Improved performance

hdpg1d/adaptation.py

@@ -2,7 +2,7 @@ import numpy as np
 from numpy import concatenate as cat
 import matplotlib.pyplot as plt
 import warnings
-from .preprocess import shape, discretization
+from .preprocess import shape, discretization, boundaryCondition
 
 plt.rc('text', usetex=True)
 plt.rc('font', family='serif')
@@ -28,24 +28,6 @@ class hdpg1d(object):
         self.estErrorList = [[], []]
         self.trueErrorList = [[], []]
 
-    def bc(self, case, t=None):
-        # boundary condition
-        if case == 0:
-            # advection-diffusion
-            bc = [0, 0]
-        if case == 1:
-            # simple convection
-            # bc = np.sin(2*np.pi*t)
-            # adjoint boundary
-            bc = [0, 1]
-        return bc
-
-    def forcing(self, x):
-        # f = np.cos(2*np.pi*x)
-        # f = 4*pi**2*sin(2*pi*x)
-        f = 1
-        return f
-
     def plotU(self, counter):
         """Plot solution u with smooth higher oredr quadrature"""
         uSmooth = np.array([])
@@ -90,9 +72,8 @@ class hdpg1d(object):
             matLocal.mesh = self.mesh
         else:
             matLocal = discretization(self.coeff, self.mesh)
-        F, L, R = matLocal.lhsGen()
-        A, B, D = matLocal.eleGen()
-        C, E, G, H = matLocal.interfaceGen()
+        matGroup = matLocal.matGroup()
+        A, B, _, C, D, E, F, G, H, L, R = matGroup
         # solve
         K = -cat((C.T, G), axis=1)\
             .dot(np.linalg.inv(np.bmat([[A, -B], [B.T, D]]))
@@ -104,7 +85,6 @@ class hdpg1d(object):
         u = np.linalg.inv(np.bmat([[A, -B], [B.T, D]]))\
             .dot(np.array([np.concatenate((R, F))]).T -
                  cat((C, E)).dot(uFace))
-        # self.u = u.A1
         return u.A1, uFace.A1
 
     def solveAdjoint(self):
@@ -116,182 +96,51 @@ class hdpg1d(object):
         else:
             matAdjoint = discretization(self.coeff, self.mesh)
         self.coeff.pOrder = self.coeff.pOrder - 1
-        F, L, R = matAdjoint.lhsGen()
-        A, B, D = matAdjoint.eleGen()
-        C, E, G, H = matAdjoint.interfaceGen()
+        matGroup = matAdjoint.matGroup()
+        A, B, _, C, D, E, F, G, H, L, R = matGroup
         # add adjoint LHS conditions
         F = np.zeros(len(F))
-        R[-1] = -self.bc(1)[1]
+        R[-1] = -boundaryCondition(1)[1]
 
         # assemble global matrix LHS
         LHS = np.bmat([[A, -B, C],
                        [B.T, D, E],
                        [C.T, G, H]])
         # solve
-        U = np.linalg.solve(LHS.T, np.concatenate(
-            (R, F, L)))
+        U = np.linalg.solve(LHS.T, cat((R, F, L)))
         return U[0:2 * len(C)], U[len(C):len(U)]
 
     def residual(self, U, hat_U, z, hat_z):
+        enrich = 1
         if 'matResidual' in locals():
-            # if matLocal exists,
-            # only change the mesh instead of initializing again
             matResidual.mesh = self.mesh
         else:
-            matResidual = discretization(self.coeff, self.mesh, 1)
-        F, L, R = matResidual.lhsGen()
-        # A, B, D = matResidual.eleGen()
-        # C, E, G, H = matResidual.interfaceGen()
-
-        numEle = self.numEle
-        p = self.numBasisFuncs + 1
-
-        p_l = p - 1
-        # order of gauss quadrature
-        gorder = 2 * p
-        # shape function and gauss quadrature
-        xi, wi = np.polynomial.legendre.leggauss(gorder)
-        shp, shpx = shape(xi, p)
-        shp_l, shpx_l = shape(xi, p_l)
-        # ---------------------------------------------------------------------
-        # advection constant
-        con = self.c
-
-        # diffusion constant
-        kappa = self.kappa
-        z_q, z_u, z_hat = np.zeros(p * numEle), \
-            np.zeros(p *
-                     numEle), np.zeros(numEle - 1)
-
-        q, u, lamba = np.zeros(p_l * numEle), \
-            np.zeros(p_l * numEle), np.zeros(numEle - 1)
-        for i in np.arange(p * numEle):
-            z_q[i] = z[i]
-            z_u[i] = z[i + p * numEle]
-
-        for i in np.arange(p_l * numEle):
-            q[i] = U[i]
-            u[i] = U[i + p_l * numEle]
-
-        for i in np.arange(numEle - 1):
-            z_hat[i] = hat_z[i]
-
-        # add boundary condtions to U_hat
-        U_hat = np.zeros(numEle + 1)
-        for i, x in enumerate(hat_U):
-            U_hat[i + 1] = x
-        U_hat[0] = self.bc(0)[0]
-        U_hat[-1] = self.bc(0)[1]
-
-        # L, easy in 1d
-        # L = np.zeros(numEle + 1)
-
-        # R, easy in 1d
-        RR = np.zeros(p * numEle)
-
-        # elemental forcing vector
-        dist = np.zeros(numEle)
-        for i in range(1, numEle + 1):
-            dist[i - 1] = self.mesh[i] - self.mesh[i - 1]
-
-        # elemental h
-        h = np.zeros((2, 2))
-        h[0, 0], h[-1, -1] = -con - self.tau_pos, con - self.tau_neg
-        # mappinng matrix
-        map_h = np.zeros((2, numEle), dtype=int)
-        map_h[:, 0] = np.arange(2)
-        for i in np.arange(1, 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((numEle + 1, numEle + 1))
-        for i in range(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:numEle][:, 1:numEle]
-
-        # elemental g
-        g = np.zeros((2, p_l))
-        g[0, 0], g[-1, -1] = self.tau_pos, self.tau_neg
-        # mapping matrix
-        map_g_x = map_h
-        map_g_y = np.arange(p_l * numEle, dtype=int).reshape(numEle, p_l).T
-        # assemble global G
-        G = np.zeros((numEle + 1, p_l * numEle))
-        for i in range(numEle):
-            for j in range(2):
-                m = map_g_x[j, i]
-                for k in range(p_l):
-                    n = map_g_y[k, i]
-                    G[m, n] += g[j, k]
-        G = G[1:numEle, :]
-
-        # elemental c
-        c = np.zeros((p_l, 2))
-        c[0, 0], c[-1, -1] = -1, 1
-
-        # mapping matrix
-        map_e_x = np.arange(p_l * numEle, dtype=int).reshape(numEle, p_l).T
-        map_e_y = map_h
-
-        # assemble global C
-        C = np.zeros((p_l * numEle, numEle + 1))
-        for i in range(numEle):
-            for j in range(p_l):
-                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:numEle]
-
-        # L, easy in 1d
-        # L = np.zeros(numEle - 1)
-
-        # residual vector
-        Residual = np.zeros(self.numEle)
+            matResidual = discretization(self.coeff, self.mesh, enrich)
+        matGroup = matResidual.matGroup()
+        A, B, BonU, C, D, E, F, G, H, L, R = matGroup
+        LHS = np.bmat([[A, -B, C],
+                       [BonU, D, E]])
+        RHS = cat((R, F))
+        residual = np.zeros(self.numEle)
+        numEnrich = self.numBasisFuncs + enrich
         for i in np.arange(self.numEle):
-            a = dist[i] / 2 * 1 / kappa * \
-                ((shp.T).T).dot(np.diag(wi).dot(shp_l.T))
-
-            b = ((shpx.T) * np.ones((gorder, p))
-                 ).T.dot(np.diag(wi).dot(shp_l.T))
-
-            b_t = ((shpx_l.T) * np.ones((gorder, p_l))
-                   ).T.dot(np.diag(wi).dot(shp.T))
-
-            d = dist[i] / 2 * shp.dot(np.diag(wi).dot(shp_l.T))
-            d[0, 0] += self.tau_pos
-            d[-1, -1] += self.tau_neg
-
-            h = np.zeros((2, 2))
-            h[0, 0], h[-1, -1] = -con - self.tau_pos, con - self.tau_neg
-
-            g = np.zeros((2, p_l))
-            g[0, 0], g[-1, -1] = self.tau_pos, self.tau_neg
-
-            e = np.zeros((p, 2))
-            e[0, 0], e[-1, -1] = -con - self.tau_pos, con - self.tau_neg
-
-            c = np.zeros((p, 2))
-            c[0, 0], c[-1, -1] = -1, 1
-
-            m = np.zeros((2, p_l))
-            m[0, 0], m[-1, -1] = -1, 1
-            # local error
-            Residual[i] = (np.concatenate((a.dot(q[p_l * i:p_l * i + p_l]) + -b.dot(u[p_l * i:p_l * i + p_l]) + c.dot(U_hat[i:i + 2]),
-                                           b_t.T.dot(q[p_l * i:p_l * i + p_l]) + d.dot(u[p_l * i:p_l * i + p_l]) + e.dot(U_hat[i:i + 2]))) - np.concatenate((R[p * i:p * i + p], F[p * i:p * i + p]))).dot(1 - np.concatenate((z_q[p * i:p * i + p], z_u[p * i:p * i + p])))
-
-        com_index = np.argsort(np.abs(Residual))
+            primalResidual = (LHS.dot(cat((U, hat_U))) - RHS).A1
+            uLength = self.numEle * numEnrich
+            stepLength = i * numEnrich
+            uDWR = primalResidual[stepLength:stepLength + numEnrich].dot(
+                (1 - z)[stepLength:stepLength + numEnrich])
+            qDWR = primalResidual[uLength + stepLength:uLength +
+                                  stepLength + numEnrich]\
+                .dot((1 - z)[uLength + stepLength:uLength +
+                             stepLength + numEnrich])
+            residual[i] = uDWR + qDWR
+        # sort residual index
+        com_index = np.argsort(np.abs(residual))
         # select \theta% elements with the large error
         theta = 0.15
-        refine_index = com_index[int(self.numEle * (1 - theta)):len(Residual)]
-
-        # global error
-        R_g = (C.T.dot(q) + G.dot(u) + H.dot(U_hat[1:-1])).dot(1 - z_hat)
-        return np.abs(np.sum(Residual) + np.sum(R_g)), refine_index + 1
+        refine_index = com_index[
+            int(self.numEle * (1 - theta)):len(residual)] + 1
+        return np.abs(np.sum(residual)), refine_index
 
     def adaptive(self):
         tol = 1e-10

hdpg1d/preprocess.py

@@ -1,4 +1,5 @@
 import numpy as np
+from collections import namedtuple
 from scipy.linalg import block_diag
 
 
@@ -18,15 +19,13 @@ def forcing(x):
     return f
 
 
-def bc(case, t=None):
+def boundaryCondition(case, t=None):
     # boundary condition
     if case == 0:
-        # advection-diffusion
+        # primal problem
         bc = [0, 0]
     if case == 1:
-        # simple convection
-        # bc = np.sin(2*np.pi*t)
-        # adjoint boundary
+        # adjoint problem
         bc = [0, 1]
     return bc
 
@@ -61,9 +60,22 @@ class discretization(object):
             dist[i - 1] = self.mesh[i] - self.mesh[i - 1]
         return dist
 
+    def matGroup(self):
+        lhsMat = self.lhsGen()
+        F, L, R = lhsMat.F, lhsMat.L, lhsMat.R
+        eleMat = self.eleGen()
+        A, B, BonU, D = eleMat.A, eleMat.B, eleMat.BonU, eleMat.D
+        faceMat = self.interfaceGen()
+        C, E, G, H = faceMat.C, faceMat.E, faceMat.G, faceMat.H
+        matGroup = namedtuple(
+            'matGroup', ['A', 'B', 'BonU', 'C', 'D', 'E', 'F', 'G', 'H', 'L', 'R'])
+        return matGroup(A, B, BonU, C, D, E, F, G, H, L, R)
+
     def lhsGen(self):
         """Generate matrices associated with left hand side"""
         if self.enrich is not None:
+            # when enrich is given
+            # provide matrices used in the DWR residual calculation
             numBasisFuncs = self.numBasisFuncs + self.enrich
             shp = self.shpEnrich
         else:
@@ -77,36 +89,57 @@ 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) * bc(0)[0]
-        F[-1] += (-self.conv + self.tau_neg) * bc(0)[1]
+        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)
         # R, easy in 1d
         R = np.zeros(numBasisFuncs * self.n_ele)
-        R[0] = bc(0)[0]
-        R[-1] = -bc(0)[1]
-        return F, L, R
+        R[0] = boundaryCondition(0)[0]
+        R[-1] = -boundaryCondition(0)[1]
+        lhsMat = namedtuple('lhsMat', ['F', 'L', 'R'])
+        return lhsMat(F, L, R)
 
     def eleGen(self):
         """Generate matrices associated with elements"""
-        a = 1 / self.kappa * self.shp.dot(np.diag(self.wi).dot(self.shp.T))
+        if self.enrich is not None:
+            numBasisFuncsEnrich = self.numBasisFuncs + self.enrich
+            shpEnrich = self.shpEnrich
+            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
+        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))
         A = np.repeat(self.dist, self.numBasisFuncs) / \
             2 * block_diag(*[a] * (self.n_ele))
-        b = (self.shpx.T * np.ones((self.gqOrder, self.numBasisFuncs))
+        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 = self.shp.dot(np.diag(self.wi).dot(self.shp.T))
+        d = shpEnrich.dot(np.diag(self.wi).dot(self.shp.T))
         # assemble D
-        d_face = np.zeros((self.numBasisFuncs, self.numBasisFuncs))
+        d_face = np.zeros((numBasisFuncsEnrich, self.numBasisFuncs))
         d_face[0, 0] = self.tau_pos
         d_face[-1, -1] = self.tau_neg
         D = np.repeat(self.dist, self.numBasisFuncs) / 2 * \
             block_diag(*[d] * (self.n_ele)) + \
             block_diag(*[d_face] * (self.n_ele))
-        return A, B, D
+        eleMat = namedtuple('eleMat', ['A', 'B', 'BonU', 'D'])
+        return eleMat(A, B, BonU, D)
 
     def interfaceGen(self):
         """Generate matrices associated with interfaces"""
+        if self.enrich is not None:
+            # when enrich is given
+            # provide matrices used in the DWR residual calculation
+            numBasisFuncs = self.numBasisFuncs + self.enrich
+        else:
+            numBasisFuncs = self.numBasisFuncs
         tau_pos, tau_neg, conv = self.tau_pos, self.tau_neg, self.conv
         # elemental h
         h = np.zeros((2, 2))
@@ -128,17 +161,17 @@ class discretization(object):
         H = H[1:self.n_ele][:, 1:self.n_ele]
 
         # elemental e
-        e = np.zeros((self.numBasisFuncs, 2))
+        e = np.zeros((numBasisFuncs, 2))
         e[0, 0], e[-1, -1] = -conv - tau_pos, conv - tau_neg
         # mapping matrix
-        map_e_x = np.arange(self.numBasisFuncs * self.n_ele,
+        map_e_x = np.arange(numBasisFuncs * self.n_ele,
                             dtype=int).reshape(self.n_ele,
-                                               self.numBasisFuncs).T
+                                               numBasisFuncs).T
         map_e_y = map_h
         # assemble global E
-        E = np.zeros((self.numBasisFuncs * self.n_ele, self.n_ele + 1))
+        E = np.zeros((numBasisFuncs * self.n_ele, self.n_ele + 1))
         for i in range(self.n_ele):
-            for j in range(self.numBasisFuncs):
+            for j in range(numBasisFuncs):
                 m = map_e_x[j, i]
                 for k in range(2):
                     n = map_e_y[k, i]
@@ -146,12 +179,12 @@ class discretization(object):
         E = E[:, 1:self.n_ele]
 
         # elemental c
-        c = np.zeros((self.numBasisFuncs, 2))
+        c = np.zeros((numBasisFuncs, 2))
         c[0, 0], c[-1, -1] = -1, 1
         # assemble global C
-        C = np.zeros((self.numBasisFuncs * self.n_ele, self.n_ele + 1))
+        C = np.zeros((numBasisFuncs * self.n_ele, self.n_ele + 1))
         for i in range(self.n_ele):
-            for j in range(self.numBasisFuncs):
+            for j in range(numBasisFuncs):
                 m = map_e_x[j, i]
                 for k in range(2):
                     n = map_e_y[k, i]
@@ -159,21 +192,21 @@ class discretization(object):
         C = C[:, 1:self.n_ele]
 
         # elemental g
-        g = np.zeros((2, self.numBasisFuncs))
+        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(self.numBasisFuncs * self.n_ele,
+        map_g_y = np.arange(numBasisFuncs * self.n_ele,
                             dtype=int).reshape(self.n_ele,
-                                               self.numBasisFuncs).T
+                                               numBasisFuncs).T
         # assemble global G
-        G = np.zeros((self.n_ele + 1, self.numBasisFuncs * self.n_ele))
+        G = np.zeros((self.n_ele + 1, numBasisFuncs * self.n_ele))
         for i in range(self.n_ele):
             for j in range(2):
                 m = map_g_x[j, i]
-                for k in range(self.numBasisFuncs):
+                for k in range(numBasisFuncs):
                     n = map_g_y[k, i]
                     G[m, n] += g[j, k]
         G = G[1:self.n_ele, :]
-
-        return C, E, G, H
+        faceMat = namedtuple('faceMat', ['C', 'E', 'G', 'H'])
+        return faceMat(C, E, G, H)