diff --git a/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp b/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp index 556b8f8fca4..1b66d7317af 100644 --- a/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp +++ b/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp @@ -161,13 +161,21 @@ void WarpX::HybridPICEvolveFields () } } - // Calculate the electron pressure at t=n+1 - m_hybrid_pic_model->CalculateElectronPressure(); - // Update the E field to t=n+1 using the extrapolated J_i^n+1 value m_hybrid_pic_model->CalculateCurrentAmpere( m_fields.get_mr_levels_alldirs(FieldType::Bfield_fp, finest_level), m_fields.get_mr_levels_alldirs(FieldType::edge_lengths, finest_level)); + + // Update Ue in electron fluid container + 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 + + // Calculate the electron pressure at t=n+1 + m_hybrid_pic_model->CalculateElectronPressure(); + m_hybrid_pic_model->HybridPICSolveE( m_fields.get_mr_levels_alldirs(FieldType::Efield_fp, finest_level), current_fp_temp, diff --git a/Source/Fluids/WarpXFluidContainer.H b/Source/Fluids/WarpXFluidContainer.H index 8b21195cdfa..a5223d941dd 100644 --- a/Source/Fluids/WarpXFluidContainer.H +++ b/Source/Fluids/WarpXFluidContainer.H @@ -131,6 +131,9 @@ 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 + void HybridUpdateUe (ablastr::fields::MultiFabRegister& m_fields, 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 d1870274a9f..3d3624dcae2 100644 --- a/Source/Fluids/WarpXFluidContainer.cpp +++ b/Source/Fluids/WarpXFluidContainer.cpp @@ -1406,3 +1406,46 @@ void WarpXFluidContainer::DepositCurrent( ); } } + +void WarpXFluidContainer::HybridUpdateUe (ablastr::fields::MultiFabRegister& fields, int lev) +{ + using ablastr::fields::Direction; + using warpx::fields::FieldType; + + ablastr::fields::ScalarField rho_fp = fields.get(FieldType::rho_fp, lev); + ablastr::fields::VectorField current_fp_ampere = fields.get_alldirs(FieldType::hybrid_current_fp_ampere, lev); + ablastr::fields::VectorField current_fp = fields.get_alldirs(FieldType::current_fp, lev); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*rho_fp, TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + Array4 const& rho = rho_fp->array(mfi); + + amrex::Array4 const & Jx = current_fp_ampere[0]->array(mfi); + amrex::Array4 const & Jy = current_fp_ampere[1]->array(mfi); + amrex::Array4 const & Jz = current_fp_ampere[2]->array(mfi); + + amrex::Array4 const & Jix = current_fp[0]->array(mfi); + amrex::Array4 const & Jiy = current_fp[1]->array(mfi); + amrex::Array4 const & Jiz = current_fp[2]->array(mfi); + + const amrex::Array4 Uex = fields.get(name_mf_NU, Direction{0}, lev)->array(mfi); + const amrex::Array4 Uey = fields.get(name_mf_NU, Direction{1}, lev)->array(mfi); + const amrex::Array4 Uez = fields.get(name_mf_NU, Direction{2}, 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 ){ + Uex(i, j, k) = (Jx(i, j, k) - Jix(i, j, k))/rho(i, j, k); + Uey(i, j, k) = (Jy(i, j, k) - Jiy(i, j, k))/rho(i, j, k); + Uez(i, j, k) = (Jz(i, j, k) - Jiz(i, j, k))/rho(i, j, k); + } + + }); + } + + +} \ No newline at end of file