try to move cython version to python version

Author: Weiwei Wang

fidimag/atomistic/lib/clib.h

@@ -79,6 +79,10 @@ void compute_guiding_center(double *spin, int nx, int ny, int nz, double *res);
 void compute_px_py_c(double *spin, int nx, int ny, int nz,
                      double *px, double *py);
 
+//======================================================================
+
+void llg_rhs_dw_c(double *m, double *h, double *dm, double *T, double *alpha, 
+                  double *mu_s_inv, int *pins, double *eta, int n, double gamma, double dt);
 
 //=========================================================
 //=========================================================

fidimag/atomistic/lib/clib.pyx

@@ -4,8 +4,8 @@ np.import_array()
 
 
 cdef extern from "clib.h":
-
-    void gauss_random_vec_with_init(double *x, int n)
+    void initial_random(int seed)
+    void gauss_random_vec(double *x, int n)
 
     double skyrmion_number(double *spin, double *charge,
                            int nx, int ny, int nz, int *ngbs)
@@ -70,6 +70,8 @@ cdef extern from "clib.h":
 			      double *alpha, int *pins, double *a_J, double beta, double gamma, int n)
 
     # used for sllg
+    void llg_rhs_dw_c(double *m, double *h, double *dm, double *T, double *alpha, double *mu_s_inv, int *pins, double *eta, int n, double gamma, double dt)
+    
     ctypedef struct ode_solver:
         pass
 
@@ -88,15 +90,6 @@ cdef extern from "clib.h":
                    double *h, double *m, double *T,
                    double *alpha, double *mu_s_inv, int *pins)
 
-def random_number(np.ndarray[double, ndim=1, mode="c"] v):
-    cdef int n = len(v)
-
-
-    print n
-
-    gauss_random_vec_with_init(&v[0], n)
-
-
 
 def compute_skyrmion_number(np.ndarray[double, ndim=1, mode="c"] spin,
                             np.ndarray[double, ndim=1, mode="c"] charge,
@@ -273,6 +266,26 @@ def normalise_spin(np.ndarray[double, ndim=1, mode="c"] spin, n):
     normalise(&spin[0], n)
 
 
+def init_random(seed):
+    initial_random(seed);
+
+def random_number_array(np.ndarray[double, ndim=1, mode="c"] v):
+    
+    cdef int n = len(v)
+
+    gauss_random_vec(&v[0], n)
+
+
+def compute_llg_rhs_dw(np.ndarray[double, ndim=1, mode="c"] dm,
+                np.ndarray[double, ndim=1, mode="c"] spin,
+                np.ndarray[double, ndim=1, mode="c"] field,
+                np.ndarray[double, ndim=1, mode="c"] T,
+                np.ndarray[double, ndim=1, mode="c"] alpha,
+                np.ndarray[double, ndim=1, mode="c"] mu_s_inv,
+                np.ndarray[double, ndim=1, mode="c"] eta,
+                np.ndarray[int, ndim=1, mode="c"] pin, n, gamma, dt):
+    llg_rhs_dw_c(&spin[0], &field[0], &dm[0], &T[0], &alpha[0],  &mu_s_inv[0], &pin[0], &eta[0], n, gamma, dt)
+
 
 cdef class RK2S(object):
     cdef ode_solver * _c_plan

fidimag/atomistic/lib/llg_random.h

@@ -43,8 +43,8 @@ static inline unsigned int int_rand(void) {
 }
 
 //temporary random generator
-void initial_random(void) {
-	unsigned int seed = (unsigned int) time(NULL);
+void initial_random(int seed) {
+	//unsigned int seed = (unsigned int) time(NULL);
 
 	int i;
 	MT[0] = seed & 0xFFFFFFFFU;
@@ -136,12 +136,4 @@ void gauss_random_vec(double *x, int n) {
 
 }
 
-void gauss_random_vec_with_init(double *x, int n) {
-
-	initial_random();
-
-	gauss_random_vec(x,n);
-
-}
-
 #endif

fidimag/atomistic/lib/sllg.c

@@ -1,6 +1,48 @@
 #include "clib.h"
 #include "llg_random.h"
 
+/*
+* n is the spin number
+* eta is the random number array
+*/
+void llg_rhs_dw_c(double *m, double *h, double *dm, double *T, double *alpha, double *mu_s_inv, int *pins, double *eta, int n, double gamma, double dt) {
+        
+        double k_B = 1.3806505e-23;
+        double Q = 2 * k_B * dt / gamma;
+
+	//#pragma omp parallel for
+	for (int id = 0; id < n; id++) {
+		int i = 3*id;
+		int j = i+1;
+		int k = j+1;
+		
+		if (pins[id]>0){
+			 dm[i] = 0;
+			 dm[j] = 0;
+			 dm[k] = 0;
+			 continue;
+		}
+
+
+                double coeff = -gamma/ (1.0 + alpha[id] * alpha[id]);
+		double q = sqrt(Q * alpha[id] * T[id] * mu_s_inv[id]);
+		
+		double hi = h[i]*dt + eta[i]*q;
+		double hj = h[j]*dt + eta[j]*q;
+		double hk = h[k]*dt + eta[k]*q;
+		
+		double mth0 = coeff * (m[j] * hk - m[k] * hj);
+		double mth1 = coeff * (m[k] * hi - m[i] * hk);
+		double mth2 = coeff * (m[i] * hj - m[j] * hi);
+
+		dm[i] = mth0 + alpha[id] * (m[j] * mth2 - m[k] * mth1);
+		dm[j] = mth1 + alpha[id] * (m[k] * mth0 - m[i] * mth2);
+		dm[k] = mth2 + alpha[id] * (m[i] * mth1 - m[j] * mth0);
+		
+	}
+}
+
+
 
 ode_solver *create_ode_plan(void) {
 
@@ -33,7 +75,7 @@ void init_solver(ode_solver *s, double k_B, double theta, int n, double dt, doub
 		s->eta[i] = 0;
 	}
 
-	initial_random();
+	initial_random(2);
 }
 
 
@@ -57,7 +99,7 @@ void llg_rhs_dw(ode_solver *s, double *m, double *h, double *dm, double *T, doub
 		j = i+1;
 		k = j+1;
 
-		if (pins[i]>0){
+		if (pins[id]>0){
 			 dm[i] = 0;
 			 dm[j] = 0;
 			 dm[k] = 0;
@@ -66,7 +108,7 @@ void llg_rhs_dw(ode_solver *s, double *m, double *h, double *dm, double *T, doub
 
 
         coeff = -s->gamma/ (1.0 + alpha[id] * alpha[id]) ;
-		q = sqrt(Q * alpha[i] * T[i] * mu_s_inv[id]);
+		q = sqrt(Q * alpha[i] * T[id] * mu_s_inv[id]);
 
 		hi = h[i]*dt + eta[i]*q;
 		hj = h[j]*dt + eta[j]*q;

fidimag/atomistic/sllg.py

@@ -19,6 +19,10 @@ def __init__(self, mesh, name='unnamed'):
         super(SLLG, self).__init__(mesh, name=name)
 
         self._T = np.zeros(self.n, dtype=np.float)
+        self.next_spin = np.zeros(3*self.n, dtype=np.float)
+        self.eta = np.zeros(3*self.n, dtype=np.float)
+        self.dm1 = np.zeros(3*self.n, dtype=np.float)
+        self.dm2 = np.zeros(3*self.n, dtype=np.float)
 
         self.set_options()
 
@@ -30,30 +34,78 @@ def set_T(self, T0):
 
     T = property(get_T, set_T)
 
-    def set_options(self, dt=1e-15, theta=1.0, gamma=const.gamma, k_B=const.k_B):
+    def set_options(self, dt=1e-15, theta=1.0, gamma=const.gamma, k_B=const.k_B, seed=100):
 
+        clib.init_random(seed)
         self.gamma = gamma
         self.k_B = k_B
-        print 'self.n=',self.n
-
-        self.vode = clib.RK2S(dt,
-                              self.n,
-                              self.gamma,
-                              self.k_B,
-                              theta,
-                              self._mu_s_inv,
-                              self.alpha,
-                              self.spin,
-                              self.field,
-                              self._T,
-                              self._pins,
-                              self.update_effective_field)
+        self.dt = dt
+        self.theta = theta
+        self.theta1 = 1-0.5/theta;
+	self.theta2 = 0.5/theta;
+
+    def run_step(self):
+        
+        clib.random_number_array(self.eta)
+        
+        #step1
+        self.update_effective_field(self.spin, self.t)
+        clib.compute_llg_rhs_dw(self.dm1,
+				self.spin,
+				self.field,
+				self._T,
+				self._alpha,
+				self._mu_s_inv,
+				self.eta,
+				self._pins, 
+				self.n, 
+				self.gamma, 
+				self.dt)
+        self.next_spin = self.spin + self.theta * self.dm1
+
+        self.step += 1
+        self.t = self.dt*self.step
+
+        #step2
+        self.update_effective_field(self.next_spin, self.t)
+        clib.compute_llg_rhs_dw(self.dm2,
+				self.next_spin,
+				self.field,
+				self._T,
+				self._alpha,
+				self._mu_s_inv,
+				self.eta,
+				self._pins, 
+				self.n, 
+				self.gamma, 
+				self.dt)
+        self.spin += (self.theta1*self.dm1+self.theta2*self.dm2)
+
+        clib.normalise_spin(self.spin, self.n)
+        
 
-    def update_effective_field(self, y, t):
 
-        self.spin[:] = y[:]
+    def update_effective_field(self, y, t):
 
         self.field[:] = 0
-        return 0
+
         for obj in self.interactions:
-            self.field += obj.compute_field(t)
+            self.field += obj.compute_field(t, spin=y)
+
+
+    def run_until(self, t):
+
+        if t <= self.t:
+            if t == self.t and self.t == 0.0:
+                self.compute_effective_field(t)
+                self.saver.save()
+            return
+
+        self.spin_last[:] = self.spin[:]
+
+        while (self.t<t):
+            self.run_step()
+
+        # update field before saving data
+        self.compute_effective_field(t)
+        self.saver.save()

fidimag/atomistic/zeeman.py

@@ -30,7 +30,7 @@ def setup(self, mesh, spin, mu_s):
         self.field = np.zeros(3 * self.n)
         self.field[:] = helper.init_vector(self.H0, self.mesh)
 
-    def compute_field(self, t=0):
+    def compute_field(self, t=0, spin=None):
         return self.field
 
     def average_field(self):