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 :
15 : using namespace libMesh;
16 :
17 : InputParameters
18 1618 : LinearAssemblySegregatedSolve::validParams()
19 : {
20 1618 : InputParameters params = SIMPLESolveBase::validParams();
21 :
22 3236 : 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 1618 : params.addParam<std::vector<Real>>("active_scalar_equation_relaxation",
29 1618 : 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 3236 : params.addParam<MultiMooseEnum>("active_scalar_petsc_options",
35 3236 : Moose::PetscSupport::getCommonPetscFlags(),
36 : "Singleton PETSc options for the active scalar equation(s)");
37 3236 : params.addParam<MultiMooseEnum>(
38 : "active_scalar_petsc_options_iname",
39 3236 : Moose::PetscSupport::getCommonPetscKeys(),
40 : "Names of PETSc name/value pairs for the active scalar equation(s)");
41 3236 : 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 1618 : params.addParam<std::vector<Real>>(
46 : "active_scalar_absolute_tolerance",
47 1618 : std::vector<Real>(),
48 : "The absolute tolerance(s) on the normalized residual(s) of the active scalar equation(s).");
49 4854 : params.addRangeCheckedParam<Real>("active_scalar_l_tol",
50 3236 : 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 4854 : params.addRangeCheckedParam<Real>("active_scalar_l_abs_tol",
55 3236 : 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 3236 : params.addParam<unsigned int>(
60 : "active_scalar_l_max_its",
61 3236 : 10000,
62 : "The maximum allowed iterations in the linear solver of the turbulence equation.");
63 :
64 3236 : 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 1618 : return params;
73 0 : }
74 :
75 809 : LinearAssemblySegregatedSolve::LinearAssemblySegregatedSolve(Executioner & ex)
76 : : SIMPLESolveBase(ex),
77 1618 : _pressure_sys_number(_problem.linearSysNum(getParam<SolverSystemName>("pressure_system"))),
78 809 : _pressure_system(_problem.getLinearSystem(_pressure_sys_number)),
79 1618 : _energy_sys_number(_has_energy_system
80 1329 : ? _problem.linearSysNum(getParam<SolverSystemName>("energy_system"))
81 : : libMesh::invalid_uint),
82 809 : _energy_system(_has_energy_system ? &_problem.getLinearSystem(_energy_sys_number) : nullptr),
83 809 : _solid_energy_sys_number(
84 809 : _has_solid_energy_system
85 845 : ? _problem.linearSysNum(getParam<SolverSystemName>("solid_energy_system"))
86 : : libMesh::invalid_uint),
87 809 : _solid_energy_system(
88 809 : _has_solid_energy_system ? &_problem.getLinearSystem(_solid_energy_sys_number) : nullptr),
89 1618 : _active_scalar_system_names(getParam<std::vector<SolverSystemName>>("active_scalar_systems")),
90 809 : _has_active_scalar_systems(!_active_scalar_system_names.empty()),
91 1618 : _active_scalar_equation_relaxation(
92 : getParam<std::vector<Real>>("active_scalar_equation_relaxation")),
93 1618 : _active_scalar_l_abs_tol(getParam<Real>("active_scalar_l_abs_tol")),
94 1618 : _active_scalar_absolute_tolerance(
95 1618 : getParam<std::vector<Real>>("active_scalar_absolute_tolerance"))
96 : {
97 : // We fetch the systems and their numbers for the momentum equations.
98 2388 : for (auto system_i : index_range(_momentum_system_names))
99 : {
100 1579 : _momentum_system_numbers.push_back(_problem.linearSysNum(_momentum_system_names[system_i]));
101 1579 : _momentum_systems.push_back(&_problem.getLinearSystem(_momentum_system_numbers[system_i]));
102 1579 : _systems_to_solve.push_back(_momentum_systems.back());
103 : }
104 :
105 809 : _systems_to_solve.push_back(&_pressure_system);
106 :
107 809 : if (_has_energy_system)
108 260 : _systems_to_solve.push_back(_energy_system);
109 :
110 809 : if (_has_solid_energy_system)
111 18 : _systems_to_solve.push_back(_solid_energy_system);
112 : // and for the turbulence surrogate equations
113 809 : if (_has_turbulence_systems)
114 132 : for (auto system_i : index_range(_turbulence_system_names))
115 : {
116 88 : _turbulence_system_numbers.push_back(
117 88 : _problem.linearSysNum(_turbulence_system_names[system_i]));
118 88 : _turbulence_systems.push_back(
119 88 : &_problem.getLinearSystem(_turbulence_system_numbers[system_i]));
120 : }
121 :
122 : // and for the passive scalar equations
123 809 : if (_has_passive_scalar_systems)
124 138 : for (auto system_i : index_range(_passive_scalar_system_names))
125 : {
126 92 : _passive_scalar_system_numbers.push_back(
127 92 : _problem.linearSysNum(_passive_scalar_system_names[system_i]));
128 92 : _passive_scalar_systems.push_back(
129 92 : &_problem.getLinearSystem(_passive_scalar_system_numbers[system_i]));
130 92 : _systems_to_solve.push_back(_passive_scalar_systems.back());
131 : }
132 :
133 : // and for the active scalar equations
134 809 : if (_has_active_scalar_systems)
135 118 : for (auto system_i : index_range(_active_scalar_system_names))
136 : {
137 59 : _active_scalar_system_numbers.push_back(
138 59 : _problem.linearSysNum(_active_scalar_system_names[system_i]));
139 59 : _active_scalar_systems.push_back(
140 59 : &_problem.getLinearSystem(_active_scalar_system_numbers[system_i]));
141 59 : _systems_to_solve.push_back(_active_scalar_systems.back());
142 :
143 : const auto & active_scalar_petsc_options =
144 59 : getParam<MultiMooseEnum>("active_scalar_petsc_options");
145 : const auto & active_scalar_petsc_pair_options = getParam<MooseEnumItem, std::string>(
146 118 : "active_scalar_petsc_options_iname", "active_scalar_petsc_options_value");
147 59 : Moose::PetscSupport::addPetscFlagsToPetscOptions(
148 : active_scalar_petsc_options, "-", *this, _active_scalar_petsc_options);
149 118 : Moose::PetscSupport::addPetscPairsToPetscOptions(active_scalar_petsc_pair_options,
150 59 : _problem.mesh().dimension(),
151 : "-",
152 : *this,
153 : _active_scalar_petsc_options);
154 :
155 59 : _active_scalar_linear_control.real_valued_data["rel_tol"] =
156 118 : getParam<Real>("active_scalar_l_tol");
157 59 : _active_scalar_linear_control.real_valued_data["abs_tol"] =
158 118 : getParam<Real>("active_scalar_l_abs_tol");
159 59 : _active_scalar_linear_control.int_valued_data["max_its"] =
160 118 : getParam<unsigned int>("active_scalar_l_max_its");
161 59 : }
162 :
163 809 : if (_active_scalar_equation_relaxation.size() != _active_scalar_system_names.size())
164 0 : 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 3626 : for (auto & system : _systems_to_solve)
170 2817 : system->system().prefix_with_name(false);
171 809 : }
172 :
173 : void
174 807 : LinearAssemblySegregatedSolve::linkRhieChowUserObject()
175 : {
176 807 : _rc_uo =
177 807 : const_cast<RhieChowMassFlux *>(&getUserObject<RhieChowMassFlux>("rhie_chow_user_object"));
178 807 : _rc_uo->linkMomentumPressureSystems(
179 807 : _momentum_systems, _pressure_system, _momentum_system_numbers);
180 :
181 : // Initialize the face velocities in the RC object
182 807 : if (!_app.isRecovering())
183 739 : _rc_uo->initFaceMassFlux();
184 807 : _rc_uo->initCouplingField();
185 807 : }
186 :
187 : std::vector<std::pair<unsigned int, Real>>
188 162117 : LinearAssemblySegregatedSolve::solveMomentumPredictor()
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 162117 : libMesh::cast_ref<LinearImplicitSystem &>(_momentum_systems[0]->system());
196 :
197 : PetscLinearSolver<Real> & momentum_solver =
198 162117 : 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 423538 : for (const auto system_i : index_range(_momentum_systems))
204 : {
205 261421 : _problem.setCurrentLinearSystem(_momentum_system_numbers[system_i]);
206 :
207 : // We will need the right hand side and the solution of the next component
208 : LinearImplicitSystem & momentum_system =
209 261421 : libMesh::cast_ref<LinearImplicitSystem &>(_momentum_systems[system_i]->system());
210 :
211 : NumericVector<Number> & solution = *(momentum_system.solution);
212 261421 : NumericVector<Number> & rhs = *(momentum_system.rhs);
213 261421 : SparseMatrix<Number> & mmat = *(momentum_system.matrix);
214 :
215 261421 : auto diff_diagonal = solution.zero_clone();
216 :
217 : // We assemble the matrix and the right hand side
218 261421 : _problem.computeLinearSystemSys(momentum_system, mmat, rhs, /*compute_grads*/ true);
219 :
220 : // Still need to relax the right hand side with the same vector
221 261421 : NS::FV::relaxMatrix(mmat, _momentum_equation_relaxation, *diff_diagonal);
222 261421 : 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 261421 : 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 261421 : 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.
233 261421 : _momentum_linear_control.real_valued_data["abs_tol"] = _momentum_l_abs_tol * norm_factor;
234 261421 : momentum_solver.set_solver_configuration(_momentum_linear_control);
235 :
236 : // We solve the equation
237 261421 : auto its_resid_pair = momentum_solver.solve(mmat, mmat, solution, rhs);
238 261421 : momentum_system.update();
239 :
240 : // We will reuse the preconditioner for every momentum system
241 261421 : if (system_i == 0)
242 162117 : momentum_solver.reuse_preconditioner(true);
243 :
244 : // Save the normalized residual
245 : its_normalized_residuals.push_back(
246 261421 : std::make_pair(its_resid_pair.first, momentum_solver.get_initial_residual() / norm_factor));
247 :
248 261421 : if (_print_fields)
249 : {
250 0 : _console << " solution after solve " << std::endl;
251 0 : solution.print();
252 0 : _console << " matrix when we solve " << std::endl;
253 0 : mmat.print();
254 0 : _console << " rhs when we solve " << std::endl;
255 0 : rhs.print();
256 0 : _console << " velocity solution component " << system_i << std::endl;
257 0 : solution.print();
258 0 : _console << "Norm factor " << norm_factor << std::endl;
259 0 : _console << Moose::stringify(momentum_solver.get_initial_residual()) << std::endl;
260 : }
261 :
262 : // Printing residuals
263 261421 : _console << " Momentum equation:"
264 : << (_momentum_systems.size() > 1
265 1115291 : ? std::string(" Component ") + std::to_string(system_i + 1) + std::string(" ")
266 261421 : : std::string(" "))
267 261421 : << COLOR_GREEN << its_normalized_residuals[system_i].second << COLOR_DEFAULT
268 261421 : << " Linear its: " << its_normalized_residuals[system_i].first << std::endl;
269 261421 : }
270 :
271 423538 : for (const auto system_i : index_range(_momentum_systems))
272 : {
273 : LinearImplicitSystem & momentum_system =
274 261421 : libMesh::cast_ref<LinearImplicitSystem &>(_momentum_systems[system_i]->system());
275 261421 : _momentum_systems[system_i]->setSolution(*(momentum_system.current_local_solution));
276 261421 : _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 162117 : momentum_solver.reuse_preconditioner(false);
282 :
283 162117 : return its_normalized_residuals;
284 0 : }
285 :
286 : std::pair<unsigned int, Real>
287 183817 : LinearAssemblySegregatedSolve::solvePressureCorrector()
288 : {
289 183817 : _problem.setCurrentLinearSystem(_pressure_sys_number);
290 :
291 : // We will need some members from the linear system
292 : LinearImplicitSystem & pressure_system =
293 183817 : 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 183817 : SparseMatrix<Number> & mmat = *(pressure_system.matrix);
299 183817 : NumericVector<Number> & rhs = *(pressure_system.rhs);
300 :
301 : // Fetch the linear solver from the system
302 : PetscLinearSolver<Real> & pressure_solver =
303 183817 : libMesh::cast_ref<PetscLinearSolver<Real> &>(*pressure_system.get_linear_solver());
304 :
305 183817 : _problem.computeLinearSystemSys(pressure_system, mmat, rhs, false);
306 :
307 183817 : if (_print_fields)
308 : {
309 0 : _console << "Pressure matrix" << std::endl;
310 0 : mmat.print();
311 : }
312 :
313 : // We compute the normalization factors based on the fluxes
314 183817 : 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 183817 : LibmeshPetscCall(KSPSetNormType(pressure_solver.ksp(), KSP_NORM_UNPRECONDITIONED));
318 :
319 : // Setting the linear tolerances and maximum iteration counts
320 183817 : _pressure_linear_control.real_valued_data["abs_tol"] = _pressure_l_abs_tol * norm_factor;
321 183817 : pressure_solver.set_solver_configuration(_pressure_linear_control);
322 :
323 183817 : if (_pin_pressure)
324 51119 : NS::FV::constrainSystem(mmat, rhs, _pressure_pin_value, _pressure_pin_dof);
325 183817 : pressure_system.update();
326 :
327 183817 : auto its_res_pair = pressure_solver.solve(mmat, mmat, solution, rhs);
328 183817 : pressure_system.update();
329 :
330 183817 : if (_print_fields)
331 : {
332 0 : _console << " rhs when we solve pressure " << std::endl;
333 0 : rhs.print();
334 0 : _console << " Pressure " << std::endl;
335 0 : solution.print();
336 0 : _console << "Norm factor " << norm_factor << std::endl;
337 : }
338 :
339 183817 : _pressure_system.setSolution(current_local_solution);
340 :
341 : const auto residuals =
342 183817 : std::make_pair(its_res_pair.first, pressure_solver.get_initial_residual() / norm_factor);
343 :
344 183817 : _console << " Pressure equation: " << COLOR_GREEN << residuals.second << COLOR_DEFAULT
345 183817 : << " Linear its: " << residuals.first << std::endl;
346 :
347 183817 : return residuals;
348 : }
349 :
350 : std::pair<unsigned int, Real>
351 5058 : LinearAssemblySegregatedSolve::solveSolidEnergy()
352 : {
353 5058 : _problem.setCurrentLinearSystem(_solid_energy_sys_number);
354 :
355 : // We will need some members from the linear system
356 : LinearImplicitSystem & system =
357 5058 : 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 5058 : SparseMatrix<Number> & mmat = *(system.matrix);
363 5058 : NumericVector<Number> & rhs = *(system.rhs);
364 :
365 : // Fetch the linear solver from the system
366 : PetscLinearSolver<Real> & solver =
367 5058 : libMesh::cast_ref<PetscLinearSolver<Real> &>(*system.get_linear_solver());
368 :
369 5058 : _problem.computeLinearSystemSys(system, mmat, rhs, false);
370 :
371 5058 : if (_print_fields)
372 : {
373 0 : _console << "Solid energy matrix" << std::endl;
374 0 : mmat.print();
375 : }
376 :
377 : // We compute the normalization factors based on the fluxes
378 5058 : 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 5058 : LibmeshPetscCall(KSPSetNormType(solver.ksp(), KSP_NORM_UNPRECONDITIONED));
382 :
383 : // Setting the linear tolerances and maximum iteration counts
384 5058 : _solid_energy_linear_control.real_valued_data["abs_tol"] = _solid_energy_l_abs_tol * norm_factor;
385 5058 : solver.set_solver_configuration(_solid_energy_linear_control);
386 :
387 5058 : auto its_res_pair = solver.solve(mmat, mmat, solution, rhs);
388 5058 : system.update();
389 :
390 5058 : if (_print_fields)
391 : {
392 0 : _console << " rhs when we solve solid energy " << std::endl;
393 0 : rhs.print();
394 0 : _console << " Solid energy " << std::endl;
395 0 : solution.print();
396 0 : _console << "Norm factor " << norm_factor << std::endl;
397 : }
398 :
399 5058 : _solid_energy_system->setSolution(current_local_solution);
400 :
401 : const auto residuals =
402 5058 : std::make_pair(its_res_pair.first, solver.get_initial_residual() / norm_factor);
403 :
404 5058 : _console << " Solid energy equation: " << COLOR_GREEN << residuals.second << COLOR_DEFAULT
405 5058 : << " Linear its: " << residuals.first << std::endl;
406 :
407 5058 : return residuals;
408 : }
409 :
410 : std::pair<unsigned int, Real>
411 183817 : 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 183817 : _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.
422 183817 : Moose::PetscSupport::petscSetOptions(_pressure_petsc_options, solver_params);
423 :
424 : // Solve the pressure corrector
425 183817 : 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 183817 : if (recompute_face_mass_flux)
430 162117 : _rc_uo->computeFaceMassFlux();
431 :
432 183817 : auto & pressure_current_solution = *(_pressure_system.system().current_local_solution.get());
433 183817 : auto & pressure_old_solution = *(_pressure_system.solutionPreviousNewton());
434 :
435 : // Relax the pressure update for the next momentum predictor
436 183817 : NS::FV::relaxSolutionUpdate(
437 183817 : pressure_current_solution, pressure_old_solution, _pressure_variable_relaxation);
438 :
439 : // Overwrite old solution
440 183817 : pressure_old_solution = pressure_current_solution;
441 183817 : _pressure_system.setSolution(pressure_current_solution);
442 :
443 : // We recompute the updated pressure gradient
444 183817 : _pressure_system.computeGradients();
445 :
446 : // Reconstruct the cell velocity as well to accelerate convergence
447 183817 : _rc_uo->computeCellVelocity();
448 :
449 183817 : return residuals;
450 : }
451 :
452 : std::pair<unsigned int, Real>
453 105108 : LinearAssemblySegregatedSolve::solveAdvectedSystem(const unsigned int system_num,
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 105108 : _problem.setCurrentLinearSystem(system_num);
462 :
463 : // We will need some members from the implicit linear system
464 105108 : 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 105108 : SparseMatrix<Number> & mmat = *(li_system.matrix);
470 105108 : NumericVector<Number> & rhs = *(li_system.rhs);
471 :
472 : // We need a vector that stores the (diagonal_relaxed-original_diagonal) vector
473 105108 : auto diff_diagonal = solution.zero_clone();
474 :
475 : // Fetch the linear solver from the system
476 : PetscLinearSolver<Real> & linear_solver =
477 105108 : libMesh::cast_ref<PetscLinearSolver<Real> &>(*li_system.get_linear_solver());
478 :
479 105108 : _problem.computeLinearSystemSys(li_system, mmat, rhs, true);
480 :
481 : // Go and relax the system matrix and the right hand side
482 105108 : NS::FV::relaxMatrix(mmat, relaxation_factor, *diff_diagonal);
483 105108 : NS::FV::relaxRightHandSide(rhs, solution, *diff_diagonal);
484 :
485 105108 : if (_print_fields)
486 : {
487 0 : _console << system.name() << " system matrix" << std::endl;
488 0 : mmat.print();
489 : }
490 :
491 : // We compute the normalization factors based on the fluxes
492 105108 : 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 105108 : LibmeshPetscCall(KSPSetNormType(linear_solver.ksp(), KSP_NORM_UNPRECONDITIONED));
496 :
497 : // Setting the linear tolerances and maximum iteration counts
498 105108 : solver_config.real_valued_data["abs_tol"] = absolute_tol * norm_factor;
499 105108 : linear_solver.set_solver_configuration(solver_config);
500 :
501 : // Solve the system and update current local solution
502 105108 : auto its_res_pair = linear_solver.solve(mmat, mmat, solution, rhs);
503 105108 : li_system.update();
504 :
505 105108 : if (_print_fields)
506 : {
507 0 : _console << " rhs when we solve " << system.name() << std::endl;
508 0 : rhs.print();
509 0 : _console << system.name() << " solution " << std::endl;
510 0 : solution.print();
511 0 : _console << " Norm factor " << norm_factor << std::endl;
512 : }
513 :
514 : // Limiting scalar solution
515 105108 : if (min_value_limiter != std::numeric_limits<Real>::min())
516 37298 : NS::FV::limitSolutionUpdate(current_local_solution, min_value_limiter);
517 :
518 : // Relax the field update for the next momentum predictor
519 105108 : if (field_relaxation != 1.0)
520 : {
521 12082 : auto & old_local_solution = *(system.solutionPreviousNewton());
522 12082 : NS::FV::relaxSolutionUpdate(current_local_solution, old_local_solution, field_relaxation);
523 :
524 : // Update old solution, only needed if relaxing the field
525 12082 : old_local_solution = current_local_solution;
526 : }
527 :
528 105108 : system.setSolution(current_local_solution);
529 :
530 : const auto residuals =
531 105108 : std::make_pair(its_res_pair.first, linear_solver.get_initial_residual() / norm_factor);
532 :
533 105108 : _console << " Advected system: " << system.name() << " " << COLOR_GREEN << residuals.second
534 105108 : << COLOR_DEFAULT << " Linear its: " << residuals.first << std::endl;
535 :
536 105108 : return residuals;
537 105108 : }
538 :
539 : bool
540 1634 : LinearAssemblySegregatedSolve::solve()
541 : {
542 : // Do not solve if problem is set not to
543 1634 : if (!_problem.shouldSolve())
544 : return true;
545 :
546 : // Dummy solver parameter file which is needed for switching petsc options
547 1634 : SolverParams solver_params;
548 1634 : solver_params._type = Moose::SolveType::ST_LINEAR;
549 1634 : solver_params._line_search = Moose::LineSearchType::LS_NONE;
550 :
551 : // Initialize the SIMPLE iteration counter
552 1634 : unsigned int simple_iteration_counter = 0;
553 :
554 : // Assign residuals to general residual vector
555 1634 : const unsigned int no_systems = _momentum_systems.size() + 1 + _has_energy_system +
556 1634 : _has_solid_energy_system + _active_scalar_systems.size() +
557 1634 : _turbulence_systems.size();
558 :
559 1634 : std::vector<std::pair<unsigned int, Real>> ns_residuals(no_systems, std::make_pair(0, 1.0));
560 1634 : std::vector<Real> ns_abs_tols(_momentum_systems.size(), _momentum_absolute_tolerance);
561 1634 : ns_abs_tols.push_back(_pressure_absolute_tolerance);
562 :
563 : // Push back energy tolerances
564 1634 : if (_has_energy_system)
565 840 : ns_abs_tols.push_back(_energy_absolute_tolerance);
566 1634 : if (_has_solid_energy_system)
567 18 : ns_abs_tols.push_back(_solid_energy_absolute_tolerance);
568 1634 : if (_has_active_scalar_systems)
569 338 : for (const auto scalar_tol : _active_scalar_absolute_tolerance)
570 169 : ns_abs_tols.push_back(scalar_tol);
571 :
572 : // Push back turbulence tolerances
573 1634 : if (_has_turbulence_systems)
574 378 : for (const auto turbulence_tol : _turbulence_absolute_tolerance)
575 252 : ns_abs_tols.push_back(turbulence_tol);
576 :
577 : bool converged = false;
578 : // Loop until converged or hit the maximum allowed iteration number
579 163751 : while (simple_iteration_counter < _num_iterations && !converged)
580 : {
581 162117 : 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
586 162117 : Moose::PetscSupport::petscSetOptions(_momentum_petsc_options, solver_params);
587 :
588 : // Initialize pressure gradients, after this we just reuse the last ones from each
589 : // iteration
590 162117 : if (simple_iteration_counter == 1)
591 1634 : _pressure_system.computeGradients();
592 :
593 162117 : _console << "Iteration " << simple_iteration_counter << " Initial residual norms:" << std::endl;
594 :
595 : // Solve the momentum predictor step
596 162117 : auto momentum_residual = solveMomentumPredictor();
597 423538 : 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 162117 : 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 162117 : if (_has_energy_system)
608 : {
609 : // We set the preconditioner/controllable parameters through petsc options. Linear
610 : // tolerances will be overridden within the solver.
611 26223 : Moose::PetscSupport::petscSetOptions(_energy_petsc_options, solver_params);
612 26223 : ns_residuals[momentum_residual.size() + _has_energy_system] =
613 26223 : solveAdvectedSystem(_energy_sys_number,
614 26223 : *_energy_system,
615 26223 : _energy_equation_relaxation,
616 : _energy_linear_control,
617 26223 : _energy_l_abs_tol);
618 : }
619 162117 : if (_has_solid_energy_system)
620 : {
621 : // We set the preconditioner/controllable parameters through petsc options. Linear
622 : // tolerances will be overridden within the solver.
623 5058 : Moose::PetscSupport::petscSetOptions(_solid_energy_petsc_options, solver_params);
624 5058 : ns_residuals[momentum_residual.size() + _has_solid_energy_system + _has_energy_system] =
625 5058 : solveSolidEnergy();
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
631 162117 : if (_has_active_scalar_systems)
632 : {
633 24023 : _problem.execute(EXEC_NONLINEAR);
634 :
635 : // We set the preconditioner/controllable parameters through petsc options. Linear
636 : // tolerances will be overridden within the solver.
637 24023 : Moose::PetscSupport::petscSetOptions(_active_scalar_petsc_options, solver_params);
638 48046 : for (const auto i : index_range(_active_scalar_system_names))
639 24023 : ns_residuals[momentum_residual.size() + 1 + _has_energy_system + _has_solid_energy_system +
640 24023 : i] = solveAdvectedSystem(_active_scalar_system_numbers[i],
641 : *_active_scalar_systems[i],
642 : _active_scalar_equation_relaxation[i],
643 : _active_scalar_linear_control,
644 24023 : _active_scalar_l_abs_tol);
645 : }
646 :
647 : // If we have turbulence equations, solve them here.
648 : // The turbulent viscosity depends on the value of the turbulence surrogate variables
649 162117 : if (_has_turbulence_systems)
650 : {
651 : // We set the preconditioner/controllable parameters through petsc options. Linear
652 : // tolerances will be overridden within the solver.
653 18649 : Moose::PetscSupport::petscSetOptions(_turbulence_petsc_options, solver_params);
654 55947 : for (const auto i : index_range(_turbulence_system_names))
655 : {
656 37298 : ns_residuals[momentum_residual.size() + 1 + _has_energy_system + _has_solid_energy_system +
657 37298 : _active_scalar_system_names.size() + i] =
658 37298 : solveAdvectedSystem(_turbulence_system_numbers[i],
659 : *_turbulence_systems[i],
660 : _turbulence_equation_relaxation[i],
661 : _turbulence_linear_control,
662 37298 : _turbulence_l_abs_tol,
663 : _turbulence_field_relaxation[i],
664 : _turbulence_field_min_limit[i]);
665 : }
666 : }
667 :
668 162117 : _problem.execute(EXEC_NONLINEAR);
669 :
670 162117 : converged = NS::FV::converged(ns_residuals, ns_abs_tols);
671 162117 : }
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
676 1634 : if (_has_passive_scalar_systems && (converged || _continue_on_max_its))
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 98 : unsigned int ps_iteration_counter = 0;
682 :
683 98 : _console << "Passive scalar iteration " << ps_iteration_counter
684 98 : << " Initial residual norms:" << std::endl;
685 :
686 8880 : while (ps_iteration_counter < _num_iterations && !passive_scalar_converged)
687 : {
688 8782 : ps_iteration_counter++;
689 : std::vector<std::pair<unsigned int, Real>> scalar_residuals(
690 8782 : _passive_scalar_system_names.size(), std::make_pair(0, 1.0));
691 : std::vector<Real> scalar_abs_tols;
692 26346 : for (const auto scalar_tol : _passive_scalar_absolute_tolerance)
693 17564 : 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.
697 8782 : Moose::PetscSupport::petscSetOptions(_passive_scalar_petsc_options, solver_params);
698 26346 : for (const auto i : index_range(_passive_scalar_system_names))
699 17564 : scalar_residuals[i] = solveAdvectedSystem(_passive_scalar_system_numbers[i],
700 : *_passive_scalar_systems[i],
701 : _passive_scalar_equation_relaxation[i],
702 : _passive_scalar_linear_control,
703 17564 : _passive_scalar_l_abs_tol);
704 :
705 8782 : passive_scalar_converged = NS::FV::converged(scalar_residuals, scalar_abs_tols);
706 8782 : }
707 :
708 : // Both flow and scalars must converge
709 98 : converged = passive_scalar_converged && converged;
710 : }
711 :
712 1634 : converged = _continue_on_max_its ? true : converged;
713 :
714 : return converged;
715 1634 : }
|
