Skip to content

Commit

Permalink
Fix indexing error in computing the advective term for
Browse files Browse the repository at this point in the history 
non-conservative tracers. In practice, this makes little difference
because it's multiplied by div(Umac) which the MAC projection just
enforced to be numerically zero to within its tolerance.

Also make the density field non-constant and add a second,
non-conservative tracer in the tracer_adv_diff_cn regression test.
  • Loading branch information
@cgilet
cgilet committed Jun 7, 2024
1 parent 0a40f11 commit b1464aa
Show file tree
Hide file tree
Showing 5 changed files with 55 additions and 23 deletions.
6 changes: 3 additions & 3 deletions src/convection/incflo_compute_advection_term.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -697,9 +697,9 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
// If convective, we define u dot grad trac = div (u trac) - trac div(u)
HydroUtils::ComputeConvectiveTerm(bx, m_ntrac, mfi,
tracer[lev]->array(mfi,0),
AMREX_D_DECL(face_x[lev].array(mfi),
face_y[lev].array(mfi),
face_z[lev].array(mfi)),
AMREX_D_DECL(face_x[lev].array(mfi,flux_comp),
face_y[lev].array(mfi,flux_comp),
face_z[lev].array(mfi,flux_comp)),
divu[lev].array(mfi),
update_arr,
get_tracer_iconserv_device_ptr(),
Expand Down
4 changes: 2 additions & 2 deletions src/prob/incflo_prob_I.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,14 +84,14 @@
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& problo,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& probhi) const;

static void init_periodic_tracer (amrex::Box const& vbx, amrex::Box const& gbx,
void init_periodic_tracer (amrex::Box const& vbx, amrex::Box const& gbx,
amrex::Array4<amrex::Real> const& vel,
amrex::Array4<amrex::Real> const& density,
amrex::Array4<amrex::Real> const& tracer,
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::GpuArray<amrex::Real, AMREX_SPACEDIM> const& probhi) const;

void init_double_shear_layer (amrex::Box const& vbx, amrex::Box const& gbx,
amrex::Array4<amrex::Real> const& vel,
Expand Down
55 changes: 41 additions & 14 deletions src/prob/prob_init_fluid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,7 +100,7 @@ void incflo::prob_init_fluid (int lev)
ld.tracer.array(mfi),
domain, dx, problo, probhi);
}
else if (12 == m_probtype)
else if (12 == m_probtype || 122 == m_probtype)
{
init_periodic_tracer(vbx, gbx,
ld.velocity.array(mfi),
Expand Down Expand Up @@ -809,30 +809,57 @@ void incflo::init_boussinesq_bubble (Box const& vbx, Box const& /*gbx*/,

void incflo::init_periodic_tracer (Box const& vbx, Box const& /*gbx*/,
Array4<Real> const& vel,
Array4<Real> const& /*density*/,
Array4<Real> const& density,
Array4<Real> const& tracer,
Box const& /*domain*/,
GpuArray<Real, AMREX_SPACEDIM> const& dx,
GpuArray<Real, AMREX_SPACEDIM> const& problo,
GpuArray<Real, AMREX_SPACEDIM> const& probhi)
GpuArray<Real, AMREX_SPACEDIM> const& probhi) const
{
Real L = probhi[0]-problo[0];
Real C = Real(2.0)*Real(3.1415926535897932) / L;
amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept

if (m_probtype == 12)
{
constexpr Real A = Real(1.0);
Real x = Real(i+0.5)*dx[0];
Real y = Real(j+0.5)*dx[1];
amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
constexpr Real A = Real(1.0);
Real x = Real(i+0.5)*dx[0];
Real y = Real(j+0.5)*dx[1];
#if (AMREX_SPACEDIM == 3)
Real z = Real(k+0.5)*dx[2];
Real z = Real(k+0.5)*dx[2];
#else
Real z = 0.0_rt;
Real z = 0.0_rt;
#endif
AMREX_D_TERM(vel(i,j,k,0) = Real(1.0);,
vel(i,j,k,1) = Real(0.1)*(std::sin(C*(x+z) - Real(0.00042)) + Real(1.0)) * std::exp(y);,
vel(i,j,k,2) = Real(0.1)*(std::sin(C*(x+y) - Real(0.00042)) + Real(1.0)) * std::exp(z););
tracer(i,j,k) = A *(std::sin(C*(y+z) - Real(0.00042)) + Real(1.0)) * std::exp(x);
});
AMREX_D_TERM(vel(i,j,k,0) = Real(1.0);,
vel(i,j,k,1) = Real(0.1)*(std::sin(C*(x+z) - Real(0.00042)) + Real(1.0)) * std::exp(y);,
vel(i,j,k,2) = Real(0.1)*(std::sin(C*(x+y) - Real(0.00042)) + Real(1.0)) * std::exp(z););
tracer(i,j,k,0) = A *(std::sin(C*(y+z) - Real(0.00042)) + Real(1.0)) * std::exp(x);
});
}
else if (m_probtype == 122)
{
Real B = Real(0.1);
amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
constexpr Real A = Real(1.0);
Real x = Real(i+0.5)*dx[0];
Real y = Real(j+0.5)*dx[1];
#if (AMREX_SPACEDIM == 3)
Real z = Real(k+0.5)*dx[2];
#else
Real z = 0.0_rt;
#endif
AMREX_D_TERM(vel(i,j,k,0) = Real(1.0);,
vel(i,j,k,1) = Real(0.1)*(std::sin(C*(x+z) - Real(0.00042)) + Real(1.0)) * std::exp(y);,
vel(i,j,k,2) = Real(0.1)*(std::sin(C*(x+y) - Real(0.00042)) + Real(1.0)) * std::exp(z););
density(i,j,k) = A + B*(x + y + z);
tracer(i,j,k,0) = A *(std::sin(C*(y+z) - Real(0.00042)) + Real(1.0)) * std::exp(x);
tracer(i,j,k,1) = A *(std::sin(C*(y+z) - Real(0.00042)) + Real(1.0)) * std::exp(x);
});
} else {
Abort("Unknow periodic tracer probtype");
}
}

void incflo::init_double_shear_layer (Box const& vbx, Box const& /*gbx*/,
Expand Down
6 changes: 4 additions & 2 deletions test_2d/benchmark.tracer_adv_diff_cn
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,9 @@ incflo.mu = 0.001 # Dynamic viscosity coefficient

incflo.constant_density = false #
incflo.advect_tracer = true #
incflo.mu_s = 0.1
incflo.ntrac = 2
incflo.mu_s = 0.1 0.1
incflo.trac_is_conservative = 1 0

incflo.diffusion_type = 1 # Crank-Nicolson (instead of Implicit)

Expand Down Expand Up @@ -71,7 +73,7 @@ incflo.write_eb_surface = false
#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# INITIAL CONDITIONS #
#.......................................#
incflo.probtype = 12
incflo.probtype = 122
incflo.test_tracer_conservation = true

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
Expand Down
7 changes: 5 additions & 2 deletions test_3d/benchmark.tracer_adv_diff_cn
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,10 @@ incflo.mu = 0.001 # Dynamic viscosity coefficient

incflo.constant_density = false #
incflo.advect_tracer = true #
incflo.mu_s = 0.1
incflo.ntrac = 2
incflo.mu_s = 0.1 0.1
incflo.trac_is_conservative = 1 0


incflo.diffusion_type = 1 # Crank-Nicolson (instead of Implicit)

Expand Down Expand Up @@ -75,7 +78,7 @@ incflo.write_eb_surface = false
#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# INITIAL CONDITIONS #
#.......................................#
incflo.probtype = 12
incflo.probtype = 122
incflo.test_tracer_conservation = true

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
Expand Down

0 comments on commit b1464aa

Please sign in to comment.