Line data Source code
1 : #ifdef MOOSE_MFEM_ENABLED
2 :
3 : #include "ComplexEquationSystem.h"
4 : #include "CoefficientManager.h"
5 : #include "libmesh/int_range.h"
6 :
7 : namespace Moose::MFEM
8 : {
9 :
10 : void
11 63 : ComplexEquationSystem::Init(GridFunctions & gridfunctions,
12 : ComplexGridFunctions & cmplx_gridfunctions,
13 : mfem::AssemblyLevel assembly_level)
14 : {
15 63 : _assembly_level = assembly_level;
16 :
17 63 : if (gridfunctions.size())
18 0 : mooseError("Mixing real and complex variables is currently not supported.");
19 :
20 105 : for (auto & test_var_name : _test_var_names)
21 : {
22 42 : if (!cmplx_gridfunctions.Has(test_var_name))
23 : {
24 0 : mooseError("MFEM complex variable ",
25 : test_var_name,
26 : " requested by equation system during initialization was "
27 : "not found in gridfunctions");
28 : }
29 : // Store pointers to test FESpaces
30 42 : _test_pfespaces.push_back(cmplx_gridfunctions.Get(test_var_name)->ParFESpace());
31 : // Create auxiliary gridfunctions for storing essential constraints from Dirichlet conditions
32 42 : _cmplx_var_ess_constraints.emplace_back(std::make_unique<mfem::ParComplexGridFunction>(
33 84 : cmplx_gridfunctions.Get(test_var_name)->ParFESpace()));
34 : }
35 :
36 : // Store pointers to FESpaces of all coupled variables
37 105 : for (auto & coupled_var_name : _coupled_var_names)
38 42 : _coupled_pfespaces.push_back(cmplx_gridfunctions.Get(coupled_var_name)->ParFESpace());
39 :
40 : // Extract which coupled variables are to be trivially eliminated and which are trial variables
41 63 : SetTrialVariableNames();
42 :
43 : // Store pointers to coupled variable ComplexGridFunctions that are to be eliminated prior to
44 : // forming the jacobian
45 63 : for (auto & eliminated_var_name : _eliminated_var_names)
46 0 : _cmplx_eliminated_variables.Register(eliminated_var_name,
47 0 : cmplx_gridfunctions.GetShared(eliminated_var_name));
48 :
49 : // Get a reference to the complex GridFunctions
50 63 : _complex_gfuncs = &cmplx_gridfunctions;
51 63 : }
52 :
53 : void
54 42 : ComplexEquationSystem::BuildEquationSystem()
55 : {
56 42 : BuildBilinearForms();
57 42 : BuildLinearForms();
58 42 : }
59 :
60 : void
61 42 : ComplexEquationSystem::BuildLinearForms()
62 : {
63 : // Register linear forms
64 84 : for (const auto i : index_range(_test_var_names))
65 : {
66 42 : auto test_var_name = _test_var_names.at(i);
67 42 : _clfs.Register(test_var_name,
68 84 : std::make_shared<mfem::ParComplexLinearForm>(_test_pfespaces.at(i)));
69 42 : _clfs.GetRef(test_var_name) = 0.0;
70 42 : }
71 : // Apply boundary conditions
72 42 : ApplyEssentialBCs();
73 :
74 84 : for (auto & test_var_name : _test_var_names)
75 : {
76 : // Apply kernels
77 42 : auto clf = _clfs.GetShared(test_var_name);
78 42 : ApplyDomainLFIntegrators(test_var_name, clf, _cmplx_kernels_map);
79 42 : ApplyBoundaryLFIntegrators(test_var_name, clf, _cmplx_integrated_bc_map);
80 42 : clf->Assemble();
81 42 : }
82 42 : }
83 :
84 : void
85 42 : ComplexEquationSystem::BuildBilinearForms()
86 : {
87 : // Register bilinear forms
88 84 : for (const auto i : index_range(_test_var_names))
89 : {
90 42 : auto test_var_name = _test_var_names.at(i);
91 42 : _slfs.Register(test_var_name,
92 84 : std::make_shared<mfem::ParSesquilinearForm>(_test_pfespaces.at(i)));
93 :
94 : // Apply kernels
95 42 : auto slf = _slfs.GetShared(test_var_name);
96 42 : slf->SetAssemblyLevel(_assembly_level);
97 42 : ApplyBoundaryBLFIntegrators<mfem::ParSesquilinearForm>(
98 42 : test_var_name, test_var_name, slf, _cmplx_integrated_bc_map);
99 42 : ApplyDomainBLFIntegrators<mfem::ParSesquilinearForm>(
100 42 : test_var_name, test_var_name, slf, _cmplx_kernels_map);
101 : // Assemble
102 42 : slf->Assemble();
103 42 : }
104 42 : }
105 :
106 : void
107 42 : ComplexEquationSystem::ApplyComplexEssentialBC(const std::string & var_name,
108 : mfem::ParComplexGridFunction & trial_gf,
109 : mfem::Array<int> & global_ess_markers)
110 : {
111 42 : if (_cmplx_essential_bc_map.Has(var_name))
112 105 : for (auto & bc : _cmplx_essential_bc_map.GetRef(var_name))
113 : {
114 : // Set constrained DoFs values on essential boundaries
115 63 : bc->ApplyBC(trial_gf);
116 : // Fetch marker array labelling essential boundaries of current BC
117 63 : mfem::Array<int> ess_bdrs(bc->getBoundaryMarkers());
118 : // Add these boundary markers to the set of markers labelling all essential boundaries
119 392 : for (const auto i : make_range(ess_bdrs.Size()))
120 329 : global_ess_markers[i] |= ess_bdrs[i];
121 63 : }
122 42 : }
123 :
124 : void
125 42 : ComplexEquationSystem::ApplyEssentialBCs()
126 : {
127 42 : _ess_tdof_lists.resize(_trial_var_names.size());
128 42 : _ess_markers.resize(_trial_var_names.size());
129 84 : for (const auto i : index_range(_trial_var_names))
130 : {
131 42 : const auto & trial_var_name = _trial_var_names.at(i);
132 42 : mfem::ParComplexGridFunction & trial_gf = *_cmplx_var_ess_constraints.at(i);
133 :
134 : // Make sure we update the size, if this mesh has changed recently for instance
135 42 : trial_gf.Update();
136 :
137 : // Initial guess for iterative solvers (initial condition or the previous time step solution)
138 42 : static_cast<mfem::Vector &>(trial_gf) = _complex_gfuncs->GetRef(trial_var_name);
139 :
140 42 : _ess_markers.at(i).SetSize(trial_gf.ParFESpace()->GetParMesh()->bdr_attributes.Max(), 0);
141 : // Set strongly constrained DoFs of trial_gf on essential boundaries and add markers for all
142 : // essential boundaries to the _ess_markers array
143 42 : ApplyComplexEssentialBC(trial_var_name, trial_gf, _ess_markers.at(i));
144 42 : trial_gf.ParFESpace()->GetEssentialTrueDofs(_ess_markers.at(i), _ess_tdof_lists.at(i));
145 : }
146 42 : }
147 :
148 : void
149 63 : ComplexEquationSystem::AddComplexKernel(std::shared_ptr<MFEMComplexKernel> kernel)
150 : {
151 63 : const auto & trial_var_name = kernel->getTrialVariableName();
152 63 : const auto & test_var_name = kernel->getTestVariableName();
153 63 : AddCoupledVariableNameIfMissing(trial_var_name);
154 63 : AddTestVariableNameIfMissing(test_var_name);
155 : // Register new complex kernels map if not present for the test variable
156 63 : if (!_cmplx_kernels_map.Has(test_var_name))
157 : {
158 : auto kernel_field_map =
159 42 : std::make_shared<NamedFieldsMap<std::vector<std::shared_ptr<MFEMComplexKernel>>>>();
160 42 : _cmplx_kernels_map.Register(test_var_name, std::move(kernel_field_map));
161 42 : }
162 : // Register new complex kernels map if not present for the test/trial variable pair
163 63 : if (!_cmplx_kernels_map.Get(test_var_name)->Has(trial_var_name))
164 : {
165 42 : auto kernels = std::make_shared<std::vector<std::shared_ptr<MFEMComplexKernel>>>();
166 42 : _cmplx_kernels_map.Get(test_var_name)->Register(trial_var_name, std::move(kernels));
167 42 : }
168 63 : _cmplx_kernels_map.GetRef(test_var_name).Get(trial_var_name)->push_back(std::move(kernel));
169 63 : }
170 :
171 : void
172 14 : ComplexEquationSystem::AddComplexIntegratedBC(std::shared_ptr<MFEMComplexIntegratedBC> bc)
173 : {
174 14 : const auto & trial_var_name = bc->getTrialVariableName();
175 14 : const auto & test_var_name = bc->getTestVariableName();
176 14 : AddCoupledVariableNameIfMissing(trial_var_name);
177 14 : AddTestVariableNameIfMissing(test_var_name);
178 : // Register new complex integrated bc map if not present for the test variable
179 14 : if (!_cmplx_integrated_bc_map.Has(test_var_name))
180 : {
181 : auto integrated_bc_field_map =
182 7 : std::make_shared<NamedFieldsMap<std::vector<std::shared_ptr<MFEMComplexIntegratedBC>>>>();
183 7 : _cmplx_integrated_bc_map.Register(test_var_name, std::move(integrated_bc_field_map));
184 7 : }
185 : // Register new complex integrated bc map if not present for the test/trial variable pair
186 14 : if (!_cmplx_integrated_bc_map.Get(test_var_name)->Has(trial_var_name))
187 : {
188 7 : auto bcs = std::make_shared<std::vector<std::shared_ptr<MFEMComplexIntegratedBC>>>();
189 7 : _cmplx_integrated_bc_map.Get(test_var_name)->Register(trial_var_name, std::move(bcs));
190 7 : }
191 14 : _cmplx_integrated_bc_map.GetRef(test_var_name).Get(trial_var_name)->push_back(std::move(bc));
192 14 : }
193 :
194 : void
195 63 : ComplexEquationSystem::AddComplexEssentialBCs(std::shared_ptr<MFEMComplexEssentialBC> bc)
196 : {
197 63 : const auto & test_var_name = bc->getTestVariableName();
198 63 : AddTestVariableNameIfMissing(test_var_name);
199 : // Register new complex essential bc map if not present for the test variable
200 63 : if (!_cmplx_essential_bc_map.Has(test_var_name))
201 : {
202 42 : auto bcs = std::make_shared<std::vector<std::shared_ptr<MFEMComplexEssentialBC>>>();
203 42 : _cmplx_essential_bc_map.Register(test_var_name, std::move(bcs));
204 42 : }
205 63 : _cmplx_essential_bc_map.GetRef(test_var_name).push_back(std::move(bc));
206 63 : }
207 :
208 : void
209 0 : ComplexEquationSystem::FormSystemOperator(mfem::OperatorHandle & op,
210 : mfem::BlockVector & trueX,
211 : mfem::BlockVector & trueRHS)
212 : {
213 0 : auto & test_var_name = _test_var_names.at(0);
214 0 : mfem::Vector aux_x, aux_rhs;
215 0 : mfem::OperatorPtr aux_a;
216 :
217 0 : auto slf = _slfs.Get(test_var_name);
218 0 : slf->FormLinearSystem(_ess_tdof_lists.at(0),
219 0 : *_cmplx_var_ess_constraints.at(0),
220 0 : *_clfs.Get(test_var_name),
221 : aux_a,
222 : aux_x,
223 : aux_rhs,
224 : /*copy_interior=*/true);
225 :
226 0 : trueX.GetBlock(0) = aux_x;
227 0 : trueRHS.GetBlock(0) = aux_rhs;
228 0 : trueX.SyncFromBlocks();
229 0 : trueRHS.SyncFromBlocks();
230 :
231 0 : op.Reset(aux_a.Ptr());
232 0 : aux_a.SetOperatorOwner(false);
233 0 : }
234 :
235 : void
236 42 : ComplexEquationSystem::FormSystemMatrix(mfem::OperatorHandle & op,
237 : mfem::BlockVector & trueX,
238 : mfem::BlockVector & trueRHS)
239 : {
240 :
241 : // Allocate block operator
242 42 : DeleteHBlocks();
243 42 : _h_blocks.SetSize(_test_var_names.size(), _trial_var_names.size());
244 42 : _h_blocks = nullptr;
245 : // Zero out RHS and sync memory
246 42 : trueRHS = 0.0;
247 42 : trueRHS.SyncToBlocks();
248 :
249 : // Form diagonal blocks.
250 84 : for (const auto i : index_range(_test_var_names))
251 : {
252 42 : auto & test_var_name = _test_var_names.at(i);
253 :
254 42 : mfem::Vector aux_x, aux_rhs;
255 42 : mfem::OperatorHandle aux_a;
256 :
257 42 : auto slf = _slfs.Get(test_var_name);
258 42 : slf->FormLinearSystem(_ess_tdof_lists.at(i),
259 42 : *_cmplx_var_ess_constraints.at(i),
260 42 : *_clfs.Get(test_var_name),
261 : aux_a,
262 : aux_x,
263 : aux_rhs,
264 : /*copy_interior=*/true);
265 42 : trueX.GetBlock(i) = aux_x;
266 42 : trueRHS.GetBlock(i) = aux_rhs;
267 42 : _h_blocks(i, i) = aux_a.As<mfem::ComplexHypreParMatrix>()->GetSystemMatrix();
268 42 : }
269 : // Sync memory
270 42 : trueX.SyncFromBlocks();
271 42 : trueRHS.SyncFromBlocks();
272 :
273 : // Create monolithic matrix
274 42 : op.Reset(mfem::HypreParMatrixFromBlocks(_h_blocks));
275 42 : }
276 :
277 : // Equation system Mult
278 : void
279 0 : ComplexEquationSystem::Mult(const mfem::Vector & x, mfem::Vector & residual) const
280 : {
281 0 : _linear_operator->Mult(x, residual);
282 0 : x.HostRead();
283 0 : residual.HostRead();
284 0 : }
285 :
286 : void
287 42 : ComplexEquationSystem::SetTrialVariablesFromTrueVectors(const mfem::BlockVector & trueX) const
288 : {
289 84 : for (const auto i : index_range(_trial_var_names))
290 : {
291 42 : auto & trial_var_name = _trial_var_names.at(i);
292 42 : trueX.GetBlock(i).SyncMemory(trueX);
293 42 : _complex_gfuncs->Get(trial_var_name)->Distribute(&(trueX.GetBlock(i)));
294 : }
295 : // Solution variables changed: stored projections of solution-dependent coefficients are stale.
296 42 : if (_coefficient_manager)
297 42 : _coefficient_manager->markSolutionChanged();
298 42 : }
299 : }
300 :
301 : #endif
|
