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