From 31eed1d853e5476a023082b932949bc2575db851 Mon Sep 17 00:00:00 2001
From: Daniel Ruprecht <d.ruprecht@leeds.ac.uk>
Date: Mon, 25 Apr 2016 15:41:35 +0100
Subject: [PATCH] Add special integrator which is defined by passing a
 stability function to the constructor. Allows to create tailored integrators
 for testing purposes

---
 src/special_integrator.py        | 17 +++++++++
 tests/test_special_integrator.py | 63 ++++++++++++++++++++++++++++++++
 2 files changed, 80 insertions(+)
 create mode 100644 src/special_integrator.py
 create mode 100644 tests/test_special_integrator.py

diff --git a/src/special_integrator.py b/src/special_integrator.py
new file mode 100644
index 0000000..f22243f
--- /dev/null
+++ b/src/special_integrator.py
@@ -0,0 +1,17 @@
+from integrator import integrator
+from solution import solution
+from solution_linear import solution_linear
+
+class special_integrator(integrator):
+
+  def __init__(self, tstart, tend, nsteps, stab_function):
+    super(special_integrator, self).__init__(tstart, tend, nsteps)
+    self.stab_function = stab_function
+
+  def run(self, u0):
+    assert isinstance(u0, solution_linear), "special_integrator can only be run for solutions of type solution_linear"
+    for i in range(0,self.nsteps):
+      u0.y = self.stab_function.dot(u0.y)
+
+  def get_update_matrix(self):
+    return self.stab_function
diff --git a/tests/test_special_integrator.py b/tests/test_special_integrator.py
new file mode 100644
index 0000000..9494eb0
--- /dev/null
+++ b/tests/test_special_integrator.py
@@ -0,0 +1,63 @@
+import sys
+sys.path.append('../src')
+
+import numpy as np
+from scipy import sparse
+from scipy import linalg
+from integrator import integrator
+from impeuler import impeuler
+from trapezoidal import trapezoidal
+from special_integrator import special_integrator
+from solution_linear import solution_linear
+from nose.tools import *
+import unittest
+import copy
+
+class TestSpecialIntegrator(unittest.TestCase):
+
+  def setUp(self):
+    self.ndof = np.random.randint(255)
+    self.A    = sparse.spdiags([ np.ones(self.ndof), -2.0*np.ones(self.ndof), np.ones(self.ndof)], [-1,0,1], self.ndof, self.ndof, format="csc")
+    self.M    = sparse.spdiags([ np.random.rand(self.ndof) ], [0], self.ndof, self.ndof, format="csc")
+    self.sol  = solution_linear(np.ones(self.ndof), self.A, self.M)
+
+  # Can instantiate object
+  def test_caninstantiate(self):
+    integ = special_integrator(0.0, 1.0, 10, sparse.eye(self.ndof))
+
+  # Throws if tend < tstart
+  def test_tendbeforetstartthrows(self):
+    with self.assertRaises(AssertionError):
+      integ = special_integrator(1.0, 0.0, 10, sparse.eye(self.ndof))
+
+  # See if run function can be called
+  def test_cancallrun(self):
+    integ = special_integrator(0.0, 1.0, 10, sparse.eye(self.ndof))
+    integ.run(self.sol)
+
+  # Make sure it provides identical solution to other integrator when stability function is given
+  def test_reproducesimpeuler(self):
+    ie = impeuler(0.0, 1.0, 25)
+    M = sparse.csc_matrix(self.sol.M)
+    A = sparse.csc_matrix(self.sol.A)
+    Rmat = sparse.linalg.inv(M - ie.dt*A)
+    Rmat = Rmat.dot(M)
+    integ = special_integrator(0.0, 1.0, 25, Rmat)
+    sol2 = copy.deepcopy(self.sol)
+    ie.run(sol2)
+    integ.run(self.sol)
+    self.sol.axpy(-1.0, sol2)
+    assert (self.sol.norm()<1e-14), "special_integrator did produce output identical to backward Euler"
+
+  def test_reproducetrapezoidal(self):
+    trap = trapezoidal(0.0, 1.0, 25)
+    M = sparse.csc_matrix(self.sol.M)
+    A = sparse.csc_matrix(self.sol.A)
+    Rmat = sparse.linalg.inv(M - 0.5*trap.dt*A)
+    Rmat = Rmat.dot(M+0.5*trap.dt*A)
+    integ = special_integrator(0.0, 1.0, 25, Rmat)
+    sol2 = copy.deepcopy(self.sol)
+    trap.run(sol2)
+    integ.run(self.sol)
+    self.sol.axpy(-1.0, sol2)
+    assert (self.sol.norm()<1e-14), "special_integrator did produce output identical to trapezoidal rule"