From 13f9d43e28564bd9a2efa4caa3c6d88229562775 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Tue, 5 Nov 2024 09:02:30 -0800 Subject: [PATCH 1/8] Implement slice structure factor for compressible_stag To use, set project_dir (normal direction of plane) and slicepoint (coordinate in project_dir to slice) in the inputs file; e.g. if you care about the k=0 plane use project_dir=2 and slicepoint=0 If you do not specify slicepoint it will take the vertical average over project_dir over the entire dataset plt_SF_prim_Flattened plt_SF_cons_Flattened contains the result --- src_compressible_stag/main_driver.cpp | 45 ++++++++++++++++++--------- 1 file changed, 31 insertions(+), 14 deletions(-) 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"); From c01b5831d4d58e09eb721cf453cc52df9acad555 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Mon, 18 Nov 2024 11:31:25 -0800 Subject: [PATCH 2/8] fix bug when hi-x wall moves in the y-direction --- src_common/MultiFabPhysBC.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 { From 6041c07a165cd28ea992c436347af05a286e4ad6 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Mon, 18 Nov 2024 11:47:47 -0800 Subject: [PATCH 3/8] fix moving wall solver in delta form --- exec/multispec/AdvanceTimestepBousq.cpp | 4 ++-- exec/multispec/AdvanceTimestepInertial.cpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) 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()); From 0cac562fdd31693f5110df4720fa2000566f10ea Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Tue, 19 Nov 2024 19:21:28 -0800 Subject: [PATCH 4/8] introduce code that forces the slice multifab to have index 0:0 in the dir direction, regardless of slicepoint the structure factor flattened code requires this --- src_common/ComputeAverages.cpp | 41 ++++++++++++++++++++++++++++++++-- 1 file changed, 39 insertions(+), 2 deletions(-) diff --git a/src_common/ComputeAverages.cpp b/src_common/ComputeAverages.cpp index 6880c70d..e63b371d 100644 --- a/src_common/ComputeAverages.cpp +++ b/src_common/ComputeAverages.cpp @@ -384,7 +384,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 & 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) = slice_tmp(i,j,k); + }); + } + if (dir == 1) { + amrex::ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + slice(i,0,k) = slice_tmp(i,j,k); + }); + } + if (dir == 2) { + amrex::ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + slice(i,j,0) = slice_tmp(i,j,k); + }); + } + } } From 81ae8acc4f0aaf00ffac22f3b21d2fdd842e42d3 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Thu, 21 Nov 2024 14:32:39 -0800 Subject: [PATCH 5/8] fix copy of slice MF components --- src_common/ComputeAverages.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src_common/ComputeAverages.cpp b/src_common/ComputeAverages.cpp index e63b371d..34791dcc 100644 --- a/src_common/ComputeAverages.cpp +++ b/src_common/ComputeAverages.cpp @@ -408,19 +408,19 @@ void ExtractSlice(const MultiFab& mf, MultiFab& mf_slice, if (dir == 0) { amrex::ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - slice(0,j,k) = slice_tmp(i,j,k); + 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) = slice_tmp(i,j,k); + 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) = slice_tmp(i,j,k); + slice(i,j,0,n) = slice_tmp(i,j,k,n); }); } } From dd17133f446417d48966050abca99d0dad658e5a Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Wed, 4 Dec 2024 15:29:27 -0800 Subject: [PATCH 6/8] add ability to write deltaY to plotfile must specify plot_deltaY = 1 project_dir = (0, 1, or 2) --- src_common/ComputeAverages.cpp | 8 +++++-- src_common/common_functions.H | 2 +- src_common/common_functions.cpp | 3 +++ src_common/common_namespace.H | 1 + src_compressible_stag/writePlotFileStag.cpp | 24 ++++++++++++++++++++- 5 files changed, 34 insertions(+), 4 deletions(-) diff --git a/src_common/ComputeAverages.cpp b/src_common/ComputeAverages.cpp index 34791dcc..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; n0) nplot += nspec_surfcov*6; } - + + if (plot_deltaY == 1) { + nplot += nspecies; + } + amrex::BoxArray ba = cuMeans.boxArray(); amrex::DistributionMapping dmap = cuMeans.DistributionMap(); @@ -257,6 +261,16 @@ void WritePlotFileStag(int step, } } + if (plot_deltaY == 1) { + MultiFab Ybar(ba, dmap, nspecies, 0); + // Yk is component 6: in prim + WriteHorizontalAverageToMF(prim,Ybar,project_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 == 1) { + x = "deltaYk_"; + for (i=0; i Date: Wed, 4 Dec 2024 17:21:52 -0800 Subject: [PATCH 7/8] fix inputs_giantflux_3d to work by ensuring bc YK's sum to 1 add error checking for this condition in the compressible_stag code --- exec/compressible/inputs_giantfluct_3d | 4 ++-- src_compressible_stag/boundaryStag.cpp | 13 +++++++++++++ 2 files changed, 15 insertions(+), 2 deletions(-) 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/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"); } } From 9f1ccd2cbfb4ed0f978a245541b33fbdf4d0c9ce Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Wed, 4 Dec 2024 18:05:26 -0800 Subject: [PATCH 8/8] use plot_deltaY_dir instead of project_dir --- src_common/common_functions.cpp | 6 +++--- src_common/common_namespace.H | 2 +- src_compressible_stag/writePlotFileStag.cpp | 8 ++++---- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src_common/common_functions.cpp b/src_common/common_functions.cpp index 8dfa7029..bbc4d429 100644 --- a/src_common/common_functions.cpp +++ b/src_common/common_functions.cpp @@ -240,7 +240,7 @@ int common::plot_means; int common::plot_vars; int common::plot_covars; int common::plot_cross; -int common::plot_deltaY; +int common::plot_deltaY_dir; int common::particle_motion; AMREX_GPU_MANAGED amrex::Real common::turb_a; @@ -609,7 +609,7 @@ void InitializeCommonNamespace() { plot_vars = 0; plot_covars = 0; plot_cross = 0; - plot_deltaY = 0; + plot_deltaY_dir = -1; particle_motion = 0; // turblent forcing parameters @@ -1137,7 +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",plot_deltaY); + 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 bddc3961..b5493d3a 100644 --- a/src_common/common_namespace.H +++ b/src_common/common_namespace.H @@ -314,7 +314,7 @@ namespace common { extern int plot_vars; extern int plot_covars; extern int plot_cross; - extern int plot_deltaY; + extern int plot_deltaY_dir; extern int particle_motion; // parameters for turbulent forcing example diff --git a/src_compressible_stag/writePlotFileStag.cpp b/src_compressible_stag/writePlotFileStag.cpp index 94f56adc..820c4cd8 100644 --- a/src_compressible_stag/writePlotFileStag.cpp +++ b/src_compressible_stag/writePlotFileStag.cpp @@ -88,7 +88,7 @@ void WritePlotFileStag(int step, if (nspec_surfcov>0) nplot += nspec_surfcov*6; } - if (plot_deltaY == 1) { + if (plot_deltaY_dir != -1) { nplot += nspecies; } @@ -261,10 +261,10 @@ void WritePlotFileStag(int step, } } - if (plot_deltaY == 1) { + if (plot_deltaY_dir != -1) { MultiFab Ybar(ba, dmap, nspecies, 0); // Yk is component 6: in prim - WriteHorizontalAverageToMF(prim,Ybar,project_dir,6,nspecies,0); + 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); @@ -459,7 +459,7 @@ void WritePlotFileStag(int step, } - if (plot_deltaY == 1) { + if (plot_deltaY_dir != 1) { x = "deltaYk_"; for (i=0; i