Precision error while calculating deformations in the Linear Spring

35 views
Skip to first unread message

HRISHIV NEUPANE

unread,
May 24, 2020, 11:22:22 AM5/24/20
to mastodo...@googlegroups.com, Christopher Wong
Hi,

I was implementing a bilinear translational spring (Ibarra Kwranklier model) and found an issue that may be insignificant for many but made my model work in the wrong way and had to spin my head for few days before catching this issue.
 
I simply added a few lines for testing in the linear spring model and used the displacement controlled problem in the X direction by using a piecewise linear function and FunctionDirichletBC that can be seen in the input file that I have attached. At, time step 1.25 I was expecting the deformation value of 2.5 based on the linear interpolation of the values that I have provided but unfortunately it calculating value in between 2.50000000001 and 2.50000000000001. Since the problem is displacement controlled I am wondering if it is the problem with interpolation in MOOSE (for which I am cc'ing Chris who is looking into MOOSE bugs this summer) or calculation of deformations. 

void
LinearSpring::computeForces()
// spring forces = Kdd * deformations
// spring moments = Krr * rotations
{
//added from here
std::cout<<" \n deformation outside if condition " << _deformations[_qp](0);
if ( _deformations[_qp](0) > ( 2.5 ))
{
std::cout<<" \n deformation inside if condition " << _deformations[_qp](0);
}​
  //added till here
  // forces
  _spring_forces_local(0) = _kx[_qp] * _deformations[_qp](0);
Time Step 25, time = 1.25, dt = 0.05
 
 deformation outside if condition 2
 deformation outside if condition 2 0 Nonlinear |R| = 5.000000e-01
 
 deformation outside if condition 2      0 Linear |R| = 5.000000e-01
      1 Linear |R| = -0.000000e+00
 
 deformation outside if condition 2.5
 deformation inside if condition 2.5

Comment: If the value was actually 2.5 it would never go inside the if condition. 


Thanks
Hrishiv
Graduate Research Assistant
University of Utah


Christopher Wong

unread,
May 24, 2020, 11:24:40 AM5/24/20
to mastodon-users
Ibarra-Medina-Krawlinker* ;)

Chandrakanth Bolisetti

unread,
May 24, 2020, 12:43:33 PM5/24/20
to mastodon-users
Hrishiv,

That is a good question, and I think it belongs to the MOOSE group instead. Have you checked if there are any posts out there on this? I would guess this is related to precision, and that would be venturing into some theoretical aspects of computation. I'll be curious to see what the MOOSE group has to say about this.

Best regards,
Chandu

--
Chandu Bolisetti, Ph.D.
Facility Risk Group Lead
Idaho National Laboratory
716-352-5107 (M)
208-526-8161 (O)


On 5/24/20, 9:26 AM, "mastodo...@googlegroups.com on behalf of Christopher Wong" <mastodo...@googlegroups.com on behalf of crswo...@gmail.com> wrote:

Ibarra-Medina-Krawlinker* ;)

--
https://gcc01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fmooseframework.org%2Fmastodon&amp;data=02%7C01%7Cchandrakanth.bolisetti%40inl.gov%7Cd07eeb596c83440e7bc708d7fff6dfe5%7C4cf464b7869a42368da2a98566485554%7C0%7C0%7C637259308101989084&amp;sdata=ET36Z%2BR%2BMmOGfRrbQfH0oyFA%2FBsijQk5HYwMRDSe8tY%3D&amp;reserved=0
https://gcc01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fidaholab%2Fmastodon&amp;data=02%7C01%7Cchandrakanth.bolisetti%40inl.gov%7Cd07eeb596c83440e7bc708d7fff6dfe5%7C4cf464b7869a42368da2a98566485554%7C0%7C0%7C637259308101989084&amp;sdata=Mj7YcbNiqJBdQS8wo1HOcarMxL2IwX4Yccc0jYGOra8%3D&amp;reserved=0
---
You received this message because you are subscribed to the Google Groups "mastodon-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mastodon-user...@googlegroups.com.
To view this discussion on the web visit https://gcc01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fmsgid%2Fmastodon-users%2F62f63bda-30d5-4da9-b8c8-fca71ecdafb0%2540googlegroups.com&amp;data=02%7C01%7Cchandrakanth.bolisetti%40inl.gov%7Cd07eeb596c83440e7bc708d7fff6dfe5%7C4cf464b7869a42368da2a98566485554%7C0%7C0%7C637259308101989084&amp;sdata=1o%2BzYHp2Vc2hJstaNiDOMf7VbY0qwTYZvpjMvh%2F4xMI%3D&amp;reserved=0.

Swetha Veeraraghavan

unread,
May 24, 2020, 12:45:18 PM5/24/20
to mastodon-users
Hi Hrishiv, 

I had noticed an error in piecewise linear interpolated functions that was causing weird spikes in acceleration. Refer to this thread 

The error in this case was in moose/framework/src/utils/LinearInterpolation.C line 79 due to tolerance issues. 

Changing line 79 to "if (MooseUtils::absoluteFuzzyGreaterEqual(x, _x[i]) && MooseUtils::absoluteFuzzyLessThan(x, _x[i + 1])) instead of if (x >= _x[i] && x < _x[i + 1])" fixes this issue. 

I haven't had a chance to submit a PR to fix this issue yet. But you can try this fix and see this works for you. 

Regards,
Swetha

On Sun, 24 May 2020 at 20:54, Christopher Wong <crswo...@gmail.com> wrote:
Ibarra-Medina-Krawlinker* ;)

--
https://mooseframework.org/mastodon
https://github.com/idaholab/mastodon

---
You received this message because you are subscribed to the Google Groups "mastodon-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mastodon-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mastodon-users/62f63bda-30d5-4da9-b8c8-fca71ecdafb0%40googlegroups.com.

HRISHIV NEUPANE

unread,
May 26, 2020, 12:05:06 AM5/26/20
to mastodon-users
I changed the code as per Swetha's suggestion and compiled MOOSE and then MASTODON and performed the simulation again but the problem that I mentioned earlier still persists.

Also, I have asked the same question in the MOOSE group but have not heard back from them.

In the meantime, I would be happy if somebody could tell me if there is a MathUtils to round off a real number to certain significant digits so that I could bypass the problem for now while this bug (probably) is being looked into by the core team. 

Thanks
Hrishiv 
From: mastodo...@googlegroups.com <mastodo...@googlegroups.com> on behalf of Swetha Veeraraghavan <swetha.vee...@gmail.com>
Sent: Sunday, May 24, 2020 10:44 AM
To: mastodon-users <mastodo...@googlegroups.com>
Subject: Re: [mastodon-users] Precision error while calculating deformations in the Linear Spring
 

Swetha Veeraraghavan

unread,
May 26, 2020, 7:05:15 AM5/26/20
to mastodon-users
Well I couldn't find a MathUtils function that rounds to a certain number of significant digits. There is a function that rounds to an integer though.

But you can write your own as a temporary work around. To round off a real number x upto say 4 significant digits, you would just have to do MathUitls::round(x*10000.0)/10000.0

Regards,
Swetha

To view this discussion on the web visit https://groups.google.com/d/msgid/mastodon-users/BY5PR11MB3910002815FA9567F12851C995B00%40BY5PR11MB3910.namprd11.prod.outlook.com.

Swetha Veeraraghavan

unread,
May 26, 2020, 7:17:06 AM5/26/20
to mastodon-users
Just of out curiosity, how is the change in the 11th significant digit affecting your algorithm? If you are comparing two real numbers x and y and checking if x > y or performing similar operations, you can set tolerances on those operations.

Check MooseUtils::absoluteFuzzyGreaterEqual and similar operations in https://github.com/idaholab/moose/blob/next/framework/include/utils/MooseUtils.h

Regards,
Swetha


HRISHIV NEUPANE

unread,
May 26, 2020, 10:47:01 AM5/26/20
to Swetha Veeraraghavan, mastodon-users
Hi Swetha,

Thanks for the suggestion. I had noticed the Utils that round to an integer value but didn't think that way around. Yes, you are correct I have many comparisons for the displacement of the spring based on which certain flags are turned or off for the calculation of spring forces and this code has been tested and works well in MATLAB.  
 
if ( _isdetpos[_qp] == false)
  {
    if ( _deformations[_qp](0) > ( _th_disp ))
    {
      _isdetpos[_qp] = true;
    }
  }

My threshold displacement was 2.5 and for the displacement of 2.5 _isdetpos[_qp] should've been false but the precision error made it true. 


I have defined these flags (eg _isdamage) as MaterialProperty because these flags are turned on and off based on certain conditions and I don't use them immediately but have to wait for more timesteps for a certain condition to meet to use them.  Can I define them as a variable so that the changed value is inherited to the next time step?


Thanks
Hrishiv
Sent: Tuesday, May 26, 2020 5:16 AM

Christopher Wong

unread,
May 28, 2020, 2:25:49 PM5/28/20
to mastodon-users
Hrishiv, 

First off, `PiecewiseLinear()` is going to invoke `PiecewiseLinearBase::value()`, which invokes `LinearInterpolationTempl<T>::sample()`, and not the `sampleDerivative()` function as Swetha had mistakenly suggested. So you should do some round offs in the way Swetha suggested, but on line 61. And actually, this is not your case, because the `if` statement on line 61 is satisfied at time = 1.25. Your precision error actually arises on the return value on line 62, so you should be rounding that number instead (or, in addition to, to be generic).

Second, let's recognize a deeper question here: why are these numerical errors even arising in the first place? Here are the values produced in `LinearInterpolationTempl<T>::sample()` on the t = 1.25 step:

y[i] = 0, y[i+1] = 3.5, x[i] = 1, x[i+1] = 1.350000000000000089, x = 1.250000000000000444, y = 2.500000000000003997

And your absissca/ordinate pairs are (1, 0) and (1.35, 3.5). Thus, x[i+1] and x (which should be x = 1.25) are slightly wrong here. Even if the math were exact, the error would still arise, because the inputs aren't. However, even the math is not exact, let's say we wrote the interpolation a little differently:

y = y[i+1] + (x - x[i+1]) * (y[i+1] - y[i]) / (x[i+1] - x[i]) = 2.500000000000003553

This is just a tad closer than what you're getting right now, even though it's an equivalent formula.

If you want PiecewiseLinear to pass the exact values that you expect to your spring code, you need to go into LinearInterpolation.C, or any other objects that the value passes through, and consider wether to do a round-off like Swetha suggested or to find the source of the numerical error. For the latter option, you'll find print_trace() a useful tool. However, I suggest you do what you need to do right in your spring model. Make a round-off to the highest precision of the user's machine by using something like `std::numeric_limits<long double>::digits10 + 1`, and then check your yield condition with that. 

For Civil Engineering problems, a value of 4e-14 has absolutely no meaning, even if your units were micrometers. Even the IMK model cannot claim such high accuracy. Still, it would be nice if, at least, the mathematics were exact, even though the theory is not, so its up to you to decide what's most important for this model. 

On Tuesday, May 26, 2020 at 8:47:01 AM UTC-6, HRISHIV NEUPANE wrote:
Hi Swetha,

Thanks for the suggestion. I had noticed the Utils that round to an integer value but didn't think that way around. Yes, you are correct I have many comparisons for the displacement of the spring based on which certain flags are turned or off for the calculation of spring forces and this code has been tested and works well in MATLAB.  
 
if ( _isdetpos[_qp] == false)
  {
    if ( _deformations[_qp](0) > ( _th_disp ))
    {
      _isdetpos[_qp] = true;
    }
  }

My threshold displacement was 2.5 and for the displacement of 2.5 _isdetpos[_qp] should've been false but the precision error made it true. 


I have defined these flags (eg _isdamage) as MaterialProperty because these flags are turned on and off based on certain conditions and I don't use them immediately but have to wait for more timesteps for a certain condition to meet to use them.  Can I define them as a variable so that the changed value is inherited to the next time step?


Thanks
Hrishiv

Sent: Tuesday, May 26, 2020 5:16 AM
To: mastodon-users <mastodo...@googlegroups.com>
Subject: Re: [mastodon-users] Precision error while calculating deformations in the Linear Spring
Just of out curiosity, how is the change in the 11th significant digit affecting your algorithm? If you are comparing two real numbers x and y and checking if x > y or performing similar operations, you can set tolerances on those operations.

Check MooseUtils::absoluteFuzzyGreaterEqual and similar operations in https://github.com/idaholab/moose/blob/next/framework/include/utils/MooseUtils.h

Regards,
Swetha



On Tue, 26 May 2020 at 16:34, Swetha Veeraraghavan <swetha.ve...@gmail.com> wrote:
Well I couldn't find a MathUtils function that rounds to a certain number of significant digits. There is a function that rounds to an integer though.

But you can write your own as a temporary work around. To round off a real number x upto say 4 significant digits, you would just have to do MathUitls::round(x*10000.0)/10000.0

Regards,
Swetha

On Tue, 26 May 2020 at 09:35, HRISHIV NEUPANE <Hrishiv...@utah.edu> wrote:
I changed the code as per Swetha's suggestion and compiled MOOSE and then MASTODON and performed the simulation again but the problem that I mentioned earlier still persists.

Also, I have asked the same question in the MOOSE group but have not heard back from them.

In the meantime, I would be happy if somebody could tell me if there is a MathUtils to round off a real number to certain significant digits so that I could bypass the problem for now while this bug (probably) is being looked into by the core team. 

Thanks
Hrishiv 

Sent: Sunday, May 24, 2020 10:44 AM
To: mastodon-users <mastodo...@googlegroups.com>
Subject: Re: [mastodon-users] Precision error while calculating deformations in the Linear Spring
Hi Hrishiv, 

I had noticed an error in piecewise linear interpolated functions that was causing weird spikes in acceleration. Refer to this thread 

The error in this case was in moose/framework/src/utils/LinearInterpolation.C line 79 due to tolerance issues. 

Changing line 79 to "if (MooseUtils::absoluteFuzzyGreaterEqual(x, _x[i]) && MooseUtils::absoluteFuzzyLessThan(x, _x[i + 1])) instead of if (x >= _x[i] && x < _x[i + 1])" fixes this issue. 

I haven't had a chance to submit a PR to fix this issue yet. But you can try this fix and see this works for you. 

Regards,
Swetha

On Sun, 24 May 2020 at 20:54, Christopher Wong <crswo...@gmail.com> wrote:
Ibarra-Medina-Krawlinker* ;)

--
https://mooseframework.org/mastodon
https://github.com/idaholab/mastodon
---
You received this message because you are subscribed to the Google Groups "mastodon-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mastodo...@googlegroups.com.

--
https://mooseframework.org/mastodon
https://github.com/idaholab/mastodon
---
You received this message because you are subscribed to the Google Groups "mastodon-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mastodo...@googlegroups.com.

--
https://mooseframework.org/mastodon
https://github.com/idaholab/mastodon
---
You received this message because you are subscribed to the Google Groups "mastodon-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mastodo...@googlegroups.com.

--
https://mooseframework.org/mastodon
https://github.com/idaholab/mastodon
---
You received this message because you are subscribed to the Google Groups "mastodon-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mastodo...@googlegroups.com.

HRISHIV NEUPANE

unread,
May 28, 2020, 5:22:54 PM5/28/20
to Christopher Wong, mastodon-users
Hi Chris,

Thanks. I used Swetha's suggestion to round off the deformation values and it worked well for my first case. 
_th_disp_neg[_qp] = MathUtils::round((std::abs(_deformations_old[_qp](0))*10000000)) /10000000;


I need to test some more cases of random vibrations next week before I am done with this model. I feel that the rounding off would only be required for my test case when I use linear interpolation function, as in real simulation the displacement would be calculated based on the accelerations applied to the basket.

Thanks
Hrishiv 


From: mastodo...@googlegroups.com <mastodo...@googlegroups.com> on behalf of Christopher Wong <crswo...@gmail.com>
Sent: Thursday, May 28, 2020 12:25 PM
To unsubscribe from this group and stop receiving emails from it, send an email to mastodon-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mastodon-users/eb2ecca0-4eaf-4c90-a882-e901ce943136%40googlegroups.com.

Christopher Wong

unread,
May 28, 2020, 6:04:28 PM5/28/20
to mastodon-users
Yes, we shouldn't be doing any rounding in the source code as this might result in over-control fluctuations, since the cause is pretty much just random noise. You should look into how you can create your tests such that rounding only occurs for those special cases, i.e, so that the csvdiff or exodiff can simply be close matches, and not exact ones, or something like that. Maybe ask the MOOSE groups if theres a way to do this. 

You could also simply remaunfacture your solution such that this precision error doesn't become a problem. Maybe avoid having the displacement cross the threshold at exactly a given time-step - just offset it a bit.

Swetha Veeraraghavan

unread,
May 30, 2020, 2:23:49 AM5/30/20
to mastodon-users
Hi Hrishiv, 

Sorry for the delayed response. 

When you are performing comparisons among two real numbers, it is advisable to use the fuzzy comparison functions to prevent errors due to precision and any changes in results from one computer architecture to another.     

In your case, if ( _deformations[_qp](0) > ( _th_disp )) would become if (MooseUtils::absoluteFuzzyGreaterThan(_deformations[_qp](0), _th_disp). The third argument to MooseUtils::absoluteFuzzyGreaterThan is the tolerance. If you don't provide any third argument, the default tolerance set in Libmesh is used, which is I think 1e-6 or 1e-8. 

Regarding your question on whether to use the flags as material properties or variables:

In my understanding, material properties are used when the parameter changes with qp, and also when you want to have a stable copy of what happened in the previous time step (like stress old) that is not affected by the changes during the multiple iterations happening in the current time step. This way you can do something like stress[_qp] = stress_old[_qp] + stress_increment[_qp]. Also, when one time step (n) is completed and you are going to n+1, then the material properties at time step n are automatically copied into the "old" material properties. So material properties are designed to only store the values at the current time step (which changes with iterations), previous time step ("old" material property - you cannot change this value) and 2 steps before the current time step ("older" material property - you cannot change this either). 

Assuming that the boolean flag for the linear spring changes only with each spring element and not with the qps within each element, you can use variables for your purpose. When the simulation starts, each element gets a copy of the boolean flag, which you would have to declare as a variable in the include file as say "bool _isdetpos". This variable is local to each spring element, so it can be edited only when actions/calculations are being performed on that spring element. So at time step 5, you can set to true based on some if condition, and again later if the condition is not satisfied at time step 10, it can be set to false then. However, the value of this flag may keep changing during the linear/nonlinear iterations as well. So, you need to decide if that is something acceptable to you.  

Hope this helps,
Swetha

 

To unsubscribe from this group and stop receiving emails from it, send an email to mastodon-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mastodon-users/f37d7c2c-04cf-45cb-84c3-012d924aafc5%40googlegroups.com.
Reply all
Reply to author
Forward
0 new messages