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 "libmesh/int_range.h"
14 :
15 : namespace Moose::MFEM
16 : {
17 :
18 353 : EquationSystem::~EquationSystem() { DeleteAllBlocks(); }
19 :
20 : void
21 1153 : EquationSystem::DeleteAllBlocks()
22 : {
23 1973 : for (const auto i : make_range(_h_blocks.NumRows()))
24 1680 : for (const auto j : make_range(_h_blocks.NumCols()))
25 860 : delete _h_blocks(i, j);
26 1153 : _h_blocks.DeleteAll();
27 1153 : }
28 :
29 : bool
30 1936 : EquationSystem::VectorContainsName(const std::vector<std::string> & the_vector,
31 : const std::string & name) const
32 : {
33 1936 : auto iter = std::find(the_vector.begin(), the_vector.end(), name);
34 1936 : return (iter != the_vector.end());
35 : }
36 :
37 : void
38 844 : EquationSystem::AddCoupledVariableNameIfMissing(const std::string & coupled_var_name)
39 : {
40 844 : if (!VectorContainsName(_coupled_var_names, coupled_var_name))
41 474 : _coupled_var_names.push_back(coupled_var_name);
42 844 : }
43 :
44 : void
45 602 : EquationSystem::AddTestVariableNameIfMissing(const std::string & test_var_name)
46 : {
47 602 : if (!VectorContainsName(_test_var_names, test_var_name))
48 317 : _test_var_names.push_back(test_var_name);
49 602 : }
50 :
51 : void
52 236 : EquationSystem::SetTrialVariableNames()
53 : {
54 : // If a coupled variable has an equation associated with it,
55 : // add it to the set of trial variables.
56 468 : for (const auto & coupled_var_name : _coupled_var_names)
57 : {
58 232 : if (VectorContainsName(_test_var_names, coupled_var_name))
59 196 : _trial_var_names.push_back(coupled_var_name);
60 : else
61 36 : _eliminated_var_names.push_back(coupled_var_name);
62 : }
63 236 : }
64 :
65 : void
66 446 : EquationSystem::AddKernel(std::shared_ptr<MFEMKernel> kernel)
67 : {
68 446 : AddTestVariableNameIfMissing(kernel->getTestVariableName());
69 446 : AddCoupledVariableNameIfMissing(kernel->getTrialVariableName());
70 446 : auto trial_var_name = kernel->getTrialVariableName();
71 446 : auto test_var_name = kernel->getTestVariableName();
72 446 : if (!_kernels_map.Has(test_var_name))
73 : {
74 : auto kernel_field_map =
75 313 : std::make_shared<Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>>();
76 313 : _kernels_map.Register(test_var_name, std::move(kernel_field_map));
77 313 : }
78 : // Register new kernels map if not present for the test/trial variable
79 : // pair
80 446 : if (!_kernels_map.Get(test_var_name)->Has(trial_var_name))
81 : {
82 357 : auto kernels = std::make_shared<std::vector<std::shared_ptr<MFEMKernel>>>();
83 357 : _kernels_map.Get(test_var_name)->Register(trial_var_name, std::move(kernels));
84 357 : }
85 446 : _kernels_map.GetRef(test_var_name).Get(trial_var_name)->push_back(std::move(kernel));
86 446 : }
87 :
88 : void
89 38 : EquationSystem::AddIntegratedBC(std::shared_ptr<MFEMIntegratedBC> bc)
90 : {
91 38 : AddTestVariableNameIfMissing(bc->getTestVariableName());
92 38 : AddCoupledVariableNameIfMissing(bc->getTrialVariableName());
93 38 : auto trial_var_name = bc->getTrialVariableName();
94 38 : auto test_var_name = bc->getTestVariableName();
95 38 : if (!_integrated_bc_map.Has(test_var_name))
96 : {
97 : auto integrated_bc_field_map = std::make_shared<
98 34 : Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>>();
99 34 : _integrated_bc_map.Register(test_var_name, std::move(integrated_bc_field_map));
100 34 : }
101 : // Register new integrated bc map if not present for the test/trial variable
102 : // pair
103 38 : if (!_integrated_bc_map.Get(test_var_name)->Has(trial_var_name))
104 : {
105 34 : auto bcs = std::make_shared<std::vector<std::shared_ptr<MFEMIntegratedBC>>>();
106 34 : _integrated_bc_map.Get(test_var_name)->Register(trial_var_name, std::move(bcs));
107 34 : }
108 38 : _integrated_bc_map.GetRef(test_var_name).Get(trial_var_name)->push_back(std::move(bc));
109 38 : }
110 :
111 : void
112 450 : EquationSystem::AddEssentialBC(std::shared_ptr<MFEMEssentialBC> bc)
113 : {
114 450 : auto test_var_name = bc->getTestVariableName();
115 450 : if (!_essential_bc_map.Has(test_var_name))
116 : {
117 281 : auto bcs = std::make_shared<std::vector<std::shared_ptr<MFEMEssentialBC>>>();
118 281 : _essential_bc_map.Register(test_var_name, std::move(bcs));
119 281 : }
120 450 : _essential_bc_map.GetRef(test_var_name).push_back(std::move(bc));
121 450 : }
122 :
123 : void
124 353 : EquationSystem::Init(Moose::MFEM::GridFunctions & gridfunctions, mfem::AssemblyLevel assembly_level)
125 : {
126 353 : _assembly_level = assembly_level;
127 :
128 670 : for (auto & test_var_name : _test_var_names)
129 : {
130 317 : if (!gridfunctions.Has(test_var_name))
131 : {
132 0 : mooseError("MFEM variable ",
133 : test_var_name,
134 : " requested by equation system during initialisation was "
135 : "not found in gridfunctions");
136 : }
137 : // Store pointers to test FESpaces
138 317 : _test_pfespaces.push_back(gridfunctions.Get(test_var_name)->ParFESpace());
139 : // Create auxiliary gridfunctions for storing essential constraints from Dirichlet conditions
140 317 : _var_ess_constraints.emplace_back(
141 634 : std::make_unique<mfem::ParGridFunction>(gridfunctions.Get(test_var_name)->ParFESpace()));
142 : }
143 :
144 : // Extract which coupled variables are to be trivially eliminated and which are trial variables
145 353 : SetTrialVariableNames();
146 :
147 : // Store pointers to FESpaces of all coupled variables
148 827 : for (auto & coupled_var_name : _coupled_var_names)
149 474 : _coupled_pfespaces.push_back(gridfunctions.Get(coupled_var_name)->ParFESpace());
150 :
151 : // Store pointers to coupled variable GridFunctions that are to be eliminated prior to forming the
152 : // jacobian
153 518 : for (auto & eliminated_var_name : _eliminated_var_names)
154 165 : _eliminated_variables.Register(eliminated_var_name,
155 330 : gridfunctions.GetShared(eliminated_var_name));
156 353 : }
157 :
158 : void
159 1768 : EquationSystem::ApplyEssentialBC(const std::string & var_name,
160 : mfem::ParGridFunction & trial_gf,
161 : mfem::Array<int> & global_ess_markers)
162 : {
163 1768 : if (_essential_bc_map.Has(var_name))
164 : {
165 934 : auto & bcs = _essential_bc_map.GetRef(var_name);
166 2608 : for (auto & bc : bcs)
167 : {
168 : // Set constrained DoFs values on essential boundaries
169 1674 : bc->ApplyBC(trial_gf);
170 : // Fetch marker array labelling essential boundaries of current BC
171 1674 : mfem::Array<int> ess_bdrs(bc->getBoundaryMarkers());
172 : // Add these boundary markers to the set of markers labelling all essential boundaries
173 7892 : for (auto it = 0; it != trial_gf.ParFESpace()->GetParMesh()->bdr_attributes.Max(); ++it)
174 6218 : global_ess_markers[it] = std::max(global_ess_markers[it], ess_bdrs[it]);
175 1674 : }
176 : }
177 1768 : }
178 :
179 : void
180 236 : EquationSystem::ApplyEssentialBCs()
181 : {
182 236 : _ess_tdof_lists.resize(_test_var_names.size());
183 432 : for (const auto i : index_range(_test_var_names))
184 : {
185 196 : const auto & test_var_name = _test_var_names.at(i);
186 196 : mfem::ParGridFunction & trial_gf = *(_var_ess_constraints.at(i));
187 196 : mfem::Array<int> global_ess_markers(trial_gf.ParFESpace()->GetParMesh()->bdr_attributes.Max());
188 196 : global_ess_markers = 0;
189 : // Set strongly constrained DoFs of trial_gf on essential boundaries and add markers for all
190 : // essential boundaries to the global_ess_markers array
191 196 : ApplyEssentialBC(test_var_name, trial_gf, global_ess_markers);
192 196 : trial_gf.ParFESpace()->GetEssentialTrueDofs(global_ess_markers, _ess_tdof_lists.at(i));
193 196 : }
194 236 : }
195 :
196 : void
197 1006 : EquationSystem::EliminateCoupledVariables()
198 : {
199 1988 : for (const auto & test_var_name : _test_var_names)
200 : {
201 982 : auto lf = _lfs.Get(test_var_name);
202 1900 : for (const auto & eliminated_var_name : _eliminated_var_names)
203 : {
204 918 : if (_mblfs.Has(test_var_name) && _mblfs.Get(test_var_name)->Has(eliminated_var_name))
205 : {
206 52 : auto mblf = _mblfs.Get(test_var_name)->Get(eliminated_var_name);
207 : // The AddMult method in mfem::BilinearForm is not defined for non-legacy assembly
208 52 : mfem::Vector lf_prev(lf->Size());
209 52 : mblf->Mult(*_eliminated_variables.Get(eliminated_var_name), lf_prev);
210 52 : *lf -= lf_prev;
211 52 : }
212 : }
213 : }
214 1006 : }
215 :
216 : void
217 841 : EquationSystem::FormLinearSystem(mfem::OperatorHandle & op,
218 : mfem::BlockVector & trueX,
219 : mfem::BlockVector & trueRHS)
220 : {
221 :
222 841 : switch (_assembly_level)
223 : {
224 800 : case mfem::AssemblyLevel::LEGACY:
225 800 : FormLegacySystem(op, trueX, trueRHS);
226 800 : break;
227 41 : default:
228 : mooseAssert(_test_var_names.size() == 1,
229 : "Non-legacy assembly is only supported for single-variable systems");
230 : mooseAssert(
231 : _test_var_names.size() == _trial_var_names.size(),
232 : "Non-legacy assembly is only supported for single test and trial variable systems");
233 41 : FormSystem(op, trueX, trueRHS);
234 : }
235 841 : }
236 :
237 : void
238 9 : EquationSystem::FormSystem(mfem::OperatorHandle & op,
239 : mfem::BlockVector & trueX,
240 : mfem::BlockVector & trueRHS)
241 : {
242 9 : auto & test_var_name = _test_var_names.at(0);
243 9 : auto blf = _blfs.Get(test_var_name);
244 9 : auto lf = _lfs.Get(test_var_name);
245 9 : mfem::BlockVector aux_x, aux_rhs;
246 9 : mfem::OperatorPtr aux_a;
247 :
248 27 : blf->FormLinearSystem(
249 9 : _ess_tdof_lists.at(0), *(_var_ess_constraints.at(0)), *lf, aux_a, aux_x, aux_rhs);
250 :
251 9 : trueX.GetBlock(0) = aux_x;
252 9 : trueRHS.GetBlock(0) = aux_rhs;
253 9 : trueX.SyncFromBlocks();
254 9 : trueRHS.SyncFromBlocks();
255 :
256 9 : op.Reset(aux_a.Ptr());
257 9 : aux_a.SetOperatorOwner(false);
258 9 : }
259 :
260 : void
261 800 : EquationSystem::AssembleJacobian(
262 : Moose::MFEM::NamedFieldsMap<mfem::ParBilinearForm> & jac_blfs,
263 : Moose::MFEM::NamedFieldsMap<Moose::MFEM::NamedFieldsMap<mfem::ParMixedBilinearForm>> &
264 : jac_mblfs,
265 : Moose::MFEM::NamedFieldsMap<mfem::ParLinearForm> & rhs_lfs,
266 : std::vector<mfem::Array<int>> & ess_tdof_lists,
267 : std::vector<std::unique_ptr<mfem::ParGridFunction>> & var_ess_constraints,
268 : mfem::OperatorHandle & op,
269 : mfem::BlockVector & trueX,
270 : mfem::BlockVector & trueRHS)
271 : {
272 : // Allocate block operator
273 800 : DeleteAllBlocks();
274 800 : _h_blocks.SetSize(_test_var_names.size(), _test_var_names.size());
275 : // Form diagonal blocks.
276 1620 : for (const auto i : index_range(_test_var_names))
277 : {
278 820 : auto & test_var_name = _test_var_names.at(i);
279 820 : auto blf = jac_blfs.Get(test_var_name);
280 820 : auto lf = rhs_lfs.Get(test_var_name);
281 820 : mfem::Vector aux_x, aux_rhs;
282 820 : mfem::HypreParMatrix * aux_a = new mfem::HypreParMatrix;
283 2460 : blf->FormLinearSystem(
284 820 : ess_tdof_lists.at(i), *(var_ess_constraints.at(i)), *lf, *aux_a, aux_x, aux_rhs);
285 820 : _h_blocks(i, i) = aux_a;
286 820 : trueX.GetBlock(i) = aux_x;
287 820 : trueRHS.GetBlock(i) = aux_rhs;
288 820 : }
289 :
290 : // Form off-diagonal blocks
291 1620 : for (const auto i : index_range(_test_var_names))
292 : {
293 820 : auto test_var_name = _test_var_names.at(i);
294 1680 : for (const auto j : index_range(_trial_var_names))
295 : {
296 860 : auto trial_var_name = _trial_var_names.at(j);
297 :
298 860 : mfem::Vector aux_x, aux_rhs;
299 860 : mfem::ParLinearForm aux_lf(_test_pfespaces.at(i));
300 860 : aux_lf = 0.0;
301 860 : if (jac_mblfs.Has(test_var_name) && jac_mblfs.Get(test_var_name)->Has(trial_var_name))
302 : {
303 40 : auto mblf = jac_mblfs.Get(test_var_name)->Get(trial_var_name);
304 40 : mfem::HypreParMatrix * aux_a = new mfem::HypreParMatrix;
305 40 : mblf->FormRectangularLinearSystem(ess_tdof_lists.at(j),
306 40 : ess_tdof_lists.at(i),
307 40 : *(var_ess_constraints.at(j)),
308 : aux_lf,
309 : *aux_a,
310 : aux_x,
311 : aux_rhs);
312 40 : _h_blocks(i, j) = aux_a;
313 40 : trueRHS.GetBlock(i) += aux_rhs;
314 : }
315 860 : }
316 820 : }
317 : // Sync memory
318 800 : trueX.SyncFromBlocks();
319 800 : trueRHS.SyncFromBlocks();
320 :
321 : // Create monolithic matrix
322 800 : op.Reset(mfem::HypreParMatrixFromBlocks(_h_blocks));
323 800 : }
324 :
325 : void
326 179 : EquationSystem::FormLegacySystem(mfem::OperatorHandle & op,
327 : mfem::BlockVector & trueX,
328 : mfem::BlockVector & trueRHS)
329 : {
330 179 : AssembleJacobian(_blfs, _mblfs, _lfs, _ess_tdof_lists, _var_ess_constraints, op, trueX, trueRHS);
331 179 : }
332 :
333 : void
334 841 : EquationSystem::BuildJacobian(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS)
335 : {
336 841 : height = trueX.Size();
337 841 : width = trueRHS.Size();
338 841 : FormLinearSystem(_jacobian, trueX, trueRHS);
339 841 : }
340 :
341 : void
342 1682 : EquationSystem::Mult(const mfem::Vector & x, mfem::Vector & residual) const
343 : {
344 1682 : _jacobian->Mult(x, residual);
345 1682 : x.HostRead();
346 1682 : residual.HostRead();
347 1682 : }
348 :
349 : mfem::Operator &
350 841 : EquationSystem::GetGradient(const mfem::Vector &) const
351 : {
352 841 : return *_jacobian;
353 : }
354 :
355 : void
356 188 : EquationSystem::RecoverFEMSolution(mfem::BlockVector & trueX,
357 : Moose::MFEM::GridFunctions & gridfunctions)
358 : {
359 384 : for (const auto i : index_range(_trial_var_names))
360 : {
361 196 : auto & trial_var_name = _trial_var_names.at(i);
362 196 : trueX.GetBlock(i).SyncAliasMemory(trueX);
363 196 : gridfunctions.Get(trial_var_name)->Distribute(&(trueX.GetBlock(i)));
364 : }
365 188 : }
366 :
367 : void
368 1006 : EquationSystem::BuildLinearForms()
369 : {
370 : // Register linear forms
371 1988 : for (const auto i : index_range(_test_var_names))
372 : {
373 982 : auto test_var_name = _test_var_names.at(i);
374 982 : _lfs.Register(test_var_name, std::make_shared<mfem::ParLinearForm>(_test_pfespaces.at(i)));
375 982 : _lfs.GetRef(test_var_name) = 0.0;
376 982 : }
377 :
378 1988 : for (auto & test_var_name : _test_var_names)
379 : {
380 : // Apply kernels
381 982 : auto lf = _lfs.GetShared(test_var_name);
382 982 : ApplyDomainLFIntegrators(test_var_name, lf, _kernels_map);
383 982 : ApplyBoundaryLFIntegrators(test_var_name, lf, _integrated_bc_map);
384 982 : lf->Assemble();
385 982 : }
386 :
387 : // Apply essential boundary conditions
388 1006 : ApplyEssentialBCs();
389 :
390 : // Eliminate trivially eliminated variables by subtracting contributions from linear forms
391 1006 : EliminateCoupledVariables();
392 1006 : }
393 :
394 : void
395 1006 : EquationSystem::BuildBilinearForms()
396 : {
397 : // Register bilinear forms
398 1988 : for (const auto i : index_range(_test_var_names))
399 : {
400 982 : auto test_var_name = _test_var_names.at(i);
401 982 : _blfs.Register(test_var_name, std::make_shared<mfem::ParBilinearForm>(_test_pfespaces.at(i)));
402 :
403 : // Apply kernels
404 982 : auto blf = _blfs.GetShared(test_var_name);
405 982 : blf->SetAssemblyLevel(_assembly_level);
406 982 : ApplyBoundaryBLFIntegrators<mfem::ParBilinearForm>(
407 982 : test_var_name, test_var_name, blf, _integrated_bc_map);
408 982 : ApplyDomainBLFIntegrators<mfem::ParBilinearForm>(
409 982 : test_var_name, test_var_name, blf, _kernels_map);
410 : // Assemble
411 982 : blf->Assemble();
412 982 : }
413 1006 : }
414 :
415 : void
416 1006 : EquationSystem::BuildMixedBilinearForms()
417 : {
418 : // Register mixed bilinear forms. Note that not all combinations may
419 : // have a kernel.
420 :
421 : // Create mblf for each test/coupled variable pair with an added kernel.
422 : // Mixed bilinear forms with coupled variables that are not trial variables are
423 : // associated with contributions from eliminated variables.
424 1988 : for (const auto i : index_range(_test_var_names))
425 : {
426 982 : auto test_var_name = _test_var_names.at(i);
427 982 : auto test_mblfs = std::make_shared<Moose::MFEM::NamedFieldsMap<mfem::ParMixedBilinearForm>>();
428 2866 : for (const auto j : index_range(_coupled_var_names))
429 : {
430 1884 : const auto & coupled_var_name = _coupled_var_names.at(j);
431 3768 : auto mblf = std::make_shared<mfem::ParMixedBilinearForm>(_coupled_pfespaces.at(j),
432 1884 : _test_pfespaces.at(i));
433 : // Register MixedBilinearForm if kernels exist for it, and assemble
434 : // kernels
435 1884 : if (_kernels_map.Has(test_var_name) &&
436 2894 : _kernels_map.Get(test_var_name)->Has(coupled_var_name) &&
437 1010 : test_var_name != coupled_var_name)
438 : {
439 68 : mblf->SetAssemblyLevel(_assembly_level);
440 : // Apply all mixed kernels with this test/trial pair
441 68 : ApplyDomainBLFIntegrators<mfem::ParMixedBilinearForm>(
442 68 : coupled_var_name, test_var_name, mblf, _kernels_map);
443 : // Assemble mixed bilinear forms
444 68 : mblf->Assemble();
445 : // Register mixed bilinear forms associated with a single trial variable
446 : // for the current test variable
447 68 : test_mblfs->Register(coupled_var_name, mblf);
448 : }
449 1884 : }
450 : // Register all mixed bilinear form sets associated with a single test
451 : // variable
452 982 : _mblfs.Register(test_var_name, test_mblfs);
453 982 : }
454 1006 : }
455 :
456 : void
457 1006 : EquationSystem::BuildEquationSystem()
458 : {
459 1006 : BuildBilinearForms();
460 1006 : BuildMixedBilinearForms();
461 1006 : BuildLinearForms();
462 1006 : }
463 :
464 117 : TimeDependentEquationSystem::TimeDependentEquationSystem(
465 117 : const Moose::MFEM::TimeDerivativeMap & time_derivative_map)
466 117 : : _dt_coef(1.0), _time_derivative_map(time_derivative_map)
467 : {
468 117 : }
469 :
470 : void
471 117 : TimeDependentEquationSystem::Init(Moose::MFEM::GridFunctions & gridfunctions,
472 : mfem::AssemblyLevel assembly_level)
473 : {
474 117 : EquationSystem::Init(gridfunctions, assembly_level);
475 238 : for (auto & test_var_name : _test_var_names)
476 121 : _td_var_ess_constraints.emplace_back(
477 242 : std::make_unique<mfem::ParGridFunction>(gridfunctions.Get(test_var_name)->ParFESpace()));
478 117 : }
479 :
480 : void
481 653 : TimeDependentEquationSystem::SetTimeStep(mfem::real_t dt)
482 : {
483 653 : if (fabs(dt - _dt_coef.constant) > 1.0e-12 * dt)
484 : {
485 208 : _dt_coef.constant = dt;
486 420 : for (auto test_var_name : _test_var_names)
487 : {
488 212 : auto blf = _blfs.Get(test_var_name);
489 212 : blf->Update();
490 212 : blf->Assemble();
491 212 : }
492 : }
493 653 : }
494 :
495 : void
496 117 : TimeDependentEquationSystem::SetTrialVariableNames()
497 : {
498 : // The TimeDependentEquationSystem operator expects to act on a vector of variable time
499 : // derivatives, so the trial variable must be the time derivative of the 'base' variable. The base
500 : // variable (test_var_name) without derivatives applied must also be coupled in implicit
501 : // timestepping schemes for the elimination of 'old' variable values from the previous timestep
502 238 : for (const auto & test_var_name : _test_var_names)
503 : {
504 121 : AddCoupledVariableNameIfMissing(test_var_name);
505 121 : AddCoupledVariableNameIfMissing(_time_derivative_map.getTimeDerivativeName(test_var_name));
506 : }
507 :
508 : // If a coupled variable does not have an equation associated with it,
509 : // add it to the set of eliminated variables.
510 238 : for (const auto & test_var_name : _test_var_names)
511 : {
512 : const auto time_derivative_test_var_name =
513 121 : _time_derivative_map.getTimeDerivativeName(test_var_name);
514 379 : for (const auto & coupled_var_name : _coupled_var_names)
515 : {
516 395 : if (time_derivative_test_var_name != coupled_var_name &&
517 137 : !VectorContainsName(_eliminated_var_names, coupled_var_name))
518 129 : _eliminated_var_names.push_back(coupled_var_name);
519 : }
520 121 : if (!VectorContainsName(_trial_var_names, time_derivative_test_var_name))
521 121 : _trial_var_names.push_back(time_derivative_test_var_name);
522 121 : }
523 117 : }
524 :
525 : void
526 235 : TimeDependentEquationSystem::AddKernel(std::shared_ptr<MFEMKernel> kernel)
527 : {
528 235 : const auto & trial_var_name = kernel->getTrialVariableName();
529 235 : const auto & test_var_name = kernel->getTestVariableName();
530 :
531 235 : if (_time_derivative_map.isTimeDerivative(trial_var_name))
532 : {
533 118 : AddTestVariableNameIfMissing(test_var_name);
534 118 : AddCoupledVariableNameIfMissing(trial_var_name);
535 :
536 : // Register new kernels map if not present for the test variable
537 118 : if (!_td_kernels_map.Has(test_var_name))
538 : {
539 : auto kernel_field_map =
540 114 : std::make_shared<Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>>();
541 114 : _td_kernels_map.Register(test_var_name, std::move(kernel_field_map));
542 114 : }
543 118 : if (!_td_kernels_map.Get(test_var_name)->Has(trial_var_name))
544 : {
545 118 : auto kernels = std::make_shared<std::vector<std::shared_ptr<MFEMKernel>>>();
546 118 : _td_kernels_map.Get(test_var_name)->Register(trial_var_name, std::move(kernels));
547 118 : }
548 118 : _td_kernels_map.GetRef(test_var_name).Get(trial_var_name)->push_back(std::move(kernel));
549 : }
550 : else
551 : {
552 117 : EquationSystem::AddKernel(kernel);
553 : }
554 235 : }
555 :
556 : void
557 770 : TimeDependentEquationSystem::BuildBilinearForms()
558 : {
559 770 : EquationSystem::BuildBilinearForms();
560 :
561 : // Build and assemble bilinear forms acting on time derivatives
562 1556 : for (const auto i : index_range(_test_var_names))
563 : {
564 786 : auto test_var_name = _test_var_names.at(i);
565 :
566 786 : _td_blfs.Register(test_var_name,
567 1572 : std::make_shared<mfem::ParBilinearForm>(_test_pfespaces.at(i)));
568 :
569 : // Apply kernels to td_blf
570 786 : auto td_blf = _td_blfs.GetShared(test_var_name);
571 786 : td_blf->SetAssemblyLevel(_assembly_level);
572 786 : ApplyBoundaryBLFIntegrators<mfem::ParBilinearForm>(
573 786 : test_var_name, test_var_name, td_blf, _integrated_bc_map);
574 786 : ApplyDomainBLFIntegrators<mfem::ParBilinearForm>(
575 786 : _time_derivative_map.getTimeDerivativeName(test_var_name),
576 : test_var_name,
577 : td_blf,
578 786 : _td_kernels_map);
579 :
580 : // Recover and scale integrators from blf. This is to apply the dt*du/dt contributions from the
581 : // operator on the trial variable in the implicit integration scheme
582 786 : auto blf = _blfs.GetShared(test_var_name);
583 786 : ScaleAndAddBLFIntegrators(blf, td_blf, _dt_coef.constant);
584 : // Assemble form
585 786 : td_blf->Assemble();
586 786 : }
587 770 : }
588 :
589 : void
590 770 : TimeDependentEquationSystem::BuildMixedBilinearForms()
591 : {
592 770 : EquationSystem::BuildMixedBilinearForms();
593 : // Register mixed bilinear forms. Note that not all combinations may
594 : // have a kernel.
595 :
596 : // Create mblf for each test/trial variable pair with an added kernel
597 1556 : for (const auto i : index_range(_test_var_names))
598 : {
599 786 : const auto & test_var_name = _test_var_names.at(i);
600 : auto test_td_mblfs =
601 786 : std::make_shared<Moose::MFEM::NamedFieldsMap<mfem::ParMixedBilinearForm>>();
602 1604 : for (const auto j : index_range(_trial_var_names))
603 : {
604 818 : const auto & trial_var_name = _trial_var_names.at(j);
605 1636 : auto td_mblf = std::make_shared<mfem::ParMixedBilinearForm>(_test_pfespaces.at(j),
606 818 : _test_pfespaces.at(i));
607 : // Register MixedBilinearForm if kernels exist for it, and assemble
608 : // kernels
609 818 : if (_td_kernels_map.Has(test_var_name) &&
610 1606 : _td_kernels_map.Get(test_var_name)->Has(trial_var_name) &&
611 788 : _time_derivative_map.getTimeDerivativeName(test_var_name) != trial_var_name)
612 : {
613 : // Apply all mixed kernels with this test/trial pair
614 16 : ApplyDomainBLFIntegrators<mfem::ParMixedBilinearForm>(
615 16 : trial_var_name, test_var_name, td_mblf, _td_kernels_map);
616 : }
617 : // Recover and scale integrators from the mblf acting on the time integral of the trial
618 : // variable corresponding to coupled_var_name. This is to apply the dt*du/dt contributions
619 : // from the operator on the trial variable in the implicit integration scheme
620 : const auto & trial_var_time_integral_name =
621 818 : _time_derivative_map.getTimeIntegralName(trial_var_name);
622 818 : if (_mblfs.Has(test_var_name) && _mblfs.Get(test_var_name)->Has(trial_var_time_integral_name))
623 : {
624 : // Recover and scale integrators from mblf. This is to apply the dt*du/dt contributions
625 : // from the operator on the trial variable in the implicit integration scheme
626 16 : auto mblf = _mblfs.Get(test_var_name)->GetShared(trial_var_time_integral_name);
627 16 : ScaleAndAddBLFIntegrators(mblf, td_mblf, _dt_coef.constant);
628 16 : }
629 : // If implicit contributions exist, scale them by timestep and add to td_mblf
630 818 : if (td_mblf->GetDBFI()->Size() || td_mblf->GetBBFI()->Size())
631 : {
632 : // Assemble mixed bilinear form
633 32 : td_mblf->SetAssemblyLevel(_assembly_level);
634 32 : td_mblf->Assemble();
635 : // Register mixed bilinear forms associated with a single trial variable
636 : // for the current test variable
637 32 : test_td_mblfs->Register(trial_var_name, td_mblf);
638 : }
639 818 : }
640 : // Register all mixed bilinear forms associated with a single test variable
641 786 : _td_mblfs.Register(test_var_name, test_td_mblfs);
642 786 : }
643 770 : }
644 :
645 : void
646 770 : TimeDependentEquationSystem::ApplyEssentialBCs()
647 : {
648 770 : _ess_tdof_lists.resize(_test_var_names.size());
649 1556 : for (const auto i : index_range(_test_var_names))
650 : {
651 786 : const auto & test_var_name = _test_var_names.at(i);
652 : const auto time_derivative_test_var_name =
653 786 : _time_derivative_map.getTimeDerivativeName(test_var_name);
654 786 : mfem::ParGridFunction & trial_gf = *(_var_ess_constraints.at(i));
655 786 : mfem::ParGridFunction & trial_gf_time_derivative = *(_td_var_ess_constraints.at(i));
656 786 : mfem::Array<int> global_ess_markers(trial_gf.ParFESpace()->GetParMesh()->bdr_attributes.Max());
657 786 : global_ess_markers = 0;
658 : // Set strongly constrained DoFs of trial_gf on essential boundaries and add markers for all
659 : // essential boundaries to the global_ess_markers array
660 786 : EquationSystem::ApplyEssentialBC(test_var_name, trial_gf, global_ess_markers);
661 : // Update solution values on Dirichlet values to be in terms of du/dt instead of u
662 786 : *_td_var_ess_constraints.at(i).get() = *(_var_ess_constraints.at(i).get());
663 786 : *_td_var_ess_constraints.at(i).get() -= *_eliminated_variables.Get(test_var_name);
664 786 : *_td_var_ess_constraints.at(i).get() /= _dt_coef.constant;
665 : // Apply any remaining Dirichlet BCs specified directly on du/dt
666 786 : EquationSystem::ApplyEssentialBC(
667 : time_derivative_test_var_name, trial_gf_time_derivative, global_ess_markers);
668 786 : trial_gf.ParFESpace()->GetEssentialTrueDofs(global_ess_markers, _ess_tdof_lists.at(i));
669 786 : }
670 770 : }
671 :
672 : void
673 770 : TimeDependentEquationSystem::EliminateCoupledVariables()
674 : {
675 : // Eliminate contributions from variables at previous timestep.
676 1556 : for (const auto i : index_range(_test_var_names))
677 : {
678 786 : auto & test_var_name = _test_var_names.at(i);
679 786 : auto blf = _blfs.Get(test_var_name);
680 786 : auto lf = _lfs.Get(test_var_name);
681 : // if implicit, add contribution to linear form from terms involving state
682 : // The AddMult method in mfem::BilinearForm is not defined for non-legacy assembly
683 786 : mfem::Vector lf_prev(lf->Size());
684 786 : blf->Mult(*_eliminated_variables.Get(test_var_name), lf_prev);
685 786 : *lf -= lf_prev;
686 786 : }
687 : // Eliminate contributions from other coupled variables.
688 770 : EquationSystem::EliminateCoupledVariables();
689 770 : }
690 :
691 : void
692 621 : TimeDependentEquationSystem::FormLegacySystem(mfem::OperatorHandle & op,
693 : mfem::BlockVector & truedXdt,
694 : mfem::BlockVector & trueRHS)
695 : {
696 : // Form linear system for operator acting on vector of du/dt
697 621 : AssembleJacobian(
698 621 : _td_blfs, _td_mblfs, _lfs, _ess_tdof_lists, _td_var_ess_constraints, op, truedXdt, trueRHS);
699 621 : }
700 :
701 : void
702 32 : TimeDependentEquationSystem::FormSystem(mfem::OperatorHandle & op,
703 : mfem::BlockVector & truedXdt,
704 : mfem::BlockVector & trueRHS)
705 : {
706 32 : auto & test_var_name = _test_var_names.at(0);
707 32 : auto td_blf = _td_blfs.Get(test_var_name);
708 32 : auto lf = _lfs.Get(test_var_name);
709 :
710 : // Form linear system for operator acting on vector of du/dt
711 32 : mfem::OperatorPtr aux_a;
712 32 : mfem::Vector aux_x, aux_rhs;
713 96 : td_blf->FormLinearSystem(
714 32 : _ess_tdof_lists.at(0), *(_td_var_ess_constraints.at(0)), *lf, aux_a, aux_x, aux_rhs);
715 :
716 32 : truedXdt.GetBlock(0) = aux_x;
717 32 : trueRHS.GetBlock(0) = aux_rhs;
718 32 : truedXdt.SyncFromBlocks();
719 32 : trueRHS.SyncFromBlocks();
720 :
721 : // Create monolithic matrix
722 32 : op.Reset(aux_a.Ptr());
723 32 : aux_a.SetOperatorOwner(false);
724 32 : }
725 :
726 : void
727 653 : TimeDependentEquationSystem::UpdateEquationSystem()
728 : {
729 653 : EquationSystem::BuildEquationSystem();
730 653 : }
731 :
732 : } // namespace Moose::MFEM
733 :
734 : #endif
|
