make guiding center more flexible

Author: Weiwei Wang

examples/micromagnetic/dw_stt_cpp/main.py

@@ -120,7 +120,7 @@ def excite_system(mesh, beta=0.0):
     # The simulation will run for 5 ns and save
     # 500 snapshots of the system in the process
     ts = np.linspace(0, 0.5e-9, 21)
-    
+
     xs=[]
     thetas=[]
 
@@ -154,15 +154,14 @@ def plot():
     z_an, theta_an = analytical(ts)
 
     fig = plt.figure(figsize=(5, 2.6))
-    
+
     gs = gridspec.GridSpec(1, 2, width_ratios=[4, 4])
     ax0 = plt.subplot(gs[0])
     ax0.plot(ts*1e9,xs,'.')
     ax0.plot(ts*1e9, z_an, '-')
     plt.ylabel(r'DW shift')
     plt.xlabel(r'Time (ns)')
 
-
     ax1 = plt.subplot(gs[1])
     ax1.plot(ts*1e9,theta,'.')
     ax1.plot(ts*1e9, theta_an, '-')
@@ -173,7 +172,6 @@ def plot():
     analytical(ts)
 
     plt.savefig('xs_theta.pdf')
-    
 
 if __name__ == '__main__':
 
@@ -189,4 +187,3 @@ def plot():
 
     excite_system(mesh, beta=0.16)
     plot()
-

fidimag/atomistic/lib/clib.h

@@ -1,97 +1,101 @@
 #ifndef __CLIB__
 #define __CLIB__
 
-#include<math.h>
-#include<complex.h>
-#include<fftw3.h>
+#include <complex.h>
+#include <fftw3.h>
+#include <math.h>
 //#include<omp.h>
 
-#include"fidimag_random.h"
+#include "fidimag_random.h"
 
 #define WIDE_PI 3.1415926535897932384626433832795L
 
 /* 3 components for the cross product calculations */
-inline double cross_x(double a0, double a1, double a2,
-                      double b0, double b1, double b2) { return a1 * b2 - a2 * b1; }
-inline double cross_y(double a0, double a1, double a2,
-                      double b0, double b1, double b2) { return a2 * b0 - a0 * b2; }
-inline double cross_z(double a0, double a1, double a2,
-                      double b0, double b1, double b2) { return a0 * b1 - a1 * b0; }
-
-void compute_exch_field(double *spin, double *field, double *energy,
-						double Jx, double Jy, double Jz,
-                        int *ngbs, int n);
-
-
-void compute_exch_field_spatial(double *spin, double *field, double *energy,double *J,
-                        int *ngbs, int n);
+inline double cross_x(double a0, double a1, double a2, double b0, double b1,
+                      double b2) {
+  return a1 * b2 - a2 * b1;
+}
+inline double cross_y(double a0, double a1, double a2, double b0, double b1,
+                      double b2) {
+  return a2 * b0 - a0 * b2;
+}
+inline double cross_z(double a0, double a1, double a2, double b0, double b1,
+                      double b2) {
+  return a0 * b1 - a1 * b0;
+}
+
+void compute_exch_field(double *spin, double *field, double *energy, double Jx,
+                        double Jy, double Jz, int *ngbs, int n);
+
+void compute_exch_field_spatial(double *spin, double *field, double *energy,
+                                double *J, int *ngbs, int n);
 
 double compute_exch_energy(double *spin, double Jx, double Jy, double Jz,
-                           int nx, int ny, int nz,
-                           int xperiodic, int yperiodic);
+                           int nx, int ny, int nz, int xperiodic,
+                           int yperiodic);
 
-void compute_anis(double *spin, double *field, double *energy,
-	              double *Ku, double *axis, int n);
+void compute_anis(double *spin, double *field, double *energy, double *Ku,
+                  double *axis, int n);
 
-void dmi_field_bulk(double *spin, double *field, double *energy,
-                    double *D, int *ngbs, int n);
+void dmi_field_bulk(double *spin, double *field, double *energy, double *D,
+                    int *ngbs, int n);
 
 void dmi_field_interfacial_atomistic(double *spin, double *field,
-                                     double *energy, double D, int *ngbs,
-                                     int n, int nneighbours,
-                                     double *DMI_vec);
+                                     double *energy, double D, int *ngbs, int n,
+                                     int nneighbours, double *DMI_vec);
 
-void demag_full(double *spin, double *field, double *energy,
-                double *coords,
+void demag_full(double *spin, double *field, double *energy, double *coords,
                 double *mu_s, double *mu_s_scale, int n);
 
-double dmi_energy(double *spin, double D, int nx, int ny, int nz,
-                  int xperiodic, int yperiodic);
+double dmi_energy(double *spin, double D, int nx, int ny, int nz, int xperiodic,
+                  int yperiodic);
 
-void llg_rhs(double * dm_dt, double * spin, double * h, double *alpha,
-		     int *pins, double gamma, int n, int do_precession,
-             double default_c);
+void llg_rhs(double *dm_dt, double *spin, double *h, double *alpha, int *pins,
+             double gamma, int n, int do_precession, double default_c);
 
-void llg_rhs_jtimes(double *jtn, double *m, double *h,
-                    double *mp, double *hp, double *alpha, int *pins,
-                    double gamma, int n, int do_precession, double default_c);
+void llg_rhs_jtimes(double *jtn, double *m, double *h, double *mp, double *hp,
+                    double *alpha, int *pins, double gamma, int n,
+                    int do_precession, double default_c);
 
-void llg_s_rhs(double * dm_dt, double * spin, double * h, double *alpha,
-             double *chi, double gamma, int n);
+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 *jz, double dx, double dy, double dz, 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,
-                 double beta, double u0, double gamma, int n);
-
-void llg_stt_cpp(double *dm_dt, double *m, double *h, double *p,
-			      double *alpha, int *pins, double *a_J, double beta, double gamma, 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);
 
+void llg_stt_cpp(double *dm_dt, double *m, double *h, double *p, double *alpha,
+                 int *pins, double *a_J, double beta, double gamma, int n);
 
 void normalise(double *m, int *pins, int n);
 
-double skyrmion_number(double *spin, double *charge,
-                       int nx, int ny, int nz, int *ngbs);
+double skyrmion_number(double *spin, double *charge, int nx, int ny, int nz,
+                       int *ngbs);
 
-double skyrmion_number_BergLuscher(double *spin, double *charge,
-                                   int nx, int ny, int nz, int *ngbs);
+double skyrmion_number_BergLuscher(double *spin, double *charge, int nx, int ny,
+                                   int nz, int *ngbs);
 
-void compute_guiding_center(double *spin, int nx, int ny, int nz, double *res);
+void compute_guiding_center(double *spin, int nx, int ny, int nz, int nx_start,
+                            int nx_stop, int ny_start, int ny_stop,
+                            double *res);
 
-void compute_px_py_c(double *spin, int nx, int ny, int nz,
-                     double *px, double *py);
+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);
+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);
 
 //======================================================================
 
-void run_step_mc(mt19937_state *state, double *spin, double *new_spin, int *ngbs, int *nngbs, double J, double J1, double D, double D1, double *h, double Kc, int n, double T, int hexagnoal_mesh);
-
+void run_step_mc(mt19937_state *state, double *spin, double *new_spin,
+                 int *ngbs, int *nngbs, double J, double J1, double D,
+                 double D1, double *h, double Kc, int n, double T,
+                 int hexagnoal_mesh);
 
 #endif

fidimag/atomistic/lib/clib.pyx

@@ -5,7 +5,7 @@ np.import_array()
 cdef extern from "fidimag_random.h":
     ctypedef struct mt19937_state:
         pass
-    
+
     mt19937_state *create_mt19937_state()
     void finalize_mt19937_state(mt19937_state *state)
 
@@ -19,7 +19,7 @@ cdef extern from "time.h":
     time_t time(time_t *timer)
 
 cdef extern from "clib.h":
-    void run_step_mc(mt19937_state *state, double *spin, double *new_spin, int *ngbs, int *nngbs, 
+    void run_step_mc(mt19937_state *state, double *spin, double *new_spin, int *ngbs, int *nngbs,
                     double J, double J1, double D, double D1, double *h, double Kc, int n, double T, int hexagnoal_mesh)
     double skyrmion_number(double *spin, double *charge,
                            int nx, int ny, int nz, int *ngbs)
@@ -28,7 +28,8 @@ cdef extern from "clib.h":
     double skyrmion_number_BergLuscher(double *spin, double *charge,
                                        int nx, int ny, int nz, int *ngbs)
 
-    void compute_guiding_center(double *spin, int nx, int ny, int nz,
+    void compute_guiding_center(double *spin, int nx, int ny, int nz, int nx_start,
+                                int nx_stop, int ny_start, int ny_stop,
                                 double *res)
 
     void compute_px_py_c(double *spin, int nx, int ny, int nz,
@@ -107,13 +108,17 @@ def compute_skyrmion_number_BergLuscher(np.ndarray[double, ndim=1, mode="c"] spi
     return skyrmion_number_BergLuscher(&spin[0], &charge[0], nx, ny, nz, &ngbs[0,0])
 
 def compute_RxRy(np.ndarray[double, ndim=1, mode="c"] spin,
-                            nx, ny, nz):
+                            nx, ny, nz, nx_start=0, nx_stop=-1, ny_start=0, ny_stop=-1):
 
     res = numpy.array([0.0,0.0])
+    if nx_stop < 0 or nx_stop > nx:
+        nx_stop = nx
+    if ny_stop < 0 or ny_stop > ny:
+        ny_stop = ny
 
     cdef np.ndarray[double, ndim=1, mode="c"] R = res
 
-    compute_guiding_center(&spin[0], nx, ny, nz, &R[0])
+    compute_guiding_center(&spin[0], nx, ny, nz, nx_start, nx_stop, ny_start, ny_stop, &R[0])
 
     return res[0], res[1]
 
@@ -258,7 +263,7 @@ def compute_llg_stt_cpp(np.ndarray[double, ndim=1, mode="c"] dm_dt,
 		            np.ndarray[int, ndim=1, mode="c"] pin,
                 np.ndarray[double, ndim=1, mode="c"] a_J,
                 beta, gamma, n):
-    llg_stt_cpp(&dm_dt[0], &spin[0], &field[0], &p[0], 
+    llg_stt_cpp(&dm_dt[0], &spin[0], &field[0], &p[0],
                 &alpha[0], &pin[0], &a_J[0], beta, gamma, n)
 
 
@@ -286,11 +291,11 @@ cdef class rng_mt19937(object):
             self.seed = int(seed)
         else:
             self.seed = time(NULL)
-                
+
         self._c_state = create_mt19937_state()
         if self._c_state is NULL:
             raise MemoryError()
-        
+
         initial_rng_mt19973(self._c_state, self.seed)
 
     def set_seed(self, seed):
@@ -302,7 +307,7 @@ cdef class rng_mt19937(object):
             return a random number in [0,1)
         """
         return random_double_half_open(self._c_state)
-    
+
     def fill_vector_gaussian(self, np.ndarray[np.float64_t, ndim=1] vector):
         gauss_random_vector(self._c_state, &vector[0], vector.shape[0])
 
@@ -315,25 +320,24 @@ cdef class rng_mt19937(object):
             self._c_state = NULL
 
 cdef class monte_carlo(object):
-    cdef mt19937_state *_c_state            
+    cdef mt19937_state *_c_state
 
     def __init__(self, seed=-1):
-                
+
         self._c_state = create_mt19937_state()
         if self._c_state is NULL:
             raise MemoryError()
-        
+
         initial_rng_mt19973(self._c_state, seed)
 
     def set_seed(self, seed):
         initial_rng_mt19973(self._c_state, seed)
- 
+
     def run_step(self,np.ndarray[double, ndim=1, mode="c"] spin,
                 np.ndarray[double, ndim=1, mode="c"] new_spin,
                 np.ndarray[int, ndim=2, mode="c"] ngbs,
                 np.ndarray[int, ndim=2, mode="c"] nngbs,
-                J, J1, D, D1, np.ndarray[double, ndim=1, mode="c"] h, 
+                J, J1, D, D1, np.ndarray[double, ndim=1, mode="c"] h,
                 Kc, n, T, hexagnoal_mesh):
 
         run_step_mc(self._c_state, &spin[0], &new_spin[0], &ngbs[0,0], &nngbs[0,0], J, J1, D, D1, &h[0], Kc, n, T, hexagnoal_mesh)
-