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

## 59 statements

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

1import numpy as np

2import scipy.sparse as sp

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

5 getMatrix,

6 getUpwindMatrix,

7 getBCLeft,

8 getBCRight,

9)

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

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'

32 h = np.zeros(2)

33 x = None

34 z = None

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]

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]

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]

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]

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

57 return xx, zz, h

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'

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)

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

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

73 return Dx, Dz

76#

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

78#

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'

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

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

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

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

97 return bl, br

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'

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

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

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

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

116 return bd, bu