Skip to content

Commit

Permalink
Clean up related to Inertial Number in rheology/incflo_rheology.cpp
Browse files Browse the repository at this point in the history
  • Loading branch information
@siddanib
siddanib committed Nov 4, 2024
1 parent 3a1e0f5 commit 34b1c36
Show file tree
Hide file tree
Showing 4 changed files with 60 additions and 51 deletions.
15 changes: 6 additions & 9 deletions src/derive/incflo_derive.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -228,23 +228,20 @@ void incflo::compute_nodal_hydrostatic_pressure_at_level (int lev,

void incflo::compute_nodal_inertial_num_at_level (int lev,
MultiFab* inertial_num,
MultiFab* vel,
MultiFab* strainrate,
MultiFab* press,
Real p_eps, Real ro_grain,
Real diam_grain,
Geometry& lev_geom,
Real time, int nghost)
int nghost)
{
// Get strainrate in inertial_num
// Strainrate calculated is two times the actual value
compute_nodal_strainrate_at_level(lev,inertial_num,vel,lev_geom,time,nghost);
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
for (MFIter mfi(*inertial_num,TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
Box const& bx = mfi.growntilebox(nghost);
Array4<Real const> const& p_nd_arr = press->const_array(mfi);
Array4<Real const> const& sr_arr = strainrate->const_array(mfi);
Array4<Real> const& inrt_num_arr = inertial_num->array(mfi);
const Real eps = p_eps;
const Real diam_scnd = diam_grain;
Expand All @@ -256,9 +253,9 @@ void incflo::compute_nodal_inertial_num_at_level (int lev,
+ eps*eps);
p_reg += p_nd_arr(i,j,k);
p_reg *= Real(0.5);
inrt_num_arr(i,j,k) *=
std::sqrt(ro_scnd/p_reg)*
diam_scnd*Real(0.5);
// Strainrate in incflo is two-times the actual value
inrt_num_arr(i,j,k) = std::sqrt(ro_scnd/p_reg)*
diam_scnd*Real(0.5)*sr_arr(i,j,k);
});
}
}
Expand Down
7 changes: 3 additions & 4 deletions src/incflo.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -318,14 +318,13 @@ public:
amrex::MultiFab* rho_cc,
amrex::Geometry& lev_geom,
int nghost);
virtual void compute_nodal_inertial_num_at_level (int lev,
void compute_nodal_inertial_num_at_level (int lev,
amrex::MultiFab* inertial_num,
amrex::MultiFab* vel,
amrex::MultiFab* strainrate,
amrex::MultiFab* press,
amrex::Real p_eps, amrex::Real ro_grain,
amrex::Real diam_grain,
amrex::Geometry& lev_geom,
amrex::Real time, int nghost);
int nghost);
void compute_nodal_second_fluid_conc (amrex::MultiFab* conc_second_nd,
amrex::MultiFab* rho_cc,
int nghost) const;
Expand Down
59 changes: 27 additions & 32 deletions src/rheology/incflo_rheology.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -228,41 +228,27 @@ void incflo::compute_nodal_viscosity_at_level (int lev,
// nodal MultiFab for hydrostatic pressure
MultiFab p_static(vel_eta->boxArray(),vel_eta->DistributionMap(),1,nghost);
compute_nodal_hydrostatic_pressure_at_level(lev,&p_static,rho,lev_geom,nghost);
// Inertial Number = diameter*strainrate*sqrt(rho_grain/p)
// NOTE: Strain-rate calculated is TWO TIMES the actual value
// The second component will carry concentration
MultiFab inertial_num(vel_eta->boxArray(),vel_eta->DistributionMap(),2,nghost);
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
for (MFIter mfi(sr_mf,TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
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& 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)
+ eps*eps);
p_reg += p_static_arr(i,j,k);
p_reg *= Real(0.5);
inrt_num_arr(i,j,k,0) =
std::sqrt(ro_scnd/p_reg)*
diam_scnd*Real(0.5)*sr_arr(i,j,k);
});
}
// Copy concentration
MultiFab::Copy(inertial_num,conc_second_nd,0,1,1,nghost);

#ifdef USE_AMREX_MPMD
if (m_fluid_model_second == FluidModel::DataDrivenMPMD) {
// Copier send inertial_num
mpmd_copiers_send_lev(inertial_num,0,2,lev);
// Inertial Number = diameter*strainrate*sqrt(rho_grain/p)
// NOTE: Strain-rate calculated is TWO TIMES the actual value
// The second component will carry concentration
MultiFab inertial_num(vel_eta->boxArray(),vel_eta->DistributionMap(),
1,nghost);
compute_nodal_inertial_num_at_level(lev,&inertial_num,
&sr_mf,&p_static,m_mu_p_eps_second,
m_ro_grain_second,m_diam_second,
nghost);

MultiFab inertial_num_mpmd(vel_eta->boxArray(),vel_eta->DistributionMap(),
2,nghost);
// Copy concentration
MultiFab::Copy(inertial_num_mpmd,inertial_num,0,0,1,nghost);
// Copy concentration
MultiFab::Copy(inertial_num_mpmd,conc_second_nd,0,1,1,nghost);
// Copier send inertial_num_mpmd
mpmd_copiers_send_lev(inertial_num_mpmd,0,2,lev);
// NOTE: Actual received quantity is stress ratio
mpmd_copiers_recv_lev(vel_eta_second,0,1,lev);
#ifdef _OPENMP
Expand All @@ -289,6 +275,15 @@ void incflo::compute_nodal_viscosity_at_level (int lev,
} else
#endif
if (m_fluid_model_second == FluidModel::Rauter) {
// Inertial Number = diameter*strainrate*sqrt(rho_grain/p)
// NOTE: Strain-rate calculated is TWO TIMES the actual value
// The second component will carry concentration
MultiFab inertial_num(vel_eta->boxArray(),vel_eta->DistributionMap(),
1,nghost);
compute_nodal_inertial_num_at_level(lev,&inertial_num,
&sr_mf,&p_static,m_mu_p_eps_second,
m_ro_grain_second,m_diam_second,
nghost);
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
Expand Down
30 changes: 24 additions & 6 deletions src/utilities/io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -728,18 +728,27 @@ void incflo::WritePlotFile()
&m_leveldata[lev]->density,
Geom(lev),0);

MultiFab strainrate(amrex::convert(mf[lev].boxArray(),
IndexType::TheNodeType().ixType()),
mf[lev].DistributionMap(),1,0);

compute_nodal_strainrate_at_level(lev,&strainrate,
&m_leveldata[lev]->velocity,
Geom(lev),
m_cur_time, 0);

MultiFab inertial_num(amrex::convert(mf[lev].boxArray(),
IndexType::TheNodeType().ixType()),
mf[lev].DistributionMap(),1,0);

compute_nodal_inertial_num_at_level(lev,
&inertial_num,
&m_leveldata[lev]->velocity,
&strainrate,
&p_static,
m_mu_p_eps_second,
m_ro_grain_second,
m_diam_second,
Geom(lev),
m_cur_time, 0);
0);
amrex::average_node_to_cellcenter(mf[lev],icomp,inertial_num,0,1);
}
pltscaVarsName.push_back("inertial_num");
Expand Down Expand Up @@ -767,15 +776,24 @@ void incflo::WritePlotFile()
compute_nodal_hydrostatic_pressure_at_level(lev,&p_static,
&m_leveldata[lev]->density,
Geom(lev),0);

MultiFab strainrate(amrex::convert(mf[lev].boxArray(),
IndexType::TheNodeType().ixType()),
mf[lev].DistributionMap(),1,0);

compute_nodal_strainrate_at_level(lev,&strainrate,
&m_leveldata[lev]->velocity,
Geom(lev),
m_cur_time, 0);

compute_nodal_inertial_num_at_level(lev,
&mu_I,
&m_leveldata[lev]->velocity,
&strainrate,
&p_static,
m_mu_p_eps_second,
m_ro_grain_second,
m_diam_second,
Geom(lev),
m_cur_time, 0);
0);
// Copier send inertial_num
mpmd_copiers_send_lev(mu_I,0,1,lev);
// Receive mu(I) = stress ratio
Expand Down

0 comments on commit 34b1c36

Please sign in to comment.