Added numerically approximated exact solution

hdpg1d/discretization.py

@@ -4,23 +4,20 @@ from scipy.linalg import block_diag
 import matplotlib.pyplot as plt
 
 
-class HDPG1d(object):
+class hdpg1d(object):
     """
     1d advection hdg solver outlined in 'an implicit HHDG method for
     confusion'. Test case: /tau = 1, convection only, linear and higher order.
     Please enter number of elements and polynomial order, i.e., HDG1d(10,2)
     """
-    # stablization parameter
-    tau_pos = 1e-6
-    tau_neg = 1e-6
-    # advection constant
-    c = 0
-    # diffusion constant
-    kappa = 1e-6
 
     def __init__(self, n_ele, p):
         self.n_ele = n_ele
         self.p = p + 1
+        self.tau_pos = 1e-6
+        self.tau_neg = 1e-6
+        self.c = 0
+        self.kappa = 1e-6
         self.estError = []
         self.trueError = []
 
@@ -66,11 +63,12 @@ class HDPG1d(object):
         dx = x_c[n_ele] - x_c[n_ele - 1]
         return x_i, dx, x_c
 
-    def uexact(self, x):
-        # Ue = self.bc(0) + np.sin(2*pi*x)
-        Ue = (np.exp(1 / self.kappa * (x - 1)) - 1) / \
-            (np.exp(-1 / self.kappa) - 1)
-        return Ue
+    def exact(self, x):
+        """solve the problem in an enriched space to simulate exact soltuion"""
+        self.n_ele = 1000
+        self.p = 3
+        x = np.linspace(0, 1, self.n_ele + 1)
+        self.exactSol = self.solve_local([], x)
 
     def matrix_gen(self, index, x):
         n_ele = self.n_ele

hdpg1d/postprocess.py

@@ -5,38 +5,50 @@ import matplotlib.pyplot as plt
 import numpy as np
 
 
-def errorL2(solution):
-    errorL2 = 0.
-    # solve the uniform case
-    x = np.linspace(0, 1, solution.n_ele + 1)
-    U, _ = solution.solve_local([], x)
-    errorL2 = np.abs(U[solution.p * solution.n_ele - 1] -
-                     np.sqrt(solution.kappa))
-    return errorL2
+class utils(object):
+    def __init__(self, solution):
+        self.solution = solution
+        exactNumEle = 200
+        exactPolyOrder = 5
+        self.solution.n_ele = exactNumEle
+        self.solution.p = exactPolyOrder
+        x = np.linspace(0, 1, exactNumEle + 1)
+        self.exactSol = self.solution.solve_local(
+            [], x)[0].A1[exactNumEle * exactPolyOrder - 1]
 
+    def errorL2(self):
+        errorL2 = 0.
+        n_ele = self.solution.n_ele
+        p = self.solution.p
+        # solve the uniform case
+        x = np.linspace(0, 1, n_ele + 1)
+        U, _ = self.solution.solve_local([], x)
+        errorL2 = np.abs(U[p * n_ele - 1] - self.exactSol)
+        return errorL2
 
-def uniConv(solution):
-    p = np.arange(2, 3)
-    n_ele = 2**np.arange(1, 9)
-    uniError = np.zeros((n_ele.size, p.size))
-    for i in range(p.size):
-        solution.p = p[i]
-        for j, n in enumerate(n_ele):
-            solution.n_ele = n
-            uniError[j, i] = errorL2(solution)
-    return n_ele, uniError
+    def uniConv(self):
+        p = np.arange(2, 3)
+        n_ele = 2**np.arange(1, 9)
+        uniError = np.zeros((n_ele.size, p.size))
+        for i in range(p.size):
+            self.solution.p = p[i]
+            for j, n in enumerate(n_ele):
+                self.solution.n_ele = n
+                uniError[j, i] = self.errorL2()
+        return n_ele, uniError
 
-
-def convHistory(solution):
-    plt.loglog(solution.trueError[0, 0:-1],
-               solution.trueError[1, 0:-1], '-ro')
-    # plt.axis([1, 250, 1e-13, 1e-2])
-    n_ele, errorL2 = uniConv(solution)
-    plt.loglog(n_ele, errorL2, '-o')
-    plt.loglog(solution.estError[0, :],
-               solution.estError[1, :], '--', color='#1f77b4')
-    plt.xlabel('Number of elements', fontsize=17)
-    plt.ylabel('Error', fontsize=17)
-    plt.grid()
-    plt.legend(('Adaptive', 'Uniform', 'Estimator'), loc=3, fontsize=15)
-    plt.show()
+    def convHistory(self):
+        trueError = self.solution.trueError
+        estError = self.solution.estError
+        plt.loglog(trueError[0, 0:-1],
+                   trueError[1, 0:-1], '-ro')
+        # plt.axis([1, 250, 1e-13, 1e-2])
+        n_ele, errorL2 = self.uniConv()
+        plt.loglog(n_ele, errorL2, '-o')
+        plt.loglog(estError[0, :],
+                   estError[1, :], '--', color='#1f77b4')
+        plt.xlabel('Number of elements', fontsize=17)
+        plt.ylabel('Error', fontsize=17)
+        plt.grid()
+        plt.legend(('Adaptive', 'Uniform', 'Estimator'), loc=3, fontsize=15)
+        plt.show()

hdpg1d/preprocess.py

@@ -0,0 +1 @@
+def mesh:

hdpg1d/solve.py

@@ -1,6 +1,6 @@
 from .coefficients import coefficients
-from .discretization import HDPG1d
-from .postprocess import convHistory
+from .discretization import hdpg1d
+from .postprocess import utils
 import sys
 
 
@@ -50,10 +50,11 @@ def menu():
     selection = input("Please Select: ")
     if selection == '1':
         hdgCoeff = getCoefficients()
-        hdgSolution = HDPG1d(hdgCoeff.nele, hdgCoeff.porder)
+        hdgSolution = hdpg1d(hdgCoeff.nele, hdgCoeff.porder)
         # solve the problem adaptively and plot convergence history
         hdgSolution.adaptive()
-        convHistory(hdgSolution)
+        post = utils(hdgSolution)
+        post.convHistory()
     elif selection == '2':
         print("In development...")
     elif selection == '3':