Merge pull request #70 from computationalmodelling/nebm_spherical

Nebm spherical

Author: rpep

fidimag/common/neb_spherical.py

@@ -365,12 +365,11 @@ def init_m(pos):
         self.coords = np.zeros(2 * self.n * self.total_image_num)
         self.last_m = np.zeros(self.coords.shape)
 
-        # For the effective field we use the extremes (to fit the energies
-        # array length). We could save some memory if we don't consider them
-        # (CHECK this in the future)
-        self.Heff = np.zeros(self.coords.shape)
+        # For the effective field we do not consider the extremes
+        # since they are fixed
+        self.Heff = np.zeros(2 * self.n * self.image_num)
         # self.Heff = np.zeros(2 * self.n * self.total_image_num)
-        self.Heff.shape = (self.total_image_num, -1)
+        self.Heff.shape = (self.image_num, -1)
 
         # Tangent components in spherical coordinates
         self.tangents = np.zeros(2 * self.n * self.image_num)
@@ -548,8 +547,22 @@ def save_vtks(self):
 
         self.coords.shape = (self.total_image_num, -1)
 
+        # We use Ms from the simulation assuming that all the
+        # images are the same
         for i in range(self.total_image_num):
-            self.vtk.save_vtk(spherical2cartesian(self.coords[i]), step=i, vtkname='m')
+            # We will try to save for the micromagnetic simulation (Ms)
+            # or an atomistic simulation (mu_s)
+            # TODO: maybe this can be done with an: isinstance
+            try:
+                self.vtk.save_vtk(spherical2cartesian(self.coords[i]).reshape(-1, 3),
+                                  self.sim.Ms,
+                                  step=i,
+                                  vtkname='m')
+            except:
+                self.vtk.save_vtk(spherical2cartesian(self.coords[i]).reshape(-1, 3),
+                                  self.sim._mu_s,
+                                  step=i,
+                                  vtkname='m')
 
         self.coords.shape = (-1, )
 
@@ -606,7 +619,7 @@ def compute_effective_field(self, y):
 
             # Set the simulation magnetisation to the (i+1)-th image
             # spin components
-            self.sim.spin = spherical2cartesian(y[i + 1])
+            self.sim.spin[:] = spherical2cartesian(y[i + 1])
 
             # Compute the effective field using Fidimag's methods.
             # (we use the time=0 since we are only using the simulation
@@ -618,7 +631,7 @@ def compute_effective_field(self, y):
             h = self.sim.field
             # Save the field components of the (i + 1)-th image
             # to the effective field array which is in spherical coords
-            self.Heff[i + 1, :] = cartesian2spherical_field(h, y[i + 1])
+            self.Heff[i, :] = cartesian2spherical_field(h, y[i + 1])
             # Compute and save the total energy for this image
             self.energy[i + 1] = self.sim.compute_energy()
 
@@ -674,7 +687,7 @@ def sundials_rhs(self, time, y, ydot):
         #   D_climb = -nabla E + 2 * [nabla E * t] t
         for i in range(self.image_num):
             # The eff field has (i + 1) since it considers the extreme images
-            h = self.Heff[i + 1]
+            h = self.Heff[i]
             # Tangents and springs has length: total images -2 (no extremes)
             t = self.tangents[i]
             sf = self.springs[i]

fidimag/common/vtk.py

@@ -51,6 +51,7 @@ def init_VtkData(self, structure, header):
     def reset_data(self):
         self.vtk_data = pyvtk.VtkData(self.structure, self.header)
 
+    @ignore_stderr
     def save_scalar(self, s, name="my_field", step=0):
         self.vtk_data.cell_data.append(pyvtk.Scalars(s, name))
 

tests/test_two_particles.py tests/test_two_particles_nebm.pytests/test_two_particles.py renamed to tests/test_two_particles_nebm.py

@@ -4,8 +4,9 @@
 # FIDIMAG:
 from fidimag.micro import Sim
 from fidimag.common import CuboidMesh
-from fidimag.micro import UniformExchange, UniaxialAnisotropy, Demag
+from fidimag.micro import UniformExchange, UniaxialAnisotropy
 from fidimag.common.neb_cartesian import NEB_Sundials
+from fidimag.common.neb_spherical import NEB_Sundials as NEB_Sundials_spherical
 import numpy as np
 
 # Material Parameters
@@ -25,7 +26,7 @@
 
 def two_part(pos):
 
-    x, y = pos[0], pos[1]
+    x = pos[0]
 
     if x > 6 or x < 3:
         return Ms
@@ -34,15 +35,16 @@ def two_part(pos):
 
 # Finite differences mesh
 mesh = CuboidMesh(nx=3,
-              ny=1,
-              nz=1,
-              dx=3, dy=3, dz=3,
-              unit_length=1e-9
-              )
+                  ny=1,
+                  nz=1,
+                  dx=3, dy=3, dz=3,
+                  unit_length=1e-9
+                  )
 
 
 # Simulation Function
-def relax_neb(k, maxst, simname, init_im, interp, save_every=10000):
+def relax_neb(k, maxst, simname, init_im, interp, save_every=10000,
+              coordinates='Cartesian'):
     """
     Execute a simulation with the NEB function of the FIDIMAG code, for an
     elongated particle (long cylinder)
@@ -78,28 +80,36 @@ def relax_neb(k, maxst, simname, init_im, interp, save_every=10000):
     # equal to 'the number of initial states specified', minus one.
     interpolations = interp
 
-    neb = NEB_Sundials(sim,
-                       init_images,
-                       interpolations=interpolations,
-                       spring=k,
-                       name=simname)
+    if coordinates == 'Cartesian':
+        neb = NEB_Sundials(sim,
+                           init_images,
+                           interpolations=interpolations,
+                           spring=k,
+                           name=simname)
+    elif coordinates == 'Spherical':
+        neb = NEB_Sundials_spherical(sim,
+                                     init_images,
+                                     interpolations=interpolations,
+                                     spring=k,
+                                     name=simname)
 
     neb.relax(max_steps=maxst,
               save_vtk_steps=save_every,
               save_npy_steps=save_every,
               stopping_dmdt=1e-2)
 
 
+def mid_m(pos):
+    if pos[0] > 4:
+        return (0.5, 0, 0.2)
+    else:
+        return (-0.5, 0, 0.2)
+
+
 # this test runs for over a minute
 @pytest.mark.slow
 def test_energy_barrier_2particles():
     # Initial images: we set here a rotation interpolating
-    def mid_m(pos):
-        if pos[0] > 4:
-            return (0.5, 0, 0.2)
-        else:
-            return (-0.5, 0, 0.2)
-
     init_im = [(-1, 0, 0), mid_m, (1, 0, 0)]
     interp = [6, 6]
 
@@ -114,12 +124,27 @@ def mid_m(pos):
                   interp,
                   save_every=5000)
 
+        # Relax the same system using spherical coordinates
+        relax_neb(float(k), 2000,
+                  'neb_2particles_k{}_10-10int_spherical'.format(k),
+                  init_im,
+                  interp,
+                  save_every=5000,
+                  coordinates='Spherical'
+                  )
+
     # Get the energies from the last state
     data = np.loadtxt('neb_2particles_k1e8_10-10int_energy.ndt')[-1][1:]
+    data_spherical = np.loadtxt(
+        'neb_2particles_k1e8_10-10int_spherical_energy.ndt')[-1][1:]
 
     ebarrier = np.abs(np.max(data) - np.min(data)) / (1.602e-19)
     print(ebarrier)
 
+    ebarrier_spherical = np.abs(np.max(data_spherical) -
+                                np.min(data_spherical)) / (1.602e-19)
+    print(ebarrier_spherical)
+
     # Analitically, the energy when a single particle rotates is:
     #   K V cos^2 theta
     # where theta is the angle of the direction of one particle with respect
@@ -133,6 +158,9 @@ def mid_m(pos):
     assert ebarrier < 0.017
     assert ebarrier > 0.005
 
+    assert ebarrier_spherical < 0.017
+    assert ebarrier_spherical > 0.005
+
 
-if __name__=='__main__':
-  test_energy_barrier_2particles()
+if __name__ == '__main__':
+    test_energy_barrier_2particles()