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 "MFEMProblem.h"
13 : #include "MFEMInitialCondition.h"
14 : #include "MFEMVariable.h"
15 : #include "MFEMSubMesh.h"
16 : #include "MFEMFunctorMaterial.h"
17 : #include "libmesh/string_to_enum.h"
18 :
19 : #include <vector>
20 : #include <algorithm>
21 :
22 : registerMooseObject("MooseApp", MFEMProblem);
23 :
24 : InputParameters
25 9158 : MFEMProblem::validParams()
26 : {
27 9158 : InputParameters params = ExternalProblem::validParams();
28 9158 : params.addClassDescription("Problem type for building and solving finite element problem using"
29 : " the MFEM finite element library.");
30 9158 : return params;
31 0 : }
32 :
33 264 : MFEMProblem::MFEMProblem(const InputParameters & params) : ExternalProblem(params)
34 : {
35 : // Initialise Hypre for all MFEM problems.
36 264 : mfem::Hypre::Init();
37 264 : }
38 :
39 : void
40 134 : MFEMProblem::initialSetup()
41 : {
42 134 : FEProblemBase::initialSetup();
43 134 : addMFEMNonlinearSolver();
44 134 : }
45 :
46 : void
47 134 : MFEMProblem::setMesh()
48 : {
49 134 : auto pmesh = mesh().getMFEMParMeshPtr();
50 134 : getProblemData().pmesh = pmesh;
51 134 : getProblemData().comm = pmesh->GetComm();
52 134 : MPI_Comm_size(pmesh->GetComm(), &(getProblemData().num_procs));
53 134 : MPI_Comm_rank(pmesh->GetComm(), &(getProblemData().myid));
54 134 : }
55 :
56 : void
57 134 : MFEMProblem::initProblemOperator()
58 : {
59 134 : setMesh();
60 134 : auto mfem_exec_ptr = dynamic_cast<MFEMExecutioner *>(_app.getExecutioner());
61 134 : if (mfem_exec_ptr != nullptr)
62 : {
63 134 : mfem_exec_ptr->constructProblemOperator();
64 : }
65 : else
66 : {
67 0 : mooseError("Executioner used that is not currently supported by MFEMProblem");
68 : }
69 134 : }
70 :
71 : void
72 136 : MFEMProblem::addMFEMPreconditioner(const std::string & user_object_name,
73 : const std::string & name,
74 : InputParameters & parameters)
75 : {
76 136 : FEProblemBase::addUserObject(user_object_name, name, parameters);
77 136 : }
78 :
79 : void
80 102 : MFEMProblem::addMFEMSolver(const std::string & user_object_name,
81 : const std::string & name,
82 : InputParameters & parameters)
83 : {
84 102 : FEProblemBase::addUserObject(user_object_name, name, parameters);
85 102 : auto object_ptr = getUserObject<MFEMSolverBase>(name).getSharedPtr();
86 :
87 102 : getProblemData().jacobian_solver = std::dynamic_pointer_cast<MFEMSolverBase>(object_ptr);
88 102 : }
89 :
90 : void
91 134 : MFEMProblem::addMFEMNonlinearSolver()
92 : {
93 134 : auto nl_solver = std::make_shared<mfem::NewtonSolver>(getProblemData().comm);
94 :
95 : // Defaults to one iteration, without further nonlinear iterations
96 134 : nl_solver->SetRelTol(0.0);
97 134 : nl_solver->SetAbsTol(0.0);
98 134 : nl_solver->SetMaxIter(1);
99 :
100 134 : getProblemData().nonlinear_solver = nl_solver;
101 134 : }
102 :
103 : void
104 170 : MFEMProblem::addBoundaryCondition(const std::string & bc_name,
105 : const std::string & name,
106 : InputParameters & parameters)
107 : {
108 170 : FEProblemBase::addUserObject(bc_name, name, parameters);
109 170 : const UserObject * mfem_bc_uo = &(getUserObjectBase(name));
110 170 : if (dynamic_cast<const MFEMIntegratedBC *>(mfem_bc_uo) != nullptr)
111 : {
112 8 : auto object_ptr = getUserObject<MFEMIntegratedBC>(name).getSharedPtr();
113 8 : auto bc = std::dynamic_pointer_cast<MFEMIntegratedBC>(object_ptr);
114 8 : if (getProblemData().eqn_system)
115 : {
116 8 : getProblemData().eqn_system->AddIntegratedBC(std::move(bc));
117 : }
118 : else
119 : {
120 0 : mooseError("Cannot add integrated BC with name '" + name +
121 : "' because there is no corresponding equation system.");
122 : }
123 8 : }
124 162 : else if (dynamic_cast<const MFEMEssentialBC *>(mfem_bc_uo) != nullptr)
125 : {
126 162 : auto object_ptr = getUserObject<MFEMEssentialBC>(name).getSharedPtr();
127 162 : auto mfem_bc = std::dynamic_pointer_cast<MFEMEssentialBC>(object_ptr);
128 162 : if (getProblemData().eqn_system)
129 : {
130 162 : getProblemData().eqn_system->AddEssentialBC(std::move(mfem_bc));
131 : }
132 : else
133 : {
134 0 : mooseError("Cannot add boundary condition with name '" + name +
135 : "' because there is no corresponding equation system.");
136 : }
137 162 : }
138 : else
139 : {
140 0 : mooseError("Unsupported bc of type '", bc_name, "' and name '", name, "' detected.");
141 : }
142 170 : }
143 :
144 : void
145 0 : MFEMProblem::addMaterial(const std::string &, const std::string &, InputParameters &)
146 : {
147 0 : mooseError(
148 : "MFEM materials must be added through the 'FunctorMaterials' block and not 'Materials'");
149 : }
150 :
151 : void
152 57 : MFEMProblem::addFunctorMaterial(const std::string & material_name,
153 : const std::string & name,
154 : InputParameters & parameters)
155 : {
156 57 : FEProblemBase::addUserObject(material_name, name, parameters);
157 49 : getUserObject<MFEMFunctorMaterial>(name);
158 49 : }
159 :
160 : void
161 212 : MFEMProblem::addFESpace(const std::string & user_object_name,
162 : const std::string & name,
163 : InputParameters & parameters)
164 : {
165 212 : FEProblemBase::addUserObject(user_object_name, name, parameters);
166 212 : MFEMFESpace & mfem_fespace(getUserObject<MFEMFESpace>(name));
167 :
168 : // Register fespace and associated fe collection.
169 212 : getProblemData().fecs.Register(name, mfem_fespace.getFEC());
170 212 : getProblemData().fespaces.Register(name, mfem_fespace.getFESpace());
171 212 : }
172 :
173 : void
174 210 : MFEMProblem::addVariable(const std::string & var_type,
175 : const std::string & var_name,
176 : InputParameters & parameters)
177 : {
178 210 : addGridFunction(var_type, var_name, parameters);
179 : // MOOSE variables store DoFs for the trial variable and its time derivatives up to second order;
180 : // MFEM GridFunctions store data for only one set of DoFs each, so we must add additional
181 : // GridFunctions for time derivatives.
182 210 : if (isTransient())
183 : {
184 40 : addGridFunction(var_type, Moose::MFEM::GetTimeDerivativeName(var_name), parameters);
185 : }
186 210 : }
187 :
188 : void
189 250 : MFEMProblem::addGridFunction(const std::string & var_type,
190 : const std::string & var_name,
191 : InputParameters & parameters)
192 : {
193 250 : if (var_type == "MFEMVariable")
194 : {
195 : // Add MFEM variable directly.
196 218 : FEProblemBase::addUserObject(var_type, var_name, parameters);
197 : }
198 : else
199 : {
200 : // Add MOOSE variable.
201 32 : FEProblemBase::addVariable(var_type, var_name, parameters);
202 :
203 : // Add MFEM variable indirectly ("gridfunction").
204 32 : InputParameters mfem_variable_params = addMFEMFESpaceFromMOOSEVariable(parameters);
205 32 : FEProblemBase::addUserObject("MFEMVariable", var_name, mfem_variable_params);
206 32 : }
207 :
208 : // Register gridfunction.
209 250 : MFEMVariable & mfem_variable = getUserObject<MFEMVariable>(var_name);
210 250 : getProblemData().gridfunctions.Register(var_name, mfem_variable.getGridFunction());
211 250 : if (mfem_variable.getFESpace().isScalar())
212 175 : getCoefficients().declareScalar<mfem::GridFunctionCoefficient>(
213 175 : var_name, mfem_variable.getGridFunction().get());
214 : else
215 75 : getCoefficients().declareVector<mfem::VectorGridFunctionCoefficient>(
216 75 : var_name, mfem_variable.getGridFunction().get());
217 250 : }
218 :
219 : void
220 74 : MFEMProblem::addAuxVariable(const std::string & var_type,
221 : const std::string & var_name,
222 : InputParameters & parameters)
223 : {
224 : // We do not handle MFEM AuxVariables separately from variables currently
225 74 : addVariable(var_type, var_name, parameters);
226 74 : }
227 :
228 : void
229 40 : MFEMProblem::addAuxKernel(const std::string & kernel_name,
230 : const std::string & name,
231 : InputParameters & parameters)
232 : {
233 40 : FEProblemBase::addUserObject(kernel_name, name, parameters);
234 40 : }
235 :
236 : void
237 162 : MFEMProblem::addKernel(const std::string & kernel_name,
238 : const std::string & name,
239 : InputParameters & parameters)
240 : {
241 162 : FEProblemBase::addUserObject(kernel_name, name, parameters);
242 162 : const UserObject * kernel_uo = &(getUserObjectBase(name));
243 :
244 162 : if (dynamic_cast<const MFEMKernel *>(kernel_uo) != nullptr)
245 : {
246 162 : auto object_ptr = getUserObject<MFEMKernel>(name).getSharedPtr();
247 162 : auto kernel = std::dynamic_pointer_cast<MFEMKernel>(object_ptr);
248 162 : if (getProblemData().eqn_system)
249 : {
250 162 : getProblemData().eqn_system->AddKernel(std::move(kernel));
251 : }
252 : else
253 : {
254 0 : mooseError("Cannot add kernel with name '" + name +
255 : "' because there is no corresponding equation system.");
256 : }
257 162 : }
258 : else
259 : {
260 0 : mooseError("Unsupported kernel of type '", kernel_name, "' and name '", name, "' detected.");
261 : }
262 162 : }
263 :
264 : libMesh::Point
265 1039138 : pointFromMFEMVector(const mfem::Vector & vec)
266 : {
267 : return libMesh::Point(
268 1039138 : vec.Elem(0), vec.Size() > 1 ? vec.Elem(1) : 0., vec.Size() > 2 ? vec.Elem(2) : 0.);
269 : }
270 :
271 : int
272 64 : vectorFunctionDim(const std::string & type, const InputParameters & parameters)
273 : {
274 64 : if (type == "LevelSetOlssonVortex")
275 : {
276 0 : return 2;
277 : }
278 64 : else if (type == "ParsedVectorFunction")
279 : {
280 64 : if (parameters.isParamSetByUser("expression_z") || parameters.isParamSetByUser("value_z"))
281 : {
282 56 : return 3;
283 : }
284 8 : else if (parameters.isParamSetByUser("expression_y") || parameters.isParamSetByUser("value_y"))
285 : {
286 7 : return 2;
287 : }
288 : else
289 : {
290 1 : return 1;
291 : }
292 : }
293 : else
294 : {
295 0 : return 3;
296 : }
297 : }
298 :
299 : const std::vector<std::string> SCALAR_FUNCS = {"Axisymmetric2D3DSolutionFunction",
300 : "BicubicSplineFunction",
301 : "CoarsenedPiecewiseLinear",
302 : "CompositeFunction",
303 : "ConstantFunction",
304 : "ImageFunction",
305 : "ParsedFunction",
306 : "ParsedGradFunction",
307 : "PeriodicFunction",
308 : "PiecewiseBilinear",
309 : "PiecewiseConstant",
310 : "PiecewiseConstantFromCSV",
311 : "PiecewiseLinear",
312 : "PiecewiseLinearFromVectorPostprocessor",
313 : "PiecewiseMultiInterpolation",
314 : "PiecewiseMulticonstant",
315 : "SolutionFunction",
316 : "SplineFunction",
317 : "FunctionSeries",
318 : "LevelSetOlssonBubble",
319 : "LevelSetOlssonPlane",
320 : "NearestReporterCoordinatesFunction",
321 : "ParameterMeshFunction",
322 : "ParsedOptimizationFunction",
323 : "FourierNoise",
324 : "MovingPlanarFront",
325 : "MultiControlDrumFunction",
326 : "Grad2ParsedFunction",
327 : "GradParsedFunction",
328 : "RichardsExcavGeom",
329 : "ScaledAbsDifferenceDRLRewardFunction",
330 : "CircularAreaHydraulicDiameterFunction",
331 : "CosineHumpFunction",
332 : "CosineTransitionFunction",
333 : "CubicTransitionFunction",
334 : "GeneralizedCircumference",
335 : "PiecewiseFunction",
336 : "TimeRampFunction"},
337 : VECTOR_FUNCS = {"ParsedVectorFunction", "LevelSetOlssonVortex"};
338 :
339 : void
340 125 : MFEMProblem::addFunction(const std::string & type,
341 : const std::string & name,
342 : InputParameters & parameters)
343 : {
344 125 : ExternalProblem::addFunction(type, name, parameters);
345 125 : auto & func = getFunction(name);
346 : // FIXME: Do we want to have optimised versions for when functions
347 : // are only of space or only of time.
348 125 : if (std::find(SCALAR_FUNCS.begin(), SCALAR_FUNCS.end(), type) != SCALAR_FUNCS.end())
349 : {
350 60 : getCoefficients().declareScalar<mfem::FunctionCoefficient>(
351 : name,
352 60 : [&func](const mfem::Vector & p, double t) -> mfem::real_t
353 642099 : { return func.value(t, pointFromMFEMVector(p)); });
354 : }
355 65 : else if (std::find(VECTOR_FUNCS.begin(), VECTOR_FUNCS.end(), type) != VECTOR_FUNCS.end())
356 : {
357 64 : int dim = vectorFunctionDim(type, parameters);
358 64 : getCoefficients().declareVector<mfem::VectorFunctionCoefficient>(
359 : name,
360 : dim,
361 64 : [&func, dim](const mfem::Vector & p, double t, mfem::Vector & u)
362 : {
363 397039 : libMesh::RealVectorValue vector_value = func.vectorValue(t, pointFromMFEMVector(p));
364 1476889 : for (int i = 0; i < dim; i++)
365 : {
366 1079850 : u[i] = vector_value(i);
367 : }
368 397039 : });
369 : }
370 : else
371 : {
372 1 : mooseWarning("Could not identify whether function ",
373 : type,
374 : " is scalar or vector; no MFEM coefficient object created.");
375 : }
376 124 : }
377 :
378 : void
379 13 : MFEMProblem::addPostprocessor(const std::string & type,
380 : const std::string & name,
381 : InputParameters & parameters)
382 : {
383 : // For some reason this isn't getting called
384 13 : ExternalProblem::addPostprocessor(type, name, parameters);
385 13 : const PostprocessorValue & val = getPostprocessorValueByName(name);
386 13 : getCoefficients().declareScalar<mfem::FunctionCoefficient>(
387 14 : name, [&val](const mfem::Vector &, double) -> mfem::real_t { return val; });
388 13 : }
389 :
390 : InputParameters
391 32 : MFEMProblem::addMFEMFESpaceFromMOOSEVariable(InputParameters & parameters)
392 : {
393 :
394 32 : InputParameters fespace_params = _factory.getValidParams("MFEMGenericFESpace");
395 32 : InputParameters mfem_variable_params = _factory.getValidParams("MFEMVariable");
396 :
397 : auto moose_fe_type =
398 64 : FEType(Utility::string_to_enum<Order>(parameters.get<MooseEnum>("order")),
399 96 : Utility::string_to_enum<FEFamily>(parameters.get<MooseEnum>("family")));
400 :
401 32 : std::string mfem_family;
402 32 : int mfem_vdim = 1;
403 :
404 32 : switch (moose_fe_type.family)
405 : {
406 8 : case FEFamily::LAGRANGE:
407 8 : mfem_family = "H1";
408 8 : mfem_vdim = 1;
409 8 : break;
410 8 : case FEFamily::LAGRANGE_VEC:
411 8 : mfem_family = "H1";
412 8 : mfem_vdim = 3;
413 8 : break;
414 8 : case FEFamily::MONOMIAL:
415 8 : mfem_family = "L2";
416 8 : mfem_vdim = 1;
417 8 : break;
418 8 : case FEFamily::MONOMIAL_VEC:
419 8 : mfem_family = "L2";
420 8 : mfem_vdim = 3;
421 8 : break;
422 0 : default:
423 0 : mooseError("Unable to set MFEM FESpace for MOOSE variable");
424 : break;
425 : }
426 :
427 : // Create fespace name. If this already exists, we will reuse this for
428 : // the mfem variable ("gridfunction").
429 64 : const std::string fespace_name = mfem_family + "_" +
430 128 : std::to_string(mesh().getMFEMParMesh().Dimension()) + "D_P" +
431 128 : std::to_string(moose_fe_type.order.get_order());
432 :
433 : // Set all fespace parameters.
434 32 : fespace_params.set<std::string>("fec_name") = fespace_name;
435 32 : fespace_params.set<int>("vdim") = mfem_vdim;
436 :
437 32 : if (!hasUserObject(fespace_name)) // Create the fespace (implicit).
438 : {
439 8 : addFESpace("MFEMGenericFESpace", fespace_name, fespace_params);
440 : }
441 :
442 32 : mfem_variable_params.set<UserObjectName>("fespace") = fespace_name;
443 :
444 64 : return mfem_variable_params;
445 32 : }
446 :
447 : void
448 110 : MFEMProblem::displaceMesh()
449 : {
450 : // Displace mesh
451 110 : if (mesh().shouldDisplace())
452 : {
453 8 : mesh().displace(static_cast<mfem::GridFunction const &>(*getMeshDisplacementGridFunction()));
454 : // TODO: update FESpaces GridFunctions etc for transient solves
455 : }
456 110 : }
457 :
458 : std::optional<std::reference_wrapper<mfem::ParGridFunction const>>
459 142 : MFEMProblem::getMeshDisplacementGridFunction()
460 : {
461 : // If C++23 transform were available this would be easier
462 142 : auto const displacement_variable = mesh().getMeshDisplacementVariable();
463 142 : if (displacement_variable)
464 : {
465 16 : return *_problem_data.gridfunctions.Get(displacement_variable.value());
466 : }
467 : else
468 : {
469 126 : return std::nullopt;
470 : }
471 : }
472 :
473 : std::vector<VariableName>
474 0 : MFEMProblem::getAuxVariableNames()
475 : {
476 0 : return systemBaseAuxiliary().getVariableNames();
477 : }
478 :
479 : MFEMMesh &
480 3706 : MFEMProblem::mesh()
481 : {
482 : mooseAssert(ExternalProblem::mesh().type() == "MFEMMesh",
483 : "Please choose the MFEMMesh mesh type for an MFEMProblem\n");
484 3706 : return static_cast<MFEMMesh &>(_mesh);
485 : }
486 :
487 : const MFEMMesh &
488 506 : MFEMProblem::mesh() const
489 : {
490 506 : return const_cast<MFEMProblem *>(this)->mesh();
491 : }
492 :
493 : void
494 12 : MFEMProblem::addSubMesh(const std::string & var_type,
495 : const std::string & var_name,
496 : InputParameters & parameters)
497 : {
498 : // Add MFEM SubMesh.
499 12 : FEProblemBase::addUserObject(var_type, var_name, parameters);
500 : // Register submesh.
501 12 : MFEMSubMesh & mfem_submesh = getUserObject<MFEMSubMesh>(var_name);
502 12 : getProblemData().submeshes.Register(var_name, mfem_submesh.getSubMesh());
503 12 : }
504 :
505 : void
506 22 : MFEMProblem::addTransfer(const std::string & transfer_name,
507 : const std::string & name,
508 : InputParameters & parameters)
509 : {
510 22 : if (parameters.get<std::string>("_moose_base") == "MFEMSubMeshTransfer")
511 4 : FEProblemBase::addUserObject(transfer_name, name, parameters);
512 : else
513 18 : FEProblemBase::addTransfer(transfer_name, name, parameters);
514 22 : }
515 :
516 : std::shared_ptr<mfem::ParGridFunction>
517 50 : MFEMProblem::getGridFunction(const std::string & name)
518 : {
519 50 : return getUserObject<MFEMVariable>(name).getGridFunction();
520 : }
521 :
522 : void
523 50 : MFEMProblem::addInitialCondition(const std::string & ic_name,
524 : const std::string & name,
525 : InputParameters & parameters)
526 : {
527 50 : FEProblemBase::addUserObject(ic_name, name, parameters);
528 50 : getUserObject<MFEMInitialCondition>(name); // error check
529 50 : }
530 :
531 : std::string
532 134 : MFEMProblem::solverTypeString(const unsigned int libmesh_dbg_var(solver_sys_num))
533 : {
534 : mooseAssert(solver_sys_num == 0, "No support for multi-system with MFEM right now");
535 134 : return MooseUtils::prettyCppType(getProblemData().jacobian_solver.get());
536 : }
537 :
538 : #endif
|
