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