OpenCMISS / iron

Source code repository for OpenCMISS-Iron
9 stars 62 forks source link

Problem_Control_Loop_Solve using Ceiling for time_loop%number_of_iterations #185

Open jaredmcollette opened 4 years ago

jaredmcollette commented 4 years ago

In working with reaction diffusion routines, I noticed that FiniteElementCalculate seemed to be called either once or twice per time step (per chemical species/per element) arbitrarily within the same simulation and would switch between being called once or twice without any pattern.

Tracing this bug to the root, I found the variable TIME_LOOP%TIME_INCREMENT would solve to values between 0.99999999999999944 and 1.0000000000000002, but because the CEILING Function is used, this would round to 1 and 2 respectively, when the value should clearly be 1 in both cases.

Switching the function CEILING to instead use NINT (line 464 in Problem_Routines) fixed the bug in my simulation, but I am unsure if it would cause issues in other cases. Was there a particular motivation for using CEILING for the TIME_LOOP%TIME_INCREMENT variable?

chrispbradley commented 4 years ago

That code was put in to try and make the time loops consistent in that due to small floating point errors when you kept adding time increment to start time it would sometime be just short of stop time and thus do an extra iteration.With CEILING it will ensure that the time is always greater than stop time on exit. With NINT there may be a case that it rounds down and so stop time is not reached when the loop exits. I can see pros and cons for each method. That said your case does seem wrong that we have 2 iterations instead of 1. Does anybody else have any thoughts?

vraj004 commented 4 years ago

Perhaps one way to deal with this is to not use NINT or CEILING but to test if adding an additional time step to the current time would end up overshooting the end time. One could even set a tolerance on how far that calculation can overshoot.

chrispbradley commented 4 years ago

I guess the most important thing is that the number of steps the time loop takes is deterministic. We really want different runs to always execute the same number of steps. With whatever choice of NINT or CEILING that you choose you may end up with different results due to rounding and other numerical errors. There is the command

controlLoop.NumberOfIterationsSet(number)

which will then redefine the time increment to be closer to what is required to achieve that number of iterations between the start time and stop time?

vraj004 commented 4 years ago

Related to your comment - and similar to mechanics control loops - maybe we should set START_TIME, END_TIME and TIME_STEPS instead of the time step size to solve this issue.