Skip to content

Commit

Permalink
Changes to leverage hydrostatic pressure in plotfile(s)
Browse files Browse the repository at this point in the history
  • Loading branch information
@siddanib
siddanib committed Nov 3, 2024
1 parent a0ac0d0 commit 26efd58
Show file tree
Hide file tree
Showing 4 changed files with 92 additions and 35 deletions.
68 changes: 68 additions & 0 deletions src/derive/incflo_derive.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,6 +158,74 @@ void incflo::compute_nodal_strainrate_at_level (int /*lev*/,
}
}

void incflo::compute_nodal_hydrostatic_pressure (int lev,
MultiFab* p_static,
MultiFab* rho_cc,
Geometry& lev_geom,
int nghost)
{
if (lev > 0) {
amrex::Abort("Hydrostatic pressure is not implemented for lev > 0");
}

MultiFab rho_nodal(p_static->boxArray(), p_static->DistributionMap(),1,nghost);
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
for (MFIter mfi(rho_nodal,TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
Box const& bx = mfi.tilebox();
Array4<Real const> const& rho_arr = rho_cc->const_array(mfi);
Array4<Real> const& rho_nodal_arr = rho_nodal.array(mfi);
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
rho_nodal_arr(i,j,k) = incflo_nodal_density(i,j,k,rho_arr);
});
}

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);

Real idx = Real(1.0) / lev_geom.CellSize(0);
Real idy = Real(1.0) / lev_geom.CellSize(1);
#if (AMREX_SPACEDIM == 3)
Real idz = Real(1.0) / lev_geom.CellSize(2);
#endif
const Dim3 dlo = amrex::lbound(lev_geom.Domain());
const Dim3 dhi = amrex::ubound(lev_geom.Domain());

#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
p_static->ParallelCopy(pencil_p_static,lev_geom.periodicity());
}

void incflo::compute_nodal_inertial_num_at_level (int lev,
MultiFab* inertial_num,
MultiFab* vel,
Expand Down
5 changes: 5 additions & 0 deletions src/incflo.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -313,6 +313,11 @@ public:
amrex::Geometry& lev_geom,
amrex::Real time, int nghost);
void compute_tracer_diff_coeff (amrex::Vector<amrex::MultiFab*> const& tra_eta, int nghost);
void compute_nodal_hydrostatic_pressure (int lev,
amrex::MultiFab* p_static,
amrex::MultiFab* rho_cc,
amrex::Geometry& lev_geom,
int nghost);
virtual void compute_nodal_inertial_num_at_level (int lev,
amrex::MultiFab* inertial_num,
amrex::MultiFab* vel,
Expand Down
34 changes: 1 addition & 33 deletions src/rheology/incflo_rheology.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -449,41 +449,9 @@ 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());
compute_nodal_hydrostatic_pressure(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
Expand Down
20 changes: 18 additions & 2 deletions src/utilities/io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -720,13 +720,21 @@ void incflo::WritePlotFile()
&& (m_fluid_model_second == FluidModel::DataDrivenMPMD
|| m_fluid_model_second == FluidModel::Rauter)) {
for (int lev = 0; lev <= finest_level; ++lev) {
MultiFab p_static(amrex::convert(mf[lev].boxArray(),
IndexType::TheNodeType().ixType()),
mf[lev].DistributionMap(),1,0);

compute_nodal_hydrostatic_pressure(lev,&p_static,
&m_leveldata[lev]->density,
Geom(lev),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,
&m_leveldata[lev]->p_nd,
&p_static,
m_mu_p_eps_second,
m_ro_grain_second,
m_diam_second,
Expand All @@ -751,10 +759,18 @@ void incflo::WritePlotFile()
MultiFab mu_I(amrex::convert(mf[lev].boxArray(),
IndexType::TheNodeType().ixType()),
mf[lev].DistributionMap(),1,0);

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

compute_nodal_hydrostatic_pressure(lev,&p_static,
&m_leveldata[lev]->density,
Geom(lev),0);
compute_nodal_inertial_num_at_level(lev,
&mu_I,
&m_leveldata[lev]->velocity,
&m_leveldata[lev]->p_nd,
&p_static,
m_mu_p_eps_second,
m_ro_grain_second,
m_diam_second,
Expand Down

0 comments on commit 26efd58

Please sign in to comment.