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 : #include "EquationSystem.h"
13 : #include "MFEMLinearSolverBase.h"
14 : #include "CoefficientManager.h"
15 : #include "libmesh/int_range.h"
16 :
17 : namespace Moose::MFEM
18 : {
19 :
20 1673 : EquationSystem::~EquationSystem()
21 : {
22 1673 : DeleteHBlocks();
23 1673 : DeleteJacobianBlocks();
24 1673 : }
25 :
26 : void
27 4266 : EquationSystem::DeleteHBlocks()
28 : {
29 6888 : for (const auto i : make_range(_h_blocks.NumRows()))
30 5306 : for (const auto j : make_range(_h_blocks.NumCols()))
31 : {
32 2684 : if (_jacobian_blocks.NumRows() && _jacobian_blocks(i, j) == _h_blocks(i, j))
33 0 : _jacobian_blocks(i, j) = nullptr;
34 2684 : delete _h_blocks(i, j);
35 : }
36 4266 : _h_blocks.DeleteAll();
37 4266 : }
38 :
39 : void
40 3372 : EquationSystem::DeleteJacobianBlocks()
41 : {
42 5069 : for (const auto i : make_range(_jacobian_blocks.NumRows()))
43 3394 : for (const auto j : make_range(_jacobian_blocks.NumCols()))
44 1697 : if (!_h_blocks.NumRows() || _jacobian_blocks(i, j) != _h_blocks(i, j))
45 1697 : delete _jacobian_blocks(i, j);
46 3372 : _jacobian_blocks.DeleteAll();
47 3372 : }
48 :
49 : bool
50 8342 : EquationSystem::VectorContainsName(const std::vector<std::string> & the_vector,
51 : const std::string & name) const
52 : {
53 8342 : return std::find(the_vector.begin(), the_vector.end(), name) != the_vector.end();
54 : }
55 :
56 : void
57 1975 : EquationSystem::AddCoupledVariableNameIfMissing(const std::string & coupled_var_name)
58 : {
59 1975 : if (!VectorContainsName(_coupled_var_names, coupled_var_name))
60 1260 : _coupled_var_names.push_back(coupled_var_name);
61 1975 : }
62 :
63 : void
64 201 : EquationSystem::AddEliminatedVariableNameIfMissing(const std::string & eliminated_var_name)
65 : {
66 201 : if (!VectorContainsName(_eliminated_var_names, eliminated_var_name))
67 195 : _eliminated_var_names.push_back(eliminated_var_name);
68 201 : }
69 :
70 : void
71 3618 : EquationSystem::AddTestVariableNameIfMissing(const std::string & test_var_name)
72 : {
73 3618 : if (!VectorContainsName(_test_var_names, test_var_name))
74 1162 : _test_var_names.push_back(test_var_name);
75 3618 : }
76 :
77 : void
78 1673 : EquationSystem::SetTrialVariableNames()
79 : {
80 : // If a coupled variable has an equation associated with it,
81 : // add it to the set of trial variables.
82 2833 : for (const auto & test_var_name : _test_var_names)
83 1160 : if (VectorContainsName(_coupled_var_names, test_var_name))
84 1160 : _trial_var_names.push_back(test_var_name);
85 :
86 : // Otherwise, add it to the set of eliminated variables.
87 2931 : for (const auto & coupled_var_name : _coupled_var_names)
88 1258 : if (!VectorContainsName(_test_var_names, coupled_var_name))
89 98 : _eliminated_var_names.push_back(coupled_var_name);
90 1673 : }
91 :
92 : void
93 1832 : EquationSystem::AddKernel(std::shared_ptr<MFEMKernel> kernel)
94 : {
95 1832 : const auto & trial_var_name = kernel->getTrialVariableName();
96 1832 : const auto & test_var_name = kernel->getTestVariableName();
97 1832 : AddCoupledVariableNameIfMissing(trial_var_name);
98 1832 : AddTestVariableNameIfMissing(test_var_name);
99 : // Register new kernels map if not present for the test variable
100 1832 : if (!_kernels_map.Has(test_var_name))
101 : {
102 : auto kernel_field_map =
103 1114 : std::make_shared<Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>>();
104 1114 : _kernels_map.Register(test_var_name, std::move(kernel_field_map));
105 1114 : }
106 : // Register new kernels map if not present for the test/trial variable pair
107 1832 : if (!_kernels_map.Get(test_var_name)->Has(trial_var_name))
108 : {
109 1225 : auto kernels = std::make_shared<std::vector<std::shared_ptr<MFEMKernel>>>();
110 1225 : _kernels_map.Get(test_var_name)->Register(trial_var_name, std::move(kernels));
111 1225 : }
112 1832 : _kernels_map.GetRef(test_var_name).Get(trial_var_name)->push_back(std::move(kernel));
113 1832 : }
114 :
115 : void
116 66 : EquationSystem::AddIntegratedBC(std::shared_ptr<MFEMIntegratedBC> bc)
117 : {
118 66 : const auto & trial_var_name = bc->getTrialVariableName();
119 66 : const auto & test_var_name = bc->getTestVariableName();
120 66 : AddCoupledVariableNameIfMissing(trial_var_name);
121 66 : AddTestVariableNameIfMissing(test_var_name);
122 : // Register new integrated bc map if not present for the test variable
123 66 : if (!_integrated_bc_map.Has(test_var_name))
124 : {
125 : auto integrated_bc_field_map = std::make_shared<
126 60 : Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>>();
127 60 : _integrated_bc_map.Register(test_var_name, std::move(integrated_bc_field_map));
128 60 : }
129 : // Register new integrated bc map if not present for the test/trial variable pair
130 66 : if (!_integrated_bc_map.Get(test_var_name)->Has(trial_var_name))
131 : {
132 60 : auto bcs = std::make_shared<std::vector<std::shared_ptr<MFEMIntegratedBC>>>();
133 60 : _integrated_bc_map.Get(test_var_name)->Register(trial_var_name, std::move(bcs));
134 60 : }
135 66 : _integrated_bc_map.GetRef(test_var_name).Get(trial_var_name)->push_back(std::move(bc));
136 66 : }
137 :
138 : void
139 1379 : EquationSystem::AddEssentialBC(std::shared_ptr<MFEMEssentialBC> bc)
140 : {
141 1379 : const auto & test_var_name = bc->getTestVariableName();
142 1379 : AddTestVariableNameIfMissing(test_var_name);
143 : // Register new essential bc map if not present for the test variable
144 1379 : if (!_essential_bc_map.Has(test_var_name))
145 : {
146 1011 : auto bcs = std::make_shared<std::vector<std::shared_ptr<MFEMEssentialBC>>>();
147 1011 : _essential_bc_map.Register(test_var_name, std::move(bcs));
148 1011 : }
149 1379 : _essential_bc_map.GetRef(test_var_name).push_back(std::move(bc));
150 1379 : }
151 :
152 : void
153 1610 : EquationSystem::Init(Moose::MFEM::GridFunctions & gridfunctions,
154 : Moose::MFEM::ComplexGridFunctions & /*cmplx_gridfunctions*/,
155 : mfem::AssemblyLevel assembly_level)
156 : {
157 1610 : _assembly_level = assembly_level;
158 :
159 : // Extract which coupled variables are to be trivially eliminated and which are trial variables
160 1610 : SetTrialVariableNames();
161 :
162 2728 : for (auto & test_var_name : _test_var_names)
163 : {
164 1118 : if (!gridfunctions.Has(test_var_name))
165 : {
166 0 : mooseError("MFEM variable ",
167 : test_var_name,
168 : " requested by equation system during initialization was "
169 : "not found in gridfunctions");
170 : }
171 : // Store pointers to test FESpaces
172 1118 : _test_pfespaces.push_back(gridfunctions.Get(test_var_name)->ParFESpace());
173 : }
174 :
175 2728 : for (auto & trial_var_name : _trial_var_names)
176 : {
177 1118 : if (!gridfunctions.Has(trial_var_name))
178 : {
179 0 : mooseError("MFEM variable ",
180 : trial_var_name,
181 : " requested by equation system during initialization was "
182 : "not found in gridfunctions");
183 : }
184 : // Create auxiliary gridfunctions for storing essential constraints from Dirichlet conditions
185 1118 : _var_ess_constraints.emplace_back(
186 2236 : std::make_unique<mfem::ParGridFunction>(gridfunctions.Get(trial_var_name)->ParFESpace()));
187 : }
188 :
189 : // Store pointers to FESpaces of all coupled variables
190 2826 : for (auto & coupled_var_name : _coupled_var_names)
191 1216 : _coupled_pfespaces.push_back(gridfunctions.Get(coupled_var_name)->ParFESpace());
192 :
193 : // Store pointers to coupled variable GridFunctions that are to be eliminated prior to forming the
194 : // jacobian
195 1901 : for (auto & eliminated_var_name : _eliminated_var_names)
196 291 : _eliminated_variables.Register(eliminated_var_name,
197 582 : gridfunctions.GetShared(eliminated_var_name));
198 :
199 : // Get a reference to the GridFunctions
200 1610 : _gfuncs = &gridfunctions;
201 1610 : }
202 :
203 : void
204 2655 : EquationSystem::ApplyEssentialBC(const std::string & var_name,
205 : mfem::ParGridFunction & trial_gf,
206 : mfem::Array<int> & global_ess_markers)
207 : {
208 2655 : if (_essential_bc_map.Has(var_name))
209 4911 : for (auto & bc : _essential_bc_map.GetRef(var_name))
210 : {
211 : // Set constrained DoFs values on essential boundaries
212 3039 : bc->ApplyBC(trial_gf);
213 : // Fetch marker array labelling essential boundaries of current BC
214 3039 : mfem::Array<int> ess_bdrs(bc->getBoundaryMarkers());
215 : // Add these boundary markers to the set of markers labelling all essential boundaries
216 15260 : for (const auto i : make_range(ess_bdrs.Size()))
217 12221 : global_ess_markers[i] |= ess_bdrs[i];
218 3039 : }
219 2655 : }
220 :
221 : void
222 2620 : EquationSystem::ApplyEssentialBCs()
223 : {
224 2620 : _ess_tdof_lists.resize(_trial_var_names.size());
225 2620 : _ess_markers.resize(_trial_var_names.size());
226 5275 : for (const auto i : index_range(_trial_var_names))
227 : {
228 2655 : const auto & trial_var_name = _trial_var_names.at(i);
229 2655 : mfem::ParGridFunction & trial_gf = *_var_ess_constraints.at(i);
230 :
231 : // Make sure we update the size, if this mesh has changed recently for instance
232 2655 : trial_gf.Update();
233 :
234 : // Initial guess for iterative solvers (initial condition or the previous time step solution)
235 2655 : trial_gf = _gfuncs->GetRef(trial_var_name);
236 :
237 2655 : _ess_markers.at(i).SetSize(trial_gf.ParFESpace()->GetParMesh()->bdr_attributes.Max(), 0);
238 : // Set strongly constrained DoFs of trial_gf on essential boundaries and add markers for all
239 : // essential boundaries to the _ess_markers array
240 2655 : ApplyEssentialBC(trial_var_name, trial_gf, _ess_markers.at(i));
241 2655 : trial_gf.ParFESpace()->GetEssentialTrueDofs(_ess_markers.at(i), _ess_tdof_lists.at(i));
242 : }
243 2620 : }
244 :
245 : void
246 2646 : EquationSystem::EliminateCoupledVariables()
247 : {
248 5327 : for (const auto & test_var_name : _test_var_names)
249 4535 : for (const auto & eliminated_var_name : _eliminated_var_names)
250 1984 : if (_mblfs.Has(test_var_name) && _mblfs.Get(test_var_name)->Has(eliminated_var_name) &&
251 130 : !VectorContainsName(_test_var_names, eliminated_var_name))
252 : {
253 94 : auto & mblf = *_mblfs.Get(test_var_name)->Get(eliminated_var_name);
254 94 : mblf.AddMult(*_eliminated_variables.Get(eliminated_var_name), *_lfs.Get(test_var_name), -1);
255 : }
256 2646 : }
257 :
258 : void
259 2652 : EquationSystem::FormLinearSystem(mfem::OperatorHandle & op,
260 : mfem::BlockVector & trueX,
261 : mfem::BlockVector & trueRHS)
262 : {
263 : mooseAssert(_test_var_names.size() == _trial_var_names.size(),
264 : "Number of test and trial variables must be the same for block matrix assembly.");
265 :
266 2652 : if (_assembly_level == mfem::AssemblyLevel::LEGACY)
267 2593 : FormSystemMatrix(op, trueX, trueRHS);
268 : else
269 59 : FormSystemOperator(op, trueX, trueRHS);
270 2652 : }
271 :
272 : void
273 59 : EquationSystem::FormSystemOperator(mfem::OperatorHandle & op,
274 : mfem::BlockVector & trueX,
275 : mfem::BlockVector & trueRHS)
276 : {
277 : mooseAssert(_test_var_names.size() == 1 && _test_var_names.size() == _trial_var_names.size(),
278 : "Non-legacy assembly is only supported for single test and trial variable systems");
279 :
280 59 : auto & test_var_name = _test_var_names.at(0);
281 59 : mfem::Vector aux_x, aux_rhs;
282 59 : mfem::OperatorPtr aux_a;
283 :
284 59 : auto blf = _blfs.Get(test_var_name);
285 59 : blf->FormLinearSystem(_ess_tdof_lists.at(0),
286 59 : *_var_ess_constraints.at(0),
287 59 : *_lfs.Get(test_var_name),
288 : aux_a,
289 : aux_x,
290 : aux_rhs,
291 : /*copy_interior=*/true);
292 :
293 59 : trueX.GetBlock(0) = aux_x;
294 59 : trueRHS.GetBlock(0) = aux_rhs;
295 59 : trueX.SyncFromBlocks();
296 59 : trueRHS.SyncFromBlocks();
297 :
298 59 : op.Reset(aux_a.Ptr());
299 59 : aux_a.SetOperatorOwner(false);
300 59 : }
301 :
302 : void
303 2551 : EquationSystem::FormSystemMatrix(mfem::OperatorHandle & op,
304 : mfem::BlockVector & trueX,
305 : mfem::BlockVector & trueRHS)
306 : {
307 : // Allocate block operator
308 2551 : DeleteHBlocks();
309 2551 : _h_blocks.SetSize(_test_var_names.size(), _trial_var_names.size());
310 2551 : _h_blocks = nullptr;
311 : // Zero out RHS and sync memory
312 2551 : trueRHS = 0.0;
313 2551 : trueRHS.SyncToBlocks();
314 :
315 5133 : for (const auto i : index_range(_test_var_names))
316 : {
317 2582 : auto test_var_name = _test_var_names.at(i);
318 :
319 5226 : for (const auto j : index_range(_trial_var_names))
320 : {
321 2644 : auto trial_var_name = _trial_var_names.at(j);
322 :
323 2644 : mfem::Vector aux_x, aux_rhs;
324 2644 : mfem::ParLinearForm aux_lf(_test_pfespaces.at(i));
325 2644 : mfem::HypreParMatrix * aux_a = new mfem::HypreParMatrix;
326 :
327 2644 : if (test_var_name == trial_var_name)
328 : {
329 : mooseAssert(i == j, "Trial and test variables must have the same ordering.");
330 2582 : auto blf = _blfs.Get(test_var_name);
331 2582 : blf->FormLinearSystem(_ess_tdof_lists.at(j),
332 2582 : *_var_ess_constraints.at(j),
333 2582 : *_lfs.Get(test_var_name),
334 : *aux_a,
335 : aux_x,
336 : aux_rhs,
337 : /*copy_interior=*/true);
338 2582 : trueX.GetBlock(j) = aux_x;
339 : }
340 62 : else if (_mblfs.Has(test_var_name) && _mblfs.Get(test_var_name)->Has(trial_var_name))
341 : {
342 62 : auto mblf = _mblfs.Get(test_var_name)->Get(trial_var_name);
343 62 : mblf->FormRectangularLinearSystem(_ess_tdof_lists.at(j),
344 62 : _ess_tdof_lists.at(i),
345 62 : *_var_ess_constraints.at(j),
346 62 : aux_lf = 0,
347 : *aux_a,
348 : aux_x,
349 : aux_rhs);
350 : }
351 : else
352 0 : continue;
353 :
354 2644 : trueRHS.GetBlock(i) += aux_rhs;
355 2644 : _h_blocks(i, j) = aux_a;
356 2644 : }
357 2582 : }
358 : // Sync memory
359 2551 : trueX.SyncFromBlocks();
360 2551 : trueRHS.SyncFromBlocks();
361 :
362 : // Create monolithic matrix
363 2551 : op.Reset(mfem::HypreParMatrixFromBlocks(_h_blocks));
364 2551 : }
365 :
366 : void
367 2654 : EquationSystem::FormSystem(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS)
368 : {
369 2654 : BuildEquationSystem();
370 2652 : height = trueX.Size();
371 2652 : width = trueRHS.Size();
372 : // Store block offsets
373 2652 : _block_true_offsets.SetSize(trueX.NumBlocks() + 1);
374 2652 : _block_true_offsets[0] = 0;
375 5335 : for (unsigned i = 0; i < _trial_var_names.size(); i++)
376 2683 : _block_true_offsets[i + 1] = trueX.BlockSize(i);
377 2652 : _block_true_offsets.PartialSum();
378 2652 : FormLinearSystem(_linear_operator, trueX, trueRHS);
379 2652 : }
380 :
381 : void
382 1697 : EquationSystem::Mult(const mfem::Vector & sol, mfem::Vector & residual) const
383 : {
384 1697 : if (_non_linear)
385 : {
386 1697 : ComputeNonlinearResidual(sol, residual);
387 1697 : _linear_operator->AddMult(sol, residual);
388 : }
389 : else
390 : {
391 0 : residual = 0.0;
392 0 : _linear_operator->Mult(sol, residual);
393 : }
394 :
395 1697 : sol.HostRead();
396 1697 : residual.HostRead();
397 1697 : }
398 :
399 : void
400 1697 : EquationSystem::ComputeNonlinearResidual(const mfem::Vector & sol, mfem::Vector & residual) const
401 : {
402 : mooseAssert(_non_linear, "Should not be calling this method if our forms are not nonlinear");
403 1697 : residual = 0.0;
404 :
405 1697 : const mfem::BlockVector block_solution(const_cast<mfem::Vector &>(sol), _block_true_offsets);
406 1697 : SetTrialVariablesFromTrueVectors(block_solution);
407 :
408 1697 : mfem::BlockVector block_residual(residual, _block_true_offsets);
409 3394 : for (unsigned int i = 0; i < _test_var_names.size(); i++)
410 : {
411 1697 : auto & test_var_name = _test_var_names.at(i);
412 1697 : auto nlf = _nlfs.GetShared(test_var_name);
413 1697 : nlf->Mult(block_solution.GetBlock(i), block_residual.GetBlock(i));
414 1697 : block_residual.GetBlock(i).SyncAliasMemory(block_residual);
415 1697 : }
416 1697 : }
417 :
418 : void
419 1699 : EquationSystem::FormJacobianMatrix(const mfem::Vector & u)
420 : {
421 1699 : DeleteJacobianBlocks();
422 1699 : _jacobian_blocks.SetSize(_test_var_names.size(), _trial_var_names.size());
423 1699 : _jacobian_blocks = nullptr;
424 :
425 1699 : const mfem::BlockVector update_vector(const_cast<mfem::Vector &>(u), _block_true_offsets);
426 3398 : for (const auto i : index_range(_test_var_names))
427 : {
428 1699 : auto test_var_name = _test_var_names.at(i);
429 1699 : if (_nlfs.Has(test_var_name))
430 : {
431 1699 : auto nlf = _nlfs.Get(test_var_name);
432 : mfem::HypreParMatrix * nlf_jac =
433 1699 : dynamic_cast<mfem::HypreParMatrix *>(&nlf->GetGradient(update_vector.GetBlock(i)));
434 : mooseAssert(nlf_jac,
435 : "Jacobian contribution of nonlinear form associated with " + test_var_name +
436 : " is not castable into a HypreParMatrix");
437 1699 : _jacobian_blocks(i, i) = mfem::ParAdd(_h_blocks(i, i), nlf_jac);
438 : }
439 : else
440 0 : _jacobian_blocks(i, i) = _h_blocks(i, i);
441 3398 : for (const auto j : index_range(_trial_var_names))
442 1699 : if (i != j) // nlf->GetGradient only contributes to on-diagonal blocks
443 0 : _jacobian_blocks(i, j) = _h_blocks(i, j);
444 1699 : }
445 : // Create monolithic matrix
446 1699 : _jacobian.Reset(mfem::HypreParMatrixFromBlocks(_jacobian_blocks));
447 1699 : }
448 :
449 : mfem::Operator &
450 3617 : EquationSystem::GetGradient(const mfem::Vector & u) const
451 : {
452 3617 : _linearization_point = &u;
453 :
454 3617 : if (_non_linear)
455 : {
456 1701 : if (_assembly_level != mfem::AssemblyLevel::LEGACY)
457 2 : mooseError("MFEM nonlinear solvers that require GetGradient() currently require legacy "
458 : "assembly in EquationSystem.");
459 1699 : const_cast<EquationSystem *>(this)->FormJacobianMatrix(u);
460 : }
461 : else
462 1916 : _jacobian = _linear_operator;
463 :
464 3615 : return *_jacobian;
465 : }
466 :
467 : void
468 2549 : EquationSystem::SetTrialVariablesFromTrueVectors(const mfem::BlockVector & trueX) const
469 : {
470 5111 : for (const auto i : index_range(_trial_var_names))
471 : {
472 2562 : auto & trial_var_name = _trial_var_names.at(i);
473 2562 : trueX.GetBlock(i).SyncMemory(trueX);
474 2562 : _gfuncs->Get(trial_var_name)->Distribute(&(trueX.GetBlock(i)));
475 : }
476 : // Solution variables changed: stored projections of solution-dependent coefficients are stale.
477 2549 : if (_coefficient_manager)
478 2549 : _coefficient_manager->markSolutionChanged();
479 2549 : }
480 :
481 : void
482 2648 : EquationSystem::BuildLinearForms()
483 : {
484 : // Register linear forms
485 5331 : for (const auto i : index_range(_test_var_names))
486 : {
487 2683 : auto test_var_name = _test_var_names.at(i);
488 2683 : _lfs.Register(test_var_name, std::make_shared<mfem::ParLinearForm>(_test_pfespaces.at(i)));
489 2683 : _lfs.GetRef(test_var_name) = 0.0;
490 2683 : }
491 :
492 5329 : for (auto & test_var_name : _test_var_names)
493 : {
494 : // Apply kernels
495 2683 : auto lf = _lfs.GetShared(test_var_name);
496 2683 : ApplyDomainLFIntegrators(test_var_name, lf, _kernels_map);
497 2683 : ApplyBoundaryLFIntegrators(test_var_name, lf, _integrated_bc_map);
498 2683 : lf->Assemble();
499 2681 : }
500 :
501 : // Apply essential boundary conditions
502 2646 : ApplyEssentialBCs();
503 :
504 : // Eliminate trivially eliminated variables by subtracting contributions from linear forms
505 2646 : EliminateCoupledVariables();
506 2646 : }
507 :
508 : void
509 918 : EquationSystem::BuildNonlinearForms()
510 : {
511 : // Register non-linear Action forms
512 1849 : for (const auto i : index_range(_test_var_names))
513 : {
514 935 : auto test_var_name = _test_var_names.at(i);
515 935 : _nlfs.Register(test_var_name, std::make_shared<mfem::ParNonlinearForm>(_test_pfespaces.at(i)));
516 : // Apply kernels
517 935 : auto nlf = _nlfs.GetShared(test_var_name);
518 935 : nlf->SetEssentialTrueDofs(_ess_tdof_lists.at(i));
519 937 : ApplyDomainNLFIntegrators(test_var_name, nlf, _kernels_map, std::nullopt);
520 935 : ApplyBoundaryNLFIntegrators(test_var_name, nlf, _integrated_bc_map, std::nullopt);
521 939 : }
522 914 : }
523 :
524 : void
525 920 : EquationSystem::BuildBilinearForms()
526 : {
527 : // Register bilinear forms
528 1857 : for (const auto i : index_range(_test_var_names))
529 : {
530 937 : auto test_var_name = _test_var_names.at(i);
531 937 : _blfs.Register(test_var_name, std::make_shared<mfem::ParBilinearForm>(_test_pfespaces.at(i)));
532 :
533 : // Apply kernels
534 937 : auto blf = _blfs.GetShared(test_var_name);
535 937 : blf->SetAssemblyLevel(_assembly_level);
536 1874 : ApplyBoundaryBLFIntegrators<mfem::ParBilinearForm>(
537 937 : test_var_name, test_var_name, blf, _integrated_bc_map);
538 1874 : ApplyDomainBLFIntegrators<mfem::ParBilinearForm>(
539 937 : test_var_name, test_var_name, blf, _kernels_map);
540 : // Assemble
541 937 : blf->Assemble();
542 937 : }
543 920 : }
544 :
545 : void
546 920 : EquationSystem::BuildMixedBilinearForms()
547 : {
548 : // Register mixed bilinear forms. Note that not all combinations may
549 : // have a kernel.
550 :
551 : // Create mblf for each test/coupled variable pair with an added kernel.
552 : // Mixed bilinear forms with coupled variables that are not trial variables are
553 : // associated with contributions from eliminated variables.
554 1857 : for (const auto i : index_range(_test_var_names))
555 : {
556 937 : auto test_var_name = _test_var_names.at(i);
557 937 : auto test_mblfs = std::make_shared<Moose::MFEM::NamedFieldsMap<mfem::ParMixedBilinearForm>>();
558 2004 : for (const auto j : index_range(_coupled_var_names))
559 : {
560 1067 : const auto & coupled_var_name = _coupled_var_names.at(j);
561 2134 : auto mblf = std::make_shared<mfem::ParMixedBilinearForm>(_coupled_pfespaces.at(j),
562 1067 : _test_pfespaces.at(i));
563 : // Register MixedBilinearForm if kernels exist for it, and assemble kernels
564 1067 : if (_kernels_map.Has(test_var_name) &&
565 2113 : _kernels_map.Get(test_var_name)->Has(coupled_var_name) &&
566 1046 : test_var_name != coupled_var_name)
567 : {
568 122 : mblf->SetAssemblyLevel(_assembly_level);
569 : // Apply all mixed kernels with this test/trial pair
570 244 : ApplyDomainBLFIntegrators<mfem::ParMixedBilinearForm>(
571 122 : coupled_var_name, test_var_name, mblf, _kernels_map);
572 : // Assemble mixed bilinear forms
573 122 : mblf->Assemble();
574 : // Register mixed bilinear forms associated with a single trial variable
575 : // for the current test variable
576 122 : test_mblfs->Register(coupled_var_name, mblf);
577 : }
578 1067 : }
579 : // Register all mixed bilinear form sets associated with a single test variable
580 937 : _mblfs.Register(test_var_name, test_mblfs);
581 937 : }
582 920 : }
583 :
584 : void
585 2648 : EquationSystem::BuildEquationSystem()
586 : {
587 2648 : BuildBilinearForms();
588 2648 : BuildMixedBilinearForms();
589 2648 : BuildLinearForms();
590 2646 : BuildNonlinearForms();
591 2642 : }
592 :
593 : void
594 2683 : EquationSystem::ApplyDomainLFIntegrators(
595 : const std::string & test_var_name,
596 : std::shared_ptr<mfem::ParLinearForm> form,
597 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map)
598 : {
599 2683 : if (kernels_map.Has(test_var_name) && kernels_map.Get(test_var_name)->Has(test_var_name))
600 : {
601 2634 : auto kernels = kernels_map.GetRef(test_var_name).GetRef(test_var_name);
602 6511 : for (auto & kernel : kernels)
603 : {
604 3877 : mfem::LinearFormIntegrator * integ = kernel->createLFIntegrator();
605 :
606 3877 : if (integ)
607 : {
608 498 : kernel->isSubdomainRestricted()
609 498 : ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
610 463 : : form->AddDomainIntegrator(std::move(integ));
611 : }
612 : }
613 2634 : }
614 2683 : }
615 :
616 : void
617 2681 : EquationSystem::ApplyDomainNLFIntegrators(
618 : const std::string & test_var_name,
619 : std::shared_ptr<mfem::ParNonlinearForm> form,
620 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map,
621 : std::optional<mfem::real_t> scale_factor)
622 : {
623 2681 : if (kernels_map.Has(test_var_name))
624 5433 : for (const auto & [trial_var_name, kernels] : kernels_map.GetRef(test_var_name))
625 6783 : for (auto & kernel : *kernels)
626 4013 : if (auto * integ = kernel->createNLIntegrator())
627 : {
628 698 : if (_gradient_required && (trial_var_name != test_var_name))
629 2 : mooseError("Support for off-diagonal MFEM nonlinear domain integrators in conjunction "
630 : "with a nonlinear solver that requires a gradient is not currently "
631 : "implemented. Kernel '",
632 2 : kernel->name(),
633 : "' contributes to test variable '",
634 : test_var_name,
635 : "' from trial variable '",
636 : trial_var_name,
637 : "'.");
638 :
639 696 : _non_linear = true;
640 696 : if (scale_factor.has_value())
641 674 : integ = new NLScaleIntegrator(integ, scale_factor.value(), true);
642 696 : kernel->isSubdomainRestricted()
643 696 : ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
644 672 : : form->AddDomainIntegrator(std::move(integ));
645 : }
646 2679 : }
647 :
648 : void
649 2683 : EquationSystem::ApplyBoundaryLFIntegrators(
650 : const std::string & test_var_name,
651 : std::shared_ptr<mfem::ParLinearForm> form,
652 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
653 : integrated_bc_map)
654 : {
655 2803 : if (integrated_bc_map.Has(test_var_name) &&
656 120 : integrated_bc_map.Get(test_var_name)->Has(test_var_name))
657 : {
658 116 : auto bcs = integrated_bc_map.GetRef(test_var_name).GetRef(test_var_name);
659 250 : for (auto & bc : bcs)
660 : {
661 134 : mfem::LinearFormIntegrator * integ = bc->createLFIntegrator();
662 :
663 134 : if (integ)
664 : {
665 92 : bc->isBoundaryRestricted()
666 92 : ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
667 13 : : form->AddBoundaryIntegrator(std::move(integ));
668 : }
669 : }
670 116 : }
671 2683 : }
672 :
673 : void
674 2679 : EquationSystem::ApplyBoundaryNLFIntegrators(
675 : const std::string & test_var_name,
676 : std::shared_ptr<mfem::ParNonlinearForm> form,
677 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
678 : integrated_bc_map,
679 : std::optional<mfem::real_t> scale_factor)
680 : {
681 2679 : if (integrated_bc_map.Has(test_var_name))
682 238 : for (const auto & [trial_var_name, bcs] : integrated_bc_map.GetRef(test_var_name))
683 256 : for (auto & bc : *bcs)
684 138 : if (auto * integ = bc->createNLIntegrator())
685 : {
686 44 : if (_gradient_required && (test_var_name != trial_var_name))
687 2 : mooseError(
688 : "Support for Off-diagonal MFEM nonlinear boundary integrators in conjunction with "
689 : "a nonlinear solver that requires a gradient is not currently "
690 : "implemented. Boundary condition '",
691 2 : bc->name(),
692 : "' contributes to test variable '",
693 : test_var_name,
694 : "' from trial variable '",
695 : trial_var_name,
696 : "'.");
697 :
698 42 : _non_linear = true;
699 42 : if (scale_factor.has_value())
700 42 : integ = new NLScaleIntegrator(integ, scale_factor.value(), true);
701 42 : bc->isBoundaryRestricted()
702 42 : ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
703 0 : : form->AddBoundaryIntegrator(std::move(integ));
704 : }
705 2677 : }
706 :
707 : const mfem::Vector &
708 0 : EquationSystem::GetLinearizationPoint() const
709 : {
710 0 : if (!_linearization_point)
711 0 : mooseError("EquationSystem::GetLinearizationPoint() called before GetGradient().");
712 0 : return *_linearization_point;
713 : }
714 :
715 : std::shared_ptr<mfem::ParBilinearForm>
716 7 : EquationSystem::BuildBilinearFormForFESpace(const std::string & var_name,
717 : mfem::ParFiniteElementSpace & fespace,
718 : mfem::AssemblyLevel assembly_level)
719 : {
720 7 : auto blf = std::make_shared<mfem::ParBilinearForm>(&fespace);
721 7 : blf->SetAssemblyLevel(assembly_level);
722 7 : ApplyBoundaryBLFIntegrators<mfem::ParBilinearForm>(var_name, var_name, blf, _integrated_bc_map);
723 7 : ApplyDomainBLFIntegrators<mfem::ParBilinearForm>(var_name, var_name, blf, _kernels_map);
724 7 : blf->Assemble();
725 7 : return blf;
726 0 : }
727 :
728 : std::shared_ptr<mfem::ParNonlinearForm>
729 0 : EquationSystem::BuildNonlinearFormForFESpace(const std::string & var_name,
730 : mfem::ParFiniteElementSpace & fespace,
731 : mfem::AssemblyLevel /*assembly_level*/)
732 : {
733 0 : auto nlf = std::make_shared<mfem::ParNonlinearForm>(&fespace);
734 0 : ApplyDomainNLFIntegrators(var_name, nlf, _kernels_map, std::nullopt);
735 0 : ApplyBoundaryNLFIntegrators(var_name, nlf, _integrated_bc_map, std::nullopt);
736 0 : return nlf;
737 0 : }
738 :
739 : mfem::Array<int> &
740 71 : EquationSystem::GetEssentialBoundaryMarkers(const std::string & var_name)
741 : {
742 71 : for (const auto i : index_range(_trial_var_names))
743 71 : if (_trial_var_names.at(i) == var_name)
744 71 : return _ess_markers.at(i);
745 :
746 0 : mooseError("No essential boundary markers found for variable '", var_name, "'.");
747 : }
748 :
749 : } // namespace Moose::MFEM
750 :
751 : #endif
|
