Line data Source code
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 :
10 : #include "LinearAssemblySegregatedSolve.h"
11 : #include "FEProblem.h"
12 : #include "SegregatedSolverUtils.h"
13 : #include "LinearSystem.h"
14 : #include "Executioner.h"
15 :
16 : using namespace libMesh;
17 :
18 : InputParameters
19 1344 : LinearAssemblySegregatedSolve::validParams()
20 : {
21 1344 : InputParameters params = SIMPLESolveBase::validParams();
22 :
23 2688 : params.addParam<std::vector<SolverSystemName>>(
24 : "active_scalar_systems", {}, "The solver system for each active scalar advection equation.");
25 :
26 : /*
27 : * Parameters to control the solution of each scalar advection system
28 : */
29 1344 : params.addParam<std::vector<Real>>("active_scalar_equation_relaxation",
30 1344 : std::vector<Real>(),
31 : "The relaxation which should be used for the active scalar "
32 : "equations. (=1 for no relaxation, "
33 : "diagonal dominance will still be enforced)");
34 :
35 2688 : params.addParam<MultiMooseEnum>("active_scalar_petsc_options",
36 2688 : Moose::PetscSupport::getCommonPetscFlags(),
37 : "Singleton PETSc options for the active scalar equation(s)");
38 2688 : params.addParam<MultiMooseEnum>(
39 : "active_scalar_petsc_options_iname",
40 2688 : Moose::PetscSupport::getCommonPetscKeys(),
41 : "Names of PETSc name/value pairs for the active scalar equation(s)");
42 2688 : params.addParam<std::vector<std::string>>(
43 : "active_scalar_petsc_options_value",
44 : "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\" for the "
45 : "active scalar equation(s)");
46 1344 : params.addParam<std::vector<Real>>(
47 : "active_scalar_absolute_tolerance",
48 1344 : std::vector<Real>(),
49 : "The absolute tolerance(s) on the normalized residual(s) of the active scalar equation(s).");
50 4032 : params.addRangeCheckedParam<Real>("active_scalar_l_tol",
51 2688 : 1e-5,
52 : "0.0<=active_scalar_l_tol & active_scalar_l_tol<1.0",
53 : "The relative tolerance on the normalized residual in the "
54 : "linear solver of the active scalar equation(s).");
55 4032 : params.addRangeCheckedParam<Real>("active_scalar_l_abs_tol",
56 2688 : 1e-10,
57 : "0.0<active_scalar_l_abs_tol",
58 : "The absolute tolerance on the normalized residual in the "
59 : "linear solver of the active scalar equation(s).");
60 2688 : params.addParam<unsigned int>(
61 : "active_scalar_l_max_its",
62 2688 : 10000,
63 : "The maximum allowed iterations in the linear solver of the turbulence equation.");
64 :
65 2688 : params.addParamNamesToGroup(
66 : "active_scalar_systems active_scalar_equation_relaxation active_scalar_petsc_options "
67 : "active_scalar_petsc_options_iname "
68 : "active_scalar_petsc_options_value active_scalar_petsc_options_value "
69 : "active_scalar_absolute_tolerance "
70 : "active_scalar_l_tol active_scalar_l_abs_tol active_scalar_l_max_its",
71 : "Active Scalars Equations");
72 :
73 : /*
74 : * Flags to optionally skip solving subsets of the thermal-hydraulics system (useful when
75 : * recovering a converged solution and only advancing scalar transport for example).
76 : */
77 2688 : params.addParam<bool>(
78 2688 : "should_solve_momentum", true, "Whether we should solve the momentum predictor/corrector.");
79 2688 : params.addParam<bool>(
80 2688 : "should_solve_pressure", true, "Whether we should solve the pressure corrector.");
81 2688 : params.addParam<bool>(
82 2688 : "should_solve_energy", true, "Whether we should solve the fluid energy equation.");
83 2688 : params.addParam<bool>(
84 2688 : "should_solve_solid_energy", true, "Whether we should solve the solid energy equation.");
85 2688 : params.addParam<bool>("should_solve_turbulence",
86 2688 : true,
87 : "Whether we should solve the turbulence surrogate equations.");
88 2688 : params.addParam<bool>(
89 2688 : "should_solve_passive_scalars", true, "Whether we should solve passive scalar equations.");
90 2688 : params.addParam<bool>(
91 2688 : "should_solve_active_scalars", true, "Whether we should solve active scalar equations.");
92 2688 : params.addParam<bool>("should_solve_pm_radiation",
93 2688 : true,
94 : "Whether we should solve participating media radiation equations.");
95 2688 : params.addParamNamesToGroup("should_solve_momentum should_solve_pressure should_solve_energy "
96 : "should_solve_solid_energy should_solve_turbulence "
97 : "should_solve_passive_scalars should_solve_active_scalars",
98 : "Solve control");
99 :
100 : /*
101 : * Parameters to control the conjugate heat transfer
102 : */
103 1344 : params += NS::FV::CHTHandler::validParams();
104 :
105 1344 : return params;
106 0 : }
107 :
108 672 : LinearAssemblySegregatedSolve::LinearAssemblySegregatedSolve(Executioner & ex)
109 : : SIMPLESolveBase(ex),
110 1344 : _pressure_sys_number(_problem.linearSysNum(getParam<SolverSystemName>("pressure_system"))),
111 672 : _pressure_system(_problem.getLinearSystem(_pressure_sys_number)),
112 1344 : _energy_sys_number(_has_energy_system
113 1116 : ? _problem.linearSysNum(getParam<SolverSystemName>("energy_system"))
114 : : libMesh::invalid_uint),
115 672 : _energy_system(_has_energy_system ? &_problem.getLinearSystem(_energy_sys_number) : nullptr),
116 672 : _solid_energy_sys_number(
117 672 : _has_solid_energy_system
118 776 : ? _problem.linearSysNum(getParam<SolverSystemName>("solid_energy_system"))
119 : : libMesh::invalid_uint),
120 672 : _solid_energy_system(
121 672 : _has_solid_energy_system ? &_problem.getLinearSystem(_solid_energy_sys_number) : nullptr),
122 1344 : _should_solve_momentum(getParam<bool>("should_solve_momentum")),
123 1344 : _should_solve_pressure(getParam<bool>("should_solve_pressure")),
124 1344 : _should_solve_energy(getParam<bool>("should_solve_energy")),
125 1344 : _should_solve_solid_energy(getParam<bool>("should_solve_solid_energy")),
126 1344 : _should_solve_turbulence(getParam<bool>("should_solve_turbulence")),
127 1344 : _should_solve_passive_scalars(getParam<bool>("should_solve_passive_scalars")),
128 1344 : _should_solve_active_scalars(getParam<bool>("should_solve_active_scalars")),
129 1344 : _should_solve_pm_radiation(getParam<bool>("should_solve_pm_radiation")),
130 1344 : _active_scalar_system_names(getParam<std::vector<SolverSystemName>>("active_scalar_systems")),
131 672 : _has_active_scalar_systems(!_active_scalar_system_names.empty()),
132 1344 : _active_scalar_equation_relaxation(
133 : getParam<std::vector<Real>>("active_scalar_equation_relaxation")),
134 1344 : _active_scalar_l_abs_tol(getParam<Real>("active_scalar_l_abs_tol")),
135 1344 : _active_scalar_absolute_tolerance(
136 : getParam<std::vector<Real>>("active_scalar_absolute_tolerance")),
137 2016 : _cht(ex.parameters())
138 : {
139 672 : if (!_should_solve_momentum && _should_solve_pressure)
140 0 : paramError("should_solve_momentum",
141 : "Pressure correction requires solving the momentum equations.");
142 672 : if (_should_solve_momentum && !_should_solve_pressure)
143 0 : paramError("should_solve_pressure",
144 : "Solving momentum without a pressure corrector is not supported.");
145 672 : if (_has_solid_energy_system && !_should_solve_energy && _should_solve_solid_energy)
146 0 : paramError("should_solve_solid_energy",
147 : "Solid energy solve cannot be enabled when the fluid energy solve is disabled.");
148 :
149 : // We fetch the systems and their numbers for the momentum equations only if we solve them
150 672 : if (_should_solve_momentum)
151 1970 : for (auto system_i : index_range(_momentum_system_names))
152 : {
153 1298 : _momentum_system_numbers.push_back(_problem.linearSysNum(_momentum_system_names[system_i]));
154 1298 : _momentum_systems.push_back(&_problem.getLinearSystem(_momentum_system_numbers[system_i]));
155 1298 : _systems_to_solve.push_back(_momentum_systems.back());
156 : }
157 :
158 672 : if (_should_solve_pressure)
159 672 : _systems_to_solve.push_back(&_pressure_system);
160 :
161 672 : if (_has_energy_system && _should_solve_energy)
162 222 : _systems_to_solve.push_back(_energy_system);
163 :
164 672 : if (_has_solid_energy_system && _should_solve_solid_energy)
165 52 : _systems_to_solve.push_back(_solid_energy_system);
166 : // and for the turbulence surrogate equations
167 672 : if (_has_turbulence_systems && _should_solve_turbulence)
168 213 : for (auto system_i : index_range(_turbulence_system_names))
169 : {
170 142 : _turbulence_system_numbers.push_back(
171 142 : _problem.linearSysNum(_turbulence_system_names[system_i]));
172 142 : _turbulence_systems.push_back(
173 142 : &_problem.getLinearSystem(_turbulence_system_numbers[system_i]));
174 : }
175 :
176 : // and for the passive scalar equations
177 672 : if (_has_passive_scalar_systems && _should_solve_passive_scalars)
178 114 : for (auto system_i : index_range(_passive_scalar_system_names))
179 : {
180 70 : _passive_scalar_system_numbers.push_back(
181 70 : _problem.linearSysNum(_passive_scalar_system_names[system_i]));
182 70 : _passive_scalar_systems.push_back(
183 70 : &_problem.getLinearSystem(_passive_scalar_system_numbers[system_i]));
184 70 : if (_should_solve_passive_scalars)
185 70 : _systems_to_solve.push_back(_passive_scalar_systems.back());
186 : }
187 :
188 : // and for the participating media radiation equations
189 672 : if (_has_pm_radiation_systems && _should_solve_pm_radiation)
190 32 : for (auto system_i : index_range(_pm_radiation_system_names))
191 : {
192 16 : _pm_radiation_system_numbers.push_back(
193 16 : _problem.linearSysNum(_pm_radiation_system_names[system_i]));
194 16 : _pm_radiation_systems.push_back(
195 16 : &_problem.getLinearSystem(_pm_radiation_system_numbers[system_i]));
196 16 : _systems_to_solve.push_back(_pm_radiation_systems.back());
197 : }
198 :
199 : // and for the active scalar equations
200 672 : if (_has_active_scalar_systems && _should_solve_active_scalars)
201 50 : for (auto system_i : index_range(_active_scalar_system_names))
202 : {
203 25 : _active_scalar_system_numbers.push_back(
204 25 : _problem.linearSysNum(_active_scalar_system_names[system_i]));
205 25 : _active_scalar_systems.push_back(
206 25 : &_problem.getLinearSystem(_active_scalar_system_numbers[system_i]));
207 25 : _systems_to_solve.push_back(_active_scalar_systems.back());
208 :
209 : const auto & active_scalar_petsc_options =
210 25 : getParam<MultiMooseEnum>("active_scalar_petsc_options");
211 : const auto & active_scalar_petsc_pair_options = getParam<MooseEnumItem, std::string>(
212 50 : "active_scalar_petsc_options_iname", "active_scalar_petsc_options_value");
213 25 : Moose::PetscSupport::addPetscFlagsToPetscOptions(
214 : active_scalar_petsc_options, "", *this, _active_scalar_petsc_options);
215 50 : Moose::PetscSupport::addPetscPairsToPetscOptions(active_scalar_petsc_pair_options,
216 25 : _problem.mesh().dimension(),
217 : "",
218 : *this,
219 : _active_scalar_petsc_options);
220 :
221 25 : _active_scalar_linear_control.real_valued_data["rel_tol"] =
222 50 : getParam<Real>("active_scalar_l_tol");
223 25 : _active_scalar_linear_control.real_valued_data["abs_tol"] =
224 50 : getParam<Real>("active_scalar_l_abs_tol");
225 25 : _active_scalar_linear_control.int_valued_data["max_its"] =
226 50 : getParam<unsigned int>("active_scalar_l_max_its");
227 25 : }
228 :
229 672 : if (_active_scalar_equation_relaxation.size() != _active_scalar_system_names.size())
230 0 : paramError("active_scalar_equation_relaxation",
231 : "Should be the same size as the number of systems");
232 :
233 : // We disable the prefix here for the time being, the segregated solvers use a different approach
234 : // for setting the petsc parameters
235 3027 : for (auto & system : _systems_to_solve)
236 2355 : system->system().prefix_with_name(false);
237 :
238 : // Link CHT objects, this will also do some error checking
239 : // Make a copy for compatibility. These could change in the future
240 : // Convert _pm_radiation_systems to std::vector<SystemBase *>
241 672 : if (_cht.enabled())
242 : {
243 40 : if (!_should_solve_energy)
244 0 : paramError("should_solve_energy",
245 : "Conjugate heat transfer requires solving the fluid energy equation.");
246 40 : if (_has_solid_energy_system && !_should_solve_solid_energy)
247 0 : paramError("should_solve_solid_energy",
248 : "Conjugate heat transfer requires solving the solid energy equation.");
249 :
250 : std::vector<SystemBase *> pm_radiation_systems_base(_pm_radiation_systems.begin(),
251 40 : _pm_radiation_systems.end());
252 :
253 40 : _cht.linkEnergySystems(_solid_energy_system, _energy_system, pm_radiation_systems_base);
254 40 : }
255 672 : }
256 :
257 : void
258 671 : LinearAssemblySegregatedSolve::linkRhieChowUserObject()
259 : {
260 671 : if (!_should_solve_momentum)
261 : return;
262 :
263 671 : _rc_uo =
264 671 : const_cast<RhieChowMassFlux *>(&getUserObject<RhieChowMassFlux>("rhie_chow_user_object"));
265 671 : _rc_uo->linkMomentumPressureSystems(
266 671 : _momentum_systems, _pressure_system, _momentum_system_numbers);
267 :
268 : // Initialize the face velocities in the RC object
269 671 : if (!_app.isRecovering())
270 652 : _rc_uo->initFaceMassFlux();
271 671 : _rc_uo->initCouplingField();
272 : }
273 :
274 : std::vector<std::pair<unsigned int, Real>>
275 119163 : LinearAssemblySegregatedSolve::solveMomentumPredictor()
276 : {
277 : // Temporary storage for the (flux-normalized) residuals from
278 : // different momentum components
279 : std::vector<std::pair<unsigned int, Real>> its_normalized_residuals;
280 :
281 : LinearImplicitSystem & momentum_system_0 =
282 119163 : libMesh::cast_ref<LinearImplicitSystem &>(_momentum_systems[0]->system());
283 :
284 : PetscLinearSolver<Real> & momentum_solver =
285 119163 : libMesh::cast_ref<PetscLinearSolver<Real> &>(*momentum_system_0.get_linear_solver());
286 :
287 : // Solve the momentum equations.
288 : // TO DO: These equations are VERY similar. If we can store the differences (things coming from
289 : // BCs for example) separately, it is enough to construct one matrix.
290 321006 : for (const auto system_i : index_range(_momentum_systems))
291 : {
292 201843 : _problem.setCurrentLinearSystem(_momentum_system_numbers[system_i]);
293 :
294 : // We will need the right hand side and the solution of the next component
295 : LinearImplicitSystem & momentum_system =
296 201843 : libMesh::cast_ref<LinearImplicitSystem &>(_momentum_systems[system_i]->system());
297 :
298 : NumericVector<Number> & solution = *(momentum_system.solution);
299 201843 : NumericVector<Number> & rhs = *(momentum_system.rhs);
300 201843 : SparseMatrix<Number> & mmat = *(momentum_system.matrix);
301 :
302 201843 : auto diff_diagonal = solution.zero_clone();
303 :
304 : // We assemble the matrix and the right hand side
305 201843 : _problem.computeLinearSystemSys(momentum_system, mmat, rhs, /*compute_grads*/ true);
306 :
307 : // Still need to relax the right hand side with the same vector
308 201843 : NS::FV::relaxMatrix(mmat, _momentum_equation_relaxation, *diff_diagonal);
309 201843 : NS::FV::relaxRightHandSide(rhs, solution, *diff_diagonal);
310 :
311 : // The normalization factor depends on the right hand side so we need to recompute it for this
312 : // component
313 201843 : Real norm_factor = NS::FV::computeNormalizationFactor(solution, mmat, rhs);
314 :
315 : // Very important, for deciding the convergence, we need the unpreconditioned
316 : // norms in the linear solve
317 201843 : LibmeshPetscCall(KSPSetNormType(momentum_solver.ksp(), KSP_NORM_UNPRECONDITIONED));
318 : // Solve this component. We don't update the ghosted solution yet, that will come at the end
319 : // of the corrector step. Also setting the linear tolerances and maximum iteration counts.
320 201843 : _momentum_linear_control.real_valued_data["abs_tol"] = _momentum_l_abs_tol * norm_factor;
321 201843 : momentum_solver.set_solver_configuration(_momentum_linear_control);
322 :
323 : // We solve the equation
324 201843 : auto its_resid_pair = momentum_solver.solve(mmat, mmat, solution, rhs);
325 201843 : momentum_system.update();
326 :
327 : // We will reuse the preconditioner for every momentum system
328 201843 : if (system_i == 0)
329 119163 : momentum_solver.reuse_preconditioner(true);
330 :
331 : // Save the normalized residual
332 : its_normalized_residuals.push_back(
333 201843 : std::make_pair(its_resid_pair.first, momentum_solver.get_initial_residual() / norm_factor));
334 :
335 201843 : if (_print_fields)
336 : {
337 0 : _console << " matrix when we solve " << std::endl;
338 0 : mmat.print();
339 0 : _console << " rhs when we solve " << std::endl;
340 0 : rhs.print();
341 0 : _console << " velocity solution component " << system_i << std::endl;
342 0 : solution.print();
343 0 : _console << "Norm factor " << norm_factor << std::endl;
344 0 : _console << Moose::stringify(momentum_solver.get_initial_residual()) << std::endl;
345 : }
346 :
347 : // Printing residuals
348 201843 : _console << " Momentum equation:"
349 : << (_momentum_systems.size() > 1
350 897516 : ? std::string(" Component ") + std::to_string(system_i + 1) + std::string(" ")
351 201843 : : std::string(" "))
352 201843 : << COLOR_GREEN << its_normalized_residuals[system_i].second << COLOR_DEFAULT
353 201843 : << " Linear its: " << its_normalized_residuals[system_i].first << std::endl;
354 201843 : }
355 :
356 321006 : for (const auto system_i : index_range(_momentum_systems))
357 : {
358 : LinearImplicitSystem & momentum_system =
359 201843 : libMesh::cast_ref<LinearImplicitSystem &>(_momentum_systems[system_i]->system());
360 201843 : _momentum_systems[system_i]->setSolution(*(momentum_system.current_local_solution));
361 201843 : _momentum_systems[system_i]->copyPreviousNonlinearSolutions();
362 : }
363 :
364 : // We reset this to ensure the preconditioner is recomputed new time we go to the momentum
365 : // predictor
366 119163 : momentum_solver.reuse_preconditioner(false);
367 :
368 119163 : return its_normalized_residuals;
369 0 : }
370 :
371 : void
372 671 : LinearAssemblySegregatedSolve::initialSetup()
373 : {
374 671 : if (_cht.enabled())
375 : {
376 40 : _cht.deduceCHTBoundaryCoupling();
377 40 : _cht.setupConjugateHeatTransferContainers();
378 : }
379 671 : }
380 :
381 : std::pair<unsigned int, Real>
382 132163 : LinearAssemblySegregatedSolve::solvePressureCorrector()
383 : {
384 132163 : _problem.setCurrentLinearSystem(_pressure_sys_number);
385 :
386 : // We will need some members from the linear system
387 : LinearImplicitSystem & pressure_system =
388 132163 : libMesh::cast_ref<LinearImplicitSystem &>(_pressure_system.system());
389 :
390 : // We will need the solution, the right hand side and the matrix
391 : NumericVector<Number> & current_local_solution = *(pressure_system.current_local_solution);
392 : NumericVector<Number> & solution = *(pressure_system.solution);
393 132163 : SparseMatrix<Number> & mmat = *(pressure_system.matrix);
394 132163 : NumericVector<Number> & rhs = *(pressure_system.rhs);
395 :
396 : // Fetch the linear solver from the system
397 : PetscLinearSolver<Real> & pressure_solver =
398 132163 : libMesh::cast_ref<PetscLinearSolver<Real> &>(*pressure_system.get_linear_solver());
399 :
400 132163 : _problem.computeLinearSystemSys(pressure_system, mmat, rhs, false);
401 :
402 132163 : if (_print_fields)
403 : {
404 0 : _console << "Pressure matrix" << std::endl;
405 0 : mmat.print();
406 : }
407 :
408 : // We compute the normalization factors based on the fluxes
409 132163 : Real norm_factor = NS::FV::computeNormalizationFactor(solution, mmat, rhs);
410 :
411 : // We need the non-preconditioned norm to be consistent with the norm factor
412 132163 : LibmeshPetscCall(KSPSetNormType(pressure_solver.ksp(), KSP_NORM_UNPRECONDITIONED));
413 :
414 : // Setting the linear tolerances and maximum iteration counts
415 132163 : _pressure_linear_control.real_valued_data["abs_tol"] = _pressure_l_abs_tol * norm_factor;
416 132163 : pressure_solver.set_solver_configuration(_pressure_linear_control);
417 :
418 132163 : if (_pin_pressure)
419 27901 : NS::FV::constrainSystem(mmat, rhs, _pressure_pin_value, _pressure_pin_dof);
420 132163 : pressure_system.update();
421 :
422 132163 : auto its_res_pair = pressure_solver.solve(mmat, mmat, solution, rhs);
423 132163 : pressure_system.update();
424 :
425 132163 : if (_print_fields)
426 : {
427 0 : _console << " rhs when we solve pressure " << std::endl;
428 0 : rhs.print();
429 0 : _console << " Pressure " << std::endl;
430 0 : solution.print();
431 0 : _console << "Norm factor " << norm_factor << std::endl;
432 : }
433 :
434 132163 : _pressure_system.setSolution(current_local_solution);
435 :
436 : const auto residuals =
437 132163 : std::make_pair(its_res_pair.first, pressure_solver.get_initial_residual() / norm_factor);
438 :
439 132163 : _console << " Pressure equation: " << COLOR_GREEN << residuals.second << COLOR_DEFAULT
440 132163 : << " Linear its: " << residuals.first << std::endl;
441 :
442 132163 : return residuals;
443 : }
444 :
445 : std::pair<unsigned int, Real>
446 10642 : LinearAssemblySegregatedSolve::solveSolidEnergy()
447 : {
448 10642 : _problem.setCurrentLinearSystem(_solid_energy_sys_number);
449 :
450 : // We will need some members from the linear system
451 : LinearImplicitSystem & system =
452 10642 : libMesh::cast_ref<LinearImplicitSystem &>(_solid_energy_system->system());
453 :
454 : // We will need the solution, the right hand side and the matrix
455 : NumericVector<Number> & current_local_solution = *(system.current_local_solution);
456 : NumericVector<Number> & solution = *(system.solution);
457 10642 : SparseMatrix<Number> & mmat = *(system.matrix);
458 10642 : NumericVector<Number> & rhs = *(system.rhs);
459 :
460 : // Fetch the linear solver from the system
461 : PetscLinearSolver<Real> & solver =
462 10642 : libMesh::cast_ref<PetscLinearSolver<Real> &>(*system.get_linear_solver());
463 :
464 10642 : _problem.computeLinearSystemSys(system, mmat, rhs, false);
465 :
466 10642 : if (_print_fields)
467 : {
468 0 : _console << "Solid energy matrix" << std::endl;
469 0 : mmat.print();
470 : }
471 :
472 : // We compute the normalization factors based on the fluxes
473 10642 : Real norm_factor = NS::FV::computeNormalizationFactor(solution, mmat, rhs);
474 :
475 : // We need the non-preconditioned norm to be consistent with the norm factor
476 10642 : LibmeshPetscCall(KSPSetNormType(solver.ksp(), KSP_NORM_UNPRECONDITIONED));
477 :
478 : // Setting the linear tolerances and maximum iteration counts
479 10642 : _solid_energy_linear_control.real_valued_data["abs_tol"] = _solid_energy_l_abs_tol * norm_factor;
480 10642 : solver.set_solver_configuration(_solid_energy_linear_control);
481 :
482 10642 : auto its_res_pair = solver.solve(mmat, mmat, solution, rhs);
483 10642 : system.update();
484 :
485 10642 : if (_print_fields)
486 : {
487 0 : _console << " rhs when we solve solid energy " << std::endl;
488 0 : rhs.print();
489 0 : _console << " Solid energy " << std::endl;
490 0 : solution.print();
491 0 : _console << "Norm factor " << norm_factor << std::endl;
492 : }
493 :
494 10642 : _solid_energy_system->setSolution(current_local_solution);
495 :
496 : const auto residuals =
497 10642 : std::make_pair(its_res_pair.first, solver.get_initial_residual() / norm_factor);
498 :
499 10642 : _console << " Solid energy equation: " << COLOR_GREEN << residuals.second << COLOR_DEFAULT
500 10642 : << " Linear its: " << residuals.first << std::endl;
501 :
502 10642 : return residuals;
503 : }
504 :
505 : std::pair<unsigned int, Real>
506 132163 : LinearAssemblySegregatedSolve::correctVelocity(const bool subtract_updated_pressure,
507 : const bool recompute_face_mass_flux,
508 : const SolverParams & solver_params)
509 : {
510 : // Compute the coupling fields between the momentum and pressure equations.
511 : // The first argument makes sure the pressure gradient is staged at the first
512 : // iteration
513 132163 : _rc_uo->computeHbyA(subtract_updated_pressure, _print_fields);
514 :
515 : // We set the preconditioner/controllable parameters for the pressure equations through
516 : // petsc options. Linear tolerances will be overridden within the solver.
517 132163 : Moose::PetscSupport::petscSetOptions(_pressure_petsc_options, solver_params);
518 :
519 : // Solve the pressure corrector
520 132163 : const auto residuals = solvePressureCorrector();
521 :
522 : // Compute the face velocity which is used in the advection terms. In certain
523 : // segregated solver algorithms (like PISO) this is only done on the last iteration.
524 132163 : if (recompute_face_mass_flux)
525 119163 : _rc_uo->computeFaceMassFlux();
526 :
527 132163 : auto & pressure_current_solution = *(_pressure_system.system().current_local_solution.get());
528 132163 : auto & pressure_old_solution = *(_pressure_system.solutionPreviousNewton());
529 :
530 : // Relax the pressure update for the next momentum predictor
531 132163 : NS::FV::relaxSolutionUpdate(
532 132163 : pressure_current_solution, pressure_old_solution, _pressure_variable_relaxation);
533 :
534 : // Overwrite old solution
535 132163 : pressure_old_solution = pressure_current_solution;
536 132163 : _pressure_system.setSolution(pressure_current_solution);
537 :
538 : // We recompute the updated pressure gradient
539 132163 : _pressure_system.computeGradients();
540 :
541 : // Reconstruct the cell velocity as well to accelerate convergence
542 132163 : _rc_uo->computeCellVelocity();
543 :
544 132163 : return residuals;
545 : }
546 :
547 : std::pair<unsigned int, Real>
548 107388 : LinearAssemblySegregatedSolve::solveAdvectedSystem(const unsigned int system_num,
549 : LinearSystem & system,
550 : const Real relaxation_factor,
551 : SolverConfiguration & solver_config,
552 : const Real absolute_tol,
553 : const Real field_relaxation,
554 : const Real min_value_limiter)
555 : {
556 107388 : _problem.setCurrentLinearSystem(system_num);
557 :
558 : // We will need some members from the implicit linear system
559 107388 : LinearImplicitSystem & li_system = libMesh::cast_ref<LinearImplicitSystem &>(system.system());
560 :
561 : // We will need the solution, the right hand side and the matrix
562 : NumericVector<Number> & current_local_solution = *(li_system.current_local_solution);
563 : NumericVector<Number> & solution = *(li_system.solution);
564 107388 : SparseMatrix<Number> & mmat = *(li_system.matrix);
565 107388 : NumericVector<Number> & rhs = *(li_system.rhs);
566 :
567 : // We need a vector that stores the (diagonal_relaxed-original_diagonal) vector
568 107388 : auto diff_diagonal = solution.zero_clone();
569 :
570 : // Fetch the linear solver from the system
571 : PetscLinearSolver<Real> & linear_solver =
572 107388 : libMesh::cast_ref<PetscLinearSolver<Real> &>(*li_system.get_linear_solver());
573 :
574 107388 : _problem.computeLinearSystemSys(li_system, mmat, rhs, true);
575 :
576 : // Go and relax the system matrix and the right hand side
577 107388 : NS::FV::relaxMatrix(mmat, relaxation_factor, *diff_diagonal);
578 107388 : NS::FV::relaxRightHandSide(rhs, solution, *diff_diagonal);
579 :
580 107388 : if (_print_fields)
581 : {
582 0 : _console << system.name() << " system matrix" << std::endl;
583 0 : mmat.print();
584 : }
585 :
586 : // We compute the normalization factors based on the fluxes
587 107388 : Real norm_factor = NS::FV::computeNormalizationFactor(solution, mmat, rhs);
588 :
589 : // We need the non-preconditioned norm to be consistent with the norm factor
590 107388 : LibmeshPetscCall(KSPSetNormType(linear_solver.ksp(), KSP_NORM_UNPRECONDITIONED));
591 :
592 : // Setting the linear tolerances and maximum iteration counts
593 107388 : solver_config.real_valued_data["abs_tol"] = absolute_tol * norm_factor;
594 107388 : linear_solver.set_solver_configuration(solver_config);
595 :
596 : // Solve the system and update current local solution
597 107388 : auto its_res_pair = linear_solver.solve(mmat, mmat, solution, rhs);
598 107388 : li_system.update();
599 :
600 107388 : if (_print_fields)
601 : {
602 0 : _console << " rhs when we solve " << system.name() << std::endl;
603 0 : rhs.print();
604 0 : _console << system.name() << " solution " << std::endl;
605 0 : solution.print();
606 0 : _console << " Norm factor " << norm_factor << std::endl;
607 : }
608 :
609 : // Limiting scalar solution
610 107388 : if (min_value_limiter != std::numeric_limits<Real>::min())
611 52996 : NS::FV::limitSolutionUpdate(current_local_solution, min_value_limiter);
612 :
613 : // Relax the field update for the next momentum predictor
614 107388 : if (field_relaxation != 1.0)
615 : {
616 35652 : auto & old_local_solution = *(system.solutionPreviousNewton());
617 35652 : NS::FV::relaxSolutionUpdate(current_local_solution, old_local_solution, field_relaxation);
618 :
619 : // Update old solution, only needed if relaxing the field
620 35652 : old_local_solution = current_local_solution;
621 : }
622 :
623 107388 : system.setSolution(current_local_solution);
624 :
625 : const auto residuals =
626 107388 : std::make_pair(its_res_pair.first, linear_solver.get_initial_residual() / norm_factor);
627 :
628 107388 : _console << " Advected system: " << system.name() << " " << COLOR_GREEN << residuals.second
629 107388 : << COLOR_DEFAULT << " Linear its: " << residuals.first << std::endl;
630 :
631 107388 : return residuals;
632 107388 : }
633 :
634 : bool
635 1183 : LinearAssemblySegregatedSolve::solve()
636 : {
637 : // Do not solve if problem is set not to
638 1183 : if (!_problem.shouldSolve())
639 : return true;
640 :
641 : // Dummy solver parameter file which is needed for switching petsc options
642 1183 : SolverParams solver_params;
643 1183 : solver_params._type = Moose::SolveType::ST_LINEAR;
644 1183 : solver_params._line_search = Moose::LineSearchType::LS_NONE;
645 :
646 : // Initialize the SIMPLE iteration counter
647 1183 : unsigned int simple_iteration_counter = 0;
648 :
649 : // We set up the residual storage and the corresponding tolerances.
650 1183 : ResidualStorage residual_storage = setupResidualStorage();
651 : auto & ns_residuals = residual_storage.ns_residuals;
652 : auto & ns_abs_tols = residual_storage.ns_abs_tols;
653 : const auto & momentum_indices = residual_storage.momentum_indices;
654 1183 : const auto pressure_index = residual_storage.pressure_index;
655 1183 : const auto energy_index = residual_storage.energy_index;
656 1183 : const auto solid_energy_index = residual_storage.solid_energy_index;
657 : const auto & active_scalar_indices = residual_storage.active_scalar_indices;
658 : const auto & turbulence_indices = residual_storage.turbulence_indices;
659 : const auto & pm_radiation_indices = residual_storage.pm_radiation_indices;
660 :
661 1183 : bool converged = residual_storage.converged;
662 :
663 : // Loop until converged or hit the maximum allowed iteration number
664 1183 : if (_cht.enabled() && _should_solve_energy)
665 38 : _cht.initializeCHTCouplingFields();
666 :
667 120346 : while (simple_iteration_counter < _num_iterations && !converged)
668 : {
669 119163 : simple_iteration_counter++;
670 :
671 : // We set the preconditioner/controllable parameters through petsc options. Linear
672 : // tolerances will be overridden within the solver. In case of a segregated momentum
673 : // solver, we assume that every velocity component uses the same preconditioner
674 119163 : if (_should_solve_momentum)
675 119163 : Moose::PetscSupport::petscSetOptions(_momentum_petsc_options, solver_params);
676 :
677 : // Initialize pressure gradients, after this we just reuse the last ones from each
678 : // iteration
679 119163 : if (_should_solve_pressure && simple_iteration_counter == 1)
680 1183 : _pressure_system.computeGradients();
681 :
682 119163 : _console << "Iteration " << simple_iteration_counter << " Initial residual norms:" << std::endl;
683 :
684 : // Solve the momentum predictor step
685 119163 : if (_should_solve_momentum)
686 : {
687 119163 : auto momentum_residual = solveMomentumPredictor();
688 321006 : for (const auto system_i : index_range(momentum_residual))
689 201843 : ns_residuals[momentum_indices[system_i]] = momentum_residual[system_i];
690 119163 : }
691 :
692 : // Now we correct the velocity, this function depends on the method, it differs for
693 : // SIMPLE/PIMPLE, this returns the pressure errors
694 119163 : if (_should_solve_pressure)
695 119163 : ns_residuals[pressure_index] = correctVelocity(true, true, solver_params);
696 :
697 : // If we have an energy equation, solve it here.We assume the material properties in the
698 : // Navier-Stokes equations depend on temperature, therefore we can not solve for temperature
699 : // outside of the velocity-pressure loop
700 119163 : if (_has_energy_system && _should_solve_energy)
701 : {
702 : // If there is no CHT specified this will just do go once through this block
703 : _cht.resetCHTConvergence();
704 56908 : while (!_cht.converged())
705 : {
706 29806 : if (_cht.enabled())
707 7276 : _cht.updateCHTBoundaryCouplingFields(NS::CHTSide::FLUID);
708 :
709 : // We set the preconditioner/controllable parameters through petsc options. Linear
710 : // tolerances will be overridden within the solver.
711 29806 : Moose::PetscSupport::petscSetOptions(_energy_petsc_options, solver_params);
712 59612 : ns_residuals[energy_index] = solveAdvectedSystem(_energy_sys_number,
713 29806 : *_energy_system,
714 29806 : _energy_equation_relaxation,
715 : _energy_linear_control,
716 29806 : _energy_l_abs_tol);
717 :
718 29806 : if (_has_pm_radiation_systems && _should_solve_pm_radiation)
719 : {
720 : // We set the preconditioner/controllable parameters through petsc options. Linear
721 : // tolerances will be overridden within the solver.
722 2212 : Moose::PetscSupport::petscSetOptions(_pm_radiation_petsc_options, solver_params);
723 4424 : for (const auto i : index_range(_pm_radiation_system_names))
724 : {
725 2212 : ns_residuals[pm_radiation_indices[i]] =
726 2212 : solveAdvectedSystem(_pm_radiation_system_numbers[i],
727 : *_pm_radiation_systems[i],
728 : _pm_radiation_equation_relaxation[i],
729 : _pm_radiation_linear_control,
730 2212 : _pm_radiation_l_abs_tol);
731 : }
732 : }
733 :
734 29806 : if (_has_solid_energy_system && _should_solve_solid_energy)
735 : {
736 : // For now we only update gradients if cht is needed, might change in the future
737 10642 : if (_cht.enabled())
738 : {
739 7276 : _energy_system->computeGradients();
740 7276 : _cht.updateCHTBoundaryCouplingFields(NS::CHTSide::SOLID);
741 : }
742 :
743 : // We set the preconditioner/controllable parameters through petsc options. Linear
744 : // tolerances will be overridden within the solver.
745 10642 : Moose::PetscSupport::petscSetOptions(_solid_energy_petsc_options, solver_params);
746 10642 : ns_residuals[solid_energy_index] = solveSolidEnergy();
747 :
748 : // For now we only update gradients if cht is needed, might change in the future
749 10642 : if (_cht.enabled())
750 7276 : _solid_energy_system->computeGradients();
751 : }
752 :
753 29806 : if (_cht.enabled())
754 : {
755 7276 : _cht.sumIntegratedFluxes();
756 7276 : _cht.printIntegratedFluxes();
757 : }
758 :
759 : _cht.incrementCHTIterators();
760 : }
761 27102 : if (_cht.enabled())
762 4572 : _cht.resetIntegratedFluxes();
763 : }
764 :
765 : // If we have active scalar equations, solve them here in case they depend on temperature
766 : // or they affect the fluid properties such that they must be solved concurrently with
767 : // pressure and velocity
768 119163 : if (_has_active_scalar_systems && _should_solve_active_scalars)
769 : {
770 12088 : _problem.execute(EXEC_NONLINEAR);
771 :
772 : // We set the preconditioner/controllable parameters through petsc options. Linear
773 : // tolerances will be overridden within the solver.
774 12088 : Moose::PetscSupport::petscSetOptions(_active_scalar_petsc_options, solver_params);
775 24176 : for (const auto i : index_range(_active_scalar_system_names))
776 12088 : ns_residuals[active_scalar_indices[i]] =
777 12088 : solveAdvectedSystem(_active_scalar_system_numbers[i],
778 : *_active_scalar_systems[i],
779 : _active_scalar_equation_relaxation[i],
780 : _active_scalar_linear_control,
781 12088 : _active_scalar_l_abs_tol);
782 : }
783 :
784 : // If we have turbulence equations, solve them here.
785 : // The turbulent viscosity depends on the value of the turbulence surrogate variables
786 119163 : if (_has_turbulence_systems && _should_solve_turbulence)
787 : {
788 : // We set the preconditioner/controllable parameters through petsc options. Linear
789 : // tolerances will be overridden within the solver.
790 26498 : Moose::PetscSupport::petscSetOptions(_turbulence_petsc_options, solver_params);
791 79494 : for (const auto i : index_range(_turbulence_system_names))
792 : {
793 52996 : ns_residuals[turbulence_indices[i]] =
794 52996 : solveAdvectedSystem(_turbulence_system_numbers[i],
795 : *_turbulence_systems[i],
796 : _turbulence_equation_relaxation[i],
797 : _turbulence_linear_control,
798 52996 : _turbulence_l_abs_tol,
799 : _turbulence_field_relaxation[i],
800 : _turbulence_field_min_limit[i]);
801 : }
802 : }
803 :
804 119163 : _problem.execute(EXEC_NONLINEAR);
805 :
806 119163 : converged = NS::FV::converged(ns_residuals, ns_abs_tols);
807 : }
808 :
809 : // If we have passive scalar equations, solve them here. We assume the material properties in
810 : // the Navier-Stokes equations do not depend on passive scalars, as they are passive, therefore
811 : // we solve outside of the velocity-pressure loop
812 1183 : if (_has_passive_scalar_systems && _should_solve_passive_scalars &&
813 47 : (converged || _continue_on_max_its))
814 : {
815 : // The reason why we need more than one iteration is due to the matrix relaxation
816 : // which can be used to stabilize the equations
817 : bool passive_scalar_converged = false;
818 72 : unsigned int ps_iteration_counter = 0;
819 :
820 72 : _console << "Passive scalar iteration " << ps_iteration_counter
821 72 : << " Initial residual norms:" << std::endl;
822 :
823 5299 : while (ps_iteration_counter < _num_iterations && !passive_scalar_converged)
824 : {
825 5227 : ps_iteration_counter++;
826 : std::vector<std::pair<unsigned int, Real>> scalar_residuals(
827 5227 : _passive_scalar_system_names.size(), std::make_pair(0, 1.0));
828 : std::vector<Real> scalar_abs_tols;
829 15513 : for (const auto scalar_tol : _passive_scalar_absolute_tolerance)
830 10286 : scalar_abs_tols.push_back(scalar_tol);
831 :
832 : // We set the preconditioner/controllable parameters through petsc options. Linear
833 : // tolerances will be overridden within the solver.
834 5227 : Moose::PetscSupport::petscSetOptions(_passive_scalar_petsc_options, solver_params);
835 15513 : for (const auto i : index_range(_passive_scalar_system_names))
836 10286 : scalar_residuals[i] = solveAdvectedSystem(_passive_scalar_system_numbers[i],
837 : *_passive_scalar_systems[i],
838 : _passive_scalar_equation_relaxation[i],
839 : _passive_scalar_linear_control,
840 10286 : _passive_scalar_l_abs_tol);
841 :
842 5227 : passive_scalar_converged = NS::FV::converged(scalar_residuals, scalar_abs_tols);
843 5227 : }
844 :
845 : // Both flow and scalars must converge
846 72 : converged = passive_scalar_converged && converged;
847 : }
848 :
849 1183 : converged = _continue_on_max_its ? true : converged;
850 :
851 : return converged;
852 1183 : }
853 :
854 : LinearAssemblySegregatedSolve::ResidualStorage
855 1183 : LinearAssemblySegregatedSolve::setupResidualStorage() const
856 : {
857 1183 : ResidualStorage storage;
858 :
859 : // Residual store: position in this vector defines the ordering used by NS::FV::converged()
860 : // Each entry holds (linear its, normalized residual) for one system
861 1183 : if (_should_solve_momentum)
862 3479 : for ([[maybe_unused]] const auto system_i : index_range(_momentum_systems))
863 : {
864 2296 : storage.momentum_indices.push_back(storage.ns_residuals.size());
865 2296 : storage.ns_residuals.push_back(std::make_pair(0, 1.0));
866 2296 : storage.ns_abs_tols.push_back(_momentum_absolute_tolerance);
867 : }
868 :
869 1183 : if (_should_solve_pressure)
870 : {
871 1183 : storage.pressure_index = storage.ns_residuals.size();
872 1183 : storage.ns_residuals.push_back(std::make_pair(0, 1.0));
873 1183 : storage.ns_abs_tols.push_back(_pressure_absolute_tolerance);
874 : }
875 :
876 1183 : if (_has_energy_system && _should_solve_energy)
877 : {
878 590 : storage.energy_index = storage.ns_residuals.size();
879 590 : storage.ns_residuals.push_back(std::make_pair(0, 1.0));
880 590 : storage.ns_abs_tols.push_back(_energy_absolute_tolerance);
881 : }
882 :
883 1183 : if (_has_solid_energy_system && _should_solve_solid_energy)
884 : {
885 50 : storage.solid_energy_index = storage.ns_residuals.size();
886 50 : storage.ns_residuals.push_back(std::make_pair(0, 1.0));
887 50 : storage.ns_abs_tols.push_back(_solid_energy_absolute_tolerance);
888 : }
889 :
890 1183 : if (_has_active_scalar_systems && _should_solve_active_scalars)
891 170 : for (const auto i : index_range(_active_scalar_system_names))
892 : {
893 85 : storage.active_scalar_indices.push_back(storage.ns_residuals.size());
894 85 : storage.ns_residuals.push_back(std::make_pair(0, 1.0));
895 85 : storage.ns_abs_tols.push_back(_active_scalar_absolute_tolerance[i]);
896 : }
897 :
898 1183 : if (_has_turbulence_systems && _should_solve_turbulence)
899 369 : for (const auto i : index_range(_turbulence_system_names))
900 : {
901 246 : storage.turbulence_indices.push_back(storage.ns_residuals.size());
902 246 : storage.ns_residuals.push_back(std::make_pair(0, 1.0));
903 246 : storage.ns_abs_tols.push_back(_turbulence_absolute_tolerance[i]);
904 : }
905 :
906 1183 : if (_has_pm_radiation_systems && _should_solve_pm_radiation)
907 28 : for (const auto i : index_range(_pm_radiation_system_names))
908 : {
909 14 : storage.pm_radiation_indices.push_back(storage.ns_residuals.size());
910 14 : storage.ns_residuals.push_back(std::make_pair(0, 1.0));
911 14 : storage.ns_abs_tols.push_back(_pm_radiation_absolute_tolerance[i]);
912 : }
913 :
914 1183 : storage.converged = storage.ns_residuals.empty();
915 1183 : return storage;
916 0 : }
|
