Skip to content

Commit

Permalink
Changes to hydrostatic pressure function - Leading to a change of O(1…
Browse files Browse the repository at this point in the history 
…e-6) compared to previous function
  • Loading branch information
@siddanib
siddanib committed Nov 4, 2024
1 parent 0d1246e commit 6eac46f
Show file tree
Hide file tree
Showing 5 changed files with 75 additions and 38 deletions.
28 changes: 22 additions & 6 deletions src/derive/incflo_derive.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -161,6 +161,7 @@ void incflo::compute_nodal_strainrate_at_level (int /*lev*/,
void incflo::compute_nodal_hydrostatic_pressure_at_level (int lev,
MultiFab* p_static,
MultiFab* rho_cc,
Real p_surface,
Geometry& lev_geom,
int nghost)
{
Expand Down Expand Up @@ -200,8 +201,6 @@ void incflo::compute_nodal_hydrostatic_pressure_at_level (int lev,
#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())
Expand All @@ -213,13 +212,30 @@ void incflo::compute_nodal_hydrostatic_pressure_at_level (int lev,
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]);

// Even for pencil_ba, this needs to be based on validbox
const Dim3 v_bxlo = amrex::lbound(mfi.validbox());
const Dim3 v_bxhi = amrex::ubound(mfi.validbox());
const int level = lev;
const Real p_srf = p_surface;
#if (AMREX_SPACEDIM == 2)
int h_end = v_bxhi.y;
#else
int h_end = v_bxhi.z;
#endif
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{

p_static_arr(i,j,k) = incflo_local_hydrostatic_pressure_nodal(
Real p_validbox_top = p_srf;
if (level > 0) {
#if (AMREX_SPACEDIM == 2)
p_validbox_top = p_static_arr(i,h_end,k);
#else
p_validbox_top = p_static_arr(i,j,h_end);
#endif
}
p_static_arr(i,j,k) = p_validbox_top;
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);
gravity,rho_nodal_arr,v_bxlo,v_bxhi);
});
}
// nodal MultiFab for hydrostatic pressure
Expand Down
50 changes: 26 additions & 24 deletions src/incflo.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -316,6 +316,7 @@ public:
void compute_nodal_hydrostatic_pressure_at_level (int lev,
amrex::MultiFab* p_static,
amrex::MultiFab* rho_cc,
amrex::Real p_surface,
amrex::Geometry& lev_geom,
int nghost);
static void compute_nodal_inertial_num_at_level (int lev,
Expand Down Expand Up @@ -621,31 +622,32 @@ private:
amrex::Vector<amrex::Real> tag_region_hi;

// Flags for saving fluid data in plot files
int m_plt_velx = 1;
int m_plt_vely = 1;
int m_plt_velz = 1;
int m_plt_gpx = 1;
int m_plt_gpy = 1;
int m_plt_gpz = 1;
int m_plt_rho = 1;
int m_plt_tracer = 1;
int m_plt_p_nd = 0;
int m_plt_macphi = 0;
int m_plt_eta = 0;
int m_plt_magvel = 1;
int m_plt_vort = 1;
int m_plt_forcing = 0;
int m_plt_strainrate = 0;
int m_plt_divu = 0;
int m_plt_vfrac = 1;
int m_plt_error_u = 0;
int m_plt_error_v = 0;
int m_plt_error_w = 0;
int m_plt_error_p = 0;
int m_plt_error_mac_p = 0;
int m_plt_inertial_num = 0;
int m_plt_velx = 1;
int m_plt_vely = 1;
int m_plt_velz = 1;
int m_plt_gpx = 1;
int m_plt_gpy = 1;
int m_plt_gpz = 1;
int m_plt_rho = 1;
int m_plt_tracer = 1;
int m_plt_p_nd = 0;
int m_plt_macphi = 0;
int m_plt_eta = 0;
int m_plt_magvel = 1;
int m_plt_vort = 1;
int m_plt_forcing = 0;
int m_plt_strainrate = 0;
int m_plt_divu = 0;
int m_plt_vfrac = 1;
int m_plt_error_u = 0;
int m_plt_error_v = 0;
int m_plt_error_w = 0;
int m_plt_error_p = 0;
int m_plt_error_mac_p = 0;
int m_plt_inertial_num = 0;
int m_plt_hydrostatic_p = 0;
#ifdef USE_AMREX_MPMD
int m_plt_mu_I = 0;
int m_plt_mu_I = 0;
#endif

#ifdef INCFLO_USE_PARTICLES
Expand Down
13 changes: 5 additions & 8 deletions src/rheology/incflo_rheology.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -227,24 +227,21 @@ void incflo::compute_nodal_viscosity_at_level (int lev,
compute_nodal_strainrate_at_level(lev,&sr_mf,vel,lev_geom,time,nghost);
// 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);
compute_nodal_hydrostatic_pressure_at_level(lev,&p_static,rho,Real(0.0),
lev_geom,nghost);

#ifdef USE_AMREX_MPMD
if (m_fluid_model_second == FluidModel::DataDrivenMPMD) {
MultiFab inertial_num_mpmd(vel_eta->boxArray(),vel_eta->DistributionMap(),
2,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(),
1,nghost);
MultiFab inertial_num(inertial_num_mpmd,amrex::make_alias,0,1);
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
Expand Down
1 change: 1 addition & 0 deletions src/setup/init.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -300,6 +300,7 @@ void incflo::ReadIOParameters()
pp.query("plt_error_p", m_plt_error_p );
pp.query("plt_error_mac_p",m_plt_error_mac_p );
pp.query("plt_inertial_num", m_plt_inertial_num);
pp.query("plt_hydrostatic_p", m_plt_hydrostatic_p);
#ifdef USE_AMREX_MPMD
pp.query("plt_mu_I", m_plt_mu_I);
#endif
Expand Down
21 changes: 21 additions & 0 deletions src/utilities/io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -428,6 +428,9 @@ void incflo::WritePlotFile()
&& m_fluid_model_second == FluidModel::DataDrivenMPMD) ++ncomp;
#endif

// hydrostatic pressure
if (m_plt_hydrostatic_p) ++ncomp;

Vector<MultiFab> mf(finest_level + 1);
for (int lev = 0; lev <= finest_level; ++lev) {
mf[lev].define(grids[lev], dmap[lev], ncomp, 0, MFInfo(), Factory(lev));
Expand Down Expand Up @@ -726,6 +729,7 @@ void incflo::WritePlotFile()

compute_nodal_hydrostatic_pressure_at_level(lev,&p_static,
&m_leveldata[lev]->density,
Real(0.0),
Geom(lev),0);

MultiFab strainrate(amrex::convert(mf[lev].boxArray(),
Expand Down Expand Up @@ -775,6 +779,7 @@ void incflo::WritePlotFile()

compute_nodal_hydrostatic_pressure_at_level(lev,&p_static,
&m_leveldata[lev]->density,
Real(0.0),
Geom(lev),0);

MultiFab strainrate(amrex::convert(mf[lev].boxArray(),
Expand Down Expand Up @@ -805,6 +810,22 @@ void incflo::WritePlotFile()
}
#endif

if (m_plt_hydrostatic_p) {
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_at_level(lev,&p_static,
&m_leveldata[lev]->density,
Real(0.0),
Geom(lev),0);
amrex::average_node_to_cellcenter(mf[lev],icomp,p_static,0,1);
}
pltscaVarsName.push_back("hydrostatic_p");
++icomp;
}

#ifdef AMREX_USE_EB
for (int lev = 0; lev <= finest_level; ++lev) {
EB_set_covered(mf[lev], 0.0);
Expand Down

0 comments on commit 6eac46f

Please sign in to comment.