diff --git a/src/convection/incflo_compute_advection_term.cpp b/src/convection/incflo_compute_advection_term.cpp index 16d2df1b..8f7fbcbf 100644 --- a/src/convection/incflo_compute_advection_term.cpp +++ b/src/convection/incflo_compute_advection_term.cpp @@ -697,9 +697,9 @@ incflo::compute_convective_term (Vector 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(), diff --git a/src/prob/incflo_prob_I.H b/src/prob/incflo_prob_I.H index 080c1736..7f035e76 100644 --- a/src/prob/incflo_prob_I.H +++ b/src/prob/incflo_prob_I.H @@ -84,14 +84,14 @@ amrex::GpuArray const& problo, amrex::GpuArray 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 const& vel, amrex::Array4 const& density, amrex::Array4 const& tracer, amrex::Box const& domain, amrex::GpuArray const& dx, amrex::GpuArray const& problo, - amrex::GpuArray const& probhi); + amrex::GpuArray const& probhi) const; void init_double_shear_layer (amrex::Box const& vbx, amrex::Box const& gbx, amrex::Array4 const& vel, diff --git a/src/prob/prob_init_fluid.cpp b/src/prob/prob_init_fluid.cpp index c6f14299..93cecb35 100644 --- a/src/prob/prob_init_fluid.cpp +++ b/src/prob/prob_init_fluid.cpp @@ -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), @@ -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 const& vel, - Array4 const& /*density*/, + Array4 const& density, Array4 const& tracer, Box const& /*domain*/, GpuArray const& dx, GpuArray const& problo, - GpuArray const& probhi) + GpuArray 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*/, diff --git a/test_2d/benchmark.tracer_adv_diff_cn b/test_2d/benchmark.tracer_adv_diff_cn index 3ee73021..29f30b09 100644 --- a/test_2d/benchmark.tracer_adv_diff_cn +++ b/test_2d/benchmark.tracer_adv_diff_cn @@ -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) @@ -71,7 +73,7 @@ incflo.write_eb_surface = false #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# # INITIAL CONDITIONS # #.......................................# -incflo.probtype = 12 +incflo.probtype = 122 incflo.test_tracer_conservation = true #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# diff --git a/test_3d/benchmark.tracer_adv_diff_cn b/test_3d/benchmark.tracer_adv_diff_cn index c6bd43df..d64727ee 100644 --- a/test_3d/benchmark.tracer_adv_diff_cn +++ b/test_3d/benchmark.tracer_adv_diff_cn @@ -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) @@ -75,7 +78,7 @@ incflo.write_eb_surface = false #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# # INITIAL CONDITIONS # #.......................................# -incflo.probtype = 12 +incflo.probtype = 122 incflo.test_tracer_conservation = true #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#