aboutsummaryrefslogtreecommitdiff
path: root/Source/WarpXEvolve.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'Source/WarpXEvolve.cpp')
-rw-r--r--Source/WarpXEvolve.cpp145
1 files changed, 98 insertions, 47 deletions
diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp
index a92f033b5..b4e1c8e95 100644
--- a/Source/WarpXEvolve.cpp
+++ b/Source/WarpXEvolve.cpp
@@ -102,9 +102,6 @@ WarpX::EvolveEM (int numsteps)
SyncCurrent();
SyncRho(rho_fp, rho_cp);
-#ifdef WARPX_USE_PSATD
- SyncRho(rho2_fp, rho2_cp);
-#endif
// Push E and B from {n} to {n+1}
// (And update guard cells immediately afterwards)
@@ -113,12 +110,17 @@ WarpX::EvolveEM (int numsteps)
FillBoundaryE();
FillBoundaryB();
#else
- EvolveF(dt[0], DtType::Full);
- EvolveB(0.5*dt[0], DtType::SecondHalf); // We now have B^{n+1/2}
+ EvolveF(0.5*dt[0], DtType::FirstHalf);
+ EvolveB(0.5*dt[0]); // We now have B^{n+1/2}
FillBoundaryB();
- EvolveE(dt[0], DtType::Full); // We now have E^{n+1}
+ EvolveE(dt[0]); // We now have E^{n+1}
FillBoundaryE();
- EvolveB(0.5*dt[0], DtType::FirstHalf); // We now have B^{n+1}
+ EvolveF(0.5*dt[0], DtType::SecondHalf);
+ EvolveB(0.5*dt[0]); // We now have B^{n+1}
+ if (do_pml) {
+ DampPML();
+ FillBoundaryE();
+ }
FillBoundaryB();
#endif
@@ -222,15 +224,15 @@ WarpX::EvolveEM (int numsteps)
}
void
-WarpX::EvolveB (Real dt, DtType typ)
+WarpX::EvolveB (Real dt)
{
for (int lev = 0; lev <= finest_level; ++lev) {
- EvolveB(lev, dt, typ);
+ EvolveB(lev, dt);
}
}
void
-WarpX::EvolveB (int lev, Real dt, DtType typ)
+WarpX::EvolveB (int lev, Real dt)
{
BL_PROFILE("WarpX::EvolveB()");
@@ -305,14 +307,10 @@ WarpX::EvolveB (int lev, Real dt, DtType typ)
if (do_pml && pml[lev]->ok())
{
- const int dttype = static_cast<int>(typ);
-
for (int ipatch = 0; ipatch < npatches; ++ipatch)
{
const auto& pml_B = (ipatch==0) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp();
const auto& pml_E = (ipatch==0) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
- const auto& sigba = (ipatch==0) ? pml[lev]->GetMultiSigmaBox_fp()
- : pml[lev]->GetMultiSigmaBox_cp();
int patch_level = (ipatch == 0) ? lev : lev-1;
const std::array<Real,3>& dx = WarpX::CellSize(patch_level);
const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]};
@@ -335,7 +333,6 @@ WarpX::EvolveB (int lev, Real dt, DtType typ)
BL_TO_FORTRAN_3D((*pml_B[0])[mfi]),
BL_TO_FORTRAN_3D((*pml_B[1])[mfi]),
BL_TO_FORTRAN_3D((*pml_B[2])[mfi]),
- WRPX_PML_SIGMA_STAR_TO_FORTRAN(sigba[mfi],dttype),
&dtsdx[0], &dtsdx[1], &dtsdx[2],
&WarpX::maxwell_fdtd_solver_id);
}
@@ -344,15 +341,15 @@ WarpX::EvolveB (int lev, Real dt, DtType typ)
}
void
-WarpX::EvolveE (Real dt, DtType typ)
+WarpX::EvolveE (Real dt)
{
for (int lev = 0; lev <= finest_level; ++lev) {
- EvolveE(lev, dt, typ);
+ EvolveE(lev, dt);
}
}
void
-WarpX::EvolveE (int lev, Real dt, DtType typ)
+WarpX::EvolveE (int lev, Real dt)
{
BL_PROFILE("WarpX::EvolveE()");
@@ -360,7 +357,7 @@ WarpX::EvolveE (int lev, Real dt, DtType typ)
const int norder = 2;
static constexpr Real c2 = PhysConst::c*PhysConst::c;
const Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * dt;
- const Real foo = (PhysConst::c*PhysConst::c) * dt;
+ const Real c2dt = (PhysConst::c*PhysConst::c) * dt;
int npatches = (lev == 0) ? 1 : 2;
@@ -368,7 +365,7 @@ WarpX::EvolveE (int lev, Real dt, DtType typ)
{
int patch_level = (ipatch == 0) ? lev : lev-1;
const std::array<Real,3>& dx = WarpX::CellSize(patch_level);
- const std::array<Real,3> dtsdx_c2 {foo/dx[0], foo/dx[1], foo/dx[2]};
+ const std::array<Real,3> dtsdx_c2 {c2dt/dx[0], c2dt/dx[1], c2dt/dx[2]};
MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz, *jx, *jy, *jz, *F;
if (ipatch == 0)
@@ -454,16 +451,14 @@ WarpX::EvolveE (int lev, Real dt, DtType typ)
{
pml[lev]->ExchangeF(F_fp[lev].get(), F_cp[lev].get());
- const int dttype = static_cast<int>(typ);
-
for (int ipatch = 0; ipatch < npatches; ++ipatch)
{
const auto& pml_B = (ipatch==0) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp();
const auto& pml_E = (ipatch==0) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
- const auto& sigba = (ipatch==0) ? pml[lev]->GetMultiSigmaBox_fp()
- : pml[lev]->GetMultiSigmaBox_cp();
- const MultiFab* pml_F = (ipatch==0) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp();
-
+ const auto& pml_F = (ipatch==0) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp();
+ int patch_level = (ipatch == 0) ? lev : lev-1;
+ const std::array<Real,3>& dx = WarpX::CellSize(patch_level);
+ const std::array<Real,3> dtsdx_c2 {c2dt/dx[0], c2dt/dx[1], c2dt/dx[2]};
#ifdef _OPENMP
#pragma omp parallel
#endif
@@ -483,7 +478,7 @@ WarpX::EvolveE (int lev, Real dt, DtType typ)
BL_TO_FORTRAN_3D((*pml_B[0])[mfi]),
BL_TO_FORTRAN_3D((*pml_B[1])[mfi]),
BL_TO_FORTRAN_3D((*pml_B[2])[mfi]),
- WRPX_PML_SIGMA_TO_FORTRAN(sigba[mfi],dttype));
+ &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]);
if (pml_F)
{
@@ -495,8 +490,7 @@ WarpX::EvolveE (int lev, Real dt, DtType typ)
BL_TO_FORTRAN_3D((*pml_E[1])[mfi]),
BL_TO_FORTRAN_3D((*pml_E[2])[mfi]),
BL_TO_FORTRAN_3D((*pml_F )[mfi]),
- WRPX_PML_SIGMA_STAR_TO_FORTRAN(sigba[mfi],dttype),
- &c2);
+ &dtsdx_c2[0]);
}
}
}
@@ -504,17 +498,17 @@ WarpX::EvolveE (int lev, Real dt, DtType typ)
}
void
-WarpX::EvolveF (Real dt, DtType typ)
+WarpX::EvolveF (Real dt, DtType dt_type)
{
if (!do_dive_cleaning) return;
for (int lev = 0; lev <= finest_level; ++lev) {
- EvolveF(lev, dt, typ);
+ EvolveF(lev, dt, dt_type);
}
}
void
-WarpX::EvolveF (int lev, Real dt, DtType typ)
+WarpX::EvolveF (int lev, Real dt, DtType dt_type)
{
if (!do_dive_cleaning) return;
@@ -529,6 +523,7 @@ WarpX::EvolveF (int lev, Real dt, DtType typ)
{
int patch_level = (ipatch == 0) ? lev : lev-1;
const auto& dx = WarpX::CellSize(patch_level);
+ const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]};
MultiFab *Ex, *Ey, *Ez, *rho, *F;
if (ipatch == 0)
@@ -548,20 +543,17 @@ WarpX::EvolveF (int lev, Real dt, DtType typ)
F = F_cp[lev].get();
}
+ const int rhocomp = (dt_type == DtType::FirstHalf) ? 0 : 1;
+
MultiFab src(rho->boxArray(), rho->DistributionMap(), 1, 0);
ComputeDivE(src, 0, {Ex,Ey,Ez}, dx);
- MultiFab::Saxpy(src, -mu_c2, *rho, 0, 0, 1, 0);
+ MultiFab::Saxpy(src, -mu_c2, *rho, rhocomp, 0, 1, 0);
MultiFab::Saxpy(*F, dt, src, 0, 0, 1, 0);
if (do_pml && pml[lev]->ok())
{
- const int dttype = static_cast<int>(typ);
-
const auto& pml_F = (ipatch==0) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp();
const auto& pml_E = (ipatch==0) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
- const auto& sigba = (ipatch==0) ? pml[lev]->GetMultiSigmaBox_fp()
- : pml[lev]->GetMultiSigmaBox_cp();
-
#ifdef _OPENMP
#pragma omp parallel
@@ -574,8 +566,73 @@ WarpX::EvolveF (int lev, Real dt, DtType typ)
BL_TO_FORTRAN_ANYD((*pml_E[0])[mfi]),
BL_TO_FORTRAN_ANYD((*pml_E[1])[mfi]),
BL_TO_FORTRAN_ANYD((*pml_E[2])[mfi]),
- WRPX_PML_SIGMA_TO_FORTRAN(sigba[mfi],dttype),
- &c2inv);
+ &dtsdx[0], &dtsdx[1], &dtsdx[2]);
+ }
+ }
+ }
+}
+
+void
+WarpX::DampPML ()
+{
+ if (!do_pml) return;
+
+ for (int lev = 0; lev <= finest_level; ++lev) {
+ DampPML(lev);
+ }
+}
+
+void
+WarpX::DampPML (int lev)
+{
+ if (!do_pml) return;
+
+ BL_PROFILE("WarpX::DampPML()");
+
+ int npatches = (lev == 0) ? 1 : 2;
+
+ for (int ipatch = 0; ipatch < npatches; ++ipatch)
+ {
+ if (pml[lev]->ok())
+ {
+ const auto& pml_E = (ipatch==0) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
+ const auto& pml_B = (ipatch==0) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp();
+ const auto& pml_F = (ipatch==0) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp();
+ const auto& sigba = (ipatch==0) ? pml[lev]->GetMultiSigmaBox_fp()
+ : pml[lev]->GetMultiSigmaBox_cp();
+
+#ifdef _OPENMP
+#pragma omp parallel
+#endif
+ for ( MFIter mfi(*pml_E[0],true); mfi.isValid(); ++mfi )
+ {
+ const Box& tex = mfi.tilebox(Ex_nodal_flag);
+ const Box& tey = mfi.tilebox(Ey_nodal_flag);
+ const Box& tez = mfi.tilebox(Ez_nodal_flag);
+ const Box& tbx = mfi.tilebox(Bx_nodal_flag);
+ const Box& tby = mfi.tilebox(By_nodal_flag);
+ const Box& tbz = mfi.tilebox(Bz_nodal_flag);
+
+ WRPX_DAMP_PML(tex.loVect(), tex.hiVect(),
+ tey.loVect(), tey.hiVect(),
+ tez.loVect(), tez.hiVect(),
+ tbx.loVect(), tbx.hiVect(),
+ tby.loVect(), tby.hiVect(),
+ tbz.loVect(), tbz.hiVect(),
+ BL_TO_FORTRAN_3D((*pml_E[0])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_E[1])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_E[2])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_B[0])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_B[1])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_B[2])[mfi]),
+ WRPX_PML_TO_FORTRAN(sigba[mfi]));
+
+ if (pml_F) {
+ const Box& tnd = mfi.nodaltilebox();
+ WRPX_DAMP_PML_F(tnd.loVect(), tnd.hiVect(),
+ BL_TO_FORTRAN_3D((*pml_F)[mfi]),
+ WRPX_PML_TO_FORTRAN(sigba[mfi]));
+ }
}
}
}
@@ -594,17 +651,11 @@ WarpX::PushParticlesandDepose (Real cur_time)
void
WarpX::PushParticlesandDepose (int lev, Real cur_time)
{
-#ifdef WARPX_USE_PSATD
- MultiFab* prho2 = rho2_fp[lev].get();
-#else
- MultiFab* prho2 = nullptr;
-#endif
-
mypc->Evolve(lev,
*Efield_aux[lev][0],*Efield_aux[lev][1],*Efield_aux[lev][2],
*Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2],
*current_fp[lev][0],*current_fp[lev][1],*current_fp[lev][2],
- rho_fp[lev].get(), prho2, cur_time, dt[lev]);
+ rho_fp[lev].get(), cur_time, dt[lev]);
}
void
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-18 14:11:33 +0000