clean the cython version of sllg class

Author: Weiwei Wang

fidimag/atomistic/lib/clib.h

@@ -84,40 +84,5 @@ void compute_px_py_c(double *spin, int nx, int ny, int nz,
 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);
 
-//=========================================================
-//=========================================================
-//used for sode
-typedef struct {
-	int n;
-
-	double dt;
-	double T;
-	double gamma;
-	double *mu_s;
-	double coeff;
-	double Q;
-
-	double theta;
-	double theta1;
-	double theta2;
-
-	double *dm1;
-	double *dm2;
-	double *eta;
-
-} ode_solver;
-
-void init_solver(ode_solver *s, double k_B, double theta,
-                 int n, double dt, double gamma);
-
-ode_solver *create_ode_plan(void);
-
-void finalize_ode_plan(ode_solver *plan);
-
-void run_step1(ode_solver *s, double *m, double *h, double *m_pred, double *T,
-		double *alpha, double *mu_s_inv, int *pins);
-
-void run_step2(ode_solver *s, double *m_pred, double *h, double *m, double *T,
-		double *alpha, double *mu_s_inv, int *pins);
 
 #endif

fidimag/atomistic/lib/clib.pyx

@@ -71,24 +71,6 @@ cdef extern from "clib.h":
 
     # 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
-
-    ode_solver *create_ode_plan()
-
-    void init_solver(ode_solver *s, double k_B, double theta,
-                     int n, double dt, double gamma)
-
-    void finalize_ode_plan(ode_solver *plan)
-
-    void run_step1(ode_solver *s, double *m, double *h,
-                   double *m_pred, double *T, double *alpha,
-                   double *mu_s_inv, int *pins)
-
-    void run_step2(ode_solver *s, double *m_pred,
-                   double *h, double *m, double *T,
-                   double *alpha, double *mu_s_inv, int *pins)
 
 
 def compute_skyrmion_number(np.ndarray[double, ndim=1, mode="c"] spin,
@@ -285,89 +267,3 @@ def compute_llg_rhs_dw(np.ndarray[double, ndim=1, mode="c"] dm,
                 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
-    cdef double dt
-    cdef update_fun
-    cdef np.ndarray pred_m
-    cdef np.ndarray field
-    cdef np.ndarray mu_s_inv
-    cdef np.ndarray T
-    cdef np.ndarray alpha
-    cdef np.ndarray pins
-
-    cdef public double t
-    cdef public int step
-    cdef public np.ndarray y
-
-    def __cinit__(self,dt,n,gamma,k_B,theta,
-                            np.ndarray[double, ndim=1, mode="c"] mu_s_inv,
-                            np.ndarray[double, ndim=1, mode="c"] alpha,
-                            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[int, ndim=1, mode="c"] pins,
-                            update_fun):
-
-        self.t = 0
-        self.step = 0
-        self.dt = dt
-
-        self.update_fun = update_fun
-        self.mu_s_inv = mu_s_inv
-        self.field = field
-        self.T = T
-        self.alpha = alpha
-        self.pins = pins
-        self.pred_m = numpy.zeros(3*n,dtype=numpy.float)
-        self.y = numpy.zeros(3*n,dtype=numpy.float)
-
-
-        self._c_plan = create_ode_plan()
-        if self._c_plan is NULL:
-            raise MemoryError()
-
-        init_solver(self._c_plan,k_B,theta,n,dt,gamma)
-
-    def __dealloc__(self):
-        if self._c_plan is not NULL:
-            finalize_ode_plan(self._c_plan)
-            self._c_plan = NULL
-    
-    def set_initial_value(self,np.ndarray[double, ndim=1, mode="c"] spin, t):
-        self.t = t
-        self.y[:] = spin[:]
-
-    def successful(self):
-        #print self.spin
-        return True
-
-    def run_step(self):
-        cdef np.ndarray[double, ndim=1, mode="c"] y=self.y
-        cdef np.ndarray[double, ndim=1, mode="c"] field=self.field
-        cdef np.ndarray[double, ndim=1, mode="c"] pred_m=self.pred_m
-        cdef np.ndarray[double, ndim=1, mode="c"] T=self.T
-        cdef np.ndarray[double, ndim=1, mode="c"] alpha=self.alpha
-        cdef np.ndarray[double, ndim=1, mode="c"] mu_s_inv=self.mu_s_inv
-        cdef np.ndarray[int, ndim=1, mode="c"] pins = self.pins
-
-        #print "from cython1", self.spin,self.field,self.pred_m
-        self.update_fun(self.y, self.t)
-        run_step1(self._c_plan, &y[0], &field[0], &pred_m[0],
-                  &T[0], &alpha[0], &mu_s_inv[0], &pins[0])
-
-        self.step += 1
-        self.t = self.step * self.dt
-
-        self.update_fun(self.pred_m, self.t)
-        run_step2(self._c_plan, &pred_m[0], &field[0],
-                  &y[0], &T[0], &alpha[0], &mu_s_inv[0], &pins[0])
-
-
-    def run_until(self, t):
-        while (self.t<t):
-            print "from cython %g"%self.t
-            self.run_step()
-        return 0

fidimag/atomistic/lib/sllg.c

@@ -42,150 +42,6 @@ void llg_rhs_dw_c(double *m, double *h, double *dm, double *T, double *alpha, do
 	}
 }
 
-
-
-ode_solver *create_ode_plan(void) {
-
-	ode_solver *plan = (ode_solver*) malloc(sizeof(ode_solver));
-
-	return plan;
-}
-
-
-void init_solver(ode_solver *s, double k_B, double theta, int n, double dt, double gamma) {
-
-	s->theta = theta;
-	s->theta1 = 1-0.5/theta;
-	s->theta2 = 0.5/theta;
-
-	s->dt = dt;
-	s->n = n;
-	s->gamma = gamma;
-	
-	s->Q = 2 * k_B * dt / gamma;
-
-	s->dm1 = (double*) malloc(3 * n * sizeof(double));
-	s->dm2 = (double*) malloc(3 * n * sizeof(double));
-	s->eta = (double*) malloc(3 * n * sizeof(double));
-
-	int i = 0;
-	for (i = 0; i < 3 * n; i++) {
-		s->dm1[i] = 0;
-		s->dm2[i] = 0;
-		s->eta[i] = 0;
-	}
-
-	initial_random(2);
-}
-
-
-void llg_rhs_dw(ode_solver *s, double *m, double *h, double *dm, double *T, double *alpha, double *mu_s_inv, int *pins) {
-
-	int i, j, k;
-
-	double mth0, mth1, mth2;
-	
-	int n = s->n;
-	double *eta = &s->eta[0];
-	double dt = s->dt;
-	double Q=s->Q;
-
-	double coeff, q;
-	double hi,hj,hk;
-
-	#pragma omp parallel for private(i,j,k,q,coeff,hi,hj,hk,mth0,mth1,mth2)
-	for (int id = 0; id < n; id++) {
-		i = 3*id;
-		j = i+1;
-		k = j+1;
-		
-		if (pins[id]>0){
-			 dm[i] = 0;
-			 dm[j] = 0;
-			 dm[k] = 0;
-			 continue;
-		}
-
-
-        coeff = -s->gamma/ (1.0 + alpha[id] * alpha[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;
-		hk = h[k]*dt + eta[k]*q;
-		
-		mth0 = coeff * (m[j] * hk - m[k] * hj);
-		mth1 = coeff * (m[k] * hi - m[i] * hk);
-		mth2 = coeff * (m[i] * hj - m[j] * hi);
-
-
-		dm[i] = mth0 + alpha[i] * (m[j] * mth2 - m[k] * mth1);
-		dm[j] = mth1 + alpha[i] * (m[k] * mth0 - m[i] * mth2);
-		dm[k] = mth2 + alpha[i] * (m[i] * mth1 - m[j] * mth0);
-		
-	}
-}
-
-void run_step1(ode_solver *s, double *m, double *h, double *m_pred, double *T, double *alpha, double *mu_s_inv, int *pins) {
-	int i;
-	int n = s->n;
-	double *dm1 = s->dm1;
-	double theta = s->theta;
-
-	gauss_random_vec(s->eta, 3 * n);
-	llg_rhs_dw(s, m, h, dm1, T, alpha, mu_s_inv, pins);
-
-	//#pragma omp parallel for private(i)
-	for (i = 0; i < 3 * s->n; i++) {
-		m_pred[i] = m[i] + theta * dm1[i];
-	}
-
-}
-
-void run_step2(ode_solver *s, double *m_pred, double *h, double *m, double *T, double *alpha, double *mu_s_inv, int *pins) {
-	int i, j, k;
-	int n = s->n;
-	double *dm1 = s->dm1;
-	double *dm2 = s->dm2;
-	double theta1 = s->theta1;
-	double theta2 = s->theta2;
-
-	llg_rhs_dw(s, m_pred, h, dm2, T, alpha, mu_s_inv, pins);
-
-	//#pragma omp parallel for private(i)
-	for (i = 0; i < 3 * n; i++) {
-		m[i] += (theta1 * dm1[i] + theta2 * dm2[i]);
-	}
-
-	double mm;
-	#pragma omp parallel for private(i,j,k,mm)
-	for (int id = 0; id < s->n; id++) {
-
-		if (pins[id]>0){
-			continue;
-		}
-
-		i= 3*id;
-		j = i + 1;
-		k = j + 1;
-
-		mm = 1.0 / sqrt(m[i] * m[i] + m[j] * m[j] + m[k] * m[k]);
-		m[i] *= mm;
-		m[j] *= mm;
-		m[k] *= mm;
-	}
-	//we have to say that this kind of method is quite inaccurate, so in future we can try other methods
-
-}
-
-void finalize_ode_plan(ode_solver *plan) {
-	free(plan->dm1);
-	free(plan->dm2);
-	free(plan->eta);
-	free(plan);
-}
-
-
 void normalise(double *m, int n){
 	int i, j, k;
 	double mm;