NEW COMPUTATIONAL TECHNOLOGIES for RELIABLE COMPUTER SIMULATIONS - - PowerPoint PPT Presentation

✬ ✫ ✩ ✪

NEW COMPUTATIONAL TECHNOLOGIES for RELIABLE COMPUTER SIMULATIONS Sergey Korotov⋆,†

⋆ Institute of Mathematics

Helsinki University of Technology, Finland

† Academy of Finland

1

✬ ✫ ✩ ✪

Main Problem in Mathematical Modelling and Numerical Analysis: RELIABLE VERIFICATION OF ACCURACY OF APPROXIMATE SOLUTIONS OBTAINED IN COMPUTER SIMULATIONS

• What are the sources of various errors affecting the

reliability of numerical solutions then ?

2

✬ ✫ ✩ ✪ U Physical process ⇓ ε1 − →

Modelling error

⇓ u Differential model Au = f ⇓ ε2 − →

Approximation error

⇓ uh Discrete model Ahuh = fh ⇓ ε3 − →

Computational error

⇓ uǫ

h

Numerical solution Aǫ

huǫ h = f ǫ h 3

✬ ✫ ✩ ✪

Two Principal Relations

• Computations on the base of a reliable (certified) model.

Here ε1 is assumed to be small and computed uǫ

h gives a

desired information on U U − uǫ

h ≤

ε1 + ε2 + ε3 (1)

• Verification of a mathematical model. Here physical

data U and uǫ

h are compared to judge on the quality of

mathematical model ε1 ≤ U − uǫ

h + ε2 + ε3

(2)

4

✬ ✫ ✩ ✪ Thus, two major problems of mathematical modelling:

• reliable computer simulation
• verification of mathematical models by comparing

physical and mathematical experiments require efficient methods able to provide COMPUTABLE AND REALISTIC estimates of ε2 + ε3

5

✬ ✫ ✩ ✪

No Error Control - No Reliability

0.2 0.4 0.6 0.8 1 0.5 1 −0.5 0.5 1 1.5 2

exact solution

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

mesh with 741 nodes error = 0.09212

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

mesh with 785 nodes error = 0.11562

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

mesh with 855 nodes error = 0.44662

6

✬ ✫ ✩ ✪

NEW TECHNOLOGIES OF ERROR CONTROL Elliptic Type BVPs

Model elliptic BVP: Find a function u such that

−div(A∇u) = f in Ω & u = 0

• n

∂Ω (3)

Ω is bounded domain in Rd with Lipschitz boundary ∂Ω, f ∈ L2(Ω), matrix of coefficients A is symmetric, with entries aij ∈ L∞(Ω), i, j = 1, . . . , d, and is such that c2|ξ|2 ≥ A(x)ξ · ξ ≥ c1|ξ|2 ∀ξ ∈ Rd ∀x ∈ Ω (4)

7

✬ ✫ ✩ ✪ Weak formulation: Find u ∈ H1

0(Ω) such that

• Ω

A∇u · ∇w dx =

• Ω

fw dx ∀w ∈ H1

0(Ω)

(5)

• Let ¯

u be any function from H1

0(Ω) (e.g., computed by some

numerical method) considered as approximation of u

• We are interested in diverse and reliable control of deviation (or

error) e := u − ¯ u during computational process, depending on the final goal(s) of calculations

• So far two possibilities for the error control are well realized and

used in mathematical and engineering communities

8

✬ ✫ ✩ ✪ Control of error in global norm: We are usually interested in estimation of gradient of deviation in energy norm | | |∇(u − ¯ u)| | | :=  

• Ω

A∇(u − ¯ u) · ∇(u − ¯ u) dx  

1/2

(6)

• Such estimation gives a general presentation on the quality of ¯

u

• Such type of control is in focus since [Babuˇ

ska, Rheinboldt, 1978] Control of local errors: Another trend in a posteriori error estimation is based on concept of local error control in addition to classical control in energy norm

• This approach is motivated by practical needs: analysts are often

interested not only in the value of the overall error, but also in controlling errors over certain parts of solution domain, or relative to some interesting characteristics (“quantities of interest”)

9

✬ ✫ ✩ ✪

• Common way for performing such type control is to introduce

linear functional ℓ associated with subdomain of interest (and/or with a “quantity of interest”) and to obtain a computable estimate for ℓ(u − ¯ u)

• For example, one can be interested in estimation of

ℓ(e) = ℓ(u − ¯ u) =

• Ω

ϕ(u − ¯ u) dx (7) where ϕ ∈ L2(Ω) and supp ϕ = ω ⊆ Ω, which provides us with info

• n error e locally, in subdomain ω
• Estimates for ℓ(u − ¯

u) can also be used for estimation of unknown “quantity of interest” ℓ(u) since ℓ(u) = ℓ(¯ u) + ℓ(u − ¯ u)

• Certain “quantities of interest” (e.g., in linear elasticity) are often

more interesting to practitioners than solution itself

10

✬ ✫ ✩ ✪

Error Control in Energy Norm

Upper estimate: Since u − ¯ u ∈ H1

0(Ω) we have

| | |∇(u − ¯ u)| | |2 =

• Ω

f(u − ¯ u) dx −

• Ω

A∇¯ u · ∇ (u − ¯ u) dx = =

• Ω

f(u− ¯ u) dx−

• Ω

(A∇¯ u−y∗)·∇(u− ¯ u) dx−

• Ω

y∗ ·∇(u− ¯ u) dx = =

• Ω

(f + div y∗)(u − ¯ u) dx −

• Ω

A(∇¯ u − A−1y∗) · ∇(u − ¯ u) dx ≤ ≤ f + div y∗ u − ¯ u + | | |∇¯ u − A−1y∗| | | | | |∇(u − ¯ u)| | | where y∗ ∈ H (div; Ω) and · is a standard norm in L2(Ω)

• In above Green formula and CBS inequality are used

11

✬ ✫ ✩ ✪ From Friedrichs inequality w ≤ cΩ∇w ∀w ∈ H1

0(Ω)

(8) and conditions on A we observe that u − ¯ u2 ≤ c2

Ω∇ (u − ¯

u) 2 ≤ c2

Ω

c1 | | |∇(u − ¯ u)| | |2 (9) Thus, we immediately get | | |∇(u − ¯ u)| | | ≤ cΩ √c1 f + div y∗ + | | |∇¯ u − A−1y∗| | | (10) After squaring and using Young inequality, we get | | |∇(u − ¯ u)| | |2 ≤

• 1 + 1

β c2

Ω

c1 f + div y∗2 + (1 + β) | | |∇¯ u − A−1y∗| | |2 (11) valid for any y∗ ∈ H(div; Ω) and any β > 0

12

✬ ✫ ✩ ✪

• Upper bound (11) was first obtained in [Repin, 1997] in a general

form using quite complicated tools of duality theory

• More general problems with mixed Dirichlet/Neumann/Robin

boundary conditions are treated in a just presented (and more simple) way in [Korotov, 2005], where (11) follows as a particular case Lower estimate: It is easy to show that | | |∇(u − v)| | |2 = 2(J(v) − J(u)) ∀v ∈ H1

0(Ω)

(12) where J(w) = 1

2

• Ω

A∇w · ∇w dx −

• Ω

fw dx (energy functional)

• Function u minimizes J =

⇒ J(u) ≤ J(w) ∀w ∈ H1

0(Ω)

= ⇒ | | |∇e| | |2 ≥ 2(J(¯ u) − J(w)) (13) where w is any function from H1

0(Ω) 13

✬ ✫ ✩ ✪

Practical Realization

• Assume that computations are performed on Th1, Th2, Th3, . . . ,

where h = h1 > h2 > h3 > . . . , and, thus, we always have several successive approximations uh1, uh2, uh3, . . . On computation of constant cΩ:

• Upper bound contains the only unknown constant cΩ
• cΩ =

1 √λΩ , where λΩ is smallest eigenvalue of Laplacian for Ω

• Other estimation techniques involve many unknown constants:
• which are very hard to estimate
• which have to be always recomputed if computational mesh has

changed

• But cΩ remains the same under any changes of meshes

14

✬ ✫ ✩ ✪ On minimization of upper bound: now ¯ u ≡ uh M ⊕(uh, β, y∗) :=

• 1 + 1

β c2

Ω

c1 f+div y∗2+(1 + β) | | |∇uh−A−1y∗| | |2

• Coarse upper bounds can be computed fast using, e.g.,

y∗ = Gµ(∇uµ) ∈ H(div, Ω), µ = h1, h2, h3, . . . , and Gµ is some gradient averaging operator

• Sharp estimates require a real minimization of M ⊕(uh, β, y∗)

with respect to “free variables” y∗ and β On computation of lower bound: We usually try to have J(uh1) > J(uh2) > J(uh3) > . . . , which suggests lower bounds | | |∇e| | |2 = | | |∇(u − uh)| | |2 ≥ 2(J(uh) − J(uµ)) > 0 (14) where µ = h2, h3, . . .

15

✬ ✫ ✩ ✪

Test No. 1

−1 1 −1 1 −1 −0.5 0.5 1 −1 −0.5 0.5 1 0.5 1 1.5 2 1 2 3 4 1 1.5 1.7514 2 2.5 3 3.5 4 4.5 5

• A ≡ I, f ≡ 10, Th with 92 nodes
• Several successive meshes are used in computations of bounds
• Upper bound decreases from 4.6426 to 2.4443
• Lower bound grows from 1.2191 to 1.7469
• |

| |∇e| | |2 ≈ 1.7514

16

✬ ✫ ✩ ✪

Error Control Via Linear Functionals

• Second technology is developed for controlling local errors

Consider e = u − ¯ u measured in terms of linear functional ℓ Example: ℓ(e) =

• Ω

ϕ (u − ¯ u) dx (15) where ϕ ∈ L2(Ω) and supp ϕ = ω ⊆ Ω Adjoint problem: Find v ∈ H1

0(Ω) such that

• Ω

A∇v · ∇w dx = ℓ(w) ∀w ∈ H1

0(Ω)

(16)

• Adjoint problem is uniquely solvable

17

✬ ✫ ✩ ✪

• Usually we only have some approximation ¯

v (e.g., computed by FEM on adjoint mesh Tτ)

• It can be shown that

ℓ(u − ¯ u) = E0(¯ u, ¯ v) + E1(e, eϕ) (17) where E0(¯ u, ¯ v) =

• Ω

f¯ v dx −

• Ω

A∇¯ v · ∇¯ u dx (18) and E1(e, eϕ) =

• Ω

A∇e · ∇eϕ dx (19)

• Deviation eϕ := v − ¯

v is computational error in adjoint problem

18

✬ ✫ ✩ ✪ Proof: From integral equality (16), we have (w := u − ¯ u) ℓ(u − ¯ u) =

• Ω

A∇v · ∇(u − ¯ u) dx = =

• Ω

A∇(v − ¯ v) · ∇(u − ¯ u) dx +

• Ω

A∇¯ v · ∇(u − ¯ u) dx = = E1(e, eϕ) −

• Ω

A∇¯ v · ∇¯ u dx +

• Ω

A∇u · ∇¯ v dx = = E1(e, eϕ) −

• Ω

A∇¯ v · ∇¯ u dx +

• Ω

f¯ v dx = = E0(¯ u, ¯ v) + E1(e, eϕ)

19

✬ ✫ ✩ ✪

• First term E0 is directly computable once ¯

u and ¯ v are computed, but term E1 contains unknown gradients ∇u and ∇v

• In order to estimate E1, we use the following identity

2E1(e, eϕ) = | | |∇(αe + 1 αeϕ)| | |2 − α2| | |∇e| | |2 − 1 α2 | | |∇eϕ| | |2 (20) valid for any α > 0

• Above identity contains errors in energy norm for both, primal

and adjoint problems, and we can use corresponding two-sided estimates M ⊖

f ≤ |

| |∇e| | |2 ≤ M ⊕

f

& M ⊖

ϕ ≤ |

| |∇eϕ| | |2 ≤ M ⊕

ϕ

(21)

20

✬ ✫ ✩ ✪

• For the first term in RHS of (20), we observe that

| | |∇(αe + 1 αeϕ)| | |2 = | | |∇(αu + 1 αv) − ∇(α¯ u + 1 α ¯ v)| | |2 (22)

• Function αu + 1

αv can be perceived as the solution of problem:

Find uα ∈ H1

0(Ω) such that

• Ω

A∇uα · ∇w dx =

• Ω

(αf + 1 αϕ)w dx ∀w ∈ H1

0(Ω)

(23) and function α¯ u + 1

α ¯

v can be considered as an approximation of uα

• Then we can again use two-sided estimates

M ⊖

α,f,ϕ ≤ |

| |∇(αe + 1 αeϕ)| | |2 ≤ M ⊕

α,f,ϕ

(24)

21

✬ ✫ ✩ ✪

• We immediately observe that

1 2(M ⊖

α,f,ϕ − α2M ⊕ f − 1

α2 M ⊕

ϕ ) ≤ E1(e, eϕ)

and E1(e, eϕ) ≤ 1 2(M ⊕

α,f,ϕ − α2M ⊖ f − 1

α2 M ⊖

ϕ )

which together with directly computable term E0(¯ u, ¯ v) =

• Ω

f¯ v dx −

• Ω

A∇¯ v · ∇¯ u dx provide with two-sided estimates for the error ℓ(e) = E0 + E1

22

✬ ✫ ✩ ✪

Test No. 2

−1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 1 2 −0.02 −0.01 0.01 0.02 0.03 0.04 0.05

• A ≡ I, f ≡ 10, Th with 249 nodes
• ℓ(w) :=
• Ω

ϕw dx, supp ϕ = ω ⊂ Ω, ϕ ∈ L2(Ω), Tτ with 541 nodes

• Upper bound decreases from 0.04364 to 0.02433
• Lower bound grows from -0.01606 to 0.00591
• ℓ(e) =
• ω

ϕ(u − uh) dx ≈ 0.01317

23

✬ ✫ ✩ ✪

Final Comments

• Techniques proposed can be adapted for the other boundary

conditions and other linear elliptic problems (e.g., in linear elasticity), also for parabolic problems.

• For both types of control the corresponding upper and lower

bounds can be made arbitrary close to the true errors

• Such closeness only depends on computer resources (memory,

velocity) available

• Technologies can be easily coded and added as block-checker to

most of existing educational and industrial software products like MATLAB, FEMLAB, ANSYS, etc

24

✬ ✫ ✩ ✪

Compatibility With Existing Codes

CHECKER

Error Distribution Error + Error Distribution Total

THIS CONSIDERED AS A "BLACK BOX"

Data file Solver Grid Generator Local Error +

25

✬ ✫ ✩ ✪

Main Target of New Technologies

To create a set of numerical algorithms and respective computer codes, which is able to verify the accuracy of computed solutions. Such programs, CHECKERS, will be able to give the following information about the accuracy of obtained approximate solutions:

• two-sided estimation of errors in terms of special “goal-oriented”

quantities that can be proposed by engineers

• guaranteed two-sided bounds of the error in the global norm
• distribution of relative errors in subdomains

26

✬ ✫ ✩ ✪

Computational Benefits

By means of the above checkers, we can

• explicitely verify the error for the concrete computed solution

and, therefore, see when the desired tolerance is achieved

• verify the exisiting codes and discover their possible drawbacks

and bugs

• create optimal (or “quasi-optimal”) meshes in the process of

the mesh-refinement

• improve the quality of computer–aided modelling
• reduce the time of computer simulation

27

✬ ✫ ✩ ✪

Background

• The corresponding software is based on own very recent

mathematical results and present “know-how” technologies

• One of the main advantages of our technologies is that the

error estimators proposed are valid for a wide spectrum of different approximations, i.e., the checking code is absolutely independent of the solver used

• The connection between a particular solver and the checker can

be organized via the exchange of the data files

28

✬ ✫ ✩ ✪ INPUT DATA: a) geometry of the body b) material constants c) approximate solution obtained on the concrete mesh and the mesh parameters

• The input data are to be delivered as files of real numbers, the

structure of such files is to be specified OUTPUT DATA: a) two-sided estimates of errors in terms of special “goal-oriented” quantities suggested by customer b) two-sided bounds of the error in the global norm c) distribution of local errors

• The output data can be presented both as tables of numbers and

in graphical form

29