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

65 

66 if not num_nodes > 0: 

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

68 if not tleft < tright: 

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

70 

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

72 try: 

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

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

75 ) 

76 except Exception as e: 

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

78 

79 # Set base attributes 

80 self.num_nodes = num_nodes 

81 self.tleft = tleft 

82 self.tright = tright 

83 self.node_type = node_type 

84 self.quad_type = quad_type 

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

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

87 

88 # Integration order 

89 self.order = self.generator.order 

90 

91 # Compute coefficients 

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

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

94 

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

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

97 self.Qmat = Q 

98 

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

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

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

102 # instead of differences of the Q matrix coefficients. 

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

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

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

106 self.Smat = self._gen_Smatrix = S 

107 

108 self.delta_m = self._gen_deltas 

109 

110 @staticmethod 

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

112 """ 

113 Evaluates the quadrature over the full interval 

114 

115 Args: 

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

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

118 

119 Returns: 

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

121 """ 

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

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

124 

125 return np.dot(weights, data) 

126 

127 @property 

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

129 """ 

130 Compute distances between the nodes 

131 

132 Returns: 

133 numpy.ndarray: distances between the nodes 

134 """ 

135 M = self.num_nodes 

136 delta = np.zeros(M) 

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

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

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

140 

141 return delta