/* Copyright 2023 The WarpX Community * * This file is part of WarpX. * * Authors: Roelof Groenewald (TAE Technologies) * * License: BSD-3-Clause-LBNL */ #ifndef WARPX_HYBRIDPICMODEL_H_ #define WARPX_HYBRIDPICMODEL_H_ #include "HybridPICModel_fwd.H" #include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H" #include "Utils/Parser/ParserUtils.H" #include "Utils/WarpXConst.H" #include "Utils/WarpXProfilerWrapper.H" #include "WarpX.H" #include #include /** * \brief This class contains the parameters needed to evaluate hybrid field * solutions (kinetic ions with fluid electrons). */ class HybridPICModel { public: HybridPICModel (int nlevs_max); // constructor /** Read user-defined model parameters. Called in constructor. */ void ReadParameters (); /** Allocate hybrid-PIC specific multifabs. Called in constructor. */ void AllocateMFs (int nlevs_max); void AllocateLevelMFs (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm, const int ncomps, const amrex::IntVect& ngJ, const amrex::IntVect& ngRho, const amrex::IntVect& jx_nodal_flag, const amrex::IntVect& jy_nodal_flag, const amrex::IntVect& jz_nodal_flag, const amrex::IntVect& rho_nodal_flag); /** Helper function to clear values from hybrid-PIC specific multifabs. */ void ClearLevel (int lev); void InitData (); /** * \brief * Function to calculate the total current based on Ampere's law while * neglecting displacement current (J = curl x B). Used in the Ohm's law * solver (kinetic-fluid hybrid model). * * \param[in] Bfield Magnetic field from which the current is calculated. * \param[in] edge_lengths Length of cell edges taking embedded boundaries into account */ void CalculateCurrentAmpere ( amrex::Vector, 3>> const& Bfield, amrex::Vector, 3>> const& edge_lengths ); void CalculateCurrentAmpere ( std::array< std::unique_ptr, 3> const& Bfield, std::array< std::unique_ptr, 3> const& edge_lengths, const int lev ); /** * \brief * Function to update the E-field using Ohm's law (hybrid-PIC model). */ void HybridPICSolveE ( amrex::Vector, 3>>& Efield, amrex::Vector, 3>> const& Jfield, amrex::Vector, 3>> const& Bfield, amrex::Vector> const& rhofield, amrex::Vector, 3>> const& edge_lengths, DtType dt_type); void HybridPICSolveE ( std::array< std::unique_ptr, 3>& Efield, std::array< std::unique_ptr, 3> const& Jfield, std::array< std::unique_ptr, 3> const& Bfield, std::unique_ptr const& rhofield, std::array< std::unique_ptr, 3> const& edge_lengths, const int lev, DtType dt_type); void HybridPICSolveE ( std::array< std::unique_ptr, 3>& Efield, std::array< std::unique_ptr, 3> const& Jfield, std::array< std::unique_ptr, 3> const& Bfield, std::unique_ptr const& rhofield, std::array< std::unique_ptr, 3> const& edge_lengths, const int lev, PatchType patch_type, DtType dt_type); /** * \brief * Function to calculate the electron pressure at a given timestep type * using the simulation charge density. Used in the Ohm's law solver * (kinetic-fluid hybrid model). */ void CalculateElectronPressure ( DtType a_dt_type); void CalculateElectronPressure (const int lev, DtType a_dt_type); /** * \brief Fill the electron pressure multifab given the kinetic particle * charge density (and assumption of quasi-neutrality) using the user * specified electron equation of state. * * \param[out] Pe scalar electron pressure MultiFab at a given level * \param[in] rhofield scalar ion chrge density Multifab at a given level */ void FillElectronPressureMF ( std::unique_ptr const& Pe, amrex::MultiFab* const& rhofield ); // Declare variables to hold hybrid-PIC model parameters /** Number of substeps to take when evolving B */ int m_substeps = 100; /** Electron temperature in eV */ amrex::Real m_elec_temp; /** Reference electron density */ amrex::Real m_n0_ref = 1.0; /** Electron pressure scaling exponent */ amrex::Real m_gamma = 5.0/3.0; /** Plasma density floor - if n < n_floor it will be set to n_floor */ amrex::Real m_n_floor = 1.0; /** Plasma resistivity */ std::string m_eta_expression = "0.0"; std::unique_ptr m_resistivity_parser; amrex::ParserExecutor<1> m_eta; // Declare multifabs specifically needed for the hybrid-PIC model amrex::Vector< std::unique_ptr > rho_fp_temp; amrex::Vector, 3 > > current_fp_temp; amrex::Vector, 3 > > current_fp_ampere; amrex::Vector< std::unique_ptr > electron_pressure_fp; // Helper functions to retrieve hybrid-PIC multifabs amrex::MultiFab * get_pointer_current_fp_ampere (int lev, int direction) const { return current_fp_ampere[lev][direction].get(); } amrex::MultiFab * get_pointer_electron_pressure_fp (int lev) const { return electron_pressure_fp[lev].get(); } /** Gpu Vector with index type of the Jx multifab */ amrex::GpuArray Jx_IndexType; /** Gpu Vector with index type of the Jy multifab */ amrex::GpuArray Jy_IndexType; /** Gpu Vector with index type of the Jz multifab */ amrex::GpuArray Jz_IndexType; /** Gpu Vector with index type of the Bx multifab */ amrex::GpuArray Bx_IndexType; /** Gpu Vector with index type of the By multifab */ amrex::GpuArray By_IndexType; /** Gpu Vector with index type of the Bz multifab */ amrex::GpuArray Bz_IndexType; /** Gpu Vector with index type of the Ex multifab */ amrex::GpuArray Ex_IndexType; /** Gpu Vector with index type of the Ey multifab */ amrex::GpuArray Ey_IndexType; /** Gpu Vector with index type of the Ez multifab */ amrex::GpuArray Ez_IndexType; }; /** * \brief * This struct contains only static functions to compute the electron pressure * using the particle density at a given point and the user provided reference * density and temperatures. */ struct ElectronPressure { AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real get_pressure (amrex::Real const n0, amrex::Real const T0, amrex::Real const gamma, amrex::Real const rho) { return n0 * T0 * pow((rho/PhysConst::q_e)/n0, gamma); } }; #endif // WARPX_HYBRIDPICMODEL_H_