Fix bugs and improve modularity

README.md

@@ -65,7 +65,7 @@ Please, read the full description at [MRAM Framework Description](./doc/README.m
 	* `fokker_plank`
 		* [README.md](./fokker_plank/README.md) for the MRAM Fokker-Plank description
 		* `fvm`
-			* `fvm_classes.py` Finite Volume Method classes, see [FVM](https://github.com/danieljfarrell/FVM)
+			* `fvm_classes.py` Finite Volume Method classes, see [FVM](https://github.com/danieljfarrell/FVM and [FVM2](https://github.com/Fratorhe/FVM/))
 			* `mtj_fp_fvm.py` MTJ Fokker-Plank FVM solver
 		* `analytical`
 			* `analytical.py` MTJ Fokker-Plank Analytical solver

src/README.md

@@ -6,8 +6,8 @@ Summary:
 * `test_sllgs_solver.py` shows you the basic s-LLGS solver config, calling (`sllgs_solver.py`)
 * `stochastic_multithread_simulation.py` (calling `sllgs_solver.py`) is the script
 that helps you launching parallel python s-LLGS simulations
-* These results can be compared against Fooker-Plank simulations (see `plot_sllgs_fpe_comparison.py` script)
-* `analytical.py` and `mtj_fp_fvm.py` contain the Fooker-Plank solvers. Analytical contains the WER/RER fitter for the problem optimization
+* These results can be compared against Fokker-Plank simulations (see `plot_sllgs_fpe_comparison.py` script)
+* `analytical.py` and `mtj_fp_fvm.py` contain the Fokker-Plank solvers. Analytical contains the WER/RER fitter for the problem optimization
 * Verilog-a compact models: run the testbenches `tb.scs` and `tb_subckt.scs`
 
 Please, read the full description at [MRAM Framework Description](./doc/README.md).

src/fokker_plank/README.md

@@ -1,6 +1,6 @@
 ## Quick Start & More info
 Summary:
-* Fooker-Plank simulations (see `plot_sllgs_fpe_comparison.py` script) `analytical.py` and `mtj_fp_fvm.py`
+* Fokker-Plank simulations (see `plot_sllgs_fpe_comparison.py` script) `analytical.py` and `mtj_fp_fvm.py`
 
 Please, read the full description at [MRAM Framework Description](./doc/README.md).
 

src/fokker_plank/fvm/fvm_lib/fvm_classes.py

@@ -331,7 +331,7 @@ def alpha_matrix(self):
         system it is equal to the identity matrix.
         """
         a1 = 0 if self.left_flux is None else 1
-        aJ = 0 if self.left_flux is None else 1
+        aJ = 0 if self.right_flux is None else 1
         diagonals = np.ones(self.mesh.J)
         diagonals[0] = a1
         diagonals[-1] = aJ

src/fokker_plank/fvm/mtj_fp_fvm.py

@@ -9,40 +9,38 @@
 MTJ Fokker-Plank Finite Volume Method Solver.
 
 Fokker-Plank or advection-diffusion for
-MTJ magnetization probability evolution. 
+MTJ magnetization probability evolution.
 """
 
+import matplotlib.cm as cm
+import matplotlib.pyplot as plt
 import numpy as np
 from scipy import sparse
 from scipy.sparse import linalg
-import matplotlib.pyplot as plt
-import matplotlib.cm as cm
 
 import fvm_lib.fvm_classes as fvm
 
 np.random.seed(seed=1)
 
 
-def solve_mtj_fp(dim_points=1000,
-                 rho_init=None,
-                 delta=60,
-                 i0=1.5,
-                 h=0.,
-                 t_step=0.001,
-                 T=10,
-                 lin_space_z=False,
-                 do_3d=True):
+def solve_mtj_fp(
+    dim_points=1000,
+    rho_init=None,
+    delta=60,
+    i0=1.5,
+    h=0.0,
+    t_step=0.001,
+    T=10,
+    lin_space_z=False,
+    discretization="exponential",  # exponential, upwind, central
+    do_3d=True,
+):
     """Solve FVM for the FPE on a perpendicular symetric MTJ."""
     # cranck-nickolson
     theta = 0.5
     # fully implicit
     # theta = 1.
 
-    # discretization
-    discretization = 'exponential'
-    # discretization = 'upwind'
-    # discretization = 'central'
-
     # Dirichlet boundary conditions
     # right_value = 1.0
     # left_value = 1.0
@@ -59,9 +57,9 @@ def solve_mtj_fp(dim_points=1000,
     mesh = fvm.Mesh(faces)
 
     # drift/advection/convection
-    U = (i0-h-mesh.cells)*(1-mesh.cells*mesh.cells)
+    U = (i0 - h - mesh.cells) * (1 - mesh.cells * mesh.cells)
     # diffusion
-    D = (1-mesh.cells*mesh.cells)/(2*delta)
+    D = (1 - mesh.cells * mesh.cells) / (2 * delta)
 
     # drift/advection/convection
     a = fvm.CellVariable(-U, mesh=mesh)
@@ -75,106 +73,95 @@ def solve_mtj_fp(dim_points=1000,
         # w_init = np.exp(
         #     -delta*(1-mesh.cells*mesh.cells))*np.heaviside(mesh.cells, 0.5)
         _theta_x = np.arccos(mesh.cells)
-        rho_init = np.exp(
-            -delta*np.sin(_theta_x)*np.sin(_theta_x))*np.heaviside(mesh.cells,
-                                                                   0.5)
+        rho_init = np.exp(-delta * np.sin(_theta_x) * np.sin(_theta_x)) * np.heaviside(
+            mesh.cells, 0.5
+        )
         rho_init /= np.trapz(rho_init, x=mesh.cells)
     # w_init = w_init[::-1]
-    print(f'\trho_init area: {np.trapz(rho_init, x=mesh.cells)}')
+    print(f"\trho_init area: {np.trapz(rho_init, x=mesh.cells)}")
 
     model = fvm.AdvectionDiffusionModel(
-        faces, a, d, t_step, discretization=discretization)
+        faces, a, d, t_step, discretization=discretization
+    )
     # model.set_boundary_conditions(left_value=1., right_value=0.)
     model.set_boundary_conditions(left_flux=left_flux, right_flux=right_flux)
     M = model.coefficient_matrix()
     alpha = model.alpha_matrix()
     beta = model.beta_vector()
     identity = sparse.identity(model.mesh.J)
-    print(f'\talpha: [{np.min(alpha.todense())}, {np.max(alpha.todense())}]')
-    print(f'\tbeta: [{np.min(beta)}, {np.max(beta)}]')
-    print(f'\tidentity: {identity.shape}')
+    print(f"\talpha: [{np.min(alpha.todense())}, {np.max(alpha.todense())}]")
+    print(f"\tbeta: [{np.min(beta)}, {np.max(beta)}]")
+    print(f"\tidentity: {identity.shape}")
 
     # Construct linear system from discretised matrices, A.x = d
-    A = identity - t_step*theta*alpha*M
-    d = (identity + t_step*(1-theta)*alpha*M)*rho_init + beta
+    A = identity - t_step * theta * alpha * M
+    d = (identity + t_step * (1 - theta) * alpha * M) * rho_init + beta
 
-    print("\tPeclet number", np.min(model.peclet_number()),
-          np.max(model.peclet_number()))
-    print("\tCFL condition", np.min(model.CFL_condition()),
-          np.max(model.CFL_condition()))
+    print(
+        "\tPeclet number", np.min(model.peclet_number()), np.max(model.peclet_number())
+    )
+    print(
+        "\tCFL condition", np.min(model.CFL_condition()), np.max(model.CFL_condition())
+    )
 
     rho = rho_init
 
-    t0 = np.linspace(0, T, int(T/t_step)+1)
+    t0 = np.linspace(0, T, int(T / t_step) + 1)
     if do_3d:
         rho = np.zeros((t0.shape[0], mesh.cells.shape[0]))
         area = np.zeros(t0.shape[0])
         rho[0] = rho_init
         area[0] = np.trapz(rho[0], x=mesh.cells)
         for i in range(1, t0.shape[0]):
-            d = (identity + t_step*(1-theta)*alpha*M)*rho[i-1] + beta
+            d = (identity + t_step * (1 - theta) * alpha * M) * rho[i - 1] + beta
             rho[i] = linalg.spsolve(A, d)
             # normalization not needed, flux is kept
         # PS/PNS
-        ps = np.trapz(rho.T[mesh.cells < 0],
-                      x=mesh.cells[mesh.cells < 0], axis=0)
+        ps = np.trapz(rho.T[mesh.cells < 0], x=mesh.cells[mesh.cells < 0], axis=0)
         t_sw = t0[np.argmax(ps > 0.5)]
     else:
         rho = np.array(rho_init)
         rho_next = np.array(rho_init)
         t_sw = 0
         for i in range(1, t0.shape[0]):
-            d = (identity + t_step*(1-theta)*alpha*M)*rho + beta
+            d = (identity + t_step * (1 - theta) * alpha * M) * rho + beta
             rho_next = linalg.spsolve(A, d)
             # normalization not needed, flux is kept
-            ps = np.trapz(rho.T[mesh.cells < 0],
-                          x=mesh.cells[mesh.cells < 0],
-                          axis=0)
+            ps = np.trapz(rho.T[mesh.cells < 0], x=mesh.cells[mesh.cells < 0], axis=0)
             if t_sw == 0 and ps > 0.5:
-                t_sw = t_step*i
+                t_sw = t_step * i
             # update variable by switching
             rho_next, rho = rho, rho_next
 
     # return
-    return {'t_sw': t_sw,
-            'rho': rho.T,
-            't0': t0,
-            'z0': mesh.cells}
+    return {"t_sw": t_sw, "rho": rho.T, "t0": t0, "z0": mesh.cells}
 
 
 def simple_test():
     """FVM simple test."""
     z_points = 500
     delta = 60
     i0 = 1.5
-    h = 0.
+    h = 0.0
     # time step
     t_step = 0.001
     T = 20
 
-    data = solve_mtj_fp(dim_points=z_points,
-                        delta=delta,
-                        i0=i0,
-                        h=h,
-                        t_step=t_step,
-                        T=T)
-    rho = data['rho']
-    t0 = data['t0']
-    z0 = data['z0']
+    data = solve_mtj_fp(
+        dim_points=z_points, delta=delta, i0=i0, h=h, t_step=t_step, T=T
+    )
+    rho = data["rho"]
+    t0 = data["t0"]
+    z0 = data["z0"]
     _theta_x = np.arccos(z0)
     t_mesh0, z_mesh0 = np.meshgrid(t0, z0)
     fig = plt.figure()
-    ax = fig.add_subplot(111, projection='3d')
+    ax = fig.add_subplot(111, projection="3d")
     # z_mesh0 = np.arccos(z_mesh0)
-    ax.plot_surface(z_mesh0,
-                    t_mesh0,
-                    rho,
-                    alpha=0.7,
-                    cmap='viridis',
-                    edgecolor='k')
+    ax.plot_surface(z_mesh0, t_mesh0, rho, alpha=0.7, cmap="viridis", edgecolor="k")
     plt.show()
 
-    print('plotting 2d')
+    print("plotting 2d")
     fig = plt.figure()
     ax_0 = fig.add_subplot(211)
     ax_1 = fig.add_subplot(212)
@@ -184,12 +171,12 @@ def simple_test():
         t0_idx = np.argmax(t0 >= tt)
         if tt == t0[-1]:
             t0_idx = -1
-        ax_0.plot(_theta_x, rho[:, t0_idx], label=f't={tt}')
-        ax_1.plot(z0, rho[:, t0_idx], label=f't={tt}')
+        ax_0.plot(_theta_x, rho[:, t0_idx], label=f"t={tt}")
+        ax_1.plot(z0, rho[:, t0_idx], label=f"t={tt}")
     ax_0.legend()
-    ax_0.set_yscale('log', base=10)
+    ax_0.set_yscale("log", base=10)
     ax_1.legend()
-    ax_1.set_yscale('log', base=10)
+    ax_1.set_yscale("log", base=10)
     # ax.set_ylim(bottom=1e-10)
     ax_0.grid()
     ax_1.grid()
@@ -204,24 +191,17 @@ def simple_test():
     # plot_3d_evolution(p_imp_0, t0, z0, plot_res=1e6,
     #                   title='ro implicit matmult')
     pt = np.trapz(rho, z0, axis=0)
-    ps = np.trapz(y=rho[z0 < 0], x=z0[z0 < 0], axis=0)/pt
-    pns = np.trapz(y=rho[z0 >= 0], x=z0[z0 >= 0], axis=0)/pt
-
-    ax.plot(t0,
-            ps,
-            '--',
-            color=colors[0],
-            alpha=0.5,
-            label=f'ps i: {i0}')
-    ax.plot(t0,
-            pns,
-            label=f'pns i: {i0}')
-    ax.set_yscale('log', base=10)
+    ps = np.trapz(y=rho[z0 < 0], x=z0[z0 < 0], axis=0) / pt
+    pns = np.trapz(y=rho[z0 >= 0], x=z0[z0 >= 0], axis=0) / pt
+
+    ax.plot(t0, ps, "--", color=colors[0], alpha=0.5, label=f"ps i: {i0}")
+    ax.plot(t0, pns, label=f"pns i: {i0}")
+    ax.set_yscale("log", base=10)
     # ax.set_ylim(bottom=1e-10)
     ax.legend()
     ax.grid()
     plt.show()
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     simple_test()

src/python_compact_model/README.md

@@ -14,8 +14,8 @@ Summary:
 * `test_sllgs_solver.py` shows you the basic s-LLGS solver config, calling (`sllgs_solver.py`)
 * `stochastic_multithread_simulation.py` (calling `sllgs_solver.py`) is the script
 that helps you launching parallel python s-LLGS simulations
-* These results can be compared against Fooker-Plank simulations (see `plot_sllgs_fpe_comparison.py` script)
-* `analytical.py` and `mtj_fp_fvm.py` contain the Fooker-Plank solvers. Analytical contains the WER/RER fitter for the problem optimization
+* These results can be compared against Fokker-Plank simulations (see `plot_sllgs_fpe_comparison.py` script)
+* `analytical.py` and `mtj_fp_fvm.py` contain the Fokker-Plank solvers. Analytical contains the WER/RER fitter for the problem optimization
 * Verilog-a compact models: run the testbenches `tb.scs` and `tb_subckt.scs`
 
 Please, read the full description at [MRAM Framework Description](./doc/README.md).

src/python_compact_model/sllgs_solver.py

@@ -382,11 +382,20 @@ def __init__(self,
         self.l_fl = l_fl
         self._area = np.pi/4*self.w_fl*self.l_fl
         self._volume = self._area*self.t_fl
+        print(f'\t[debug] t_fl    {self.t_fl}')
+        print(f'\t[debug] t_ox    {self.t_ox}')
+        print(f'\t[debug] w_fl    {self.w_fl}')
+        print(f'\t[debug] l_fl    {self.l_fl}')
+        print(f'\t[debug] _area   {self._area}')
+        print(f'\t[debug] _volume {self._volume}')
 
         # magnetic properties
         self.alpha = alpha
         self._U0_gamma_0_alpha = c_U0 * c_gamma_0 / (1 + self.alpha*self.alpha)
+        print(f'U0gamma0alpha: {self._U0_gamma_0_alpha}')
         self.ms = ms
+        print(f'\t[debug] alpha: {self.alpha}')
+        print(f'\t[debug] ms: {self.ms}')
         # Full anisotropy:
         # Watanabe, K.,
         # Shape anisotropy revisited in single-digit nanometer magnetic
@@ -414,8 +423,10 @@ def __init__(self,
                 z_0 = 1
             else:
                 z_0 = -1
-        print(f'z_0: {z_0}')
+        print(f'\t[debug] z_0: {z_0}')
+        print(f'\t[debug] _volume: {self._volume}')
         h_k = 2 * self.k_i * np.abs(z_0) / (self.t_fl * c_U0 * self.ms)
+        print(f'\t[debug] hk_mode: {self.hk_mode}')
         if self.hk_mode == 0:
             h_k += -self.ms * \
                 (self.shape_ani_demag_n[2] -
@@ -573,6 +584,7 @@ def _init_stt_constants(self,
                 _lamb_FL2*np.sqrt((_lamb_PL2-1)/(_lamb_FL2-1))
         elif self.stt_mode == 'stt_oommf_simple':
             self.lambda_s = lambda_s
+            print(f'\t[debug] lambda_s: {self.lambda_s}')
             self.p = p
             self._eps_simple_num = p*lambda_s*lambda_s
             self._eps_simple_den0 = lambda_s*lambda_s + 1
@@ -1223,7 +1235,7 @@ def _f_cart(self, t, m_cart, dt=np.inf):
         h_eff_cart += self.get_h_demag(m_cart)
         h_eff_cart += self.get_h_exchange()
         h_eff_cart += self.get_h_una(m_cart)
-        h_eff_cart += self.get_h_vcma(v_mtj)
+        h_eff_cart += self.get_h_vcma(v_mtj, m_cart)
 
         # wienner process
         # similar to what SPICE solvers do
@@ -1309,7 +1321,7 @@ def solve_ode(self,
             saved_max_step = max_step
 
         # max_step/tstep checks
-        if not scipy_ivp and max_step is None or max_step == np.inf:
+        if not scipy_ivp and (max_step is None or max_step == np.inf):
             print('[error] Max step should be specified when not using '
                   'scipy_ivp mode')
             return None
@@ -1330,8 +1342,11 @@ def solve_ode(self,
             _g = self._g_cart
             y0 = self.m_cart_init
             method_info = ', cartesian coordinates'
+        if scipy_ivp and (max_step == np.inf or max_step is None):
+            save_every = 1
+        else:
+            save_every = int(saved_max_step/max_step)
 
-        save_every = int(saved_max_step/max_step)
         # normalization
         # do time normalization
         if self.solve_normalized:
@@ -1415,6 +1430,7 @@ def solve_ode(self,
             2. * self.alpha * self.temperature * c_KB /
             (c_U0 * th_gamma * self.ms * self._volume))
         v_cart = np.array([_sig, _sig, _sig])
+        print(f'[debug] v_cart: {v_cart}')
 
         # different options
         while t_idx < n_t: