added WarpXFluidContainer::HybridUpdateUe and called in WarpXPushFiel…

…dsHybridPIC between CalculateCurrentAmpere and HybridPICSolveE

moved CalculateElectronPressure after HybridUpdateUe

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,

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;}
 

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<Real> const& rho = rho_fp->array(mfi);
+
+            amrex::Array4<amrex::Real> const & Jx = current_fp_ampere[0]->array(mfi);
+            amrex::Array4<amrex::Real> const & Jy = current_fp_ampere[1]->array(mfi);
+            amrex::Array4<amrex::Real> const & Jz = current_fp_ampere[2]->array(mfi);
+
+            amrex::Array4<amrex::Real> const & Jix = current_fp[0]->array(mfi);
+            amrex::Array4<amrex::Real> const & Jiy = current_fp[1]->array(mfi);
+            amrex::Array4<amrex::Real> const & Jiz = current_fp[2]->array(mfi);
+
+            const amrex::Array4<amrex::Real> Uex = fields.get(name_mf_NU, Direction{0}, lev)->array(mfi);
+            const amrex::Array4<amrex::Real> Uey = fields.get(name_mf_NU, Direction{1}, lev)->array(mfi);
+            const amrex::Array4<amrex::Real> 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);
+                }
+
+            });
+        }
+
+
+}