Dynamic tolerance

src/OMSimulatorLib/AlgLoop.cpp

@@ -29,14 +29,24 @@
  *
  */
 
+#include <sstream>
+#include <cmath>
+#include <algorithm>
+
 #include "AlgLoop.h"
 
 #include "Component.h"
 #include "DirectedGraph.h"
 #include "Flags.h"
 #include "System.h"
 
-#include <sstream>
+#ifdef min
+#undef min
+#endif
+
+#ifdef max
+#undef max
+#endif
 
 /**
  * @brief Check flag returned by KINSOL function and log error
@@ -292,7 +302,7 @@ oms::KinsolSolver::~KinsolSolver()
  * @param absoluteTolerance     Tolerance used for solving the loop
  * @return oms::KinsolSolver*   Retruns pointer to KinsolSolver object
  */
-oms::KinsolSolver* oms::KinsolSolver::NewKinsolSolver(const int algLoopNum, const unsigned int size, double absoluteTolerance, const bool useDirectionalDerivative)
+oms::KinsolSolver* oms::KinsolSolver::NewKinsolSolver(const int algLoopNum, const unsigned int size, double absoluteTolerance, double relativeTolerance, const bool useDirectionalDerivative)
 {
   int flag;
   int printLevel;
@@ -301,6 +311,7 @@ oms::KinsolSolver* oms::KinsolSolver::NewKinsolSolver(const int algLoopNum, cons
   logDebug("Create new KinsolSolver object for algebraic loop number " + std::to_string(algLoopNum));
 
   kinsolSolver->size = size;
+  kinsolSolver->firstSolution = true;
 
   /* Allocate memory */
   kinsolSolver->initialGuess = N_VNew_Serial(kinsolSolver->size);
@@ -364,10 +375,12 @@ oms::KinsolSolver* oms::KinsolSolver::NewKinsolSolver(const int algLoopNum, cons
   if (!checkFlag(flag, "KINSetJacFn")) return NULL;
 
   /* Set function-norm stopping tolerance */
-  kinsolSolver->fnormtol = absoluteTolerance;
+  kinsolSolver->fnormtol = 0.01 * absoluteTolerance;
   flag = KINSetFuncNormTol(kinsolSolver->kinsolMemory, kinsolSolver->fnormtol);
   if (!checkFlag(flag, "KINSetFuncNormTol")) return NULL;
 
+  kinsolSolver->freltol = relativeTolerance * 0.01;
+
   /* Set scaled-step stopping tolerance */
   flag = KINSetScaledStepTol(kinsolSolver->kinsolMemory, 0.0);
   if (!checkFlag(flag, "KINSetScaledStepTol")) return NULL;
@@ -400,7 +413,7 @@ oms::KinsolSolver* oms::KinsolSolver::NewKinsolSolver(const int algLoopNum, cons
  *                            `oms_status_warning` if solving was computed, but solution is not within tolerance
  *                            and `oms_status_error` if an error occured.
  */
-oms_status_enu_t oms::KinsolSolver::kinsolSolve(System& syst, DirectedGraph& graph)
+oms_status_enu_t oms::KinsolSolver::kinsolSolve(System& syst, DirectedGraph& graph, double tolerance /*= 0.0*/)
 {
   /* Update user data */
   KINSOL_USER_DATA* kinsolUserData = (KINSOL_USER_DATA*) user_data;
@@ -413,7 +426,7 @@ oms_status_enu_t oms::KinsolSolver::kinsolSolve(System& syst, DirectedGraph& gra
   int flag;
   double fNormValue;
 
-  logDebug("Solving system " + std::to_string(kinsolUserData->algLoopNumber));
+  double tol = tolerance != 0.0 ? tolerance : fnormtol;
 
   if (SCC.connections.size() != size)
   {
@@ -432,9 +445,29 @@ oms_status_enu_t oms::KinsolSolver::kinsolSolve(System& syst, DirectedGraph& gra
     }
   }
 
+  /* Apply relative tolerance to magnitude of initial guess */
+  fNormValue = std::sqrt(N_VDotProd_Serial(initialGuess, initialGuess));
+  tol = std::max(tol, fNormValue * freltol * tolerance / fnormtol);
+
+  /* Check residual for initial guess */
+  nlsKinsolResiduals(initialGuess, fTmp, user_data);
+  fNormValue = std::sqrt(N_VDotProd_Serial(fTmp, fTmp));
+  if (fNormValue > 0.0)
+    tol = std::min(fNormValue, tol); // We already know this is achievable, but let's get in at least one iteration
+  else
+    return oms_status_ok;
+
+  if (!firstSolution) {
+    // Predict a better initial guess
+    SUNLinSolSolve(linSol, J, y, fTmp, tol);
+    N_VLinearSum(1.0, initialGuess, -1.0, y, initialGuess);
+  }
+
   /* u and f scaling */
   // TODO: Add scaling that is not only constant ones
 
+  KINSetFuncNormTol(kinsolMemory, tol);
+
   /* Solve algebraic loop with KINSol() */
   flag = KINSol(kinsolMemory,   /* KINSol memory block */
                 initialGuess,   /* initial guess on input; solution vector */
@@ -444,15 +477,19 @@ oms_status_enu_t oms::KinsolSolver::kinsolSolve(System& syst, DirectedGraph& gra
   if (!checkFlag(flag, "KINSol")) return oms_status_error;
 
   /* Check solution */
-  flag = nlsKinsolResiduals(initialGuess, fTmp, user_data);
-  fNormValue = N_VWL2Norm(fTmp, fTmp);
-  if ( fNormValue > fnormtol )
+  KINGetFuncNorm(kinsolMemory, &fNormValue);
+  if ( fNormValue > tol )
   {
     logWarning("Solution of algebraic loop " + std::to_string(((KINSOL_USER_DATA *)user_data)->algLoopNumber) + "not within precission given by fnormtol: " + std::to_string(fnormtol));
     logDebug("2-norm of residual of solution: " + std::to_string(fNormValue));
     return oms_status_warning;
   }
 
+  if (flag == KIN_SUCCESS) {
+    firstSolution = false;
+    KINSetNoInitSetup(kinsolMemory, SUNTRUE);
+  }
+
   logDebug("Solved system " + std::to_string(kinsolUserData->algLoopNumber) + " successfully");
   return oms_status_ok;
 }
@@ -467,7 +504,7 @@ oms_status_enu_t oms::KinsolSolver::kinsolSolve(System& syst, DirectedGraph& gra
  * @param scc     Strong Connected Componten, a vector of connected
  * @param number
  */
-oms::AlgLoop::AlgLoop(oms_alg_solver_enu_t method, double absTol, scc_t scc, const int number, const bool useDirectionalDerivative): absoluteTolerance(absTol), SCC(scc), systNumber(number)
+oms::AlgLoop::AlgLoop(oms_alg_solver_enu_t method, double absTol, double relTol, scc_t scc, const int number, const bool useDirectionalDerivative): absoluteTolerance(absTol), SCC(scc), systNumber(number)
 {
   switch (method)
   {
@@ -482,7 +519,7 @@ oms::AlgLoop::AlgLoop(oms_alg_solver_enu_t method, double absTol, scc_t scc, con
 
   if (method == oms_alg_solver_kinsol)
   {
-    kinsolData = KinsolSolver::NewKinsolSolver(systNumber, SCC.connections.size(), absoluteTolerance, useDirectionalDerivative);
+    kinsolData = KinsolSolver::NewKinsolSolver(systNumber, SCC.connections.size(), absoluteTolerance, relTol, useDirectionalDerivative);
     if (kinsolData==NULL)
     {
       logError("NewKinsolSolver() failed. Aborting!");
@@ -501,17 +538,17 @@ oms::AlgLoop::AlgLoop(oms_alg_solver_enu_t method, double absTol, scc_t scc, con
  * @param graph               Reference to directed graph
  * @return oms_status_enu_t   Returns `oms_status_ok` on success
  */
-oms_status_enu_t oms::AlgLoop::solveAlgLoop(System& syst, DirectedGraph& graph)
+oms_status_enu_t oms::AlgLoop::solveAlgLoop(System& syst, DirectedGraph& graph, double tolerance)
 {
   logDebug("Solving algebraic loop formed by connections\n" + dumpLoopVars(graph));
   logDebug("Using solver " + getAlgSolverName());
 
   switch (algSolverMethod)
   {
   case oms_alg_solver_fixedpoint:
-    return fixPointIteration(syst, graph);
+    return fixPointIteration(syst, graph, tolerance);
   case oms_alg_solver_kinsol:
-    return kinsolData->kinsolSolve(syst, graph);
+    return kinsolData->kinsolSolve(syst, graph, tolerance);
   default:
     logError("Invalid algebraic solver method!");
     return oms_status_error;
@@ -525,14 +562,17 @@ oms_status_enu_t oms::AlgLoop::solveAlgLoop(System& syst, DirectedGraph& graph)
  * @param graph
  * @return oms_status_enu_t
  */
-oms_status_enu_t oms::AlgLoop::fixPointIteration(System& syst, DirectedGraph& graph)
+oms_status_enu_t oms::AlgLoop::fixPointIteration(System& syst, DirectedGraph& graph, double tolerance)
 {
   const int size = SCC.connections.size();
   const int maxIterations = Flags::MaxLoopIteration();
   int it=0;
   double maxRes;
   double *res = new double[size]();
 
+  if (tolerance == 0.0 || tolerance > absoluteTolerance)
+    tolerance = absoluteTolerance;
+
   do
   {
     std::stringstream ss;
@@ -641,7 +681,7 @@ oms_status_enu_t oms::AlgLoop::fixPointIteration(System& syst, DirectedGraph& gr
       logInfo(ss.str());
     }
 
-  } while(maxRes > absoluteTolerance && it < maxIterations);
+  } while(maxRes > tolerance && it < maxIterations);
 
   delete[] res;
 

src/OMSimulatorLib/AlgLoop.h

@@ -57,12 +57,13 @@ namespace oms
   {
   public:
     ~KinsolSolver();
-    static KinsolSolver* NewKinsolSolver(const int algLoopNum, const unsigned int size, double absoluteTolerance, const bool useDirectionalDerivative);
-    oms_status_enu_t kinsolSolve(System& syst, DirectedGraph& graph);
+    static KinsolSolver* NewKinsolSolver(const int algLoopNum, const unsigned int size, double absoluteTolerance, double relativeTolerance, const bool useDirectionalDerivative);
+    oms_status_enu_t kinsolSolve(System& syst, DirectedGraph& graph, double tolerance = 0.0);
 
   private:
     /* tolerances */
     double fnormtol;        /* function tolerance */
+    double freltol;         /* relative function tolerance */
 
     /* work arrays */
     N_Vector initialGuess;
@@ -80,6 +81,8 @@ namespace oms
     N_Vector y;             /* Template for cloning vectors needed inside linear solver */
     SUNMatrix J;            /* (Non-)Sparse matrix template for cloning matrices needed within linear solver */
 
+    bool firstSolution;
+
     /* member function */
     static int nlsKinsolJac(N_Vector u, N_Vector fu, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2);
     static int nlsKinsolResiduals(N_Vector u, N_Vector fval, void *user_data);
@@ -90,16 +93,16 @@ namespace oms
   class AlgLoop
   {
   public:
-    AlgLoop(oms_alg_solver_enu_t method, double absTol, scc_t SCC, const int systNumber, const bool useDirectionalDerivative);
+    AlgLoop(oms_alg_solver_enu_t method, double absTol, double relTol, scc_t SCC, const int systNumber, const bool useDirectionalDerivative);
 
     scc_t getSCC() {return SCC;}
-    oms_status_enu_t solveAlgLoop(System& syst, DirectedGraph& graph);
+    oms_status_enu_t solveAlgLoop(System& syst, DirectedGraph& graph, double tolerance);
     std::string getAlgSolverName();
     std::string dumpLoopVars(DirectedGraph& graph);
 
   private:
     oms_alg_solver_enu_t algSolverMethod;
-    oms_status_enu_t fixPointIteration(System& syst, DirectedGraph& graph);
+    oms_status_enu_t fixPointIteration(System& syst, DirectedGraph& graph, double tolerance);
 
     KinsolSolver* kinsolData;
 

src/OMSimulatorLib/System.cpp

@@ -3066,7 +3066,7 @@ oms_status_enu_t oms::System::addAlgLoop(scc_t SCC, const int algLoopNum, Direct
     loopsNeedUpdate = false;
   }
 
-  algLoops.push_back( AlgLoop(Flags::AlgLoopSolver(), absoluteTolerance, SCC, algLoopNum, supportsDirectionalDerivatives));
+  algLoops.push_back( AlgLoop(Flags::AlgLoopSolver(), absoluteTolerance, relativeTolerance, SCC, algLoopNum, supportsDirectionalDerivatives));
 
   return oms_status_ok;
 }
@@ -3110,9 +3110,9 @@ oms_status_enu_t oms::System::updateAlgebraicLoops(const std::vector<scc_t>& sor
 }
 
 
-oms_status_enu_t oms::System::solveAlgLoop(DirectedGraph& graph, int loopNumber)
+oms_status_enu_t oms::System::solveAlgLoop(DirectedGraph& graph, int loopNumber, double tolerance)
 {
-  return algLoops[loopNumber].solveAlgLoop(*this, graph);
+  return algLoops[loopNumber].solveAlgLoop(*this, graph, tolerance);
 }
 
 oms_status_enu_t oms::System::rename(const oms::ComRef& newCref)

src/OMSimulatorLib/System.h

@@ -188,7 +188,7 @@ namespace oms
     AlgLoop* getAlgLoop(const int systemNumber);
     oms_status_enu_t addAlgLoop(scc_t SCC, const int algLoopNum, DirectedGraph& graph, bool supportsDirectionalDerivatives);
     oms_status_enu_t updateAlgebraicLoops(const std::vector< scc_t >& sortedConnections, DirectedGraph& graph);
-    oms_status_enu_t solveAlgLoop(DirectedGraph& graph, int loopNumber);
+    oms_status_enu_t solveAlgLoop(DirectedGraph& graph, int loopNumber, double tolerance = 0.0);
 
     bool useThreadPool();
     ctpl::thread_pool& getThreadPool();