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 1491 : EquationSystem::~EquationSystem()
20 : {
21 1491 : DeleteHBlocks();
22 1491 : DeleteJacobianBlocks();
23 1491 : }
24 :
25 : void
26 3608 : EquationSystem::DeleteHBlocks()
27 : {
28 5756 : for (const auto i : make_range(_h_blocks.NumRows()))
29 4358 : for (const auto j : make_range(_h_blocks.NumCols()))
30 : {
31 2210 : if (_jacobian_blocks.NumRows() && _jacobian_blocks(i, j) == _h_blocks(i, j))
32 0 : _jacobian_blocks(i, j) = nullptr;
33 2210 : delete _h_blocks(i, j);
34 : }
35 3608 : _h_blocks.DeleteAll();
36 3608 : }
37 :
38 : void
39 2091 : EquationSystem::DeleteJacobianBlocks()
40 : {
41 2691 : 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 2091 : _jacobian_blocks.DeleteAll();
46 2091 : }
47 :
48 : bool
49 7623 : EquationSystem::VectorContainsName(const std::vector<std::string> & the_vector,
50 : const std::string & name) const
51 : {
52 7623 : return std::find(the_vector.begin(), the_vector.end(), name) != the_vector.end();
53 : }
54 :
55 : void
56 1794 : EquationSystem::AddCoupledVariableNameIfMissing(const std::string & coupled_var_name)
57 : {
58 1794 : if (!VectorContainsName(_coupled_var_names, coupled_var_name))
59 1156 : _coupled_var_names.push_back(coupled_var_name);
60 1794 : }
61 :
62 : void
63 177 : EquationSystem::AddEliminatedVariableNameIfMissing(const std::string & eliminated_var_name)
64 : {
65 177 : if (!VectorContainsName(_eliminated_var_names, eliminated_var_name))
66 171 : _eliminated_var_names.push_back(eliminated_var_name);
67 177 : }
68 :
69 : void
70 3310 : EquationSystem::AddTestVariableNameIfMissing(const std::string & test_var_name)
71 : {
72 3310 : if (!VectorContainsName(_test_var_names, test_var_name))
73 1062 : _test_var_names.push_back(test_var_name);
74 3310 : }
75 :
76 : void
77 1491 : EquationSystem::SetTrialVariableNames()
78 : {
79 : // If a coupled variable has an equation associated with it,
80 : // add it to the set of trial variables.
81 2551 : for (const auto & test_var_name : _test_var_names)
82 1060 : if (VectorContainsName(_coupled_var_names, test_var_name))
83 1060 : _trial_var_names.push_back(test_var_name);
84 :
85 : // Otherwise, add it to the set of eliminated variables.
86 2645 : for (const auto & coupled_var_name : _coupled_var_names)
87 1154 : if (!VectorContainsName(_test_var_names, coupled_var_name))
88 94 : _eliminated_var_names.push_back(coupled_var_name);
89 1491 : }
90 :
91 : void
92 1655 : EquationSystem::AddKernel(std::shared_ptr<MFEMKernel> kernel)
93 : {
94 1655 : const auto & trial_var_name = kernel->getTrialVariableName();
95 1655 : const auto & test_var_name = kernel->getTestVariableName();
96 1655 : AddCoupledVariableNameIfMissing(trial_var_name);
97 1655 : AddTestVariableNameIfMissing(test_var_name);
98 : // Register new kernels map if not present for the test variable
99 1655 : if (!_kernels_map.Has(test_var_name))
100 : {
101 : auto kernel_field_map =
102 1014 : std::make_shared<Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>>();
103 1014 : _kernels_map.Register(test_var_name, std::move(kernel_field_map));
104 1014 : }
105 : // Register new kernels map if not present for the test/trial variable pair
106 1655 : if (!_kernels_map.Get(test_var_name)->Has(trial_var_name))
107 : {
108 1121 : auto kernels = std::make_shared<std::vector<std::shared_ptr<MFEMKernel>>>();
109 1121 : _kernels_map.Get(test_var_name)->Register(trial_var_name, std::move(kernels));
110 1121 : }
111 1655 : _kernels_map.GetRef(test_var_name).Get(trial_var_name)->push_back(std::move(kernel));
112 1655 : }
113 :
114 : void
115 62 : EquationSystem::AddIntegratedBC(std::shared_ptr<MFEMIntegratedBC> bc)
116 : {
117 62 : const auto & trial_var_name = bc->getTrialVariableName();
118 62 : const auto & test_var_name = bc->getTestVariableName();
119 62 : AddCoupledVariableNameIfMissing(trial_var_name);
120 62 : AddTestVariableNameIfMissing(test_var_name);
121 : // Register new integrated bc map if not present for the test variable
122 62 : if (!_integrated_bc_map.Has(test_var_name))
123 : {
124 : auto integrated_bc_field_map = std::make_shared<
125 56 : Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>>();
126 56 : _integrated_bc_map.Register(test_var_name, std::move(integrated_bc_field_map));
127 56 : }
128 : // Register new integrated bc map if not present for the test/trial variable pair
129 62 : if (!_integrated_bc_map.Get(test_var_name)->Has(trial_var_name))
130 : {
131 56 : auto bcs = std::make_shared<std::vector<std::shared_ptr<MFEMIntegratedBC>>>();
132 56 : _integrated_bc_map.Get(test_var_name)->Register(trial_var_name, std::move(bcs));
133 56 : }
134 62 : _integrated_bc_map.GetRef(test_var_name).Get(trial_var_name)->push_back(std::move(bc));
135 62 : }
136 :
137 : void
138 1276 : EquationSystem::AddEssentialBC(std::shared_ptr<MFEMEssentialBC> bc)
139 : {
140 1276 : const auto & test_var_name = bc->getTestVariableName();
141 1276 : AddTestVariableNameIfMissing(test_var_name);
142 : // Register new essential bc map if not present for the test variable
143 1276 : if (!_essential_bc_map.Has(test_var_name))
144 : {
145 923 : auto bcs = std::make_shared<std::vector<std::shared_ptr<MFEMEssentialBC>>>();
146 923 : _essential_bc_map.Register(test_var_name, std::move(bcs));
147 923 : }
148 1276 : _essential_bc_map.GetRef(test_var_name).push_back(std::move(bc));
149 1276 : }
150 :
151 : void
152 1428 : EquationSystem::Init(Moose::MFEM::GridFunctions & gridfunctions,
153 : Moose::MFEM::ComplexGridFunctions & cmplx_gridfunctions,
154 : mfem::AssemblyLevel assembly_level)
155 : {
156 1428 : _assembly_level = assembly_level;
157 :
158 1428 : 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 1428 : SetTrialVariableNames();
164 :
165 2446 : for (auto & test_var_name : _test_var_names)
166 : {
167 1018 : 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 1018 : _test_pfespaces.push_back(gridfunctions.Get(test_var_name)->ParFESpace());
176 : }
177 :
178 2446 : for (auto & trial_var_name : _trial_var_names)
179 : {
180 1018 : 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 1018 : _var_ess_constraints.emplace_back(
189 2036 : std::make_unique<mfem::ParGridFunction>(gridfunctions.Get(trial_var_name)->ParFESpace()));
190 : }
191 :
192 : // Store pointers to FESpaces of all coupled variables
193 2540 : for (auto & coupled_var_name : _coupled_var_names)
194 1112 : _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 1691 : 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 1428 : _gfuncs = &gridfunctions;
204 1428 : }
205 :
206 : void
207 2231 : EquationSystem::ApplyEssentialBC(const std::string & var_name,
208 : mfem::ParGridFunction & trial_gf,
209 : mfem::Array<int> & global_ess_markers)
210 : {
211 2231 : if (_essential_bc_map.Has(var_name))
212 : {
213 1804 : auto & bcs = _essential_bc_map.GetRef(var_name);
214 4759 : for (auto & bc : bcs)
215 : {
216 : // Set constrained DoFs values on essential boundaries
217 2955 : bc->ApplyBC(trial_gf);
218 : // Fetch marker array labelling essential boundaries of current BC
219 2955 : mfem::Array<int> ess_bdrs(bc->getBoundaryMarkers());
220 : // Add these boundary markers to the set of markers labelling all essential boundaries
221 14972 : for (const auto i : make_range(trial_gf.ParFESpace()->GetParMesh()->bdr_attributes.Max()))
222 12017 : global_ess_markers[i] = std::max(global_ess_markers[i], ess_bdrs[i]);
223 2955 : }
224 : }
225 2231 : }
226 :
227 : void
228 2144 : EquationSystem::ApplyEssentialBCs()
229 : {
230 2144 : _ess_tdof_lists.resize(_trial_var_names.size());
231 4323 : for (const auto i : index_range(_trial_var_names))
232 : {
233 2179 : const auto & trial_var_name = _trial_var_names.at(i);
234 2179 : 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 2179 : trial_gf.Update();
238 :
239 : // Initial guess for iterative solvers (initial condition or the previous time step solution)
240 2179 : trial_gf = _gfuncs->GetRef(trial_var_name);
241 :
242 2179 : mfem::Array<int> global_ess_markers(trial_gf.ParFESpace()->GetParMesh()->bdr_attributes.Max());
243 2179 : 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 2179 : ApplyEssentialBC(trial_var_name, trial_gf, global_ess_markers);
247 2179 : trial_gf.ParFESpace()->GetEssentialTrueDofs(global_ess_markers, _ess_tdof_lists.at(i));
248 2179 : }
249 2144 : }
250 :
251 : void
252 2170 : EquationSystem::EliminateCoupledVariables()
253 : {
254 4375 : for (const auto & test_var_name : _test_var_names)
255 3653 : 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 2170 : }
263 :
264 : void
265 2176 : 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 2176 : if (_assembly_level == mfem::AssemblyLevel::LEGACY)
273 2117 : 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 59 : FormSystemOperator(op, trueX, trueRHS);
279 : }
280 2176 : }
281 :
282 : void
283 59 : EquationSystem::FormSystemOperator(mfem::OperatorHandle & op,
284 : mfem::BlockVector & trueX,
285 : mfem::BlockVector & trueRHS)
286 : {
287 59 : auto & test_var_name = _test_var_names.at(0);
288 59 : mfem::Vector aux_x, aux_rhs;
289 59 : mfem::OperatorPtr aux_a;
290 :
291 59 : auto blf = _blfs.Get(test_var_name);
292 59 : blf->FormLinearSystem(_ess_tdof_lists.at(0),
293 59 : *_var_ess_constraints.at(0),
294 59 : *_lfs.Get(test_var_name),
295 : aux_a,
296 : aux_x,
297 : aux_rhs,
298 : /*copy_interior=*/true);
299 :
300 59 : trueX.GetBlock(0) = aux_x;
301 59 : trueRHS.GetBlock(0) = aux_rhs;
302 59 : trueX.SyncFromBlocks();
303 59 : trueRHS.SyncFromBlocks();
304 :
305 59 : op.Reset(aux_a.Ptr());
306 59 : aux_a.SetOperatorOwner(false);
307 59 : }
308 :
309 : void
310 2075 : EquationSystem::FormSystemMatrix(mfem::OperatorHandle & op,
311 : mfem::BlockVector & trueX,
312 : mfem::BlockVector & trueRHS)
313 : {
314 : // Allocate block operator
315 2075 : DeleteHBlocks();
316 2075 : _h_blocks.SetSize(_test_var_names.size(), _trial_var_names.size());
317 2075 : _h_blocks = nullptr;
318 : // Zero out RHS and sync memory
319 2075 : trueRHS = 0.0;
320 2075 : trueRHS.SyncToBlocks();
321 :
322 4181 : for (const auto i : index_range(_test_var_names))
323 : {
324 2106 : auto test_var_name = _test_var_names.at(i);
325 :
326 4274 : for (const auto j : index_range(_trial_var_names))
327 : {
328 2168 : auto trial_var_name = _trial_var_names.at(j);
329 :
330 2168 : mfem::Vector aux_x, aux_rhs;
331 2168 : mfem::ParLinearForm aux_lf(_test_pfespaces.at(i));
332 2168 : mfem::HypreParMatrix * aux_a = new mfem::HypreParMatrix;
333 :
334 2168 : if (test_var_name == trial_var_name)
335 : {
336 : mooseAssert(i == j, "Trial and test variables must have the same ordering.");
337 2106 : auto blf = _blfs.Get(test_var_name);
338 2106 : blf->FormLinearSystem(_ess_tdof_lists.at(j),
339 2106 : *_var_ess_constraints.at(j),
340 2106 : *_lfs.Get(test_var_name),
341 : *aux_a,
342 : aux_x,
343 : aux_rhs,
344 : /*copy_interior=*/true);
345 2106 : 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 2168 : trueRHS.GetBlock(i) += aux_rhs;
362 2168 : _h_blocks(i, j) = aux_a;
363 2168 : }
364 2106 : }
365 : // Sync memory
366 2075 : trueX.SyncFromBlocks();
367 2075 : trueRHS.SyncFromBlocks();
368 :
369 : // Create monolithic matrix
370 2075 : op.Reset(mfem::HypreParMatrixFromBlocks(_h_blocks));
371 2075 : }
372 :
373 : void
374 2176 : EquationSystem::FormSystem(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS)
375 : {
376 2176 : height = trueX.Size();
377 2176 : width = trueRHS.Size();
378 : // Store block offsets
379 2176 : _block_true_offsets.SetSize(trueX.NumBlocks() + 1);
380 2176 : _block_true_offsets[0] = 0;
381 4383 : for (unsigned i = 0; i < _trial_var_names.size(); i++)
382 2207 : _block_true_offsets[i + 1] = trueX.BlockSize(i);
383 2176 : _block_true_offsets.PartialSum();
384 2176 : FormLinearSystem(_linear_operator, trueX, trueRHS);
385 2176 : }
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 1788 : EquationSystem::SetTrialVariablesFromTrueVectors(const mfem::BlockVector & trueX) const
473 : {
474 3589 : for (const auto i : index_range(_trial_var_names))
475 : {
476 1801 : auto & trial_var_name = _trial_var_names.at(i);
477 1801 : trueX.GetBlock(i).SyncMemory(trueX);
478 1801 : _gfuncs->Get(trial_var_name)->Distribute(&(trueX.GetBlock(i)));
479 : }
480 1788 : }
481 :
482 : void
483 2170 : EquationSystem::BuildLinearForms()
484 : {
485 : // Register linear forms
486 4375 : for (const auto i : index_range(_test_var_names))
487 : {
488 2205 : auto test_var_name = _test_var_names.at(i);
489 2205 : _lfs.Register(test_var_name, std::make_shared<mfem::ParLinearForm>(_test_pfespaces.at(i)));
490 2205 : _lfs.GetRef(test_var_name) = 0.0;
491 2205 : }
492 :
493 4375 : for (auto & test_var_name : _test_var_names)
494 : {
495 : // Apply kernels
496 2205 : auto lf = _lfs.GetShared(test_var_name);
497 2205 : ApplyDomainLFIntegrators(test_var_name, lf, _kernels_map);
498 2205 : ApplyBoundaryLFIntegrators(test_var_name, lf, _integrated_bc_map);
499 2205 : lf->Assemble();
500 2205 : }
501 :
502 : // Apply essential boundary conditions
503 2170 : ApplyEssentialBCs();
504 :
505 : // Eliminate trivially eliminated variables by subtracting contributions from linear forms
506 2170 : EliminateCoupledVariables();
507 2170 : }
508 :
509 : void
510 846 : EquationSystem::BuildNonlinearForms()
511 : {
512 : // Register non-linear Action forms
513 1705 : for (const auto i : index_range(_test_var_names))
514 : {
515 863 : auto test_var_name = _test_var_names.at(i);
516 863 : _nlfs.Register(test_var_name, std::make_shared<mfem::ParNonlinearForm>(_test_pfespaces.at(i)));
517 : // Apply kernels
518 863 : auto nlf = _nlfs.GetShared(test_var_name);
519 863 : nlf->SetEssentialTrueDofs(_ess_tdof_lists.at(i));
520 865 : ApplyDomainNLFIntegrators(test_var_name, nlf, _kernels_map, std::nullopt);
521 863 : ApplyBoundaryNLFIntegrators(test_var_name, nlf, _integrated_bc_map, std::nullopt);
522 867 : }
523 842 : }
524 :
525 : void
526 846 : EquationSystem::BuildBilinearForms()
527 : {
528 : // Register bilinear forms
529 1709 : for (const auto i : index_range(_test_var_names))
530 : {
531 863 : auto test_var_name = _test_var_names.at(i);
532 863 : _blfs.Register(test_var_name, std::make_shared<mfem::ParBilinearForm>(_test_pfespaces.at(i)));
533 :
534 : // Apply kernels
535 863 : auto blf = _blfs.GetShared(test_var_name);
536 863 : blf->SetAssemblyLevel(_assembly_level);
537 1726 : ApplyBoundaryBLFIntegrators<mfem::ParBilinearForm>(
538 863 : test_var_name, test_var_name, blf, _integrated_bc_map);
539 1726 : ApplyDomainBLFIntegrators<mfem::ParBilinearForm>(
540 863 : test_var_name, test_var_name, blf, _kernels_map);
541 : // Assemble
542 863 : blf->Assemble();
543 863 : }
544 846 : }
545 :
546 : void
547 846 : 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 1709 : for (const auto i : index_range(_test_var_names))
556 : {
557 863 : auto test_var_name = _test_var_names.at(i);
558 863 : auto test_mblfs = std::make_shared<Moose::MFEM::NamedFieldsMap<mfem::ParMixedBilinearForm>>();
559 1854 : for (const auto j : index_range(_coupled_var_names))
560 : {
561 991 : const auto & coupled_var_name = _coupled_var_names.at(j);
562 1982 : auto mblf = std::make_shared<mfem::ParMixedBilinearForm>(_coupled_pfespaces.at(j),
563 991 : _test_pfespaces.at(i));
564 : // Register MixedBilinearForm if kernels exist for it, and assemble kernels
565 991 : if (_kernels_map.Has(test_var_name) &&
566 1961 : _kernels_map.Get(test_var_name)->Has(coupled_var_name) &&
567 970 : 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 991 : }
580 : // Register all mixed bilinear form sets associated with a single test variable
581 863 : _mblfs.Register(test_var_name, test_mblfs);
582 863 : }
583 846 : }
584 :
585 : void
586 2170 : EquationSystem::BuildEquationSystem()
587 : {
588 2170 : BuildBilinearForms();
589 2170 : BuildMixedBilinearForms();
590 2170 : BuildLinearForms();
591 2170 : BuildNonlinearForms();
592 2166 : }
593 :
594 : void
595 2205 : 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 2205 : if (kernels_map.Has(test_var_name) && kernels_map.Get(test_var_name)->Has(test_var_name))
601 : {
602 2156 : auto kernels = kernels_map.GetRef(test_var_name).GetRef(test_var_name);
603 5140 : for (auto & kernel : kernels)
604 : {
605 2984 : mfem::LinearFormIntegrator * integ = kernel->createLFIntegrator();
606 :
607 2984 : if (integ)
608 : {
609 433 : kernel->isSubdomainRestricted()
610 433 : ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
611 413 : : form->AddDomainIntegrator(std::move(integ));
612 : }
613 : }
614 2156 : }
615 2205 : }
616 :
617 : void
618 2205 : 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 2205 : if (kernels_map.Has(test_var_name))
625 4479 : for (const auto & [trial_var_name, kernels] : kernels_map.GetRef(test_var_name))
626 5414 : for (auto & kernel : *kernels)
627 3122 : 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 2203 : }
648 :
649 : void
650 2205 : 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 2319 : 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 2205 : }
673 :
674 : void
675 2203 : 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 2203 : 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 2201 : }
707 :
708 : void
709 2169 : EquationSystem::PrepareLinearSolver(LinearSolverBase & solver)
710 : {
711 2169 : if (solver.IsLOR())
712 : {
713 52 : if (Complex())
714 0 : mooseError("LOR solve is not supported for complex equation systems.");
715 52 : if (_test_var_names.size() > 1)
716 0 : mooseError("LOR solve is only supported for single-variable systems");
717 :
718 52 : const auto & test_var_name = _test_var_names.at(0);
719 52 : const auto & trial_var_name = _trial_var_names.at(0);
720 52 : mfem::ParGridFunction & trial_gf = _gfuncs->GetRef(trial_var_name);
721 52 : mfem::Array<int> global_ess_markers(trial_gf.ParFESpace()->GetParMesh()->bdr_attributes.Max());
722 52 : global_ess_markers = 0;
723 52 : ApplyEssentialBC(trial_var_name, trial_gf, global_ess_markers);
724 52 : solver.SetupLOR(*_blfs.Get(test_var_name), global_ess_markers);
725 52 : }
726 :
727 : mooseAssert(_linear_operator.Ptr(),
728 : "If we are preparing a linear solver, we better have a linear operator");
729 2169 : solver.SetOperator(_linear_operator);
730 2169 : }
731 :
732 : } // namespace Moose::MFEM
733 :
734 : #endif
|
