Optimized algorithm

hdpg1d/adaptation.py

@@ -90,12 +90,12 @@ class hdpg1d(object):
             .dot(np.linalg.inv(np.bmat([[A, -B], [B.T, D]])))\
             .dot(np.array([cat((R, F))]).T)
 
-        def M_x(vec):
+        def invRHS(vec):
             """Construct preconditioner"""
             matVec = spla.spsolve(sK, vec)
             return matVec
         n = len(F_hat)
-        preconditioner = spla.LinearOperator((n, n), M_x)
+        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(
@@ -115,7 +115,7 @@ 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],
@@ -124,14 +124,14 @@ class hdpg1d(object):
         RHS = cat((R, F, L))
 
         # solve in one shoot using GMRES
-
-        def M_x(vec):
+        def invRHS(vec):
             """Construct preconditioner"""
             matVec = spla.spsolve(sLHS, vec)
             return matVec
         n = len(RHS)
-        preconditioner = spla.LinearOperator((n, n), M_x)
+        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):

hdpg1d/coefficients.py

@@ -1,34 +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:
-                print("Something is wrong. Exit.")
-                break
-        return cls(diff, conv, reaction, pOrder, numEle, tauPlus, tauMinus)

hdpg1d/config.json

@@ -1,19 +1,17 @@
 {
-    "maps": [
-        {
-            "id": "blabla",
-            "iscategorical": "0"
-        },
-        {
-            "id": "blabla",
-            "iscategorical": "0"
-        }
-    ],
-    "masks": {
-        "id": "valore"
-    },
-    "om_points": "value",
-    "parameters": {
-        "id": "valore"
-    }
+	"coefficients": {
+		"diffusion": 1e-6,
+		"convection": 0,
+		"reaction": 1,
+		"pOrder": 2,
+		"numEle": 2,
+		"tauPlus": 1,
+		"tauMinus": 1
+	},
+	"forcing": "x",
+	"flux": "x",
+	"boundary": {
+		"left": 0,
+		"right": 0
+	}
 }

hdpg1d/preprocess.py

@@ -1,16 +1,83 @@
+import os
 import json
+import ast
+import operator as op
 import numpy as np
-import os
 from collections import namedtuple
 from scipy.linalg import block_diag
 
 for loc in os.curdir, os.path.expanduser("~"), "/etc/hdpg1d":
     try:
         with open(os.path.join(loc, "config.json")) as source:
-            data = json.load(source)
+            configdata = json.load(source)
     except IOError:
         pass
 
+# supported 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"]
+        coeff = coefficients(diffDefault, convDefault, reactionDefault,
+                             pOrderDefault, numEleDefault,
+                             tauPlusDefault, tauMinusDefault)
+    else:
+        coeff = coefficients.fromInput()
+    return coeff
+
 
 def shape(x, p):
     """generate p shape functions and its first order derivative
@@ -24,20 +91,61 @@ def shape(x, p):
 
 
 def forcing(x):
-    f = 1
+    f = np.zeros(len(x))
+    for i, forcingItem in enumerate(x):
+        forcingExpr = configdata["forcing"]
+        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):
+        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:
+                print("Something is wrong. Exit.")
+                break
+        return cls(diff, conv, reaction, pOrder, numEle, tauPlus, tauMinus)
+
+
 class discretization(object):
     """Given the problem statement and current mesh,
     construct the discretization matrices"""
@@ -98,14 +206,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

tests/test_solve_test_problems.py

@@ -3,7 +3,7 @@ import pytest
 import math
 import matplotlib
 matplotlib.use('Agg')           # supress figures in the following modules
-from hdpg1d.coefficients import coefficients as coeff
+from hdpg1d.preprocess import coefficients as coeff
 from hdpg1d.adaptation import hdpg1d