Merge branch 'master' of github.com:computationalmodelling/fidimag

davidcorteso · davidcorteso · commit 74cbb147b405 · 2016-07-12T12:32:56.000+01:00
diff --git a/doc/index.rst b/doc/index.rst
@@ -4,7 +4,7 @@
    contain the root `toctree` directive.
 
 Welcome to fidimag's documentation!
-================================
+===================================
 
 The official webpage is
 http://computationalmodelling.github.io/fidimag/, the github
@@ -17,6 +17,7 @@ Contents:
    :maxdepth: 2
 
    install
+   ipynb/tutorial-docker-container
    core_eqs
    extended_eqs
    monte_carlo
diff --git a/doc/ipynb/tutorial-docker-container.ipynb b/doc/ipynb/tutorial-docker-container.ipynb
@@ -28,10 +28,10 @@
    "source": [
     "## Setup fidimag docker container with Jupyter Notebook\n",
     "\n",
-    "Pull the `fidimag notebook` container:\n",
+    "Pull the fidimag notebook container:\n",
     "\n",
     "```\n",
-    "docker pull fidimag\n",
+    "docker pull fidimag/notebook\n",
     "    ```"
    ]
   },
@@ -60,9 +60,9 @@
     "```\n",
     "docker-machine ip\n",
     "```\n",
-    "For example, if `docker-machine ip` returns `192.168.99.100`, then the right URL to paste into the browser on the host system is `http://192.168.99.100/30000`.\n",
+    "For example, if `docker-machine ip` returns `192.168.99.100`, then the right URL to paste into the browser on the host system is `http://192.168.99.100:30000`.\n",
     "\n",
-    "[How does this work on Windows? Pull requests welcome.]"
+    "[*How does this work on Windows? Pull requests welcome.*]"
    ]
   },
   {
@@ -84,10 +84,13 @@
     "## Show active Docker containers\n",
     "```docker ps``` lists all running containers.\n",
     "\n",
-    "To only show the `id`s, we can use ```docker ps -q```\n",
+    "To only show the `id`s, we can use \n",
+    "\n",
+    "    docker ps -q\n",
     "\n",
     "To only show the containers that was last started, we can use the `-l` flag:\n",
-    "``` docker ps -l```"
+    "\n",
+    "    docker ps -l"
    ]
   },
   {
@@ -163,7 +166,7 @@
     "We can use the container like this to achieve that:\n",
     "\n",
     "```\n",
-    "docker run -v `pwd`:/io -ti fidimag/notebook python hello.py```\n",
+    "docker run -v `pwd`:/io -ti fidimag/notebook python run.py```\n",
     "\n",
     "The ```-v `pwd`:/io``` tells the docker container to take the current working directory on the host (``` `pwd` ```) and mount it to the path `/io` on the container. The container is set up so that the default working directory is `/io`."
    ]
diff --git a/fidimag/atomistic/lib/clib.h b/fidimag/atomistic/lib/clib.h
@@ -60,8 +60,8 @@ void llg_s_rhs(double * dm_dt, double * spin, double * h, double *alpha,
              double *chi, double gamma, int n);
 
 
-void compute_stt_field_c(double *spin, double *field, double *jx, double *jy,
-		double dx, double dy, int *ngbs, int n);
+void compute_stt_field_c(double *spin, double *field, double *jx, double *jy, double *jz,
+		double dx, double dy, double dz, int *ngbs, int n);
 
 void llg_stt_rhs(double *dm_dt, double *m, double *h,
                  double *h_stt, double *alpha,
diff --git a/fidimag/atomistic/lib/clib.pyx b/fidimag/atomistic/lib/clib.pyx
@@ -78,8 +78,8 @@ cdef extern from "clib.h":
     void normalise(double *m, int *pins, int n)
 
     void compute_stt_field_c(double *spin, double *field,
-                             double *jx, double *jy,
-                             double dx, double dy, int *ngbs, int n)
+                             double *jx, double *jy, double *jz,
+                             double dx, double dy, double dz, int *ngbs, int n)
 
     void llg_stt_rhs(double *dm_dt, double *m, double *h, double *h_stt,
                      double *alpha,double beta, double u0, double gamma, int n)
@@ -241,12 +241,13 @@ def compute_stt_field(np.ndarray[double, ndim=1, mode="c"] spin,
                       np.ndarray[double, ndim=1, mode="c"] field,
                       np.ndarray[double, ndim=1, mode="c"] jx,
                       np.ndarray[double, ndim=1, mode="c"] jy,
-                      dx, dy,
+                      np.ndarray[double, ndim=1, mode="c"] jz,
+                      dx, dy, dz,
                       np.ndarray[int, ndim=2, mode="c"] ngbs,
                       n
                       ):
-    compute_stt_field_c(&spin[0], &field[0], &jx[0], &jy[0],
-                        dx, dy, &ngbs[0, 0], n)
+    compute_stt_field_c(&spin[0], &field[0], &jx[0], &jy[0],&jz[0],
+                        dx, dy, dz, &ngbs[0, 0], n)
 
 def compute_llg_stt_rhs(np.ndarray[double, ndim=1, mode="c"] dm_dt,
                 np.ndarray[double, ndim=1, mode="c"] spin,
diff --git a/fidimag/atomistic/lib/stt.c b/fidimag/atomistic/lib/stt.c
@@ -4,7 +4,7 @@
  * is
  *          ( \vec{j}_s \cdot \nabla ) \vec{m}
  *
- * where j_s is the current vector in 2D
+ * where j_s is the current vector in 3D
  * For this computation we use the neighbours matrix to make easier the
  * derivative calculation at the boundaries.
  * ** For now this only works with the CUBOID mesh
@@ -17,6 +17,7 @@
  *  the derivatives are simply:
  *          f(x + 1, y, z) - f(x - 1, y, z) / 2 dx
  *          f(x, y + 1, z) - f(x, y - 1, z) / 2 dy
+ *          f(x, y, z + 1) - f(x, y, z - 1) / 2 dz
  *
  *  where x+-1 is the neighbour to the left or right.
  *
@@ -31,7 +32,7 @@
  * i-th site and nn_x2 for the NN to the right. These variables are simply
  * the i-th index (in the corresponding cases) when there is a single NN
  *
- * Similar for y
+ * Similar for y and z
  *
  * n is the number of spins or lattice sites in the system
  *
@@ -50,25 +51,24 @@
  *  Neighbouring sites where there is no material, has index -1
  *
  */
-void compute_stt_field_c(double *spin, double *field, double *jx, double *jy,
-		double dx, double dy, int *ngbs, int n) {
+void compute_stt_field_c(double *spin, double *field, double *jx, double *jy, double *jz,
+		double dx, double dy, double dz, int *ngbs, int n) {
 
-    int i, j;
-    float factor_x, factor_y;
-    int nn_i;
-    int nn_x1, nn_x2, nn_y1, nn_y2;
-
-	for (i = 0; i < 3 * n; i++) {
+    //#pragma omp parallel for
+	for (int i = 0; i < 3 * n; i++) {
 		field[i] = 0;
 	}
-
+    
+    #pragma omp parallel for
     /* Iterate through every lattice site */
-    for (i = 0; i < n; i++){
+    for (int i = 0; i < n; i++){
 
         /* Starting index for the NNs of the i-th site
          * i+0, i+1, i+2, i+3 ...  --> -x, +x, -y, +y ...
          */
-        nn_i = 6 * i;
+        int nn_i = 6 * i;
+        double factor_x, factor_y, factor_z;
+        int nn_x1, nn_x2, nn_y1, nn_y2, nn_z1, nn_z2;
 
         /* Here we distinguish if there are 2 NNs, no NN in the
          * -x direction, or no NN in the +x direction,
@@ -100,7 +100,7 @@ void compute_stt_field_c(double *spin, double *field, double *jx, double *jy,
          * This calculation is:  jx[i] * d m[i] / dx
          * */
         if (factor_x){
-            for(j = 0; j < 3; j++){
+            for(int j = 0; j < 3; j++){
                 field[3 * i + j] += jx[i] * (spin[3 * nn_x2 + j]
                                              - spin[3 * nn_x1 + j]) / (factor_x * dx);
             }
@@ -124,11 +124,37 @@ void compute_stt_field_c(double *spin, double *field, double *jx, double *jy,
         }
 
         if (factor_y){
-            for(j = 0; j < 3; j++){
+            for(int j = 0; j < 3; j++){
                 field[3 * i + j] += jy[i] * (spin[3 * nn_y2 + j]
                                              - spin[3 * nn_y1 + j]) / (factor_y * dy);
             }
         }
+        
+        
+        // We do the same along the z direction
+        if(ngbs[nn_i + 4] >= 0 && ngbs[nn_i + 5] >= 0) {
+            factor_z = 2;
+            nn_z1 = ngbs[nn_i + 4];
+            nn_z2 = ngbs[nn_i + 5];
+        } else if(ngbs[nn_i + 4] >= 0 && ngbs[nn_i + 5] < 0){
+            factor_z = 1;
+            nn_z1 = ngbs[nn_i + 4];
+            nn_z2 = i;
+        } else if(ngbs[nn_i + 4] < 0 && ngbs[nn_i + 5] >= 0 ){
+            factor_z = 1;
+            nn_z1 = i;
+            nn_z2 = ngbs[nn_i + 5];
+        } else {
+            factor_z = 0;
+        }
+        
+        if (factor_z){
+            for(int j = 0; j < 3; j++){
+                field[3 * i + j] += jz[i] * (spin[3 * nn_z2 + j]
+                                             - spin[3 * nn_z1 + j]) / (factor_z * dz);
+            }
+        }
+        
     }
 }
 
diff --git a/fidimag/atomistic/llg.py b/fidimag/atomistic/llg.py
@@ -86,6 +86,7 @@ def set_default_options(self, gamma=1, mu_s=1, alpha=0.1):
         self.default_c = -1
         self._alpha[:] = alpha
         self._mu_s[:] = mu_s
+        self.mu_s_const = mu_s
         self.gamma = gamma
         self.do_precession = True
 
@@ -156,6 +157,8 @@ def set_mu_s(self, value):
                 # Set the neighbour index to -1 for sites with mu_s = 0
                 self.mesh.neighbours[self.mesh.neighbours == i] = -1
 
+        self.mu_s_const = np.max(self._mu_s)
+
     mu_s = property(get_mu_s, set_mu_s)
 
     def add(self, interaction, save_field=False):
diff --git a/fidimag/atomistic/llg_stt.py b/fidimag/atomistic/llg_stt.py
@@ -0,0 +1,103 @@
+from __future__ import division
+
+import fidimag.extensions.clib as clib
+import numpy as np
+
+from .llg import LLG
+import fidimag.common.helper as helper
+import fidimag.common.constant as const
+
+
+class LLG_STT(LLG):
+
+    def __init__(self, mesh, name='unnamed'):
+        """Simulation object.
+
+        *Arguments*
+
+          name : the Simulation name (used for writing data files, for examples)
+
+        """
+        super(LLG_STT, self).__init__(mesh, name=name)
+
+        self.field_stt = np.zeros(3 * self.n)
+
+        self._jx = np.zeros(self.n, dtype=np.float)
+        self._jy = np.zeros(self.n, dtype=np.float)
+        self._jz = np.zeros(self.n, dtype=np.float)
+
+        self.p = 0.5
+        self.beta = 0
+        self.update_j_fun = None
+
+        # FIXME: change the u0 to spatial
+        v = self.mesh.dx* self.mesh.dy * self.mesh.dz * (self.mesh.unit_length**3)
+        self.u0 = const.g_e * const.mu_B / (2 * const.c_e) * v
+
+    def get_jx(self):
+        return self._jx
+
+    def set_jx(self, value):
+        self._jx[:] = helper.init_scalar(value, self.mesh)
+
+    jx = property(get_jx, set_jx)
+
+    def get_jy(self):
+        return self._jy
+
+    def set_jy(self, value):
+        self._jy[:] = helper.init_scalar(value, self.mesh)
+
+    jy = property(get_jy, set_jy)
+
+    def get_jz(self):
+        return self._jz
+
+    def set_jz(self, value):
+        self._jz[:] = helper.init_scalar(value, self.mesh)
+
+    jz = property(get_jz, set_jz)
+
+    def sundials_rhs(self, t, y, ydot):
+
+        self.t = t
+
+        # already synchronized when call this funciton
+        # self.spin[:]=y[:]
+
+        self.compute_effective_field(t)
+
+        if self.update_j_fun is not None:
+            clib.compute_stt_field(self.spin,
+                                   self.field_stt,
+                                   self._jx * self.update_j_fun(t),
+                                   self._jy * self.update_j_fun(t),
+                                   self._jz * self.update_j_fun(t),
+                                   self.mesh.dx * self.mesh.unit_length,
+                                   self.mesh.dy * self.mesh.unit_length,
+                                   self.mesh.dz * self.mesh.unit_length,
+                                   self.mesh.neighbours,
+                                   self.n
+                                   )
+        else:
+            clib.compute_stt_field(self.spin,
+                                   self.field_stt,
+                                   self._jx,
+                                   self._jy,
+                                   self._jz,
+                                   self.mesh.dx * self.mesh.unit_length,
+                                   self.mesh.dy * self.mesh.unit_length,
+                                   self.mesh.dz * self.mesh.unit_length,
+                                   self.mesh.neighbours,
+                                   self.n
+                                   )
+
+        clib.compute_llg_stt_rhs(ydot,
+                                 self.spin,
+                                 self.field,
+                                 self.field_stt,
+                                 self.alpha,
+                                 self.beta,
+                                 self.u0 * self.p / self.mu_s_const,
+                                 self.gamma,
+                                 self.n)
diff --git a/fidimag/atomistic/llg_stt_cpp.py b/fidimag/atomistic/llg_stt_cpp.py
@@ -9,7 +9,7 @@
 import fidimag.common.helper as helper
 
 
-class LLG_STT_Slonczewski(LLG):
+class LLG_STT_CPP(LLG):
 
     def __init__(self, mesh, name='unnamed'):
         """Simulation object.
@@ -19,7 +19,7 @@ def __init__(self, mesh, name='unnamed'):
           name : the Simulation name (used for writing data files, for examples)
 
         """
-        super(LLG_STT_Slonczewski, self).__init__(mesh, name=name)
+        super(LLG_STT_CPP, self).__init__(mesh, name=name)
 
         self._p = np.zeros(3 * self.n, dtype=np.float)
         self._a_J = np.zeros(self.n, dtype=np.float)
diff --git a/fidimag/atomistic/sim.py b/fidimag/atomistic/sim.py
@@ -1,11 +1,13 @@
 from .llg import LLG
 from .sllg import SLLG
-from .llg_stt_slonczewski_type import LLG_STT_Slonczewski
+from .llg_stt import LLG_STT
+from .llg_stt_cpp import LLG_STT_CPP
 
 
 KNOWN_DRIVERS = {'llg': LLG,
                  'sllg': SLLG,
-                 'slonczewski': LLG_STT_Slonczewski
+                 'llg_stt': LLG_STT,
+                 'llg_stt_cpp': LLG_STT_CPP
 		}
 
 
diff --git a/fidimag/micro/llg_stt.py b/fidimag/micro/llg_stt.py
diff --git a/tests/test_stt.py b/tests/test_stt.py
diff --git a/tests/test_stt_dw_atomistic.py b/tests/test_stt_dw_atomistic.py
diff --git a/tests/test_stt_slonczewski.py b/tests/test_stt_slonczewski.py
