diff options
Diffstat (limited to 'Source/WarpXEvolve.cpp')
| -rw-r--r-- | Source/WarpXEvolve.cpp | 662 |
1 files changed, 330 insertions, 332 deletions
diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index be5152fc2..17037e576 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -49,7 +49,7 @@ WarpX::EvolveEM (int numsteps) { Real walltime_beg_step = amrex::second(); - // Start loop on time steps + // Start loop on time steps amrex::Print() << "\nSTEP " << step+1 << " starts ...\n"; #ifdef WARPX_USE_PY if (warpx_py_beforestep) warpx_py_beforestep(); @@ -63,16 +63,16 @@ WarpX::EvolveEM (int numsteps) if (step > 0 && (step+1) % load_balance_int == 0) { LoadBalance(); - // Reset the costs to 0 - for (int lev = 0; lev <= finest_level; ++lev) { - costs[lev]->setVal(0.0); - } + // Reset the costs to 0 + for (int lev = 0; lev <= finest_level; ++lev) { + costs[lev]->setVal(0.0); + } } for (int lev = 0; lev <= finest_level; ++lev) { - // Perform running average of the costs - // (Giving more importance to most recent costs) - (*costs[lev].get()).mult( (1. - 2./load_balance_int) ); + // Perform running average of the costs + // (Giving more importance to most recent costs) + (*costs[lev].get()).mult( (1. - 2./load_balance_int) ); } } @@ -89,12 +89,12 @@ WarpX::EvolveEM (int numsteps) *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}. UpdateAuxilaryData(); - } + } // Push particle from x^{n} to x^{n+1} // from p^{n-1/2} to p^{n+1/2} @@ -146,7 +146,7 @@ WarpX::EvolveEM (int numsteps) mypc->PushP(lev, 0.5*dt[0], *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 @@ -157,14 +157,14 @@ WarpX::EvolveEM (int numsteps) ++istep[lev]; } - cur_time += dt[0]; + cur_time += dt[0]; bool to_make_plot = (plot_int > 0) && ((step+1) % plot_int == 0); bool move_j = is_synchronized || to_make_plot; // If is_synchronized we need to shift j too so that next step we can evolve E by dt/2. // We might need to move j because we are going to make a plotfile. - MoveWindow(move_j); + MoveWindow(move_j); if (max_level == 0) { mypc->RedistributeLocal(); @@ -181,10 +181,10 @@ WarpX::EvolveEM (int numsteps) << " s; This step = " << walltime_end_step-walltime_beg_step << " s; Avg. per step = " << walltime/(step+1) << " s\n"; - // sync up time - for (int i = 0; i <= max_level; ++i) { - t_new[i] = cur_time; - } + // sync up time + for (int i = 0; i <= max_level; ++i) { + t_new[i] = cur_time; + } if (do_boosted_frame_diagnostic) { std::unique_ptr<MultiFab> cell_centered_data = nullptr; @@ -194,7 +194,7 @@ WarpX::EvolveEM (int numsteps) myBFD->writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); } - if (to_make_plot) + if (to_make_plot) { FillBoundaryE(); FillBoundaryB(); @@ -206,24 +206,24 @@ WarpX::EvolveEM (int numsteps) *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2]); } - last_plot_file_step = step+1; - WritePlotFile(); - } + last_plot_file_step = step+1; + WritePlotFile(); + } - if (check_int > 0 && (step+1) % check_int == 0) { - last_check_file_step = step+1; - WriteCheckPointFile(); - } + if (check_int > 0 && (step+1) % check_int == 0) { + last_check_file_step = step+1; + WriteCheckPointFile(); + } - if (cur_time >= stop_time - 1.e-3*dt[0]) { - max_time_reached = true; - break; - } + if (cur_time >= stop_time - 1.e-3*dt[0]) { + max_time_reached = true; + break; + } #ifdef WARPX_USE_PY if (warpx_py_afterstep) warpx_py_afterstep(); #endif - // End loop on time steps + // End loop on time steps } if (plot_int > 0 && istep[0] > last_plot_file_step && (max_time_reached || istep[0] >= max_step)) @@ -238,11 +238,11 @@ WarpX::EvolveEM (int numsteps) *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2]); } - WritePlotFile(); + WritePlotFile(); } if (check_int > 0 && istep[0] > last_check_file_step && (max_time_reached || istep[0] >= max_step)) { - WriteCheckPointFile(); + WriteCheckPointFile(); } if (do_boosted_frame_diagnostic) { @@ -262,107 +262,103 @@ void WarpX::EvolveB (int lev, Real dt) { BL_PROFILE("WarpX::EvolveB()"); - - // Parameters of the solver: order and mesh spacing - const int norder = 2; - - int npatches = (lev == 0) ? 1 : 2; - - for (int ipatch = 0; ipatch < npatches; ++ipatch) + EvolveB(lev, PatchType::fine, dt); + if (lev > 0) { - 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]}; + EvolveB(lev, PatchType::coarse, dt); + } +} - 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(); - } +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(); + } - 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(*Bx,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( - 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); - } + for ( MFIter mfi(*Bx,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( + 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); } } 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]}; + 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 ) - { - 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); - } + 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); } } } @@ -370,7 +366,8 @@ WarpX::EvolveB (int lev, Real dt) void WarpX::EvolveE (Real dt) { - for (int lev = 0; lev <= finest_level; ++lev) { + for (int lev = 0; lev <= finest_level; ++lev) + { EvolveE(lev, dt); } } @@ -379,151 +376,147 @@ 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 ) + 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) { - 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) { - - // 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); - - } + 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); - } + 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( - 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); - } + 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); } } } @@ -534,7 +527,8 @@ WarpX::EvolveF (Real dt, DtType dt_type) { if (!do_dive_cleaning) return; - for (int lev = 0; lev <= finest_level; ++lev) { + for (int lev = 0; lev <= finest_level; ++lev) + { EvolveF(lev, dt, dt_type); } } @@ -544,62 +538,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]); } } } @@ -607,8 +605,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); } @@ -617,54 +613,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])); } } } @@ -701,18 +699,18 @@ WarpX::ComputeDt () Real deltat = 0.; if (maxwell_fdtd_solver_id == 0) { - // CFL time step Yee solver - deltat = cfl * 1./( std::sqrt(AMREX_D_TERM( 1./(dx[0]*dx[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 ); } else { - // CFL time step CKC solver + // CFL time step CKC solver #if (BL_SPACEDIM == 3) - const Real delta = std::min(dx[0],std::min(dx[1],dx[2])); + const Real delta = std::min(dx[0],std::min(dx[1],dx[2])); #elif (BL_SPACEDIM == 2) - const Real delta = std::min(dx[0],dx[1]); + const Real delta = std::min(dx[0],dx[1]); #endif - deltat = cfl*delta/PhysConst::c; + deltat = cfl*delta/PhysConst::c; } dt.resize(0); dt.resize(max_level+1,deltat); |