From dedbed9fca7f9b325d1a1e2768d1f275af3a5bdd Mon Sep 17 00:00:00 2001 From: Edoardo Zoni Date: Thu, 6 Mar 2025 15:31:28 -0800 Subject: [PATCH] Apply clang-format to WarpX.H, WarpX.cpp --- .pre-commit-config.yaml | 2 +- Source/WarpX.H | 1404 ++++++++++------- Source/WarpX.cpp | 3181 ++++++++++++++++++++++----------------- 3 files changed, 2603 insertions(+), 1984 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 078d4802b7f..a5459c67f46 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -68,7 +68,7 @@ repos: rev: v19.1.7 hooks: - id: clang-format - files: '^Source/main.cpp' + files: '(Source/main.cpp|Source/WarpX.H|Source/WarpX.cpp)' # Python: Ruff linter & formatter # https://docs.astral.sh/ruff/ diff --git a/Source/WarpX.H b/Source/WarpX.H index 56dfa3e2f10..4f965dcc949 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -19,31 +19,31 @@ #include "EmbeddedBoundary/WarpXFaceInfoBox_fwd.H" #include "FieldSolver/ElectrostaticSolvers/ElectrostaticSolver_fwd.H" #include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver_fwd.H" -#include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties_fwd.H" #include "FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel_fwd.H" +#include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties_fwd.H" #include "Filter/NCIGodfreyFilter_fwd.H" +#include "Fluids/MultiFluidContainer_fwd.H" +#include "Fluids/WarpXFluidContainer_fwd.H" #include "Initialization/ExternalField_fwd.H" -#include "Particles/ParticleBoundaryBuffer_fwd.H" #include "Particles/MultiParticleContainer_fwd.H" +#include "Particles/ParticleBoundaryBuffer_fwd.H" #include "Particles/WarpXParticleContainer_fwd.H" -#include "Fluids/MultiFluidContainer_fwd.H" -#include "Fluids/WarpXFluidContainer_fwd.H" #ifdef WARPX_USE_FFT -# ifdef WARPX_DIM_RZ -# include "FieldSolver/SpectralSolver/SpectralSolverRZ_fwd.H" -# include "BoundaryConditions/PML_RZ_fwd.H" -# else -# include "FieldSolver/SpectralSolver/SpectralSolver_fwd.H" -# endif +#ifdef WARPX_DIM_RZ +#include "BoundaryConditions/PML_RZ_fwd.H" +#include "FieldSolver/SpectralSolver/SpectralSolverRZ_fwd.H" +#else +#include "FieldSolver/SpectralSolver/SpectralSolver_fwd.H" +#endif #endif #include "AcceleratorLattice/AcceleratorLattice.H" #include "Evolve/WarpXDtType.H" #include "Evolve/WarpXPushType.H" -#include "Fields.H" -#include "FieldSolver/MagnetostaticSolver/MagnetostaticSolver.H" #include "FieldSolver/ImplicitSolvers/ImplicitSolver.H" #include "FieldSolver/ImplicitSolvers/WarpXSolverVec.H" +#include "FieldSolver/MagnetostaticSolver/MagnetostaticSolver.H" +#include "Fields.H" #include "Filter/BilinearFilter.H" #include "Parallelization/GuardCellManager.H" #include "Utils/Parser/IntervalsParser.H" @@ -58,7 +58,7 @@ #include #include #ifdef AMREX_USE_EB -# include +#include #endif #include #include @@ -70,8 +70,8 @@ #include #include -#include #include +#include #include #include @@ -82,45 +82,47 @@ #include #include -class WARPX_EXPORT WarpX - : public amrex::AmrCore -{ -public: +class WARPX_EXPORT WarpX : public amrex::AmrCore { + public: static WarpX& GetInstance (); static void ResetInstance (); /** * \brief - * This method has to be called at the end of the simulation. It deletes the WarpX instance. + * This method has to be called at the end of the simulation. It deletes the + * WarpX instance. */ - static void Finalize(); + static void Finalize (); /** Destructor */ ~WarpX () override; /** Copy constructor */ - WarpX ( WarpX const &) = delete; + WarpX (WarpX const&) = delete; /** Copy operator */ - WarpX& operator= ( WarpX const & ) = delete; + WarpX& operator=(WarpX const&) = delete; /** Move constructor */ - WarpX ( WarpX && ) = delete; + WarpX (WarpX&&) = delete; /** Move operator */ - WarpX& operator= ( WarpX && ) = delete; + WarpX& operator=(WarpX&&) = delete; - static std::string Version (); //!< Version of WarpX executable + static std::string Version (); //!< Version of WarpX executable static std::string PicsarVersion (); //!< Version of PICSAR dependency - [[nodiscard]] int Verbose () const { return verbose; } + [[nodiscard]] int + Verbose () const { + return verbose; + } - [[nodiscard]] const amrex::Array& GetFieldBoundaryLo () const - { + [[nodiscard]] const amrex::Array& + GetFieldBoundaryLo () const { return field_boundary_lo; } - [[nodiscard]] const amrex::Array& GetFieldBoundaryHi () const - { + [[nodiscard]] const amrex::Array& + GetFieldBoundaryHi () const { return field_boundary_hi; } @@ -136,112 +138,160 @@ public: // // Functions used by implicit solvers // - void ImplicitPreRHSOp ( amrex::Real cur_time, - amrex::Real a_full_dt, - int a_nl_iter, - bool a_from_jacobian ); + void ImplicitPreRHSOp (amrex::Real cur_time, amrex::Real a_full_dt, + int a_nl_iter, bool a_from_jacobian); void SaveParticlesAtImplicitStepStart (); void FinishImplicitParticleUpdate (); - void SetElectricFieldAndApplyBCs ( const WarpXSolverVec& a_E, amrex::Real a_time ); - void UpdateMagneticFieldAndApplyBCs ( ablastr::fields::MultiLevelVectorField const& a_Bn, - amrex::Real a_thetadt, amrex::Real start_time ); - void SpectralSourceFreeFieldAdvance ( amrex::Real start_time); - void FinishMagneticFieldAndApplyBCs ( ablastr::fields::MultiLevelVectorField const& a_Bn, - amrex::Real a_theta, amrex::Real a_time ); - void FinishImplicitField ( const ablastr::fields::MultiLevelVectorField& Field_fp, - const ablastr::fields::MultiLevelVectorField& Field_n, - amrex::Real theta ); - void ImplicitComputeRHSE ( amrex::Real dt, WarpXSolverVec& a_Erhs_vec); - void ImplicitComputeRHSE (int lev, amrex::Real dt, WarpXSolverVec& a_Erhs_vec); - void ImplicitComputeRHSE (int lev, PatchType patch_type, amrex::Real dt, WarpXSolverVec& a_Erhs_vec); - - MultiParticleContainer& GetPartContainer () { return *mypc; } - MultiFluidContainer& GetFluidContainer () { return *myfl; } - ElectrostaticSolver& GetElectrostaticSolver () {return *m_electrostatic_solver;} - [[nodiscard]] HybridPICModel * get_pointer_HybridPICModel () const { return m_hybrid_pic_model.get(); } - MultiDiagnostics& GetMultiDiags () {return *multi_diags;} - ParticleBoundaryBuffer& GetParticleBoundaryBuffer () { return *m_particle_boundary_buffer; } - amrex::Vector,3 > >& GetEBUpdateEFlag() { return m_eb_update_E; } - amrex::Vector,3 > >& GetEBUpdateBFlag() { return m_eb_update_B; } - amrex::Vector< std::unique_ptr > const & GetEBReduceParticleShapeFlag() const { return m_eb_reduce_particle_shape; } + void SetElectricFieldAndApplyBCs (const WarpXSolverVec& a_E, + amrex::Real a_time); + void UpdateMagneticFieldAndApplyBCs ( + ablastr::fields::MultiLevelVectorField const& a_Bn, + amrex::Real a_thetadt, amrex::Real start_time); + void SpectralSourceFreeFieldAdvance (amrex::Real start_time); + void FinishMagneticFieldAndApplyBCs ( + ablastr::fields::MultiLevelVectorField const& a_Bn, amrex::Real a_theta, + amrex::Real a_time); + void + FinishImplicitField (const ablastr::fields::MultiLevelVectorField& Field_fp, + const ablastr::fields::MultiLevelVectorField& Field_n, + amrex::Real theta); + void ImplicitComputeRHSE (amrex::Real dt, WarpXSolverVec& a_Erhs_vec); + void ImplicitComputeRHSE (int lev, amrex::Real dt, + WarpXSolverVec& a_Erhs_vec); + void ImplicitComputeRHSE (int lev, PatchType patch_type, amrex::Real dt, + WarpXSolverVec& a_Erhs_vec); + + MultiParticleContainer& + GetPartContainer () { + return *mypc; + } + MultiFluidContainer& + GetFluidContainer () { + return *myfl; + } + ElectrostaticSolver& + GetElectrostaticSolver () { + return *m_electrostatic_solver; + } + [[nodiscard]] HybridPICModel* + get_pointer_HybridPICModel () const { + return m_hybrid_pic_model.get(); + } + MultiDiagnostics& + GetMultiDiags () { + return *multi_diags; + } + ParticleBoundaryBuffer& + GetParticleBoundaryBuffer () { + return *m_particle_boundary_buffer; + } + amrex::Vector, 3>>& + GetEBUpdateEFlag () { + return m_eb_update_E; + } + amrex::Vector, 3>>& + GetEBUpdateBFlag () { + return m_eb_update_B; + } + amrex::Vector> const& + GetEBReduceParticleShapeFlag () const { + return m_eb_reduce_particle_shape; + } /** * \brief - * If an authors' string is specified in the inputfile, this method returns that string. - * Otherwise, it returns an empty string. + * If an authors' string is specified in the inputfile, this method returns + * that string. Otherwise, it returns an empty string. */ - [[nodiscard]] std::string GetAuthors () const { return m_authors; } + [[nodiscard]] std::string + GetAuthors () const { + return m_authors; + } - /** Maximum level up to which the externally defined electric and magnetic fields are initialized. - * The default value is set to the max levels in the simulation. - * if lev > maxlevel_extEMfield_init, the fields on those levels will have a default value of 0 + /** Maximum level up to which the externally defined electric and magnetic + * fields are initialized. The default value is set to the max levels in the + * simulation. if lev > maxlevel_extEMfield_init, the fields on those levels + * will have a default value of 0 */ int maxlevel_extEMfield_init; // Algorithms - //! Integer that corresponds to the current deposition algorithm (Esirkepov, direct, Vay, Villasenor) + //! Integer that corresponds to the current deposition algorithm (Esirkepov, + //! direct, Vay, Villasenor) static inline auto current_deposition_algo = CurrentDepositionAlgo::Default; - //! Integer that corresponds to the charge deposition algorithm (only standard deposition) + //! Integer that corresponds to the charge deposition algorithm (only + //! standard deposition) static inline auto charge_deposition_algo = ChargeDepositionAlgo::Default; - //! Integer that corresponds to the field gathering algorithm (energy-conserving, momentum-conserving) + //! Integer that corresponds to the field gathering algorithm + //! (energy-conserving, momentum-conserving) static inline auto field_gathering_algo = GatheringAlgo::Default; - //! Integer that corresponds to the particle push algorithm (Boris, Vay, Higuera-Cary) + //! Integer that corresponds to the particle push algorithm (Boris, Vay, + //! Higuera-Cary) static inline auto particle_pusher_algo = ParticlePusherAlgo::Default; - //! Integer that corresponds to the type of Maxwell solver (Yee, CKC, PSATD, ECT) - static inline auto electromagnetic_solver_id = ElectromagneticSolverAlgo::Default; - //! Integer that corresponds to the evolve scheme (explicit, semi_implicit_em, theta_implicit_em) + //! Integer that corresponds to the type of Maxwell solver (Yee, CKC, PSATD, + //! ECT) + static inline auto electromagnetic_solver_id = + ElectromagneticSolverAlgo::Default; + //! Integer that corresponds to the evolve scheme (explicit, + //! semi_implicit_em, theta_implicit_em) EvolveScheme evolve_scheme = EvolveScheme::Default; - //! Maximum iterations used for self-consistent particle update in implicit particle-suppressed evolve schemes + //! Maximum iterations used for self-consistent particle update in implicit + //! particle-suppressed evolve schemes static int max_particle_its_in_implicit_scheme; - //! Relative tolerance used for self-consistent particle update in implicit particle-suppressed evolve schemes + //! Relative tolerance used for self-consistent particle update in implicit + //! particle-suppressed evolve schemes static amrex::ParticleReal particle_tol_in_implicit_scheme; /** Records a number corresponding to the load balance cost update strategy * being used (0 or 1 corresponding to timers or heuristic). */ - static inline auto load_balance_costs_update_algo = LoadBalanceCostsUpdateAlgo::Default; + static inline auto load_balance_costs_update_algo = + LoadBalanceCostsUpdateAlgo::Default; /** Integer that correspond to macroscopic Maxwell solver algorithm * (BackwardEuler - 0, Lax-Wendroff - 1) */ static inline auto macroscopic_solver_algo = MacroscopicSolverAlgo::Default; /** Boundary conditions applied to fields at the lower domain boundaries - * (Possible values PML, Periodic, PEC, PMC, Neumann, Damped, Absorbing Silver-Mueller, None) + * (Possible values PML, Periodic, PEC, PMC, Neumann, Damped, Absorbing + * Silver-Mueller, None) */ - static inline amrex::Array - field_boundary_lo {AMREX_D_DECL(FieldBoundaryType::Default, - FieldBoundaryType::Default, - FieldBoundaryType::Default)}; + static inline amrex::Array + field_boundary_lo{AMREX_D_DECL(FieldBoundaryType::Default, + FieldBoundaryType::Default, + FieldBoundaryType::Default)}; /** Boundary conditions applied to fields at the upper domain boundaries - * (Possible values PML, Periodic, PEC, PMC, Neumann, Damped, Absorbing Silver-Mueller, None) - */ - static inline amrex::Array - field_boundary_hi {AMREX_D_DECL(FieldBoundaryType::Default, - FieldBoundaryType::Default, - FieldBoundaryType::Default)}; - /** Integers that correspond to boundary condition applied to particles at the - * lower domain boundaries - * (0 to 4 correspond to Absorbing, Open, Reflecting, Periodic, Thermal) - */ - static inline amrex::Array - particle_boundary_lo {AMREX_D_DECL(ParticleBoundaryType::Default, - ParticleBoundaryType::Default, - ParticleBoundaryType::Default)}; - /** Integers that correspond to boundary condition applied to particles at the - * upper domain boundaries - * (0 to 4 correspond to Absorbing, Open, Reflecting, Periodic, Thermal) - */ - static inline amrex::Array - particle_boundary_hi {AMREX_D_DECL(ParticleBoundaryType::Default, - ParticleBoundaryType::Default, - ParticleBoundaryType::Default)}; + * (Possible values PML, Periodic, PEC, PMC, Neumann, Damped, Absorbing + * Silver-Mueller, None) + */ + static inline amrex::Array + field_boundary_hi{AMREX_D_DECL(FieldBoundaryType::Default, + FieldBoundaryType::Default, + FieldBoundaryType::Default)}; + /** Integers that correspond to boundary condition applied to particles at + * the lower domain boundaries (0 to 4 correspond to Absorbing, Open, + * Reflecting, Periodic, Thermal) + */ + static inline amrex::Array + particle_boundary_lo{AMREX_D_DECL(ParticleBoundaryType::Default, + ParticleBoundaryType::Default, + ParticleBoundaryType::Default)}; + /** Integers that correspond to boundary condition applied to particles at + * the upper domain boundaries (0 to 4 correspond to Absorbing, Open, + * Reflecting, Periodic, Thermal) + */ + static inline amrex::Array + particle_boundary_hi{AMREX_D_DECL(ParticleBoundaryType::Default, + ParticleBoundaryType::Default, + ParticleBoundaryType::Default)}; //! Integers that correspond to the time dependency of J (constant, linear) //! and rho (linear, quadratic) for the PSATD algorithm static inline auto J_in_time = JInTime::Default; static inline auto rho_in_time = RhoInTime::Default; - //! If true, the current is deposited on a nodal grid and then centered onto a staggered grid - //! using finite centering of order given by #current_centering_nox, #current_centering_noy, - //! and #current_centering_noz + //! If true, the current is deposited on a nodal grid and then centered onto + //! a staggered grid using finite centering of order given by + //! #current_centering_nox, #current_centering_noy, and + //! #current_centering_noz bool do_current_centering = false; //! If true, a correction is applied to the current in Fourier space, @@ -273,10 +323,12 @@ public: //! Whether to fill guard cells when computing inverse FFTs of currents static amrex::IntVect m_fill_guards_current; - //! Solve additional Maxwell equation for F in order to control errors in Gauss' law - //! (useful when using current deposition algorithms that are not charge-conserving) + //! Solve additional Maxwell equation for F in order to control errors in + //! Gauss' law (useful when using current deposition algorithms that are not + //! charge-conserving) static bool do_dive_cleaning; - //! Solve additional Maxwell equation for G in order to control errors in magnetic Gauss' law + //! Solve additional Maxwell equation for G in order to control errors in + //! magnetic Gauss' law static bool do_divb_cleaning; //! Order of the particle shape factors (splines) along x @@ -286,64 +338,79 @@ public: //! Order of the particle shape factors (splines) along z static int noz; - //! Order of finite centering of fields (from staggered grid to nodal grid), along x + //! Order of finite centering of fields (from staggered grid to nodal grid), + //! along x static int field_centering_nox; - //! Order of finite centering of fields (from staggered grid to nodal grid), along y + //! Order of finite centering of fields (from staggered grid to nodal grid), + //! along y static int field_centering_noy; - //! Order of finite centering of fields (from staggered grid to nodal grid), along z + //! Order of finite centering of fields (from staggered grid to nodal grid), + //! along z static int field_centering_noz; - //! Order of finite centering of currents (from nodal grid to staggered grid), along x + //! Order of finite centering of currents (from nodal grid to staggered + //! grid), along x static int current_centering_nox; - //! Order of finite centering of currents (from nodal grid to staggered grid), along y + //! Order of finite centering of currents (from nodal grid to staggered + //! grid), along y static int current_centering_noy; - //! Order of finite centering of currents (from nodal grid to staggered grid), along z + //! Order of finite centering of currents (from nodal grid to staggered + //! grid), along z static int current_centering_noz; //! Number of modes for the RZ multi-mode version static int n_rz_azimuthal_modes; //! Number of MultiFab components - //! (with the RZ multi-mode version, each mode has a real and imaginary component, - //! except mode 0 which is purely real: component 0 is mode 0, odd components are - //! the real parts, even components are the imaginary parts) + //! (with the RZ multi-mode version, each mode has a real and imaginary + //! component, except mode 0 which is purely real: component 0 is mode 0, + //! odd components are the real parts, even components are the imaginary + //! parts) static int ncomps; //! If true, a Numerical Cherenkov Instability (NCI) corrector is applied //! (for simulations using the FDTD Maxwell solver) static bool use_fdtd_nci_corr; - //! If true (Galerkin method): The fields are interpolated directly from the staggered - //! grid to the particle positions (with reduced interpolation order in the parallel - //! direction). This scheme is energy conserving in the limit of infinitely small time - //! steps. - //! Otherwise, "momentum conserving" (in the same limit): The fields are first - //! interpolated from the staggered grid points to the corners of each cell, and then - //! from the cell corners to the particle position (with the same order of interpolation - //! in all directions). This scheme is momentum conserving in the limit of infinitely - //! small time steps. + //! If true (Galerkin method): The fields are interpolated directly from the + //! staggered grid to the particle positions (with reduced interpolation + //! order in the parallel direction). This scheme is energy conserving in + //! the limit of infinitely small time steps. Otherwise, "momentum + //! conserving" (in the same limit): The fields are first interpolated from + //! the staggered grid points to the corners of each cell, and then from the + //! cell corners to the particle position (with the same order of + //! interpolation in all directions). This scheme is momentum conserving in + //! the limit of infinitely small time steps. static bool galerkin_interpolation; //! If true, a bilinear filter is used to smooth charge and currents static bool use_filter; - //! If true, the bilinear filtering of charge and currents is done in Fourier space + //! If true, the bilinear filtering of charge and currents is done in + //! Fourier space static bool use_kspace_filter; - //! If true, a compensation step is added to the bilinear filtering of charge and currents + //! If true, a compensation step is added to the bilinear filtering of + //! charge and currents static bool use_filter_compensation; - //! If true, the initial conditions from random number generators are serialized (useful for reproducible testing with OpenMP) + //! If true, the initial conditions from random number generators are + //! serialized (useful for reproducible testing with OpenMP) static bool serialize_initial_conditions; - //! Lorentz factor of the boosted frame in which a boosted-frame simulation is run + //! Lorentz factor of the boosted frame in which a boosted-frame simulation + //! is run static amrex::Real gamma_boost; - //! Beta value corresponding to the Lorentz factor of the boosted frame of the simulation + //! Beta value corresponding to the Lorentz factor of the boosted frame of + //! the simulation static amrex::Real beta_boost; - //! Direction of the Lorentz transform that defines the boosted frame of the simulation + //! Direction of the Lorentz transform that defines the boosted frame of the + //! simulation static amrex::Vector boost_direction; - //! store initial value of zmin_domain_boost because WarpX::computeMaxStepBoostAccelerator - //! needs the initial value of zmin_domain_boost, even if restarting from a checkpoint file + //! store initial value of zmin_domain_boost because + //! WarpX::computeMaxStepBoostAccelerator needs the initial value of + //! zmin_domain_boost, even if restarting from a checkpoint file static amrex::Real zmin_domain_boost_step_0; - //! If true, the code will compute max_step from the back transformed diagnostics + //! If true, the code will compute max_step from the back transformed + //! diagnostics static bool compute_max_step_from_btd; static bool do_dynamic_scheduling; @@ -352,35 +419,40 @@ public: static utils::parser::IntervalsParser sort_intervals; static amrex::IntVect sort_bin_size; - //! If true, particles will be sorted in the order x -> y -> z -> ppc for faster deposition + //! If true, particles will be sorted in the order x -> y -> z -> ppc for + //! faster deposition static bool sort_particles_for_deposition; - //! Specifies the type of grid used for the above sorting, i.e. cell-centered, nodal, or mixed + //! Specifies the type of grid used for the above sorting, i.e. + //! cell-centered, nodal, or mixed static amrex::IntVect sort_idx_type; static bool do_multi_J; static int do_multi_J_n_depositions; - //! With mesh refinement, particles located inside a refinement patch, but within - //! #n_field_gather_buffer cells of the edge of the patch, will gather the fields - //! from the lower refinement level instead of the refinement patch itself + //! With mesh refinement, particles located inside a refinement patch, but + //! within #n_field_gather_buffer cells of the edge of the patch, will + //! gather the fields from the lower refinement level instead of the + //! refinement patch itself static int n_field_gather_buffer; - //! With mesh refinement, particles located inside a refinement patch, but within - //! #n_current_deposition_buffer cells of the edge of the patch, will deposit their charge - //! and current onto the lower refinement level instead of the refinement patch itself + //! With mesh refinement, particles located inside a refinement patch, but + //! within #n_current_deposition_buffer cells of the edge of the patch, will + //! deposit their charge and current onto the lower refinement level instead + //! of the refinement patch itself static int n_current_deposition_buffer; //! Integer that corresponds to the type of grid used in the simulation //! (collocated, staggered, hybrid) static inline auto grid_type = ablastr::utils::enums::GridType::Default; - // Global rho nodal flag to know about rho index type when rho MultiFab is not allocated + // Global rho nodal flag to know about rho index type when rho MultiFab is + // not allocated amrex::IntVect m_rho_nodal_flag; /** * \brief - * Allocate and optionally initialize the iMultiFab. This also adds the iMultiFab - * to the map of MultiFabs (used to ease the access to MultiFabs from the Python - * interface + * Allocate and optionally initialize the iMultiFab. This also adds the + * iMultiFab to the map of MultiFabs (used to ease the access to MultiFabs + * from the Python interface * * \param[out] mf The iMultiFab unique pointer to be allocated * \param[in] ba The BoxArray describing the iMultiFab @@ -391,37 +463,46 @@ public: * \param[in] name The name of the iMultiFab to use in the map * \param[in] initial_value The optional initial value */ - void AllocInitMultiFab ( - std::unique_ptr& mf, - const amrex::BoxArray& ba, - const amrex::DistributionMapping& dm, - int ncomp, - const amrex::IntVect& ngrow, - int level, - const std::string& name, - std::optional initial_value = {}); + void AllocInitMultiFab (std::unique_ptr& mf, + const amrex::BoxArray& ba, + const amrex::DistributionMapping& dm, int ncomp, + const amrex::IntVect& ngrow, int level, + const std::string& name, + std::optional initial_value = {}); - // Maps of all of the iMultiFabs used (this can include MFs from other classes) - // This is a convenience for the Python interface, allowing all iMultiFabs - // to be easily referenced from Python. - std::map imultifab_map; + // Maps of all of the iMultiFabs used (this can include MFs from other + // classes) This is a convenience for the Python interface, allowing all + // iMultiFabs to be easily referenced from Python. + std::map imultifab_map; /** * \brief - * Get pointer to the amrex::MultiFab containing the dotMask for the specified field + * Get pointer to the amrex::MultiFab containing the dotMask for the + * specified field */ [[nodiscard]] const amrex::iMultiFab* - getFieldDotMaskPointer (warpx::fields::FieldType field_type, int lev, ablastr::fields::Direction dir) const; + getFieldDotMaskPointer (warpx::fields::FieldType field_type, int lev, + ablastr::fields::Direction dir) const; - [[nodiscard]] bool DoPML () const {return do_pml;} - [[nodiscard]] bool DoFluidSpecies () const {return do_fluid_species;} + [[nodiscard]] bool + DoPML () const { + return do_pml; + } + [[nodiscard]] bool + DoFluidSpecies () const { + return do_fluid_species; + } #if (defined WARPX_DIM_RZ) && (defined WARPX_USE_FFT) - const PML_RZ* getPMLRZ() {return pml_rz[0].get();} + const PML_RZ* + getPMLRZ () { + return pml_rz[0].get(); + } #endif - /** get low-high-low-high-... vector for each direction indicating if mother grid PMLs are enabled */ - [[nodiscard]] std::vector getPMLdirections() const; + /** get low-high-low-high-... vector for each direction indicating if mother + * grid PMLs are enabled */ + [[nodiscard]] std::vector getPMLdirections () const; static amrex::LayoutData* getCosts (int lev); @@ -431,33 +512,38 @@ public: static amrex::IntVect filter_npass_each_dir; BilinearFilter bilinear_filter; - amrex::Vector< std::unique_ptr > nci_godfrey_filter_exeybz; - amrex::Vector< std::unique_ptr > nci_godfrey_filter_bxbyez; + amrex::Vector> nci_godfrey_filter_exeybz; + amrex::Vector> nci_godfrey_filter_bxbyez; amrex::Real time_of_last_gal_shift = 0; - amrex::Vector m_v_galilean = amrex::Vector(3, amrex::Real(0.)); - amrex::Array m_galilean_shift = {{0}}; + amrex::Vector m_v_galilean = + amrex::Vector(3, amrex::Real(0.)); + amrex::Array m_galilean_shift = {{0}}; - amrex::Vector m_v_comoving = amrex::Vector(3, amrex::Real(0.)); + amrex::Vector m_v_comoving = + amrex::Vector(3, amrex::Real(0.)); - /// object with all reduced diagnostics, similar to MultiParticleContainer for species. + /// object with all reduced diagnostics, similar to MultiParticleContainer + /// for species. std::unique_ptr reduced_diags; - void applyMirrors(amrex::Real time); + void applyMirrors (amrex::Real time); /** Determine the timestep of the simulation. */ void ComputeDt (); /** * Determine the simulation timestep from the maximum speed of all particles - * Sets timestep so that a particle can only cross cfl*dx cells per timestep. + * Sets timestep so that a particle can only cross cfl*dx cells per + * timestep. */ void UpdateDtFromParticleSpeeds (); /** Print main PIC parameters to stdout */ void PrintMainPICparameters (); - /** Write a file that record all inputs: inputs file + command line options */ + /** Write a file that record all inputs: inputs file + command line options + */ void WriteUsedInputsFile () const; /** @@ -467,18 +553,20 @@ public: */ void ComputeMaxStep (); // Compute max_step automatically for simulations in a boosted frame. - void computeMaxStepBoostAccelerator(); + void computeMaxStepBoostAccelerator (); /** \brief Move the moving window * \param step Time step - * \param move_j whether the current (and the charge, if allocated) is shifted or not + * \param move_j whether the current (and the charge, if allocated) is + * shifted or not */ - int MoveWindow (int step, bool move_j); + int MoveWindow (int step, bool move_j); /** * \brief * This function shifts the boundary of the grid by 'm_v_galilean*dt'. - * In doding so, only positions attributes are changed while fields remain unchanged. + * In doding so, only positions attributes are changed while fields remain + * unchanged. */ void ShiftGalileanBoundary (); @@ -489,22 +577,28 @@ public: void UpdateInjectionPosition (amrex::Real dt); void ResetProbDomain (const amrex::RealBox& rb); - void EvolveE ( amrex::Real dt, amrex::Real start_time); + void EvolveE (amrex::Real dt, amrex::Real start_time); void EvolveE (int lev, amrex::Real dt, amrex::Real start_time); - void EvolveB ( amrex::Real dt, DtType dt_type, amrex::Real start_time); - void EvolveB (int lev, amrex::Real dt, DtType dt_type, amrex::Real start_time); - void EvolveF ( amrex::Real dt, DtType dt_type); + void EvolveB (amrex::Real dt, DtType dt_type, amrex::Real start_time); + void EvolveB (int lev, amrex::Real dt, DtType dt_type, + amrex::Real start_time); + void EvolveF (amrex::Real dt, DtType dt_type); void EvolveF (int lev, amrex::Real dt, DtType dt_type); - void EvolveG ( amrex::Real dt, DtType dt_type); + void EvolveG (amrex::Real dt, DtType dt_type); void EvolveG (int lev, amrex::Real dt, DtType dt_type); - void EvolveB (int lev, PatchType patch_type, amrex::Real dt, DtType dt_type, amrex::Real start_time); - void EvolveE (int lev, PatchType patch_type, amrex::Real dt, amrex::Real start_time); - void EvolveF (int lev, PatchType patch_type, amrex::Real dt, DtType dt_type); - void EvolveG (int lev, PatchType patch_type, amrex::Real dt, DtType dt_type); - - void MacroscopicEvolveE ( amrex::Real dt, amrex::Real start_time); + void EvolveB (int lev, PatchType patch_type, amrex::Real dt, DtType dt_type, + amrex::Real start_time); + void EvolveE (int lev, PatchType patch_type, amrex::Real dt, + amrex::Real start_time); + void EvolveF (int lev, PatchType patch_type, amrex::Real dt, + DtType dt_type); + void EvolveG (int lev, PatchType patch_type, amrex::Real dt, + DtType dt_type); + + void MacroscopicEvolveE (amrex::Real dt, amrex::Real start_time); void MacroscopicEvolveE (int lev, amrex::Real dt, amrex::Real start_time); - void MacroscopicEvolveE (int lev, PatchType patch_type, amrex::Real dt, amrex::Real start_time); + void MacroscopicEvolveE (int lev, PatchType patch_type, amrex::Real dt, + amrex::Real start_time); /** * \brief Hybrid-PIC field evolve function. @@ -526,7 +620,7 @@ public: * * \param dt vector of time steps (for all levels) */ - void Hybrid_QED_Push ( amrex::Vector dt); + void Hybrid_QED_Push (amrex::Vector dt); /** apply QED correction on electric field for level lev * @@ -535,7 +629,8 @@ public: */ void Hybrid_QED_Push (int lev, amrex::Real dt); - /** apply QED correction on electric field for level lev and patch type patch_type + /** apply QED correction on electric field for level lev and patch type + * patch_type * * \param lev mesh refinement level * \param patch_type which MR patch: PatchType::fine or PatchType::coarse @@ -549,7 +644,8 @@ public: */ void CheckLoadBalance (int step); - /** \brief perform load balance; compute and communicate new `amrex::DistributionMapping` + /** \brief perform load balance; compute and communicate new + * `amrex::DistributionMapping` */ void LoadBalance (); @@ -567,8 +663,8 @@ public: /** \brief returns the load balance interval */ - [[nodiscard]] utils::parser::IntervalsParser get_load_balance_intervals () const - { + [[nodiscard]] utils::parser::IntervalsParser + get_load_balance_intervals () const { return load_balance_intervals; } @@ -593,13 +689,13 @@ public: void DampFieldsInGuards (int lev, amrex::MultiFab* mf); #ifdef WARPX_DIM_RZ - void ApplyInverseVolumeScalingToCurrentDensity(amrex::MultiFab* Jx, - amrex::MultiFab* Jy, - amrex::MultiFab* Jz, - int lev) const; + void ApplyInverseVolumeScalingToCurrentDensity (amrex::MultiFab* Jx, + amrex::MultiFab* Jy, + amrex::MultiFab* Jz, + int lev) const; - void ApplyInverseVolumeScalingToChargeDensity(amrex::MultiFab* Rho, - int lev) const; + void ApplyInverseVolumeScalingToChargeDensity (amrex::MultiFab* Rho, + int lev) const; #endif /** @@ -615,16 +711,18 @@ public: * density deposited in the guard cells have to be reflected back into * the simulation domain. This function performs that logic. */ - void ApplyJfieldBoundary (int lev, amrex::MultiFab* Jx, - amrex::MultiFab* Jy, amrex::MultiFab* Jz, - PatchType patch_type); + void ApplyJfieldBoundary (int lev, amrex::MultiFab* Jx, amrex::MultiFab* Jy, + amrex::MultiFab* Jz, PatchType patch_type); - void ApplyEfieldBoundary (int lev, PatchType patch_type, amrex::Real cur_time); - void ApplyBfieldBoundary (int lev, PatchType patch_type, DtType dt_type, amrex::Real cur_time); + void ApplyEfieldBoundary (int lev, PatchType patch_type, + amrex::Real cur_time); + void ApplyBfieldBoundary (int lev, PatchType patch_type, DtType dt_type, + amrex::Real cur_time); #ifdef WARPX_DIM_RZ // Applies the boundary conditions that are specific to the axis when in RZ. - void ApplyFieldBoundaryOnAxis (amrex::MultiFab* Er, amrex::MultiFab* Et, amrex::MultiFab* Ez, int lev) const; + void ApplyFieldBoundaryOnAxis (amrex::MultiFab* Er, amrex::MultiFab* Et, + amrex::MultiFab* Ez, int lev) const; #endif /** @@ -649,7 +747,7 @@ public: void CopyJPML (); /** True if any of the particle boundary condition type is Thermal */ - static bool isAnyParticleBoundaryThermal(); + static bool isAnyParticleBoundaryThermal (); PML* GetPML (int lev); #if (defined WARPX_DIM_RZ) && (defined WARPX_USE_FFT) @@ -664,10 +762,13 @@ public: void doQEDEvents (); #endif - void PushParticlesandDeposit (int lev, amrex::Real cur_time, DtType a_dt_type=DtType::Full, bool skip_current=false, - PushType push_type=PushType::Explicit); - void PushParticlesandDeposit (amrex::Real cur_time, bool skip_current=false, - PushType push_type=PushType::Explicit); + void PushParticlesandDeposit (int lev, amrex::Real cur_time, + DtType a_dt_type = DtType::Full, + bool skip_current = false, + PushType push_type = PushType::Explicit); + void PushParticlesandDeposit (amrex::Real cur_time, + bool skip_current = false, + PushType push_type = PushType::Explicit); // This function does aux(lev) = fp(lev) + I(aux(lev-1)-cp(lev)). // Caller must make sure fp and cp have ghost cells filled. @@ -680,27 +781,38 @@ public: * it centers the currents from a nodal grid to a staggered grid (Yee) using * finite-order interpolation based on the Fornberg coefficients. * - * \param[in,out] dst destination \c MultiFab where the results of the finite-order centering are stored - * \param[in] src source \c MultiFab that contains the values of the nodal current to be centered + * \param[in,out] dst destination \c MultiFab where the results of the + * finite-order centering are stored + * \param[in] src source \c MultiFab that contains the values of the nodal + * current to be centered */ - void UpdateCurrentNodalToStag (amrex::MultiFab& dst, amrex::MultiFab const& src); + void UpdateCurrentNodalToStag (amrex::MultiFab& dst, + amrex::MultiFab const& src); // Fill boundary cells including coarse/fine boundaries - void FillBoundaryB (amrex::IntVect ng, std::optional nodal_sync = std::nullopt); - void FillBoundaryE (amrex::IntVect ng, std::optional nodal_sync = std::nullopt); - void FillBoundaryB_avg (amrex::IntVect ng); - void FillBoundaryE_avg (amrex::IntVect ng); - - void FillBoundaryF (amrex::IntVect ng, std::optional nodal_sync = std::nullopt); - void FillBoundaryG (amrex::IntVect ng, std::optional nodal_sync = std::nullopt); + void FillBoundaryB (amrex::IntVect ng, + std::optional nodal_sync = std::nullopt); + void FillBoundaryE (amrex::IntVect ng, + std::optional nodal_sync = std::nullopt); + void FillBoundaryB_avg (amrex::IntVect ng); + void FillBoundaryE_avg (amrex::IntVect ng); + + void FillBoundaryF (amrex::IntVect ng, + std::optional nodal_sync = std::nullopt); + void FillBoundaryG (amrex::IntVect ng, + std::optional nodal_sync = std::nullopt); void FillBoundaryAux (amrex::IntVect ng); - void FillBoundaryE (int lev, amrex::IntVect ng, std::optional nodal_sync = std::nullopt); - void FillBoundaryB (int lev, amrex::IntVect ng, std::optional nodal_sync = std::nullopt); - void FillBoundaryE_avg (int lev, amrex::IntVect ng); - void FillBoundaryB_avg (int lev, amrex::IntVect ng); - - void FillBoundaryF (int lev, amrex::IntVect ng, std::optional nodal_sync = std::nullopt); - void FillBoundaryG (int lev, amrex::IntVect ng, std::optional nodal_sync = std::nullopt); + void FillBoundaryE (int lev, amrex::IntVect ng, + std::optional nodal_sync = std::nullopt); + void FillBoundaryB (int lev, amrex::IntVect ng, + std::optional nodal_sync = std::nullopt); + void FillBoundaryE_avg (int lev, amrex::IntVect ng); + void FillBoundaryB_avg (int lev, amrex::IntVect ng); + + void FillBoundaryF (int lev, amrex::IntVect ng, + std::optional nodal_sync = std::nullopt); + void FillBoundaryG (int lev, amrex::IntVect ng, + std::optional nodal_sync = std::nullopt); void FillBoundaryAux (int lev, amrex::IntVect ng); /** @@ -725,73 +837,134 @@ public: void SyncRho (); - void SyncRho ( - const ablastr::fields::MultiLevelScalarField& charge_fp, - const ablastr::fields::MultiLevelScalarField& charge_cp, - ablastr::fields::MultiLevelScalarField const & charge_buffer); - - [[nodiscard]] amrex::Vector getnsubsteps () const {return nsubsteps;} - [[nodiscard]] int getnsubsteps (int lev) const {return nsubsteps[lev];} - [[nodiscard]] amrex::Vector getistep () const {return istep;} - [[nodiscard]] int getistep (int lev) const {return istep[lev];} - void setistep (int lev, int ii) {istep[lev] = ii;} - [[nodiscard]] amrex::Vector gett_old () const {return t_old;} - [[nodiscard]] amrex::Real gett_old (int lev) const {return t_old[lev];} - [[nodiscard]] amrex::Vector gett_new () const {return t_new;} - [[nodiscard]] amrex::Real gett_new (int lev) const {return t_new[lev];} - void sett_new (int lev, amrex::Real time) {t_new[lev] = time;} - [[nodiscard]] amrex::Vector getdt () const {return dt;} - [[nodiscard]] amrex::Real getdt (int lev) const {return dt.at(lev);} - [[nodiscard]] int getdo_moving_window() const {return do_moving_window;} - [[nodiscard]] amrex::Real getmoving_window_x() const {return moving_window_x;} - [[nodiscard]] bool getis_synchronized() const {return is_synchronized;} - - [[nodiscard]] int maxStep () const {return max_step;} - void updateMaxStep (const int new_max_step) {max_step = new_max_step;} - [[nodiscard]] amrex::Real stopTime () const {return stop_time;} - void updateStopTime (const amrex::Real new_stop_time) {stop_time = new_stop_time;} - - static std::array CellSize (int lev); + void SyncRho (const ablastr::fields::MultiLevelScalarField& charge_fp, + const ablastr::fields::MultiLevelScalarField& charge_cp, + ablastr::fields::MultiLevelScalarField const& charge_buffer); + + [[nodiscard]] amrex::Vector + getnsubsteps () const { + return nsubsteps; + } + [[nodiscard]] int + getnsubsteps (int lev) const { + return nsubsteps[lev]; + } + [[nodiscard]] amrex::Vector + getistep () const { + return istep; + } + [[nodiscard]] int + getistep (int lev) const { + return istep[lev]; + } + void + setistep (int lev, int ii) { + istep[lev] = ii; + } + [[nodiscard]] amrex::Vector + gett_old () const { + return t_old; + } + [[nodiscard]] amrex::Real + gett_old (int lev) const { + return t_old[lev]; + } + [[nodiscard]] amrex::Vector + gett_new () const { + return t_new; + } + [[nodiscard]] amrex::Real + gett_new (int lev) const { + return t_new[lev]; + } + void + sett_new (int lev, amrex::Real time) { + t_new[lev] = time; + } + [[nodiscard]] amrex::Vector + getdt () const { + return dt; + } + [[nodiscard]] amrex::Real + getdt (int lev) const { + return dt.at(lev); + } + [[nodiscard]] int + getdo_moving_window () const { + return do_moving_window; + } + [[nodiscard]] amrex::Real + getmoving_window_x () const { + return moving_window_x; + } + [[nodiscard]] bool + getis_synchronized () const { + return is_synchronized; + } + + [[nodiscard]] int + maxStep () const { + return max_step; + } + void + updateMaxStep (const int new_max_step) { + max_step = new_max_step; + } + [[nodiscard]] amrex::Real + stopTime () const { + return stop_time; + } + void + updateStopTime (const amrex::Real new_stop_time) { + stop_time = new_stop_time; + } + + static std::array CellSize (int lev); static amrex::XDim3 InvCellSize (int lev); - static amrex::RealBox getRealBox(const amrex::Box& bx, int lev); + static amrex::RealBox getRealBox (const amrex::Box& bx, int lev); /** * \brief Return the lower corner of the box in real units. * \param bx The box * \param lev The refinement level of the box - * \param time_shift_delta The time relative to the current time at which to calculate the position - * (when v_galilean is not zero) + * \param time_shift_delta The time relative to the current time at which to + * calculate the position (when v_galilean is not zero) * \return An array of the position coordinates */ - static amrex::XDim3 LowerCorner (const amrex::Box& bx, int lev, amrex::Real time_shift_delta); + static amrex::XDim3 LowerCorner (const amrex::Box& bx, int lev, + amrex::Real time_shift_delta); /** * \brief Return the upper corner of the box in real units. * \param bx The box * \param lev The refinement level of the box - * \param time_shift_delta The time relative to the current time at which to calculate the position - * (when v_galilean is not zero) + * \param time_shift_delta The time relative to the current time at which to + * calculate the position (when v_galilean is not zero) * \return An array of the position coordinates */ - static amrex::XDim3 UpperCorner (const amrex::Box& bx, int lev, amrex::Real time_shift_delta); + static amrex::XDim3 UpperCorner (const amrex::Box& bx, int lev, + amrex::Real time_shift_delta); static amrex::IntVect RefRatio (int lev); static const amrex::iMultiFab* CurrentBufferMasks (int lev); static const amrex::iMultiFab* GatherBufferMasks (int lev); - static inline auto electrostatic_solver_id = ElectrostaticSolverAlgo::Default; + static inline auto electrostatic_solver_id = + ElectrostaticSolverAlgo::Default; static inline auto poisson_solver_id = PoissonSolverAlgo::Default; - static int do_moving_window; // boolean + static int do_moving_window; // boolean static int start_moving_window_step; // the first step to move window - static int end_moving_window_step; // the last step to move window + static int end_moving_window_step; // the last step to move window /** Returns true if the moving window is active for the provided step * * @param step time step * @return true if active, else false */ - static int moving_window_active (int const step) { - bool const step_before_end = (step < end_moving_window_step) || (end_moving_window_step < 0); + static int + moving_window_active (int const step) { + bool const step_before_end = + (step < end_moving_window_step) || (end_moving_window_step < 0); bool const step_after_start = (step >= start_moving_window_step); return do_moving_window && step_before_end && step_after_start; } @@ -801,51 +974,76 @@ public: // these should be private, but can't due to Cuda limitations static void ComputeDivB (amrex::MultiFab& divB, int dcomp, - ablastr::fields::VectorField const & B, - const std::array& dx); + ablastr::fields::VectorField const& B, + const std::array& dx); static void ComputeDivB (amrex::MultiFab& divB, int dcomp, - ablastr::fields::VectorField const & B, - const std::array& dx, amrex::IntVect ngrow); + ablastr::fields::VectorField const& B, + const std::array& dx, + amrex::IntVect ngrow); - void ComputeDivE(amrex::MultiFab& divE, int lev); + void ComputeDivE (amrex::MultiFab& divE, int lev); void ProjectionCleanDivB (); void CalculateExternalCurlA (); - [[nodiscard]] amrex::IntVect getngEB() const { return guard_cells.ng_alloc_EB; } - [[nodiscard]] amrex::IntVect getngF() const { return guard_cells.ng_alloc_F; } - [[nodiscard]] amrex::IntVect getngUpdateAux() const { return guard_cells.ng_UpdateAux; } - [[nodiscard]] amrex::IntVect get_ng_depos_J() const {return guard_cells.ng_depos_J;} - [[nodiscard]] amrex::IntVect get_ng_depos_rho() const {return guard_cells.ng_depos_rho;} - [[nodiscard]] amrex::IntVect get_ng_fieldgather () const {return guard_cells.ng_FieldGather;} + [[nodiscard]] amrex::IntVect + getngEB () const { + return guard_cells.ng_alloc_EB; + } + [[nodiscard]] amrex::IntVect + getngF () const { + return guard_cells.ng_alloc_F; + } + [[nodiscard]] amrex::IntVect + getngUpdateAux () const { + return guard_cells.ng_UpdateAux; + } + [[nodiscard]] amrex::IntVect + get_ng_depos_J () const { + return guard_cells.ng_depos_J; + } + [[nodiscard]] amrex::IntVect + get_ng_depos_rho () const { + return guard_cells.ng_depos_rho; + } + [[nodiscard]] amrex::IntVect + get_ng_fieldgather () const { + return guard_cells.ng_FieldGather; + } /** Coarsest-level Domain Decomposition * * If specified, the domain will be chopped into the exact number * of pieces in each dimension as specified by this parameter. * - * @return the number of MPI processes per dimension if specified, otherwise a 0-vector + * @return the number of MPI processes per dimension if specified, otherwise + * a 0-vector */ - [[nodiscard]] amrex::IntVect get_numprocs() const {return numprocs;} + [[nodiscard]] amrex::IntVect + get_numprocs () const { + return numprocs; + } /** Electrostatic solve call */ void ComputeSpaceChargeField (bool reset_fields); // Magnetostatic Solver Interface - MagnetostaticSolver::VectorPoissonBoundaryHandler m_vector_poisson_boundary_handler; + MagnetostaticSolver::VectorPoissonBoundaryHandler + m_vector_poisson_boundary_handler; int magnetostatic_solver_max_iters = 200; int magnetostatic_solver_verbosity = 2; void ComputeMagnetostaticField (); void AddMagnetostaticFieldLabFrame (); - void computeVectorPotential (ablastr::fields::MultiLevelVectorField const& curr, - ablastr::fields::MultiLevelVectorField const& A, - amrex::Real required_precision=amrex::Real(1.e-11), - amrex::Real absolute_tolerance=amrex::Real(0.0), - int max_iters=200, - int verbosity=2); // const; + void computeVectorPotential ( + ablastr::fields::MultiLevelVectorField const& curr, + ablastr::fields::MultiLevelVectorField const& A, + amrex::Real required_precision = amrex::Real(1.e-11), + amrex::Real absolute_tolerance = amrex::Real(0.0), int max_iters = 200, + int verbosity = 2); // const; - void setVectorPotentialBC (ablastr::fields::MultiLevelVectorField const& A) const; + void setVectorPotentialBC ( + ablastr::fields::MultiLevelVectorField const& A) const; /** * \brief @@ -861,28 +1059,33 @@ public: * \param[in] fy_parser parser function to initialize y-field * \param[in] fz_parser parser function to initialize z-field * \param[in] lev level of the Multifabs that is initialized - * \param[in] patch_type PatchType on which the field is initialized (fine or coarse) - * \param[in] eb_update_field flag indicating which gridpoints should be modified by this functions - * \param[in] use_eb_flags (default:true) flag indicating if eb points should be excluded or not + * \param[in] patch_type PatchType on which the field is initialized (fine + * or coarse) + * \param[in] eb_update_field flag indicating which gridpoints should be + * modified by this functions + * \param[in] use_eb_flags (default:true) flag indicating if eb points + * should be excluded or not */ void ComputeExternalFieldOnGridUsingParser ( warpx::fields::FieldType field, amrex::ParserExecutor<4> const& fx_parser, amrex::ParserExecutor<4> const& fy_parser, - amrex::ParserExecutor<4> const& fz_parser, - int lev, PatchType patch_type, - amrex::Vector,3 > > const& eb_update_field, + amrex::ParserExecutor<4> const& fz_parser, int lev, + PatchType patch_type, + amrex::Vector, 3>> const& + eb_update_field, bool use_eb_flags); void ComputeExternalFieldOnGridUsingParser ( warpx::fields::FieldType field, amrex::ParserExecutor<4> const& fx_parser, amrex::ParserExecutor<4> const& fy_parser, - amrex::ParserExecutor<4> const& fz_parser, - int lev, PatchType patch_type, - amrex::Vector,3 > > const& eb_update_field); + amrex::ParserExecutor<4> const& fz_parser, int lev, + PatchType patch_type, + amrex::Vector, 3>> const& + eb_update_field); - /** + /** * \brief * This function computes the E, B, and J fields on each level * using the parser and the user-defined function for the external fields. @@ -899,26 +1102,29 @@ public: * \param[in] face_areas face areas information * \param[in] topology flag indicating if field is edge-based or face-based * \param[in] lev level of the Multifabs that is initialized - * \param[in] patch_type PatchType on which the field is initialized (fine or coarse) - * \param[in] eb_update_field flag indicating which gridpoints should be modified by this functions - * \param[in] use_eb_flags (default:true) flag indicating if eb points should be excluded or not + * \param[in] patch_type PatchType on which the field is initialized (fine + * or coarse) + * \param[in] eb_update_field flag indicating which gridpoints should be + * modified by this functions + * \param[in] use_eb_flags (default:true) flag indicating if eb points + * should be excluded or not */ void ComputeExternalFieldOnGridUsingParser ( - std::string const& field, - amrex::ParserExecutor<4> const& fx_parser, + std::string const& field, amrex::ParserExecutor<4> const& fx_parser, amrex::ParserExecutor<4> const& fy_parser, - amrex::ParserExecutor<4> const& fz_parser, - int lev, PatchType patch_type, - amrex::Vector< std::array< std::unique_ptr,3> > const& eb_update_field, + amrex::ParserExecutor<4> const& fz_parser, int lev, + PatchType patch_type, + amrex::Vector, 3>> const& + eb_update_field, bool use_eb_flags); void ComputeExternalFieldOnGridUsingParser ( - std::string const& field, - amrex::ParserExecutor<4> const& fx_parser, + std::string const& field, amrex::ParserExecutor<4> const& fx_parser, amrex::ParserExecutor<4> const& fy_parser, - amrex::ParserExecutor<4> const& fz_parser, - int lev, PatchType patch_type, - amrex::Vector< std::array< std::unique_ptr,3> > const& eb_update_field); + amrex::ParserExecutor<4> const& fz_parser, int lev, + PatchType patch_type, + amrex::Vector, 3>> const& + eb_update_field); /** * \brief Load field values from a user-specified openPMD file, @@ -930,9 +1136,10 @@ public: * \brief Load field values from a user-specified openPMD file * for a specific field (specified by `F_name`) */ - void ReadExternalFieldFromFile ( - const std::string& read_fields_from_path, amrex::MultiFab* mf, - const std::string& F_name, const std::string& F_component); + void ReadExternalFieldFromFile (const std::string& read_fields_from_path, + amrex::MultiFab* mf, + const std::string& F_name, + const std::string& F_component); /** * \brief @@ -942,115 +1149,143 @@ public: * * \param[in] lev level of the Multifabs that is initialized */ - void InitializeEBGridData(int lev); - - /** \brief adds particle and cell contributions in cells to compute heuristic - * cost in each box on each level, and records in `costs` - * @param[in] costs vector of (`unique_ptr` to) vectors; expected to be initialized - * to correct number of boxes and boxes per level - */ - void ComputeCostsHeuristic (amrex::Vector > >& costs); - - void ApplyFilterandSumBoundaryRho (int lev, int glev, amrex::MultiFab& rho, int icomp, int ncomp); - - // Device vectors of stencil coefficients used for finite-order centering of fields - amrex::Gpu::DeviceVector device_field_centering_stencil_coeffs_x; - amrex::Gpu::DeviceVector device_field_centering_stencil_coeffs_y; - amrex::Gpu::DeviceVector device_field_centering_stencil_coeffs_z; - - // Device vectors of stencil coefficients used for finite-order centering of currents - amrex::Gpu::DeviceVector device_current_centering_stencil_coeffs_x; - amrex::Gpu::DeviceVector device_current_centering_stencil_coeffs_y; - amrex::Gpu::DeviceVector device_current_centering_stencil_coeffs_z; + void InitializeEBGridData (int lev); + + /** \brief adds particle and cell contributions in cells to compute + * heuristic cost in each box on each level, and records in `costs` + * @param[in] costs vector of (`unique_ptr` to) vectors; expected to be + * initialized to correct number of boxes and boxes per level + */ + void ComputeCostsHeuristic ( + amrex::Vector>>& costs); + + void ApplyFilterandSumBoundaryRho (int lev, int glev, amrex::MultiFab& rho, + int icomp, int ncomp); + + // Device vectors of stencil coefficients used for finite-order centering of + // fields + amrex::Gpu::DeviceVector + device_field_centering_stencil_coeffs_x; + amrex::Gpu::DeviceVector + device_field_centering_stencil_coeffs_y; + amrex::Gpu::DeviceVector + device_field_centering_stencil_coeffs_z; + + // Device vectors of stencil coefficients used for finite-order centering of + // currents + amrex::Gpu::DeviceVector + device_current_centering_stencil_coeffs_x; + amrex::Gpu::DeviceVector + device_current_centering_stencil_coeffs_y; + amrex::Gpu::DeviceVector + device_current_centering_stencil_coeffs_z; // This needs to be public for CUDA. //! Tagging cells for refinement - void ErrorEst (int lev, amrex::TagBoxArray& tags, amrex::Real time, int /*ngrow*/) final; - - // Return the accelerator lattice instance defined at the given refinement level - const AcceleratorLattice& get_accelerator_lattice (int lev) {return *(m_accelerator_lattice[lev]);} + void ErrorEst (int lev, amrex::TagBoxArray& tags, amrex::Real time, + int /*ngrow*/) final; + + // Return the accelerator lattice instance defined at the given refinement + // level + const AcceleratorLattice& + get_accelerator_lattice (int lev) { + return *(m_accelerator_lattice[lev]); + } // for cuda - void BuildBufferMasksInBox ( amrex::Box tbx, amrex::IArrayBox &buffer_mask, - const amrex::IArrayBox &guard_mask, int ng ); + void BuildBufferMasksInBox (amrex::Box tbx, amrex::IArrayBox& buffer_mask, + const amrex::IArrayBox& guard_mask, int ng); #ifdef AMREX_USE_EB - [[nodiscard]] amrex::EBFArrayBoxFactory const& fieldEBFactory (int lev) const noexcept { - return static_cast(*m_field_factory[lev]); + [[nodiscard]] amrex::EBFArrayBoxFactory const& + fieldEBFactory (int lev) const noexcept { + return static_cast( + *m_field_factory[lev]); } #endif void InitEB (); /** - * \brief Compute the level set function used for particle-boundary interaction. - */ + * \brief Compute the level set function used for particle-boundary + * interaction. + */ void ComputeDistanceToEB (); /** - * \brief Auxiliary function to count the amount of faces which still need to be extended - */ - amrex::Array1D CountExtFaces(); + * \brief Auxiliary function to count the amount of faces which still need + * to be extended + */ + amrex::Array1D CountExtFaces (); /** - * \brief Main function computing the cell extension. Where possible it computes one-way - * extensions and, when this is not possible, it does eight-ways extensions. - */ - void ComputeFaceExtensions(); + * \brief Main function computing the cell extension. Where possible it + * computes one-way extensions and, when this is not possible, it does + * eight-ways extensions. + */ + void ComputeFaceExtensions (); /** - * \brief Initialize the memory for the FaceInfoBoxes - */ - void InitBorrowing(); + * \brief Initialize the memory for the FaceInfoBoxes + */ + void InitBorrowing (); /** - * \brief Shrink the vectors in the FaceInfoBoxes - */ - void ShrinkBorrowing(); + * \brief Shrink the vectors in the FaceInfoBoxes + */ + void ShrinkBorrowing (); /** - * \brief Do the one-way extension - */ - void ComputeOneWayExtensions(); + * \brief Do the one-way extension + */ + void ComputeOneWayExtensions (); /** - * \brief Do the eight-ways extension - */ - void ComputeEightWaysExtensions(); + * \brief Do the eight-ways extension + */ + void ComputeEightWaysExtensions (); /** - * \brief Whenever an unstable cell cannot be extended we increase its area to be the minimal for stability. - * This is the method Benkler-Chavannes-Kuster method and it is less accurate than the regular ECT but it - * still works better than staircasing. (see https://ieeexplore.ieee.org/document/1638381) - * - * @param idim Integer indicating the dimension (x->0, y->1, z->2) for which the BCK correction is done - * - */ - void ApplyBCKCorrection(int idim); + * \brief Whenever an unstable cell cannot be extended we increase its area + * to be the minimal for stability. This is the method + * Benkler-Chavannes-Kuster method and it is less accurate than the regular + * ECT but it still works better than staircasing. (see + * https://ieeexplore.ieee.org/document/1638381) + * + * @param idim Integer indicating the dimension (x->0, y->1, z->2) for which + * the BCK correction is done + * + */ + void ApplyBCKCorrection (int idim); /** - * \brief Subtract the average of the cumulative sums of the preliminary current D - * from the current J (computed from D according to the Vay deposition scheme) + * \brief Subtract the average of the cumulative sums of the preliminary + * current D from the current J (computed from D according to the Vay + * deposition scheme) */ void PSATDSubtractCurrentPartialSumsAvg (); #ifdef WARPX_USE_FFT - auto& get_spectral_solver_fp (int lev) {return *spectral_solver_fp[lev];} + auto& + get_spectral_solver_fp (int lev) { + return *spectral_solver_fp[lev]; + } #endif - FiniteDifferenceSolver * get_pointer_fdtd_solver_fp (int lev) { return m_fdtd_solver_fp[lev].get(); } + FiniteDifferenceSolver* + get_pointer_fdtd_solver_fp (int lev) { + return m_fdtd_solver_fp[lev].get(); + } // Field container ablastr::fields::MultiFabRegister m_fields; -protected: - + protected: /** * \brief * This function initializes E, B, rho, and F, at all the levels * of the multifab. rho and F are initialized with 0. * The E and B fields are initialized using user-defined inputs. * The initialization type is set using "B_ext_grid_init_style" - * and "E_ext_grid_init_style". The initialization style is set to "default" - * if not explicitly defined by the user, and the E and B fields are - * initialized with E_external_grid and B_external_grid, respectively, each with - * a default value of 0. - * If the initialization type for the E and B field is "constant", - * then, the E and B fields at all the levels are initialized with - * user-defined values for E_external_grid and B_external_grid. - * If the initialization type for B-field is set to + * and "E_ext_grid_init_style". The initialization style is set to + * "default" if not explicitly defined by the user, and the E and B fields + * are initialized with E_external_grid and B_external_grid, respectively, + * each with a default value of 0. If the initialization type for the E and + * B field is "constant", then, the E and B fields at all the levels are + * initialized with user-defined values for E_external_grid and + * B_external_grid. If the initialization type for B-field is set to * "parse_ext_grid_function", then, the parser is used to read * Bx_external_grid_function(x,y,z), By_external_grid_function(x,y,z), * and Bz_external_grid_function(x,y,z). @@ -1060,7 +1295,8 @@ protected: * and Ex_external_grid_function(x,y,z). The parser for the E and B * initialization assumes that the function has three independent * variables, at max, namely, x, y, z. However, any number of constants - * can be used in the function used to define the E and B fields on the grid. + * can be used in the function used to define the E and B fields on the + * grid. */ void InitLevelData (int lev, amrex::Real time); @@ -1075,31 +1311,35 @@ protected: //! Make a new level from scratch using provided BoxArray and //! DistributionMapping. Only used during initialization. Called //! by AmrCoreInitFromScratch. - void MakeNewLevelFromScratch (int lev, amrex::Real time, const amrex::BoxArray& new_grids, - const amrex::DistributionMapping& new_dmap) final; + void + MakeNewLevelFromScratch (int lev, amrex::Real time, + const amrex::BoxArray& new_grids, + const amrex::DistributionMapping& new_dmap) final; //! Make a new level using provided BoxArray and //! DistributionMapping and fill with interpolated coarse level //! data. Called by AmrCore::regrid. - void MakeNewLevelFromCoarse (int /*lev*/, amrex::Real /*time*/, const amrex::BoxArray& /*ba*/, - const amrex::DistributionMapping& /*dm*/) final; + void + MakeNewLevelFromCoarse (int /*lev*/, amrex::Real /*time*/, + const amrex::BoxArray& /*ba*/, + const amrex::DistributionMapping& /*dm*/) final; //! Remake an existing level using provided BoxArray and //! DistributionMapping and fill with existing fine and coarse //! data. Called by AmrCore::regrid. void RemakeLevel (int lev, amrex::Real time, const amrex::BoxArray& ba, - const amrex::DistributionMapping& dm) final; + const amrex::DistributionMapping& dm) final; //! Delete level data. Called by AmrCore::regrid. void ClearLevel (int lev) final; -private: - + private: /** * \brief - * WarpX constructor. This method should not be called directly, but rather through - * the static member function MakeWarpX(). MakeWarpX() is called by GetInstance () - * if an instance of the WarpX class does not currently exist. + * WarpX constructor. This method should not be called directly, but rather + * through the static member function MakeWarpX(). MakeWarpX() is called by + * GetInstance () if an instance of the WarpX class does not currently + * exist. */ WarpX (); @@ -1112,13 +1352,18 @@ private: // Singleton is used when the code is run from python static WarpX* m_instance; - //! Complete the asynchronous broadcast of signal flags, and initiate a checkpoint if requested + //! Complete the asynchronous broadcast of signal flags, and initiate a + //! checkpoint if requested void HandleSignals (); - void FillBoundaryB (int lev, PatchType patch_type, amrex::IntVect ng, std::optional nodal_sync = std::nullopt); - void FillBoundaryE (int lev, PatchType patch_type, amrex::IntVect ng, std::optional nodal_sync = std::nullopt); - void FillBoundaryF (int lev, PatchType patch_type, amrex::IntVect ng, std::optional nodal_sync = std::nullopt); - void FillBoundaryG (int lev, PatchType patch_type, amrex::IntVect ng, std::optional nodal_sync = std::nullopt); + void FillBoundaryB (int lev, PatchType patch_type, amrex::IntVect ng, + std::optional nodal_sync = std::nullopt); + void FillBoundaryE (int lev, PatchType patch_type, amrex::IntVect ng, + std::optional nodal_sync = std::nullopt); + void FillBoundaryF (int lev, PatchType patch_type, amrex::IntVect ng, + std::optional nodal_sync = std::nullopt); + void FillBoundaryG (int lev, PatchType patch_type, amrex::IntVect ng, + std::optional nodal_sync = std::nullopt); void FillBoundaryB_avg (int lev, PatchType patch_type, amrex::IntVect ng); void FillBoundaryE_avg (int lev, PatchType patch_type, amrex::IntVect ng); @@ -1129,53 +1374,39 @@ private: void OneStep_sub1 (amrex::Real cur_time); /** - * \brief Perform one PIC iteration, with the multiple J deposition per time step + * \brief Perform one PIC iteration, with the multiple J deposition per time + * step */ void OneStep_multiJ (amrex::Real cur_time); void RestrictCurrentFromFineToCoarsePatch ( const ablastr::fields::MultiLevelVectorField& J_fp, - const ablastr::fields::MultiLevelVectorField& J_cp, - int lev); + const ablastr::fields::MultiLevelVectorField& J_cp, int lev); void AddCurrentFromFineLevelandSumBoundary ( const ablastr::fields::MultiLevelVectorField& J_fp, const ablastr::fields::MultiLevelVectorField& J_cp, - const ablastr::fields::MultiLevelVectorField& J_buffer, - int lev); + const ablastr::fields::MultiLevelVectorField& J_buffer, int lev); void StoreCurrent (int lev); void RestoreCurrent (int lev); - void ApplyFilterMF ( - const ablastr::fields::MultiLevelVectorField& mfvec, - int lev, - int idim); - void ApplyFilterMF ( - const ablastr::fields::MultiLevelVectorField& mfvec, - int lev); - void SumBoundaryJ ( - const ablastr::fields::MultiLevelVectorField& current, - int lev, - int idim, - const amrex::Periodicity& period); - void SumBoundaryJ ( - const ablastr::fields::MultiLevelVectorField& current, - int lev, - const amrex::Periodicity& period); - - void RestrictRhoFromFineToCoarsePatch (int lev ); + void ApplyFilterMF (const ablastr::fields::MultiLevelVectorField& mfvec, + int lev, int idim); + void ApplyFilterMF (const ablastr::fields::MultiLevelVectorField& mfvec, + int lev); + void SumBoundaryJ (const ablastr::fields::MultiLevelVectorField& current, + int lev, int idim, const amrex::Periodicity& period); + void SumBoundaryJ (const ablastr::fields::MultiLevelVectorField& current, + int lev, const amrex::Periodicity& period); + + void RestrictRhoFromFineToCoarsePatch (int lev); void ApplyFilterandSumBoundaryRho ( const ablastr::fields::MultiLevelScalarField& charge_fp, - const ablastr::fields::MultiLevelScalarField& charge_cp, - int lev, - PatchType patch_type, - int icomp, - int ncomp); + const ablastr::fields::MultiLevelScalarField& charge_cp, int lev, + PatchType patch_type, int icomp, int ncomp); void AddRhoFromFineLevelandSumBoundary ( const ablastr::fields::MultiLevelScalarField& charge_fp, const ablastr::fields::MultiLevelScalarField& charge_cp, - ablastr::fields::MultiLevelScalarField const & charge_buffer, - int lev, - int icomp, - int ncomp); + ablastr::fields::MultiLevelScalarField const& charge_buffer, int lev, + int icomp, int ncomp); void ReadParameters (); @@ -1189,7 +1420,8 @@ private: const amrex::DistributionMapping& dm); [[nodiscard]] amrex::DistributionMapping - GetRestartDMap (const std::string& chkfile, const amrex::BoxArray& ba, int lev) const; + GetRestartDMap (const std::string& chkfile, const amrex::BoxArray& ba, + int lev) const; void InitFromCheckpoint (); void PostRestart (); @@ -1204,79 +1436,90 @@ private: void InitNCICorrector (); /** - * \brief Check that the number of guard cells is smaller than the number of valid cells, - * for all available MultiFabs, and abort otherwise. + * \brief Check that the number of guard cells is smaller than the number of + * valid cells, for all available MultiFabs, and abort otherwise. */ - void CheckGuardCells(); + void CheckGuardCells (); /** - * \brief Checks for known numerical issues involving different WarpX modules + * \brief Checks for known numerical issues involving different WarpX + * modules */ - void CheckKnownIssues(); + void CheckKnownIssues (); void BuildBufferMasks (); - [[nodiscard]] const amrex::iMultiFab* getCurrentBufferMasks (int lev) const { + [[nodiscard]] const amrex::iMultiFab* + getCurrentBufferMasks (int lev) const { return current_buffer_masks[lev].get(); } - [[nodiscard]] const amrex::iMultiFab* getGatherBufferMasks (int lev) const - { + [[nodiscard]] const amrex::iMultiFab* + getGatherBufferMasks (int lev) const { return gather_buffer_masks[lev].get(); } /** - * \brief Allocates and initializes the stencil coefficients used for the finite-order centering - * of fields and currents, and stores them in the given device vectors. + * \brief Allocates and initializes the stencil coefficients used for the + * finite-order centering of fields and currents, and stores them in the + * given device vectors. * - * \param[in,out] device_centering_stencil_coeffs_x device vector where the stencil coefficients along x will be stored - * \param[in,out] device_centering_stencil_coeffs_y device vector where the stencil coefficients along y will be stored - * \param[in,out] device_centering_stencil_coeffs_z device vector where the stencil coefficients along z will be stored + * \param[in,out] device_centering_stencil_coeffs_x device vector where the + * stencil coefficients along x will be stored + * \param[in,out] device_centering_stencil_coeffs_y device vector where the + * stencil coefficients along y will be stored + * \param[in,out] device_centering_stencil_coeffs_z device vector where the + * stencil coefficients along z will be stored * \param[in] centering_nox order of the finite-order centering along x * \param[in] centering_noy order of the finite-order centering along y * \param[in] centering_noz order of the finite-order centering along z * \param[in] a_grid_type type of grid (collocated or not) */ - void AllocateCenteringCoefficients (amrex::Gpu::DeviceVector& device_centering_stencil_coeffs_x, - amrex::Gpu::DeviceVector& device_centering_stencil_coeffs_y, - amrex::Gpu::DeviceVector& device_centering_stencil_coeffs_z, - int centering_nox, - int centering_noy, - int centering_noz, - ablastr::utils::enums::GridType a_grid_type); - - void AllocLevelMFs (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm, + void + AllocateCenteringCoefficients (amrex::Gpu::DeviceVector& + device_centering_stencil_coeffs_x, + amrex::Gpu::DeviceVector& + device_centering_stencil_coeffs_y, + amrex::Gpu::DeviceVector& + device_centering_stencil_coeffs_z, + int centering_nox, int centering_noy, + int centering_noz, + ablastr::utils::enums::GridType a_grid_type); + + void AllocLevelMFs (int lev, const amrex::BoxArray& ba, + const amrex::DistributionMapping& dm, const amrex::IntVect& ngEB, amrex::IntVect& ngJ, const amrex::IntVect& ngRho, const amrex::IntVect& ngF, const amrex::IntVect& ngG, bool aux_is_nodal); #ifdef WARPX_USE_FFT -# ifdef WARPX_DIM_RZ - void AllocLevelSpectralSolverRZ (amrex::Vector>& spectral_solver, - int lev, - const amrex::BoxArray& realspace_ba, - const amrex::DistributionMapping& dm, - const std::array& dx); -# else - void AllocLevelSpectralSolver (amrex::Vector>& spectral_solver, - int lev, - const amrex::BoxArray& realspace_ba, - const amrex::DistributionMapping& dm, - const std::array& dx, - bool pml_flag=false); -# endif +#ifdef WARPX_DIM_RZ + void AllocLevelSpectralSolverRZ ( + amrex::Vector>& spectral_solver, + int lev, const amrex::BoxArray& realspace_ba, + const amrex::DistributionMapping& dm, + const std::array& dx); +#else + void AllocLevelSpectralSolver ( + amrex::Vector>& spectral_solver, + int lev, const amrex::BoxArray& realspace_ba, + const amrex::DistributionMapping& dm, + const std::array& dx, bool pml_flag = false); +#endif #endif //! Author of an input file / simulation setup std::string m_authors; - amrex::Vector istep; // which step? - amrex::Vector nsubsteps; // how many substeps on each level? + amrex::Vector istep; // which step? + amrex::Vector nsubsteps; // how many substeps on each level? amrex::Vector t_new; amrex::Vector t_old; amrex::Vector dt; - utils::parser::IntervalsParser m_dt_update_interval = utils::parser::IntervalsParser{}; // How often to update the timestep when using adaptive timestepping + utils::parser::IntervalsParser m_dt_update_interval = + utils::parser::IntervalsParser{}; // How often to update the timestep + // when using adaptive timestepping bool m_safe_guard_cells = false; @@ -1288,32 +1531,42 @@ private: bool do_fluid_species = false; std::unique_ptr myfl; - //! Integer that corresponds to electromagnetic Maxwell solver (vacuum - 0, macroscopic - 1) + //! Integer that corresponds to electromagnetic Maxwell solver (vacuum - 0, + //! macroscopic - 1) MediumForEM m_em_solver_medium = MediumForEM::Default; // // Fields: First array for level, second for direction // - // Masks for computing dot product and global moments of fields when using grids that - // have shared locations across different ranks (e.g., a Yee grid) - mutable amrex::Vector,3 > > Efield_dotMask; - mutable amrex::Vector,3 > > Bfield_dotMask; - mutable amrex::Vector,3 > > Afield_dotMask; - mutable amrex::Vector< std::unique_ptr > phi_dotMask; - - /** EB: Flag to indicate whether a gridpoint is inside the embedded boundary and therefore - * whether the E or B should not be updated. (One array per level and per direction, due to staggering) - */ - amrex::Vector,3 > > m_eb_update_E; - amrex::Vector,3 > > m_eb_update_B; - - /** EB: Mask that indicates whether a particle should use its regular shape factor (mask set to 0) - * or a reduced, order-1 shape factor instead (mask set to 1) in a given cell, when depositing charge/current. - * The flag is typically set to 1 in cells that are close to the embedded boundary, in order to avoid - * errors in charge conservation when a particle is too close to the embedded boundary. - */ - amrex::Vector< std::unique_ptr > m_eb_reduce_particle_shape; + // Masks for computing dot product and global moments of fields when using + // grids that have shared locations across different ranks (e.g., a Yee + // grid) + mutable amrex::Vector, 3>> + Efield_dotMask; + mutable amrex::Vector, 3>> + Bfield_dotMask; + mutable amrex::Vector, 3>> + Afield_dotMask; + mutable amrex::Vector> phi_dotMask; + + /** EB: Flag to indicate whether a gridpoint is inside the embedded boundary + * and therefore whether the E or B should not be updated. (One array per + * level and per direction, due to staggering) + */ + amrex::Vector, 3>> + m_eb_update_E; + amrex::Vector, 3>> + m_eb_update_B; + + /** EB: Mask that indicates whether a particle should use its regular shape + * factor (mask set to 0) or a reduced, order-1 shape factor instead (mask + * set to 1) in a given cell, when depositing charge/current. The flag is + * typically set to 1 in cells that are close to the embedded boundary, in + * order to avoid errors in charge conservation when a particle is too close + * to the embedded boundary. + */ + amrex::Vector> m_eb_reduce_particle_shape; /** EB: for every mesh face flag_info_face contains a: * * 0 if the face needs to be extended @@ -1321,23 +1574,28 @@ private: * * 2 if the face is actually intruded by other face * It is initialized in WarpX::MarkExtensionCells * This is only used for the ECT solver.*/ - amrex::Vector, 3 > > m_flag_info_face; + amrex::Vector, 3>> + m_flag_info_face; /** EB: for every mesh face face flag_ext_face contains a: * * 1 if the face needs to be extended * * 0 otherwise - * It is initialized in WarpX::MarkExtensionCells and then modified in WarpX::ComputeOneWayExtensions - * and in WarpX::ComputeEightWaysExtensions + * It is initialized in WarpX::MarkExtensionCells and then modified in + * WarpX::ComputeOneWayExtensions and in WarpX::ComputeEightWaysExtensions * This is only used for the ECT solver.*/ - amrex::Vector, 3 > > m_flag_ext_face; + amrex::Vector, 3>> + m_flag_ext_face; - /** EB: m_borrowing contains the info about the enlarged cells, i.e. for every enlarged cell it - * contains the info of which neighbors are being intruded (and the amount of borrowed area). - * This is only used for the ECT solver.*/ - amrex::Vector >, 3 > > m_borrowing; + /** EB: m_borrowing contains the info about the enlarged cells, i.e. for + * every enlarged cell it contains the info of which neighbors are being + * intruded (and the amount of borrowed area). This is only used for the ECT + * solver.*/ + amrex::Vector< + std::array>, 3>> + m_borrowing; // Copy of the coarse aux - amrex::Vector > current_buffer_masks; - amrex::Vector > gather_buffer_masks; + amrex::Vector> current_buffer_masks; + amrex::Vector> gather_buffer_masks; // PML int do_pml = 0; @@ -1352,9 +1610,9 @@ private: bool do_pml_divb_cleaning; // default set in WarpX.cpp amrex::Vector do_pml_Lo; amrex::Vector do_pml_Hi; - amrex::Vector > pml; + amrex::Vector> pml; #if (defined WARPX_DIM_RZ) && (defined WARPX_USE_FFT) - amrex::Vector > pml_rz; + amrex::Vector> pml_rz; #endif amrex::Real v_particle_pml; @@ -1384,8 +1642,8 @@ private: // whether to use subcycling bool m_do_subcycling = false; - //! Flag whether the Verboncoeur correction is applied to the current and charge density - //! on the axis when using RZ. + //! Flag whether the Verboncoeur correction is applied to the current and + //! charge density on the axis when using RZ. bool m_verboncoeur_axis_correction = true; // Macroscopic properties @@ -1398,38 +1656,40 @@ private: std::unique_ptr m_hybrid_pic_model; // Load balancing - /** Load balancing intervals that reads the "load_balance_intervals" string int the input file - * for getting steps at which load balancing is performed */ + /** Load balancing intervals that reads the "load_balance_intervals" string + * int the input file for getting steps at which load balancing is performed + */ utils::parser::IntervalsParser load_balance_intervals; /** Collection of LayoutData to keep track of weights used in load balancing - * routines. Contains timer-based or heuristic-based costs depending on input option */ - amrex::Vector > > costs; + * routines. Contains timer-based or heuristic-based costs depending on + * input option */ + amrex::Vector>> costs; /** Load balance with 'space filling curve' strategy. */ int load_balance_with_sfc = 0; - /** Controls the maximum number of boxes that can be assigned to a rank during - * load balance via the 'knapsack' strategy; e.g., if there are 4 boxes per rank, - * `load_balance_knapsack_factor=2` limits the maximum number of boxes that can - * be assigned to a rank to 8. */ + /** Controls the maximum number of boxes that can be assigned to a rank + * during load balance via the 'knapsack' strategy; e.g., if there are 4 + * boxes per rank, `load_balance_knapsack_factor=2` limits the maximum + * number of boxes that can be assigned to a rank to 8. */ amrex::Real load_balance_knapsack_factor = amrex::Real(1.24); /** Threshold value that controls whether to adopt the proposed distribution * mapping during load balancing. The new distribution mapping is adopted * if the ratio of proposed distribution mapping efficiency to current - * distribution mapping efficiency is larger than the threshold; 'efficiency' - * here means the average cost per MPI rank. */ + * distribution mapping efficiency is larger than the threshold; + * 'efficiency' here means the average cost per MPI rank. */ amrex::Real load_balance_efficiency_ratio_threshold = amrex::Real(1.1); /** Current load balance efficiency for each level. */ amrex::Vector load_balance_efficiency; /** Weight factor for cells in `Heuristic` costs update. * Default values on GPU are determined from single-GPU tests on Summit. * The problem setup for these tests is an empty (i.e. no particles) domain - * of size 256 by 256 by 256 cells, from which the average time per iteration - * per cell is computed. */ + * of size 256 by 256 by 256 cells, from which the average time per + * iteration per cell is computed. */ amrex::Real costs_heuristic_cells_wt = amrex::Real(0); /** Weight factor for particles in `Heuristic` costs update. * Default values on GPU are determined from single-GPU tests on Summit. * The problem setup for these tests is a high-ppc (27 particles per cell) - * uniform plasma on a domain of size 128 by 128 by 128, from which the approximate - * time per iteration per particle is computed. */ + * uniform plasma on a domain of size 128 by 128 by 128, from which the + * approximate time per iteration per particle is computed. */ amrex::Real costs_heuristic_particles_wt = amrex::Real(0); // Determines timesteps for override sync @@ -1440,12 +1700,12 @@ private: bool use_hybrid_QED = false; - int max_step = std::numeric_limits::max(); + int max_step = std::numeric_limits::max(); amrex::Real stop_time = std::numeric_limits::max(); - //! If specified, the maximum number of iterations is computed automatically so that - //! the lower end of the simulation domain along z reaches #zmax_plasma_to_compute_max_step - //! in the boosted frame + //! If specified, the maximum number of iterations is computed automatically + //! so that the lower end of the simulation domain along z reaches + //! #zmax_plasma_to_compute_max_step in the boosted frame std::optional m_zmax_plasma_to_compute_max_step = std::nullopt; int regrid_int = -1; @@ -1454,11 +1714,14 @@ private: std::string restart_chkfile; - /** When `true`, write the diagnostics after restart at the time of the restart. */ + /** When `true`, write the diagnostics after restart at the time of the + * restart. */ bool write_diagnostics_on_restart = false; - amrex::VisMF::Header::Version plotfile_headerversion = amrex::VisMF::Header::Version_v1; - amrex::VisMF::Header::Version slice_plotfile_headerversion = amrex::VisMF::Header::Version_v1; + amrex::VisMF::Header::Version plotfile_headerversion = + amrex::VisMF::Header::Version_v1; + amrex::VisMF::Header::Version slice_plotfile_headerversion = + amrex::VisMF::Header::Version_v1; bool synchronize_velocity_for_diagnostics = false; @@ -1491,8 +1754,9 @@ private: int noy_fft = 16; int noz_fft = 16; - //! Solve Poisson equation when loading an external magnetic field to clean divergence - //! This is useful to remove errors that could lead to non-zero B field divergence + //! Solve Poisson equation when loading an external magnetic field to clean + //! divergence This is useful to remove errors that could lead to non-zero B + //! field divergence bool m_do_divb_cleaning_external = false; //! Domain decomposition on Level 0 @@ -1502,22 +1766,24 @@ private: std::unique_ptr m_particle_boundary_buffer; // Accelerator lattice elements - amrex::Vector< std::unique_ptr > m_accelerator_lattice; + amrex::Vector> m_accelerator_lattice; // // Embedded Boundary // // Factory for field data - amrex::Vector > > m_field_factory; + amrex::Vector>> + m_field_factory; [[nodiscard]] - amrex::FabFactory const& fieldFactory (int lev) const noexcept - { + amrex::FabFactory const& + fieldFactory (int lev) const noexcept { return *m_field_factory[lev]; } - /** Stop the simulation at the end of the current step due to a received Unix signal? + /** Stop the simulation at the end of the current step due to a received + * Unix signal? */ bool m_exit_loop_due_to_interrupt_signal = false; @@ -1534,7 +1800,8 @@ private: * @param cur_time current time * @param num_moved number of cells the moving window moved */ - void HandleParticlesAtBoundaries (int step, amrex::Real cur_time, int num_moved); + void HandleParticlesAtBoundaries (int step, amrex::Real cur_time, + int num_moved); /** Update the E and B fields in the explicit em PIC scheme. * @@ -1567,13 +1834,17 @@ private: * \brief Backward FFT of averaged E,B on all mesh refinement levels * * \param E_avg_fp Vector of three-dimensional arrays (for each level) - * storing the fine patch averaged electric field to be transformed + * storing the fine patch averaged electric field to be + * transformed * \param B_avg_fp Vector of three-dimensional arrays (for each level) - * storing the fine patch averaged magnetic field to be transformed + * storing the fine patch averaged magnetic field to be + * transformed * \param E_avg_cp Vector of three-dimensional arrays (for each level) - * storing the coarse patch averaged electric field to be transformed + * storing the coarse patch averaged electric field to be + * transformed * \param B_avg_cp Vector of three-dimensional arrays (for each level) - * storing the coarse patch averaged magnetic field to be transformed + * storing the coarse patch averaged magnetic field to be + * transformed */ void PSATDBackwardTransformEBavg ( ablastr::fields::MultiLevelVectorField const& E_avg_fp, @@ -1590,12 +1861,12 @@ private: * \param J_cp Vector of three-dimensional arrays (for each level) * storing the coarse patch current to be transformed * \param[in] apply_kspace_filter Control whether to apply filtering - * (only used in RZ geometry to avoid double filtering) + * (only used in RZ geometry to avoid double + * filtering) */ - void PSATDForwardTransformJ ( - std::string const & J_fp_string, - std::string const & J_cp_string, - bool apply_kspace_filter=true); + void PSATDForwardTransformJ (std::string const& J_fp_string, + std::string const& J_cp_string, + bool apply_kspace_filter = true); /** * \brief Backward FFT of J on all mesh refinement levels @@ -1605,28 +1876,32 @@ private: * \param J_cp Vector of three-dimensional arrays (for each level) * storing the coarse patch current to be transformed */ - void PSATDBackwardTransformJ ( - std::string const & J_fp_string, - std::string const & J_cp_string); + void PSATDBackwardTransformJ (std::string const& J_fp_string, + std::string const& J_cp_string); /** * \brief Forward FFT of rho on all mesh refinement levels, * with k-space filtering (if needed) * - * \param charge_fp Vector (for each level) storing the fine patch charge to be transformed - * \param charge_cp Vector (for each level) storing the coarse patch charge to be transformed + * \param charge_fp Vector (for each level) storing the fine patch charge to + * be transformed + * \param charge_cp Vector (for each level) storing the coarse patch charge + * to be transformed * \param[in] icomp index of fourth component (0 for rho_old, 1 for rho_new) - * \param[in] dcomp index of spectral component (0 for rho_old, 1 for rho_new) + * \param[in] dcomp index of spectral component (0 for rho_old, 1 for + * rho_new) * \param[in] apply_kspace_filter Control whether to apply filtering - * (only used in RZ geometry to avoid double filtering) + * (only used in RZ geometry to avoid double + * filtering) */ - void PSATDForwardTransformRho ( - std::string const & charge_fp_string, - std::string const & charge_cp_string, - int icomp, int dcomp, bool apply_kspace_filter=true); + void PSATDForwardTransformRho (std::string const& charge_fp_string, + std::string const& charge_cp_string, + int icomp, int dcomp, + bool apply_kspace_filter = true); /** - * \brief Copy rho_new to rho_old in spectral space (when rho is linear in time) + * \brief Copy rho_new to rho_old in spectral space (when rho is linear in + * time) */ void PSATDMoveRhoNewToRhoOld (); @@ -1672,13 +1947,13 @@ private: */ void PSATDEraseAverageFields (); -# ifdef WARPX_DIM_RZ - amrex::Vector> spectral_solver_fp; - amrex::Vector> spectral_solver_cp; -# else - amrex::Vector> spectral_solver_fp; - amrex::Vector> spectral_solver_cp; -# endif +#ifdef WARPX_DIM_RZ + amrex::Vector> spectral_solver_fp; + amrex::Vector> spectral_solver_cp; +#else + amrex::Vector> spectral_solver_fp; + amrex::Vector> spectral_solver_cp; +#endif #endif @@ -1687,7 +1962,6 @@ private: // implicit solver object std::unique_ptr m_implicit_solver; - }; #endif diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index f129987b991..072d2e13168 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -17,31 +17,31 @@ #include "Diagnostics/ReducedDiags/MultiReducedDiags.H" #include "EmbeddedBoundary/Enabled.H" #include "EmbeddedBoundary/WarpXFaceInfoBox.H" +#include "FieldSolver/ElectrostaticSolvers/EffectivePotentialES.H" #include "FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.H" #include "FieldSolver/ElectrostaticSolvers/LabFrameExplicitES.H" #include "FieldSolver/ElectrostaticSolvers/RelativisticExplicitES.H" -#include "FieldSolver/ElectrostaticSolvers/EffectivePotentialES.H" #include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H" -#include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H" #include "FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H" +#include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H" #ifdef WARPX_USE_FFT -# include "FieldSolver/SpectralSolver/SpectralKSpace.H" -# ifdef WARPX_DIM_RZ -# include "FieldSolver/SpectralSolver/SpectralSolverRZ.H" -# include "BoundaryConditions/PML_RZ.H" -# else -# include "FieldSolver/SpectralSolver/SpectralSolver.H" -# endif // RZ ifdef +#include "FieldSolver/SpectralSolver/SpectralKSpace.H" +#ifdef WARPX_DIM_RZ +#include "BoundaryConditions/PML_RZ.H" +#include "FieldSolver/SpectralSolver/SpectralSolverRZ.H" +#else +#include "FieldSolver/SpectralSolver/SpectralSolver.H" +#endif // RZ ifdef #endif // use PSATD ifdef +#include "AcceleratorLattice/AcceleratorLattice.H" #include "FieldSolver/WarpX_FDTD.H" #include "Filter/NCIGodfreyFilter.H" +#include "Fluids/MultiFluidContainer.H" +#include "Fluids/WarpXFluidContainer.H" #include "Initialization/ExternalField.H" #include "Initialization/WarpXInit.H" #include "Particles/MultiParticleContainer.H" -#include "Fluids/MultiFluidContainer.H" -#include "Fluids/WarpXFluidContainer.H" #include "Particles/ParticleBoundaryBuffer.H" -#include "AcceleratorLattice/AcceleratorLattice.H" #include "Utils/TextMsg.H" #include "Utils/WarpXAlgorithmSelection.H" #include "Utils/WarpXConst.H" @@ -55,7 +55,7 @@ #include #ifdef AMREX_USE_SENSEI_INSITU -# include +#include #endif #include #include @@ -63,8 +63,8 @@ #include #include #ifdef AMREX_USE_EB -# include -# include +#include +#include #endif #include #include @@ -106,12 +106,12 @@ Real WarpX::moving_window_v = std::numeric_limits::max(); bool WarpX::fft_do_time_averaging = false; -amrex::IntVect WarpX::m_fill_guards_fields = amrex::IntVect(0); +amrex::IntVect WarpX::m_fill_guards_fields = amrex::IntVect(0); amrex::IntVect WarpX::m_fill_guards_current = amrex::IntVect(0); Real WarpX::gamma_boost = 1._rt; Real WarpX::beta_boost = 0._rt; -Vector WarpX::boost_direction = {0,0,0}; +Vector WarpX::boost_direction = {0, 0, 0}; bool WarpX::compute_max_step_from_btd = false; Real WarpX::zmin_domain_boost_step_0 = 0._rt; @@ -124,12 +124,13 @@ bool WarpX::do_single_precision_comms = false; bool WarpX::do_shared_mem_charge_deposition = false; bool WarpX::do_shared_mem_current_deposition = false; #if defined(WARPX_DIM_3D) -amrex::IntVect WarpX::shared_tilesize(AMREX_D_DECL(6,6,8)); +amrex::IntVect WarpX::shared_tilesize(AMREX_D_DECL(6, 6, 8)); #elif (AMREX_SPACEDIM == 2) -amrex::IntVect WarpX::shared_tilesize(AMREX_D_DECL(14,14,0)); +amrex::IntVect WarpX::shared_tilesize(AMREX_D_DECL(14, 14, 0)); #else -//Have not experimented with good tilesize here because expect use case to be low -amrex::IntVect WarpX::shared_tilesize(AMREX_D_DECL(1,1,1)); +// Have not experimented with good tilesize here because expect use case to be +// low +amrex::IntVect WarpX::shared_tilesize(AMREX_D_DECL(1, 1, 1)); #endif int WarpX::shared_mem_current_tpb = 128; @@ -155,14 +156,14 @@ bool WarpX::use_fdtd_nci_corr = false; bool WarpX::galerkin_interpolation = true; bool WarpX::use_filter = true; -bool WarpX::use_kspace_filter = true; +bool WarpX::use_kspace_filter = true; bool WarpX::use_filter_compensation = false; bool WarpX::serialize_initial_conditions = false; -bool WarpX::refine_plasma = false; +bool WarpX::refine_plasma = false; utils::parser::IntervalsParser WarpX::sort_intervals; -amrex::IntVect WarpX::sort_bin_size(AMREX_D_DECL(1,1,1)); +amrex::IntVect WarpX::sort_bin_size(AMREX_D_DECL(1, 1, 1)); #if defined(AMREX_USE_CUDA) bool WarpX::sort_particles_for_deposition = true; @@ -170,7 +171,7 @@ bool WarpX::sort_particles_for_deposition = true; bool WarpX::sort_particles_for_deposition = false; #endif -amrex::IntVect WarpX::sort_idx_type(AMREX_D_DECL(0,0,0)); +amrex::IntVect WarpX::sort_idx_type(AMREX_D_DECL(0, 0, 0)); bool WarpX::do_dynamic_scheduling = true; @@ -186,50 +187,54 @@ amrex::IntVect m_rho_nodal_flag; WarpX* WarpX::m_instance = nullptr; -namespace -{ - - [[nodiscard]] bool - isAnyBoundaryPML( - const amrex::Array& field_boundary_lo, - const amrex::Array& field_boundary_hi) - { - constexpr auto is_pml = [](const FieldBoundaryType fbt) {return (fbt == FieldBoundaryType::PML);}; - const auto is_any_pml = - std::any_of(field_boundary_lo.begin(), field_boundary_lo.end(), is_pml) || - std::any_of(field_boundary_hi.begin(), field_boundary_hi.end(), is_pml); - return is_any_pml; - } +namespace { + +[[nodiscard]] bool +isAnyBoundaryPML ( + const amrex::Array& field_boundary_lo, + const amrex::Array& field_boundary_hi) { + constexpr auto is_pml = [] (const FieldBoundaryType fbt) { + return (fbt == FieldBoundaryType::PML); + }; + const auto is_any_pml = + std::any_of(field_boundary_lo.begin(), field_boundary_lo.end(), + is_pml) || + std::any_of(field_boundary_hi.begin(), field_boundary_hi.end(), is_pml); + return is_any_pml; +} - /** - * \brief - * Set the dotMask container - */ - void SetDotMask( std::unique_ptr& field_dotMask, - ablastr::fields::ConstScalarField const& field, - amrex::Periodicity const& periodicity) +/** + * \brief + * Set the dotMask container + */ +void +SetDotMask (std::unique_ptr& field_dotMask, + ablastr::fields::ConstScalarField const& field, + amrex::Periodicity const& periodicity) - { - // Define the dot mask for this field_type needed to properly compute dotProduct() - // for field values that have shared locations on different MPI ranks - if (field_dotMask != nullptr) { return; } +{ + // Define the dot mask for this field_type needed to properly compute + // dotProduct() for field values that have shared locations on different MPI + // ranks + if (field_dotMask != nullptr) { + return; + } - const auto& this_ba = field->boxArray(); - const auto tmp = amrex::MultiFab{ - this_ba, field->DistributionMap(), - 1, 0, amrex::MFInfo().SetAlloc(false)}; + const auto& this_ba = field->boxArray(); + const auto tmp = amrex::MultiFab{this_ba, field->DistributionMap(), 1, 0, + amrex::MFInfo().SetAlloc(false)}; - field_dotMask = tmp.OwnerMask(periodicity); - } + field_dotMask = tmp.OwnerMask(periodicity); } +} // namespace -void WarpX::MakeWarpX () -{ +void +WarpX::MakeWarpX () { warpx::initialization::check_dims(); - ReadMovingWindowParameters( - do_moving_window, start_moving_window_step, end_moving_window_step, - moving_window_dir, moving_window_v); + ReadMovingWindowParameters(do_moving_window, start_moving_window_step, + end_moving_window_step, moving_window_dir, + moving_window_v); ConvertLabParamsToBoost(); ReadBCParams(); @@ -242,8 +247,7 @@ void WarpX::MakeWarpX () } WarpX& -WarpX::GetInstance () -{ +WarpX::GetInstance () { if (!m_instance) { MakeWarpX(); } @@ -251,29 +255,28 @@ WarpX::GetInstance () } void -WarpX::ResetInstance () -{ - if (m_instance){ +WarpX::ResetInstance () { + if (m_instance) { delete m_instance; m_instance = nullptr; } } void -WarpX::Finalize() -{ +WarpX::Finalize () { WarpX::ResetInstance(); } -WarpX::WarpX () -{ +WarpX::WarpX () { warpx::initialization::initialize_warning_manager(); ReadParameters(); BackwardCompatibility(); - if (EB::enabled()) { InitEB(); } + if (EB::enabled()) { + InitEB(); + } ablastr::utils::SignalHandling::InitSignalHandling(); @@ -294,24 +297,23 @@ WarpX::WarpX () // Loop over species (particles and lasers) // and set current injection position per species - if (do_moving_window){ + if (do_moving_window) { const int n_containers = mypc->nContainers(); - for (int i=0; iGetParticleContainer(i); // Storing injection position for all species, regardless of whether // they are continuously injected, since it makes looping over the - // elements of current_injection_position easier elsewhere in the code. - if (moving_window_v > 0._rt) - { + // elements of current_injection_position easier elsewhere in the + // code. + if (moving_window_v > 0._rt) { // Inject particles continuously from the right end of the box - pc.m_current_injection_position = geom[0].ProbHi(moving_window_dir); - } - else if (moving_window_v < 0._rt) - { + pc.m_current_injection_position = + geom[0].ProbHi(moving_window_dir); + } else if (moving_window_v < 0._rt) { // Inject particles continuously from the left end of the box - pc.m_current_injection_position = geom[0].ProbLo(moving_window_dir); + pc.m_current_injection_position = + geom[0].ProbLo(moving_window_dir); } } } @@ -338,23 +340,24 @@ WarpX::WarpX () m_borrowing.resize(nlevs_max); // Create Electrostatic Solver object if needed - if ((WarpX::electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame) - || (WarpX::electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic)) - { - m_electrostatic_solver = std::make_unique(nlevs_max); + if ((WarpX::electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame) || + (WarpX::electrostatic_solver_id == + ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic)) { + m_electrostatic_solver = + std::make_unique(nlevs_max); } // Initialize the effective potential electrostatic solver if required - else if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameEffectivePotential) - { - m_electrostatic_solver = std::make_unique(nlevs_max); - } - else - { - m_electrostatic_solver = std::make_unique(nlevs_max); + else if (electrostatic_solver_id == + ElectrostaticSolverAlgo::LabFrameEffectivePotential) { + m_electrostatic_solver = + std::make_unique(nlevs_max); + } else { + m_electrostatic_solver = + std::make_unique(nlevs_max); } - if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC) - { + if (WarpX::electromagnetic_solver_id == + ElectromagneticSolverAlgo::HybridPIC) { // Create hybrid-PIC model object if needed m_hybrid_pic_model = std::make_unique(); } @@ -383,54 +386,53 @@ WarpX::WarpX () // Set default values for particle and cell weights for costs update; // Default values listed here for the case AMREX_USE_GPU are determined // from single-GPU tests on Summit. - if (costs_heuristic_cells_wt<=0. && costs_heuristic_particles_wt<=0. - && WarpX::load_balance_costs_update_algo==LoadBalanceCostsUpdateAlgo::Heuristic) - { + if (costs_heuristic_cells_wt <= 0. && costs_heuristic_particles_wt <= 0. && + WarpX::load_balance_costs_update_algo == + LoadBalanceCostsUpdateAlgo::Heuristic) { #ifdef AMREX_USE_GPU - if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { - switch (WarpX::nox) - { - case 1: - costs_heuristic_cells_wt = 0.575_rt; - costs_heuristic_particles_wt = 0.425_rt; - break; - case 2: - costs_heuristic_cells_wt = 0.405_rt; - costs_heuristic_particles_wt = 0.595_rt; - break; - case 3: - costs_heuristic_cells_wt = 0.250_rt; - costs_heuristic_particles_wt = 0.750_rt; - break; - case 4: - // this is only a guess - costs_heuristic_cells_wt = 0.200_rt; - costs_heuristic_particles_wt = 0.800_rt; - break; + if (WarpX::electromagnetic_solver_id == + ElectromagneticSolverAlgo::PSATD) { + switch (WarpX::nox) { + case 1: + costs_heuristic_cells_wt = 0.575_rt; + costs_heuristic_particles_wt = 0.425_rt; + break; + case 2: + costs_heuristic_cells_wt = 0.405_rt; + costs_heuristic_particles_wt = 0.595_rt; + break; + case 3: + costs_heuristic_cells_wt = 0.250_rt; + costs_heuristic_particles_wt = 0.750_rt; + break; + case 4: + // this is only a guess + costs_heuristic_cells_wt = 0.200_rt; + costs_heuristic_particles_wt = 0.800_rt; + break; } } else { // FDTD - switch (WarpX::nox) - { - case 1: - costs_heuristic_cells_wt = 0.401_rt; - costs_heuristic_particles_wt = 0.599_rt; - break; - case 2: - costs_heuristic_cells_wt = 0.268_rt; - costs_heuristic_particles_wt = 0.732_rt; - break; - case 3: - costs_heuristic_cells_wt = 0.145_rt; - costs_heuristic_particles_wt = 0.855_rt; - break; - case 4: - // this is only a guess - costs_heuristic_cells_wt = 0.100_rt; - costs_heuristic_particles_wt = 0.900_rt; - break; + switch (WarpX::nox) { + case 1: + costs_heuristic_cells_wt = 0.401_rt; + costs_heuristic_particles_wt = 0.599_rt; + break; + case 2: + costs_heuristic_cells_wt = 0.268_rt; + costs_heuristic_particles_wt = 0.732_rt; + break; + case 3: + costs_heuristic_cells_wt = 0.145_rt; + costs_heuristic_particles_wt = 0.855_rt; + break; + case 4: + // this is only a guess + costs_heuristic_cells_wt = 0.100_rt; + costs_heuristic_particles_wt = 0.900_rt; + break; } } -#else // CPU +#else // CPU costs_heuristic_cells_wt = 0.1_rt; costs_heuristic_particles_wt = 0.9_rt; #endif // AMREX_USE_GPU @@ -468,22 +470,20 @@ WarpX::WarpX () } m_accelerator_lattice.resize(nlevs_max); - } -WarpX::~WarpX () -{ - const int nlevs_max = maxLevel() +1; +WarpX::~WarpX () { + const int nlevs_max = maxLevel() + 1; for (int lev = 0; lev < nlevs_max; ++lev) { ClearLevel(lev); } } void -WarpX::ReadParameters () -{ +WarpX::ReadParameters () { { - const ParmParse pp;// Traditionally, max_step and stop_time do not have prefix. + const ParmParse + pp; // Traditionally, max_step and stop_time do not have prefix. utils::parser::queryWithParser(pp, "max_step", max_step); utils::parser::queryWithParser(pp, "stop_time", stop_time); pp.query("authors", m_authors); @@ -497,20 +497,24 @@ WarpX::ReadParameters () { const ParmParse pp_algo("algo"); - pp_algo.query_enum_sloppy("maxwell_solver", electromagnetic_solver_id, "-_"); - if (electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT && !EB::enabled()) { - throw std::runtime_error("ECP Solver requires to enable embedded boundaries at runtime."); + pp_algo.query_enum_sloppy("maxwell_solver", electromagnetic_solver_id, + "-_"); + if (electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT && + !EB::enabled()) { + throw std::runtime_error("ECP Solver requires to enable embedded " + "boundaries at runtime."); } #ifdef WARPX_DIM_RZ - if (electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) - { - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(Geom(0).ProbLo(0) == 0., - "Lower bound of radial coordinate (prob_lo[0]) with RZ PSATD solver must be zero"); - } - else - { - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(Geom(0).ProbLo(0) >= 0., - "Lower bound of radial coordinate (prob_lo[0]) with RZ FDTD solver must be non-negative"); + if (electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + Geom(0).ProbLo(0) == 0., + "Lower bound of radial coordinate (prob_lo[0]) with RZ PSATD " + "solver must be zero"); + } else { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + Geom(0).ProbLo(0) >= 0., + "Lower bound of radial coordinate (prob_lo[0]) with RZ FDTD " + "solver must be non-negative"); } #endif @@ -521,24 +525,28 @@ WarpX::ReadParameters () ParmParse const pp_warpx("warpx"); std::vector numprocs_in; - utils::parser::queryArrWithParser( - pp_warpx, "numprocs", numprocs_in, 0, AMREX_SPACEDIM); + utils::parser::queryArrWithParser(pp_warpx, "numprocs", numprocs_in, 0, + AMREX_SPACEDIM); if (not numprocs_in.empty()) { #ifdef WARPX_DIM_RZ if (electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(numprocs_in[0] == 1, - "Domain decomposition in RZ with spectral solvers works only along z direction"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + numprocs_in[0] == 1, + "Domain decomposition in RZ with spectral solvers works " + "only along z direction"); } #endif - WARPX_ALWAYS_ASSERT_WITH_MESSAGE - (numprocs_in.size() == AMREX_SPACEDIM, - "warpx.numprocs, if specified, must have AMREX_SPACEDIM numbers"); - WARPX_ALWAYS_ASSERT_WITH_MESSAGE - (ParallelDescriptor::NProcs() == AMREX_D_TERM(numprocs_in[0], + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + numprocs_in.size() == AMREX_SPACEDIM, + "warpx.numprocs, if specified, must have AMREX_SPACEDIM " + "numbers"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + ParallelDescriptor::NProcs() == AMREX_D_TERM(numprocs_in[0], *numprocs_in[1], *numprocs_in[2]), - "warpx.numprocs, if specified, its product must be equal to the number of processes"); + "warpx.numprocs, if specified, its product must be equal to " + "the number of processes"); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { numprocs[idim] = numprocs_in[idim]; } @@ -549,14 +557,16 @@ WarpX::ReadParameters () pp_warpx.queryarr("break_signals", signals_in); #if defined(__linux__) || defined(__APPLE__) - for (const std::string &str : signals_in) { + for (const std::string& str : signals_in) { const int sig = SignalHandling::parseSignalNameToNumber(str); - SignalHandling::signal_conf_requests[SignalHandling::SIGNAL_REQUESTS_BREAK][sig] = true; + SignalHandling::signal_conf_requests + [SignalHandling::SIGNAL_REQUESTS_BREAK][sig] = true; } signals_in.clear(); #else - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(signals_in.empty(), - "Signal handling requested in input, but is not supported on this platform"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + signals_in.empty(), "Signal handling requested in input, but is " + "not supported on this platform"); #endif bool have_checkpoint_diagnostic = false; @@ -565,7 +575,7 @@ WarpX::ReadParameters () std::vector diags_names; pp.queryarr("diags_names", diags_names); - for (const auto &diag : diags_names) { + for (const auto& diag : diags_names) { const ParmParse dd(diag); std::string format; dd.query("format", format); @@ -575,33 +585,38 @@ WarpX::ReadParameters () } } - pp_warpx.query("write_diagnostics_on_restart", write_diagnostics_on_restart); + pp_warpx.query("write_diagnostics_on_restart", + write_diagnostics_on_restart); pp_warpx.queryarr("checkpoint_signals", signals_in); #if defined(__linux__) || defined(__APPLE__) - for (const std::string &str : signals_in) { + for (const std::string& str : signals_in) { const int sig = SignalHandling::parseSignalNameToNumber(str); - SignalHandling::signal_conf_requests[SignalHandling::SIGNAL_REQUESTS_CHECKPOINT][sig] = true; - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(have_checkpoint_diagnostic, - "Signal handling was requested to checkpoint, but no checkpoint diagnostic is configured"); + SignalHandling::signal_conf_requests + [SignalHandling::SIGNAL_REQUESTS_CHECKPOINT][sig] = true; + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + have_checkpoint_diagnostic, + "Signal handling was requested to checkpoint, but no " + "checkpoint diagnostic is configured"); } #else - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(signals_in.empty(), - "Signal handling requested in input, but is not supported on this platform"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + signals_in.empty(), "Signal handling requested in input, but is " + "not supported on this platform"); #endif // set random seed std::string random_seed = "default"; pp_warpx.query("random_seed", random_seed); - if ( random_seed != "default" ) { + if (random_seed != "default") { const unsigned long myproc_1 = ParallelDescriptor::MyProc() + 1; - if ( random_seed == "random" ) { + if (random_seed == "random") { std::random_device rd; std::uniform_int_distribution dist(2, INT_MAX); const unsigned long cpu_seed = myproc_1 * dist(rd); const unsigned long gpu_seed = myproc_1 * dist(rd); ResetRandomSeed(cpu_seed, gpu_seed); - } else if ( std::stoi(random_seed) > 0 ) { + } else if (std::stoi(random_seed) > 0) { const unsigned long nprocs = ParallelDescriptor::NProcs(); const unsigned long seed_long = std::stoul(random_seed); const unsigned long cpu_seed = myproc_1 * seed_long; @@ -609,7 +624,8 @@ WarpX::ReadParameters () ResetRandomSeed(cpu_seed, gpu_seed); } else { WARPX_ABORT_WITH_MESSAGE( - "warpx.random_seed must be \"default\", \"random\" or an integer > 0."); + "warpx.random_seed must be \"default\", \"random\" or an " + "integer > 0."); } } @@ -618,50 +634,56 @@ WarpX::ReadParameters () utils::parser::queryWithParser(pp_warpx, "regrid_int", regrid_int); pp_warpx.query("do_subcycling", m_do_subcycling); pp_warpx.query("do_multi_J", do_multi_J); - if (do_multi_J) - { - utils::parser::getWithParser( - pp_warpx, "do_multi_J_n_depositions", do_multi_J_n_depositions); + if (do_multi_J) { + utils::parser::getWithParser(pp_warpx, "do_multi_J_n_depositions", + do_multi_J_n_depositions); } pp_warpx.query("use_hybrid_QED", use_hybrid_QED); pp_warpx.query("safe_guard_cells", m_safe_guard_cells); std::vector override_sync_intervals_string_vec = {"1"}; - pp_warpx.queryarr("override_sync_intervals", override_sync_intervals_string_vec); + pp_warpx.queryarr("override_sync_intervals", + override_sync_intervals_string_vec); override_sync_intervals = utils::parser::IntervalsParser(override_sync_intervals_string_vec); - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_do_subcycling != 1 || max_level <= 1, - "Subcycling method 1 only works for 2 levels."); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + m_do_subcycling != 1 || max_level <= 1, + "Subcycling method 1 only works for 2 levels."); ReadBoostedFrameParameters(gamma_boost, beta_boost, boost_direction); - // queryWithParser returns 1 if argument zmax_plasma_to_compute_max_step is - // specified by the user, 0 otherwise. - if(auto temp = 0.0_rt; utils::parser::queryWithParser(pp_warpx, "zmax_plasma_to_compute_max_step",temp)){ + // queryWithParser returns 1 if argument zmax_plasma_to_compute_max_step + // is specified by the user, 0 otherwise. + if (auto temp = 0.0_rt; utils::parser::queryWithParser( + pp_warpx, "zmax_plasma_to_compute_max_step", temp)) { m_zmax_plasma_to_compute_max_step = temp; } - pp_warpx.query("compute_max_step_from_btd", - compute_max_step_from_btd); + pp_warpx.query("compute_max_step_from_btd", compute_max_step_from_btd); if (do_moving_window) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( Geom(0).isPeriodic(moving_window_dir) == 0, - "The problem must be non-periodic in the moving window direction"); + "The problem must be non-periodic in the moving window " + "direction"); moving_window_x = geom[0].ProbLo(moving_window_dir); } m_p_ext_field_params = std::make_unique(pp_warpx); - if (m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::read_from_file || - m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::read_from_file){ - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(max_level == 0, - "External field reading is not implemented for more than one level"); + if (m_p_ext_field_params->B_ext_grid_type == + ExternalFieldType::read_from_file || + m_p_ext_field_params->E_ext_grid_type == + ExternalFieldType::read_from_file) { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + max_level == 0, "External field reading is not implemented for " + "more than one level"); } maxlevel_extEMfield_init = maxLevel(); pp_warpx.query("maxlevel_extEMfield_init", maxlevel_extEMfield_init); - pp_warpx.query_enum_sloppy("do_electrostatic", electrostatic_solver_id, "-_"); + pp_warpx.query_enum_sloppy("do_electrostatic", electrostatic_solver_id, + "-_"); // if an electrostatic solver is used, set the Maxwell solver to None if (electrostatic_solver_id != ElectrostaticSolverAlgo::None) { electromagnetic_solver_id = ElectromagneticSolverAlgo::None; @@ -670,36 +692,37 @@ WarpX::ReadParameters () pp_warpx.query_enum_sloppy("poisson_solver", poisson_solver_id, "-_"); #ifndef WARPX_DIM_3D WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - poisson_solver_id!=PoissonSolverAlgo::IntegratedGreenFunction, - "The FFT Poisson solver only works in 3D."); + poisson_solver_id != PoissonSolverAlgo::IntegratedGreenFunction, + "The FFT Poisson solver only works in 3D."); #endif #ifndef WARPX_USE_FFT WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - poisson_solver_id!=PoissonSolverAlgo::IntegratedGreenFunction, - "To use the FFT Poisson solver, compile with WARPX_USE_FFT=ON."); + poisson_solver_id != PoissonSolverAlgo::IntegratedGreenFunction, + "To use the FFT Poisson solver, compile with WARPX_USE_FFT=ON."); #endif - utils::parser::queryWithParser(pp_warpx, "self_fields_max_iters", magnetostatic_solver_max_iters); - utils::parser::queryWithParser(pp_warpx, "self_fields_verbosity", magnetostatic_solver_verbosity); + utils::parser::queryWithParser(pp_warpx, "self_fields_max_iters", + magnetostatic_solver_max_iters); + utils::parser::queryWithParser(pp_warpx, "self_fields_verbosity", + magnetostatic_solver_verbosity); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - ( - electrostatic_solver_id!=ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic || - poisson_solver_id!=PoissonSolverAlgo::IntegratedGreenFunction - ), - "The FFT Poisson solver is not implemented in labframe-electromagnetostatic mode yet." - ); + (electrostatic_solver_id != + ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic || + poisson_solver_id != PoissonSolverAlgo::IntegratedGreenFunction), + "The FFT Poisson solver is not implemented in " + "labframe-electromagnetostatic mode yet."); [[maybe_unused]] bool const eb_enabled = EB::enabled(); #if !defined(AMREX_USE_EB) WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - !eb_enabled, - "Embedded boundaries are requested via warpx.eb_enabled but were not compiled!" - ); + !eb_enabled, "Embedded boundaries are requested via " + "warpx.eb_enabled but were not compiled!"); #endif #ifdef WARPX_DIM_RZ const ParmParse pp_boundary("boundary"); - pp_boundary.query("verboncoeur_axis_correction", m_verboncoeur_axis_correction); + pp_boundary.query("verboncoeur_axis_correction", + m_verboncoeur_axis_correction); #endif // Read timestepping options @@ -709,24 +732,29 @@ WarpX::ReadParameters () pp_warpx.queryarr("dt_update_interval", dt_interval_vec); m_dt_update_interval = utils::parser::IntervalsParser(dt_interval_vec); - // Filter defaults to true for the explicit scheme, and false for the implicit schemes + // Filter defaults to true for the explicit scheme, and false for the + // implicit schemes if (evolve_scheme != EvolveScheme::Explicit) { use_filter = false; } - // Filter currently not working with FDTD solver in RZ geometry: turn OFF by default - // (see https://github.com/ECP-WarpX/WarpX/issues/1943) + // Filter currently not working with FDTD solver in RZ geometry: turn + // OFF by default (see https://github.com/ECP-WarpX/WarpX/issues/1943) #ifdef WARPX_DIM_RZ - if (WarpX::electromagnetic_solver_id != ElectromagneticSolverAlgo::PSATD) { WarpX::use_filter = false; } + if (WarpX::electromagnetic_solver_id != + ElectromagneticSolverAlgo::PSATD) { + WarpX::use_filter = false; + } #endif // Read filter and fill IntVect filter_npass_each_dir with // proper size for AMREX_SPACEDIM pp_warpx.query("use_filter", use_filter); pp_warpx.query("use_filter_compensation", use_filter_compensation); - Vector parse_filter_npass_each_dir(AMREX_SPACEDIM,1); - utils::parser::queryArrWithParser( - pp_warpx, "filter_npass_each_dir", parse_filter_npass_each_dir, 0, AMREX_SPACEDIM); + Vector parse_filter_npass_each_dir(AMREX_SPACEDIM, 1); + utils::parser::queryArrWithParser(pp_warpx, "filter_npass_each_dir", + parse_filter_npass_each_dir, 0, + AMREX_SPACEDIM); filter_npass_each_dir[0] = parse_filter_npass_each_dir[0]; #if (AMREX_SPACEDIM >= 2) filter_npass_each_dir[1] = parse_filter_npass_each_dir[1]; @@ -735,45 +763,50 @@ WarpX::ReadParameters () filter_npass_each_dir[2] = parse_filter_npass_each_dir[2]; #endif - // TODO When k-space filtering will be implemented also for Cartesian geometries, - // this code block will have to be applied in all cases (remove #ifdef condition) + // TODO When k-space filtering will be implemented also for Cartesian + // geometries, this code block will have to be applied in all cases + // (remove #ifdef condition) #ifdef WARPX_DIM_RZ - if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { + if (WarpX::electromagnetic_solver_id == + ElectromagneticSolverAlgo::PSATD) { // With RZ spectral, only use k-space filtering use_kspace_filter = use_filter; use_filter = false; - } - else - { - if (WarpX::electromagnetic_solver_id != ElectromagneticSolverAlgo::HybridPIC) { - // Filter currently not working with FDTD solver in RZ geometry along R - // (see https://github.com/ECP-WarpX/WarpX/issues/1943) - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(!use_filter || filter_npass_each_dir[0] == 0, - "In RZ geometry with FDTD, filtering can only be applied along z. This can be controlled by setting warpx.filter_npass_each_dir"); + } else { + if (WarpX::electromagnetic_solver_id != + ElectromagneticSolverAlgo::HybridPIC) { + // Filter currently not working with FDTD solver in RZ geometry + // along R (see https://github.com/ECP-WarpX/WarpX/issues/1943) + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + !use_filter || filter_npass_each_dir[0] == 0, + "In RZ geometry with FDTD, filtering can only be applied " + "along z. This can be controlled by setting " + "warpx.filter_npass_each_dir"); } else { if (use_filter && filter_npass_each_dir[0] > 0) { ablastr::warn_manager::WMRecordWarning( "HybridPIC ElectromagneticSolver", - "Radial Filtering in RZ is not currently using radial geometric weighting to conserve charge. Use at your own risk.", - ablastr::warn_manager::WarnPriority::low - ); + "Radial Filtering in RZ is not currently using radial " + "geometric weighting to conserve charge. Use at your " + "own risk.", + ablastr::warn_manager::WarnPriority::low); } } } #endif - utils::parser::queryWithParser( - pp_warpx, "num_mirrors", m_num_mirrors); - if (m_num_mirrors>0){ + utils::parser::queryWithParser(pp_warpx, "num_mirrors", m_num_mirrors); + if (m_num_mirrors > 0) { m_mirror_z.resize(m_num_mirrors); - utils::parser::getArrWithParser( - pp_warpx, "mirror_z", m_mirror_z, 0, m_num_mirrors); + utils::parser::getArrWithParser(pp_warpx, "mirror_z", m_mirror_z, 0, + m_num_mirrors); m_mirror_z_width.resize(m_num_mirrors); - utils::parser::getArrWithParser( - pp_warpx, "mirror_z_width", m_mirror_z_width, 0, m_num_mirrors); + utils::parser::getArrWithParser(pp_warpx, "mirror_z_width", + m_mirror_z_width, 0, m_num_mirrors); m_mirror_z_npoints.resize(m_num_mirrors); - utils::parser::getArrWithParser( - pp_warpx, "mirror_z_npoints", m_mirror_z_npoints, 0, m_num_mirrors); + utils::parser::getArrWithParser(pp_warpx, "mirror_z_npoints", + m_mirror_z_npoints, 0, + m_num_mirrors); } pp_warpx.query("do_single_precision_comms", do_single_precision_comms); @@ -782,73 +815,94 @@ WarpX::ReadParameters () do_single_precision_comms = false; ablastr::warn_manager::WMRecordWarning( "comms", - "Overwrote warpx.do_single_precision_comms to be 0, since WarpX was built in single precision.", + "Overwrote warpx.do_single_precision_comms to be 0, since " + "WarpX was built in single precision.", ablastr::warn_manager::WarnPriority::low); } #endif - pp_warpx.query("do_shared_mem_charge_deposition", do_shared_mem_charge_deposition); - pp_warpx.query("do_shared_mem_current_deposition", do_shared_mem_current_deposition); + pp_warpx.query("do_shared_mem_charge_deposition", + do_shared_mem_charge_deposition); + pp_warpx.query("do_shared_mem_current_deposition", + do_shared_mem_current_deposition); #if !(defined(AMREX_USE_HIP) || defined(AMREX_USE_CUDA)) - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(!do_shared_mem_current_deposition, - "requested shared memory for current deposition, but shared memory is only available for CUDA or HIP"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + !do_shared_mem_current_deposition, + "requested shared memory for current deposition, but shared memory " + "is only available for CUDA or HIP"); #endif pp_warpx.query("shared_mem_current_tpb", shared_mem_current_tpb); // initialize the shared tilesize Vector vect_shared_tilesize(AMREX_SPACEDIM, 1); - const bool shared_tilesize_is_specified = utils::parser::queryArrWithParser(pp_warpx, "shared_tilesize", - vect_shared_tilesize, 0, AMREX_SPACEDIM); - if (shared_tilesize_is_specified){ - for (int i=0; i(quantum_xi * PhysConst::c * PhysConst::c); + m_quantum_xi_c2 = static_cast( + quantum_xi * PhysConst::c * PhysConst::c); } const auto at_least_one_boundary_is_pml = - (std::any_of(WarpX::field_boundary_lo.begin(), WarpX::field_boundary_lo.end(), - [](const auto& cc){return cc == FieldBoundaryType::PML;}) - || - std::any_of(WarpX::field_boundary_hi.begin(), WarpX::field_boundary_hi.end(), - [](const auto& cc){return cc == FieldBoundaryType::PML;}) - ); + (std::any_of(WarpX::field_boundary_lo.begin(), + WarpX::field_boundary_lo.end(), + [] (const auto& cc) { + return cc == FieldBoundaryType::PML; + }) || + std::any_of(WarpX::field_boundary_hi.begin(), + WarpX::field_boundary_hi.end(), [] (const auto& cc) { + return cc == FieldBoundaryType::PML; + })); const auto at_least_one_boundary_is_silver_mueller = - (std::any_of(WarpX::field_boundary_lo.begin(), WarpX::field_boundary_lo.end(), - [](const auto& cc){return cc == FieldBoundaryType::Absorbing_SilverMueller;}) - || - std::any_of(WarpX::field_boundary_hi.begin(), WarpX::field_boundary_hi.end(), - [](const auto& cc){return cc == FieldBoundaryType::Absorbing_SilverMueller;}) - ); + (std::any_of(WarpX::field_boundary_lo.begin(), + WarpX::field_boundary_lo.end(), + [] (const auto& cc) { + return cc == + FieldBoundaryType::Absorbing_SilverMueller; + }) || + std::any_of(WarpX::field_boundary_hi.begin(), + WarpX::field_boundary_hi.end(), [] (const auto& cc) { + return cc == + FieldBoundaryType::Absorbing_SilverMueller; + })); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - !(at_least_one_boundary_is_pml && at_least_one_boundary_is_silver_mueller), - "PML and Silver-Mueller boundary conditions cannot be activated at the same time."); + !(at_least_one_boundary_is_pml && + at_least_one_boundary_is_silver_mueller), + "PML and Silver-Mueller boundary conditions cannot be activated at " + "the same time."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( (!at_least_one_boundary_is_silver_mueller) || - (electromagnetic_solver_id == ElectromagneticSolverAlgo::Yee), - "The Silver-Mueller boundary condition can only be used with the Yee solver."); + (electromagnetic_solver_id == ElectromagneticSolverAlgo::Yee), + "The Silver-Mueller boundary condition can only be used with the " + "Yee solver."); utils::parser::queryWithParser(pp_warpx, "pml_ncell", pml_ncell); utils::parser::queryWithParser(pp_warpx, "pml_delta", pml_delta); @@ -858,31 +912,42 @@ WarpX::ReadParameters () pp_warpx.query("do_similar_dm_pml", do_similar_dm_pml); // Read `v_particle_pml` in units of the speed of light v_particle_pml = 1._rt; - utils::parser::queryWithParser(pp_warpx, "v_particle_pml", v_particle_pml); - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(0._rt < v_particle_pml && v_particle_pml <= 1._rt, - "Input value for the velocity warpx.v_particle_pml of the macroparticle must be in (0,1] (in units of c)."); + utils::parser::queryWithParser(pp_warpx, "v_particle_pml", + v_particle_pml); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + 0._rt < v_particle_pml && v_particle_pml <= 1._rt, + "Input value for the velocity warpx.v_particle_pml of the " + "macroparticle must be in (0,1] (in units of c)."); // Scale by the speed of light v_particle_pml = v_particle_pml * PhysConst::c; - // Default values of WarpX::do_pml_dive_cleaning and WarpX::do_pml_divb_cleaning: - // true for Cartesian PSATD solver, false otherwise + // Default values of WarpX::do_pml_dive_cleaning and + // WarpX::do_pml_divb_cleaning: true for Cartesian PSATD solver, false + // otherwise do_pml_dive_cleaning = false; do_pml_divb_cleaning = false; #ifndef WARPX_DIM_RZ - if (electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) - { + if (electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { do_pml_dive_cleaning = true; do_pml_divb_cleaning = true; } - // If WarpX::do_dive_cleaning = true, set also WarpX::do_pml_dive_cleaning = true - // (possibly overwritten by users in the input file, see query below) - if (do_dive_cleaning) { do_pml_dive_cleaning = true; } + // If WarpX::do_dive_cleaning = true, set also + // WarpX::do_pml_dive_cleaning = true (possibly overwritten by users in + // the input file, see query below) + if (do_dive_cleaning) { + do_pml_dive_cleaning = true; + } - // If WarpX::do_divb_cleaning = true, set also WarpX::do_pml_divb_cleaning = true - // (possibly overwritten by users in the input file, see query below) - // TODO Implement div(B) cleaning in PML with FDTD and remove second if condition - if (do_divb_cleaning && electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { do_pml_divb_cleaning = true; } + // If WarpX::do_divb_cleaning = true, set also + // WarpX::do_pml_divb_cleaning = true (possibly overwritten by users in + // the input file, see query below) + // TODO Implement div(B) cleaning in PML with FDTD and remove second if + // condition + if (do_divb_cleaning && + electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { + do_pml_divb_cleaning = true; + } #endif // Query input parameters to use div(E) and div(B) cleaning in PMLs @@ -890,43 +955,52 @@ WarpX::ReadParameters () pp_warpx.query("do_pml_divb_cleaning", do_pml_divb_cleaning); // TODO Implement div(B) cleaning in PML with FDTD and remove ASSERT - if (electromagnetic_solver_id != ElectromagneticSolverAlgo::PSATD) - { + if (electromagnetic_solver_id != ElectromagneticSolverAlgo::PSATD) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( do_pml_divb_cleaning == false, - "warpx.do_pml_divb_cleaning = true not implemented for FDTD solver"); + "warpx.do_pml_divb_cleaning = true not implemented for FDTD " + "solver"); } // Divergence cleaning in PMLs for PSATD solver implemented only // for both div(E) and div(B) cleaning - if (electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) - { + if (electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( do_pml_dive_cleaning == do_pml_divb_cleaning, - "warpx.do_pml_dive_cleaning = " - + std::to_string(do_pml_dive_cleaning) - +" and warpx.do_pml_divb_cleaning = " - + std::to_string(do_pml_divb_cleaning) - + ": this case is not implemented yet," - + " please set both parameters to the same value" - ); + "warpx.do_pml_dive_cleaning = " + + std::to_string(do_pml_dive_cleaning) + + " and warpx.do_pml_divb_cleaning = " + + std::to_string(do_pml_divb_cleaning) + + ": this case is not implemented yet," + + " please set both parameters to the same value"); } #ifdef WARPX_DIM_RZ - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( ::isAnyBoundaryPML(field_boundary_lo, field_boundary_hi) == false || electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD, - "PML are not implemented in RZ geometry with FDTD; please set a different boundary condition using boundary.field_lo and boundary.field_hi."); - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( field_boundary_lo[1] != FieldBoundaryType::PML && field_boundary_hi[1] != FieldBoundaryType::PML, - "PML are not implemented in RZ geometry along z; please set a different boundary condition using boundary.field_lo and boundary.field_hi."); - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( (do_pml_dive_cleaning == false && do_pml_divb_cleaning == false), - "do_pml_dive_cleaning and do_pml_divb_cleaning are not implemented in RZ geometry." ); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + ::isAnyBoundaryPML(field_boundary_lo, field_boundary_hi) == false || + electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD, + "PML are not implemented in RZ geometry with FDTD; please set a " + "different boundary condition using boundary.field_lo and " + "boundary.field_hi."); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + field_boundary_lo[1] != FieldBoundaryType::PML && + field_boundary_hi[1] != FieldBoundaryType::PML, + "PML are not implemented in RZ geometry along z; please set a " + "different boundary condition using boundary.field_lo and " + "boundary.field_hi."); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + (do_pml_dive_cleaning == false && do_pml_divb_cleaning == false), + "do_pml_dive_cleaning and do_pml_divb_cleaning are not implemented " + "in RZ geometry."); #endif WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - (do_pml_j_damping==0)||(do_pml_in_domain==1), - "J-damping can only be done when PML are inside simulation domain (do_pml_in_domain=1)" - ); + (do_pml_j_damping == 0) || (do_pml_in_domain == 1), + "J-damping can only be done when PML are inside simulation domain " + "(do_pml_in_domain=1)"); - pp_warpx.query("synchronize_velocity_for_diagnostics", synchronize_velocity_for_diagnostics); + pp_warpx.query("synchronize_velocity_for_diagnostics", + synchronize_velocity_for_diagnostics); { // Parameters below control all plotfile diagnostics @@ -942,29 +1016,38 @@ WarpX::ReadParameters () ParmParse pp_vismf("vismf"); pp_vismf.add("usesingleread", use_single_read); pp_vismf.add("usesinglewrite", use_single_write); - utils::parser::queryWithParser(pp_warpx, "mffile_nstreams", mffile_nstreams); + utils::parser::queryWithParser(pp_warpx, "mffile_nstreams", + mffile_nstreams); VisMF::SetMFFileInStreams(mffile_nstreams); - utils::parser::queryWithParser(pp_warpx, "field_io_nfiles", field_io_nfiles); + utils::parser::queryWithParser(pp_warpx, "field_io_nfiles", + field_io_nfiles); VisMF::SetNOutFiles(field_io_nfiles); - utils::parser::queryWithParser(pp_warpx, "particle_io_nfiles", particle_io_nfiles); + utils::parser::queryWithParser(pp_warpx, "particle_io_nfiles", + particle_io_nfiles); ParmParse pp_particles("particles"); pp_particles.add("particles_nfiles", particle_io_nfiles); } if (maxLevel() > 0) { Vector lo, hi; - const bool fine_tag_lo_specified = utils::parser::queryArrWithParser(pp_warpx, "fine_tag_lo", lo); - const bool fine_tag_hi_specified = utils::parser::queryArrWithParser(pp_warpx, "fine_tag_hi", hi); + const bool fine_tag_lo_specified = + utils::parser::queryArrWithParser(pp_warpx, "fine_tag_lo", lo); + const bool fine_tag_hi_specified = + utils::parser::queryArrWithParser(pp_warpx, "fine_tag_hi", hi); std::string ref_patch_function; - const bool parser_specified = pp_warpx.query("ref_patch_function(x,y,z)",ref_patch_function); - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( ((fine_tag_lo_specified && fine_tag_hi_specified) || - parser_specified ), - "For max_level > 0, you need to either set\ + const bool parser_specified = + pp_warpx.query("ref_patch_function(x,y,z)", ref_patch_function); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + ((fine_tag_lo_specified && fine_tag_hi_specified) || + parser_specified), + "For max_level > 0, you need to either set\ warpx.fine_tag_lo and warpx.fine_tag_hi\ or warpx.ref_patch_function(x,y,z)"); - if ( (fine_tag_lo_specified && fine_tag_hi_specified) && parser_specified) { - ablastr::warn_manager::WMRecordWarning("Refined patch", "Both fine_tag_lo,fine_tag_hi\ + if ((fine_tag_lo_specified && fine_tag_hi_specified) && + parser_specified) { + ablastr::warn_manager::WMRecordWarning( + "Refined patch", "Both fine_tag_lo,fine_tag_hi\ and ref_patch_function(x,y,z) are provided. Note that fine_tag_lo/fine_tag_hi will\ override the ref_patch_function(x,y,z) for defining the refinement patches"); } @@ -972,9 +1055,11 @@ WarpX::ReadParameters () fine_tag_lo = RealVect{lo}; fine_tag_hi = RealVect{hi}; } else { - utils::parser::Store_parserString(pp_warpx, "ref_patch_function(x,y,z)", ref_patch_function); - ref_patch_parser = std::make_unique( - utils::parser::makeParser(ref_patch_function,{"x","y","z"})); + utils::parser::Store_parserString( + pp_warpx, "ref_patch_function(x,y,z)", ref_patch_function); + ref_patch_parser = + std::make_unique(utils::parser::makeParser( + ref_patch_function, {"x", "y", "z"})); } } @@ -985,13 +1070,18 @@ WarpX::ReadParameters () pp_warpx.query_enum_sloppy("grid_type", grid_type, "-_"); // Use same shape factors in all directions, for gathering - if (grid_type == GridType::Collocated) { galerkin_interpolation = false; } + if (grid_type == GridType::Collocated) { + galerkin_interpolation = false; + } #ifdef WARPX_DIM_RZ // Only needs to be set with WARPX_DIM_RZ, otherwise defaults to 1 - utils::parser::queryWithParser(pp_warpx, "n_rz_azimuthal_modes", n_rz_azimuthal_modes); - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( n_rz_azimuthal_modes > 0, - "The number of azimuthal modes (n_rz_azimuthal_modes) must be at least 1"); + utils::parser::queryWithParser(pp_warpx, "n_rz_azimuthal_modes", + n_rz_azimuthal_modes); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + n_rz_azimuthal_modes > 0, + "The number of azimuthal modes (n_rz_azimuthal_modes) must be at " + "least 1"); #endif // Check whether fluid species will be used @@ -1001,21 +1091,25 @@ WarpX::ReadParameters () pp_fluids.queryarr("species_names", fluid_species_names); do_fluid_species = !fluid_species_names.empty(); if (do_fluid_species) { - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(max_level <= 1, - "Fluid species cannot currently be used with mesh refinement."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - electrostatic_solver_id != ElectrostaticSolverAlgo::Relativistic, - "Fluid species cannot currently be used with the relativistic electrostatic solver."); + max_level <= 1, "Fluid species cannot currently be used " + "with mesh refinement."); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + electrostatic_solver_id != + ElectrostaticSolverAlgo::Relativistic, + "Fluid species cannot currently be used with the " + "relativistic electrostatic solver."); #ifdef WARPX_DIM_RZ - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( n_rz_azimuthal_modes <= 1, - "Fluid species cannot be used with more than 1 azimuthal mode."); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + n_rz_azimuthal_modes <= 1, + "Fluid species cannot be used with more than 1 azimuthal " + "mode."); #endif } } // Set default parameters with hybrid grid (parsed later below) - if (grid_type == GridType::Hybrid) - { + if (grid_type == GridType::Hybrid) { // Finite-order centering of fields (staggered to nodal) // Default field gathering algorithm will be set below field_centering_nox = 8; @@ -1034,21 +1128,28 @@ WarpX::ReadParameters () "warpx.grid_type=hybrid is not implemented in RZ geometry"); #endif - // Update default to external projection divb cleaner if external fields are loaded, - // the grids are staggered, and the solver is compatible with the cleaner - if (!do_divb_cleaning - && m_p_ext_field_params->B_ext_grid_type != ExternalFieldType::default_zero - && m_p_ext_field_params->B_ext_grid_type != ExternalFieldType::constant - && grid_type != GridType::Collocated - && (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::Yee - || WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC - || ( (WarpX::electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame - || WarpX::electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) - && WarpX::poisson_solver_id == PoissonSolverAlgo::Multigrid))) - { + // Update default to external projection divb cleaner if external fields + // are loaded, the grids are staggered, and the solver is compatible + // with the cleaner + if (!do_divb_cleaning && + m_p_ext_field_params->B_ext_grid_type != + ExternalFieldType::default_zero && + m_p_ext_field_params->B_ext_grid_type != + ExternalFieldType::constant && + grid_type != GridType::Collocated && + (WarpX::electromagnetic_solver_id == + ElectromagneticSolverAlgo::Yee || + WarpX::electromagnetic_solver_id == + ElectromagneticSolverAlgo::HybridPIC || + ((WarpX::electrostatic_solver_id == + ElectrostaticSolverAlgo::LabFrame || + WarpX::electrostatic_solver_id == + ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) && + WarpX::poisson_solver_id == PoissonSolverAlgo::Multigrid))) { m_do_divb_cleaning_external = true; } - pp_warpx.query("do_divb_cleaning_external", m_do_divb_cleaning_external); + pp_warpx.query("do_divb_cleaning_external", + m_do_divb_cleaning_external); // If true, the current is deposited on a nodal grid and centered onto // a staggered grid. Setting warpx.do_current_centering=1 makes sense @@ -1056,54 +1157,58 @@ WarpX::ReadParameters () // warpx.grid_type=staggered, Maxwell's equations are solved either on a // collocated grid or on a staggered grid without current centering. pp_warpx.query("do_current_centering", do_current_centering); - if (do_current_centering) - { + if (do_current_centering) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( grid_type == GridType::Hybrid, - "warpx.do_current_centering=1 can be used only with warpx.grid_type=hybrid"); + "warpx.do_current_centering=1 can be used only with " + "warpx.grid_type=hybrid"); - utils::parser::queryWithParser( - pp_warpx, "current_centering_nox", current_centering_nox); - utils::parser::queryWithParser( - pp_warpx, "current_centering_noy", current_centering_noy); - utils::parser::queryWithParser( - pp_warpx, "current_centering_noz", current_centering_noz); + utils::parser::queryWithParser(pp_warpx, "current_centering_nox", + current_centering_nox); + utils::parser::queryWithParser(pp_warpx, "current_centering_noy", + current_centering_noy); + utils::parser::queryWithParser(pp_warpx, "current_centering_noz", + current_centering_noz); - AllocateCenteringCoefficients(device_current_centering_stencil_coeffs_x, - device_current_centering_stencil_coeffs_y, - device_current_centering_stencil_coeffs_z, - current_centering_nox, - current_centering_noy, - current_centering_noz, - grid_type); + AllocateCenteringCoefficients( + device_current_centering_stencil_coeffs_x, + device_current_centering_stencil_coeffs_y, + device_current_centering_stencil_coeffs_z, + current_centering_nox, current_centering_noy, + current_centering_noz, grid_type); } WARPX_ALWAYS_ASSERT_WITH_MESSAGE( maxLevel() == 0 || !do_current_centering, - "Finite-order centering of currents is not implemented with mesh refinement" - ); + "Finite-order centering of currents is not implemented with mesh " + "refinement"); } { const ParmParse pp_algo("algo"); #ifdef WARPX_DIM_RZ - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( electromagnetic_solver_id != ElectromagneticSolverAlgo::CKC, + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + electromagnetic_solver_id != ElectromagneticSolverAlgo::CKC, "algo.maxwell_solver = ckc is not (yet) available for RZ geometry"); #endif #ifndef WARPX_USE_FFT - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( electromagnetic_solver_id != ElectromagneticSolverAlgo::PSATD, - "algo.maxwell_solver = psatd is not supported because WarpX was built without spectral solvers"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + electromagnetic_solver_id != ElectromagneticSolverAlgo::PSATD, + "algo.maxwell_solver = psatd is not supported because WarpX was " + "built without spectral solvers"); #endif #if defined(WARPX_DIM_1D_Z) && defined(WARPX_USE_FFT) - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( electromagnetic_solver_id != ElectromagneticSolverAlgo::PSATD, + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + electromagnetic_solver_id != ElectromagneticSolverAlgo::PSATD, "algo.maxwell_solver = psatd is not available for 1D geometry"); #endif #ifdef WARPX_DIM_RZ - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(Geom(0).isPeriodic(0) == 0, + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + Geom(0).isPeriodic(0) == 0, "The problem must not be periodic in the radial direction"); // Ensure code aborts if PEC is specified at r=0 for RZ - if (Geom(0).ProbLo(0) == 0){ + if (Geom(0).ProbLo(0) == 0) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( WarpX::field_boundary_lo[0] == FieldBoundaryType::None, "Error : Field boundary at r=0 must be ``none``. \n"); @@ -1111,10 +1216,10 @@ WarpX::ReadParameters () const ParmParse pp_boundary("boundary"); if (pp_boundary.contains("particle_lo")) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - WarpX::particle_boundary_lo[0] == ParticleBoundaryType::None, + WarpX::particle_boundary_lo[0] == + ParticleBoundaryType::None, "Error : Particle boundary at r=0 must be ``none``. \n"); } - } if (electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { @@ -1125,56 +1230,63 @@ WarpX::ReadParameters () } #endif - // note: current_deposition must be set after maxwell_solver (electromagnetic_solver_id) or - // do_electrostatic (electrostatic_solver_id) are already determined, - // because its default depends on the solver selection + // note: current_deposition must be set after maxwell_solver + // (electromagnetic_solver_id) or + // do_electrostatic (electrostatic_solver_id) are already + // determined, because its default depends on the solver selection if (electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD || electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC || electrostatic_solver_id != ElectrostaticSolverAlgo::None) { current_deposition_algo = CurrentDepositionAlgo::Direct; } - pp_algo.query_enum_sloppy("current_deposition", current_deposition_algo, "-_"); - pp_algo.query_enum_sloppy("charge_deposition", charge_deposition_algo, "-_"); - pp_algo.query_enum_sloppy("particle_pusher", particle_pusher_algo, "-_"); + pp_algo.query_enum_sloppy("current_deposition", current_deposition_algo, + "-_"); + pp_algo.query_enum_sloppy("charge_deposition", charge_deposition_algo, + "-_"); + pp_algo.query_enum_sloppy("particle_pusher", particle_pusher_algo, + "-_"); // check for implicit evolve scheme if (evolve_scheme == EvolveScheme::SemiImplicitEM) { m_implicit_solver = std::make_unique(); - } - else if (evolve_scheme == EvolveScheme::ThetaImplicitEM) { + } else if (evolve_scheme == EvolveScheme::ThetaImplicitEM) { m_implicit_solver = std::make_unique(); - } - else if (evolve_scheme == EvolveScheme::StrangImplicitSpectralEM) { + } else if (evolve_scheme == EvolveScheme::StrangImplicitSpectralEM) { m_implicit_solver = std::make_unique(); } // implicit evolve schemes not setup to use mirrors if (evolve_scheme == EvolveScheme::SemiImplicitEM || evolve_scheme == EvolveScheme::ThetaImplicitEM) { - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( m_num_mirrors == 0, + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + m_num_mirrors == 0, "Mirrors cannot be used with Implicit evolve schemes."); } WARPX_ALWAYS_ASSERT_WITH_MESSAGE( current_deposition_algo != CurrentDepositionAlgo::Esirkepov || - !do_current_centering, - "Current centering (nodal deposition) cannot be used with Esirkepov deposition." - "Please set warpx.do_current_centering = 0 or algo.current_deposition = direct."); + !do_current_centering, + "Current centering (nodal deposition) cannot be used with " + "Esirkepov deposition." + "Please set warpx.do_current_centering = 0 or " + "algo.current_deposition = direct."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( current_deposition_algo != CurrentDepositionAlgo::Villasenor || - !do_current_centering, - "Current centering (nodal deposition) cannot be used with Villasenor deposition." - "Please set warpx.do_current_centering = 0 or algo.current_deposition = direct."); + !do_current_centering, + "Current centering (nodal deposition) cannot be used with " + "Villasenor deposition." + "Please set warpx.do_current_centering = 0 or " + "algo.current_deposition = direct."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( WarpX::current_deposition_algo != CurrentDepositionAlgo::Vay || - !do_current_centering, + !do_current_centering, "Vay deposition not implemented with current centering"); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( WarpX::current_deposition_algo != CurrentDepositionAlgo::Vay || - maxLevel() <= 0, + maxLevel() <= 0, "Vay deposition not implemented with mesh refinement"); if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) { @@ -1190,31 +1302,31 @@ WarpX::ReadParameters () } // Query algo.field_gathering from input, set field_gathering_algo to - // "default" if not found (default defined in Utils/WarpXAlgorithmSelection.cpp) - pp_algo.query_enum_sloppy("field_gathering", field_gathering_algo, "-_"); + // "default" if not found (default defined in + // Utils/WarpXAlgorithmSelection.cpp) + pp_algo.query_enum_sloppy("field_gathering", field_gathering_algo, + "-_"); - // Set default field gathering algorithm for hybrid grids (momentum-conserving) + // Set default field gathering algorithm for hybrid grids + // (momentum-conserving) std::string tmp_algo; // - algo.field_gathering not found in the input // - field_gathering_algo set to "default" above // (default defined in Utils/WarpXAlgorithmSelection.cpp) // - reset default value here for hybrid grids - if (!pp_algo.query("field_gathering", tmp_algo)) - { - if (grid_type == GridType::Hybrid) - { + if (!pp_algo.query("field_gathering", tmp_algo)) { + if (grid_type == GridType::Hybrid) { field_gathering_algo = GatheringAlgo::MomentumConserving; } } // - algo.field_gathering found in the input // - field_gathering_algo read above and set to user-defined value - else - { - if (grid_type == GridType::Hybrid) - { + else { + if (grid_type == GridType::Hybrid) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( field_gathering_algo == GatheringAlgo::MomentumConserving, - "Hybrid grid (warpx.grid_type=hybrid) should be used only with " + "Hybrid grid (warpx.grid_type=hybrid) should be used only " + "with " "momentum-conserving field gathering algorithm " "(algo.field_gathering=momentum-conserving)"); } @@ -1222,10 +1334,12 @@ WarpX::ReadParameters () // Use same shape factors in all directions // - with momentum-conserving field gathering - if (field_gathering_algo == GatheringAlgo::MomentumConserving) {galerkin_interpolation = false;} + if (field_gathering_algo == GatheringAlgo::MomentumConserving) { + galerkin_interpolation = false; + } // - with direct current deposition and the EM solver - if( electromagnetic_solver_id != ElectromagneticSolverAlgo::None && - electromagnetic_solver_id != ElectromagneticSolverAlgo::HybridPIC ) { + if (electromagnetic_solver_id != ElectromagneticSolverAlgo::None && + electromagnetic_solver_id != ElectromagneticSolverAlgo::HybridPIC) { if (current_deposition_algo == CurrentDepositionAlgo::Direct) { galerkin_interpolation = false; } @@ -1233,7 +1347,7 @@ WarpX::ReadParameters () { const ParmParse pp_interpolation("interpolation"); - pp_interpolation.query("galerkin_scheme",galerkin_interpolation); + pp_interpolation.query("galerkin_scheme", galerkin_interpolation); } // With the PSATD solver, momentum-conserving field gathering @@ -1241,15 +1355,14 @@ WarpX::ReadParameters () // TODO Needs debugging if (electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD && field_gathering_algo == GatheringAlgo::MomentumConserving && - maxLevel() > 0) - { + maxLevel() > 0) { WARPX_ABORT_WITH_MESSAGE( "With the PSATD solver, momentum-conserving field gathering" " combined with mesh refinement is currently not implemented"); } pp_algo.query_enum_sloppy("em_solver_medium", m_em_solver_medium, "-_"); - if (m_em_solver_medium == MediumForEM::Macroscopic ) { + if (m_em_solver_medium == MediumForEM::Macroscopic) { pp_algo.query_enum_sloppy("macroscopic_sigma_method", macroscopic_solver_algo, "-_"); } @@ -1260,49 +1373,64 @@ WarpX::ReadParameters () WARPX_ALWAYS_ASSERT_WITH_MESSAGE( current_deposition_algo == CurrentDepositionAlgo::Esirkepov || - current_deposition_algo == CurrentDepositionAlgo::Villasenor || - current_deposition_algo == CurrentDepositionAlgo::Direct, - "Only Esirkepov, Villasenor, or Direct current deposition supported with the implicit and semi-implicit schemes"); + current_deposition_algo == + CurrentDepositionAlgo::Villasenor || + current_deposition_algo == CurrentDepositionAlgo::Direct, + "Only Esirkepov, Villasenor, or Direct current deposition " + "supported with the implicit and semi-implicit schemes"); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( electromagnetic_solver_id == ElectromagneticSolverAlgo::Yee || - electromagnetic_solver_id == ElectromagneticSolverAlgo::CKC || - electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD, - "Only the Yee, CKC, and PSATD EM solvers are supported with the implicit and semi-implicit schemes"); + electromagnetic_solver_id == + ElectromagneticSolverAlgo::CKC || + electromagnetic_solver_id == + ElectromagneticSolverAlgo::PSATD, + "Only the Yee, CKC, and PSATD EM solvers are supported with " + "the implicit and semi-implicit schemes"); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( particle_pusher_algo == ParticlePusherAlgo::Boris || - particle_pusher_algo == ParticlePusherAlgo::HigueraCary, - "Only the Boris and Higuera particle pushers are supported with the implicit and semi-implicit schemes"); + particle_pusher_algo == ParticlePusherAlgo::HigueraCary, + "Only the Boris and Higuera particle pushers are supported " + "with the implicit and semi-implicit schemes"); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( field_gathering_algo != GatheringAlgo::MomentumConserving, - "With implicit and semi-implicit schemes, the momentum conserving field gather is not supported as it would not conserve energy"); + "With implicit and semi-implicit schemes, the momentum " + "conserving field gather is not supported as it would not " + "conserve energy"); } if (evolve_scheme == EvolveScheme::StrangImplicitSpectralEM) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD, - "With the strang_implicit_spectral_em evolve scheme, the algo.maxwell_solver must be psatd"); + "With the strang_implicit_spectral_em evolve scheme, the " + "algo.maxwell_solver must be psatd"); } // Load balancing parameters std::vector load_balance_intervals_string_vec = {"0"}; - pp_algo.queryarr("load_balance_intervals", load_balance_intervals_string_vec); - load_balance_intervals = utils::parser::IntervalsParser( - load_balance_intervals_string_vec); + pp_algo.queryarr("load_balance_intervals", + load_balance_intervals_string_vec); + load_balance_intervals = + utils::parser::IntervalsParser(load_balance_intervals_string_vec); pp_algo.query("load_balance_with_sfc", load_balance_with_sfc); // Knapsack factor only used with non-SFC strategy if (!load_balance_with_sfc) { - pp_algo.query("load_balance_knapsack_factor", load_balance_knapsack_factor); + pp_algo.query("load_balance_knapsack_factor", + load_balance_knapsack_factor); } - utils::parser::queryWithParser(pp_algo, "load_balance_efficiency_ratio_threshold", - load_balance_efficiency_ratio_threshold); - pp_algo.query_enum_sloppy("load_balance_costs_update", load_balance_costs_update_algo, "-_"); - if (WarpX::load_balance_costs_update_algo==LoadBalanceCostsUpdateAlgo::Heuristic) { - utils::parser::queryWithParser( - pp_algo, "costs_heuristic_cells_wt", costs_heuristic_cells_wt); - utils::parser::queryWithParser( - pp_algo, "costs_heuristic_particles_wt", costs_heuristic_particles_wt); + utils::parser::queryWithParser( + pp_algo, "load_balance_efficiency_ratio_threshold", + load_balance_efficiency_ratio_threshold); + pp_algo.query_enum_sloppy("load_balance_costs_update", + load_balance_costs_update_algo, "-_"); + if (WarpX::load_balance_costs_update_algo == + LoadBalanceCostsUpdateAlgo::Heuristic) { + utils::parser::queryWithParser(pp_algo, "costs_heuristic_cells_wt", + costs_heuristic_cells_wt); + utils::parser::queryWithParser(pp_algo, + "costs_heuristic_particles_wt", + costs_heuristic_particles_wt); } // Parse algo.particle_shape and check that input is acceptable @@ -1322,40 +1450,44 @@ WarpX::ReadParameters () // representation of a laser pulse in cylindrical coordinates // requires at least 2 azimuthal modes if (!lasers_names.empty() && n_rz_azimuthal_modes < 2) { - ablastr::warn_manager::WMRecordWarning("Laser", - "Laser pulse representation in RZ requires at least 2 azimuthal modes", - ablastr::warn_manager::WarnPriority::high); + ablastr::warn_manager::WMRecordWarning( + "Laser", + "Laser pulse representation in RZ requires at least 2 " + "azimuthal modes", + ablastr::warn_manager::WarnPriority::high); } #endif std::vector sort_intervals_string_vec = {"-1"}; int particle_shape; if (!species_names.empty() || !lasers_names.empty()) { - if (utils::parser::queryWithParser(pp_algo, "particle_shape", particle_shape)){ + if (utils::parser::queryWithParser(pp_algo, "particle_shape", + particle_shape)) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - (particle_shape >= 1) && (particle_shape <=4), - "algo.particle_shape can be only 1, 2, 3, or 4" - ); + (particle_shape >= 1) && (particle_shape <= 4), + "algo.particle_shape can be only 1, 2, 3, or 4"); nox = particle_shape; noy = particle_shape; noz = particle_shape; - } - else{ + } else { WARPX_ABORT_WITH_MESSAGE( "algo.particle_shape must be set in the input file:" " please set algo.particle_shape to 1, 2, 3, or 4"); } - if ((maxLevel() > 0) && (particle_shape > 1) && (do_pml_j_damping == 1)) - { - ablastr::warn_manager::WMRecordWarning("Particles", - "When algo.particle_shape > 1," - "some numerical artifact will be present at the interface between coarse and fine patch." - "We recommend setting algo.particle_shape = 1 in order to avoid this issue"); + if ((maxLevel() > 0) && (particle_shape > 1) && + (do_pml_j_damping == 1)) { + ablastr::warn_manager::WMRecordWarning( + "Particles", "When algo.particle_shape > 1," + "some numerical artifact will be present at " + "the interface between coarse and fine patch." + "We recommend setting algo.particle_shape = 1 " + "in order to avoid this issue"); } - // default sort interval for particles if species or lasers vector is not empty + // default sort interval for particles if species or lasers vector + // is not empty #ifdef AMREX_USE_GPU sort_intervals_string_vec = {"4"}; #else @@ -1365,31 +1497,32 @@ WarpX::ReadParameters () const amrex::ParmParse pp_warpx("warpx"); pp_warpx.queryarr("sort_intervals", sort_intervals_string_vec); - sort_intervals = utils::parser::IntervalsParser(sort_intervals_string_vec); + sort_intervals = + utils::parser::IntervalsParser(sort_intervals_string_vec); - Vector vect_sort_bin_size(AMREX_SPACEDIM,1); + Vector vect_sort_bin_size(AMREX_SPACEDIM, 1); const bool sort_bin_size_is_specified = - utils::parser::queryArrWithParser( - pp_warpx, "sort_bin_size", - vect_sort_bin_size, 0, AMREX_SPACEDIM); - if (sort_bin_size_is_specified){ - for (int i=0; i vect_sort_idx_type(AMREX_SPACEDIM,0); + pp_warpx.query("sort_particles_for_deposition", + sort_particles_for_deposition); + Vector vect_sort_idx_type(AMREX_SPACEDIM, 0); const bool sort_idx_type_is_specified = - utils::parser::queryArrWithParser( - pp_warpx, "sort_idx_type", - vect_sort_idx_type, 0, AMREX_SPACEDIM); - if (sort_idx_type_is_specified){ - for (int i=0; i 0 && WarpX::grid_type != GridType::Collocated) - { + // (note that when warpx.grid_type=collocated, finite-order centering is + // not used anyways) + if (maxLevel() > 0 && WarpX::grid_type != GridType::Collocated) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - field_centering_nox == 2 && field_centering_noy == 2 && field_centering_noz == 2, - "High-order centering of fields (order > 2) is not implemented with mesh refinement"); + field_centering_nox == 2 && field_centering_noy == 2 && + field_centering_noz == 2, + "High-order centering of fields (order > 2) is not implemented " + "with mesh refinement"); } } - if (electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) - { + if (electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { const ParmParse pp_psatd("psatd"); pp_psatd.query("periodic_single_box_fft", fft_periodic_single_box); @@ -1444,68 +1576,78 @@ WarpX::ReadParameters () pp_psatd.query("noy", noy_str); pp_psatd.query("noz", noz_str); - if(nox_str == "inf") { + if (nox_str == "inf") { nox_fft = -1; } else { utils::parser::queryWithParser(pp_psatd, "nox", nox_fft); } - if(noy_str == "inf") { + if (noy_str == "inf") { noy_fft = -1; } else { utils::parser::queryWithParser(pp_psatd, "noy", noy_fft); } - if(noz_str == "inf") { + if (noz_str == "inf") { noz_fft = -1; } else { utils::parser::queryWithParser(pp_psatd, "noz", noz_fft); } if (!fft_periodic_single_box) { - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(nox_fft > 0, "PSATD order must be finite unless psatd.periodic_single_box_fft is used"); - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(noy_fft > 0, "PSATD order must be finite unless psatd.periodic_single_box_fft is used"); - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(noz_fft > 0, "PSATD order must be finite unless psatd.periodic_single_box_fft is used"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + nox_fft > 0, "PSATD order must be finite unless " + "psatd.periodic_single_box_fft is used"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + noy_fft > 0, "PSATD order must be finite unless " + "psatd.periodic_single_box_fft is used"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + noz_fft > 0, "PSATD order must be finite unless " + "psatd.periodic_single_box_fft is used"); } // Integer that corresponds to the order of the PSATD solution // (whether the PSATD equations are derived from first-order or // second-order solution) - pp_psatd.query_enum_sloppy("solution_type", m_psatd_solution_type, "-_"); + pp_psatd.query_enum_sloppy("solution_type", m_psatd_solution_type, + "-_"); - // Integers that correspond to the time dependency of J (constant, linear) - // and rho (linear, quadratic) for the PSATD algorithm + // Integers that correspond to the time dependency of J (constant, + // linear) and rho (linear, quadratic) for the PSATD algorithm pp_psatd.query_enum_sloppy("J_in_time", J_in_time, "-_"); pp_psatd.query_enum_sloppy("rho_in_time", rho_in_time, "-_"); - if (m_psatd_solution_type != PSATDSolutionType::FirstOrder || !do_multi_J) - { + if (m_psatd_solution_type != PSATDSolutionType::FirstOrder || + !do_multi_J) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( rho_in_time == RhoInTime::Linear, "psatd.rho_in_time=constant not yet implemented, " - "except for psatd.solution_type=first-order and warpx.do_multi_J=1"); + "except for psatd.solution_type=first-order and " + "warpx.do_multi_J=1"); } // Current correction activated by default, unless a charge-conserving // current deposition (Esirkepov, Vay) or the div(E) cleaning scheme // are used - if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Esirkepov || - WarpX::current_deposition_algo == CurrentDepositionAlgo::Villasenor || + if (WarpX::current_deposition_algo == + CurrentDepositionAlgo::Esirkepov || + WarpX::current_deposition_algo == + CurrentDepositionAlgo::Villasenor || WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay || - WarpX::do_dive_cleaning) - { + WarpX::do_dive_cleaning) { current_correction = false; } // TODO Remove this default when current correction will // be implemented for the multi-J algorithm as well. - if (do_multi_J) { current_correction = false; } + if (do_multi_J) { + current_correction = false; + } pp_psatd.query("current_correction", current_correction); if (!current_correction && current_deposition_algo != CurrentDepositionAlgo::Esirkepov && current_deposition_algo != CurrentDepositionAlgo::Villasenor && - current_deposition_algo != CurrentDepositionAlgo::Vay) - { + current_deposition_algo != CurrentDepositionAlgo::Vay) { ablastr::warn_manager::WMRecordWarning( "Algorithms", "The chosen current deposition algorithm does not guarantee" @@ -1518,21 +1660,22 @@ WarpX::ReadParameters () pp_psatd.query("do_time_averaging", fft_do_time_averaging); - if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) - { + if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !fft_periodic_single_box, - "Option algo.current_deposition=vay must be used with psatd.periodic_single_box_fft=0."); + "Option algo.current_deposition=vay must be used with " + "psatd.periodic_single_box_fft=0."); } - if (current_deposition_algo == CurrentDepositionAlgo::Vay) - { + if (current_deposition_algo == CurrentDepositionAlgo::Vay) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !current_correction, - "Options algo.current_deposition=vay and psatd.current_correction=1 cannot be combined together."); + "Options algo.current_deposition=vay and " + "psatd.current_correction=1 cannot be combined together."); } - // Auxiliary: boosted_frame = true if WarpX::gamma_boost is set in the inputs + // Auxiliary: boosted_frame = true if WarpX::gamma_boost is set in the + // inputs const amrex::ParmParse pp_warpx("warpx"); const bool boosted_frame = pp_warpx.query("gamma_boost", gamma_boost); @@ -1542,17 +1685,15 @@ WarpX::ReadParameters () WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !use_default_v_galilean || boosted_frame, - "psatd.use_default_v_galilean = 1 can be used only if WarpX::gamma_boost is also set" - ); + "psatd.use_default_v_galilean = 1 can be used only if " + "WarpX::gamma_boost is also set"); - if (use_default_v_galilean && boosted_frame) - { - m_v_galilean[2] = -std::sqrt(1._rt - 1._rt / (gamma_boost * gamma_boost)); - } - else - { - utils::parser::queryArrWithParser( - pp_psatd, "v_galilean", m_v_galilean, 0, 3); + if (use_default_v_galilean && boosted_frame) { + m_v_galilean[2] = + -std::sqrt(1._rt - 1._rt / (gamma_boost * gamma_boost)); + } else { + utils::parser::queryArrWithParser(pp_psatd, "v_galilean", + m_v_galilean, 0, 3); } // Check whether the default comoving velocity should be used @@ -1561,94 +1702,93 @@ WarpX::ReadParameters () WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !use_default_v_comoving || boosted_frame, - "psatd.use_default_v_comoving = 1 can be used only if WarpX::gamma_boost is also set" - ); + "psatd.use_default_v_comoving = 1 can be used only if " + "WarpX::gamma_boost is also set"); - if (use_default_v_comoving && boosted_frame) - { - m_v_comoving[2] = -std::sqrt(1._rt - 1._rt / (gamma_boost * gamma_boost)); - } - else - { - utils::parser::queryArrWithParser( - pp_psatd, "v_comoving", m_v_comoving, 0, 3); + if (use_default_v_comoving && boosted_frame) { + m_v_comoving[2] = + -std::sqrt(1._rt - 1._rt / (gamma_boost * gamma_boost)); + } else { + utils::parser::queryArrWithParser(pp_psatd, "v_comoving", + m_v_comoving, 0, 3); } // Scale the Galilean/comoving velocity by the speed of light - for (auto& vv : m_v_galilean) { vv*= PhysConst::c; } - for (auto& vv : m_v_comoving) { vv*= PhysConst::c; } + for (auto& vv : m_v_galilean) { + vv *= PhysConst::c; + } + for (auto& vv : m_v_comoving) { + vv *= PhysConst::c; + } const auto v_galilean_is_zero = std::all_of(m_v_galilean.begin(), m_v_galilean.end(), - [](const auto& val){return val == 0.;}); + [] (const auto& val) { return val == 0.; }); const auto v_comoving_is_zero = std::all_of(m_v_comoving.begin(), m_v_comoving.end(), - [](const auto& val){return val == 0.;}); - + [] (const auto& val) { return val == 0.; }); // Galilean and comoving algorithms should not be used together WARPX_ALWAYS_ASSERT_WITH_MESSAGE( v_galilean_is_zero || v_comoving_is_zero, - "Galilean and comoving algorithms should not be used together" - ); - + "Galilean and comoving algorithms should not be used together"); if (current_deposition_algo == CurrentDepositionAlgo::Esirkepov || current_deposition_algo == CurrentDepositionAlgo::Villasenor) { - // The comoving PSATD algorithm is not implemented nor tested with Esirkepov current deposition - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(v_comoving_is_zero, - "charge-conserving current depositions (Esirkepov and Villasenor) cannot be used with the comoving PSATD algorithm"); + // The comoving PSATD algorithm is not implemented nor tested with + // Esirkepov current deposition + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + v_comoving_is_zero, + "charge-conserving current depositions (Esirkepov and " + "Villasenor) cannot be used with the comoving PSATD algorithm"); - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(v_galilean_is_zero, - "charge-conserving current depositions (Esirkepov and Villasenor) cannot be used with the Galilean algorithm."); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + v_galilean_is_zero, + "charge-conserving current depositions (Esirkepov and " + "Villasenor) cannot be used with the Galilean algorithm."); } WARPX_ALWAYS_ASSERT_WITH_MESSAGE( (current_deposition_algo != CurrentDepositionAlgo::Vay) || - v_galilean_is_zero, - "Vay current deposition not implemented for Galilean algorithms" - ); + v_galilean_is_zero, + "Vay current deposition not implemented for Galilean algorithms"); -# ifdef WARPX_DIM_RZ +#ifdef WARPX_DIM_RZ update_with_rho = true; -# else - if (m_v_galilean[0] == 0. && m_v_galilean[1] == 0. && m_v_galilean[2] == 0. && - m_v_comoving[0] == 0. && m_v_comoving[1] == 0. && m_v_comoving[2] == 0.) { +#else + if (m_v_galilean[0] == 0. && m_v_galilean[1] == 0. && + m_v_galilean[2] == 0. && m_v_comoving[0] == 0. && + m_v_comoving[1] == 0. && m_v_comoving[2] == 0.) { update_with_rho = do_dive_cleaning; // standard PSATD + } else { + update_with_rho = true; // Galilean PSATD or comoving PSATD } - else { - update_with_rho = true; // Galilean PSATD or comoving PSATD - } -# endif +#endif // Overwrite update_with_rho with value set in input file pp_psatd.query("update_with_rho", update_with_rho); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( (!do_dive_cleaning) || update_with_rho, - "warpx.do_dive_cleaning = 1 not implemented with psatd.update_with_rho = 0" - ); + "warpx.do_dive_cleaning = 1 not implemented with " + "psatd.update_with_rho = 0"); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( v_comoving_is_zero || update_with_rho, - "psatd.update_with_rho must be equal to 1 for comoving PSATD" - ); + "psatd.update_with_rho must be equal to 1 for comoving PSATD"); - if (do_multi_J) - { + if (do_multi_J) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( v_galilean_is_zero, - "Multi-J algorithm not implemented with Galilean PSATD" - ); + "Multi-J algorithm not implemented with Galilean PSATD"); } - if (J_in_time == JInTime::Linear) - { + if (J_in_time == JInTime::Linear) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - update_with_rho, - "psatd.update_with_rho must be set to 1 when psatd.J_in_time=linear"); + update_with_rho, "psatd.update_with_rho must be set to 1 when " + "psatd.J_in_time=linear"); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( v_galilean_is_zero, @@ -1659,54 +1799,51 @@ WarpX::ReadParameters () "psatd.J_in_time=linear not implemented with comoving PSATD"); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - !current_correction, - "psatd.current_correction=1 not implemented with psatd.J_in_time=linear"); + !current_correction, "psatd.current_correction=1 not " + "implemented with psatd.J_in_time=linear"); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( current_deposition_algo != CurrentDepositionAlgo::Vay, - "algo.current_deposition=vay not implemented with psatd.J_in_time=linear"); + "algo.current_deposition=vay not implemented with " + "psatd.J_in_time=linear"); } - for (int dir = 0; dir < AMREX_SPACEDIM; dir++) - { + for (int dir = 0; dir < AMREX_SPACEDIM; dir++) { if (WarpX::field_boundary_lo[dir] == FieldBoundaryType::Damped || - WarpX::field_boundary_hi[dir] == FieldBoundaryType::Damped ) { + WarpX::field_boundary_hi[dir] == FieldBoundaryType::Damped) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - WarpX::field_boundary_lo[dir] == WarpX::field_boundary_hi[dir], - "field boundary in both lo and hi must be set to Damped for PSATD" - ); + WarpX::field_boundary_lo[dir] == + WarpX::field_boundary_hi[dir], + "field boundary in both lo and hi must be set to Damped " + "for PSATD"); } } // Fill guard cells with backward FFTs in directions with field damping - for (int dir = 0; dir < AMREX_SPACEDIM; dir++) - { + for (int dir = 0; dir < AMREX_SPACEDIM; dir++) { if (WarpX::field_boundary_lo[dir] == FieldBoundaryType::Damped || - WarpX::field_boundary_hi[dir] == FieldBoundaryType::Damped) - { + WarpX::field_boundary_hi[dir] == FieldBoundaryType::Damped) { WarpX::m_fill_guards_fields[dir] = 1; } } // Without periodic single box, fill guard cells with backward FFTs, // with current correction or Vay deposition - if (!fft_periodic_single_box) - { + if (!fft_periodic_single_box) { if (current_correction || - current_deposition_algo == CurrentDepositionAlgo::Vay) - { + current_deposition_algo == CurrentDepositionAlgo::Vay) { WarpX::m_fill_guards_current = amrex::IntVect(1); } } } - if (electromagnetic_solver_id != ElectromagneticSolverAlgo::PSATD ) { + if (electromagnetic_solver_id != ElectromagneticSolverAlgo::PSATD) { for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( (WarpX::field_boundary_lo[idim] != FieldBoundaryType::Damped) && - (WarpX::field_boundary_hi[idim] != FieldBoundaryType::Damped), - "FieldBoundaryType::Damped is only supported for PSATD" - ); + (WarpX::field_boundary_hi[idim] != + FieldBoundaryType::Damped), + "FieldBoundaryType::Damped is only supported for PSATD"); } } @@ -1720,24 +1857,21 @@ WarpX::ReadParameters () amrex::Vector slice_hi(AMREX_SPACEDIM); Vector slice_crse_ratio(AMREX_SPACEDIM); // set default slice_crse_ratio // - for (int idim=0; idim < AMREX_SPACEDIM; ++idim ) - { + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { slice_crse_ratio[idim] = 1; } - utils::parser::queryArrWithParser( - pp_slice, "dom_lo", slice_lo, 0, AMREX_SPACEDIM); - utils::parser::queryArrWithParser( - pp_slice, "dom_hi", slice_hi, 0, AMREX_SPACEDIM); - utils::parser::queryArrWithParser( - pp_slice, "coarsening_ratio",slice_crse_ratio,0,AMREX_SPACEDIM); - utils::parser::queryWithParser( - pp_slice, "plot_int",slice_plot_int); + utils::parser::queryArrWithParser(pp_slice, "dom_lo", slice_lo, 0, + AMREX_SPACEDIM); + utils::parser::queryArrWithParser(pp_slice, "dom_hi", slice_hi, 0, + AMREX_SPACEDIM); + utils::parser::queryArrWithParser(pp_slice, "coarsening_ratio", + slice_crse_ratio, 0, AMREX_SPACEDIM); + utils::parser::queryWithParser(pp_slice, "plot_int", slice_plot_int); slice_realbox.setLo(slice_lo); slice_realbox.setHi(slice_hi); - slice_cr_ratio = IntVect(AMREX_D_DECL(1,1,1)); - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) - { - if (slice_crse_ratio[idim] > 1 ) { + slice_cr_ratio = IntVect(AMREX_D_DECL(1, 1, 1)); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if (slice_crse_ratio[idim] > 1) { slice_cr_ratio[idim] = slice_crse_ratio[idim]; } } @@ -1745,8 +1879,7 @@ WarpX::ReadParameters () } void -WarpX::BackwardCompatibility () -{ +WarpX::BackwardCompatibility () { // Auxiliary variables int backward_int; bool backward_bool; @@ -1756,174 +1889,164 @@ WarpX::BackwardCompatibility () const ParmParse pp_amr("amr"); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_amr.query("plot_int", backward_int), - "amr.plot_int is not supported anymore. Please use the new syntax for diagnostics:\n" + "amr.plot_int is not supported anymore. Please use the new syntax for " + "diagnostics:\n" "diagnostics.diags_names = my_diag\n" "my_diag.intervals = 10\n" - "for output every 10 iterations. See documentation for more information" - ); + "for output every 10 iterations. See documentation for more " + "information"); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_amr.query("plot_file", backward_str), "amr.plot_file is not supported anymore. " - "Please use the new syntax for diagnostics, see documentation." - ); + "Please use the new syntax for diagnostics, see documentation."); const ParmParse pp_warpx("warpx"); std::vector backward_strings; WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_warpx.queryarr("fields_to_plot", backward_strings), "warpx.fields_to_plot is not supported anymore. " - "Please use the new syntax for diagnostics, see documentation." - ); + "Please use the new syntax for diagnostics, see documentation."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_warpx.query("plot_finepatch", backward_int), "warpx.plot_finepatch is not supported anymore. " - "Please use the new syntax for diagnostics, see documentation." - ); + "Please use the new syntax for diagnostics, see documentation."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_warpx.query("plot_crsepatch", backward_int), "warpx.plot_crsepatch is not supported anymore. " - "Please use the new syntax for diagnostics, see documentation." - ); + "Please use the new syntax for diagnostics, see documentation."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_warpx.queryarr("load_balance_int", backward_strings), "warpx.load_balance_int is no longer a valid option. " - "Please use the renamed option algo.load_balance_intervals instead." - ); + "Please use the renamed option algo.load_balance_intervals instead."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_warpx.queryarr("load_balance_intervals", backward_strings), "warpx.load_balance_intervals is no longer a valid option. " - "Please use the renamed option algo.load_balance_intervals instead." - ); + "Please use the renamed option algo.load_balance_intervals instead."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - !pp_warpx.query("load_balance_efficiency_ratio_threshold", backward_Real), - "warpx.load_balance_efficiency_ratio_threshold is not supported anymore. " - "Please use the renamed option algo.load_balance_efficiency_ratio_threshold." - ); + !pp_warpx.query("load_balance_efficiency_ratio_threshold", + backward_Real), + "warpx.load_balance_efficiency_ratio_threshold is not supported " + "anymore. " + "Please use the renamed option " + "algo.load_balance_efficiency_ratio_threshold."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_warpx.query("load_balance_with_sfc", backward_int), "warpx.load_balance_with_sfc is not supported anymore. " - "Please use the renamed option algo.load_balance_with_sfc." - ); + "Please use the renamed option algo.load_balance_with_sfc."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_warpx.query("load_balance_knapsack_factor", backward_Real), "warpx.load_balance_knapsack_factor is not supported anymore. " - "Please use the renamed option algo.load_balance_knapsack_factor." - ); + "Please use the renamed option algo.load_balance_knapsack_factor."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_warpx.queryarr("override_sync_int", backward_strings), "warpx.override_sync_int is no longer a valid option. " - "Please use the renamed option warpx.override_sync_intervals instead." - ); + "Please use the renamed option warpx.override_sync_intervals instead."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_warpx.queryarr("sort_int", backward_strings), "warpx.sort_int is no longer a valid option. " - "Please use the renamed option warpx.sort_intervals instead." - ); + "Please use the renamed option warpx.sort_intervals instead."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_warpx.query("do_nodal", backward_int), "warpx.do_nodal is not supported anymore. " - "Please use the flag warpx.grid_type instead." - ); + "Please use the flag warpx.grid_type instead."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_warpx.query("use_kspace_filter", backward_int), "warpx.use_kspace_filter is not supported anymore. " - "Please use the flag use_filter, see documentation." - ); + "Please use the flag use_filter, see documentation."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_warpx.query("do_pml", backward_int), - "do_pml is not supported anymore. Please use boundary.field_lo and boundary.field_hi" - " to set the boundary conditions." - ); + "do_pml is not supported anymore. Please use boundary.field_lo and " + "boundary.field_hi" + " to set the boundary conditions."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_warpx.query("serialize_ics", backward_bool), "warpx.serialize_ics is no longer a valid option. " - "Please use the renamed option warpx.serialize_initial_conditions instead." - ); + "Please use the renamed option warpx.serialize_initial_conditions " + "instead."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_warpx.query("do_back_transformed_diagnostics", backward_int), "Legacy back-transformed diagnostics are not supported anymore. " - "Please use the new syntax for back-transformed diagnostics, see documentation." - ); + "Please use the new syntax for back-transformed diagnostics, see " + "documentation."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_warpx.query("lab_data_directory", backward_str), "Legacy back-transformed diagnostics are not supported anymore. " - "Please use the new syntax for back-transformed diagnostics, see documentation." - ); + "Please use the new syntax for back-transformed diagnostics, see " + "documentation."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_warpx.query("num_snapshots_lab", backward_int), "Legacy back-transformed diagnostics are not supported anymore. " - "Please use the new syntax for back-transformed diagnostics, see documentation." - ); + "Please use the new syntax for back-transformed diagnostics, see " + "documentation."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_warpx.query("dt_snapshots_lab", backward_Real), "Legacy back-transformed diagnostics are not supported anymore. " - "Please use the new syntax for back-transformed diagnostics, see documentation." - ); + "Please use the new syntax for back-transformed diagnostics, see " + "documentation."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_warpx.query("dz_snapshots_lab", backward_Real), "Legacy back-transformed diagnostics are not supported anymore. " - "Please use the new syntax for back-transformed diagnostics, see documentation." - ); + "Please use the new syntax for back-transformed diagnostics, see " + "documentation."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_warpx.query("do_back_transformed_fields", backward_int), "Legacy back-transformed diagnostics are not supported anymore. " - "Please use the new syntax for back-transformed diagnostics, see documentation." - ); + "Please use the new syntax for back-transformed diagnostics, see " + "documentation."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_warpx.query("buffer_size", backward_int), "Legacy back-transformed diagnostics are not supported anymore. " - "Please use the new syntax for back-transformed diagnostics, see documentation." - ); + "Please use the new syntax for back-transformed diagnostics, see " + "documentation."); const ParmParse pp_slice("slice"); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_slice.query("num_slice_snapshots_lab", backward_int), "Legacy back-transformed diagnostics are not supported anymore. " - "Please use the new syntax for back-transformed diagnostics, see documentation." - ); + "Please use the new syntax for back-transformed diagnostics, see " + "documentation."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_slice.query("dt_slice_snapshots_lab", backward_Real), "Legacy back-transformed diagnostics are not supported anymore. " - "Please use the new syntax for back-transformed diagnostics, see documentation." - ); + "Please use the new syntax for back-transformed diagnostics, see " + "documentation."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_slice.query("particle_slice_width_lab", backward_Real), "Legacy back-transformed diagnostics are not supported anymore. " - "Please use the new syntax for back-transformed diagnostics, see documentation." - ); + "Please use the new syntax for back-transformed diagnostics, see " + "documentation."); const ParmParse pp_interpolation("interpolation"); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_interpolation.query("nox", backward_int) && - !pp_interpolation.query("noy", backward_int) && - !pp_interpolation.query("noz", backward_int), + !pp_interpolation.query("noy", backward_int) && + !pp_interpolation.query("noz", backward_int), "interpolation.nox (as well as .noy, .noz) are not supported anymore." - " Please use the new syntax algo.particle_shape instead" - ); + " Please use the new syntax algo.particle_shape instead"); const ParmParse pp_algo("algo"); int backward_mw_solver; @@ -1934,15 +2057,17 @@ WarpX::BackwardCompatibility () const ParmParse pp_particles("particles"); int nspecies; - if (pp_particles.query("nspecies", nspecies)){ - ablastr::warn_manager::WMRecordWarning("Species", - "particles.nspecies is ignored. Just use particles.species_names please.", + if (pp_particles.query("nspecies", nspecies)) { + ablastr::warn_manager::WMRecordWarning( + "Species", + "particles.nspecies is ignored. Just use particles.species_names " + "please.", ablastr::warn_manager::WarnPriority::low); } std::vector backward_sp_names; pp_particles.queryarr("species_names", backward_sp_names); - for(const std::string& speciesiter : backward_sp_names){ + for (const std::string& speciesiter : backward_sp_names) { const ParmParse pp_species(speciesiter); std::vector backward_vel; std::stringstream ssspecies; @@ -1953,8 +2078,9 @@ WarpX::BackwardCompatibility () ssspecies << "'" << speciesiter << ".multiple_particles_u' ."; WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_species.queryarr("multiple_particles_vel_x", backward_vel) && - !pp_species.queryarr("multiple_particles_vel_y", backward_vel) && - !pp_species.queryarr("multiple_particles_vel_z", backward_vel), + !pp_species.queryarr("multiple_particles_vel_y", + backward_vel) && + !pp_species.queryarr("multiple_particles_vel_z", backward_vel), ssspecies.str()); ssspecies.str(""); @@ -1970,17 +2096,19 @@ WarpX::BackwardCompatibility () const ParmParse pp_collisions("collisions"); int ncollisions; - if (pp_collisions.query("ncollisions", ncollisions)){ - ablastr::warn_manager::WMRecordWarning("Collisions", - "collisions.ncollisions is ignored. Just use particles.collision_names please.", + if (pp_collisions.query("ncollisions", ncollisions)) { + ablastr::warn_manager::WMRecordWarning( + "Collisions", + "collisions.ncollisions is ignored. Just use " + "particles.collision_names please.", ablastr::warn_manager::WarnPriority::low); } const ParmParse pp_lasers("lasers"); int nlasers; - if (pp_lasers.query("nlasers", nlasers)){ - ablastr::warn_manager::WMRecordWarning("Laser", - "lasers.nlasers is ignored. Just use lasers.names please.", + if (pp_lasers.query("nlasers", nlasers)) { + ablastr::warn_manager::WMRecordWarning( + "Laser", "lasers.nlasers is ignored. Just use lasers.names please.", ablastr::warn_manager::WarnPriority::low); } } @@ -1988,29 +2116,27 @@ WarpX::BackwardCompatibility () // This is a virtual function. void WarpX::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& new_grids, - const DistributionMapping& new_dmap) -{ + const DistributionMapping& new_dmap) { AllocLevelData(lev, new_grids, new_dmap); InitLevelData(lev, time); } // This is a virtual function. void -WarpX::MakeNewLevelFromCoarse (int /*lev*/, amrex::Real /*time*/, const amrex::BoxArray& /*ba*/, - const amrex::DistributionMapping& /*dm*/) -{ +WarpX::MakeNewLevelFromCoarse (int /*lev*/, amrex::Real /*time*/, + const amrex::BoxArray& /*ba*/, + const amrex::DistributionMapping& /*dm*/) { WARPX_ABORT_WITH_MESSAGE("MakeNewLevelFromCoarse: To be implemented"); } void -WarpX::ClearLevel (int lev) -{ +WarpX::ClearLevel (int lev) { m_fields.clear_level(lev); for (int i = 0; i < 3; ++i) { - Efield_dotMask [lev][i].reset(); - Bfield_dotMask [lev][i].reset(); - Afield_dotMask [lev][i].reset(); + Efield_dotMask[lev][i].reset(); + Bfield_dotMask[lev][i].reset(); + Afield_dotMask[lev][i].reset(); } phi_dotMask[lev].reset(); @@ -2027,58 +2153,42 @@ WarpX::ClearLevel (int lev) } void -WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& dm) -{ - const bool aux_is_nodal = (field_gathering_algo == GatheringAlgo::MomentumConserving); +WarpX::AllocLevelData (int lev, const BoxArray& ba, + const DistributionMapping& dm) { + const bool aux_is_nodal = + (field_gathering_algo == GatheringAlgo::MomentumConserving); const Real* dx = Geom(lev).CellSize(); // Initialize filter before guard cells manager // (needs info on length of filter's stencil) - if (use_filter) - { + if (use_filter) { InitFilter(); } guard_cells.Init( - dt[lev], - dx, - m_do_subcycling, - WarpX::use_fdtd_nci_corr, - grid_type, - do_moving_window, - moving_window_dir, - WarpX::nox, - nox_fft, noy_fft, noz_fft, - NCIGodfreyFilter::m_stencil_width, - electromagnetic_solver_id, - maxLevel(), - WarpX::m_v_galilean, - WarpX::m_v_comoving, - m_safe_guard_cells, - WarpX::do_multi_J, - WarpX::fft_do_time_averaging, + dt[lev], dx, m_do_subcycling, WarpX::use_fdtd_nci_corr, grid_type, + do_moving_window, moving_window_dir, WarpX::nox, nox_fft, noy_fft, + noz_fft, NCIGodfreyFilter::m_stencil_width, electromagnetic_solver_id, + maxLevel(), WarpX::m_v_galilean, WarpX::m_v_comoving, + m_safe_guard_cells, WarpX::do_multi_J, WarpX::fft_do_time_averaging, ::isAnyBoundaryPML(field_boundary_lo, field_boundary_hi), - WarpX::do_pml_in_domain, - WarpX::pml_ncell, - this->refRatio(), - use_filter, + WarpX::do_pml_in_domain, WarpX::pml_ncell, this->refRatio(), use_filter, bilinear_filter.stencil_length_each_dir); #ifdef AMREX_USE_EB bool const eb_enabled = EB::enabled(); if (eb_enabled) { int const max_guard = guard_cells.ng_FieldSolver.max(); - m_field_factory[lev] = amrex::makeEBFabFactory(Geom(lev), ba, dm, - {max_guard, max_guard, max_guard}, - amrex::EBSupport::full); + m_field_factory[lev] = amrex::makeEBFabFactory( + Geom(lev), ba, dm, {max_guard, max_guard, max_guard}, + amrex::EBSupport::full); } else #endif { m_field_factory[lev] = std::make_unique(); } - if (mypc->nSpeciesDepositOnMainGrid() && n_current_deposition_buffer == 0) { n_current_deposition_buffer = 1; // This forces the allocation of buffers and allows the code associated @@ -2096,18 +2206,18 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d } 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_alloc_G, aux_is_nodal); + guard_cells.ng_alloc_Rho, guard_cells.ng_alloc_F, + guard_cells.ng_alloc_G, aux_is_nodal); m_accelerator_lattice[lev] = std::make_unique(); m_accelerator_lattice[lev]->InitElementFinder(lev, gamma_boost, ba, dm); - } void -WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm, - const IntVect& ngEB, IntVect& ngJ, const IntVect& ngRho, - const IntVect& ngF, const IntVect& ngG, const bool aux_is_nodal) -{ +WarpX::AllocLevelMFs (int lev, const BoxArray& ba, + const DistributionMapping& dm, const IntVect& ngEB, + IntVect& ngJ, const IntVect& ngRho, const IntVect& ngF, + const IntVect& ngG, const bool aux_is_nodal) { using ablastr::fields::Direction; // Declare nodal flags @@ -2119,8 +2229,9 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm amrex::IntVect F_nodal_flag, G_nodal_flag; // Set nodal flags -#if defined(WARPX_DIM_1D_Z) - // AMReX convention: x = missing dimension, y = missing dimension, z = only dimension +#if defined(WARPX_DIM_1D_Z) + // AMReX convention: x = missing dimension, y = missing dimension, z = only + // dimension Ex_nodal_flag = IntVect(1); Ey_nodal_flag = IntVect(1); Ez_nodal_flag = IntVect(0); @@ -2130,66 +2241,67 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm jx_nodal_flag = IntVect(1); jy_nodal_flag = IntVect(1); jz_nodal_flag = IntVect(0); -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - // AMReX convention: x = first dimension, y = missing dimension, z = second dimension - Ex_nodal_flag = IntVect(0,1); - Ey_nodal_flag = IntVect(1,1); - Ez_nodal_flag = IntVect(1,0); - Bx_nodal_flag = IntVect(1,0); - By_nodal_flag = IntVect(0,0); - Bz_nodal_flag = IntVect(0,1); - jx_nodal_flag = IntVect(0,1); - jy_nodal_flag = IntVect(1,1); - jz_nodal_flag = IntVect(1,0); +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + // AMReX convention: x = first dimension, y = missing dimension, z = second + // dimension + Ex_nodal_flag = IntVect(0, 1); + Ey_nodal_flag = IntVect(1, 1); + Ez_nodal_flag = IntVect(1, 0); + Bx_nodal_flag = IntVect(1, 0); + By_nodal_flag = IntVect(0, 0); + Bz_nodal_flag = IntVect(0, 1); + jx_nodal_flag = IntVect(0, 1); + jy_nodal_flag = IntVect(1, 1); + jz_nodal_flag = IntVect(1, 0); #elif defined(WARPX_DIM_3D) - Ex_nodal_flag = IntVect(0,1,1); - Ey_nodal_flag = IntVect(1,0,1); - Ez_nodal_flag = IntVect(1,1,0); - Bx_nodal_flag = IntVect(1,0,0); - By_nodal_flag = IntVect(0,1,0); - Bz_nodal_flag = IntVect(0,0,1); - jx_nodal_flag = IntVect(0,1,1); - jy_nodal_flag = IntVect(1,0,1); - jz_nodal_flag = IntVect(1,1,0); + Ex_nodal_flag = IntVect(0, 1, 1); + Ey_nodal_flag = IntVect(1, 0, 1); + Ez_nodal_flag = IntVect(1, 1, 0); + Bx_nodal_flag = IntVect(1, 0, 0); + By_nodal_flag = IntVect(0, 1, 0); + Bz_nodal_flag = IntVect(0, 0, 1); + jx_nodal_flag = IntVect(0, 1, 1); + jy_nodal_flag = IntVect(1, 0, 1); + jz_nodal_flag = IntVect(1, 1, 0); #endif - if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) - { - jx_nodal_flag = IntVect::TheNodeVector(); - jy_nodal_flag = IntVect::TheNodeVector(); - jz_nodal_flag = IntVect::TheNodeVector(); + if (electrostatic_solver_id == + ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) { + jx_nodal_flag = IntVect::TheNodeVector(); + jy_nodal_flag = IntVect::TheNodeVector(); + jz_nodal_flag = IntVect::TheNodeVector(); ngJ = ngRho; } - rho_nodal_flag = IntVect( AMREX_D_DECL(1,1,1) ); + rho_nodal_flag = IntVect(AMREX_D_DECL(1, 1, 1)); phi_nodal_flag = IntVect::TheNodeVector(); F_nodal_flag = amrex::IntVect::TheNodeVector(); G_nodal_flag = amrex::IntVect::TheCellVector(); // Overwrite nodal flags if necessary if (grid_type == GridType::Collocated) { - Ex_nodal_flag = IntVect::TheNodeVector(); - Ey_nodal_flag = IntVect::TheNodeVector(); - Ez_nodal_flag = IntVect::TheNodeVector(); - Bx_nodal_flag = IntVect::TheNodeVector(); - By_nodal_flag = IntVect::TheNodeVector(); - Bz_nodal_flag = IntVect::TheNodeVector(); - jx_nodal_flag = IntVect::TheNodeVector(); - jy_nodal_flag = IntVect::TheNodeVector(); - jz_nodal_flag = IntVect::TheNodeVector(); + Ex_nodal_flag = IntVect::TheNodeVector(); + Ey_nodal_flag = IntVect::TheNodeVector(); + Ez_nodal_flag = IntVect::TheNodeVector(); + Bx_nodal_flag = IntVect::TheNodeVector(); + By_nodal_flag = IntVect::TheNodeVector(); + Bz_nodal_flag = IntVect::TheNodeVector(); + jx_nodal_flag = IntVect::TheNodeVector(); + jy_nodal_flag = IntVect::TheNodeVector(); + jz_nodal_flag = IntVect::TheNodeVector(); rho_nodal_flag = IntVect::TheNodeVector(); G_nodal_flag = amrex::IntVect::TheNodeVector(); } #ifdef WARPX_DIM_RZ if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { // Force cell-centered IndexType in r and z - Ex_nodal_flag = IntVect::TheCellVector(); - Ey_nodal_flag = IntVect::TheCellVector(); - Ez_nodal_flag = IntVect::TheCellVector(); - Bx_nodal_flag = IntVect::TheCellVector(); - By_nodal_flag = IntVect::TheCellVector(); - Bz_nodal_flag = IntVect::TheCellVector(); - jx_nodal_flag = IntVect::TheCellVector(); - jy_nodal_flag = IntVect::TheCellVector(); - jz_nodal_flag = IntVect::TheCellVector(); + Ex_nodal_flag = IntVect::TheCellVector(); + Ey_nodal_flag = IntVect::TheCellVector(); + Ez_nodal_flag = IntVect::TheCellVector(); + Bx_nodal_flag = IntVect::TheCellVector(); + By_nodal_flag = IntVect::TheCellVector(); + Bz_nodal_flag = IntVect::TheCellVector(); + jx_nodal_flag = IntVect::TheCellVector(); + jy_nodal_flag = IntVect::TheCellVector(); + jz_nodal_flag = IntVect::TheCellVector(); rho_nodal_flag = IntVect::TheCellVector(); F_nodal_flag = IntVect::TheCellVector(); G_nodal_flag = IntVect::TheCellVector(); @@ -2200,626 +2312,874 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm // Component 0 is mode 0. // Odd components are the real parts. // Even components are the imaginary parts. - ncomps = n_rz_azimuthal_modes*2 - 1; + ncomps = n_rz_azimuthal_modes * 2 - 1; #endif - // Set global rho nodal flag to know about rho index type when rho MultiFab is not allocated + // Set global rho nodal flag to know about rho index type when rho MultiFab + // is not allocated m_rho_nodal_flag = rho_nodal_flag; // // The fine patch // - const std::array dx = CellSize(lev); - - m_fields.alloc_init(FieldType::Bfield_fp, Direction{0}, lev, amrex::convert(ba, Bx_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Bfield_fp, Direction{1}, lev, amrex::convert(ba, By_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Bfield_fp, Direction{2}, lev, amrex::convert(ba, Bz_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - - m_fields.alloc_init(FieldType::Efield_fp, Direction{0}, lev, amrex::convert(ba, Ex_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Efield_fp, Direction{1}, lev, amrex::convert(ba, Ey_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Efield_fp, Direction{2}, lev, amrex::convert(ba, Ez_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - - m_fields.alloc_init(FieldType::current_fp, Direction{0}, lev, amrex::convert(ba, jx_nodal_flag), dm, ncomps, ngJ, 0.0_rt); - m_fields.alloc_init(FieldType::current_fp, Direction{1}, lev, amrex::convert(ba, jy_nodal_flag), dm, ncomps, ngJ, 0.0_rt); - m_fields.alloc_init(FieldType::current_fp, Direction{2}, lev, amrex::convert(ba, jz_nodal_flag), dm, ncomps, ngJ, 0.0_rt); - - if (do_current_centering) - { - amrex::BoxArray const& nodal_ba = amrex::convert(ba, amrex::IntVect::TheNodeVector()); - m_fields.alloc_init(FieldType::current_fp_nodal, Direction{0}, lev, nodal_ba, dm, ncomps, ngJ, 0.0_rt); - m_fields.alloc_init(FieldType::current_fp_nodal, Direction{1}, lev, nodal_ba, dm, ncomps, ngJ, 0.0_rt); - m_fields.alloc_init(FieldType::current_fp_nodal, Direction{2}, lev, nodal_ba, dm, ncomps, ngJ, 0.0_rt); - } - - if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) - { - m_fields.alloc_init(FieldType::current_fp_vay, Direction{0}, lev, amrex::convert(ba, rho_nodal_flag), dm, ncomps, ngJ, 0.0_rt); - m_fields.alloc_init(FieldType::current_fp_vay, Direction{1}, lev, amrex::convert(ba, rho_nodal_flag), dm, ncomps, ngJ, 0.0_rt); - m_fields.alloc_init(FieldType::current_fp_vay, Direction{2}, lev, amrex::convert(ba, rho_nodal_flag), dm, ncomps, ngJ, 0.0_rt); - } - - if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) - { - m_fields.alloc_init(FieldType::vector_potential_fp_nodal, Direction{0}, lev, amrex::convert(ba, rho_nodal_flag), dm, ncomps, ngRho, 0.0_rt); - m_fields.alloc_init(FieldType::vector_potential_fp_nodal, Direction{1}, lev, amrex::convert(ba, rho_nodal_flag), dm, ncomps, ngRho, 0.0_rt); - m_fields.alloc_init(FieldType::vector_potential_fp_nodal, Direction{2}, lev, amrex::convert(ba, rho_nodal_flag), dm, ncomps, ngRho, 0.0_rt); + const std::array dx = CellSize(lev); + + m_fields.alloc_init(FieldType::Bfield_fp, Direction{0}, lev, + amrex::convert(ba, Bx_nodal_flag), dm, ncomps, ngEB, + 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_fp, Direction{1}, lev, + amrex::convert(ba, By_nodal_flag), dm, ncomps, ngEB, + 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_fp, Direction{2}, lev, + amrex::convert(ba, Bz_nodal_flag), dm, ncomps, ngEB, + 0.0_rt); + + m_fields.alloc_init(FieldType::Efield_fp, Direction{0}, lev, + amrex::convert(ba, Ex_nodal_flag), dm, ncomps, ngEB, + 0.0_rt); + m_fields.alloc_init(FieldType::Efield_fp, Direction{1}, lev, + amrex::convert(ba, Ey_nodal_flag), dm, ncomps, ngEB, + 0.0_rt); + m_fields.alloc_init(FieldType::Efield_fp, Direction{2}, lev, + amrex::convert(ba, Ez_nodal_flag), dm, ncomps, ngEB, + 0.0_rt); + + m_fields.alloc_init(FieldType::current_fp, Direction{0}, lev, + amrex::convert(ba, jx_nodal_flag), dm, ncomps, ngJ, + 0.0_rt); + m_fields.alloc_init(FieldType::current_fp, Direction{1}, lev, + amrex::convert(ba, jy_nodal_flag), dm, ncomps, ngJ, + 0.0_rt); + m_fields.alloc_init(FieldType::current_fp, Direction{2}, lev, + amrex::convert(ba, jz_nodal_flag), dm, ncomps, ngJ, + 0.0_rt); + + if (do_current_centering) { + amrex::BoxArray const& nodal_ba = + amrex::convert(ba, amrex::IntVect::TheNodeVector()); + m_fields.alloc_init(FieldType::current_fp_nodal, Direction{0}, lev, + nodal_ba, dm, ncomps, ngJ, 0.0_rt); + m_fields.alloc_init(FieldType::current_fp_nodal, Direction{1}, lev, + nodal_ba, dm, ncomps, ngJ, 0.0_rt); + m_fields.alloc_init(FieldType::current_fp_nodal, Direction{2}, lev, + nodal_ba, dm, ncomps, ngJ, 0.0_rt); + } + + if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) { + m_fields.alloc_init(FieldType::current_fp_vay, Direction{0}, lev, + amrex::convert(ba, rho_nodal_flag), dm, ncomps, ngJ, + 0.0_rt); + m_fields.alloc_init(FieldType::current_fp_vay, Direction{1}, lev, + amrex::convert(ba, rho_nodal_flag), dm, ncomps, ngJ, + 0.0_rt); + m_fields.alloc_init(FieldType::current_fp_vay, Direction{2}, lev, + amrex::convert(ba, rho_nodal_flag), dm, ncomps, ngJ, + 0.0_rt); + } + + if (electrostatic_solver_id == + ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) { + m_fields.alloc_init(FieldType::vector_potential_fp_nodal, Direction{0}, + lev, amrex::convert(ba, rho_nodal_flag), dm, ncomps, + ngRho, 0.0_rt); + m_fields.alloc_init(FieldType::vector_potential_fp_nodal, Direction{1}, + lev, amrex::convert(ba, rho_nodal_flag), dm, ncomps, + ngRho, 0.0_rt); + m_fields.alloc_init(FieldType::vector_potential_fp_nodal, Direction{2}, + lev, amrex::convert(ba, rho_nodal_flag), dm, ncomps, + ngRho, 0.0_rt); // Memory buffers for computing magnetostatic fields - // Vector Potential A and previous step. Time buffer needed for computing dA/dt to first order - m_fields.alloc_init(FieldType::vector_potential_grad_buf_e_stag, Direction{0}, lev, amrex::convert(ba, Ex_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::vector_potential_grad_buf_e_stag, Direction{1}, lev, amrex::convert(ba, Ey_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::vector_potential_grad_buf_e_stag, Direction{2}, lev, amrex::convert(ba, Ez_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - - m_fields.alloc_init(FieldType::vector_potential_grad_buf_b_stag, Direction{0}, lev, amrex::convert(ba, Bx_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::vector_potential_grad_buf_b_stag, Direction{1}, lev, amrex::convert(ba, By_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::vector_potential_grad_buf_b_stag, Direction{2}, lev, amrex::convert(ba, Bz_nodal_flag), dm, ncomps, ngEB, 0.0_rt); + // Vector Potential A and previous step. Time buffer needed for + // computing dA/dt to first order + m_fields.alloc_init( + FieldType::vector_potential_grad_buf_e_stag, Direction{0}, lev, + amrex::convert(ba, Ex_nodal_flag), dm, ncomps, ngEB, 0.0_rt); + m_fields.alloc_init( + FieldType::vector_potential_grad_buf_e_stag, Direction{1}, lev, + amrex::convert(ba, Ey_nodal_flag), dm, ncomps, ngEB, 0.0_rt); + m_fields.alloc_init( + FieldType::vector_potential_grad_buf_e_stag, Direction{2}, lev, + amrex::convert(ba, Ez_nodal_flag), dm, ncomps, ngEB, 0.0_rt); + + m_fields.alloc_init( + FieldType::vector_potential_grad_buf_b_stag, Direction{0}, lev, + amrex::convert(ba, Bx_nodal_flag), dm, ncomps, ngEB, 0.0_rt); + m_fields.alloc_init( + FieldType::vector_potential_grad_buf_b_stag, Direction{1}, lev, + amrex::convert(ba, By_nodal_flag), dm, ncomps, ngEB, 0.0_rt); + m_fields.alloc_init( + FieldType::vector_potential_grad_buf_b_stag, Direction{2}, lev, + amrex::convert(ba, Bz_nodal_flag), dm, ncomps, ngEB, 0.0_rt); } // Allocate extra multifabs needed by the kinetic-fluid hybrid algorithm. - if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC) - { + if (WarpX::electromagnetic_solver_id == + ElectromagneticSolverAlgo::HybridPIC) { m_hybrid_pic_model->AllocateLevelMFs( - m_fields, - lev, ba, dm, ncomps, ngJ, ngRho, ngEB, jx_nodal_flag, jy_nodal_flag, - jz_nodal_flag, rho_nodal_flag, Ex_nodal_flag, Ey_nodal_flag, Ez_nodal_flag, - Bx_nodal_flag, By_nodal_flag, Bz_nodal_flag - ); + m_fields, lev, ba, dm, ncomps, ngJ, ngRho, ngEB, jx_nodal_flag, + jy_nodal_flag, jz_nodal_flag, rho_nodal_flag, Ex_nodal_flag, + Ey_nodal_flag, Ez_nodal_flag, Bx_nodal_flag, By_nodal_flag, + Bz_nodal_flag); } // Allocate extra multifabs needed for fluids if (do_fluid_species) { myfl->AllocateLevelMFs(m_fields, ba, dm, lev); - auto & warpx = GetInstance(); + auto& warpx = GetInstance(); const amrex::Real cur_time = warpx.gett_new(lev); myfl->InitData(m_fields, geom[lev].Domain(), cur_time, lev); } // Allocate extra multifabs for macroscopic properties of the medium if (m_em_solver_medium == MediumForEM::Macroscopic) { - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( lev==0, + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + lev == 0, "Macroscopic properties are not supported with mesh refinement."); m_macroscopic_properties->AllocateLevelMFs(ba, dm, ngEB); } - if (fft_do_time_averaging) - { - m_fields.alloc_init(FieldType::Bfield_avg_fp, Direction{0}, lev, amrex::convert(ba, Bx_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Bfield_avg_fp, Direction{1}, lev, amrex::convert(ba, By_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Bfield_avg_fp, Direction{2}, lev, amrex::convert(ba, Bz_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - - m_fields.alloc_init(FieldType::Efield_avg_fp, Direction{0}, lev, amrex::convert(ba, Ex_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Efield_avg_fp, Direction{1}, lev, amrex::convert(ba, Ey_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Efield_avg_fp, Direction{2}, lev, amrex::convert(ba, Ez_nodal_flag), dm, ncomps, ngEB, 0.0_rt); + if (fft_do_time_averaging) { + m_fields.alloc_init(FieldType::Bfield_avg_fp, Direction{0}, lev, + amrex::convert(ba, Bx_nodal_flag), dm, ncomps, ngEB, + 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_avg_fp, Direction{1}, lev, + amrex::convert(ba, By_nodal_flag), dm, ncomps, ngEB, + 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_avg_fp, Direction{2}, lev, + amrex::convert(ba, Bz_nodal_flag), dm, ncomps, ngEB, + 0.0_rt); + + m_fields.alloc_init(FieldType::Efield_avg_fp, Direction{0}, lev, + amrex::convert(ba, Ex_nodal_flag), dm, ncomps, ngEB, + 0.0_rt); + m_fields.alloc_init(FieldType::Efield_avg_fp, Direction{1}, lev, + amrex::convert(ba, Ey_nodal_flag), dm, ncomps, ngEB, + 0.0_rt); + m_fields.alloc_init(FieldType::Efield_avg_fp, Direction{2}, lev, + amrex::convert(ba, Ez_nodal_flag), dm, ncomps, ngEB, + 0.0_rt); } if (EB::enabled()) { constexpr int nc_ls = 1; amrex::IntVect const ng_ls(2); - //EB level set - m_fields.alloc_init(FieldType::distance_to_eb, lev, amrex::convert(ba, IntVect::TheNodeVector()), dm, nc_ls, ng_ls, 0.0_rt); + // EB level set + m_fields.alloc_init(FieldType::distance_to_eb, lev, + amrex::convert(ba, IntVect::TheNodeVector()), dm, + nc_ls, ng_ls, 0.0_rt); // Whether to reduce the particle shape to order 1 when close to the EB - AllocInitMultiFab(m_eb_reduce_particle_shape[lev], amrex::convert(ba, IntVect::TheCellVector()), dm, ncomps, - ngRho, lev, "m_eb_reduce_particle_shape"); + AllocInitMultiFab(m_eb_reduce_particle_shape[lev], + amrex::convert(ba, IntVect::TheCellVector()), dm, + ncomps, ngRho, lev, "m_eb_reduce_particle_shape"); // EB info are needed only at the finest level if (lev == maxLevel()) { - if (WarpX::electromagnetic_solver_id != ElectromagneticSolverAlgo::PSATD) { - - AllocInitMultiFab(m_eb_update_E[lev][0], amrex::convert(ba, Ex_nodal_flag), dm, ncomps, - guard_cells.ng_FieldSolver, lev, "m_eb_update_E[x]"); - AllocInitMultiFab(m_eb_update_E[lev][1], amrex::convert(ba, Ey_nodal_flag), dm, ncomps, - guard_cells.ng_FieldSolver, lev, "m_eb_update_E[y]"); - AllocInitMultiFab(m_eb_update_E[lev][2], amrex::convert(ba, Ez_nodal_flag), dm, ncomps, - guard_cells.ng_FieldSolver, lev, "m_eb_update_E[z]"); - - AllocInitMultiFab(m_eb_update_B[lev][0], amrex::convert(ba, Bx_nodal_flag), dm, ncomps, - guard_cells.ng_FieldSolver, lev, "m_eb_update_B[x]"); - AllocInitMultiFab(m_eb_update_B[lev][1], amrex::convert(ba, By_nodal_flag), dm, ncomps, - guard_cells.ng_FieldSolver, lev, "m_eb_update_B[y]"); - AllocInitMultiFab(m_eb_update_B[lev][2], amrex::convert(ba, Bz_nodal_flag), dm, ncomps, - guard_cells.ng_FieldSolver, lev, "m_eb_update_B[z]"); + if (WarpX::electromagnetic_solver_id != + ElectromagneticSolverAlgo::PSATD) { + + AllocInitMultiFab(m_eb_update_E[lev][0], + amrex::convert(ba, Ex_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, + "m_eb_update_E[x]"); + AllocInitMultiFab(m_eb_update_E[lev][1], + amrex::convert(ba, Ey_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, + "m_eb_update_E[y]"); + AllocInitMultiFab(m_eb_update_E[lev][2], + amrex::convert(ba, Ez_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, + "m_eb_update_E[z]"); + + AllocInitMultiFab(m_eb_update_B[lev][0], + amrex::convert(ba, Bx_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, + "m_eb_update_B[x]"); + AllocInitMultiFab(m_eb_update_B[lev][1], + amrex::convert(ba, By_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, + "m_eb_update_B[y]"); + AllocInitMultiFab(m_eb_update_B[lev][2], + amrex::convert(ba, Bz_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, + "m_eb_update_B[z]"); } - if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT) { + if (WarpX::electromagnetic_solver_id == + ElectromagneticSolverAlgo::ECT) { //! EB: Lengths of the mesh edges - m_fields.alloc_init(FieldType::edge_lengths, Direction{0}, lev, amrex::convert(ba, Ex_nodal_flag), - dm, ncomps, guard_cells.ng_FieldSolver, 0.0_rt); - m_fields.alloc_init(FieldType::edge_lengths, Direction{1}, lev, amrex::convert(ba, Ey_nodal_flag), - dm, ncomps, guard_cells.ng_FieldSolver, 0.0_rt); - m_fields.alloc_init(FieldType::edge_lengths, Direction{2}, lev, amrex::convert(ba, Ez_nodal_flag), - dm, ncomps, guard_cells.ng_FieldSolver, 0.0_rt); + m_fields.alloc_init(FieldType::edge_lengths, Direction{0}, lev, + amrex::convert(ba, Ex_nodal_flag), dm, + ncomps, guard_cells.ng_FieldSolver, 0.0_rt); + m_fields.alloc_init(FieldType::edge_lengths, Direction{1}, lev, + amrex::convert(ba, Ey_nodal_flag), dm, + ncomps, guard_cells.ng_FieldSolver, 0.0_rt); + m_fields.alloc_init(FieldType::edge_lengths, Direction{2}, lev, + amrex::convert(ba, Ez_nodal_flag), dm, + ncomps, guard_cells.ng_FieldSolver, 0.0_rt); //! EB: Areas of the mesh faces - m_fields.alloc_init(FieldType::face_areas, Direction{0}, lev, amrex::convert(ba, Bx_nodal_flag), - dm, ncomps, guard_cells.ng_FieldSolver, 0.0_rt); - m_fields.alloc_init(FieldType::face_areas, Direction{1}, lev, amrex::convert(ba, By_nodal_flag), - dm, ncomps, guard_cells.ng_FieldSolver, 0.0_rt); - m_fields.alloc_init(FieldType::face_areas, Direction{2}, lev, amrex::convert(ba, Bz_nodal_flag), - dm, ncomps, guard_cells.ng_FieldSolver, 0.0_rt); - - AllocInitMultiFab(m_flag_info_face[lev][0], amrex::convert(ba, Bx_nodal_flag), dm, ncomps, - guard_cells.ng_FieldSolver, lev, "m_flag_info_face[x]"); - AllocInitMultiFab(m_flag_info_face[lev][1], amrex::convert(ba, By_nodal_flag), dm, ncomps, - guard_cells.ng_FieldSolver, lev, "m_flag_info_face[y]"); - AllocInitMultiFab(m_flag_info_face[lev][2], amrex::convert(ba, Bz_nodal_flag), dm, ncomps, - guard_cells.ng_FieldSolver, lev, "m_flag_info_face[z]"); - AllocInitMultiFab(m_flag_ext_face[lev][0], amrex::convert(ba, Bx_nodal_flag), dm, ncomps, - guard_cells.ng_FieldSolver, lev, "m_flag_ext_face[x]"); - AllocInitMultiFab(m_flag_ext_face[lev][1], amrex::convert(ba, By_nodal_flag), dm, ncomps, - guard_cells.ng_FieldSolver, lev, "m_flag_ext_face[y]"); - AllocInitMultiFab(m_flag_ext_face[lev][2], amrex::convert(ba, Bz_nodal_flag), dm, ncomps, - guard_cells.ng_FieldSolver, lev, "m_flag_ext_face[z]"); - - /** EB: area_mod contains the modified areas of the mesh faces, i.e. if a face is enlarged it - * contains the area of the enlarged face - * This is only used for the ECT solver.*/ - m_fields.alloc_init(FieldType::area_mod, Direction{0}, lev, amrex::convert(ba, Bx_nodal_flag), - dm, ncomps, guard_cells.ng_FieldSolver, 0.0_rt); - m_fields.alloc_init(FieldType::area_mod, Direction{1}, lev, amrex::convert(ba, By_nodal_flag), - dm, ncomps, guard_cells.ng_FieldSolver, 0.0_rt); - m_fields.alloc_init(FieldType::area_mod, Direction{2}, lev, amrex::convert(ba, Bz_nodal_flag), - dm, ncomps, guard_cells.ng_FieldSolver, 0.0_rt); - - m_borrowing[lev][0] = std::make_unique>( + m_fields.alloc_init(FieldType::face_areas, Direction{0}, lev, + amrex::convert(ba, Bx_nodal_flag), dm, + ncomps, guard_cells.ng_FieldSolver, 0.0_rt); + m_fields.alloc_init(FieldType::face_areas, Direction{1}, lev, + amrex::convert(ba, By_nodal_flag), dm, + ncomps, guard_cells.ng_FieldSolver, 0.0_rt); + m_fields.alloc_init(FieldType::face_areas, Direction{2}, lev, + amrex::convert(ba, Bz_nodal_flag), dm, + ncomps, guard_cells.ng_FieldSolver, 0.0_rt); + + AllocInitMultiFab(m_flag_info_face[lev][0], + amrex::convert(ba, Bx_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, + "m_flag_info_face[x]"); + AllocInitMultiFab(m_flag_info_face[lev][1], + amrex::convert(ba, By_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, + "m_flag_info_face[y]"); + AllocInitMultiFab(m_flag_info_face[lev][2], + amrex::convert(ba, Bz_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, + "m_flag_info_face[z]"); + AllocInitMultiFab(m_flag_ext_face[lev][0], + amrex::convert(ba, Bx_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, + "m_flag_ext_face[x]"); + AllocInitMultiFab(m_flag_ext_face[lev][1], + amrex::convert(ba, By_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, + "m_flag_ext_face[y]"); + AllocInitMultiFab(m_flag_ext_face[lev][2], + amrex::convert(ba, Bz_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, + "m_flag_ext_face[z]"); + + /** EB: area_mod contains the modified areas of the mesh faces, + * i.e. if a face is enlarged it contains the area of the + * enlarged face This is only used for the ECT solver.*/ + m_fields.alloc_init(FieldType::area_mod, Direction{0}, lev, + amrex::convert(ba, Bx_nodal_flag), dm, + ncomps, guard_cells.ng_FieldSolver, 0.0_rt); + m_fields.alloc_init(FieldType::area_mod, Direction{1}, lev, + amrex::convert(ba, By_nodal_flag), dm, + ncomps, guard_cells.ng_FieldSolver, 0.0_rt); + m_fields.alloc_init(FieldType::area_mod, Direction{2}, lev, + amrex::convert(ba, Bz_nodal_flag), dm, + ncomps, guard_cells.ng_FieldSolver, 0.0_rt); + + m_borrowing[lev][0] = + std::make_unique>( amrex::convert(ba, Bx_nodal_flag), dm); - m_borrowing[lev][1] = std::make_unique>( + m_borrowing[lev][1] = + std::make_unique>( amrex::convert(ba, By_nodal_flag), dm); - m_borrowing[lev][2] = std::make_unique>( + m_borrowing[lev][2] = + std::make_unique>( amrex::convert(ba, Bz_nodal_flag), dm); - /** Venl contains the electromotive force for every mesh face, i.e. every entry is - * the corresponding entry in ECTRhofield multiplied by the total area (possibly with enlargement) - * This is only used for the ECT solver.*/ - m_fields.alloc_init(FieldType::Venl, Direction{0}, lev, amrex::convert(ba, Bx_nodal_flag), - dm, ncomps, guard_cells.ng_FieldSolver, 0.0_rt); - m_fields.alloc_init(FieldType::Venl, Direction{1}, lev, amrex::convert(ba, By_nodal_flag), - dm, ncomps, guard_cells.ng_FieldSolver, 0.0_rt); - m_fields.alloc_init(FieldType::Venl, Direction{2}, lev, amrex::convert(ba, Bz_nodal_flag), - dm, ncomps, guard_cells.ng_FieldSolver, 0.0_rt); + /** Venl contains the electromotive force for every mesh face, + * i.e. every entry is the corresponding entry in ECTRhofield + * multiplied by the total area (possibly with enlargement) This + * is only used for the ECT solver.*/ + m_fields.alloc_init(FieldType::Venl, Direction{0}, lev, + amrex::convert(ba, Bx_nodal_flag), dm, + ncomps, guard_cells.ng_FieldSolver, 0.0_rt); + m_fields.alloc_init(FieldType::Venl, Direction{1}, lev, + amrex::convert(ba, By_nodal_flag), dm, + ncomps, guard_cells.ng_FieldSolver, 0.0_rt); + m_fields.alloc_init(FieldType::Venl, Direction{2}, lev, + amrex::convert(ba, Bz_nodal_flag), dm, + ncomps, guard_cells.ng_FieldSolver, 0.0_rt); /** ECTRhofield is needed only by the ect - * solver and it contains the electromotive force density for every mesh face. - * The name ECTRhofield has been used to comply with the notation of the paper - * https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=4463918 (page 9, equation 4 - * and below). - * Although it's called rho it has nothing to do with the charge density! - * This is only used for the ECT solver.*/ - m_fields.alloc_init(FieldType::ECTRhofield, Direction{0}, lev, amrex::convert(ba, Bx_nodal_flag), - dm, ncomps, guard_cells.ng_FieldSolver, 0.0_rt); - m_fields.alloc_init(FieldType::ECTRhofield, Direction{1}, lev, amrex::convert(ba, By_nodal_flag), - dm, ncomps, guard_cells.ng_FieldSolver, 0.0_rt); - m_fields.alloc_init(FieldType::ECTRhofield, Direction{2}, lev, amrex::convert(ba, Bz_nodal_flag), - dm, ncomps, guard_cells.ng_FieldSolver, 0.0_rt); + * solver and it contains the electromotive force density for + * every mesh face. The name ECTRhofield has been used to comply + * with the notation of the paper + * https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=4463918 + * (page 9, equation 4 and below). Although it's called rho it + * has nothing to do with the charge density! This is only used + * for the ECT solver.*/ + m_fields.alloc_init(FieldType::ECTRhofield, Direction{0}, lev, + amrex::convert(ba, Bx_nodal_flag), dm, + ncomps, guard_cells.ng_FieldSolver, 0.0_rt); + m_fields.alloc_init(FieldType::ECTRhofield, Direction{1}, lev, + amrex::convert(ba, By_nodal_flag), dm, + ncomps, guard_cells.ng_FieldSolver, 0.0_rt); + m_fields.alloc_init(FieldType::ECTRhofield, Direction{2}, lev, + amrex::convert(ba, Bz_nodal_flag), dm, + ncomps, guard_cells.ng_FieldSolver, 0.0_rt); } } } int rho_ncomps = 0; - if( (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame) || - (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) || - (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameEffectivePotential) || - (electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC) ) { + if ((electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame) || + (electrostatic_solver_id == + ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) || + (electrostatic_solver_id == + ElectrostaticSolverAlgo::LabFrameEffectivePotential) || + (electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC)) { rho_ncomps = ncomps; } if (do_dive_cleaning) { - rho_ncomps = 2*ncomps; + rho_ncomps = 2 * ncomps; } if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { if (do_dive_cleaning || update_with_rho || current_correction) { - // For the multi-J algorithm we can allocate only one rho component (no distinction between old and new) - rho_ncomps = (WarpX::do_multi_J) ? ncomps : 2*ncomps; + // For the multi-J algorithm we can allocate only one rho component + // (no distinction between old and new) + rho_ncomps = (WarpX::do_multi_J) ? ncomps : 2 * ncomps; } } - if (rho_ncomps > 0) - { - m_fields.alloc_init(FieldType::rho_fp, - lev, amrex::convert(ba, rho_nodal_flag), dm, - rho_ncomps, ngRho, 0.0_rt); + if (rho_ncomps > 0) { + m_fields.alloc_init(FieldType::rho_fp, lev, + amrex::convert(ba, rho_nodal_flag), dm, rho_ncomps, + ngRho, 0.0_rt); } if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame || - electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic || - electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameEffectivePotential ) - { - const IntVect ngPhi = IntVect( AMREX_D_DECL(1,1,1) ); - m_fields.alloc_init(FieldType::phi_fp, lev, amrex::convert(ba, phi_nodal_flag), dm, - ncomps, ngPhi, 0.0_rt ); + electrostatic_solver_id == + ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic || + electrostatic_solver_id == + ElectrostaticSolverAlgo::LabFrameEffectivePotential) { + const IntVect ngPhi = IntVect(AMREX_D_DECL(1, 1, 1)); + m_fields.alloc_init(FieldType::phi_fp, lev, + amrex::convert(ba, phi_nodal_flag), dm, ncomps, + ngPhi, 0.0_rt); + } + + if (m_do_subcycling && lev == 0) { + m_fields.alloc_init(FieldType::current_store, Direction{0}, lev, + amrex::convert(ba, jx_nodal_flag), dm, ncomps, ngJ, + 0.0_rt); + m_fields.alloc_init(FieldType::current_store, Direction{1}, lev, + amrex::convert(ba, jy_nodal_flag), dm, ncomps, ngJ, + 0.0_rt); + m_fields.alloc_init(FieldType::current_store, Direction{2}, lev, + amrex::convert(ba, jz_nodal_flag), dm, ncomps, ngJ, + 0.0_rt); } - if (m_do_subcycling && lev == 0) - { - m_fields.alloc_init(FieldType::current_store, Direction{0}, lev, amrex::convert(ba,jx_nodal_flag), dm, ncomps, ngJ, 0.0_rt); - m_fields.alloc_init(FieldType::current_store, Direction{1}, lev, amrex::convert(ba,jy_nodal_flag), dm, ncomps, ngJ, 0.0_rt); - m_fields.alloc_init(FieldType::current_store, Direction{2}, lev, amrex::convert(ba,jz_nodal_flag), dm, ncomps, ngJ, 0.0_rt); - } - - if (do_dive_cleaning) - { - m_fields.alloc_init(FieldType::F_fp, - lev, amrex::convert(ba, F_nodal_flag), dm, - ncomps, ngF, 0.0_rt); + if (do_dive_cleaning) { + m_fields.alloc_init(FieldType::F_fp, lev, + amrex::convert(ba, F_nodal_flag), dm, ncomps, ngF, + 0.0_rt); } - if (do_divb_cleaning) - { - m_fields.alloc_init(FieldType::G_fp, - lev, amrex::convert(ba, G_nodal_flag), dm, - ncomps, ngG, 0.0_rt); + if (do_divb_cleaning) { + m_fields.alloc_init(FieldType::G_fp, lev, + amrex::convert(ba, G_nodal_flag), dm, ncomps, ngG, + 0.0_rt); } - if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) - { + if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { // Allocate and initialize the spectral solver #ifndef WARPX_USE_FFT - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( false, - "WarpX::AllocLevelMFs: PSATD solver requires WarpX build with spectral solver support."); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + false, "WarpX::AllocLevelMFs: PSATD solver requires WarpX build " + "with spectral solver support."); #else // Check whether the option periodic, single box is valid here if (fft_periodic_single_box) { -# ifdef WARPX_DIM_RZ +#ifdef WARPX_DIM_RZ WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - geom[0].isPeriodic(1) // domain is periodic in z - && ba.size() == 1 && lev == 0, // domain is decomposed in a single box - "The option `psatd.periodic_single_box_fft` can only be used for a periodic domain, decomposed in a single box"); -# else + geom[0].isPeriodic(1) // domain is periodic in z + && ba.size() == 1 && + lev == 0, // domain is decomposed in a single box + "The option `psatd.periodic_single_box_fft` can only be used " + "for a periodic domain, decomposed in a single box"); +#else WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - geom[0].isAllPeriodic() // domain is periodic in all directions - && ba.size() == 1 && lev == 0, // domain is decomposed in a single box - "The option `psatd.periodic_single_box_fft` can only be used for a periodic domain, decomposed in a single box"); -# endif + geom[0].isAllPeriodic() // domain is periodic in all directions + && ba.size() == 1 && + lev == 0, // domain is decomposed in a single box + "The option `psatd.periodic_single_box_fft` can only be used " + "for a periodic domain, decomposed in a single box"); +#endif } // Get the cell-centered box - BoxArray realspace_ba = ba; // Copy box + BoxArray realspace_ba = ba; // Copy box realspace_ba.enclosedCells(); // Make it cell-centered - // Define spectral solver -# ifdef WARPX_DIM_RZ - if ( !fft_periodic_single_box ) { + // Define spectral solver +#ifdef WARPX_DIM_RZ + if (!fft_periodic_single_box) { realspace_ba.grow(1, ngEB[1]); // add guard cells only in z } - if (field_boundary_hi[0] == FieldBoundaryType::PML && !do_pml_in_domain) { + if (field_boundary_hi[0] == FieldBoundaryType::PML && + !do_pml_in_domain) { // Extend region that is solved for to include the guard cells // which is where the PML boundary is applied. realspace_ba.growHi(0, pml_ncell); } - AllocLevelSpectralSolverRZ(spectral_solver_fp, - lev, - realspace_ba, - dm, + AllocLevelSpectralSolverRZ(spectral_solver_fp, lev, realspace_ba, dm, dx); -# else - if ( !fft_periodic_single_box ) { - realspace_ba.grow(ngEB); // add guard cells +#else + if (!fft_periodic_single_box) { + realspace_ba.grow(ngEB); // add guard cells } bool const pml_flag_false = false; - AllocLevelSpectralSolver(spectral_solver_fp, - lev, - realspace_ba, - dm, - dx, + AllocLevelSpectralSolver(spectral_solver_fp, lev, realspace_ba, dm, dx, pml_flag_false); -# endif +#endif #endif } // ElectromagneticSolverAlgo::PSATD else { - m_fdtd_solver_fp[lev] = std::make_unique(electromagnetic_solver_id, dx, grid_type); + m_fdtd_solver_fp[lev] = std::make_unique( + electromagnetic_solver_id, dx, grid_type); } // // The Aux patch (i.e., the full solution) // - if (aux_is_nodal and grid_type != GridType::Collocated) - { + if (aux_is_nodal and grid_type != GridType::Collocated) { // Create aux multifabs on Nodal Box Array - BoxArray const nba = amrex::convert(ba,IntVect::TheNodeVector()); - - m_fields.alloc_init(FieldType::Bfield_aux, Direction{0}, lev, nba, dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Bfield_aux, Direction{1}, lev, nba, dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Bfield_aux, Direction{2}, lev, nba, dm, ncomps, ngEB, 0.0_rt); - - m_fields.alloc_init(FieldType::Efield_aux, Direction{0}, lev, nba, dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Efield_aux, Direction{1}, lev, nba, dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Efield_aux, Direction{2}, lev, nba, dm, ncomps, ngEB, 0.0_rt); + BoxArray const nba = amrex::convert(ba, IntVect::TheNodeVector()); + + m_fields.alloc_init(FieldType::Bfield_aux, Direction{0}, lev, nba, dm, + ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_aux, Direction{1}, lev, nba, dm, + ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_aux, Direction{2}, lev, nba, dm, + ncomps, ngEB, 0.0_rt); + + m_fields.alloc_init(FieldType::Efield_aux, Direction{0}, lev, nba, dm, + ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Efield_aux, Direction{1}, lev, nba, dm, + ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Efield_aux, Direction{2}, lev, nba, dm, + ncomps, ngEB, 0.0_rt); } else if (lev == 0) { if (WarpX::fft_do_time_averaging) { - m_fields.alias_init(FieldType::Bfield_aux, FieldType::Bfield_avg_fp, Direction{0}, lev, 0.0_rt); - m_fields.alias_init(FieldType::Bfield_aux, FieldType::Bfield_avg_fp, Direction{1}, lev, 0.0_rt); - m_fields.alias_init(FieldType::Bfield_aux, FieldType::Bfield_avg_fp, Direction{2}, lev, 0.0_rt); - - m_fields.alias_init(FieldType::Efield_aux, FieldType::Efield_avg_fp, Direction{0}, lev, 0.0_rt); - m_fields.alias_init(FieldType::Efield_aux, FieldType::Efield_avg_fp, Direction{1}, lev, 0.0_rt); - m_fields.alias_init(FieldType::Efield_aux, FieldType::Efield_avg_fp, Direction{2}, lev, 0.0_rt); + m_fields.alias_init(FieldType::Bfield_aux, FieldType::Bfield_avg_fp, + Direction{0}, lev, 0.0_rt); + m_fields.alias_init(FieldType::Bfield_aux, FieldType::Bfield_avg_fp, + Direction{1}, lev, 0.0_rt); + m_fields.alias_init(FieldType::Bfield_aux, FieldType::Bfield_avg_fp, + Direction{2}, lev, 0.0_rt); + + m_fields.alias_init(FieldType::Efield_aux, FieldType::Efield_avg_fp, + Direction{0}, lev, 0.0_rt); + m_fields.alias_init(FieldType::Efield_aux, FieldType::Efield_avg_fp, + Direction{1}, lev, 0.0_rt); + m_fields.alias_init(FieldType::Efield_aux, FieldType::Efield_avg_fp, + Direction{2}, lev, 0.0_rt); } else { if (mypc->m_B_ext_particle_s == "read_from_file") { - m_fields.alloc_init(FieldType::Bfield_aux, Direction{0}, lev, amrex::convert(ba, Bx_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Bfield_aux, Direction{1}, lev, amrex::convert(ba, By_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Bfield_aux, Direction{2}, lev, amrex::convert(ba, Bz_nodal_flag), dm, ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_aux, Direction{0}, lev, + amrex::convert(ba, Bx_nodal_flag), dm, + ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_aux, Direction{1}, lev, + amrex::convert(ba, By_nodal_flag), dm, + ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_aux, Direction{2}, lev, + amrex::convert(ba, Bz_nodal_flag), dm, + ncomps, ngEB, 0.0_rt); } else { - // In this case, the aux grid is simply an alias of the fp grid (most common case in WarpX) - m_fields.alias_init(FieldType::Bfield_aux, FieldType::Bfield_fp, Direction{0}, lev, 0.0_rt); - m_fields.alias_init(FieldType::Bfield_aux, FieldType::Bfield_fp, Direction{1}, lev, 0.0_rt); - m_fields.alias_init(FieldType::Bfield_aux, FieldType::Bfield_fp, Direction{2}, lev, 0.0_rt); + // In this case, the aux grid is simply an alias of the fp grid + // (most common case in WarpX) + m_fields.alias_init(FieldType::Bfield_aux, FieldType::Bfield_fp, + Direction{0}, lev, 0.0_rt); + m_fields.alias_init(FieldType::Bfield_aux, FieldType::Bfield_fp, + Direction{1}, lev, 0.0_rt); + m_fields.alias_init(FieldType::Bfield_aux, FieldType::Bfield_fp, + Direction{2}, lev, 0.0_rt); } if (mypc->m_E_ext_particle_s == "read_from_file") { - m_fields.alloc_init(FieldType::Efield_aux, Direction{0}, lev, amrex::convert(ba, Ex_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Efield_aux, Direction{1}, lev, amrex::convert(ba, Ey_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Efield_aux, Direction{2}, lev, amrex::convert(ba, Ez_nodal_flag), dm, ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Efield_aux, Direction{0}, lev, + amrex::convert(ba, Ex_nodal_flag), dm, + ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Efield_aux, Direction{1}, lev, + amrex::convert(ba, Ey_nodal_flag), dm, + ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Efield_aux, Direction{2}, lev, + amrex::convert(ba, Ez_nodal_flag), dm, + ncomps, ngEB, 0.0_rt); } else { - // In this case, the aux grid is simply an alias of the fp grid (most common case in WarpX) - m_fields.alias_init(FieldType::Efield_aux, FieldType::Efield_fp, Direction{0}, lev, 0.0_rt); - m_fields.alias_init(FieldType::Efield_aux, FieldType::Efield_fp, Direction{1}, lev, 0.0_rt); - m_fields.alias_init(FieldType::Efield_aux, FieldType::Efield_fp, Direction{2}, lev, 0.0_rt); + // In this case, the aux grid is simply an alias of the fp grid + // (most common case in WarpX) + m_fields.alias_init(FieldType::Efield_aux, FieldType::Efield_fp, + Direction{0}, lev, 0.0_rt); + m_fields.alias_init(FieldType::Efield_aux, FieldType::Efield_fp, + Direction{1}, lev, 0.0_rt); + m_fields.alias_init(FieldType::Efield_aux, FieldType::Efield_fp, + Direction{2}, lev, 0.0_rt); } } } else { - m_fields.alloc_init(FieldType::Bfield_aux, Direction{0}, lev, amrex::convert(ba, Bx_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Bfield_aux, Direction{1}, lev, amrex::convert(ba, By_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Bfield_aux, Direction{2}, lev, amrex::convert(ba, Bz_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - - m_fields.alloc_init(FieldType::Efield_aux, Direction{0}, lev, amrex::convert(ba, Ex_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Efield_aux, Direction{1}, lev, amrex::convert(ba, Ey_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Efield_aux, Direction{2}, lev, amrex::convert(ba, Ez_nodal_flag), dm, ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_aux, Direction{0}, lev, + amrex::convert(ba, Bx_nodal_flag), dm, ncomps, ngEB, + 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_aux, Direction{1}, lev, + amrex::convert(ba, By_nodal_flag), dm, ncomps, ngEB, + 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_aux, Direction{2}, lev, + amrex::convert(ba, Bz_nodal_flag), dm, ncomps, ngEB, + 0.0_rt); + + m_fields.alloc_init(FieldType::Efield_aux, Direction{0}, lev, + amrex::convert(ba, Ex_nodal_flag), dm, ncomps, ngEB, + 0.0_rt); + m_fields.alloc_init(FieldType::Efield_aux, Direction{1}, lev, + amrex::convert(ba, Ey_nodal_flag), dm, ncomps, ngEB, + 0.0_rt); + m_fields.alloc_init(FieldType::Efield_aux, Direction{2}, lev, + amrex::convert(ba, Ez_nodal_flag), dm, ncomps, ngEB, + 0.0_rt); } // The external fields that are read from file - if (m_p_ext_field_params->B_ext_grid_type != ExternalFieldType::default_zero && m_p_ext_field_params->B_ext_grid_type != ExternalFieldType::constant) { - // These fields will be added directly to the grid, i.e. to fp, and need to match the index type - m_fields.alloc_init(FieldType::Bfield_fp_external, Direction{0}, lev, - amrex::convert(ba, m_fields.get(FieldType::Bfield_fp,Direction{0},lev)->ixType()), + if (m_p_ext_field_params->B_ext_grid_type != + ExternalFieldType::default_zero && + m_p_ext_field_params->B_ext_grid_type != ExternalFieldType::constant) { + // These fields will be added directly to the grid, i.e. to fp, and need + // to match the index type + m_fields.alloc_init( + FieldType::Bfield_fp_external, Direction{0}, lev, + amrex::convert( + ba, m_fields.get(FieldType::Bfield_fp, Direction{0}, lev) + ->ixType()), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Bfield_fp_external, Direction{1}, lev, - amrex::convert(ba, m_fields.get(FieldType::Bfield_fp,Direction{1},lev)->ixType()), + m_fields.alloc_init( + FieldType::Bfield_fp_external, Direction{1}, lev, + amrex::convert( + ba, m_fields.get(FieldType::Bfield_fp, Direction{1}, lev) + ->ixType()), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Bfield_fp_external, Direction{2}, lev, - amrex::convert(ba, m_fields.get(FieldType::Bfield_fp,Direction{2},lev)->ixType()), + m_fields.alloc_init( + FieldType::Bfield_fp_external, Direction{2}, lev, + amrex::convert( + ba, m_fields.get(FieldType::Bfield_fp, Direction{2}, lev) + ->ixType()), dm, ncomps, ngEB, 0.0_rt); } if (mypc->m_B_ext_particle_s == "read_from_file") { - // These fields will be added to the fields that the particles see, and need to match the index type - auto *Bfield_aux_levl_0 = m_fields.get(FieldType::Bfield_aux, Direction{0}, lev); - auto *Bfield_aux_levl_1 = m_fields.get(FieldType::Bfield_aux, Direction{1}, lev); - auto *Bfield_aux_levl_2 = m_fields.get(FieldType::Bfield_aux, Direction{2}, lev); + // These fields will be added to the fields that the particles see, and + // need to match the index type + auto* Bfield_aux_levl_0 = + m_fields.get(FieldType::Bfield_aux, Direction{0}, lev); + auto* Bfield_aux_levl_1 = + m_fields.get(FieldType::Bfield_aux, Direction{1}, lev); + auto* Bfield_aux_levl_2 = + m_fields.get(FieldType::Bfield_aux, Direction{2}, lev); // Same as Bfield_fp for reading external field data - m_fields.alloc_init(FieldType::B_external_particle_field, Direction{0}, lev, amrex::convert(ba, Bfield_aux_levl_0->ixType()), - dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::B_external_particle_field, Direction{1}, lev, amrex::convert(ba, Bfield_aux_levl_1->ixType()), - dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::B_external_particle_field, Direction{2}, lev, amrex::convert(ba, Bfield_aux_levl_2->ixType()), + m_fields.alloc_init(FieldType::B_external_particle_field, Direction{0}, + lev, + amrex::convert(ba, Bfield_aux_levl_0->ixType()), dm, + ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::B_external_particle_field, Direction{1}, + lev, + amrex::convert(ba, Bfield_aux_levl_1->ixType()), dm, + ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::B_external_particle_field, Direction{2}, + lev, + amrex::convert(ba, Bfield_aux_levl_2->ixType()), dm, + ncomps, ngEB, 0.0_rt); + } + if (m_p_ext_field_params->E_ext_grid_type != + ExternalFieldType::default_zero && + m_p_ext_field_params->E_ext_grid_type != ExternalFieldType::constant) { + // These fields will be added directly to the grid, i.e. to fp, and need + // to match the index type + m_fields.alloc_init( + FieldType::Efield_fp_external, Direction{0}, lev, + amrex::convert( + ba, m_fields.get(FieldType::Efield_fp, Direction{0}, lev) + ->ixType()), dm, ncomps, ngEB, 0.0_rt); - } - if (m_p_ext_field_params->E_ext_grid_type != ExternalFieldType::default_zero && m_p_ext_field_params->E_ext_grid_type != ExternalFieldType::constant) { - // These fields will be added directly to the grid, i.e. to fp, and need to match the index type - m_fields.alloc_init(FieldType::Efield_fp_external, Direction{0}, lev, amrex::convert(ba, m_fields.get(FieldType::Efield_fp, Direction{0}, lev)->ixType()), - dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Efield_fp_external, Direction{1}, lev, amrex::convert(ba, m_fields.get(FieldType::Efield_fp, Direction{1}, lev)->ixType()), + m_fields.alloc_init( + FieldType::Efield_fp_external, Direction{1}, lev, + amrex::convert( + ba, m_fields.get(FieldType::Efield_fp, Direction{1}, lev) + ->ixType()), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Efield_fp_external, Direction{2}, lev, amrex::convert(ba, m_fields.get(FieldType::Efield_fp, Direction{2}, lev)->ixType()), + m_fields.alloc_init( + FieldType::Efield_fp_external, Direction{2}, lev, + amrex::convert( + ba, m_fields.get(FieldType::Efield_fp, Direction{2}, lev) + ->ixType()), dm, ncomps, ngEB, 0.0_rt); } if (mypc->m_E_ext_particle_s == "read_from_file") { - // These fields will be added to the fields that the particles see, and need to match the index type - auto *Efield_aux_levl_0 = m_fields.get(FieldType::Efield_aux, Direction{0}, lev); - auto *Efield_aux_levl_1 = m_fields.get(FieldType::Efield_aux, Direction{1}, lev); - auto *Efield_aux_levl_2 = m_fields.get(FieldType::Efield_aux, Direction{2}, lev); + // These fields will be added to the fields that the particles see, and + // need to match the index type + auto* Efield_aux_levl_0 = + m_fields.get(FieldType::Efield_aux, Direction{0}, lev); + auto* Efield_aux_levl_1 = + m_fields.get(FieldType::Efield_aux, Direction{1}, lev); + auto* Efield_aux_levl_2 = + m_fields.get(FieldType::Efield_aux, Direction{2}, lev); // Same as Efield_fp for reading external field data - m_fields.alloc_init(FieldType::E_external_particle_field, Direction{0}, lev, amrex::convert(ba, Efield_aux_levl_0->ixType()), - dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::E_external_particle_field, Direction{1}, lev, amrex::convert(ba, Efield_aux_levl_1->ixType()), - dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::E_external_particle_field, Direction{2}, lev, amrex::convert(ba, Efield_aux_levl_2->ixType()), - dm, ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::E_external_particle_field, Direction{0}, + lev, + amrex::convert(ba, Efield_aux_levl_0->ixType()), dm, + ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::E_external_particle_field, Direction{1}, + lev, + amrex::convert(ba, Efield_aux_levl_1->ixType()), dm, + ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::E_external_particle_field, Direction{2}, + lev, + amrex::convert(ba, Efield_aux_levl_2->ixType()), dm, + ncomps, ngEB, 0.0_rt); } // // The coarse patch // - if (lev > 0) - { + if (lev > 0) { BoxArray cba = ba; - cba.coarsen(refRatio(lev-1)); - const std::array cdx = CellSize(lev-1); + cba.coarsen(refRatio(lev - 1)); + const std::array cdx = CellSize(lev - 1); // Create the MultiFabs for B - m_fields.alloc_init(FieldType::Bfield_cp, Direction{0}, lev, amrex::convert(cba, Bx_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Bfield_cp, Direction{1}, lev, amrex::convert(cba, By_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Bfield_cp, Direction{2}, lev, amrex::convert(cba, Bz_nodal_flag), dm, ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_cp, Direction{0}, lev, + amrex::convert(cba, Bx_nodal_flag), dm, ncomps, + ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_cp, Direction{1}, lev, + amrex::convert(cba, By_nodal_flag), dm, ncomps, + ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_cp, Direction{2}, lev, + amrex::convert(cba, Bz_nodal_flag), dm, ncomps, + ngEB, 0.0_rt); // Create the MultiFabs for E - m_fields.alloc_init(FieldType::Efield_cp, Direction{0}, lev, amrex::convert(cba, Ex_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Efield_cp, Direction{1}, lev, amrex::convert(cba, Ey_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Efield_cp, Direction{2}, lev, amrex::convert(cba, Ez_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - - if (fft_do_time_averaging) - { - m_fields.alloc_init(FieldType::Bfield_avg_cp, Direction{0}, lev, amrex::convert(cba, Bx_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Bfield_avg_cp, Direction{1}, lev, amrex::convert(cba, By_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Bfield_avg_cp, Direction{2}, lev, amrex::convert(cba, Bz_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - - m_fields.alloc_init(FieldType::Efield_avg_cp, Direction{0}, lev, amrex::convert(cba, Ex_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Efield_avg_cp, Direction{1}, lev, amrex::convert(cba, Ey_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Efield_avg_cp, Direction{2}, lev, amrex::convert(cba, Ez_nodal_flag), dm, ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Efield_cp, Direction{0}, lev, + amrex::convert(cba, Ex_nodal_flag), dm, ncomps, + ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Efield_cp, Direction{1}, lev, + amrex::convert(cba, Ey_nodal_flag), dm, ncomps, + ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Efield_cp, Direction{2}, lev, + amrex::convert(cba, Ez_nodal_flag), dm, ncomps, + ngEB, 0.0_rt); + + if (fft_do_time_averaging) { + m_fields.alloc_init(FieldType::Bfield_avg_cp, Direction{0}, lev, + amrex::convert(cba, Bx_nodal_flag), dm, ncomps, + ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_avg_cp, Direction{1}, lev, + amrex::convert(cba, By_nodal_flag), dm, ncomps, + ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_avg_cp, Direction{2}, lev, + amrex::convert(cba, Bz_nodal_flag), dm, ncomps, + ngEB, 0.0_rt); + + m_fields.alloc_init(FieldType::Efield_avg_cp, Direction{0}, lev, + amrex::convert(cba, Ex_nodal_flag), dm, ncomps, + ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Efield_avg_cp, Direction{1}, lev, + amrex::convert(cba, Ey_nodal_flag), dm, ncomps, + ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Efield_avg_cp, Direction{2}, lev, + amrex::convert(cba, Ez_nodal_flag), dm, ncomps, + ngEB, 0.0_rt); } // Create the MultiFabs for the current - m_fields.alloc_init(FieldType::current_cp, Direction{0}, lev, amrex::convert(cba, jx_nodal_flag), dm, ncomps, ngJ, 0.0_rt); - m_fields.alloc_init(FieldType::current_cp, Direction{1}, lev, amrex::convert(cba, jy_nodal_flag), dm, ncomps, ngJ, 0.0_rt); - m_fields.alloc_init(FieldType::current_cp, Direction{2}, lev, amrex::convert(cba, jz_nodal_flag), dm, ncomps, ngJ, 0.0_rt); + m_fields.alloc_init(FieldType::current_cp, Direction{0}, lev, + amrex::convert(cba, jx_nodal_flag), dm, ncomps, ngJ, + 0.0_rt); + m_fields.alloc_init(FieldType::current_cp, Direction{1}, lev, + amrex::convert(cba, jy_nodal_flag), dm, ncomps, ngJ, + 0.0_rt); + m_fields.alloc_init(FieldType::current_cp, Direction{2}, lev, + amrex::convert(cba, jz_nodal_flag), dm, ncomps, ngJ, + 0.0_rt); if (rho_ncomps > 0) { - m_fields.alloc_init(FieldType::rho_cp, - lev, amrex::convert(cba, rho_nodal_flag), dm, - rho_ncomps, ngRho, 0.0_rt); + m_fields.alloc_init(FieldType::rho_cp, lev, + amrex::convert(cba, rho_nodal_flag), dm, + rho_ncomps, ngRho, 0.0_rt); } - if (do_dive_cleaning) - { - m_fields.alloc_init(FieldType::F_cp, - lev, amrex::convert(cba, IntVect::TheUnitVector()), dm, - ncomps, ngF, 0.0_rt); + if (do_dive_cleaning) { + m_fields.alloc_init(FieldType::F_cp, lev, + amrex::convert(cba, IntVect::TheUnitVector()), + dm, ncomps, ngF, 0.0_rt); } - if (do_divb_cleaning) - { - if (grid_type == GridType::Collocated) - { - m_fields.alloc_init(FieldType::G_cp, - lev, amrex::convert(cba, IntVect::TheUnitVector()), dm, - ncomps, ngG, 0.0_rt); - } - else // grid_type=staggered or grid_type=hybrid + if (do_divb_cleaning) { + if (grid_type == GridType::Collocated) { + m_fields.alloc_init( + FieldType::G_cp, lev, + amrex::convert(cba, IntVect::TheUnitVector()), dm, ncomps, + ngG, 0.0_rt); + } else // grid_type=staggered or grid_type=hybrid { - m_fields.alloc_init(FieldType::G_cp, - lev, amrex::convert(cba, IntVect::TheZeroVector()), dm, - ncomps, ngG, 0.0_rt); + m_fields.alloc_init( + FieldType::G_cp, lev, + amrex::convert(cba, IntVect::TheZeroVector()), dm, ncomps, + ngG, 0.0_rt); } } - if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) - { + if (WarpX::electromagnetic_solver_id == + ElectromagneticSolverAlgo::PSATD) { // Allocate and initialize the spectral solver #ifndef WARPX_USE_FFT - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( false, - "WarpX::AllocLevelMFs: PSATD solver requires WarpX build with spectral solver support."); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + false, "WarpX::AllocLevelMFs: PSATD solver requires WarpX " + "build with spectral solver support."); #else // Get the cell-centered box, with guard cells - BoxArray c_realspace_ba = cba;// Copy box + BoxArray c_realspace_ba = cba; // Copy box c_realspace_ba.enclosedCells(); // Make it cell-centered - // Define spectral solver + // Define spectral solver #ifdef WARPX_DIM_RZ c_realspace_ba.grow(1, ngEB[1]); // add guard cells only in z - if (field_boundary_hi[0] == FieldBoundaryType::PML && !do_pml_in_domain) { + if (field_boundary_hi[0] == FieldBoundaryType::PML && + !do_pml_in_domain) { // Extend region that is solved for to include the guard cells // which is where the PML boundary is applied. c_realspace_ba.growHi(0, pml_ncell); } - AllocLevelSpectralSolverRZ(spectral_solver_cp, - lev, - c_realspace_ba, - dm, - cdx); -# else + AllocLevelSpectralSolverRZ(spectral_solver_cp, lev, c_realspace_ba, + dm, cdx); +#else c_realspace_ba.grow(ngEB); bool const pml_flag_false = false; - AllocLevelSpectralSolver(spectral_solver_cp, - lev, - c_realspace_ba, - dm, - cdx, - pml_flag_false); -# endif + AllocLevelSpectralSolver(spectral_solver_cp, lev, c_realspace_ba, + dm, cdx, pml_flag_false); +#endif #endif } // ElectromagneticSolverAlgo::PSATD else { - m_fdtd_solver_cp[lev] = std::make_unique(electromagnetic_solver_id, cdx, - grid_type); + m_fdtd_solver_cp[lev] = std::make_unique( + electromagnetic_solver_id, cdx, grid_type); } } // // Copy of the coarse aux // - if (lev > 0 && (n_field_gather_buffer > 0 || n_current_deposition_buffer > 0 || - mypc->nSpeciesGatherFromMainGrid() > 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)); + cba.coarsen(refRatio(lev - 1)); - if (n_field_gather_buffer > 0 || mypc->nSpeciesGatherFromMainGrid() > 0) { + if (n_field_gather_buffer > 0 || + mypc->nSpeciesGatherFromMainGrid() > 0) { if (aux_is_nodal) { - BoxArray const& cnba = amrex::convert(cba,IntVect::TheNodeVector()); + BoxArray const& cnba = + amrex::convert(cba, IntVect::TheNodeVector()); // Create the MultiFabs for B - m_fields.alloc_init(FieldType::Bfield_cax, Direction{0}, lev, cnba, dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Bfield_cax, Direction{1}, lev, cnba, dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Bfield_cax, Direction{2}, lev, cnba, dm, ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_cax, Direction{0}, lev, + cnba, dm, ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_cax, Direction{1}, lev, + cnba, dm, ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_cax, Direction{2}, lev, + cnba, dm, ncomps, ngEB, 0.0_rt); // Create the MultiFabs for E - m_fields.alloc_init(FieldType::Efield_cax, Direction{0}, lev, cnba, dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Efield_cax, Direction{1}, lev, cnba, dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Efield_cax, Direction{2}, lev, cnba, dm, ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Efield_cax, Direction{0}, lev, + cnba, dm, ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Efield_cax, Direction{1}, lev, + cnba, dm, ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Efield_cax, Direction{2}, lev, + cnba, dm, ncomps, ngEB, 0.0_rt); } else { // Create the MultiFabs for B - m_fields.alloc_init(FieldType::Bfield_cax, Direction{0}, lev, amrex::convert(cba, Bx_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Bfield_cax, Direction{1}, lev, amrex::convert(cba, By_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Bfield_cax, Direction{2}, lev, amrex::convert(cba, Bz_nodal_flag), dm, ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_cax, Direction{0}, lev, + amrex::convert(cba, Bx_nodal_flag), dm, + ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_cax, Direction{1}, lev, + amrex::convert(cba, By_nodal_flag), dm, + ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Bfield_cax, Direction{2}, lev, + amrex::convert(cba, Bz_nodal_flag), dm, + ncomps, ngEB, 0.0_rt); // Create the MultiFabs for E - m_fields.alloc_init(FieldType::Efield_cax, Direction{0}, lev, amrex::convert(cba,Ex_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Efield_cax, Direction{1}, lev, amrex::convert(cba,Ey_nodal_flag), dm, ncomps, ngEB, 0.0_rt); - m_fields.alloc_init(FieldType::Efield_cax, Direction{2}, lev, amrex::convert(cba,Ez_nodal_flag), dm, ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Efield_cax, Direction{0}, lev, + amrex::convert(cba, Ex_nodal_flag), dm, + ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Efield_cax, Direction{1}, lev, + amrex::convert(cba, Ey_nodal_flag), dm, + ncomps, ngEB, 0.0_rt); + m_fields.alloc_init(FieldType::Efield_cax, Direction{2}, lev, + amrex::convert(cba, Ez_nodal_flag), dm, + ncomps, ngEB, 0.0_rt); } - AllocInitMultiFab(gather_buffer_masks[lev], ba, dm, ncomps, amrex::IntVect(1), lev, "gather_buffer_masks"); + AllocInitMultiFab(gather_buffer_masks[lev], ba, dm, ncomps, + amrex::IntVect(1), lev, "gather_buffer_masks"); // Gather buffer masks have 1 ghost cell, because of the fact - // that particles may move by more than one cell when using subcycling. + // that particles may move by more than one cell when using + // subcycling. } if (n_current_deposition_buffer > 0) { - m_fields.alloc_init(FieldType::current_buf, Direction{0}, lev, amrex::convert(cba,jx_nodal_flag), dm, ncomps, ngJ, 0.0_rt); - m_fields.alloc_init(FieldType::current_buf, Direction{1}, lev, amrex::convert(cba,jy_nodal_flag), dm, ncomps, ngJ, 0.0_rt); - m_fields.alloc_init(FieldType::current_buf, Direction{2}, lev, amrex::convert(cba,jz_nodal_flag), dm, ncomps, ngJ, 0.0_rt); + m_fields.alloc_init(FieldType::current_buf, Direction{0}, lev, + amrex::convert(cba, jx_nodal_flag), dm, ncomps, + ngJ, 0.0_rt); + m_fields.alloc_init(FieldType::current_buf, Direction{1}, lev, + amrex::convert(cba, jy_nodal_flag), dm, ncomps, + ngJ, 0.0_rt); + m_fields.alloc_init(FieldType::current_buf, Direction{2}, lev, + amrex::convert(cba, jz_nodal_flag), dm, ncomps, + ngJ, 0.0_rt); if (m_fields.has(FieldType::rho_cp, lev)) { - m_fields.alloc_init(FieldType::rho_buf, lev, amrex::convert(cba,rho_nodal_flag), dm, 2*ncomps, ngRho, 0.0_rt); + m_fields.alloc_init(FieldType::rho_buf, lev, + amrex::convert(cba, rho_nodal_flag), dm, + 2 * ncomps, ngRho, 0.0_rt); } - AllocInitMultiFab(current_buffer_masks[lev], ba, dm, ncomps, amrex::IntVect(1), lev, "current_buffer_masks"); + AllocInitMultiFab(current_buffer_masks[lev], ba, dm, ncomps, + amrex::IntVect(1), lev, "current_buffer_masks"); // Current buffer masks have 1 ghost cell, because of the fact - // that particles may move by more than one cell when using subcycling. + // that particles may move by more than one cell when using + // subcycling. } } - if (load_balance_intervals.isActivated()) - { + if (load_balance_intervals.isActivated()) { costs[lev] = std::make_unique>(ba, dm); load_balance_efficiency[lev] = -1; } } #ifdef WARPX_USE_FFT -# ifdef WARPX_DIM_RZ +#ifdef WARPX_DIM_RZ /* \brief Allocate spectral Maxwell solver (RZ dimensions) at a level * - * \param[in, out] spectral_solver Vector of pointer to SpectralSolver, to point to allocated spectral Maxwell - * solver at a given level - * \param[in] lev Level at which to allocate spectral Maxwell solver - * \param[in] realspace_ba Box array that corresponds to the decomposition of the fields in real space - * (cell-centered; includes guard cells) - * \param[in] dm Indicates which MPI proc owns which box, in realspace_ba + * \param[in, out] spectral_solver Vector of pointer to SpectralSolver, to point + * to allocated spectral Maxwell solver at a given level + * \param[in] lev Level at which to allocate spectral Maxwell + * solver + * \param[in] realspace_ba Box array that corresponds to the + * decomposition of the fields in real space (cell-centered; includes guard + * cells) + * \param[in] dm Indicates which MPI proc owns which box, in + * realspace_ba * \param[in] dx Cell size along each dimension */ -void WarpX::AllocLevelSpectralSolverRZ (amrex::Vector>& spectral_solver, - const int lev, - const amrex::BoxArray& realspace_ba, - const amrex::DistributionMapping& dm, - const std::array& dx) -{ +void +WarpX::AllocLevelSpectralSolverRZ ( + amrex::Vector>& spectral_solver, + const int lev, const amrex::BoxArray& realspace_ba, + const amrex::DistributionMapping& dm, const std::array& dx) { const RealVect dx_vect(dx[0], dx[2]); amrex::Real solver_dt = dt[lev]; - if (WarpX::do_multi_J) { solver_dt /= static_cast(WarpX::do_multi_J_n_depositions); } + if (WarpX::do_multi_J) { + solver_dt /= static_cast(WarpX::do_multi_J_n_depositions); + } if (evolve_scheme == EvolveScheme::StrangImplicitSpectralEM) { // The step is Strang split into two half steps solver_dt /= 2.; } - auto pss = std::make_unique(lev, - realspace_ba, - dm, - n_rz_azimuthal_modes, - noz_fft, - grid_type, - m_v_galilean, - dx_vect, - solver_dt, - ::isAnyBoundaryPML(field_boundary_lo, field_boundary_hi), - update_with_rho, - fft_do_time_averaging, - J_in_time, - rho_in_time, - do_dive_cleaning, - do_divb_cleaning); + auto pss = std::make_unique( + lev, realspace_ba, dm, n_rz_azimuthal_modes, noz_fft, grid_type, + m_v_galilean, dx_vect, solver_dt, + ::isAnyBoundaryPML(field_boundary_lo, field_boundary_hi), + update_with_rho, fft_do_time_averaging, J_in_time, rho_in_time, + do_dive_cleaning, do_divb_cleaning); spectral_solver[lev] = std::move(pss); if (use_kspace_filter) { @@ -2827,25 +3187,28 @@ void WarpX::AllocLevelSpectralSolverRZ (amrex::Vector>& spectral_solver, - const int lev, - const amrex::BoxArray& realspace_ba, - const amrex::DistributionMapping& dm, - const std::array& dx, - const bool pml_flag) -{ +void +WarpX::AllocLevelSpectralSolver ( + amrex::Vector>& spectral_solver, + const int lev, const amrex::BoxArray& realspace_ba, + const amrex::DistributionMapping& dm, const std::array& dx, + const bool pml_flag) { #if defined(WARPX_DIM_3D) const RealVect dx_vect(dx[0], dx[1], dx[2]); #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) @@ -2855,139 +3218,133 @@ void WarpX::AllocLevelSpectralSolver (amrex::Vector(WarpX::do_multi_J_n_depositions); } + if (WarpX::do_multi_J) { + solver_dt /= static_cast(WarpX::do_multi_J_n_depositions); + } if (evolve_scheme == EvolveScheme::StrangImplicitSpectralEM) { // The step is Strang split into two half steps solver_dt /= 2.; } - auto pss = std::make_unique(lev, - realspace_ba, - dm, - nox_fft, - noy_fft, - noz_fft, - grid_type, - m_v_galilean, - m_v_comoving, - dx_vect, - solver_dt, - pml_flag, - fft_periodic_single_box, - update_with_rho, - fft_do_time_averaging, - m_psatd_solution_type, - J_in_time, - rho_in_time, - do_dive_cleaning, - do_divb_cleaning); + auto pss = std::make_unique( + lev, realspace_ba, dm, nox_fft, noy_fft, noz_fft, grid_type, + m_v_galilean, m_v_comoving, dx_vect, solver_dt, pml_flag, + fft_periodic_single_box, update_with_rho, fft_do_time_averaging, + m_psatd_solution_type, J_in_time, rho_in_time, do_dive_cleaning, + do_divb_cleaning); spectral_solver[lev] = std::move(pss); } -# endif +#endif #endif -std::array -WarpX::CellSize (int lev) -{ +std::array +WarpX::CellSize (int lev) { const amrex::Geometry& gm = GetInstance().Geom(lev); const Real* dx = gm.CellSize(); #if defined(WARPX_DIM_3D) - return { dx[0], dx[1], dx[2] }; + return {dx[0], dx[1], dx[2]}; #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - return { dx[0], 1.0, dx[1] }; + return {dx[0], 1.0, dx[1]}; #else - return { 1.0, 1.0, dx[0] }; + return {1.0, 1.0, dx[0]}; #endif } amrex::XDim3 -WarpX::InvCellSize (int lev) -{ - std::array dx = WarpX::CellSize(lev); - return {1._rt/dx[0], 1._rt/dx[1], 1._rt/dx[2]}; +WarpX::InvCellSize (int lev) { + std::array dx = WarpX::CellSize(lev); + return {1._rt / dx[0], 1._rt / dx[1], 1._rt / dx[2]}; } amrex::RealBox -WarpX::getRealBox(const Box& bx, int lev) -{ +WarpX::getRealBox (const Box& bx, int lev) { const amrex::Geometry& gm = GetInstance().Geom(lev); const RealBox grid_box{bx, gm.CellSize(), gm.ProbLo()}; - return( grid_box ); + return (grid_box); } amrex::XDim3 -WarpX::LowerCorner(const Box& bx, const int lev, const amrex::Real time_shift_delta) -{ - auto & warpx = GetInstance(); - const RealBox grid_box = getRealBox( bx, lev ); +WarpX::LowerCorner (const Box& bx, const int lev, + const amrex::Real time_shift_delta) { + auto& warpx = GetInstance(); + const RealBox grid_box = getRealBox(bx, lev); const Real* grid_min = grid_box.lo(); const amrex::Real cur_time = warpx.gett_new(lev); - const amrex::Real time_shift = (cur_time + time_shift_delta - warpx.time_of_last_gal_shift); - amrex::Array galilean_shift = { warpx.m_v_galilean[0]*time_shift, - warpx.m_v_galilean[1]*time_shift, - warpx.m_v_galilean[2]*time_shift }; + const amrex::Real time_shift = + (cur_time + time_shift_delta - warpx.time_of_last_gal_shift); + amrex::Array galilean_shift = { + warpx.m_v_galilean[0] * time_shift, warpx.m_v_galilean[1] * time_shift, + warpx.m_v_galilean[2] * time_shift}; #if defined(WARPX_DIM_3D) - return { grid_min[0] + galilean_shift[0], grid_min[1] + galilean_shift[1], grid_min[2] + galilean_shift[2] }; + return {grid_min[0] + galilean_shift[0], grid_min[1] + galilean_shift[1], + grid_min[2] + galilean_shift[2]}; #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - return { grid_min[0] + galilean_shift[0], std::numeric_limits::lowest(), grid_min[1] + galilean_shift[2] }; + return {grid_min[0] + galilean_shift[0], + std::numeric_limits::lowest(), + grid_min[1] + galilean_shift[2]}; #elif defined(WARPX_DIM_1D_Z) - return { std::numeric_limits::lowest(), std::numeric_limits::lowest(), grid_min[0] + galilean_shift[2] }; + return {std::numeric_limits::lowest(), + std::numeric_limits::lowest(), + grid_min[0] + galilean_shift[2]}; #endif } amrex::XDim3 -WarpX::UpperCorner(const Box& bx, const int lev, const amrex::Real time_shift_delta) -{ - auto & warpx = GetInstance(); - const RealBox grid_box = getRealBox( bx, lev ); +WarpX::UpperCorner (const Box& bx, const int lev, + const amrex::Real time_shift_delta) { + auto& warpx = GetInstance(); + const RealBox grid_box = getRealBox(bx, lev); const Real* grid_max = grid_box.hi(); const amrex::Real cur_time = warpx.gett_new(lev); - const amrex::Real time_shift = (cur_time + time_shift_delta - warpx.time_of_last_gal_shift); - amrex::Array galilean_shift = { warpx.m_v_galilean[0]*time_shift, - warpx.m_v_galilean[1]*time_shift, - warpx.m_v_galilean[2]*time_shift }; + const amrex::Real time_shift = + (cur_time + time_shift_delta - warpx.time_of_last_gal_shift); + amrex::Array galilean_shift = { + warpx.m_v_galilean[0] * time_shift, warpx.m_v_galilean[1] * time_shift, + warpx.m_v_galilean[2] * time_shift}; #if defined(WARPX_DIM_3D) - return { grid_max[0] + galilean_shift[0], grid_max[1] + galilean_shift[1], grid_max[2] + galilean_shift[2] }; + return {grid_max[0] + galilean_shift[0], grid_max[1] + galilean_shift[1], + grid_max[2] + galilean_shift[2]}; #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - return { grid_max[0] + galilean_shift[0], std::numeric_limits::max(), grid_max[1] + galilean_shift[1] }; + return {grid_max[0] + galilean_shift[0], std::numeric_limits::max(), + grid_max[1] + galilean_shift[1]}; #elif defined(WARPX_DIM_1D_Z) - return { std::numeric_limits::max(), std::numeric_limits::max(), grid_max[0] + galilean_shift[0] }; + return {std::numeric_limits::max(), std::numeric_limits::max(), + grid_max[0] + galilean_shift[0]}; #endif } IntVect -WarpX::RefRatio (int lev) -{ +WarpX::RefRatio (int lev) { return GetInstance().refRatio(lev); } void WarpX::ComputeDivB (amrex::MultiFab& divB, int const dcomp, ablastr::fields::VectorField const& B, - const std::array& dx) -{ + const std::array& dx) { ComputeDivB(divB, dcomp, B, dx, IntVect::TheZeroVector()); } void WarpX::ComputeDivB (amrex::MultiFab& divB, int const dcomp, ablastr::fields::VectorField const& B, - const std::array& dx, IntVect const ngrow) -{ - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(grid_type != GridType::Collocated, + const std::array& dx, IntVect const ngrow) { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + grid_type != GridType::Collocated, "ComputeDivB not implemented with warpx.grid_type=collocated."); - const Real dxinv = 1._rt/dx[0], dyinv = 1._rt/dx[1], dzinv = 1._rt/dx[2]; + const Real dxinv = 1._rt / dx[0], dyinv = 1._rt / dx[1], + dzinv = 1._rt / dx[2]; #ifdef WARPX_DIM_RZ const Real rmin = GetInstance().Geom(0).ProbLo(0); @@ -2996,47 +3353,46 @@ WarpX::ComputeDivB (amrex::MultiFab& divB, int const dcomp, #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - for (MFIter mfi(divB, TilingIfNotGPU()); mfi.isValid(); ++mfi) - { + for (MFIter mfi(divB, TilingIfNotGPU()); mfi.isValid(); ++mfi) { const Box bx = mfi.growntilebox(ngrow); amrex::Array4 const& Bxfab = B[0]->array(mfi); amrex::Array4 const& Byfab = B[1]->array(mfi); amrex::Array4 const& Bzfab = B[2]->array(mfi); amrex::Array4 const& divBfab = divB.array(mfi); - ParallelFor(bx, - [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - warpx_computedivb(i, j, k, dcomp, divBfab, Bxfab, Byfab, Bzfab, dxinv, dyinv, dzinv + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + warpx_computedivb(i, j, k, dcomp, divBfab, Bxfab, Byfab, Bzfab, + dxinv, dyinv, dzinv #ifdef WARPX_DIM_RZ - ,rmin + , + rmin #endif - ); + ); }); } } void -WarpX::ComputeDivE(amrex::MultiFab& divE, const int lev) -{ - if ( WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD ) { +WarpX::ComputeDivE (amrex::MultiFab& divE, const int lev) { + if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { #ifdef WARPX_USE_FFT - const ablastr::fields::VectorField Efield_aux_lev = m_fields.get_alldirs(FieldType::Efield_aux, lev); + const ablastr::fields::VectorField Efield_aux_lev = + m_fields.get_alldirs(FieldType::Efield_aux, lev); spectral_solver_fp[lev]->ComputeSpectralDivE(lev, Efield_aux_lev, divE); #else WARPX_ABORT_WITH_MESSAGE( "ComputeDivE: PSATD requested but not compiled"); #endif } else { - const ablastr::fields::VectorField Efield_aux_lev = m_fields.get_alldirs(FieldType::Efield_aux, lev); + const ablastr::fields::VectorField Efield_aux_lev = + m_fields.get_alldirs(FieldType::Efield_aux, lev); m_fdtd_solver_fp[lev]->ComputeDivE(Efield_aux_lev, divE); } } #if (defined WARPX_DIM_RZ) && (defined WARPX_USE_FFT) PML_RZ* -WarpX::GetPML_RZ (int lev) -{ +WarpX::GetPML_RZ (int lev) { if (pml_rz[lev]) { // This should check if pml was initialized return pml_rz[lev].get(); @@ -3047,8 +3403,7 @@ WarpX::GetPML_RZ (int lev) #endif PML* -WarpX::GetPML (int lev) -{ +WarpX::GetPML (int lev) { if (do_pml) { // This should check if pml was initialized return pml[lev].get(); @@ -3057,87 +3412,75 @@ WarpX::GetPML (int lev) } } -std::vector< bool > -WarpX::getPMLdirections() const -{ - std::vector< bool > dirsWithPML( 6, false ); -#if AMREX_SPACEDIM!=3 - dirsWithPML.resize( 4 ); +std::vector +WarpX::getPMLdirections () const { + std::vector dirsWithPML(6, false); +#if AMREX_SPACEDIM != 3 + dirsWithPML.resize(4); #endif - if( do_pml ) - { - for( int i = 0; i < static_cast(dirsWithPML.size()) / 2; ++i ) - { - dirsWithPML.at( 2u*i ) = bool(do_pml_Lo[0][i]); // on level 0 - dirsWithPML.at( 2u*i + 1u ) = bool(do_pml_Hi[0][i]); // on level 0 + if (do_pml) { + for (int i = 0; i < static_cast(dirsWithPML.size()) / 2; ++i) { + dirsWithPML.at(2u * i) = bool(do_pml_Lo[0][i]); // on level 0 + dirsWithPML.at(2u * i + 1u) = bool(do_pml_Hi[0][i]); // on level 0 } } return dirsWithPML; } amrex::LayoutData* -WarpX::getCosts (int lev) -{ - if (m_instance) - { +WarpX::getCosts (int lev) { + if (m_instance) { return m_instance->costs[lev].get(); - } else - { + } else { return nullptr; } } void -WarpX::setLoadBalanceEfficiency (const int lev, const amrex::Real efficiency) -{ - if (m_instance) - { +WarpX::setLoadBalanceEfficiency (const int lev, const amrex::Real efficiency) { + if (m_instance) { m_instance->load_balance_efficiency[lev] = efficiency; - } else - { + } else { return; } } amrex::Real -WarpX::getLoadBalanceEfficiency (const int lev) -{ - if (m_instance) - { +WarpX::getLoadBalanceEfficiency (const int lev) { + if (m_instance) { return m_instance->load_balance_efficiency[lev]; - } else - { + } else { return -1; } } void -WarpX::BuildBufferMasks () -{ - for (int lev = 1; lev <= maxLevel(); ++lev) - { - for (int ipass = 0; ipass < 2; ++ipass) - { - const int ngbuffer = (ipass == 0) ? n_current_deposition_buffer : n_field_gather_buffer; - iMultiFab* bmasks = (ipass == 0) ? current_buffer_masks[lev].get() : gather_buffer_masks[lev].get(); - if (bmasks) - { +WarpX::BuildBufferMasks () { + for (int lev = 1; lev <= maxLevel(); ++lev) { + for (int ipass = 0; ipass < 2; ++ipass) { + const int ngbuffer = (ipass == 0) ? n_current_deposition_buffer + : n_field_gather_buffer; + iMultiFab* bmasks = (ipass == 0) ? current_buffer_masks[lev].get() + : gather_buffer_masks[lev].get(); + if (bmasks) { const IntVect ngtmp = ngbuffer + bmasks->nGrowVect(); - iMultiFab tmp(bmasks->boxArray(), bmasks->DistributionMap(), 1, ngtmp); + iMultiFab tmp(bmasks->boxArray(), bmasks->DistributionMap(), 1, + ngtmp); const int covered = 1; const int notcovered = 0; const int physbnd = 1; const int interior = 1; const Box& dom = Geom(lev).Domain(); const amrex::Periodicity& period = Geom(lev).periodicity(); - tmp.BuildMask(dom, period, covered, notcovered, physbnd, interior); + tmp.BuildMask(dom, period, covered, notcovered, physbnd, + interior); #ifdef AMREX_USE_OMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif - for (MFIter mfi(*bmasks, true); mfi.isValid(); ++mfi) - { + for (MFIter mfi(*bmasks, true); mfi.isValid(); ++mfi) { const Box tbx = mfi.growntilebox(); - BuildBufferMasksInBox( tbx, (*bmasks)[mfi], tmp[mfi], ngbuffer ); + BuildBufferMasksInBox(tbx, (*bmasks)[mfi], tmp[mfi], + ngbuffer); } } } @@ -3153,36 +3496,35 @@ WarpX::BuildBufferMasks () * \param ng Number of guard cells */ void -WarpX::BuildBufferMasksInBox ( const amrex::Box tbx, amrex::IArrayBox &buffer_mask, - const amrex::IArrayBox &guard_mask, const int ng ) -{ +WarpX::BuildBufferMasksInBox (const amrex::Box tbx, + amrex::IArrayBox& buffer_mask, + const amrex::IArrayBox& guard_mask, + const int ng) { auto const& msk = buffer_mask.array(); auto const& gmsk = guard_mask.const_array(); const amrex::Dim3 ng3 = amrex::IntVect(ng).dim3(); - amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - for (int kk = k-ng3.z; kk <= k+ng3.z; ++kk) { - for (int jj = j-ng3.y; jj <= j+ng3.y; ++jj) { - for (int ii = i-ng3.x; ii <= i+ng3.x; ++ii) { - if ( gmsk(ii,jj,kk) == 0 ) { - msk(i,j,k) = 0; + amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + for (int kk = k - ng3.z; kk <= k + ng3.z; ++kk) { + for (int jj = j - ng3.y; jj <= j + ng3.y; ++jj) { + for (int ii = i - ng3.x; ii <= i + ng3.x; ++ii) { + if (gmsk(ii, jj, kk) == 0) { + msk(i, j, k) = 0; return; } } } } - msk(i,j,k) = 1; + msk(i, j, k) = 1; }); } -void WarpX::AllocateCenteringCoefficients (amrex::Gpu::DeviceVector& device_centering_stencil_coeffs_x, - amrex::Gpu::DeviceVector& device_centering_stencil_coeffs_y, - amrex::Gpu::DeviceVector& device_centering_stencil_coeffs_z, - const int centering_nox, - const int centering_noy, - const int centering_noz, - ablastr::utils::enums::GridType a_grid_type) -{ +void +WarpX::AllocateCenteringCoefficients ( + amrex::Gpu::DeviceVector& device_centering_stencil_coeffs_x, + amrex::Gpu::DeviceVector& device_centering_stencil_coeffs_y, + amrex::Gpu::DeviceVector& device_centering_stencil_coeffs_z, + const int centering_nox, const int centering_noy, const int centering_noz, + ablastr::utils::enums::GridType a_grid_type) { // Vectors of Fornberg stencil coefficients amrex::Vector Fornberg_stencil_coeffs_x; amrex::Vector Fornberg_stencil_coeffs_y; @@ -3193,9 +3535,12 @@ void WarpX::AllocateCenteringCoefficients (amrex::Gpu::DeviceVector amrex::Vector host_centering_stencil_coeffs_y; amrex::Vector host_centering_stencil_coeffs_z; - Fornberg_stencil_coeffs_x = ablastr::math::getFornbergStencilCoefficients(centering_nox, a_grid_type); - Fornberg_stencil_coeffs_y = ablastr::math::getFornbergStencilCoefficients(centering_noy, a_grid_type); - Fornberg_stencil_coeffs_z = ablastr::math::getFornbergStencilCoefficients(centering_noz, a_grid_type); + Fornberg_stencil_coeffs_x = ablastr::math::getFornbergStencilCoefficients( + centering_nox, a_grid_type); + Fornberg_stencil_coeffs_y = ablastr::math::getFornbergStencilCoefficients( + centering_noy, a_grid_type); + Fornberg_stencil_coeffs_z = ablastr::math::getFornbergStencilCoefficients( + centering_noz, a_grid_type); host_centering_stencil_coeffs_x.resize(centering_nox); host_centering_stencil_coeffs_y.resize(centering_noy); @@ -3203,37 +3548,34 @@ void WarpX::AllocateCenteringCoefficients (amrex::Gpu::DeviceVector // Re-order Fornberg stencil coefficients: // example for order 6: (c_0,c_1,c_2) becomes (c_2,c_1,c_0,c_0,c_1,c_2) - ablastr::math::ReorderFornbergCoefficients( - host_centering_stencil_coeffs_x, - Fornberg_stencil_coeffs_x, - centering_nox - ); - ablastr::math::ReorderFornbergCoefficients( - host_centering_stencil_coeffs_y, - Fornberg_stencil_coeffs_y, - centering_noy - ); - ablastr::math::ReorderFornbergCoefficients( - host_centering_stencil_coeffs_z, - Fornberg_stencil_coeffs_z, - centering_noz - ); + ablastr::math::ReorderFornbergCoefficients(host_centering_stencil_coeffs_x, + Fornberg_stencil_coeffs_x, + centering_nox); + ablastr::math::ReorderFornbergCoefficients(host_centering_stencil_coeffs_y, + Fornberg_stencil_coeffs_y, + centering_noy); + ablastr::math::ReorderFornbergCoefficients(host_centering_stencil_coeffs_z, + Fornberg_stencil_coeffs_z, + centering_noz); // Device vectors of stencil coefficients used for finite-order centering - device_centering_stencil_coeffs_x.resize(host_centering_stencil_coeffs_x.size()); + device_centering_stencil_coeffs_x.resize( + host_centering_stencil_coeffs_x.size()); amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, host_centering_stencil_coeffs_x.begin(), host_centering_stencil_coeffs_x.end(), device_centering_stencil_coeffs_x.begin()); - device_centering_stencil_coeffs_y.resize(host_centering_stencil_coeffs_y.size()); + device_centering_stencil_coeffs_y.resize( + host_centering_stencil_coeffs_y.size()); amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, host_centering_stencil_coeffs_y.begin(), host_centering_stencil_coeffs_y.end(), device_centering_stencil_coeffs_y.begin()); - device_centering_stencil_coeffs_z.resize(host_centering_stencil_coeffs_z.size()); + device_centering_stencil_coeffs_z.resize( + host_centering_stencil_coeffs_z.size()); amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, host_centering_stencil_coeffs_z.begin(), host_centering_stencil_coeffs_z.end(), @@ -3243,33 +3585,32 @@ void WarpX::AllocateCenteringCoefficients (amrex::Gpu::DeviceVector } const iMultiFab* -WarpX::CurrentBufferMasks (int lev) -{ +WarpX::CurrentBufferMasks (int lev) { return GetInstance().getCurrentBufferMasks(lev); } const iMultiFab* -WarpX::GatherBufferMasks (int lev) -{ +WarpX::GatherBufferMasks (int lev) { return GetInstance().getGatherBufferMasks(lev); } void -WarpX::StoreCurrent (int lev) -{ +WarpX::StoreCurrent (int lev) { using ablastr::fields::Direction; for (int idim = 0; idim < 3; ++idim) { - if (m_fields.has(FieldType::current_store, Direction{idim},lev)) { - MultiFab::Copy(*m_fields.get(FieldType::current_store, Direction{idim}, lev), - *m_fields.get(FieldType::current_fp, Direction{idim}, lev), - 0, 0, 1, m_fields.get(FieldType::current_store, Direction{idim}, lev)->nGrowVect()); + if (m_fields.has(FieldType::current_store, Direction{idim}, lev)) { + MultiFab::Copy( + *m_fields.get(FieldType::current_store, Direction{idim}, lev), + *m_fields.get(FieldType::current_fp, Direction{idim}, lev), 0, + 0, 1, + m_fields.get(FieldType::current_store, Direction{idim}, lev) + ->nGrowVect()); } } } void -WarpX::RestoreCurrent (int lev) -{ +WarpX::RestoreCurrent (int lev) { using ablastr::fields::Direction; using warpx::fields::FieldType; @@ -3277,35 +3618,36 @@ WarpX::RestoreCurrent (int lev) if (m_fields.has(FieldType::current_store, Direction{idim}, lev)) { std::swap( *m_fields.get(FieldType::current_fp, Direction{idim}, lev), - *m_fields.get(FieldType::current_store, Direction{idim}, lev) - ); + *m_fields.get(FieldType::current_store, Direction{idim}, lev)); } } } bool -WarpX::isAnyParticleBoundaryThermal () -{ +WarpX::isAnyParticleBoundaryThermal () { for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - if (WarpX::particle_boundary_lo[idim] == ParticleBoundaryType::Thermal) {return true;} - if (WarpX::particle_boundary_hi[idim] == ParticleBoundaryType::Thermal) {return true;} + if (WarpX::particle_boundary_lo[idim] == + ParticleBoundaryType::Thermal) { + return true; + } + if (WarpX::particle_boundary_hi[idim] == + ParticleBoundaryType::Thermal) { + return true; + } } return false; } void -WarpX::AllocInitMultiFab ( - std::unique_ptr& mf, - const amrex::BoxArray& ba, - const amrex::DistributionMapping& dm, - const int ncomp, - const amrex::IntVect& ngrow, - const int level, - const std::string& name, - std::optional initial_value) -{ +WarpX::AllocInitMultiFab (std::unique_ptr& mf, + const amrex::BoxArray& ba, + const amrex::DistributionMapping& dm, const int ncomp, + const amrex::IntVect& ngrow, const int level, + const std::string& name, + std::optional initial_value) { // Add the suffix "[level=level]" - const auto name_with_suffix = name + "[level=" + std::to_string(level) + "]"; + const auto name_with_suffix = + name + "[level=" + std::to_string(level) + "]"; const auto tag = amrex::MFInfo().SetTag(name_with_suffix); mf = std::make_unique(ba, dm, ncomp, ngrow, tag); if (initial_value) { @@ -3315,8 +3657,7 @@ WarpX::AllocInitMultiFab ( } amrex::DistributionMapping -WarpX::MakeDistributionMap (int lev, amrex::BoxArray const& ba) -{ +WarpX::MakeDistributionMap (int lev, amrex::BoxArray const& ba) { bool roundrobin_sfc = false; const ParmParse pp("warpx"); pp.query("roundrobin_sfc", roundrobin_sfc); @@ -3341,25 +3682,29 @@ WarpX::MakeDistributionMap (int lev, amrex::BoxArray const& ba) } const amrex::iMultiFab* -WarpX::getFieldDotMaskPointer ( FieldType field_type, int lev, ablastr::fields::Direction dir ) const -{ +WarpX::getFieldDotMaskPointer (FieldType field_type, int lev, + ablastr::fields::Direction dir) const { const auto periodicity = Geom(lev).periodicity(); - switch(field_type) - { - case FieldType::Efield_fp : - ::SetDotMask( Efield_dotMask[lev][dir], m_fields.get("Efield_fp", dir, lev), periodicity); - return Efield_dotMask[lev][dir].get(); - case FieldType::Bfield_fp : - ::SetDotMask( Bfield_dotMask[lev][dir], m_fields.get("Bfield_fp", dir, lev), periodicity); - return Bfield_dotMask[lev][dir].get(); - case FieldType::vector_potential_fp : - ::SetDotMask( Afield_dotMask[lev][dir], m_fields.get("vector_potential_fp", dir, lev), periodicity); - return Afield_dotMask[lev][dir].get(); - case FieldType::phi_fp : - ::SetDotMask( phi_dotMask[lev], m_fields.get("phi_fp", dir, lev), periodicity); - return phi_dotMask[lev].get(); - default: - WARPX_ABORT_WITH_MESSAGE("Invalid field type for dotMask"); - return Efield_dotMask[lev][dir].get(); + switch (field_type) { + case FieldType::Efield_fp: + ::SetDotMask(Efield_dotMask[lev][dir], + m_fields.get("Efield_fp", dir, lev), periodicity); + return Efield_dotMask[lev][dir].get(); + case FieldType::Bfield_fp: + ::SetDotMask(Bfield_dotMask[lev][dir], + m_fields.get("Bfield_fp", dir, lev), periodicity); + return Bfield_dotMask[lev][dir].get(); + case FieldType::vector_potential_fp: + ::SetDotMask(Afield_dotMask[lev][dir], + m_fields.get("vector_potential_fp", dir, lev), + periodicity); + return Afield_dotMask[lev][dir].get(); + case FieldType::phi_fp: + ::SetDotMask(phi_dotMask[lev], m_fields.get("phi_fp", dir, lev), + periodicity); + return phi_dotMask[lev].get(); + default: + WARPX_ABORT_WITH_MESSAGE("Invalid field type for dotMask"); + return Efield_dotMask[lev][dir].get(); } }