Merge pull request #10 from DamynChipman/v1-joss-paper

v1.0.2 - JOSS Paper

.github/workflows/draft-pdf.yml

@@ -0,0 +1,23 @@
+on: [push]
+
+jobs:
+  paper:
+    runs-on: ubuntu-latest
+    name: Paper Draft
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v4
+      - name: Build draft PDF
+        uses: openjournals/openjournals-draft-action@master
+        with:
+          journal: joss
+          # This should be the path to the paper within your repo.
+          paper-path: paper.md
+      - name: Upload
+        uses: actions/upload-artifact@v1
+        with:
+          name: paper
+          # This is the output path where Pandoc will write the compiled
+          # PDF. Note, this should be the same directory as the input
+          # paper.md
+          path: paper.pdf

examples/elliptic-single/main.cpp

@@ -276,4 +276,4 @@ int main(int argc, char** argv) {
 
     // All clean up is done in destructors
     return EXIT_SUCCESS;
-}
+}

paper.bib

@@ -0,0 +1,151 @@
+@article{chipman2024fast,
+	title={A Fast Direct Solver for Elliptic PDEs on a Hierarchy of Adaptively Refined Quadtrees},
+	author={Chipman, D and Calhoun, D and Burstedde, C},
+	journal={arXiv preprint arXiv:2402.14936},
+	year={2024},
+	doi={10.48550/arXiv.2402.14936}
+}
+
+@article{chipman2024parallel,
+	author = {Chipman, D and Calhoun, D and Burstedde, C},
+	journal = {[Forthcoming]},
+	title = {{MPI Implementation of the Quadtree-Adaptive Hierarchical Poincar{\'e}-Steklov Method}},
+	year = {2024}
+}
+
+@article{gillman2014direct,
+    author = {Gillman, A. and Martinsson, P. G.},
+    title = {{A Direct Solver with $O(N)$ Complexity for Variable Coefficient Elliptic PDEs Discretized via a High-Order Composite Spectral Collocation Method}},
+    journal = {SIAM Journal on Scientific Computing},
+    volume = {36},
+    number = {4},
+    pages = {A2023-A2046},
+    year = {2014},
+    doi = {10.1137/130918988},
+    url = {https://doi.org/10.1137/130918988},
+    eprint = {https://doi.org/10.1137/130918988},
+}
+
+@book{martinsson2019fast,
+	author = {Martinsson, PG},
+	publisher = {SIAM},
+	title = {{Fast direct solvers for elliptic PDEs}},
+	year = {2019},
+	doi={10.1137/1.9781611976045}
+}
+
+@article{swarztrauber1999fishpack,
+	author = {Swarztrauber, P and Sweet, R and Adams, J},
+	journal = {UCAR Publication, July},
+	title = {FISHPACK: Efficient FORTRAN subprograms for the solution of elliptic partial differential equations},
+	year = {1999},
+	doi={10.1145/800207.806417}
+}
+
+@article{adams2016fishpack90,
+	author = {Adams, J C and Swarztrauber, P and Sweet, R},
+	journal = {Astrophysics Source Code Library},
+	pages = {ascl--1609},
+	title = {FISHPACK90: Efficient fortran subprograms for the solution of separable elliptic partial differential equations},
+	year = {2016},
+	url={https://ui.adsabs.harvard.edu/abs/2016ascl.soft09005A/abstract}
+}
+
+@manual{mpi41,
+    author = "{Message Passing Interface Forum}",
+    title  = "{MPI}: A Message-Passing Interface Standard Version 4.1",
+    url    = "https://www.mpi-forum.org/docs/mpi-4.1/mpi41-report.pdf",
+    year   = 2023,
+    month  = nov
+}
+
+@article{burstedde2011p4est,
+	author = {Burstedde, C and Wilcox, L and Ghattas, O},
+	journal = {SIAM Journal on Scientific Computing},
+	number = {3},
+	pages = {1103-1133},
+	title = {{\texttt{p4est}}: Scalable Algorithms for Parallel Adaptive Mesh Refinement on Forests of Octrees},
+	volume = {33},
+	year = {2011},
+	doi={10.1137/100791634}
+}
+
+@incollection{paraview,
+	author = {Ahrens, J and Geveci, B and Law, C},
+	booktitle = {Visualization Handbook},
+	publisher = {Elesvier},
+	title = {{ParaView}: An End-User Tool for Large Data Visualization},
+	year = {2005},
+	doi = {10.1016/b978-012387582-2/50038-1}
+}
+
+@Misc{            anl2023petsc,
+  author        = {Satish Balay and Shrirang Abhyankar and Mark~F. Adams and Steven Benson and Jed
+                  Brown and Peter Brune and Kris Buschelman and Emil~M. Constantinescu and Lisandro
+                  Dalcin and Alp Dener and Victor Eijkhout and Jacob Faibussowitsch and William~D.
+                  Gropp and V\'{a}clav Hapla and Tobin Isaac and Pierre Jolivet and Dmitry Karpeev
+                  and Dinesh Kaushik and Matthew~G. Knepley and Fande Kong and Scott Kruger and
+                  Dave~A. May and Lois Curfman McInnes and Richard Tran Mills and Lawrence Mitchell
+                  and Todd Munson and Jose~E. Roman and Karl Rupp and Patrick Sanan and Jason Sarich
+                  and Barry~F. Smith and Stefano Zampini and Hong Zhang and Hong Zhang and Junchao
+                  Zhang},
+  title         = {{PETS}c {W}eb page},
+  url           = {https://petsc.org/},
+  year          = {2024}
+}
+
+@article{martinsson2015hierarchical,
+	author = {Martinsson, PG},
+	journal = {arXiv preprint arXiv:1506.01308},
+	title = {{The hierarchical Poincar{\'e}-Steklov (HPS) solver for elliptic PDEs: A tutorial}},
+	year = {2015},
+	doi={10.48550/arXiv.1506.01308}
+}
+
+@incollection{quarteroni1991theory,
+	author = {Quarteroni, A and Valli, A},
+	booktitle = {Applied and Industrial Mathematics},
+	pages = {179--203},
+	publisher = {Springer},
+	title = {{Theory and application of Steklov-Poincar{\'e} operators for boundary-value problems}},
+	year = {1991},
+	doi = {10.1007/978-94-009-1908-2_14}
+}
+
+@incollection{HPV:VisIt,
+	author = {Hank Childs and Eric Brugger and Brad Whitlock and Jeremy Meredith and Sean Ahern and David Pugmire and Kathleen Biagas and Mark Miller and Cyrus Harrison and Gunther H. Weber and Hari Krishnan and Thomas Fogal and Allen Sanderson and Christoph Garth and E. Wes Bethel and David Camp and Oliver R\"{u}bel and Marc Durant and Jean M. Favre and Paul Navr\'{a}til},
+	booktitle = {High Performance Visualization--Enabling Extreme-Scale Scientific Insight},
+	month = {October},
+	pages = {357-372},
+	title = {VisIt: An End-User Tool For Visualizing and Analyzing Very Large Data},
+	year = {2012},
+	doi={10.1201/b12985-29}
+}
+
+@misc{ultraSEM,
+	author = {Fortunato, D and Hale, N and Townsend, A},
+	title = {ultraSEM: The ultraspherical spectral element method},
+	year = {2022},
+	publisher = {GitHub},
+	journal = {GitHub repository},
+	url = {https://github.com/danfortunato/ultraSEM}
+}
+
+@misc{HPS_Demos,
+	author = {Gillman, A},
+	title = {HPS_Demos: A collection of codes applying the HPS method},
+	year = {2023},
+	publisher = {GitHub},
+	journal = {GitHub repository},
+	url = {https://github.com/agillman20/HPS_Demos}
+}
+
+@misc{Streamer_HPS,
+	author = {Semenov, I},
+	title = {Streamer_HPS_DGSEM: Implementation of the spectral element method for modelling streamer discharges},
+	year = {2023},
+	publisher = {GitHub},
+	journal = {GitHub repository},
+	url = {https://github.com/igsemenov/Streamer_HPS_DGSEM}
+}
+

paper.md

@@ -0,0 +1,88 @@
+---
+title: 'EllipticForest: A Direct Solver Library for Elliptic Partial Differential Equations on Adaptive Meshes'
+tags:
+  - C++
+  - numerical linear algebra
+  - partial differential equations
+  - adaptive mesh refinement
+authors:
+  - name: Damyn Chipman
+    orcid: 0000-0001-6600-3720
+    affiliation: 1
+    corresponding: true
+affiliations:
+ - name: Boise State University, USA
+   index: 1
+date: 1 January 2024
+bibliography: paper.bib
+---
+
+# Summary
+
+EllipticForest is a software library with utilities to solve elliptic partial differential equations (PDEs) with adaptive mesh refinement (AMR) using a direct matrix factorization. It implements a quadtree-adaptive variation of the Hierarchical Poincaré-Steklov (HPS) method [@gillman2014direct]. The HPS method is a direct method for solving elliptic PDEs based on the recursive merging of Poincaré-Steklov operators [@quarteroni1991theory]. EllipticForest is built on top of the parallel and highly efficient mesh library `p4est` [@burstedde2011p4est] for mesh adaptivity and mesh management. Distributed memory parallelism is implemented through the Message Passing Interface (MPI) [@mpi41]. EllipticForest wraps the fast, cyclic-reduction methods found in the FISHPACK [@swarztrauber1999fishpack] library and updated in the FISHPACK90 [@adams2016fishpack90] library at the lowest grid level (called leaf patches). In addition, for more general elliptic problems, EllipticForest wraps solvers from the PDE solver library PETSc [@anl2023petsc]. The numerical methods used in EllipticForest are detailed in [@chipman2024fast]. A key feature of EllipticForest is the ability for users to extend the solver interface classes to implement custom solvers on leaf patches. EllipticForest is an implementation of the HPS method to be used as a software library, either as a standalone to solve elliptic PDEs or for coupling with other scientific libraries for broader applications.
+
+# Statement of Need
+
+Elliptic PDEs arise in a wide-range of physics and engineering applications, including fluid modeling, electromagnetism, astrophysics, heat transfer, and more. Solving elliptic PDEs is often one of the most computationally expensive steps in numerical algorithms due to the need to solve large systems of equations. Parallel algorithms are desirable in order solve larger systems at scale on small to large computing clusters. Communication patterns for elliptic solvers make implementing parallel solvers difficult due to to the global nature of the underlying mathematics. In addition, adaptive mesh refinement adds coarse-fine interfaces and more complex meshes that make development and scalability difficult. The solvers implemented in EllipticForest address these complexities through proven numerical methods and efficient software implementations.
+
+The general form of elliptic PDEs that EllipticForest is tailored to solve is the following:
+
+$$
+\alpha(x,y) \nabla \cdot \Big[ \beta(x,y) \nabla u(x,y) \Big] + \lambda(x,y) u(x,y) = f(x,y)
+$$
+
+where $\alpha(x,y)$, $\beta(x,y)$, $\lambda(x,y)$, and $f(x,y)$ are known functions in $x$ and $y$ and the goal is to solve for $u(x,y)$. Currently, EllipticForest solves the above problem in a rectangular domain $\Omega = [x_L, x_U] \times [y_L, y_U]$. The above PDE is discretized using a finite-volume approach using a standard five-point stencil yielding a second-order accurate solution. This leads to a standard linear system of equations of the form
+
+$$
+\textbf{A} \textbf{u} = \textbf{f}
+$$
+
+which is solved via the HPS method, a direct matrix factorization method.
+
+Similar to other direct methods, the HPS method is comprised of two stages: a build stage and a solve stage. In the build stage, a set of solution operators are formed that act as the factorization of the system matrix corresponding to the discretization stencil. This is done with $\mathcal{O}(N^{3/2})$ complexity, where $N$ is the size of the system matrix. In the solve stage, the factorization is applied to boundary and non-homogeneous data to solve for the solution vector with linear complexity $\mathcal{O}(N)$. The build and the solve stages are recursive applications of a merge and a split algorithm, respectively. The advantages of this approach over iterative methods such as conjugate gradient and multi-grid methods include the ability to apply the factorization to multiple right-hand side vectors.
+
+In addition, another advantage of the quadtree-adaptive HPS method as implemented in EllipticForest is the ability to adapt the matrix factorization to a changing grid. When the mesh changes due to a refining/coarsening criteria, traditional matrix factorizations must be recomputed. The quadtree-adaptive HPS method can adapt the factorization locally, eliminating the need to recompute the factorization. This works as the HPS method builds a set of solution operators that act like a global solution operator. When the mesh changes, the set can be updated by locally applying the merging and splitting algorithms. This is especially practical for implicit time-dependent problems that require a full linear solve each time step.
+
+EllipticForest builds upon the `p4est` mesh library [@burstedde2011p4est]. The quadtree-adaptive HPS method is uniquely suited for quadtree meshes. `p4est`, as a parallel and highly efficient mesh library, provides routines for creating, adapting, and iterating over quadtree meshes. The routines in EllipticForest wrap or extend the capabilities in `p4est`. A primary extension is the development of a *path-indexed* quadtree. This is in contrast to the *leaf-indexed* quadtree implemented in `p4est`. A *path-indexed* quadtree is a data structure that stores data at all nodes in a quadtree, as opposed to just the leaf nodes (see \autoref{fig:parallel_quadtree}). The *path-indexed* quadtree data structure is designed to store the various data and operators required in the quadtree-adaptive HPS method.
+
+The novelty of EllipticForest as software is the implementation of the HPS method for coupling with other scientific software as well as user extension. Currently, other implementations of the HPS method are MATLAB or Python codes designed by research groups and used in-house for solving specific problems [@ultraSEM; @HPS_Demos; @Streamer_HPS]. EllipticForest is designed to be extended and coupled with external libraries. This paper highlights the software details including the user-friendly interface to the HPS method and the ability for users to extend the solver interface using modern object-oriented programming (OOP) paradigms.
+
+# Software Overview
+
+Below, we outline various components of the software implemented in EllipticForest. These classes and utilities allow the user to create and refine meshes tailored for their use case, initialize the solver for the elliptic PDE, and visualize the output solution. A user may also extend the functionality of EllipticForest through inheritance of the `Patch` classes for user-defined solvers at the leaf level.
+
+## Quadtree
+
+The underlying data structure that encodes the mesh is a path-indexed quadtree. The `Quadtree` object is a class that implements a path-indexed quadtree using a `NodeMap`, which is equivalent to `std::map<std::string, Node<T>*>`. The template parameter `T` refers to the type of data that is stored on quadtree nodes. The `Quadtree` implemented in EllipticForest wraps the `p4est` leaf-indexed quadtree to create, iterate, and operate on the path-indexed quadtree. Functions to iterate over the quadtree include `traversePreOrder`, `traversePostOrder`, `merge`, and `split`. The `traversePreOrder` and `traversePostOrder` functions iterate over the tree in a pre- and post-order fashion, respectively, and provide the user with access to the node or node data via a provided callback function. The `merge` and `split` functions iterate over the tree in a post- and pre-order fashion, respectively, and provide the user with access to a family of nodes, or a group of four siblings and their parent node.
+
+![A path-indexed quadtree representation of a mesh. Colors indicate which rank owns that node. The nodes colored by gradient indicate they are owned by multiple ranks.\label{fig:parallel_quadtree}](assets/parallel_path_indexed_tree.png)
+
+## Mesh
+
+The user interfaces with the domain discretization through the `Mesh` class. The `Mesh` class has an instance of the `Quadtree` detailed above. `Mesh` provides functions to iterate over patches or cells via `iteratePatches` or `iterateCells`.
+
+`Mesh` also provides the user with an interface to the visualization features of EllipticForest. A user may add mesh functions via `addMeshFunction`, which are functions in $x$ and $y$ that are defined over the entire mesh. This can either be a mathematical function $f(x,y)$ that is provided via a `std::function<double(double x, double y)>`, or as a `Vector<double>` that has the value of $f(x,y)$ at each cell in the domain, ordered by patch and then by the ordering of patch grid. Once a mesh function is added to the mesh, the user may call `toVTK`, which writes the mesh to a parallel, unstructured VTK file format. See the section below on output and visualization for more information.
+
+## Patches
+
+The fundamental building block of the mesh and quadtree structures are the patches. A `Patch` is a class that contains data matrices and vectors that store the solution data and operators needed in the HPS method. A `Patch` also has an instance of a `PatchGrid` which represents the discretization of the problem. Each node in the path-indexed quadtree stores a pointer to a `Patch`.
+
+In EllipticForest, the patch, patch grid, patch solver, and patch node factory interfaces are implemented as a pure virtual interface for the user to extend. Internally, EllipticForest uses these interfaces to call the implemented solver or discretization. By default, EllipticForest implements a 2nd-order, finite volume discretization and solver. This implementation is found under `src/Patches/FiniteVolume` and each class therein implements the pure virtual interface of the `Patch`, `PatchGrid`, `PatchSolver`, and `AbstractNodeFactory` classes. Users may use the finite volume implementation shipped with EllipticForest, or they may implement different solvers to be used in the HPS method.
+
+## HPS Solver
+
+Once the mesh has been created and refined, and the patch solver has been initialized, solving the elliptic problem on the input mesh is done by creating an instance of the `HPSAlgorithm` class. The `HPSAlgorithm` class has member functions that perform the setup, build, upwards, and solve stages of the HPS method. As the HPS method is a direct method, once the build stage has been completed, the upwards and solve stages can be called without rebuilding the matrix factorization.
+
+## Output and Visualization
+
+Once the problem has been solved over the entire mesh, each leaf patch in the mesh has the solution stored in one of its data vectors, `vectorU`. This is a discrete representation of the solution to the PDE.
+
+The user may choose to output the mesh and solution in an unstructured PVTK format using the VTK functionality built-in. To output to VTK files, the user first adds mesh functions to the mesh. This includes the solution stored in `vectorU` after the HPS solve. Then, the user calls the `toVTK` member function of the `Mesh` class. This will write a `.pvtu` file for the mesh and a `.pvtu` file for the quadtree. An example of this output for a Poisson equation is shown in \autoref{fig:poisson_solution}.
+
+![Solution of Poisson equation on a quadtree mesh using EllipticForest. The mesh and data are output in an unstructured PVTK format and visualized with VisIt [@HPV:VisIt].\label{fig:poisson_solution}](examples/elliptic-single/output.png)
+
+# Acknowledgements
+
+The development of EllipticForest has been funded by the National Science Foundation (NSF-DMS #1819257) and the Boise State University School of Computing. The author acknowledges the assistance and guidance of Dr. Donna Calhoun and Dr. Carsten Burstedde through discussions and direction.
+
+# References