diff --git a/src/derive/incflo_derive.cpp b/src/derive/incflo_derive.cpp index 2311125f..0e574821 100644 --- a/src/derive/incflo_derive.cpp +++ b/src/derive/incflo_derive.cpp @@ -178,9 +178,18 @@ void incflo::compute_nodal_hydrostatic_pressure_at_level (int lev, Box const& bx = mfi.tilebox(); Array4 const& rho_arr = rho_cc->const_array(mfi); Array4 const& rho_nodal_arr = rho_nodal.array(mfi); + const int prob_534 = (m_probtype == 534) ? 1 : 0; + const Real rho_1 = m_ro_0; 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); + // In inclined_plane_granular, the free surface needs to have + // hydrostatic pressure of zero + if (prob_534) { + if (rho_nodal_arr(i,j,k) == rho_1) { + rho_nodal_arr(i,j,k) -= rho_1; + } + } }); } diff --git a/src/prob/incflo_prob_usr_I.H b/src/prob/incflo_prob_usr_I.H index 578c0825..ca7ef205 100644 --- a/src/prob/incflo_prob_usr_I.H +++ b/src/prob/incflo_prob_usr_I.H @@ -159,4 +159,33 @@ amrex::GpuArray const& problo, amrex::GpuArray const& probhi) const; #endif + void single_layer_inclined_plane_granular (amrex::Box const& vbx, amrex::Box const& nbx, + amrex::Array4 const& density, + amrex::Array4 const& tracer, + amrex::Array4 const& pressure, + amrex::Array4 const& velocity, + amrex::Box const& domain, + amrex::GpuArray const& dx, + amrex::GpuArray const& problo, + amrex::GpuArray const& probhi) const; + + void double_layer_inclined_plane_granular (amrex::Box const& vbx, amrex::Box const& nbx, + amrex::Array4 const& density, + amrex::Array4 const& tracer, + amrex::Array4 const& pressure, + amrex::Array4 const& velocity, + amrex::Box const& domain, + amrex::GpuArray const& dx, + amrex::GpuArray const& problo, + amrex::GpuArray const& probhi) const; + + void smooth_double_layer_inclined_plane_granular (amrex::Box const& vbx, amrex::Box const& nbx, + amrex::Array4 const& density, + amrex::Array4 const& tracer, + amrex::Array4 const& pressure, + amrex::Array4 const& velocity, + amrex::Box const& domain, + amrex::GpuArray const& dx, + amrex::GpuArray const& problo, + amrex::GpuArray const& probhi) const; #endif diff --git a/src/prob/prob_init_fluid.cpp b/src/prob/prob_init_fluid.cpp index 76846a44..6c3df575 100644 --- a/src/prob/prob_init_fluid.cpp +++ b/src/prob/prob_init_fluid.cpp @@ -218,6 +218,41 @@ void incflo::prob_init_fluid (int lev) domain, dx, problo, probhi); } #endif + // 534 corresponds to inclined_plane_granular, + // but the initialization is the same as column_collapse_granular + // 534 ENFORCES ZERO hydrostatic pressure for free surface + else if (534 == m_probtype) + { + column_collapse_granular(vbx, nbx, + ld.density.array(mfi), + ld.tracer.array(mfi), + ld.p_nd.array(mfi), + domain, dx, problo, probhi); + } + else if (535 == m_probtype) { + single_layer_inclined_plane_granular(vbx, nbx, + ld.density.array(mfi), + ld.tracer.array(mfi), + ld.p_nd.array(mfi), + ld.velocity.array(mfi), + domain, dx, problo, probhi); + } + else if (536 == m_probtype) { + double_layer_inclined_plane_granular(vbx, nbx, + ld.density.array(mfi), + ld.tracer.array(mfi), + ld.p_nd.array(mfi), + ld.velocity.array(mfi), + domain, dx, problo, probhi); + } + else if (537 == m_probtype) { + smooth_double_layer_inclined_plane_granular(vbx, nbx, + ld.density.array(mfi), + ld.tracer.array(mfi), + ld.p_nd.array(mfi), + ld.velocity.array(mfi), + domain, dx, problo, probhi); + } else { amrex::Abort("prob_init_fluid: unknown m_probtype"); diff --git a/src/prob/prob_init_fluid_usr.cpp b/src/prob/prob_init_fluid_usr.cpp index ac3356df..42cbd86b 100644 --- a/src/prob/prob_init_fluid_usr.cpp +++ b/src/prob/prob_init_fluid_usr.cpp @@ -281,3 +281,282 @@ void incflo::cylinder_collapse_granular (Box const& vbx, Box const& nbx, amrex::Abort("Include non-zero initial_iterations as pressure is not set"); } #endif + +void incflo::single_layer_inclined_plane_granular (Box const& vbx, Box const& nbx, + Array4 const& density, + Array4 const& tracer, + Array4 const& pressure, + Array4 const& velocity, + Box const& /*domain*/, + GpuArray const& dx, + GpuArray const& problo, + GpuArray const& probhi) const +{ + amrex::ignore_unused(nbx); + amrex::ignore_unused>(pressure); + // Ensure it is set to two_fluid + if (!m_two_fluid) amrex::Abort("probtype 535 leverages two_fluid infrastructure"); + Real grav_mag = m_gravity[0]*m_gravity[0] + m_gravity[1]*m_gravity[1]+ + m_gravity[2]*m_gravity[2]; + grav_mag = std::sqrt(grav_mag); + if (grav_mag == Real(0.)) { + amrex::Abort("Gravity is required in probtype 535"); + } + Real zeta_0; // Inclination angle used in initialization + Real I_zeta_0; // Inertial number corresponding to initialization + ParmParse pp; + pp.get("zeta_0", zeta_0); + pp.get("I_zeta_0",I_zeta_0); + amrex::Print() << "zeta_0: "< const& density, + Array4 const& tracer, + Array4 const& pressure, + Array4 const& velocity, + Box const& /*domain*/, + GpuArray const& dx, + GpuArray const& problo, + GpuArray const& probhi) const +{ + amrex::ignore_unused(nbx); + amrex::ignore_unused>(pressure); + // Ensure it is set to two_fluid + if (!m_two_fluid) amrex::Abort("probtype 536 leverages two_fluid infrastructure"); + Real grav_mag = m_gravity[0]*m_gravity[0] + m_gravity[1]*m_gravity[1]+ + m_gravity[2]*m_gravity[2]; + grav_mag = std::sqrt(grav_mag); + if (grav_mag == Real(0.)) { + amrex::Abort("Gravity is required in probtype 536"); + } + Real zeta_0; // Inclination angle used in initialization + Real I_zeta_0; // Inertial number corresponding to initialization + ParmParse pp; + pp.get("zeta_0", zeta_0); + pp.get("I_zeta_0",I_zeta_0); + amrex::Print() << "zeta_0: "< granlen_vec{AMREX_D_DECL(Real(0.),Real(0.),Real(0.))}; + 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 needs to be equal to setup"); +#if (AMREX_SPACEDIM == 2) + if (granLen[1] >= probhi[1]) + amrex::Abort("Granular length along y needs to be smaller than setup"); +#else + if (granLen[1] != probhi[1]) + amrex::Abort("Granular length along y needs to be equal to setup"); + + granLen[2] += problo[2]; + if (granLen[2] >= probhi[2]) + amrex::Abort("Granular length along z needs to be smaller than setup"); +#endif + Real rho_1 = m_ro_0; Real rho_2 = m_ro_0_second; + Real rho_p = m_ro_grain_second; // Granular particle density + Real diam_p = m_diam_second; // Granular particle mean diameter + Real H_c = granLen[1]; // Column Height + Real H_d = probhi[1]; // Domain Height +#if (AMREX_SPACEDIM == 3) + H_c = granLen[2]; + H_d = probhi[2]; +#endif + // Create a variable to store max Bagnold velocity + Real u_bag_max = std::sqrt(grav_mag*std::cos(zeta_0)*rho_2/rho_p); + u_bag_max *= (Real(4.0)*I_zeta_0/(Real(3.0)*diam_p)); + u_bag_max *= std::pow(H_c,Real(1.5)); + + 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; + tracer(i,j,k,0) = Real(1.0); + // Bagnold profile based on zeta_0 + velocity(i,j,k,0) = std::sqrt(grav_mag*std::cos(zeta_0)*rho_2/rho_p); + velocity(i,j,k,0) *= (Real(4.0)*I_zeta_0/(Real(3.0)*diam_p)); +#if (AMREX_SPACEDIM == 2) + velocity(i,j,k,0) *= (std::pow(H_c,Real(1.5)) - std::pow(H_c-y,Real(1.5))); +#else + velocity(i,j,k,0) *= (std::pow(H_c,Real(1.5)) - std::pow(H_c-z,Real(1.5))); +#endif + } else { + density(i,j,k) = rho_1; + tracer(i,j,k,0) = Real(0.0); + // Linear profile for lighter fluid +#if (AMREX_SPACEDIM == 2) + velocity(i,j,k,0) = (H_d - y)/(H_d - H_c); +#else + velocity(i,j,k,0) = (H_d - z)/(H_d - H_c); +#endif + velocity(i,j,k,0) *= u_bag_max; + } + + velocity(i,j,k,1) = Real(0.0); +#if (AMREX_SPACEDIM == 3) + velocity(i,j,k,2) = Real(0.0); +#endif + }); + + if (m_initial_iterations == 0 ) + amrex::Abort("Include non-zero initial_iterations as pressure is not set"); +} + +void incflo::smooth_double_layer_inclined_plane_granular (Box const& vbx, Box const& nbx, + Array4 const& density, + Array4 const& tracer, + Array4 const& pressure, + Array4 const& velocity, + Box const& /*domain*/, + GpuArray const& dx, + GpuArray const& problo, + GpuArray const& probhi) const +{ + amrex::ignore_unused(nbx); + amrex::ignore_unused>(pressure); + // Ensure it is set to two_fluid + if (!m_two_fluid) amrex::Abort("probtype 537 requires two_fluid"); + Real grav_mag = m_gravity[0]*m_gravity[0] + m_gravity[1]*m_gravity[1]+ + m_gravity[2]*m_gravity[2]; + grav_mag = std::sqrt(grav_mag); + if (grav_mag == Real(0.)) { + amrex::Abort("Gravity is required in probtype 537"); + } + Real zeta_0; // Inclination angle used in initialization + Real I_zeta_0; // Inertial number corresponding to initialization + ParmParse pp; + pp.get("zeta_0", zeta_0); + pp.get("I_zeta_0",I_zeta_0); + amrex::Print() << "zeta_0: "< granlen_vec{AMREX_D_DECL(Real(0.),Real(0.),Real(0.))}; + 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 needs to be equal to setup"); +#if (AMREX_SPACEDIM == 2) + if (granLen[1] >= probhi[1]) + amrex::Abort("Granular length along y needs to be smaller than setup"); +#else + if (granLen[1] != probhi[1]) + amrex::Abort("Granular length along y needs to be equal to setup"); + + granLen[2] += problo[2]; + if (granLen[2] >= probhi[2]) + amrex::Abort("Granular length along z needs to be smaller than setup"); +#endif + Real rho_1 = m_ro_0; Real rho_2 = m_ro_0_second; + Real rho_p = m_ro_grain_second; // Granular particle density + Real diam_p = m_diam_second; // Granular particle mean diameter + Real H_c = granLen[1]; // Column Height + Real H_d = probhi[1]; // Domain Height +#if (AMREX_SPACEDIM == 3) + H_c = granLen[2]; + H_d = probhi[2]; +#endif + // Create a variable to store max Bagnold velocity + Real u_bag_max = std::sqrt(grav_mag*std::cos(zeta_0)*rho_2/rho_p); + u_bag_max *= (Real(4.0)*I_zeta_0/(Real(3.0)*diam_p)); + u_bag_max *= std::pow(H_c,Real(1.5)); + + amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + // Tracer is 1 in heavier fluid + 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 == 2) + tracer(i,j,k,0) = std::tanh((granLen[1]-y)/(smoothing_factor*dx[1])) + + Real(1.0); +#else + Real z = problo[2] + (Real(k)+Real(0.5))*dx[2]; + tracer(i,j,k,0) = std::tanh((granLen[2]-z)/(smoothing_factor*dx[2])) + + Real(1.0); + +#endif + tracer(i,j,k,0) *= Real(0.5); + density(i,j,k) = tracer(i,j,k,0)*(rho_2-rho_1); + density(i,j,k) += rho_1; + + if ( (x <= granLen[0]) and (y <= granLen[1]) +#if (AMREX_SPACEDIM == 3) + and (z <= granLen[2]) +#endif + ) { + // Bagnold profile based on zeta_0 + velocity(i,j,k,0) = std::sqrt(grav_mag*std::cos(zeta_0)*rho_2/rho_p); + velocity(i,j,k,0) *= (Real(4.0)*I_zeta_0/(Real(3.0)*diam_p)); +#if (AMREX_SPACEDIM == 2) + velocity(i,j,k,0) *= (std::pow(H_c,Real(1.5)) - std::pow(H_c-y,Real(1.5))); +#else + velocity(i,j,k,0) *= (std::pow(H_c,Real(1.5)) - std::pow(H_c-z,Real(1.5))); +#endif + } else { + // Linear profile for lighter fluid +#if (AMREX_SPACEDIM == 2) + velocity(i,j,k,0) = (H_d - y)/(H_d - H_c); +#else + velocity(i,j,k,0) = (H_d - z)/(H_d - H_c); +#endif + velocity(i,j,k,0) *= u_bag_max; + } + + velocity(i,j,k,1) = Real(0.0); +#if (AMREX_SPACEDIM == 3) + velocity(i,j,k,2) = Real(0.0); +#endif + }); + + if (m_initial_iterations == 0 ) + amrex::Abort("Include non-zero initial_iterations as pressure is not set"); +} diff --git a/src/setup/set_background_pressure.cpp b/src/setup/set_background_pressure.cpp index e80218c5..a98561bd 100644 --- a/src/setup/set_background_pressure.cpp +++ b/src/setup/set_background_pressure.cpp @@ -15,10 +15,11 @@ void incflo::set_background_pressure () GpuArray problen{AMREX_D_DECL(probhi[0]-problo[0], probhi[1]-problo[1], probhi[2]-problo[2])}; + int delp_dir; // There are 3 exclusive sources for background pressure gradient. - // (1) incflo.delp in inputs - int delp_dir = -1; for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { + delp_dir = -1; + // (1) incflo.delp in inputs if (std::abs(m_delp[dir]) > std::numeric_limits::epsilon()) { if (delp_dir == -1) { delp_dir = dir; @@ -27,9 +28,7 @@ void incflo::set_background_pressure () amrex::Abort("set_background_pressure: how did this happen?"); } } - } - // (2) pressure inflow and pressure outflow - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { + // (2) pressure inflow and pressure outflow if ((m_bc_type[Orientation(dir,Orientation::low)] == BC::pressure_inflow and m_bc_type[Orientation(dir,Orientation::high)] == BC::pressure_outflow) or (m_bc_type[Orientation(dir,Orientation::high)] == BC::pressure_inflow and @@ -43,9 +42,7 @@ void incflo::set_background_pressure () amrex::Abort("set_background_pressure: how did this happen?"); } } - } - // (3) gravity - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { + // (3) gravity Real dpdx = m_gravity[dir] * m_ro_0; if (std::abs(dpdx) > std::numeric_limits::epsilon()) { if (delp_dir == -1) {