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