Double checked nebm in spherical coordinates. It appears to work much slower than the Cartesian version. The test using two particles is working now for both coordinate systems giving similar results, thus Spherical coordinates are working. Other tests are too slow to implement for now, like the elongated particle, this may indicate that we need to optimise the nebm code somehow

Author: davidcorteso

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)
@@ -620,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
@@ -632,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()
 
@@ -688,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()