From 1d0fb7e398c8ee5fcb4720a2fc0db2cbc7ea8c6c Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Fri, 3 Jan 2025 19:54:34 -0800
Subject: [PATCH] more SF cleanup

---
 src_compressible_stag/main_driver.cpp | 428 ++++++++++++--------------
 1 file changed, 200 insertions(+), 228 deletions(-)

diff --git a/src_compressible_stag/main_driver.cpp b/src_compressible_stag/main_driver.cpp
index 38856d65..59fc217d 100644
--- a/src_compressible_stag/main_driver.cpp
+++ b/src_compressible_stag/main_driver.cpp
@@ -209,16 +209,22 @@ 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 (project_dir >= 0) {
+    if (struct_fact_int >0 and 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 < membrane_cell < n_cells[project_dir] - 1");
         }
+        if (do_1D) {
+            Abort("Projected structure factors (project_dir) does not work for do_1D case");
+        }
         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");
+        if (do_2D and slicepoint >= 0) {
+            Abort("Cannot use do_2D and slicepoint");
+        }
+        if (do_2D and project_dir != 2) {
+            Abort("Structure factors with do_2D requires project_dir == 2");
+        }
     }
     if ((all_correl > 1) or (all_correl < 0)) {
         Abort("all_correl can be 0 or 1");
@@ -285,26 +291,41 @@ void main_driver(const char* argv)
     // Standard 3D structure factors
     StructFact structFactPrim;
     StructFact structFactCons;
+    // MultiFabs to copy data into for snapshots for full 3D data
     MultiFab structFactPrimMF;
     MultiFab structFactConsMF;
 
     // Structure factor for vertically-averaged or sliced data
+    // these are enabled if project_dir >= 0 AND do_slab_sf == 0
     StructFact structFactPrimFlattened;
     StructFact structFactConsFlattened;
-
     // Structure factor for 2D averaged data (across a membrane)
+    // these are enabled if project_dir >= 0 AND do_slab_sf == 1
     StructFact structFactPrimVerticalAverageMembraneLo;
     StructFact structFactPrimVerticalAverageMembraneHi;
     StructFact structFactConsVerticalAverageMembraneLo;
     StructFact structFactConsVerticalAverageMembraneHi;
-    MultiFab primFlattenedMaster;
-    MultiFab consFlattenedMaster;
-
-    // Vector of structure factors for 2D simulation
+    // Vector of structure factors of slices for 2D simulation
+    // these are enabled if do_2D (this mode assumes z slices; project_dir must equal 2)
     Vector < StructFact > structFactPrimArray;
     Vector < StructFact > structFactConsArray;
-    MultiFab prim2DFlattenedMaster;
-    MultiFab cons2DFlattenedMaster;
+    // MultiFabs to copy data into for snapshots for flattened
+    MultiFab structFactPrimFlattenedMF;
+    MultiFab structFactConsFlattenedMF;
+    
+    Geometry geom_flat;
+    BoxArray ba_flat;
+    DistributionMapping dmap_flat;
+
+    // Structure factor for surface coverage slice
+    // these are enabled if n_ads_spec > 0 and assumes the k=0 plane is the slice of interest
+    StructFact structFactSurfCov;
+    // MultiFab to copy data into for snapshop
+    MultiFab structFactSurfCovMF;
+    
+    Geometry geom_surfcov;
+    BoxArray ba_surfcov;
+    DistributionMapping dmap_surfcov;
 
 #if defined(TURB)
     // Structure factor for compressible turbulence
@@ -313,20 +334,6 @@ void main_driver(const char* argv)
     StructFact turbStructFactScalar; // scalars 
 #endif
 
-    // surface coverage structure factor
-    StructFact surfcovStructFact;
-    MultiFab surfcovFlattenedMaster;
-    Geometry surfcov_geom_flat;
-    BoxArray surfcov_ba_flat;
-    DistributionMapping surfcov_dmap_flat;
-    
-    Geometry geom_flat;
-    Geometry geom_flat_2D;
-    BoxArray ba_flat;
-    BoxArray ba_flat_2D;
-    DistributionMapping dmap_flat;
-    DistributionMapping dmap_flat_2D;
-
     // "primitive" variable structure factor will contain
     // rho
     // vel (shifted)
@@ -758,165 +765,131 @@ void main_driver(const char* argv)
     ///////////////////////////////////////////
 
     if (struct_fact_int > 0) {
-        structFactPrimMF.define(ba, dmap, structVarsPrim, 0);
-        structFactPrim.define(ba,dmap,prim_var_names,var_scaling_prim);
-            
-        structFactConsMF.define(ba, dmap, structVarsCons, 0);
-        structFactCons.define(ba,dmap,cons_var_names,var_scaling_cons);
-        
+
+        if ((do_1D==0) and (do_2D==0)) {
+            structFactPrim.define(ba,dmap,prim_var_names,var_scaling_prim);
+            structFactCons.define(ba,dmap,cons_var_names,var_scaling_cons);
+            structFactConsMF.define(ba,dmap,structVarsCons,0);
+            structFactPrimMF.define(ba,dmap,structVarsPrim,0);
+        }
+
         // structure factor class for vertically-averaged dataset
         if (project_dir >= 0) {
+
             MultiFab Flattened; // flattened multifab define below
 
-            // we are only calling ComputeVerticalAverage or ExtractSlice here to obtain
+            // we are only calling ExtractSlice here to obtain
             // a built version of Flattened so can obtain what we need to build the
             // structure factor and geometry objects for flattened data
-            if (slicepoint < 0) {
-                ComputeVerticalAverage(prim, Flattened, geom, project_dir, 0, nprimvars);
-            } else {
-                ExtractSlice(prim, Flattened, geom, project_dir, slicepoint, 0, 1);
-            }
+            ExtractSlice(prim, Flattened, geom, project_dir, 0, 0, 1);
+
             ba_flat = Flattened.boxArray();
             dmap_flat = Flattened.DistributionMap();
-            primFlattenedMaster.define(ba_flat,dmap_flat,structVarsPrim,0);
-            consFlattenedMaster.define(ba_flat,dmap_flat,structVarsCons,0);
-            {
-                Box domain_flat = primFlattenedMaster.boxArray().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);
+            structFactPrimFlattenedMF.define(ba_flat,dmap_flat,structVarsPrim,0);
+            structFactConsFlattenedMF.define(ba_flat,dmap_flat,structVarsCons,0);
 
-                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];
+            Box domain_flat = structFactPrimFlattenedMF.boxArray().minimalBox();
 
-                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());
-            }
-           
+            // 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);
 
-            if (do_slab_sf == 0) {
-                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 {
-                structFactPrimVerticalAverageMembraneLo.define(ba_flat,dmap_flat,prim_var_names,var_scaling_prim);
-                structFactPrimVerticalAverageMembraneHi.define(ba_flat,dmap_flat,prim_var_names,var_scaling_prim);
-                structFactConsVerticalAverageMembraneLo.define(ba_flat,dmap_flat,cons_var_names,var_scaling_cons);
-                structFactConsVerticalAverageMembraneHi.define(ba_flat,dmap_flat,cons_var_names,var_scaling_cons);
+            for (int d=0; d<AMREX_SPACEDIM; ++d) {
+                projected_lo[d] = prob_lo[d];
+                projected_hi[d] = prob_hi[d];
             }
-    
-        }
-
-        if (do_2D) { // 2D is coded only for XY plane
-
-            {
-                MultiFab X;
-                ExtractSlice(prim, X, geom, 2, 0, 0, nprimvars);
-                ba_flat_2D = X.boxArray();
-                dmap_flat_2D = X.DistributionMap();
-                prim2DFlattenedMaster.define(ba_flat_2D,dmap_flat_2D,structVarsPrim,0);
-                cons2DFlattenedMaster.define(ba_flat_2D,dmap_flat_2D,structVarsCons,0);
-
-                IntVect dom_lo_flat(AMREX_D_DECL(0,0,0));
-                IntVect dom_hi_flat;
-                dom_hi_flat[0] = n_cells[0]-1;
-                dom_hi_flat[1] = n_cells[1]-1;
-                dom_hi_flat[2] = 0;
-                Box domain_flat(dom_lo_flat, dom_hi_flat);
+            projected_lo[project_dir] = 0.;
+            projected_hi[project_dir] = dx[project_dir];
 
-                // This defines the physical box
-                Vector<Real> projected_hi(AMREX_SPACEDIM);
-                for (int d=0; d<AMREX_SPACEDIM; d++) {
-                    projected_hi[d] = prob_hi[d];
-                }
-                projected_hi[AMREX_SPACEDIM-1] = prob_hi[2] / n_cells[2];
-
-                RealBox real_box_flat({AMREX_D_DECL(     prob_lo[0],     prob_lo[1],     prob_lo[2])},
-                                      {AMREX_D_DECL(projected_hi[0],projected_hi[1],projected_hi[2])});
+            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_2D.define(domain_flat,&real_box_flat,CoordSys::cartesian,is_periodic.data());
+            // This defines a Geometry object
+            geom_flat.define(domain_flat,&real_box_flat,CoordSys::cartesian,is_periodic.data());
 
-            }
+            if (do_2D) {
 
-            structFactPrimArray.resize(n_cells[2]);
-            structFactConsArray.resize(n_cells[2]);
+                structFactPrimArray.resize(n_cells[project_dir]);
+                structFactConsArray.resize(n_cells[project_dir]);
 
-            for (int i = 0; i < n_cells[2]; ++i) { 
-                structFactPrimArray[i].define(ba_flat_2D,dmap_flat_2D,prim_var_names,var_scaling_prim,2);
-                structFactConsArray[i].define(ba_flat_2D,dmap_flat_2D,cons_var_names,var_scaling_cons,2);
-            }
+                for (int i = 0; i < n_cells[project_dir]; ++i) { 
+                    structFactPrimArray[i].define(ba_flat,dmap_flat,prim_var_names,var_scaling_prim);
+                    structFactConsArray[i].define(ba_flat,dmap_flat,cons_var_names,var_scaling_cons);
+                }
 
+            } else {
+                if (do_slab_sf == 0) {
+                    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 {
+                    structFactPrimVerticalAverageMembraneLo.define(ba_flat,dmap_flat,prim_var_names,var_scaling_prim);
+                    structFactPrimVerticalAverageMembraneHi.define(ba_flat,dmap_flat,prim_var_names,var_scaling_prim);
+                    structFactConsVerticalAverageMembraneLo.define(ba_flat,dmap_flat,cons_var_names,var_scaling_cons);
+                    structFactConsVerticalAverageMembraneHi.define(ba_flat,dmap_flat,cons_var_names,var_scaling_cons);
+                }
+            }
+    
         }
-    }
 
-    if (n_ads_spec>0) {
+        if (n_ads_spec>0) {
 
-        MultiFab Flattened;  // flattened multifab defined below
-
-        // we are only calling ExtractSlice here to obtain
-        // a built version of Flattened so can obtain what we need to build the
-        // structure factor and geometry objects for flattened data
-        // assume surface covered is stored in the "k" direction in the k=0 coordinate.
-        int surfcov_dir = 2;
-        int surfcov_plane = 0;
-        int surfcov_structVars = n_ads_spec;
-        int surfcov_nPairs = surfcov_structVars*(surfcov_structVars+1)/2;
-
-        Vector< std::string > surfcov_var_names;
-        surfcov_var_names.resize(surfcov_structVars);
-        for (int d=0; d<surfcov_structVars; d++) {
-            x = "surfCov";
-            x += (48+d);
-            surfcov_var_names[d] = x;
-        }
+            MultiFab Flattened;  // flattened multifab defined below
 
-        Vector<Real> surfcov_var_scaling(surfcov_nPairs);
-        for (int d=0; d<surfcov_var_scaling.size(); ++d) {
-            surfcov_var_scaling[d] = 1.;
-        }
+            // we are only calling ExtractSlice here to obtain
+            // a built version of Flattened so can obtain what we need to build the
+            // structure factor and geometry objects for flattened data
+            // assume surface covered is stored in the "k" direction in the k=0 coordinate.
+            int surfcov_dir = 2;
+            int surfcov_plane = 0;
+            int surfcov_structVars = n_ads_spec;
+            int surfcov_nPairs = surfcov_structVars*(surfcov_structVars+1)/2;
+
+            Vector< std::string > surfcov_var_names;
+            surfcov_var_names.resize(surfcov_structVars);
+            for (int d=0; d<surfcov_structVars; d++) {
+                x = "surfCov";
+                x += (48+d);
+                surfcov_var_names[d] = x;
+            }
+
+            Vector<Real> surfcov_var_scaling(surfcov_nPairs);
+            for (int d=0; d<surfcov_var_scaling.size(); ++d) {
+                surfcov_var_scaling[d] = 1.;
+            }
       
-        ExtractSlice(surfcov, Flattened, geom, surfcov_dir, surfcov_plane, 0, surfcov_structVars);
-        BoxArray surfcov_ba_flat = Flattened.boxArray();
-        const DistributionMapping& dmap_flat = Flattened.DistributionMap();
-        surfcovFlattenedMaster.define(surfcov_ba_flat,dmap_flat,surfcov_structVars,0);
-        {
-            Box domain_flat = surfcovFlattenedMaster.boxArray().minimalBox();
+            ExtractSlice(surfcov, Flattened, geom, surfcov_dir, surfcov_plane, 0, surfcov_structVars);
+            BoxArray ba_surfcov = Flattened.boxArray();
+            const DistributionMapping& dmap_surfcov = Flattened.DistributionMap();
+            structFactSurfCovMF.define(ba_surfcov,dmap_surfcov,surfcov_structVars,0);
+            {
+                Box domain_surfcov = structFactSurfCovMF.boxArray().minimalBox();
         
-            // This defines the physical box
-            // we retain prob_lo and prob_hi in all directions except surfcov_dir,
-            // where the physical size is 0 to dx[surfcov_dir]
-            Vector<Real> projected_lo(AMREX_SPACEDIM);
-            Vector<Real> projected_hi(AMREX_SPACEDIM);
+                // This defines the physical box
+                // we retain prob_lo and prob_hi in all directions except surfcov_dir,
+                // where the physical size is 0 to dx[surfcov_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[surfcov_dir] = 0.;
-            projected_hi[surfcov_dir] = (prob_hi[surfcov_dir] - prob_lo[surfcov_dir]) / n_cells[surfcov_dir];
+                for (int d=0; d<AMREX_SPACEDIM; ++d) {
+                    projected_lo[d] = prob_lo[d];
+                    projected_hi[d] = prob_hi[d];
+                }
+                projected_lo[surfcov_dir] = 0.;
+                projected_hi[surfcov_dir] = dx[surfcov_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());
+                RealBox real_box_surfcov({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
-            surfcov_geom_flat.define(domain_flat,&real_box_flat,CoordSys::cartesian,is_periodic.data());
-        }
+                // This defines a Geometry object
+                geom_surfcov.define(domain_surfcov,&real_box_surfcov,CoordSys::cartesian,is_periodic.data());
+            }
 
-        surfcovStructFact.define(surfcov_ba_flat,dmap_flat,surfcov_var_names,surfcov_var_scaling);
+            structFactSurfCov.define(ba_surfcov,dmap_surfcov,surfcov_var_names,surfcov_var_scaling);
+        }
     }
 
     /////////////////////////////////////////////////
@@ -1361,89 +1334,88 @@ void main_driver(const char* argv)
 
             if (project_dir >= 0) {
 
-                if (do_slab_sf == 0) {
-                    
-                    {
-                        MultiFab X;
+                if (do_2D) {
 
-                        if (slicepoint < 0) {
-                            ComputeVerticalAverage(structFactPrimMF, X, geom, project_dir, 0, structVarsPrim);
-                        } else {
-                            ExtractSlice(structFactPrimMF, X, geom, project_dir, slicepoint, 0, structVarsPrim);
-                        }
-                        primFlattenedMaster.ParallelCopy(X, 0, 0, structVarsPrim); 
-                        structFactPrimFlattened.FortStructure(primFlattenedMaster);
-                    }
+                    for (int i=0; i<n_cells[project_dir]; ++i) {
 
-                    {
-                        MultiFab X;
+                        {
+                            MultiFab Flattened;
 
-                        if (slicepoint < 0) {
-                            ComputeVerticalAverage(structFactConsMF, X, geom, project_dir, 0, structVarsCons);
-                        } else {
-                            ExtractSlice(structFactConsMF, X, geom, project_dir, slicepoint, 0, structVarsCons);
+                            ExtractSlice(structFactPrimMF, Flattened, geom, project_dir, i, 0, structVarsPrim);
+                            structFactPrimFlattenedMF.ParallelCopy(Flattened, 0, 0, structVarsPrim); 
+                            structFactPrimArray[i].FortStructure(structFactPrimFlattenedMF);
                         }
-                        consFlattenedMaster.ParallelCopy(X, 0, 0, structVarsCons);
-                        structFactConsFlattened.FortStructure(consFlattenedMaster);
-                    }
-
-                }
-                else {
-                    
-                    {
-                        MultiFab X;
 
-                        ComputeVerticalAverage(structFactPrimMF, X, geom, project_dir, 0, structVarsPrim, 0, membrane_cell-1);
-                        primFlattenedMaster.ParallelCopy(X, 0, 0, structVarsPrim);
-                        structFactPrimVerticalAverageMembraneLo.FortStructure(primFlattenedMaster);
-                    }
+                        {
+                            MultiFab Flattened;
 
-                    {
-                        MultiFab X;
+                            ExtractSlice(structFactConsMF, Flattened, geom, project_dir, i, 0, structVarsCons);
+                            structFactConsFlattenedMF.ParallelCopy(Flattened, 0, 0, structVarsCons); 
+                            structFactConsArray[i].FortStructure(structFactConsFlattenedMF);
+                        }
 
-                        ComputeVerticalAverage(structFactPrimMF, X, geom, project_dir, 0, structVarsPrim, membrane_cell, n_cells[project_dir]-1);
-                        primFlattenedMaster.ParallelCopy(X, 0, 0, structVarsPrim); 
-                        structFactPrimVerticalAverageMembraneHi.FortStructure(primFlattenedMaster);
                     }
+                } else {
 
-                    {
-                        MultiFab X;
+                    if (do_slab_sf == 0) {
+                    
+                        {
+                            MultiFab Flattened;
+
+                            if (slicepoint < 0) {
+                                ComputeVerticalAverage(structFactPrimMF, Flattened, geom, project_dir, 0, structVarsPrim);
+                            } else {
+                                ExtractSlice(structFactPrimMF, Flattened, geom, project_dir, slicepoint, 0, structVarsPrim);
+                            }
+                            structFactPrimFlattenedMF.ParallelCopy(Flattened, 0, 0, structVarsPrim); 
+                            structFactPrimFlattened.FortStructure(structFactPrimFlattenedMF);
+                        }
 
-                        ComputeVerticalAverage(structFactConsMF, X, geom, project_dir, 0, structVarsCons, 0, membrane_cell-1);
-                        consFlattenedMaster.ParallelCopy(X, 0, 0, structVarsCons); 
-                        structFactConsVerticalAverageMembraneLo.FortStructure(consFlattenedMaster);
-                    }
+                        {
+                            MultiFab Flattened;
 
-                    {
-                        MultiFab X;
+                            if (slicepoint < 0) {
+                                ComputeVerticalAverage(structFactConsMF, Flattened, geom, project_dir, 0, structVarsCons);
+                            } else {
+                                ExtractSlice(structFactConsMF, Flattened, geom, project_dir, slicepoint, 0, structVarsCons);
+                            }
+                            structFactConsFlattenedMF.ParallelCopy(Flattened, 0, 0, structVarsCons);
+                            structFactConsFlattened.FortStructure(structFactConsFlattenedMF);
+                        }
+                    } else {
+                    
+                        {
+                            MultiFab Flattened;
 
-                        ComputeVerticalAverage(structFactConsMF, X, geom, project_dir, 0, structVarsCons, membrane_cell, n_cells[project_dir]-1);
-                        consFlattenedMaster.ParallelCopy(X, 0, 0, structVarsCons); 
-                        structFactConsVerticalAverageMembraneHi.FortStructure(consFlattenedMaster);
-                    }
-                }
-            }
+                            ComputeVerticalAverage(structFactPrimMF, Flattened, geom, project_dir, 0, structVarsPrim, 0, membrane_cell-1);
+                            structFactPrimFlattenedMF.ParallelCopy(Flattened, 0, 0, structVarsPrim);
+                            structFactPrimVerticalAverageMembraneLo.FortStructure(structFactPrimFlattenedMF);
+                        }
 
-            if (do_2D) {
+                        {
+                            MultiFab Flattened;
 
-                for (int i=0; i<n_cells[2]; ++i) {
+                            ComputeVerticalAverage(structFactPrimMF, Flattened, geom, project_dir, 0, structVarsPrim, membrane_cell, n_cells[project_dir]-1);
+                            structFactPrimFlattenedMF.ParallelCopy(Flattened, 0, 0, structVarsPrim); 
+                            structFactPrimVerticalAverageMembraneHi.FortStructure(structFactPrimFlattenedMF);
+                        }
 
-                    {
-                        MultiFab X;
+                        {
+                            MultiFab Flattened;
 
-                        ExtractSlice(structFactPrimMF, X, geom, 2, i, 0, structVarsPrim);
-                        prim2DFlattenedMaster.ParallelCopy(X, 0, 0, structVarsPrim); 
-                        structFactPrimArray[i].FortStructure(prim2DFlattenedMaster);
-                    }
+                            ComputeVerticalAverage(structFactConsMF, Flattened, geom, project_dir, 0, structVarsCons, 0, membrane_cell-1);
+                            structFactConsFlattenedMF.ParallelCopy(Flattened, 0, 0, structVarsCons); 
+                            structFactConsVerticalAverageMembraneLo.FortStructure(structFactConsFlattenedMF);
+                        }
 
-                    {
-                        MultiFab X;
+                        {
+                            MultiFab Flattened;
 
-                        ExtractSlice(structFactConsMF, X, geom, 2, i, 0, structVarsCons);
-                        cons2DFlattenedMaster.ParallelCopy(X, 0, 0, structVarsCons); 
-                        structFactConsArray[i].FortStructure(cons2DFlattenedMaster);
+                            ComputeVerticalAverage(structFactConsMF, Flattened, geom, project_dir, 0, structVarsCons, membrane_cell, n_cells[project_dir]-1);
+                            structFactConsFlattenedMF.ParallelCopy(Flattened, 0, 0, structVarsCons); 
+                            structFactConsVerticalAverageMembraneHi.FortStructure(structFactConsFlattenedMF);
+                        }
                     }
-
                 }
             }
 
@@ -1453,8 +1425,8 @@ void main_driver(const char* argv)
                 int surfcov_structVars = n_ads_spec;
                 MultiFab Flattened;  // flattened multifab defined below
                 ExtractSlice(surfcov, Flattened, geom, surfcov_dir, surfcov_plane, 0, surfcov_structVars);
-                surfcovFlattenedMaster.ParallelCopy(Flattened,0,0,surfcov_structVars);
-                surfcovStructFact.FortStructure(surfcovFlattenedMaster);
+                structFactSurfCovMF.ParallelCopy(Flattened,0,0,surfcov_structVars);
+                structFactSurfCov.FortStructure(structFactSurfCovMF);
             }
 
         }
@@ -1486,10 +1458,10 @@ void main_driver(const char* argv)
                     
                 MultiFab prim_mag, prim_realimag, cons_mag, cons_realimag;
 
-                prim_mag.define(ba_flat_2D,dmap_flat_2D,structFactPrimArray[0].get_ncov(),0);
-                prim_realimag.define(ba_flat_2D,dmap_flat_2D,2*structFactPrimArray[0].get_ncov(),0);
-                cons_mag.define(ba_flat_2D,dmap_flat_2D,structFactConsArray[0].get_ncov(),0);
-                cons_realimag.define(ba_flat_2D,dmap_flat_2D,2*structFactConsArray[0].get_ncov(),0);
+                prim_mag.define(ba_flat,dmap_flat,structFactPrimArray[0].get_ncov(),0);
+                prim_realimag.define(ba_flat,dmap_flat,2*structFactPrimArray[0].get_ncov(),0);
+                cons_mag.define(ba_flat,dmap_flat,structFactConsArray[0].get_ncov(),0);
+                cons_realimag.define(ba_flat,dmap_flat,2*structFactConsArray[0].get_ncov(),0);
 
                 prim_mag.setVal(0.0);
                 cons_mag.setVal(0.0);
@@ -1507,15 +1479,15 @@ void main_driver(const char* argv)
                 prim_realimag.mult(ncellsinv);
                 cons_realimag.mult(ncellsinv);
 
-                WritePlotFilesSF_2D(prim_mag,prim_realimag,geom_flat_2D,step,time,
+                WritePlotFilesSF_2D(prim_mag,prim_realimag,geom_flat,step,time,
                                     structFactPrimArray[0].get_names(),"plt_SF_prim_2D");
-                WritePlotFilesSF_2D(cons_mag,cons_realimag,geom_flat_2D,step,time,
+                WritePlotFilesSF_2D(cons_mag,cons_realimag,geom_flat,step,time,
                                     structFactConsArray[0].get_names(),"plt_SF_cons_2D");
 
             }
 
             if (n_ads_spec > 0) {
-                surfcovStructFact.WritePlotFile(step,time,surfcov_geom_flat,"plt_SF_surfcov");
+                structFactSurfCov.WritePlotFile(step,time,geom_surfcov,"plt_SF_surfcov");
             }
         }