Coverage for pySDC/implementations/problem_classes/AdvectionEquation_ND_FD.py: 67%

1import numpy as np 

2 

3from pySDC.implementations.problem_classes.generic_ND_FD import GenericNDimFinDiff 

4 

5 

6# noinspection PyUnusedLocal 

7class advectionNd(GenericNDimFinDiff): 

8 r""" 

9 Example implementing the unforced ND advection equation with periodic 

10 or Dirichlet boundary conditions in :math:`[0,1]^N` 

11 

12 .. math:: 

13 \frac{\partial u}{\partial t} = -c \frac{\partial u}{\partial x}, 

14 

15 and initial solution of the form 

16 

17 .. math:: 

18 u({\bf x},0) = \prod_{i=1}^N \sin(f\pi x_i), 

19 

20 with :math:`x_i` the coordinate in :math:`i^{th}` dimension. 

21 Discretization uses central finite differences. 

22 

23 Parameters 

24 ---------- 

25 nvars : int of tuple, optional 

26 Spatial resolution (same in all dimensions). Using a tuple allows to 

27 consider several dimensions, e.g ``nvars=(16,16)`` for a 2D problem. 

28 c : float, optional 

29 Advection speed (same in all dimensions). 

30 freq : int of tuple, optional 

31 Spatial frequency :math:`f` of the initial conditions, can be tuple. 

32 stencil_type : str, optional 

33 Type of the finite difference stencil. 

34 order : int, optional 

35 Order of the finite difference discretization. 

36 lintol : float, optional 

37 Tolerance for spatial solver (GMRES). 

38 liniter : int, optional 

39 Max. iterations number for GMRES. 

40 solver_type : str, optional 

41 Solve the linear system directly or using GMRES or CG 

42 bc : str, optional 

43 Boundary conditions, either ``'periodic'`` or ``'dirichlet'``. 

44 sigma : float, optional 

45 If ``freq=-1`` and ``ndim=1``, uses a Gaussian initial solution of the form 

46 

47 .. math:: 

48 u(x,0) = e^{ 

49 \frac{\displaystyle 1}{\displaystyle 2} 

50 \left( 

51 \frac{\displaystyle x-1/2}{\displaystyle \sigma} 

52 \right)^2 

53 } 

54 

55 Attributes 

56 ---------- 

57 A : sparse matrix (CSC) 

58 FD discretization matrix of the ND grad operator. 

59 Id : sparse matrix (CSC) 

60 Identity matrix of the same dimension as A. 

61 

62 Note 

63 ---- 

64 Args can be set as values or as tuples, which will increase the dimension. 

65 Do, however, take care that all spatial parameters have the same dimension. 

66 """ 

67 

68 def __init__( 

69 self, 

70 nvars=512, 

71 c=1.0, 

72 freq=2, 

73 stencil_type='center', 

74 order=2, 

75 lintol=1e-12, 

76 liniter=10000, 

77 solver_type='direct', 

78 bc='periodic', 

79 sigma=6e-2, 

80 ): 

81 super().__init__(nvars, -c, 1, freq, stencil_type, order, lintol, liniter, solver_type, bc) 

82 

83 if solver_type == 'CG': # pragma: no cover 

84 self.logger.warning('CG is not usually used for advection equation') 

85 self._makeAttributeAndRegister('c', localVars=locals(), readOnly=True) 

86 self._makeAttributeAndRegister('sigma', localVars=locals()) 

87 

88 def u_exact(self, t, **kwargs): 

89 r""" 

90 Routine to compute the exact solution at time :math:`t`. 

91 

92 Parameters 

93 ---------- 

94 t : float 

95 Time of the exact solution. 

96 **kwargs : dict 

97 Additional arguments (that won't be used). 

98 

99 Returns 

100 ------- 

101 sol : dtype_u 

102 The exact solution. 

103 """ 

104 if 'u_init' in kwargs.keys() or 't_init' in kwargs.keys(): 

105 self.logger.warning( 

106 f'{type(self).__name__} uses an analytic exact solution from t=0. If you try to compute the local error, you will get the global error instead!' 

107 ) 

108 

109 # Initialize pointers and variables 

110 ndim, freq, c, sigma, sol = self.ndim, self.freq, self.c, self.sigma, self.u_init 

111 

112 if ndim == 1: 

113 x = self.grids 

114 if freq[0] >= 0: 

115 sol[:] = np.sin(np.pi * freq[0] * (x - c * t)) 

116 elif freq[0] == -1: 

117 # Gaussian initial solution 

118 sol[:] = np.exp(-0.5 * (((x - (c * t)) % 1.0 - 0.5) / sigma) ** 2) 

119 

120 elif ndim == 2: 

121 x, y = self.grids 

122 sol[:] = np.sin(np.pi * freq[0] * (x - c * t)) * np.sin(np.pi * freq[1] * (y - c * t)) 

123 

124 elif ndim == 3: 

125 x, y, z = self.grids 

126 sol[:] = ( 

127 np.sin(np.pi * freq[0] * (x - c * t)) 

128 * np.sin(np.pi * freq[1] * (y - c * t)) 

129 * np.sin(np.pi * freq[2] * (z - c * t)) 

130 ) 

131 

132 return sol