Merge pull request #59 from ubermag/time-dependence

Time-dependent fields and currents.

oommfc/drivers/driver.py

@@ -1,4 +1,5 @@
 import os
+import numpy as np
 import sys
 import abc
 import glob
@@ -170,6 +171,12 @@ def drive(self, system, /, dirname='.', append=True, fixed_subregions=None,
         if not os.path.exists(workingdir):
             os.makedirs(workingdir)
 
+        # compute tlist for time-dependent field (current)
+        for term in system.energy:
+            if (hasattr(term, 'time_dependence')
+                    and callable(term.time_dependence)):
+                self._time_dependence(term=term, **kwargs)
+
         # Change directory to workingdir
         with _changedir(workingdir):
             # Generate the necessary filenames.
@@ -185,6 +192,14 @@ def drive(self, system, /, dirname='.', append=True, fixed_subregions=None,
             with open(miffilename, 'w') as miffile:
                 miffile.write(mif)
 
+            # free memory
+            for term in system.energy:
+                if (hasattr(term, 'time_dependence')
+                        and callable(term.time_dependence)):
+                    del term.tlist
+                    del term.dtlist
+                    del mif
+
             # Generate and save json info file for a drive (not compute).
             if compute is None:
                 info = {}
@@ -238,3 +253,20 @@ def drive(self, system, /, dirname='.', append=True, fixed_subregions=None,
         # remove information about fixed cells for subsequent runs
         if hasattr(self.evolver, 'fixed_spins'):
             del self.evolver.fixed_spins
+
+    def _time_dependence(self, term, **kwargs):
+        try:
+            tmax = kwargs['t']
+        except KeyError:
+            msg = (f'Time-dependent term {term.__class__.__name__=} must be '
+                   'used with time driver.')
+            raise RuntimeError(msg)
+        ts = np.arange(0, tmax + term.tstep, term.tstep)
+        try:  # vector output from term.time_dependence
+            tlist = [list(term.time_dependence(t)) for t in ts]
+            dtlist = (np.gradient(tlist)[0] / term.tstep).tolist()
+        except TypeError:
+            tlist = [term.time_dependence(t) for t in ts]
+            dtlist = list(np.gradient(tlist) / term.tstep)
+        term.tlist = tlist
+        term.dtlist = dtlist

oommfc/evolvers/spintevolver.py

@@ -40,5 +40,7 @@ class SpinTEvolver(mm.Evolver):
                            'gamma_G',
                            'do_precess',
                            'u',
+                           'u_profile',
+                           'u_profile_args',
                            'beta',
                            'method']

oommfc/scripts/driver.py

@@ -143,6 +143,10 @@ def driver_script(driver, system, fixed_subregions=None, compute=None,
         if mm.ZhangLi() in system.dynamics:
             driver.evolver.u = system.dynamics.zhangli.u
             driver.evolver.beta = system.dynamics.zhangli.beta
+            for arg in ['time_dependence', 'tstep', 'tcl_strings']:
+                if hasattr(system.dynamics.zhangli, arg):
+                    setattr(driver.evolver, arg,
+                            getattr(system.dynamics.zhangli, arg))
         if mm.Slonczewski() in system.dynamics:
             driver.evolver.J = system.dynamics.slonczewski.J
             driver.evolver.mp = system.dynamics.slonczewski.mp
@@ -154,6 +158,13 @@ def driver_script(driver, system, fixed_subregions=None, compute=None,
             else:
                 driver.evolver.eps_prime = \
                     system.dynamics.slonczewski.eps_prime
+            for arg in ['time_dependence', 'tstep', 'tcl_strings']:
+                if (hasattr(system.dynamics.slonczewski, arg)
+                        and not isinstance(
+                            getattr(system.dynamics.slonczewski, arg),
+                            ts.Descriptor)):
+                    setattr(driver.evolver, arg,
+                            getattr(system.dynamics.slonczewski, arg))
         if isinstance(driver.evolver, (oc.UHH_ThetaEvolver,
                                        oc.Xf_ThermHeunEvolver,
                                        oc.Xf_ThermSpinXferEvolver)):
@@ -169,7 +180,7 @@ def driver_script(driver, system, fixed_subregions=None, compute=None,
             resstr = f'{{main_atlas {" ".join(fixed_subregions)}}}'
             driver.evolver.fixed_spins = resstr
 
-        mif += oc.scripts.evolver_script(driver.evolver)
+        mif += oc.scripts.evolver_script(driver.evolver, **kwargs)
 
         # Extract time and number of steps.
         t, n = kwargs['t'], kwargs['n']

oommfc/scripts/energy.py

@@ -144,6 +144,74 @@ def zeeman_script(term, system):
             mif += f'  script_args total_time\n'
             mif += f'  script TimeFunction\n'
             mif += '}\n\n'
+    elif hasattr(term, 'tlist'):
+        if isinstance(term.H, (df.Field, dict)):
+            mif += 'proc TimeFunction { total_time } {\n'
+            mif += f'  set tstep {term.tstep}\n'
+            mif += '  set index [expr round($total_time/$tstep)]\n'
+            tstr = ' '.join(map(str, term.tlist))
+            for char, replacement in [[',', ''], ['[', '{ '], [']', ' }']]:
+                tstr = tstr.replace(char, replacement)
+            mif += f'  set H_t_fac {{ {tstr} }}\n'
+            dtstr = ' '.join(map(str, term.dtlist))
+            for char, replacement in [[',', ''], ['[', '{ '], [']', ' }']]:
+                dtstr = dtstr.replace(char, replacement)
+            mif += f'  set dH_t_fac {{ {dtstr} }}\n'
+            mif += '  set H [lindex $H_t_fac $index]\n'
+            mif += '  set dH [lindex $dH_t_fac $index]\n'
+            if isinstance(term.tlist[0], list):
+                mif += (f'  return [list'
+                        f' {" ".join([f"[lindex $H {i}]" for i in range(9)])}'
+                        f' {" ".join([f"[lindex $dH {i}]" for i in range(9)])}'
+                        ']\n')
+            else:
+                mif += '  return [list $H $H $H $dH $dH $dH]\n'
+            mif += '}\n\n'
+
+            mif += '# TransformZeeman\n'
+            mif += f'Specify Oxs_TransformZeeman:{term.name} {{\n'
+            if isinstance(term.tlist[0], list):
+                mif += '  type general\n'
+            else:
+                mif += '  type diagonal\n'
+            mif += '  script_args total_time\n'
+            mif += '  script TimeFunction\n'
+            mif += f'  field {Hname}\n'
+            mif += '}\n\n'
+        else:
+            mif += 'proc TimeFunction { total_time } {\n'
+            mif += f'  set tstep {term.tstep}\n'
+            mif += '  set index [expr round($total_time/$tstep)]\n'
+            mif += f'  set H_t_fac {{ {" ".join(map(str, term.tlist))} }}\n'
+            mif += f'  set dH_t_fac {{ {" ".join(map(str, term.dtlist))} }}\n'
+            mif += '  set H_fac [lindex $H_t_fac $index]\n'
+            mif += '  set dH_fac [lindex $dH_t_fac $index]\n'
+            mif += f'  set Hx [expr {{ {term.H[0]}*$H_fac }}]\n'
+            mif += f'  set Hy [expr {{ {term.H[1]}*$H_fac }}]\n'
+            mif += f'  set Hz [expr {{ {term.H[2]}*$H_fac }}]\n'
+            mif += f'  set dHx [expr {{ {term.H[0]}*$dH_fac }}]\n'
+            mif += f'  set dHy [expr {{ {term.H[1]}*$dH_fac }}]\n'
+            mif += f'  set dHz [expr {{ {term.H[2]}*$dH_fac }}]\n'
+            mif += '  return [list $Hx $Hy $Hz $dHx $dHy $dHz]\n'
+            mif += '}\n\n'
+
+            mif += '# ScriptUZeeman\n'
+            mif += f'Specify Oxs_ScriptUZeeman:{term.name} {{\n'
+            mif += f'  script_args total_time\n'
+            mif += f'  script TimeFunction\n'
+            mif += '}\n\n'
+    elif isinstance(term.tcl_strings, dict):
+        mif += term.tcl_strings['proc']
+        mif += f'\n# {term.tcl_strings["energy"][4:]}\n'  # 3.9 removeprefix
+        mif += f'Specify {term.tcl_strings["energy"]}:{term.name} {{\n'
+        for key in ['type', 'script_args', 'script']:
+            try:
+                mif += f'  {key} {term.tcl_strings[key]}\n'
+            except KeyError:
+                pass
+        if term.tcl_strings['energy'] == 'Oxs_TransformZeeman':
+            mif += f'  field {Hname}\n'
+        mif += '}\n\n'
     else:
         mif += '# FixedZeeman\n'
         mif += f'Specify Oxs_FixedZeeman:{term.name} {{\n'

oommfc/scripts/evolver.py

@@ -1,7 +1,9 @@
 import oommfc as oc
+import numbers
+import numpy as np
 
 
-def evolver_script(evolver):
+def evolver_script(evolver, **kwargs):
     mif = ''
 
     # Prepare parameters depending on what attributes are defined in evolver.
@@ -35,8 +37,14 @@ def evolver_script(evolver):
         evolver.P = Pname
         mif += Pmif
     if hasattr(evolver, 'Lambda'):
+        Lambda = evolver.Lambda
+        if isinstance(evolver.Lambda, dict):
+            # the automatically added default value 0 is not allowed
+            Lambda = evolver.Lambda.copy()
+            if 'default' not in Lambda:
+                Lambda['default'] = 1
         Lambdamif, Lambdaname = oc.scripts.setup_scalar_parameter(
-            evolver.Lambda, 'sl_Lambda')
+            Lambda, 'sl_Lambda')
         evolver.Lambda = Lambdaname
         mif += Lambdamif
     if hasattr(evolver, 'eps_prime'):
@@ -45,6 +53,39 @@ def evolver_script(evolver):
         evolver.eps_prime = eps_primename
         mif += eps_primemif
 
+    if hasattr(evolver, 'tstep') and isinstance(evolver.tstep, numbers.Real):
+        ts = np.arange(0, kwargs['t'] + evolver.tstep, evolver.tstep)
+        tlist = [evolver.time_dependence(t) for t in ts]
+
+        mif += 'proc TimeFunction { total_time } {\n'
+        mif += f'  set tstep {evolver.tstep}\n'
+        mif += '  set index [expr round($total_time/$tstep)]\n'
+        mif += f'  set profile {{ {" ".join(map(str, tlist))} }}\n'
+        mif += '  set factor [lindex $profile $index]\n'
+        mif += '  return $factor\n'
+        mif += '}\n\n'
+
+        if isinstance(evolver, (oc.SpinXferEvolver)):
+            # oc.Xf_TherrmSpinXferEvolver)):
+            setattr(evolver, 'J_profile', 'TimeFunction')
+            setattr(evolver, 'J_profile_args', 'total_time')
+        elif isinstance(evolver, oc.SpinTEvolver):
+            setattr(evolver, 'u_profile', 'TimeFunction')
+            setattr(evolver, 'u_profile_args', 'total_time')
+    if hasattr(evolver, 'tcl_strings') and isinstance(evolver.tcl_strings,
+                                                      dict):
+        print(evolver.tcl_strings)
+        mif += evolver.tcl_strings['proc']
+        if isinstance(evolver, (oc.SpinXferEvolver)):
+            # oc.Xf_ThermSpinXferEvolver)):
+            setattr(evolver, 'J_profile', evolver.tcl_strings['proc_name'])
+            setattr(evolver, 'J_profile_args',
+                    evolver.tcl_strings['proc_args'])
+        elif isinstance(evolver, oc.SpinTEvolver):
+            setattr(evolver, 'u_profile', evolver.tcl_strings['proc_name'])
+            setattr(evolver, 'u_profile_args',
+                    evolver.tcl_strings['proc_args'])
+
     # temperature cannot spacially vary
 
     # Scripts for a specific evolver.
@@ -55,9 +96,11 @@ def evolver_script(evolver):
         mif += '# RungeKuttaEvolver\n'
         mif += 'Specify Oxs_RungeKuttaEvolve:evolver {\n'
     elif isinstance(evolver, oc.SpinTEvolver):
+        # time dependence
         mif += '# Zhang-Li evolver\n'
         mif += 'Specify Anv_SpinTEvolve:evolver {\n'
     elif isinstance(evolver, oc.SpinXferEvolver):
+        # time dependence
         mif += '# Slonczewski evolver\n'
         mif += 'Specify Oxs_SpinXferEvolve:evolver {\n'
     elif isinstance(evolver, oc.CGEvolver):
@@ -75,7 +118,8 @@ def evolver_script(evolver):
 
     # Define all other parameters.
     for attr, value in evolver:
-        mif += f'  {attr} {value}\n'
+        if attr not in ['time_dependence', 'tstep', 'tcl_strings']:
+            mif += f'  {attr} {value}\n'
     mif += '}\n\n'
 
     return mif

oommfc/scripts/util.py

@@ -82,9 +82,12 @@ def setup_scalar_parameter(parameter, name):
         return mif, f'{name}_norm'
 
     elif isinstance(parameter, dict):
-        if 'default' not in parameter.keys():
-            parameter['default'] = 0
-        mif = atlas_scalar_field(parameter, f'{name}')
+        # to avoid changing the dictionary used in the respective term
+        # in the system (from micromagneticmodel)
+        param = parameter.copy()
+        if 'default' not in param.keys():
+            param['default'] = 0
+        mif = atlas_scalar_field(param, f'{name}')
         return mif, f'{name}'
 
     elif isinstance(parameter, numbers.Real):