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

## 145 statements

, created at 2024-09-09 14:59 +0000

1import numpy as np

2import scipy.linalg as la

3import scipy.sparse as sp

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

13 elif order == 2:

14 stencil = [1.0, -4.0, 3.0]

15 coeff = 1.0 / 2.0

16 zero_pos = 3

18 elif order == 3:

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

20 coeff = 1.0 / 6.0

21 zero_pos = 3

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

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.")

35 first_col = np.zeros(N)

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)

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)

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

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

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

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

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

63 A = sp.lil_matrix(A)

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"

71 if bc_left in ['periodic']:

72 if order == 2:

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

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]

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]

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]

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

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

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

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

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

144 return sp.csc_matrix(A)

147#

148#

149#

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

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

153 if order == 2:

154 coeff = 1.0 / 2.0

155 elif order == 4:

156 coeff = 1.0 / 12.0

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

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

173 return coeff * (1.0 / dx) * b

176#

177#

178#

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

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

182 if order == 2:

183 coeff = 1.0 / 2.0

184 elif order == 4:

185 coeff = 1.0 / 12.0

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

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

202 return coeff * (1.0 / dx) * b