diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index 47060e0ff24..e067513fa06 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -1265,7 +1265,9 @@ Numerics and algorithms The medium for evaluating the Maxwell solver. Available options are : - ``vacuum``: vacuum properties are used in the Maxwell solver. - - ``macroscopic``: macroscopic Maxwell equation is evaluated. If this option is selected, then the corresponding properties of the medium must be provided using ``macroscopic.sigma``, ``macroscopic.epsilon``, and ``macroscopic.mu`` for each case where the initialization style is ``constant``. Otherwise if the initialization style uses the parser, ``macroscopic.sigma_function(x,y,z)``, ``macroscopic.epsilon_function(x,y,z)`` and/or ``macroscopic.mu_function(x,y,z)`` must be provided using the parser initialization style for spatially varying macroscopic properties. + - ``macroscopic``: macroscopic Maxwell equation is evaluated. If this option is selected, then the corresponding properties of the medium must be provided using ``macroscopic.sigma``, + ``macroscopic.epsilon``, and ``macroscopic.mu`` for each case where the initialization style is ``constant``. Otherwise if the initialization style uses the parser, ``macroscopic.sigma_function(x,y,z)``, ``macroscopic.epsilon_function(x,y,z)`` + and/or ``macroscopic.mu_function(x,y,z)`` must be provided using the parser initialization style for spatially varying macroscopic properties. If ``algo.em_solver_medium`` is not specified, ``vacuum`` is the default. @@ -1285,7 +1287,10 @@ Numerics and algorithms * ``macroscopic.sigma``, ``macroscopic.epsilon``, ``macroscopic.mu`` (`double`) To initialize a constant conductivity, permittivity, and permeability of the computational medium, respectively. The default values are the corresponding values - in vacuum. + in vacuum. + In particular, if `USE_LLG=TRUE` in the GNUMakefile and `mag_Ms` is non-zero, ``macroscopic.mu`` must be set to the vacuum value (``mu0``, i.e. 1.25663706212e-06). + If `USE_LLG=TRUE` in the GNUMakefile and `mag_Ms == 0`, LLG evolution is turned off and the macroscopic magnetic permeability, i.e. ``mur * mu0`` should be parsed in ``macroscopic.mu``, with + the dimensionless relative permeability ``mur >=1``. * ``macroscopic.mag_normalized_error`` (`double`; default: `0.1`) The maximum relative amount we let M deviate from Ms before aborting for the LLG equation for saturated cases, i.e., `mag_M_normalization>0`. @@ -2300,6 +2305,9 @@ Solving magnetization using LLG equation ``macroscopic.mag_Ms_function(x,y,z)`` to initialize the saturation magnetization. If ``algo.em_solver_medium`` is set to macroscopic, and ``USE_LLG = TRUE``, then this input property must be provided. + In addition, if `USE_LLG=TRUE` in the GNUMakefile and `mag_Ms` is non-zero, ``macroscopic.mu`` must be set to the vacuum value (``mu0``, i.e. 1.25663706212e-06). + If `USE_LLG=TRUE` in the GNUMakefile and `mag_Ms == 0`, LLG evolution is turned off and the macroscopic magnetic permeability, i.e. ``mur * mu0`` should be parsed in ``macroscopic.mu``, with + the dimensionless relative permeability ``mur >=1``. * ``macroscopic.mag_alpha_init_style`` (string) optional (default is "default") This parameter determines the type of initialization for the Gilbert damping factor diff --git a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary new file mode 100644 index 00000000000..b8b667f45df --- /dev/null +++ b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary @@ -0,0 +1,73 @@ +# This is a E+H (use USE_LLG=TRUE) simulation of a +# plane wave in a periodic domain +################################ +####### GENERAL PARAMETERS ###### +################################# +max_step = 1000 +amr.n_cell = 16 16 128 +amr.max_grid_size = 32 +amr.blocking_factor = 16 +geometry.coord_sys = 0 +geometry.is_periodic = 1 1 0 +geometry.prob_lo = -16.e-6 -16.e-6 -128.e-6 +geometry.prob_hi = 16.e-6 16.e-6 128.e-6 +amr.max_level = 0 + +################################# +############ NUMERICS ########### +################################# +warpx.verbose = 0 +warpx.use_filter = 0 +warpx.cfl = 0.9 +warpx.do_pml = 1 +algo.em_solver_medium = macroscopic # vacuum/macroscopic +algo.macroscopic_sigma_method = laxwendroff # laxwendroff or backwardeuler +macroscopic.sigma_init_style = "parse_sigma_function" # parse or "constant" +macroscopic.sigma_function(x,y,z) = "0.0" +macroscopic.epsilon_init_style = "parse_epsilon_function" # parse or "constant" +macroscopic.epsilon_function(x,y,z) = "epsilon_r*8.8541878128e-12" +macroscopic.mu_init_style = "parse_mu_function" # parse or "constant" +macroscopic.mu_function(x,y,z) = "mu_r*1.25663706212e-06" +particles.nspecies = 0 + +################################# +############ FIELDS ############# +################################# +my_constants.mu_r = 4 +my_constants.epsilon_r = 1 +my_constants.pi = 3.14159265359 +my_constants.L = 141.4213562373095e-6 +my_constants.c = 299792458. +my_constants.wavelength = 64.e-6 + +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) = "1.e7*exp(-z**2/L**2)*cos(2*pi*z/wavelength * sqrt(epsilon_r*mu_r))" +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) = "-1.e7*exp(-z**2/L**2)*cos(2*pi*z/wavelength * sqrt(epsilon_r*mu_r)) / sqrt( mu_r*1.25663706212e-06 / (epsilon_r*8.8541878128e-12) )" +warpx.Hy_external_grid_function(x,y,z) = 0. +warpx.Hz_external_grid_function(x,y,z) = 0. + +# If you want to use a USE_LLG build with Ms=0 everywhere, uncomment these settings +warpx.mag_M_normalization = 1 +macroscopic.mag_Ms_init_style = parse_mag_Ms_function +macroscopic.mag_Ms_function(x,y,z) = "1.4e5 * (z>64e-6)" + +macroscopic.mag_alpha_init_style = "parse_mag_alpha_function" # parse or "constant" +macroscopic.mag_alpha_function(x,y,z) = "0.0058 * (z>64e-6)" # 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 * (z>64e-6)" # gyromagnetic ratio is constant for electrons in all materials + +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.4e5 * (z>64e-6)" +warpx.Mz_external_grid_function(x,y,z) = 0. + +# Diagnostics +diagnostics.diags_names = plt +plt.period = 10 +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.diag_type = Full diff --git a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_noMs b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_noMs index 58bb683488a..60c7ba3884e 100644 --- a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_noMs +++ b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_noMs @@ -1,10 +1,9 @@ # This is a E+H (use USE_LLG=TRUE) simulation of a # plane wave in a periodic domain - ################################ ####### GENERAL PARAMETERS ###### ################################# -max_step = 500 +max_step = 1000 amr.n_cell = 16 16 128 amr.max_grid_size = 32 amr.blocking_factor = 16 @@ -21,35 +20,34 @@ warpx.verbose = 0 warpx.use_filter = 0 warpx.cfl = 0.9 warpx.do_pml = 1 - algo.em_solver_medium = macroscopic # vacuum/macroscopic algo.macroscopic_sigma_method = laxwendroff # laxwendroff or backwardeuler macroscopic.sigma_init_style = "parse_sigma_function" # parse or "constant" macroscopic.sigma_function(x,y,z) = "0.0" macroscopic.epsilon_init_style = "parse_epsilon_function" # parse or "constant" -macroscopic.epsilon_function(x,y,z) = "8.8541878128e-12" +macroscopic.epsilon_function(x,y,z) = "epsilon_r*8.8541878128e-12" macroscopic.mu_init_style = "parse_mu_function" # parse or "constant" -macroscopic.mu_function(x,y,z) = "1.25663706212e-06" - +macroscopic.mu_function(x,y,z) = "mu_r*1.25663706212e-06" particles.nspecies = 0 ################################# ############ FIELDS ############# ################################# - +my_constants.mu_r = 1 +my_constants.epsilon_r = 1 my_constants.pi = 3.14159265359 my_constants.L = 141.4213562373095e-6 my_constants.c = 299792458. my_constants.wavelength = 64.e-6 warpx.E_ext_grid_init_style = parse_E_ext_grid_function -warpx.Ez_external_grid_function(x,y,z) = 0. warpx.Ex_external_grid_function(x,y,z) = 0. -warpx.Ey_external_grid_function(x,y,z) = "1.e5*exp(-z**2/L**2)*cos(2*pi*z/wavelength)" +warpx.Ey_external_grid_function(x,y,z) = "1.e5*exp(-z**2/L**2)*cos(2*pi*z/wavelength * sqrt(epsilon_r*mu_r))" +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)= "-1.e5*exp(-z**2/L**2)*cos(2*pi*z/wavelength)/(c*1.25663706212e-06)" -warpx.Hy_external_grid_function(x,y,z)= 0. +warpx.Hx_external_grid_function(x,y,z) = "-1.e5*exp(-z**2/L**2)*cos(2*pi*z/wavelength * sqrt(epsilon_r*mu_r)) / sqrt( mu_r*1.25663706212e-06 / (epsilon_r*8.8541878128e-12) )" +warpx.Hy_external_grid_function(x,y,z) = 0. warpx.Hz_external_grid_function(x,y,z) = 0. # If you want to use a USE_LLG build with Ms=0 everywhere, uncomment these settings diff --git a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs index d9393c00681..32cf5ef0e8e 100644 --- a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs +++ b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs @@ -1,10 +1,9 @@ -# This is a E+B (use USE_LLG=FALSE, ie., no M field) simulation of a +# This is a E+B (use USE_LLG=FALSE) simulation of a # plane wave in a periodic domain - ################################ ####### GENERAL PARAMETERS ###### ################################# -max_step = 500 +max_step = 1000 amr.n_cell = 16 16 128 amr.max_grid_size = 32 amr.blocking_factor = 16 @@ -21,35 +20,34 @@ warpx.verbose = 0 warpx.use_filter = 0 warpx.cfl = 0.9 warpx.do_pml = 1 - algo.em_solver_medium = macroscopic # vacuum/macroscopic algo.macroscopic_sigma_method = laxwendroff # laxwendroff or backwardeuler macroscopic.sigma_init_style = "parse_sigma_function" # parse or "constant" macroscopic.sigma_function(x,y,z) = "0.0" macroscopic.epsilon_init_style = "parse_epsilon_function" # parse or "constant" -macroscopic.epsilon_function(x,y,z) = "8.8541878128e-12" +macroscopic.epsilon_function(x,y,z) = "epsilon_r*8.8541878128e-12" macroscopic.mu_init_style = "parse_mu_function" # parse or "constant" -macroscopic.mu_function(x,y,z) = "1.25663706212e-06" - +macroscopic.mu_function(x,y,z) = "mu_r*1.25663706212e-06" particles.nspecies = 0 ################################# ############ FIELDS ############# ################################# - +my_constants.mu_r = 1 +my_constants.epsilon_r = 1 my_constants.pi = 3.14159265359 my_constants.L = 141.4213562373095e-6 my_constants.c = 299792458. my_constants.wavelength = 64.e-6 warpx.E_ext_grid_init_style = parse_E_ext_grid_function -warpx.Ez_external_grid_function(x,y,z) = 0. warpx.Ex_external_grid_function(x,y,z) = 0. -warpx.Ey_external_grid_function(x,y,z) = "1.e5*exp(-z**2/L**2)*cos(2*pi*z/wavelength)" +warpx.Ey_external_grid_function(x,y,z) = "1.e5*exp(-z**2/L**2)*cos(2*pi*z/wavelength * sqrt(epsilon_r*mu_r))" +warpx.Ez_external_grid_function(x,y,z) = 0. warpx.B_ext_grid_init_style = parse_B_ext_grid_function -warpx.Bx_external_grid_function(x,y,z)= "-1.e5*exp(-z**2/L**2)*cos(2*pi*z/wavelength)/c" -warpx.By_external_grid_function(x,y,z)= 0. +warpx.Bx_external_grid_function(x,y,z) = "-1.e5*exp(-z**2/L**2)*cos(2*pi*z/wavelength * sqrt(epsilon_r*mu_r))/c * sqrt(epsilon_r*mu_r)" +warpx.By_external_grid_function(x,y,z) = 0. warpx.Bz_external_grid_function(x,y,z) = 0. # Diagnostics diff --git a/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt b/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt index a792fbbfa14..34f01cd8e03 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt +++ b/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt @@ -19,7 +19,7 @@ if(WarpX_MAG_LLG) MacroscopicEvolveHM.cpp MacroscopicEvolveM_2nd.cpp MacroscopicEvolveM.cpp - EvolveHPML.cpp + MacroscopicEvolveHPML.cpp ) endif() add_subdirectory(MacroscopicProperties) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index a454816f79a..faa31bb2ffe 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -146,9 +146,11 @@ class FiniteDifferenceSolver #ifdef WARPX_MAG_LLG - void EvolveHPML ( std::array< amrex::MultiFab*, 3 > Hfield, + void MacroscopicEvolveHPML ( std::array< amrex::MultiFab*, 3 > Hfield, std::array< amrex::MultiFab*, 3 > const Efield, amrex::Real const dt, + std::unique_ptr const& macroscopic_properties, + amrex::MultiFab* const mu_mf, const bool dive_cleaning); #endif @@ -303,10 +305,12 @@ class FiniteDifferenceSolver #ifdef WARPX_MAG_LLG template< typename T_Algo > - void EvolveHPMLCartesian ( - std::array< amrex::MultiFab*, 3 > Bfield, + void MacroscopicEvolveHPMLCartesian ( + std::array< amrex::MultiFab*, 3 > Hfield, std::array< amrex::MultiFab*, 3 > const Efield, amrex::Real const dt, + std::unique_ptr const& macroscopic_properties, + amrex::MultiFab* const mu_mf, const bool dive_cleaning); #endif diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveHPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp similarity index 70% rename from Source/FieldSolver/FiniteDifferenceSolver/EvolveHPML.cpp rename to Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp index 71f3aa0567f..96d62a0ec62 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveHPML.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp @@ -1,9 +1,3 @@ -/* Copyright 2020 Remi Lehe - * - * This file is part of WarpX. - * - * License: BSD-3-Clause-LBNL - */ #include "Utils/WarpXAlgorithmSelection.H" #include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H" @@ -17,36 +11,39 @@ #include "BoundaryConditions/PMLComponent.H" #include #include +#include "Utils/CoarsenIO.H" using namespace amrex; #ifdef WARPX_MAG_LLG /** - * \brief Update the H field, over one timestep + * \brief Update Hfield in PML region */ -void FiniteDifferenceSolver::EvolveHPML ( +void FiniteDifferenceSolver::MacroscopicEvolveHPML ( std::array< amrex::MultiFab*, 3 > Hfield, std::array< amrex::MultiFab*, 3 > const Efield, amrex::Real const dt, + std::unique_ptr const& macroscopic_properties, + amrex::MultiFab* const mu_mf, const bool dive_cleaning) { // Select algorithm (The choice of algorithm is a runtime option, // but we compile code for each algorithm, using templates) #ifdef WARPX_DIM_RZ - amrex::ignore_unused(Hfield, Efield, dt); + amrex::ignore_unused(Hfield, Efield, dt, macroscopic_properties, mu_mf); amrex::Abort("PML are not implemented in cylindrical geometry."); #else if (m_do_nodal) { - - EvolveHPMLCartesian (Hfield, Efield, dt, dive_cleaning); + + MacroscopicEvolveHPMLCartesian ( Hfield, Efield, dt, macroscopic_properties, mu_mf, dive_cleaning); } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) { - EvolveHPMLCartesian (Hfield, Efield, dt, dive_cleaning); + MacroscopicEvolveHPMLCartesian ( Hfield, Efield, dt, macroscopic_properties, mu_mf, dive_cleaning); } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) { - EvolveHPMLCartesian (Hfield, Efield, dt, dive_cleaning); - + MacroscopicEvolveHPMLCartesian ( Hfield, Efield, dt, macroscopic_properties, mu_mf, dive_cleaning); + } else { amrex::Abort("EvolveHPML: Unknown algorithm"); } @@ -57,12 +54,20 @@ void FiniteDifferenceSolver::EvolveHPML ( #ifndef WARPX_DIM_RZ template -void FiniteDifferenceSolver::EvolveHPMLCartesian ( +void FiniteDifferenceSolver::MacroscopicEvolveHPMLCartesian ( std::array< amrex::MultiFab*, 3 > Hfield, std::array< amrex::MultiFab*, 3 > const Efield, amrex::Real const dt, + std::unique_ptr const& macroscopic_properties, + amrex::MultiFab* const mu_mf, const bool dive_cleaning) { + amrex::GpuArray const& mu_stag = macroscopic_properties->mu_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; + // Loop through the grids, and over the tiles within each grid #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) @@ -90,12 +95,17 @@ void FiniteDifferenceSolver::EvolveHPMLCartesian ( Box const& tby = mfi.tilebox(Hfield[1]->ixType().ixType()); Box const& tbz = mfi.tilebox(Hfield[2]->ixType().ixType()); - amrex::Real mu0_inv = 1._rt/PhysConst::mu0; + // starting component to interpolate macro properties to Hx, Hy, Hz locations + const int scomp = 0; + // mu_mf will be imported but will only be called at grids where Ms == 0 + 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){ + + Real mu_inv = 1._rt/CoarsenIO::Interp( mu_arr, mu_stag, Hx_stag, macro_cr, i, j, k, scomp); amrex::Real UpwardDz_Ey_yy = 0._rt; amrex::Real UpwardDy_Ez_zz = 0._rt; @@ -105,19 +115,21 @@ void FiniteDifferenceSolver::EvolveHPMLCartesian ( UpwardDy_Ez_zz = T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k, PMLComp::zz); } - Hx(i, j, k, PMLComp::xz) += mu0_inv * dt * ( + Hx(i, j, k, PMLComp::xz) += mu_inv * dt * ( T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, PMLComp::yx) + T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, PMLComp::yz) + UpwardDz_Ey_yy); - Hx(i, j, k, PMLComp::xy) -= mu0_inv * dt * ( + Hx(i, j, k, PMLComp::xy) -= mu_inv * dt * ( T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k, PMLComp::zx) + T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k, PMLComp::zy) + UpwardDy_Ez_zz); }, [=] AMREX_GPU_DEVICE (int i, int j, int k){ - + + Real mu_inv = 1._rt/CoarsenIO::Interp( mu_arr, mu_stag, Hy_stag, macro_cr, i, j, k, scomp); + amrex::Real UpwardDx_Ez_zz = 0._rt; amrex::Real UpwardDz_Ex_xx = 0._rt; if (dive_cleaning) @@ -126,12 +138,12 @@ void FiniteDifferenceSolver::EvolveHPMLCartesian ( UpwardDz_Ex_xx = T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k, PMLComp::xx); } - Hy(i, j, k, PMLComp::yx) += mu0_inv * dt * ( + Hy(i, j, k, PMLComp::yx) += mu_inv * dt * ( T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k, PMLComp::zx) + T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k, PMLComp::zy) + UpwardDx_Ez_zz); - Hy(i, j, k, PMLComp::yz) -= mu0_inv * dt * ( + Hy(i, j, k, PMLComp::yz) -= mu_inv * dt * ( UpwardDz_Ex_xx + T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k, PMLComp::xy) + T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k, PMLComp::xz)); @@ -139,6 +151,8 @@ void FiniteDifferenceSolver::EvolveHPMLCartesian ( [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Real mu_inv = 1._rt/CoarsenIO::Interp( mu_arr, mu_stag, Hz_stag, macro_cr, i, j, k, scomp); + amrex::Real UpwardDy_Ex_xx = 0._rt; amrex::Real UpwardDx_Ey_yy = 0._rt; if (dive_cleaning) @@ -147,12 +161,12 @@ void FiniteDifferenceSolver::EvolveHPMLCartesian ( UpwardDx_Ey_yy = T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k, PMLComp::yy); } - Hz(i, j, k, PMLComp::zy) += mu0_inv * dt * ( + Hz(i, j, k, PMLComp::zy) += mu_inv * dt * ( UpwardDy_Ex_xx + T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k, PMLComp::xy) + T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k, PMLComp::xz) ); - Hz(i, j, k, PMLComp::zx) -= mu0_inv * dt * ( + Hz(i, j, k, PMLComp::zx) -= mu_inv * dt * ( T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k, PMLComp::yx) + T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k, PMLComp::yz) + UpwardDx_Ey_yy); diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H index 6c4795957cc..c5dd42f04b9 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H @@ -58,6 +58,12 @@ public: amrex::GpuArray My_IndexType; /** Gpu Vector with index type of the Mz multifab */ amrex::GpuArray Mz_IndexType; + /** Gpu Vector with index type of the Hx multifab */ + amrex::GpuArray Hx_IndexType; + /** Gpu Vector with index type of the Hy multifab */ + amrex::GpuArray Hy_IndexType; + /** Gpu Vector with index type of the Hz multifab */ + amrex::GpuArray Hz_IndexType; #endif /** Gpu Vector with index type of coarsening ratio with default value (1,1,1) */ amrex::GpuArray macro_cr_ratio; diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp index 8a58e1ea9f5..2621c1d52f3 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp @@ -237,6 +237,9 @@ MacroscopicProperties::InitData () IntVect Mx_stag = warpx.getMfield_fp(0,0).ixType().toIntVect(); // face-centered IntVect My_stag = warpx.getMfield_fp(0,1).ixType().toIntVect(); IntVect Mz_stag = warpx.getMfield_fp(0,2).ixType().toIntVect(); + IntVect Hx_stag = warpx.getHfield_fp(0,0).ixType().toIntVect(); // face-centered + IntVect Hy_stag = warpx.getHfield_fp(0,1).ixType().toIntVect(); + IntVect Hz_stag = warpx.getHfield_fp(0,2).ixType().toIntVect(); #endif for ( int idim = 0; idim < AMREX_SPACEDIM; ++idim) { sigma_IndexType[idim] = sigma_stag[idim]; @@ -252,6 +255,9 @@ MacroscopicProperties::InitData () Mx_IndexType[idim] = Mx_stag[idim]; My_IndexType[idim] = My_stag[idim]; Mz_IndexType[idim] = Mz_stag[idim]; + Hx_IndexType[idim] = Hx_stag[idim]; + Hy_IndexType[idim] = Hy_stag[idim]; + Hz_IndexType[idim] = Hz_stag[idim]; #endif macro_cr_ratio[idim] = 1; } @@ -269,6 +275,9 @@ MacroscopicProperties::InitData () Mx_IndexType[2] = 0; My_IndexType[2] = 0; Mz_IndexType[2] = 0; + Hx_IndexType[2] = 0; + Hy_IndexType[2] = 0; + Hz_IndexType[2] = 0; #endif macro_cr_ratio[2] = 1; #endif diff --git a/Source/FieldSolver/FiniteDifferenceSolver/Make.package b/Source/FieldSolver/FiniteDifferenceSolver/Make.package index d269911a527..c20f4531c72 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/Make.package +++ b/Source/FieldSolver/FiniteDifferenceSolver/Make.package @@ -8,7 +8,7 @@ CEXE_sources += MacroscopicEvolveE.cpp #ifdef WARPX_MAG_LLG CEXE_sources += MacroscopicEvolveHM.cpp CEXE_sources += MacroscopicEvolveHM_2nd.cpp -CEXE_sources += EvolveHPML.cpp +CEXE_sources += MacroscopicEvolveHPML.cpp #endif CEXE_sources += EvolveBPML.cpp diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index d0c474a7a85..9822595b558 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -427,11 +427,21 @@ WarpX::MacroscopicEvolveHM (int lev, PatchType patch_type, amrex::Real a_dt) { // 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); + m_fdtd_solver_fp[lev]->MacroscopicEvolveHPML( + pml[lev]->GetH_fp(), + pml[lev]->GetE_fp(), + a_dt, + m_macroscopic_properties, + pml[lev]->Getmu_fp(), + 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 ); + m_fdtd_solver_cp[lev]->MacroscopicEvolveHPML( + pml[lev]->GetH_cp(), + pml[lev]->GetE_cp(), + a_dt, + m_macroscopic_properties, + pml[lev]->Getmu_cp(), + WarpX::do_dive_cleaning); } } } @@ -470,11 +480,21 @@ WarpX::MacroscopicEvolveHM_2nd (int lev, PatchType patch_type, amrex::Real a_dt) // 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 ); + m_fdtd_solver_fp[lev]->MacroscopicEvolveHPML( + pml[lev]->GetH_fp(), + pml[lev]->GetE_fp(), + a_dt, + m_macroscopic_properties, + pml[lev]->Getmu_fp(), + 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 ); + m_fdtd_solver_cp[lev]->MacroscopicEvolveHPML( + pml[lev]->GetH_cp(), + pml[lev]->GetE_cp(), + a_dt, + m_macroscopic_properties, + pml[lev]->Getmu_cp(), + WarpX::do_dive_cleaning); } } }