Merge pull request #114 from computationalmodelling/current-time-dependence

Current time dependence

Author: rpep

examples/micromagnetic/dw_stt_time_dependence/main.py

@@ -0,0 +1,199 @@
+import matplotlib as mpl
+mpl.use("Agg")
+import matplotlib.pyplot as plt
+
+import numpy as np
+from fidimag.micro import Sim
+from fidimag.common import CuboidMesh
+
+# The energies, we can use DMI in a future simulation
+from fidimag.micro import UniformExchange
+from fidimag.micro import UniaxialAnisotropy
+# from micro import DMI
+
+mu0 = 4 * np.pi * 1e-7
+
+# Initial State, a rough DW ina 1D chain
+
+
+def init_m(pos):
+
+    x = pos[0]
+
+    if x < 400:
+        return (-1, 0, 0)
+    elif 400 <= x < 500:
+        return (0, 1, 1)
+    else:
+        return (1, 0, 0)
+
+
+def relax_system(mesh):
+
+    # Only relaxation
+    sim = Sim(mesh, name='relax')
+
+    # Simulation parameters
+    sim.driver.set_tols(rtol=1e-8, atol=1e-10)
+    sim.driver.alpha = 0.5
+    sim.driver.gamma = 2.211e5
+    sim.Ms = 8.6e5
+    sim.driver.do_precession = False
+
+    # The initial state passed as a function
+    sim.set_m(init_m)
+    # sim.set_m(np.load('m0.npy'))
+    # Energies
+    A = 1.3e-11
+    exch = UniformExchange(A=A)
+    sim.add(exch)
+    anis = UniaxialAnisotropy(5e4)
+    sim.add(anis)
+    # Start relaxation and save the state in m0.npy
+    sim.relax(dt=1e-14, stopping_dmdt=0.01, max_steps=5000,
+              save_m_steps=None, save_vtk_steps=50)
+
+    np.save('m0.npy', sim.spin)
+
+
+def deal_plot(_list, output_file):
+    """
+    The list contain any number of list with two entries:
+    the _file path and the component of the magnetisation
+    in the formats:
+
+            mx, my, mz
+
+    e.g.
+            _list = [ ['m0.npy', 'mz'], ... ]
+
+    output_file     :: output file name
+
+    """
+
+    plt.figure()
+
+    comp = {'mx': 0, 'my': 1, 'mz': 2}
+
+    for element in _list:
+        # element[0]  --> _file
+        # element[1]  --> m component: mx, my or mz
+
+        data = np.load(element[0])
+
+        data.shape = (3, -1)
+
+        mx = data[comp[element[1]]]
+
+        plt.plot(mx, '--',
+                 label=element[1] + ' / %s' % element[0],
+                 dashes=(2, 2))
+
+    plt.legend()
+    # plt.xlim([0, 0.012])
+    # plt.ylim([-5, 100])
+    plt.xlabel(r'x  [ nm ]')
+    plt.savefig(output_file)
+
+
+# THIS NEEDS REVISION
+def deal_plot_dynamics(spin, m_component='mz', output_file='output.pdf'):
+    """
+    The dynamics of the m_component of the i-th spin
+    of the chain, from the magnetisation snapshot
+    m_{i} files, located at dyn_npys
+
+    m_component     :: 'mx', 'my' or 'mz'
+
+    """
+
+    comp = {'mx': 0, 'my': 1, 'mz': 2}
+
+    plt.figure()
+
+    _all = []
+    for i in range(500):
+        data = np.load('dyn2_npys/m_%d.npy' % i)
+        data.shape = (3, -1)
+        _all.append(data[comp[m_component]][spin])
+
+    plt.plot(_all, '--', label=m_component)
+
+    plt.xlabel(r't')
+    # plt.ylabel('Susceptibility')
+    plt.savefig(output_file)
+
+
+# This function excites the system with a
+# current in the x-direction using Spin Transfer Torque
+# formalism
+def excite_system(mesh):
+
+    # Specify the stt dynamics in the simulation
+    sim = Sim(mesh, name='dyn2', driver='llg_stt')
+    sim.driver.set_tols(rtol=1e-12, atol=1e-14)
+    sim.driver.alpha = 0.2
+    sim.driver.gamma = 2.211e5
+    sim.Ms = 8.6e5
+
+    # sim.set_m(init_m)
+    sim.set_m(np.load('m0.npy'))
+
+    # Energies
+    A = 1.3e-11
+    exch = UniformExchange(A=A)
+    sim.add(exch)
+
+    anis = UniaxialAnisotropy(5e4)
+    sim.add(anis)
+
+    # dmi = DMI(D=8e-4)
+    # sim.add(dmi)
+
+    # Set the current in the x direction, in A / m
+    # beta is the parameter in the STT torque
+    def jx_func(pos, t):
+        T = 1e-9
+        return (-1e12 +  -0.1e12 * np.sin(t/T))
+
+    sim.driver.jx_function = jx_func
+    sim.driver.beta = 0.01
+
+    # The simulation will run for 5 ns and save
+    # 500 snapshots of the system in the process
+    ts = np.linspace(0, 10e-9, 1001)
+
+    for t in ts:
+        print('time', t)
+        sim.driver.run_until(t)
+        print('j = {}'.format(sim.driver._jx[0]))
+        sim.save_vtk()
+        sim.save_m()
+
+if __name__ == '__main__':
+    # We will crate a mesh with 1000 elements of elements
+    # in the x direction, and 1 along y and z
+    # (so we have a 1D system)
+    mesh = CuboidMesh(nx=1000, ny=1, nz=1,
+                  dx=2, dy=2, dz=2.0,
+                  unit_length=1e-9)
+
+    # Relax the initial state
+    relax_system(mesh)
+    deal_plot([['m0.npy', 'mx']], 'initial_state.pdf')
+
+    # Now we excite the system with the current
+    excite_system(mesh)
+
+    # Plot the dynamics of a single spin
+    deal_plot_dynamics(0)
+
+    # We can plot the m_z component for a number snapshots
+    # to observe the DW motion
+    # We will plot the 200th and 499th file (the last one
+    # is the system at ~5 ns)
+    deal_plot([['m0.npy', 'mx'],
+               ['dyn2_npys/m_200.npy', 'mx'],
+               ['dyn2_npys/m_499.npy', 'mx']
+               ],
+              'dw_motion_stt.pdf')

fidimag/atomistic/llg_stt.py

@@ -23,12 +23,12 @@ class LLG_STT(AtomisticDriver):
           dt             2
                   ( 1 + a  )
 
-                        u  
+                        u
                   +  --------  s X [ (1 + a * b) (j . grad) s  - (b - a) s X (j . grad) s ]
                            2
                     ( 1 + a  )
 
-        
+
         with         g mu_B P
                 u =  --------
                      2 mu_s e
@@ -67,10 +67,13 @@ def __init__(self, mesh, spin, mu_s, mu_s_inv, field, pins,
 
         self.p = 0.5
         self.beta = 0
-        self.update_j_fun = None
+        self.jx_function = None
+        self.jy_function = None
+        self.jz_function = None
 
         # FIXME: change the u0 to spatial
-        v = self.mesh.dx * self.mesh.dy * self.mesh.dz * (self.mesh.unit_length ** 3)
+        v = self.mesh.dx * self.mesh.dy * \
+            self.mesh.dz * (self.mesh.unit_length ** 3)
         self.u0 = const.g_e * const.mu_B / (2 * const.c_e) * v
 
     def get_jx(self):
@@ -98,38 +101,28 @@ def set_jz(self, value):
     jz = property(get_jz, set_jz)
 
     def sundials_rhs(self, t, y, ydot):
-
         self.t = t
-
         # already synchronized when call this funciton
         # self.spin[:]=y[:]
-
         self.compute_effective_field(t)
-
-        if self.update_j_fun is not None:
-            clib.compute_stt_field(self.spin,
-                                   self.field_stt,
-                                   self._jx * self.update_j_fun(t),
-                                   self._jy * self.update_j_fun(t),
-                                   self._jz * self.update_j_fun(t),
-                                   self.mesh.dx * self.mesh.unit_length,
-                                   self.mesh.dy * self.mesh.unit_length,
-                                   self.mesh.dz * self.mesh.unit_length,
-                                   self.mesh.neighbours,
-                                   self.n
-                                   )
-        else:
-            clib.compute_stt_field(self.spin,
-                                   self.field_stt,
-                                   self._jx,
-                                   self._jy,
-                                   self._jz,
-                                   self.mesh.dx * self.mesh.unit_length,
-                                   self.mesh.dy * self.mesh.unit_length,
-                                   self.mesh.dz * self.mesh.unit_length,
-                                   self.mesh.neighbours,
-                                   self.n
-                                   )
+        if self.jx_function:
+            self.set_jx(self.jx_function, t)
+        if self.jy_function:
+            self.set_jy(self.jy_function, t)
+        if self.jz_function:
+            self.set_jz(self.jz_function, t)
+
+        clib.compute_stt_field(self.spin,
+                               self.field_stt,
+                               self._jx,
+                               self._jy,
+                               self._jz,
+                               self.mesh.dx * self.mesh.unit_length,
+                               self.mesh.dy * self.mesh.unit_length,
+                               self.mesh.dz * self.mesh.unit_length,
+                               self.mesh.neighbours,
+                               self.n
+                               )
 
         clib.compute_llg_stt_rhs(ydot,
                                  self.spin,

fidimag/atomistic/llg_stt_cpp.py

@@ -47,13 +47,14 @@ def __init__(self, mesh, spin, mu_s, mu_s_inv, field, pins,
 
         self.a_J = 1
         self.beta = 0
-        self.J_time_fun = None
+        self.j_function = None
 
     def get_p(self):
         return self._p
 
     def set_p(self, value):
         self._p[:] = helper.init_vector(value, self.mesh)
+
     p = property(get_p, set_p)
 
     def get_a_J(self):
@@ -73,25 +74,16 @@ def sundials_rhs(self, t, y, ydot):
 
         self.compute_effective_field(t)
 
-        if self.J_time_fun is not None:
-            clib.compute_llg_stt_cpp(ydot,
-                                     self.spin,
-                                     self.field,
-                                     self._p,
-                                     self._alpha,
-                                     self._pins,
-                                     self._a_J*self.J_time_fun(t),
-                                     self.beta,
-                                     self.gamma,
-                                     self.n)
-        else:
-            clib.compute_llg_stt_cpp(ydot,
-                                     self.spin,
-                                     self.field,
-                                     self._p,
-                                     self._alpha,
-                                     self._pins,
-                                     self._a_J,
-                                     self.beta,
-                                     self.gamma,
-                                     self.n)
+        if self.j_function is not None:
+            self.set_a_J(self.j_function)
+
+        clib.compute_llg_stt_cpp(ydot,
+                                 self.spin,
+                                 self.field,
+                                 self._p,
+                                 self._alpha,
+                                 self._pins,
+                                 self._a_J,
+                                 self.beta,
+                                 self.gamma,
+                                 self.n)

fidimag/common/__init__.py

@@ -4,3 +4,4 @@
 from .citation import citation
 from .cuboid_mesh import CuboidMesh
 #from .neb_cartesian import NEB_Sundials
+from .helper import init_scalar, init_vector

fidimag/common/helper.py

@@ -24,7 +24,7 @@ def normalise(a):
     a.shape = (-1,)
 
 
-def init_vector(m0, mesh, norm=False):
+def init_vector(m0, mesh, norm=False, *args):
 
     n = mesh.n
 
@@ -36,7 +36,7 @@ def init_vector(m0, mesh, norm=False):
 
     elif hasattr(m0, '__call__'):
         for i in range(n):
-            v = m0(mesh.coordinates[i])
+            v = m0(mesh.coordinates[i], *args)
             if len(v) != 3:
                 raise Exception(
                     'The length of the value in init_vector method must be 3.')
@@ -58,7 +58,7 @@ def init_vector(m0, mesh, norm=False):
     return spin
 
 
-def init_scalar(value, mesh):
+def init_scalar(value, mesh, *args):
 
     n = mesh.n
 
@@ -68,7 +68,7 @@ def init_scalar(value, mesh):
         mesh_v[:] = value
     elif hasattr(value, '__call__'):
         for i in range(n):
-            mesh_v[i] = value(mesh.coordinates[i])
+            mesh_v[i] = value(mesh.coordinates[i], *args)
 
     elif isinstance(value, np.ndarray):
         if value.shape == mesh_v.shape: