24 statements

, created at 2024-09-20 16:55 +0000

1import numpy as np

3from pySDC.implementations.problem_classes.generic_ND_FD import GenericNDimFinDiff

6# noinspection PyUnusedLocal

8 r"""

9 Example implementing the unforced ND advection equation with periodic

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

12 .. math::

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

15 and initial solution of the form

17 .. math::

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

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

21 Discretization uses central finite differences.

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

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 }

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.

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

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)

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

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

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

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

89 r"""

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

92 Parameters

93 ----------

94 t : float

95 Time of the exact solution.

96 **kwargs : dict

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

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 )

109 # Initialize pointers and variables

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

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)

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

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 )

132 return sol