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