Corrected convection, reaction constant.

hdpg1d/adaptation.py

@@ -21,6 +21,7 @@ class hdpg1d(object):
         self.numBasisFuncs = coeff.pOrder + 1
         self.coeff = coeff
         self.mesh = np.linspace(0, 1, self.numEle + 1)
+        self.enrichOrder = 1
         self.primalSoln = None
         self.adjointSoln = None
         self.estErrorList = [[], []]
@@ -34,8 +35,8 @@ class hdpg1d(object):
 
     def plotState(self, counter):
         """Plot solution u with smooth higher oredr quadrature"""
-        uSmooth = np.array([])
-        uNode = np.zeros(self.numEle + 1)
+        stateSmooth = np.array([])
+        stateNode = np.zeros(self.numEle + 1)
         xSmooth = np.array([])
         gradState, _ = self.separateSoln(self.primalSoln)
         halfLenState = int(len(gradState) / 2)
@@ -47,19 +48,18 @@ class hdpg1d(object):
         for j in range(1, self.numEle + 1):
             xSmooth = np.hstack((xSmooth, (self.mesh[(j - 1)] + self.mesh[j]) / 2 + (
                 self.mesh[j] - self.mesh[j - 1]) / 2 * xi))
-            uSmooth = np.hstack(
-                (uSmooth, shp.T.dot(state[(j - 1) * self.numBasisFuncs:j * self.numBasisFuncs])))
-            uNode[j - 1] = state[(j - 1) * self.numBasisFuncs]
-        uNode[-1] = state[-1]
+            stateSmooth = np.hstack(
+                (stateSmooth, shp.T.dot(state[(j - 1) * self.numBasisFuncs:j * self.numBasisFuncs])))
+            stateNode[j - 1] = state[(j - 1) * self.numBasisFuncs]
+        stateNode[-1] = state[-1]
         plt.figure(1)
-        plt.plot(xSmooth, uSmooth, '-', color='C3')
-        plt.plot(self.mesh, uNode, 'C3.')
+        plt.plot(xSmooth, stateSmooth, '-', color='C3')
+        plt.plot(self.mesh, stateNode, 'C3.')
         plt.xlabel('$x$', fontsize=17)
         plt.ylabel('$u$', fontsize=17)
         # plt.axis([-0.05, 1.05, 0, 1.3])
         plt.grid()
         plt.pause(5e-1)
-        plt.clf()
 
     def meshAdapt(self, index):
         """Given the index list, adapt the mesh"""
@@ -79,7 +79,7 @@ class hdpg1d(object):
             matLocal = discretization(self.coeff, self.mesh)
         matGroup = matLocal.matGroup()
         A, B, _, C, D, E, F, G, H, L, R = matGroup
-        # solve
+        # solve by exploiting the local global separation
         K = -cat((C.T, G), axis=1)\
             .dot(np.linalg.inv(np.bmat([[A, -B], [B.T, D]]))
                  .dot(cat((C, E)))) + H
@@ -110,24 +110,23 @@ class hdpg1d(object):
         LHS = np.bmat([[A, -B, C],
                        [B.T, D, E],
                        [C.T, G, H]])
-        # solve
+        # solve in one shoot
         soln = np.linalg.solve(LHS.T, cat((R, F, L)))
         self.adjointSoln = soln
 
-    def residual(self):
-        enrich = 1
+    def DWResidual(self):
         if 'matResidual' in locals():
             matResidual.mesh = self.mesh
         else:
-            matResidual = discretization(self.coeff, self.mesh, enrich)
+            matResidual = discretization(
+                self.coeff, self.mesh, self.enrichOrder)
         matGroup = matResidual.matGroup()
-        A, B, BonU, C, D, E, F, G, H, L, R = matGroup
+        A, B, BonQ, C, D, E, F, G, H, L, R = matGroup
         LHS = np.bmat([[A, -B, C],
-                       [BonU, D, E]])
+                       [BonQ, D, E]])
         RHS = cat((R, F))
         residual = np.zeros(self.numEle)
-        numEnrich = self.numBasisFuncs + enrich
-        primalGradState, primalStateFace = self.separateSoln(self.primalSoln)
+        numEnrich = self.numBasisFuncs + self.enrichOrder
         adjointGradState, adjointStateFace = self.separateSoln(
             self.adjointSoln)
         for i in np.arange(self.numEle):
@@ -150,28 +149,36 @@ class hdpg1d(object):
         return np.abs(np.sum(residual)), refineIndex
 
     def adaptive(self):
-        tol = 1e-10
+        TOL = 1e-10
         estError = 10
-        counter = 0
-        ceilCounter = 30
-        while estError > tol and counter < ceilCounter:
-            print("Iteration {}. Target function error {:.3e}.".format(
-                counter, estError))
+        nodeCount = 0
+        maxCount = 30
+        while estError > TOL and nodeCount < maxCount:
             # solve
             self.solveLocal()
             self.solveAdjoint()
             # plot the solution at certain counter
-            if counter in [0, 4, 9, 19]:
-                self.plotState(counter)
+            if nodeCount in [0, 4, 9, 19, maxCount]:
+                plt.clf()
+                self.plotState(nodeCount)
             # record error
             self.trueErrorList[0].append(self.numEle)
             self.trueErrorList[1].append(
                 self.primalSoln[self.numEle * self.numBasisFuncs - 1])
-            estError, index = self.residual()
+            estError, index = self.DWResidual()
             self.estErrorList[0].append(self.numEle)
             self.estErrorList[1].append(estError)
             # adapt
             index = index.tolist()
             self.meshAdapt(index)
             self.numEle = self.numEle + len(index)
-            counter += 1
+            nodeCount += 1
+            print("Iteration {}. 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,28 +1,33 @@
 class coefficients:
-    def __init__(self, diff, conv, flux, pOrder, numEle, tauPlus, tauMinus):
-        self.diffusion = diff
-        self.convection = conv
-        self.flux = flux
+    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
+        self.TAUPLUS = tauPlus
+        self.TAUMINUS = tauMinus
 
     @classmethod
     def fromInput(cls):
         while True:
             try:
                 print("Please provide the following coefficients.")
-                diff = float(input("Diffusion coefficient: "))
-                conv = float(input("Covection coefficient: "))
-                flux = float(input("Flux: "))
-                pOrder = int(input("Order of polynomials: "))
-                numEle = int(input("Number of elements: "))
-                tauPlus = float(input("Stablization parameter plus: "))
-                tauMinus = float(input("Stablization parameter minus: "))
+                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, flux, pOrder, numEle, tauPlus, tauMinus)
+        return cls(diff, conv, reaction, pOrder, numEle, tauPlus, tauMinus)

hdpg1d/postprocess.py

@@ -6,17 +6,25 @@ import numpy as np
 
 
 class utils(object):
+    exactNumEle = 300
+    exactBasisFuncs = 5
+
     def __init__(self, solution):
         self.solution = solution
-        exactNumEle = 500
-        exactBasisFuncs = 5
-        self.solution.coeff.numEle = exactNumEle
-        self.solution.coeff.pOrder = exactBasisFuncs - 1
-        self.solution.mesh = np.linspace(0, 1, exactNumEle + 1)
+        self.numEle = self.exactNumEle
+        self.exactSoln = self.solveExact()
+
+    def solveExact(self):
+        self.solution.coeff.numEle = self.exactNumEle
+        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.exactSol = self.solution.solveLocal()[0][exactNumEle * exactBasisFuncs - 1]
+        self.solution.solveLocal()
+        exactSoln = self.solution.separateSoln(self.solution.primalSoln)[0][
+            self.exactNumEle * self.exactBasisFuncs - 1]
         # for the reaction diffusion test problem, we know the exact solution
-        self.exactSol = np.sqrt(self.solution.coeff.diffusion)
+        # self.exactSol = np.sqrt(self.solution.coeff.diffusion)
+        return exactSoln
 
     def errorL2(self):
         errorL2 = 0.
@@ -27,7 +35,8 @@ class utils(object):
         self.solution.mesh = np.linspace(0, 1, numEle + 1)
         self.solution.solveLocal()
         gradState, _ = self.solution.separateSoln(self.solution.primalSoln)
-        errorL2 = np.abs(gradState[numBasisFuncs * numEle - 1] - self.exactSol)
+        errorL2 = np.abs(
+            gradState[numBasisFuncs * numEle - 1] - self.exactSoln)
         return errorL2
 
     def uniConv(self):
@@ -44,9 +53,12 @@ class utils(object):
 
     def convHistory(self):
         """Plot the uniform and adaptive convergence history"""
+        print("Please note that the error is calculated using {} "
+              "elements with polynomial order {}."
+              .format(self.exactNumEle, self.exactBasisFuncs))
         plt.figure(2)
         trueErrorList = self.solution.trueErrorList
-        trueErrorList[1] = np.abs(trueErrorList[1] - self.exactSol)
+        trueErrorList[1] = np.abs(trueErrorList[1] - self.exactSoln)
         estErrorList = self.solution.estErrorList
         plt.loglog(trueErrorList[0],
                    trueErrorList[1], '-ro')

hdpg1d/preprocess.py

@@ -5,7 +5,7 @@ from scipy.linalg import block_diag
 
 def shape(x, p):
     """generate p shape functions and its first order derivative
-    (order p-1) at the given location x"""
+    at the given location x. x can be an array"""
     A = np.array([np.linspace(-1, 1, p)]).T**np.arange(p)
     C = np.linalg.inv(A).T
     x = np.array([x]).T
@@ -19,8 +19,7 @@ def forcing(x):
     return f
 
 
-def boundaryCondition(case, t=None):
-    # boundary condition
+def boundaryCondition(case):
     if case == 0:
         # primal problem
         bc = [0, 0]
@@ -31,7 +30,8 @@ def boundaryCondition(case, t=None):
 
 
 class discretization(object):
-    """Given the problem statement, construct the discretization matrices"""
+    """Given the problem statement and current mesh,
+    construct the discretization matrices"""
 
     def __init__(self, coeff, mesh, enrich=None):
         self.mesh = mesh
@@ -39,11 +39,11 @@ class discretization(object):
         self.enrich = enrich
         # the following init are for the sake of simplicity
         self.numBasisFuncs = coeff.pOrder + 1
-        self.tau_pos = coeff.tauPlus
-        self.tau_neg = coeff.tauMinus
-        self.conv = coeff.convection
-        self.kappa = coeff.diffusion
-        self.n_ele = len(mesh) - 1
+        self.tau_pos = coeff.TAUPLUS
+        self.tau_neg = coeff.TAUMINUS
+        self.conv = coeff.CONVECTION
+        self.kappa = coeff.DIFFUSION
+        self.numEle = len(mesh) - 1
         self.dist = self.distGen()
         # shape function and gauss quadrature
         self.gqOrder = 5 * self.numBasisFuncs
@@ -55,8 +55,8 @@ class discretization(object):
                 self.xi, self.numBasisFuncs + enrich)
 
     def distGen(self):
-        dist = np.zeros(self.n_ele)
-        for i in range(1, self.n_ele + 1):
+        dist = np.zeros(self.numEle)
+        for i in range(1, self.numEle + 1):
             dist[i - 1] = self.mesh[i] - self.mesh[i - 1]
         return dist
 
@@ -82,8 +82,8 @@ class discretization(object):
             numBasisFuncs = self.numBasisFuncs
             shp = self.shp
         # forcing vector F
-        F = np.zeros(numBasisFuncs * self.n_ele)
-        for i in range(1, self.n_ele + 1):
+        F = np.zeros(numBasisFuncs * self.numEle)
+        for i in range(1, self.numEle + 1):
             f = self.dist[i - 1] / 2 \
                 * shp.dot(self.wi * forcing(self.mesh[i - 1] +
                                             1 / 2 * (1 + self.xi) *
@@ -92,9 +92,9 @@ class discretization(object):
         F[0] += (self.conv + self.tau_pos) * boundaryCondition(0)[0]
         F[-1] += (-self.conv + self.tau_neg) * boundaryCondition(0)[1]
         # L, easy in 1d
-        L = np.zeros(self.n_ele - 1)
+        L = np.zeros(self.numEle - 1)
         # R, easy in 1d
-        R = np.zeros(numBasisFuncs * self.n_ele)
+        R = np.zeros(numBasisFuncs * self.numEle)
         R[0] = boundaryCondition(0)[0]
         R[-1] = -boundaryCondition(0)[1]
         lhsMat = namedtuple('lhsMat', ['F', 'L', 'R'])
@@ -108,29 +108,35 @@ class discretization(object):
             shpxEnrich = self.shpxEnrich
             b = (self.shpx.T * np.ones((self.gqOrder, self.numBasisFuncs))
                  ).T.dot(np.diag(self.wi).dot(shpEnrich.T))
-            # BonU is only used in calculating DWR residual
-            BonU = block_diag(*[b] * (self.n_ele)).T
+            # BonQ is only used in calculating DWR residual
+            BonQ = block_diag(*[b] * (self.numEle)).T
         else:
             numBasisFuncsEnrich = self.numBasisFuncs
             shpEnrich = self.shp
             shpxEnrich = self.shpx
-            BonU = None
-        a = 1 / self.kappa * shpEnrich.dot(np.diag(self.wi).dot(self.shp.T))
+            BonQ = None
+        a = 1 / self.coeff.DIFFUSION * \
+            shpEnrich.dot(np.diag(self.wi).dot(self.shp.T))
         A = np.repeat(self.dist, self.numBasisFuncs) / \
-            2 * block_diag(*[a] * (self.n_ele))
+            2 * block_diag(*[a] * (self.numEle))
         b = (shpxEnrich.T * np.ones((self.gqOrder, numBasisFuncsEnrich))
              ).T.dot(np.diag(self.wi).dot(self.shp.T))
-        B = block_diag(*[b] * (self.n_ele))
-        d = shpEnrich.dot(np.diag(self.wi).dot(self.shp.T))
+        B = block_diag(*[b] * (self.numEle))
+        d = self.coeff.REACTION * \
+            shpEnrich.dot(np.diag(self.wi).dot(self.shp.T))
         # assemble D
-        d_face = np.zeros((numBasisFuncsEnrich, self.numBasisFuncs))
-        d_face[0, 0] = self.tau_pos
-        d_face[-1, -1] = self.tau_neg
+        dFace = np.zeros((numBasisFuncsEnrich, self.numBasisFuncs))
+        dFace[0, 0] = self.tau_pos
+        dFace[-1, -1] = self.tau_neg
+        dConv = -self.conv * (shpxEnrich.T * np.ones((self.gqOrder,
+                                                      numBasisFuncsEnrich)))\
+            .T.dot(np.diag(self.wi).dot(self.shp.T))
         D = np.repeat(self.dist, self.numBasisFuncs) / 2 * \
-            block_diag(*[d] * (self.n_ele)) + \
-            block_diag(*[d_face] * (self.n_ele))
+            block_diag(*[d] * (self.numEle)) + \
+            block_diag(*[dFace] * (self.numEle)) +\
+            block_diag(*[dConv] * (self.numEle))
         eleMat = namedtuple('eleMat', ['A', 'B', 'BonU', 'D'])
-        return eleMat(A, B, BonU, D)
+        return eleMat(A, B, BonQ, D)
 
     def interfaceGen(self):
         """Generate matrices associated with interfaces"""
@@ -145,68 +151,68 @@ class discretization(object):
         h = np.zeros((2, 2))
         h[0, 0], h[-1, -1] = -conv - tau_pos, conv - tau_neg
         # mappinng matrix
-        map_h = np.zeros((2, self.n_ele), dtype=int)
+        map_h = np.zeros((2, self.numEle), dtype=int)
         map_h[:, 0] = np.arange(2)
-        for i in np.arange(1, self.n_ele):
+        for i in np.arange(1, self.numEle):
             map_h[:, i] = np.arange(
                 map_h[2 - 1, i - 1], map_h[2 - 1, i - 1] + 2)
         # assemble H and eliminate boundaries
-        H = np.zeros((self.n_ele + 1, self.n_ele + 1))
-        for i in range(self.n_ele):
+        H = np.zeros((self.numEle + 1, self.numEle + 1))
+        for i in range(self.numEle):
             for j in range(2):
                 m = map_h[j, i]
                 for k in range(2):
                     n = map_h[k, i]
                     H[m, n] += h[j, k]
-        H = H[1:self.n_ele][:, 1:self.n_ele]
+        H = H[1:self.numEle][:, 1:self.numEle]
 
         # elemental e
         e = np.zeros((numBasisFuncs, 2))
         e[0, 0], e[-1, -1] = -conv - tau_pos, conv - tau_neg
         # mapping matrix
-        map_e_x = np.arange(numBasisFuncs * self.n_ele,
-                            dtype=int).reshape(self.n_ele,
+        map_e_x = np.arange(numBasisFuncs * self.numEle,
+                            dtype=int).reshape(self.numEle,
                                                numBasisFuncs).T
         map_e_y = map_h
         # assemble global E
-        E = np.zeros((numBasisFuncs * self.n_ele, self.n_ele + 1))
-        for i in range(self.n_ele):
+        E = np.zeros((numBasisFuncs * self.numEle, self.numEle + 1))
+        for i in range(self.numEle):
             for j in range(numBasisFuncs):
                 m = map_e_x[j, i]
                 for k in range(2):
                     n = map_e_y[k, i]
                     E[m, n] += e[j, k]
-        E = E[:, 1:self.n_ele]
+        E = E[:, 1:self.numEle]
 
         # elemental c
         c = np.zeros((numBasisFuncs, 2))
         c[0, 0], c[-1, -1] = -1, 1
         # assemble global C
-        C = np.zeros((numBasisFuncs * self.n_ele, self.n_ele + 1))
-        for i in range(self.n_ele):
+        C = np.zeros((numBasisFuncs * self.numEle, self.numEle + 1))
+        for i in range(self.numEle):
             for j in range(numBasisFuncs):
                 m = map_e_x[j, i]
                 for k in range(2):
                     n = map_e_y[k, i]
                     C[m, n] += c[j, k]
-        C = C[:, 1:self.n_ele]
+        C = C[:, 1:self.numEle]
 
         # elemental g
         g = np.zeros((2, numBasisFuncs))
         g[0, 0], g[-1, -1] = tau_pos, tau_neg
         # mapping matrix
         map_g_x = map_h
-        map_g_y = np.arange(numBasisFuncs * self.n_ele,
-                            dtype=int).reshape(self.n_ele,
+        map_g_y = np.arange(numBasisFuncs * self.numEle,
+                            dtype=int).reshape(self.numEle,
                                                numBasisFuncs).T
         # assemble global G
-        G = np.zeros((self.n_ele + 1, numBasisFuncs * self.n_ele))
-        for i in range(self.n_ele):
+        G = np.zeros((self.numEle + 1, numBasisFuncs * self.numEle))
+        for i in range(self.numEle):
             for j in range(2):
                 m = map_g_x[j, i]
                 for k in range(numBasisFuncs):
                     n = map_g_y[k, i]
                     G[m, n] += g[j, k]
-        G = G[1:self.n_ele, :]
+        G = G[1:self.numEle, :]
         faceMat = namedtuple('faceMat', ['C', 'E', 'G', 'H'])
         return faceMat(C, E, G, H)

hdpg1d/solve.py

@@ -1,7 +1,6 @@
 from .coefficients import coefficients
 from .adaptation import hdpg1d
 from .postprocess import utils
-import sys
 
 
 def queryYesNo(question, default="yes"):
@@ -17,22 +16,22 @@ def queryYesNo(question, default="yes"):
         raise ValueError("invalid default answer: '%s'" % default)
 
     while True:
-        sys.stdout.write(question + prompt)
+        print(question + prompt)
         choice = input().lower()
         if default is not None and choice == '':
             return valid[default]
         elif choice in valid:
             return valid[choice]
         else:
-            sys.stdout.write("Please respond with 'yes' or 'no' "
-                             "(or 'y' or 'n').\n")
+            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, 0, 2, 2, 1e-6, 1e-6)
+        Coeff = coefficients(1e-6, 0, 1, 2, 2, 1, 1)
     else:
         Coeff = coefficients.fromInput()
     return Coeff
@@ -61,11 +60,15 @@ def hdgSolve():
 def runInteractive():
     menu()
     selection = input("Please Select: ")
-    if selection == '1':
-        hdgSolve()
-    elif selection == '2':
-        print("In development...")
-    elif selection == '3':
-        print("Bye.")
-    else:
-        print("Unknown Option Selected!")
+    while True:
+        if selection == '1':
+            hdgSolve()
+            break
+        elif selection == '2':
+            print("In development...")
+        elif selection == '3':
+            print("Bye.")
+            break
+        else:
+            print("Unknown Option Selected!")
+            continue