From 34b1c36e1cfa2df353cc2604f5448696c73c7201 Mon Sep 17 00:00:00 2001 From: Bhargav Sriram Siddani Date: Sun, 3 Nov 2024 23:37:19 -0800 Subject: [PATCH] Clean up related to Inertial Number in rheology/incflo_rheology.cpp --- src/derive/incflo_derive.cpp | 15 ++++---- src/incflo.H | 7 ++-- src/rheology/incflo_rheology.cpp | 59 +++++++++++++++----------------- src/utilities/io.cpp | 30 ++++++++++++---- 4 files changed, 60 insertions(+), 51 deletions(-) diff --git a/src/derive/incflo_derive.cpp b/src/derive/incflo_derive.cpp index 4fb8baa6..67786e18 100644 --- a/src/derive/incflo_derive.cpp +++ b/src/derive/incflo_derive.cpp @@ -228,16 +228,12 @@ void incflo::compute_nodal_hydrostatic_pressure_at_level (int lev, void incflo::compute_nodal_inertial_num_at_level (int lev, MultiFab* inertial_num, - MultiFab* vel, + MultiFab* strainrate, MultiFab* press, Real p_eps, Real ro_grain, Real diam_grain, - Geometry& lev_geom, - Real time, int nghost) + int nghost) { - // Get strainrate in inertial_num - // Strainrate calculated is two times the actual value - compute_nodal_strainrate_at_level(lev,inertial_num,vel,lev_geom,time,nghost); #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -245,6 +241,7 @@ void incflo::compute_nodal_inertial_num_at_level (int lev, { Box const& bx = mfi.growntilebox(nghost); Array4 const& p_nd_arr = press->const_array(mfi); + Array4 const& sr_arr = strainrate->const_array(mfi); Array4 const& inrt_num_arr = inertial_num->array(mfi); const Real eps = p_eps; const Real diam_scnd = diam_grain; @@ -256,9 +253,9 @@ void incflo::compute_nodal_inertial_num_at_level (int lev, + eps*eps); p_reg += p_nd_arr(i,j,k); p_reg *= Real(0.5); - inrt_num_arr(i,j,k) *= - std::sqrt(ro_scnd/p_reg)* - diam_scnd*Real(0.5); + // Strainrate in incflo is two-times the actual value + inrt_num_arr(i,j,k) = std::sqrt(ro_scnd/p_reg)* + diam_scnd*Real(0.5)*sr_arr(i,j,k); }); } } diff --git a/src/incflo.H b/src/incflo.H index b6a6716c..2b03731f 100644 --- a/src/incflo.H +++ b/src/incflo.H @@ -318,14 +318,13 @@ public: amrex::MultiFab* rho_cc, amrex::Geometry& lev_geom, int nghost); - virtual void compute_nodal_inertial_num_at_level (int lev, + void compute_nodal_inertial_num_at_level (int lev, amrex::MultiFab* inertial_num, - amrex::MultiFab* vel, + amrex::MultiFab* strainrate, amrex::MultiFab* press, amrex::Real p_eps, amrex::Real ro_grain, amrex::Real diam_grain, - amrex::Geometry& lev_geom, - amrex::Real time, int nghost); + int nghost); void compute_nodal_second_fluid_conc (amrex::MultiFab* conc_second_nd, amrex::MultiFab* rho_cc, int nghost) const; diff --git a/src/rheology/incflo_rheology.cpp b/src/rheology/incflo_rheology.cpp index f1312a0a..4a14c94f 100644 --- a/src/rheology/incflo_rheology.cpp +++ b/src/rheology/incflo_rheology.cpp @@ -228,41 +228,27 @@ void incflo::compute_nodal_viscosity_at_level (int lev, // 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); - // 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(),2,nghost); -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(sr_mf,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - 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& inrt_num_arr = inertial_num.array(mfi); - const Real eps = m_mu_p_eps_second; - const Real diam_scnd = m_diam_second; - const Real ro_scnd = m_ro_grain_second; - 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) - + eps*eps); - p_reg += p_static_arr(i,j,k); - p_reg *= Real(0.5); - inrt_num_arr(i,j,k,0) = - std::sqrt(ro_scnd/p_reg)* - diam_scnd*Real(0.5)*sr_arr(i,j,k); - }); - } - // Copy concentration - MultiFab::Copy(inertial_num,conc_second_nd,0,1,1,nghost); #ifdef USE_AMREX_MPMD if (m_fluid_model_second == FluidModel::DataDrivenMPMD) { - // Copier send inertial_num - mpmd_copiers_send_lev(inertial_num,0,2,lev); + // 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); + 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 + mpmd_copiers_send_lev(inertial_num_mpmd,0,2,lev); // NOTE: Actual received quantity is stress ratio mpmd_copiers_recv_lev(vel_eta_second,0,1,lev); #ifdef _OPENMP @@ -289,6 +275,15 @@ void incflo::compute_nodal_viscosity_at_level (int lev, } else #endif if (m_fluid_model_second == FluidModel::Rauter) { + // 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); + 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); #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif diff --git a/src/utilities/io.cpp b/src/utilities/io.cpp index 4cc36547..b9011caa 100644 --- a/src/utilities/io.cpp +++ b/src/utilities/io.cpp @@ -728,18 +728,27 @@ void incflo::WritePlotFile() &m_leveldata[lev]->density, Geom(lev),0); + MultiFab strainrate(amrex::convert(mf[lev].boxArray(), + IndexType::TheNodeType().ixType()), + mf[lev].DistributionMap(),1,0); + + compute_nodal_strainrate_at_level(lev,&strainrate, + &m_leveldata[lev]->velocity, + Geom(lev), + m_cur_time, 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, + &strainrate, &p_static, m_mu_p_eps_second, m_ro_grain_second, m_diam_second, - Geom(lev), - m_cur_time, 0); + 0); amrex::average_node_to_cellcenter(mf[lev],icomp,inertial_num,0,1); } pltscaVarsName.push_back("inertial_num"); @@ -767,15 +776,24 @@ void incflo::WritePlotFile() compute_nodal_hydrostatic_pressure_at_level(lev,&p_static, &m_leveldata[lev]->density, Geom(lev),0); + + MultiFab strainrate(amrex::convert(mf[lev].boxArray(), + IndexType::TheNodeType().ixType()), + mf[lev].DistributionMap(),1,0); + + compute_nodal_strainrate_at_level(lev,&strainrate, + &m_leveldata[lev]->velocity, + Geom(lev), + m_cur_time, 0); + compute_nodal_inertial_num_at_level(lev, &mu_I, - &m_leveldata[lev]->velocity, + &strainrate, &p_static, m_mu_p_eps_second, m_ro_grain_second, m_diam_second, - Geom(lev), - m_cur_time, 0); + 0); // Copier send inertial_num mpmd_copiers_send_lev(mu_I,0,1,lev); // Receive mu(I) = stress ratio