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