Fixed indexing error in PID controller #4855

erikhide · 2023-09-08T13:12:42Z

The PID controller formula (see articles referred to in the comments in the OPM code, e.g. "D. Kuzmin and S. Turek, Numerical simulation of turbulent bubbly flows (2004)") as implemented in the current version of OPM contains an indexing error.

Comparing the implementation and PID controller formula, errors_[0] correspond to $e_{n-2}$, errors_[1] to $e_{n-1}$, and errors_[2] to $e_n$.

bska · 2023-09-08T13:16:45Z

jenkins build this please

bska · 2023-09-08T13:39:38Z

Hm, if this is indeed wrong–I don't have the reference on hand for easy verification–then it's an error that's been with us for a very long time. As far as I can tell this code was first introduced in OPM/opm-core#662 and has remained more or less unchanged, though moved and renamed a number of times, since then.

erikhide · 2023-09-08T13:55:07Z

After looking at it quite a few times, I'm more or less certain that the indexing was wrong in the original implementation. And yes, the PID controller in OPM has remained unchanged since at least 2015. The time stepping algorithm seems to have been adapted to the current PID controller implementation, so there might be a need to update other parts and/or default parameter values to make it work with the new implementation.

That said, there is nothing inherently wrong with using the existing version, but as it stands now it is not a correct implementation of the PID controller as derived in the literature.

totto82 · 2023-09-08T14:04:31Z

Note that the PID controller has in my opinion never worked as intended. It was therefore removed as the default time-stepper long time ago. Thanks for looking into this. If you can make it to work, it would be great.

alfbr · 2023-09-13T12:57:37Z

benchmark please

atgeirr · 2023-09-13T13:32:14Z

Note that the PID controller has in my opinion never worked as intended. It was therefore removed as the default time-stepper long time ago. Thanks for looking into this. If you can make it to work, it would be great.

Our current default is "pid+newtoniteration" which still includes the erroneous pid implementation. I am very much in favour of merging this to correct it, but we must deal with the fallout including failed tests etc. and potentially also change other things as suggested by @erikhide:

so there might be a need to update other parts and/or default parameter values to make it work with the new implementation.

Just a hunch: I cannot remember to have ever seen the verbose messages from the PID function. It might of course be activated in some places, but is it possible that it never was, and only after the current change is it relevant? If so, we could try to make the default TimeStepControl just "newtoniteration" and see if that is easier to deal with.

ytelses · 2023-09-14T02:58:53Z

Benchmark result overview:

Test	Configuration	Relative
opm-git	OPM Benchmark: drogon - Threads: 1	0.998
opm-git	OPM Benchmark: drogon - Threads: 8	1.169
opm-git	OPM Benchmark: punqs3 - Threads: 1	0.979
opm-git	OPM Benchmark: punqs3 - Threads: 8	1.026
opm-git	OPM Benchmark: smeaheia - Threads: 1	1.009
opm-git	OPM Benchmark: smeaheia - Threads: 8	1.123
opm-git	OPM Benchmark: spe10_model_1 - Threads: 1	1.059
opm-git	OPM Benchmark: spe10_model_1 - Threads: 8	1
opm-git	OPM Benchmark: flow_mpi_extra - Threads: 1	1.019
opm-git	OPM Benchmark: flow_mpi_extra - Threads: 8	1.062
opm-git	OPM Benchmark: flow_mpi_norne - Threads: 1	1.026
opm-git	OPM Benchmark: flow_mpi_norne - Threads: 8	1.012
opm-git	OPM Benchmark: flow_mpi_norne_4c_msw - Threads: 1	1.021
opm-git	OPM Benchmark: flow_mpi_norne_4c_msw - Threads: 8	1.007

Speed-up = Total time master / Total time pull request. Above 1.0 is an improvement. *

View result details @ https://www.ytelses.com/opm/?page=result&id=2228

totto82 · 2023-09-14T07:25:11Z

The benchmark results looks very promising.

Our current default is "pid+newtoniteration" which still includes the erroneous pid implementation

My memory is getting worse. I remember looking into this 5 years ago and concluded that PID didn't work as expected but I never noticed the bug. I agree with @atgeirr that this should be merged and with @erikhide that we need to test this and possibly also tune the timestepping differently. Have you done any testing on only using "pid". The "newton" part have some issues at the moment. For some cases it can get stuck on very small timesteps due to some unlucky combinations of NUPCOL and group/well switching.

erikhide · 2023-09-14T10:54:11Z

I have tried to use only "pid" as input to --time-step-control, but so far I cannot say anything conclusive about what works best. My understanding is that the "target iteration approach" was added to ensure that the time steps remained low, because the PID controller tended to favour larger time steps than what might be preferable. The time stepping algorithm as a whole (as currently implemented) finds several different time step proposals, and always chooses the smallest one. Hence, it might become rather conservative in some cases and possibly also fail to increase the time steps when it is both possible and desirable.

I'm also working on some other alternatives (other heuristics than the PID controller, as well as other versions of the PID controller itself), and I'm trying to decide how best to compare different time stepping algorithms to ensure both fast computations and robust behaviour. This might take some time, but hopefully possible improvements will come sometime in the near future.

bska · 2023-09-14T11:03:22Z

I'm also working on some other alternatives (other heuristics than the PID controller, as well as other versions of the PID controller itself), and I'm trying to decide how best to compare different time stepping algorithms to ensure both fast computations and robust behaviour. This might take some time, but hopefully possible improvements will come sometime in the near future.

This is very exciting news! Thank you so much for attacking this area 👍

bska · 2023-10-16T11:13:31Z

opm/simulators/timestepping/TimeStepControl.cpp

@@ -206,7 +206,7 @@ namespace Opm
            const double kD = 0.01 ;
            const double newDt = (dt * std::pow( errors_[ 1 ] / errors_[ 2 ], kP ) *
                                 std::pow( tol_         / errors_[ 2 ], kI ) *
-                                 std::pow( errors_[0]*errors_[0]/errors_[ 1 ]/errors_[ 2 ], kD ));
+                                 std::pow( errors_[1]*errors_[1]/errors_[ 0 ]/errors_[ 2 ], kD ));


I've looked a little bit closer at this. Suppose errors_[0] == 0, which it theoretically could be since the conditional above only checks errors_[1] and errors_[2]. Then newDt will be non-finite, typically inf. Is that reasonable here?

I guess the errors_[1] and errors[2] were included in the conditional because they appeared in the denominators of the PID formula, and hence it was only relevant to check those. However, having changed indexing, errors_[0] (and not errors_[1]) appears in the denominator. It might therefore be reasonable to add errors[0] to (and possibly just remove errors_[1] from) the conditional.

It might therefore be reasonable to add errors[0] to [...] the conditional.

Yes, that's probably needed.

and possibly just remove errors_[1] from the conditional

I'd be wary of doing that. If errors_[1] = 0 we'll get newDt = 0 as well and that seems unexpected to me. Maybe I'm just misunderstanding something.

You're probably right. It appears that that might have been a problem earlier, though. Regardless, it certainly shouldn't hurt to add errors_[0] == 0 to the conditional.

I added it now.

It appears that that (errors_[0] = 0 giving newDt = 0) might have been a problem earlier, though

Agreed. I only remarked on this possibility because we're in the area now so it's an opportunity to fix it. For what it's worth, I tried this PR on a field case and ran into that exact problem. Otherwise I might not have noticed it.

bska · 2023-10-17T13:31:19Z

jenkins build this please

bska

Thank you for the update. I think this is starting to look good so I'm running a build check to assess the impact of the change.

bska · 2023-10-17T13:34:55Z

opm/simulators/timestepping/TimeStepControl.cpp

@@ -206,7 +206,7 @@ namespace Opm
            const double kD = 0.01 ;
            const double newDt = (dt * std::pow( errors_[ 1 ] / errors_[ 2 ], kP ) *
                                 std::pow( tol_         / errors_[ 2 ], kI ) *
-                                 std::pow( errors_[0]*errors_[0]/errors_[ 1 ]/errors_[ 2 ], kD ));
+                                 std::pow( errors_[1]*errors_[1]/errors_[ 0 ]/errors_[ 2 ], kD ));


It appears that that (errors_[0] = 0 giving newDt = 0) might have been a problem earlier, though

Agreed. I only remarked on this possibility because we're in the area now so it's an opportunity to fix it. For what it's worth, I tried this PR on a field case and ran into that exact problem. Otherwise I might not have noticed it.

alfbr · 2024-04-25T11:13:27Z

Any chance we can get this in?

bska · 2024-04-25T11:42:50Z

Any chance we can get this in?

There's absolutely a chance we can get it in. We're a little constrained by the way we handle updates to reference solutions in the time leading up to a new official release, so it'll definitely be easier to get it in once the new release is finalised.

bska · 2024-04-25T13:47:54Z

jenkins build this please

bska · 2024-04-25T13:53:38Z

jenkins build this please

Note to @erikhide : The above build request is just to refresh my memory of the impact of this PR. We are going to merge this regardless.

Fixed indexing error in PID controller

c160416

bska reviewed Oct 16, 2023

View reviewed changes

Added errors_[0] == 0 to conditional

2b7f391

bska reviewed Oct 17, 2023

View reviewed changes

Merge branch 'OPM:master' into master

366e558

Merge branch 'OPM:master' into master

a9da973

bska mentioned this pull request Apr 26, 2024

Include rejection of time step if tolerance test fails #5317

Draft

erikhide mentioned this pull request Apr 26, 2024

Changed the relative change used in adaptive time stepping #4856

Draft

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fixed indexing error in PID controller #4855

Fixed indexing error in PID controller #4855

erikhide commented Sep 8, 2023 •

edited by bska

bska commented Sep 8, 2023

bska commented Sep 8, 2023

erikhide commented Sep 8, 2023 •

edited

totto82 commented Sep 8, 2023

alfbr commented Sep 13, 2023

atgeirr commented Sep 13, 2023

ytelses commented Sep 14, 2023

totto82 commented Sep 14, 2023

erikhide commented Sep 14, 2023

bska commented Sep 14, 2023

bska Oct 16, 2023

erikhide Oct 16, 2023 •

edited

bska Oct 16, 2023

erikhide Oct 17, 2023

bska Oct 17, 2023

bska commented Oct 17, 2023

bska left a comment

bska Oct 17, 2023

alfbr commented Apr 25, 2024

bska commented Apr 25, 2024

bska commented Apr 25, 2024

bska commented Apr 25, 2024

Fixed indexing error in PID controller #4855

Are you sure you want to change the base?

Fixed indexing error in PID controller #4855

Conversation

erikhide commented Sep 8, 2023 • edited by bska

bska commented Sep 8, 2023

bska commented Sep 8, 2023

erikhide commented Sep 8, 2023 • edited

totto82 commented Sep 8, 2023

alfbr commented Sep 13, 2023

atgeirr commented Sep 13, 2023

ytelses commented Sep 14, 2023

totto82 commented Sep 14, 2023

erikhide commented Sep 14, 2023

bska commented Sep 14, 2023

bska Oct 16, 2023

Choose a reason for hiding this comment

erikhide Oct 16, 2023 • edited

Choose a reason for hiding this comment

bska Oct 16, 2023

Choose a reason for hiding this comment

erikhide Oct 17, 2023

Choose a reason for hiding this comment

bska Oct 17, 2023

Choose a reason for hiding this comment

bska commented Oct 17, 2023

bska left a comment

Choose a reason for hiding this comment

bska Oct 17, 2023

Choose a reason for hiding this comment

alfbr commented Apr 25, 2024

bska commented Apr 25, 2024

bska commented Apr 25, 2024

bska commented Apr 25, 2024

erikhide commented Sep 8, 2023 •

edited by bska

erikhide commented Sep 8, 2023 •

edited

erikhide Oct 16, 2023 •

edited