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

1import logging 

2import numpy as np 

3from qmat import Q_GENERATORS 

4 

5from pySDC.core.errors import CollocationError 

6 

7 

8class CollBase(object): 

9 """ 

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

11 intervals and between nodes. 

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

13 distribution (awesome!). 

14 

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

16 

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

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

19 

20 Current implementation provides the following available parameter values 

21 for node_type: 

22 

23 - EQUID: equidistant node distribution 

24 - LEGENDRE: distribution from Legendre polynomials 

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

26 

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

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

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

30 

31 All coefficients are generated using 

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

33 

34 Attributes: 

35 num_nodes (int): number of collocation nodes 

36 tleft (float): left interval point 

37 tright (float): right interval point 

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

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

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

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

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

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

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

45 """ 

46 

47 def __init__(self, num_nodes=None, tleft=0, tright=1, node_type='LEGENDRE', quad_type=None, **kwargs): 

48 """ 

49 Initialization routine for a collocation object 

50 

51 Args: 

52 num_nodes (int): number of collocation nodes 

53 tleft (float): left interval point 

54 tright (float): right interval point 

55 """ 

56 

57 if not num_nodes > 0: 

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

59 if not tleft < tright: 

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

61 

62 self.logger = logging.getLogger('collocation') 

63 try: 

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

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

66 ) 

67 except Exception as e: 

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

69 

70 # Set base attributes 

71 self.num_nodes = num_nodes 

72 self.tleft = tleft 

73 self.tright = tright 

74 self.node_type = node_type 

75 self.quad_type = quad_type 

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

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

78 

79 # Integration order 

80 self.order = self.generator.order 

81 

82 # Compute coefficients 

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

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

85 

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

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

88 self.Qmat = Q 

89 

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

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

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

93 # instead of differences of the Q matrix coefficients. 

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

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

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

97 self.Smat = self._gen_Smatrix = S 

98 

99 self.delta_m = self._gen_deltas 

100 

101 @staticmethod 

102 def evaluate(weights, data): 

103 """ 

104 Evaluates the quadrature over the full interval 

105 

106 Args: 

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

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

109 

110 Returns: 

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

112 """ 

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

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

115 

116 return np.dot(weights, data) 

117 

118 @property 

119 def _gen_deltas(self): 

120 """ 

121 Compute distances between the nodes 

122 

123 Returns: 

124 numpy.ndarray: distances between the nodes 

125 """ 

126 M = self.num_nodes 

127 delta = np.zeros(M) 

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

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

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

131 

132 return delta