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 : #include "libmesh/ignore_warnings.h"
14 : #include "mfem/miniapps/common/pfem_extras.hpp"
15 : #include "libmesh/restore_warnings.h"
16 : #include "MFEMIntegratedBC.h"
17 : #include "MFEMEssentialBC.h"
18 : #include "MFEMContainers.h"
19 : #include "MFEMKernel.h"
20 : #include "MFEMMixedBilinearFormKernel.h"
21 : #include "ScaleIntegrator.h"
22 :
23 : namespace Moose::MFEM
24 : {
25 :
26 : /**
27 : * Class to store weak form components (bilinear and linear forms, and optionally
28 : * mixed and nonlinear forms) and build methods
29 : */
30 : class EquationSystem : public mfem::Operator
31 : {
32 :
33 : public:
34 : friend class EquationSystemProblemOperator;
35 : friend class TimeDomainEquationSystemProblemOperator;
36 :
37 353 : EquationSystem() = default;
38 : ~EquationSystem() override;
39 :
40 : /// Add test variable to EquationSystem.
41 : virtual void AddTestVariableNameIfMissing(const std::string & test_var_name);
42 : /// Add coupled variable to EquationSystem.
43 : virtual void AddCoupledVariableNameIfMissing(const std::string & coupled_var_name);
44 :
45 : /// Add kernels.
46 : virtual void AddKernel(std::shared_ptr<MFEMKernel> kernel);
47 : virtual void AddIntegratedBC(std::shared_ptr<MFEMIntegratedBC> kernel);
48 :
49 : /// Add BC associated with essentially constrained DoFs on boundaries.
50 : virtual void AddEssentialBC(std::shared_ptr<MFEMEssentialBC> bc);
51 :
52 : /// Initialise
53 : virtual void Init(Moose::MFEM::GridFunctions & gridfunctions, mfem::AssemblyLevel assembly_level);
54 :
55 : /// Build linear forms and eliminate constrained DoFs
56 : virtual void BuildLinearForms();
57 :
58 : /// Apply essential BC(s) associated with test_var_name to set true DoFs of trial_gf and update
59 : /// markers of all essential boundaries
60 : virtual void ApplyEssentialBC(const std::string & test_var_name,
61 : mfem::ParGridFunction & trial_gf,
62 : mfem::Array<int> & global_ess_markers);
63 : /// Update all essentially constrained true DoF markers and values on boundaries
64 : virtual void ApplyEssentialBCs();
65 :
66 : /// Perform trivial eliminations of coupled variables lacking corresponding test variables
67 : virtual void EliminateCoupledVariables();
68 :
69 : /// Build bilinear forms (diagonal Jacobian contributions)
70 : virtual void BuildBilinearForms();
71 : /// Build mixed bilinear forms (off-diagonal Jacobian contributions)
72 : virtual void BuildMixedBilinearForms();
73 : /// Build all forms comprising this EquationSystem
74 : virtual void BuildEquationSystem();
75 :
76 : /// Form Jacobian operator based on on- and off-diagonal bilinear form contributions, populate
77 : /// solution and RHS vectors of true DoFs, and apply constraints
78 : void AssembleJacobian(
79 : Moose::MFEM::NamedFieldsMap<mfem::ParBilinearForm> & jac_blfs,
80 : Moose::MFEM::NamedFieldsMap<Moose::MFEM::NamedFieldsMap<mfem::ParMixedBilinearForm>> &
81 : jac_mblfs,
82 : Moose::MFEM::NamedFieldsMap<mfem::ParLinearForm> & rhs_lfs,
83 : std::vector<mfem::Array<int>> & ess_tdof_lists,
84 : std::vector<std::unique_ptr<mfem::ParGridFunction>> & var_ess_constraints,
85 : mfem::OperatorHandle & op,
86 : mfem::BlockVector & trueX,
87 : mfem::BlockVector & trueRHS);
88 :
89 : /// Form linear system, with essential boundary conditions accounted for
90 : virtual void FormLinearSystem(mfem::OperatorHandle & op,
91 : mfem::BlockVector & trueX,
92 : mfem::BlockVector & trueRHS);
93 :
94 : virtual void
95 : FormSystem(mfem::OperatorHandle & op, mfem::BlockVector & trueX, mfem::BlockVector & trueRHS);
96 : virtual void FormLegacySystem(mfem::OperatorHandle & op,
97 : mfem::BlockVector & trueX,
98 : mfem::BlockVector & trueRHS);
99 :
100 : /// Build linear system, with essential boundary conditions accounted for
101 : virtual void BuildJacobian(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS);
102 :
103 : /// Compute residual y = Mu
104 : void Mult(const mfem::Vector & u, mfem::Vector & residual) const override;
105 :
106 : /// Compute J = M + grad_H(u)
107 : mfem::Operator & GetGradient(const mfem::Vector & u) const override;
108 :
109 : /// Update variable from solution vector after solve
110 : virtual void RecoverFEMSolution(mfem::BlockVector & trueX,
111 : Moose::MFEM::GridFunctions & gridfunctions);
112 :
113 : std::vector<mfem::Array<int>> _ess_tdof_lists;
114 :
115 353 : const std::vector<std::string> & TrialVarNames() const { return _trial_var_names; }
116 353 : const std::vector<std::string> & TestVarNames() const { return _test_var_names; }
117 :
118 : private:
119 : /// Disallowed inherited method
120 : using mfem::Operator::RecoverFEMSolution;
121 :
122 : /// Set trial variable names from subset of coupled variables that have an associated test variable.
123 : virtual void SetTrialVariableNames();
124 :
125 : protected:
126 : /// Deletes the HypreParMatrix associated with any pointer stored in _h_blocks,
127 : /// and then proceeds to delete all dynamically allocated memory for _h_blocks
128 : /// itself, resetting all dimensions to zero.
129 : void DeleteAllBlocks();
130 :
131 : bool VectorContainsName(const std::vector<std::string> & the_vector,
132 : const std::string & name) const;
133 :
134 : // Test variables are associated with LinearForms,
135 : // whereas trial variables are associated with gridfunctions.
136 :
137 : /// Names of all trial variables of kernels and boundary conditions
138 : /// added to this EquationSystem.
139 : std::vector<std::string> _coupled_var_names;
140 : /// Subset of _coupled_var_names of all variables corresponding to gridfunctions with degrees of
141 : /// freedom that comprise the state vector of this EquationSystem. This will differ from
142 : /// _coupled_var_names when time derivatives or other eliminated variables are present.
143 : std::vector<std::string> _trial_var_names;
144 : /// Names of all coupled variables without a corresponding test variable.
145 : std::vector<std::string> _eliminated_var_names;
146 : /// Pointers to coupled variables not part of the reduced EquationSystem.
147 : Moose::MFEM::GridFunctions _eliminated_variables;
148 : /// Names of all test variables corresponding to linear forms in this equation system
149 : std::vector<std::string> _test_var_names;
150 : /// Pointers to finite element spaces associated with test variables.
151 : std::vector<mfem::ParFiniteElementSpace *> _test_pfespaces;
152 : /// Pointers to finite element spaces associated with coupled variables.
153 : std::vector<mfem::ParFiniteElementSpace *> _coupled_pfespaces;
154 :
155 : // Components of weak form. // Named according to test variable
156 : Moose::MFEM::NamedFieldsMap<mfem::ParBilinearForm> _blfs;
157 : Moose::MFEM::NamedFieldsMap<mfem::ParLinearForm> _lfs;
158 : Moose::MFEM::NamedFieldsMap<mfem::ParNonlinearForm> _nlfs;
159 : Moose::MFEM::NamedFieldsMap<Moose::MFEM::NamedFieldsMap<mfem::ParMixedBilinearForm>>
160 : _mblfs; // named according to trial variable
161 :
162 : /**
163 : * Template method for applying BilinearFormIntegrators on domains from kernels to a BilinearForm,
164 : * or MixedBilinearForm
165 : */
166 : template <class FormType>
167 : void ApplyDomainBLFIntegrators(
168 : const std::string & trial_var_name,
169 : const std::string & test_var_name,
170 : std::shared_ptr<FormType> form,
171 : Moose::MFEM::NamedFieldsMap<
172 : Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map);
173 :
174 : void ApplyDomainLFIntegrators(
175 : const std::string & test_var_name,
176 : std::shared_ptr<mfem::ParLinearForm> form,
177 : Moose::MFEM::NamedFieldsMap<
178 : Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map);
179 :
180 : template <class FormType>
181 : void ApplyBoundaryBLFIntegrators(
182 : const std::string & trial_var_name,
183 : const std::string & test_var_name,
184 : std::shared_ptr<FormType> form,
185 : Moose::MFEM::NamedFieldsMap<
186 : Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
187 : integrated_bc_map);
188 :
189 : void ApplyBoundaryLFIntegrators(
190 : const std::string & test_var_name,
191 : std::shared_ptr<mfem::ParLinearForm> form,
192 : Moose::MFEM::NamedFieldsMap<
193 : Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
194 : integrated_bc_map);
195 :
196 : /// Gridfunctions holding essential constraints from Dirichlet BCs
197 : std::vector<std::unique_ptr<mfem::ParGridFunction>> _var_ess_constraints;
198 :
199 : mfem::Array2D<const mfem::HypreParMatrix *> _h_blocks;
200 :
201 : /// Arrays to store kernels to act on each component of weak form.
202 : /// Named according to test and trial variables.
203 : Moose::MFEM::NamedFieldsMap<Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>>
204 : _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 : Moose::MFEM::NamedFieldsMap<
208 : Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>>
209 : _integrated_bc_map;
210 : /// Arrays to store essential BCs to act on each component of weak form.
211 : /// Named according to test variable.
212 : Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMEssentialBC>>> _essential_bc_map;
213 :
214 : mutable mfem::OperatorHandle _jacobian;
215 :
216 : mfem::AssemblyLevel _assembly_level;
217 : };
218 :
219 : template <class FormType>
220 : void
221 1852 : EquationSystem::ApplyDomainBLFIntegrators(
222 : const std::string & trial_var_name,
223 : const std::string & test_var_name,
224 : std::shared_ptr<FormType> form,
225 : Moose::MFEM::NamedFieldsMap<
226 : Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map)
227 : {
228 1852 : if (kernels_map.Has(test_var_name) && kernels_map.Get(test_var_name)->Has(trial_var_name))
229 : {
230 1798 : auto kernels = kernels_map.GetRef(test_var_name).GetRef(trial_var_name);
231 3685 : for (auto & kernel : kernels)
232 : {
233 1887 : mfem::BilinearFormIntegrator * integ = kernel->createBFIntegrator();
234 1887 : if (integ != nullptr)
235 : {
236 1836 : kernel->isSubdomainRestricted()
237 1836 : ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
238 1800 : : form->AddDomainIntegrator(std::move(integ));
239 : }
240 : }
241 1798 : }
242 1852 : }
243 :
244 : inline void
245 982 : EquationSystem::ApplyDomainLFIntegrators(
246 : const std::string & test_var_name,
247 : std::shared_ptr<mfem::ParLinearForm> form,
248 : Moose::MFEM::NamedFieldsMap<
249 : Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map)
250 : {
251 982 : if (kernels_map.Has(test_var_name) && kernels_map.Get(test_var_name)->Has(test_var_name))
252 : {
253 942 : auto kernels = kernels_map.GetRef(test_var_name).GetRef(test_var_name);
254 1973 : for (auto & kernel : kernels)
255 : {
256 1031 : mfem::LinearFormIntegrator * integ = kernel->createLFIntegrator();
257 1031 : if (integ != nullptr)
258 : {
259 51 : kernel->isSubdomainRestricted()
260 51 : ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
261 51 : : form->AddDomainIntegrator(std::move(integ));
262 : }
263 : }
264 942 : }
265 982 : }
266 :
267 : template <class FormType>
268 : void
269 1768 : EquationSystem::ApplyBoundaryBLFIntegrators(
270 : const std::string & trial_var_name,
271 : const std::string & test_var_name,
272 : std::shared_ptr<FormType> form,
273 : Moose::MFEM::NamedFieldsMap<
274 : Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
275 : integrated_bc_map)
276 : {
277 1928 : if (integrated_bc_map.Has(test_var_name) &&
278 160 : integrated_bc_map.Get(test_var_name)->Has(trial_var_name))
279 : {
280 160 : auto bcs = integrated_bc_map.GetRef(test_var_name).GetRef(trial_var_name);
281 352 : for (auto & bc : bcs)
282 : {
283 192 : mfem::BilinearFormIntegrator * integ = bc->createBFIntegrator();
284 192 : if (integ != nullptr)
285 : {
286 112 : bc->isBoundaryRestricted()
287 112 : ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
288 0 : : form->AddBoundaryIntegrator(std::move(integ));
289 : }
290 : }
291 160 : }
292 1768 : }
293 :
294 : inline void
295 982 : EquationSystem::ApplyBoundaryLFIntegrators(
296 : const std::string & test_var_name,
297 : std::shared_ptr<mfem::ParLinearForm> form,
298 : Moose::MFEM::NamedFieldsMap<
299 : Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
300 : integrated_bc_map)
301 : {
302 1070 : if (integrated_bc_map.Has(test_var_name) &&
303 88 : integrated_bc_map.Get(test_var_name)->Has(test_var_name))
304 : {
305 88 : auto bcs = integrated_bc_map.GetRef(test_var_name).GetRef(test_var_name);
306 192 : for (auto & bc : bcs)
307 : {
308 104 : mfem::LinearFormIntegrator * integ = bc->createLFIntegrator();
309 104 : if (integ != nullptr)
310 : {
311 104 : bc->isBoundaryRestricted()
312 104 : ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
313 8 : : form->AddBoundaryIntegrator(std::move(integ));
314 : }
315 : }
316 88 : }
317 982 : }
318 :
319 : /**
320 : * Class to store weak form components for time dependent PDEs
321 : */
322 : class TimeDependentEquationSystem : public EquationSystem
323 : {
324 : public:
325 : TimeDependentEquationSystem(const Moose::MFEM::TimeDerivativeMap & time_derivative_map);
326 :
327 : /// Initialise
328 : virtual void Init(Moose::MFEM::GridFunctions & gridfunctions,
329 : mfem::AssemblyLevel assembly_level) override;
330 :
331 : virtual void SetTimeStep(mfem::real_t dt);
332 : virtual void UpdateEquationSystem();
333 :
334 : virtual void AddKernel(std::shared_ptr<MFEMKernel> kernel) override;
335 : virtual void BuildBilinearForms() override;
336 : virtual void BuildMixedBilinearForms() override;
337 : virtual void ApplyEssentialBCs() override;
338 : virtual void EliminateCoupledVariables() override;
339 : virtual void FormLegacySystem(mfem::OperatorHandle & op,
340 : mfem::BlockVector & truedXdt,
341 : mfem::BlockVector & trueRHS) override;
342 : virtual void FormSystem(mfem::OperatorHandle & op,
343 : mfem::BlockVector & truedXdt,
344 : mfem::BlockVector & trueRHS) override;
345 :
346 : /// Fetch all integrators on a source bilinear form, scale them by a real factor, and add to a second target bilienar form.
347 : /// Useful for scaling bilinear form integrators by timesteps.
348 : template <class FormType>
349 802 : void ScaleAndAddBLFIntegrators(std::shared_ptr<FormType> source_blf,
350 : std::shared_ptr<FormType> target_blf,
351 : mfem::real_t scale_factor)
352 : {
353 : // Add and scale all domain integrators existing on source_blf to target_blf
354 802 : auto domain_integrators = source_blf->GetDBFI();
355 802 : auto domain_markers = source_blf->GetDBFI_Marker();
356 1572 : for (int i = 0; i < domain_integrators->Size(); ++i)
357 1540 : target_blf->AddDomainIntegrator(
358 770 : new ScaleIntegrator(*domain_integrators[i], scale_factor, false), *(*domain_markers[i]));
359 : // Add and scale all boundary integrators existing on source_blf to target_blf
360 802 : auto boundary_integrators = source_blf->GetBBFI();
361 802 : auto boundary_markers = source_blf->GetBBFI_Marker();
362 858 : for (int i = 0; i < boundary_integrators->Size(); ++i)
363 112 : target_blf->AddBoundaryIntegrator(
364 56 : new ScaleIntegrator(*boundary_integrators[i], scale_factor, false),
365 56 : *(*boundary_markers[i]));
366 802 : }
367 :
368 : protected:
369 : /// Coefficient for timestep scaling
370 : mfem::ConstantCoefficient _dt_coef;
371 :
372 : Moose::MFEM::NamedFieldsMap<Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>>
373 : _td_kernels_map;
374 : /// Containers to store contributions to weak form of the form (F du/dt, v)
375 : Moose::MFEM::NamedFieldsMap<mfem::ParBilinearForm> _td_blfs;
376 : Moose::MFEM::NamedFieldsMap<Moose::MFEM::NamedFieldsMap<mfem::ParMixedBilinearForm>>
377 : _td_mblfs; // named according to trial variable
378 :
379 : /// Gridfunctions holding essential constraints from Dirichlet BCs
380 : std::vector<std::unique_ptr<mfem::ParGridFunction>> _td_var_ess_constraints;
381 :
382 : /// Map between variable names and their time derivatives
383 : const Moose::MFEM::TimeDerivativeMap & _time_derivative_map;
384 :
385 : private:
386 : /// Set trial variable names from subset of coupled variables that have an associated test variable.
387 : virtual void SetTrialVariableNames() override;
388 : };
389 :
390 : } // namespace Moose::MFEM
391 :
392 : #endif
|
