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