Implemented GMRES iterative method

hdpg1d/adaptation.py

@@ -1,5 +1,7 @@
 import numpy as np
 from numpy import concatenate as cat
+from scipy.sparse import csr_matrix
+import scipy.sparse.linalg as spla
 from copy import copy
 import matplotlib.pyplot as plt
 import warnings
@@ -83,14 +85,22 @@ class hdpg1d(object):
         K = -cat((C.T, G), axis=1)\
             .dot(np.linalg.inv(np.bmat([[A, -B], [B.T, D]]))
                  .dot(cat((C, E)))) + H
+        sK = csr_matrix(K)
         F_hat = np.array([L]).T - cat((C.T, G), axis=1)\
             .dot(np.linalg.inv(np.bmat([[A, -B], [B.T, D]])))\
             .dot(np.array([cat((R, F))]).T)
-        stateFace = np.linalg.solve(K, F_hat)
-        gradState = np.linalg.inv(np.bmat([[A, -B], [B.T, D]]))\
-            .dot(np.array([np.concatenate((R, F))]).T -
-                 cat((C, E)).dot(stateFace))
-        self.primalSoln = cat((gradState.A1, stateFace.A1))
+
+        def M_x(vec):
+            """Construct preconditioner"""
+            matVec = spla.spsolve(sK, vec)
+            return matVec
+        n = len(F_hat)
+        preconditioner = spla.LinearOperator((n, n), M_x)
+        stateFace = spla.gmres(sK, F_hat, M=preconditioner)[0]
+        # stateFace = np.linalg.solve(K, F_hat)
+        gradState = np.linalg.inv(np.asarray(np.bmat([[A, -B], [B.T, D]]))).dot(
+            cat((R, F)) - cat((C, E)).dot(stateFace))
+        self.primalSoln = cat((gradState, stateFace))
 
     def solveAdjoint(self):
         """Solve the adjoint problem"""
@@ -110,8 +120,18 @@ class hdpg1d(object):
         LHS = np.bmat([[A, -B, C],
                        [B.T, D, E],
                        [C.T, G, H]])
-        # solve in one shoot
-        soln = np.linalg.solve(LHS.T, cat((R, F, L)))
+        sLHS = csr_matrix(LHS)
+        RHS = cat((R, F, L))
+
+        # solve in one shoot using GMRES
+
+        def M_x(vec):
+            """Construct preconditioner"""
+            matVec = spla.spsolve(sLHS, vec)
+            return matVec
+        n = len(RHS)
+        preconditioner = spla.LinearOperator((n, n), M_x)
+        soln = spla.gmres(sLHS, RHS, M=preconditioner)[0]
         self.adjointSoln = soln
 
     def DWResidual(self):
@@ -152,7 +172,7 @@ class hdpg1d(object):
         TOL = 1e-10
         estError = 10
         nodeCount = 0
-        maxCount = 30
+        maxCount = 40
         while estError > TOL and nodeCount < maxCount:
             # solve
             self.solvePrimal()

tests/test_solve_test_problems.py

@@ -7,18 +7,26 @@ from hdpg1d.coefficients import coefficients as coeff
 from hdpg1d.adaptation import hdpg1d
 
 
+testData = [
+    ([1e-4, 0, 1, 2, 2, 1, 1], 1e-2),  # diffusion reaction
+    ([0, 1, 0, 2, 2, 1, 1], 0),        # convection
+    # ([1, 1, 0, 2, 2, 1, 1], 1)         # diffusion convection)
+]
+
+
 class TestClass(object):
+    @pytest.fixture(scope="module", params=testData)
+    def coeffGen(self, request):
+        coeffTest = coeff(*request.param[0])
+        expected = request.param[1]
+        yield coeffTest, expected  # teardown
+
     def test_zeroDivision(self, monkeypatch):
         coeffTest = coeff(*([0] * 7))
         assert coeffTest.DIFFUSION != 0
 
-    @pytest.mark.parametrize("coeffInput, expected", [
-        ([1e-4, 0, 1, 2, 2, 1, 1], 1e-2),  # diffusion reaction
-        ([0, 1, 0, 2, 2, 1, 1], 0),        # convection
-        # ([1, 1, 0, 2, 2, 1, 1], 1)         # diffusion convection
-    ])
-    def test_solveAdaptive(self, coeffInput, expected):
-        coeffTest = coeff(*coeffInput)
+    def test_solveAdaptive(self, coeffGen):
+        coeffTest, expected = coeffGen
         hdpgTest = hdpg1d(coeffTest)
         hdpgTest.adaptive()
         # get the target function value
@@ -26,16 +34,11 @@ class TestClass(object):
         soln = hdpgTest.trueErrorList[1][-1]
         assert math.isclose(soln, expected, rel_tol=1e-5, abs_tol=1e-10)
 
-    @pytest.mark.parametrize("coeffInput, expected", [
-        ([1e-4, 0, 1, 2, 2, 1, 1], 1e-2),  # diffusion reaction
-        ([0, 1, 0, 2, 2, 1, 1], 0),        # convection
-        # ([1, 1, 0, 2, 2, 1, 1], 1)         # diffusion convection
-    ])
-    def test_solvePrimal(self, coeffInput, expected):
-        # test the primal solver with a refined mesh and poly order
-        coeffInput[3] = 5
-        coeffInput[4] = 200
-        coeffTest = coeff(*coeffInput)
+    def test_solvePrimal(self, coeffGen):
+        coeffTest, expected = coeffGen
+        # test the primal solver on a refined mesh and higher poly order
+        coeffTest.pOrder = 5
+        coeffTest.numEle = 300
         hdpgTest = hdpg1d(coeffTest)
         hdpgTest.solvePrimal()
         soln = hdpgTest.primalSoln[hdpgTest.numEle *