Merge remote-tracking branch 'origin' into feature/mblomquist/WriteSF…

…ToASCII

exec/compressible/inputs_giantfluct_3d

@@ -98,7 +98,7 @@
 
   # Xk and Yk at the wall for Dirichlet (concentrations) - set one to zero
   # Ordering: (species 1, x-dir), (species 2, x-dir), ... (species 1, y-dir), ...
-  bc_Yk_y_lo = 0.2 0.09316672 0.70683296     # lo BC
-  bc_Yk_y_hi = 0.3 0.40683328 0.29316704     # hi BC
+  bc_Yk_y_lo = 0.2 0.09316672 0.70683328     # lo BC
+  bc_Yk_y_hi = 0.3 0.40683328 0.29316672     # hi BC
 
 

exec/multispec/AdvanceTimestepBousq.cpp

@@ -273,7 +273,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
         // set normal velocity of physical domain boundaries
         MultiFabPhysBCDomainVel(umac[i],geom,i);
         // set transverse velocity behind physical boundaries
-        int is_inhomogeneous = 1;
+        int is_inhomogeneous = 0;
         MultiFabPhysBCMacVel(umac[i],geom,i,is_inhomogeneous);
         // fill periodic and interior ghost cells
         umac[i].FillBoundary(geom.periodicity());
@@ -707,7 +707,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
         // set normal velocity of physical domain boundaries
         MultiFabPhysBCDomainVel(umac[i],geom,i);
         // set transverse velocity behind physical boundaries
-        int is_inhomogeneous = 1;
+        int is_inhomogeneous = 0;
         MultiFabPhysBCMacVel(umac[i],geom,i,is_inhomogeneous);
         // fill periodic and interior ghost cells
         umac[i].FillBoundary(geom.periodicity());

exec/multispec/AdvanceTimestepInertial.cpp

@@ -309,7 +309,7 @@ void AdvanceTimestepInertial(std::array< MultiFab, AMREX_SPACEDIM >& umac,
         // set normal velocity of physical domain boundaries
         MultiFabPhysBCDomainVel(umac[i],geom,i);
         // set transverse velocity behind physical boundaries
-        int is_inhomogeneous = 1;
+        int is_inhomogeneous = 0;
         MultiFabPhysBCMacVel(umac[i],geom,i,is_inhomogeneous);
         // fill periodic and interior ghost cells
         umac[i].FillBoundary(geom.periodicity());
@@ -598,7 +598,7 @@ void AdvanceTimestepInertial(std::array< MultiFab, AMREX_SPACEDIM >& umac,
         // set normal velocity of physical domain boundaries
         MultiFabPhysBCDomainVel(umac[i],geom,i);
         // set transverse velocity behind physical boundaries
-        int is_inhomogeneous = 1;
+        int is_inhomogeneous = 0;
         MultiFabPhysBCMacVel(umac[i],geom,i,is_inhomogeneous);
         // fill periodic and interior ghost cells
         umac[i].FillBoundary(geom.periodicity());

src_common/ComputeAverages.cpp

@@ -95,8 +95,12 @@ void WriteHorizontalAverage(const MultiFab& mf_in, const int& dir, const int& in
 
 void WriteHorizontalAverageToMF(const MultiFab& mf_in, MultiFab& mf_out,
                                 const int& dir, const int& incomp,
-                                const int& ncomp)
+                                const int& ncomp, int outcomp)
 {
+    if (outcomp == -1) {
+        outcomp = incomp; // default outcomp is incomp unless specified
+    }
+
     // number of points in the averaging direction
     int npts = n_cells[dir];
 
@@ -166,7 +170,7 @@ void WriteHorizontalAverageToMF(const MultiFab& mf_in, MultiFab& mf_out,
         const Array4<Real> mf = mf_out.array(mfi);
 
         for (auto n=0; n<ncomp; ++n) {
-            comp = incomp+n;
+            comp = outcomp+n;
             for (auto k = lo.z; k <= hi.z; ++k) {
             for (auto j = lo.y; j <= hi.y; ++j) {
             for (auto i = lo.x; i <= hi.x; ++i) {
@@ -384,7 +388,44 @@ void ExtractSlice(const MultiFab& mf, MultiFab& mf_slice,
 
     // create a new DistributionMapping and define the MultiFab
     DistributionMapping dmap_slice(ba_slice);
-    mf_slice.define(ba_slice,dmap_slice,ncomp,0);
+    MultiFab mf_slice_tmp(ba_slice,dmap_slice,ncomp,0);
 
-    mf_slice.ParallelCopy(mf, incomp, 0, ncomp);
+    mf_slice_tmp.ParallelCopy(mf, incomp, 0, ncomp);
+
+    // now copy this into a multifab with index zero in the dir direction rather than slicepoint
+    // (structure factor code requires this)
+    dom_lo[dir] = 0;
+    dom_hi[dir] = 0;
+
+    Box domain_slice2(dom_lo,dom_hi);
+    BoxArray ba_slice2(domain_slice2);
+    ba_slice2.maxSize(IntVect(max_grid_slice));
+    mf_slice.define(ba_slice2,dmap_slice,ncomp,0);
+
+    for ( MFIter mfi(mf_slice_tmp,TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
+
+        const Box& bx = mfi.tilebox();
+
+        const Array4<Real> & slice = mf_slice.array(mfi);
+        const Array4<Real> & slice_tmp = mf_slice_tmp.array(mfi);
+
+        if (dir == 0) {
+            amrex::ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
+            {
+                slice(0,j,k,n) = slice_tmp(i,j,k,n);
+            });
+        }
+        if (dir == 1) {
+            amrex::ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
+            {
+                slice(i,0,k,n) = slice_tmp(i,j,k,n);
+            });
+        }
+        if (dir == 2) {
+            amrex::ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
+            {
+                slice(i,j,0,n) = slice_tmp(i,j,k,n);
+            });
+        }
+    }
 }

src_common/MultiFabPhysBC.cpp

@@ -390,7 +390,7 @@ void MultiFabPhysBCMacVel(MultiFab& vel, const Geometry& geom, int dim, int is_i
             }
         }
 
-        if ((dim != 0) && (bc_vel_lo[0] == 1 || bc_vel_hi[0] == 2) && (bx.bigEnd(0) > dom.bigEnd(0))) {
+        if ((dim != 0) && (bc_vel_hi[0] == 1 || bc_vel_hi[0] == 2) && (bx.bigEnd(0) > dom.bigEnd(0))) {
             if (bc_vel_hi[0] == 1) { // slip
                 amrex::ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {

src_common/common_functions.H

@@ -148,7 +148,7 @@ void WriteHorizontalAverage(const MultiFab& mf_in, const int& dir, const int& in
 
 void WriteHorizontalAverageToMF(const MultiFab& mf_in, MultiFab& mf_out,
                                 const int& dir, const int& incomp,
-                                const int& ncomp);
+                                const int& ncomp, int outcomp=-1);
 
 void ComputeVerticalAverage(const MultiFab & mf, MultiFab & mf_flat, const Geometry & geom,
                             const int& dir, const int& incomp, const int& ncomp,

src_common/common_functions.cpp

@@ -240,6 +240,7 @@ int                        common::plot_means;
 int                        common::plot_vars;
 int                        common::plot_covars;
 int                        common::plot_cross;
+int                        common::plot_deltaY_dir;
 int                        common::particle_motion;
 
 AMREX_GPU_MANAGED amrex::Real common::turb_a;
@@ -608,6 +609,7 @@ void InitializeCommonNamespace() {
     plot_vars = 0;
     plot_covars = 0;
     plot_cross = 0;
+    plot_deltaY_dir = -1;
     particle_motion = 0;
 
     // turblent forcing parameters
@@ -1135,6 +1137,7 @@ void InitializeCommonNamespace() {
     pp.query("plot_vars",plot_vars);
     pp.query("plot_covars",plot_covars);
     pp.query("plot_cross",plot_cross);
+    pp.query("plot_deltaY_dir",plot_deltaY_dir);
     pp.query("particle_motion",particle_motion);
     pp.query("turb_a",turb_a);
     pp.query("turb_b",turb_b);

src_common/common_namespace.H

@@ -314,6 +314,7 @@ namespace common {
     extern int                        plot_vars;
     extern int                        plot_covars;
     extern int                        plot_cross;
+    extern int                        plot_deltaY_dir;
     extern int                        particle_motion;
 
     // parameters for turbulent forcing example

src_compressible_stag/boundaryStag.cpp

@@ -32,6 +32,9 @@ void SetupCWallStag() {
        } else if (amrex::Math::abs(sumy-1) < 1.e-10) {
            GetMolfrac(bc_Yk_x_lo,bc_Xk_x_lo);
        }
+       else {
+           Abort("SetupCWallStag: lo-x; mass or mole fractions do not sum to 1");
+       }
     }
 
     if (bc_mass_lo[0] >= 3) {
@@ -69,6 +72,8 @@ void SetupCWallStag() {
           GetMassfrac(bc_Xk_x_hi,bc_Yk_x_hi);
        } else if (amrex::Math::abs(sumy-1) < 1.e-10) {
           GetMolfrac(bc_Yk_x_hi,bc_Xk_x_hi);
+       } else {
+           Abort("SetupCWallStag: hi-x; mass or mole fractions do not sum to 1");
        }
     }
 
@@ -108,6 +113,8 @@ void SetupCWallStag() {
           GetMassfrac(bc_Xk_y_lo,bc_Yk_y_lo);
        } else if (amrex::Math::abs(sumy-1) < 1.e-10) {
           GetMolfrac(bc_Yk_y_lo,bc_Xk_y_lo);
+       } else {
+           Abort("SetupCWallStag: lo-y; mass or mole fractions do not sum to 1");
        }
     }
 
@@ -146,6 +153,8 @@ void SetupCWallStag() {
           GetMassfrac(bc_Xk_y_hi,bc_Yk_y_hi);
        } else if (amrex::Math::abs(sumy-1) < 1.e-10) {
           GetMolfrac(bc_Yk_y_hi,bc_Xk_y_hi);
+       } else {
+           Abort("SetupCWallStag: hi-y; mass or mole fractions do not sum to 1");
        }
     }
 
@@ -185,6 +194,8 @@ void SetupCWallStag() {
           GetMassfrac(bc_Xk_z_lo,bc_Yk_z_lo);
        } else if (amrex::Math::abs(sumy-1) < 1.e-10) {
           GetMolfrac(bc_Yk_z_lo,bc_Xk_z_lo);
+       } else {
+           Abort("SetupCWallStag: lo-z; mass or mole fractions do not sum to 1");
        }
     }
 
@@ -223,6 +234,8 @@ void SetupCWallStag() {
           GetMassfrac(bc_Xk_z_hi,bc_Yk_z_hi);
        } else if (amrex::Math::abs(sumy-1) < 1.e-10) {
           GetMolfrac(bc_Yk_z_hi,bc_Xk_z_hi);
+       } else {
+           Abort("SetupCWallStag: hi-z; mass or mole fractions do not sum to 1");
        }
     }
 

src_compressible_stag/main_driver.cpp

@@ -208,8 +208,13 @@ void main_driver(const char* argv)
     if ((plot_cross) and ((cross_cell < 0) or (cross_cell > n_cells[0]-1))) {
         Abort("Cross cell needs to be within the domain: 0 <= cross_cell <= n_cells[0] - 1");
     }
-    if ((do_slab_sf) and ((membrane_cell <= 0) or (membrane_cell >= n_cells[0]-1))) {
-        Abort("Slab structure factor needs a membrane cell within the domain: 0 < cross_cell < n_cells[0] - 1");
+    if (project_dir >= 0) {
+        if (do_slab_sf and ((membrane_cell <= 0) or (membrane_cell >= n_cells[project_dir]-1))) {
+            Abort("Slab structure factor needs a membrane cell within the domain: 0 < cross_cell < n_cells[project_dir] - 1");
+        }
+        if (do_slab_sf and slicepoint >= 0) {
+            Abort("Cannot use do_slab_sf and slicepoint");
+        }
     }
     if ((project_dir >= 0) and ((do_1D) or (do_2D))) {
         Abort("Projected structure factors (project_dir) works only for 3D case");
@@ -276,9 +281,9 @@ void main_driver(const char* argv)
     MultiFab structFactPrimMF;
     MultiFab structFactConsMF;
 
-    // Structure factor for 2D averaged data
-    StructFact structFactPrimVerticalAverage;
-    StructFact structFactConsVerticalAverage;
+    // Structure factor for vertically-averaged or sliced data
+    StructFact structFactPrimFlattened;
+    StructFact structFactConsFlattened;
 
     // Structure factor for 2D averaged data (across a membrane)
     StructFact structFactPrimVerticalAverage0;
@@ -781,7 +786,11 @@ void main_driver(const char* argv)
 
             {
                 MultiFab X, XRot;
-                ComputeVerticalAverage(prim, X, geom, project_dir, 0, nprimvars);
+                if (slicepoint < 0) {
+                    ComputeVerticalAverage(prim, X, geom, project_dir, 0, nprimvars);
+                } else {
+                    ExtractSlice(prim, X, geom, project_dir, slicepoint, 0, 1);
+                }
                 XRot = RotateFlattenedMF(X);
                 ba_flat = XRot.boxArray();
                 dmap_flat = XRot.DistributionMap();
@@ -845,8 +854,8 @@ void main_driver(const char* argv)
             }
 
             if (do_slab_sf == 0) {
-                structFactPrimVerticalAverage.define(ba_flat,dmap_flat,prim_var_names,var_scaling_prim,2);
-                structFactConsVerticalAverage.define(ba_flat,dmap_flat,cons_var_names,var_scaling_cons,2);
+                structFactPrimFlattened.define(ba_flat,dmap_flat,prim_var_names,var_scaling_prim);
+                structFactConsFlattened.define(ba_flat,dmap_flat,cons_var_names,var_scaling_cons);
             }
             else {
                 structFactPrimVerticalAverage0.define(ba_flat,dmap_flat,prim_var_names,var_scaling_prim);
@@ -1468,19 +1477,27 @@ void main_driver(const char* argv)
                     {
                         MultiFab X, XRot;
 
-                        ComputeVerticalAverage(structFactPrimMF, X, geom, project_dir, 0, structVarsPrim);
+                        if (slicepoint < 0) {
+                            ComputeVerticalAverage(structFactPrimMF, X, geom, project_dir, 0, structVarsPrim);
+                        } else {
+                            ExtractSlice(structFactPrimMF, X, geom, project_dir, slicepoint, 0, structVarsPrim);
+                        }
                         XRot = RotateFlattenedMF(X);
                         master_project_rot_prim.ParallelCopy(XRot, 0, 0, structVarsPrim); 
-                        structFactPrimVerticalAverage.FortStructure(master_project_rot_prim,geom_flat);
+                        structFactPrimFlattened.FortStructure(master_project_rot_prim,geom_flat);
                     }
 
                     {
                         MultiFab X, XRot;
 
-                        ComputeVerticalAverage(structFactConsMF, X, geom, project_dir, 0, structVarsCons);
+                        if (slicepoint < 0) {
+                            ComputeVerticalAverage(structFactConsMF, X, geom, project_dir, 0, structVarsCons);
+                        } else {
+                            ExtractSlice(structFactConsMF, X, geom, project_dir, slicepoint, 0, structVarsCons);
+                        }
                         XRot = RotateFlattenedMF(X);
                         master_project_rot_cons.ParallelCopy(XRot, 0, 0, structVarsCons);
-                        structFactConsVerticalAverage.FortStructure(master_project_rot_cons,geom_flat);
+                        structFactConsFlattened.FortStructure(master_project_rot_cons,geom_flat);
                     }
 
                 }
@@ -1576,8 +1593,8 @@ void main_driver(const char* argv)
 
             if (project_dir >= 0) {
                 if (do_slab_sf == 0) {
-                    structFactPrimVerticalAverage.WritePlotFile(step,time,geom_flat,"plt_SF_prim_VerticalAverage");
-                    structFactConsVerticalAverage.WritePlotFile(step,time,geom_flat,"plt_SF_cons_VerticalAverage");
+                    structFactPrimFlattened.WritePlotFile(step,time,geom_flat,"plt_SF_prim_Flattened");
+                    structFactConsFlattened.WritePlotFile(step,time,geom_flat,"plt_SF_cons_Flattened");
                 }
                 else {
                     structFactPrimVerticalAverage0.WritePlotFile(step,time,geom_flat,"plt_SF_prim_VerticalAverageSlab0");

src_compressible_stag/writePlotFileStag.cpp

@@ -87,7 +87,11 @@ void WritePlotFileStag(int step,
 
 	if (nspec_surfcov>0) nplot += nspec_surfcov*6;
     }
-
+
+    if (plot_deltaY_dir != -1) {
+        nplot += nspecies;
+    }
+
     amrex::BoxArray ba = cuMeans.boxArray();
     amrex::DistributionMapping dmap = cuMeans.DistributionMap();
 
@@ -257,6 +261,16 @@ void WritePlotFileStag(int step,
 	}
     }
 
+    if (plot_deltaY_dir != -1) {
+        MultiFab Ybar(ba, dmap, nspecies, 0);
+        // Yk is component 6: in prim
+        WriteHorizontalAverageToMF(prim,Ybar,plot_deltaY_dir,6,nspecies,0);
+        Ybar.mult(-1.);
+        amrex::MultiFab::Add(Ybar,prim,6,0,nspecies,0);
+        amrex::MultiFab::Copy(plotfile,Ybar,0,cnt,nspecies,0);
+        cnt+= nspecies;
+    }
+
     // Set variable names
     cnt = 0;
 
@@ -445,6 +459,14 @@ void WritePlotFileStag(int step,
 
     }
 
+    if (plot_deltaY_dir != 1) {
+        x = "deltaYk_";
+        for (i=0; i<nspecies; i++) {
+            varNames[cnt] = x;
+            varNames[cnt++] += 48+i;
+        }
+    }
+
     AMREX_ASSERT(cnt==nplot);
 
     // write a plotfile