https://mooseframework.inl.gov
LinearSystem.C
Go to the documentation of this file.
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"
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"
29 #include "Moose.h"
30 #include "ConsoleStream.h"
31 #include "MooseError.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 compute_linear_system(libMesh::EquationSystems & es, const std::string & system_name)
66 {
67  FEProblemBase * p = es.parameters.get<FEProblemBase *>("_fe_problem_base");
68  auto & sys = p->getLinearSystem(p->linearSysNum(system_name));
69  auto & lin_sys = sys.linearImplicitSystem();
70  auto & matrix = *(sys.linearImplicitSystem().matrix);
71  auto & rhs = *(sys.linearImplicitSystem().rhs);
72  p->computeLinearSystemSys(lin_sys, matrix, rhs);
73 }
74 }
75 
76 LinearSystem::LinearSystem(FEProblemBase & fe_problem, const std::string & name)
77  : SolverSystem(fe_problem, fe_problem, name, Moose::VAR_SOLVER),
78  PerfGraphInterface(fe_problem.getMooseApp().perfGraph(), "LinearSystem"),
79  _sys(fe_problem.es().add_system<LinearImplicitSystem>(name)),
80  _rhs_time_tag(-1),
81  _rhs_time(NULL),
82  _rhs_non_time_tag(-1),
83  _rhs_non_time(NULL),
84  _n_linear_iters(0),
85  _converged(false),
86  _linear_implicit_system(fe_problem.es().get_system<LinearImplicitSystem>(name))
87 {
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.
93 
94  // We create a tag for the right hand side, the vector is already in the libmesh system
97 
99 }
100 
101 LinearSystem::~LinearSystem() = default;
102 
103 void
105 {
107  // Checking if somebody accidentally assigned nonlinear variables to this system
108  const auto & var_names = _vars[0].names();
109  for (const auto & name : var_names)
110  if (!dynamic_cast<MooseLinearVariableFVReal *>(_vars[0].getVariable(name)))
111  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  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
117  {
118  std::vector<LinearFVElementalKernel *> fv_elemental_kernels;
120  .query()
121  .template condition<AttribSystem>("LinearFVElementalKernel")
122  .template condition<AttribThread>(tid)
123  .queryInto(fv_elemental_kernels);
124 
125  for (auto * fv_kernel : fv_elemental_kernels)
126  fv_kernel->initialSetup();
127 
128  std::vector<LinearFVFluxKernel *> fv_flux_kernels;
130  .query()
131  .template condition<AttribSystem>("LinearFVFluxKernel")
132  .template condition<AttribThread>(tid)
133  .queryInto(fv_flux_kernels);
134 
135  for (auto * fv_kernel : fv_flux_kernels)
136  fv_kernel->initialSetup();
137  }
138 }
139 
140 void
141 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  TIME_SECTION("LinearSystem::computeLinearSystemTags", 5);
148 
150 
152 
153  try
154  {
155  computeLinearSystemInternal(vector_tags, matrix_tags, compute_gradients);
156  }
157  catch (MooseException & e)
158  {
159  _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  }
164 }
165 
166 void
168 {
169  _new_gradient.clear();
170  for (auto & vec : _raw_grad_container)
171  _new_gradient.push_back(vec->zero_clone());
172 
173  TIME_SECTION("LinearVariableFV_Gradients", 3 /*, "Computing Linear FV variable gradients"*/);
174 
175  PARALLEL_TRY
176  {
178  FaceInfoRange face_info_range(_fe_problem.mesh().ownedFaceInfoBegin(),
180 
181  ComputeLinearFVGreenGaussGradientFaceThread gradient_face_thread(
183  Threads::parallel_reduce(face_info_range, gradient_face_thread);
184  }
185  PARALLEL_CATCH;
186 
187  for (auto & vec : _new_gradient)
188  vec->close();
189 
190  PARALLEL_TRY
191  {
193  ElemInfoRange elem_info_range(_fe_problem.mesh().ownedElemInfoBegin(),
195 
196  ComputeLinearFVGreenGaussGradientVolumeThread gradient_volume_thread(
198  Threads::parallel_reduce(elem_info_range, gradient_volume_thread);
199  }
200  PARALLEL_CATCH;
201 
202  for (const auto i : index_range(_raw_grad_container))
203  {
204  _new_gradient[i]->close();
205  _raw_grad_container[i] = std::move(_new_gradient[i]);
206  }
207 }
208 
209 void
210 LinearSystem::computeLinearSystemInternal(const std::set<TagID> & vector_tags,
211  const std::set<TagID> & matrix_tags,
212  const bool compute_gradients)
213 {
214  TIME_SECTION("computeLinearSystemInternal", 3);
215 
216  // Before we assemble we clear up the matrix and the vector
219 
220  // Make matrix ready to use
222 
223  for (auto tag : matrix_tags)
224  {
225  auto & matrix = getMatrix(tag);
226  // Necessary for speed
227  if (auto petsc_matrix = dynamic_cast<PetscMatrix<Number> *>(&matrix))
228  {
229  LibmeshPetscCall(MatSetOption(petsc_matrix->mat(),
230  MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
231  PETSC_TRUE));
233  LibmeshPetscCall(
234  MatSetOption(petsc_matrix->mat(), MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
235  }
236  }
237 
238  if (compute_gradients)
240 
241  // linear contributions from the domain
242  PARALLEL_TRY
243  {
244  TIME_SECTION("LinearFVKernels_FullSystem", 3 /*, "Computing LinearFVKernels"*/);
245 
247  ElemInfoRange elem_info_range(_fe_problem.mesh().ownedElemInfoBegin(),
249 
251  FaceInfoRange face_info_range(_fe_problem.mesh().ownedFaceInfoBegin(),
253 
254  ComputeLinearFVElementalThread elem_thread(
255  _fe_problem, this->number(), vector_tags, matrix_tags);
256  Threads::parallel_reduce(elem_info_range, elem_thread);
257 
259  this->number(),
261  vector_tags,
262  matrix_tags);
263  Threads::parallel_reduce(face_info_range, face_thread);
264  }
265  PARALLEL_CATCH;
266 
267  closeTaggedMatrices(matrix_tags);
268 
271 
272  // Accumulate the occurrence of solution invalid warnings for the current iteration cumulative
273  // counters
276 }
277 
280 {
281  return *_rhs_time;
282 }
283 
286 {
287  return *_rhs_non_time;
288 }
289 
290 void
292  std::vector<dof_id_type> & /*n_nz*/,
293  std::vector<dof_id_type> & /*n_oz*/)
294 {
295  mooseError("LinearSystem does not support AugmentSparsity!");
296 }
297 
298 void
300 {
301  TIME_SECTION("LinearSystem::solve", 2, "Solving linear system");
302 
303  // Clear the iteration counters
304  _current_l_its = 0;
305 
306  system().solve();
307 
308  // store info about the solve
310 
311  auto & linear_solver =
312  libMesh::cast_ref<PetscLinearSolver<Real> &>(*_linear_implicit_system.get_linear_solver());
313  _initial_linear_residual = linear_solver.get_initial_residual();
315  _converged = linear_solver.get_converged_reason() > 0;
316 
317  _console << "System: " << this->name() << " Initial residual: " << _initial_linear_residual
318  << " Final residual: " << _final_linear_residual << " Num. of Iter. " << _n_linear_iters
319  << std::endl;
320 
321  // determine whether solution invalid occurs in the converged solution
323 }
324 
325 void
326 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  closeTaggedVectors(vector_tags_to_close);
332 }
333 
334 bool
336 {
337  // Right now, FV kernels are in TheWarehouse so we have to use that.
338  std::vector<LinearFVKernel *> kernels;
339  auto base_query = _fe_problem.theWarehouse()
340  .query()
341  .template condition<AttribSysNum>(this->number())
342  .template condition<AttribSystem>("LinearFVKernel")
343  .queryInto(kernels);
344 
345  bool contains_time_kernel = false;
346  for (const auto kernel : kernels)
347  {
348  contains_time_kernel = dynamic_cast<LinearFVTimeDerivative *>(kernel);
349  if (contains_time_kernel)
350  break;
351  }
352 
353  return contains_time_kernel;
354 }
355 
356 void
358 {
359  // Linear systems have their own time derivative computation machinery
360 }
std::string name(const ElemQuality q)
virtual void stopSolve(const ExecFlagType &exec_flag, const std::set< TagID > &vector_tags_to_close) override
Quit the current solve as soon as possible.
Definition: LinearSystem.C:326
The gradient in a volume using Green Gauss theorem and a cell-centered finite-volume approximation ca...
virtual const char * what() const
Get out the error message.
unsigned int n_threads()
void computeGradients()
Compute the Green-Gauss gradients.
Definition: LinearSystem.C:167
unsigned int _n_linear_iters
Number of linear iterations.
Definition: LinearSystem.h:190
void checkInvalidSolution()
Definition: SolverSystem.C:111
face_info_iterator ownedFaceInfoBegin()
Iterators to owned faceInfo objects.
Definition: MooseMesh.C:1507
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
virtual bool containsTimeKernel() override
If the system has a kernel that corresponds to a time derivative.
Definition: LinearSystem.C:335
The gradient in a volume using Green Gauss theorem and a cell-centered finite-volume approximation ca...
NumericVector< Number > & getRightHandSideNonTimeVector()
Return a numeric vector that is associated with the nontime tag.
Definition: LinearSystem.C:285
virtual TagID addVectorTag(const TagName &tag_name, const Moose::VectorTagType type=Moose::VECTOR_TAG_RESIDUAL)
Create a Tag.
Definition: SubProblem.C:92
std::vector< T * > & queryInto(std::vector< T *> &results, Args &&... args)
queryInto executes the query and stores the results in the given vector.
Definition: TheWarehouse.h:311
virtual void associateVectorToTag(NumericVector< Number > &vec, TagID tag)
Associate a vector for a given tag.
Definition: SystemBase.C:964
void computeLinearSystemInternal(const std::set< TagID > &vector_tags, const std::set< TagID > &matrix_tags, const bool compute_gradients=true)
Compute the right hand side and system matrix for given tags.
Definition: LinearSystem.C:210
virtual void compute(ExecFlagType type) override
Compute time derivatives, auxiliary variables, etc.
Definition: LinearSystem.C:357
void compute_linear_system(libMesh::EquationSystems &es, const std::string &system_name)
Definition: LinearSystem.C:65
NumericVector< Number > * rhs
virtual LinearSolver< Number > * get_linear_solver() const override
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
unsigned int _current_l_its
The linear iterations needed for convergence.
Definition: LinearSystem.h:160
std::vector< std::unique_ptr< NumericVector< Number > > > _raw_grad_container
A cache for storing gradients at dof locations.
Definition: SystemBase.h:1065
Real _initial_linear_residual
The initial linear residual.
Definition: LinearSystem.h:193
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
std::vector< std::unique_ptr< NumericVector< Number > > > _new_gradient
Vectors to store the new gradients during the computation.
Definition: LinearSystem.h:208
Scope guard for starting and stopping Floating Point Exception Trapping.
virtual void zero()=0
elem_info_iterator ownedElemInfoBegin()
Iterators to owned faceInfo objects.
Definition: MooseMesh.C:1525
Adds contributions from volumetric terms discretized using the finite volume method to the matrix and...
virtual void initialSetup() override
Setup Functions.
Definition: LinearSystem.C:104
virtual void computeLinearSystemSys(libMesh::LinearImplicitSystem &sys, libMesh::SparseMatrix< libMesh::Number > &system_matrix, NumericVector< libMesh::Number > &rhs, const bool compute_gradients=true)
Assemble both the right hand side and the system matrix of a given linear system. ...
virtual const std::string & name() const
Definition: SystemBase.C:1330
void closeTaggedMatrices(const std::set< TagID > &tags)
Close all matrices associated the tags.
Definition: SystemBase.C:1043
void solutionInvalidAccumulation()
Pass the number of solution invalid occurrences from current iteration to cumulative counters...
void syncIteration()
Sync iteration counts to main processor.
libMesh::LinearImplicitSystem & _linear_implicit_system
Base class reference to the linear implicit system in libmesh.
Definition: LinearSystem.h:202
NumericVector< Number > * _rhs_time
right hand side vector for time contributions
Definition: LinearSystem.h:172
const T & get(std::string_view) const
virtual void zero()=0
TheWarehouse & theWarehouse() const
unsigned int n_linear_iterations() const
Real _final_linear_residual
The final linear residual.
Definition: LinearSystem.h:196
SolutionInvalidity & solutionInvalidity()
Get the SolutionInvalidity for this app.
Definition: MooseApp.h:178
virtual TagID addMatrixTag(TagName tag_name)
Create a Tag.
Definition: SubProblem.C:311
TagID _rhs_tag
Used for the right hand side vector from PETSc.
Definition: LinearSystem.h:181
bool errorOnJacobianNonzeroReallocation() const
Will return True if the user wants to get an error when a nonzero is reallocated in the Jacobian by P...
virtual void activeAllMatrixTags()
Make all exsiting matrices ative.
Definition: SystemBase.C:1137
unsigned int number() const
Gets the number of this system.
Definition: SystemBase.C:1159
void attach_assemble_function(void fptr(EquationSystems &es, const std::string &name))
virtual void solve() override
Solve the system (using libMesh magic)
Definition: LinearSystem.C:299
void closeTaggedVectors(const std::set< TagID > &tags)
Close all vectors for given tags.
Definition: SystemBase.C:659
Interface for objects interacting with the PerfGraph.
virtual void solve()
virtual void close()=0
LinearSystem & getLinearSystem(unsigned int sys_num)
Get non-constant reference to a linear system.
Kernel that adds contributions from a time derivative term to a linear system populated using the fin...
TagID _system_matrix_tag
Tag for every contribution to system matrix.
Definition: LinearSystem.h:187
MooseApp & _app
Definition: SystemBase.h:982
FEProblemBase & _fe_problem
the governing finite element/volume problem
Definition: SystemBase.h:980
bool _converged
If the solve on the linear system converged.
Definition: LinearSystem.h:199
std::vector< VariableWarehouse > _vars
Variable warehouses (one for each thread)
Definition: SystemBase.h:990
Provides a way for users to bail out of the current solve.
NumericVector< Number > * _rhs_non_time
right hand side vector for non-time contributions
Definition: LinearSystem.h:178
virtual void close()=0
Class for containing MooseEnum item information.
Definition: MooseEnumItem.h:18
NumericVector< Number > & getRightHandSideTimeVector()
Return a numeric vector that is associated with the time tag.
Definition: LinearSystem.C:279
virtual libMesh::SparseMatrix< Number > & getMatrix(TagID tag)
Get a raw SparseMatrix.
Definition: SystemBase.C:1007
SparseMatrix< Number > * matrix
void setCurrentLinearSystem(unsigned int sys_num)
Set the current linear system pointer.
Query query()
query creates and returns an initialized a query object for querying objects from the warehouse...
Definition: TheWarehouse.h:466
virtual MooseMesh & mesh() override
unsigned int linearSysNum(const LinearSystemName &linear_sys_name) const override
virtual void augmentSparsity(SparsityPattern::Graph &sparsity, std::vector< dof_id_type > &n_nz, std::vector< dof_id_type > &n_oz) override
Will modify the sparsity pattern to add logical geometric connections.
Definition: LinearSystem.C:291
Adds contributions from face terms discretized using the finite volume method to the matrix and right...
elem_info_iterator ownedElemInfoEnd()
Definition: MooseMesh.C:1533
MOOSE now contains C++17 code, so give a reasonable error message stating what the user can do to add...
const ConsoleStream _console
An instance of helper class to write streams to the Console objects.
face_info_iterator ownedFaceInfoEnd()
Definition: MooseMesh.C:1516
MooseVariableFieldBase & getVariable(THREAD_ID tid, const std::string &var_name) const
Gets a reference to a variable of with specified name.
Definition: SystemBase.C:89
virtual ~LinearSystem()
virtual void initialSetup()
Setup Functions.
Definition: SystemBase.C:1525
LinearSystem(FEProblemBase &problem, const std::string &name)
Definition: LinearSystem.C:76
auto index_range(const T &sizable)
virtual System & system() override
Get the reference to the libMesh system.
Definition: LinearSystem.h:107
void computeLinearSystemTags(const std::set< TagID > &vector_tags, const std::set< TagID > &matrix_tags, const bool compute_gradients=true)
Compute the right hand side and the system matrix of the system for given tags.
Definition: LinearSystem.C:141
unsigned int THREAD_ID
Definition: MooseTypes.h:209