From 46fb348476d07cd4d6d3d9950465e21b2ef206b9 Mon Sep 17 00:00:00 2001 From: davidcortesortuno <david.cortes.o@gmail.com> Date: Sat, 5 Oct 2024 13:35:09 -0300 Subject: [PATCH] Debugging action calculation --- fidimag/common/chain_method_integrators.py | 12 ++++++++---- fidimag/common/nebm_FS.py | 21 +++++++++++++-------- tests/test_two_particles_nebm-fs.py | 9 +++++---- 3 files changed, 26 insertions(+), 16 deletions(-) diff --git a/fidimag/common/chain_method_integrators.py b/fidimag/common/chain_method_integrators.py index 823b8400..c86c087b 100644 --- a/fidimag/common/chain_method_integrators.py +++ b/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,7 +271,8 @@ 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 @@ -279,11 +280,14 @@ def run_for(self, n_steps): # 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) diff --git a/fidimag/common/nebm_FS.py b/fidimag/common/nebm_FS.py index 958c4c38..fad86ae1 100644 --- a/fidimag/common/nebm_FS.py +++ b/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, diff --git a/tests/test_two_particles_nebm-fs.py b/tests/test_two_particles_nebm-fs.py index 9c46382d..fade321a 100644 --- a/tests/test_two_particles_nebm-fs.py +++ b/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 = []