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 349 : EquationSystem::~EquationSystem() { DeleteAllBlocks(); }
19 :
20 : void
21 1137 : EquationSystem::DeleteAllBlocks()
22 : {
23 1933 : for (const auto i : make_range(_h_blocks.NumRows()))
24 1608 : for (const auto j : make_range(_h_blocks.NumCols()))
25 812 : delete _h_blocks(i, j);
26 1137 : _h_blocks.DeleteAll();
27 1137 : }
28 :
29 : bool
30 2293 : EquationSystem::VectorContainsName(const std::vector<std::string> & the_vector,
31 : const std::string & name) const
32 : {
33 :
34 2293 : auto iter = std::find(the_vector.begin(), the_vector.end(), name);
35 :
36 2293 : return (iter != the_vector.end());
37 : }
38 :
39 : void
40 811 : EquationSystem::AddCoupledVariableNameIfMissing(const std::string & coupled_var_name)
41 : {
42 811 : if (!VectorContainsName(_coupled_var_names, coupled_var_name))
43 458 : _coupled_var_names.push_back(coupled_var_name);
44 811 : }
45 :
46 : void
47 1024 : EquationSystem::AddTestVariableNameIfMissing(const std::string & test_var_name)
48 : {
49 1024 : if (!VectorContainsName(_test_var_names, test_var_name))
50 309 : _test_var_names.push_back(test_var_name);
51 1024 : }
52 :
53 : void
54 236 : EquationSystem::SetTrialVariableNames()
55 : {
56 : // If a coupled variable has an equation associated with it,
57 : // add it to the set of trial variables.
58 468 : for (const auto & coupled_var_name : _coupled_var_names)
59 : {
60 232 : if (VectorContainsName(_test_var_names, coupled_var_name))
61 196 : _trial_var_names.push_back(coupled_var_name);
62 : else
63 36 : _eliminated_var_names.push_back(coupled_var_name);
64 : }
65 236 : }
66 :
67 : void
68 442 : EquationSystem::AddKernel(std::shared_ptr<MFEMKernel> kernel)
69 : {
70 442 : AddTestVariableNameIfMissing(kernel->getTestVariableName());
71 442 : AddCoupledVariableNameIfMissing(kernel->getTrialVariableName());
72 442 : auto trial_var_name = kernel->getTrialVariableName();
73 442 : auto test_var_name = kernel->getTestVariableName();
74 442 : if (!_kernels_map.Has(test_var_name))
75 : {
76 : auto kernel_field_map =
77 309 : std::make_shared<Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>>();
78 309 : _kernels_map.Register(test_var_name, std::move(kernel_field_map));
79 309 : }
80 : // Register new kernels map if not present for the test/trial variable
81 : // pair
82 442 : if (!_kernels_map.Get(test_var_name)->Has(trial_var_name))
83 : {
84 353 : auto kernels = std::make_shared<std::vector<std::shared_ptr<MFEMKernel>>>();
85 353 : _kernels_map.Get(test_var_name)->Register(trial_var_name, std::move(kernels));
86 353 : }
87 442 : _kernels_map.GetRef(test_var_name).Get(trial_var_name)->push_back(std::move(kernel));
88 442 : }
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 446 : EquationSystem::AddEssentialBC(std::shared_ptr<MFEMEssentialBC> bc)
115 : {
116 446 : AddTestVariableNameIfMissing(bc->getTestVariableName());
117 446 : auto test_var_name = bc->getTestVariableName();
118 446 : if (!_essential_bc_map.Has(test_var_name))
119 : {
120 277 : auto bcs = std::make_shared<std::vector<std::shared_ptr<MFEMEssentialBC>>>();
121 277 : _essential_bc_map.Register(test_var_name, std::move(bcs));
122 277 : }
123 446 : _essential_bc_map.GetRef(test_var_name).push_back(std::move(bc));
124 446 : }
125 :
126 : void
127 349 : EquationSystem::Init(Moose::MFEM::GridFunctions & gridfunctions,
128 : const Moose::MFEM::FESpaces & /*fespaces*/,
129 : mfem::AssemblyLevel assembly_level)
130 : {
131 349 : _assembly_level = assembly_level;
132 :
133 658 : for (auto & test_var_name : _test_var_names)
134 : {
135 309 : 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 309 : _test_pfespaces.push_back(gridfunctions.Get(test_var_name)->ParFESpace());
144 : // Create auxiliary gridfunctions for storing essential constraints from Dirichlet conditions
145 309 : _var_ess_constraints.emplace_back(
146 618 : std::make_unique<mfem::ParGridFunction>(gridfunctions.Get(test_var_name)->ParFESpace()));
147 : }
148 :
149 : // Store pointers to FESpaces of all coupled variables
150 807 : for (auto & coupled_var_name : _coupled_var_names)
151 458 : _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 349 : SetTrialVariableNames();
155 :
156 : // Store pointers to coupled variable GridFunctions that are to be eliminated prior to forming the
157 : // jacobian
158 498 : for (auto & eliminated_var_name : _eliminated_var_names)
159 149 : _eliminated_variables.Register(eliminated_var_name,
160 298 : gridfunctions.GetShared(eliminated_var_name));
161 349 : }
162 :
163 : void
164 990 : EquationSystem::ApplyEssentialBCs()
165 : {
166 990 : _ess_tdof_lists.resize(_test_var_names.size());
167 1940 : for (const auto i : index_range(_test_var_names))
168 : {
169 950 : auto test_var_name = _test_var_names.at(i);
170 950 : if (!_essential_bc_map.Has(test_var_name))
171 32 : continue;
172 :
173 : // Set default value of gridfunction used in essential BC. Values
174 : // overwritten in applyEssentialBCs
175 918 : mfem::ParGridFunction & trial_gf(*(_var_ess_constraints.at(i)));
176 918 : auto * const pmesh = _test_pfespaces.at(i)->GetParMesh();
177 : mooseAssert(pmesh, "parallel mesh is null");
178 :
179 918 : auto bcs = _essential_bc_map.GetRef(test_var_name);
180 918 : mfem::Array<int> global_ess_markers(pmesh->bdr_attributes.Max());
181 918 : global_ess_markers = 0;
182 2576 : for (auto & bc : bcs)
183 : {
184 1658 : bc->ApplyBC(trial_gf);
185 :
186 1658 : mfem::Array<int> ess_bdrs(bc->getBoundaryMarkers());
187 7812 : for (auto it = 0; it != pmesh->bdr_attributes.Max(); ++it)
188 : {
189 6154 : global_ess_markers[it] = std::max(global_ess_markers[it], ess_bdrs[it]);
190 : }
191 1658 : }
192 918 : trial_gf.FESpace()->GetEssentialTrueDofs(global_ess_markers, _ess_tdof_lists.at(i));
193 950 : }
194 990 : }
195 :
196 : void
197 990 : EquationSystem::EliminateCoupledVariables()
198 : {
199 1940 : for (const auto & test_var_name : _test_var_names)
200 : {
201 950 : auto lf = _lfs.Get(test_var_name);
202 1740 : for (const auto & eliminated_var_name : _eliminated_var_names)
203 : {
204 790 : if (_mblfs.Has(test_var_name) && _mblfs.Get(test_var_name)->Has(eliminated_var_name))
205 : {
206 36 : auto mblf = _mblfs.Get(test_var_name)->Get(eliminated_var_name);
207 36 : mblf->AddMult(*_eliminated_variables.Get(eliminated_var_name), *lf, -1.0);
208 : }
209 : }
210 : }
211 990 : }
212 :
213 : void
214 829 : EquationSystem::FormLinearSystem(mfem::OperatorHandle & op,
215 : mfem::BlockVector & trueX,
216 : mfem::BlockVector & trueRHS)
217 : {
218 :
219 829 : switch (_assembly_level)
220 : {
221 788 : case mfem::AssemblyLevel::LEGACY:
222 788 : FormLegacySystem(op, trueX, trueRHS);
223 788 : 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 829 : }
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 179 : EquationSystem::FormLegacySystem(mfem::OperatorHandle & op,
259 : mfem::BlockVector & trueX,
260 : mfem::BlockVector & trueRHS)
261 : {
262 :
263 : // Allocate block operator
264 179 : DeleteAllBlocks();
265 179 : _h_blocks.SetSize(_test_var_names.size(), _test_var_names.size());
266 : // Form diagonal blocks.
267 366 : for (const auto i : index_range(_test_var_names))
268 : {
269 187 : auto & test_var_name = _test_var_names.at(i);
270 187 : auto blf = _blfs.Get(test_var_name);
271 187 : auto lf = _lfs.Get(test_var_name);
272 187 : mfem::Vector aux_x, aux_rhs;
273 187 : mfem::HypreParMatrix * aux_a = new mfem::HypreParMatrix;
274 561 : blf->FormLinearSystem(
275 187 : _ess_tdof_lists.at(i), *(_var_ess_constraints.at(i)), *lf, *aux_a, aux_x, aux_rhs);
276 187 : _h_blocks(i, i) = aux_a;
277 187 : trueX.GetBlock(i) = aux_x;
278 187 : trueRHS.GetBlock(i) = aux_rhs;
279 187 : }
280 :
281 : // Form off-diagonal blocks
282 366 : for (const auto i : index_range(_test_var_names))
283 : {
284 187 : auto test_var_name = _test_var_names.at(i);
285 390 : for (const auto j : index_range(_trial_var_names))
286 : {
287 203 : auto trial_var_name = _trial_var_names.at(j);
288 :
289 203 : mfem::Vector aux_x, aux_rhs;
290 203 : mfem::ParLinearForm aux_lf(_test_pfespaces.at(i));
291 203 : aux_lf = 0.0;
292 203 : 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 203 : }
307 187 : }
308 : // Sync memory
309 366 : for (const auto i : index_range(_test_var_names))
310 : {
311 187 : trueX.GetBlock(i).SyncAliasMemory(trueX);
312 187 : trueRHS.GetBlock(i).SyncAliasMemory(trueRHS);
313 : }
314 :
315 : // Create monolithic matrix
316 179 : op.Reset(mfem::HypreParMatrixFromBlocks(_h_blocks));
317 179 : }
318 :
319 : void
320 829 : EquationSystem::BuildJacobian(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS)
321 : {
322 829 : height = trueX.Size();
323 829 : width = trueRHS.Size();
324 829 : FormLinearSystem(_jacobian, trueX, trueRHS);
325 829 : }
326 :
327 : void
328 1658 : EquationSystem::Mult(const mfem::Vector & x, mfem::Vector & residual) const
329 : {
330 1658 : _jacobian->Mult(x, residual);
331 1658 : x.HostRead();
332 1658 : residual.HostRead();
333 1658 : }
334 :
335 : mfem::Operator &
336 829 : EquationSystem::GetGradient(const mfem::Vector &) const
337 : {
338 829 : return *_jacobian;
339 : }
340 :
341 : void
342 188 : EquationSystem::RecoverFEMSolution(mfem::BlockVector & trueX,
343 : Moose::MFEM::GridFunctions & gridfunctions)
344 : {
345 384 : for (const auto i : index_range(_trial_var_names))
346 : {
347 196 : auto & trial_var_name = _trial_var_names.at(i);
348 196 : trueX.GetBlock(i).SyncAliasMemory(trueX);
349 196 : gridfunctions.Get(trial_var_name)->Distribute(&(trueX.GetBlock(i)));
350 : }
351 188 : }
352 :
353 : void
354 990 : EquationSystem::BuildLinearForms()
355 : {
356 : // Register linear forms
357 1940 : for (const auto i : index_range(_test_var_names))
358 : {
359 950 : auto test_var_name = _test_var_names.at(i);
360 950 : _lfs.Register(test_var_name, std::make_shared<mfem::ParLinearForm>(_test_pfespaces.at(i)));
361 950 : _lfs.GetRef(test_var_name) = 0.0;
362 950 : }
363 :
364 : // Apply boundary conditions
365 990 : ApplyEssentialBCs();
366 :
367 1940 : for (auto & test_var_name : _test_var_names)
368 : {
369 : // Apply kernels
370 950 : auto lf = _lfs.GetShared(test_var_name);
371 950 : ApplyDomainLFIntegrators(test_var_name, lf, _kernels_map);
372 950 : ApplyBoundaryLFIntegrators(test_var_name, lf, _integrated_bc_map);
373 950 : lf->Assemble();
374 950 : }
375 :
376 : // Eliminate trivially eliminated variables by subtracting contributions from linear forms
377 990 : EliminateCoupledVariables();
378 990 : }
379 :
380 : void
381 990 : EquationSystem::BuildBilinearForms()
382 : {
383 : // Register bilinear forms
384 1940 : for (const auto i : index_range(_test_var_names))
385 : {
386 950 : auto test_var_name = _test_var_names.at(i);
387 950 : _blfs.Register(test_var_name, std::make_shared<mfem::ParBilinearForm>(_test_pfespaces.at(i)));
388 :
389 : // Apply kernels
390 950 : auto blf = _blfs.GetShared(test_var_name);
391 950 : blf->SetAssemblyLevel(_assembly_level);
392 950 : ApplyBoundaryBLFIntegrators<mfem::ParBilinearForm>(
393 950 : test_var_name, test_var_name, blf, _integrated_bc_map);
394 950 : ApplyDomainBLFIntegrators<mfem::ParBilinearForm>(
395 950 : test_var_name, test_var_name, blf, _kernels_map);
396 : // Assemble
397 950 : blf->Assemble();
398 950 : }
399 990 : }
400 :
401 : void
402 990 : 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 1940 : for (const auto i : index_range(_test_var_names))
409 : {
410 950 : auto test_var_name = _test_var_names.at(i);
411 950 : auto test_mblfs = std::make_shared<Moose::MFEM::NamedFieldsMap<mfem::ParMixedBilinearForm>>();
412 2706 : for (const auto j : index_range(_coupled_var_names))
413 : {
414 1756 : const auto & coupled_var_name = _coupled_var_names.at(j);
415 3512 : auto mblf = std::make_shared<mfem::ParMixedBilinearForm>(_coupled_pfespaces.at(j),
416 1756 : _test_pfespaces.at(i));
417 : // Register MixedBilinearForm if kernels exist for it, and assemble
418 : // kernels
419 1756 : if (_kernels_map.Has(test_var_name) &&
420 2750 : _kernels_map.Get(test_var_name)->Has(coupled_var_name) &&
421 994 : test_var_name != coupled_var_name)
422 : {
423 52 : mblf->SetAssemblyLevel(_assembly_level);
424 : // Apply all mixed kernels with this test/trial pair
425 52 : ApplyDomainBLFIntegrators<mfem::ParMixedBilinearForm>(
426 52 : coupled_var_name, test_var_name, mblf, _kernels_map);
427 : // Assemble mixed bilinear forms
428 52 : mblf->Assemble();
429 : // Register mixed bilinear forms associated with a single trial variable
430 : // for the current test variable
431 52 : test_mblfs->Register(coupled_var_name, mblf);
432 : }
433 1756 : }
434 : // Register all mixed bilinear form sets associated with a single test
435 : // variable
436 950 : _mblfs.Register(test_var_name, test_mblfs);
437 950 : }
438 990 : }
439 :
440 : void
441 990 : EquationSystem::BuildEquationSystem()
442 : {
443 990 : BuildBilinearForms();
444 990 : BuildMixedBilinearForms();
445 990 : BuildLinearForms();
446 990 : }
447 :
448 113 : TimeDependentEquationSystem::TimeDependentEquationSystem() : _dt_coef(1.0) {}
449 :
450 : void
451 233 : 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 233 : 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 233 : EquationSystem::AddCoupledVariableNameIfMissing(coupled_var_name);
459 233 : }
460 :
461 : void
462 641 : TimeDependentEquationSystem::SetTimeStep(double dt)
463 : {
464 641 : if (fabs(dt - _dt_coef.constant) > 1.0e-12 * dt)
465 : {
466 204 : _dt_coef.constant = dt;
467 408 : for (auto test_var_name : _test_var_names)
468 : {
469 204 : auto blf = _blfs.Get(test_var_name);
470 204 : blf->Update();
471 204 : blf->Assemble();
472 204 : }
473 : }
474 641 : }
475 :
476 : void
477 113 : TimeDependentEquationSystem::SetTrialVariableNames()
478 : {
479 : // If a coupled variable has an equation associated with it,
480 : // add it to the set of trial variables.
481 339 : for (const auto & coupled_var_name : _coupled_var_names)
482 : {
483 452 : for (const auto & test_var_name : _test_var_names)
484 : {
485 226 : const auto time_derivative_test_var_name = GetTimeDerivativeName(test_var_name);
486 226 : if (time_derivative_test_var_name == coupled_var_name)
487 : {
488 113 : if (!VectorContainsName(_trial_var_names, coupled_var_name))
489 113 : _trial_var_names.push_back(coupled_var_name);
490 : }
491 : else
492 : {
493 113 : if (!VectorContainsName(_eliminated_var_names, coupled_var_name))
494 113 : _eliminated_var_names.push_back(coupled_var_name);
495 : }
496 226 : }
497 : }
498 113 : }
499 :
500 : void
501 219 : TimeDependentEquationSystem::AddKernel(std::shared_ptr<MFEMKernel> kernel)
502 : {
503 219 : if (kernel->getTrialVariableName() == GetTimeDerivativeName(kernel->getTestVariableName()))
504 : {
505 106 : auto trial_var_name = kernel->getTrialVariableName();
506 106 : auto test_var_name = kernel->getTestVariableName();
507 106 : AddTestVariableNameIfMissing(test_var_name);
508 106 : AddCoupledVariableNameIfMissing(test_var_name);
509 :
510 106 : if (!_td_kernels_map.Has(test_var_name))
511 : {
512 : auto kernel_field_map =
513 106 : std::make_shared<Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>>();
514 106 : _td_kernels_map.Register(test_var_name, std::move(kernel_field_map));
515 106 : }
516 : // Register new kernels map if not present for the test variable
517 106 : if (!_td_kernels_map.Get(test_var_name)->Has(test_var_name))
518 : {
519 106 : auto kernels = std::make_shared<std::vector<std::shared_ptr<MFEMKernel>>>();
520 106 : _td_kernels_map.Get(test_var_name)->Register(test_var_name, std::move(kernels));
521 106 : }
522 106 : _td_kernels_map.GetRef(test_var_name).Get(test_var_name)->push_back(std::move(kernel));
523 106 : }
524 : else
525 : {
526 113 : EquationSystem::AddKernel(kernel);
527 : }
528 219 : }
529 :
530 : void
531 754 : TimeDependentEquationSystem::BuildBilinearForms()
532 : {
533 754 : EquationSystem::BuildBilinearForms();
534 :
535 : // Build and assemble bilinear forms acting on time derivatives
536 1508 : for (const auto i : index_range(_test_var_names))
537 : {
538 754 : auto test_var_name = _test_var_names.at(i);
539 :
540 754 : _td_blfs.Register(test_var_name,
541 1508 : std::make_shared<mfem::ParBilinearForm>(_test_pfespaces.at(i)));
542 :
543 : // Apply kernels to td_blf
544 754 : auto td_blf = _td_blfs.GetShared(test_var_name);
545 754 : td_blf->SetAssemblyLevel(_assembly_level);
546 754 : ApplyBoundaryBLFIntegrators<mfem::ParBilinearForm>(
547 754 : test_var_name, test_var_name, td_blf, _integrated_bc_map);
548 754 : ApplyDomainBLFIntegrators<mfem::ParBilinearForm>(
549 754 : 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 754 : auto blf = _blfs.Get(test_var_name);
554 754 : auto integs = blf->GetDBFI();
555 754 : auto b_integs = blf->GetBBFI();
556 754 : auto markers = blf->GetBBFI_Marker();
557 :
558 : // If implicit contributions exist, scale them by timestep and add to td_blf
559 754 : if (integs->Size() || b_integs->Size())
560 : {
561 754 : mfem::SumIntegrator * sum = new mfem::SumIntegrator(false);
562 754 : ScaleIntegrator * scaled_sum = new ScaleIntegrator(sum, _dt_coef.constant, true);
563 :
564 1508 : for (int i = 0; i < integs->Size(); ++i)
565 754 : sum->AddIntegrator(*integs[i]);
566 :
567 810 : 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 754 : td_blf->AddDomainIntegrator(scaled_sum);
573 : }
574 :
575 : // Assemble form
576 754 : td_blf->Assemble();
577 754 : }
578 754 : }
579 :
580 : void
581 609 : TimeDependentEquationSystem::FormLegacySystem(mfem::OperatorHandle & op,
582 : mfem::BlockVector & truedXdt,
583 : mfem::BlockVector & trueRHS)
584 : {
585 :
586 : // Allocate block operator
587 609 : DeleteAllBlocks();
588 609 : _h_blocks.SetSize(_test_var_names.size(), _test_var_names.size());
589 : // Form diagonal blocks.
590 1218 : for (const auto i : index_range(_test_var_names))
591 : {
592 609 : auto & test_var_name = _test_var_names.at(i);
593 609 : auto td_blf = _td_blfs.Get(test_var_name);
594 609 : auto blf = _blfs.Get(test_var_name);
595 609 : 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 609 : blf->AddMult(*_eliminated_variables.Get(test_var_name), *lf, -1.0);
599 : // }
600 609 : mfem::Vector aux_x, aux_rhs;
601 : // Update solution values on Dirichlet values to be in terms of du/dt instead of u
602 609 : mfem::Vector bc_x = *(_var_ess_constraints.at(i).get());
603 609 : bc_x -= *_eliminated_variables.Get(test_var_name);
604 609 : bc_x /= _dt_coef.constant;
605 :
606 : // Form linear system for operator acting on vector of du/dt
607 609 : mfem::HypreParMatrix * aux_a = new mfem::HypreParMatrix;
608 609 : td_blf->FormLinearSystem(_ess_tdof_lists.at(i), bc_x, *lf, *aux_a, aux_x, aux_rhs);
609 609 : _h_blocks(i, i) = aux_a;
610 609 : truedXdt.GetBlock(i) = aux_x;
611 609 : trueRHS.GetBlock(i) = aux_rhs;
612 609 : }
613 :
614 609 : truedXdt.SyncFromBlocks();
615 609 : trueRHS.SyncFromBlocks();
616 :
617 : // Create monolithic matrix
618 609 : op.Reset(mfem::HypreParMatrixFromBlocks(_h_blocks));
619 609 : }
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 641 : TimeDependentEquationSystem::UpdateEquationSystem()
660 : {
661 641 : EquationSystem::BuildEquationSystem();
662 641 : }
663 :
664 : } // namespace Moose::MFEM
665 :
666 : #endif
|
