diff --git a/Source/PeleLMeX.H b/Source/PeleLMeX.H index 975d352e..c6ff4fa3 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; + int m_plotExtSource = 0; 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_Forces.cpp b/Source/PeleLMeX_Forces.cpp index d5170821..82b5a502 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 + return; +#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..243e1c63 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/Submodules/PelePhysics b/Submodules/PelePhysics index 760e2feb..c3c955f6 160000 --- a/Submodules/PelePhysics +++ b/Submodules/PelePhysics @@ -1 +1 @@ -Subproject commit 760e2febb1573bc9f1a3dd9b70ac5a36ebc8ce5c +Subproject commit c3c955f65d730689e1ee7065312c02b6541f8a1a