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 "LinearSystem.h"
11 : #include "AuxiliarySystem.h"
12 : #include "Problem.h"
13 : #include "FEProblem.h"
14 : #include "PetscSupport.h"
15 : #include "Factory.h"
16 : #include "ParallelUniqueId.h"
17 : #include "ComputeLinearFVGreenGaussGradientFaceThread.h"
18 : #include "ComputeLinearFVGreenGaussGradientVolumeThread.h"
19 : #include "ComputeLinearFVElementalThread.h"
20 : #include "ComputeLinearFVFaceThread.h"
21 : #include "DisplacedProblem.h"
22 : #include "Parser.h"
23 : #include "MooseMesh.h"
24 : #include "MooseUtils.h"
25 : #include "MooseApp.h"
26 : #include "TimeIntegrator.h"
27 : #include "Assembly.h"
28 : #include "FloatingPointExceptionGuard.h"
29 : #include "Moose.h"
30 : #include "ConsoleStream.h"
31 : #include "MooseError.h"
32 : #include "LinearFVKernel.h"
33 : #include "UserObject.h"
34 : #include "SolutionInvalidity.h"
35 : #include "MooseLinearVariableFV.h"
36 : #include "LinearFVTimeDerivative.h"
37 :
38 : // libMesh
39 : #include "libmesh/linear_solver.h"
40 : #include "libmesh/quadrature_gauss.h"
41 : #include "libmesh/dense_vector.h"
42 : #include "libmesh/boundary_info.h"
43 : #include "libmesh/petsc_matrix.h"
44 : #include "libmesh/petsc_vector.h"
45 : #include "libmesh/petsc_nonlinear_solver.h"
46 : #include "libmesh/numeric_vector.h"
47 : #include "libmesh/mesh.h"
48 : #include "libmesh/dense_subvector.h"
49 : #include "libmesh/dense_submatrix.h"
50 : #include "libmesh/dof_map.h"
51 : #include "libmesh/sparse_matrix.h"
52 : #include "libmesh/petsc_matrix.h"
53 : #include "libmesh/default_coupling.h"
54 : #include "libmesh/diagonal_matrix.h"
55 : #include "libmesh/petsc_solver_exception.h"
56 :
57 : #include <ios>
58 :
59 : using namespace libMesh;
60 :
61 : namespace Moose
62 : {
63 : void
64 4897 : compute_linear_system(libMesh::EquationSystems & es, const std::string & system_name)
65 : {
66 4897 : FEProblemBase * p = es.parameters.get<FEProblemBase *>("_fe_problem_base");
67 4897 : auto & sys = p->getLinearSystem(p->linearSysNum(system_name));
68 4897 : auto & lin_sys = sys.linearImplicitSystem();
69 4897 : auto & matrix = *(sys.linearImplicitSystem().matrix);
70 4897 : auto & rhs = *(sys.linearImplicitSystem().rhs);
71 4897 : p->computeLinearSystemSys(lin_sys, matrix, rhs);
72 4897 : }
73 : }
74 :
75 1119 : LinearSystem::LinearSystem(FEProblemBase & fe_problem, const std::string & name)
76 : : SolverSystem(fe_problem, fe_problem, name, Moose::VAR_SOLVER),
77 1119 : PerfGraphInterface(fe_problem.getMooseApp().perfGraph(), "LinearSystem"),
78 1119 : _sys(fe_problem.es().add_system<LinearImplicitSystem>(name)),
79 1119 : _rhs_time_tag(-1),
80 1119 : _rhs_time(NULL),
81 1119 : _rhs_non_time_tag(-1),
82 1119 : _rhs_non_time(NULL),
83 1119 : _n_linear_iters(0),
84 1119 : _converged(false),
85 3357 : _linear_implicit_system(fe_problem.es().get_system<LinearImplicitSystem>(name))
86 : {
87 1119 : getRightHandSideNonTimeVector();
88 : // Don't need to add the matrix - it already exists. Well, technically it will exist
89 : // after the initialization. Right now it is just a nullpointer. We will just make sure
90 : // we associate the tag with the system matrix for now.
91 1119 : _system_matrix_tag = _fe_problem.addMatrixTag("SYSTEM");
92 :
93 : // We create a tag for the right hand side, the vector is already in the libmesh system
94 1119 : _rhs_tag = _fe_problem.addVectorTag("RHS");
95 1119 : associateVectorToTag(*_linear_implicit_system.rhs, _rhs_tag);
96 :
97 1119 : _linear_implicit_system.attach_assemble_function(Moose::compute_linear_system);
98 1119 : }
99 :
100 1119 : LinearSystem::~LinearSystem() = default;
101 :
102 : void
103 1119 : LinearSystem::initialSetup()
104 : {
105 1119 : SystemBase::initialSetup();
106 : // Checking if somebody accidentally assigned nonlinear variables to this system
107 1119 : const auto & var_names = _vars[0].names();
108 2238 : for (const auto & name : var_names)
109 1119 : if (!dynamic_cast<MooseLinearVariableFVReal *>(_vars[0].getVariable(name)))
110 0 : mooseError("You are trying to add a nonlinear variable to a linear system! The variable "
111 : "which is assigned to the wrong system: ",
112 : name);
113 1119 : }
114 :
115 : void
116 4957 : LinearSystem::computeLinearSystemTags(const std::set<TagID> & vector_tags,
117 : const std::set<TagID> & matrix_tags,
118 : const bool compute_gradients)
119 : {
120 : parallel_object_only();
121 :
122 4957 : TIME_SECTION("LinearSystem::computeLinearSystemTags", 5);
123 :
124 4957 : _fe_problem.setCurrentLinearSystem(_fe_problem.linearSysNum(name()));
125 :
126 4957 : FloatingPointExceptionGuard fpe_guard(_app);
127 :
128 : try
129 : {
130 4957 : computeLinearSystemInternal(vector_tags, matrix_tags, compute_gradients);
131 : }
132 0 : catch (MooseException & e)
133 : {
134 0 : _console << "Exception detected " << e.what() << std::endl;
135 : // The buck stops here, we have already handled the exception by
136 : // calling stopSolve(), it is now up to PETSc to return a
137 : // "diverged" reason during the next solve.
138 0 : }
139 4957 : }
140 :
141 : void
142 9865 : LinearSystem::computeGradients()
143 : {
144 9865 : _new_gradient.clear();
145 23741 : for (auto & vec : _raw_grad_container)
146 13876 : _new_gradient.push_back(vec->zero_clone());
147 :
148 9865 : TIME_SECTION("LinearVariableFV_Gradients", 3 /*, "Computing Linear FV variable gradients"*/);
149 :
150 : PARALLEL_TRY
151 : {
152 : using FaceInfoRange = StoredRange<MooseMesh::const_face_info_iterator, const FaceInfo *>;
153 19730 : FaceInfoRange face_info_range(_fe_problem.mesh().ownedFaceInfoBegin(),
154 29595 : _fe_problem.mesh().ownedFaceInfoEnd());
155 :
156 : ComputeLinearFVGreenGaussGradientFaceThread gradient_face_thread(
157 9865 : _fe_problem, _fe_problem.linearSysNum(name()));
158 9865 : Threads::parallel_reduce(face_info_range, gradient_face_thread);
159 9865 : }
160 9865 : PARALLEL_CATCH;
161 :
162 23741 : for (auto & vec : _new_gradient)
163 13876 : vec->close();
164 :
165 : PARALLEL_TRY
166 : {
167 : using ElemInfoRange = StoredRange<MooseMesh::const_elem_info_iterator, const ElemInfo *>;
168 19730 : ElemInfoRange elem_info_range(_fe_problem.mesh().ownedElemInfoBegin(),
169 29595 : _fe_problem.mesh().ownedElemInfoEnd());
170 :
171 : ComputeLinearFVGreenGaussGradientVolumeThread gradient_volume_thread(
172 9865 : _fe_problem, _fe_problem.linearSysNum(name()));
173 9865 : Threads::parallel_reduce(elem_info_range, gradient_volume_thread);
174 9865 : }
175 9865 : PARALLEL_CATCH;
176 :
177 23741 : for (const auto i : index_range(_raw_grad_container))
178 : {
179 13876 : _new_gradient[i]->close();
180 13876 : _raw_grad_container[i] = std::move(_new_gradient[i]);
181 : }
182 9865 : }
183 :
184 : void
185 4957 : LinearSystem::computeLinearSystemInternal(const std::set<TagID> & vector_tags,
186 : const std::set<TagID> & matrix_tags,
187 : const bool compute_gradients)
188 : {
189 4957 : TIME_SECTION("computeLinearSystemInternal", 3);
190 :
191 : // Before we assemble we clear up the matrix and the vector
192 4957 : _linear_implicit_system.matrix->zero();
193 4957 : _linear_implicit_system.rhs->zero();
194 :
195 : // Make matrix ready to use
196 4957 : activeAllMatrixTags();
197 :
198 9936 : for (auto tag : matrix_tags)
199 : {
200 4979 : auto & matrix = getMatrix(tag);
201 : // Necessary for speed
202 4979 : if (auto petsc_matrix = dynamic_cast<PetscMatrix<Number> *>(&matrix))
203 : {
204 4979 : LibmeshPetscCall(MatSetOption(petsc_matrix->mat(),
205 : MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
206 : PETSC_TRUE));
207 4979 : if (!_fe_problem.errorOnJacobianNonzeroReallocation())
208 0 : LibmeshPetscCall(
209 : MatSetOption(petsc_matrix->mat(), MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
210 : }
211 : }
212 :
213 4957 : if (compute_gradients)
214 4957 : computeGradients();
215 :
216 : // linear contributions from the domain
217 : PARALLEL_TRY
218 : {
219 4957 : TIME_SECTION("LinearFVKernels_FullSystem", 3 /*, "Computing LinearFVKernels"*/);
220 :
221 : using ElemInfoRange = StoredRange<MooseMesh::const_elem_info_iterator, const ElemInfo *>;
222 9914 : ElemInfoRange elem_info_range(_fe_problem.mesh().ownedElemInfoBegin(),
223 14871 : _fe_problem.mesh().ownedElemInfoEnd());
224 :
225 : using FaceInfoRange = StoredRange<MooseMesh::const_face_info_iterator, const FaceInfo *>;
226 9914 : FaceInfoRange face_info_range(_fe_problem.mesh().ownedFaceInfoBegin(),
227 14871 : _fe_problem.mesh().ownedFaceInfoEnd());
228 :
229 : ComputeLinearFVElementalThread elem_thread(
230 4957 : _fe_problem, this->number(), vector_tags, matrix_tags);
231 4957 : Threads::parallel_reduce(elem_info_range, elem_thread);
232 :
233 : ComputeLinearFVFaceThread face_thread(_fe_problem,
234 : this->number(),
235 : Moose::FV::LinearFVComputationMode::FullSystem,
236 : vector_tags,
237 4957 : matrix_tags);
238 4957 : Threads::parallel_reduce(face_info_range, face_thread);
239 4957 : }
240 4957 : PARALLEL_CATCH;
241 :
242 4957 : closeTaggedMatrices(matrix_tags);
243 :
244 4957 : _linear_implicit_system.matrix->close();
245 4957 : _linear_implicit_system.rhs->close();
246 :
247 : // Accumulate the occurrence of solution invalid warnings for the current iteration cumulative
248 : // counters
249 4957 : _app.solutionInvalidity().syncIteration();
250 4957 : _app.solutionInvalidity().solutionInvalidAccumulation();
251 4957 : }
252 :
253 : NumericVector<Number> &
254 0 : LinearSystem::getRightHandSideTimeVector()
255 : {
256 0 : return *_rhs_time;
257 : }
258 :
259 : NumericVector<Number> &
260 1119 : LinearSystem::getRightHandSideNonTimeVector()
261 : {
262 1119 : return *_rhs_non_time;
263 : }
264 :
265 : void
266 0 : LinearSystem::augmentSparsity(SparsityPattern::Graph & /*sparsity*/,
267 : std::vector<dof_id_type> & /*n_nz*/,
268 : std::vector<dof_id_type> & /*n_oz*/)
269 : {
270 0 : mooseError("LinearSystem does not support AugmentSparsity!");
271 : }
272 :
273 : void
274 4897 : LinearSystem::solve()
275 : {
276 4897 : TIME_SECTION("LinearSystem::solve", 2, "Solving linear system");
277 :
278 : // Clear the iteration counters
279 4897 : _current_l_its = 0;
280 :
281 4897 : system().solve();
282 :
283 : // store info about the solve
284 4897 : _n_linear_iters = _linear_implicit_system.n_linear_iterations();
285 :
286 : auto & linear_solver =
287 4897 : libMesh::cast_ref<PetscLinearSolver<Real> &>(*_linear_implicit_system.get_linear_solver());
288 4897 : _initial_linear_residual = linear_solver.get_initial_residual();
289 4897 : _final_linear_residual = _linear_implicit_system.final_linear_residual();
290 4897 : _converged = linear_solver.get_converged_reason() > 0;
291 :
292 4897 : _console << "System: " << this->name() << " Initial residual: " << _initial_linear_residual
293 4897 : << " Final residual: " << _final_linear_residual << " Num. of Iter. " << _n_linear_iters
294 4897 : << std::endl;
295 :
296 : // determine whether solution invalid occurs in the converged solution
297 4897 : checkInvalidSolution();
298 4897 : }
299 :
300 : void
301 0 : LinearSystem::stopSolve(const ExecFlagType & /*exec_flag*/,
302 : const std::set<TagID> & vector_tags_to_close)
303 : {
304 : // We close the containers in case the solve restarts from a failed iteration
305 0 : closeTaggedVectors(vector_tags_to_close);
306 0 : _linear_implicit_system.matrix->close();
307 0 : }
308 :
309 : bool
310 1095 : LinearSystem::containsTimeKernel()
311 : {
312 : // Right now, FV kernels are in TheWarehouse so we have to use that.
313 1095 : std::vector<LinearFVKernel *> kernels;
314 1095 : auto base_query = _fe_problem.theWarehouse()
315 1095 : .query()
316 2190 : .template condition<AttribSysNum>(this->number())
317 1095 : .template condition<AttribSystem>("LinearFVKernel")
318 1095 : .queryInto(kernels);
319 :
320 1095 : bool contains_time_kernel = false;
321 1095 : for (const auto kernel : kernels)
322 : {
323 0 : contains_time_kernel = dynamic_cast<LinearFVTimeDerivative *>(kernel);
324 0 : if (contains_time_kernel)
325 0 : break;
326 : }
327 :
328 1095 : return contains_time_kernel;
329 1095 : }
330 :
331 : void
332 14806 : LinearSystem::compute(ExecFlagType)
333 : {
334 : // Linear systems have their own time derivative computation machinery
335 14806 : }
|
