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