From e3381c21155589505ec6464d29426208fc9bea6e Mon Sep 17 00:00:00 2001 From: LemurPwned Date: Sun, 15 Sep 2024 15:12:45 +0200 Subject: [PATCH 1/2] updates to the domain -- fixing the user upload, adding help and better start params --- view/domain.py | 182 ++++++++++++++++++++++++++++++----------- view/simulation_fns.py | 17 ++-- view/streamlit_app.py | 139 ++++++++++++++++++------------- view/utils.py | 11 +++ 4 files changed, 239 insertions(+), 110 deletions(-) diff --git a/view/domain.py b/view/domain.py index c2f6ffc..8574b59 100644 --- a/view/domain.py +++ b/view/domain.py @@ -5,22 +5,33 @@ import streamlit as st from helpers import read_mh_data from simulation_fns import create_single_domain, get_axis_angles -from utils import GENERIC_BOUNDS, GENERIC_UNITS +from utils import GENERIC_BOUNDS, GENERIC_UNITS, get_init_kval from cmtj import AxialDriver, Junction from cmtj.utils import FieldScan _lock = RLock() + with st.container(): st.write("# Domain fitting") - - st.write( - "Fit M(H) to multidomain model. " - "Upload your file with experimental data: with columns H, mx, my, mz.\n" - "First column is always H in A/m, the rest are the components of the magnetisation in range (-1, 1)." - "If you upload just two columns, the script will assume that the data is in the form H, mx." - ) + with st.expander("Explanation"): + st.write( + "#### Fit M(H) to multidomain model. \n" + "##### Experimental data\n" + "Upload your file with experimental data: with columns H, mx, my, mz.\n" + "First column is always H in A/m, the rest are the components of the magnetisation in range (-1, 1)." + "If you upload just two columns, the script will assume that the data is in the form H, mx.\n" + "##### Simulation\n" + "Add new domains in the sidebar panel. Area of each domain is a weight. The resulting $m_\mathrm{mixed}(h)$ plot is produces with\n\n" + r"$m_\mathrm{mixed} = (1/{\sum_i a_i})\sum_i a_i m_i$" + "\n\n where $a_i$ is the area of $i$th domain and $m_i$ is $m(h)$ for that domain.\n" + "##### Coordinate system\n" + r"$\theta$ is the polar angle and $\phi$ is the azimuthal angle." + r" $\theta=90^\circ$ denotes fully in-plane magnetisation, " + r"$\theta=0^\circ$ denotes out-of-plane magnetisation (positive z-direction)" + r"$\phi=0^\circ$ denotes magnetisation pointing in the positive x direction." + ) progress_bar = st.progress(0) @@ -46,7 +57,8 @@ def simulate_domains(): weights = [st.session_state[f"area{i}"] * 1e-18 for i in range(st.session_state.N)] wsum = sum(weights) weights = [w / wsum for w in weights] - Mmixed = np.zeros((len(Hscan), 3)) + Mmixed = np.zeros((len(Hscan), 3), dtype=np.float16) + Mdomains = np.zeros((st.session_state.N, len(Hscan), 3), np.float16) for i, H in enumerate(Hvecs): j.clearLog() j.setLayerExternalFieldDriver("all", AxialDriver(*H)) @@ -60,11 +72,13 @@ def simulate_domains(): mx = np.mean(log[f"domain_{l}_mx"][-st.session_state.last_N :]) my = np.mean(log[f"domain_{l}_my"][-st.session_state.last_N :]) mz = np.mean(log[f"domain_{l}_mz"][-st.session_state.last_N :]) - Mmixed[i] += np.array([mx, my, mz]) * weights[l] + Mdomains[l, i] = np.array([mx, my, mz]) + Mmixed[i] += Mdomains[l, i] * weights[l] # for each layer take last N values progress_bar.progress(int((i / maxlen) * 100)) st.session_state["Mmixed"] = Mmixed st.session_state["Hscan"] = Hscan + st.session_state["Mdomains"] = Mdomains with st.sidebar: @@ -77,6 +91,52 @@ def simulate_domains(): accept_multiple_files=False, key="upload", ) + st.checkbox("Show domains in polar angles", key="domain_use_angle", value=False) + domain_params = st.expander("Domain parameters", expanded=True) + with domain_params: + N = st.number_input( + "Number of domains", min_value=1, max_value=10, value=1, key="N" + ) + for i in range(N): + st.markdown(f"#### Domain {i+1}") + unit_k = GENERIC_UNITS["K"] + st.number_input( + f"K ({i+1}) " r" ($\mathrm{" f"{unit_k}" "}$)", + min_value=GENERIC_BOUNDS["K"][0], + max_value=GENERIC_BOUNDS["K"][1], + value=1.2e3, + step=10.0, + key=f"K{i}", + help="Uniaxial anisotropy constant of the domain. Does not include shape anisotropy.", + ) + + st.number_input( + f"area ({i+1}) " r"($\mathrm{nm^2}$)", + min_value=1.0, + max_value=500.0, + value=10.0, + key=f"area{i}", + help="Area of the domain. In fact, this is a weight of the domain.", + ) + + st.number_input( + r"$\theta_\mathrm{K}$ (deg.)", + min_value=0.0, + max_value=180.0, + value=get_init_kval(i), + key=f"theta_K{i}", + help="Polar angle of the anisotropy axis", + ) + st.number_input( + r"$\phi_\mathrm{K}$ (deg.)", + min_value=0.0, + max_value=180.0, + value=get_init_kval(i), + key=f"phi_K{i}", + help="Azimuthal angle of the anisotropy axis", + ) + st.markdown("-----\n") + with st.expander("Simulation parameters", expanded=True): st.selectbox( "H axis", options=["x", "y", "z", "xy", "xz", "yz"], key="H_axis", index=0 @@ -125,6 +185,7 @@ def simulate_domains(): value=1.8, step=0.01, key="Ms_shared", + help="Saturation magnetisation of the system. Affects the shape anisotropy", ) st.number_input( "thickness (nm)", @@ -141,40 +202,9 @@ def simulate_domains(): key="alpha_shared", format="%.3f", ) - domain_params = st.expander("Domain parameters", expanded=True) - with domain_params: - N = st.number_input( - "Number of domains", min_value=1, max_value=10, value=1, key="N" - ) - for i in range(N): - st.markdown(f"#### Domain {i+1}") - unit_k = GENERIC_UNITS["K"] - st.number_input( - f"K ({i+1}) " r" ($\mathrm{" f"{unit_k}" "}$)", - min_value=GENERIC_BOUNDS["K"][0], - max_value=GENERIC_BOUNDS["K"][1], - value=1.2e3, - step=10.0, - key=f"K{i}", - ) - st.number_input( - f"area ({i+1}) " r"($\mathrm{nm^2}$)", - min_value=1.0, - max_value=500.0, - value=10.0, - key=f"area{i}", - ) - st.radio( - f"anisotropy axis ({i+1})", - options=["x", "y", "z"], - key=f"anisotropy_axis{i}", - index=2, - ) - st.markdown("-----\n") - -def render(Hscan, Mmixed): +def render(Hscan, Mmixed, Mdomains=None): if len(Hscan) <= 0 or len(Mmixed) <= 0: return with _lock: @@ -187,21 +217,35 @@ def render(Hscan, Mmixed): fields, mh = read_mh_data() render_from_exp(ax, fields=fields, mh=mh) render_from_sim(ax, Hscan, Mmixed) - + fig.suptitle("Mixed Domain model") + st.pyplot(fig) + n = 3 + if st.session_state["domain_use_angle"]: + n = 2 + w, h = plt.figaspect(st.session_state.N / n) + fig, ax = plt.subplots( + st.session_state.N, n, dpi=400, figsize=(w, h), sharex=True, sharey=True + ) + if st.session_state.N <= 1: + ax = np.array([ax]) + try: + render_individual_domains(fig, ax, Hscan, Mdomains) + except IndexError: + pass st.pyplot(fig) def render_from_exp(ax, fields, mh): if len(fields) <= 0 or len(mh) <= 0: st.toast("No data to render") - + fields = np.asarray(fields) m = np.asarray(mh) for i, l in zip(range(m.shape[1]), "xyz"): ax[i].plot( - fields, + fields / 1e3, # A/m --> kA/m m[:, i], label=f"$m_{l}$", - color="black", + color="yellow", marker="x", alpha=0.5, linestyle="--", @@ -221,6 +265,46 @@ def render_from_sim(ax, Hscan, Mmixed): ax[2].set_xlabel("H (kA/m)") +def render_individual_domains(fig, ax, Hscan, M: list[list[float]]): + n = 3 + if st.session_state["domain_use_angle"]: + n = 2 + ax[0, 0].set_title(r"$\theta$") + ax[0, 1].set_title(r"$\phi$") + for i, domain_m in enumerate(M): + # convert to polar + theta_m, phi_m, _ = FieldScan.vector2angle( + domain_m[:, 0], domain_m[:, 1], domain_m[:, 2] + ) + ax[i, 0].plot( + Hscan / 1e3, theta_m, label=rf"$\theta_{i+1}$", color="crimson" + ) + ax[i, 1].plot( + Hscan / 1e3, phi_m, label=rf"$\phi_{i+1}$", color="forestgreen" + ) + ax[i, 0].set_ylabel(rf"$\theta_{i+1}$ (deg.)") + ax[i, 1].set_ylabel(rf"$\phi_{i+1}$ (deg.)") + else: + n = 3 + ax[0, 0].set_title(r"$m_\mathrm{x}$") + ax[0, 1].set_title(r"$m_\mathrm{y}$") + ax[0, 2].set_title(r"$m_\mathrm{z}$") + for i, domain_m in enumerate(M): + ax[i, 0].plot( + Hscan / 1e3, domain_m[:, 0], label=rf"$m_{i+1}$", color="crimson" + ) + ax[i, 1].plot( + Hscan / 1e3, domain_m[:, 1], label=rf"$m_{i+1}$", color="forestgreen" + ) + ax[i, 2].plot( + Hscan / 1e3, domain_m[:, 2], label=rf"$m_{i+1}$", color="royalblue" + ) + ax[i, 0].set_ylabel(rf"$m_{i+1}$") + for j in range(n): + ax[-1, j].set_xlabel("H (kA/m)") + fig.suptitle("Individual Domains") # fig.legend() + + simulate_btn = st.button( "Simulate", key="simulate", @@ -228,4 +312,8 @@ def render_from_sim(ax, Hscan, Mmixed): type="primary", ) -render(st.session_state.get("Hscan", []), st.session_state.get("Mmixed", [])) +render( + st.session_state.get("Hscan", []), + st.session_state.get("Mmixed", []), + Mdomains=st.session_state.get("Mdomains", []), +) diff --git a/view/simulation_fns.py b/view/simulation_fns.py index 70516e4..0be0c48 100644 --- a/view/simulation_fns.py +++ b/view/simulation_fns.py @@ -14,11 +14,14 @@ def create_single_domain(id_: str) -> Layer: demag = [CVector(0, 0, 0), CVector(0, 0, 0), CVector(0, 0, 1)] - Kdir1 = get_axis_cvector(st.session_state[f"anisotropy_axis{id_}"]) + # Kdir1 = get_axis_cvector(st.session_state[f"anisotropy_axis{id_}"]) + Kdir = FieldScan.angle2vector( + theta=st.session_state[f"theta_K{id_}"], phi=st.session_state[f"phi_K{id_}"] + ) layer = Layer( id=f"domain_{id_}", - mag=Kdir1, - anis=Kdir1, + mag=Kdir, + anis=Kdir, Ms=st.session_state["Ms_shared"], thickness=st.session_state["thickness_shared"] * 1e-9, cellSurface=1e-16, @@ -34,11 +37,13 @@ def create_single_domain(id_: str) -> Layer: def create_single_layer(id_: str) -> tuple: """Do not forget to rescale the units!""" demag = [CVector(0, 0, 0), CVector(0, 0, 0), CVector(0, 0, 1)] - Kdir1 = get_axis_cvector(st.session_state[f"anisotropy_axis{id_}"]) + Kdir = FieldScan.angle2vector( + theta=st.session_state[f"theta_K{id_}"], phi=st.session_state[f"phi_K{id_}"] + ) layer = Layer( id=f"layer_{id_}", - mag=Kdir1, - anis=Kdir1, + mag=Kdir, + anis=Kdir, Ms=st.session_state[f"Ms{id_}"], thickness=st.session_state[f"thickness{id_}"] * 1e-9, cellSurface=1e-16, diff --git a/view/streamlit_app.py b/view/streamlit_app.py index c43c30d..08092c4 100644 --- a/view/streamlit_app.py +++ b/view/streamlit_app.py @@ -6,29 +6,24 @@ from helpers import simulate_pimm, simulate_vsd from utils import GENERIC_BOUNDS, GENERIC_UNITS -apptitle = "CMTJ simulator" - -st.title(apptitle) - -container = st.container() - -container.markdown( - """ - ## Data Upload - If you want to upload data, to overlay it on the plot, please upload - a file with \t separated columns: H, f1, f2, ... Put H in (A/m) and f in (Hz). - They will be rescaled to (kA/m) and (GHz) automatically. - """ -) -container.file_uploader( - "Upload your data here", - help="Upload your data here. Must be `\t` separated values and have H and f headers.", - type=["txt", "dat"], - accept_multiple_files=False, - key="upload", -) +with st.expander("# Read me"): + st.write( + """ + ### Data Upload + If you want to upload data, to overlay it on the plot, please upload + a file with \t separated columns: H, f1, f2, ... Put H in (A/m) and f in (Hz). + They will be rescaled to (kA/m) and (GHz) automatically. + """ + ) with st.sidebar: + st.file_uploader( + "Upload your data here", + help="Upload your data here. Must be `\t` separated values and have H and f headers.", + type=["txt", "dat"], + accept_multiple_files=False, + key="upload", + ) N = st.number_input( "Number of layers", min_value=1, max_value=10, value=1, key="N", format="%d" ) @@ -42,6 +37,7 @@ value=0.52, step=0.01, key=f"Ms{i}", + help="Saturation magnetisation (T). Affectes shape aniostropy.", ) st.number_input( f"K ({i+1}) ({GENERIC_UNITS['K']})", @@ -58,6 +54,7 @@ value=1e-3, key=f"alpha{i}", format="%.3f", + help="Gilbert damping parameter", ) st.number_input( f"thickness ({i+1}) (nm)", @@ -72,6 +69,7 @@ max_value=500.0, value=10.0, key=f"width{i}", + help="Width of the sample. For resistance calculations", ) st.number_input( f"length ({i+1}) (um)", @@ -79,12 +77,23 @@ max_value=500.0, value=10.0, key=f"length{i}", + help="Length of the sample. For resistance calculations", ) - st.radio( - f"anisotropy axis ({i+1})", - options=["x", "y", "z"], - key=f"anisotropy_axis{i}", - index=2, + st.number_input( + r"$\theta_\mathrm{K}$ (deg.)", + min_value=0.0, + max_value=180.0, + value=0.0, + key=f"theta_K{i}", + help="Polar angle of the anisotropy axis", + ) + st.number_input( + r"$\phi_\mathrm{K}$ (deg.)", + min_value=0.0, + max_value=180.0, + value=0.0, + key=f"phi_K{i}", + help="Azimuthal angle of the anisotropy axis", ) st.markdown("-----\n") @@ -97,6 +106,7 @@ value=0.0, key=f"J{j}", format="%.3f", + help="Interlayer exchange coupling constant", ) with st.expander("Simulation & control parameters"): st.selectbox( @@ -118,6 +128,7 @@ value=1e-12, key="int_step", format="%.1e", + help="Integration step for the simulation", ) st.number_input( "sim_time (ns)", @@ -126,6 +137,7 @@ value=16, key="sim_time", format="%d", + help="Simulation time in nanoseconds", ) st.number_input( "max_freq (GHz)", @@ -134,32 +146,34 @@ value=50, key="max_freq", format="%d", + help="Maximum frequency (cutoff) visible in plot", ) pimm_tab, vsd_tab, opt_tab = st.tabs(["PIMM", "VSD", "Optimization"]) with opt_tab: - st.markdown( - """ - ## Optimization + with st.expander("Optimization parameters"): + st.markdown( + """ + ### Optimization - Run Bayesian optimisation -- fitting data is source from file upload. - Select the number of iterations. - The only optimised values are: Ms, and K and J. - All other parameters are treated as constants. - Narrow the bounds to improve the optimisation. + Run Bayesian optimisation -- fitting data is source from file upload. + Select the number of iterations. + The only optimised values are: Ms, and K and J. + All other parameters are treated as constants. + Narrow the bounds to improve the optimisation. - ### Fixed parameters - Parameters marked as fixed will not be optimised. - In the panel below, __the lower bound value will be fixed__ + #### Fixed parameters + Parameters marked as fixed will not be optimised. + In the panel below, __the lower bound value will be fixed__ - ### Values to set on the left hand side panel - - H axis - - K axis for each layer - - alpha for each layer - - thickness for each layer - """ - ) + #### Values to set on the left hand side panel + - H axis + - K axis for each layer + - alpha for each layer + - thickness for each layer + """ + ) st.number_input( "Number of iterations", min_value=1, @@ -211,6 +225,14 @@ grid[3].toggle(f"fix {param_name} ({i+1})", key=f"check_{param_name}{i}") with vsd_tab: + with st.expander("VSD simulation parameters"): + st.markdown( + """### Simulation info + + This is Voltage Spin Diode experiment. + The frequency is the sinusoidal excitation frequency fed into the system. + """ + ) st.number_input( "Frequency min (GHz)", min_value=0.0, max_value=50.0, value=0.0, key="fmin" ) @@ -233,32 +255,35 @@ options=["Rx", "Ry"], key="res_type", index=1, + help="Resistance type for the VSD simulation. PRL, Kim et al. 2016", ) st.selectbox("excitation axis", options=["x", "y", "z"], key="Hoex", index=1) - st.markdown( - """### Simulation info - - This is Voltage Spin Diode experiment (WIP). - """ - ) st.button("Simulate VSD", on_click=simulate_vsd, key="VSD_btn") with pimm_tab: fn = simulate_pimm - st.markdown( - """### Simulation info - - This app simulates PIMM experiment. - """ - ) + with st.expander("PIMM simulation parameters"): + st.write( + """### Simulation info\n + This app simulates PIMM experiment wherin we apply a step impulse. + Here, the impulse is field impulse (in the experiment that could be Oersted field impulse). + This impulse is called $\mathbf{H}_{oe}$. + """ + ) st.button("Simulate PIMM", on_click=simulate_pimm, key="PIMM_btn") st.number_input( - "Hoe (kA/m)", min_value=0.05, max_value=50.0, value=0.05, key="Hoe_mag" + "Hoe (kA/m)", + min_value=0.05, + max_value=50.0, + value=0.05, + key="Hoe_mag", + help="Magnitude of the Oersted field impulse (PIMM excitation)", ) st.radio( "Hoe axis", options=["x", "y", "z"], key="Hoeaxis", index=1, + help="Direction of the Oersted field impulse", ) diff --git a/view/utils.py b/view/utils.py index aa2596e..bcf4a80 100644 --- a/view/utils.py +++ b/view/utils.py @@ -1,3 +1,4 @@ +import os from typing import List import numpy as np @@ -14,6 +15,10 @@ "J": "uJ/m^2", } + +PREGENERATED_ANGLE_VALUES = [float(np.random.randint(0, 20)) for _ in range(30)] + + def extract_max_resonance_lines( spectrum, h_vals: List[float], @@ -43,3 +48,9 @@ def extract_max_resonance_lines( freqs = frequencies[indices].ravel() + [None] * n_diff linespectra[h_vals[i]] = freqs return linespectra + + +def get_init_kval(i: int = 0): + if i < len(PREGENERATED_ANGLE_VALUES): + return PREGENERATED_ANGLE_VALUES[i] + return 0.0 From a3b032472b4dbdc2e1ffaf8fa4c1ef7aa20c7617 Mon Sep 17 00:00:00 2001 From: Jakub Date: Sun, 15 Sep 2024 15:50:10 +0200 Subject: [PATCH 2/2] Update Streamlit app links --- README.md | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 8555abb..2a4c38b 100644 --- a/README.md +++ b/README.md @@ -6,7 +6,7 @@ [![pages-build-deployment](https://github.com/LemurPwned/cmtj/actions/workflows/pages/pages-build-deployment/badge.svg?branch=gh-pages)](https://github.com/LemurPwned/cmtj/actions/workflows/pages/pages-build-deployment) [![Version](https://img.shields.io/pypi/v/cmtj)](https://pypi.org/project/cmtj/) [![License](https://img.shields.io/pypi/l/cmtj.svg)](https://github.com/LemurPwned/cmtj/blob/master/LICENSE) -[![Streamlit](https://static.streamlit.io/badges/streamlit_badge_black_white.svg)](http://cmtj-simulations.streamlit.app/) +[![Streamlit](https://static.streamlit.io/badges/streamlit_badge_black_white.svg)](https://cmtj-app.streamlit.app/spectrum) ![Downloads](https://img.shields.io/pypi/dm/cmtj.svg) ## Table of contents @@ -36,7 +36,14 @@ It is also possible to connect devices in parallel or in series to have electric ## Web GUI -Check out the [streamlit hosted demo here](http://cmtj-simulations.streamlit.app/). You can simulate PIMM spectra and Spin-Diode spectra there. Let us know if you have any issues with the demo. +Check out the [streamlit hosted demo here](https://cmtj-app.streamlit.app/spectrum). +You can simulate: + +* PIMM spectra and Spin-Diode spectra +* Try some optimization fitting +* Fit multi-domain or multi-level M(H) or R(H) loops in [Domain mode](https://cmtj-app.streamlit.app) + +Let us know if you have any issues with the demo. ## Quickstart