Merge branch 'main' of github.com:AMReX-FHD/FHDeX into main

exec/DSMC/main_driver.cpp

@@ -415,7 +415,7 @@ void main_driver(const char* argv)
 			writePlotFile(cuInst,cuMeans,cuVars,primInst,primMeans,primVars,
 				coVars,spatialCross1D,particles,geom,time,ncross,istep);
 
-			structFactPrim.WritePlotFile(istep,time,geom,"plt_SF_prim");
+			structFactPrim.WritePlotFile(istep,time,"plt_SF_prim");
 		}
 
         if ((n_steps_skip > 0 && istep == n_steps_skip) || (n_steps_skip < 0 && istep%n_steps_skip == 0) ) {

exec/DSMC_granular/main_driver.cpp

@@ -409,7 +409,7 @@ void main_driver(const char* argv)
 			struct_fact_int > 0 && plot_int > 0 &&
 			istep%plot_int == 0)
 		{
-			structFactPrim.WritePlotFile(istep,time,geom,"plt_SF_prim");
+			structFactPrim.WritePlotFile(istep,time,"plt_SF_prim");
 		}
 
 		//////////////////////////////////////

exec/cellbdytest_new/main_driver.cpp

@@ -1819,8 +1819,8 @@ void main_driver(const char* argv)
 
             // plot structure factor on plot_int
             if (istep%plot_int == 0) {
-                structFact_charge.WritePlotFile(istep,time,geomP,"plt_SF_charge");
-                structFact_vel   .WritePlotFile(istep,time,geom ,"plt_SF_vel");
+                structFact_charge.WritePlotFile(istep,time,"plt_SF_charge");
+                structFact_vel   .WritePlotFile(istep,time,"plt_SF_vel");
             }
         }
 

exec/hydro/main_driver.cpp

@@ -351,8 +351,6 @@ void main_driver(const char* argv)
 
     StructFact structFactFlattened;
 
-    Geometry geom_flat;
-
     if(project_dir >= 0){
       MultiFab Flattened;  // flattened multifab defined below
 
@@ -364,30 +362,9 @@ void main_driver(const char* argv)
       } else {
           ExtractSlice(structFactMF, Flattened, project_dir, slicepoint, 0, 1);
       }
+
       BoxArray ba_flat = Flattened.boxArray();
       const DistributionMapping& dmap_flat = Flattened.DistributionMap();
-      {
-        Box domain_flat = ba_flat.minimalBox();
-
-        // This defines the physical box
-        // we retain prob_lo and prob_hi in all directions except project_dir,
-        // where the physical size is 0 to dx[project_dir]
-        Vector<Real> projected_lo(AMREX_SPACEDIM);
-        Vector<Real> projected_hi(AMREX_SPACEDIM);
-
-        for (int d=0; d<AMREX_SPACEDIM; ++d) {
-            projected_lo[d] = prob_lo[d];
-            projected_hi[d] = prob_hi[d];
-        }
-        projected_lo[project_dir] = 0.;
-        projected_hi[project_dir] = (prob_hi[project_dir] - prob_lo[project_dir]) / n_cells[project_dir];
-
-        RealBox real_box_flat({AMREX_D_DECL(projected_lo[0],projected_lo[1],projected_lo[2])},
-                              {AMREX_D_DECL(projected_hi[0],projected_hi[1],projected_hi[2])});
-
-        // This defines a Geometry object
-        geom_flat.define(domain_flat,&real_box_flat,CoordSys::cartesian,is_periodic.data());
-      }
 
       structFactFlattened.define(ba_flat,dmap_flat,var_names,var_scaling);
     }
@@ -445,9 +422,9 @@ void main_driver(const char* argv)
         if (plot_int > 0) {
             WritePlotFile(0,time,geom,umac,pres);
             if (n_steps_skip == 0 && struct_fact_int > 0) {
-                structFact.WritePlotFile(0,0.,geom,"plt_SF");
+                structFact.WritePlotFile(0,0.,"plt_SF");
                 if(project_dir >= 0) {
-                    structFactFlattened.WritePlotFile(0,time,geom_flat,"plt_SF_Flattened");
+                    structFactFlattened.WritePlotFile(0,time,"plt_SF_Flattened");
                 }
             }
         }
@@ -547,9 +524,9 @@ void main_driver(const char* argv)
 
             // write out structure factor to plotfile
             if (step > n_steps_skip && struct_fact_int > 0) {
-                structFact.WritePlotFile(step,time,geom,"plt_SF");
+                structFact.WritePlotFile(step,time,"plt_SF");
                 if(project_dir >= 0) {
-                    structFactFlattened.WritePlotFile(step,time,geom_flat,"plt_SF_Flattened");
+                    structFactFlattened.WritePlotFile(step,time,"plt_SF_Flattened");
                 }
             }
 

exec/immersedIons/main_driver.cpp

@@ -1589,8 +1589,8 @@ void main_driver(const char* argv)
 
             // plot structure factor on plot_int
             if (istep%plot_int == 0) {
-                structFact_charge.WritePlotFile(istep,time,geomP,"plt_SF_charge");
-                structFact_vel   .WritePlotFile(istep,time,geom ,"plt_SF_vel");
+                structFact_charge.WritePlotFile(istep,time,"plt_SF_charge");
+                structFact_vel   .WritePlotFile(istep,time,"plt_SF_vel");
             }
         }
 

exec/multispec/main_driver.cpp

@@ -508,7 +508,7 @@ void main_driver(const char* argv)
         if (plot_int > 0) {
             WritePlotFile(0,0.,geom,umac,rhotot_old,rho_old,pi,charge_old,Epot);
             if (n_steps_skip == 0 && struct_fact_int > 0) {
-                structFact.WritePlotFile(0,0.,geom,"plt_SF");
+                structFact.WritePlotFile(0,0.,"plt_SF");
             }
         }
 
@@ -593,7 +593,7 @@ void main_driver(const char* argv)
         if (plot_int > 0 && istep%plot_int == 0) {
             WritePlotFile(istep,time,geom,umac,rhotot_new,rho_new,pi,charge_new,Epot);
             if (istep > n_steps_skip && struct_fact_int > 0) {
-                structFact.WritePlotFile(istep,time,geom,"plt_SF");
+                structFact.WritePlotFile(istep,time,"plt_SF");
             }
         }
 

exec/structFactTest/main_driver.cpp

@@ -150,7 +150,7 @@ void main_driver(const char* argv)
 
     structFact.FortStructure(struct_cc);
 
-    structFact.WritePlotFile(0,0.,geom,"plt_SF");
+    structFact.WritePlotFile(0,0.,"plt_SF");
 
     // Call the timer again and compute the maximum difference between the start time
     // and stop time over all processors

src_MFsurfchem/MFsurfchem_functions.H

@@ -20,8 +20,8 @@ void InitializeMFSurfchemNamespace();
 
 void init_surfcov(MultiFab& surfcov, const amrex::Geometry& geom);
 
-void sample_MFsurfchem(MultiFab& cu, MultiFab& prim, MultiFab& surfcov, MultiFab& dNadsdes, const amrex::Geometry& geom, const amrex::Real dt);
+void sample_MFsurfchem(MultiFab& cu, MultiFab& prim, MultiFab& surfcov, MultiFab& dNadsdes, MultiFab& dNads, MultiFab& dNdes, const amrex::Geometry& geom, const amrex::Real dt);
 
-void update_MFsurfchem(MultiFab& cu, MultiFab& prim, MultiFab& surfcov, MultiFab& dNadsdes, const amrex::Geometry& geom);
+void update_MFsurfchem(MultiFab& cu, MultiFab& prim, MultiFab& surfcov, MultiFab& dNadsdes, MultiFab& dNads, MultiFab& dNdes, const amrex::Geometry& geom);
 
 #endif

src_MFsurfchem/MFsurfchem_functions.cpp

@@ -2,6 +2,7 @@
 #include "AMReX_ParmParse.H"
 
 AMREX_GPU_MANAGED int MFsurfchem::n_ads_spec;
+AMREX_GPU_MANAGED int MFsurfchem::ads_wall_dir;
 AMREX_GPU_MANAGED GpuArray<amrex::Real, MAX_SPECIES> MFsurfchem::surfcov0;
 
 AMREX_GPU_MANAGED amrex::Real MFsurfchem::surf_site_num_dens;
@@ -16,6 +17,7 @@ AMREX_GPU_MANAGED amrex::Real MFsurfchem::k_beta;
 AMREX_GPU_MANAGED amrex::Real MFsurfchem::e_beta;
 
 AMREX_GPU_MANAGED int MFsurfchem::splitting_MFsurfchem;
+AMREX_GPU_MANAGED int MFsurfchem::conversion_MFsurfchem;
 
 void InitializeMFSurfchemNamespace()
 {
@@ -29,6 +31,9 @@ void InitializeMFSurfchemNamespace()
     // if n_ads_spec is set to 0 or not defined in the inputs file, quit the routine
     if (n_ads_spec==0) return;
 
+    ads_wall_dir = 2;
+    pp.query("ads_wall_dir",ads_wall_dir);
+
     // load default values to surfcov0 array
     for (int m=0;m<n_ads_spec;m++) surfcov0[m] = 0.;
 
@@ -75,6 +80,13 @@ void InitializeMFSurfchemNamespace()
     // get splitting type: first order (0), strang (1)
     splitting_MFsurfchem = 0; // default value
     pp.query("splitting_MFsurfchem",splitting_MFsurfchem);
+
+    // ads/des of species 1 -> ads of species 1 + desorption of species (1+n)
+    conversion_MFsurfchem = 0; // default value
+    pp.query("conversion_MFsurfchem",conversion_MFsurfchem);
+    if ( (n_ads_spec + conversion_MFsurfchem) > nspecies) {
+        Abort("ERROR: desorption species is not included in nspecies");
+    }
     return;
 }
 
@@ -96,7 +108,7 @@ void init_surfcov(MultiFab& surfcov, const amrex::Geometry& geom)
 
         amrex::ParallelForRNG(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k, RandomEngine const& engine) noexcept
         {
-            if (k==0) {
+            if ( (ads_wall_dir == 0 && i == 0) || (ads_wall_dir == 1 && j == 0) || (ads_wall_dir == 2 && k == 0) ) {
                 if (stoch_surfcov0==1) {
                     amrex::Real Ntot = rint(surf_site_num_dens*dx[0]*dx[1]);  // total number of reactive sites
                     GpuArray<int,MAX_SPECIES> Nocc;
@@ -126,7 +138,7 @@ void init_surfcov(MultiFab& surfcov, const amrex::Geometry& geom)
     return;
 }
 
-void sample_MFsurfchem(MultiFab& cu, MultiFab& prim, MultiFab& surfcov, MultiFab& dNadsdes,
+void sample_MFsurfchem(MultiFab& cu, MultiFab& prim, MultiFab& surfcov, MultiFab& dNadsdes, MultiFab& dNads, MultiFab& dNdes,
                        const amrex::Geometry& geom, const amrex::Real dt)
 {
 
@@ -139,12 +151,14 @@ void sample_MFsurfchem(MultiFab& cu, MultiFab& prim, MultiFab& surfcov, MultiFab
         const Array4<Real> & prim_arr = prim.array(mfi);
         const Array4<Real> & surfcov_arr = surfcov.array(mfi);
         const Array4<Real> & dNadsdes_arr = dNadsdes.array(mfi);
+        const Array4<Real> & dNads_arr = dNads.array(mfi);
+        const Array4<Real> & dNdes_arr = dNdes.array(mfi);
 
         amrex::Real Ntot = surf_site_num_dens*dx[0]*dx[1];  // total number of reactive sites
 
         amrex::ParallelForRNG(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k, RandomEngine const& engine) noexcept
         {
-            if (k==0) {
+            if ( (ads_wall_dir == 0 && i == 0) || (ads_wall_dir == 1 && j == 0) || (ads_wall_dir == 2 && k == 0) ) {
                 amrex::Real sumtheta = 0.;
                 for (int m=0;m<n_ads_spec;m++) {
                     sumtheta += surfcov_arr(i,j,k,m);
@@ -173,7 +187,13 @@ void sample_MFsurfchem(MultiFab& cu, MultiFab& prim, MultiFab& surfcov, MultiFab
                         Ndes = RandomPoisson(meanNdes,engine);
                     }
 
-                    dNadsdes_arr(i,j,k,m) = Nads-Ndes;
+		    if (conversion_MFsurfchem > 0) {
+		        dNads_arr(i,j,k,m) = Nads;
+			dNdes_arr(i,j,k,m) = Ndes;
+		    }
+		    else {
+		        dNadsdes_arr(i,j,k,m) = Nads-Ndes;
+		    }
                 }
             }
         });
@@ -182,7 +202,7 @@ void sample_MFsurfchem(MultiFab& cu, MultiFab& prim, MultiFab& surfcov, MultiFab
     return;
 }
 
-void update_MFsurfchem(MultiFab& cu, MultiFab& prim, MultiFab& surfcov, MultiFab& dNadsdes,
+void update_MFsurfchem(MultiFab& cu, MultiFab& prim, MultiFab& surfcov, MultiFab& dNadsdes, MultiFab& dNads, MultiFab& dNdes,
                        const amrex::Geometry& geom)
 {
     const GpuArray<Real, AMREX_SPACEDIM> dx = geom.CellSizeArray();
@@ -194,22 +214,40 @@ void update_MFsurfchem(MultiFab& cu, MultiFab& prim, MultiFab& surfcov, MultiFab
         const Array4<Real> & prim_arr = prim.array(mfi);
         const Array4<Real> & surfcov_arr = surfcov.array(mfi);
         const Array4<Real> & dNadsdes_arr = dNadsdes.array(mfi);
+        const Array4<Real> & dNads_arr = dNads.array(mfi);
+        const Array4<Real> & dNdes_arr = dNdes.array(mfi);
 
         amrex::Real Ntot = surf_site_num_dens*dx[0]*dx[1];  // total number of reactive sites
         amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
-            if (k==0) {
+            if ( (ads_wall_dir == 0 && i == 0) || (ads_wall_dir == 1 && j == 0) || (ads_wall_dir == 2 && k == 0) ) {
 
                 amrex::Real T_inst = prim_arr(i,j,k,4);
 
                 for (int m=0;m<n_ads_spec;m++) {
-                    amrex::Real dN = dNadsdes_arr(i,j,k,m);
-                    amrex::Real factor1 = molmass[m]/AVONUM/(dx[0]*dx[1]*dx[2]);
-		    amrex::Real factor2 = (e_beta*k_B*T_inst+(e0[m]+hcv[m]*T_inst)*molmass[m]/AVONUM)/(dx[0]*dx[1]*dx[2]);
-                    surfcov_arr(i,j,k,m) += dN/Ntot;
-                    cu_arr(i,j,k,0) -= factor1*dN;
-                    cu_arr(i,j,k,5+m) -= factor1*dN;
-                    cu_arr(i,j,k,4) -= factor2*dN;
+		    if (conversion_MFsurfchem > 0) {
+                        amrex::Real mconv = m + conversion_MFsurfchem;
+                        amrex::Real dNads = dNads_arr(i,j,k,m);
+                        amrex::Real dNdes = dNdes_arr(i,j,k,m);
+                        amrex::Real factor1 = molmass[m]/AVONUM/(dx[0]*dx[1]*dx[2]);
+                        amrex::Real factor1conv = molmass[mconv]/AVONUM/(dx[0]*dx[1]*dx[2]);
+                        amrex::Real factor2 = (e_beta*k_B*T_inst+(e0[m]+hcv[m]*T_inst)*molmass[m]/AVONUM)/(dx[0]*dx[1]*dx[2]);
+                        amrex::Real factor2conv = (e_beta*k_B*T_inst+(e0[mconv]+hcv[mconv]*T_inst)*molmass[mconv]/AVONUM)/(dx[0]*dx[1]*dx[2]);
+                        surfcov_arr(i,j,k,m) += (dNads-dNdes)/Ntot;
+                        cu_arr(i,j,k,0) -= factor1*dNads - factor1conv*dNdes;
+                        cu_arr(i,j,k,5+m) -= factor1*dNads;
+                        cu_arr(i,j,k,5+mconv) += factor1conv*dNdes;
+                        cu_arr(i,j,k,4) -= factor2*dNads - factor2conv*dNdes;
+		    }
+		    else {
+                        amrex::Real dN = dNadsdes_arr(i,j,k,m);
+                        amrex::Real factor1 = molmass[m]/AVONUM/(dx[0]*dx[1]*dx[2]);
+                        amrex::Real factor2 = (e_beta*k_B*T_inst+(e0[m]+hcv[m]*T_inst)*molmass[m]/AVONUM)/(dx[0]*dx[1]*dx[2]);
+                        surfcov_arr(i,j,k,m) += dN/Ntot;
+                        cu_arr(i,j,k,0) -= factor1*dN;
+                        cu_arr(i,j,k,5+m) -= factor1*dN;
+                        cu_arr(i,j,k,4) -= factor2*dN;
+		    }
                 }
             }
         });

src_MFsurfchem/MFsurfchem_namespace.H

@@ -1,6 +1,7 @@
 namespace MFsurfchem {
 
     extern AMREX_GPU_MANAGED int n_ads_spec;
+    extern AMREX_GPU_MANAGED int ads_wall_dir;
 
     extern AMREX_GPU_MANAGED GpuArray<amrex::Real, MAX_SPECIES> surfcov0;
     extern AMREX_GPU_MANAGED amrex::Real surf_site_num_dens;
@@ -15,5 +16,6 @@ namespace MFsurfchem {
     extern AMREX_GPU_MANAGED amrex::Real e_beta;
 
     extern AMREX_GPU_MANAGED int splitting_MFsurfchem;
+    extern AMREX_GPU_MANAGED int conversion_MFsurfchem;
 
 }

src_analysis/StructFact.H

@@ -37,6 +37,8 @@ class StructFact {
     // Define vector of unique selected variables
     amrex::Vector< int > var_u;
 
+    Geometry geom_sf;
+
 public:
 
     // Vector containing running sums of real and imaginary components
@@ -86,8 +88,7 @@ public:
                     amrex::MultiFab&,
                     bool unpack=true);
 
-    void WritePlotFile(const int, const amrex::Real, const amrex::Geometry&, 
-                       std::string, const int& zero_avg=1);
+    void WritePlotFile(const int, const amrex::Real, std::string, const int& zero_avg=1);
 
     void Finalize(amrex::MultiFab&, amrex::MultiFab&,
                   const int& zero_avg=1);

src_analysis/StructFact.cpp

@@ -203,6 +203,36 @@ void StructFact::define(const BoxArray& ba_in,
     cov_names[cnt] = x;
     cnt++;
   }
+
+
+  Box domain = ba_in.minimalBox();
+
+  Vector<Real> kspace_lo(AMREX_SPACEDIM);
+  Vector<Real> kspace_hi(AMREX_SPACEDIM);
+
+  for (int d=0; d<AMREX_SPACEDIM; ++d) {
+
+      if (domain.length(d) % 2 == 0) {
+          // even number of cells
+          kspace_lo[d] = -domain.length(d) / 2. - 0.5;
+          kspace_hi[d] =  domain.length(d) / 2. - 0.5;
+      } else {
+          // odd number of cells
+          kspace_lo[d] = -domain.length(d) / 2.;
+          kspace_hi[d] =  domain.length(d) / 2.;
+      }
+  }
+
+  RealBox kspace({AMREX_D_DECL(kspace_lo[0],kspace_lo[1],kspace_lo[2])},
+                 {AMREX_D_DECL(kspace_hi[0],kspace_hi[1],kspace_hi[2])});
+
+  // required only to define geom object
+  Vector<int> is_periodic(AMREX_SPACEDIM,1);
+
+  geom_sf.define(domain,&kspace,CoordSys::cartesian,is_periodic.data());
+
+
+
 }
 
 void StructFact::FortStructure(const MultiFab& variables,
@@ -505,7 +535,7 @@ void StructFact::ComputeFFT(const MultiFab& variables,
     }
 }
 
-void StructFact::WritePlotFile(const int step, const Real time, const Geometry& geom,
+void StructFact::WritePlotFile(const int step, const Real time,
                                std::string plotfile_base,
                                const int& zero_avg) {
 
@@ -546,7 +576,7 @@ void StructFact::WritePlotFile(const int step, const Real time, const Geometry&
   MultiFab::Copy(plotfile, cov_mag, 0, 0, NCOV, 0); // copy structure factor into plotfile
 
   // write a plotfile
-  WriteSingleLevelPlotfile(plotfilename1,plotfile,varNames,geom,time,step);
+  WriteSingleLevelPlotfile(plotfilename1,plotfile,varNames,geom_sf,time,step);
 
   //////////////////////////////////////////////////////////////////////////////////
   // Write out real and imaginary components of structure factor to plot file
@@ -579,7 +609,7 @@ void StructFact::WritePlotFile(const int step, const Real time, const Geometry&
   MultiFab::Copy(plotfile,cov_imag_temp,0,NCOV,NCOV,0);
 
   // write a plotfile
-  WriteSingleLevelPlotfile(plotfilename2,plotfile,varNames,geom,time,step);
+  WriteSingleLevelPlotfile(plotfilename2,plotfile,varNames,geom_sf,time,step);
 }
 
 void StructFact::Finalize(MultiFab& cov_real_in, MultiFab& cov_imag_in,