the error is too large

Issue #852 invalid
BeeBreeze created an issue

I have three functions, f1, f2 and f3, satisfied f1 + f3 = f2. I need to get the projections of all three into linear basis fem space. The mesh is defined by UnitCubeMesh(nx,nx,nx), where nx = 8,16,32. I get the L2 projection of f1+f3 and calculate the error between f2 and the projection. The question is that as nx goes from 8 to 32, the error goes from 61 to 49. However, as we all know, the order of error about nx should be 2. Below is my code and it is easy to check the relationship among the three functions with Mathematica or maple.

from dolfin import *
nx_list = [8,16]
for nx in  nx_list:
    mesh = UnitCubeMesh(nx,nx,nx)
    V = FunctionSpace(mesh,'CG',1)  

    A = TrialFunction(V)
    v = TestFunction(V)

    f1 = Expression('-2*(2*cos(x[0])*cos(x[0])*cos(x[0])*cos(2*x[1])*cos(x[2])*cos(x[2])*sin(x[0]) \
        + cos(2*x[1])*sin(2*x[0])*(exp(2*t)*cos(x[1])*cos(x[1])*cos(x[2])*cos(x[2])) \
        - (1/2)*exp(t)*cos(x[2])*cos(x[2])*(-2*pi*sin(x[0])*sin(x[1])+exp(t)*sin(2*x[0])*sin(2*x[1])*(sin(2*x[1]) - sin(2*x[2])))\
        )',t = 1.0,degree = 6)
    f2 = Expression('(-4*cos(x[0])*cos(x[0])*cos(x[0])*cos(2*x[1])*cos(x[2])*cos(x[2])*sin(x[0]) \
        - 2*cos(2*x[1])*sin(2*x[0])*(exp(2*t)*cos(x[1])*cos(x[1])*cos(x[2])*cos(x[2])) \
        + exp(t)*cos(x[2])*cos(x[2])*(-2*pi*sin(x[0])*sin(x[1])+exp(t)*sin(2*x[0])*sin(2*x[1])*(sin(2*x[1]) - sin(2*x[2])))\
        )',t = 1.0,degree = 6)

    s2 = Function(V)

    parameters["form_compiler"]["quadrature_degree"] = 14
    parameters["form_compiler"]["optimize"] = True
    parameters["form_compiler"]["representation"] = 'quadrature'
    A = assemble(A*v*dx)
    b2 = assemble((f2)*v*dx)

    solver = LinearSolver('cg','none')
    solver.parameters['relative_tolerance'] = 1E-16
    solver.parameters['absolute_tolerance'] = 1E-11 
    solver.parameters['maximum_iterations'] = 100000000000
    solver.parameters['nonzero_initial_guess'] = True
    solver.parameters['monitor_convergence'] = False
    solver.solve(A,s2.vector(),b2)

    print errornorm(f1,s2)

Comments (12)

  1. BeeBreeze reporter

    I try to make a simple example. The relationship is that f1 = f2. You can see that the only different is that I convert f1 = -2 * (g1 + g2 + g3 ) into f2 = (-2g1 -2g2 -2*g3). The error is about 2.67 no matter what nx is. If the demo is still too long, I am sorry I have tried my best and let us just ignore it. Thanks for your comments.

  2. Jan Blechta

    Garth is right. This is complex numerical problem (rather than a bug report) which needs disassembling into smaller pieces. Also note that the snippet contains a lot of unrelated code. I don't recommend writing MWEs this way. Moreover is it relevant that CG is used?

  3. BeeBreeze reporter

    My research area is numerical PDE. Based on my knowledge, these two functions are smooth enough, so that CG can work well and P6 is also good enough. The problem is much more complicated actually. In my research, all functions are some kinds of numerical solution and they all come from real life. The function f1 or f2, (sorry I am new here and don't know how to write it like you), is just one term of the exact expression of one numerical solution.

  4. BeeBreeze reporter

    Maybe bigger nx can work better. However, as nx goes from 8 to 16, the error nearly keeps the same. That's very strange.

  5. Jan Blechta

    It's not really important that you're a PhD student. That's not a good excuse for providing a MWE with dozens of boilerplate lines - first one is import numpy and then it continues on and on.

  6. Michal Habera

    Dear @BeeBreeze , you have an integer division in the first expression! Note 1/2 == 0 but 1.0/2 == 0.5. With the proper division your error goes approx. 0.04, 0.01, 0.0026, which is what you, I hope, expect.

    Please reconsider filing a bug next time, as guys recommended you above.

  7. BeeBreeze reporter

    Thank you so much. You are right. In my research, the expression is so long that I may make this kind of mistake. The expression I provide is just copied from it. Thanks a lot.

  8. Log in to comment