Skip to content

Commit

Permalink
WIP - Temperature added to LevelData, boundary conditions read,
Browse files Browse the repository at this point in the history 
some progress on updates elsewhere...
  • Loading branch information
@cgilet
cgilet committed Sep 11, 2024
1 parent f3f129b commit 7834f46
Show file tree
Hide file tree
Showing 13 changed files with 261 additions and 56 deletions.
28 changes: 15 additions & 13 deletions src/boundary_conditions/README.org
View file
Original file line number Diff line number Diff line change
@@ -1,17 +1,19 @@
* fillpatch

| | pi | po | mi | nsw | sw | dir_dep |
|--------+----------+----------+----------+----------------------+---------------------------------------|-----------------------
| v_n | foextrap | foextrap | ext_dir | ext_dir (0) | ext_dir (0) | ext_dir if inflowing |
| | | | | | | foextrap otherwise |
| v_t | foextrap | foextrap | ext_dir | ext_dir (0) | hoextrap | ext_dir if inflowing |
| | | | | | | foextrap otherwise |
| rho | foextrap | foextrap | ext_dir | foextrap | if (advection_type == "BDS") foextrap | ext_dir if inflowing |
| | | | | | else hoextrap | foextrap otherwise |
| tracer | foextrap | foextrap | ext_dir | ext_dir if specified | ext_dir if specified, otherwise | ext_dir if inflowing |
| | | | | foextrap otherwise | if (advection_type == "BDS") foextrap | foextrap otherwise |
| | | | | | else hoextrap | |
| force | foextrap | foextrap | foextrap | foextrap | foextrap | foextrap |
| | pi | po | mi | nsw | sw | dir_dep |
|-------------+----------+----------+----------+------------------------+---------------------------------------|-----------------------
| v_n | foextrap | foextrap | ext_dir | ext_dir (0) | ext_dir (0) | ext_dir if inflowing |
| | | | | | | foextrap otherwise |
| v_t | foextrap | foextrap | ext_dir | ext_dir (0) | hoextrap | ext_dir if inflowing |
| | | | | | | foextrap otherwise |
| rho | foextrap | foextrap | ext_dir | foextrap | if (advection_type == "BDS") foextrap | ext_dir if inflowing |
| | | | | | else hoextrap | foextrap otherwise |
| tracer | foextrap | foextrap | ext_dir | ext_dir if specified | ext_dir if specified, otherwise | ext_dir if inflowing |
| | | | | foextrap otherwise | if (advection_type == "BDS") foextrap | foextrap otherwise |
| | | | | | else hoextrap | |
| temperature | foextrap | foextrap | ext_dir | ext_dir if specified | ext_dir if specified, otherwise | ext_dir if inflowing |
| | | | | reflect_even otherwise | reflect_even | foextrap otherwise |
| force | foextrap | foextrap | foextrap | foextrap | foextrap | foextrap |

* projection

Expand All @@ -32,5 +34,5 @@
|---------+---------+---------+-----------+------------------------+------------------------|------------
| v_n | Neumann | Neumann | Dirichlet | Dirichlet | Dirichlet | Dirichlet |
| v_t | Neumann | Neumann | Dirichlet | Dirichlet | Neumann | Dirichlet |
| tracer | Neumann | Neumann | Dirichlet | Dirchelet if specified | Dirichlet if specified | Dirichlet |
| scalars | Neumann | Neumann | Dirichlet | Dirchelet if specified | Dirichlet if specified | Dirichlet |
| | | | | Neumann otherwise | Neumann otherwise | |
37 changes: 36 additions & 1 deletion src/boundary_conditions/boundary_conditions.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@ void incflo::init_bcs ()
m_bcrec_velocity.resize(AMREX_SPACEDIM);
m_bcrec_density.resize(1);
if (m_ntrac > 0) { m_bcrec_tracer.resize(m_ntrac); }
m_bcrec_temperature.resize(1);

auto f = [this] (std::string const& bcid, Orientation ori)
{
Expand All @@ -21,6 +22,8 @@ void incflo::init_bcs ()
m_bc_velocity[ori][1] = 0.0;,
m_bc_velocity[ori][2] = 0.0;);
m_bc_tracer[ori].resize(m_ntrac,0.0);
// FIXME? do i want something reasonable here, or force users to set inputs?
m_bc_temperature[ori] = 1.0;

ParmParse pp(bcid);
std::string bc_type_in = "null";
Expand All @@ -41,6 +44,7 @@ void incflo::init_bcs ()
m_bcrec_velocity[2].set(ori, BCType::foextrap););
m_bcrec_density[0].set(ori, BCType::foextrap);
for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::foextrap); }
m_bcrec_temperature[0].set(ori, BCType::foextrap);
}
else if (bc_type == "pressure_outflow" || bc_type == "po")
{
Expand All @@ -56,6 +60,7 @@ void incflo::init_bcs ()
m_bcrec_velocity[2].set(ori, BCType::foextrap););
m_bcrec_density[0].set(ori, BCType::foextrap);
for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::foextrap); }
m_bcrec_temperature[0].set(ori, BCType::foextrap);
}
else if (bc_type == "mass_inflow" || bc_type == "mi")
{
Expand All @@ -72,13 +77,15 @@ void incflo::init_bcs ()

pp.query("density", m_bc_density[ori]);
pp.queryarr("tracer", m_bc_tracer[ori], 0, m_ntrac);
pp.query("temperature", m_bc_temperature[ori]);

// Set mathematical BCs
AMREX_D_TERM(m_bcrec_velocity[0].set(ori, BCType::ext_dir);,
m_bcrec_velocity[1].set(ori, BCType::ext_dir);,
m_bcrec_velocity[2].set(ori, BCType::ext_dir););
m_bcrec_density[0].set(ori, BCType::ext_dir);
for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::ext_dir); }
m_bcrec_temperature[0].set(ori, BCType::ext_dir);
}
else if (bc_type == "direction_dependent" || bc_type == "dd" )
{
Expand All @@ -97,12 +104,14 @@ void incflo::init_bcs ()

pp.query("density", m_bc_density[ori]);
pp.queryarr("tracer", m_bc_tracer[ori], 0, m_ntrac);
pp.query("temperature", m_bc_temperature[ori]);

AMREX_D_TERM(m_bcrec_velocity[0].set(ori, BCType::direction_dependent);,
m_bcrec_velocity[1].set(ori, BCType::direction_dependent);,
m_bcrec_velocity[2].set(ori, BCType::direction_dependent););
m_bcrec_density[0].set(ori, BCType::direction_dependent);
for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::direction_dependent); }
m_bcrec_temperature[0].set(ori, BCType::direction_dependent);
}
else if (bc_type == "no_slip_wall" || bc_type == "nsw")
{
Expand All @@ -126,6 +135,7 @@ void incflo::init_bcs ()
// We potentially read in values at no-slip walls in the event that the
// tracer has Dirichlet bcs
pp.queryarr("tracer", m_bc_tracer[ori], 0, m_ntrac);
pp.query("temperature", m_bc_temperature[ori], 0, 1);

// Set mathematical BCs
AMREX_D_TERM(m_bcrec_velocity[0].set(ori, BCType::ext_dir);,
Expand All @@ -137,6 +147,11 @@ void incflo::init_bcs ()
} else {
for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::foextrap); }
}
if ( pp.contains("temperature") ) {
for (auto& b : m_bcrec_temperature) { b.set(ori, BCType::ext_dir); }
} else {
for (auto& b : m_bcrec_temperature) { b.set(ori, BCType::reflect_even); }
}
}
else if (bc_type == "slip_wall" || bc_type == "sw")
{
Expand All @@ -151,6 +166,7 @@ void incflo::init_bcs ()
// We potentially read in values at slip walls in the event that the
// tracer has Dirichlet bcs
pp.queryarr("tracer", m_bc_tracer[ori], 0, m_ntrac);
pp.query("temperature", m_bc_temperature[ori], 0, 1);

// Tangential directions have hoextrap
AMREX_D_TERM(m_bcrec_velocity[0].set(ori, BCType::hoextrap);,
Expand All @@ -175,6 +191,11 @@ void incflo::init_bcs ()
}
}
}
if ( pp.contains("temperature") ) {
for (auto& b : m_bcrec_temperature) { b.set(ori, BCType::ext_dir); }
} else {
for (auto& b : m_bcrec_temperature) { b.set(ori, BCType::reflect_even); }
}
}
else if (bc_type == "mixed" )
{
Expand Down Expand Up @@ -208,13 +229,15 @@ void incflo::init_bcs ()

pp.query("density", m_bc_density[ori]);
pp.queryarr("tracer", m_bc_tracer[ori], 0, m_ntrac);
pp.query("temperature", m_bc_temperature[ori]);

// Set mathematical BCs. BC_mask will handle Dirichlet part.
AMREX_D_TERM(m_bcrec_velocity[0].set(ori, BCType::foextrap);,
m_bcrec_velocity[1].set(ori, BCType::foextrap);,
m_bcrec_velocity[2].set(ori, BCType::foextrap););
m_bcrec_density[0].set(ori, BCType::foextrap);
for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::foextrap); }
m_bcrec_temperature[0].set(ori, BCType::foextrap);
}
else
{
Expand All @@ -231,6 +254,7 @@ void incflo::init_bcs ()
m_bcrec_velocity[2].set(ori, BCType::int_dir););
m_bcrec_density[0].set(ori, BCType::int_dir);
for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::int_dir); }
m_bcrec_temperature[0].set(ori, BCType::int_dir);
} else {
amrex::Abort("Wrong BC type for periodic boundary");
}
Expand Down Expand Up @@ -285,7 +309,8 @@ void incflo::init_bcs ()
#else
std::memcpy
#endif
(m_bcrec_density_d.data(), m_bcrec_density.data(), sizeof(BCRec));
(m_bcrec_density_d.data(), m_bcrec_density.data(), sizeof(BCRec));
}

if (m_ntrac > 0)
{
Expand All @@ -298,6 +323,16 @@ void incflo::init_bcs ()
(m_bcrec_tracer_d.data(), m_bcrec_tracer.data(), sizeof(BCRec)*m_ntrac);
}

if (m_use_temperature) {
m_bcrec_temperature_d.resize(1);
#ifdef AMREX_USE_GPU
Gpu::htod_memcpy
#else
std::memcpy
#endif
(m_bcrec_temperature_d.data(), m_bcrec_temperature.data(), sizeof(BCRec));
}

// force
{
const int ncomp = std::max(m_ntrac, AMREX_SPACEDIM);
Expand Down
67 changes: 65 additions & 2 deletions src/convection/incflo_compute_advection_term.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,15 +46,18 @@ void incflo::init_advection ()
void
incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
Vector<MultiFab*> const& conv_r,
Vector<MultiFab*> const& conv_t,
Vector<MultiFab*> const& conv_tra,
Vector<MultiFab*> const& conv_tem,
Vector<MultiFab const*> const& vel,
Vector<MultiFab const*> const& density,
Vector<MultiFab const*> const& tracer,
Vector<MultiFab const*> const& temperature,
AMREX_D_DECL(Vector<MultiFab*> const& u_mac,
Vector<MultiFab*> const& v_mac,
Vector<MultiFab*> const& w_mac),
Vector<MultiFab* > const& vel_forces,
Vector<MultiFab* > const& tra_forces,
Vector<MultiFab* > const& tem_forces,
Real /*time*/)
{
bool fluxes_are_area_weighted = false;
Expand Down Expand Up @@ -191,6 +194,17 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
fillpatch_force(m_cur_time, tra_forces, nghost_force());
}

if (m_use_temperature)
{
compute_T_forces(tem_forces, get_density_old_const());
if (m_godunov_include_diff_in_forcing)
for (int lev = 0; lev <= finest_level; ++lev)
// FIXME? Saxpy? with 1/rhocp - could be as low as point function or by level?
MultiFab::Add(*tem_forces[lev], m_leveldata[lev]->laps_T_o, 0, 0, 1, 0);
if (nghost_force() > 0)
fillpatch_force(m_cur_time, tem_forces, nghost_force());
}

} // end m_advection_type

for (int lev = 0; lev <= finest_level; ++lev)
Expand Down Expand Up @@ -300,7 +314,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
// *************************************************************************************
// Define domain boundary conditions at half-time to be used for fluxes if using Godunov
// *************************************************************************************
//
// For time-dependent Dirichlet BCs
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);
Expand Down Expand Up @@ -354,6 +368,13 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
tracBC_MF = make_BC_MF(lev, m_bcrec_tracer_d, "tracer");
}
}
// FIXME? do all the scalars really need individual BC MFs?
// Need to think about general mixed BC vs inflow-outflow only?
std::unique_ptr<iMultiFab> tempBC_MF;
if (m_use_temperature && m_has_mixedBC) {
Abort("Temperature equation with mixed BC not completed yet");
tempBC_MF = make_BC_MF(lev, m_bcrec_temperature_d, "temperature");
}

// ************************************************************************
// Compute advective fluxes
Expand Down Expand Up @@ -551,6 +572,48 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
m_advection_type, PPM::default_limiter,
allow_inflow_on_outflow, tracbc_arr);
}

// ************************************************************************
// Temperature
// ************************************************************************
if (m_use_temperature) {
// Temperature adveciton is always non-conservative

// FIXME- check this!!!
face_comp += m_ntrac;
ncomp = 1;
is_velocity = false;
allow_inflow_on_outflow = false;
Array4<int const> const& tempbc_arr = tempBC_MF ? (*tempBC_MF).const_array(mfi)
: Array4<int const>{};
HydroUtils::ComputeFluxesOnBoxFromState( bx, ncomp, mfi,
temperature[lev]->const_array(mfi),
temp_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)),
AMREX_D_DECL(face_x[lev].array(mfi,face_comp),
face_y[lev].array(mfi,face_comp),
face_z[lev].array(mfi,face_comp)),
knownFaceStates,
AMREX_D_DECL(u_mac[lev]->const_array(mfi),
v_mac[lev]->const_array(mfi),
w_mac[lev]->const_array(mfi)),
divu_arr,
(!tem_forces.empty()) ? tem_forces[lev]->const_array(mfi) : Array4<Real const>{},
geom[lev], m_dt,
get_temperature_bcrec(),
get_temperature_bcrec_device_ptr(),
get_temperature_iconserv_device_ptr(),
#ifdef AMREX_USE_EB
ebfact,
m_eb_flow.enabled ? get_temperature_eb()[lev]->const_array(mfi) : Array4<Real const>{},
#endif
m_godunov_ppm, m_godunov_use_forces_in_trans,
is_velocity, fluxes_are_area_weighted,
m_advection_type, PPM::default_limiter,
allow_inflow_on_outflow, tempbc_arr);
}
} // mfi
} // lev

Expand Down
3 changes: 2 additions & 1 deletion src/diffusion/DiffusionScalarOp.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,9 +15,10 @@ class DiffusionScalarOp
public:
DiffusionScalarOp (incflo* a_incflo);

void diffuse_scalar (amrex::Vector<amrex::MultiFab*> const& tracer,
void diffuse_scalar (amrex::Vector<amrex::MultiFab*> const& a_scalar,
amrex::Vector<amrex::MultiFab*> const& density,
amrex::Vector<amrex::MultiFab const*> const& eta,
amrex::Vector<int> const& iconserv,
amrex::Real dt);

void diffuse_vel_components (
Expand Down
Loading

0 comments on commit 7834f46

Please sign in to comment.