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 1473 : EquationSystem::~EquationSystem()
20 : {
21 1473 : DeleteHBlocks();
22 1473 : DeleteJacobianBlocks();
23 1473 : }
24 :
25 : void
26 3572 : EquationSystem::DeleteHBlocks()
27 : {
28 5702 : 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 3572 : _h_blocks.DeleteAll();
36 3572 : }
37 :
38 : void
39 2073 : EquationSystem::DeleteJacobianBlocks()
40 : {
41 2673 : 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 2073 : _jacobian_blocks.DeleteAll();
46 2073 : }
47 :
48 : bool
49 7501 : EquationSystem::VectorContainsName(const std::vector<std::string> & the_vector,
50 : const std::string & name) const
51 : {
52 7501 : return std::find(the_vector.begin(), the_vector.end(), name) != the_vector.end();
53 : }
54 :
55 : void
56 1772 : EquationSystem::AddCoupledVariableNameIfMissing(const std::string & coupled_var_name)
57 : {
58 1772 : if (!VectorContainsName(_coupled_var_names, coupled_var_name))
59 1136 : _coupled_var_names.push_back(coupled_var_name);
60 1772 : }
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 3248 : EquationSystem::AddTestVariableNameIfMissing(const std::string & test_var_name)
71 : {
72 3248 : if (!VectorContainsName(_test_var_names, test_var_name))
73 1042 : _test_var_names.push_back(test_var_name);
74 3248 : }
75 :
76 : void
77 1473 : EquationSystem::SetTrialVariableNames()
78 : {
79 : // If a coupled variable has an equation associated with it,
80 : // add it to the set of trial variables.
81 2515 : for (const auto & test_var_name : _test_var_names)
82 1042 : if (VectorContainsName(_coupled_var_names, test_var_name))
83 1042 : _trial_var_names.push_back(test_var_name);
84 :
85 : // Otherwise, add it to the set of eliminated variables.
86 2609 : for (const auto & coupled_var_name : _coupled_var_names)
87 1136 : if (!VectorContainsName(_test_var_names, coupled_var_name))
88 94 : _eliminated_var_names.push_back(coupled_var_name);
89 1473 : }
90 :
91 : void
92 1635 : EquationSystem::AddKernel(std::shared_ptr<MFEMKernel> kernel)
93 : {
94 1635 : const auto & trial_var_name = kernel->getTrialVariableName();
95 1635 : const auto & test_var_name = kernel->getTestVariableName();
96 1635 : AddCoupledVariableNameIfMissing(trial_var_name);
97 1635 : AddTestVariableNameIfMissing(test_var_name);
98 : // Register new kernels map if not present for the test variable
99 1635 : if (!_kernels_map.Has(test_var_name))
100 : {
101 : auto kernel_field_map =
102 994 : std::make_shared<Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>>();
103 994 : _kernels_map.Register(test_var_name, std::move(kernel_field_map));
104 994 : }
105 : // Register new kernels map if not present for the test/trial variable pair
106 1635 : if (!_kernels_map.Get(test_var_name)->Has(trial_var_name))
107 : {
108 1101 : auto kernels = std::make_shared<std::vector<std::shared_ptr<MFEMKernel>>>();
109 1101 : _kernels_map.Get(test_var_name)->Register(trial_var_name, std::move(kernels));
110 1101 : }
111 1635 : _kernels_map.GetRef(test_var_name).Get(trial_var_name)->push_back(std::move(kernel));
112 1635 : }
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 1238 : EquationSystem::AddEssentialBC(std::shared_ptr<MFEMEssentialBC> bc)
139 : {
140 1238 : const auto & test_var_name = bc->getTestVariableName();
141 1238 : AddTestVariableNameIfMissing(test_var_name);
142 : // Register new essential bc map if not present for the test variable
143 1238 : if (!_essential_bc_map.Has(test_var_name))
144 : {
145 903 : auto bcs = std::make_shared<std::vector<std::shared_ptr<MFEMEssentialBC>>>();
146 903 : _essential_bc_map.Register(test_var_name, std::move(bcs));
147 903 : }
148 1238 : _essential_bc_map.GetRef(test_var_name).push_back(std::move(bc));
149 1238 : }
150 :
151 : void
152 1410 : EquationSystem::Init(Moose::MFEM::GridFunctions & gridfunctions,
153 : Moose::MFEM::ComplexGridFunctions & cmplx_gridfunctions,
154 : mfem::AssemblyLevel assembly_level)
155 : {
156 1410 : _assembly_level = assembly_level;
157 :
158 1410 : 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 1410 : SetTrialVariableNames();
164 :
165 2410 : for (auto & test_var_name : _test_var_names)
166 : {
167 1000 : 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 1000 : _test_pfespaces.push_back(gridfunctions.Get(test_var_name)->ParFESpace());
176 : }
177 :
178 2410 : for (auto & trial_var_name : _trial_var_names)
179 : {
180 1000 : 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 1000 : _var_ess_constraints.emplace_back(
189 2000 : std::make_unique<mfem::ParGridFunction>(gridfunctions.Get(trial_var_name)->ParFESpace()));
190 : }
191 :
192 : // Store pointers to FESpaces of all coupled variables
193 2504 : for (auto & coupled_var_name : _coupled_var_names)
194 1094 : _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 1673 : 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 1410 : _gfuncs = &gridfunctions;
204 1410 : }
205 :
206 : void
207 2213 : EquationSystem::ApplyEssentialBC(const std::string & var_name,
208 : mfem::ParGridFunction & trial_gf,
209 : mfem::Array<int> & global_ess_markers)
210 : {
211 2213 : if (_essential_bc_map.Has(var_name))
212 : {
213 1786 : auto & bcs = _essential_bc_map.GetRef(var_name);
214 4705 : for (auto & bc : bcs)
215 : {
216 : // Set constrained DoFs values on essential boundaries
217 2919 : bc->ApplyBC(trial_gf);
218 : // Fetch marker array labelling essential boundaries of current BC
219 2919 : mfem::Array<int> ess_bdrs(bc->getBoundaryMarkers());
220 : // Add these boundary markers to the set of markers labelling all essential boundaries
221 14864 : for (const auto i : make_range(trial_gf.ParFESpace()->GetParMesh()->bdr_attributes.Max()))
222 11945 : global_ess_markers[i] = std::max(global_ess_markers[i], ess_bdrs[i]);
223 2919 : }
224 : }
225 2213 : }
226 :
227 : void
228 2126 : EquationSystem::ApplyEssentialBCs()
229 : {
230 2126 : _ess_tdof_lists.resize(_trial_var_names.size());
231 4287 : for (const auto i : index_range(_trial_var_names))
232 : {
233 2161 : const auto & trial_var_name = _trial_var_names.at(i);
234 2161 : 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 2161 : trial_gf.Update();
238 :
239 : // Initial guess for iterative solvers (initial condition or the previous time step solution)
240 2161 : trial_gf = _gfuncs->GetRef(trial_var_name);
241 :
242 2161 : mfem::Array<int> global_ess_markers(trial_gf.ParFESpace()->GetParMesh()->bdr_attributes.Max());
243 2161 : 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 2161 : ApplyEssentialBC(trial_var_name, trial_gf, global_ess_markers);
247 2161 : trial_gf.ParFESpace()->GetEssentialTrueDofs(global_ess_markers, _ess_tdof_lists.at(i));
248 2161 : }
249 2126 : }
250 :
251 : void
252 2152 : EquationSystem::EliminateCoupledVariables()
253 : {
254 4339 : for (const auto & test_var_name : _test_var_names)
255 3635 : 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 2152 : }
263 :
264 : void
265 2158 : 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 2158 : 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 59 : FormSystemOperator(op, trueX, trueRHS);
279 : }
280 2158 : }
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 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 2158 : EquationSystem::FormSystem(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS)
375 : {
376 2158 : height = trueX.Size();
377 2158 : width = trueRHS.Size();
378 : // Store block offsets
379 2158 : _block_true_offsets.SetSize(trueX.NumBlocks() + 1);
380 2158 : _block_true_offsets[0] = 0;
381 4347 : for (unsigned i = 0; i < _trial_var_names.size(); i++)
382 2189 : _block_true_offsets[i + 1] = trueX.BlockSize(i);
383 2158 : _block_true_offsets.PartialSum();
384 2158 : FormLinearSystem(_linear_operator, trueX, trueRHS);
385 2158 : }
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 1770 : EquationSystem::SetTrialVariablesFromTrueVectors(const mfem::BlockVector & trueX) const
473 : {
474 3553 : for (const auto i : index_range(_trial_var_names))
475 : {
476 1783 : auto & trial_var_name = _trial_var_names.at(i);
477 1783 : trueX.GetBlock(i).SyncMemory(trueX);
478 1783 : _gfuncs->Get(trial_var_name)->Distribute(&(trueX.GetBlock(i)));
479 : }
480 1770 : }
481 :
482 : void
483 2152 : EquationSystem::BuildLinearForms()
484 : {
485 : // Register linear forms
486 4339 : for (const auto i : index_range(_test_var_names))
487 : {
488 2187 : auto test_var_name = _test_var_names.at(i);
489 2187 : _lfs.Register(test_var_name, std::make_shared<mfem::ParLinearForm>(_test_pfespaces.at(i)));
490 2187 : _lfs.GetRef(test_var_name) = 0.0;
491 2187 : }
492 :
493 4339 : for (auto & test_var_name : _test_var_names)
494 : {
495 : // Apply kernels
496 2187 : auto lf = _lfs.GetShared(test_var_name);
497 2187 : ApplyDomainLFIntegrators(test_var_name, lf, _kernels_map);
498 2187 : ApplyBoundaryLFIntegrators(test_var_name, lf, _integrated_bc_map);
499 2187 : lf->Assemble();
500 2187 : }
501 :
502 : // Apply essential boundary conditions
503 2152 : ApplyEssentialBCs();
504 :
505 : // Eliminate trivially eliminated variables by subtracting contributions from linear forms
506 2152 : EliminateCoupledVariables();
507 2152 : }
508 :
509 : void
510 828 : EquationSystem::BuildNonlinearForms()
511 : {
512 : // Register non-linear Action forms
513 1669 : for (const auto i : index_range(_test_var_names))
514 : {
515 845 : auto test_var_name = _test_var_names.at(i);
516 845 : _nlfs.Register(test_var_name, std::make_shared<mfem::ParNonlinearForm>(_test_pfespaces.at(i)));
517 : // Apply kernels
518 845 : auto nlf = _nlfs.GetShared(test_var_name);
519 845 : nlf->SetEssentialTrueDofs(_ess_tdof_lists.at(i));
520 847 : ApplyDomainNLFIntegrators(test_var_name, nlf, _kernels_map, std::nullopt);
521 845 : ApplyBoundaryNLFIntegrators(test_var_name, nlf, _integrated_bc_map, std::nullopt);
522 849 : }
523 824 : }
524 :
525 : void
526 828 : EquationSystem::BuildBilinearForms()
527 : {
528 : // Register bilinear forms
529 1673 : for (const auto i : index_range(_test_var_names))
530 : {
531 845 : auto test_var_name = _test_var_names.at(i);
532 845 : _blfs.Register(test_var_name, std::make_shared<mfem::ParBilinearForm>(_test_pfespaces.at(i)));
533 :
534 : // Apply kernels
535 845 : auto blf = _blfs.GetShared(test_var_name);
536 845 : blf->SetAssemblyLevel(_assembly_level);
537 1690 : ApplyBoundaryBLFIntegrators<mfem::ParBilinearForm>(
538 845 : test_var_name, test_var_name, blf, _integrated_bc_map);
539 1690 : ApplyDomainBLFIntegrators<mfem::ParBilinearForm>(
540 845 : test_var_name, test_var_name, blf, _kernels_map);
541 : // Assemble
542 845 : blf->Assemble();
543 845 : }
544 828 : }
545 :
546 : void
547 828 : 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 1673 : for (const auto i : index_range(_test_var_names))
556 : {
557 845 : auto test_var_name = _test_var_names.at(i);
558 845 : auto test_mblfs = std::make_shared<Moose::MFEM::NamedFieldsMap<mfem::ParMixedBilinearForm>>();
559 1818 : for (const auto j : index_range(_coupled_var_names))
560 : {
561 973 : const auto & coupled_var_name = _coupled_var_names.at(j);
562 1946 : auto mblf = std::make_shared<mfem::ParMixedBilinearForm>(_coupled_pfespaces.at(j),
563 973 : _test_pfespaces.at(i));
564 : // Register MixedBilinearForm if kernels exist for it, and assemble kernels
565 973 : if (_kernels_map.Has(test_var_name) &&
566 1925 : _kernels_map.Get(test_var_name)->Has(coupled_var_name) &&
567 952 : 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 973 : }
580 : // Register all mixed bilinear form sets associated with a single test variable
581 845 : _mblfs.Register(test_var_name, test_mblfs);
582 845 : }
583 828 : }
584 :
585 : void
586 2152 : EquationSystem::BuildEquationSystem()
587 : {
588 2152 : BuildBilinearForms();
589 2152 : BuildMixedBilinearForms();
590 2152 : BuildLinearForms();
591 2152 : BuildNonlinearForms();
592 2148 : }
593 :
594 : void
595 2187 : 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 2187 : if (kernels_map.Has(test_var_name) && kernels_map.Get(test_var_name)->Has(test_var_name))
601 : {
602 2138 : auto kernels = kernels_map.GetRef(test_var_name).GetRef(test_var_name);
603 5104 : for (auto & kernel : kernels)
604 : {
605 2966 : mfem::LinearFormIntegrator * integ = kernel->createLFIntegrator();
606 :
607 2966 : 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 2138 : }
615 2187 : }
616 :
617 : void
618 2187 : 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 2187 : if (kernels_map.Has(test_var_name))
625 4443 : for (const auto & [trial_var_name, kernels] : kernels_map.GetRef(test_var_name))
626 5378 : for (auto & kernel : *kernels)
627 3104 : 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 2185 : }
648 :
649 : void
650 2187 : 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 2301 : 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 2187 : }
673 :
674 : void
675 2185 : 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 2185 : 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 2183 : }
707 :
708 : void
709 2151 : EquationSystem::PrepareLinearSolver(LinearSolverBase & solver)
710 : {
711 2151 : 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 2151 : solver.SetOperator(_linear_operator);
730 2151 : }
731 :
732 : } // namespace Moose::MFEM
733 :
734 : #endif
|
