Commit 9faf5fe7 authored by HEMERY Mathieu's avatar HEMERY Mathieu
Browse files

Implement a dichotomic event handling and plug rosenbrock with the ODE_solver options

parent 0cecac1a
......@@ -7,6 +7,7 @@ option(
error_epsilon_relative: 1e-6,
initial_step_size: 1e-6,
maximum_step_size: 1e-2,
minimum_step_size: 1e-4,
precision: 6, %5
nusmv_initial_states: all,
nusmv_counter_example: no,
......
......@@ -143,6 +143,10 @@ numerical_simulation :-
maximum_step_size, number, _MaximumStepSize,
'maximum step size for the numerical solver, as a fraction of the total time'
),
option(
minimum_step_size, number, _MinimumStepSize,
'minimum step size, as a fraction of the total time (used to trigger events)'
),
option(
precision, number, _Precision,
'precision for the numerical solver'
......@@ -273,6 +277,7 @@ prepare_numerical_simulation_options(Options) :-
get_option(error_epsilon_relative, ErrorEpsilonRelative),
get_option(initial_step_size, InitialStepSize),
get_option(maximum_step_size, MaximumStepSize),
get_option(minimum_step_size, MinimumStepSize),
get_option(precision, Precision),
get_option(filter, Filter),
set_counter(conditional_events,0),
......@@ -311,6 +316,7 @@ prepare_numerical_simulation_options(Options) :-
error_epsilon_relative: ErrorEpsilonRelative,
initial_step_size: InitialStepSize,
maximum_step_size: MaximumStepSize,
minimum_step_size: MinimumStepSize,
precision: Precision,
filter: Filter,
time_initial: 0,
......
......@@ -15,11 +15,11 @@
% Initialize coefficients for the rosenbrock method, the errors and the minimal step size
% for events handling
rosenbrock_init :-
nb_setval(not_zero_error,0.0001),
rosenbrock_init(Epsilon_abs) :-
nb_setval(not_zero_error, Epsilon_abs),
nb_setval(k_method,stiff),
nb_setval(event_step_size,0.01),
nb_setval(step_doubling_error,1e-6), % Control the method's precision
nb_setval(event_step_size, 1e-3),
nb_setval(step_doubling_error, Epsilon_abs), % 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,
GAM is 1.0/2.0, A21=2.0, A31 is 48.0/25.0, A32 is 6.0/25.0, C21 = -8.0,
......@@ -68,10 +68,11 @@ ode_solver([
initial_parameter_values: InitialParameters,
events: Events,
method: _,
error_epsilon_absolute: _,
error_epsilon_absolute: Epsilon_abs,
error_epsilon_relative: _,
initial_step_size: _,
maximum_step_size: _,
maximum_step_size: MaxStSz,
minimum_step_size: MinStSz,
precision: _,
filter: _,
time_initial: InitialTime,
......@@ -116,8 +117,8 @@ ode_solver([
assertz(saved_row([InitialTime|FullState])),
TrueInitialTime is InitialTime
),
rosenbrock_init,
rosenbrock(Equations,TrueInitialTime,Duration,Jacobian),
rosenbrock_init(Epsilon_abs),
rosenbrock(Equations,TrueInitialTime,Duration,MaxStSz,MinStSz,Jacobian),
retractall(k_fail).
......@@ -355,9 +356,11 @@ get_last_row(Time,CurrentState):-
% The current state of the simulation if obtained in the final row of the global
% variable numerical_table
rosenbrock(Equations,InitialTime,Duration,Jacobian) :-
rosenbrock(Equations, InitialTime, Duration, MaxStSz, MinStSz, Jacobian) :-
FinalTime is InitialTime + Duration,
StepSizeMax is Duration/20,
StepSizeMax is Duration*MaxStSz,
StepSizeMin is Duration*MinStSz,
nb_setval(event_step_size, StepSizeMin),
(
StepSizeMax =:= 0.0,
!,
......@@ -393,20 +396,13 @@ rosenbrock(Equations,InitialTime,Duration,Jacobian) :-
FinalStep is (FinalTime - Time),
MaxStep is min(EventStep, FinalStep),
% Here is the true rosenbrock step
% format("Time:~w,\tMaxStep:~w~n",[Time, MaxStep]),
rosenbrock_step(Equations, Jacobian, Time, CurrentState, Parameters, MaxStep, Next_state),
nb_getval(hdid,NewStepSize),
Time2 is Time+NewStepSize,
( % Store the point when not handling event
nb_getval(event_handling, false)
->
add_functions_values(Next_state, Time2, Next_state_full),
NewLine =.. [row,Time2|Next_state_full],
nb_setval(last_row, NewLine),
assertz(saved_row([Time2|Next_state_full]))
;
true
),
add_functions_values(Next_state, Time2, Next_state_full),
NewLine =.. [row,Time2|Next_state_full],
nb_setval(last_row, NewLine),
assertz(saved_row([Time2|Next_state_full])),
nb_setval(last_time,Time2),
Time2 >= FinalTime,!,
%Collect the data for numerical table
......@@ -598,14 +594,10 @@ rosenbrock_step(Equations,Jacobian,Time,CurrentState,Parameters,MaxStep,NewCurre
%%% EVENTS HANDLING PART %%%
(
NewTime is Time+StepSize,
(
nb_getval(event_handling, true)
;
check_events(NewCurrentState, NewTime),
check_conditional_events(NewCurrentState, NewTime),
nb_getval(events_number,NumberOfEvents),
NumberOfEvents > 0
)
check_events(NewCurrentState, NewTime),
check_conditional_events(NewCurrentState, NewTime),
nb_getval(events_number,NumberOfEvents),
NumberOfEvents > 0
->
nb_setval(event_handling, true),
(
......@@ -623,7 +615,15 @@ rosenbrock_step(Equations,Jacobian,Time,CurrentState,Parameters,MaxStep,NewCurre
fail
)
;
true
(
nb_getval(event_handling, true)
->
nb_getval(event_step_size, ESS),
StepSizeRed is max(SHRINK*StepSize, ESS),
nb_setval(h, StepSizeRed)
;
true
)
),
!,
true.
......
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