diff --git a/Exec/DevTests/EB_Test/amrvis.defaults b/Exec/DevTests/EB_Test/amrvis.defaults new file mode 100644 index 000000000..71eb95eab --- /dev/null +++ b/Exec/DevTests/EB_Test/amrvis.defaults @@ -0,0 +1,10 @@ +palette /home/almgren/bin/Palette +initialderived volfrac +initialscale 8 +expansion 1 +numberformat %8.6e +maxpixmapsize 1000000 +reservesystemcolors 35 +showboxes TRUE +windowheight 650 +windowwidth 1100 diff --git a/Exec/DevTests/EB_Test/inputs b/Exec/DevTests/EB_Test/inputs index b811a2675..61cccb18b 100644 --- a/Exec/DevTests/EB_Test/inputs +++ b/Exec/DevTests/EB_Test/inputs @@ -2,6 +2,8 @@ max_step = 100 max_step = 0 +amr.max_grid_size = 256 256 256 + eb2.geometry = terrain amrex.fpe_trap_invalid = 1 @@ -9,16 +11,27 @@ amrex.fpe_trap_invalid = 1 fabarray.mfiter_tile_size = 1024 1024 1024 # PROBLEM SIZE & GEOMETRY -geometry.prob_extent = 16. 16. 16. -amr.n_cell = 16 16 16 +geometry.prob_lo = 0. 0. 0. +geometry.prob_hi = 10. 1. 4. + +amr.n_cell = 256 8 64 # dx=dy=dz=100 m, Straka et al 1993 / Xue et al 2000 + +geometry.is_periodic = 0 1 0 -geometry.is_periodic = 1 1 0 +xlo.type = "Inflow" +xhi.type = "Outflow" + +xlo.velocity = 1. 0. 0. +xlo.density = 1.16 +xlo.theta = 300. +xlo.scalar = 0. zlo.type = "SlipWall" zhi.type = "SlipWall" # TIME STEP CONTROL -erf.cfl = 0.9 # cfl number for hyperbolic system +erf.no_substepping = 1 +erf.fixed_dt = 1.e-5 # DIAGNOSTICS & VERBOSITY erf.sum_interval = 1 # timesteps between computing mass @@ -38,23 +51,22 @@ erf.plot_int_1 = 20 # number of timesteps between plotfiles erf.plot_vars_1 = density rhoadv_0 x_velocity y_velocity z_velocity pressure scalar volfrac # SOLVER CHOICE -erf.alpha_T = 0.0 -erf.alpha_C = 0.0 -erf.use_gravity = false +erf.use_gravity = true +erf.use_coriolis = false +erf.les_type = "None" erf.les_type = "None" erf.molec_diff_type = "None" erf.dynamicViscosity = 0.0 +erf.alpha_T = 0.0 +erf.alpha_C = 0.0 + erf.init_type = "uniform" # PROBLEM PARAMETERS -prob.rho_0 = 1.0 -prob.T_0 = 1.0 -prob.A_0 = 1.0 -prob.u_0 = 10.0 -prob.v_0 = 5.0 -prob.rad_0 = 0.125 -prob.uRef = 0.0 - -prob.prob_type = 11 +prob.u_0 = 1.0 + +erf.use_rayleigh_damping = true +prob.dampcoef = 2000. # ==> time scale = 0.0005 s (for testing) +prob.zdamp = 1.0 diff --git a/Exec/DevTests/EB_Test/prob.cpp b/Exec/DevTests/EB_Test/prob.cpp index 1c2ccd558..ad4a5f000 100644 --- a/Exec/DevTests/EB_Test/prob.cpp +++ b/Exec/DevTests/EB_Test/prob.cpp @@ -1,4 +1,5 @@ #include "prob.H" +#include "TerrainMetrics.H" using namespace amrex; @@ -66,98 +67,60 @@ Problem::init_custom_pert( Array4 const& z_vel_pert, Array4 const& /*r_hse*/, Array4 const& /*p_hse*/, - Array4 const& /*z_nd*/, - Array4 const& /*z_cc*/, + Array4 const& z_nd, + Array4 const& z_cc, GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, Array4 const& /*mf_v*/, const SolverChoice& sc) { + const int khi = geomdata.Domain().bigEnd()[2]; + const bool use_moisture = (sc.moisture_type != MoistureType::None); - amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - // Geometry - const Real* prob_lo = geomdata.ProbLo(); - const Real* prob_hi = geomdata.ProbHi(); - const Real* dx = geomdata.CellSize(); - const Real x = prob_lo[0] + (i + 0.5) * dx[0]; - const Real y = prob_lo[1] + (j + 0.5) * dx[1]; - const Real z = prob_lo[2] + (k + 0.5) * dx[2]; - - // Define a point (xc,yc,zc) at the center of the domain - const Real xc = parms.xc_frac * (prob_lo[0] + prob_hi[0]); - const Real yc = parms.yc_frac * (prob_lo[1] + prob_hi[1]); - const Real zc = parms.zc_frac * (prob_lo[2] + prob_hi[2]); - - // Define ellipse parameters - const Real r0 = parms.rad_0 * (prob_hi[0] - prob_lo[0]); - - const Real r3d = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc) + (z-zc)*(z-zc)); - const Real r2d_xy = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc)); - const Real r2d_xz = std::sqrt((x-xc)*(x-xc) + (z-zc)*(z-zc)); - - if (parms.prob_type == 10) - { - // Set scalar = A_0*exp(-10r^2), where r is distance from center of domain, - // + B_0*sin(x) - state_pert(i, j, k, RhoScalar_comp) = parms.A_0 * exp(-10.*r3d*r3d) + parms.B_0*sin(x); - - } else if (parms.prob_type == 11) { - state_pert(i, j, k, RhoScalar_comp) = parms.A_0 * 0.25 * (1.0 + std::cos(PI * std::min(r2d_xy, r0) / r0)); - } else if (parms.prob_type == 12) { - state_pert(i, j, k, RhoScalar_comp) = parms.A_0 * 0.25 * (1.0 + std::cos(PI * std::min(r2d_xz, r0) / r0)); - } else if (parms.prob_type == 13) { - const Real r0_z = parms.rad_0 * (prob_hi[2] - prob_lo[2]); - //ellipse for mapfac shear validation - const Real r2d_xz_ell = std::sqrt((x-xc)*(x-xc)/(r0*r0) + (z-zc)*(z-zc)/(r0_z*r0_z)); - state_pert(i, j, k, RhoScalar_comp) = parms.A_0 * 0.25 * (1.0 + std::cos(PI * std::min(r2d_xz_ell, r0_z) / r0_z)); - } else if (parms.prob_type == 14) { - state_pert(i, j, k, RhoScalar_comp) = std::cos(PI*x); - } else { - // Set scalar = A_0 in a ball of radius r0 and 0 elsewhere - if (r3d < r0) { - state_pert(i, j, k, RhoScalar_comp) = parms.A_0; - } else { - state_pert(i, j, k, RhoScalar_comp) = 0.0; - } - } + AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); - state_pert(i, j, k, RhoScalar_comp) *= parms.rho_0; + ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Set scalar = 0 everywhere + state_pert(i, j, k, RhoScalar_comp) = 0.0; - if (use_moisture) { - state_pert(i, j, k, RhoQ1_comp) = 0.0; - state_pert(i, j, k, RhoQ2_comp) = 0.0; - } - }); + if (use_moisture) { + state_pert(i, j, k, RhoQ1_comp) = 0.0; + state_pert(i, j, k, RhoQ2_comp) = 0.0; + } + }); - // Set the x-velocity - amrex::ParallelFor(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - x_vel_pert(i, j, k) = parms.u_0; + // Set the x-velocity + ParallelFor(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + x_vel_pert(i, j, k) = parms.u_0; + }); - const Real* prob_lo = geomdata.ProbLo(); - const Real* dx = geomdata.CellSize(); - const Real z = prob_lo[2] + (k + 0.5) * dx[2]; + // Set the y-velocity + ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + y_vel_pert(i, j, k) = 0.0; + }); - // Set the x-velocity - x_vel_pert(i, j, k) = parms.u_0 + parms.uRef * - std::log((z + parms.z0)/parms.z0)/ - std::log((parms.zRef +parms.z0)/parms.z0); - }); + const auto dx = geomdata.CellSize(); + GpuArray dxInv; + dxInv[0] = 1. / dx[0]; + dxInv[1] = 1. / dx[1]; + dxInv[2] = 1. / dx[2]; - // Set the y-velocity - amrex::ParallelFor(ybx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - y_vel_pert(i, j, k) = parms.v_0; - }); + // Set the z-velocity from impenetrable condition + ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { +#ifdef ERF_USE_EB + z_vel_pert(i, j, k) = 0.0; +#else + z_vel_pert(i, j, k) = WFromOmega(i, j, k, 0.0, x_vel_pert, y_vel_pert, z_nd, dxInv); +#endif + }); - // Set the z-velocity - amrex::ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - z_vel_pert(i, j, k) = 0.0; - }); + amrex::Gpu::streamSynchronize(); } /** @@ -167,32 +130,47 @@ Problem::init_custom_pert( * */ void - Problem::init_custom_terrain (const amrex::Geometry& /*geom*/, + Problem::init_custom_terrain (const amrex::Geometry& geom, amrex::FArrayBox& z_phys_nd, const amrex::Real& /*time*/) { - // Note that this only sets the terrain value at the ground IF k=0 is in the box - amrex::Print() << "Initializing flat terrain at z=0" << std::endl; - // Bottom of domain int k0 = 0; + // Domain valid box (z_nd is nodal) + const amrex::Box& domain = geom.Domain(); + int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; + + const Real* prob_lo = geom.ProbLo(); + const Real* prob_hi = geom.ProbHi(); + + const Real* dx = geom.CellSize(); + + // User function parameters + Real a = 0.5; + Real num = 8 * a * a * a; + Real xcen = 0.5 * (prob_lo[0] + prob_hi[0]); + // Grown box with no z range amrex::Box bx = z_phys_nd.box(); + bx.setRange(2,0); + amrex::Array4 const& z_arr = z_phys_nd.array(); - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int) { - z_arr(i,j,k0) = 2.5; + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int) + { + // Clip indices for ghost-cells + int ii = amrex::min(amrex::max(i,domlo_x),domhi_x); + + // Location of nodes + Real x = (ii * dx[0] - xcen); + + // WoA Hill in x-direction + Real height = num / (x*x + 4 * a * a); + + // Populate terrain height + z_arr(i,j,k0) = height; + + // z_arr(i,j,k0) = 2.5; }); } - -#if 0 -AMREX_GPU_DEVICE -Real -dhdt(int /*i*/, int /*j*/, - const GpuArray /*dx*/, - const Real /*time_mt*/, const Real /*delta_t*/) -{ - return 0.; -} -#endif diff --git a/Source/Advection/Advection.H b/Source/Advection/Advection.H index c5e9e4334..448c392e8 100644 --- a/Source/Advection/Advection.H +++ b/Source/Advection/Advection.H @@ -19,18 +19,14 @@ void AdvectionSrcForRho (const amrex::Box& bx, const amrex::Array4< amrex::Real>& avg_xmom, // These are being defined const amrex::Array4< amrex::Real>& avg_ymom, // from the rho fluxes const amrex::Array4< amrex::Real>& avg_zmom, - const amrex::Array4& z_nd, + const amrex::Array4& ax_arr, + const amrex::Array4& ay_arr, + const amrex::Array4& az_arr, const amrex::Array4& detJ, const amrex::GpuArray& cellSizeInv, const amrex::Array4& mf_m, const amrex::Array4& mf_u, const amrex::Array4& mf_v, -#ifdef ERF_USE_EB - const amrex::Array4& ax_arr, - const amrex::Array4& ay_arr, - const amrex::Array4& az_arr, -#endif - const bool use_terrain, const amrex::GpuArray, AMREX_SPACEDIM>& flx_arr); /** Compute advection tendency for all scalars other than density and potential temperature */ @@ -46,7 +42,6 @@ void AdvectionSrcForScalars (const amrex::Box& bx, const amrex::Array4& mf_m, const AdvType horiz_adv_type, const AdvType vert_adv_type, const amrex::Real horiz_upw_frac, const amrex::Real vert_upw_frac, - const bool use_terrain, const amrex::GpuArray, AMREX_SPACEDIM>& flx_arr, const amrex::Box& domain, const amrex::BCRec* bc_ptr_h); @@ -60,7 +55,10 @@ void AdvectionSrcForMom (const amrex::Box& bxx, const amrex::Box& bxy, const amr const amrex::Array4& w , const amrex::Array4& rho_u , const amrex::Array4& rho_v, const amrex::Array4& Omega , - const amrex::Array4& z_nd , const amrex::Array4& detJ, + const amrex::Array4& ax, + const amrex::Array4& ay, + const amrex::Array4& az, + const amrex::Array4& detJ, const amrex::GpuArray& cellSizeInv, const amrex::Array4& mf_m, const amrex::Array4& mf_u, @@ -90,10 +88,10 @@ AdvectionSrcForOpenBC_Tangent_Xmom (const amrex::Box& bxx, const amrex::Array4& rho_u, const amrex::Array4& rho_v, const amrex::Array4& Omega, - const amrex::Array4& z_nd, + const amrex::Array4& ax, + const amrex::Array4& az, const amrex::Array4& detJ, const amrex::GpuArray& cellSizeInv, - const bool use_terrain, const bool do_lo=false); /** Compute advection tendencies for v momentum tangent to open bc*/ @@ -105,10 +103,10 @@ AdvectionSrcForOpenBC_Tangent_Ymom (const amrex::Box& bxy, const amrex::Array4& rho_u, const amrex::Array4& rho_v, const amrex::Array4& Omega, - const amrex::Array4& z_nd, + const amrex::Array4& ay, + const amrex::Array4& az, const amrex::Array4& detJ, const amrex::GpuArray& cellSizeInv, - const bool use_terrain, const bool do_lo=false); /** Compute advection tendencies for w momentum tangent to open bc*/ @@ -120,10 +118,12 @@ AdvectionSrcForOpenBC_Tangent_Zmom (const amrex::Box& bxz, const amrex::Array4& rho_u, const amrex::Array4& rho_v, const amrex::Array4& Omega, - const amrex::Array4& z_nd, + const amrex::Array4& ax, + const amrex::Array4& ay, + const amrex::Array4& az, const amrex::Array4& detJ, const amrex::GpuArray& cellSizeInv, - const bool use_terrain, const int domhi_z, + const int domhi_z, const bool do_lo=false); /** Compute advection tendencies for cons tangential to BC (2nd order)*/ @@ -139,7 +139,6 @@ AdvectionSrcForOpenBC_Tangent_Cons (const amrex::Box& bx, const amrex::Array4& avg_zmom, const amrex::Array4& detJ, const amrex::GpuArray& cellSizeInv, - const bool use_terrain, const bool do_lo=false); /** Compute advection tendencies in normal dir for vars tangent to open bc */ diff --git a/Source/Advection/AdvectionSrcForMom.cpp b/Source/Advection/AdvectionSrcForMom.cpp index 3875c40a9..0735f66b3 100644 --- a/Source/Advection/AdvectionSrcForMom.cpp +++ b/Source/Advection/AdvectionSrcForMom.cpp @@ -24,7 +24,9 @@ using namespace amrex; * @param[in] rho_u x-component of the momentum * @param[in] rho_v y-component of the momentum * @param[in] Omega component of the momentum normal to the z-coordinate surface - * @param[in] z_nd height coordinate at nodes + * @param[in] ax Area fraction of x-faces + * @param[in] ay Area fraction of y-faces + * @param[in] az Area fraction of z-faces * @param[in] detJ Jacobian of the metric transformation (= 1 if use_terrain is false) * @param[in] cellSizeInv inverse of the mesh spacing * @param[in] mf_m map factor at cell centers @@ -46,7 +48,9 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz, const Array4& rho_u, const Array4& rho_v, const Array4& Omega, - const Array4& z_nd, + const Array4& ax, + const Array4& ay, + const Array4& az, const Array4& detJ, const GpuArray& cellSizeInv, const Array4& mf_m, @@ -85,6 +89,9 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz, mf_v_inv(i,j,0) = 1. / mf_v(i,j,0); }); +#ifdef ERF_USE_EB + amrex::ignore_unused(use_terrain); +#else if (!use_terrain) { // Inline with 2nd order for efficiency if (horiz_adv_type == AdvType::Centered_2nd && vert_adv_type == AdvType::Centered_2nd) @@ -194,27 +201,30 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz, } } // end of use_terrain == false else - { // now do use_terrain == true +#endif + { // now do use_terrain = true (or ERF_USE_EB) // Inline with 2nd order for efficiency if (horiz_adv_type == AdvType::Centered_2nd && vert_adv_type == AdvType::Centered_2nd) { ParallelFor(bxx, bxy, bxz, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real met_h_zeta_xhi = Compute_h_zeta_AtCellCenter(i ,j ,k ,cellSizeInv,z_nd); - Real xflux_hi = 0.25 * (rho_u(i, j , k) * mf_u_inv(i,j,0) + rho_u(i+1, j , k) * mf_u_inv(i+1,j,0)) * (u(i+1,j,k) + u(i,j,k)) * met_h_zeta_xhi; + Real xflux_hi = 0.25 * (rho_u(i,j,k) * mf_u_inv(i,j,0) + rho_u(i+1,j,k) * mf_u_inv(i+1,j,0)) * + (u(i+1,j,k) + u(i,j,k)) * 0.5 * (ax(i,j,k) + ax(i+1,j,k)); - Real met_h_zeta_xlo = Compute_h_zeta_AtCellCenter(i-1,j ,k ,cellSizeInv,z_nd); - Real xflux_lo = 0.25 * (rho_u(i, j , k) * mf_u_inv(i,j,0) + rho_u(i-1, j , k) * mf_u_inv(i-1,j,0)) * (u(i-1,j,k) + u(i,j,k)) * met_h_zeta_xlo; + Real xflux_lo = 0.25 * (rho_u(i,j,k) * mf_u_inv(i,j,0) + rho_u(i-1,j,k) * mf_u_inv(i-1,j,0)) * + (u(i-1,j,k) + u(i,j,k)) * 0.5 * (ax(i,j,k) + ax(i-1,j,k)); - Real met_h_zeta_yhi = Compute_h_zeta_AtEdgeCenterK(i ,j+1,k ,cellSizeInv,z_nd); - Real yflux_hi = 0.25 * (rho_v(i, j+1, k) * mf_v_inv(i ,j+1,0) + rho_v(i-1, j+1, k) * mf_v_inv(i-1,j+1,0)) * (u(i,j+1,k) + u(i,j,k)) * met_h_zeta_yhi; + Real yflux_hi = 0.25 * (rho_v(i,j+1,k) * mf_v_inv(i,j+1,0) + rho_v(i-1,j+1,k) * mf_v_inv(i-1,j+1,0)) * + (u(i,j+1,k) + u(i,j,k)) * 0.5 * (ay(i,j,k) + ay(i,j+1,k)); - Real met_h_zeta_ylo = Compute_h_zeta_AtEdgeCenterK(i ,j ,k ,cellSizeInv,z_nd); - Real yflux_lo = 0.25 * (rho_v(i, j , k) * mf_v_inv(i ,j ,0) + rho_v(i-1, j , k) * mf_v_inv(i-1,j ,0)) * (u(i,j-1,k) + u(i,j,k)) * met_h_zeta_ylo; + Real yflux_lo = 0.25 * (rho_v(i,j,k) * mf_v_inv(i,j,0) + rho_v(i-1,j,k) * mf_v_inv(i-1,j,0)) * + (u(i,j-1,k) + u(i,j,k)) * 0.5 * (ay(i,j,k) + ay(i,j-1,k)); - Real zflux_hi = 0.25 * (Omega(i, j, k+1) + Omega(i-1, j, k+1)) * (u(i,j,k+1) + u(i,j,k)); - Real zflux_lo = 0.25 * (Omega(i, j, k ) + Omega(i-1, j, k )) * (u(i,j,k-1) + u(i,j,k)); + Real zflux_hi = 0.25 * (Omega(i,j,k+1) + Omega(i-1,j,k+1)) * (u(i,j,k+1) + u(i,j,k)) * + 0.5 * (az(i,j,k+1) + az(i-1,j,k+1)); + Real zflux_lo = 0.25 * (Omega(i,j,k ) + Omega(i-1,j,k )) * (u(i,j,k-1) + u(i,j,k)) * + 0.5 * (az(i,j,k ) + az(i-1,j,k )); Real mfsq = mf_u(i,j,0) * mf_u(i,j,0); @@ -226,20 +236,22 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz, }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real met_h_zeta_xhi = Compute_h_zeta_AtEdgeCenterK(i+1,j ,k ,cellSizeInv,z_nd); - Real xflux_hi = 0.25 * (rho_u(i+1,j ,k) * mf_u_inv(i+1,j ,0) + rho_u(i+1,j-1, k) * mf_u_inv(i+1,j-1,0)) * (v(i+1,j,k) + v(i,j,k)) * met_h_zeta_xhi; + Real xflux_hi = 0.25 * (rho_u(i+1,j,k) * mf_u_inv(i+1,j,0) + rho_u(i+1,j-1, k) * mf_u_inv(i+1,j-1,0)) * + (v(i+1,j,k) + v(i,j,k)) * 0.5 * (ax(i,j,k) + ax(i+1,j,k)); - Real met_h_zeta_xlo = Compute_h_zeta_AtEdgeCenterK(i ,j ,k ,cellSizeInv,z_nd); - Real xflux_lo = 0.25 * (rho_u(i, j , k) * mf_u_inv(i ,j ,0) + rho_u(i ,j-1, k) * mf_u_inv(i-1,j ,0)) * (v(i-1,j,k) + v(i,j,k)) * met_h_zeta_xlo; + Real xflux_lo = 0.25 * (rho_u(i,j,k) * mf_u_inv(i,j,0) + rho_u(i ,j-1, k) * mf_u_inv(i-1,j ,0)) * + (v(i-1,j,k) + v(i,j,k)) * 0.5 * (ax(i,j,k) + ax(i-1,j,k)); - Real met_h_zeta_yhi = Compute_h_zeta_AtCellCenter(i ,j ,k ,cellSizeInv,z_nd); - Real yflux_hi = 0.25 * (rho_v(i ,j+1, k) * mf_v_inv(i ,j+1,0) + rho_v(i ,j ,k) * mf_v_inv(i ,j ,0)) * (v(i,j+1,k) + v(i,j,k)) * met_h_zeta_yhi; + Real yflux_hi = 0.25 * (rho_v(i,j+1,k) * mf_v_inv(i,j+1,0) + rho_v(i ,j ,k) * mf_v_inv(i ,j ,0)) * + (v(i,j+1,k) + v(i,j,k)) * 0.5 * (ay(i,j,k) + ay(i,j+1,k)); - Real met_h_zeta_ylo = Compute_h_zeta_AtCellCenter(i ,j-1,k ,cellSizeInv,z_nd); - Real yflux_lo = 0.25 * (rho_v(i ,j ,k) * mf_v_inv(i ,j ,0) + rho_v(i , j-1, k) * mf_v_inv(i ,j-1,0)) * (v(i,j-1,k) + v(i,j,k)) * met_h_zeta_ylo; + Real yflux_lo = 0.25 * (rho_v(i,j,k) * mf_v_inv(i,j,0) + rho_v(i , j-1, k) * mf_v_inv(i ,j-1,0)) * + (v(i,j-1,k) + v(i,j,k)) * 0.5 * (ay(i,j,k) + ay(i,j-1,k)); - Real zflux_hi = 0.25 * (Omega(i, j, k+1) + Omega(i, j-1, k+1)) * (v(i,j,k+1) + v(i,j,k)); - Real zflux_lo = 0.25 * (Omega(i, j, k ) + Omega(i, j-1, k )) * (v(i,j,k-1) + v(i,j,k)); + Real zflux_hi = 0.25 * (Omega(i,j,k+1) + Omega(i, j-1, k+1)) * (v(i,j,k+1) + v(i,j,k)) * + 0.5 * (az(i,j,k+1) + az(i,j-1,k+1)); + Real zflux_lo = 0.25 * (Omega(i,j,k ) + Omega(i, j-1, k )) * (v(i,j,k-1) + v(i,j,k)) * + 0.5 * (az(i,j,k ) + az(i,j-1,k )); Real mfsq = mf_v(i,j,0) * mf_v(i,j,0); @@ -251,22 +263,23 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz, }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real met_h_zeta_xhi = Compute_h_zeta_AtEdgeCenterJ(i+1,j ,k ,cellSizeInv,z_nd); - Real xflux_hi = 0.25*(rho_u(i+1,j ,k) + rho_u(i+1, j, k-1)) * mf_u_inv(i+1,j ,0) * (w(i+1,j,k) + w(i,j,k)) * met_h_zeta_xhi; + Real xflux_hi = 0.25*(rho_u(i+1,j,k) + rho_u(i+1,j,k-1)) * mf_u_inv(i+1,j ,0) * + (w(i+1,j,k) + w(i,j,k)) * 0.5 * (ax(i+1,j,k) + ax(i+1,j,k-1)); - Real met_h_zeta_xlo = Compute_h_zeta_AtEdgeCenterJ(i ,j ,k ,cellSizeInv,z_nd); - Real xflux_lo = 0.25*(rho_u(i ,j ,k) + rho_u(i , j, k-1)) * mf_u_inv(i ,j ,0) * (w(i-1,j,k) + w(i,j,k)) * met_h_zeta_xlo; + Real xflux_lo = 0.25*(rho_u(i,j,k) + rho_u(i,j,k-1)) * mf_u_inv(i ,j ,0) * + (w(i-1,j,k) + w(i,j,k)) * 0.5 * (ax(i,j,k) + ax(i,j,k-1)); - Real met_h_zeta_yhi = Compute_h_zeta_AtEdgeCenterI(i ,j+1,k ,cellSizeInv,z_nd); - Real yflux_hi = 0.25*(rho_v(i ,j+1,k) + rho_v(i, j+1, k-1)) * mf_v_inv(i ,j+1,0) * (w(i,j+1,k) + w(i,j,k)) * met_h_zeta_yhi; + Real yflux_hi = 0.25*(rho_v(i ,j+1,k) + rho_v(i, j+1, k-1)) * mf_v_inv(i ,j+1,0) * + (w(i,j+1,k) + w(i,j,k)) * 0.5 * (ay(i,j+1,k) + ay(i,j+1,k-1)); - Real met_h_zeta_ylo = Compute_h_zeta_AtEdgeCenterI(i ,j ,k ,cellSizeInv,z_nd); - Real yflux_lo = 0.25*(rho_v(i ,j ,k) + rho_v(i, j , k-1)) * mf_v_inv(i ,j ,0) * (w(i,j-1,k) + w(i,j,k)) * met_h_zeta_ylo; + Real yflux_lo = 0.25*(rho_v(i ,j ,k) + rho_v(i, j , k-1)) * mf_v_inv(i ,j ,0) * + (w(i,j-1,k) + w(i,j,k)) * 0.5 * (ay(i,j,k) + ay(i,j,k-1)); Real zflux_lo = 0.25 * (Omega(i,j,k) + Omega(i,j,k-1)) * (w(i,j,k) + w(i,j,k-1)); - Real zflux_hi = (k == hi_z_face) ? Omega(i,j,k) * w(i,j,k) : - 0.25 * (Omega(i,j,k) + Omega(i,j,k+1)) * (w(i,j,k) + w(i,j,k+1)); + Real zflux_hi = (k == hi_z_face) ? Omega(i,j,k) * w(i,j,k) * az(i,j,k): + 0.25 * (Omega(i,j,k) + Omega(i,j,k+1)) * (w(i,j,k) + w(i,j,k+1)) * + 0.5 * (az(i,j,k) + az(i,j,k+1)); Real mfsq = mf_m(i,j,0) * mf_m(i,j,0); @@ -279,45 +292,45 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz, // Template higher order methods } else { if (horiz_adv_type == AdvType::Centered_2nd) { - AdvectionSrcForMomVert_T(bxx, bxy, bxz, - rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, z_nd, detJ, - cellSizeInv, mf_m, mf_u_inv, mf_v_inv, - horiz_upw_frac, vert_upw_frac, - vert_adv_type, lo_z_face, hi_z_face); + AdvectionSrcForMomVert(bxx, bxy, bxz, + rho_u_rhs, rho_v_rhs, rho_w_rhs, + rho_u, rho_v, Omega, u, v, w, ax, ay, az, detJ, + cellSizeInv, mf_m, mf_u_inv, mf_v_inv, + horiz_upw_frac, vert_upw_frac, + vert_adv_type, lo_z_face, hi_z_face); } else if (horiz_adv_type == AdvType::Upwind_3rd) { - AdvectionSrcForMomVert_T(bxx, bxy, bxz, - rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, z_nd, detJ, - cellSizeInv, mf_m, mf_u_inv, mf_v_inv, - horiz_upw_frac, vert_upw_frac, - vert_adv_type, lo_z_face, hi_z_face); + AdvectionSrcForMomVert(bxx, bxy, bxz, + rho_u_rhs, rho_v_rhs, rho_w_rhs, + rho_u, rho_v, Omega, u, v, w, ax, ay, az, detJ, + cellSizeInv, mf_m, mf_u_inv, mf_v_inv, + horiz_upw_frac, vert_upw_frac, + vert_adv_type, lo_z_face, hi_z_face); } else if (horiz_adv_type == AdvType::Centered_4th) { - AdvectionSrcForMomVert_T(bxx, bxy, bxz, - rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, z_nd, detJ, - cellSizeInv, mf_m, mf_u_inv, mf_v_inv, - horiz_upw_frac, vert_upw_frac, - vert_adv_type, lo_z_face, hi_z_face); + AdvectionSrcForMomVert(bxx, bxy, bxz, + rho_u_rhs, rho_v_rhs, rho_w_rhs, + rho_u, rho_v, Omega, u, v, w, ax, ay, az, detJ, + cellSizeInv, mf_m, mf_u_inv, mf_v_inv, + horiz_upw_frac, vert_upw_frac, + vert_adv_type, lo_z_face, hi_z_face); } else if (horiz_adv_type == AdvType::Upwind_5th) { - AdvectionSrcForMomVert_T(bxx, bxy, bxz, - rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, z_nd, detJ, - cellSizeInv, mf_m, mf_u_inv, mf_v_inv, - horiz_upw_frac, vert_upw_frac, - vert_adv_type, lo_z_face, hi_z_face); + AdvectionSrcForMomVert(bxx, bxy, bxz, + rho_u_rhs, rho_v_rhs, rho_w_rhs, + rho_u, rho_v, Omega, u, v, w, ax, ay, az, detJ, + cellSizeInv, mf_m, mf_u_inv, mf_v_inv, + horiz_upw_frac, vert_upw_frac, + vert_adv_type, lo_z_face, hi_z_face); } else if (horiz_adv_type == AdvType::Centered_6th) { - AdvectionSrcForMomVert_T(bxx, bxy, bxz, - rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, z_nd, detJ, - cellSizeInv, mf_m, mf_u_inv, mf_v_inv, - horiz_upw_frac, vert_upw_frac, - vert_adv_type, lo_z_face, hi_z_face); + AdvectionSrcForMomVert(bxx, bxy, bxz, + rho_u_rhs, rho_v_rhs, rho_w_rhs, + rho_u, rho_v, Omega, u, v, w, ax, ay, az, detJ, + cellSizeInv, mf_m, mf_u_inv, mf_v_inv, + horiz_upw_frac, vert_upw_frac, + vert_adv_type, lo_z_face, hi_z_face); } else { - AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!"); + AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!"); } - } - } + } // higher order + } // terrain // Open bc will be imposed upon all vars (we only access cons here for simplicity) const bool xlo_open = (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::open); @@ -331,24 +344,24 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz, Box tby_xlo, tby_xhi, tby_ylo, tby_yhi; Box tbz_xlo, tbz_xhi, tbz_ylo, tbz_yhi; if (xlo_open) { - if (bxx.smallEnd(0) == domain.smallEnd(0)) { tbx_xlo = makeSlab(bxx,0,domain.smallEnd(0)); tbx_xlo.growLo(0,-1); } - if (bxy.smallEnd(0) == domain.smallEnd(0)) { tby_xlo = makeSlab(bxy,0,domain.smallEnd(0)); } - if (bxz.smallEnd(0) == domain.smallEnd(0)) { tbz_xlo = makeSlab(bxz,0,domain.smallEnd(0)); } + if (bxx.smallEnd(0) == domain.smallEnd(0)) { tbx_xlo = makeSlab(bxx,0,domain.smallEnd(0));} + if (bxy.smallEnd(0) == domain.smallEnd(0)) { tby_xlo = makeSlab(bxy,0,domain.smallEnd(0));} + if (bxz.smallEnd(0) == domain.smallEnd(0)) { tbz_xlo = makeSlab(bxz,0,domain.smallEnd(0));} } if (xhi_open) { - if (bxx.bigEnd(0) == domain.bigEnd(0)+1) { tbx_xhi = makeSlab(bxx,0,domain.bigEnd(0)+1); tbx_xhi.growHi(0,-1); } - if (bxy.bigEnd(0) == domain.bigEnd(0)) { tby_xhi = makeSlab(bxy,0,domain.bigEnd(0) ); } - if (bxz.bigEnd(0) == domain.bigEnd(0)) { tbz_xhi = makeSlab(bxz,0,domain.bigEnd(0) ); } + if (bxx.bigEnd(0) == domain.bigEnd(0)+1) { tbx_xhi = makeSlab(bxx,0,domain.bigEnd(0)+1);} + if (bxy.bigEnd(0) == domain.bigEnd(0)) { tby_xhi = makeSlab(bxy,0,domain.bigEnd(0) );} + if (bxz.bigEnd(0) == domain.bigEnd(0)) { tbz_xhi = makeSlab(bxz,0,domain.bigEnd(0) );} } if (ylo_open) { - if (bxx.smallEnd(1) == domain.smallEnd(1)) { tbx_ylo = makeSlab(bxx,1,domain.smallEnd(1)); } - if (bxy.smallEnd(1) == domain.smallEnd(1)) { tby_ylo = makeSlab(bxy,1,domain.smallEnd(1)); tby_ylo.growLo(1,-1); } - if (bxz.smallEnd(1) == domain.smallEnd(1)) { tbz_ylo = makeSlab(bxz,1,domain.smallEnd(1)); } + if (bxx.smallEnd(1) == domain.smallEnd(1)) { tbx_ylo = makeSlab(bxx,1,domain.smallEnd(1));} + if (bxy.smallEnd(1) == domain.smallEnd(1)) { tby_ylo = makeSlab(bxy,1,domain.smallEnd(1));} + if (bxz.smallEnd(1) == domain.smallEnd(1)) { tbz_ylo = makeSlab(bxz,1,domain.smallEnd(1));} } if (yhi_open) { - if (bxx.bigEnd(1) == domain.bigEnd(1)) { tbx_yhi = makeSlab(bxx,1,domain.bigEnd(1) ); } - if (bxy.bigEnd(1) == domain.bigEnd(1)+1) { tby_yhi = makeSlab(bxy,1,domain.bigEnd(1)+1); tby_yhi.growHi(1,-1); } - if (bxz.bigEnd(1) == domain.bigEnd(1)) { tbz_yhi = makeSlab(bxz,1,domain.bigEnd(1) ); } + if (bxx.bigEnd(1) == domain.bigEnd(1)) { tbx_yhi = makeSlab(bxx,1,domain.bigEnd(1) );} + if (bxy.bigEnd(1) == domain.bigEnd(1)+1) { tby_yhi = makeSlab(bxy,1,domain.bigEnd(1)+1);} + if (bxz.bigEnd(1) == domain.bigEnd(1)) { tbz_yhi = makeSlab(bxz,1,domain.bigEnd(1) );} } // Special advection operator for open BC (bndry normal/tangent operations) @@ -357,46 +370,45 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz, AdvectionSrcForOpenBC_Normal(tbx_xlo, 0, rho_u_rhs, u, cell_data, cellSizeInv, do_lo); AdvectionSrcForOpenBC_Tangent_Ymom(tby_xlo, 0, rho_v_rhs, v, rho_u, rho_v, Omega, - z_nd, detJ, cellSizeInv, - use_terrain, do_lo); + ay, az, detJ, cellSizeInv, + do_lo); AdvectionSrcForOpenBC_Tangent_Zmom(tbz_xlo, 0, rho_w_rhs, w, rho_u, rho_v, Omega, - z_nd, detJ, cellSizeInv, - use_terrain, domhi_z, do_lo); + ax, ay, az, detJ, cellSizeInv, + domhi_z, do_lo); } if (xhi_open) { AdvectionSrcForOpenBC_Normal(tbx_xhi, 0, rho_u_rhs, u, cell_data, cellSizeInv); AdvectionSrcForOpenBC_Tangent_Ymom(tby_xhi, 0, rho_v_rhs, v, rho_u, rho_v, Omega, - z_nd, detJ, cellSizeInv, - use_terrain); + ay, az, detJ, cellSizeInv); + AdvectionSrcForOpenBC_Tangent_Zmom(tbz_xhi, 0, rho_w_rhs, w, rho_u, rho_v, Omega, - z_nd, detJ, cellSizeInv, - use_terrain, domhi_z); + ax, ay, az, detJ, cellSizeInv, + domhi_z); } if (ylo_open) { bool do_lo = true; AdvectionSrcForOpenBC_Tangent_Xmom(tbx_ylo, 1, rho_u_rhs, u, rho_u, rho_v, Omega, - z_nd, detJ, cellSizeInv, - use_terrain, do_lo); + ax, az, detJ, cellSizeInv, + do_lo); AdvectionSrcForOpenBC_Normal(tby_ylo, 1, rho_v_rhs, v, cell_data, cellSizeInv, do_lo); AdvectionSrcForOpenBC_Tangent_Zmom(tbz_ylo, 1, rho_w_rhs, w, rho_u, rho_v, Omega, - z_nd, detJ, cellSizeInv, - use_terrain, domhi_z, do_lo); + ax, ay, az, detJ, cellSizeInv, + domhi_z, do_lo); } if (yhi_open) { AdvectionSrcForOpenBC_Tangent_Xmom(tbx_yhi, 1, rho_u_rhs, u, rho_u, rho_v, Omega, - z_nd, detJ, cellSizeInv, - use_terrain); + ax, az, detJ, cellSizeInv); AdvectionSrcForOpenBC_Normal(tby_yhi, 1, rho_v_rhs, v, cell_data, cellSizeInv); AdvectionSrcForOpenBC_Tangent_Zmom(tbz_yhi, 1, rho_w_rhs, w, rho_u, rho_v, Omega, - z_nd, detJ, cellSizeInv, - use_terrain, domhi_z); + ax, ay, az, detJ, cellSizeInv, + domhi_z); } } diff --git a/Source/Advection/AdvectionSrcForMom_T.H b/Source/Advection/AdvectionSrcForMom_T.H index 05d4c6df2..fdbbf6fca 100644 --- a/Source/Advection/AdvectionSrcForMom_T.H +++ b/Source/Advection/AdvectionSrcForMom_T.H @@ -11,7 +11,9 @@ * @param[in] rho_v y-component of momentum * @param[in] Omega component of the momentum normal to the z-coordinate surface * @param[in] u x-component of velocity - * @param[in] z_nd height coordinate at nodes + * @param[in] ax Area fractions on x-faces + * @param[in] ay Area fractions on y-faces + * @param[in] az Area fractions on z-faces * @param[in] detJ Jacobian of the metric transformation (= 1 if use_terrain is false) * @param[in] cellSizeInv inverse of the mesh spacing * @param[in] mf_u map factor on x-faces @@ -21,19 +23,21 @@ template AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real -AdvectionSrcForXMom_T (int i, int j, int k, - const amrex::Array4& rho_u, - const amrex::Array4& rho_v, - const amrex::Array4& Omega, - const amrex::Array4& z_nd, - const amrex::Array4& detJ, - InterpType_H interp_u_h, - InterpType_V interp_u_v, - const amrex::Real upw_frac_h, - const amrex::Real upw_frac_v, - const amrex::GpuArray& cellSizeInv, - const amrex::Array4& mf_u_inv, - const amrex::Array4& mf_v_inv) +AdvectionSrcForXMom (int i, int j, int k, + const amrex::Array4& rho_u, + const amrex::Array4& rho_v, + const amrex::Array4& Omega, + const amrex::Array4& ax, + const amrex::Array4& ay, + const amrex::Array4& az, + const amrex::Array4& detJ, + InterpType_H interp_u_h, + InterpType_V interp_u_v, + const amrex::Real upw_frac_h, + const amrex::Real upw_frac_v, + const amrex::GpuArray& cellSizeInv, + const amrex::Array4& mf_u_inv, + const amrex::Array4& mf_v_inv) { amrex::Real advectionSrc; auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2]; @@ -54,8 +58,8 @@ AdvectionSrcForXMom_T (int i, int j, int k, rho_u_avg_lo = 0.5 * (rho_u(i-1, j, k) * mf_u_inv(i-1,j ,0) + rho_u(i, j, k) * mf_u_inv(i ,j ,0)); interp_u_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo,upw_frac_h); - amrex::Real centFluxXXNext = rho_u_avg_hi * Compute_h_zeta_AtCellCenter(i ,j ,k ,cellSizeInv,z_nd) * interp_hi; - amrex::Real centFluxXXPrev = rho_u_avg_lo * Compute_h_zeta_AtCellCenter(i-1,j ,k ,cellSizeInv,z_nd) * interp_lo; + amrex::Real centFluxXXNext = rho_u_avg_hi * interp_hi * 0.5 * (ax(i,j,k) + ax(i+1,j,k)); + amrex::Real centFluxXXPrev = rho_u_avg_lo * interp_lo * 0.5 * (ax(i,j,k) + ax(i-1,j,k)); // **************************************************************************************** // Y-fluxes (at edges in k-direction) @@ -66,14 +70,14 @@ AdvectionSrcForXMom_T (int i, int j, int k, rho_v_avg_lo = 0.5 * (rho_v(i, j , k) * mf_v_inv(i ,j ,0) + rho_v(i-1, j , k) * mf_v_inv(i-1,j ,0)); interp_u_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo,upw_frac_h); - amrex::Real edgeFluxXYNext = rho_v_avg_hi * Compute_h_zeta_AtEdgeCenterK(i ,j+1,k ,cellSizeInv,z_nd) * interp_hi; - amrex::Real edgeFluxXYPrev = rho_v_avg_lo * Compute_h_zeta_AtEdgeCenterK(i ,j ,k ,cellSizeInv,z_nd) * interp_lo; + amrex::Real edgeFluxXYNext = rho_v_avg_hi * interp_hi * 0.5 * (ay(i,j+1,k) + ay(i-1,j+1,k)); + amrex::Real edgeFluxXYPrev = rho_v_avg_lo * interp_lo * 0.5 * (ay(i,j ,k) + ay(i-1,j ,k)); // **************************************************************************************** // Z-fluxes (at edges in j-direction) // **************************************************************************************** - Omega_avg_hi = 0.5 * (Omega(i, j, k+1) + Omega(i-1, j, k+1)); - Omega_avg_lo = 0.5 * (Omega(i, j, k ) + Omega(i-1, j, k )); + Omega_avg_hi = 0.5 * (Omega(i, j, k+1) + Omega(i-1, j, k+1)) * 0.5 * (az(i,j,k+1) + az(i-1,j,k+1)); + Omega_avg_lo = 0.5 * (Omega(i, j, k ) + Omega(i-1, j, k )) * 0.5 * (az(i,j,k ) + az(i-1,j,k )); interp_u_v.InterpolateInZ(i,j,k+1,0,interp_hi,Omega_avg_hi,upw_frac_v); interp_u_v.InterpolateInZ(i,j,k ,0,interp_lo,Omega_avg_lo,upw_frac_v); @@ -102,7 +106,9 @@ AdvectionSrcForXMom_T (int i, int j, int k, * @param[in] rho_v y-component of momentum * @param[in] Omega component of the momentum normal to the z-coordinate surface * @param[in] v y-component of velocity - * @param[in] z_nd height coordinate at nodes + * @param[in] ax Area fractions on x-faces + * @param[in] ay Area fractions on y-faces + * @param[in] az Area fractions on z-faces * @param[in] detJ Jacobian of the metric transformation (= 1 if use_terrain is false) * @param[in] cellSizeInv inverse of the mesh spacing * @param[in] mf_u map factor on x-faces @@ -112,19 +118,21 @@ template AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real -AdvectionSrcForYMom_T (int i, int j, int k, - const amrex::Array4& rho_u, - const amrex::Array4& rho_v, - const amrex::Array4& Omega, - const amrex::Array4& z_nd, - const amrex::Array4& detJ, - InterpType_H interp_v_h, - InterpType_V interp_v_v, - const amrex::Real upw_frac_h, - const amrex::Real upw_frac_v, - const amrex::GpuArray& cellSizeInv, - const amrex::Array4& mf_u_inv, - const amrex::Array4& mf_v_inv) +AdvectionSrcForYMom (int i, int j, int k, + const amrex::Array4& rho_u, + const amrex::Array4& rho_v, + const amrex::Array4& Omega, + const amrex::Array4& ax, + const amrex::Array4& ay, + const amrex::Array4& az, + const amrex::Array4& detJ, + InterpType_H interp_v_h, + InterpType_V interp_v_v, + const amrex::Real upw_frac_h, + const amrex::Real upw_frac_v, + const amrex::GpuArray& cellSizeInv, + const amrex::Array4& mf_u_inv, + const amrex::Array4& mf_v_inv) { amrex::Real advectionSrc; auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2]; @@ -138,32 +146,32 @@ AdvectionSrcForYMom_T (int i, int j, int k, // **************************************************************************************** // x-fluxes (at edges in k-direction) // **************************************************************************************** - rho_u_avg_hi = 0.5 * (rho_u(i+1, j, k) * mf_u_inv(i+1,j ,0) + rho_u(i+1, j-1, k) * mf_u_inv(i+1,j-1,0)); - interp_v_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi,upw_frac_h); + rho_u_avg_hi = 0.5 * (rho_u(i+1,j,k) * mf_u_inv(i+1,j,0) + rho_u(i+1,j-1,k) * mf_u_inv(i+1,j-1,0)); + rho_u_avg_lo = 0.5 * (rho_u(i ,j,k) * mf_u_inv(i ,j,0) + rho_u(i ,j-1,k) * mf_u_inv(i ,j-1,0)); - rho_u_avg_lo = 0.5 * (rho_u(i , j, k) * mf_u_inv(i ,j ,0) + rho_u(i , j-1, k) * mf_u_inv(i ,j-1,0)); - interp_v_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo,upw_frac_h); + interp_v_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi,upw_frac_h); + interp_v_h.InterpolateInX(i ,j,k,0,interp_lo,rho_u_avg_lo,upw_frac_h); - amrex::Real edgeFluxYXNext = rho_u_avg_hi * Compute_h_zeta_AtEdgeCenterK(i+1,j ,k ,cellSizeInv,z_nd) * interp_hi; - amrex::Real edgeFluxYXPrev = rho_u_avg_lo * Compute_h_zeta_AtEdgeCenterK(i ,j ,k ,cellSizeInv,z_nd) * interp_lo; + amrex::Real edgeFluxYXNext = rho_u_avg_hi * 0.5 * (ax(i+1,j,k) + ax(i+1,j-1,k)) * interp_hi; + amrex::Real edgeFluxYXPrev = rho_u_avg_lo * 0.5 * (ax(i ,j,k) + ax(i ,j-1,k)) * interp_lo; // **************************************************************************************** // y-fluxes (at cell centers) // **************************************************************************************** - rho_v_avg_hi = 0.5 * (rho_v(i, j, k) * mf_v_inv(i ,j ,0) + rho_v(i, j+1, k) * mf_v_inv(i ,j+1,0)); - interp_v_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi,upw_frac_h); + rho_v_avg_hi = 0.5 * (rho_v(i,j,k) * mf_v_inv(i,j,0) + rho_v(i,j+1,k) * mf_v_inv(i,j+1,0)); + rho_v_avg_lo = 0.5 * (rho_v(i,j,k) * mf_v_inv(i,j,0) + rho_v(i,j-1,k) * mf_v_inv(i,j-1,0)); - rho_v_avg_lo = 0.5 * (rho_v(i, j, k) * mf_v_inv(i ,j ,0) + rho_v(i, j-1, k) * mf_v_inv(i ,j-1,0)); - interp_v_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo,upw_frac_h); + interp_v_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi,upw_frac_h); + interp_v_h.InterpolateInY(i,j ,k,0,interp_lo,rho_v_avg_lo,upw_frac_h); - amrex::Real centFluxYYNext = rho_v_avg_hi * Compute_h_zeta_AtCellCenter(i ,j ,k ,cellSizeInv,z_nd) * interp_hi; - amrex::Real centFluxYYPrev = rho_v_avg_lo * Compute_h_zeta_AtCellCenter(i ,j-1,k ,cellSizeInv,z_nd) * interp_lo; + amrex::Real centFluxYYNext = rho_v_avg_hi * 0.5 * (ay(i,j,k) + ay(i,j+1,k)) * interp_hi; + amrex::Real centFluxYYPrev = rho_v_avg_lo * 0.5 * (ay(i,j,k) + ay(i,j-1,k)) * interp_lo; // **************************************************************************************** // Z-fluxes (at edges in j-direction) // **************************************************************************************** - Omega_avg_hi = 0.5 * (Omega(i, j, k+1) + Omega(i, j-1, k+1)); - Omega_avg_lo = 0.5 * (Omega(i, j, k ) + Omega(i, j-1, k )); + Omega_avg_hi = 0.5 * (Omega(i, j, k+1) + Omega(i, j-1, k+1)) * 0.5 * (az(i,j,k+1) + az(i,j-1,k+1)); + Omega_avg_lo = 0.5 * (Omega(i, j, k ) + Omega(i, j-1, k )) * 0.5 * (az(i,j,k ) + az(i,j-1,k )); interp_v_v.InterpolateInZ(i,j,k+1,0,interp_hi,Omega_avg_hi,upw_frac_v); interp_v_v.InterpolateInZ(i,j,k ,0,interp_lo,Omega_avg_lo,upw_frac_v); @@ -192,7 +200,9 @@ AdvectionSrcForYMom_T (int i, int j, int k, * @param[in] rho_v y-component of momentum * @param[in] Omega component of the momentum normal to the z-coordinate surface * @param[in] w z-component of velocity - * @param[in] z_nd height coordinate at nodes + * @param[in] ax Area fractions on x-faces + * @param[in] ay Area fractions on y-faces + * @param[in] az Area fractions on z-faces * @param[in] detJ Jacobian of the metric transformation (= 1 if use_terrain is false) * @param[in] cellSizeInv inverse of the mesh spacing * @param[in] mf_m map factor on cell centers @@ -206,24 +216,26 @@ template AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real -AdvectionSrcForZMom_T (int i, int j, int k, - const amrex::Array4& rho_u, - const amrex::Array4& rho_v, - const amrex::Array4& Omega, - const amrex::Array4& w, - const amrex::Array4& z_nd, - const amrex::Array4& detJ, - InterpType_H interp_omega_h, - InterpType_V interp_omega_v, - WallInterpType interp_omega_wall, - const amrex::Real upw_frac_h, - const amrex::Real upw_frac_v, - const amrex::GpuArray& cellSizeInv, - const amrex::Array4& mf_m, - const amrex::Array4& mf_u_inv, - const amrex::Array4& mf_v_inv, - const AdvType vert_adv_type, - const int lo_z_face, const int hi_z_face) +AdvectionSrcForZMom (int i, int j, int k, + const amrex::Array4& rho_u, + const amrex::Array4& rho_v, + const amrex::Array4& Omega, + const amrex::Array4& w, + const amrex::Array4& ax, + const amrex::Array4& ay, + const amrex::Array4& az, + const amrex::Array4& detJ, + InterpType_H interp_omega_h, + InterpType_V interp_omega_v, + WallInterpType interp_omega_wall, + const amrex::Real upw_frac_h, + const amrex::Real upw_frac_v, + const amrex::GpuArray& cellSizeInv, + const amrex::Array4& mf_m, + const amrex::Array4& mf_u_inv, + const amrex::Array4& mf_v_inv, + const AdvType vert_adv_type, + const int lo_z_face, const int hi_z_face) { amrex::Real advectionSrc; auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2]; @@ -238,13 +250,13 @@ AdvectionSrcForZMom_T (int i, int j, int k, // x-fluxes (at edges in j-direction) // **************************************************************************************** rho_u_avg_hi = 0.5 * (rho_u(i+1, j, k) + rho_u(i+1, j, k-1)) * mf_u_inv(i+1,j ,0); - interp_omega_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi,upw_frac_h); - rho_u_avg_lo = 0.5 * (rho_u(i , j, k) + rho_u(i , j, k-1)) * mf_u_inv(i ,j ,0); + + interp_omega_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi,upw_frac_h); interp_omega_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo,upw_frac_h); - amrex::Real edgeFluxZXNext = rho_u_avg_hi * Compute_h_zeta_AtEdgeCenterJ(i+1,j ,k ,cellSizeInv,z_nd) * interp_hi; - amrex::Real edgeFluxZXPrev = rho_u_avg_lo * Compute_h_zeta_AtEdgeCenterJ(i ,j ,k ,cellSizeInv,z_nd) * interp_lo; + amrex::Real edgeFluxZXNext = rho_u_avg_hi * 0.5 * (ax(i+1,j,k) + ax(i+1,j,k-1)) * interp_hi; + amrex::Real edgeFluxZXPrev = rho_u_avg_lo * 0.5 * (ax(i ,j,k) + ax(i ,j,k-1)) * interp_lo; // **************************************************************************************** // y-fluxes (at edges in i-direction) @@ -255,14 +267,15 @@ AdvectionSrcForZMom_T (int i, int j, int k, rho_v_avg_lo = 0.5 * (rho_v(i, j , k) + rho_v(i, j , k-1)) * mf_v_inv(i ,j ,0); interp_omega_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo,upw_frac_h); - amrex::Real edgeFluxZYNext = rho_v_avg_hi * Compute_h_zeta_AtEdgeCenterI(i ,j+1,k ,cellSizeInv,z_nd) * interp_hi; - amrex::Real edgeFluxZYPrev = rho_v_avg_lo * Compute_h_zeta_AtEdgeCenterI(i ,j ,k ,cellSizeInv,z_nd) * interp_lo; + amrex::Real edgeFluxZYNext = rho_v_avg_hi * 0.5 * (ay(i,j+1,k) + ay(i,j+1,k-1)) * interp_hi; + amrex::Real edgeFluxZYPrev = rho_v_avg_lo * 0.5 * (ay(i,j ,k) + ay(i,j ,k-1)) * interp_lo; // **************************************************************************************** // z-fluxes (at cell centers) // **************************************************************************************** - Omega_avg_hi = (k == hi_z_face) ? Omega(i,j,k) : 0.5 * (Omega(i,j,k) + Omega(i,j,k+1)); + Omega_avg_hi = (k == hi_z_face) ? Omega(i,j,k) * az(i,j,k) : + 0.25 * (Omega(i,j,k) + Omega(i,j,k+1)) * (az(i,j,k) + az(i,j,k+1)); amrex::Real centFluxZZNext = Omega_avg_hi; // int l_spatial_order_hi = std::min(std::min(vert_spatial_order, 2*(hi_z_face-k)), 2*(k+1)); @@ -289,7 +302,8 @@ AdvectionSrcForZMom_T (int i, int j, int k, // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * - Omega_avg_lo = (k == 0) ? Omega(i,j,k) : 0.5 * (Omega(i,j,k) + Omega(i,j,k-1)); + Omega_avg_lo = (k == 0) ? Omega(i,j,k) * az(i,j,k) : + 0.25 * (Omega(i,j,k) + Omega(i,j,k-1)) * (az(i,j,k) + az(i,j,k-1)); amrex::Real centFluxZZPrev = Omega_avg_lo; // int l_spatial_order_lo = std::min(std::min(vert_spatial_order, 2*(hi_z_face+1-k)), 2*k); @@ -332,26 +346,28 @@ AdvectionSrcForZMom_T (int i, int j, int k, */ template void -AdvectionSrcForMomWrapper_T (const amrex::Box& bxx, const amrex::Box& bxy, const amrex::Box& bxz, - const amrex::Array4& rho_u_rhs, - const amrex::Array4& rho_v_rhs, - const amrex::Array4& rho_w_rhs, - const amrex::Array4& rho_u, - const amrex::Array4& rho_v, - const amrex::Array4& Omega, - const amrex::Array4& u, - const amrex::Array4& v, - const amrex::Array4& w, - const amrex::Array4& z_nd, - const amrex::Array4& detJ, - const amrex::GpuArray& cellSizeInv, - const amrex::Array4& mf_m, - const amrex::Array4& mf_u_inv, - const amrex::Array4& mf_v_inv, - const amrex::Real upw_frac_h, - const amrex::Real upw_frac_v, - const AdvType vert_adv_type, - const int lo_z_face, const int hi_z_face) +AdvectionSrcForMomWrapper (const amrex::Box& bxx, const amrex::Box& bxy, const amrex::Box& bxz, + const amrex::Array4& rho_u_rhs, + const amrex::Array4& rho_v_rhs, + const amrex::Array4& rho_w_rhs, + const amrex::Array4& rho_u, + const amrex::Array4& rho_v, + const amrex::Array4& Omega, + const amrex::Array4& u, + const amrex::Array4& v, + const amrex::Array4& w, + const amrex::Array4& ax, + const amrex::Array4& ay, + const amrex::Array4& az, + const amrex::Array4& detJ, + const amrex::GpuArray& cellSizeInv, + const amrex::Array4& mf_m, + const amrex::Array4& mf_u_inv, + const amrex::Array4& mf_v_inv, + const amrex::Real upw_frac_h, + const amrex::Real upw_frac_v, + const AdvType vert_adv_type, + const int lo_z_face, const int hi_z_face) { // Instantiate the appropriate structs InterpType_H interp_u_h(u); InterpType_V interp_u_v(u); // X-MOM @@ -362,25 +378,25 @@ AdvectionSrcForMomWrapper_T (const amrex::Box& bxx, const amrex::Box& bxy, const amrex::ParallelFor(bxx, bxy, bxz, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - rho_u_rhs(i, j, k) = -AdvectionSrcForXMom_T(i, j, k, rho_u, rho_v, Omega, z_nd, detJ, - interp_u_h, interp_u_v, - upw_frac_h, upw_frac_v, + rho_u_rhs(i, j, k) = -AdvectionSrcForXMom(i, j, k, rho_u, rho_v, Omega, ax, ay, az, detJ, + interp_u_h, interp_u_v, + upw_frac_h, upw_frac_v, cellSizeInv, mf_u_inv, mf_v_inv); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - rho_v_rhs(i, j, k) = -AdvectionSrcForYMom_T(i, j, k, rho_u, rho_v, Omega, z_nd, detJ, - interp_v_h, interp_v_v, - upw_frac_h, upw_frac_v, - cellSizeInv, mf_u_inv, mf_v_inv); + rho_v_rhs(i, j, k) = -AdvectionSrcForYMom(i, j, k, rho_u, rho_v, Omega, ax, ay, az, detJ, + interp_v_h, interp_v_v, + upw_frac_h, upw_frac_v, + cellSizeInv, mf_u_inv, mf_v_inv); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - rho_w_rhs(i, j, k) = -AdvectionSrcForZMom_T(i, j, k, rho_u, rho_v, Omega, w, z_nd, detJ, - interp_w_h, interp_w_v, interp_w_wall, - upw_frac_h, upw_frac_v, - cellSizeInv, mf_m, mf_u_inv, mf_v_inv, - vert_adv_type, lo_z_face, hi_z_face); + rho_w_rhs(i, j, k) = -AdvectionSrcForZMom(i, j, k, rho_u, rho_v, Omega, w, ax, ay, az, detJ, + interp_w_h, interp_w_v, interp_w_wall, + upw_frac_h, upw_frac_v, + cellSizeInv, mf_m, mf_u_inv, mf_v_inv, + vert_adv_type, lo_z_face, hi_z_face); }); } @@ -389,67 +405,69 @@ AdvectionSrcForMomWrapper_T (const amrex::Box& bxx, const amrex::Box& bxy, const */ template void -AdvectionSrcForMomVert_T (const amrex::Box& bxx, const amrex::Box& bxy, const amrex::Box& bxz, - const amrex::Array4& rho_u_rhs, - const amrex::Array4& rho_v_rhs, - const amrex::Array4& rho_w_rhs, - const amrex::Array4& rho_u, - const amrex::Array4& rho_v, - const amrex::Array4& Omega, - const amrex::Array4& u, - const amrex::Array4& v, - const amrex::Array4& w, - const amrex::Array4& z_nd, - const amrex::Array4& detJ, - const amrex::GpuArray& cellSizeInv, - const amrex::Array4& mf_m, - const amrex::Array4& mf_u_inv, - const amrex::Array4& mf_v_inv, - const amrex::Real upw_frac_h, - const amrex::Real upw_frac_v, - const AdvType vert_adv_type, - const int lo_z_face, const int hi_z_face) +AdvectionSrcForMomVert (const amrex::Box& bxx, const amrex::Box& bxy, const amrex::Box& bxz, + const amrex::Array4& rho_u_rhs, + const amrex::Array4& rho_v_rhs, + const amrex::Array4& rho_w_rhs, + const amrex::Array4& rho_u, + const amrex::Array4& rho_v, + const amrex::Array4& Omega, + const amrex::Array4& u, + const amrex::Array4& v, + const amrex::Array4& w, + const amrex::Array4& ax, + const amrex::Array4& ay, + const amrex::Array4& az, + const amrex::Array4& detJ, + const amrex::GpuArray& cellSizeInv, + const amrex::Array4& mf_m, + const amrex::Array4& mf_u_inv, + const amrex::Array4& mf_v_inv, + const amrex::Real upw_frac_h, + const amrex::Real upw_frac_v, + const AdvType vert_adv_type, + const int lo_z_face, const int hi_z_face) { if (vert_adv_type == AdvType::Centered_2nd) { - AdvectionSrcForMomWrapper_T(bxx, bxy, bxz, - rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, z_nd, detJ, - cellSizeInv, mf_m, - mf_u_inv, mf_v_inv, - upw_frac_h, upw_frac_v, - vert_adv_type, lo_z_face, hi_z_face); - } else if (vert_adv_type == AdvType::Upwind_3rd) { - AdvectionSrcForMomWrapper_T(bxx, bxy, bxz, + AdvectionSrcForMomWrapper(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, z_nd, detJ, + rho_u, rho_v, Omega, u, v, w, ax, ay, az, detJ, cellSizeInv, mf_m, mf_u_inv, mf_v_inv, upw_frac_h, upw_frac_v, vert_adv_type, lo_z_face, hi_z_face); + } else if (vert_adv_type == AdvType::Upwind_3rd) { + AdvectionSrcForMomWrapper(bxx, bxy, bxz, + rho_u_rhs, rho_v_rhs, rho_w_rhs, + rho_u, rho_v, Omega, u, v, w, ax, ay, az, detJ, + cellSizeInv, mf_m, + mf_u_inv, mf_v_inv, + upw_frac_h, upw_frac_v, + vert_adv_type, lo_z_face, hi_z_face); } else if (vert_adv_type == AdvType::Centered_4th) { - AdvectionSrcForMomWrapper_T(bxx, bxy, bxz, - rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, z_nd, detJ, - cellSizeInv, mf_m, - mf_u_inv, mf_v_inv, - upw_frac_h, upw_frac_v, - vert_adv_type, lo_z_face, hi_z_face); - } else if (vert_adv_type == AdvType::Upwind_5th) { - AdvectionSrcForMomWrapper_T(bxx, bxy, bxz, + AdvectionSrcForMomWrapper(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, z_nd, detJ, + rho_u, rho_v, Omega, u, v, w, ax, ay, az, detJ, cellSizeInv, mf_m, mf_u_inv, mf_v_inv, upw_frac_h, upw_frac_v, vert_adv_type, lo_z_face, hi_z_face); + } else if (vert_adv_type == AdvType::Upwind_5th) { + AdvectionSrcForMomWrapper(bxx, bxy, bxz, + rho_u_rhs, rho_v_rhs, rho_w_rhs, + rho_u, rho_v, Omega, u, v, w, ax, ay, az, detJ, + cellSizeInv, mf_m, + mf_u_inv, mf_v_inv, + upw_frac_h, upw_frac_v, + vert_adv_type, lo_z_face, hi_z_face); } else if (vert_adv_type == AdvType::Centered_6th) { - AdvectionSrcForMomWrapper_T(bxx, bxy, bxz, - rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, z_nd, detJ, - cellSizeInv, mf_m, - mf_u_inv, mf_v_inv, - upw_frac_h, upw_frac_v, - vert_adv_type, lo_z_face, hi_z_face); + AdvectionSrcForMomWrapper(bxx, bxy, bxz, + rho_u_rhs, rho_v_rhs, rho_w_rhs, + rho_u, rho_v, Omega, u, v, w, ax, ay, az, detJ, + cellSizeInv, mf_m, + mf_u_inv, mf_v_inv, + upw_frac_h, upw_frac_v, + vert_adv_type, lo_z_face, hi_z_face); } else { AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!"); } diff --git a/Source/Advection/AdvectionSrcForOpenBC.cpp b/Source/Advection/AdvectionSrcForOpenBC.cpp index ef7ac4f90..fc7b5a442 100644 --- a/Source/Advection/AdvectionSrcForOpenBC.cpp +++ b/Source/Advection/AdvectionSrcForOpenBC.cpp @@ -44,43 +44,24 @@ AdvectionSrcForOpenBC_Tangent_Xmom (const Box& bxx, const Array4& rho_u, const Array4& rho_v, const Array4& Omega, - const Array4& z_nd, + const Array4& ax, + const Array4& az, const Array4& detJ, const GpuArray& cellSizeInv, - const bool use_terrain, const bool do_lo) { AMREX_ALWAYS_ASSERT(dir==1); auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2]; - if (!use_terrain) { - ParallelFor(bxx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real xflux_hi = 0.25 * (rho_u(i, j , k) + rho_u(i+1, j , k)) * (u(i+1,j,k) + u(i,j,k)); - Real xflux_lo = 0.25 * (rho_u(i, j , k) + rho_u(i-1, j , k)) * (u(i-1,j,k) + u(i,j,k)); - - Real zflux_hi = 0.25 * (Omega(i, j, k+1) + Omega(i-1, j, k+1)) * (u(i,j,k+1) + u(i,j,k)); - Real zflux_lo = 0.25 * (Omega(i, j, k ) + Omega(i-1, j, k )) * (u(i,j,k-1) + u(i,j,k)); - - Real x_src = (xflux_hi - xflux_lo) * dxInv; - Real y_src = AdvectionSrcForOpenBC_Tangent(i, j, k, 0, dir, u, rho_v, dyInv, do_lo); - Real z_src = (zflux_hi - zflux_lo) * dzInv; - Real advectionSrc = x_src + y_src + z_src; - - rho_u_rhs(i, j, k) = -advectionSrc; - }); - } else { - ParallelFor(bxx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real met_h_zeta_xhi = Compute_h_zeta_AtCellCenter(i ,j ,k ,cellSizeInv,z_nd); - Real xflux_hi = 0.25 * (rho_u(i, j , k) + rho_u(i+1, j , k)) * (u(i+1,j,k) + u(i,j,k)) * met_h_zeta_xhi; - - Real met_h_zeta_xlo = Compute_h_zeta_AtCellCenter(i-1,j ,k ,cellSizeInv,z_nd); - Real xflux_lo = 0.25 * (rho_u(i, j , k) + rho_u(i-1, j , k)) * (u(i-1,j,k) + u(i,j,k)) * met_h_zeta_xlo; + ParallelFor(bxx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (ax(i,j,k) > 0.) { + Real xflux_hi = 0.25 * (rho_u(i, j , k) + rho_u(i+1, j , k)) * (u(i+1,j,k) + u(i,j,k)) * 0.5 * (ax(i,j,k) + ax(i+1,j,k)); + Real xflux_lo = 0.25 * (rho_u(i, j , k) + rho_u(i-1, j , k)) * (u(i-1,j,k) + u(i,j,k)) * 0.5 * (ax(i,j,k) + ax(i-1,j,k)); - Real zflux_hi = 0.25 * (Omega(i, j, k+1) + Omega(i-1, j, k+1)) * (u(i,j,k+1) + u(i,j,k)); - Real zflux_lo = 0.25 * (Omega(i, j, k ) + Omega(i-1, j, k )) * (u(i,j,k-1) + u(i,j,k)); + Real zflux_hi = 0.25 * (Omega(i, j, k+1) + Omega(i-1, j, k+1)) * (u(i,j,k+1) + u(i,j,k)) * az(i,j,k+1); + Real zflux_lo = 0.25 * (Omega(i, j, k ) + Omega(i-1, j, k )) * (u(i,j,k-1) + u(i,j,k)) * az(i,j,k ); Real x_src = (xflux_hi - xflux_lo) * dxInv; Real y_src = AdvectionSrcForOpenBC_Tangent(i, j, k, 0, dir, u, rho_v, dyInv, do_lo); @@ -88,8 +69,10 @@ AdvectionSrcForOpenBC_Tangent_Xmom (const Box& bxx, Real advectionSrc = x_src + y_src + z_src; rho_u_rhs(i, j, k) = -advectionSrc / (0.5 * (detJ(i,j,k) + detJ(i-1,j,k))); - }); - } + } else { + rho_u_rhs(i, j, k) = 0.0; + } + }); } /** Compute advection tendencies for x momentum tangential to BC (2nd order)*/ @@ -101,43 +84,24 @@ AdvectionSrcForOpenBC_Tangent_Ymom (const Box& bxy, const Array4& rho_u, const Array4& rho_v, const Array4& Omega, - const Array4& z_nd, + const Array4& ay, + const Array4& az, const Array4& detJ, const GpuArray& cellSizeInv, - const bool use_terrain, const bool do_lo) { AMREX_ALWAYS_ASSERT(dir==0); auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2]; - if (!use_terrain) { - ParallelFor(bxy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real yflux_hi = 0.25 * (rho_v(i ,j+1,k) + rho_v(i ,j ,k)) * (v(i,j+1,k) + v(i,j,k)); - Real yflux_lo = 0.25 * (rho_v(i ,j ,k) + rho_v(i ,j-1,k)) * (v(i,j-1,k) + v(i,j,k)); - - Real zflux_hi = 0.25 * (Omega(i, j, k+1) + Omega(i, j-1, k+1)) * (v(i,j,k+1) + v(i,j,k)); - Real zflux_lo = 0.25 * (Omega(i, j, k ) + Omega(i, j-1, k )) * (v(i,j,k-1) + v(i,j,k)); - - Real x_src = AdvectionSrcForOpenBC_Tangent(i, j, k, 0, dir, v, rho_u, dxInv, do_lo); - Real y_src = (yflux_hi - yflux_lo) * dyInv; - Real z_src = (zflux_hi - zflux_lo) * dzInv; - Real advectionSrc = x_src + y_src + z_src; - - rho_v_rhs(i, j, k) = -advectionSrc; - }); - } else { - ParallelFor(bxy, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real met_h_zeta_yhi = Compute_h_zeta_AtCellCenter(i ,j ,k ,cellSizeInv,z_nd); - Real yflux_hi = 0.25 * (rho_v(i ,j+1, k) + rho_v(i ,j ,k)) * (v(i,j+1,k) + v(i,j,k)) * met_h_zeta_yhi; - - Real met_h_zeta_ylo = Compute_h_zeta_AtCellCenter(i ,j-1,k ,cellSizeInv,z_nd); - Real yflux_lo = 0.25 * (rho_v(i ,j ,k) + rho_v(i , j-1, k)) * (v(i,j-1,k) + v(i,j,k)) * met_h_zeta_ylo; + ParallelFor(bxy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (ay(i,j,k) > 0.) { + Real yflux_hi = 0.25 * (rho_v(i ,j+1,k) + rho_v(i ,j ,k)) * (v(i,j+1,k) + v(i,j,k)) * 0.5 * (ay(i,j,k) + ay(i,j+1,k)); + Real yflux_lo = 0.25 * (rho_v(i ,j ,k) + rho_v(i ,j-1,k)) * (v(i,j-1,k) + v(i,j,k)) * 0.5 * (ay(i,j,k) + ay(i,j-1,k)); - Real zflux_hi = 0.25 * (Omega(i, j, k+1) + Omega(i, j-1, k+1)) * (v(i,j,k+1) + v(i,j,k)); - Real zflux_lo = 0.25 * (Omega(i, j, k ) + Omega(i, j-1, k )) * (v(i,j,k-1) + v(i,j,k)); + Real zflux_hi = 0.25 * (Omega(i, j, k+1) + Omega(i, j-1, k+1)) * (v(i,j,k+1) + v(i,j,k)) * az(i,j,k+1); + Real zflux_lo = 0.25 * (Omega(i, j, k ) + Omega(i, j-1, k )) * (v(i,j,k-1) + v(i,j,k)) * az(i,j,k ); Real x_src = AdvectionSrcForOpenBC_Tangent(i, j, k, 0, dir, v, rho_u, dxInv, do_lo); Real y_src = (yflux_hi - yflux_lo) * dyInv; @@ -145,8 +109,10 @@ AdvectionSrcForOpenBC_Tangent_Ymom (const Box& bxy, Real advectionSrc = x_src + y_src + z_src; rho_v_rhs(i, j, k) = -advectionSrc / (0.5 * (detJ(i,j,k) + detJ(i,j-1,k))); - }); - } + } else { + rho_v_rhs(i, j, k) = 0.0; + } + }); } /** Compute advection tendencies for x momentum tangential to BC (2nd order)*/ @@ -158,10 +124,11 @@ AdvectionSrcForOpenBC_Tangent_Zmom (const Box& bxz, const Array4& rho_u, const Array4& rho_v, const Array4& Omega, - const Array4& z_nd, + const Array4& ax, + const Array4& ay, + const Array4& az, const Array4& detJ, const GpuArray& cellSizeInv, - const bool use_terrain, const int domhi_z, const bool do_lo) { @@ -172,48 +139,27 @@ AdvectionSrcForOpenBC_Tangent_Zmom (const Box& bxz, auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2]; - if (!use_terrain) { - ParallelFor(bxz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real xflux_hi = 0.25*(rho_u(i+1,j ,k) + rho_u(i+1, j, k-1)) * (w(i+1,j,k) + w(i,j,k)); - Real xflux_lo = 0.25*(rho_u(i ,j ,k) + rho_u(i , j, k-1)) * (w(i-1,j,k) + w(i,j,k)); - - Real yflux_hi = 0.25*(rho_v(i ,j+1,k) + rho_v(i, j+1, k-1)) * (w(i,j+1,k) + w(i,j,k)); - Real yflux_lo = 0.25*(rho_v(i ,j ,k) + rho_v(i, j , k-1)) * (w(i,j-1,k) + w(i,j,k)); - - Real zflux_lo = 0.25 * (Omega(i,j,k) + Omega(i,j,k-1)) * (w(i,j,k) + w(i,j,k-1)); - - Real zflux_hi = (k == domhi_z+1) ? Omega(i,j,k) * w(i,j,k) : - 0.25 * (Omega(i,j,k) + Omega(i,j,k+1)) * (w(i,j,k) + w(i,j,k+1)); - - Real x_src = (xopen) ? AdvectionSrcForOpenBC_Tangent(i, j, k, 0, dir, w, rho_u, dxInv, do_lo) : - (xflux_hi - xflux_lo) * dxInv; - Real y_src = (yopen) ? AdvectionSrcForOpenBC_Tangent(i, j, k, 0, dir, w, rho_v, dyInv, do_lo) : - (yflux_hi - yflux_lo) * dyInv; - Real z_src = (zflux_hi - zflux_lo) * dzInv; - Real advectionSrc = x_src + y_src + z_src; - - rho_w_rhs(i, j, k) = -advectionSrc; - }); - } else { - ParallelFor(bxz, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real met_h_zeta_xhi = Compute_h_zeta_AtEdgeCenterJ(i+1,j ,k ,cellSizeInv,z_nd); - Real xflux_hi = 0.25*(rho_u(i+1,j ,k) + rho_u(i+1, j, k-1)) * (w(i+1,j,k) + w(i,j,k)) * met_h_zeta_xhi; + ParallelFor(bxz, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (az(i,j,k) > 0.) { + Real xflux_hi = 0.25*(rho_u(i+1,j ,k) + rho_u(i+1, j, k-1)) * (w(i+1,j,k) + w(i,j,k)) * + 0.5 *(ax(i+1,j,k) + ax(i+1,j,k-1)); - Real met_h_zeta_xlo = Compute_h_zeta_AtEdgeCenterJ(i ,j ,k ,cellSizeInv,z_nd); - Real xflux_lo = 0.25*(rho_u(i ,j ,k) + rho_u(i , j, k-1)) * (w(i-1,j,k) + w(i,j,k)) * met_h_zeta_xlo; + Real xflux_lo = 0.25*(rho_u(i ,j ,k) + rho_u(i , j, k-1)) * (w(i-1,j,k) + w(i,j,k)) * + 0.5 *(ax(i,j,k) + ax(i,j,k-1)); - Real met_h_zeta_yhi = Compute_h_zeta_AtEdgeCenterI(i ,j+1,k ,cellSizeInv,z_nd); - Real yflux_hi = 0.25*(rho_v(i ,j+1,k) + rho_v(i, j+1, k-1)) * (w(i,j+1,k) + w(i,j,k)) * met_h_zeta_yhi; + Real yflux_hi = 0.25*(rho_v(i ,j+1,k) + rho_v(i, j+1, k-1)) * (w(i,j+1,k) + w(i,j,k)) * + 0.5 *(ay(i,j+1,k) + ay(i,j+1,k-1)); - Real met_h_zeta_ylo = Compute_h_zeta_AtEdgeCenterI(i ,j ,k ,cellSizeInv,z_nd); - Real yflux_lo = 0.25*(rho_v(i ,j ,k) + rho_v(i, j , k-1)) * (w(i,j-1,k) + w(i,j,k)) * met_h_zeta_ylo; + Real yflux_lo = 0.25*(rho_v(i ,j ,k) + rho_v(i, j , k-1)) * (w(i,j-1,k) + w(i,j,k)) * + 0.5 *(ay(i,j,k) + ay(i,j,k-1)); - Real zflux_lo = 0.25 * (Omega(i,j,k) + Omega(i,j,k-1)) * (w(i,j,k) + w(i,j,k-1)); + Real zflux_lo = 0.25 * (Omega(i,j,k) + Omega(i,j,k-1)) * (w(i,j,k) + w(i,j,k-1)) * + 0.5 *(az(i,j,k) + az(i,j,k-1)); - Real zflux_hi = (k == domhi_z+1) ? Omega(i,j,k) * w(i,j,k) : - 0.25 * (Omega(i,j,k) + Omega(i,j,k+1)) * (w(i,j,k) + w(i,j,k+1)); + Real zflux_hi = (k == domhi_z+1) ? Omega(i,j,k) * w(i,j,k) * az(i,j,k): + 0.25 * (Omega(i,j,k) + Omega(i,j,k+1)) * (w(i,j,k) + w(i,j,k+1)) * + 0.5 * (az(i,j,k) + az(i,j,k-1)); Real x_src = (xopen) ? AdvectionSrcForOpenBC_Tangent(i, j, k, 0, dir, w, rho_u, dxInv, do_lo) : (xflux_hi - xflux_lo) * dxInv; @@ -223,8 +169,10 @@ AdvectionSrcForOpenBC_Tangent_Zmom (const Box& bxz, Real advectionSrc = x_src + y_src + z_src; rho_w_rhs(i, j, k) = -advectionSrc / (0.5*(detJ(i,j,k) + detJ(i,j,k-1))); - }); - } + } else { + rho_w_rhs(i, j, k) = 0.0; + } + }); } /** Compute advection tendencies for x momentum tangential to BC (2nd order)*/ @@ -240,7 +188,6 @@ AdvectionSrcForOpenBC_Tangent_Cons (const Box& bx, const Array4& avg_zmom, const Array4& detJ, const GpuArray& cellSizeInv, - const bool use_terrain, const bool do_lo) { AMREX_ALWAYS_ASSERT(dir!=2 && icomp>0); @@ -252,33 +199,33 @@ AdvectionSrcForOpenBC_Tangent_Cons (const Box& bx, ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - const int cons_index = icomp + n; - const int prim_index = cons_index - 1; - - Real prim_xlo = 0.5 * (cell_prim(i,j,k,prim_index) + cell_prim(i-1,j,k,prim_index)); - Real prim_xhi = 0.5 * (cell_prim(i,j,k,prim_index) + cell_prim(i+1,j,k,prim_index)); - Real xflux_lo = avg_xmom(i ,j,k) * prim_xlo; - Real xflux_hi = avg_xmom(i+1,j,k) * prim_xhi; - - Real prim_ylo = 0.5 * (cell_prim(i,j,k,prim_index) + cell_prim(i,j-1,k,prim_index)); - Real prim_yhi = 0.5 * (cell_prim(i,j,k,prim_index) + cell_prim(i,j+1,k,prim_index)); - Real yflux_lo = avg_ymom(i,j ,k) * prim_ylo; - Real yflux_hi = avg_ymom(i,j+1,k) * prim_yhi; - - Real prim_zlo = 0.5 * (cell_prim(i,j,k,prim_index) + cell_prim(i,j,k-1,prim_index)); - Real prim_zhi = 0.5 * (cell_prim(i,j,k,prim_index) + cell_prim(i,j,k+1,prim_index)); - Real zflux_lo = avg_zmom(i,j,k ) * prim_zlo; - Real zflux_hi = avg_zmom(i,j,k+1) * prim_zhi; - - Real invdetJ = (use_terrain) ? 1. / detJ(i,j,k) : 1.; - - Real x_src = (xopen) ? AdvectionSrcForOpenBC_Tangent(i, j, k, prim_index, dir, cell_prim, avg_xmom, dxInv, do_lo) : - (xflux_hi - xflux_lo) * dxInv; - Real y_src = (yopen) ? AdvectionSrcForOpenBC_Tangent(i, j, k, prim_index, dir, cell_prim, avg_ymom, dyInv, do_lo) : - (yflux_hi - yflux_lo) * dyInv; - Real z_src = (zflux_hi - zflux_lo) * dzInv; - Real advectionSrc = x_src + y_src + z_src; - cell_rhs(i,j,k,cons_index) = -advectionSrc * invdetJ; + if (detJ(i,j,k) > 0.) { + const int cons_index = icomp + n; + const int prim_index = cons_index - 1; + + Real prim_xlo = 0.5 * (cell_prim(i,j,k,prim_index) + cell_prim(i-1,j,k,prim_index)); + Real prim_xhi = 0.5 * (cell_prim(i,j,k,prim_index) + cell_prim(i+1,j,k,prim_index)); + Real xflux_lo = avg_xmom(i ,j,k) * prim_xlo; + Real xflux_hi = avg_xmom(i+1,j,k) * prim_xhi; + + Real prim_ylo = 0.5 * (cell_prim(i,j,k,prim_index) + cell_prim(i,j-1,k,prim_index)); + Real prim_yhi = 0.5 * (cell_prim(i,j,k,prim_index) + cell_prim(i,j+1,k,prim_index)); + Real yflux_lo = avg_ymom(i,j ,k) * prim_ylo; + Real yflux_hi = avg_ymom(i,j+1,k) * prim_yhi; + + Real prim_zlo = 0.5 * (cell_prim(i,j,k,prim_index) + cell_prim(i,j,k-1,prim_index)); + Real prim_zhi = 0.5 * (cell_prim(i,j,k,prim_index) + cell_prim(i,j,k+1,prim_index)); + Real zflux_lo = avg_zmom(i,j,k ) * prim_zlo; + Real zflux_hi = avg_zmom(i,j,k+1) * prim_zhi; + + Real x_src = (xopen) ? AdvectionSrcForOpenBC_Tangent(i, j, k, prim_index, dir, cell_prim, avg_xmom, dxInv, do_lo) : + (xflux_hi - xflux_lo) * dxInv; + Real y_src = (yopen) ? AdvectionSrcForOpenBC_Tangent(i, j, k, prim_index, dir, cell_prim, avg_ymom, dyInv, do_lo) : + (yflux_hi - yflux_lo) * dyInv; + Real z_src = (zflux_hi - zflux_lo) * dzInv; + Real advectionSrc = x_src + y_src + z_src; + cell_rhs(i,j,k,cons_index) = -advectionSrc / detJ(i,j,k); + } }); } diff --git a/Source/Advection/AdvectionSrcForState.cpp b/Source/Advection/AdvectionSrcForState.cpp index 263927312..9aff1bab5 100644 --- a/Source/Advection/AdvectionSrcForState.cpp +++ b/Source/Advection/AdvectionSrcForState.cpp @@ -19,17 +19,11 @@ using namespace amrex; * @param[out] avg_xmom x-component of time-averaged momentum defined in this routine * @param[out] avg_ymom y-component of time-averaged momentum defined in this routine * @param[out] avg_zmom z-component of time-averaged momentum defined in this routine - * @param[in] z_nd height coordinate at nodes * @param[in] detJ Jacobian of the metric transformation (= 1 if use_terrain is false) * @param[in] cellSizeInv inverse of the mesh spacing * @param[in] mf_m map factor at cell centers * @param[in] mf_u map factor at x-faces * @param[in] mf_v map factor at y-faces - * @param[in] horiz_adv_type advection scheme to be used in horiz. directions for dry scalars - * @param[in] vert_adv_type advection scheme to be used in vert. directions for dry scalars - * @param[in] horiz_upw_frac upwinding fraction to be used in horiz. directions for dry scalars (for Blended schemes only) - * @param[in] vert_upw_frac upwinding fraction to be used in vert. directions for dry scalars (for Blended schemes only) - * @param[in] use_terrain if true, use the terrain-aware derivatives (with metric terms) */ void @@ -41,18 +35,14 @@ AdvectionSrcForRho (const Box& bx, const Array4< Real>& avg_xmom, const Array4< Real>& avg_ymom, const Array4< Real>& avg_zmom, - const Array4& z_nd, + const Array4& ax_arr, + const Array4& ay_arr, + const Array4& az_arr, const Array4& detJ, const GpuArray& cellSizeInv, const Array4& mf_m, const Array4& mf_u, const Array4& mf_v, -#ifdef ERF_USE_EB - const Array4& ax_arr, - const Array4& ay_arr, - const Array4& az_arr, -#endif - const bool use_terrain, const GpuArray, AMREX_SPACEDIM>& flx_arr) { BL_PROFILE_VAR("AdvectionSrcForRhoAndTheta", AdvectionSrcForRhoAndTheta); @@ -64,68 +54,32 @@ AdvectionSrcForRho (const Box& bx, ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - (flx_arr[0])(i,j,k,0) = rho_u(i,j,k) / mf_u(i,j,0); + (flx_arr[0])(i,j,k,0) = ax_arr(i,j,k) * rho_u(i,j,k) / mf_u(i,j,0); avg_xmom(i,j,k) = (flx_arr[0])(i,j,k,0); }); ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - (flx_arr[1])(i,j,k,0) = rho_v(i,j,k) / mf_v(i,j,0); + (flx_arr[1])(i,j,k,0) = ay_arr(i,j,k) * rho_v(i,j,k) / mf_v(i,j,0); avg_ymom(i,j,k) = (flx_arr[1])(i,j,k,0); }); ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real mfsq = mf_m(i,j,0) * mf_m(i,j,0); - (flx_arr[2])(i,j,k,0) = Omega(i,j,k) / mfsq; + (flx_arr[2])(i,j,k,0) = az_arr(i,j,k) * Omega(i,j,k) / mfsq; avg_zmom(i,j,k ) = (flx_arr[2])(i,j,k,0); }); -#ifdef ERF_USE_EB - ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - (flx_arr[0])(i,j,k,0) *= ax_arr(i,j,k); - avg_xmom(i,j,k) *= ax_arr(i,j,k); - }); - ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - (flx_arr[1])(i,j,k,0) *= ay_arr(i,j,k); - avg_ymom(i,j,k) *= ay_arr(i,j,k); - }); - ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - (flx_arr[2])(i,j,k,0) *= az_arr(i,j,k); - avg_zmom(i,j,k) *= az_arr(i,j,k); - }); -#else - if (use_terrain) { - ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real h_zeta = Compute_h_zeta_AtIface(i,j,k,cellSizeInv,z_nd); - (flx_arr[0])(i,j,k,0) *= h_zeta; - avg_xmom(i,j,k) *= h_zeta; - }); - ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real h_zeta = Compute_h_zeta_AtJface(i,j,k,cellSizeInv,z_nd); - (flx_arr[1])(i,j,k,0) *= h_zeta; - avg_ymom(i,j,k) *= h_zeta; - }); - } -#endif - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { -#ifdef ERF_USE_EB - Real invdetJ = (detJ(i,j,k) > 0.) ? 1. / detJ(i,j,k) : 1.; -#else - Real invdetJ = (use_terrain) ? 1. / detJ(i,j,k) : 1.; -#endif - - Real mfsq = mf_m(i,j,0) * mf_m(i,j,0); - - advectionSrc(i,j,k,0) = - invdetJ * mfsq * ( - ( (flx_arr[0])(i+1,j,k,0) - (flx_arr[0])(i ,j,k,0) ) * dxInv + - ( (flx_arr[1])(i,j+1,k,0) - (flx_arr[1])(i,j ,k,0) ) * dyInv + - ( (flx_arr[2])(i,j,k+1,0) - (flx_arr[2])(i,j,k ,0) ) * dzInv ); + if (detJ(i,j,k) > 0.) { + Real mfsq = mf_m(i,j,0) * mf_m(i,j,0); + advectionSrc(i,j,k,0) = - mfsq / detJ(i,j,k) * ( + ( (flx_arr[0])(i+1,j,k,0) - (flx_arr[0])(i ,j,k,0) ) * dxInv + + ( (flx_arr[1])(i,j+1,k,0) - (flx_arr[1])(i,j ,k,0) ) * dyInv + + ( (flx_arr[2])(i,j,k+1,0) - (flx_arr[2])(i,j,k ,0) ) * dzInv ); + } else { + advectionSrc(i,j,k,0) = 0.; + } }); } @@ -150,7 +104,6 @@ AdvectionSrcForRho (const Box& bx, * @param[in] vert_adv_type advection scheme to be used in vert. directions for dry scalars * @param[in] horiz_upw_frac upwinding fraction to be used in horiz. directions for dry scalars (for Blended schemes only) * @param[in] vert_upw_frac upwinding fraction to be used in vert. directions for dry scalars (for Blended schemes only) - * @param[in] use_terrain if true, use the terrain-aware derivatives (with metric terms) */ void @@ -167,7 +120,6 @@ AdvectionSrcForScalars (const Box& bx, const int icomp, const int ncomp, const AdvType vert_adv_type, const Real horiz_upw_frac, const Real vert_upw_frac, - const bool use_terrain, const GpuArray, AMREX_SPACEDIM>& flx_arr, const Box& domain, const BCRec* bc_ptr_h) @@ -288,11 +240,7 @@ AdvectionSrcForScalars (const Box& bx, const int icomp, const int ncomp, ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { -#ifdef ERF_USE_EB Real invdetJ = (detJ(i,j,k) > 0.) ? 1. / detJ(i,j,k) : 1.; -#else - Real invdetJ = (use_terrain) ? 1. / detJ(i,j,k) : 1.; -#endif Real mfsq = mf_m(i,j,0) * mf_m(i,j,0); @@ -308,22 +256,22 @@ AdvectionSrcForScalars (const Box& bx, const int icomp, const int ncomp, bool do_lo = true; AdvectionSrcForOpenBC_Tangent_Cons(bx_xlo, 0, icomp, ncomp, advectionSrc, cell_prim, avg_xmom, avg_ymom, avg_zmom, - detJ, cellSizeInv, use_terrain, do_lo); + detJ, cellSizeInv, do_lo); } if (xhi_open) { AdvectionSrcForOpenBC_Tangent_Cons(bx_xhi, 0, icomp, ncomp, advectionSrc, cell_prim, avg_xmom, avg_ymom, avg_zmom, - detJ, cellSizeInv, use_terrain); + detJ, cellSizeInv); } if (ylo_open) { bool do_lo = true; AdvectionSrcForOpenBC_Tangent_Cons(bx_ylo, 1, icomp, ncomp, advectionSrc, cell_prim, avg_xmom, avg_ymom, avg_zmom, - detJ, cellSizeInv, use_terrain, do_lo); + detJ, cellSizeInv, do_lo); } if (yhi_open) { AdvectionSrcForOpenBC_Tangent_Cons(bx_yhi, 1, icomp, ncomp, advectionSrc, cell_prim, avg_xmom, avg_ymom, avg_zmom, - detJ, cellSizeInv, use_terrain); + detJ, cellSizeInv); } } diff --git a/Source/BoundaryConditions/MOSTAverage.cpp b/Source/BoundaryConditions/MOSTAverage.cpp index e2dfc0d99..b2a5f72e8 100644 --- a/Source/BoundaryConditions/MOSTAverage.cpp +++ b/Source/BoundaryConditions/MOSTAverage.cpp @@ -270,7 +270,7 @@ MOSTAverage::set_k_indices_N() AMREX_ASSERT_WITH_MESSAGE(m_zref <= m_zhi - 0.5 * m_dz, "Query point must be below the last z-cell!"); - int lk = static_cast(floor((m_zref - m_zlo) / m_dz)); + int lk = static_cast(floor((m_zref - m_zlo) / m_dz - 0.5)); m_zref = (lk + 0.5) * m_dz + m_zlo; diff --git a/Source/Diffusion/ComputeStrain_T.cpp b/Source/Diffusion/ComputeStrain_T.cpp index 6dafd0587..bc069d232 100644 --- a/Source/Diffusion/ComputeStrain_T.cpp +++ b/Source/Diffusion/ComputeStrain_T.cpp @@ -35,7 +35,8 @@ ComputeStrain_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain, Array4& tau12, Array4& tau13, Array4& tau21, Array4& tau23, Array4& tau31, Array4& tau32, - const Array4& z_nd , + const Array4& z_nd, + const Array4& detJ, const BCRec* bc_ptr, const GpuArray& dxInv, const Array4& /*mf_m*/, const Array4& mf_u, @@ -337,7 +338,7 @@ ComputeStrain_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain, Real met_h_xi,met_h_eta,met_h_zeta; met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_nd); met_h_eta = Compute_h_eta_AtCellCenter (i,j,k,dxInv,z_nd); - met_h_zeta = Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_nd); + met_h_zeta = detJ(i,j,k); tau11(i,j,k) = ( (u(i+1, j, k)/mf_u(i+1,j,0) - u(i, j, k)/mf_u(i,j,0))*dxInv[0] - (met_h_xi/met_h_zeta)*GradUz ) * mf_u(i,j,0)*mf_u(i,j,0); @@ -357,7 +358,8 @@ ComputeStrain_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain, Real met_h_xi,met_h_eta,met_h_zeta; met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd); met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd); - met_h_zeta = Compute_h_zeta_AtEdgeCenterK(i,j,k,dxInv,z_nd); + // met_h_zeta = Compute_h_zeta_AtEdgeCenterK(i,j,k,dxInv,z_nd); + met_h_zeta = 0.5 * (detJ(i,j,k) + detJ(i,j,k-1)); tau12(i,j,k) = 0.5 * ( (u(i, j, k)/mf_u(i,j,0) - u(i , j-1, k)/mf_u(i,j-1,0))*dxInv[1] + (v(i, j, k)/mf_v(i,j,0) - v(i-1, j , k)/mf_v(i-1,j,0))*dxInv[0] @@ -443,7 +445,7 @@ ComputeStrain_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain, Real met_h_xi,met_h_eta,met_h_zeta; met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_nd); met_h_eta = Compute_h_eta_AtCellCenter (i,j,k,dxInv,z_nd); - met_h_zeta = Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_nd); + met_h_zeta = detJ(i,j,k); tau11(i,j,k) = ( (u(i+1, j, k)/mf_u(i+1,j,0) - u(i, j, k)/mf_u(i,j,0))*dxInv[0] - (met_h_xi/met_h_zeta)*GradUz ) * mf_u(i,j,0)*mf_u(i,j,0); diff --git a/Source/Diffusion/ComputeStress_T.cpp b/Source/Diffusion/ComputeStress_T.cpp index d1f6f5628..265bf7ca9 100644 --- a/Source/Diffusion/ComputeStress_T.cpp +++ b/Source/Diffusion/ComputeStress_T.cpp @@ -33,7 +33,8 @@ ComputeStressConsVisc_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff, Array4& tau21, Array4& tau23, Array4& tau31, Array4& tau32, const Array4& er_arr, - const Array4& z_nd , + const Array4& z_nd, + const Array4& detJ, const GpuArray& dxInv) { // Handle constant alpha case, in which the provided mu_eff is actually @@ -257,8 +258,7 @@ ComputeStressConsVisc_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff, ParallelFor(bxcc,tbxxy, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real met_h_zeta; - met_h_zeta = Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_nd); + Real met_h_zeta = detJ(i,j,k); Real mu_tot = rhoAlpha(i,j,k); tau11(i,j,k) *= -mu_tot*met_h_zeta; @@ -266,8 +266,7 @@ ComputeStressConsVisc_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff, }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real met_h_zeta; - met_h_zeta = Compute_h_zeta_AtEdgeCenterK(i,j,k,dxInv,z_nd); + Real met_h_zeta = Compute_h_zeta_AtEdgeCenterK(i,j,k,dxInv,z_nd); Real mu_tot = 0.25*( rhoAlpha(i-1, j , k) + rhoAlpha(i, j , k) + rhoAlpha(i-1, j-1, k) + rhoAlpha(i, j-1, k) ); @@ -309,7 +308,8 @@ ComputeStressVarVisc_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff, Array4& tau21, Array4& tau23, Array4& tau31, Array4& tau32, const Array4& er_arr, - const Array4& z_nd , + const Array4& z_nd, + const Array4& detJ, const GpuArray& dxInv) { // Handle constant alpha case, in which the provided mu_eff is actually @@ -554,8 +554,7 @@ ComputeStressVarVisc_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff, ParallelFor(bxcc,tbxxy, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real met_h_zeta; - met_h_zeta = Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_nd); + Real met_h_zeta = detJ(i,j,k); Real mu_tot = rhoAlpha(i,j,k) + 2.0*mu_turb(i, j, k, EddyDiff::Mom_h); @@ -564,8 +563,7 @@ ComputeStressVarVisc_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff, }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real met_h_zeta; - met_h_zeta = Compute_h_zeta_AtEdgeCenterK(i,j,k,dxInv,z_nd); + Real met_h_zeta = Compute_h_zeta_AtEdgeCenterK(i,j,k,dxInv,z_nd); Real mu_bar = 0.25*( mu_turb(i-1, j , k, EddyDiff::Mom_h) + mu_turb(i, j , k, EddyDiff::Mom_h) + mu_turb(i-1, j-1, k, EddyDiff::Mom_h) + mu_turb(i, j-1, k, EddyDiff::Mom_h) ); diff --git a/Source/Diffusion/Diffusion.H b/Source/Diffusion/Diffusion.H index 578a80809..ce4aefc81 100644 --- a/Source/Diffusion/Diffusion.H +++ b/Source/Diffusion/Diffusion.H @@ -76,6 +76,9 @@ void DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain, const amrex::Array4& yflux, const amrex::Array4& zflux, const amrex::Array4& z_nd, + const amrex::Array4& ax, + const amrex::Array4& ay, + const amrex::Array4& az, const amrex::Array4& detJ, const amrex::GpuArray& dxInv, const amrex::Array4& SmnSmn_a, @@ -107,7 +110,8 @@ void ComputeStressConsVisc_T (amrex::Box bxcc, amrex::Box tbxxy, amrex::Box tbxx amrex::Array4& tau21, amrex::Array4& tau23, amrex::Array4& tau31, amrex::Array4& tau32, const amrex::Array4& er_arr, - const amrex::Array4& z_nd , + const amrex::Array4& z_nd, + const amrex::Array4& detJ, const amrex::GpuArray& dxInv); @@ -127,7 +131,8 @@ void ComputeStressVarVisc_T (amrex::Box bxcc, amrex::Box tbxxy, amrex::Box tbxxz amrex::Array4& tau21, amrex::Array4& tau23, amrex::Array4& tau31, amrex::Array4& tau32, const amrex::Array4& er_arr, - const amrex::Array4& z_nd , + const amrex::Array4& z_nd, + const amrex::Array4& detJ, const amrex::GpuArray& dxInv); @@ -149,7 +154,8 @@ void ComputeStrain_T (amrex::Box bxcc, amrex::Box tbxxy, amrex::Box tbxxz, amrex amrex::Array4& tau12, amrex::Array4& tau13, amrex::Array4& tau21, amrex::Array4& tau23, amrex::Array4& tau31, amrex::Array4& tau32, - const amrex::Array4& z_nd , + const amrex::Array4& z_nd, + const amrex::Array4& detJ, const amrex::BCRec* bc_ptr, const amrex::GpuArray& dxInv, const amrex::Array4& mf_m, const amrex::Array4& mf_u, const amrex::Array4& mf_v); #endif diff --git a/Source/Diffusion/DiffusionSrcForState_T.cpp b/Source/Diffusion/DiffusionSrcForState_T.cpp index 1d7c79048..ebd653050 100644 --- a/Source/Diffusion/DiffusionSrcForState_T.cpp +++ b/Source/Diffusion/DiffusionSrcForState_T.cpp @@ -49,6 +49,9 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, const Array4& yflux, const Array4& zflux, const Array4& z_nd, + const Array4& ax, + const Array4& ay, + const Array4& az, const Array4& detJ, const GpuArray& cellSizeInv, const Array4& SmnSmn_a, @@ -203,9 +206,8 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, rhoAlpha += 0.5 * ( mu_turb(i , j, k, d_eddy_diff_idx[prim_index]) + mu_turb(i-1, j, k, d_eddy_diff_idx[prim_index]) ); - Real met_h_xi,met_h_zeta; - met_h_xi = Compute_h_xi_AtIface (i,j,k,cellSizeInv,z_nd); - met_h_zeta = Compute_h_zeta_AtIface(i,j,k,cellSizeInv,z_nd); + Real met_h_xi = Compute_h_xi_AtIface (i,j,k,cellSizeInv,z_nd); + Real met_h_zeta = ax(i,j,k); Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i-1, j, k+1, prim_index) - cell_prim(i, j, k-1, prim_index) - cell_prim(i-1, j, k-1, prim_index) ); @@ -223,9 +225,8 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, rhoAlpha += 0.5 * ( mu_turb(i, j , k, d_eddy_diff_idy[prim_index]) + mu_turb(i, j-1, k, d_eddy_diff_idy[prim_index]) ); - Real met_h_eta,met_h_zeta; - met_h_eta = Compute_h_eta_AtJface (i,j,k,cellSizeInv,z_nd); - met_h_zeta = Compute_h_zeta_AtJface(i,j,k,cellSizeInv,z_nd); + Real met_h_eta = Compute_h_eta_AtJface (i,j,k,cellSizeInv,z_nd); + Real met_h_zeta = ay(i,j,k); Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i, j-1, k+1, prim_index) - cell_prim(i, j, k-1, prim_index) - cell_prim(i, j-1, k-1, prim_index) ); @@ -243,8 +244,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_index]) + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_index]) ); - Real met_h_zeta; - met_h_zeta = Compute_h_zeta_AtKface(i,j,k,cellSizeInv,z_nd); + Real met_h_zeta = az(i,j,k); Real GradCz; bool ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir) && k == 0); @@ -285,9 +285,8 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, rhoAlpha += 0.5 * ( mu_turb(i , j, k, d_eddy_diff_idx[prim_index]) + mu_turb(i-1, j, k, d_eddy_diff_idx[prim_index]) ); - Real met_h_xi,met_h_zeta; - met_h_xi = Compute_h_xi_AtIface (i,j,k,cellSizeInv,z_nd); - met_h_zeta = Compute_h_zeta_AtIface(i,j,k,cellSizeInv,z_nd); + Real met_h_xi = Compute_h_xi_AtIface (i,j,k,cellSizeInv,z_nd); + Real met_h_zeta = ax(i,j,k); Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i-1, j, k+1, prim_index) - cell_prim(i, j, k-1, prim_index) - cell_prim(i-1, j, k-1, prim_index) ); @@ -304,9 +303,8 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, rhoAlpha += 0.5 * ( mu_turb(i, j , k, d_eddy_diff_idy[prim_index]) + mu_turb(i, j-1, k, d_eddy_diff_idy[prim_index]) ); - Real met_h_eta,met_h_zeta; - met_h_eta = Compute_h_eta_AtJface (i,j,k,cellSizeInv,z_nd); - met_h_zeta = Compute_h_zeta_AtJface(i,j,k,cellSizeInv,z_nd); + Real met_h_eta = Compute_h_eta_AtJface (i,j,k,cellSizeInv,z_nd); + Real met_h_zeta = ay(i,j,k); Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i, j-1, k+1, prim_index) - cell_prim(i, j, k-1, prim_index) - cell_prim(i, j-1, k-1, prim_index) ); @@ -324,8 +322,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_index]) + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_index]) ); - Real met_h_zeta; - met_h_zeta = Compute_h_zeta_AtKface(i,j,k,cellSizeInv,z_nd); + Real met_h_zeta = az(i,j,k); Real GradCz; bool ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir) && k == 0); @@ -366,9 +363,8 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) ); Real rhoAlpha = rhoFace * d_alpha_eff[prim_index]; - Real met_h_xi,met_h_zeta; - met_h_xi = Compute_h_xi_AtIface (i,j,k,cellSizeInv,z_nd); - met_h_zeta = Compute_h_zeta_AtIface(i,j,k,cellSizeInv,z_nd); + Real met_h_xi = Compute_h_xi_AtIface (i,j,k,cellSizeInv,z_nd); + Real met_h_zeta = ax(i,j,k); Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i-1, j, k+1, prim_index) - cell_prim(i, j, k-1, prim_index) - cell_prim(i-1, j, k-1, prim_index) ); @@ -384,9 +380,8 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) ); Real rhoAlpha = rhoFace * d_alpha_eff[prim_index]; - Real met_h_eta,met_h_zeta; - met_h_eta = Compute_h_eta_AtJface (i,j,k,cellSizeInv,z_nd); - met_h_zeta = Compute_h_zeta_AtJface(i,j,k,cellSizeInv,z_nd); + Real met_h_eta = Compute_h_eta_AtJface (i,j,k,cellSizeInv,z_nd); + Real met_h_zeta = ay(i,j,k); Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i, j-1, k+1, prim_index) - cell_prim(i, j, k-1, prim_index) - cell_prim(i, j-1, k-1, prim_index) ); @@ -402,8 +397,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) ); Real rhoAlpha = rhoFace * d_alpha_eff[prim_index]; - Real met_h_zeta; - met_h_zeta = Compute_h_zeta_AtKface(i,j,k,cellSizeInv,z_nd); + Real met_h_zeta = az(i,j,k); Real GradCz; bool ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir) && k == 0); @@ -653,7 +647,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, bool v_ext_dir_on_zhi = ( (bc_ptr[BCVars::yvel_bc].lo(5) == ERFBCType::ext_dir) ); ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - const Real met_h_zeta = Compute_h_zeta_AtCellCenter(i,j,k,cellSizeInv,z_nd); + const Real met_h_zeta = detJ(i,j,k); cell_rhs(i, j, k, qty_index) += ComputeQKESourceTerms(i,j,k,u,v,cell_data,cell_prim, mu_turb,cellSizeInv,domain,pbl_mynn_B1_l,tm_arr(i,j,0), c_ext_dir_on_zlo, c_ext_dir_on_zhi, diff --git a/Source/EB/TerrainIF.H b/Source/EB/TerrainIF.H index 89001e218..066a6be45 100644 --- a/Source/EB/TerrainIF.H +++ b/Source/EB/TerrainIF.H @@ -17,7 +17,9 @@ public: TerrainIF (amrex::FArrayBox& a_z_terrain, amrex::Geometry& a_geom) : m_terr(a_z_terrain), m_geom(a_geom) - {} + { + amrex::Print() << " EB type = Terrain " << std::endl; + } AMREX_GPU_HOST_DEVICE inline amrex::Real operator() (AMREX_D_DECL(amrex::Real x, amrex::Real y, amrex::Real z)) @@ -31,6 +33,8 @@ public: amrex::Array4 const& terr_arr = m_terr.const_array(); + if (i == 128 and j == 3) amrex::Print() <<" TERRARR " << amrex::IntVect(i,j,0) << " " << z << " " << terr_arr(i,j,0) << std::endl; + return -(z - terr_arr(i,j,0)); } diff --git a/Source/ERF.H b/Source/ERF.H index 5497b9c22..8433ce9bc 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -635,13 +635,23 @@ private: amrex::Vector zlevels_stag; // nominal height levels amrex::Vector> z_phys_nd; amrex::Vector> z_phys_cc; + amrex::Vector> detJ_cc; + amrex::Vector> ax; + amrex::Vector> ay; + amrex::Vector> az; amrex::Vector> z_phys_nd_src; amrex::Vector> detJ_cc_src; + amrex::Vector> ax_src; + amrex::Vector> ay_src; + amrex::Vector> az_src; amrex::Vector> z_phys_nd_new; amrex::Vector> detJ_cc_new; + amrex::Vector> ax_new; + amrex::Vector> ay_new; + amrex::Vector> az_new; amrex::Vector> z_t_rk; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 2cc6cea9a..b74cac1d9 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -221,10 +221,22 @@ ERF::ERF () z_phys_nd.resize(nlevs_max); z_phys_cc.resize(nlevs_max); detJ_cc.resize(nlevs_max); + ax.resize(nlevs_max); + ay.resize(nlevs_max); + az.resize(nlevs_max); + z_phys_nd_new.resize(nlevs_max); detJ_cc_new.resize(nlevs_max); + ax_new.resize(nlevs_max); + ay_new.resize(nlevs_max); + az_new.resize(nlevs_max); + z_phys_nd_src.resize(nlevs_max); detJ_cc_src.resize(nlevs_max); + ax_src.resize(nlevs_max); + ay_src.resize(nlevs_max); + az_src.resize(nlevs_max); + z_t_rk.resize(nlevs_max); // Mapfactors diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index 29155bee6..316217410 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/Source/ERF_make_new_arrays.cpp @@ -56,11 +56,20 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, // ******************************************************************************************** if (solverChoice.use_terrain) { z_phys_cc[lev] = std::make_unique(ba,dm,1,1); - detJ_cc[lev] = std::make_unique(ba,dm,1,1); - if (solverChoice.terrain_type != TerrainType::Static) { + if (solverChoice.terrain_type != TerrainType::Static) + { detJ_cc_new[lev] = std::make_unique(ba,dm,1,1); detJ_cc_src[lev] = std::make_unique(ba,dm,1,1); + + ax_src[lev] = std::make_unique(convert(ba,IntVect(1,0,0)),dm,1,1); + ay_src[lev] = std::make_unique(convert(ba,IntVect(0,1,0)),dm,1,1); + az_src[lev] = std::make_unique(convert(ba,IntVect(0,0,1)),dm,1,1); + + ax_new[lev] = std::make_unique(convert(ba,IntVect(1,0,0)),dm,1,1); + ay_new[lev] = std::make_unique(convert(ba,IntVect(0,1,0)),dm,1,1); + az_new[lev] = std::make_unique(convert(ba,IntVect(0,0,1)),dm,1,1); + z_t_rk[lev] = std::make_unique( convert(ba, IntVect(0,0,1)), dm, 1, 1 ); } @@ -78,7 +87,6 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, } else { z_phys_nd[lev] = nullptr; z_phys_cc[lev] = nullptr; - detJ_cc[lev] = nullptr; z_phys_nd_new[lev] = nullptr; detJ_cc_new[lev] = nullptr; @@ -89,6 +97,17 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, z_t_rk[lev] = nullptr; } + // We use these area arrays regardless of terrain, EB or none of the above + detJ_cc[lev] = std::make_unique(ba,dm,1,1); + ax[lev] = std::make_unique(convert(ba,IntVect(1,0,0)),dm,1,1); + ay[lev] = std::make_unique(convert(ba,IntVect(0,1,0)),dm,1,1); + az[lev] = std::make_unique(convert(ba,IntVect(0,0,1)),dm,1,1); + + detJ_cc[lev]->setVal(1.0); + ax[lev]->setVal(1.0); + ay[lev]->setVal(1.0); + az[lev]->setVal(1.0); + // ******************************************************************************************** // These are the persistent containers for the old and new data // ******************************************************************************************** @@ -319,6 +338,8 @@ ERF::update_terrain_arrays (int lev, Real time) } make_J(geom[lev],*z_phys_nd[lev],*detJ_cc[lev]); + make_areas(geom[lev],*z_phys_nd[lev],*ax[lev],*ay[lev],*az[lev]); + make_zcc(geom[lev],*z_phys_nd[lev],*z_phys_cc[lev]); } } diff --git a/Source/Initialization/ERF_init_from_metgrid.cpp b/Source/Initialization/ERF_init_from_metgrid.cpp index 576d223fd..2a0a6a505 100644 --- a/Source/Initialization/ERF_init_from_metgrid.cpp +++ b/Source/Initialization/ERF_init_from_metgrid.cpp @@ -206,6 +206,7 @@ ERF::init_from_metgrid (int lev) // This makes the Jacobian. make_J(geom[lev],*z_phys, *detJ_cc[lev]); + make_areas(geom[lev],*z_phys,*ax[lev],*ay[lev],*az[lev]); // This defines z at w-cell faces. make_zcc(geom[lev],*z_phys,*z_phys_cc[lev]); diff --git a/Source/Initialization/ERF_init_from_wrfinput.cpp b/Source/Initialization/ERF_init_from_wrfinput.cpp index 14757f8f2..9e59f2adf 100644 --- a/Source/Initialization/ERF_init_from_wrfinput.cpp +++ b/Source/Initialization/ERF_init_from_wrfinput.cpp @@ -216,6 +216,7 @@ ERF::init_from_wrfinput (int lev) } // mf make_J (geom[lev],*z_phys_nd[lev],* detJ_cc[lev]); + make_areas(geom[lev],*z_phys_nd[lev],*ax[lev],*ay[lev],*az[lev]); make_zcc(geom[lev],*z_phys_nd[lev],*z_phys_cc[lev]); } // use_terrain diff --git a/Source/TimeIntegration/ERF_ComputeTimestep.cpp b/Source/TimeIntegration/ERF_ComputeTimestep.cpp index 4bd1c485f..53807cb5b 100644 --- a/Source/TimeIntegration/ERF_ComputeTimestep.cpp +++ b/Source/TimeIntegration/ERF_ComputeTimestep.cpp @@ -62,40 +62,58 @@ ERF::estTimeStep (int level, long& dt_fast_ratio) const average_face_to_cellcenter(ccvel,0, Array{&vars_new[level][Vars::xvel], - &vars_new[level][Vars::yvel], - &vars_new[level][Vars::zvel]}); + &vars_new[level][Vars::yvel], + &vars_new[level][Vars::zvel]}); int l_no_substepping = solverChoice.no_substepping; +#ifdef ERF_USE_EB + EBFArrayBoxFactory ebfact = EBFactory(level); + const MultiFab& detJ = ebfact.getVolFrac(); +#endif + +#ifdef ERF_USE_EB + Real estdt_comp_inv = ReduceMax(S_new, ccvel, detJ, 0, + [=] AMREX_GPU_HOST_DEVICE (Box const& b, + Array4 const& s, + Array4 const& vf, + Array4 const& u) -> Real +#else Real estdt_comp_inv = ReduceMax(S_new, ccvel, 0, [=] AMREX_GPU_HOST_DEVICE (Box const& b, Array4 const& s, Array4 const& u) -> Real +#endif { Real new_comp_dt = -1.e100; amrex::Loop(b, [=,&new_comp_dt] (int i, int j, int k) noexcept { - const Real rho = s(i, j, k, Rho_comp); - const Real rhotheta = s(i, j, k, RhoTheta_comp); - - // NOTE: even when moisture is present, - // we only use the partial pressure of the dry air - // to compute the soundspeed - Real pressure = getPgivenRTh(rhotheta); - Real c = std::sqrt(Gamma * pressure / rho); - - // If we are not doing the acoustic substepping, then the z-direction contributes - // to the computation of the time step - if (l_no_substepping) { - new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0))+c)*dxinv[0]), - ((amrex::Math::abs(u(i,j,k,1))+c)*dxinv[1]), - ((amrex::Math::abs(u(i,j,k,2))+c)*dzinv ), new_comp_dt); - - // If we are doing the acoustic substepping, then the z-direction does not contribute - // to the computation of the time step - } else { - new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0))+c)*dxinv[0]), - ((amrex::Math::abs(u(i,j,k,1))+c)*dxinv[1]), new_comp_dt); +#ifdef ERF_USE_EB + if (vf(i,j,k) > 0.) +#endif + { + const Real rho = s(i, j, k, Rho_comp); + const Real rhotheta = s(i, j, k, RhoTheta_comp); + + // NOTE: even when moisture is present, + // we only use the partial pressure of the dry air + // to compute the soundspeed + Real pressure = getPgivenRTh(rhotheta); + Real c = std::sqrt(Gamma * pressure / rho); + + // If we are not doing the acoustic substepping, then the z-direction contributes + // to the computation of the time step + if (l_no_substepping) { + new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0))+c)*dxinv[0]), + ((amrex::Math::abs(u(i,j,k,1))+c)*dxinv[1]), + ((amrex::Math::abs(u(i,j,k,2))+c)*dzinv ), new_comp_dt); + + // If we are doing the acoustic substepping, then the z-direction does not contribute + // to the computation of the time step + } else { + new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0))+c)*dxinv[0]), + ((amrex::Math::abs(u(i,j,k,1))+c)*dxinv[1]), new_comp_dt); + } } }); return new_comp_dt; diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index 464904398..f0f0c3a99 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -141,7 +141,7 @@ void ERF::advance_dycore(int level, tau12, tau13, tau21, tau23, tau31, tau32, - z_nd, bc_ptr_h, dxInv, + z_nd, detJ_cc[level]->const_array(mfi), bc_ptr_h, dxInv, mf_m, mf_u, mf_v); } else { ComputeStrain_N(bxcc, tbxxy, tbxxz, tbxyz, domain, diff --git a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp index 2bc27505e..413757c60 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp @@ -293,7 +293,7 @@ void erf_slow_rhs_inc (int level, int nrk, s12, s13, s21, s23, s31, s32, - z_nd, bc_ptr_h, dxInv, + z_nd, detJ_cc[level]->const_array(mfi), bc_ptr_h, dxInv, mf_m, mf_u, mf_v); } // profile diff --git a/Source/TimeIntegration/ERF_slow_rhs_post.cpp b/Source/TimeIntegration/ERF_slow_rhs_post.cpp index a0a6a3b2a..7cac3812f 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_post.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_post.cpp @@ -36,6 +36,9 @@ using namespace amrex; * @param[in] most Pointer to MOST class for Monin-Obukhov Similarity Theory boundary condition * @param[in] domain_bcs_type_d device vector for domain boundary conditions * @param[in] z_phys_nd height coordinate at nodes + * @param[in] ax area fractions on x-faces + * @param[in] ay area fractions on y-faces + * @param[in] az area fractions on z-faces * @param[in] detJ Jacobian of the metric transformation at start of time step (= 1 if use_terrain is false) * @param[in] detJ_new Jacobian of the metric transformation at new RK stage time (= 1 if use_terrain is false) * @param[in] mapfac_m map factor at cell centers @@ -67,6 +70,9 @@ void erf_slow_rhs_post (int level, int finest_level, const Gpu::DeviceVector& domain_bcs_type_d, const Vector& domain_bcs_type_h, std::unique_ptr& z_phys_nd, + std::unique_ptr& ax, + std::unique_ptr& ay, + std::unique_ptr& az, std::unique_ptr& detJ, std::unique_ptr& detJ_new, std::unique_ptr& mapfac_m, @@ -222,9 +228,7 @@ void erf_slow_rhs_post (int level, int finest_level, const Array4& mu_turb = l_use_turb ? eddyDiffs->const_array(mfi) : Array4{}; - // Metric terms const Array4& z_nd = l_use_terrain ? z_phys_nd->const_array(mfi) : Array4{}; - const Array4& detJ_arr = l_use_terrain ? detJ->const_array(mfi) : Array4{}; const Array4& detJ_new_arr = l_moving_terrain ? detJ_new->const_array(mfi) : Array4{}; // Map factors @@ -274,8 +278,17 @@ void erf_slow_rhs_post (int level, int finest_level, // ************************************************************************** // Define updates in the RHS of continuity, temperature, and scalar equations // ************************************************************************** + // Metric terms #ifdef ERF_USE_EB - const auto& vf_arr = ebfact.getVolFrac().const_array(mfi); + auto const& ax_arr = ebfact.getAreaFrac()[0]->const_array(mfi); + auto const& ay_arr = ebfact.getAreaFrac()[1]->const_array(mfi); + auto const& az_arr = ebfact.getAreaFrac()[2]->const_array(mfi); + const auto& detJ_arr = ebfact.getVolFrac().const_array(mfi); +#else + auto const& ax_arr = ax->const_array(mfi); + auto const& ay_arr = ay->const_array(mfi); + auto const& az_arr = az->const_array(mfi); + auto const& detJ_arr = detJ->const_array(mfi); #endif AdvType horiz_adv_type, vert_adv_type; @@ -330,15 +343,11 @@ void erf_slow_rhs_post (int level, int finest_level, AdvectionSrcForScalars(tbx, start_comp, num_comp, avg_xmom, avg_ymom, avg_zmom, cur_prim, cell_rhs, -#ifdef ERF_USE_EB - vf_arr, -#else detJ_arr, -#endif dxInv, mf_m, horiz_adv_type, vert_adv_type, horiz_upw_frac, vert_upw_frac, - l_use_terrain, flx_arr, domain, bc_ptr_h); + flx_arr, domain, bc_ptr_h); if (l_use_diff) { @@ -347,7 +356,8 @@ void erf_slow_rhs_post (int level, int finest_level, if (l_use_terrain) { DiffusionSrcForState_T(tbx, domain, start_comp, num_comp, u, v, cur_cons, cur_prim, cell_rhs, - diffflux_x, diffflux_y, diffflux_z, z_nd, detJ_arr, + diffflux_x, diffflux_y, diffflux_z, + z_nd, ax_arr, ay_arr, az_arr, detJ_arr, dxInv, SmnSmn_a, mf_m, mf_u, mf_v, hfx_z, diss, mu_turb, dc, tc, tm_arr, grav_gpu, bc_ptr_d, use_most); diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index 53aefd170..9db14dc2d 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -47,6 +47,9 @@ using namespace amrex; * @param[in] domain_bcs_type_d device vector for domain boundary conditions * @param[in] domain_bcs_type_h host vector for domain boundary conditions * @param[in] z_phys_nd height coordinate at nodes + * @param[in] ax area fractions on x-faces + * @param[in] ay area fractions on y-faces + * @param[in] az area fractions on z-faces * @param[in] detJ Jacobian of the metric transformation (= 1 if use_terrain is false) * @param[in] p0 Reference (hydrostatically stratified) pressure * @param[in] mapfac_m map factor at cell centers @@ -85,7 +88,11 @@ void erf_slow_rhs_pre (int level, int finest_level, std::unique_ptr& most, const Gpu::DeviceVector& domain_bcs_type_d, const Vector& domain_bcs_type_h, - std::unique_ptr& z_phys_nd, std::unique_ptr& detJ, + std::unique_ptr& z_phys_nd, + std::unique_ptr& ax, + std::unique_ptr& ay, + std::unique_ptr& az, + std::unique_ptr& detJ, const MultiFab* p0, std::unique_ptr& mapfac_m, std::unique_ptr& mapfac_u, @@ -102,6 +109,10 @@ void erf_slow_rhs_pre (int level, int finest_level, { BL_PROFILE_REGION("erf_slow_rhs_pre()"); +#ifdef ERF_USE_EB + amrex::ignore_unused(ax,ay,az,detJ); +#endif + const BCRec* bc_ptr_d = domain_bcs_type_d.data(); const BCRec* bc_ptr_h = domain_bcs_type_h.data(); @@ -256,7 +267,7 @@ void erf_slow_rhs_pre (int level, int finest_level, // Terrain metrics const Array4& z_nd = l_use_terrain ? z_phys_nd->const_array(mfi) : Array4{}; - const Array4& detJ_arr = l_use_terrain ? detJ->const_array(mfi) : Array4{}; + const Array4& detJ_arr = detJ->const_array(mfi); //------------------------------------------------------------------------------- // NOTE: Tile boxes with terrain are not intuitive. The linear combination of @@ -356,7 +367,7 @@ void erf_slow_rhs_pre (int level, int finest_level, s12, s13, s21, s23, s31, s32, - z_nd, bc_ptr_h, dxInv, + z_nd, detJ_arr, bc_ptr_h, dxInv, mf_m, mf_u, mf_v); } // profile @@ -400,7 +411,7 @@ void erf_slow_rhs_pre (int level, int finest_level, s12, s13, s21, s23, s31, s32, - er_arr, z_nd, dxInv); + er_arr, z_nd, detJ_arr, dxInv); } else { ComputeStressVarVisc_T(bxcc, tbxxy, tbxxz, tbxyz, mu_eff, mu_turb, cell_data, @@ -408,7 +419,7 @@ void erf_slow_rhs_pre (int level, int finest_level, s12, s13, s21, s23, s31, s32, - er_arr, z_nd, dxInv); + er_arr, z_nd, detJ_arr, dxInv); } // Remove halo cells from tau_ii but extend across valid_box bdry @@ -617,7 +628,6 @@ void erf_slow_rhs_pre (int level, int finest_level, // Terrain metrics const Array4& z_nd = l_use_terrain ? z_phys_nd->const_array(mfi) : Array4{}; - const Array4& detJ_arr = l_use_terrain ? detJ->const_array(mfi) : Array4{}; // Base state const Array4& p0_arr = p0->const_array(mfi); @@ -739,38 +749,29 @@ void erf_slow_rhs_pre (int level, int finest_level, auto const& ax_arr = ebfact.getAreaFrac()[0]->const_array(mfi); auto const& ay_arr = ebfact.getAreaFrac()[1]->const_array(mfi); auto const& az_arr = ebfact.getAreaFrac()[2]->const_array(mfi); - - const auto& vf_arr = ebfact.getVolFrac().const_array(mfi); + const auto& detJ_arr = ebfact.getVolFrac().const_array(mfi); +#else + auto const& ax_arr = ax->const_array(mfi); + auto const& ay_arr = ay->const_array(mfi); + auto const& az_arr = az->const_array(mfi); + auto const& detJ_arr = detJ->const_array(mfi); #endif AdvectionSrcForRho(bx, cell_rhs, rho_u, rho_v, omega_arr, // these are being used to build the fluxes avg_xmom, avg_ymom, avg_zmom, // these are being defined from the fluxes - z_nd, -#ifdef ERF_USE_EB - vf_arr, -#else - detJ_arr, -#endif + ax_arr, ay_arr, az_arr, detJ_arr, dxInv, mf_m, mf_u, mf_v, -#ifdef ERF_USE_EB - ax_arr, ay_arr, az_arr, -#endif - l_use_terrain, flx_arr); + flx_arr); int icomp = RhoTheta_comp; int ncomp = 1; AdvectionSrcForScalars(bx, icomp, ncomp, avg_xmom, avg_ymom, avg_zmom, - cell_prim, cell_rhs, -#ifdef ERF_USE_EB - vf_arr, -#else - detJ_arr, -#endif + cell_prim, cell_rhs, detJ_arr, dxInv, mf_m, l_horiz_adv_type, l_vert_adv_type, l_horiz_upw_frac, l_vert_upw_frac, - l_use_terrain, flx_arr, domain, bc_ptr_h); + flx_arr, domain, bc_ptr_h); if (l_use_diff) { Array4 diffflux_x = dflux_x->array(mfi); @@ -789,7 +790,8 @@ void erf_slow_rhs_pre (int level, int finest_level, if (l_use_terrain) { DiffusionSrcForState_T(bx, domain, n_start, n_comp, u, v, cell_data, cell_prim, cell_rhs, - diffflux_x, diffflux_y, diffflux_z, z_nd, detJ_arr, + diffflux_x, diffflux_y, diffflux_z, + z_nd, ax_arr, ay_arr, az_arr, detJ_arr, dxInv, SmnSmn_a, mf_m, mf_u, mf_v, hfx_z, diss, mu_turb, dc, tc, tm_arr, grav_gpu, bc_ptr_d, use_most); @@ -836,7 +838,8 @@ void erf_slow_rhs_pre (int level, int finest_level, rho_u_rhs, rho_v_rhs, rho_w_rhs, cell_data, u, v, w, rho_u, rho_v, omega_arr, - z_nd, detJ_arr, dxInv, mf_m, mf_u, mf_v, + ax_arr, ay_arr, az_arr, detJ_arr, + dxInv, mf_m, mf_u, mf_v, l_horiz_adv_type, l_vert_adv_type, l_horiz_upw_frac, l_vert_upw_frac, l_use_terrain, lo_z_face, hi_z_face, diff --git a/Source/TimeIntegration/TI_fast_rhs_fun.H b/Source/TimeIntegration/TI_fast_rhs_fun.H index 5ac82f630..ec4414648 100644 --- a/Source/TimeIntegration/TI_fast_rhs_fun.H +++ b/Source/TimeIntegration/TI_fast_rhs_fun.H @@ -55,6 +55,7 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk, prob->init_custom_terrain(fine_geom,*z_phys_nd_new[level],new_substep_time); init_terrain_grid (level,fine_geom,*z_phys_nd_new[level], zlevels_stag); make_J (fine_geom,*z_phys_nd_new[level], *detJ_cc_new[level]); + make_areas (fine_geom,*z_phys_nd_new[level], *ax_new[level], *ay_new[level], *az_new[level]); Real inv_dt = 1./dtau; diff --git a/Source/TimeIntegration/TI_slow_headers.H b/Source/TimeIntegration/TI_slow_headers.H index 031871904..3c8d4cf46 100644 --- a/Source/TimeIntegration/TI_slow_headers.H +++ b/Source/TimeIntegration/TI_slow_headers.H @@ -57,6 +57,9 @@ void erf_slow_rhs_pre (int level, int finest_level, int nrk, const amrex::Gpu::DeviceVector& domain_bcs_type_d, const amrex::Vector& domain_bcs_type, std::unique_ptr& z_phys_nd, + std::unique_ptr& ax, + std::unique_ptr& ay, + std::unique_ptr& az, std::unique_ptr& dJ, const amrex::MultiFab* p0, std::unique_ptr& mapfac_m, @@ -98,6 +101,9 @@ void erf_slow_rhs_post (int level, int finest_level, int nrk, const amrex::Gpu::DeviceVector& domain_bcs_type_d, const amrex::Vector& domain_bcs_type, std::unique_ptr& z_phys_nd, + std::unique_ptr& ax, + std::unique_ptr& ay, + std::unique_ptr& az, std::unique_ptr& dJ_old, std::unique_ptr& dJ_new, std::unique_ptr& mapfac_m, diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index 3f10fbced..978315c11 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -61,16 +61,19 @@ prob->init_custom_terrain(fine_geom,*z_phys_nd[level],old_step_time); init_terrain_grid (level,fine_geom,*z_phys_nd[level], zlevels_stag); make_J (fine_geom,*z_phys_nd[level], *detJ_cc[level]); + make_areas (fine_geom,*z_phys_nd[level], *ax[level], *ay[level], *az[level]); if (verbose) Print() << "Making src geometry at old_stage_time: " << old_stage_time << std::endl; prob->init_custom_terrain(fine_geom,*z_phys_nd_src[level],old_stage_time); init_terrain_grid (level,fine_geom,*z_phys_nd_src[level], zlevels_stag); make_J (fine_geom,*z_phys_nd_src[level], *detJ_cc_src[level]); + make_areas (fine_geom,*z_phys_nd_src[level], *ax_src[level], *ay_src[level], *az_src[level]); if (verbose) Print() << "Making new geometry at new_stage_time: " << new_stage_time << std::endl; prob->init_custom_terrain(fine_geom,*z_phys_nd_new[level],new_stage_time); init_terrain_grid (level,fine_geom,*z_phys_nd_new[level], zlevels_stag); make_J (fine_geom,*z_phys_nd_new[level], *detJ_cc_new[level]); + make_areas (fine_geom,*z_phys_nd_new[level], *ax_new[level], *ay_new[level], *az_new[level]); Real inv_dt = 1./slow_dt; @@ -114,7 +117,7 @@ Tau13_lev[level].get(), Tau21_lev[level].get(), Tau23_lev[level].get(), Tau31_lev[level].get(), Tau32_lev[level].get(), SmnSmn, eddyDiffs, Hfx3, Diss, fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type, - z_phys_nd_src[level], detJ_cc_src[level], p0_new, + z_phys_nd_src[level], ax_src[level], ay_src[level], az_src[level], detJ_cc_src[level], p0_new, mapfac_m[level], mapfac_u[level], mapfac_v[level], #ifdef ERF_USE_EB EBFactory(level), @@ -206,7 +209,7 @@ Tau13_lev[level].get(), Tau21_lev[level].get(), Tau23_lev[level].get(), Tau31_lev[level].get(), Tau32_lev[level].get(), SmnSmn, eddyDiffs, Hfx3, Diss, fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type, - z_phys_nd[level], detJ_cc[level], p0, + z_phys_nd[level], ax[level], ay[level], az[level], detJ_cc[level], p0, mapfac_m[level], mapfac_u[level], mapfac_v[level], #ifdef ERF_USE_EB EBFactory(level), @@ -315,7 +318,7 @@ source, SmnSmn, eddyDiffs, Hfx3, Diss, fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type, - z_phys_nd_src[level], detJ_cc[level], detJ_cc_new[level], + z_phys_nd[level], ax[level], ay[level], az[level], detJ_cc[level], detJ_cc_new[level], mapfac_m[level], mapfac_u[level], mapfac_v[level], #ifdef ERF_USE_EB EBFactory(level), @@ -333,7 +336,7 @@ source, SmnSmn, eddyDiffs, Hfx3, Diss, fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type, - z_phys_nd[level], detJ_cc[level], detJ_cc[level], + z_phys_nd[level], ax[level], ay[level], az[level], detJ_cc[level], detJ_cc[level], mapfac_m[level], mapfac_u[level], mapfac_v[level], #ifdef ERF_USE_EB EBFactory(level), diff --git a/Source/Utils/TerrainMetrics.cpp b/Source/Utils/TerrainMetrics.cpp index b17dfd189..c57fdaaf0 100644 --- a/Source/Utils/TerrainMetrics.cpp +++ b/Source/Utils/TerrainMetrics.cpp @@ -484,6 +484,71 @@ make_J (const Geometry& geom, detJ_cc.FillBoundary(geom.periodicity()); } +/** + * Computation of area fractions on faces + */ +void +make_areas (const Geometry& geom, + MultiFab& z_phys_nd, MultiFab& ax, + MultiFab& ay, MultiFab& az) +{ + const auto* dx = geom.CellSize(); + Real dzInv = 1.0/dx[2]; + + // Domain valid box (z_nd is nodal) + const Box& domain = geom.Domain(); + int domlo_z = domain.smallEnd(2); + + // The z-faces are always full when using terrain-fitted coordinates + az.setVal(1.0); + + // + // x-areas + // +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(ax, TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + Box gbx = mfi.growntilebox(ax.nGrow()); + if (gbx.smallEnd(2) < domlo_z) { + gbx.setSmall(2,domlo_z); + } + + Array4 z_nd = z_phys_nd.const_array(mfi); + Array4 ax_arr = ax.array(mfi); + ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + ax_arr(i, j, k) = .5 * dzInv * ( + z_nd(i,j,k+1) + z_nd(i,j+1,k+1) - z_nd(i,j,k) - z_nd(i,j+1,k)); + }); + } + + // + // y-areas + // +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(ay, TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + Box gbx = mfi.growntilebox(ay.nGrow()); + if (gbx.smallEnd(2) < domlo_z) { + gbx.setSmall(2,domlo_z); + } + + Array4 z_nd = z_phys_nd.const_array(mfi); + Array4 ay_arr = ay.array(mfi); + ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + ay_arr(i, j, k) = .5 * dzInv * ( + z_nd(i,j,k+1) + z_nd(i+1,j,k+1) - z_nd(i,j,k) - z_nd(i+1,j,k)); + }); + } + + ax.FillBoundary(geom.periodicity()); + ay.FillBoundary(geom.periodicity()); + az.FillBoundary(geom.periodicity()); +} + /** * Computation of z_phys at cell-center */ diff --git a/Source/Utils/Utils.H b/Source/Utils/Utils.H index d25773f51..282b77ddc 100644 --- a/Source/Utils/Utils.H +++ b/Source/Utils/Utils.H @@ -16,6 +16,12 @@ void make_J (const amrex::Geometry& geom, amrex::MultiFab& z_phys_nd, amrex::MultiFab& detJ_cc); +void make_areas (const amrex::Geometry& geom, + amrex::MultiFab& z_phys_nd, + amrex::MultiFab& ax, + amrex::MultiFab& ay, + amrex::MultiFab& az); + /* * Average z_phys_nd on nodes to cell centers */