A POSTERIORI ERROR ESTIMATION AND ADAPTIVITY ... - CiteSeerX
A POSTERIORI ERROR ESTIMATION AND ADAPTIVITY FOR DEGENERATE PARABOLIC PROBLEMS R.H. NOCHETTO, A. SCHMIDT, and C. VERDI Abstract. Two explicit error representation formulas are derived for degenerate
parabolic PDEs, which are based on evaluating a parabolic residual in negative norms. The resulting upper bounds are valid for any numerical method, and rely on regularity properties of solutions of a dual parabolic problem in nondivergence form with vanishing di usion coecient. They are applied to a practical space-time discretization consisting of C 0 piecewise linear nite elements over highly graded unstructured meshes, and backward nite di erences with varying time-steps. Two rigorous a posteriori error estimates are derived for this scheme, and used in designing an ecient adaptive algorithm, which equidistributes space and time discretization errors via re nement/coarsening. A simulation nally compares the behavior of the rigorous a posteriori error estimators with a heuristic approach, and hints at the potentials and reliability of the proposed method.
1. Introduction
A posteriori error estimates are a fundamental component in the design of reliable and ecient adaptive algorithms for the numerical solution of PDEs. Even though rigorous results are available for linear and mildly nonlinear parabolic PDEs [5], [6], [7], the theory is much less satisfactory for strongly nonlinear PDEs. There are no results applicable to degenerate parabolic PDEs, which in turn exhibit lack of regularity across interfaces and corresponding numerical pollution e ects. The use of highly graded meshes and varying time-steps is thus motivated by the nonlinear structure of the PDE, as opposed to domain geometry, and is a vehicle for resolving small scale features with optimal computational complexity. In this paper we introduce a rigorous theory of a posteriori error estimation for degenerate parabolic problems of the form (1.1)
@t u ? (u) = f in Q := (0; T );
where is nondecreasing and Lipschitz. A typical example of industrial interest is that of solidi cation (classical Stefan problem), for which (1.2)
(s) := ? min(s; 0) + + max(s ? L; 0):
1991 Mathematics Subject Classi cation. 65N15, 65N30, 65N50, 80A22, 35K65, 35R35. Key words and phrases. degenerate parabolic equations, Stefan problem, nite elements, parabolic duality, a posteriori estimates, adaptivity. This work was partially supported by NSF Grants DMS-9305935 and DMS-9623394, EU Grant HCM \Phase Transitions and Surface Tension", MURST, CNR Contract 95.00735.01 and TMR contract \Viscosity Solutions and their Applications". Typeset by AMS-TEX 1
2
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
This corresponds to an ideal material with constant thermal coecients ? ; + > 0 and latent heat L. Any approximation U of u satis es (1.3) @t U ? (U ) = f ? R in Q; where R, an oscillatory distribution of singular character, is the so-called parabolic residual. In spite of the simple parabolic structure of (1.1), (1.2), its degenerate nature makes the theory of [5] fail in that it exploits the regularizing e ect of a linear parabolic dual problem. The corresponding dual PDE in this context is the nonstrictly parabolic equation in nondivergence form (1.4) @t + b = in Q; with vanishing and rough di usion coecient 0 b max( ? ; +) [8], [9], [10]. Such a PDE does not exhibit a regularizing mechanism, and is not computable in that b is discontinuous and depends on both u and U . Problem (1.4) measures the error accumulation in time, and is thus crucial. Our objective is to prove global regularity properties of in x2, and use them in x3 to derive two representation formulas for the errors u ? U and (u) ? (U ) in energy norms. These formulas are valid for any numerical method, evaluate R in two distinct negative norms, and lead to rigorous a posteriori upper error bounds. Since negative norms entail averaging, they are appropriate to quantify the oscillatory character of R. We next apply these ideas to a practical scheme consisting of C 0 piecewise linear nite elements over highly graded unstructured meshes and backward nite di erences with varying time-steps. The method uses mass lumping and evaluates (U ) solely at the nodes, which makes it easy to implement and solve iteratively. We discuss the method in x4 and derive in x5 two rigorous a posteriori error estimates for it of the form (Approaches I and II): (1.5) ku ? U kL2 (Q) + k (u) ? (U )kL1 (0;T ;H ?1( )) E (u0 ; f; T; ; U; h; ): The estimator E is computable in terms of data u0 = u(; 0); f; T; , computed solution U , meshsize h, and time step , but entails L1 or L2 norms in time. This is impractical in that the entire evolution history would be needed to control E . We thus resort to an L1 norm in time, and introduce an adaptive algorithm which equidistributes space discretization errors for a uniform error distribution in time. Such errors are estimated via local a posteriori error indicators, and further equidistributed via a re nement/coarsening strategy based on bisection. This yields compatible consecutive meshes. These issues are fully discussed in x6. We conclude in x7 with simulations illustrating the viability of our Approaches I and II, as well as a heuristic Approach III based on using local regularity of in (1.4) and heat estimators away from discrete interfaces. We clearly show that they are all able to detect the presence of interfaces, and re ne accordingly, and that Approaches II and III perform best. There is no need to compute the interface explicitly for mesh design, which is a major improvement with respect to [12]. Further simulations, comparisons of several nonlinear solvers, and a detailed description of the adaptive algorithm will be presented elsewhere [14]. Even though our error estimates (1.5) are rigorous, they do not necessarily imply E (u0 ; f; T; ; U; h; ) ! 0 as h; ! 0 because E depends on discrete quantities that change with h and . Stability and convergence will be assessed in [13].
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE
3
2. Setting and Parabolic Duality
Let Rd (d > 1) be a bounded convex polyhedral domain, let T > 0 be xed, and set Q := (0; T ). Let 2 W 1;1(R) satisfy (s) = 0 for all s 2 (0; L) and 0 < a 0(s) A for a.e. s 2= (0; L). Let u0 indicate the initial enthalpy, let 0 := (u0) denote the initial temperature, and let F0 := fx 2 : 0(x) = 0g be the initial interface. They satisfy 0 2 W01;1( ); F0 is a Lipschitz curve: Therefore, u0 2 W 1;1( nF0 ) and u0 has a jump discontinuity across F0. Finally let f be suciently smooth. The continuous problem then reads as follows. Continuous problem. Find u and such that 2 L2(0; T ; H01 ( )); u 2 L1(0; T ; L2 ( )) \ H 1(0; T ; H ?1 ( )); (x; t) = (u(x; t)) a.e. (x; t) 2 Q; u(; 0) = u0; and for a.e. t 2 (0; T ) and all 2 H01 ( ) the following equation holds (2.1) h@t u; i + hr; ri = hf; i: Hereafter, h; i stands for either the inner product in L2( ) or the duality pairing between H ?1 ( ) and H01( ). It is to be observed that the vanishing Dirichlet boundary condition on is assumed only for simplicity and so that the interface F (t) := fx 2 : (x; t) = 0g does not include @ . Existence and uniqueness for this problem are known [8], [9]. To motivate the dual problem (1.4), already studied in [8], [9], [10], we subtract (1.3) from (1.1) and integrate by parts over Q against a smooth test function vanishing on @ (0; T ). The error eu := u ? U satis es Z
T
heu ; ijt=T = heu ; ijt=0 + heu; @t + b i + R( ); 0 where 0 b(x; t) A is the discontinuous function ( (u(x;t))? (U (x;t)) if u(x; t) = 6 U (x; t) u(x;t)?U (x;t) b(x; t) :=
(2.2)
A otherwise: We could thus represent norms of eu(; T ) or their integrals over Q by making judicious choices of (; T ) and @t + b . Evaluation of R( ) depends on regularity of , which we now investigate. Given a regularization parameter > 0 to be chosen later, we consider two backward parabolic problems (2.3) J ( ) = ?b1=2 in Q; (; T ) = 0 in ; (2.4) J ( ) = 0 in Q; (; T ) = in ; with operator J in nondivergence form (2.5) J ( ) := @t + (b + ); ; = 0 on @ (0; T ), and 2 L2(Q), 2 H01 ( ). Existence of unique solutions ; 2 H 1;2 (Q) follows directly from the theory of nonlinear strictly parabolic problems [9].
4
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
Lemma 2.1. The following a priori bounds are valid for all 0 t T (2.6) 2kr (; t)k2L2 ( ); A+1 k@t k2L2 (Q); 4k k2L2 (Q) k k2L2 (Q): Proof. We multiply (2.3) by and integrate by parts in (t; T ) to obtain Z
1 2
jr
(; t)j2 +
Z Z
t
T
bj
Z Z
T
(b + )j
j2 = ?
t Z Z T 2 1 j j2: j +4
t
Z Z
t
T
b1=2
This yields the a priori bounds for r and . Finally, we multiply (2.3) by @t =(b + ) and integrate by parts in (t; T ) to get Z Z
T
Z
Z Z
T
b1=2 @s 1 j@s j2 + 1 jr (; t)j2 = ? b + 2
t b+
t Z Z T Z Z T 2 1 1 1 j@s j + 2 j j2; 2 b +
t
t
because > 0. This implies the a priori bounds for @t because b A. Lemma 2.2. The following a priori bounds are valid for all 0 t T (2.7)
kr (; t)k2L2 ( ); A+2 k@t k2L2 (Q); 2k k2L2(Q) krk2L2( ):
Proof. We multiply (2.4) by and integrate by parts in (t; T ) to get Z
Z Z
1 jr (; t)j2 + 2
t
T
Z
(b + )j j2 = 21 jrj2:
The a priori bound for thus follows from b 0. In view of (2.4) we further have
k@t k2L2(Q) = k(b + ) k2L2(Q) (A + )k(b + )1=2 k2L2(Q) A2+ krk2L2( ); where we have used that b A. This completes the proof. Corollary 2.1. The following L2t Hx2 a priori bounds are valid (2.8)
Z
0
T
j
j2 2
H ( )
2 1 4 k kL2 (Q) ;
Z
0
T
j j2H 2( ) 21 krk2L2( ):
Proof. It suces to invoke the well-known estimate [5], [9] for convex
jjH 2( ) kkL2( ) 8 2 H01 ( ) \ H 2 ( ); in conjunction with (2.6) and (2.7).
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE
5
The above dual problems will be instrumental in this paper. To see why, we point out that taking = in (2.2) yields an estimate for k (u) ? (U )kL2 (Q), whereas taking = in (2.2) gives rise to an estimate for ku ? U kL1 (0;T ;H ?1( )); see x3. Based on this argument, we will present in x5 two rigorous a posteriori error estimates for problem (2.1) in these natural energy norms (Approaches I and II). In contrast to the heat equation, problems (2.3) and (2.4) do not exhibit any regularizing e ect; in fact as ! 0 the information on second derivatives is lost (e.g., if b = 0 in Q then, in the limit as ! 0, (; t) = () 2 H01 ( ) for all 0 t T ). Since the regularity of the dual problem dictates the weights (powers of meshsize and time step) of the a posteriori error estimators, this is an early indication of the striking di erence between degenerate and strictly parabolic problems. A simple calculation shows, however, that b a := a=(L + ) in the open set
Q := f(x; t) 2 Q : jU (x; t) ? L2 j > L2 + g; for any > 0. The proofs of Lemma 2.1, with splitting 2jb1=2 j bj j2 + j j2, and Lemma 2.2 yield Z
(2.9)
Q
j2 1
j
a
k k2 2
L (Q) ;
Z
Q
j j2 2a1 krk2L2( ):
Then L2t Hx2 seminorms of and in any compact subset Q of Q depend on the inverse of the parabolic distance between the parabolic boundaries of Q and Q . This lack of uniform regularity in L2t Hx2 prevents the construction of rigorous error indicators based on the heat equation away from discrete interfaces. However, in xx5 and 6 we present an empirical estimator (Approach III), which utilizes heat estimators in Q0 , and compare its performance with Approaches I and II in x7. 3. Error Representation Formulas
The purpose of this derivation is to obtain formulas for the errors
e (u) := (u) ? (U ); eu := u ? U; where u is the true solution and U 2 L2(Q) is any other function. We do not assume that U is any speci c approximation of u, and so the resulting formulas are quite general. In particular we do not require stability of U and (U ) in the energy p 2 1 2 norms L1 t Lx and Lt Hx , respectively. As a by-product we rederive the usual O( ") rate of convergence for a vanishing viscosity approximation U of u. For any function 2 C 0([0; T ]; H ?1 ( )), we denote t() := (; t) for all t 2 [0; T ]. We multiply the PDE operator J ( ) in (2.5) by eu(; t), which is in L2 ( ) for a.e. t 2 (0; T ), and use property beu = e (u) to write Z
Q
euJ ( ) =
Z
Q
?
u@t + (u) ?
Z
Q
?
U@t + (U ) +
Z
Q
eu:
Let U 0 ; U T 2 H ?1 ( ) be given, and set e0u := u0 ? U 0 and eTu := uT ? U T . At this stage, both U 0 and U T are arbitrary, but they will later be the initial value U (; 0)
6
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
and nal value U (; T ) of U which make no sense in the present context. We then integrate by parts in space and time the rst term on the right-hand side, employ (2.1), and add and subtract hU 0 ; 0i ? hU T ; T i, to arrive at (3.1)
heTu ; T i ?
Z
Q
euJ ( ) = he0 ; 0i + R( ) ?
Z
u
Q
eu;
where R( ), the parabolic residual, is the distribution (3.2)
R( ) := hU 0 ; 0i ? hU T ; T i +
Z
Q
?
f + U@t + (U ) :
Together with the initial error he0u ; 0i, R( ) is a measure of the amount by which U misses to be a solution of (2.1) and must be evaluated in negative norms. We show rst an L2t L2x error estimate for (u). To this end, let = be the solution of the dual problem (2.3) and note that T = 0 and
eu b1=2 = (eu e (u))1=2 A11=2 je (u)j: From (3.1) and (3.2), we easily get R
(3.3)
e
Q (u) ke (u)kL2(Q) = sup 2L2 (Q) k kL2 (Q) ke0ukH ?1( )kr (; 0)kL2 ( ) A1=2 sup k kL2 (Q) 2L2 (Q) jR( )j + A1=2 sup 2L2 (Q) k kL2 (Q) keukL2 (Q)k1=2 kL2 (Q) : + 1=2A1=2 sup k kL2 (Q) 2L2 (Q)
It is thus apparent from (2.6) with ! 0 that the representation formula (3.3) hinges on a negative norm of the residual R( ) involving rst derivatives of ,
jR( )j ; ?1 := sup 2L2 (Q) k kL2 (Q) and leads to the following a posteriori error estimates for (u). Lemma 3.1. Let be the solution of (2.3) with arbitrary 2 L2 (Q). Then
ke (u)kL2 (Q) A1=2
?
p12 ke0u kH ?1( ) + ?1 :
Proof. It remains to deal with the last term in (3.3). Since
jeuj L + a1 je (u)j;
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE
7
combining (2.6) with (3.3) implies ? ? (3.4) 1 ? A21a=2 1=2 ke (u)kL2 (Q) A1=2 p12 ke0u kH ?1( ) + ?1 + 21 LjQj1=2 1=2 : The assertion then follows upon taking ! 0. ?1 To derive an L1 t Hx error estimate for u, let = be the solution of the dual problem (2.4). Therefore, according with (3.1) and (3.2), since eu J ( ) = 0 we get R T T keu kH ?1( ) = sup1 kr keu2 L ( ) 2H0 ( ) 0 ( ) kr (; 0)kL2 ( ) sup1 keu kH ?1kr kL2( ) 2H0 ( ) (3.5) + sup1 krjRk( 2)j L ( ) 2H0 ( ) k eukL2 (Q)k1=2 kL2(Q) 1 = 2 + sup : krkL2( ) 2H01( ) Again, in view of (2.7) with ! 0, this representation formula hinges on a negative norm of the residual R( ) involving rst derivatives of , that is ?1 := sup1 krjRk( 2)j : L ( ) 2H0 ( ) Lemma 3.2. Let be the solution of (2.4) with an arbitrary 2 H01( ). Then keu(; T )kH ?1 ( ) ke0ukH ?1( ) + ?1: Proof. We make use of (2.7) and argue as in Lemma 3.1 to derive from (3.5) (3.6) keTu kH ?1( ) ke0ukH ?1( ) + ?1 + p12a 1=2 ke (u)kL2 (Q) + p12 LjQj1=2 1=2 : The assertion then follows upon taking ! 0. Lemmas 3.1 and 3.2 lead to Approach I below, and are pessimistic in that they deal with the worst scenario in terms of regularity of and , namely one space derivative. An alternative and fruitful avenue consists of keeping > 0, thereby allowing H 2 space regularity of and , and optimizing later without sending it to 0; this yields Approach II below and works best. To this end, we set
R?2 ; C := 1 ?A1=2 + p1 ; 0 := C aL 0 a 2 jQj1=2 0
p
where R?2 := A1=2 ?2 + ?2= 2 involves the following negative norms of the residuals R( ) and R( ) with two space derivatives of and jR( )j ; := sup jR( )j : ?2 := sup ?2 2L2 (Q) k kL2 (Q) 2H01( ) k kL2 (Q) We expect 0 to be small because it involves R( ) and R( ). However this cannot be guaranteed a priori and is re ected in the statement of next result.
8
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
Lemma 3.3. Let be the solution of (2.3) with an arbitrary 2 L2(Q) and let be the solution of (2.4) with an arbitrary 2 H01 ( ). Then p ke (u)kL2 (Q) + keu (; T )kH ?1( ) 2C0ake0ukH ?1( ) 8 ? < 2 C0 aLjQj1=2 R?2 1=2 if 0 12 C0 +:
if 0 > C102 :
2C0R?2
Proof. We add twice (3.3) with (3.5), and argue as with (3.4) and (3.6). Since jR( )j 1 jR( )j ; jR( )j p 1 jR( )j ; 1 = 2 1 = 2 k kL2 (Q) 2 k kL2 (Q) krkL2( ) 2 k kL2(Q) we readily get p ? (3.7) 2 ? C01=2 ke (u)kL2 (Q) + keu(; T )kH ?1( ) 2C0ake0u kH ?1( ) + q(); where q() := q? ?1=2 + q+ 1=2 with q? := R?2 and q+ := C0aLjQj1=2 . We now observe that = 0 minimizes q(). If C001=2 1, then the rst assertion follows trivially from q(0 ) = 2(q? q+ )1=2 . Otherwise, if C001=2 > 1, then q+ < C02q? and q(1=C02) < 2C0q? . This concludes the argument. If U is a nite element solution, then the last term in Lemmas 3.1{3.3 (parabolic residual) can be further evaluated via Galerkin orthogonality; this is accomplished in x5. We stress that U need not be a discrete solution as the following application of Lemma 3.3 illustrates. Let U be the solution of (2.1) with " (s) := (s) + "s in place of (s), namely U is the usual vanishing viscosity approximation of u. The proof below is di erent from the original one in [10], easily extends to other perturbations of such as that in [10], and is valid under minimal regularity of both u0 and f which precludes compactness. Corollary 3.1. The following perturbation estimate is valid for " small ke (u)kL2 (Q) + keukL1 (0;T ;H ?1( )) C"1=2; where constant C depends on ku0kH ?1( ); j j; T; L; a; A, and kf kL1 (0;T ;H ?1( )). R Proof. Since R( ) = ?" Q U , the assertion results from Lemma 3.3 upon realizing that R?2 aC0"kU kL2 (Q) C". To derive a bound on kU kL2 (Q) under minimal regularity, we consider (2.1) with " instead of and strong data U 0 2 L2 ( ) and f 2 L2(0; T ; H ?1 ( )). We take = GU (; t) 2 H01 ( ), where the Green's operator G is the inverse of ? in H01 ( ), and obtain 1 (3.8) 2 dt hGU; U i + h " (U ); U i = hGf; U i: We then integrate (3.8) in time for 0 t s T , the maximum of kU (; t)kH ?1( ) being attained at t = s. We easily deduce
(3.9)
2 1 1 2 kU (; s)kH ?1 ( ) + A
Z
0
s
k " (U )(; t)k2L2 ( )dt Z
1 0 2 2 kU kH ?1 ( ) + kU (; s)kH ?1 ( ) 0
s
kGf (; t)kH01 ( )dt;
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE
9
which yields a bound for kU (; s)kH ?1 ( ) without use of Gronwall's inequality. We insert this bound in (3.9), now for arbitrary s, and use that "0 (s) 0(s) a > 0 and is bounded to arrive at the desired estimate for kU kL2 (Q) for strong data. We nally observe that the mapping (U 0 ; f ) ! (U; " (U )) from H ?1 ( ) L1(0; T ; H ?1 ( )) ! L1(0; T ; H ?1 ( )) L2 (Q) is Lipschitz, which is easily seen with the technique leading to (3.9). This allows us to regularize the weak data U 0 2 H ?1 ( ) and f 2 L1(0; T ; H ?1 ( )) via U0 2 L2( ) and f 2 L2 (0; T ; H ?1 ( )) and eventually pass to the limit in the regularization parameter . 4. Finite Element Discretization
We now introduce the fully discrete problem, which combines continuous piecewise linear nite elements in space with backward di erences in time.Pn We denote by n the time step at the n-th step and set tn := i=1 i. Let N be the total number of time steps, that is tN T . For any function 2 C 0((tn?1; tn ]; H ?1 ( )), we denote n() := (; tn). We denote by Mn a uniformly regular partition of into simplices [3]. Mesh Mn is obtained from Mn?1 by re ning/coarsening, and thus Mn and Mn?1 are compatible. Given a triangle S 2 Mn, hS stands for its diameter and S for its sphericity and they satisfy hS 2S = sin( S =2), where S is the minimum angle of S . Uniform regularity of the family of triangulations is equivalent to S > 0, with independent of n. We also denote by Bn the collection of interior interelement boundaries e of Mn in ; he stands for the size of e 2 Bn. the usual space of piecewise linear nite elements Let Vn H01( ) indicate n over Mn. Let fxnk gKk=1 denote the interior nodes of Mn. Let I n : C00( ) ! Vn be the usual Lagrange interpolation operator, namely (I n)(xnk ) = (xnk) for all 1 k K n. Finally, let the discrete inner product h; in be de ned by
h'; in :=
Z
1 I n(')dx = d+1
X
S 2Mn
jS j
X
xnk 2S
'(xnk )(xnk ) 8 '; 2 Vn:
This corresponds to the vertex quadrature rule, which can be easily evaluated element by element and leads to mass lumping [3]. We observe the validity of the local estimates k' ? I n(')kL1 (S) 81 h2S jr'jjrj 8 '; 2 Vn; kr'kL2(S) 13 ?S 1k'kL2(S) 8 ' 2 Vn: (4.1) R Consequently, the element inner product h'; inS := S I n(') satis es [15] k'kL2 (S) (h'; 'inS )1=2 =: k'knS C k'kL2(S) 8 ' 2 Vn; R and the discrepancy between h'; iS := S ' and h'; inS can be bounded by jh'; iS ? h'; inS j 18 h2S kr'kL2(S)krkL2(S) (4.2) 241 h2S ?S 1k'kL2(S)krkL2(S) 8 '; 2 Vn: The discrete initial enthalpy U 0 2 V0 is de ned at nodes x0k of M0 = M1 to be (4.3) U 0 (x0k ) := u0(x0k ) 8 x0k 2 nF0; U 0 (x0k ) := 0 8 x0k 2 F0: Hence, U 0 is easy to evaluate in practice. Then we set 0 := I 0 (U 0 ).
10
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
Discrete problem. Given U n?1 ; n?1 2 Vn?1, then Mn?1 and n?1 are modi ed as described below to get Mn and n and thereafter U n ; n 2 Vn computed according to n = I n (U n ) and
(4.4)
n hU 1
n ? I n U n?1 ; 'in + hrn ; r'i = hI n f n ; 'in
8 ' 2 Vn:
In view of the constitutive relation n = I n (U n ) being enforced only at the nodes, and the use of mass lumping, (4.4) is easy to implement and yields a monotone operator in RKn . This problem is solved below via an optimized nonlinear SOR [12]. However, these computational tricks introduce further consistency errors that are apparent from (3.2). Whether these devices preserve optimal accuracy is still to be explored. We nally introduce further notation. The following sets will be used later
F n := [fS 2 Mn : n(x) = 0 for some x 2 S g discrete free boundary; T n := [fS 2 Mn : n(x)n?1 (x) 0 for some x 2 S g transition region; C n := [fS 2 Mn : S is coarsened from Mn?1g coarsening set: We point out that the compatibility of Mn and Mn?1 yields I nU n?1 6= U n?1 only in C n, in which case Vn?1 \ Vn 6= Vn?1. Let the jump Jen of rn across e 2 Bn be
Jen := [ rn] e e = (rnjS1 ? rnjS2 ) e: If the unit normal vector e to e always points from S2 to S1, then Jen is well de ned. For any element S 2 Mn, JSn stands for the jumps of rn across @S n@ . Let U be the piecewise constant extension of fU ng de ned by U (; 0) = U 0 () and U (; t) := U n () for all tn?1 < t tn with n 1. Let n () ? I n U n?1 () U n?1 < t tn ; n 1; 8 t (4.5) Ut(; t) := n and let the interior residual Rn be
Rn := I nf n ? Ut (; tn): 5. A Posteriori Error Analysis
We now state and prove two rigorous a posteriori error estimates. Theorem 5.1 is based on H 1 regularity of and , solutions of the dual problems (2.3) and (2.4), and leads to Approach I. Exploiting H 2 space regularity of and in the spirit of Lemma 3.3 yields Theorem 5.2 and corresponding Approach II. The derivation below parallels, and in fact extends, that in x3 but exploits Galerkin orthogonality to express negative norms of the residuals R( ) and R( ) in terms of computable quantities.
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE 11
Theorem 5.1 (Approach I). There exist constants Ci > 0 depending pon the minimum angle of Mn and the space dimension d, and C := A1=2 = 2 + 1 and
C] := A1=2 C0a, such that the following a posteriori error estimates holds
ke (u)kL2(Q) + keukL1 (0;T ;H ?1( )) E I (u0; f; T; ; U; h; ) := C
X
i6=4
Ei + C] E4 ;
where the error indicators Ei = EiI for 0 i 8 are given by
E0 := ku0 ? U 0 kH ?1( )
N X 1=2 X I hekJenk2L2 (e) E1 := C1 n n=1 e2Bn N X 1=2 X I h2S kRn k2L2(S) E2 := C2 n n=1 S 2Mn N
X E3I := nkr (U n) ? rI n (U n)kL2 ( ) n=1 N 1=2 X nkU n ? I nU n?1 k2L2( ) E4 := n=1
E5 :=
N X
n=1
E6 := C6 E7 := E8 :=
kI nU n?1 ? U n?1 kH ?1( )
N X
n=1
N X
n
S 2Mn
h4S krRnk2L2 (S)
nkf n ? I nf n kH ?1( )
n=1 N Z tn X n=1
X
tn?1
kf ? f n kH ?1( )
initial error; jump residual; interior residual; constitutive relation; time residual; coarsening;
1=2
quadrature; interpolation; time discretization :
All indicators Ei can be evaluated explicitly in terms of the computed solution U , initial datum u0, and source term f . Indicators E0 ; E1I ; E2I ; E4 , and E5 are essential, and are also present for the heat equation but with di erent weights and cumulative e ect in time [5]. The error accumulation is measured here in L1 for E1I ; E2I ; E5 and in L2 for E4 , whereas it is in L1 for the heat equation; the latter exhibits a weaker dependence on T for T 1. The powers of meshsize in E1I and E2I are smaller than those for the heat equation, namely h3e and h4S , respectively, thereby re ecting the degenerate nature of (1.1), or equivalently the lack of H 2 space regularity of and . The remainding indicators E3I ; E6 ; E7 ; E8 are not essential and could in principle be removed at the expense of complicating the implementation of (4.4). In particular, we note that (U n ) 6= I n (U n ) only for S F n and njS 6= 0, provided is piecewise linear, and that E6 C E2I as results from (4.1).
12
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
The following theorem yields the same weights for E1I and E2I as the heat equation but yet with a worse error accumulation in time. Note that this improvement comes at the expense of a smaller outermost power, namely 1=4 instead of 1=2. 2 )1=2 (2A1=2 + 1=p2); Cz := Theorem 5.2 (Approach II) . Let C := ( A + 1 =C y 0 p 2C0a, II II II 1 := E1 L+jQE2j1=+2 E3 ;
and
E1II := C1 E2II := C2 E3II
N X n=1 N X n=1
N X
:=
n=1
n n
X
e2Bn
h3e kJenk2L2 (e)
X
S 2Mn
1=2
h4S kRn k2L2(S)
1=2
jump residual; interior residual;
1=2 n n n 2 constitutive relation: nk (U ) ? I (U )kL2 ( )
The following a posteriori error estimate holds
ke (u)kL2 (Q) + keukL1 (0;T ;H ?1( )) E II (u0; f; T; ; U; h; ) := Cz E0 + p
8
C102 :
8 X
i=5
Ei
5.1. Residuals. We rst express the residual R( ), for a generic function 2 H 1;2(Q), in terms of computable quantities. We notice that, since U is piecewise constant in time, summation by parts yields
N X n n n ? 1 N N 0 0 hU; @t i = hU ; ? i = hU ? i ? hU ; i ? hU n ? U n?1; n?1i: 0 n=1 n=1
Z
N X
T
Moreover, in light of (4.5), the last term can be written equivalently as N X
hU n ? U n?1 ; n?1i =
n=1
Z
0
T
hUt ; i
N X n ? 1 h U ; ? i + hI nU n?1 ? U n?1; n?1i: + t n=1 n=1 tn?1 N Z tn X
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE 13
Rearranging terms, the residual R( ) becomes
R( ) = (5.1)
N Z tn ? X n=1 tn?1 N Z tn X
hRn ; i ? hrI n (U n); r i
N X n ? 1 hUt ; ? i ? hI nU n?1 ? U n?1 ; n?1i ? n ? 1 n=1 n=1 t n Z N Z tn N t X X n f n ; i: n n n h f ? I hr I ( U ) ? r ( U ) ; r i + + n=1 tn?1 n=1 tn?1
We next rewrite the discrete problem (4.4) as follows
hRn; 'i ? hrI n (U n ); r'i = hRn ; 'i ? hRn ; 'in 8 ' 2 Vn; and subtract this expression from the right hand side of the residual R( ) given by (5.1). This crucial step is usually referred to as Galerkin orthogonality. We nally decompose the integral hrn; r( ? ')i over all elements S 2 Mn and next integrate by parts to obtain the equivalent expression
?hrI n (U n); r( ? ')i =
X
e2Bn
hhJen; ? 'iie 8 ' 2 Vn;
where hh; iie denotes the L2 -scalar product on e 2 Bn. Thus we easily arrive at N Z tn X n ; ? 'i n h R hh J ; ? ' ii + R( ) = e e n=1 tn?1 n=1 tn?1 e2Bn N Z tn N Z tn X X n?1 ? i n n n h U ; hr I ( U ) ? r ( U ) ; r i ? + t n=1 tn?1 n=1 tn?1 N Z tn ? N X X n ; 'i ? hRn ; 'in n n ? 1 n ? 1 n ? 1 h R ? hI U ? U ; i + n=1 tn?1 n=1 N Z tn X hf ? I nf n ; i =: I + + VII: + n ? 1 n=1 t N Z tn X
X
Since we need to approximate under minimal regularity assumptions, we resort to the Clement interpolation operator n : L2( ) ! Vn, which satis es the following local approximation properties [4], for all 2 H k ( ) and k = 1; 2, (5.2)
~ kS jjH k(S~); k ? nkL2(S) + hS kr( ? n)kL2 (S) Ch ~ ek?1=2jjH k(S~); k ? nkL2(e) Ch
where S~ is the union of all elements surrounding S 2 Mn or e 2 Bn. Constants C~ depend solely on the minimum angle of the mesh Mn. An important by-product of
14
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
uniform mesh regularity is that the number of adjacent simplices to a given element is bounded by a constant M independent of n, meshsizes, and time steps. Hence X kk2L2(S~) M kk2L2 ( ) 8 2 L2( ): (5.3) S 2Mn
This, in conjunction with (5.2) for k = 1, yields ~ 1=2 )krkL2( ) 8 2 H01 ( ): (5.4) krnkL2 ( ) (1 + CM We now estimate each term I to VII separately to derive Theorems 5.1 and 5.2. We argue with = , solution of (2.3) with 2 L2(Q). Similar estimates are valid for = , solution of (2.4) with 2 H01 ( ). We make extensive use of the a priori estimates of Lemma 2.1 (and Lemma 2.2). Selecting '(; t) = n (; t) for tn?1 < t tn and using (2.6), (5.2), and (5.3), terms I and II can be bounded as follows
jIj C1 jIIj C2
N Z tn X X tn?1
n=1 N Z tn X
e2Bn
hekJenk2L2(e)
X
n=1 tn?1 S 2Mn
1=2
h2 kRnk2 2 S
L (S )
kr (; t)kL2 ( ) p12 E1I k kL2 (Q);
1=2
kr (; t)kL2 ( ) p12 E2I k kL2 (Q):
Moreover, (2.6) also yields
jIIIj
N Z tn X n=1
tn?1
kr (U n) ? rI n (U n )kL2 ( )kr (; t)kL2 ( ) p12 E3I k kL2(Q):
Estimators E1I ; E2I , and E3I are those de ned in Theorem 5.1; constants C1 and C2 depend upon C~ in (5.2) and M in (5.3). Alternatively, if E1II ; E2II , and E3II are the estimators in Theorem 5.2, in light of (2.8) we can also write
jIj C1 jIIj C2
N Z tn X X tn?1
n=1 N Z tn X
n=1
e2Bn
X
tn?1
S 2Mn
h3e kJenk2L2 (e)
1=2
h4S kRn k2L2(S)
j (; t)jH 2( ) 21 ?1=2 E1II k kL2 (Q);
1=2
j (; t)jH 2( ) 21 ?1=2 E2II k kL2 (Q);
and, upon integration by parts,
jIIIj
N Z tn X
n=1
tn?1
k (U n ) ? I n (U n )kL2 ( )k (; t)kL2 ( ) 21 ?1=2 E3II k kL2 (Q):
R With the aid of (2.6) and the fact that ? n?1 = ttn?1 @s , we readily obtain
jIVj
N X
n=1
1=2 Z T n n n ? 1 2 k@t nkU ? I U kL2 ( )
(A + )1=2 E4 k kL2 (Q):
0
k2L2 ( )
1=2
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE 15
Using (2.6) again, now combined with (4.2) and (5.4), we easily deduce
jVIj C6
N Z X n=1
tn
tn?1
X
S 2Mn
h4S krRnk2L2(S)
1=2
kr (; t)kL2 ( ) p12 E6k kL2 (Q);
where C6 depends on C~ and M . The remaining bounds follow directly from (2.6)
jVj jVIIj
N X
kI nU n?1 ? U n?1kH ?1( )kr
n=1 N X Z tn
n=1
tn?1
n?1 k
L2 ( ) p12 E5 k kL2 (Q) ;
kf ? I nf n kH ?1( )kr (; t)kL2 ( ) p12 (E7 + E8)k kL2 (Q):
Collecting all the previous estimates we get the following bound for R( ) (5.5)
8 X
jR( )j p1 E + (A + )1=2 E + 4 k kL2(Q) 2 i=5 i
(
p12 (E1I + E2I + E3I ) 1 ?1=2 ?E II + E II + E II : 1 2 3 2
Similarly, using (2.7) we obtain (
(5.6)
8 I I I ? jR( )j X 1=2 E + E1 + E2 +? E3 1 p A + ) E + 4 p12 ?1=2 E1II + E2II + E3II : krkL2( ) i=5 i 2
5.2. Proofs of Theorems 5.1 and 5.2. Upon adding A1=2 times (5.5) with (5.6) and taking ! 0, Theorem 5.1 is an easy consequence of Lemmas 3.1 and 3.2, because E I is increasing in T and tN T . In order to prove Theorem 5.2 we proceed as in Lemma 3.3. We take 1=C02 and resort to (3.7), which in the present context becomes (5.7)
8
X ke (u)kL2 (Q) + keNu kH ?1( ) Cz E0 + Ei + CyE4 + p12 Czq() i=5
with
q() := q? ?1=2 + q+ 1=2; q? := E1II + E2II + E3II ; q+ := LjQj1=2 : We then realize that = 1 optimizes q(), and argue as in Lemma 3.3. 5.3. Approach III. The local treatment of terms I, II, and III above, in conjunction with (2.9), suggests removing the factor ?1=2 away from discrete interfaces. This idea, heuristic in the sense that L2t Hx2 regularity of and is not necessarily uniform, yields a method somewhat in between Approaches I and II. We set the estimators E1III ; E2III ; E3III equal to those of Approach I near discrete interfaces and
16
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
equal to those of Approach II (heat estimators), but with factor D=C instead of 1 away from discrete interfaces, where D := a1=2 C0, and obtain ke (u)kL2 (Q) + keukL1(0;T ;H ?1( )) E III (u0; f; T; ; U; h; ) 8 X
:= C E0 +
i=5
?
Ei + C] E4 + C E1III + E2III + E3III :
To this end, we add A1=2 times (5.5) with (5.6), let ! 0, replace kD2 kL2(Q0 ) by k kL2(Q0) for = ; and use (2.9), and ultimately pretend that the lower bound of b is a instead of a0 = 0. If and were computable, we could keep these Hx2 norms and evaluate them locally, in which case they could be viewed as weights. 6. Adaptive Algorithm
The error estimators E (u0 ; f; T; ; U; h; ) of x5 entail an L1 or L2 norm in time, which is impractical in that the entire evolution history would be needed to control E . We resort here to an L1 norm in time, and explain our error equidistribution strategy, which is based upon minimizing the spatial degrees of freedom for a uniform error distribution in time. A similar strategy is derived in [1] for a linear elliptic problem. We then fully discuss element error indicators, and all tests necessary for mesh and time step admissibility. 6.1. Mesh Design. Let M(t) be a time-dependent mesh with variable meshsize h(x; t) and let (t) be the underlying variable time-step. Since hd is proportional to the volume of a generic space-time nite element, then the computational complexity of (4.4) can be measured by the total number of degrees of freedom (6.1)
M :=
Z
Q
where 0 < ? (x; t) + and
Z
M (t) :=
(6.2)
T
Z
h(x; t)?d (t)?1 (x; t) dxdt =
0
(t)?1 M (t) dt;
h(x; t)?d (x; t) dx
stands for the cardinality of M(t). Let E be a generic a priori error of the form (6.3)
E=
Z
0
T
E (t)dt =
Z
Q
?
h(x; t) + (t) E (x; t) dxdt:
Given an error tolerance ", we then pose the following question: optimize h and for an error distribution (6.4) E = ": We have to minimize (6.1) subject to the constraint (6.4). This constrained optimization problem is equivalent to seeking a saddle point of the functional
L(h; ; ) :=
Z
Q
h?d ?1 ? " ?
Z
Q
(h + )E ;
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE 17
with Lagrange multiplier . Di erentiation with respect to h and yield 0=
Z
Q
0=
@h L(h; ; ) =
Z
Q
Z
0
@ L(h; ; ) =
TZ Z
T
0
(?dh?(d+1) ?1 + h ?1E ); Z
(?h?d ?2 + ?1E );
where = (x; t) and = (t) are smooth functions. Hence Z d 1 ?1M (t): ? d ? 1 (6.5) h E = h ; E =
These relations, in conjunction with (6.1) and (6.3), yield Z Z T d + M: 1 ? d ? 1 " = h + ?1M (t) = d 0 Q We thus deduce the following expression for the Lagrange multiplier = d + M" : We conclude from (6.5) that the optimization procedure consists of equidistributing the local space-time errors according to the recipes +1 Z d M ? 1 d + E = + d ": (6.6) M h E = + d "; M (t)
A serious diculty of (6.6) is that M is never available for all times so as to make the optimization feasible. A di erent and more stringent objective would be to equidistribute pointwise discretization errors in (6.3). This leads to the local requirements
T j jh E = 1 ";
(6.7)
2
T
Z
E = 21 ";
which corresponds to having a constant integrant in (6.3) and yields (6.4). This method does not require M , and can thus be implemented. If we further impose = T=N to be constant, we end up with the method chosen in [11] for a priori mesh design. Upon combining (6.6) and (6.7), we get a third strategy that can still be implemented: we optimize the spatial degrees of freedom for a uniform error distribution E (t) in time. Minimizing (6.2) for each 0 < t T , subject to the constraints Z
h(x; t) E (x; t) dx = 1 ";
T (t)
Z
T E (x; t) dx = 12 "; 2
results in the following restrictions on h and , which are intermediate between (6.6) and (6.7), Z ? 1 d + 1 (6.9) TM (t)(x; t) h(x; t) E (x; t) = 2 "; T (t) E (x; t) dx = 12 ":
(6.8)
We adopt such a viewpoint here, and discuss its implementation next. Note rst that E would depend on the discrete solution, and so implicitly on h and , if E in (6.3) were an a posteriori error estimator. So the above analysis is a priori.
18
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
6.2. Equidistribution Strategy. The equidistribution strategy (ES) is an iterative procedure that in the kth step improves upon a mesh density hk (x; t) by means of the following two opposite operations, based on the rst constraint in (6.9). Let R ?d 0 < < 1 be a re nement factor and Mk = hk from (6.2) with h = hk . (a) Re nement. Set hk+1 = hk for all elements satisfying (6.10)
k := 2TMk ?1 hkd+ E > ":
(b) Coarsening. Set hk+1 = ?1hk for all elements verifying
k < 2d+ ":
(6.11)
In practice, more elements will be re ned and fewer will be coarsened to preserve mesh conformity. In fact, every local mesh modi cation involves also adjacent elements. Since Z
Z
d ?d h?d = ?d Mk ; Mk+1 = h?k+1 k
we see that k+1 k for re ned elements whereas k+1 ?2d? k < " for coarsened elements. We conclude that coarsened elements will not be candidates for re nement in the (k +1)th step, because k+1 < " and so test (a) fails, and that ?
k+1 min "; k : The latter guarantees convergence of ES in a nite number of steps. 6.3. Element Indicators. We now introduce error indicators of the form (6.8) with density E (x; t) computable element by element and denoted by E n(S ) for S 2 Mn and t = tn. We rst replace the global norm H ?1 ( ) in E0 ; E5 ; E7 ; E8 by the L2 ( ) norm scaled with the smallest diameter d of , which is the Poincare constant. We next show how to determine computable global indicators E0 (initial error), E ;1 ; E ;2 (temporal errors) and Eh;1; Eh;k 2 (spatial errors) satisfying (6.12)
E0 E0; E4 E ;1; E8 E ;2; E5 + E6 + E7 Eh;1 ; E1k + E2k + E3k Eh;k 2 ;
where k = I; II stands for the approach. Each of these errors Ep is obtained from element indicators Epn(S ) via
Ep = 1max nN
X
S 2Mn
Epn(S )
1=2
;
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE 19
where the max is super uous for p = 0. The element indicators are given by
E0(S ) := d2 jS jku0 ? U 0 k2L1 (S) E ;n 1(S ) := T kU n ? I nU n?1 k2L2(S) E ;n 2(S ) := T 2d2 kf ? f nk2L1 (tn?1;tn;L2 (S)) Eh;n 1(S ) := 3T 2 d2 n?2kI nU n?1 ? U n?1k2L2 (S) + C62h4S krRnk2L2 (S) + d2 kf n ? I nf n k2L2(S) I ;n 2 Eh;2(S ) := 3T 21 C12hS kJSnk2L2 (@S) + C22h2S kRnk2L2(S) + kr (U n) ? rI n (U n )k2L2 (S) Eh;II;n2 (S ) := 3T 21 C12h3S kJSnk2L2(@S) + C22h4S kRnk2L2(S) n n n 2 + k (U ) ? I (U )kL2 (S)
local initial error; local time residual; local time discretization; local coarsening local quadrature local interpolation; local jump residual I local interior residual I local constitutive relation I; local jump residual II local interior residual II local constitutive relation II:
Since u0 2 W 1;1( nF0), interpolation gives ku0 ? U 0 kL1 (S) ChS kru0kL1 (S) = O(hS ) away from F0 . If S intersects F0, instead, then u0 has a jump of size L which yields ku0 ? U 0 kL1 (S) = O(L). In both instances, the ensuing orders match that of the local jump residual I, which indicates that M0 will not be drastically modi ed. If in addition u0 2 W 2;1( nF0), then ku0?U 0 kL1(S) Ch2S kD2 u0kL1 (S) = O(h2S ) d away from F0 , which yields the order O(h4+ S ) of the local jump residual II. According to x5.3, we introduce the local indicators for heuristic Approach III which is based on combining Eh;I;n2(S ) and Eh;II;n2 (S ). With D = a1=2 C0 we set
Eh;III2;n(S ) :=
(
Eh;I;n2(S ) if S T n D22 E II;n (S ) if S 6 T n : C h;2
A simple but tedious calculation shows that (6.12) is ful lled, and the following inequality is valid for Approaches k = I and III (heuristic for Approach III) (6.13)
E k (u0; f; T; ; U; h; ) E k := C (E0 + E ;2 + Eh;1 + Eh;k 2 ) + C]E ;1:
On the other hand, if 1=C02 and q() := Eh;II 2?1=2 + LT 1=2j j1=2 1=2, then (5.7) and (6.12) lead to (6.14) E II (u0 ; f; T; ; U; h; ) E II := Cz(E0 + E ;2 + Eh;1) + CyE ;1 + p12 Czq():
20
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
6.4. Adaptive Strategy. Given an error tolerance ", the objective is to adap-
tively select time steps and mesh densities h so that the spacial degrees of freedom are optimized for a uniform error distribution in time and (6.15) E k(u0; f; T; ; U; h; ) ": To this end we equidistribute the local estimators of x6.3 according to (6.9) via ES, and use bisection to perform re nement/coarsening operations (6.10) and (6.11). Bisection creates compatible consecutive meshes, extends naturally from 2D to 3D, and is handy for combined re nement/coarsening operations; we refer to [2]. Given re nement parameters ? > 0 and coarsening parameters > 0 satisfying ?0 + ? + ?h 1; < ? ; h < ?h ; then time steps and mesh densities for Approaches k = I and III are reduced until CE0 ?0 "; 1 1 (6.16) 2 " C] E ;1 ; C E ;2 2 ? "; 1 h " C Eh;1 ; C E k 1 ?h ": h;2 2 2 In fact, from (6.13) we readily obtain (6.15). Achieving (6.15) for Approach II is more problematic in that the optimal choice of in (6.14) turns out to be Eh;II 2 2 := LT 1=2j j1=2 ; which is not a computable quantity as it involves a maximum over all time steps. To overcome this diculty we restrict the meshsize and time step as follows CzE0 ?0 "; 1 " CyE ;1 ; CzE ;2 1 ? "; 2 2 (6.17) 1 h " Cz Eh;1 1 ?h "; 2 2 1 G h " C 2aE II 1 G?h "; 0 h;2 4 4 with G := min 1; 4aLT?1h=2"j j1=2 : We claim this yields (6.15). In fact, if G < 1 then 2 2 Eh;II 2 2 = LT 1=2j j1=2 16C 2?ah2 L" 2T j j = CG022 < C102 ; 0 whence II )1=2 1 ?h ": p12 Czq (2 ) = 2C0 aL1=2 T 1=4 j j1=4 (Eh; 2 2 2 On the other hand, if G = 1, then with = 1=C0 in (6.14) we deduce ? II + aLT 1=2 j j1=2 1 ? " + 1 ? " = 1 ? ": p12 Cz q C12 = C02 aEh; 2 4 h 4 h 2 h 0 We now describe an algorithm based on (6.16) for Approach I (and III). The implementation of Approach II is similar, but based on (6.17) instead. Let M n denote the cardinality of Mn at any step of ES.
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE 21
6.4.1. Initial Mesh. Given a coarse mesh, ES bisects all S 2 M0 such that 2 "2 ? 0 2 C E0(S ) > M 0 ;
as suggested by (6.10) without T . We interpolate u0 by means of (4.3) to nd U 0 (and 0) and estimate E0(S ) as E0(S ) h2S jS jkrU 0k2L1 (S). 6.4.2. Time Step Selection. Starting with n = n?1, the algorithm checks whether
E n := 4 max C]2
X
S 2Mn
E ;n 1(S ); C2
X
S 2Mn
E ;n 2(S ) > ? 2 "2 ; E n < 2 "2 :
In the rst case n is reduced, whereas in the second one (corresponding to n being too small) n is accepted but the initial guess for the next time step size is enlarged. 6.4.3. Mesh Size Selection. Starting from Mn = Mn?1, ES checks whether 2 2 2 2 ? EhI;n(S ) := 4C2 max Eh;n 1 (S ); Eh;I;n2(S ) > ?Mh "n ; EhI;n(S ) < Mh "n ;
as suggested by (6.10) and (6.11). Then re nement and coarsening operations are performed accordingly, with the precaution of choosing h ?h properly to prevent ES from alternating such operations over the same elements. 6.4.4. Flow Chart. Figure 6.1 shows a ow diagram of the adaptive algorithm for Approach I. After each iteration of ES, both U n and n are recalculated on the new mesh using the new time step size. To minimize the overall computational cost, a compromise is reached between the minimization of the number of mesh elements (degrees of freedom) and iterations of ES. ES stops iterating as soon as X
S 2Mn
EhI;n(S ) ?h2 "2
is ful lled, thereby allowing some elements to violate the local error tolerance. Consequently, discretization errors might not be equidistributed. As implemented, and shown in the ow diagram, at least one mesh modi cation per time step is performed in order to permit elements near the moving interface to be re ned, even if the global error bound is already ful lled by the old mesh. 6.5. Convergence. We elaborate here on convergence of the algorithm of x6.4; a more comprehensive study will be given elsewhere [13]. For all 1 n N , we assume the following discrete a priori estimates, which are slightly stronger than the natural bounds
krnkL2 ( ) = O(1); jF nnZ nj = o(1); X kU n ? I nU n?1k2L2 ( ) + n?1kn ? I nn?1k2L2 ( ) + hS kJSnk2L2(@S) = o(1); S nF n
22
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
Figure 6.1. Flow diagram of the adaptive Approach I. where Z n := [fS 2 Mn : njS = 0g is the numerical mush. We now examine a few
relevant terms, using the hyperbolic and parabolic relations h2S C ; rn := max h2S C : (6.18) rHn := Smax P H P S nT n n T n n2 We consider rst the jump residual, for which it suces to estimate the contribution
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE 23
of F nnZ n, namely X
eF nnZ n
hekJenk2L2 (e) C krnk2L2 (F nnZ n) = o(1):
The same calculation applies to the constitutive relation. For the interior residual, we decompose the integral over T n and the complement, where we use U n ? I nU n?1 = n ? I nn?1, to deduce
h2S kU n ? I nU n?1 k2 rn kU n ? I nU n?1 k2 H L2 (S ) L2 ( ) 2 S 2Mn n n + r P kn ? I nn?1k2L2( ) = o(1): n X
Finally, kU n ? I nU n?1 kL2( ) = o(1) also controls the time residual. This demonstrates that the a posteriori estimators of x6.4 tend to 0 as both meshsize and time-step approach 0. Therefore the goal (6.15) is achievable. On the other hand, (6.18) re ects the hyperbolic structure of the interface as well as the parabolic structure of the problem elsewhere. We refer to [11] for a similar observation. 7. Simulations
We present results of the adaptive method for the model problem from [12] with known exact solution (
?
0:75 6r2 ? 61 if r < 61 ; u(x; y; t) := ? 1 + 1:5 ? 0 (t) y? r (t) (r ? 16 ) if r 61 ;
(7.1)
where r := (x ? 13 )2 + (y ? (t))2 1=2. The interface is a circle of constant radius r = 61 centered at ( 13 ; (t)) with (t) := 0:5 + 0:1 sin(12:5t), and thus oscillates in the vertical direction with highly varying normal velocity. Initial datum, boundary values, and heat source f are directly computed from (7.1). The domain is = (0:1; 0:7) (0:2; 0:8) and the nal time is T = 1. Note that f is discontinuous across the true interface. Therefore a simple minded use of I nf leads to large errors which dominate the local indicators. However a discontinuous f is not realistic in practice, being just the cost of having a simple exact solution. To examine the essense of the proposed methodology, a special quadrature is used for those triangles crossed by the exact free boundary: their intersection is determined rst and then separate quadrature used in the resulting quadrilateral and triangle. The nonhomogeneous Dirichlet boundary condition is handled in the standard fashion since the interface does not touch the xed boundary. Several simulations were carried out with error tolerances " = 20; 14; 10; 7; 5. The following parameters were used for partitioning the total error into initial, temporal, and spatial components ?
?0 = ? = 0:2; ?h = 0:6; = 0:155; h = 0:2683 (I); 0:1897 (II, III):
24
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
Figure 7.1. Approaches I with " = 14(10) and II and III with " =
10(7); Triangle count: I 1016 (2747), II 700 (1592), III 786 (2129)
Figure 7.2. Meshes and isotherms for Approach II with " = 7
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE 25
Figure p p 7.3. E ect of coarsening for Approach I and " = 10, h =?h =
0:2( 0:1).
Figure 7.4. Total element count vs true error for Approaches I, II,
and III.
Since the expected local coarsening error is proportional to the power of meshsize, and such a power is smaller for Approach I than for II and III, more coarsening (larger h) is allowed for Approach I. The results are summarized in Figures 7.1 to 7.4. Figure 7.1 shows the meshes at time t = 0:3 produced by Approaches I with " = 14 (10), and II and III with " = 10 (7). Meshes from left to right correspond to Approaches I, II, and III. The true errors for all three pictures on top and botton are similar. Approach I leads to more triangles than II and III because of the lower power of meshsize in the estimators. Approaches II and III are comparable. Figure 7.2 shows three meshes generated by Approach II for quite di erent interface velocities. In the left and middle pictures the velocity is highest, and the
26
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
interface is moving upwards (t = 0:5) and downwards (t = 0:765) as re ected in the isotherms (the interface is enhanced with a bold line). The rightmost picture corresponds to vanishing velocity and so to circular isotherms (t = 0:876). Higher interface velocities lead to higher re nement near the free boundary, because velocity is proportional to the temperature gradient jump. Higher velocities yield also larger temperature gradients in the liquid phase and corresponding additional mesh re nement. The triangle counts are 1670, 1734, and 1366, respectively. Figure 7.3 illustrates the e ect of h for Approach I with the interface moving downwards. The ratio ?h = h cannot be close to 1 nor too small. The former situation leads to mesh oscillations due to repeated coarsening/re nement operations over the same elements. A small ratio instead leads to meshes with unnecessary triangles, which thereby re ect the evolution history. This is an undesirable event. In Figure 7.4 we provide a quantitative comparison: a plot of the total element count (added over all time steps), which is proportional to the total amount of work, as a function of the true error. Approach II performs slightly better in the range tested. References 1. I. Babuska and W.C. Rheinboldt, Error estimates for adaptive nite element computations, SIAM J. Numer. Anal. 15 (1978), 736{754. 2. E. Bansch, Local mesh re nement in 2 and 3 dimensions, IMPACT Comput. Sci. Engrg. 3 (1991), 181{191. 3. P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. 4. Ph. Clement, Approximation by nite element functions using local regularization, RAIRO Anal. Numer. 9 (1975), 77{84. 5. K. Eriksson and C. Johnson, Adaptive nite element methods for parabolic problems I: a linear model problem, SIAM J. Numer. Anal. 28 (1991), 43{77. 6. , Adaptive nite element methods for parabolic problems II: optimal error estimates in L1 L2 and L1 L1 , SIAM J. Numer. Anal. 32 (1995), 706{740. 7. , Adaptive nite element methods for parabolic problems IV: nonlinear problems, SIAM J. Numer. Anal. 32 (1995), 1750{1763. 8. A. Friedman, Variational Principles and Free Boundary Problems, Wiley, New York, 1982. 9. O.A. Ladyzenskaja, V. Solonnikov, and N. Ural'ceva, Linear and Quasilinear Equations of Parabolic Type, vol. TMM 23, AMS, Providence, 1968. 10. R.H. Nochetto, Error estimates for multidimensional singular parabolic problems, Japan J. Indust. Appl. Math. 4 (1987), 111{138. 11. R.H. Nochetto, M. Paolini, and C. Verdi, An adaptive nite elements method for two-phase Stefan problems in two space dimensions. Part I: stability and error estimates. Supplement, Math. Comp. 57 (1991), 73{108, S1{S11. 12. , An adaptive nite elements method for two-phase Stefan problems in two space dimensions. Part II: implementation and numerical experiments, SIAM J. Sci. Statist. Comput. 12 (1991), 1207{1244. 13. R.H. Nochetto, A. Schmidt, and C. Verdi, Mesh and time step modi cation for degenerate parabolic problems, in preparation. , Adaptive algorithm and simulations for Stefan problems in two and three dimensions, 14. in preparation. 15. P.A. Raviart, The use of numerical integration in nite element methods for solving parabolic equations, Topics in Numerical Analysis (J.J.H. Miller ed.), Academic Press, London, 1973, pp. 233{264. Department of Mathematics and Institute for Physical Science and Technology, University of Maryland, College Park, MD 20742, USA.
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE 27 E-mail address: [email protected] Institut fur Angewandte Mathematik, Universitat Freiburg, 79106 Freiburg, Germany.
E-mail address: [email protected]
Dipartimento di Matematica, Universita di Milano, 20133 Milano, Italy.
E-mail address: [email protected]
parabolic PDEs, which are based on evaluating a parabolic residual in negative norms. The resulting upper bounds are valid for any numerical method, and rely on regularity properties of solutions of a dual parabolic problem in nondivergence form with vanishing di usion coecient. They are applied to a practical space-time discretization consisting of C 0 piecewise linear nite elements over highly graded unstructured meshes, and backward nite di erences with varying time-steps. Two rigorous a posteriori error estimates are derived for this scheme, and used in designing an ecient adaptive algorithm, which equidistributes space and time discretization errors via re nement/coarsening. A simulation nally compares the behavior of the rigorous a posteriori error estimators with a heuristic approach, and hints at the potentials and reliability of the proposed method.
1. Introduction
A posteriori error estimates are a fundamental component in the design of reliable and ecient adaptive algorithms for the numerical solution of PDEs. Even though rigorous results are available for linear and mildly nonlinear parabolic PDEs [5], [6], [7], the theory is much less satisfactory for strongly nonlinear PDEs. There are no results applicable to degenerate parabolic PDEs, which in turn exhibit lack of regularity across interfaces and corresponding numerical pollution e ects. The use of highly graded meshes and varying time-steps is thus motivated by the nonlinear structure of the PDE, as opposed to domain geometry, and is a vehicle for resolving small scale features with optimal computational complexity. In this paper we introduce a rigorous theory of a posteriori error estimation for degenerate parabolic problems of the form (1.1)
@t u ? (u) = f in Q := (0; T );
where is nondecreasing and Lipschitz. A typical example of industrial interest is that of solidi cation (classical Stefan problem), for which (1.2)
(s) := ? min(s; 0) + + max(s ? L; 0):
1991 Mathematics Subject Classi cation. 65N15, 65N30, 65N50, 80A22, 35K65, 35R35. Key words and phrases. degenerate parabolic equations, Stefan problem, nite elements, parabolic duality, a posteriori estimates, adaptivity. This work was partially supported by NSF Grants DMS-9305935 and DMS-9623394, EU Grant HCM \Phase Transitions and Surface Tension", MURST, CNR Contract 95.00735.01 and TMR contract \Viscosity Solutions and their Applications". Typeset by AMS-TEX 1
2
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
This corresponds to an ideal material with constant thermal coecients ? ; + > 0 and latent heat L. Any approximation U of u satis es (1.3) @t U ? (U ) = f ? R in Q; where R, an oscillatory distribution of singular character, is the so-called parabolic residual. In spite of the simple parabolic structure of (1.1), (1.2), its degenerate nature makes the theory of [5] fail in that it exploits the regularizing e ect of a linear parabolic dual problem. The corresponding dual PDE in this context is the nonstrictly parabolic equation in nondivergence form (1.4) @t + b = in Q; with vanishing and rough di usion coecient 0 b max( ? ; +) [8], [9], [10]. Such a PDE does not exhibit a regularizing mechanism, and is not computable in that b is discontinuous and depends on both u and U . Problem (1.4) measures the error accumulation in time, and is thus crucial. Our objective is to prove global regularity properties of in x2, and use them in x3 to derive two representation formulas for the errors u ? U and (u) ? (U ) in energy norms. These formulas are valid for any numerical method, evaluate R in two distinct negative norms, and lead to rigorous a posteriori upper error bounds. Since negative norms entail averaging, they are appropriate to quantify the oscillatory character of R. We next apply these ideas to a practical scheme consisting of C 0 piecewise linear nite elements over highly graded unstructured meshes and backward nite di erences with varying time-steps. The method uses mass lumping and evaluates (U ) solely at the nodes, which makes it easy to implement and solve iteratively. We discuss the method in x4 and derive in x5 two rigorous a posteriori error estimates for it of the form (Approaches I and II): (1.5) ku ? U kL2 (Q) + k (u) ? (U )kL1 (0;T ;H ?1( )) E (u0 ; f; T; ; U; h; ): The estimator E is computable in terms of data u0 = u(; 0); f; T; , computed solution U , meshsize h, and time step , but entails L1 or L2 norms in time. This is impractical in that the entire evolution history would be needed to control E . We thus resort to an L1 norm in time, and introduce an adaptive algorithm which equidistributes space discretization errors for a uniform error distribution in time. Such errors are estimated via local a posteriori error indicators, and further equidistributed via a re nement/coarsening strategy based on bisection. This yields compatible consecutive meshes. These issues are fully discussed in x6. We conclude in x7 with simulations illustrating the viability of our Approaches I and II, as well as a heuristic Approach III based on using local regularity of in (1.4) and heat estimators away from discrete interfaces. We clearly show that they are all able to detect the presence of interfaces, and re ne accordingly, and that Approaches II and III perform best. There is no need to compute the interface explicitly for mesh design, which is a major improvement with respect to [12]. Further simulations, comparisons of several nonlinear solvers, and a detailed description of the adaptive algorithm will be presented elsewhere [14]. Even though our error estimates (1.5) are rigorous, they do not necessarily imply E (u0 ; f; T; ; U; h; ) ! 0 as h; ! 0 because E depends on discrete quantities that change with h and . Stability and convergence will be assessed in [13].
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE
3
2. Setting and Parabolic Duality
Let Rd (d > 1) be a bounded convex polyhedral domain, let T > 0 be xed, and set Q := (0; T ). Let 2 W 1;1(R) satisfy (s) = 0 for all s 2 (0; L) and 0 < a 0(s) A for a.e. s 2= (0; L). Let u0 indicate the initial enthalpy, let 0 := (u0) denote the initial temperature, and let F0 := fx 2 : 0(x) = 0g be the initial interface. They satisfy 0 2 W01;1( ); F0 is a Lipschitz curve: Therefore, u0 2 W 1;1( nF0 ) and u0 has a jump discontinuity across F0. Finally let f be suciently smooth. The continuous problem then reads as follows. Continuous problem. Find u and such that 2 L2(0; T ; H01 ( )); u 2 L1(0; T ; L2 ( )) \ H 1(0; T ; H ?1 ( )); (x; t) = (u(x; t)) a.e. (x; t) 2 Q; u(; 0) = u0; and for a.e. t 2 (0; T ) and all 2 H01 ( ) the following equation holds (2.1) h@t u; i + hr; ri = hf; i: Hereafter, h; i stands for either the inner product in L2( ) or the duality pairing between H ?1 ( ) and H01( ). It is to be observed that the vanishing Dirichlet boundary condition on is assumed only for simplicity and so that the interface F (t) := fx 2 : (x; t) = 0g does not include @ . Existence and uniqueness for this problem are known [8], [9]. To motivate the dual problem (1.4), already studied in [8], [9], [10], we subtract (1.3) from (1.1) and integrate by parts over Q against a smooth test function vanishing on @ (0; T ). The error eu := u ? U satis es Z
T
heu ; ijt=T = heu ; ijt=0 + heu; @t + b i + R( ); 0 where 0 b(x; t) A is the discontinuous function ( (u(x;t))? (U (x;t)) if u(x; t) = 6 U (x; t) u(x;t)?U (x;t) b(x; t) :=
(2.2)
A otherwise: We could thus represent norms of eu(; T ) or their integrals over Q by making judicious choices of (; T ) and @t + b . Evaluation of R( ) depends on regularity of , which we now investigate. Given a regularization parameter > 0 to be chosen later, we consider two backward parabolic problems (2.3) J ( ) = ?b1=2 in Q; (; T ) = 0 in ; (2.4) J ( ) = 0 in Q; (; T ) = in ; with operator J in nondivergence form (2.5) J ( ) := @t + (b + ); ; = 0 on @ (0; T ), and 2 L2(Q), 2 H01 ( ). Existence of unique solutions ; 2 H 1;2 (Q) follows directly from the theory of nonlinear strictly parabolic problems [9].
4
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
Lemma 2.1. The following a priori bounds are valid for all 0 t T (2.6) 2kr (; t)k2L2 ( ); A+1 k@t k2L2 (Q); 4k k2L2 (Q) k k2L2 (Q): Proof. We multiply (2.3) by and integrate by parts in (t; T ) to obtain Z
1 2
jr
(; t)j2 +
Z Z
t
T
bj
Z Z
T
(b + )j
j2 = ?
t Z Z T 2 1 j j2: j +4
t
Z Z
t
T
b1=2
This yields the a priori bounds for r and . Finally, we multiply (2.3) by @t =(b + ) and integrate by parts in (t; T ) to get Z Z
T
Z
Z Z
T
b1=2 @s 1 j@s j2 + 1 jr (; t)j2 = ? b + 2
t b+
t Z Z T Z Z T 2 1 1 1 j@s j + 2 j j2; 2 b +
t
t
because > 0. This implies the a priori bounds for @t because b A. Lemma 2.2. The following a priori bounds are valid for all 0 t T (2.7)
kr (; t)k2L2 ( ); A+2 k@t k2L2 (Q); 2k k2L2(Q) krk2L2( ):
Proof. We multiply (2.4) by and integrate by parts in (t; T ) to get Z
Z Z
1 jr (; t)j2 + 2
t
T
Z
(b + )j j2 = 21 jrj2:
The a priori bound for thus follows from b 0. In view of (2.4) we further have
k@t k2L2(Q) = k(b + ) k2L2(Q) (A + )k(b + )1=2 k2L2(Q) A2+ krk2L2( ); where we have used that b A. This completes the proof. Corollary 2.1. The following L2t Hx2 a priori bounds are valid (2.8)
Z
0
T
j
j2 2
H ( )
2 1 4 k kL2 (Q) ;
Z
0
T
j j2H 2( ) 21 krk2L2( ):
Proof. It suces to invoke the well-known estimate [5], [9] for convex
jjH 2( ) kkL2( ) 8 2 H01 ( ) \ H 2 ( ); in conjunction with (2.6) and (2.7).
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE
5
The above dual problems will be instrumental in this paper. To see why, we point out that taking = in (2.2) yields an estimate for k (u) ? (U )kL2 (Q), whereas taking = in (2.2) gives rise to an estimate for ku ? U kL1 (0;T ;H ?1( )); see x3. Based on this argument, we will present in x5 two rigorous a posteriori error estimates for problem (2.1) in these natural energy norms (Approaches I and II). In contrast to the heat equation, problems (2.3) and (2.4) do not exhibit any regularizing e ect; in fact as ! 0 the information on second derivatives is lost (e.g., if b = 0 in Q then, in the limit as ! 0, (; t) = () 2 H01 ( ) for all 0 t T ). Since the regularity of the dual problem dictates the weights (powers of meshsize and time step) of the a posteriori error estimators, this is an early indication of the striking di erence between degenerate and strictly parabolic problems. A simple calculation shows, however, that b a := a=(L + ) in the open set
Q := f(x; t) 2 Q : jU (x; t) ? L2 j > L2 + g; for any > 0. The proofs of Lemma 2.1, with splitting 2jb1=2 j bj j2 + j j2, and Lemma 2.2 yield Z
(2.9)
Q
j2 1
j
a
k k2 2
L (Q) ;
Z
Q
j j2 2a1 krk2L2( ):
Then L2t Hx2 seminorms of and in any compact subset Q of Q depend on the inverse of the parabolic distance between the parabolic boundaries of Q and Q . This lack of uniform regularity in L2t Hx2 prevents the construction of rigorous error indicators based on the heat equation away from discrete interfaces. However, in xx5 and 6 we present an empirical estimator (Approach III), which utilizes heat estimators in Q0 , and compare its performance with Approaches I and II in x7. 3. Error Representation Formulas
The purpose of this derivation is to obtain formulas for the errors
e (u) := (u) ? (U ); eu := u ? U; where u is the true solution and U 2 L2(Q) is any other function. We do not assume that U is any speci c approximation of u, and so the resulting formulas are quite general. In particular we do not require stability of U and (U ) in the energy p 2 1 2 norms L1 t Lx and Lt Hx , respectively. As a by-product we rederive the usual O( ") rate of convergence for a vanishing viscosity approximation U of u. For any function 2 C 0([0; T ]; H ?1 ( )), we denote t() := (; t) for all t 2 [0; T ]. We multiply the PDE operator J ( ) in (2.5) by eu(; t), which is in L2 ( ) for a.e. t 2 (0; T ), and use property beu = e (u) to write Z
Q
euJ ( ) =
Z
Q
?
u@t + (u) ?
Z
Q
?
U@t + (U ) +
Z
Q
eu:
Let U 0 ; U T 2 H ?1 ( ) be given, and set e0u := u0 ? U 0 and eTu := uT ? U T . At this stage, both U 0 and U T are arbitrary, but they will later be the initial value U (; 0)
6
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
and nal value U (; T ) of U which make no sense in the present context. We then integrate by parts in space and time the rst term on the right-hand side, employ (2.1), and add and subtract hU 0 ; 0i ? hU T ; T i, to arrive at (3.1)
heTu ; T i ?
Z
Q
euJ ( ) = he0 ; 0i + R( ) ?
Z
u
Q
eu;
where R( ), the parabolic residual, is the distribution (3.2)
R( ) := hU 0 ; 0i ? hU T ; T i +
Z
Q
?
f + U@t + (U ) :
Together with the initial error he0u ; 0i, R( ) is a measure of the amount by which U misses to be a solution of (2.1) and must be evaluated in negative norms. We show rst an L2t L2x error estimate for (u). To this end, let = be the solution of the dual problem (2.3) and note that T = 0 and
eu b1=2 = (eu e (u))1=2 A11=2 je (u)j: From (3.1) and (3.2), we easily get R
(3.3)
e
Q (u) ke (u)kL2(Q) = sup 2L2 (Q) k kL2 (Q) ke0ukH ?1( )kr (; 0)kL2 ( ) A1=2 sup k kL2 (Q) 2L2 (Q) jR( )j + A1=2 sup 2L2 (Q) k kL2 (Q) keukL2 (Q)k1=2 kL2 (Q) : + 1=2A1=2 sup k kL2 (Q) 2L2 (Q)
It is thus apparent from (2.6) with ! 0 that the representation formula (3.3) hinges on a negative norm of the residual R( ) involving rst derivatives of ,
jR( )j ; ?1 := sup 2L2 (Q) k kL2 (Q) and leads to the following a posteriori error estimates for (u). Lemma 3.1. Let be the solution of (2.3) with arbitrary 2 L2 (Q). Then
ke (u)kL2 (Q) A1=2
?
p12 ke0u kH ?1( ) + ?1 :
Proof. It remains to deal with the last term in (3.3). Since
jeuj L + a1 je (u)j;
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE
7
combining (2.6) with (3.3) implies ? ? (3.4) 1 ? A21a=2 1=2 ke (u)kL2 (Q) A1=2 p12 ke0u kH ?1( ) + ?1 + 21 LjQj1=2 1=2 : The assertion then follows upon taking ! 0. ?1 To derive an L1 t Hx error estimate for u, let = be the solution of the dual problem (2.4). Therefore, according with (3.1) and (3.2), since eu J ( ) = 0 we get R T T keu kH ?1( ) = sup1 kr keu2 L ( ) 2H0 ( ) 0 ( ) kr (; 0)kL2 ( ) sup1 keu kH ?1kr kL2( ) 2H0 ( ) (3.5) + sup1 krjRk( 2)j L ( ) 2H0 ( ) k eukL2 (Q)k1=2 kL2(Q) 1 = 2 + sup : krkL2( ) 2H01( ) Again, in view of (2.7) with ! 0, this representation formula hinges on a negative norm of the residual R( ) involving rst derivatives of , that is ?1 := sup1 krjRk( 2)j : L ( ) 2H0 ( ) Lemma 3.2. Let be the solution of (2.4) with an arbitrary 2 H01( ). Then keu(; T )kH ?1 ( ) ke0ukH ?1( ) + ?1: Proof. We make use of (2.7) and argue as in Lemma 3.1 to derive from (3.5) (3.6) keTu kH ?1( ) ke0ukH ?1( ) + ?1 + p12a 1=2 ke (u)kL2 (Q) + p12 LjQj1=2 1=2 : The assertion then follows upon taking ! 0. Lemmas 3.1 and 3.2 lead to Approach I below, and are pessimistic in that they deal with the worst scenario in terms of regularity of and , namely one space derivative. An alternative and fruitful avenue consists of keeping > 0, thereby allowing H 2 space regularity of and , and optimizing later without sending it to 0; this yields Approach II below and works best. To this end, we set
R?2 ; C := 1 ?A1=2 + p1 ; 0 := C aL 0 a 2 jQj1=2 0
p
where R?2 := A1=2 ?2 + ?2= 2 involves the following negative norms of the residuals R( ) and R( ) with two space derivatives of and jR( )j ; := sup jR( )j : ?2 := sup ?2 2L2 (Q) k kL2 (Q) 2H01( ) k kL2 (Q) We expect 0 to be small because it involves R( ) and R( ). However this cannot be guaranteed a priori and is re ected in the statement of next result.
8
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
Lemma 3.3. Let be the solution of (2.3) with an arbitrary 2 L2(Q) and let be the solution of (2.4) with an arbitrary 2 H01 ( ). Then p ke (u)kL2 (Q) + keu (; T )kH ?1( ) 2C0ake0ukH ?1( ) 8 ? < 2 C0 aLjQj1=2 R?2 1=2 if 0 12 C0 +:
if 0 > C102 :
2C0R?2
Proof. We add twice (3.3) with (3.5), and argue as with (3.4) and (3.6). Since jR( )j 1 jR( )j ; jR( )j p 1 jR( )j ; 1 = 2 1 = 2 k kL2 (Q) 2 k kL2 (Q) krkL2( ) 2 k kL2(Q) we readily get p ? (3.7) 2 ? C01=2 ke (u)kL2 (Q) + keu(; T )kH ?1( ) 2C0ake0u kH ?1( ) + q(); where q() := q? ?1=2 + q+ 1=2 with q? := R?2 and q+ := C0aLjQj1=2 . We now observe that = 0 minimizes q(). If C001=2 1, then the rst assertion follows trivially from q(0 ) = 2(q? q+ )1=2 . Otherwise, if C001=2 > 1, then q+ < C02q? and q(1=C02) < 2C0q? . This concludes the argument. If U is a nite element solution, then the last term in Lemmas 3.1{3.3 (parabolic residual) can be further evaluated via Galerkin orthogonality; this is accomplished in x5. We stress that U need not be a discrete solution as the following application of Lemma 3.3 illustrates. Let U be the solution of (2.1) with " (s) := (s) + "s in place of (s), namely U is the usual vanishing viscosity approximation of u. The proof below is di erent from the original one in [10], easily extends to other perturbations of such as that in [10], and is valid under minimal regularity of both u0 and f which precludes compactness. Corollary 3.1. The following perturbation estimate is valid for " small ke (u)kL2 (Q) + keukL1 (0;T ;H ?1( )) C"1=2; where constant C depends on ku0kH ?1( ); j j; T; L; a; A, and kf kL1 (0;T ;H ?1( )). R Proof. Since R( ) = ?" Q U , the assertion results from Lemma 3.3 upon realizing that R?2 aC0"kU kL2 (Q) C". To derive a bound on kU kL2 (Q) under minimal regularity, we consider (2.1) with " instead of and strong data U 0 2 L2 ( ) and f 2 L2(0; T ; H ?1 ( )). We take = GU (; t) 2 H01 ( ), where the Green's operator G is the inverse of ? in H01 ( ), and obtain 1 (3.8) 2 dt hGU; U i + h " (U ); U i = hGf; U i: We then integrate (3.8) in time for 0 t s T , the maximum of kU (; t)kH ?1( ) being attained at t = s. We easily deduce
(3.9)
2 1 1 2 kU (; s)kH ?1 ( ) + A
Z
0
s
k " (U )(; t)k2L2 ( )dt Z
1 0 2 2 kU kH ?1 ( ) + kU (; s)kH ?1 ( ) 0
s
kGf (; t)kH01 ( )dt;
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE
9
which yields a bound for kU (; s)kH ?1 ( ) without use of Gronwall's inequality. We insert this bound in (3.9), now for arbitrary s, and use that "0 (s) 0(s) a > 0 and is bounded to arrive at the desired estimate for kU kL2 (Q) for strong data. We nally observe that the mapping (U 0 ; f ) ! (U; " (U )) from H ?1 ( ) L1(0; T ; H ?1 ( )) ! L1(0; T ; H ?1 ( )) L2 (Q) is Lipschitz, which is easily seen with the technique leading to (3.9). This allows us to regularize the weak data U 0 2 H ?1 ( ) and f 2 L1(0; T ; H ?1 ( )) via U0 2 L2( ) and f 2 L2 (0; T ; H ?1 ( )) and eventually pass to the limit in the regularization parameter . 4. Finite Element Discretization
We now introduce the fully discrete problem, which combines continuous piecewise linear nite elements in space with backward di erences in time.Pn We denote by n the time step at the n-th step and set tn := i=1 i. Let N be the total number of time steps, that is tN T . For any function 2 C 0((tn?1; tn ]; H ?1 ( )), we denote n() := (; tn). We denote by Mn a uniformly regular partition of into simplices [3]. Mesh Mn is obtained from Mn?1 by re ning/coarsening, and thus Mn and Mn?1 are compatible. Given a triangle S 2 Mn, hS stands for its diameter and S for its sphericity and they satisfy hS 2S = sin( S =2), where S is the minimum angle of S . Uniform regularity of the family of triangulations is equivalent to S > 0, with independent of n. We also denote by Bn the collection of interior interelement boundaries e of Mn in ; he stands for the size of e 2 Bn. the usual space of piecewise linear nite elements Let Vn H01( ) indicate n over Mn. Let fxnk gKk=1 denote the interior nodes of Mn. Let I n : C00( ) ! Vn be the usual Lagrange interpolation operator, namely (I n)(xnk ) = (xnk) for all 1 k K n. Finally, let the discrete inner product h; in be de ned by
h'; in :=
Z
1 I n(')dx = d+1
X
S 2Mn
jS j
X
xnk 2S
'(xnk )(xnk ) 8 '; 2 Vn:
This corresponds to the vertex quadrature rule, which can be easily evaluated element by element and leads to mass lumping [3]. We observe the validity of the local estimates k' ? I n(')kL1 (S) 81 h2S jr'jjrj 8 '; 2 Vn; kr'kL2(S) 13 ?S 1k'kL2(S) 8 ' 2 Vn: (4.1) R Consequently, the element inner product h'; inS := S I n(') satis es [15] k'kL2 (S) (h'; 'inS )1=2 =: k'knS C k'kL2(S) 8 ' 2 Vn; R and the discrepancy between h'; iS := S ' and h'; inS can be bounded by jh'; iS ? h'; inS j 18 h2S kr'kL2(S)krkL2(S) (4.2) 241 h2S ?S 1k'kL2(S)krkL2(S) 8 '; 2 Vn: The discrete initial enthalpy U 0 2 V0 is de ned at nodes x0k of M0 = M1 to be (4.3) U 0 (x0k ) := u0(x0k ) 8 x0k 2 nF0; U 0 (x0k ) := 0 8 x0k 2 F0: Hence, U 0 is easy to evaluate in practice. Then we set 0 := I 0 (U 0 ).
10
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
Discrete problem. Given U n?1 ; n?1 2 Vn?1, then Mn?1 and n?1 are modi ed as described below to get Mn and n and thereafter U n ; n 2 Vn computed according to n = I n (U n ) and
(4.4)
n hU 1
n ? I n U n?1 ; 'in + hrn ; r'i = hI n f n ; 'in
8 ' 2 Vn:
In view of the constitutive relation n = I n (U n ) being enforced only at the nodes, and the use of mass lumping, (4.4) is easy to implement and yields a monotone operator in RKn . This problem is solved below via an optimized nonlinear SOR [12]. However, these computational tricks introduce further consistency errors that are apparent from (3.2). Whether these devices preserve optimal accuracy is still to be explored. We nally introduce further notation. The following sets will be used later
F n := [fS 2 Mn : n(x) = 0 for some x 2 S g discrete free boundary; T n := [fS 2 Mn : n(x)n?1 (x) 0 for some x 2 S g transition region; C n := [fS 2 Mn : S is coarsened from Mn?1g coarsening set: We point out that the compatibility of Mn and Mn?1 yields I nU n?1 6= U n?1 only in C n, in which case Vn?1 \ Vn 6= Vn?1. Let the jump Jen of rn across e 2 Bn be
Jen := [ rn] e e = (rnjS1 ? rnjS2 ) e: If the unit normal vector e to e always points from S2 to S1, then Jen is well de ned. For any element S 2 Mn, JSn stands for the jumps of rn across @S n@ . Let U be the piecewise constant extension of fU ng de ned by U (; 0) = U 0 () and U (; t) := U n () for all tn?1 < t tn with n 1. Let n () ? I n U n?1 () U n?1 < t tn ; n 1; 8 t (4.5) Ut(; t) := n and let the interior residual Rn be
Rn := I nf n ? Ut (; tn): 5. A Posteriori Error Analysis
We now state and prove two rigorous a posteriori error estimates. Theorem 5.1 is based on H 1 regularity of and , solutions of the dual problems (2.3) and (2.4), and leads to Approach I. Exploiting H 2 space regularity of and in the spirit of Lemma 3.3 yields Theorem 5.2 and corresponding Approach II. The derivation below parallels, and in fact extends, that in x3 but exploits Galerkin orthogonality to express negative norms of the residuals R( ) and R( ) in terms of computable quantities.
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE 11
Theorem 5.1 (Approach I). There exist constants Ci > 0 depending pon the minimum angle of Mn and the space dimension d, and C := A1=2 = 2 + 1 and
C] := A1=2 C0a, such that the following a posteriori error estimates holds
ke (u)kL2(Q) + keukL1 (0;T ;H ?1( )) E I (u0; f; T; ; U; h; ) := C
X
i6=4
Ei + C] E4 ;
where the error indicators Ei = EiI for 0 i 8 are given by
E0 := ku0 ? U 0 kH ?1( )
N X 1=2 X I hekJenk2L2 (e) E1 := C1 n n=1 e2Bn N X 1=2 X I h2S kRn k2L2(S) E2 := C2 n n=1 S 2Mn N
X E3I := nkr (U n) ? rI n (U n)kL2 ( ) n=1 N 1=2 X nkU n ? I nU n?1 k2L2( ) E4 := n=1
E5 :=
N X
n=1
E6 := C6 E7 := E8 :=
kI nU n?1 ? U n?1 kH ?1( )
N X
n=1
N X
n
S 2Mn
h4S krRnk2L2 (S)
nkf n ? I nf n kH ?1( )
n=1 N Z tn X n=1
X
tn?1
kf ? f n kH ?1( )
initial error; jump residual; interior residual; constitutive relation; time residual; coarsening;
1=2
quadrature; interpolation; time discretization :
All indicators Ei can be evaluated explicitly in terms of the computed solution U , initial datum u0, and source term f . Indicators E0 ; E1I ; E2I ; E4 , and E5 are essential, and are also present for the heat equation but with di erent weights and cumulative e ect in time [5]. The error accumulation is measured here in L1 for E1I ; E2I ; E5 and in L2 for E4 , whereas it is in L1 for the heat equation; the latter exhibits a weaker dependence on T for T 1. The powers of meshsize in E1I and E2I are smaller than those for the heat equation, namely h3e and h4S , respectively, thereby re ecting the degenerate nature of (1.1), or equivalently the lack of H 2 space regularity of and . The remainding indicators E3I ; E6 ; E7 ; E8 are not essential and could in principle be removed at the expense of complicating the implementation of (4.4). In particular, we note that (U n ) 6= I n (U n ) only for S F n and njS 6= 0, provided is piecewise linear, and that E6 C E2I as results from (4.1).
12
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
The following theorem yields the same weights for E1I and E2I as the heat equation but yet with a worse error accumulation in time. Note that this improvement comes at the expense of a smaller outermost power, namely 1=4 instead of 1=2. 2 )1=2 (2A1=2 + 1=p2); Cz := Theorem 5.2 (Approach II) . Let C := ( A + 1 =C y 0 p 2C0a, II II II 1 := E1 L+jQE2j1=+2 E3 ;
and
E1II := C1 E2II := C2 E3II
N X n=1 N X n=1
N X
:=
n=1
n n
X
e2Bn
h3e kJenk2L2 (e)
X
S 2Mn
1=2
h4S kRn k2L2(S)
1=2
jump residual; interior residual;
1=2 n n n 2 constitutive relation: nk (U ) ? I (U )kL2 ( )
The following a posteriori error estimate holds
ke (u)kL2 (Q) + keukL1 (0;T ;H ?1( )) E II (u0; f; T; ; U; h; ) := Cz E0 + p
8
C102 :
8 X
i=5
Ei
5.1. Residuals. We rst express the residual R( ), for a generic function 2 H 1;2(Q), in terms of computable quantities. We notice that, since U is piecewise constant in time, summation by parts yields
N X n n n ? 1 N N 0 0 hU; @t i = hU ; ? i = hU ? i ? hU ; i ? hU n ? U n?1; n?1i: 0 n=1 n=1
Z
N X
T
Moreover, in light of (4.5), the last term can be written equivalently as N X
hU n ? U n?1 ; n?1i =
n=1
Z
0
T
hUt ; i
N X n ? 1 h U ; ? i + hI nU n?1 ? U n?1; n?1i: + t n=1 n=1 tn?1 N Z tn X
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE 13
Rearranging terms, the residual R( ) becomes
R( ) = (5.1)
N Z tn ? X n=1 tn?1 N Z tn X
hRn ; i ? hrI n (U n); r i
N X n ? 1 hUt ; ? i ? hI nU n?1 ? U n?1 ; n?1i ? n ? 1 n=1 n=1 t n Z N Z tn N t X X n f n ; i: n n n h f ? I hr I ( U ) ? r ( U ) ; r i + + n=1 tn?1 n=1 tn?1
We next rewrite the discrete problem (4.4) as follows
hRn; 'i ? hrI n (U n ); r'i = hRn ; 'i ? hRn ; 'in 8 ' 2 Vn; and subtract this expression from the right hand side of the residual R( ) given by (5.1). This crucial step is usually referred to as Galerkin orthogonality. We nally decompose the integral hrn; r( ? ')i over all elements S 2 Mn and next integrate by parts to obtain the equivalent expression
?hrI n (U n); r( ? ')i =
X
e2Bn
hhJen; ? 'iie 8 ' 2 Vn;
where hh; iie denotes the L2 -scalar product on e 2 Bn. Thus we easily arrive at N Z tn X n ; ? 'i n h R hh J ; ? ' ii + R( ) = e e n=1 tn?1 n=1 tn?1 e2Bn N Z tn N Z tn X X n?1 ? i n n n h U ; hr I ( U ) ? r ( U ) ; r i ? + t n=1 tn?1 n=1 tn?1 N Z tn ? N X X n ; 'i ? hRn ; 'in n n ? 1 n ? 1 n ? 1 h R ? hI U ? U ; i + n=1 tn?1 n=1 N Z tn X hf ? I nf n ; i =: I + + VII: + n ? 1 n=1 t N Z tn X
X
Since we need to approximate under minimal regularity assumptions, we resort to the Clement interpolation operator n : L2( ) ! Vn, which satis es the following local approximation properties [4], for all 2 H k ( ) and k = 1; 2, (5.2)
~ kS jjH k(S~); k ? nkL2(S) + hS kr( ? n)kL2 (S) Ch ~ ek?1=2jjH k(S~); k ? nkL2(e) Ch
where S~ is the union of all elements surrounding S 2 Mn or e 2 Bn. Constants C~ depend solely on the minimum angle of the mesh Mn. An important by-product of
14
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
uniform mesh regularity is that the number of adjacent simplices to a given element is bounded by a constant M independent of n, meshsizes, and time steps. Hence X kk2L2(S~) M kk2L2 ( ) 8 2 L2( ): (5.3) S 2Mn
This, in conjunction with (5.2) for k = 1, yields ~ 1=2 )krkL2( ) 8 2 H01 ( ): (5.4) krnkL2 ( ) (1 + CM We now estimate each term I to VII separately to derive Theorems 5.1 and 5.2. We argue with = , solution of (2.3) with 2 L2(Q). Similar estimates are valid for = , solution of (2.4) with 2 H01 ( ). We make extensive use of the a priori estimates of Lemma 2.1 (and Lemma 2.2). Selecting '(; t) = n (; t) for tn?1 < t tn and using (2.6), (5.2), and (5.3), terms I and II can be bounded as follows
jIj C1 jIIj C2
N Z tn X X tn?1
n=1 N Z tn X
e2Bn
hekJenk2L2(e)
X
n=1 tn?1 S 2Mn
1=2
h2 kRnk2 2 S
L (S )
kr (; t)kL2 ( ) p12 E1I k kL2 (Q);
1=2
kr (; t)kL2 ( ) p12 E2I k kL2 (Q):
Moreover, (2.6) also yields
jIIIj
N Z tn X n=1
tn?1
kr (U n) ? rI n (U n )kL2 ( )kr (; t)kL2 ( ) p12 E3I k kL2(Q):
Estimators E1I ; E2I , and E3I are those de ned in Theorem 5.1; constants C1 and C2 depend upon C~ in (5.2) and M in (5.3). Alternatively, if E1II ; E2II , and E3II are the estimators in Theorem 5.2, in light of (2.8) we can also write
jIj C1 jIIj C2
N Z tn X X tn?1
n=1 N Z tn X
n=1
e2Bn
X
tn?1
S 2Mn
h3e kJenk2L2 (e)
1=2
h4S kRn k2L2(S)
j (; t)jH 2( ) 21 ?1=2 E1II k kL2 (Q);
1=2
j (; t)jH 2( ) 21 ?1=2 E2II k kL2 (Q);
and, upon integration by parts,
jIIIj
N Z tn X
n=1
tn?1
k (U n ) ? I n (U n )kL2 ( )k (; t)kL2 ( ) 21 ?1=2 E3II k kL2 (Q):
R With the aid of (2.6) and the fact that ? n?1 = ttn?1 @s , we readily obtain
jIVj
N X
n=1
1=2 Z T n n n ? 1 2 k@t nkU ? I U kL2 ( )
(A + )1=2 E4 k kL2 (Q):
0
k2L2 ( )
1=2
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE 15
Using (2.6) again, now combined with (4.2) and (5.4), we easily deduce
jVIj C6
N Z X n=1
tn
tn?1
X
S 2Mn
h4S krRnk2L2(S)
1=2
kr (; t)kL2 ( ) p12 E6k kL2 (Q);
where C6 depends on C~ and M . The remaining bounds follow directly from (2.6)
jVj jVIIj
N X
kI nU n?1 ? U n?1kH ?1( )kr
n=1 N X Z tn
n=1
tn?1
n?1 k
L2 ( ) p12 E5 k kL2 (Q) ;
kf ? I nf n kH ?1( )kr (; t)kL2 ( ) p12 (E7 + E8)k kL2 (Q):
Collecting all the previous estimates we get the following bound for R( ) (5.5)
8 X
jR( )j p1 E + (A + )1=2 E + 4 k kL2(Q) 2 i=5 i
(
p12 (E1I + E2I + E3I ) 1 ?1=2 ?E II + E II + E II : 1 2 3 2
Similarly, using (2.7) we obtain (
(5.6)
8 I I I ? jR( )j X 1=2 E + E1 + E2 +? E3 1 p A + ) E + 4 p12 ?1=2 E1II + E2II + E3II : krkL2( ) i=5 i 2
5.2. Proofs of Theorems 5.1 and 5.2. Upon adding A1=2 times (5.5) with (5.6) and taking ! 0, Theorem 5.1 is an easy consequence of Lemmas 3.1 and 3.2, because E I is increasing in T and tN T . In order to prove Theorem 5.2 we proceed as in Lemma 3.3. We take 1=C02 and resort to (3.7), which in the present context becomes (5.7)
8
X ke (u)kL2 (Q) + keNu kH ?1( ) Cz E0 + Ei + CyE4 + p12 Czq() i=5
with
q() := q? ?1=2 + q+ 1=2; q? := E1II + E2II + E3II ; q+ := LjQj1=2 : We then realize that = 1 optimizes q(), and argue as in Lemma 3.3. 5.3. Approach III. The local treatment of terms I, II, and III above, in conjunction with (2.9), suggests removing the factor ?1=2 away from discrete interfaces. This idea, heuristic in the sense that L2t Hx2 regularity of and is not necessarily uniform, yields a method somewhat in between Approaches I and II. We set the estimators E1III ; E2III ; E3III equal to those of Approach I near discrete interfaces and
16
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
equal to those of Approach II (heat estimators), but with factor D=C instead of 1 away from discrete interfaces, where D := a1=2 C0, and obtain ke (u)kL2 (Q) + keukL1(0;T ;H ?1( )) E III (u0; f; T; ; U; h; ) 8 X
:= C E0 +
i=5
?
Ei + C] E4 + C E1III + E2III + E3III :
To this end, we add A1=2 times (5.5) with (5.6), let ! 0, replace kD2 kL2(Q0 ) by k kL2(Q0) for = ; and use (2.9), and ultimately pretend that the lower bound of b is a instead of a0 = 0. If and were computable, we could keep these Hx2 norms and evaluate them locally, in which case they could be viewed as weights. 6. Adaptive Algorithm
The error estimators E (u0 ; f; T; ; U; h; ) of x5 entail an L1 or L2 norm in time, which is impractical in that the entire evolution history would be needed to control E . We resort here to an L1 norm in time, and explain our error equidistribution strategy, which is based upon minimizing the spatial degrees of freedom for a uniform error distribution in time. A similar strategy is derived in [1] for a linear elliptic problem. We then fully discuss element error indicators, and all tests necessary for mesh and time step admissibility. 6.1. Mesh Design. Let M(t) be a time-dependent mesh with variable meshsize h(x; t) and let (t) be the underlying variable time-step. Since hd is proportional to the volume of a generic space-time nite element, then the computational complexity of (4.4) can be measured by the total number of degrees of freedom (6.1)
M :=
Z
Q
where 0 < ? (x; t) + and
Z
M (t) :=
(6.2)
T
Z
h(x; t)?d (t)?1 (x; t) dxdt =
0
(t)?1 M (t) dt;
h(x; t)?d (x; t) dx
stands for the cardinality of M(t). Let E be a generic a priori error of the form (6.3)
E=
Z
0
T
E (t)dt =
Z
Q
?
h(x; t) + (t) E (x; t) dxdt:
Given an error tolerance ", we then pose the following question: optimize h and for an error distribution (6.4) E = ": We have to minimize (6.1) subject to the constraint (6.4). This constrained optimization problem is equivalent to seeking a saddle point of the functional
L(h; ; ) :=
Z
Q
h?d ?1 ? " ?
Z
Q
(h + )E ;
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE 17
with Lagrange multiplier . Di erentiation with respect to h and yield 0=
Z
Q
0=
@h L(h; ; ) =
Z
Q
Z
0
@ L(h; ; ) =
TZ Z
T
0
(?dh?(d+1) ?1 + h ?1E ); Z
(?h?d ?2 + ?1E );
where = (x; t) and = (t) are smooth functions. Hence Z d 1 ?1M (t): ? d ? 1 (6.5) h E = h ; E =
These relations, in conjunction with (6.1) and (6.3), yield Z Z T d + M: 1 ? d ? 1 " = h + ?1M (t) = d 0 Q We thus deduce the following expression for the Lagrange multiplier = d + M" : We conclude from (6.5) that the optimization procedure consists of equidistributing the local space-time errors according to the recipes +1 Z d M ? 1 d + E = + d ": (6.6) M h E = + d "; M (t)
A serious diculty of (6.6) is that M is never available for all times so as to make the optimization feasible. A di erent and more stringent objective would be to equidistribute pointwise discretization errors in (6.3). This leads to the local requirements
T j jh E = 1 ";
(6.7)
2
T
Z
E = 21 ";
which corresponds to having a constant integrant in (6.3) and yields (6.4). This method does not require M , and can thus be implemented. If we further impose = T=N to be constant, we end up with the method chosen in [11] for a priori mesh design. Upon combining (6.6) and (6.7), we get a third strategy that can still be implemented: we optimize the spatial degrees of freedom for a uniform error distribution E (t) in time. Minimizing (6.2) for each 0 < t T , subject to the constraints Z
h(x; t) E (x; t) dx = 1 ";
T (t)
Z
T E (x; t) dx = 12 "; 2
results in the following restrictions on h and , which are intermediate between (6.6) and (6.7), Z ? 1 d + 1 (6.9) TM (t)(x; t) h(x; t) E (x; t) = 2 "; T (t) E (x; t) dx = 12 ":
(6.8)
We adopt such a viewpoint here, and discuss its implementation next. Note rst that E would depend on the discrete solution, and so implicitly on h and , if E in (6.3) were an a posteriori error estimator. So the above analysis is a priori.
18
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
6.2. Equidistribution Strategy. The equidistribution strategy (ES) is an iterative procedure that in the kth step improves upon a mesh density hk (x; t) by means of the following two opposite operations, based on the rst constraint in (6.9). Let R ?d 0 < < 1 be a re nement factor and Mk = hk from (6.2) with h = hk . (a) Re nement. Set hk+1 = hk for all elements satisfying (6.10)
k := 2TMk ?1 hkd+ E > ":
(b) Coarsening. Set hk+1 = ?1hk for all elements verifying
k < 2d+ ":
(6.11)
In practice, more elements will be re ned and fewer will be coarsened to preserve mesh conformity. In fact, every local mesh modi cation involves also adjacent elements. Since Z
Z
d ?d h?d = ?d Mk ; Mk+1 = h?k+1 k
we see that k+1 k for re ned elements whereas k+1 ?2d? k < " for coarsened elements. We conclude that coarsened elements will not be candidates for re nement in the (k +1)th step, because k+1 < " and so test (a) fails, and that ?
k+1 min "; k : The latter guarantees convergence of ES in a nite number of steps. 6.3. Element Indicators. We now introduce error indicators of the form (6.8) with density E (x; t) computable element by element and denoted by E n(S ) for S 2 Mn and t = tn. We rst replace the global norm H ?1 ( ) in E0 ; E5 ; E7 ; E8 by the L2 ( ) norm scaled with the smallest diameter d of , which is the Poincare constant. We next show how to determine computable global indicators E0 (initial error), E ;1 ; E ;2 (temporal errors) and Eh;1; Eh;k 2 (spatial errors) satisfying (6.12)
E0 E0; E4 E ;1; E8 E ;2; E5 + E6 + E7 Eh;1 ; E1k + E2k + E3k Eh;k 2 ;
where k = I; II stands for the approach. Each of these errors Ep is obtained from element indicators Epn(S ) via
Ep = 1max nN
X
S 2Mn
Epn(S )
1=2
;
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE 19
where the max is super uous for p = 0. The element indicators are given by
E0(S ) := d2 jS jku0 ? U 0 k2L1 (S) E ;n 1(S ) := T kU n ? I nU n?1 k2L2(S) E ;n 2(S ) := T 2d2 kf ? f nk2L1 (tn?1;tn;L2 (S)) Eh;n 1(S ) := 3T 2 d2 n?2kI nU n?1 ? U n?1k2L2 (S) + C62h4S krRnk2L2 (S) + d2 kf n ? I nf n k2L2(S) I ;n 2 Eh;2(S ) := 3T 21 C12hS kJSnk2L2 (@S) + C22h2S kRnk2L2(S) + kr (U n) ? rI n (U n )k2L2 (S) Eh;II;n2 (S ) := 3T 21 C12h3S kJSnk2L2(@S) + C22h4S kRnk2L2(S) n n n 2 + k (U ) ? I (U )kL2 (S)
local initial error; local time residual; local time discretization; local coarsening local quadrature local interpolation; local jump residual I local interior residual I local constitutive relation I; local jump residual II local interior residual II local constitutive relation II:
Since u0 2 W 1;1( nF0), interpolation gives ku0 ? U 0 kL1 (S) ChS kru0kL1 (S) = O(hS ) away from F0 . If S intersects F0, instead, then u0 has a jump of size L which yields ku0 ? U 0 kL1 (S) = O(L). In both instances, the ensuing orders match that of the local jump residual I, which indicates that M0 will not be drastically modi ed. If in addition u0 2 W 2;1( nF0), then ku0?U 0 kL1(S) Ch2S kD2 u0kL1 (S) = O(h2S ) d away from F0 , which yields the order O(h4+ S ) of the local jump residual II. According to x5.3, we introduce the local indicators for heuristic Approach III which is based on combining Eh;I;n2(S ) and Eh;II;n2 (S ). With D = a1=2 C0 we set
Eh;III2;n(S ) :=
(
Eh;I;n2(S ) if S T n D22 E II;n (S ) if S 6 T n : C h;2
A simple but tedious calculation shows that (6.12) is ful lled, and the following inequality is valid for Approaches k = I and III (heuristic for Approach III) (6.13)
E k (u0; f; T; ; U; h; ) E k := C (E0 + E ;2 + Eh;1 + Eh;k 2 ) + C]E ;1:
On the other hand, if 1=C02 and q() := Eh;II 2?1=2 + LT 1=2j j1=2 1=2, then (5.7) and (6.12) lead to (6.14) E II (u0 ; f; T; ; U; h; ) E II := Cz(E0 + E ;2 + Eh;1) + CyE ;1 + p12 Czq():
20
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
6.4. Adaptive Strategy. Given an error tolerance ", the objective is to adap-
tively select time steps and mesh densities h so that the spacial degrees of freedom are optimized for a uniform error distribution in time and (6.15) E k(u0; f; T; ; U; h; ) ": To this end we equidistribute the local estimators of x6.3 according to (6.9) via ES, and use bisection to perform re nement/coarsening operations (6.10) and (6.11). Bisection creates compatible consecutive meshes, extends naturally from 2D to 3D, and is handy for combined re nement/coarsening operations; we refer to [2]. Given re nement parameters ? > 0 and coarsening parameters > 0 satisfying ?0 + ? + ?h 1; < ? ; h < ?h ; then time steps and mesh densities for Approaches k = I and III are reduced until CE0 ?0 "; 1 1 (6.16) 2 " C] E ;1 ; C E ;2 2 ? "; 1 h " C Eh;1 ; C E k 1 ?h ": h;2 2 2 In fact, from (6.13) we readily obtain (6.15). Achieving (6.15) for Approach II is more problematic in that the optimal choice of in (6.14) turns out to be Eh;II 2 2 := LT 1=2j j1=2 ; which is not a computable quantity as it involves a maximum over all time steps. To overcome this diculty we restrict the meshsize and time step as follows CzE0 ?0 "; 1 " CyE ;1 ; CzE ;2 1 ? "; 2 2 (6.17) 1 h " Cz Eh;1 1 ?h "; 2 2 1 G h " C 2aE II 1 G?h "; 0 h;2 4 4 with G := min 1; 4aLT?1h=2"j j1=2 : We claim this yields (6.15). In fact, if G < 1 then 2 2 Eh;II 2 2 = LT 1=2j j1=2 16C 2?ah2 L" 2T j j = CG022 < C102 ; 0 whence II )1=2 1 ?h ": p12 Czq (2 ) = 2C0 aL1=2 T 1=4 j j1=4 (Eh; 2 2 2 On the other hand, if G = 1, then with = 1=C0 in (6.14) we deduce ? II + aLT 1=2 j j1=2 1 ? " + 1 ? " = 1 ? ": p12 Cz q C12 = C02 aEh; 2 4 h 4 h 2 h 0 We now describe an algorithm based on (6.16) for Approach I (and III). The implementation of Approach II is similar, but based on (6.17) instead. Let M n denote the cardinality of Mn at any step of ES.
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE 21
6.4.1. Initial Mesh. Given a coarse mesh, ES bisects all S 2 M0 such that 2 "2 ? 0 2 C E0(S ) > M 0 ;
as suggested by (6.10) without T . We interpolate u0 by means of (4.3) to nd U 0 (and 0) and estimate E0(S ) as E0(S ) h2S jS jkrU 0k2L1 (S). 6.4.2. Time Step Selection. Starting with n = n?1, the algorithm checks whether
E n := 4 max C]2
X
S 2Mn
E ;n 1(S ); C2
X
S 2Mn
E ;n 2(S ) > ? 2 "2 ; E n < 2 "2 :
In the rst case n is reduced, whereas in the second one (corresponding to n being too small) n is accepted but the initial guess for the next time step size is enlarged. 6.4.3. Mesh Size Selection. Starting from Mn = Mn?1, ES checks whether 2 2 2 2 ? EhI;n(S ) := 4C2 max Eh;n 1 (S ); Eh;I;n2(S ) > ?Mh "n ; EhI;n(S ) < Mh "n ;
as suggested by (6.10) and (6.11). Then re nement and coarsening operations are performed accordingly, with the precaution of choosing h ?h properly to prevent ES from alternating such operations over the same elements. 6.4.4. Flow Chart. Figure 6.1 shows a ow diagram of the adaptive algorithm for Approach I. After each iteration of ES, both U n and n are recalculated on the new mesh using the new time step size. To minimize the overall computational cost, a compromise is reached between the minimization of the number of mesh elements (degrees of freedom) and iterations of ES. ES stops iterating as soon as X
S 2Mn
EhI;n(S ) ?h2 "2
is ful lled, thereby allowing some elements to violate the local error tolerance. Consequently, discretization errors might not be equidistributed. As implemented, and shown in the ow diagram, at least one mesh modi cation per time step is performed in order to permit elements near the moving interface to be re ned, even if the global error bound is already ful lled by the old mesh. 6.5. Convergence. We elaborate here on convergence of the algorithm of x6.4; a more comprehensive study will be given elsewhere [13]. For all 1 n N , we assume the following discrete a priori estimates, which are slightly stronger than the natural bounds
krnkL2 ( ) = O(1); jF nnZ nj = o(1); X kU n ? I nU n?1k2L2 ( ) + n?1kn ? I nn?1k2L2 ( ) + hS kJSnk2L2(@S) = o(1); S nF n
22
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
Figure 6.1. Flow diagram of the adaptive Approach I. where Z n := [fS 2 Mn : njS = 0g is the numerical mush. We now examine a few
relevant terms, using the hyperbolic and parabolic relations h2S C ; rn := max h2S C : (6.18) rHn := Smax P H P S nT n n T n n2 We consider rst the jump residual, for which it suces to estimate the contribution
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE 23
of F nnZ n, namely X
eF nnZ n
hekJenk2L2 (e) C krnk2L2 (F nnZ n) = o(1):
The same calculation applies to the constitutive relation. For the interior residual, we decompose the integral over T n and the complement, where we use U n ? I nU n?1 = n ? I nn?1, to deduce
h2S kU n ? I nU n?1 k2 rn kU n ? I nU n?1 k2 H L2 (S ) L2 ( ) 2 S 2Mn n n + r P kn ? I nn?1k2L2( ) = o(1): n X
Finally, kU n ? I nU n?1 kL2( ) = o(1) also controls the time residual. This demonstrates that the a posteriori estimators of x6.4 tend to 0 as both meshsize and time-step approach 0. Therefore the goal (6.15) is achievable. On the other hand, (6.18) re ects the hyperbolic structure of the interface as well as the parabolic structure of the problem elsewhere. We refer to [11] for a similar observation. 7. Simulations
We present results of the adaptive method for the model problem from [12] with known exact solution (
?
0:75 6r2 ? 61 if r < 61 ; u(x; y; t) := ? 1 + 1:5 ? 0 (t) y? r (t) (r ? 16 ) if r 61 ;
(7.1)
where r := (x ? 13 )2 + (y ? (t))2 1=2. The interface is a circle of constant radius r = 61 centered at ( 13 ; (t)) with (t) := 0:5 + 0:1 sin(12:5t), and thus oscillates in the vertical direction with highly varying normal velocity. Initial datum, boundary values, and heat source f are directly computed from (7.1). The domain is = (0:1; 0:7) (0:2; 0:8) and the nal time is T = 1. Note that f is discontinuous across the true interface. Therefore a simple minded use of I nf leads to large errors which dominate the local indicators. However a discontinuous f is not realistic in practice, being just the cost of having a simple exact solution. To examine the essense of the proposed methodology, a special quadrature is used for those triangles crossed by the exact free boundary: their intersection is determined rst and then separate quadrature used in the resulting quadrilateral and triangle. The nonhomogeneous Dirichlet boundary condition is handled in the standard fashion since the interface does not touch the xed boundary. Several simulations were carried out with error tolerances " = 20; 14; 10; 7; 5. The following parameters were used for partitioning the total error into initial, temporal, and spatial components ?
?0 = ? = 0:2; ?h = 0:6; = 0:155; h = 0:2683 (I); 0:1897 (II, III):
24
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
Figure 7.1. Approaches I with " = 14(10) and II and III with " =
10(7); Triangle count: I 1016 (2747), II 700 (1592), III 786 (2129)
Figure 7.2. Meshes and isotherms for Approach II with " = 7
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE 25
Figure p p 7.3. E ect of coarsening for Approach I and " = 10, h =?h =
0:2( 0:1).
Figure 7.4. Total element count vs true error for Approaches I, II,
and III.
Since the expected local coarsening error is proportional to the power of meshsize, and such a power is smaller for Approach I than for II and III, more coarsening (larger h) is allowed for Approach I. The results are summarized in Figures 7.1 to 7.4. Figure 7.1 shows the meshes at time t = 0:3 produced by Approaches I with " = 14 (10), and II and III with " = 10 (7). Meshes from left to right correspond to Approaches I, II, and III. The true errors for all three pictures on top and botton are similar. Approach I leads to more triangles than II and III because of the lower power of meshsize in the estimators. Approaches II and III are comparable. Figure 7.2 shows three meshes generated by Approach II for quite di erent interface velocities. In the left and middle pictures the velocity is highest, and the
26
R.H. NOCHETTO, A. SCHMIDT, AND C. VERDI
interface is moving upwards (t = 0:5) and downwards (t = 0:765) as re ected in the isotherms (the interface is enhanced with a bold line). The rightmost picture corresponds to vanishing velocity and so to circular isotherms (t = 0:876). Higher interface velocities lead to higher re nement near the free boundary, because velocity is proportional to the temperature gradient jump. Higher velocities yield also larger temperature gradients in the liquid phase and corresponding additional mesh re nement. The triangle counts are 1670, 1734, and 1366, respectively. Figure 7.3 illustrates the e ect of h for Approach I with the interface moving downwards. The ratio ?h = h cannot be close to 1 nor too small. The former situation leads to mesh oscillations due to repeated coarsening/re nement operations over the same elements. A small ratio instead leads to meshes with unnecessary triangles, which thereby re ect the evolution history. This is an undesirable event. In Figure 7.4 we provide a quantitative comparison: a plot of the total element count (added over all time steps), which is proportional to the total amount of work, as a function of the true error. Approach II performs slightly better in the range tested. References 1. I. Babuska and W.C. Rheinboldt, Error estimates for adaptive nite element computations, SIAM J. Numer. Anal. 15 (1978), 736{754. 2. E. Bansch, Local mesh re nement in 2 and 3 dimensions, IMPACT Comput. Sci. Engrg. 3 (1991), 181{191. 3. P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. 4. Ph. Clement, Approximation by nite element functions using local regularization, RAIRO Anal. Numer. 9 (1975), 77{84. 5. K. Eriksson and C. Johnson, Adaptive nite element methods for parabolic problems I: a linear model problem, SIAM J. Numer. Anal. 28 (1991), 43{77. 6. , Adaptive nite element methods for parabolic problems II: optimal error estimates in L1 L2 and L1 L1 , SIAM J. Numer. Anal. 32 (1995), 706{740. 7. , Adaptive nite element methods for parabolic problems IV: nonlinear problems, SIAM J. Numer. Anal. 32 (1995), 1750{1763. 8. A. Friedman, Variational Principles and Free Boundary Problems, Wiley, New York, 1982. 9. O.A. Ladyzenskaja, V. Solonnikov, and N. Ural'ceva, Linear and Quasilinear Equations of Parabolic Type, vol. TMM 23, AMS, Providence, 1968. 10. R.H. Nochetto, Error estimates for multidimensional singular parabolic problems, Japan J. Indust. Appl. Math. 4 (1987), 111{138. 11. R.H. Nochetto, M. Paolini, and C. Verdi, An adaptive nite elements method for two-phase Stefan problems in two space dimensions. Part I: stability and error estimates. Supplement, Math. Comp. 57 (1991), 73{108, S1{S11. 12. , An adaptive nite elements method for two-phase Stefan problems in two space dimensions. Part II: implementation and numerical experiments, SIAM J. Sci. Statist. Comput. 12 (1991), 1207{1244. 13. R.H. Nochetto, A. Schmidt, and C. Verdi, Mesh and time step modi cation for degenerate parabolic problems, in preparation. , Adaptive algorithm and simulations for Stefan problems in two and three dimensions, 14. in preparation. 15. P.A. Raviart, The use of numerical integration in nite element methods for solving parabolic equations, Topics in Numerical Analysis (J.J.H. Miller ed.), Academic Press, London, 1973, pp. 233{264. Department of Mathematics and Institute for Physical Science and Technology, University of Maryland, College Park, MD 20742, USA.
A POSTERIORI ERROR ESTIMATION FOR DEGENERATE PARABOLIC PDE 27 E-mail address: [email protected] Institut fur Angewandte Mathematik, Universitat Freiburg, 79106 Freiburg, Germany.
E-mail address: [email protected]
Dipartimento di Matematica, Universita di Milano, 20133 Milano, Italy.
E-mail address: [email protected]
