diff options
Diffstat (limited to 'Source/WarpX.cpp')
| -rw-r--r-- | Source/WarpX.cpp | 464 |
1 files changed, 248 insertions, 216 deletions
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 07f124820..5a51d7d13 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -1,17 +1,3 @@ - -#include <limits> -#include <algorithm> -#include <cctype> -#include <cmath> -#include <numeric> - -#ifdef _OPENMP -#include <omp.h> -#endif - -#include <AMReX_ParmParse.H> -#include <AMReX_MultiFabUtil.H> - #include <WarpX.H> #include <WarpX_f.H> #include <WarpXConst.H> @@ -20,16 +6,29 @@ #include <WarpXAlgorithmSelection.H> #include <WarpX_FDTD.H> +#include <AMReX_ParmParse.H> +#include <AMReX_MultiFabUtil.H> #ifdef BL_USE_SENSEI_INSITU -#include <AMReX_AmrMeshInSituBridge.H> +# include <AMReX_AmrMeshInSituBridge.H> #endif +#ifdef _OPENMP +# include <omp.h> +#endif + +#include <limits> +#include <algorithm> +#include <cctype> +#include <cmath> +#include <numeric> + using namespace amrex; Vector<Real> WarpX::B_external(3, 0.0); int WarpX::do_moving_window = 0; int WarpX::moving_window_dir = -1; +Real WarpX::moving_window_v = std::numeric_limits<amrex::Real>::max(); Real WarpX::gamma_boost = 1.; Real WarpX::beta_boost = 0.; @@ -37,13 +36,15 @@ Vector<int> WarpX::boost_direction = {0,0,0}; int WarpX::do_compute_max_step_from_zmax = 0; Real WarpX::zmax_plasma_to_compute_max_step = 0.; -long WarpX::use_picsar_deposition = 1; long WarpX::current_deposition_algo; long WarpX::charge_deposition_algo; long WarpX::field_gathering_algo; long WarpX::particle_pusher_algo; int WarpX::maxwell_fdtd_solver_id; +long WarpX::n_rz_azimuthal_modes = 1; +long WarpX::ncomps = 1; + long WarpX::nox = 1; long WarpX::noy = 1; long WarpX::noz = 1; @@ -68,6 +69,8 @@ bool WarpX::do_boosted_frame_particles = true; bool WarpX::do_dynamic_scheduling = true; +int WarpX::do_subcycling = 0; + #if (AMREX_SPACEDIM == 3) IntVect WarpX::Bx_nodal_flag(1,0,0); IntVect WarpX::By_nodal_flag(0,1,0); @@ -100,7 +103,7 @@ IntVect WarpX::jz_nodal_flag(1,0); // z is the second dimension to 2D AMReX IntVect WarpX::filter_npass_each_dir(1); -int WarpX::n_field_gather_buffer = 0; +int WarpX::n_field_gather_buffer = -1; int WarpX::n_current_deposition_buffer = -1; int WarpX::do_nodal = false; @@ -111,7 +114,7 @@ WarpX& WarpX::GetInstance () { if (!m_instance) { - m_instance = new WarpX(); + m_instance = new WarpX(); } return *m_instance; } @@ -134,20 +137,20 @@ WarpX::WarpX () // No valid BoxArray and DistributionMapping have been defined. // But the arrays for them have been resized. - int nlevs_max = maxLevel() + 1; + const int nlevs_max = maxLevel() + 1; istep.resize(nlevs_max, 0); nsubsteps.resize(nlevs_max, 1); #if 0 // no subcycling yet - for (int lev = 1; lev <= maxLevel(); ++lev) { - nsubsteps[lev] = MaxRefRatio(lev-1); + for (int lev = 1; lev < nlevs_max; ++lev) { + nsubsteps[lev] = MaxRefRatio(lev-1); } #endif t_new.resize(nlevs_max, 0.0); - t_old.resize(nlevs_max, -1.e100); - dt.resize(nlevs_max, 1.e100); + t_old.resize(nlevs_max, std::numeric_limits<Real>::lowest()); + dt.resize(nlevs_max, std::numeric_limits<Real>::max()); // Particle Container mypc = std::unique_ptr<MultiParticleContainer> (new MultiParticleContainer(this)); @@ -187,11 +190,6 @@ WarpX::WarpX () current_buf.resize(nlevs_max); charge_buf.resize(nlevs_max); - current_fp_owner_masks.resize(nlevs_max); - current_cp_owner_masks.resize(nlevs_max); - rho_fp_owner_masks.resize(nlevs_max); - rho_cp_owner_masks.resize(nlevs_max); - pml.resize(nlevs_max); #ifdef WARPX_DO_ELECTROSTATIC @@ -241,7 +239,7 @@ WarpX::WarpX () WarpX::~WarpX () { - int nlevs_max = maxLevel() +1; + const int nlevs_max = maxLevel() +1; for (int lev = 0; lev < nlevs_max; ++lev) { ClearLevel(lev); } @@ -255,133 +253,149 @@ void 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); + ParmParse pp; // Traditionally, max_step and stop_time do not have prefix. + pp.query("max_step", max_step); + pp.query("stop_time", stop_time); } { - ParmParse pp("amr"); // Traditionally, these have prefix, amr. + ParmParse pp("amr"); // Traditionally, these have prefix, amr. - pp.query("check_file", check_file); - pp.query("check_int", check_int); + pp.query("check_file", check_file); + pp.query("check_int", check_int); - pp.query("plot_file", plot_file); - pp.query("plot_int", plot_int); + pp.query("plot_file", plot_file); + pp.query("plot_int", plot_int); - pp.query("restart", restart_chkfile); + pp.query("restart", restart_chkfile); } { - ParmParse pp("warpx"); - - pp.query("cfl", cfl); - pp.query("verbose", verbose); - pp.query("regrid_int", regrid_int); - pp.query("do_subcycling", do_subcycling); - - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_subcycling != 1 || max_level <= 1, - "Subcycling method 1 only works for 2 levels."); - - ReadBoostedFrameParameters(gamma_boost, beta_boost, boost_direction); - - // pp.query returns 1 if argument zmax_plasma_to_compute_max_step is - // specified by the user, 0 otherwise. - do_compute_max_step_from_zmax = - pp.query("zmax_plasma_to_compute_max_step", - zmax_plasma_to_compute_max_step); - - pp.queryarr("B_external", B_external); - - pp.query("do_moving_window", do_moving_window); - if (do_moving_window) - { - std::string s; - pp.get("moving_window_dir", s); - if (s == "x" || s == "X") { - moving_window_dir = 0; - } + ParmParse pp("warpx"); + + pp.query("cfl", cfl); + pp.query("verbose", verbose); + pp.query("regrid_int", regrid_int); + pp.query("do_subcycling", do_subcycling); + pp.query("override_sync_int", override_sync_int); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_subcycling != 1 || max_level <= 1, + "Subcycling method 1 only works for 2 levels."); + + ReadBoostedFrameParameters(gamma_boost, beta_boost, boost_direction); + + // pp.query returns 1 if argument zmax_plasma_to_compute_max_step is + // specified by the user, 0 otherwise. + do_compute_max_step_from_zmax = + pp.query("zmax_plasma_to_compute_max_step", + zmax_plasma_to_compute_max_step); + + pp.queryarr("B_external", B_external); + + pp.query("do_moving_window", do_moving_window); + if (do_moving_window) + { + std::string s; + pp.get("moving_window_dir", s); + if (s == "x" || s == "X") { + moving_window_dir = 0; + } #if (AMREX_SPACEDIM == 3) - else if (s == "y" || s == "Y") { - moving_window_dir = 1; - } + else if (s == "y" || s == "Y") { + moving_window_dir = 1; + } #endif - else if (s == "z" || s == "Z") { - moving_window_dir = AMREX_SPACEDIM-1; - } - else { - const std::string msg = "Unknown moving_window_dir: "+s; - amrex::Abort(msg.c_str()); - } + else if (s == "z" || s == "Z") { + moving_window_dir = AMREX_SPACEDIM-1; + } + else { + const std::string msg = "Unknown moving_window_dir: "+s; + amrex::Abort(msg.c_str()); + } - moving_window_x = geom[0].ProbLo(moving_window_dir); + moving_window_x = geom[0].ProbLo(moving_window_dir); - pp.get("moving_window_v", moving_window_v); - moving_window_v *= PhysConst::c; - } + pp.get("moving_window_v", moving_window_v); + moving_window_v *= PhysConst::c; + } - pp.query("do_boosted_frame_diagnostic", do_boosted_frame_diagnostic); - if (do_boosted_frame_diagnostic) { + pp.query("do_boosted_frame_diagnostic", do_boosted_frame_diagnostic); + if (do_boosted_frame_diagnostic) { - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(gamma_boost > 1.0, - "gamma_boost must be > 1 to use the boosted frame diagnostic."); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(gamma_boost > 1.0, + "gamma_boost must be > 1 to use the boosted frame diagnostic."); - pp.query("lab_data_directory", lab_data_directory); + pp.query("lab_data_directory", lab_data_directory); - std::string s; - pp.get("boost_direction", s); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( (s == "z" || s == "Z"), - "The boosted frame diagnostic currently only works if the boost is in the z direction."); + std::string s; + pp.get("boost_direction", s); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( (s == "z" || s == "Z"), + "The boosted frame diagnostic currently only works if the boost is in the z direction."); - pp.get("num_snapshots_lab", num_snapshots_lab); - pp.get("dt_snapshots_lab", dt_snapshots_lab); - pp.get("gamma_boost", gamma_boost); + pp.get("num_snapshots_lab", num_snapshots_lab); - pp.query("do_boosted_frame_fields", do_boosted_frame_fields); + // Read either dz_snapshots_lab or dt_snapshots_lab + bool snapshot_interval_is_specified = 0; + Real dz_snapshots_lab = 0; + snapshot_interval_is_specified += pp.query("dt_snapshots_lab", dt_snapshots_lab); + if ( pp.query("dz_snapshots_lab", dz_snapshots_lab) ){ + dt_snapshots_lab = dz_snapshots_lab/PhysConst::c; + snapshot_interval_is_specified = 1; + } + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + snapshot_interval_is_specified, + "When using back-transformed diagnostics, user should specify either dz_snapshots_lab or dt_snapshots_lab."); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_moving_window, - "The moving window should be on if using the boosted frame diagnostic."); + pp.get("gamma_boost", gamma_boost); - pp.get("moving_window_dir", s); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( (s == "z" || s == "Z"), - "The boosted frame diagnostic currently only works if the moving window is in the z direction."); - } + pp.query("do_boosted_frame_fields", do_boosted_frame_fields); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_moving_window, + "The moving window should be on if using the boosted frame diagnostic."); + + pp.get("moving_window_dir", s); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( (s == "z" || s == "Z"), + "The boosted frame diagnostic currently only works if the moving window is in the z direction."); + } - pp.query("do_electrostatic", do_electrostatic); - pp.query("n_buffer", n_buffer); - pp.query("const_dt", const_dt); - - // Read filter and fill IntVect filter_npass_each_dir with - // proper size for AMREX_SPACEDIM - pp.query("use_filter", use_filter); - Vector<int> parse_filter_npass_each_dir(AMREX_SPACEDIM,1); - pp.queryarr("filter_npass_each_dir", parse_filter_npass_each_dir); - filter_npass_each_dir[0] = parse_filter_npass_each_dir[0]; - filter_npass_each_dir[1] = parse_filter_npass_each_dir[1]; + pp.query("do_electrostatic", do_electrostatic); + pp.query("n_buffer", n_buffer); + pp.query("const_dt", const_dt); + + // Read filter and fill IntVect filter_npass_each_dir with + // proper size for AMREX_SPACEDIM + pp.query("use_filter", use_filter); + Vector<int> parse_filter_npass_each_dir(AMREX_SPACEDIM,1); + pp.queryarr("filter_npass_each_dir", parse_filter_npass_each_dir); + filter_npass_each_dir[0] = parse_filter_npass_each_dir[0]; + filter_npass_each_dir[1] = parse_filter_npass_each_dir[1]; #if (AMREX_SPACEDIM == 3) - filter_npass_each_dir[2] = parse_filter_npass_each_dir[2]; + filter_npass_each_dir[2] = parse_filter_npass_each_dir[2]; #endif - pp.query("num_mirrors", num_mirrors); - if (num_mirrors>0){ - mirror_z.resize(num_mirrors); - pp.getarr("mirror_z", mirror_z, 0, num_mirrors); - mirror_z_width.resize(num_mirrors); - pp.getarr("mirror_z_width", mirror_z_width, 0, num_mirrors); - mirror_z_npoints.resize(num_mirrors); - pp.getarr("mirror_z_npoints", mirror_z_npoints, 0, num_mirrors); - } + pp.query("num_mirrors", num_mirrors); + if (num_mirrors>0){ + mirror_z.resize(num_mirrors); + pp.getarr("mirror_z", mirror_z, 0, num_mirrors); + mirror_z_width.resize(num_mirrors); + pp.getarr("mirror_z_width", mirror_z_width, 0, num_mirrors); + mirror_z_npoints.resize(num_mirrors); + pp.getarr("mirror_z_npoints", mirror_z_npoints, 0, num_mirrors); + } - pp.query("serialize_ics", serialize_ics); - pp.query("refine_plasma", refine_plasma); + pp.query("serialize_ics", serialize_ics); + pp.query("refine_plasma", refine_plasma); pp.query("do_dive_cleaning", do_dive_cleaning); pp.query("n_field_gather_buffer", n_field_gather_buffer); pp.query("n_current_deposition_buffer", n_current_deposition_buffer); - pp.query("sort_int", sort_int); + pp.query("sort_int", sort_int); pp.query("do_pml", do_pml); pp.query("pml_ncell", pml_ncell); pp.query("pml_delta", pml_delta); + pp.query("pml_has_particles", pml_has_particles); + pp.query("do_pml_j_damping", do_pml_j_damping); + pp.query("do_pml_in_domain", do_pml_in_domain); Vector<int> parse_do_pml_Lo(AMREX_SPACEDIM,1); pp.queryarr("do_pml_Lo", parse_do_pml_Lo); @@ -398,8 +412,12 @@ WarpX::ReadParameters () do_pml_Hi[2] = parse_do_pml_Hi[2]; #endif + if ( (do_pml_j_damping==1)&&(do_pml_in_domain==0) ){ + amrex::Abort("J-damping can only be done when PML are inside simulation domain (do_pml_in_domain=1)"); + } pp.query("dump_openpmd", dump_openpmd); + pp.query("openpmd_backend", openpmd_backend); pp.query("dump_plotfiles", dump_plotfiles); pp.query("plot_raw_fields", plot_raw_fields); pp.query("plot_raw_fields_guards", plot_raw_fields_guards); @@ -434,7 +452,8 @@ WarpX::ReadParameters () } // Check that the coarsening_ratio can divide the blocking factor - for (int lev=0; lev<maxLevel(); lev++){ + const int nlevs_max = maxLevel(); + for (int lev=0; lev<nlevs_max; lev++){ for (int comp=0; comp<AMREX_SPACEDIM; comp++){ if ( blockingFactor(lev)[comp] % plot_coarsening_ratio != 0 ){ amrex::Abort("plot_coarsening_ratio should be an integer " @@ -498,24 +517,23 @@ WarpX::ReadParameters () // Use same shape factors in all directions, for gathering l_lower_order_in_v = false; } + + // Only needs to be set with WARPX_DIM_RZ, otherwise defaults to 1. + pp.query("n_rz_azimuthal_modes", n_rz_azimuthal_modes); } { - ParmParse pp("interpolation"); - pp.query("nox", nox); - pp.query("noy", noy); - pp.query("noz", noz); + ParmParse pp("interpolation"); + pp.query("nox", nox); + pp.query("noy", noy); + pp.query("noz", noz); AMREX_ALWAYS_ASSERT_WITH_MESSAGE( nox == noy and nox == noz , - "warpx.nox, noy and noz must be equal"); + "warpx.nox, noy and noz must be equal"); AMREX_ALWAYS_ASSERT_WITH_MESSAGE( nox >= 1, "warpx.nox must >= 1"); } { ParmParse pp("algo"); - // If not in RZ mode, read use_picsar_deposition - // In RZ mode, use_picsar_deposition is on, as the C++ version - // of the deposition does not support RZ - pp.query("use_picsar_deposition", use_picsar_deposition); current_deposition_algo = GetAlgorithmInteger(pp, "current_deposition"); charge_deposition_algo = GetAlgorithmInteger(pp, "charge_deposition"); field_gathering_algo = GetAlgorithmInteger(pp, "field_gathering"); @@ -601,30 +619,24 @@ void WarpX::ClearLevel (int lev) { for (int i = 0; i < 3; ++i) { - Efield_aux[lev][i].reset(); - Bfield_aux[lev][i].reset(); + Efield_aux[lev][i].reset(); + Bfield_aux[lev][i].reset(); - current_fp[lev][i].reset(); - Efield_fp [lev][i].reset(); - Bfield_fp [lev][i].reset(); + current_fp[lev][i].reset(); + Efield_fp [lev][i].reset(); + Bfield_fp [lev][i].reset(); current_store[lev][i].reset(); - current_cp[lev][i].reset(); - Efield_cp [lev][i].reset(); - Bfield_cp [lev][i].reset(); + current_cp[lev][i].reset(); + Efield_cp [lev][i].reset(); + Bfield_cp [lev][i].reset(); - Efield_cax[lev][i].reset(); - Bfield_cax[lev][i].reset(); + Efield_cax[lev][i].reset(); + Bfield_cax[lev][i].reset(); current_buf[lev][i].reset(); - - current_fp_owner_masks[lev][i].reset(); - current_cp_owner_masks[lev][i].reset(); } - rho_fp_owner_masks[lev].reset(); - rho_cp_owner_masks[lev].reset(); - charge_buf[lev].reset(); current_buffer_masks[lev].reset(); @@ -682,7 +694,7 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d 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::stencil_width; + int ng = ngz_tmp + NCIGodfreyFilter::m_stencil_width; ngz = (ng % 2) ? ng+1 : ng; } else { ngz = ngz_nonci; @@ -719,11 +731,19 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d if (mypc->nSpeciesDepositOnMainGrid() && n_current_deposition_buffer == 0) { n_current_deposition_buffer = 1; + // This forces the allocation of buffers and allows the code associated + // with buffers to run. But the buffer size of `1` is in fact not used, + // `deposit_on_main_grid` forces all particles (whether or not they + // are in buffers) to deposition on the main grid. } if (n_current_deposition_buffer < 0) { n_current_deposition_buffer = ngJ.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 @@ -767,48 +787,51 @@ void WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm, const IntVect& ngE, const IntVect& ngJ, const IntVect& ngRho, int ngF) { + +#if defined WARPX_DIM_RZ + // With RZ multimode, there is a real and imaginary component + // for each mode, except mode 0 which is purely real + // Component 0 is mode 0. + // Odd components are the real parts. + // Even components are the imaginary parts. + ncomps = n_rz_azimuthal_modes*2 - 1; +#endif + // // The fine patch // - Bfield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Bx_nodal_flag),dm,1,ngE)); - Bfield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,By_nodal_flag),dm,1,ngE)); - Bfield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Bz_nodal_flag),dm,1,ngE)); + Bfield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Bx_nodal_flag),dm,ncomps,ngE)); + Bfield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,By_nodal_flag),dm,ncomps,ngE)); + Bfield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Bz_nodal_flag),dm,ncomps,ngE)); - Efield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Ex_nodal_flag),dm,1,ngE)); - Efield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,Ey_nodal_flag),dm,1,ngE)); - Efield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Ez_nodal_flag),dm,1,ngE)); + Efield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Ex_nodal_flag),dm,ncomps,ngE)); + Efield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,Ey_nodal_flag),dm,ncomps,ngE)); + Efield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Ez_nodal_flag),dm,ncomps,ngE)); - current_fp[lev][0].reset( new MultiFab(amrex::convert(ba,jx_nodal_flag),dm,1,ngJ)); - current_fp[lev][1].reset( new MultiFab(amrex::convert(ba,jy_nodal_flag),dm,1,ngJ)); - current_fp[lev][2].reset( new MultiFab(amrex::convert(ba,jz_nodal_flag),dm,1,ngJ)); - - const auto& period = Geom(lev).periodicity(); - current_fp_owner_masks[lev][0] = std::move(current_fp[lev][0]->OwnerMask(period)); - current_fp_owner_masks[lev][1] = std::move(current_fp[lev][1]->OwnerMask(period)); - current_fp_owner_masks[lev][2] = std::move(current_fp[lev][2]->OwnerMask(period)); + current_fp[lev][0].reset( new MultiFab(amrex::convert(ba,jx_nodal_flag),dm,ncomps,ngJ)); + current_fp[lev][1].reset( new MultiFab(amrex::convert(ba,jy_nodal_flag),dm,ncomps,ngJ)); + current_fp[lev][2].reset( new MultiFab(amrex::convert(ba,jz_nodal_flag),dm,ncomps,ngJ)); if (do_dive_cleaning || plot_rho) { - rho_fp[lev].reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,2,ngRho)); - rho_fp_owner_masks[lev] = std::move(rho_fp[lev]->OwnerMask(period)); + rho_fp[lev].reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,2*ncomps,ngRho)); } if (do_subcycling == 1 && lev == 0) { - current_store[lev][0].reset( new MultiFab(amrex::convert(ba,jx_nodal_flag),dm,1,ngJ)); - current_store[lev][1].reset( new MultiFab(amrex::convert(ba,jy_nodal_flag),dm,1,ngJ)); - current_store[lev][2].reset( new MultiFab(amrex::convert(ba,jz_nodal_flag),dm,1,ngJ)); + current_store[lev][0].reset( new MultiFab(amrex::convert(ba,jx_nodal_flag),dm,ncomps,ngJ)); + current_store[lev][1].reset( new MultiFab(amrex::convert(ba,jy_nodal_flag),dm,ncomps,ngJ)); + current_store[lev][2].reset( new MultiFab(amrex::convert(ba,jz_nodal_flag),dm,ncomps,ngJ)); } if (do_dive_cleaning) { - F_fp[lev].reset (new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,1, ngF)); + F_fp[lev].reset (new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,ncomps, ngF)); } #ifdef WARPX_USE_PSATD else { - rho_fp[lev].reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,2,ngRho)); - rho_fp_owner_masks[lev] = std::move(rho_fp[lev]->OwnerMask(period)); + rho_fp[lev].reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,2*ncomps,ngRho)); } if (fft_hybrid_mpi_decomposition == false){ // Allocate and initialize the spectral solver @@ -833,19 +856,19 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm if (lev == 0) { for (int idir = 0; idir < 3; ++idir) { - Efield_aux[lev][idir].reset(new MultiFab(*Efield_fp[lev][idir], amrex::make_alias, 0, 1)); - Bfield_aux[lev][idir].reset(new MultiFab(*Bfield_fp[lev][idir], amrex::make_alias, 0, 1)); + Efield_aux[lev][idir].reset(new MultiFab(*Efield_fp[lev][idir], amrex::make_alias, 0, ncomps)); + Bfield_aux[lev][idir].reset(new MultiFab(*Bfield_fp[lev][idir], amrex::make_alias, 0, ncomps)); } } else { - Bfield_aux[lev][0].reset( new MultiFab(amrex::convert(ba,Bx_nodal_flag),dm,1,ngE)); - Bfield_aux[lev][1].reset( new MultiFab(amrex::convert(ba,By_nodal_flag),dm,1,ngE)); - Bfield_aux[lev][2].reset( new MultiFab(amrex::convert(ba,Bz_nodal_flag),dm,1,ngE)); + Bfield_aux[lev][0].reset( new MultiFab(amrex::convert(ba,Bx_nodal_flag),dm,ncomps,ngE)); + Bfield_aux[lev][1].reset( new MultiFab(amrex::convert(ba,By_nodal_flag),dm,ncomps,ngE)); + Bfield_aux[lev][2].reset( new MultiFab(amrex::convert(ba,Bz_nodal_flag),dm,ncomps,ngE)); - Efield_aux[lev][0].reset( new MultiFab(amrex::convert(ba,Ex_nodal_flag),dm,1,ngE)); - Efield_aux[lev][1].reset( new MultiFab(amrex::convert(ba,Ey_nodal_flag),dm,1,ngE)); - Efield_aux[lev][2].reset( new MultiFab(amrex::convert(ba,Ez_nodal_flag),dm,1,ngE)); + Efield_aux[lev][0].reset( new MultiFab(amrex::convert(ba,Ex_nodal_flag),dm,ncomps,ngE)); + Efield_aux[lev][1].reset( new MultiFab(amrex::convert(ba,Ey_nodal_flag),dm,ncomps,ngE)); + Efield_aux[lev][2].reset( new MultiFab(amrex::convert(ba,Ez_nodal_flag),dm,ncomps,ngE)); } // @@ -857,38 +880,46 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm cba.coarsen(refRatio(lev-1)); // Create the MultiFabs for B - Bfield_cp[lev][0].reset( new MultiFab(amrex::convert(cba,Bx_nodal_flag),dm,1,ngE)); - Bfield_cp[lev][1].reset( new MultiFab(amrex::convert(cba,By_nodal_flag),dm,1,ngE)); - Bfield_cp[lev][2].reset( new MultiFab(amrex::convert(cba,Bz_nodal_flag),dm,1,ngE)); + Bfield_cp[lev][0].reset( new MultiFab(amrex::convert(cba,Bx_nodal_flag),dm,ncomps,ngE)); + Bfield_cp[lev][1].reset( new MultiFab(amrex::convert(cba,By_nodal_flag),dm,ncomps,ngE)); + Bfield_cp[lev][2].reset( new MultiFab(amrex::convert(cba,Bz_nodal_flag),dm,ncomps,ngE)); // Create the MultiFabs for E - Efield_cp[lev][0].reset( new MultiFab(amrex::convert(cba,Ex_nodal_flag),dm,1,ngE)); - Efield_cp[lev][1].reset( new MultiFab(amrex::convert(cba,Ey_nodal_flag),dm,1,ngE)); - Efield_cp[lev][2].reset( new MultiFab(amrex::convert(cba,Ez_nodal_flag),dm,1,ngE)); + Efield_cp[lev][0].reset( new MultiFab(amrex::convert(cba,Ex_nodal_flag),dm,ncomps,ngE)); + Efield_cp[lev][1].reset( new MultiFab(amrex::convert(cba,Ey_nodal_flag),dm,ncomps,ngE)); + Efield_cp[lev][2].reset( new MultiFab(amrex::convert(cba,Ez_nodal_flag),dm,ncomps,ngE)); // Create the MultiFabs for the current - current_cp[lev][0].reset( new MultiFab(amrex::convert(cba,jx_nodal_flag),dm,1,ngJ)); - current_cp[lev][1].reset( new MultiFab(amrex::convert(cba,jy_nodal_flag),dm,1,ngJ)); - current_cp[lev][2].reset( new MultiFab(amrex::convert(cba,jz_nodal_flag),dm,1,ngJ)); - - const auto& cperiod = Geom(lev-1).periodicity(); - current_cp_owner_masks[lev][0] = std::move(current_cp[lev][0]->OwnerMask(cperiod)); - current_cp_owner_masks[lev][1] = std::move(current_cp[lev][1]->OwnerMask(cperiod)); - current_cp_owner_masks[lev][2] = std::move(current_cp[lev][2]->OwnerMask(cperiod)); + current_cp[lev][0].reset( new MultiFab(amrex::convert(cba,jx_nodal_flag),dm,ncomps,ngJ)); + current_cp[lev][1].reset( new MultiFab(amrex::convert(cba,jy_nodal_flag),dm,ncomps,ngJ)); + current_cp[lev][2].reset( new MultiFab(amrex::convert(cba,jz_nodal_flag),dm,ncomps,ngJ)); if (do_dive_cleaning || plot_rho){ - rho_cp[lev].reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2,ngRho)); - rho_cp_owner_masks[lev] = std::move(rho_cp[lev]->OwnerMask(cperiod)); + rho_cp[lev].reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2*ncomps,ngRho)); } if (do_dive_cleaning) { - F_cp[lev].reset (new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,1, ngF)); + F_cp[lev].reset (new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,ncomps, ngF)); } #ifdef WARPX_USE_PSATD else { - rho_cp[lev].reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2,ngRho)); - rho_cp_owner_masks[lev] = std::move(rho_cp[lev]->OwnerMask(cperiod)); + rho_cp[lev].reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2*ncomps,ngRho)); + } + if (fft_hybrid_mpi_decomposition == false){ + // Allocate and initialize the spectral solver + std::array<Real,3> cdx = CellSize(lev-1); + #if (AMREX_SPACEDIM == 3) + RealVect cdx_vect(cdx[0], cdx[1], cdx[2]); + #elif (AMREX_SPACEDIM == 2) + RealVect cdx_vect(cdx[0], cdx[2]); + #endif + // Get the cell-centered box, with guard cells + BoxArray realspace_ba = cba; // Copy box + realspace_ba.enclosedCells().grow(ngE); // cell-centered + guard cells + // Define spectral solver + spectral_solver_cp[lev].reset( new SpectralSolver( realspace_ba, dm, + nox_fft, noy_fft, noz_fft, do_nodal, cdx_vect, dt[lev] ) ); } #endif } @@ -896,35 +927,36 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm // // Copy of the coarse aux // - if (lev > 0 && (n_field_gather_buffer > 0 || n_current_deposition_buffer > 0)) + if (lev > 0 && (n_field_gather_buffer > 0 || n_current_deposition_buffer > 0 || + mypc->nSpeciesGatherFromMainGrid() > 0)) { BoxArray cba = ba; cba.coarsen(refRatio(lev-1)); - if (n_field_gather_buffer > 0) { + if (n_field_gather_buffer > 0 || mypc->nSpeciesGatherFromMainGrid() > 0) { // Create the MultiFabs for B - Bfield_cax[lev][0].reset( new MultiFab(amrex::convert(cba,Bx_nodal_flag),dm,1,ngE)); - Bfield_cax[lev][1].reset( new MultiFab(amrex::convert(cba,By_nodal_flag),dm,1,ngE)); - Bfield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Bz_nodal_flag),dm,1,ngE)); + Bfield_cax[lev][0].reset( new MultiFab(amrex::convert(cba,Bx_nodal_flag),dm,ncomps,ngE)); + Bfield_cax[lev][1].reset( new MultiFab(amrex::convert(cba,By_nodal_flag),dm,ncomps,ngE)); + Bfield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Bz_nodal_flag),dm,ncomps,ngE)); // Create the MultiFabs for E - Efield_cax[lev][0].reset( new MultiFab(amrex::convert(cba,Ex_nodal_flag),dm,1,ngE)); - Efield_cax[lev][1].reset( new MultiFab(amrex::convert(cba,Ey_nodal_flag),dm,1,ngE)); - Efield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Ez_nodal_flag),dm,1,ngE)); + Efield_cax[lev][0].reset( new MultiFab(amrex::convert(cba,Ex_nodal_flag),dm,ncomps,ngE)); + Efield_cax[lev][1].reset( new MultiFab(amrex::convert(cba,Ey_nodal_flag),dm,ncomps,ngE)); + Efield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Ez_nodal_flag),dm,ncomps,ngE)); - gather_buffer_masks[lev].reset( new iMultiFab(ba, dm, 1, 1) ); + gather_buffer_masks[lev].reset( new iMultiFab(ba, dm, ncomps, 1) ); // Gather buffer masks have 1 ghost cell, because of the fact // that particles may move by more than one cell when using subcycling. } if (n_current_deposition_buffer > 0) { - current_buf[lev][0].reset( new MultiFab(amrex::convert(cba,jx_nodal_flag),dm,1,ngJ)); - current_buf[lev][1].reset( new MultiFab(amrex::convert(cba,jy_nodal_flag),dm,1,ngJ)); - current_buf[lev][2].reset( new MultiFab(amrex::convert(cba,jz_nodal_flag),dm,1,ngJ)); - if (do_dive_cleaning || plot_rho) { - charge_buf[lev].reset( new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2,ngRho)); + current_buf[lev][0].reset( new MultiFab(amrex::convert(cba,jx_nodal_flag),dm,ncomps,ngJ)); + current_buf[lev][1].reset( new MultiFab(amrex::convert(cba,jy_nodal_flag),dm,ncomps,ngJ)); + current_buf[lev][2].reset( new MultiFab(amrex::convert(cba,jz_nodal_flag),dm,ncomps,ngJ)); + if (rho_cp[lev]) { + charge_buf[lev].reset( new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2*ncomps,ngRho)); } - current_buffer_masks[lev].reset( new iMultiFab(ba, dm, 1, 1) ); + current_buffer_masks[lev].reset( new iMultiFab(ba, dm, ncomps, 1) ); // Current buffer masks have 1 ghost cell, because of the fact // that particles may move by more than one cell when using subcycling. } @@ -965,7 +997,7 @@ WarpX::LowerCorner(const Box& bx, int lev) #if (AMREX_SPACEDIM == 3) return { xyzmin[0], xyzmin[1], xyzmin[2] }; #elif (AMREX_SPACEDIM == 2) - return { xyzmin[0], -1.e100, xyzmin[1] }; + return { xyzmin[0], std::numeric_limits<Real>::lowest(), xyzmin[1] }; #endif } @@ -977,7 +1009,7 @@ WarpX::UpperCorner(const Box& bx, int lev) #if (AMREX_SPACEDIM == 3) return { xyzmax[0], xyzmax[1], xyzmax[2] }; #elif (AMREX_SPACEDIM == 2) - return { xyzmax[0], 1.e100, xyzmax[1] }; + return { xyzmax[0], std::numeric_limits<Real>::max(), xyzmax[1] }; #endif } |