Skip to content

Commit

Permalink
Leveraging solved pressure instead of hydrostatic pressure in mu(I)
Browse files Browse the repository at this point in the history
  • Loading branch information
@siddanib
siddanib committed Nov 1, 2024
1 parent 151252d commit bb1d35b
Show file tree
Hide file tree
Showing 6 changed files with 45 additions and 48 deletions.
1 change: 1 addition & 0 deletions src/incflo.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -861,6 +861,7 @@ private:
amrex::Vector<amrex::MultiFab*> get_divtau_new () noexcept;
amrex::Vector<amrex::MultiFab*> get_laps_old () noexcept;
amrex::Vector<amrex::MultiFab*> get_laps_new () noexcept;
amrex::Vector<amrex::MultiFab*> get_pressure () noexcept;
//
[[nodiscard]] amrex::Vector<amrex::MultiFab const*> get_velocity_old_const () const noexcept;
[[nodiscard]] amrex::Vector<amrex::MultiFab const*> get_velocity_new_const () const noexcept;
Expand Down
3 changes: 2 additions & 1 deletion src/incflo_apply_corrector.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,7 +139,8 @@ void incflo::ApplyCorrector()
// *************************************************************************************
// Compute viscosity / diffusive coefficients
// *************************************************************************************
compute_viscosity(GetVecOfPtrs(vel_eta), get_density_new(), get_velocity_new(), new_time, 1);
compute_viscosity(GetVecOfPtrs(vel_eta), get_density_new(), get_velocity_new(),
get_pressure(), new_time, 1);

// Here we create divtau of the (n+1,*) state that was computed in the predictor
if ( (m_diff_type == DiffusionType::Explicit) || use_tensor_correction )
Expand Down
1 change: 1 addition & 0 deletions src/incflo_apply_predictor.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -130,6 +130,7 @@ void incflo::ApplyPredictor (bool incremental_projection)
// *************************************************************************************
compute_viscosity(GetVecOfPtrs(vel_eta),
get_density_old(), get_velocity_old(),
get_pressure(),
m_cur_time, 1);

// *************************************************************************************
Expand Down
10 changes: 10 additions & 0 deletions src/incflo_utils.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -364,6 +364,16 @@ void incflo::copy_from_old_to_new_tracer (int lev, IntVect const& ng)
}
}

Vector<MultiFab*> incflo::get_pressure () noexcept
{
Vector<MultiFab*> r;
r.reserve(finest_level+1);
for (int lev = 0; lev <= finest_level; ++lev) {
r.push_back(&(m_leveldata[lev]->p_nd));
}
return r;
}

#ifdef USE_AMREX_MPMD
void incflo::mpmd_copiers_send_lev (amrex::MultiFab& send_mf,
int icomp, int ncomp, int lev){
Expand Down
77 changes: 30 additions & 47 deletions src/rheology/incflo_rheology.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,7 @@ struct NonNewtonianViscosity
void incflo::compute_viscosity (Vector<MultiFab*> const& vel_eta,
Vector<MultiFab*> const& rho,
Vector<MultiFab*> const& vel,
Vector<MultiFab*> const& press,
Real time, int nghost)
{
#ifdef USE_AMREX_MPMD
Expand All @@ -65,12 +66,16 @@ void incflo::compute_viscosity (Vector<MultiFab*> const& vel_eta,
{
if (m_nodal_vel_eta) {
compute_nodal_viscosity_at_level(lev, vel_eta[lev], rho[lev],
vel[lev], geom[lev], time, 0);
vel[lev],press[lev], geom[lev], time, 0);
} else {
compute_viscosity_at_level(lev, vel_eta[lev], rho[lev],
vel[lev], geom[lev], time, nghost);
}
}

if (m_nodal_vel_eta == 0) {
amrex::ignore_unused<Vector<MultiFab*>>(press);
}
}

#ifdef AMREX_USE_EB
Expand Down Expand Up @@ -157,6 +162,7 @@ void incflo::compute_nodal_viscosity_at_level (int /*lev*/,
MultiFab* vel_eta,
MultiFab* rho,
MultiFab* vel,
MultiFab* press,
Geometry& lev_geom,
Real /*time*/, int nghost)
{
Expand Down Expand Up @@ -449,41 +455,6 @@ void incflo::compute_nodal_viscosity_at_level (int /*lev*/,
}
}

BoxArray pencil_ba(surroundingNodes(lev_geom.Domain()));
#if (AMREX_SPACEDIM==2)
IntVect pencil_iv(8,1048576);
#else
IntVect pencil_iv(8,8,1048576);
#endif
pencil_ba.maxSize(pencil_iv);
DistributionMapping pencil_dm{pencil_ba};
MultiFab pencil_rho_nodal(pencil_ba,pencil_dm,1,nghost);
pencil_rho_nodal.ParallelCopy(rho_nodal,lev_geom.periodicity());
MultiFab pencil_p_static(pencil_ba,pencil_dm,1,nghost);

#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
for (MFIter mfi(pencil_rho_nodal,TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
// Ensure that static pressure is calculated based on validbox
Box const& bx = mfi.tilebox();
Array4<Real> const& p_static_arr = pencil_p_static.array(mfi);
Array4<Real const> const& rho_nodal_arr = pencil_rho_nodal.const_array(mfi);
const Real gravity = std::abs(m_gravity[AMREX_SPACEDIM-1]);

amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{

p_static_arr(i,j,k) = incflo_local_hydrostatic_pressure_nodal(
i,j,k,AMREX_D_DECL(idx,idy,idz),
gravity,rho_nodal_arr,dlo,dhi);
});

}
// nodal MultiFab for hydrostatic pressure
MultiFab p_static(vel_eta->boxArray(),vel_eta->DistributionMap(),1,nghost);
p_static.ParallelCopy(pencil_p_static,lev_geom.periodicity());
// Inertial Number = diameter*strainrate*sqrt(rho_grain/p)
// NOTE: Strain-rate calculated is TWO TIMES the actual value
// The second component will carry concentration
Expand All @@ -495,17 +466,17 @@ void incflo::compute_nodal_viscosity_at_level (int /*lev*/,
{
Box const& bx = mfi.growntilebox(nghost);
Array4<Real const> const& sr_arr = sr_mf.const_array(mfi);
Array4<Real const> const& p_static_arr = p_static.const_array(mfi);
Array4<Real const> const& p_nd_arr = press->const_array(mfi);
Array4<Real> const& inrt_num_arr = inertial_num.array(mfi);
const Real eps = m_mu_p_eps_second;
const Real diam_scnd = m_diam_second;
const Real ro_scnd = m_ro_grain_second;
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
// Regularized Pressure
Real p_reg = std::sqrt(p_static_arr(i,j,k)*p_static_arr(i,j,k)
Real p_reg = std::sqrt(p_nd_arr(i,j,k)*p_nd_arr(i,j,k)
+ eps*eps);
p_reg += p_static_arr(i,j,k);
p_reg += p_nd_arr(i,j,k);
p_reg *= Real(0.5);
inrt_num_arr(i,j,k,0) =
std::sqrt(ro_scnd/p_reg)*
Expand All @@ -528,16 +499,22 @@ void incflo::compute_nodal_viscosity_at_level (int /*lev*/,
{
Box const& bx = mfi.growntilebox(nghost);
Array4<Real const> const& sr_arr = sr_mf.const_array(mfi);
Array4<Real const> const& p_static_arr = p_static.const_array(mfi);
Array4<Real const> const& p_nd_arr = press->const_array(mfi);
Array4<Real> const& vel_eta_snd_arr = vel_eta_second.array(mfi);
const Real eps = m_mu_sr_eps_second;
const Real eps_p = m_mu_p_eps_second;
const Real eps_sr = m_mu_sr_eps_second;
// Note: sr_mf contains TWO TIMES strain rate
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
// Regularized Pressure
Real p_reg = std::sqrt(p_nd_arr(i,j,k)*p_nd_arr(i,j,k)
+ eps_p*eps_p);
p_reg += p_nd_arr(i,j,k);
p_reg *= Real(0.5);
// Regularized strain rate
Real sr_reg = std::sqrt(Real(0.25)*sr_arr(i,j,k)*sr_arr(i,j,k)
+ eps*eps);
vel_eta_snd_arr(i,j,k) *= p_static_arr(i,j,k);
+ eps_sr*eps_sr);
vel_eta_snd_arr(i,j,k) *= p_reg;
vel_eta_snd_arr(i,j,k) /= (Real(2.0)*sr_reg);
});
}
Expand All @@ -552,13 +529,14 @@ void incflo::compute_nodal_viscosity_at_level (int /*lev*/,
{
Box const& bx = mfi.growntilebox(nghost);
Array4<Real const> const& sr_arr = sr_mf.const_array(mfi);
Array4<Real const> const& p_static_arr = p_static.const_array(mfi);
Array4<Real const> const& p_nd_arr = press->const_array(mfi);
Array4<Real const> const& inrt_num_arr = inertial_num.const_array(mfi);
Array4<Real> const& vel_eta_snd_arr = vel_eta_second.array(mfi);
const Real mu_1_scnd = m_mu_1_second;
const Real mu_2_scnd = m_mu_2_second;
const Real I_0_scnd = m_I_0_second;
const Real eps = m_mu_sr_eps_second;
const Real eps_p = m_mu_p_eps_second;
const Real eps_sr = m_mu_sr_eps_second;
// Note: sr_mf contains TWO TIMES strain rate
// Note: Inertial number in Rauter 2021 (Eq. 2.29)
// has an extra factor of 2
Expand All @@ -569,10 +547,15 @@ void incflo::compute_nodal_viscosity_at_level (int /*lev*/,
vel_eta_snd_arr(i,j,k) *= (mu_2_scnd-mu_1_scnd);
vel_eta_snd_arr(i,j,k) += mu_1_scnd;
// The above value is stress ratio
// Regularized Pressure
Real p_reg = std::sqrt(p_nd_arr(i,j,k)*p_nd_arr(i,j,k)
+ eps_p*eps_p);
p_reg += p_nd_arr(i,j,k);
p_reg *= Real(0.5);
// Regularized strain rate
Real sr_reg = std::sqrt(Real(0.25)*sr_arr(i,j,k)*sr_arr(i,j,k)
+ eps*eps);
vel_eta_snd_arr(i,j,k) *= p_static_arr(i,j,k);
+ eps_sr*eps_sr);
vel_eta_snd_arr(i,j,k) *= p_reg;
vel_eta_snd_arr(i,j,k) /= (Real(2.0)*sr_reg);
});
}
Expand Down
1 change: 1 addition & 0 deletions src/utilities/io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -611,6 +611,7 @@ void incflo::WritePlotFile()
&vel_eta,
&m_leveldata[lev]->density,
&m_leveldata[lev]->velocity,
&m_leveldata[lev]->p_nd,
Geom(lev),
m_cur_time, 0);
amrex::average_node_to_cellcenter(mf[lev],icomp,vel_eta,0,1);
Expand Down

0 comments on commit bb1d35b

Please sign in to comment.