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 "SIMPLESolveNonlinearAssembly.h"
11 : #include "FEProblem.h"
12 : #include "SegregatedSolverUtils.h"
13 : #include "NonlinearSystemBase.h"
14 :
15 : #include "libmesh/nonlinear_implicit_system.h"
16 :
17 : using namespace libMesh;
18 :
19 : InputParameters
20 1484 : SIMPLESolveNonlinearAssembly::validParams()
21 : {
22 1484 : InputParameters params = SIMPLESolveBase::validParams();
23 :
24 2968 : params.addParam<TagName>("pressure_gradient_tag",
25 : "pressure_momentum_kernels",
26 : "The name of the tags associated with the kernels in the momentum "
27 : "equations which are not related to the pressure gradient.");
28 :
29 1484 : return params;
30 0 : }
31 :
32 742 : SIMPLESolveNonlinearAssembly::SIMPLESolveNonlinearAssembly(Executioner & ex)
33 : : SIMPLESolveBase(ex),
34 1484 : _pressure_sys_number(_problem.nlSysNum(getParam<SolverSystemName>("pressure_system"))),
35 742 : _pressure_system(_problem.getNonlinearSystemBase(_pressure_sys_number)),
36 1484 : _has_turbulence_systems(!getParam<std::vector<SolverSystemName>>("turbulence_systems").empty()),
37 1484 : _energy_sys_number(_has_energy_system
38 1106 : ? _problem.nlSysNum(getParam<SolverSystemName>("energy_system"))
39 : : libMesh::invalid_uint),
40 742 : _energy_system(_has_energy_system ? &_problem.getNonlinearSystemBase(_energy_sys_number)
41 : : nullptr),
42 742 : _solid_energy_sys_number(
43 742 : _has_solid_energy_system
44 760 : ? _problem.nlSysNum(getParam<SolverSystemName>("solid_energy_system"))
45 : : libMesh::invalid_uint),
46 742 : _solid_energy_system(_has_solid_energy_system
47 742 : ? &_problem.getNonlinearSystemBase(_solid_energy_sys_number)
48 : : nullptr),
49 1484 : _turbulence_system_names(getParam<std::vector<SolverSystemName>>("turbulence_systems")),
50 1484 : _turbulence_equation_relaxation(getParam<std::vector<Real>>("turbulence_equation_relaxation")),
51 1484 : _turbulence_field_min_limit(getParam<std::vector<Real>>("turbulence_field_min_limit")),
52 1484 : _turbulence_l_abs_tol(getParam<Real>("turbulence_l_abs_tol")),
53 1484 : _turbulence_absolute_tolerance(getParam<std::vector<Real>>("turbulence_absolute_tolerance")),
54 742 : _pressure_tag_name(getParam<TagName>("pressure_gradient_tag")),
55 2226 : _pressure_tag_id(_problem.addVectorTag(_pressure_tag_name))
56 : {
57 : // We disable this considering that this object passes petsc options a little differently
58 742 : _pressure_system.system().prefix_with_name(false);
59 :
60 : // We fetch the system numbers for the momentum components plus add vectors
61 : // for removing the contribution from the pressure gradient terms.
62 2283 : for (auto system_i : index_range(_momentum_system_names))
63 : {
64 1541 : _momentum_system_numbers.push_back(_problem.nlSysNum(_momentum_system_names[system_i]));
65 1541 : _momentum_systems.push_back(
66 1541 : &_problem.getNonlinearSystemBase(_momentum_system_numbers[system_i]));
67 1541 : _momentum_systems[system_i]->addVector(_pressure_tag_id, false, ParallelType::PARALLEL);
68 :
69 : // We disable this considering that this object passes petsc options a little differently
70 1541 : _momentum_systems[system_i]->system().prefix_with_name(false);
71 : }
72 :
73 742 : if (_has_passive_scalar_systems)
74 156 : for (auto system_i : index_range(_passive_scalar_system_names))
75 : {
76 101 : _passive_scalar_system_numbers.push_back(
77 101 : _problem.nlSysNum(_passive_scalar_system_names[system_i]));
78 101 : _passive_scalar_systems.push_back(
79 101 : &_problem.getNonlinearSystemBase(_passive_scalar_system_numbers[system_i]));
80 :
81 : // We disable this considering that this object passes petsc options a little differently
82 101 : _passive_scalar_systems[system_i]->system().prefix_with_name(false);
83 : }
84 :
85 742 : if (_has_turbulence_systems)
86 : {
87 882 : for (auto system_i : index_range(_turbulence_system_names))
88 : {
89 588 : _turbulence_system_numbers.push_back(_problem.nlSysNum(_turbulence_system_names[system_i]));
90 588 : _turbulence_systems.push_back(
91 588 : &_problem.getNonlinearSystemBase(_turbulence_system_numbers[system_i]));
92 :
93 : // We disable this considering that this object passes petsc options a little differently
94 588 : _turbulence_systems[system_i]->system().prefix_with_name(false);
95 : }
96 :
97 : // We check for input errors with regards to the turbulence equations. At the same time, we
98 : // set up the corresponding system numbers
99 294 : if (_turbulence_system_names.size() != _turbulence_equation_relaxation.size())
100 0 : paramError("turbulence_equation_relaxation",
101 : "The number of equation relaxation parameters does not match the number of "
102 : "turbulence scalar equations!");
103 294 : if (_turbulence_system_names.size() != _turbulence_absolute_tolerance.size())
104 0 : paramError("turbulence_absolute_tolerance",
105 : "The number of absolute tolerances does not match the number of "
106 : "turbulence equations!");
107 294 : if (_turbulence_field_min_limit.empty())
108 : // If no minimum bounds are given, initialize to default value 1e-8
109 294 : _turbulence_field_min_limit.resize(_turbulence_system_names.size(), 1e-8);
110 0 : else if (_turbulence_system_names.size() != _turbulence_field_min_limit.size())
111 0 : paramError("turbulence_field_min_limit",
112 : "The number of lower bounds for turbulent quantities does not match the "
113 : "number of turbulence equations!");
114 : }
115 :
116 2226 : if (isParamValid("solid_energy_system") && !_has_energy_system)
117 0 : paramError(
118 : "solid_energy_system",
119 : "We cannot solve a solid energy system without solving for the fluid energy as well!");
120 :
121 742 : if (_has_turbulence_systems)
122 : {
123 294 : const auto & turbulence_petsc_options = getParam<MultiMooseEnum>("turbulence_petsc_options");
124 : const auto & turbulence_petsc_pair_options = getParam<MooseEnumItem, std::string>(
125 588 : "turbulence_petsc_options_iname", "turbulence_petsc_options_value");
126 294 : Moose::PetscSupport::addPetscFlagsToPetscOptions(
127 : turbulence_petsc_options, "-", *this, _turbulence_petsc_options);
128 588 : Moose::PetscSupport::addPetscPairsToPetscOptions(turbulence_petsc_pair_options,
129 294 : _problem.mesh().dimension(),
130 : "-",
131 : *this,
132 : _turbulence_petsc_options);
133 :
134 588 : _turbulence_linear_control.real_valued_data["rel_tol"] = getParam<Real>("turbulence_l_tol");
135 588 : _turbulence_linear_control.real_valued_data["abs_tol"] = getParam<Real>("turbulence_l_abs_tol");
136 294 : _turbulence_linear_control.int_valued_data["max_its"] =
137 588 : getParam<unsigned int>("turbulence_l_max_its");
138 294 : }
139 : else
140 896 : checkDependentParameterError("turbulence_system",
141 : {"turbulence_petsc_options",
142 : "turbulence_petsc_options_iname",
143 : "turbulence_petsc_options_value",
144 : "turbulence_l_tol",
145 : "turbulence_l_abs_tol",
146 : "turbulence_l_max_its",
147 : "turbulence_equation_relaxation",
148 : "turbulence_absolute_tolerance"},
149 : false);
150 742 : }
151 :
152 : void
153 742 : SIMPLESolveNonlinearAssembly::linkRhieChowUserObject()
154 : {
155 : // Fetch the segregated rhie-chow object and transfer some information about the momentum
156 : // system(s)
157 742 : _rc_uo = const_cast<INSFVRhieChowInterpolatorSegregated *>(
158 742 : &getUserObject<INSFVRhieChowInterpolatorSegregated>("rhie_chow_user_object"));
159 742 : _rc_uo->linkMomentumSystem(_momentum_systems, _momentum_system_numbers, _pressure_tag_id);
160 :
161 : // Initialize the face velocities in the RC object
162 742 : _rc_uo->initFaceVelocities();
163 742 : }
164 :
165 : std::vector<std::pair<unsigned int, Real>>
166 43563 : SIMPLESolveNonlinearAssembly::solveMomentumPredictor()
167 : {
168 : // Temporary storage for the (flux-normalized) residuals form
169 : // different momentum components
170 : std::vector<std::pair<unsigned int, Real>> its_normalized_residuals;
171 :
172 : // We can create this here with the assumption that every momentum component has the same number
173 : // of dofs
174 43563 : auto zero_solution = _momentum_systems[0]->system().current_local_solution->zero_clone();
175 :
176 : // Solve the momentum equations.
177 : // TO DO: These equations are VERY similar. If we can store the differences (things coming from
178 : // BCs for example) separately, it is enough to construct one matrix.
179 135771 : for (const auto system_i : index_range(_momentum_systems))
180 : {
181 92208 : _problem.setCurrentNonlinearSystem(_momentum_system_numbers[system_i]);
182 :
183 : // We will need the right hand side and the solution of the next component
184 : NonlinearImplicitSystem & momentum_system =
185 92208 : libMesh::cast_ref<NonlinearImplicitSystem &>(_momentum_systems[system_i]->system());
186 :
187 : PetscLinearSolver<Real> & momentum_solver =
188 92208 : libMesh::cast_ref<PetscLinearSolver<Real> &>(*momentum_system.get_linear_solver());
189 :
190 : NumericVector<Number> & solution = *(momentum_system.solution);
191 92208 : NumericVector<Number> & rhs = *(momentum_system.rhs);
192 92208 : SparseMatrix<Number> & mmat = *(momentum_system.matrix);
193 :
194 92208 : auto diff_diagonal = solution.zero_clone();
195 :
196 : // We plug zero in this to get the system matrix and the right hand side of the linear problem
197 92208 : _problem.computeResidualAndJacobian(*zero_solution, rhs, mmat);
198 : // Sadly, this returns -b so we multiply with -1
199 92208 : rhs.scale(-1.0);
200 :
201 : // Still need to relax the right hand side with the same vector
202 92208 : NS::FV::relaxMatrix(mmat, _momentum_equation_relaxation, *diff_diagonal);
203 92208 : NS::FV::relaxRightHandSide(rhs, solution, *diff_diagonal);
204 :
205 : // The normalization factor depends on the right hand side so we need to recompute it for this
206 : // component
207 92208 : Real norm_factor = NS::FV::computeNormalizationFactor(solution, mmat, rhs);
208 :
209 : // Very important, for deciding the convergence, we need the unpreconditioned
210 : // norms in the linear solve
211 92208 : LibmeshPetscCall(KSPSetNormType(momentum_solver.ksp(), KSP_NORM_UNPRECONDITIONED));
212 : // Solve this component. We don't update the ghosted solution yet, that will come at the end
213 : // of the corrector step. Also setting the linear tolerances and maximum iteration counts.
214 92208 : _momentum_linear_control.real_valued_data["abs_tol"] = _momentum_l_abs_tol * norm_factor;
215 92208 : momentum_solver.set_solver_configuration(_momentum_linear_control);
216 :
217 : // We solve the equation
218 92208 : auto its_resid_pair = momentum_solver.solve(mmat, mmat, solution, rhs);
219 92208 : momentum_system.update();
220 :
221 : // Save the normalized residual
222 : its_normalized_residuals.push_back(
223 92208 : std::make_pair(its_resid_pair.first, momentum_solver.get_initial_residual() / norm_factor));
224 :
225 92208 : if (_print_fields)
226 : {
227 0 : _console << " matrix when we solve " << std::endl;
228 0 : mmat.print();
229 0 : _console << " rhs when we solve " << std::endl;
230 0 : rhs.print();
231 0 : _console << " velocity solution component " << system_i << std::endl;
232 0 : solution.print();
233 0 : _console << "Norm factor " << norm_factor << std::endl;
234 0 : _console << Moose::stringify(momentum_solver.get_initial_residual()) << std::endl;
235 : }
236 :
237 92208 : _momentum_systems[system_i]->setSolution(*(momentum_system.current_local_solution));
238 92208 : _momentum_systems[system_i]->copyPreviousNonlinearSolutions();
239 92208 : }
240 :
241 43563 : return its_normalized_residuals;
242 43563 : }
243 :
244 : std::pair<unsigned int, Real>
245 43563 : SIMPLESolveNonlinearAssembly::solvePressureCorrector()
246 : {
247 43563 : _problem.setCurrentNonlinearSystem(_pressure_sys_number);
248 :
249 : // We will need some members from the implicit nonlinear system
250 : NonlinearImplicitSystem & pressure_system =
251 43563 : libMesh::cast_ref<NonlinearImplicitSystem &>(_pressure_system.system());
252 :
253 : // We will need the solution, the right hand side and the matrix
254 : NumericVector<Number> & current_local_solution = *(pressure_system.current_local_solution);
255 : NumericVector<Number> & solution = *(pressure_system.solution);
256 43563 : SparseMatrix<Number> & mmat = *(pressure_system.matrix);
257 43563 : NumericVector<Number> & rhs = *(pressure_system.rhs);
258 :
259 : // Fetch the linear solver from the system
260 : PetscLinearSolver<Real> & pressure_solver =
261 43563 : libMesh::cast_ref<PetscLinearSolver<Real> &>(*pressure_system.get_linear_solver());
262 :
263 : // We need a zero vector to be able to emulate the Ax=b system by evaluating the
264 : // residual and jacobian. Unfortunately, this will leave us with the -b on the right hand side
265 : // so we correct it by multiplying it with (-1)
266 43563 : auto zero_solution = current_local_solution.zero_clone();
267 43563 : _problem.computeResidualAndJacobian(*zero_solution, rhs, mmat);
268 43563 : rhs.scale(-1.0);
269 :
270 43563 : if (_print_fields)
271 : {
272 0 : _console << "Pressure matrix" << std::endl;
273 0 : mmat.print();
274 : }
275 :
276 : // We compute the normalization factors based on the fluxes
277 43563 : Real norm_factor = NS::FV::computeNormalizationFactor(solution, mmat, rhs);
278 :
279 : // We need the non-preconditioned norm to be consistent with the norm factor
280 43563 : LibmeshPetscCall(KSPSetNormType(pressure_solver.ksp(), KSP_NORM_UNPRECONDITIONED));
281 :
282 : // Setting the linear tolerances and maximum iteration counts
283 43563 : _pressure_linear_control.real_valued_data["abs_tol"] = _pressure_l_abs_tol * norm_factor;
284 43563 : pressure_solver.set_solver_configuration(_pressure_linear_control);
285 :
286 43563 : if (_pin_pressure)
287 16646 : NS::FV::constrainSystem(mmat, rhs, _pressure_pin_value, _pressure_pin_dof);
288 :
289 43563 : auto its_res_pair = pressure_solver.solve(mmat, mmat, solution, rhs);
290 43563 : pressure_system.update();
291 :
292 43563 : if (_print_fields)
293 : {
294 0 : _console << " rhs when we solve pressure " << std::endl;
295 0 : rhs.print();
296 0 : _console << " Pressure " << std::endl;
297 0 : solution.print();
298 0 : _console << "Norm factor " << norm_factor << std::endl;
299 : }
300 :
301 43563 : _pressure_system.setSolution(current_local_solution);
302 :
303 43563 : return std::make_pair(its_res_pair.first, pressure_solver.get_initial_residual() / norm_factor);
304 43563 : }
305 :
306 : std::pair<unsigned int, Real>
307 35667 : SIMPLESolveNonlinearAssembly::solveAdvectedSystem(const unsigned int system_num,
308 : NonlinearSystemBase & system,
309 : const Real relaxation_factor,
310 : libMesh::SolverConfiguration & solver_config,
311 : const Real absolute_tol)
312 : {
313 35667 : _problem.setCurrentNonlinearSystem(system_num);
314 :
315 : // We will need some members from the implicit nonlinear system
316 : NonlinearImplicitSystem & ni_system =
317 35667 : libMesh::cast_ref<NonlinearImplicitSystem &>(system.system());
318 :
319 : // We will need the solution, the right hand side and the matrix
320 : NumericVector<Number> & current_local_solution = *(ni_system.current_local_solution);
321 : NumericVector<Number> & solution = *(ni_system.solution);
322 35667 : SparseMatrix<Number> & mmat = *(ni_system.matrix);
323 35667 : NumericVector<Number> & rhs = *(ni_system.rhs);
324 :
325 : // We need a vector that stores the (diagonal_relaxed-original_diagonal) vector
326 35667 : auto diff_diagonal = solution.zero_clone();
327 :
328 : // Fetch the linear solver from the system
329 : PetscLinearSolver<Real> & linear_solver =
330 35667 : libMesh::cast_ref<PetscLinearSolver<Real> &>(*ni_system.get_linear_solver());
331 :
332 : // We need a zero vector to be able to emulate the Ax=b system by evaluating the
333 : // residual and jacobian. Unfortunately, this will leave us with the -b on the right hand side
334 : // so we correct it by multiplying it with (-1)
335 35667 : auto zero_solution = current_local_solution.zero_clone();
336 35667 : _problem.computeResidualAndJacobian(*zero_solution, rhs, mmat);
337 35667 : rhs.scale(-1.0);
338 :
339 : // Go and relax the system matrix and the right hand side
340 35667 : NS::FV::relaxMatrix(mmat, relaxation_factor, *diff_diagonal);
341 35667 : NS::FV::relaxRightHandSide(rhs, solution, *diff_diagonal);
342 :
343 35667 : if (_print_fields)
344 : {
345 0 : _console << system.name() << " system matrix" << std::endl;
346 0 : mmat.print();
347 0 : _console << system.name() << " RHS vector" << std::endl;
348 0 : rhs.print();
349 : }
350 :
351 : // We compute the normalization factors based on the fluxes
352 35667 : Real norm_factor = NS::FV::computeNormalizationFactor(solution, mmat, rhs);
353 :
354 : // We need the non-preconditioned norm to be consistent with the norm factor
355 35667 : LibmeshPetscCall(KSPSetNormType(linear_solver.ksp(), KSP_NORM_UNPRECONDITIONED));
356 :
357 : // Setting the linear tolerances and maximum iteration counts
358 35667 : solver_config.real_valued_data["abs_tol"] = absolute_tol * norm_factor;
359 35667 : linear_solver.set_solver_configuration(solver_config);
360 :
361 : // Solve the system and update current local solution
362 35667 : auto its_res_pair = linear_solver.solve(mmat, mmat, solution, rhs);
363 35667 : ni_system.update();
364 :
365 35667 : if (_print_fields)
366 : {
367 0 : _console << " rhs when we solve " << system.name() << std::endl;
368 0 : rhs.print();
369 0 : _console << system.name() << " solution " << std::endl;
370 0 : solution.print();
371 0 : _console << " Norm factor " << norm_factor << std::endl;
372 : }
373 :
374 35667 : system.setSolution(current_local_solution);
375 :
376 71334 : return std::make_pair(its_res_pair.first, linear_solver.get_initial_residual() / norm_factor);
377 35667 : }
378 :
379 : std::pair<unsigned int, Real>
380 666 : SIMPLESolveNonlinearAssembly::solveSolidEnergySystem()
381 : {
382 666 : _problem.setCurrentNonlinearSystem(_solid_energy_sys_number);
383 :
384 : // We will need some members from the implicit nonlinear system
385 : NonlinearImplicitSystem & se_system =
386 666 : libMesh::cast_ref<NonlinearImplicitSystem &>(_solid_energy_system->system());
387 :
388 : // We will need the solution, the right hand side and the matrix
389 : NumericVector<Number> & current_local_solution = *(se_system.current_local_solution);
390 : NumericVector<Number> & solution = *(se_system.solution);
391 666 : SparseMatrix<Number> & mat = *(se_system.matrix);
392 666 : NumericVector<Number> & rhs = *(se_system.rhs);
393 :
394 : // Fetch the linear solver from the system
395 : PetscLinearSolver<Real> & se_solver =
396 666 : libMesh::cast_ref<PetscLinearSolver<Real> &>(*se_system.get_linear_solver());
397 :
398 : // We need a zero vector to be able to emulate the Ax=b system by evaluating the
399 : // residual and jacobian. Unfortunately, this will leave us with the -b on the righ hand side
400 : // so we correct it by multiplying it with (-1)
401 666 : auto zero_solution = current_local_solution.zero_clone();
402 666 : _problem.computeResidualAndJacobian(*zero_solution, rhs, mat);
403 666 : rhs.scale(-1.0);
404 :
405 666 : if (_print_fields)
406 : {
407 0 : _console << "Solid energy matrix" << std::endl;
408 0 : mat.print();
409 : }
410 :
411 : // We compute the normalization factors based on the fluxes
412 666 : Real norm_factor = NS::FV::computeNormalizationFactor(solution, mat, rhs);
413 :
414 : // We need the non-preconditioned norm to be consistent with the norm factor
415 666 : LibmeshPetscCall(KSPSetNormType(se_solver.ksp(), KSP_NORM_UNPRECONDITIONED));
416 :
417 : // Setting the linear tolerances and maximum iteration counts
418 666 : _solid_energy_linear_control.real_valued_data["abs_tol"] = _solid_energy_l_abs_tol * norm_factor;
419 666 : se_solver.set_solver_configuration(_solid_energy_linear_control);
420 :
421 666 : auto its_res_pair = se_solver.solve(mat, mat, solution, rhs);
422 666 : se_system.update();
423 :
424 666 : if (_print_fields)
425 : {
426 0 : _console << " Solid energy rhs " << std::endl;
427 0 : rhs.print();
428 0 : _console << " Solid temperature " << std::endl;
429 0 : solution.print();
430 0 : _console << "Norm factor " << norm_factor << std::endl;
431 : }
432 :
433 666 : _solid_energy_system->setSolution(current_local_solution);
434 :
435 666 : return std::make_pair(its_res_pair.first, se_solver.get_initial_residual() / norm_factor);
436 666 : }
437 :
438 : bool
439 647 : SIMPLESolveNonlinearAssembly::solve()
440 : {
441 : // Dummy solver parameter file which is needed for switching petsc options
442 647 : SolverParams solver_params;
443 647 : solver_params._type = Moose::SolveType::ST_LINEAR;
444 647 : solver_params._line_search = Moose::LineSearchType::LS_NONE;
445 :
446 : // Initialize the quantities which matter in terms of the iteration
447 647 : unsigned int iteration_counter = 0;
448 :
449 : // Assign residuals to general residual vector
450 : unsigned int no_systems =
451 647 : _momentum_systems.size() + 1 + _has_energy_system + _has_solid_energy_system;
452 647 : if (_has_turbulence_systems)
453 219 : no_systems += _turbulence_systems.size();
454 647 : std::vector<std::pair<unsigned int, Real>> ns_its_residuals(no_systems, std::make_pair(0, 1.0));
455 647 : std::vector<Real> ns_abs_tols(_momentum_systems.size(), _momentum_absolute_tolerance);
456 647 : ns_abs_tols.push_back(_pressure_absolute_tolerance);
457 647 : if (_has_energy_system)
458 : {
459 152 : ns_abs_tols.push_back(_energy_absolute_tolerance);
460 152 : if (_has_solid_energy_system)
461 9 : ns_abs_tols.push_back(_solid_energy_absolute_tolerance);
462 : }
463 647 : if (_has_turbulence_systems)
464 657 : for (auto system_i : index_range(_turbulence_absolute_tolerance))
465 438 : ns_abs_tols.push_back(_turbulence_absolute_tolerance[system_i]);
466 :
467 : bool converged = false;
468 : // Loop until converged or hit the maximum allowed iteration number
469 44210 : while (iteration_counter < _num_iterations && !converged)
470 : {
471 43563 : iteration_counter++;
472 : // Resdiual index
473 : size_t residual_index = 0;
474 :
475 : // Execute all objects tagged as nonlinear
476 : // This will execute everything in the problem at nonlinear, including the aux kernels.
477 : // This way we compute the aux kernels before the momentum equations are solved.
478 43563 : _problem.execute(EXEC_NONLINEAR);
479 :
480 : // We clear the caches in the momentum and pressure variables
481 135771 : for (auto system_i : index_range(_momentum_systems))
482 92208 : _momentum_systems[system_i]->residualSetup();
483 43563 : _pressure_system.residualSetup();
484 :
485 : // If we solve for energy, we clear the caches there too
486 43563 : if (_has_energy_system)
487 : {
488 9463 : _energy_system->residualSetup();
489 9463 : if (_has_solid_energy_system)
490 666 : _solid_energy_system->residualSetup();
491 : }
492 :
493 : // If we solve for turbulence, we clear the caches there too
494 43563 : if (_has_turbulence_systems)
495 27495 : for (auto system_i : index_range(_turbulence_systems))
496 18330 : _turbulence_systems[system_i]->residualSetup();
497 :
498 : // We set the preconditioner/controllable parameters through petsc options. Linear
499 : // tolerances will be overridden within the solver. In case of a segregated momentum
500 : // solver, we assume that every velocity component uses the same preconditioner
501 43563 : Moose::PetscSupport::petscSetOptions(_momentum_petsc_options, solver_params);
502 :
503 : // Solve the momentum predictor step
504 43563 : auto momentum_residual = solveMomentumPredictor();
505 135771 : for (const auto system_i : index_range(momentum_residual))
506 : ns_its_residuals[system_i] = momentum_residual[system_i];
507 :
508 : // Compute the coupling fields between the momentum and pressure equations
509 43563 : _rc_uo->computeHbyA(_print_fields);
510 :
511 : // We set the preconditioner/controllable parameters for the pressure equations through
512 : // petsc options. Linear tolerances will be overridden within the solver.
513 43563 : Moose::PetscSupport::petscSetOptions(_pressure_petsc_options, solver_params);
514 :
515 : // Solve the pressure corrector
516 43563 : ns_its_residuals[momentum_residual.size()] = solvePressureCorrector();
517 : // We need this to make sure we evaluate cell gradients for the nonorthogonal correction in
518 : // the face velocity update
519 43563 : _pressure_system.residualSetup();
520 :
521 : // Compute the face velocity which is used in the advection terms
522 43563 : _rc_uo->computeFaceVelocity();
523 :
524 43563 : auto & pressure_current_solution = *(_pressure_system.system().current_local_solution.get());
525 43563 : auto & pressure_old_solution = *(_pressure_system.solutionPreviousNewton());
526 : // Relax the pressure update for the next momentum predictor
527 43563 : NS::FV::relaxSolutionUpdate(
528 43563 : pressure_current_solution, pressure_old_solution, _pressure_variable_relaxation);
529 :
530 : // Overwrite old solution
531 43563 : pressure_old_solution = pressure_current_solution;
532 43563 : _pressure_system.setSolution(pressure_current_solution);
533 :
534 : // We clear out the caches so that the gradients can be computed with the relaxed solution
535 43563 : _pressure_system.residualSetup();
536 :
537 : // Reconstruct the cell velocity as well to accelerate convergence
538 43563 : _rc_uo->computeCellVelocity();
539 :
540 : // Update residual index
541 : residual_index = momentum_residual.size();
542 :
543 : // If we have an energy equation, solve it here. We assume the material properties in the
544 : // Navier-Stokes equations depend on temperature, therefore we can not solve for temperature
545 : // outside of the velocity-pressure loop
546 43563 : if (_has_energy_system)
547 : {
548 : // We set the preconditioner/controllable parameters through petsc options. Linear
549 : // tolerances will be overridden within the solver.
550 9463 : Moose::PetscSupport::petscSetOptions(_energy_petsc_options, solver_params);
551 9463 : residual_index += 1;
552 18926 : ns_its_residuals[residual_index] = solveAdvectedSystem(_energy_sys_number,
553 9463 : *_energy_system,
554 9463 : _energy_equation_relaxation,
555 : _energy_linear_control,
556 9463 : _energy_l_abs_tol);
557 :
558 9463 : if (_has_solid_energy_system)
559 : {
560 : // We set the preconditioner/controllable parameters through petsc options. Linear
561 : // tolerances will be overridden within the solver.
562 666 : Moose::PetscSupport::petscSetOptions(_energy_petsc_options, solver_params);
563 666 : residual_index += 1;
564 666 : ns_its_residuals[residual_index] = solveSolidEnergySystem();
565 : }
566 : }
567 :
568 : // If we have an turbulence equations, we solve it here. We solve it inside the
569 : // momentum-pressure loop because it affects the turbulent viscosity
570 43563 : if (_has_turbulence_systems)
571 : {
572 9165 : Moose::PetscSupport::petscSetOptions(_turbulence_petsc_options, solver_params);
573 :
574 27495 : for (auto system_i : index_range(_turbulence_systems))
575 : {
576 18330 : residual_index += 1;
577 : ns_its_residuals[residual_index] =
578 36660 : solveAdvectedSystem(_turbulence_system_numbers[system_i],
579 : *_turbulence_systems[system_i],
580 : _turbulence_equation_relaxation[system_i],
581 : _turbulence_linear_control,
582 18330 : _turbulence_l_abs_tol);
583 :
584 : auto & current_solution =
585 18330 : *(_turbulence_systems[system_i]->system().current_local_solution.get());
586 18330 : NS::FV::limitSolutionUpdate(current_solution, _turbulence_field_min_limit[system_i]);
587 :
588 : // Relax the turbulence update for the next momentum predictor
589 18330 : auto & old_solution = *(_turbulence_systems[system_i]->solutionPreviousNewton());
590 :
591 : // Relax the pressure update for the next momentum predictor
592 18330 : NS::FV::relaxSolutionUpdate(
593 : current_solution, old_solution, _turbulence_equation_relaxation[system_i]);
594 :
595 : // Overwrite old solution
596 18330 : old_solution = current_solution;
597 18330 : _turbulence_systems[system_i]->setSolution(current_solution);
598 :
599 : // We clear out the caches so that the gradients can be computed with the relaxed solution
600 18330 : _turbulence_systems[system_i]->residualSetup();
601 : }
602 : }
603 :
604 : // Printing residuals
605 : residual_index = 0;
606 43563 : _console << "Iteration " << iteration_counter << " Initial residual norms:" << std::endl;
607 135771 : for (auto system_i : index_range(_momentum_systems))
608 92208 : _console << " Momentum equation:"
609 : << (_momentum_systems.size() > 1
610 461040 : ? std::string(" Component ") + std::to_string(system_i + 1) +
611 276624 : std::string(" ")
612 92208 : : std::string(" "))
613 92208 : << COLOR_GREEN << ns_its_residuals[system_i].second << COLOR_DEFAULT << std::endl;
614 43563 : _console << " Pressure equation: " << COLOR_GREEN
615 43563 : << ns_its_residuals[momentum_residual.size()].second << COLOR_DEFAULT << std::endl;
616 : residual_index = momentum_residual.size();
617 :
618 43563 : if (_has_energy_system)
619 : {
620 9463 : residual_index += 1;
621 9463 : _console << " Energy equation: " << COLOR_GREEN << ns_its_residuals[residual_index].second
622 9463 : << COLOR_DEFAULT << std::endl;
623 9463 : if (_has_solid_energy_system)
624 : {
625 666 : residual_index += 1;
626 666 : _console << " Solid energy equation: " << COLOR_GREEN
627 666 : << ns_its_residuals[residual_index].second << COLOR_DEFAULT << std::endl;
628 : }
629 : }
630 :
631 43563 : if (_has_turbulence_systems)
632 : {
633 9165 : _console << "Turbulence Iteration " << std::endl;
634 27495 : for (auto system_i : index_range(_turbulence_systems))
635 : {
636 18330 : residual_index += 1;
637 18330 : _console << _turbulence_systems[system_i]->name() << " " << COLOR_GREEN
638 18330 : << ns_its_residuals[residual_index].second << COLOR_DEFAULT << std::endl;
639 : }
640 : }
641 :
642 43563 : converged = NS::FV::converged(ns_its_residuals, ns_abs_tols);
643 43563 : }
644 :
645 647 : converged = _continue_on_max_its ? true : converged;
646 :
647 : // Now we solve for the passive scalar equations, they should not influence the solution of the
648 : // system above. The reason why we need more than one iteration is due to the matrix relaxation
649 : // which can be used to stabilize the equations
650 647 : if (_has_passive_scalar_systems)
651 : {
652 50 : _console << " Passive Scalar Iteration " << iteration_counter << std::endl;
653 :
654 : // We set the options used by Petsc (preconditioners etc). We assume that every passive
655 : // scalar equation uses the same options for now.
656 50 : Moose::PetscSupport::petscSetOptions(_passive_scalar_petsc_options, solver_params);
657 :
658 50 : iteration_counter = 0;
659 : std::vector<std::pair<unsigned int, Real>> passive_scalar_residuals(
660 50 : _passive_scalar_systems.size(), std::make_pair(0, 1.0));
661 :
662 : bool passive_scalar_converged =
663 50 : NS::FV::converged(passive_scalar_residuals, _passive_scalar_absolute_tolerance);
664 4185 : while (iteration_counter < _num_iterations && !passive_scalar_converged)
665 : {
666 : // We clear the caches in the passive scalar variables
667 12009 : for (auto system_i : index_range(_passive_scalar_systems))
668 7874 : _passive_scalar_systems[system_i]->residualSetup();
669 :
670 4135 : iteration_counter++;
671 :
672 : // Solve the passive scalar equations
673 12009 : for (auto system_i : index_range(_passive_scalar_systems))
674 : passive_scalar_residuals[system_i] =
675 7874 : solveAdvectedSystem(_passive_scalar_system_numbers[system_i],
676 : *_passive_scalar_systems[system_i],
677 : _passive_scalar_equation_relaxation[system_i],
678 : _passive_scalar_linear_control,
679 7874 : _passive_scalar_l_abs_tol);
680 :
681 4135 : _console << "Iteration " << iteration_counter << " Initial residual norms:" << std::endl;
682 12009 : for (auto system_i : index_range(_passive_scalar_systems))
683 7874 : _console << _passive_scalar_systems[system_i]->name() << " " << COLOR_GREEN
684 7874 : << passive_scalar_residuals[system_i].second << COLOR_DEFAULT << std::endl;
685 :
686 : passive_scalar_converged =
687 4135 : NS::FV::converged(passive_scalar_residuals, _passive_scalar_absolute_tolerance);
688 : }
689 :
690 50 : converged = _continue_on_max_its ? true : passive_scalar_converged;
691 50 : }
692 :
693 647 : return converged;
694 647 : }
695 :
696 : void
697 742 : SIMPLESolveNonlinearAssembly::checkIntegrity()
698 : {
699 : // check to make sure that we don't have any time kernels in this simulation (Steady State)
700 2283 : for (const auto system : _momentum_systems)
701 1541 : checkTimeKernels(*system);
702 :
703 742 : checkTimeKernels(_pressure_system);
704 :
705 742 : if (_has_energy_system)
706 : {
707 182 : checkTimeKernels(*_energy_system);
708 182 : if (_has_solid_energy_system)
709 9 : checkTimeKernels(*_solid_energy_system);
710 : }
711 :
712 742 : if (_has_passive_scalar_systems)
713 156 : for (const auto system : _passive_scalar_systems)
714 101 : checkTimeKernels(*system);
715 :
716 742 : if (_has_turbulence_systems)
717 882 : for (const auto system : _turbulence_systems)
718 588 : checkTimeKernels(*system);
719 742 : }
720 :
721 : void
722 3163 : SIMPLESolveNonlinearAssembly::checkTimeKernels(NonlinearSystemBase & system)
723 : {
724 : // check to make sure that we don't have any time kernels in this simulation (Steady State)
725 3163 : if (system.containsTimeKernel())
726 0 : mooseError("You have specified time kernels in your steady state simulation in system",
727 0 : system.name());
728 3163 : }
|
