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 : {
113 42 : auto & bcs = _cmplx_essential_bc_map.GetRef(var_name);
114 105 : for (auto & bc : bcs)
115 : {
116 : // Set constrained DoFs values on essential boundaries
117 63 : bc->ApplyBC(trial_gf);
118 : // Fetch marker array labelling essential boundaries of current BC
119 63 : mfem::Array<int> ess_bdrs(bc->getBoundaryMarkers());
120 : // Add these boundary markers to the set of markers labelling all essential boundaries
121 392 : for (const auto i : make_range(trial_gf.ParFESpace()->GetParMesh()->bdr_attributes.Max()))
122 329 : global_ess_markers[i] = std::max(global_ess_markers[i], ess_bdrs[i]);
123 63 : }
124 : }
125 42 : }
126 :
127 : void
128 42 : ComplexEquationSystem::ApplyEssentialBCs()
129 : {
130 42 : _ess_tdof_lists.resize(_trial_var_names.size());
131 84 : for (const auto i : index_range(_trial_var_names))
132 : {
133 42 : const auto & trial_var_name = _trial_var_names.at(i);
134 42 : mfem::ParComplexGridFunction & trial_gf = *_cmplx_var_ess_constraints.at(i);
135 :
136 : // Make sure we update the size, if this mesh has changed recently for instance
137 42 : trial_gf.Update();
138 :
139 : // Initial guess for iterative solvers (initial condition or the previous time step solution)
140 42 : static_cast<mfem::Vector &>(trial_gf) = _complex_gfuncs->GetRef(trial_var_name);
141 :
142 42 : mfem::Array<int> global_ess_markers(trial_gf.ParFESpace()->GetParMesh()->bdr_attributes.Max());
143 42 : global_ess_markers = 0;
144 : // Set strongly constrained DoFs of trial_gf on essential boundaries and add markers for all
145 : // essential boundaries to the global_ess_markers array
146 42 : ApplyComplexEssentialBC(trial_var_name, trial_gf, global_ess_markers);
147 42 : trial_gf.FESpace()->GetEssentialTrueDofs(global_ess_markers, _ess_tdof_lists.at(i));
148 42 : }
149 42 : }
150 :
151 : void
152 63 : ComplexEquationSystem::AddComplexKernel(std::shared_ptr<MFEMComplexKernel> kernel)
153 : {
154 63 : const auto & trial_var_name = kernel->getTrialVariableName();
155 63 : const auto & test_var_name = kernel->getTestVariableName();
156 63 : AddCoupledVariableNameIfMissing(trial_var_name);
157 63 : AddTestVariableNameIfMissing(test_var_name);
158 : // Register new complex kernels map if not present for the test variable
159 63 : if (!_cmplx_kernels_map.Has(test_var_name))
160 : {
161 : auto kernel_field_map =
162 42 : std::make_shared<NamedFieldsMap<std::vector<std::shared_ptr<MFEMComplexKernel>>>>();
163 42 : _cmplx_kernels_map.Register(test_var_name, std::move(kernel_field_map));
164 42 : }
165 : // Register new complex kernels map if not present for the test/trial variable pair
166 63 : if (!_cmplx_kernels_map.Get(test_var_name)->Has(trial_var_name))
167 : {
168 42 : auto kernels = std::make_shared<std::vector<std::shared_ptr<MFEMComplexKernel>>>();
169 42 : _cmplx_kernels_map.Get(test_var_name)->Register(trial_var_name, std::move(kernels));
170 42 : }
171 63 : _cmplx_kernels_map.GetRef(test_var_name).Get(trial_var_name)->push_back(std::move(kernel));
172 63 : }
173 :
174 : void
175 14 : ComplexEquationSystem::AddComplexIntegratedBC(std::shared_ptr<MFEMComplexIntegratedBC> bc)
176 : {
177 14 : const auto & trial_var_name = bc->getTrialVariableName();
178 14 : const auto & test_var_name = bc->getTestVariableName();
179 14 : AddCoupledVariableNameIfMissing(trial_var_name);
180 14 : AddTestVariableNameIfMissing(test_var_name);
181 : // Register new complex integrated bc map if not present for the test variable
182 14 : if (!_cmplx_integrated_bc_map.Has(test_var_name))
183 : {
184 : auto integrated_bc_field_map =
185 7 : std::make_shared<NamedFieldsMap<std::vector<std::shared_ptr<MFEMComplexIntegratedBC>>>>();
186 7 : _cmplx_integrated_bc_map.Register(test_var_name, std::move(integrated_bc_field_map));
187 7 : }
188 : // Register new complex integrated bc map if not present for the test/trial variable pair
189 14 : if (!_cmplx_integrated_bc_map.Get(test_var_name)->Has(trial_var_name))
190 : {
191 7 : auto bcs = std::make_shared<std::vector<std::shared_ptr<MFEMComplexIntegratedBC>>>();
192 7 : _cmplx_integrated_bc_map.Get(test_var_name)->Register(trial_var_name, std::move(bcs));
193 7 : }
194 14 : _cmplx_integrated_bc_map.GetRef(test_var_name).Get(trial_var_name)->push_back(std::move(bc));
195 14 : }
196 :
197 : void
198 63 : ComplexEquationSystem::AddComplexEssentialBCs(std::shared_ptr<MFEMComplexEssentialBC> bc)
199 : {
200 63 : const auto & test_var_name = bc->getTestVariableName();
201 63 : AddTestVariableNameIfMissing(test_var_name);
202 : // Register new complex essential bc map if not present for the test variable
203 63 : if (!_cmplx_essential_bc_map.Has(test_var_name))
204 : {
205 42 : auto bcs = std::make_shared<std::vector<std::shared_ptr<MFEMComplexEssentialBC>>>();
206 42 : _cmplx_essential_bc_map.Register(test_var_name, std::move(bcs));
207 42 : }
208 63 : _cmplx_essential_bc_map.GetRef(test_var_name).push_back(std::move(bc));
209 63 : }
210 :
211 : void
212 0 : ComplexEquationSystem::FormSystemOperator(mfem::OperatorHandle & op,
213 : mfem::BlockVector & trueX,
214 : mfem::BlockVector & trueRHS)
215 : {
216 0 : auto & test_var_name = _test_var_names.at(0);
217 0 : mfem::Vector aux_x, aux_rhs;
218 0 : mfem::OperatorPtr aux_a;
219 :
220 0 : auto slf = _slfs.Get(test_var_name);
221 0 : slf->FormLinearSystem(_ess_tdof_lists.at(0),
222 0 : *_cmplx_var_ess_constraints.at(0),
223 0 : *_clfs.Get(test_var_name),
224 : aux_a,
225 : aux_x,
226 : aux_rhs,
227 : /*copy_interior=*/true);
228 :
229 0 : trueX.GetBlock(0) = aux_x;
230 0 : trueRHS.GetBlock(0) = aux_rhs;
231 0 : trueX.SyncFromBlocks();
232 0 : trueRHS.SyncFromBlocks();
233 :
234 0 : op.Reset(aux_a.Ptr());
235 0 : aux_a.SetOperatorOwner(false);
236 0 : }
237 :
238 : void
239 42 : ComplexEquationSystem::FormSystemMatrix(mfem::OperatorHandle & op,
240 : mfem::BlockVector & trueX,
241 : mfem::BlockVector & trueRHS)
242 : {
243 :
244 : // Allocate block operator
245 42 : DeleteHBlocks();
246 42 : _h_blocks.SetSize(_test_var_names.size(), _trial_var_names.size());
247 42 : _h_blocks = nullptr;
248 : // Zero out RHS and sync memory
249 42 : trueRHS = 0.0;
250 42 : trueRHS.SyncToBlocks();
251 :
252 : // Form diagonal blocks.
253 84 : for (const auto i : index_range(_test_var_names))
254 : {
255 42 : auto & test_var_name = _test_var_names.at(i);
256 :
257 42 : mfem::Vector aux_x, aux_rhs;
258 42 : mfem::OperatorHandle aux_a;
259 :
260 42 : auto slf = _slfs.Get(test_var_name);
261 42 : slf->FormLinearSystem(_ess_tdof_lists.at(i),
262 42 : *_cmplx_var_ess_constraints.at(i),
263 42 : *_clfs.Get(test_var_name),
264 : aux_a,
265 : aux_x,
266 : aux_rhs,
267 : /*copy_interior=*/true);
268 42 : trueX.GetBlock(i) = aux_x;
269 42 : trueRHS.GetBlock(i) = aux_rhs;
270 42 : _h_blocks(i, i) = aux_a.As<mfem::ComplexHypreParMatrix>()->GetSystemMatrix();
271 42 : }
272 : // Sync memory
273 42 : trueX.SyncFromBlocks();
274 42 : trueRHS.SyncFromBlocks();
275 :
276 : // Create monolithic matrix
277 42 : op.Reset(mfem::HypreParMatrixFromBlocks(_h_blocks));
278 42 : }
279 :
280 : // Equation system Mult
281 : void
282 0 : ComplexEquationSystem::Mult(const mfem::Vector & x, mfem::Vector & residual) const
283 : {
284 0 : _linear_operator->Mult(x, residual);
285 0 : x.HostRead();
286 0 : residual.HostRead();
287 0 : }
288 :
289 : void
290 42 : ComplexEquationSystem::SetTrialVariablesFromTrueVectors(const mfem::BlockVector & trueX) const
291 : {
292 84 : for (const auto i : index_range(_trial_var_names))
293 : {
294 42 : auto & trial_var_name = _trial_var_names.at(i);
295 42 : trueX.GetBlock(i).SyncMemory(trueX);
296 42 : _complex_gfuncs->Get(trial_var_name)->Distribute(&(trueX.GetBlock(i)));
297 : }
298 : // Solution variables changed: stored projections of solution-dependent coefficients are stale.
299 42 : if (_coefficient_manager)
300 42 : _coefficient_manager->markSolutionChanged();
301 42 : }
302 : }
303 :
304 : #endif
|
