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