grt: add analog_solver (work in progress)

4 files changed, 197 insertions, 9 deletions

diff --git a/src/grt/grt-analog_solver.adb b/src/grt/grt-analog_solver.adb
new file mode 100644
index 000000000..c4d2c3b39
--- /dev/null
+++ b/src/grt/grt-analog_solver.adb
@@ -0,0 +1,84 @@
+--  GHDL Run Time (GRT) -  analog solver for sundials.
+--  Copyright (C) 2022 Tristan Gingold
+--
+--  This program is free software: you can redistribute it and/or modify
+--  it under the terms of the GNU General Public License as published by
+--  the Free Software Foundation, either version 2 of the License, or
+--  (at your option) any later version.
+--
+--  This program is distributed in the hope that it will be useful,
+--  but WITHOUT ANY WARRANTY; without even the implied warranty of
+--  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+--  GNU General Public License for more details.
+--
+--  You should have received a copy of the GNU General Public License
+--  along with this program.  If not, see <gnu.org/licenses>.
+--
+--  As a special exception, if other files instantiate generics from this
+--  unit, or you link this unit with other files to produce an executable,
+--  this unit does not by itself cause the resulting executable to be
+--  covered by the GNU General Public License. This exception does not
+--  however invalidate any other reasons why the executable file might be
+--  covered by the GNU Public License.
+
+with Grt.Errors; use Grt.Errors;
+
+package body Grt.Analog_Solver is
+   function Sundials_Init (Sz : Ghdl_I32) return Ghdl_I32;
+   pragma Import (C, Sundials_Init, "grt_sundials_init");
+
+   function Sundials_Start return Ghdl_I32;
+   pragma Import (C, Sundials_Start, "grt_sundials_start");
+
+   procedure Sundials_Solve (T0 : Ghdl_F64; Tn : in out Ghdl_F64;
+                                            Res : out Integer);
+   pragma Import (C, Sundials_Solve, "grt_sundials_solve");
+
+   procedure Init (Size : Ghdl_I32)
+   is
+      Res : Ghdl_I32;
+   begin
+      Res := Sundials_Init (Size);
+      if Res < 0 then
+         Internal_Error ("sundials initialization failure");
+      end if;
+   end Init;
+
+   function Sundials_Get_Yy_Vec return F64_C_Arr_Ptr;
+   pragma Import (C, Sundials_Get_Yy_Vec, "grt_sundials_get_yy_vec");
+
+   function Get_Init_Val_Ptr return F64_C_Arr_Ptr is
+   begin
+      return Sundials_Get_Yy_Vec;
+   end Get_Init_Val_Ptr;
+
+   function Sundials_Get_Yp_Vec return F64_C_Arr_Ptr;
+   pragma Import (C, Sundials_Get_Yp_Vec, "grt_sundials_get_yp_vec");
+
+   function Get_Init_Der_Ptr return F64_C_Arr_Ptr is
+   begin
+      return Sundials_Get_Yp_Vec;
+   end Get_Init_Der_Ptr;
+
+   procedure Start is
+   begin
+      if Sundials_Start /= 0 then
+         Internal_Error ("sundials start");
+      end if;
+   end Start;
+
+   procedure Set_Root_Size (Size : Ghdl_I32) is
+   begin
+      Internal_Error ("sundials set_root_size");
+   end Set_Root_Size;
+
+   procedure Solve (T : Ghdl_F64; Tn : in out Ghdl_F64; Res : out Integer) is
+   begin
+      Sundials_Solve (T, Tn, Res);
+   end Solve;
+
+   procedure Finish is
+   begin
+      Internal_Error ("sundials finish");
+   end Finish;
+end Grt.Analog_Solver;
diff --git a/src/grt/grt-analog_solver.ads b/src/grt/grt-analog_solver.ads
new file mode 100644
index 000000000..84490f16f
--- /dev/null
+++ b/src/grt/grt-analog_solver.ads
@@ -0,0 +1,69 @@
+--  GHDL Run Time (GRT) -  analog solver.
+--  Copyright (C) 2022 Tristan Gingold
+--
+--  This program is free software: you can redistribute it and/or modify
+--  it under the terms of the GNU General Public License as published by
+--  the Free Software Foundation, either version 2 of the License, or
+--  (at your option) any later version.
+--
+--  This program is distributed in the hope that it will be useful,
+--  but WITHOUT ANY WARRANTY; without even the implied warranty of
+--  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+--  GNU General Public License for more details.
+--
+--  You should have received a copy of the GNU General Public License
+--  along with this program.  If not, see <gnu.org/licenses>.
+--
+--  As a special exception, if other files instantiate generics from this
+--  unit, or you link this unit with other files to produce an executable,
+--  this unit does not by itself cause the resulting executable to be
+--  covered by the GNU General Public License. This exception does not
+--  however invalidate any other reasons why the executable file might be
+--  covered by the GNU Public License.
+
+with Grt.Types; use Grt.Types;
+
+package Grt.Analog_Solver is
+   type F64_Array is array (Natural range <>) of Ghdl_F64;
+   subtype F64_Fat_Array is F64_Array (Natural);
+
+   --  A pointer to an F64 array.
+   type F64_C_Arr_Ptr is access all F64_Fat_Array;
+   pragma Convention (C, F64_C_Arr_Ptr);
+
+   --  Initialize the analog solver, SIZE is the number of scalar quantities.
+   procedure Init (Size : Ghdl_I32);
+
+   --  Return the address of the initial values vector and derivative vector.
+   function Get_Init_Val_Ptr return F64_C_Arr_Ptr;
+   function Get_Init_Der_Ptr return F64_C_Arr_Ptr;
+
+   --  Finish initialization.
+   procedure Start;
+
+   procedure Residues (T : Ghdl_F64;
+                       Y : F64_C_Arr_Ptr;
+                       Yp : F64_C_Arr_Ptr;
+                       Res : F64_C_Arr_Ptr);
+   pragma Import (C, Residues, "grt__analog_solver__residues");
+
+   procedure Set_Root_Size (Size : Ghdl_I32);
+
+   procedure Roots (T : Ghdl_F64;
+                    Y : F64_C_Arr_Ptr;
+                    Yp: F64_C_Arr_Ptr;
+                    Res : F64_C_Arr_Ptr);
+   pragma Import (C, Roots, "grt__analog_solver__roots");
+
+   procedure Set_Values (Y : F64_C_Arr_Ptr;
+                         Yp: F64_C_Arr_Ptr);
+   pragma Import (C, Set_Values, "grt__analog_solver__set_values");
+
+   --  Return status:
+   --  0: Ok, Tn reached
+   --  1: Stopped due to zero crossing
+   --  2: failure
+   procedure Solve (T : Ghdl_F64; Tn : in out Ghdl_F64; Res : out Integer);
+
+   procedure Finish;
+end Grt.Analog_Solver;
diff --git a/src/grt/grt-processes.adb b/src/grt/grt-processes.adb
index 9eba3b677..6133343d8 100644
--- a/src/grt/grt-processes.adb
+++ b/src/grt/grt-processes.adb
@@ -22,8 +22,6 @@
 --  covered by the GNU Public License.
 with Grt.Table;
 with Ada.Unchecked_Deallocation;
-with System.Storage_Elements; --  Work around GNAT bug.
-pragma Unreferenced (System.Storage_Elements);
 with Grt.Disp;
 with Grt.Astdio;
 with Grt.Astdio.Vhdl; use Grt.Astdio.Vhdl;
@@ -39,6 +37,7 @@ with Grt.Stats;
 with Grt.Threads; use Grt.Threads;
 pragma Elaborate_All (Grt.Table);
 with Grt.Stdio;
+with Grt.Analog_Solver;
 
 package body Grt.Processes is
    Last_Time : constant Std_Time := Std_Time'Last;
@@ -880,12 +879,20 @@ package body Grt.Processes is
       --  - The time of the next simulation cycle (which in this case is the
       --    first simulation cycle), Tn, is calculated according to the rules
       --    of step f of the simulation cycle, below.
-      Next_Time := Compute_Next_Time;
-      if Next_Time /= 0 then
-         if Has_Callbacks (Hooks.Cb_Last_Known_Delta) then
-            Call_Callbacks (Hooks.Cb_Last_Known_Delta);
-            Flush_Active_Chain;
-            Next_Time := Compute_Next_Time;
+
+      --  LRM 1076.1-2007
+      --  - The time of the next simulation cycle (which in this case is the
+      --    first simulation cycle), Tn, is set to 0.0
+      if Flag_AMS then
+         Next_Time := 0;
+      else
+         Next_Time := Compute_Next_Time;
+         if Next_Time /= 0 then
+            if Has_Callbacks (Hooks.Cb_Last_Known_Delta) then
+               Call_Callbacks (Hooks.Cb_Last_Known_Delta);
+               Flush_Active_Chain;
+               Next_Time := Compute_Next_Time;
+            end if;
          end if;
       end if;
 
@@ -896,12 +903,26 @@ package body Grt.Processes is
    --  Launch a simulation cycle.
    function Simulation_Cycle return Integer
    is
+      use Grt.Options;
       Tn : Std_Time;
+      Tn_AMS : Ghdl_F64;
       Status : Integer;
    begin
       --  LRM08 14.7.5.3 Simulation cycle (ex LRM93 12.6.4)
       --  A simulation cycle consists of the following steps:
       --
+
+      --  LRM 1076.1-2007 12.6.4 Simulation cycle
+      --  a) The analog solver is executed
+      if Flag_AMS and Next_Time > Current_Time then
+         Current_Time_AMS := Ghdl_F64 (Current_Time) * Time_Phys_To_Real;
+         Tn_AMS := Ghdl_F64 (Next_Time) * Time_Phys_To_Real;
+         Grt.Analog_Solver.Solve (Current_Time_AMS, Tn_AMS, Status);
+         if Status /= 0 then
+            Internal_Error ("simulation_cycle - analog_solver");
+         end if;
+      end if;
+
       --  a) The current time, Tc is set equal to Tn.  Simulation is complete
       --     when Tn = TIME'HIGH and there are no active drivers or process
       --     resumptions at Tn.
@@ -1013,7 +1034,11 @@ package body Grt.Processes is
       if Options.Flag_Stats then
          Stats.Start_Next_Time;
       end if;
-      Tn := Compute_Next_Time;
+      if Flag_AMS and Break_Flag then
+         Tn := Current_Time;
+      else
+         Tn := Compute_Next_Time;
+      end if;
 
       --  h) If the next simulation cycle will be a delta cycle, the remainder
       --     of the step is skipped. Otherwise the following actions occur
diff --git a/src/grt/grt-processes.ads b/src/grt/grt-processes.ads
index 80df93641..da820a8b0 100644
--- a/src/grt/grt-processes.ads
+++ b/src/grt/grt-processes.ads
@@ -70,6 +70,16 @@ package Grt.Processes is
    --  Number of non-delta cycles.
    Nbr_Cycles : Ghdl_I64;
 
+   --  If True, do AMS simulation
+   Flag_AMS : Boolean := False;
+
+   --  The break flag of AMS simulation
+   --
+   --  LRM 1076.1-2017 14.7.5.2 Initialization
+   --  At the beginning of initialization, [...] and the break flag is assumed
+   --  to be cleared.
+   Break_Flag : Boolean := False;
+
    type Process_Type is private;
    --  type Process_Acc is access all Process_Type;