diff --git a/Exec/RegTests/FlameSheet/GNUmakefile b/Exec/RegTests/FlameSheet/GNUmakefile index b19495cf..6ed5bc6b 100644 --- a/Exec/RegTests/FlameSheet/GNUmakefile +++ b/Exec/RegTests/FlameSheet/GNUmakefile @@ -29,12 +29,12 @@ THREAD_SANITIZER = FALSE USE_EFIELD = FALSE # PelePhysics -USE_MANIFOLD = FALSE +USE_MANIFOLD = TRUE ifeq ($(USE_MANIFOLD), TRUE) Eos_Model = Manifold Chemistry_Model = Null Transport_Model = Manifold - Manifold_Dim = 1 + Manifold_Dim = 2 # Manifold_Type = Table else Chemistry_Model = drm19 diff --git a/Exec/RegTests/FlameSheet/input.cvar b/Exec/RegTests/FlameSheet/input.cvar new file mode 100644 index 00000000..549bf258 --- /dev/null +++ b/Exec/RegTests/FlameSheet/input.cvar @@ -0,0 +1,101 @@ +#----------------------DOMAIN DEFINITION------------------------ +geometry.is_periodic = 0 0 # For each dir, 0: non-perio, 1: periodic +geometry.coord_sys = 0 # 0 => cart, 1 => RZ +geometry.prob_lo = 0.0 0.0 0.0 # x_lo y_lo (z_lo) +geometry.prob_hi = 0.003 0.012 0.016 # x_hi y_hi (z_hi) + +# >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< +# Interior, Inflow, Outflow, Symmetry, +# SlipWallAdiab, NoSlipWallAdiab, SlipWallIsotherm, NoSlipWallIsotherm +peleLM.lo_bc = Interior Inflow +peleLM.hi_bc = Interior Outflow +peleLM.lo_bc = SlipWallAdiab Inflow +peleLM.hi_bc = SlipWallAdiab Outflow + + +#-------------------------AMR CONTROL---------------------------- +amr.n_cell = 64 256 32 # Level 0 number of cells in each direction +amr.v = 1 # AMR verbose +amr.max_level = 0 # maximum level number allowed +amr.ref_ratio = 2 2 2 2 # refinement ratio +amr.regrid_int = 5 # how often to regrid +amr.n_error_buf = 1 1 2 2 # number of buffer cells in error est +amr.grid_eff = 0.7 # what constitutes an efficient grid +amr.blocking_factor = 16 # block factor in grid generation (min box size) +amr.max_grid_size = 64 # max box size + + +#--------------------------- Problem ------------------------------- +prob.P_mean = 101325.0 +prob.standoff = -.03 +prob.pertmag = 0.0000 +pmf.datafile = "prem_drm19_phi1_p1_t298_mani_filter2.344e-03.dat" +#pmf.datafile = "prem_drm19_phi1_p1_t298_mani.dat" + +#-------------------------PeleLM CONTROL---------------------------- +peleLM.v = 3 +peleLM.incompressible = 0 +peleLM.sdc_iterMax = 1 +peleLM.floor_species = 0 + +#peleLM.num_init_iter = 0 +#peleLM.num_divu_iter = 0 + +peleLM.do_temporals = 1 +peleLM.temporal_int = 2 +peleLM.do_mass_balance = 1 + +#amr.restart = chk00005 +#amr.check_int = 20 +amr.plot_int = 20 +amr.max_step = 1000 +amr.dt_shrink = 0.01 +amr.stop_time = 0.01 +#amr.stop_time = 1.00 +amr.cfl = 0.5 +amr.derive_plot_vars = avg_pressure mag_vort mass_fractions maniout progress_variable + +peleLM.progressVariable.format = Cantera +peleLM.progressVariable.weights = "PROG:1" +peleLM.progressVariable.coldState = "PROG:0" +peleLM.progressVariable.hotState = "PROG:1" + +# ------------------- INPUTS DERIVED DIAGS ------------------ +peleLM.fuel_name = CH4 + +# --------------- INPUTS TO CHEMISTRY REACTOR --------------- +peleLM.chem_integrator = "ReactorRK64" + +mac_proj.verbose = 0 +nodal_proj.verbose = 0 + +#--------------------REFINEMENT CONTROL------------------------ +amr.refinement_indicators = temp +amr.temp.max_level = 3 +amr.temp.adjacent_difference_greater = 10 +amr.temp.field_name = temp + +#amrex.fpe_trap_invalid = 1 +#amrex.fpe_trap_zero = 1 +#amrex.fpe_trap_overflow = 1 + +#----------------------- Manifold Stuff ---------------------------- +manifold.model = Table +manifold.v = 1 +manifold.table.v = 2 +manifold.table.filename = prem_drm19_phi1_p1_t298.ctb +manifold.nominal_pressure_cgs = 1013250.0 +manifold.compute_temperature = 1 +manifold.density_lookup_type = log +peleLM.use_wbar = 0 +peleLM.chi_correction_type = DivuFirstIter +#peleLM.chi_correction_type = NoDivu +peleLM.print_chi_convergence = 1 + +peleLM.deltaT_iterMax = -1 +peleLM.les_model = "Smagorinsky" +peleLM.les_cs_smag = 0.18 +peleLM.les_v = 1 +peleLM.plot_les = 1 +amr.plot_extSource = 1 +peleLM.les_c_chi = 2.0 \ No newline at end of file diff --git a/Source/PeleLMeX.H b/Source/PeleLMeX.H index 975d352e..49999a93 100644 --- a/Source/PeleLMeX.H +++ b/Source/PeleLMeX.H @@ -827,6 +827,8 @@ public: void addSpark(const PeleLM::TimeStamp& a_time); + void addScalarDissipation(const PeleLM::TimeStamp& a_time); + //----------------------------------------------------------------------------- #ifdef PELE_USE_SPRAY @@ -1802,6 +1804,7 @@ public: int m_plotChemDiag = 0; int m_plotHeatRelease = 0; int m_plotStateSpec = 0; + bool m_plotExtSource = false; int m_plot_int = 0; bool m_plot_overwrite = false; int m_plot_zeroEBcovered = 1; @@ -1911,7 +1914,7 @@ public: amrex::Real m_les_cs_smag = 0.18; amrex::Real m_les_cm_wale = 0.60; amrex::Real m_les_cs_sigma = 1.35; - amrex::Vector m_turb_visc_time; + amrex::Real m_les_c_chi = 20.0; // switch on/off different processes // TODO: actually implement that diff --git a/Source/PeleLMeX_Advance.cpp b/Source/PeleLMeX_Advance.cpp index b011dbbc..eb77f394 100644 --- a/Source/PeleLMeX_Advance.cpp +++ b/Source/PeleLMeX_Advance.cpp @@ -105,6 +105,9 @@ PeleLM::Advance(int is_initIter) poissonSolveEF(AmrOldTime); #endif } + if (m_do_les) { + calcTurbViscosity(AmrOldTime); + } //---------------------------------------------------------------- BL_PROFILE_VAR_STOP(PLM_SETUP); diff --git a/Source/PeleLMeX_Evaluate.cpp b/Source/PeleLMeX_Evaluate.cpp index 11b07f28..53090904 100644 --- a/Source/PeleLMeX_Evaluate.cpp +++ b/Source/PeleLMeX_Evaluate.cpp @@ -171,6 +171,10 @@ PeleLM::MLevaluate( finest_level, grids, dmap, m_factory, m_nGrowAdv, m_use_wbar, m_use_soret); calcDiffusivity(AmrNewTime); + // If doing LES, need to be able to add in turbulent component + if (m_do_les) { + calcTurbViscosity(AmrNewTime); + } computeDifferentialDiffusionTerms(AmrNewTime, diffData); for (int lev = 0; lev <= finest_level; ++lev) { MultiFab::Copy( @@ -215,6 +219,7 @@ PeleLM::MLevaluate( } else if (a_var == "transportCC") { // Cell-centered transport coefficients functions go through the level // data container. Simply copy once the later has been filled. + // No LES component here - this is just moledular transport coeffs calcViscosity(AmrNewTime); calcDiffusivity(AmrNewTime); for (int lev = 0; lev <= finest_level; ++lev) { @@ -290,6 +295,9 @@ PeleLM::evaluateChemExtForces( // compute t^{n} data calcViscosity(AmrOldTime); calcDiffusivity(AmrOldTime); + if (m_do_les) { + calcTurbViscosity(AmrOldTime); + } floorSpecies(AmrOldTime); setThermoPress(AmrOldTime); @@ -392,6 +400,9 @@ PeleLM::evaluateAdvectionTerms( // compute t^{n} data calcViscosity(AmrOldTime); calcDiffusivity(AmrOldTime); + if (m_do_les) { + calcTurbViscosity(AmrOldTime); + } floorSpecies(AmrOldTime); setThermoPress(AmrOldTime); diff --git a/Source/PeleLMeX_Forces.cpp b/Source/PeleLMeX_Forces.cpp index 502fda58..bbfd3e7e 100644 --- a/Source/PeleLMeX_Forces.cpp +++ b/Source/PeleLMeX_Forces.cpp @@ -237,6 +237,82 @@ PeleLM::addSpark(const TimeStamp& a_timestamp) } } +// Manifold model - dissipation rate sources for variances +void +PeleLM::addScalarDissipation(const TimeStamp& a_timestamp) +{ + BL_PROFILE("PeleLM::addScalarDissipationRate"); + // no scalar dissipation sources if not using a manifold model +#ifndef USE_MANIFOLD_EOS + amrex::ignore_unused(a_timestamp); +#else + + // Dissipation source term for subfilter variances + for (int lev = 0; lev <= finest_level; lev++) { + + auto* ldata_p = getLevelDataPtr(lev, a_timestamp); + auto const& leosparm = eos_parms.host_parm(); + + // For now - we require the turbulent viscosity to be pre-computed + // it always is stored at AmrOldTime, so we just use that + // We need cell-centered mu_t but have it at faces + // The simple interpolation below probably isn't valid for EB +#ifdef AMREX_USE_EB + amrex::Abort( + "PeleLM::addScalarDissipation(): this is not supported with EB"); +#endif + + for (int n = 0; n < MANIFOLD_DIM; ++n) { + if (leosparm.is_variance_of[n] >= 0) { + if (!m_do_les) { + amrex::Abort("PeleLM::addScalarDissipation(): cannot add a " + "scalarDissipation without an active LES model"); + } + + constexpr amrex::Real fact = 0.5 / AMREX_SPACEDIM; + const amrex::Real C_chi = m_les_c_chi; + + AMREX_D_TERM(auto const& mut_arr_x = + m_leveldata_old[lev]->visc_turb_fc[0].const_arrays(); + , auto const& mut_arr_y = + m_leveldata_old[lev]->visc_turb_fc[1].const_arrays(); + , auto const& mut_arr_z = + m_leveldata_old[lev]->visc_turb_fc[2].const_arrays();) + auto extma = m_extSource[lev]->arrays(); + auto statema = ldata_p->state.const_arrays(); + + // l_scale will also need modification for EB + const amrex::Real vol = AMREX_D_TERM( + geom[lev].CellSize(0), *geom[lev].CellSize(1), + *geom[lev].CellSize(2)); + const amrex::Real l_scale = + (AMREX_SPACEDIM == 2) ? std::sqrt(vol) : std::cbrt(vol); + const amrex::Real inv_l_scale2 = 1.0 / (l_scale * l_scale); + + amrex::ParallelFor( + *m_extSource[lev], + [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { + amrex::Real mu_t = + fact * + (AMREX_D_TERM( + mut_arr_x[box_no](i, j, k) + mut_arr_x[box_no](i + 1, j, k), + +mut_arr_y[box_no](i, j, k) + mut_arr_y[box_no](i, j + 1, k), + +mut_arr_z[box_no](i, j, k) + mut_arr_z[box_no](i, j, k + 1))); + + // Linear Relaxation model + // rho chi_sgs = C_chi * mu_t / Delta^2 * Variance + extma[box_no](i, j, k, FIRSTSPEC + n) -= + C_chi * mu_t * inv_l_scale2 * + statema[box_no](i, j, k, FIRSTSPEC + n); + }); + } + } + Gpu::streamSynchronize(); + } + +#endif +} + // Calculate additional external sources (soot, radiation, user defined, etc.) void PeleLM::getExternalSources( @@ -268,6 +344,8 @@ PeleLM::getExternalSources( } #endif + addScalarDissipation(a_timestamp_old); + // User defined external sources if (m_user_defined_ext_sources) { for (int lev = 0; lev <= finest_level; lev++) { @@ -279,4 +357,4 @@ PeleLM::getExternalSources( ldata_p_new->state, ext_src, geom[lev].data(), *prob_parm_d); } } -} \ No newline at end of file +} diff --git a/Source/PeleLMeX_Plot.cpp b/Source/PeleLMeX_Plot.cpp index 51a3b02c..a2e1cb57 100644 --- a/Source/PeleLMeX_Plot.cpp +++ b/Source/PeleLMeX_Plot.cpp @@ -165,6 +165,10 @@ PeleLM::WritePlotFile() ncomp += 1; } + if (m_plotExtSource) { + ncomp += NVAR; + } + //---------------------------------------------------------------- // Plot MultiFabs Vector mf_plt(finest_level + 1); @@ -282,6 +286,12 @@ PeleLM::WritePlotFile() plt_VarsName.push_back("viscturb"); } + if (m_plotExtSource) { + for (int ivar = 0; ivar < NVAR; ++ivar) { + plt_VarsName.push_back("extsource_" + stateVariableName(ivar)); + } + } + #if NUM_ODE > 0 for (int n = 0; n < NUM_ODE; n++) { plt_VarsName.push_back(m_ode_names[n]); @@ -437,6 +447,12 @@ PeleLM::WritePlotFile() }); Gpu::streamSynchronize(); } + cnt += 1; + + if (m_plotExtSource) { + MultiFab::Copy(mf_plt[lev], *m_extSource[lev], 0, cnt, NVAR, 0); + // cnt += NVAR; + } #ifdef AMREX_USE_EB if (m_plot_zeroEBcovered != 0) { diff --git a/Source/PeleLMeX_Setup.cpp b/Source/PeleLMeX_Setup.cpp index d77a4dc5..d812519e 100644 --- a/Source/PeleLMeX_Setup.cpp +++ b/Source/PeleLMeX_Setup.cpp @@ -413,9 +413,7 @@ PeleLM::readParameters() m_les_verbose = m_verbose; pp.query("plot_les", m_plot_les); pp.query("les_v", m_les_verbose); - for (int lev = 0; lev <= max_level; ++lev) { - m_turb_visc_time.push_back(-1.0E200); - } + pp.query("les_c_chi", m_les_c_chi); #ifdef PELE_USE_EFIELD amrex::Abort("LES implementation is not yet compatible with efield/ions"); #endif @@ -837,6 +835,7 @@ PeleLM::readIOParameters() } pp.query("plot_zeroEBcovered", m_plot_zeroEBcovered); pp.query("plot_speciesState", m_plotStateSpec); + pp.query("plot_extSource", m_plotExtSource); m_initial_grid_file = ""; m_regrid_file = ""; pp.query("initial_grid_file", m_initial_grid_file); diff --git a/Source/PeleLMeX_TransportProp.cpp b/Source/PeleLMeX_TransportProp.cpp index 23a77f96..cde9af38 100644 --- a/Source/PeleLMeX_TransportProp.cpp +++ b/Source/PeleLMeX_TransportProp.cpp @@ -414,24 +414,9 @@ PeleLM::getDiffusivity( // supported if ((addTurbContrib != 0) and m_do_les) { - // If initializing the simulation, always recompute turbulent viscosity - // otherwise, only recompute once per level per timestep (at old time) - // calcTurbViscosity computes for all levels, so only call from the base - // level - TimeStamp tstamp; - if (getTime(lev, AmrNewTime) == 0.0) { - tstamp = AmrNewTime; - if (lev == 0) { - calcTurbViscosity(tstamp); - } - } else if (lev == 0 and getTime(lev, AmrOldTime) > m_turb_visc_time[lev]) { - tstamp = AmrOldTime; - calcTurbViscosity(tstamp); - m_turb_visc_time[lev] = getTime(lev, AmrOldTime); - } else { - tstamp = AmrOldTime; - } - auto* ldata_p = getLevelDataPtr(lev, tstamp); + // Turbulent viscosity - always in old data + // (not updated because velocity doesn't get updated until end of advance) + auto* ldata_p = getLevelDataPtr(lev, AmrOldTime); // Identify and add the correct turbulent contribution for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { diff --git a/Submodules/PelePhysics b/Submodules/PelePhysics index 760e2feb..f523711a 160000 --- a/Submodules/PelePhysics +++ b/Submodules/PelePhysics @@ -1 +1 @@ -Subproject commit 760e2febb1573bc9f1a3dd9b70ac5a36ebc8ce5c +Subproject commit f523711a3bf579b9a8c82c53009ac8ddd428a7ef