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 10334 : MFEMProblem::validParams()
26 : {
27 10334 : InputParameters params = ExternalProblem::validParams();
28 10334 : params.addClassDescription("Problem type for building and solving finite element problem using"
29 : " the MFEM finite element library.");
30 10334 : 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 611 : setMesh();
38 611 : }
39 :
40 : void
41 349 : MFEMProblem::initialSetup()
42 : {
43 349 : FEProblemBase::initialSetup();
44 349 : addMFEMNonlinearSolver();
45 349 : }
46 :
47 : void
48 611 : MFEMProblem::setMesh()
49 : {
50 611 : auto pmesh = mesh().getMFEMParMeshPtr();
51 611 : getProblemData().pmesh = pmesh;
52 611 : getProblemData().comm = pmesh->GetComm();
53 611 : MPI_Comm_size(pmesh->GetComm(), &(getProblemData().num_procs));
54 611 : MPI_Comm_rank(pmesh->GetComm(), &(getProblemData().myid));
55 611 : }
56 :
57 : void
58 353 : MFEMProblem::addMFEMPreconditioner(const std::string & user_object_name,
59 : const std::string & name,
60 : InputParameters & parameters)
61 : {
62 353 : FEProblemBase::addUserObject(user_object_name, name, parameters);
63 353 : }
64 :
65 : void
66 301 : MFEMProblem::addMFEMSolver(const std::string & user_object_name,
67 : const std::string & name,
68 : InputParameters & parameters)
69 : {
70 301 : FEProblemBase::addUserObject(user_object_name, name, parameters);
71 301 : auto object_ptr = getUserObject<MFEMSolverBase>(name).getSharedPtr();
72 :
73 301 : getProblemData().jacobian_solver = std::dynamic_pointer_cast<MFEMSolverBase>(object_ptr);
74 301 : }
75 :
76 : void
77 349 : MFEMProblem::addMFEMNonlinearSolver()
78 : {
79 349 : auto nl_solver = std::make_shared<mfem::NewtonSolver>(getProblemData().comm);
80 :
81 : // Defaults to one iteration, without further nonlinear iterations
82 349 : nl_solver->SetRelTol(0.0);
83 349 : nl_solver->SetAbsTol(0.0);
84 349 : nl_solver->SetMaxIter(1);
85 :
86 349 : getProblemData().nonlinear_solver = nl_solver;
87 349 : }
88 :
89 : void
90 476 : MFEMProblem::addBoundaryCondition(const std::string & bc_name,
91 : const std::string & name,
92 : InputParameters & parameters)
93 : {
94 476 : FEProblemBase::addUserObject(bc_name, name, parameters);
95 476 : const UserObject * mfem_bc_uo = &(getUserObjectBase(name));
96 476 : 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 446 : else if (dynamic_cast<const MFEMEssentialBC *>(mfem_bc_uo) != nullptr)
111 : {
112 446 : auto object_ptr = getUserObject<MFEMEssentialBC>(name).getSharedPtr();
113 446 : auto mfem_bc = std::dynamic_pointer_cast<MFEMEssentialBC>(object_ptr);
114 446 : if (getProblemData().eqn_system)
115 : {
116 446 : 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 446 : }
124 : else
125 : {
126 0 : mooseError("Unsupported bc of type '", bc_name, "' and name '", name, "' detected.");
127 : }
128 476 : }
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 137 : MFEMProblem::addFunctorMaterial(const std::string & material_name,
139 : const std::string & name,
140 : InputParameters & parameters)
141 : {
142 137 : FEProblemBase::addUserObject(material_name, name, parameters);
143 121 : getUserObject<MFEMFunctorMaterial>(name);
144 121 : }
145 :
146 : void
147 661 : MFEMProblem::addFESpace(const std::string & user_object_name,
148 : const std::string & name,
149 : InputParameters & parameters)
150 : {
151 661 : FEProblemBase::addUserObject(user_object_name, name, parameters);
152 661 : MFEMFESpace & mfem_fespace(getUserObject<MFEMFESpace>(name));
153 :
154 : // Register fespace and associated fe collection.
155 661 : getProblemData().fecs.Register(name, mfem_fespace.getFEC());
156 661 : getProblemData().fespaces.Register(name, mfem_fespace.getFESpace());
157 661 : }
158 :
159 : void
160 386 : MFEMProblem::addVariable(const std::string & var_type,
161 : const std::string & var_name,
162 : InputParameters & parameters)
163 : {
164 386 : 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 386 : if (isTransient())
169 120 : addGridFunction(var_type, Moose::MFEM::GetTimeDerivativeName(var_name), parameters);
170 386 : }
171 :
172 : void
173 863 : MFEMProblem::addGridFunction(const std::string & var_type,
174 : const std::string & var_name,
175 : InputParameters & parameters)
176 : {
177 863 : if (var_type == "MFEMVariable")
178 : {
179 : // Add MFEM variable directly.
180 821 : 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 863 : MFEMVariable & mfem_variable = getUserObject<MFEMVariable>(var_name);
194 863 : getProblemData().gridfunctions.Register(var_name, mfem_variable.getGridFunction());
195 863 : if (mfem_variable.getFESpace().isScalar())
196 528 : getCoefficients().declareScalar<mfem::GridFunctionCoefficient>(
197 528 : var_name, mfem_variable.getGridFunction().get());
198 : else
199 335 : getCoefficients().declareVector<mfem::VectorGridFunctionCoefficient>(
200 335 : var_name, mfem_variable.getGridFunction().get());
201 863 : }
202 :
203 : void
204 357 : 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 357 : addGridFunction(var_type, var_name, parameters);
211 357 : }
212 :
213 : void
214 190 : MFEMProblem::addAuxKernel(const std::string & kernel_name,
215 : const std::string & name,
216 : InputParameters & parameters)
217 : {
218 190 : FEProblemBase::addUserObject(kernel_name, name, parameters);
219 190 : }
220 :
221 : void
222 548 : MFEMProblem::addKernel(const std::string & kernel_name,
223 : const std::string & name,
224 : InputParameters & parameters)
225 : {
226 548 : FEProblemBase::addUserObject(kernel_name, name, parameters);
227 548 : const UserObject * kernel_uo = &(getUserObjectBase(name));
228 :
229 548 : if (dynamic_cast<const MFEMKernel *>(kernel_uo) != nullptr)
230 : {
231 548 : auto object_ptr = getUserObject<MFEMKernel>(name).getSharedPtr();
232 548 : auto kernel = std::dynamic_pointer_cast<MFEMKernel>(object_ptr);
233 548 : if (getProblemData().eqn_system)
234 : {
235 548 : 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 548 : }
243 : else
244 : {
245 0 : mooseError("Unsupported kernel of type '", kernel_name, "' and name '", name, "' detected.");
246 : }
247 548 : }
248 :
249 : libMesh::Point
250 2281620 : pointFromMFEMVector(const mfem::Vector & vec)
251 : {
252 : return libMesh::Point(
253 2281620 : vec.Elem(0), vec.Size() > 1 ? vec.Elem(1) : 0., vec.Size() > 2 ? vec.Elem(2) : 0.);
254 : }
255 :
256 : int
257 141 : vectorFunctionDim(const std::string & type, const InputParameters & parameters)
258 : {
259 141 : if (type == "LevelSetOlssonVortex")
260 : {
261 0 : return 2;
262 : }
263 141 : else if (type == "ParsedVectorFunction")
264 : {
265 463 : if (parameters.isParamSetByUser("expression_z") || parameters.isParamSetByUser("value_z"))
266 : {
267 121 : return 3;
268 : }
269 64 : else if (parameters.isParamSetByUser("expression_y") || parameters.isParamSetByUser("value_y"))
270 : {
271 18 : return 2;
272 : }
273 : else
274 : {
275 2 : 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 241 : MFEMProblem::addFunction(const std::string & type,
326 : const std::string & name,
327 : InputParameters & parameters)
328 : {
329 241 : ExternalProblem::addFunction(type, name, parameters);
330 241 : 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 241 : if (std::find(SCALAR_FUNCS.begin(), SCALAR_FUNCS.end(), type) != SCALAR_FUNCS.end())
334 : {
335 98 : getCoefficients().declareScalar<mfem::FunctionCoefficient>(
336 : name,
337 98 : [&func](const mfem::Vector & p, mfem::real_t t) -> mfem::real_t
338 711278 : { return func.value(t, pointFromMFEMVector(p)); });
339 : }
340 143 : else if (std::find(VECTOR_FUNCS.begin(), VECTOR_FUNCS.end(), type) != VECTOR_FUNCS.end())
341 : {
342 141 : int dim = vectorFunctionDim(type, parameters);
343 141 : getCoefficients().declareVector<mfem::VectorFunctionCoefficient>(
344 : name,
345 : dim,
346 141 : [&func, dim](const mfem::Vector & p, mfem::real_t t, mfem::Vector & u)
347 : {
348 1570342 : libMesh::RealVectorValue vector_value = func.vectorValue(t, pointFromMFEMVector(p));
349 5346466 : for (int i = 0; i < dim; i++)
350 : {
351 3776124 : u[i] = vector_value(i);
352 : }
353 1570342 : });
354 : }
355 : else
356 : {
357 2 : mooseWarning("Could not identify whether function ",
358 : type,
359 : " is scalar or vector; no MFEM coefficient object created.");
360 : }
361 239 : }
362 :
363 : void
364 78 : MFEMProblem::addPostprocessor(const std::string & type,
365 : const std::string & name,
366 : InputParameters & parameters)
367 : {
368 78 : ExternalProblem::addPostprocessor(type, name, parameters);
369 78 : const PostprocessorValue & val = getPostprocessorValueByName(name);
370 78 : getCoefficients().declareScalar<mfem::FunctionCoefficient>(
371 80 : name, [&val](const mfem::Vector &) -> mfem::real_t { return val; });
372 78 : }
373 :
374 : InputParameters
375 42 : MFEMProblem::addMFEMFESpaceFromMOOSEVariable(InputParameters & parameters)
376 : {
377 :
378 84 : InputParameters fespace_params = _factory.getValidParams("MFEMGenericFESpace");
379 84 : InputParameters mfem_variable_params = _factory.getValidParams("MFEMVariable");
380 :
381 : auto moose_fe_type =
382 84 : FEType(Utility::string_to_enum<Order>(parameters.get<MooseEnum>("order")),
383 126 : Utility::string_to_enum<FEFamily>(parameters.get<MooseEnum>("family")));
384 :
385 42 : std::string mfem_family;
386 42 : int mfem_vdim = 1;
387 :
388 42 : switch (moose_fe_type.family)
389 : {
390 14 : case FEFamily::LAGRANGE:
391 14 : mfem_family = "H1";
392 14 : mfem_vdim = 1;
393 14 : break;
394 14 : case FEFamily::LAGRANGE_VEC:
395 14 : mfem_family = "H1";
396 14 : mfem_vdim = 3;
397 14 : break;
398 7 : case FEFamily::MONOMIAL:
399 7 : mfem_family = "L2";
400 7 : mfem_vdim = 1;
401 7 : break;
402 7 : case FEFamily::MONOMIAL_VEC:
403 7 : mfem_family = "L2";
404 7 : mfem_vdim = 3;
405 7 : break;
406 0 : default:
407 0 : mooseError("Unable to set MFEM FESpace for MOOSE variable");
408 : break;
409 : }
410 :
411 : // Create fespace name. If this already exists, we will reuse this for
412 : // the mfem variable ("gridfunction").
413 84 : const std::string fespace_name = mfem_family + "_" +
414 168 : std::to_string(mesh().getMFEMParMesh().Dimension()) + "D_P" +
415 168 : std::to_string(moose_fe_type.order.get_order());
416 :
417 : // Set all fespace parameters.
418 42 : fespace_params.set<std::string>("fec_name") = fespace_name;
419 84 : fespace_params.set<int>("vdim") = mfem_vdim;
420 :
421 42 : if (!hasUserObject(fespace_name)) // Create the fespace (implicit).
422 : {
423 28 : addFESpace("MFEMGenericFESpace", fespace_name, fespace_params);
424 : }
425 :
426 126 : mfem_variable_params.set<UserObjectName>("fespace") = fespace_name;
427 :
428 84 : return mfem_variable_params;
429 42 : }
430 :
431 : void
432 877 : MFEMProblem::displaceMesh()
433 : {
434 : // Displace mesh
435 877 : if (mesh().shouldDisplace())
436 : {
437 13 : mesh().displace(static_cast<mfem::GridFunction const &>(*getMeshDisplacementGridFunction()));
438 : // TODO: update FESpaces GridFunctions etc for transient solves
439 : }
440 877 : }
441 :
442 : std::optional<std::reference_wrapper<mfem::ParGridFunction const>>
443 362 : MFEMProblem::getMeshDisplacementGridFunction()
444 : {
445 : // If C++23 transform were available this would be easier
446 362 : auto const displacement_variable = mesh().getMeshDisplacementVariable();
447 362 : if (displacement_variable)
448 : {
449 26 : return *_problem_data.gridfunctions.Get(displacement_variable.value());
450 : }
451 : else
452 : {
453 336 : return std::nullopt;
454 : }
455 : }
456 :
457 : std::vector<VariableName>
458 0 : MFEMProblem::getAuxVariableNames()
459 : {
460 0 : return systemBaseAuxiliary().getVariableNames();
461 : }
462 :
463 : MFEMMesh &
464 11456 : MFEMProblem::mesh()
465 : {
466 : mooseAssert(ExternalProblem::mesh().type() == "MFEMMesh",
467 : "Please choose the MFEMMesh mesh type for an MFEMProblem\n");
468 11456 : return static_cast<MFEMMesh &>(_mesh);
469 : }
470 :
471 : const MFEMMesh &
472 1417 : MFEMProblem::mesh() const
473 : {
474 1417 : return const_cast<MFEMProblem *>(this)->mesh();
475 : }
476 :
477 : void
478 82 : MFEMProblem::addSubMesh(const std::string & var_type,
479 : const std::string & var_name,
480 : InputParameters & parameters)
481 : {
482 : // Add MFEM SubMesh.
483 82 : FEProblemBase::addUserObject(var_type, var_name, parameters);
484 : // Register submesh.
485 82 : MFEMSubMesh & mfem_submesh = getUserObject<MFEMSubMesh>(var_name);
486 82 : getProblemData().submeshes.Register(var_name, mfem_submesh.getSubMesh());
487 82 : }
488 :
489 : void
490 140 : MFEMProblem::addTransfer(const std::string & transfer_name,
491 : const std::string & name,
492 : InputParameters & parameters)
493 : {
494 140 : if (parameters.getBase() == "MFEMSubMeshTransfer")
495 96 : FEProblemBase::addUserObject(transfer_name, name, parameters);
496 : else
497 44 : FEProblemBase::addTransfer(transfer_name, name, parameters);
498 140 : }
499 :
500 : std::shared_ptr<mfem::ParGridFunction>
501 85 : MFEMProblem::getGridFunction(const std::string & name)
502 : {
503 85 : return getUserObject<MFEMVariable>(name).getGridFunction();
504 : }
505 :
506 : void
507 85 : MFEMProblem::addInitialCondition(const std::string & ic_name,
508 : const std::string & name,
509 : InputParameters & parameters)
510 : {
511 85 : FEProblemBase::addUserObject(ic_name, name, parameters);
512 85 : getUserObject<MFEMInitialCondition>(name); // error check
513 85 : }
514 :
515 : std::string
516 349 : MFEMProblem::solverTypeString(const unsigned int libmesh_dbg_var(solver_sys_num))
517 : {
518 : mooseAssert(solver_sys_num == 0, "No support for multi-system with MFEM right now");
519 349 : return MooseUtils::prettyCppType(getProblemData().jacobian_solver.get());
520 : }
521 :
522 : #endif
|
