From c276af9fe6b7d56e25c375d27db02d2915819184 Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Tue, 22 Aug 2023 19:55:44 -0700 Subject: [PATCH] Streamline FPr BA (#1199) * Fixes for use with amrex patch for face linear interp. * not sure why this is needed with substepping. * This requires amrex patch 3483. Also, do relaxation on all RK stages; from tests with 2d cylinder. --------- Co-authored-by: Ann Almgren --- .../DensityCurrent/inputs_periodic_twolevel | 68 ++++++++++++++ Exec/RegTests/Terrain3d_Hemisphere/inputs.amr | 94 +++++++++++++++++++ Source/BoundaryConditions/ERF_FillPatcher.H | 2 +- Source/BoundaryConditions/ERF_FillPatcher.cpp | 36 ++----- Source/TimeIntegration/TI_slow_rhs_fun.H | 8 +- Source/Utils/InteriorGhostCells.cpp | 3 + 6 files changed, 179 insertions(+), 32 deletions(-) create mode 100644 Exec/RegTests/DensityCurrent/inputs_periodic_twolevel create mode 100644 Exec/RegTests/Terrain3d_Hemisphere/inputs.amr diff --git a/Exec/RegTests/DensityCurrent/inputs_periodic_twolevel b/Exec/RegTests/DensityCurrent/inputs_periodic_twolevel new file mode 100644 index 000000000..ec107d2d0 --- /dev/null +++ b/Exec/RegTests/DensityCurrent/inputs_periodic_twolevel @@ -0,0 +1,68 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 999999 +stop_time = 900.0 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_lo = -25600. 0. 0. +geometry.prob_hi = 25600. 400. 6400. + +amr.n_cell = 512 12 64 # dx=dy=dz=100 m, Straka et al 1993 / Xue et al 2000 + +# periodic in x to match WRF setup +# - as an alternative, could use symmetry at x=0 and outflow at x=25600 +geometry.is_periodic = 1 1 0 +zlo.type = "SlipWall" +zhi.type = "SlipWall" + +# TIME STEP CONTROL +erf.fixed_dt = 0.1 # fixed time step [s] -- Straka et al 1993 +erf.fixed_fast_dt = 0.025 # fixed time step [s] -- Straka et al 1993 + +# DIAGNOSTICS & VERBOSITY +erf.sum_interval = 1 # timesteps between computing mass +erf.v = 1 # verbosity in ERF.cpp +amr.v = 1 # verbosity in Amr.cpp + +# MULTILEVEL FLAGS +amr.max_level = 1 +amr.ref_ratio_vect = 2 2 1 +erf.refinement_indicators = box1 +erf.box1.max_level = 1 +erf.box1.in_box_lo = -6400. 0. +erf.box1.in_box_hi = 6400. 400. +erf.coupling_type = "OneWay" +erf.cf_width = 5 +erf.cf_set_width = 1 + +# CHECKPOINT FILES +erf.check_file = chk # root name of checkpoint file +erf.check_int = -57600 # number of timesteps between checkpoints + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 1 # number of timesteps between plotfiles +erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta pres_hse dens_hse pert_pres pert_dens + +# SOLVER CHOICE +erf.use_gravity = true +erf.use_coriolis = false +erf.use_rayleigh_damping = false + +# +# diffusion coefficient from Straka, K = 75 m^2/s +# +erf.les_type = "None" +erf.molec_diff_type = "ConstantAlpha" +erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities +erf.dynamicViscosity = 75.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s +erf.alpha_T = 75.0 # [m^2/s] + +erf.c_p = 1004.5 + +# PROBLEM PARAMETERS (optional) +prob.T_0 = 300.0 +prob.U_0 = 10.0 \ No newline at end of file diff --git a/Exec/RegTests/Terrain3d_Hemisphere/inputs.amr b/Exec/RegTests/Terrain3d_Hemisphere/inputs.amr new file mode 100644 index 000000000..a10dd8858 --- /dev/null +++ b/Exec/RegTests/Terrain3d_Hemisphere/inputs.amr @@ -0,0 +1,94 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 4000 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_lo = 0. 0. 0. +geometry.prob_hi = 10 10 10 + +amr.n_cell = 100 100 100 # dx=dy=dz=100 m, Straka et al 1993 / Xue et al 2000 + +geometry.is_periodic = 0 1 0 + +xlo.type = "Outflow" +xhi.type = "Outflow" +xlo.velocity = 1. 0. 0. +xlo.density = 1.16 +xlo.theta = 300. +xlo.scalar = 0. + +zlo.type = "SlipWall" +zhi.type = "Outflow" + +erf.sponge_strength = 100.0 +erf.use_xlo_sponge_damping = true +erf.xlo_sponge_end = 2.0 +erf.use_xhi_sponge_damping = true +erf.xhi_sponge_start = 8.0 +erf.use_ylo_sponge_damping = true +erf.ylo_sponge_end = 2.0 +erf.use_yhi_sponge_damping = true +erf.yhi_sponge_start = 8.0 +erf.use_zhi_sponge_damping = true +erf.zhi_sponge_start = 8.0 + +erf.sponge_density = 1.2 +erf.sponge_x_velocity = 10.0 +erf.sponge_y_velocity = 0.0 +erf.sponge_z_velocity = 0.0 + + +# TIME STEP CONTROL +erf.no_substepping = 1 +#erf.fixed_dt = 1E-5 + +erf.cfl = 0.3 + +# DIAGNOSTICS & VERBOSITY +erf.sum_interval = 1 # timesteps between computing mass +erf.v = 1 # verbosity in ERF.cpp +amr.v = 1 # verbosity in Amr.cpp + +# REFINEMENT / REGRIDDING +amr.max_level = 1 +amr.ref_ratio_vect = 2 2 1 +erf.refinement_indicators = box1 +erf.box1.max_level = 1 +erf.box1.in_box_lo = 3. 3. +erf.box1.in_box_hi = 7. 7. +erf.coupling_type = "OneWay" +erf.cf_width = 0 #5 +erf.cf_set_width = 0 #1 + +# CHECKPOINT FILES +erf.check_file = chk # root name of checkpoint file +erf.check_int = -57600 # number of timesteps between checkpoints + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 100 # number of timesteps between plotfiles +erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta pres_hse dens_hse pert_pres pert_dens z_phys detJ dpdx dpdy pres_hse_x pres_hse_y + +# SOLVER CHOICE +erf.use_gravity = true +erf.use_coriolis = false +erf.use_rayleigh_damping = false +erf.les_type = "None" + +# TERRRAIN GRID TYPE +erf.use_terrain = true +erf.terrain_smoothing = 0 + +# Diffusion coefficient from Straka, K = 75 m^2/s +erf.molec_diff_type = "ConstantAlpha" +erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities +erf.dynamicViscosity = 0.1 # [kg/(m-s)] ==> nu = 75.0 m^2/s +erf.alpha_T = 0.0 # [m^2/s] + +# PROBLEM PARAMETERS (optional) +prob.T_0 = 300.0 +prob.U_0 = 0.0 +prob.rho_0 = 1.16 diff --git a/Source/BoundaryConditions/ERF_FillPatcher.H b/Source/BoundaryConditions/ERF_FillPatcher.H index 48ec6e002..b52afe349 100644 --- a/Source/BoundaryConditions/ERF_FillPatcher.H +++ b/Source/BoundaryConditions/ERF_FillPatcher.H @@ -67,7 +67,7 @@ void ERFFillPatcher::fill (amrex::MultiFab& mf, amrex::Real time, BC& cbc, amrex::Vector const& bcs, bool fill_subset) { - constexpr amrex::Real eps = std::numeric_limits::epsilon(); + constexpr amrex::Real eps = std::numeric_limits::epsilon(); AMREX_ALWAYS_ASSERT((time >= m_crse_times[0]-eps) && (time <= m_crse_times[1]+eps)); // Time interpolation factors diff --git a/Source/BoundaryConditions/ERF_FillPatcher.cpp b/Source/BoundaryConditions/ERF_FillPatcher.cpp index e028ffc32..381e68778 100644 --- a/Source/BoundaryConditions/ERF_FillPatcher.cpp +++ b/Source/BoundaryConditions/ERF_FillPatcher.cpp @@ -116,49 +116,31 @@ void ERFFillPatcher::Define (BoxArray const& fba, DistributionMapping fdm, } } - // Grown fine box list (expand interior halo to touch crse face) - BoxList gbl; - gbl.set(m_ixt); - for (auto bl_iter = bl.begin(); bl_iter < bl.end(); bl_iter++) { - const int* lo = (*bl_iter).loVect(); - const int* hi = (*bl_iter).hiVect(); - IntVect glo(0); IntVect ghi(0); - for (int idim(0); idim < AMREX_SPACEDIM; ++idim) { - int rlo = lo[idim]%m_ratio[idim]; - int rhi = hi[idim]%m_ratio[idim]; - glo[idim] = lo[idim]; ghi[idim] = hi[idim]; - if (rlo > 0) glo[idim] = lo[idim] - rlo; - if (rhi > 0) ghi[idim] = hi[idim] - rhi + m_ratio[idim]; - } - Box gbx(glo,ghi,m_ixt); - gbl.push_back(gbx); - } - // Coarse box list (coincides with grown fine box list) BoxList cbl; cbl.set(m_ixt); - cbl.reserve(gbl.size()); - for (auto const& b : gbl) { + cbl.reserve(bl.size()); + for (auto const& b : bl) { cbl.push_back(interp->CoarseBox(b, m_ratio)); } // Box arrays for fine and coarse - BoxArray gcf_fba(std::move(gbl)); - BoxArray gcf_cba(std::move(cbl)); + BoxArray cf_fba(std::move(bl)); + BoxArray cf_cba(std::move(cbl)); BoxArray cf_fba_s(std::move(bl_subset)); - DistributionMapping gcf_dm(gcf_fba); + DistributionMapping cf_dm(cf_fba); DistributionMapping cf_dm_s(cf_fba_s); // Fine patch to hold the time-interpolated state - m_cf_fine_data = new MultiFab (gcf_fba, gcf_dm, m_ncomp, 0); + m_cf_fine_data = new MultiFab (cf_fba, cf_dm, m_ncomp, 0); // Fine subset patch to hold the time-interpolated state - m_cf_fine_subset_data = new MultiFab (cf_fba_s, gcf_dm, m_ncomp, 0); + m_cf_fine_subset_data = new MultiFab (cf_fba_s, cf_dm, m_ncomp, 0); // Two coarse patches to hold the data to be interpolated - m_cf_crse_data[0] = new MultiFab (gcf_cba, gcf_dm, m_ncomp, 0); - m_cf_crse_data[1] = new MultiFab (gcf_cba, gcf_dm, m_ncomp, 0); + m_cf_crse_data[0] = new MultiFab (cf_cba, cf_dm, m_ncomp, 0); + m_cf_crse_data[1] = new MultiFab (cf_cba, cf_dm, m_ncomp, 0); } diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index be8a41e3e..9f1661199 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -194,7 +194,7 @@ #ifdef ERF_USE_NETCDF // Populate RHS for relaxation zones - if (init_type=="real" && level==0 && nrk==0) { + if (init_type=="real" && level==0) { wrfbdy_compute_interior_ghost_RHS(bdy_time_interval, start_bdy_time, new_stage_time, slow_dt, wrfbdy_width-1, wrfbdy_set_width, fine_geom, S_rhs, S_data, bdy_data_xlo, bdy_data_xhi, @@ -203,7 +203,7 @@ #endif // Compute RHS for fine interior ghost - if (level>0 && coupling_type=="OneWay" && cf_width>0 && nrk==0) { + if (level>0 && coupling_type=="OneWay" && cf_width>0) { fine_compute_interior_ghost_RHS(new_stage_time, slow_dt, cf_width-1, cf_set_width, fine_geom, &FPr_c[level-1], &FPr_u[level-1], &FPr_v[level-1], &FPr_w[level-1], @@ -342,7 +342,7 @@ #ifdef ERF_USE_NETCDF // Populate RHS for relaxation zones - if (init_type=="real" && level==0 && nrk==0) { + if (init_type=="real" && level==0) { wrfbdy_compute_interior_ghost_RHS(bdy_time_interval, start_bdy_time, new_stage_time, slow_dt, wrfbdy_width-1, wrfbdy_set_width, fine_geom, S_rhs, S_data, bdy_data_xlo, bdy_data_xhi, @@ -351,7 +351,7 @@ #endif // Compute RHS for fine interior ghost - if (level>0 && coupling_type=="OneWay" && cf_width>0 && nrk==0) { + if (level>0 && coupling_type=="OneWay" && cf_width>0) { fine_compute_interior_ghost_RHS(new_stage_time, slow_dt, cf_width-1, cf_set_width, &FPr_c[level-1], &FPr_u[level-1], &FPr_v[level-1], &FPr_w[level-1], diff --git a/Source/Utils/InteriorGhostCells.cpp b/Source/Utils/InteriorGhostCells.cpp index 138a6ca6e..902b81844 100644 --- a/Source/Utils/InteriorGhostCells.cpp +++ b/Source/Utils/InteriorGhostCells.cpp @@ -577,6 +577,7 @@ fine_compute_interior_ghost_RHS(const Real& time, Box domain = boxes_at_level[g_ind]; domain.convert(fmf.boxArray().ixType()); + // NOTE:: width+1 to get halo interior cell from fill Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; compute_interior_ghost_bxs_xy(tbx, domain, width+1, 0, bx_xlo, bx_xhi, @@ -625,6 +626,7 @@ fine_compute_interior_ghost_RHS(const Real& time, Box domain = boxes_at_level[g_ind]; domain.convert(fmf.boxArray().ixType()); + // NOTE:: width+1 to get halo interior cell from fill Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; compute_interior_ghost_bxs_xy(tbx, domain, width+1, 0, bx_xlo, bx_xhi, @@ -672,6 +674,7 @@ fine_compute_interior_ghost_RHS(const Real& time, Box domain = boxes_at_level[g_ind]; domain.convert(fmf.boxArray().ixType()); + // NOTE:: width+1 to get halo interior cell from fill Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; compute_interior_ghost_bxs_xy(tbx, domain, width+1, 0, bx_xlo, bx_xhi,