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