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 1784 : LinearAssemblySegregatedSolve::validParams()
20 : {
21 1784 : InputParameters params = SIMPLESolveBase::validParams();
22 :
23 3568 : 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 1784 : params.addParam<std::vector<Real>>("active_scalar_equation_relaxation",
30 1784 : 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 3568 : params.addParam<MultiMooseEnum>("active_scalar_petsc_options",
36 3568 : Moose::PetscSupport::getCommonPetscFlags(),
37 : "Singleton PETSc options for the active scalar equation(s)");
38 3568 : params.addParam<MultiMooseEnum>(
39 : "active_scalar_petsc_options_iname",
40 3568 : Moose::PetscSupport::getCommonPetscKeys(),
41 : "Names of PETSc name/value pairs for the active scalar equation(s)");
42 3568 : 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 1784 : params.addParam<std::vector<Real>>(
47 : "active_scalar_absolute_tolerance",
48 1784 : std::vector<Real>(),
49 : "The absolute tolerance(s) on the normalized residual(s) of the active scalar equation(s).");
50 5352 : params.addRangeCheckedParam<Real>("active_scalar_l_tol",
51 3568 : 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 5352 : params.addRangeCheckedParam<Real>("active_scalar_l_abs_tol",
56 3568 : 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 3568 : params.addParam<unsigned int>(
61 : "active_scalar_l_max_its",
62 3568 : 10000,
63 : "The maximum allowed iterations in the linear solver of the turbulence equation.");
64 :
65 3568 : 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 : * Parameters to control the conjugate heat transfer
75 : */
76 1784 : params += NS::FV::CHTHandler::validParams();
77 :
78 1784 : return params;
79 0 : }
80 :
81 892 : LinearAssemblySegregatedSolve::LinearAssemblySegregatedSolve(Executioner & ex)
82 : : SIMPLESolveBase(ex),
83 1784 : _pressure_sys_number(_problem.linearSysNum(getParam<SolverSystemName>("pressure_system"))),
84 892 : _pressure_system(_problem.getLinearSystem(_pressure_sys_number)),
85 1784 : _energy_sys_number(_has_energy_system
86 1560 : ? _problem.linearSysNum(getParam<SolverSystemName>("energy_system"))
87 : : libMesh::invalid_uint),
88 892 : _energy_system(_has_energy_system ? &_problem.getLinearSystem(_energy_sys_number) : nullptr),
89 892 : _solid_energy_sys_number(
90 892 : _has_solid_energy_system
91 1000 : ? _problem.linearSysNum(getParam<SolverSystemName>("solid_energy_system"))
92 : : libMesh::invalid_uint),
93 892 : _solid_energy_system(
94 892 : _has_solid_energy_system ? &_problem.getLinearSystem(_solid_energy_sys_number) : nullptr),
95 1784 : _active_scalar_system_names(getParam<std::vector<SolverSystemName>>("active_scalar_systems")),
96 892 : _has_active_scalar_systems(!_active_scalar_system_names.empty()),
97 1784 : _active_scalar_equation_relaxation(
98 : getParam<std::vector<Real>>("active_scalar_equation_relaxation")),
99 1784 : _active_scalar_l_abs_tol(getParam<Real>("active_scalar_l_abs_tol")),
100 1784 : _active_scalar_absolute_tolerance(
101 : getParam<std::vector<Real>>("active_scalar_absolute_tolerance")),
102 2676 : _cht(ex.parameters())
103 : {
104 : // We fetch the systems and their numbers for the momentum equations.
105 2637 : for (auto system_i : index_range(_momentum_system_names))
106 : {
107 1745 : _momentum_system_numbers.push_back(_problem.linearSysNum(_momentum_system_names[system_i]));
108 1745 : _momentum_systems.push_back(&_problem.getLinearSystem(_momentum_system_numbers[system_i]));
109 1745 : _systems_to_solve.push_back(_momentum_systems.back());
110 : }
111 :
112 892 : _systems_to_solve.push_back(&_pressure_system);
113 :
114 892 : if (_has_energy_system)
115 334 : _systems_to_solve.push_back(_energy_system);
116 :
117 892 : if (_has_solid_energy_system)
118 54 : _systems_to_solve.push_back(_solid_energy_system);
119 : // and for the turbulence surrogate equations
120 892 : if (_has_turbulence_systems)
121 159 : for (auto system_i : index_range(_turbulence_system_names))
122 : {
123 106 : _turbulence_system_numbers.push_back(
124 106 : _problem.linearSysNum(_turbulence_system_names[system_i]));
125 106 : _turbulence_systems.push_back(
126 106 : &_problem.getLinearSystem(_turbulence_system_numbers[system_i]));
127 : }
128 :
129 : // and for the passive scalar equations
130 892 : if (_has_passive_scalar_systems)
131 138 : for (auto system_i : index_range(_passive_scalar_system_names))
132 : {
133 92 : _passive_scalar_system_numbers.push_back(
134 92 : _problem.linearSysNum(_passive_scalar_system_names[system_i]));
135 92 : _passive_scalar_systems.push_back(
136 92 : &_problem.getLinearSystem(_passive_scalar_system_numbers[system_i]));
137 92 : _systems_to_solve.push_back(_passive_scalar_systems.back());
138 : }
139 :
140 : // and for the active scalar equations
141 892 : if (_has_active_scalar_systems)
142 118 : for (auto system_i : index_range(_active_scalar_system_names))
143 : {
144 59 : _active_scalar_system_numbers.push_back(
145 59 : _problem.linearSysNum(_active_scalar_system_names[system_i]));
146 59 : _active_scalar_systems.push_back(
147 59 : &_problem.getLinearSystem(_active_scalar_system_numbers[system_i]));
148 59 : _systems_to_solve.push_back(_active_scalar_systems.back());
149 :
150 : const auto & active_scalar_petsc_options =
151 59 : getParam<MultiMooseEnum>("active_scalar_petsc_options");
152 : const auto & active_scalar_petsc_pair_options = getParam<MooseEnumItem, std::string>(
153 118 : "active_scalar_petsc_options_iname", "active_scalar_petsc_options_value");
154 59 : Moose::PetscSupport::addPetscFlagsToPetscOptions(
155 : active_scalar_petsc_options, "", *this, _active_scalar_petsc_options);
156 118 : Moose::PetscSupport::addPetscPairsToPetscOptions(active_scalar_petsc_pair_options,
157 59 : _problem.mesh().dimension(),
158 : "",
159 : *this,
160 : _active_scalar_petsc_options);
161 :
162 59 : _active_scalar_linear_control.real_valued_data["rel_tol"] =
163 118 : getParam<Real>("active_scalar_l_tol");
164 59 : _active_scalar_linear_control.real_valued_data["abs_tol"] =
165 118 : getParam<Real>("active_scalar_l_abs_tol");
166 59 : _active_scalar_linear_control.int_valued_data["max_its"] =
167 118 : getParam<unsigned int>("active_scalar_l_max_its");
168 59 : }
169 :
170 892 : if (_active_scalar_equation_relaxation.size() != _active_scalar_system_names.size())
171 0 : paramError("active_scalar_equation_relaxation",
172 : "Should be the same size as the number of systems");
173 :
174 : // We disable the prefix here for the time being, the segregated solvers use a different approach
175 : // for setting the petsc parameters
176 4068 : for (auto & system : _systems_to_solve)
177 3176 : system->system().prefix_with_name(false);
178 :
179 : // Link CHT objects, this will also do some error checking
180 892 : if (_cht.enabled())
181 36 : _cht.linkEnergySystems(_solid_energy_system, _energy_system);
182 892 : }
183 :
184 : void
185 890 : LinearAssemblySegregatedSolve::linkRhieChowUserObject()
186 : {
187 890 : _rc_uo =
188 890 : const_cast<RhieChowMassFlux *>(&getUserObject<RhieChowMassFlux>("rhie_chow_user_object"));
189 890 : _rc_uo->linkMomentumPressureSystems(
190 890 : _momentum_systems, _pressure_system, _momentum_system_numbers);
191 :
192 : // Initialize the face velocities in the RC object
193 890 : if (!_app.isRecovering())
194 812 : _rc_uo->initFaceMassFlux();
195 890 : _rc_uo->initCouplingField();
196 890 : }
197 :
198 : std::vector<std::pair<unsigned int, Real>>
199 171133 : LinearAssemblySegregatedSolve::solveMomentumPredictor()
200 : {
201 : // Temporary storage for the (flux-normalized) residuals from
202 : // different momentum components
203 : std::vector<std::pair<unsigned int, Real>> its_normalized_residuals;
204 :
205 : LinearImplicitSystem & momentum_system_0 =
206 171133 : libMesh::cast_ref<LinearImplicitSystem &>(_momentum_systems[0]->system());
207 :
208 : PetscLinearSolver<Real> & momentum_solver =
209 171133 : libMesh::cast_ref<PetscLinearSolver<Real> &>(*momentum_system_0.get_linear_solver());
210 :
211 : // Solve the momentum equations.
212 : // TO DO: These equations are VERY similar. If we can store the differences (things coming from
213 : // BCs for example) separately, it is enough to construct one matrix.
214 450586 : for (const auto system_i : index_range(_momentum_systems))
215 : {
216 279453 : _problem.setCurrentLinearSystem(_momentum_system_numbers[system_i]);
217 :
218 : // We will need the right hand side and the solution of the next component
219 : LinearImplicitSystem & momentum_system =
220 279453 : libMesh::cast_ref<LinearImplicitSystem &>(_momentum_systems[system_i]->system());
221 :
222 : NumericVector<Number> & solution = *(momentum_system.solution);
223 279453 : NumericVector<Number> & rhs = *(momentum_system.rhs);
224 279453 : SparseMatrix<Number> & mmat = *(momentum_system.matrix);
225 :
226 279453 : auto diff_diagonal = solution.zero_clone();
227 :
228 : // We assemble the matrix and the right hand side
229 279453 : _problem.computeLinearSystemSys(momentum_system, mmat, rhs, /*compute_grads*/ true);
230 :
231 : // Still need to relax the right hand side with the same vector
232 279453 : NS::FV::relaxMatrix(mmat, _momentum_equation_relaxation, *diff_diagonal);
233 279453 : NS::FV::relaxRightHandSide(rhs, solution, *diff_diagonal);
234 :
235 : // The normalization factor depends on the right hand side so we need to recompute it for this
236 : // component
237 279453 : Real norm_factor = NS::FV::computeNormalizationFactor(solution, mmat, rhs);
238 :
239 : // Very important, for deciding the convergence, we need the unpreconditioned
240 : // norms in the linear solve
241 279453 : LibmeshPetscCall(KSPSetNormType(momentum_solver.ksp(), KSP_NORM_UNPRECONDITIONED));
242 : // Solve this component. We don't update the ghosted solution yet, that will come at the end
243 : // of the corrector step. Also setting the linear tolerances and maximum iteration counts.
244 279453 : _momentum_linear_control.real_valued_data["abs_tol"] = _momentum_l_abs_tol * norm_factor;
245 279453 : momentum_solver.set_solver_configuration(_momentum_linear_control);
246 :
247 : // We solve the equation
248 279453 : auto its_resid_pair = momentum_solver.solve(mmat, mmat, solution, rhs);
249 279453 : momentum_system.update();
250 :
251 : // We will reuse the preconditioner for every momentum system
252 279453 : if (system_i == 0)
253 171133 : momentum_solver.reuse_preconditioner(true);
254 :
255 : // Save the normalized residual
256 : its_normalized_residuals.push_back(
257 279453 : std::make_pair(its_resid_pair.first, momentum_solver.get_initial_residual() / norm_factor));
258 :
259 279453 : if (_print_fields)
260 : {
261 0 : _console << " solution after solve " << std::endl;
262 0 : solution.print();
263 0 : _console << " matrix when we solve " << std::endl;
264 0 : mmat.print();
265 0 : _console << " rhs when we solve " << std::endl;
266 0 : rhs.print();
267 0 : _console << " velocity solution component " << system_i << std::endl;
268 0 : solution.print();
269 0 : _console << "Norm factor " << norm_factor << std::endl;
270 0 : _console << Moose::stringify(momentum_solver.get_initial_residual()) << std::endl;
271 : }
272 :
273 : // Printing residuals
274 279453 : _console << " Momentum equation:"
275 : << (_momentum_systems.size() > 1
276 1205451 : ? std::string(" Component ") + std::to_string(system_i + 1) + std::string(" ")
277 279453 : : std::string(" "))
278 394869 : << COLOR_GREEN << its_normalized_residuals[system_i].second << COLOR_DEFAULT
279 279453 : << " Linear its: " << its_normalized_residuals[system_i].first << std::endl;
280 279453 : }
281 :
282 450586 : for (const auto system_i : index_range(_momentum_systems))
283 : {
284 : LinearImplicitSystem & momentum_system =
285 279453 : libMesh::cast_ref<LinearImplicitSystem &>(_momentum_systems[system_i]->system());
286 279453 : _momentum_systems[system_i]->setSolution(*(momentum_system.current_local_solution));
287 279453 : _momentum_systems[system_i]->copyPreviousNonlinearSolutions();
288 : }
289 :
290 : // We reset this to ensure the preconditioner is recomputed new time we go to the momentum
291 : // predictor
292 171133 : momentum_solver.reuse_preconditioner(false);
293 :
294 171133 : return its_normalized_residuals;
295 0 : }
296 :
297 : void
298 749 : LinearAssemblySegregatedSolve::initialSetup()
299 : {
300 749 : if (_cht.enabled())
301 : {
302 36 : _cht.deduceCHTBoundaryCoupling();
303 36 : _cht.setupConjugateHeatTransferContainers();
304 : }
305 749 : }
306 :
307 : std::pair<unsigned int, Real>
308 192833 : LinearAssemblySegregatedSolve::solvePressureCorrector()
309 : {
310 192833 : _problem.setCurrentLinearSystem(_pressure_sys_number);
311 :
312 : // We will need some members from the linear system
313 : LinearImplicitSystem & pressure_system =
314 192833 : libMesh::cast_ref<LinearImplicitSystem &>(_pressure_system.system());
315 :
316 : // We will need the solution, the right hand side and the matrix
317 : NumericVector<Number> & current_local_solution = *(pressure_system.current_local_solution);
318 : NumericVector<Number> & solution = *(pressure_system.solution);
319 192833 : SparseMatrix<Number> & mmat = *(pressure_system.matrix);
320 192833 : NumericVector<Number> & rhs = *(pressure_system.rhs);
321 :
322 : // Fetch the linear solver from the system
323 : PetscLinearSolver<Real> & pressure_solver =
324 192833 : libMesh::cast_ref<PetscLinearSolver<Real> &>(*pressure_system.get_linear_solver());
325 :
326 192833 : _problem.computeLinearSystemSys(pressure_system, mmat, rhs, false);
327 :
328 192833 : if (_print_fields)
329 : {
330 0 : _console << "Pressure matrix" << std::endl;
331 0 : mmat.print();
332 : }
333 :
334 : // We compute the normalization factors based on the fluxes
335 192833 : Real norm_factor = NS::FV::computeNormalizationFactor(solution, mmat, rhs);
336 :
337 : // We need the non-preconditioned norm to be consistent with the norm factor
338 192833 : LibmeshPetscCall(KSPSetNormType(pressure_solver.ksp(), KSP_NORM_UNPRECONDITIONED));
339 :
340 : // Setting the linear tolerances and maximum iteration counts
341 192833 : _pressure_linear_control.real_valued_data["abs_tol"] = _pressure_l_abs_tol * norm_factor;
342 192833 : pressure_solver.set_solver_configuration(_pressure_linear_control);
343 :
344 192833 : if (_pin_pressure)
345 52925 : NS::FV::constrainSystem(mmat, rhs, _pressure_pin_value, _pressure_pin_dof);
346 192833 : pressure_system.update();
347 :
348 192833 : auto its_res_pair = pressure_solver.solve(mmat, mmat, solution, rhs);
349 192833 : pressure_system.update();
350 :
351 192833 : if (_print_fields)
352 : {
353 0 : _console << " rhs when we solve pressure " << std::endl;
354 0 : rhs.print();
355 0 : _console << " Pressure " << std::endl;
356 0 : solution.print();
357 0 : _console << "Norm factor " << norm_factor << std::endl;
358 : }
359 :
360 192833 : _pressure_system.setSolution(current_local_solution);
361 :
362 : const auto residuals =
363 192833 : std::make_pair(its_res_pair.first, pressure_solver.get_initial_residual() / norm_factor);
364 :
365 276467 : _console << " Pressure equation: " << COLOR_GREEN << residuals.second << COLOR_DEFAULT
366 192833 : << " Linear its: " << residuals.first << std::endl;
367 :
368 192833 : return residuals;
369 : }
370 :
371 : std::pair<unsigned int, Real>
372 12564 : LinearAssemblySegregatedSolve::solveSolidEnergy()
373 : {
374 12564 : _problem.setCurrentLinearSystem(_solid_energy_sys_number);
375 :
376 : // We will need some members from the linear system
377 : LinearImplicitSystem & system =
378 12564 : libMesh::cast_ref<LinearImplicitSystem &>(_solid_energy_system->system());
379 :
380 : // We will need the solution, the right hand side and the matrix
381 : NumericVector<Number> & current_local_solution = *(system.current_local_solution);
382 : NumericVector<Number> & solution = *(system.solution);
383 12564 : SparseMatrix<Number> & mmat = *(system.matrix);
384 12564 : NumericVector<Number> & rhs = *(system.rhs);
385 :
386 : // Fetch the linear solver from the system
387 : PetscLinearSolver<Real> & solver =
388 12564 : libMesh::cast_ref<PetscLinearSolver<Real> &>(*system.get_linear_solver());
389 :
390 12564 : _problem.computeLinearSystemSys(system, mmat, rhs, false);
391 :
392 12564 : if (_print_fields)
393 : {
394 0 : _console << "Solid energy matrix" << std::endl;
395 0 : mmat.print();
396 : }
397 :
398 : // We compute the normalization factors based on the fluxes
399 12564 : Real norm_factor = NS::FV::computeNormalizationFactor(solution, mmat, rhs);
400 :
401 : // We need the non-preconditioned norm to be consistent with the norm factor
402 12564 : LibmeshPetscCall(KSPSetNormType(solver.ksp(), KSP_NORM_UNPRECONDITIONED));
403 :
404 : // Setting the linear tolerances and maximum iteration counts
405 12564 : _solid_energy_linear_control.real_valued_data["abs_tol"] = _solid_energy_l_abs_tol * norm_factor;
406 12564 : solver.set_solver_configuration(_solid_energy_linear_control);
407 :
408 12564 : auto its_res_pair = solver.solve(mmat, mmat, solution, rhs);
409 12564 : system.update();
410 :
411 12564 : if (_print_fields)
412 : {
413 0 : _console << " rhs when we solve solid energy " << std::endl;
414 0 : rhs.print();
415 0 : _console << " Solid energy " << std::endl;
416 0 : solution.print();
417 0 : _console << "Norm factor " << norm_factor << std::endl;
418 : }
419 :
420 12564 : _solid_energy_system->setSolution(current_local_solution);
421 :
422 : const auto residuals =
423 12564 : std::make_pair(its_res_pair.first, solver.get_initial_residual() / norm_factor);
424 :
425 20940 : _console << " Solid energy equation: " << COLOR_GREEN << residuals.second << COLOR_DEFAULT
426 12564 : << " Linear its: " << residuals.first << std::endl;
427 :
428 12564 : return residuals;
429 : }
430 :
431 : std::pair<unsigned int, Real>
432 192833 : LinearAssemblySegregatedSolve::correctVelocity(const bool subtract_updated_pressure,
433 : const bool recompute_face_mass_flux,
434 : const SolverParams & solver_params)
435 : {
436 : // Compute the coupling fields between the momentum and pressure equations.
437 : // The first argument makes sure the pressure gradient is staged at the first
438 : // iteration
439 192833 : _rc_uo->computeHbyA(subtract_updated_pressure, _print_fields);
440 :
441 : // We set the preconditioner/controllable parameters for the pressure equations through
442 : // petsc options. Linear tolerances will be overridden within the solver.
443 192833 : Moose::PetscSupport::petscSetOptions(_pressure_petsc_options, solver_params);
444 :
445 : // Solve the pressure corrector
446 192833 : const auto residuals = solvePressureCorrector();
447 :
448 : // Compute the face velocity which is used in the advection terms. In certain
449 : // segregated solver algorithms (like PISO) this is only done on the last iteration.
450 192833 : if (recompute_face_mass_flux)
451 171133 : _rc_uo->computeFaceMassFlux();
452 :
453 192833 : auto & pressure_current_solution = *(_pressure_system.system().current_local_solution.get());
454 192833 : auto & pressure_old_solution = *(_pressure_system.solutionPreviousNewton());
455 :
456 : // Relax the pressure update for the next momentum predictor
457 192833 : NS::FV::relaxSolutionUpdate(
458 192833 : pressure_current_solution, pressure_old_solution, _pressure_variable_relaxation);
459 :
460 : // Overwrite old solution
461 192833 : pressure_old_solution = pressure_current_solution;
462 192833 : _pressure_system.setSolution(pressure_current_solution);
463 :
464 : // We recompute the updated pressure gradient
465 192833 : _pressure_system.computeGradients();
466 :
467 : // Reconstruct the cell velocity as well to accelerate convergence
468 192833 : _rc_uo->computeCellVelocity();
469 :
470 192833 : return residuals;
471 : }
472 :
473 : std::pair<unsigned int, Real>
474 120422 : LinearAssemblySegregatedSolve::solveAdvectedSystem(const unsigned int system_num,
475 : LinearSystem & system,
476 : const Real relaxation_factor,
477 : SolverConfiguration & solver_config,
478 : const Real absolute_tol,
479 : const Real field_relaxation,
480 : const Real min_value_limiter)
481 : {
482 120422 : _problem.setCurrentLinearSystem(system_num);
483 :
484 : // We will need some members from the implicit linear system
485 120422 : LinearImplicitSystem & li_system = libMesh::cast_ref<LinearImplicitSystem &>(system.system());
486 :
487 : // We will need the solution, the right hand side and the matrix
488 : NumericVector<Number> & current_local_solution = *(li_system.current_local_solution);
489 : NumericVector<Number> & solution = *(li_system.solution);
490 120422 : SparseMatrix<Number> & mmat = *(li_system.matrix);
491 120422 : NumericVector<Number> & rhs = *(li_system.rhs);
492 :
493 : // We need a vector that stores the (diagonal_relaxed-original_diagonal) vector
494 120422 : auto diff_diagonal = solution.zero_clone();
495 :
496 : // Fetch the linear solver from the system
497 : PetscLinearSolver<Real> & linear_solver =
498 120422 : libMesh::cast_ref<PetscLinearSolver<Real> &>(*li_system.get_linear_solver());
499 :
500 120422 : _problem.computeLinearSystemSys(li_system, mmat, rhs, true);
501 :
502 : // Go and relax the system matrix and the right hand side
503 120422 : NS::FV::relaxMatrix(mmat, relaxation_factor, *diff_diagonal);
504 120422 : NS::FV::relaxRightHandSide(rhs, solution, *diff_diagonal);
505 :
506 120422 : if (_print_fields)
507 : {
508 0 : _console << system.name() << " system matrix" << std::endl;
509 0 : mmat.print();
510 : }
511 :
512 : // We compute the normalization factors based on the fluxes
513 120422 : Real norm_factor = NS::FV::computeNormalizationFactor(solution, mmat, rhs);
514 :
515 : // We need the non-preconditioned norm to be consistent with the norm factor
516 120422 : LibmeshPetscCall(KSPSetNormType(linear_solver.ksp(), KSP_NORM_UNPRECONDITIONED));
517 :
518 : // Setting the linear tolerances and maximum iteration counts
519 120422 : solver_config.real_valued_data["abs_tol"] = absolute_tol * norm_factor;
520 120422 : linear_solver.set_solver_configuration(solver_config);
521 :
522 : // Solve the system and update current local solution
523 120422 : auto its_res_pair = linear_solver.solve(mmat, mmat, solution, rhs);
524 120422 : li_system.update();
525 :
526 120422 : if (_print_fields)
527 : {
528 0 : _console << " rhs when we solve " << system.name() << std::endl;
529 0 : rhs.print();
530 0 : _console << system.name() << " solution " << std::endl;
531 0 : solution.print();
532 0 : _console << " Norm factor " << norm_factor << std::endl;
533 : }
534 :
535 : // Limiting scalar solution
536 120422 : if (min_value_limiter != std::numeric_limits<Real>::min())
537 43306 : NS::FV::limitSolutionUpdate(current_local_solution, min_value_limiter);
538 :
539 : // Relax the field update for the next momentum predictor
540 120422 : if (field_relaxation != 1.0)
541 : {
542 18090 : auto & old_local_solution = *(system.solutionPreviousNewton());
543 18090 : NS::FV::relaxSolutionUpdate(current_local_solution, old_local_solution, field_relaxation);
544 :
545 : // Update old solution, only needed if relaxing the field
546 18090 : old_local_solution = current_local_solution;
547 : }
548 :
549 120422 : system.setSolution(current_local_solution);
550 :
551 : const auto residuals =
552 120422 : std::make_pair(its_res_pair.first, linear_solver.get_initial_residual() / norm_factor);
553 :
554 153802 : _console << " Advected system: " << system.name() << " " << COLOR_GREEN << residuals.second
555 153802 : << COLOR_DEFAULT << " Linear its: " << residuals.first << std::endl;
556 :
557 120422 : return residuals;
558 120422 : }
559 :
560 : bool
561 1707 : LinearAssemblySegregatedSolve::solve()
562 : {
563 : // Do not solve if problem is set not to
564 1707 : if (!_problem.shouldSolve())
565 : return true;
566 :
567 : // Dummy solver parameter file which is needed for switching petsc options
568 1707 : SolverParams solver_params;
569 1707 : solver_params._type = Moose::SolveType::ST_LINEAR;
570 1707 : solver_params._line_search = Moose::LineSearchType::LS_NONE;
571 :
572 : // Initialize the SIMPLE iteration counter
573 1707 : unsigned int simple_iteration_counter = 0;
574 :
575 : // Assign residuals to general residual vector
576 1707 : const unsigned int no_systems = _momentum_systems.size() + 1 + _has_energy_system +
577 1707 : _has_solid_energy_system + _active_scalar_systems.size() +
578 1707 : _turbulence_systems.size();
579 :
580 1707 : std::vector<std::pair<unsigned int, Real>> ns_residuals(no_systems, std::make_pair(0, 1.0));
581 1707 : std::vector<Real> ns_abs_tols(_momentum_systems.size(), _momentum_absolute_tolerance);
582 1707 : ns_abs_tols.push_back(_pressure_absolute_tolerance);
583 :
584 : // Push back energy tolerances
585 1707 : if (_has_energy_system)
586 904 : ns_abs_tols.push_back(_energy_absolute_tolerance);
587 1707 : if (_has_solid_energy_system)
588 54 : ns_abs_tols.push_back(_solid_energy_absolute_tolerance);
589 1707 : if (_has_active_scalar_systems)
590 338 : for (const auto scalar_tol : _active_scalar_absolute_tolerance)
591 169 : ns_abs_tols.push_back(scalar_tol);
592 :
593 : // Push back turbulence tolerances
594 1707 : if (_has_turbulence_systems)
595 405 : for (const auto turbulence_tol : _turbulence_absolute_tolerance)
596 270 : ns_abs_tols.push_back(turbulence_tol);
597 :
598 : bool converged = false;
599 : // Loop until converged or hit the maximum allowed iteration number
600 1707 : if (_cht.enabled())
601 36 : _cht.initializeCHTCouplingFields();
602 :
603 172840 : while (simple_iteration_counter < _num_iterations && !converged)
604 : {
605 171133 : simple_iteration_counter++;
606 :
607 : // We set the preconditioner/controllable parameters through petsc options. Linear
608 : // tolerances will be overridden within the solver. In case of a segregated momentum
609 : // solver, we assume that every velocity component uses the same preconditioner
610 171133 : Moose::PetscSupport::petscSetOptions(_momentum_petsc_options, solver_params);
611 :
612 : // Initialize pressure gradients, after this we just reuse the last ones from each
613 : // iteration
614 171133 : if (simple_iteration_counter == 1)
615 1707 : _pressure_system.computeGradients();
616 :
617 171133 : _console << "Iteration " << simple_iteration_counter << " Initial residual norms:" << std::endl;
618 :
619 : // Solve the momentum predictor step
620 171133 : auto momentum_residual = solveMomentumPredictor();
621 450586 : for (const auto system_i : index_range(momentum_residual))
622 : ns_residuals[system_i] = momentum_residual[system_i];
623 :
624 : // Now we correct the velocity, this function depends on the method, it differs for
625 : // SIMPLE/PIMPLE, this returns the pressure errors
626 171133 : ns_residuals[momentum_residual.size()] = correctVelocity(true, true, solver_params);
627 :
628 : // If we have an energy equation, solve it here.We assume the material properties in the
629 : // Navier-Stokes equations depend on temperature, therefore we can not solve for temperature
630 : // outside of the velocity-pressure loop
631 171133 : if (_has_energy_system)
632 : {
633 : // If there is no CHT specified this will just do go once through this block
634 : _cht.resetCHTConvergence();
635 67764 : while (!_cht.converged())
636 : {
637 35529 : if (_cht.enabled())
638 7506 : _cht.updateCHTBoundaryCouplingFields(NS::CHTSide::FLUID);
639 :
640 : // We set the preconditioner/controllable parameters through petsc options. Linear
641 : // tolerances will be overridden within the solver.
642 35529 : Moose::PetscSupport::petscSetOptions(_energy_petsc_options, solver_params);
643 35529 : ns_residuals[momentum_residual.size() + _has_energy_system] =
644 71058 : solveAdvectedSystem(_energy_sys_number,
645 35529 : *_energy_system,
646 35529 : _energy_equation_relaxation,
647 : _energy_linear_control,
648 35529 : _energy_l_abs_tol);
649 :
650 35529 : if (_has_solid_energy_system)
651 : {
652 : // For now we only update gradients if cht is needed, might change in the future
653 12564 : if (_cht.enabled())
654 : {
655 7506 : _energy_system->computeGradients();
656 7506 : _cht.updateCHTBoundaryCouplingFields(NS::CHTSide::SOLID);
657 : }
658 :
659 : // We set the preconditioner/controllable parameters through petsc options. Linear
660 : // tolerances will be overridden within the solver.
661 12564 : Moose::PetscSupport::petscSetOptions(_solid_energy_petsc_options, solver_params);
662 12564 : ns_residuals[momentum_residual.size() + _has_solid_energy_system + _has_energy_system] =
663 12564 : solveSolidEnergy();
664 :
665 : // For now we only update gradients if cht is needed, might change in the future
666 12564 : if (_cht.enabled())
667 7506 : _solid_energy_system->computeGradients();
668 : }
669 :
670 35529 : if (_cht.enabled())
671 : {
672 7506 : _cht.sumIntegratedFluxes();
673 7506 : _cht.printIntegratedFluxes();
674 : }
675 :
676 : _cht.incrementCHTIterators();
677 : }
678 32235 : if (_cht.enabled())
679 4212 : _cht.resetIntegratedFluxes();
680 : }
681 :
682 : // If we have active scalar equations, solve them here in case they depend on temperature
683 : // or they affect the fluid properties such that they must be solved concurrently with
684 : // pressure and velocity
685 171133 : if (_has_active_scalar_systems)
686 : {
687 24023 : _problem.execute(EXEC_NONLINEAR);
688 :
689 : // We set the preconditioner/controllable parameters through petsc options. Linear
690 : // tolerances will be overridden within the solver.
691 24023 : Moose::PetscSupport::petscSetOptions(_active_scalar_petsc_options, solver_params);
692 48046 : for (const auto i : index_range(_active_scalar_system_names))
693 24023 : ns_residuals[momentum_residual.size() + 1 + _has_energy_system + _has_solid_energy_system +
694 24023 : i] = solveAdvectedSystem(_active_scalar_system_numbers[i],
695 : *_active_scalar_systems[i],
696 : _active_scalar_equation_relaxation[i],
697 : _active_scalar_linear_control,
698 24023 : _active_scalar_l_abs_tol);
699 : }
700 :
701 : // If we have turbulence equations, solve them here.
702 : // The turbulent viscosity depends on the value of the turbulence surrogate variables
703 171133 : if (_has_turbulence_systems)
704 : {
705 : // We set the preconditioner/controllable parameters through petsc options. Linear
706 : // tolerances will be overridden within the solver.
707 21653 : Moose::PetscSupport::petscSetOptions(_turbulence_petsc_options, solver_params);
708 64959 : for (const auto i : index_range(_turbulence_system_names))
709 : {
710 43306 : ns_residuals[momentum_residual.size() + 1 + _has_energy_system + _has_solid_energy_system +
711 43306 : _active_scalar_system_names.size() + i] =
712 43306 : solveAdvectedSystem(_turbulence_system_numbers[i],
713 : *_turbulence_systems[i],
714 : _turbulence_equation_relaxation[i],
715 : _turbulence_linear_control,
716 43306 : _turbulence_l_abs_tol,
717 : _turbulence_field_relaxation[i],
718 : _turbulence_field_min_limit[i]);
719 : }
720 : }
721 :
722 171133 : _problem.execute(EXEC_NONLINEAR);
723 :
724 171133 : converged = NS::FV::converged(ns_residuals, ns_abs_tols);
725 171133 : }
726 :
727 : // If we have passive scalar equations, solve them here. We assume the material properties in
728 : // the Navier-Stokes equations do not depend on passive scalars, as they are passive, therefore
729 : // we solve outside of the velocity-pressure loop
730 1707 : if (_has_passive_scalar_systems && (converged || _continue_on_max_its))
731 : {
732 : // The reason why we need more than one iteration is due to the matrix relaxation
733 : // which can be used to stabilize the equations
734 : bool passive_scalar_converged = false;
735 98 : unsigned int ps_iteration_counter = 0;
736 :
737 98 : _console << "Passive scalar iteration " << ps_iteration_counter
738 98 : << " Initial residual norms:" << std::endl;
739 :
740 8880 : while (ps_iteration_counter < _num_iterations && !passive_scalar_converged)
741 : {
742 8782 : ps_iteration_counter++;
743 : std::vector<std::pair<unsigned int, Real>> scalar_residuals(
744 8782 : _passive_scalar_system_names.size(), std::make_pair(0, 1.0));
745 : std::vector<Real> scalar_abs_tols;
746 26346 : for (const auto scalar_tol : _passive_scalar_absolute_tolerance)
747 17564 : scalar_abs_tols.push_back(scalar_tol);
748 :
749 : // We set the preconditioner/controllable parameters through petsc options. Linear
750 : // tolerances will be overridden within the solver.
751 8782 : Moose::PetscSupport::petscSetOptions(_passive_scalar_petsc_options, solver_params);
752 26346 : for (const auto i : index_range(_passive_scalar_system_names))
753 17564 : scalar_residuals[i] = solveAdvectedSystem(_passive_scalar_system_numbers[i],
754 : *_passive_scalar_systems[i],
755 : _passive_scalar_equation_relaxation[i],
756 : _passive_scalar_linear_control,
757 17564 : _passive_scalar_l_abs_tol);
758 :
759 8782 : passive_scalar_converged = NS::FV::converged(scalar_residuals, scalar_abs_tols);
760 8782 : }
761 :
762 : // Both flow and scalars must converge
763 98 : converged = passive_scalar_converged && converged;
764 : }
765 :
766 1707 : converged = _continue_on_max_its ? true : converged;
767 :
768 : return converged;
769 1707 : }
|
