Feature/add tests

examples/boussinesq_2d_imex/playground.py

@@ -28,22 +28,22 @@
 
     # This comes as read-in for the level class
     lparams = {}
-    lparams['restol'] = 1E-14
+    lparams['restol'] = 1E-8
 
     swparams = {}
     swparams['collocation_class'] = collclass.CollGaussLobatto
     swparams['num_nodes'] = 3
     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 = {}

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)

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)
+    assert np.all(a3.values==300.0)