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