Skip to content
GitLab
Commit 8d8fe677 authored by HEMERY Mathieu's avatar HEMERY Mathieu
Browse files

Reflow and cleaning

parent bbe7783c
Hide whitespace changes
Inline Side-by-side
rosenbrock.pl
View file @ 8d8fe677
......@@ -5,9 +5,10 @@
).
% Only for separate compilation and linting
:- use_module(arithmetic_rules).
:- use_module(counters).
:- use_module(doc).
:- use_module(util).
:- use_module(arithmetic_rules).
%! rosenbrock_init
%
......@@ -348,306 +349,316 @@ rw_coeff(sin(Element),CurrentState,Parameters,Time,sin(ElementEvaluated)) :-
eval_state(State, Parameters, Evaluated) :-
eval_row(State, [], Parameters, Evaluated).
/*
eval_state([],_,[]).
eval_state([p(ParameterIndex)|Elements],Parameters,[NewElement|NewElements]) :-
nth0(ParameterIndex,Parameters,NewElement),!,
eval_state(Elements,Parameters,NewElements).
eval_state([Element|Elements],Parameters,[Element|NewElements]) :-
eval_state(Elements,Parameters,NewElements).
*/
/* This a tool for the rosenbrock method */
%! update_state(+Old, +StepSize, +Grad, -New)
%
% Compute New = Old + StepSize x Grad
% This a tool for the rosenbrock method
update_state([],_,[],[]).
update_state([Element|CurrentState],Constant,[Value|Values],[NewElement|NewCurrentState]) :-
NewElement is Element + Constant * Value,
update_state(CurrentState,Constant,Values,NewCurrentState).
NewElement is Element + Constant * Value,
update_state(CurrentState,Constant,Values,NewCurrentState).
get_last_row(Data,Time,CurrentState):-
last(Data, LastRow),
LastRow=..[row,Time|CurrentState].
/* Can count how many rosenbrock steps were done */
%! get_last_row(+Data, -Time, -CurrentState)
%
%
get_last_row(Data,Time,CurrentState):-
last(Data, LastRow),
LastRow=..[row,Time|CurrentState].
counter(Times) :-
nb_getval(counterr,Time),
Times is Time + 1,
nb_setval(counterr,Times).
/* Linear interpolation of between the two final points of the simulation, in order to not go further the final time desired */
%! 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).
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).
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 */
/* Parameters are stored within the global variable parameters_list
The current state of the simulation if obtained in the final row of the global variable numerical_table */
%! rosenbrock(+Equations, +InitialTime, +Duration, +Jacobian)
%
% Parameters are stored within the global variable parameters_list
% The current state of the simulation if obtained in the final row of the global
% variable numerical_table
rosenbrock(Equations,InitialTime,Duration,Jacobian) :-
set_counter(counterr,0),
FinalTime is InitialTime + Duration,
StepSizeMax is Duration/20,
(
StepSizeMax =:= 0.0,
!,
debug(error, "Simulation length cannot be 0", []),
fail
;
true
),
nb_setval(hmax,StepSizeMax),
(
nb_current(h,_)
->
true
;
StepSize is StepSizeMax/5,
nb_setval(h,StepSize),
nb_setval(k_time_units,2)
),
repeat,
% counter(TTTTT), write(TTTTT), nl,
nb_getval(parameters_list,Parameters),
nb_getval(numerical_table,CurrentTable),
get_last_row(CurrentTable,Time,CurrentStateInitial),
debug(debug,"entering step for ~p",[Time]),
length(Equations,NumberTrueVariable), % It prevents the functions values to be used in the simulation step
take(NumberTrueVariable,CurrentStateInitial,CurrentState),
rosenbrock_step(Equations,Jacobian,CurrentState,Parameters,CurrentState2),
nb_getval(hdid,NewStepSize),
Time2 is Time+NewStepSize,
(
check_events(CurrentState2,Time2),
check_conditional_events(CurrentState2,Time2),
set_counter(rosenbrock_step,0),
FinalTime is InitialTime + Duration,
StepSizeMax is Duration/20,
(
StepSizeMax =:= 0.0,
!,
debug(error, "Simulation length cannot be 0", []),
fail
;
true
),
nb_setval(hmax,StepSizeMax),
(
nb_current(h, _) % same as nb_getval but fail if 'h' do not exist instead of crashing
->
true
;
StepSize is StepSizeMax/5,
nb_setval(h,StepSize),
nb_setval(k_time_units,2)
),
repeat,
% count(rosenbrock_step, STEP), write_ln(STEP),
nb_getval(parameters_list, Parameters),
nb_getval(numerical_table, CurrentTable),
get_last_row(CurrentTable, Time, State_full),
debug(rosenbrock, "entering step for ~p", [Time]),
length(Equations, NumberTrueVariable), % Cut the function values at the end of State_full
take(NumberTrueVariable, State_full, CurrentState),
% Here is the true rosenbrock step
rosenbrock_step(Equations, Jacobian, CurrentState, Parameters, Next_state),
nb_getval(hdid,NewStepSize),
Time2 is Time+NewStepSize,
(
check_events(Next_state,Time2),
check_conditional_events(Next_state,Time2),
nb_getval(events_number,NumberOfEvents),
NumberOfEvents > 0
->
->
euler(Equations,Time,CurrentState,CurrentTable)
;
;
(
nb_current(functions_list,_)
nb_current(functions_list,_)
->
update_functions_state(CurrentState2,Time2)
update_functions_state(Next_state,Time2)
;
NewLine =.. [row,Time2|CurrentState2],
append(CurrentTable,[NewLine],NewTable),
nb_setval(numerical_table,NewTable)
NewLine =.. [row,Time2|Next_state],
append(CurrentTable,[NewLine],NewTable),
nb_setval(numerical_table,NewTable)
),
nb_setval(last_time,Time2)
),
nb_getval(last_time,LastTime),
LastTime >= FinalTime,!,
final_interpolation(FinalTime).
nb_setval(last_time,Time2)
),
nb_getval(last_time,LastTime),
LastTime >= FinalTime,!,
final_interpolation(FinalTime).
%! rosenbrock_step(+Equations, +Jacobian, +CurrentState, +Parameters, -NextState)
%
% perform one step of rosenbrock integration, take also care of adjusting the step_size
% (stored in global variable nb_getval(hdid, _))
rosenbrock_step(Equations,Jacobian,CurrentState,Parameters,NewCurrentState) :-
nb_getval(k_rbk_safe,SAFE),
nb_getval(k_rbk_grow,GROW),
nb_getval(k_rbk_pgrow,PGROW),
nb_getval(k_rbk_shrink,SHRINK),
nb_getval(k_rbk_pshrink,PSHRINK),
nb_getval(k_rbk_errcon,ERRCON),
nb_getval(k_rbk_gam,GAM),
nb_getval(k_rbk_a21,A21),
nb_getval(k_rbk_a31,A31),
nb_getval(k_rbk_a32,A32),
nb_getval(k_rbk_c21,C21),
nb_getval(k_rbk_c31,C31),
nb_getval(k_rbk_c32,C32),
nb_getval(k_rbk_c41,C41),
nb_getval(k_rbk_c42,C42),
nb_getval(k_rbk_c43,C43),
nb_getval(k_rbk_b1,B1),
nb_getval(k_rbk_b2,B2),
nb_getval(k_rbk_b3,B3),
nb_getval(k_rbk_b4,B4),
nb_getval(k_rbk_e1,E1),
nb_getval(k_rbk_e2,E2),
nb_getval(k_rbk_e3,E3),
nb_getval(k_rbk_e4,E4),
repeat,
eval_state(CurrentState,Parameters,TrueCurrentState),
eval_matrix(Jacobian,TrueCurrentState,Parameters,JacobianCurrentState), % Create a choicepoint which I can't find a way to prevent
nb_getval(h,StepSize),
nb_setval(hdid,StepSize),
length(Equations,Length),
% evaluate Jacobian at current time point
debug(debug,"Jacobian evaluated",[]), %%%%%%% j'en suis là
debug(debug,"Jacobian = ~w CurrentState = ~w",[Jacobian, TrueCurrentState]),
debug(debug,"JacobianCurrentState = ~w",[JacobianCurrentState]),
% fill the left hand side of AX=B
GGAM is 1.0/(GAM*StepSize),
create_diag_matrix(Length,GGAM,_),
sub_matrix(JacobianCurrentState),
% LU decompose A
ludcmp_matrix,
debug(debug,"Left side evaluated",[]),
% fill the right hand side, G1
eval_row(Equations,TrueCurrentState,Parameters,G1),
debug(debug,"Right side evaluated",[]),
debug(debug,"Equations = ~w CurrentState = ~w",[Equations, TrueCurrentState]),
debug(debug,"G1 = ~w",[G1]),
% compute by backsubstitution
lusubst(G1,G1X),
debug(debug,"G1=~p",[G1X]),
% implement inc_vector/4 and inc_val_list/4 (which will be the same predicate I think)
% get the new values at time + A2X * StepSize
update_state(TrueCurrentState,A21,G1X,CurrentState1),
debug(debug,"CurrentState1 = ~w",[CurrentState1]),
eval_row(Equations,CurrentState1,Parameters,Equations1),
debug(debug,"Equations at +A2X = ~p",[Equations1]),
% fill rhs again for G2
C21H is C21/StepSize,
update_state(Equations1,C21H,G1X,G2),
debug(debug,"Right side evaluated",[]),
% compute by backsubstitution
lusubst(G2,G2X),
debug(debug,"G2=~p",[G2X]),
% get the new values at time + (A2X+A3X) * StepSize
update_state(TrueCurrentState,A31,G1X,CurrentState21),
update_state(CurrentState21,A32,G2X,CurrentState2),
eval_row(Equations,CurrentState2,Parameters,Equations2),
debug(debug,"Equations at +A2X+A3X = ~p",[Equations2]),
% fill rhs again for G3
C31H is C31/StepSize,
C32H is C32/StepSize,
update_state(Equations2,C31H,G1X,G31),
update_state(G31,C32H,G2X,G3),
debug(debug,"Right side evaluated",[]),
% compute by backsubstitution
lusubst(G3,G3X),
debug(debug,"G3=~p",[G3X]),
% for G4 no need to recompute CurrentState and Equations
C41H is C41/StepSize,
C42H is C42/StepSize,
C43H is C43/StepSize,
update_state(Equations2,C41H,G1X,G41),
update_state(G41,C42H,G2X,G42),
update_state(G42,C43H,G3X,G4),
debug(debug,"Right side evaluated",[]),
% compute by backsubstitution
lusubst(G4,G4X),
debug(debug,"G4=~p",[G4X]),
% now get 4th order estimate for CurrentState
update_state(TrueCurrentState,B1,G1X,CurrentState41),
update_state(CurrentState41,B2,G2X,CurrentState42),
update_state(CurrentState42,B3,G3X,CurrentState43),
update_state(CurrentState43,B4,G4X,NewCurrentState),
debug(debug,"CurrentState 4th order evaluated",[]),
% Estimate error
constant_list(Length,0,ErrorListInit), nb_setval(error,ErrorListInit),
nb_getval(error,ErrorList0),
update_state(ErrorList0,E1,G1X,ErrorList1),
update_state(ErrorList1,E2,G2X,ErrorList2),
update_state(ErrorList2,E3,G3X,ErrorList3),
update_state(ErrorList3,E4,G4X,ErrorList4),
scale_error(ErrorList4,TrueCurrentState,ErrorList),
max_list(ErrorList,Erro),
debug(debug,"Error computed",[]),
debug(debug,"Erro = ~w",[Erro]),
% check error
nb_getval(step_doubling_error,ErrMax),
Error is Erro/ErrMax,
debug(debug,"Error = ~w",[Error]),
nb_getval(k_rbk_safe,SAFE),
nb_getval(k_rbk_grow,GROW),
nb_getval(k_rbk_pgrow,PGROW),
nb_getval(k_rbk_shrink,SHRINK),
nb_getval(k_rbk_pshrink,PSHRINK),
nb_getval(k_rbk_errcon,ERRCON),
nb_getval(k_rbk_gam,GAM),
nb_getval(k_rbk_a21,A21),
nb_getval(k_rbk_a31,A31),
nb_getval(k_rbk_a32,A32),
nb_getval(k_rbk_c21,C21),
nb_getval(k_rbk_c31,C31),
nb_getval(k_rbk_c32,C32),
nb_getval(k_rbk_c41,C41),
nb_getval(k_rbk_c42,C42),
nb_getval(k_rbk_c43,C43),
nb_getval(k_rbk_b1,B1),
nb_getval(k_rbk_b2,B2),
nb_getval(k_rbk_b3,B3),
nb_getval(k_rbk_b4,B4),
nb_getval(k_rbk_e1,E1),
nb_getval(k_rbk_e2,E2),
nb_getval(k_rbk_e3,E3),
nb_getval(k_rbk_e4,E4),
repeat,
% evaluate Jacobian at current time point
eval_state(CurrentState,Parameters,TrueCurrentState),
eval_matrix(Jacobian,TrueCurrentState,Parameters,JacobianCurrentState), % Create a choicepoint which I can't find a way to prevent
debug(rosenbrock,"Jacobian evaluated",[]),
debug(rosenbrock,"Jacobian = ~w CurrentState = ~w",[Jacobian, TrueCurrentState]),
debug(rosenbrock,"JacobianCurrentState = ~w",[JacobianCurrentState]),
nb_getval(h,StepSize),
nb_setval(hdid,StepSize),
length(Equations,Length),
% fill the left hand side of AX=B
GGAM is 1.0/(GAM*StepSize),
create_diag_matrix(Length,GGAM,_),
sub_matrix(JacobianCurrentState),
% LU decompose A
ludcmp_matrix,
debug(rosenbrock,"Left side evaluated",[]),
% fill the right hand side, G1
eval_row(Equations,TrueCurrentState,Parameters,G1),
debug(rosenbrock,"Right side evaluated",[]),
debug(rosenbrock,"Equations = ~w CurrentState = ~w",[Equations, TrueCurrentState]),
debug(rosenbrock,"G1 = ~w",[G1]),
% compute by backsubstitution
lusubst(G1,G1X),
debug(rosenbrock,"G1=~p",[G1X]),
% implement inc_vector/4 and inc_val_list/4 (which will be the same predicate I think)
% get the new values at time + A2X * StepSize
update_state(TrueCurrentState,A21,G1X,CurrentState1),
debug(rosenbrock,"CurrentState1 = ~w",[CurrentState1]),
eval_row(Equations,CurrentState1,Parameters,Equations1),
debug(rosenbrock,"Equations at +A2X = ~p",[Equations1]),
% fill rhs again for G2
C21H is C21/StepSize,
update_state(Equations1,C21H,G1X,G2),
debug(rosenbrock,"Right side evaluated",[]),
% compute by backsubstitution
lusubst(G2,G2X),
debug(rosenbrock,"G2=~p",[G2X]),
% get the new values at time + (A2X+A3X) * StepSize
update_state(TrueCurrentState,A31,G1X,CurrentState21),
update_state(CurrentState21,A32,G2X,CurrentState2),
eval_row(Equations,CurrentState2,Parameters,Equations2),
debug(rosenbrock,"Equations at +A2X+A3X = ~p",[Equations2]),
% fill rhs again for G3
C31H is C31/StepSize,
C32H is C32/StepSize,
update_state(Equations2,C31H,G1X,G31),
update_state(G31,C32H,G2X,G3),
debug(rosenbrock,"Right side evaluated",[]),
% compute by backsubstitution
lusubst(G3,G3X),
debug(rosenbrock,"G3=~p",[G3X]),
% for G4 no need to recompute CurrentState and Equations
C41H is C41/StepSize,
C42H is C42/StepSize,
C43H is C43/StepSize,
update_state(Equations2,C41H,G1X,G41),
update_state(G41,C42H,G2X,G42),
update_state(G42,C43H,G3X,G4),
debug(rosenbrock,"Right side evaluated",[]),
% compute by backsubstitution
lusubst(G4,G4X),
debug(rosenbrock,"G4=~p",[G4X]),
% now get 4th order estimate for CurrentState
update_state(TrueCurrentState,B1,G1X,CurrentState41),
update_state(CurrentState41,B2,G2X,CurrentState42),
update_state(CurrentState42,B3,G3X,CurrentState43),
update_state(CurrentState43,B4,G4X,NewCurrentState),
debug(rosenbrock,"CurrentState 4th order evaluated",[]),
% Estimate error
constant_list(Length,0,ErrorListInit), nb_setval(error,ErrorListInit),
nb_getval(error,ErrorList0),
update_state(ErrorList0,E1,G1X,ErrorList1),
update_state(ErrorList1,E2,G2X,ErrorList2),
update_state(ErrorList2,E3,G3X,ErrorList3),
update_state(ErrorList3,E4,G4X,ErrorList4),
scale_error(ErrorList4,TrueCurrentState,ErrorList),
max_list(ErrorList,Erro),
debug(rosenbrock,"Error computed",[]),
debug(rosenbrock,"Erro = ~w",[Erro]),
% check error
nb_getval(step_doubling_error,ErrMax),
Error is Erro/ErrMax,
debug(rosenbrock,"Error = ~w",[Error]),
(
(
(
Error > 1
;
(Error = nan) % catch the case where Error=nan
)
Error > 1
;
(Error = nan) % catch the case where Error=nan
)
->
% failure: reduce the step
% TODO if not yet at maxtry
(
StepSize < ErrMax/1000
->
% failure: reduce the step
% TODO if not yet at maxtry
(
StepSize < ErrMax/1000
->
(
k_fail
->
true
;
assertz(k_fail),
debug(warning, "Failure with minimal step, error will not be guaranteed.", [])
)
;
StepSize1 is SAFE*StepSize*(Error**PSHRINK),
StepSize2 is SHRINK*StepSize,
max_list([StepSize1,StepSize2],StepSizeNext),
nb_setval(h,StepSizeNext),
debug(debug,"error ~p, candidates ~p ~p",[Erro,StepSize1,StepSize2]),
debug(debug,"step size ~p failed, trying ~p...",[StepSize,StepSizeNext]),
fail
)
(
k_fail
->
true
;
assertz(k_fail),
debug(warning, "Failure with minimal step, error will not be guaranteed.", [])
)
;
% success, try increasing the step size
(
Error > ERRCON*ErrMax
->
StepSizeNext is SAFE*StepSize*(Error**PGROW)
;
StepSizeNext is StepSize*GROW
),
debug(debug,"step size ~p ok, trying ~p...",[StepSize,StepSizeNext]),
nb_getval(hmax,StepSizeMAX),
(
StepSizeNext < StepSizeMAX
->
nb_setval(h,StepSizeNext)
;
true
)
StepSize1 is SAFE*StepSize*(Error**PSHRINK),
StepSize2 is SHRINK*StepSize,
max_list([StepSize1,StepSize2],StepSizeNext),
nb_setval(h,StepSizeNext),
debug(rosenbrock,"error ~p, candidates ~p ~p",[Erro,StepSize1,StepSize2]),
debug(rosenbrock,"step size ~p failed, trying ~p...",[StepSize,StepSizeNext]),
fail
)
;
% success, try increasing the step size
(
Error > ERRCON*ErrMax
->
StepSizeNext is SAFE*StepSize*(Error**PGROW)
;
StepSizeNext is StepSize*GROW
),
debug(rosenbrock,"step size ~p ok, trying ~p...",[StepSize,StepSizeNext]),
nb_getval(hmax,StepSizeMAX),
(
StepSizeNext < StepSizeMAX
->
nb_setval(h,StepSizeNext)
;
true
)
),
!,
true.
!,
true.
/* Scale error vector with absolute value of variable (in Name,Value list)
if that |value| is > 1 */
%! scale_error(+Errors, +Values, -Scaled_errors)
%
% Scale error vector with absolute value of variable (if |value| is > 1)
scale_error([],[],[]).
scale_error([E|Err],[V|L],[SE|SErr]) :-
(
abs(V) < 1
->
SE is abs(E)
;
SE is abs(E/V)
),
scale_error(Err,L,SErr).
(
abs(V) < 1
->
SE is abs(E)
;
SE is abs(E/V)
),
scale_error(Err,L,SErr).
/* Events handling */
......@@ -750,20 +761,6 @@ check_conditional_events(State,Time) :-
/* Same for traditionnal events */
/*check_events(State,Time) :-
get_events_lists(Events,_),
set_counter(events_number,0),
forall(
(member(Event, Events),Event =.. [->,Condition|_]),
( test_condition(State,Time,Condition,Result),
Result
->
count(events_number,_)
;
true
)
).*/
check_events(State,Time) :-
get_events_lists(Events,_),
set_counter(events_number,0),
......
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