5
$\begingroup$

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.

Improve this question
$\endgroup$
5
  • 4
    $\begingroup$ 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? $\endgroup$
    – Brian Borchers
    Nov 23, 2013 at 14:01
  • $\begingroup$ 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 :-) $\endgroup$
    – Wolfgang Bangerth
    Nov 23, 2013 at 17:11
  • $\begingroup$ @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 $\endgroup$
    – Justin Dong
    Nov 23, 2013 at 17:29
  • $\begingroup$ 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 :) $\endgroup$
    – Justin Dong
    Nov 23, 2013 at 17:36
  • $\begingroup$ We've all done that at one point or other. I've made that mistake more than once. $\endgroup$
    – Wolfgang Bangerth
    Nov 23, 2013 at 20:03

2 Answers 2

Reset to default
1
$\begingroup$

Are you solving your linear systems iteratively, using solvers such as CG or GMRES? If so, you need to set the tolerance within these solvers appropriately or your overall accuracy will be limited by the accuracy you get from the linear solver.

Improve this answer
$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – Justin Dong
    Nov 23, 2013 at 17:07
  • $\begingroup$ I don't understand. In order to compute the solution, don't you have to solve a linear system? Or what is it that you're comparing between this code and the analytic solution? $\endgroup$
    – Wolfgang Bangerth
    Nov 23, 2013 at 17:10
  • $\begingroup$ 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. $\endgroup$
    – Justin Dong
    Nov 23, 2013 at 17:15
  • $\begingroup$ Yes, 0.000002 instead of 0 is a good indication that something was done in single precision in the process of setting up these matrices. You seem to be on the right track. $\endgroup$
    – Brian Borchers
    Nov 23, 2013 at 17:37
1
$\begingroup$

For archival purposes, note that the resolution to the problem is found in the comments below the question. Data was being read from a file containing lower-precision data.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.