diff --git a/.gitignore b/.gitignore index 07bd0d1..d76bac3 100644 --- a/.gitignore +++ b/.gitignore @@ -7,3 +7,9 @@ *.snm *.toc *.vrb +*.fdb_latexmk +*.fls +*.pdf +*.spl +*.bbl +*.blg diff --git a/papers/2024_softwarex/Makefile b/papers/2024_softwarex/Makefile new file mode 100644 index 0000000..99d1a9e --- /dev/null +++ b/papers/2024_softwarex/Makefile @@ -0,0 +1,14 @@ +NAME=paper + +all: $(NAME).pdf + +# We have a custom .latexmkrc file in this directory which tells latexmk how to make the glossary. +# https://tex.stackexchange.com/questions/1226/how-to-make-latexmk-use-makeglossaries +$(NAME).pdf: $(NAME).tex + latexmk -g -pdf -logfilewarnings $< + +# Removes all temporary files including output files (dvi/ps/pdf). +# Note: "latexmk -c" doesn't remove output files. +clean: + latexmk -C + rm -f *~ .*~ diff --git a/papers/2024_softwarex/figures/netgen_hole_demo.png b/papers/2024_softwarex/figures/netgen_hole_demo.png new file mode 100644 index 0000000..d6b8f89 Binary files /dev/null and b/papers/2024_softwarex/figures/netgen_hole_demo.png differ diff --git a/papers/2024_softwarex/figures/vector_fe_ex9_streamlines.png b/papers/2024_softwarex/figures/vector_fe_ex9_streamlines.png new file mode 100644 index 0000000..7cd3abe Binary files /dev/null and b/papers/2024_softwarex/figures/vector_fe_ex9_streamlines.png differ diff --git a/papers/2024_softwarex/paper.tex b/papers/2024_softwarex/paper.tex new file mode 100644 index 0000000..c84b650 --- /dev/null +++ b/papers/2024_softwarex/paper.tex @@ -0,0 +1,388 @@ +%% +%% Copyright 2026 Elsevier Ltd +%% +%% This file is part of the 'Elsarticle Bundle'. +%% --------------------------------------------- +%% +%% It may be distributed under the conditions of the LaTeX Project Public +%% License, either version 1.2 of this license or (at your option) any +%% later version. The latest version of this license is in +%% http://www.latex-project.org/lppl.txt +%% and version 1.2 or later is part of all distributions of LaTeX +%% version 1999/12/01 or later. +%% +%% The list of all files belonging to the 'Elsarticle Bundle' is +%% given in the file `manifest.txt'. +%% + +%% Template article for Elsevier's document class `elsarticle' +%% with numbered style bibliographic references +%% SP 2008/03/01 + +\documentclass[preprint,12pt, a4paper]{elsarticle} + +%% Use the option review to obtain double line spacing +%% \documentclass[authoryear,preprint,review,12pt]{elsarticle} + +%% For including figures, graphicx.sty has been loaded in +%% elsarticle.cls. If you prefer to use the old commands +%% please give \usepackage{epsfig} + +%% The amssymb package provides various useful mathematical symbols +\usepackage{amssymb} +\usepackage{hyperref} +\setlength{\parindent}{0pt} +%% The amsthm package provides extended theorem environments +%% \usepackage{amsthm} + +\journal{SoftwareX} + +\begin{document} +\renewcommand{\labelenumii}{\arabic{enumi}.\arabic{enumii}} + +\begin{frontmatter} + +\title{libMesh: A C++ Library for Finite Element and Finite Volume Based Multiphysics} + +\author[inl]{Roy H. Stogner} +\author[akselos]{John W. Peterson} +\author[inl]{Alexander D. Lindsay} +\author[akselos]{David J. Knezevic} +\author[hartree]{Nuno Nobre} +\author[inl]{Logan Harbour} +\author[inl]{Guillaume L. Giudicelli} +\author[inl]{Derek R. Gaston} +\author[inl]{Cody J. Permann} +\author[ub]{Paul T. Bauman} +\author[inlcmm]{Daniel Schwen} +\author[inl]{Patrick Behne} +\author[ut]{Benjamin Kirk} + +\address[inl]{Computational Frameworks, Idaho National Laboratory, Idaho Falls, ID 83415, USA} +\address[akselos]{Akselos} +\address[hartree]{Hartree Centre, Science and Technology Facilities Council (STFC), UK} +\address[ub]{University at Buffalo, Department of Mechanical and Aerospace Engineering, Buffalo, NY, USA} +\address[inlcmm]{Computational Mechanics and Materials, Idaho National Laboratory, Idaho Falls, ID 83415, USA} +\address[ut]{Department of Aerospace Engineering and Engineering Mechanics (ASE), University of Texas at Austin, Austin, TX, USA} + +\begin{abstract} +libMesh is an open-source C++ library for parallel, spatially adaptive finite element and finite +volume simulation, previously described in a 2006 publication. In the two decades since, the +library has grown substantially: it now supports hybridizable discontinuous Galerkin and other +statically condensed discretizations, isogeometric analysis with rational finite element bases and +arbitrary trimmed-mesh extraction operators, native polyhedral elements, and a standalone parallel +communication submodule (TIMPI) implementing modern nonblocking-consensus algorithms. Automatic +differentiation, certified reduced-basis model order reduction, and adjoint-driven $hp$-adaptivity +round out a substantially expanded numerical toolkit. This article surveys these capabilities and +the multiphysics applications they enable. +\end{abstract} + +\begin{keyword} +finite elements \sep finite volume \sep parallel computing \sep isogeometric analysis \sep multiphysics \sep adaptive mesh refinement +\end{keyword} + +\end{frontmatter} + +\section*{Required Metadata} + +\section*{Current code version} + +\begin{table}[!h] +\begin{tabular}{|l|p{6.5cm}|p{6.5cm}|} +\hline +\textbf{Nr.} & \textbf{Code metadata description} & \textbf{Metadata} \\ +\hline +C1 & Current code version & v1.9.0 development series (confirm exact tag at submission) \\ +\hline +C2 & Permanent GitHub link to code/repository used for this code version & \url{https://github.com/libMesh/libmesh} \\ +\hline +C3 & Legal Code License & GNU Lesser General Public License v2.1 (LGPL-2.1) \\ +\hline +C4 & Code versioning system used & git (hosted on GitHub) \\ +\hline +C5 & Software code languages, tools, and services used & C++ (C++17); MPI; PETSc/SLEPc (optional); TIMPI (bundled, required); MetaPhysicL (bundled, optional); Eigen; TetGen, Netgen, Triangle (bundled); ExodusII/SEACAS; NetCDF; VTK \\ +\hline +C6 & Compilation requirements, operating environments \& dependencies & Autoconf/Automake build (\texttt{configure} \& \texttt{make}); C++17 compiler; MPI implementation (no supported build configuration omits it in practice); TIMPI submodule (required); PETSc strongly recommended for solvers and static condensation \\ +\hline +C7 & If available Link to developer documentation/manual & \url{https://libmesh.github.io/} \\ +\hline +C8 & Support email for questions & \url{https://github.com/libmesh/libmesh/discussions} \\ +\hline +\end{tabular} +\caption{Code metadata (mandatory)} +\label{tab:metadata} +\end{table} + +\begin{enumerate} + +\item Motivation and significance + +libMesh is an open-source C++ library that provides a general, physics-independent +infrastructure for parallel, spatially adaptive finite element -- and, increasingly, finite +volume -- simulation. An earlier publication described the library's architecture as of 2006: a +mesh and degree-of-freedom data model, $h$-adaptive refinement and coarsening with hanging-node +constraints, and a domain-decomposition-based parallel model built atop third-party solver +libraries such as PETSc~\cite{kirk2006libmesh}. That description, while still broadly applicable +to the library's core design, predates nearly all of the capabilities surveyed in this article. + +In the twenty years since, libMesh has remained under active, continuous development -- the +codebase now spans tens of thousands of commits from a broad developer community -- and has grown +well beyond its original scope. Applications built on libMesh increasingly require discretizations +and geometric representations unavailable in 2006: hybridizable and other statically condensed +formulations well suited to advection-dominated transport problems, isogeometric meshes derived +from trimmed or piecewise-rational-spline CAD geometry, and general polyhedral cells that are a +much better fit for finite-volume-style discretizations than the fixed simplicial and +tensor-product element types available previously. The library's parallel communication layer has +also been substantially reworked, and new numerical infrastructure -- automatic differentiation, +certified reduced-order modeling, and goal-oriented adaptivity -- has been added. + +These capabilities are motivated by, and used in, two largely independent application +communities. The MOOSE multiphysics framework~\cite{permann2020moose}, developed at Idaho National +Laboratory (INL), relies on libMesh for both finite element and, through the discontinuous/monomial +discretizations and generalized ghosting infrastructure described in Section 2, finite-volume-style +capabilities across INL's fission and fusion reactor multiphysics work; a comparable MOOSE-based +fusion-device multiphysics effort is under way in the UK, where the Hartree Centre (STFC) +collaborates closely with the UK Atomic Energy Authority. Independently, Akselos builds commercial +reduced-order-modeling and digital-twin engineering-analysis products directly on libMesh's +certified reduced-basis infrastructure, a use case with no dependency on MOOSE. This article +surveys the library-level capabilities underlying both communities. + +\item Software description + +\begin{enumerate} + +\item Software architecture + +libMesh is organized around a small set of core abstractions -- \texttt{Mesh}, \texttt{Elem}, +\texttt{DofMap}, \texttt{System}, and finite element (\texttt{FEBase}) classes -- that are +deliberately physics-independent: an application supplies the equations to be solved, and libMesh +supplies mesh management, degree-of-freedom numbering, adaptivity, and parallel communication +around that specification. The capabilities described below extend this same architecture rather +than replacing it: new finite element families, element types, and mesh constructs plug into the +existing \texttt{FEBase}/\texttt{Elem} hierarchy, new solver and preconditioning machinery plugs +into the existing \texttt{System}/\texttt{DofMap} interfaces, and new parallel primitives are +exposed through the bundled TIMPI submodule rather than through a separate communication layer. + +\item Software functionalities + +\textbf{Discretizations for finite element and finite volume multiphysics.} libMesh implements +hybridization and static condensation through the \texttt{StaticCondensation}, +\texttt{StaticCondensationDofMap}, and \texttt{StaticCondensationPreconditioner} classes, which +assemble a full elemental system and eliminate element-local degrees of freedom via a local Schur +complement (using Eigen for the local dense factorizations) before the global solve. Paired with +the \texttt{SIDE\_HIERARCHIC} finite element family, which represents variables on element facets +rather than volumetrically, this supports hybridizable discontinuous Galerkin (HDG) discretizations +for general transport problems, including the advection-dominated regimes in which HDG is used +downstream in MOOSE-based applications. General polyhedral cells are supported natively through the +\texttt{Polyhedron} and \texttt{C0Polyhedron} classes, whose topology (number of faces, sides, and +nodes) is determined at run time rather than fixed at compile time as with libMesh's traditional +hexahedral, tetrahedral, and prismatic elements; polyhedral cells -- including coarser cells built +by merging groups of finer elements -- are a substantially better match for finite-volume-style +stencils than fixed simplicial elements. Finite-volume-style discretizations more broadly are +enabled by libMesh's discontinuous/monomial finite element spaces together with a generalized +ghosting-functor framework (\texttt{GhostingFunctor} and its specializations for point-neighbor, +overlap, non-manifold, and sibling coupling), which lets an application declare exactly which +off-element data a cell-centered discretization needs without a dedicated finite-volume module +inside libMesh itself. + +\textbf{Isogeometric analysis.} A \texttt{RATIONAL\_BERNSTEIN} finite element family reweights the +underlying Bernstein/B\'ezier shape functions with per-node NURBS weights, giving the classic +rational-basis interpolation used in isogeometric analysis (IGA)~\cite{hughes2005iga}; a +corresponding element mapping type lets the geometric map itself, not only the solution field, use +the same rational representation, so that meshes derived from piecewise-rational-spline CAD +geometry can be represented exactly. Arbitrary B\'ezier extraction operators -- needed when a +traditional, structured IGA mesh is impractical to generate automatically, as with trimmed +geometry -- are supported by reusing libMesh's existing generic linear-constraint infrastructure: +extraction coefficients read from ExodusII files are stored as constraint rows relating spline-basis +degrees of freedom to element-local nodes, requiring no separate, IGA-specific parallel code +path~\cite{borden2011extraction}. + +\textbf{Meshing and I/O.} A \texttt{MeshTetInterface} abstract base class unifies libMesh's +interfaces to the bundled TetGen and Netgen tetrahedralization libraries, giving applications a +common API for generating volume meshes from boundary surfaces regardless of which third-party +mesher is used. Mesh I/O has grown substantially beyond the formats available in 2006, with +extended ExodusII/Nemesis, VTK, Abaqus, and NetCDF support, and the element library itself has grown +to include additional high-order and shell element types together with generic element-quality +metrics (aspect ratio, minimum/maximum angle, and scaled Jacobian). + +\textbf{Parallel infrastructure.} Parallel communication is now implemented in TIMPI, a submodule +that began as libMesh-internal code and was split into its own independently versioned and tested +repository in 2019. TIMPI exposes multiple interchangeable synchronization strategies for +parallel push/pull-style data exchange, selectable through a \texttt{SyncType} enumeration +(\texttt{NBX}, \texttt{ALLTOALL\_COUNTS}, \texttt{SENDRECEIVE}); its \texttt{NBX} strategy +implements the nonblocking-consensus dynamic sparse data exchange algorithm of +Hoefler et al.~\cite{hoefler2010nbx}, which avoids the all-to-all communication and a priori +knowledge of message counts that earlier, simpler synchronization schemes required. TIMPI is a +required dependency of the library, not an optional one. + +\textbf{Numerics and solvers.} libMesh's core generic numeric types -- \texttt{TypeTensor}, +\texttt{TypeVector}, \texttt{VectorValue}, \texttt{TensorValue}, and \texttt{DenseMatrix}/ +\texttt{DenseVector} -- are templated to work transparently with the dual-number type provided by +the bundled MetaPhysicL submodule. This is what allows downstream codes, notably +MOOSE~\cite{lindsay2021metaphysicl}, to build exact-derivative Jacobian assembly on top of libMesh +without any automatic-differentiation-specific code path inside libMesh itself. A certified +reduced-basis model-order-reduction capability (offline/online reduced-basis construction, the +empirical interpolation method for nonaffine operators, and the successive constraint method for +stability-factor lower bounds) provides the numerical core of Akselos's commercial +reduced-order-modeling products. Adjoint-based, goal-oriented error estimation and $hp$-adaptive +refinement strategies extend the simpler flux-jump and patch-recovery indicators available in 2006 +to quantity-of-interest-driven adaptivity. + +\item Sample code snippets analysis + +The following excerpt, adapted from libMesh's HDG example driver, illustrates how a hybridized +discretization is declared: an element-interior vector unknown, its scalar components, a facet +Lagrange multiplier represented with \texttt{SIDE\_HIERARCHIC}, and a global scalar constraint are +added as ordinary libMesh \texttt{System} variables, after which the \texttt{StaticCondensation} +matrix class is attached to eliminate the element-interior unknowns before the global solve: + +\begin{verbatim} +system.add_variable("qu", FIRST, L2_LAGRANGE_VEC); +system.add_variable("vel_x", FIRST, L2_LAGRANGE); +system.add_variable("lm_u", FIRST, SIDE_HIERARCHIC); +system.add_variable("lm_v", FIRST, SIDE_HIERARCHIC); +system.add_variable("pressure", FIRST, L2_LAGRANGE); +system.add_variable("global_lm", FIRST, SCALAR); +\end{verbatim} + +No new solver, assembly, or parallel-communication code is required for this to run in parallel: +the same \texttt{DofMap}, ghosting, and TIMPI-based communication infrastructure used by +conventional continuous finite element systems apply directly to the hybridized system. + +\end{enumerate} + +\item Illustrative examples + +libMesh's \texttt{examples/} tree includes runnable drivers that exercise the capabilities above +end to end, and its test suite exercises capabilities for which a dedicated example does not yet +exist. Table~\ref{tab:examples} summarizes, for each capability surveyed in Section 2, a +representative runnable driver or test. + +\texttt{examples/vector\_fe/vector\_fe\_ex9} solves a Stokes-like vector transport problem using +the HDG formulation described in Section 2.2: an element-interior velocity/pressure system is +coupled to facet Lagrange multipliers through \texttt{SIDE\_HIERARCHIC} variables, with +\texttt{StaticCondensation} eliminating the interior unknowns before the global (facet-only) system +is solved. Figure~\ref{fig:hdg} shows velocity streamlines for the resulting lid-driven-cavity flow, +colored by speed, exhibiting the expected primary recirculation vortex. + +\begin{figure}[!h] +\centering +\includegraphics[width=0.55\textwidth]{figures/vector_fe_ex9_streamlines.png} +\caption{Velocity streamlines, colored by speed, for the HDG cavity-flow problem solved by +\texttt{examples/vector\_fe/vector\_fe\_ex9}, after static condensation of the element-interior +velocity/pressure unknowns.} +\label{fig:hdg} +\end{figure} + +The unified tetrahedralization interface is illustrated with a minimal driver, in the same spirit +as \texttt{examples/miscellaneous/miscellaneous\_ex6}, that attaches a cube-shaped hole to a cubic +domain and hands both to the \texttt{MeshTetInterface}-derived \texttt{NetGenMeshInterface} to +produce a constrained-quality volume mesh; Figure~\ref{fig:mesh} shows a clipped view of the +resulting tetrahedralization, with the interior hole and the fill between it and the outer boundary +both visible. The same abstraction accepts a \texttt{TetGenMeshInterface} in place of the Netgen +backend without any change to the calling code, which is the uniformity across bundled mesh +generators that \texttt{MeshTetInterface} is intended to provide. + +\begin{figure}[!h] +\centering +\includegraphics[width=0.55\textwidth]{figures/netgen_hole_demo.png} +\caption{Clipped view of a tetrahedralized cubic domain with a cube-shaped interior hole, generated +via \texttt{MeshTetInterface}'s Netgen backend.} +\label{fig:mesh} +\end{figure} + +\begin{table}[!h] +\begin{tabular}{|p{4.3cm}|p{8.7cm}|} +\hline +\textbf{Capability} & \textbf{Runnable example or test} \\ +\hline +Hybridization / static condensation (HDG) & \texttt{examples/vector\_fe/vector\_fe\_ex9} \\ +\hline +Native polyhedral elements & \texttt{tests/geom/volume\_test.C} +(\texttt{testC0PolyhedronCube}, \texttt{testC0PolyhedronHexagonalPrism}) \\ +\hline +Isogeometric rational bases \& B\'ezier extraction & \texttt{tests/fe/fe\_rational\_test.C}, +\texttt{tests/fe/fe\_rational\_map.C}, and the two-element IGA regression mesh +\texttt{tests/meshes/two\_element\_iga\_in.e}; a dedicated example is planned \\ +\hline +Unified tetrahedralization (TetGen/Netgen) & \texttt{examples/miscellaneous/miscellaneous\_ex6} \\ +\hline +TIMPI synchronization strategies (\texttt{NBX}, \texttt{ALLTOALL\_COUNTS}, +\texttt{SENDRECEIVE}) & \texttt{contrib/timpi/test/parallel\_sync\_unit.C} \\ +\hline +Certified reduced-basis model order reduction & \texttt{examples/reduced\_basis/reduced\_basis\_ex1} \\ +\hline +Automatic differentiation (dual-number-compatible numerics) & No dedicated libMesh-side example; +exercised downstream by AD-enabled Jacobian assembly in MOOSE \\ +\hline +\end{tabular} +\caption{Runnable examples and tests exercising the capabilities surveyed in Section 2.} +\label{tab:examples} +\end{table} + +\item Impact + +libMesh's expanded capabilities enable classes of problems that were not reachable with the 2006 +codebase: hybridized discretizations for advection-dominated transport, trimmed and +piecewise-rational-spline isogeometric geometries, finite-volume-style multiphysics on general +polyhedral meshes, and exact-derivative Jacobian assembly via dual-number-compatible numerics. +These capabilities underlie the two application communities described in Section 1, at comparable +weight. Decoupling TIMPI into its own repository has also had a maintenance-level impact: TIMPI is +now independently versioned and tested, and is reusable by other parallel scientific codes without +a dependency on the rest of libMesh. + +\item Conclusions + +libMesh's original identity as a parallel, adaptive, physics-independent finite element library is +intact, but the library today spans considerably more ground than its 2006 description: finite +element and finite-volume multiphysics on general polyhedral meshes, isogeometric analysis on both +structured and trimmed geometry, hybridized discretizations for advection-dominated transport, and +supporting parallel, automatic-differentiation, and model-order-reduction infrastructure. These +capabilities are the product of two decades of contributions from a broad developer community; +several of them -- native polyhedral elements and continued refinement of the B\'ezier extraction +machinery -- were added or substantially revised within the past year, evidence that libMesh remains +an actively developed project responding to the needs of its national-laboratory, academic, and +commercial user base. + +\end{enumerate} + +\section*{Acknowledgements} + +\begin{thebibliography}{00} + +\bibitem{kirk2006libmesh} B.~S. Kirk, J.~W. Peterson, R.~H. Stogner, G.~F. Carey, libMesh: a C++ +library for parallel adaptive mesh refinement/coarsening simulations, Engineering with Computers +22 (2006) 237--254. + +\bibitem{permann2020moose} C.~J. Permann, D.~R. Gaston, D.~Andr\v{s}, R.~W. Carlsen, F.~Kong, +A.~D. Lindsay, J.~M. Miller, J.~W. Peterson, A.~E. Slaughter, R.~H. Stogner, R.~C. Martineau, MOOSE: +Enabling massively parallel multiphysics simulation, SoftwareX 11 (2020) 100430. + +\bibitem{hoefler2010nbx} T.~Hoefler, C.~Siebert, A.~Lumsdaine, Scalable communication protocols for +dynamic sparse data exchange, in: Proceedings of the 15th ACM SIGPLAN Symposium on Principles and +Practice of Parallel Programming (PPoPP '10), ACM, 2010, pp. 159--168. + +\bibitem{hughes2005iga} T.~J.~R. Hughes, J.~A. Cottrell, Y.~Bazilevs, Isogeometric analysis: CAD, +finite elements, NURBS, exact geometry and mesh refinement, Computer Methods in Applied Mechanics +and Engineering 194 (39--41) (2005) 4135--4195. + +\bibitem{borden2011extraction} M.~J. Borden, M.~A. Scott, J.~A. Evans, T.~J.~R. Hughes, Isogeometric +finite element data structures based on B\'ezier extraction of NURBS, International Journal for +Numerical Methods in Engineering 87 (1--5) (2011) 15--47. + +\bibitem{lindsay2021metaphysicl} A.~D. Lindsay, R.~H. Stogner, D.~R. Gaston, D.~Schwen, +C.~Matthews, W.~Jiang, L.~K. Aagesen, R.~Carlsen, F.~Kong, A.~E. Slaughter, et al., Automatic +differentiation in MetaPhysicL and its applications in MOOSE, Nuclear Technology 207 (7) (2021) +905--922. + +\end{thebibliography} + +\end{document} +\endinput +%% +%% End of file `SoftwareX_article_template.tex'. + +%%% Local Variables: +%%% mode: latex +%%% TeX-master: t +%%% End: