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 "MFEMProblem.h"
13 : #include "MFEMVariable.h"
14 : #include "MFEMIndicator.h"
15 : #include "MFEMSubMesh.h"
16 : #include "MFEMFunctorMaterial.h"
17 : #include "MFEMExecutedObject.h"
18 : #include "MFEMVectorUtils.h"
19 : #include "MFEMFESpaceHierarchy.h"
20 : #include "Postprocessor.h"
21 : #include "VectorPostprocessor.h"
22 : #include "MFEMNonlinearSolverBase.h"
23 : #include "DependencyResolver.h"
24 : #include "MooseUtils.h"
25 : #include "DataIO.h"
26 :
27 : #include "libmesh/string_to_enum.h"
28 :
29 : #include <vector>
30 : #include <algorithm>
31 : #include <map>
32 : #include <deque>
33 : #include <sstream>
34 :
35 : registerMooseObject("MooseApp", MFEMProblem);
36 :
37 : namespace
38 : {
39 : std::vector<MFEMSolverName>
40 2164 : getMFEMSolverDependencies(const InputParameters & parameters)
41 : {
42 2164 : std::vector<MFEMSolverName> dependencies;
43 :
44 54210 : for (const auto & [param_name, _] : parameters)
45 : {
46 52046 : if (parameters.isPrivate(param_name))
47 32493 : continue;
48 :
49 19553 : if (auto * name = parameters.queryParam<MFEMSolverName>(param_name))
50 1001 : dependencies.push_back(*name);
51 18552 : else if (auto * names = parameters.queryParam<std::vector<MFEMSolverName>>(param_name))
52 11 : dependencies.insert(dependencies.end(), names->begin(), names->end());
53 : }
54 :
55 2164 : return dependencies;
56 0 : }
57 : }
58 :
59 : InputParameters
60 8360 : MFEMProblem::validParams()
61 : {
62 8360 : InputParameters params = ExternalProblem::validParams();
63 16720 : params.addClassDescription("Problem type for building and solving the finite element problem "
64 : "using the MFEM finite element library.");
65 33440 : MooseEnum numeric_types("real complex", "real");
66 25080 : params.addParam<MooseEnum>("numeric_type", numeric_types, "Number type used for the problem");
67 :
68 16720 : return params;
69 8360 : }
70 :
71 2050 : MFEMProblem::MFEMProblem(const InputParameters & params)
72 : : ExternalProblem(params),
73 4100 : _num_type{static_cast<int>(getParam<MooseEnum>("numeric_type"))},
74 2050 : _solution_state_data(declareRestartableDataWithContext<Moose::MFEM::SolutionState>(
75 4100 : "mfem_solution_state", &_problem_data))
76 : {
77 : // Initialise Hypre for all MFEM problems.
78 2050 : mfem::Hypre::Init();
79 : // Disable multithreading for all MFEM problems (including any libMesh or MFEM subapps).
80 2050 : libMesh::libMeshPrivateData::_n_threads = 1;
81 : #ifdef LIBMESH_HAVE_OPENMP
82 2050 : omp_set_num_threads(1);
83 : #endif
84 2050 : setMesh();
85 2050 : }
86 :
87 : void
88 1694 : MFEMProblem::initialSetup()
89 : {
90 1694 : ExternalProblem::initialSetup();
91 :
92 1694 : std::vector<MFEMExecutedObject *> objects;
93 1694 : theWarehouse()
94 1694 : .query()
95 1694 : .condition<AttribSystem>("MFEMExecutedObject")
96 3388 : .condition<AttribThread>(0)
97 1694 : .queryInto(objects);
98 4985 : for (auto * const object : objects)
99 3309 : object->initialSetup();
100 :
101 : // MFEM indicators create their estimators during addIndicator(); markers still need an explicit
102 : // setup pass because they are no longer initialized through the libMesh/MOOSE user-object path.
103 1676 : std::vector<MFEMRefinementMarker *> markers;
104 1676 : theWarehouse().query().condition<AttribSystem>("Marker").queryInto(markers);
105 1702 : for (auto marker : markers)
106 26 : marker->initialSetup();
107 1676 : }
108 :
109 : void
110 11750 : MFEMProblem::execute(const ExecFlagType & exec_type)
111 : {
112 11750 : setCurrentExecuteOnFlag(exec_type);
113 11750 : executeMFEMObjects(exec_type);
114 :
115 11750 : ExternalProblem::execute(exec_type);
116 11750 : }
117 :
118 : void
119 2050 : MFEMProblem::setMesh()
120 : {
121 2050 : auto pmesh = mesh().getMFEMParMeshPtr();
122 2050 : getProblemData().pmesh = pmesh;
123 2050 : getProblemData().comm = pmesh->GetComm();
124 2050 : getProblemData().num_procs = pmesh->GetNRanks();
125 2050 : getProblemData().myid = pmesh->GetMyRank();
126 2050 : }
127 :
128 : void
129 26 : MFEMProblem::addIndicator(const std::string & indicator_type,
130 : const std::string & name,
131 : InputParameters & parameters)
132 : {
133 52 : auto estimator = addObject<MFEMIndicator>(indicator_type, name, parameters).front();
134 :
135 : // construct the estimator itself
136 26 : estimator->createEstimator();
137 26 : }
138 :
139 : void
140 26 : MFEMProblem::addMarker(const std::string & marker_type,
141 : const std::string & name,
142 : InputParameters & parameters)
143 : {
144 52 : getProblemData().refiner = addObject<MFEMRefinementMarker>(marker_type, name, parameters).front();
145 26 : }
146 :
147 : void
148 2164 : MFEMProblem::addMFEMSolver(const std::string & solver_type,
149 : const std::string & name,
150 : InputParameters & parameters)
151 : {
152 : mooseAssert(!_mfem_solver_definitions.count(name), "Multiple MFEM solvers named '" + name + "'.");
153 2164 : _mfem_solver_definitions.emplace(name, MFEMSolverDefinition{solver_type, ¶meters});
154 2164 : }
155 :
156 : void
157 1696 : MFEMProblem::resolveMFEMSolvers()
158 : {
159 1696 : if (_mfem_solver_definitions.empty())
160 593 : return;
161 :
162 1103 : DependencyResolver<std::string> resolver;
163 :
164 3267 : for (auto & [solver_name, definition] : _mfem_solver_definitions)
165 : {
166 2164 : const auto dependencies = getMFEMSolverDependencies(*definition.parameters);
167 2164 : if (dependencies.empty())
168 1163 : resolver.addNode(solver_name);
169 :
170 3176 : for (const auto & dependency_name : dependencies)
171 : {
172 1012 : auto dependency_it = _mfem_solver_definitions.find(dependency_name);
173 1012 : if (dependency_it == _mfem_solver_definitions.end())
174 0 : mooseError("MFEM solver '",
175 : solver_name,
176 : "' references MFEM solver '",
177 : dependency_name,
178 : "', but no solver with that name was provided in the [Solvers] block.");
179 :
180 1012 : dependency_it->second.referenced = true;
181 1012 : resolver.addEdge(dependency_name, solver_name);
182 : }
183 2164 : }
184 :
185 1103 : const std::vector<std::string> * sorted_solver_names = nullptr;
186 : try
187 : {
188 1103 : sorted_solver_names = &resolver.getSortedValues();
189 : }
190 0 : catch (CyclicDependencyException<std::string> & e)
191 : {
192 0 : mooseError("Cyclic MFEM solver dependency detected: ",
193 0 : MooseUtils::join(e.getCyclicDependencies(), " <- "));
194 0 : }
195 :
196 1103 : auto & problem_data = getProblemData();
197 : mooseAssert(!problem_data.jacobian_solver, "MFEM linear solver driver already assigned");
198 : mooseAssert(!problem_data.nonlinear_solver, "MFEM nonlinear solver driver already assigned");
199 :
200 3263 : for (const auto & solver_name : *sorted_solver_names)
201 : {
202 2162 : auto & definition = libmesh_map_find(_mfem_solver_definitions, solver_name);
203 : auto solver =
204 6484 : addObject<Moose::MFEM::SolverBase>(definition.type, solver_name, *definition.parameters)
205 2160 : .front();
206 :
207 2160 : if (definition.referenced)
208 1010 : continue;
209 :
210 1150 : if (auto lin_solver = std::dynamic_pointer_cast<Moose::MFEM::LinearSolverBase>(solver))
211 : {
212 1094 : if (problem_data.jacobian_solver)
213 0 : mooseError("Multiple MFEM linear solver drivers provided. '",
214 0 : problem_data.jacobian_solver->name(),
215 : "' and '",
216 0 : lin_solver->name(),
217 : "' are not referenced by another MFEM solver.");
218 1094 : problem_data.jacobian_solver = lin_solver;
219 : }
220 56 : else if (auto nonlinear_solver =
221 56 : std::dynamic_pointer_cast<Moose::MFEM::NonlinearSolverBase>(solver);
222 56 : nonlinear_solver)
223 : {
224 56 : if (problem_data.nonlinear_solver)
225 0 : mooseError("Multiple MFEM nonlinear solver drivers provided. '",
226 0 : problem_data.nonlinear_solver->name(),
227 : "' and '",
228 0 : nonlinear_solver->name(),
229 : "' are not referenced by another MFEM solver.");
230 56 : problem_data.nonlinear_solver = nonlinear_solver;
231 : }
232 : else
233 0 : mooseError("Unsupported MFEM solver object type '",
234 0 : solver->type(),
235 : "' for solver '",
236 0 : solver->name(),
237 1206 : "'.");
238 2160 : }
239 :
240 1101 : _mfem_solver_definitions.clear();
241 1101 : }
242 :
243 : void
244 1490 : MFEMProblem::addBoundaryCondition(const std::string & bc_name,
245 : const std::string & name,
246 : InputParameters & parameters)
247 : {
248 2980 : auto bc = addObject<MFEMBoundaryCondition>(bc_name, name, parameters).front();
249 1490 : const auto & mfem_bc = *bc;
250 :
251 1490 : if (dynamic_cast<const MFEMIntegratedBC *>(&mfem_bc))
252 : {
253 62 : auto integrated_bc = std::dynamic_pointer_cast<MFEMIntegratedBC>(bc);
254 : auto eqsys =
255 62 : std::dynamic_pointer_cast<Moose::MFEM::EquationSystem>(getProblemData().eqn_system);
256 62 : if (eqsys)
257 62 : eqsys->AddIntegratedBC(std::move(integrated_bc));
258 : else
259 0 : mooseError("Cannot add integrated BC with name '" + name +
260 : "' because there is no corresponding equation system.");
261 62 : }
262 1428 : else if (dynamic_cast<const MFEMComplexIntegratedBC *>(&mfem_bc))
263 : {
264 14 : auto integrated_bc = std::dynamic_pointer_cast<MFEMComplexIntegratedBC>(bc);
265 : auto eqsys =
266 14 : std::dynamic_pointer_cast<Moose::MFEM::ComplexEquationSystem>(getProblemData().eqn_system);
267 14 : if (eqsys)
268 14 : eqsys->AddComplexIntegratedBC(std::move(integrated_bc));
269 : else
270 0 : mooseError("Cannot add complex integrated BC with name '" + name +
271 : "' because there is no corresponding equation system.");
272 14 : }
273 1414 : else if (dynamic_cast<const MFEMComplexEssentialBC *>(&mfem_bc))
274 : {
275 63 : auto essential_bc = std::dynamic_pointer_cast<MFEMComplexEssentialBC>(bc);
276 : auto eqsys =
277 63 : std::dynamic_pointer_cast<Moose::MFEM::ComplexEquationSystem>(getProblemData().eqn_system);
278 63 : if (eqsys)
279 63 : eqsys->AddComplexEssentialBCs(std::move(essential_bc));
280 : else
281 0 : mooseError("Cannot add boundary condition with name '" + name +
282 : "' because there is no corresponding equation system.");
283 63 : }
284 1351 : else if (dynamic_cast<const MFEMEssentialBC *>(&mfem_bc))
285 : {
286 1351 : auto essential_bc = std::dynamic_pointer_cast<MFEMEssentialBC>(bc);
287 : auto eqsys =
288 1351 : std::dynamic_pointer_cast<Moose::MFEM::EquationSystem>(getProblemData().eqn_system);
289 1351 : if (eqsys)
290 1351 : eqsys->AddEssentialBC(std::move(essential_bc));
291 : else
292 0 : mooseError("Cannot add boundary condition with name '" + name +
293 : "' because there is no corresponding equation system.");
294 1351 : }
295 : else
296 : {
297 0 : mooseError("Unsupported bc of type '", bc_name, "' and name '", name, "' detected.");
298 : }
299 1490 : }
300 :
301 : void
302 0 : MFEMProblem::addMaterial(const std::string &, const std::string &, InputParameters &)
303 : {
304 0 : mooseError(
305 : "MFEM materials must be added through the 'FunctorMaterials' block and not 'Materials'");
306 : }
307 :
308 : void
309 378 : MFEMProblem::addFunctorMaterial(const std::string & material_name,
310 : const std::string & name,
311 : InputParameters & parameters)
312 : {
313 772 : addObject<MFEMFunctorMaterial>(material_name, name, parameters);
314 362 : }
315 :
316 : void
317 3009 : MFEMProblem::addFESpace(const std::string & type,
318 : const std::string & name,
319 : InputParameters & parameters)
320 : {
321 3009 : if (getProblemData().fespace_hierarchies.Has(name))
322 0 : mooseError("Cannot add FESpace '",
323 : name,
324 : "': an MFEMFESpaceHierarchy with the same name already exists. "
325 : "FESpaces and FESpaceHierarchies share the fespaces namespace.");
326 :
327 6018 : auto & mfem_fespace = *addObject<MFEMFESpace>(type, name, parameters).front();
328 :
329 : // Register fespace and associated fe collection.
330 3009 : getProblemData().fecs.Register(name, mfem_fespace.getFEC());
331 3009 : getProblemData().fespaces.Register(name, mfem_fespace.getFESpace());
332 3009 : }
333 :
334 : void
335 11 : MFEMProblem::addFESpaceHierarchy(const std::string & type,
336 : const std::string & name,
337 : InputParameters & parameters)
338 : {
339 11 : if (getProblemData().fespaces.Has(name))
340 0 : mooseError("Cannot add MFEMFESpaceHierarchy '",
341 : name,
342 : "': a FESpace with the same name already exists. "
343 : "FESpaces and FESpaceHierarchies share the fespaces namespace.");
344 :
345 22 : auto hierarchy_obj = addObject<MFEMFESpaceHierarchy>(type, name, parameters).front();
346 11 : auto hierarchy_shared = hierarchy_obj->getHierarchyShared();
347 : // Register the hierarchy for co-ownership by solvers.
348 11 : getProblemData().fespace_hierarchies.Register(name, hierarchy_shared);
349 : // Register the finest-level FESpace in fespaces under the hierarchy name so that
350 : // variables can say `fespace = <hierarchy_name>` without a separate FESpace definition.
351 : // The aliasing shared_ptr keeps the hierarchy alive as long as this entry lives.
352 : auto finest = std::shared_ptr<mfem::ParFiniteElementSpace>(
353 11 : hierarchy_shared, &hierarchy_obj->getHierarchy().GetFinestFESpace());
354 11 : getProblemData().fespaces.Register(name, finest);
355 11 : }
356 :
357 : void
358 2117 : MFEMProblem::validateVariableNumericType(const std::string & var_type,
359 : const std::string & var_name) const
360 : {
361 2117 : const bool variable_is_complex = var_type == "MFEMComplexVariable";
362 2117 : const bool problem_is_complex = _num_type == NumericType::COMPLEX;
363 2117 : if (variable_is_complex != problem_is_complex)
364 8 : paramError("numeric_type",
365 : "The problem numeric type does not match primary MFEM variable '",
366 : var_name,
367 : "', which is ",
368 : variable_is_complex ? "complex." : "real.");
369 2113 : }
370 :
371 : void
372 2091 : MFEMProblem::addVariable(const std::string & var_type,
373 : const std::string & var_name,
374 : InputParameters & parameters)
375 : {
376 2091 : validateVariableNumericType(var_type, var_name);
377 2087 : addGridFunction(var_type, var_name, parameters);
378 : // MOOSE variables store DoFs for the trial variable and its time derivatives up to second order;
379 : // MFEM GridFunctions store data for only one set of DoFs each, so we must add additional
380 : // GridFunctions for time derivatives.
381 2087 : if (isTransient())
382 : {
383 : const auto time_derivative_var_name =
384 299 : getMFEMObject<MFEMVariable>("MooseVariableBase", var_name).getTimeDerivativeName();
385 299 : getProblemData().time_derivative_map.addTimeDerivativeAssociation(var_name,
386 : time_derivative_var_name);
387 299 : addGridFunction(var_type, time_derivative_var_name, parameters);
388 299 : }
389 2087 : }
390 :
391 : void
392 4208 : MFEMProblem::addGridFunction(const std::string & var_type,
393 : const std::string & var_name,
394 : InputParameters & parameters)
395 : {
396 :
397 4208 : if (var_type == "MFEMVariable" || var_type == "MFEMComplexVariable")
398 : {
399 : // Add MFEM variable directly.
400 3991 : if (var_type == "MFEMComplexVariable")
401 231 : addObject<MFEMComplexVariable>(var_type, var_name, parameters);
402 : else
403 11742 : addObject<MFEMVariable>(var_type, var_name, parameters);
404 : }
405 : else
406 : {
407 : // Add MOOSE variable.
408 217 : ExternalProblem::addVariable(var_type, var_name, parameters);
409 :
410 : // Add MFEM variable indirectly ("gridfunction").
411 217 : InputParameters mfem_variable_params = addMFEMFESpaceFromMOOSEVariable(parameters);
412 868 : addObject<MFEMVariable>("MFEMVariable", var_name, mfem_variable_params);
413 217 : }
414 :
415 : // Register gridfunction.
416 4208 : if (var_type == "MFEMComplexVariable")
417 : {
418 : MFEMComplexVariable & mfem_variable =
419 77 : getMFEMObject<MFEMComplexVariable>("MooseVariableBase", var_name);
420 77 : getProblemData().cmplx_gridfunctions.Register(var_name, mfem_variable.getComplexGridFunction());
421 77 : mfem_variable.declareCoefficients();
422 : }
423 : else // must be real, but may have been set up indirectly from a MOOSE variable
424 : {
425 4131 : MFEMVariable & mfem_variable = getMFEMObject<MFEMVariable>("MooseVariableBase", var_name);
426 4131 : getProblemData().gridfunctions.Register(var_name, mfem_variable.getGridFunction());
427 4131 : mfem_variable.declareCoefficients();
428 : }
429 4208 : }
430 :
431 : void
432 1666 : MFEMProblem::addAuxVariable(const std::string & var_type,
433 : const std::string & var_name,
434 : InputParameters & parameters)
435 : {
436 : // We handle MFEM AuxVariables just like MFEM Variables, except
437 : // we do not add additional GridFunctions for time derivatives.
438 1666 : addGridFunction(var_type, var_name, parameters);
439 1666 : }
440 :
441 : void
442 542 : MFEMProblem::addAuxKernel(const std::string & kernel_name,
443 : const std::string & name,
444 : InputParameters & parameters)
445 : {
446 1084 : addObject<MFEMExecutedObject>(kernel_name, name, parameters);
447 542 : }
448 :
449 : void
450 2020 : MFEMProblem::addKernel(const std::string & kernel_name,
451 : const std::string & name,
452 : InputParameters & parameters)
453 : {
454 4040 : auto kernel = addObject<MFEMKernel>(kernel_name, name, parameters).front();
455 2020 : const auto & kernel_object = *kernel;
456 :
457 2020 : if (dynamic_cast<const MFEMComplexKernel *>(&kernel_object))
458 : {
459 63 : auto complex_kernel = std::dynamic_pointer_cast<MFEMComplexKernel>(kernel);
460 : auto eqsys =
461 63 : std::dynamic_pointer_cast<Moose::MFEM::ComplexEquationSystem>(getProblemData().eqn_system);
462 63 : if (eqsys)
463 63 : eqsys->AddComplexKernel(std::move(complex_kernel));
464 : else
465 0 : mooseError("Cannot add complex kernel with name '" + name +
466 : "' because there is no corresponding equation system.");
467 63 : }
468 : else
469 : {
470 : auto eqsys =
471 1957 : std::dynamic_pointer_cast<Moose::MFEM::EquationSystem>(getProblemData().eqn_system);
472 1957 : if (eqsys)
473 1957 : eqsys->AddKernel(std::move(kernel));
474 : else
475 0 : mooseError("Cannot add kernel with name '" + name +
476 : "' because there is no corresponding equation system.");
477 1957 : }
478 2020 : }
479 :
480 : void
481 63 : MFEMProblem::addRealComponentToKernel(const std::string & kernel_name,
482 : const std::string & name,
483 : InputParameters & parameters)
484 : {
485 : auto parent_ptr = std::dynamic_pointer_cast<MFEMComplexKernel>(
486 63 : getMFEMObject<MFEMComplexKernel>("Kernel", name).getSharedPtr());
487 252 : parameters.set<VariableName>("variable") = parent_ptr->getParam<VariableName>("variable");
488 126 : auto kernel_ptr = addObject<MFEMKernel>(kernel_name, name + "_real", parameters).front();
489 63 : parent_ptr->setRealKernel(kernel_ptr);
490 63 : }
491 :
492 : void
493 56 : MFEMProblem::addImagComponentToKernel(const std::string & kernel_name,
494 : const std::string & name,
495 : InputParameters & parameters)
496 : {
497 : auto parent_ptr = std::dynamic_pointer_cast<MFEMComplexKernel>(
498 56 : getMFEMObject<MFEMComplexKernel>("Kernel", name).getSharedPtr());
499 224 : parameters.set<VariableName>("variable") = parent_ptr->getParam<VariableName>("variable");
500 112 : auto kernel_ptr = addObject<MFEMKernel>(kernel_name, name + "_imag", parameters).front();
501 56 : parent_ptr->setImagKernel(kernel_ptr);
502 56 : }
503 :
504 : void
505 0 : MFEMProblem::addRealComponentToBC(const std::string & kernel_name,
506 : const std::string & name,
507 : InputParameters & parameters)
508 : {
509 : auto parent_ptr = std::dynamic_pointer_cast<MFEMComplexIntegratedBC>(
510 0 : getMFEMObject<MFEMComplexIntegratedBC>("BoundaryCondition", name).getSharedPtr());
511 0 : parameters.set<VariableName>("variable") = parent_ptr->getParam<VariableName>("variable");
512 0 : parameters.set<std::vector<BoundaryName>>("boundary") =
513 0 : parent_ptr->getParam<std::vector<BoundaryName>>("boundary");
514 : auto bc_ptr = std::dynamic_pointer_cast<MFEMIntegratedBC>(
515 0 : addObject<MFEMBoundaryCondition>(kernel_name, name + "_real", parameters).front());
516 0 : parent_ptr->setRealBC(bc_ptr);
517 0 : }
518 :
519 : void
520 0 : MFEMProblem::addImagComponentToBC(const std::string & kernel_name,
521 : const std::string & name,
522 : InputParameters & parameters)
523 : {
524 : auto parent_ptr = std::dynamic_pointer_cast<MFEMComplexIntegratedBC>(
525 0 : getMFEMObject<MFEMComplexIntegratedBC>("BoundaryCondition", name).getSharedPtr());
526 0 : parameters.set<VariableName>("variable") = parent_ptr->getParam<VariableName>("variable");
527 0 : parameters.set<std::vector<BoundaryName>>("boundary") =
528 0 : parent_ptr->getParam<std::vector<BoundaryName>>("boundary");
529 : auto bc_ptr = std::dynamic_pointer_cast<MFEMIntegratedBC>(
530 0 : addObject<MFEMBoundaryCondition>(kernel_name, name + "_imag", parameters).front());
531 0 : parent_ptr->setImagBC(bc_ptr);
532 0 : }
533 :
534 : int
535 396 : vectorFunctionDim(const std::string & type, const InputParameters & parameters)
536 : {
537 792 : if (parameters.isParamSetByUser("expression_z"))
538 366 : return 3;
539 90 : if (parameters.isParamSetByUser("expression_y") || type == "LevelSetOlssonVortex")
540 28 : return 2;
541 4 : if (parameters.isParamSetByUser("expression_x"))
542 2 : return 1;
543 :
544 0 : return 3;
545 : }
546 :
547 : const std::vector<std::string> SCALAR_FUNCS = {"Axisymmetric2D3DSolutionFunction",
548 : "BicubicSplineFunction",
549 : "CoarsenedPiecewiseLinear",
550 : "CompositeFunction",
551 : "ConstantFunction",
552 : "ImageFunction",
553 : "ParsedFunction",
554 : "ParsedGradFunction",
555 : "PeriodicFunction",
556 : "PiecewiseBilinear",
557 : "PiecewiseConstant",
558 : "PiecewiseConstantFromCSV",
559 : "PiecewiseLinear",
560 : "PiecewiseLinearFromVectorPostprocessor",
561 : "PiecewiseMultiInterpolation",
562 : "PiecewiseMulticonstant",
563 : "SolutionFunction",
564 : "SplineFunction",
565 : "FunctionSeries",
566 : "LevelSetOlssonBubble",
567 : "LevelSetOlssonPlane",
568 : "NearestReporterCoordinatesFunction",
569 : "ParameterMeshFunction",
570 : "ParsedOptimizationFunction",
571 : "FourierNoise",
572 : "MovingPlanarFront",
573 : "MultiControlDrumFunction",
574 : "Grad2ParsedFunction",
575 : "GradParsedFunction",
576 : "ScaledAbsDifferenceDRLRewardFunction",
577 : "CircularAreaHydraulicDiameterFunction",
578 : "CosineHumpFunction",
579 : "CosineTransitionFunction",
580 : "CubicTransitionFunction",
581 : "GeneralizedCircumference",
582 : "PiecewiseFunction",
583 : "TimeRampFunction"},
584 : VECTOR_FUNCS = {"ParsedVectorFunction", "LevelSetOlssonVortex"};
585 :
586 : void
587 1697 : MFEMProblem::addFunction(const std::string & type,
588 : const std::string & name,
589 : InputParameters & parameters)
590 : {
591 1697 : ExternalProblem::addFunction(type, name, parameters);
592 1697 : auto & func = getFunction(name);
593 : // FIXME: Do we want to have optimised versions for when functions
594 : // are only of space or only of time.
595 1697 : if (std::find(SCALAR_FUNCS.begin(), SCALAR_FUNCS.end(), type) != SCALAR_FUNCS.end())
596 : {
597 1144 : getCoefficients().declareScalar<mfem::FunctionCoefficient>(
598 : name,
599 1144 : [&func](const mfem::Vector & p, mfem::real_t t) -> mfem::real_t
600 2327744 : { return func.value(t, Moose::MFEM::libMeshPointFromMFEMVector(p)); });
601 : }
602 553 : else if (std::find(VECTOR_FUNCS.begin(), VECTOR_FUNCS.end(), type) != VECTOR_FUNCS.end())
603 : {
604 396 : int dim = vectorFunctionDim(type, parameters);
605 396 : getCoefficients().declareVector<mfem::VectorFunctionCoefficient>(
606 : name,
607 : dim,
608 396 : [&func, dim](const mfem::Vector & p, mfem::real_t t, mfem::Vector & u)
609 : {
610 : libMesh::RealVectorValue vector_value =
611 2875414 : func.vectorValue(t, Moose::MFEM::libMeshPointFromMFEMVector(p));
612 11183410 : for (int i = 0; i < dim; i++)
613 : {
614 8307996 : u[i] = vector_value(i);
615 : }
616 2875414 : });
617 : }
618 157 : else if ("MFEMParsedFunction" != type)
619 : {
620 2 : mooseWarning("Could not identify whether function ",
621 : type,
622 : " is scalar or vector; no MFEM coefficient object created.");
623 : }
624 1695 : }
625 :
626 : void
627 904 : MFEMProblem::addPostprocessor(const std::string & type,
628 : const std::string & name,
629 : InputParameters & parameters)
630 : {
631 904 : if (parameters.getSystemAttributeName() == "MFEMExecutedObject")
632 : {
633 1600 : checkUserObjectNameCollision(name, "Postprocessor");
634 1600 : addObject<MFEMExecutedObject>(type, name, parameters);
635 800 : const PostprocessorValue & val = getPostprocessorValueByName(name);
636 800 : getCoefficients().declareScalar<mfem::FunctionCoefficient>(
637 802 : name, [&val](const mfem::Vector &) -> mfem::real_t { return val; });
638 : }
639 : else
640 104 : ExternalProblem::addPostprocessor(type, name, parameters);
641 904 : }
642 :
643 : void
644 452 : MFEMProblem::addVectorPostprocessor(const std::string & type,
645 : const std::string & name,
646 : InputParameters & parameters)
647 : {
648 452 : if (parameters.getSystemAttributeName() == "MFEMExecutedObject")
649 : {
650 904 : checkUserObjectNameCollision(name, "VectorPostprocessor");
651 1354 : addObject<MFEMExecutedObject>(type, name, parameters);
652 : }
653 : else
654 0 : ExternalProblem::addVectorPostprocessor(type, name, parameters);
655 450 : }
656 :
657 : InputParameters
658 217 : MFEMProblem::addMFEMFESpaceFromMOOSEVariable(InputParameters & parameters)
659 : {
660 :
661 434 : InputParameters fespace_params = _factory.getValidParams("MFEMGenericFESpace");
662 434 : InputParameters variable_params = _factory.getValidParams("MFEMVariable");
663 :
664 217 : const auto family = Utility::string_to_enum<FEFamily>(parameters.get<MooseEnum>("family"));
665 217 : auto order = static_cast<int>(parameters.get<MooseEnum>("order"));
666 217 : const auto dim = mesh().dimension();
667 :
668 217 : std::string space;
669 217 : int vdim = 1;
670 :
671 217 : switch (family)
672 : {
673 35 : case FEFamily::LAGRANGE:
674 35 : space = "H1";
675 35 : break;
676 35 : case FEFamily::NEDELEC_ONE:
677 35 : space = "ND";
678 35 : break;
679 35 : case FEFamily::RAVIART_THOMAS:
680 35 : space = "RT";
681 35 : --order;
682 35 : break;
683 42 : case FEFamily::MONOMIAL:
684 : case FEFamily::L2_LAGRANGE:
685 42 : space = "L2";
686 42 : break;
687 28 : case FEFamily::LAGRANGE_VEC:
688 28 : space = "H1";
689 28 : vdim = dim;
690 28 : break;
691 42 : case FEFamily::MONOMIAL_VEC:
692 : case FEFamily::L2_LAGRANGE_VEC:
693 42 : space = "L2";
694 42 : vdim = dim;
695 42 : break;
696 0 : default:
697 0 : mooseError("Unable to set MFEM FESpace for MOOSE variable");
698 : break;
699 : }
700 :
701 : // Create fespace name. If this already exists, we will reuse this for
702 : // the mfem variable ("gridfunction"). If using AMR, this implies all
703 : // variables sharing the fespace are affected.
704 217 : const auto fec_name = space + "_" + std::to_string(dim) + "D_P" + std::to_string(order);
705 217 : const auto fes_name = fec_name + "_X" + std::to_string(vdim);
706 :
707 : // Set all fespace parameters.
708 217 : fespace_params.set<std::string>("fec_name") = fec_name;
709 651 : fespace_params.set<int>("vdim") = vdim;
710 :
711 434 : if (!hasMFEMObject("MFEMFESpace", fes_name))
712 266 : addFESpace("MFEMGenericFESpace", fes_name, fespace_params);
713 :
714 651 : variable_params.set<MFEMFESpaceName>("fespace") = fes_name;
715 :
716 434 : return variable_params;
717 217 : }
718 :
719 : void
720 3272 : MFEMProblem::displaceMesh()
721 : {
722 : // Displace mesh
723 3272 : if (mesh().shouldDisplace())
724 : {
725 38 : mesh().displace(static_cast<mfem::GridFunction const &>(*getMeshDisplacementGridFunction()));
726 : // TODO: update FESpaces GridFunctions etc for transient solves
727 : }
728 3272 : }
729 :
730 : std::optional<std::reference_wrapper<mfem::ParGridFunction const>>
731 1732 : MFEMProblem::getMeshDisplacementGridFunction()
732 : {
733 : // If C++23 transform were available this would be easier
734 1732 : auto const displacement_variable = mesh().getMeshDisplacementVariable();
735 1732 : if (displacement_variable)
736 : {
737 70 : return *_problem_data.gridfunctions.Get(displacement_variable.value());
738 : }
739 : else
740 : {
741 1662 : return std::nullopt;
742 : }
743 : }
744 :
745 : void
746 0 : MFEMProblem::rebalanceMesh(mfem::ParMesh & pmesh)
747 : {
748 0 : if (pmesh.Nonconforming())
749 : {
750 0 : pmesh.Rebalance();
751 0 : updateFESpaces();
752 0 : updateGridFunctions();
753 : }
754 0 : }
755 :
756 : void
757 13 : MFEMProblem::updateFESpaces()
758 : {
759 39 : for (const auto & fe_space_pair : _problem_data.fespaces)
760 26 : fe_space_pair.second->Update();
761 13 : }
762 :
763 : void
764 26 : MFEMProblem::updateGridFunctions()
765 : {
766 78 : for (const auto & gridfunction_pair : _problem_data.gridfunctions)
767 52 : gridfunction_pair.second->Update();
768 26 : }
769 :
770 : std::vector<VariableName>
771 0 : MFEMProblem::getAuxVariableNames()
772 : {
773 0 : return systemBaseAuxiliary().getVariableNames();
774 : }
775 :
776 : MFEMMesh &
777 75549 : MFEMProblem::mesh()
778 : {
779 : mooseAssert(ExternalProblem::mesh().type() == "MFEMMesh",
780 : "Please choose the MFEMMesh mesh type for an MFEMProblem\n");
781 75549 : return static_cast<MFEMMesh &>(_mesh);
782 : }
783 :
784 : const MFEMMesh &
785 3693 : MFEMProblem::mesh() const
786 : {
787 3693 : return const_cast<MFEMProblem *>(this)->mesh();
788 : }
789 :
790 : void
791 206 : MFEMProblem::addSubMesh(const std::string & var_type,
792 : const std::string & var_name,
793 : InputParameters & parameters)
794 : {
795 412 : auto & mfem_submesh = *addObject<MFEMSubMesh>(var_type, var_name, parameters).front();
796 : // Register submesh.
797 206 : getProblemData().submeshes.Register(var_name, mfem_submesh.getSubMesh());
798 206 : }
799 :
800 : void
801 36 : MFEMProblem::addQuadratureFunction(const std::string & type,
802 : const std::string & name,
803 : InputParameters & parameters)
804 : {
805 : // The object declares its coefficient with the CoefficientManager on construction.
806 72 : addObject<MFEMObject>(type, name, parameters);
807 36 : }
808 :
809 : void
810 688 : MFEMProblem::addTransfer(const std::string & transfer_name,
811 : const std::string & name,
812 : InputParameters & parameters)
813 : {
814 688 : if (parameters.getBase() == "MFEMSubMeshTransfer")
815 747 : addObject<MFEMExecutedObject>(transfer_name, name, parameters);
816 : else
817 439 : ExternalProblem::addTransfer(transfer_name, name, parameters);
818 688 : }
819 :
820 : void
821 1272 : MFEMProblem::addInitialCondition(const std::string & ic_name,
822 : const std::string & name,
823 : InputParameters & parameters)
824 : {
825 2544 : addObject<MFEMExecutedObject>(ic_name, name, parameters);
826 1272 : }
827 :
828 : void
829 11756 : MFEMProblem::executeMFEMObjects(const ExecFlagType & exec_type)
830 : {
831 11756 : std::vector<MFEMExecutedObject *> objects;
832 11756 : theWarehouse()
833 11756 : .query()
834 11756 : .condition<AttribSystem>("MFEMExecutedObject")
835 11756 : .condition<AttribExecOns>(exec_type)
836 23512 : .condition<AttribThread>(0)
837 11756 : .queryInto(objects);
838 :
839 11756 : std::map<std::string, const MFEMExecutedObject *> suppliers;
840 16382 : for (auto * const object : objects)
841 9254 : for (const auto & item : object->getSuppliedItems())
842 : {
843 4628 : const auto [it, inserted] = suppliers.emplace(item, object);
844 4628 : if (!inserted && it->second != object)
845 2 : mooseError("MFEM executed-object dependency ambiguity on ",
846 : exec_type,
847 : ": both '",
848 2 : it->second->name(),
849 : "' and '",
850 2 : object->name(),
851 : "' supply '",
852 : item,
853 : "'.");
854 : }
855 :
856 16378 : for (auto * const object : objects)
857 : {
858 4624 : object->initialize();
859 4624 : object->execute();
860 4624 : object->finalize();
861 :
862 4624 : if (auto * const pp = dynamic_cast<const Postprocessor *>(object))
863 : {
864 1328 : _reporter_data.finalize(pp->PPName());
865 1328 : setPostprocessorValueByName(pp->PPName(), pp->getValue());
866 : }
867 :
868 4624 : if (auto * const vpp = dynamic_cast<VectorPostprocessor *>(object))
869 1096 : _reporter_data.finalize(vpp->PPName());
870 : }
871 11758 : }
872 :
873 : std::string
874 1676 : MFEMProblem::solverTypeString(const unsigned int libmesh_dbg_var(solver_sys_num))
875 : {
876 : mooseAssert(solver_sys_num == 0, "No support for multi-system with MFEM right now");
877 :
878 1676 : std::vector<std::string> solvers;
879 :
880 1676 : if (getProblemData().nonlinear_solver)
881 56 : solvers.push_back(MooseUtils::prettyCppType(getProblemData().nonlinear_solver.get()));
882 :
883 1676 : if (getProblemData().jacobian_solver)
884 : {
885 1094 : solvers.push_back(MooseUtils::prettyCppType(getProblemData().jacobian_solver.get()));
886 1094 : if (const auto * prec = getProblemData().jacobian_solver->GetPreconditioner())
887 977 : solvers.push_back(MooseUtils::prettyCppType(prec));
888 : }
889 :
890 6129 : return solvers.empty() ? "None" : MooseUtils::stringJoin(solvers);
891 1676 : }
892 :
893 : bool
894 217 : MFEMProblem::hasMFEMObject(const std::string & system, const std::string & name) const
895 : {
896 217 : std::vector<MooseObject *> objs;
897 217 : theWarehouse()
898 217 : .query()
899 217 : .condition<AttribSystem>(system)
900 434 : .condition<AttribThread>(0)
901 217 : .condition<AttribName>(name)
902 217 : .queryInto(objs);
903 434 : return !objs.empty();
904 217 : }
905 :
906 : #endif
|
