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