diff --git a/src/derive/incflo_derive.cpp b/src/derive/incflo_derive.cpp index ad61596e..48a72622 100644 --- a/src/derive/incflo_derive.cpp +++ b/src/derive/incflo_derive.cpp @@ -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 const& rho_arr = rho_cc->const_array(mfi); + Array4 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 const& p_static_arr = pencil_p_static.array(mfi); + Array4 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, diff --git a/src/incflo.H b/src/incflo.H index cda8d6ce..831bc7e2 100644 --- a/src/incflo.H +++ b/src/incflo.H @@ -313,6 +313,11 @@ public: amrex::Geometry& lev_geom, amrex::Real time, int nghost); void compute_tracer_diff_coeff (amrex::Vector 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, diff --git a/src/rheology/incflo_rheology.cpp b/src/rheology/incflo_rheology.cpp index 133a1795..848d750c 100644 --- a/src/rheology/incflo_rheology.cpp +++ b/src/rheology/incflo_rheology.cpp @@ -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 const& p_static_arr = pencil_p_static.array(mfi); - Array4 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 diff --git a/src/utilities/io.cpp b/src/utilities/io.cpp index f74f57ed..a6a92bab 100644 --- a/src/utilities/io.cpp +++ b/src/utilities/io.cpp @@ -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, @@ -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,