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

1import scipy.sparse as sp 

2 

3from pySDC.helpers.problem_helper import ( 

4 get_finite_difference_matrix, 

5 get_finite_difference_stencil, 

6) 

7 

8 

9def getUpwindMatrix(N, dx, order): 

10 A, _ = get_finite_difference_matrix( 

11 derivative=1, order=order, stencil_type='upwind', dx=dx, size=N, dim=1, bc='periodic' 

12 ) 

13 return sp.csc_matrix(A) 

14 

15 

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

17 r""" 

18 Centered first derivative matrix for homogeneous boundary conditions. 

19 

20 For periodic boundaries this is just the generic finite difference matrix. For Neumann and 

21 Dirichlet boundaries only the interior stencil is taken from the generic helper, while the 

22 rows next to the boundary keep the closures that were derived for this problem. 

23 

24 These closures are *not* interchangeable with the ones the generic helper builds, even though 

25 the generic ones are more accurate: the Boussinesq operator differentiates the pressure with 

26 the Neumann closure and the vertical velocity with the Dirichlet closure, and neutral stability 

27 of the resulting wave operator depends on how those two closures relate to each other. Building 

28 each of them in isolation, as the generic helper does, gives a discretisation whose spectrum has 

29 positive real parts and which therefore grows exponentially in time. See #233. 

30 """ 

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

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

33 

34 if bc_left == 'periodic' or bc_right == 'periodic': 

35 assert bc_left == bc_right, "Periodic BC can only be selected for both sides simultaneously" 

36 A, _ = get_finite_difference_matrix( 

37 derivative=1, order=order, stencil_type='center', dx=dx, size=N, dim=1, bc='periodic' 

38 ) 

39 return sp.csc_matrix(A) 

40 

41 assert order in [2, 4], "Neumann and Dirichlet closures are only available for order 2 and 4" 

42 

43 coeff, steps = get_finite_difference_stencil(derivative=1, order=order, stencil_type='center') 

44 A = sp.lil_matrix(sp.diags(coeff, steps, shape=(N, N))) 

45 

46 # Neumann boundary conditions 

47 if bc_left == 'neumann': 

48 A[0, :] = 0.0 

49 if order == 2: 

50 A[0, 0] = -2.0 / 3.0 

51 A[0, 1] = 2.0 / 3.0 

52 elif order == 4: 

53 A[0, 0] = -2.0 / 3.0 

54 A[0, 1] = 2.0 / 3.0 

55 A[1, 0] = -5.0 / 9.0 

56 A[1, 1] = -1.0 / 36.0 

57 

58 if bc_right == 'neumann': 

59 A[N - 1, :] = 0.0 

60 if order == 2: 

61 A[N - 1, N - 2] = -2.0 / 3.0 

62 A[N - 1, N - 1] = 2.0 / 3.0 

63 elif order == 4: 

64 A[N - 2, N - 1] = 5.0 / 9.0 

65 A[N - 2, N - 2] = 1.0 / 36.0 

66 A[N - 1, N - 1] = 2.0 / 3.0 

67 A[N - 1, N - 2] = -2.0 / 3.0 

68 

69 # Dirichlet boundary conditions. For order 2 the ghost value drops out of the stencil, so 

70 # there is nothing to do. 

71 if bc_left == 'dirichlet' and order == 4: 

72 A[0, :] = 0.0 

73 A[0, 1] = 1.0 / 2.0 

74 

75 if bc_right == 'dirichlet' and order == 4: 

76 A[N - 1, :] = 0.0 

77 A[N - 1, N - 2] = -1.0 / 2.0 

78 

79 return sp.csc_matrix(A / dx)