Hari/gradlimiter (#5)

* adding potential gradient limiter --------- Co-authored-by: Hariswaran Sitaraman <[email protected]> Co-authored-by: Hariswaran Sitaraman <[email protected]>
NREL · Oct 25, 2024 · 5d1a50a · 5d1a50a
1 parent d5d859e
commit 5d1a50a
Show file tree

Hide file tree

Showing 6 changed files with 112 additions and 2 deletions.
diff --git a/Source/Evolve.cpp b/Source/Evolve.cpp
@@ -202,6 +202,20 @@ void Vidyut::Evolve()
                 }
             }
 
+            if(efield_limiter)
+            {
+                potential_gradlimiter(Sborder); 
+
+                //fillpatching here to get the latest efields 
+                //in sborder so that it can be used in drift vel calcs
+                //may be there is a clever way to improve performance 
+                for(int lev=0;lev<=finest_level;lev++)
+                {
+                    Sborder[lev].setVal(0.0);
+                    FillPatch(lev, cur_time+dt_common, Sborder[lev], 0, Sborder[lev].nComp());
+                }
+            }
+
             // Calculate the reactive source terms for all species/levels
             if(do_reactions)
             {

diff --git a/Source/HelperFuncs.H b/Source/HelperFuncs.H
@@ -90,4 +90,41 @@ amrex::Real get_applied_potential(Real current_time, int domain_end, int vprof,
 
     return voltage;
 }
+AMREX_GPU_DEVICE
+    AMREX_FORCE_INLINE
+amrex::Real minmod_limiter(amrex::Real left_difference,amrex::Real right_difference)
+{
+    amrex::Real phi=1.0;
+    if(amrex::Math::abs(right_difference) > 0.0)
+    {
+        amrex::Real r = left_difference/right_difference;
+        phi = std::max(0.0,std::min(1.0,r));
+    }
+    else
+    {
+        phi = 1.0;
+    }
+    return(phi);
+}
+
+AMREX_GPU_DEVICE
+AMREX_FORCE_INLINE
+amrex::Real efieldlimiter(int i,int j,int k, int sweepdir,
+                          Array4<Real> const& state_arr)
+{
+    IntVect iv{AMREX_D_DECL(i,j,k)};
+    IntVect ivl{AMREX_D_DECL(i,j,k)};
+    IntVect ivr{AMREX_D_DECL(i,j,k)};
+
+    ivl[sweepdir]-=1;
+    ivr[sweepdir]+=1;
+
+    amrex::Real ldiff=state_arr(iv,POT_ID)-state_arr(ivl,POT_ID);
+    amrex::Real rdiff=state_arr(ivr,POT_ID)-state_arr(iv,POT_ID);
+
+    amrex::Real lim=minmod_limiter(ldiff,rdiff);
+
+    return(lim);
+}
+
 #endif
diff --git a/Source/PotentialSolve.cpp b/Source/PotentialSolve.cpp
@@ -570,3 +570,59 @@ void Vidyut::update_cs_technique_potential()
         }
     }
 }
+
+void Vidyut::potential_gradlimiter(Vector<MultiFab>& Sborder)
+{
+    for (int ilev = 0; ilev <= finest_level; ilev++)
+    {
+        auto prob_lo = geom[ilev].ProbLoArray();
+        auto prob_hi = geom[ilev].ProbHiArray();
+        const auto dx = geom[ilev].CellSizeArray();
+        const int* domlo_arr = geom[ilev].Domain().loVect();
+        const int* domhi_arr = geom[ilev].Domain().hiVect();
+
+        GpuArray<int,AMREX_SPACEDIM> domlo={AMREX_D_DECL(domlo_arr[0], domlo_arr[1], domlo_arr[2])};
+        GpuArray<int,AMREX_SPACEDIM> domhi={AMREX_D_DECL(domhi_arr[0], domhi_arr[1], domhi_arr[2])};
+
+        for (MFIter mfi(phi_new[ilev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
+        {
+            const Box& bx = mfi.tilebox();
+            Box bx_x = convert(bx, {AMREX_D_DECL(1, 0, 0)});
+#if AMREX_SPACEDIM > 1
+            Box bx_y = convert(bx, {AMREX_D_DECL(0, 1, 0)});
+#if AMREX_SPACEDIM == 3
+            Box bx_z = convert(bx, {0, 0, 1});
+#endif
+#endif
+            Array4<Real> phi_arr   = phi_new[ilev].array(mfi);
+            Array4<Real> sb_arr   = Sborder[ilev].array(mfi);
+
+            //x sweep
+            amrex::ParallelFor(bx_x, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+
+                amrex::Real lim=efieldlimiter(i, j, k, 0, sb_arr);
+                phi_arr(i,j,k,EFX_ID)*=lim;
+            });
+
+#if AMREX_SPACEDIM > 1
+            //y sweep
+            amrex::ParallelFor(bx_y, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+
+                amrex::Real lim=efieldlimiter(i, j, k, 1, sb_arr);
+                phi_arr(i,j,k,EFY_ID)*=lim;
+
+            });
+
+#if AMREX_SPACEDIM == 3
+            //z sweep
+            amrex::ParallelFor(bx_z, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+
+                amrex::Real lim=efieldlimiter(i, j, k, 2, sb_arr);
+                phi_arr(i,j,k,EFZ_ID)*=lim;
+
+            });
+#endif
+#endif
+        }
+    }
+}
diff --git a/Source/ScalarSolve.cpp b/Source/ScalarSolve.cpp
@@ -152,7 +152,7 @@ void Vidyut::update_rxnsrc_at_all_levels(Vector<MultiFab>& Sborder,
                 // Get molar production rates
                 CKWC(captured_gastemp, spec_C, spec_wdot, Te, EN, &ener_exch);
 
-                // Convert from mol/cm3-s to 1/m3-s and add to scalar react source MF
+                // Convert from mol/m3-s to 1/m3-s and add to scalar react source MF
                 for(int sp = 0; sp<NUM_SPECIES; sp++) rxn_arr(i,j,k,sp) = spec_wdot[sp] * N_A;
                 rxn_arr(i,j,k,NUM_SPECIES) = ener_exch;
 

diff --git a/Source/Vidyut.H b/Source/Vidyut.H
@@ -244,6 +244,7 @@ public:
 
     void update_cc_efields(amrex::Vector<MultiFab>& Sborder); 
     void update_cs_technique_potential(); 
+    void potential_gradlimiter(amrex::Vector<MultiFab>& Sborder);
 
 
     // how often each level regrids the higher levels of refinement
@@ -269,8 +270,9 @@ public:
 
     // Electron species index
     int E_IDX = 0;
+    int efield_limiter=0;
 
-    // Moitor file options
+    // Monitor file options
     int monitor_file_int = -1;
 
 private:

diff --git a/Source/Vidyut.cpp b/Source/Vidyut.cpp
@@ -364,6 +364,7 @@ void Vidyut::ReadParameters()
 
         pp.query("monitor_file_int", monitor_file_int);
         pp.query("num_timestep_correctors",num_timestep_correctors);
+        pp.query("efield_limiter",efield_limiter);
 
         pp.query("cs_technique",cs_technique);
         if(cs_technique)