diff options
| -rw-r--r-- | Source/WarpXEvolve.cpp | 152 | ||||
| -rw-r--r-- | Source/WarpXPML.cpp | 120 |
2 files changed, 136 insertions, 136 deletions
diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index d142e175b..46dd2f25a 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -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) { @@ -605,65 +605,65 @@ WarpX::EvolveF (int lev, PatchType patch_type, Real dt, DtType dt_type) void WarpX::DampPML () { - for (int lev = 0; lev <= finest_level; ++lev) { - DampPML(lev); - } + for (int lev = 0; lev <= finest_level; ++lev) { + DampPML(lev); + } } void WarpX::DampPML (int lev) { - DampPML(lev, PatchType::fine); - if (lev > 0) DampPML(lev, PatchType::coarse); + DampPML(lev, PatchType::fine); + if (lev > 0) DampPML(lev, PatchType::coarse); } void WarpX::DampPML (int lev, PatchType patch_type) { - if (!do_pml) return; + if (!do_pml) return; - BL_PROFILE("WarpX::DampPML()"); + BL_PROFILE("WarpX::DampPML()"); - if (pml[lev]->ok()) + if (pml[lev]->ok()) { - 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(); + 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 ) + 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])); + 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); - if (pml_F) { - const Box& tnd = mfi.nodaltilebox(); - WRPX_DAMP_PML_F(tnd.loVect(), tnd.hiVect(), - BL_TO_FORTRAN_3D((*pml_F)[mfi]), + 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])); + } } } } @@ -699,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); diff --git a/Source/WarpXPML.cpp b/Source/WarpXPML.cpp index f9b924725..ac8fda790 100644 --- a/Source/WarpXPML.cpp +++ b/Source/WarpXPML.cpp @@ -528,17 +528,17 @@ void PML::ExchangeB (PatchType patch_type, const std::array<amrex::MultiFab*,3>& Bp) { - if (patch_type == PatchType::fine && pml_B_fp[0] && Bp[0]) + if (patch_type == PatchType::fine && pml_B_fp[0] && Bp[0]) { - Exchange(*pml_B_fp[0], *Bp[0], *m_geom); - Exchange(*pml_B_fp[1], *Bp[1], *m_geom); - Exchange(*pml_B_fp[2], *Bp[2], *m_geom); + Exchange(*pml_B_fp[0], *Bp[0], *m_geom); + Exchange(*pml_B_fp[1], *Bp[1], *m_geom); + Exchange(*pml_B_fp[2], *Bp[2], *m_geom); } - else if (patch_type == PatchType::coarse && pml_B_cp[0] && Bp[0]) + else if (patch_type == PatchType::coarse && pml_B_cp[0] && Bp[0]) { - Exchange(*pml_B_cp[0], *Bp[0], *m_cgeom); - Exchange(*pml_B_cp[1], *Bp[1], *m_cgeom); - Exchange(*pml_B_cp[2], *Bp[2], *m_cgeom); + Exchange(*pml_B_cp[0], *Bp[0], *m_cgeom); + Exchange(*pml_B_cp[1], *Bp[1], *m_cgeom); + Exchange(*pml_B_cp[2], *Bp[2], *m_cgeom); } } @@ -546,43 +546,43 @@ void PML::ExchangeE (const std::array<amrex::MultiFab*,3>& E_fp, const std::array<amrex::MultiFab*,3>& E_cp) { - ExchangeB(PatchType::fine, E_fp); - ExchangeB(PatchType::coarse, E_cp); + ExchangeB(PatchType::fine, E_fp); + ExchangeB(PatchType::coarse, E_cp); } void PML::ExchangeE (PatchType patch_type, const std::array<amrex::MultiFab*,3>& Ep) { - if (patch_type == PatchType::fine && pml_E_fp[0] && Ep[0]) + if (patch_type == PatchType::fine && pml_E_fp[0] && Ep[0]) { - Exchange(*pml_E_fp[0], *Ep[0], *m_geom); - Exchange(*pml_E_fp[1], *Ep[1], *m_geom); - Exchange(*pml_E_fp[2], *Ep[2], *m_geom); + Exchange(*pml_E_fp[0], *Ep[0], *m_geom); + Exchange(*pml_E_fp[1], *Ep[1], *m_geom); + Exchange(*pml_E_fp[2], *Ep[2], *m_geom); } - else if (patch_type == PatchType::coarse && pml_E_cp[0] && Ep[0]) + else if (patch_type == PatchType::coarse && pml_E_cp[0] && Ep[0]) { - Exchange(*pml_E_cp[0], *Ep[0], *m_cgeom); - Exchange(*pml_E_cp[1], *Ep[1], *m_cgeom); - Exchange(*pml_E_cp[2], *Ep[2], *m_cgeom); + Exchange(*pml_E_cp[0], *Ep[0], *m_cgeom); + Exchange(*pml_E_cp[1], *Ep[1], *m_cgeom); + Exchange(*pml_E_cp[2], *Ep[2], *m_cgeom); } } void PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp) { - ExchangeF(PatchType::fine, F_fp); - ExchangeF(PatchType::coarse, F_cp); + ExchangeF(PatchType::fine, F_fp); + ExchangeF(PatchType::coarse, F_cp); } void PML::ExchangeF (PatchType patch_type, MultiFab* Fp) { - if (patch_type == PatchType::fine && pml_F_fp && Fp) { - Exchange(*pml_F_fp, *Fp, *m_geom); - } else if (patch_type == PatchType::coarse && pml_F_cp && Fp) { - Exchange(*pml_F_cp, *Fp, *m_cgeom); - } + if (patch_type == PatchType::fine && pml_F_fp && Fp) { + Exchange(*pml_F_fp, *Fp, *m_geom); + } else if (patch_type == PatchType::coarse && pml_F_cp && Fp) { + Exchange(*pml_F_cp, *Fp, *m_cgeom); + } } void @@ -639,74 +639,74 @@ PML::FillBoundary () void PML::FillBoundaryE () { - FillBoundaryE(PatchType::fine); - FillBoundaryE(PatchType::coarse); + FillBoundaryE(PatchType::fine); + FillBoundaryE(PatchType::coarse); } void PML::FillBoundaryE (PatchType patch_type) { - if (patch_type == PatchType::fine && pml_E_fp[0] && pml_E_fp[0]->nGrowVect().max() > 0) + if (patch_type == PatchType::fine && pml_E_fp[0] && pml_E_fp[0]->nGrowVect().max() > 0) { - const auto& period = m_geom->periodicity(); - pml_E_fp[0]->FillBoundary(period); - pml_E_fp[1]->FillBoundary(period); - pml_E_fp[2]->FillBoundary(period); - } - else if (patch_type == PatchType::coarse && pml_E_cp[0] && pml_E_cp[0]->nGrowVect().max() > 0) + const auto& period = m_geom->periodicity(); + pml_E_fp[0]->FillBoundary(period); + pml_E_fp[1]->FillBoundary(period); + pml_E_fp[2]->FillBoundary(period); + } + else if (patch_type == PatchType::coarse && pml_E_cp[0] && pml_E_cp[0]->nGrowVect().max() > 0) { - const auto& period = m_cgeom->periodicity(); - pml_E_cp[0]->FillBoundary(period); - pml_E_cp[1]->FillBoundary(period); - pml_E_cp[2]->FillBoundary(period); + const auto& period = m_cgeom->periodicity(); + pml_E_cp[0]->FillBoundary(period); + pml_E_cp[1]->FillBoundary(period); + pml_E_cp[2]->FillBoundary(period); } } void PML::FillBoundaryB () { - FillBoundaryB(PatchType::fine); - FillBoundaryB(PatchType::coarse); + FillBoundaryB(PatchType::fine); + FillBoundaryB(PatchType::coarse); } void PML::FillBoundaryB (PatchType patch_type) { - if (patch_type == PatchType::fine && pml_B_fp[0]) + if (patch_type == PatchType::fine && pml_B_fp[0]) { - const auto& period = m_geom->periodicity(); - pml_B_fp[0]->FillBoundary(period); - pml_B_fp[1]->FillBoundary(period); - pml_B_fp[2]->FillBoundary(period); - } - else if (patch_type == PatchType::coarse && pml_B_cp[0]) + const auto& period = m_geom->periodicity(); + pml_B_fp[0]->FillBoundary(period); + pml_B_fp[1]->FillBoundary(period); + pml_B_fp[2]->FillBoundary(period); + } + else if (patch_type == PatchType::coarse && pml_B_cp[0]) { - const auto& period = m_cgeom->periodicity(); - pml_B_cp[0]->FillBoundary(period); - pml_B_cp[1]->FillBoundary(period); - pml_B_cp[2]->FillBoundary(period); + const auto& period = m_cgeom->periodicity(); + pml_B_cp[0]->FillBoundary(period); + pml_B_cp[1]->FillBoundary(period); + pml_B_cp[2]->FillBoundary(period); } } void PML::FillBoundaryF () { - FillBoundaryF(PatchType::fine); - FillBoundaryF(PatchType::coarse); + FillBoundaryF(PatchType::fine); + FillBoundaryF(PatchType::coarse); } void PML::FillBoundaryF (PatchType patch_type) { - if (patch_type == PatchType::fine && pml_F_fp && pml_F_fp->nGrowVect().max() > 0) + if (patch_type == PatchType::fine && pml_F_fp && pml_F_fp->nGrowVect().max() > 0) { - const auto& period = m_geom->periodicity(); - pml_F_fp->FillBoundary(period); - } - else if (patch_type == PatchType::coarse && pml_F_cp && pml_F_cp->nGrowVect().max() > 0) + const auto& period = m_geom->periodicity(); + pml_F_fp->FillBoundary(period); + } + else if (patch_type == PatchType::coarse && pml_F_cp && pml_F_cp->nGrowVect().max() > 0) { - const auto& period = m_cgeom->periodicity(); - pml_F_cp->FillBoundary(period); + const auto& period = m_cgeom->periodicity(); + pml_F_cp->FillBoundary(period); } } |