# Odd accuracy barrier in C/PETSc regarding finite elements

I’m implementing a finite element code (translating from a working MATLAB version, so I have results to compare to) and for some odd reason, some of my computations are only accurate to around 6 decimal places.

I’m working with doubles and PetscScalars and I’ve double and quadruple checked my code and am fairly certain there is nothing wrong with my implementation. I went back and made sure all of my multiplications are with doubles or PetscScalars, though I’m not sure if that changes anything.

I test on an exact solution where I should get machine precision error. On a refinement study, the coarsest mesh (two-element triangular mesh on $[0,1]^{2}$) I get machine precision error but for every refinement after that, the error bumps up to around $10^{-7}$ and completely flatlines there. For exact solutions that can’t be exactly expressed in terms of the function space, I’m obtaining proper convergence rates for basis degree 1 and 2, and then the rates rapidly drop off for any high degrees.

I’ve been trying to hunt down the error for a while now. This is my first time programming such a large project in C so I’m sure there are some nuances that I just don’t know about, although I think I’m experience enough with C to have not made any errors when translating the working code. Do I have to specify that the PetscScalar’s are doubles? I reconfigured --with-precision=double but it didn’t change the results.

• I believe you earlier asked a question about using files for input/output with PETSc. Are you writing data out to a file and then passing that data into the solver? If so, what format are you using for writing out the coefficients of your system of equations, and does it give the required precision? – Brian Borchers Nov 23 '13 at 14:01
• Ha, good point. We've all made the mistake to write solutions to a file with standard C precision of 6 digits and then wonder why we can't get more accurate :-) – Wolfgang Bangerth Nov 23 '13 at 17:11
• @BrianBorchers I’m actually not doing that now. But I AM reading my meshes from files though since there were a few I had from MATLAB that I wanted to use. In particular, the coordinates of the nodes and the normal vectors to the edges of the mesh are read in from a .dat file (with default precision of 5 digits, yikes). Will change that now and see – Justin Dong Nov 23 '13 at 17:29
• Wow, I can’t believe I overlooked that. Reading the mesh data in with default precision from MATLAB turned out to be the issue. It wasn’t that big of an issue for the nodes but definitely was for the normals. Resaved the mesh data with 15 digits and that made all the difference. Thanks to everyone for the quick help :) – Justin Dong Nov 23 '13 at 17:36
• We've all done that at one point or other. I've made that mistake more than once. – Wolfgang Bangerth Nov 23 '13 at 20:03

• for now, I am just solving using PCLU. The issues I was mentioning above are with the local matrices, which don’t involve any linear solves yet. – Justin Dong Nov 23 '13 at 17:07
• yes, I have to solve a linear system to compute the solution. The computed solutions are inaccurate as mentioned in my original post, but I was viewing the local matrices and just comparing them to expected results and noting that even the local matrices seem to be suffering from some nontrivial accuracy issues. For example, in some of the local matrices, the computed entries should be 0. In MATLAB, those entries are close to machine precision in value, but in C I’m getting values of 0.000002, which I think is a problem. Sorry if I was unclear before. I will append my original post. – Justin Dong Nov 23 '13 at 17:15