From 2e46b508f89205c6d7048f0060e8ceb4bdc346e0 Mon Sep 17 00:00:00 2001 From: RevathiJambunathan Date: Mon, 7 Dec 2020 17:55:33 -0800 Subject: [PATCH 01/30] initialize and allocate macro pml mfs --- Source/BoundaryConditions/PML.H | 19 ++++- Source/BoundaryConditions/PML.cpp | 115 +++++++++++++++++++++++++++++- 2 files changed, 132 insertions(+), 2 deletions(-) diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index ec54f19aa56..96e5af865ed 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -107,7 +107,8 @@ public: int do_dive_cleaning, int do_moving_window, int pml_has_particles, int do_pml_in_domain, const amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(), - const amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector()); + const amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector(), + const int lev = 0); void ComputePMLFactors (amrex::Real dt); @@ -125,6 +126,14 @@ public: amrex::MultiFab* GetF_fp (); amrex::MultiFab* GetF_cp (); + // return macroscopic pml multifabs + amrex::MultiFab* Geteps_fp(); + amrex::MultiFab* Geteps_cp(); + amrex::MultiFab* Getmu_fp(); + amrex::MultiFab* Getmu_cp(); + amrex::MultiFab* Getsigma_fp(); + amrex::MultiFab* Getsigma_cp(); + const MultiSigmaBox& GetMultiSigmaBox_fp () const { return *sigba_fp; } @@ -201,6 +210,14 @@ private: std::array,3> pml_H_cp; #endif + // macroproperties in PML, eps-permittivity, mu-permeability, sigma-condictivity + std::unique_ptr pml_eps_fp; + std::unique_ptr pml_mu_fp; + std::unique_ptr pml_sigma_fp; + std::unique_ptr pml_eps_cp; + std::unique_ptr pml_mu_cp; + std::unique_ptr pml_sigma_cp; + std::unique_ptr sigba_fp; std::unique_ptr sigba_cp; diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 3ff7b647184..daee9760bbf 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -11,6 +11,7 @@ #include "Utils/WarpXConst.H" #include "Utils/WarpXAlgorithmSelection.H" #include "WarpX.H" +//#include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H" #include #include @@ -432,7 +433,8 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& /*grid_dm*/, #endif int do_dive_cleaning, int do_moving_window, int /*pml_has_particles*/, int do_pml_in_domain, - const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi) + const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi, + const int lev) : m_geom(geom), m_cgeom(cgeom) { @@ -514,6 +516,42 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& /*grid_dm*/, pml_H_fp[2] = std::make_unique(amrex::convert( ba, WarpX::GetInstance().getHfield_fp(0,2).ixType().toIntVect() ), dm, 2, ngb ); #endif + + if (WarpX::em_solver_medium == MediumForEM::Macroscopic) { + // Allocating macroproperties in pml at cell-centers + pml_eps_fp = std::make_unique(ba, dm, 1, nge); + pml_mu_fp = std::make_unique(ba, dm, 1, nge); + pml_sigma_fp = std::make_unique(ba, dm, 1, nge); + + // Initializing macroparameter multifab // + auto& warpx = WarpX::GetInstance(); + auto& macroscopic_properties = warpx.m_macroscopic_properties; + + // Initialize sigma, conductivity + if (macroscopic_properties->m_sigma_s == "constant") { + pml_sigma_fp->setVal(macroscopic_properties->m_sigma); + } else if (macroscopic_properties->m_sigma_s == "parse_sigma_function") { + macroscopic_properties->InitializeMacroMultiFabUsingParser(pml_sigma_fp.get(), + getParser(macroscopic_properties->m_sigma_parser), lev); + } + + // Initialize epsilon, permittivity + if (macroscopic_properties->m_epsilon_s == "constant") { + pml_eps_fp->setVal(macroscopic_properties->m_epsilon); + } else if (macroscopic_properties->m_epsilon_s == "parse_epsilon_function") { + macroscopic_properties->InitializeMacroMultiFabUsingParser(pml_eps_fp.get(), + getParser(macroscopic_properties->m_epsilon_parser), lev); + } + + // Initialize mu, permeability + if (macroscopic_properties->m_mu_s == "constant") { + pml_mu_fp->setVal(macroscopic_properties->m_mu); + } else if (macroscopic_properties->m_sigma_s == "parse_mu_function") { + macroscopic_properties->InitializeMacroMultiFabUsingParser(pml_mu_fp.get(), + getParser(macroscopic_properties->m_mu_parser), lev); + } + + } pml_E_fp[0]->setVal(0.0); pml_E_fp[1]->setVal(0.0); @@ -602,6 +640,44 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& /*grid_dm*/, WarpX::GetInstance().getHfield_cp(1,2).ixType().toIntVect() ), cdm, 2, ngb ); #endif + + // Allocating macroproperties in pml at cell-centers + if (WarpX::em_solver_medium == MediumForEM::Macroscopic) { + pml_eps_cp = std::make_unique(cba, dm, 1, nge); + pml_mu_cp = std::make_unique(cba, dm, 1, nge); + pml_sigma_cp = std::make_unique(cba, dm, 1, nge); + + // Initializing macroparameter multifab // + auto& warpx = WarpX::GetInstance(); + auto& macroscopic_properties = warpx.m_macroscopic_properties; + + // Initialize sigma, conductivity + if (macroscopic_properties->m_sigma_s == "constant") { + pml_sigma_cp->setVal(macroscopic_properties->m_sigma); + } else if (macroscopic_properties->m_sigma_s == "parse_sigma_function") { + macroscopic_properties->InitializeMacroMultiFabUsingParser(pml_sigma_cp.get(), + getParser(macroscopic_properties->m_sigma_parser), lev); + } + + // Initialize epsilon, permittivity + if (macroscopic_properties->m_epsilon_s == "constant") { + pml_eps_cp->setVal(macroscopic_properties->m_epsilon); + } else if (macroscopic_properties->m_epsilon_s == "parse_epsilon_function") { + macroscopic_properties->InitializeMacroMultiFabUsingParser(pml_eps_cp.get(), + getParser(macroscopic_properties->m_epsilon_parser), lev); + } + + // Initialize mu, permeability + if (macroscopic_properties->m_mu_s == "constant") { + pml_mu_cp->setVal(macroscopic_properties->m_mu); + } else if (macroscopic_properties->m_sigma_s == "parse_mu_function") { + macroscopic_properties->InitializeMacroMultiFabUsingParser(pml_mu_cp.get(), + getParser(macroscopic_properties->m_mu_parser), lev); + } + + + } + pml_E_cp[0]->setVal(0.0); pml_E_cp[1]->setVal(0.0); pml_E_cp[2]->setVal(0.0); @@ -805,6 +881,43 @@ PML::GetF_cp () return pml_F_cp.get(); } +// Return macroscopic pml multifabs +amrex::MultiFab* +PML::Geteps_fp() +{ + return pml_eps_fp.get(); +} + +amrex::MultiFab* +PML::Getmu_fp() +{ + return pml_mu_fp.get(); +} + +amrex::MultiFab* +PML::Getsigma_fp() +{ + return pml_sigma_fp.get(); +} + +amrex::MultiFab* +PML::Geteps_cp() +{ + return pml_eps_cp.get(); +} + +amrex::MultiFab* +PML::Getmu_cp() +{ + return pml_mu_cp.get(); +} + +amrex::MultiFab* +PML::Getsigma_cp() +{ + return pml_sigma_cp.get(); +} + void PML::ExchangeB (const std::array& B_fp, const std::array& B_cp, From 97abb25149279a369ad8de3fab60f2d809b682da Mon Sep 17 00:00:00 2001 From: RevathiJambunathan Date: Mon, 7 Dec 2020 17:56:15 -0800 Subject: [PATCH 02/30] add level when initializing PML constructor --- Source/Initialization/WarpXInitData.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 706ce246653..a3a8c259d0a 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -176,7 +176,7 @@ WarpX::InitPML () #endif do_dive_cleaning, do_moving_window, pml_has_particles, do_pml_in_domain, - do_pml_Lo_corrected, do_pml_Hi); + do_pml_Lo_corrected, do_pml_Hi,0); for (int lev = 1; lev <= finest_level; ++lev) { amrex::IntVect do_pml_Lo_MR = amrex::IntVect::TheUnitVector(); @@ -194,7 +194,7 @@ WarpX::InitPML () #endif do_dive_cleaning, do_moving_window, pml_has_particles, do_pml_in_domain, - do_pml_Lo_MR, amrex::IntVect::TheUnitVector()); + do_pml_Lo_MR, amrex::IntVect::TheUnitVector(), lev); } } } From 63bbf84f34caa239b1da3b607554ee4e7a65f5ed Mon Sep 17 00:00:00 2001 From: RevathiJambunathan Date: Mon, 7 Dec 2020 17:57:29 -0800 Subject: [PATCH 03/30] change to public so we can call from PML constructor --- .../MacroscopicProperties.H | 28 ++++++++++++------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H index 4308b75cc27..6c4795957cc 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H @@ -184,7 +184,6 @@ public: }; #endif //closes ifdef MAG_LLG -private: /** Conductivity, sigma, of the medium */ amrex::Real m_sigma = 0.0; /** Permittivity, epsilon, of the medium */ @@ -192,6 +191,24 @@ private: /** Permeability, mu, of the medium */ amrex::Real m_mu = PhysConst::mu0; + /** Stores initialization type for conductivity : constant or parser */ + std::string m_sigma_s = "constant"; + /** Stores initialization type for permittivity : constant or parser */ + std::string m_epsilon_s = "constant"; + /** Stores initialization type for permeability : constant or parser */ + std::string m_mu_s = "constant"; + + /** Parser Wrappers */ + // The ParserWrapper struct is constructed to safely use the GpuParser class + // that can essentially be though of as a raw pointer. The GpuParser does + // not have an explicit destructor and the AddPlasma subroutines handle the GpuParser + // in a safe way. The ParserWrapper struct is used to avoid memory leak + // in the EB parser functions. + std::unique_ptr > m_sigma_parser; + std::unique_ptr > m_epsilon_parser; + std::unique_ptr > m_mu_parser; +private: + #ifdef WARPX_MAG_LLG // preferred to use this tag multiple times for different variable types to keep formatting consistent /** Saturation magnetization, only applies for magnetic materials */ amrex::Real m_mag_Ms; @@ -228,12 +245,6 @@ private: std::unique_ptr m_mag_gamma_mf; #endif - /** Stores initialization type for conductivity : constant or parser */ - std::string m_sigma_s = "constant"; - /** Stores initialization type for permittivity : constant or parser */ - std::string m_epsilon_s = "constant"; - /** Stores initialization type for permeability : constant or parser */ - std::string m_mu_s = "constant"; #ifdef WARPX_MAG_LLG std::string m_mag_Ms_s; @@ -257,9 +268,6 @@ private: // not have an explicit destructor and the AddPlasma subroutines handle the GpuParser // in a safe way. The ParserWrapper struct is used to avoid memory leak // in the EB parser functions. - std::unique_ptr > m_sigma_parser; - std::unique_ptr > m_epsilon_parser; - std::unique_ptr > m_mu_parser; #ifdef WARPX_MAG_LLG std::unique_ptr > m_mag_Ms_parser; std::unique_ptr > m_mag_alpha_parser; From 955e67db28078dbed013e60345e01b2c0fb0d637 Mon Sep 17 00:00:00 2001 From: RevathiJambunathan Date: Mon, 7 Dec 2020 17:58:29 -0800 Subject: [PATCH 04/30] Add new file for macroscopic evolveE pml account for material --- .../FiniteDifferenceSolver/CMakeLists.txt | 1 + .../MacroscopicEvolveEPML.cpp | 302 ++++++++++++++++++ .../FiniteDifferenceSolver/Make.package | 1 + 3 files changed, 304 insertions(+) create mode 100644 Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveEPML.cpp diff --git a/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt b/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt index ffcfe4271d1..a792fbbfa14 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt +++ b/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt @@ -9,6 +9,7 @@ target_sources(WarpX EvolveFPML.cpp FiniteDifferenceSolver.cpp MacroscopicEvolveE.cpp + MacroscopicEvolveEPML.cpp ) if(WarpX_MAG_LLG) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveEPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveEPML.cpp new file mode 100644 index 00000000000..b3a6db30bb2 --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveEPML.cpp @@ -0,0 +1,302 @@ +#include "Utils/WarpXAlgorithmSelection.H" +#include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H" +#ifdef WARPX_DIM_RZ +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H" +#else +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H" +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H" +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H" +#endif +#include "BoundaryConditions/PML.H" +#include "BoundaryConditions/PML_current.H" +#include "BoundaryConditions/PMLComponent.H" +#include "Utils/CoarsenIO.H" +#include "Utils/WarpXConst.H" +#include +#include + +using namespace amrex; + +/** + * \brief Update the E field, over one timestep + */ +void FiniteDifferenceSolver::MacroscopicEvolveEPML ( + std::array< amrex::MultiFab*, 3 > Efield, +#ifndef WARPX_MAG_LLG + std::array< amrex::MultiFab*, 3 > const Bfield, +#else + std::array< amrex::MultiFab*, 3 > const Hfield, +#endif + std::array< amrex::MultiFab*, 3 > const Jfield, + amrex::MultiFab* const Ffield, + MultiSigmaBox const& sigba, + amrex::Real const dt, bool pml_has_particles, + std::unique_ptr const& macroscopic_properties, + amrex::MultiFab* const eps_mf, + amrex::MultiFab* const mu_mf, + amrex::MultiFab* const sigma_mf) { + + // 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(Efield, Bfield, Jfield, Ffield, sigba, dt, pml_has_particles); + amrex::Abort("PML are not implemented in cylindrical geometry."); +#else + if (m_do_nodal) { + + amrex::Abort("Macro E-push is not implemented for nodal, yet."); + + } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) { + + if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::LaxWendroff) { + MacroscopicEvolveEPMLCartesian ( + Efield, +#ifndef WARPX_MAG_LLG + Bfield, +#else + Hfield, +#endif + Jfield, Ffield, sigba, dt, pml_has_particles, + macroscopic_properties, eps_mf, mu_mf, sigma_mf ); + } + else if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::BackwardEuler) { + MacroscopicEvolveEPMLCartesian ( + Efield, +#ifndef WARPX_MAG_LLG + Bfield, +#else + Hfield, +#endif + Jfield, Ffield, sigba, dt, pml_has_particles, + macroscopic_properties, eps_mf, mu_mf, sigma_mf ); + } + + } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) { + // Note :: Macroscopic Evolve E for PML is the same for CKC and Yee + if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::LaxWendroff) { + MacroscopicEvolveEPMLCartesian ( + Efield, +#ifndef WARPX_MAG_LLG + Bfield, +#else + Hfield, +#endif + Jfield, Ffield, sigba, dt, pml_has_particles, + macroscopic_properties, eps_mf, mu_mf, sigma_mf ); + } + else if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::BackwardEuler) { + MacroscopicEvolveEPMLCartesian ( + Efield, +#ifndef WARPX_MAG_LLG + Bfield, +#else + Hfield, +#endif + Jfield, Ffield, sigba, dt, pml_has_particles, + macroscopic_properties, eps_mf, mu_mf, sigma_mf ); + } + + } else { + amrex::Abort("Unknown algorithm"); + } +#endif +} + + +#ifndef WARPX_DIM_RZ + +template +void FiniteDifferenceSolver::MacroscopicEvolveEPMLCartesian ( + std::array< amrex::MultiFab*, 3 > Efield, +#ifndef WARPX_MAG_LLG + std::array< amrex::MultiFab*, 3 > const Bfield, +#else + std::array< amrex::MultiFab*, 3 > const Hfield, +#endif + std::array< amrex::MultiFab*, 3 > const Jfield, + amrex::MultiFab* const Ffield, + MultiSigmaBox const& sigba, + amrex::Real const dt, bool pml_has_particles, + std::unique_ptr const& macroscopic_properties, + amrex::MultiFab* const eps_mf, + amrex::MultiFab* const mu_mf, + amrex::MultiFab* const sigma_mf ) { + + Real c2 = PhysConst::c * PhysConst::c; + +#ifdef WARPX_MAG_LLG + c2 *= PhysConst::mu0; +#endif + + // Index type required for calling CoarsenIO::Interp to interpolate macroscopic + // properties from their respective staggering to the Ex, Ey, Ez locations + amrex::GpuArray const& sigma_stag = macroscopic_properties->sigma_IndexType; + amrex::GpuArray const& epsilon_stag = macroscopic_properties->epsilon_IndexType; + amrex::GpuArray const& Ex_stag = macroscopic_properties->Ex_IndexType; + amrex::GpuArray const& Ey_stag = macroscopic_properties->Ey_IndexType; + amrex::GpuArray const& Ez_stag = macroscopic_properties->Ez_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()) +#endif + for ( MFIter mfi(*Efield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + + // Extract field data for this grid/tile + Array4 const& Ex = Efield[0]->array(mfi); + Array4 const& Ey = Efield[1]->array(mfi); + Array4 const& Ez = Efield[2]->array(mfi); +#ifndef WARPX_MAG_LLG + Array4 const& Bx = Bfield[0]->array(mfi); + Array4 const& By = Bfield[1]->array(mfi); + Array4 const& Bz = Bfield[2]->array(mfi); +#endif + + // material macroscopic properties + Array4 const& sigma_arr = sigma_mf->array(mfi); + Array4 const& eps_arr = eps_mf->array(mfi); + Array4 const& mu_arr = mu_mf->array(mfi); + + + // Extract stencil coefficients + Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr(); + int const n_coefs_x = m_stencil_coefs_x.size(); + Real const * const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr(); + int const n_coefs_y = m_stencil_coefs_y.size(); + Real const * const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr(); + int const n_coefs_z = m_stencil_coefs_z.size(); + +#ifndef WARPX_MAG_LLG + FieldAccessorMacroscopic const Hx(Bx, mu_arr); + FieldAccessorMacroscopic const Hy(By, mu_arr); + FieldAccessorMacroscopic const Hz(Bz, mu_arr); +#else + Array4 const Hx = Hfield[0]->array(mfi); + Array4 const Hy = Hfield[1]->array(mfi); + Array4 const Hz = Hfield[2]->array(mfi); +#endif + + // Extract tileboxes for which to loop + Box const& tex = mfi.tilebox(Efield[0]->ixType().toIntVect()); + Box const& tey = mfi.tilebox(Efield[1]->ixType().toIntVect()); + Box const& tez = mfi.tilebox(Efield[2]->ixType().toIntVect()); + // starting component to interpolate macro properties to Ex, Ey, Ez locations + const int scomp = 0; + + // Loop over the cells and update the fields + amrex::ParallelFor(tex, tey, tez, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + // Interpolate conductivity, sigma, to Ex position on the grid + amrex::Real const sigma_interp = CoarsenIO::Interp( sigma_arr, sigma_stag, + Ex_stag, macro_cr, i, j, k, scomp); + // Interpolated permittivity, epsilon, to Ex position on the grid + amrex::Real const epsilon_interp = CoarsenIO::Interp( eps_arr, epsilon_stag, + Ex_stag, macro_cr, i, j, k, scomp); + amrex::Real alpha = T_MacroAlgo::alpha( sigma_interp, epsilon_interp, dt); + amrex::Real beta = T_MacroAlgo::beta( sigma_interp, epsilon_interp, dt); + + Ex(i, j, k, PMLComp::xz) = alpha * Ex(i, j, k, PMLComp::xz) - beta * ( + T_Algo::DownwardDz(Hy, coefs_z, n_coefs_z, i, j, k, PMLComp::yx) + + T_Algo::DownwardDz(Hy, coefs_z, n_coefs_z, i, j, k, PMLComp::yz) ); + Ex(i, j, k, PMLComp::xy) = alpha * Ex(i, j, k, PMLComp::xy) + beta * ( + T_Algo::DownwardDy(Hz, coefs_y, n_coefs_y, i, j, k, PMLComp::zx) + + T_Algo::DownwardDy(Hz, coefs_y, n_coefs_y, i, j, k, PMLComp::zy) ); + }, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + amrex::Real const sigma_interp = CoarsenIO::Interp( sigma_arr, sigma_stag, + Ey_stag, macro_cr, i, j, k, scomp); + amrex::Real const epsilon_interp = CoarsenIO::Interp( eps_arr, epsilon_stag, + Ey_stag, macro_cr, i, j, k, scomp); + amrex::Real alpha = T_MacroAlgo::alpha( sigma_interp, epsilon_interp, dt); + amrex::Real beta = T_MacroAlgo::beta( sigma_interp, epsilon_interp, dt); + + Ey(i, j, k, PMLComp::yx) = alpha * Ey(i, j, k, PMLComp::yx) - beta * ( + T_Algo::DownwardDx(Hz, coefs_x, n_coefs_x, i, j, k, PMLComp::zx) + + T_Algo::DownwardDx(Hz, coefs_x, n_coefs_x, i, j, k, PMLComp::zy) ); + Ey(i, j, k, PMLComp::yz) = alpha * Ey(i, j, k, PMLComp::yz) + beta * ( + T_Algo::DownwardDz(Hx, coefs_z, n_coefs_z, i, j, k, PMLComp::xy) + + T_Algo::DownwardDz(Hx, coefs_z, n_coefs_z, i, j, k, PMLComp::xz) ); + }, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + amrex::Real const sigma_interp = CoarsenIO::Interp( sigma_arr, sigma_stag, + Ez_stag, macro_cr, i, j, k, scomp); + amrex::Real const epsilon_interp = CoarsenIO::Interp( eps_arr, epsilon_stag, + Ez_stag, macro_cr, i, j, k, scomp); + amrex::Real alpha = T_MacroAlgo::alpha( sigma_interp, epsilon_interp, dt); + amrex::Real beta = T_MacroAlgo::beta( sigma_interp, epsilon_interp, dt); + + Ez(i, j, k, PMLComp::zy) = alpha * Ez(i, j, k, PMLComp::zy) - beta * ( + T_Algo::DownwardDy(Hx, coefs_y, n_coefs_y, i, j, k, PMLComp::xy) + + T_Algo::DownwardDy(Hx, coefs_y, n_coefs_y, i, j, k, PMLComp::xz) ); + Ez(i, j, k, PMLComp::zx) = alpha * Ez(i, j, k, PMLComp::zx) + beta * ( + T_Algo::DownwardDx(Hy, coefs_x, n_coefs_x, i, j, k, PMLComp::yx) + + T_Algo::DownwardDx(Hy, coefs_x, n_coefs_x, i, j, k, PMLComp::yz) ); + } + + ); + + // Update the E field in the PML, using the current + // deposited by the particles in the PML + if (pml_has_particles) { + + // Extract field data for this grid/tile + Array4 const& Jx = Jfield[0]->array(mfi); + Array4 const& Jy = Jfield[1]->array(mfi); + Array4 const& Jz = Jfield[2]->array(mfi); + const Real* sigmaj_x = sigba[mfi].sigma[0].data(); + const Real* sigmaj_y = sigba[mfi].sigma[1].data(); + const Real* sigmaj_z = sigba[mfi].sigma[2].data(); + int const x_lo = sigba[mfi].sigma[0].lo(); +#if (AMREX_SPACEDIM == 3) + int const y_lo = sigba[mfi].sigma[1].lo(); + int const z_lo = sigba[mfi].sigma[2].lo(); +#else + int const y_lo = 0; + int const z_lo = sigba[mfi].sigma[1].lo(); +#endif + + amrex::ParallelFor( tex, tey, tez, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + // Interpolate conductivity, sigma, to Ex position on the grid + amrex::Real const sigma_interp = CoarsenIO::Interp( sigma_arr, sigma_stag, + Ex_stag, macro_cr, i, j, k, scomp); + // Interpolated permittivity, epsilon, to Ex position on the grid + amrex::Real const epsilon_interp = CoarsenIO::Interp( eps_arr, epsilon_stag, + Ex_stag, macro_cr, i, j, k, scomp); + amrex::Real beta = T_MacroAlgo::beta( sigma_interp, epsilon_interp, dt); + + push_ex_pml_current(i, j, k, Ex, Jx, + sigmaj_y, sigmaj_z, y_lo, z_lo, beta); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + amrex::Real const sigma_interp = CoarsenIO::Interp( sigma_arr, sigma_stag, + Ey_stag, macro_cr, i, j, k, scomp); + amrex::Real const epsilon_interp = CoarsenIO::Interp( eps_arr, epsilon_stag, + Ey_stag, macro_cr, i, j, k, scomp); + amrex::Real beta = T_MacroAlgo::beta( sigma_interp, epsilon_interp, dt); + + push_ey_pml_current(i, j, k, Ey, Jy, + sigmaj_x, sigmaj_z, x_lo, z_lo, beta); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + amrex::Real const sigma_interp = CoarsenIO::Interp( sigma_arr, sigma_stag, + Ez_stag, macro_cr, i, j, k, scomp); + amrex::Real const epsilon_interp = CoarsenIO::Interp( eps_arr, epsilon_stag, + Ez_stag, macro_cr, i, j, k, scomp); + amrex::Real beta = T_MacroAlgo::beta( sigma_interp, epsilon_interp, dt); + + push_ez_pml_current(i, j, k, Ez, Jz, + sigmaj_x, sigmaj_y, x_lo, y_lo, beta); + } + ); + } + + } + +} + +#endif // corresponds to ifndef WARPX_DIM_RZ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/Make.package b/Source/FieldSolver/FiniteDifferenceSolver/Make.package index 62f00ffc3a6..d269911a527 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/Make.package +++ b/Source/FieldSolver/FiniteDifferenceSolver/Make.package @@ -14,6 +14,7 @@ CEXE_sources += EvolveHPML.cpp CEXE_sources += EvolveBPML.cpp CEXE_sources += EvolveEPML.cpp CEXE_sources += EvolveFPML.cpp +CEXE_sources += MacroscopicEvolveEPML.cpp include $(WARPX_HOME)/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/Make.package From 3f5e7b02de018308f4c5ddeabc97a164bf844892 Mon Sep 17 00:00:00 2001 From: RevathiJambunathan Date: Mon, 7 Dec 2020 17:59:05 -0800 Subject: [PATCH 05/30] initialize macro PML and call macro PML --- .../FiniteDifferenceSolver.H | 36 +++++++++++++++++++ Source/FieldSolver/WarpXPushFieldsEM.cpp | 28 +++++++++------ 2 files changed, 54 insertions(+), 10 deletions(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index 48bf5bbbaef..c205d268b7a 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -128,12 +128,30 @@ class FiniteDifferenceSolver std::array< amrex::MultiFab*, 3 > const Efield, amrex::Real const dt ); + void MacroscopicEvolveEPML ( std::array< amrex::MultiFab*, 3 > Efield, +#ifndef WARPX_MAG_LLG + std::array< amrex::MultiFab*, 3 > const Bfield, +#else + std::array< amrex::MultiFab*, 3 > const Hfield, +#endif + std::array< amrex::MultiFab*, 3 > const Jfield, + amrex::MultiFab* const Ffield, + MultiSigmaBox const& sigba, + amrex::Real const dt, bool pml_has_particles, + std::unique_ptr const& macroscopic_properties, + amrex::MultiFab* const eps_mf, + amrex::MultiFab* const mu_mf, + amrex::MultiFab* const sigma_mf); + + #ifdef WARPX_MAG_LLG void EvolveHPML ( std::array< amrex::MultiFab*, 3 > Hfield, std::array< amrex::MultiFab*, 3 > const Efield, amrex::Real const dt ); #endif + + private: int m_fdtd_algo; @@ -264,6 +282,24 @@ class FiniteDifferenceSolver std::array< amrex::MultiFab*, 3 > const Efield, amrex::Real const dt ); + template< typename T_Algo, typename T_MacroAlgo > + void MacroscopicEvolveEPMLCartesian ( + std::array< amrex::MultiFab*, 3 > Efield, +#ifndef WARPX_MAG_LLG + std::array< amrex::MultiFab*, 3 > const Bfield, +#else + std::array< amrex::MultiFab*, 3 > const Hfield, +#endif + std::array< amrex::MultiFab*, 3 > const Jfield, + amrex::MultiFab* const Ffield, + MultiSigmaBox const& sigba, + amrex::Real const dt, bool pml_has_particles, + std::unique_ptr const& macroscopic_properties, + amrex::MultiFab* const eps_mf, + amrex::MultiFab* const mu_mf, + amrex::MultiFab* const sigma_mf ); + + #ifdef WARPX_MAG_LLG template< typename T_Algo > void EvolveHPMLCartesian ( diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index f1415f6812b..370893f44a3 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -346,27 +346,35 @@ WarpX::MacroscopicEvolveE (int lev, PatchType patch_type, amrex::Real a_dt) { // Evolve E field in PML cells if (do_pml && pml[lev]->ok()) { if (patch_type == PatchType::fine) { - m_fdtd_solver_fp[lev]->EvolveEPML( + m_fdtd_solver_fp[lev]->MacroscopicEvolveEPML( pml[lev]->GetE_fp(), -#ifdef WARPX_MAG_LLG - pml[lev]->GetH_fp(), -#else +#ifndef WARPX_MAG_LLG pml[lev]->GetB_fp(), +#else + pml[lev]->GetH_fp(), #endif pml[lev]->Getj_fp(), pml[lev]->GetF_fp(), pml[lev]->GetMultiSigmaBox_fp(), - a_dt, pml_has_particles ); + a_dt, pml_has_particles, + m_macroscopic_properties, + pml[lev]->Geteps_fp(), + pml[lev]->Getmu_fp(), + pml[lev]->Getsigma_fp() ); } else { - m_fdtd_solver_cp[lev]->EvolveEPML( + m_fdtd_solver_cp[lev]->MacroscopicEvolveEPML( pml[lev]->GetE_cp(), -#ifdef WARPX_MAG_LLG - pml[lev]->GetH_cp(), -#else +#ifndef WARPX_MAG_LLG pml[lev]->GetB_cp(), +#else + pml[lev]->GetH_cp(), #endif pml[lev]->Getj_cp(), pml[lev]->GetF_cp(), pml[lev]->GetMultiSigmaBox_cp(), - a_dt, pml_has_particles ); + a_dt, pml_has_particles, + m_macroscopic_properties, + pml[lev]->Geteps_cp(), + pml[lev]->Getmu_cp(), + pml[lev]->Getsigma_cp() ); } } } From 86c452fa735ae18a1d756951fffbd8c4c335006d Mon Sep 17 00:00:00 2001 From: Revathi Jambunathan Date: Mon, 7 Dec 2020 19:15:51 -0800 Subject: [PATCH 06/30] fix typo --- Source/BoundaryConditions/PML.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index daee9760bbf..91bfe82bc32 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -546,7 +546,7 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& /*grid_dm*/, // Initialize mu, permeability if (macroscopic_properties->m_mu_s == "constant") { pml_mu_fp->setVal(macroscopic_properties->m_mu); - } else if (macroscopic_properties->m_sigma_s == "parse_mu_function") { + } else if (macroscopic_properties->m_mu_s == "parse_mu_function") { macroscopic_properties->InitializeMacroMultiFabUsingParser(pml_mu_fp.get(), getParser(macroscopic_properties->m_mu_parser), lev); } From c786179d63ff8f2f83fa032972ac0493c938f377 Mon Sep 17 00:00:00 2001 From: Revathi Jambunathan Date: Mon, 7 Dec 2020 19:18:03 -0800 Subject: [PATCH 07/30] add relevant header --- .../FiniteDifferenceSolver/MacroscopicEvolveEPML.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveEPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveEPML.cpp index b3a6db30bb2..0f7762b8778 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveEPML.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveEPML.cpp @@ -5,7 +5,7 @@ #else # include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H" # include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H" -# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H" +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/FieldAccessorFunctors.H" #endif #include "BoundaryConditions/PML.H" #include "BoundaryConditions/PML_current.H" From 0946ea083a4b1b0bb57e2c728ac49c46cb6a7e02 Mon Sep 17 00:00:00 2001 From: Revathi Jambunathan Date: Mon, 7 Dec 2020 21:06:06 -0800 Subject: [PATCH 08/30] fix eol whitespace --- Source/BoundaryConditions/PML.cpp | 8 ++++---- .../FiniteDifferenceSolver/FiniteDifferenceSolver.H | 2 +- .../FiniteDifferenceSolver/MacroscopicEvolveEPML.cpp | 2 +- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 91bfe82bc32..ebaa0320bb0 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -516,17 +516,17 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& /*grid_dm*/, pml_H_fp[2] = std::make_unique(amrex::convert( ba, WarpX::GetInstance().getHfield_fp(0,2).ixType().toIntVect() ), dm, 2, ngb ); #endif - + if (WarpX::em_solver_medium == MediumForEM::Macroscopic) { // Allocating macroproperties in pml at cell-centers pml_eps_fp = std::make_unique(ba, dm, 1, nge); pml_mu_fp = std::make_unique(ba, dm, 1, nge); pml_sigma_fp = std::make_unique(ba, dm, 1, nge); - + // Initializing macroparameter multifab // auto& warpx = WarpX::GetInstance(); auto& macroscopic_properties = warpx.m_macroscopic_properties; - + // Initialize sigma, conductivity if (macroscopic_properties->m_sigma_s == "constant") { pml_sigma_fp->setVal(macroscopic_properties->m_sigma); @@ -650,7 +650,7 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& /*grid_dm*/, // Initializing macroparameter multifab // auto& warpx = WarpX::GetInstance(); auto& macroscopic_properties = warpx.m_macroscopic_properties; - + // Initialize sigma, conductivity if (macroscopic_properties->m_sigma_s == "constant") { pml_sigma_cp->setVal(macroscopic_properties->m_sigma); diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index c205d268b7a..4b36ed1b2cf 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -150,7 +150,7 @@ class FiniteDifferenceSolver amrex::Real const dt ); #endif - + private: diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveEPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveEPML.cpp index 0f7762b8778..f34c3213df5 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveEPML.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveEPML.cpp @@ -157,7 +157,7 @@ void FiniteDifferenceSolver::MacroscopicEvolveEPMLCartesian ( Array4 const& sigma_arr = sigma_mf->array(mfi); Array4 const& eps_arr = eps_mf->array(mfi); Array4 const& mu_arr = mu_mf->array(mfi); - + // Extract stencil coefficients Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr(); From 397f387f9278eb3b7b71ca6865a61a70fd97fbde Mon Sep 17 00:00:00 2001 From: Revathi Jambunathan <41089244+RevathiJambunathan@users.noreply.github.com> Date: Tue, 8 Dec 2020 09:46:39 -0800 Subject: [PATCH 09/30] Apply suggestions from code review Co-authored-by: Andy Nonaka --- Source/BoundaryConditions/PML.cpp | 1 - .../FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H | 2 -- 2 files changed, 3 deletions(-) diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index ebaa0320bb0..620609a43b0 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -11,7 +11,6 @@ #include "Utils/WarpXConst.H" #include "Utils/WarpXAlgorithmSelection.H" #include "WarpX.H" -//#include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H" #include #include diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index 4b36ed1b2cf..c16fc2eab67 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -150,8 +150,6 @@ class FiniteDifferenceSolver amrex::Real const dt ); #endif - - private: int m_fdtd_algo; From 01510f5d7cbff64e5eb5cfbaee39bcd20971c6f7 Mon Sep 17 00:00:00 2001 From: Revathi Jambunathan <41089244+RevathiJambunathan@users.noreply.github.com> Date: Tue, 8 Dec 2020 09:56:35 -0800 Subject: [PATCH 10/30] Update Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveEPML.cpp Co-authored-by: Andy Nonaka --- .../FiniteDifferenceSolver/MacroscopicEvolveEPML.cpp | 6 ------ 1 file changed, 6 deletions(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveEPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveEPML.cpp index f34c3213df5..6a5336e6dcf 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveEPML.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveEPML.cpp @@ -122,12 +122,6 @@ void FiniteDifferenceSolver::MacroscopicEvolveEPMLCartesian ( amrex::MultiFab* const mu_mf, amrex::MultiFab* const sigma_mf ) { - Real c2 = PhysConst::c * PhysConst::c; - -#ifdef WARPX_MAG_LLG - c2 *= PhysConst::mu0; -#endif - // Index type required for calling CoarsenIO::Interp to interpolate macroscopic // properties from their respective staggering to the Ex, Ey, Ez locations amrex::GpuArray const& sigma_stag = macroscopic_properties->sigma_IndexType; From c852c398ee5f055e2b8eba7cd2820de9b875547f Mon Sep 17 00:00:00 2001 From: Jackie Yao Date: Fri, 11 Dec 2020 09:17:21 -0800 Subject: [PATCH 11/30] added macroscopic mu in MacroscopicEvolveHPML --- .../Macroscopic_Maxwell/inputs_3d_LLG_noMs | 18 +- .../Tests/Macroscopic_Maxwell/inputs_3d_noMs | 21 +-- .../FiniteDifferenceSolver/CMakeLists.txt | 1 + .../FiniteDifferenceSolver.H | 16 +- .../MacroscopicEvolveHPML.cpp | 167 ++++++++++++++++++ .../MacroscopicProperties.H | 6 + .../FiniteDifferenceSolver/Make.package | 1 + Source/FieldSolver/WarpXPushFieldsEM.cpp | 32 +++- 8 files changed, 234 insertions(+), 28 deletions(-) create mode 100644 Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp diff --git a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_noMs b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_noMs index 58bb683488a..924e4ecfe54 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,22 +20,23 @@ 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 = 4 +my_constants.epsilon_r = 1 my_constants.pi = 3.14159265359 my_constants.L = 141.4213562373095e-6 my_constants.c = 299792458. @@ -45,10 +45,10 @@ 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)" +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.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.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. @@ -65,4 +65,4 @@ macroscopic.mag_gamma = 0. diagnostics.diags_names = plt plt.period = 10 plt.fields_to_plot = Ex Ey Ez Hx Hy Hz Bx By Bz -plt.diag_type = Full +plt.diag_type = Full \ No newline at end of file diff --git a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs index d9393c00681..7b8ab242b2f 100644 --- a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs +++ b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs @@ -1,15 +1,14 @@ -# This is a E+B (use USE_LLG=FALSE, ie., no M field) simulation of a +# This is a E+H (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 geometry.coord_sys = 0 -geometry.is_periodic = 1 1 0 +geometry.is_periodic = 1 1 1 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 @@ -21,22 +20,23 @@ 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. @@ -45,10 +45,11 @@ 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)" +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.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.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. @@ -56,4 +57,4 @@ warpx.Bz_external_grid_function(x,y,z) = 0. diagnostics.diags_names = plt plt.period = 10 plt.fields_to_plot = Ex Ey Ez Bx By Bz -plt.diag_type = Full +plt.diag_type = Full \ No newline at end of file diff --git a/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt b/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt index a792fbbfa14..67bdeb85835 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt +++ b/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt @@ -20,6 +20,7 @@ if(WarpX_MAG_LLG) 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 c16fc2eab67..e77c9a53f4f 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -148,6 +148,12 @@ class FiniteDifferenceSolver void EvolveHPML ( std::array< amrex::MultiFab*, 3 > Hfield, std::array< amrex::MultiFab*, 3 > const Efield, amrex::Real const dt ); + + 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 ); #endif private: @@ -301,9 +307,17 @@ class FiniteDifferenceSolver #ifdef WARPX_MAG_LLG template< typename T_Algo > void EvolveHPMLCartesian ( - std::array< amrex::MultiFab*, 3 > Bfield, + std::array< amrex::MultiFab*, 3 > Hfield, std::array< amrex::MultiFab*, 3 > const Efield, amrex::Real const dt ); + + template< typename T_Algo > + 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 ); #endif #endif diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp new file mode 100644 index 00000000000..8ad35084a56 --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp @@ -0,0 +1,167 @@ +/* Copyright 2020 Remi Lehe + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#include "Utils/WarpXAlgorithmSelection.H" +#include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H" +#ifdef WARPX_DIM_RZ +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H" +#else +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H" +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H" +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H" +#endif +#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 + */ +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 ) { + + // 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, macroscopic_properties, mu_mf); + amrex::Abort("PML are not implemented in cylindrical geometry."); +#else + if (m_do_nodal) { + + MacroscopicEvolveHPMLCartesian ( Hfield, Efield, dt, macroscopic_properties, mu_mf ); + + } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) { + + MacroscopicEvolveHPMLCartesian ( Hfield, Efield, dt, macroscopic_properties, mu_mf ); + + } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) { + + MacroscopicEvolveHPMLCartesian ( Hfield, Efield, dt, macroscopic_properties, mu_mf ); + + } else { + amrex::Abort("Unknown algorithm"); + } +#endif +} + + +#ifndef WARPX_DIM_RZ + +template +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 ) { + + 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()) +#endif + + for ( MFIter mfi(*Hfield[0], TilingIfNotGPU()); mfi.isValid(); ++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& Ex = Efield[0]->array(mfi); + Array4 const& Ey = Efield[1]->array(mfi); + Array4 const& Ez = Efield[2]->array(mfi); + + // Extract stencil coefficients + Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr(); + int const n_coefs_x = m_stencil_coefs_x.size(); + Real const * const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr(); + int const n_coefs_y = m_stencil_coefs_y.size(); + 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 + Box const& tbx = mfi.tilebox(Hfield[0]->ixType().ixType()); + Box const& tby = mfi.tilebox(Hfield[1]->ixType().ixType()); + Box const& tbz = mfi.tilebox(Hfield[2]->ixType().ixType()); + + // 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_arrx = CoarsenIO::Interp( mu_arr, mu_stag, Hx_stag, macro_cr, i, j, k, 0); + + if (i ==0 && j == 0 && k == 0){ + std::cout << "i:" << i << "\n j:" << j << "\n k:" << k << std::endl; + std::cout << "mux: " << mu_arrx << std::endl; + std::cout << "Hxz: " << Hx(i, j, k, PMLComp::xz) << std::endl; + std::cout << "Hxy: " << Hx(i, j, k, PMLComp::xy) << std::endl;} + + Hx(i, j, k, PMLComp::xz) += 1._rt / mu_arrx * 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::yy) + + T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, PMLComp::yz) ); + Hx(i, j, k, PMLComp::xy) -= 1._rt / mu_arrx * 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) + + T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k, PMLComp::zz) ); + }, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + + Real mu_arry = CoarsenIO::Interp( mu_arr, mu_stag, Hy_stag, macro_cr, i, j, k, 0); + + Hy(i, j, k, PMLComp::yx) += 1._rt / mu_arry * 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) + + T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k, PMLComp::zz) ); + Hy(i, j, k, PMLComp::yz) -= 1._rt / mu_arry * dt * ( + T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k, PMLComp::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) ); + }, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + + Real mu_arrz = CoarsenIO::Interp( mu_arr, mu_stag, Hz_stag, macro_cr, i, j, k, 0); + + Hz(i, j, k, PMLComp::zy) += 1._rt / mu_arrz * dt * ( + T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k, PMLComp::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) -= 1._rt / mu_arrz * 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::yy) + + T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k, PMLComp::yz) ); + } + + ); + + } + +} + +#endif // corresponds to ifndef WARPX_DIM_RZ + +#endif // #ifdef WARPX_MAG_LLG 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/Make.package b/Source/FieldSolver/FiniteDifferenceSolver/Make.package index d269911a527..65424badbbb 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/Make.package +++ b/Source/FieldSolver/FiniteDifferenceSolver/Make.package @@ -9,6 +9,7 @@ CEXE_sources += MacroscopicEvolveE.cpp 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 370893f44a3..7015ed0d40e 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -414,11 +414,19 @@ 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 ); + m_fdtd_solver_fp[lev]->MacroscopicEvolveHPML( + pml[lev]->GetH_fp(), + pml[lev]->GetE_fp(), + a_dt, + m_macroscopic_properties, + pml[lev]->Getmu_fp() ); } else { - m_fdtd_solver_cp[lev]->EvolveHPML( - pml[lev]->GetH_cp(), pml[lev]->GetE_cp(), a_dt ); + m_fdtd_solver_cp[lev]->MacroscopicEvolveHPML( + pml[lev]->GetH_cp(), + pml[lev]->GetE_cp(), + a_dt, + m_macroscopic_properties, + pml[lev]->Getmu_cp() ); } } } @@ -457,11 +465,19 @@ 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 ); + m_fdtd_solver_fp[lev]->MacroscopicEvolveHPML( + pml[lev]->GetH_fp(), + pml[lev]->GetE_fp(), + a_dt, + m_macroscopic_properties, + pml[lev]->Getmu_fp() ); } else { - m_fdtd_solver_cp[lev]->EvolveHPML( - pml[lev]->GetH_cp(), pml[lev]->GetE_cp(), a_dt ); + m_fdtd_solver_cp[lev]->MacroscopicEvolveHPML( + pml[lev]->GetH_cp(), + pml[lev]->GetE_cp(), + a_dt, + m_macroscopic_properties, + pml[lev]->Getmu_cp() ); } } } From c3af5e0d1ef7abdc30f6ea8b752754f5ce09ca50 Mon Sep 17 00:00:00 2001 From: Jackie Yao Date: Fri, 11 Dec 2020 09:23:55 -0800 Subject: [PATCH 12/30] removed printing residual in Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp --- .../FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp | 6 ------ 1 file changed, 6 deletions(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp index 8ad35084a56..b3d2bfdcc6c 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp @@ -111,12 +111,6 @@ void FiniteDifferenceSolver::MacroscopicEvolveHPMLCartesian ( [=] AMREX_GPU_DEVICE (int i, int j, int k){ Real mu_arrx = CoarsenIO::Interp( mu_arr, mu_stag, Hx_stag, macro_cr, i, j, k, 0); - - if (i ==0 && j == 0 && k == 0){ - std::cout << "i:" << i << "\n j:" << j << "\n k:" << k << std::endl; - std::cout << "mux: " << mu_arrx << std::endl; - std::cout << "Hxz: " << Hx(i, j, k, PMLComp::xz) << std::endl; - std::cout << "Hxy: " << Hx(i, j, k, PMLComp::xy) << std::endl;} Hx(i, j, k, PMLComp::xz) += 1._rt / mu_arrx * dt * ( T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, PMLComp::yx) From 7689f472ea8fc2f387a91a912bc5ef29013d03da Mon Sep 17 00:00:00 2001 From: Jackie Yao Date: Fri, 11 Dec 2020 13:49:25 -0800 Subject: [PATCH 13/30] updated ncomp in Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp --- .../FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp index b3d2bfdcc6c..aa81a0b40df 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp @@ -101,7 +101,8 @@ void FiniteDifferenceSolver::MacroscopicEvolveHPMLCartesian ( Box const& tbx = mfi.tilebox(Hfield[0]->ixType().ixType()); Box const& tby = mfi.tilebox(Hfield[1]->ixType().ixType()); Box const& tbz = mfi.tilebox(Hfield[2]->ixType().ixType()); - + // 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); @@ -110,7 +111,7 @@ void FiniteDifferenceSolver::MacroscopicEvolveHPMLCartesian ( [=] AMREX_GPU_DEVICE (int i, int j, int k){ - Real mu_arrx = CoarsenIO::Interp( mu_arr, mu_stag, Hx_stag, macro_cr, i, j, k, 0); + Real mu_arrx = CoarsenIO::Interp( mu_arr, mu_stag, Hx_stag, macro_cr, i, j, k, scomp); Hx(i, j, k, PMLComp::xz) += 1._rt / mu_arrx * dt * ( T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, PMLComp::yx) @@ -124,7 +125,7 @@ void FiniteDifferenceSolver::MacroscopicEvolveHPMLCartesian ( [=] AMREX_GPU_DEVICE (int i, int j, int k){ - Real mu_arry = CoarsenIO::Interp( mu_arr, mu_stag, Hy_stag, macro_cr, i, j, k, 0); + Real mu_arry = CoarsenIO::Interp( mu_arr, mu_stag, Hy_stag, macro_cr, i, j, k, scomp); Hy(i, j, k, PMLComp::yx) += 1._rt / mu_arry * dt * ( T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k, PMLComp::zx) @@ -138,7 +139,7 @@ void FiniteDifferenceSolver::MacroscopicEvolveHPMLCartesian ( [=] AMREX_GPU_DEVICE (int i, int j, int k){ - Real mu_arrz = CoarsenIO::Interp( mu_arr, mu_stag, Hz_stag, macro_cr, i, j, k, 0); + Real mu_arrz = CoarsenIO::Interp( mu_arr, mu_stag, Hz_stag, macro_cr, i, j, k, scomp); Hz(i, j, k, PMLComp::zy) += 1._rt / mu_arrz * dt * ( T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k, PMLComp::xx) From d479e63b582d07ee803ea64eee6a333e001c412e Mon Sep 17 00:00:00 2001 From: Jackie Yao Date: Fri, 11 Dec 2020 14:01:58 -0800 Subject: [PATCH 14/30] EOL fix --- .../MacroscopicEvolveHPML.cpp | 2 +- Source/FieldSolver/WarpXPushFieldsEM.cpp | 32 +++++++++---------- 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp index aa81a0b40df..bbd67f2e0e1 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp @@ -126,7 +126,7 @@ void FiniteDifferenceSolver::MacroscopicEvolveHPMLCartesian ( [=] AMREX_GPU_DEVICE (int i, int j, int k){ Real mu_arry = CoarsenIO::Interp( mu_arr, mu_stag, Hy_stag, macro_cr, i, j, k, scomp); - + Hy(i, j, k, PMLComp::yx) += 1._rt / mu_arry * 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) diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 7015ed0d40e..3d127efee0b 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -415,17 +415,17 @@ WarpX::MacroscopicEvolveHM (int lev, PatchType patch_type, amrex::Real a_dt) { if (do_pml && pml[lev]->ok()) { if (patch_type == PatchType::fine) { m_fdtd_solver_fp[lev]->MacroscopicEvolveHPML( - pml[lev]->GetH_fp(), - pml[lev]->GetE_fp(), - a_dt, - m_macroscopic_properties, + pml[lev]->GetH_fp(), + pml[lev]->GetE_fp(), + a_dt, + m_macroscopic_properties, pml[lev]->Getmu_fp() ); } else { m_fdtd_solver_cp[lev]->MacroscopicEvolveHPML( - pml[lev]->GetH_cp(), - pml[lev]->GetE_cp(), - a_dt, - m_macroscopic_properties, + pml[lev]->GetH_cp(), + pml[lev]->GetE_cp(), + a_dt, + m_macroscopic_properties, pml[lev]->Getmu_cp() ); } } @@ -466,17 +466,17 @@ WarpX::MacroscopicEvolveHM_2nd (int lev, PatchType patch_type, amrex::Real a_dt) if (do_pml && pml[lev]->ok()) { if (patch_type == PatchType::fine) { m_fdtd_solver_fp[lev]->MacroscopicEvolveHPML( - pml[lev]->GetH_fp(), - pml[lev]->GetE_fp(), - a_dt, - m_macroscopic_properties, + pml[lev]->GetH_fp(), + pml[lev]->GetE_fp(), + a_dt, + m_macroscopic_properties, pml[lev]->Getmu_fp() ); } else { m_fdtd_solver_cp[lev]->MacroscopicEvolveHPML( - pml[lev]->GetH_cp(), - pml[lev]->GetE_cp(), - a_dt, - m_macroscopic_properties, + pml[lev]->GetH_cp(), + pml[lev]->GetE_cp(), + a_dt, + m_macroscopic_properties, pml[lev]->Getmu_cp() ); } } From 7cab0d7800391ada5cb5f7eb02ee21c2bcb339b8 Mon Sep 17 00:00:00 2001 From: Jackie Yao Date: Fri, 11 Dec 2020 14:04:43 -0800 Subject: [PATCH 15/30] fixed pml setup in Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs --- Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs index 7b8ab242b2f..b1fc054a3d8 100644 --- a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs +++ b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs @@ -8,7 +8,7 @@ amr.n_cell = 16 16 128 amr.max_grid_size = 32 amr.blocking_factor = 16 geometry.coord_sys = 0 -geometry.is_periodic = 1 1 1 +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 From cd52411ec6f886b9537e97a850c43ee222c13f9e Mon Sep 17 00:00:00 2001 From: jackiezy <58234082+jackieyao0114@users.noreply.github.com> Date: Fri, 11 Dec 2020 14:06:15 -0800 Subject: [PATCH 16/30] Apply suggestions from code review (annotations in code) Co-authored-by: Revathi Jambunathan <41089244+RevathiJambunathan@users.noreply.github.com> --- .../FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp index bbd67f2e0e1..60202976be9 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.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" @@ -24,7 +18,7 @@ using namespace amrex; #ifdef WARPX_MAG_LLG /** - * \brief Update the H field, over one timestep + * \brief Update Hfield in PML region */ void FiniteDifferenceSolver::MacroscopicEvolveHPML ( std::array< amrex::MultiFab*, 3 > Hfield, From edac43ea6cbe3920eb33567a8f00c66d858a48c2 Mon Sep 17 00:00:00 2001 From: Jackie Yao Date: Fri, 11 Dec 2020 15:00:05 -0800 Subject: [PATCH 17/30] clarified mur and Ms requirement in Docs/source/running_cpp/parameters.rst --- Docs/source/running_cpp/parameters.rst | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index 970b5e9c10b..80262b0df6e 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -1216,7 +1216,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. @@ -1236,7 +1238,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 unitless 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`. @@ -2219,6 +2224,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 unitless 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 From 643f04d6b240f7ac33a2b7b312f31914bbc02b6f Mon Sep 17 00:00:00 2001 From: jackiezy <58234082+jackieyao0114@users.noreply.github.com> Date: Tue, 15 Dec 2020 13:55:10 -0500 Subject: [PATCH 18/30] Update Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs Co-authored-by: Andy Nonaka --- Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs index b1fc054a3d8..e463deef278 100644 --- a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs +++ b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs @@ -1,4 +1,4 @@ -# This is a E+H (use USE_LLG=FALSE) simulation of a +# This is a E+B (use USE_LLG=FALSE) simulation of a # plane wave in a periodic domain ################################ ####### GENERAL PARAMETERS ###### @@ -57,4 +57,4 @@ warpx.Bz_external_grid_function(x,y,z) = 0. diagnostics.diags_names = plt plt.period = 10 plt.fields_to_plot = Ex Ey Ez Bx By Bz -plt.diag_type = Full \ No newline at end of file +plt.diag_type = Full From 3034b3f355cad2e18e37a5c06bd05724c748b549 Mon Sep 17 00:00:00 2001 From: jackiezy <58234082+jackieyao0114@users.noreply.github.com> Date: Thu, 17 Dec 2020 11:32:17 -0500 Subject: [PATCH 19/30] Update Docs/source/running_cpp/parameters.rst Co-authored-by: Revathi Jambunathan <41089244+RevathiJambunathan@users.noreply.github.com> --- Docs/source/running_cpp/parameters.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index 80262b0df6e..7cb21c79669 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -1241,7 +1241,7 @@ Numerics and algorithms 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 unitless relative permeability ``mur >=1``. + 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`. From a9aba04be07cf6c40347bcb50347f4c399518ba9 Mon Sep 17 00:00:00 2001 From: jackiezy <58234082+jackieyao0114@users.noreply.github.com> Date: Thu, 17 Dec 2020 14:49:00 -0500 Subject: [PATCH 20/30] removed part of EvolveHPML residual Co-authored-by: Revathi Jambunathan <41089244+RevathiJambunathan@users.noreply.github.com> --- Docs/source/running_cpp/parameters.rst | 2 +- Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt | 1 - Source/FieldSolver/FiniteDifferenceSolver/Make.package | 1 - 3 files changed, 1 insertion(+), 3 deletions(-) diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index 7cb21c79669..9570d043fce 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -2226,7 +2226,7 @@ Solving magnetization using LLG equation 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 unitless relative permeability ``mur >=1``. + 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/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt b/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt index 67bdeb85835..34f01cd8e03 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt +++ b/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt @@ -19,7 +19,6 @@ if(WarpX_MAG_LLG) MacroscopicEvolveHM.cpp MacroscopicEvolveM_2nd.cpp MacroscopicEvolveM.cpp - EvolveHPML.cpp MacroscopicEvolveHPML.cpp ) endif() diff --git a/Source/FieldSolver/FiniteDifferenceSolver/Make.package b/Source/FieldSolver/FiniteDifferenceSolver/Make.package index 65424badbbb..c20f4531c72 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/Make.package +++ b/Source/FieldSolver/FiniteDifferenceSolver/Make.package @@ -8,7 +8,6 @@ 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 From 993e8342d25bf44c0ff047d47f518704f560be24 Mon Sep 17 00:00:00 2001 From: Jackie Yao Date: Thu, 17 Dec 2020 12:35:45 -0800 Subject: [PATCH 21/30] removed more EvolveHPML residual --- .../FiniteDifferenceSolver/FiniteDifferenceSolver.H | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index e77c9a53f4f..814639e0c72 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -145,10 +145,6 @@ class FiniteDifferenceSolver #ifdef WARPX_MAG_LLG - void EvolveHPML ( std::array< amrex::MultiFab*, 3 > Hfield, - std::array< amrex::MultiFab*, 3 > const Efield, - amrex::Real const dt ); - void MacroscopicEvolveHPML ( std::array< amrex::MultiFab*, 3 > Hfield, std::array< amrex::MultiFab*, 3 > const Efield, amrex::Real const dt, @@ -305,12 +301,6 @@ class FiniteDifferenceSolver #ifdef WARPX_MAG_LLG - template< typename T_Algo > - void EvolveHPMLCartesian ( - std::array< amrex::MultiFab*, 3 > Hfield, - std::array< amrex::MultiFab*, 3 > const Efield, - amrex::Real const dt ); - template< typename T_Algo > void MacroscopicEvolveHPMLCartesian ( std::array< amrex::MultiFab*, 3 > Hfield, From 61fe0dc4d4c75e34402a25bc4af41139d7555d57 Mon Sep 17 00:00:00 2001 From: Jackie Yao Date: Thu, 17 Dec 2020 12:37:50 -0800 Subject: [PATCH 22/30] removed Source/FieldSolver/FiniteDifferenceSolver/EvolveHPML.cpp --- .../FiniteDifferenceSolver/EvolveHPML.cpp | 139 ------------------ 1 file changed, 139 deletions(-) delete mode 100644 Source/FieldSolver/FiniteDifferenceSolver/EvolveHPML.cpp diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveHPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveHPML.cpp deleted file mode 100644 index 470dce0beba..00000000000 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveHPML.cpp +++ /dev/null @@ -1,139 +0,0 @@ -/* Copyright 2020 Remi Lehe - * - * This file is part of WarpX. - * - * License: BSD-3-Clause-LBNL - */ - -#include "Utils/WarpXAlgorithmSelection.H" -#include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H" -#ifdef WARPX_DIM_RZ -# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H" -#else -# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H" -# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H" -# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H" -#endif -#include "BoundaryConditions/PMLComponent.H" -#include -#include - -using namespace amrex; - -#ifdef WARPX_MAG_LLG - -/** - * \brief Update the H field, over one timestep - */ -void FiniteDifferenceSolver::EvolveHPML ( - std::array< amrex::MultiFab*, 3 > Hfield, - std::array< amrex::MultiFab*, 3 > const Efield, - amrex::Real const dt ) { - - // 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::Abort("PML are not implemented in cylindrical geometry."); -#else - if (m_do_nodal) { - - EvolveHPMLCartesian ( Hfield, Efield, dt ); - - } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) { - - EvolveHPMLCartesian ( Hfield, Efield, dt ); - - } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) { - - EvolveHPMLCartesian ( Hfield, Efield, dt ); - - } else { - amrex::Abort("Unknown algorithm"); - } -#endif -} - - -#ifndef WARPX_DIM_RZ - -template -void FiniteDifferenceSolver::EvolveHPMLCartesian ( - std::array< amrex::MultiFab*, 3 > Hfield, - std::array< amrex::MultiFab*, 3 > const Efield, - amrex::Real const dt ) { - - // Loop through the grids, and over the tiles within each grid -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(*Hfield[0], TilingIfNotGPU()); mfi.isValid(); ++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& Ex = Efield[0]->array(mfi); - Array4 const& Ey = Efield[1]->array(mfi); - Array4 const& Ez = Efield[2]->array(mfi); - - // Extract stencil coefficients - Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr(); - int const n_coefs_x = m_stencil_coefs_x.size(); - Real const * const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr(); - int const n_coefs_y = m_stencil_coefs_y.size(); - 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 - Box const& tbx = mfi.tilebox(Hfield[0]->ixType().ixType()); - 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; - - // Loop over the cells and update the fields - amrex::ParallelFor(tbx, tby, tbz, - - [=] AMREX_GPU_DEVICE (int i, int j, int k){ - Hx(i, j, k, PMLComp::xz) += mu0_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::yy) - + T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, PMLComp::yz) ); - Hx(i, j, k, PMLComp::xy) -= mu0_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) - + T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k, PMLComp::zz) ); - }, - - [=] AMREX_GPU_DEVICE (int i, int j, int k){ - Hy(i, j, k, PMLComp::yx) += mu0_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) - + T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k, PMLComp::zz) ); - Hy(i, j, k, PMLComp::yz) -= mu0_inv * dt * ( - T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k, PMLComp::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) ); - }, - - [=] AMREX_GPU_DEVICE (int i, int j, int k){ - Hz(i, j, k, PMLComp::zy) += mu0_inv * dt * ( - T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k, PMLComp::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 * ( - 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::yy) - + T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k, PMLComp::yz) ); - } - - ); - - } - -} - -#endif // corresponds to ifndef WARPX_DIM_RZ - -#endif // #ifdef WARPX_MAG_LLG From 117a3a0f14024d35e0ce469ed03f2672c2c2b541 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Thu, 17 Dec 2020 19:31:04 -0800 Subject: [PATCH 23/30] initialize H_indexType --- .../MacroscopicProperties/MacroscopicProperties.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp index a74f9fa268b..ad201a3a8df 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp @@ -235,6 +235,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]; @@ -250,6 +253,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; } @@ -267,6 +273,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 From 8f9716fa0020fa31dce6ba4fc046ac18300af36e Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Mon, 21 Dec 2020 10:58:32 -0800 Subject: [PATCH 24/30] Update Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp --- .../FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp index 60202976be9..173ce663a31 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp @@ -72,7 +72,6 @@ void FiniteDifferenceSolver::MacroscopicEvolveHPMLCartesian ( #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif - for ( MFIter mfi(*Hfield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { // Extract field data for this grid/tile From 537495b78e97443a332da5c4549e6d879b30d153 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Mon, 21 Dec 2020 11:15:36 -0800 Subject: [PATCH 25/30] cleanup only --- .../MacroscopicEvolveHPML.cpp | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp index 173ce663a31..8e9c7095683 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp @@ -94,6 +94,7 @@ void FiniteDifferenceSolver::MacroscopicEvolveHPMLCartesian ( Box const& tbx = mfi.tilebox(Hfield[0]->ixType().ixType()); Box const& tby = mfi.tilebox(Hfield[1]->ixType().ixType()); Box const& tbz = mfi.tilebox(Hfield[2]->ixType().ixType()); + // 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 @@ -104,13 +105,13 @@ void FiniteDifferenceSolver::MacroscopicEvolveHPMLCartesian ( [=] AMREX_GPU_DEVICE (int i, int j, int k){ - Real mu_arrx = CoarsenIO::Interp( mu_arr, mu_stag, Hx_stag, macro_cr, i, j, k, scomp); + Real mu_inv = 1._rt/CoarsenIO::Interp( mu_arr, mu_stag, Hx_stag, macro_cr, i, j, k, scomp); - Hx(i, j, k, PMLComp::xz) += 1._rt / mu_arrx * 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::yy) + T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, PMLComp::yz) ); - Hx(i, j, k, PMLComp::xy) -= 1._rt / mu_arrx * 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) + T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k, PMLComp::zz) ); @@ -118,13 +119,13 @@ void FiniteDifferenceSolver::MacroscopicEvolveHPMLCartesian ( [=] AMREX_GPU_DEVICE (int i, int j, int k){ - Real mu_arry = CoarsenIO::Interp( mu_arr, mu_stag, Hy_stag, macro_cr, i, j, k, scomp); + Real mu_inv = 1._rt/CoarsenIO::Interp( mu_arr, mu_stag, Hy_stag, macro_cr, i, j, k, scomp); - Hy(i, j, k, PMLComp::yx) += 1._rt / mu_arry * 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) + T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k, PMLComp::zz) ); - Hy(i, j, k, PMLComp::yz) -= 1._rt / mu_arry * dt * ( + Hy(i, j, k, PMLComp::yz) -= mu_inv * dt * ( T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k, PMLComp::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) ); @@ -132,13 +133,13 @@ void FiniteDifferenceSolver::MacroscopicEvolveHPMLCartesian ( [=] AMREX_GPU_DEVICE (int i, int j, int k){ - Real mu_arrz = CoarsenIO::Interp( mu_arr, mu_stag, Hz_stag, macro_cr, i, j, k, scomp); + Real mu_inv = 1._rt/CoarsenIO::Interp( mu_arr, mu_stag, Hz_stag, macro_cr, i, j, k, scomp); - Hz(i, j, k, PMLComp::zy) += 1._rt / mu_arrz * dt * ( + Hz(i, j, k, PMLComp::zy) += mu_inv * dt * ( T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k, PMLComp::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) -= 1._rt / mu_arrz * 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::yy) + T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k, PMLComp::yz) ); From 13ad49e73128a76c3ac139d5c7e15f198f488b35 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Mon, 21 Dec 2020 11:45:05 -0800 Subject: [PATCH 26/30] clean up inputs files (no functional changes) --- .../Tests/Macroscopic_Maxwell/inputs_3d_LLG_noMs | 12 ++++++------ Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs | 7 +++---- 2 files changed, 9 insertions(+), 10 deletions(-) diff --git a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_noMs b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_noMs index 924e4ecfe54..a1923861a8d 100644 --- a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_noMs +++ b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_noMs @@ -35,7 +35,7 @@ particles.nspecies = 0 ################################# ############ FIELDS ############# ################################# -my_constants.mu_r = 4 +my_constants.mu_r = 1 my_constants.epsilon_r = 1 my_constants.pi = 3.14159265359 my_constants.L = 141.4213562373095e-6 @@ -43,13 +43,13 @@ 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 * 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.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 @@ -65,4 +65,4 @@ macroscopic.mag_gamma = 0. diagnostics.diags_names = plt plt.period = 10 plt.fields_to_plot = Ex Ey Ez Hx Hy Hz Bx By Bz -plt.diag_type = Full \ No newline at end of file +plt.diag_type = Full diff --git a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs index e463deef278..35e881e1b2b 100644 --- a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs +++ b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs @@ -43,14 +43,13 @@ 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 * sqrt(epsilon_r*mu_r))/c * sqrt(epsilon_r*mu_r)" -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 From 94e089fe35e5e20e7bba2828b06ad0fb53f58663 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Mon, 21 Dec 2020 11:53:13 -0800 Subject: [PATCH 27/30] add Ms at boundary test. sync up inputs files --- .../inputs_3d_LLG_MsBoundary | 71 +++++++++++++++++++ .../Macroscopic_Maxwell/inputs_3d_LLG_noMs | 2 - .../Tests/Macroscopic_Maxwell/inputs_3d_noMs | 2 - 3 files changed, 71 insertions(+), 4 deletions(-) create mode 100644 Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary 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..5d509c678d4 --- /dev/null +++ b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary @@ -0,0 +1,71 @@ +# 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 = 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.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 * 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 * 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 = constant +macroscopic.mag_alpha = 0. +macroscopic.mag_gamma_init_style = constant +macroscopic.mag_gamma = 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.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 a1923861a8d..60c7ba3884e 100644 --- a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_noMs +++ b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_noMs @@ -25,10 +25,8 @@ 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 diff --git a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs index 35e881e1b2b..32cf5ef0e8e 100644 --- a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs +++ b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs @@ -25,10 +25,8 @@ 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 From f1ce428cb3ac4068a8b9b77a5676d7c5afd5d10b Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Mon, 21 Dec 2020 12:10:37 -0800 Subject: [PATCH 28/30] set alpha and gamma --- .../Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary index 5d509c678d4..8b725d996c1 100644 --- a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary +++ b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary @@ -54,10 +54,12 @@ warpx.Hz_external_grid_function(x,y,z) = 0. 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 = constant -macroscopic.mag_alpha = 0. -macroscopic.mag_gamma_init_style = constant -macroscopic.mag_gamma = 0. + +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. From fa242357e472753a1ece760b7e22ca3805addf25 Mon Sep 17 00:00:00 2001 From: Jackie Yao Date: Tue, 5 Jan 2021 16:25:12 -0800 Subject: [PATCH 29/30] added Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary --- Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary index 8b725d996c1..b8b667f45df 100644 --- a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary +++ b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary @@ -33,7 +33,7 @@ particles.nspecies = 0 ################################# ############ FIELDS ############# ################################# -my_constants.mu_r = 1 +my_constants.mu_r = 4 my_constants.epsilon_r = 1 my_constants.pi = 3.14159265359 my_constants.L = 141.4213562373095e-6 @@ -42,11 +42,11 @@ 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.e5*exp(-z**2/L**2)*cos(2*pi*z/wavelength * sqrt(epsilon_r*mu_r))" +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.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.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. From 461796ecdf442e67f81f46e0432694833c726ca7 Mon Sep 17 00:00:00 2001 From: Jackie Yao Date: Tue, 12 Jan 2021 08:24:30 -0800 Subject: [PATCH 30/30] added Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary --- Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary index 8b725d996c1..b8b667f45df 100644 --- a/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary +++ b/Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_MsBoundary @@ -33,7 +33,7 @@ particles.nspecies = 0 ################################# ############ FIELDS ############# ################################# -my_constants.mu_r = 1 +my_constants.mu_r = 4 my_constants.epsilon_r = 1 my_constants.pi = 3.14159265359 my_constants.L = 141.4213562373095e-6 @@ -42,11 +42,11 @@ 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.e5*exp(-z**2/L**2)*cos(2*pi*z/wavelength * sqrt(epsilon_r*mu_r))" +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.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.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.