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

1from pySDC.implementations.problem_classes.generic_ND_FD import GenericNDimFinDiff 

2 

3 

4# noinspection PyUnusedLocal 

5class advectionNd(GenericNDimFinDiff): 

6 r""" 

7 Example implementing the unforced ND advection equation with periodic 

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

9 

10 .. math:: 

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

12 

13 and initial solution of the form 

14 

15 .. math:: 

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

17 

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

19 Discretization uses central finite differences. 

20 

21 Parameters 

22 ---------- 

23 nvars : int of tuple, optional 

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

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

26 c : float, optional 

27 Advection speed (same in all dimensions). 

28 freq : int of tuple, optional 

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

30 stencil_type : str, optional 

31 Type of the finite difference stencil. 

32 order : int, optional 

33 Order of the finite difference discretization. 

34 lintol : float, optional 

35 Tolerance for spatial solver (GMRES). 

36 liniter : int, optional 

37 Max. iterations number for GMRES. 

38 solver_type : str, optional 

39 Solve the linear system directly or using GMRES or CG 

40 bc : str, optional 

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

42 dtype : dtype-like, optional 

43 Precision the state is stored at, ``float64`` by default. See :class:`GenericNDimFinDiff`. 

44 useGPU : bool, optional 

45 Run on the GPU with CuPy instead of on the CPU with NumPy. 

46 solve_precision : dtype-like or None, optional 

47 Precision the implicit solve runs at, independently of ``dtype``. See :class:`GenericNDimFinDiff`. 

48 sigma : float, optional 

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

50 

51 .. math:: 

52 u(x,0) = e^{ 

53 \frac{\displaystyle 1}{\displaystyle 2} 

54 \left( 

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

56 \right)^2 

57 } 

58 

59 Attributes 

60 ---------- 

61 A : sparse matrix (CSC) 

62 FD discretization matrix of the ND grad operator. 

63 Id : sparse matrix (CSC) 

64 Identity matrix of the same dimension as A. 

65 

66 Note 

67 ---- 

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

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

70 """ 

71 

72 def __init__( 

73 self, 

74 nvars=512, 

75 c=1.0, 

76 freq=2, 

77 stencil_type='center', 

78 order=2, 

79 lintol=1e-12, 

80 liniter=10000, 

81 solver_type='direct', 

82 bc='periodic', 

83 sigma=6e-2, 

84 dtype='float64', 

85 useGPU=False, 

86 solve_precision=None, 

87 ): 

88 """Initialization routine""" 

89 super().__init__( 

90 nvars, 

91 -c, 

92 1, 

93 freq, 

94 stencil_type, 

95 order, 

96 lintol, 

97 liniter, 

98 solver_type, 

99 bc, 

100 dtype=dtype, 

101 useGPU=useGPU, 

102 solve_precision=solve_precision, 

103 ) 

104 

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

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

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

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

109 

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

111 r""" 

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

113 

114 Parameters 

115 ---------- 

116 t : float 

117 Time of the exact solution. 

118 **kwargs : dict 

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

120 

121 Returns 

122 ------- 

123 sol : dtype_u 

124 The exact solution. 

125 """ 

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

127 self.logger.warning( 

128 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!' 

129 ) 

130 

131 # Initialize pointers and variables 

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

133 

134 if ndim == 1: 

135 x = self.grids 

136 if freq[0] >= 0: 

137 sol[:] = self.xp.sin(self.xp.pi * freq[0] * (x - c * t)) 

138 elif freq[0] == -1: 

139 # Gaussian initial solution 

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

141 

142 elif ndim == 2: 

143 x, y = self.grids 

144 sol[:] = self.xp.sin(self.xp.pi * freq[0] * (x - c * t)) * self.xp.sin(self.xp.pi * freq[1] * (y - c * t)) 

145 

146 elif ndim == 3: 

147 x, y, z = self.grids 

148 sol[:] = ( 

149 self.xp.sin(self.xp.pi * freq[0] * (x - c * t)) 

150 * self.xp.sin(self.xp.pi * freq[1] * (y - c * t)) 

151 * self.xp.sin(self.xp.pi * freq[2] * (z - c * t)) 

152 ) 

153 

154 return sol