bug in CTimeDiscretization::CalcTimeStep

Dear all,

we found an issue with the "Critical time match" in CTimeDiscretization::CalcTimeStep(double current_time).

As I understand this functionality, it is intended to make sure that no critical time for model output is missed in simulations with automatic time step control. Generally, this is a very good feature, especially if it can be extended to include critical times for transient BCs or STs:)
In the case of a model with a finite number of explicitely defined - but variable - time step lengths (i.e. no automatic control) like
\$TIME_STEPS
1 500
1 10000
5 10
the current implementation however may cause a reduction of a single time step ahead of a critical time for model output, which is not compensated for by prolongation of the following time step, as the (maximum) length of each time step is defined and thus fixed in the time step list. This will result in a loss of total simulation time, and worse, the variable predefined time step lengths no longer match transient BC/ST changes defined by time curves or functions.

I think the main problem is with the "else if" statement below:

CTimeDiscretization::CalcTimeStep(double current_time){
...
// WW. Critical time match (JT2012 modified)
for(int i = 0; i < (int)critical_time.size(); i++){
if(current_time < critical_time[i]) {

next = current_time + time_step_length;
// JT2012. A tiny increase in dt is better
// than a miniscule dt on the next step
tval = next + time_step_length/1.0e3;

if(tval > critical_time[i]){
// Critical time is hit
if(next != critical_time[i]){
// otherwise, match is already exact
time_step_length = (critical_time[i] - current_time);
}
break;
}
else if(tval + time_step_length > critical_time[i]){
// We can hit the critical time in 2 time steps,
// smooth the transition
if(next + time_step_length != critical_time[i]){
// otherwise, match is already exact
time_step_length = (critical_time[i] - current_time)/2.0;
}
break;
}
break;
}
}
...
next_active_time = current_time + time_step_length;
return time_step_length;
}

With the example of the list above:
After the first time step of 500 s we have a long time step of 10000 s followed by 5 short time steps of 10 s and model output requested at the end of the simulation at t = 10550 s.
i.e.
current time = 500
time_step_length = 10000
critical_time = 10550

next = current_time + time_step_length;
= 10500
tval = next + time_step_length/1.0e3;
= 10510

In the "if" statement tval < critical_time
so nothing happens here.
In the "else if" statement we compare
tval + time_step_length > critical_time
10510 + 10000 > 10550
which is true. Hence, after the following "if":
next + time_step_length != critical_time[i]
10500 + 10000 != 10550
--> time_step_length = (critical_time[i] - current_time)/2.0;
5025 = (10550 - 500)/2

Accordingly, the simulation will be advanced by 5025 s.
At the new time step, the numbers now are:
current time = 5525
time_step_length = 10
critical_time = 10550 (still)
next = 5535
tval = 5535.01

Accordingly
"if": tval < critical_time
and
"else if": tval + time_step_length < critical_time

The time_step_length hence will stay at 10 s, also for the following 4 steps and the simulation will not reach the intended final time of 10550 s.

A quick fix would be to restrict the critical time match functionality to automatic time stepping simulations and leave it to the responsibility of the user to make sure that simulation and output time steps (as well as transient BC/ST) are synchronized.

What do you think?
Christof

Dear Christof,

I support your idea to limit the critical time match feature to automatic time stepping. I think, that the user may take care to get output at the desired time steps, if she/he defines the stepping manually.
However, did you test, whether this issue may occur with the automatic stepping, too?

Best,
Marc

···

On 22.10.2015 13:17, Christof Beyer wrote:

Dear all,

we found an issue with the "Critical time match" in
CTimeDiscretization::CalcTimeStep(double current_time).

As I understand this functionality, it is intended to make sure that no
critical time for model output is missed in simulations with automatic
time step control. Generally, this is a very good feature, especially if
it can be extended to include critical times for transient BCs or STs:)
In the case of a model with a finite number of explicitely defined - but
variable - time step lengths (i.e. no automatic control) like
\$TIME_STEPS
1 500
1 10000
5 10
the current implementation however may cause a reduction of a single
time step ahead of a critical time for model output, which is not
compensated for by prolongation of the following time step, as the
(maximum) length of each time step is defined and thus fixed in the time
step list. This will result in a loss of total simulation time, and
worse, the variable predefined time step lengths no longer match
transient BC/ST changes defined by time curves or functions.

I think the main problem is with the "else if" statement below:

CTimeDiscretization::CalcTimeStep(double current_time){
...
// WW. Critical time match (JT2012 modified)
for(int i = 0; i < (int)critical_time.size(); i++){
if(current_time < critical_time[i]) {

next = current_time + time_step_length;
// JT2012. A tiny increase in dt is better
// than a miniscule dt on the next step
tval = next + time_step_length/1.0e3;

if(tval > critical_time[i]){
// Critical time is hit
if(next != critical_time[i]){
// otherwise, match is already exact
time_step_length = (critical_time[i] - current_time);
}
break;
}
else if(tval + time_step_length > critical_time[i]){
// We can hit the critical time in 2 time steps,
// smooth the transition
if(next + time_step_length != critical_time[i]){
// otherwise, match is already exact
time_step_length = (critical_time[i] - current_time)/2.0;
}
break;
}
break;
}
}
...
next_active_time = current_time + time_step_length;
return time_step_length;
}

With the example of the list above:
After the first time step of 500 s we have a long time step of 10000 s
followed by 5 short time steps of 10 s and model output requested at the
end of the simulation at t = 10550 s.
i.e.
current time = 500
time_step_length = 10000
critical_time = 10550

next = current_time + time_step_length;
= 10500
tval = next + time_step_length/1.0e3;
= 10510

In the "if" statement tval < critical_time
so nothing happens here.
In the "else if" statement we compare
tval + time_step_length > critical_time
10510 + 10000 > 10550
which is true. Hence, after the following "if":
next + time_step_length != critical_time[i]
10500 + 10000 != 10550
--> time_step_length = (critical_time[i] - current_time)/2.0;
5025 = (10550 - 500)/2

Accordingly, the simulation will be advanced by 5025 s.
At the new time step, the numbers now are:
current time = 5525
time_step_length = 10
critical_time = 10550 (still)
next = 5535
tval = 5535.01

Accordingly
"if": tval < critical_time
and
"else if": tval + time_step_length < critical_time

The time_step_length hence will stay at 10 s, also for the following 4
steps and the simulation will not reach the intended final time of 10550 s.

A quick fix would be to restrict the critical time match functionality
to automatic time stepping simulations and leave it to the
responsibility of the user to make sure that simulation and output time
steps (as well as transient BC/ST) are synchronized.

What do you think?
Christof

--
- - - - - - - - - - - - - - - - - -

Marc Walther
Hydrologist
JProf Dr. rer. nat.

Helmholtz Centre for Environmental Research - UFZ
Department of Environmental Informatics
Permoserstraße 15, 04318 Leipzig, Germany

Technische Universität Dresden
Faculty of Environmental Sciences
Department of Hydrosciences
Institute for Groundwater Management
01062 Dresden

Permoserstraße 15, 04318 Leipzig
Building 7.2, Room E110

Bergstr. 66, 01069 Dresden
Neubau Chemie 2. Bauabschnitt, CHE-BA2, Room E60

Phone +49 341 235 10 54
Fax +49 341 235 19 39
Mobile +49 178 334 18 19

Helmholtz-Zentrum für Umweltforschung GmbH - UFZ
Registered Office: Leipzig
Registration Office: Amtsgericht Leipzig
Chairman of the Supervisory Board: MinR Wilfried Kraus
Scientific Director: Prof. Dr. Georg Teutsch

- - - - - - - - - - - - - - - - - -

Yes, we have to avoid using 'Critical time' for the manually specified (or fixed) time steps. The reason is that we can adjust the step sizes manually for the times when we want to get solutions.

···

On 10/22/2015 01:17 PM, Christof Beyer wrote:

Dear all,

we found an issue with the "Critical time match" in CTimeDiscretization::CalcTimeStep(double current_time).

As I understand this functionality, it is intended to make sure that no critical time for model output is missed in simulations with automatic time step control. Generally, this is a very good feature, especially if it can be extended to include critical times for transient BCs or STs:)
In the case of a model with a finite number of explicitely defined - but variable - time step lengths (i.e. no automatic control) like
\$TIME_STEPS
1 500
1 10000
5 10
the current implementation however may cause a reduction of a single time step ahead of a critical time for model output, which is not compensated for by prolongation of the following time step, as the (maximum) length of each time step is defined and thus fixed in the time step list. This will result in a loss of total simulation time, and worse, the variable predefined time step lengths no longer match transient BC/ST changes defined by time curves or functions.

I think the main problem is with the "else if" statement below:

CTimeDiscretization::CalcTimeStep(double current_time){
...
// WW. Critical time match (JT2012 modified)
for(int i = 0; i < (int)critical_time.size(); i++){
if(current_time < critical_time[i]) {

next = current_time + time_step_length;
// JT2012. A tiny increase in dt is better
// than a miniscule dt on the next step
tval = next + time_step_length/1.0e3;

if(tval > critical_time[i]){
// Critical time is hit
if(next != critical_time[i]){
// otherwise, match is already exact
time_step_length = (critical_time[i] - current_time);
}
break;
}
else if(tval + time_step_length > critical_time[i]){
// We can hit the critical time in 2 time steps,
// smooth the transition
if(next + time_step_length != critical_time[i]){
// otherwise, match is already exact
time_step_length = (critical_time[i] - current_time)/2.0;
}
break;
}
break;
}
}
...
next_active_time = current_time + time_step_length;
return time_step_length;
}

With the example of the list above:
After the first time step of 500 s we have a long time step of 10000 s followed by 5 short time steps of 10 s and model output requested at the end of the simulation at t = 10550 s.
i.e.
current time = 500
time_step_length = 10000
critical_time = 10550

next = current_time + time_step_length;
= 10500
tval = next + time_step_length/1.0e3;
= 10510

In the "if" statement tval < critical_time
so nothing happens here.
In the "else if" statement we compare
tval + time_step_length > critical_time
10510 + 10000 > 10550
which is true. Hence, after the following "if":
next + time_step_length != critical_time[i]
10500 + 10000 != 10550
--> time_step_length = (critical_time[i] - current_time)/2.0;
5025 = (10550 - 500)/2

Accordingly, the simulation will be advanced by 5025 s.
At the new time step, the numbers now are:
current time = 5525
time_step_length = 10
critical_time = 10550 (still)
next = 5535
tval = 5535.01

Accordingly
"if": tval < critical_time
and
"else if": tval + time_step_length < critical_time

The time_step_length hence will stay at 10 s, also for the following 4 steps and the simulation will not reach the intended final time of 10550 s.

A quick fix would be to restrict the critical time match functionality to automatic time stepping simulations and leave it to the responsibility of the user to make sure that simulation and output time steps (as well as transient BC/ST) are synchronized.

What do you think?
Christof