From 15abd99dee05926a02fd29bd548827d375d78ad0 Mon Sep 17 00:00:00 2001 From: Fernando Garcia Redondo Date: Thu, 7 Jul 2022 08:52:16 +0100 Subject: [PATCH 1/4] fix typo --- src/python_compact_model/sllgs_solver.py | 2 +- .../stochastic_multithread_simulation.py | 17 +++++++++++++++++ 2 files changed, 18 insertions(+), 1 deletion(-) diff --git a/src/python_compact_model/sllgs_solver.py b/src/python_compact_model/sllgs_solver.py index 0a0badf..939ee29 100644 --- a/src/python_compact_model/sllgs_solver.py +++ b/src/python_compact_model/sllgs_solver.py @@ -1223,7 +1223,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 diff --git a/src/python_compact_model/stochastic_multithread_simulation.py b/src/python_compact_model/stochastic_multithread_simulation.py index ff099ab..91289aa 100644 --- a/src/python_compact_model/stochastic_multithread_simulation.py +++ b/src/python_compact_model/stochastic_multithread_simulation.py @@ -317,3 +317,20 @@ def run_walks( print('\t[experiment] reproduce results with:') for k in args.__dict__: print('\t\t', k, ': ', args.__dict__[k]) + + +def test(): + + run_walks( + solve_normalized=True, + method='stratonovich_heun', + suffix='rand_walk', + max_step=1e-13, + I0=40e-6, + alpha=0.01, + t_delay=3e-9, + t_pulse=60e-9, + total_th_sims=int(1e2), + cores=2, + plot_walks=True + ) From c9473e6296fe949940227d39fc633890771078ea Mon Sep 17 00:00:00 2001 From: Fernando Garcia Redondo Date: Thu, 7 Jul 2022 08:53:55 +0100 Subject: [PATCH 2/4] Plot sllgs vs fp results --- .../plot_sllgs_fpe_comparison.py | 259 +++++++++--------- 1 file changed, 137 insertions(+), 122 deletions(-) diff --git a/src/sllgs_fpe_comparison/plot_sllgs_fpe_comparison.py b/src/sllgs_fpe_comparison/plot_sllgs_fpe_comparison.py index 48278e9..e6b7d7f 100644 --- a/src/sllgs_fpe_comparison/plot_sllgs_fpe_comparison.py +++ b/src/sllgs_fpe_comparison/plot_sllgs_fpe_comparison.py @@ -26,43 +26,47 @@ # these should match the ones you used for the s-LLGS simulation ###################################################################### # Define a macrospin mesh (i.e. one discretisation cell). -z0 = 1 -initial_m = None # np.array([0.01, 0.01, z0]) # vector in x direction -diam = 50e-9 -t_fl = 1.e-9 -temperature = 300 -k_i_0 = 1.0e-3 -alpha = 0.01 # Gilbert damping -Ms = 1.2e6 # magnetisation saturation (A/m) +Z0 = 1 +INITIAL_M = None # np.array([0.01, 0.01, z0]) # vector in x direction +DIAM = 50e-9 +T_FL = 1.e-9 +TEMP = 300 +K_I_0 = 1.0e-3 +ALPHA = 0.01 # Gilbert damping +MS = 1.2e6 # magnetisation saturation (A/m) # Zeeman (external field) validation -Hext = (0., 0., 0.) # low external magnetic field (A/m) (test anisotropy) +HEXT = (0., 0., 0.) # low external magnetic field (A/m) (test anisotropy) # stt params -eps_prime = 0.0 # [FITTING] constant +EPS_PRIME = 0.0 # [FITTING] constant P = 0.75 -Lambda = 1. +LAMBDA = 1. ################################### -sllgs_files = [ +SLLGS_FILES = [ 'path_to_csv_file_exported_by_stochastic_multithread_simulation.py', ] -sllgs_time_files = [ +SLLGS_TIME_FILES = [ 'path_to_time_csv_file_exported_by_stochastic_multithread_simulation.py', ] -I0_uA = [50]*len(sllgs_files) -t_delays_ns = [5]*len(sllgs_files) +I0_UA = [50]*len(SLLGS_FILES) +T_DELAYS_NS = [5]*len(SLLGS_FILES) -ps_sslgs = [None] * len(t_delays_ns) -time_sslgs = [None] * len(t_delays_ns) -ps_fp = [None] * len(t_delays_ns) -time_fp = [None] * len(t_delays_ns) -compute_analytical_manual = False -plot_intermediate = False -colors = itertools.cycle(mcolors.TABLEAU_COLORS) +ps_sslgs = [None] * len(T_DELAYS_NS) +time_sslgs = [None] * len(T_DELAYS_NS) +ps_fp = [None] * len(T_DELAYS_NS) +time_fp = [None] * len(T_DELAYS_NS) +COMPUTE_ANALYTICAL_MANUAL = False +PLOT_INTERMEDIATE = False + + +COLORS = itertools.cycle(mcolors.TABLEAU_COLORS) +MARKERS = itertools.cycle( + ('o', '*', 'd', '1', ',', '+', '.', 's', 'X', 'x')) def test_h_ext(t): """Set H_ext.""" - return np.array(Hext) + return np.array(HEXT) def test_current(t): @@ -71,30 +75,32 @@ def test_current(t): # get Ic, tau_d etc -llg_o = sllgs_solver.LLG(w_fl=diam, l_fl=diam, t_fl=t_fl, - ms=Ms, - alpha=alpha, - k_i_0=k_i_0, - # thermal_stability_0=thermal_stability_0, - temperature=temperature, - stt_mode='stt_oommf_simple', - p=P, - lambda_s=Lambda, - eps_prime=eps_prime, - # not include temperature noise, - # and do not force its effects on the theta_0 - do_thermal=False, - do_fake_thermal=False, - do_theta_windowing=False, - m_init=initial_m, - # shape_ani_demag_mode=0, - # theta_init=0.09, - theta_pl=0.0, # [rad] pinned layer theta - phi_pl=0.0, # [rad] pinned layer phi - h_ext_cart=test_h_ext, - i_amp_fn=test_current) - -print(f'P: {P}, P/2: {P/2}, eps: {llg_o.get_epsilon_stt(1)}') +LLG_O = sllgs_solver.LLG(w_fl=DIAM, l_fl=DIAM, t_fl=T_FL, + ms=MS, + alpha=ALPHA, + k_i_0=K_I_0, + # thermal_stability_0=thermal_stability_0, + temperature=TEMP, + stt_mode='stt_oommf_simple', + p=P, + lambda_s=LAMBDA, + eps_prime=EPS_PRIME, + # not include temperature noise, + # and do not force its effects on the theta_0 + do_thermal=False, + do_fake_thermal=False, + do_theta_windowing=False, + m_init=INITIAL_M, + # shape_ani_demag_mode=0, + # theta_init=0.09, + # [rad] pinned layer theta + theta_pl=0.0, + # [rad] pinned layer phi + phi_pl=0.0, + h_ext_cart=test_h_ext, + i_amp_fn=test_current) + +print(f'P: {P}, P/2: {P/2}, eps: {LLG_O.get_epsilon_stt(1)}') ############################################################################## # importing raw data @@ -102,10 +108,12 @@ def test_current(t): ############################################################################## -def analyze_s_llgs_tolerances(data_file, time_file, t_delay, I0_uA, - compute_analytical_manual=False, - plot_intermediate=True, - dim_points=1000): +def analyze_s_llgs_tolerances( + data_file, time_file, t_delay, I0_uA, + llg_o=LLG_O, + compute_analytical_manual=False, + plot_intermediate=True, + dim_points=1000): """Analyze s-LLGS.""" # get all data, note it is not delayed anymore sllgs_time, sllgs_data, sllgs_sw_idx, sllgs_sw_ps = sllgs_i.process_sllgs_data( @@ -121,7 +129,7 @@ def analyze_s_llgs_tolerances(data_file, time_file, t_delay, I0_uA, # } ######################################################################## # Initializations - I0 = np.sign(z0)*I0_uA*1e-6 + I0 = np.sign(Z0)*I0_uA*1e-6 lin_space_z = False L0_max = 200 @@ -154,12 +162,13 @@ def analyze_s_llgs_tolerances(data_file, time_file, t_delay, I0_uA, print(f't_final: {xlim}') print(f'tau: {tau} i: {i}') lin_space_z = False - rho_0_at_pi = z0 < 0 + rho_0_at_pi = llg_o.theta_init > np.pi/2 + # sllgs_data has (data, time) s_rho_0_fit = np.polynomial.legendre.legfit( x=np.cos(theta_axis), # sllgs is normalized to 1 already - y=sllgs_data[:, 0][::-1], + y=sllgs_data[0][:, 0][::-1], deg=L0_max) _, rho_0_a, s_rho_0_a = analytical.get_state_rho_0( delta=delta, @@ -192,7 +201,7 @@ def analyze_s_llgs_tolerances(data_file, time_file, t_delay, I0_uA, fig, axs = plt.subplots(1, 1, figsize=(8, 6)) axs.plot(theta_axis, - sllgs_data[:, 0][::-1], + sllgs_data[0][:, 0][::-1], '+-', label='sllgs') axs.plot(fp_theta, @@ -279,10 +288,10 @@ def analyze_s_llgs_tolerances(data_file, time_file, t_delay, I0_uA, # cmap = plt.cm.PuBu_r # s-LLGS - v_min = np.max([1e-8, np.min(sllgs_data)]) - v_max = np.max(sllgs_data) + v_min = np.max([1e-8, np.min(sllgs_data[0])]) + v_max = np.max(sllgs_data[0]) - sllgs_data[sllgs_data < v_min] = v_min + sllgs_data[0][sllgs_data[0] < v_min] = v_min fp_data_a[fp_data_a < v_min] = v_min fp_data_b[fp_data_b < v_min] = v_min @@ -290,7 +299,7 @@ def analyze_s_llgs_tolerances(data_file, time_file, t_delay, I0_uA, ax0 = axs[0] ax0.pcolormesh(1e9*sllgs_time, np.pi-fp_theta, - sllgs_data, + sllgs_data[0], # vmin=v_min, norm=mcolors.LogNorm(vmin=v_min, vmax=v_max), @@ -417,73 +426,79 @@ def analyze_s_llgs_tolerances(data_file, time_file, t_delay, I0_uA, 'fp_ps_manual': fp_a_ps_manual, 'sllgs_time': sllgs_time, 'sllgs_theta': theta_axis, - 'sllgs_data': sllgs_data, + 'sllgs_data': sllgs_data[0], 'sllgs_ps': sllgs_sw_ps } -################################ -# comparison -################################ - -if compute_analytical_manual: - ps_fp_manual = [None] * len(t_delays_ns) - -plot_idx = 3 -for idx in range(len(sllgs_files)): - # if idx%4 != plot_idx: - # continue - print('---------------------', idx) - data = analyze_s_llgs_tolerances( - sllgs_files[idx], - sllgs_time_files[idx], - t_delays_ns[idx]*1e-9, - I0_uA[idx], - compute_analytical_manual=compute_analytical_manual, - plot_intermediate=plot_intermediate) - ps_sslgs[idx] = data['sllgs_ps'] - time_sslgs[idx] = data['sllgs_time'] - ps_fp[idx] = data['fp_ps'] - time_fp[idx] = data['fp_time'] - if compute_analytical_manual: - ps_fp_manual[idx] = data['fp_ps_manual'] - -fig, ax = plt.subplots(1, 1, - sharex=True, - figsize=(8, 6)) -markers = itertools.cycle( - ('o', '*', 'd', '1', ',', '+', '.', 's', 'X', 'x')) -for idx in range(len(sllgs_files)): - # if idx%4 != plot_idx: - # continue - print('---------------------', idx) - color = next(colors) - marker = next(markers) - np.savetxt(f'dummy_sllgs_wer_{idx}.csv', - np.array([time_sslgs[idx], 1-ps_sslgs[idx]])) - ax.plot(1e9*time_sslgs[idx], 1-ps_sslgs[idx], - # label=sllgs_files[idx].split('/')[-1], - label='10000 s-LLGS walks', - marker=marker, - color=color) - ax.plot(1e9*time_fp[idx], - 1-ps_fp[idx], - '--', - # label=f'fp {sllgs_files[idx]}', - label='Fokker-Plank', - color=color) +def plot_results(sllgs_files=SLLGS_FILES, + sllgs_time_files=SLLGS_TIME_FILES, + t_delays_ns=T_DELAYS_NS, + plot_intermediate=PLOT_INTERMEDIATE, + I0s_uA=I0_UA, + llg_os=LLG_O, + compute_analytical_manual=COMPUTE_ANALYTICAL_MANUAL, + ): + """Plot results for comparison.""" if compute_analytical_manual: + ps_fp_manual = [None] * len(t_delays_ns) + + for idx in range(len(sllgs_files)): + print('---------------------', idx) + data = analyze_s_llgs_tolerances( + data_file=sllgs_files[idx], + time_file=sllgs_time_files[idx], + t_delay=t_delays_ns[idx]*1e-9, + I0_uA=I0s_uA[idx], + llg_o=llg_os[idx], + compute_analytical_manual=compute_analytical_manual, + plot_intermediate=plot_intermediate) + ps_sslgs[idx] = data['sllgs_ps'] + time_sslgs[idx] = data['sllgs_time'] + ps_fp[idx] = data['fp_ps'] + time_fp[idx] = data['fp_time'] + if compute_analytical_manual: + ps_fp_manual[idx] = data['fp_ps_manual'] + + fig, ax = plt.subplots(1, 1, + sharex=True, + figsize=(8, 6)) + for idx in range(len(sllgs_files)): + print('---------------------', idx) + color = next(COLORS) + marker = next(MARKERS) + np.savetxt(f'dummy_sllgs_wer_{idx}.csv', + np.array([time_sslgs[idx], 1-ps_sslgs[idx]])) + ax.plot(1e9*time_sslgs[idx], 1-ps_sslgs[idx], + # label=sllgs_files[idx].split('/')[-1], + label='10000 s-LLGS walks', + marker=marker, + color=color) ax.plot(1e9*time_fp[idx], - 1-ps_fp_manual[idx], - ':', - # label=sllgs_files[idx], + 1-ps_fp[idx], + '--', + # label=f'fp {sllgs_files[idx]}', + label='Fokker-Plank', color=color) -# ax.plot(1e9*fp_time, 1-fp_b_ps, label='fpb') -# ax.plot(1e9*time, 1-sw_ps, '--', label='sllgs') -ax.set_title('WER over time') -ax.set_xlabel('time [a.u.]') -ax.set_ylabel('WER') -ax.set_yscale('log', base=10) -ax.legend() -ax.grid() -plt.show() + if compute_analytical_manual: + ax.plot(1e9*time_fp[idx], + 1-ps_fp_manual[idx], + ':', + # label=sllgs_files[idx], + color=color) + # ax.plot(1e9*fp_time, 1-fp_b_ps, label='fpb') + # ax.plot(1e9*time, 1-sw_ps, '--', label='sllgs') + ax.set_title('WER over time') + ax.set_xlabel('time [a.u.]') + ax.set_ylabel('WER') + ax.set_yscale('log', base=10) + ax.legend() + ax.grid() + plt.show() + + +################################ +# comparison +################################ +if __name__ == "__main__": + plot_results() From a9bc2aa8f885163266a534589c5fa067eec50090 Mon Sep 17 00:00:00 2001 From: Fernando Garcia Redondo Date: Wed, 18 Feb 2026 08:38:37 +0000 Subject: [PATCH 3/4] Sync with internal repo Including improvements from https://github.com/Fratorhe/FVM/ --- src/README.md | 4 +- src/fokker_plank/README.md | 2 +- src/fokker_plank/fvm/fvm_lib/fvm_classes.py | 2 +- src/fokker_plank/fvm/mtj_fp_fvm.py | 144 ++++++++---------- src/python_compact_model/README.md | 4 +- src/python_compact_model/sllgs_solver.py | 22 ++- .../stochastic_multithread_simulation.py | 65 +++++++- .../plot_sllgs_fpe_comparison.py | 34 +++-- .../mram_lib/mtj_subcircuit.scs | 2 +- 9 files changed, 167 insertions(+), 112 deletions(-) diff --git a/src/README.md b/src/README.md index 8310c1d..a669412 100644 --- a/src/README.md +++ b/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). diff --git a/src/fokker_plank/README.md b/src/fokker_plank/README.md index 83d7baf..48d8e0e 100644 --- a/src/fokker_plank/README.md +++ b/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). diff --git a/src/fokker_plank/fvm/fvm_lib/fvm_classes.py b/src/fokker_plank/fvm/fvm_lib/fvm_classes.py index 4b69550..869454b 100644 --- a/src/fokker_plank/fvm/fvm_lib/fvm_classes.py +++ b/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 diff --git a/src/fokker_plank/fvm/mtj_fp_fvm.py b/src/fokker_plank/fvm/mtj_fp_fvm.py index 0a6dbfd..694d85e 100644 --- a/src/fokker_plank/fvm/mtj_fp_fvm.py +++ b/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,71 +73,68 @@ 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(): @@ -147,34 +142,26 @@ def 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() diff --git a/src/python_compact_model/README.md b/src/python_compact_model/README.md index 56875dd..a5d1c61 100644 --- a/src/python_compact_model/README.md +++ b/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). diff --git a/src/python_compact_model/sllgs_solver.py b/src/python_compact_model/sllgs_solver.py index 939ee29..2de0adf 100644 --- a/src/python_compact_model/sllgs_solver.py +++ b/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 @@ -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: diff --git a/src/python_compact_model/stochastic_multithread_simulation.py b/src/python_compact_model/stochastic_multithread_simulation.py index 91289aa..726b978 100644 --- a/src/python_compact_model/stochastic_multithread_simulation.py +++ b/src/python_compact_model/stochastic_multithread_simulation.py @@ -132,6 +132,7 @@ def test_current(t): h_ext_cart=test_h_ext, i_amp_fn=test_current, seed=idx+2000) + llg_b._get_Ic(debug=True) llg_sol = llg_b.solve_and_plot(final_t=t_final, max_step=max_step, solve_spherical=False, @@ -159,6 +160,7 @@ def run_walks( alpha, max_step, t_delay, + initial_idx=0, t_pulse=60e-9, total_th_sims=10000, cores=32, @@ -172,7 +174,7 @@ def run_walks( results_a = [pool.apply_async( worker_stt_simplest, args=( - idx, + initial_idx+idx, solve_normalized, method, I0, @@ -205,11 +207,11 @@ def run_walks( print('[info][parallel_solver] all threads computed') np.savetxt( - f'{suffix}_{t_delay}_del_{I0}A_interp_{total_th_sims}_sims.csv', + f'{suffix}_{t_delay}_del_{I0}A_interp_{total_th_sims}_sims_{initial_idx}.csv', np.array(theta_i_list)) np.savetxt( f'{suffix}_{t_delay}_del_{I0}A_time_interp_' - f'{total_th_sims}_sims.csv', + f'{total_th_sims}_sims_{initial_idx}.csv', t_i) if not plot_walks: @@ -291,6 +293,11 @@ def run_walks( type=str2bool, default=False, help='Plot distribution') + parser.add_argument( + '--initial_idx', + type=int, + default=int(0), + help='initial idx (for seeding)') args = parser.parse_args() @@ -310,6 +317,7 @@ def run_walks( t_delay=args.t_delay, t_pulse=args.t_pulse, total_th_sims=args.total_th_sims, + initial_idx=args.initial_idx, cores=args.cores, plot_walks=args.plot_walks ) @@ -324,13 +332,56 @@ def test(): run_walks( solve_normalized=True, method='stratonovich_heun', - suffix='rand_walk', + suffix='rand_walk_1Em13', max_step=1e-13, I0=40e-6, alpha=0.01, t_delay=3e-9, - t_pulse=60e-9, - total_th_sims=int(1e2), - cores=2, + t_pulse=15e-9, + total_th_sims=int(1), + initial_idx=1000, + cores=19, + plot_walks=True + ) + +def test_tstep(): + + run_walks( + solve_normalized=True, + method='stratonovich_heun', + suffix='rand_walk_1Em14', + max_step=1e-14, + I0=40e-6, + alpha=0.01, + t_delay=3e-9, + t_pulse=12e-9, + total_th_sims=int(1000), + cores=19, + plot_walks=True + ) + run_walks( + solve_normalized=True, + method='stratonovich_heun', + suffix='rand_walk_0p5Em13', + max_step=0.5e-13, + I0=40e-6, + alpha=0.01, + t_delay=3e-9, + t_pulse=12e-9, + total_th_sims=int(1000), + cores=19, + plot_walks=True + ) + run_walks( + solve_normalized=True, + method='stratonovich_heun', + suffix='rand_walk_1Em13', + max_step=1e-13, + I0=40e-6, + alpha=0.01, + t_delay=3e-9, + t_pulse=12e-9, + total_th_sims=int(1000), + cores=19, plot_walks=True ) diff --git a/src/sllgs_fpe_comparison/plot_sllgs_fpe_comparison.py b/src/sllgs_fpe_comparison/plot_sllgs_fpe_comparison.py index e6b7d7f..d32125e 100644 --- a/src/sllgs_fpe_comparison/plot_sllgs_fpe_comparison.py +++ b/src/sllgs_fpe_comparison/plot_sllgs_fpe_comparison.py @@ -17,7 +17,7 @@ import itertools import matplotlib.colors as mcolors -import fooker_plank.analytical as analytical +import fokker_plank.analytical.analytical as analytical import python_compact_model.sllgs_solver as sllgs_solver import sllgs_importer as sllgs_i @@ -51,17 +51,13 @@ I0_UA = [50]*len(SLLGS_FILES) T_DELAYS_NS = [5]*len(SLLGS_FILES) -ps_sslgs = [None] * len(T_DELAYS_NS) -time_sslgs = [None] * len(T_DELAYS_NS) -ps_fp = [None] * len(T_DELAYS_NS) -time_fp = [None] * len(T_DELAYS_NS) COMPUTE_ANALYTICAL_MANUAL = False PLOT_INTERMEDIATE = False COLORS = itertools.cycle(mcolors.TABLEAU_COLORS) MARKERS = itertools.cycle( - ('o', '*', 'd', '1', ',', '+', '.', 's', 'X', 'x')) + ('o', '*', 'd', '1', ',', '+', '.', 's', 'X', 'x')) def test_h_ext(t): @@ -116,8 +112,9 @@ def analyze_s_llgs_tolerances( dim_points=1000): """Analyze s-LLGS.""" # get all data, note it is not delayed anymore - sllgs_time, sllgs_data, sllgs_sw_idx, sllgs_sw_ps = sllgs_i.process_sllgs_data( - data_file, time_file, t_delay, dim_points) + (sllgs_time, sllgs_data, + sllgs_sw_idx, sllgs_sw_ps) = sllgs_i.process_sllgs_data( + data_file, time_file, 0, dim_points) theta_axis = np.linspace(np.pi, 0, dim_points) # debug sllgs @@ -164,11 +161,16 @@ def analyze_s_llgs_tolerances( lin_space_z = False rho_0_at_pi = llg_o.theta_init > np.pi/2 + rho_0_time_idx = np.searchsorted(sllgs_time, t_delay+t_delay/50) + print(f'\n->fitting rho at time {sllgs_time[rho_0_time_idx]}\n') + print(f'\n->mean: {np.mean(sllgs_data[0][:, rho_0_time_idx])}') + print(f'\n->std: {np.std(sllgs_data[0][:, rho_0_time_idx])}') + # sllgs_data has (data, time) s_rho_0_fit = np.polynomial.legendre.legfit( x=np.cos(theta_axis), # sllgs is normalized to 1 already - y=sllgs_data[0][:, 0][::-1], + y=sllgs_data[0][:, rho_0_time_idx][::-1], deg=L0_max) _, rho_0_a, s_rho_0_a = analytical.get_state_rho_0( delta=delta, @@ -201,7 +203,7 @@ def analyze_s_llgs_tolerances( fig, axs = plt.subplots(1, 1, figsize=(8, 6)) axs.plot(theta_axis, - sllgs_data[0][:, 0][::-1], + sllgs_data[0][:, rho_0_time_idx][::-1], '+-', label='sllgs') axs.plot(fp_theta, @@ -443,6 +445,11 @@ def plot_results(sllgs_files=SLLGS_FILES, if compute_analytical_manual: ps_fp_manual = [None] * len(t_delays_ns) + ps_sslgs = [None] * len(t_delays_ns) + time_sslgs = [None] * len(t_delays_ns) + ps_fp = [None] * len(t_delays_ns) + time_fp = [None] * len(t_delays_ns) + for idx in range(len(sllgs_files)): print('---------------------', idx) data = analyze_s_llgs_tolerances( @@ -469,9 +476,9 @@ def plot_results(sllgs_files=SLLGS_FILES, marker = next(MARKERS) np.savetxt(f'dummy_sllgs_wer_{idx}.csv', np.array([time_sslgs[idx], 1-ps_sslgs[idx]])) - ax.plot(1e9*time_sslgs[idx], 1-ps_sslgs[idx], - # label=sllgs_files[idx].split('/')[-1], - label='10000 s-LLGS walks', + ax.plot(1e9*(time_sslgs[idx])-t_delays_ns[idx], 1-ps_sslgs[idx], + label=sllgs_files[idx].split('/')[-1], + # label='10000 s-LLGS walks', marker=marker, color=color) ax.plot(1e9*time_fp[idx], @@ -494,6 +501,7 @@ def plot_results(sllgs_files=SLLGS_FILES, ax.set_yscale('log', base=10) ax.legend() ax.grid() + # plt.savefig('wer.png') plt.show() diff --git a/src/verilog_a_compact_model/mram_lib/mtj_subcircuit.scs b/src/verilog_a_compact_model/mram_lib/mtj_subcircuit.scs index 3792d35..3c6f613 100644 --- a/src/verilog_a_compact_model/mram_lib/mtj_subcircuit.scs +++ b/src/verilog_a_compact_model/mram_lib/mtj_subcircuit.scs @@ -34,7 +34,7 @@ v_aux (aux_1 0) bsource v=v(fl, pl) r_aux (aux_0 0) resistor r=1 // simple TMR based conduction -r_mtj_module (fl pl) bsource r=p_r_p*(1 + p_tmr/(p_tmr + 2))/(1 - p_tmr/(p_tmr + 2)*V(mz)) +r_mtj_module (fl pl) bsource r=p_r_p*(1 + p_tmr/(p_tmr + 2))/(1 + p_tmr/(p_tmr + 2)*V(mz)) ends mtj_subcircuit From 4d21f22e3ec63f3d15b2b081530c6ee4aebae8bc Mon Sep 17 00:00:00 2001 From: Fernando Garcia Redondo Date: Wed, 18 Feb 2026 08:40:49 +0000 Subject: [PATCH 4/4] Readme --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 3544de0..7a28fe3 100644 --- a/README.md +++ b/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