diff --git a/src/diffusion/DiffusionScalarOp.H b/src/diffusion/DiffusionScalarOp.H index ef843825..7a17c69f 100644 --- a/src/diffusion/DiffusionScalarOp.H +++ b/src/diffusion/DiffusionScalarOp.H @@ -32,6 +32,7 @@ public: void compute_laps (amrex::Vector const& laps, amrex::Vector const& a_scalar, amrex::Vector const& eta, + amrex::Vector const& eb_dirichlet, amrex::Vector bcrec); void compute_divtau (amrex::Vector const& a_divtau, diff --git a/src/diffusion/DiffusionScalarOp.cpp b/src/diffusion/DiffusionScalarOp.cpp index 19db4d9e..c78f31ca 100644 --- a/src/diffusion/DiffusionScalarOp.cpp +++ b/src/diffusion/DiffusionScalarOp.cpp @@ -226,8 +226,11 @@ DiffusionScalarOp::diffuse_scalar (Vector 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 @@ -488,6 +491,7 @@ DiffusionScalarOp::diffuse_vel_components (Vector const& vel, void DiffusionScalarOp::compute_laps (Vector const& a_laps, Vector const& a_scalar, Vector const& a_eta, + Vector const& eb_dirichlet, amrex::Vector bcrec) { BL_PROFILE("DiffusionScalarOp::compute_laps"); @@ -534,6 +538,19 @@ void DiffusionScalarOp::compute_laps (Vector 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(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 b = m_incflo->average_scalar_eta_to_faces(lev, eta_comp, *a_eta[lev]); @@ -572,6 +589,8 @@ void DiffusionScalarOp::compute_laps (Vector 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); @@ -635,9 +654,6 @@ void DiffusionScalarOp::compute_divtau (Vector 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); @@ -665,6 +681,16 @@ void DiffusionScalarOp::compute_divtau (Vector 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]); diff --git a/src/diffusion/incflo_diffusion.cpp b/src/diffusion/incflo_diffusion.cpp index cd8c6020..c559023b 100644 --- a/src/diffusion/incflo_diffusion.cpp +++ b/src/diffusion/incflo_diffusion.cpp @@ -13,17 +13,25 @@ incflo::compute_divtau(Vector 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 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 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) @@ -37,7 +45,9 @@ incflo::compute_divtau(Vector 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"; @@ -59,7 +69,7 @@ incflo::compute_laps(Vector const& laps, Vector const& scalar, Vector 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()); } @@ -69,7 +79,7 @@ incflo::compute_laps_T(Vector const& laps, Vector const& scalar, Vector 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()); } diff --git a/src/incflo.H b/src/incflo.H index a2a96850..3950fa87 100644 --- a/src/incflo.H +++ b/src/incflo.H @@ -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; @@ -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 (); @@ -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; } diff --git a/src/rheology/incflo_read_rheology_parameters.cpp b/src/rheology/incflo_read_rheology_parameters.cpp index 358c1fbf..e27963d7 100644 --- a/src/rheology/incflo_read_rheology_parameters.cpp +++ b/src/rheology/incflo_read_rheology_parameters.cpp @@ -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") { diff --git a/src/rheology/incflo_rheology.cpp b/src/rheology/incflo_rheology.cpp index 88d3dea6..fdeef121 100644 --- a/src/rheology/incflo_rheology.cpp +++ b/src/rheology/incflo_rheology.cpp @@ -15,7 +15,7 @@ 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 { @@ -23,7 +23,9 @@ struct NonNewtonianViscosity { 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: { @@ -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); diff --git a/src/setup/init.cpp b/src/setup/init.cpp index f2536970..c7fccb1d 100644 --- a/src/setup/init.cpp +++ b/src/setup/init.cpp @@ -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(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"); }