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