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