diff --git a/exec/compressible/inputs_giantfluct_3d b/exec/compressible/inputs_giantfluct_3d index 575927ee..9a18b240 100644 --- a/exec/compressible/inputs_giantfluct_3d +++ b/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 diff --git a/exec/multispec/AdvanceTimestepBousq.cpp b/exec/multispec/AdvanceTimestepBousq.cpp index 8d9f3931..bed33cd1 100644 --- a/exec/multispec/AdvanceTimestepBousq.cpp +++ b/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()); diff --git a/exec/multispec/AdvanceTimestepInertial.cpp b/exec/multispec/AdvanceTimestepInertial.cpp index bbfbc263..b89d1373 100644 --- a/exec/multispec/AdvanceTimestepInertial.cpp +++ b/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()); diff --git a/src_common/ComputeAverages.cpp b/src_common/ComputeAverages.cpp index 6880c70d..2b474cc5 100644 --- a/src_common/ComputeAverages.cpp +++ b/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 mf = mf_out.array(mfi); for (auto n=0; n & slice = mf_slice.array(mfi); + const Array4 & 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); + }); + } + } } diff --git a/src_common/MultiFabPhysBC.cpp b/src_common/MultiFabPhysBC.cpp index 529efaf5..06948555 100644 --- a/src_common/MultiFabPhysBC.cpp +++ b/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 { diff --git a/src_common/common_functions.H b/src_common/common_functions.H index 353d38ef..7cfeba66 100644 --- a/src_common/common_functions.H +++ b/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, diff --git a/src_common/common_functions.cpp b/src_common/common_functions.cpp index e9ca0ad0..bbc4d429 100644 --- a/src_common/common_functions.cpp +++ b/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); diff --git a/src_common/common_namespace.H b/src_common/common_namespace.H index 7e91aa22..b5493d3a 100644 --- a/src_common/common_namespace.H +++ 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 diff --git a/src_compressible_stag/boundaryStag.cpp b/src_compressible_stag/boundaryStag.cpp index b4c10039..96568857 100644 --- a/src_compressible_stag/boundaryStag.cpp +++ b/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"); } } diff --git a/src_compressible_stag/main_driver.cpp b/src_compressible_stag/main_driver.cpp index 98267a1e..da87869e 100644 --- a/src_compressible_stag/main_driver.cpp +++ b/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"); diff --git a/src_compressible_stag/writePlotFileStag.cpp b/src_compressible_stag/writePlotFileStag.cpp index c606b196..820c4cd8 100644 --- a/src_compressible_stag/writePlotFileStag.cpp +++ b/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