Merge branch 'release/v3.0'

This is a major update that includes multiple bug fixes, performance improvements, and optimization.

.travis.yml

@@ -0,0 +1,12 @@
+language: python
+python:
+  - "3.5"
+  - "3.6"
+
+# command to install dependencies
+install:
+  - pip install -r requirements.txt
+  - python setup.py install
+
+# command to run tests
+script: pytest

MANIFEST.in

@@ -0,0 +1,4 @@
+include hdpg1d/config/config.json
+include LICENSE
+include README.md
+include requirements.txt

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
@@ -69,7 +71,7 @@ class hdpg1d(object):
                           self.mesh[index[i] - 1]) / 2
         self.mesh = np.sort(np.insert(self.mesh, 0, inValue))
 
-    def solveLocal(self):
+    def solvePrimal(self):
         """Solve the primal problem"""
         if 'matLocal' in locals():
             # if matLocal exists,
@@ -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 invRHS(vec):
+            """Construct preconditioner"""
+            matVec = spla.spsolve(sK, vec)
+            return matVec
+        n = len(F_hat)
+        preconditioner = spla.LinearOperator((n, n), invRHS)
+        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"""
@@ -105,13 +115,23 @@ class hdpg1d(object):
         A, B, _, C, D, E, F, G, H, L, R = matGroup
         # add adjoint LHS conditions
         F = np.zeros(len(F))
-        R[-1] = -boundaryCondition(1)[1]
+        R[-1] = -boundaryCondition('adjoint')[1]
         # assemble global matrix LHS
         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 invRHS(vec):
+            """Construct preconditioner"""
+            matVec = spla.spsolve(sLHS, vec)
+            return matVec
+        n = len(RHS)
+        preconditioner = spla.LinearOperator((n, n), invRHS)
+        soln = spla.gmres(sLHS, RHS, M=preconditioner)[0]
+        # soln = np.linalg.solve(LHS.T, RHS)
         self.adjointSoln = soln
 
     def DWResidual(self):
@@ -149,13 +169,13 @@ class hdpg1d(object):
         return np.abs(np.sum(residual)), refineIndex
 
     def adaptive(self):
-        TOL = 1e-10
+        TOL = self.coeff.TOL
         estError = 10
         nodeCount = 0
-        maxCount = 30
+        maxCount = self.coeff.MAXIT
         while estError > TOL and nodeCount < maxCount:
             # solve
-            self.solveLocal()
+            self.solvePrimal()
             self.solveAdjoint()
             # plot the solution at certain counter
             if nodeCount in [0, 4, 9, 19, maxCount]:
@@ -173,12 +193,10 @@ class hdpg1d(object):
             self.meshAdapt(index)
             self.numEle = self.numEle + len(index)
             nodeCount += 1
-            print("Iteration {}. Target function error {:.3e}.".format(
-                nodeCount, estError))
+            print("Iteration {}. Estimated target function error {:.3e}."
+                  .format(nodeCount, estError))
             if nodeCount == maxCount:
                 print("Max iteration number is reached "
                       "while the convergence criterion is not satisfied.\n"
                       "Check the problem statement or "
                       "raise the max iteration number, then try again.\n")
-                from .solve import runInteractive
-                runInteractive()

hdpg1d/coefficients.py

@@ -1,33 +0,0 @@
-class coefficients:
-    def __init__(self, diff, conv, reaction, pOrder, numEle, tauPlus, tauMinus):
-        if diff == 0:
-            # set the diffusion constant to a small number
-            # to avoid division by zero error
-            diff = 1e-16
-        self.DIFFUSION = diff
-        self.CONVECTION = conv
-        self.REACTION = reaction
-        self.pOrder = pOrder
-        self.numEle = numEle
-        self.TAUPLUS = tauPlus
-        self.TAUMINUS = tauMinus
-
-    @classmethod
-    def fromInput(cls):
-        while True:
-            try:
-                print("Please provide the following coefficients.")
-                diff = float(input("Diffusion constant (float): "))
-                conv = float(input("Covection constant (float): "))
-                reaction = float(input("Reaction constant (float): "))
-                pOrder = int(input("Order of polynomials (int): "))
-                numEle = int(input("Number of elements (int): "))
-                tauPlus = float(input("Stablization parameter plus (float): "))
-                tauMinus = float(
-                    input("Stablization parameter minus (float): "))
-            except ValueError:
-                print("Sorry, wrong data type.")
-                continue
-            else:
-                break
-        return cls(diff, conv, reaction, pOrder, numEle, tauPlus, tauMinus)

hdpg1d/config/config.json

@@ -0,0 +1,18 @@
+{
+	"coefficients": {
+		"diffusion": 1e-6,
+		"convection": 0,
+		"reaction": 1,
+		"pOrder": 2,
+		"numEle": 2,
+		"tauPlus": 1,
+		"tauMinus": 1,
+        "tol": 1e-10,
+        "maxIt": 30
+	},
+	"forcing": "1",
+	"boundary": {
+		"left": 0,
+		"right": 0
+	}
+}

hdpg1d/postprocess.py

@@ -19,7 +19,7 @@ class utils(object):
         self.solution.coeff.pOrder = self.exactBasisFuncs - 1
         self.solution.mesh = np.linspace(0, 1, self.exactNumEle + 1)
         # approximate the exact solution for general problems
-        self.solution.solveLocal()
+        self.solution.solvePrimal()
         exactSoln = self.solution.separateSoln(self.solution.primalSoln)[0][
             self.exactNumEle * self.exactBasisFuncs - 1]
         # for the reaction diffusion test problem, we know the exact solution
@@ -33,7 +33,7 @@ class utils(object):
         numBasisFuncs = self.solution.coeff.pOrder + 1
         # solve on the uniform mesh
         self.solution.mesh = np.linspace(0, 1, numEle + 1)
-        self.solution.solveLocal()
+        self.solution.solvePrimal()
         gradState, _ = self.solution.separateSoln(self.solution.primalSoln)
         errorL2 = np.abs(
             gradState[numBasisFuncs * numEle - 1] - self.exactSoln)

hdpg1d/preprocess.py

@@ -1,8 +1,90 @@
+import os
+import json
+import ast
+import operator as op
 import numpy as np
 from collections import namedtuple
 from scipy.linalg import block_diag
 
 
+# load cfg file
+installDir = os.path.split(__file__)[0]
+cfgPath = os.path.join(installDir, "config")
+for loc in cfgPath, os.curdir, os.path.expanduser("~"):
+    try:
+        with open(os.path.join(loc, "config.json")) as source:
+            configdata = json.load(source)
+    except IOError:
+        pass
+
+# evaluate the input json function with only these math operators
+operators = {ast.Add: op.add, ast.Sub: op.sub, ast.Mult: op.mul,
+             ast.Div: op.truediv, ast.Pow: op.pow, ast.BitXor: op.xor,
+             ast.USub: op.neg}
+
+
+def eval_expr(expr):
+    return eval_(ast.parse(expr, mode='eval').body)
+
+
+def eval_(node):
+    if isinstance(node, ast.Num):  # <number>
+        return node.n
+    elif isinstance(node, "x"):
+        return node.n
+    elif isinstance(node, ast.BinOp):  # <left> <operator> <right>
+        return operators[type(node.op)](eval_(node.left), eval_(node.right))
+    elif isinstance(node, ast.UnaryOp):  # <operator> <operand> e.g., -1
+        return operators[type(node.op)](eval_(node.operand))
+    else:
+        raise TypeError(node)
+
+
+def queryYesNo(question, default="yes"):
+    valid = {"yes": True, "y": True, "ye": True,
+             "no": False, "n": False}
+    if default is None:
+        prompt = " [y/n] "
+    elif default == "yes":
+        prompt = " [Y/n] "
+    elif default == "no":
+        prompt = " [y/N] "
+    else:
+        raise ValueError("invalid default answer: '%s'" % default)
+
+    while True:
+        print(question + prompt, end='')
+        choice = input().lower()
+        if default is not None and choice == '':
+            return valid[default]
+        elif choice in valid:
+            return valid[choice]
+        else:
+            print("Please respond with 'yes' or 'no' "
+                  "(or 'y' or 'n').\n")
+
+
+def setDefaultCoefficients():
+    question = "Do you want to use the default parameters?"
+    isDefault = queryYesNo(question, "yes")
+    if (isDefault):
+        diffDefault = configdata["coefficients"]["diffusion"]
+        convDefault = configdata["coefficients"]["convection"]
+        reactionDefault = configdata["coefficients"]["reaction"]
+        pOrderDefault = configdata["coefficients"]["pOrder"]
+        numEleDefault = configdata["coefficients"]["numEle"]
+        tauPlusDefault = configdata["coefficients"]["tauPlus"]
+        tauMinusDefault = configdata["coefficients"]["tauMinus"]
+        tol = configdata["coefficients"]["tol"]
+        maxIt = configdata["coefficients"]["maxIt"]
+        coeff = coefficients(diffDefault, convDefault, reactionDefault,
+                             pOrderDefault, numEleDefault, tauPlusDefault,
+                             tauMinusDefault, tol, maxIt)
+    else:
+        coeff = coefficients.fromInput()
+    return coeff
+
+
 def shape(x, p):
     """generate p shape functions and its first order derivative
     at the given location x. x can be an array"""
@@ -15,20 +97,68 @@ def shape(x, p):
 
 
 def forcing(x):
-    f = 1
+    f = np.zeros(len(x))
+    for i, forcingItem in enumerate(x):
+        forcingExpr = configdata["forcing"]
+        # replace the 'x' in the json file with the function parameters
+        f[i] = eval_expr(forcingExpr.replace("x", str(forcingItem)))
     return f
 
 
 def boundaryCondition(case):
-    if case == 0:
+    if case == 'primal':
         # primal problem
-        bc = [0, 0]
-    if case == 1:
+        bcLeft = configdata["boundary"]["left"]
+        bcRight = configdata["boundary"]["right"]
+        bc = [bcLeft, bcRight]
+    if case == 'adjoint':
         # adjoint problem
         bc = [0, 1]
     return bc
 
 
+class coefficients:
+    def __init__(self, diff, conv, reaction, pOrder, numEle,
+                 tauPlus, tauMinus, tol, maxIt):
+        if diff == 0:
+            # set the diffusion constant to a small number
+            # to avoid division by zero error
+            diff = 1e-16
+        self.DIFFUSION = diff
+        self.CONVECTION = conv
+        self.REACTION = reaction
+        self.pOrder = pOrder
+        self.numEle = numEle
+        self.TAUPLUS = tauPlus
+        self.TAUMINUS = tauMinus
+        self.TOL = tol
+        self.MAXIT = maxIt
+
+    @classmethod
+    def fromInput(cls):
+        while True:
+            try:
+                print("Please provide the following coefficients.")
+                diff = float(input("Diffusion constant (float): "))
+                conv = float(input("Covection constant (float): "))
+                reaction = float(input("Reaction constant (float): "))
+                pOrder = int(input("Order of polynomials (int): "))
+                numEle = int(input("Number of elements (int): "))
+                tauPlus = float(input("Stablization parameter plus (float): "))
+                tauMinus = float(
+                    input("Stablization parameter minus (float): "))
+                tol = float(input("Error tolerance (float): "))
+                maxIt = int(input("Max adaptive iterations (int): "))
+            except ValueError:
+                print("Sorry, wrong data type.")
+                continue
+            else:
+                print("Something is wrong. Exit.")
+                break
+        return cls(diff, conv, reaction, pOrder, numEle,
+                   tauPlus, tauMinus, tol, maxIt)
+
+
 class discretization(object):
     """Given the problem statement and current mesh,
     construct the discretization matrices"""
@@ -89,14 +219,14 @@ class discretization(object):
                                             1 / 2 * (1 + self.xi) *
                                             self.dist[i - 1]))
             F[(i - 1) * numBasisFuncs:i * numBasisFuncs] = f
-        F[0] += (self.conv + self.tau_pos) * boundaryCondition(0)[0]
-        F[-1] += (-self.conv + self.tau_neg) * boundaryCondition(0)[1]
+        F[0] += (self.conv + self.tau_pos) * boundaryCondition('primal')[0]
+        F[-1] += (-self.conv + self.tau_neg) * boundaryCondition('primal')[1]
         # L, easy in 1d
         L = np.zeros(self.numEle - 1)
         # R, easy in 1d
         R = np.zeros(numBasisFuncs * self.numEle)
-        R[0] = boundaryCondition(0)[0]
-        R[-1] = -boundaryCondition(0)[1]
+        R[0] = boundaryCondition('primal')[0]
+        R[-1] = -boundaryCondition('primal')[1]
         lhsMat = namedtuple('lhsMat', ['F', 'L', 'R'])
         return lhsMat(F, L, R)
 

hdpg1d/solve.py

@@ -1,42 +1,8 @@
-from .coefficients import coefficients
+from .preprocess import setDefaultCoefficients
 from .adaptation import hdpg1d
 from .postprocess import utils
 
 
-def queryYesNo(question, default="yes"):
-    valid = {"yes": True, "y": True, "ye": True,
-             "no": False, "n": False}
-    if default is None:
-        prompt = " [y/n] "
-    elif default == "yes":
-        prompt = " [Y/n] "
-    elif default == "no":
-        prompt = " [y/N] "
-    else:
-        raise ValueError("invalid default answer: '%s'" % default)
-
-    while True:
-        print(question + prompt)
-        choice = input().lower()
-        if default is not None and choice == '':
-            return valid[default]
-        elif choice in valid:
-            return valid[choice]
-        else:
-            print("Please respond with 'yes' or 'no' "
-                  "(or 'y' or 'n').\n")
-
-
-def getCoefficients():
-    question = 'Do you want to use the default parameters?'
-    isDefault = queryYesNo(question, "yes")
-    if (isDefault):
-        Coeff = coefficients(1e-6, 0, 1, 2, 2, 1, 1)
-    else:
-        Coeff = coefficients.fromInput()
-    return Coeff
-
-
 def menu():
     menu = {}
     menu['1.'] = "Solve with HDG."
@@ -48,7 +14,7 @@ def menu():
 
 
 def hdgSolve():
-    hdgCoeff = getCoefficients()
+    hdgCoeff = setDefaultCoefficients()
     print("Solving...")
     hdgSolution = hdpg1d(hdgCoeff)
     # solve the problem adaptively and plot convergence history
@@ -58,9 +24,9 @@ def hdgSolve():
 
 
 def runInteractive():
-    menu()
-    selection = input("Please Select: ")
     while True:
+        menu()
+        selection = input("Please Select: ")
         if selection == '1':
             hdgSolve()
             break
@@ -70,5 +36,5 @@ def runInteractive():
             print("Bye.")
             break
         else:
-            print("Unknown Option Selected!")
+            print("Unknown Option Selected!\n")
             continue

setup.py

@@ -1,20 +1,19 @@
 #!/usr/bin/env python
 
 from setuptools import setup, find_packages
-import os
-
-# here = os.path.abspath(os.path.dirname(__file__))
 
 
 setup(name='hdpg1d',
-      version='2.0',
+      version='3.0',
       description='An 1D finite element solver using hybridizable discontinuous\
       Petrov-Galerkin method',
       author='Keyi Ni',
       author_email='[email protected]',
-      url='test',
+      url='https://github.com/asdf123101/HDPG1D',
       license='MIT',
       packages=find_packages(),
+      data_files=[('config', ['hdpg1d/config/config.json'])],
+      include_package_data=True,
       entry_points={
           'console_scripts': [
               'PGsolve = hdpg1d.cmd:main'

test/test_solve_test_problems.py

@@ -0,0 +1,46 @@
+# content of test_module.py
+import pytest
+import math
+import matplotlib
+matplotlib.use('Agg')           # supress figures in the following modules
+from hdpg1d.preprocess 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
+
+    def test_solveAdaptive(self, coeffGen):
+        coeffTest, expected = coeffGen
+        hdpgTest = hdpg1d(coeffTest)
+        hdpgTest.adaptive()
+        # get the target function value
+        # and compare to the expected value
+        soln = hdpgTest.trueErrorList[1][-1]
+        assert math.isclose(soln, expected, rel_tol=1e-5, abs_tol=1e-10)
+
+    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 *
+                                   hdpgTest.numBasisFuncs - 1]
+        assert math.isclose(soln, expected, rel_tol=1e-5, abs_tol=1e-10)