nebm fs: implemented calc of action

fidimag/common/chain_method_integrators.py

@@ -227,7 +227,7 @@ def run_for(self, n_steps):
                 # TODO: remove time from chain method rhs
                 #       make a specific function to update G??
                 self.rhs(t, self.y)
-                self.action = self.action_fun()
+                self.action = self.ChainObj.compute_action()
 
                 # self.step += 1
                 self.forces_old[:] = self.forces  # Scale field??
@@ -259,15 +259,20 @@ def run_for(self, n_steps):
 
                 # Getting averages of forces from the INNER images in the band (no extrema)
                 # (forces are given by vector G in the chain method code)
-                # TODO: we might use all band images, not only inner ones
-                G_norms = np.linalg.norm(self.forces[INNER_DOFS].reshape(-1, 3), axis=1)
+                # 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(G_norms.reshape(self.n_images - 2, -1), axis=1))
+                rms_G_norms_per_image = np.sqrt(np.mean(Gnorms2.reshape(self.n_images - 2, -1), axis=1))
                 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 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}')
+
                 # 10 seems like a magic number; we set here a minimum number of evaulations
                 if (nStart > self.nTrail * 10) and (deltaAction < self.actionTol):
                     print('Change in action is negligible')
@@ -289,7 +294,7 @@ def run_for(self, n_steps):
                     if (eta < 1e-3):
                         print('')
                         resetCount += 1
-                        bestAction = self.action_old
+                        # bestAction = self.action_old
                         self.refine_path(self.ChainObj.distances, self.band)  # Resets the path to equidistant structures (smoothing  kinks?)
                         # PathChanged[:] = True
 
@@ -329,13 +334,13 @@ def refine_path(self, distances, band):
             cs = si.CubicSpline(distances, bandrs[:, i])
             bandrs[:, i] = cs(new_dist)
 
-    def set_options(self):
-        pass
+    # def set_options(self):
+    #     pass
 
-    def _step(self, t, y, h, f):
-        """
-        """
-        pass
+    # def _step(self, t, y, h, f):
+    #     """
+    #     """
+    #     pass
 
 def normalise_spins(y):
     # Normalise an array of spins y with 3 * N elements

fidimag/common/nebm_FS.py

@@ -7,13 +7,16 @@
 from .chain_method_tools import m_to_zero_nomaterial
 from .chain_method_base import ChainMethodBase
 
+from .chain_method_integrators import FSIntegrator
+
 import scipy.integrate as spi
 
 import logging
 logging.basicConfig(level=logging.DEBUG)
 log = logging.getLogger(name="fidimag")
 
 
+# TODO: we can just inherit from the geodesic nebm class!
 class NEBM_FS(ChainMethodBase):
     """
     ARGUMENTS -----------------------------------------------------------------
@@ -153,16 +156,9 @@ def initialise_integrator(self):
         self.iterations = 0
         self.ode_count = 1
 
-        self.integrator = FSIntegrator(self.band,    # y
-                                       self.G,       # forces
-                                       self.step_RHS,
-                                       self.n_images,
-                                       self.n_dofs_image,
-                                       stepsize=1e-4)
-        self.integrator.set_options()
-        # In Verlet algorithm we only use the total force G and not YxYxG:
-        self._llg_evolve = False
-
+        # Use default integrator parameters. Pass this class object to the integrator
+        self.integrator = FSIntegrator(self)
+
     def generate_initial_band(self, method='linear'):
         """
         method      :: linear, rotation
@@ -258,7 +254,9 @@ def compute_effective_field_and_energy(self, y):
             self.energies[i] = self.sim.compute_energy()
 
         # Compute the gradient norm per every image
-        self.gradientENorm[:] = np.linalg.norm(self.gradientE, axis=1)
+        Gnorms2 = np.sum(self.gradientE.reshape(-1, 3)**2, axis=1)
+        # Compute the root mean square per image
+        self.gradientENorm[:] = np.sqrt(np.mean(Gnorms2.reshape(self.n_images, -1), axis=1))
 
         y = y.reshape(-1)
         self.gradientE = self.gradientE.reshape(-1)
@@ -314,9 +312,39 @@ def compute_spring_force(self, y):
                                        )
 
     def compute_action(self):
+        """
+        """
+        # Check that action is computed AFTER calculating and projecting the forces
+
+        # nebm_clib.image_distances_GreatCircle(self.distances,
+        #                                       self.path_distances,
+        #                                       y,
+        #                                       self.n_images,
+        #                                       self.n_dofs_image,
+        #                                       self._material_int,
+        #                                       self.n_dofs_image_material
+        #                                       )
+
         # TODO: we can use a better quadrature such as Gaussian
+        # notice that the gradient norm here is using the RMS
         action = spi.trapezoid(self.gradientENorm, self.distances)
-        # TODO: Add spring fore term to the action
+
+        # The spring term in the action is added as |F_k|^2 / (2 * self.k) = self.k * x^2 / 2
+        # (CHECK) This assumes the spring force is orthogonal to the force gradient (after projection)
+        # Knorms2 = np.sum(self.spring_force.reshape(-1, 3)**2, axis=1)
+        # # Compute the root mean square per image
+        # springF_norms = np.sqrt(np.mean(Knorms2.reshape(self.n_images, -1), axis=1))
+
+        # Norm of the spring force per image (assuming tangents as unit vectors)
+        # These are the norms of the inner images
+        dist_plus_norm = self.distances[1:]
+        dist_minus_norm = self.distances[:-1]
+        # dY_plus_norm = distances[i];
+        # dY_minus_norm = distances[i - 1];
+        springF2 = self.k * ((dist_plus_norm - dist_minus_norm)**2)
+        # CHECK: do we need to scale?
+        action += np.sum(springF2) / (self.n_images - 2)
+
         return action
 
     def compute_min_action(self):