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 25708 : compute_linear_system(libMesh::EquationSystems & es, const std::string & system_name)
68 : {
69 25708 : FEProblemBase * p = es.parameters.get<FEProblemBase *>("_fe_problem_base");
70 25708 : auto & sys = p->getLinearSystem(p->linearSysNum(system_name));
71 25708 : auto & lin_sys = sys.linearImplicitSystem();
72 25708 : auto & matrix = *(sys.linearImplicitSystem().matrix);
73 25708 : auto & rhs = *(sys.linearImplicitSystem().rhs);
74 25708 : p->computeLinearSystemSys(lin_sys, matrix, rhs);
75 25708 : }
76 : }
77 :
78 1264 : LinearSystem::LinearSystem(FEProblemBase & fe_problem, const std::string & name)
79 : : SolverSystem(fe_problem, fe_problem, name, Moose::VAR_SOLVER),
80 1264 : PerfGraphInterface(fe_problem.getMooseApp().perfGraph(), "LinearSystem"),
81 : LinearFVGradientInterface(static_cast<SystemBase &>(*this)),
82 1264 : _sys(fe_problem.es().add_system<LinearImplicitSystem>(name)),
83 1264 : _rhs_time_tag(-1),
84 1264 : _rhs_time(NULL),
85 1264 : _rhs_non_time_tag(-1),
86 1264 : _rhs_non_time(NULL),
87 1264 : _n_linear_iters(0),
88 1264 : _converged(false),
89 6320 : _linear_implicit_system(fe_problem.es().get_system<LinearImplicitSystem>(name))
90 : {
91 1264 : 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 1264 : _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 1264 : _rhs_tag = _fe_problem.addVectorTag("RHS");
99 1264 : associateVectorToTag(*_linear_implicit_system.rhs, _rhs_tag);
100 :
101 1264 : _linear_implicit_system.attach_assemble_function(Moose::compute_linear_system);
102 1264 : }
103 :
104 1262 : LinearSystem::~LinearSystem() = default;
105 :
106 : void
107 1228 : LinearSystem::preInit()
108 : {
109 1228 : SolverSystem::preInit();
110 :
111 : #ifdef MOOSE_KOKKOS_ENABLED
112 899 : if (_fe_problem.hasKokkosResidualObjects())
113 181 : _sys.get_dof_map().full_sparsity_pattern_needed();
114 : #endif
115 1228 : }
116 :
117 : void
118 1226 : LinearSystem::initialSetup()
119 : {
120 1226 : SystemBase::initialSetup();
121 1226 : _current_solution = system().current_local_solution.get();
122 : // Checking if somebody accidentally assigned nonlinear variables to this system
123 1226 : const auto & var_names = _vars[0].names();
124 2472 : for (const auto & name : var_names)
125 1246 : if (!dynamic_cast<MooseLinearVariableFVReal *>(_vars[0].getVariable(name)))
126 0 : mooseError("You are trying to add a nonlinear variable to a linear system! The variable "
127 : "which is assigned to the wrong system: ",
128 : name);
129 :
130 1226 : LinearFVGradientInterface::rebuildLinearFVGradientStorage();
131 :
132 : // Calling initial setup for the linear kernels
133 2452 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
134 : {
135 1226 : std::vector<LinearFVElementalKernel *> fv_elemental_kernels;
136 1226 : _fe_problem.theWarehouse()
137 1226 : .query()
138 1226 : .template condition<AttribSysNum>(number())
139 1226 : .template condition<AttribSystem>("LinearFVElementalKernel")
140 2452 : .template condition<AttribKokkos>(false)
141 1226 : .template condition<AttribThread>(tid)
142 1226 : .queryInto(fv_elemental_kernels);
143 :
144 2265 : for (auto * fv_kernel : fv_elemental_kernels)
145 1039 : fv_kernel->initialSetup();
146 :
147 1226 : std::vector<LinearFVFluxKernel *> fv_flux_kernels;
148 1226 : _fe_problem.theWarehouse()
149 1226 : .query()
150 1226 : .template condition<AttribSysNum>(number())
151 1226 : .template condition<AttribSystem>("LinearFVFluxKernel")
152 2452 : .template condition<AttribKokkos>(false)
153 1226 : .template condition<AttribThread>(tid)
154 1226 : .queryInto(fv_flux_kernels);
155 :
156 2383 : for (auto * fv_kernel : fv_flux_kernels)
157 1157 : fv_kernel->initialSetup();
158 :
159 1226 : std::vector<LinearFVBoundaryCondition *> fv_bcs;
160 1226 : _fe_problem.theWarehouse()
161 1226 : .query()
162 1226 : .template condition<AttribSysNum>(number())
163 1226 : .template condition<AttribSystem>("LinearFVBoundaryCondition")
164 2452 : .template condition<AttribKokkos>(false)
165 1226 : .template condition<AttribThread>(tid)
166 1226 : .queryInto(fv_bcs);
167 :
168 3290 : for (auto * fv_bc : fv_bcs)
169 2064 : fv_bc->initialSetup();
170 1226 : }
171 :
172 : #ifdef MOOSE_KOKKOS_ENABLED
173 898 : initialSetupKokkosLinearFV();
174 : #endif
175 1226 : }
176 :
177 : void
178 10 : LinearSystem::reinit()
179 : {
180 10 : _current_solution = system().current_local_solution.get();
181 10 : LinearFVGradientInterface::rebuildLinearFVGradientStorage();
182 10 : }
183 :
184 : void
185 25792 : LinearSystem::computeLinearSystemTags(const std::set<TagID> & vector_tags,
186 : const std::set<TagID> & matrix_tags,
187 : const bool compute_gradients)
188 : {
189 : parallel_object_only();
190 :
191 77376 : TIME_SECTION("LinearSystem::computeLinearSystemTags", 5);
192 :
193 25792 : _fe_problem.setCurrentLinearSystem(_fe_problem.linearSysNum(name()));
194 :
195 25792 : FloatingPointExceptionGuard fpe_guard(_app);
196 :
197 : try
198 : {
199 25792 : computeLinearSystemInternal(vector_tags, matrix_tags, compute_gradients);
200 : }
201 0 : catch (MooseException & e)
202 : {
203 0 : _console << "Exception detected " << e.what() << std::endl;
204 : // The buck stops here, we have already handled the exception by
205 : // calling stopSolve(), it is now up to PETSc to return a
206 : // "diverged" reason during the next solve.
207 0 : }
208 25792 : }
209 :
210 : void
211 25792 : LinearSystem::computeLinearSystemInternal(const std::set<TagID> & vector_tags,
212 : const std::set<TagID> & matrix_tags,
213 : const bool compute_gradients)
214 : {
215 77376 : TIME_SECTION("computeLinearSystemInternal", 3);
216 :
217 : // Before we assemble we clear up the matrix and the vector
218 25792 : _linear_implicit_system.matrix->zero();
219 25792 : _linear_implicit_system.rhs->zero();
220 :
221 : // Make matrix ready to use
222 25792 : activateAllMatrixTags();
223 :
224 51604 : for (auto tag : matrix_tags)
225 : {
226 25812 : auto & matrix = getMatrix(tag);
227 : // Necessary for speed
228 25812 : if (auto petsc_matrix = dynamic_cast<PetscMatrix<Number> *>(&matrix))
229 : {
230 25812 : LibmeshPetscCall(MatSetOption(petsc_matrix->mat(),
231 : MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
232 : PETSC_TRUE));
233 25812 : if (!_fe_problem.errorOnJacobianNonzeroReallocation())
234 0 : LibmeshPetscCall(
235 : MatSetOption(petsc_matrix->mat(), MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
236 : }
237 : }
238 :
239 25792 : if (compute_gradients)
240 25792 : computeGradients();
241 :
242 : #ifdef MOOSE_KOKKOS_ENABLED
243 : // Kokkos assembly runs first: it accumulates via device-side atomics and syncs
244 : // back to PETSc before the CPU path writes. The CPU queries below filter out
245 : // Kokkos objects via AttribKokkos(false) so there is no double-counting.
246 17842 : if (_fe_problem.hasKokkosResidualObjects())
247 292 : computeKokkosLinearSystem(vector_tags, matrix_tags);
248 : #endif
249 :
250 : // linear contributions from the domain
251 : PARALLEL_TRY
252 : {
253 77376 : TIME_SECTION("LinearFVKernels_FullSystem", 3 /*, "Computing LinearFVKernels"*/);
254 :
255 : using ElemInfoRange = StoredRange<MooseMesh::const_elem_info_iterator, const ElemInfo *>;
256 51584 : ElemInfoRange elem_info_range(_fe_problem.mesh().ownedElemInfoBegin(),
257 77376 : _fe_problem.mesh().ownedElemInfoEnd());
258 :
259 : using FaceInfoRange = StoredRange<MooseMesh::const_face_info_iterator, const FaceInfo *>;
260 51584 : FaceInfoRange face_info_range(_fe_problem.mesh().ownedFaceInfoBegin(),
261 77376 : _fe_problem.mesh().ownedFaceInfoEnd());
262 :
263 : ComputeLinearFVElementalThread elem_thread(
264 25792 : _fe_problem, this->number(), vector_tags, matrix_tags);
265 25792 : Threads::parallel_reduce(elem_info_range, elem_thread);
266 :
267 : ComputeLinearFVFaceThread face_thread(_fe_problem,
268 : this->number(),
269 : Moose::FV::LinearFVComputationMode::FullSystem,
270 : vector_tags,
271 25792 : matrix_tags);
272 25792 : Threads::parallel_reduce(face_info_range, face_thread);
273 25792 : }
274 25792 : PARALLEL_CATCH;
275 :
276 25792 : closeTaggedMatrices(matrix_tags);
277 :
278 25792 : _linear_implicit_system.matrix->close();
279 25792 : _linear_implicit_system.rhs->close();
280 :
281 : // Accumulate the occurrence of solution invalid warnings for the current iteration cumulative
282 : // counters
283 25792 : _app.solutionInvalidity().syncIteration();
284 25792 : _app.solutionInvalidity().accumulateIterationIntoTimeStepOccurences();
285 25792 : }
286 :
287 : NumericVector<Number> &
288 0 : LinearSystem::getRightHandSideTimeVector()
289 : {
290 0 : return *_rhs_time;
291 : }
292 :
293 : NumericVector<Number> &
294 1264 : LinearSystem::getRightHandSideNonTimeVector()
295 : {
296 1264 : return *_rhs_non_time;
297 : }
298 :
299 : void
300 0 : LinearSystem::augmentSparsity(SparsityPattern::Graph & /*sparsity*/,
301 : std::vector<dof_id_type> & /*n_nz*/,
302 : std::vector<dof_id_type> & /*n_oz*/)
303 : {
304 0 : mooseError("LinearSystem does not support AugmentSparsity!");
305 : }
306 :
307 : void
308 25708 : LinearSystem::solve()
309 : {
310 128540 : TIME_SECTION("LinearSystem::solve", 2, "Solving linear system");
311 :
312 : // Clear the iteration counters
313 25708 : _current_l_its = 0;
314 :
315 25708 : system().solve();
316 :
317 : // store info about the solve
318 25708 : _n_linear_iters = _linear_implicit_system.n_linear_iterations();
319 :
320 : auto & linear_solver =
321 25708 : libMesh::cast_ref<PetscLinearSolver<Real> &>(*_linear_implicit_system.get_linear_solver());
322 25708 : _initial_linear_residual = linear_solver.get_initial_residual();
323 25708 : _final_linear_residual = _linear_implicit_system.final_linear_residual();
324 25708 : _converged = linear_solver.get_converged_reason() > 0;
325 :
326 25708 : _console << "System: " << this->name() << " Initial residual: " << _initial_linear_residual
327 25708 : << " Final residual: " << _final_linear_residual << " Num. of Iter. " << _n_linear_iters
328 25708 : << std::endl;
329 :
330 : // determine whether solution invalid occurs in the converged solution
331 25708 : checkInvalidSolution();
332 25708 : }
333 :
334 : void
335 0 : LinearSystem::stopSolve(const ExecFlagType & /*exec_flag*/,
336 : const std::set<TagID> & vector_tags_to_close)
337 : {
338 : // We close the containers in case the solve restarts from a failed iteration
339 0 : closeTaggedVectors(vector_tags_to_close);
340 0 : _linear_implicit_system.matrix->close();
341 0 : }
342 :
343 : bool
344 1197 : LinearSystem::containsTimeKernel()
345 : {
346 : // Right now, FV kernels are in TheWarehouse so we have to use that.
347 1197 : std::vector<LinearFVKernel *> kernels;
348 1197 : auto base_query = _fe_problem.theWarehouse()
349 1197 : .query()
350 2394 : .template condition<AttribSysNum>(this->number())
351 1197 : .template condition<AttribSystem>("LinearFVKernel")
352 1197 : .queryInto(kernels);
353 :
354 1197 : bool contains_time_kernel = false;
355 1197 : for (const auto kernel : kernels)
356 : {
357 0 : contains_time_kernel = dynamic_cast<LinearFVTimeDerivative *>(kernel);
358 0 : if (contains_time_kernel)
359 0 : break;
360 : }
361 :
362 1197 : return contains_time_kernel;
363 1197 : }
364 :
365 : void
366 67031 : LinearSystem::compute(const ExecFlagType type)
367 : {
368 : // - Linear system assembly is associated with EXEC_NONLINEAR
369 : // - Avoid division by 0 dt
370 67031 : if (type == EXEC_NONLINEAR && _fe_problem.dt() > 0.)
371 390 : for (auto & ti : _time_integrators)
372 : // Do things like compute integration weights
373 195 : ti->preStep();
374 67031 : }
|
