Skip to content

Commit

Permalink
Merge branch 'AMReX-Fluids:development' into vof
Browse files Browse the repository at this point in the history
  • Loading branch information
@IFDL-WSU
IFDL-WSU authored Aug 17, 2024
2 parents f8911d1 + 94aa2ba commit eabae52
Show file tree
Hide file tree
Showing 8 changed files with 297 additions and 12 deletions.
2 changes: 2 additions & 0 deletions CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 )

Expand Down Expand Up @@ -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 "")
Expand Down
8 changes: 6 additions & 2 deletions src/boundary_conditions/incflo_set_bcs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,13 +59,16 @@ incflo::make_BC_MF(int lev, amrex::Gpu::DeviceVector<amrex::BCRec> const& bcs,
std::string const& field)
{
auto ncomp = static_cast<int>(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<iMultiFab> 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();
Expand All @@ -86,6 +89,7 @@ incflo::make_BC_MF(int lev, amrex::Gpu::DeviceVector<amrex::BCRec> const& bcs,
Box bhi = mfi.growntilebox() & dhi;
Array4<int> 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 {
Expand Down
78 changes: 75 additions & 3 deletions src/convection/incflo_compute_MAC_projected_velocities.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
#include <hydro_utils.H>
#include <hydro_bcs_K.H>
#include <incflo.H>
#include <AMReX_MultiFabUtil.H>

Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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";
}

Expand Down Expand Up @@ -184,6 +186,77 @@ 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<iMultiFab> 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);
time_dep_inflow_vel.setVal(0.);
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<int const> const& velbc_arr = velBC_MF ? (*velBC_MF).const_array(mfi)
: Array4<int const>{};

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
Expand All @@ -205,7 +278,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);

Expand Down
50 changes: 47 additions & 3 deletions src/convection/incflo_compute_advection_term.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -297,12 +297,49 @@ incflo::compute_convective_term (Vector<MultiFab*> 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();
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<iMultiFab> velBC_MF;
if (m_has_mixedBC) {
velBC_MF = make_BC_MF(lev, m_bcrec_velocity_d, "velocity");
Expand All @@ -318,6 +355,10 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
}
}

// ************************************************************************
// Compute advective fluxes
// ************************************************************************
//
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
Expand Down Expand Up @@ -374,6 +415,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
: Array4<int const>{};
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)),
Expand Down Expand Up @@ -416,6 +458,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
: Array4<int const>{};
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)),
Expand Down Expand Up @@ -482,6 +525,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
: Array4<int const>{};
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)),
Expand Down
14 changes: 11 additions & 3 deletions src/prob/prob_bc.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ struct IncfloVelFill
AMREX_GPU_DEVICE
void operator() (const amrex::IntVect& iv, amrex::Array4<amrex::Real> 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
{
Expand All @@ -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)
Expand Down Expand Up @@ -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));
Expand All @@ -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
}
}
Expand Down
13 changes: 12 additions & 1 deletion src/prob/prob_init_fluid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Expand Down Expand Up @@ -1023,6 +1023,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;
Expand Down
72 changes: 72 additions & 0 deletions test_no_eb_2d/benchmark.time_dep_inflow
View file
Original file line number Diff line number Diff line change
@@ -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
Loading

0 comments on commit eabae52

Please sign in to comment.