diff --git a/Python/pywarpx/fields.py b/Python/pywarpx/fields.py index 0100f64f261..b82977e0ab4 100644 --- a/Python/pywarpx/fields.py +++ b/Python/pywarpx/fields.py @@ -1030,3 +1030,45 @@ def GCPPMLWrapper(level=0, include_ghosts=False): return _MultiFABWrapper( mf_name="pml_G_cp", level=level, include_ghosts=include_ghosts ) + + +# remove +def Ne_HybridWrapper(level=0, include_ghosts=False): + return _MultiFABWrapper( + mf_name="fluid_density_electrons_hybrid", level=level, include_ghosts=include_ghosts + ) + + +def NUex_HybridWrapper(level=0, include_ghosts=False): + return _MultiFABWrapper( + mf_name="fluid_momentum_density_electrons_hybrid", idir=0, level=level, include_ghosts=include_ghosts + ) + + +def NUey_HybridWrapper(level=0, include_ghosts=False): + return _MultiFABWrapper( + mf_name="fluid_momentum_density_electrons_hybrid", idir=1, level=level, include_ghosts=include_ghosts + ) + + +def NUez_HybridWrapper(level=0, include_ghosts=False): + return _MultiFABWrapper( + mf_name="fluid_momentum_density_electrons_hybrid", idir=2, level=level, include_ghosts=include_ghosts + ) + + +def Te_HybridWrapper(level=0, include_ghosts=False): + return _MultiFABWrapper( + mf_name="fluid_temperature_electrons_hybrid", level=level, include_ghosts=include_ghosts + ) + + +def Ke_HybridWrapper(level=0, include_ghosts=False): + return _MultiFABWrapper( + mf_name="fluid_entropy_electrons_hybrid", level=level, include_ghosts=include_ghosts + ) + +def Hybrid_Pe_Wrapper(level=0, include_ghosts=False): + return _MultiFABWrapper( + mf_name="hybrid_electron_pressure_fp", level=level, include_ghosts=include_ghosts + ) \ No newline at end of file diff --git a/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp b/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp index 52843ac0aea..1b51c0cc6fa 100644 --- a/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp +++ b/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp @@ -166,15 +166,24 @@ void WarpX::HybridPICEvolveFields () m_fields.get_mr_levels_alldirs(FieldType::Bfield_fp, finest_level), m_fields.get_mr_levels_alldirs(FieldType::edge_lengths, finest_level)); + // all the qdsmc solver functions should be in a ElectronEnergyEquationSolver class as well as other solvers like Layer method if(m_hybrid_pic_model->m_solve_electron_energy_equation){ + // Initialize Ke = ne^(1-gamma)*Te + // Move up, shouldn't Ke be defined with both ne and Te at n instead of ne at n+1 ? + hybrid_electron_fl->HybridInitializeKe(m_fields, m_hybrid_pic_model->m_gamma , finest_level); + // Update Ue in electron fluid container // check at what step this should be calculated, here Ue is at n+1 + // maybe move up to be consistent with Ke initialization hybrid_electron_fl->HybridUpdateUe(m_fields, finest_level); // Solve electron energy equation without sources or sinks // and update electron temperature before calculating Pe - // hybrid_electron_fl->SolveEEQqdsmc(); // should take rho, Je + // hybrid_electron_fl->SolveEEEqQDSMC(); // should take rho, Je + + // Update Te after QDSMC solver is used + hybrid_electron_fl->HybridUpdateTe(m_fields, m_hybrid_pic_model->m_gamma, finest_level); } diff --git a/Source/Fluids/WarpXFluidContainer.H b/Source/Fluids/WarpXFluidContainer.H index 0d12adf9711..dbcc768039a 100644 --- a/Source/Fluids/WarpXFluidContainer.H +++ b/Source/Fluids/WarpXFluidContainer.H @@ -131,9 +131,16 @@ public: */ void DepositCharge (ablastr::fields::MultiFabRegister& m_fields, amrex::MultiFab &rho, int lev, int icomp = 0); - // Function that calculates Ue from Je = J_ampere - Ji and saves in NU for hybrid solver + // HybridUpdateUe calculates Ue from Je = J_ampere - Ji and saves in NU for hybrid solver void HybridUpdateUe (ablastr::fields::MultiFabRegister& m_fields, int lev); + // Sets Ke in fluid container MF to be used by qdsmc particles + void HybridInitializeKe (ablastr::fields::MultiFabRegister& m_fields, amrex::Real gamma, int lev); + + // UpdateTe takes the already interpolated Ke values from the grid after pushing the fictitious particles + // and calculates Te at n+1 + void HybridUpdateTe (ablastr::fields::MultiFabRegister& m_fields, amrex::Real gamma, int lev); + [[nodiscard]] amrex::Real getCharge () const {return charge;} [[nodiscard]] amrex::Real getMass () const {return mass;} diff --git a/Source/Fluids/WarpXFluidContainer.cpp b/Source/Fluids/WarpXFluidContainer.cpp index 131ac6b2ad5..e8bf41262e2 100644 --- a/Source/Fluids/WarpXFluidContainer.cpp +++ b/Source/Fluids/WarpXFluidContainer.cpp @@ -1447,6 +1447,54 @@ void WarpXFluidContainer::HybridUpdateUe (ablastr::fields::MultiFabRegister& fie }); } +} + +void WarpXFluidContainer::HybridInitializeKe (ablastr::fields::MultiFabRegister& fields, amrex::Real gamma, int lev) +{ + using warpx::fields::FieldType; + ablastr::fields::ScalarField rho_fp = fields.get(FieldType::rho_fp, lev); +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*rho_fp, TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + amrex::Array4 const & rho = rho_fp->array(mfi); + amrex::Array4 const & Te = fields.get(name_mf_T, lev)->array(mfi); + const amrex::Array4 Ke = fields.get(name_mf_K, lev)->array(mfi); + + const Box& tilebox = mfi.tilebox(); + + ParallelFor(tilebox, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + if( rho(i, j, k) > 0.0_rt ){ + Ke(i, j, k) = Te(i, j, k)*std::pow((rho(i, j, k)/PhysConst::q_e), 1-gamma)/PhysConst::q_e; + } + }); + } } + + +void WarpXFluidContainer::HybridUpdateTe (ablastr::fields::MultiFabRegister& fields, amrex::Real gamma, int lev) +{ + using warpx::fields::FieldType; + ablastr::fields::ScalarField rho_fp = fields.get(FieldType::rho_fp, lev); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*rho_fp, TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + amrex::Array4 const & rho = rho_fp->array(mfi); + amrex::Array4 const & Ke = fields.get(name_mf_K, lev)->array(mfi); + const amrex::Array4 Te = fields.get(name_mf_T, lev)->array(mfi); + + const Box& tilebox = mfi.tilebox(); + + ParallelFor(tilebox, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + if( rho(i, j, k) > 0.0_rt ){ + Te(i, j, k) = Ke(i, j, k)*std::pow((rho(i, j, k)/PhysConst::q_e), gamma-1)*PhysConst::q_e; + } + }); + } +} \ No newline at end of file