isoGeometric Residual Minimization Method (iGRM)

(1)

isoGeometric Residual Minimization Method (iGRM)

Maciej Paszyński¹

M. Łoś¹ I. Muga² Q. Deng³ V. M. Calo³

1Department of Computer Science

AGH University of Science and Technology, Kraków, Poland home.agh.edu.pl/paszynsk

2The University of Basque Country, Bilbao, Spain

3Curtin University, Perth, Western Australia

National Science Centre, Poland grant HARMONIA 2017/26/M/ST1/00281

1 / 45

(2)

Outline

Unstable problems

Factorization of Kronecker product matrices Residual minimization method

Factorization of residual minimization problem matrix Iterative solver

isoGeometric Residual Minimization method (iGRM) Numerical results: Manufactured solution problem, Eriksson-Johnson model problem

Conclusions

2 / 45

(3)

Eriksson-Johnson problem

Given Ω = (0, 1)², β = (1, 0)^T, we seek the solution of the advection-diffusion problem

∂u

∂x −  ∂²u

∂x² +∂²u

∂y²

!

= 0 (1)

with Dirichlet boundary conditions

u = 0 for x ∈ (0, 1), y ∈ {0, 1} u = sin(Πy ) for x = 0

The problem is driven by the inflow Dirichlet boundary condition.

It develops a boundary layer of width at the outflow x = 1.

3 / 45

(4)

Eriksson-Johnson problem

Figure:Galerkin FEM solution to Erikkson-Johnson problem, Pe = 10⁶

Figure:Cross-section and zoom

4 / 45

(5)

Higher continuity of isogeometric analysis

Figure:Isogeometric analysis provided smooth approximations with lower number of degrees of freedom

5 / 45

(6)

Higher continuity of isogeometric analysis

Figure:Recursive definition of B-splines

6 / 45

(7)

Higher continuity of isogeometric analysis

Figure:Tensor product structure of the 3D mesh

Isogeometric basis functions:

1D B-splines basis along x axis B_1,p^x (x ), . . . , B_N^x

x,p(x ) 1D B-splines basis along y axis B_1,p^y (y ), . . . , B_N^y_y_,p(y ) 1D B-splines basis along z axis B_1,p^z (z), . . . , B^z_N_z_,p(z) In 2D we take tensor product basis

{B_{i ,p}^x (x )B_j,p^y (y )}_{i =1,...,N}_x_;j=1,...,N_y In 3D we take tensor product basis

{B_{i ,p}^x (x )B_j,p^y (y )B_k,p^z (z)}_{i =1,...,N}_x_;j=1,...,N_y_;1,...,N_z

7 / 45

(8)

Mass and stiffness matrices over 2D domain Ω = Ω

_x

× Ω

_y

M= (Bij, Bkl)_L2 = Z

Ω

BijBkldΩ =

Z

Ω

B^x_i(x )B_j^y(y )B_k^x(x )B_l^y(y ) dΩ = Z

Ω

(B_iB_k)(x ) (B_jB_l)(y ) dΩ

=

Z

Ωx

BiBkdx

Z

Ωy

BjBldy

!

=M^x⊗ M^y

S= (∇Bij, ∇Bkl)L² = Z

Ω

∇Bij· ∇BkldΩ =

Z

Ω

∂(B_i^x(x )B_j^y(y ))

∂x

∂(B_k^x(x )B_l^y(y ))

∂x +∂(B_i^x(x )B_j^y(y ))

∂y

∂(B^x_k(x )B_l^y(y ))

∂y dΩ

= Z

Ω

∂B_i^x(x )

∂x B^y_j(y )∂B_k^x(x )

∂x B_l^y(y ) + B_i^x(x )∂B_j^y(y )

∂y B_k^x(x )∂B_l^y(y ))

∂y dΩ

= Z

Ωx

∂Bi

∂x

∂Bk

∂x dx Z

Ωy

BjBldy + Z

Ωx

BiBkdx Z

Ωy

∂Bj

∂y

∂Bl

∂y dy

=S^x⊗ M^y+ M^x⊗ S^y

8 / 45

(9)

Derivation of Spatial Direction Splitting

Idea exploit Kronecker product structure of Mx = b

with M = A ⊗ B, where A is n × n, B is m × m Definition of Kronecker (tensor) product:

M = A ⊗ B =







A B11 A B12 · · · A B1m

A B₂₁ A B₂₂ · · · A B_2m ... ... . .. ... A B_m1 A B_m2 · · · A B_mm







9 / 45

(10)

Derivation of Spatial Direction Splitting

RHS and solution are partitioned into m blocks of size n each xi = (xi 1, . . . , xin)^T

bi = (bi 1, . . . , bin)^T

We can rewrite the system as a block matrix equation:











AB₁₁x₁+ AB₁₂x₂ + · · · + AB_1mx_m = b₁ AB₂₁x₁+ AB₂₂x₂ + · · · + AB_2mx_m = b₂

... ... ... ... AB_m1x₁+ AB_m2x₂+ · · · + AB_mmx_m= b_m

10 / 45

(11)

Derivation of Spatial Direction Splitting

Factor out A:











A B11x1+ B12x2 + · · · + B1mxm

= b1

A B₂₁x₁+ B₂₂x₂ + · · · + B_2mx_m = b₂ ... ... ... ... A Bm1x1+ Bm2x2+ · · · + Bmmxm

= bm

Wy multiply by A⁻¹ and define yⁱ = A⁻¹bⁱ

(we have one 1D problem here A yⁱ = bⁱ with multiple RHS)











B₁₁x₁+ B₁₂x₂ + · · · + B_1mx_m = y₁ B₂₁x₁+ B₂₂x₂ + · · · + B_2mx_m = y₂

... ... ... ... B_m1x₁+ B_m2x₂+ · · · + B_mmx_m = y_m

11 / 45

(12)

Derivation of Spatial Direction Splitting

Consider each component of xi and yi ⇒ family of linear systems











B₁₁x¹ⁱ + B₁₂x²ⁱ + · · · + B_1mx^mi = y_1i B₂₁x¹ⁱ + B₂₂x²ⁱ + · · · + B_2mx^mi = y_2i

... ... ... ... B_m1x¹ⁱ + B_m2x²ⁱ+ · · · + B_mmx^mi = y_mi for each i = 1, . . . , n

⇒ linear systems with matrix B (We have another 1D problem here with multiple RHS B xⁱ = yⁱ )

12 / 45

(13)

Residual minimization method

b(u, v ) = l (v ) ∀v ∈ V (2)

For our weak problem (2) we define the operator B : U → V⁰ such as < Bu, v >_V⁰×V= b (u, v ).

B : U → V⁰ (3)

such that

hBu, v i_V⁰_×V = b(u, v ) (4) so we can reformulate the problem as

Bu − l = 0 (5)

We wish to minimize the residual u_h= argmin_w_h∈U_h

1

2kBw_h− lk²_V0 (6)

13 / 45

(14)

Residual minimization method

We introduce the Riesz operator being the isometric isomorphism R_V: V 3 v → (v , .) ∈ V⁰ (7) We can project the problem back to V

uh= argminw_h∈U_h

1

2kR_V⁻¹(Bwh− l)k²_V (8) The minimum is attained at u_h when the Gâteaux derivative is equal to 0 in all directions:

hR_V⁻¹(Bu_h− l), R_V⁻¹(B w_h)i_V = 0 ∀ w_h∈ U_h (9) We define the residual r = R_V⁻¹(Bu_h− l) and we get

hr , R_V⁻¹(B w_h)i = 0 ∀ w_h∈ U_h (10) which is equivalent to

hBw_h, r i = 0 ∀w_h∈ U_h. (11) From the definition of the residual we have

(r , v )_V = hBu_h− l, v i ∀v ∈ V . (12)

14 / 45

(15)

Residual minimization method with semi-infinite problem

Find (r , u_h)_{V ×U}_h such as

(r , v )_V − hBu_h− l, v i = 0 ∀v ∈ V

hBw_h, r i = 0 ∀w_h∈ U_h (13) We discretize the test space V_m ∈ V to get the discrete problem:

Find (rm, uh)Vm×U_h such as

(rm, vm)Vm− hBu_h− l, v_mi = 0 ∀v_m ∈ V_m

hBw_h, rmi = 0 ∀w_h∈ U_h (14) where (∗, ∗)_V_m is an inner product in Vm, hBu_h, vmi = b (u_h, vm), hBw_h, r_mi = b (w_h, r_m).

Remark

We define the discrete test space V_m in such a way that it is as close as possible to the abstract V space, to ensure stability, in a sense that the discrete inf-sup condition is satisfied. In our method it is possible to gain stability enriching the test space V_m without changing the trial space Uh.

15 / 45

(16)

Discretization of the residual minimization method

We approximate the solution with tensor product of one dimensional B-splines basis functions of order p

u_h=^X

i ,j

u_{i ,j}B_{i ;p}^x (x )B_j;p^y (y ). (15) We test with tensor product of one dimensional B-splines basis functions, where we enrich the order in the direction of the x axis from p to r (r ≥ p, and we enrich the space only in the direction of the alternating splitting)

v_m← B_{i ;r}^x (x )B_j;p^y (y ). (16) We approximate the residual with tensor product of one dimensional B-splines basis functions of order p

r_m =^X

s,t

r_s,tB_s;r^x (x )B_t;p^y (y ), (17) and we test again with tensor product of 1D B-spline basis functions of order r and p, in the corresponding directions

wh← B_k;p^x (x )B_{l ;p}^y (y ). (18) Remark

We perform the enrichment of the test space also in the alternating directions manner. In this way, when we solve the problem with derivatives along the x direction, we enrich the test space by increasing the B-splines order in the x direction, but we keep the B-splines order along y constant (same as in the trial space). By doing that, we preserve the Kronecker product structure of the matrix, to ensure that we can apply the alternating direction solver.

16 / 45

(17)

Decomposition into Kronecker product structure

A = A_y⊗ A_x; B = B_x ⊗ B_y; B^T = B_y^T⊗ B^T_x; A_y = B_y (19)

A B

B^T 0

!

= Ay 0

0 A^T_y

! Ax Bx

B_x^T 0

!

= AxAy BxAy

B_x^TA^T_y 0

! . (20) Both matrices Ay 0

0 A^T_y

!

and Ax Bx

B_x^T 0

!

can be factorized in a linear O(N) computational cost.

Figure:Factorization of second block. _{17 / 45}

(18)

Towards iterative solver

"

A B

B^T 0

# "

r u

#

=

"

F 0

#

(21) A = M + ηK

M = M_x ⊗ M_y,

K = Kx⊗ M_y+ Mx⊗ K_y. A= M + ηK

= (Mx + ηKx) ⊗ (My + ηKy) − η²Kx ⊗ K_y =A − ˜˜ K

"

A˜ B B^T 0

# "

r u

#

=

"

F +K˜r 0 ∗ u

#

(22) so we iterate

"

r^k+1 u^k+1

#

=

"

A˜ B B^T 0

#−1"

F +K˜r^k 0 ∗ u^k

#

(23)

18 / 45

(19)

Towards iterative solver

Firstly, consider the residual error

r^k = F − Aw^k − Bu^k = F + ˜K w^k− ˜Aw^k− Bu^k,

s^k = −B^Tw^k. (24)

To minimize these residual-type error, we build a dual problem, which resembles the features of preconditioner.

"

A˜ B B^T 0

# "

d^k c^k

#

=

"

r^k s^k

#

, (25)

Using (24), solving these equations gives

d^k = Ã⁻¹(F + ˜K w^k − Bu^k) − w^k− Ã⁻¹Bc^k := ˜d^k − w^k− Ã⁻¹Bc^k,

B^TA˜⁻¹Bc^k = −B^Td˜^k + B^Tw^k+ s^k = −B^T˜d^k,

(26)

where the first equation is a primal update and the second equation is to be solved by a Conjugate-Gradient (CG) type method.

19 / 45

(20)

Inner loop

Initialize {q⁽⁰⁾ = p⁽⁰⁾ = ˜c^k; u⁽⁰⁾ = u^k} for j = 1 until convergence

θ^(j)= Bp^(j); δ^(j)= ˜A⁻¹θ^(j); α^(j)= p^(j), q^(j)

θ^(j), δ^(j); u^(j+1)= u^(j)+ α^(j)p^(j); q^(j+1)= q^(j)− α^(j)B^Tδ^(j); β^(j+1)= q^(j+1), q^(j+1)

q^(j), q^(j) ; p^(j+1)= q^(j+1)+ β^(j+1)p^(j);

Algorithm 1: Inner CG

20 / 45

(21)

Iterative solver

To determine the convergence, for the inner loop,

we iterate until α^(j+1) ≤ tolerance. We denote this j as j_c. The outer iteration calculates

c^k = u^(j^c⁾− u⁽⁰⁾; u^k+1= u^(j^c⁾;

w^k+1= ˜d^k+ ˜A⁻¹Bc^k.

(27)

The outer iteration stops at c^k+1 ≤ tolerance.

21 / 45

(22)

Towards iterative solver

Initialize {u⁽⁰⁾ = 0; w⁽⁰⁾= 0}

for k = 1, ..., N until convergence Inner loop

c^k = u^(N^k⁾− u⁽⁰⁾; u^k+1= u^(N^k⁾;

w^k+1= ˜d^k+ ˜A⁻¹Bc^k; k = k + 1;

d˜^k+1=A˜⁻¹(F + ˜K w^k+1− Bu^k+1);

˜

c^k+1= −B^Td˜^k+1.

Algorithm 2: Inner-Outer CG

22 / 45

(23)

isogeometric Residual Minimization Method

isoGeometric Residual Minimization Method (iGRM) Residual minimization

(unconditional stability in space)

Discretization with B-spline basis functions (higher continuity smooth solutions)

Kronecker product structure of the inner product matrix (linear cost O(N) preconditioner for iterative solver) Symmetric positive definite system

(convergence of the conjugated gradient method )

23 / 45

(24)

A manufactured solution problem: strong form

We focus on a model problem with a manufactured solution. For a unitary square domain Ω = (0, 1)², the advection vector

β = (1, 1)^T, and Pe = 100, = 1/Pe we seek the solution of the advection-diffusion equation

∂u

∂x + ∂u

∂y −  ∂²u

∂x² +∂²u

∂y²

!

= f (28)

with Dirichlet boundary conditions u = g on the whole of Γ = ∂Ω.

We utilize a manufactured solution u(x , y ) = (x + e^Pe∗x − 1

1 − e^Pe )(y + e^Pe∗y − 1 1 − e^Pe )

enforced by the right-hand side, and we use homogeneous Dirichlet boundary conditions on ∂Ω.

24 / 45

(25)

A manufactured solution problem: weak form

b(u, v ) = l (v ) ∀v ∈ V (29) b(u, v ) =

∂u

∂x, v

Ω

+

∂u

∂y, v

Ω

+ 

∂u

∂x,∂v

∂x

Ω

+ 

∂u

∂y,∂v

∂y

Ω

−

∂u

∂xnx, v

Γ

−

∂u

∂yny, v

Γ

− (u, ∇v · n)_Γ− (u, β · nv )_Γ−u, 3p²/hv

Γ

n = (nx, ny) is versor normal to Γ, and h is element diameter, l (v ) = (f , v )_Ω− (g , ∇v · n)_Γ− (g , β · nv )_Γ−g , 3p²/hv

Γ (30) redterms correspond to weak imposition of the Dirichlet b.c. on Γ with g = 0, f is the manufactured solution, blueterms are the integration by parts,grayterms the penalty terms. We seek the solution in space U = V = H¹(Ω). The inner product in V is

(u, v )_V = (u, v )_L

2+

∂u

∂x,∂v

∂x

L2

+

∂u

∂y,∂v

∂y

L2

(31)

25 / 45

(26)

A manufactured solution results

n trial(2,1) trial(3,2) trial(4,3) trial(5,4)

test(2,0) test(2,0) test(2,0) test(2,0)

#DOF 389 410 433 458

L2 192 151 78 28

H1 101 74 44 32

8 × 8

#DOF 1413 1450 1489 1530

L2 80 16 3.29 1.48

H1 59 29 18 10

16 × 16

#DOF 5381 5450 5521 5594

L2 32 1.33 0.27 0.056

H1 31 9.77 3.16 0.82

32 × 32

#DOF 20997 21130 21265 21402

L2 7.66 0.07 0.01 0.003

H1 9.86 1.67 0.26 0.068

64 × 64 26 / 45

(27)

Eriksson-Johnson problem strong form

Given Ω = (0, 1)², β = (1, 0)^T, we seek the solution of the advection-diffusion problem

∂u

∂x −  ∂²u

∂x² +∂²u

∂y²

!

= 0 (32)

with Dirichlet boundary conditions

u = 0 for x ∈ (0, 1), y ∈ {0, 1} u = sin(Πy ) for x = 0

The problem is driven by the inflow Dirichlet boundary condition.

It develops a boundary layer of width at the outflow x = 1.

27 / 45

(28)

Eriksson-Johnson problem weak form

We introduce first the weak formulation for the Eriksson-Johnson problem

b(u, v ) = l (v ) ∀v ∈ V (33) b(u, v ) =

∂u

∂x, v

+ 

∂u

∂x,∂v

∂x

+ 

∂u

∂y,∂v

∂y

−

∂u

∂xnx, v

Γ

−

∂u

∂yny, v

Γ

− (u, ∇v · n)_Γ− (u, β · nv )_Γ−u, 3p²/hv

Γ

where theblue,red, andgrayterms correspond to the weak

imposition of the Dirichlet b.c. on the boundary Γ, and n = (n_x, n_y) is the versor normal to the boundary,

l (v ) =− (g , ∇v · n)_Γ−− (g , β · nv )_Γ−−g , 3p²/hv

Γ⁻ (34) where thered andgrayterms correspond to the weak introduction of the Dirichlet b.c. on the boundary Γ.

We plug the weak form and the inner product into the iGRM setup and we use the preconditioned CG solver. ^{28 / 45}

(29)

Residual minimization method for Eriksson-Johnson problem

Find (r_m, u_h)_V_m_×U_h such as (rm, vm)Vm−

∂uh

∂x , vm

− 

∂uh

∂x ,∂vm

∂x +∂uh

∂y ,∂vm

∂y

= (f , vm)

∀v_m∈ V_m

∂wh

∂x , rm

+ 

∂wh

∂x ,∂rm

∂x +∂wh

∂y ,∂rm

∂y

= 0

∀w_h∈ U_h where (r_m, v_m)_V_m= (r_m, v_m) + (^∂r_∂x^m,^∂v_∂x^m) + (^∂r_∂y^m,^∂v_∂y^m)

is the H¹ norm induced inner product.

Remark

We will use trial space U_h as quadratic B-splines with C¹ continuity and test space V_h as cubic B-splines with C² continuity

29 / 45

(30)

Numerical results for the Eriksson-Johnson problem

Figure:Left panel: Solution to the Erikkson-Johnsson problem with Galerkin method with = 10⁻², with trial = test = quadratic B-splines on a uniform mesh of 20 × 20 elements. Right panel: Cross-section at y = 0.5.

30 / 45

(31)

Numerical results for the Eriksson-Johnson problem

Figure:Solution to the Erikkson-Johnsson problem withresidual minimization methodwith = 10⁻², with trial = quadratic B-splines, test = cubic B-splines on a uniform mesh of 20 × 20 elements.

31 / 45

(32)

Numerical results for the Eriksson-Johnson problem

Figure:Cross-section at y = 0.5 of the solution to the Erikkson-Johnsson problem withresidual minimization methodwith = 10⁻², with trial = quadratic B-splines, test = cubic B-splines

on a uniform mesh of 20 × 20 elements.

32 / 45

(33)

Numerical results for the Eriksson-Johnson problem

Figure:Zoom of the cross-section at y = 0.5 of the solution at (0.99,1.0) of the Erikkson-Johnsson problem withresidual minimization methodwith

 = 10⁻², with trial = quadratic B-splines, test = cubic B-splines on a uniform mesh of 20 × 20 elements.

33 / 45

(34)

Eriksson-Johnson problem weak form, SUPG method

b(u, v ) +(R(u), τ β · ∇v )= l (v ) ∀v ∈ V (35) where R(u) = ^∂u_∂x + ∆u, and^β_h^x

x +^β_h^y

y

+ 3_h2¹ x+h²_y,

and in our case diffusion term = 10⁻⁶, and convection term β = (1, 0), and h_x and h_y are dimensions of an element.

b_SUPG(u, v ) = l (v ) ∀v ∈ V (36)

b_SUPG(u, v ) =

∂u

∂x, v

+ 

∂u

∂x,∂v

∂x

+ 

∂u

∂y,∂v

∂y

−

∂u

∂xnx, v

Γ

−

∂u

∂yny, v

Γ

− (u, ∇v · n)_Γ− (u, β · nv )_Γ−u, 3p²/hv

Γ

+





∂u

∂x + ∆u, 1 hx

+ 3 1 h²_x + h_y²

!2

∂v

∂x





34 / 45

(35)

Numerical results for the Eriksson-Johnson problem

0 0.2 0.4 0.6 0.8 1 1.2

0 0.2 0.4 0.6 0.8 1

(a)iGRM (2,1), (3,0) (b)iGRM

0 0.2 0.4 0.6 0.8 1 1.2

0.999990 0.999995 1

(c)iGRM zoom

0 0.2 0.4 0.6 0.8 1 1.2

0 0.2 0.4 0.6 0.8 1

(d)SUPG trial (3,2) (e)SUPG

0 0.2 0.4 0.6 0.8 1 1.2

0.999990 0.999995 1

(f)SUPG zoom Figure:Comparison of solutions of the Erikkson-Johnsson by using iGRM and SUPG methods for Pe=1,000,000 on 2x2 mesh,

using (2,1) for trial and (3,0) for testing.

35 / 45

(36)

Numerical results for the Eriksson-Johnson problem

0 0.2 0.4 0.6 0.8 1 1.2

0 0.2 0.4 0.6 0.8 1

(a)iGRM (2,1), (3,0) (b)iGRM

0 0.2 0.4 0.6 0.8 1 1.2

0.999990 0.999995 1

(c)iGRM zoom

0 0.2 0.4 0.6 0.8 1 1.2

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1 1.2

0.999990 0.999995 1

(f)SUPG zoom Figure:Comparison of solutions of the Erikkson-Johnsson by using iGRM and SUPG methods for Pe=10000 on 4x4 mesh,

36 / 45

(37)

Numerical results for the Eriksson-Johnson problem

0 0.2 0.4 0.6 0.8 1 1.2

0 0.2 0.4 0.6 0.8 1

(a)iGRM (2,1), (3,0) (b)iGRM

0 0.2 0.4 0.6 0.8 1 1.2

0.999990 0.999995 1

(c)iGRM zoom

0 0.2 0.4 0.6 0.8 1 1.2

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1 1.2

0.999990 0.999995 1

37 / 45

(38)

Numerical results for the Eriksson-Johnson problem

0 0.2 0.4 0.6 0.8 1 1.2

0 0.2 0.4 0.6 0.8 1

(a)iGRM (2,1), (3,0) (b)iGRM

0 0.2 0.4 0.6 0.8 1 1.2

0.999990 0.999995 1

(c)iGRM zoom

0 0.2 0.4 0.6 0.8 1 1.2

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1 1.2

0.999990 0.999995 1

38 / 45

(39)

Numerical results for the Eriksson-Johnson problem

0 0.2 0.4 0.6 0.8 1 1.2

0 0.2 0.4 0.6 0.8 1

(a)iGRM (2,1), (3,0) (b)iGRM

0 0.2 0.4 0.6 0.8 1 1.2

0.999990 0.999995 1

(c)iGRM zoom

0 0.2 0.4 0.6 0.8 1 1.2

0 0.2 0.4 0.6 0.8 1

(d)SUPG trial (3,2) (e)SUPG crossection

0 0.2 0.4 0.6 0.8 1 1.2

0.999990 0.999995 1

(f)SUPG crossection zoom

Figure:Comparison of solutions of the Erikkson-Johnsson by using iGRM and SUPG methods for Pe=1,000,000 on 32x32 mesh,

39 / 45

(40)

Erikkson-Johnson iterations

2 4 8 16 32

n 0

5 10 15 20

#ofouteriterations

2 4 8 16 32

n 0

50 100 150 200

#ofCGiterations

Figure:Left panel: Number of iterations of outer loop for. Right panel:

Number of iterations of inner loop (CG).

2 4 8 16 32

n 0.0

0.5 1.0 1.5 2.0 2.5 3.0

#ofouteriterations

Figure:Number of iterations of outer loop, with inner loop treated by

MUMPS solver. _{40 / 45}

(41)

Circular wind problem

β_x∂u

∂x + β_y∂u

∂y −  ∂²u

∂x² +∂²u

∂y²

!

= 0

over the rectangular domain Ω = (0, 1) × (−1, 1), with zero right-hand side f = 0, Pe = 1, 000, 000, the advection vector β(x , y ) = (β_x(x , y ), β_y(x , y )) = ψ( ^−y

(x²+y²)¹², ^x

(x²+y²)¹²) modeling the circular wind, where ψ is the wind force coefficient.

Γ₁ = {(x , y ) : x = 0, 0.5 ≤ y ≤ 1.0}, Γ2 = {(x , y ) : x = 0, 0.0 ≤ y ≤ 0.5}, Γ₃ = {(x , y ) : x = 0, −0.5 ≤ y ≤ 0.0}, Γ₄ = {(x , y ) : x = 0, −1.0 ≤ y ≤ −0.5},

We utilize the Dirichlet boundary conditions u = g on Γ = ∂Ω where g = 1

2

tanh

(|y | − 0.35)b

+ 1

, for x ∈ Γ2∪ Γ3

g = 1 2

0.65 − tanh

(|y |)b

+ 1

, for x ∈ Γ1∪ Γ4

g = 0, for x ∈ Γ \ Γ₁∪ Γ2∪ Γ3∪ Γ4 _{41 / 45}

(42)

Circular wind problem

The weak formulation

b(u, v ) = l (v ) ∀v ∈ V

b(u, v ) =

β_x∂u

∂x, v

Ω

+

β_y∂u

∂y, v

Ω

+ 

∂u

∂x,∂v

∂x

Ω

+ 

∂u

∂y,∂v

∂y

Ω

−

∂u

∂xn_x, v

Γ

−

∂u

∂yn_y, v

Γ

− (u, ∇v · n)_Γ− (u, β · nv )_Γ−u, 3p²/hv

Γ

where n = (n_x, n_y) is the versor normal to Γ,

l (v ) =− (g , ∇v · n)_Γ− (g , β · nv )_Γ−g , 3p²/hv

Γ

n = (n_x, n_y) is the versor normal to the boundary, and the right-hand side forcing is equal to 0.

42 / 45

(43)

Circular wind problem

Figure:Solution to the circular wind problem on the mesh of 128 × 128 elements with trial(2,1),test(2,0), for Pecklet number Pe = 1, 000, 000, wind force b = 1. Horizontal cross-section at x = 0.

43 / 45

(44)

Conclusions and open questions

Conclusions

We use B-spline basis functions on tensor product patches of elements and we do not break the test spaces

We employ the residual minimization method with fixed trial space, and we enrich the test space to improve the stability To preserve the Kronecker product structure of the residual minimization system, we enrich the test space in the alternating direction manner (e.g. (3,2) for x direction, and (2, 3) for y direction) to preserve the Kronecker product structure of the matrix, and linear cost O(N) factorization For a linear cost O(N) Kronecker product precondtioner we can enrich the test space in arbitrary way, since we only need a Kronecker product structure of the Gramm matrix

Future work

Stokes, Oseen, Navier-Stokes and Maxwell equations This work is supported by National Science Centre, Poland grant

HARMONIA 2017/26/M/ST1/00281. _{44 / 45}

(45)

Papers

M. Los, Q.Deng, I. Muga, V.M.Calo, M. Paszynski, Isogeometric Residual Minimization Method (iGRM) with Direction Splitting Preconditoner for Stationary Advection-Diffusion Problems, submitted to Computer Methods in Applied Mechanics and Engineering (2019) IF: 4.441

Code based on:

M. Los, M. Wozniak, M. Paszynski, A. Lenharth, K. Pingali, IGA-ADS : Isogeometric Analysis FEM using ADS solver, Computer

& Physics Communications, 217 (2017), pp. 99-116 IF: 3.748

45 / 45