|
1 | 1 | import numpy as np |
2 | 2 | import scipy.sparse as sp |
3 | 3 | from pySDC.implementations.problem_classes.boussinesq_helpers.build2DFDMatrix import get2DMatrix, getBCHorizontal, \ |
4 | | - getBCVertical, get2DUpwindMatrix |
| 4 | + get2DUpwindMatrix |
| 5 | + |
5 | 6 |
|
6 | 7 | def getBoussinesq2DUpwindMatrix(N, dx, u_adv, order): |
| 8 | + Dx = get2DUpwindMatrix(N, dx, order) |
| 9 | + |
| 10 | + # Note: In the equations it is u_t + u_adv* D_x u = ... so in order to comply with the form u_t = M u, |
| 11 | + # add a minus sign in front of u_adv |
| 12 | + |
| 13 | + Zero = np.zeros((N[0] * N[1], N[0] * N[1])) |
| 14 | + M1 = sp.hstack((-u_adv * Dx, Zero, Zero, Zero), format="csr") |
| 15 | + M2 = sp.hstack((Zero, -u_adv * Dx, Zero, Zero), format="csr") |
| 16 | + M3 = sp.hstack((Zero, Zero, -u_adv * Dx, Zero), format="csr") |
| 17 | + M4 = sp.hstack((Zero, Zero, Zero, -u_adv * Dx), format="csr") |
| 18 | + M = sp.vstack((M1, M2, M3, M4), format="csr") |
| 19 | + |
| 20 | + return sp.csc_matrix(M) |
| 21 | + |
7 | 22 |
|
8 | | - Dx = get2DUpwindMatrix(N, dx, order) |
9 | | - |
10 | | - # Note: In the equations it is u_t + u_adv* D_x u = ... so in order to comply with the form u_t = M u, |
11 | | - # add a minus sign in front of u_adv |
12 | | - |
13 | | - Zero = np.zeros((N[0]*N[1], N[0]*N[1])) |
14 | | - M1 = sp.hstack((-u_adv*Dx, Zero, Zero, Zero), format="csr") |
15 | | - M2 = sp.hstack(( Zero, -u_adv*Dx, Zero, Zero), format="csr") |
16 | | - M3 = sp.hstack(( Zero, Zero, -u_adv*Dx, Zero), format="csr") |
17 | | - M4 = sp.hstack(( Zero, Zero, Zero, -u_adv*Dx), format="csr") |
18 | | - M = sp.vstack((M1,M2,M3,M4), format="csr") |
19 | | - |
20 | | - return sp.csc_matrix(M) |
21 | | - |
22 | 23 | def getBoussinesq2DMatrix(N, h, bc_hor, bc_ver, c_s, Nfreq, order): |
| 24 | + Dx_u, Dz_u = get2DMatrix(N, h, bc_hor[0], bc_ver[0], order) |
| 25 | + Dx_w, Dz_w = get2DMatrix(N, h, bc_hor[1], bc_ver[1], order) |
| 26 | + # Dx_b, Dz_b = get2DMatrix(N, h, bc_hor[2], bc_ver[2], order) |
| 27 | + Dx_p, Dz_p = get2DMatrix(N, h, bc_hor[3], bc_ver[3], order) |
23 | 28 |
|
24 | | - Dx_u, Dz_u = get2DMatrix(N, h, bc_hor[0], bc_ver[0], order) |
25 | | - Dx_w, Dz_w = get2DMatrix(N, h, bc_hor[1], bc_ver[1], order) |
26 | | - Dx_b, Dz_b = get2DMatrix(N, h, bc_hor[2], bc_ver[2], order) |
27 | | - Dx_p, Dz_p = get2DMatrix(N, h, bc_hor[3], bc_ver[3], order) |
| 29 | + # Id_N = sp.eye(N[0] * N[1]) |
28 | 30 |
|
29 | | - Id_N = sp.eye(N[0]*N[1]) |
| 31 | + Zero = np.zeros((N[0] * N[1], N[0] * N[1])) |
| 32 | + Id_w = sp.eye(N[0] * N[1]) |
30 | 33 |
|
31 | | - Zero = np.zeros((N[0]*N[1],N[0]*N[1])) |
32 | | - Id_w = sp.eye(N[0]*N[1]) |
| 34 | + # Note: Bring all terms to right hand side, therefore a couple of minus signs |
| 35 | + # are needed |
33 | 36 |
|
34 | | - # Note: Bring all terms to right hand side, therefore a couple of minus signs |
35 | | - # are needed |
| 37 | + M1 = sp.hstack((Zero, Zero, Zero, -Dx_p), format="csr") |
| 38 | + M2 = sp.hstack((Zero, Zero, Id_w, -Dz_p), format="csr") |
| 39 | + M3 = sp.hstack((Zero, -Nfreq ** 2 * Id_w, Zero, Zero), format="csr") |
| 40 | + M4 = sp.hstack((-c_s ** 2 * Dx_u, -c_s ** 2 * Dz_w, Zero, Zero), format="csr") |
| 41 | + M = sp.vstack((M1, M2, M3, M4), format="csr") |
36 | 42 |
|
37 | | - M1 = sp.hstack(( Zero, Zero, Zero, -Dx_p), format="csr") |
38 | | - M2 = sp.hstack(( Zero, Zero, Id_w, -Dz_p), format="csr") |
39 | | - M3 = sp.hstack(( Zero, -Nfreq**2*Id_w, Zero, Zero), format="csr") |
40 | | - M4 = sp.hstack((-c_s**2*Dx_u, -c_s**2*Dz_w, Zero, Zero), format="csr") |
41 | | - M = sp.vstack((M1,M2,M3,M4), format="csr") |
| 43 | + Id = sp.eye(4 * N[0] * N[1]) |
42 | 44 |
|
43 | | - Id = sp.eye(4*N[0]*N[1]) |
| 45 | + return sp.csc_matrix(Id), sp.csc_matrix(M) |
44 | 46 |
|
45 | | - return sp.csc_matrix(Id), sp.csc_matrix(M) |
46 | 47 |
|
47 | 48 | def getBoussinesqBCHorizontal(value, N, dx, bc_hor): |
48 | | - |
49 | | - bu_left, bu_right = getBCHorizontal( value[0], N, dx, bc_hor[0] ) |
50 | | - bw_left, bw_right = getBCHorizontal( value[1], N, dx, bc_hor[1] ) |
51 | | - bb_left, bb_right = getBCHorizontal( value[2], N, dx, bc_hor[2] ) |
52 | | - bp_left, bp_right = getBCHorizontal( value[3], N, dx, bc_hor[3] ) |
53 | | - |
54 | | - b_left = np.concatenate(( bp_left, bp_left, bu_left + bw_left)) |
55 | | - b_right = np.concatenate(( bp_right, bp_right, bu_right + bw_right)) |
56 | | - return b_left, b_right |
| 49 | + bu_left, bu_right = getBCHorizontal(value[0], N, dx, bc_hor[0]) |
| 50 | + bw_left, bw_right = getBCHorizontal(value[1], N, dx, bc_hor[1]) |
| 51 | + # bb_left, bb_right = getBCHorizontal(value[2], N, dx, bc_hor[2]) |
| 52 | + bp_left, bp_right = getBCHorizontal(value[3], N, dx, bc_hor[3]) |
| 53 | + |
| 54 | + b_left = np.concatenate((bp_left, bp_left, bu_left + bw_left)) |
| 55 | + b_right = np.concatenate((bp_right, bp_right, bu_right + bw_right)) |
| 56 | + return b_left, b_right |
| 57 | + |
57 | 58 |
|
58 | 59 | def getBoussinesqBCVertical(): |
59 | | - return 0.0 |
| 60 | + return 0.0 |
0 commit comments