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

## 129 statements

, created at 2024-09-20 17:10 +0000

1import sys

3import numpy as np

4import scipy.linalg as la

5import scipy.sparse as sp

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

15 elif order == 2:

16 stencil = [1.0, -4.0, 3.0]

17 coeff = 1.0 / 2.0

18 zero_pos = 3

20 elif order == 3:

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

22 coeff = 1.0 / 6.0

23 zero_pos = 3

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

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')

37 first_col = np.zeros(N)

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)

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)

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

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

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

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')

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 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]

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

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

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

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

128 A *= coeff * (1.0 / dx)

130 return sp.csc_matrix(A)

133#

134#

135#

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

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

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)

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

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

161 return coeff * (1.0 / dx) * b

164#

165#

166#

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

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

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)

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

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

192 return coeff * (1.0 / dx) * b