2
$\begingroup$

I'm currently attempting to turn my code for solving the laplace equation using finite element approximations into an adaptive one using the dual weighted residual as my error estimator: i.e. my indicator for element i: $η_i=(f[i]−∑_0^N A_{ij}U_h[j])∗z[i]$. i.e. I set up f, A and z over a very fine mesh and solve U_h over the normal mesh then interpolate it so i can perform the above operation.

I am solving for the dual solution z relatively simply; i.e. B(v,z) = Q(v) where B is the stiffness matrix and Q is a vector of 0s with 1 at the quantity of interest.

However the errors do not behave how they should; as an example i got the mesh to just increase from 10 in uniform multiples of 10 until 150. However these are the error estimates (as solved above):

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

They seem to drop then grow, drop then grow repeatedly and i have no idea why, does anybody have any clues what could be causing this?

The following is my code

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]

I understand theres alot of code here but i have no idea where else to turn. I'm almost certain my functions for solving the finite element method are correct (i.e. the Stiffness matrix, nodal quadrature and the solver). Any help would be appreciated alot. Thanks, James

EDIT, this is my plot of the error indicator which must be incorrect but i cannot figure out what is wrong with it. Error estimate over mesh

$\endgroup$
  • 1
    $\begingroup$ Did you verify that the solution on a uniform mesh is correct? Did you verify that the dual solution is correct? If they are, did you verify that the error estimator (as a function over the mesh) looks like you expect it to look like? You can't expect to write a larger code and assume that it will work correctly on first try. You need to verify each individual step or your end result will be wrong (not may be wrong, but will be wrong). $\endgroup$ – Wolfgang Bangerth Sep 5 '15 at 20:10
  • $\begingroup$ Hi Wolfgang, Thanks for reading; Both the finite element solution over a uniform mesh and the dual solution are correct. You are correct my error estimator as a function over the mesh is incredibly strange but I have analysed it hundreds of times now. Not quite sure what's wrong with it but i'm sure i'm missing something stupid. I have edited in the plot of my error estimate over the mesh $\endgroup$ – malonej Sep 5 '15 at 20:21
  • $\begingroup$ I have no idea either, obviously, other than saying that this does not look correct. My experience is that you can do almost everything on a coarse mesh with a solution that you have manufactured and consequently know. In your case, I would first verify that both the primal and dual solutions converge for a simple case (say, computing the error in the average, which yields a pleasant dual problem), and then make sure that your error estimator converges under global mesh refinement. $\endgroup$ – Wolfgang Bangerth Sep 8 '15 at 13:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.