From a89aed65bc2890eb4fc2520bbac2238dc5aa649b Mon Sep 17 00:00:00 2001 From: Bhargav Sriram Siddani Date: Wed, 30 Oct 2024 16:10:25 -0700 Subject: [PATCH] Leveraging solved pressure instead of hydrostatic pressure in mu(I) --- src/incflo.H | 1 + src/incflo_apply_corrector.cpp | 3 +- src/incflo_apply_predictor.cpp | 1 + src/incflo_utils.cpp | 10 +++++ src/rheology/incflo_rheology.cpp | 77 +++++++++++++------------------- src/utilities/io.cpp | 1 + 6 files changed, 45 insertions(+), 48 deletions(-) diff --git a/src/incflo.H b/src/incflo.H index 8396827d..3a9cd4ab 100644 --- a/src/incflo.H +++ b/src/incflo.H @@ -861,6 +861,7 @@ private: amrex::Vector get_divtau_new () noexcept; amrex::Vector get_laps_old () noexcept; amrex::Vector get_laps_new () noexcept; + amrex::Vector get_pressure () noexcept; // [[nodiscard]] amrex::Vector get_velocity_old_const () const noexcept; [[nodiscard]] amrex::Vector get_velocity_new_const () const noexcept; diff --git a/src/incflo_apply_corrector.cpp b/src/incflo_apply_corrector.cpp index f7cb706b..e02dc62b 100644 --- a/src/incflo_apply_corrector.cpp +++ b/src/incflo_apply_corrector.cpp @@ -139,7 +139,8 @@ void incflo::ApplyCorrector() // ************************************************************************************* // Compute viscosity / diffusive coefficients // ************************************************************************************* - compute_viscosity(GetVecOfPtrs(vel_eta), get_density_new(), get_velocity_new(), new_time, 1); + compute_viscosity(GetVecOfPtrs(vel_eta), get_density_new(), get_velocity_new(), + get_pressure(), new_time, 1); // Here we create divtau of the (n+1,*) state that was computed in the predictor if ( (m_diff_type == DiffusionType::Explicit) || use_tensor_correction ) diff --git a/src/incflo_apply_predictor.cpp b/src/incflo_apply_predictor.cpp index c53d68c3..437cb33f 100644 --- a/src/incflo_apply_predictor.cpp +++ b/src/incflo_apply_predictor.cpp @@ -130,6 +130,7 @@ void incflo::ApplyPredictor (bool incremental_projection) // ************************************************************************************* compute_viscosity(GetVecOfPtrs(vel_eta), get_density_old(), get_velocity_old(), + get_pressure(), m_cur_time, 1); // ************************************************************************************* diff --git a/src/incflo_utils.cpp b/src/incflo_utils.cpp index 437a0b0e..d4c35235 100644 --- a/src/incflo_utils.cpp +++ b/src/incflo_utils.cpp @@ -364,6 +364,16 @@ void incflo::copy_from_old_to_new_tracer (int lev, IntVect const& ng) } } +Vector incflo::get_pressure () noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->p_nd)); + } + return r; +} + #ifdef USE_AMREX_MPMD void incflo::mpmd_copiers_send_lev (amrex::MultiFab& send_mf, int icomp, int ncomp, int lev){ diff --git a/src/rheology/incflo_rheology.cpp b/src/rheology/incflo_rheology.cpp index 133a1795..1a71eb17 100644 --- a/src/rheology/incflo_rheology.cpp +++ b/src/rheology/incflo_rheology.cpp @@ -50,6 +50,7 @@ struct NonNewtonianViscosity void incflo::compute_viscosity (Vector const& vel_eta, Vector const& rho, Vector const& vel, + Vector const& press, Real time, int nghost) { #ifdef USE_AMREX_MPMD @@ -65,12 +66,16 @@ void incflo::compute_viscosity (Vector const& vel_eta, { if (m_nodal_vel_eta) { compute_nodal_viscosity_at_level(lev, vel_eta[lev], rho[lev], - vel[lev], geom[lev], time, 0); + vel[lev],press[lev], geom[lev], time, 0); } else { compute_viscosity_at_level(lev, vel_eta[lev], rho[lev], vel[lev], geom[lev], time, nghost); } } + + if (m_nodal_vel_eta == 0) { + amrex::ignore_unused>(press); + } } #ifdef AMREX_USE_EB @@ -157,6 +162,7 @@ void incflo::compute_nodal_viscosity_at_level (int /*lev*/, MultiFab* vel_eta, MultiFab* rho, MultiFab* vel, + MultiFab* press, Geometry& lev_geom, Real /*time*/, int nghost) { @@ -449,41 +455,6 @@ 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()); // Inertial Number = diameter*strainrate*sqrt(rho_grain/p) // NOTE: Strain-rate calculated is TWO TIMES the actual value // The second component will carry concentration @@ -495,7 +466,7 @@ void incflo::compute_nodal_viscosity_at_level (int /*lev*/, { Box const& bx = mfi.growntilebox(nghost); Array4 const& sr_arr = sr_mf.const_array(mfi); - Array4 const& p_static_arr = p_static.const_array(mfi); + Array4 const& p_nd_arr = press->const_array(mfi); Array4 const& inrt_num_arr = inertial_num.array(mfi); const Real eps = m_mu_p_eps_second; const Real diam_scnd = m_diam_second; @@ -503,9 +474,9 @@ void incflo::compute_nodal_viscosity_at_level (int /*lev*/, 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) + Real p_reg = std::sqrt(p_nd_arr(i,j,k)*p_nd_arr(i,j,k) + eps*eps); - p_reg += p_static_arr(i,j,k); + p_reg += p_nd_arr(i,j,k); p_reg *= Real(0.5); inrt_num_arr(i,j,k,0) = std::sqrt(ro_scnd/p_reg)* @@ -528,16 +499,22 @@ void incflo::compute_nodal_viscosity_at_level (int /*lev*/, { Box const& bx = mfi.growntilebox(nghost); Array4 const& sr_arr = sr_mf.const_array(mfi); - Array4 const& p_static_arr = p_static.const_array(mfi); + Array4 const& p_nd_arr = press->const_array(mfi); Array4 const& vel_eta_snd_arr = vel_eta_second.array(mfi); - const Real eps = m_mu_sr_eps_second; + const Real eps_p = m_mu_p_eps_second; + const Real eps_sr = m_mu_sr_eps_second; // Note: sr_mf contains TWO TIMES strain rate amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + // Regularized Pressure + Real p_reg = std::sqrt(p_nd_arr(i,j,k)*p_nd_arr(i,j,k) + + eps_p*eps_p); + p_reg += p_nd_arr(i,j,k); + p_reg *= Real(0.5); // Regularized strain rate Real sr_reg = std::sqrt(Real(0.25)*sr_arr(i,j,k)*sr_arr(i,j,k) - + eps*eps); - vel_eta_snd_arr(i,j,k) *= p_static_arr(i,j,k); + + eps_sr*eps_sr); + vel_eta_snd_arr(i,j,k) *= p_reg; vel_eta_snd_arr(i,j,k) /= (Real(2.0)*sr_reg); }); } @@ -552,13 +529,14 @@ void incflo::compute_nodal_viscosity_at_level (int /*lev*/, { Box const& bx = mfi.growntilebox(nghost); Array4 const& sr_arr = sr_mf.const_array(mfi); - Array4 const& p_static_arr = p_static.const_array(mfi); + Array4 const& p_nd_arr = press->const_array(mfi); Array4 const& inrt_num_arr = inertial_num.const_array(mfi); Array4 const& vel_eta_snd_arr = vel_eta_second.array(mfi); const Real mu_1_scnd = m_mu_1_second; const Real mu_2_scnd = m_mu_2_second; const Real I_0_scnd = m_I_0_second; - const Real eps = m_mu_sr_eps_second; + const Real eps_p = m_mu_p_eps_second; + const Real eps_sr = m_mu_sr_eps_second; // Note: sr_mf contains TWO TIMES strain rate // Note: Inertial number in Rauter 2021 (Eq. 2.29) // has an extra factor of 2 @@ -569,10 +547,15 @@ void incflo::compute_nodal_viscosity_at_level (int /*lev*/, vel_eta_snd_arr(i,j,k) *= (mu_2_scnd-mu_1_scnd); vel_eta_snd_arr(i,j,k) += mu_1_scnd; // The above value is stress ratio + // Regularized Pressure + Real p_reg = std::sqrt(p_nd_arr(i,j,k)*p_nd_arr(i,j,k) + + eps_p*eps_p); + p_reg += p_nd_arr(i,j,k); + p_reg *= Real(0.5); // Regularized strain rate Real sr_reg = std::sqrt(Real(0.25)*sr_arr(i,j,k)*sr_arr(i,j,k) - + eps*eps); - vel_eta_snd_arr(i,j,k) *= p_static_arr(i,j,k); + + eps_sr*eps_sr); + vel_eta_snd_arr(i,j,k) *= p_reg; vel_eta_snd_arr(i,j,k) /= (Real(2.0)*sr_reg); }); } diff --git a/src/utilities/io.cpp b/src/utilities/io.cpp index f74f57ed..3cefa1b4 100644 --- a/src/utilities/io.cpp +++ b/src/utilities/io.cpp @@ -611,6 +611,7 @@ void incflo::WritePlotFile() &vel_eta, &m_leveldata[lev]->density, &m_leveldata[lev]->velocity, + &m_leveldata[lev]->p_nd, Geom(lev), m_cur_time, 0); amrex::average_node_to_cellcenter(mf[lev],icomp,vel_eta,0,1);