WIAS-PDELib · chmerdon · Jul 8, 2025 · Jul 8, 2025
diff --git a/docs/src/coefficients.md b/docs/src/coefficients.md
@@ -1,9 +1,19 @@
 # Stochastic Coefficients
 
-Coefficients are assumed to be represented by a Karhunen-Loeve expansion (KLE)
-that have the following general structure and API.
+Stochastic coefficients play a central role in uncertainty quantification and stochastic finite element methods. In this package, coefficients are represented using a Karhunen-Loève expansion (KLE), which allows for the efficient representation of random fields with prescribed covariance structure.
+
+## Overview
+
+The Karhunen-Loève expansion expresses a stochastic process as a series of orthogonal functions weighted by uncorrelated random variables:
+
+$$
+    a(x, \omega) = a_0(x) + \sum_{n=1}^N \sqrt{\lambda_n} \phi_n(x) \xi_n(\omega)
+$$
+
+where $a_0(x)$ is the mean, $\lambda_n$ and $\phi_n(x)$ are the eigenvalues and eigenfunctions of the covariance operator, and $\xi_n$ are independent standard normal random variables.
 
 ## API
+
 ```@autodocs
 Modules = [ExtendableASGFEM]
 Pages = ["coefficients/coefficients.jl"]

diff --git a/docs/src/estimators.md b/docs/src/estimators.md
@@ -1,14 +1,9 @@
 # Adaptivity and error control
 
 ## Residual-based a posteriori error estimation
-Spatial error estimation refers to classical residual-based error estimation
-for the zero-th multi index that refers to the mean value.
-Stochastic error control has to estimate which stochastic mode needs to be refined
-in the sense that either the polynomial degree is increased or neighbouring
-stochastic modes are activated. Both is represented by the multi-indices.
-Then unified error control allows to perform residual-based error estimation for
-the subresiduals that are associated to each multi-index and depend on the model
-problem (see references on the main page for details).
+Spatial error estimation refers to classical residual-based error estimation for the zero-th multi-index, which corresponds to the mean value.
+Stochastic error control determines which stochastic mode should be refined, either by increasing the polynomial degree or by activating neighboring stochastic modes. Both are represented by multi-indices.
+Unified error control enables residual-based error estimation for the subresiduals associated with each multi-index, depending on the model problem (see references on the main page for details).
 
 
 ```@autodocs
@@ -17,12 +12,10 @@ Pages = ["estimate.jl"]
 Order   = [:type, :function]
 ```
 
-## Multi index management
+## Multi-index management
 
-Depending on the model problem and stochastic coefficient the
-amount of multi indices that should be added to the error estimation
-varies.
-Here are some methods that help with enriching the set of multi-indices.
+Depending on the model problem and stochastic coefficient, the number of multi-indices to be added for error estimation varies.
+The following methods assist in enriching the set of multi-indices.
 
 
 ```@autodocs
@@ -31,12 +24,10 @@ Pages = ["mopcontrol.jl"]
 Order   = [:type, :function]
 ```
 
-## Monte carlo sampling estimator
+## Monte Carlo sampling estimator
+
+A hierarchical Monte Carlo error estimator is also available. It compares the solution with a higher-order discrete solution for sampled deterministic problems. This is mainly intended to compute a reference error to assess the efficiency of the residual-based error estimator.
 
-There is also a hierarchical Monte carlo error estimator available that
-compares the solution with a higher order discrete solution for sampled
-deterministic problems. This is merely intended as a way to compute the
-reference error to assess the efficiency of the residual-based error estimator.
 
 ```@autodocs
 Modules = [ExtendableASGFEM]

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -1,26 +1,22 @@
 # ExtendableASGFEM
 
-This package implements the stochastic Galerkin finite element method for certain two dimensional
-model problems involving KLE of stochastic coefficients. The rather huge systems have a tensorized
-structure and are solved by iterative solvers. A posteriori error estimators steer the spatial
-and stochastic refinement.
-
-The spatial discretization is based on the finite element package [ExtendableFEM.jl](https://github.com/WIAS-PDELib/ExtendableFEM.jl)/[ExtendableFEMBase.jl](https://github.com/WIAS-PDELib/ExtendableFEMBase.jl)
+This package provides an implementation of the stochastic Galerkin finite element method (SGFEM) for selected two-dimensional model problems involving Karhunen-Loève expansions (KLE) of stochastic coefficients. The resulting large-scale systems exhibit a tensorized structure and are efficiently solved using iterative solvers. Adaptive a posteriori error estimators guide both spatial and stochastic refinement to ensure accuracy and efficiency.
 
+Spatial discretization is performed using the finite element packages [ExtendableFEM.jl](https://github.com/WIAS-PDELib/ExtendableFEM.jl) and [ExtendableFEMBase.jl](https://github.com/WIAS-PDELib/ExtendableFEMBase.jl).
 
 ## References
 
-- [1]   "Adaptive stochastic Galerkin FEM"\
-        M. Eigel, C.J. Gittelson, C. Schwab, E. Zander\
-        CMAME 270, 1 (2014), 247--269\
-        [>Journal-Link<](https://www.doi.org/10.1016/J.CMA.2013.11.015)
-        [>Preprint-Link<](https://www.research-collection.ethz.ch/handle/20.500.11850/154941)
-- [2]   "A posteriori error control for stochastic Galerkin FEM with high-dimensional random parametric PDEs",\
-        M. Eigel, C. Merdon\
-        to be published in: Error Control, Adaptive Discretizations, and Applications, Part 3, Academic Press\
-        [>Preprint-Link<](https://www.wias-berlin.de/publications/wias-publ/run.jsp?template=abstract&type=Preprint&year=&number=3174)
-- [3]   "Local equilibration error estimators for guaranteed error control in adaptive higher-order stochastic  Galerkin finite element methods"\
-        M. Eigel and C. Merdon\
-        SIAM/ASA J. Uncertainty Quantification 4(1) (2016), 1372--1397"\
-        [>Journal-Link<](https://epubs.siam.org/doi/10.1137/15M102188X)
-        [>Preprint-Link<](http://www.wias-berlin.de/publications/wias-publ/run.jsp?template=abstract&type=Preprint&year=2014&number=1997)
+- [1]   "Adaptive stochastic Galerkin FEM"  
+        M. Eigel, C.J. Gittelson, C. Schwab, E. Zander  
+        CMAME 270, 1 (2014), 247–269  
+        [Journal-Link](https://www.doi.org/10.1016/J.CMA.2013.11.015)  
+        [Preprint-Link](https://www.research-collection.ethz.ch/handle/20.500.11850/154941)
+- [2]   "A posteriori error control for stochastic Galerkin FEM with high-dimensional random parametric PDEs"  
+        M. Eigel, C. Merdon  
+        To appear in: Error Control, Adaptive Discretizations, and Applications, Part 3, Academic Press  
+        [Preprint-Link](https://www.wias-berlin.de/publications/wias-publ/run.jsp?template=abstract&type=Preprint&year=&number=3174)
+- [3]   "Local equilibration error estimators for guaranteed error control in adaptive higher-order stochastic Galerkin finite element methods"  
+        M. Eigel and C. Merdon  
+        SIAM/ASA J. Uncertainty Quantification 4(1) (2016), 1372–1397  
+        [Journal-Link](https://epubs.siam.org/doi/10.1137/15M102188X)  
+        [Preprint-Link](http://www.wias-berlin.de/publications/wias-publ/run.jsp?template=abstract&type=Preprint&year=2014&number=1997)
diff --git a/docs/src/index.pdf b/docs/src/index.pdf
diff --git a/docs/src/onbasis.md b/docs/src/onbasis.md
@@ -1,10 +1,6 @@
 # ONBasis
 
-An ONBasis (=orthonormal basis) stores information for the orthogonal polynomials
-of the disctribution, like norms, quadrature rules and cached evaluations at quadrature points.
-It is the main building brick for the tensorized basis associated to the multi-indices
-of the stochastic discretization.
-
+An ONBasis (orthonormal basis) encapsulates information about the orthogonal polynomials associated with a given probability distribution. This includes norms, quadrature rules, and cached evaluations at quadrature points. The ONBasis serves as a fundamental building block for constructing the tensorized basis linked to the multi-indices in the stochastic discretization.
 
 ```@autodocs
 Modules = [ExtendableASGFEM]

diff --git a/docs/src/onbasis.pdf b/docs/src/onbasis.pdf
diff --git a/docs/src/orthogonal_polynomials.md b/docs/src/orthogonal_polynomials.md
@@ -1,41 +1,34 @@
 # Orthogonal Polynomials
 
-The stochastic discretization of random variables involves global
-polynomials that are orthogonal with respect to the involved
-random distribution of the `y_m`. These polynomials
-can be generated by some recurrence relation with coefficients
-that depend on the distribution.
+In the stochastic discretization of random variables, global polynomials that are orthogonal with respect to the probability distribution of each random variable $y_m$ are used. These orthogonal polynomials can be generated via recurrence relations, with coefficients determined by the underlying distribution.
+
+## Recurrence Relations
+Orthogonal polynomials $H_n$ with respect to a weight function $\omega$ satisfy
 
-## Recurrence relations
-Orthogonal polynomials
-$H_n$ with respect to some weight function $\omega$, i.e.,
 ```math
-\int_\Gamma \omega(y) H_{n}(y) H_m(y) dy = N^2_{nm}\delta_{nm}
+\int_\Gamma \omega(y) H_{n}(y) H_m(y) \,dy = N^2_{nm}\delta_{nm}
 ```
-where 
+where
 ```math
-  N_{nn} := \| H_n \|_{\omega}^2 := \int_\Gamma \omega(y) H_{n}(y) H_n(y) dy
+N_{nn} := \| H_n \|_{\omega}^2 := \int_\Gamma \omega(y) H_{n}(y) H_n(y) \,dy
 ```
-satisfy the three-term recurrence relation
+The polynomials satisfy the three-term recurrence relation:
 ```math
 \begin{aligned}
-  H_{n+2}(y) & = (a_n y-b_n) H_{n+1}(y) - c_n H_{n}(y)
-  \end{aligned}
+  H_{n+2}(y) & = (a_n y - b_n) H_{n+1}(y) - c_n H_{n}(y)
+\end{aligned}
 ```
-initialized by $H_0  = 0$ and $H_1 = 1$.
-
+with initial values $H_0 = 0$ and $H_1 = 1$.
 
 ```@autodocs
 Modules = [ExtendableASGFEM]
 Pages = ["orthogonal_polynomials/orthogonal_polynomials.jl"]
 Order   = [:type, :function]
 ```
 
-## Legendre Polynomials (uniform distribution)
+## Legendre Polynomials (Uniform Distribution)
 
-For the weight function $\omega(y) = 1/2$ in the interval $[-1,1]$ (uniform distribution),
-take $a_n = (2n+1)/(n+1)$, $b_n = 0$ and $c_n = n/(n+1)$. The norms
-of the resulting Legendre polynomials are given by
+For the weight function $\omega(y) = 1/2$ on the interval $[-1,1]$ (uniform distribution), the recurrence coefficients are $a_n = (2n+1)/(n+1)$, $b_n = 0$, and $c_n = n/(n+1)$. The norms of the resulting Legendre polynomials are
 ```math
     \| H_n \|^2_\omega = \frac{2}{2n+1}
 ```
@@ -46,21 +39,20 @@ Pages = ["orthogonal_polynomials/Legendre_uniform.jl"]
 Order   = [:type, :function]
 ```
 
+## Hermite Polynomials (Normal Distribution)
 
-## Hermite Polynomials (normal distribution)
-
-For the weight function $\omega(y) = \exp(-y^2/2)/(2\pi)$ (normal distribution), take $a_n = 1$, $b_n = 0$ and $c_n = n$. Then, the first six polynomials read
+For the weight function $\omega(y) = \exp(-y^2/2)/(2\pi)$ (normal distribution), the recurrence coefficients are $a_n = 1$, $b_n = 0$, and $c_n = n$. The first six Hermite polynomials are:
 ```math
 \begin{aligned}
 H_0 & = 0\\
 H_1 & = 1\\
 H_2 & = y\\
 H_3 & = y^2 - 1\\
 H_4 & = y^3 - 3y\\
-H_5 & = y^4 - 6y^2 +3\\
+H_5 & = y^4 - 6y^2 + 3\\
 \end{aligned}
 ```
-and their norms are given by
+Their norms are given by
 ```math
     \| H_n \|^2_\omega = n!
 ```

diff --git a/docs/src/orthogonal_polynomials.pdf b/docs/src/orthogonal_polynomials.pdf
diff --git a/docs/src/solvers.md b/docs/src/solvers.md
@@ -1,30 +1,21 @@
 # Solver
 
-Solving requires a spatial and stochastic discretization.
-Both are connected in a special vector structure
-that is passed to a solve function that runs a special
-iterative solver for each model problem.
+Solving a problem requires both spatial and stochastic discretization. These are combined into a specialized vector structure, which is then passed to a solver function that executes an iterative algorithm tailored to each model problem.
 
 ## SGFEVector
 
-The spatial discretization is represented by
-s single finite element space from [ExtendableFEM.jl](https://github.com/WIAS-PDELib/ExtendableFEM.jl),
-while the stochastic discretization is represented by a tensorized basis
-for the parameter space of the stochastic coefficient. Both have to be prepared in
-advance.
+The spatial discretization is defined by a single finite element space from [ExtendableFEM.jl](https://github.com/WIAS-PDELib/ExtendableFEM.jl), while the stochastic discretization uses a tensorized basis for the parameter space of the stochastic coefficient. Both components must be set up in advance.
 
 !!! note
-
-    Currently it is not possible to use different finite element spaces for different multi-indices,
-    but this feature might be added in the future.
+    Currently, it is not possible to use different finite element spaces for different multi-indices. This feature may be added in the future.
 
 ```@autodocs
 Modules = [ExtendableASGFEM]
 Pages = ["sgfevector.jl"]
 Order   = [:type, :function]
 ```
 
-## Solve function
+## Solve Function
 
 ```@autodocs
 Modules = [ExtendableASGFEM]

diff --git a/docs/src/tonbasis.md b/docs/src/tonbasis.md
@@ -1,19 +1,8 @@
 # TensorizedBasis
 
-Each multi-index ``\mu = [\mu_1,\mu_2,\ldots,\mu_M]``
-encodes a tensorized basis function for the parameter space
-of the form ``H_\mu = \prod_{k=1}^M H_k`` where the
-``H_k`` are the orthogonal polynomials.
-The TensorizedBasis collects all information necessary
-to evaluate those basis functions, i.e. the set of multi-indices
-and the triple products of the form ``(y_mH_\mu, H_\lambda)``
-for each ``m`` and ``\mu, \lambda`` in the set of multi-indices
-as a sparse matrix. There are analytic formulas to evaluate
-these triple products in terms of recurrence coefficients, but it makes
-sense to store them for faster evaluation times.
-
-
+Each multi-index $\mu = [\mu_1, \mu_2, \ldots, \mu_M]$ defines a tensorized basis function for the parameter space of the form $H_\mu = \prod_{k=1}^M H_k$, where each $H_k$ is an orthogonal polynomial.
 
+The `TensorizedBasis` object collects all information required to evaluate these basis functions, including the set of multi-indices and the triple products of the form $(y_m H_\mu, H_\lambda)$ for each $m$ and $\mu, \lambda$ in the set of multi-indices, stored as a sparse matrix. While analytic formulas exist to compute these triple products using recurrence coefficients, storing them can significantly speed up evaluations.
 
 ```@autodocs
 Modules = [ExtendableASGFEM]

diff --git a/docs/src/tonbasis.pdf b/docs/src/tonbasis.pdf
diff --git a/src/ExtendableASGFEM.jl b/src/ExtendableASGFEM.jl
@@ -44,7 +44,7 @@ export get_neighbourhood_relation_matrix
 include("orthogonal_polynomials/orthogonal_polynomials.jl")
 export OrthogonalPolynomialType
 export HermitePolynomials, LegendrePolynomials #, LaguerrePolynomials, ChebyshevTPolynomials, ChebyshevUPolynomials
-export evaluate, gauss_rule
+export evaluate, evaluate!, gauss_rule
 export norm, norms
 export distribution
 export normalise_recurrence_coefficients

diff --git a/src/estimate.jl b/src/estimate.jl
@@ -38,12 +38,28 @@ end
 """
 $(TYPEDSIGNATURES)
 
-computes the residual-based a posteriori error estimator for the current solution vector `sol`
-for the given model problem and stochastic coefficient `C` and returns:
-- `eta4modes` = array with total error estimators for each multi-index (corresponding to the enriched set of multi-indices with the current active modes coming first)
-- `eta4cell` = array with total error estimators for each cell in grid (for spatial refinement)
-- `multi_indices_extended` = enriched set of multi-indices used for the computations
-
+Compute the residual-based a posteriori error estimator for a stochastic Galerkin solution vector `sol`.
+
+# Arguments
+- `::Type{AbstractModelProblem}`: The model problem type for which the estimator is called. Used for dispatching to the appropriate estimator implementation.
+- `sol::SGFEVector`: The current stochastic Galerkin solution vector.
+- `C::AbstractStochasticCoefficient`: The stochastic coefficient (random field or parameterization).
+- `rhs`: (Optional) Right-hand side function for the PDE (default: `nothing`).
+- `bonus_quadorder`: (Optional) Additional quadrature order for integration (default: 1).
+- `tail_extension`: (Optional) Number of additional boundary modes to include in the multi-index set (default: 5).
+- `kwargs...`: Additional keyword arguments passed to the estimator.
+
+# Returns
+A tuple containing:
+- `eta4modes::Vector{Float64}`: Total error estimator for each multi-index (stochastic mode), corresponding to the enriched set of multi-indices (with current active modes first).
+- `eta4cell::Matrix{Float64}`: Error estimator for each cell in the spatial grid (for spatial refinement), for each multi-index.
+- `multi_indices_extended`: The enriched set of multi-indices used in the computation (including boundary extensions).
+
+# Description
+This function computes a residual-based a posteriori error estimator for the current SGFEM solution. It supports both spatial and stochastic adaptivity by providing error indicators for each cell and each stochastic mode. The estimator is tailored to the model problem and the stochastic coefficient, and can be extended to include additional boundary modes for improved reliability.
+
+If no specialized estimator is available for the given model problem type, an error is raised and empty arrays are returned.
+```
 """
 function estimate(::Type{AbstractModelProblem}, sol::SGFEVector, C::AbstractStochasticCoefficient; kwargs...)
     @error "no error estimator for the model problem type available"

diff --git a/src/modelproblems/modelproblems.jl b/src/modelproblems/modelproblems.jl
@@ -3,15 +3,15 @@ abstract type AbstractModelProblem end
 """
 $(TYPEDSIGNATURES)
 
-solves the specified model problem with the given stochastic coefficient `C` and right-hand side `rhs`
-and writes the solution into `sol`. Via this `sol` vector the spatial and stochastic discretization
-is communicated as well as initial data for the iterative solver.
-The boolean `use_iterative_solver` (default is true) determines if the iterative solver is used
-or if the full matrix is assembled and solved by a direct solver (very slow for larger systems).
-The parameters `bonus_quadorder_f` (default is 0) and `bonus_quadorder_a` (default is 2) can be used
-to increase the quadrature order in terms that involve the rhs or the stochastic coefficient, respectively.
+Solves the specified model problem using the provided stochastic coefficient `C` and right-hand side `rhs`, writing the solution into `sol`.
 
+The `sol` vector communicates both the spatial and stochastic discretization, as well as any initial data required for the iterative solver.
 
+- If `use_iterative_solver` (default: `true`) is set, an iterative solver is used. Otherwise, the full system matrix is assembled and solved directly (note: this is very slow for large systems).
+- The parameters `bonus_quadorder_f` (default: `0`) and `bonus_quadorder_a` (default: `2`) allow you to increase the quadrature order for terms involving the right-hand side or the stochastic coefficient, respectively.
+- Additional keyword arguments can be passed via `kwargs`.
+
+If no solver is implemented for the given model problem, an error is thrown.
 """
 function solve!(
         ::Type{AbstractModelProblem},