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