Convex Optimization-Based Structure-Preserving Filter for Multidimensional Finite Element Simulations]Structure-preserving convex optimization for multidimensional finite element simulations^††thanks: Accepted by editors (JCP) 07/08/2023.

Many sophisticated approaches have been considered in the literature to tackle different aspects of the structure-preservation problem. We make a broad classification into ‘intrusive’ and ‘nonintrusive’ classes of solutions. Many popular approaches generally involve changing either the solution or the PDE and require the imposition of stringent conditions on solver parameters, such as the step-size of the numerical method or the domain shape and granularity. Such methods can be classified as intrusive. On the other hand, some methods constrain the solution obtained from the solver with minimal, often superficial, change to the numerical scheme, i.e., in a nonintrusive way.

Prominent among the intrusive class of approaches is the transformation of structure preservation into a PDE-constrained optimization problem, methods that modify the spatial discretization [30], and limiters derived from Karush-Kahan-Tucker (KKT) optimality conditions [29]. Another example of the intrusive approach [21] proposes an operator-splitting positivity-preservation method based on the energy dissipation law.

The strategies for constraint satisfaction in [3] consider a version of the problem that is a specialization of the one solved by the solution proposed in this manuscript, which falls under the nonintrusive class. In [2, 1] and extensions [3, 23], the authors explore approaches for positivity preservation on the basis functions derived from Bernstein polynomials such that they are non-negative and possess a partition of unity property. Therefore, the interpolated solution respects the original bounds at any point in the domain, nonintrusively. Other nonintrusive methods prescribe solutions such as limiters or truncation and modification of the domain (e.g., curvilinear mesh adaptation) to adhere to the desired solution structures. One such method is established in [4] where the authors present a new approach for multi-material arbitrary Lagrangian–Eulerian (ALE) hydrodynamics simulations based on high-order finite elements posed on high-order curvilinear meshes in which conservative fields are remapped and additional synchronization steps are introduced for structural preservation. Additionally, many nonintrusive methods perform pre- or postprocessing of intermediate solutions without modifying the underlying problem to incorporate the constraints or making changes to the domain. Some nonintrusive approaches to structure preservation can be found in [20, 34, 25, 22, 27]. In [33], the authors present a positivity-preserving limiter by maintaining cell averages within a certain tolerance. The work in [16] imposes positivity on a discrete grid for kinetic equations via a nonlinear filtering procedure. Satisfying constraints up to a numerical tolerance is sufficient for many applications. However some applications require the solution to adhere to a strict non-negative structure. For example, in [13], the authors enforce strict, global non-negativity of the diffusion propagator by formulating constraints specific to the propagator model. The algorithm in [28] employs alternating least squares procedures to implement general linear equality and inequality constraints by reducing it to a non-negativity-constrained least squares (NNLS) problem. The non-negativity constraints based on square root representations have also been of interest in many other works [14, 10]. We propose a nonintrusive method based on a general filtering approach in [31, 32].

Many of the existing techniques prescribe approaches to preserve the structure that work only for low-order problems. Additionally, an increase in the dimensionality of the problem further exacerbates the issue of accuracy and convergence. Many approaches struggle to maintain high-order accuracy and convergence robustly, including the filter’s ability to preserve structure without any additional dependence or restrictions imposed on a per-case basis. This includes problem-specific modifications as well as varying solutions based on other aspects of the problem. For example, [15, 26] presents redistribution-based local and global repair methods depending on a sweep through cells to locate and rectify the violating structures. The final solution depends on the sweep order, producing different solutions for different sweep orders. In addition, this repair procedure preserves element-averaged bound/positivity constraints at particular finite difference nodes. It does not preserve them pointwise, e.g., at every spatial location. The proposed method tries to preserve the constraints at every spatial location. An ideal approach would achieve a balance between robust structure preservation, maintaining accuracy as well as convergence, and providing usability by deterministic termination within an acceptable time.

We pose the problem of preserving the structure in polynomial-based numerical methods as a convex optimization problem. Our previous works discuss a filtering algorithm to solve this problem and demonstrate the solution for function approximations [31] and PDE solutions in 1D [32]. The robust design of the problem enables the preservation of different properties of the solution simultaneously (e.g., positivity, monotonicity, and boundedness). It transforms the problem into a composite constraint satisfaction problem that can be solved by convex minimization. To address this problem, we propose a filter that can be applied as a postprocessing step to deal with structural inconsistencies. The novel idea in the design of the filter is the efficiency achieved by applying the filter only to the violating regions in a subdomain and structure-preservation continuously throughout the domain, and not just a subset of points within the domain. This guarantee of structure-preservation depends on the granularity of the global minimization procedure in 2D and 3D, which is expensive and less reliable compared to its 1D counterpart. In 1D, the minimum can be found accurately using the confederate matrix approach to find zeros of the derivative of a polynomial. For finding the critical points in higher dimensions, an efficient gradient-based method is required which solves the nonconvex problem of finding the global minimum. In the case of the proposed filter, the global minimum represents a structure violation. Gradient-based methods are prone to get stuck at a local minimum, causing it to miss the point of larger structural inconsistency, thus failing to provide strong guarantees of preservation. Another important aspect that affects the outcome of the gradient-based method is the choice of a starting point. To address this aspect, we propose an investigation of structural inconsistencies on a lattice of points as the initial step in the gradient-based method. The point with the most inconsistency is chosen as a starting point. The accuracy of the minimum found depends on the pattern and granularity of the initial lattice used, which is detailed in Section 2.1.

To streamline the minima-finding stage and reduce the time and space complexity, we adopt the Barycentric interpolation approach from [6], expanded upon by [17]. Using these techniques, a considerable speed-up in time taken by minimization and overall filter computations is achieved.

This paper focuses on designing a structure-preserving filter using a convex optimization approach for 2D and 3D problems on different element types. We will present the formulation derived from [31] in Section 2 and extend the concepts to more complex problems tackled in this paper. The remainder of the paper is organized as follows: Section 2 discusses the setup for applications of the filter to time-dependent PDEs in 2D and 3D. Section 3 introduces the notations, details the design of different building blocks of the filter, and summarizes the 1D problem formulation from [31, 32]. A procedure for the filtering process developed for multidimensional applications is presented in Section 4. Section 5 describes the numerical results to demonstrate the filter’s efficacy in preserving the desired structures in different application scenarios. This section is divided into subsections describing the process of choosing the parameters to run the experiments, the 2D and 3D canonical experiments, and an advection problem on different homogeneous meshes using the discontinuous Galerkin (dG) formulation. We conclude with a demonstration of the use of the proposed filter on a real-world application: the model of platelet aggregation and blood coagulation problem.

If the polynomial projections lose structural conformity, they are likely to do so at or between the quadrature points. We hope to capture the structural inconsistencies by checking the values on a lattice of quadrature points and the centroids formed by adjacent quadrature points. Here, the arbitrary choice of the centroids is one of the many possible options to choose the lattice. A quadrilateral is used as a sample element to describe the setup; however, the same applies to any canonical element type in 2D and 3D. Let us call such a lattice of points $(X)$. For a quadrilateral, let $Q$ represent the number of quadrature points in one dimension. $X$ is defined by a combination of the quadrature grid on the quadrilateral (total $Q^{2}$ points and $(Q-1)^{2}$ centroids formed by the midpoints of the quadrature grid) as defined in Equation 2.1.1.

For a triangle, we can follow the same idea of combining the quadrature points and the centroids to construct the lattice. An example lattice on a mesh with quadrilaterals and triangles is shown in Figure 2.1.1.

Consider a function in 2D defined by Equation 2.1.2.

where

Projecting Equation 2.1.2 using polynomial order $N=5$ on a quadrilateral in the domain $\Omega=[0,1]\times[0,1]$, we get the coefficients $\boldsymbol{\tilde{v}}=\{\tilde{v}_{j}\}_{j=0}^{P-1}$. The projection is shown in the left subfigure of Figure 2.1.2. In this case, $P=(N+1)^{d}$, where $d$ is dimension. Since $d=2$, $P=36$. The original function is non-negative, so the projection should preserve positivity. To this end, we employ the optimization outlined in Section 4 as a postprocessing filter. The process behaves as a spectral filter based on the empirical evidence presented in [31, Section 5.2] and further theoretical explanation provided in [31, Proof 5.1]. We will hereafter refer to this optimization as a nonlinear filter ($\mathcal{F}$). Let $\boldsymbol{{v}}$ be the filtered version of $\boldsymbol{\tilde{v}}$. The filtered projection is shown in the right subfigure of Figure 2.1.2.

The advection equation for a quantity described by a scalar field $u$ is expressed mathematically by a continuity equation Equation 2.2.1 in the domain $\Omega\subset{\mathbbm{R}}^{d}$.

where $\boldsymbol{a}$ represents the velocity of advection. Assume proper initial and boundary conditions are specified.

Galerkin-type methods assume an ansatz for $u$ as a time-varying element of a fixed $P$-dimensional linear subspace $V$ that contains the polynomial functions for a fixed degree $P$, where frequently $V\subset L^{2}(\Omega)$.

where $x\in\Omega$, and $\{\phi_{j}\}_{j=0}^{P-1}$ represents the traditional FEM basis, which is not necessarily orthogonal. Let $\{\psi_{j}\}_{j=0}^{P-1}$ represent a collection of orthonormal basis functions We can transform the vector of coefficients $\boldsymbol{\widehat{v}}$ into its orthonormal form $\boldsymbol{\tilde{v}}$.

Consider a partition of $\Omega$ consisting of $E$ subdomains defined by $\tau(\Omega)=\{e_{0},e_{1},\cdots,e_{E-1}\}$. The discontinuous Galerkin formulation assumes that $V$ is comprised of functions that are polynomials of a fixed degree $P$ on each element on $\Omega$; where discontinuities in derivatives are allowed only at partition boundaries. The semidiscrete form for Equation 2.2.1 is derived in the standard Galerkin way by using the ansatz Equation 2.2.2 and forcing the residual to be $L^{2}$-orthogonal to $V$. Usually, integration by parts is performed in the residual orthogonalization step, and often depending on the equation and spatial discretization, numerical flux and/or stabilization terms are included in the weak formulation.

The result is a system of ordinary differential equations prescribing time-evolution of the discrete degrees of freedom represented by vector $\boldsymbol{\tilde{v}}=\{\tilde{v}_{0},\ldots,\tilde{v}_{(P-1)}\}.$

where $\boldsymbol{M}and\boldsymbol{A}$ are the $P\times P$ mass and advection matrices, respectively, defined as

and F is a general term to compensate for any numerical fluxes or stabilization terms. Finally, Equation 2.2.3 is transformed into a fully discrete system that can be solved using an appropriate numerical integration scheme.

Consider the case when the solution defined by $\widehat{v}^{n}$ on a particular partition $e$ at a timestep $n$ does not satisfy the desired structural properties. To preserve the structure, we employ the optimization defined in Equation 3.1.2 as a postprocessing filter. Let $\boldsymbol{v}$ represent the filtered solution that satisfies the structure. The process can be illustrated in Equation 2.2.4. The proposed filter works by simply augmenting this scheme nonintrusively.

Since DG discretizations have degrees of freedom that are decoupled across the subdomains, the total degrees of freedom in the optimization step need to be the only ones over the subdomain. In addition, the optimization needs to be applied only in the subdomains that have structural violations discovered at the lattice points, which leads to a parallelizable set of independent optimizations each in $P$ dimensions. We can filter the selected subdomains simultaneously, which further adds to the procedure’s efficiency.

For a mesh with the number of elements $E>1$, the process of enforcing positivity per element leads to changes in solution properties; in particular, elementwise boundary values and the mean of the discrete solution are not preserved. Assume that numerical fluxes are computed as explicit functions of element boundary values. A shift in element boundary values by the filter causes changes in the corresponding fluxes, thus adding errors to the simulation. We can resolve this issue by imposing additional equality constraints (i.e., function values at element boundaries) in the filter. However, for the flux conservation, we cannot make an assertion that the resulting filtered solutions should satisfy the original flux/jump conditions or some updated conditions, based on the new filtered basis coefficients. In the case of mass (integral) conservation, additional equality constraints can be imposed. The values of fluxes are often used to ensure the properties of numerical schemes therefore by preserving fluxes one can attest that the properties of the original scheme have been retained. Preserving fluxes in this context corresponds to corrections in one element need not have an impact on neighboring elements, and one way to promote isolation of the optimization effect to the element under consideration is to ensure that interelement communication through fluxes is undisturbed. The formulation for flux conservation in [32, Section 3.2 and 3.3] is unique to one-dimensional scenarios. Generalization to higher dimensions is not immediate. Additional analysis is needed to incorporate different elements and, therefore, beyond the scope of this paper. The robust design of the filter allows for the incorporation of an arbitrary number of structural constraints into one optimization problem. However, increasing the number of constraints reduces the degrees of freedom for the filter accordingly. In particular, if the filter has $P$ degrees of freedom and we have $q$ linear equality constraints, then the dimension of the optimization problem can be reduced to $(P-q)$ using the procedure described in [32, Section 3.2].

Let $u(x)\in V$ be a function, where $V$ is a finite-dimensional subspace of a Hilbert space of real-valued functions on $\Omega\in\mathbbm{R}^{d}$. Let V contain a collection of orthonormal basis functions $\{\psi_{n}\}_{n=1}^{N}$ for some $N\in\mathbbm{N}$.

Therefore, we have

where $\left\langle\cdot,\cdot\right\rangle$ is the inner product on $V$, and $\delta_{ij}$, the Kronecker delta function.

The value of $d$ is set to 1 for simplicity. The formulation follows the same steps for higher dimensions. A comprehensive constraint-satisfaction problem can be constructed as detailed in this section.

Let $K$ be the collection of $u$ that satisfies the constraints. Therefore, $K$ is the feasible subset of $V$. An example of such a family of linear constraints is {positivity, boundedness}. Consider a family of linear constraints such as the one in Equation 3.1.1, where ${L}$ is a linear operator bounded on $V$, and $\ell$ is a function on the domain $\Omega$.

Note that the framework in [31] allows a finite number of such constraints to be considered simultaneously.

According to the Riesz representation theorem, if $V^{*}$ is a dual of $V$, then a functional $L\in V^{*}$ can be associated with a unique V-representor $\ell\in V$ satisfying

Furthermore, this $L\leftrightarrow\ell$ identification is an isometry. We will use these facts in what follows. Given $L$ that identifies $\ell$, we consider the coordinates $\widehat{\ell}_{j}$ of $\ell$ in a $V$-orthonormal basis,

Then we have the following relations:

where $\|\cdot\|_{2}$ is the Euclidean norm on vectors in $\mathbbm{R}^{N}$.

The problem is essentially that of developing a map $\mathcal{F}$ from a function $u\in V$, to a unique function $\mathcal{F}(u)\in K$. Initially, $u$ may not satisfy the structural constraints but its mapped version $\mathcal{F}(u)$ does. In [31] we introduce a strategy to solve the following well-posed convex feasibility problem.

Under certain conditions, such as if $0\in K$, $\mathcal{F}$ is norm-contractive; therefore, the operation can be called a filter. For a brief discussion of the norm-contractive nature of Equation 3.1.2, see [31, Proposition 5.1].

We first transform the continuous problem Equation 3.1.2 to a discrete version to construct a feasible solution. Let $\{\psi_{j}\}_{j=0}^{P-1}$ be a collection of orthonormal basis functions in $V$, different from the standard FEM basis ($\phi$), which are not orthonormal.

Consider $C$ as the affine conic region in $\mathbbm{R}^{P}$ corresponding to $K\subset V$. From previously established notations, $\boldsymbol{\tilde{v}}$ is a vector of expansion coefficients of $u$ that does not satisfy that desired constraints, and $\boldsymbol{v}$ is a vector of filtered coefficients of $u$ that satisfies the desired constraints. Representing the Euclidean 2-norm on vectors as $\|\cdot\|_{2}$, and the fact that $C\subset V$, we get the equivalent of the optimization problem Equation 3.1.2 as follows:

For a fixed $x\in\Omega$, let $H_{x}$ be a $(P-1)$-dimensional planar surface defined by an equality constraint corresponding to Equation 3.1.1 for a fixed $x$. Therefore, $L(u)=\ell(x)$. Writing $u$ as an expansion in degrees of freedom $\tilde{v}$ that corresponds to a single linear equality constraint for the $\tilde{v}$, i.e., a $(P-1)$-dimensional plane. In a geometric sense, $H_{x}$ is a surface representing the constraint boundary in the space $\mathbbm{R}^{P}$, dividing it into two hyperspaces as shown in Figure 3.2.1. One hyperspace represents the region that satisfies a linear inequality constraint and the other that does not. For a single linear constraint, $C$ can be defined in terms of hyperspaces as Equation 3.1.1.

To obtain the full feasible set $C$, the filter iteratively projects $\tilde{\boldsymbol{v}}$ onto the intersection of feasible halfspaces defined by individual hyperplanes $\{H_{x}\}$. To describe the feasible regions further, let us consider a geometric interpretation. Let the current state vector of the projection defined by $\tilde{\boldsymbol{v}}$ be a point shown by the blue dot in Figure 3.2.1. The filter works by applying correction Equation 3.2.7 to $\tilde{\boldsymbol{v}}$ by computationally inspecting the signed distance function Equation 3.2.4.

For the positivity-preservation example, the process of “inspection” means the ability to compute the global minimum of $s(x)$ to determine a region where $s$ is negative. Based on this inspection, the algorithm projects the state vector of the current iterate $\boldsymbol{\tilde{v}}$ onto $H_{y}$ for some $y\in\Omega$. This projection can cause a particular $H_{y1}$ with a positive $s(y1)$ to go negative. Therefore, we need to repeat the projection step until the solution $\boldsymbol{\tilde{v}}$ lies completely in (or on the boundary of) $C$. [31, Algorithm 1] summarizes all the steps of the filter.

At each iteration, the spatial point $x$ is found that minimizes the signed distance function given by Equation 3.2.5. Note that for more than one constraint we need to keep track of $x$ and the constraint specific to the minimum signed distance function.

This greedy procedure is an infinite-constraint generalization of Motzkin’s relaxation method [24]. The major computational work in the filtering procedure is to minimize the objective defined by Equation 3.2.6 at each iteration.

Next, update $\boldsymbol{\tilde{v}}$ by projecting it onto the hyperplane corresponding to $x^{\ast}$:

where $\boldsymbol{h}(x,)$ is the normal vector corresponding to hyperplane $H_{x}$ of a single constraint family that points toward $C$. This vector is readily computable from the orthonormal basis, see [31] for details. This procedure is repeated until $s(x^{\ast})$ vanishes to within a numerical tolerance.

Finding the violating hyperplane that is farthest away corresponds to finding the global minimum of $s(x)$ on $\Omega$. Unlike the 1D procedure described in [31], multidimensional global minimization requires a more expensive approach such as gradient descent (GD). In this paper, we implement 2D and 3D global minimization using line-search-based GD with backtracking. We further provide a strategy to investigate the times taken for convergence and the accuracy of this GD approach by varying the values of the parameters. Then, we choose the best values for the parameters for the accuracy of the method. Section 4.1 provides details of the GD variant chosen for filter implementation and corresponding parameter selection to efficiently solve the minimization problem.

An important and the most expensive part of the optimization process is the global minimization Equation 3.2.5. For efficiency and usability of the filter, we need to fine-tune the most expensive step, i.e., minimizing the scaled distance function Equation 3.2.6, which quantifies the difference between the current state of the coefficients and the coefficients that satisfy the constraints at a critical point of the current iteration. The success of the proposed filter depends on accurately finding the points in the domain that violate the structure. Since the problem formulation uses dG, it can be solved element by element on a domain with $E$ elements. Applying the filter at each timestep is essentially solving $E$ optimization problems at every timestep. Each optimization problem is over a single element. We can reduce the number of optimization problems by testing the solution for structural inconsistencies at a lattice of points defined by Equation 2.1.1. For example, in the case of positivity preservation, if the minimum value of an element is nonpositive, the filter is applied to that element.

In the 1D case, the filter finds the minimum value by evaluating the solution at the critical points on the element. If the subspace $V$ restricted to an element is a space of polynomials, then the critical points of Equation 3.2.6 are defined by the roots of the derivative of Equation 3.2.6. For the root-finding problem, we employ the confederate matrix-generation approach. Relative to 1D, minimization in higher dimensions is more complex and expensive. Many variants of the gradient descent-based (GD-based) methods have been proposed to perform this operation efficiently and accurately. These methods operate only up to a certain tolerance and, like all the GD-based methods, have the potential for stagnation.

Since there is no multidimensional version of the confederate linearization approach that exactly solves the first-order optimality conditions, the problem of finding the critical points on the multidimensional domain is nonconvex. Therefore, we cannot ensure that all such points will be flagged for correction by any particular minimization algorithm. From our analysis, one suitable GD approach for this job is to use adaptive descent size and the ability to retrace the steps. For this reason, we choose the method of steepest descent with a backtracking line search strategy. The backtracking line search starts with a relatively large step size and repeatedly shrinks it by a factor $\gamma$ until the Armijo–Goldstein condition Equation 4.1.1 is fulfilled. This widely used algorithm makes an informed choice about step size and direction at each iteration to efficiently arrive at the optimum value. It avoids stagnation by tracing back its steps if the difference between descent values falls below a fixed tolerance. An essential part of this algorithm is choosing the parameters denoted by $c\in(0,1)$ and $\gamma\in(0,1)$ that determine the step size and backtracking performed by the algorithm. For an objective function $f$, with a starting position $x$, given parameters $c$ and $\gamma$ the Armijo–Goldstein condition is defined as Equation 4.1.1.

where $\nabla f(x)^{T}$ is the slope in the direction of descent $p$ and $\gamma$ is calculated as

$c$ and $\gamma$ play a vital role in testing the condition [5], which determines whether a step-wise movement from a current position to a modified position achieves an adequately corresponding decrease in the objective function. Although many recommendations exist [5, 11, 12, 7, 8] for an optimal way to choose $c$ and $\gamma$, these parameters must be tuned for maximum efficiency.

From the discussion in Section 2, the optimization that enforces a constrained solution is applied as a postprocessing part inside the timestepper, independent of the choice of the timestepping routine. As predicted by the ‘curse of dimensionality’, structure preservation and optimization for higher dimensional problems are more complex than their 1D counterparts. Since the central idea of the filter is agnostic to the dimensions of the problem, the Algorithm 1 remains the same as described in [31]. Similarly, the applications to PDE solutions follow the procedures established in [32] for 1D problems.

Conversion from an input coefficient vector $\widehat{\boldsymbol{v}}$ to corresponding coefficients $\tilde{\boldsymbol{v}}$ in an orthonormal basis is the first stage of the filtering procedure. Algorithm 1 presents a summary of steps taken by the PDE solver for the filter application.

A significant chunk of work in the iterative correction stage is done by evaluating the solution at the lattice points using basis interpolation, an expensive operation that involves processing and storage of large interpolation matrices. To make the processing computationally efficient, we use the Barycentric interpolation method from [6], expanded upon by [17]. This method reduces the number of calculations and storage required to evaluate the basis functions at arbitrary points in the domain, thereby reducing the cost of the iterative correction stage. Although the implementation details of Barycentric interpolation method are beyond the scope of this paper, for the minimization part of the numerical experiments in Section 5, we observe a significant speed-up using the Barycentric approach.

To fix the values for $c$ and $\gamma$ Equation 4.1.1 for the numerical results in Section 5, we first choose a set of functions. The choice of the functions depends on factors such as the constraint type and the domain geometry. For example, for an application where positivity is the primary structural concern, a list of functions with values close to zero is a good choice. For applications involving discontinuities in the function or the domain, discontinuous functions would be a good choice. The feasible space for $c$ and $\gamma$ is $(0,1)$. A non-comprehensive sampling on $(0,1)$ is performed for both parameters, and the GD is called with the parameters set to permutations of the sampled values. GD_linesearch denotes the call to the GD with line search. As its output, the GD_linesearch routine provides the number of iterations (denoted by $niter$) taken by the algorithm. The error refers to the difference between the minimum value found and the known minimum value.

As a large number of GD iterations per filter iteration makes the filtering process more expensive, choosing the parameter values for which GD converges in the least number of iterations is desirable. Another important quantity to consider here is the error. It is desirable to have GD return with a minimum value as close to the known (or golden) minimum value as possible. With these quantities ($niter$ and $err$) as the metrics, we prescribe a procedure to select the ideal value of $c$ and $\gamma$ in Algorithm A.1 [Appendix Appendix A]. The selection criteria narrow the samples to a range of $c$ and $\gamma$ with the minimal number of iterations and then pick from that subset the ones with the smallest error. In the case of a tie, all the $c$ and $\gamma$ values are included in the return variables. The procedure can be summarized as taking a list of functions that fairly represent the complexity of space and returning two ranges of ideal values: one for each of the parameters $c$ and $\gamma$, respectively. Since the quadrilateral and hexahedron represent richer space than their other canonical 2D and 3D counterparts, Algorithm A.1 [Appendix Appendix A] is applied and analyzed only on a quadrilateral for 2D and a hexahedron for 3D.

Future developments in the approaches for finding optima in multiple dimensions may improve the efficacy and efficiency of the filter. To this end and to provide robustness, we employ the object-oriented principle by modularizing the minimization part of the filter.

We consider different example functions Section 5.1 on the 2D domain $[-1,1]\times[-1,1]$ and Section 5.1 on the 3D domain $[-1,1]\times[-1,1]\times[-1,1]$, which have varying complexities and projection orders to set the GD parameters comprehensively and fairly.

Using the orthogonal basis for Equation 3.2.7, $\boldsymbol{\tilde{{v}}}$ is obtained, which is then used by Algorithm A.1 to find the values of parameters $c$ and $\gamma$. For functions that have an exact analytical minimum, calculating the error between the exact minimum and the optimum returned by the GD is straightforward. For the discontinuous functions such as $f_{2}$ and $f_{5}$, an approximate golden minimum is calculated by projecting the function using polynomial order $N=8$ on a dense grid of points in $\Omega$ and finding the least value of the projected polynomial. In all our test cases, GD converges to the exact (or golden) minimum within an acceptable tolerance (numerical 0 set to $10^{-7}$). Therefore, the selection of the parameters is made based on the fewest iterations taken by GD ($niter$) to converge.

Consider $k$ discrete equispaced samples of $\in(0,1)$. For $k=9$, by applying Algorithm A.1, we get the results shown in Table 5.1.1.

Based on the results of these tests, we infer that for the numerical experiments, the appropriate choices for $c$ and $\gamma$ in 2D are 0.7 and 0.7, respectively. Following a similar procedure, Algorithm A.1 on a hexahedron for functions defined in Section 5.1, the choices of $c$ and $\gamma$ for 3D are 0.2 and 0.7, respectively. These values are used for all the numerical experiments in the numerical results.

The C++ implementation of the filter in Nektar++ [9] supports change in $c$ and $\gamma$ as parameters in the configuration file.

Consider a 2D function Equation 5.2.1 and a 3D function Equation 5.2.2 that are both discontinuous clamped versions of a smooth sinusoidal function. The initial projections of $f(x,y)$ on a quadrilateral and $f(x,y,z)$ on a hexahedron are shown in Figure 5.2.1. The discontinuity produces oscillations similar to Gibb’s phenomenon upon projection, which is interesting for the application and analysis of the filter.

where

where

It is evident from the left subfigures in Figure 5.2.1 that the projection of discontinuous non-negative functions leads to the negative values as shown in the highlighted regions. As seen in the right subfigures of Figure 5.2.1, the application of the structure-preserving filter restores the non-negative structures of Equation 5.2.1 and Equation 5.2.2, respectively. To analyze the p-convergence, the experiment is repeated for different polynomial orders, as shown in Figure 5.2.2.

Figure 5.2.2 presents $L_{2}$ error comparison between the filtered and unfiltered version of projections of nonsmooth functions in 2D and 3D. The choice of a nonsmooth function is to demonstrate the effect of filter application in worst-case scenarios.

Consider a dG solution to the advection problem Equation 5.3.1 on a 2D domain as shown in the first row of Figure 5.3.1.

Consider a function $f$ defined on $[-1,1]\times[-1,1]$ as Equation 5.3.2.

For a given initial condition to be $f(x,y)$, $\textbf{a}=[1,1]$, and the periodic boundary conditions, we formulate the problem in the discontinuous Galerkin framework. The first step is projecting $f$ on a set of $E$ elements $\{e_{0},e_{1}\cdots e_{E-1}\}\in\Omega$ using the typical nonorthogonal (hats and bubbles) basis $\phi$. Locally, for an element $e$, we have

where $\Omega$ is the mesh shown in Figure 5.3.1 consisting of a mix of quadrilaterals and triangles. The filter steps do not change for different element types as emphasized by the choice of composite mesh.

The solution for the advection problem using the timestep $\Delta t=1e-3$, polynomial order 4, RK-4 integration scheme, and upwind flux calculations is shown in Figure 5.3.1. Row 2 of the figure shows the simulation state at a particular timestep and highlights the negative values in the domain. In row 3 the non-negative structure is restored after the application of the filter, The filter is applied at each timestep to preserve the structure (positivity) of the solution on a lattice of points of interest. The lattice is a set of points defined in Equation 2.1.1. If a structure violation is found at any point in the lattice, the parent element is flagged for filtering. Since the boundaries are periodic, the final state of the simulation looks similar to the initial state.

For a similar analysis of the 3D advection problem, consider the initial state as a smooth continuous function Equation 5.3.3 on a cube mesh of 64 hexahedron elements defined in the domain $[-1,-1,-1]\times[1,1,1]$. The advection velocity is defined by a $=[1,1,1]$, and all the boundary conditions are periodic. The integration method used is RK-4 with a timestep of $1e{-3}$. The timestepping is performed for a total of 2000 timesteps and the flux calculation is upwind. Following a similar selection procedure as in 2D, we discover the elements that violate the structure at each timestep and apply the structure-preserving filter to those elements.

Figure 5.3.2 shows an instance during the advection process where a loss of structure is encountered, i.e., negative intermediate values are found, as shown by the highlighted region in the left and middle subfigures of Figure 5.3.2. The process of discovering the structurally nonconformal elements and application of the filter to those elements adds to the total time cost of the advection solution, which is an important aspect to consider when deciding on the criteria for applying the filter.

The percentage increase in total time taken by the experiment by the application of structure-preserving filter is shown in Figure 5.3.3. Note that for 2D, at low orders (N=2,3,4), the cost of filtering is between 50-70$\%$ of the total time. This is due, in part, to the efficiency at these orders of the unfiltered solver such as caching effects, as well as up-front costs associated with the optimization. We find that for 2D the ratio of filtered to unfiltered is lowest at N=5, and then proceeds to climb linearly as shown in the left subfigure of Figure 5.3.3. The time taken by the filter in 3D is notably higher than the time taken by the filter in 2D because of the higher complexity of finding the global minimum in 3D. The cost of overall filtering per timestep depends on the number of elements that the filter operates upon. Therefore, the procedure of selective application of the filter becomes increasingly important.

We now investigate the application of the filter on the solid body rotation experiment using a discontinuous initial condition defined as a combination of a notched cylinder, a cone, and a smooth hump. Consider the initial data shown in Figure 5.4.1. This example is extensively used in advection and structure preservation literature [19]. The parameters to reproduce the initial state: the cylinder, cone, and hump have a radius of 0.3 each. The centers of the cylinder, cone, and hump are (0, 0.5), (0,-0.5), and (-0.6,0) respectively. Without changing the domain and boundary conditions used for the 2D advection experiment shown in Figure 5.3.1, the advection velocity is changed to circular such that in one time period, the solution returns to its original state. For all tests in this section, the domain has 144 quadrilateral elements and the timestep $(\Delta t)$ is $5e-4$.

Since the initial condition for the experiment is defined by a complicated discontinuous function, the solution suffers from a larger loss of structure as compared to the case with a smooth initial condition. Figure 5.4.2 shows the initial and final states during rotation of the body shown in Figure 5.4.1, with and without the filter application.

Figures 5.4.3 and 5.4.4 show different cross sections of the solution with and without filtering at the initial and final time. The total time taken by the filter as a function of the entire simulation time is shown in Figure 5.4.5.

As evident from the top left panel in Figure 5.4.5, the filter takes proportionally longer to restore the structure as the loss of structure increases. In many practical applications, such as the one presented in Section 5.6, we observe that paying the extra filtering cost results in the successful termination of the simulation, which otherwise fails due to the presence of structural inconsistencies. In such cases, the tradeoff between the time cost vs guarantee of structure preservation leans towards the latter. The filtered and unfiltered $L_{2}$ errors after one rotation of the solid body are shown in the top right panel of Figure 5.4.5. Whereas filtering changes the final state of the solution, thus contributing to the errors, the difference between filtered and unfiltered errors is still comparable. The filtered version of the tests using different polynomial order follows the same p-convergence as the unfiltered counterpart. A comparison between the number of times the filter is invoked and the number of iterations taken per invocation is shown in the bottom left panel of Figure 5.4.5. We notice that the number of iterations taken per call to the filter increases with the increase in polynomial order because of the larger magnitude of structural inconsistencies observed as a result of the oscillations in the solution. The largest contributor to the cost of applying the filter occurs from the GD line-search to find the global minimum. At each iteration of the filter, we find the minimum. Each call to find the minimum, in turn, takes a few iterations of GD. On the bottom right panel of Figure 5.4.5 is a comparison between the average number of GD iterations per iteration of the filter. Also shown is the average number of filter iterations per call to the filter.

The procedure to filter the elements that exhibit loss of structure remains agnostic to the type of element. To demonstrate this, consider different tesselations of the domain $[-1,-1,-1]\times[1,1,1]$ and solve an advection problem with the nontrivial initial state Equation 5.5.1. The setup remains the same as the previous 3D example, i.e., a $=[1,1,1]$ in all directions, periodic boundaries, RK-4 integration scheme using timestep $1e{-3}$, and for a total of 2000 timesteps.

The initial state looks like a torus inside a cube. The advection experiment is repeated on different homogeneous meshes with element types hexahedra, tetrahedra, and prisms individually. Figure 5.5.1 shows the results before and after the filter application at a particular timestep for all these meshes.

The analysis of the filtering effects on the accuracy and convergence of the advection process, especially on the meshes with hexahedron and prism elements, reveals that the difference between the $L_{2}$ errors is too small to warrant a convergence comparison. However, for the mesh with tetrahedral elements, the $L_{2}$ errors perceptibly vary at each timestep, as shown by a longer run of the same experiment (total 4000 timesteps) in Figure 5.5.2.

A prominent scenario in which the problem of structure-preservation is encountered is in the mathematical-biology domain. A detailed mathematical model for platelet aggregation and blood coagulation (PAC) by [18] is considered for this experiment. The model describes the process of evolution, interaction, and decay of the chemical species involved in the process of thrombosis.

Although the details of the individual species in the PAC process are beyond the scope of this paper, [18] has a detailed explanation of the nature and the evolution of the chemical species. The model can be summarized as a comprehensive collection of ODEs and PDEs that track all the chemical species in various phases during the process. The results in [18] use a finite-difference method with a specified limiter to truncate nonpositive concentration values to zero. In order to implement this model using the FEM without loss of structure, an alternative approach to truncation, such as the one presented in this paper, is needed. For our experiments, we use the structure-preserving filter instead of a truncation limiter. The focus is limited to the system of equations to track the chemical species that advect, diffuse, and react. One such species is fluid-phase thrombin (FP $e_{2}$), an essential component of the thrombosis process. For this section, we track the evolution of thrombin by advection, diffusion, and decay in the fluid-phase. For the detailed results of the evolution of all the chemical species, including thrombin, under various circumstances in the original model, refer to [18].

We set up a version of the PAC model that solves the species evolution problem by a combination of advection, diffusion, and reaction (ADR) PDEs. Unlike the original model, the velocity $\boldsymbol{u}$ of the fluid medium is not reported by a modified Navier-Stokes solver. Instead, the velocity is constant $\boldsymbol{u}=[u_{x},u_{y}]=[5,0]$. To solve this system of PDEs on a sample blood vessel represented by a rectangular domain $\Omega=[0,300]\times[0,20]$, we employ the dG FEM. The domain is tessellated using 2048 quadrilaterals. The bottom wall has an injury site spanning from $x=20$ to $x=60$. We use polynomial order $=4$ and timestep of $\Delta t=1e{-2}$ to solve the problem.

The total number of species tracked in our implementation of the model is 56. Since the focus is on the behavior of thrombin, we present the PDE describing the evolution of fluid-phase thrombin (FP $e_{2}$) Equation 5.6.1. For details on the other chemical species involved in Equation 5.6.1, refer to [18]. The boundary conditions vary depending on the chemical being tracked. For FP $e_{2}$, we have a Dirichlet zero boundary on the left, and a no-flux boundary on the right and top. The bottom boundary has two kinds of conditions: a robin boundary condition on the injury site and Dirichlet zero everywhere else.

Figure 5.6.1 shows filtered and unfiltered runs of our implementation of this model and the resulting concentration profile of FP $e_{2}$. At timestep = 2, we get the invalid (negative) values around the injury site because of the sharp concentration changes (refer to the slice at $y=5$ in the bottom left part of Figure 5.6.1). Therefore, the fluid phase thrombin concentration does not adhere to the desired structure (positivity), and the simulation does not succeed. If this experiment is allowed to continue, the effect of invalid values will compound over time, rendering the simulation results nonphysical and useless. At timestep 38, the simulation starts becoming unstable because of structural discrepancies. If it is allowed to run further, the accumulation and propagation of negative values results in the simulation blow-up at timestep 87, at which point the code reports invalid values such as NaN (datatype: Not a Number).

When the filter is applied, it runs for significantly more timesteps (total 674) and terminates naturally. The right side of Figure 5.6.1 shows the state of FP $e_{2}$ at the depletion point of the chemicals involved.