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 :
30 : /**
31 : * Owns the weak-form mathematics of a MOOSE MFEM problem.
32 : *
33 : * An EquationSystem stores weak-form components (bilinear, linear, mixed-bilinear,
34 : * and nonlinear forms) contributed by kernels and boundary conditions. It forms the
35 : * constrained linear part, keeps nonlinear action forms for residual/Jacobian
36 : * evaluation, and exposes the solve interface as an mfem::Operator. It is responsible
37 : * for applying essential-DoF constraints and propagating either itself or the assembled
38 : * linear operator (and bilinear form, for LOR) to the configured solver tree.
39 : *
40 : * EquationSystem is *not* responsible for grid-function bookkeeping, time stepping, or
41 : * solver selection - those belong to the ProblemOperator layer.
42 : *
43 : * @see ProblemOperatorBase for the conceptual split between the two layers.
44 : */
45 : class EquationSystem : public mfem::Operator
46 : {
47 :
48 : public:
49 1594 : EquationSystem() = default;
50 : ~EquationSystem() override;
51 :
52 : /// Add kernels.
53 : virtual void AddKernel(std::shared_ptr<MFEMKernel> kernel);
54 : virtual void AddIntegratedBC(std::shared_ptr<MFEMIntegratedBC> kernel);
55 : /// Add BC associated with essentially constrained DoFs on boundaries.
56 : virtual void AddEssentialBC(std::shared_ptr<MFEMEssentialBC> bc);
57 :
58 : /// Initialise
59 : virtual void Init(GridFunctions & gridfunctions,
60 : ComplexGridFunctions & cmplx_gridfunctions,
61 : mfem::AssemblyLevel assembly_level);
62 : /**
63 : * Build all weak-form components via BuildEquationSystem(), form the constrained linear part of
64 : * the system, and populate the true-DoF vectors used by the solve.
65 : *
66 : * For nonlinear problems, nonlinear forms are registered here but are not folded into the
67 : * assembled linear operator. Their residual action is evaluated by Mult(), and any Jacobian
68 : * contribution is formed at the current nonlinear iterate by GetGradient().
69 : *
70 : * This is the single public entry point for callers that want the equation system prepared for a
71 : * solve. Subclasses customize the form-building step by overriding BuildEquationSystem().
72 : */
73 : void FormSystem(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS);
74 : /// Compute residual y = Mu
75 : void Mult(const mfem::Vector & u, mfem::Vector & residual) const override;
76 : /// Compute the contribution to the residual from nonlinear forms only.
77 : virtual void ComputeNonlinearResidual(const mfem::Vector & u, mfem::Vector & residual) const;
78 : /// Get Jacobian at the provided vector of true DoFs of trial variables
79 : mfem::Operator & GetGradient(const mfem::Vector & u) const override;
80 :
81 : /// Update variable from solution vector after solve
82 : virtual void SetTrialVariablesFromTrueVectors(const mfem::BlockVector & trueX) const;
83 :
84 : /// Set whether the nonlinear solver driving this equation system requires Jacobian information.
85 56 : void SetSolverRequiresGradient(bool requires_gradient)
86 : {
87 56 : _solver_requires_gradient = requires_gradient;
88 56 : }
89 :
90 : /// Set the coefficient manager to notify when trial variables are updated, so that stored
91 : /// projections of solution-dependent coefficients are invalidated.
92 1576 : void SetCoefficientManager(CoefficientManager & coefficients)
93 : {
94 1576 : _coefficient_manager = &coefficients;
95 1576 : }
96 :
97 : // Test variables are associated with linear forms,
98 : // whereas trial variables are associated with gridfunctions.
99 1602 : const std::vector<std::string> & GetTrialVarNames() const { return _trial_var_names; }
100 1628 : const std::vector<std::string> & GetTestVarNames() const { return _test_var_names; }
101 :
102 : /**
103 : * @returns Whether nonlinear integrators are present
104 : */
105 2520 : bool Nonlinear() const { return _non_linear; }
106 :
107 : /**
108 : * Prepare the provided linear solver. First calls SetupLOR on the solver if it's using a Low
109 : * Order Refined methodology and then calls SetOperator on the solver with the assembled linear
110 : * operator
111 : */
112 : void PrepareLinearSolver(LinearSolverBase & solver);
113 :
114 : /// The true-DoF vector used for the most recent Jacobian linearization.
115 : const mfem::Vector & GetLinearizationPoint() const;
116 :
117 : /**
118 : * Build a fresh ParBilinearForm on the given FESpace using the same kernels as the main
119 : * system's bilinear form for var_name. Caller owns the returned form.
120 : */
121 : std::shared_ptr<mfem::ParBilinearForm>
122 : BuildBilinearFormForFESpace(const std::string & var_name,
123 : mfem::ParFiniteElementSpace & fespace,
124 : mfem::AssemblyLevel assembly_level);
125 :
126 : /**
127 : * Build a fresh ParNonlinearForm on the given FESpace using the same kernels as the main
128 : * system's nonlinear form for var_name. Caller owns the returned form.
129 : */
130 : std::shared_ptr<mfem::ParNonlinearForm>
131 : BuildNonlinearFormForFESpace(const std::string & var_name,
132 : mfem::ParFiniteElementSpace & fespace,
133 : mfem::AssemblyLevel assembly_level);
134 :
135 : /**
136 : * Returns true if the equation system has active mixed bilinear form contributions for var_name.
137 : */
138 : bool HasMixedBilinearForms(const std::string & var_name) const;
139 :
140 : /**
141 : * Build and return the essential boundary attribute marker array for a given trial variable.
142 : * The returned array has size == pmesh.bdr_attributes.Max() with 1 at essential boundaries.
143 : */
144 : mfem::Array<int> BuildEssentialBoundaryMarkers(const std::string & var_name) const;
145 :
146 : protected:
147 : /// Add coupled variable to EquationSystem.
148 : virtual void AddCoupledVariableNameIfMissing(const std::string & coupled_var_name);
149 : /// Add eliminated variable to EquationSystem.
150 : virtual void AddEliminatedVariableNameIfMissing(const std::string & eliminated_var_name);
151 : /// Add test variable to EquationSystem.
152 : virtual void AddTestVariableNameIfMissing(const std::string & test_var_name);
153 : /// Set trial variable names from subset of coupled variables that have an associated test variable.
154 : virtual void SetTrialVariableNames();
155 :
156 : /// Deletes the HypreParMatrix associated with any pointer stored in _h_blocks,
157 : /// and then proceeds to delete all dynamically allocated memory for _h_blocks
158 : /// itself, resetting all dimensions to zero.
159 : void DeleteHBlocks();
160 :
161 : /// Deletes the HypreParMatrix associated with any pointer stored in _jacobian_blocks,
162 : /// and then proceeds to delete all dynamically allocated memory for _jacobian_blocks
163 : /// itself, resetting all dimensions to zero.
164 : void DeleteJacobianBlocks();
165 :
166 : bool VectorContainsName(const std::vector<std::string> & the_vector,
167 : const std::string & name) const;
168 :
169 : /// Apply essential BC(s) associated with var_name to set true DoFs of trial_gf and update
170 : /// markers of all essential boundaries
171 : virtual void ApplyEssentialBC(const std::string & var_name,
172 : mfem::ParGridFunction & trial_gf,
173 : mfem::Array<int> & global_ess_markers);
174 : /// Update all essentially constrained true DoF markers and values on boundaries
175 : virtual void ApplyEssentialBCs();
176 : /// Perform trivial eliminations of coupled variables lacking corresponding test variables
177 : virtual void EliminateCoupledVariables();
178 : /// Build linear forms and eliminate constrained DoFs
179 : virtual void BuildLinearForms();
180 : /// Build non-linear action forms
181 : virtual void BuildNonlinearForms();
182 : /// Build bilinear forms (diagonal Jacobian contributions)
183 : virtual void BuildBilinearForms();
184 : /// Build mixed bilinear forms (off-diagonal Jacobian contributions)
185 : virtual void BuildMixedBilinearForms();
186 : /// Build all forms comprising this EquationSystem
187 : virtual void BuildEquationSystem();
188 :
189 : /// Form linear components of system based on on- and off-diagonal bilinear form
190 : /// contributions, populate solution and RHS vectors of true DoFs, and apply constraints.
191 : virtual void FormLinearSystem(mfem::OperatorHandle & op,
192 : mfem::BlockVector & trueX,
193 : mfem::BlockVector & trueRHS);
194 : using mfem::Operator::FormSystemOperator;
195 : /// Form matrix-free representation of linear components of system operator.
196 : /// Used when EquationSystem assembly level is set to 'FULL', 'ELEMENT', 'PARTIAL', or 'NONE'.
197 : virtual void FormSystemOperator(mfem::OperatorHandle & op,
198 : mfem::BlockVector & trueX,
199 : mfem::BlockVector & trueRHS);
200 : /// Form matrix representation of linear components of system operator as a HypreParMatrix.
201 : /// Used when EquationSystem assembly level is set to 'LEGACY'.
202 : virtual void FormSystemMatrix(mfem::OperatorHandle & op,
203 : mfem::BlockVector & trueX,
204 : mfem::BlockVector & trueRHS);
205 : /// Compute Jacobian matrix at the provided vector of true DoFs of trial variables
206 : void FormJacobianMatrix(const mfem::Vector & u);
207 :
208 : /**
209 : * Template method for applying BilinearFormIntegrators on domains from kernels to a BilinearForm,
210 : * or MixedBilinearForm
211 : */
212 : template <class FormType>
213 : void ApplyDomainBLFIntegrators(
214 : const std::string & trial_var_name,
215 : const std::string & test_var_name,
216 : std::shared_ptr<FormType> form,
217 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map,
218 : std::optional<mfem::real_t> scale_factor = std::nullopt);
219 :
220 : /**
221 : * Apply domain LinearFormIntegrators from kernels to the linear form associated with the
222 : * supplied test variable.
223 : */
224 : void ApplyDomainLFIntegrators(
225 : const std::string & test_var_name,
226 : std::shared_ptr<mfem::ParLinearForm> form,
227 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map);
228 :
229 : /**
230 : * Apply domain NonlinearFormIntegrators from kernels to the nonlinear form associated with the
231 : * supplied test variable.
232 : */
233 : void ApplyDomainNLFIntegrators(
234 : const std::string & test_var_name,
235 : std::shared_ptr<mfem::ParNonlinearForm> form,
236 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map,
237 : std::optional<mfem::real_t> scale_factor = std::nullopt);
238 :
239 : /**
240 : * Template method for applying BilinearFormIntegrators on boundaries from integrated boundary
241 : * conditions to a BilinearForm, or MixedBilinearForm.
242 : */
243 : template <class FormType>
244 : void ApplyBoundaryBLFIntegrators(
245 : const std::string & trial_var_name,
246 : const std::string & test_var_name,
247 : std::shared_ptr<FormType> form,
248 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
249 : integrated_bc_map,
250 : std::optional<mfem::real_t> scale_factor = std::nullopt);
251 :
252 : /**
253 : * Apply boundary LinearFormIntegrators from integrated boundary conditions to the linear form
254 : * associated with the supplied test variable.
255 : */
256 : void ApplyBoundaryLFIntegrators(
257 : const std::string & test_var_name,
258 : std::shared_ptr<mfem::ParLinearForm> form,
259 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
260 : integrated_bc_map);
261 :
262 : /**
263 : * Apply boundary NonlinearFormIntegrators from integrated boundary conditions to the nonlinear
264 : * form associated with the supplied test variable.
265 : */
266 : void ApplyBoundaryNLFIntegrators(
267 : const std::string & test_var_name,
268 : std::shared_ptr<mfem::ParNonlinearForm> form,
269 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
270 : integrated_bc_map,
271 : std::optional<mfem::real_t> scale_factor = std::nullopt);
272 :
273 : /**
274 : * Whether this a complex equation system
275 : */
276 52 : virtual bool Complex() const { return false; }
277 :
278 : /// Names of all trial variables of kernels and boundary conditions
279 : /// added to this EquationSystem.
280 : std::vector<std::string> _coupled_var_names;
281 : /// Subset of _coupled_var_names of all variables corresponding to gridfunctions with degrees of
282 : /// freedom that comprise the state vector of this EquationSystem. This will differ from
283 : /// _coupled_var_names when time derivatives or other eliminated variables are present.
284 : std::vector<std::string> _trial_var_names;
285 : /// Names of all coupled variables without a corresponding test variable.
286 : std::vector<std::string> _eliminated_var_names;
287 : /// Pointers to coupled variables not part of the reduced EquationSystem.
288 : Moose::MFEM::GridFunctions _eliminated_variables;
289 : /// Names of all test variables corresponding to linear forms in this equation system
290 : std::vector<std::string> _test_var_names;
291 : /// Pointers to finite element spaces associated with test variables.
292 : std::vector<mfem::ParFiniteElementSpace *> _test_pfespaces;
293 : /// Pointers to finite element spaces associated with coupled variables.
294 : std::vector<mfem::ParFiniteElementSpace *> _coupled_pfespaces;
295 :
296 : // Components of weak form, named according to test variable
297 : NamedFieldsMap<mfem::ParBilinearForm> _blfs;
298 : NamedFieldsMap<mfem::ParLinearForm> _lfs;
299 : NamedFieldsMap<mfem::ParNonlinearForm> _nlfs;
300 : NamedFieldsMap<NamedFieldsMap<mfem::ParMixedBilinearForm>> _mblfs; // named according to trial var
301 :
302 : /// Gridfunctions holding essential constraints from Dirichlet BCs
303 : std::vector<std::unique_ptr<mfem::ParGridFunction>> _var_ess_constraints;
304 : std::vector<mfem::Array<int>> _ess_tdof_lists;
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
|
