Fast and smooth simulation of space-time problems

(1)

Fast and smooth simulation of space-time problems

Day 1

Department of Computer Science

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

1 / 43

(2)

Department of Computer Science AGH University, Kraków, Poland

2 / 43

(3)

Outline

Isogeometric finite element method

Alternating Directions Implicit (ADI) method Isogeometric L2 projections

Explicit dynamics

Example 1: Heat transfer Installation of IGA-ADS solver

Parallel distributed memory explicit dynamics Parallel shared memory explicit dynamics

Example 2: Non-linear flow in heterogenous media Implicit dynamics

Example 3: Implicit heat transfer Example 4: Linear elasticity Example 5: Pollution problem Labs with implict dynamics

3 / 43

(4)

Software

Program Title: IGA-ADS

Code: git clone https://github.com/marcinlos/iga-ads Licensing provisions: MIT license (MIT)

Programming language: C++

Nature of problem: Solving non-stationary problems in 1D, 2D and 3D

Solution method: Alternating direction solver with isogeometric finite element method

If you use this software in your work, please cite

Marcin Łoś, Maciej Woźniak, Maciej Paszyński, Andrew Lenharth, Keshav Pingali IGA-ADS : Isogeometric Analysis FEM using ADS solver, Computer & Physics Communications 217 (2017) 99-116 (available on researchgate.org)

4 / 43

(5)

Isogeometric finite element method

J.A. Cottrel, T.J.R. Hughes, Y. Bazilevs, Isogeometric Analysis.

Toward Integration of CAD and FEA, Wiley, (2009).

5 / 43

(6)

Isogeometric finite element method

Original recursive definition of B-spline basis functions

Figure:Recursive formulae for B-spline basis functions

6 / 43

(7)

Isogeometric finite element method

How to remember this formulae graphically

Figure:Practical implementation of the recursive formulae for B-spline basis functions

7 / 43

(8)

Isogeometric finite element method

How these B-spline basis functions look like

Figure:Basis functions of order 0,1,2 for uniform knot vector {0,1,2,3,4,5}

8 / 43

(9)

Isogeometric finite element method

Representation of B-splines by knot vectors

Figure:B-spline basis functions represented by knot vector {0,0,0,1,2,3,4,4,5,5,5}

9 / 43

(10)

Isogeometric finite element method

Quadratic B-spline basis functions represented by knot vector {0,0,0,1,2,3,4,4,5,5,5}

B-spline curve:

N1,2*(0,1)+N2,2*(1,0)+N3,2*(2,0)+N4,2*(2,2)+N5,2*(4,3)+

N_6,2*(4,4)+N_7,2*(2,4)+N_8,2*(1,2)

10 / 43

(11)

Alternating Direction Implicit (ADI) method

The Alternating Direction Implicit (ADI) method

G. Birkhoff, R.S. Varga, D. Young, Alternating direction implicit methods, Advanced Computing (1962)

du

dt − L_xu − L_yu = f du

dt −ui −1,j − 2u_{i ,j}+ ui +1,j

h² −ui ,j−1− 2u_{i ,j} + ui ,j+1

h² = f

u^t+0.5_{i ,j} − u_{i ,j}^t

dt −u^t+0.5_{i −1,j} − 2u_{i ,j}^t+0.5+ u^t+0.5_{i +1,j}

h² = u^t_{i ,j−1}− 2u_{i ,j}^t + u_{i ,j+1}^t h² +f_{i ,j}^t u_{i ,j}^t+1− u^t+0.5_{i ,j}

dt −u^t+1_{i ,j−1}− 2u_{i ,j}^t+1+ u_{i ,j+1}^t+1

h² =

u^t+0.5_{i −1,j} − 2u_{i ,j}^t+0.5+ u_{i +1,j}^t+0.5

h² + f_{i ,j}^t+0.5

11 / 43

(12)

Alternating Direction Implicit (ADI) method

u_{i −1,j}^t+0.5[−2dt

h² ] + u_{i ,j}^t+0.5[1 + 2dt

h² ] + u_{i +1,j}^t+0.5[−2dt h² ] = dtu^t_{i ,j−1}− 2u_{i ,j}^t + u^t_{i ,j+1}

h² + dtf_{i ,j}^t for i = 1, ..., N_x, j = 1, ..., N_y.

u_{i ,j−1}^t [−2dt

h² ] + u^t_{i ,j}[1 + 2dt

h² ] + u_{i ,j+1}^t [−2dt h² ] = dtu_{i −1,j}^t+0.5− 2u_{i ,j}^t+0.5+ u_{i −1,j}^t+0.5

h² + dtf_{i ,j}^t+0.5 for j = 1, ..., N_y, i = 1, ..., N_x.

12 / 43

(13)

Alternating Direction Implicit (ADI) method

u_{i −1,j}^t+0.5[−2dt] + u_{i ,j}^t+0.5[h²+ 2dt] + u_{i +1,j}^t+0.5[−2dt] = dtu_{i ,j−1}^t − 2u^t_{i ,j}+ u_{i ,j+1}^t + h²dtf_{i ,j}^t for i = 1, ..., N_x, j = 1, ..., N_y.

u_{i ,j−1}^t [−2dt] + u^t_{i ,j}[h²+ 2dt] + u_{i ,j+1}^t [−2dt] = u^t+0.5_{i −1,j} − 2u_{i ,j}^t+0.5+ u^t+0.5_{i −1,j} + h²dtf_{i ,j}^t+0.5 for j = 1, ..., Ny, i = 1, ..., Nx.

13 / 43

(14)

Alternating Direction Implicit (ADI) method







h²+ 2dt −2dt 0 · · · · · · · · · · · · 0

−2dt h²+ 2dt −2dt 0 · · · · · · · · · 0

0 −2dt h²+ 2dt −2dt 0 · · · · · · 0

.. .

0 · · · · · · 0 −2dt h²+ 2dt −2dt

0 · · · · · · · · · 0 −2dt h²+ 2dt













u^t+0.5_1,1 u^t+0.5 1,2 u^t+0.5_1,3

.. . u^t+0.5

Nx ,Ny −1 u^t+0.5

Nx ,Ny







=







−2u^t_1,1+ u_1,2^t + h²dtf_1,1^t u^t_1,1− 2u^t_1,2+ u^t_1,3+ h²dtf_{i ,j}^t

.. . u^t

Nx ,Ny −2− 2u^t

Nx ,Ny −1+ u^t

Nx ,Ny+ h²dtf^t Nx ,Ny −1 u^t

Nx ,Ny −1− 2u^t

Nx ,Ny+ h²dtf^t Nx ,Ny







14 / 43

(15)

Isogeometric L2 projections

Longfei Gao, Kronecker Products on Preconditioning, PhD. Thesis, KAUST (supervised by Victor Calo), 2013.

Isogeometric basis functions:

1D B-splines basis B₁(x ), . . . , B_n(x ) higher dimensions: tensor product basis Bi1···i_d(x1, . . . , x_d) ≡ B_i^x₁¹(x1) · · · B_i^x^d

d (x_d)

15 / 43

(16)

Isogeometric L2 projections

Gram matrix of B-spline basis on 2D domain Ω = Ωx× Ω_y: M_ijkl = (B_ij, B_kl)_L2 =

Z

Ω

B_ijB_kldΩ

= Z

Ω

B_i^x(x )B_j^y(y )B_k^x(x )B_l^y(y ) dΩ

= Z

Ω

(BiB_k)(x ) (BjB_l)(y ) dΩ

=

Z

Ωx

BiBkdx

Z

Ωy

BjBldy

!

= M^x_ikM^y_jl M = M^x ⊗ M^y (Kronecker product)

16 / 43

(17)

Isogeometric L2 projections

B-spline basis functions have local support (over p + 1 elements) M^x, M^y, . . . – banded structure

M^x_ij = 0 ⇐⇒ |i − j| > 2p + 1

Exemplary basis functions and matrix for cubics







(B₁, B₁)_L2 (B₁, B₂)_L2 (B₁, B₃)_L2 (B₁, B₄)_L2 0 0 · · · 0

(B₂, B₁)_L2 (B₂, B₂)_L2 (B₂, B₃)_L2 (B₂, B₄)_L2 (B₂, B₅)_L2 0 · · · 0 (B₃, B₁)_L2 (B₃, B₂)_L2 (B₃, B₃)_L2 (B₃, B₄)_L2 (B₃, B₅)_L2 (B₃, B₆)_L2 · · · 0

.. .

.. . 0 0 . . . (Bn, B_n−3)_L2 (Bn, B_n−2)_L2 (Bn, B_n−1)_L2 (Bn, Bn)_L2







17 / 43

(18)

Isogeometric L2 projections

Two steps – solving systems with A and B in different directions







A11 A12 · · · 0 A₂₁ A₂₂ · · · 0 ... ... . .. ... 0 0 · · · A_nn













y11 y21 · · · y_m1 y₁₂ y₂₂ · · · y_m1 ... ... . .. ... y_1n y_2n · · · y_mn







=







b11 b21 · · · b_m1 b₁₂ b₂₂ · · · b_m2 ... ... . .. ... b_1n b_2n · · · b_mn













B11 B12 · · · 0 B21 B22 · · · 0 ... ... . .. ... 0 0 · · · B_mm













x11 · · · x1n

x21 · · · x2n

... . .. ... xm1 · · · x_mn







=







y11 y12 · · · y1n

y21 y22 · · · y2n

... ... . .. ... ym1 ym2 · · · y_mn







Two 1D problems with multiple RHS, linear cost O(N) n × n with m right hand sides → O(n ∗ m) = O(N) m × m with n right hand sides → O(m ∗ n) = O(N)

18 / 43

(19)

Derivation of Spatial Direction Splitting

Idea exploit Kronecker product structure of M = M^x⊗ M^y Generally, consider

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 Bm1 A Bm2 · · · A B_mm







19 / 43

(20)

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

20 / 43

(21)

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

21 / 43

(22)

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ⁱ )

22 / 43

(23)

Explicit dynamics

Applications to time-dependent problems (Fortran sequential) M. Łoś, M. Woźniak, M. Paszyński, L. Dalcin, V.M. Calo, Dynamics with Matrices Possessing Kronecker Product Structure, Procedia Computer Science 51 (2015) 286-295

In general: non-stationary problem of the form

∂_tu − L(u) = f (x , t) with some initial state u₀ and boundary conditions L – well-posed linear spatial partial differential operator Discretization:

spatial discretization: isogeometric FEM Basis functions: tensor product B-splines u(x , y ) ≈^P_{i ,j}u_{i ,j}B_{i ,p}^x (x )B^y_j,p(y )

23 / 43

(24)

Explicit dynamics

spatial discretization: isogeometric FEM Basis functions: tensor product B-splines u(x , y ) ≈^P_{i ,j}u_k,lB_{i ,p}^x (x )B_j,p^y (y )

time discretization with explicit method

ut+1−ut

dt = Lu_t+ f_t → u_t+1= u_t+ dtLu_t

implies isogeometric L2 projections in every time step (ut+1, v )L2= (ut+ dtLut, v )L2

24 / 43

(25)

Explicit dynamics

implies isogeometric L2 projections in every time step (u_t+1, v )_L2= (u_t+ dtLu_t, v )_L2

ut+1≈^P_{i ,j}u^{i ,j}_t+1B_{i ,p}^x (x )B_j,p^y (y ), v ← B_k,p^x (x )B_{l ,p}^y (y ) u_t≈^P_{i ,j}u_t^{i ,j}B_{i ,p}^x (x )B_j,p^y (y ))

so the system looks like P

i ,ju_t+1^{i ,j} (B_{i ,p}^x (x )B_j,p^y (y ), B_k,p^x (x )B_{l ,p}^y (y ))_L2 = P

i ,ju_t^{i ,j}(B_{i ,p}^x (x )B_j,p^y (y )) +

dt^P_{i ,j}u_t^{i ,j}L(B_{i ,p}^x (x )B_j,p^y (y )), v )_L2 ∀k, l

25 / 43

(26)

Explicit dynamics

sequence of isogeometric L2 projections P

i ,ju^{i ,j}_t+1(B^x_{i ,p}(x )B_j,p^y (y ), B^x_k,p(x )B_{l ,p}^y (y ))_L2=

P

i ,ju^{i ,j}_t (B_{i ,p}^x (x )B^y_j,p(y )) + dtP

i ,ju^{i ,j}_t L(B_{i ,p}^x (x )B_j,p^y (y )), B^x_k,pB^y l ,p)_L2 ∀k, l







(B_1,p^x B^y 1,p, B^x_1,pB^y

1,p)_L2 (B_1,p^x B^y 1,p, B_2,p^x B^y

1,p)_L2 · · · (B^x_1,pB^y 1,p, B^x

Nx ,pB^y Ny ,p)_L2 (B_2,p^x B_1,p^y , B^x_1,pB_1,p^y )_L2 (B_2,p^x B_1,p^y , B_2,p^x B_1,p^y )_L2 · · · (B^x_2,pB^y_1,p, B^x

Nx ,pB^y Ny ,p)_L2 ..

.

.. .

.. . (B^x

Nx ,pB^y

Ny ,p, B_1,p^x B_1,p^y )_L2 (B^x Nx ,pB^y

Ny ,p, B_2,p^x B_1,p^y )_L2 · · · (B^x Nx ,pB^y

Ny ,p, B^x Nx ,pB^y

Ny ,p)_L2













u^1,1 t+1 u_t+1^2,1 .. . u_t+1^{Nx ,Ny}







=





 P

i ,ju_t^{i ,j}(B_{i ,p}^x (x )B^y

j,p(y )) + dtP

i ,ju^{i ,j}_t L(B^x_{i ,p}(x )B^y

j,p(y )), B_1,p^x B^y 1,p)_L2

P

i ,ju_t^{i ,j}(B_{i ,p}^x (x )B_j,p^y (y )) + dtP

i ,ju^{i ,j}_t L(B^x_{i ,p}(x )B_j,p^y (y )), B_2,p^x B_1,p^y )_L2 ..

.

P

i ,ju_t^{i ,j}(B_{i ,p}^x (x )B^y_j,p(y )) + dtP

i ,ju^{i ,j}_t L(B^x_{i ,p}(x )B_j,p^y (y )), B^x_{Nx ,p}B^y Ny ,p)_L2







26 / 43

(27)

Example 1: Heat transfer equation

We seek the temperature scalar field u : Ω → R such as:











∂u

∂t = ∆u + f (x) on Ω × [0, T ]

∇u · ˆn = 0 on ∂ Ω × [0, T ] u(x, 0) = u₀(x) on Ω

(1)

where Ω = [0, 1]², ˆ

n is a normal vector of the domain boundary, T is a length of the time interval for the simulation, and u₀ is an initial state.

f = 0 (no heat source)

27 / 43

(28)

Example 1: Heat transfer equation

The corresponding weak formulation is obtained by multiplying (1) by test function w ∈ H¹(Ω), integrating by parts over Ω, and imposing the boundary conditions.

Find u ∈ C¹ [0, T ] , H¹(Ω)such that for each t ∈ [0, T ] (∂u

∂t, w )_L2 = −(∇u, ∇w )_L2+ (f , w )_L2 ∀w ∈ H¹(Ω) (2) where (·, ·)_L2 stands for the L²(Ω) scalar product

We utilize Euler time integration scheme

(u_t+1, w )_L2 = (u_t, w )_L₂− dt ∗ ∇u_t, ∇w )_L2 ∀w ∈ H¹(Ω) (3)

28 / 43

(29)

Example 1: Heat transfer equation

Click in the middle

29 / 43

(30)

Code for Example 1 (Heat transfer equation)

"problems/heat/heat_3d.cpp"

#include "problems/heat/heat_3d.hpp"

using namespace ads;

using namespace ads::problems;

pilot for the simulation int main() {

quadratic B-splnes, 12 elements along axis dim_config dim{ 2, 12 };

5000 time steps, time step size 10⁻⁷

timesteps_config steps{ 5000, 1e-7 };

we will need to compute first derivatives during the computations int ders = 1;

some auxiliary objects for configuration and simulation config_3d c{dim, dim, dim, steps, ders};

heat_3d sim{c};

run the simulation sim.run();

}

30 / 43

(31)

Code for Example 1 (Heat transfer equation)

"problems/heat/heat_3d.hpp"

#include "ads/simulation.hpp"

using namespace ads;

using namespace problems;

class heat_3d : public simulation_3d { ...

implementation of the initial state

double init_state(double x, double y, double z) executed once before the simulation starts

void before() override

executed before every simulation step void before_step() override implementation of the simulation step void step() override

executed after every simulation step void after_step() override implementation of generation of RHS void compute_rhs() override executed once after the simulation ends

void after() override ^{31 / 43}

(32)

Code for Example 1 (Heat transfer equation)

this functions is called from before at the beginning of the simulation

the function returns the value of u0 = u(x , y , z)|t=0) computed at point (x , y , z)

double init_state(double x, double y, double z) { double dx = x - 0.5;

double dy = y - 0.5;

double dz = z - 0.5;

double r2 = std::min(8*(dx*dx+dy*dy+dz*dz),1.0);

return (r2 - 1) * (r2 - 1) * (r2 + 1) * (r2 + 1);

};

32 / 43

(33)

Code for Example 1 (Heat transfer equation)

this function is called once before the simulation starts void before() override {

performs LU factorization of three 1D systems, representing B-splines along x , y and z axes

prepare_matrices();

pointer to init_state function

auto init = [this](double x, double y, double z) { return init_state(x, y, z); };

preparation of the initial state projection(u, init);

forward and backward substitutions with multiple RHS solve(u);

}

33 / 43

(34)

Code for Example 1 (Heat transfer equation)

this function is called before every time step

void before_step(int /*iter*/, double /*t*/) override {

using std::swap;

swap u_t and u_t−1 swap(u, u_prev);

}

this function implements every time step

void step(int /*iter*/, double /*t*/) override { generate new RHS using u_prev

compute_rhs();

forward and backward substitutions with multiple RHS solve(u);

}

34 / 43

(35)

Example 1: Heat transfer equation

(u_t+1, w )_L2 = (u_t, w )_L₂− dt ∗ (∇u_t, ∇w )_L2 ∀w ∈ H¹(Ω) (4)

value of test function a over element e at Gauss point q value_type v = eval_basis(e, q, a);

value of ut at Gauss point

value_type u = eval_fun(u_prev, e, q);

computations of double gradient

double gradient = u.dx*v.dx+u.dy*v.dy+u.dz*v.dz;

RHS = u_t− dt∇u_t· ∇v

double val = u.val*v.val - steps.dt * gradient;

scale by Jacobian and weight

rhs(a[0],a[1],a[2])+=val*w*J;

35 / 43

(36)

Code for Example 1 (Heat transfer equation)

void compute_rhs() { auto& rhs = u; zero(rhs);

for (auto e : elements()) { loop through elements double J = jacobian(e);compute Jacobian

for (auto q:quad_points()){ loop through Gauss points double w = weigth(q); Gauss weight

for (auto a : dofs_on_element(e)){loop through dofs value of basis function q over element e at Gauss point q

value_type v = eval_basis(e, q, a);

value of u_t at Gauss point

this also computes derivatives and stored at *.dx value_type u = eval_fun(u_prev, e, q);

computations of double gradient

double gradient = u.dx*v.dx+u.dy*v.dy+u.dz*v.dz;

RHS = ut− dt∇u · ∇v

double val = u.val*v.val - steps.dt * gradient;

scale by Jacobian and weight

rhs(a[0],a[1],a[2])+=val*w*J;

} } } }

36 / 43

(37)

Example 2: Non-linear flow in heterogenous media

Hydraulic fracturing - oil/gas extraction technique consisting in high-pressure fluid injection into the deposit

37 / 43

(38)

Example 2: Non-linear flow in heterogenous media

Hydraulic fracturing - oil/gas extraction technique consisting in high-pressure fluid injection into the deposit

Spatial domain = Ω = [0, 1]³











∂u

∂t − ∇ · (κ(x, u) ∇u) = h(x, t) in Ω × [0, T ]

∇u · ˆn = 0 on ∂ Ω × [0, T ] u(x , 0) = u₀ in Ω

u – pressure

zero Neumann boundary conditions initial state u₀

κ – permeability

h – forcing (induced by extraction method)

M. Alotaibi, V.M. Calo, Y. Efendiev, J. Galvis, M. Ghommem, Global-Local Nonlinear Model Reduction for Flows in Heterogeneous Porous Media arXiv:1407.0782 [math.NA]

38 / 43

(39)

Example 2: Non-linear flow in heterogenous media

κ(x, u) = Kq(x ) b(u) b(u) = e^µu

µ = 10

Kq(x) – property of the terrain (example below)

39 / 43

(40)

Example 2: Non-linear flow in heterogenous media

Extraction process modeled by pumps and sinks pump/sink has a location x ∈ Ω

pumps locally increase the pressure u

sinks locally decrease u (the higher, the faster) h(x , t) = ^X

p∈P

φ (kx_p− x k)−^X

s∈S

u(x , t)φ (kx_s− x k)

P, S – sets of pump and sinks xp, xs – location of pump p/sink s φ – cut-off function (r = 0.15)

φ(t) = ( _t

r − 1² ^t_r + 1² for t ≤ r

0 for t > r

0 r = 0.15 0.4

0 1

40 / 43

(41)

Example 2: Non-linear flow in heterogenous media

Initial state is derived from the permeability of the material K_q K˜q(x) = (Kq(x) − 1)/(1000 − 1)

u0(x) = 0.1 ˜Kq(x) θ0.2,0.3(kx − ck) c = (0.5, 0.5, 0.5)

0 r = 0.2 R = 0.3 0.4

0 1

Figure:θr ,R

41 / 43

(42)

Example 2: Non-linear flow in heterogenous media

We utilize Euler time integration scheme

(u_t+1, w )_L2 = (u_t−dt∗K_q(x ) e^10∗u^t, u_t)+(∇u_t+h(u_t), ∇w )_L2 ∀w ∈ H¹(Ω) where Kq(x , t) does not change with time, and it is given by the

permeability map,

h(x , t) are pumps and sinks h(x , t) = ^X

p∈P

φ (kx_p− x k) −^X

s∈S

u(x , t)φ (kx_s− x k)

42 / 43

(43)

Example 2: Non-linear flow in heterogenous media

Click in the middle

43 / 43