Coverage for pySDC/implementations/problem_classes/boussinesq_helpers/build2DFDMatrix.py: 0%

1import numpy as np 

2import scipy.sparse as sp 

3 

4from pySDC.implementations.problem_classes.boussinesq_helpers.buildFDMatrix import ( 

5 getMatrix, 

6 getUpwindMatrix, 

7 getBCLeft, 

8 getBCRight, 

9) 

10 

11 

12# 

13# 

14# 

15def get2DUpwindMatrix(N, dx, order): 

16 Dx = getUpwindMatrix(N[0], dx, order) 

17 return sp.kron(Dx, sp.eye(N[1]), format="csr") 

18 

19 

20# 

21# 

22# 

23def get2DMesh(N, x_b, z_b, bc_hor, bc_ver): 

24 assert np.size(N) == 2, 'N needs to be an array with two entries: N[0]=Nx and N[1]=Nz' 

25 assert ( 

26 np.size(x_b) == 2 

27 ), 'x_b needs to be an array with two entries: x_b[0] = left boundary, x_b[1] = right boundary' 

28 assert ( 

29 np.size(z_b) == 2 

30 ), 'z_b needs to be an array with two entries: z_b[0] = lower boundary, z_b[1] = upper boundary' 

31 

32 h = np.zeros(2) 

33 x = None 

34 z = None 

35 

36 if bc_hor[0] in ['periodic']: 

37 assert bc_hor[1] in ['periodic'], 'Periodic boundary conditions must be prescribed at both boundaries' 

38 x = np.linspace(x_b[0], x_b[1], N[0], endpoint=False) 

39 h[0] = x[1] - x[0] 

40 

41 if bc_hor[0] in ['dirichlet', 'neumann']: 

42 x = np.linspace(x_b[0], x_b[1], N[0] + 2, endpoint=True) 

43 x = x[1 : N[0] + 1] 

44 h[0] = x[1] - x[0] 

45 

46 if bc_ver[0] in ['periodic']: 

47 assert bc_ver[1] in ['periodic'], 'Periodic boundary conditions must be prescribed at both boundaries' 

48 z = np.linspace(z_b[0], z_b[1], N[1], endpoint=False) 

49 h[1] = z[1] - z[0] 

50 

51 if bc_ver[0] in ['dirichlet', 'neumann']: 

52 z = np.linspace(z_b[0], z_b[1], N[1] + 2, endpoint=True) 

53 z = z[1 : N[1] + 1] 

54 h[1] = z[1] - z[0] 

55 

56 xx, zz = np.meshgrid(x, z, indexing="ij") 

57 return xx, zz, h 

58 

59 

60# 

61# 

62# 

63def get2DMatrix(N, h, bc_hor, bc_ver, order): 

64 assert np.size(N) == 2, 'N needs to be an array with two entries: N[0]=Nx and N[1]=Nz' 

65 assert np.size(h) == 2, 'h needs to be an array with two entries: h[0]=dx and h[1]=dz' 

66 

67 Ax = getMatrix(N[0], h[0], bc_hor[0], bc_hor[1], order) 

68 Az = getMatrix(N[1], h[1], bc_ver[0], bc_ver[1], order) 

69 

70 Dx = sp.kron(Ax, sp.eye(N[1]), format="csr") 

71 Dz = sp.kron(sp.eye(N[0]), Az, format="csr") 

72 

73 return Dx, Dz 

74 

75 

76# 

77# NOTE: So far only constant dirichlet values can be prescribed, i.e. one fixed value for a whole segment 

78# 

79 

80 

81def getBCHorizontal(value, N, dx, bc_hor): 

82 assert ( 

83 np.size(value) == 2 

84 ), 'Value needs to be an array with two entries: value[0] for the left and value[1] for the right boundary' 

85 assert np.size(N) == 2, 'N needs to be an array with two entries: N[0]=Nx and N[1]=Nz' 

86 assert np.size(dx) == 1, 'dx must be a scalar' 

87 assert ( 

88 np.size(bc_hor) == 2 

89 ), 'bc_hor must have two entries, bc_hor[0] specifying the BC at the left, bc_hor[1] at the right boundary' 

90 

91 bl = getBCLeft(value[0], N[0], dx, bc_hor[0]) 

92 bl = np.kron(bl, np.ones(N[1])) 

93 

94 br = getBCRight(value[1], N[0], dx, bc_hor[1]) 

95 br = np.kron(br, np.ones(N[1])) 

96 

97 return bl, br 

98 

99 

100def getBCVertical(value, N, dz, bc_ver): 

101 assert ( 

102 np.size(value) == 2 

103 ), 'Value needs to be an array with two entries: value[0] for the left and value[1] for the right boundary' 

104 assert np.size(N) == 2, 'N needs to be an array with two entries: N[0]=Nx and N[1]=Nz' 

105 assert np.size(dz) == 1, 'dx must be a scalar' 

106 assert ( 

107 np.size(bc_ver) == 2 

108 ), 'bc_hor must have two entries, bc_hor[0] specifying the BC at the left, bc_hor[1] at the right boundary' 

109 

110 bd = getBCLeft(value[0], N[1], dz, bc_ver[0]) 

111 bd = np.kron(np.ones(N[0]), bd) 

112 

113 bu = getBCRight(value[1], N[1], dz, bc_ver[1]) 

114 bu = np.kron(np.ones(N[0]), bu) 

115 

116 return bd, bu