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