diff --git a/src/derive/incflo_derive.cpp b/src/derive/incflo_derive.cpp index cc931ec9..935d227b 100644 --- a/src/derive/incflo_derive.cpp +++ b/src/derive/incflo_derive.cpp @@ -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) { @@ -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()) @@ -213,13 +212,30 @@ void incflo::compute_nodal_hydrostatic_pressure_at_level (int lev, 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]); - + // 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 diff --git a/src/incflo.H b/src/incflo.H index 422177b9..a3de11ce 100644 --- a/src/incflo.H +++ b/src/incflo.H @@ -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, @@ -621,31 +622,32 @@ private: amrex::Vector 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 diff --git a/src/rheology/incflo_rheology.cpp b/src/rheology/incflo_rheology.cpp index 4a14c94f..4a97b535 100644 --- a/src/rheology/incflo_rheology.cpp +++ b/src/rheology/incflo_rheology.cpp @@ -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 diff --git a/src/setup/init.cpp b/src/setup/init.cpp index 2fd329ac..dc6010a4 100644 --- a/src/setup/init.cpp +++ b/src/setup/init.cpp @@ -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 diff --git a/src/utilities/io.cpp b/src/utilities/io.cpp index b9011caa..adf94af7 100644 --- a/src/utilities/io.cpp +++ b/src/utilities/io.cpp @@ -428,6 +428,9 @@ void incflo::WritePlotFile() && m_fluid_model_second == FluidModel::DataDrivenMPMD) ++ncomp; #endif + // hydrostatic pressure + if (m_plt_hydrostatic_p) ++ncomp; + Vector mf(finest_level + 1); for (int lev = 0; lev <= finest_level; ++lev) { mf[lev].define(grids[lev], dmap[lev], ncomp, 0, MFInfo(), Factory(lev)); @@ -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(), @@ -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(), @@ -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);