diff --git a/Examples/Tests/Macroscopic_Maxwell/PSSW_halfCos_exci/analysis_getline.bash b/Examples/Tests/Macroscopic_Maxwell/PSSW_halfCos_exci/analysis_getline.bash new file mode 100644 index 000000000..7733d89f5 --- /dev/null +++ b/Examples/Tests/Macroscopic_Maxwell/PSSW_halfCos_exci/analysis_getline.bash @@ -0,0 +1,11 @@ +#!/bin/bash +FILE="$1" +LINE_NO=$2 +i=0 +while read line; do + i=$(( i + 1 )) + case $i in $LINE_NO) echo "$line"; break;; esac +done <"$FILE" + +### how to run this script: +### $ bash squashTime.bash FILE LINE_NO \ No newline at end of file diff --git a/Examples/Tests/Macroscopic_Maxwell/PSSW_halfCos_exci/analysis_readfield.bash b/Examples/Tests/Macroscopic_Maxwell/PSSW_halfCos_exci/analysis_readfield.bash new file mode 100644 index 000000000..a66156be4 --- /dev/null +++ b/Examples/Tests/Macroscopic_Maxwell/PSSW_halfCos_exci/analysis_readfield.bash @@ -0,0 +1,6 @@ +#!/bin/bash +for timestep in $(seq 0 500 5000000 ); do + newnumber='000000000'${timestep} # get number, pack with zeros + newnumber=${newnumber:(-9)} # the last seven characters + ./WritePlotfileToASCII3d.gnu.ex infile=diags/plt${newnumber} | tee raw_data/${newnumber}.txt +done diff --git a/Examples/Tests/Macroscopic_Maxwell/PSSW_halfCos_exci/analysis_squashTime.bash b/Examples/Tests/Macroscopic_Maxwell/PSSW_halfCos_exci/analysis_squashTime.bash new file mode 100644 index 000000000..e57827115 --- /dev/null +++ b/Examples/Tests/Macroscopic_Maxwell/PSSW_halfCos_exci/analysis_squashTime.bash @@ -0,0 +1,6 @@ +#!/bin/bash +for timestep in $(seq 0 500 5000000); do + newnumber='000000000'${timestep} # get number, pack with zeros + newnumber=${newnumber:(-9)} # the last seven characters + bash getLine.bash raw_data/${newnumber}.txt 128 | tee -a Mx_xface_center_128.txt +done diff --git a/Examples/Tests/Macroscopic_Maxwell/PSSW_halfCos_exci/inputs_LLG_PSSW_asymm_exci b/Examples/Tests/Macroscopic_Maxwell/PSSW_halfCos_exci/inputs_LLG_PSSW_asymm_exci new file mode 100644 index 000000000..d76f3cd3b --- /dev/null +++ b/Examples/Tests/Macroscopic_Maxwell/PSSW_halfCos_exci/inputs_LLG_PSSW_asymm_exci @@ -0,0 +1,129 @@ +# set up to duplicate perpendicular standing spin wave for exchange boundary condition demo +# max step +max_step = 5000000 + +# number of grid points +amr.n_cell = n_cellx n_celly n_cellz + +# Maximum allowable size of each subdomain +amr.max_grid_size = grid_number +amr.blocking_factor = 8 + +amr.max_level = 0 + +# Geometry +geometry.coord_sys = 0 # 0: Cartesian +geometry.prob_lo = -width/2 -length/2 -thickness +geometry.prob_hi = width/2 length/2 thickness + +# Boundary condition +boundary.field_lo = periodic periodic pec +boundary.field_hi = periodic periodic pec +warpx.serialize_ics = 1 + +# Verbosity +warpx.verbose = 1 +warpx.use_filter = 0 +warpx.cfl = 1000.0 # We can use a large cfl number here because the Maxwell equations are not included. +warpx.mag_time_scheme_order = 2 # default 1 +warpx.mag_M_normalization = 1 +warpx.mag_LLG_coupling = 0 +warpx.mag_LLG_exchange_coupling = 1 +warpx.mag_LLG_anisotropy_coupling = 0 + +# Algorithms +algo.current_deposition = esirkepov +algo.em_solver_medium = macroscopic # vacuum/macroscopic + +my_constants.thickness = 345.0e-9 +my_constants.width = 6900.0e-9 +my_constants.length = 6900.0e-9 +my_constants.n_cellx = 8 +my_constants.n_celly = 8 +my_constants.n_cellz = 256 +my_constants.grid_number = 256 + +my_constants.epsilon_r = 13.0 +my_constants.pi = 3.14159265359 +my_constants.mu0 = 1.25663706212e-6 +my_constants.epsilon0 = 8.84e-12 +my_constants.Ms_ga = 1750. + +my_constants.flag_none = 0 +my_constants.flag_hs = 1 +my_constants.flag_ss = 2 + +algo.macroscopic_sigma_method = laxwendroff # laxwendroff or backwardeuler +macroscopic.sigma_function(x,y,z) = "0.0" +macroscopic.epsilon_function(x,y,z) = "epsilon0 * epsilon_r" +macroscopic.mu_function(x,y,z) = "mu0" + +#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) = "(z < thickness/2 + 1e-12) * (z > -thickness/2 - 1e-12) * Ms_ga*1000/4/pi" # in unit A/m + +macroscopic.mag_alpha_init_style = "parse_mag_alpha_function" # parse or "constant" +macroscopic.mag_alpha_function(x,y,z) = "(z < thickness/2 + 1e-12) * (z > -thickness/2 - 1e-12) * 0.005" # alpha is unitless, typical values range from 1e-3 ~ 1e-5 + +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) = "(z < thickness/2 + 1e-12) * (z > -thickness/2 - 1e-12) * 3.1e-12" # YIG + +#macroscopic.mag_anisotropy_init_style = "parse_mag_anisotropy_function" # parse or "constant" +#macroscopic.mag_anisotropy_function(x,y,z) = "0.0" # Ku = mu0 * Ms * H_anis / 2.0, all fields in SI units, H_anis = 20 Oe + +macroscopic.mag_max_iter = 100 # maximum number of M iteration in each time step +macroscopic.mag_tol = 1.e-6 # 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 +#macroscopic.mag_LLG_anisotropy_axis = 0.0 1.0 0.0 + +################################# +############ FIELDS ############# +################################# +my_constants.Hy_bias_oe = 300 # in Oe +my_constants.Hx_rf_oe = 3 # in Oe +my_constants.c = 299792458. +my_constants.frequency = 2.5e9 # frequency of input microwave H +my_constants.TP = 1/frequency +my_constants.flag_none = 0 +my_constants.flag_hs = 1 +my_constants.flag_ss = 2 + +#warpx.E_excitation_on_grid_style = "parse_E_excitation_grid_function" +#warpx.Ex_excitation_grid_function(x,y,z,t) = "0.0" +#warpx.Ey_excitation_grid_function(x,y,z,t) = "0.0" +#warpx.Ez_excitation_grid_function(x,y,z,t) = "0.0" +#warpx.Ex_excitation_flag_function(x,y,z) = "flag_none" +# warpx.Ey_excitation_flag_function(x,y,z) = "flag_ss * (z > -domain_size/2.0 + thickness - domain_size/grid_number/2.0) * (z <= -domain_size/2.0 + thickness + domain_size/grid_number/2.0)" +#warpx.Ey_excitation_flag_function(x,y,z) = "flag_none" +#warpx.Ez_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) = "Hx_rf_oe * 1000/4/pi * (exp(-(t-3*TP)**2/(2*TP**2))*cos(2*pi*frequency*t))" +warpx.Hx_bias_excitation_grid_function(x,y,z,t) = "Hx_rf_oe * 1000/4/pi * (exp(-(t-3*TP)**2/(2*TP**2))*cos(2*pi*frequency*t)) * cos(z/thickness * pi)" + +#warpx.Hx_bias_excitation_grid_function(x,y,z,t) = "0.0" +warpx.Hy_bias_excitation_grid_function(x,y,z,t) = "Hy_bias_oe * 1000/4/pi" +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_hs" +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)= "0.0" +warpx.My_external_grid_function(x,y,z)= " (z < thickness/2 + 1e-12) * (z > -thickness/2 - 1e-12) * Ms_ga*1000/4/pi " +warpx.Mz_external_grid_function(x,y,z)= "0.0" + + +# Diagnostics +diagnostics.diags_names = plt +plt.intervals = 500 +plt.diag_type = Full +plt.fields_to_plot = Mx_xface Mx_yface Mx_zface My_xface My_yface My_zface Mz_xface Mz_yface Mz_zface +plt.file_min_digits = 9 diff --git a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_exchangeBC b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_exchangeBC new file mode 100644 index 000000000..dc6244912 --- /dev/null +++ b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_exchangeBC @@ -0,0 +1,88 @@ +################################ +####### GENERAL PARAMETERS ###### +################################# +max_step = 6000 +amr.n_cell = 64 64 64 # number of cells spanning the domain in each coordinate direction at level 0 +amr.max_grid_size = 32 # maximum size of each AMReX box, used to decompose the domain +amr.blocking_factor = 16 +geometry.coord_sys = 0 + +geometry.prob_lo = -10e-9 -10e-9 -10.0e-9 +geometry.prob_hi = 10e-9 10e-9 10.0e-9 +boundary.field_lo = periodic periodic periodic +boundary.field_hi = periodic periodic periodic +amr.max_level = 0 + +################################# +############ NUMERICS ########### +################################# +warpx.verbose = 1 +warpx.use_filter = 0 +warpx.cfl = 1000 +warpx.mag_time_scheme_order = 2 # default 1 +warpx.mag_M_normalization = 1 # 1 is saturated +warpx.mag_LLG_coupling = 0 +warpx.mag_LLG_exchange_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" +macroscopic.mu_function(x,y,z) = "1.25663706212e-06" + +my_constants.mag_lo_x = -5.e-9 +my_constants.mag_hi_x = 5.e-9 +my_constants.mag_lo_y = -5.e-9 +my_constants.mag_hi_y = 5.e-9 +my_constants.mag_lo_z = -5.e-9 +my_constants.mag_hi_z = 5.e-9 + +#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.4e5 * (x >= mag_lo_x - 1e-12) * (x <= mag_hi_x + 1e-12) * (y >= mag_lo_y - 1e-12) * (y <= mag_hi_y + 1e-12) * (z >= mag_lo_z - 1e-12) * (z <= mag_hi_z + 1e-12)" # in unit A/m, equal to 1750 Gauss; Ms must be nonzero for LLG + +macroscopic.mag_alpha_init_style = "parse_mag_alpha_function" # parse or "constant" +macroscopic.mag_alpha_function(x,y,z) = "0.058" # * (x >= mag_lo_x - 1e-12) * (x <= mag_hi_x + 1e-12) * (y >= mag_lo_y - 1e-12) * (y <= mag_hi_y + 1e-12) * (z >= mag_lo_z - 1e-12) * (z <= mag_hi_z + 1e-12)" # alpha is unitless, calculated from linewidth Delta_H = 40 Oersted + +macroscopic.mag_gamma_init_style = "parse_mag_gamma_function" # parse or "constant" +macroscopic.mag_gamma_function(x,y,z) = "-1.759e11" # * (x >= mag_lo_x - 1e-12) * (x <= mag_hi_x + 1e-12) * (y >= mag_lo_y - 1e-12) * (y <= mag_hi_y + 1e-12) * (z >= mag_lo_z - 1e-12) * (z <= mag_hi_z + 1e-12)" # 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) = "3.76e-12" # * (x >= mag_lo_x) * (x <= mag_hi_x) * (y >= mag_lo_y - 1e-12) * (y <= mag_hi_y + 1e-12) * (z >= mag_lo_z - 1e-12) * (z <= mag_hi_z + 1e-12)" # exchange coupling constanit; Must be non-zero when exchange coupling is ON + +macroscopic.mag_max_iter = 100 # maximum number of M iteration in each time step +macroscopic.mag_tol = 1.e-6 # 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.E_ext_grid_init_style = parse_E_ext_grid_function +warpx.Ex_external_grid_function(x,y,z) = 0. +warpx.Ey_external_grid_function(x,y,z) = 0. +warpx.Ez_external_grid_function(x,y,z) = 0. + +warpx.H_ext_grid_init_style = parse_H_ext_grid_function +warpx.Hx_external_grid_function(x,y,z)= 0. +warpx.Hy_external_grid_function(x,y,z) = 0. +warpx.Hz_external_grid_function(x,y,z) = 0. + +#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)= "3.7e4" # in A/m, equal to 464 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)= "1.4e5 * (x >= mag_lo_x - 1e-12) * (x <= mag_hi_x + 1e-12) * (y >= mag_lo_y - 1e-12) * (y < 0.) * (z >= mag_lo_z - 1e-12) * (z <= mag_hi_z + 1e-12)" +warpx.My_external_grid_function(x,y,z) = 0. +warpx.Mz_external_grid_function(x,y,z) = "1.4e5 * (x >= mag_lo_x - 1e-12) * (x <= mag_hi_x + 1e-12) * (y >= 0.) * (y <= mag_hi_y + 1e-12) * (z >= mag_lo_z - 1e-12) * (z <= mag_hi_z + 1e-12)" + +#Diagnostics +diagnostics.diags_names = plt +plt.intervals = 100 +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 +plt.plot_raw_fields = 0 diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H index ad0445e85..f3081f733 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H @@ -15,6 +15,9 @@ #include #include #include +#include "Utils/WarpXUtil.H" +#include "WarpX.H" +#include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H" #include #include @@ -239,6 +242,109 @@ struct CartesianYeeAlgorithm { #endif +#ifdef WARPX_MAG_LLG + + /** + * Perform divergence of gradient along x on M field when exchange coupling is on */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real LaplacianDx_Mag ( + amrex::Array4 const& F, + amrex::Real const * const coefs_x, int const n_coefs_x, amrex::Real const Ms_lo_x, amrex::Real const Ms_hi_x, + int const i, int const j, int const k, int const ncomp=0, int const nodality=0) { + + amrex::Real const inv_dx = coefs_x[0]; + if (nodality == 0){ // at x face (normal face). dM/dx = 0 + if (Ms_hi_x == 0.){ + return 0.5 * inv_dx * inv_dx * (8. * F(i-1, j, k, ncomp) - F(i-2, j, k, ncomp) - 7. * F(i, j, k, ncomp)); + } else if (Ms_lo_x == 0){ + return 0.5 * inv_dx * inv_dx * (8. * F(i+1, j, k, ncomp) - F(i+2, j, k, ncomp) - 7. * F(i, j, k, ncomp)); + } else { + return inv_dx*(UpwardDx(F, coefs_x, n_coefs_x, i, j, k, ncomp) - DownwardDx(F, coefs_x, n_coefs_x, i, j, k, ncomp)); + } + } else { // at y or z faces + if (Ms_hi_x == 0.){ + return inv_dx*(0. - DownwardDx(F, coefs_x, n_coefs_x, i, j, k, ncomp)); + } else if (Ms_lo_x == 0.){ + return inv_dx*(UpwardDx(F, coefs_x, n_coefs_x, i, j, k, ncomp) - 0.); + } else { + return inv_dx*(UpwardDx(F, coefs_x, n_coefs_x, i, j, k, ncomp) - DownwardDx(F, coefs_x, n_coefs_x, i, j, k, ncomp)); + } + } + } + + /** + * Perform divergence of gradient along y on M field when exchange coupling is on*/ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real LaplacianDy_Mag ( + amrex::Array4 const& F, + amrex::Real const * const coefs_y, int const n_coefs_y, amrex::Real const Ms_lo_y, amrex::Real const Ms_hi_y, + int const i, int const j, int const k, int const ncomp=0, int const nodality=0) { + + amrex::Real const inv_dy = coefs_y[0]; + if (nodality == 1){ // y face (normal face), dM/dy = 0 + if (Ms_hi_y == 0.){ + return 0.5 * inv_dy * inv_dy * (8. * F(i, j-1, k, ncomp) - F(i, j-2, k, ncomp) - 7. * F(i, j, k, ncomp)); + } else if (Ms_lo_y == 0.) { + return 0.5 * inv_dy * inv_dy * (8. * F(i, j+1, k, ncomp) - F(i, j+2, k, ncomp) - 7. * F(i, j, k, ncomp)); + } else { + return inv_dy*(UpwardDy(F, coefs_y, n_coefs_y, i, j, k, ncomp) - DownwardDy(F, coefs_y, n_coefs_y, i, j, k, ncomp)); + } + + } else { // x or z faces + if (Ms_hi_y == 0.){ + return inv_dy*(0. - DownwardDy(F, coefs_y, n_coefs_y, i, j, k, ncomp)); + } else if (Ms_lo_y == 0.) { + return inv_dy*(UpwardDy(F, coefs_y, n_coefs_y, i, j, k, ncomp) - 0.); + } else { + return inv_dy*(UpwardDy(F, coefs_y, n_coefs_y, i, j, k, ncomp) - DownwardDy(F, coefs_y, n_coefs_y, i, j, k, ncomp)); + } + } + } + + /** + * Perform divergence of gradient along z on M field when exchange coupling is on*/ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real LaplacianDz_Mag ( + amrex::Array4 const& F, + amrex::Real const * const coefs_z, int const n_coefs_z, amrex::Real const Ms_lo_z, amrex::Real const Ms_hi_z, + int const i, int const j, int const k, int const ncomp=0, int const nodality=0) { + + amrex::Real const inv_dz = coefs_z[0]; + if (nodality == 2){ // z face (normal face) dM/dz = 0 + if ( Ms_hi_z == 0.) { + return 0.5 * inv_dz * inv_dz * (8. * F(i, j, k-1, ncomp) - F(i, j, k-2, ncomp) - 7. * F(i, j, k, ncomp)); + } else if ( Ms_lo_z == 0.){ + return 0.5 * inv_dz * inv_dz * (8. * F(i, j, k+1, ncomp) - F(i, j, k+2, ncomp) - 7. * F(i, j, k, ncomp)); + } else { + return inv_dz*(UpwardDz(F, coefs_z, n_coefs_z, i, j, k, ncomp) - DownwardDz(F, coefs_z, n_coefs_z, i, j, k, ncomp)); + } + + } else { // x or y face + if ( Ms_hi_z == 0.) { + return inv_dz*(0. - DownwardDz(F, coefs_z, n_coefs_z, i, j, k, ncomp)); + } else if ( Ms_lo_z == 0.){ + return inv_dz*(UpwardDz(F, coefs_z, n_coefs_z, i, j, k, ncomp) - 0.); + } else { + return inv_dz*(UpwardDz(F, coefs_z, n_coefs_z, i, j, k, ncomp) - DownwardDz(F, coefs_z, n_coefs_z, i, j, k, ncomp)); + } + } + } + + /** + * Compute the sum to get Laplacian of M field when exchange coupling is on*/ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real Laplacian_Mag ( + amrex::Array4 const& F, + amrex::Real const * const coefs_x, amrex::Real const * const coefs_y, amrex::Real const * const coefs_z, + int const n_coefs_x, int const n_coefs_y, int const n_coefs_z, + amrex::Real const Ms_lo_x, amrex::Real const Ms_hi_x, amrex::Real const Ms_lo_y, amrex::Real const Ms_hi_y, amrex::Real const Ms_lo_z, amrex::Real const Ms_hi_z, + int const i, int const j, int const k, int const ncomp=0, int const nodality=0) { + + return LaplacianDx_Mag(F, coefs_x, n_coefs_x, Ms_lo_x, Ms_hi_x, i, j, k, ncomp, nodality) + LaplacianDy_Mag(F, coefs_y, n_coefs_y, Ms_lo_y, Ms_hi_y, i, j, k, ncomp, nodality) + LaplacianDz_Mag(F, coefs_z, n_coefs_z, Ms_lo_z, Ms_hi_z, i, j, k, ncomp, nodality); + } + +#endif + }; #endif // WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_YEE_H_ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index 88aa868b4..77c778170 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -142,6 +142,7 @@ class FiniteDifferenceSolver std::unique_ptr const ¯oscopic_properties); void MacroscopicEvolveHM_2nd ( + int lev, std::array, 3> &Mfield, std::array, 3> &Hfield, // H Maxwell std::array< std::unique_ptr, 3>& Bfield, @@ -327,6 +328,7 @@ class FiniteDifferenceSolver template< typename T_Algo > void MacroscopicEvolveHMCartesian_2nd( + int lev, std::array, 3> &Mfield, std::array, 3> &Hfield, // H Maxwell std::array< std::unique_ptr, 3 >& Bfield, diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHM.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHM.cpp index ed751b951..ab3d2dc1c 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHM.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHM.cpp @@ -154,9 +154,12 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian( Box const &tbz = mfi.tilebox(Hfield[2]->ixType().toIntVect()); // 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(); - amrex::Real const * const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr(); - amrex::Real const * const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr(); + 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, @@ -188,9 +191,17 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian( // 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); - Hx_eff += H_exchange_coeff * T_Algo::Laplacian(M_old_xface, coefs_x, coefs_y, coefs_z, i, j, k, 0); - Hy_eff += H_exchange_coeff * T_Algo::Laplacian(M_old_xface, coefs_x, coefs_y, coefs_z, i, j, k, 1); - Hz_eff += H_exchange_coeff * T_Algo::Laplacian(M_old_xface, coefs_x, coefs_y, coefs_z, i, j, k, 2); + + 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); + + Hx_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_old_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 + Hy_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_old_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 + Hz_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_old_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 } if (mag_anisotropy_coupling == 1){ @@ -297,9 +308,17 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian( // 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); - Hx_eff += H_exchange_coeff * T_Algo::Laplacian(M_old_yface, coefs_x, coefs_y, coefs_z, i, j, k, 0); - Hy_eff += H_exchange_coeff * T_Algo::Laplacian(M_old_yface, coefs_x, coefs_y, coefs_z, i, j, k, 1); - Hz_eff += H_exchange_coeff * T_Algo::Laplacian(M_old_yface, coefs_x, coefs_y, coefs_z, i, j, k, 2); + + 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); + + Hx_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_old_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 + Hy_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_old_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 + Hz_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_old_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 } if (mag_anisotropy_coupling == 1){ @@ -407,9 +426,17 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian( // 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); - Hx_eff += H_exchange_coeff * T_Algo::Laplacian(M_old_zface, coefs_x, coefs_y, coefs_z, i, j, k, 0); - Hy_eff += H_exchange_coeff * T_Algo::Laplacian(M_old_zface, coefs_x, coefs_y, coefs_z, i, j, k, 1); - Hz_eff += H_exchange_coeff * T_Algo::Laplacian(M_old_zface, coefs_x, coefs_y, coefs_z, i, j, k, 2); + + 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); + + Hx_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_old_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 + Hy_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_old_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 + Hz_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_old_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 } if (mag_anisotropy_coupling == 1){ @@ -450,7 +477,7 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian( // z component on z-faces of grid M_zface(i, j, k, 2) += dt * (PhysConst::mu0 * mag_gammaL) * (M_old_zface(i, j, k, 0) * Hy_eff - M_old_zface(i, j, k, 1) * Hx_eff) + dt * Gil_damp * (M_old_zface(i, j, k, 0) * (M_old_zface(i, j, k, 2) * Hx_eff - M_old_zface(i, j, k, 0) * Hz_eff) - - M_old_zface(i, j, k, 1) * (M_old_zface(i, j, k, 1) * Hz_eff - M_old_yface(i, j, k, 2) * Hy_eff)); + - M_old_zface(i, j, k, 1) * (M_old_zface(i, j, k, 1) * Hz_eff - M_old_zface(i, j, k, 2) * Hy_eff)); // 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 diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHM_2nd.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHM_2nd.cpp index ad861c85d..622d97400 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHM_2nd.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHM_2nd.cpp @@ -29,6 +29,7 @@ using namespace amrex; void FiniteDifferenceSolver::MacroscopicEvolveHM_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> &Mfield, // Mfield contains three components MultiFab std::array, 3> &Hfield, std::array, 3> &Bfield, @@ -38,7 +39,7 @@ void FiniteDifferenceSolver::MacroscopicEvolveHM_2nd( std::unique_ptr const ¯oscopic_properties) { if (m_fdtd_algo == MaxwellSolverAlgo::Yee){ - MacroscopicEvolveHMCartesian_2nd(Mfield, Hfield, Bfield, H_biasfield, Efield, dt, macroscopic_properties); + MacroscopicEvolveHMCartesian_2nd(lev, Mfield, Hfield, Bfield, H_biasfield, Efield, dt, macroscopic_properties); } else { amrex::Abort("Only yee algorithm is compatible for M updates."); } @@ -47,6 +48,7 @@ void FiniteDifferenceSolver::MacroscopicEvolveHM_2nd( #ifdef WARPX_MAG_LLG template void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian_2nd( + int lev, std::array, 3> &Mfield, std::array, 3> &Hfield, std::array, 3> &Bfield, @@ -167,9 +169,12 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian_2nd( 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(); - amrex::Real const * const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr(); - amrex::Real const * const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr(); + 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, @@ -201,9 +206,17 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian_2nd( // 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); - Hx_eff += H_exchange_coeff * T_Algo::Laplacian(M_xface, coefs_x, coefs_y, coefs_z, i, j, k, 0); - Hy_eff += H_exchange_coeff * T_Algo::Laplacian(M_xface, coefs_x, coefs_y, coefs_z, i, j, k, 1); - Hz_eff += H_exchange_coeff * T_Algo::Laplacian(M_xface, coefs_x, coefs_y, coefs_z, i, j, k, 2); + + 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); + + Hx_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_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 + Hy_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_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 + Hz_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_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 } if (mag_anisotropy_coupling == 1){ @@ -278,9 +291,17 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian_2nd( // 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); - Hx_eff += H_exchange_coeff * T_Algo::Laplacian(M_yface, coefs_x, coefs_y, coefs_z, i, j, k, 0); - Hy_eff += H_exchange_coeff * T_Algo::Laplacian(M_yface, coefs_x, coefs_y, coefs_z, i, j, k, 1); - Hz_eff += H_exchange_coeff * T_Algo::Laplacian(M_yface, coefs_x, coefs_y, coefs_z, i, j, k, 2); + + 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); + + Hx_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_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 + Hy_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_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 + Hz_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_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 } if (mag_anisotropy_coupling == 1){ @@ -355,9 +376,17 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian_2nd( // 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); - Hx_eff += H_exchange_coeff * T_Algo::Laplacian(M_zface, coefs_x, coefs_y, coefs_z, i, j, k, 0); - Hy_eff += H_exchange_coeff * T_Algo::Laplacian(M_zface, coefs_x, coefs_y, coefs_z, i, j, k, 1); - Hz_eff += H_exchange_coeff * T_Algo::Laplacian(M_zface, coefs_x, coefs_y, coefs_z, i, j, k, 2); + + 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); + + Hx_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_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 + Hy_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_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 + Hz_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_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 } if (mag_anisotropy_coupling == 1){ @@ -491,9 +520,12 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian_2nd( 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(); - amrex::Real const * const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr(); - amrex::Real const * const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr(); + 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, @@ -525,9 +557,17 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian_2nd( // 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); - Hx_eff += H_exchange_coeff * T_Algo::Laplacian(M_prev_xface, coefs_x, coefs_y, coefs_z, i, j, k, 0); - Hy_eff += H_exchange_coeff * T_Algo::Laplacian(M_prev_xface, coefs_x, coefs_y, coefs_z, i, j, k, 1); - Hz_eff += H_exchange_coeff * T_Algo::Laplacian(M_prev_xface, coefs_x, coefs_y, coefs_z, i, j, k, 2); + + 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); + + Hx_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_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 + Hy_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_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 + Hz_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_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 } if (mag_anisotropy_coupling == 1){ @@ -632,9 +672,17 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian_2nd( // 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); - Hx_eff += H_exchange_coeff * T_Algo::Laplacian(M_prev_yface, coefs_x, coefs_y, coefs_z, i, j, k, 0); - Hy_eff += H_exchange_coeff * T_Algo::Laplacian(M_prev_yface, coefs_x, coefs_y, coefs_z, i, j, k, 1); - Hz_eff += H_exchange_coeff * T_Algo::Laplacian(M_prev_yface, coefs_x, coefs_y, coefs_z, i, j, k, 2); + + 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); + + Hx_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_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 + Hy_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_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 + Hz_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_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 } if (mag_anisotropy_coupling == 1){ @@ -740,9 +788,17 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian_2nd( // 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); - Hx_eff += H_exchange_coeff * T_Algo::Laplacian(M_prev_zface, coefs_x, coefs_y, coefs_z, i, j, k, 0); - Hy_eff += H_exchange_coeff * T_Algo::Laplacian(M_prev_zface, coefs_x, coefs_y, coefs_z, i, j, k, 1); - Hz_eff += H_exchange_coeff * T_Algo::Laplacian(M_prev_zface, coefs_x, coefs_y, coefs_z, i, j, k, 2); + + 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); + + Hx_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_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 + Hy_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_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 + Hz_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_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 } if (mag_anisotropy_coupling == 1){ @@ -1027,9 +1083,11 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian_2nd( } } else{ - // Copy Mfield to Mfield_previous + 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); } } diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 88a10b801..3b42afa89 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -819,7 +819,7 @@ WarpX::MacroscopicEvolveHM_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_2nd(Mfield_fp[lev], Hfield_fp[lev], Bfield_fp[lev], H_biasfield_fp[lev], Efield_fp[lev], + m_fdtd_solver_fp[lev]->MacroscopicEvolveHM_2nd(lev, Mfield_fp[lev], Hfield_fp[lev], Bfield_fp[lev], H_biasfield_fp[lev], Efield_fp[lev], a_dt, m_macroscopic_properties); } else {