diff options
Diffstat (limited to 'Source/WarpXEvolve.cpp')
| -rw-r--r-- | Source/WarpXEvolve.cpp | 318 |
1 files changed, 161 insertions, 157 deletions
diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index f1a0faa15..d142e175b 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -263,7 +263,8 @@ WarpX::EvolveB (int lev, Real dt) { BL_PROFILE("WarpX::EvolveB()"); EvolveB(lev, PatchType::fine, dt); - if (lev > 0) { + if (lev > 0) + { EvolveB(lev, PatchType::coarse, dt); } } @@ -271,47 +272,47 @@ WarpX::EvolveB (int lev, Real dt) 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]}; + 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) + 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(); + 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 + 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(); + 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(); } - MultiFab* cost = costs[lev].get(); - const IntVect& rr = (lev > 0) ? refRatio(lev-1) : 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(*Bx,true); mfi.isValid(); ++mfi ) + for ( MFIter mfi(*Bx,true); mfi.isValid(); ++mfi ) { - Real wt = amrex::second(); + 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); + 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( + // Call picsar routine for each tile + warpx_push_bvec( tbx.loVect(), tbx.hiVect(), tby.loVect(), tby.hiVect(), tbz.loVect(), tbz.hiVect(), @@ -324,29 +325,29 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) &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); - } + 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 (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(); + 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(); #ifdef _OPENMP #pragma omp parallel #endif - for ( MFIter mfi(*pml_B[0],true); mfi.isValid(); ++mfi ) + 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( + 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(), @@ -365,78 +366,80 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) void WarpX::EvolveE (Real dt) { - for (int lev = 0; lev <= finest_level; ++lev) { - EvolveE(lev, dt); - } + for (int lev = 0; lev <= finest_level; ++lev) + { + EvolveE(lev, dt); + } } void WarpX::EvolveE (int lev, Real dt) { - BL_PROFILE("WarpX::EvolveE()"); - EvolveE(lev, PatchType::fine, dt); - if (lev > 0) { - EvolveE(lev, PatchType::coarse, 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 - 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; + // Parameters of the solver: order and mesh spacing + 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 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]}; + 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]}; - MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz, *jx, *jy, *jz, *F; - if (patch_type == PatchType::fine) + MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz, *jx, *jy, *jz, *F; + 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(); - jx = current_fp[lev][0].get(); - jy = current_fp[lev][1].get(); - jz = current_fp[lev][2].get(); - F = F_fp[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) + 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(); + 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 > 0) ? refRatio(lev-1) : 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 ) + 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( + 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(), @@ -451,9 +454,10 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) BL_TO_FORTRAN_3D((*jz)[mfi]), &mu_c2_dt, &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]); - - if (F) { - warpx_push_evec_f( + + if (F) + { + warpx_push_evec_f( tex.loVect(), tex.hiVect(), tey.loVect(), tey.hiVect(), tez.loVect(), tez.hiVect(), @@ -463,33 +467,33 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) 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 (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 (do_pml && pml[lev]->ok()) { - if (F) pml[lev]->ExchangeF(patch_type, F); + if (F) pml[lev]->ExchangeF(patch_type, F); - 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(); + 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( + 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(), @@ -501,9 +505,9 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) BL_TO_FORTRAN_3D((*pml_B[2])[mfi]), &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]); - if (pml_F) + if (pml_F) { - WRPX_PUSH_PML_EVEC_F( + WRPX_PUSH_PML_EVEC_F( tex.loVect(), tex.hiVect(), tey.loVect(), tey.hiVect(), tez.loVect(), tez.hiVect(), @@ -518,77 +522,77 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) } } - void WarpX::EvolveF (Real dt, DtType dt_type) { - if (!do_dive_cleaning) return; + if (!do_dive_cleaning) return; - for (int lev = 0; lev <= finest_level; ++lev) { - EvolveF(lev, dt, dt_type); - } + for (int lev = 0; lev <= finest_level; ++lev) + { + EvolveF(lev, dt, dt_type); + } } void WarpX::EvolveF (int lev, Real dt, DtType dt_type) { - if (!do_dive_cleaning) return; + if (!do_dive_cleaning) return; - EvolveF(lev, PatchType::fine, dt, dt_type); - if (lev > 0) EvolveF(lev, PatchType::coarse, dt, dt_type); + 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; + if (!do_dive_cleaning) return; - BL_PROFILE("WarpX::EvolveF()"); + 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; + static constexpr Real c2inv = 1.0/(PhysConst::c*PhysConst::c); + static constexpr Real mu_c2 = PhysConst::mu0*PhysConst::c*PhysConst::c; - 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]}; + 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]}; - MultiFab *Ex, *Ey, *Ez, *rho, *F; - if (patch_type == PatchType::fine) + MultiFab *Ex, *Ey, *Ez, *rho, *F; + if (patch_type == PatchType::fine) { - 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(); + 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 + 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(); + 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; + 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); + 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()) + 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(); + 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 ) + for ( MFIter mfi(*pml_F,true); mfi.isValid(); ++mfi ) { - const Box& bx = mfi.tilebox(); - WRPX_PUSH_PML_F(bx.loVect(), bx.hiVect(), + 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]), |