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