diff options
Diffstat (limited to 'Source/Initialization')
-rw-r--r-- | Source/Initialization/CustomDensityProb.H | 2 | ||||
-rw-r--r-- | Source/Initialization/InjectorDensity.H | 17 | ||||
-rw-r--r-- | Source/Initialization/InjectorDensity.cpp | 3 | ||||
-rw-r--r-- | Source/Initialization/InjectorMomentum.H | 15 | ||||
-rw-r--r-- | Source/Initialization/InjectorMomentum.cpp | 1 | ||||
-rw-r--r-- | Source/Initialization/InjectorPosition.H | 14 | ||||
-rw-r--r-- | Source/Initialization/PlasmaInjector.H | 16 | ||||
-rw-r--r-- | Source/Initialization/PlasmaInjector.cpp | 17 | ||||
-rw-r--r-- | Source/Initialization/WarpXInitData.cpp | 51 |
9 files changed, 76 insertions, 60 deletions
diff --git a/Source/Initialization/CustomDensityProb.H b/Source/Initialization/CustomDensityProb.H index b00830e6c..c5c159ee8 100644 --- a/Source/Initialization/CustomDensityProb.H +++ b/Source/Initialization/CustomDensityProb.H @@ -27,7 +27,7 @@ struct InjectorDensityCustom } } - // Return density at given position, using user-defined parameters + // Return density at given position, using user-defined parameters // stored in p. AMREX_GPU_HOST_DEVICE amrex::Real diff --git a/Source/Initialization/InjectorDensity.H b/Source/Initialization/InjectorDensity.H index b7f5c26eb..7e61ae27d 100644 --- a/Source/Initialization/InjectorDensity.H +++ b/Source/Initialization/InjectorDensity.H @@ -1,12 +1,13 @@ #ifndef INJECTOR_DENSITY_H_ #define INJECTOR_DENSITY_H_ -#include <AMReX_Gpu.H> -#include <AMReX_Dim3.H> #include <GpuParser.H> #include <CustomDensityProb.H> #include <WarpXConst.H> +#include <AMReX_Gpu.H> +#include <AMReX_Dim3.H> + // struct whose getDensity returns constant density. struct InjectorDensityConstant { @@ -67,18 +68,20 @@ struct InjectorDensityPredefined amrex::Real kp = PhysConst::q_e/PhysConst::c *std::sqrt( n0/(PhysConst::m_e*PhysConst::ep0) ); + // Longitudinal profile, normalized to 1 if ((z-z_start)>=0 and (z-z_start)<ramp_up ) { - n = (z-z_start)/ramp_up; + n = 0.5*(1.-std::cos(MathConst::pi*(z-z_start)/ramp_up)); } else if ((z-z_start)>=ramp_up and (z-z_start)< ramp_up+plateau ) { n = 1.; } else if ((z-z_start)>=ramp_up+plateau and (z-z_start)< ramp_up+plateau+ramp_down) { - n = 1.-((z-z_start)-ramp_up-plateau)/ramp_down; + n = 0.5*(1.+std::cos(MathConst::pi*((z-z_start)-ramp_up-plateau)/ramp_down)); } else { n = 0.; } + // Multiply by transverse profile, and physical density n *= n0*(1.+4.*(x*x+y*y)/(kp*kp*rc*rc*rc*rc)); return n; } @@ -94,9 +97,9 @@ private: amrex::Real* p; }; -// Base struct for density injector. -// InjectorDensity contains a union (called Object) that holds any one -// instance of: +// Base struct for density injector. +// InjectorDensity contains a union (called Object) that holds any one +// instance of: // - InjectorDensityConstant : to generate constant density; // - InjectorDensityParser : to generate density from parser; // - InjectorDensityCustom : to generate density from custom profile; diff --git a/Source/Initialization/InjectorDensity.cpp b/Source/Initialization/InjectorDensity.cpp index 54df4b14d..ddf6d1a18 100644 --- a/Source/Initialization/InjectorDensity.cpp +++ b/Source/Initialization/InjectorDensity.cpp @@ -1,3 +1,4 @@ +#include <InjectorDensity.H> #include <PlasmaInjector.H> using namespace amrex; @@ -48,7 +49,7 @@ InjectorDensityPredefined::InjectorDensityPredefined ( ParmParse pp(a_species_name); std::vector<amrex::Real> v; - // Read parameters for the predefined plasma profile, + // Read parameters for the predefined plasma profile, // and store them in managed memory pp.getarr("predefined_profile_params", v); p = static_cast<amrex::Real*> diff --git a/Source/Initialization/InjectorMomentum.H b/Source/Initialization/InjectorMomentum.H index 399ee7759..ba0d26fc4 100644 --- a/Source/Initialization/InjectorMomentum.H +++ b/Source/Initialization/InjectorMomentum.H @@ -1,10 +1,11 @@ #ifndef INJECTOR_MOMENTUM_H_ #define INJECTOR_MOMENTUM_H_ +#include <CustomMomentumProb.H> +#include <GpuParser.H> + #include <AMReX_Gpu.H> #include <AMReX_Dim3.H> -#include <GpuParser.H> -#include <CustomMomentumProb.H> // struct whose getMomentum returns constant momentum. struct InjectorMomentumConstant @@ -22,7 +23,7 @@ private: amrex::Real m_ux, m_uy, m_uz; }; -// struct whose getMomentum returns momentum for 1 particle, from random +// struct whose getMomentum returns momentum for 1 particle, from random // gaussian distribution. struct InjectorMomentumGaussian { @@ -84,9 +85,9 @@ struct InjectorMomentumParser GpuParser m_ux_parser, m_uy_parser, m_uz_parser; }; -// Base struct for momentum injector. -// InjectorMomentum contains a union (called Object) that holds any one -// instance of: +// Base struct for momentum injector. +// InjectorMomentum contains a union (called Object) that holds any one +// instance of: // - InjectorMomentumConstant : to generate constant density; // - InjectorMomentumGaussian : to generate gaussian distribution; // - InjectorMomentumRadialExpansion: to generate radial expansion; @@ -189,7 +190,7 @@ private: Type type; // An instance of union Object constructs and stores any one of - // the objects declared (constant or custom or gaussian or + // the objects declared (constant or custom or gaussian or // radial_expansion or parser). union Object { Object (InjectorMomentumConstant*, diff --git a/Source/Initialization/InjectorMomentum.cpp b/Source/Initialization/InjectorMomentum.cpp index a197b5bef..8fadf0c4b 100644 --- a/Source/Initialization/InjectorMomentum.cpp +++ b/Source/Initialization/InjectorMomentum.cpp @@ -1,3 +1,4 @@ +#include <InjectorMomentum.H> #include <PlasmaInjector.H> using namespace amrex; diff --git a/Source/Initialization/InjectorPosition.H b/Source/Initialization/InjectorPosition.H index 19bb092dd..6ecae93e0 100644 --- a/Source/Initialization/InjectorPosition.H +++ b/Source/Initialization/InjectorPosition.H @@ -25,7 +25,7 @@ struct InjectorPositionRegular // i_part: particle number within the cell, required to evenly space // particles within the cell. - // ref_fac: the number of particles evenly-spaced within a cell + // ref_fac: the number of particles evenly-spaced within a cell // is a_ppc*(ref_fac**AMREX_SPACEDIM). AMREX_GPU_HOST_DEVICE amrex::XDim3 @@ -33,7 +33,7 @@ struct InjectorPositionRegular { int nx = ref_fac*ppc.x; int ny = ref_fac*ppc.y; -#if (AMREX_SPACEDIM == 3) +#if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) int nz = ref_fac*ppc.z; #else int nz = 1; @@ -41,15 +41,17 @@ struct InjectorPositionRegular int ix_part = i_part/(ny*nz); // written this way backward compatibility int iz_part = (i_part-ix_part*(ny*nz)) / ny; int iy_part = (i_part-ix_part*(ny*nz)) - ny*iz_part; - return amrex::XDim3{(0.5+ix_part)/nx, (0.5+iy_part)/ny, (0.5+iz_part) / nz}; + return amrex::XDim3{(amrex::Real(0.5)+ix_part)/nx, + (amrex::Real(0.5)+iy_part)/ny, + (amrex::Real(0.5)+iz_part) / nz}; } private: amrex::Dim3 ppc; }; -// Base struct for position injector. -// InjectorPosition contains a union (called Object) that holds any one -// instance of: +// Base struct for position injector. +// InjectorPosition contains a union (called Object) that holds any one +// instance of: // - InjectorPositionRandom : to generate random distribution; // - InjectorPositionRegular: to generate regular distribution. // The choice is made at runtime, depending in the constructor called. diff --git a/Source/Initialization/PlasmaInjector.H b/Source/Initialization/PlasmaInjector.H index f7e86bff5..56b32c827 100644 --- a/Source/Initialization/PlasmaInjector.H +++ b/Source/Initialization/PlasmaInjector.H @@ -5,13 +5,15 @@ #include <InjectorDensity.H> #include <InjectorMomentum.H> -#include <array> -#include <AMReX_Vector.H> #include <WarpXConst.H> #include <WarpXParser.H> + +#include <AMReX_Vector.H> #include <AMReX_ParmParse.H> #include <AMReX_Utility.H> +#include <array> + /// /// The PlasmaInjector class parses and stores information about the plasma /// type used in the particle container. This information is used to create the @@ -41,9 +43,9 @@ public: bool doInjection () const noexcept { return inj_pos != NULL;} bool add_single_particle = false; - amrex::Vector<amrex::Real> single_particle_pos; - amrex::Vector<amrex::Real> single_particle_vel; - amrex::Real single_particle_weight; + amrex::Vector<amrex::ParticleReal> single_particle_pos; + amrex::Vector<amrex::ParticleReal> single_particle_vel; + amrex::ParticleReal single_particle_weight; bool gaussian_beam = false; amrex::Real x_m; @@ -94,9 +96,9 @@ protected: std::unique_ptr<InjectorPosition> inj_pos; std::unique_ptr<InjectorDensity > inj_rho; std::unique_ptr<InjectorMomentum> inj_mom; - + void parseDensity (amrex::ParmParse& pp); - void parseMomentum (amrex::ParmParse& pp); + void parseMomentum (amrex::ParmParse& pp); }; #endif diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index 541999789..af04c6b31 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -1,13 +1,14 @@ #include "PlasmaInjector.H" -#include <sstream> -#include <functional> - #include <WarpXConst.H> #include <WarpX_f.H> -#include <AMReX.H> #include <WarpX.H> +#include <AMReX.H> + +#include <sstream> +#include <functional> + using namespace amrex; namespace { @@ -45,7 +46,7 @@ namespace { } else if (name == "m_p"){ return PhysConst::m_p; } else if (name == "inf"){ - return std::numeric_limits<double>::infinity(); + return std::numeric_limits<double>::infinity(); } else if (pp.query("mass", result)) { return result; } else { @@ -143,9 +144,11 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) parseDensity(pp); parseMomentum(pp); } else if (part_pos_s == "nuniformpercell") { - num_particles_per_cell_each_dim.resize(3); + // Note that for RZ, three numbers are expected, r, theta, and z. + // For 2D, only two are expected. The third is overwritten with 1. + num_particles_per_cell_each_dim.assign(3, 1); pp.getarr("num_particles_per_cell_each_dim", num_particles_per_cell_each_dim); -#if ( AMREX_SPACEDIM == 2 ) +#if WARPX_DIM_XZ num_particles_per_cell_each_dim[2] = 1; #endif // Construct InjectorPosition from InjectorPositionRegular. diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 590c11b84..385993f78 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -1,12 +1,12 @@ -#include <AMReX_ParallelDescriptor.H> -#include <AMReX_ParmParse.H> - #include <WarpX.H> #include <WarpX_f.H> #include <BilinearFilter.H> #include <NCIGodfreyFilter.H> +#include <AMReX_ParallelDescriptor.H> +#include <AMReX_ParmParse.H> + #ifdef BL_USE_SENSEI_INSITU #include <AMReX_AmrMeshInSituBridge.H> #endif @@ -25,11 +25,11 @@ WarpX::InitData () } else { - InitFromCheckpoint(); + InitFromCheckpoint(); if (is_synchronized) { ComputeDt(); } - PostRestart(); + PostRestart(); } ComputePMLFactors(); @@ -87,12 +87,12 @@ WarpX::InitDiagnostics () { const Real* current_hi = geom[0].ProbHi(); Real dt_boost = dt[0]; - // Find the positions of the lab-frame box that corresponds to the boosted-frame box at t=0 - Real zmin_lab = current_lo[moving_window_dir]/( (1.+beta_boost)*gamma_boost ); - Real zmax_lab = current_hi[moving_window_dir]/( (1.+beta_boost)*gamma_boost ); + // Find the positions of the lab-frame box that corresponds to the boosted-frame box at t=0 + Real zmin_lab = current_lo[moving_window_dir]/( (1.+beta_boost)*gamma_boost ); + Real zmax_lab = current_hi[moving_window_dir]/( (1.+beta_boost)*gamma_boost ); myBFD.reset(new BoostedFrameDiagnostic(zmin_lab, - zmax_lab, + zmax_lab, moving_window_v, dt_snapshots_lab, num_snapshots_lab, gamma_boost, t_new[0], dt_boost, @@ -142,6 +142,7 @@ WarpX::InitPML () dt[0], nox_fft, noy_fft, noz_fft, do_nodal, #endif do_dive_cleaning, do_moving_window, + pml_has_particles, do_pml_in_domain, do_pml_Lo_corrected, do_pml_Hi)); for (int lev = 1; lev <= finest_level; ++lev) { @@ -159,6 +160,7 @@ WarpX::InitPML () dt[lev], nox_fft, noy_fft, noz_fft, do_nodal, #endif do_dive_cleaning, do_moving_window, + pml_has_particles, do_pml_in_domain, do_pml_Lo_MR, amrex::IntVect::TheUnitVector())); } } @@ -195,9 +197,10 @@ WarpX::InitNCICorrector () // Initialize Godfrey filters // Same filter for fields Ex, Ey and Bz - nci_godfrey_filter_exeybz[lev].reset( new NCIGodfreyFilter(godfrey_coeff_set::Ex_Ey_Bz, cdtodz, WarpX::l_lower_order_in_v) ); + const bool nodal_gather = (l_lower_order_in_v == 0); + nci_godfrey_filter_exeybz[lev].reset( new NCIGodfreyFilter(godfrey_coeff_set::Ex_Ey_Bz, cdtodz, nodal_gather) ); // Same filter for fields Bx, By and Ez - nci_godfrey_filter_bxbyez[lev].reset( new NCIGodfreyFilter(godfrey_coeff_set::Bx_By_Ez, cdtodz, WarpX::l_lower_order_in_v) ); + nci_godfrey_filter_bxbyez[lev].reset( new NCIGodfreyFilter(godfrey_coeff_set::Bx_By_Ez, cdtodz, nodal_gather) ); // Compute Godfrey filters stencils nci_godfrey_filter_exeybz[lev]->ComputeStencils(); nci_godfrey_filter_bxbyez[lev]->ComputeStencils(); @@ -249,9 +252,9 @@ WarpX::InitOpenbc () BoxList bl{IndexType::TheNodeType()}; for (int i = 0; i < nprocs; ++i) { - bl.push_back(Box(IntVect(alllohi[6*i ],alllohi[6*i+1],alllohi[6*i+2]), - IntVect(alllohi[6*i+3],alllohi[6*i+4],alllohi[6*i+5]), - IndexType::TheNodeType())); + bl.push_back(Box(IntVect(alllohi[6*i ],alllohi[6*i+1],alllohi[6*i+2]), + IntVect(alllohi[6*i+3],alllohi[6*i+4],alllohi[6*i+5]), + IndexType::TheNodeType())); } BoxArray ba{bl}; @@ -284,13 +287,13 @@ WarpX::InitOpenbc () #endif for (MFIter mfi(phi); mfi.isValid(); ++mfi) { - const Box& bx = mfi.validbox(); - warpx_compute_E(bx.loVect(), bx.hiVect(), - BL_TO_FORTRAN_3D(phi[mfi]), - BL_TO_FORTRAN_3D((*Efield[lev][0])[mfi]), - BL_TO_FORTRAN_3D((*Efield[lev][1])[mfi]), - BL_TO_FORTRAN_3D((*Efield[lev][2])[mfi]), - dx); + const Box& bx = mfi.validbox(); + warpx_compute_E(bx.loVect(), bx.hiVect(), + BL_TO_FORTRAN_3D(phi[mfi]), + BL_TO_FORTRAN_3D((*Efield[lev][0])[mfi]), + BL_TO_FORTRAN_3D((*Efield[lev][1])[mfi]), + BL_TO_FORTRAN_3D((*Efield[lev][2])[mfi]), + dx); } } #endif @@ -299,9 +302,9 @@ void WarpX::InitLevelData (int lev, Real time) { for (int i = 0; i < 3; ++i) { - current_fp[lev][i]->setVal(0.0); - Efield_fp[lev][i]->setVal(0.0); - Bfield_fp[lev][i]->setVal(0.0); + current_fp[lev][i]->setVal(0.0); + Efield_fp[lev][i]->setVal(0.0); + Bfield_fp[lev][i]->setVal(0.0); } if (lev > 0) { |