{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6eccf3f7",
   "metadata": {
    "tags": [
     "remove-input"
    ]
   },
   "outputs": [],
   "source": [
    "%config InlineBackend.figure_formats = ['svg']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34ad4fbd",
   "metadata": {},
   "source": [
    "# Stochastic thermal field\n",
    "## Verification experiment\n",
    "\n",
    "In this notebook we set up the flatspin runs and analyse the data for the verification of the stochastic thermal field, as presented in the paper.\n",
    "In short, we compare flatspin simulations with varying fields and temperatures to experiments where the results are known analytically (or numerically through micromagnetic simulations).\n",
    "\n",
    "We consider a typical ensemble of uncoupled spins with magnetization direction parallel to an external magnetic field direction.\n",
    "The ensembles are subject to both temperature and an external field influence.\n",
    "In the lack of other influences and spin-spin interactions, the ensemble net magnetization is completely defined by temperature and the external field $M(H,T)$.\n",
    "\n",
    "The experiment has a low coercivity and a high coercivity scenario.\n",
    "The low coercivity scenario, $h_k$ = 0.001 mT, can be described analytically, by $\\langle m_x \\rangle = \\tanh (A\\mu_0 H)$, where $A = M_SV/K_BT)$. The high coercivity scenario, $h_k$ = 0.020 mT, we model by micromagnetic simulations."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39b41d6d",
   "metadata": {},
   "source": [
    "### The low coercivity scenario\n",
    "\n",
    "First, we look at the analytical description of the low coercivity scenario.\n",
    "We simply plot the following expression for $m_x$:\n",
    "\n",
    "$\\langle m_x \\rangle = \\tanh (A\\mu_0 H)$, where $A = M_SV/K_BT)$,\n",
    "\n",
    "using the volume and saturation magnetization that we will use in the experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "972905a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "kB = 1.38064852e-23\n",
    "msat = 860e3\n",
    "volume = 10e-9*10e-9*10e-9\n",
    "temperatures = [100, 300, 500]\n",
    "\n",
    "def hyptan(msat, V, T, h):\n",
    "    A = msat * V / (kB * T)\n",
    "    return np.tanh(A * h)\n",
    "\n",
    "h_range = np.linspace(-0.025, 0.025, 200) # Unit Tesla, i.e. mu0 * H\n",
    "for T in temperatures:\n",
    "    plt.plot(h_range, hyptan(msat, volume, T, h_range), label=f'$T={T}$ K')\n",
    "    \n",
    "plt.xlabel('$\\mu_0H$ [T]')\n",
    "plt.ylabel('$m_x$')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8eefe80e",
   "metadata": {},
   "source": [
    "Note that this analytical expression does not take into account any coercivty, as it assumes no coercive field.\n",
    "Therefore, this only applies to the low coercivity scenario."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0760c8c2",
   "metadata": {},
   "source": [
    "### Generate flatspin data\n",
    "\n",
    "To set up the experiment in flatspin, we use the `IsingSpinIce` class, as all the spins are parallel.\n",
    "`alpha` is set to 0 to isolate each spin, and we use a 16 by 16 spin ensemble, initialized in the negative direction.\n",
    "\n",
    "Below we build the flatspin-run command by going through the relevant parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4a53382",
   "metadata": {},
   "outputs": [],
   "source": [
    "from flatspin.model import IsingSpinIce\n",
    "\n",
    "run_command = 'flatspin-run-sweep -m IsingSpinIce'\n",
    "\n",
    "alpha = 0\n",
    "size = (16, 16)\n",
    "\n",
    "model = IsingSpinIce(alpha=alpha, size=size, init=-1, spin_angle=0)\n",
    "model.plot()\n",
    "plt.title('Initial state');\n",
    "\n",
    "run_command += f\" -p alpha={alpha} -p size='{size}' -p init=-1 -p spin_angle=0\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be04d812",
   "metadata": {},
   "source": [
    "We use a quasi-static field protocol, i.e., we increase the field by a small amount and simulate many timesteps at each field value to reach a steady state. We also save the state at each timestep (i.e. `spp = timesteps`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ae1d69f",
   "metadata": {},
   "outputs": [],
   "source": [
    "H_max = 0.1 # Tesla, plenty enough to saturate the ensembles\n",
    "timesteps_per_H  = 200 # We repeat each field value this many times\n",
    "\n",
    "H_range = H_max * np.linspace(-1,1,1001) # The number of inputs\n",
    "H_profile = np.repeat(H_range, timesteps_per_H)\n",
    "\n",
    "plt.plot(H_profile)\n",
    "plt.ylabel('$H$ [T]')\n",
    "plt.xlabel('Step');\n",
    "\n",
    "run_command += f\" -e Constant -p H={H_max} -p 'input=np.linspace(-1,1,1001)'\"\n",
    "run_command += f\" -p timesteps={timesteps_per_H} -p spp={timesteps_per_H}\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a73a7ff5",
   "metadata": {},
   "outputs": [],
   "source": [
    "hc_low = 0.001\n",
    "hc_high = 0.020\n",
    "delta_t_low = 1e-10\n",
    "delta_t_high = 1e-9\n",
    "field_angles = [0]\n",
    "\n",
    "# Add temperature sweep and other parameters\n",
    "run_command += f\" -s 'temperature={temperatures}' -p msat={msat} -p volume={volume} -p 'm_therm=msat*volume' -s 'phi={field_angles}'\"\n",
    "# Astroid parameters\n",
    "run_command += \" -p sw_b=1 -p sw_c=1 -p sw_gamma=3 -p sw_beta=3\"\n",
    "\n",
    "# Specifiy distributed running on GPUs\n",
    "run_command += \" -p use_opencl=1 -r dist\"\n",
    "\n",
    "# Create one command for the low hc and one for the high hc\n",
    "run_command_low = run_command + f\" -s 'therm_timescale=[{delta_t_low}]' -p hc={hc_low} -o OUTPUT_FILE_LOW\" \n",
    "run_command_high = run_command + f\" -s 'therm_timescale=[{delta_t_high}]' -p hc={hc_high} -o OUTPUT_FILE_HIGH\" \n",
    "\n",
    "print(run_command_low)\n",
    "print()\n",
    "print(run_command_high)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4c6b97d",
   "metadata": {},
   "source": [
    "The above commands can be used to run and generate {download}`the relevant datasets<temp-verification.zip>`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7374e85",
   "metadata": {},
   "source": [
    "### Results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9368bf3b",
   "metadata": {},
   "source": [
    "With our data generated, we read the datasets and take a closer look"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e72737f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from flatspin.data import Dataset, read_table\n",
    "from numpy.linalg import norm\n",
    "\n",
    "# Selected parameters\n",
    "phi = 0\n",
    "temperatures = [100,300,500]\n",
    "\n",
    "# Read data\n",
    "flatspin_ds_high_hc = Dataset.read('/data/flatspin/temp-verification/temp-verification-high-hc')\n",
    "flatspin_ds_low_hc = Dataset.read('/data/flatspin/temp-verification/temp-verification-low-hc')\n",
    "\n",
    "# If the dataset includes extra swept parameters, we filter these out:\n",
    "flatspin_ds_high_hc = flatspin_ds_high_hc.filter(phi=round(phi), temperature=temperatures, therm_timescale=[1e-9]) \n",
    "flatspin_ds_low_hc = flatspin_ds_low_hc.filter(phi=round(phi), temperature=temperatures, therm_timescale=[1e-10])\n",
    "\n",
    "\n",
    "# Function to return a dataframe with the relevant field and magnetization from a single dataset\n",
    "def read_hyst_data(ds):\n",
    "    df = read_table(ds.tablefile(\"h_ext\"), index_col=\"t\")\n",
    "    df['h'] = np.sign(df[\"h_extx\"]) * norm(df[[\"h_extx\", \"h_exty\"]].values, axis=1)\n",
    "    \n",
    "    # Average magnetization\n",
    "    spin_states = read_table(ds.tablefile(\"spin\"), index_col=\"t\")\n",
    "    df[\"m\"] = spin_states.agg('mean', axis=1)\n",
    "    \n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "212057eb",
   "metadata": {},
   "source": [
    "We also read in some data that we have generated with micromagnetic simulations.\n",
    "The data is saved as a csv with the relevant temperature and coercivity (hc). \n",
    "The hc field has two values, 'high' and 'low', indicating a coercivity equivalent to hc=0.020 mT and hc=0.001 mT respectively. \n",
    "For each external field value, 'B', the average of the binned spin states of each micromagnetic cell, 'sx', is recorded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e0a56dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "mx = pd.read_csv('mumax-temp-vs-field.csv', usecols=['T', 'hc', 'B', 'sx'])\n",
    "mx = mx[mx['hc'] == 'high'] # The comparison is only valid for the high hc scenario"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "335c2bc9",
   "metadata": {},
   "source": [
    "### Plotting everything\n",
    "\n",
    "We average the magnetization at each field value, and plot everything together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ef88888",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(7.08*1.135, 4.395))\n",
    "\n",
    "H_axis = mx[mx['T']==temperatures[0]]['B']\n",
    "### high hc subplot ###\n",
    "plt.subplot(1,2,2)\n",
    "\n",
    "hc = flatspin_ds_high_hc.params['hc']*1000\n",
    "plt.title(f'$h_c={hc}$ mT')\n",
    "plt.xlim(-.025, .025)\n",
    "plt.xlabel(\"$\\mu_0H$\")\n",
    "\n",
    "for i, T in enumerate(temperatures):\n",
    "    # Read the temperature dataset\n",
    "    df = read_hyst_data(flatspin_ds_high_hc.filter(temperature=T))\n",
    "    # Take the mean for each field value\n",
    "    df = df.groupby('h').agg('mean').reset_index()\n",
    "    \n",
    "    # Plot flatpsin results\n",
    "    plt.plot(df['h'], df['m'], label=f'$T={T}$ K, flatspin', c=f'C{i}', ls=':' )\n",
    "    \n",
    "    # Plot MuMax3 results\n",
    "    plt.plot(mx[mx['T']==T]['B'], \n",
    "             mx[mx['T']==T]['sx'], label=f'$T={T}$ K, MuMax3', c=f'C{i}', ls='--')\n",
    "    \n",
    "    # Plot analytical\n",
    "    ht = hyptan(msat, volume, T, H_axis)\n",
    "    plt.plot(H_axis, ht, label=f'$T={T}$ K, analytical', ls='-',alpha=0.5, linewidth=2 )\n",
    "\n",
    "ax = plt.gca()\n",
    "h, l = ax.get_legend_handles_labels()\n",
    "\n",
    "### low hc subplot ###\n",
    "plt.subplot(1,2,1)\n",
    "\n",
    "hc = flatspin_ds_low_hc.params['hc']*1000\n",
    "plt.title(f'$h_c={hc}$ mT')\n",
    "plt.xlim(-.025, .025)\n",
    "plt.xlabel(\"$\\mu_0H$\")\n",
    "plt.ylabel('$m_x$')\n",
    "\n",
    "for i, T in enumerate(temperatures):\n",
    "    df = read_hyst_data(flatspin_ds_low_hc.filter(temperature=T))\n",
    "    df = df.groupby('h').agg('mean').reset_index()\n",
    "    \n",
    "    # Plot flatspin results\n",
    "    plt.plot(df['h'], df['m'], label=f'$T={T}$ K, flatspin', c=f'C{i}', ls=':')\n",
    "\n",
    "    # Plot analytical\n",
    "    ht = hyptan(msat, volume, T, H_axis)\n",
    "    plt.plot(H_axis, ht, label=f'$T={T}$ K, analytical', ls='-',alpha=0.5, linewidth=2 )\n",
    "\n",
    "\n",
    "fig = plt.gcf()\n",
    "fig.subplots_adjust(bottom=0.25, wspace=0.25)\n",
    "fig.legend(h, l, loc='lower center', bbox_to_anchor=(0.5, -0.05), ncol=3)\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,md:myst",
   "text_representation": {
    "extension": ".md",
    "format_name": "myst",
    "format_version": 0.13,
    "jupytext_version": "1.11.5"
   }
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  },
  "source_map": [
   15,
   19,
   34,
   45,
   65,
   70,
   79,
   92,
   96,
   111,
   133,
   137,
   141,
   145,
   173,
   180,
   183,
   189
  ]
 },
 "nbformat": 4,
 "nbformat_minor": 5
}