diff --git a/Examples/Artemis_Tutorial_Cases/README.md b/Examples/Artemis_Tutorial_Cases/README.md new file mode 100644 index 000000000..5b875d2ea --- /dev/null +++ b/Examples/Artemis_Tutorial_Cases/README.md @@ -0,0 +1,18 @@ +# Test Cases +This file contains two test cases: an air-filled waveguide and a magnetically tunable filter. The simulation parameters are defined in an input file, "/inputs_*". To simulate the air-filled waveguide using NERSC's Perlmutter, run +> sbatch perl.run +The simulation will generate a folder of diagnostic data titled "/diags". Next, to convert the diagnostic data into readable txt files, run +> python3 batch_convert.py +The converted data will be stored in "/raw_data*". Finally, to extract the phase constant and attenuation spectrums, run +> python3 alphaExtraction.py + +This folders contains the following files: +- alphaExtraction.py : Extracts the components of the wavevector, k, from a "/raw_data*" file +- batch_convert.py : Converts diagnostic data, "/diags", into readable data, "raw_data*". +- inputs_3d_empty_X_band : Input file defining simulation parameters for air-filled waveguide case. +- inputs_3d_LLG_filter : Input file defining simulation parameters for the magnetically tunable filter. +- main3d.gnu.TPROF.MTMPI.CUDA.ex : Executable file used for air-filled waveguide; the non-LLG case. +- main3d.gnu.TPROF.MTMPI.CUDA.LLG.ex : Executable file used for the magnetically tunable filter; the LLG case. +- perl_empty.run : Run file for NERSC's Perlmutter for the non-LLG case. +- perl_llg.run : Run file for NERSC's Perlmutter for LLG case. + diff --git a/Examples/Artemis_Tutorial_Cases/alphaExtraction.py b/Examples/Artemis_Tutorial_Cases/alphaExtraction.py new file mode 100644 index 000000000..1f59f47a3 --- /dev/null +++ b/Examples/Artemis_Tutorial_Cases/alphaExtraction.py @@ -0,0 +1,184 @@ +import glob +import math +import re +import numpy as np +import matplotlib.pyplot as plt +from scipy.interpolate import interp1d +from scipy.linalg import hankel, svd + + +def esprit(x, p, n_sources): + """1D Total Least Squares ESPRIT (TLS-ESPRIT) algorithm.""" + N = len(x) + H = hankel(x[: N - p + 1], x[N - p :]) + + U, _, _ = svd(H, full_matrices=False) + Us = U[:, :n_sources] + + Us1, Us2 = Us[:-1, :], Us[1:, :] + + K = np.hstack((Us1, Us2)) + _, _, Vh_k = svd(K) + V = Vh_k.conj().T + + V12 = V[:n_sources, n_sources:] + V22 = V[n_sources:, n_sources:] + Psi = -V12 @ np.linalg.inv(V22) + + z = np.linalg.eigvals(Psi) + wm = np.log(z) + + return wm + + +# --- Parameters & Grid Setup --- +# Physical Constants & Waveguide Setup +EP0 = 8.854e-12 +MU0 = 4 * math.pi * 1e-7 +CC = 3e8 + +# GENERAL PARAMETERS +freq = 10.5e9 # frequency +XX, YY, ZZ = 14.95e-3, 10.16e-3, 500.0e-3 # dimensions +NX, NY, NZ = 128, 4, 128 # number of cells + +DX, DY, DZ = XX / NX, YY / NY, ZZ / NZ + +DT = 3.11275047e-13 # time step in seconds +TP = 1.0 / freq +Tstep_end = 1000000 # maximum time steps +Tstep_int = 100 # interval between time steps + +file_pattern = "./raw_data_empty/*.txt" # data location + +k0 = 1 +k1 = NZ - 1 +dk = 1 +NZ_record = int(math.floor(k1 / dk)) + +Tstep_record = math.floor((4.7 / 6.3 * Tstep_end) / Tstep_int) +Tstep_trunc = Tstep_record + +Time_axis = np.arange(Tstep_int, Tstep_end + 1, Tstep_int) + +# Spectral Setup +fs = 1.0 / (Tstep_int * DT) +Tstep_zeros = 0 +Tstep_ext = Tstep_record + Tstep_zeros +NFFT_ext = 1 << (Tstep_ext - 1).bit_length() + +f_ext = (fs / 2.0) * np.linspace(0, 1, (NFFT_ext // 2) + 1) +f_ext_GHz = f_ext / 1e9 +f_ext_step = np.mean(np.diff(f_ext)) + +# --- Data Import --- +txt_files = sorted( + glob.glob(file_pattern), + key=lambda x: int(re.search(r"(\d+)\.txt$", x).group(1)) if re.search(r"(\d+)\.txt$", x) else 0, +) + +Ey_obs_ret = np.zeros((NZ_record, Tstep_record)) +for tcount, fname in enumerate(txt_files[:Tstep_record]): + data = np.loadtxt(fname, skiprows=0) + Ey_obs_ret[:, tcount] = data[:NZ_record, 1] + +ZZ_axis = DZ * np.arange(k0, k1 + 1, dk) + +# --- FFT Processing --- +ey_ext_rec = np.hstack( + ( + Ey_obs_ret[:, :Tstep_trunc], + np.zeros((NZ_record, Tstep_record - Tstep_trunc)), + np.zeros((NZ_record, Tstep_zeros)), + ) +) + +ey_extff_rec = np.fft.fft(ey_ext_rec, n=NFFT_ext, axis=1) / Tstep_ext + +# --- ESPRIT Processing (Integrated into Main Script) --- +ko1_m = int(math.floor((ZZ - DZ) / (DZ * dk))) +ko2_m = int(math.floor(DZ / (DZ * dk))) + +num_freqs = int(NFFT_ext // 2) + 1 +beta_esprit = np.zeros(num_freqs) +alpha_esprit = np.zeros(num_freqs) +spatial_step = DZ * dk + +for fpoint in range(num_freqs): + slice_m = ey_extff_rec[ko2_m:ko1_m, fpoint] + wm = esprit(slice_m, p=3, n_sources=1) + + beta_esprit[fpoint] = np.abs(np.imag(wm[0])) / spatial_step + alpha_esprit[fpoint] = np.abs(np.real(wm[0])) / spatial_step + +N_m = len(beta_esprit) + +# --- Data Exports --- +alpha_m_db = 10 * np.log10(np.exp(2)) * (alpha_esprit / 100) + +alpha_data = np.column_stack((f_ext_GHz[:N_m], alpha_m_db)) +np.savetxt("alpha.txt", alpha_data, header="Frequency Forward", fmt="%.6f", comments="") + +beta_data = np.column_stack((f_ext_GHz[:N_m], beta_esprit)) +np.savetxt("beta.txt", beta_data, header="Frequency Forward", fmt="%.6f", comments="") + +Nf_ey = int(round(10.2e9 / f_ext_step)) +ey_data = np.column_stack((ZZ_axis, np.abs(ey_extff_rec[:, Nf_ey]))) +np.savetxt("Ey.txt", ey_data, header="Z position |Ey|", fmt="%.6f", comments="") + +# --- Target Plots --- + +# 1. beta.png +plt.figure() +plt.plot(f_ext_GHz[:N_m], beta_esprit, label="Forward") +plt.xlim([8, 12]) +plt.xlabel("Frequency [GHz]") +plt.ylabel(r"Phase constant $\beta$ [rad/m]") +plt.title("Dispersion relation") +plt.savefig("fig1_dispersion.png", dpi=300) +plt.close() + +# 2. alpha.png +plt.figure() +plt.plot(f_ext_GHz[:N_m], alpha_esprit, label="Forward") +plt.xlim([8.5, 12]) +plt.xlabel("Frequency [GHz]") +plt.ylabel(r"Attenuation rate $\alpha$ [Np/m]") +plt.title("Attenuation Spectrum") +plt.savefig("fig2_attenuation.png", dpi=300) +plt.close() + +# 3. ey.png +plt.figure() +plt.plot(ZZ_axis * 1000, np.abs(ey_extff_rec[:, Nf_ey])) +plt.ylabel("|Ey| [V/m]") +plt.xlabel("z position [mm]") +plt.title(f"{f_ext_GHz[Nf_ey]:.2f} GHz") +plt.savefig("fig3_ey.png", dpi=300) +plt.close() + +# --- Interpolation and Export --- +f_start, f_end = 8.5e9, 12e9 +Nf_start = int(round(f_start / f_ext_step)) +Nf_end = int(round(f_end / f_ext_step)) + +f_interp = np.linspace(f_ext_GHz[Nf_start], f_ext_GHz[Nf_end], 150) +f_source = f_ext_GHz[Nf_start : Nf_end + 1] +alpha_source = alpha_m_db[Nf_start : Nf_end + 1] + +spline_interp = interp1d(f_source, alpha_source, kind="cubic", fill_value="extrapolate") +alpha_total = 4.9 * spline_interp(f_interp) + +# 4. alpha_total.png +plt.figure() +plt.plot(f_interp, alpha_total, "k") +plt.xlim([8.5, 12]) +plt.ylim([0, 90]) +plt.xlabel("Frequency [GHz]") +plt.ylabel("Attenuation over 49mm [dB]") +plt.title("Attenuation Spectrum") +plt.savefig("fig4_atten_total.png", dpi=300) +plt.close() + +alpha_tot_data = np.column_stack((f_interp, alpha_total)) +np.savetxt("alpha_total.txt", alpha_tot_data, header="Frequency Attenuation", fmt="%.6f", comments="") diff --git a/Examples/Artemis_Tutorial_Cases/batch_convert.py b/Examples/Artemis_Tutorial_Cases/batch_convert.py new file mode 100644 index 000000000..a682f80d9 --- /dev/null +++ b/Examples/Artemis_Tutorial_Cases/batch_convert.py @@ -0,0 +1,43 @@ +import glob +import os +import re +import numpy as np +import yt + +# Directory containing your plotfile folders (e.g., plt00001000, plt00002000, etc.) +input_dir = "./diags_LLG" +output_dir = "./raw_data_LLG" + +os.makedirs(output_dir, exist_ok=True) + +# Find and sort all plotfile folders numerically by timestep +plt_folders = sorted( + [d for d in glob.glob(os.path.join(input_dir, "plt*")) if os.path.isdir(d)], + key=lambda x: int(re.search(r'plt(\d+)', x).group(1)) +) + +print(f"Found {len(plt_folders)} plotfiles to process.") + +for plt_path in plt_folders: + # Extract integer timestep from folder name (e.g., 'plt00001000' -> 1000) + step_num = int(re.search(r'plt(\d+)', plt_path).group(1)) + + # Load dataset and extract covering grid + ds = yt.load(plt_path) + ad0 = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions) + + # Extract 3D Ey array and slice 1D profile at transverse cell (0, 0) + E_array = ad0['Ey'].to_ndarray() + ey_1d = E_array[0, 0, :] + + # Generate 0-based k indices + k_indices = np.arange(len(ey_1d)) + + # Format 2-column data matrix: [k_index, Ey] + output_data = np.column_stack((k_indices, ey_1d)) + + # Save output text file matching your naming scheme (e.g., 1000.txt) + output_filename = os.path.join(output_dir, f"{step_num}.txt") + np.savetxt(output_filename, output_data, fmt="%d %12.8e") + +print(f"Batch conversion complete! Output saved to '{output_dir}'.") diff --git a/Examples/Waveguide/inputs_3d_LLG_filter b/Examples/Artemis_Tutorial_Cases/inputs_3d_LLG_filter similarity index 95% rename from Examples/Waveguide/inputs_3d_LLG_filter rename to Examples/Artemis_Tutorial_Cases/inputs_3d_LLG_filter index 19215a0d7..e5684150a 100644 --- a/Examples/Waveguide/inputs_3d_LLG_filter +++ b/Examples/Artemis_Tutorial_Cases/inputs_3d_LLG_filter @@ -101,10 +101,10 @@ warpx.My_external_grid_function(x,y,z)= "1.3926e5 * (x<=thickness-width/2) + 0 * warpx.Mz_external_grid_function(x,y,z) = "0." #Diagnostics -diagnostics.diags_names = zline_BernadoFilter0111_ +diagnostics.diags_names = plt -zline_BernadoFilter0111_.intervals = 2 -zline_BernadoFilter0111_.diag_lo = 0.0 0.0 -250.e-3 -zline_BernadoFilter0111_.diag_hi = 0.0 0.0 250.e-3 -zline_BernadoFilter0111_.diag_type = Full -zline_BernadoFilter0111_.fields_to_plot = Ey Hx Hz +plt.intervals = 20 +plt.diag_lo = 0.0 0.0 -250.e-3 +plt.diag_hi = 0.0 0.0 250.e-3 +plt.diag_type = Full +plt.fields_to_plot = Ey Hx Hz diff --git a/Examples/Waveguide/inputs_3d_empty_X_band b/Examples/Artemis_Tutorial_Cases/inputs_3d_empty_X_band similarity index 72% rename from Examples/Waveguide/inputs_3d_empty_X_band rename to Examples/Artemis_Tutorial_Cases/inputs_3d_empty_X_band index 5d4ecec41..047c7bddd 100644 --- a/Examples/Waveguide/inputs_3d_empty_X_band +++ b/Examples/Artemis_Tutorial_Cases/inputs_3d_empty_X_band @@ -8,18 +8,18 @@ ####### GENERAL PARAMETERS ###### ################################# max_step = 1000000 -amr.n_cell = 128 4 128 # number of cells spanning the domain in each coordinate direction at level 0 +amr.n_cell = Nx Ny Nz # number of cells spanning the domain in each coordinate direction at level 0 amr.max_grid_size = 64 # maximum size of each AMReX box, used to decompose the domain amr.blocking_factor = 4 # only meaningful for AMR geometry.dims = 3 -boundary.field_lo = pec pec pml # PEC on side walls; PML at -z end -boundary.field_hi = pec pec pec # PEC on side walls; PEC at +z end +boundary.field_lo = pec pec pml # PEC on side walls; PML at -z end +boundary.field_hi = pec pec pec # PEC on side walls; PEC at +z end # waveguide width 14.95mm, height is 11.43mm, and length is 500mm # 14.95mm to make the wave guide work on fundamental mode at 10.5GHz # 11.43mm < 14.95mm so that the fundamental mode is TE10 -geometry.prob_lo = -7.475e-3 -5.715e-3 -250.0e-3 # must be consistent with my_constants.length and .diag_hi/lo -geometry.prob_hi = 7.475e-3 5.715e-3 250.0e-3 +geometry.prob_lo = -width/2 -height/2 -length/2 # must be consistent with my_constants.length and .diag_hi/lo +geometry.prob_hi = width/2 height/2 length/2 amr.max_level = 0 @@ -28,14 +28,23 @@ my_constants.c = 299792458. my_constants.width = 14.95e-3 # waveguide width is 14.95mm my_constants.height = 11.43e-3 # waveguide height is 10.16mm my_constants.length = 500.0e-3 # waveguide length is 400mm -my_constants.rjx = 1.0e-4 # the x dimension of current source cross-section -my_constants.rjz = 10.0e-4 # the z dimension of current source cross-section; should be just larger than 2*dz my_constants.wavelength = 0.0286 # frequency is 10.5 GHz my_constants.TP = 9.5238e-11 # Gaussian pulse width, 1 x time period of excitation my_constants.flag_none = 0 # no source flag my_constants.flag_hs = 1 # hard source flag my_constants.flag_ss = 2 # soft source flag +my_constants.Nx = 128 +my_constants.Ny = 4 +my_constants.Nz = 128 + +my_constants.dx = width/Nx +my_constants.dy = height/Ny +my_constants.dz = length/Nz + +my_constants.rjx = dx/100 +my_constants.rjz = dz/100 + ################################# ############ NUMERICS ########### ################################# @@ -58,19 +67,27 @@ macroscopic.mu_function(x,y,z) = "1.25663706212e-06" ################################# warpx.B_excitation_on_grid_style = "parse_B_excitation_grid_function" -warpx.Bx_excitation_grid_function(x,y,z,t) = "1.25663706212e-06 * 2.5e-5 * (exp(-(t-3*TP)**2/(2*TP**2))*cos(2*pi*c/wavelength*t)) * cos(x/(width/2)*(pi/2)) * (z > - rjz/2 + length/2)" +warpx.Bx_excitation_grid_function(x,y,z,t) = "1.25663706212e-06 * 2.5e-5 * (exp(-(t-3*TP)**2/(2*TP**2))*cos(2*pi*c/wavelength*t)) * cos(x/(width/2)*(pi/2))" # current line source Jy = Bx at the +z end of waveguide; spanning over entire y dimension warpx.By_excitation_grid_function(x,y,z,t) = "0.0" warpx.Bz_excitation_grid_function(x,y,z,t) = "0.0" -warpx.Bx_excitation_flag_function(x,y,z) = "flag_ss * (z > - rjz/2 + length/2) + flag_none * (z <= - rjz/2 + length/2)" +warpx.Bx_excitation_flag_function(x,y,z) = "flag_ss * (z > - rjz + length/2 - dz) + flag_none * (z <= - rjz + length/2 - dz)" warpx.By_excitation_flag_function(x,y,z) = "flag_none" warpx.Bz_excitation_flag_function(x,y,z) = "flag_none" #Diagnostics -diagnostics.diags_names = zline_BernadoFilter0126_ +diagnostics.diags_names = plt chk + +plt.intervals = 100 +plt.diag_lo = 0.0 0.0 -length/2 +plt.diag_hi = 0.0 0.0 length/2 +plt.diag_type = Full +plt.fields_to_plot = Ey Bx + +# checkpoint +chk.intervals = 250000 +chk.diag_type = Full +chk.format = checkpoint +chk.file_min_digits = 7 -zline_BernadoFilter0126_.intervals = 40 -zline_BernadoFilter0126_.diag_lo = 0.0 0.0 -250.e-3 -zline_BernadoFilter0126_.diag_hi = 0.0 0.0 250.e-3 -zline_BernadoFilter0126_.diag_type = Full -zline_BernadoFilter0126_.fields_to_plot = Ey Bx +# amr.restart = diags/chk0500000 diff --git a/Examples/Artemis_Tutorial_Cases/perl_empty.run b/Examples/Artemis_Tutorial_Cases/perl_empty.run new file mode 100644 index 000000000..112792d2a --- /dev/null +++ b/Examples/Artemis_Tutorial_Cases/perl_empty.run @@ -0,0 +1,81 @@ +#!/bin/bash +#SBATCH -N 1 +#SBATCH -C gpu +#SBATCH -G 4 +#SBATCH -q regular +#SBATCH -J waveguide_test +#SBATCH -A m3766 +#SBATCH -t 02:00:00 + +# OpenMP settings: +export OMP_NUM_THREADS=1 +export OMP_PLACES=threads +export OMP_PROC_BIND=spread + +export MPICH_GPU_SUPPORT_ENABLED=1 +export AMREX_CUDA_SHARE_PARAM_MEM=1 + +#run the application: +# applications may perform better with --gpu-bind=none instead of --gpu-bind=single:1 +srun -n 4 -c 32 --cpu_bind=cores -G 4 ./main3d.gnu.TPROF.MTMPI.CUDA.ex inputs_3d_empty_X_band + +#!/bin/bash + +rm -r raw_data +mkdir raw_data +rm -r FieldInTime +mkdir FieldInTime + +for timestep in $(seq 0 500 1700000); do + newnumber='000000000'${timestep} # get number, pack with zeros + newnumber=${newnumber:(-9)} # the last eight characters + + ./WritePlotfileToASCII3d.gnu.x86-milan.ex infile=diags/CPW_surf${newnumber} | tee raw_data/raw_${newnumber}.txt + + for i in $(seq 53 1 53); do # Hz + for j in $(seq 125 1 125); do + string="^$i $j 501 5" # Hz + prefix="^$i " + if [[ "$string" =~ ^"$prefix" ]]; then + ALLNAME=${i}_${j}_${newnumber}_Hz + grep "$string" raw_data/raw_${newnumber}.txt | tee -a raw_data/$ALLNAME.txt # + + XYNAME=${i}_${j}_Hz + grep "$string" raw_data/$ALLNAME.txt | tee -a FieldInTime/$XYNAME.txt # begin y + fi + + done + done + + for i in $(seq 60 1 60); do # Ez + for j in $(seq 95 1 95); do + string="^$i $j 505 2" # Ez + prefix="^$i " + if [[ "$string" =~ ^"$prefix" ]]; then + ALLNAME=${i}_${j}_${newnumber}_Ez + grep "$string" raw_data/raw_${newnumber}.txt | tee -a raw_data/$ALLNAME.txt # + + XYNAME=${i}_${j}_Ez + grep "$string" raw_data/$ALLNAME.txt | tee -a FieldInTime/$XYNAME.txt # begin y + fi + + done + done + + for i in $(seq 60 1 60); do # Hx + for j in $(seq 100 1 100); do + string="^$i $j 508 3" # Hx + prefix="^$i " + if [[ "$string" =~ ^"$prefix" ]]; then + ALLNAME=${i}_${j}_${newnumber}_Hx + grep "$string" raw_data/raw_${newnumber}.txt | tee -a raw_data/$ALLNAME.txt # + + XYNAME=${i}_${j}_Hx + grep "$string" raw_data/$ALLNAME.txt | tee -a FieldInTime/$XYNAME.txt # begin y + fi + + done + done + +done + diff --git a/Examples/Artemis_Tutorial_Cases/perl_llg.run b/Examples/Artemis_Tutorial_Cases/perl_llg.run new file mode 100644 index 000000000..65d24d45e --- /dev/null +++ b/Examples/Artemis_Tutorial_Cases/perl_llg.run @@ -0,0 +1,20 @@ +#!/bin/bash +#SBATCH -N 2 +#SBATCH -C gpu +#SBATCH -G 8 +#SBATCH -q regular +#SBATCH -J wg_LLG_test +#SBATCH -A m3766 +#SBATCH -t 03:00:00 + +# OpenMP settings: +export OMP_NUM_THREADS=1 +export OMP_PLACES=threads +export OMP_PROC_BIND=spread + +export MPICH_GPU_SUPPORT_ENABLED=1 +export AMREX_CUDA_SHARE_PARAM_MEM=1 + +#run the application: +# applications may perform better with --gpu-bind=none instead of --gpu-bind=single:1 +srun -n 4 -c 32 --cpu_bind=cores -G 8 ./main3d.gnu.TPROF.MTMPI.CUDA.LLG.ex inputs_3d_LLG_filter \ No newline at end of file diff --git a/Examples/Tests/afm_testing/inputs_NiO b/Examples/Tests/afm_testing/inputs_NiO new file mode 100644 index 000000000..7181bc000 --- /dev/null +++ b/Examples/Tests/afm_testing/inputs_NiO @@ -0,0 +1,130 @@ +#################################################################################################### +## This input file simulates thin-film NiO +## The plane wave excitation is the time-dependent modified Gaussian pulse + +## L = (M1 - M2)/2 +## M = (M1 + M2)/2 +#################################################################################################### + +################################ +####### GENERAL PARAMETERS ###### +################################# +max_step = 5 +amr.n_cell = Nx Ny Nz # number of cells spanning the domain in each coordinate direction at level 0 +amr.max_grid_size = 16 # maximum size of each AMReX box, used to decompose the domain +amr.blocking_factor = 4 # only meaningful for AMR +geometry.dims = 3 +boundary.field_lo = periodic periodic pml +boundary.field_hi = periodic periodic pml + +geometry.prob_lo = -length/2 -width/2 -height/2 # must be consistent with my_constants.length and .diag_hi/lo +geometry.prob_hi = length/2 width/2 height/2 + +amr.max_level = 0 + +my_constants.pi = 3.14159265359 +my_constants.c = 299792458. +my_constants.length = 10.0e-7 +my_constants.width = 10.0e-7 +my_constants.height = 45.0e-9 +my_constants.wavelength = 0.0003 # frequency is 1 THz +my_constants.TP = 1e-12 # Gaussian pulse width, 1 x time period of excitation +my_constants.flag_none = 0 # no source flag +my_constants.flag_hs = 1 # hard source flag +my_constants.flag_ss = 2 # soft source flag + +my_constants.Nx = 4 +my_constants.Ny = 4 +my_constants.Nz = 32 + +my_constants.dx = length/Nx +my_constants.dy = width/Ny +my_constants.dz = height/Nz + +my_constants.rjx = dx/100 +my_constants.rjz = dz/100 + +my_constants.thickness = height + +################################# +############ NUMERICS ########### +################################# +warpx.verbose = 1 +warpx.use_filter = 0 +warpx.cfl = 0.8 +warpx.mag_time_scheme_order = 2 # default 1 +warpx.mag_M_normalization = 1 # 1 is saturated +warpx.mag_LLG_coupling = 1 +warpx.mag_LLG_exchange_coupling = 1 +warpx.mag_LLG_afm_coupling = 1 + +algo.em_solver_medium = macroscopic # vacuum/macroscopic + +algo.macroscopic_sigma_method = laxwendroff # laxwendroff or backwardeuler + +macroscopic.sigma_function(x,y,z) = "0.0" + +macroscopic.epsilon_function(x,y,z) = "8.8541878128e-12 " # air-filled + +macroscopic.mu_function(x,y,z) = "1.25663706212e-06" + +#unit conversion: 1 Gauss = (1000/4pi) A/m +macroscopic.mag_Ms_init_style = "parse_mag_Ms_function" # parse or "constant" +macroscopic.mag_Ms_function(x,y,z) = "3.82e5" # in unit A/m + +macroscopic.mag_alpha_init_style = "parse_mag_alpha_function" # parse or "constant" +macroscopic.mag_alpha_function(x,y,z) = "5e-5" # alpha is unitless + +macroscopic.mag_gamma_init_style = "parse_mag_gamma_function" # parse or "constant" +macroscopic.mag_gamma_function(x,y,z) = "-1.759e11" # gyromagnetic ratio is constant for electrons in all materials + +macroscopic.mag_exchange_init_style = "parse_mag_exchange_function" # parse or "constant" +macroscopic.mag_exchange_function(x,y,z) = "1.15e-13" # + +macroscopic.mag_afm_init_style = "parse_mag_afm_function" # parse or "constant" +macroscopic.mag_afm_function(x,y,z) = "-5.539e7" # J/m + +macroscopic.mag_max_iter = 100 # maximum number of M iteration in each time step +macroscopic.mag_tol = 1.e-7 # M magnitude relative error tolerance compared to previous iteration +macroscopic.mag_normalized_error = 0.1 # if M magnitude relatively changes more than this value, raise a red flag + +################################# +############ FIELDS ############# +################################# + +warpx.H_excitation_on_grid_style = "parse_H_excitation_grid_function" +warpx.Hx_excitation_grid_function(x,y,z,t) = "0.0" # plane source +warpx.Hy_excitation_grid_function(x,y,z,t) = "0.0" +warpx.Hz_excitation_grid_function(x,y,z,t) = "0.0" +warpx.Hx_excitation_flag_function(x,y,z) = "flag_none" # plane source +warpx.Hy_excitation_flag_function(x,y,z) = "flag_none" +warpx.Hz_excitation_flag_function(x,y,z) = "flag_none" + +#unit conversion: 1 Gauss = 1 Oersted = (1000/4pi) A/m +#calculation of H_bias: H_bias (oe) = frequency / 2.8e6 + +warpx.H_bias_excitation_on_grid_style = parse_h_bias_excitation_grid_function +#warpx.Hx_bias_excitation_grid_function(x,y,z,t) = "0.707 * 1e5 * (exp(-(t-3*TP)**2/(2*TP**2))*cos(2*pi*c/wavelength*t)) * cos(x/(width/2)*(pi/2))" +#warpx.Hy_bias_excitation_grid_function(x,y,z,t) = "0.0" +#warpx.Hz_bias_excitation_grid_function(x,y,z,t) = "0.707 * 1e5 * (exp(-(t-3*TP)**2/(2*TP**2))*cos(2*pi*c/wavelength*t)) * cos(x/(width/2)*(pi/2))" +warpx.Hx_bias_excitation_grid_function(x,y,z,t) = "1.0" +warpx.Hy_bias_excitation_grid_function(x,y,z,t) = "0.0" +warpx.Hz_bias_excitation_grid_function(x,y,z,t) = "0.0" +warpx.Hx_bias_excitation_flag_function(x,y,z) = "flag_hs" +warpx.Hy_bias_excitation_flag_function(x,y,z) = "flag_none" +warpx.Hz_bias_excitation_flag_function(x,y,z) = "flag_none" + +warpx.M_ext_grid_init_style = parse_M_ext_grid_function +warpx.Mx_external_grid_function(x,y,z)= "3.82e5 * 0.408" +warpx.My_external_grid_function(x,y,z)= "3.82e5 * 0.408" # in unit A/m, equal to 1750 Gauss; Ms=0 triggers off LLG +warpx.Mz_external_grid_function(x,y,z) = "-3.82e5 * 0.817" + +#Diagnostics +diagnostics.diags_names = plt + +plt.intervals = 1 +#plt.diag_lo = 0.0 -height/2 0.0 +#plt.diag_hi = 0.0 height/2 0.0 +plt.diag_type = Full +plt.fields_to_plot = Lx_xface Lx_yface Lx_zface Ly_xface Ly_yface Ly_zface Lz_xface Lz_yface Lz_zface Mx_xface My_xface Mz_xface +#plt.fields_to_plot = My_xface My_yface diff --git a/Examples/Waveguide/inputs_3d_LLG_filter_restarttest b/Examples/Waveguide/inputs_3d_LLG_filter_restarttest deleted file mode 100644 index 444eb444d..000000000 --- a/Examples/Waveguide/inputs_3d_LLG_filter_restarttest +++ /dev/null @@ -1,113 +0,0 @@ -#################################################################################################### -## This input file simulates the thin-film ferromagnetic inserted waveguide -## PEC applied on x, y and +z boundaries -## PML applied on -z boundary -## The plane wave excitation is the time-dependent modified Gaussian pulse -## This input file requires USE_LLG=TRUE in the GNUMakefile. -#################################################################################################### - -################################ -####### GENERAL PARAMETERS ###### -################################# -max_step = 100 -amr.n_cell = 512 4 256 # number of cells spanning the domain in each coordinate direction at level 0 -amr.max_grid_size = 128 # maximum size of each AMReX box, used to decompose the domain -amr.blocking_factor = 4 # only meaningful for AMR -geometry.dims = 3 -boundary.field_lo = pec pec pml # PEC on side walls; PML at -z end -boundary.field_hi = pec pec pec # PEC on side walls; PEC at +z end - -# waveguide width 14.95mm, height is 11.43mm, and length is 280mm -# 14.95mm to make the wave guide work on fundamental mode at 10.5GHz -# 11.43mm < 14.95mm so that the fundamental mode is TE10 -geometry.prob_lo = -7.475e-3 -5.715e-3 -250.0e-3 # must be consistent with my_constants.length and .diag_hi/lo -geometry.prob_hi = 7.475e-3 5.715e-3 250.0e-3 - -amr.max_level = 0 - -my_constants.pi = 3.14159265359 -my_constants.c = 299792458. -my_constants.thickness = 0.45e-3 # thicness of the film is 0.45mm -my_constants.width = 14.95e-3 # waveguide width is 14.95mm -my_constants.height = 11.43e-3 # waveguide height is 10.16mm -my_constants.length = 500.0e-3 # waveguide length is 400mm -my_constants.rjx = 1.0e-4 # the x dimension of current source cross-section -my_constants.rjz = 10.0e-4 # the z dimension of current source cross-section; should be just larger than 2*dz -my_constants.wavelength = 0.0286 # frequency is 10.5 GHz -my_constants.TP = 9.5238e-11 # Gaussian pulse width, 1 x time period of excitation -my_constants.flag_none = 0 # no source flag -my_constants.flag_hs = 1 # hard source flag -my_constants.flag_ss = 2 # soft source flag -my_constants.epr = 13 # relative permittivity of ferrite slab - -################################# -############ NUMERICS ########### -################################# -warpx.verbose = 1 -warpx.use_filter = 0 -warpx.cfl = 0.8 -warpx.mag_time_scheme_order = 2 # default 1 -warpx.mag_M_normalization = 1 # 1 is saturated -warpx.mag_LLG_coupling = 1 - -algo.em_solver_medium = macroscopic # vacuum/macroscopic - -algo.macroscopic_sigma_method = laxwendroff # laxwendroff or backwardeuler - -macroscopic.sigma_function(x,y,z) = "0.0" - -macroscopic.epsilon_function(x,y,z) = "epr * 8.8541878128e-12 * (x<=thickness-width/2) + 8.8541878128e-12 * (x>thickness-width/2)" # EPr is 13 of the ferrite slab - -macroscopic.mu_function(x,y,z) = "1.25663706212e-06" # MUr is not predefined in ferrite materials - -#unit conversion: 1 Gauss = (1000/4pi) A/m -macroscopic.mag_Ms_init_style = "parse_mag_Ms_function" # parse or "constant" -macroscopic.mag_Ms_function(x,y,z) = "1.3926e5 * (x<=thickness-width/2) + 0 * (x>thickness-width/2)" # in unit A/m, equal to 1750 Gauss; Ms=0 triggers off LLG - -macroscopic.mag_alpha_init_style = "parse_mag_alpha_function" # parse or "constant" -macroscopic.mag_alpha_function(x,y,z) = "0.0051 * (x<=thickness-width/2) + 0 * (x>thickness-width/2)" # alpha is unitless, calculated from linewidth Delta_H = 35 Oersted - -macroscopic.mag_gamma_init_style = "parse_mag_gamma_function" # parse or "constant" -macroscopic.mag_gamma_function(x,y,z) = "-1.759e11 * (x<=thickness-width/2) + 0 * (x>thickness-width/2)" # gyromagnetic ratio is constant for electrons in all materials - -macroscopic.mag_max_iter = 100 # maximum number of M iteration in each time step -macroscopic.mag_tol = 1.e-7 # M magnitude relative error tolerance compared to previous iteration -macroscopic.mag_normalized_error = 0.1 # if M magnitude relatively changes more than this value, raise a red flag - -################################# -############ FIELDS ############# -################################# - -warpx.H_excitation_on_grid_style = "parse_H_excitation_grid_function" -warpx.Hx_excitation_grid_function(x,y,z,t) = "2.5e-5 * (exp(-(t-3*TP)**2/(2*TP**2))*cos(2*pi*c/wavelength*t)) * cos(x/(width/2)*(pi/2)) * (z > - rjz/2 + length/2)" # plane source -# magnetic current line source at the +z end of waveguide; spanning over entire y dimension -warpx.Hy_excitation_grid_function(x,y,z,t) = "0.0" -warpx.Hz_excitation_grid_function(x,y,z,t) = "0.0" -warpx.Hx_excitation_flag_function(x,y,z) = "flag_ss * (z > - rjz/2 + length/2) + flag_none * (z <= - rjz/2 + length/2)" # plane source -warpx.Hy_excitation_flag_function(x,y,z) = "flag_none" -warpx.Hz_excitation_flag_function(x,y,z) = "flag_none" - -#unit conversion: 1 Gauss = 1 Oersted = (1000/4pi) A/m -#calculation of H_bias: H_bias (oe) = frequency / 2.8e6 - -warpx.H_bias_ext_grid_init_style = parse_H_bias_ext_grid_function -warpx.Hx_bias_external_grid_function(x,y,z)= "0." -warpx.Hy_bias_external_grid_function(x,y,z)= "2.3475e+05 * (x<=thickness-width/2) + 0 * (x>thickness-width/2)" # in A/m, equal to 2950 Oersted -warpx.Hz_bias_external_grid_function(x,y,z)= "0." - -warpx.M_ext_grid_init_style = parse_M_ext_grid_function -warpx.Mx_external_grid_function(x,y,z)= "0." -warpx.My_external_grid_function(x,y,z)= "1.3926e5 * (x<=thickness-width/2) + 0 * (x>thickness-width/2)" # in unit A/m, equal to 1750 Gauss; Ms=0 triggers off LLG -warpx.Mz_external_grid_function(x,y,z) = "0." - -# Diagnostics -diagnostics.diags_names = plt chk -plt.intervals = 10 -plt.diag_type = Full -plt.fields_to_plot = Ex Ey Ez Hx Hy Hz Bx By Bz Mx_xface My_xface Mz_xface Mx_yface My_yface Mz_yface Mx_zface My_zface Mz_zface - -chk.format = checkpoint -chk.intervals = 10 -chk.diag_type = Full - -#amr.restart = diags/chk000090 diff --git a/Examples/Waveguide/inputs_regression_3d_LLG_filter b/Examples/Waveguide/inputs_regression_3d_LLG_filter deleted file mode 100644 index 3c4457872..000000000 --- a/Examples/Waveguide/inputs_regression_3d_LLG_filter +++ /dev/null @@ -1,107 +0,0 @@ -#################################################################################################### -## This input file simulates the thin-film ferromagnetic inserted waveguide -## PEC applied on x, y and +z boundaries -## PML applied on -z boundary -## The plane wave excitation is the time-dependent modified Gaussian pulse -## This input file requires USE_LLG=TRUE in the GNUMakefile. -#################################################################################################### - -################################ -####### GENERAL PARAMETERS ###### -################################# -max_step = 3 -amr.n_cell = 1024 4 512 # number of cells spanning the domain in each coordinate direction at level 0 -amr.max_grid_size = 1024 # maximum size of each AMReX box, used to decompose the domain -amr.blocking_factor = 4 # only meaningful for AMR -geometry.dims = 3 -boundary.field_lo = pec pec pml # PEC on side walls; PML at -z end -boundary.field_hi = pec pec pec # PEC on side walls; PEC at +z end - -# waveguide width 14.95mm, height is 11.43mm, and length is 280mm -# 14.95mm to make the wave guide work on fundamental mode at 10.5GHz -# 11.43mm < 14.95mm so that the fundamental mode is TE10 -geometry.prob_lo = -7.475e-3 -5.715e-3 -250.0e-3 # must be consistent with my_constants.length and .diag_hi/lo -geometry.prob_hi = 7.475e-3 5.715e-3 250.0e-3 - -amr.max_level = 0 - -my_constants.pi = 3.14159265359 -my_constants.c = 299792458. -my_constants.thickness = 0.45e-3 # thicness of the film is 0.45mm -my_constants.width = 14.95e-3 # waveguide width is 14.95mm -my_constants.height = 11.43e-3 # waveguide height is 10.16mm -my_constants.length = 500.0e-3 # waveguide length is 400mm -my_constants.rjx = 1.0e-4 # the x dimension of current source cross-section -my_constants.rjz = 10.0e-4 # the z dimension of current source cross-section; should be just larger than 2*dz -my_constants.wavelength = 0.0286 # frequency is 10.5 GHz -my_constants.TP = 9.5238e-11 # Gaussian pulse width, 1 x time period of excitation -my_constants.flag_none = 0 # no source flag -my_constants.flag_hs = 1 # hard source flag -my_constants.flag_ss = 2 # soft source flag -my_constants.epr = 13 # relative permittivity of ferrite slab - -################################# -############ NUMERICS ########### -################################# -warpx.verbose = 1 -warpx.use_filter = 0 -warpx.cfl = 0.8 -warpx.mag_time_scheme_order = 2 # default 1 -warpx.mag_M_normalization = 1 # 1 is saturated -warpx.mag_LLG_coupling = 1 - -algo.em_solver_medium = macroscopic # vacuum/macroscopic - -algo.macroscopic_sigma_method = laxwendroff # laxwendroff or backwardeuler - -macroscopic.sigma_function(x,y,z) = "0.0" - -macroscopic.epsilon_function(x,y,z) = "epr * 8.8541878128e-12 * (x<=thickness-width/2) + 8.8541878128e-12 * (x>thickness-width/2)" # EPr is 13 of the ferrite slab - -macroscopic.mu_function(x,y,z) = "1.25663706212e-06" # MUr is not predefined in ferrite materials - -#unit conversion: 1 Gauss = (1000/4pi) A/m -macroscopic.mag_Ms_init_style = "parse_mag_Ms_function" # parse or "constant" -macroscopic.mag_Ms_function(x,y,z) = "1.3926e5 * (x<=thickness-width/2) + 0 * (x>thickness-width/2)" # in unit A/m, equal to 1750 Gauss; Ms=0 triggers off LLG - -macroscopic.mag_alpha_init_style = "parse_mag_alpha_function" # parse or "constant" -macroscopic.mag_alpha_function(x,y,z) = "0.0051 * (x<=thickness-width/2) + 0 * (x>thickness-width/2)" # alpha is unitless, calculated from linewidth Delta_H = 35 Oersted - -macroscopic.mag_gamma_init_style = "parse_mag_gamma_function" # parse or "constant" -macroscopic.mag_gamma_function(x,y,z) = "-1.759e11 * (x<=thickness-width/2) + 0 * (x>thickness-width/2)" # gyromagnetic ratio is constant for electrons in all materials - -macroscopic.mag_max_iter = 100 # maximum number of M iteration in each time step -macroscopic.mag_tol = 1.e-7 # M magnitude relative error tolerance compared to previous iteration -macroscopic.mag_normalized_error = 0.1 # if M magnitude relatively changes more than this value, raise a red flag - -################################# -############ FIELDS ############# -################################# - -warpx.H_excitation_on_grid_style = "parse_H_excitation_grid_function" -warpx.Hx_excitation_grid_function(x,y,z,t) = "2.5e-5 * (exp(-(t-3*TP)**2/(2*TP**2))*cos(2*pi*c/wavelength*t)) * cos(x/(width/2)*(pi/2)) * (z > - rjz/2 + length/2)" # plane source -# magnetic current line source at the +z end of waveguide; spanning over entire y dimension -warpx.Hy_excitation_grid_function(x,y,z,t) = "0.0" -warpx.Hz_excitation_grid_function(x,y,z,t) = "0.0" -warpx.Hx_excitation_flag_function(x,y,z) = "flag_ss * (z > - rjz/2 + length/2) + flag_none * (z <= - rjz/2 + length/2)" # plane source -warpx.Hy_excitation_flag_function(x,y,z) = "flag_none" -warpx.Hz_excitation_flag_function(x,y,z) = "flag_none" - -#unit conversion: 1 Gauss = 1 Oersted = (1000/4pi) A/m -#calculation of H_bias: H_bias (oe) = frequency / 2.8e6 - -warpx.H_bias_ext_grid_init_style = parse_H_bias_ext_grid_function -warpx.Hx_bias_external_grid_function(x,y,z)= "0." -warpx.Hy_bias_external_grid_function(x,y,z)= "2.3475e+05 * (x<=thickness-width/2) + 0 * (x>thickness-width/2)" # in A/m, equal to 2950 Oersted -warpx.Hz_bias_external_grid_function(x,y,z)= "0." - -warpx.M_ext_grid_init_style = parse_M_ext_grid_function -warpx.Mx_external_grid_function(x,y,z)= "0." -warpx.My_external_grid_function(x,y,z)= "1.3926e5 * (x<=thickness-width/2) + 0 * (x>thickness-width/2)" # in unit A/m, equal to 1750 Gauss; Ms=0 triggers off LLG -warpx.Mz_external_grid_function(x,y,z) = "0." - -# Diagnostics -diagnostics.diags_names = plt -plt.intervals = 3 -plt.diag_type = Full -plt.fields_to_plot = Ex Ey Ez Hx Hy Hz Bx By Bz Mx_xface My_xface Mz_xface Mx_yface My_yface Mz_yface Mx_zface My_zface Mz_zface diff --git a/Source/Diagnostics/FullDiagnostics.cpp b/Source/Diagnostics/FullDiagnostics.cpp index f9814c84a..d499bfdcc 100644 --- a/Source/Diagnostics/FullDiagnostics.cpp +++ b/Source/Diagnostics/FullDiagnostics.cpp @@ -127,6 +127,7 @@ FullDiagnostics::BackwardCompatibility () void FullDiagnostics::Flush ( int i_buffer ) { + // This function should be moved to Diagnostics when plotfiles/openpmd format // is supported for BackTransformed Diagnostics, in BTDiagnostics class. auto & warpx = WarpX::GetInstance(); @@ -591,7 +592,57 @@ FullDiagnostics::InitializeFieldFunctors (int lev) m_all_field_functors[lev].resize(ntot); // Fill vector of functors for all components except individual cylindrical modes. + #ifdef WARPX_MAG_LLG + // AFM diagnostics: M1, M2, L (3 components stored on each face, like M) + // and H1, H2 (one value per direction, like H). + // Returns nullptr if the name is not an AFM variable. + const std::array xyz{"x","y","z"}; + auto make_afm_functor = [&](const std::string& name) -> std::unique_ptr + { + for (int idir = 0; idir < 3; ++idir) { + // H1x, H1y, H1z, H2x, ... (same pattern as Hx, Hy, Hz) + if (name == "H1" + xyz[idir]) { + return std::make_unique( + warpx.get_pointer_H1field_fp(lev, idir), lev, m_crse_ratio); + } + if (name == "H2" + xyz[idir]) { + return std::make_unique( + warpx.get_pointer_H2field_fp(lev, idir), lev, m_crse_ratio); + } + // M1x_xface, M1y_xface, ..., M2z_zface, Lx_xface, ... Lz_zface + // idir selects the face (MultiFab), icomp selects the vector component. + for (int icomp = 0; icomp < 3; ++icomp) { + const std::string suffix = xyz[icomp] + "_" + xyz[idir] + "face"; + if (name == "M1" + suffix) { + amrex::Print() << "functor build M1 dir " << idir << " ptr=" + << static_cast(warpx.get_pointer_M1field_fp(lev, idir)) << "\n"; + return std::make_unique( + warpx.get_pointer_M1field_fp(lev, idir), lev, m_crse_ratio, true, 1, icomp); + } + if (name == "M2" + suffix) { + return std::make_unique( + warpx.get_pointer_M2field_fp(lev, idir), lev, m_crse_ratio, true, 1, icomp); + } + if (name == "L" + suffix) { + return std::make_unique( + warpx.get_pointer_Lfield_fp(lev, idir), lev, m_crse_ratio, true, 1, icomp); + } + + } + } + return nullptr; + }; +#endif for (int comp=0; comp(warpx.get_pointer_Efield_aux(lev, 2), lev, m_crse_ratio); } else if ( m_varnames[comp] == "Bz" ){ @@ -797,6 +848,7 @@ FullDiagnostics::PrepareFieldDataForOutput () m_output_species[i_buffer][i].m_diag_domain = m_geom_output[i_buffer][0].ProbDomain(); } } + } void diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 71d387c2d..a103d91e1 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -499,7 +499,12 @@ WarpX::OneStep_nosub (Real cur_time) if (mag_time_scheme_order==1){ MacroscopicEvolveHM(0.5*dt[0]); // we now have M^{n+1/2} and H^{n+1/2} } else if (mag_time_scheme_order==2){ - MacroscopicEvolveHM_2nd(0.5*dt[0]); // we now have M^{n+1/2} and H^{n+1/2} + if (mag_LLG_afm_coupling==1){ //if afm is required + MacroscopicEvolveHM_AFM_2nd(0.5*dt[0]); + } + else{ + MacroscopicEvolveHM_2nd(0.5*dt[0]); // we now have M^{n+1/2} and H^{n+1/2} + } } else { amrex::Abort("unsupported mag_time_scheme_order for M field"); } @@ -559,7 +564,12 @@ WarpX::OneStep_nosub (Real cur_time) if (mag_time_scheme_order==1){ MacroscopicEvolveHM(0.5*dt[0]); // we now have M^{n+1} and H^{n+1} } else if (mag_time_scheme_order==2){ - MacroscopicEvolveHM_2nd(0.5*dt[0]); // we now have M^{n+1} and H^{n+1} + if (mag_LLG_afm_coupling==1){ //if afm is required + MacroscopicEvolveHM_AFM_2nd(0.5*dt[0]); + } + else{ + MacroscopicEvolveHM_2nd(0.5*dt[0]); // we now have M^{n+1/2} and H^{n+1/2} + } } else { amrex::Abort("unsupported mag_time_scheme_order for M field"); } diff --git a/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt b/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt index c95e646ba..2199f6b75 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt +++ b/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt @@ -20,6 +20,7 @@ if(WarpX_MAG_LLG) target_sources(WarpX PRIVATE MacroscopicEvolveHM_2nd.cpp + MacroscopicEvolveHM_AFM_2nd.cpp MacroscopicEvolveHM.cpp EvolveHPML.cpp ) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index c1d24d74c..84eed0cee 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -163,6 +163,21 @@ class FiniteDifferenceSolver amrex::Real const dt, std::unique_ptr const ¯oscopic_properties); + void MacroscopicEvolveHM_AFM_2nd( + int lev, + std::array, 3> &M1field, + std::array, 3> &M2field, + std::array, 3> &Mfield, // Mfield contains three components MultiFab + std::array, 3> &Hfield, + std::array, 3> &H1field, + std::array, 3> &H2field, + std::array, 3> &Bfield, + std::array, 3> &Lfield, + std::array, 3> const &H_biasfield, // H bias + std::array, 3> const &Efield, + amrex::Real const dt, + std::unique_ptr const ¯oscopic_properties); + #endif #endif // ifndef WARPX_DIM_RZ @@ -356,6 +371,22 @@ class FiniteDifferenceSolver std::array, 3> const &Efield, amrex::Real const dt, std::unique_ptr const ¯oscopic_properties); + + template< typename T_Algo > + void MacroscopicEvolveHMCartesian_AFM_2nd( + int lev, + std::array, 3> &M1field, + std::array, 3> &M2field, + std::array, 3> &Mfield, // Mfield contains three components MultiFab + std::array, 3> &Hfield, + std::array, 3> &H1field, + std::array, 3> &H2field, + std::array, 3> &Bfield, + std::array, 3> &Lfield, + std::array, 3> const &H_biasfield, // H bias + std::array, 3> const &Efield, + amrex::Real const dt, + std::unique_ptr const ¯oscopic_properties); #endif template< typename T_Algo > diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHM_2nd.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHM_2nd.cpp index fd9bef015..3bedba24f 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHM_2nd.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHM_2nd.cpp @@ -15,7 +15,7 @@ blank #include "Utils/WarpXConst.H" #include #include "Utils/WarpXUtil.H" -#include + using namespace amrex; diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHM_AFM_2nd.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHM_AFM_2nd.cpp new file mode 100644 index 000000000..ac232f33b --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHM_AFM_2nd.cpp @@ -0,0 +1,1647 @@ +/* copyright +blank +*/ + +#include "WarpX.H" +#include "Utils/WarpXAlgorithmSelection.H" +#include "FiniteDifferenceSolver.H" +#ifdef WARPX_DIM_RZ +#include "FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H" +#else +#include "FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H" +#endif +#include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H" + +#include "Utils/WarpXConst.H" +#include +#include "Utils/WarpXUtil.H" +#include + +using namespace amrex; + +/** + * \brief Update H and M fields with iterative correction, over one timestep + */ + +#ifndef WARPX_DIM_RZ +#ifdef WARPX_MAG_LLG + +void FiniteDifferenceSolver::MacroscopicEvolveHM_AFM_2nd( + // The MField here is a vector of three multifabs, with M on each face, and each multifab is a three-component multifab. + // Each M-multifab has three components, one for each component in x, y, z. (All multifabs are four dimensional, (i,j,k,n)), where, n=1 for E, B, but, n=3 for M_xface, M_yface, M_zface + int lev, + std::array, 3> &M1field, + std::array, 3> &M2field, + std::array, 3> &Mfield, // Mfield contains three components MultiFab + std::array, 3> &Hfield, + std::array, 3> &H1field, + std::array, 3> &H2field, + std::array, 3> &Bfield, + std::array, 3> &Lfield, + std::array, 3> const &H_biasfield, // H bias + std::array, 3> const &Efield, + amrex::Real const dt, + std::unique_ptr const ¯oscopic_properties) { + + if (m_fdtd_algo == ElectromagneticSolverAlgo::Yee){ + MacroscopicEvolveHMCartesian_AFM_2nd(lev, M1field, M2field, Mfield, Hfield, H1field, H2field, Bfield, Lfield, H_biasfield, Efield, dt, macroscopic_properties); + } else { + amrex::Abort("Only yee algorithm is compatible for M updates."); + } +} // closes function MacroscopicEvolveHM_AMF_2nd +#endif +#ifdef WARPX_MAG_LLG +template +void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian_AFM_2nd( + int lev, + std::array, 3> &M1field, + std::array, 3> &M2field, + std::array, 3> &Mfield, + std::array, 3> &Hfield, + std::array, 3> &H1field, + std::array, 3> &H2field, + std::array, 3> &Bfield, + std::array, 3> &Lfield, + std::array, 3> const &H_biasfield, // H bias + std::array, 3> const &Efield, + amrex::Real const dt, + std::unique_ptr const ¯oscopic_properties) { + + // obtain the maximum relative amount we let M deviate from Ms before aborting + amrex::Real mag_normalized_error = macroscopic_properties->getmag_normalized_error(); + + auto &warpx = WarpX::GetInstance(); + int coupling = warpx.mag_LLG_coupling; + int M_normalization = warpx.mag_M_normalization; + int mag_exchange_coupling = warpx.mag_LLG_exchange_coupling; + int mag_anisotropy_coupling = warpx.mag_LLG_anisotropy_coupling; + int mag_afm_coupling = warpx.mag_LLG_afm_coupling; + + // build temporary vector Mfield_prev, Mfield_error, a_temp, a_temp_static, b_temp_static + std::array, 3> Hfield_old; // H^(old_time) before the current time step + std::array, 3> H1field_old; // H1^(old_time) before the current time step + std::array, 3> H2field_old; // H2^(old_time) before the current time step + std::array, 3> Mfield_old; // M^(old_time) before the current time step + std::array, 3> M1field_old; // M1^(old_time) before the current time step + std::array, 3> M2field_old; // M2^(old_time) before the current time step + std::array, 3> Mfield_prev; // M^(new_time) of the (r-1)th iteration + std::array, 3> M1field_prev; // M1^(new_time) of the (r-1)th iteration + std::array, 3> M2field_prev; // M2^(new_time) of the (r-1)th iteration + std::array, 3> Mfield_error; // The error of the M field between the two consecutive iterations + std::array, 3> M1field_error; // The error of the M1 field between the two consecutive iterations + std::array, 3> M2field_error; // The error of the M2 field between the two consecutive iterations + std::array, 3> a1_temp; // right-hand side of vector a, see the documentation + std::array, 3> a2_temp; // right-hand side of vector a, see the documentation + std::array, 3> a1_temp_static; // α M^(old_time)/|M| in the right-hand side of vector a, see the documentation + std::array, 3> a2_temp_static; // α M^(old_time)/|M| in the right-hand side of vector a, see the documentation + std::array, 3> b1_temp_static; // right-hand side of vector b, see the documentation + std::array, 3> b2_temp_static; // right-hand side of vector b, see the documentation + + amrex::GpuArray const& mu_stag = macroscopic_properties->mu_IndexType; + amrex::GpuArray const& Bx_stag = macroscopic_properties->Bx_IndexType; + amrex::GpuArray const& By_stag = macroscopic_properties->By_IndexType; + amrex::GpuArray const& Bz_stag = macroscopic_properties->Bz_IndexType; + amrex::GpuArray const& Hx_stag = macroscopic_properties->Hx_IndexType; + amrex::GpuArray const& Hy_stag = macroscopic_properties->Hy_IndexType; + amrex::GpuArray const& Hz_stag = macroscopic_properties->Hz_IndexType; + amrex::GpuArray const& macro_cr = macroscopic_properties->macro_cr_ratio; + amrex::GpuArray const& anisotropy_axis = macroscopic_properties->mag_LLG_anisotropy_axis; + + // Initialize Hfield_old (H^(old_time)), Mfield_old (M^(old_time)), Mfield_prev (M^[(new_time),r-1]), Mfield_error + for (int i = 0; i < 3; i++){ + Hfield_old[i].reset(new MultiFab(Hfield[i]->boxArray(), Hfield[i]->DistributionMap(), 1, Hfield[i]->nGrow())); + H1field_old[i].reset(new MultiFab(H1field[i]->boxArray(), H1field[i]->DistributionMap(), 1, H1field[i]->nGrow())); + H2field_old[i].reset(new MultiFab(H2field[i]->boxArray(), H2field[i]->DistributionMap(), 1, H2field[i]->nGrow())); + Mfield_old[i].reset(new MultiFab(Mfield[i]->boxArray(), Mfield[i]->DistributionMap(), 3, Mfield[i]->nGrow())); + M1field_old[i].reset(new MultiFab(M1field[i]->boxArray(), M1field[i]->DistributionMap(), 3, M1field[i]->nGrow())); + M2field_old[i].reset(new MultiFab(M2field[i]->boxArray(), M2field[i]->DistributionMap(), 3, M2field[i]->nGrow())); + Mfield_prev[i].reset(new MultiFab(Mfield[i]->boxArray(), Mfield[i]->DistributionMap(), 3, Mfield[i]->nGrow())); + M1field_prev[i].reset(new MultiFab(M1field[i]->boxArray(), M1field[i]->DistributionMap(), 3, M1field[i]->nGrow())); + M2field_prev[i].reset(new MultiFab(M2field[i]->boxArray(), M2field[i]->DistributionMap(), 3, M2field[i]->nGrow())); + Mfield_error[i].reset(new MultiFab(Mfield[i]->boxArray(), Mfield[i]->DistributionMap(), 3, Mfield[i]->nGrow())); + M1field_error[i].reset(new MultiFab(M1field[i]->boxArray(), M1field[i]->DistributionMap(), 3, M1field[i]->nGrow())); + M2field_error[i].reset(new MultiFab(M2field[i]->boxArray(), M2field[i]->DistributionMap(), 3, M2field[i]->nGrow())); + Mfield_error[i]->setVal(0.); // reset Mfield_error to zero + M1field_error[i]->setVal(0.); // reset M1field_error to zero + M2field_error[i]->setVal(0.); // reset M2field_error to zero + MultiFab::Copy(*Hfield_old[i], *Hfield[i], 0, 0, 1, Hfield[i]->nGrow()); + MultiFab::Copy(*H1field_old[i], *H1field[i], 0, 0, 1, H1field[i]->nGrow()); + MultiFab::Copy(*H2field_old[i], *H2field[i], 0, 0, 1, H2field[i]->nGrow()); + MultiFab::Copy(*Mfield_old[i], *Mfield[i], 0, 0, 3, Mfield[i]->nGrow()); + MultiFab::Copy(*M1field_old[i], *M1field[i], 0, 0, 3, M1field[i]->nGrow()); + MultiFab::Copy(*M2field_old[i], *M2field[i], 0, 0, 3, M2field[i]->nGrow()); + MultiFab::Copy(*Mfield_prev[i], *Mfield[i], 0, 0, 3, Mfield[i]->nGrow()); + MultiFab::Copy(*M1field_prev[i], *M1field[i], 0, 0, 3, M1field[i]->nGrow()); + MultiFab::Copy(*M2field_prev[i], *M2field[i], 0, 0, 3, M2field[i]->nGrow()); + } + // initialize a_temp, a_temp_static, b_temp_static + for (int i = 0; i < 3; i++){ + a1_temp[i].reset(new MultiFab(M1field[i]->boxArray(), M1field[i]->DistributionMap(), 3, M1field[i]->nGrow())); + a2_temp[i].reset(new MultiFab(M2field[i]->boxArray(), M2field[i]->DistributionMap(), 3, M2field[i]->nGrow())); + a1_temp_static[i].reset(new MultiFab(M1field[i]->boxArray(), M1field[i]->DistributionMap(), 3, M1field[i]->nGrow())); + a2_temp_static[i].reset(new MultiFab(M2field[i]->boxArray(), M2field[i]->DistributionMap(), 3, M2field[i]->nGrow())); + b1_temp_static[i].reset(new MultiFab(M1field[i]->boxArray(), M1field[i]->DistributionMap(), 3, M1field[i]->nGrow())); + b2_temp_static[i].reset(new MultiFab(M2field[i]->boxArray(), M2field[i]->DistributionMap(), 3, M2field[i]->nGrow())); + } + + amrex::MultiFab& mu_mf = macroscopic_properties->getmu_mf(); + + // calculate the b_temp_static, a_temp_static + for (MFIter mfi(*a1_temp_static[0], TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + auto& mag_Ms_xface_mf = macroscopic_properties->getmag_Ms_mf(0); + auto& mag_Ms_yface_mf = macroscopic_properties->getmag_Ms_mf(1); + auto& mag_Ms_zface_mf = macroscopic_properties->getmag_Ms_mf(2); + auto& mag_alpha_xface_mf = macroscopic_properties->getmag_alpha_mf(0); + auto& mag_alpha_yface_mf = macroscopic_properties->getmag_alpha_mf(1); + auto& mag_alpha_zface_mf = macroscopic_properties->getmag_alpha_mf(2); + auto& mag_gamma_xface_mf = macroscopic_properties->getmag_gamma_mf(0); + auto& mag_gamma_yface_mf = macroscopic_properties->getmag_gamma_mf(1); + auto& mag_gamma_zface_mf = macroscopic_properties->getmag_gamma_mf(2); + auto& mag_exchange_xface_mf = macroscopic_properties->getmag_exchange_mf(0); + auto& mag_exchange_yface_mf = macroscopic_properties->getmag_exchange_mf(1); + auto& mag_exchange_zface_mf = macroscopic_properties->getmag_exchange_mf(2); + auto& mag_anisotropy_xface_mf = macroscopic_properties->getmag_anisotropy_mf(0); + auto& mag_anisotropy_yface_mf = macroscopic_properties->getmag_anisotropy_mf(1); + auto& mag_anisotropy_zface_mf = macroscopic_properties->getmag_anisotropy_mf(2); + auto& mag_afm_xface_mf = macroscopic_properties->getmag_afm_mf(0); + auto& mag_afm_yface_mf = macroscopic_properties->getmag_afm_mf(1); + auto& mag_afm_zface_mf = macroscopic_properties->getmag_afm_mf(2); + + // extract material properties + Array4 const& mag_Ms_xface_arr = mag_Ms_xface_mf.array(mfi); + Array4 const& mag_Ms_yface_arr = mag_Ms_yface_mf.array(mfi); + Array4 const& mag_Ms_zface_arr = mag_Ms_zface_mf.array(mfi); + Array4 const& mag_alpha_xface_arr = mag_alpha_xface_mf.array(mfi); + Array4 const& mag_alpha_yface_arr = mag_alpha_yface_mf.array(mfi); + Array4 const& mag_alpha_zface_arr = mag_alpha_zface_mf.array(mfi); + Array4 const& mag_gamma_xface_arr = mag_gamma_xface_mf.array(mfi); + Array4 const& mag_gamma_yface_arr = mag_gamma_yface_mf.array(mfi); + Array4 const& mag_gamma_zface_arr = mag_gamma_zface_mf.array(mfi); + Array4 const& mag_exchange_xface_arr = mag_exchange_xface_mf.array(mfi); + Array4 const& mag_exchange_yface_arr = mag_exchange_yface_mf.array(mfi); + Array4 const& mag_exchange_zface_arr = mag_exchange_zface_mf.array(mfi); + Array4 const& mag_anisotropy_xface_arr = mag_anisotropy_xface_mf.array(mfi); + Array4 const& mag_anisotropy_yface_arr = mag_anisotropy_yface_mf.array(mfi); + Array4 const& mag_anisotropy_zface_arr = mag_anisotropy_zface_mf.array(mfi); + Array4 const& mag_afm_xface_arr = mag_afm_xface_mf.array(mfi); + Array4 const& mag_afm_yface_arr = mag_afm_yface_mf.array(mfi); + Array4 const& mag_afm_zface_arr = mag_afm_zface_mf.array(mfi); + + // extract field data + Array4 const &M_xface = Mfield[0]->array(mfi); // note M_xface include x,y,z components at |_x faces + Array4 const &M_yface = Mfield[1]->array(mfi); // note M_yface include x,y,z components at |_y faces + Array4 const &M_zface = Mfield[2]->array(mfi); // note M_zface include x,y,z components at |_z faces + Array4 const &M1_xface = M1field[0]->array(mfi); // note M_xface include x,y,z components at |_x faces + Array4 const &M1_yface = M1field[1]->array(mfi); // note M_yface include x,y,z components at |_y faces + Array4 const &M1_zface = M1field[2]->array(mfi); // note M_zface include x,y,z components at |_z faces + Array4 const &M2_xface = M2field[0]->array(mfi); // note M_xface include x,y,z components at |_x faces + Array4 const &M2_yface = M2field[1]->array(mfi); // note M_yface include x,y,z components at |_y faces + Array4 const &M2_zface = M2field[2]->array(mfi); // note M_zface include x,y,z components at |_z faces + Array4 const &Hx_bias = H_biasfield[0]->array(mfi); // Hx_bias is the x component at |_x faces + Array4 const &Hy_bias = H_biasfield[1]->array(mfi); // Hy_bias is the y component at |_y faces + Array4 const &Hz_bias = H_biasfield[2]->array(mfi); // Hz_bias is the z component at |_z faces + Array4 const &Hx_old = Hfield_old[0]->array(mfi); // Hx_old is the x component at |_x faces + Array4 const &Hy_old = Hfield_old[1]->array(mfi); // Hy_old is the y component at |_y faces + Array4 const &Hz_old = Hfield_old[2]->array(mfi); // Hz_old is the z component at |_z faces + Array4 const &Hx1_old = H1field_old[0]->array(mfi); // Hx1_old is the x component at |_x faces + Array4 const &Hy1_old = H1field_old[1]->array(mfi); // Hy1_old is the y component at |_y faces + Array4 const &Hz1_old = H1field_old[2]->array(mfi); // Hz1_old is the z component at |_z faces + Array4 const &Hx2_old = H2field_old[0]->array(mfi); // Hx2_old is the x component at |_x faces + Array4 const &Hy2_old = H2field_old[1]->array(mfi); // Hy2_old is the y component at |_y faces + Array4 const &Hz2_old = H2field_old[2]->array(mfi); // Hz2_old is the z component at |_z faces + + // extract field data of a_temp_static and b_temp_static + Array4 const &a1_temp_static_xface = a1_temp_static[0]->array(mfi); + Array4 const &a1_temp_static_yface = a1_temp_static[1]->array(mfi); + Array4 const &a1_temp_static_zface = a1_temp_static[2]->array(mfi); + Array4 const &b1_temp_static_xface = b1_temp_static[0]->array(mfi); + Array4 const &b1_temp_static_yface = b1_temp_static[1]->array(mfi); + Array4 const &b1_temp_static_zface = b1_temp_static[2]->array(mfi); + Array4 const &a2_temp_static_xface = a2_temp_static[0]->array(mfi); + Array4 const &a2_temp_static_yface = a2_temp_static[1]->array(mfi); + Array4 const &a2_temp_static_zface = a2_temp_static[2]->array(mfi); + Array4 const &b2_temp_static_xface = b2_temp_static[0]->array(mfi); + Array4 const &b2_temp_static_yface = b2_temp_static[1]->array(mfi); + Array4 const &b2_temp_static_zface = b2_temp_static[2]->array(mfi); + + // extract tileboxes for which to loop + amrex::IntVect Mxface_stag = Mfield[0]->ixType().toIntVect(); + amrex::IntVect Myface_stag = Mfield[1]->ixType().toIntVect(); + amrex::IntVect Mzface_stag = Mfield[2]->ixType().toIntVect(); + Box const &tbx = mfi.tilebox(Mxface_stag); /* just define which grid type */ + Box const &tby = mfi.tilebox(Myface_stag); + Box const &tbz = mfi.tilebox(Mzface_stag); + + // Extract stencil coefficients for calculating the exchange field H_exchange and the anisotropy field H_anisotropy + amrex::Real const *const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr(); + int const n_coefs_x = m_stencil_coefs_x.size(); + amrex::Real const *const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr(); + int const n_coefs_y = m_stencil_coefs_y.size(); + amrex::Real const *const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr(); + int const n_coefs_z = m_stencil_coefs_z.size(); + + // loop over cells and update fields + amrex::ParallelFor(tbx, + [=] AMREX_GPU_DEVICE(int i, int j, int k) { + + // determine if the material is nonmagnetic or not + if (mag_Ms_xface_arr(i,j,k) > 0._rt){ + + // when working on M_xface(i,j,k, 0:2) we have direct access to M_xface(i,j,k,0:2) and Hx(i,j,k) + // Hy and Hz can be acquired by interpolation + + // H_bias + amrex::Real Hx1_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mxface_stag, Mxface_stag, Hx_bias); + amrex::Real Hy1_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Myface_stag, Mxface_stag, Hy_bias); + amrex::Real Hz1_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mzface_stag, Mxface_stag, Hz_bias); + + amrex::Real Hx2_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mxface_stag, Mxface_stag, Hx_bias); + amrex::Real Hy2_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Myface_stag, Mxface_stag, Hy_bias); + amrex::Real Hz2_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mzface_stag, Mxface_stag, Hz_bias); + + if (coupling == 1){ + // H_eff = H_maxwell + H_bias + H_exchange + H_anisotropy + H_afm + + // H_maxwell - use H^(old_time) + Hx1_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mxface_stag, Mxface_stag, Hx_old); + Hy1_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Myface_stag, Mxface_stag, Hy_old); + Hz1_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mzface_stag, Mxface_stag, Hz_old); + + Hx2_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mxface_stag, Mxface_stag, Hx_old); + Hy2_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Myface_stag, Mxface_stag, Hy_old); + Hz2_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mzface_stag, Mxface_stag, Hz_old); + } + + if (mag_exchange_coupling == 1){ + + if (mag_exchange_xface_arr(i,j,k) == 0._rt) amrex::Abort("The mag_exchange_xface_arr(i,j,k) is 0.0 while including the exchange coupling term H_exchange for H_eff"); + + // H_exchange - use M^(old_time) + amrex::Real const H_exchange_coeff = 2.0 * mag_exchange_xface_arr(i,j,k) / PhysConst::mu0 / mag_Ms_xface_arr(i,j,k) / mag_Ms_xface_arr(i,j,k); + + amrex::Real Ms_lo_x = mag_Ms_xface_arr(i-1, j, k); + amrex::Real Ms_hi_x = mag_Ms_xface_arr(i+1, j, k); + amrex::Real Ms_lo_y = mag_Ms_xface_arr(i, j-1, k); + amrex::Real Ms_hi_y = mag_Ms_xface_arr(i, j+1, k); + amrex::Real Ms_lo_z = mag_Ms_xface_arr(i, j, k-1); + amrex::Real Ms_hi_z = mag_Ms_xface_arr(i, j, k+1); + + Hx1_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M1_xface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 0, 0); //Last argument is nodality -- xface = 0 + Hy1_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M1_xface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 1, 0); //Last argument is nodality -- xface = 0 + Hz1_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M1_xface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 2, 0); //Last argument is nodality -- xface = 0 + + Hx2_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M2_xface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 0, 0); //Last argument is nodality -- xface = 0 + Hy2_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M2_xface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 1, 0); //Last argument is nodality -- xface = 0 + Hz2_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M2_xface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 2, 0); //Last argument is nodality -- xface = 0 + } + + // always turn off anisotrophy for now + //mag_anisotropy_coupling = 0; + if (mag_anisotropy_coupling == 1){ + + if (mag_anisotropy_xface_arr(i,j,k) == 0._rt) amrex::Abort("The mag_anisotropy_xface_arr(i,j,k) is 0.0 while including the anisotropy coupling term H_anisotropy for H_eff"); + + // H_anisotropy - use M^(old_time) + amrex::Real M_dot_anisotropy_axis = 0.0; + for (int comp=0; comp<3; ++comp) { + M_dot_anisotropy_axis += M_xface(i, j, k, comp) * anisotropy_axis[comp]; + } + amrex::Real const H_anisotropy_coeff = - 2.0 * mag_anisotropy_xface_arr(i,j,k) / PhysConst::mu0 / mag_Ms_xface_arr(i,j,k) / mag_Ms_xface_arr(i,j,k); + Hx1_eff += H_anisotropy_coeff * M_dot_anisotropy_axis * anisotropy_axis[0]; + Hy1_eff += H_anisotropy_coeff * M_dot_anisotropy_axis * anisotropy_axis[1]; + Hz1_eff += H_anisotropy_coeff * M_dot_anisotropy_axis * anisotropy_axis[2]; + } + + if(mag_afm_coupling == 1){ + if (mag_afm_xface_arr(i,j,k) == 0._rt) amrex::Abort("The mag_afm_xface_arr(i,j,k) is 0.0 while including the afm coupling term H_afm for H_eff"); + + // H_afm - use M^(old_time) + amrex::Real const H_afm_coeff = 4.0 * mag_afm_xface_arr(i,j,k) / PhysConst::mu0 / mag_Ms_xface_arr(i,j,k) / mag_Ms_xface_arr(i,j,k); + Hx1_eff += H_afm_coeff * M2_xface(i,j,k,0); + Hy1_eff += H_afm_coeff * M2_xface(i,j,k,1); + Hz1_eff += H_afm_coeff * M2_xface(i,j,k,2); + + Hx2_eff += H_afm_coeff * M1_xface(i,j,k,0); + Hy2_eff += H_afm_coeff * M1_xface(i,j,k,1); + Hz2_eff += H_afm_coeff * M1_xface(i,j,k,2); + } + + // 0 = unsaturated; compute |M| locally. 1 = saturated; use M_s + amrex::Real M1_magnitude = (M_normalization == 0) ? std::sqrt(std::pow(M1_xface(i, j, k, 0), 2._rt) + std::pow(M1_xface(i, j, k, 1), 2._rt) + std::pow(M1_xface(i, j, k, 2), 2._rt)) + : mag_Ms_xface_arr(i,j,k); + amrex::Real M2_magnitude = (M_normalization == 0) ? std::sqrt(std::pow(M2_xface(i, j, k, 0), 2._rt) + std::pow(M2_xface(i, j, k, 1), 2._rt) + std::pow(M2_xface(i, j, k, 2), 2._rt)) + : mag_Ms_xface_arr(i,j,k); + // a_temp_static_coeff does not change in the current step for SATURATED materials; but it does change for UNSATURATED ones + amrex::Real a1_temp_static_coeff = mag_alpha_xface_arr(i,j,k) / M1_magnitude; + amrex::Real a2_temp_static_coeff = mag_alpha_xface_arr(i,j,k) / M2_magnitude; + + // calculate the b_temp_static_coeff (it is divided by 2.0 because the derivation is based on an interger dt, + // while in real simulations, the input dt is actually dt/2.0) + amrex::Real b_temp_static_coeff = - PhysConst::mu0 * amrex::Math::abs(mag_gamma_xface_arr(i,j,k)) / 2._rt; //not dependent on sublattice + + for (int comp=0; comp<3; ++comp) { + // calculate a_temp_static_xface + // all components on x-faces of grid + a1_temp_static_xface(i, j, k, comp) = a1_temp_static_coeff * M1_xface(i, j, k, comp); + a2_temp_static_xface(i, j, k, comp) = a2_temp_static_coeff * M2_xface(i, j, k, comp); + } + + // calculate b_temp_static_xface + // x component on x-faces of grid + b1_temp_static_xface(i, j, k, 0) = M1_xface(i, j, k, 0) + dt * b_temp_static_coeff * (M1_xface(i, j, k, 1) * Hz1_eff - M1_xface(i, j, k, 2) * Hy1_eff); + + // y component on x-faces of grid + b1_temp_static_xface(i, j, k, 1) = M1_xface(i, j, k, 1) + dt * b_temp_static_coeff * (M1_xface(i, j, k, 2) * Hx1_eff - M1_xface(i, j, k, 0) * Hz1_eff); + + // z component on x-faces of grid + b1_temp_static_xface(i, j, k, 2) = M1_xface(i, j, k, 2) + dt * b_temp_static_coeff * (M1_xface(i, j, k, 0) * Hy1_eff - M1_xface(i, j, k, 1) * Hx1_eff); + + //x component on x-faces of grid + b2_temp_static_xface(i, j, k, 0) = M2_xface(i, j, k, 0) + dt * b_temp_static_coeff * (M2_xface(i, j, k, 1) * Hz2_eff - M2_xface(i, j, k, 2) * Hy2_eff); + + // y component on x-faces of grid + b2_temp_static_xface(i, j, k, 1) = M2_xface(i, j, k, 1) + dt * b_temp_static_coeff * (M2_xface(i, j, k, 2) * Hx2_eff - M2_xface(i, j, k, 0) * Hz2_eff); + + // z component on x-faces of grid + b2_temp_static_xface(i, j, k, 2) = M2_xface(i, j, k, 2) + dt * b_temp_static_coeff * (M2_xface(i, j, k, 0) * Hy2_eff - M2_xface(i, j, k, 1) * Hx2_eff); + } + }); + + amrex::ParallelFor(tby, + [=] AMREX_GPU_DEVICE(int i, int j, int k) { + + // determine if the material is nonmagnetic or not + if (mag_Ms_yface_arr(i,j,k) > 0._rt){ + + // when working on M_yface(i,j,k,0:2) we have direct access to M_yface(i,j,k,0:2) and Hy(i,j,k) + // Hy and Hz can be acquired by interpolation + + // H_bias + amrex::Real Hx1_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mxface_stag, Myface_stag, Hx_bias); + amrex::Real Hy1_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Myface_stag, Myface_stag, Hy_bias); + amrex::Real Hz1_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mzface_stag, Myface_stag, Hz_bias); + + amrex::Real Hx2_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mxface_stag, Myface_stag, Hx_bias); + amrex::Real Hy2_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Myface_stag, Myface_stag, Hy_bias); + amrex::Real Hz2_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mzface_stag, Myface_stag, Hz_bias); + + if (coupling == 1){ + // H_eff = H_maxwell + H_bias + H_exchange + H_anisotropy + + // H_maxwell - use H^(old_time) + Hx1_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mxface_stag, Myface_stag, Hx_old); + Hy1_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Myface_stag, Myface_stag, Hy_old); + Hz1_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mzface_stag, Myface_stag, Hz_old); + + Hx2_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mxface_stag, Myface_stag, Hx_old); + Hy2_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Myface_stag, Myface_stag, Hy_old); + Hz2_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mzface_stag, Myface_stag, Hz_old); + } + + if (mag_exchange_coupling == 1){ + + if (mag_exchange_yface_arr(i,j,k) == 0._rt) amrex::Abort("The mag_exchange_yface_arr(i,j,k) is 0.0 while including the exchange coupling term H_exchange for H_eff"); + + // H_exchange - use M^(old_time) + amrex::Real const H_exchange_coeff = 2.0 * mag_exchange_yface_arr(i,j,k) / PhysConst::mu0 / mag_Ms_yface_arr(i,j,k) / mag_Ms_yface_arr(i,j,k); + + amrex::Real Ms_lo_x = mag_Ms_yface_arr(i-1, j, k); + amrex::Real Ms_hi_x = mag_Ms_yface_arr(i+1, j, k); + amrex::Real Ms_lo_y = mag_Ms_yface_arr(i, j-1, k); + amrex::Real Ms_hi_y = mag_Ms_yface_arr(i, j+1, k); + amrex::Real Ms_lo_z = mag_Ms_yface_arr(i, j, k-1); + amrex::Real Ms_hi_z = mag_Ms_yface_arr(i, j, k+1); + + Hx1_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M1_yface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 0, 1); //Last argument is nodality -- yface = 1 + Hy1_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M1_yface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 1, 1); //Last argument is nodality -- yface = 1 + Hz1_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M1_yface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 2, 1); //Last argument is nodality -- yface = 1 + + Hx2_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M2_yface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 0, 1); //Last argument is nodality -- yface = 1 + Hy2_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M2_yface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 1, 1); //Last argument is nodality -- yface = 1 + Hz2_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M2_yface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 2, 1); //Last argument is nodality -- yface = 1 + } + + //mag_anisotropy_coupling = 0; + if (mag_anisotropy_coupling == 1){ + + if (mag_anisotropy_yface_arr(i,j,k) == 0._rt) amrex::Abort("The mag_anisotropy_yface_arr(i,j,k) is 0.0 while including the anisotropy coupling term H_anisotropy for H_eff"); + + // H_anisotropy - use M^(old_time) + amrex::Real M_dot_anisotropy_axis = 0.0; + for (int comp=0; comp<3; ++comp) { + M_dot_anisotropy_axis += M_yface(i, j, k, comp) * anisotropy_axis[comp]; + } + amrex::Real const H_anisotropy_coeff = - 2.0 * mag_anisotropy_yface_arr(i,j,k) / PhysConst::mu0 / mag_Ms_yface_arr(i,j,k) / mag_Ms_yface_arr(i,j,k); + Hx1_eff += H_anisotropy_coeff * M_dot_anisotropy_axis * anisotropy_axis[0]; + Hy1_eff += H_anisotropy_coeff * M_dot_anisotropy_axis * anisotropy_axis[1]; + Hz1_eff += H_anisotropy_coeff * M_dot_anisotropy_axis * anisotropy_axis[2]; + } + + if(mag_afm_coupling == 1){ + if (mag_afm_yface_arr(i,j,k) == 0._rt) amrex::Abort("The mag_afm_yface_arr(i,j,k) is 0.0 while including the afm coupling term H_afm for H_eff"); + + // H_afm - use M^(old_time) + amrex::Real const H_afm_coeff = 4.0 * mag_afm_yface_arr(i,j,k) / PhysConst::mu0 / mag_Ms_yface_arr(i,j,k) / mag_Ms_yface_arr(i,j,k); + Hx1_eff += H_afm_coeff * M2_yface(i,j,k,0); + Hy1_eff += H_afm_coeff * M2_yface(i,j,k,1); + Hz1_eff += H_afm_coeff * M2_yface(i,j,k,2); + + Hx2_eff += H_afm_coeff * M1_yface(i,j,k,0); + Hy2_eff += H_afm_coeff * M1_yface(i,j,k,1); + Hz2_eff += H_afm_coeff * M1_yface(i,j,k,2); + } + + // 0 = unsaturated; compute |M| locally. 1 = saturated; use M_s + // note the unsaturated case is less usefull in real devices + amrex::Real M1_magnitude = (M_normalization == 0) ? std::sqrt(std::pow(M1_yface(i, j, k, 0), 2._rt) + std::pow(M1_yface(i, j, k, 1), 2._rt) + std::pow(M1_yface(i, j, k, 2), 2._rt)) + : mag_Ms_yface_arr(i,j,k); + amrex::Real M2_magnitude = (M_normalization == 0) ? std::sqrt(std::pow(M2_yface(i, j, k, 0), 2._rt) + std::pow(M2_yface(i, j, k, 1), 2._rt) + std::pow(M2_yface(i, j, k, 2), 2._rt)) + : mag_Ms_yface_arr(i,j,k); + amrex::Real a1_temp_static_coeff = mag_alpha_yface_arr(i,j,k) / M1_magnitude; + amrex::Real a2_temp_static_coeff = mag_alpha_yface_arr(i,j,k) / M2_magnitude; + + // calculate the b_temp_static_coeff (it is divided by 2.0 because the derivation is based on an interger dt, + // while in real simulations, the input dt is actually dt/2.0) + amrex::Real b_temp_static_coeff = - PhysConst::mu0 * amrex::Math::abs(mag_gamma_yface_arr(i,j,k)) / 2._rt; + + for (int comp=0; comp<3; ++comp) { + // calculate a_temp_static_yface + // all component on y-faces of grid + a1_temp_static_yface(i, j, k, comp) = a1_temp_static_coeff * M1_yface(i, j, k, comp); + a2_temp_static_yface(i, j, k, comp) = a2_temp_static_coeff * M2_yface(i, j, k, comp); + } + + // calculate b_temp_static_yface + // x component on y-faces of grid + b1_temp_static_yface(i, j, k, 0) = M1_yface(i, j, k, 0) + dt * b_temp_static_coeff * (M1_yface(i, j, k, 1) * Hz1_eff - M1_yface(i, j, k, 2) * Hy1_eff); + b2_temp_static_yface(i, j, k, 0) = M2_yface(i, j, k, 0) + dt * b_temp_static_coeff * (M2_yface(i, j, k, 1) * Hz2_eff - M2_yface(i, j, k, 2) * Hy2_eff); + + // y component on y-faces of grid + b1_temp_static_yface(i, j, k, 1) = M1_yface(i, j, k, 1) + dt * b_temp_static_coeff * (M1_yface(i, j, k, 2) * Hx1_eff - M1_yface(i, j, k, 0) * Hz1_eff); + b2_temp_static_yface(i, j, k, 1) = M2_yface(i, j, k, 1) + dt * b_temp_static_coeff * (M2_yface(i, j, k, 2) * Hx2_eff - M2_yface(i, j, k, 0) * Hz2_eff); + + // z component on y-faces of grid + b1_temp_static_yface(i, j, k, 2) = M1_yface(i, j, k, 2) + dt * b_temp_static_coeff * (M1_yface(i, j, k, 0) * Hy1_eff - M1_yface(i, j, k, 1) * Hx1_eff); + b2_temp_static_yface(i, j, k, 2) = M2_yface(i, j, k, 2) + dt * b_temp_static_coeff * (M2_yface(i, j, k, 0) * Hy2_eff - M2_yface(i, j, k, 1) * Hx2_eff); + } + }); + + amrex::ParallelFor(tbz, + [=] AMREX_GPU_DEVICE(int i, int j, int k) { + + // determine if the material is nonmagnetic or not + if (mag_Ms_zface_arr(i,j,k) > 0._rt){ + + // when working on M_zface(i,j,k,0:2) we have direct access to M_zface(i,j,k,0:2) and Hz(i,j,k) + // Hy and Hz can be acquired by interpolation + + // H_bias + amrex::Real Hx1_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mxface_stag, Mzface_stag, Hx_bias); + amrex::Real Hy1_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Myface_stag, Mzface_stag, Hy_bias); + amrex::Real Hz1_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mzface_stag, Mzface_stag, Hz_bias); + + amrex::Real Hx2_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mxface_stag, Mzface_stag, Hx_bias); + amrex::Real Hy2_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Myface_stag, Mzface_stag, Hy_bias); + amrex::Real Hz2_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mzface_stag, Mzface_stag, Hz_bias); + + if (coupling == 1){ + // H_eff = H_maxwell + H_bias + H_exchange + H_anisotropy + + // H_maxwell - use H^(old_time) + Hx1_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mxface_stag, Mzface_stag, Hx_old); + Hy1_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Myface_stag, Mzface_stag, Hy_old); + Hz1_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mzface_stag, Mzface_stag, Hz_old); + + Hx2_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mxface_stag, Mzface_stag, Hx_old); + Hy2_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Myface_stag, Mzface_stag, Hy_old); + Hz2_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Mzface_stag, Mzface_stag, Hz_old); + } + + if (mag_exchange_coupling == 1){ + + if (mag_exchange_zface_arr(i,j,k) == 0._rt) amrex::Abort("The mag_exchange_zface_arr(i,j,k) is 0.0 while including the exchange coupling term H_exchange for H_eff"); + + // H_exchange - use M^(old_time) + amrex::Real const H_exchange_coeff = 2.0 * mag_exchange_zface_arr(i,j,k) / PhysConst::mu0 / mag_Ms_zface_arr(i,j,k) / mag_Ms_zface_arr(i,j,k); + + amrex::Real Ms_lo_x = mag_Ms_zface_arr(i-1, j, k); + amrex::Real Ms_hi_x = mag_Ms_zface_arr(i+1, j, k); + amrex::Real Ms_lo_y = mag_Ms_zface_arr(i, j-1, k); + amrex::Real Ms_hi_y = mag_Ms_zface_arr(i, j+1, k); + amrex::Real Ms_lo_z = mag_Ms_zface_arr(i, j, k-1); + amrex::Real Ms_hi_z = mag_Ms_zface_arr(i, j, k+1); + + Hx1_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M1_zface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 0, 2); //Last argument is nodality -- zface = 2 + Hy1_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M1_zface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 1, 2); //Last argument is nodality -- zface = 2 + Hz1_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M1_zface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 2, 2); //Last argument is nodality -- zface = 2 + + Hx2_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M2_zface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 0, 2); //Last argument is nodality -- zface = 2 + Hy2_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M2_zface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 1, 2); //Last argument is nodality -- zface = 2 + Hz2_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M2_zface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 2, 2); //Last argument is nodality -- zface = 2 + } + + //mag_anisotropy_coupling = 0; + if (mag_anisotropy_coupling == 1){ + + if (mag_anisotropy_zface_arr(i,j,k) == 0._rt) amrex::Abort("The mag_anisotropy_zface_arr(i,j,k) is 0.0 while including the anisotropy coupling term H_anisotropy for H_eff"); + + // H_anisotropy - use M^(old_time) + amrex::Real M_dot_anisotropy_axis = 0.0; + for (int comp=0; comp<3; ++comp) { + M_dot_anisotropy_axis += M_zface(i, j, k, comp) * anisotropy_axis[comp]; + } + amrex::Real const H_anisotropy_coeff = - 2.0 * mag_anisotropy_zface_arr(i,j,k) / PhysConst::mu0 / mag_Ms_zface_arr(i,j,k) / mag_Ms_zface_arr(i,j,k); + Hx1_eff += H_anisotropy_coeff * M_dot_anisotropy_axis * anisotropy_axis[0]; + Hy1_eff += H_anisotropy_coeff * M_dot_anisotropy_axis * anisotropy_axis[1]; + Hz1_eff += H_anisotropy_coeff * M_dot_anisotropy_axis * anisotropy_axis[2]; + } + + if(mag_afm_coupling == 1){ + if (mag_afm_zface_arr(i,j,k) == 0._rt) amrex::Abort("The mag_afm_zface_arr(i,j,k) is 0.0 while including the afm coupling term H_afm for H_eff"); + + // H_afm - use M^(old_time) + amrex::Real const H_afm_coeff = 4.0 * mag_afm_zface_arr(i,j,k) / PhysConst::mu0 / mag_Ms_zface_arr(i,j,k) / mag_Ms_zface_arr(i,j,k); + Hx1_eff += H_afm_coeff * M2_zface(i,j,k,0); + Hy1_eff += H_afm_coeff * M2_zface(i,j,k,1); + Hz1_eff += H_afm_coeff * M2_zface(i,j,k,2); + + Hx2_eff += H_afm_coeff * M1_zface(i,j,k,0); + Hy2_eff += H_afm_coeff * M1_zface(i,j,k,1); + Hz2_eff += H_afm_coeff * M1_zface(i,j,k,2); + } + + // 0 = unsaturated; compute |M| locally. 1 = saturated; use M_s + amrex::Real M1_magnitude = (M_normalization == 0) ? std::sqrt(std::pow(M1_zface(i, j, k, 0), 2._rt) + std::pow(M1_zface(i, j, k, 1), 2._rt) + std::pow(M1_zface(i, j, k, 2), 2._rt)) + : mag_Ms_zface_arr(i,j,k); + amrex::Real M2_magnitude = (M_normalization == 0) ? std::sqrt(std::pow(M2_zface(i, j, k, 0), 2._rt) + std::pow(M2_zface(i, j, k, 1), 2._rt) + std::pow(M2_zface(i, j, k, 2), 2._rt)) + : mag_Ms_zface_arr(i,j,k); + amrex::Real a1_temp_static_coeff = mag_alpha_zface_arr(i,j,k) / M1_magnitude; + amrex::Real a2_temp_static_coeff = mag_alpha_zface_arr(i,j,k) / M2_magnitude; + + // calculate the b_temp_static_coeff (it is divided by 2.0 because the derivation is based on an interger dt, + // while in real simulations, the input dt is actually dt/2.0) + amrex::Real b_temp_static_coeff = - PhysConst::mu0 * amrex::Math::abs(mag_gamma_zface_arr(i,j,k)) / 2._rt; + + for (int comp=0; comp<3; ++comp) { + // calculate a_temp_static_zface + // all components on z-faces of grid + a1_temp_static_zface(i, j, k, comp) = a1_temp_static_coeff * M1_zface(i, j, k, comp); + a2_temp_static_zface(i, j, k, comp) = a2_temp_static_coeff * M2_zface(i, j, k, comp); + } + + // calculate b_temp_static_zface + // x component on z-faces of grid + b1_temp_static_zface(i, j, k, 0) = M1_zface(i, j, k, 0) + dt * b_temp_static_coeff * (M1_zface(i, j, k, 1) * Hz1_eff - M1_zface(i, j, k, 2) * Hy1_eff); + b2_temp_static_zface(i, j, k, 0) = M2_zface(i, j, k, 0) + dt * b_temp_static_coeff * (M2_zface(i, j, k, 1) * Hz2_eff - M2_zface(i, j, k, 2) * Hy2_eff); + + // y component on z-faces of grid + b1_temp_static_zface(i, j, k, 1) = M1_zface(i, j, k, 1) + dt * b_temp_static_coeff * (M1_zface(i, j, k, 2) * Hx1_eff - M1_zface(i, j, k, 0) * Hz1_eff); + b2_temp_static_zface(i, j, k, 1) = M2_zface(i, j, k, 1) + dt * b_temp_static_coeff * (M2_zface(i, j, k, 2) * Hx2_eff - M2_zface(i, j, k, 0) * Hz2_eff); + + // z component on z-faces of grid + b1_temp_static_zface(i, j, k, 2) = M1_zface(i, j, k, 2) + dt * b_temp_static_coeff * (M1_zface(i, j, k, 0) * Hy1_eff - M1_zface(i, j, k, 1) * Hx1_eff); + b2_temp_static_zface(i, j, k, 2) = M2_zface(i, j, k, 2) + dt * b_temp_static_coeff * (M2_zface(i, j, k, 0) * Hy2_eff - M2_zface(i, j, k, 1) * Hx2_eff); + } + }); + } + + // initialize M_max_iter, M_iter, M_tol, M_iter_error + // maximum number of iterations allowed + int M_max_iter = macroscopic_properties->getmag_max_iter(); + int M_iter = 0; + // relative tolerance stopping criteria for 2nd-order iterative algorithm + amrex::Real M_tol = macroscopic_properties->getmag_tol(); + int stop_iter = 0; + + // begin the iteration + while (!stop_iter){ + + warpx.FillBoundaryH(warpx.getngEB()); + + for (MFIter mfi(*Mfield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi){ + + auto& mag_Ms_xface_mf = macroscopic_properties->getmag_Ms_mf(0); + auto& mag_Ms_yface_mf = macroscopic_properties->getmag_Ms_mf(1); + auto& mag_Ms_zface_mf = macroscopic_properties->getmag_Ms_mf(2); + auto& mag_alpha_xface_mf = macroscopic_properties->getmag_alpha_mf(0); + auto& mag_alpha_yface_mf = macroscopic_properties->getmag_alpha_mf(1); + auto& mag_alpha_zface_mf = macroscopic_properties->getmag_alpha_mf(2); + auto& mag_gamma_xface_mf = macroscopic_properties->getmag_gamma_mf(0); + auto& mag_gamma_yface_mf = macroscopic_properties->getmag_gamma_mf(1); + auto& mag_gamma_zface_mf = macroscopic_properties->getmag_gamma_mf(2); + auto& mag_exchange_xface_mf = macroscopic_properties->getmag_exchange_mf(0); + auto& mag_exchange_yface_mf = macroscopic_properties->getmag_exchange_mf(1); + auto& mag_exchange_zface_mf = macroscopic_properties->getmag_exchange_mf(2); + auto& mag_anisotropy_xface_mf = macroscopic_properties->getmag_anisotropy_mf(0); + auto& mag_anisotropy_yface_mf = macroscopic_properties->getmag_anisotropy_mf(1); + auto& mag_anisotropy_zface_mf = macroscopic_properties->getmag_anisotropy_mf(2); + auto& mag_afm_xface_mf = macroscopic_properties->getmag_afm_mf(0); + auto& mag_afm_yface_mf = macroscopic_properties->getmag_afm_mf(1); + auto& mag_afm_zface_mf = macroscopic_properties->getmag_afm_mf(2); + + // extract material properties + Array4 const& mag_Ms_xface_arr = mag_Ms_xface_mf.array(mfi); + Array4 const& mag_Ms_yface_arr = mag_Ms_yface_mf.array(mfi); + Array4 const& mag_Ms_zface_arr = mag_Ms_zface_mf.array(mfi); + Array4 const& mag_alpha_xface_arr = mag_alpha_xface_mf.array(mfi); + Array4 const& mag_alpha_yface_arr = mag_alpha_yface_mf.array(mfi); + Array4 const& mag_alpha_zface_arr = mag_alpha_zface_mf.array(mfi); + Array4 const& mag_gamma_xface_arr = mag_gamma_xface_mf.array(mfi); + Array4 const& mag_gamma_yface_arr = mag_gamma_yface_mf.array(mfi); + Array4 const& mag_gamma_zface_arr = mag_gamma_zface_mf.array(mfi); + Array4 const& mag_exchange_xface_arr = mag_exchange_xface_mf.array(mfi); + Array4 const& mag_exchange_yface_arr = mag_exchange_yface_mf.array(mfi); + Array4 const& mag_exchange_zface_arr = mag_exchange_zface_mf.array(mfi); + Array4 const& mag_anisotropy_xface_arr = mag_anisotropy_xface_mf.array(mfi); + Array4 const& mag_anisotropy_yface_arr = mag_anisotropy_yface_mf.array(mfi); + Array4 const& mag_anisotropy_zface_arr = mag_anisotropy_zface_mf.array(mfi); + Array4 const& mag_afm_xface_arr = mag_afm_xface_mf.array(mfi); + Array4 const& mag_afm_yface_arr = mag_afm_yface_mf.array(mfi); + Array4 const& mag_afm_zface_arr = mag_afm_zface_mf.array(mfi); + + // extract field data + Array4 const &M_xface = Mfield[0]->array(mfi); // note M_xface include x,y,z components at |_x faces + Array4 const &M_yface = Mfield[1]->array(mfi); // note M_yface include x,y,z components at |_y faces + Array4 const &M_zface = Mfield[2]->array(mfi); // note M_zface include x,y,z components at |_z faces + Array4 const &M1_xface = M1field[0]->array(mfi); // note M_xface include x,y,z components at |_x faces + Array4 const &M1_yface = M1field[1]->array(mfi); // note M_yface include x,y,z components at |_y faces + Array4 const &M1_zface = M1field[2]->array(mfi); // note M_zface include x,y,z components at |_z faces + Array4 const &M2_xface = M2field[0]->array(mfi); // note M_xface include x,y,z components at |_x faces + Array4 const &M2_yface = M2field[1]->array(mfi); // note M_yface include x,y,z components at |_y faces + Array4 const &M2_zface = M2field[2]->array(mfi); // note M_zface include x,y,z components at |_z faces + Array4 const &Hx_bias = H_biasfield[0]->array(mfi); // Hx_bias is the x component at |_x faces + Array4 const &Hy_bias = H_biasfield[1]->array(mfi); // Hy_bias is the y component at |_y faces + Array4 const &Hz_bias = H_biasfield[2]->array(mfi); // Hz_bias is the z component at |_z faces + Array4 const &Hx = Hfield[0]->array(mfi); // Hx is the x component at |_x faces + Array4 const &Hy = Hfield[1]->array(mfi); // Hy is the y component at |_y faces + Array4 const &Hz = Hfield[2]->array(mfi); // Hz is the z component at |_z faces + + // extract field data of Mfield_prev, Mfield_error, a_temp, a_temp_static, and b_temp_static + Array4 const &M_prev_xface = Mfield_prev[0]->array(mfi); + Array4 const &M_prev_yface = Mfield_prev[1]->array(mfi); + Array4 const &M_prev_zface = Mfield_prev[2]->array(mfi); + Array4 const &M1_prev_xface = M1field_prev[0]->array(mfi); + Array4 const &M1_prev_yface = M1field_prev[1]->array(mfi); + Array4 const &M1_prev_zface = M1field_prev[2]->array(mfi); + Array4 const &M2_prev_xface = M2field_prev[0]->array(mfi); + Array4 const &M2_prev_yface = M2field_prev[1]->array(mfi); + Array4 const &M2_prev_zface = M2field_prev[2]->array(mfi); + Array4 const &M_old_xface = Mfield_old[0]->array(mfi); + Array4 const &M_old_yface = Mfield_old[1]->array(mfi); + Array4 const &M_old_zface = Mfield_old[2]->array(mfi); + Array4 const &M1_old_xface = M1field_old[0]->array(mfi); + Array4 const &M1_old_yface = M1field_old[1]->array(mfi); + Array4 const &M1_old_zface = M1field_old[2]->array(mfi); + Array4 const &M2_old_xface = M2field_old[0]->array(mfi); + Array4 const &M2_old_yface = M2field_old[1]->array(mfi); + Array4 const &M2_old_zface = M2field_old[2]->array(mfi); + Array4 const &M_error_xface = Mfield_error[0]->array(mfi); + Array4 const &M_error_yface = Mfield_error[1]->array(mfi); + Array4 const &M_error_zface = Mfield_error[2]->array(mfi); + Array4 const &M2_error_xface = M2field_error[0]->array(mfi); + Array4 const &M2_error_yface = M2field_error[1]->array(mfi); + Array4 const &M2_error_zface = M2field_error[2]->array(mfi); + Array4 const &M1_error_xface = M1field_error[0]->array(mfi); + Array4 const &M1_error_yface = M1field_error[1]->array(mfi); + Array4 const &M1_error_zface = M1field_error[2]->array(mfi); + Array4 const &a1_temp_xface = a1_temp[0]->array(mfi); + Array4 const &a1_temp_yface = a1_temp[1]->array(mfi); + Array4 const &a1_temp_zface = a1_temp[2]->array(mfi); + Array4 const &a2_temp_xface = a2_temp[0]->array(mfi); + Array4 const &a2_temp_yface = a2_temp[1]->array(mfi); + Array4 const &a2_temp_zface = a2_temp[2]->array(mfi); + Array4 const &a1_temp_static_xface = a1_temp_static[0]->array(mfi); + Array4 const &a1_temp_static_yface = a1_temp_static[1]->array(mfi); + Array4 const &a1_temp_static_zface = a1_temp_static[2]->array(mfi); + Array4 const &a2_temp_static_xface = a2_temp_static[0]->array(mfi); + Array4 const &a2_temp_static_yface = a2_temp_static[1]->array(mfi); + Array4 const &a2_temp_static_zface = a2_temp_static[2]->array(mfi); + Array4 const &b1_temp_static_xface = b1_temp_static[0]->array(mfi); + Array4 const &b1_temp_static_yface = b1_temp_static[1]->array(mfi); + Array4 const &b1_temp_static_zface = b1_temp_static[2]->array(mfi); + Array4 const &b2_temp_static_xface = b2_temp_static[0]->array(mfi); + Array4 const &b2_temp_static_yface = b2_temp_static[1]->array(mfi); + Array4 const &b2_temp_static_zface = b2_temp_static[2]->array(mfi); + + // extract tileboxes for which to loop + amrex::IntVect Hxnodal = Hfield[0]->ixType().toIntVect(); + amrex::IntVect Hynodal = Hfield[1]->ixType().toIntVect(); + amrex::IntVect Hznodal = Hfield[2]->ixType().toIntVect(); + Box const &tbx = mfi.tilebox(Hxnodal); /* just define which grid type */ + Box const &tby = mfi.tilebox(Hynodal); + Box const &tbz = mfi.tilebox(Hznodal); + + // Extract stencil coefficients for calculating the exchange field H_exchange and the anisotropy field H_anisotropy + amrex::Real const *const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr(); + int const n_coefs_x = m_stencil_coefs_x.size(); + amrex::Real const *const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr(); + int const n_coefs_y = m_stencil_coefs_y.size(); + amrex::Real const *const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr(); + int const n_coefs_z = m_stencil_coefs_z.size(); + + // loop over cells and update fields + amrex::ParallelFor(tbx, + [=] AMREX_GPU_DEVICE(int i, int j, int k) { + + // determine if the material is nonmagnetic or not + if (mag_Ms_xface_arr(i,j,k) > 0._rt){ + + // when working on M_xface(i,j,k, 0:2) we have direct access to M_xface(i,j,k,0:2) and Hx(i,j,k) + // Hy and Hz can be acquired by interpolation + + // H_bias + amrex::Real Hx1_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hxnodal, Hxnodal, Hx_bias); + amrex::Real Hy1_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hynodal, Hxnodal, Hy_bias); + amrex::Real Hz1_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hznodal, Hxnodal, Hz_bias); + + amrex::Real Hx2_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hxnodal, Hxnodal, Hx_bias); + amrex::Real Hy2_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hynodal, Hxnodal, Hy_bias); + amrex::Real Hz2_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hznodal, Hxnodal, Hz_bias); + + if (coupling == 1){ + // H_eff = H_maxwell + H_bias + H_exchange + H_anisotropy + + // H_maxwell - use H^[(new_time),r-1] + Hx1_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hxnodal, Hxnodal, Hx); + Hy1_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hynodal, Hxnodal, Hy); + Hz1_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hznodal, Hxnodal, Hz); + + Hx2_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hxnodal, Hxnodal, Hx); + Hy2_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hynodal, Hxnodal, Hy); + Hz2_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hznodal, Hxnodal, Hz); + } + + if (mag_exchange_coupling == 1){ + + if (mag_exchange_xface_arr(i,j,k) == 0._rt) amrex::Abort("The mag_exchange_xface_arr(i,j,k) is 0.0 while including the exchange coupling term H_exchange for H_eff"); + + // H_exchange - use M^[(new_time),r-1] + amrex::Real const H_exchange_coeff = 2.0 * mag_exchange_xface_arr(i,j,k) / PhysConst::mu0 / mag_Ms_xface_arr(i,j,k) / mag_Ms_xface_arr(i,j,k); + + amrex::Real Ms_lo_x = mag_Ms_xface_arr(i-1, j, k); + amrex::Real Ms_hi_x = mag_Ms_xface_arr(i+1, j, k); + amrex::Real Ms_lo_y = mag_Ms_xface_arr(i, j-1, k); + amrex::Real Ms_hi_y = mag_Ms_xface_arr(i, j+1, k); + amrex::Real Ms_lo_z = mag_Ms_xface_arr(i, j, k-1); + amrex::Real Ms_hi_z = mag_Ms_xface_arr(i, j, k+1); + + Hx1_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M1_prev_xface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 0, 0); //Last argument is nodality -- xface = 0 + Hy1_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M1_prev_xface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 1, 0); //Last argument is nodality -- xface = 0 + Hz1_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M1_prev_xface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 2, 0); //Last argument is nodality -- xface = 0 + + Hx2_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M2_prev_xface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 0, 0); //Last argument is nodality -- xface = 0 + Hy2_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M2_prev_xface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 1, 0); //Last argument is nodality -- xface = 0 + Hz2_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M2_prev_xface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 2, 0); //Last argument is nodality -- xface = 0 + } + + + //mag_anisotropy_coupling = 0; + if (mag_anisotropy_coupling == 1){ + + if (mag_anisotropy_xface_arr(i,j,k) == 0._rt) amrex::Abort("The mag_anisotropy_xface_arr(i,j,k) is 0.0 while including the anisotropy coupling term H_anisotropy for H_eff"); + + // H_anisotropy - use M^[(new_time),r-1] + amrex::Real M_dot_anisotropy_axis = 0.0; + for (int comp=0; comp<3; ++comp) { + M_dot_anisotropy_axis += M_xface(i, j, k, comp) * anisotropy_axis[comp]; + } + amrex::Real const H_anisotropy_coeff = - 2.0 * mag_anisotropy_xface_arr(i,j,k) / PhysConst::mu0 / mag_Ms_xface_arr(i,j,k) / mag_Ms_xface_arr(i,j,k); + Hx1_eff += H_anisotropy_coeff * M_dot_anisotropy_axis * anisotropy_axis[0]; + Hy1_eff += H_anisotropy_coeff * M_dot_anisotropy_axis * anisotropy_axis[1]; + Hz1_eff += H_anisotropy_coeff * M_dot_anisotropy_axis * anisotropy_axis[2]; + } + + if(mag_afm_coupling == 1){ + if (mag_afm_xface_arr(i,j,k) == 0._rt) amrex::Abort("The mag_afm_xface_arr(i,j,k) is 0.0 while including the afm coupling term H_afm for H_eff"); + + // H_afm - use M^(old_time) + amrex::Real const H_afm_coeff = 4.0 * mag_afm_xface_arr(i,j,k) / PhysConst::mu0 / mag_Ms_xface_arr(i,j,k) / mag_Ms_xface_arr(i,j,k); + Hx1_eff += H_afm_coeff * M2_xface(i,j,k,0); + Hy1_eff += H_afm_coeff * M2_xface(i,j,k,1); + Hz1_eff += H_afm_coeff * M2_xface(i,j,k,2); + + Hx2_eff += H_afm_coeff * M1_xface(i,j,k,0); + Hy2_eff += H_afm_coeff * M1_xface(i,j,k,1); + Hz2_eff += H_afm_coeff * M1_xface(i,j,k,2); + } + + // calculate the a_temp_dynamic_coeff (it is divided by 2.0 because the derivation is based on an interger dt, + // while in real simulations, the input dt is actually dt/2.0) + amrex::Real a_temp_dynamic_coeff = PhysConst::mu0 * amrex::Math::abs(mag_gamma_xface_arr(i,j,k)) / 2._rt; //not dependent on sublattice + + amrex::GpuArray H1_eff; + H1_eff[0] = Hx1_eff; + H1_eff[1] = Hy1_eff; + H1_eff[2] = Hz1_eff; + + amrex::GpuArray H2_eff; + H2_eff[0] = Hx2_eff; + H2_eff[1] = Hy2_eff; + H2_eff[2] = Hz2_eff; + + for (int comp=0; comp<3; ++comp) { + // calculate a_temp_xface + // all components on x-faces of grid + a1_temp_xface(i, j, k, comp) = (M_normalization != 0) ? -(dt * a_temp_dynamic_coeff * H1_eff[comp] + a1_temp_static_xface(i, j, k, comp)) + : -(dt * a_temp_dynamic_coeff * H1_eff[comp] + 0.5 * a1_temp_static_xface(i, j, k, comp) + + 0.5 * mag_alpha_xface_arr(i,j,k) * 1. / std::sqrt(std::pow(M1_xface(i, j, k, 0), 2._rt) + std::pow(M1_xface(i, j, k, 1), 2._rt) + std::pow(M1_xface(i, j, k, 2), 2._rt)) * M1_old_xface(i, j, k, comp)); + a2_temp_xface(i, j, k, comp) = (M_normalization != 0) ? -(dt * a_temp_dynamic_coeff * H2_eff[comp] + a2_temp_static_xface(i, j, k, comp)) + : -(dt * a_temp_dynamic_coeff * H2_eff[comp] + 0.5 * a2_temp_static_xface(i, j, k, comp) + + 0.5 * mag_alpha_xface_arr(i,j,k) * 1. / std::sqrt(std::pow(M2_xface(i, j, k, 0), 2._rt) + std::pow(M2_xface(i, j, k, 1), 2._rt) + std::pow(M2_xface(i, j, k, 2), 2._rt)) * M2_old_xface(i, j, k, comp)); + } + + for (int comp=0; comp<3; ++comp) { + // update M_xface from a and b using the updateM_field + // all components on x-faces of grid + M1_xface(i, j, k, comp) = MacroscopicProperties::updateM_field(i, j, k, comp, a1_temp_xface, b1_temp_static_xface); + M2_xface(i, j, k, comp) = MacroscopicProperties::updateM_field(i, j, k, comp, a2_temp_xface, b2_temp_static_xface); + } + + // temporary normalized magnitude of M_xface field at the fixed point + // re-investigate the way we do Ms interp, in case we encounter the case where Ms changes across two adjacent cells that you are doing interp + amrex::Real M1_magnitude_normalized = std::sqrt(std::pow(M1_xface(i, j, k, 0), 2._rt) + std::pow(M1_xface(i, j, k, 1), 2._rt) + std::pow(M1_xface(i, j, k, 2), 2._rt)) / mag_Ms_xface_arr(i,j,k); + amrex::Real M2_magnitude_normalized = std::sqrt(std::pow(M2_xface(i, j, k, 0), 2._rt) + std::pow(M2_xface(i, j, k, 1), 2._rt) + std::pow(M2_xface(i, j, k, 2), 2._rt)) / mag_Ms_xface_arr(i,j,k); + if (M_normalization == 1){ + // saturated case; if |M| has drifted from M_s too much, abort. Otherwise, normalize + // check the normalized error + if (amrex::Math::abs(1._rt - M1_magnitude_normalized) > mag_normalized_error){ + amrex::Abort("Exceed the normalized error of the M_xface field"); + } + if (amrex::Math::abs(1._rt - M2_magnitude_normalized) > mag_normalized_error){ + amrex::Abort("Exceed the normalized error of the M_xface field"); + } + // normalize the M_xface field + M1_xface(i, j, k, 0) /= M1_magnitude_normalized; + M1_xface(i, j, k, 1) /= M1_magnitude_normalized; + M1_xface(i, j, k, 2) /= M1_magnitude_normalized; + + M2_xface(i, j, k, 0) /= M2_magnitude_normalized; + M2_xface(i, j, k, 1) /= M2_magnitude_normalized; + M2_xface(i, j, k, 2) /= M2_magnitude_normalized; + } + else if (M_normalization == 0){ + // check the normalized error + if (M1_magnitude_normalized > (1._rt + mag_normalized_error)){ + amrex::Abort("Caution: Unsaturated material has M_xface exceeding the saturation magnetization"); + } + else if (M1_magnitude_normalized > 1._rt && M1_magnitude_normalized <= 1._rt + mag_normalized_error){ + // normalize the M_xface field + M1_xface(i, j, k, 0) /= M1_magnitude_normalized; + M1_xface(i, j, k, 1) /= M1_magnitude_normalized; + M1_xface(i, j, k, 2) /= M1_magnitude_normalized; + } + if (M2_magnitude_normalized > (1._rt + mag_normalized_error)){ + amrex::Abort("Caution: Unsaturated material has M_xface exceeding the saturation magnetization"); + } + else if (M2_magnitude_normalized > 1._rt && M2_magnitude_normalized <= 1._rt + mag_normalized_error){ + // normalize the M_xface field + M2_xface(i, j, k, 0) /= M2_magnitude_normalized; + M2_xface(i, j, k, 1) /= M2_magnitude_normalized; + M2_xface(i, j, k, 2) /= M2_magnitude_normalized; + } + } + + // calculate M_error_xface + // x,y,z component on M-error on x-faces of grid + for (int icomp = 0; icomp < 3; ++icomp) { + M1_error_xface(i, j, k, icomp) = amrex::Math::abs((M1_xface(i, j, k, icomp) - M1_prev_xface(i, j, k, icomp))) / mag_Ms_xface_arr(i,j,k); + M2_error_xface(i, j, k, icomp) = amrex::Math::abs((M2_xface(i, j, k, icomp) - M2_prev_xface(i, j, k, icomp))) / mag_Ms_xface_arr(i,j,k); + } + } + }); + + amrex::ParallelFor(tby, + [=] AMREX_GPU_DEVICE(int i, int j, int k) { + + // determine if the material is nonmagnetic or not + if (mag_Ms_yface_arr(i,j,k) > 0._rt){ + + // when working on M_yface(i,j,k,0:2) we have direct access to M_yface(i,j,k,0:2) and Hy(i,j,k) + // Hy and Hz can be acquired by interpolation + + // H_bias + amrex::Real Hx1_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hxnodal, Hynodal, Hx_bias); + amrex::Real Hy1_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hynodal, Hynodal, Hy_bias); + amrex::Real Hz1_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hznodal, Hynodal, Hz_bias); + + amrex::Real Hx2_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hxnodal, Hynodal, Hx_bias); + amrex::Real Hy2_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hynodal, Hynodal, Hy_bias); + amrex::Real Hz2_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hznodal, Hynodal, Hz_bias); + + if (coupling == 1){ + // H_eff = H_maxwell + H_bias + H_exchange + H_anisotropy + + // H_maxwell - use H^[(new_time),r-1] + Hx1_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hxnodal, Hynodal, Hx); + Hy1_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hynodal, Hynodal, Hy); + Hz1_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hznodal, Hynodal, Hz); + + Hx2_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hxnodal, Hynodal, Hx); + Hy2_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hynodal, Hynodal, Hy); + Hz2_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hznodal, Hynodal, Hz); + } + + if (mag_exchange_coupling == 1){ + + if (mag_exchange_yface_arr(i,j,k) == 0._rt) amrex::Abort("The mag_exchange_yface_arr(i,j,k) is 0.0 while including the exchange coupling term H_exchange for H_eff"); + + // H_exchange - use M^[(new_time),r-1] + amrex::Real const H_exchange_coeff = 2.0 * mag_exchange_yface_arr(i,j,k) / PhysConst::mu0 / mag_Ms_yface_arr(i,j,k) / mag_Ms_yface_arr(i,j,k); + + amrex::Real Ms_lo_x = mag_Ms_yface_arr(i-1, j, k); + amrex::Real Ms_hi_x = mag_Ms_yface_arr(i+1, j, k); + amrex::Real Ms_lo_y = mag_Ms_yface_arr(i, j-1, k); + amrex::Real Ms_hi_y = mag_Ms_yface_arr(i, j+1, k); + amrex::Real Ms_lo_z = mag_Ms_yface_arr(i, j, k-1); + amrex::Real Ms_hi_z = mag_Ms_yface_arr(i, j, k+1); + + Hx1_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M1_prev_yface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 0, 1); //Last argument is nodality -- yface = 1 + Hy1_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M1_prev_yface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 1, 1); //Last argument is nodality -- yface = 1 + Hz1_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M1_prev_yface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 2, 1); //Last argument is nodality -- yface = 1 + + Hx2_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M2_prev_yface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 0, 1); //Last argument is nodality -- yface = 1 + Hy2_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M2_prev_yface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 1, 1); //Last argument is nodality -- yface = 1 + Hz2_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M2_prev_yface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 2, 1); //Last argument is nodality -- yface = 1 + } + + //mag_anisotropy_coupling = 0; + if (mag_anisotropy_coupling == 1){ + + if (mag_anisotropy_yface_arr(i,j,k) == 0._rt) amrex::Abort("The mag_anisotropy_yface_arr(i,j,k) is 0.0 while including the anisotropy coupling term H_anisotropy for H_eff"); + + // H_anisotropy - use M^[(new_time),r-1] + amrex::Real M_dot_anisotropy_axis = 0.0; + for (int comp=0; comp<3; ++comp) { + M_dot_anisotropy_axis += M_yface(i, j, k, comp) * anisotropy_axis[comp]; + } + amrex::Real const H_anisotropy_coeff = - 2.0 * mag_anisotropy_yface_arr(i,j,k) / PhysConst::mu0 / mag_Ms_yface_arr(i,j,k) / mag_Ms_yface_arr(i,j,k); + Hx1_eff += H_anisotropy_coeff * M_dot_anisotropy_axis * anisotropy_axis[0]; + Hy1_eff += H_anisotropy_coeff * M_dot_anisotropy_axis * anisotropy_axis[1]; + Hz1_eff += H_anisotropy_coeff * M_dot_anisotropy_axis * anisotropy_axis[2]; + } + + if(mag_afm_coupling == 1){ + if (mag_afm_yface_arr(i,j,k) == 0._rt) amrex::Abort("The mag_afm_yface_arr(i,j,k) is 0.0 while including the afm coupling term H_afm for H_eff"); + + // H_afm - use M^(old_time) + amrex::Real const H_afm_coeff = 4.0 * mag_afm_yface_arr(i,j,k) / PhysConst::mu0 / mag_Ms_yface_arr(i,j,k) / mag_Ms_yface_arr(i,j,k); + Hx1_eff += H_afm_coeff * M2_yface(i,j,k,0); + Hy1_eff += H_afm_coeff * M2_yface(i,j,k,1); + Hz1_eff += H_afm_coeff * M2_yface(i,j,k,2); + + Hx2_eff += H_afm_coeff * M1_yface(i,j,k,0); + Hy2_eff += H_afm_coeff * M1_yface(i,j,k,1); + Hz2_eff += H_afm_coeff * M1_yface(i,j,k,2); + } + + // calculate the a_temp_dynamic_coeff (it is divided by 2.0 because the derivation is based on an interger dt, + // while in real simulations, the input dt is actually dt/2.0) + amrex::Real a_temp_dynamic_coeff = PhysConst::mu0 * amrex::Math::abs(mag_gamma_yface_arr(i,j,k)) / 2._rt; + + amrex::GpuArray H1_eff; + H1_eff[0] = Hx1_eff; + H1_eff[1] = Hy1_eff; + H1_eff[2] = Hz1_eff; + + amrex::GpuArray H2_eff; + H2_eff[0] = Hx2_eff; + H2_eff[1] = Hy2_eff; + H2_eff[2] = Hz2_eff; + + for (int comp=0; comp<3; ++comp) { + // calculate a_temp_yface + // all components on y-faces of grid + a1_temp_yface(i, j, k, comp) = (M_normalization != 0) ? -(dt * a_temp_dynamic_coeff * H1_eff[comp] + a1_temp_static_yface(i, j, k, comp)) + : -(dt * a_temp_dynamic_coeff * H1_eff[comp] + 0.5 * a1_temp_static_yface(i, j, k, comp) + + 0.5 * mag_alpha_yface_arr(i,j,k) * 1. / std::sqrt(std::pow(M1_yface(i, j, k, 0), 2._rt) + std::pow(M1_yface(i, j, k, 1), 2._rt) + std::pow(M1_yface(i, j, k, 2), 2._rt)) * M1_old_yface(i, j, k, comp)); + a2_temp_yface(i, j, k, comp) = (M_normalization != 0) ? -(dt * a_temp_dynamic_coeff * H2_eff[comp] + a2_temp_static_yface(i, j, k, comp)) + : -(dt * a_temp_dynamic_coeff * H2_eff[comp] + 0.5 * a2_temp_static_yface(i, j, k, comp) + + 0.5 * mag_alpha_yface_arr(i,j,k) * 1. / std::sqrt(std::pow(M2_yface(i, j, k, 0), 2._rt) + std::pow(M2_yface(i, j, k, 1), 2._rt) + std::pow(M2_yface(i, j, k, 2), 2._rt)) * M2_old_yface(i, j, k, comp)); + } + + + for (int comp=0; comp<3; ++comp) { + // update M_yface from a and b using the updateM_field + // all components on y-faces of grid + M1_yface(i, j, k, comp) = MacroscopicProperties::updateM_field(i, j, k, comp, a1_temp_yface, b1_temp_static_yface); + M2_yface(i, j, k, comp) = MacroscopicProperties::updateM_field(i, j, k, comp, a2_temp_yface, b2_temp_static_yface); + } + + // temporary normalized magnitude of M_yface field at the fixed point + // re-investigate the way we do Ms interp, in case we encounter the case where Ms changes across two adjacent cells that you are doing interp + amrex::Real M1_magnitude_normalized = std::sqrt(std::pow(M1_yface(i, j, k, 0), 2._rt) + std::pow(M1_yface(i, j, k, 1), 2._rt) + std::pow(M1_yface(i, j, k, 2), 2._rt)) / mag_Ms_yface_arr(i,j,k); + amrex::Real M2_magnitude_normalized = std::sqrt(std::pow(M2_yface(i, j, k, 0), 2._rt) + std::pow(M2_yface(i, j, k, 1), 2._rt) + std::pow(M2_yface(i, j, k, 2), 2._rt)) / mag_Ms_yface_arr(i,j,k); + + if (M_normalization == 1){ + // saturated case; if |M| has drifted from M_s too much, abort. Otherwise, normalize + // check the normalized error + if (amrex::Math::abs(1._rt - M1_magnitude_normalized) > mag_normalized_error){ + amrex::Abort("Exceed the normalized error of the M_yface field"); + } + if (amrex::Math::abs(1._rt - M2_magnitude_normalized) > mag_normalized_error){ + amrex::Abort("Exceed the normalized error of the M_yface field"); + } + // normalize the M_yface field + M1_yface(i, j, k, 0) /= M1_magnitude_normalized; + M1_yface(i, j, k, 1) /= M1_magnitude_normalized; + M1_yface(i, j, k, 2) /= M1_magnitude_normalized; + + M2_yface(i, j, k, 0) /= M2_magnitude_normalized; + M2_yface(i, j, k, 1) /= M2_magnitude_normalized; + M2_yface(i, j, k, 2) /= M2_magnitude_normalized; + } + else if (M_normalization == 0){ + // check the normalized error + if (M1_magnitude_normalized > 1._rt + mag_normalized_error){ + amrex::Abort("Caution: Unsaturated material has M_yface exceeding the saturation magnetization"); + } + else if (M1_magnitude_normalized > 1._rt && M1_magnitude_normalized <= 1._rt + mag_normalized_error){ + // normalize the M_yface field + M1_yface(i, j, k, 0) /= M1_magnitude_normalized; + M1_yface(i, j, k, 1) /= M1_magnitude_normalized; + M1_yface(i, j, k, 2) /= M1_magnitude_normalized; + } + if (M2_magnitude_normalized > 1._rt + mag_normalized_error){ + amrex::Abort("Caution: Unsaturated material has M_yface exceeding the saturation magnetization"); + } + else if (M2_magnitude_normalized > 1._rt && M2_magnitude_normalized <= 1._rt + mag_normalized_error){ + // normalize the M_yface field + M2_yface(i, j, k, 0) /= M2_magnitude_normalized; + M2_yface(i, j, k, 1) /= M2_magnitude_normalized; + M2_yface(i, j, k, 2) /= M2_magnitude_normalized; + } + } + + // calculate M_error_yface + // x,y,z component on y-faces of grid + for (int icomp = 0; icomp < 3; ++icomp) { + M1_error_yface(i, j, k, icomp) = amrex::Math::abs((M1_yface(i, j, k, icomp) - M1_prev_yface(i, j, k, icomp))) / mag_Ms_yface_arr(i,j,k); + M2_error_yface(i, j, k, icomp) = amrex::Math::abs((M2_yface(i, j, k, icomp) - M2_prev_yface(i, j, k, icomp))) / mag_Ms_yface_arr(i,j,k); + } + } + }); + + amrex::ParallelFor(tbz, + [=] AMREX_GPU_DEVICE(int i, int j, int k) { + + // determine if the material is nonmagnetic or not + if (mag_Ms_zface_arr(i,j,k) > 0._rt){ + + // when working on M_zface(i,j,k,0:2) we have direct access to M_zface(i,j,k,0:2) and Hz(i,j,k) + // Hy and Hz can be acquired by interpolation + + // H_bias + amrex::Real Hx1_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hxnodal, Hznodal, Hx_bias); + amrex::Real Hy1_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hynodal, Hznodal, Hy_bias); + amrex::Real Hz1_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hznodal, Hznodal, Hz_bias); + + amrex::Real Hx2_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hxnodal, Hznodal, Hx_bias); + amrex::Real Hy2_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hynodal, Hznodal, Hy_bias); + amrex::Real Hz2_eff = MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hznodal, Hznodal, Hz_bias); + + if (coupling == 1){ + // H_eff = H_maxwell + H_bias + H_exchange + H_anisotropy + + // H_maxwell - use H^[(new_time),r-1] + Hx1_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hxnodal, Hznodal, Hx); + Hy1_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hynodal, Hznodal, Hy); + Hz1_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hznodal, Hznodal, Hz); + + Hx2_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hxnodal, Hznodal, Hx); + Hy2_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hynodal, Hznodal, Hy); + Hz2_eff += MacroscopicProperties::face_avg_to_face(i, j, k, 0, Hznodal, Hznodal, Hz); + } + + if (mag_exchange_coupling == 1){ + + if (mag_exchange_zface_arr(i,j,k) == 0._rt) amrex::Abort("The mag_exchange_zface_arr(i,j,k) is 0.0 while including the exchange coupling term H_exchange for H_eff"); + + // H_exchange - use M^[(new_time),r-1] + amrex::Real const H_exchange_coeff = 2.0 * mag_exchange_zface_arr(i,j,k) / PhysConst::mu0 / mag_Ms_zface_arr(i,j,k) / mag_Ms_zface_arr(i,j,k); + + amrex::Real Ms_lo_x = mag_Ms_zface_arr(i-1, j, k); + amrex::Real Ms_hi_x = mag_Ms_zface_arr(i+1, j, k); + amrex::Real Ms_lo_y = mag_Ms_zface_arr(i, j-1, k); + amrex::Real Ms_hi_y = mag_Ms_zface_arr(i, j+1, k); + amrex::Real Ms_lo_z = mag_Ms_zface_arr(i, j, k-1); + amrex::Real Ms_hi_z = mag_Ms_zface_arr(i, j, k+1); + + Hx1_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M1_prev_zface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 0, 2); //Last argument is nodality -- zface = 2 + Hy1_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M1_prev_zface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 1, 2); //Last argument is nodality -- zface = 2 + Hz1_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M1_prev_zface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 2, 2); //Last argument is nodality -- zface = 2 + + Hx2_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M2_prev_zface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 0, 2); //Last argument is nodality -- zface = 2 + Hy2_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M2_prev_zface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 1, 2); //Last argument is nodality -- zface = 2 + Hz2_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M2_prev_zface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 2, 2); //Last argument is nodality -- zface = 2 + } + + //mag_anisotropy_coupling = 0; + if (mag_anisotropy_coupling == 1){ + + if (mag_anisotropy_zface_arr(i,j,k) == 0._rt) amrex::Abort("The mag_anisotropy_zface_arr(i,j,k) is 0.0 while including the anisotropy coupling term H_anisotropy for H_eff"); + + // H_anisotropy - use M^[(new_time),r-1] + amrex::Real M_dot_anisotropy_axis = 0.0; + for (int comp=0; comp<3; ++comp) { + M_dot_anisotropy_axis += M_zface(i, j, k, comp) * anisotropy_axis[comp]; + } + amrex::Real const H_anisotropy_coeff = - 2.0 * mag_anisotropy_zface_arr(i,j,k) / PhysConst::mu0 / mag_Ms_zface_arr(i,j,k) / mag_Ms_zface_arr(i,j,k); + Hx1_eff += H_anisotropy_coeff * M_dot_anisotropy_axis * anisotropy_axis[0]; + Hy1_eff += H_anisotropy_coeff * M_dot_anisotropy_axis * anisotropy_axis[1]; + Hz1_eff += H_anisotropy_coeff * M_dot_anisotropy_axis * anisotropy_axis[2]; + } + + if(mag_afm_coupling == 1){ + if (mag_afm_zface_arr(i,j,k) == 0._rt) amrex::Abort("The mag_afm_zface_arr(i,j,k) is 0.0 while including the afm coupling term H_afm for H_eff"); + + // H_afm - use M^(old_time) + amrex::Real const H_afm_coeff = 4.0 * mag_afm_zface_arr(i,j,k) / PhysConst::mu0 / mag_Ms_zface_arr(i,j,k) / mag_Ms_zface_arr(i,j,k); + Hx1_eff += H_afm_coeff * M2_zface(i,j,k,0); + Hy1_eff += H_afm_coeff * M2_zface(i,j,k,1); + Hz1_eff += H_afm_coeff * M2_zface(i,j,k,2); + + Hx2_eff += H_afm_coeff * M1_zface(i,j,k,0); + Hy2_eff += H_afm_coeff * M1_zface(i,j,k,1); + Hz2_eff += H_afm_coeff * M1_zface(i,j,k,2); + } + + // calculate the a_temp_dynamic_coeff (it is divided by 2.0 because the derivation is based on an interger dt, + // while in real simulations, the input dt is actually dt/2.0) + amrex::Real a_temp_dynamic_coeff = PhysConst::mu0 * amrex::Math::abs(mag_gamma_zface_arr(i,j,k)) / 2._rt; + + amrex::GpuArray H1_eff; + H1_eff[0] = Hx1_eff; + H1_eff[1] = Hy1_eff; + H1_eff[2] = Hz1_eff; + + amrex::GpuArray H2_eff; + H2_eff[0] = Hx2_eff; + H2_eff[1] = Hy2_eff; + H2_eff[2] = Hz2_eff; + + for (int comp=0; comp<3; ++comp) { + // calculate a_temp_zface + // all components on z-faces of grid + a1_temp_zface(i, j, k, comp) = (M_normalization != 0) ? -(dt * a_temp_dynamic_coeff * H1_eff[comp] + a1_temp_static_zface(i, j, k, comp)) + : -(dt * a_temp_dynamic_coeff * H1_eff[comp] + 0.5 * a1_temp_static_zface(i, j, k, comp) + + 0.5 * mag_alpha_zface_arr(i,j,k) * 1. / std::sqrt(std::pow(M1_zface(i, j, k, 0), 2._rt) + std::pow(M1_zface(i, j, k, 1), 2._rt) + std::pow(M1_zface(i, j, k, 2), 2._rt)) * M1_old_zface(i, j, k, comp)); + a2_temp_zface(i, j, k, comp) = (M_normalization != 0) ? -(dt * a_temp_dynamic_coeff * H2_eff[comp] + a2_temp_static_zface(i, j, k, comp)) + : -(dt * a_temp_dynamic_coeff * H2_eff[comp] + 0.5 * a2_temp_static_zface(i, j, k, comp) + + 0.5 * mag_alpha_zface_arr(i,j,k) * 1. / std::sqrt(std::pow(M2_zface(i, j, k, 0), 2._rt) + std::pow(M2_zface(i, j, k, 1), 2._rt) + std::pow(M2_zface(i, j, k, 2), 2._rt)) * M2_old_zface(i, j, k, comp)); + } + + for (int comp=0; comp<3; ++comp) { + // update M_zface from a and b using the updateM_field + // all components on z-faces of grid + M1_zface(i, j, k, comp) = MacroscopicProperties::updateM_field(i, j, k, comp, a1_temp_zface, b1_temp_static_zface); + M2_zface(i, j, k, comp) = MacroscopicProperties::updateM_field(i, j, k, comp, a2_temp_zface, b2_temp_static_zface); + } + + // temporary normalized magnitude of M_zface field at the fixed point + // re-investigate the way we do Ms interp, in case we encounter the case where Ms changes across two adjacent cells that you are doing interp + amrex::Real M1_magnitude_normalized = std::sqrt(std::pow(M1_zface(i, j, k, 0), 2._rt) + std::pow(M1_zface(i, j, k, 1), 2._rt) + std::pow(M1_zface(i, j, k, 2), 2._rt)) / mag_Ms_zface_arr(i,j,k); + amrex::Real M2_magnitude_normalized = std::sqrt(std::pow(M2_zface(i, j, k, 0), 2._rt) + std::pow(M2_zface(i, j, k, 1), 2._rt) + std::pow(M2_zface(i, j, k, 2), 2._rt)) / mag_Ms_zface_arr(i,j,k); + + if (M_normalization == 1){ + // saturated case; if |M| has drifted from M_s too much, abort. Otherwise, normalize + // check the normalized error + if (amrex::Math::abs(1. - M1_magnitude_normalized) > mag_normalized_error){ + amrex::Abort("Exceed the normalized error of the M_zface field"); + } + if (amrex::Math::abs(1. - M2_magnitude_normalized) > mag_normalized_error){ + amrex::Abort("Exceed the normalized error of the M_zface field"); + } + // normalize the M_zface field + M1_zface(i, j, k, 0) /= M1_magnitude_normalized; + M1_zface(i, j, k, 1) /= M1_magnitude_normalized; + M1_zface(i, j, k, 2) /= M1_magnitude_normalized; + + M2_zface(i, j, k, 0) /= M2_magnitude_normalized; + M2_zface(i, j, k, 1) /= M2_magnitude_normalized; + M2_zface(i, j, k, 2) /= M2_magnitude_normalized; + } + else if (M_normalization == 0){ + // check the normalized error + if (M1_magnitude_normalized > 1._rt + mag_normalized_error){ + amrex::Abort("Caution: Unsaturated material has M_zface exceeding the saturation magnetization"); + } + else if (M1_magnitude_normalized > 1._rt && M1_magnitude_normalized <= 1._rt + mag_normalized_error){ + // normalize the M_zface field + M1_zface(i, j, k, 0) /= M1_magnitude_normalized; + M1_zface(i, j, k, 1) /= M1_magnitude_normalized; + M1_zface(i, j, k, 2) /= M1_magnitude_normalized; + } + if (M2_magnitude_normalized > 1._rt + mag_normalized_error){ + amrex::Abort("Caution: Unsaturated material has M_zface exceeding the saturation magnetization"); + } + else if (M2_magnitude_normalized > 1._rt && M2_magnitude_normalized <= 1._rt + mag_normalized_error){ + // normalize the M_zface field + M2_zface(i, j, k, 0) /= M2_magnitude_normalized; + M2_zface(i, j, k, 1) /= M2_magnitude_normalized; + M2_zface(i, j, k, 2) /= M2_magnitude_normalized; + } + } + + // calculate M_error_zface + // x,y,z component on z-faces of grid + for (int icomp = 0; icomp < 3; ++icomp) { + M1_error_zface(i, j, k, icomp) = amrex::Math::abs((M1_zface(i, j, k, icomp) - M1_prev_zface(i, j, k, icomp))) / mag_Ms_zface_arr(i,j,k); + M2_error_zface(i, j, k, icomp) = amrex::Math::abs((M2_zface(i, j, k, icomp) - M2_prev_zface(i, j, k, icomp))) / mag_Ms_zface_arr(i,j,k); + } + } + }); + } + + // update H + for (MFIter mfi(*Hfield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi){ + + auto& mag_Ms_xface_mf = macroscopic_properties->getmag_Ms_mf(0); + auto& mag_Ms_yface_mf = macroscopic_properties->getmag_Ms_mf(1); + auto& mag_Ms_zface_mf = macroscopic_properties->getmag_Ms_mf(2); + + // extract material properties + Array4 const& mag_Ms_xface_arr = mag_Ms_xface_mf.array(mfi); + Array4 const& mag_Ms_yface_arr = mag_Ms_yface_mf.array(mfi); + Array4 const& mag_Ms_zface_arr = mag_Ms_zface_mf.array(mfi); + + // Extract field data for this grid/tile + Array4 const &Hx = Hfield[0]->array(mfi); + Array4 const &Hy = Hfield[1]->array(mfi); + Array4 const &Hz = Hfield[2]->array(mfi); + Array4 const &Hx_old = Hfield_old[0]->array(mfi); + Array4 const &Hy_old = Hfield_old[1]->array(mfi); + Array4 const &Hz_old = Hfield_old[2]->array(mfi); + Array4 const &Ex = Efield[0]->array(mfi); + Array4 const &Ey = Efield[1]->array(mfi); + Array4 const &Ez = Efield[2]->array(mfi); + Array4 const &M1_xface = M1field[0]->array(mfi); // note M_xface include x,y,z components at |_x faces + Array4 const &M1_yface = M1field[1]->array(mfi); // note M_yface include x,y,z components at |_y faces + Array4 const &M1_zface = M1field[2]->array(mfi); // note M_zface include x,y,z components at |_z faces + Array4 const &M1_xface_old = M1field_old[0]->array(mfi); // note M_xface_old include x,y,z components at |_x faces + Array4 const &M1_yface_old = M1field_old[1]->array(mfi); // note M_yface_old include x,y,z components at |_y faces + Array4 const &M1_zface_old = M1field_old[2]->array(mfi); // note M_zface_old include x,y,z components at |_z faces + Array4 const &M2_xface = M2field[0]->array(mfi); // note M_xface include x,y,z components at |_x faces + Array4 const &M2_yface = M2field[1]->array(mfi); // note M_yface include x,y,z components at |_y faces + Array4 const &M2_zface = M2field[2]->array(mfi); // note M_zface include x,y,z components at |_z faces + Array4 const &M2_xface_old = M2field_old[0]->array(mfi); // note M_xface_old include x,y,z components at |_x faces + Array4 const &M2_yface_old = M2field_old[1]->array(mfi); // note M_yface_old include x,y,z components at |_y faces + Array4 const &M2_zface_old = M2field_old[2]->array(mfi); // note M_zface_old include x,y,z components at |_z faces + Array4 const &M_xface = Mfield[0]->array(mfi); // note M_xface include x,y,z components at |_x faces + Array4 const &M_yface = Mfield[1]->array(mfi); // note M_yface include x,y,z components at |_y faces + Array4 const &M_zface = Mfield[2]->array(mfi); // note M_zface include x,y,z components at |_z faces + Array4 const &M_xface_old = Mfield_old[0]->array(mfi); // note M_xface_old include x,y,z components at |_x faces + Array4 const &M_yface_old = Mfield_old[1]->array(mfi); // note M_yface_old include x,y,z components at |_y faces + Array4 const &M_zface_old = Mfield_old[2]->array(mfi); // note M_zface_old include x,y,z components at |_z faces + + // Extract stencil coefficients + amrex::Real const *const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr(); + int const n_coefs_x = m_stencil_coefs_x.size(); + amrex::Real const *const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr(); + int const n_coefs_y = m_stencil_coefs_y.size(); + amrex::Real const *const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr(); + int const n_coefs_z = m_stencil_coefs_z.size(); + + // Extract tileboxes for which to loop + amrex::IntVect Hxnodal = Hfield[0]->ixType().toIntVect(); + amrex::IntVect Hynodal = Hfield[1]->ixType().toIntVect(); + amrex::IntVect Hznodal = Hfield[2]->ixType().toIntVect(); + Box const &tbx = mfi.tilebox(Hxnodal); + Box const &tby = mfi.tilebox(Hynodal); + Box const &tbz = mfi.tilebox(Hznodal); + + amrex::Array4 const& mu_arr = mu_mf.array(mfi); + + amrex::Real const mu0_inv = 1. / PhysConst::mu0; + + // Loop over the cells and update the fields + amrex::ParallelFor(tbx, tby, tbz, + + [=] AMREX_GPU_DEVICE(int i, int j, int k) { + + + for(int n=0; n<3; ++n){ + + M_xface(i,j,k,n) = (M1_xface(i,j,k,n) + M2_xface(i,j,k,n))/2.0 ; + + } + + if (mag_Ms_xface_arr(i,j,k) == 0._rt){ // nonmagnetic region + amrex::Real mu_arrx = ablastr::coarsen::sample::Interp( mu_arr, mu_stag, Hx_stag, macro_cr, i, j, k, 0); + Hx(i, j, k) = Hx_old(i, j, k) + 1. / mu_arrx * dt * (T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k) + - T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k)); + } else if (mag_Ms_xface_arr(i,j,k) > 0){ // magnetic region + Hx(i, j, k) = Hx_old(i, j, k) + mu0_inv * dt * (T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k) + - T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k)); + if (coupling == 1) { + Hx(i, j, k) += - M_xface(i, j, k, 0) + M_xface_old(i, j, k, 0); + } + } + }, + + [=] AMREX_GPU_DEVICE(int i, int j, int k) { + + for(int n=0; n<3; ++n){ + M_yface(i,j,k,n) = (M1_yface(i,j,k,n) + M2_yface(i,j,k,n))/2.0; + } + + if (mag_Ms_yface_arr(i,j,k) == 0._rt){ // nonmagnetic region + amrex::Real mu_arry = ablastr::coarsen::sample::Interp( mu_arr, mu_stag, Hy_stag, macro_cr, i, j, k, 0); + Hy(i, j, k) = Hy_old(i, j, k) + 1. / mu_arry * dt * (T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k) + - T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k)); + } else if (mag_Ms_yface_arr(i,j,k) > 0){ // magnetic region + Hy(i, j, k) = Hy_old(i, j, k) + mu0_inv * dt * (T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k) + - T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k)); + if (coupling == 1){ + Hy(i, j, k) += - M_yface(i, j, k, 1) + M_yface_old(i, j, k, 1); + } + } + }, + + [=] AMREX_GPU_DEVICE(int i, int j, int k) { + + for(int n=0; n<3; ++n){ + M_zface(i,j,k,n) = (M1_zface(i,j,k,n) + M2_zface(i,j,k,n))/2.0; + } + + if (mag_Ms_zface_arr(i,j,k) == 0._rt){ // nonmagnetic region + amrex::Real mu_arrz = ablastr::coarsen::sample::Interp( mu_arr, mu_stag, Hz_stag, macro_cr, i, j, k, 0); + Hz(i, j, k) = Hz_old(i, j, k) + 1. / mu_arrz * dt * (T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k) + - T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k)); + } else if (mag_Ms_zface_arr(i,j,k) > 0){ // magnetic region + Hz(i, j, k) = Hz_old(i, j, k) + mu0_inv * dt * (T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k) + - T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k)); + if (coupling == 1){ + Hz(i, j, k) += - M_zface(i, j, k, 2) + M_zface_old(i, j, k, 2); + } + } + } + + ); + } + + // Check the error between Mfield and Mfield_prev and decide whether another iteration is needed + amrex::Real M_iter_maxerror = -1._rt; + for (int iface = 0; iface < 3; iface++){ + for (int jcomp = 0; jcomp < 3; jcomp++){ + Real M1_iter_error = M1field_error[iface]->norm0(jcomp); + Real M2_iter_error = M2field_error[iface]->norm0(jcomp); + if (M1_iter_error >= M_iter_maxerror || M2_iter_error >= M_iter_maxerror){ + if(M2_iter_error >= M1_iter_error){ + M_iter_maxerror = M2_iter_error; + } + else{ + M_iter_maxerror = M1_iter_error; + } + } + } + } + + if (M_iter_maxerror <= M_tol){ + + stop_iter = 1; + + // normalize M + if (M_normalization == 2){ + + for (MFIter mfi(*Mfield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi){ + + auto& mag_Ms_xface_mf = macroscopic_properties->getmag_Ms_mf(0); + auto& mag_Ms_yface_mf = macroscopic_properties->getmag_Ms_mf(1); + auto& mag_Ms_zface_mf = macroscopic_properties->getmag_Ms_mf(2); + + // extract material properties + Array4 const& mag_Ms_xface_arr = mag_Ms_xface_mf.array(mfi); + Array4 const& mag_Ms_yface_arr = mag_Ms_yface_mf.array(mfi); + Array4 const& mag_Ms_zface_arr = mag_Ms_zface_mf.array(mfi); + + // extract field data + Array4 const &M1_xface = M1field[0]->array(mfi); // note M_xface include x,y,z components at |_x faces + Array4 const &M1_yface = M1field[1]->array(mfi); // note M_yface include x,y,z components at |_y faces + Array4 const &M1_zface = M1field[2]->array(mfi); // note M_zface include x,y,z components at |_z faces + + Array4 const &M2_xface = M2field[0]->array(mfi); // note M_xface include x,y,z components at |_x faces + Array4 const &M2_yface = M2field[1]->array(mfi); // note M_yface include x,y,z components at |_y faces + Array4 const &M2_zface = M2field[2]->array(mfi); // note M_zface include x,y,z components at |_z faces + + // extract tileboxes for which to loop + amrex::IntVect Mxface_stag = Mfield[0]->ixType().toIntVect(); + amrex::IntVect Myface_stag = Mfield[1]->ixType().toIntVect(); + amrex::IntVect Mzface_stag = Mfield[2]->ixType().toIntVect(); + Box const &tbx = mfi.tilebox(Mxface_stag); /* just define which grid type */ + Box const &tby = mfi.tilebox(Myface_stag); + Box const &tbz = mfi.tilebox(Mzface_stag); + + // loop over cells and update fields + amrex::ParallelFor(tbx, tby, tbz, + [=] AMREX_GPU_DEVICE(int i, int j, int k) { + + if (mag_Ms_xface_arr(i,j,k) > 0._rt){ + // temporary normalized magnitude of M_xface field at the fixed point + // re-investigate the way we do Ms interp, in case we encounter the case where Ms changes across two adjacent cells that you are doing interp + amrex::Real M1_magnitude_normalized = std::sqrt(std::pow(M1_xface(i, j, k, 0), 2._rt) + std::pow(M1_xface(i, j, k, 1), 2._rt) + + std::pow(M1_xface(i, j, k, 2), 2._rt)) / + mag_Ms_xface_arr(i,j,k); + amrex::Real M2_magnitude_normalized = std::sqrt(std::pow(M2_xface(i, j, k, 0), 2._rt) + std::pow(M2_xface(i, j, k, 1), 2._rt) + + std::pow(M2_xface(i, j, k, 2), 2._rt)) / + mag_Ms_xface_arr(i,j,k); + + // check the normalized error + if (amrex::Math::abs(1._rt - M1_magnitude_normalized) > mag_normalized_error){ + amrex::Abort("Exceed the normalized error of the M_xface field"); + } + if (amrex::Math::abs(1._rt - M2_magnitude_normalized) > mag_normalized_error){ + amrex::Abort("Exceed the normalized error of the M_xface field"); + } + // normalize the M_xface field + M1_xface(i, j, k, 0) /= M1_magnitude_normalized; + M1_xface(i, j, k, 1) /= M1_magnitude_normalized; + M1_xface(i, j, k, 2) /= M1_magnitude_normalized; + + M2_xface(i, j, k, 0) /= M2_magnitude_normalized; + M2_xface(i, j, k, 1) /= M2_magnitude_normalized; + M2_xface(i, j, k, 2) /= M2_magnitude_normalized; + } + }, + + [=] AMREX_GPU_DEVICE(int i, int j, int k) { + + if (mag_Ms_yface_arr(i,j,k) > 0._rt){ + // temporary normalized magnitude of M_yface field at the fixed point + // re-investigate the way we do Ms interp, in case we encounter the case where Ms changes across two adjacent cells that you are doing interp + amrex::Real M1_magnitude_normalized = std::sqrt(std::pow(M1_yface(i, j, k, 0), 2._rt) + std::pow(M1_yface(i, j, k, 1), 2._rt) + + std::pow(M1_yface(i, j, k, 2), 2._rt)) / + mag_Ms_yface_arr(i,j,k); + amrex::Real M2_magnitude_normalized = std::sqrt(std::pow(M2_yface(i, j, k, 0), 2._rt) + std::pow(M2_yface(i, j, k, 1), 2._rt) + + std::pow(M2_yface(i, j, k, 2), 2._rt)) / + mag_Ms_yface_arr(i,j,k); + + // check the normalized error + if (amrex::Math::abs(1._rt - M1_magnitude_normalized) > mag_normalized_error){ + amrex::Abort("Exceed the normalized error of the M_yface field"); + } + if (amrex::Math::abs(1._rt - M2_magnitude_normalized) > mag_normalized_error){ + amrex::Abort("Exceed the normalized error of the M_yface field"); + } + // normalize the M_yface field + M1_yface(i, j, k, 0) /= M1_magnitude_normalized; + M1_yface(i, j, k, 1) /= M1_magnitude_normalized; + M1_yface(i, j, k, 2) /= M1_magnitude_normalized; + + M2_yface(i, j, k, 0) /= M2_magnitude_normalized; + M2_yface(i, j, k, 1) /= M2_magnitude_normalized; + M2_yface(i, j, k, 2) /= M2_magnitude_normalized; + } + }, + + [=] AMREX_GPU_DEVICE(int i, int j, int k) { + + if (mag_Ms_zface_arr(i,j,k) > 0._rt){ + // temporary normalized magnitude of M_zface field at the fixed point + // re-investigate the way we do Ms interp, in case we encounter the case where Ms changes across two adjacent cells that you are doing interp + amrex::Real M1_magnitude_normalized = std::sqrt(std::pow(M1_zface(i, j, k, 0), 2._rt) + std::pow(M1_zface(i, j, k, 1), 2._rt) + + std::pow(M1_zface(i, j, k, 2), 2._rt)) / + mag_Ms_zface_arr(i,j,k); + amrex::Real M2_magnitude_normalized = std::sqrt(std::pow(M2_zface(i, j, k, 0), 2._rt) + std::pow(M2_zface(i, j, k, 1), 2._rt) + + std::pow(M2_zface(i, j, k, 2), 2._rt)) / + mag_Ms_zface_arr(i,j,k); + + // check the normalized error + if (amrex::Math::abs(1. - M1_magnitude_normalized) > mag_normalized_error){ + amrex::Abort("Exceed the normalized error of the M_zface field"); + } + if (amrex::Math::abs(1. - M2_magnitude_normalized) > mag_normalized_error){ + amrex::Abort("Exceed the normalized error of the M_zface field"); + } + // normalize the M_zface field + M1_zface(i, j, k, 0) /= M1_magnitude_normalized; + M1_zface(i, j, k, 1) /= M1_magnitude_normalized; + M1_zface(i, j, k, 2) /= M1_magnitude_normalized; + + M2_zface(i, j, k, 0) /= M2_magnitude_normalized; + M2_zface(i, j, k, 1) /= M2_magnitude_normalized; + M2_zface(i, j, k, 2) /= M2_magnitude_normalized; + } + }); + } + } + } + else{ + + const auto& period = warpx.Geom(lev).periodicity(); + // Copy Mfield to Mfield_previous and fill periodic/interior ghost cells + for (int i = 0; i < 3; i++){ + MultiFab::Copy(*Mfield_prev[i], *Mfield[i], 0, 0, 3, Mfield[i]->nGrow()); + (*Mfield_prev[i]).FillBoundary(Mfield[i]->nGrowVect(), period); + + MultiFab::Copy(*M1field_prev[i], *M1field[i], 0, 0, 3, M1field[i]->nGrow()); + (*M1field_prev[i]).FillBoundary(M1field[i]->nGrowVect(), period); + + MultiFab::Copy(*M2field_prev[i], *M2field[i], 0, 0, 3, M2field[i]->nGrow()); + (*M2field_prev[i]).FillBoundary(M2field[i]->nGrowVect(), period); + } + } + + if (M_iter >= M_max_iter){ + amrex::Abort("The M_iter exceeds the M_max_iter"); + amrex::Print() << "The M_iter = " << M_iter << " exceeds the M_max_iter = " << M_max_iter << std::endl; + } + else{ + M_iter++; + amrex::Print() << "Finish " << M_iter << " times iteration with M_iter_maxerror = " << M_iter_maxerror << " and M_tol = " << M_tol << std::endl; + } + + } // end the iteration + + // update B + for (MFIter mfi(*Bfield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi){ + + auto& mag_Ms_xface_mf = macroscopic_properties->getmag_Ms_mf(0); + auto& mag_Ms_yface_mf = macroscopic_properties->getmag_Ms_mf(1); + auto& mag_Ms_zface_mf = macroscopic_properties->getmag_Ms_mf(2); + + // Extract field data for this grid/tile + Array4 const& mag_Ms_xface_arr = mag_Ms_xface_mf.array(mfi); + Array4 const& mag_Ms_yface_arr = mag_Ms_yface_mf.array(mfi); + Array4 const& mag_Ms_zface_arr = mag_Ms_zface_mf.array(mfi); + Array4 const &Hx = Hfield[0]->array(mfi); + Array4 const &Hy = Hfield[1]->array(mfi); + Array4 const &Hz = Hfield[2]->array(mfi); + Array4 const &Bx = Bfield[0]->array(mfi); + Array4 const &By = Bfield[1]->array(mfi); + Array4 const &Bz = Bfield[2]->array(mfi); + Array4 const &L_xface = Lfield[0]->array(mfi); + Array4 const &L_yface = Lfield[1]->array(mfi); + Array4 const &L_zface = Lfield[2]->array(mfi); + Array4 const &M_xface = Mfield[0]->array(mfi); // note M_xface include x,y,z components at |_x faces + Array4 const &M_yface = Mfield[1]->array(mfi); // note M_yface include x,y,z components at |_y faces + Array4 const &M_zface = Mfield[2]->array(mfi); // note M_zface include x,y,z components at |_z faces + Array4 const &M1_xface = M1field[0]->array(mfi); // note M_xface include x,y,z components at |_x faces + Array4 const &M1_yface = M1field[1]->array(mfi); // note M_yface include x,y,z components at |_y faces + Array4 const &M1_zface = M1field[2]->array(mfi); // note M_zface include x,y,z components at |_z faces + Array4 const &M2_xface = M2field[0]->array(mfi); // note M_xface include x,y,z components at |_x faces + Array4 const &M2_yface = M2field[1]->array(mfi); // note M_yface include x,y,z components at |_y faces + Array4 const &M2_zface = M2field[2]->array(mfi); // note M_zface include x,y,z components at |_z faces + + // Extract tileboxes for which to loop + amrex::IntVect Bxnodal = Bfield[0]->ixType().toIntVect(); + amrex::IntVect Bynodal = Bfield[1]->ixType().toIntVect(); + amrex::IntVect Bznodal = Bfield[2]->ixType().toIntVect(); + Box const &tbx = mfi.tilebox(Bxnodal); + Box const &tby = mfi.tilebox(Bynodal); + Box const &tbz = mfi.tilebox(Bznodal); + + amrex::Array4 const& mu_arr = mu_mf.array(mfi); + + // Loop over the cells and update the fields + amrex::ParallelFor(tbx, tby, tbz, + + [=] AMREX_GPU_DEVICE(int i, int j, int k) { + + if (mag_Ms_xface_arr(i,j,k) == 0._rt){ // nonmagnetic region + amrex::Real mu_arrx = ablastr::coarsen::sample::Interp( mu_arr, mu_stag, Bx_stag, macro_cr, i, j, k, 0); + Bx(i, j, k) = mu_arrx * Hx(i, j, k); + } else if (mag_Ms_xface_arr(i,j,k) > 0){ + Bx(i, j, k) = PhysConst::mu0 * (M_xface(i, j, k, 0) + Hx(i, j, k)); + } + + for(int n=0; n<3; ++n){ + L_xface(i,j,k,n) = (M1_xface(i,j,k,n) - M2_xface(i,j,k,n))/2.0; + } + + }, + + [=] AMREX_GPU_DEVICE(int i, int j, int k) { + + if (mag_Ms_yface_arr(i,j,k) == 0._rt){ // nonmagnetic region + amrex::Real mu_arry = ablastr::coarsen::sample::Interp( mu_arr, mu_stag, By_stag, macro_cr, i, j, k, 0); + By(i, j, k) = mu_arry * Hy(i, j, k); + } else if (mag_Ms_yface_arr(i,j,k) > 0){ + By(i, j, k) = PhysConst::mu0 * (M_yface(i, j, k, 1) + Hy(i, j, k)); + } + + for(int n=0; n<3; ++n){ + L_yface(i,j,k,n) = (M1_yface(i,j,k,n) - M2_yface(i,j,k,n))/2.0; + } + }, + + [=] AMREX_GPU_DEVICE(int i, int j, int k) { + + if (mag_Ms_zface_arr(i,j,k) == 0._rt){ // nonmagnetic region + amrex::Real mu_arrz = ablastr::coarsen::sample::Interp( mu_arr, mu_stag, Bz_stag, macro_cr, i, j, k, 0); + Bz(i, j, k) = mu_arrz * Hz(i, j, k); + } else if (mag_Ms_zface_arr(i,j,k) > 0){ + Bz(i, j, k) = PhysConst::mu0 * (M_zface(i, j, k, 2) + Hz(i, j, k)); + } + + for(int n=0; n<3; ++n){ + L_zface(i,j,k,n) = (M1_zface(i,j,k,n) - M2_zface(i,j,k,n))/2; + } + } + + ); + } + +} +#endif // ifdef WARPX_MAG_LLG +#endif // ifndef WARPX_DIM_RZ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H index c8784aafd..6dd33d5e3 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H @@ -183,6 +183,8 @@ public: amrex::MultiFab * getmag_pointer_exchange (int dir) {return m_mag_exchange_mf[dir].get();} amrex::MultiFab& getmag_anisotropy_mf(int dir) {return (*m_mag_anisotropy_mf[dir]);} amrex::MultiFab * getmag_pointer_anisotropy (int dir) {return m_mag_anisotropy_mf[dir].get();} + amrex::MultiFab& getmag_afm_mf(int dir) {return (*m_mag_afm_mf[dir]);} + amrex::MultiFab * getmag_pointer_afm (int dir) {return m_mag_afm_mf[dir].get();} amrex::Real getmag_normalized_error () {return m_mag_normalized_error;} int getmag_max_iter () {return m_mag_max_iter;} @@ -312,6 +314,8 @@ public: amrex::Real m_mag_exchange; /** The coefficient of the anisotropy coupling term, only applies for magnetic materials */ amrex::Real m_mag_anisotropy; + /** The coefficient of the antiferromagnetic coupling term, only applies for magnetic materials */ + amrex::Real m_mag_afm; // If the magnitude of the magnetization deviates by more than this amount relative // to the user-defined Ms, abort. Default 0.1. @@ -333,6 +337,8 @@ public: std::array, 3> m_mag_exchange_mf; /** Multifabs storing spatially varying coefficient of the anisotropy coupling term on three faces */ std::array, 3> m_mag_anisotropy_mf; + /** Multifabs storing spatially varying coefficient of the antiferromagnetic coupling term on three faces */ + std::array, 3> m_mag_afm_mf; // these store the type of initialization, e.g., "constant", "parse_X_function", etc. std::string m_mag_Ms_s; @@ -340,6 +346,7 @@ public: std::string m_mag_gamma_s; std::string m_mag_exchange_s; std::string m_mag_anisotropy_s; + std::string m_mag_afm_s; // these store the parsed expression for the material properties std::string m_str_mag_Ms_function; @@ -347,12 +354,14 @@ public: std::string m_str_mag_gamma_function; std::string m_str_mag_exchange_function; std::string m_str_mag_anisotropy_function; + std::string m_str_mag_afm_function; std::unique_ptr m_mag_Ms_parser; std::unique_ptr m_mag_alpha_parser; std::unique_ptr m_mag_gamma_parser; std::unique_ptr m_mag_exchange_parser; std::unique_ptr m_mag_anisotropy_parser; + std::unique_ptr m_mag_afm_parser; #endif private: diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp index b2388c12a..5f21a91b8 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp @@ -250,6 +250,18 @@ MacroscopicProperties::ReadParameters () } } + if (warpx.mag_LLG_afm_coupling == 1) { // magnetic crystal is considered to have one sublattice by default + pp_macroscopic.get("mag_afm_init_style", m_mag_afm_s); + if (m_mag_afm_s == "constant") pp_macroscopic.get("mag_afm", m_mag_afm); + // _mag_ such that it's clear the exch variable is only meaningful for magnetic materials + //initialization with parser + if (m_mag_afm_s == "parse_mag_afm_function") { + utils::parser::Store_parserString(pp_macroscopic, "mag_afm_function(x,y,z)", m_str_mag_afm_function); + m_mag_afm_parser = std::make_unique( + utils::parser::makeParser(m_str_mag_afm_function,{"x","y","z"})); + } + } + m_mag_normalized_error = 0.1; pp_macroscopic.query("mag_normalized_error",m_mag_normalized_error); @@ -425,6 +437,7 @@ MacroscopicProperties::InitData () m_mag_gamma_mf[i] = std::make_unique(amrex::convert(ba,IntVect::TheDimensionVector(i)), dmap, 1, ng_EB_alloc); m_mag_exchange_mf[i] = std::make_unique(amrex::convert(ba,IntVect::TheDimensionVector(i)), dmap, 1, ng_EB_alloc); m_mag_anisotropy_mf[i] = std::make_unique(amrex::convert(ba,IntVect::TheDimensionVector(i)), dmap, 1, ng_EB_alloc); + m_mag_afm_mf[i] = std::make_unique(amrex::convert(ba,IntVect::TheDimensionVector(i)), dmap, 1, ng_EB_alloc); } // mag_Ms - defined at cell centers @@ -509,6 +522,18 @@ MacroscopicProperties::InitData () InitializeMacroMultiFabUsingParser(m_mag_anisotropy_mf[1].get(), m_mag_anisotropy_parser->compile<3>(), lev); InitializeMacroMultiFabUsingParser(m_mag_anisotropy_mf[2].get(), m_mag_anisotropy_parser->compile<3>(), lev); } + + // mag_afm - defined at faces + if (m_mag_afm_s == "constant") { + m_mag_afm_mf[0]->setVal(m_mag_afm); + m_mag_afm_mf[1]->setVal(m_mag_afm); + m_mag_afm_mf[2]->setVal(m_mag_afm); + } + else if (m_mag_afm_s == "parse_mag_afm_function"){ + InitializeMacroMultiFabUsingParser(m_mag_afm_mf[0].get(), m_mag_afm_parser->compile<3>(), lev); + InitializeMacroMultiFabUsingParser(m_mag_afm_mf[1].get(), m_mag_afm_parser->compile<3>(), lev); + InitializeMacroMultiFabUsingParser(m_mag_afm_mf[2].get(), m_mag_afm_parser->compile<3>(), lev); + } #endif @@ -788,4 +813,4 @@ MacroscopicProperties::InitializePECFromSigma (amrex::MultiFab* sigma_mf, }); } -} +} \ No newline at end of file diff --git a/Source/FieldSolver/FiniteDifferenceSolver/Make.package b/Source/FieldSolver/FiniteDifferenceSolver/Make.package index 0f96e6bb5..0fcef8471 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/Make.package +++ b/Source/FieldSolver/FiniteDifferenceSolver/Make.package @@ -10,6 +10,7 @@ CEXE_sources += MacroscopicEvolveE.cpp #ifdef WARPX_MAG_LLG CEXE_sources += MacroscopicEvolveHM.cpp CEXE_sources += MacroscopicEvolveHM_2nd.cpp +CEXE_sources += MacroscopicEvolveHM_AFM_2nd.cpp CEXE_sources += EvolveHPML.cpp #endif diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index aa7e6a768..50e414880 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -1233,6 +1233,49 @@ WarpX::MacroscopicEvolveHM_2nd (int lev, PatchType patch_type, amrex::Real a_dt) } } +// define WarpX::MacroscopicEvolveHM_AFM_2nd +void +WarpX::MacroscopicEvolveHM_AFM_2nd (amrex::Real a_dt) +{ + for (int lev = 0; lev <= finest_level; ++lev ) { + MacroscopicEvolveHM_AFM_2nd(lev, a_dt); + } +} + +void +WarpX::MacroscopicEvolveHM_AFM_2nd (int lev, amrex::Real a_dt) { + + WARPX_PROFILE("WarpX::MacroscopicEvolveHM_AFM_2nd()"); + MacroscopicEvolveHM_AFM_2nd(lev, PatchType::fine, a_dt); + if (lev > 0) { + amrex::Abort("Macroscopic EvolveHM_AFM_2nd is not implemented for lev>0, yet."); + } +} + +void +WarpX::MacroscopicEvolveHM_AFM_2nd (int lev, PatchType patch_type, amrex::Real a_dt) { + + // Evolve H field in regular cells + if (patch_type == PatchType::fine) { + m_fdtd_solver_fp[lev]->MacroscopicEvolveHM_AFM_2nd(lev, M1field_fp[lev], M2field_fp[lev], Mfield_fp[lev], Hfield_fp[lev], H1field_fp[lev], H2field_fp[lev], Bfield_fp[lev], Lfield_fp[lev], H_biasfield_fp[lev], Efield_fp[lev], + a_dt, m_macroscopic_properties); + } + else { + amrex::Abort("Macroscopic EvolveHM_AFM_2nd is not implemented for lev > 0 yet"); + } + + // Evolve H field in PML cells + if (do_pml && pml[lev]->ok()) { + if (patch_type == PatchType::fine) { + m_fdtd_solver_fp[lev]->EvolveHPML( + pml[lev]->GetH_fp(), pml[lev]->GetE_fp(), a_dt, WarpX::do_dive_cleaning ); + } else { + m_fdtd_solver_cp[lev]->EvolveHPML( + pml[lev]->GetH_cp(), pml[lev]->GetE_cp(), a_dt, WarpX::do_dive_cleaning ); + } + } +} + #endif #endif // ifndef WARPX_DIM_RZ diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index ade792947..48033b63a 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -1295,12 +1295,283 @@ WarpX::InitLevelData (int lev, Real /*time*/) 'H', lev, PatchType::coarse); } + if (mag_LLG_afm_coupling==1){ + H1xfield_parser = std::make_unique( + utils::parser::makeParser(str_Hx_ext_grid_function,{"x","y","z"})); + H1yfield_parser = std::make_unique( + utils::parser::makeParser(str_Hy_ext_grid_function,{"x","y","z"})); + H1zfield_parser = std::make_unique( + utils::parser::makeParser(str_Hz_ext_grid_function,{"x","y","z"})); + + // Initialize H1field_fp with external function + InitializeExternalFieldsOnGridUsingParser(H1field_fp[lev][0].get(), + H1field_fp[lev][1].get(), + H1field_fp[lev][2].get(), + H1xfield_parser->compile<3>(), + H1yfield_parser->compile<3>(), + H1zfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'S', //H1 field + lev, PatchType::fine); + if (lev > 0) { + InitializeExternalFieldsOnGridUsingParser(H1field_aux[lev][0].get(), + H1field_aux[lev][1].get(), + H1field_aux[lev][2].get(), + H1xfield_parser->compile<3>(), + H1yfield_parser->compile<3>(), + H1zfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'S', //H1 field + lev, PatchType::coarse); + + InitializeExternalFieldsOnGridUsingParser(H1field_cp[lev][0].get(), + H1field_cp[lev][1].get(), + H1field_cp[lev][2].get(), + H1xfield_parser->compile<3>(), + H1yfield_parser->compile<3>(), + H1zfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'S', //H1 field + lev, PatchType::coarse); + } + H2xfield_parser = std::make_unique( + utils::parser::makeParser(str_Hx_ext_grid_function,{"x","y","z"})); + H2yfield_parser = std::make_unique( + utils::parser::makeParser(str_Hy_ext_grid_function,{"x","y","z"})); + H2zfield_parser = std::make_unique( + utils::parser::makeParser(str_Hz_ext_grid_function,{"x","y","z"})); + + // Initialize H2field_fp with external function + InitializeExternalFieldsOnGridUsingParser(H2field_fp[lev][0].get(), + H2field_fp[lev][1].get(), + H2field_fp[lev][2].get(), + H2xfield_parser->compile<3>(), + H2yfield_parser->compile<3>(), + H2zfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'C', //H2 field + lev, PatchType::fine); + if (lev > 0) { + InitializeExternalFieldsOnGridUsingParser(H2field_aux[lev][0].get(), + H2field_aux[lev][1].get(), + H2field_aux[lev][2].get(), + H2xfield_parser->compile<3>(), + H2yfield_parser->compile<3>(), + H2zfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'C', //H2 field + lev, PatchType::coarse); + + InitializeExternalFieldsOnGridUsingParser(H2field_cp[lev][0].get(), + H2field_cp[lev][1].get(), + H2field_cp[lev][2].get(), + H2xfield_parser->compile<3>(), + H2yfield_parser->compile<3>(), + H2zfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'C', //H2 field + lev, PatchType::coarse); + } + } } if (M_ext_grid_s == "parse_m_ext_grid_function") { #ifdef WARPX_DIM_RZ amrex::Abort("M-field parser for external fields does not work with RZ"); #endif + if(mag_LLG_afm_coupling == 1){ + utils::parser::Store_parserString(pp_warpx, "Mx_external_grid_function(x,y,z)", + str_Mx_ext_grid_function); + utils::parser::Store_parserString(pp_warpx, "My_external_grid_function(x,y,z)", + str_My_ext_grid_function); + utils::parser::Store_parserString(pp_warpx, "Mz_external_grid_function(x,y,z)", + str_Mz_ext_grid_function); + + M1xfield_parser = std::make_unique( + utils::parser::makeParser(str_Mx_ext_grid_function,{"x","y","z"})); + M1yfield_parser = std::make_unique( + utils::parser::makeParser(str_My_ext_grid_function,{"x","y","z"})); + M1zfield_parser = std::make_unique( + utils::parser::makeParser(str_Mz_ext_grid_function,{"x","y","z"})); + + // Initialize M1field_fp with external function directly on the faces + +amrex::Print() << "M1field_fp.size()=" << M1field_fp.size() << " lev=" << lev << "\n"; +if (static_cast(M1field_fp.size()) <= lev || !M1field_fp[lev][0]) { + amrex::Abort("M1field_fp not allocated at this level"); +} + + InitializeExternalFieldsOnGridUsingParser(M1field_fp[lev][0].get(), + M1field_fp[lev][1].get(), + M1field_fp[lev][2].get(), + M1xfield_parser->compile<3>(), + M1yfield_parser->compile<3>(), + M1zfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'U', //M1 field + lev, PatchType::fine); + if (lev > 0) { + InitializeExternalFieldsOnGridUsingParser(M1field_aux[lev][0].get(), + M1field_aux[lev][1].get(), + M1field_aux[lev][2].get(), + M1xfield_parser->compile<3>(), + M1yfield_parser->compile<3>(), + M1zfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'U', //M1 field + lev, PatchType::coarse); + + InitializeExternalFieldsOnGridUsingParser(M1field_cp[lev][0].get(), + M1field_cp[lev][1].get(), + M1field_cp[lev][2].get(), + M1xfield_parser->compile<3>(), + M1yfield_parser->compile<3>(), + M1zfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'U', //M1 field + lev, PatchType::coarse); + } + // Construct M2 parsers using the negated string from M1 + M2xfield_parser = std::make_unique( + utils::parser::makeParser("-(" + str_Mx_ext_grid_function + ")", {"x","y","z"})); + M2yfield_parser = std::make_unique( + utils::parser::makeParser("-(" + str_My_ext_grid_function + ")", {"x","y","z"})); + M2zfield_parser = std::make_unique( + utils::parser::makeParser("-(" + str_Mz_ext_grid_function + ")", {"x","y","z"})); + + // Initialize M2field_fp with the negated M1 function directly on the faces + InitializeExternalFieldsOnGridUsingParser(M2field_fp[lev][0].get(), + M2field_fp[lev][1].get(), + M2field_fp[lev][2].get(), + M2xfield_parser->compile<3>(), + M2yfield_parser->compile<3>(), + M2zfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'D', //M2 field + lev, PatchType::fine); + if (lev > 0) { + InitializeExternalFieldsOnGridUsingParser(M2field_aux[lev][0].get(), + M2field_aux[lev][1].get(), + M2field_aux[lev][2].get(), + M2xfield_parser->compile<3>(), + M2yfield_parser->compile<3>(), + M2zfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'D', //M2 field + lev, PatchType::coarse); + + InitializeExternalFieldsOnGridUsingParser(M2field_cp[lev][0].get(), + M2field_cp[lev][1].get(), + M2field_cp[lev][2].get(), + M2xfield_parser->compile<3>(), + M2yfield_parser->compile<3>(), + M2zfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'D', //M2 field + lev, PatchType::coarse); + } + // Construct M parsers by adding the M1 and M2 string expressions + // (Note: M2 is strictly -M1, this evaluates to 0) + Mxfield_parser = std::make_unique( + utils::parser::makeParser("0.0", {"x","y","z"})); + Myfield_parser = std::make_unique( + utils::parser::makeParser("0.0", {"x","y","z"})); + Mzfield_parser = std::make_unique( + utils::parser::makeParser("0.0", {"x","y","z"})); + + // Initialize Mfield_fp directly on the faces + InitializeExternalFieldsOnGridUsingParser(Mfield_fp[lev][0].get(), + Mfield_fp[lev][1].get(), + Mfield_fp[lev][2].get(), + Mxfield_parser->compile<3>(), + Myfield_parser->compile<3>(), + Mzfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'M', + lev, PatchType::fine); + if (lev > 0) { + InitializeExternalFieldsOnGridUsingParser(Mfield_aux[lev][0].get(), + Mfield_aux[lev][1].get(), + Mfield_aux[lev][2].get(), + Mxfield_parser->compile<3>(), + Myfield_parser->compile<3>(), + Mzfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'M', + lev, PatchType::coarse); + + InitializeExternalFieldsOnGridUsingParser(Mfield_cp[lev][0].get(), + Mfield_cp[lev][1].get(), + Mfield_cp[lev][2].get(), + Mxfield_parser->compile<3>(), + Myfield_parser->compile<3>(), + Mzfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'M', + lev, PatchType::coarse); + } + // Construct L parsers by subtracting the M1 and M2 string expressions + // (Note: M2 is strictly -M1, this evaluates to 2*M1) + Lxfield_parser = std::make_unique( + utils::parser::makeParser("(" + str_Mx_ext_grid_function + ")", {"x","y","z"})); + Lyfield_parser = std::make_unique( + utils::parser::makeParser("(" + str_My_ext_grid_function + ")", {"x","y","z"})); + Lzfield_parser = std::make_unique( + utils::parser::makeParser("(" + str_Mz_ext_grid_function + ")", {"x","y","z"})); + + // Initialize Mfield_fp directly on the faces + InitializeExternalFieldsOnGridUsingParser(Lfield_fp[lev][0].get(), + Lfield_fp[lev][1].get(), + Lfield_fp[lev][2].get(), + Lxfield_parser->compile<3>(), + Lyfield_parser->compile<3>(), + Lzfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'L', + lev, PatchType::fine); + if (lev > 0) { + InitializeExternalFieldsOnGridUsingParser(Lfield_aux[lev][0].get(), + Lfield_aux[lev][1].get(), + Lfield_aux[lev][2].get(), + Lxfield_parser->compile<3>(), + Lyfield_parser->compile<3>(), + Lzfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'L', + lev, PatchType::coarse); + + InitializeExternalFieldsOnGridUsingParser(Lfield_cp[lev][0].get(), + Lfield_cp[lev][1].get(), + Lfield_cp[lev][2].get(), + Lxfield_parser->compile<3>(), + Lyfield_parser->compile<3>(), + Lzfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'L', + lev, PatchType::coarse); + } + } + + + else{ utils::parser::Store_parserString(pp_warpx, "Mx_external_grid_function(x,y,z)", str_Mx_ext_grid_function); utils::parser::Store_parserString(pp_warpx, "My_external_grid_function(x,y,z)", @@ -1350,6 +1621,7 @@ WarpX::InitLevelData (int lev, Real /*time*/) lev, PatchType::coarse); } } + } #endif //closes #ifdef WARPX_MAG_LLG @@ -2031,4 +2303,4 @@ WarpX::ReadExternalFieldFromFile (std::string , amrex::MultiFab* ,std::string, s WARPX_ABORT_WITH_MESSAGE("OpenPMD field reading requires OpenPMD support to be enabled"); #endif } -#endif // WARPX_USE_OPENPMD +#endif // WARPX_USE_OPENPMD \ No newline at end of file diff --git a/Source/WarpX.H b/Source/WarpX.H index 1cfe69217..7a910f42a 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -125,6 +125,7 @@ public: static amrex::Vector H_external_grid; static amrex::Vector M_external_grid; static amrex::Vector H_bias_external_grid; + static amrex::Vector M1_external_grid; #endif //! Initialization type for external magnetic field on the grid @@ -143,6 +144,8 @@ public: #ifdef WARPX_MAG_LLG static std::string H_ext_grid_s; // added for H; static std::string M_ext_grid_s; // added for M; + static std::string M1_ext_grid_s; + static std::string M2_ext_grid_s; static std::string H_bias_ext_grid_s; #endif static int ApplyExcitationInPML; @@ -265,10 +268,25 @@ public: std::unique_ptr Hxfield_parser; std::unique_ptr Hyfield_parser; std::unique_ptr Hzfield_parser; + std::unique_ptr H1xfield_parser; + std::unique_ptr H1yfield_parser; + std::unique_ptr H1zfield_parser; + std::unique_ptr H2xfield_parser; + std::unique_ptr H2yfield_parser; + std::unique_ptr H2zfield_parser; // Parser for M_external on the grid std::unique_ptr Mxfield_parser; std::unique_ptr Myfield_parser; std::unique_ptr Mzfield_parser; + std::unique_ptr M1xfield_parser; + std::unique_ptr M1yfield_parser; + std::unique_ptr M1zfield_parser; + std::unique_ptr M2xfield_parser; + std::unique_ptr M2yfield_parser; + std::unique_ptr M2zfield_parser; + std::unique_ptr Lxfield_parser; + std::unique_ptr Lyfield_parser; + std::unique_ptr Lzfield_parser; // Parser for H_bias_external on the grid std::unique_ptr Hx_biasfield_parser; std::unique_ptr Hy_biasfield_parser; @@ -342,6 +360,8 @@ public: int mag_LLG_exchange_coupling = 0; // turn off the anisotropy coupling term H_anisotropy in H_eff for the LLG updates int mag_LLG_anisotropy_coupling = 0; + // turn off the antiferromagnetic coupling term H_afm in H_eff for the LGG updates + int mag_LLG_afm_coupling = 0; #endif //! If true, the current is deposited on a nodal grid and then centered onto a staggered grid //! using finite centering of order given by #current_centering_nox, #current_centering_noy, @@ -641,10 +661,36 @@ public: amrex::MultiFab * get_pointer_Bfield_sc_fp (int lev, int direction) const { return Bfield_sc_fp[lev][direction].get(); } #ifdef WARPX_MAG_LLG - amrex::MultiFab * get_pointer_Hfield_fp (int lev, int direction) const { return Hfield_fp[lev][direction].get();} + amrex::MultiFab * get_pointer_Hfield_fp (int lev, int direction) const { return Hfield_fp[lev][direction].get(); } // note "direction" of M means face. For M, each face stores all 3 vector components of M - amrex::MultiFab * get_pointer_Mfield_fp (int lev, int direction) const { return Mfield_fp[lev][direction].get();} - amrex::MultiFab * get_pointer_H_biasfield_fp (int lev, int direction) const { return H_biasfield_fp[lev][direction].get();} + amrex::MultiFab * get_pointer_Mfield_fp (int lev, int direction) const { return Mfield_fp[lev][direction].get(); } + amrex::MultiFab * get_pointer_H_biasfield_fp (int lev, int direction) const { return H_biasfield_fp[lev][direction].get(); } + + amrex::MultiFab * get_pointer_H1field_fp (int lev, int direction) const { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(mag_LLG_afm_coupling == 1, + "get_pointer_H1field_fp called but mag_LLG_afm_coupling is not enabled!"); + return H1field_fp[lev][direction].get(); + } + amrex::MultiFab * get_pointer_H2field_fp (int lev, int direction) const { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(mag_LLG_afm_coupling == 1, + "get_pointer_H2field_fp called but mag_LLG_afm_coupling is not enabled!"); + return H2field_fp[lev][direction].get(); + } + amrex::MultiFab * get_pointer_M1field_fp (int lev, int direction) const { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(mag_LLG_afm_coupling == 1, + "get_pointer_M1field_fp called but mag_LLG_afm_coupling is not enabled!"); + return M1field_fp[lev][direction].get(); + } + amrex::MultiFab * get_pointer_M2field_fp (int lev, int direction) const { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(mag_LLG_afm_coupling == 1, + "get_pointer_M2field_fp called but mag_LLG_afm_coupling is not enabled!"); + return M2field_fp[lev][direction].get(); + } + amrex::MultiFab * get_pointer_Lfield_fp (int lev, int direction) const { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(mag_LLG_afm_coupling == 1, + "get_pointer_Lfield_fp called but mag_LLG_afm_coupling is not enabled!"); + return Lfield_fp[lev][direction].get(); + } #endif amrex::MultiFab * get_pointer_current_fp (int lev, int direction) const { return current_fp[lev][direction].get(); } amrex::MultiFab * get_pointer_current_fp_nodal (int lev, int direction) const { return current_fp_nodal[lev][direction].get(); } @@ -697,6 +743,26 @@ public: const amrex::MultiFab& getHfield_fp (int lev, int direction) {return *Hfield_fp[lev][direction];} const amrex::MultiFab& getMfield_fp (int lev, int direction) {return *Mfield_fp[lev][direction];} const amrex::MultiFab& getH_biasfield_fp (int lev, int direction) {return *H_biasfield_fp[lev][direction];} + const amrex::MultiFab& getH1field_fp (int lev, int direction) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(mag_LLG_afm_coupling == 1, + "getH1field_fp called but mag_LLG_afm_coupling is not enabled!"); + return *H1field_fp[lev][direction]; + } + const amrex::MultiFab& getM1field_fp (int lev, int direction) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(mag_LLG_afm_coupling == 1, + "getM1field_fp called but mag_LLG_afm_coupling is not enabled!"); + return *M1field_fp[lev][direction]; + } + const amrex::MultiFab& getH2field_fp (int lev, int direction) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(mag_LLG_afm_coupling == 1, + "getH2field_fp called but mag_LLG_afm_coupling is not enabled!"); + return *H2field_fp[lev][direction]; + } + const amrex::MultiFab& getM2field_fp (int lev, int direction) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(mag_LLG_afm_coupling == 1, + "getM2field_fp called but mag_LLG_afm_coupling is not enabled!"); + return *M2field_fp[lev][direction]; + } #endif const amrex::MultiFab& getrho_fp (int lev) {return *rho_fp[lev];} const amrex::MultiFab& getphi_fp (int lev) {return *phi_fp[lev];} @@ -825,6 +891,10 @@ public: void MacroscopicEvolveHM_2nd ( amrex::Real dt); void MacroscopicEvolveHM_2nd (int lev, amrex::Real dt); void MacroscopicEvolveHM_2nd (int lev, PatchType patch_type, amrex::Real dt); + + void MacroscopicEvolveHM_AFM_2nd ( amrex::Real dt); + void MacroscopicEvolveHM_AFM_2nd (int lev, amrex::Real dt); + void MacroscopicEvolveHM_AFM_2nd (int lev, PatchType patch_type, amrex::Real dt); #endif /** apply QED correction on electric field @@ -1596,6 +1666,13 @@ private: amrex::Vector, 3 > > Mfield_aux; amrex::Vector, 3 > > Hfield_aux; amrex::Vector, 3 > > H_biasfield_aux; + + // Auxiliary AFM coupling fields (allocated at runtime if mag_LLG_afm_coupling == 1) + amrex::Vector, 3 > > M1field_aux; + amrex::Vector, 3 > > M2field_aux; + amrex::Vector, 3 > > H1field_aux; + amrex::Vector, 3 > > H2field_aux; + amrex::Vector, 3 > > Lfield_aux; #endif amrex::Vector, 3 > > Efield_avg_aux; amrex::Vector, 3 > > Bfield_avg_aux; @@ -1614,6 +1691,13 @@ private: amrex::Vector, 3 > > Mfield_fp; amrex::Vector, 3 > > Hfield_fp; amrex::Vector, 3 > > H_biasfield_fp; + + // AFM coupling fields (allocated at runtime if mag_LLG_afm_coupling == 1) + amrex::Vector, 3 > > M1field_fp; + amrex::Vector, 3 > > M2field_fp; + amrex::Vector, 3 > > H1field_fp; + amrex::Vector, 3 > > H2field_fp; + amrex::Vector, 3 > > Lfield_fp; #endif amrex::Vector, 3 > > Efield_avg_fp; amrex::Vector, 3 > > Bfield_avg_fp; @@ -1690,6 +1774,13 @@ private: amrex::Vector, 3 > > Mfield_cp; amrex::Vector, 3 > > Hfield_cp; amrex::Vector, 3 > > H_biasfield_cp; + + // Coarse patch AFM coupling fields (allocated at runtime if mag_LLG_afm_coupling == 1) + amrex::Vector, 3 > > M1field_cp; + amrex::Vector, 3 > > M2field_cp; + amrex::Vector, 3 > > H1field_cp; + amrex::Vector, 3 > > H2field_cp; + amrex::Vector, 3 > > Lfield_cp; #endif amrex::Vector, 3 > > Efield_avg_cp; amrex::Vector, 3 > > Bfield_avg_cp; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 11a95f6c5..16da4f94f 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -397,6 +397,13 @@ WarpX::WarpX () Mfield_fp.resize(nlevs_max); Hfield_fp.resize(nlevs_max); H_biasfield_fp.resize(nlevs_max); + if (mag_LLG_afm_coupling == 1) { + M1field_fp.resize(nlevs_max); + M2field_fp.resize(nlevs_max); + H1field_fp.resize(nlevs_max); + H2field_fp.resize(nlevs_max); + Lfield_fp.resize(nlevs_max); + } #endif // Only allocate vector potential arrays when using the Magnetostatic Solver @@ -1053,6 +1060,8 @@ WarpX::ReadParameters () pp_warpx.query("mag_LLG_exchange_coupling",mag_LLG_exchange_coupling); // turn on the anisotropy coupling term H_anisotropy for H_eff in the LLG equation pp_warpx.query("mag_LLG_anisotropy_coupling",mag_LLG_anisotropy_coupling); + // turn on the antiferromagnetic coupling term H_afm for H_eff in the LLG equation + pp_warpx.query("mag_LLG_afm_coupling",mag_LLG_afm_coupling); #endif #ifdef WARPX_DIM_RZ @@ -2012,6 +2021,13 @@ WarpX::ClearLevel (int lev) Mfield_fp [lev][i].reset(); Hfield_fp [lev][i].reset(); H_biasfield_fp [lev][i].reset(); + if (mag_LLG_afm_coupling == 1) { + M1field_fp[lev][i].reset(); + M2field_fp[lev][i].reset(); + H1field_fp[lev][i].reset(); + H2field_fp[lev][i].reset(); + Lfield_fp[lev][i].reset(); + } #endif current_store[lev][i].reset(); @@ -2300,6 +2316,24 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm tagname.append("[l=").append(std::to_string(lev)).append("]"); return tagname; }; + if (mag_LLG_afm_coupling == 1) { + AllocInitMultiFab(M1field_fp[lev][0], amrex::convert(ba,Mx_nodal_flag), dm, 3, ngEB, tag("M1field_fp[x]")); + AllocInitMultiFab(M1field_fp[lev][1], amrex::convert(ba,My_nodal_flag), dm, 3, ngEB, tag("M1field_fp[y]")); + AllocInitMultiFab(M1field_fp[lev][2], amrex::convert(ba,Mz_nodal_flag), dm, 3, ngEB, tag("M1field_fp[z]")); + AllocInitMultiFab(M2field_fp[lev][0], amrex::convert(ba,Mx_nodal_flag), dm, 3, ngEB, tag("M2field_fp[x]")); + AllocInitMultiFab(M2field_fp[lev][1], amrex::convert(ba,My_nodal_flag), dm, 3, ngEB, tag("M2field_fp[y]")); + AllocInitMultiFab(M2field_fp[lev][2], amrex::convert(ba,Mz_nodal_flag), dm, 3, ngEB, tag("M2field_fp[z]")); + AllocInitMultiFab(Lfield_fp[lev][0], amrex::convert(ba,Mx_nodal_flag), dm, 3, ngEB, tag("Lfield_fp[x]")); + AllocInitMultiFab(Lfield_fp[lev][1], amrex::convert(ba,My_nodal_flag), dm, 3, ngEB, tag("Lfield_fp[y]")); + AllocInitMultiFab(Lfield_fp[lev][2], amrex::convert(ba,Mz_nodal_flag), dm, 3, ngEB, tag("Lfield_fp[z]")); + + AllocInitMultiFab(H1field_fp[lev][0], amrex::convert(ba,Hx_nodal_flag), dm, 3, ngEB, tag("H1field_fp[x]")); + AllocInitMultiFab(H1field_fp[lev][1], amrex::convert(ba,Hy_nodal_flag), dm, 3, ngEB, tag("H1field_fp[y]")); + AllocInitMultiFab(H1field_fp[lev][2], amrex::convert(ba,Hz_nodal_flag), dm, 3, ngEB, tag("H1field_fp[z]")); + AllocInitMultiFab(H2field_fp[lev][0], amrex::convert(ba,Hx_nodal_flag), dm, 3, ngEB, tag("H2field_fp[x]")); + AllocInitMultiFab(H2field_fp[lev][1], amrex::convert(ba,Hy_nodal_flag), dm, 3, ngEB, tag("H2field_fp[y]")); + AllocInitMultiFab(H2field_fp[lev][2], amrex::convert(ba,Hz_nodal_flag), dm, 3, ngEB, tag("H2field_fp[z]")); + } // // The fine patch