diff --git a/src/prob/incflo_prob_I.H b/src/prob/incflo_prob_I.H index 0cb105c9..7f035e76 100644 --- a/src/prob/incflo_prob_I.H +++ b/src/prob/incflo_prob_I.H @@ -164,12 +164,4 @@ amrex::GpuArray const& dx, amrex::GpuArray const& problo, amrex::GpuArray const& probhi) const; - - void column_collapse_granular (amrex::Box const& vbx, amrex::Box const& nbx, - amrex::Array4 const& density, - amrex::Array4 const& pressure, - amrex::Box const& domain, - amrex::GpuArray const& dx, - amrex::GpuArray const& problo, - amrex::GpuArray const& probhi); #endif diff --git a/src/prob/incflo_prob_usr_I.H b/src/prob/incflo_prob_usr_I.H index fcc1025f..77d3b090 100644 --- a/src/prob/incflo_prob_usr_I.H +++ b/src/prob/incflo_prob_usr_I.H @@ -131,4 +131,21 @@ amrex::GpuArray const& dx, amrex::GpuArray const& problo, amrex::GpuArray const& probhi); + + void column_collapse_granular (amrex::Box const& vbx, amrex::Box const& nbx, + amrex::Array4 const& density, + amrex::Array4 const& pressure, + amrex::Box const& domain, + amrex::GpuArray const& dx, + amrex::GpuArray const& problo, + amrex::GpuArray const& probhi); + + void smooth_column_collapse_granular (amrex::Box const& vbx, amrex::Box const& nbx, + amrex::Array4 const& density, + amrex::Array4 const& pressure, + amrex::Box const& domain, + amrex::GpuArray const& dx, + amrex::GpuArray const& problo, + amrex::GpuArray const& probhi, + amrex::Real smoothing_factor); #endif diff --git a/src/prob/prob_init_fluid.cpp b/src/prob/prob_init_fluid.cpp index 80dc02b6..d15a831f 100644 --- a/src/prob/prob_init_fluid.cpp +++ b/src/prob/prob_init_fluid.cpp @@ -193,6 +193,18 @@ void incflo::prob_init_fluid (int lev) ld.p_nd.array(mfi), domain, dx, problo, probhi); } + else if (532 == m_probtype) + { + Real smoothing_factor; + amrex::ParmParse pp_tf("two_fluid"); + pp_tf.get("init_smoothing_factor", smoothing_factor); + + smooth_column_collapse_granular(vbx, nbx, + ld.density.array(mfi), + ld.p_nd.array(mfi), + domain, dx, problo, probhi, + smoothing_factor); + } else { amrex::Abort("prob_init_fluid: unknown m_probtype"); @@ -1128,92 +1140,3 @@ void incflo::init_burggraf (Box const& vbx, Box const& /*gbx*/, #endif }); } - -void incflo::column_collapse_granular (Box const& vbx, Box const& nbx, - Array4 const& density, - Array4 const& pressure, - Box const& /*domain*/, - GpuArray const& dx, - GpuArray const& problo, - GpuArray const& probhi) -{ - // Ensure it is set to two_fluid - if (!m_two_fluid) amrex::Abort("probtype 531 involves two fluids"); - - Vector granlen_vec{AMREX_D_DECL(Real(0.),Real(0.),Real(0.))}; - ParmParse pp; - pp.getarr("granular_length",granlen_vec,0,AMREX_SPACEDIM); - GpuArray granLen{AMREX_D_DECL( - granlen_vec[0],granlen_vec[1],granlen_vec[2])}; - granLen[0] += problo[0]; - granLen[1] += problo[1]; - if (granLen[0] > probhi[0]) - amrex::Abort("Granular length along x is larger than setup"); - if (granLen[1] > probhi[1]) - amrex::Abort("Granular length along y is larger than setup"); -#if (AMREX_SPACEDIM==3) - granLen[2] += problo[2]; - if (granLen[2] > probhi[2]) - amrex::Abort("Granular length along z is larger than setup"); -#endif - Real rho_1 = m_ro_0; Real rho_2 = m_ro_0_second; - if (rho_1 > rho_2) - amrex::Abort("Primary fluid must be lighter than second fluid"); - // Density - amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real x = problo[0] + (Real(i)+Real(0.5))*dx[0]; - Real y = problo[1] + (Real(j)+Real(0.5))*dx[1]; -#if (AMREX_SPACEDIM == 3) - Real z = problo[2] + (Real(k)+Real(0.5))*dx[2]; -#endif - if ( (x <= granLen[0]) and (y <= granLen[1]) -#if (AMREX_SPACEDIM == 3) - and (z <= granLen[2]) -#endif - ) { - density(i,j,k) = rho_2; - } else { - density(i,j,k) = rho_1; - } - }); - // Pressure - GpuArray grav{0.,0.,0.}; - if (m_gravity[0]*m_gravity[0] + m_gravity[1]*m_gravity[1]+ - m_gravity[2]*m_gravity[2] > Real(0.)) { - grav[0] = std::abs(m_gravity[0]); - grav[1] = std::abs(m_gravity[1]); - grav[2] = std::abs(m_gravity[2]); - } else if (m_gp0[0]*m_gp0[0] + m_gp0[1]*m_gp0[1]+ - m_gp0[2]*m_gp0[2] > Real(0.)) { - grav[0] = std::abs(m_gp0[0]); - grav[1] = std::abs(m_gp0[1]); - grav[2] = std::abs(m_gp0[2]); - } - else { - amrex::Abort("Body force required in probtype 531"); - } - - amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - // Taking pressure based on the last dimension - Real x = problo[0] + Real(i)*dx[0]; - Real y = problo[1] + Real(j)*dx[1]; - Real p_h = y; -#if (AMREX_SPACEDIM == 3) - Real z = problo[2] + Real(k)*dx[2]; - p_h = z; -#endif - int p_dim = AMREX_SPACEDIM - 1; - if ( (x <= granLen[0]) and (y <= granLen[1]) -#if (AMREX_SPACEDIM == 3) - and (z <= granLen[2]) -#endif - ) { - pressure(i,j,k) = rho_1*grav[p_dim]*(probhi[p_dim]-granLen[p_dim]) - + rho_2*grav[p_dim]*(granLen[p_dim]-p_h); - } else { - pressure(i,j,k) = rho_1*grav[p_dim]*(probhi[p_dim]-p_h); - } - }); -} diff --git a/src/prob/prob_init_fluid_usr.cpp b/src/prob/prob_init_fluid_usr.cpp index bb486826..eeb7d725 100644 --- a/src/prob/prob_init_fluid_usr.cpp +++ b/src/prob/prob_init_fluid_usr.cpp @@ -2,4 +2,193 @@ using namespace amrex; -// This is just a placeholder file +void incflo::column_collapse_granular (Box const& vbx, Box const& nbx, + Array4 const& density, + Array4 const& pressure, + Box const& /*domain*/, + GpuArray const& dx, + GpuArray const& problo, + GpuArray const& probhi) +{ + // Ensure it is set to two_fluid + if (!m_two_fluid) amrex::Abort("probtype 531 involves two fluids"); + + Vector granlen_vec{AMREX_D_DECL(Real(0.),Real(0.),Real(0.))}; + ParmParse pp; + pp.getarr("granular_length",granlen_vec,0,AMREX_SPACEDIM); + GpuArray granLen{AMREX_D_DECL( + granlen_vec[0],granlen_vec[1],granlen_vec[2])}; + granLen[0] += problo[0]; + granLen[1] += problo[1]; + if (granLen[0] > probhi[0]) + amrex::Abort("Granular length along x is larger than setup"); + if (granLen[1] > probhi[1]) + amrex::Abort("Granular length along y is larger than setup"); +#if (AMREX_SPACEDIM==3) + granLen[2] += problo[2]; + if (granLen[2] > probhi[2]) + amrex::Abort("Granular length along z is larger than setup"); +#endif + Real rho_1 = m_ro_0; Real rho_2 = m_ro_0_second; + if (rho_1 > rho_2) + amrex::Abort("Primary fluid must be lighter than second fluid"); + // Density + amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real x = problo[0] + (Real(i)+Real(0.5))*dx[0]; + Real y = problo[1] + (Real(j)+Real(0.5))*dx[1]; +#if (AMREX_SPACEDIM == 3) + Real z = problo[2] + (Real(k)+Real(0.5))*dx[2]; +#endif + if ( (x <= granLen[0]) and (y <= granLen[1]) +#if (AMREX_SPACEDIM == 3) + and (z <= granLen[2]) +#endif + ) { + density(i,j,k) = rho_2; + } else { + density(i,j,k) = rho_1; + } + }); + // Pressure + GpuArray grav{0.,0.,0.}; + if (m_gravity[0]*m_gravity[0] + m_gravity[1]*m_gravity[1]+ + m_gravity[2]*m_gravity[2] > Real(0.)) { + grav[0] = std::abs(m_gravity[0]); + grav[1] = std::abs(m_gravity[1]); + grav[2] = std::abs(m_gravity[2]); + } else if (m_gp0[0]*m_gp0[0] + m_gp0[1]*m_gp0[1]+ + m_gp0[2]*m_gp0[2] > Real(0.)) { + grav[0] = std::abs(m_gp0[0]); + grav[1] = std::abs(m_gp0[1]); + grav[2] = std::abs(m_gp0[2]); + } + else { + amrex::Abort("Body force required in probtype 531"); + } + + amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + // Taking pressure based on the last dimension + Real x = problo[0] + Real(i)*dx[0]; + Real y = problo[1] + Real(j)*dx[1]; + Real p_h = y; +#if (AMREX_SPACEDIM == 3) + Real z = problo[2] + Real(k)*dx[2]; + p_h = z; +#endif + int p_dim = AMREX_SPACEDIM - 1; + if ( (x <= granLen[0]) and (y <= granLen[1]) +#if (AMREX_SPACEDIM == 3) + and (z <= granLen[2]) +#endif + ) { + pressure(i,j,k) = rho_1*grav[p_dim]*(probhi[p_dim]-granLen[p_dim]) + + rho_2*grav[p_dim]*(granLen[p_dim]-p_h); + } else { + pressure(i,j,k) = rho_1*grav[p_dim]*(probhi[p_dim]-p_h); + } + }); +} + +void incflo::smooth_column_collapse_granular (Box const& vbx, Box const& nbx, + Array4 const& density, + Array4 const& pressure, + Box const& /*domain*/, + GpuArray const& dx, + GpuArray const& problo, + GpuArray const& probhi, + Real smoothing_factor) +{ + // Ensure it is set to two_fluid + if (!m_two_fluid) amrex::Abort("probtype 531 involves two fluids"); + + Vector granlen_vec{AMREX_D_DECL(Real(0.),Real(0.),Real(0.))}; + ParmParse pp; + pp.getarr("granular_length",granlen_vec,0,AMREX_SPACEDIM); + GpuArray granLen{AMREX_D_DECL( + granlen_vec[0],granlen_vec[1],granlen_vec[2])}; + granLen[0] += problo[0]; + granLen[1] += problo[1]; + if (granLen[0] > probhi[0]) + amrex::Abort("Granular length along x is larger than setup"); + if (granLen[1] > probhi[1]) + amrex::Abort("Granular length along y is larger than setup"); +#if (AMREX_SPACEDIM==3) + granLen[2] += problo[2]; + if (granLen[2] > probhi[2]) + amrex::Abort("Granular length along z is larger than setup"); +#endif + Real rho_1 = m_ro_0; Real rho_2 = m_ro_0_second; + if (rho_1 > rho_2) + amrex::Abort("Primary fluid must be lighter than second fluid"); + // Density + amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real x = problo[0] + (Real(i)+Real(0.5))*dx[0]; + Real y = problo[1] + (Real(j)+Real(0.5))*dx[1]; +#if (AMREX_SPACEDIM == 3) + Real z = problo[2] + (Real(k)+Real(0.5))*dx[2]; +#endif + +#if (AMREX_SPACEDIM == 2) + density(i,j,k) = std::tanh((granLen[1]-y)/(smoothing_factor*dx[1])) + + Real(1.0); +#else + density(i,j,k) = std::tanh((granLen[2]-z)/(smoothing_factor*dx[2])) + + Real(1.0); + +#endif + density(i,j,k) *= (std::tanh((granLen[0]-x)/(smoothing_factor*dx[0])) + + Real(1.0)); + density(i,j,k) *= Real(0.25)*(rho_2-rho_1); + density(i,j,k) += rho_1; + }); + // Pressure + GpuArray grav{0.,0.,0.}; + if (m_gravity[0]*m_gravity[0] + m_gravity[1]*m_gravity[1]+ + m_gravity[2]*m_gravity[2] > Real(0.)) { + grav[0] = std::abs(m_gravity[0]); + grav[1] = std::abs(m_gravity[1]); + grav[2] = std::abs(m_gravity[2]); + } else if (m_gp0[0]*m_gp0[0] + m_gp0[1]*m_gp0[1]+ + m_gp0[2]*m_gp0[2] > Real(0.)) { + grav[0] = std::abs(m_gp0[0]); + grav[1] = std::abs(m_gp0[1]); + grav[2] = std::abs(m_gp0[2]); + } + else { + amrex::Abort("Body force required in probtype 531"); + } + + amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + // Taking pressure based on the last dimension + Real x = problo[0] + Real(i)*dx[0]; + Real y = problo[1] + Real(j)*dx[1]; + Real p_grav = grav[1]; + Real p_h = y; + Real p_h_max = probhi[1]; + Real p_h_gran = granLen[1]; + // This is a multiplication factor along gravity + Real p_fctr = smoothing_factor*dx[1]; +#if (AMREX_SPACEDIM == 3) + Real z = problo[2] + Real(k)*dx[2]; + p_grav = grav[2]; + p_h = z; + p_h_max = probhi[2]; + p_h_gran = granLen[2]; + p_fctr = smoothing_factor*dx[2]; +#endif + Real fctr_x = std::tanh((granLen[0]-x)/(smoothing_factor*dx[0])) + + Real(1.0); + fctr_x *= ((rho_2-rho_1)*p_grav*Real(0.25)); + + pressure(i,j,k) = std::cosh((p_h_gran-p_h)/p_fctr); + pressure(i,j,k) /= std::cosh((p_h_gran-p_h_max)/p_fctr); + pressure(i,j,k) = std::log(pressure(i,j,k)); + pressure(i,j,k) *= (fctr_x*p_fctr); + pressure(i,j,k) += (fctr_x*(p_h_max-p_h)); + pressure(i,j,k) += (rho_1*p_grav*(p_h_max-p_h)); + }); +} diff --git a/src/rheology/incflo_rheology.cpp b/src/rheology/incflo_rheology.cpp index 2794e295..6591e594 100644 --- a/src/rheology/incflo_rheology.cpp +++ b/src/rheology/incflo_rheology.cpp @@ -314,16 +314,17 @@ void incflo::compute_nodal_viscosity_at_level (int /*lev*/, } // Obtain concentration of the second fluid, based on nodal density MultiFab rho_nodal(vel_eta->boxArray(),vel_eta->DistributionMap(),1,nghost); - MultiFab conc_second(vel_eta->boxArray(),vel_eta->DistributionMap(),1,nghost); + // Nodal second fluid concentration MultiFab + MultiFab conc_second_nd(vel_eta->boxArray(),vel_eta->DistributionMap(),1,nghost); #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - for (MFIter mfi(conc_second,TilingIfNotGPU()); mfi.isValid(); ++mfi) + for (MFIter mfi(conc_second_nd,TilingIfNotGPU()); mfi.isValid(); ++mfi) { Box const& bx = mfi.growntilebox(nghost); Array4 const& rho_arr = rho->const_array(mfi); Array4 const& rho_nodal_arr = rho_nodal.array(mfi); - Array4 const& conc_second_arr = conc_second.array(mfi); + Array4 const& conc_second_arr = conc_second_nd.array(mfi); const Real rho_first = m_ro_0; const Real rho_second = m_ro_0_second; amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -428,7 +429,7 @@ void incflo::compute_nodal_viscosity_at_level (int /*lev*/, }); } // Copy concentration - MultiFab::Copy(inertial_num,conc_second,0,1,1,nghost); + MultiFab::Copy(inertial_num,conc_second_nd,0,1,1,nghost); #ifdef USE_AMREX_MPMD if (m_fluid_model_second == FluidModel::DataDrivenMPMD) { @@ -530,10 +531,10 @@ void incflo::compute_nodal_viscosity_at_level (int /*lev*/, #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - for (MFIter mfi(conc_second,TilingIfNotGPU()); mfi.isValid(); ++mfi) + for (MFIter mfi(conc_second_nd,TilingIfNotGPU()); mfi.isValid(); ++mfi) { Box const& bx = mfi.growntilebox(nghost); - Array4 const& conc_second_arr = conc_second.array(mfi); + Array4 const& conc_second_arr = conc_second_nd.array(mfi); Array4 const& eta_arr_second = vel_eta_second.const_array(mfi); Array4 const& eta_arr = vel_eta->array(mfi); const Real min_conc_scnd = m_min_conc_second;