From 8248e8fc4c239db7b7e9c05b55e946811b051846 Mon Sep 17 00:00:00 2001 From: Guillaume Giudicelli Date: Wed, 25 Dec 2024 00:11:28 +0100 Subject: [PATCH] Add support for active scalars in SIMPLE --- .../include/executioners/SIMPLESolve.h | 31 +++++++ .../src/executioners/SIMPLESolve.C | 93 ++++++++++++++++++- 2 files changed, 122 insertions(+), 2 deletions(-) diff --git a/modules/navier_stokes/include/executioners/SIMPLESolve.h b/modules/navier_stokes/include/executioners/SIMPLESolve.h index 3931ef4d898c..0787b966e053 100644 --- a/modules/navier_stokes/include/executioners/SIMPLESolve.h +++ b/modules/navier_stokes/include/executioners/SIMPLESolve.h @@ -71,6 +71,37 @@ class SIMPLESolve : public SIMPLESolveBase /// Pointer(s) to the system(s) corresponding to the passive scalar equation(s) std::vector _passive_scalar_systems; + /// Pointer(s) to the system(s) corresponding to the active scalar equation(s) + std::vector _active_scalar_systems; + /// Pointer to the segregated RhieChow interpolation object RhieChowMassFlux * _rc_uo; + + + // ************************ Active Scalar Variables ************************ // + + /// The names of the active scalar systems + const std::vector & _active_scalar_system_names; + + /// Boolean for easy check if a active scalar systems shall be solved or not + const bool _has_active_scalar_systems; + + // The number(s) of the system(s) corresponding to the active scalar equation(s) + std::vector _active_scalar_system_numbers; + + /// The user-defined relaxation parameter(s) for the active scalar equation(s) + const std::vector _active_scalar_equation_relaxation; + + /// Options which hold the petsc settings for the active scalar equation(s) + Moose::PetscSupport::PetscOptions _active_scalar_petsc_options; + + /// Options for the linear solver of the active scalar equation(s) + SIMPLESolverConfiguration _active_scalar_linear_control; + + /// Absolute linear tolerance for the active scalar equation(s). We need to store this, because + /// it needs to be scaled with a representative flux. + const Real _active_scalar_l_abs_tol; + + /// The user-defined absolute tolerance for determining the convergence in active scalars + const std::vector _active_scalar_absolute_tolerance; }; diff --git a/modules/navier_stokes/src/executioners/SIMPLESolve.C b/modules/navier_stokes/src/executioners/SIMPLESolve.C index 5c94e83bcda3..44f94d669829 100644 --- a/modules/navier_stokes/src/executioners/SIMPLESolve.C +++ b/modules/navier_stokes/src/executioners/SIMPLESolve.C @@ -18,6 +18,54 @@ InputParameters SIMPLESolve::validParams() { InputParameters params = SIMPLESolveBase::validParams(); + + /* + * Parameters to control the solution of each scalar advection system + */ + params.addParam>("active_scalar_equation_relaxation", + std::vector(), + "The relaxation which should be used for the active scalar " + "equations. (=1 for no relaxation, " + "diagonal dominance will still be enforced)"); + + params.addParam("active_scalar_petsc_options", + Moose::PetscSupport::getCommonPetscFlags(), + "Singleton PETSc options for the active scalar equation(s)"); + params.addParam( + "active_scalar_petsc_options_iname", + Moose::PetscSupport::getCommonPetscKeys(), + "Names of PETSc name/value pairs for the active scalar equation(s)"); + params.addParam>( + "active_scalar_petsc_options_value", + "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\" for the " + "active scalar equation(s)"); + params.addParam>( + "active_scalar_absolute_tolerance", + std::vector(), + "The absolute tolerance(s) on the normalized residual(s) of the active scalar equation(s)."); + params.addRangeCheckedParam("active_scalar_l_tol", + 1e-5, + "0.0<=active_scalar_l_tol & active_scalar_l_tol<1.0", + "The relative tolerance on the normalized residual in the " + "linear solver of the active scalar equation(s)."); + params.addRangeCheckedParam("active_scalar_l_abs_tol", + 1e-10, + "0.0( + "active_scalar_l_max_its", + 10000, + "The maximum allowed iterations in the linear solver of the turbulence equation."); + + params.addParamNamesToGroup( + "active_scalar_systems active_scalar_equation_relaxation active_scalar_petsc_options " + "active_scalar_petsc_options_iname " + "active_scalar_petsc_options_value active_scalar_petsc_options_value " + "active_scalar_absolute_tolerance " + "active_scalar_l_tol active_scalar_l_abs_tol active_scalar_l_max_its", + "Active Scalars Equations"); + return params; } @@ -28,7 +76,14 @@ SIMPLESolve::SIMPLESolve(Executioner & ex) _energy_sys_number(_has_energy_system ? _problem.linearSysNum(getParam("energy_system")) : libMesh::invalid_uint), - _energy_system(_has_energy_system ? &_problem.getLinearSystem(_energy_sys_number) : nullptr) + _energy_system(_has_energy_system ? &_problem.getLinearSystem(_energy_sys_number) : nullptr), + _active_scalar_system_names(getParam>("active_scalar_systems")), + _has_active_scalar_systems(!_active_scalar_system_names.empty()), + _active_scalar_equation_relaxation( + getParam>("active_scalar_equation_relaxation")), + _active_scalar_l_abs_tol(getParam("active_scalar_l_abs_tol")), + _active_scalar_absolute_tolerance( + getParam>("active_scalar_absolute_tolerance")) { // We fetch the systems and their numbers for the momentum equations. for (auto system_i : index_range(_momentum_system_names)) @@ -45,6 +100,15 @@ SIMPLESolve::SIMPLESolve(Executioner & ex) _passive_scalar_systems.push_back( &_problem.getLinearSystem(_passive_scalar_system_numbers[system_i])); } + // and for the active scalar equations + if (_has_active_scalar_systems) + for (auto system_i : index_range(_active_scalar_system_names)) + { + _active_scalar_system_numbers.push_back( + _problem.linearSysNum(_active_scalar_system_names[system_i])); + _active_scalar_systems.push_back( + &_problem.getLinearSystem(_active_scalar_system_numbers[system_i])); + } } void @@ -284,12 +348,15 @@ SIMPLESolve::solve() unsigned int iteration_counter = 0; // Assign residuals to general residual vector - const unsigned int no_systems = _momentum_systems.size() + 1 + _has_energy_system; + const unsigned int no_systems = _momentum_systems.size() + 1 + _has_energy_system + _active_scalar_absolute_tolerance.size(); std::vector> ns_residuals(no_systems, std::make_pair(0, 1.0)); std::vector ns_abs_tols(_momentum_systems.size(), _momentum_absolute_tolerance); ns_abs_tols.push_back(_pressure_absolute_tolerance); if (_has_energy_system) ns_abs_tols.push_back(_energy_absolute_tolerance); + if (_has_active_scalar_systems) + for (const auto scalar_tol : _active_scalar_absolute_tolerance) + ns_abs_tols.push_back(scalar_tol); bool converged = false; // Loop until converged or hit the maximum allowed iteration number @@ -356,6 +423,21 @@ SIMPLESolve::solve() _energy_linear_control, _energy_l_abs_tol); } + + // If we have active scalar equations, solve them here in case they depend on temperature + if (_has_active_scalar_systems && converged) + { + // We set the preconditioner/controllable parameters through petsc options. Linear + // tolerances will be overridden within the solver. + Moose::PetscSupport::petscSetOptions(_active_scalar_petsc_options, solver_params); + for (const auto i : index_range(_active_scalar_system_names)) + ns_residuals[momentum_residual.size() + 1 + _has_energy_system + i] = solveAdvectedSystem(_active_scalar_system_numbers[i], + *_active_scalar_systems[i], + _active_scalar_equation_relaxation[i], + _active_scalar_linear_control, + _active_scalar_l_abs_tol); + } + _problem.execute(EXEC_NONLINEAR); // Printing residuals @@ -377,6 +459,13 @@ SIMPLESolve::solve() << ns_residuals[momentum_residual.size() + 1].second << COLOR_DEFAULT << " Linear its: " << ns_residuals[momentum_residual.size() + 1].first << std::endl; + if (_has_active_scalar_systems) + _console << " Active scalar equations:\n"; + for (const auto i : index_range(_active_scalar_system_names)) + _console << COLOR_GREEN << " " << _active_scalar_system_names[i] << ": " + << ns_residuals[momentum_residual.size() + 1 + _has_energy_system + i].second << COLOR_DEFAULT + << " Linear its: " << ns_residuals[momentum_residual.size() + 1 + _has_energy_system + i].first << std::endl; + converged = NS::FV::converged(ns_residuals, ns_abs_tols); }