{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "755d81e8",
   "metadata": {},
   "source": [
    "# RKKY\n",
    "\n",
    "Example and benchmark problem according to \"Accurate finite-difference micromagnetics of magnets including RKKY interaction\", Suess, Dieter, et al.,  arXiv preprint arXiv:2206.11063 (2022).\n",
    "\n",
    "The analytical saturation field Hx = 5*Hk"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54df1504",
   "metadata": {},
   "source": [
    "## Google Colab Link\n",
    "\n",
    "The demo can be run on Google Colab without any local installation.\n",
    "Use the following [link](https://colab.research.google.com/drive/1SIdiiz8plOI0SG3HhxNJYOxbknG178Qo) to try it out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "47f0b5e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install -q magnumnp numpy==1.22.4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80a432d8",
   "metadata": {},
   "source": [
    "## Run Demo:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "bed27907",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2023-05-08 12:26:30  magnum.np:INFO \u001b[1;37;32m[State] running on device: cpu (dtype = float64)\u001b[0m\n",
      "2023-05-08 12:26:30  magnum.np:INFO \u001b[1;37;32m[Mesh] 1x1x400 (size= 2e-09 x 2e-09 x 2e-09)\u001b[0m\n",
      "2023-05-08 12:26:30  magnum.np:INFO \u001b[1;37;32m[LLGSolver] using RKF45 solver (atol = 1e-06)\u001b[0m\n"
     ]
    },
    {
     "ename": "IndexError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "File \u001b[0;32m~/git/magnum.np/magnumnp/common/timer.py:236\u001b[0m, in \u001b[0;36mtimedmethod.<locals>.timed\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m    235\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m Timer(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m.\u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;241m%\u001b[39m (args[\u001b[38;5;241m0\u001b[39m]\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__class__\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__name__\u001b[39m, method\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__name__\u001b[39m)):\n\u001b[0;32m--> 236\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mmethod\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/git/magnum.np/magnumnp/field_terms/external.py:51\u001b[0m, in \u001b[0;36mExternalField.h\u001b[0;34m(self, state)\u001b[0m\n\u001b[1;32m     49\u001b[0m \u001b[38;5;129m@timedmethod\u001b[39m\n\u001b[1;32m     50\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mh\u001b[39m(\u001b[38;5;28mself\u001b[39m, state):\n\u001b[0;32m---> 51\u001b[0m     h \u001b[38;5;241m=\u001b[39m \u001b[43mstate\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mTensor\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_h\u001b[49m\u001b[43m)\u001b[49m\u001b[43m(\u001b[49m\u001b[43mstate\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mt\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     52\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(h\u001b[38;5;241m.\u001b[39mshape) \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m1\u001b[39m:\n",
      "File \u001b[0;32m~/git/magnum.np/magnumnp/common/state.py:115\u001b[0m, in \u001b[0;36mState.Tensor.<locals>.<lambda>\u001b[0;34m(t)\u001b[0m\n\u001b[1;32m    114\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m callable(data):\n\u001b[0;32m--> 115\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;01mlambda\u001b[39;00m t: \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mTensor(\u001b[43mdata\u001b[49m\u001b[43m(\u001b[49m\u001b[43mt\u001b[49m\u001b[43m)\u001b[49m)\n\u001b[1;32m    116\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n",
      "File \u001b[0;32m~/git/magnum.np/magnumnp/common/time_interpolator.py:36\u001b[0m, in \u001b[0;36mTimeInterpolator.__call__\u001b[0;34m(self, t)\u001b[0m\n\u001b[1;32m     35\u001b[0m fp \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_fp\n\u001b[0;32m---> 36\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m fp[i\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m] \u001b[38;5;241m+\u001b[39m (t\u001b[38;5;241m-\u001b[39mtp[i\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m]) \u001b[38;5;241m/\u001b[39m (\u001b[43mtp\u001b[49m\u001b[43m[\u001b[49m\u001b[43mi\u001b[49m\u001b[43m]\u001b[49m\u001b[38;5;241m-\u001b[39mtp[i\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m]) \u001b[38;5;241m*\u001b[39m (fp[i] \u001b[38;5;241m-\u001b[39m fp[i\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m])\n",
      "File \u001b[0;32m~/git/magnum.np/magnumnp/common/decorated_tensor.py:62\u001b[0m, in \u001b[0;36mDecoratedTensor.__getitem__\u001b[0;34m(self, idx)\u001b[0m\n\u001b[1;32m     61\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m: \u001b[38;5;66;03m# apply fancy indexing\u001b[39;00m\n\u001b[0;32m---> 62\u001b[0m     item\u001b[38;5;241m.\u001b[39m_cell_volumes \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_cell_volumes\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[38;5;21;43m__getitem__\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43midx\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     64\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m item\n",
      "File \u001b[0;32m~/git/magnum.np/venv/lib/python3.8/site-packages/torch/_tensor.py:1295\u001b[0m, in \u001b[0;36mTensor.__torch_function__\u001b[0;34m(cls, func, types, args, kwargs)\u001b[0m\n\u001b[1;32m   1294\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m _C\u001b[38;5;241m.\u001b[39mDisableTorchFunctionSubclass():\n\u001b[0;32m-> 1295\u001b[0m     ret \u001b[38;5;241m=\u001b[39m \u001b[43mfunc\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m   1296\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m func \u001b[38;5;129;01min\u001b[39;00m get_default_nowrap_functions():\n",
      "\u001b[0;31mIndexError\u001b[0m: index 1 is out of bounds for dimension 0 with size 1",
      "\nThe above exception was the direct cause of the following exception:\n",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[12], line 52\u001b[0m\n\u001b[1;32m     50\u001b[0m logger \u001b[38;5;241m=\u001b[39m Logger(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mdata\u001b[39m\u001b[38;5;124m\"\u001b[39m, [\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mt\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mm\u001b[39m\u001b[38;5;124m'\u001b[39m, zeeman\u001b[38;5;241m.\u001b[39mh], [\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mm\u001b[39m\u001b[38;5;124m'\u001b[39m], fields_every \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m100\u001b[39m)\n\u001b[1;32m     51\u001b[0m \u001b[38;5;28;01mwhile\u001b[39;00m state\u001b[38;5;241m.\u001b[39mt \u001b[38;5;241m<\u001b[39m tfinal:\n\u001b[0;32m---> 52\u001b[0m     \u001b[43mlogger\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m<<\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mstate\u001b[49m\n\u001b[1;32m     53\u001b[0m     llg\u001b[38;5;241m.\u001b[39mstep(state, \u001b[38;5;241m1e-9\u001b[39m)\n\u001b[1;32m     55\u001b[0m Timer\u001b[38;5;241m.\u001b[39mprint_report()\n",
      "File \u001b[0;32m~/git/magnum.np/magnumnp/loggers/logger.py:71\u001b[0m, in \u001b[0;36mLogger.__lshift__\u001b[0;34m(self, state)\u001b[0m\n\u001b[1;32m     70\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__lshift__\u001b[39m(\u001b[38;5;28mself\u001b[39m, state):\n\u001b[0;32m---> 71\u001b[0m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlog\u001b[49m\u001b[43m(\u001b[49m\u001b[43mstate\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/git/magnum.np/magnumnp/loggers/logger.py:68\u001b[0m, in \u001b[0;36mLogger.log\u001b[0;34m(self, state)\u001b[0m\n\u001b[1;32m     66\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m\n\u001b[1;32m     67\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m logger \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mloggers\u001b[38;5;241m.\u001b[39mvalues():\n\u001b[0;32m---> 68\u001b[0m     \u001b[43mlogger\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m<<\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mstate\u001b[49m\n",
      "File \u001b[0;32m~/git/magnum.np/magnumnp/loggers/scalar_logger.py:113\u001b[0m, in \u001b[0;36mScalarLogger.__lshift__\u001b[0;34m(self, state)\u001b[0m\n\u001b[1;32m    112\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__lshift__\u001b[39m(\u001b[38;5;28mself\u001b[39m, state):\n\u001b[0;32m--> 113\u001b[0m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlog\u001b[49m\u001b[43m(\u001b[49m\u001b[43mstate\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/git/magnum.np/magnumnp/loggers/scalar_logger.py:92\u001b[0m, in \u001b[0;36mScalarLogger.log\u001b[0;34m(self, state)\u001b[0m\n\u001b[1;32m     90\u001b[0m     \u001b[38;5;28;01mexcept\u001b[39;00m:\n\u001b[1;32m     91\u001b[0m         name \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124munnamed\u001b[39m\u001b[38;5;124m'\u001b[39m\n\u001b[0;32m---> 92\u001b[0m     raw_value \u001b[38;5;241m=\u001b[39m \u001b[43mcolumn\u001b[49m\u001b[43m(\u001b[49m\u001b[43mstate\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     93\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(column, \u001b[38;5;28mtuple\u001b[39m) \u001b[38;5;129;01mor\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(column, \u001b[38;5;28mlist\u001b[39m):\n\u001b[1;32m     94\u001b[0m     name \u001b[38;5;241m=\u001b[39m column[\u001b[38;5;241m0\u001b[39m]\n",
      "File \u001b[0;32m~/git/magnum.np/magnumnp/common/timer.py:236\u001b[0m, in \u001b[0;36mtimedmethod.<locals>.timed\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m    233\u001b[0m \u001b[38;5;129m@wraps\u001b[39m(method)\n\u001b[1;32m    234\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mtimed\u001b[39m(\u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[1;32m    235\u001b[0m     \u001b[38;5;28;01mwith\u001b[39;00m Timer(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m.\u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;241m%\u001b[39m (args[\u001b[38;5;241m0\u001b[39m]\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__class__\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__name__\u001b[39m, method\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__name__\u001b[39m)):\n\u001b[0;32m--> 236\u001b[0m         \u001b[38;5;28;01mreturn\u001b[39;00m method(\u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n",
      "File \u001b[0;32m~/git/magnum.np/magnumnp/common/timer.py:103\u001b[0m, in \u001b[0;36mTimer.__exit__\u001b[0;34m(self, e_typ, e_val, trcbak)\u001b[0m\n\u001b[1;32m    101\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__exit__\u001b[39m(\u001b[38;5;28mself\u001b[39m, e_typ, e_val, trcbak):\n\u001b[1;32m    102\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mall\u001b[39m((e_typ, e_val, trcbak)):\n\u001b[0;32m--> 103\u001b[0m         \u001b[38;5;28;01mraise\u001b[39;00m e_typ \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01me_val\u001b[39;00m\n\u001b[1;32m    105\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m Timer\u001b[38;5;241m.\u001b[39m_options[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mactive\u001b[39m\u001b[38;5;124m'\u001b[39m]: \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\n\u001b[1;32m    107\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m Timer\u001b[38;5;241m.\u001b[39m_options[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mskip\u001b[39m\u001b[38;5;124m'\u001b[39m] \u001b[38;5;241m<\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_data[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mcalls\u001b[39m\u001b[38;5;124m'\u001b[39m]:\n",
      "\u001b[0;31mIndexError\u001b[0m: "
     ]
    }
   ],
   "source": [
    "from magnumnp import *\n",
    "import torch\n",
    "import numpy as np\n",
    "\n",
    "Timer.enable()\n",
    "A = 1e-11\n",
    "Js = 1.0\n",
    "Ms = Js/constants.mu_0\n",
    "Ku = 1e5\n",
    "tfinal = 20e-9\n",
    "J_rkky = -2.0e-3\n",
    "Hxmin = 4.8 * 2*Ku/Js\n",
    "Hxmax = 5.0 * 2*Ku/Js\n",
    "\n",
    "# initialize mesh\n",
    "n  = (1, 1, 400)\n",
    "dx = (2e-9, 2e-9, 2e-9)\n",
    "mesh = Mesh(n, dx)\n",
    "state = State(mesh)\n",
    "\n",
    "state.material = {\n",
    "        \"Ms\": Ms,\n",
    "        \"A\": A,\n",
    "        \"Ku\": Ku,\n",
    "        \"Ku_axis\": state.Tensor([0,1,0]),\n",
    "        \"alpha\": 1.0\n",
    "        }\n",
    "\n",
    "domain1 = state._zeros(n, dtype=torch.bool)\n",
    "domain1[:,:,n[2]//2:] = True\n",
    "\n",
    "domain2 = state._zeros(n, dtype=torch.bool)\n",
    "domain2[:,:,:-n[2]//2] = True\n",
    "\n",
    "# Set initial magnetization\n",
    "state.m = state.Constant([0, 0, 0]) \n",
    "state.m[domain1] = state.Tensor([0, -1, 0])\n",
    "state.m[domain2] = state.Tensor([0, 1, 0])\n",
    "\n",
    "# define field terms\n",
    "exchange1 = ExchangeField(domain1) # Two seperate exchange regions are required, \n",
    "exchange2 = ExchangeField(domain2) # so that at the RKKY interface the bulk exchange is zero\n",
    "rkky      = RKKYField(J_rkky, \"z\", n[2]//2-1, n[2]//2, order=2)\n",
    "aniso     = UniaxialAnisotropyField()\n",
    "zeeman    = ExternalField(TimeInterpolator(state, {0e-0:   [0, 0, Hxmin],\n",
    "                                                   tfinal: [0, 0, Hxmax]}))\n",
    "\n",
    "# integrate\n",
    "llg = LLGSolver([aniso, exchange1, exchange2, rkky, zeeman], atol = 1e-6)\n",
    "logger = Logger(\"data\", ['t', 'm', zeeman.h], ['m'], fields_every = 100)\n",
    "while state.t < tfinal:\n",
    "    logger << state\n",
    "    llg.step(state, 1e-9)\n",
    "\n",
    "Timer.print_report()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f08871ce",
   "metadata": {},
   "source": [
    "## Plot Results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3dc1640",
   "metadata": {},
   "outputs": [],
   "source": [
    "from os import path\n",
    "if not path.isfile(\"ref/log.dat\"):\n",
    "    !mkdir ref\n",
    "    !wget -P ref https://gitlab.com/magnum.np/magnum.np/raw/main/demos/rkky/ref/log.dat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f07a8913",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "data = np.loadtxt(\"data/log.dat\")\n",
    "ref = np.loadtxt(\"ref/log.dat\")\n",
    "\n",
    "Hk = 2 * 1e5 / 1.\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10,5))\n",
    "cycle = plt.rcParams['axes.prop_cycle'].by_key()['color']\n",
    "\n",
    "ax.plot(data[:,6]/Hk, data[:,3], '-', color = cycle[2], label = \"magnum.np\")\n",
    "ax.plot(ref[:,6]/Hk, ref[:,3], '-', color = cycle[2], linewidth = 6, alpha = 0.4, label = \"reference\")\n",
    "\n",
    "ax.set_xlim([4.82,5.0])\n",
    "ax.set_ylim([0.99994,1.0])\n",
    "ax.set_yticklabels(['0.99994', '0.99995', '0.99996', '0.99997', '0.99998', '0.99999', '1.00000'])\n",
    "ax.set_title(\"RKKY Standard Problem\")\n",
    "ax.set_xlabel(\"External Field H$_{ext}$ / H$_k$ [1]\")\n",
    "ax.set_ylabel(\"Magnetization m\")\n",
    "ax.legend(ncol=3)\n",
    "ax.grid()\n",
    "fig.savefig(\"data/results.png\")\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}