diff options
Diffstat (limited to 'Source/WarpX.cpp')
| -rw-r--r-- | Source/WarpX.cpp | 169 |
1 files changed, 67 insertions, 102 deletions
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index c9084cd2a..d49ead1bc 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -24,12 +24,10 @@ using namespace amrex; -Vector<Real> WarpX::B_external_particle(3, 0.0); -Vector<Real> WarpX::E_external_particle(3, 0.0); - Vector<Real> WarpX::E_external_grid(3, 0.0); Vector<Real> WarpX::B_external_grid(3, 0.0); +std::string WarpX::authors = ""; std::string WarpX::B_ext_grid_s = "default"; std::string WarpX::E_ext_grid_s = "default"; @@ -57,6 +55,7 @@ long WarpX::charge_deposition_algo; long WarpX::field_gathering_algo; long WarpX::particle_pusher_algo; int WarpX::maxwell_fdtd_solver_id; +int WarpX::do_dive_cleaning = 0; long WarpX::n_rz_azimuthal_modes = 1; long WarpX::ncomps = 1; @@ -90,6 +89,7 @@ Real WarpX::particle_slice_width_lab = 0.0; bool WarpX::do_dynamic_scheduling = true; int WarpX::do_subcycling = 0; +bool WarpX::exchange_all_guard_cells = 0; #if (AMREX_SPACEDIM == 3) IntVect WarpX::Bx_nodal_flag(1,0,0); @@ -153,7 +153,7 @@ WarpX::WarpX () ReadParameters(); #ifdef WARPX_USE_OPENPMD - m_OpenPMDPlotWriter = new WarpXOpenPMDPlot(openpmd_tspf, openpmd_backend); + m_OpenPMDPlotWriter = new WarpXOpenPMDPlot(openpmd_tspf, openpmd_backend, WarpX::getPMLdirections()); #endif // Geometry on all levels has been defined already. @@ -263,6 +263,25 @@ WarpX::WarpX () // at different levels (the stencil depends on c*dt/dz) nci_godfrey_filter_exeybz.resize(nlevs_max); nci_godfrey_filter_bxbyez.resize(nlevs_max); + + // Sanity checks. Must be done after calling the MultiParticleContainer + // constructor, as it reads additional parameters + // (e.g., use_fdtd_nci_corr) + +#ifndef WARPX_USE_PSATD + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + not ( do_pml && do_nodal ), + "PML + do_nodal for finite-difference not implemented" + ); +#endif + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + not ( do_dive_cleaning && do_nodal ), + "divE cleaning + do_nodal not implemented" + ); +#ifdef WARPX_USE_PSATD + AMREX_ALWAYS_ASSERT(use_fdtd_nci_corr == 0); + AMREX_ALWAYS_ASSERT(do_subcycling == 0); +#endif } WarpX::~WarpX () @@ -288,6 +307,7 @@ WarpX::ReadParameters () ParmParse pp;// Traditionally, max_step and stop_time do not have prefix. pp.query("max_step", max_step); pp.query("stop_time", stop_time); + pp.query("authors", authors); } { @@ -309,6 +329,7 @@ WarpX::ReadParameters () pp.query("verbose", verbose); pp.query("regrid_int", regrid_int); pp.query("do_subcycling", do_subcycling); + pp.query("exchange_all_guard_cells", exchange_all_guard_cells); pp.query("override_sync_int", override_sync_int); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_subcycling != 1 || max_level <= 1, @@ -322,9 +343,6 @@ WarpX::ReadParameters () pp.query("zmax_plasma_to_compute_max_step", zmax_plasma_to_compute_max_step); - pp.queryarr("B_external_particle", B_external_particle); - pp.queryarr("E_external_particle", E_external_particle); - pp.query("do_moving_window", do_moving_window); if (do_moving_window) { @@ -646,7 +664,6 @@ WarpX::ReadParameters () } } - } // This is a virtual function. @@ -730,58 +747,23 @@ WarpX::ClearLevel (int lev) void WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& dm) { - // When using subcycling, the particles on the fine level perform two pushes - // before being redistributed ; therefore, we need one extra guard cell - // (the particles may move by 2*c*dt) - const int ngx_tmp = (maxLevel() > 0 && do_subcycling == 1) ? WarpX::nox+1 : WarpX::nox; - const int ngy_tmp = (maxLevel() > 0 && do_subcycling == 1) ? WarpX::noy+1 : WarpX::noy; - const int ngz_tmp = (maxLevel() > 0 && do_subcycling == 1) ? WarpX::noz+1 : WarpX::noz; - - // Ex, Ey, Ez, Bx, By, and Bz have the same number of ghost cells. - // jx, jy, jz and rho have the same number of ghost cells. - // E and B have the same number of ghost cells as j and rho if NCI filter is not used, - // but different number of ghost cells in z-direction if NCI filter is used. - // The number of cells should be even, in order to easily perform the - // interpolation from coarse grid to fine grid. - int ngx = (ngx_tmp % 2) ? ngx_tmp+1 : ngx_tmp; // Always even number - int ngy = (ngy_tmp % 2) ? ngy_tmp+1 : ngy_tmp; // Always even number - int ngz_nonci = (ngz_tmp % 2) ? ngz_tmp+1 : ngz_tmp; // Always even number - int ngz; - if (WarpX::use_fdtd_nci_corr) { - int ng = ngz_tmp + NCIGodfreyFilter::m_stencil_width; - ngz = (ng % 2) ? ng+1 : ng; - } else { - ngz = ngz_nonci; - } - // J is only interpolated from fine to coarse (not coarse to fine) - // and therefore does not need to be even. - int ngJx = ngx_tmp; - int ngJy = ngy_tmp; - int ngJz = ngz_tmp; - - // When calling the moving window (with one level of refinement), we shift - // the fine grid by 2 cells ; therefore, we need at least 2 guard cells - // on level 1. This may not be necessary for level 0. - if (do_moving_window) { - ngx = std::max(ngx,2); - ngy = std::max(ngy,2); - ngz = std::max(ngz,2); - ngJx = std::max(ngJx,2); - ngJy = std::max(ngJy,2); - ngJz = std::max(ngJz,2); - } - -#if (AMREX_SPACEDIM == 3) - IntVect ngE(ngx,ngy,ngz); - IntVect ngJ(ngJx,ngJy,ngJz); -#elif (AMREX_SPACEDIM == 2) - IntVect ngE(ngx,ngz); - IntVect ngJ(ngJx,ngJz); -#endif + bool aux_is_nodal = (field_gathering_algo == GatheringAlgo::MomentumConserving); - IntVect ngRho = ngJ+1; //One extra ghost cell, so that it's safe to deposit charge density - // after pushing particle. + guard_cells.Init( + do_subcycling, + WarpX::use_fdtd_nci_corr, + do_nodal, + do_moving_window, + fft_hybrid_mpi_decomposition, + aux_is_nodal, + moving_window_dir, + WarpX::nox, + nox_fft, noy_fft, noz_fft, + NCIGodfreyFilter::m_stencil_width, + maxwell_fdtd_solver_id, + maxLevel(), + exchange_all_guard_cells); if (mypc->nSpeciesDepositOnMainGrid() && n_current_deposition_buffer == 0) { n_current_deposition_buffer = 1; @@ -792,54 +774,22 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d } if (n_current_deposition_buffer < 0) { - n_current_deposition_buffer = ngJ.max(); + n_current_deposition_buffer = guard_cells.ng_alloc_J.max(); } if (n_field_gather_buffer < 0) { // Field gather buffer should be larger than current deposition buffers n_field_gather_buffer = n_current_deposition_buffer + 1; } - int ngF = (do_moving_window) ? 2 : 0; - // CKC solver requires one additional guard cell - if (maxwell_fdtd_solver_id == 1) ngF = std::max( ngF, 1 ); - -#ifdef WARPX_USE_PSATD - if (fft_hybrid_mpi_decomposition == false){ - // All boxes should have the same number of guard cells - // (to avoid temporary parallel copies) - // Thus take the max of the required number of guards for each field - // Also: the number of guard cell should be enough to contain - // the stencil of the FFT solver. Here, this number (`ngFFT`) - // is determined *empirically* to be the order of the solver - // for nodal, and half the order of the solver for staggered. - IntVect ngFFT; - if (do_nodal) { - ngFFT = IntVect(AMREX_D_DECL(nox_fft, noy_fft, noz_fft)); - } else { - ngFFT = IntVect(AMREX_D_DECL(nox_fft/2, noy_fft/2, noz_fft/2)); - } - for (int i_dim=0; i_dim<AMREX_SPACEDIM; i_dim++ ){ - int ng_required = ngFFT[i_dim]; - // Get the max - ng_required = std::max( ng_required, ngE[i_dim] ); - ng_required = std::max( ng_required, ngJ[i_dim] ); - ng_required = std::max( ng_required, ngRho[i_dim] ); - ng_required = std::max( ng_required, ngF ); - // Set the guard cells to this max - ngE[i_dim] = ng_required; - ngJ[i_dim] = ng_required; - ngRho[i_dim] = ng_required; - ngF = ng_required; - } - } -#endif - - AllocLevelMFs(lev, ba, dm, ngE, ngJ, ngRho, ngF); + AllocLevelMFs(lev, ba, dm, guard_cells.ng_alloc_EB, guard_cells.ng_alloc_J, + guard_cells.ng_alloc_Rho, guard_cells.ng_alloc_F, + guard_cells.ng_Extra, aux_is_nodal); } void WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm, - const IntVect& ngE, const IntVect& ngJ, const IntVect& ngRho, int ngF) + const IntVect& ngE, const IntVect& ngJ, const IntVect& ngRho, + const IntVect& ngF, const IntVect& ngextra, const bool aux_is_nodal) { #if defined WARPX_DIM_RZ @@ -851,9 +801,6 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm ncomps = n_rz_azimuthal_modes*2 - 1; #endif - bool aux_is_nodal = (field_gathering_algo == GatheringAlgo::MomentumConserving); - IntVect ngextra(static_cast<int>(aux_is_nodal and !do_nodal)); - // // The fine patch // @@ -883,7 +830,7 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm if (do_dive_cleaning) { - F_fp[lev].reset (new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,ncomps, ngF)); + F_fp[lev].reset (new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,ncomps, ngF.max())); } #ifdef WARPX_USE_PSATD else @@ -970,7 +917,7 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm } if (do_dive_cleaning) { - F_cp[lev].reset (new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,ncomps, ngF)); + F_cp[lev].reset (new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,ncomps, ngF.max())); } #ifdef WARPX_USE_PSATD else @@ -1280,6 +1227,24 @@ WarpX::GetPML (int lev) } } +std::vector< bool > +WarpX::getPMLdirections() const +{ + std::vector< bool > dirsWithPML( 6, false ); +#if AMREX_SPACEDIM!=3 + dirsWithPML.resize( 4 ); +#endif + if( do_pml ) + { + for( auto i = 0u; i < dirsWithPML.size() / 2u; ++i ) + { + dirsWithPML.at( 2u*i ) = bool(do_pml_Lo[i]); + dirsWithPML.at( 2u*i + 1u ) = bool(do_pml_Hi[i]); + } + } + return dirsWithPML; +} + void WarpX::BuildBufferMasks () { @@ -1362,7 +1327,7 @@ WarpX::Version () std::string WarpX::PicsarVersion () { -#ifdef WARPX_GIT_VERSION +#ifdef PICSAR_GIT_VERSION return std::string(PICSAR_GIT_VERSION); #else return std::string("Unknown"); |