- changed status to invalid
the error is too large
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)
-
-
reporter - edited description
-
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.
-
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?
-
Why should it be sufficient to interpolate
f1
andf2
to P6? -
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.
-
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.
-
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. -
reporter Well, sorry for that. I will edit my comment.
-
reporter - edited description
-
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.
-
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.
- Log in to comment
Please use fenics-support@fenicsproject.org.
The problem needs to be simplified to be a bug report that can be digested. We can't spend time verifying the very long strings.