Skip to content

address issues 175,182,183,185,194 #214

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
asalmgren merged 1 commit into from
Aug 21, 2026
Merged

address issues 175,182,183,185,194 #214

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions src/diffusion/DiffusionScalarOp.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@ public:
void compute_laps (amrex::Vector<amrex::MultiFab*> const& laps,
amrex::Vector<amrex::MultiFab const*> const& a_scalar,
amrex::Vector<amrex::MultiFab const*> const& eta,
amrex::Vector<amrex::MultiFab*> const& eb_dirichlet,
amrex::Vector<amrex::BCRec> bcrec);

void compute_divtau (amrex::Vector<amrex::MultiFab*> const& a_divtau,
Expand Down
36 changes: 31 additions & 5 deletions src/diffusion/DiffusionScalarOp.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -226,8 +226,11 @@ DiffusionScalarOp::diffuse_scalar (Vector<MultiFab*> const& a_scalar,
}

if (!eb_dirichlet[lev]->empty()) {
MultiFab phi(*eb_dirichlet[lev], amrex::make_alias, comp, 1);
m_eb_scal_solve_op->setEBDirichlet(lev, phi, *eta[lev]);
// This op has one component, so the EB beta must be the
// single component of eta that we are solving for
MultiFab phi (*eb_dirichlet[lev], amrex::make_alias, comp, 1);
MultiFab beta(*eta[lev] , amrex::make_alias, comp, 1);
m_eb_scal_solve_op->setEBDirichlet(lev, phi, beta);
} // else use default homogeneous Neumann on EB


Expand Down Expand Up @@ -488,6 +491,7 @@ DiffusionScalarOp::diffuse_vel_components (Vector<MultiFab*> const& vel,
void DiffusionScalarOp::compute_laps (Vector<MultiFab*> const& a_laps,
Vector<MultiFab const*> const& a_scalar,
Vector<MultiFab const*> const& a_eta,
Vector<MultiFab*> const& eb_dirichlet,
amrex::Vector<amrex::BCRec> bcrec)
{
BL_PROFILE("DiffusionScalarOp::compute_laps");
Expand Down Expand Up @@ -534,6 +538,19 @@ void DiffusionScalarOp::compute_laps (Vector<MultiFab*> const& a_laps,
laps_comp.emplace_back(laps_tmp[lev],amrex::make_alias,comp,1);
scalar_comp.emplace_back(*a_scalar[lev],amrex::make_alias,comp,1);

// Use the same EB BC as the implicit solve in diffuse_scalar,
// otherwise the explicit and implicit diffusion terms are
// inconsistent at the EB. NOTE that, exactly as for the solve
// op, this op is shared by the tracer and temperature passes and
// MLEBABecLap has no way to revert an EB Dirichlet BC back to
// Neumann, so if only one of tracer_eb/temperature_eb is defined
// the other pass will see the EB phi left by the first.
if (lev < static_cast<int>(eb_dirichlet.size()) && !eb_dirichlet[lev]->empty()) {
MultiFab phi_eb (*eb_dirichlet[lev], amrex::make_alias, comp , 1);
MultiFab beta_eb(*a_eta[lev] , amrex::make_alias, eta_comp, 1);
m_eb_scal_apply_op->setEBDirichlet(lev, phi_eb, beta_eb);
} // else use default homogeneous Neumann on EB

Array<MultiFab,AMREX_SPACEDIM>
b = m_incflo->average_scalar_eta_to_faces(lev, eta_comp, *a_eta[lev]);

Expand Down Expand Up @@ -572,6 +589,8 @@ void DiffusionScalarOp::compute_laps (Vector<MultiFab*> const& a_laps,
else
#endif
{
amrex::ignore_unused(eb_dirichlet);

// We want to return div (mu grad)) phi
m_reg_scal_apply_op->setScalars(0.0, -1.0);

Expand Down Expand Up @@ -635,9 +654,6 @@ void DiffusionScalarOp::compute_divtau (Vector<MultiFab*> const& a_divtau,
divtau_tmp[lev].setVal(0.0);
}

for (int lev = 0; lev <= finest_level; ++lev) {
m_eb_vel_apply_op->setEBHomogDirichlet(lev, *a_eta[lev]);
}
// We want to return div (mu grad)) phi
m_eb_vel_apply_op->setScalars(0.0, -1.0);

Expand Down Expand Up @@ -665,6 +681,16 @@ void DiffusionScalarOp::compute_divtau (Vector<MultiFab*> const& a_divtau,
divtau_single.emplace_back(divtau_tmp[lev],amrex::make_alias,comp,1);
vel_single.emplace_back( vel[lev],amrex::make_alias,comp,1);

// Match the implicit solve in diffuse_vel_components; the EB
// Dirichlet value is per component, so this must be set inside
// the component loop
if (m_incflo->hasEBFlow()) {
MultiFab phi(*m_incflo->get_velocity_eb()[lev], amrex::make_alias, comp, 1);
m_eb_vel_apply_op->setEBDirichlet(lev, phi, *a_eta[lev]);
} else {
m_eb_vel_apply_op->setEBHomogDirichlet(lev, *a_eta[lev]);
}

if ( m_incflo->m_has_mixedBC ) {
auto const robin = m_incflo->make_robinBC_MFs(lev, &vel_single[lev]);

Expand Down
30 changes: 20 additions & 10 deletions src/diffusion/incflo_diffusion.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,17 +13,25 @@ incflo::compute_divtau(Vector<MultiFab *> const& divtau,

get_diffusion_tensor_op()->compute_divtau(divtau, vel, density, eta);
#ifdef AMREX_USE_EB
EB_set_covered(*divtau[0] , 0.0);
for (int lev = 0; lev <= finest_level; ++lev) {
EB_set_covered(*divtau[lev], 0.0);
}
#endif

Vector<MultiFab*> divtau_scal;
divtau_scal.push_back(new MultiFab(grids[0], dmap[0], divtau[0]->nComp(),
divtau[0]->nGrow(),MFInfo(),*m_factory[0]));
divtau_scal[0]->setVal(0.);
// Note that the diffusion ops loop over all levels, so divtau_scal must be
// defined on all levels, not just level 0.
Vector<MultiFab> divtau_scal(finest_level+1);
for (int lev = 0; lev <= finest_level; ++lev) {
divtau_scal[lev].define(grids[lev], dmap[lev], divtau[lev]->nComp(),
divtau[lev]->nGrow(), MFInfo(), *m_factory[lev]);
divtau_scal[lev].setVal(0.);
}

get_diffusion_scalar_op()->compute_divtau({divtau_scal}, vel, density, eta);
get_diffusion_scalar_op()->compute_divtau(GetVecOfPtrs(divtau_scal), vel, density, eta);
#ifdef AMREX_USE_EB
EB_set_covered(*divtau_scal[0], 0.0);
for (int lev = 0; lev <= finest_level; ++lev) {
EB_set_covered(divtau_scal[lev], 0.0);
}
#endif

// Define divtau to be (divtau_full - divtau_separate)
Expand All @@ -37,7 +45,9 @@ incflo::compute_divtau(Vector<MultiFab *> const& divtau,
// amrex::Print() << "Z-comp: Norm of tensor apply vs scalar apply " <<
// divtau[0]->norm0(2) << " " << divtau_scal[0]->norm0(2) << "\n";

MultiFab::Saxpy(*divtau[0], -1.0, *divtau_scal[0], 0, 0, AMREX_SPACEDIM, 0);
for (int lev = 0; lev <= finest_level; ++lev) {
MultiFab::Saxpy(*divtau[lev], -1.0, divtau_scal[lev], 0, 0, AMREX_SPACEDIM, 0);
}

// amrex::Print() << "X-comp: Norm of difference of tensor apply vs scalar apply " <<
// divtau[0]->norm0(0) << "\n";
Expand All @@ -59,7 +69,7 @@ incflo::compute_laps(Vector<MultiFab *> const& laps,
Vector<MultiFab const*> const& scalar,
Vector<MultiFab const*> const& eta)
{
get_diffusion_scalar_op()->compute_laps(laps, scalar, eta,
get_diffusion_scalar_op()->compute_laps(laps, scalar, eta, get_tracer_eb(),
get_tracer_bcrec());

}
Expand All @@ -69,7 +79,7 @@ incflo::compute_laps_T(Vector<MultiFab *> const& laps,
Vector<MultiFab const*> const& scalar,
Vector<MultiFab const*> const& eta)
{
get_diffusion_scalar_op()->compute_laps(laps, scalar, eta,
get_diffusion_scalar_op()->compute_laps(laps, scalar, eta, get_temperature_eb(),
get_temperature_bcrec());
}

Expand Down
12 changes: 11 additions & 1 deletion src/incflo.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -598,6 +598,9 @@ private:
amrex::Real m_n_0 = 0.0;
amrex::Real m_tau_0 = 0.0;
amrex::Real m_papa_reg = 0.0;
//! Lower bound on the strain rate used by the power-law viscosity, so that
//! eta = mu*sr^(n-1) stays finite (and non-zero) at zero strain rate
amrex::Real m_sr_floor = 1.e-9;
amrex::Real m_eta_0 = 0.0;

int m_plot_int = -1;
Expand Down Expand Up @@ -651,7 +654,7 @@ private:
ParticleData particleData;

// variables and functions for tracers particles
bool m_use_tracer_particles; /*!< tracer particles that advect with flow */
bool m_use_tracer_particles = false; /*!< tracer particles that advect with flow */

/*! Read tracer particles parameters */
void readTracerParticlesParams ();
Expand Down Expand Up @@ -805,6 +808,13 @@ private:
[[nodiscard]] int nghost_mac () const {
#ifdef AMREX_USE_EB
if (!EBFactory(0).isAllRegular()) return (m_advection_type == "MOL") ? 3 : 4;
#endif
#ifdef INCFLO_USE_PARTICLES
// incflo_PC::AdvectWithFlow interpolates from the faces bracketing each
// particle, so it needs at least one ghost face even with MOL, which
// otherwise needs none. Returning a non-zero value here also makes
// compute_convective_term FillPatch those ghost faces.
if (m_use_tracer_particles) return 1;
#endif
return (m_advection_type == "MOL") ? 0 : 1;
}
Expand Down
11 changes: 10 additions & 1 deletion src/rheology/incflo_read_rheology_parameters.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,9 +22,18 @@ void incflo::ReadRheologyParameters()
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_n_0 != 1.0,
"No point in using power-law rheology with n = 1");

// The strain rate is floored at this value before evaluating
// mu*sr^(n-1), which would otherwise be infinite at sr = 0 for n < 1
// (and zero at sr = 0 for n > 1). Note that sr has units of 1/time,
// so the appropriate value is problem-dependent.
pp.query("sr_floor", m_sr_floor);
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_sr_floor > 0.0,
"Power-law strain-rate floor must be positive");

amrex::Print() << "Power-law fluid with"
<< " mu = " << m_mu
<< ", n = " << m_n_0 << "\n";
<< ", n = " << m_n_0
<< ", sr_floor = " << m_sr_floor << "\n";
}
else if(fluid_model_s == "bingham")
{
Expand Down
7 changes: 5 additions & 2 deletions src/rheology/incflo_rheology.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,15 +15,17 @@ amrex::Real expterm (amrex::Real nu) noexcept
struct NonNewtonianViscosity
{
incflo::FluidModel fluid_model;
amrex::Real mu, n_flow, tau_0, eta_0, papa_reg;
amrex::Real mu, n_flow, tau_0, eta_0, papa_reg, sr_floor;

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real operator() (amrex::Real sr) const noexcept {
switch (fluid_model)
{
case incflo::FluidModel::powerlaw:
{
return mu * std::pow(sr,n_flow-Real(1.0));
// Regularize the strain rate so that eta stays finite as sr -> 0 for
// shear-thinning fluids (n < 1), and stays non-zero for n > 1.
return mu * std::pow(amrex::max(sr,sr_floor),n_flow-Real(1.0));
}
case incflo::FluidModel::Bingham:
{
Expand Down Expand Up @@ -82,6 +84,7 @@ void incflo::compute_viscosity_at_level (int /*lev*/,
non_newtonian_viscosity.tau_0 = m_tau_0;
non_newtonian_viscosity.eta_0 = m_eta_0;
non_newtonian_viscosity.papa_reg = m_papa_reg;
non_newtonian_viscosity.sr_floor = m_sr_floor;

#ifdef AMREX_USE_EB
auto const& fact = EBFactory(lev);
Expand Down
7 changes: 7 additions & 0 deletions src/setup/init.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -234,6 +234,13 @@ void incflo::ReadParameters ()
if (m_advect_tracer && m_eb_flow.enabled && m_eb_flow.tracer.empty()) {
Abort("Must specify tracer EB value for flow through EB");
}
// set_eb_tracer reads m_eb_flow.tracer[n] for every n < m_ntrac, and is
// called whenever the list is non-empty -- whether or not eb_flow itself
// is enabled -- so the list must have exactly ntrac entries.
if (!m_eb_flow.tracer.empty() &&
static_cast<int>(m_eb_flow.tracer.size()) != m_ntrac) {
Abort("eb_flow.tracer must have exactly ntrac values");
}
if (m_use_temperature && m_eb_flow.enabled && m_eb_flow.temperature.empty()) {
Abort("Must specify temperature EB value for flow through EB");
}
Expand Down
Loading