https://mooseframework.inl.gov
LinearAssemblySegregatedSolve.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://mooseframework.inl.gov
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
11 #include "FEProblem.h"
12 #include "SegregatedSolverUtils.h"
13 #include "LinearSystem.h"
14 
15 using namespace libMesh;
16 
19 {
21 
22  params.addParam<std::vector<SolverSystemName>>(
23  "active_scalar_systems", {}, "The solver system for each active scalar advection equation.");
24 
25  /*
26  * Parameters to control the solution of each scalar advection system
27  */
28  params.addParam<std::vector<Real>>("active_scalar_equation_relaxation",
29  std::vector<Real>(),
30  "The relaxation which should be used for the active scalar "
31  "equations. (=1 for no relaxation, "
32  "diagonal dominance will still be enforced)");
33 
34  params.addParam<MultiMooseEnum>("active_scalar_petsc_options",
36  "Singleton PETSc options for the active scalar equation(s)");
37  params.addParam<MultiMooseEnum>(
38  "active_scalar_petsc_options_iname",
40  "Names of PETSc name/value pairs for the active scalar equation(s)");
41  params.addParam<std::vector<std::string>>(
42  "active_scalar_petsc_options_value",
43  "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\" for the "
44  "active scalar equation(s)");
45  params.addParam<std::vector<Real>>(
46  "active_scalar_absolute_tolerance",
47  std::vector<Real>(),
48  "The absolute tolerance(s) on the normalized residual(s) of the active scalar equation(s).");
49  params.addRangeCheckedParam<Real>("active_scalar_l_tol",
50  1e-5,
51  "0.0<=active_scalar_l_tol & active_scalar_l_tol<1.0",
52  "The relative tolerance on the normalized residual in the "
53  "linear solver of the active scalar equation(s).");
54  params.addRangeCheckedParam<Real>("active_scalar_l_abs_tol",
55  1e-10,
56  "0.0<active_scalar_l_abs_tol",
57  "The absolute tolerance on the normalized residual in the "
58  "linear solver of the active scalar equation(s).");
59  params.addParam<unsigned int>(
60  "active_scalar_l_max_its",
61  10000,
62  "The maximum allowed iterations in the linear solver of the turbulence equation.");
63 
64  params.addParamNamesToGroup(
65  "active_scalar_systems active_scalar_equation_relaxation active_scalar_petsc_options "
66  "active_scalar_petsc_options_iname "
67  "active_scalar_petsc_options_value active_scalar_petsc_options_value "
68  "active_scalar_absolute_tolerance "
69  "active_scalar_l_tol active_scalar_l_abs_tol active_scalar_l_max_its",
70  "Active Scalars Equations");
71 
72  return params;
73 }
74 
76  : SIMPLESolveBase(ex),
77  _pressure_sys_number(_problem.linearSysNum(getParam<SolverSystemName>("pressure_system"))),
78  _pressure_system(_problem.getLinearSystem(_pressure_sys_number)),
79  _energy_sys_number(_has_energy_system
80  ? _problem.linearSysNum(getParam<SolverSystemName>("energy_system"))
82  _energy_system(_has_energy_system ? &_problem.getLinearSystem(_energy_sys_number) : nullptr),
83  _solid_energy_sys_number(
84  _has_solid_energy_system
85  ? _problem.linearSysNum(getParam<SolverSystemName>("solid_energy_system"))
87  _solid_energy_system(
88  _has_solid_energy_system ? &_problem.getLinearSystem(_solid_energy_sys_number) : nullptr),
89  _active_scalar_system_names(getParam<std::vector<SolverSystemName>>("active_scalar_systems")),
90  _has_active_scalar_systems(!_active_scalar_system_names.empty()),
91  _active_scalar_equation_relaxation(
92  getParam<std::vector<Real>>("active_scalar_equation_relaxation")),
93  _active_scalar_l_abs_tol(getParam<Real>("active_scalar_l_abs_tol")),
94  _active_scalar_absolute_tolerance(
95  getParam<std::vector<Real>>("active_scalar_absolute_tolerance"))
96 {
97  // We fetch the systems and their numbers for the momentum equations.
98  for (auto system_i : index_range(_momentum_system_names))
99  {
102  _systems_to_solve.push_back(_momentum_systems.back());
103  }
104 
106 
107  if (_has_energy_system)
109 
112  // and for the turbulence surrogate equations
114  for (auto system_i : index_range(_turbulence_system_names))
115  {
116  _turbulence_system_numbers.push_back(
118  _turbulence_systems.push_back(
120  }
121 
122  // and for the passive scalar equations
124  for (auto system_i : index_range(_passive_scalar_system_names))
125  {
128  _passive_scalar_systems.push_back(
130  _systems_to_solve.push_back(_passive_scalar_systems.back());
131  }
132 
133  // and for the active scalar equations
135  for (auto system_i : index_range(_active_scalar_system_names))
136  {
139  _active_scalar_systems.push_back(
141  _systems_to_solve.push_back(_active_scalar_systems.back());
142 
143  const auto & active_scalar_petsc_options =
144  getParam<MultiMooseEnum>("active_scalar_petsc_options");
145  const auto & active_scalar_petsc_pair_options = getParam<MooseEnumItem, std::string>(
146  "active_scalar_petsc_options_iname", "active_scalar_petsc_options_value");
148  active_scalar_petsc_options, "-", *this, _active_scalar_petsc_options);
149  Moose::PetscSupport::addPetscPairsToPetscOptions(active_scalar_petsc_pair_options,
150  _problem.mesh().dimension(),
151  "-",
152  *this,
154 
156  getParam<Real>("active_scalar_l_tol");
158  getParam<Real>("active_scalar_l_abs_tol");
160  getParam<unsigned int>("active_scalar_l_max_its");
161  }
162 
164  paramError("active_scalar_equation_relaxation",
165  "Should be the same size as the number of systems");
166 
167  // We disable the prefix here for the time being, the segregated solvers use a different approach
168  // for setting the petsc parameters
169  for (auto & system : _systems_to_solve)
170  system->system().prefix_with_name(false);
171 }
172 
173 void
175 {
176  _rc_uo =
177  const_cast<RhieChowMassFlux *>(&getUserObject<RhieChowMassFlux>("rhie_chow_user_object"));
180 
181  // Initialize the face velocities in the RC object
182  if (!_app.isRecovering())
185 }
186 
187 std::vector<std::pair<unsigned int, Real>>
189 {
190  // Temporary storage for the (flux-normalized) residuals from
191  // different momentum components
192  std::vector<std::pair<unsigned int, Real>> its_normalized_residuals;
193 
194  LinearImplicitSystem & momentum_system_0 =
195  libMesh::cast_ref<LinearImplicitSystem &>(_momentum_systems[0]->system());
196 
197  PetscLinearSolver<Real> & momentum_solver =
198  libMesh::cast_ref<PetscLinearSolver<Real> &>(*momentum_system_0.get_linear_solver());
199 
200  // Solve the momentum equations.
201  // TO DO: These equations are VERY similar. If we can store the differences (things coming from
202  // BCs for example) separately, it is enough to construct one matrix.
203  for (const auto system_i : index_range(_momentum_systems))
204  {
206 
207  // We will need the right hand side and the solution of the next component
208  LinearImplicitSystem & momentum_system =
209  libMesh::cast_ref<LinearImplicitSystem &>(_momentum_systems[system_i]->system());
210 
211  NumericVector<Number> & solution = *(momentum_system.solution);
212  NumericVector<Number> & rhs = *(momentum_system.rhs);
213  SparseMatrix<Number> & mmat = *(momentum_system.matrix);
214 
215  auto diff_diagonal = solution.zero_clone();
216 
217  // We assemble the matrix and the right hand side
218  _problem.computeLinearSystemSys(momentum_system, mmat, rhs, /*compute_grads*/ true);
219 
220  // Still need to relax the right hand side with the same vector
221  NS::FV::relaxMatrix(mmat, _momentum_equation_relaxation, *diff_diagonal);
222  NS::FV::relaxRightHandSide(rhs, solution, *diff_diagonal);
223 
224  // The normalization factor depends on the right hand side so we need to recompute it for this
225  // component
226  Real norm_factor = NS::FV::computeNormalizationFactor(solution, mmat, rhs);
227 
228  // Very important, for deciding the convergence, we need the unpreconditioned
229  // norms in the linear solve
230  LibmeshPetscCall(KSPSetNormType(momentum_solver.ksp(), KSP_NORM_UNPRECONDITIONED));
231  // Solve this component. We don't update the ghosted solution yet, that will come at the end
232  // of the corrector step. Also setting the linear tolerances and maximum iteration counts.
234  momentum_solver.set_solver_configuration(_momentum_linear_control);
235 
236  // We solve the equation
237  auto its_resid_pair = momentum_solver.solve(mmat, mmat, solution, rhs);
238  momentum_system.update();
239 
240  // We will reuse the preconditioner for every momentum system
241  if (system_i == 0)
242  momentum_solver.reuse_preconditioner(true);
243 
244  // Save the normalized residual
245  its_normalized_residuals.push_back(
246  std::make_pair(its_resid_pair.first, momentum_solver.get_initial_residual() / norm_factor));
247 
248  if (_print_fields)
249  {
250  _console << " solution after solve " << std::endl;
251  solution.print();
252  _console << " matrix when we solve " << std::endl;
253  mmat.print();
254  _console << " rhs when we solve " << std::endl;
255  rhs.print();
256  _console << " velocity solution component " << system_i << std::endl;
257  solution.print();
258  _console << "Norm factor " << norm_factor << std::endl;
259  _console << Moose::stringify(momentum_solver.get_initial_residual()) << std::endl;
260  }
261 
262  // Printing residuals
263  _console << " Momentum equation:"
264  << (_momentum_systems.size() > 1
265  ? std::string(" Component ") + std::to_string(system_i + 1) + std::string(" ")
266  : std::string(" "))
267  << COLOR_GREEN << its_normalized_residuals[system_i].second << COLOR_DEFAULT
268  << " Linear its: " << its_normalized_residuals[system_i].first << std::endl;
269  }
270 
271  for (const auto system_i : index_range(_momentum_systems))
272  {
273  LinearImplicitSystem & momentum_system =
274  libMesh::cast_ref<LinearImplicitSystem &>(_momentum_systems[system_i]->system());
275  _momentum_systems[system_i]->setSolution(*(momentum_system.current_local_solution));
276  _momentum_systems[system_i]->copyPreviousNonlinearSolutions();
277  }
278 
279  // We reset this to ensure the preconditioner is recomputed new time we go to the momentum
280  // predictor
281  momentum_solver.reuse_preconditioner(false);
282 
283  return its_normalized_residuals;
284 }
285 
286 std::pair<unsigned int, Real>
288 {
290 
291  // We will need some members from the linear system
292  LinearImplicitSystem & pressure_system =
293  libMesh::cast_ref<LinearImplicitSystem &>(_pressure_system.system());
294 
295  // We will need the solution, the right hand side and the matrix
296  NumericVector<Number> & current_local_solution = *(pressure_system.current_local_solution);
297  NumericVector<Number> & solution = *(pressure_system.solution);
298  SparseMatrix<Number> & mmat = *(pressure_system.matrix);
299  NumericVector<Number> & rhs = *(pressure_system.rhs);
300 
301  // Fetch the linear solver from the system
302  PetscLinearSolver<Real> & pressure_solver =
303  libMesh::cast_ref<PetscLinearSolver<Real> &>(*pressure_system.get_linear_solver());
304 
305  _problem.computeLinearSystemSys(pressure_system, mmat, rhs, false);
306 
307  if (_print_fields)
308  {
309  _console << "Pressure matrix" << std::endl;
310  mmat.print();
311  }
312 
313  // We compute the normalization factors based on the fluxes
314  Real norm_factor = NS::FV::computeNormalizationFactor(solution, mmat, rhs);
315 
316  // We need the non-preconditioned norm to be consistent with the norm factor
317  LibmeshPetscCall(KSPSetNormType(pressure_solver.ksp(), KSP_NORM_UNPRECONDITIONED));
318 
319  // Setting the linear tolerances and maximum iteration counts
322 
323  if (_pin_pressure)
325  pressure_system.update();
326 
327  auto its_res_pair = pressure_solver.solve(mmat, mmat, solution, rhs);
328  pressure_system.update();
329 
330  if (_print_fields)
331  {
332  _console << " rhs when we solve pressure " << std::endl;
333  rhs.print();
334  _console << " Pressure " << std::endl;
335  solution.print();
336  _console << "Norm factor " << norm_factor << std::endl;
337  }
338 
339  _pressure_system.setSolution(current_local_solution);
340 
341  const auto residuals =
342  std::make_pair(its_res_pair.first, pressure_solver.get_initial_residual() / norm_factor);
343 
344  _console << " Pressure equation: " << COLOR_GREEN << residuals.second << COLOR_DEFAULT
345  << " Linear its: " << residuals.first << std::endl;
346 
347  return residuals;
348 }
349 
350 std::pair<unsigned int, Real>
352 {
354 
355  // We will need some members from the linear system
356  LinearImplicitSystem & system =
357  libMesh::cast_ref<LinearImplicitSystem &>(_solid_energy_system->system());
358 
359  // We will need the solution, the right hand side and the matrix
360  NumericVector<Number> & current_local_solution = *(system.current_local_solution);
361  NumericVector<Number> & solution = *(system.solution);
362  SparseMatrix<Number> & mmat = *(system.matrix);
363  NumericVector<Number> & rhs = *(system.rhs);
364 
365  // Fetch the linear solver from the system
366  PetscLinearSolver<Real> & solver =
367  libMesh::cast_ref<PetscLinearSolver<Real> &>(*system.get_linear_solver());
368 
369  _problem.computeLinearSystemSys(system, mmat, rhs, false);
370 
371  if (_print_fields)
372  {
373  _console << "Solid energy matrix" << std::endl;
374  mmat.print();
375  }
376 
377  // We compute the normalization factors based on the fluxes
378  Real norm_factor = NS::FV::computeNormalizationFactor(solution, mmat, rhs);
379 
380  // We need the non-preconditioned norm to be consistent with the norm factor
381  LibmeshPetscCall(KSPSetNormType(solver.ksp(), KSP_NORM_UNPRECONDITIONED));
382 
383  // Setting the linear tolerances and maximum iteration counts
386 
387  auto its_res_pair = solver.solve(mmat, mmat, solution, rhs);
388  system.update();
389 
390  if (_print_fields)
391  {
392  _console << " rhs when we solve solid energy " << std::endl;
393  rhs.print();
394  _console << " Solid energy " << std::endl;
395  solution.print();
396  _console << "Norm factor " << norm_factor << std::endl;
397  }
398 
399  _solid_energy_system->setSolution(current_local_solution);
400 
401  const auto residuals =
402  std::make_pair(its_res_pair.first, solver.get_initial_residual() / norm_factor);
403 
404  _console << " Solid energy equation: " << COLOR_GREEN << residuals.second << COLOR_DEFAULT
405  << " Linear its: " << residuals.first << std::endl;
406 
407  return residuals;
408 }
409 
410 std::pair<unsigned int, Real>
411 LinearAssemblySegregatedSolve::correctVelocity(const bool subtract_updated_pressure,
412  const bool recompute_face_mass_flux,
413  const SolverParams & solver_params)
414 {
415  // Compute the coupling fields between the momentum and pressure equations.
416  // The first argument makes sure the pressure gradient is staged at the first
417  // iteration
418  _rc_uo->computeHbyA(subtract_updated_pressure, _print_fields);
419 
420  // We set the preconditioner/controllable parameters for the pressure equations through
421  // petsc options. Linear tolerances will be overridden within the solver.
423 
424  // Solve the pressure corrector
425  const auto residuals = solvePressureCorrector();
426 
427  // Compute the face velocity which is used in the advection terms. In certain
428  // segregated solver algorithms (like PISO) this is only done on the last iteration.
429  if (recompute_face_mass_flux)
431 
432  auto & pressure_current_solution = *(_pressure_system.system().current_local_solution.get());
433  auto & pressure_old_solution = *(_pressure_system.solutionPreviousNewton());
434 
435  // Relax the pressure update for the next momentum predictor
437  pressure_current_solution, pressure_old_solution, _pressure_variable_relaxation);
438 
439  // Overwrite old solution
440  pressure_old_solution = pressure_current_solution;
441  _pressure_system.setSolution(pressure_current_solution);
442 
443  // We recompute the updated pressure gradient
445 
446  // Reconstruct the cell velocity as well to accelerate convergence
448 
449  return residuals;
450 }
451 
452 std::pair<unsigned int, Real>
454  LinearSystem & system,
455  const Real relaxation_factor,
456  SolverConfiguration & solver_config,
457  const Real absolute_tol,
458  const Real field_relaxation,
459  const Real min_value_limiter)
460 {
461  _problem.setCurrentLinearSystem(system_num);
462 
463  // We will need some members from the implicit linear system
464  LinearImplicitSystem & li_system = libMesh::cast_ref<LinearImplicitSystem &>(system.system());
465 
466  // We will need the solution, the right hand side and the matrix
467  NumericVector<Number> & current_local_solution = *(li_system.current_local_solution);
468  NumericVector<Number> & solution = *(li_system.solution);
469  SparseMatrix<Number> & mmat = *(li_system.matrix);
470  NumericVector<Number> & rhs = *(li_system.rhs);
471 
472  // We need a vector that stores the (diagonal_relaxed-original_diagonal) vector
473  auto diff_diagonal = solution.zero_clone();
474 
475  // Fetch the linear solver from the system
476  PetscLinearSolver<Real> & linear_solver =
477  libMesh::cast_ref<PetscLinearSolver<Real> &>(*li_system.get_linear_solver());
478 
479  _problem.computeLinearSystemSys(li_system, mmat, rhs, true);
480 
481  // Go and relax the system matrix and the right hand side
482  NS::FV::relaxMatrix(mmat, relaxation_factor, *diff_diagonal);
483  NS::FV::relaxRightHandSide(rhs, solution, *diff_diagonal);
484 
485  if (_print_fields)
486  {
487  _console << system.name() << " system matrix" << std::endl;
488  mmat.print();
489  }
490 
491  // We compute the normalization factors based on the fluxes
492  Real norm_factor = NS::FV::computeNormalizationFactor(solution, mmat, rhs);
493 
494  // We need the non-preconditioned norm to be consistent with the norm factor
495  LibmeshPetscCall(KSPSetNormType(linear_solver.ksp(), KSP_NORM_UNPRECONDITIONED));
496 
497  // Setting the linear tolerances and maximum iteration counts
498  solver_config.real_valued_data["abs_tol"] = absolute_tol * norm_factor;
499  linear_solver.set_solver_configuration(solver_config);
500 
501  // Solve the system and update current local solution
502  auto its_res_pair = linear_solver.solve(mmat, mmat, solution, rhs);
503  li_system.update();
504 
505  if (_print_fields)
506  {
507  _console << " rhs when we solve " << system.name() << std::endl;
508  rhs.print();
509  _console << system.name() << " solution " << std::endl;
510  solution.print();
511  _console << " Norm factor " << norm_factor << std::endl;
512  }
513 
514  // Limiting scalar solution
515  if (min_value_limiter != std::numeric_limits<Real>::min())
516  NS::FV::limitSolutionUpdate(current_local_solution, min_value_limiter);
517 
518  // Relax the field update for the next momentum predictor
519  if (field_relaxation != 1.0)
520  {
521  auto & old_local_solution = *(system.solutionPreviousNewton());
522  NS::FV::relaxSolutionUpdate(current_local_solution, old_local_solution, field_relaxation);
523 
524  // Update old solution, only needed if relaxing the field
525  old_local_solution = current_local_solution;
526  }
527 
528  system.setSolution(current_local_solution);
529 
530  const auto residuals =
531  std::make_pair(its_res_pair.first, linear_solver.get_initial_residual() / norm_factor);
532 
533  _console << " Advected system: " << system.name() << " " << COLOR_GREEN << residuals.second
534  << COLOR_DEFAULT << " Linear its: " << residuals.first << std::endl;
535 
536  return residuals;
537 }
538 
539 bool
541 {
542  // Do not solve if problem is set not to
543  if (!_problem.shouldSolve())
544  return true;
545 
546  // Dummy solver parameter file which is needed for switching petsc options
547  SolverParams solver_params;
548  solver_params._type = Moose::SolveType::ST_LINEAR;
549  solver_params._line_search = Moose::LineSearchType::LS_NONE;
550 
551  // Initialize the SIMPLE iteration counter
552  unsigned int simple_iteration_counter = 0;
553 
554  // Assign residuals to general residual vector
555  const unsigned int no_systems = _momentum_systems.size() + 1 + _has_energy_system +
557  _turbulence_systems.size();
558 
559  std::vector<std::pair<unsigned int, Real>> ns_residuals(no_systems, std::make_pair(0, 1.0));
560  std::vector<Real> ns_abs_tols(_momentum_systems.size(), _momentum_absolute_tolerance);
561  ns_abs_tols.push_back(_pressure_absolute_tolerance);
562 
563  // Push back energy tolerances
564  if (_has_energy_system)
565  ns_abs_tols.push_back(_energy_absolute_tolerance);
567  ns_abs_tols.push_back(_solid_energy_absolute_tolerance);
569  for (const auto scalar_tol : _active_scalar_absolute_tolerance)
570  ns_abs_tols.push_back(scalar_tol);
571 
572  // Push back turbulence tolerances
574  for (const auto turbulence_tol : _turbulence_absolute_tolerance)
575  ns_abs_tols.push_back(turbulence_tol);
576 
577  bool converged = false;
578  // Loop until converged or hit the maximum allowed iteration number
579  while (simple_iteration_counter < _num_iterations && !converged)
580  {
581  simple_iteration_counter++;
582 
583  // We set the preconditioner/controllable parameters through petsc options. Linear
584  // tolerances will be overridden within the solver. In case of a segregated momentum
585  // solver, we assume that every velocity component uses the same preconditioner
587 
588  // Initialize pressure gradients, after this we just reuse the last ones from each
589  // iteration
590  if (simple_iteration_counter == 1)
592 
593  _console << "Iteration " << simple_iteration_counter << " Initial residual norms:" << std::endl;
594 
595  // Solve the momentum predictor step
596  auto momentum_residual = solveMomentumPredictor();
597  for (const auto system_i : index_range(momentum_residual))
598  ns_residuals[system_i] = momentum_residual[system_i];
599 
600  // Now we correct the velocity, this function depends on the method, it differs for
601  // SIMPLE/PIMPLE, this returns the pressure errors
602  ns_residuals[momentum_residual.size()] = correctVelocity(true, true, solver_params);
603 
604  // If we have an energy equation, solve it here.We assume the material properties in the
605  // Navier-Stokes equations depend on temperature, therefore we can not solve for temperature
606  // outside of the velocity-pressure loop
607  if (_has_energy_system)
608  {
609  // We set the preconditioner/controllable parameters through petsc options. Linear
610  // tolerances will be overridden within the solver.
612  ns_residuals[momentum_residual.size() + _has_energy_system] =
618  }
620  {
621  // We set the preconditioner/controllable parameters through petsc options. Linear
622  // tolerances will be overridden within the solver.
624  ns_residuals[momentum_residual.size() + _has_solid_energy_system + _has_energy_system] =
626  }
627 
628  // If we have active scalar equations, solve them here in case they depend on temperature
629  // or they affect the fluid properties such that they must be solved concurrently with pressure
630  // and velocity
632  {
634 
635  // We set the preconditioner/controllable parameters through petsc options. Linear
636  // tolerances will be overridden within the solver.
638  for (const auto i : index_range(_active_scalar_system_names))
639  ns_residuals[momentum_residual.size() + 1 + _has_energy_system + _has_solid_energy_system +
645  }
646 
647  // If we have turbulence equations, solve them here.
648  // The turbulent viscosity depends on the value of the turbulence surrogate variables
650  {
651  // We set the preconditioner/controllable parameters through petsc options. Linear
652  // tolerances will be overridden within the solver.
654  for (const auto i : index_range(_turbulence_system_names))
655  {
656  ns_residuals[momentum_residual.size() + 1 + _has_energy_system + _has_solid_energy_system +
657  _active_scalar_system_names.size() + i] =
665  }
666  }
667 
669 
670  converged = NS::FV::converged(ns_residuals, ns_abs_tols);
671  }
672 
673  // If we have passive scalar equations, solve them here. We assume the material properties in the
674  // Navier-Stokes equations do not depend on passive scalars, as they are passive, therefore we
675  // solve outside of the velocity-pressure loop
677  {
678  // The reason why we need more than one iteration is due to the matrix relaxation
679  // which can be used to stabilize the equations
680  bool passive_scalar_converged = false;
681  unsigned int ps_iteration_counter = 0;
682 
683  _console << "Passive scalar iteration " << ps_iteration_counter
684  << " Initial residual norms:" << std::endl;
685 
686  while (ps_iteration_counter < _num_iterations && !passive_scalar_converged)
687  {
688  ps_iteration_counter++;
689  std::vector<std::pair<unsigned int, Real>> scalar_residuals(
690  _passive_scalar_system_names.size(), std::make_pair(0, 1.0));
691  std::vector<Real> scalar_abs_tols;
692  for (const auto scalar_tol : _passive_scalar_absolute_tolerance)
693  scalar_abs_tols.push_back(scalar_tol);
694 
695  // We set the preconditioner/controllable parameters through petsc options. Linear
696  // tolerances will be overridden within the solver.
698  for (const auto i : index_range(_passive_scalar_system_names))
699  scalar_residuals[i] = solveAdvectedSystem(_passive_scalar_system_numbers[i],
704 
705  passive_scalar_converged = NS::FV::converged(scalar_residuals, scalar_abs_tols);
706  }
707 
708  // Both flow and scalars must converge
709  converged = passive_scalar_converged && converged;
710  }
711 
713 
714  return converged;
715 }
MultiMooseEnum getCommonPetscKeys()
void relaxMatrix(SparseMatrix< Number > &matrix_in, const Real relaxation_parameter, NumericVector< Number > &diff_diagonal)
Relax the matrix to ensure diagonal dominance, we hold onto the difference in diagonals for later use...
bool shouldSolve() const
const std::vector< Real > _active_scalar_equation_relaxation
The user-defined relaxation parameter(s) for the active scalar equation(s)
FEProblemBase & _problem
SIMPLESolverConfiguration _momentum_linear_control
Options for the linear solver of the momentum equation.
SIMPLESolverConfiguration _turbulence_linear_control
Options for the linear solver of the turbulence equation(s)
const Real _energy_absolute_tolerance
The user-defined absolute tolerance for determining the convergence in energy.
const unsigned int _num_iterations
The maximum number of momentum-pressure iterations.
User object responsible for determining the face fluxes using the Rhie-Chow interpolation in a segreg...
SIMPLESolverConfiguration _energy_linear_control
Options for the linear solver of the energy equation.
const unsigned int invalid_uint
void computeGradients()
const bool _print_fields
Debug parameter which allows printing the coupling and solution vectors/matrices. ...
Moose::PetscSupport::PetscOptions _turbulence_petsc_options
Options which hold the petsc settings for the turbulence equation(s)
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
const Real _solid_energy_l_abs_tol
Absolute linear tolerance for the energy equations.
const bool _has_energy_system
Boolean for easy check if a fluid energy system shall be solved or not.
Moose::LineSearchType _line_search
Real computeNormalizationFactor(const NumericVector< Number > &solution, const SparseMatrix< Number > &mat, const NumericVector< Number > &rhs)
Compute a normalization factor which is applied to the linear residual to determine convergence...
SIMPLESolverConfiguration _active_scalar_linear_control
Options for the linear solver of the active scalar equation(s)
std::vector< Real > _turbulence_field_relaxation
The user-defined relaxation parameter(s) for the turbulence field(s)
const Real _passive_scalar_l_abs_tol
Absolute linear tolerance for the passive scalar equation(s).
std::pair< unsigned int, Real > solveSolidEnergy()
Solve an equation which contains the solid energy conservation.
std::vector< LinearSystem * > _passive_scalar_systems
Pointer(s) to the system(s) corresponding to the passive scalar equation(s)
Moose::PetscSupport::PetscOptions _solid_energy_petsc_options
Options which hold the petsc settings for the fluid energy equation.
const Real _momentum_equation_relaxation
The user-defined relaxation parameter for the momentum equation.
virtual std::unique_ptr< NumericVector< T > > zero_clone() const=0
const bool _has_turbulence_systems
Boolean for easy check if a turbulence scalar systems shall be solved or not.
void computeFaceMassFlux()
Update the values of the face velocities in the containers.
std::vector< unsigned int > _turbulence_system_numbers
static InputParameters validParams()
LinearSystem * _solid_energy_system
Pointer to the nonlinear system corresponding to the solid energy equation.
std::pair< unsigned int, Real > solveAdvectedSystem(const unsigned int system_num, LinearSystem &system, const Real relaxation_factor, libMesh::SolverConfiguration &solver_config, const Real abs_tol, const Real field_relaxation=1.0, const Real min_value_limiter=std::numeric_limits< Real >::min())
Solve an equation which contains an advection term that depends on the solution of the segregated Nav...
NumericVector< Number > * rhs
const Real _momentum_l_abs_tol
Absolute linear tolerance for the momentum equation(s).
void setSolution(const NumericVector< Number > &soln)
const bool _has_active_scalar_systems
Boolean for easy check if a active scalar systems shall be solved or not.
virtual std::pair< unsigned int, Real > solve(SparseMatrix< T > &matrix_in, NumericVector< T > &solution_in, NumericVector< T > &rhs_in, const std::optional< double > tol=std::nullopt, const std::optional< unsigned int > m_its=std::nullopt) override
const Real _active_scalar_l_abs_tol
Absolute linear tolerance for the active scalar equation(s).
virtual LinearSolver< Number > * get_linear_solver() const override
RhieChowMassFlux * _rc_uo
Pointer to the segregated RhieChow interpolation object.
The following methods are specializations for using the Parallel::packed_range_* routines for a vecto...
void relaxRightHandSide(NumericVector< Number > &rhs_in, const NumericVector< Number > &solution_in, const NumericVector< Number > &diff_diagonal)
Relax the right hand side of an equation, this needs to be called once and the system matrix has been...
std::map< std::string, Real > real_valued_data
virtual std::pair< unsigned int, Real > correctVelocity(const bool subtract_updated_pressure, const bool recompute_face_mass_flux, const SolverParams &solver_params)
Computes new velocity field based on computed pressure gradients.
bool converged(const std::vector< std::pair< unsigned int, Real >> &residuals, const std::vector< Real > &abs_tolerances)
Based on the residuals, determine if the iterative process converged or not.
const Real _pressure_absolute_tolerance
The user-defined absolute tolerance for determining the convergence in pressure.
const Real _pressure_variable_relaxation
The user-defined relaxation parameter for the pressure variable.
virtual void computeLinearSystemSys(libMesh::LinearImplicitSystem &sys, libMesh::SparseMatrix< libMesh::Number > &system_matrix, NumericVector< libMesh::Number > &rhs, const bool compute_gradients=true)
const Real _momentum_absolute_tolerance
The user-defined absolute tolerance for determining the convergence in momentum.
virtual const std::string & name() const
const std::vector< SolverSystemName > & _passive_scalar_system_names
The names of the passive scalar systems.
const unsigned int _solid_energy_sys_number
The number of the system corresponding to the solid energy equation.
std::vector< LinearSystem * > _active_scalar_systems
Pointer(s) to the system(s) corresponding to the active scalar equation(s)
SIMPLESolverConfiguration _passive_scalar_linear_control
Options for the linear solver of the passive scalar equation(s)
void computeHbyA(const bool with_updated_pressure, const bool verbose)
Computes the inverse of the diagonal (1/A) of the system matrix plus the H/A components for the press...
virtual void execute(const ExecFlagType &exec_type)
Moose::PetscSupport::PetscOptions _energy_petsc_options
Options which hold the petsc settings for the fluid energy equation.
LinearSystem & _pressure_system
Reference to the nonlinear system corresponding to the pressure equation.
std::vector< unsigned int > _momentum_system_numbers
The number(s) of the system(s) corresponding to the momentum equation(s)
SIMPLESolverConfiguration _pressure_linear_control
Options for the linear solver of the pressure equation.
MultiMooseEnum getCommonPetscFlags()
LinearSystem * _energy_system
Pointer to the nonlinear system corresponding to the fluid energy equation.
const unsigned int _energy_sys_number
The number of the system corresponding to the energy equation.
Moose::PetscSupport::PetscOptions _passive_scalar_petsc_options
Options which hold the petsc settings for the passive scalar equation(s)
virtual unsigned int dimension() const
std::map< std::string, int > int_valued_data
virtual std::unique_ptr< SparseMatrix< T > > zero_clone() const=0
void linkMomentumPressureSystems(const std::vector< LinearSystem *> &momentum_systems, const LinearSystem &pressure_system, const std::vector< unsigned int > &momentum_system_numbers)
Update the momentum system-related information.
Moose::PetscSupport::PetscOptions _momentum_petsc_options
Options which hold the petsc settings for the momentum equation.
const Real _pressure_l_abs_tol
Absolute linear tolerance for the pressure equation.
const std::vector< Real > _passive_scalar_absolute_tolerance
The user-defined absolute tolerance for determining the convergence in passive scalars.
std::unique_ptr< NumericVector< Number > > solution
dof_id_type _pressure_pin_dof
The dof ID where the pressure needs to be pinned.
void computeCellVelocity()
Update the cell values of the velocity variables.
Moose::SolveType _type
const std::vector< SolverSystemName > & _turbulence_system_names
The names of the turbulence systems.
virtual bool solve() override
Performs the momentum pressure coupling.
virtual void print(std::ostream &os=libMesh::out) const
const bool _pin_pressure
If the pressure needs to be pinned.
void paramError(const std::string &param, Args... args) const
const Real _energy_l_abs_tol
Absolute linear tolerance for the energy equations.
std::string stringify(const T &t)
void petscSetOptions(const PetscOptions &po, const SolverParams &solver_params, FEProblemBase *const problem=nullptr)
virtual void linkRhieChowUserObject() override
Fetch the Rhie Chow user object that is reponsible for determining face velocities and mass flux...
LinearSystem & getLinearSystem(unsigned int sys_num)
void set_solver_configuration(SolverConfiguration &solver_configuration)
const std::vector< Real > _turbulence_equation_relaxation
The user-defined relaxation parameter(s) for the turbulence equation(s)
const std::vector< SolverSystemName > & _momentum_system_names
The names of the momentum systems.
std::vector< LinearSystem * > _turbulence_systems
Pointer(s) to the system(s) corresponding to the turbulence equation(s)
const ExecFlagType EXEC_NONLINEAR
const std::vector< Real > _passive_scalar_equation_relaxation
The user-defined relaxation parameter(s) for the passive scalar equation(s)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const std::vector< SolverSystemName > & _active_scalar_system_names
The names of the active scalar systems.
MooseApp & _app
void addPetscPairsToPetscOptions(const std::vector< std::pair< MooseEnumItem, std::string >> &petsc_pair_options, const unsigned int mesh_dimension, const std::string &prefix, const ParallelParamObject &param_object, PetscOptions &petsc_options)
const Real _energy_equation_relaxation
The user-defined relaxation parameter for the energy equation.
SparseMatrix< Number > * matrix
void setCurrentLinearSystem(unsigned int sys_num)
std::vector< LinearSystem * > _systems_to_solve
Shortcut to every linear system that we solve for here.
const Real _pressure_pin_value
The value we want to enforce for pressure.
virtual std::vector< std::pair< unsigned int, Real > > solveMomentumPredictor() override
Solve a momentum predictor step with a fixed pressure field.
const Real _turbulence_l_abs_tol
Absolute linear tolerance for the turbulence equation(s).
virtual MooseMesh & mesh() override
unsigned int linearSysNum(const LinearSystemName &linear_sys_name) const override
virtual std::pair< unsigned int, Real > solvePressureCorrector() override
Solve a pressure corrector step.
Moose::PetscSupport::PetscOptions _active_scalar_petsc_options
Options which hold the petsc settings for the active scalar equation(s)
std::unique_ptr< NumericVector< Number > > current_local_solution
virtual const NumericVector< Number > * solutionPreviousNewton() const
Solve class serving as a base class for the two SIMPLE solvers that operate with different assembly a...
const bool _has_passive_scalar_systems
Boolean for easy check if a passive scalar systems shall be solved or not.
const std::vector< Real > _active_scalar_absolute_tolerance
The user-defined absolute tolerance for determining the convergence in active scalars.
void addPetscFlagsToPetscOptions(const MultiMooseEnum &petsc_flags, const std::string &prefix, const ParallelParamObject &param_object, PetscOptions &petsc_options)
void addRangeCheckedParam(const std::string &name, const T &value, const std::string &parsed_function, const std::string &doc_string)
const ConsoleStream _console
const Real _solid_energy_absolute_tolerance
The user-defined absolute tolerance for determining the convergence in solid energy.
std::vector< unsigned int > _active_scalar_system_numbers
SIMPLESolverConfiguration _solid_energy_linear_control
Options for the linear solver of the energy equation.
std::vector< LinearSystem * > _momentum_systems
Pointer(s) to the system(s) corresponding to the momentum equation(s)
const bool _continue_on_max_its
If solve should continue if maximum number of iterations is hit.
std::vector< Real > _turbulence_field_min_limit
The user-defined lower limit for turbulent quantities e.g. k, eps/omega, etc..
void constrainSystem(SparseMatrix< Number > &mx, NumericVector< Number > &rhs, const Real desired_value, const dof_id_type dof_id)
Implicitly constrain the system by adding a factor*(u-u_desired) to it at a desired dof value...
void limitSolutionUpdate(NumericVector< Number > &solution, const Real min_limit=std::numeric_limits< Real >::epsilon(), const Real max_limit=1e10)
Limit a solution to its minimum and maximum bounds: $u = min(max(u, min_limit), max_limit)$.
bool isRecovering() const
void relaxSolutionUpdate(NumericVector< Number > &vec_new, const NumericVector< Number > &vec_old, const Real relaxation_factor)
Relax the update on a solution field using the following approach: $u = u_{old}+ (u - u_{old})$...
const std::vector< Real > _turbulence_absolute_tolerance
The user-defined absolute tolerance for determining the convergence turbulence variables.
auto index_range(const T &sizable)
void print(std::ostream &os=libMesh::out, const bool sparse=false) const
std::vector< unsigned int > _passive_scalar_system_numbers
void initCouplingField()
Initialize the coupling fields (HbyA and Ainv)
Moose::PetscSupport::PetscOptions _pressure_petsc_options
Options which hold the petsc settings for the pressure equation.
virtual System & system() override
const bool _has_solid_energy_system
Boolean for easy check if a solid energy system shall be solved or not.
void initFaceMassFlux()
Initialize the container for face velocities.
void addParamNamesToGroup(const std::string &space_delim_names, const std::string group_name)
const unsigned int _pressure_sys_number
The number of the system corresponding to the pressure equation.