Hellow, I have found two bugs in the codes, and I listed them in the attachment. I have a question about why the backward differentiation formula starts after m_pressureCalls>2, in function void Nektar::Extrapolate::AccelerationBDF(Array< OneD, Array< OneD, NekDouble > > & array). Since when m_pressureCalls=2, we can already get a first order accuracy result. thanks Ankang
Hi Ankang, Thanks for sending the bug reports! I've only looked at bug (1) so far. The interpolation error is quite interesting and I agree not necessarily accurate in the vicinity of the interpolation points in its present form. In your example, the interpolation is quite close indeed and perhaps closer than we would usually do this, which may explain why this has not come up before. I was therefore wondering what lead you to find this inaccuracy, just so I can better understand the context a bit? However the 'traditional implementation' of the Lagrange interpolant as you've written it is also numerically unstable (I believe) for larger numbers of interpolation points. We might consider using a barycentric form, as explained in the following two articles: https://people.maths.ox.ac.uk/trefethen/mythspaper.pdf http://www.maths.manchester.ac.uk/~higham/narep/narep440.pdf Don't know if someone with more knowledge of this might be able to weigh in on this aspect! On the m_pressureCalls issue -- I agree one can get a first-order approximation from the first timestep onwards. Not sure why the original implementation doesn't have this. Dave
On 14 Apr 2017, at 13:06, Ankang Gao <gaoak@pku.edu.cn> wrote:
Hellow,
I have found two bugs in the codes, and I listed them in the attachment.
I have a question about why the backward differentiation formula starts after m_pressureCalls>2, in function void Nektar::Extrapolate::AccelerationBDF(Array< OneD, Array< OneD, NekDouble > > & array). Since when m_pressureCalls=2, we can already get a first order accuracy result.
thanks
Ankang <reportbug_question.pdf>_______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/nektar-users
Hi Dave Thanks for your fast reply. I encountered the interpolation error problem when I used function PhysEvaluate(const Array< OneD, const NekDouble > & coords, const Array< OneD, const NekDouble > & physvals). I found if I set coords as the integral points(these were obtained by using function GetCoords) in phyEvaluate, sometimes I got a very poor accuracy. Further debug led me to find this inaccuracy in Lagrange interpolation functions.
-----原始邮件----- 发件人: "David Moxey" <d.moxey@imperial.ac.uk> 发送时间: 2017-04-15 05:47:20 (星期六) 收件人: "Ankang Gao" <gaoak@pku.edu.cn> 抄送: nektar-users <nektar-users@imperial.ac.uk> 主题: Re: [Nektar-users] 2 bugs and 1 question
Hi Ankang,
Thanks for sending the bug reports! I've only looked at bug (1) so far.
The interpolation error is quite interesting and I agree not necessarily accurate in the vicinity of the interpolation points in its present form. In your example, the interpolation is quite close indeed and perhaps closer than we would usually do this, which may explain why this has not come up before. I was therefore wondering what lead you to find this inaccuracy, just so I can better understand the context a bit?
However the 'traditional implementation' of the Lagrange interpolant as you've written it is also numerically unstable (I believe) for larger numbers of interpolation points. We might consider using a barycentric form, as explained in the following two articles:
https://people.maths.ox.ac.uk/trefethen/mythspaper.pdf http://www.maths.manchester.ac.uk/~higham/narep/narep440.pdf
Don't know if someone with more knowledge of this might be able to weigh in on this aspect!
On the m_pressureCalls issue -- I agree one can get a first-order approximation from the first timestep onwards. Not sure why the original implementation doesn't have this.
Dave
On 14 Apr 2017, at 13:06, Ankang Gao <gaoak@pku.edu.cn> wrote:
Hellow,
I have found two bugs in the codes, and I listed them in the attachment.
I have a question about why the backward differentiation formula starts after m_pressureCalls>2, in function void Nektar::Extrapolate::AccelerationBDF(Array< OneD, Array< OneD, NekDouble > > & array). Since when m_pressureCalls=2, we can already get a first order accuracy result.
thanks
Ankang <reportbug_question.pdf>_______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/nektar-users
HI Ankang, Thanks for your comments which are very insightful. I need to get together and discuss this with Dave to see what we might do about this. Cheers, Spencer. PS I have not had a close look at m_pressureCalls yet! On 15 Apr 2017, at 01:53, Ankang Gao <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> wrote: Hi Dave Thanks for your fast reply. I encountered the interpolation error problem when I used function PhysEvaluate(const Array< OneD, const NekDouble > & coords, const Array< OneD, const NekDouble > & physvals). I found if I set coords as the integral points(these were obtained by using function GetCoords) in phyEvaluate, sometimes I got a very poor accuracy. Further debug led me to find this inaccuracy in Lagrange interpolation functions. -----原始邮件----- 发件人: "David Moxey" <d.moxey@imperial.ac.uk<mailto:d.moxey@imperial.ac.uk>> 发送时间: 2017-04-15 05:47:20 (星期六) 收件人: "Ankang Gao" <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> 抄送: nektar-users <nektar-users@imperial.ac.uk<mailto:nektar-users@imperial.ac.uk>> 主题: Re: [Nektar-users] 2 bugs and 1 question Hi Ankang, Thanks for sending the bug reports! I've only looked at bug (1) so far. The interpolation error is quite interesting and I agree not necessarily accurate in the vicinity of the interpolation points in its present form. In your example, the interpolation is quite close indeed and perhaps closer than we would usually do this, which may explain why this has not come up before. I was therefore wondering what lead you to find this inaccuracy, just so I can better understand the context a bit? However the 'traditional implementation' of the Lagrange interpolant as you've written it is also numerically unstable (I believe) for larger numbers of interpolation points. We might consider using a barycentric form, as explained in the following two articles: https://people.maths.ox.ac.uk/trefethen/mythspaper.pdf http://www.maths.manchester.ac.uk/~higham/narep/narep440.pdf Don't know if someone with more knowledge of this might be able to weigh in on this aspect! On the m_pressureCalls issue -- I agree one can get a first-order approximation from the first timestep onwards. Not sure why the original implementation doesn't have this. Dave On 14 Apr 2017, at 13:06, Ankang Gao <gaoak@pku.edu.cn> wrote: Hellow, I have found two bugs in the codes, and I listed them in the attachment. I have a question about why the backward differentiation formula starts after m_pressureCalls>2, in function void Nektar::Extrapolate::AccelerationBDF(Array< OneD, Array< OneD, NekDouble > > & array). Since when m_pressureCalls=2, we can already get a first order accuracy result. thanks Ankang <reportbug_question.pdf>_______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/nektar-users _______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk<mailto:Nektar-users@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/nektar-users Spencer Sherwin McLaren Racing/Royal Academy of Engineering Research Chair, Professor of Computational Fluid Mechanics, Department of Aeronautics, Imperial College London South Kensington Campus London SW7 2AZ s.sherwin@imperial.ac.uk<mailto:s.sherwin@imperial.ac.uk> +44 (0) 20 759 45052
Hi Ankang, You are right about the second bug. Thanks for spotting it. I have pushed a fix for this bug to the repository, which I think will be in version 4.4.1. Cheers, Douglas ________________________________ From: nektar-users-bounces@imperial.ac.uk <nektar-users-bounces@imperial.ac.uk> on behalf of Sherwin, Spencer J <s.sherwin@imperial.ac.uk> Sent: 16 April 2017 22:21:08 To: Ankang Gao Cc: nektar-users Subject: Re: [Nektar-users] 2 bugs and 1 question HI Ankang, Thanks for your comments which are very insightful. I need to get together and discuss this with Dave to see what we might do about this. Cheers, Spencer. PS I have not had a close look at m_pressureCalls yet! On 15 Apr 2017, at 01:53, Ankang Gao <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> wrote: Hi Dave Thanks for your fast reply. I encountered the interpolation error problem when I used function PhysEvaluate(const Array< OneD, const NekDouble > & coords, const Array< OneD, const NekDouble > & physvals). I found if I set coords as the integral points(these were obtained by using function GetCoords) in phyEvaluate, sometimes I got a very poor accuracy. Further debug led me to find this inaccuracy in Lagrange interpolation functions. -----原始邮件----- 发件人: "David Moxey" <d.moxey@imperial.ac.uk<mailto:d.moxey@imperial.ac.uk>> 发送时间: 2017-04-15 05:47:20 (星期六) 收件人: "Ankang Gao" <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> 抄送: nektar-users <nektar-users@imperial.ac.uk<mailto:nektar-users@imperial.ac.uk>> 主题: Re: [Nektar-users] 2 bugs and 1 question Hi Ankang, Thanks for sending the bug reports! I've only looked at bug (1) so far. The interpolation error is quite interesting and I agree not necessarily accurate in the vicinity of the interpolation points in its present form. In your example, the interpolation is quite close indeed and perhaps closer than we would usually do this, which may explain why this has not come up before. I was therefore wondering what lead you to find this inaccuracy, just so I can better understand the context a bit? However the 'traditional implementation' of the Lagrange interpolant as you've written it is also numerically unstable (I believe) for larger numbers of interpolation points. We might consider using a barycentric form, as explained in the following two articles: https://people.maths.ox.ac.uk/trefethen/mythspaper.pdf http://www.maths.manchester.ac.uk/~higham/narep/narep440.pdf Don't know if someone with more knowledge of this might be able to weigh in on this aspect! On the m_pressureCalls issue -- I agree one can get a first-order approximation from the first timestep onwards. Not sure why the original implementation doesn't have this. Dave On 14 Apr 2017, at 13:06, Ankang Gao <gaoak@pku.edu.cn> wrote: Hellow, I have found two bugs in the codes, and I listed them in the attachment. I have a question about why the backward differentiation formula starts after m_pressureCalls>2, in function void Nektar::Extrapolate::AccelerationBDF(Array< OneD, Array< OneD, NekDouble > > & array). Since when m_pressureCalls=2, we can already get a first order accuracy result. thanks Ankang <reportbug_question.pdf>_______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/nektar-users _______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk<mailto:Nektar-users@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/nektar-users Spencer Sherwin McLaren Racing/Royal Academy of Engineering Research Chair, Professor of Computational Fluid Mechanics, Department of Aeronautics, Imperial College London South Kensington Campus London SW7 2AZ s.sherwin@imperial.ac.uk<mailto:s.sherwin@imperial.ac.uk> +44 (0) 20 759 45052
Hi Ankang, Sorry it has taken rather long to get back to you. Thanks for your clear explanation which I think is very interesting and we should do something about this. Do you have a modified polylib.cpp we can include. I am not sure if the papers that Dave forwarded offer a more robust way of evaluating the lagging polynomial near the quadrature points since that is the only other modification I can think we might consider. I am cc’ing Mike so he can also perhaps comments on this point. Cheers, Spencer.
On 14 Apr 2017, at 13:06, Ankang Gao <gaoak@pku.edu.cn> wrote:
Hellow,
I have found two bugs in the codes, and I listed them in the attachment.
I have a question about why the backward differentiation formula starts after m_pressureCalls>2, in function void Nektar::Extrapolate::AccelerationBDF(Array< OneD, Array< OneD, NekDouble > > & array). Since when m_pressureCalls=2, we can already get a first order accuracy result.
thanks
Ankang <reportbug_question.pdf>_______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/nektar-users
Spencer Sherwin McLaren Racing/Royal Academy of Engineering Research Chair, Professor of Computational Fluid Mechanics, Department of Aeronautics, Imperial College London South Kensington Campus London SW7 2AZ s.sherwin@imperial.ac.uk +44 (0) 20 759 45052
Hi, I use the same Lagrange interpolation polynomial as Nektar::LibUtilities::GaussPoints::LagrangePoly or Nektar::LibUtilities::PolyEPoints::LagrangePoly. The Polylib.cpp file I modified is attached, and the differences between this file and the original one are shown in another attachment 'difference.txt' which is obtained by command "diff" in linux OS. After this modification, the main numerical error of funtion PhysEvaluate is contributed by NewtonIterationForLocCoord, which is 2.e-6 maximum in my tests with 12 points each direction in a quad geometry. I'm glad to see that my report is of use to improve nektar++. In my opinion, Nektar++ is a wonderful code, and it has become the main numerical tool in my research. I am looking forward to nektar++ getting more and more powerful. Ankang -----原始邮件----- 发件人:"Sherwin, Spencer J" <s.sherwin@imperial.ac.uk> 发送时间:2017-04-21 22:51:00 (星期五) 收件人: "Ankang Gao" <gaoak@pku.edu.cn> 抄送: nektar-users <nektar-users@imperial.ac.uk> 主题: Re: [Nektar-users] 2 bugs and 1 question Hi Ankang, Sorry it has taken rather long to get back to you. Thanks for your clear explanation which I think is very interesting and we should do something about this. Do you have a modified polylib.cpp we can include. I am not sure if the papers that Dave forwarded offer a more robust way of evaluating the lagging polynomial near the quadrature points since that is the only other modification I can think we might consider. I am cc’ing Mike so he can also perhaps comments on this point. Cheers, Spencer.
On 14 Apr 2017, at 13:06, Ankang Gao <gaoak@pku.edu.cn> wrote:
Hellow,
I have found two bugs in the codes, and I listed them in the attachment.
I have a question about why the backward differentiation formula starts after m_pressureCalls>2, in function void Nektar::Extrapolate::AccelerationBDF(Array< OneD, Array< OneD, NekDouble > > & array). Since when m_pressureCalls=2, we can already get a first order accuracy result.
thanks
Ankang <reportbug_question.pdf>_______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/nektar-users
Spencer Sherwin McLaren Racing/Royal Academy of Engineering Research Chair, Professor of Computational Fluid Mechanics, Department of Aeronautics, Imperial College London South Kensington Campus London SW7 2AZ s.sherwin@imperial.ac.uk +44 (0) 20 759 45052
Hi All, From Ankang’s report about polylib.cpp, it was said that “Problem happens when z is very close to zj but their distance is still greater than EPS, since at this situation the formula in hj contains subtraction of two very close (double) float numbers, which causes heavy loss of accuracy." It seemed that the double-precision floating-point format was being challenged. Double precision only gives 15-17 significant decimal digits precision. (See wikipedia) Sign bit: 1 bit Exponent: 11 bits Significant precision: 53 bits (52 explicitly stored) -> 15-17 significant decimal digits However, we could have a more professional way to circumvent/reduce "subtractive cancellation". See the example given below and it can be a way to do some modifications. The codes: #include <iostream> #include <limits.h> #include <float.h> #include <iomanip> #include <cmath> using namespace std; typedef std::numeric_limits< double > dbl; double scale=1e16; int main(){ double a = 1.111111112355343; double b = 1.111111112344554; cout<< setprecision(30)<< a << endl << setprecision(30) << b <<endl; double expn = floor(log10(abs(a-b))); // Assuming that the two numbers are very close to each other and both are greater than 1. The difference is approaching to the limitation of significant digits. double a0 = a*pow(10,-expn-1)-floor(a*pow(10,-expn-1)); double b0 = b*pow(10,-expn-1)-floor(b*pow(10,-expn-1)); double a1 = round(a0*pow(10,15+expn+1)); // get rid of rubbish double b1 = round(b0*pow(10,15+expn+1)); // get rid of rubbish cout<< setprecision(30)<< a1 << endl << setprecision(30) << b1 <<endl; cout<<"The difference between two number : "<<endl<< setprecision(20)<< a-b <<endl; cout<<"Reference Value :"<<endl <<" 92686996014.459171378255630735" <<endl; cout<<"Calculations : "<<endl; cout<<"Direct Calc : "<<setprecision(30)<<1/(a-b)<<endl; cout<<"Optimised Calc : "<<setprecision(30)<<1e15/(a1-b1)<<endl; return 0; } The output: 1.11111111235534298913307793555 1.11111111234455406382437558932 55343 44554 The difference between two number : 1.0788925308702346229e-11 Reference Value: 92686996014.459171378255630735 Calculations : Direct Calc : 92687637682.8190765380859375 Optimised Calc : 92686996014.45916748046875 Hope it could be useful. Cheers, Hui On 21 Apr 2017, at 19:47, Ankang Gao <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> wrote: Hi, I use the same Lagrange interpolation polynomial as Nektar::LibUtilities::GaussPoints::LagrangePoly or Nektar::LibUtilities::PolyEPoints::LagrangePoly. The Polylib.cpp file I modified is attached, and the differences between this file and the original one are shown in another attachment 'difference.txt' which is obtained by command "diff" in linux OS. After this modification, the main numerical error of funtion PhysEvaluate is contributed by NewtonIterationForLocCoord, which is 2.e-6 maximum in my tests with 12 points each direction in a quad geometry. I'm glad to see that my report is of use to improve nektar++. In my opinion, Nektar++ is a wonderful code, and it has become the main numerical tool in my research. I am looking forward to nektar++ getting more and more powerful. Ankang -----原始邮件----- 发件人:"Sherwin, Spencer J" <s.sherwin@imperial.ac.uk<mailto:s.sherwin@imperial.ac.uk>> 发送时间:2017-04-21 22:51:00 (星期五) 收件人: "Ankang Gao" <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> 抄送: nektar-users <nektar-users@imperial.ac.uk<mailto:nektar-users@imperial.ac.uk>> 主题: Re: [Nektar-users] 2 bugs and 1 question Hi Ankang, Sorry it has taken rather long to get back to you. Thanks for your clear explanation which I think is very interesting and we should do something about this. Do you have a modified polylib.cpp we can include. I am not sure if the papers that Dave forwarded offer a more robust way of evaluating the lagging polynomial near the quadrature points since that is the only other modification I can think we might consider. I am cc’ing Mike so he can also perhaps comments on this point. Cheers, Spencer.
On 14 Apr 2017, at 13:06, Ankang Gao <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> wrote:
Hellow,
I have found two bugs in the codes, and I listed them in the attachment.
I have a question about why the backward differentiation formula starts after m_pressureCalls>2, in function void Nektar::Extrapolate::AccelerationBDF(Array< OneD, Array< OneD, NekDouble > > & array). Since when m_pressureCalls=2, we can already get a first order accuracy result.
thanks
Ankang <reportbug_question.pdf>_______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk<mailto:Nektar-users@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/nektar-users
Spencer Sherwin McLaren Racing/Royal Academy of Engineering Research Chair, Professor of Computational Fluid Mechanics, Department of Aeronautics, Imperial College London South Kensington Campus London SW7 2AZ s.sherwin@imperial.ac.uk<mailto:s.sherwin@imperial.ac.uk> +44 (0) 20 759 45052 <Polylib.cpp><difference.txt>_______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk<mailto:Nektar-users@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/nektar-users
Hi Hui, Thanks for this code snippet which is very interesting. However I think the issue here is that we are using a threshold where we set the value to 1 when a points z is near to zj within a prescribed tolerance EPS. However depending on this value of EPS we are giving a value of 1 when it should be 0.999999xxx. In the code from Ankang I believe he is swapping to the classical form of the Lagrange polynomial when we are within this tolerance (@Ankang can you confirm I am correct). So although we could swap to the arithmetic you suggest below, will we not be left with the issue that we still need to evaluate the function when z = zj. Will the approach below allow for this case? If not we might prefer to swap to Anyang’s suggested modification. Cheers, Spencer On 22 Apr 2017, at 15:28, Xu, Hui <hui.xu@imperial.ac.uk<mailto:hui.xu@imperial.ac.uk>> wrote: Hi All, From Ankang’s report about polylib.cpp, it was said that “Problem happens when z is very close to zj but their distance is still greater than EPS, since at this situation the formula in hj contains subtraction of two very close (double) float numbers, which causes heavy loss of accuracy." It seemed that the double-precision floating-point format was being challenged. Double precision only gives 15-17 significant decimal digits precision. (See wikipedia) Sign bit: 1 bit Exponent: 11 bits Significant precision: 53 bits (52 explicitly stored) -> 15-17 significant decimal digits However, we could have a more professional way to circumvent/reduce "subtractive cancellation". See the example given below and it can be a way to do some modifications. The codes: #include <iostream> #include <limits.h> #include <float.h> #include <iomanip> #include <cmath> using namespace std; typedef std::numeric_limits< double > dbl; double scale=1e16; int main(){ double a = 1.111111112355343; double b = 1.111111112344554; cout<< setprecision(30)<< a << endl << setprecision(30) << b <<endl; double expn = floor(log10(abs(a-b))); // Assuming that the two numbers are very close to each other and both are greater than 1. The difference is approaching to the limitation of significant digits. double a0 = a*pow(10,-expn-1)-floor(a*pow(10,-expn-1)); double b0 = b*pow(10,-expn-1)-floor(b*pow(10,-expn-1)); double a1 = round(a0*pow(10,15+expn+1)); // get rid of rubbish double b1 = round(b0*pow(10,15+expn+1)); // get rid of rubbish cout<< setprecision(30)<< a1 << endl << setprecision(30) << b1 <<endl; cout<<"The difference between two number : "<<endl<< setprecision(20)<< a-b <<endl; cout<<"Reference Value :"<<endl <<" 92686996014.459171378255630735" <<endl; cout<<"Calculations : "<<endl; cout<<"Direct Calc : "<<setprecision(30)<<1/(a-b)<<endl; cout<<"Optimised Calc : "<<setprecision(30)<<1e15/(a1-b1)<<endl; return 0; } The output: 1.11111111235534298913307793555 1.11111111234455406382437558932 55343 44554 The difference between two number : 1.0788925308702346229e-11 Reference Value: 92686996014.459171378255630735 Calculations : Direct Calc : 92687637682.8190765380859375 Optimised Calc : 92686996014.45916748046875 Hope it could be useful. Cheers, Hui On 21 Apr 2017, at 19:47, Ankang Gao <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> wrote: Hi, I use the same Lagrange interpolation polynomial as Nektar::LibUtilities::GaussPoints::LagrangePoly or Nektar::LibUtilities::PolyEPoints::LagrangePoly. The Polylib.cpp file I modified is attached, and the differences between this file and the original one are shown in another attachment 'difference.txt' which is obtained by command "diff" in linux OS. After this modification, the main numerical error of funtion PhysEvaluate is contributed by NewtonIterationForLocCoord, which is 2.e-6 maximum in my tests with 12 points each direction in a quad geometry. I'm glad to see that my report is of use to improve nektar++. In my opinion, Nektar++ is a wonderful code, and it has become the main numerical tool in my research. I am looking forward to nektar++ getting more and more powerful. Ankang -----原始邮件----- 发件人:"Sherwin, Spencer J" <s.sherwin@imperial.ac.uk<mailto:s.sherwin@imperial.ac.uk>> 发送时间:2017-04-21 22:51:00 (星期五) 收件人: "Ankang Gao" <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> 抄送: nektar-users <nektar-users@imperial.ac.uk<mailto:nektar-users@imperial.ac.uk>> 主题: Re: [Nektar-users] 2 bugs and 1 question Hi Ankang, Sorry it has taken rather long to get back to you. Thanks for your clear explanation which I think is very interesting and we should do something about this. Do you have a modified polylib.cpp we can include. I am not sure if the papers that Dave forwarded offer a more robust way of evaluating the lagging polynomial near the quadrature points since that is the only other modification I can think we might consider. I am cc’ing Mike so he can also perhaps comments on this point. Cheers, Spencer.
On 14 Apr 2017, at 13:06, Ankang Gao <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> wrote:
Hellow,
I have found two bugs in the codes, and I listed them in the attachment.
I have a question about why the backward differentiation formula starts after m_pressureCalls>2, in function void Nektar::Extrapolate::AccelerationBDF(Array< OneD, Array< OneD, NekDouble > > & array). Since when m_pressureCalls=2, we can already get a first order accuracy result.
thanks
Ankang <reportbug_question.pdf>_______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk<mailto:Nektar-users@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/nektar-users
Spencer Sherwin McLaren Racing/Royal Academy of Engineering Research Chair, Professor of Computational Fluid Mechanics, Department of Aeronautics, Imperial College London South Kensington Campus London SW7 2AZ s.sherwin@imperial.ac.uk<mailto:s.sherwin@imperial.ac.uk> +44 (0) 20 759 45052 <Polylib.cpp><difference.txt>_______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk<mailto:Nektar-users@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/nektar-users Spencer Sherwin McLaren Racing/Royal Academy of Engineering Research Chair, Professor of Computational Fluid Mechanics, Department of Aeronautics, Imperial College London South Kensington Campus London SW7 2AZ s.sherwin@imperial.ac.uk<mailto:s.sherwin@imperial.ac.uk> +44 (0) 20 759 45052
Hi, I think it's not exactly the case Spencer said. Firstly, EPS=1.E-12 is small enough, and in addition the polynomial's first derivative almost vanishes there, so the numerical error caused by giving value 1 when fabs(dz)<EPS is very small. I use a much larger threshold value(EPSERR=3.E-3 compared with EPS=2.E-12) when swapping to the classical form of the Lagrange polynomial, since double-precision floating-point format is be challenged in the original evaluation method if fabs(dz) is too small (like 1.E-11) but still larger than EPS. In the classical Lagrange polynomial, the challenging to double-precision format is avoided since no subtraction of to very close float numbers appears there. Hui's code snippet is amazing(at least for me) and very accurate in calculating 1/(a-b) when 'a' and 'b' are very close, but I'm afraid it is too costly since log10 and pow functions are called, and at the same time its performance is not very stable when fabs(a-b) is not small enough, for example when a=1.1, and b=1.+ 1.E-11: direct calculate of 1/(a-b) gives 10.0000000009999912009561739978 which is acceptable, while Hui's code gives -100000000000 which is wrong. best regard Ankang -----原始邮件----- 发件人:"Sherwin, Spencer J" <s.sherwin@imperial.ac.uk> 发送时间:2017-04-24 17:31:20 (星期一) 收件人: "Xu, Hui" <hui.xu@imperial.ac.uk> 抄送: nektar-users <nektar-users@imperial.ac.uk>, "Ankang Gao" <gaoak@pku.edu.cn> 主题: Re: [Nektar-users] 2 bugs and 1 question Hi Hui, Thanks for this code snippet which is very interesting. However I think the issue here is that we are using a threshold where we set the value to 1 when a points z is near to zj within a prescribed tolerance EPS. However depending on this value of EPS we are giving a value of 1 when it should be 0.999999xxx. In the code from Ankang I believe he is swapping to the classical form of the Lagrange polynomial when we are within this tolerance (@Ankang can you confirm I am correct). So although we could swap to the arithmetic you suggest below, will we not be left with the issue that we still need to evaluate the function when z = zj. Will the approach below allow for this case? If not we might prefer to swap to Anyang’s suggested modification. Cheers, Spencer On 22 Apr 2017, at 15:28, Xu, Hui <hui.xu@imperial.ac.uk> wrote: Hi All, From Ankang’s report about polylib.cpp, it was said that “Problem happens when z is very close to zj but their distance is still greater than EPS, since at this situation the formula in hj contains subtraction of two very close (double) float numbers, which causes heavy loss of accuracy." It seemed that the double-precision floating-point format was being challenged. Double precision only gives 15-17 significant decimal digits precision. (See wikipedia) Sign bit: 1 bit Exponent: 11 bits Significant precision: 53 bits (52 explicitly stored) -> 15-17 significant decimal digits However, we could have a more professional way to circumvent/reduce "subtractive cancellation". See the example given below and it can be a way to do some modifications. The codes: #include <iostream> #include <limits.h> #include <float.h> #include <iomanip> #include <cmath> using namespace std; typedef std::numeric_limits< double > dbl; double scale=1e16; int main(){ double a = 1.111111112355343; double b = 1.111111112344554; cout<< setprecision(30)<< a << endl << setprecision(30) << b <<endl; double expn = floor(log10(abs(a-b))); // Assuming that the two numbers are very close to each other and both are greater than 1. The difference is approaching to the limitation of significant digits. double a0 = a*pow(10,-expn-1)-floor(a*pow(10,-expn-1)); double b0 = b*pow(10,-expn-1)-floor(b*pow(10,-expn-1)); double a1 = round(a0*pow(10,15+expn+1)); // get rid of rubbish double b1 = round(b0*pow(10,15+expn+1)); // get rid of rubbish cout<< setprecision(30)<< a1 << endl << setprecision(30) << b1 <<endl; cout<<"The difference between two number : "<<endl<< setprecision(20)<< a-b <<endl; cout<<"Reference Value :"<<endl <<" 92686996014.459171378255630735" <<endl; cout<<"Calculations : "<<endl; cout<<"Direct Calc : "<<setprecision(30)<<1/(a-b)<<endl; cout<<"Optimised Calc : "<<setprecision(30)<<1e15/(a1-b1)<<endl; return 0; } The output: 1.11111111235534298913307793555 1.11111111234455406382437558932 55343 44554 The difference between two number : 1.0788925308702346229e-11 Reference Value: 92686996014.459171378255630735 Calculations : Direct Calc : 92687637682.8190765380859375 Optimised Calc : 92686996014.45916748046875 Hope it could be useful. Cheers, Hui On 21 Apr 2017, at 19:47, Ankang Gao <gaoak@pku.edu.cn> wrote: Hi, I use the same Lagrange interpolation polynomial as Nektar::LibUtilities::GaussPoints::LagrangePoly or Nektar::LibUtilities::PolyEPoints::LagrangePoly. The Polylib.cpp file I modified is attached, and the differences between this file and the original one are shown in another attachment 'difference.txt' which is obtained by command "diff" in linux OS. After this modification, the main numerical error of funtion PhysEvaluate is contributed by NewtonIterationForLocCoord, which is 2.e-6 maximum in my tests with 12 points each direction in a quad geometry. I'm glad to see that my report is of use to improve nektar++. In my opinion, Nektar++ is a wonderful code, and it has become the main numerical tool in my research. I am looking forward to nektar++ getting more and more powerful. Ankang -----原始邮件----- 发件人:"Sherwin, Spencer J" <s.sherwin@imperial.ac.uk> 发送时间:2017-04-21 22:51:00 (星期五) 收件人: "Ankang Gao" <gaoak@pku.edu.cn> 抄送: nektar-users <nektar-users@imperial.ac.uk> 主题: Re: [Nektar-users] 2 bugs and 1 question Hi Ankang, Sorry it has taken rather long to get back to you. Thanks for your clear explanation which I think is very interesting and we should do something about this. Do you have a modified polylib.cpp we can include. I am not sure if the papers that Dave forwarded offer a more robust way of evaluating the lagging polynomial near the quadrature points since that is the only other modification I can think we might consider. I am cc’ing Mike so he can also perhaps comments on this point. Cheers, Spencer.
On 14 Apr 2017, at 13:06, Ankang Gao <gaoak@pku.edu.cn> wrote:
Hellow,
I have found two bugs in the codes, and I listed them in the attachment.
I have a question about why the backward differentiation formula starts after m_pressureCalls>2, in function void Nektar::Extrapolate::AccelerationBDF(Array< OneD, Array< OneD, NekDouble > > & array). Since when m_pressureCalls=2, we can already get a first order accuracy result.
thanks
Ankang <reportbug_question.pdf>_______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/nektar-users
Spencer Sherwin McLaren Racing/Royal Academy of Engineering Research Chair, Professor of Computational Fluid Mechanics, Department of Aeronautics, Imperial College London South Kensington Campus London SW7 2AZ s.sherwin@imperial.ac.uk +44 (0) 20 759 45052 <Polylib.cpp><difference.txt>_______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/nektar-users Spencer Sherwin McLaren Racing/Royal Academy of Engineering Research Chair, Professor of Computational Fluid Mechanics, Department of Aeronautics, Imperial College London South Kensington Campus London SW7 2AZ s.sherwin@imperial.ac.uk +44 (0) 20 759 45052
Hi, I understood what you meant here about this extreme event. In the case, you do not have enough significant digits to do calculations (the machine precision was still being challenged). From my codes and the digits in the xml, we should know a consistent tolerance should be 1e-7 (or 1e-8). (digit number 15/2) and then you can use what I said to improve. This was what I would like to indicate there. See if it works and gives you a precision about 1e-7. What do you think? Cheers, Hui On 24 Apr 2017, at 13:00, Ankang Gao <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> wrote: Hi, I think it's not exactly the case Spencer said. Firstly, EPS=1.E-12 is small enough, and in addition the polynomial's first derivative almost vanishes there, so the numerical error caused by giving value 1 when fabs(dz)<EPS is very small. I use a much larger threshold value(EPSERR=3.E-3 compared with EPS=2.E-12) when swapping to the classical form of the Lagrange polynomial, since double-precision floating-point format is be challenged in the original evaluation method if fabs(dz) is too small (like 1.E-11) but still larger than EPS. In the classical Lagrange polynomial, the challenging to double-precision format is avoided since no subtraction of to very close float numbers appears there. Hui's code snippet is amazing(at least for me) and very accurate in calculating 1/(a-b) when 'a' and 'b' are very close, but I'm afraid it is too costly since log10 and pow functions are called, and at the same time its performance is not very stable when fabs(a-b) is not small enough, for example when a=1.1, and b=1.+ 1.E-11: direct calculate of 1/(a-b) gives 10.0000000009999912009561739978 which is acceptable, while Hui's code gives -100000000000 which is wrong. best regard Ankang -----原始邮件----- 发件人:"Sherwin, Spencer J" <s.sherwin@imperial.ac.uk<mailto:s.sherwin@imperial.ac.uk>> 发送时间:2017-04-24 17:31:20 (星期一) 收件人: "Xu, Hui" <hui.xu@imperial.ac.uk<mailto:hui.xu@imperial.ac.uk>> 抄送: nektar-users <nektar-users@imperial.ac.uk<mailto:nektar-users@imperial.ac.uk>>, "Ankang Gao" <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> 主题: Re: [Nektar-users] 2 bugs and 1 question Hi Hui, Thanks for this code snippet which is very interesting. However I think the issue here is that we are using a threshold where we set the value to 1 when a points z is near to zj within a prescribed tolerance EPS. However depending on this value of EPS we are giving a value of 1 when it should be 0.999999xxx. In the code from Ankang I believe he is swapping to the classical form of the Lagrange polynomial when we are within this tolerance (@Ankang can you confirm I am correct). So although we could swap to the arithmetic you suggest below, will we not be left with the issue that we still need to evaluate the function when z = zj. Will the approach below allow for this case? If not we might prefer to swap to Anyang’s suggested modification. Cheers, Spencer On 22 Apr 2017, at 15:28, Xu, Hui <hui.xu@imperial.ac.uk<mailto:hui.xu@imperial.ac.uk>> wrote: Hi All, From Ankang’s report about polylib.cpp, it was said that “Problem happens when z is very close to zj but their distance is still greater than EPS, since at this situation the formula in hj contains subtraction of two very close (double) float numbers, which causes heavy loss of accuracy." It seemed that the double-precision floating-point format was being challenged. Double precision only gives 15-17 significant decimal digits precision. (See wikipedia) Sign bit: 1 bit Exponent: 11 bits Significant precision: 53 bits (52 explicitly stored) -> 15-17 significant decimal digits However, we could have a more professional way to circumvent/reduce "subtractive cancellation". See the example given below and it can be a way to do some modifications. The codes: #include <iostream> #include <limits.h> #include <float.h> #include <iomanip> #include <cmath> using namespace std; typedef std::numeric_limits< double > dbl; double scale=1e16; int main(){ double a = 1.111111112355343; double b = 1.111111112344554; cout<< setprecision(30)<< a << endl << setprecision(30) << b <<endl; double expn = floor(log10(abs(a-b))); // Assuming that the two numbers are very close to each other and both are greater than 1. The difference is approaching to the limitation of significant digits. double a0 = a*pow(10,-expn-1)-floor(a*pow(10,-expn-1)); double b0 = b*pow(10,-expn-1)-floor(b*pow(10,-expn-1)); double a1 = round(a0*pow(10,15+expn+1)); // get rid of rubbish double b1 = round(b0*pow(10,15+expn+1)); // get rid of rubbish cout<< setprecision(30)<< a1 << endl << setprecision(30) << b1 <<endl; cout<<"The difference between two number : "<<endl<< setprecision(20)<< a-b <<endl; cout<<"Reference Value :"<<endl <<" 92686996014.459171378255630735" <<endl; cout<<"Calculations : "<<endl; cout<<"Direct Calc : "<<setprecision(30)<<1/(a-b)<<endl; cout<<"Optimised Calc : "<<setprecision(30)<<1e15/(a1-b1)<<endl; return 0; } The output: 1.11111111235534298913307793555 1.11111111234455406382437558932 55343 44554 The difference between two number : 1.0788925308702346229e-11 Reference Value: 92686996014.459171378255630735 Calculations : Direct Calc : 92687637682.8190765380859375 Optimised Calc : 92686996014.45916748046875 Hope it could be useful. Cheers, Hui On 21 Apr 2017, at 19:47, Ankang Gao <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> wrote: Hi, I use the same Lagrange interpolation polynomial as Nektar::LibUtilities::GaussPoints::LagrangePoly or Nektar::LibUtilities::PolyEPoints::LagrangePoly. The Polylib.cpp file I modified is attached, and the differences between this file and the original one are shown in another attachment 'difference.txt' which is obtained by command "diff" in linux OS. After this modification, the main numerical error of funtion PhysEvaluate is contributed by NewtonIterationForLocCoord, which is 2.e-6 maximum in my tests with 12 points each direction in a quad geometry. I'm glad to see that my report is of use to improve nektar++. In my opinion, Nektar++ is a wonderful code, and it has become the main numerical tool in my research. I am looking forward to nektar++ getting more and more powerful. Ankang -----原始邮件----- 发件人:"Sherwin, Spencer J" <s.sherwin@imperial.ac.uk<mailto:s.sherwin@imperial.ac.uk>> 发送时间:2017-04-21 22:51:00 (星期五) 收件人: "Ankang Gao" <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> 抄送: nektar-users <nektar-users@imperial.ac.uk<mailto:nektar-users@imperial.ac.uk>> 主题: Re: [Nektar-users] 2 bugs and 1 question Hi Ankang, Sorry it has taken rather long to get back to you. Thanks for your clear explanation which I think is very interesting and we should do something about this. Do you have a modified polylib.cpp we can include. I am not sure if the papers that Dave forwarded offer a more robust way of evaluating the lagging polynomial near the quadrature points since that is the only other modification I can think we might consider. I am cc’ing Mike so he can also perhaps comments on this point. Cheers, Spencer.
On 14 Apr 2017, at 13:06, Ankang Gao <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> wrote:
Hellow,
I have found two bugs in the codes, and I listed them in the attachment.
I have a question about why the backward differentiation formula starts after m_pressureCalls>2, in function void Nektar::Extrapolate::AccelerationBDF(Array< OneD, Array< OneD, NekDouble > > & array). Since when m_pressureCalls=2, we can already get a first order accuracy result.
thanks
Ankang <reportbug_question.pdf>_______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk<mailto:Nektar-users@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/nektar-users
Spencer Sherwin McLaren Racing/Royal Academy of Engineering Research Chair, Professor of Computational Fluid Mechanics, Department of Aeronautics, Imperial College London South Kensington Campus London SW7 2AZ s.sherwin@imperial.ac.uk<mailto:s.sherwin@imperial.ac.uk> +44 (0) 20 759 45052 <Polylib.cpp><difference.txt>_______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk<mailto:Nektar-users@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/nektar-users Spencer Sherwin McLaren Racing/Royal Academy of Engineering Research Chair, Professor of Computational Fluid Mechanics, Department of Aeronautics, Imperial College London South Kensington Campus London SW7 2AZ s.sherwin@imperial.ac.uk<mailto:s.sherwin@imperial.ac.uk> +44 (0) 20 759 45052
Hi Hui, I have tried the cases from a = 1.+1.E-1, 1.+1.E-2, 1.+1.E-3, ..., to 1.E-10, with fixed b = 1. + 1.E-11, your code gives the same result -1.E11 in all cases which is wrong. I have found where the problem is, that is because expn = floor(log10(abs(a-b))) is a little smaller than it should be, I change it to expn = 1 + floor(log10(abs(a-b))) i.e. plus 1, then high accuracy is restored for all cases. thanks Ankang -----原始邮件----- 发件人:"Xu, Hui" <hui.xu@imperial.ac.uk> 发送时间:2017-04-24 20:46:32 (星期一) 收件人: "Ankang Gao" <gaoak@pku.edu.cn> 抄送: "Sherwin, Spencer J" <s.sherwin@imperial.ac.uk>, nektar-users <nektar-users@imperial.ac.uk> 主题: Re: [Nektar-users] 2 bugs and 1 question Hi, I understood what you meant here about this extreme event. In the case, you do not have enough significant digits to do calculations (the machine precision was still being challenged). From my codes and the digits in the xml, we should know a consistent tolerance should be 1e-7 (or 1e-8). (digit number 15/2) and then you can use what I said to improve. This was what I would like to indicate there. See if it works and gives you a precision about 1e-7. What do you think? Cheers, Hui On 24 Apr 2017, at 13:00, Ankang Gao <gaoak@pku.edu.cn> wrote: Hi, I think it's not exactly the case Spencer said. Firstly, EPS=1.E-12 is small enough, and in addition the polynomial's first derivative almost vanishes there, so the numerical error caused by giving value 1 when fabs(dz)<EPS is very small. I use a much larger threshold value(EPSERR=3.E-3 compared with EPS=2.E-12) when swapping to the classical form of the Lagrange polynomial, since double-precision floating-point format is be challenged in the original evaluation method if fabs(dz) is too small (like 1.E-11) but still larger than EPS. In the classical Lagrange polynomial, the challenging to double-precision format is avoided since no subtraction of to very close float numbers appears there. Hui's code snippet is amazing(at least for me) and very accurate in calculating 1/(a-b) when 'a' and 'b' are very close, but I'm afraid it is too costly since log10 and pow functions are called, and at the same time its performance is not very stable when fabs(a-b) is not small enough, for example when a=1.1, and b=1.+ 1.E-11: direct calculate of 1/(a-b) gives 10.0000000009999912009561739978 which is acceptable, while Hui's code gives -100000000000 which is wrong. best regard Ankang -----原始邮件----- 发件人:"Sherwin, Spencer J" <s.sherwin@imperial.ac.uk> 发送时间:2017-04-24 17:31:20 (星期一) 收件人: "Xu, Hui" <hui.xu@imperial.ac.uk> 抄送: nektar-users <nektar-users@imperial.ac.uk>, "Ankang Gao" <gaoak@pku.edu.cn> 主题: Re: [Nektar-users] 2 bugs and 1 question Hi Hui, Thanks for this code snippet which is very interesting. However I think the issue here is that we are using a threshold where we set the value to 1 when a points z is near to zj within a prescribed tolerance EPS. However depending on this value of EPS we are giving a value of 1 when it should be 0.999999xxx. In the code from Ankang I believe he is swapping to the classical form of the Lagrange polynomial when we are within this tolerance (@Ankang can you confirm I am correct). So although we could swap to the arithmetic you suggest below, will we not be left with the issue that we still need to evaluate the function when z = zj. Will the approach below allow for this case? If not we might prefer to swap to Anyang’s suggested modification. Cheers, Spencer On 22 Apr 2017, at 15:28, Xu, Hui <hui.xu@imperial.ac.uk> wrote: Hi All, From Ankang’s report about polylib.cpp, it was said that “Problem happens when z is very close to zj but their distance is still greater than EPS, since at this situation the formula in hj contains subtraction of two very close (double) float numbers, which causes heavy loss of accuracy." It seemed that the double-precision floating-point format was being challenged. Double precision only gives 15-17 significant decimal digits precision. (See wikipedia) Sign bit: 1 bit Exponent: 11 bits Significant precision: 53 bits (52 explicitly stored) -> 15-17 significant decimal digits However, we could have a more professional way to circumvent/reduce "subtractive cancellation". See the example given below and it can be a way to do some modifications. The codes: #include <iostream> #include <limits.h> #include <float.h> #include <iomanip> #include <cmath> using namespace std; typedef std::numeric_limits< double > dbl; double scale=1e16; int main(){ double a = 1.111111112355343; double b = 1.111111112344554; cout<< setprecision(30)<< a << endl << setprecision(30) << b <<endl; double expn = floor(log10(abs(a-b))); // Assuming that the two numbers are very close to each other and both are greater than 1. The difference is approaching to the limitation of significant digits. double a0 = a*pow(10,-expn-1)-floor(a*pow(10,-expn-1)); double b0 = b*pow(10,-expn-1)-floor(b*pow(10,-expn-1)); double a1 = round(a0*pow(10,15+expn+1)); // get rid of rubbish double b1 = round(b0*pow(10,15+expn+1)); // get rid of rubbish cout<< setprecision(30)<< a1 << endl << setprecision(30) << b1 <<endl; cout<<"The difference between two number : "<<endl<< setprecision(20)<< a-b <<endl; cout<<"Reference Value :"<<endl <<" 92686996014.459171378255630735" <<endl; cout<<"Calculations : "<<endl; cout<<"Direct Calc : "<<setprecision(30)<<1/(a-b)<<endl; cout<<"Optimised Calc : "<<setprecision(30)<<1e15/(a1-b1)<<endl; return 0; } The output: 1.11111111235534298913307793555 1.11111111234455406382437558932 55343 44554 The difference between two number : 1.0788925308702346229e-11 Reference Value: 92686996014.459171378255630735 Calculations : Direct Calc : 92687637682.8190765380859375 Optimised Calc : 92686996014.45916748046875 Hope it could be useful. Cheers, Hui On 21 Apr 2017, at 19:47, Ankang Gao <gaoak@pku.edu.cn> wrote: Hi, I use the same Lagrange interpolation polynomial as Nektar::LibUtilities::GaussPoints::LagrangePoly or Nektar::LibUtilities::PolyEPoints::LagrangePoly. The Polylib.cpp file I modified is attached, and the differences between this file and the original one are shown in another attachment 'difference.txt' which is obtained by command "diff" in linux OS. After this modification, the main numerical error of funtion PhysEvaluate is contributed by NewtonIterationForLocCoord, which is 2.e-6 maximum in my tests with 12 points each direction in a quad geometry. I'm glad to see that my report is of use to improve nektar++. In my opinion, Nektar++ is a wonderful code, and it has become the main numerical tool in my research. I am looking forward to nektar++ getting more and more powerful. Ankang -----原始邮件----- 发件人:"Sherwin, Spencer J" <s.sherwin@imperial.ac.uk> 发送时间:2017-04-21 22:51:00 (星期五) 收件人: "Ankang Gao" <gaoak@pku.edu.cn> 抄送: nektar-users <nektar-users@imperial.ac.uk> 主题: Re: [Nektar-users] 2 bugs and 1 question Hi Ankang, Sorry it has taken rather long to get back to you. Thanks for your clear explanation which I think is very interesting and we should do something about this. Do you have a modified polylib.cpp we can include. I am not sure if the papers that Dave forwarded offer a more robust way of evaluating the lagging polynomial near the quadrature points since that is the only other modification I can think we might consider. I am cc’ing Mike so he can also perhaps comments on this point. Cheers, Spencer.
On 14 Apr 2017, at 13:06, Ankang Gao <gaoak@pku.edu.cn> wrote:
Hellow,
I have found two bugs in the codes, and I listed them in the attachment.
I have a question about why the backward differentiation formula starts after m_pressureCalls>2, in function void Nektar::Extrapolate::AccelerationBDF(Array< OneD, Array< OneD, NekDouble > > & array). Since when m_pressureCalls=2, we can already get a first order accuracy result.
thanks
Ankang <reportbug_question.pdf>_______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/nektar-users
Spencer Sherwin McLaren Racing/Royal Academy of Engineering Research Chair, Professor of Computational Fluid Mechanics, Department of Aeronautics, Imperial College London South Kensington Campus London SW7 2AZ s.sherwin@imperial.ac.uk +44 (0) 20 759 45052 <Polylib.cpp><difference.txt>_______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/nektar-users Spencer Sherwin McLaren Racing/Royal Academy of Engineering Research Chair, Professor of Computational Fluid Mechanics, Department of Aeronautics, Imperial College London South Kensington Campus London SW7 2AZ s.sherwin@imperial.ac.uk +44 (0) 20 759 45052
Hi Ankang, Great !!! Cheers, Hui On 24 Apr 2017, at 14:09, Ankang Gao <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> wrote: Hi Hui, I have tried the cases from a = 1.+1.E-1, 1.+1.E-2, 1.+1.E-3, ..., to 1.E-10, with fixed b = 1. + 1.E-11, your code gives the same result -1.E11 in all cases which is wrong. I have found where the problem is, that is because expn = floor(log10(abs(a-b))) is a little smaller than it should be, I change it to expn = 1 + floor(log10(abs(a-b))) i.e. plus 1, then high accuracy is restored for all cases. thanks Ankang -----原始邮件----- 发件人:"Xu, Hui" <hui.xu@imperial.ac.uk<mailto:hui.xu@imperial.ac.uk>> 发送时间:2017-04-24 20:46:32 (星期一) 收件人: "Ankang Gao" <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> 抄送: "Sherwin, Spencer J" <s.sherwin@imperial.ac.uk<mailto:s.sherwin@imperial.ac.uk>>, nektar-users <nektar-users@imperial.ac.uk<mailto:nektar-users@imperial.ac.uk>> 主题: Re: [Nektar-users] 2 bugs and 1 question Hi, I understood what you meant here about this extreme event. In the case, you do not have enough significant digits to do calculations (the machine precision was still being challenged). From my codes and the digits in the xml, we should know a consistent tolerance should be 1e-7 (or 1e-8). (digit number 15/2) and then you can use what I said to improve. This was what I would like to indicate there. See if it works and gives you a precision about 1e-7. What do you think? Cheers, Hui On 24 Apr 2017, at 13:00, Ankang Gao <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> wrote: Hi, I think it's not exactly the case Spencer said. Firstly, EPS=1.E-12 is small enough, and in addition the polynomial's first derivative almost vanishes there, so the numerical error caused by giving value 1 when fabs(dz)<EPS is very small. I use a much larger threshold value(EPSERR=3.E-3 compared with EPS=2.E-12) when swapping to the classical form of the Lagrange polynomial, since double-precision floating-point format is be challenged in the original evaluation method if fabs(dz) is too small (like 1.E-11) but still larger than EPS. In the classical Lagrange polynomial, the challenging to double-precision format is avoided since no subtraction of to very close float numbers appears there. Hui's code snippet is amazing(at least for me) and very accurate in calculating 1/(a-b) when 'a' and 'b' are very close, but I'm afraid it is too costly since log10 and pow functions are called, and at the same time its performance is not very stable when fabs(a-b) is not small enough, for example when a=1.1, and b=1.+ 1.E-11: direct calculate of 1/(a-b) gives 10.0000000009999912009561739978 which is acceptable, while Hui's code gives -100000000000 which is wrong. best regard Ankang -----原始邮件----- 发件人:"Sherwin, Spencer J" <s.sherwin@imperial.ac.uk<mailto:s.sherwin@imperial.ac.uk>> 发送时间:2017-04-24 17:31:20 (星期一) 收件人: "Xu, Hui" <hui.xu@imperial.ac.uk<mailto:hui.xu@imperial.ac.uk>> 抄送: nektar-users <nektar-users@imperial.ac.uk<mailto:nektar-users@imperial.ac.uk>>, "Ankang Gao" <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> 主题: Re: [Nektar-users] 2 bugs and 1 question Hi Hui, Thanks for this code snippet which is very interesting. However I think the issue here is that we are using a threshold where we set the value to 1 when a points z is near to zj within a prescribed tolerance EPS. However depending on this value of EPS we are giving a value of 1 when it should be 0.999999xxx. In the code from Ankang I believe he is swapping to the classical form of the Lagrange polynomial when we are within this tolerance (@Ankang can you confirm I am correct). So although we could swap to the arithmetic you suggest below, will we not be left with the issue that we still need to evaluate the function when z = zj. Will the approach below allow for this case? If not we might prefer to swap to Anyang’s suggested modification. Cheers, Spencer On 22 Apr 2017, at 15:28, Xu, Hui <hui.xu@imperial.ac.uk<mailto:hui.xu@imperial.ac.uk>> wrote: Hi All, From Ankang’s report about polylib.cpp, it was said that “Problem happens when z is very close to zj but their distance is still greater than EPS, since at this situation the formula in hj contains subtraction of two very close (double) float numbers, which causes heavy loss of accuracy." It seemed that the double-precision floating-point format was being challenged. Double precision only gives 15-17 significant decimal digits precision. (See wikipedia) Sign bit: 1 bit Exponent: 11 bits Significant precision: 53 bits (52 explicitly stored) -> 15-17 significant decimal digits However, we could have a more professional way to circumvent/reduce "subtractive cancellation". See the example given below and it can be a way to do some modifications. The codes: #include <iostream> #include <limits.h> #include <float.h> #include <iomanip> #include <cmath> using namespace std; typedef std::numeric_limits< double > dbl; double scale=1e16; int main(){ double a = 1.111111112355343; double b = 1.111111112344554; cout<< setprecision(30)<< a << endl << setprecision(30) << b <<endl; double expn = floor(log10(abs(a-b))); // Assuming that the two numbers are very close to each other and both are greater than 1. The difference is approaching to the limitation of significant digits. double a0 = a*pow(10,-expn-1)-floor(a*pow(10,-expn-1)); double b0 = b*pow(10,-expn-1)-floor(b*pow(10,-expn-1)); double a1 = round(a0*pow(10,15+expn+1)); // get rid of rubbish double b1 = round(b0*pow(10,15+expn+1)); // get rid of rubbish cout<< setprecision(30)<< a1 << endl << setprecision(30) << b1 <<endl; cout<<"The difference between two number : "<<endl<< setprecision(20)<< a-b <<endl; cout<<"Reference Value :"<<endl <<" 92686996014.459171378255630735" <<endl; cout<<"Calculations : "<<endl; cout<<"Direct Calc : "<<setprecision(30)<<1/(a-b)<<endl; cout<<"Optimised Calc : "<<setprecision(30)<<1e15/(a1-b1)<<endl; return 0; } The output: 1.11111111235534298913307793555 1.11111111234455406382437558932 55343 44554 The difference between two number : 1.0788925308702346229e-11 Reference Value: 92686996014.459171378255630735 Calculations : Direct Calc : 92687637682.8190765380859375 Optimised Calc : 92686996014.45916748046875 Hope it could be useful. Cheers, Hui On 21 Apr 2017, at 19:47, Ankang Gao <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> wrote: Hi, I use the same Lagrange interpolation polynomial as Nektar::LibUtilities::GaussPoints::LagrangePoly or Nektar::LibUtilities::PolyEPoints::LagrangePoly. The Polylib.cpp file I modified is attached, and the differences between this file and the original one are shown in another attachment 'difference.txt' which is obtained by command "diff" in linux OS. After this modification, the main numerical error of funtion PhysEvaluate is contributed by NewtonIterationForLocCoord, which is 2.e-6 maximum in my tests with 12 points each direction in a quad geometry. I'm glad to see that my report is of use to improve nektar++. In my opinion, Nektar++ is a wonderful code, and it has become the main numerical tool in my research. I am looking forward to nektar++ getting more and more powerful. Ankang -----原始邮件----- 发件人:"Sherwin, Spencer J" <s.sherwin@imperial.ac.uk<mailto:s.sherwin@imperial.ac.uk>> 发送时间:2017-04-21 22:51:00 (星期五) 收件人: "Ankang Gao" <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> 抄送: nektar-users <nektar-users@imperial.ac.uk<mailto:nektar-users@imperial.ac.uk>> 主题: Re: [Nektar-users] 2 bugs and 1 question Hi Ankang, Sorry it has taken rather long to get back to you. Thanks for your clear explanation which I think is very interesting and we should do something about this. Do you have a modified polylib.cpp we can include. I am not sure if the papers that Dave forwarded offer a more robust way of evaluating the lagging polynomial near the quadrature points since that is the only other modification I can think we might consider. I am cc’ing Mike so he can also perhaps comments on this point. Cheers, Spencer.
On 14 Apr 2017, at 13:06, Ankang Gao <gaoak@pku.edu.cn<mailto:gaoak@pku.edu.cn>> wrote:
Hellow,
I have found two bugs in the codes, and I listed them in the attachment.
I have a question about why the backward differentiation formula starts after m_pressureCalls>2, in function void Nektar::Extrapolate::AccelerationBDF(Array< OneD, Array< OneD, NekDouble > > & array). Since when m_pressureCalls=2, we can already get a first order accuracy result.
thanks
Ankang <reportbug_question.pdf>_______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk<mailto:Nektar-users@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/nektar-users
Spencer Sherwin McLaren Racing/Royal Academy of Engineering Research Chair, Professor of Computational Fluid Mechanics, Department of Aeronautics, Imperial College London South Kensington Campus London SW7 2AZ s.sherwin@imperial.ac.uk<mailto:s.sherwin@imperial.ac.uk> +44 (0) 20 759 45052 <Polylib.cpp><difference.txt>_______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk<mailto:Nektar-users@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/nektar-users Spencer Sherwin McLaren Racing/Royal Academy of Engineering Research Chair, Professor of Computational Fluid Mechanics, Department of Aeronautics, Imperial College London South Kensington Campus London SW7 2AZ s.sherwin@imperial.ac.uk<mailto:s.sherwin@imperial.ac.uk> +44 (0) 20 759 45052
participants (5)
-
Ankang Gao
-
David Moxey
-
Serson, Douglas
-
Sherwin, Spencer J
-
Xu, Hui