Debugging action calculation

fidimag/common/chain_method_integrators.py

@@ -160,7 +160,7 @@ def __init__(self, ChainObj,
                  # band, forces, distances, rhs_fun, action_fun,
                  # n_images, n_dofs_image,
                  maxSteps=1000,
-                 maxCreep=6, actionTol=1e-10, forcesTol=1e-6,
+                 maxCreep=6, actionTol=1e-2, forcesTol=1e-6,
                  etaScale=1e-6, dEta=4., minEta=1e-6,
                  # perturbSeed=42, perturbFactor=0.1,
                  nTrail=13, resetMax=20
@@ -249,7 +249,7 @@ def run_for(self, n_steps):
             # Creep stage: minimise with a fixed eta
             while creepCount < self.maxCreep:
                 # Update spin. Avoid pinned or zero-Ms sites
-                self.band[:] = self.band_old + eta * self.etaScale * self.forces_old
+                self.band[:] = self.band_old - eta * self.etaScale * self.forces_old
                 normalise_spins(self.band)
 
                 self.trailAction
@@ -271,19 +271,23 @@ def run_for(self, n_steps):
                 # TODO: we might use all band images, not only inner ones, although G is zero at the extrema
                 Gnorms2 = np.sum(self.forces[INNER_DOFS].reshape(-1, 3)**2, axis=1)
                 # Compute the root mean square per image
-                rms_G_norms_per_image = np.sqrt(np.mean(Gnorms2.reshape(self.n_images - 2, -1), axis=1))
+                rms_G_norms_per_image = np.sum(Gnorms2.reshape(self.n_images - 2, -1), axis=1) / self.ChainObj.n_spins
+                rms_G_norms_per_image = np.sqrt(rms_G_norms_per_image)
                 mean_rms_G_norms_per_image = np.mean(rms_G_norms_per_image)
 
                 # Average step difference between trailing action and new action
                 deltaAction = (np.abs(self.trailAction[nStart] - self.action)) / self.nTrail
 
                 # print('trail Actions', self.trailAction)
 
+                ma = self.ChainObj.compute_min_action()
+
                 # Print log
                 print(f'Step {self.i_step}  ⟨RMS(G)〉= {mean_rms_G_norms_per_image:.5e}  ',
                       f'deltaAction = {deltaAction:.5e}  Creep n = {creepCount:>3}  resetC = {resetCount:>3}  ',
                       f'eta = {eta:>5.4e}  '
-                      f'action = {self.action:>5.4e}  action_old = {self.action_old:>5.4e}'
+                      f'action = {self.action:>5.4e}  action_old = {self.action_old:>5.4e} '
+                      f'MIN action = {ma:>5.4e}'
                       )
                 # print(self.forces)
 

fidimag/common/nebm_FS.py

@@ -263,9 +263,9 @@ def compute_tangents(self, y):
         nebm_clib.project_images(self.tangents, y,
                                  self.n_images, self.n_dofs_image
                                  )
-        nebm_clib.normalise_images(self.tangents,
-                                   self.n_images, self.n_dofs_image
-                                   )
+        # nebm_clib.normalise_images(self.tangents,
+        #                            self.n_images, self.n_dofs_image
+        #                            )
 
     def compute_spring_force(self, y):
         """
@@ -322,7 +322,9 @@ def compute_action(self):
 
         # NOTE: Gradient here is projected in the S2^N tangent space
         self.gradientE.shape = (self.n_images, -1)
-        Gnorms2 = np.sum(self.gradientE**2, axis=1) / self.n_images
+        # NOTE: HEre we have to divide by the number of spins per image,not n_images:
+        Gnorms2 = np.sum(self.gradientE**2, axis=1) / self.n_spins
+
         # Compute the root mean square per image
         self.gradientENorm[:] = np.sqrt(Gnorms2)
         self.gradientE.shape = (-1)
@@ -333,6 +335,10 @@ def compute_action(self):
         # TODO: we can use a better quadrature such as Gaussian
         # notice that the gradient norm here is using the RMS
         action = spi.simpson(self.gradientENorm, x=self.path_distances)
+        # print('E', self.energies / (self.mesh.dx * self.mesh.dy * self.mesh.dz * self.mesh.unit_length**3))
+        # print('gradE norm', self.gradientENorm)
+        # print('Path distance', self.path_distances)
+        print('Images', self.band.reshape(-1, 3).reshape(self.n_images, -1))
 
         # DEBUG:
         # print('action from gradE', action)
@@ -361,15 +367,14 @@ def compute_action(self):
     def compute_min_action(self):
         dE = self.energies[-1] - self.energies[0]
         minAction = np.sum(np.abs(self.energies[1:] - self.energies[:-1]))
-        return 2 * (dE + minAction)
+        return 2 * (dE + minAction) / (self.mesh.dx * self.mesh.dy * self.mesh.dz * self.mesh.unit_length**3) / self.path_distances[-1]
 
     def nebm_step(self, y, ensure_zero_extrema=False):
 
         self.compute_effective_field_and_energy(y)
-        nebm_clib.project_images(self.gradientE, y,
-                                 self.n_images, self.n_dofs_image
-                                 )
+        nebm_clib.project_images(self.gradientE, y, self.n_images, self.n_dofs_image)
         self.compute_tangents(y)
+        nebm_clib.normalise_spins(self.tangents, self.n_images, self.n_dofs_image)
         self.compute_spring_force(y)
 
         nebm_clib.compute_effective_force(self.G,

tests/test_two_particles_nebm-fs.py

@@ -3,9 +3,10 @@
 # FIDIMAG:
 from fidimag.micro import Sim
 from fidimag.common import CuboidMesh
-from fidimag.micro import UniformExchange, UniaxialAnisotropy
+from fidimag.micro import UniformExchange, UniaxialAnisotropy, Demag
 from fidimag.common.nebm_FS import NEBM_FS
 import numpy as np
+import matplotlib.pyplot as plt
 
 # Material Parameters
 # Parameters
@@ -61,12 +62,12 @@ def relax_string(maxst, simname, init_im, interp, save_every=10000):
     interpolations = interp
 
     nebm = NEBM_FS(sim, init_im, interpolations=interpolations, name=simname,
-                   interpolation_method='linear')
+                   interpolation_method='rotation', spring_constant=1e5)
 
     # dt = integrator.stepsize means after every integrator step, the images
     # are rescaled. We can run more integrator steps if we decrease the
     # stepsize, e.g. dt=1e-3 and integrator.stepsize=1e-4
-    nebm.integrator.maxSteps = 33
+    nebm.integrator.maxSteps = 400
     nebm.integrator.run_for(maxst,
                             # save_vtks_every=save_every,
                             # save_npys_every=save_every,
@@ -84,7 +85,7 @@ def mid_m(pos):
 
 def test_energy_barrier_2particles_string():
     # Initial images: we set here a rotation interpolating
-    init_im = [(-1, 0, 0), (0.0, 0.0, 1.0), (1, 0, 0)]
+    init_im = [(-1, 0, 0), (0.0, 0.2, 1.0), (1, 0, 0)]
     interp = [6, 6]
 
     barriers = []