diff options
Diffstat (limited to 'Source/WarpXEvolve.cpp')
| -rw-r--r-- | Source/WarpXEvolve.cpp | 792 |
1 files changed, 464 insertions, 328 deletions
diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 54c3db78f..a0fe641a3 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -60,6 +60,7 @@ WarpX::EvolveEM (int numsteps) #ifdef WARPX_USE_PSATD amrex::Abort("LoadBalance for PSATD: TODO"); #endif + if (step > 0 && (step+1) % load_balance_int == 0) { LoadBalance(); @@ -85,58 +86,27 @@ WarpX::EvolveEM (int numsteps) UpdateAuxilaryData(); // on first step, push p by -0.5*dt for (int lev = 0; lev <= finest_level; ++lev) { - mypc->PushP(lev, -0.5*dt[0], + mypc->PushP(lev, -0.5*dt[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]); } - is_synchronized = false; + is_synchronized = false; } else { // Beyond one step, we have E^{n} and B^{n}. // Particles have p^{n-1/2} and x^{n}. - FillBoundaryE(); - FillBoundaryB(); - UpdateAuxilaryData(); - } - - // Push particle from x^{n} to x^{n+1} - // from p^{n-1/2} to p^{n+1/2} - // Deposit current j^{n+1/2} - // Deposit charge density rho^{n} -#ifdef WARPX_USE_PY - if (warpx_py_particleinjection) warpx_py_particleinjection(); - if (warpx_py_particlescraper) warpx_py_particlescraper(); - if (warpx_py_beforedeposition) warpx_py_beforedeposition(); -#endif - PushParticlesandDepose(cur_time); -#ifdef WARPX_USE_PY - if (warpx_py_afterdeposition) warpx_py_afterdeposition(); -#endif - - SyncCurrent(); - - SyncRho(rho_fp, rho_cp); - - // Push E and B from {n} to {n+1} - // (And update guard cells immediately afterwards) -#ifdef WARPX_USE_PSATD - PushPSATD(dt[0]); - FillBoundaryE(); - FillBoundaryB(); -#else - EvolveF(0.5*dt[0], DtType::FirstHalf); - FillBoundaryF(); - 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(); FillBoundaryE(); + FillBoundaryB(); + UpdateAuxilaryData(); + } + + if (do_subcycling == 0 || finest_level == 0) { + OneStep_nosub(cur_time); + } else if (do_subcycling == 1 && finest_level == 1) { + OneStep_sub1(cur_time); + } else { + amrex::Print() << "Error: do_subcycling = " << do_subcycling << std::endl; + amrex::Abort("Unsupported do_subcycling type"); } - FillBoundaryB(); -#endif #ifdef WARPX_USE_PY if (warpx_py_beforeEsolve) warpx_py_beforeEsolve(); @@ -145,10 +115,10 @@ WarpX::EvolveEM (int numsteps) // At the end of last step, push p by 0.5*dt to synchronize UpdateAuxilaryData(); for (int lev = 0; lev <= finest_level; ++lev) { - mypc->PushP(lev, 0.5*dt[0], + mypc->PushP(lev, 0.5*dt[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]); - } + } is_synchronized = true; } #ifdef WARPX_USE_PY @@ -253,6 +223,173 @@ WarpX::EvolveEM (int numsteps) } void +WarpX::OneStep_nosub (Real cur_time) +{ + // Push particle from x^{n} to x^{n+1} + // from p^{n-1/2} to p^{n+1/2} + // Deposit current j^{n+1/2} + // Deposit charge density rho^{n} +#ifdef WARPX_USE_PY + if (warpx_py_particleinjection) warpx_py_particleinjection(); + if (warpx_py_particlescraper) warpx_py_particlescraper(); + if (warpx_py_beforedeposition) warpx_py_beforedeposition(); +#endif + PushParticlesandDepose(cur_time); +#ifdef WARPX_USE_PY + if (warpx_py_afterdeposition) warpx_py_afterdeposition(); +#endif + + SyncCurrent(); + + SyncRho(rho_fp, rho_cp); + + // Push E and B from {n} to {n+1} + // (And update guard cells immediately afterwards) +#ifdef WARPX_USE_PSATD + PushPSATD(dt[0]); + FillBoundaryE(); + FillBoundaryB(); +#else + EvolveF(0.5*dt[0], DtType::FirstHalf); + FillBoundaryF(); + 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(); + FillBoundaryE(); + } + FillBoundaryB(); +#endif +} + +void +WarpX::OneStep_sub1 (Real curtime) +{ + // TODO: we could save some charge depositions + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(finest_level == 1, "Must have exactly two levels"); + const int fine_lev = 1; + const int coarse_lev = 0; + + // i) + PushParticlesandDepose(fine_lev, curtime); + RestrictCurrentFromFineToCoarsePatch(fine_lev); + RestrictRhoFromFineToCoarsePatch(fine_lev); + ApplyFilterandSumBoundaryJ(fine_lev, PatchType::fine); + NodalSyncJ(fine_lev, PatchType::fine); + ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, 2); + NodalSyncRho(fine_lev, PatchType::fine, 0, 2); + + EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); + EvolveF(fine_lev, PatchType::fine, 0.5*dt[fine_lev], DtType::FirstHalf); + FillBoundaryB(fine_lev, PatchType::fine); + FillBoundaryF(fine_lev, PatchType::fine); + + EvolveE(fine_lev, PatchType::fine, dt[fine_lev]); + FillBoundaryE(fine_lev, PatchType::fine); + + EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); + EvolveF(fine_lev, PatchType::fine, 0.5*dt[fine_lev], DtType::SecondHalf); + + if (do_pml) { + DampPML(fine_lev, PatchType::fine); + FillBoundaryE(fine_lev, PatchType::fine); + } + + FillBoundaryB(fine_lev, PatchType::fine); + + // ii) + PushParticlesandDepose(coarse_lev, curtime); + StoreCurrent(coarse_lev); + AddCurrentFromFineLevelandSumBoundary(coarse_lev); + AddRhoFromFineLevelandSumBoundary(coarse_lev, 0, 1); + + EvolveB(fine_lev, PatchType::coarse, dt[fine_lev]); + EvolveF(fine_lev, PatchType::coarse, dt[fine_lev], DtType::FirstHalf); + FillBoundaryB(fine_lev, PatchType::coarse); + FillBoundaryF(fine_lev, PatchType::coarse); + + EvolveE(fine_lev, PatchType::coarse, dt[fine_lev]); + FillBoundaryE(fine_lev, PatchType::coarse); + + EvolveB(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev]); + EvolveF(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev], DtType::FirstHalf); + FillBoundaryB(coarse_lev, PatchType::fine); + FillBoundaryF(coarse_lev, PatchType::fine); + + EvolveE(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev]); + FillBoundaryE(coarse_lev, PatchType::fine); + + // iii) + UpdateAuxilaryData(); + + // iv) + PushParticlesandDepose(fine_lev, curtime+dt[fine_lev]); + RestrictCurrentFromFineToCoarsePatch(fine_lev); + RestrictRhoFromFineToCoarsePatch(fine_lev); + ApplyFilterandSumBoundaryJ(fine_lev, PatchType::fine); + NodalSyncJ(fine_lev, PatchType::fine); + ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, 2); + NodalSyncRho(fine_lev, PatchType::fine, 0, 2); + + EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); + EvolveF(fine_lev, PatchType::fine, 0.5*dt[fine_lev], DtType::FirstHalf); + FillBoundaryB(fine_lev, PatchType::fine); + FillBoundaryF(fine_lev, PatchType::fine); + + EvolveE(fine_lev, PatchType::fine, dt[fine_lev]); + FillBoundaryE(fine_lev, PatchType::fine); + + EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); + EvolveF(fine_lev, PatchType::fine, 0.5*dt[fine_lev], DtType::SecondHalf); + + if (do_pml) { + DampPML(fine_lev, PatchType::fine); + FillBoundaryE(fine_lev, PatchType::fine); + } + + FillBoundaryB(fine_lev, PatchType::fine); + FillBoundaryF(fine_lev, PatchType::fine); + + // v) + RestoreCurrent(coarse_lev); + AddCurrentFromFineLevelandSumBoundary(coarse_lev); + AddRhoFromFineLevelandSumBoundary(coarse_lev, 1, 1); + + EvolveE(fine_lev, PatchType::coarse, dt[fine_lev]); + FillBoundaryE(fine_lev, PatchType::coarse); + + EvolveB(fine_lev, PatchType::coarse, dt[fine_lev]); + EvolveF(fine_lev, PatchType::coarse, dt[fine_lev], DtType::SecondHalf); + + if (do_pml) { + DampPML(fine_lev, PatchType::coarse); // do it twice + DampPML(fine_lev, PatchType::coarse); + FillBoundaryE(fine_lev, PatchType::coarse); + } + + FillBoundaryB(fine_lev, PatchType::coarse); + FillBoundaryF(fine_lev, PatchType::coarse); + + EvolveE(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev]); + FillBoundaryE(coarse_lev, PatchType::fine); + + EvolveB(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev]); + EvolveF(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev], DtType::SecondHalf); + + if (do_pml) { + DampPML(coarse_lev, PatchType::fine); + FillBoundaryE(coarse_lev, PatchType::fine); + } + + FillBoundaryB(coarse_lev, PatchType::fine); +} + +void WarpX::EvolveB (Real dt) { for (int lev = 0; lev <= finest_level; ++lev) { @@ -264,107 +401,102 @@ void WarpX::EvolveB (int lev, Real dt) { BL_PROFILE("WarpX::EvolveB()"); + EvolveB(lev, PatchType::fine, dt); + if (lev > 0) { + EvolveB(lev, PatchType::coarse, dt); + } +} - // Parameters of the solver: order and mesh spacing - const int norder = 2; +void +WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) +{ + const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; + const std::array<Real,3>& dx = WarpX::CellSize(patch_level); + const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]}; + + MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz; + if (patch_type == PatchType::fine) + { + Ex = Efield_fp[lev][0].get(); + Ey = Efield_fp[lev][1].get(); + Ez = Efield_fp[lev][2].get(); + Bx = Bfield_fp[lev][0].get(); + By = Bfield_fp[lev][1].get(); + Bz = Bfield_fp[lev][2].get(); + } + else + { + Ex = Efield_cp[lev][0].get(); + Ey = Efield_cp[lev][1].get(); + Ez = Efield_cp[lev][2].get(); + Bx = Bfield_cp[lev][0].get(); + By = Bfield_cp[lev][1].get(); + Bz = Bfield_cp[lev][2].get(); + } - int npatches = (lev == 0) ? 1 : 2; + MultiFab* cost = costs[lev].get(); + const IntVect& rr = (lev > 0) ? refRatio(lev-1) : IntVect::TheUnitVector(); - for (int ipatch = 0; ipatch < npatches; ++ipatch) + // Loop through the grids, and over the tiles within each grid +#ifdef _OPENMP +#pragma omp parallel +#endif + for ( MFIter mfi(*Bx,true); mfi.isValid(); ++mfi ) { - int patch_level = (ipatch == 0) ? lev : lev-1; - const std::array<Real,3>& dx = WarpX::CellSize(patch_level); - const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]}; - - MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz; - if (ipatch == 0) - { - Ex = Efield_fp[lev][0].get(); - Ey = Efield_fp[lev][1].get(); - Ez = Efield_fp[lev][2].get(); - Bx = Bfield_fp[lev][0].get(); - By = Bfield_fp[lev][1].get(); - Bz = Bfield_fp[lev][2].get(); - } - else - { - Ex = Efield_cp[lev][0].get(); - Ey = Efield_cp[lev][1].get(); - Ez = Efield_cp[lev][2].get(); - Bx = Bfield_cp[lev][0].get(); - By = Bfield_cp[lev][1].get(); - Bz = Bfield_cp[lev][2].get(); + Real wt = amrex::second(); + + const Box& tbx = mfi.tilebox(Bx_nodal_flag); + const Box& tby = mfi.tilebox(By_nodal_flag); + const Box& tbz = mfi.tilebox(Bz_nodal_flag); + + // Call picsar routine for each tile + warpx_push_bvec( + tbx.loVect(), tbx.hiVect(), + tby.loVect(), tby.hiVect(), + tbz.loVect(), tbz.hiVect(), + BL_TO_FORTRAN_3D((*Ex)[mfi]), + BL_TO_FORTRAN_3D((*Ey)[mfi]), + BL_TO_FORTRAN_3D((*Ez)[mfi]), + BL_TO_FORTRAN_3D((*Bx)[mfi]), + BL_TO_FORTRAN_3D((*By)[mfi]), + BL_TO_FORTRAN_3D((*Bz)[mfi]), + &dtsdx[0], &dtsdx[1], &dtsdx[2], + &WarpX::maxwell_fdtd_solver_id); + + if (cost) { + Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)}); + if (patch_type == PatchType::coarse) cbx.refine(rr); + wt = (amrex::second() - wt) / cbx.d_numPts(); + (*cost)[mfi].plus(wt, cbx); } + } - MultiFab* cost = costs[lev].get(); - const IntVect& rr = (lev < finestLevel()) ? refRatio(lev) : IntVect::TheUnitVector(); + if (do_pml && pml[lev]->ok()) + { + const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); + const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); - // Loop through the grids, and over the tiles within each grid #ifdef _OPENMP #pragma omp parallel #endif - for ( MFIter mfi(*Bx,true); mfi.isValid(); ++mfi ) + for ( MFIter mfi(*pml_B[0],true); mfi.isValid(); ++mfi ) { - Real wt = amrex::second(); - const Box& tbx = mfi.tilebox(Bx_nodal_flag); const Box& tby = mfi.tilebox(By_nodal_flag); const Box& tbz = mfi.tilebox(Bz_nodal_flag); - - // Call picsar routine for each tile - warpx_push_bvec( + + WRPX_PUSH_PML_BVEC( tbx.loVect(), tbx.hiVect(), tby.loVect(), tby.hiVect(), tbz.loVect(), tbz.hiVect(), - BL_TO_FORTRAN_3D((*Ex)[mfi]), - BL_TO_FORTRAN_3D((*Ey)[mfi]), - BL_TO_FORTRAN_3D((*Ez)[mfi]), - BL_TO_FORTRAN_3D((*Bx)[mfi]), - BL_TO_FORTRAN_3D((*By)[mfi]), - BL_TO_FORTRAN_3D((*Bz)[mfi]), - &dtsdx[0], &dtsdx[1], &dtsdx[2], - &WarpX::maxwell_fdtd_solver_id); - - if (cost) { - Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)}); - if (ipatch == 1) cbx.refine(rr); - wt = (amrex::second() - wt) / cbx.d_numPts(); - (*cost)[mfi].plus(wt, cbx); - } - } - } - - if (do_pml && pml[lev]->ok()) - { - 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(); - int patch_level = (ipatch == 0) ? lev : lev-1; - const std::array<Real,3>& dx = WarpX::CellSize(patch_level); - const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]}; -#ifdef _OPENMP -#pragma omp parallel -#endif - for ( MFIter mfi(*pml_B[0],true); mfi.isValid(); ++mfi ) - { - 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_PUSH_PML_BVEC( - 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]), - &dtsdx[0], &dtsdx[1], &dtsdx[2], - &WarpX::maxwell_fdtd_solver_id); - } + 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]), + &dtsdx[0], &dtsdx[1], &dtsdx[2], + &WarpX::maxwell_fdtd_solver_id); } } } @@ -381,151 +513,145 @@ void WarpX::EvolveE (int lev, Real dt) { BL_PROFILE("WarpX::EvolveE()"); + EvolveE(lev, PatchType::fine, dt); + if (lev > 0) { + EvolveE(lev, PatchType::coarse, dt); + } +} +void +WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) +{ // Parameters of the solver: order and mesh spacing - 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 c2dt = (PhysConst::c*PhysConst::c) * dt; - int npatches = (lev == 0) ? 1 : 2; + int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; + const std::array<Real,3>& dx = WarpX::CellSize(patch_level); + const std::array<Real,3> dtsdx_c2 {c2dt/dx[0], c2dt/dx[1], c2dt/dx[2]}; - for (int ipatch = 0; ipatch < npatches; ++ipatch) + MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz, *jx, *jy, *jz, *F; + if (patch_type == PatchType::fine) { - int patch_level = (ipatch == 0) ? lev : lev-1; - const std::array<Real,3>& dx = WarpX::CellSize(patch_level); - const std::array<Real,3> dtsdx_c2 {c2dt/dx[0], c2dt/dx[1], c2dt/dx[2]}; - - MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz, *jx, *jy, *jz, *F; - if (ipatch == 0) - { - Ex = Efield_fp[lev][0].get(); - Ey = Efield_fp[lev][1].get(); - Ez = Efield_fp[lev][2].get(); - Bx = Bfield_fp[lev][0].get(); - By = Bfield_fp[lev][1].get(); - Bz = Bfield_fp[lev][2].get(); - jx = current_fp[lev][0].get(); - jy = current_fp[lev][1].get(); - jz = current_fp[lev][2].get(); - F = F_fp[lev].get(); - } - else - { - Ex = Efield_cp[lev][0].get(); - Ey = Efield_cp[lev][1].get(); - Ez = Efield_cp[lev][2].get(); - Bx = Bfield_cp[lev][0].get(); - By = Bfield_cp[lev][1].get(); - Bz = Bfield_cp[lev][2].get(); - jx = current_cp[lev][0].get(); - jy = current_cp[lev][1].get(); - jz = current_cp[lev][2].get(); - F = F_cp[lev].get(); - } + Ex = Efield_fp[lev][0].get(); + Ey = Efield_fp[lev][1].get(); + Ez = Efield_fp[lev][2].get(); + Bx = Bfield_fp[lev][0].get(); + By = Bfield_fp[lev][1].get(); + Bz = Bfield_fp[lev][2].get(); + jx = current_fp[lev][0].get(); + jy = current_fp[lev][1].get(); + jz = current_fp[lev][2].get(); + F = F_fp[lev].get(); + } + else if (patch_type == PatchType::coarse) + { + Ex = Efield_cp[lev][0].get(); + Ey = Efield_cp[lev][1].get(); + Ez = Efield_cp[lev][2].get(); + Bx = Bfield_cp[lev][0].get(); + By = Bfield_cp[lev][1].get(); + Bz = Bfield_cp[lev][2].get(); + jx = current_cp[lev][0].get(); + jy = current_cp[lev][1].get(); + jz = current_cp[lev][2].get(); + F = F_cp[lev].get(); + } - MultiFab* cost = costs[lev].get(); - const IntVect& rr = (lev < finestLevel()) ? refRatio(lev) : IntVect::TheUnitVector(); + MultiFab* cost = costs[lev].get(); + const IntVect& rr = (lev > 0) ? refRatio(lev-1) : IntVect::TheUnitVector(); - // Loop through the grids, and over the tiles within each grid + // Loop through the grids, and over the tiles within each grid #ifdef _OPENMP #pragma omp parallel #endif - for ( MFIter mfi(*Ex,true); mfi.isValid(); ++mfi ) - { - Real wt = amrex::second(); - - const Box& tex = mfi.tilebox(Ex_nodal_flag); - const Box& tey = mfi.tilebox(Ey_nodal_flag); - const Box& tez = mfi.tilebox(Ez_nodal_flag); - - // Call picsar routine for each tile - warpx_push_evec( + for ( MFIter mfi(*Ex,true); mfi.isValid(); ++mfi ) + { + Real wt = amrex::second(); + + const Box& tex = mfi.tilebox(Ex_nodal_flag); + const Box& tey = mfi.tilebox(Ey_nodal_flag); + const Box& tez = mfi.tilebox(Ez_nodal_flag); + + // Call picsar routine for each tile + warpx_push_evec( + tex.loVect(), tex.hiVect(), + tey.loVect(), tey.hiVect(), + tez.loVect(), tez.hiVect(), + BL_TO_FORTRAN_3D((*Ex)[mfi]), + BL_TO_FORTRAN_3D((*Ey)[mfi]), + BL_TO_FORTRAN_3D((*Ez)[mfi]), + BL_TO_FORTRAN_3D((*Bx)[mfi]), + BL_TO_FORTRAN_3D((*By)[mfi]), + BL_TO_FORTRAN_3D((*Bz)[mfi]), + BL_TO_FORTRAN_3D((*jx)[mfi]), + BL_TO_FORTRAN_3D((*jy)[mfi]), + BL_TO_FORTRAN_3D((*jz)[mfi]), + &mu_c2_dt, + &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]); + + if (F) { + warpx_push_evec_f( tex.loVect(), tex.hiVect(), tey.loVect(), tey.hiVect(), tez.loVect(), tez.hiVect(), BL_TO_FORTRAN_3D((*Ex)[mfi]), BL_TO_FORTRAN_3D((*Ey)[mfi]), BL_TO_FORTRAN_3D((*Ez)[mfi]), - BL_TO_FORTRAN_3D((*Bx)[mfi]), - BL_TO_FORTRAN_3D((*By)[mfi]), - BL_TO_FORTRAN_3D((*Bz)[mfi]), - BL_TO_FORTRAN_3D((*jx)[mfi]), - BL_TO_FORTRAN_3D((*jy)[mfi]), - BL_TO_FORTRAN_3D((*jz)[mfi]), - &mu_c2_dt, - &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]); - - if (F) { - - // Call picsar routine for each tile - warpx_push_evec_f( - tex.loVect(), tex.hiVect(), - tey.loVect(), tey.hiVect(), - tez.loVect(), tez.hiVect(), - BL_TO_FORTRAN_3D((*Ex)[mfi]), - BL_TO_FORTRAN_3D((*Ey)[mfi]), - BL_TO_FORTRAN_3D((*Ez)[mfi]), - BL_TO_FORTRAN_3D((*F)[mfi]), - &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2], - &WarpX::maxwell_fdtd_solver_id); - - } - - if (cost) { - Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)}); - if (ipatch == 1) cbx.refine(rr); - wt = (amrex::second() - wt) / cbx.d_numPts(); - (*cost)[mfi].plus(wt, cbx); - } + BL_TO_FORTRAN_3D((*F)[mfi]), + &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2], + &WarpX::maxwell_fdtd_solver_id); + } + + if (cost) { + Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)}); + if (patch_type == PatchType::coarse) cbx.refine(rr); + wt = (amrex::second() - wt) / cbx.d_numPts(); + (*cost)[mfi].plus(wt, cbx); } } if (do_pml && pml[lev]->ok()) { + if (F) pml[lev]->ExchangeF(patch_type, F); - 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& pml_F = (ipatch==0) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); - int patch_level = (ipatch == 0) ? lev : lev-1; - const std::array<Real,3>& dx = WarpX::CellSize(patch_level); - const std::array<Real,3> dtsdx_c2 {c2dt/dx[0], c2dt/dx[1], c2dt/dx[2]}; + const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); + const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); + const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); #ifdef _OPENMP #pragma omp parallel #endif - for ( MFIter mfi(*pml_E[0],true); mfi.isValid(); ++mfi ) + 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); + + WRPX_PUSH_PML_EVEC( + tex.loVect(), tex.hiVect(), + tey.loVect(), tey.hiVect(), + tez.loVect(), tez.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]), + &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]); + + if (pml_F) { - const Box& tex = mfi.tilebox(Ex_nodal_flag); - const Box& tey = mfi.tilebox(Ey_nodal_flag); - const Box& tez = mfi.tilebox(Ez_nodal_flag); - - WRPX_PUSH_PML_EVEC( + WRPX_PUSH_PML_EVEC_F( tex.loVect(), tex.hiVect(), tey.loVect(), tey.hiVect(), tez.loVect(), tez.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]), - &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]); - - if (pml_F) - { - WRPX_PUSH_PML_EVEC_F( - tex.loVect(), tex.hiVect(), - tey.loVect(), tey.hiVect(), - tez.loVect(), tez.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_F )[mfi]), - &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2], - &WarpX::maxwell_fdtd_solver_id); - } + BL_TO_FORTRAN_3D((*pml_F )[mfi]), + &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2], + &WarpX::maxwell_fdtd_solver_id); } } } @@ -546,62 +672,66 @@ WarpX::EvolveF (int lev, Real dt, DtType dt_type) { if (!do_dive_cleaning) return; + EvolveF(lev, PatchType::fine, dt, dt_type); + if (lev > 0) EvolveF(lev, PatchType::coarse, dt, dt_type); +} + +void +WarpX::EvolveF (int lev, PatchType patch_type, Real dt, DtType dt_type) +{ + if (!do_dive_cleaning) return; + BL_PROFILE("WarpX::EvolveF()"); static constexpr Real c2inv = 1.0/(PhysConst::c*PhysConst::c); static constexpr Real mu_c2 = PhysConst::mu0*PhysConst::c*PhysConst::c; - int npatches = (lev == 0) ? 1 : 2; + int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; + const auto& dx = WarpX::CellSize(patch_level); + const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]}; - for (int ipatch = 0; ipatch < npatches; ++ipatch) + MultiFab *Ex, *Ey, *Ez, *rho, *F; + if (patch_type == PatchType::fine) { - int patch_level = (ipatch == 0) ? lev : lev-1; - const auto& dx = WarpX::CellSize(patch_level); - const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]}; - - MultiFab *Ex, *Ey, *Ez, *rho, *F; - if (ipatch == 0) - { - Ex = Efield_fp[lev][0].get(); - Ey = Efield_fp[lev][1].get(); - Ez = Efield_fp[lev][2].get(); - rho = rho_fp[lev].get(); - F = F_fp[lev].get(); - } - else - { - Ex = Efield_cp[lev][0].get(); - Ey = Efield_cp[lev][1].get(); - Ez = Efield_cp[lev][2].get(); - rho = rho_cp[lev].get(); - 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, rhocomp, 0, 1, 0); - MultiFab::Saxpy(*F, dt, src, 0, 0, 1, 0); - - if (do_pml && pml[lev]->ok()) - { - 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(); + Ex = Efield_fp[lev][0].get(); + Ey = Efield_fp[lev][1].get(); + Ez = Efield_fp[lev][2].get(); + rho = rho_fp[lev].get(); + F = F_fp[lev].get(); + } + else + { + Ex = Efield_cp[lev][0].get(); + Ey = Efield_cp[lev][1].get(); + Ez = Efield_cp[lev][2].get(); + rho = rho_cp[lev].get(); + 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, rhocomp, 0, 1, 0); + MultiFab::Saxpy(*F, dt, src, 0, 0, 1, 0); + + if (do_pml && pml[lev]->ok()) + { + const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); + const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); #ifdef _OPENMP #pragma omp parallel #endif - for ( MFIter mfi(*pml_F,true); mfi.isValid(); ++mfi ) - { - const Box& bx = mfi.tilebox(); - WRPX_PUSH_PML_F(bx.loVect(), bx.hiVect(), - BL_TO_FORTRAN_ANYD((*pml_F )[mfi]), - BL_TO_FORTRAN_ANYD((*pml_E[0])[mfi]), - BL_TO_FORTRAN_ANYD((*pml_E[1])[mfi]), - BL_TO_FORTRAN_ANYD((*pml_E[2])[mfi]), - &dtsdx[0], &dtsdx[1], &dtsdx[2]); - } + for ( MFIter mfi(*pml_F,true); mfi.isValid(); ++mfi ) + { + const Box& bx = mfi.tilebox(); + WRPX_PUSH_PML_F(bx.loVect(), bx.hiVect(), + BL_TO_FORTRAN_ANYD((*pml_F )[mfi]), + BL_TO_FORTRAN_ANYD((*pml_E[0])[mfi]), + BL_TO_FORTRAN_ANYD((*pml_E[1])[mfi]), + BL_TO_FORTRAN_ANYD((*pml_E[2])[mfi]), + &dtsdx[0], &dtsdx[1], &dtsdx[2]); } } } @@ -609,8 +739,6 @@ WarpX::EvolveF (int lev, Real dt, DtType dt_type) void WarpX::DampPML () { - if (!do_pml) return; - for (int lev = 0; lev <= finest_level; ++lev) { DampPML(lev); } @@ -619,54 +747,56 @@ WarpX::DampPML () void WarpX::DampPML (int lev) { + DampPML(lev, PatchType::fine); + if (lev > 0) DampPML(lev, PatchType::coarse); +} + +void +WarpX::DampPML (int lev, PatchType patch_type) +{ 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()) { - if (pml[lev]->ok()) - { - const auto& pml_E = (ipatch==0) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); - const auto& pml_B = (ipatch==0) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); - const auto& pml_F = (ipatch==0) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); - const auto& sigba = (ipatch==0) ? pml[lev]->GetMultiSigmaBox_fp() - : pml[lev]->GetMultiSigmaBox_cp(); + const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); + const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); + const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); + const auto& sigba = (patch_type == PatchType::fine) ? pml[lev]->GetMultiSigmaBox_fp() + : pml[lev]->GetMultiSigmaBox_cp(); #ifdef _OPENMP #pragma omp parallel #endif - for ( MFIter mfi(*pml_E[0],true); mfi.isValid(); ++mfi ) - { - const Box& tex = mfi.tilebox(Ex_nodal_flag); - const Box& tey = mfi.tilebox(Ey_nodal_flag); - const Box& tez = mfi.tilebox(Ez_nodal_flag); - const Box& tbx = mfi.tilebox(Bx_nodal_flag); - const Box& tby = mfi.tilebox(By_nodal_flag); - const Box& tbz = mfi.tilebox(Bz_nodal_flag); - - WRPX_DAMP_PML(tex.loVect(), tex.hiVect(), - tey.loVect(), tey.hiVect(), - tez.loVect(), tez.hiVect(), - tbx.loVect(), tbx.hiVect(), - tby.loVect(), tby.hiVect(), - tbz.loVect(), tbz.hiVect(), - BL_TO_FORTRAN_3D((*pml_E[0])[mfi]), - BL_TO_FORTRAN_3D((*pml_E[1])[mfi]), - BL_TO_FORTRAN_3D((*pml_E[2])[mfi]), - BL_TO_FORTRAN_3D((*pml_B[0])[mfi]), - BL_TO_FORTRAN_3D((*pml_B[1])[mfi]), - BL_TO_FORTRAN_3D((*pml_B[2])[mfi]), - WRPX_PML_TO_FORTRAN(sigba[mfi])); - - if (pml_F) { - const Box& tnd = mfi.nodaltilebox(); - WRPX_DAMP_PML_F(tnd.loVect(), tnd.hiVect(), - BL_TO_FORTRAN_3D((*pml_F)[mfi]), - WRPX_PML_TO_FORTRAN(sigba[mfi])); - } + for ( MFIter mfi(*pml_E[0],true); mfi.isValid(); ++mfi ) + { + const Box& tex = mfi.tilebox(Ex_nodal_flag); + const Box& tey = mfi.tilebox(Ey_nodal_flag); + const Box& tez = mfi.tilebox(Ez_nodal_flag); + const Box& tbx = mfi.tilebox(Bx_nodal_flag); + const Box& tby = mfi.tilebox(By_nodal_flag); + const Box& tbz = mfi.tilebox(Bz_nodal_flag); + + WRPX_DAMP_PML(tex.loVect(), tex.hiVect(), + tey.loVect(), tey.hiVect(), + tez.loVect(), tez.hiVect(), + tbx.loVect(), tbx.hiVect(), + tby.loVect(), tby.hiVect(), + tbz.loVect(), tbz.hiVect(), + BL_TO_FORTRAN_3D((*pml_E[0])[mfi]), + BL_TO_FORTRAN_3D((*pml_E[1])[mfi]), + BL_TO_FORTRAN_3D((*pml_E[2])[mfi]), + BL_TO_FORTRAN_3D((*pml_B[0])[mfi]), + BL_TO_FORTRAN_3D((*pml_B[1])[mfi]), + BL_TO_FORTRAN_3D((*pml_B[2])[mfi]), + WRPX_PML_TO_FORTRAN(sigba[mfi])); + + if (pml_F) { + const Box& tnd = mfi.nodaltilebox(); + WRPX_DAMP_PML_F(tnd.loVect(), tnd.hiVect(), + BL_TO_FORTRAN_3D((*pml_F)[mfi]), + WRPX_PML_TO_FORTRAN(sigba[mfi])); } } } @@ -705,8 +835,8 @@ WarpX::ComputeDt () if (maxwell_fdtd_solver_id == 0) { // CFL time step Yee solver deltat = cfl * 1./( std::sqrt(AMREX_D_TERM( 1./(dx[0]*dx[0]), - + 1./(dx[1]*dx[1]), - + 1./(dx[2]*dx[2]))) * PhysConst::c ); + + 1./(dx[1]*dx[1]), + + 1./(dx[2]*dx[2]))) * PhysConst::c ); } else { // CFL time step CKC solver #if (BL_SPACEDIM == 3) @@ -719,6 +849,12 @@ WarpX::ComputeDt () dt.resize(0); dt.resize(max_level+1,deltat); + if (do_subcycling) { + for (int lev = max_level-1; lev >= 0; --lev) { + dt[lev] = dt[lev+1] * refRatio(lev)[0]; + } + } + if (do_electrostatic) { dt[0] = const_dt; } |