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 CoefficientManager;
28 :
29 : /**
30 : * Owns the weak-form mathematics of a MOOSE MFEM problem.
31 : *
32 : * An EquationSystem stores weak-form components (bilinear, linear, mixed-bilinear,
33 : * and nonlinear forms) contributed by kernels and boundary conditions. It forms the
34 : * constrained linear part, keeps nonlinear action forms for residual/Jacobian
35 : * evaluation, and exposes the solve interface as an mfem::Operator. It is responsible
36 : * for applying essential-DoF constraints and propagating either itself or the assembled
37 : * linear operator (and bilinear form, for LOR) to the configured solver tree.
38 : *
39 : * EquationSystem is *not* responsible for grid-function bookkeeping, time stepping, or
40 : * solver selection - those belong to the ProblemOperator layer.
41 : *
42 : * @see ProblemOperatorBase for the conceptual split between the two layers.
43 : */
44 : class EquationSystem : public mfem::Operator
45 : {
46 :
47 : public:
48 1705 : EquationSystem() = default;
49 : ~EquationSystem() override;
50 :
51 : /// Add kernels.
52 : virtual void AddKernel(std::shared_ptr<MFEMKernel> kernel);
53 : virtual void AddIntegratedBC(std::shared_ptr<MFEMIntegratedBC> kernel);
54 : /// Add BC associated with essentially constrained DoFs on boundaries.
55 : virtual void AddEssentialBC(std::shared_ptr<MFEMEssentialBC> bc);
56 :
57 : /// Initialise
58 : virtual void Init(GridFunctions & gridfunctions,
59 : ComplexGridFunctions & cmplx_gridfunctions,
60 : mfem::AssemblyLevel assembly_level);
61 : /**
62 : * Build all weak-form components via BuildEquationSystem(), form the constrained linear part of
63 : * the system, and populate the true-DoF vectors used by the solve.
64 : *
65 : * For nonlinear problems, nonlinear forms are registered here but are not folded into the
66 : * assembled linear operator. Their residual action is evaluated by Mult(), and any Jacobian
67 : * contribution is formed at the current nonlinear iterate by GetGradient().
68 : *
69 : * This is the single public entry point for callers that want the equation system prepared for a
70 : * solve. Subclasses customize the form-building step by overriding BuildEquationSystem().
71 : */
72 : void FormSystem(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS);
73 : /// Compute residual y = Mu
74 : void Mult(const mfem::Vector & u, mfem::Vector & residual) const override;
75 : /// Compute the contribution to the residual from nonlinear forms only.
76 : virtual void ComputeNonlinearResidual(const mfem::Vector & u, mfem::Vector & residual) const;
77 : /// Get Jacobian at the provided vector of true DoFs of trial variables
78 : mfem::Operator & GetGradient(const mfem::Vector & u) const override;
79 : /// Get operator handle for linear component of system operator
80 26 : mfem::OperatorHandle & GetLinearOperator() const { return _linear_operator; };
81 :
82 : /// Update variable from solution vector after solve
83 : virtual void SetTrialVariablesFromTrueVectors(const mfem::BlockVector & trueX) const;
84 :
85 : /// Set whether an external object (such as a nonlinear solver) requires Jacobian information for this EquationSystem.
86 62 : void SetGradientRequired(bool requires_gradient) { _gradient_required = requires_gradient; }
87 :
88 : /// Set the coefficient manager to notify when trial variables are updated, so that stored
89 : /// projections of solution-dependent coefficients are invalidated.
90 1635 : void SetCoefficientManager(CoefficientManager & coefficients)
91 : {
92 1635 : _coefficient_manager = &coefficients;
93 1635 : }
94 :
95 : // Test variables are associated with linear forms,
96 : // whereas trial variables are associated with gridfunctions.
97 1661 : const std::vector<std::string> & GetTrialVarNames() const { return _trial_var_names; }
98 1751 : const std::vector<std::string> & GetTestVarNames() const { return _test_var_names; }
99 :
100 : /// Getter for block true offsets associated with the EquationSystem operator
101 : const mfem::Array<int> & GetBlockOffsets() const { return _block_true_offsets; }
102 :
103 : /**
104 : * @returns a reference to the MFEM ParBilinearForm corresponding to test_var_name
105 : */
106 64 : mfem::ParBilinearForm & GetBilinearForm(const std::string & test_var_name)
107 : {
108 64 : return _blfs.GetRef(test_var_name);
109 : }
110 :
111 : /**
112 : * @returns a reference to the MFEM ParGridFunction corresponding to trial_var_name
113 : */
114 : mfem::ParGridFunction & GetGridFunction(const std::string & trial_var_name)
115 : {
116 : return _gfuncs->GetRef(trial_var_name);
117 : }
118 :
119 : /**
120 : * @returns Whether nonlinear integrators are present
121 : */
122 : bool Nonlinear() const { return _non_linear; }
123 :
124 : /// The true-DoF vector used for the most recent Jacobian linearization.
125 : const mfem::Vector & GetLinearizationPoint() const;
126 :
127 : /**
128 : * Build a fresh ParBilinearForm on the given FESpace using the same kernels as the main
129 : * system's bilinear form for var_name. Caller owns the returned form.
130 : */
131 : std::shared_ptr<mfem::ParBilinearForm>
132 : BuildBilinearFormForFESpace(const std::string & var_name,
133 : mfem::ParFiniteElementSpace & fespace,
134 : mfem::AssemblyLevel assembly_level);
135 :
136 : /**
137 : * Build a fresh ParNonlinearForm on the given FESpace using the same kernels as the main
138 : * system's nonlinear form for var_name. Caller owns the returned form.
139 : */
140 : std::shared_ptr<mfem::ParNonlinearForm>
141 : BuildNonlinearFormForFESpace(const std::string & var_name,
142 : mfem::ParFiniteElementSpace & fespace,
143 : mfem::AssemblyLevel assembly_level);
144 :
145 : /**
146 : * Return the essential boundary attribute marker array for a given trial variable.
147 : * The returned array has size == pmesh.bdr_attributes.Max() with 1 at essential boundaries.
148 : */
149 : mfem::Array<int> & GetEssentialBoundaryMarkers(const std::string & var_name);
150 :
151 : /// @returns Whether this EquationSystem includes complex components
152 73 : virtual bool IsComplex() const { return false; }
153 : /// @returns Whether this EquationSystem represents an eigenproblem
154 11 : virtual bool IsEigen() const { return false; }
155 : /// @returns Whether this EquationSystem has time-dependent components
156 0 : virtual bool IsTimeDependent() const { return false; }
157 : /// @returns Whether this is a multivariate (maybe mixed) equation system
158 71 : bool IsMultivariate() const { return _test_var_names.size() > 1; }
159 : /// @returns Whether nonlinear integrators are present in the equation system
160 9 : bool IsNonlinear() const { return _non_linear; }
161 :
162 : /// Build all forms comprising this EquationSystem
163 : virtual void BuildEquationSystem();
164 :
165 : protected:
166 : /// Add coupled variable to EquationSystem.
167 : virtual void AddCoupledVariableNameIfMissing(const std::string & coupled_var_name);
168 : /// Add eliminated variable to EquationSystem.
169 : virtual void AddEliminatedVariableNameIfMissing(const std::string & eliminated_var_name);
170 : /// Add test variable to EquationSystem.
171 : virtual void AddTestVariableNameIfMissing(const std::string & test_var_name);
172 : /// Set trial variable names from subset of coupled variables that have an associated test variable.
173 : virtual void SetTrialVariableNames();
174 :
175 : /// Deletes the HypreParMatrix associated with any pointer stored in _h_blocks,
176 : /// and then proceeds to delete all dynamically allocated memory for _h_blocks
177 : /// itself, resetting all dimensions to zero.
178 : void DeleteHBlocks();
179 :
180 : /// Deletes the HypreParMatrix associated with any pointer stored in _jacobian_blocks,
181 : /// and then proceeds to delete all dynamically allocated memory for _jacobian_blocks
182 : /// itself, resetting all dimensions to zero.
183 : void DeleteJacobianBlocks();
184 :
185 : bool VectorContainsName(const std::vector<std::string> & the_vector,
186 : const std::string & name) const;
187 :
188 : /// Apply essential BC(s) associated with var_name to set true DoFs of trial_gf and update
189 : /// markers of all essential boundaries
190 : virtual void ApplyEssentialBC(const std::string & var_name,
191 : mfem::ParGridFunction & trial_gf,
192 : mfem::Array<int> & global_ess_markers);
193 : /// Update all essentially constrained true DoF markers and values on boundaries
194 : virtual void ApplyEssentialBCs();
195 : /// Perform trivial eliminations of coupled variables lacking corresponding test variables
196 : virtual void EliminateCoupledVariables();
197 : /// Build linear forms and eliminate constrained DoFs
198 : virtual void BuildLinearForms();
199 : /// Build non-linear action forms
200 : virtual void BuildNonlinearForms();
201 : /// Build bilinear forms (diagonal Jacobian contributions)
202 : virtual void BuildBilinearForms();
203 : /// Build mixed bilinear forms (off-diagonal Jacobian contributions)
204 : virtual void BuildMixedBilinearForms();
205 :
206 : /// Form linear components of system based on on- and off-diagonal bilinear form
207 : /// contributions, populate solution and RHS vectors of true DoFs, and apply constraints.
208 : virtual void FormLinearSystem(mfem::OperatorHandle & op,
209 : mfem::BlockVector & trueX,
210 : mfem::BlockVector & trueRHS);
211 : using mfem::Operator::FormSystemOperator;
212 : /// Form matrix-free representation of linear components of system operator.
213 : /// Used when EquationSystem assembly level is set to 'FULL', 'ELEMENT', 'PARTIAL', or 'NONE'.
214 : virtual void FormSystemOperator(mfem::OperatorHandle & op,
215 : mfem::BlockVector & trueX,
216 : mfem::BlockVector & trueRHS);
217 : /// Form matrix representation of linear components of system operator as a HypreParMatrix.
218 : /// Used when EquationSystem assembly level is set to 'LEGACY'.
219 : virtual void FormSystemMatrix(mfem::OperatorHandle & op,
220 : mfem::BlockVector & trueX,
221 : mfem::BlockVector & trueRHS);
222 : /// Compute Jacobian matrix at the provided vector of true DoFs of trial variables
223 : void FormJacobianMatrix(const mfem::Vector & u);
224 :
225 : /**
226 : * Template method for applying BilinearFormIntegrators on domains from kernels to a BilinearForm,
227 : * or MixedBilinearForm
228 : */
229 : template <class FormType>
230 : void ApplyDomainBLFIntegrators(
231 : const std::string & trial_var_name,
232 : const std::string & test_var_name,
233 : std::shared_ptr<FormType> form,
234 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map,
235 : std::optional<mfem::real_t> scale_factor = std::nullopt);
236 :
237 : /**
238 : * Apply domain LinearFormIntegrators from kernels to the linear form associated with the
239 : * supplied test variable.
240 : */
241 : void ApplyDomainLFIntegrators(
242 : const std::string & test_var_name,
243 : std::shared_ptr<mfem::ParLinearForm> form,
244 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map);
245 :
246 : /**
247 : * Apply domain NonlinearFormIntegrators from kernels to the nonlinear form associated with the
248 : * supplied test variable.
249 : */
250 : void ApplyDomainNLFIntegrators(
251 : const std::string & test_var_name,
252 : std::shared_ptr<mfem::ParNonlinearForm> form,
253 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map,
254 : std::optional<mfem::real_t> scale_factor = std::nullopt);
255 :
256 : /**
257 : * Template method for applying BilinearFormIntegrators on boundaries from integrated boundary
258 : * conditions to a BilinearForm, or MixedBilinearForm.
259 : */
260 : template <class FormType>
261 : void ApplyBoundaryBLFIntegrators(
262 : const std::string & trial_var_name,
263 : const std::string & test_var_name,
264 : std::shared_ptr<FormType> form,
265 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
266 : integrated_bc_map,
267 : std::optional<mfem::real_t> scale_factor = std::nullopt);
268 :
269 : /**
270 : * Apply boundary LinearFormIntegrators from integrated boundary conditions to the linear form
271 : * associated with the supplied test variable.
272 : */
273 : void ApplyBoundaryLFIntegrators(
274 : const std::string & test_var_name,
275 : std::shared_ptr<mfem::ParLinearForm> form,
276 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
277 : integrated_bc_map);
278 :
279 : /**
280 : * Apply boundary NonlinearFormIntegrators from integrated boundary conditions to the nonlinear
281 : * form associated with the supplied test variable.
282 : */
283 : void ApplyBoundaryNLFIntegrators(
284 : const std::string & test_var_name,
285 : std::shared_ptr<mfem::ParNonlinearForm> form,
286 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
287 : integrated_bc_map,
288 : std::optional<mfem::real_t> scale_factor = std::nullopt);
289 :
290 : /// Names of all trial variables of kernels and boundary conditions
291 : /// added to this EquationSystem.
292 : std::vector<std::string> _coupled_var_names;
293 : /// Subset of _coupled_var_names of all variables corresponding to gridfunctions with degrees of
294 : /// freedom that comprise the state vector of this EquationSystem. This will differ from
295 : /// _coupled_var_names when time derivatives or other eliminated variables are present.
296 : std::vector<std::string> _trial_var_names;
297 : /// Names of all coupled variables without a corresponding test variable.
298 : std::vector<std::string> _eliminated_var_names;
299 : /// Pointers to coupled variables not part of the reduced EquationSystem.
300 : Moose::MFEM::GridFunctions _eliminated_variables;
301 : /// Names of all test variables corresponding to linear forms in this equation system
302 : std::vector<std::string> _test_var_names;
303 : /// Pointers to finite element spaces associated with test variables.
304 : std::vector<mfem::ParFiniteElementSpace *> _test_pfespaces;
305 : /// Pointers to finite element spaces associated with coupled variables.
306 : std::vector<mfem::ParFiniteElementSpace *> _coupled_pfespaces;
307 :
308 : // Components of weak form, named according to test variable
309 : NamedFieldsMap<mfem::ParBilinearForm> _blfs;
310 : NamedFieldsMap<mfem::ParLinearForm> _lfs;
311 : NamedFieldsMap<mfem::ParNonlinearForm> _nlfs;
312 : NamedFieldsMap<NamedFieldsMap<mfem::ParMixedBilinearForm>> _mblfs; // named according to trial var
313 :
314 : /// Gridfunctions holding essential constraints from Dirichlet BCs
315 : std::vector<std::unique_ptr<mfem::ParGridFunction>> _var_ess_constraints;
316 : std::vector<mfem::Array<int>> _ess_tdof_lists;
317 : std::vector<mfem::Array<int>> _ess_markers;
318 :
319 : mfem::Array2D<const mfem::HypreParMatrix *> _h_blocks, _jacobian_blocks;
320 : /// Arrays to store kernels to act on each component of weak form.
321 : /// Named according to test and trial variables.
322 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> _kernels_map;
323 : /// Arrays to store integrated BCs to act on each component of weak form.
324 : /// Named according to test and trial variables.
325 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> _integrated_bc_map;
326 : /// Arrays to store essential BCs to act on each component of weak form.
327 : /// Named according to test variable.
328 : NamedFieldsMap<std::vector<std::shared_ptr<MFEMEssentialBC>>> _essential_bc_map;
329 :
330 : // Operator handle for the jacobian
331 : mutable mfem::OperatorHandle _jacobian;
332 : // Operator handle for the linear components of the system operator
333 : mutable mfem::OperatorHandle _linear_operator;
334 : mfem::AssemblyLevel _assembly_level;
335 :
336 : // Pointer to GridFunctions to enable updates during nonlinear iterations
337 : Moose::MFEM::GridFunctions * _gfuncs;
338 : // Array storing block offsets of solution and residual vector
339 : mfem::Array<int> _block_true_offsets;
340 : // Non-owning pointer to the vector passed into the most recent GetGradient() call.
341 : mutable const mfem::Vector * _linearization_point = nullptr;
342 : // Boolean indicating if EquationSystem contains nonlinear integrators
343 : bool _non_linear = false;
344 : // Whether an external object (e.g. solver) requires Jacobian/gradient information.
345 : bool _gradient_required = false;
346 : // Coefficient manager notified when trial variables are updated, so that stored projections
347 : // of solution-dependent coefficients are invalidated.
348 : CoefficientManager * _coefficient_manager = nullptr;
349 :
350 : private:
351 : /// Disallowed inherited method
352 : using mfem::Operator::RecoverFEMSolution;
353 : };
354 :
355 : template <class FormType>
356 : void
357 6412 : EquationSystem::ApplyDomainBLFIntegrators(
358 : const std::string & trial_var_name,
359 : const std::string & test_var_name,
360 : std::shared_ptr<FormType> form,
361 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map,
362 : std::optional<mfem::real_t> scale_factor)
363 : {
364 6412 : if (kernels_map.Has(test_var_name) && kernels_map.Get(test_var_name)->Has(trial_var_name))
365 : {
366 6261 : auto kernels = kernels_map.GetRef(test_var_name).GetRef(trial_var_name);
367 13772 : for (auto & kernel : kernels)
368 : {
369 7511 : mfem::BilinearFormIntegrator * integ = kernel->createBFIntegrator();
370 :
371 7511 : if (integ)
372 : {
373 6308 : if (scale_factor.has_value())
374 1704 : integ = new ScaleIntegrator(integ, scale_factor.value(), true);
375 6308 : kernel->isSubdomainRestricted()
376 6308 : ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
377 6216 : : form->AddDomainIntegrator(std::move(integ));
378 : }
379 : }
380 6261 : }
381 6412 : }
382 :
383 : template <class FormType>
384 : void
385 2726 : EquationSystem::ApplyBoundaryBLFIntegrators(
386 : const std::string & trial_var_name,
387 : const std::string & test_var_name,
388 : std::shared_ptr<FormType> form,
389 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
390 : integrated_bc_map,
391 : std::optional<mfem::real_t> scale_factor)
392 : {
393 2864 : if (integrated_bc_map.Has(test_var_name) &&
394 138 : integrated_bc_map.Get(test_var_name)->Has(trial_var_name))
395 : {
396 116 : auto bcs = integrated_bc_map.GetRef(test_var_name).GetRef(trial_var_name);
397 250 : for (auto & bc : bcs)
398 : {
399 134 : mfem::BilinearFormIntegrator * integ = bc->createBFIntegrator();
400 :
401 134 : if (integ)
402 : {
403 36 : if (scale_factor.has_value())
404 36 : integ = new ScaleIntegrator(integ, scale_factor.value(), true);
405 36 : bc->isBoundaryRestricted()
406 36 : ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
407 0 : : form->AddBoundaryIntegrator(std::move(integ));
408 : }
409 : }
410 116 : }
411 2726 : }
412 :
413 : } // namespace Moose::MFEM
414 :
415 : #endif
|
