Updated demag C lib to allow calc of energy density per mesh site

Author: davidcorteso

fidimag/atomistic/demag.py

@@ -67,6 +67,6 @@ def compute_exact(self):
     def compute_energy(self):
 
         energy = self.demag.compute_energy(
-            self.spin, self.mu_s_scale, self.field)
+            self.spin, self.mu_s_scale, self.field, self.energy)
 
         return energy / self.scale

fidimag/atomistic/demag_hexagonal.py

@@ -81,6 +81,9 @@ def setup(self, mesh, spin, mu_s, mu_s_inv):
         self.field = np.zeros(3 * self.n, dtype=np.float)
         self.field_c = np.zeros(3 * self.n_c, dtype=np.float)
 
+        self.energy = np.zeros(self.n, dtype=np.float)
+        self.energy_c = np.zeros(self.n_c, dtype=np.float)
+
         unit_length = mesh.unit_length
         self.mu_s_scale = np.zeros(mesh.n, dtype=np.float)
         self.mu_s_scale_c = np.zeros(2 * self.n, dtype=np.float)
@@ -144,7 +147,7 @@ def compute_energy(self):
         self.vector2cuboid(self.spin, self.spin_c)
 
         energy = self.demag.compute_energy(
-            self.spin_c, self.mu_s_scale_c, self.field_c)
+            self.spin_c, self.mu_s_scale_c, self.field_c, self.energy_c)
 
         return energy / self.scale
 

fidimag/common/dipolar/demag.c

@@ -417,22 +417,28 @@ void exact_compute(fft_demag_plan *restrict plan, double *restrict spin,  double
 
 }
 
-double compute_demag_energy(fft_demag_plan *restrict plan, double *restrict spin, double *restrict mu_s, double *restrict field) {
+double compute_demag_energy(fft_demag_plan *restrict plan,
+                            double *restrict spin,
+                            double *restrict mu_s,
+                            double *restrict field,
+                            double *restrict energy
+                            ) {
 
 	int i,j;
 
 	int nxyz = plan->nx * plan->ny * plan->nz;
 
-	double energy = 0;
+	double total_energy = 0;
 
 	for (i = 0; i < nxyz; i++) {
 		j = 3*i;
-		energy += mu_s[i]*(spin[j]*field[j]+spin[j+1]*field[j+1]+spin[j+2]*field[j+2]);
+		energy[i] = -0.5 * mu_s[i] * (spin[j]   * field[j]   + 
+                                      spin[j+1] * field[j+1] + 
+                                      spin[j+2] * field[j+2] );
+        total_energy += energy[i];
 	}
 
-	energy = -energy / 2.0;
-
-	return energy;
+	return total_energy;
 
 }
 

fidimag/common/dipolar/dipolar.h

@@ -78,7 +78,7 @@ void fill_demag_tensors_c(fft_demag_plan *restrict plan, double *restrict tensor
 
 void compute_fields(fft_demag_plan *restrict plan, double *restrict spin, double *mu_s, double *restrict field);
 void exact_compute(fft_demag_plan *restrict plan, double *restrict spin, double *mu_s, double *restrict field);
-double compute_demag_energy(fft_demag_plan *restrict plan, double *restrict spin, double *restrict mu_s, double *restrict field);
+double compute_demag_energy(fft_demag_plan *restrict plan, double *restrict spin, double *restrict mu_s, double *restrict field, double *restrict energy);
 
 
 

fidimag/common/dipolar/dipolar.pyx

@@ -41,7 +41,7 @@ cdef extern from "dipolar.h":
     void init_plan(fft_demag_plan * plan, double dx, double dy, double dz, int nx,int ny, int nz)
     void compute_fields(fft_demag_plan * plan, double *spin, double *mu_s, double *field)
     void exact_compute(fft_demag_plan * plan, double *spin, double *mu_s, double *field)
-    double compute_demag_energy(fft_demag_plan *plan, double *spin, double *mu_s, double *field)
+    double compute_demag_energy(fft_demag_plan *plan, double *spin, double *mu_s, double *field, double *energy)
     void compute_dipolar_tensors(fft_demag_plan *plan)
     void compute_demag_tensors(fft_demag_plan *plan)
     void create_fftw_plan(fft_demag_plan *plan)
@@ -177,9 +177,11 @@ cdef class FFTDemag(object):
     def compute_energy(self,
                       np.ndarray[double, ndim=1, mode="c"] spin,
                       np.ndarray[double, ndim=1, mode="c"] mu_s,
-                      np.ndarray[double, ndim=1, mode="c"] field):
+                      np.ndarray[double, ndim=1, mode="c"] field,
+                      np.ndarray[double, ndim=1, mode="c"] energy):
 
-        return compute_demag_energy(self._c_plan, &spin[0], &mu_s[0], &field[0])
+        return compute_demag_energy(self._c_plan, &spin[0], &mu_s[0],
+                                    &field[0], &energy[0])
 
 
 cdef extern from "demagcoef.h":

fidimag/micro/demag.py

@@ -105,7 +105,13 @@ def compute_exact(self):
     def compute_energy(self):
 
         self.compute_field()
-        energy = self.demag.compute_energy(self.spin, self.Ms, self.field)
+        energy = self.demag.compute_energy(self.spin, self.Ms,
+                                           self.field, self.energy)
+
+        self.energy *= mu_0 * (self.mesh.dx *
+                               self.mesh.dy *
+                               self.mesh.dz *
+                               self.mesh.unit_length ** 3.)
 
         return energy * mu_0 * (self.mesh.dx *
                                 self.mesh.dy *