Next Article in Journal
Temperature Patterns in TSA for Different Frequencies and Material Properties: A FEM Approach
Previous Article in Journal
An Experimental Study of Grouping Mutation Operators for the Unrelated Parallel-Machine Scheduling Problem
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

An Adaptive in Space, Stabilized Finite Element Method via Residual Minimization for Linear and Nonlinear Unsteady Advection–Diffusion–Reaction Equations

1
School of Electrical Engineering, Computing and Mathematical Sciences, Curtin University, GPO Box U1987, Perth, WA 6845, Australia
2
Mineral Resources, Commonwealth Scientific and Industrial Research Organisation (CSIRO), Kensington, Perth, WA 6152, Australia
*
Author to whom correspondence should be addressed.
Math. Comput. Appl. 2023, 28(1), 7; https://doi.org/10.3390/mca28010007
Submission received: 7 October 2022 / Revised: 22 December 2022 / Accepted: 26 December 2022 / Published: 6 January 2023
(This article belongs to the Special Issue Discontinuous Galerkin Methods)

Abstract

:
We construct a stabilized finite element method for linear and nonlinear unsteady advection–diffusion–reaction equations using the method of lines. We propose a residual minimization strategy that uses an ad-hoc modified discrete system that couples a time-marching schema and a semi-discrete discontinuous Galerkin formulation in space. This combination delivers a stable continuous solution and an on-the-fly error estimate that robustly guides adaptivity at every discrete time. We show the performance of advection-dominated problems to demonstrate stability in the solution and efficiency in the adaptivity strategy. We also present the method’s robustness in the nonlinear Bratu equation in two dimensions.

1. Introduction

The unsteady advection–diffusion–reaction model system poses distinct challenges for its numerical approximation. A limit case of interest arises when the equation becomes advection-dominated, showing sharp internal or boundary layers. Classical numerical methods (e.g., standard FEM) lead to numerical instabilities when the mesh is not sufficiently fine to capture the fine scales associated with these flow features, leading to unphysical oscillation.
Stabilized methods, such as Stream Upwind Petrov Galerkin (SUPG) [1] or Galerkin least squares (GLS) [2], overcome this problem by including artificial diffusion in the governing equation’s variational form. The variational multi-scale method (VMS) [3,4] captures fine-scale features by modifying the variational form [5,6,7,8]. Other technics such as the discontinuous Galerkin (dG) methods [9,10,11,12,13] stabilize the solution by providing local conservation and adding inner penalization across the element interfaces.
Methods based on residual minimization, including Least-Squares Finite Element Methods (LSFEM) and the Discontinuous Petrov–Galerkin method (DPG), seek stability by minimizing the discrete residual in dual norms [14,15,16]. Extensions to these ideas regarding parabolic problems include: [17,18,19,20] for dG and [14,21,22,23,24] for DPG methods.
Although these methods show stability for advection-dominated problems, the lack of a priori localization of the inner or boundary layers in the exact solution leads to expensive simulations on quasi-uniform meshes. Additionally, the lack of robust refinement strategies is critical for unsteady problems where the solution varies in space and time. Thus, we reduce this computational cost by using adaptive methods that rely on a posteriori error estimators to refine solution singularities. Posteriori error estimators for unsteady diffusion advection reaction methods are described in [25,26] and for unsteady dG implementations in [27,28,29].
Recently, Calo et al. [30] introduced a new class of adaptive stabilized conforming finite elements via residual minimization for steady problems. The method combines residual minimization ideas with the stability of the discontinuous Galerkin formulations. As a result, the method delivers a stable continuous solution and a robust error representation to perform on-the-fly adaptivity. The authors introduce the framework in a series of papers for linear and nonlinear applications (e.g., advection–diffusion problems with heterogeneous and highly anisotropic diffusion [31], its use combined with isogeometric analysis [32], nonlinear weak constraint enforcement [33], goal-oriented adaptivity [34], for incompressible flows [35,36], flow in porous media [37], and dynamic fracture propagation [38]).
In this paper, we extend [30] for unsteady advection–diffusion–reaction problems using the method of lines. Our method offers robust spatial refinements for a user-selected time marching method. We first approximate the spatial derivatives using a space semi-discrete scheme and then solve the resulting system using a time-marching discretization. As particular examples, this paper uses implicit first- and second-order time-stepping (BDF1 and BDF2) discretizations [39,40,41].
Compared to other techniques, the main advantage of this method relies on the non-conformity of the starting dG formulation, which allows us to work with stronger norms from the dG theory with a continuous trial space. Moreover, the refinement strategy and its efficiency in obtaining high-resolution approximations from coarse meshes allow us to overcome the computational cost resulting from the implicit temporal schemes and the extra degrees of freedom in the saddle point formulation.
The paper’s outline follows: Section 2 introduces the model problem. Section 3 describes some preliminary concepts for dG discretizations for time marching. Then, we present the well-posedness of the dG method combined with the backward differentiating formula for time marching. Section 4 describes the residual minimization problem and introduces our adaptive stabilized finite element method for parabolic problems. Finally, Section 5 contains some numerical examples showing uniform and adaptive refinement cases for two-dimension linear and non-linear problems, followed by some concluding remarks.

2. Model Problem

Let Ω R d , with d = 2 , 3 be an open, bounded Lipschitz polygon with boundary Γ . We denote n as the outward normal vector to Γ . For a given open and bounded domain K, we represent its L 2 inner product and L 2 norm as ( · , · ) 0 , K and · 0 , K , respectively. We set ( · , · ) 0 : = ( · , · ) 0 , Ω and · 0 : = · 0 , Ω for convenience. We define the well-known Hilbert space H 1 ( Ω ) : = v L 2 ( Ω ) : v L 2 ( Ω ) with the inner product on Ω denoted by ( · , · ) 1 and the space H 0 1 ( Ω ) : = v H 1 ( Ω ) : v = 0 on Γ .
For any T > 0 , l > 0 , and let V be a Hilbert space, we denote by C l ( V ) : = C l ( 0 , T ; V ) the l-times continuously-differentiable function space in [ 0 , T ] . Thus, C 0 ( V ) and C 1 ( V ) represent the continuous and continuously differentiable space function in [ 0 , T ] , respectively.
The inflow (−) and outflow(+) subsets of the boundary Γ are defined by
Γ : = x Γ | β · n < 0 , Γ + : = x Γ | β · n 0 .
where β L i p ( Ω ) (Lipschitz continuous) represents a velocity vector field.
We denote by Γ D and Γ N the Dirichlet and Nuemann boundary, respectively, such that Γ = Γ N Γ D ¯ . Thus, we define the inner and outer part of the Neumann boundary as follows:
Γ N : = Γ N Γ , Γ N + : = Γ N Γ + .
We consider the time evolution of the advection–diffusion–reaction solution defined in the space–time cylinder Ω × ( 0 , T ] for T > 0 . The governing equations in strong form read:
t u · κ u + β · u + μ u = f in Ω × ( 0 , T ] , u = g D on Γ D × ( 0 , T ] , ( κ u β u ) · n = g N on Γ N × ( 0 , T ] , κ u · n = g N on Γ N + × ( 0 , T ] , u ( · , t = 0 ) = u 0 ( x ) in Ω ,
where μ L represents the reaction coefficient, κ ( x ) > 0 the diffusivity, f C 0 ( L 2 ( Ω ) ) the source term, and g D C 0 ( H 1 / 2 ( Γ D ) ) and g N C 0 ( L 2 ( Γ N ) ) the Dirichlet and Neumann boundary values. We assume that β , κ and μ are time-independent, and that ( β · v + μ v , v ) 0 0 for all v H 0 1 ( Ω ) . Denoting L β , the Lipschitz modulus of β , we consider a reference velocity β c and a reference time τ c defined respectively as: β c : = β and τ c : = max μ , L β 1 .
Introducing the billinear form
a ( u , v ) : = ( κ u , v ) 0 + ( β · u , v ) 0 + ( μ u , v ) 0 ,
the weak form of (1) then reads: find u L 2 ( 0 , T ; H 1 ( Ω ) ) H 1 ( 0 , T ; L 2 ( Ω ) ) such that, for each t ( 0 , T ] ,
( t u , v ) 0 + a ( u , v ) = ( f , v ) 0 , v H 0 1 ( Ω ) .

3. Discontinuous Galerkin-Based Time Marching Discretization

3.1. Discrete Setting

We set T h as the triangulation of Ω , and K an element of T h . We define the finite dimensional spaces:
V h ( T h ) : = v L 2 ( T ) : v | K P b ( K ) , K T h
and
U h ( T h ) : = V h ( T h ) C 0 ( Ω ) ,
where P b ( K ) denotes the set of functions with degree lower or equal than b on K. Given K 1 and K 2 T h two disjoint adjacent elements in T h , sharing an internal face F = K 1 K 2 , we define n F as the normal vector on the face F from K 1 to K 2 (see Figure 1). We define the set of all faces as S h : = K T h F and the internal and boundary faces set by S h 0 : = S h \ Γ and S h = S h Γ , respectively. Moreover, we denote by S h D : = S h Γ D the set of Dirichlet boundary faces and by S h N : = S h Γ N the set of Neumann boundary faces. Let h K be the element diameter of K T h and h F be the face diameter of F S h . Given a face F S h 0 , we define the jump and average of v across F by
[ [ v ] ] F ( x ) : = v | K 1 ( x ) v | K 2 ( x ) x F
and
{ v } F ( x ) : = 1 2 ( v | K 1 ( x ) + v | K 2 ( x ) ) x F .
If F S h , we set [ [ v ] ] F ( x ) = { u } F ( x ) : = v | K ( x ) x F .

3.2. Space Semi-Discretization

We formulate the space semi-discretization by combining the Symmetric Interior Penalty (SIP) and the upwind dG formulations (UPW) for the steady advection–diffusion-reaction equation. We set that h K β c min ( T , τ c ) to avoid strong reaction regimes, to allow the mesh to resolve the spatial variation of the velocity field, and to guarantee that a particle at speed β c crosses at least one mesh element over the time interval ( 0 , T ) . Let the semi-discrete dG approximation (1) in the space V h be: For t ( 0 , T ] , find θ h , such that
( t θ h , v h ) 0 + a h ( θ h , v h ) = h ( v h ) v h V h ,
with θ h ( 0 ) = θ 0 . We define the advection–diffusion-reaction billinear form a h as follows:
a h ( u , v ) : = a h ( u , v ) SIP + a h ( u , v ) UPW ,
with
a h ( u , v ) SIP : = K T ( κ h u , h v ) K + F S h 0 n e κ ( [ [ u ] ] , [ [ v ] ] ) F F S h 0 ( κ h u · n F , [ [ v ] ] ) F + ( [ [ u ] ] , κ h v · n F ) F F S h D ( κ h u · n F , v ) F + ( u , κ h v · n F ) F n e κ ( u , v ) F
and
a h ( u , v ) UPW : = K T μ u + β · h u , v K + F S h Γ ( β · n F u , v ) F F S h 0 ( β · n F ) [ [ u ] ] , v + n a 2 | β · n F | [ [ u ] ] , [ [ v ] ] F ,
where n e and n a are two positive penalty coefficients for the diffusion and advection billinear forms. We define n e explicitly as (see [30]):
n e : = n o ( p + 1 ) ( p + d ) d 1 2 A ( K 1 ) V ( K 1 ) + A ( K 2 ) V ( K 2 ) , if F = K 1 K 2 A ( K ) V ( K ) , if F = K Γ ,
where n o > 0 is a user-defined constant, p is the polynomial degree of the test space and V and A represent the volume and area of an element in 3D, and its length and area in 2D, respectively. Moreover, n a modifies the numerical flux associated with the upwinding billinear form. The centered fluxes correspond to n a 0 while the upwind fluxes correspond to n a 1 . We set n a = 1 and n o = 1 for this research. In the general case when weakly non-homogeneous boundary conditions are enforced, the linear form h ( v h ) for the discrete problem (5) reads:
h ( v ) : = ( f h ( t ) , v ) + F S h D n e κ ( g D , v ) F ( g D , κ h v · n F ) F + F S h D Γ β · n F g D , v F + F S h N ( g N , v ) F .
We set f h ( t ) = π h f ( t ) t [ 0 , T ] , where π h is the L 2 projection onto V h . Next, we manipulate functions of the form ( u ( t ) v h ) in the space V h : = H 1 ( Ω ) + V h . Thus, we can write an equivalent form of (5) in terms of the discrete differential operator A h : V h V h , such that, for all ( u , v h ) V * h × V h ,
( A h u , v h ) 0 : = a h ( u , v h ) .
We use the discrete operator A h to formulate the space semi-discrete problem (12) in the form: for each t ( 0 , T ] , then
( t θ h ( t ) , v ) 0 + ( A h θ h ( t ) , v ) 0 = h ( v ) in v H 0 1 ( Ω ) ,
with the initial condition θ h ( 0 ) = π h u ( 0 ) . We endow V h with the norm:
v V 2 : = v SIP 2 + v UPW 2 ,
where v SIP 2 and v UPW 2 correspond to the symmetric interior penalty (SIP) and upwinding (UPW) norms defined as follow:
v SIP 2 : = κ v 0 2 + F S h n e κ v 0 , F ,
and
v UPW 2 : = τ c 1 v 0 2 + 1 2 | β · n F | ( v , v ) 0 , Γ + F S h 0 1 2 | β · n F | ( v , v ) 0 , F + K T β c 1 h K β · v 0 2 .
Additionally, we define its extended norm as:
v V , * 2 : = v V 2 + K T β c v 0 , Γ 2 + h K 1 v 0 2 + K T h K κ v · n 0 , Γ . 2
We introduce the discrete properties of the operator A h following [§ 3–4] [42]
Theorem 1 
(Discrete operator A h properties).
1. 
Consistency: The exact solution u of (1) satisfies
t u ( t ) + A h u ( t ) = h ( t ) t ( 0 , T ] .
2. 
Boundedness: There is a constant C b n d < , independent of h and τ, such that
A h v , w h 0 C b n d v V , * w h V ( v , w h ) V h * × V h .
3. 
Discrete inf-sup stability: There is a constant C s t a > 0 , such that
C s t a v h V sup w h V h { 0 } a h ( v h , w h ) w h V v h V h .

3.3. Backward Euler Time Discretization

We first consider the Backward Euler method (BDF1) for time marching and implement the second-order Backward differentiation formula (BDF2) in Section 3.4. Herein, the ≲ symbol denotes less or equal to a mesh-independent constant. We define τ : = T / N as the time step, where T is the final time, and N is a positive integer. We set τ min ( T , τ c ) . We use the following first-order approximation of the time derivative:
δ t ( 1 ) v n + 1 : = v n + 1 v n τ V n 0 , N .
Thus, the fully discrete problem is: for n = 0 , , n 1 , find θ h n + 1 V h , such that
( δ t ( 1 ) θ h n + 1 , v h n + 1 ) 0 + ( A h θ h n + 1 , v h n + 1 ) 0 = ( h n + 1 , v h n + 1 ) 0 v h n + 1 V h ,
where θ h 0 = π h u 0 and h n + 1 denotes the discrete linear form on V h * (10) at time n + 1 . We define the discrete-time operator A h , τ : V * h V h , such that, for all ( u , w h ) V h * × V h ,
( A h , τ u , w h ) = a h , τ ( u , w h ) : = ( u , w h ) 0 + τ a h ( u , w h ) .
Thus, we rewrite problem (18) in terms of the new operator A h , τ as:
Given θ h n , find θ h n + 1 V h such that : ( A h , τ θ h n + 1 , v h ) = ( l h dG , v h ) 0 v h V h ,
with
l h dG : = θ h n + τ h n + 1 .
We endow V h with the time-step dependent norm:
w h τ 2 : = w h 0 2 + τ w h V 2
and its extension:
w h τ , * 2 : = w h 0 2 + τ w h V , * 2 .
The operator A h , τ satisfies the inf-sup condition in terms of the norm w h τ 2 by extending Theorem 1. Moreover, the operator A h , τ is bounded in terms of the above norm and its extension w h τ , * 2 as:
( A h , τ ( u ) , v h ) u τ , * v h τ , ( v , w h ) V * h × V h .

3.4. Second-Order Backward Differencing Formula (BDF2)

As above, we use the second-order backward differencing formula as a time marching method to obtain a fully discrete solution,
δ t ( 2 ) v n + 1 : = 3 v n + 1 4 v n + v n 1 2 τ V n 1 , N .
For n = 1 , , k 1 , find θ h n + 1 V h , such that
( δ t ( 2 ) θ h n + 1 , v h n + 1 ) 0 + ( A h θ h n + 1 , v h n + 1 ) 0 = ( h n + 1 , v h n + 1 ) 0 v h n + 1 V h ,
for this case, we redefine the discrete-time operator A h , τ , as well as the billinear form a h , τ as: A h , τ : V * h V h , such that, for all ( u , w h ) V h * × V h ,
( A h , τ u , w h ) = a h , τ ( u , w h ) : = ( u , w h ) 0 + 3 2 τ a h ( u , w h ) .
We now write problem (26) following the derivation of (20) with
l h dG : = 2 3 τ h n + 1 + 4 3 θ h n 1 3 θ h n 1
and a given initial condition θ h 0 = π h u 0 . We compute θ h 1 , if necessary, with a first-order method. This operator satisfies the stability properties described in Section 3.3. Thus, updating (20) with the definitions (27) and (28), the operator A h , τ is well-posed.

4. Fully Discrete Residual Minimization

This section describes the stabilized finite element formulation via residual minimization on dual discontinuous Galerkin norms formulated in [30] for unsteady problems. We develop a method that delivers a stabilized discrete solution in a continuous space by minimizing the residual in a dual discontinuous norm at each time step. Thus, we choose V h as a broken polynomial space described in (3) and U h as its H 1 -conforming subspace. Following the formulation (20) in V h , we chose a trial conforming subspace U h V h to solve the following residual minimization problem:
Given u h n , find u h n + 1 U h V h , such that : u h n + 1 = arg min z h U h 1 2 l h A h , τ z h τ , * 2 = arg min z h U h 1 2 R τ 1 ( l h A h , τ z h ) τ 2 ,
where u h 0 = π h u 0 and l h is defined as l h : = u h n + τ h n + 1 for BDF1 and l h : = 2 3 τ h n + 1 + 4 3 u h n 1 3 u h n 1 for BDF2. R τ 1 denotes the inverse of the Riesz map:
R τ : V h V h * ( R τ y h , v h ) V h * × V h : = ( y h , v h ) τ v h V h .
Problem (29) is equivalent to the following saddle-point problem:
Given u h n , find ( ε h n + 1 , u h n + 1 ) V h × U h , such that : ( ε h n + 1 , v h ) τ + ( A h , τ u h n + 1 , v h ) = ( l h , v h ) 0 v h V h , ( A h , τ z h , ε h n + 1 ) = 0 , z h U h ,
where the residual representation function ε h n + 1 is defined by:
ε h n + 1 : = R τ 1 ( l h A h , τ ) V h ,
We write (31) in the dual space,
Given u h n , find ( ε h n + 1 , u h n + 1 ) V h × U h , such that : R τ ε h n + 1 + A h , τ u h n + 1 = l h , in V h * , A h , τ ε h n + 1 = 0 , in U h * .
Remark 1. 
Substituting the source term ( h n + 1 ) from (20) into the first identity in (33), we obtain, for BDF1, that:
R τ ε h n + 1 + A h , τ u h n + 1 = u h n + A h , τ θ h n + 1 θ h n .
Rearranging and defining the spatial error at time step i by ξ i : = θ h i u h i ; then, (34) implies that:
ε h n + 1 = R τ 1 ( A h , τ ξ n + 1 ξ n ) .
Or, equivalently, for the BDF2 implementation:
ε h n + 1 = R τ 1 ( A h , τ ξ n + 1 4 3 ξ n + 1 3 ξ n 1 ) .
Hence, we can alternatively define ε h n + 1 as an error measure distance from the continuous to discontinuous approximation at the n + 1 time step with the k previous time-step spatial error contributions (for a k-order BDF method).

Adaptive Mesh Refinement

This section describes the adaptive refinement procedure with the following steps. First, we solve the saddle-point problem (33) to obtain an error representation ( ε n + 1 V h ) in the norm ( ε n + 1 τ 2 ). Then, we construct a local version of the time-dependent norm (22), having an error indicator per cell ( E K ), such that:
E K 2 = ε n + 1 0 , K 2 + τ ε n + 1 V , L o c , 2
where
ε n + 1 V , L o c 2 : = κ ε n + 1 0 , K 2 + β c 1 h K β · ε n + 1 0 , K 2 + F S h n e κ + 1 2 | β · n F | ( ε n + 1 , ε n + 1 ) 0 , F .
Here, we use an extension of the Dörfler bulk-chasing criterion [43] to mark the cells with the highest E K values based on an accumulative error in a cell loop. We first organize the cells in the order of decreasing error per cell. Then, the algorithm marks the elements in two cases: when the accumulative error in a first loop reaches a user-defined fraction of the error ε n + 1 τ 2 , and when the error of the remaining cells in the first loop is larger than a chosen fraction of the last refined element. By refining all elements with comparable errors in an iteration, we guarantee refinement in the elements close to the cutoff, which the original strategy did not mark; this combined strategy reduces the computational cost per iteration. Let η r e f be 0.25 in 2D and 0.125 in 3D (see [33]) and ν = 0.2 in all cases. Then, we refine the marked cells using bisection. Algorithm 1 summarizes the marking strategy.
Algorithm 1 Marking strategy
Input: T h , ε n + 1 τ 2 , η r e f , N , ν
1:
Compute E K all K T h from (36)
2:
Sort and store in sortK all K T h from highest to lowest E K values
3:
Initialize cell to mark Kmarked = sortK[0]
4:
Initialize the local error of the marked mesh cell E K m = E K [ K m a r k e d ]
5:
Initialize sum = 0 , i = 0 , flag = True and E c u t = 0
6:
while (sum < η r e f 2 ε n + 1 τ 2 or E K m ( 1 ν ) E c u t ) and i < N do
7:
  Mark Kmarked
8:
  if sum < η r e f 2 ε τ 2 then
9:
    sum ← sum + E K m
10:
  else
11:
    if flag then
12:
       E c u t E K m
13:
      flag F a l s e
14:
  ii + 1
15:
  KmarkedsortK[i]
16:
   E K m E K [ K m a r k e d ]
The stopping criterion for the refinement algorithm is as follows. Starting with a coarse mesh, we refine while the total estimated error in the norm ε n + 1 τ remains above a time-step dependent tolerance E t o l = τ C t o l , where C t o l is a user-defined constant. For the numerical examples in this paper, we use C t o l = 1 × 10 5 . Algorithm 2 details the implementation of BDF1.
Algorithm 2 Algorithm for BDF1
Input: T h 0 u h 0 , T , τ , E t o l , n = 0
1:
u h n = u h 0
2:
for t { 0 T } do
3:
  while E E t o l do
4:
    Project u h n on T h n
5:
    Solve u h n + 1 and ε h n + 1 using the saddle point problem from (33)
6:
    Compute E using ε h n + 1
7:
    if E < E T o l then
8:
      Use ε h n + 1 and the marking criteria described in the Algorithm 1 to obtain the refined mesh T
9:
       T h n T
10:
   t t + τ
11:
   n n + 1

5. Numerical Examples

This section presents four numerical examples to show the performance properties of our adaptive stabilized finite element method. First, we solve the heat equation problem, obtaining optimal space and time convergences for uniform refinements. We use the classic Eriksson–Johnson problem in the second case to test the adaptive refinement strategy and its convergence in space for different polynomial degrees and Péclet numbers. The third example shows the stability in two dimensions for the unsteady pure-advection problem where the mesh moves in time. Here, we compare the computational time of the stabilized finite element method using adaptivity for a uniform mesh with the dG method. Finally, we show the performance of our procedure in a nonlinear unsteady reaction–diffusion problem with two-branched numerical solutions.
Since we minimized the residual in the energy norm ( τ ), we focused this research on the spatial convergence study in this norm. For the following numerical examples, we implement the iterative algorithm described in [30,44] to solve the resulting saddle point system (33) and use FEniCS [45] as a platform to perform all the numerical simulations.

5.1. Heat Equation (2D)

We start the method’s performance analysis by solving the 2D heat equation while refining the spatial domain uniformly. Although this case does not present any particular challenge to classical methods, it is a standard benchmark problem to test space and time convergences of parabolic problems.
Let the domain Ω be [ 0 , 1 ] 2 ; we consider the problem:
t u Δ u = f in Ω × ( 0 , T ] , u = 0 on Γ D × ( 0 , T ] , u ( · , t = 0 ) = u 0 in Ω ,
with the initial condition:
u 0 = sin ( π x ) sin ( π y ) ,
and the source term f that satisfies the exact solution:
u ( ( x , y ) , t ) = exp ( π 2 t ) sin ( π x ) sin ( π y ) .
To formulate the fully discrete problem, we combine the SIP billinear form (for κ = 1 ) with the BDF1/BDF2 time marching scheme. We show the convergence plots for linear and quadratic polynomials in space (Figure 2 and Figure 3) and time (Figure 4). To perform the spatial convergence test, we compute the errors u u h τ (in black), u θ h τ (in blue), θ h u h τ (in green) and ε h τ (in red) for different mesh sizes ( Δ x ). We denote Δ x equal to h K for uniform meshes. To perform the time convergence, we show u u h 0 varying with the time step τ for Δ x = 0.01 . As a result, we recover space optimality for the continuous approximation (from dG formulation) and the first- and second-order time convergence for BDF1 and BDF2, respectively. We show that the saturation assumption in [30] holds in our formulation (i.e., u ( T ) θ h ( T ) τ u ( T ) u h ( T ) τ ) and the residual representation is efficient until the error dominates (i.e., ε h τ u ( T ) u h ( T ) τ in [30]). To illustrate, Figure 3b shows that for p = 2 and high DoF, the temporal error is no longer negligible to the spatial error; however, the error estimator continues decaying since it does not consider the temporal error contribution. We will seek to prove these properties in future work.

5.2. Advection–Diffusion Problem

We complement the spatial convergence study of the previous example, using adaptive refinement for the unsteady advection–dominated Eriksson–Johnson problem. Let the domain Ω be [ 0 , 1 ] × [ 0.5 , 0.5 ] ; we consider the exact solution
u ( ( x , y ) , t ) = exp ( l t ) [ exp ( λ 1 x ) exp ( λ 2 x ) ] + c o s ( π y ) exp ( s 1 x ) exp ( r 1 x ) exp ( s 1 ) exp ( r 1 ) ,
for f = 0 and l = 2 , λ 1 , 2 = 1 ± 1 4 κ l 2 κ , r 1 = 1 + 1 + 4 κ 2 π 2 2 κ and s 1 = 1 1 + 4 κ 2 π 2 2 κ . Here, we set β = [ 1 , 0 ] and μ = 0 for different diffusion coefficient values. Based on the exact solution, we apply Neumann boundary conditions at x = 1 and t = 0 ; meanwhile, we impose Dirichlet boundary conditions at x = 0 , y = 0.5 and y = 0.5 at time t = 0 .
The problem’s main challenge is capturing the boundary layer, especially for high Péclet numbers. Figure 5 shows how the error estimator drives spatial adaptivity to smooth the regions with sharp gradients in each time step. Figure 6 shows the errors u u h τ , u θ h τ and ε h τ versus the square of total degrees of freedom ( D o F 1 / 2 ) (i.e., dim ( U h ) + dim ( V h ) ); these plots verify the optimal spatial convergence in the fully discrete energy norm using BDF1 and BDF2 time integrators for linear and quadratic polynomial trial functions at T = 0.1 and τ = 0.005 . In Figure 7, we verify our method’s convergence for higher Péclet numbers by setting the diffusivity in 10 3 and 10 4 . Figure 8 shows the evolution of our transient solution to the analytical steady-state Eriksson–Johnson. Similarly to the uniform refinement case, we preserve the efficiency of the residual representative and the saturation assumptions as stated for the adaptive steady state case.

5.3. Rotating Flow Transporting a Gaussian Profile

In this example, we analyze the performance of our method in a convective transport problem with a localized disturbance. Thus, we test the adaptive algorithm in the case where the region of interest moves within the domain as time passes. We study the solution of a 2D convection–diffusion transport by a rotatory flow of a Gaussian profile. We set Ω = [ 2 , 2 ] 2 , T = π , β = [ y , x ] , κ = 10 5 , μ = 0 and f = 0 . The initial condition is
u 0 = exp ( 64 ( x 0.5 ) 2 ) exp ( 64 y 2 ) ,
we impose Dirichlet boundary conditions from the exact solution:
u ( ( x , y ) , t ) = 1 1 + 256 κ t exp 64 ( x 0.5 cos ( t ) ) 2 1 + 256 κ t exp 64 ( y + 0.5 sin ( t ) ) 2 1 + 256 κ t .
Figure 9 shows the profiles of the solutions and the corresponding adaptively refined meshes at different time steps. These results demonstrate the continuous solution stability and consistency with the physical phenomena, even for low diffusivities. Moreover, the mesh nodes concentrate where the solution varies largely, showing the robustness of the error estimator and the efficiency of the marking strategy when adding new degrees of freedom. Regarding computational cost, our stabilized finite element formulation using adaptivity is competitive with the dG methodology using uniform refinement. Solving the saddle point formulation requires an extra cost due to the additional degrees of freedom. However, adaptivity compensates for the excess due to the solution’s stability in coarse meshes and the robustness of the error estimator. Figure 10 shows a comparison between the total computational cost required to obtain a solution with the adaptive stabilized method (blue line) and the computational cost using a regular mesh in the dG method (red line). Besides, the figure shows that the adaptivity can reduce the computational cost by up to one order of magnitude to get a resolution of 1 × 10 5 in the energy norm.

5.4. Unsteady Bratu Equation: Non-Linear Diffusion-Reaction Equation

We conclude the method’s performance analysis by studying the solution of a nonlinear diffusion-reaction equation. Let λ be a positive real constant and Ω = [ 0 , 1 ] 2 ; we solve the unsteady version of Bratu’s problem in the following form:
Find u such that , for T > 0 , t u = Δ u + λ exp ( u ) in Ω × ( 0 , T ] , u = 0 on Γ × ( 0 , T ] , u ( · , t = 0 ) = u 0 ( x ) in Ω ,
The two-dimensional steady version of (39) has a branched solution for λ < λ c (lower and upper branches) and a unique solution when λ = λ c , with λ c 6.8081 as a critical point. The problem’s main challenges are the lack of stable solutions in the upper branch and close to the critical point λ c , leading to classical techniques converging only to the stable lower branch. We test our method’s robustness, accuracy and performance in this transient bifurcation problem; we compare the solutions obtained in (39) when t , with the 2D steady Bratu’s approach obtained in [33]. We formulate the space semi-discretization of (39) as follows:
Find θ h V h , such that : ( t θ h , v h ) 0 + η h ( θ h ; v h ) = h ( v h ) , v h V h ,
where η h ( u h ; v h ) denotes the nonlinear form, including the SIP formulation in (7) with a non-linear reactive contribution. We define it as:
η h ( u h ; v h ) : = K T ( h u h · h v h ) K K S h ( λ exp ( u h ) , v h ) K F S h ( h u h · n F , [ [ v h ] ] ) F + ( [ [ u h ] ] , h v h · n F ) F n e κ ( [ [ u h ] ] , [ [ v h ] ] ) F .
For BDF1 time marching, we define the discrete-time nonlinear form as:
η h , τ ( u h ; v h ) : = ( u h , v h ) K + τ η h ( u h ; v h ) .
Thus, the full-discrete formulation for problem (39) is:
Given θ h n , find θ h n + 1 V h such that : η h , τ ( θ h n + 1 ; v h ) = ( l h dG , v h ) 0 v h V h ,
with l h dG : = θ h n + τ h n + 1 . We use a Newton–Raphson iteration scheme combined with the residual minimization strategy described in Section 4 to solve (43). We seek a solution at every Newton step increment w h for each time step n + 1 by using the linearized form:
η h , τ ( u h ; w h , v h ) : = ( w h , v h ) K + τ ( K T ( h w h · h v h ) K K S h ( λ exp ( u h ) w h , v h ) K F S h ( h w h · n F , [ [ v h ] ] ) F + ( [ [ w h ] ] , h v h · n F ) F n e κ ( [ [ w h ] ] , [ [ v h ] ] ) F )
Since (44) takes the form of a diffusion-reaction problem, we use a time-step dependent norm (22), with the SIP contribution to the V h norm, to minimize the discrete residual of the linearized system. The norm · τ is enforced with an L 2 contribution to measure the nonlinear reactive term. Starting with an initial guess ( ε h , 0 n + 1 , u h , 0 n + 1 ) and given ( ε h , i n + 1 , u h , i n + 1 ) , we find:
( δ ε h n + 1 , δ u h n + 1 ) V h × U h , such that : ( z h , v h ) V h , × U h ( δ ε h n + 1 , v h ) τ + η h , τ ( u h , i n + 1 ; δ u h n + 1 , v h ) = ( l h , v h ) 0 ( ε h , i n + 1 , v h ) τ η h , τ ( u h , i n + 1 ; v h ) η h , τ ( u h , i n + 1 ; z h , δ ε h n + 1 ) = η h , τ ( u h , i n + 1 ; z h , ε h , i n + 1 )
u h , i n + 1 and ε h , i n + 1 are updated at every i-th increment as follows:
u h , i + 1 n + 1 = u h , i n + 1 + k δ u h n + 1 , ε h , i + 1 n + 1 = ε h , i n + 1 + k δ ε h n + 1 ,
where k denotes a relaxation parameter from the Damped Newton’s method [46], and it is detailed to our formulation’s context in [33]. For the time step n = 1 , we set the initial guess ( ε h , 0 n + 1 , u h , 0 n + 1 ) = ( 0 , u IG ) , where u IG varies depending on the solution branch we want to capture. Here, we assume u IG equal to the initial solution (i.e., u IG = u 0 ) with u 0 = 0 for the stable lower branch and u 0 = u up for the upper branch. Since the lower branch is stable, many different initial guesses converge to it; however, we only use one option. The unstable upper branch is more restrictive; therefore, we follow [47] and use:
u up ( x , y ) = 50 ( 2 + λ ) λ ( x x 2 ) ( y y 2 ) .
Figure 11 and Figure 12 show the two branch solutions obtained for a time step increment τ = 0.1 with an initial mesh of 4 × 4 elements and a final time T = 1.0 . Figure 11 shows a time sequence for lower and upper solutions from t = 0 to t = T at λ = 2 . Figure 12 shows the classical bifurcation diagram for Bratu’s problem evaluating the maximum value u max at the time T for different λ values from 0 to λ c . At this time, we guarantee a convergent solution over time to approach the steady state of this problem. We demonstrate the robustness of our approach and the efficient refinement strategy to capture all stable and unstable branches, even close to the critical point ( λ c ). We test the accuracy of the results by successfully validating our bifurcation map at time T with results obtained from different authors at the steady state [31,47].

6. Discussion

This paper proposes an adaptive-stabilized finite element method based on residual minimization for unsteady advection–diffusion–reaction problems using the method of lines. The method provides a stable solution and a robust error representation to guide adaptivity at every discrete time. We demonstrate the method’s performance for challenging linear and nonlinear problems with optimal spatial and temporal convergence and the efficiency of the adaptive refinement strategy to capture sharp inner and boundary layers. We present evidence that the adaptive refinement process could overcome the required computational cost to solve the saddle-point problem compared to the uniformly refined schemas on discontinuous Galerkin approximation. Other time-marching approaches, including explicit time-marching schemas with time adaptivity and space–time formulation, will be described in future publications.

Author Contributions

Conceptualization, J.F.G. and V.M.C.; methodology, J.F.G. and V.M.C.; software, J.F.G.; validation, J.F.G. and V.M.C.; formal analysis, J.F.G. and V.M.C.; investigation, J.F.G. and V.M.C.; resources, J.F.G. and V.M.C.; writing—original draft preparation, J.F.G. and V.M.C.; writing—review and editing, V.M.C.; supervision, V.M.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

This work is supported by the close collaboration between Curtin University and CSIRO under the CSIRO DEI FSP postgraduate top-up scholarship (Grant no. 50068868). J.G. gratefully acknowledges Roberto Rocca Education Program for its support.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Brooks, A.N.; Hughes, T.J. Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations. Comput. Methods Appl. Mech. Eng. 1982, 32, 199–259. [Google Scholar] [CrossRef]
  2. Hughes, T.J.R.; Franca, L.P.; Hulbert, G.M. A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/Least-Squares method for advection–diffusive equations. Comput. Methods Appl. Mech. Eng. 1989, 73, 173–189. [Google Scholar] [CrossRef]
  3. Hughes, T.J.; Feijóo, G.R.; Mazzei, L.; Quincy, J.B. The variational multiscale method—A paradigm for computational mechanics. Comput. Methods Appl. Mech. Eng. 1998, 166, 3–24. [Google Scholar] [CrossRef]
  4. Hughes, T.J.R.; Scovazzi, G.; Franca, L.P. Multiscale and Stabilized Methods. In Encyclopedia of Computational Mechanics, 2nd ed.; American Cancer Society: Kennesaw, GA, USA, 2017; pp. 1–64. [Google Scholar]
  5. Bazilevs, Y.; Calo, V.; Cottrell, J.; Hughes, T.; Reali, A.; Scovazzi, G. Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows. Comput. Methods Appl. Mech. Eng. 2007, 197, 173–201. [Google Scholar] [CrossRef]
  6. Bazilevs, Y.; Michler, C.; Calo, V.; Hughes, T. Isogeometric variational multiscale modeling of wall-bounded turbulent flows with weakly enforced boundary conditions on unstretched meshes. Comput. Methods Appl. Mech. Eng. 2010, 199, 780–790. [Google Scholar] [CrossRef]
  7. Chang, K.; Hughes, T.; Calo, V. Isogeometric variational multiscale large-eddy simulation of fully-developed turbulent flow over a wavy wall. Comput. Fluids 2012, 68, 94–104. [Google Scholar] [CrossRef]
  8. Ghaffari Motlagh, Y.; Ahn, H.T.; Hughes, T.J.; Calo, V.M. Simulation of laminar and turbulent concentric pipe flows with the isogeometric variational multiscale method. Comput. Fluids 2013, 71, 146–155. [Google Scholar] [CrossRef]
  9. Arnold, D.N.; Brezzi, F.; Cockburn, B.; Donatella Marini, L. Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM J. Numer. Anal. 2001, 39, 1749–1779. [Google Scholar] [CrossRef]
  10. Burman, E.; Hansbo, P. Edge stabilization for Galerkin approximations of convection-diffusion-reaction problems. Comput. Methods Appl. Mech. Eng. 2004, 193, 1437–1453. [Google Scholar] [CrossRef] [Green Version]
  11. Brezzi, F.; Marini, L.D.; Süli, E. Discontinuous Galerkin methods for first-order hyperbolic problems. Math. Model. Methods Appl. Sci. 2004, 14, 1893–1903. [Google Scholar] [CrossRef]
  12. Ayuso, B.; Marini, L.D. Discontinuous Galerkin methods for advection–diffusion–reaction problems. SIAM J. Numer. Anal. 2009, 47, 1391–1420. [Google Scholar] [CrossRef]
  13. Cockburn, B.; Karniadakis, G.E.; Shu, C.W. Discontinuous Galerkin Methods: Theory, Computation and Applications; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2012; Volume 11. [Google Scholar]
  14. Demwkowicz, L.; Gopalakrishnan, J. Analysis of the DPG Method for the Poisson Equation. SIAM J. Numer. Anal. 2011, 49, 1788–1809. [Google Scholar] [CrossRef] [Green Version]
  15. Zitelli, J.; Muga, I.; Demkowicz, L.; Gopalakrishnan, J.; Pardo, D.; Calo, V. A class of discontinuous Petrov–Galerkin methods. Part IV: The optimal test norm and time-harmonic wave propagation in 1D. J. Comput. Phys. 2011, 230, 2406–2432. [Google Scholar] [CrossRef] [Green Version]
  16. Niemi, A.H.; Collier, N.O.; Calo, V.M. Automatically stable discontinuous Petrov–Galerkin methods for stationary transport problems: Quasi-optimal test space norm. Comput. Math. Appl. 2013, 66, 2096–2113. [Google Scholar] [CrossRef]
  17. Cockburn, B.; Shu, C.W. The local discontinuous Galerkin method for time-dependent convection–diffusion systems. SIAM J. Numer. Anal. 1998, 35, 2440–2463. [Google Scholar] [CrossRef]
  18. Borker, R.; Farhat, C.; Tezaur, R. A high-order discontinuous Galerkin method for unsteady advection–diffusion problems. J. Comput. Phys. 2017, 332, 520–537. [Google Scholar] [CrossRef]
  19. Rivière, B. Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations: Theory and Implementation; SIAM: Philadelphia, PA, USA, 2008. [Google Scholar]
  20. Di Pietro, D.A.; Ern, A. Mathematical Aspects of Discontinuous Galerkin Methods; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2011; Volume 69. [Google Scholar]
  21. Demkowicz, L.; Gopalakrishnan, J. A class of discontinuous Petrov–Galerkin methods. Part I: The transport equation. Comput. Methods Appl. Mech. Eng. 2010, 199, 1558–1572. [Google Scholar] [CrossRef] [Green Version]
  22. Führer, T.; Heuer, N.; Karkulik, M. Analysis of backward Euler primal DPG methods. Comput. Methods Appl. Math. 2021, 21, 811–826. [Google Scholar] [CrossRef]
  23. Barrenechea, G.R.; Brezzi, F.; Cangiani, A.; Georgoulis, E.H. Building Bridges: Connections and Challenges in Modern Approaches to Numerical Partial Differential Equations; Springer: Berlin/Heidelberg, Germany, 2016. [Google Scholar]
  24. Roberts, N.V.; Henneking, S. Time-stepping DPG formulations for the heat equation. Comput. Math. Appl. 2020, 95, 242–255. [Google Scholar] [CrossRef]
  25. Ern, A.; Vohralík, M. A posteriori error estimation based on potential and flux reconstruction for the heat equation. SIAM J. Numer. Anal. 2010, 48, 198–223. [Google Scholar] [CrossRef]
  26. Zhu, L.; Schötzau, D. A robust a posteriori error estimate for hp-adaptive DG methods for convection–diffusion equations. IMA J. Numer. Anal. 2010, 31, 971–1005. [Google Scholar] [CrossRef] [Green Version]
  27. Ern, A.; Proft, J. A posteriori discontinuous Galerkin error estimates for transient convection-diffusion equations. Appl. Math. Lett. 2005, 18, 833–841. [Google Scholar] [CrossRef] [Green Version]
  28. Araya, R.; Venegas, P. An a posteriori error estimator for an unsteady advection-diffusion- reaction problem. Comput. Math. Appl. 2014, 66, 2456–2476. [Google Scholar] [CrossRef]
  29. Cangiani, A.; Georgoulis, E.H.; Metcalfe, S. Adaptive discontinuous Galerkin methods for nonstationary convection-diffusion problems. IMA J. Numer. Anal. 2014, 34, 1578–1597. [Google Scholar] [CrossRef] [Green Version]
  30. Calo, V.M.; Ern, A.; Muga, I.; Rojas, S. An adaptive stabilized conforming finite element method via residual minimization on dual discontinuous Galerkin norms. Comput. Methods Appl. Mech. Eng. 2020, 363, 112891. [Google Scholar] [CrossRef] [Green Version]
  31. Cier, R.J.; Rojas, S.; Calo, V.M. Automatically adaptive, stabilized finite element method via residual minimization for heterogeneous, anisotropic advection–diffusion–reaction problems. Comput. Methods Appl. Mech. Eng. 2021, 385, 114027. [Google Scholar] [CrossRef]
  32. Calo, V.; Łoś, M.; Deng, Q.; Muga, I.; Paszyński, M. Isogeometric Residual Minimization Method (iGRM) with direction splitting preconditioner for stationary advection-dominated diffusion problems. Comput. Methods Appl. Mech. Eng. 2021, 373, 113214. [Google Scholar] [CrossRef]
  33. Cier, R.J.; Rojas, S.; Calo, V.M. A nonlinear weak constraint enforcement method for advection-dominated diffusion problems. Mech. Res. Commun. 2021, 112, 103602. [Google Scholar] [CrossRef]
  34. Rojas, S.; Pardo, D.; Behnoudfar, P.; Calo, V.M. Goal-oriented adaptivity for a conforming residual minimization method in a dual discontinuous Galerkin norm. Comput. Methods Appl. Mech. Eng. 2021, 377, 113686. [Google Scholar] [CrossRef]
  35. Kyburg, F.E.; Rojas, S.; Calo, V.M. Incompressible flow modeling using an adaptive stabilized finite element method based on residual minimization. Int. J. Numer. Methods Eng. 2022, 123, 1717–1735. [Google Scholar] [CrossRef]
  36. Łoś, M.; Rojas, S.; Paszyński, M.; Muga, I.; Calo, V.M. DGIRM: Discontinuous Galerkin based isogeometric residual minimization for the Stokes problem. J. Comput. Sci. 2021, 50, 101306. [Google Scholar] [CrossRef]
  37. Poulet, T.; Giraldo, J.F.; Ramanaidou, E.; Piechocka, A.; Calo, V.M. Paleo-stratigraphic permeability anisotropy controls supergene mimetic martite goethite deposits. Basin Res. 2022. [Google Scholar] [CrossRef]
  38. Labanda, N.A.; Espath, L.; Calo, V.M. A spatio-temporal adaptive phase-field fracture method. Comput. Methods Appl. Mech. Eng. 2022, 392, 114675. [Google Scholar] [CrossRef]
  39. Bramble, J.H.; Thomée, V. Semidiscrete least-squares methods for a parabolic boundary value problem. Math. Comput. 1972, 26, 633–648. [Google Scholar]
  40. Becker, J. A second order backward difference method with variable steps for a parabolic problem. BIT Numer. Math. 1998, 38, 644–662. [Google Scholar] [CrossRef]
  41. Crouzeix, M.; Lisbona, F. The convergence of variable-stepsize, variable-formula, multistep methods. SIAM J. Numer. Anal. 1984, 21, 512–534. [Google Scholar] [CrossRef]
  42. Ern, A.; Stephansen, A.F.; Zunino, P. A discontinuous Galerkin method with weighted averages for advection–diffusion equations with locally small and anisotropic diffusivity. IMA J. Numer. Anal. 2009, 29, 235–256. [Google Scholar] [CrossRef]
  43. Dörfler, W. A convergent adaptive algorithm for Poisson’s equation. SIAM J. Numer. Anal. 1996, 33, 1106–1124. [Google Scholar] [CrossRef]
  44. Bank, R.E.; Welfert, B.D.; Yserentant, H. A class of iterative methods for solving saddle point problems. Numer. Math. 1989, 56, 645–666. [Google Scholar] [CrossRef]
  45. Alnæs, M.; Blechta, J.; Hake, J.; Johansson, A.; Kehlet, B.; Logg, A.; Richardson, C.; Ring, J.; Rognes, M.E.; Wells, G.N. The FEniCS project version 1.5. Arch. Numer. Softw. 2015, 3. [Google Scholar] [CrossRef]
  46. Bank, R.E.; Rose, D.J. Global approximate Newton methods. Numer. Math. 1981, 37, 279–295. [Google Scholar] [CrossRef]
  47. Hajipour, M.; Jajarmi, A.; Baleanu, D. On the accurate discretization of a highly nonlinear boundary value problem. Numer. Algorithms 2018, 79, 679–695. [Google Scholar] [CrossRef]
Figure 1. Notation of the element interface.
Figure 1. Notation of the element interface.
Mca 28 00007 g001
Figure 2. BDF1 spatial convergence using fixed time step and uniform meshes.
Figure 2. BDF1 spatial convergence using fixed time step and uniform meshes.
Mca 28 00007 g002
Figure 3. BDF2 spatial convergence using fixed time step and uniform meshes.
Figure 3. BDF2 spatial convergence using fixed time step and uniform meshes.
Mca 28 00007 g003aMca 28 00007 g003b
Figure 4. Time convergence for BDF1 and BDF2 time integrators using a fixed mesh.
Figure 4. Time convergence for BDF1 and BDF2 time integrators using a fixed mesh.
Mca 28 00007 g004
Figure 5. Mesh refinement τ = 0.005 , T = 0.1 , p = 1 .
Figure 5. Mesh refinement τ = 0.005 , T = 0.1 , p = 1 .
Mca 28 00007 g005
Figure 6. Spatial convergence for adaptive refinement (BDF1 & BDF2: κ = 10 2 , T = 0.1 , τ = 0.005 ).
Figure 6. Spatial convergence for adaptive refinement (BDF1 & BDF2: κ = 10 2 , T = 0.1 , τ = 0.005 ).
Mca 28 00007 g006
Figure 7. Spatial convergence for adaptive refinement using BDF2. p = 2 , T = 0.1 , τ = 0.005 .
Figure 7. Spatial convergence for adaptive refinement using BDF2. p = 2 , T = 0.1 , τ = 0.005 .
Mca 28 00007 g007
Figure 8. Solution convergence to the steady Eriksson–Johnson solution.
Figure 8. Solution convergence to the steady Eriksson–Johnson solution.
Mca 28 00007 g008
Figure 9. Time evolution p = 1 ( T = π , τ = π / 512 ).
Figure 9. Time evolution p = 1 ( T = π , τ = π / 512 ).
Mca 28 00007 g009
Figure 10. Computational cost [s] vs. total degrees of freedom.
Figure 10. Computational cost [s] vs. total degrees of freedom.
Mca 28 00007 g010
Figure 11. Solution’s temporal evolution for λ = 2 for the lower and upper branches.
Figure 11. Solution’s temporal evolution for λ = 2 for the lower and upper branches.
Mca 28 00007 g011
Figure 12. Bratu’s bifurcation diagram for T = 1.0 and τ = 0.1 .
Figure 12. Bratu’s bifurcation diagram for T = 1.0 and τ = 0.1 .
Mca 28 00007 g012
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Giraldo, J.F.; Calo, V.M. An Adaptive in Space, Stabilized Finite Element Method via Residual Minimization for Linear and Nonlinear Unsteady Advection–Diffusion–Reaction Equations. Math. Comput. Appl. 2023, 28, 7. https://doi.org/10.3390/mca28010007

AMA Style

Giraldo JF, Calo VM. An Adaptive in Space, Stabilized Finite Element Method via Residual Minimization for Linear and Nonlinear Unsteady Advection–Diffusion–Reaction Equations. Mathematical and Computational Applications. 2023; 28(1):7. https://doi.org/10.3390/mca28010007

Chicago/Turabian Style

Giraldo, Juan F., and Victor M. Calo. 2023. "An Adaptive in Space, Stabilized Finite Element Method via Residual Minimization for Linear and Nonlinear Unsteady Advection–Diffusion–Reaction Equations" Mathematical and Computational Applications 28, no. 1: 7. https://doi.org/10.3390/mca28010007

APA Style

Giraldo, J. F., & Calo, V. M. (2023). An Adaptive in Space, Stabilized Finite Element Method via Residual Minimization for Linear and Nonlinear Unsteady Advection–Diffusion–Reaction Equations. Mathematical and Computational Applications, 28(1), 7. https://doi.org/10.3390/mca28010007

Article Metrics

Back to TopTop