Improved stability of algebraic loop solutions

src/OMSimulatorLib/SystemSC.h

@@ -65,14 +65,18 @@ namespace oms
 
     oms_status_enu_t updateInputs(DirectedGraph& graph);
 
+    oms_status_enu_t saveInputs(DirectedGraph& graph, std::vector<double> &saved_inputs);
+
+    oms_status_enu_t restoreInputs(DirectedGraph& graph, const std::vector<double>& saved_outputs);
+
     std::string getSolverName() const;
     oms_status_enu_t setSolverMethod(std::string);
 
     oms_status_enu_t setSolver(oms_solver_enu_t solver) {if (solver > oms_solver_sc_min && solver < oms_solver_sc_max) {solverMethod=solver; return oms_status_ok;} return oms_status_error;}
 
   private:
-    oms_status_enu_t doStepEuler();
-    oms_status_enu_t doStepCVODE();
+    oms_status_enu_t doStepEuler(double stopTime);
+    oms_status_enu_t doStepCVODE(double stopTime);
 
   protected:
     SystemSC(const ComRef& cref, Model* parentModel, System* parentSystem);
@@ -95,6 +99,8 @@ namespace oms
     std::vector<double*> event_indicators;
     std::vector<double*> event_indicators_prev;
 
+    std::vector<double> initial_guess;  // Initial guess for all algebraic loops
+
     struct SolverDataEuler_t
     {
     };

testsuite/reference-fmus/Dahlquist.lua

@@ -58,7 +58,7 @@ end
 -- info:    maximum step size for 'model.root': 0.001000
 -- info:    Result file: Dahlquist-me.mat (bufferSize=1)
 -- info:    Final Statistics for 'model.root':
---          NumSteps = 10001 NumRhsEvals  = 10002 NumLinSolvSetups = 501
+--          NumSteps = 10001 NumRhsEvals  = 10002 NumLinSolvSetups = 502
 --          NumNonlinSolvIters = 10001 NumNonlinSolvConvFails = 0 NumErrTestFails = 0
 -- signal x is equal
 -- endResult

testsuite/simulation/DualMassOscillator.lua

@@ -50,9 +50,9 @@ oms_delete("DualMassOscillator")
 -- info:      system1.x1: 0.0
 -- info:      system2.x2: 0.5
 -- info:    Simulation
--- info:      system1.x1: 0.051037644866066
--- info:      system2.x2: 0.031639500249976
+-- info:      system1.x1: 0.051037644963887
+-- info:      system2.x2: 0.031639500306724
 -- info:    Final Statistics for 'DualMassOscillator.root':
---          NumSteps = 10007 NumRhsEvals  = 10008 NumLinSolvSetups = 507
+--          NumSteps = 10007 NumRhsEvals  = 10008 NumLinSolvSetups = 508
 --          NumNonlinSolvIters = 10007 NumNonlinSolvConvFails = 0 NumErrTestFails = 0
 -- endResult

testsuite/simulation/EventTest.lua

@@ -79,16 +79,15 @@ readFile("EventTest_lua.csv")
 -- 0, 0.3, -1
 -- 0.30000001, -1.00000002723e-08, -1
 -- 0.30000001, 0.3, -1
--- 0.60000002, -1.00000154823e-08, -1
+-- 0.60000002, -1.00000004943e-08, -1
 -- 0.60000002, 0.3, -1
--- 0.90000003, -1.00000214776e-08, -1
+-- 0.90000003, -1.00000004943e-08, -1
 -- 0.90000003, 0.3, -1
--- 1, 0.20000003, -1
--- 1.20000004, -1.0000003936e-08, -1
+-- 1.20000004, -1.00000006054e-08, -1
 -- 1.20000004, 0.3, -1
--- 1.50000005, -1.00000143721e-08, -1
+-- 1.50000005, -1.00000290271e-08, -1
 -- 1.50000005, 0.3, -1
--- 1.80000006, -1.00000390191e-08, -1
+-- 1.80000006, -1.0000053674e-08, -1
 -- 1.80000006, 0.3, -1
 -- 2, 0.10000006, -1
 -- 2, 0.10000006, -1

testsuite/simulation/Makefile

@@ -50,6 +50,7 @@ importPartialSnapshot.lua \
 importStartValues.lua \
 importStartValues1.lua \
 importSuppressUnitConversion.lua \
+integratorTest.lua \
 listVariants1.lua \
 multipleConnections.lua \
 multipleInstance.lua \

testsuite/simulation/integratorTest.lua

@@ -0,0 +1,79 @@
+-- status: correct
+-- teardown_command: rm -rf integratorTest_lua/
+-- linux: yes
+-- ucrt64: yes
+-- win: yes
+-- mac: no
+
+function readFile(filename)
+  local f = assert(io.open(filename, "r"))
+  local content = f:read("*all")
+  print(content)
+  f:close()
+end
+
+oms_setCommandLineOption("--suppressPath=true")
+oms_setTempDirectory("./integratorTest_lua/")
+oms_setWorkingDirectory("./integratorTest_lua/")
+
+oms_newModel("IntegratorTest")
+oms_addSystem("IntegratorTest.root", oms_system_sc)
+oms_addSubModel("IntegratorTest.root.integrator", "../../resources/Modelica.Blocks.Continuous.Integrator.fmu")
+
+-- simulation settings
+oms_setSolver("IntegratorTest.root", oms_solver_sc_cvode)
+oms_setResultFile("IntegratorTest", "integratorTest_lua.csv")
+oms_setStopTime("IntegratorTest", 10.0)
+oms_setVariableStepSize("IntegratorTest.root", 1.0, 1e-8, 1.0)
+
+oms_instantiate("IntegratorTest")
+oms_initialize("IntegratorTest")
+
+oms_setReal("IntegratorTest.root.integrator.u", 0.0)
+print(string.format("t = 0.0: y = %g", oms_getReal("IntegratorTest.root.integrator.y")))
+oms_stepUntil("IntegratorTest", 2.5)
+oms_setReal("IntegratorTest.root.integrator.u", 10.0)
+
+print(string.format("t = 2.5: y = %g", oms_getReal("IntegratorTest.root.integrator.y")))
+oms_stepUntil("IntegratorTest", 5.0)
+print(string.format("t = 5.0: y = %g", oms_getReal("IntegratorTest.root.integrator.y")))
+oms_stepUntil("IntegratorTest", 10.0)
+print(string.format("t = 10.0: y = %g", oms_getReal("IntegratorTest.root.integrator.y")))
+
+oms_delete("IntegratorTest")
+
+readFile("integratorTest_lua.csv")
+
+-- Result:
+-- info:    maximum step size for 'IntegratorTest.root': 1.000000
+-- info:    Result file: integratorTest_lua.csv (bufferSize=1)
+-- t = 0.0: y = 0
+-- t = 2.5: y = 0
+-- t = 5.0: y = 25
+-- t = 10.0: y = 75
+-- info:    Final Statistics for 'IntegratorTest.root':
+--          NumSteps = 16 NumRhsEvals  = 26 NumLinSolvSetups = 13
+--          NumNonlinSolvIters = 24 NumNonlinSolvConvFails = 0 NumErrTestFails = 4
+-- time,IntegratorTest.root.integrator._D_outputAlias_y,IntegratorTest.root.integrator.der(_D_outputAlias_y),IntegratorTest.root.integrator.local_set,IntegratorTest.root.integrator.u,IntegratorTest.root.integrator.y,IntegratorTest.root.integrator.initType,IntegratorTest.root.integrator.local_reset,IntegratorTest.root.integrator.use_reset,IntegratorTest.root.integrator.use_set
+-- 0, 0, 0, 0, 0, 0, 3, 0, 0, 0
+-- 1, 0, 0, 0, 0, 0, 3, 0, 0, 0
+-- 2, 0, 0, 0, 0, 0, 3, 0, 0, 0
+-- 2.5, 0, 0, 0, 0, 0, 3, 0, 0, 0
+-- 2.5, 0, 0, 0, 0, 0, 3, 0, 0, 0
+-- 2.5001, 0.001, 10, 0, 10, 0.001, 3, 0, 0, 0
+-- 2.5002, 0.002, 10, 0, 10, 0.002, 3, 0, 0, 0
+-- 2.5012, 0.012, 10, 0, 10, 0.012, 3, 0, 0, 0
+-- 2.5112, 0.112, 10, 0, 10, 0.112, 3, 0, 0, 0
+-- 2.6112, 1.112, 10, 0, 10, 1.112, 3, 0, 0, 0
+-- 3.6112, 11.112, 10, 0, 10, 11.112, 3, 0, 0, 0
+-- 4.6112, 21.112, 10, 0, 10, 21.112, 3, 0, 0, 0
+-- 5, 25, 10, 0, 10, 25, 3, 0, 0, 0
+-- 5, 25, 10, 0, 10, 25, 3, 0, 0, 0
+-- 6, 35, 10, 0, 10, 35, 3, 0, 0, 0
+-- 7, 45, 10, 0, 10, 45, 3, 0, 0, 0
+-- 8, 55, 10, 0, 10, 55, 3, 0, 0, 0
+-- 9, 65, 10, 0, 10, 65, 3, 0, 0, 0
+-- 10, 75, 10, 0, 10, 75, 3, 0, 0, 0
+-- 10, 75, 10, 0, 10, 75, 3, 0, 0, 0
+--
+-- endResult