Very high-order finite volume scheme for the 2D linear - - PowerPoint PPT Presentation

Very high-order finite volume scheme for the 2D linear convection-diffusion problem

• S. Clain1,3, G.J. Machado1, J.M. Nobrega2, R.M.S. Pereira1,
• A. Boularas4

1Mathematical Centre, University of Minho, Portugal 2Institute for Polymers and Composites/I3N, University of Minho, Portugal 3Mathematical Institute, Paul Sabatier University, France 4Laplace centre, Paul Sabatier University, France

Ofir, 2014 april-28th may-2nd

Outline

➓ 2D linear convection-diffusion with finite volume ➓ Discretization and residual formulation ➓ Polynomial reconstructions ➓ Scheme design ➓ Numerical tests

Model

Find φ ✏ ♣φ1, φ2q on the bounded open domain Ω such that ∇.♣V1φ1 ✁ k1∇φ1q ✏ f1, in Ω1, (1a) ∇.♣V2φ2 ✁ k2∇φ2q ✏ f2, in Ω2, (1b) k1∇φ1.nΓ ✏ k2∇φ2.nΓ,

• n Γ,

(1c) φ ✏ φD,

• n ΓD,

(1d) φ1 ✏ φ2 on Γ

• r

k1∇φ1.nΓ ✏ h♣φ2 ✁ φ1q .

Discretization

➺

❇ci

♣V.nφ ✁ k∇φ.nqds ✁ ➺

ci

f dx ✏ 0. (divergence theorem) ➳

jPν♣iq

⑤eij⑤

R

➳

r✏1

ζr ✑ V♣qij,rq.nijφ♣qij,rq ✁ k♣qij,rq∇φ♣qij,rq.nij ✙ ✁⑤ci⑤fi ✏ O♣h2R

i q. (Gauss Points)

Residual formulation

Based on the previous expression: the residual formulation is Gi ✏ ➳

jPν♣iq

⑤eij⑤ ⑤ci⑤

R

➳

r✏1

ζrFij,r ✁ fi, where Fij,r ✓ V♣qij,rq.nijφ♣qij,rq ✁ k♣qij,rq∇φ♣qij,rq.nij, approximation of the flux at the Gauss point qij,r. ☞ Sixth-order approximation for Fij,r.

Polynomial reconstructions

Conservative reconstruction for cells

ci: cell of mesh Th, d: the polynomial degree, S♣ci, dq: the associated stencil, φi: mean value on ci. ♣ φi♣x; dq ✏ φi ➳

1↕⑤α⑤↕d

Rd,α

i

✦ ♣x ✁ biqα ✁ Mα

i

✮ α ✏ ♣α1, α2q, ⑤α⑤ ✏ α1 α2, x ✏ ♣x1, x2q, bi the centroid of cell ci ☞ Set Mα

i ✏ 1

⑤ci⑤ ➺

ci

♣x ✁ biqα dx to provide the conservation.

Coefficients for ♣ φi

Rd

i vector gathering coefficients Rd,α i

Assume mean values φℓ on cells cℓ, ℓ P S♣ci, dq are known, ♣ Rd

i minimizes the functional

Ei♣Rd

i ; dq ✏

➳

ℓPS♣ci,dq

✑ 1 ⑤cℓ⑤ ➺

cℓ

♣ φi♣x; dq dx ✁ φℓ ✙2 , ☞ Lead to an over-determined linear system Ad

i Rd i ✏ bS i where

bS

i represents the variations φℓ ✁ φi, ℓ P S♣ci; dq

Preconditioning and solving

Determine the Moore-Penrose pseudo-inverse matrix for system ♣Ad

i Pd i q♣Pd i q✁1Rd i ✏ bS i with the diagonal matrix

Pd

i ✏ diag♣⑤ci⑤✁⑤α⑤④2q1↕⑤α⑤↕d.

Motivation: the Ad

i matrix coefficients are

1 ⑤cℓ⑤ ➺

cℓ

♣x ✁ biqαdx ✁ 1 ⑤ci⑤ ➺

ci

♣x ✁ biqαdx. ☞ Strongly reduces the effect of the power α. We compute the pseudo inverse matrix ♣Ad

i Pd i q✿ and get

Rd

i ✏ Pd i ♣Ad i Pd i q✿bS i .

Conservative reconstruction for Γ

eij ⑨ Γ: edge on the interface, cj the cell on the Ω2 side, d: the polynomial degree, S♣cj, dq: the associated stencil. q φj♣x; dq ✏ φj ➳

1↕⑤α⑤↕d

Rd,α

j

✦ ♣x ✁ bjqα ✁ Mα

j

✮ . bj the centroid of cell cj. ☞ Set Mα

j ✏ 1

⑤cj⑤ ➺

cj

♣x ✁ bjqα dx to provide the conservation. ☞ Only use the cells on the ”j” side.

Coefficients for q φj

Rd

j vector gathering coefficients Rd,α j

Assume mean values φℓ on cells cℓ, ℓ P S♣eij, dq and φij the mean value on eij are known, q Rd

j minimizes the functional

Ej♣Rd

j ; dq ✏

➳

ℓPS♣cj,dq

✑ 1 ⑤cℓ⑤ ➺

cℓ

q φj♣x; dq dx ✁ φℓ ✙2 ωij ✑ 1 ⑤eij⑤ ➺

eij

q φj♣x; dq ds ✁ φij ✙2 with ωij a positive weight.

Conservative reconstruction for ΓD

eiD: edge on the boundary ΓD, d: the polynomial degree, S♣eiD, dq: the associated stencil, φiD ✏ 1 ⑤eiD⑤ ➺

eiD

φD♣sq ds. ♣ φiD♣x; dq ✏ φiD ➳

1↕⑤α⑤↕d

Rd,α

iD

✦ ♣x ✁ miDqα ✁ Mα

iD

✮ , miD the centroid of edge eiD ☞ Set Mα

iD ✏

1 ⑤eiD⑤ ➺

eiD

♣x ✁ miDqα ds to provide the conservation.

Coefficients for ♣ φiD

Rd

iD vector gathering coefficients Rd,α iD

Assume mean values φℓ on cells cℓ, ℓ P S♣eiD, dq are known, ♣ Rd

iD minimizes the functional

EiD♣Rd

iD; dq ✏

➳

ℓPS♣eiD,dq

ωiD,ℓ ✑ 1 ⑤cℓ⑤ ➺

cℓ

♣ φiD♣x; dq dx ✁ φℓ ✙2 , with ωiD,ℓ positive weights. ☞ Coefficients ωiD,ℓ are very important to provide ”good” properties.

Non conservative reconstruction for inner edges

eij: inner edge of the mesh, d: the polynomial degree, S♣eij, dq: the associated stencil, no value associated to eij. r φij♣x; dq ✏ ➳

0↕⑤α⑤↕d

Rd,α

ij ♣x ✁ mijqα

mij the centroid of edge eij ☞ Coefficient Rd,α

ij

for ⑤α⑤ ✏ 0 is also unknown.

Coefficients for r φij

Rd

ij vector gathering coefficients Rd,α ij

Assume mean values φℓ on cells cℓ, ℓ P S♣eij, dq are known, r Rd

ij minimizes the functional

Eij♣Rd

ij; dq ✏

➳

ℓPS♣eij,dq

ωij,ℓ ✑ 1 ⑤cℓ⑤ ➺

cℓ

r φij♣x; dq dx ✁ φℓ ✙2 , where ωij,ℓ are positive weights. ☞ One more time: coefficients ωij,ℓ are very important.

Polynomial Reconstruction Operators

✌ Φ ✏ ♣φiqiPC the mean values vector. ✌ Operators Φ Ñ ♣ φi, q φj, ♣ φiD, and r φij are linear. ✌ φ polynomial function of degree d, φi the exact mean values. d-exact reconstruction if ♣ φi♣x; dq ✏ q φj♣x; dq ✏ ♣ φiD♣x; dq ✏ r φij♣x; dq ✏ φ♣xq, x P R2. ☞ The finite volume method associated to the polynomial reconstruction is a d 1th-order method.

The flux on edge (except Γ)

• 1. eij is an inner edge (not on Γ),

Fij,r ✏ rV♣qij,rq.nijs♣ φi♣qij,r; dq rV♣qij,rq.nijs✁♣ φj♣qij,r; dq ✁k♣qij,rq∇r φij♣qij,r; dq.nij.

• 2. eiD belongs to ΓD,

FiD,r ✏ rV♣qiD,rq.niDs♣ φi♣qiD,r; dq rV♣qiD,rq.niDs✁φD♣qiD,rq ✁k♣qiD,rq∇♣ φiD♣qiD,r; dq.niD.

The flux on edge of eij ⑨ Γ

✌ Transfer condition: Fij,r ✏ h♣qij,rqr♣ φi♣qij,r; dq ✁ ♣ φj♣qij,r; dqs. ✌ Continuity condition: we perform three steps Step 1: compute the reconstructions ♣ φi♣x; dq for all ci ⑨ Ω1. Step 2: compute φij ✏ 1 ⑤eij⑤ ➺

eij

♣ φi♣x; dqds for all eij ⑨ Γ. Step 3: Fij,r ✏ k2♣qij,rq∇q φj♣qij,r; dq.nij.

Resolution

① The polynomial reconstruction operators are linear. ② The flux computations are linear. ③ The residual expression is linear: Φ Ñ Gi♣Φq. We get a linear operator Φ Ñ G♣Φq ✏ ♣G1♣Φq, ..., GI♣Φqq. Problem: Find Φ such that G♣Φq ✏ 0. ☛ Matrix-free problem: use GMRES method. Preconditioning is very very important: P preconditioning matrix substitute Φ Ñ G♣Φq by Φ Ñ PG♣Φq.

Preconditioning matrix

☞ G♣Φq ✏ AΦ ✁ b but we do not have matrix A: ILU not possible. Diagonal preconditioning matrix P ✏ D✁1

P

with DP♣i, iq ✏ 1 ⑤ci⑤ ➳

jPν♣iq

⑤eij⑤ ✒k♣biq ⑤bibj⑤ rV♣mijq.nijs ✚ . More sophisticated preconditioning matrix, AP ✏ DP for the diagonal coefficients and AP♣i, jq ✏ ⑤eij⑤ ⑤ci⑤ ✒ ✁k♣mijq ⑤bibj⑤ rV♣mijq.nijs✁ ✚ , j P ν♣iq.

Incomplete inverse of AP

Preconditioning matrix is supposed to be P ✏ A✁1

P

Substitute with the incomplete inverse A✿

P with the same

non-null entries of AP. Taking advantage of the structure of AP and A✿

P provides explicit

construction of A✿

P

A✿

P♣i, jq ✏ ✁AP♣i, jqA✿ P♣i, iq

AP♣j, jq, j P ν♣iq with A✿

P♣i, iq ✏

1 AP♣i, iq ✁ ➳

jPν♣iq

AP♣i, jqAP♣j, iq AP♣j, jq .

Curved boundary treatment

☞ Problem: we set the boundary condition on the edge while it is prescribed on the curve. ✌ qiD,r Gauss points on eiD, ✌ piD,r Gauss points on the curve, ✌ ♣ φiD is evaluated using φD♣qiD,rq and not φD♣piD,rq. ☞ Idea: modify the mean value φiD ✏

R

➳

r✏1

ζrφD♣qiD,rq ds but φiD is still associated to edge eiD.

Curved boundary treatment

Find φiD approximation on edge eiD ⑨ ΓD such that φD♣piD,rq ✁ ♣ φiD♣piD,rq is minimal.

• 1. Initialize φ0

iD ✏ ➦R r✏1 ζrφD♣qiD,rq and evaluate ♣

φ0

iD

• 2. Do

➓ Compute δk

iD,r ✏ φD♣piD,rq ✁ ♣

φk

iD♣piD,rq (boundary default),

➓ Update φiD with φk1

iD

:✏ φk

iD R

➳

r✏1

ζrδk

iD,r,

➓ Evaluate ♣

φk1

iD

with new mean value φk1

iD ,

While(⑤φk1

iD

✁ φk

iD⑤ ➔ Tol)

• 3. Compute the flux on boundary with the update ♣

φiD.

Examples and numerical simulations

M-matrix

✌ The underlying matrix A must be an M-matrix (stability and positivity preserving). ✌ Number of non negative coefficients: w ✏ 2 gives 46%, w ✏ 2.5 gives 25%, w ✏ 3. gives 0%. ☞ I have to check if aij ↕ 0 with i ✘ j.

Convection Diffusion

Smooth solution with low and high Peclet number.

1.0e-10 1.0e-08 1.0e-06 1.0e-04 1.0e-02 5.0e-03 1.0e-02 2.0e-02 4.0e-02 8.0e-02

error

L1-norm convergence curves: low Peclet number

P1-reconstruction P3-reconstruction P5-reconstruction 1.0e-06 1.0e-05 1.0e-04 1.0e-03 5.0e-03 1.0e-02

error

L1-norm convergence curves: high Peclet number

P1-reconstruction P3-reconstruction P5-reconstruction

Pure convection

Revolution of a smooth pattern. Velocity V ✏ ♣✁y, xq Dirichlet condition for the inflow boundary.

1.0e-08 1.0e-06 1.0e-04 1.0e-02 5.0e-03 1.0e-02 2.0e-02 4.0e-02

error

L1-norm convergence curves: circular shape

P1-reconstruction P3-reconstruction P5-reconstruction

Pure diffusion with discontinuous coefficients

k1∇φ1 ✏ h♣φ2 ✁ φ1q φ1 ✏ φ2

1.0e-08 1.0e-06 1.0e-04 1.0e-02 2.0e-02 4.0e-02 8.0e-02

error

L1-norm convergence curves: transfer condition

P1-reconstruction P3-reconstruction P5-reconstruction 1.0e-09 1.0e-07 1.0e-05 1.0e-03 1.0e-01 2.0e-02 4.0e-02 8.0e-02

error

L1-norm convergence curves: continuity condition

P1-reconstruction P3-reconstruction P5-reconstruction

Convection diffusion with discontinous coefficients

k1∇φ1 ✏ h♣φ2 ✁ φ1q φ1 ✏ φ2

1.0e-08 1.0e-06 1.0e-04 1.0e-02 2.0e-02 4.0e-02 8.0e-02

error

L1-norm convergence curves: transfer condition

P1-reconstruction P3-reconstruction P5-reconstruction 1.0e-09 1.0e-07 1.0e-05 1.0e-03 1.0e-01 2.0e-02 4.0e-02 8.0e-02

error

L1-norm convergence curves: continuity condition

P1-reconstruction P3-reconstruction P5-reconstruction

The ring problem

∆φ ✏ 0, φ ✏ 25 at R ✏ 1000e ✁ 09, φ ✏ 0 at 100.e ✁ 09, Exact solution a b ln♣rq.

1.0e-02 5.0e-03 1.0e-02 2.0e-02

error

L1-norm convergence curves: rough Dirichlet condition

P1-reconstruction P3-reconstruction P5-reconstruction 1.0e-09 1.0e-07 1.0e-05 1.0e-03 1.0e-01 5.0e-03 1.0e-02 2.0e-02

error

L1-norm convergence curves: acccurate Dirichlet condition

P1-reconstruction P3-reconstruction P5-reconstruction

The cooler

• 5.0e-03

2.2e-02 5.0e-02 0e+00 1e-02 x y

Preconditioning

Residual histogram for the three preconditioning matrices versus the order method. Simulations of the cooler with a 23616 triangles mesh .

Conclusions

✌ Very high-order is achieved even for discontinuous coefficients or discontinuous function at the interface. ✌ Control of the M-matrix via the weights ✌ no spurious oscillation, high Peclet number. ✌ Curved boundaries seems resolved. ✌ 3D is in progress (curved boundary?) ✌ Non-stationary problem in progress. ✌ Next steps: non linear case, Stokes, Navier-Stokes (MOOD)