https://mooseframework.inl.gov
EquationSystem.h
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 #ifdef MOOSE_MFEM_ENABLED
11 
12 #pragma once
13 
14 #include "libmesh/ignore_warnings.h"
15 #include "mfem/miniapps/common/pfem_extras.hpp"
16 #include "libmesh/restore_warnings.h"
17 #include "MFEMIntegratedBC.h"
18 #include "MFEMEssentialBC.h"
19 #include "MFEMContainers.h"
20 #include "MFEMKernel.h"
22 #include "ScaleIntegrator.h"
23 #include "NLScaleIntegrator.h"
24 
25 namespace Moose::MFEM
26 {
27 
32 class EquationSystem : public mfem::Operator
33 {
34 
35 public:
36  EquationSystem() = default;
37  ~EquationSystem() override;
38 
40  virtual void AddKernel(std::shared_ptr<MFEMKernel> kernel);
41  virtual void AddIntegratedBC(std::shared_ptr<MFEMIntegratedBC> kernel);
43  virtual void AddEssentialBC(std::shared_ptr<MFEMEssentialBC> bc);
44 
46  virtual void Init(GridFunctions & gridfunctions,
47  ComplexGridFunctions & cmplx_gridfunctions,
48  mfem::AssemblyLevel assembly_level);
53  void FormSystem(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS);
55  void Mult(const mfem::Vector & u, mfem::Vector & residual) const override;
57  mfem::Operator & GetGradient(const mfem::Vector & u) const override;
58 
60  virtual void SetTrialVariablesFromTrueVectors(const mfem::BlockVector & trueX) const;
61 
62  // Test variables are associated with linear forms,
63  // whereas trial variables are associated with gridfunctions.
64  const std::vector<std::string> & GetTrialVarNames() const { return _trial_var_names; }
65  const std::vector<std::string> & GetTestVarNames() const { return _test_var_names; }
66 
67 protected:
69  virtual void AddCoupledVariableNameIfMissing(const std::string & coupled_var_name);
71  virtual void AddEliminatedVariableNameIfMissing(const std::string & eliminated_var_name);
73  virtual void AddTestVariableNameIfMissing(const std::string & test_var_name);
75  virtual void SetTrialVariableNames();
76 
80  void DeleteHBlocks();
81 
85  void DeleteJacobianBlocks();
86 
87  bool VectorContainsName(const std::vector<std::string> & the_vector,
88  const std::string & name) const;
89 
92  virtual void ApplyEssentialBC(const std::string & var_name,
93  mfem::ParGridFunction & trial_gf,
94  mfem::Array<int> & global_ess_markers);
96  virtual void ApplyEssentialBCs();
98  virtual void EliminateCoupledVariables();
100  virtual void BuildLinearForms();
102  virtual void BuildNonlinearForms();
104  virtual void BuildBilinearForms();
106  virtual void BuildMixedBilinearForms();
108  virtual void BuildEquationSystem();
109 
112  virtual void FormLinearSystem(mfem::OperatorHandle & op,
113  mfem::BlockVector & trueX,
114  mfem::BlockVector & trueRHS);
117  virtual void FormSystemOperator(mfem::OperatorHandle & op,
118  mfem::BlockVector & trueX,
119  mfem::BlockVector & trueRHS);
122  virtual void FormSystemMatrix(mfem::OperatorHandle & op,
123  mfem::BlockVector & trueX,
124  mfem::BlockVector & trueRHS);
126  void FormJacobianMatrix(const mfem::Vector & u);
127 
132  template <class FormType>
134  const std::string & trial_var_name,
135  const std::string & test_var_name,
136  std::shared_ptr<FormType> form,
137  NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map,
138  std::optional<mfem::real_t> scale_factor = std::nullopt);
139 
141  const std::string & test_var_name,
142  std::shared_ptr<mfem::ParLinearForm> form,
143  NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map);
144 
146  const std::string & test_var_name,
147  std::shared_ptr<mfem::ParNonlinearForm> form,
148  NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map,
149  std::optional<mfem::real_t> scale_factor = std::nullopt);
150 
151  template <class FormType>
153  const std::string & trial_var_name,
154  const std::string & test_var_name,
155  std::shared_ptr<FormType> form,
156  NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
157  integrated_bc_map,
158  std::optional<mfem::real_t> scale_factor = std::nullopt);
159 
161  const std::string & test_var_name,
162  std::shared_ptr<mfem::ParLinearForm> form,
163  NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
164  integrated_bc_map);
165 
167  const std::string & test_var_name,
168  std::shared_ptr<mfem::ParNonlinearForm> form,
169  NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
170  integrated_bc_map,
171  std::optional<mfem::real_t> scale_factor = std::nullopt);
172 
175  std::vector<std::string> _coupled_var_names;
179  std::vector<std::string> _trial_var_names;
181  std::vector<std::string> _eliminated_var_names;
185  std::vector<std::string> _test_var_names;
187  std::vector<mfem::ParFiniteElementSpace *> _test_pfespaces;
189  std::vector<mfem::ParFiniteElementSpace *> _coupled_pfespaces;
190 
191  // Components of weak form, named according to test variable
196 
198  std::vector<std::unique_ptr<mfem::ParGridFunction>> _var_ess_constraints;
199  std::vector<mfem::Array<int>> _ess_tdof_lists;
200 
201  mfem::Array2D<const mfem::HypreParMatrix *> _h_blocks, _jacobian_blocks;
211 
212  // Operator handle for the jacobian
213  mutable mfem::OperatorHandle _jacobian;
214  // Operator handle for the linear components of the system operator
215  mutable mfem::OperatorHandle _linear_operator;
216  mfem::AssemblyLevel _assembly_level;
217 
218  // Pointer to GridFunctions to enable updates during nonlinear iterations
220  // Array storing block offsets of solution and residual vector
221  mfem::Array<int> _block_true_offsets;
222  // Boolean indicating if EquationSystem contains nonlinear integrators
223  bool _non_linear = false;
224 
225 private:
227  friend class ::MFEMProblemSolve;
229  using mfem::Operator::RecoverFEMSolution;
230 };
231 
232 template <class FormType>
233 void
235  const std::string & trial_var_name,
236  const std::string & test_var_name,
237  std::shared_ptr<FormType> form,
238  NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map,
239  std::optional<mfem::real_t> scale_factor)
240 {
241  if (kernels_map.Has(test_var_name) && kernels_map.Get(test_var_name)->Has(trial_var_name))
242  {
243  auto kernels = kernels_map.GetRef(test_var_name).GetRef(trial_var_name);
244  for (auto & kernel : kernels)
245  {
246  mfem::BilinearFormIntegrator * integ = kernel->createBFIntegrator();
247 
248  if (integ)
249  {
250  if (scale_factor.has_value())
251  integ = new ScaleIntegrator(integ, scale_factor.value(), true);
252  kernel->isSubdomainRestricted()
253  ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
254  : form->AddDomainIntegrator(std::move(integ));
255  }
256  }
257  }
258 }
259 
260 inline void
262  const std::string & test_var_name,
263  std::shared_ptr<mfem::ParLinearForm> form,
264  NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map)
265 {
266  if (kernels_map.Has(test_var_name) && kernels_map.Get(test_var_name)->Has(test_var_name))
267  {
268  auto kernels = kernels_map.GetRef(test_var_name).GetRef(test_var_name);
269  for (auto & kernel : kernels)
270  {
271  mfem::LinearFormIntegrator * integ = kernel->createLFIntegrator();
272 
273  if (integ)
274  {
275  kernel->isSubdomainRestricted()
276  ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
277  : form->AddDomainIntegrator(std::move(integ));
278  }
279  }
280  }
281 }
282 
283 inline void
285  const std::string & test_var_name,
286  std::shared_ptr<mfem::ParNonlinearForm> form,
287  NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map,
288  std::optional<mfem::real_t> scale_factor)
289 {
290  if (kernels_map.Has(test_var_name) && kernels_map.Get(test_var_name)->Has(test_var_name))
291  {
292  auto kernels = kernels_map.GetRef(test_var_name).GetRef(test_var_name);
293  for (auto & kernel : kernels)
294  {
295  mfem::NonlinearFormIntegrator * integ = kernel->createNLIntegrator();
296  if (integ)
297  {
298  _non_linear = true;
299  if (scale_factor.has_value())
300  integ = new NLScaleIntegrator(integ, scale_factor.value(), true);
301  kernel->isSubdomainRestricted()
302  ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
303  : form->AddDomainIntegrator(std::move(integ));
304  }
305  }
306  }
307 }
308 
309 template <class FormType>
310 void
312  const std::string & trial_var_name,
313  const std::string & test_var_name,
314  std::shared_ptr<FormType> form,
315  NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
316  integrated_bc_map,
317  std::optional<mfem::real_t> scale_factor)
318 {
319  if (integrated_bc_map.Has(test_var_name) &&
320  integrated_bc_map.Get(test_var_name)->Has(trial_var_name))
321  {
322  auto bcs = integrated_bc_map.GetRef(test_var_name).GetRef(trial_var_name);
323  for (auto & bc : bcs)
324  {
325  mfem::BilinearFormIntegrator * integ = bc->createBFIntegrator();
326 
327  if (integ)
328  {
329  if (scale_factor.has_value())
330  integ = new ScaleIntegrator(integ, scale_factor.value(), true);
331  bc->isBoundaryRestricted()
332  ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
333  : form->AddBoundaryIntegrator(std::move(integ));
334  }
335  }
336  }
337 }
338 
339 inline void
341  const std::string & test_var_name,
342  std::shared_ptr<mfem::ParLinearForm> form,
343  NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
344  integrated_bc_map)
345 {
346  if (integrated_bc_map.Has(test_var_name) &&
347  integrated_bc_map.Get(test_var_name)->Has(test_var_name))
348  {
349  auto bcs = integrated_bc_map.GetRef(test_var_name).GetRef(test_var_name);
350  for (auto & bc : bcs)
351  {
352  mfem::LinearFormIntegrator * integ = bc->createLFIntegrator();
353 
354  if (integ)
355  {
356  bc->isBoundaryRestricted()
357  ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
358  : form->AddBoundaryIntegrator(std::move(integ));
359  }
360  }
361  }
362 }
363 
364 inline void
366  const std::string & test_var_name,
367  std::shared_ptr<mfem::ParNonlinearForm> form,
368  NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
369  integrated_bc_map,
370  std::optional<mfem::real_t> scale_factor)
371 {
372  if (integrated_bc_map.Has(test_var_name) &&
373  integrated_bc_map.Get(test_var_name)->Has(test_var_name))
374  {
375  auto bcs = integrated_bc_map.GetRef(test_var_name).GetRef(test_var_name);
376  for (auto & bc : bcs)
377  {
378  mfem::NonlinearFormIntegrator * integ = bc->createNLIntegrator();
379  if (integ)
380  {
381  _non_linear = true;
382  if (scale_factor.has_value())
383  integ = new NLScaleIntegrator(integ, scale_factor.value(), true);
384  bc->isBoundaryRestricted()
385  ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
386  : form->AddBoundaryIntegrator(std::move(integ));
387  }
388  }
389  }
390 }
391 
392 } // namespace Moose::MFEM
393 
394 #endif
void ApplyDomainNLFIntegrators(const std::string &test_var_name, std::shared_ptr< mfem::ParNonlinearForm > form, NamedFieldsMap< NamedFieldsMap< std::vector< std::shared_ptr< MFEMKernel >>>> &kernels_map, std::optional< mfem::real_t > scale_factor=std::nullopt)
NamedFieldsMap< mfem::ParBilinearForm > _blfs
virtual void EliminateCoupledVariables()
Perform trivial eliminations of coupled variables lacking corresponding test variables.
NamedFieldsMap< NamedFieldsMap< mfem::ParMixedBilinearForm > > _mblfs
virtual void AddTestVariableNameIfMissing(const std::string &test_var_name)
Add test variable to EquationSystem.
virtual void AddKernel(std::shared_ptr< MFEMKernel > kernel)
Add kernels.
virtual void BuildBilinearForms()
Build bilinear forms (diagonal Jacobian contributions)
std::vector< mfem::ParFiniteElementSpace * > _coupled_pfespaces
Pointers to finite element spaces associated with coupled variables.
virtual void SetTrialVariablesFromTrueVectors(const mfem::BlockVector &trueX) const
Update variable from solution vector after solve.
void DeleteJacobianBlocks()
Deletes the HypreParMatrix associated with any pointer stored in _jacobian_blocks, and then proceeds to delete all dynamically allocated memory for _jacobian_blocks itself, resetting all dimensions to zero.
bool VectorContainsName(const std::vector< std::string > &the_vector, const std::string &name) const
NamedFieldsMap< NamedFieldsMap< std::vector< std::shared_ptr< MFEMIntegratedBC > > > > _integrated_bc_map
Arrays to store integrated BCs to act on each component of weak form.
std::vector< std::string > _eliminated_var_names
Names of all coupled variables without a corresponding test variable.
void FormJacobianMatrix(const mfem::Vector &u)
Compute Jacobian matrix at the provided vector of true DoFs of trial variables.
NamedFieldsMap< mfem::ParNonlinearForm > _nlfs
Class to store weak form components (bilinear and linear forms, and optionally mixed and nonlinear fo...
Moose::MFEM::GridFunctions _eliminated_variables
Pointers to coupled variables not part of the reduced EquationSystem.
virtual void AddEssentialBC(std::shared_ptr< MFEMEssentialBC > bc)
Add BC associated with essentially constrained DoFs on boundaries.
NonlinearFormIntegrator which scales its results by a constant value.
std::vector< mfem::Array< int > > _ess_tdof_lists
mfem::AssemblyLevel _assembly_level
void ApplyBoundaryBLFIntegrators(const std::string &trial_var_name, const std::string &test_var_name, std::shared_ptr< FormType > form, NamedFieldsMap< NamedFieldsMap< std::vector< std::shared_ptr< MFEMIntegratedBC >>>> &integrated_bc_map, std::optional< mfem::real_t > scale_factor=std::nullopt)
std::vector< mfem::ParFiniteElementSpace * > _test_pfespaces
Pointers to finite element spaces associated with test variables.
void FormSystem(mfem::BlockVector &trueX, mfem::BlockVector &trueRHS)
Assemble the linear part of the operator, assemble the right-hand side, apply essential and eliminate...
void Mult(const mfem::Vector &u, mfem::Vector &residual) const override
Compute residual y = Mu.
mfem::Array< int > _block_true_offsets
std::vector< std::string > _trial_var_names
Subset of _coupled_var_names of all variables corresponding to gridfunctions with degrees of freedom ...
NamedFieldsMap< std::vector< std::shared_ptr< MFEMEssentialBC > > > _essential_bc_map
Arrays to store essential BCs to act on each component of weak form.
void ApplyDomainBLFIntegrators(const std::string &trial_var_name, const std::string &test_var_name, std::shared_ptr< FormType > form, NamedFieldsMap< NamedFieldsMap< std::vector< std::shared_ptr< MFEMKernel >>>> &kernels_map, std::optional< mfem::real_t > scale_factor=std::nullopt)
Template method for applying BilinearFormIntegrators on domains from kernels to a BilinearForm...
mfem::OperatorHandle _jacobian
NamedFieldsMap< NamedFieldsMap< std::vector< std::shared_ptr< MFEMKernel > > > > _kernels_map
Arrays to store kernels to act on each component of weak form.
std::vector< std::string > _test_var_names
Names of all test variables corresponding to linear forms in this equation system.
Steady-state problem operator with an equation system.
virtual void BuildEquationSystem()
Build all forms comprising this EquationSystem.
std::vector< std::unique_ptr< mfem::ParGridFunction > > _var_ess_constraints
Gridfunctions holding essential constraints from Dirichlet BCs.
void ApplyBoundaryLFIntegrators(const std::string &test_var_name, std::shared_ptr< mfem::ParLinearForm > form, NamedFieldsMap< NamedFieldsMap< std::vector< std::shared_ptr< MFEMIntegratedBC >>>> &integrated_bc_map)
virtual void BuildMixedBilinearForms()
Build mixed bilinear forms (off-diagonal Jacobian contributions)
virtual void FormSystemOperator(mfem::OperatorHandle &op, mfem::BlockVector &trueX, mfem::BlockVector &trueRHS)
Form matrix-free representation of linear components of system operator.
Integrator which scales its results by a constant value.
mfem::Array2D< const mfem::HypreParMatrix * > _h_blocks
virtual void SetTrialVariableNames()
Set trial variable names from subset of coupled variables that have an associated test variable...
virtual void FormSystemMatrix(mfem::OperatorHandle &op, mfem::BlockVector &trueX, mfem::BlockVector &trueRHS)
Form matrix representation of linear components of system operator as a HypreParMatrix.
virtual void AddCoupledVariableNameIfMissing(const std::string &coupled_var_name)
Add coupled variable to EquationSystem.
virtual void FormLinearSystem(mfem::OperatorHandle &op, mfem::BlockVector &trueX, mfem::BlockVector &trueRHS)
Form linear components of system based on on- and off-diagonal bilinear form contributions, populate solution and RHS vectors of true DoFs, and apply constraints.
void DeleteHBlocks()
Deletes the HypreParMatrix associated with any pointer stored in _h_blocks, and then proceeds to dele...
std::vector< std::string > _coupled_var_names
Names of all trial variables of kernels and boundary conditions added to this EquationSystem.
virtual void BuildNonlinearForms()
Build non-linear action forms.
NamedFieldsMap< mfem::ParLinearForm > _lfs
virtual void BuildLinearForms()
Build linear forms and eliminate constrained DoFs.
const std::vector< std::string > & GetTrialVarNames() const
const std::vector< std::string > & GetTestVarNames() const
void ApplyBoundaryNLFIntegrators(const std::string &test_var_name, std::shared_ptr< mfem::ParNonlinearForm > form, NamedFieldsMap< NamedFieldsMap< std::vector< std::shared_ptr< MFEMIntegratedBC >>>> &integrated_bc_map, std::optional< mfem::real_t > scale_factor=std::nullopt)
virtual void ApplyEssentialBCs()
Update all essentially constrained true DoF markers and values on boundaries.
void ApplyDomainLFIntegrators(const std::string &test_var_name, std::shared_ptr< mfem::ParLinearForm > form, NamedFieldsMap< NamedFieldsMap< std::vector< std::shared_ptr< MFEMKernel >>>> &kernels_map)
Utilities for converting between vector(s) of libMesh Points and MFEM Vector(s).
mfem::OperatorHandle _linear_operator
virtual void AddEliminatedVariableNameIfMissing(const std::string &eliminated_var_name)
Add eliminated variable to EquationSystem.
virtual void Init(GridFunctions &gridfunctions, ComplexGridFunctions &cmplx_gridfunctions, mfem::AssemblyLevel assembly_level)
Initialise.
mfem::Operator & GetGradient(const mfem::Vector &u) const override
Get Jacobian at the provided vector of true DoFs of trial variables.
virtual void ApplyEssentialBC(const std::string &var_name, mfem::ParGridFunction &trial_gf, mfem::Array< int > &global_ess_markers)
Apply essential BC(s) associated with var_name to set true DoFs of trial_gf and update markers of all...
virtual void AddIntegratedBC(std::shared_ptr< MFEMIntegratedBC > kernel)
mfem::Array2D< const mfem::HypreParMatrix * > _jacobian_blocks
Moose::MFEM::GridFunctions * _gfuncs