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 : #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"
21 : #include "MFEMMixedBilinearFormKernel.h"
22 : #include "ScaleIntegrator.h"
23 : #include "NLScaleIntegrator.h"
24 :
25 : namespace Moose::MFEM
26 : {
27 : class LinearSolverBase;
28 :
29 : /**
30 : * Class to store weak form components (bilinear and linear forms, and optionally
31 : * mixed and nonlinear forms) and build methods
32 : */
33 : class EquationSystem : public mfem::Operator
34 : {
35 :
36 : public:
37 1450 : EquationSystem() = default;
38 : ~EquationSystem() override;
39 :
40 : /// Add kernels.
41 : virtual void AddKernel(std::shared_ptr<MFEMKernel> kernel);
42 : virtual void AddIntegratedBC(std::shared_ptr<MFEMIntegratedBC> kernel);
43 : /// Add BC associated with essentially constrained DoFs on boundaries.
44 : virtual void AddEssentialBC(std::shared_ptr<MFEMEssentialBC> bc);
45 :
46 : /// Initialise
47 : virtual void Init(GridFunctions & gridfunctions,
48 : ComplexGridFunctions & cmplx_gridfunctions,
49 : mfem::AssemblyLevel assembly_level);
50 : /**
51 : * Assemble the linear part of the operator, assemble the right-hand side, apply essential and
52 : * eliminated-variable constraints, and populate the true-DoF vectors used by the solve.
53 : */
54 : void FormSystem(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS);
55 : /// Compute residual y = Mu
56 : void Mult(const mfem::Vector & u, mfem::Vector & residual) const override;
57 : /// Compute the contribution to the residual from nonlinear forms only.
58 : virtual void ComputeNonlinearResidual(const mfem::Vector & u, mfem::Vector & residual) const;
59 : /// Get Jacobian at the provided vector of true DoFs of trial variables
60 : mfem::Operator & GetGradient(const mfem::Vector & u) const override;
61 :
62 : /// Update variable from solution vector after solve
63 : virtual void SetTrialVariablesFromTrueVectors(const mfem::BlockVector & trueX) const;
64 :
65 : /// Set whether the nonlinear solver driving this equation system requires Jacobian information.
66 48 : void SetSolverRequiresGradient(bool requires_gradient)
67 : {
68 48 : _solver_requires_gradient = requires_gradient;
69 48 : }
70 :
71 : // Test variables are associated with linear forms,
72 : // whereas trial variables are associated with gridfunctions.
73 1466 : const std::vector<std::string> & GetTrialVarNames() const { return _trial_var_names; }
74 1492 : const std::vector<std::string> & GetTestVarNames() const { return _test_var_names; }
75 :
76 : /**
77 : * @returns Whether nonlinear integrators are present
78 : */
79 2135 : bool Nonlinear() const { return _non_linear; }
80 :
81 : /**
82 : * Prepare the provided linear solver. First calls SetupLOR on the solver if it's using a Low
83 : * Order Refined methodology and then calls SetOperator on the solver with the assembled linear
84 : * operator
85 : */
86 : void PrepareLinearSolver(LinearSolverBase & solver);
87 :
88 : protected:
89 : /// Add coupled variable to EquationSystem.
90 : virtual void AddCoupledVariableNameIfMissing(const std::string & coupled_var_name);
91 : /// Add eliminated variable to EquationSystem.
92 : virtual void AddEliminatedVariableNameIfMissing(const std::string & eliminated_var_name);
93 : /// Add test variable to EquationSystem.
94 : virtual void AddTestVariableNameIfMissing(const std::string & test_var_name);
95 : /// Set trial variable names from subset of coupled variables that have an associated test variable.
96 : virtual void SetTrialVariableNames();
97 :
98 : /// Deletes the HypreParMatrix associated with any pointer stored in _h_blocks,
99 : /// and then proceeds to delete all dynamically allocated memory for _h_blocks
100 : /// itself, resetting all dimensions to zero.
101 : void DeleteHBlocks();
102 :
103 : /// Deletes the HypreParMatrix associated with any pointer stored in _jacobian_blocks,
104 : /// and then proceeds to delete all dynamically allocated memory for _jacobian_blocks
105 : /// itself, resetting all dimensions to zero.
106 : void DeleteJacobianBlocks();
107 :
108 : bool VectorContainsName(const std::vector<std::string> & the_vector,
109 : const std::string & name) const;
110 :
111 : /// Apply essential BC(s) associated with var_name to set true DoFs of trial_gf and update
112 : /// markers of all essential boundaries
113 : virtual void ApplyEssentialBC(const std::string & var_name,
114 : mfem::ParGridFunction & trial_gf,
115 : mfem::Array<int> & global_ess_markers);
116 : /// Update all essentially constrained true DoF markers and values on boundaries
117 : virtual void ApplyEssentialBCs();
118 : /// Perform trivial eliminations of coupled variables lacking corresponding test variables
119 : virtual void EliminateCoupledVariables();
120 : /// Build linear forms and eliminate constrained DoFs
121 : virtual void BuildLinearForms();
122 : /// Build non-linear action forms
123 : virtual void BuildNonlinearForms();
124 : /// Build bilinear forms (diagonal Jacobian contributions)
125 : virtual void BuildBilinearForms();
126 : /// Build mixed bilinear forms (off-diagonal Jacobian contributions)
127 : virtual void BuildMixedBilinearForms();
128 : /// Build all forms comprising this EquationSystem
129 : virtual void BuildEquationSystem();
130 :
131 : /// Form linear components of system based on on- and off-diagonal bilinear form
132 : /// contributions, populate solution and RHS vectors of true DoFs, and apply constraints.
133 : virtual void FormLinearSystem(mfem::OperatorHandle & op,
134 : mfem::BlockVector & trueX,
135 : mfem::BlockVector & trueRHS);
136 : using mfem::Operator::FormSystemOperator;
137 : /// Form matrix-free representation of linear components of system operator.
138 : /// Used when EquationSystem assembly level is set to 'FULL', 'ELEMENT', 'PARTIAL', or 'NONE'.
139 : virtual void FormSystemOperator(mfem::OperatorHandle & op,
140 : mfem::BlockVector & trueX,
141 : mfem::BlockVector & trueRHS);
142 : /// Form matrix representation of linear components of system operator as a HypreParMatrix.
143 : /// Used when EquationSystem assembly level is set to 'LEGACY'.
144 : virtual void FormSystemMatrix(mfem::OperatorHandle & op,
145 : mfem::BlockVector & trueX,
146 : mfem::BlockVector & trueRHS);
147 : /// Compute Jacobian matrix at the provided vector of true DoFs of trial variables
148 : void FormJacobianMatrix(const mfem::Vector & u);
149 :
150 : /**
151 : * Template method for applying BilinearFormIntegrators on domains from kernels to a BilinearForm,
152 : * or MixedBilinearForm
153 : */
154 : template <class FormType>
155 : void ApplyDomainBLFIntegrators(
156 : const std::string & trial_var_name,
157 : const std::string & test_var_name,
158 : std::shared_ptr<FormType> form,
159 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map,
160 : std::optional<mfem::real_t> scale_factor = std::nullopt);
161 :
162 : /**
163 : * Apply domain LinearFormIntegrators from kernels to the linear form associated with the
164 : * supplied test variable.
165 : */
166 : void ApplyDomainLFIntegrators(
167 : const std::string & test_var_name,
168 : std::shared_ptr<mfem::ParLinearForm> form,
169 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map);
170 :
171 : /**
172 : * Apply domain NonlinearFormIntegrators from kernels to the nonlinear form associated with the
173 : * supplied test variable.
174 : */
175 : void ApplyDomainNLFIntegrators(
176 : const std::string & test_var_name,
177 : std::shared_ptr<mfem::ParNonlinearForm> form,
178 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map,
179 : std::optional<mfem::real_t> scale_factor = std::nullopt);
180 :
181 : /**
182 : * Template method for applying BilinearFormIntegrators on boundaries from integrated boundary
183 : * conditions to a BilinearForm, or MixedBilinearForm.
184 : */
185 : template <class FormType>
186 : void ApplyBoundaryBLFIntegrators(
187 : const std::string & trial_var_name,
188 : const std::string & test_var_name,
189 : std::shared_ptr<FormType> form,
190 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
191 : integrated_bc_map,
192 : std::optional<mfem::real_t> scale_factor = std::nullopt);
193 :
194 : /**
195 : * Apply boundary LinearFormIntegrators from integrated boundary conditions to the linear form
196 : * associated with the supplied test variable.
197 : */
198 : void ApplyBoundaryLFIntegrators(
199 : const std::string & test_var_name,
200 : std::shared_ptr<mfem::ParLinearForm> form,
201 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
202 : integrated_bc_map);
203 :
204 : /**
205 : * Apply boundary NonlinearFormIntegrators from integrated boundary conditions to the nonlinear
206 : * form associated with the supplied test variable.
207 : */
208 : void ApplyBoundaryNLFIntegrators(
209 : const std::string & test_var_name,
210 : std::shared_ptr<mfem::ParNonlinearForm> form,
211 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
212 : integrated_bc_map,
213 : std::optional<mfem::real_t> scale_factor = std::nullopt);
214 :
215 : /**
216 : * Whether this a complex equation system
217 : */
218 45 : virtual bool Complex() const { return false; }
219 :
220 : /// Names of all trial variables of kernels and boundary conditions
221 : /// added to this EquationSystem.
222 : std::vector<std::string> _coupled_var_names;
223 : /// Subset of _coupled_var_names of all variables corresponding to gridfunctions with degrees of
224 : /// freedom that comprise the state vector of this EquationSystem. This will differ from
225 : /// _coupled_var_names when time derivatives or other eliminated variables are present.
226 : std::vector<std::string> _trial_var_names;
227 : /// Names of all coupled variables without a corresponding test variable.
228 : std::vector<std::string> _eliminated_var_names;
229 : /// Pointers to coupled variables not part of the reduced EquationSystem.
230 : Moose::MFEM::GridFunctions _eliminated_variables;
231 : /// Names of all test variables corresponding to linear forms in this equation system
232 : std::vector<std::string> _test_var_names;
233 : /// Pointers to finite element spaces associated with test variables.
234 : std::vector<mfem::ParFiniteElementSpace *> _test_pfespaces;
235 : /// Pointers to finite element spaces associated with coupled variables.
236 : std::vector<mfem::ParFiniteElementSpace *> _coupled_pfespaces;
237 :
238 : // Components of weak form, named according to test variable
239 : NamedFieldsMap<mfem::ParBilinearForm> _blfs;
240 : NamedFieldsMap<mfem::ParLinearForm> _lfs;
241 : NamedFieldsMap<mfem::ParNonlinearForm> _nlfs;
242 : NamedFieldsMap<NamedFieldsMap<mfem::ParMixedBilinearForm>> _mblfs; // named according to trial var
243 :
244 : /// Gridfunctions holding essential constraints from Dirichlet BCs
245 : std::vector<std::unique_ptr<mfem::ParGridFunction>> _var_ess_constraints;
246 : std::vector<mfem::Array<int>> _ess_tdof_lists;
247 :
248 : mfem::Array2D<const mfem::HypreParMatrix *> _h_blocks, _jacobian_blocks;
249 : /// Arrays to store kernels to act on each component of weak form.
250 : /// Named according to test and trial variables.
251 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> _kernels_map;
252 : /// Arrays to store integrated BCs to act on each component of weak form.
253 : /// Named according to test and trial variables.
254 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> _integrated_bc_map;
255 : /// Arrays to store essential BCs to act on each component of weak form.
256 : /// Named according to test variable.
257 : NamedFieldsMap<std::vector<std::shared_ptr<MFEMEssentialBC>>> _essential_bc_map;
258 :
259 : // Operator handle for the jacobian
260 : mutable mfem::OperatorHandle _jacobian;
261 : // Operator handle for the linear components of the system operator
262 : mutable mfem::OperatorHandle _linear_operator;
263 : mfem::AssemblyLevel _assembly_level;
264 :
265 : // Pointer to GridFunctions to enable updates during nonlinear iterations
266 : Moose::MFEM::GridFunctions * _gfuncs;
267 : // Array storing block offsets of solution and residual vector
268 : mfem::Array<int> _block_true_offsets;
269 : // Boolean indicating if EquationSystem contains nonlinear integrators
270 : bool _non_linear = false;
271 : // Whether a nonlinear solver exists and whether it requires Jacobian/gradient information.
272 : bool _solver_requires_gradient = false;
273 :
274 : private:
275 : friend class EquationSystemProblemOperator;
276 : friend class ::MFEMProblemSolve;
277 : /// Disallowed inherited method
278 : using mfem::Operator::RecoverFEMSolution;
279 : };
280 :
281 : template <class FormType>
282 : void
283 5076 : EquationSystem::ApplyDomainBLFIntegrators(
284 : const std::string & trial_var_name,
285 : const std::string & test_var_name,
286 : std::shared_ptr<FormType> form,
287 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map,
288 : std::optional<mfem::real_t> scale_factor)
289 : {
290 5076 : if (kernels_map.Has(test_var_name) && kernels_map.Get(test_var_name)->Has(trial_var_name))
291 : {
292 4925 : auto kernels = kernels_map.GetRef(test_var_name).GetRef(trial_var_name);
293 10664 : for (auto & kernel : kernels)
294 : {
295 5739 : mfem::BilinearFormIntegrator * integ = kernel->createBFIntegrator();
296 :
297 5739 : if (integ)
298 : {
299 4967 : if (scale_factor.has_value())
300 1300 : integ = new ScaleIntegrator(integ, scale_factor.value(), true);
301 4967 : kernel->isSubdomainRestricted()
302 4967 : ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
303 4877 : : form->AddDomainIntegrator(std::move(integ));
304 : }
305 : }
306 4925 : }
307 5076 : }
308 :
309 : template <class FormType>
310 : void
311 2200 : EquationSystem::ApplyBoundaryBLFIntegrators(
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 2332 : if (integrated_bc_map.Has(test_var_name) &&
320 132 : integrated_bc_map.Get(test_var_name)->Has(trial_var_name))
321 : {
322 110 : auto bcs = integrated_bc_map.GetRef(test_var_name).GetRef(trial_var_name);
323 238 : for (auto & bc : bcs)
324 : {
325 128 : mfem::BilinearFormIntegrator * integ = bc->createBFIntegrator();
326 :
327 128 : if (integ)
328 : {
329 36 : if (scale_factor.has_value())
330 36 : integ = new ScaleIntegrator(integ, scale_factor.value(), true);
331 36 : bc->isBoundaryRestricted()
332 36 : ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
333 0 : : form->AddBoundaryIntegrator(std::move(integ));
334 : }
335 : }
336 110 : }
337 2200 : }
338 :
339 : } // namespace Moose::MFEM
340 :
341 : #endif
|
