|
| 1 | +import pySDC.Collocation |
| 2 | +from pySDC.CollocationClasses import * |
| 3 | +import numpy as np |
| 4 | +import pytest |
| 5 | + |
| 6 | +# py.test excludes classes with a constructor by default, so define parameter here |
| 7 | + |
| 8 | +# generate random boundaries for the time slice with 0.0 <= t_start < 0.2 and 0.8 <= t_end < 1.0 |
| 9 | +t_start = np.random.rand(1)*0.2 |
| 10 | +t_end = 0.8 + np.random.rand(1)*0.2 |
| 11 | + |
| 12 | +classes = [ ["CollGaussLegendre", 2, 12], ["CollGaussLobatto", 2, 12], ["CollGaussRadau_Right", 2, 12] ] |
| 13 | + |
| 14 | +class TestCollocation: |
| 15 | + |
| 16 | + # TEST 1: |
| 17 | + # Check that the quadrature rule integrates polynomials up to order p-1 exactly |
| 18 | + # ----------------------------------------------------------------------------- |
| 19 | + def test_1(self): |
| 20 | + for type in classes: |
| 21 | + for M in range(type[1],type[2]+1): |
| 22 | + coll = getattr(pySDC.CollocationClasses, type[0])(M, t_start, t_end) |
| 23 | + |
| 24 | + # some basic consistency tests |
| 25 | + assert np.size(coll.nodes)==np.size(coll.weights), "For node type " + type[0] + ", number of entries in nodes and weights is different" |
| 26 | + assert np.size(coll.nodes)==M, "For node type " + type[0] + ", requesting M nodes did not produce M entries in nodes and weights" |
| 27 | + |
| 28 | + |
| 29 | + # generate random set of polynomial coefficients |
| 30 | + poly_coeff = np.random.rand(coll.order-1) |
| 31 | + # evaluate polynomial at collocation nodes |
| 32 | + poly_vals = np.polyval(poly_coeff, coll.nodes) |
| 33 | + # use python's polyint function to compute anti-derivative of polynomial |
| 34 | + poly_int_coeff = np.polyint(poly_coeff) |
| 35 | + # Compute integral from 0.0 to 1.0 |
| 36 | + int_ex = np.polyval(poly_int_coeff, t_end) - np.polyval(poly_int_coeff, t_start) |
| 37 | + # use quadrature rule to compute integral |
| 38 | + int_coll = coll.evaluate(coll.weights, poly_vals) |
| 39 | + # For large values of M, substantial differences from different round of error have to be considered |
| 40 | + assert abs(int_ex - int_coll) < 1e-10, "For node type " + type[0] + ", failed to integrate polynomial of degree " + str(coll.order-1) + " exactly. Error: %5.3e" % abs(int_ex - int_coll) |
| 41 | + |
| 42 | + |
| 43 | + # TEST 2: |
| 44 | + # Check that the Qmat entries are equal to the sum of Smat entries |
| 45 | + # ---------------------------------------------------------------- |
| 46 | + def test_2(self): |
| 47 | + for type in classes: |
| 48 | + for M in range(type[1],type[2]+1): |
| 49 | + coll = getattr(pySDC.CollocationClasses, type[0])(M, t_start, t_end) |
| 50 | + Q = coll.Qmat[1:,1:] |
| 51 | + S = coll.Smat[1:,1:] |
| 52 | + assert np.shape(Q) == np.shape(S), "For node type " + type[0] + ", Qmat and Smat have different shape" |
| 53 | + shape = np.shape(Q) |
| 54 | + assert shape[0] == shape[1], "For node type " + type[0] + ", Qmat / Smat are not quadratic" |
| 55 | + SSum = np.cumsum(S[:,:],axis=0) |
| 56 | + for i in range(0,M): |
| 57 | + assert np.linalg.norm( Q[i,:] - SSum[i,:] ) < 1e-15, "For node type " + type[0] + ", Qmat and Smat did not satisfy the expected summation property." |
| 58 | + |
| 59 | + # TEST 3: |
| 60 | + # Check that the partial quadrature rules from Qmat entries have order equal to number of nodes M |
| 61 | + # ----------------------------------------------------------------------------------------------- |
| 62 | + def test_3(self): |
| 63 | + for type in classes: |
| 64 | + for M in range(type[1],type[2]+1): |
| 65 | + coll = getattr(pySDC.CollocationClasses, type[0])(M, t_start, t_end) |
| 66 | + Q = coll.Qmat[1:,1:] |
| 67 | + # as in TEST 1, create and integrate a polynomial with random coefficients, but now of degree M-1 |
| 68 | + poly_coeff = np.random.rand(M-1) |
| 69 | + poly_vals = np.polyval(poly_coeff, coll.nodes) |
| 70 | + poly_int_coeff = np.polyint(poly_coeff) |
| 71 | + for i in range(0,M): |
| 72 | + int_ex = np.polyval(poly_int_coeff, coll.nodes[i]) - np.polyval(poly_int_coeff, t_start) |
| 73 | + int_coll = np.dot(poly_vals, Q[i,:]) |
| 74 | + assert abs(int_ex - int_coll)<1e-12, "For node type " + type[0] + ", partial quadrature from Qmat rule failed to integrate polynomial of degree M-1 exactly for M = " + str(M) |
| 75 | + |
| 76 | + def test_4(self): |
| 77 | + for type in classes: |
| 78 | + for M in range(type[1],type[2]+1): |
| 79 | + coll = getattr(pySDC.CollocationClasses, type[0])(M, t_start, t_end) |
| 80 | + S = coll.Smat[1:,1:] |
| 81 | + # as in TEST 1, create and integrate a polynomial with random coefficients, but now of degree M-1 |
| 82 | + poly_coeff = np.random.rand(M-1) |
| 83 | + poly_vals = np.polyval(poly_coeff, coll.nodes) |
| 84 | + poly_int_coeff = np.polyint(poly_coeff) |
| 85 | + for i in range(1,M): |
| 86 | + int_ex = np.polyval(poly_int_coeff, coll.nodes[i]) - np.polyval(poly_int_coeff, coll.nodes[i-1]) |
| 87 | + int_coll = np.dot(poly_vals, S[i,:]) |
| 88 | + assert abs(int_ex - int_coll)<1e-12, "For node type " + type[0] + ", partial quadrature rule from Smat failed to integrate polynomial of degree M-1 exactly for M = " + str(M) |
0 commit comments