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