Skip to content

Commit

Permalink
A diffuse/smooth version of column_collapse_granular added
Browse files Browse the repository at this point in the history
  • Loading branch information
@siddanib
siddanib committed Jul 20, 2024
1 parent d4c2cbb commit 99c166d
Show file tree
Hide file tree
Showing 5 changed files with 226 additions and 104 deletions.
8 changes: 0 additions & 8 deletions src/prob/incflo_prob_I.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,12 +164,4 @@
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& dx,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& problo,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& probhi) const;

void column_collapse_granular (amrex::Box const& vbx, amrex::Box const& nbx,
amrex::Array4<amrex::Real> const& density,
amrex::Array4<amrex::Real> const& pressure,
amrex::Box const& domain,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& dx,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& problo,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& probhi);
#endif
17 changes: 17 additions & 0 deletions src/prob/incflo_prob_usr_I.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -131,4 +131,21 @@
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& dx,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& problo,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& probhi);

void column_collapse_granular (amrex::Box const& vbx, amrex::Box const& nbx,
amrex::Array4<amrex::Real> const& density,
amrex::Array4<amrex::Real> const& pressure,
amrex::Box const& domain,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& dx,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& problo,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& probhi);

void smooth_column_collapse_granular (amrex::Box const& vbx, amrex::Box const& nbx,
amrex::Array4<amrex::Real> const& density,
amrex::Array4<amrex::Real> const& pressure,
amrex::Box const& domain,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& dx,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& problo,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& probhi,
amrex::Real smoothing_factor);
#endif
101 changes: 12 additions & 89 deletions src/prob/prob_init_fluid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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");
Expand Down Expand Up @@ -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<Real> const& density,
Array4<Real> const& pressure,
Box const& /*domain*/,
GpuArray<Real, AMREX_SPACEDIM> const& dx,
GpuArray<Real, AMREX_SPACEDIM> const& problo,
GpuArray<Real, AMREX_SPACEDIM> const& probhi)
{
// Ensure it is set to two_fluid
if (!m_two_fluid) amrex::Abort("probtype 531 involves two fluids");

Vector<Real> granlen_vec{AMREX_D_DECL(Real(0.),Real(0.),Real(0.))};
ParmParse pp;
pp.getarr("granular_length",granlen_vec,0,AMREX_SPACEDIM);
GpuArray<Real,AMREX_SPACEDIM> 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<Real,3> 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);
}
});
}
191 changes: 190 additions & 1 deletion src/prob/prob_init_fluid_usr.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<Real> const& density,
Array4<Real> const& pressure,
Box const& /*domain*/,
GpuArray<Real, AMREX_SPACEDIM> const& dx,
GpuArray<Real, AMREX_SPACEDIM> const& problo,
GpuArray<Real, AMREX_SPACEDIM> const& probhi)
{
// Ensure it is set to two_fluid
if (!m_two_fluid) amrex::Abort("probtype 531 involves two fluids");

Vector<Real> granlen_vec{AMREX_D_DECL(Real(0.),Real(0.),Real(0.))};
ParmParse pp;
pp.getarr("granular_length",granlen_vec,0,AMREX_SPACEDIM);
GpuArray<Real,AMREX_SPACEDIM> 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<Real,3> 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<Real> const& density,
Array4<Real> const& pressure,
Box const& /*domain*/,
GpuArray<Real, AMREX_SPACEDIM> const& dx,
GpuArray<Real, AMREX_SPACEDIM> const& problo,
GpuArray<Real, AMREX_SPACEDIM> const& probhi,
Real smoothing_factor)
{
// Ensure it is set to two_fluid
if (!m_two_fluid) amrex::Abort("probtype 531 involves two fluids");

Vector<Real> granlen_vec{AMREX_D_DECL(Real(0.),Real(0.),Real(0.))};
ParmParse pp;
pp.getarr("granular_length",granlen_vec,0,AMREX_SPACEDIM);
GpuArray<Real,AMREX_SPACEDIM> 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");
if (granLen[1] != probhi[1])
amrex::Abort(
"For 3D case, domain length and granular_length along y need to be same");
#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];
#if (AMREX_SPACEDIM == 2)
Real y = problo[1] + (Real(j)+Real(0.5))*dx[1];
density(i,j,k) = std::tanh((granLen[1]-y)/(smoothing_factor*dx[1]))
+ Real(1.0);
#else
Real z = problo[2] + (Real(k)+Real(0.5))*dx[2];
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<Real,3> 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));
});
}
13 changes: 7 additions & 6 deletions src/rheology/incflo_rheology.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<Real const> const& rho_arr = rho->const_array(mfi);
Array4<Real> const& rho_nodal_arr = rho_nodal.array(mfi);
Array4<Real> const& conc_second_arr = conc_second.array(mfi);
Array4<Real> 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
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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<Real const> const& conc_second_arr = conc_second.array(mfi);
Array4<Real const> const& conc_second_arr = conc_second_nd.array(mfi);
Array4<Real const> const& eta_arr_second = vel_eta_second.const_array(mfi);
Array4<Real> const& eta_arr = vel_eta->array(mfi);
const Real min_conc_scnd = m_min_conc_second;
Expand Down

0 comments on commit 99c166d

Please sign in to comment.