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 1582 : EquationSystem::~EquationSystem()
21 : {
22 1582 : DeleteHBlocks();
23 1582 : DeleteJacobianBlocks();
24 1582 : }
25 :
26 : void
27 4034 : EquationSystem::DeleteHBlocks()
28 : {
29 6515 : for (const auto i : make_range(_h_blocks.NumRows()))
30 5024 : for (const auto j : make_range(_h_blocks.NumCols()))
31 : {
32 2543 : if (_jacobian_blocks.NumRows() && _jacobian_blocks(i, j) == _h_blocks(i, j))
33 0 : _jacobian_blocks(i, j) = nullptr;
34 2543 : delete _h_blocks(i, j);
35 : }
36 4034 : _h_blocks.DeleteAll();
37 4034 : }
38 :
39 : void
40 2488 : EquationSystem::DeleteJacobianBlocks()
41 : {
42 3394 : 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 2488 : _jacobian_blocks.DeleteAll();
47 2488 : }
48 :
49 : bool
50 7940 : EquationSystem::VectorContainsName(const std::vector<std::string> & the_vector,
51 : const std::string & name) const
52 : {
53 7940 : return std::find(the_vector.begin(), the_vector.end(), name) != the_vector.end();
54 : }
55 :
56 : void
57 1884 : EquationSystem::AddCoupledVariableNameIfMissing(const std::string & coupled_var_name)
58 : {
59 1884 : if (!VectorContainsName(_coupled_var_names, coupled_var_name))
60 1201 : _coupled_var_names.push_back(coupled_var_name);
61 1884 : }
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 3445 : EquationSystem::AddTestVariableNameIfMissing(const std::string & test_var_name)
72 : {
73 3445 : if (!VectorContainsName(_test_var_names, test_var_name))
74 1107 : _test_var_names.push_back(test_var_name);
75 3445 : }
76 :
77 : void
78 1586 : EquationSystem::SetTrialVariableNames()
79 : {
80 : // If a coupled variable has an equation associated with it,
81 : // add it to the set of trial variables.
82 2689 : for (const auto & test_var_name : _test_var_names)
83 1103 : if (VectorContainsName(_coupled_var_names, test_var_name))
84 1103 : _trial_var_names.push_back(test_var_name);
85 :
86 : // Otherwise, add it to the set of eliminated variables.
87 2783 : for (const auto & coupled_var_name : _coupled_var_names)
88 1197 : if (!VectorContainsName(_test_var_names, coupled_var_name))
89 94 : _eliminated_var_names.push_back(coupled_var_name);
90 1586 : }
91 :
92 : void
93 1745 : EquationSystem::AddKernel(std::shared_ptr<MFEMKernel> kernel)
94 : {
95 1745 : const auto & trial_var_name = kernel->getTrialVariableName();
96 1745 : const auto & test_var_name = kernel->getTestVariableName();
97 1745 : AddCoupledVariableNameIfMissing(trial_var_name);
98 1745 : AddTestVariableNameIfMissing(test_var_name);
99 : // Register new kernels map if not present for the test variable
100 1745 : if (!_kernels_map.Has(test_var_name))
101 : {
102 : auto kernel_field_map =
103 1059 : std::make_shared<Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>>();
104 1059 : _kernels_map.Register(test_var_name, std::move(kernel_field_map));
105 1059 : }
106 : // Register new kernels map if not present for the test/trial variable pair
107 1745 : if (!_kernels_map.Get(test_var_name)->Has(trial_var_name))
108 : {
109 1166 : auto kernels = std::make_shared<std::vector<std::shared_ptr<MFEMKernel>>>();
110 1166 : _kernels_map.Get(test_var_name)->Register(trial_var_name, std::move(kernels));
111 1166 : }
112 1745 : _kernels_map.GetRef(test_var_name).Get(trial_var_name)->push_back(std::move(kernel));
113 1745 : }
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 1315 : EquationSystem::AddEssentialBC(std::shared_ptr<MFEMEssentialBC> bc)
140 : {
141 1315 : const auto & test_var_name = bc->getTestVariableName();
142 1315 : AddTestVariableNameIfMissing(test_var_name);
143 : // Register new essential bc map if not present for the test variable
144 1315 : if (!_essential_bc_map.Has(test_var_name))
145 : {
146 962 : auto bcs = std::make_shared<std::vector<std::shared_ptr<MFEMEssentialBC>>>();
147 962 : _essential_bc_map.Register(test_var_name, std::move(bcs));
148 962 : }
149 1315 : _essential_bc_map.GetRef(test_var_name).push_back(std::move(bc));
150 1315 : }
151 :
152 : void
153 1523 : EquationSystem::Init(Moose::MFEM::GridFunctions & gridfunctions,
154 : Moose::MFEM::ComplexGridFunctions & cmplx_gridfunctions,
155 : mfem::AssemblyLevel assembly_level)
156 : {
157 1523 : _assembly_level = assembly_level;
158 :
159 1523 : 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 1523 : SetTrialVariableNames();
165 :
166 2584 : for (auto & test_var_name : _test_var_names)
167 : {
168 1061 : 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 1061 : _test_pfespaces.push_back(gridfunctions.Get(test_var_name)->ParFESpace());
177 : }
178 :
179 2584 : for (auto & trial_var_name : _trial_var_names)
180 : {
181 1061 : 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 1061 : _var_ess_constraints.emplace_back(
190 2122 : std::make_unique<mfem::ParGridFunction>(gridfunctions.Get(trial_var_name)->ParFESpace()));
191 : }
192 :
193 : // Store pointers to FESpaces of all coupled variables
194 2678 : for (auto & coupled_var_name : _coupled_var_names)
195 1155 : _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 1792 : 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 1523 : _gfuncs = &gridfunctions;
205 1523 : }
206 :
207 : void
208 2514 : EquationSystem::ApplyEssentialBC(const std::string & var_name,
209 : mfem::ParGridFunction & trial_gf,
210 : mfem::Array<int> & global_ess_markers)
211 : {
212 2514 : if (_essential_bc_map.Has(var_name))
213 4694 : for (auto & bc : _essential_bc_map.GetRef(var_name))
214 : {
215 : // Set constrained DoFs values on essential boundaries
216 2907 : bc->ApplyBC(trial_gf);
217 : // Fetch marker array labelling essential boundaries of current BC
218 2907 : mfem::Array<int> ess_bdrs(bc->getBoundaryMarkers());
219 : // Add these boundary markers to the set of markers labelling all essential boundaries
220 14544 : for (const auto i : make_range(ess_bdrs.Size()))
221 11637 : global_ess_markers[i] |= ess_bdrs[i];
222 2907 : }
223 2514 : }
224 :
225 : void
226 2479 : EquationSystem::ApplyEssentialBCs()
227 : {
228 2479 : _ess_tdof_lists.resize(_trial_var_names.size());
229 2479 : _ess_markers.resize(_trial_var_names.size());
230 4993 : for (const auto i : index_range(_trial_var_names))
231 : {
232 2514 : const auto & trial_var_name = _trial_var_names.at(i);
233 2514 : mfem::ParGridFunction & trial_gf = *_var_ess_constraints.at(i);
234 :
235 : // Make sure we update the size, if this mesh has changed recently for instance
236 2514 : trial_gf.Update();
237 :
238 : // Initial guess for iterative solvers (initial condition or the previous time step solution)
239 2514 : trial_gf = _gfuncs->GetRef(trial_var_name);
240 :
241 2514 : _ess_markers.at(i).SetSize(trial_gf.ParFESpace()->GetParMesh()->bdr_attributes.Max(), 0);
242 : // Set strongly constrained DoFs of trial_gf on essential boundaries and add markers for all
243 : // essential boundaries to the _ess_markers array
244 2514 : ApplyEssentialBC(trial_var_name, trial_gf, _ess_markers.at(i));
245 2514 : trial_gf.ParFESpace()->GetEssentialTrueDofs(_ess_markers.at(i), _ess_tdof_lists.at(i));
246 : }
247 2479 : }
248 :
249 : void
250 2505 : EquationSystem::EliminateCoupledVariables()
251 : {
252 5045 : for (const auto & test_var_name : _test_var_names)
253 4288 : for (const auto & eliminated_var_name : _eliminated_var_names)
254 1876 : if (_mblfs.Has(test_var_name) && _mblfs.Get(test_var_name)->Has(eliminated_var_name) &&
255 128 : !VectorContainsName(_test_var_names, eliminated_var_name))
256 : {
257 92 : auto & mblf = *_mblfs.Get(test_var_name)->Get(eliminated_var_name);
258 92 : mblf.AddMult(*_eliminated_variables.Get(eliminated_var_name), *_lfs.Get(test_var_name), -1);
259 : }
260 2505 : }
261 :
262 : void
263 2511 : EquationSystem::FormLinearSystem(mfem::OperatorHandle & op,
264 : mfem::BlockVector & trueX,
265 : mfem::BlockVector & trueRHS)
266 : {
267 : mooseAssert(_test_var_names.size() == _trial_var_names.size(),
268 : "Number of test and trial variables must be the same for block matrix assembly.");
269 :
270 2511 : if (_assembly_level == mfem::AssemblyLevel::LEGACY)
271 2452 : FormSystemMatrix(op, trueX, trueRHS);
272 : else
273 59 : FormSystemOperator(op, trueX, trueRHS);
274 2511 : }
275 :
276 : void
277 59 : EquationSystem::FormSystemOperator(mfem::OperatorHandle & op,
278 : mfem::BlockVector & trueX,
279 : mfem::BlockVector & trueRHS)
280 : {
281 : mooseAssert(_test_var_names.size() == 1 && _test_var_names.size() == _trial_var_names.size(),
282 : "Non-legacy assembly is only supported for single test and trial variable systems");
283 :
284 59 : auto & test_var_name = _test_var_names.at(0);
285 59 : mfem::Vector aux_x, aux_rhs;
286 59 : mfem::OperatorPtr aux_a;
287 :
288 59 : auto blf = _blfs.Get(test_var_name);
289 59 : blf->FormLinearSystem(_ess_tdof_lists.at(0),
290 59 : *_var_ess_constraints.at(0),
291 59 : *_lfs.Get(test_var_name),
292 : aux_a,
293 : aux_x,
294 : aux_rhs,
295 : /*copy_interior=*/true);
296 :
297 59 : trueX.GetBlock(0) = aux_x;
298 59 : trueRHS.GetBlock(0) = aux_rhs;
299 59 : trueX.SyncFromBlocks();
300 59 : trueRHS.SyncFromBlocks();
301 :
302 59 : op.Reset(aux_a.Ptr());
303 59 : aux_a.SetOperatorOwner(false);
304 59 : }
305 :
306 : void
307 2410 : EquationSystem::FormSystemMatrix(mfem::OperatorHandle & op,
308 : mfem::BlockVector & trueX,
309 : mfem::BlockVector & trueRHS)
310 : {
311 : // Allocate block operator
312 2410 : DeleteHBlocks();
313 2410 : _h_blocks.SetSize(_test_var_names.size(), _trial_var_names.size());
314 2410 : _h_blocks = nullptr;
315 : // Zero out RHS and sync memory
316 2410 : trueRHS = 0.0;
317 2410 : trueRHS.SyncToBlocks();
318 :
319 4851 : for (const auto i : index_range(_test_var_names))
320 : {
321 2441 : auto test_var_name = _test_var_names.at(i);
322 :
323 4944 : for (const auto j : index_range(_trial_var_names))
324 : {
325 2503 : auto trial_var_name = _trial_var_names.at(j);
326 :
327 2503 : mfem::Vector aux_x, aux_rhs;
328 2503 : mfem::ParLinearForm aux_lf(_test_pfespaces.at(i));
329 2503 : mfem::HypreParMatrix * aux_a = new mfem::HypreParMatrix;
330 :
331 2503 : if (test_var_name == trial_var_name)
332 : {
333 : mooseAssert(i == j, "Trial and test variables must have the same ordering.");
334 2441 : auto blf = _blfs.Get(test_var_name);
335 2441 : blf->FormLinearSystem(_ess_tdof_lists.at(j),
336 2441 : *_var_ess_constraints.at(j),
337 2441 : *_lfs.Get(test_var_name),
338 : *aux_a,
339 : aux_x,
340 : aux_rhs,
341 : /*copy_interior=*/true);
342 2441 : trueX.GetBlock(j) = aux_x;
343 : }
344 62 : else if (_mblfs.Has(test_var_name) && _mblfs.Get(test_var_name)->Has(trial_var_name))
345 : {
346 62 : auto mblf = _mblfs.Get(test_var_name)->Get(trial_var_name);
347 62 : mblf->FormRectangularLinearSystem(_ess_tdof_lists.at(j),
348 62 : _ess_tdof_lists.at(i),
349 62 : *_var_ess_constraints.at(j),
350 62 : aux_lf = 0,
351 : *aux_a,
352 : aux_x,
353 : aux_rhs);
354 : }
355 : else
356 0 : continue;
357 :
358 2503 : trueRHS.GetBlock(i) += aux_rhs;
359 2503 : _h_blocks(i, j) = aux_a;
360 2503 : }
361 2441 : }
362 : // Sync memory
363 2410 : trueX.SyncFromBlocks();
364 2410 : trueRHS.SyncFromBlocks();
365 :
366 : // Create monolithic matrix
367 2410 : op.Reset(mfem::HypreParMatrixFromBlocks(_h_blocks));
368 2410 : }
369 :
370 : void
371 2513 : EquationSystem::FormSystem(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS)
372 : {
373 2513 : BuildEquationSystem();
374 2511 : height = trueX.Size();
375 2511 : width = trueRHS.Size();
376 : // Store block offsets
377 2511 : _block_true_offsets.SetSize(trueX.NumBlocks() + 1);
378 2511 : _block_true_offsets[0] = 0;
379 5053 : for (unsigned i = 0; i < _trial_var_names.size(); i++)
380 2542 : _block_true_offsets[i + 1] = trueX.BlockSize(i);
381 2511 : _block_true_offsets.PartialSum();
382 2511 : FormLinearSystem(_linear_operator, trueX, trueRHS);
383 2511 : }
384 :
385 : void
386 1584 : EquationSystem::Mult(const mfem::Vector & sol, mfem::Vector & residual) const
387 : {
388 1584 : if (_non_linear)
389 : {
390 1584 : ComputeNonlinearResidual(sol, residual);
391 1584 : _linear_operator->AddMult(sol, residual);
392 : }
393 : else
394 : {
395 0 : residual = 0.0;
396 0 : _linear_operator->Mult(sol, residual);
397 : }
398 :
399 1584 : sol.HostRead();
400 1584 : residual.HostRead();
401 1584 : }
402 :
403 : void
404 1584 : EquationSystem::ComputeNonlinearResidual(const mfem::Vector & sol, mfem::Vector & residual) const
405 : {
406 : mooseAssert(_non_linear, "Should not be calling this method if our forms are not nonlinear");
407 1584 : residual = 0.0;
408 :
409 1584 : const mfem::BlockVector block_solution(const_cast<mfem::Vector &>(sol), _block_true_offsets);
410 1584 : SetTrialVariablesFromTrueVectors(block_solution);
411 :
412 1584 : mfem::BlockVector block_residual(residual, _block_true_offsets);
413 3168 : for (unsigned int i = 0; i < _test_var_names.size(); i++)
414 : {
415 1584 : auto & test_var_name = _test_var_names.at(i);
416 1584 : auto nlf = _nlfs.GetShared(test_var_name);
417 1584 : nlf->Mult(block_solution.GetBlock(i), block_residual.GetBlock(i));
418 1584 : block_residual.GetBlock(i).SyncAliasMemory(block_residual);
419 1584 : }
420 1584 : }
421 :
422 : void
423 906 : EquationSystem::FormJacobianMatrix(const mfem::Vector & u)
424 : {
425 906 : DeleteJacobianBlocks();
426 906 : _jacobian_blocks.SetSize(_test_var_names.size(), _trial_var_names.size());
427 906 : _jacobian_blocks = nullptr;
428 :
429 906 : const mfem::BlockVector update_vector(const_cast<mfem::Vector &>(u), _block_true_offsets);
430 1812 : for (const auto i : index_range(_test_var_names))
431 : {
432 906 : auto test_var_name = _test_var_names.at(i);
433 906 : if (_nlfs.Has(test_var_name))
434 : {
435 906 : auto nlf = _nlfs.Get(test_var_name);
436 : mfem::HypreParMatrix * nlf_jac =
437 906 : dynamic_cast<mfem::HypreParMatrix *>(&nlf->GetGradient(update_vector.GetBlock(i)));
438 : mooseAssert(nlf_jac,
439 : "Jacobian contribution of nonlinear form associated with " + test_var_name +
440 : " is not castable into a HypreParMatrix");
441 906 : _jacobian_blocks(i, i) = mfem::ParAdd(_h_blocks(i, i), nlf_jac);
442 : }
443 : else
444 0 : _jacobian_blocks(i, i) = _h_blocks(i, i);
445 1812 : for (const auto j : index_range(_trial_var_names))
446 906 : if (i != j) // nlf->GetGradient only contributes to on-diagonal blocks
447 0 : _jacobian_blocks(i, j) = _h_blocks(i, j);
448 906 : }
449 : // Create monolithic matrix
450 906 : _jacobian.Reset(mfem::HypreParMatrixFromBlocks(_jacobian_blocks));
451 906 : }
452 :
453 : mfem::Operator &
454 908 : EquationSystem::GetGradient(const mfem::Vector & u) const
455 : {
456 908 : _linearization_point = &u;
457 :
458 908 : if (_non_linear)
459 : {
460 908 : if (_assembly_level != mfem::AssemblyLevel::LEGACY)
461 2 : mooseError("MFEM nonlinear solvers that require GetGradient() currently require legacy "
462 : "assembly in EquationSystem.");
463 906 : const_cast<EquationSystem *>(this)->FormJacobianMatrix(u);
464 : }
465 : else
466 0 : _jacobian = _linear_operator;
467 :
468 906 : return *_jacobian;
469 : }
470 :
471 : void
472 2427 : EquationSystem::SetTrialVariablesFromTrueVectors(const mfem::BlockVector & trueX) const
473 : {
474 4867 : for (const auto i : index_range(_trial_var_names))
475 : {
476 2440 : auto & trial_var_name = _trial_var_names.at(i);
477 2440 : trueX.GetBlock(i).SyncMemory(trueX);
478 2440 : _gfuncs->Get(trial_var_name)->Distribute(&(trueX.GetBlock(i)));
479 : }
480 : // Solution variables changed: stored projections of solution-dependent coefficients are stale.
481 2427 : if (_coefficient_manager)
482 2427 : _coefficient_manager->markSolutionChanged();
483 2427 : }
484 :
485 : void
486 2507 : EquationSystem::BuildLinearForms()
487 : {
488 : // Register linear forms
489 5049 : for (const auto i : index_range(_test_var_names))
490 : {
491 2542 : auto test_var_name = _test_var_names.at(i);
492 2542 : _lfs.Register(test_var_name, std::make_shared<mfem::ParLinearForm>(_test_pfespaces.at(i)));
493 2542 : _lfs.GetRef(test_var_name) = 0.0;
494 2542 : }
495 :
496 5047 : for (auto & test_var_name : _test_var_names)
497 : {
498 : // Apply kernels
499 2542 : auto lf = _lfs.GetShared(test_var_name);
500 2542 : ApplyDomainLFIntegrators(test_var_name, lf, _kernels_map);
501 2542 : ApplyBoundaryLFIntegrators(test_var_name, lf, _integrated_bc_map);
502 2542 : lf->Assemble();
503 2540 : }
504 :
505 : // Apply essential boundary conditions
506 2505 : ApplyEssentialBCs();
507 :
508 : // Eliminate trivially eliminated variables by subtracting contributions from linear forms
509 2505 : EliminateCoupledVariables();
510 2505 : }
511 :
512 : void
513 881 : EquationSystem::BuildNonlinearForms()
514 : {
515 : // Register non-linear Action forms
516 1775 : for (const auto i : index_range(_test_var_names))
517 : {
518 898 : auto test_var_name = _test_var_names.at(i);
519 898 : _nlfs.Register(test_var_name, std::make_shared<mfem::ParNonlinearForm>(_test_pfespaces.at(i)));
520 : // Apply kernels
521 898 : auto nlf = _nlfs.GetShared(test_var_name);
522 898 : nlf->SetEssentialTrueDofs(_ess_tdof_lists.at(i));
523 900 : ApplyDomainNLFIntegrators(test_var_name, nlf, _kernels_map, std::nullopt);
524 898 : ApplyBoundaryNLFIntegrators(test_var_name, nlf, _integrated_bc_map, std::nullopt);
525 902 : }
526 877 : }
527 :
528 : void
529 883 : EquationSystem::BuildBilinearForms()
530 : {
531 : // Register bilinear forms
532 1783 : for (const auto i : index_range(_test_var_names))
533 : {
534 900 : auto test_var_name = _test_var_names.at(i);
535 900 : _blfs.Register(test_var_name, std::make_shared<mfem::ParBilinearForm>(_test_pfespaces.at(i)));
536 :
537 : // Apply kernels
538 900 : auto blf = _blfs.GetShared(test_var_name);
539 900 : blf->SetAssemblyLevel(_assembly_level);
540 1800 : ApplyBoundaryBLFIntegrators<mfem::ParBilinearForm>(
541 900 : test_var_name, test_var_name, blf, _integrated_bc_map);
542 1800 : ApplyDomainBLFIntegrators<mfem::ParBilinearForm>(
543 900 : test_var_name, test_var_name, blf, _kernels_map);
544 : // Assemble
545 900 : blf->Assemble();
546 900 : }
547 883 : }
548 :
549 : void
550 883 : EquationSystem::BuildMixedBilinearForms()
551 : {
552 : // Register mixed bilinear forms. Note that not all combinations may
553 : // have a kernel.
554 :
555 : // Create mblf for each test/coupled variable pair with an added kernel.
556 : // Mixed bilinear forms with coupled variables that are not trial variables are
557 : // associated with contributions from eliminated variables.
558 1783 : for (const auto i : index_range(_test_var_names))
559 : {
560 900 : auto test_var_name = _test_var_names.at(i);
561 900 : auto test_mblfs = std::make_shared<Moose::MFEM::NamedFieldsMap<mfem::ParMixedBilinearForm>>();
562 1928 : for (const auto j : index_range(_coupled_var_names))
563 : {
564 1028 : const auto & coupled_var_name = _coupled_var_names.at(j);
565 2056 : auto mblf = std::make_shared<mfem::ParMixedBilinearForm>(_coupled_pfespaces.at(j),
566 1028 : _test_pfespaces.at(i));
567 : // Register MixedBilinearForm if kernels exist for it, and assemble kernels
568 1028 : if (_kernels_map.Has(test_var_name) &&
569 2035 : _kernels_map.Get(test_var_name)->Has(coupled_var_name) &&
570 1007 : test_var_name != coupled_var_name)
571 : {
572 120 : mblf->SetAssemblyLevel(_assembly_level);
573 : // Apply all mixed kernels with this test/trial pair
574 240 : ApplyDomainBLFIntegrators<mfem::ParMixedBilinearForm>(
575 120 : coupled_var_name, test_var_name, mblf, _kernels_map);
576 : // Assemble mixed bilinear forms
577 120 : mblf->Assemble();
578 : // Register mixed bilinear forms associated with a single trial variable
579 : // for the current test variable
580 120 : test_mblfs->Register(coupled_var_name, mblf);
581 : }
582 1028 : }
583 : // Register all mixed bilinear form sets associated with a single test variable
584 900 : _mblfs.Register(test_var_name, test_mblfs);
585 900 : }
586 883 : }
587 :
588 : void
589 2507 : EquationSystem::BuildEquationSystem()
590 : {
591 2507 : BuildBilinearForms();
592 2507 : BuildMixedBilinearForms();
593 2507 : BuildLinearForms();
594 2505 : BuildNonlinearForms();
595 2501 : }
596 :
597 : void
598 2542 : EquationSystem::ApplyDomainLFIntegrators(
599 : const std::string & test_var_name,
600 : std::shared_ptr<mfem::ParLinearForm> form,
601 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map)
602 : {
603 2542 : if (kernels_map.Has(test_var_name) && kernels_map.Get(test_var_name)->Has(test_var_name))
604 : {
605 2493 : auto kernels = kernels_map.GetRef(test_var_name).GetRef(test_var_name);
606 6151 : for (auto & kernel : kernels)
607 : {
608 3658 : mfem::LinearFormIntegrator * integ = kernel->createLFIntegrator();
609 :
610 3658 : if (integ)
611 : {
612 470 : kernel->isSubdomainRestricted()
613 470 : ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
614 435 : : form->AddDomainIntegrator(std::move(integ));
615 : }
616 : }
617 2493 : }
618 2542 : }
619 :
620 : void
621 2540 : EquationSystem::ApplyDomainNLFIntegrators(
622 : const std::string & test_var_name,
623 : std::shared_ptr<mfem::ParNonlinearForm> form,
624 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map,
625 : std::optional<mfem::real_t> scale_factor)
626 : {
627 2540 : if (kernels_map.Has(test_var_name))
628 5149 : for (const auto & [trial_var_name, kernels] : kernels_map.GetRef(test_var_name))
629 6419 : for (auto & kernel : *kernels)
630 3792 : if (auto * integ = kernel->createNLIntegrator())
631 : {
632 648 : if (_solver_requires_gradient && (trial_var_name != test_var_name))
633 2 : mooseError("Support for off-diagonal MFEM nonlinear domain integrators in conjunction "
634 : "with a nonlinear solver that requires a gradient is not currently "
635 : "implemented. Kernel '",
636 2 : kernel->name(),
637 : "' contributes to test variable '",
638 : test_var_name,
639 : "' from trial variable '",
640 : trial_var_name,
641 : "'.");
642 :
643 646 : _non_linear = true;
644 646 : if (scale_factor.has_value())
645 624 : integ = new NLScaleIntegrator(integ, scale_factor.value(), true);
646 646 : kernel->isSubdomainRestricted()
647 646 : ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
648 622 : : form->AddDomainIntegrator(std::move(integ));
649 : }
650 2538 : }
651 :
652 : void
653 2542 : EquationSystem::ApplyBoundaryLFIntegrators(
654 : const std::string & test_var_name,
655 : std::shared_ptr<mfem::ParLinearForm> form,
656 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
657 : integrated_bc_map)
658 : {
659 2656 : if (integrated_bc_map.Has(test_var_name) &&
660 114 : integrated_bc_map.Get(test_var_name)->Has(test_var_name))
661 : {
662 110 : auto bcs = integrated_bc_map.GetRef(test_var_name).GetRef(test_var_name);
663 238 : for (auto & bc : bcs)
664 : {
665 128 : mfem::LinearFormIntegrator * integ = bc->createLFIntegrator();
666 :
667 128 : if (integ)
668 : {
669 92 : bc->isBoundaryRestricted()
670 92 : ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
671 13 : : form->AddBoundaryIntegrator(std::move(integ));
672 : }
673 : }
674 110 : }
675 2542 : }
676 :
677 : void
678 2538 : EquationSystem::ApplyBoundaryNLFIntegrators(
679 : const std::string & test_var_name,
680 : std::shared_ptr<mfem::ParNonlinearForm> form,
681 : NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
682 : integrated_bc_map,
683 : std::optional<mfem::real_t> scale_factor)
684 : {
685 2538 : if (integrated_bc_map.Has(test_var_name))
686 226 : for (const auto & [trial_var_name, bcs] : integrated_bc_map.GetRef(test_var_name))
687 244 : for (auto & bc : *bcs)
688 132 : if (auto * integ = bc->createNLIntegrator())
689 : {
690 38 : if (_solver_requires_gradient && (test_var_name != trial_var_name))
691 2 : mooseError(
692 : "Support for Off-diagonal MFEM nonlinear boundary integrators in conjunction with "
693 : "a nonlinear solver that requires a gradient is not currently "
694 : "implemented. Boundary condition '",
695 2 : bc->name(),
696 : "' contributes to test variable '",
697 : test_var_name,
698 : "' from trial variable '",
699 : trial_var_name,
700 : "'.");
701 :
702 36 : _non_linear = true;
703 36 : if (scale_factor.has_value())
704 36 : integ = new NLScaleIntegrator(integ, scale_factor.value(), true);
705 36 : bc->isBoundaryRestricted()
706 36 : ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
707 0 : : form->AddBoundaryIntegrator(std::move(integ));
708 : }
709 2536 : }
710 :
711 : void
712 2504 : EquationSystem::PrepareLinearSolver(LinearSolverBase & solver)
713 : {
714 2504 : if (solver.IsLOR())
715 : {
716 52 : if (IsComplex())
717 0 : mooseError("LOR solve is not supported for complex equation systems.");
718 52 : if (IsMultivariate())
719 0 : mooseError("LOR solve is only supported for single-variable systems");
720 :
721 52 : const auto & test_var_name = _test_var_names.at(0);
722 52 : solver.SetupLOR(*_blfs.Get(test_var_name), GetEssentialBoundaryMarkers(test_var_name));
723 : }
724 :
725 : mooseAssert(_linear_operator.Ptr(),
726 : "If we are preparing a linear solver, we better have a linear operator");
727 2504 : solver.SetOperator(*_linear_operator);
728 2502 : }
729 :
730 : const mfem::Vector &
731 0 : EquationSystem::GetLinearizationPoint() const
732 : {
733 0 : if (!_linearization_point)
734 0 : mooseError("EquationSystem::GetLinearizationPoint() called before GetGradient().");
735 0 : return *_linearization_point;
736 : }
737 :
738 : std::shared_ptr<mfem::ParBilinearForm>
739 7 : EquationSystem::BuildBilinearFormForFESpace(const std::string & var_name,
740 : mfem::ParFiniteElementSpace & fespace,
741 : mfem::AssemblyLevel assembly_level)
742 : {
743 7 : auto blf = std::make_shared<mfem::ParBilinearForm>(&fespace);
744 7 : blf->SetAssemblyLevel(assembly_level);
745 7 : ApplyBoundaryBLFIntegrators<mfem::ParBilinearForm>(var_name, var_name, blf, _integrated_bc_map);
746 7 : ApplyDomainBLFIntegrators<mfem::ParBilinearForm>(var_name, var_name, blf, _kernels_map);
747 7 : blf->Assemble();
748 7 : return blf;
749 0 : }
750 :
751 : std::shared_ptr<mfem::ParNonlinearForm>
752 0 : EquationSystem::BuildNonlinearFormForFESpace(const std::string & var_name,
753 : mfem::ParFiniteElementSpace & fespace,
754 : mfem::AssemblyLevel /*assembly_level*/)
755 : {
756 0 : auto nlf = std::make_shared<mfem::ParNonlinearForm>(&fespace);
757 0 : ApplyDomainNLFIntegrators(var_name, nlf, _kernels_map, std::nullopt);
758 0 : ApplyBoundaryNLFIntegrators(var_name, nlf, _integrated_bc_map, std::nullopt);
759 0 : return nlf;
760 0 : }
761 :
762 : mfem::Array<int> &
763 59 : EquationSystem::GetEssentialBoundaryMarkers(const std::string & var_name)
764 : {
765 59 : for (const auto i : index_range(_trial_var_names))
766 59 : if (_trial_var_names.at(i) == var_name)
767 59 : return _ess_markers.at(i);
768 :
769 0 : mooseError("No essential boundary markers found for variable '", var_name, "'.");
770 : }
771 :
772 : } // namespace Moose::MFEM
773 :
774 : #endif
|
