From 23446df27f0ca209da5003bc75e4a5f9b40f6332 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Mon, 6 Aug 2018 17:09:15 -0700 Subject: WIP: start PML-2SC --- Source/WarpXEvolve.cpp | 35 ++++++++++++++++++++++------------- 1 file changed, 22 insertions(+), 13 deletions(-) (limited to 'Source/WarpXEvolve.cpp') diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index a92f033b5..209e48018 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -113,12 +113,12 @@ 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(dt[0]); + 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} + EvolveB(0.5*dt[0]); // We now have B^{n+1} FillBoundaryB(); #endif @@ -222,15 +222,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()"); @@ -303,6 +303,8 @@ WarpX::EvolveB (int lev, Real dt, DtType typ) } } +// xxxxx +#if 0 if (do_pml && pml[lev]->ok()) { const int dttype = static_cast(typ); @@ -341,18 +343,19 @@ WarpX::EvolveB (int lev, Real dt, DtType typ) } } } +#endif } 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()"); @@ -450,6 +453,7 @@ WarpX::EvolveE (int lev, Real dt, DtType typ) } } +#if 0 if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeF(F_fp[lev].get(), F_cp[lev].get()); @@ -501,20 +505,21 @@ WarpX::EvolveE (int lev, Real dt, DtType typ) } } } +#endif } void -WarpX::EvolveF (Real dt, DtType typ) +WarpX::EvolveF (Real dt) { if (!do_dive_cleaning) return; for (int lev = 0; lev <= finest_level; ++lev) { - EvolveF(lev, dt, typ); + EvolveF(lev, dt); } } void -WarpX::EvolveF (int lev, Real dt, DtType typ) +WarpX::EvolveF (int lev, Real dt) { if (!do_dive_cleaning) return; @@ -553,6 +558,8 @@ WarpX::EvolveF (int lev, Real dt, DtType typ) MultiFab::Saxpy(src, -mu_c2, *rho, 0, 0, 1, 0); MultiFab::Saxpy(*F, dt, src, 0, 0, 1, 0); +// xxxxx +#if 0 if (do_pml && pml[lev]->ok()) { const int dttype = static_cast(typ); @@ -578,6 +585,8 @@ WarpX::EvolveF (int lev, Real dt, DtType typ) &c2inv); } } +#endif + } } -- cgit v1.2.3 From f91570b82488a200df6b320942d19fe828feb6a7 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Tue, 7 Aug 2018 15:52:35 -0700 Subject: WIP: first pass of PMl_2SC without div e cleaning --- Source/WarpX.H | 3 + Source/WarpXEvolve.cpp | 92 ++++++-- Source/WarpXPML.H | 38 ++-- Source/WarpX_f.H | 44 ++-- Source/WarpX_pml.F90 | 575 +++++++++++++++++++++++++++++-------------------- 5 files changed, 463 insertions(+), 289 deletions(-) (limited to 'Source/WarpXEvolve.cpp') diff --git a/Source/WarpX.H b/Source/WarpX.H index 40cf38a64..f604e2bcd 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -115,6 +115,9 @@ public: void EvolveF ( amrex::Real dt); void EvolveF (int lev, amrex::Real dt); + void DampPML (); + void DampPML (int lev); + void PushParticlesandDepose (int lev, amrex::Real cur_time); void PushParticlesandDepose ( amrex::Real cur_time); diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 209e48018..b6d26d3a5 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -119,6 +119,10 @@ WarpX::EvolveEM (int numsteps) EvolveE(dt[0]); // We now have E^{n+1} FillBoundaryE(); EvolveB(0.5*dt[0]); // We now have B^{n+1} + if (do_pml) { + DampPML(); + FillBoundaryE(); + } FillBoundaryB(); #endif @@ -303,18 +307,12 @@ WarpX::EvolveB (int lev, Real dt) } } -// xxxxx -#if 0 if (do_pml && pml[lev]->ok()) { - const int dttype = static_cast(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& dx = WarpX::CellSize(patch_level); const std::array dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]}; @@ -337,13 +335,11 @@ WarpX::EvolveB (int lev, Real dt) 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); } } } -#endif } void @@ -363,7 +359,7 @@ WarpX::EvolveE (int lev, Real dt) 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; @@ -371,7 +367,7 @@ WarpX::EvolveE (int lev, Real dt) { int patch_level = (ipatch == 0) ? lev : lev-1; const std::array& dx = WarpX::CellSize(patch_level); - const std::array dtsdx_c2 {foo/dx[0], foo/dx[1], foo/dx[2]}; + const std::array 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) @@ -453,21 +449,18 @@ WarpX::EvolveE (int lev, Real dt) } } -#if 0 if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeF(F_fp[lev].get(), F_cp[lev].get()); - const int dttype = static_cast(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& dx = WarpX::CellSize(patch_level); + const std::array dtsdx_c2 {c2dt/dx[0], c2dt/dx[1], c2dt/dx[2]}; #ifdef _OPENMP #pragma omp parallel #endif @@ -487,10 +480,12 @@ WarpX::EvolveE (int lev, Real dt) 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) { + amrex::Abort("pml_F"); +#if 0 WRPX_PUSH_PML_EVEC_F( tex.loVect(), tex.hiVect(), tey.loVect(), tey.hiVect(), @@ -501,11 +496,11 @@ WarpX::EvolveE (int lev, Real dt) BL_TO_FORTRAN_3D((*pml_F )[mfi]), WRPX_PML_SIGMA_STAR_TO_FORTRAN(sigba[mfi],dttype), &c2); +#endif } } } } -#endif } void @@ -590,6 +585,65 @@ WarpX::EvolveF (int lev, Real dt) } } +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(); + // xxxxx 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])); + } + } + } +} + void WarpX::PushParticlesandDepose (Real cur_time) { diff --git a/Source/WarpXPML.H b/Source/WarpXPML.H index 900d94822..44088e736 100644 --- a/Source/WarpXPML.H +++ b/Source/WarpXPML.H @@ -8,35 +8,21 @@ #if (AMREX_SPACEDIM == 3) -#define WRPX_PML_SIGMA_STAR_TO_FORTRAN(x,y) \ - (x).sigma_star_fac1[y][0].data(), (x).sigma_star_fac1[y][0].m_lo, (x).sigma_star_fac1[y][0].m_hi, \ - (x).sigma_star_fac2[y][0].data(), (x).sigma_star_fac2[y][0].m_lo, (x).sigma_star_fac2[y][0].m_hi, \ - (x).sigma_star_fac1[y][1].data(), (x).sigma_star_fac1[y][1].m_lo, (x).sigma_star_fac1[y][1].m_hi, \ - (x).sigma_star_fac2[y][1].data(), (x).sigma_star_fac2[y][1].m_lo, (x).sigma_star_fac2[y][1].m_hi, \ - (x).sigma_star_fac1[y][2].data(), (x).sigma_star_fac1[y][2].m_lo, (x).sigma_star_fac1[y][2].m_hi, \ - (x).sigma_star_fac2[y][2].data(), (x).sigma_star_fac2[y][2].m_lo, (x).sigma_star_fac2[y][2].m_hi - -#define WRPX_PML_SIGMA_TO_FORTRAN(x,y) \ - (x).sigma_fac1[y][0].data(), (x).sigma_fac1[y][0].m_lo, (x).sigma_fac1[y][0].m_hi, \ - (x).sigma_fac2[y][0].data(), (x).sigma_fac2[y][0].m_lo, (x).sigma_fac2[y][0].m_hi, \ - (x).sigma_fac1[y][1].data(), (x).sigma_fac1[y][1].m_lo, (x).sigma_fac1[y][1].m_hi, \ - (x).sigma_fac2[y][1].data(), (x).sigma_fac2[y][1].m_lo, (x).sigma_fac2[y][1].m_hi, \ - (x).sigma_fac1[y][2].data(), (x).sigma_fac1[y][2].m_lo, (x).sigma_fac1[y][2].m_hi, \ - (x).sigma_fac2[y][2].data(), (x).sigma_fac2[y][2].m_lo, (x).sigma_fac2[y][2].m_hi +#define WRPX_PML_TO_FORTRAN(x) \ + (x).sigma_fac[0].data(), (x).sigma_fac[0].m_lo, (x).sigma_fac[0].m_hi, \ + (x).sigma_fac[1].data(), (x).sigma_fac[1].m_lo, (x).sigma_fac[1].m_hi, \ + (x).sigma_fac[2].data(), (x).sigma_fac[2].m_lo, (x).sigma_fac[2].m_hi, \ + (x).sigma_star_fac[0].data(), (x).sigma_star_fac[0].m_lo, (x).sigma_star_fac[0].m_hi, \ + (x).sigma_star_fac[1].data(), (x).sigma_star_fac[1].m_lo, (x).sigma_star_fac[1].m_hi, \ + (x).sigma_star_fac[2].data(), (x).sigma_star_fac[2].m_lo, (x).sigma_star_fac[2].m_hi #else -#define WRPX_PML_SIGMA_STAR_TO_FORTRAN(x,y) \ - (x).sigma_star_fac1[y][0].data(), (x).sigma_star_fac1[y][0].m_lo, (x).sigma_star_fac1[y][0].m_hi, \ - (x).sigma_star_fac2[y][0].data(), (x).sigma_star_fac2[y][0].m_lo, (x).sigma_star_fac2[y][0].m_hi, \ - (x).sigma_star_fac1[y][1].data(), (x).sigma_star_fac1[y][1].m_lo, (x).sigma_star_fac1[y][1].m_hi, \ - (x).sigma_star_fac2[y][1].data(), (x).sigma_star_fac2[y][1].m_lo, (x).sigma_star_fac2[y][1].m_hi - -#define WRPX_PML_SIGMA_TO_FORTRAN(x,y) \ - (x).sigma_fac1[y][0].data(), (x).sigma_fac1[y][0].m_lo, (x).sigma_fac1[y][0].m_hi, \ - (x).sigma_fac2[y][0].data(), (x).sigma_fac2[y][0].m_lo, (x).sigma_fac2[y][0].m_hi, \ - (x).sigma_fac1[y][1].data(), (x).sigma_fac1[y][1].m_lo, (x).sigma_fac1[y][1].m_hi, \ - (x).sigma_fac2[y][1].data(), (x).sigma_fac2[y][1].m_lo, (x).sigma_fac2[y][1].m_hi +#define WRPX_PML_TO_FORTRAN(x) \ + (x).sigma_fac[0].data(), (x).sigma_fac[0].m_lo, (x).sigma_fac[0].m_hi, \ + (x).sigma_fac[1].data(), (x).sigma_fac[1].m_lo, (x).sigma_fac[1].m_hi, \ + (x).sigma_star_fac[0].data(), (x).sigma_star_fac[0].m_lo, (x).sigma_star_fac[0].m_hi, \ + (x).sigma_star_fac[1].data(), (x).sigma_star_fac[1].m_lo, (x).sigma_star_fac[1].m_hi #endif diff --git a/Source/WarpX_f.H b/Source/WarpX_f.H index 377d8fd11..bdebffc68 100644 --- a/Source/WarpX_f.H +++ b/Source/WarpX_f.H @@ -26,6 +26,7 @@ #define WRPX_PUSH_PML_EVEC warpx_push_pml_evec_3d #define WRPX_PUSH_PML_EVEC_F warpx_push_pml_evec_f_3d #define WRPX_PUSH_PML_F warpx_push_pml_f_3d +#define WRPX_DAMP_PML warpx_damp_pml_3d #define WRPX_SUM_FINE_TO_CRSE_NODAL warpx_sum_fine_to_crse_nodal_3d #define WRPX_ZERO_OUT_BNDRY warpx_zero_out_bndry_3d @@ -56,6 +57,7 @@ #define WRPX_PUSH_PML_EVEC warpx_push_pml_evec_2d #define WRPX_PUSH_PML_EVEC_F warpx_push_pml_evec_f_2d #define WRPX_PUSH_PML_F warpx_push_pml_f_2d +#define WRPX_DAMP_PML warpx_damp_pml_2d #define WRPX_SUM_FINE_TO_CRSE_NODAL warpx_sum_fine_to_crse_nodal_2d #define WRPX_ZERO_OUT_BNDRY warpx_zero_out_bndry_2d @@ -333,14 +335,6 @@ extern "C" BL_FORT_FAB_ARG_3D(bx), BL_FORT_FAB_ARG_3D(by), BL_FORT_FAB_ARG_3D(bz), - const amrex::Real* sigx1, int sigx1_lo, int sigx1_hi, - const amrex::Real* sigx2, int sigx2_lo, int sigx2_hi, -#if (AMREX_SPACEDIM == 3) - const amrex::Real* sigy1, int sigy1_lo, int sigy1_hi, - const amrex::Real* sigy2, int sigy2_lo, int sigy2_hi, -#endif - const amrex::Real* sigz1, int sigz1_lo, int sigz1_hi, - const amrex::Real* sigz2, int sigz2_lo, int sigz2_hi, const amrex::Real* dtsdx, const amrex::Real* dtsdy, const amrex::Real* dtsdz, @@ -356,14 +350,9 @@ extern "C" const BL_FORT_FAB_ARG_3D(bx), const BL_FORT_FAB_ARG_3D(by), const BL_FORT_FAB_ARG_3D(bz), - const amrex::Real* sigx1, int sigx1_lo, int sigx1_hi, - const amrex::Real* sigx2, int sigx2_lo, int sigx2_hi, -#if (AMREX_SPACEDIM == 3) - const amrex::Real* sigy1, int sigy1_lo, int sigy1_hi, - const amrex::Real* sigy2, int sigy2_lo, int sigy2_hi, -#endif - const amrex::Real* sigz1, int sigz1_lo, int sigz1_hi, - const amrex::Real* sigz2, int sigz2_lo, int sigz2_hi); + const amrex::Real* dtsdx, + const amrex::Real* dtsdy, + const amrex::Real* dtsdz); void WRPX_PUSH_PML_EVEC_F(const int* xlo, const int* xhi, const int* ylo, const int* yhi, @@ -397,6 +386,29 @@ extern "C" const amrex::Real* sigz2, int sigz2_lo, int sigz2_hi, const amrex::Real* c2inv); + void WRPX_DAMP_PML (const int* texlo, const int* texhi, + const int* teylo, const int* teyhi, + const int* tezlo, const int* tezhi, + const int* tbxlo, const int* tbxhi, + const int* tbylo, const int* tbyhi, + const int* tbzlo, const int* tbzhi, + amrex::Real* ex, const int* exlo, const int* exhi, + amrex::Real* ey, const int* eylo, const int* eyhi, + amrex::Real* ez, const int* ezlo, const int* ezhi, + amrex::Real* bx, const int* bxlo, const int* bxhi, + amrex::Real* by, const int* bylo, const int* byhi, + amrex::Real* bz, const int* bzlo, const int* bzhi, + const amrex::Real* sigex, int sigex_lo, int sigex_hi, +#if (AMREX_SPACEDIM == 3) + const amrex::Real* sigey, int sigey_lo, int sigey_hi, +#endif + const amrex::Real* sigez, int sigez_lo, int sigez_hi, + const amrex::Real* sigbx, int sigbx_lo, int sigbx_hi, +#if (AMREX_SPACEDIM == 3) + const amrex::Real* sigby, int sigby_lo, int sigby_hi, +#endif + const amrex::Real* sigbz, int sigbz_lo, int sigbz_hi); + void WRPX_SYNC_CURRENT (const int* lo, const int* hi, BL_FORT_FAB_ARG_ANYD(crse), const BL_FORT_FAB_ARG_ANYD(fine), diff --git a/Source/WarpX_pml.F90 b/Source/WarpX_pml.F90 index c87ba1442..8d929aa6c 100644 --- a/Source/WarpX_pml.F90 +++ b/Source/WarpX_pml.F90 @@ -14,32 +14,17 @@ contains & Bx, Bxlo, Bxhi, & & By, Bylo, Byhi, & & Bz, Bzlo, Bzhi, & - & sigx1, sigx1_lo, sigx1_hi, & - & sigx2, sigx2_lo, sigx2_hi, & - & sigy1, sigy1_lo, sigy1_hi, & - & sigy2, sigy2_lo, sigy2_hi, & - & sigz1, sigz1_lo, sigz1_hi, & - & sigz2, sigz2_lo, sigz2_hi, & & dtsdx, dtsdy, dtsdz, solver_type) & bind(c,name='warpx_push_pml_bvec_3d') integer, intent(in) :: xlo(3), xhi(3), ylo(3), yhi(3), zlo(3), zhi(3), & Exlo(3), Exhi(3), Eylo(3), Eyhi(3), Ezlo(3), Ezhi(3), & Bxlo(3), Bxhi(3), Bylo(3), Byhi(3), Bzlo(3), Bzhi(3), solver_type - integer, intent(in), value :: sigx1_lo, sigx1_hi, sigx2_lo, sigx2_hi - integer, intent(in), value :: sigy1_lo, sigy1_hi, sigy2_lo, sigy2_hi - integer, intent(in), value :: sigz1_lo, sigz1_hi, sigz2_lo, sigz2_hi real(amrex_real), intent(in ) :: Ex (Exlo(1):Exhi(1),Exlo(2):Exhi(2),Exlo(3):Exhi(3),3) real(amrex_real), intent(in ) :: Ey (Eylo(1):Eyhi(1),Eylo(2):Eyhi(2),Eylo(3):Eyhi(3),3) real(amrex_real), intent(in ) :: Ez (Ezlo(1):Ezhi(1),Ezlo(2):Ezhi(2),Ezlo(3):Ezhi(3),3) real(amrex_real), intent(inout) :: Bx (Bxlo(1):Bxhi(1),Bxlo(2):Bxhi(2),Bxlo(3):Bxhi(3),2) real(amrex_real), intent(inout) :: By (Bylo(1):Byhi(1),Bylo(2):Byhi(2),Bylo(3):Byhi(3),2) real(amrex_real), intent(inout) :: Bz (Bzlo(1):Bzhi(1),Bzlo(2):Bzhi(2),Bzlo(3):Bzhi(3),2) - real(amrex_real), intent(in) :: sigx1(sigx1_lo:sigx1_hi) - real(amrex_real), intent(in) :: sigx2(sigx2_lo:sigx2_hi) - real(amrex_real), intent(in) :: sigy1(sigy1_lo:sigy1_hi) - real(amrex_real), intent(in) :: sigy2(sigy2_lo:sigy2_hi) - real(amrex_real), intent(in) :: sigz1(sigz1_lo:sigz1_hi) - real(amrex_real), intent(in) :: sigz2(sigz2_lo:sigz2_hi) real(amrex_real), intent(in) :: dtsdx, dtsdy, dtsdz integer :: i, j, k @@ -56,10 +41,10 @@ contains do k = xlo(3), xhi(3) do j = xlo(2), xhi(2) do i = xlo(1), xhi(1) - Bx(i,j,k,1) = sigy1(j)*Bx(i,j,k,1) - sigy2(j)*(Ez(i,j+1,k ,1)+Ez(i,j+1,k ,2)+Ez(i,j+1,k ,3) & - & -Ez(i,j ,k ,1)-Ez(i,j ,k ,2)-Ez(i,j ,k ,3)) - Bx(i,j,k,2) = sigz1(k)*Bx(i,j,k,2) + sigz2(k)*(Ey(i,j ,k+1,1)+Ey(i,j ,k+1,2)+Ey(i,j ,k+1,3) & - & -Ey(i,j ,k ,1)-Ey(i,j ,k ,2)-Ey(i,j ,k ,3)) + Bx(i,j,k,1) = Bx(i,j,k,1) - dtsdy*(Ez(i,j+1,k ,1)+Ez(i,j+1,k ,2)+Ez(i,j+1,k ,3) & + & -Ez(i,j ,k ,1)-Ez(i,j ,k ,2)-Ez(i,j ,k ,3)) + Bx(i,j,k,2) = Bx(i,j,k,2) + dtsdz*(Ey(i,j ,k+1,1)+Ey(i,j ,k+1,2)+Ey(i,j ,k+1,3) & + & -Ey(i,j ,k ,1)-Ey(i,j ,k ,2)-Ey(i,j ,k ,3)) end do end do end do @@ -67,10 +52,10 @@ contains do k = ylo(3), yhi(3) do j = ylo(2), yhi(2) do i = ylo(1), yhi(1) - By(i,j,k,1) = sigz1(k)*By(i,j,k,1) - sigz2(k)*(Ex(i ,j,k+1,1)+Ex(i ,j,k+1,2)+Ex(i ,j,k+1,3) & - & -Ex(i ,j,k ,1)-Ex(i ,j,k ,2)-Ex(i ,j,k ,3)) - By(i,j,k,2) = sigx1(i)*By(i,j,k,2) + sigx2(i)*(Ez(i+1,j,k ,1)+Ez(i+1,j,k ,2)+Ez(i+1,j,k ,3) & - & -Ez(i ,j,k ,1)-Ez(i ,j,k ,2)-Ez(i ,j,k ,3)) + By(i,j,k,1) = By(i,j,k,1) - dtsdz*(Ex(i ,j,k+1,1)+Ex(i ,j,k+1,2)+Ex(i ,j,k+1,3) & + & -Ex(i ,j,k ,1)-Ex(i ,j,k ,2)-Ex(i ,j,k ,3)) + By(i,j,k,2) = By(i,j,k,2) + dtsdx*(Ez(i+1,j,k ,1)+Ez(i+1,j,k ,2)+Ez(i+1,j,k ,3) & + & -Ez(i ,j,k ,1)-Ez(i ,j,k ,2)-Ez(i ,j,k ,3)) end do end do end do @@ -78,10 +63,10 @@ contains do k = zlo(3), zhi(3) do j = zlo(2), zhi(2) do i = zlo(1), zhi(1) - Bz(i,j,k,1) = sigx1(i)*Bz(i,j,k,1) - sigx2(i)*(Ey(i+1,j ,k,1)+Ey(i+1,j ,k,2)+Ey(i+1,j ,k,3) & - & -Ey(i ,j ,k,1)-Ey(i ,j ,k,2)-Ey(i ,j ,k,3)) - Bz(i,j,k,2) = sigy1(j)*Bz(i,j,k,2) + sigy2(j)*(Ex(i ,j+1,k,1)+Ex(i ,j+1,k,2)+Ex(i ,j+1,k,3) & - & -Ex(i ,j ,k,1)-Ex(i ,j ,k,2)-Ex(i ,j ,k,3)) + Bz(i,j,k,1) = Bz(i,j,k,1) - dtsdx*(Ey(i+1,j ,k,1)+Ey(i+1,j ,k,2)+Ey(i+1,j ,k,3) & + & -Ey(i ,j ,k,1)-Ey(i ,j ,k,2)-Ey(i ,j ,k,3)) + Bz(i,j,k,2) = Bz(i,j,k,2) + dtsdy*(Ex(i ,j+1,k,1)+Ex(i ,j+1,k,2)+Ex(i ,j+1,k,3) & + & -Ex(i ,j ,k,1)-Ex(i ,j ,k,2)-Ex(i ,j ,k,3)) end do end do end do @@ -109,47 +94,60 @@ contains alphay = 1. - 2.*betayx - 2.* betayz - 4.*gammay alphaz = 1. - 2.*betazx - 2.* betazy - 4.*gammaz + betaxy = dtsdx*betaxy + betaxz = dtsdx*betaxz + betayx = dtsdy*betayx + betayz = dtsdy*betayz + betazx = dtsdz*betazx + betazy = dtsdz*betazy + alphax = dtsdx*alphax + alphay = dtsdy*alphay + alphaz = dtsdz*alphaz + gammax = dtsdx*gammax + gammay = dtsdy*gammay + gammaz = dtsdz*gammaz + do k = xlo(3), xhi(3) do j = xlo(2), xhi(2) do i = xlo(1), xhi(1) - Bx(i,j,k,1) = sigy1(j)*Bx(i,j,k,1) - sigy2(j)*(alphay*(Ez(i ,j+1,k ,1)+Ez(i ,j+1,k ,2)+Ez(i ,j+1,k ,3) & - -Ez(i ,j ,k ,1)-Ez(i ,j ,k ,2)-Ez(i ,j ,k ,3)) & - +betayx*(Ez(i+1,j+1,k ,1)+Ez(i+1,j+1,k ,2)+Ez(i+1,j+1,k ,3) & - -Ez(i+1,j ,k ,1)-Ez(i+1,j ,k ,2)-Ez(i+1,j ,k ,3) & - +Ez(i-1,j+1,k ,1)+Ez(i-1,j+1,k ,2)+Ez(i-1,j+1,k ,3) & - -Ez(i-1,j ,k ,1)-Ez(i-1,j ,k ,2)-Ez(i-1,j ,k ,3)) & - +betayz*(Ez(i ,j+1,k+1,1)+Ez(i ,j+1,k+1,2)+Ez(i ,j+1,k+1,3) & - -Ez(i ,j ,k+1,1)-Ez(i ,j ,k+1,2)-Ez(i ,j ,k+1,3) & - +Ez(i ,j+1,k-1,1)+Ez(i ,j+1,k-1,2)+Ez(i ,j+1,k-1,3) & - -Ez(i ,j ,k-1,1)-Ez(i ,j ,k-1,2)-Ez(i ,j ,k-1,3)) & - +gammay*(Ez(i+1,j+1,k+1,1)+Ez(i+1,j+1,k+1,2)+Ez(i+1,j+1,k+1,3) & - -Ez(i+1,j ,k+1,1)-Ez(i+1,j ,k+1,2)-Ez(i+1,j ,k+1,3) & - +Ez(i-1,j+1,k+1,1)+Ez(i-1,j+1,k+1,2)+Ez(i-1,j+1,k+1,3) & - -Ez(i-1,j ,k+1,1)-Ez(i-1,j ,k+1,2)-Ez(i-1,j ,k+1,3) & - +Ez(i+1,j+1,k-1,1)+Ez(i+1,j+1,k-1,2)+Ez(i+1,j+1,k-1,3) & - -Ez(i+1,j ,k-1,1)-Ez(i+1,j ,k-1,2)-Ez(i+1,j ,k-1,3) & - +Ez(i-1,j+1,k-1,1)+Ez(i-1,j+1,k-1,2)+Ez(i-1,j+1,k-1,3) & - -Ez(i-1,j ,k-1,1)-Ez(i-1,j ,k-1,2)-Ez(i-1,j ,k-1,3))) - - - Bx(i,j,k,2) = sigz1(k)*Bx(i,j,k,2) + sigz2(k)*(alphaz*(Ey(i ,j ,k+1,1)+Ey(i ,j ,k+1,2)+Ey(i ,j ,k+1,3) & - -Ey(i ,j ,k ,1)-Ey(i ,j ,k ,2)-Ey(i ,j ,k ,3)) & - +betazx*(Ey(i+1,j ,k+1,1)+Ey(i+1,j ,k+1,2)+Ey(i+1,j ,k+1,3) & - -Ey(i+1,j ,k ,1)-Ey(i+1,j ,k ,2)-Ey(i+1,j ,k ,3) & - +Ey(i-1,j ,k+1,1)+Ey(i-1,j ,k+1,2)+Ey(i-1,j ,k+1,3) & - -Ey(i-1,j ,k ,1)-Ey(i-1,j ,k ,2)-Ey(i-1,j ,k ,3)) & - +betazy*(Ey(i ,j+1,k+1,1)+Ey(i ,j+1,k+1,2)+Ey(i ,j+1,k+1,3) & - -Ey(i ,j+1,k ,1)-Ey(i ,j+1,k ,2)-Ey(i ,j+1,k ,3) & - +Ey(i ,j-1,k+1,1)+Ey(i ,j-1,k+1,2)+Ey(i ,j-1,k+1,3) & - -Ey(i ,j-1,k ,1)-Ey(i ,j-1,k ,2)-Ey(i ,j-1,k ,3)) & - +gammaz*(Ey(i+1,j+1,k+1,1)+Ey(i+1,j+1,k+1,2)+Ey(i+1,j+1,k+1,3) & - -Ey(i+1,j+1,k ,1)-Ey(i+1,j+1,k ,2)-Ey(i+1,j+1,k ,3) & - +Ey(i-1,j+1,k+1,1)+Ey(i-1,j+1,k+1,2)+Ey(i-1,j+1,k+1,3) & - -Ey(i-1,j+1,k ,1)-Ey(i-1,j+1,k ,2)-Ey(i-1,j+1,k ,3) & - +Ey(i+1,j-1,k+1,1)+Ey(i+1,j-1,k+1,2)+Ey(i+1,j-1,k+1,3) & - -Ey(i+1,j-1,k ,1)-Ey(i+1,j-1,k ,2)-Ey(i+1,j-1,k ,3) & - +Ey(i-1,j-1,k+1,1)+Ey(i-1,j-1,k+1,2)+Ey(i-1,j-1,k+1,3) & - -Ey(i-1,j-1,k ,1)-Ey(i-1,j-1,k ,2)-Ey(i-1,j-1,k ,3))) + Bx(i,j,k,1) = Bx(i,j,k,1) - (alphay*(Ez(i ,j+1,k ,1)+Ez(i ,j+1,k ,2)+Ez(i ,j+1,k ,3) & + -Ez(i ,j ,k ,1)-Ez(i ,j ,k ,2)-Ez(i ,j ,k ,3)) & + +betayx*(Ez(i+1,j+1,k ,1)+Ez(i+1,j+1,k ,2)+Ez(i+1,j+1,k ,3) & + -Ez(i+1,j ,k ,1)-Ez(i+1,j ,k ,2)-Ez(i+1,j ,k ,3) & + +Ez(i-1,j+1,k ,1)+Ez(i-1,j+1,k ,2)+Ez(i-1,j+1,k ,3) & + -Ez(i-1,j ,k ,1)-Ez(i-1,j ,k ,2)-Ez(i-1,j ,k ,3)) & + +betayz*(Ez(i ,j+1,k+1,1)+Ez(i ,j+1,k+1,2)+Ez(i ,j+1,k+1,3) & + -Ez(i ,j ,k+1,1)-Ez(i ,j ,k+1,2)-Ez(i ,j ,k+1,3) & + +Ez(i ,j+1,k-1,1)+Ez(i ,j+1,k-1,2)+Ez(i ,j+1,k-1,3) & + -Ez(i ,j ,k-1,1)-Ez(i ,j ,k-1,2)-Ez(i ,j ,k-1,3)) & + +gammay*(Ez(i+1,j+1,k+1,1)+Ez(i+1,j+1,k+1,2)+Ez(i+1,j+1,k+1,3) & + -Ez(i+1,j ,k+1,1)-Ez(i+1,j ,k+1,2)-Ez(i+1,j ,k+1,3) & + +Ez(i-1,j+1,k+1,1)+Ez(i-1,j+1,k+1,2)+Ez(i-1,j+1,k+1,3) & + -Ez(i-1,j ,k+1,1)-Ez(i-1,j ,k+1,2)-Ez(i-1,j ,k+1,3) & + +Ez(i+1,j+1,k-1,1)+Ez(i+1,j+1,k-1,2)+Ez(i+1,j+1,k-1,3) & + -Ez(i+1,j ,k-1,1)-Ez(i+1,j ,k-1,2)-Ez(i+1,j ,k-1,3) & + +Ez(i-1,j+1,k-1,1)+Ez(i-1,j+1,k-1,2)+Ez(i-1,j+1,k-1,3) & + -Ez(i-1,j ,k-1,1)-Ez(i-1,j ,k-1,2)-Ez(i-1,j ,k-1,3))) + + + Bx(i,j,k,2) = Bx(i,j,k,2) + (alphaz*(Ey(i ,j ,k+1,1)+Ey(i ,j ,k+1,2)+Ey(i ,j ,k+1,3) & + -Ey(i ,j ,k ,1)-Ey(i ,j ,k ,2)-Ey(i ,j ,k ,3)) & + +betazx*(Ey(i+1,j ,k+1,1)+Ey(i+1,j ,k+1,2)+Ey(i+1,j ,k+1,3) & + -Ey(i+1,j ,k ,1)-Ey(i+1,j ,k ,2)-Ey(i+1,j ,k ,3) & + +Ey(i-1,j ,k+1,1)+Ey(i-1,j ,k+1,2)+Ey(i-1,j ,k+1,3) & + -Ey(i-1,j ,k ,1)-Ey(i-1,j ,k ,2)-Ey(i-1,j ,k ,3)) & + +betazy*(Ey(i ,j+1,k+1,1)+Ey(i ,j+1,k+1,2)+Ey(i ,j+1,k+1,3) & + -Ey(i ,j+1,k ,1)-Ey(i ,j+1,k ,2)-Ey(i ,j+1,k ,3) & + +Ey(i ,j-1,k+1,1)+Ey(i ,j-1,k+1,2)+Ey(i ,j-1,k+1,3) & + -Ey(i ,j-1,k ,1)-Ey(i ,j-1,k ,2)-Ey(i ,j-1,k ,3)) & + +gammaz*(Ey(i+1,j+1,k+1,1)+Ey(i+1,j+1,k+1,2)+Ey(i+1,j+1,k+1,3) & + -Ey(i+1,j+1,k ,1)-Ey(i+1,j+1,k ,2)-Ey(i+1,j+1,k ,3) & + +Ey(i-1,j+1,k+1,1)+Ey(i-1,j+1,k+1,2)+Ey(i-1,j+1,k+1,3) & + -Ey(i-1,j+1,k ,1)-Ey(i-1,j+1,k ,2)-Ey(i-1,j+1,k ,3) & + +Ey(i+1,j-1,k+1,1)+Ey(i+1,j-1,k+1,2)+Ey(i+1,j-1,k+1,3) & + -Ey(i+1,j-1,k ,1)-Ey(i+1,j-1,k ,2)-Ey(i+1,j-1,k ,3) & + +Ey(i-1,j-1,k+1,1)+Ey(i-1,j-1,k+1,2)+Ey(i-1,j-1,k+1,3) & + -Ey(i-1,j-1,k ,1)-Ey(i-1,j-1,k ,2)-Ey(i-1,j-1,k ,3))) end do end do end do @@ -157,44 +155,44 @@ contains do k = ylo(3), yhi(3) do j = ylo(2), yhi(2) do i = ylo(1), yhi(1) - By(i,j,k,1) = sigz1(k)*By(i,j,k,1) - sigz2(k)*(alphaz*(Ex(i ,j ,k+1,1)+Ex(i ,j ,k+1,2)+Ex(i ,j ,k+1,3) & - -Ex(i ,j ,k ,1)-Ex(i ,j ,k ,2)-Ex(i ,j ,k ,3)) & - +betazx*(Ex(i+1,j ,k+1,1)+Ex(i+1,j ,k+1,2)+Ex(i+1,j ,k+1,3) & - -Ex(i+1,j ,k ,1)-Ex(i+1,j ,k ,2)-Ex(i+1,j ,k ,3) & - +Ex(i-1,j ,k+1,1)+Ex(i-1,j ,k+1,2)+Ex(i-1,j ,k+1,3) & - -Ex(i-1,j ,k ,1)-Ex(i-1,j ,k ,2)-Ex(i-1,j ,k ,3)) & - +betaxy*(Ex(i ,j+1,k+1,1)+Ex(i ,j+1,k+1,2)+Ex(i ,j+1,k+1,3) & - -Ex(i ,j+1,k ,1)-Ex(i ,j+1,k ,2)-Ex(i ,j+1,k ,3) & - +Ex(i ,j-1,k+1,1)+Ex(i ,j-1,k+1,2)+Ex(i ,j-1,k+1,3) & - -Ex(i ,j-1,k ,1)-Ex(i ,j-1,k ,2)-Ex(i ,j-1,k ,3)) & - +gammaz*(Ex(i+1,j+1,k+1,1)+Ex(i+1,j+1,k+1,2)+Ex(i+1,j+1,k+1,3) & - -Ex(i+1,j+1,k ,1)-Ex(i+1,j+1,k ,2)-Ex(i+1,j+1,k ,3) & - +Ex(i-1,j+1,k+1,1)+Ex(i-1,j+1,k+1,2)+Ex(i-1,j+1,k+1,3) & - -Ex(i-1,j+1,k ,1)-Ex(i-1,j+1,k ,2)-Ex(i-1,j+1,k ,3) & - +Ex(i+1,j-1,k+1,1)+Ex(i+1,j-1,k+1,2)+Ex(i+1,j-1,k+1,3) & - -Ex(i+1,j-1,k ,1)-Ex(i+1,j-1,k ,2)-Ex(i+1,j-1,k ,3) & - +Ex(i-1,j-1,k+1,1)+Ex(i-1,j-1,k+1,2)+Ex(i-1,j-1,k+1,3) & - -Ex(i-1,j-1,k ,1)-Ex(i-1,j-1,k ,2)-Ex(i-1,j-1,k ,3))) - - - By(i,j,k,2) = sigx1(i)*By(i,j,k,2) + sigx2(i)*(alphax*(Ez(i+1,j ,k ,1)+Ez(i+1,j ,k ,2)+Ez(i+1,j ,k ,3) & - -Ez(i ,j ,k ,1)-Ez(i ,j ,k ,2)-Ez(i ,j ,k ,3)) & - +betaxy*(Ez(i+1,j+1,k ,1)+Ez(i+1,j+1,k ,2)+Ez(i+1,j+1,k ,3) & - -Ez(i ,j+1,k ,1)-Ez(i ,j+1,k ,2)-Ez(i ,j+1,k ,3) & - +Ez(i+1,j-1,k ,1)+Ez(i+1,j-1,k ,2)+Ez(i+1,j-1,k ,3) & - -Ez(i ,j-1,k ,1)-Ez(i ,j-1,k ,2)-Ez(i ,j-1,k ,3)) & - +betaxz*(Ez(i+1,j ,k+1,1)+Ez(i+1,j ,k+1,2)+Ez(i+1,j ,k+1,3) & - -Ez(i ,j ,k+1,1)-Ez(i ,j ,k+1,2)-Ez(i ,j ,k+1,3) & - +Ez(i+1,j ,k-1,1)+Ez(i+1,j ,k-1,2)+Ez(i+1,j ,k-1,3) & - -Ez(i ,j ,k-1,1)-Ez(i ,j ,k-1,2)-Ez(i ,j ,k-1,3)) & - +gammax*(Ez(i+1,j+1,k+1,1)+Ez(i+1,j+1,k+1,2)+Ez(i+1,j+1,k+1,3) & - -Ez(i ,j+1,k+1,1)-Ez(i ,j+1,k+1,2)-Ez(i ,j+1,k+1,3) & - +Ez(i+1,j-1,k+1,1)+Ez(i+1,j-1,k+1,2)+Ez(i+1,j-1,k+1,3) & - -Ez(i ,j-1,k+1,1)-Ez(i ,j-1,k+1,2)-Ez(i ,j-1,k+1,3) & - +Ez(i+1,j+1,k-1,1)+Ez(i+1,j+1,k-1,2)+Ez(i+1,j+1,k-1,3) & - -Ez(i ,j+1,k-1,1)-Ez(i ,j+1,k-1,2)-Ez(i ,j+1,k-1,3) & - +Ez(i+1,j-1,k-1,1)+Ez(i+1,j-1,k-1,2)+Ez(i+1,j-1,k-1,3) & - -Ez(i ,j-1,k-1,1)-Ez(i ,j-1,k-1,2)-Ez(i ,j-1,k-1,3))) + By(i,j,k,1) = By(i,j,k,1) - (alphaz*(Ex(i ,j ,k+1,1)+Ex(i ,j ,k+1,2)+Ex(i ,j ,k+1,3) & + -Ex(i ,j ,k ,1)-Ex(i ,j ,k ,2)-Ex(i ,j ,k ,3)) & + +betazx*(Ex(i+1,j ,k+1,1)+Ex(i+1,j ,k+1,2)+Ex(i+1,j ,k+1,3) & + -Ex(i+1,j ,k ,1)-Ex(i+1,j ,k ,2)-Ex(i+1,j ,k ,3) & + +Ex(i-1,j ,k+1,1)+Ex(i-1,j ,k+1,2)+Ex(i-1,j ,k+1,3) & + -Ex(i-1,j ,k ,1)-Ex(i-1,j ,k ,2)-Ex(i-1,j ,k ,3)) & + +betaxy*(Ex(i ,j+1,k+1,1)+Ex(i ,j+1,k+1,2)+Ex(i ,j+1,k+1,3) & + -Ex(i ,j+1,k ,1)-Ex(i ,j+1,k ,2)-Ex(i ,j+1,k ,3) & + +Ex(i ,j-1,k+1,1)+Ex(i ,j-1,k+1,2)+Ex(i ,j-1,k+1,3) & + -Ex(i ,j-1,k ,1)-Ex(i ,j-1,k ,2)-Ex(i ,j-1,k ,3)) & + +gammaz*(Ex(i+1,j+1,k+1,1)+Ex(i+1,j+1,k+1,2)+Ex(i+1,j+1,k+1,3) & + -Ex(i+1,j+1,k ,1)-Ex(i+1,j+1,k ,2)-Ex(i+1,j+1,k ,3) & + +Ex(i-1,j+1,k+1,1)+Ex(i-1,j+1,k+1,2)+Ex(i-1,j+1,k+1,3) & + -Ex(i-1,j+1,k ,1)-Ex(i-1,j+1,k ,2)-Ex(i-1,j+1,k ,3) & + +Ex(i+1,j-1,k+1,1)+Ex(i+1,j-1,k+1,2)+Ex(i+1,j-1,k+1,3) & + -Ex(i+1,j-1,k ,1)-Ex(i+1,j-1,k ,2)-Ex(i+1,j-1,k ,3) & + +Ex(i-1,j-1,k+1,1)+Ex(i-1,j-1,k+1,2)+Ex(i-1,j-1,k+1,3) & + -Ex(i-1,j-1,k ,1)-Ex(i-1,j-1,k ,2)-Ex(i-1,j-1,k ,3))) + + + By(i,j,k,2) = By(i,j,k,2) + (alphax*(Ez(i+1,j ,k ,1)+Ez(i+1,j ,k ,2)+Ez(i+1,j ,k ,3) & + -Ez(i ,j ,k ,1)-Ez(i ,j ,k ,2)-Ez(i ,j ,k ,3)) & + +betaxy*(Ez(i+1,j+1,k ,1)+Ez(i+1,j+1,k ,2)+Ez(i+1,j+1,k ,3) & + -Ez(i ,j+1,k ,1)-Ez(i ,j+1,k ,2)-Ez(i ,j+1,k ,3) & + +Ez(i+1,j-1,k ,1)+Ez(i+1,j-1,k ,2)+Ez(i+1,j-1,k ,3) & + -Ez(i ,j-1,k ,1)-Ez(i ,j-1,k ,2)-Ez(i ,j-1,k ,3)) & + +betaxz*(Ez(i+1,j ,k+1,1)+Ez(i+1,j ,k+1,2)+Ez(i+1,j ,k+1,3) & + -Ez(i ,j ,k+1,1)-Ez(i ,j ,k+1,2)-Ez(i ,j ,k+1,3) & + +Ez(i+1,j ,k-1,1)+Ez(i+1,j ,k-1,2)+Ez(i+1,j ,k-1,3) & + -Ez(i ,j ,k-1,1)-Ez(i ,j ,k-1,2)-Ez(i ,j ,k-1,3)) & + +gammax*(Ez(i+1,j+1,k+1,1)+Ez(i+1,j+1,k+1,2)+Ez(i+1,j+1,k+1,3) & + -Ez(i ,j+1,k+1,1)-Ez(i ,j+1,k+1,2)-Ez(i ,j+1,k+1,3) & + +Ez(i+1,j-1,k+1,1)+Ez(i+1,j-1,k+1,2)+Ez(i+1,j-1,k+1,3) & + -Ez(i ,j-1,k+1,1)-Ez(i ,j-1,k+1,2)-Ez(i ,j-1,k+1,3) & + +Ez(i+1,j+1,k-1,1)+Ez(i+1,j+1,k-1,2)+Ez(i+1,j+1,k-1,3) & + -Ez(i ,j+1,k-1,1)-Ez(i ,j+1,k-1,2)-Ez(i ,j+1,k-1,3) & + +Ez(i+1,j-1,k-1,1)+Ez(i+1,j-1,k-1,2)+Ez(i+1,j-1,k-1,3) & + -Ez(i ,j-1,k-1,1)-Ez(i ,j-1,k-1,2)-Ez(i ,j-1,k-1,3))) end do end do end do @@ -202,43 +200,43 @@ contains do k = zlo(3), zhi(3) do j = zlo(2), zhi(2) do i = zlo(1), zhi(1) - Bz(i,j,k,1) = sigx1(i)*Bz(i,j,k,1) - sigx2(i)*(alphax*(Ey(i+1,j ,k ,1)+Ey(i+1,j ,k ,2)+Ey(i+1,j ,k ,3) & - -Ey(i ,j ,k ,1)-Ey(i ,j ,k ,2)-Ey(i ,j ,k ,3)) & - +betaxy*(Ey(i+1,j+1,k ,1)+Ey(i+1,j+1,k ,2)+Ey(i+1,j+1,k ,3) & - -Ey(i ,j+1,k ,1)-Ey(i ,j+1,k ,2)-Ey(i ,j+1,k ,3) & - +Ey(i+1,j-1,k ,1)+Ey(i+1,j-1,k ,2)+Ey(i+1,j-1,k ,3) & - -Ey(i ,j-1,k ,1)-Ey(i ,j-1,k ,2)-Ey(i ,j-1,k ,3)) & - +betaxz*(Ey(i+1,j ,k+1,1)+Ey(i+1,j ,k+1,2)+Ey(i+1,j ,k+1,3) & - -Ey(i ,j ,k+1,1)-Ey(i ,j ,k+1,2)-Ey(i ,j ,k+1,3) & - +Ey(i+1,j ,k-1,1)+Ey(i+1,j ,k-1,2)+Ey(i+1,j ,k-1,3) & - -Ey(i ,j ,k-1,1)-Ey(i ,j ,k-1,2)-Ey(i ,j ,k-1,3)) & - +gammax*(Ey(i+1,j+1,k+1,1)+Ey(i+1,j+1,k+1,2)+Ey(i+1,j+1,k+1,3) & - -Ey(i ,j+1,k+1,1)-Ey(i ,j+1,k+1,2)-Ey(i ,j+1,k+1,3) & - +Ey(i+1,j-1,k+1,1)+Ey(i+1,j-1,k+1,2)+Ey(i+1,j-1,k+1,3) & - -Ey(i ,j-1,k+1,1)-Ey(i ,j-1,k+1,2)-Ey(i ,j-1,k+1,3) & - +Ey(i+1,j+1,k-1,1)+Ey(i+1,j+1,k-1,2)+Ey(i+1,j+1,k-1,3) & - -Ey(i ,j+1,k-1,1)-Ey(i ,j+1,k-1,2)-Ey(i ,j+1,k-1,3) & - +Ey(i+1,j-1,k-1,1)+Ey(i+1,j-1,k-1,2)+Ey(i+1,j-1,k-1,3) & - -Ey(i ,j-1,k-1,1)-Ey(i ,j-1,k-1,2)-Ey(i ,j-1,k-1,3))) - - Bz(i,j,k,2) = sigy1(j)*Bz(i,j,k,2) + sigy2(j)*(alphay*(Ex(i ,j+1,k ,1)+Ex(i ,j+1,k ,2)+Ex(i ,j+1,k ,3) & - -Ex(i ,j ,k ,1)-Ex(i ,j ,k ,2)-Ex(i ,j ,k ,3)) & - +betayx*(Ex(i+1,j+1,k ,1)+Ex(i+1,j+1,k ,2)+Ex(i+1,j+1,k ,3) & - -Ex(i+1,j ,k ,1)-Ex(i+1,j ,k ,2)-Ex(i+1,j ,k ,3) & - +Ex(i-1,j+1,k ,1)+Ex(i-1,j+1,k ,2)+Ex(i-1,j+1,k ,3) & - -Ex(i-1,j ,k ,1)-Ex(i-1,j ,k ,2)-Ex(i-1,j ,k ,3)) & - +betayz*(Ex(i ,j+1,k+1,1)+Ex(i ,j+1,k+1,2)+Ex(i ,j+1,k+1,3) & - -Ex(i ,j ,k+1,1)-Ex(i ,j ,k+1,2)-Ex(i ,j ,k+1,3) & - +Ex(i ,j+1,k-1,1)+Ex(i ,j+1,k-1,2)+Ex(i ,j+1,k-1,3) & - -Ex(i ,j ,k-1,1)-Ex(i ,j ,k-1,2)-Ex(i ,j ,k-1,3)) & - +gammay*(Ex(i+1,j+1,k+1,1)+Ex(i+1,j+1,k+1,2)+Ex(i+1,j+1,k+1,3) & - -Ex(i+1,j ,k+1,1)-Ex(i+1,j ,k+1,2)-Ex(i+1,j ,k+1,3) & - +Ex(i-1,j+1,k+1,1)+Ex(i-1,j+1,k+1,2)+Ex(i-1,j+1,k+1,3) & - -Ex(i-1,j ,k+1,1)-Ex(i-1,j ,k+1,2)-Ex(i-1,j ,k+1,3) & - +Ex(i+1,j+1,k-1,1)+Ex(i+1,j+1,k-1,2)+Ex(i+1,j+1,k-1,3) & - -Ex(i+1,j ,k-1,1)-Ex(i+1,j ,k-1,2)-Ex(i+1,j ,k-1,3) & - +Ex(i-1,j+1,k-1,1)+Ex(i-1,j+1,k-1,2)+Ex(i-1,j+1,k-1,3) & - -Ex(i-1,j ,k-1,1)-Ex(i-1,j ,k-1,2)-Ex(i-1,j ,k-1,3))) + Bz(i,j,k,1) = Bz(i,j,k,1) - (alphax*(Ey(i+1,j ,k ,1)+Ey(i+1,j ,k ,2)+Ey(i+1,j ,k ,3) & + -Ey(i ,j ,k ,1)-Ey(i ,j ,k ,2)-Ey(i ,j ,k ,3)) & + +betaxy*(Ey(i+1,j+1,k ,1)+Ey(i+1,j+1,k ,2)+Ey(i+1,j+1,k ,3) & + -Ey(i ,j+1,k ,1)-Ey(i ,j+1,k ,2)-Ey(i ,j+1,k ,3) & + +Ey(i+1,j-1,k ,1)+Ey(i+1,j-1,k ,2)+Ey(i+1,j-1,k ,3) & + -Ey(i ,j-1,k ,1)-Ey(i ,j-1,k ,2)-Ey(i ,j-1,k ,3)) & + +betaxz*(Ey(i+1,j ,k+1,1)+Ey(i+1,j ,k+1,2)+Ey(i+1,j ,k+1,3) & + -Ey(i ,j ,k+1,1)-Ey(i ,j ,k+1,2)-Ey(i ,j ,k+1,3) & + +Ey(i+1,j ,k-1,1)+Ey(i+1,j ,k-1,2)+Ey(i+1,j ,k-1,3) & + -Ey(i ,j ,k-1,1)-Ey(i ,j ,k-1,2)-Ey(i ,j ,k-1,3)) & + +gammax*(Ey(i+1,j+1,k+1,1)+Ey(i+1,j+1,k+1,2)+Ey(i+1,j+1,k+1,3) & + -Ey(i ,j+1,k+1,1)-Ey(i ,j+1,k+1,2)-Ey(i ,j+1,k+1,3) & + +Ey(i+1,j-1,k+1,1)+Ey(i+1,j-1,k+1,2)+Ey(i+1,j-1,k+1,3) & + -Ey(i ,j-1,k+1,1)-Ey(i ,j-1,k+1,2)-Ey(i ,j-1,k+1,3) & + +Ey(i+1,j+1,k-1,1)+Ey(i+1,j+1,k-1,2)+Ey(i+1,j+1,k-1,3) & + -Ey(i ,j+1,k-1,1)-Ey(i ,j+1,k-1,2)-Ey(i ,j+1,k-1,3) & + +Ey(i+1,j-1,k-1,1)+Ey(i+1,j-1,k-1,2)+Ey(i+1,j-1,k-1,3) & + -Ey(i ,j-1,k-1,1)-Ey(i ,j-1,k-1,2)-Ey(i ,j-1,k-1,3))) + + Bz(i,j,k,2) = Bz(i,j,k,2) + (alphay*(Ex(i ,j+1,k ,1)+Ex(i ,j+1,k ,2)+Ex(i ,j+1,k ,3) & + -Ex(i ,j ,k ,1)-Ex(i ,j ,k ,2)-Ex(i ,j ,k ,3)) & + +betayx*(Ex(i+1,j+1,k ,1)+Ex(i+1,j+1,k ,2)+Ex(i+1,j+1,k ,3) & + -Ex(i+1,j ,k ,1)-Ex(i+1,j ,k ,2)-Ex(i+1,j ,k ,3) & + +Ex(i-1,j+1,k ,1)+Ex(i-1,j+1,k ,2)+Ex(i-1,j+1,k ,3) & + -Ex(i-1,j ,k ,1)-Ex(i-1,j ,k ,2)-Ex(i-1,j ,k ,3)) & + +betayz*(Ex(i ,j+1,k+1,1)+Ex(i ,j+1,k+1,2)+Ex(i ,j+1,k+1,3) & + -Ex(i ,j ,k+1,1)-Ex(i ,j ,k+1,2)-Ex(i ,j ,k+1,3) & + +Ex(i ,j+1,k-1,1)+Ex(i ,j+1,k-1,2)+Ex(i ,j+1,k-1,3) & + -Ex(i ,j ,k-1,1)-Ex(i ,j ,k-1,2)-Ex(i ,j ,k-1,3)) & + +gammay*(Ex(i+1,j+1,k+1,1)+Ex(i+1,j+1,k+1,2)+Ex(i+1,j+1,k+1,3) & + -Ex(i+1,j ,k+1,1)-Ex(i+1,j ,k+1,2)-Ex(i+1,j ,k+1,3) & + +Ex(i-1,j+1,k+1,1)+Ex(i-1,j+1,k+1,2)+Ex(i-1,j+1,k+1,3) & + -Ex(i-1,j ,k+1,1)-Ex(i-1,j ,k+1,2)-Ex(i-1,j ,k+1,3) & + +Ex(i+1,j+1,k-1,1)+Ex(i+1,j+1,k-1,2)+Ex(i+1,j+1,k-1,3) & + -Ex(i+1,j ,k-1,1)-Ex(i+1,j ,k-1,2)-Ex(i+1,j ,k-1,3) & + +Ex(i-1,j+1,k-1,1)+Ex(i-1,j+1,k-1,2)+Ex(i-1,j+1,k-1,3) & + -Ex(i-1,j ,k-1,1)-Ex(i-1,j ,k-1,2)-Ex(i-1,j ,k-1,3))) end do end do end do @@ -254,41 +252,28 @@ contains & Bx, Bxlo, Bxhi, & & By, Bylo, Byhi, & & Bz, Bzlo, Bzhi, & - & sigx1, sigx1_lo, sigx1_hi, & - & sigx2, sigx2_lo, sigx2_hi, & - & sigy1, sigy1_lo, sigy1_hi, & - & sigy2, sigy2_lo, sigy2_hi, & - & sigz1, sigz1_lo, sigz1_hi, & - & sigz2, sigz2_lo, sigz2_hi) & + & dtsdx, dtsdy, dtsdz) & bind(c,name='warpx_push_pml_evec_3d') integer, intent(in) :: xlo(3), xhi(3), ylo(3), yhi(3), zlo(3), zhi(3), & Exlo(3), Exhi(3), Eylo(3), Eyhi(3), Ezlo(3), Ezhi(3), & Bxlo(3), Bxhi(3), Bylo(3), Byhi(3), Bzlo(3), Bzhi(3) - integer, intent(in), value :: sigx1_lo, sigx1_hi, sigx2_lo, sigx2_hi - integer, intent(in), value :: sigy1_lo, sigy1_hi, sigy2_lo, sigy2_hi - integer, intent(in), value :: sigz1_lo, sigz1_hi, sigz2_lo, sigz2_hi + real(amrex_real), intent(in) :: dtsdx, dtsdy, dtsdz real(amrex_real), intent(inout) :: Ex (Exlo(1):Exhi(1),Exlo(2):Exhi(2),Exlo(3):Exhi(3),2) real(amrex_real), intent(inout) :: Ey (Eylo(1):Eyhi(1),Eylo(2):Eyhi(2),Eylo(3):Eyhi(3),2) real(amrex_real), intent(inout) :: Ez (Ezlo(1):Ezhi(1),Ezlo(2):Ezhi(2),Ezlo(3):Ezhi(3),2) real(amrex_real), intent(in ) :: Bx (Bxlo(1):Bxhi(1),Bxlo(2):Bxhi(2),Bxlo(3):Bxhi(3),2) real(amrex_real), intent(in ) :: By (Bylo(1):Byhi(1),Bylo(2):Byhi(2),Bylo(3):Byhi(3),2) real(amrex_real), intent(in ) :: Bz (Bzlo(1):Bzhi(1),Bzlo(2):Bzhi(2),Bzlo(3):Bzhi(3),2) - real(amrex_real), intent(in) :: sigx1(sigx1_lo:sigx1_hi) - real(amrex_real), intent(in) :: sigx2(sigx2_lo:sigx2_hi) - real(amrex_real), intent(in) :: sigy1(sigy1_lo:sigy1_hi) - real(amrex_real), intent(in) :: sigy2(sigy2_lo:sigy2_hi) - real(amrex_real), intent(in) :: sigz1(sigz1_lo:sigz1_hi) - real(amrex_real), intent(in) :: sigz2(sigz2_lo:sigz2_hi) integer :: i, j, k do k = xlo(3), xhi(3) do j = xlo(2), xhi(2) do i = xlo(1), xhi(1) - Ex(i,j,k,1) = sigy1(j)*Ex(i,j,k,1) + sigy2(j)*(Bz(i,j ,k ,1)+Bz(i,j ,k ,2) & - & -Bz(i,j-1,k ,1)-Bz(i,j-1,k ,2)) - Ex(i,j,k,2) = sigz1(k)*Ex(i,j,k,2) - sigz2(k)*(By(i,j ,k ,1)+By(i,j ,k ,2) & - & -By(i,j ,k-1,1)-By(i,j ,k-1,2)) + Ex(i,j,k,1) = Ex(i,j,k,1) + dtsdy*(Bz(i,j ,k ,1)+Bz(i,j ,k ,2) & + & -Bz(i,j-1,k ,1)-Bz(i,j-1,k ,2)) + Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & + & -By(i,j ,k-1,1)-By(i,j ,k-1,2)) end do end do end do @@ -296,10 +281,10 @@ contains do k = ylo(3), yhi(3) do j = ylo(2), yhi(2) do i = ylo(1), yhi(1) - Ey(i,j,k,1) = sigz1(k)*Ey(i,j,k,1) + sigz2(k)*(Bx(i ,j,k ,1)+Bx(i ,j,k ,2) & - & -Bx(i ,j,k-1,1)-Bx(i ,j,k-1,2)) - Ey(i,j,k,2) = sigx1(i)*Ey(i,j,k,2) - sigx2(i)*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & - & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2)) + Ey(i,j,k,1) = Ey(i,j,k,1) + dtsdz*(Bx(i ,j,k ,1)+Bx(i ,j,k ,2) & + & -Bx(i ,j,k-1,1)-Bx(i ,j,k-1,2)) + Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & + & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2)) end do end do end do @@ -307,10 +292,10 @@ contains do k = zlo(3), zhi(3) do j = zlo(2), zhi(2) do i = zlo(1), zhi(1) - Ez(i,j,k,1) = sigx1(i)*Ez(i,j,k,1) + sigx2(i)*(By(i ,j ,k,1)+By(i ,j ,k,2) & - & -By(i-1,j ,k,1)-By(i-1,j ,k,2)) - Ez(i,j,k,2) = sigy1(j)*Ez(i,j,k,2) - sigy2(j)*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & - & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2)) + Ez(i,j,k,1) = Ez(i,j,k,1) + dtsdx*(By(i ,j ,k,1)+By(i ,j ,k,2) & + & -By(i-1,j ,k,1)-By(i-1,j ,k,2)) + Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & + & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2)) end do end do end do @@ -324,27 +309,17 @@ contains & Bx, Bxlo, Bxhi, & & By, Bylo, Byhi, & & Bz, Bzlo, Bzhi, & - & sigx1, sigx1_lo, sigx1_hi, & - & sigx2, sigx2_lo, sigx2_hi, & - & sigz1, sigz1_lo, sigz1_hi, & - & sigz2, sigz2_lo, sigz2_hi, & & dtsdx, dtsdy, dtsdz, solver_type) & bind(c,name='warpx_push_pml_bvec_2d') integer, intent(in) :: xlo(2), xhi(2), ylo(2), yhi(2), zlo(2), zhi(2), & Exlo(2), Exhi(2), Eylo(2), Eyhi(2), Ezlo(2), Ezhi(2), & Bxlo(2), Bxhi(2), Bylo(2), Byhi(2), Bzlo(2), Bzhi(2), solver_type - integer, intent(in), value :: sigx1_lo, sigx1_hi, sigx2_lo, sigx2_hi - integer, intent(in), value :: sigz1_lo, sigz1_hi, sigz2_lo, sigz2_hi real(amrex_real), intent(in ) :: Ex (Exlo(1):Exhi(1),Exlo(2):Exhi(2),3) real(amrex_real), intent(in ) :: Ey (Eylo(1):Eyhi(1),Eylo(2):Eyhi(2),3) real(amrex_real), intent(in ) :: Ez (Ezlo(1):Ezhi(1),Ezlo(2):Ezhi(2),3) real(amrex_real), intent(inout) :: Bx (Bxlo(1):Bxhi(1),Bxlo(2):Bxhi(2),2) real(amrex_real), intent(inout) :: By (Bylo(1):Byhi(1),Bylo(2):Byhi(2),2) real(amrex_real), intent(inout) :: Bz (Bzlo(1):Bzhi(1),Bzlo(2):Bzhi(2),2) - real(amrex_real), intent(in) :: sigx1(sigx1_lo:sigx1_hi) - real(amrex_real), intent(in) :: sigx2(sigx2_lo:sigx2_hi) - real(amrex_real), intent(in) :: sigz1(sigz1_lo:sigz1_hi) - real(amrex_real), intent(in) :: sigz2(sigz2_lo:sigz2_hi) real(amrex_real), intent(in) :: dtsdx, dtsdy, dtsdz integer :: i, k @@ -360,24 +335,24 @@ contains do k = xlo(2), xhi(2) do i = xlo(1), xhi(1) - Bx(i,k,2) = sigz1(k)*Bx(i,k,2) + sigz2(k)*(Ey(i,k+1,1)+Ey(i,k+1,2)+Ey(i,k+1,3) & - & -Ey(i,k ,1)-Ey(i,k ,2)-Ey(i,k ,3)) + Bx(i,k,2) = Bx(i,k,2) + dtsdz*(Ey(i,k+1,1)+Ey(i,k+1,2)+Ey(i,k+1,3) & + & -Ey(i,k ,1)-Ey(i,k ,2)-Ey(i,k ,3)) end do end do do k = ylo(2), yhi(2) do i = ylo(1), yhi(1) - By(i,k,1) = sigz1(k)*By(i,k,1) - sigz2(k)*(Ex(i ,k+1,1)+Ex(i ,k+1,2)+Ex(i ,k+1,3) & - & -Ex(i ,k ,1)-Ex(i ,k ,2)-Ex(i ,k ,3)) - By(i,k,2) = sigx1(i)*By(i,k,2) + sigx2(i)*(Ez(i+1,k ,1)+Ez(i+1,k ,2)+Ez(i+1,k ,3) & - & -Ez(i ,k ,1)-Ez(i ,k ,2)-Ez(i ,k ,3)) + By(i,k,1) = By(i,k,1) - dtsdz*(Ex(i ,k+1,1)+Ex(i ,k+1,2)+Ex(i ,k+1,3) & + & -Ex(i ,k ,1)-Ex(i ,k ,2)-Ex(i ,k ,3)) + By(i,k,2) = By(i,k,2) + dtsdx*(Ez(i+1,k ,1)+Ez(i+1,k ,2)+Ez(i+1,k ,3) & + & -Ez(i ,k ,1)-Ez(i ,k ,2)-Ez(i ,k ,3)) end do end do do k = zlo(2), zhi(2) do i = zlo(1), zhi(1) - Bz(i,k,1) = sigx1(i)*Bz(i,k,1) - sigx2(i)*(Ey(i+1,k,1)+Ey(i+1,k,2)+Ey(i+1,k,3) & - & -Ey(i ,k,1)-Ey(i ,k,2)-Ey(i ,k,3)) + Bz(i,k,1) = Bz(i,k,1) - dtsdx*(Ey(i+1,k,1)+Ey(i+1,k,2)+Ey(i+1,k,3) & + & -Ey(i ,k,1)-Ey(i ,k,2)-Ey(i ,k,3)) end do end do @@ -394,45 +369,49 @@ contains alphax = 1. - 2.*betaxz alphaz = 1. - 2.*betazx + betaxz = dtsdx*betaxz + betazx = dtsdz*betazx + alphax = dtsdx*alphax + alphaz = dtsdz*alphaz do k = xlo(2), xhi(2) do i = xlo(1), xhi(1) - Bx(i,k,2) = sigz1(k)*Bx(i,k,2) + sigz2(k)*( alphaz*(Ey(i ,k+1,1)+Ey(i ,k+1,2)+Ey(i ,k+1,3) & - -Ey(i ,k ,1)-Ey(i ,k ,2)-Ey(i ,k ,3)) & - + betazx*(Ey(i+1,k+1,1)+Ey(i+1,k+1,2)+Ey(i+1,k+1,3) & - -Ey(i+1,k ,1)-Ey(i+1,k ,2)-Ey(i+1,k ,3) & - +Ey(i-1,k+1,1)+Ey(i-1,k+1,2)+Ey(i-1,k+1,3) & - -Ey(i-1,k ,1)-Ey(i-1,k ,2)-Ey(i-1,k ,3))) + Bx(i,k,2) = Bx(i,k,2) + ( alphaz*(Ey(i ,k+1,1)+Ey(i ,k+1,2)+Ey(i ,k+1,3) & + -Ey(i ,k ,1)-Ey(i ,k ,2)-Ey(i ,k ,3)) & + + betazx*(Ey(i+1,k+1,1)+Ey(i+1,k+1,2)+Ey(i+1,k+1,3) & + -Ey(i+1,k ,1)-Ey(i+1,k ,2)-Ey(i+1,k ,3) & + +Ey(i-1,k+1,1)+Ey(i-1,k+1,2)+Ey(i-1,k+1,3) & + -Ey(i-1,k ,1)-Ey(i-1,k ,2)-Ey(i-1,k ,3))) end do end do do k = ylo(2), yhi(2) do i = ylo(1), yhi(1) - By(i,k,1) = sigz1(k)*By(i,k,1) - sigz2(k)*( alphaz*(Ex(i ,k+1,1)+Ex(i ,k+1,2)+Ex(i ,k+1,3) & - -Ex(i ,k ,1)-Ex(i ,k ,2)-Ex(i ,k ,3)) & - + betazx*(Ex(i+1,k+1,1)+Ex(i+1,k+1,2)+Ex(i+1,k+1,3) & - -Ex(i+1,k ,1)-Ex(i+1,k ,2)-Ex(i+1,k ,3) & - +Ex(i-1,k+1,1)+Ex(i-1,k+1,2)+Ex(i-1,k+1,3) & - -Ex(i-1,k ,1)-Ex(i-1,k ,2)-Ex(i-1,k ,3))) - - By(i,k,2) = sigx1(i)*By(i,k,2) + sigx2(i)*( alphax*(Ez(i+1,k ,1)+Ez(i+1,k ,2)+Ez(i+1,k ,3) & - -Ez(i ,k ,1)-Ez(i ,k ,2)-Ez(i ,k ,3)) & - + betaxz*(Ez(i+1,k+1,1)+Ez(i+1,k+1,2)+Ez(i+1,k+1,3) & - -Ez(i ,k+1,1)-Ez(i ,k+1,2)-Ez(i ,k+1,3) & - +Ez(i+1,k-1,1)+Ez(i+1,k-1,2)+Ez(i+1,k-1,3) & - -Ez(i ,k-1,1)-Ez(i ,k-1,2)-Ez(i ,k-1,3))) + By(i,k,1) = By(i,k,1) - ( alphaz*(Ex(i ,k+1,1)+Ex(i ,k+1,2)+Ex(i ,k+1,3) & + -Ex(i ,k ,1)-Ex(i ,k ,2)-Ex(i ,k ,3)) & + + betazx*(Ex(i+1,k+1,1)+Ex(i+1,k+1,2)+Ex(i+1,k+1,3) & + -Ex(i+1,k ,1)-Ex(i+1,k ,2)-Ex(i+1,k ,3) & + +Ex(i-1,k+1,1)+Ex(i-1,k+1,2)+Ex(i-1,k+1,3) & + -Ex(i-1,k ,1)-Ex(i-1,k ,2)-Ex(i-1,k ,3))) + + By(i,k,2) = By(i,k,2) + ( alphax*(Ez(i+1,k ,1)+Ez(i+1,k ,2)+Ez(i+1,k ,3) & + -Ez(i ,k ,1)-Ez(i ,k ,2)-Ez(i ,k ,3)) & + + betaxz*(Ez(i+1,k+1,1)+Ez(i+1,k+1,2)+Ez(i+1,k+1,3) & + -Ez(i ,k+1,1)-Ez(i ,k+1,2)-Ez(i ,k+1,3) & + +Ez(i+1,k-1,1)+Ez(i+1,k-1,2)+Ez(i+1,k-1,3) & + -Ez(i ,k-1,1)-Ez(i ,k-1,2)-Ez(i ,k-1,3))) end do end do do k = zlo(2), zhi(2) do i = zlo(1), zhi(1) - Bz(i,k,1) = sigx1(i)*Bz(i,k,1) - sigx2(i)*( alphax*(Ey(i+1,k ,1)+Ey(i+1,k ,2)+Ey(i+1,k ,3) & - -Ey(i ,k ,1)-Ey(i ,k ,2)-Ey(i ,k ,3)) & - + betaxz*(Ey(i+1,k+1,1)+Ey(i+1,k+1,2)+Ey(i+1,k+1,3) & - -Ey(i ,k+1,1)-Ey(i ,k+1,2)-Ey(i ,k+1,3) & - +Ey(i+1,k-1,1)+Ey(i+1,k-1,2)+Ey(i+1,k-1,3) & - -Ey(i ,k-1,1)-Ey(i ,k-1,2)-Ey(i ,k-1,3))) + Bz(i,k,1) = Bz(i,k,1) - ( alphax*(Ey(i+1,k ,1)+Ey(i+1,k ,2)+Ey(i+1,k ,3) & + -Ey(i ,k ,1)-Ey(i ,k ,2)-Ey(i ,k ,3)) & + + betaxz*(Ey(i+1,k+1,1)+Ey(i+1,k+1,2)+Ey(i+1,k+1,3) & + -Ey(i ,k+1,1)-Ey(i ,k+1,2)-Ey(i ,k+1,3) & + +Ey(i+1,k-1,1)+Ey(i+1,k-1,2)+Ey(i+1,k-1,3) & + -Ey(i ,k-1,1)-Ey(i ,k-1,2)-Ey(i ,k-1,3))) end do end do @@ -448,49 +427,41 @@ contains & Bx, Bxlo, Bxhi, & & By, Bylo, Byhi, & & Bz, Bzlo, Bzhi, & - & sigx1, sigx1_lo, sigx1_hi, & - & sigx2, sigx2_lo, sigx2_hi, & - & sigz1, sigz1_lo, sigz1_hi, & - & sigz2, sigz2_lo, sigz2_hi) & + & dtsdx, dtsdy, dtsdz) & bind(c,name='warpx_push_pml_evec_2d') integer, intent(in) :: xlo(2), xhi(2), ylo(2), yhi(2), zlo(2), zhi(2), & Exlo(2), Exhi(2), Eylo(2), Eyhi(2), Ezlo(2), Ezhi(2), & Bxlo(2), Bxhi(2), Bylo(2), Byhi(2), Bzlo(2), Bzhi(2) - integer, intent(in), value :: sigx1_lo, sigx1_hi, sigx2_lo, sigx2_hi - integer, intent(in), value :: sigz1_lo, sigz1_hi, sigz2_lo, sigz2_hi + real(amrex_real), intent(in) :: dtsdx, dtsdy, dtsdz real(amrex_real), intent(inout) :: Ex (Exlo(1):Exhi(1),Exlo(2):Exhi(2),2) real(amrex_real), intent(inout) :: Ey (Eylo(1):Eyhi(1),Eylo(2):Eyhi(2),2) real(amrex_real), intent(inout) :: Ez (Ezlo(1):Ezhi(1),Ezlo(2):Ezhi(2),2) real(amrex_real), intent(in ) :: Bx (Bxlo(1):Bxhi(1),Bxlo(2):Bxhi(2),2) real(amrex_real), intent(in ) :: By (Bylo(1):Byhi(1),Bylo(2):Byhi(2),2) real(amrex_real), intent(in ) :: Bz (Bzlo(1):Bzhi(1),Bzlo(2):Bzhi(2),2) - real(amrex_real), intent(in) :: sigx1(sigx1_lo:sigx1_hi) - real(amrex_real), intent(in) :: sigx2(sigx2_lo:sigx2_hi) - real(amrex_real), intent(in) :: sigz1(sigz1_lo:sigz1_hi) - real(amrex_real), intent(in) :: sigz2(sigz2_lo:sigz2_hi) integer :: i, k do k = xlo(2), xhi(2) do i = xlo(1), xhi(1) - Ex(i,k,2) = sigz1(k)*Ex(i,k,2) - sigz2(k)*(By(i,k ,1)+By(i,k ,2) & - & -By(i,k-1,1)-By(i,k-1,2)) + Ex(i,k,2) = Ex(i,k,2) - dtsdz*(By(i,k ,1)+By(i,k ,2) & + & -By(i,k-1,1)-By(i,k-1,2)) end do end do do k = ylo(2), yhi(2) do i = ylo(1), yhi(1) - Ey(i,k,1) = sigz1(k)*Ey(i,k,1) + sigz2(k)*(Bx(i ,k ,1)+Bx(i ,k ,2) & - & -Bx(i ,k-1,1)-Bx(i ,k-1,2)) - Ey(i,k,2) = sigx1(i)*Ey(i,k,2) - sigx2(i)*(Bz(i ,k ,1)+Bz(i ,k ,2) & - & -Bz(i-1,k ,1)-Bz(i-1,k ,2)) + Ey(i,k,1) = Ey(i,k,1) + dtsdz*(Bx(i ,k ,1)+Bx(i ,k ,2) & + & -Bx(i ,k-1,1)-Bx(i ,k-1,2)) + Ey(i,k,2) = Ey(i,k,2) - dtsdx*(Bz(i ,k ,1)+Bz(i ,k ,2) & + & -Bz(i-1,k ,1)-Bz(i-1,k ,2)) end do end do do k = zlo(2), zhi(2) do i = zlo(1), zhi(1) - Ez(i,k,1) = sigx1(i)*Ez(i,k,1) + sigx2(i)*(By(i ,k,1)+By(i ,k,2) & - & -By(i-1,k,1)-By(i-1,k,2)) + Ez(i,k,1) = Ez(i,k,1) + dtsdx*(By(i ,k,1)+By(i ,k,2) & + & -By(i-1,k,1)-By(i-1,k,2)) end do end do @@ -703,4 +674,152 @@ contains end subroutine warpx_push_pml_evec_f_2d + + subroutine warpx_damp_pml_2d (texlo, texhi, teylo, teyhi, tezlo, tezhi, & + & tbxlo, tbxhi, tbylo, tbyhi, tbzlo, tbzhi, & + & ex, exlo, exhi, ey, eylo, eyhi, ez, ezlo, ezhi, & + & bx, bxlo, bxhi, by, bylo, byhi, bz, bzlo, bzhi, & + & sigex, sexlo, sexhi, sigez, sezlo, sezhi, & + & sigbx, sbxlo, sbxhi, sigbz, sbzlo, sbzhi) & + bind(c,name='warpx_damp_pml_2d') + integer, dimension(2), intent(in) :: texlo, texhi, teylo, teyhi, tezlo, tezhi, & + tbxlo, tbxhi, tbylo, tbyhi, tbzlo, tbzhi, & + exlo, exhi, eylo, eyhi, ezlo, ezhi, bxlo, bxhi, bylo, byhi, bzlo, bzhi + integer, intent(in), value :: sexlo, sexhi, sezlo, sezhi, & + & sbxlo, sbxhi, sbzlo, sbzhi + real(amrex_real), intent(inout) :: ex(exlo(1):exhi(1),exlo(2):exhi(2),3) + real(amrex_real), intent(inout) :: ey(eylo(1):eyhi(1),eylo(2):eyhi(2),3) + real(amrex_real), intent(inout) :: ez(ezlo(1):ezhi(1),ezlo(2):ezhi(2),3) + real(amrex_real), intent(inout) :: bx(bxlo(1):bxhi(1),bxlo(2):bxhi(2),2) + real(amrex_real), intent(inout) :: by(bylo(1):byhi(1),bylo(2):byhi(2),2) + real(amrex_real), intent(inout) :: bz(bzlo(1):bzhi(1),bzlo(2):bzhi(2),2) + real(amrex_real), intent(in) :: sigex(sexlo:sexhi) + real(amrex_real), intent(in) :: sigez(sezlo:sezhi) + real(amrex_real), intent(in) :: sigbx(sbxlo:sbxhi) + real(amrex_real), intent(in) :: sigbz(sbzlo:sbzhi) + + integer :: i,k + + do k = texlo(2), texhi(2) + do i = texlo(1), texhi(1) + ex(i,k,2) = ex(i,k,2) * sigez(k) + end do + end do + + do k = teylo(2), teyhi(2) + do i = teylo(1), teyhi(1) + ey(i,k,1) = ey(i,k,1) * sigez(k) + ey(i,k,2) = ey(i,k,2) * sigex(i) + end do + end do + + do k = tezlo(2), tezhi(2) + do i = tezlo(1), tezhi(1) + ez(i,k,1) = ez(i,k,1) * sigex(i) + end do + end do + + do k = tbxlo(2), tbxhi(2) + do i = tbxlo(1), tbxhi(1) + bx(i,k,2) = bx(i,k,2) * sigbz(k) + end do + end do + + do k = tbylo(2), tbyhi(2) + do i = tbylo(1), tbyhi(1) + by(i,k,1) = by(i,k,1) * sigbz(k) + by(i,k,2) = by(i,k,2) * sigbx(i) + end do + end do + + do k = tbzlo(2), tbzhi(2) + do i = tbzlo(1), tbzhi(1) + bz(i,k,1) = bz(i,k,1) * sigbx(i) + end do + end do + end subroutine warpx_damp_pml_2d + + + subroutine warpx_damp_pml_3d (texlo, texhi, teylo, teyhi, tezlo, tezhi, & + & tbxlo, tbxhi, tbylo, tbyhi, tbzlo, tbzhi, & + & ex, exlo, exhi, ey, eylo, eyhi, ez, ezlo, ezhi, & + & bx, bxlo, bxhi, by, bylo, byhi, bz, bzlo, bzhi, & + & sigex, sexlo, sexhi, sigey, seylo, seyhi, sigez, sezlo, sezhi, & + & sigbx, sbxlo, sbxhi, sigby, sbylo, sbyhi, sigbz, sbzlo, sbzhi) & + bind(c,name='warpx_damp_pml_3d') + integer, dimension(3), intent(in) :: texlo, texhi, teylo, teyhi, tezlo, tezhi, & + tbxlo, tbxhi, tbylo, tbyhi, tbzlo, tbzhi, & + exlo, exhi, eylo, eyhi, ezlo, ezhi, bxlo, bxhi, bylo, byhi, bzlo, bzhi + integer, intent(in), value :: sexlo, sexhi, seylo, seyhi, sezlo, sezhi, & + & sbxlo, sbxhi, sbylo, sbyhi, sbzlo, sbzhi + real(amrex_real), intent(inout) :: ex(exlo(1):exhi(1),exlo(2):exhi(2),exlo(3):exhi(3),3) + real(amrex_real), intent(inout) :: ey(eylo(1):eyhi(1),eylo(2):eyhi(2),eylo(3):eyhi(3),3) + real(amrex_real), intent(inout) :: ez(ezlo(1):ezhi(1),ezlo(2):ezhi(2),ezlo(3):ezhi(3),3) + real(amrex_real), intent(inout) :: bx(bxlo(1):bxhi(1),bxlo(2):bxhi(2),bxlo(3):bxhi(3),2) + real(amrex_real), intent(inout) :: by(bylo(1):byhi(1),bylo(2):byhi(2),bylo(3):byhi(3),2) + real(amrex_real), intent(inout) :: bz(bzlo(1):bzhi(1),bzlo(2):bzhi(2),bzlo(3):bzhi(3),2) + real(amrex_real), intent(in) :: sigex(sexlo:sexhi) + real(amrex_real), intent(in) :: sigey(seylo:seyhi) + real(amrex_real), intent(in) :: sigez(sezlo:sezhi) + real(amrex_real), intent(in) :: sigbx(sbxlo:sbxhi) + real(amrex_real), intent(in) :: sigby(sbylo:sbyhi) + real(amrex_real), intent(in) :: sigbz(sbzlo:sbzhi) + + integer :: i,j,k + + do k = texlo(3), texhi(3) + do j = texlo(2), texhi(2) + do i = texlo(1), texhi(1) + ex(i,j,k,1) = ex(i,j,k,1) * sigey(j) + ex(i,j,k,2) = ex(i,j,k,2) * sigez(k) + end do + end do + end do + + do k = teylo(3), teyhi(3) + do j = teylo(2), teyhi(2) + do i = teylo(1), teyhi(1) + ey(i,j,k,1) = ey(i,j,k,1) * sigez(k) + ey(i,j,k,2) = ey(i,j,k,2) * sigex(i) + end do + end do + end do + + do k = tezlo(3), tezhi(3) + do j = tezlo(2), tezhi(2) + do i = tezlo(1), tezhi(1) + ez(i,j,k,1) = ez(i,j,k,1) * sigex(i) + ez(i,j,k,2) = ez(i,j,k,2) * sigey(j) + end do + end do + end do + + do k = tbxlo(3), tbxhi(3) + do j = tbxlo(2), tbxhi(2) + do i = tbxlo(1), tbxhi(1) + bx(i,j,k,1) = bx(i,j,k,1) * sigby(j) + bx(i,j,k,2) = bx(i,j,k,2) * sigbz(k) + end do + end do + end do + + do k = tbylo(3), tbyhi(3) + do j = tbylo(2), tbyhi(2) + do i = tbylo(1), tbyhi(1) + by(i,j,k,1) = by(i,j,k,1) * sigbz(k) + by(i,j,k,2) = by(i,j,k,2) * sigbx(i) + end do + end do + end do + + do k = tbzlo(3), tbzhi(3) + do j = tbzlo(2), tbzhi(2) + do i = tbzlo(1), tbzhi(1) + bz(i,j,k,1) = bz(i,j,k,1) * sigbx(i) + bz(i,j,k,2) = bz(i,j,k,2) * sigby(j) + end do + end do + end do + end subroutine warpx_damp_pml_3d + end module warpx_pml_module -- cgit v1.2.3 From 2dca5b21cdcba99100aaab94096925105b611e8a Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Tue, 7 Aug 2018 18:35:47 -0700 Subject: change the number of components of rho to 2 --- Source/LaserParticleContainer.H | 3 +- Source/LaserParticleContainer.cpp | 13 +++-- Source/ParticleContainer.H | 2 +- Source/ParticleContainer.cpp | 5 +- Source/PhysicalParticleContainer.H | 1 - Source/PhysicalParticleContainer.cpp | 13 +++-- Source/RigidInjectedParticleContainer.H | 1 - Source/RigidInjectedParticleContainer.cpp | 6 +-- Source/WarpX.H | 9 +--- Source/WarpX.cpp | 35 +++++++------ Source/WarpXComm.cpp | 13 +++-- Source/WarpXEvolve.cpp | 11 +---- Source/WarpXFFT.cpp | 44 +++-------------- Source/WarpXMove.cpp | 2 +- Source/WarpXParticleContainer.H | 3 +- Source/WarpXProb.cpp | 6 +-- Source/WarpXRegrid.cpp | 10 ++-- Source/WarpX_f.F90 | 82 ++++++++++++++++--------------- Source/WarpX_f.H | 6 ++- Source/main.cpp | 2 +- 20 files changed, 113 insertions(+), 154 deletions(-) (limited to 'Source/WarpXEvolve.cpp') diff --git a/Source/LaserParticleContainer.H b/Source/LaserParticleContainer.H index 2a6c90eff..79fc5e4fa 100644 --- a/Source/LaserParticleContainer.H +++ b/Source/LaserParticleContainer.H @@ -27,8 +27,7 @@ public: const amrex::MultiFab&, const amrex::MultiFab&, const amrex::MultiFab&, const amrex::MultiFab&, const amrex::MultiFab&, const amrex::MultiFab&, amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, - amrex::MultiFab* rho, amrex::MultiFab* rho2, - amrex::Real t, amrex::Real dt) final; + amrex::MultiFab* rho, amrex::Real t, amrex::Real dt) final; virtual void PushP (int lev, amrex::Real dt, const amrex::MultiFab& , diff --git a/Source/LaserParticleContainer.cpp b/Source/LaserParticleContainer.cpp index abda0cecf..971879d7c 100644 --- a/Source/LaserParticleContainer.cpp +++ b/Source/LaserParticleContainer.cpp @@ -276,8 +276,7 @@ LaserParticleContainer::Evolve (int lev, const MultiFab&, const MultiFab&, const MultiFab&, const MultiFab&, const MultiFab&, const MultiFab&, MultiFab& jx, MultiFab& jy, MultiFab& jz, - MultiFab* rho, MultiFab* rho2, - Real t, Real dt) + MultiFab* rho, Real t, Real dt) { BL_PROFILE("Laser::Evolve()"); BL_PROFILE_VAR_NS("Laser::Evolve::Copy", blp_copy); @@ -364,7 +363,7 @@ LaserParticleContainer::Evolve (int lev, long lvect = 8; - auto depositCharge = [&] (MultiFab* rhomf) + auto depositCharge = [&] (MultiFab* rhomf, int icomp) { long ngRho = rhomf->nGrow(); long ngRhoDeposit = (WarpX::use_filter) ? ngRho +1 : ngRho; @@ -420,16 +419,16 @@ LaserParticleContainer::Evolve (int lev, ncomp); amrex_atomic_accumulate_fab(BL_TO_FORTRAN_3D(filtered_rho), - BL_TO_FORTRAN_3D(rhofab), ncomp); + BL_TO_FORTRAN_N_3D(rhofab,icomp), ncomp); } else { amrex_atomic_accumulate_fab(BL_TO_FORTRAN_3D(local_rho), - BL_TO_FORTRAN_3D(rhofab), ncomp); + BL_TO_FORTRAN_N_3D(rhofab,icomp), ncomp); } }; - if (rho) depositCharge(rho); + if (rho) depositCharge(rho,0); // // Particle Push @@ -609,7 +608,7 @@ LaserParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_pxr_cd); - if (rho2) depositCharge(rho2); + if (rho) depositCharge(rho,1); // // copy particle data back diff --git a/Source/ParticleContainer.H b/Source/ParticleContainer.H index a92ded4ea..1f582f677 100644 --- a/Source/ParticleContainer.H +++ b/Source/ParticleContainer.H @@ -97,7 +97,7 @@ public: const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz, amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, - amrex::MultiFab* rho, amrex::MultiFab* rho2, amrex::Real t, amrex::Real dt); + amrex::MultiFab* rho, amrex::Real t, amrex::Real dt); /// /// This pushes the particle positions by one half time step for all the species in the diff --git a/Source/ParticleContainer.cpp b/Source/ParticleContainer.cpp index e0539e799..eeed71fce 100644 --- a/Source/ParticleContainer.cpp +++ b/Source/ParticleContainer.cpp @@ -172,15 +172,14 @@ MultiParticleContainer::Evolve (int lev, const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz, MultiFab& jx, MultiFab& jy, MultiFab& jz, - MultiFab* rho, MultiFab* rho2, Real t, Real dt) + MultiFab* rho, Real t, Real dt) { jx.setVal(0.0); jy.setVal(0.0); jz.setVal(0.0); if (rho) rho->setVal(0.0); - if (rho2) rho2->setVal(0.0); for (auto& pc : allcontainers) { - pc->Evolve(lev, Ex, Ey, Ez, Bx, By, Bz, jx, jy, jz, rho, rho2, t, dt); + pc->Evolve(lev, Ex, Ey, Ez, Bx, By, Bz, jx, jy, jz, rho, t, dt); } } diff --git a/Source/PhysicalParticleContainer.H b/Source/PhysicalParticleContainer.H index 183d4e43a..d9af3f20c 100644 --- a/Source/PhysicalParticleContainer.H +++ b/Source/PhysicalParticleContainer.H @@ -45,7 +45,6 @@ public: amrex::MultiFab& jy, amrex::MultiFab& jz, amrex::MultiFab* rho, - amrex::MultiFab* rho2, amrex::Real t, amrex::Real dt) override; diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 2e4ff73ba..4251f2db1 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -638,8 +638,7 @@ PhysicalParticleContainer::Evolve (int lev, const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz, MultiFab& jx, MultiFab& jy, MultiFab& jz, - MultiFab* rho, MultiFab* rho2, - Real t, Real dt) + MultiFab* rho, Real t, Real dt) { BL_PROFILE("PPC::Evolve()"); BL_PROFILE_VAR_NS("PPC::Evolve::Copy", blp_copy); @@ -789,7 +788,7 @@ PhysicalParticleContainer::Evolve (int lev, long lvect = 8; - auto depositCharge = [&] (MultiFab* rhomf) + auto depositCharge = [&] (MultiFab* rhomf, int icomp) { long ngRho = rhomf->nGrow(); long ngRhoDeposit = (WarpX::use_filter) ? ngRho +1 : ngRho; @@ -845,16 +844,16 @@ PhysicalParticleContainer::Evolve (int lev, ncomp); amrex_atomic_accumulate_fab(BL_TO_FORTRAN_3D(filtered_rho), - BL_TO_FORTRAN_3D(rhofab), ncomp); + BL_TO_FORTRAN_N_3D(rhofab,icomp), ncomp); } else { amrex_atomic_accumulate_fab(BL_TO_FORTRAN_3D(local_rho), - BL_TO_FORTRAN_3D(rhofab), ncomp); + BL_TO_FORTRAN_N_3D(rhofab,icomp), ncomp); } }; - if (rho) depositCharge(rho); + if (rho) depositCharge(rho,0); if (! do_not_push) { @@ -1023,7 +1022,7 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_copy); } - if (rho2) depositCharge(rho2); + if (rho) depositCharge(rho,1); if (cost) { const Box& tbx = pti.tilebox(); diff --git a/Source/RigidInjectedParticleContainer.H b/Source/RigidInjectedParticleContainer.H index 93dda0aab..cce4dda99 100644 --- a/Source/RigidInjectedParticleContainer.H +++ b/Source/RigidInjectedParticleContainer.H @@ -29,7 +29,6 @@ public: amrex::MultiFab& jy, amrex::MultiFab& jz, amrex::MultiFab* rho, - amrex::MultiFab* rho2, amrex::Real t, amrex::Real dt) override; diff --git a/Source/RigidInjectedParticleContainer.cpp b/Source/RigidInjectedParticleContainer.cpp index eb9d6ee4f..3997eba59 100644 --- a/Source/RigidInjectedParticleContainer.cpp +++ b/Source/RigidInjectedParticleContainer.cpp @@ -303,8 +303,7 @@ RigidInjectedParticleContainer::Evolve (int lev, const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz, MultiFab& jx, MultiFab& jy, MultiFab& jz, - MultiFab* rho, MultiFab* rho2, - Real t, Real dt) + MultiFab* rho, Real t, Real dt) { // Update location of injection plane in the boosted frame @@ -323,8 +322,7 @@ RigidInjectedParticleContainer::Evolve (int lev, Ex, Ey, Ez, Bx, By, Bz, jx, jy, jz, - rho, rho2, - t, dt); + rho, t, dt); // Check if all done_injecting_temp are still true. done_injecting = std::all_of(done_injecting_temp.begin(), done_injecting_temp.end(), diff --git a/Source/WarpX.H b/Source/WarpX.H index f604e2bcd..1465219b6 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -411,22 +411,17 @@ private: bool is_synchronized = true; #ifdef WARPX_USE_PSATD - amrex::Vector< std::unique_ptr > rho2_fp; - amrex::Vector< std::unique_ptr > rho2_cp; - // Store fields in real space on the dual grid (i.e. the grid for the FFT push of the fields) // This includes data for the FFT guard cells (between FFT groups) amrex::Vector, 3 > > Efield_fp_fft; amrex::Vector, 3 > > Bfield_fp_fft; amrex::Vector, 3 > > current_fp_fft; - amrex::Vector< std::unique_ptr > rho_prev_fp_fft; - amrex::Vector< std::unique_ptr > rho_next_fp_fft; + amrex::Vector< std::unique_ptr > rho_fp_fft; amrex::Vector, 3 > > Efield_cp_fft; amrex::Vector, 3 > > Bfield_cp_fft; amrex::Vector, 3 > > current_cp_fft; - amrex::Vector< std::unique_ptr > rho_prev_cp_fft; - amrex::Vector< std::unique_ptr > rho_next_cp_fft; + amrex::Vector< std::unique_ptr > rho_cp_fft; public: // FFTData is a stuct containing a 22 pointers to auxiliary arrays diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index f9ebfc8a6..3aae523eb 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -151,19 +151,15 @@ WarpX::WarpX () costs.resize(nlevs_max); #ifdef WARPX_USE_PSATD - rho2_fp.resize(nlevs_max); - rho2_cp.resize(nlevs_max); - Efield_fp_fft.resize(nlevs_max); Bfield_fp_fft.resize(nlevs_max); current_fp_fft.resize(nlevs_max); - rho_prev_fp_fft.resize(nlevs_max); - rho_next_fp_fft.resize(nlevs_max); + rho_fp_fft.resize(nlevs_max); Efield_cp_fft.resize(nlevs_max); Bfield_cp_fft.resize(nlevs_max); current_cp_fft.resize(nlevs_max); - rho_prev_cp_fft.resize(nlevs_max); + rho_cp_fft.resize(nlevs_max); rho_next_cp_fft.resize(nlevs_max); dataptr_fp_fft.resize(nlevs_max); @@ -470,14 +466,8 @@ WarpX::ClearLevel (int lev) current_cp_fft[lev][i].reset(); } - rho2_fp[lev].reset(); - rho2_cp[lev].reset(); - - rho_prev_fp_fft[lev].reset(); - rho_next_fp_fft[lev].reset(); - - rho_prev_cp_fft[lev].reset(); - rho_next_cp_fft[lev].reset(); + rho_fp_fft[lev].reset(); + rho_cp_fft[lev].reset(); dataptr_fp_fft[lev].reset(); dataptr_cp_fft[lev].reset(); @@ -510,11 +500,12 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d #if (AMREX_SPACEDIM == 3) IntVect ngE(ngx,ngy,ngz); IntVect ngJ(ngx,ngy,ngz_nonci); - IntVect ngRho = ngJ; + IntVect ngRho = ngJ + 1; // One extra ghost cell, so that it's safe to deposit charge density + // after pushing particle. #elif (AMREX_SPACEDIM == 2) IntVect ngE(ngx,ngz); IntVect ngJ(ngx,ngz_nonci); - IntVect ngRho = ngJ; + IntVect ngRho = ngJ + 1; #endif int ngF = (do_moving_window) ? 2 : 0; @@ -539,6 +530,12 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d F_fp[lev].reset (new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,1, ngF)); rho_fp[lev].reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,1,ngRho)); } +#ifdef WARPX_USE_PSATD + else + { + rho_fp[lev].reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,1,ngRho)); + } +#endif // // The Aux patch (i.e., the full solution) @@ -589,6 +586,12 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d F_cp[lev].reset (new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,1, ngF)); rho_cp[lev].reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,1,ngRho)); } +#ifdef WARPX_USE_PSATD + else + { + rho_cp[lev].reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,1,ngRho)); + } +#endif } if (load_balance_int > 0) { diff --git a/Source/WarpXComm.cpp b/Source/WarpXComm.cpp index f71d3501c..501b4aedd 100644 --- a/Source/WarpXComm.cpp +++ b/Source/WarpXComm.cpp @@ -378,9 +378,10 @@ WarpX::SyncRho (const amrex::Vector >& rhof, for (int lev = 0; lev < finest_level; ++lev) { const auto& period = Geom(lev).periodicity(); + const int ncomp = rhoc[lev+1]->nComp(); const IntVect& ngsrc = rhoc[lev+1]->nGrowVect(); const IntVect ngdst = IntVect::TheZeroVector(); - rhof[lev]->copy(*rhoc[lev+1],0,0,1,ngsrc,ngdst,period,FabArrayBase::ADD); + rhof[lev]->copy(*rhoc[lev+1],0,0,ncomp,ngsrc,ngdst,period,FabArrayBase::ADD); } // Sum up coarse patch @@ -408,6 +409,7 @@ WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int ref_ratio) { BL_ASSERT(ref_ratio == 2); const IntVect& ng = fine.nGrowVect()/ref_ratio; + const int nc = fine.nComp(); #ifdef _OPEMP #pragma omp parallel @@ -418,13 +420,14 @@ WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int ref_ratio) { const Box& bx = mfi.growntilebox(ng); Box fbx = amrex::grow(amrex::refine(bx,ref_ratio),1); - ffab.resize(fbx); + ffab.resize(fbx, nc); fbx &= fine[mfi].box(); ffab.setVal(0.0); - ffab.copy(fine[mfi], fbx, 0, fbx, 0, 1); + ffab.copy(fine[mfi], fbx, 0, fbx, 0, nc); WRPX_SYNC_RHO(bx.loVect(), bx.hiVect(), - BL_TO_FORTRAN_ANYD(crse[mfi]), - BL_TO_FORTRAN_ANYD(ffab)); + BL_TO_FORTRAN_ANYD(crse[mfi]), + BL_TO_FORTRAN_ANYD(ffab), + &nc); } } } diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index b6d26d3a5..022a6e806 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) @@ -657,17 +654,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 diff --git a/Source/WarpXFFT.cpp b/Source/WarpXFFT.cpp index 0c1a5b619..04333bf20 100644 --- a/Source/WarpXFFT.cpp +++ b/Source/WarpXFFT.cpp @@ -129,17 +129,6 @@ WarpX::AllocLevelDataFFT (int lev) FFTDomainDecomposition(lev, ba_fp_fft, dm_fp_fft, ba_valid_fp_fft[lev], domain_fp_fft[lev], geom[lev].Domain()); - int ngRho = Efield_fp[lev][0]->nGrow(); - if (rho_fp[lev] == nullptr) - { - const BoxArray& ba = Efield_fp[lev][0]->boxArray(); - const DistributionMapping& dm = Efield_fp[lev][0]->DistributionMap(); - rho_fp[lev].reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,1,ngRho)); - } - - rho2_fp[lev].reset(new MultiFab(rho_fp[lev]->boxArray(), - rho_fp[lev]->DistributionMap(), - 1, ngRho+1)); // rho2 has one extra ghost cell, so that it's safe to deposit charge density after // pushing particle. @@ -161,10 +150,8 @@ WarpX::AllocLevelDataFFT (int lev) dm_fp_fft, 1, 0)); current_fp_fft[lev][2].reset(new MultiFab(amrex::convert(ba_fp_fft,jz_nodal_flag), dm_fp_fft, 1, 0)); - rho_prev_fp_fft[lev].reset(new MultiFab(amrex::convert(ba_fp_fft,IntVect::TheNodeVector()), - dm_fp_fft, 1, 0)); - rho_next_fp_fft[lev].reset(new MultiFab(amrex::convert(ba_fp_fft,IntVect::TheNodeVector()), - dm_fp_fft, 1, 0)); + rho_fp_fft[lev].reset(new MultiFab(amrex::convert(ba_fp_fft,IntVect::TheNodeVector()), + dm_fp_fft, 2, 0)); dataptr_fp_fft[lev].reset(new LayoutData(ba_fp_fft, dm_fp_fft)); @@ -175,20 +162,6 @@ WarpX::AllocLevelDataFFT (int lev) FFTDomainDecomposition(lev, ba_cp_fft, dm_cp_fft, ba_valid_cp_fft[lev], domain_cp_fft[lev], amrex::coarsen(geom[lev].Domain(),2)); - int ngRho = Efield_cp[lev][0]->nGrow(); - if (rho_cp[lev] == nullptr) - { - const BoxArray& ba = Efield_cp[lev][0]->boxArray(); - const DistributionMapping& dm = Efield_cp[lev][0]->DistributionMap(); - rho_cp[lev].reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,1,ngRho)); - } - - rho2_cp[lev].reset(new MultiFab(rho_cp[lev]->boxArray(), - rho_cp[lev]->DistributionMap(), - 1, ngRho)); - // rho2 has one extra ghost cell, so that it's safe to deposit charge density after - // pushing particle. - Efield_cp_fft[lev][0].reset(new MultiFab(amrex::convert(ba_cp_fft,Ex_nodal_flag), dm_cp_fft, 1, 0)); Efield_cp_fft[lev][1].reset(new MultiFab(amrex::convert(ba_cp_fft,Ey_nodal_flag), @@ -207,10 +180,8 @@ WarpX::AllocLevelDataFFT (int lev) dm_cp_fft, 1, 0)); current_cp_fft[lev][2].reset(new MultiFab(amrex::convert(ba_cp_fft,jz_nodal_flag), dm_cp_fft, 1, 0)); - rho_prev_cp_fft[lev].reset(new MultiFab(amrex::convert(ba_cp_fft,IntVect::TheNodeVector()), - dm_cp_fft, 1, 0)); - rho_next_cp_fft[lev].reset(new MultiFab(amrex::convert(ba_cp_fft,IntVect::TheNodeVector()), - dm_cp_fft, 1, 0)); + rho_cp_fft[lev].reset(new MultiFab(amrex::convert(ba_cp_fft,IntVect::TheNodeVector()), + dm_cp_fft, 2, 0)); dataptr_cp_fft[lev].reset(new LayoutData(ba_cp_fft, dm_cp_fft)); } @@ -420,8 +391,7 @@ WarpX::PushPSATD (int lev, amrex::Real /* dt */) current_fp_fft[lev][0]->ParallelCopy(*current_fp[lev][0], 0, 0, 1, 0, 0, period_fp); current_fp_fft[lev][1]->ParallelCopy(*current_fp[lev][1], 0, 0, 1, 0, 0, period_fp); current_fp_fft[lev][2]->ParallelCopy(*current_fp[lev][2], 0, 0, 1, 0, 0, period_fp); - rho_prev_fp_fft[lev]->ParallelCopy(*rho_fp[lev], 0, 0, 1, 0, 0, period_fp); - rho_next_fp_fft[lev]->ParallelCopy(*rho2_fp[lev], 0, 0, 1, 0, 0, period_fp); + rho_fp_fft[lev]->ParallelCopy(*rho_fp[lev], 0, 0, 1, 0, 0, period_fp); BL_PROFILE_VAR_STOP(blp_copy); BL_PROFILE_VAR_START(blp_push_eb); @@ -439,8 +409,8 @@ WarpX::PushPSATD (int lev, amrex::Real /* dt */) WARPX_TO_FORTRAN_ANYD((*current_fp_fft[lev][0])[mfi]), WARPX_TO_FORTRAN_ANYD((*current_fp_fft[lev][1])[mfi]), WARPX_TO_FORTRAN_ANYD((*current_fp_fft[lev][2])[mfi]), - WARPX_TO_FORTRAN_ANYD((*rho_prev_fp_fft[lev])[mfi]), - WARPX_TO_FORTRAN_ANYD((*rho_next_fp_fft[lev])[mfi])); + WARPX_TO_FORTRAN_N_ANYD((*rho_fp_fft[lev])[mfi],0), + WARPX_TO_FORTRAN_N_ANYD((*rho_fp_fft[lev])[mfi],1)); } } else if (Efield_fp_fft[lev][0]->local_size() == 0) diff --git a/Source/WarpXMove.cpp b/Source/WarpXMove.cpp index d1d40e2ff..7a4e346a3 100644 --- a/Source/WarpXMove.cpp +++ b/Source/WarpXMove.cpp @@ -165,7 +165,7 @@ WarpX::MoveWindow (bool move_j) } void -WarpX::shiftMF(MultiFab& mf, const Geometry& geom, int num_shift, int dir) +WarpX::shiftMF (MultiFab& mf, const Geometry& geom, int num_shift, int dir) { const BoxArray& ba = mf.boxArray(); const DistributionMapping& dm = mf.DistributionMap(); diff --git a/Source/WarpXParticleContainer.H b/Source/WarpXParticleContainer.H index f97c5cafe..8c904f5ba 100644 --- a/Source/WarpXParticleContainer.H +++ b/Source/WarpXParticleContainer.H @@ -94,8 +94,7 @@ public: const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz, amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, - amrex::MultiFab* rho, amrex::MultiFab* rho2, - amrex::Real t, amrex::Real dt) = 0; + amrex::MultiFab* rho, amrex::Real t, amrex::Real dt) = 0; virtual void PostRestart () = 0; diff --git a/Source/WarpXProb.cpp b/Source/WarpXProb.cpp index 3e5c248b8..ecbf1afaf 100644 --- a/Source/WarpXProb.cpp +++ b/Source/WarpXProb.cpp @@ -58,8 +58,7 @@ WarpX::InitLevelDataFFT (int lev, Real time) current_fp_fft[lev][0]->setVal(0.0); current_fp_fft[lev][1]->setVal(0.0); current_fp_fft[lev][2]->setVal(0.0); - rho_prev_fp_fft[lev]->setVal(0.0); - rho_next_fp_fft[lev]->setVal(0.0); + rho_fp_fft[lev]->setVal(0.0); if (lev > 0) { @@ -72,8 +71,7 @@ WarpX::InitLevelDataFFT (int lev, Real time) current_cp_fft[lev][0]->setVal(0.0); current_cp_fft[lev][1]->setVal(0.0); current_cp_fft[lev][2]->setVal(0.0); - rho_prev_cp_fft[lev]->setVal(0.0); - rho_next_cp_fft[lev]->setVal(0.0); + rho_cp_fft[lev]->setVal(0.0); } } diff --git a/Source/WarpXRegrid.cpp b/Source/WarpXRegrid.cpp index 489e471f4..c7aeb80bd 100644 --- a/Source/WarpXRegrid.cpp +++ b/Source/WarpXRegrid.cpp @@ -74,10 +74,11 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa } if (rho_fp[lev] != nullptr) { + const int nc = rho_fp[lev]->nComp(); const IntVect& ng = rho_fp[lev]->nGrowVect(); auto pmf = std::unique_ptr(new MultiFab(rho_fp[lev]->boxArray(), - dm, 1, ng)); - // pmf->Redistribute(*rho_fp[lev], 0, 0, 1, ng); + dm, nc, ng)); + // pmf->Redistribute(*rho_fp[lev], 0, 0, nc, ng); rho_fp[lev] = std::move(pmf); } @@ -145,10 +146,11 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa } if (rho_cp[lev] != nullptr) { + const int nc = rho_cp[lev]->nComp(); const IntVect& ng = rho_cp[lev]->nGrowVect(); auto pmf = std::unique_ptr(new MultiFab(rho_cp[lev]->boxArray(), - dm, 1, ng)); - // pmf->Redistribute(*rho_cp[lev], 0, 0, 1, ng); + dm, nc, ng)); + // pmf->Redistribute(*rho_cp[lev], 0, 0, nc, ng); rho_cp[lev] = std::move(pmf); } } diff --git a/Source/WarpX_f.F90 b/Source/WarpX_f.F90 index b1cbd26e1..46502a4d5 100644 --- a/Source/WarpX_f.F90 +++ b/Source/WarpX_f.F90 @@ -371,55 +371,59 @@ contains end subroutine warpx_clean_evec_3d - subroutine warpx_sync_rho_2d (lo, hi, crse, clo, chi, fine, flo, fhi) & + subroutine warpx_sync_rho_2d (lo, hi, crse, clo, chi, fine, flo, fhi, nc) & bind(c, name='warpx_sync_rho_2d') - integer, intent(in) :: lo(2), hi(2), flo(2), fhi(2), clo(2), chi(2) - real(amrex_real), intent(in ) :: fine(flo(1):fhi(1),flo(2):fhi(2)) - real(amrex_real), intent(inout) :: crse(clo(1):chi(1),clo(2):chi(2)) + integer, intent(in) :: lo(2), hi(2), flo(2), fhi(2), clo(2), chi(2), nc + real(amrex_real), intent(in ) :: fine(flo(1):fhi(1),flo(2):fhi(2),nc) + real(amrex_real), intent(inout) :: crse(clo(1):chi(1),clo(2):chi(2),nc) - integer :: i,j,ii,jj + integer :: i,j,ii,jj,m - do j = lo(2), hi(2) - jj = j*2 - do i = lo(1), hi(1) - ii = i*2 - crse(i,j) = 0.25d0 * & - ( fine(ii,jj) + 0.5d0 *(fine(ii-1,jj )+fine(ii+1,jj ) & - & + fine(ii ,jj-1)+fine(ii ,jj+1)) & - & + 0.25d0*(fine(ii-1,jj-1)+fine(ii+1,jj-1) & - & + fine(ii-1,jj+1)+fine(ii+1,jj+1)) ) + do m = 1, nc + do j = lo(2), hi(2) + jj = j*2 + do i = lo(1), hi(1) + ii = i*2 + crse(i,j,m) = 0.25d0 * & + ( fine(ii,jj,m) + 0.5d0 *(fine(ii-1,jj ,m)+fine(ii+1,jj ,m) & + & + fine(ii ,jj-1,m)+fine(ii ,jj+1,m)) & + & + 0.25d0*(fine(ii-1,jj-1,m)+fine(ii+1,jj-1,m) & + & + fine(ii-1,jj+1,m)+fine(ii+1,jj+1,m)) ) + end do end do end do end subroutine warpx_sync_rho_2d - subroutine warpx_sync_rho_3d (lo, hi, crse, clo, chi, fine, flo, fhi) & + subroutine warpx_sync_rho_3d (lo, hi, crse, clo, chi, fine, flo, fhi, nc) & bind(c, name='warpx_sync_rho_3d') - integer, intent(in) :: lo(3), hi(3), flo(3), fhi(3), clo(3), chi(3) - real(amrex_real), intent(in ) :: fine(flo(1):fhi(1),flo(2):fhi(2),flo(3):fhi(3)) - real(amrex_real), intent(inout) :: crse(clo(1):chi(1),clo(2):chi(2),clo(3):chi(3)) + integer, intent(in) :: lo(3), hi(3), flo(3), fhi(3), clo(3), chi(3), nc + real(amrex_real), intent(in ) :: fine(flo(1):fhi(1),flo(2):fhi(2),flo(3):fhi(3),nc) + real(amrex_real), intent(inout) :: crse(clo(1):chi(1),clo(2):chi(2),clo(3):chi(3),nc) - integer :: i,j,k,ii,jj,kk + integer :: i,j,k,ii,jj,kk,m - do k = lo(3), hi(3) - kk = k*2 - do j = lo(2), hi(2) - jj = j*2 - do i = lo(1), hi(1) - ii = i*2 - crse(i,j,k) = 0.125d0 * & - (fine(ii,jj,kk) + 0.5d0 *(fine(ii-1,jj ,kk )+fine(ii+1,jj ,kk ) & - & + fine(ii ,jj-1,kk )+fine(ii ,jj+1,kk ) & - & + fine(ii ,jj ,kk-1)+fine(ii ,jj ,kk+1)) & - & + 0.25d0 *(fine(ii-1,jj-1,kk )+fine(ii+1,jj-1,kk ) & - & + fine(ii-1,jj+1,kk )+fine(ii+1,jj+1,kk ) & - & + fine(ii-1,jj ,kk-1)+fine(ii+1,jj ,kk-1) & - & + fine(ii-1,jj ,kk+1)+fine(ii+1,jj ,kk+1) & - & + fine(ii ,jj-1,kk-1)+fine(ii ,jj+1,kk-1) & - & + fine(ii ,jj-1,kk+1)+fine(ii ,jj+1,kk+1)) & - & + 0.125d0*(fine(ii-1,jj-1,kk-1)+fine(ii-1,jj-1,kk+1) & - & + fine(ii-1,jj+1,kk-1)+fine(ii-1,jj+1,kk+1) & - & + fine(ii+1,jj-1,kk-1)+fine(ii+1,jj-1,kk+1) & - & + fine(ii+1,jj+1,kk-1)+fine(ii+1,jj+1,kk+1))) + do m = 1, nc + do k = lo(3), hi(3) + kk = k*2 + do j = lo(2), hi(2) + jj = j*2 + do i = lo(1), hi(1) + ii = i*2 + crse(i,j,k,m) = 0.125d0 * & + (fine(ii,jj,kk,m) + 0.5d0 *(fine(ii-1,jj ,kk ,m)+fine(ii+1,jj ,kk ,m) & + & + fine(ii ,jj-1,kk ,m)+fine(ii ,jj+1,kk ,m) & + & + fine(ii ,jj ,kk-1,m)+fine(ii ,jj ,kk+1,m)) & + & + 0.25d0 *(fine(ii-1,jj-1,kk ,m)+fine(ii+1,jj-1,kk ,m) & + & + fine(ii-1,jj+1,kk ,m)+fine(ii+1,jj+1,kk ,m) & + & + fine(ii-1,jj ,kk-1,m)+fine(ii+1,jj ,kk-1,m) & + & + fine(ii-1,jj ,kk+1,m)+fine(ii+1,jj ,kk+1,m) & + & + fine(ii ,jj-1,kk-1,m)+fine(ii ,jj+1,kk-1,m) & + & + fine(ii ,jj-1,kk+1,m)+fine(ii ,jj+1,kk+1,m)) & + & + 0.125d0*(fine(ii-1,jj-1,kk-1,m)+fine(ii-1,jj-1,kk+1,m) & + & + fine(ii-1,jj+1,kk-1,m)+fine(ii-1,jj+1,kk+1,m) & + & + fine(ii+1,jj-1,kk-1,m)+fine(ii+1,jj-1,kk+1,m) & + & + fine(ii+1,jj+1,kk-1,m)+fine(ii+1,jj+1,kk+1,m))) + end do end do end do end do diff --git a/Source/WarpX_f.H b/Source/WarpX_f.H index bdebffc68..5c6bb0b92 100644 --- a/Source/WarpX_f.H +++ b/Source/WarpX_f.H @@ -11,6 +11,7 @@ #define WARPX_TO_FORTRAN_BOX(x) WARPX_ARLIM_ANYD((x).loVect()), WARPX_ARLIM_ANYD((x).hiVect()) #define WARPX_TO_FORTRAN_ANYD(x) (x).dataPtr(), WARPX_ARLIM_ANYD((x).loVect()), WARPX_ARLIM_ANYD((x).hiVect()) +#define WARPX_TO_FORTRAN_N_ANYD(x,n) (x).dataPtr(n), WARPX_ARLIM_ANYD((x).loVect()), WARPX_ARLIM_ANYD((x).hiVect()) #endif @@ -415,8 +416,9 @@ extern "C" const int* dir); void WRPX_SYNC_RHO (const int* lo, const int* hi, - BL_FORT_FAB_ARG_ANYD(crse), - const BL_FORT_FAB_ARG_ANYD(fine)); + BL_FORT_FAB_ARG_ANYD(crse), + const BL_FORT_FAB_ARG_ANYD(fine), + const int* ncomp); void WRPX_CLEAN_EVEC (const int* xlo, const int* xhi, const int* ylo, const int* yhi, diff --git a/Source/main.cpp b/Source/main.cpp index 757c29e1f..d3926ef07 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -19,7 +19,7 @@ int main(int argc, char* argv[]) MPI_Init(&argc, &argv); #endif - amrex::Initialize(argc,argv,MPI_COMM_WORLD); + amrex::Initialize(argc,argv); BL_PROFILE_VAR("main()", pmain); -- cgit v1.2.3 From 429f550fc48f3022238f4989a0a2f22a66ae6bb4 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Wed, 8 Aug 2018 17:36:53 -0700 Subject: div E cleaning with PML --- Source/WarpX.H | 11 ++- Source/WarpX.cpp | 8 +- Source/WarpXEvolve.cpp | 41 ++++----- Source/WarpXFFT.cpp | 2 +- Source/WarpX_f.H | 37 ++++---- Source/WarpX_pml.F90 | 230 +++++++++++++++++++++++++------------------------ 6 files changed, 168 insertions(+), 161 deletions(-) (limited to 'Source/WarpXEvolve.cpp') diff --git a/Source/WarpX.H b/Source/WarpX.H index 1465219b6..b5b537209 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -26,6 +26,13 @@ #include #endif +enum struct DtType : int +{ + Full = 0, + FirstHalf, + SecondHalf +}; + class WarpX : public amrex::AmrCore { @@ -112,8 +119,8 @@ public: void EvolveE (int lev, amrex::Real dt); void EvolveB ( amrex::Real dt); void EvolveB (int lev, amrex::Real dt); - void EvolveF ( amrex::Real dt); - void EvolveF (int lev, amrex::Real dt); + void EvolveF ( amrex::Real dt, DtType dt_type); + void EvolveF (int lev, amrex::Real dt, DtType dt_type); void DampPML (); void DampPML (int lev); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 3aae523eb..d2c31b420 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -528,12 +528,12 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d if (do_dive_cleaning) { F_fp[lev].reset (new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,1, ngF)); - rho_fp[lev].reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,1,ngRho)); + rho_fp[lev].reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,2,ngRho)); } #ifdef WARPX_USE_PSATD else { - rho_fp[lev].reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,1,ngRho)); + rho_fp[lev].reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,2,ngRho)); } #endif @@ -584,12 +584,12 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d if (do_dive_cleaning) { F_cp[lev].reset (new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,1, ngF)); - rho_cp[lev].reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,1,ngRho)); + rho_cp[lev].reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2,ngRho)); } #ifdef WARPX_USE_PSATD else { - rho_cp[lev].reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,1,ngRho)); + rho_cp[lev].reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2,ngRho)); } #endif } diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 022a6e806..b4e1c8e95 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -110,11 +110,12 @@ WarpX::EvolveEM (int numsteps) FillBoundaryE(); FillBoundaryB(); #else - EvolveF(dt[0]); + EvolveF(0.5*dt[0], DtType::FirstHalf); EvolveB(0.5*dt[0]); // We now have B^{n+1/2} FillBoundaryB(); EvolveE(dt[0]); // We now have E^{n+1} FillBoundaryE(); + EvolveF(0.5*dt[0], DtType::SecondHalf); EvolveB(0.5*dt[0]); // We now have B^{n+1} if (do_pml) { DampPML(); @@ -481,8 +482,6 @@ WarpX::EvolveE (int lev, Real dt) if (pml_F) { - amrex::Abort("pml_F"); -#if 0 WRPX_PUSH_PML_EVEC_F( tex.loVect(), tex.hiVect(), tey.loVect(), tey.hiVect(), @@ -491,9 +490,7 @@ WarpX::EvolveE (int lev, Real dt) 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); -#endif + &dtsdx_c2[0]); } } } @@ -501,17 +498,17 @@ WarpX::EvolveE (int lev, Real dt) } void -WarpX::EvolveF (Real dt) +WarpX::EvolveF (Real dt, DtType dt_type) { if (!do_dive_cleaning) return; for (int lev = 0; lev <= finest_level; ++lev) { - EvolveF(lev, dt); + EvolveF(lev, dt, dt_type); } } void -WarpX::EvolveF (int lev, Real dt) +WarpX::EvolveF (int lev, Real dt, DtType dt_type) { if (!do_dive_cleaning) return; @@ -526,6 +523,7 @@ WarpX::EvolveF (int lev, Real dt) { int patch_level = (ipatch == 0) ? lev : lev-1; const auto& dx = WarpX::CellSize(patch_level); + const std::array dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]}; MultiFab *Ex, *Ey, *Ez, *rho, *F; if (ipatch == 0) @@ -545,22 +543,17 @@ WarpX::EvolveF (int lev, Real dt) 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); -// xxxxx -#if 0 if (do_pml && pml[lev]->ok()) { - const int dttype = static_cast(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 @@ -573,12 +566,9 @@ WarpX::EvolveF (int lev, Real dt) 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]); } } -#endif - } } @@ -607,7 +597,7 @@ WarpX::DampPML (int lev) { 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(); - // xxxxx const auto& 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(); const auto& sigba = (ipatch==0) ? pml[lev]->GetMultiSigmaBox_fp() : pml[lev]->GetMultiSigmaBox_cp(); @@ -636,6 +626,13 @@ WarpX::DampPML (int lev) 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])); + } } } } diff --git a/Source/WarpXFFT.cpp b/Source/WarpXFFT.cpp index 04333bf20..55dacf7bb 100644 --- a/Source/WarpXFFT.cpp +++ b/Source/WarpXFFT.cpp @@ -391,7 +391,7 @@ WarpX::PushPSATD (int lev, amrex::Real /* dt */) current_fp_fft[lev][0]->ParallelCopy(*current_fp[lev][0], 0, 0, 1, 0, 0, period_fp); current_fp_fft[lev][1]->ParallelCopy(*current_fp[lev][1], 0, 0, 1, 0, 0, period_fp); current_fp_fft[lev][2]->ParallelCopy(*current_fp[lev][2], 0, 0, 1, 0, 0, period_fp); - rho_fp_fft[lev]->ParallelCopy(*rho_fp[lev], 0, 0, 1, 0, 0, period_fp); + rho_fp_fft[lev]->ParallelCopy(*rho_fp[lev], 0, 0, 2, 0, 0, period_fp); BL_PROFILE_VAR_STOP(blp_copy); BL_PROFILE_VAR_START(blp_push_eb); diff --git a/Source/WarpX_f.H b/Source/WarpX_f.H index 5c6bb0b92..4e73e9433 100644 --- a/Source/WarpX_f.H +++ b/Source/WarpX_f.H @@ -28,6 +28,7 @@ #define WRPX_PUSH_PML_EVEC_F warpx_push_pml_evec_f_3d #define WRPX_PUSH_PML_F warpx_push_pml_f_3d #define WRPX_DAMP_PML warpx_damp_pml_3d +#define WRPX_DAMP_PML_F warpx_damp_pml_f_3d #define WRPX_SUM_FINE_TO_CRSE_NODAL warpx_sum_fine_to_crse_nodal_3d #define WRPX_ZERO_OUT_BNDRY warpx_zero_out_bndry_3d @@ -59,6 +60,7 @@ #define WRPX_PUSH_PML_EVEC_F warpx_push_pml_evec_f_2d #define WRPX_PUSH_PML_F warpx_push_pml_f_2d #define WRPX_DAMP_PML warpx_damp_pml_2d +#define WRPX_DAMP_PML_F warpx_damp_pml_f_2d #define WRPX_SUM_FINE_TO_CRSE_NODAL warpx_sum_fine_to_crse_nodal_2d #define WRPX_ZERO_OUT_BNDRY warpx_zero_out_bndry_2d @@ -362,30 +364,16 @@ extern "C" BL_FORT_FAB_ARG_3D(ey), BL_FORT_FAB_ARG_3D(ez), const BL_FORT_FAB_ARG_3D(f), - const amrex::Real* sigx1, int sigx1_lo, int sigx1_hi, - const amrex::Real* sigx2, int sigx2_lo, int sigx2_hi, -#if (AMREX_SPACEDIM == 3) - const amrex::Real* sigy1, int sigy1_lo, int sigy1_hi, - const amrex::Real* sigy2, int sigy2_lo, int sigy2_hi, -#endif - const amrex::Real* sigz1, int sigz1_lo, int sigz1_hi, - const amrex::Real* sigz2, int sigz2_lo, int sigz2_hi, - const amrex::Real* c2); + const amrex::Real* dtdx); void WRPX_PUSH_PML_F(const int* lo, const int* hi, BL_FORT_FAB_ARG_3D(f), const BL_FORT_FAB_ARG_3D(ex), const BL_FORT_FAB_ARG_3D(ey), const BL_FORT_FAB_ARG_3D(ez), - const amrex::Real* sigx1, int sigx1_lo, int sigx1_hi, - const amrex::Real* sigx2, int sigx2_lo, int sigx2_hi, -#if (AMREX_SPACEDIM == 3) - const amrex::Real* sigy1, int sigy1_lo, int sigy1_hi, - const amrex::Real* sigy2, int sigy2_lo, int sigy2_hi, -#endif - const amrex::Real* sigz1, int sigz1_lo, int sigz1_hi, - const amrex::Real* sigz2, int sigz2_lo, int sigz2_hi, - const amrex::Real* c2inv); + const amrex::Real* dtsdx, + const amrex::Real* dtsdy, + const amrex::Real* dtsdz); void WRPX_DAMP_PML (const int* texlo, const int* texhi, const int* teylo, const int* teyhi, @@ -410,6 +398,19 @@ extern "C" #endif const amrex::Real* sigbz, int sigbz_lo, int sigbz_hi); + void WRPX_DAMP_PML_F (const int* tndlo, const int* tndhi, + amrex::Real* F, const int* flo, const int* fhi, + const amrex::Real* sigex, int sigex_lo, int sigex_hi, +#if (AMREX_SPACEDIM == 3) + const amrex::Real* sigey, int sigey_lo, int sigey_hi, +#endif + const amrex::Real* sigez, int sigez_lo, int sigez_hi, + const amrex::Real* sigbx, int sigbx_lo, int sigbx_hi, +#if (AMREX_SPACEDIM == 3) + const amrex::Real* sigby, int sigby_lo, int sigby_hi, +#endif + const amrex::Real* sigbz, int sigbz_lo, int sigbz_hi); + void WRPX_SYNC_CURRENT (const int* lo, const int* hi, BL_FORT_FAB_ARG_ANYD(crse), const BL_FORT_FAB_ARG_ANYD(fine), diff --git a/Source/WarpX_pml.F90 b/Source/WarpX_pml.F90 index 8d929aa6c..e6428f398 100644 --- a/Source/WarpX_pml.F90 +++ b/Source/WarpX_pml.F90 @@ -473,51 +473,30 @@ contains & Ex, Exlo, Exhi, & & Ey, Eylo, Eyhi, & & Ez, Ezlo, Ezhi, & - & sigx1, sigx1_lo, sigx1_hi, & - & sigx2, sigx2_lo, sigx2_hi, & - & sigy1, sigy1_lo, sigy1_hi, & - & sigy2, sigy2_lo, sigy2_hi, & - & sigz1, sigz1_lo, sigz1_hi, & - & sigz2, sigz2_lo, sigz2_hi, c2inv) & + & dtdx, dtdy, dtdz) & bind(c,name='warpx_push_pml_f_3d') integer, intent(in) :: lo(3), hi(3), Exlo(3), Exhi(3), Eylo(3), Eyhi(3), Ezlo(3), Ezhi(3), & flo(3), fhi(3) - integer, intent(in), value :: sigx1_lo, sigx1_hi, sigx2_lo, sigx2_hi - integer, intent(in), value :: sigy1_lo, sigy1_hi, sigy2_lo, sigy2_hi - integer, intent(in), value :: sigz1_lo, sigz1_hi, sigz2_lo, sigz2_hi real(amrex_real), intent(inout) :: f ( flo(1): fhi(1), flo(2): fhi(2), flo(3): fhi(3),3) real(amrex_real), intent(in ) :: Ex (Exlo(1):Exhi(1),Exlo(2):Exhi(2),Exlo(3):Exhi(3),3) real(amrex_real), intent(in ) :: Ey (Eylo(1):Eyhi(1),Eylo(2):Eyhi(2),Eylo(3):Eyhi(3),3) real(amrex_real), intent(in ) :: Ez (Ezlo(1):Ezhi(1),Ezlo(2):Ezhi(2),Ezlo(3):Ezhi(3),3) - real(amrex_real), intent(in) :: sigx1(sigx1_lo:sigx1_hi) - real(amrex_real), intent(in) :: sigx2(sigx2_lo:sigx2_hi) - real(amrex_real), intent(in) :: sigy1(sigy1_lo:sigy1_hi) - real(amrex_real), intent(in) :: sigy2(sigy2_lo:sigy2_hi) - real(amrex_real), intent(in) :: sigz1(sigz1_lo:sigz1_hi) - real(amrex_real), intent(in) :: sigz2(sigz2_lo:sigz2_hi) - real(amrex_real), intent(in) :: c2inv + real(amrex_real), intent(in) :: dtdx, dtdy, dtdz integer :: i, j, k do k = lo(3), hi(3) do j = lo(2), hi(2) do i = lo(1), hi(1) - f(i,j,k,1) = sigx1(i)*f(i,j,k,1) + sigx2(i)*c2inv* & - & ((Ex(i,j,k,1)-Ex(i-1,j,k,1)) & - & + (Ex(i,j,k,2)-Ex(i-1,j,k,2)) & - & + (Ex(i,j,k,3)-Ex(i-1,j,k,3))) - end do - do i = lo(1), hi(1) - f(i,j,k,2) = sigy1(j)*f(i,j,k,2) + sigy2(j)*c2inv* & - & ((Ey(i,j,k,1)-Ey(i,j-1,k,1)) & - & + (Ey(i,j,k,2)-Ey(i,j-1,k,2)) & - & + (Ey(i,j,k,3)-Ey(i,j-1,k,3))) - end do - do i = lo(1), hi(1) - f(i,j,k,3) = sigz1(k)*f(i,j,k,3) + sigz2(k)*c2inv* & - & ((Ez(i,j,k,1)-Ez(i,j,k-1,1)) & - & + (Ez(i,j,k,2)-Ez(i,j,k-1,2)) & - & + (Ez(i,j,k,3)-Ez(i,j,k-1,3))) + f(i,j,k,1) = f(i,j,k,1) + dtdx*((Ex(i,j,k,1)-Ex(i-1,j,k,1)) & + & + (Ex(i,j,k,2)-Ex(i-1,j,k,2)) & + & + (Ex(i,j,k,3)-Ex(i-1,j,k,3))) + f(i,j,k,2) = f(i,j,k,2) + dtdy*((Ey(i,j,k,1)-Ey(i,j-1,k,1)) & + & + (Ey(i,j,k,2)-Ey(i,j-1,k,2)) & + & + (Ey(i,j,k,3)-Ey(i,j-1,k,3))) + f(i,j,k,3) = f(i,j,k,3) + dtdz*((Ez(i,j,k,1)-Ez(i,j,k-1,1)) & + & + (Ez(i,j,k,2)-Ez(i,j,k-1,2)) & + & + (Ez(i,j,k,3)-Ez(i,j,k-1,3))) end do end do end do @@ -528,39 +507,26 @@ contains & Ex, Exlo, Exhi, & & Ey, Eylo, Eyhi, & & Ez, Ezlo, Ezhi, & - & sigx1, sigx1_lo, sigx1_hi, & - & sigx2, sigx2_lo, sigx2_hi, & - & sigz1, sigz1_lo, sigz1_hi, & - & sigz2, sigz2_lo, sigz2_hi, c2inv) & + & dtdx, dtdy, dtdz) & bind(c,name='warpx_push_pml_f_2d') integer, intent(in) :: lo(2), hi(2), Exlo(2), Exhi(2), Eylo(2), Eyhi(2), Ezlo(2), Ezhi(2), & flo(2), fhi(2) - integer, intent(in), value :: sigx1_lo, sigx1_hi, sigx2_lo, sigx2_hi - integer, intent(in), value :: sigz1_lo, sigz1_hi, sigz2_lo, sigz2_hi real(amrex_real), intent(inout) :: f ( flo(1): fhi(1), flo(2): fhi(2),3) real(amrex_real), intent(in ) :: Ex (Exlo(1):Exhi(1),Exlo(2):Exhi(2),3) real(amrex_real), intent(in ) :: Ey (Eylo(1):Eyhi(1),Eylo(2):Eyhi(2),3) real(amrex_real), intent(in ) :: Ez (Ezlo(1):Ezhi(1),Ezlo(2):Ezhi(2),3) - real(amrex_real), intent(in) :: sigx1(sigx1_lo:sigx1_hi) - real(amrex_real), intent(in) :: sigx2(sigx2_lo:sigx2_hi) - real(amrex_real), intent(in) :: sigz1(sigz1_lo:sigz1_hi) - real(amrex_real), intent(in) :: sigz2(sigz2_lo:sigz2_hi) - real(amrex_real), intent(in) :: c2inv + real(amrex_real), intent(in) :: dtdx, dtdy, dtdz integer :: i, k do k = lo(2), hi(2) do i = lo(1), hi(1) - f(i,k,1) = sigx1(i)*f(i,k,1) + sigx2(i)*c2inv* & - & ((Ex(i,k,1)-Ex(i-1,k,1)) & - & + (Ex(i,k,2)-Ex(i-1,k,2)) & - & + (Ex(i,k,3)-Ex(i-1,k,3))) - end do - do i = lo(1), hi(1) - f(i,k,3) = sigz1(k)*f(i,k,3) + sigz2(k)*c2inv* & - & ((Ez(i,k,1)-Ez(i,k-1,1)) & - & + (Ez(i,k,2)-Ez(i,k-1,2)) & - & + (Ez(i,k,3)-Ez(i,k-1,3))) + f(i,k,1) = f(i,k,1) + dtdx*((Ex(i,k,1)-Ex(i-1,k,1)) & + & + (Ex(i,k,2)-Ex(i-1,k,2)) & + & + (Ex(i,k,3)-Ex(i-1,k,3))) + f(i,k,3) = f(i,k,3) + dtdz*((Ez(i,k,1)-Ez(i,k-1,1)) & + & + (Ez(i,k,2)-Ez(i,k-1,2)) & + & + (Ez(i,k,3)-Ez(i,k-1,3))) end do end do end subroutine warpx_push_pml_f_2d @@ -571,38 +537,24 @@ contains & Ey, Eylo, Eyhi, & & Ez, Ezlo, Ezhi, & & f, flo, fhi, & - & sigx1, sigx1_lo, sigx1_hi, & - & sigx2, sigx2_lo, sigx2_hi, & - & sigy1, sigy1_lo, sigy1_hi, & - & sigy2, sigy2_lo, sigy2_hi, & - & sigz1, sigz1_lo, sigz1_hi, & - & sigz2, sigz2_lo, sigz2_hi, c2) & + & dtdx) & bind(c,name='warpx_push_pml_evec_f_3d') integer, intent(in) :: xlo(3), xhi(3), ylo(3), yhi(3), zlo(3), zhi(3), & Exlo(3), Exhi(3), Eylo(3), Eyhi(3), Ezlo(3), Ezhi(3), flo(3), fhi(3) - integer, intent(in), value :: sigx1_lo, sigx1_hi, sigx2_lo, sigx2_hi - integer, intent(in), value :: sigy1_lo, sigy1_hi, sigy2_lo, sigy2_hi - integer, intent(in), value :: sigz1_lo, sigz1_hi, sigz2_lo, sigz2_hi real(amrex_real), intent(inout) :: Ex (Exlo(1):Exhi(1),Exlo(2):Exhi(2),Exlo(3):Exhi(3),3) real(amrex_real), intent(inout) :: Ey (Eylo(1):Eyhi(1),Eylo(2):Eyhi(2),Eylo(3):Eyhi(3),3) real(amrex_real), intent(inout) :: Ez (Ezlo(1):Ezhi(1),Ezlo(2):Ezhi(2),Ezlo(3):Ezhi(3),3) real(amrex_real), intent(in ) :: f ( flo(1): fhi(1), flo(2): fhi(2), flo(3): fhi(3),3) - real(amrex_real), intent(in) :: sigx1(sigx1_lo:sigx1_hi) - real(amrex_real), intent(in) :: sigx2(sigx2_lo:sigx2_hi) - real(amrex_real), intent(in) :: sigy1(sigy1_lo:sigy1_hi) - real(amrex_real), intent(in) :: sigy2(sigy2_lo:sigy2_hi) - real(amrex_real), intent(in) :: sigz1(sigz1_lo:sigz1_hi) - real(amrex_real), intent(in) :: sigz2(sigz2_lo:sigz2_hi) - real(amrex_real), intent(in) :: c2 + real(amrex_real), intent(in) :: dtdx(3) integer :: i, j, k do k = xlo(3), xhi(3) do j = xlo(2), xhi(2) do i = xlo(1), xhi(1) - Ex(i,j,k,3) = sigx1(i)*Ex(i,j,k,3) + sigx2(i)*c2*((f(i+1,j,k,1)-f(i,j,k,1)) & - & + (f(i+1,j,k,2)-f(i,j,k,2)) & - & + (f(i+1,j,k,3)-f(i,j,k,3))) + Ex(i,j,k,3) = Ex(i,j,k,3) + dtdx(1)*((f(i+1,j,k,1)-f(i,j,k,1)) & + & + (f(i+1,j,k,2)-f(i,j,k,2)) & + & + (f(i+1,j,k,3)-f(i,j,k,3))) end do end do end do @@ -610,9 +562,9 @@ contains do k = ylo(3), yhi(3) do j = ylo(2), yhi(2) do i = ylo(1), yhi(1) - Ey(i,j,k,3) = sigy1(j)*Ey(i,j,k,3) + sigy2(j)*c2*((f(i,j+1,k,1)-f(i,j,k,1)) & - & + (f(i,j+1,k,2)-f(i,j,k,2)) & - & + (f(i,j+1,k,3)-f(i,j,k,3))) + Ey(i,j,k,3) = Ey(i,j,k,3) + dtdx(2)*((f(i,j+1,k,1)-f(i,j,k,1)) & + & + (f(i,j+1,k,2)-f(i,j,k,2)) & + & + (f(i,j+1,k,3)-f(i,j,k,3))) end do end do end do @@ -620,9 +572,9 @@ contains do k = zlo(3), zhi(3) do j = zlo(2), zhi(2) do i = zlo(1), zhi(1) - Ez(i,j,k,3) = sigz1(k)*Ez(i,j,k,3) + sigz2(k)*c2*((f(i,j,k+1,1)-f(i,j,k,1)) & - & + (f(i,j,k+1,2)-f(i,j,k,2)) & - & + (f(i,j,k+1,3)-f(i,j,k,3))) + Ez(i,j,k,3) = Ez(i,j,k,3) + dtdx(3)*((f(i,j,k+1,1)-f(i,j,k,1)) & + & + (f(i,j,k+1,2)-f(i,j,k,2)) & + & + (f(i,j,k+1,3)-f(i,j,k,3))) end do end do end do @@ -635,40 +587,31 @@ contains & Ey, Eylo, Eyhi, & & Ez, Ezlo, Ezhi, & & f, flo, fhi, & - & sigx1, sigx1_lo, sigx1_hi, & - & sigx2, sigx2_lo, sigx2_hi, & - & sigz1, sigz1_lo, sigz1_hi, & - & sigz2, sigz2_lo, sigz2_hi, c2) & + & dtdx) & bind(c,name='warpx_push_pml_evec_f_2d') integer, intent(in) :: xlo(2), xhi(2), ylo(2), yhi(2), zlo(2), zhi(2), & Exlo(2), Exhi(2), Eylo(2), Eyhi(2), Ezlo(2), Ezhi(2), flo(2), fhi(2) - integer, intent(in), value :: sigx1_lo, sigx1_hi, sigx2_lo, sigx2_hi - integer, intent(in), value :: sigz1_lo, sigz1_hi, sigz2_lo, sigz2_hi real(amrex_real), intent(inout) :: Ex (Exlo(1):Exhi(1),Exlo(2):Exhi(2),3) real(amrex_real), intent(inout) :: Ey (Eylo(1):Eyhi(1),Eylo(2):Eyhi(2),3) real(amrex_real), intent(inout) :: Ez (Ezlo(1):Ezhi(1),Ezlo(2):Ezhi(2),3) real(amrex_real), intent(in ) :: f ( flo(1): fhi(1), flo(2): fhi(2),3) - real(amrex_real), intent(in) :: sigx1(sigx1_lo:sigx1_hi) - real(amrex_real), intent(in) :: sigx2(sigx2_lo:sigx2_hi) - real(amrex_real), intent(in) :: sigz1(sigz1_lo:sigz1_hi) - real(amrex_real), intent(in) :: sigz2(sigz2_lo:sigz2_hi) - real(amrex_real), intent(in) :: c2 + real(amrex_real), intent(in) :: dtdx(3) integer :: i, k do k = xlo(2), xhi(2) do i = xlo(1), xhi(1) - Ex(i,k,3) = sigx1(i)*Ex(i,k,3) + sigx2(i)*c2*((f(i+1,k,1)-f(i,k,1)) & - & + (f(i+1,k,2)-f(i,k,2)) & - & + (f(i+1,k,3)-f(i,k,3))) + Ex(i,k,3) = Ex(i,k,3) + dtdx(1)*((f(i+1,k,1)-f(i,k,1)) & + & + (f(i+1,k,2)-f(i,k,2)) & + & + (f(i+1,k,3)-f(i,k,3))) end do end do do k = zlo(2), zhi(2) do i = zlo(1), zhi(1) - Ez(i,k,3) = sigz1(k)*Ez(i,k,3) + sigz2(k)*c2*((f(i,k+1,1)-f(i,k,1)) & - & + (f(i,k+1,2)-f(i,k,2)) & - & + (f(i,k+1,3)-f(i,k,3))) + Ez(i,k,3) = Ez(i,k,3) + dtdx(3)*((f(i,k+1,1)-f(i,k,1)) & + & + (f(i,k+1,2)-f(i,k,2)) & + & + (f(i,k+1,3)-f(i,k,3))) end do end do @@ -680,13 +623,13 @@ contains & ex, exlo, exhi, ey, eylo, eyhi, ez, ezlo, ezhi, & & bx, bxlo, bxhi, by, bylo, byhi, bz, bzlo, bzhi, & & sigex, sexlo, sexhi, sigez, sezlo, sezhi, & - & sigbx, sbxlo, sbxhi, sigbz, sbzlo, sbzhi) & + & sigcx, scxlo, scxhi, sigcz, sczlo, sczhi) & bind(c,name='warpx_damp_pml_2d') integer, dimension(2), intent(in) :: texlo, texhi, teylo, teyhi, tezlo, tezhi, & tbxlo, tbxhi, tbylo, tbyhi, tbzlo, tbzhi, & exlo, exhi, eylo, eyhi, ezlo, ezhi, bxlo, bxhi, bylo, byhi, bzlo, bzhi integer, intent(in), value :: sexlo, sexhi, sezlo, sezhi, & - & sbxlo, sbxhi, sbzlo, sbzhi + & scxlo, scxhi, sczlo, sczhi real(amrex_real), intent(inout) :: ex(exlo(1):exhi(1),exlo(2):exhi(2),3) real(amrex_real), intent(inout) :: ey(eylo(1):eyhi(1),eylo(2):eyhi(2),3) real(amrex_real), intent(inout) :: ez(ezlo(1):ezhi(1),ezlo(2):ezhi(2),3) @@ -695,14 +638,15 @@ contains real(amrex_real), intent(inout) :: bz(bzlo(1):bzhi(1),bzlo(2):bzhi(2),2) real(amrex_real), intent(in) :: sigex(sexlo:sexhi) real(amrex_real), intent(in) :: sigez(sezlo:sezhi) - real(amrex_real), intent(in) :: sigbx(sbxlo:sbxhi) - real(amrex_real), intent(in) :: sigbz(sbzlo:sbzhi) + real(amrex_real), intent(in) :: sigcx(scxlo:scxhi) + real(amrex_real), intent(in) :: sigcz(sczlo:sczhi) integer :: i,k do k = texlo(2), texhi(2) do i = texlo(1), texhi(1) ex(i,k,2) = ex(i,k,2) * sigez(k) + ex(i,k,3) = ex(i,k,3) * sigcx(i) end do end do @@ -716,27 +660,29 @@ contains do k = tezlo(2), tezhi(2) do i = tezlo(1), tezhi(1) ez(i,k,1) = ez(i,k,1) * sigex(i) + ez(i,k,3) = ez(i,k,3) * sigcz(k) end do end do do k = tbxlo(2), tbxhi(2) do i = tbxlo(1), tbxhi(1) - bx(i,k,2) = bx(i,k,2) * sigbz(k) + bx(i,k,2) = bx(i,k,2) * sigcz(k) end do end do do k = tbylo(2), tbyhi(2) do i = tbylo(1), tbyhi(1) - by(i,k,1) = by(i,k,1) * sigbz(k) - by(i,k,2) = by(i,k,2) * sigbx(i) + by(i,k,1) = by(i,k,1) * sigcz(k) + by(i,k,2) = by(i,k,2) * sigcx(i) end do end do do k = tbzlo(2), tbzhi(2) do i = tbzlo(1), tbzhi(1) - bz(i,k,1) = bz(i,k,1) * sigbx(i) + bz(i,k,1) = bz(i,k,1) * sigcx(i) end do end do + end subroutine warpx_damp_pml_2d @@ -745,13 +691,13 @@ contains & ex, exlo, exhi, ey, eylo, eyhi, ez, ezlo, ezhi, & & bx, bxlo, bxhi, by, bylo, byhi, bz, bzlo, bzhi, & & sigex, sexlo, sexhi, sigey, seylo, seyhi, sigez, sezlo, sezhi, & - & sigbx, sbxlo, sbxhi, sigby, sbylo, sbyhi, sigbz, sbzlo, sbzhi) & + & sigcx, scxlo, scxhi, sigcy, scylo, scyhi, sigcz, sczlo, sczhi) & bind(c,name='warpx_damp_pml_3d') integer, dimension(3), intent(in) :: texlo, texhi, teylo, teyhi, tezlo, tezhi, & tbxlo, tbxhi, tbylo, tbyhi, tbzlo, tbzhi, & exlo, exhi, eylo, eyhi, ezlo, ezhi, bxlo, bxhi, bylo, byhi, bzlo, bzhi integer, intent(in), value :: sexlo, sexhi, seylo, seyhi, sezlo, sezhi, & - & sbxlo, sbxhi, sbylo, sbyhi, sbzlo, sbzhi + & scxlo, scxhi, scylo, scyhi, sczlo, sczhi real(amrex_real), intent(inout) :: ex(exlo(1):exhi(1),exlo(2):exhi(2),exlo(3):exhi(3),3) real(amrex_real), intent(inout) :: ey(eylo(1):eyhi(1),eylo(2):eyhi(2),eylo(3):eyhi(3),3) real(amrex_real), intent(inout) :: ez(ezlo(1):ezhi(1),ezlo(2):ezhi(2),ezlo(3):ezhi(3),3) @@ -761,9 +707,9 @@ contains real(amrex_real), intent(in) :: sigex(sexlo:sexhi) real(amrex_real), intent(in) :: sigey(seylo:seyhi) real(amrex_real), intent(in) :: sigez(sezlo:sezhi) - real(amrex_real), intent(in) :: sigbx(sbxlo:sbxhi) - real(amrex_real), intent(in) :: sigby(sbylo:sbyhi) - real(amrex_real), intent(in) :: sigbz(sbzlo:sbzhi) + real(amrex_real), intent(in) :: sigcx(scxlo:scxhi) + real(amrex_real), intent(in) :: sigcy(scylo:scyhi) + real(amrex_real), intent(in) :: sigcz(sczlo:sczhi) integer :: i,j,k @@ -772,6 +718,7 @@ contains do i = texlo(1), texhi(1) ex(i,j,k,1) = ex(i,j,k,1) * sigey(j) ex(i,j,k,2) = ex(i,j,k,2) * sigez(k) + ex(i,j,k,3) = ex(i,j,k,3) * sigcx(i) end do end do end do @@ -781,6 +728,7 @@ contains do i = teylo(1), teyhi(1) ey(i,j,k,1) = ey(i,j,k,1) * sigez(k) ey(i,j,k,2) = ey(i,j,k,2) * sigex(i) + ey(i,j,k,3) = ey(i,j,k,3) * sigcy(j) end do end do end do @@ -790,6 +738,7 @@ contains do i = tezlo(1), tezhi(1) ez(i,j,k,1) = ez(i,j,k,1) * sigex(i) ez(i,j,k,2) = ez(i,j,k,2) * sigey(j) + ez(i,j,k,3) = ez(i,j,k,3) * sigcz(k) end do end do end do @@ -797,8 +746,8 @@ contains do k = tbxlo(3), tbxhi(3) do j = tbxlo(2), tbxhi(2) do i = tbxlo(1), tbxhi(1) - bx(i,j,k,1) = bx(i,j,k,1) * sigby(j) - bx(i,j,k,2) = bx(i,j,k,2) * sigbz(k) + bx(i,j,k,1) = bx(i,j,k,1) * sigcy(j) + bx(i,j,k,2) = bx(i,j,k,2) * sigcz(k) end do end do end do @@ -806,8 +755,8 @@ contains do k = tbylo(3), tbyhi(3) do j = tbylo(2), tbyhi(2) do i = tbylo(1), tbyhi(1) - by(i,j,k,1) = by(i,j,k,1) * sigbz(k) - by(i,j,k,2) = by(i,j,k,2) * sigbx(i) + by(i,j,k,1) = by(i,j,k,1) * sigcz(k) + by(i,j,k,2) = by(i,j,k,2) * sigcx(i) end do end do end do @@ -815,11 +764,64 @@ contains do k = tbzlo(3), tbzhi(3) do j = tbzlo(2), tbzhi(2) do i = tbzlo(1), tbzhi(1) - bz(i,j,k,1) = bz(i,j,k,1) * sigbx(i) - bz(i,j,k,2) = bz(i,j,k,2) * sigby(j) + bz(i,j,k,1) = bz(i,j,k,1) * sigcx(i) + bz(i,j,k,2) = bz(i,j,k,2) * sigcy(j) end do end do end do + end subroutine warpx_damp_pml_3d + subroutine warpx_damp_pml_f_2d (tndlo, tndhi, f, flo, fhi,& + & sigex, sexlo, sexhi, sigez, sezlo, sezhi, & + & sigcx, scxlo, scxhi, sigcz, sczlo, sczhi) & + bind(c,name='warpx_damp_pml_f_2d') + integer, dimension(2), intent(in) :: tndlo, tndhi, flo, fhi + integer, intent(in), value :: sexlo, sexhi, sezlo, sezhi, & + & scxlo, scxhi, sczlo, sczhi + real(amrex_real), intent(inout) :: f ( flo(1): fhi(1), flo(2): fhi(2),3) + real(amrex_real), intent(in) :: sigex(sexlo:sexhi) + real(amrex_real), intent(in) :: sigez(sezlo:sezhi) + real(amrex_real), intent(in) :: sigcx(scxlo:scxhi) + real(amrex_real), intent(in) :: sigcz(sczlo:sczhi) + + integer :: i,k + + do k = tndlo(2), tndhi(2) + do i = tndlo(1), tndhi(1) + f(i,k,1) = f(i,k,1) * sigex(i) + f(i,k,3) = f(i,k,3) * sigez(k) + end do + end do + end subroutine warpx_damp_pml_f_2d + + + subroutine warpx_damp_pml_f_3d (tndlo, tndhi, f, flo, fhi,& + & sigex, sexlo, sexhi, sigey, seylo, seyhi, sigez, sezlo, sezhi, & + & sigcx, scxlo, scxhi, sigcy, scylo, scyhi, sigcz, sczlo, sczhi) & + bind(c,name='warpx_damp_pml_f_3d') + integer, dimension(3), intent(in) :: tndlo, tndhi, flo, fhi + integer, intent(in), value :: sexlo, sexhi, seylo, seyhi, sezlo, sezhi, & + & scxlo, scxhi, scylo, scyhi, sczlo, sczhi + real(amrex_real), intent(inout) :: f ( flo(1): fhi(1), flo(2): fhi(2), flo(3): fhi(3),3) + real(amrex_real), intent(in) :: sigex(sexlo:sexhi) + real(amrex_real), intent(in) :: sigey(seylo:seyhi) + real(amrex_real), intent(in) :: sigez(sezlo:sezhi) + real(amrex_real), intent(in) :: sigcx(scxlo:scxhi) + real(amrex_real), intent(in) :: sigcy(scylo:scyhi) + real(amrex_real), intent(in) :: sigcz(sczlo:sczhi) + + integer :: i,j,k + + do k = tndlo(3), tndhi(3) + do j = tndlo(2), tndhi(2) + do i = tndlo(1), tndhi(1) + f(i,j,k,1) = f(i,j,k,1) * sigex(i) + f(i,j,k,2) = f(i,j,k,2) * sigey(j) + f(i,j,k,3) = f(i,j,k,3) * sigez(k) + end do + end do + end do + end subroutine warpx_damp_pml_f_3d + end module warpx_pml_module -- cgit v1.2.3