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

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 sigma : float, optional 

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

44 

45 .. math:: 

46 u(x,0) = e^{ 

47 \frac{\displaystyle 1}{\displaystyle 2} 

48 \left( 

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

50 \right)^2 

51 } 

52 

53 Attributes 

54 ---------- 

55 A : sparse matrix (CSC) 

56 FD discretization matrix of the ND grad operator. 

57 Id : sparse matrix (CSC) 

58 Identity matrix of the same dimension as A. 

59 

60 Note 

61 ---- 

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

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

64 """ 

65 

66 def __init__( 

67 self, 

68 nvars=512, 

69 c=1.0, 

70 freq=2, 

71 stencil_type='center', 

72 order=2, 

73 lintol=1e-12, 

74 liniter=10000, 

75 solver_type='direct', 

76 bc='periodic', 

77 sigma=6e-2, 

78 ): 

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

80 

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

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

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

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

85 

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

87 r""" 

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

89 

90 Parameters 

91 ---------- 

92 t : float 

93 Time of the exact solution. 

94 **kwargs : dict 

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

96 

97 Returns 

98 ------- 

99 sol : dtype_u 

100 The exact solution. 

101 """ 

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

103 self.logger.warning( 

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

105 ) 

106 

107 # Initialize pointers and variables 

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

109 

110 if ndim == 1: 

111 x = self.grids 

112 if freq[0] >= 0: 

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

114 elif freq[0] == -1: 

115 # Gaussian initial solution 

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

117 

118 elif ndim == 2: 

119 x, y = self.grids 

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

121 

122 elif ndim == 3: 

123 x, y, z = self.grids 

124 sol[:] = ( 

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

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

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

128 ) 

129 

130 return sol