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

Start to optimize rosenbrock, core works but asides need now to be tuned

parent abaa7df5
......@@ -19,17 +19,14 @@ rosenbrock_init :-
nb_setval(not_zero_error,0.0001),
nb_setval(k_method,stiff),
nb_setval(euler_step_size,0.01),
nb_setval(step_doubling_error,0.0001),
SAFE = 0.9, GROW = 1.5, PGROW = -0.25,
nb_setval(step_doubling_error,1e-6), % Control the method's precision
SAFE = 0.9, GROW = 2.5, PGROW = -0.25,
SHRINK = 0.5, PSHRINK is -1.0/3.0, ERRCON = 0.1296,
%MAXTRY = 40,
GAM is 1.0/2.0, A21=2.0, A31 is 48.0/25.0, A32 is 6.0/25.0, C21 = -8.0,
C31 is 372.0/25.0, C32 is 12.0/5.0, C41 is -112.0/125.0,
C42 is -54.0/125.0, C43 is -2.0/5.0,
B1 is 19.0/9.0, B2 is 1.0/2.0, B3 is 25.0/108.0, B4 is 125.0/108.0,
E1 is 17.0/54.0, E2 is 7.0/36.0, E3 is 0.0, E4 is 125.0/108.0,
% C1X, C2X, C3X, C4X,
% A2X is 1.0, A3X is 3.0/5.0,
nb_setval(k_rbk_safe,SAFE),
nb_setval(k_rbk_grow,GROW),
nb_setval(k_rbk_pgrow,PGROW),
......@@ -82,7 +79,7 @@ ode_solver([
jacobian: Jacobian
]) :-
(
nb_current(rsbk_continue,yes)
nb_current(rsbk_continue,yes) %BRANCH TO BE CLEANED
->
nb_getval(numerical_table,CurrentTable),
get_last_row(CurrentTable,InitialTime2,TrueInitialState),
......@@ -104,14 +101,14 @@ ode_solver([
)
),
nb_getval(variable_list,VariableList),
FirstRow =.. [row,'#''Time'|VariableList],
assertz(saved_row(['#''Time'|VariableList])), %initialize the label of the numerical table
eval_state(InitialState,InitialParameters,TrueInitialState),
nb_setval(events_list,Events),
prepare_rosenbrock_events(TrueInitialState,InitialTime),
prepare_conditional_events(TrueInitialState,InitialTime),
SecondRow =.. [row,InitialTime|TrueInitialState],
InitialTable = [FirstRow,SecondRow],
nb_setval(numerical_table,InitialTable),
nb_setval(last_row, SecondRow),
assertz(saved_row([InitialTime|TrueInitialState])),
initialize_functions_list,
TrueInitialTime is InitialTime
),
......@@ -360,8 +357,8 @@ update_state(Old, Stepsize, Grad, New) :-
%
%
get_last_row(Data,Time,CurrentState):-
last(Data, LastRow),
get_last_row(Time,CurrentState):-
nb_getval(last_row, LastRow),
LastRow=..[row,Time|CurrentState].
......@@ -376,7 +373,6 @@ get_last_row(Data,Time,CurrentState):-
% variable numerical_table
rosenbrock(Equations,InitialTime,Duration,Jacobian) :-
set_counter(rosenbrock_step,0),
FinalTime is InitialTime + Duration,
StepSizeMax is Duration/20,
(
......@@ -398,10 +394,8 @@ rosenbrock(Equations,InitialTime,Duration,Jacobian) :-
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),
get_last_row(Time, State_full),
debug(rosenbrock, "entering step for ~p", [Time]),
length(Equations, NumberTrueVariable), % Cut the function values at the end of State_full
separate_list(State_full, NumberTrueVariable, CurrentState, _FunctionState),
......@@ -416,7 +410,7 @@ rosenbrock(Equations,InitialTime,Duration,Jacobian) :-
nb_getval(events_number,NumberOfEvents),
NumberOfEvents > 0
->
euler(Equations,Time,CurrentState,CurrentTable)
euler(Equations,Time,CurrentState,_CurrentTable)
;
(
nb_current(functions_list,_)
......@@ -424,13 +418,17 @@ rosenbrock(Equations,InitialTime,Duration,Jacobian) :-
update_functions_state(Next_state,Time2)
;
NewLine =.. [row,Time2|Next_state],
append(CurrentTable,[NewLine],NewTable),
nb_setval(numerical_table,NewTable)
nb_setval(last_row, NewLine),
assertz(saved_row([Time2|Next_state]))
),
nb_setval(last_time,Time2)
),
nb_getval(last_time,LastTime),
LastTime >= FinalTime,!.
LastTime >= FinalTime,!,
%Collect the data for numerical table
findall(Row, (saved_row(Data), Row =.. [row|Data]), RowList),
retractall(saved_row(Data)),
nb_setval(numerical_table, RowList).
%! rosenbrock_step(+Equations, +Jacobian, +CurrentState, +Parameters, +MaxStep, -NextState)
......
......@@ -42,7 +42,7 @@ test(simple_integration, [cleanup(command(clear_model))]) :-
command(present(a,1.0)),
command(numerical_simulation(method:rsbk)),
get_table_data(D),
check_integration_1(D, 1e-3).
check_integration_1(D, 1e-5).
test(cosinus_integration, [cleanup(command(clear_model))]) :-
command(1.0*a for "_" => b),
......@@ -50,6 +50,6 @@ test(cosinus_integration, [cleanup(command(clear_model))]) :-
command(present(a,1.0)),
command(numerical_simulation(method:rsbk)),
get_table_data(D),
check_integration_2(D, 1e-3).
check_integration_2(D, 1e-5).
:- end_tests(rosenbrock).
Markdown is supported
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