计算科学 - 自适应有限元法 - 拉普拉斯 - 吾爱随笔录

我目前正在尝试将使用有限元近似求解拉普拉斯方程的代码转换为使用双重加权残差作为误差估计器的自适应代码：即元素 i 的指标：。即我在一个非常精细的网格上设置 f、A 和 z，并在正常网格上求解 U_h，然后对其进行插值，这样我就可以执行上述操作。 $η_i=(f[i]−∑_0^N A_{ij}U_h[j])∗z[i]$

我正在相对简单地求解对偶解 z；即 B(v,z) = Q(v) 其中 B 是刚度矩阵，Q 是 0 的向量，其中 1 为感兴趣的量。

但是，错误的行为不应该如此；作为一个例子，我让网格从 10 以 10 的均匀倍数增加到 150。但是这些是误差估计（如上解决）：

0.0335462065038
0.0335462065038
0.0335462065038
0.0335462065038
0.0335462065038
0.0333608591916
0.0335462065038
0.0334614504804
0.0333158830798
0.0335462065038
0.0334317537862
0.0333344104201
0.0335441544821
0.0334327061457

它们似乎先下降然后增长，下降然后重复增长，我不知道为什么，有人有任何线索可能导致这种情况吗？

以下是我的代码

def rhs(x0):
    rhsfunction = np.empty(len(x0))
    for i in range(len(x0)):
        rhsfunction[i] = math.sin(math.pi*x0[i])
    return rhsfunction

def Poisson_Stiffness(x0):
    """Finds the Poisson equation stiffness matrix with any non uniform mesh x0"""

    x0 = np.array(x0)
    N = len(x0) - 1 # The amount of elements; x0, x1, ..., xN

    h = x0[1:] - x0[:-1]

    a = np.zeros(N+1)
    a[0] = 1 #BOUNDARY CONDITIONS
    a[1:-1] = 1/h[1:] + 1/h[:-1]
    a[-1] = 1/h[-1]
    a[N] = 1 #BOUNDARY CONDITIONS

    b = -1/h
    b[0] = 0 #BOUNDARY CONDITIONS

    c = -1/h
    c[N-1] = 0 #BOUNDARY CONDITIONS: DIRICHLET

    data = [a.tolist(), b.tolist(), c.tolist()]
    Positions = [0, 1, -1]
    Stiffness_Matrix = diags(data, Positions, (N+1,N+1))

return Stiffness_Matrix

def NodalQuadrature(x0):
    """Finds the Nodal Quadrature Approximation of sin(pi x)"""

    x0 = np.array(x0)
    h = x0[1:] - x0[:-1]
    N = len(x0) - 1

    approx = np.zeros(len(x0))
    approx[0] = 0 #BOUNDARY CONDITIONS

    for i in range(1,N):
        approx[i] = math.sin(math.pi*x0[i])
        approx[i] = (approx[i]*h[i-1] + approx[i]*h[i])/2

    approx[N] = 0 #BOUNDARY CONDITIONS

return approx

def Solver(x0):

    Stiff_Matrix = Poisson_Stiffness(x0)

    NodalApproximation = NodalQuadrature(x0)
    NodalApproximation[0] = 0

    U = scipy.sparse.linalg.spsolve(Stiff_Matrix, NodalApproximation)

return U

def Dualsolution(rich_mesh,qoi): #BOUNDARY CONDITIONS?
    """Find Z from stiffness matrix Z = K^-1 Q over richer mesh"""

    K = Poisson_Stiffness(rich_mesh)
    Q = np.zeros(len(rich_mesh))
    if qoi in rich_mesh:

        qoi_rich_node = rich_mesh.tolist().index(qoi)
        Q[qoi_rich_node] = 1.0
    else:

        nearest = find_nearest(rich_mesh,qoi)
        if nearest < qoi:

            qoi_lower_node = rich_mesh.tolist().index(nearest)
            qoi_upper_node = qoi_lower_node+1

        else:

            qoi_upper_node = rich_mesh.tolist().index(nearest)
            qoi_lower_node = qoi_upper_node-1

        Q[qoi_upper_node] = (qoi - rich_mesh[qoi_lower_node])/(rich_mesh[qoi_upper_node]-rich_mesh[qoi_lower_node])
        Q[qoi_lower_node] = (rich_mesh[qoi_upper_node] - qoi)/(rich_mesh[qoi_upper_node]-rich_mesh[qoi_lower_node])
    Z = scipy.sparse.linalg.spsolve(K,Q)
    return Z

def Indicators(richx0,basex0,qoi):
"""interpolate U, eta_i = r_i z_i = (f_i - sum(AijUj)Zi"""
    U = Solver(basex0)
    U_inter = interp1d(basex0,U)
    U_rich = U_inter(richx0)
    f = rhs(richx0)
    A = Poisson_Stiffness(richx0).tocsr()
    Z = Dualsolution(richx0,qoi)

    eta = np.empty(len(richx0))
    for i in range(len(eta)):
        eta[i] = f[i]
        for j in range(len(richx0)):
            eta[i] = eta[i] - A[i,j]*U_rich[j]
        eta[i] = eta[i] * Z[i]

我了解这里有很多代码，但我不知道还能去哪里。我几乎可以肯定我求解有限元方法的函数是正确的（即刚度矩阵、节点求积和求解器）。任何帮助将不胜感激。谢谢，詹姆斯

编辑，这是我的错误指示器图，它一定是不正确的，但我无法弄清楚它有什么问题。