Coverage for pySDC/core/Collocation.py: 100%
84 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-04-29 09:02 +0000
1import logging
2import numpy as np
4from pySDC.core.Nodes import NodesGenerator
5from pySDC.core.Errors import CollocationError
6from pySDC.core.Lagrange import LagrangeApproximation
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!).
16 It is based on the two main parameters that define the nodes :
18 - node_type : the node distribution used for the collocation method
19 - quad_type : the type of quadrature used (inclusion of not of boundary)
21 Current implementation provides the following available parameter values
22 for node_type :
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
28 The type of quadrature cann 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).
32 Here is the equivalency table with the (old) original classes implemented
33 in pySDC :
35 +-------------------------+-----------+-------------+
36 | Original Class | node_type | quad_type |
37 +=========================+===========+=============+
38 | Equidistant | EQUID | LOBATTO |
39 +-------------------------+-----------+-------------+
40 | EquidistantInner | EQUID | GAUSS |
41 +-------------------------+-----------+-------------+
42 | EquidistantNoLeft | EQUID | RADAU-RIGHT |
43 +-------------------------+-----------+-------------+
44 | CollGaussLegendre | LEGENDRE | GAUSS |
45 +-------------------------+-----------+-------------+
46 | CollGaussLobatto | LEGENDRE | LOBATTO |
47 +-------------------------+-----------+-------------+
48 | CollGaussRadau_Left | LEGENDRE | RADAU-LEFT |
49 +-------------------------+-----------+-------------+
50 | CollGaussRadau_Right | LEGENDRE | RADAU-RIGHT |
51 +-------------------------+-----------+-------------+
53 Attributes:
54 num_nodes (int): number of collocation nodes
55 tleft (float): left interval point
56 tright (float): right interval point
57 nodes (numpy.ndarray): array of quadrature nodes
58 weights (numpy.ndarray): array of quadrature weights for the full interval
59 Qmat (numpy.ndarray): matrix containing the weights for tleft to node
60 Smat (numpy.ndarray): matrix containing the weights for node to node
61 delta_m (numpy.ndarray): array of distances between nodes
62 right_is_node (bool): flag to indicate whether right point is collocation node
63 left_is_node (bool): flag to indicate whether left point is collocation node
64 """
66 def __init__(self, num_nodes=None, tleft=0, tright=1, node_type='LEGENDRE', quad_type=None, **kwargs):
67 """
68 Initialization routine for a collocation object
70 Args:
71 num_nodes (int): number of collocation nodes
72 tleft (float): left interval point
73 tright (float): right interval point
74 """
76 if not num_nodes > 0:
77 raise CollocationError('At least one quadrature node required, got %s' % num_nodes)
78 if not tleft < tright:
79 raise CollocationError('Interval boundaries are corrupt, got %s and %s' % (tleft, tright))
81 self.logger = logging.getLogger('collocation')
83 # Set number of nodes, left and right interval boundaries
84 self.num_nodes = num_nodes
85 self.tleft = tleft
86 self.tright = tright
88 self.node_type = node_type
89 self.quad_type = quad_type
91 # Instantiate attributes
92 self.nodeGenerator = NodesGenerator(self.node_type, self.quad_type)
93 if self.node_type == 'EQUID':
94 self.order = num_nodes
95 else:
96 if self.quad_type == 'GAUSS':
97 self.order = 2 * num_nodes
98 elif self.quad_type.startswith('RADAU'):
99 self.order = 2 * num_nodes - 1
100 elif self.quad_type == 'LOBATTO':
101 self.order = 2 * num_nodes - 2
103 self.left_is_node = self.quad_type in ['LOBATTO', 'RADAU-LEFT']
104 self.right_is_node = self.quad_type in ['LOBATTO', 'RADAU-RIGHT']
106 self.nodes = self._getNodes
107 self.weights = self._getWeights(tleft, tright)
108 self.Qmat = self._gen_Qmatrix
109 self.Smat = self._gen_Smatrix
110 self.delta_m = self._gen_deltas
112 @staticmethod
113 def evaluate(weights, data):
114 """
115 Evaluates the quadrature over the full interval
117 Args:
118 weights (numpy.ndarray): array of quadrature weights for the full interval
119 data (numpy.ndarray): f(x) to be integrated
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))
127 return np.dot(weights, data)
129 def _getWeights(self, a, b):
130 """
131 Computes weights using barycentric interpolation
133 Args:
134 a (float): left interval boundary
135 b (float): right interval boundary
137 Returns:
138 numpy.ndarray: weights of the collocation formula given by the nodes
139 """
140 if self.nodes is None:
141 raise CollocationError(f"Need nodes before computing weights, got {self.nodes}")
143 # Instantiate the Lagrange interpolator object
144 approx = LagrangeApproximation(self.nodes)
146 # Compute weights
147 tLeft = np.ravel(self.tleft)[0]
148 tRight = np.ravel(self.tright)[0]
149 weights = approx.getIntegrationMatrix([(tLeft, tRight)], numQuad='FEJER')
151 return np.ravel(weights)
153 @property
154 def _getNodes(self):
155 """
156 Computes nodes using an internal NodesGenerator object
158 Returns:
159 np.ndarray: array of Gauss-Legendre nodes
160 """
161 # Generate normalized nodes in [-1, 1]
162 nodes = self.nodeGenerator.getNodes(self.num_nodes)
164 # Scale nodes to [tleft, tright]
165 a = self.tleft
166 b = self.tright
167 nodes += 1.0
168 nodes /= 2.0
169 nodes *= b - a
170 nodes += a
172 if self.left_is_node:
173 nodes[0] = self.tleft
174 if self.right_is_node:
175 nodes[-1] = self.tright
177 return nodes
179 @property
180 def _gen_Qmatrix(self):
181 """
182 Compute tleft-to-node integration matrix for later use in collocation formulation
184 Returns:
185 numpy.ndarray: matrix containing the weights for tleft to node
186 """
187 if self.nodes is None:
188 raise CollocationError(f"Need nodes before computing weights, got {self.nodes}")
189 M = self.num_nodes
190 Q = np.zeros([M + 1, M + 1])
192 # Instantiate the Lagrange interpolator object
193 approx = LagrangeApproximation(self.nodes)
195 # Compute tleft-to-node integration matrix
196 tLeft = np.ravel(self.tleft)[0]
197 intervals = [(tLeft, tau) for tau in self.nodes]
198 intQ = approx.getIntegrationMatrix(intervals, numQuad='FEJER')
200 # Store into Q matrix
201 Q[1:, 1:] = intQ
203 return Q
205 @property
206 def _gen_Smatrix(self):
207 """
208 Compute node-to-node integration matrix for later use in collocation formulation
210 Returns:
211 numpy.ndarray: matrix containing the weights for node to node
212 """
213 M = self.num_nodes
214 Q = self.Qmat
215 S = np.zeros([M + 1, M + 1])
217 S[1, :] = Q[1, :]
218 for m in np.arange(2, M + 1):
219 S[m, :] = Q[m, :] - Q[m - 1, :]
221 return S
223 @property
224 def _gen_deltas(self):
225 """
226 Compute distances between the nodes
228 Returns:
229 numpy.ndarray: distances between the nodes
230 """
231 M = self.num_nodes
232 delta = np.zeros(M)
233 delta[0] = self.nodes[0] - self.tleft
234 for m in np.arange(1, M):
235 delta[m] = self.nodes[m] - self.nodes[m - 1]
237 return delta