standard_problem3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Micromagnetic standard problem 3\n",
    "\n",
    "![title](imgs/stdprob3.png)\n",
    "\n",
    "## Problem specification\n",
    "\n",
    "This problem is to calculate a single domain limit of a cubic magnetic particle. This is the size $L$ of equal energy for the so-called flower state (which one may also call a splayed state or a modified single-domain state) on the one hand, and the vortex or curling state on the other hand.\n",
    "\n",
    "Geometry:\n",
    "\n",
    "A cube with edge length, $L$, expressed in units of the intrinsic length scale, $l_\\text{ex} = \\sqrt{A/K_\\text{m}}$, where $K_\\text{m}$ is a magnetostatic energy density, $K_\\text{m} = \\frac{1}{2}\\mu_{0}M_\\text{s}^{2}$.\n",
    "\n",
    "Material parameters: \n",
    "\n",
    "- uniaxial anisotropy $K_\\text{u}$ with $K_\\text{u} = 0.1 K_\\text{m}$, and with the easy axis directed parallel to a principal axis of the cube (0, 0, 1),\n",
    "- exchange energy constant is $A = \\frac{1}{2}\\mu_{0}M_\\text{s}^{2}l_\\text{ex}^{2}$.\n",
    "\n",
    "More details about the standard problem 3 can be found in Ref. 1.\n",
    "\n",
    "## Simulation\n",
    "\n",
    "Firstly, we import all necessary modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import discretisedfield as df\n",
    "import oommfc as oc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following two functions are used for initialising the system's magnetisation [1]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Function for initialising the flower state.\n",
    "def m_init_flower(pos):\n",
    "    x, y, z = pos[0]/1e-9, pos[1]/1e-9, pos[2]/1e-9\n",
    "    mx = 0\n",
    "    my = -2*z + 1\n",
    "    mz = 2*y - 1\n",
    "    norm_squared = mx**2 + my**2 + mz**2\n",
    "    if norm_squared <= 0.05:\n",
    "        return (1, 0, 0)\n",
    "    else:\n",
    "        return (mx, my, mz)\n",
    "\n",
    "# Function for initialising the vortex state.\n",
    "def m_init_vortex(pos):\n",
    "    x, y, z = pos[0]/1e-9, pos[1]/1e-9, pos[2]/1e-9\n",
    "    mx = 0\n",
    "    my = np.sin(np.pi/2 * (-x-0.5))\n",
    "    mz = np.cos(np.pi/2 * (-x-0.5))\n",
    "    \n",
    "    return (mx, my, mz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function is used for convenience. It takes two arguments:\n",
    "\n",
    "- $L$ - the cube edge length in units of $l_\\text{ex}$, and\n",
    "- the function for initialising the system's magnetisation.\n",
    "\n",
    "It returns the relaxed system object.\n",
    "\n",
    "Please refer to other tutorials for more details on how to create system objects and drive them using specific drivers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def minimise_system_energy(L, m_init):\n",
    "    \"\"\"Given an edge lengch L, and a function m_init for the initial configuration,\n",
    "    relax the system and return the system object (from which the energy can be \n",
    "    calculated)\"\"\"\n",
    "    print(\"L={:9}, {} \".format(L, m_init.__name__), end=\"\")\n",
    "    N = 16  # discretisation in one dimension\n",
    "    cubesize = 100e-9  # cube edge length (m)\n",
    "    cellsize = cubesize/N  # discretisation in all three dimensions.\n",
    "    lex = cubesize/L  # exchange length.\n",
    "    \n",
    "    Km = 1e6  # magnetostatic energy density (J/m**3)\n",
    "    Ms = np.sqrt(2*Km/oc.mu0)  # magnetisation saturation (A/m)\n",
    "    A = 0.5 * oc.mu0 * Ms**2 * lex**2  # exchange energy constant\n",
    "    K = 0.1*Km  # Uniaxial anisotropy constant\n",
    "    u = (0, 0, 1)  # Uniaxial anisotropy easy-axis\n",
    "\n",
    "    p1 = (0, 0, 0)  # Minimum sample coordinate.\n",
    "    p2 = (cubesize, cubesize, cubesize)  # Maximum sample coordinate.\n",
    "    cell = (cellsize, cellsize, cellsize)  # Discretisation.\n",
    "    mesh = oc.Mesh(p1=(0, 0, 0), p2=(cubesize, cubesize, cubesize),\n",
    "                   cell=(cellsize, cellsize, cellsize))  # Create a mesh object.\n",
    "\n",
    "    system = oc.System(name=\"stdprob3\")\n",
    "    system.hamiltonian = oc.Exchange(A) + oc.UniaxialAnisotropy(K, u) + oc.Demag()\n",
    "    system.m = df.Field(mesh, value=m_init, norm=Ms)\n",
    "\n",
    "    md = oc.MinDriver()\n",
    "    md.drive(system)\n",
    "    \n",
    "    return system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Relaxed magnetisation states: Vortex state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "system = minimise_system_energy(8., m_init_vortex)\n",
    "system.m.plot_plane(\"y\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Relaxed magnetisation states: Flower state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "system = minimise_system_energy(8, m_init_flower)\n",
    "system.m.plot_plane(\"y\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Energy crossing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "system.total_energy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can plot the energies of both vortex and flower states as a function of cube edge length. This will give us an idea where the state transition occurrs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L_array = np.linspace(8, 9, 5)  # values of L for which the system is relaxed.\n",
    "\n",
    "vortex_energies = []\n",
    "flower_energies = []\n",
    "\n",
    "for L in L_array:\n",
    "    vortex = minimise_system_energy(L, m_init_vortex)\n",
    "    flower = minimise_system_energy(L, m_init_flower)\n",
    "    \n",
    "    vortex_energies.append(vortex.total_energy())\n",
    "    flower_energies.append(flower.total_energy())\n",
    "    \n",
    "# Plot the energy dependences.\n",
    "import matplotlib.pyplot as plt\n",
    "plt.plot(L_array, vortex_energies, 'o-', label='vortex')\n",
    "plt.plot(L_array, flower_energies, 'o-', label='flower')\n",
    "plt.xlabel('L (lex)')\n",
    "plt.ylabel('E')\n",
    "plt.xlim([8.0, 9.0])\n",
    "plt.grid()\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now know that the energy crossing occurrs between $8l_\\text{ex}$ and $9l_\\text{ex}$, so a bisection algorithm can be used to find the exact crossing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import bisect\n",
    "\n",
    "def energy_difference(L):\n",
    "    vortex = minimise_system_energy(L, m_init_vortex)\n",
    "    flower = minimise_system_energy(L, m_init_flower)\n",
    "    \n",
    "    return vortex.total_energy() - flower.total_energy()\n",
    "\n",
    "crossing = bisect(energy_difference, 8, 9, xtol=0.1)\n",
    "\n",
    "print(\"The transition between vortex and flower states occurs at {}*lex\".format(crossing))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "[1] µMAG Site Directory http://www.ctcms.nist.gov/~rdm/mumag.org.html"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}