From 402f76cba3f324f0a70797ec5327b15a338757a8 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sun, 11 Aug 2024 19:21:08 -0700 Subject: [PATCH 1/2] Add time dep vels 2 (#118) --- CMakeLists.txt | 2 + src/boundary_conditions/incflo_set_bcs.cpp | 8 +- ...ncflo_compute_MAC_projected_velocities.cpp | 77 ++++++++++++++++++- .../incflo_compute_advection_term.cpp | 50 +++++++++++- src/prob/prob_bc.H | 14 +++- src/prob/prob_init_fluid.cpp | 13 +++- test_no_eb_2d/benchmark.time_dep_inflow | 72 +++++++++++++++++ test_no_eb_3d/benchmark.time_dep_inflow | 72 +++++++++++++++++ 8 files changed, 296 insertions(+), 12 deletions(-) create mode 100644 test_no_eb_2d/benchmark.time_dep_inflow create mode 100644 test_no_eb_3d/benchmark.time_dep_inflow diff --git a/CMakeLists.txt b/CMakeLists.txt index 7a2b155df..f8074ceea 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -54,6 +54,7 @@ option( INCFLO_CUDA "Enable CUDA" NO ) option( INCFLO_HIP "Enable HIP" NO ) option( INCFLO_SYCL "Enable SYCL" NO ) option( INCFLO_EB "Build Embedded Boundary support" NO ) +option( INCFLO_PARTICLES "Build particle support" NO ) option( INCFLO_HYPRE "Enable HYPRE" NO ) option( INCFLO_FPE "Enable Floating Point Exceptions checks" NO ) @@ -108,6 +109,7 @@ if (AMREX_HOME) # SUPERBUILD MODE set(AMReX_GPU_BACKEND SYCL CACHE INTERNAL "") endif () set(AMReX_EB ${INCFLO_EB} CACHE INTERNAL "") + set(AMReX_PARTICLES ${INCFLO_PARTICLES} CACHE INTERNAL "") set(AMREX_FORTRAN OFF CACHE INTERNAL "") set(AMReX_LINEAR_SOLVERS ON CACHE INTERNAL "") set(AMReX_BUILD_TUTORIALS OFF CACHE INTERNAL "") diff --git a/src/boundary_conditions/incflo_set_bcs.cpp b/src/boundary_conditions/incflo_set_bcs.cpp index efe172595..b9f16aa8a 100644 --- a/src/boundary_conditions/incflo_set_bcs.cpp +++ b/src/boundary_conditions/incflo_set_bcs.cpp @@ -59,13 +59,16 @@ incflo::make_BC_MF(int lev, amrex::Gpu::DeviceVector const& bcs, std::string const& field) { auto ncomp = static_cast(bcs.size()); + // The advection routines expect that the BC type is stored in the first ghost // cell of a cell-centered MF. std::unique_ptr BC_MF(new iMultiFab(grids[lev], dmap[lev], ncomp, 1)); - // initialize to 0 == BCType::int_dir, not sure this is needed really + + // Initialize to 0 == BCType::int_dir, not sure this is needed really *BC_MF = 0; + // For advection, we use FOEXTRAP for outflow regions per incflo::init_bcs() - int inflow = BCType::ext_dir; + int inflow = BCType::ext_dir; int outflow = BCType::foextrap; Box const& domain = geom[lev].Domain(); @@ -86,6 +89,7 @@ incflo::make_BC_MF(int lev, amrex::Gpu::DeviceVector const& bcs, Box bhi = mfi.growntilebox() & dhi; Array4 const& bc_arr = BC_MF->array(mfi); BCRec const* const bc_ptr = bcs.data(); + if (m_bc_type[olo] == BC::mixed) { prob_set_BC_MF(olo, blo, bc_arr, lev, inflow, outflow, field); } else { diff --git a/src/convection/incflo_compute_MAC_projected_velocities.cpp b/src/convection/incflo_compute_MAC_projected_velocities.cpp index 88025f201..c80edc5f3 100644 --- a/src/convection/incflo_compute_MAC_projected_velocities.cpp +++ b/src/convection/incflo_compute_MAC_projected_velocities.cpp @@ -1,4 +1,5 @@ #include +#include #include #include @@ -58,8 +59,9 @@ incflo::compute_MAC_projected_velocities ( } // end m_godunov_include_diff_in_forcing - if (nghost_force() > 0) + if (nghost_force() > 0) { fillpatch_force(m_cur_time, vel_forces, nghost_force()); + } } // end m_advection_type @@ -143,7 +145,7 @@ incflo::compute_MAC_projected_velocities ( // not this velocity extrapolation. So we'll use Godunov here instead. // std::string l_advection_type = m_advection_type; - if ( l_advection_type == "BDS" ) { + if ( m_advection_type == "BDS" ) { l_advection_type = "Godunov"; } @@ -184,6 +186,76 @@ incflo::compute_MAC_projected_velocities ( mac_vec[lev][2] = w_mac[lev];); } + const auto *bc_vel_d = get_velocity_bcrec_device_ptr(); + + std::unique_ptr velBC_MF; + + // + // This loop is only necessary if there is time-dependent inflow but we don't have a good test for that + // so we do it if there is any inflow face + // + for (int lev=0; lev <= finest_level; ++lev) + { + MultiFab time_dep_inflow_vel(vel[lev]->boxArray(),vel[lev]->DistributionMap(),AMREX_SPACEDIM,1); + fillphysbc_velocity(lev, m_cur_time+0.5*l_dt, time_dep_inflow_vel, 1); + + Box domain(geom[lev].Domain()); + const auto dlo = lbound(domain); + const auto dhi = ubound(domain); + + if (m_has_mixedBC) { + velBC_MF = make_BC_MF(lev, m_bcrec_velocity_d, "velocity"); + } + + for (MFIter mfi(time_dep_inflow_vel,false); mfi.isValid(); ++mfi) + { + Box const& bx = mfi.validbox(); + auto cc_arr = time_dep_inflow_vel.const_array(mfi); + + Array4 const& velbc_arr = velBC_MF ? (*velBC_MF).const_array(mfi) + : Array4{}; + + auto umac_arr = mac_vec[lev][0]->array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + int n = 0; + const auto bc = HydroBC::getBC(i, j, k, n, domain, bc_vel_d, velbc_arr); + if (i == dlo.x && bc.lo(0) == BCType::ext_dir) { + umac_arr(i,j,k) = cc_arr(i-1,j,k,0); + } + if (i == dhi.x && bc.hi(0) == BCType::ext_dir) { + umac_arr(i+1,j,k) = cc_arr(i+1,j,k,0); + } + }); + auto vmac_arr = mac_vec[lev][1]->array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + int n = 1; + const auto bc = HydroBC::getBC(i, j, k, n, domain, bc_vel_d, velbc_arr); + if (j == dlo.y && bc.lo(1) == BCType::ext_dir) { + vmac_arr(i,j,k) = cc_arr(i,j-1,k,1); + } + if (j == dhi.y && bc.hi(1) == BCType::ext_dir) { + vmac_arr(i,j+1,k) = cc_arr(i,j+1,k,1); + } + }); +#if (AMREX_SPACEDIM > 2) + auto wmac_arr = mac_vec[lev][2]->array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + int n = 2; + const auto bc = HydroBC::getBC(i, j, k, n, domain, bc_vel_d, velbc_arr); + if (k == dlo.z && bc.lo(2) == BCType::ext_dir) { + wmac_arr(i,j,k) = cc_arr(i,j,k-1,2); + } + if (k == dhi.z && bc.hi(2) == BCType::ext_dir) { + wmac_arr(i,j,k+1) = cc_arr(i,j,k+1,2); + } + }); +#endif + } // mfi + } // lev + macproj->setUMAC(mac_vec); #ifdef AMREX_USE_EB @@ -205,7 +277,6 @@ incflo::compute_MAC_projected_velocities ( // to change that so we reduce the cost of the MAC projection for (int lev=0; lev <= finest_level; ++lev) mac_phi[lev]->setVal(0.); - //mac_phi[lev]->mult(0.5,0,1,1); macproj->project(mac_phi,m_mac_mg_rtol,m_mac_mg_atol); diff --git a/src/convection/incflo_compute_advection_term.cpp b/src/convection/incflo_compute_advection_term.cpp index 18540c7ed..8f8b8aa39 100644 --- a/src/convection/incflo_compute_advection_term.cpp +++ b/src/convection/incflo_compute_advection_term.cpp @@ -297,12 +297,49 @@ incflo::compute_convective_term (Vector const& conv_u, divu[lev].FillBoundary(geom[lev].periodicity()); + // ************************************************************************************* + // Define domain boundary conditions at half-time to be used for fluxes if using Godunov + // ************************************************************************************* + // + MultiFab vel_nph( vel[lev]->boxArray(), vel[lev]->DistributionMap(),AMREX_SPACEDIM,1); + MultiFab rho_nph(density[lev]->boxArray(),density[lev]->DistributionMap(),1,1); + MultiFab trac_nph( tracer[lev]->boxArray(), tracer[lev]->DistributionMap(),m_ntrac,1); + + if (m_advection_type != "MOL") { + vel_nph.setVal(0.); + fillphysbc_velocity(lev, m_cur_time+0.5*m_dt, vel_nph, 1); + + if ( !m_constant_density || m_advect_momentum || + (m_advect_tracer && m_ntrac > 0) ) + { + rho_nph.setVal(0.); + fillphysbc_density(lev, m_cur_time+0.5*m_dt, rho_nph, 1); + } + + if ( m_advect_momentum ) { + for (int n = 0; n < AMREX_SPACEDIM; n++) { + Multiply(vel_nph, rho_nph, 0, n, 1, 1); + } + } + + if (m_advect_tracer && (m_ntrac>0)) { + trac_nph.setVal(0.); + fillphysbc_tracer(lev, m_cur_time+0.5*m_dt, trac_nph, 1); + auto const* iconserv = get_tracer_iconserv_device_ptr(); + for (int n = 0; n < m_ntrac; n++) { + if ( iconserv[n] ){ + Multiply(trac_nph, rho_nph, 0, n, 1, 1); + } + } + } + } + // ************************************************************************ - // Compute advective fluxes + // Define mixed boundary conditions if relevant // ************************************************************************ // - // Create BC MF - // FIXME -- would it be worth it to save this from vel extrap??? + // Create BC MF first (to hold bc's that vary in space along the face) + // std::unique_ptr velBC_MF; if (m_has_mixedBC) { velBC_MF = make_BC_MF(lev, m_bcrec_velocity_d, "velocity"); @@ -318,6 +355,10 @@ incflo::compute_convective_term (Vector const& conv_u, } } + // ************************************************************************ + // Compute advective fluxes + // ************************************************************************ + // #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -374,6 +415,7 @@ incflo::compute_convective_term (Vector const& conv_u, : Array4{}; HydroUtils::ComputeFluxesOnBoxFromState( bx, ncomp, mfi, (m_advect_momentum) ? rhovel[lev].array(mfi) : vel[lev]->const_array(mfi), + vel_nph.const_array(mfi), AMREX_D_DECL(flux_x[lev].array(mfi,face_comp), flux_y[lev].array(mfi,face_comp), flux_z[lev].array(mfi,face_comp)), @@ -416,6 +458,7 @@ incflo::compute_convective_term (Vector const& conv_u, : Array4{}; HydroUtils::ComputeFluxesOnBoxFromState( bx, ncomp, mfi, density[lev]->const_array(mfi), + rho_nph.const_array(mfi), AMREX_D_DECL(flux_x[lev].array(mfi,face_comp), flux_y[lev].array(mfi,face_comp), flux_z[lev].array(mfi,face_comp)), @@ -482,6 +525,7 @@ incflo::compute_convective_term (Vector const& conv_u, : Array4{}; HydroUtils::ComputeFluxesOnBoxFromState( bx, ncomp, mfi, any_conserv_trac ? trac_tmp : tracer[lev]->const_array(mfi), + trac_nph.const_array(mfi), AMREX_D_DECL(flux_x[lev].array(mfi,face_comp), flux_y[lev].array(mfi,face_comp), flux_z[lev].array(mfi,face_comp)), diff --git a/src/prob/prob_bc.H b/src/prob/prob_bc.H index 8d54d20f4..8cb08cf03 100644 --- a/src/prob/prob_bc.H +++ b/src/prob/prob_bc.H @@ -19,7 +19,7 @@ struct IncfloVelFill AMREX_GPU_DEVICE void operator() (const amrex::IntVect& iv, amrex::Array4 const& vel, const int dcomp, const int num_comp, - amrex::GeometryData const& geom, const amrex::Real /*time*/, + amrex::GeometryData const& geom, const amrex::Real time, const amrex::BCRec* bcr, const int bcomp, const int orig_comp) const { @@ -33,7 +33,6 @@ struct IncfloVelFill #else const int k = 0; #endif - const amrex::Box& domain_box = geom.Domain(); for (int nc = 0; nc < num_comp; ++nc) @@ -80,7 +79,12 @@ struct IncfloVelFill amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1.0)/domain_box.length(1)); vel(i,j,k,dcomp+nc) *= amrex::Real(6.) * y * (amrex::Real(1.0)-y); } -#if (AMREX_SPACEDIM == 3) +#if (AMREX_SPACEDIM == 2) + else if (42 == probtype) + { + vel(i,j,k,dcomp+nc) = time; + } +#elif (AMREX_SPACEDIM == 3) else if (311 == probtype) { amrex::Real z = amrex::Real(k+0.5)*(amrex::Real(1.0)/domain_box.length(2)); @@ -91,6 +95,10 @@ struct IncfloVelFill amrex::Real z = amrex::Real(k+0.5)*(amrex::Real(1.0)/domain_box.length(2)); vel(i,j,k,dcomp+nc) = amrex::Real(0.5) * z; } + else if (42 == probtype) + { + vel(i,j,k,dcomp+nc) = time; + } #endif } } diff --git a/src/prob/prob_init_fluid.cpp b/src/prob/prob_init_fluid.cpp index 08634fa7c..09d307ed6 100644 --- a/src/prob/prob_init_fluid.cpp +++ b/src/prob/prob_init_fluid.cpp @@ -158,7 +158,7 @@ void incflo::prob_init_fluid (int lev) } else if (31 == m_probtype || 32 == m_probtype || 33 == m_probtype || 311 == m_probtype || 322 == m_probtype || 333 == m_probtype || - 41 == m_probtype) + 41 == m_probtype || 42 == m_probtype) { init_plane_poiseuille(vbx, gbx, ld.velocity.array(mfi), @@ -1010,6 +1010,17 @@ void incflo::init_plane_poiseuille (Box const& vbx, Box const& /*gbx*/, if (nt > 2 && i <= dhi.x*3/4) tracer(i,j,k,2) = 3.0; }); } + else if (42 == m_probtype) + { + ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + vel(i,j,k,0) = 0.0; + const int nt = tracer.nComp(); + for (int n = 0; n < nt; ++n) { + tracer(i,j,k,n) = 0.0; + } + }); + } else if (32 == m_probtype) { Real v = m_ic_v; diff --git a/test_no_eb_2d/benchmark.time_dep_inflow b/test_no_eb_2d/benchmark.time_dep_inflow new file mode 100644 index 000000000..17c59848f --- /dev/null +++ b/test_no_eb_2d/benchmark.time_dep_inflow @@ -0,0 +1,72 @@ +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SIMULATION STOP # +#.......................................# +max_step = 10 # Max number of time steps +steady_state = 0 # Steady-state solver? + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TIME STEP COMPUTATION # +#.......................................# +incflo.fixed_dt = 0.01 # Use this constant dt if > 0 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# INPUT AND OUTPUT # +#.......................................# +amr.plot_int = 50 # Steps between plot files +amr.check_int = 1000 # Steps between checkpoint files +amr.restart = "" # Checkpoint to restart from + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# PHYSICS # +#.......................................# +incflo.gravity = 0. 0. 0. # Gravitational force (3D) +incflo.ro_0 = 1. # Reference density + +incflo.fluid_model = "newtonian" # Fluid model (rheology) +incflo.mu = 0. # Dynamic viscosity coefficient + +incflo.advect_tracer = 0 + +incflo.advection_type = "Godunov" + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# ADAPTIVE MESH REFINEMENT # +#.......................................# +amr.n_cell = 64 16 32 # Grid cells at coarsest AMRlevel +amr.max_level = 0 # Max AMR level in hierarchy + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# GEOMETRY # +#.......................................# +geometry.prob_lo = 0. 0. 0. # Lo corner coordinates +geometry.prob_hi = 2. 0.5 1. # Hi corner coordinates + +geometry.is_periodic = 0 1 0 # Periodicity x y z (0/1) + +incflo.delp = 0. 0. 0. # Prescribed (cyclic) pressure gradient + +incflo.probtype = 42 + +# Boundary conditions +xlo.type = "mi" # mass inflow +xlo.tracer = 1.0 # inflow value of tracer +xlo.velocity = 0. 0. 0. # inflow value of tracer + +xhi.type = "po" # pressure outflow +xhi.pressure = 0.0 + +zlo.type = "sw" # No-slip wall +zhi.type = "sw" # Slip wall + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# VERBOSITY # +#.......................................# +amr.verbose = 3 # incflo_level +incflo.verbose = 3 # incflo_level +mac_proj.verbose = 1 # MAC Projector +nodal_proj.verbose = 1 # Nodal Projector + +scalar_diffusion.verbose = 0 # DiffusionEquation +tensor_diffusion.verbose = 0 # DiffusionEquation + +amr.plt_ccse_regtest = 1 diff --git a/test_no_eb_3d/benchmark.time_dep_inflow b/test_no_eb_3d/benchmark.time_dep_inflow new file mode 100644 index 000000000..17c59848f --- /dev/null +++ b/test_no_eb_3d/benchmark.time_dep_inflow @@ -0,0 +1,72 @@ +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SIMULATION STOP # +#.......................................# +max_step = 10 # Max number of time steps +steady_state = 0 # Steady-state solver? + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TIME STEP COMPUTATION # +#.......................................# +incflo.fixed_dt = 0.01 # Use this constant dt if > 0 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# INPUT AND OUTPUT # +#.......................................# +amr.plot_int = 50 # Steps between plot files +amr.check_int = 1000 # Steps between checkpoint files +amr.restart = "" # Checkpoint to restart from + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# PHYSICS # +#.......................................# +incflo.gravity = 0. 0. 0. # Gravitational force (3D) +incflo.ro_0 = 1. # Reference density + +incflo.fluid_model = "newtonian" # Fluid model (rheology) +incflo.mu = 0. # Dynamic viscosity coefficient + +incflo.advect_tracer = 0 + +incflo.advection_type = "Godunov" + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# ADAPTIVE MESH REFINEMENT # +#.......................................# +amr.n_cell = 64 16 32 # Grid cells at coarsest AMRlevel +amr.max_level = 0 # Max AMR level in hierarchy + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# GEOMETRY # +#.......................................# +geometry.prob_lo = 0. 0. 0. # Lo corner coordinates +geometry.prob_hi = 2. 0.5 1. # Hi corner coordinates + +geometry.is_periodic = 0 1 0 # Periodicity x y z (0/1) + +incflo.delp = 0. 0. 0. # Prescribed (cyclic) pressure gradient + +incflo.probtype = 42 + +# Boundary conditions +xlo.type = "mi" # mass inflow +xlo.tracer = 1.0 # inflow value of tracer +xlo.velocity = 0. 0. 0. # inflow value of tracer + +xhi.type = "po" # pressure outflow +xhi.pressure = 0.0 + +zlo.type = "sw" # No-slip wall +zhi.type = "sw" # Slip wall + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# VERBOSITY # +#.......................................# +amr.verbose = 3 # incflo_level +incflo.verbose = 3 # incflo_level +mac_proj.verbose = 1 # MAC Projector +nodal_proj.verbose = 1 # Nodal Projector + +scalar_diffusion.verbose = 0 # DiffusionEquation +tensor_diffusion.verbose = 0 # DiffusionEquation + +amr.plt_ccse_regtest = 1 From 1fbb115360467c29e38f775a17930ea31c889e5c Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 12 Aug 2024 15:03:11 -0700 Subject: [PATCH 2/2] fix ptr (#119) --- src/convection/incflo_compute_advection_term.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/convection/incflo_compute_advection_term.cpp b/src/convection/incflo_compute_advection_term.cpp index 8f8b8aa39..07a64715b 100644 --- a/src/convection/incflo_compute_advection_term.cpp +++ b/src/convection/incflo_compute_advection_term.cpp @@ -325,7 +325,7 @@ incflo::compute_convective_term (Vector const& conv_u, if (m_advect_tracer && (m_ntrac>0)) { trac_nph.setVal(0.); fillphysbc_tracer(lev, m_cur_time+0.5*m_dt, trac_nph, 1); - auto const* iconserv = get_tracer_iconserv_device_ptr(); + auto const& iconserv = get_tracer_iconserv(); for (int n = 0; n < m_ntrac; n++) { if ( iconserv[n] ){ Multiply(trac_nph, rho_nph, 0, n, 1, 1);