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

1import sys 

2 

3import numpy as np 

4import scipy.linalg as la 

5import scipy.sparse as sp 

6 

7 

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

9def getUpwindMatrix(N, dx, order): 

10 if order == 1: 

11 stencil = [-1.0, 1.0] 

12 zero_pos = 2 

13 coeff = 1.0 

14 

15 elif order == 2: 

16 stencil = [1.0, -4.0, 3.0] 

17 coeff = 1.0 / 2.0 

18 zero_pos = 3 

19 

20 elif order == 3: 

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

22 coeff = 1.0 / 6.0 

23 zero_pos = 3 

24 

25 elif order == 4: 

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

27 coeff = 1.0 / 60.0 

28 zero_pos = 4 

29 

30 elif order == 5: 

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

32 coeff = 1.0 / 60.0 

33 zero_pos = 5 

34 else: 

35 sys.exit('Order ' + str(order) + ' not implemented') 

36 

37 first_col = np.zeros(N) 

38 

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

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

41 

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

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

44 

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

46 

47 

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

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

50 

51 if order == 2: 

52 stencil = [-1.0, 0.0, 1.0] 

53 range = [-1, 0, 1] 

54 coeff = 1.0 / 2.0 

55 elif order == 4: 

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

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

58 coeff = 1.0 / 12.0 

59 else: 

60 sys.exit('Order ' + str(order) + ' not implemented') 

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 if bc_right in ['periodic']: 

81 if order == 2: 

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

83 elif order == 4: 

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

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

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

87 

88 # 

89 # Neumann boundary conditions 

90 # 

91 if bc_left in ['neumann']: 

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

93 if order == 2: 

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

95 A[0, 1] = 4.0 / 3.0 

96 elif order == 4: 

97 A[0, 0] = -8.0 

98 A[0, 1] = 8.0 

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

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

101 

102 if bc_right in ['neumann']: 

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

104 if order == 2: 

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

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

107 elif order == 4: 

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

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

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

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

112 

113 # 

114 # Dirichlet boundary conditions 

115 # 

116 if bc_left in ['dirichlet']: 

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

118 if order == 4: 

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

120 A[0, 1] = 6.0 

121 

122 if bc_right in ['dirichlet']: 

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

124 if order == 4: 

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

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

127 

128 A *= coeff * (1.0 / dx) 

129 

130 return sp.csc_matrix(A) 

131 

132 

133# 

134# 

135# 

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

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

138 

139 if order == 2: 

140 coeff = 1.0 / 2.0 

141 elif order == 4: 

142 coeff = 1.0 / 12.0 

143 else: 

144 raise NotImplementedError('wrong order, got %s' % order) 

145 

146 b = np.zeros(N) 

147 if type in ['dirichlet']: 

148 if order == 2: 

149 b[0] = -value 

150 elif order == 4: 

151 b[0] = -6.0 * value 

152 b[1] = 1.0 * value 

153 

154 if type in ['neumann']: 

155 if order == 2: 

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

157 elif order == 4: 

158 b[0] = 4.0 * dx * value 

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

160 

161 return coeff * (1.0 / dx) * b 

162 

163 

164# 

165# 

166# 

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

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

169 

170 if order == 2: 

171 coeff = 1.0 / 2.0 

172 elif order == 4: 

173 coeff = 1.0 / 12.0 

174 else: 

175 raise NotImplementedError('wrong order, got %s' % order) 

176 

177 b = np.zeros(N) 

178 if type in ['dirichlet']: 

179 if order == 2: 

180 b[N - 1] = value 

181 elif order == 4: 

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

183 b[N - 1] = 6.0 * value 

184 

185 if type in ['neumann']: 

186 if order == 2: 

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

188 elif order == 4: 

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

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

191 

192 return coeff * (1.0 / dx) * b