added HybridInitializeKe and HybridUpdateTe functions to fluidcontain…

…er and calls to WarpXPushFieldsHybridPIC.cpp
marcoacc95 · Sep 24, 2024 · a7df10a · a7df10a
1 parent 75f2892
commit a7df10a
Show file tree

Hide file tree

Showing 4 changed files with 108 additions and 2 deletions.
diff --git 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
+    )
diff --git 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
@@ -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
@@ -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<amrex::Real> const & rho = rho_fp->array(mfi);
+            amrex::Array4<amrex::Real> const & Te = fields.get(name_mf_T, lev)->array(mfi);
+            const amrex::Array4<amrex::Real> 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<amrex::Real> const & rho = rho_fp->array(mfi);
+            amrex::Array4<amrex::Real> const & Ke = fields.get(name_mf_K, lev)->array(mfi);
+            const amrex::Array4<amrex::Real> 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;
+                }
+            });
+        }
+}