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