diff --git a/examples/advection_1d_implicit/testFDOrder.py b/examples/advection_1d_implicit/plotFDconvergence.py similarity index 100% rename from examples/advection_1d_implicit/testFDOrder.py rename to examples/advection_1d_implicit/plotFDconvergence.py diff --git a/examples/boussinesq_2d_imex/playground.py b/examples/boussinesq_2d_imex/playground.py index 69401f0f2c..f75816e03b 100644 --- a/examples/boussinesq_2d_imex/playground.py +++ b/examples/boussinesq_2d_imex/playground.py @@ -28,7 +28,7 @@ # This comes as read-in for the level class lparams = {} - lparams['restol'] = 1E-14 + lparams['restol'] = 1E-8 swparams = {} swparams['collocation_class'] = collclass.CollGaussLobatto @@ -36,14 +36,14 @@ swparams['do_LU'] = True sparams = {} - sparams['maxiter'] = 4 + sparams['maxiter'] = 12 # setup parameters "in time" t0 = 0 - #Tend = 3000 - #Nsteps = 500 - Tend = 30 - Nsteps = 5 + Tend = 3000 + Nsteps = 500 + #Tend = 30 + #Nsteps = 5 dt = Tend/float(Nsteps) # This comes as read-in for the problem class @@ -58,7 +58,7 @@ pparams['order_upw'] = [5, 1] pparams['gmres_maxiter'] = [50, 50] pparams['gmres_restart'] = [20, 20] - pparams['gmres_tol'] = [1e-6, 1e-6] + pparams['gmres_tol'] = [1e-8, 1e-8] # This comes as read-in for the transfer operations tparams = {} diff --git a/tests/test_collocation.py b/tests/test_collocation.py new file mode 100644 index 0000000000..efa7cfd687 --- /dev/null +++ b/tests/test_collocation.py @@ -0,0 +1,88 @@ +import pySDC.Collocation +from pySDC.CollocationClasses import * +import numpy as np +import pytest + +# py.test excludes classes with a constructor by default, so define parameter here + +# generate random boundaries for the time slice with 0.0 <= t_start < 0.2 and 0.8 <= t_end < 1.0 +t_start = np.random.rand(1)*0.2 +t_end = 0.8 + np.random.rand(1)*0.2 + +classes = [ ["CollGaussLegendre", 2, 12], ["CollGaussLobatto", 2, 12], ["CollGaussRadau_Right", 2, 12] ] + +class TestCollocation: + + # TEST 1: + # Check that the quadrature rule integrates polynomials up to order p-1 exactly + # ----------------------------------------------------------------------------- + def test_1(self): + for type in classes: + for M in range(type[1],type[2]+1): + coll = getattr(pySDC.CollocationClasses, type[0])(M, t_start, t_end) + + # some basic consistency tests + assert np.size(coll.nodes)==np.size(coll.weights), "For node type " + type[0] + ", number of entries in nodes and weights is different" + assert np.size(coll.nodes)==M, "For node type " + type[0] + ", requesting M nodes did not produce M entries in nodes and weights" + + + # generate random set of polynomial coefficients + poly_coeff = np.random.rand(coll.order-1) + # evaluate polynomial at collocation nodes + poly_vals = np.polyval(poly_coeff, coll.nodes) + # use python's polyint function to compute anti-derivative of polynomial + poly_int_coeff = np.polyint(poly_coeff) + # Compute integral from 0.0 to 1.0 + int_ex = np.polyval(poly_int_coeff, t_end) - np.polyval(poly_int_coeff, t_start) + # use quadrature rule to compute integral + int_coll = coll.evaluate(coll.weights, poly_vals) + # For large values of M, substantial differences from different round of error have to be considered + 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) + + + # TEST 2: + # Check that the Qmat entries are equal to the sum of Smat entries + # ---------------------------------------------------------------- + def test_2(self): + for type in classes: + for M in range(type[1],type[2]+1): + coll = getattr(pySDC.CollocationClasses, type[0])(M, t_start, t_end) + Q = coll.Qmat[1:,1:] + S = coll.Smat[1:,1:] + assert np.shape(Q) == np.shape(S), "For node type " + type[0] + ", Qmat and Smat have different shape" + shape = np.shape(Q) + assert shape[0] == shape[1], "For node type " + type[0] + ", Qmat / Smat are not quadratic" + SSum = np.cumsum(S[:,:],axis=0) + for i in range(0,M): + 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." + + # TEST 3: + # Check that the partial quadrature rules from Qmat entries have order equal to number of nodes M + # ----------------------------------------------------------------------------------------------- + def test_3(self): + for type in classes: + for M in range(type[1],type[2]+1): + coll = getattr(pySDC.CollocationClasses, type[0])(M, t_start, t_end) + Q = coll.Qmat[1:,1:] + # as in TEST 1, create and integrate a polynomial with random coefficients, but now of degree M-1 + poly_coeff = np.random.rand(M-1) + poly_vals = np.polyval(poly_coeff, coll.nodes) + poly_int_coeff = np.polyint(poly_coeff) + for i in range(0,M): + int_ex = np.polyval(poly_int_coeff, coll.nodes[i]) - np.polyval(poly_int_coeff, t_start) + int_coll = np.dot(poly_vals, Q[i,:]) + 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) + + def test_4(self): + for type in classes: + for M in range(type[1],type[2]+1): + coll = getattr(pySDC.CollocationClasses, type[0])(M, t_start, t_end) + S = coll.Smat[1:,1:] + # as in TEST 1, create and integrate a polynomial with random coefficients, but now of degree M-1 + poly_coeff = np.random.rand(M-1) + poly_vals = np.polyval(poly_coeff, coll.nodes) + poly_int_coeff = np.polyint(poly_coeff) + for i in range(1,M): + int_ex = np.polyval(poly_int_coeff, coll.nodes[i]) - np.polyval(poly_int_coeff, coll.nodes[i-1]) + int_coll = np.dot(poly_vals, S[i,:]) + 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) diff --git a/tests/tests_core.py b/tests/tests_core.py index e2d5bae227..a680e2b337 100644 --- a/tests/tests_core.py +++ b/tests/tests_core.py @@ -1,9 +1,9 @@ import numpy as np -def test_python_version(): - import sys - assert sys.version_info >= (3, 3), "Need Python 3.3 at least" +#def test_python_version(): +# import sys +# assert sys.version_info >= (3, 3), "Need Python 3.3 at least" # def test_return_types(): # import pySDC.problem as pr @@ -177,4 +177,4 @@ def check_datatypes_particles(init): assert p7 >= 0 assert np.all(p8.pos.values==1.0) assert np.all(p8.vel.values==10.0) - assert np.all(a3.values==300.0) \ No newline at end of file + assert np.all(a3.values==300.0)