Coverage for pySDC/core/collocation.py: 98%

1import logging 

2from typing import Optional, Any 

3import numpy as np 

4from qmat import Q_GENERATORS 

5 

6from pySDC.core.errors import CollocationError 

7 

8 

9class CollBase(object): 

10 """ 

11 Generic collocation class, that contains everything to do integration over 

12 intervals and between nodes. 

13 It can be used to produce many kind of quadrature nodes from various 

14 distribution (awesome!). 

15 

16 It is based on the two main parameters that define the nodes: 

17 

18 - node_type: the node distribution used for the collocation method 

19 - quad_type: the type of quadrature used (inclusion of not of boundary) 

20 

21 Current implementation provides the following available parameter values 

22 for node_type: 

23 

24 - EQUID: equidistant node distribution 

25 - LEGENDRE: distribution from Legendre polynomials 

26 - CHEBY-{1,2,3,4}: distribution from Chebyshev polynomials of a given kind 

27 

28 The type of quadrature can be GAUSS (only inner nodes), RADAU-LEFT 

29 (inclusion of the left boundary), RADAU-RIGHT (inclusion of the right 

30 boundary) and LOBATTO (inclusion of left and right boundary). 

31 

32 All coefficients are generated using 

33 `qmat <https://qmat.readthedocs.io/en/latest/api/qmat/qcoeff/collocation/index.html>`_. 

34 

35 Attributes: 

36 num_nodes (int): number of collocation nodes 

37 tleft (float): left interval point 

38 tright (float): right interval point 

39 nodes (numpy.ndarray): array of quadrature nodes 

40 weights (numpy.ndarray): array of quadrature weights for the full interval 

41 Qmat (numpy.ndarray): matrix containing the weights for tleft to node 

42 Smat (numpy.ndarray): matrix containing the weights for node to node 

43 delta_m (numpy.ndarray): array of distances between nodes 

44 right_is_node (bool): flag to indicate whether right point is collocation node 

45 left_is_node (bool): flag to indicate whether left point is collocation node 

46 """ 

47 

48 def __init__( 

49 self, 

50 num_nodes: Optional[int] = None, 

51 tleft: float = 0, 

52 tright: float = 1, 

53 node_type: str = 'LEGENDRE', 

54 quad_type: Optional[str] = None, 

55 **kwargs: Any, 

56 ) -> None: 

57 """ 

58 Initialization routine for a collocation object 

59 

60 Args: 

61 num_nodes (int): number of collocation nodes 

62 tleft (float): left interval point 

63 tright (float): right interval point 

64 node_type (str): type of the node distribution, see above 

65 quad_type (str): type of quadrature, see above; there is no default, it has to be given 

66 """ 

67 

68 if not num_nodes > 0: 

69 raise CollocationError('at least one quadrature node required, got %s' % num_nodes) 

70 if not tleft < tright: 

71 raise CollocationError('interval boundaries are corrupt, got %s and %s' % (tleft, tright)) 

72 

73 self.logger: logging.Logger = logging.getLogger('collocation') 

74 try: 

75 self.generator = Q_GENERATORS["Collocation"]( 

76 nNodes=num_nodes, nodeType=node_type, quadType=quad_type, tLeft=tleft, tRight=tright 

77 ) 

78 except Exception as e: 

79 raise CollocationError(f"could not instantiate qmat generator, got error: {e}") from e 

80 

81 # Set base attributes 

82 self.num_nodes = num_nodes 

83 self.tleft = tleft 

84 self.tright = tright 

85 self.node_type = node_type 

86 self.quad_type = quad_type 

87 self.left_is_node = self.quad_type in ['LOBATTO', 'RADAU-LEFT'] 

88 self.right_is_node = self.quad_type in ['LOBATTO', 'RADAU-RIGHT'] 

89 

90 # Integration order 

91 self.order = self.generator.order 

92 

93 # Compute coefficients 

94 self.nodes = self._getNodes = self.generator.nodes.copy() 

95 self.weights = self.generator.weights.copy() 

96 

97 Q = np.zeros([num_nodes + 1, num_nodes + 1], dtype=float) 

98 Q[1:, 1:] = self.generator.Q 

99 self.Qmat = Q 

100 

101 S = np.zeros([num_nodes + 1, num_nodes + 1], dtype=float) 

102 S[1:, 1:] = super(self.generator.__class__, self.generator).S 

103 # Note: qmat redefines the S matrix for collocation with integrals, 

104 # instead of differences of the Q matrix coefficients. 

105 # This does not passes the pySDC tests ... however the default S computation 

106 # in qmat uses Q matrix coefficients differences, and that's what we 

107 # use by using the parent property from the generator object. 

108 self.Smat = self._gen_Smatrix = S 

109 

110 self.delta_m = self._gen_deltas 

111 

112 @staticmethod 

113 def evaluate(weights: np.ndarray, data: np.ndarray) -> np.ndarray: 

114 """ 

115 Evaluates the quadrature over the full interval 

116 

117 Args: 

118 weights (numpy.ndarray): array of quadrature weights for the full interval 

119 data (numpy.ndarray): f(x) to be integrated 

120 

121 Returns: 

122 numpy.ndarray: integral over f(x) between tleft and tright 

123 """ 

124 if not np.size(weights) == np.size(data): 

125 raise CollocationError("Input size does not match number of weights, but is %s" % np.size(data)) 

126 

127 return np.dot(weights, data) 

128 

129 @property 

130 def _gen_deltas(self) -> np.ndarray: 

131 """ 

132 Compute distances between the nodes 

133 

134 Returns: 

135 numpy.ndarray: distances between the nodes 

136 """ 

137 M = self.num_nodes 

138 delta = np.zeros(M) 

139 delta[0] = self.nodes[0] - self.tleft 

140 for m in np.arange(1, M): 

141 delta[m] = self.nodes[m] - self.nodes[m - 1] 

142 

143 return delta