Coverage for pySDC/implementations/problem_classes/acoustic_helpers/buildFDMatrix.py: 39%

1import numpy as np 

2import scipy.linalg as la 

3import scipy.sparse as sp 

4 

5 

6# Only for periodic BC because we have advection only in x direction 

7def getHorizontalDx(N, dx, order): 

8 if order == 1: 

9 stencil = [-1.0, 1.0] 

10 zero_pos = 2 

11 coeff = 1.0 

12 

13 elif order == 2: 

14 stencil = [1.0, -4.0, 3.0] 

15 coeff = 1.0 / 2.0 

16 zero_pos = 3 

17 

18 elif order == 3: 

19 stencil = [1.0, -6.0, 3.0, 2.0] 

20 coeff = 1.0 / 6.0 

21 zero_pos = 3 

22 

23 elif order == 4: 

24 stencil = [-5.0, 30.0, -90.0, 50.0, 15.0] 

25 coeff = 1.0 / 60.0 

26 zero_pos = 4 

27 

28 elif order == 5: 

29 stencil = [3.0, -20.0, 60.0, -120.0, 65.0, 12.0] 

30 coeff = 1.0 / 60.0 

31 zero_pos = 5 

32 else: 

33 print("Order " + order + " not implemented.") 

34 

35 first_col = np.zeros(N) 

36 

37 # Because we need to specific first column (not row) in circulant, flip stencil array 

38 first_col[0 : np.size(stencil)] = np.flipud(stencil) 

39 

40 # Circulant shift of coefficient column so that entry number zero_pos becomes first entry 

41 first_col = np.roll(first_col, -np.size(stencil) + zero_pos, axis=0) 

42 

43 return sp.csc_matrix(coeff * (1.0 / dx) * la.circulant(first_col)) 

44 

45 

46def getMatrix(N, dx, bc_left, bc_right, order): 

47 assert bc_left in ['periodic', 'neumann', 'dirichlet'], "Unknown type of BC" 

48 

49 if order == 2: 

50 stencil = [-1.0, 0.0, 1.0] 

51 range = [-1, 0, 1] 

52 coeff = 1.0 / 2.0 

53 elif order == 4: 

54 stencil = [1.0, -8.0, 0.0, 8.0, -1.0] 

55 range = [-2, -1, 0, 1, 2] 

56 coeff = 1.0 / 12.0 

57 elif order == 6: 

58 stencil = [-1.0, 9.0, -45.0, 0.0, 45.0, -9.0, 1.0] 

59 range = [-3, -2, -1, 0, 1, 2, 3] 

60 coeff = 1.0 / 60.0 

61 

62 A = sp.diags(stencil, range, shape=(N, N)) 

63 A = sp.lil_matrix(A) 

64 

65 # 

66 # Periodic boundary conditions 

67 # 

68 if bc_left in ['periodic']: 

69 assert bc_right in ['periodic'], "Periodic BC can only be selected for both sides simultaneously" 

70 

71 if bc_left in ['periodic']: 

72 if order == 2: 

73 A[0, N - 1] = stencil[0] 

74 

75 elif order == 4: 

76 A[0, N - 2] = stencil[0] 

77 A[0, N - 1] = stencil[1] 

78 A[1, N - 1] = stencil[0] 

79 

80 elif order == 6: 

81 A[0, N - 3] = stencil[0] 

82 A[0, N - 2] = stencil[1] 

83 A[0, N - 1] = stencil[2] 

84 A[1, N - 2] = stencil[0] 

85 A[1, N - 1] = stencil[1] 

86 A[2, N - 1] = stencil[0] 

87 

88 if bc_right in ['periodic']: 

89 if order == 2: 

90 A[N - 1, 0] = stencil[2] 

91 elif order == 4: 

92 A[N - 2, 0] = stencil[4] 

93 A[N - 1, 0] = stencil[3] 

94 A[N - 1, 1] = stencil[4] 

95 elif order == 6: 

96 A[N - 3, 0] = stencil[6] 

97 A[N - 2, 0] = stencil[5] 

98 A[N - 2, 1] = stencil[6] 

99 A[N - 1, 0] = stencil[4] 

100 A[N - 1, 1] = stencil[5] 

101 A[N - 1, 2] = stencil[6] 

102 

103 # 

104 # Neumann boundary conditions 

105 # 

106 if bc_left in ['neumann']: 

107 A[0, :] = np.zeros(N) 

108 if order == 2: 

109 A[0, 0] = -4.0 / 3.0 

110 A[0, 1] = 4.0 / 3.0 

111 elif order == 4: 

112 A[0, 0] = -8.0 

113 A[0, 1] = 8.0 

114 A[1, 0] = -8.0 + 4.0 / 3.0 

115 A[1, 1] = -1.0 / 3.0 

116 

117 if bc_right in ['neumann']: 

118 A[N - 1, :] = np.zeros(N) 

119 if order == 2: 

120 A[N - 1, N - 2] = -4.0 / 3.0 

121 A[N - 1, N - 1] = 4.0 / 3.0 

122 elif order == 4: 

123 A[N - 2, N - 1] = 8.0 - 4.0 / 3.0 

124 A[N - 2, N - 2] = 1.0 / 3.0 

125 A[N - 1, N - 1] = 8.0 

126 A[N - 1, N - 2] = -8.0 

127 

128 # 

129 # Dirichlet boundary conditions 

130 # 

131 if bc_left in ['dirichlet']: 

132 # For order==2, nothing to do here 

133 if order == 4: 

134 A[0, :] = np.zeros(N) 

135 A[0, 1] = 6.0 

136 

137 if bc_right in ['dirichlet']: 

138 # For order==2, nothing to do here 

139 if order == 4: 

140 A[N - 1, :] = np.zeros(N) 

141 A[N - 1, N - 2] = -6.0 

142 

143 A = coeff * (1.0 / dx) * A 

144 return sp.csc_matrix(A) 

145 

146 

147# 

148# 

149# 

150def getBCLeft(value, N, dx, type, order): 

151 assert type in ['periodic', 'neumann', 'dirichlet'], "Unknown type of BC" 

152 

153 if order == 2: 

154 coeff = 1.0 / 2.0 

155 elif order == 4: 

156 coeff = 1.0 / 12.0 

157 

158 b = np.zeros(N) 

159 if type in ['dirichlet']: 

160 if order == 2: 

161 b[0] = -value 

162 elif order == 4: 

163 b[0] = -6.0 * value 

164 b[1] = 1.0 * value 

165 

166 if type in ['neumann']: 

167 if order == 2: 

168 b[0] = (2.0 / 3.0) * dx * value 

169 elif order == 4: 

170 b[0] = 4.0 * dx * value 

171 b[1] = -(2.0 / 3.0) * dx * value 

172 

173 return coeff * (1.0 / dx) * b 

174 

175 

176# 

177# 

178# 

179def getBCRight(value, N, dx, type, order): 

180 assert type in ['periodic', 'neumann', 'dirichlet'], "Unknown type of BC" 

181 

182 if order == 2: 

183 coeff = 1.0 / 2.0 

184 elif order == 4: 

185 coeff = 1.0 / 12.0 

186 

187 b = np.zeros(N) 

188 if type in ['dirichlet']: 

189 if order == 2: 

190 b[N - 1] = value 

191 elif order == 4: 

192 b[N - 2] = -1.0 * value 

193 b[N - 1] = 6.0 * value 

194 

195 if type in ['neumann']: 

196 if order == 2: 

197 b[N - 1] = (2.0 / 3.0) * dx * value 

198 elif order == 4: 

199 b[N - 2] = -(2.0 / 3.0) * dx * value 

200 b[N - 1] = 4.0 * dx * value 

201 

202 return coeff * (1.0 / dx) * b