From 99a66902f02f674bdf7ecbcf1361292d71e17d38 Mon Sep 17 00:00:00 2001 From: Soonpil Kang <109235650+skang67@users.noreply.github.com> Date: Thu, 12 Dec 2024 07:28:56 -0800 Subject: [PATCH] Setup EB Scalar DevTest. (#2011) * Eliminated ERF_USE_EB preprocessor in AdvectionSrcForMom, because it directed to use_terrain part. * Added a scalar bubble to EB_Test. * Added ERF_redistribute to ERF class. * Fixed small things in EB_Test. * Remove whitespaces. --------- Co-authored-by: Soonpil Kang --- Exec/DevTests/EB_Test/ERF_Prob.H | 3 + Exec/DevTests/EB_Test/ERF_Prob.cpp | 66 +++++++++++++++++-- Exec/DevTests/EB_Test/inputs | 19 ++++-- ..._redistribute.cpp => ERF_Redistribute.cpp} | 36 +++++----- Source/EB/Make.package | 1 - Source/ERF.H | 12 ++++ .../TimeIntegration/ERF_TI_no_substep_fun.H | 8 ++- 7 files changed, 115 insertions(+), 30 deletions(-) rename Source/EB/{ERF_redistribute.cpp => ERF_Redistribute.cpp} (77%) diff --git a/Exec/DevTests/EB_Test/ERF_Prob.H b/Exec/DevTests/EB_Test/ERF_Prob.H index 8d670810a..9dc736590 100644 --- a/Exec/DevTests/EB_Test/ERF_Prob.H +++ b/Exec/DevTests/EB_Test/ERF_Prob.H @@ -23,6 +23,9 @@ struct ProbParm : ProbParmDefaults { amrex::Real yc_frac = 0.5; // Location of "center" of scalar (multiplies domain length) amrex::Real zc_frac = 0.5; // Location of "center" of scalar (multiplies domain length) + amrex::Real xradius = 10.0; // x-radius of scalar bubble + amrex::Real zradius = 10.0; // z-radius of scalar bubble + int prob_type = -1; }; // namespace ProbParm diff --git a/Exec/DevTests/EB_Test/ERF_Prob.cpp b/Exec/DevTests/EB_Test/ERF_Prob.cpp index 1c4ec28cc..e15bd4222 100644 --- a/Exec/DevTests/EB_Test/ERF_Prob.cpp +++ b/Exec/DevTests/EB_Test/ERF_Prob.cpp @@ -32,6 +32,9 @@ Problem::Problem() pp.query("prob_type", parms.prob_type); + pp.query("xradius", parms.xradius); + pp.query("zradius", parms.zradius); + init_base_parms(parms.rho_0, parms.T_0); } @@ -81,16 +84,71 @@ Problem::init_custom_pert( AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + // ParallelFor(bx, [=] 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; + // } + // }); + + // Set the state_pert + ParallelFor(bx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - // Set scalar = 0 everywhere - state_pert(i, j, k, RhoScalar_comp) = 0.0; + // 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_d.xc_frac * (prob_lo[0] + prob_hi[0]); + const Real yc = parms_d.yc_frac * (prob_lo[1] + prob_hi[1]); + const Real zc = parms_d.zc_frac * (prob_lo[2] + prob_hi[2]); + + // Define ellipse parameters + const Real r0 = parms_d.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)); + const Real r2d_xz_nd = std::sqrt((x-xc)*(x-xc)/parms_d.xradius/parms_d.xradius + + (z-zc)*(z-zc)/parms_d.zradius/parms_d.zradius); + + if (parms_d.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_d.A_0 * exp(-10.*r3d*r3d) + parms_d.B_0*sin(x); + state_pert(i, j, k, RhoScalar_comp) = parms_d.A_0 * exp(-0.1*r2d_xz*r2d_xz) + parms_d.B_0*sin(x); + + } else if (parms_d.prob_type == 11) { + if (r2d_xz_nd < 1.0) + { + state_pert(i, j, k, RhoScalar_comp) = 0.5 * parms_d.A_0 * (1.0 + std::cos(PI*r2d_xz_nd)); + } else { + state_pert(i, j, k, RhoScalar_comp) = 0.0; + } + } 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_d.A_0; + } else { + state_pert(i, j, k, RhoScalar_comp) = 0.0; + } + } + + state_pert(i, j, k, RhoScalar_comp) *= parms_d.rho_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 ParallelFor(xbx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept diff --git a/Exec/DevTests/EB_Test/inputs b/Exec/DevTests/EB_Test/inputs index 6cf5daf48..75e3e23f7 100644 --- a/Exec/DevTests/EB_Test/inputs +++ b/Exec/DevTests/EB_Test/inputs @@ -1,6 +1,5 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -max_step = 100 -max_step = 10 +max_step = 800 amr.max_grid_size = 256 256 256 @@ -17,16 +16,16 @@ fabarray.mfiter_tile_size = 1024 1024 1024 # PROBLEM SIZE & GEOMETRY geometry.prob_lo = 0. 0. 0. -geometry.prob_hi = 16. 1. 16. +geometry.prob_hi = 16. 4. 8. -amr.n_cell = 64 4 64 +amr.n_cell = 64 4 32 geometry.is_periodic = 0 1 0 xlo.type = Inflow xhi.type = Outflow -xlo.velocity = 1. 0. 0. +xlo.velocity = 10. 0. 0. xlo.density = 1.0 xlo.theta = 1.0 xlo.scalar = 0. @@ -36,7 +35,7 @@ zhi.type = SlipWall # TIME STEP CONTROL erf.substepping_type = None -erf.fixed_dt = 1.e-5 +erf.fixed_dt = 1.e-3 # DIAGNOSTICS & VERBOSITY erf.sum_interval = 1 # timesteps between computing mass @@ -72,4 +71,10 @@ erf.init_type = "uniform" # PROBLEM PARAMETERS prob.rho_0 = 1.0 prob.T_0 = 1.0 -prob.u_0 = 1.0 +prob.u_0 = 10.0 +prob.A_0 = 1.0 +prob.prob_type = 11 +prob.xradius = 3.0 +prob.zradius = 1.5 +prob.xc_frac = 0.25 +prob.zc_frac = 0.15 diff --git a/Source/EB/ERF_redistribute.cpp b/Source/EB/ERF_Redistribute.cpp similarity index 77% rename from Source/EB/ERF_redistribute.cpp rename to Source/EB/ERF_Redistribute.cpp index 184d73ad1..aa98941c1 100644 --- a/Source/EB/ERF_redistribute.cpp +++ b/Source/EB/ERF_Redistribute.cpp @@ -9,10 +9,11 @@ using namespace amrex; void ERF::redistribute_term ( int lev, - MultiFab& result, - MultiFab& result_tmp, // Saves doing a MF::copy. does this matter??? - MultiFab const& state, - BCRec const* bc) // this is bc for the state (needed for SRD slopes) + MultiFab& result, + MultiFab& result_tmp, // Saves doing a MF::copy. does this matter??? + MultiFab const& state, + BCRec const* bc, // this is bc for the state (needed for SRD slopes) + Real const dt) { // ************************************************************************ // Redistribute result_tmp and pass out result @@ -26,16 +27,17 @@ ERF::redistribute_term ( int lev, #endif for (MFIter mfi(state,TilingIfNotGPU()); mfi.isValid(); ++mfi) { - redistribute_term(mfi, result, result_tmp, state, bc, lev); + redistribute_term(mfi, lev, result, result_tmp, state, bc, dt); } } void ERF::redistribute_term ( MFIter const& mfi, int lev, - MultiFab& result, - MultiFab& result_tmp, - MultiFab const& state, - BCRec const* bc) // this is bc for the state (needed for SRD slopes) + MultiFab& result, + MultiFab& result_tmp, + MultiFab const& state, + BCRec const* bc, // this is bc for the state (needed for SRD slopes) + Real const dt) { AMREX_ASSERT(result.nComp() == state.nComp()); @@ -57,13 +59,13 @@ ERF::redistribute_term ( MFIter const& mfi, int lev, auto const& vfrac = ebfact.getVolFrac().const_array(mfi); auto const& ccc = ebfact.getCentroid().const_array(mfi); - auto const& apx = ebfact.getAreaFrac()[0]->const_array(mfi);, - auto const& apy = ebfact.getAreaFrac()[1]->const_array(mfi);, - auto const& apz = ebfact.getAreaFrac()[2]->const_array(mfi);); + auto const& apx = ebfact.getAreaFrac()[0]->const_array(mfi); + auto const& apy = ebfact.getAreaFrac()[1]->const_array(mfi); + auto const& apz = ebfact.getAreaFrac()[2]->const_array(mfi); - auto const& fcx = ebfact.getFaceCent()[0]->const_array(mfi);, - auto const& fcy = ebfact.getFaceCent()[1]->const_array(mfi);, - auto const& fcz = ebfact.getFaceCent()[2]->const_array(mfi);); + auto const& fcx = ebfact.getFaceCent()[0]->const_array(mfi); + auto const& fcy = ebfact.getFaceCent()[1]->const_array(mfi); + auto const& fcz = ebfact.getFaceCent()[2]->const_array(mfi); Box gbx = bx; gbx.grow(3); @@ -77,7 +79,7 @@ ERF::redistribute_term ( MFIter const& mfi, int lev, // scratch(i,j,k) = 1.; // }); - std::string redistribution_type = "StateRedistribution"; + std::string redistribution_type = "StateRedist"; // State redist acts on the state. Array4 state_arr = state.const_array(mfi); @@ -85,7 +87,7 @@ ERF::redistribute_term ( MFIter const& mfi, int lev, scratch, flag, apx, apy, apz, vfrac, fcx, fcy, fcz, ccc, - bc, geom[lev], m_dt, edistribution_type); + bc, geom[lev], dt, redistribution_type); } else { diff --git a/Source/EB/Make.package b/Source/EB/Make.package index a9e8662af..8bae36093 100644 --- a/Source/EB/Make.package +++ b/Source/EB/Make.package @@ -6,6 +6,5 @@ CEXE_sources += ERF_EBCylinder.cpp CEXE_sources += ERF_EBRegular.cpp CEXE_sources += ERF_Redistribute.cpp CEXE_sources += ERF_WriteEBSurface.cpp - CEXE_headers += ERF_EBIF.H CEXE_headers += ERF_TerrainIF.H diff --git a/Source/ERF.H b/Source/ERF.H index d49fffe32..89163e184 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -475,6 +475,18 @@ public: void make_eb_box (); void make_eb_cylinder (); void make_eb_regular (); + void redistribute_term ( int lev, + amrex::MultiFab& result, + amrex::MultiFab& result_tmp, + amrex::MultiFab const& state, + amrex::BCRec const* bc, + amrex::Real const dt); + void redistribute_term ( amrex::MFIter const& mfi, int lev, + amrex::MultiFab& result, + amrex::MultiFab& result_tmp, + amrex::MultiFab const& state, + amrex::BCRec const* bc, + amrex::Real const dt); #endif // more flexible version of AverageDown() that lets you average down across multiple levels diff --git a/Source/TimeIntegration/ERF_TI_no_substep_fun.H b/Source/TimeIntegration/ERF_TI_no_substep_fun.H index 13f6e2e3f..074c0c2c7 100644 --- a/Source/TimeIntegration/ERF_TI_no_substep_fun.H +++ b/Source/TimeIntegration/ERF_TI_no_substep_fun.H @@ -151,7 +151,13 @@ } } // mfi #ifdef ERF_USE_EB - // call to redistribute_term -- pass in ssum[IntVars::cons] which is MF + // // call to redistribute_term -- pass in ssum[IntVars::cons] which is MF + // MultiFab dUdt_tmp (ba, dm, n_data, 3, MFInfo(), Factory(level)); + // dUdt_tmp.setVal(0.); + // dUdt_tmp.FillBoundary(fine_geom.periodicity()); + // const BCRec* bc_ptr_h = domain_bcs_type.data(); // Should be host or device or both? + // redistribute_term ( level, F_slow[IntVars::cons], dUdt_tmp, + // S_sum[IntVars::cons], bc_ptr_h, slow_dt); #endif } // omp