Commit be9b1b66 authored by Mathieu Hemery's avatar Mathieu Hemery
Browse files

Avoid final_interpolation trick by adjusting rosenbrock on the fly

parent c4edc110
...@@ -362,29 +362,6 @@ get_last_row(Data,Time,CurrentState):- ...@@ -362,29 +362,6 @@ get_last_row(Data,Time,CurrentState):-
LastRow=..[row,Time|CurrentState]. LastRow=..[row,Time|CurrentState].
%! final_interpolation(+Final_time)
%
% Linear interpolation of between the two final points of the simulation,
% in order to not go further the final time desired.
% /!\ Check if we can not just ask to rosenbrock not to overtake the final time
final_interpolation(_,[],_,[],_,[]).
final_interpolation(Time1,[Value1|Data1],Time2,[Value2|Data2],FinalTime,[NewValue|NewData]) :-
NewValue is (Value2 * (FinalTime - Time1) + Value1 * (Time2 - FinalTime)) / (Time2 - Time1),
final_interpolation(Time1,Data1,Time2,Data2,FinalTime,NewData).
final_interpolation(FinalTime) :-
nb_getval(numerical_table,Data),
get_last_row(Data,Time1,State1),
last(Data,LastRow),
append(Data2,[LastRow],Data),
get_last_row(Data2,Time2,State2),
final_interpolation(Time1,State1,Time2,State2,FinalTime,FinalValues),
FinalRow =.. [row,FinalTime|FinalValues],
append(Data2,[FinalRow],FinalData),
nb_setval(numerical_table,FinalData).
/*********************/ /*********************/
/* Rosenbrock method */ /* Rosenbrock method */
/*********************/ /*********************/
...@@ -426,7 +403,8 @@ rosenbrock(Equations,InitialTime,Duration,Jacobian) :- ...@@ -426,7 +403,8 @@ rosenbrock(Equations,InitialTime,Duration,Jacobian) :-
length(Equations, NumberTrueVariable), % Cut the function values at the end of State_full length(Equations, NumberTrueVariable), % Cut the function values at the end of State_full
separate_list(State_full, NumberTrueVariable, CurrentState, _FunctionState), separate_list(State_full, NumberTrueVariable, CurrentState, _FunctionState),
% Here is the true rosenbrock step % Here is the true rosenbrock step
rosenbrock_step(Equations, Jacobian, CurrentState, Parameters, Next_state), MaxStep is FinalTime - Time,
rosenbrock_step(Equations, Jacobian, CurrentState, Parameters, MaxStep, Next_state),
nb_getval(hdid,NewStepSize), nb_getval(hdid,NewStepSize),
Time2 is Time+NewStepSize, Time2 is Time+NewStepSize,
( (
...@@ -449,16 +427,16 @@ rosenbrock(Equations,InitialTime,Duration,Jacobian) :- ...@@ -449,16 +427,16 @@ rosenbrock(Equations,InitialTime,Duration,Jacobian) :-
nb_setval(last_time,Time2) nb_setval(last_time,Time2)
), ),
nb_getval(last_time,LastTime), nb_getval(last_time,LastTime),
LastTime >= FinalTime,!, LastTime >= FinalTime,!.
final_interpolation(FinalTime).
%! rosenbrock_step(+Equations, +Jacobian, +CurrentState, +Parameters, -NextState) %! rosenbrock_step(+Equations, +Jacobian, +CurrentState, +Parameters, +MaxStep, -NextState)
% %
% perform one step of rosenbrock integration, take also care of adjusting the step_size % perform one step of rosenbrock integration, take also care of adjusting the step_size
% (stored in global variable nb_getval(hdid, _)) % (stored in global variable nb_getval(hdid, _))
% MaxStep avoid overtaking the final step
rosenbrock_step(Equations,Jacobian,CurrentState,Parameters,NewCurrentState) :- rosenbrock_step(Equations,Jacobian,CurrentState,Parameters,MaxStep,NewCurrentState) :-
nb_getval(k_rbk_safe,SAFE), nb_getval(k_rbk_safe,SAFE),
nb_getval(k_rbk_grow,GROW), nb_getval(k_rbk_grow,GROW),
nb_getval(k_rbk_pgrow,PGROW), nb_getval(k_rbk_pgrow,PGROW),
...@@ -491,7 +469,8 @@ rosenbrock_step(Equations,Jacobian,CurrentState,Parameters,NewCurrentState) :- ...@@ -491,7 +469,8 @@ rosenbrock_step(Equations,Jacobian,CurrentState,Parameters,NewCurrentState) :-
debug(rosenbrock,"Jacobian = ~w CurrentState = ~w",[Jacobian, TrueCurrentState]), debug(rosenbrock,"Jacobian = ~w CurrentState = ~w",[Jacobian, TrueCurrentState]),
debug(rosenbrock,"JacobianCurrentState = ~w",[JacobianCurrentState]), debug(rosenbrock,"JacobianCurrentState = ~w",[JacobianCurrentState]),
nb_getval(h,StepSize), nb_getval(h,StepSizeRaw),
StepSize is min(StepSizeRaw, MaxStep),
nb_setval(hdid,StepSize), nb_setval(hdid,StepSize),
length(Equations,Length), length(Equations,Length),
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment