From a6029a999d97f5173b1826c3995d2b0a2417e7a8 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 27 Nov 2023 17:49:59 -0800 Subject: [PATCH] fix for different regridding options --- Exec/inputs.mesh_refinement | 3 +- Source/BoundaryConditions/ERF_FillPatch.cpp | 2 +- Source/ERF.H | 2 +- Source/ERF.cpp | 43 ++++--- Source/ERF_make_new_level.cpp | 122 ++++++++++++++++++-- Source/IO/Plotfile.cpp | 9 +- Source/TimeIntegration/ERF_TimeStep.cpp | 35 ++---- Source/TimeIntegration/TI_slow_rhs_fun.H | 6 +- 8 files changed, 161 insertions(+), 61 deletions(-) diff --git a/Exec/inputs.mesh_refinement b/Exec/inputs.mesh_refinement index 9e5406bb9..55da0226c 100644 --- a/Exec/inputs.mesh_refinement +++ b/Exec/inputs.mesh_refinement @@ -5,9 +5,10 @@ # REFINEMENT / REGRIDDING amr.max_level = 1 # maximum level number allowed amr.ref_ratio_vect = 2 2 1 3 3 1 # refinement ratio written as lev0_x lev0_y lev0_z lev1_x lev1_y lev1_z ... -amr.regrid_int = 2 2 2 2 # how often to regrid amr.n_error_buf = 2 2 2 2 # number of buffer cells in error est +erf.regrid_int = 2 2 2 2 # how often to regrid + # How to specify a region for static refinement erf.refinement_indicators = box1 erf.box1.in_box_lo = .15 .25 diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp index d4c863c1e..9b1b6dd37 100644 --- a/Source/BoundaryConditions/ERF_FillPatch.cpp +++ b/Source/BoundaryConditions/ERF_FillPatch.cpp @@ -73,7 +73,7 @@ ERF::FillPatch (int lev, Real time, const Vector& mfs, bool fillset) } // var_idx // Coarse-Fine set region - if (lev>0 && solverChoice.coupling_type == CouplingType::OneWay && cf_set_width>0 && fillset) { + if (lev>0 && cf_set_width>0 && fillset) { FPr_c[lev-1].FillSet(*mfs[Vars::cons], time, null_bc, domain_bcs_type); FPr_u[lev-1].FillSet(*mfs[Vars::xvel], time, null_bc, domain_bcs_type); FPr_v[lev-1].FillSet(*mfs[Vars::yvel], time, null_bc, domain_bcs_type); diff --git a/Source/ERF.H b/Source/ERF.H index e94eaa87d..10020f20a 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -336,7 +336,7 @@ private: // set covered coarse cells to be the average of overlying fine cells void AverageDown (int scomp, int ncomp); - void update_arrays (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm); + void update_diffusive_arrays (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm); void update_terrain_arrays (int lev, amrex::Real time); diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 3c0a73467..ff4a2281f 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -681,21 +681,24 @@ ERF::InitData () // Fill ghost cells/faces for (int lev = 0; lev <= finest_level; ++lev) { - // - // NOTE: we must set up the FillPatcher object before calling FillPatch at a fine level - // - if (solverChoice.coupling_type != CouplingType::TwoWay && cf_width>0 && lev>0) { + if (lev > 0 && cf_width > 0) { Construct_ERFFillPatchers(lev); - Register_ERFFillPatchers(lev); } - auto& lev_new = vars_new[lev]; + // + // We don't use the FillPatcher in this call because + // we don't need to fill the interior data at this point. + // + bool fillset = false; + auto& lev_new = vars_new[lev]; FillPatch(lev, t_new[lev], - {&lev_new[Vars::cons],&lev_new[Vars::xvel],&lev_new[Vars::yvel],&lev_new[Vars::zvel]}); + {&lev_new[Vars::cons],&lev_new[Vars::xvel],&lev_new[Vars::yvel],&lev_new[Vars::zvel]}, + fillset); - // We need to fill the ghost cell values of the base state in case it wasn't - // done in the initialization + // + // We fill the ghost cell values of the base state in case it wasn't done in the initialization + // base_state[lev].FillBoundary(geom[lev].periodicity()); // For moving terrain only @@ -1102,9 +1105,9 @@ ERF::ReadParameters () // Query the set and total widths for crse-fine interior ghost cells pp.query("cf_width", cf_width); pp.query("cf_set_width", cf_set_width); - AMREX_ALWAYS_ASSERT(cf_width >= 0); - AMREX_ALWAYS_ASSERT(cf_set_width >= 0); - AMREX_ALWAYS_ASSERT(cf_width >= cf_set_width); + if (cf_width < 0 || cf_set_width < 0 || cf_width < cf_set_width) { + amrex::Abort("You must set cf_width >= cf_set_width >= 0"); + } // AmrMesh iterate on grids? bool iterate(true); @@ -1124,6 +1127,10 @@ ERF::ReadParameters () if (verbose > 0) { solverChoice.display(); } + + if (solverChoice.coupling_type == CouplingType::TwoWay && cf_width > 0) { + amrex::Abort("For two-way coupling you must set cf_width = 0"); + } } // Create horizontal average quantities for 5 variables: @@ -1402,12 +1409,12 @@ ERF::Define_ERFFillPatchers (int lev) void ERF::Register_ERFFillPatchers (int lev) { - auto& lev_new = vars_new[lev-1]; - auto& lev_old = vars_old[lev-1]; - FPr_c[lev-1].RegisterCoarseData({&lev_old[Vars::cons], &lev_new[Vars::cons]}, {t_old[lev-1], t_new[lev-1]}); - FPr_u[lev-1].RegisterCoarseData({&lev_old[Vars::xvel], &lev_new[Vars::xvel]}, {t_old[lev-1], t_new[lev-1]}); - FPr_v[lev-1].RegisterCoarseData({&lev_old[Vars::yvel], &lev_new[Vars::yvel]}, {t_old[lev-1], t_new[lev-1]}); - FPr_w[lev-1].RegisterCoarseData({&lev_old[Vars::zvel], &lev_new[Vars::zvel]}, {t_old[lev-1], t_new[lev-1]}); + auto& lev_new = vars_new[lev]; + auto& lev_old = vars_old[lev]; + FPr_c[lev].RegisterCoarseData({&lev_old[Vars::cons], &lev_new[Vars::cons]}, {t_old[lev], t_new[lev]}); + FPr_u[lev].RegisterCoarseData({&lev_old[Vars::xvel], &lev_new[Vars::xvel]}, {t_old[lev], t_new[lev]}); + FPr_v[lev].RegisterCoarseData({&lev_old[Vars::yvel], &lev_new[Vars::yvel]}, {t_old[lev], t_new[lev]}); + FPr_w[lev].RegisterCoarseData({&lev_old[Vars::zvel], &lev_new[Vars::zvel]}, {t_old[lev], t_new[lev]}); } #ifdef ERF_USE_MULTIBLOCK diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index 1f9bbfe4a..b9833011a 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -63,8 +63,10 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba, qmoist[lev].define(ba, dm, 6, ngrow_state); // qv, qc, qi, qr, qs, qg #endif + // ******************************************************************************************** // Build the data structures for calculating diffusive/turbulent terms - update_arrays(lev, ba, dm); + // ******************************************************************************************** + update_diffusive_arrays(lev, ba, dm); // ******************************************************************************************** // Base state holds r_0, pres_0, pi_0 (in that order) @@ -115,6 +117,21 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba, sst_lev[lev].resize(1); sst_lev[lev][0] = nullptr; lmask_lev[lev].resize(1); lmask_lev[lev][0] = nullptr; + // ******************************************************************************************** + // Create the ERFFillPatcher object + // ******************************************************************************************** + if (cf_width > 0 && lev > 0) { + Construct_ERFFillPatchers(lev); + Define_ERFFillPatchers(lev); + } + + // ******************************************************************************************** + // Initialize the boundary conditions + // ******************************************************************************************** + physbcs[lev] = std::make_unique (lev, geom[lev], domain_bcs_type, domain_bcs_type_d, + solverChoice.terrain_type, m_bc_extdir_vals, m_bc_neumann_vals, + z_phys_nd[lev], detJ_cc[lev]); + // ******************************************************************************************** // Initialize the integrator class // ******************************************************************************************** @@ -175,21 +192,73 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, lev_new[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, crse_new[Vars::zvel].nGrowVect()); lev_old[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, crse_new[Vars::zvel].nGrowVect()); + //******************************************************************************************** + // Microphysics + // ******************************************************************************************* +#if defined(ERF_USE_MOISTURE) + int ngrow_state = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_NumDiff) + 1; + qmoist[lev].define(ba, dm, 6, ngrow_state); // qv, qc, qi, qr, qs, qg +#endif + init_stuff(lev, ba, dm); t_new[lev] = time; t_old[lev] = time - 1.e200; + // ******************************************************************************************** // Build the data structures for terrain-related quantities + // ******************************************************************************************** update_terrain_arrays(lev, time); + // ******************************************************************************************** + // Update the base state at this level + // ******************************************************************************************** + base_state[lev].define(ba,dm,3,1); + base_state[lev].setVal(0.); + initHSE(lev); + + // ******************************************************************************************** // Build the data structures for calculating diffusive/turbulent terms - update_arrays(lev, ba, dm); + // ******************************************************************************************** + update_diffusive_arrays(lev, ba, dm); + + // ******************************************************************************************** + // Initialize flux register at new level + // ******************************************************************************************** + if (solverChoice.coupling_type == CouplingType::TwoWay) { + advflux_reg.resize(lev+1); + advflux_reg[lev] = new YAFluxRegister(ba, grids[lev-1], + dm, dmap[lev-1], + geom[lev], geom[lev-1], + ref_ratio[lev-1], lev, NVAR); + } + + // ***************************************************************************************************** + // Initialize the boundary conditions (after initializing the terrain but before calling FillCoarsePatch + // ***************************************************************************************************** + physbcs[lev] = std::make_unique (lev, geom[lev], domain_bcs_type, domain_bcs_type_d, + solverChoice.terrain_type, m_bc_extdir_vals, m_bc_neumann_vals, + z_phys_nd[lev], detJ_cc[lev]); + // ******************************************************************************************** + // Fill data at the new level by interpolation from the coarser level + // ******************************************************************************************** FillCoarsePatch(lev, time, {&lev_new[Vars::cons],&lev_new[Vars::xvel], &lev_new[Vars::yvel],&lev_new[Vars::zvel]}); + // ******************************************************************************************** + // Initialize the integrator class + // ******************************************************************************************** + dt_mri_ratio[lev] = dt_mri_ratio[lev-1]; initialize_integrator(lev, lev_new[Vars::cons], lev_new[Vars::xvel]); + + // ******************************************************************************************** + // If we are making a new level then the FillPatcher for this level hasn't been allocated yet + // ******************************************************************************************** + if (cf_width > 0) { + Construct_ERFFillPatchers(lev); + Define_ERFFillPatchers(lev); + } } // Remake an existing level using provided BoxArray and DistributionMapping and @@ -217,8 +286,29 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp temp_lev_new[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, IntVect(ngrow_vels,ngrow_vels,0)); temp_lev_old[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, IntVect(ngrow_vels,ngrow_vels,0)); + // ******************************************************************************************** + // Remake the flux registers + // ******************************************************************************************** + if (solverChoice.coupling_type == CouplingType::TwoWay) { + delete advflux_reg[lev]; + advflux_reg[lev] = new YAFluxRegister(ba, grids[lev-1], + dm, dmap[lev-1], + geom[lev], geom[lev-1], + ref_ratio[lev-1], lev, NVAR); + } + + //******************************************************************************************** + // Microphysics + // ******************************************************************************************* +#if defined(ERF_USE_MOISTURE) + qmoist[lev].define(ba, dm, 6, ngrow_state); // qv, qc, qi, qr, qs, qg +#endif + init_stuff(lev,ba,dm); + BoxArray ba_old(vars_new[lev][Vars::cons].boxArray()); + DistributionMapping dm_old(vars_new[lev][Vars::cons].DistributionMap()); + // ******************************************************************************************** // This will fill the temporary MultiFabs with data from vars_new // ******************************************************************************************** @@ -245,7 +335,7 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp t_old[lev] = time - 1.e200; // Build the data structures for calculating diffusive/turbulent terms - update_arrays(lev, ba, dm); + update_diffusive_arrays(lev, ba, dm); // Build the data structures for terrain-related quantities update_terrain_arrays(lev, time); @@ -257,12 +347,32 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp base_state[lev].setVal(0.); initHSE(lev); + // ******************************************************************************************** + // Initialize the boundary conditions + // ******************************************************************************************** + physbcs[lev] = std::make_unique (lev, geom[lev], domain_bcs_type, domain_bcs_type_d, + solverChoice.terrain_type, m_bc_extdir_vals, m_bc_neumann_vals, + z_phys_nd[lev], detJ_cc[lev]); + + // ******************************************************************************************** + // Initialize the integrator class + // ******************************************************************************************** + initialize_integrator(lev, vars_new[lev][Vars::cons], vars_new[lev][Vars::xvel]); + + // We need to re-define the FillPatcher if the grids have changed + if (lev > 0 && cf_width > 0) { + bool ba_changed = (ba != ba_old); + bool dm_changed = (dm != dm_old); + if (ba_changed || dm_changed) { + Define_ERFFillPatchers(lev); + } + } } void -ERF::update_arrays (int lev, const BoxArray& ba, const DistributionMapping& dm) +ERF::update_diffusive_arrays (int lev, const BoxArray& ba, const DistributionMapping& dm) { // ******************************************************************************************** // Diffusive terms @@ -374,10 +484,6 @@ ERF::initialize_integrator (int lev, MultiFab& cons_mf, MultiFab& vel_mf) mri_integrator_mem[lev]->setNoSubstepping(solverChoice.no_substepping); mri_integrator_mem[lev]->setIncompressible(solverChoice.incompressible); mri_integrator_mem[lev]->setForceFirstStageSingleSubstep(solverChoice.force_stage1_single_substep); - - physbcs[lev] = std::make_unique (lev, geom[lev], domain_bcs_type, domain_bcs_type_d, - solverChoice.terrain_type, m_bc_extdir_vals, m_bc_neumann_vals, - z_phys_nd[lev], detJ_cc[lev]); } void ERF::init_stuff(int lev, const BoxArray& ba, const DistributionMapping& dm) diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index 95a75e179..b2312d467 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -104,11 +104,14 @@ ERF::WritePlotFile (int which, Vector plot_var_names) if (ncomp_mf == 0) return; - // We fillpatch here because some of the derived quantities require derivatives - // which require ghost cells to be filled + // We Fillpatch here because some of the derived quantities require derivatives + // which require ghost cells to be filled. We do not need to call FillPatcher + // because we don't need to set interior fine points. for (int lev = 0; lev <= finest_level; ++lev) { + bool fillset = false; FillPatch(lev, t_new[lev], {&vars_new[lev][Vars::cons], &vars_new[lev][Vars::xvel], - &vars_new[lev][Vars::yvel], &vars_new[lev][Vars::zvel]}); + &vars_new[lev][Vars::yvel], &vars_new[lev][Vars::zvel]}, + fillset); } Vector mf(finest_level+1); diff --git a/Source/TimeIntegration/ERF_TimeStep.cpp b/Source/TimeIntegration/ERF_TimeStep.cpp index ec7aca82b..69093ce1b 100644 --- a/Source/TimeIntegration/ERF_TimeStep.cpp +++ b/Source/TimeIntegration/ERF_TimeStep.cpp @@ -1,5 +1,4 @@ #include -#include #include using namespace amrex; @@ -25,26 +24,15 @@ ERF::timeStep (int lev, Real time, int iteration) // regrid changes level "lev+1" so we don't regrid on max_level // also make sure we don't regrid fine levels again if // it was taken care of during a coarser regrid - if (lev < max_level && istep[lev] > last_regrid_step[lev]) + if (lev < max_level) { - if (istep[lev] % regrid_int == 0) + if ( (istep[lev] % regrid_int == 0) && (istep[lev] > last_regrid_step[lev]) ) { // regrid could add newly refine levels (if finest_level < max_level) // so we save the previous finest level index int old_finest = finest_level; - regrid(lev, time); - // NOTE: Define & Register index lev backwards (so we add 1 here) - // Redefine & register the ERFFillpatcher objects - if (solverChoice.coupling_type != CouplingType::TwoWay && cf_width > 0) - { - Define_ERFFillPatchers(lev+1); - Register_ERFFillPatchers(lev+1); - if (lev < max_level-1) { - Define_ERFFillPatchers(lev+2); - Register_ERFFillPatchers(lev+2); - } - } + regrid(lev, time); // mark that we have regridded this level already for (int k = lev; k <= finest_level; ++k) { @@ -55,8 +43,8 @@ ERF::timeStep (int lev, Real time, int iteration) for (int k = old_finest+1; k <= finest_level; ++k) { dt[k] = dt[k-1] / MaxRefRatio(k-1); } - } - } + } // if + } // lev } // Update what we call "old" and "new" time @@ -72,9 +60,10 @@ ERF::timeStep (int lev, Real time, int iteration) // Advance a single level for a single time step Advance(lev, time, dt[lev], iteration, nsubsteps[lev]); -#ifdef ERF_USE_PARTICLES - particleData.Redistribute(); -#endif + // We store the old and new data at level "lev" once we have finished the full advance at this level + if (finest_level > 0 && lev < finest_level && cf_set_width > 0) { + Register_ERFFillPatchers(lev); + } ++istep[lev]; @@ -92,12 +81,6 @@ ERF::timeStep (int lev, Real time, int iteration) Real strt_time_for_fine = time + (i-1)*dt[lev+1]; timeStep(lev+1, strt_time_for_fine, i); } - - if (solverChoice.coupling_type == CouplingType::TwoWay) { - int scomp = 0; - int ncomp = NVAR; - AverageDownTo(lev, scomp, ncomp); // average lev+1 down to lev - } } } diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index e9ea41747..3e28792fe 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -228,7 +228,7 @@ #endif // Compute RHS for fine interior ghost - if (level>0 && solverChoice.coupling_type != CouplingType::TwoWay && cf_width>0) { + if (level > 0 && cf_width > 0) { fine_compute_interior_ghost_rhs(new_stage_time, slow_dt, cf_width, cf_set_width, fine_geom, &FPr_c[level-1], &FPr_u[level-1], &FPr_v[level-1], &FPr_w[level-1], @@ -293,7 +293,7 @@ #endif // ************************************************************************* - // Set up flux registers if using mixed or two_way coupling + // Set up flux registers if using two_way coupling // ************************************************************************* YAFluxRegister* fr_as_crse = nullptr; YAFluxRegister* fr_as_fine = nullptr; @@ -401,7 +401,7 @@ #endif // Compute RHS for fine interior ghost - if (level>0 && solverChoice.coupling_type != CouplingType::TwoWay && cf_width>0) { + if (level > 0 && cf_width>0) { fine_compute_interior_ghost_rhs(new_stage_time, slow_dt, cf_width, cf_set_width, &FPr_c[level-1], &FPr_u[level-1], &FPr_v[level-1], &FPr_w[level-1],