[SimCode/C] update jacobian handling (#9200)

* [SimCode/C] update jacobian handling

 - use symbolical jacobian, sparsity pattern and coloring if available
 - if the flag -jacobian= is used, try to use the best possible option and issue a warning if not possible

OMCompiler/Compiler/SimCode/SimCodeUtil.mo

@@ -1440,6 +1440,18 @@ algorithm
   outJac.partitionIndex := index;
 end rewriteJacPartIdx;
 
+public function jacobianColumnsAreEmpty
+  input list<SimCode.JacobianColumn> columns;
+  output Boolean b = true;
+algorithm
+  for col in columns loop
+    if not (listEmpty(col.columnEqns) and listEmpty(col.constantEqns)) then
+      b := false;
+      return;
+    end if;
+  end for;
+end jacobianColumnsAreEmpty;
+
 // =============================================================================
 // section to create SimCode.Equations from BackendDAE.Equation
 //

OMCompiler/Compiler/Template/CodegenC.tpl

@@ -5230,6 +5230,7 @@ match sparsepattern
     {
       TRACE_PUSH
       TRACE_POP
+      jacobian->availability = JACOBIAN_NOT_AVAILABLE;
       return 1;
     }
     >>
@@ -5240,6 +5241,7 @@ match sparsepattern
       let colPtr = genSPCRSPtr(listLength(sparsepattern), sparsepattern, "colPtrIndex")
       let rowIndex = genSPCRSRows(lengthListElements(unzipSecond(sparsepattern)), sparsepattern, "rowIndex")
       let colorString = genSPColors(colorList, "jacobian->sparsePattern->colorCols")
+      let availability = if SimCodeUtil.jacobianColumnsAreEmpty(jacobianColumn) then 'JACOBIAN_ONLY_SPARSITY' else 'JACOBIAN_AVAILABLE'
       let indexColumn = (jacobianColumn |> JAC_COLUMN(numberOfResultVars=n) => '<%n%>';separator="\n")
       let tmpvarsSize = (jacobianColumn |> JAC_COLUMN(columnVars=vars) => listLength(vars);separator="\n")
       let constantEqns = (jacobianColumn |> JAC_COLUMN(constantEqns=constantEqns) =>
@@ -5257,6 +5259,7 @@ match sparsepattern
 
         initAnalyticJacobian(jacobian, <%index_%>, <%indexColumn%>, <%tmpvarsSize%>, <%constantEqns%>, jacobian->sparsePattern);
         jacobian->sparsePattern = allocSparsePattern(<%sizeleadindex%>, <%sp_size_index%>, <%maxColor%>);
+        jacobian->availability = <%availability%>;
 
         /* write lead index of compressed sparse column */
         memcpy(jacobian->sparsePattern->leadindex, colPtrIndex, (<%sizeleadindex%>+1)*sizeof(unsigned int));

OMCompiler/Compiler/Template/SimCodeTV.mo

@@ -1339,6 +1339,11 @@ package SimCodeUtil
     output list<SimCode.SimEqSystem> eqs;
   end selectNLEqSys;
 
+  function jacobianColumnsAreEmpty
+    input list<SimCode.JacobianColumn> columns;
+    output Boolean b ;
+  end jacobianColumnsAreEmpty;
+
 end SimCodeUtil;
 
 package SimCodeFunctionUtil

OMCompiler/SimulationRuntime/c/simulation/solver/dassl.c

@@ -42,6 +42,7 @@
 
 #include "gc/omc_gc.h"
 #include "util/context.h"
+#include "util/jacobian_util.h"
 #include "util/omc_error.h"
 #include "util/parallel_helper.h"
 
@@ -345,24 +346,24 @@ int dassl_initial(DATA* data, threadData_t *threadData,
     dasslData->dasslJacobian = COLOREDNUMJAC;
   }
 
-  /* selects the calculation method of the jacobian */
-  if(dasslData->dasslJacobian == COLOREDNUMJAC ||
-     dasslData->dasslJacobian == COLOREDSYMJAC ||
-     dasslData->dasslJacobian == SYMJAC)
-  {
-    ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A]);
-    if (data->callback->initialAnalyticJacobianA(data, threadData, jacobian))
-    {
-      infoStreamPrint(LOG_STDOUT, 0, "Jacobian or SparsePattern is not generated or failed to initialize! Switch back to normal.");
-      dasslData->dasslJacobian = INTERNALNUMJAC;
-    } else {
-      ANALYTIC_JACOBIAN* jac = &data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A];
-      infoStreamPrint(LOG_SIMULATION, 1, "Initialized colored Jacobian:");
-      infoStreamPrint(LOG_SIMULATION, 0, "columns: %d rows: %d", jac->sizeCols, jac->sizeRows);
-      infoStreamPrint(LOG_SIMULATION, 0, "NNZ:  %d colors: %d", jac->sparsePattern->numberOfNonZeros, jac->sparsePattern->maxColors);
-      messageClose(LOG_SIMULATION);
-    }
+  ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A]);
+  data->callback->initialAnalyticJacobianA(data, threadData, jacobian);
+  if(jacobian->availability == JACOBIAN_AVAILABLE || jacobian->availability == JACOBIAN_ONLY_SPARSITY) {
+    infoStreamPrint(LOG_SIMULATION, 1, "Initialized Jacobian:");
+    infoStreamPrint(LOG_SIMULATION, 0, "columns: %d rows: %d", jacobian->sizeCols, jacobian->sizeRows);
+    infoStreamPrint(LOG_SIMULATION, 0, "NNZ:  %d colors: %d", jacobian->sparsePattern->numberOfNonZeros, jacobian->sparsePattern->maxColors);
+    messageClose(LOG_SIMULATION);
   }
+
+  // Compare user flag to availabe Jacobian methods
+  const char* flagValue;
+  if(omc_flag[FLAG_JACOBIAN]){
+    flagValue = omc_flagValue[FLAG_JACOBIAN];
+  } else {
+    flagValue = NULL;
+  }
+  dasslData->dasslJacobian = setJacobianMethod(threadData, jacobian->availability, flagValue);
+
   /* default use a user sub-routine for JAC */
   dasslData->info[4] = 1;
 

OMCompiler/SimulationRuntime/c/simulation/solver/ida_solver.c

@@ -57,6 +57,7 @@
 #include "epsilon.h"
 #include "external_input.h"
 #include "jacobianSymbolical.h"
+#include "util/jacobian_util.h"
 #include "model_help.h"
 #include "omc_math.h"
 #include "simulation/options.h"
@@ -334,38 +335,6 @@ int ida_solver_initial(DATA* data, threadData_t *threadData,
      infoStreamPrint(LOG_SOLVER, 0, "as the output frequency time step control is used: %f", idaData->stepsTime);
   }
 
-
-  /* if FLAG_JACOBIAN is set, choose jacobian calculation method */
-  if (omc_flag[FLAG_JACOBIAN]) {
-    for (i=1; i< JAC_MAX;i++) {
-      if (!strcmp((const char*)omc_flagValue[FLAG_JACOBIAN], JACOBIAN_METHOD[i])) {
-        idaData->jacobianMethod = (enum JACOBIAN_METHOD)i;
-        break;
-      }
-    }
-    if(idaData->jacobianMethod == SYMJAC) {
-      warningStreamPrint(LOG_STDOUT, 1, "Symbolic Jacobians without coloring are currently not supported by IDA. Colored symbolical Jacobian will be used.");
-      idaData->jacobianMethod = COLOREDSYMJAC;
-    }
-    if(idaData->jacobianMethod == NUMJAC) {
-      warningStreamPrint(LOG_STDOUT, 1, "Numerical Jacobians without coloring are currently not supported by IDA. Colored numerical Jacobian will be used.");
-      idaData->jacobianMethod = COLOREDNUMJAC;
-    }
-    if(idaData->jacobianMethod == JAC_UNKNOWN){
-      if (ACTIVE_WARNING_STREAM(LOG_SOLVER)) {
-        warningStreamPrint(LOG_SOLVER, 1, "unrecognized jacobian calculation method %s, current options are:", (const char*)omc_flagValue[FLAG_JACOBIAN]);
-        for(i=1; i < JAC_MAX; ++i) {
-          warningStreamPrint(LOG_SOLVER, 0, "%-15s [%s]", JACOBIAN_METHOD[i], JACOBIAN_METHOD_DESC[i]);
-        }
-        messageClose(LOG_SOLVER);
-      }
-      throwStreamPrint(threadData,"unrecognized jacobian calculation method %s", (const char*)omc_flagValue[FLAG_JACOBIAN]);
-    }
-    /* default case colored numerical jacobian */
-  } else {
-    idaData->jacobianMethod = COLOREDNUMJAC;
-  }
-
   /* if FLAG_IDA_LS is set, choose ida linear solver method */
   if (omc_flag[FLAG_IDA_LS]) {
     for (i=1; i< IDA_LS_MAX; i++) {
@@ -388,29 +357,37 @@ int ida_solver_initial(DATA* data, threadData_t *threadData,
     idaData->linearSolverMethod = IDA_LS_KLU;
   }
 
-  /* selects the calculation method of the jacobian */
-  /* in daeMode with numerical jacobian sparse pattern is already initialized in DAEMODE_DATA */
-  if((idaData->jacobianMethod == COLOREDNUMJAC && !idaData->daeMode)||
-     idaData->jacobianMethod == COLOREDSYMJAC ||
-     idaData->jacobianMethod == SYMJAC)
-  {
-    ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A]);
-    if (data->callback->initialAnalyticJacobianA(data, threadData, jacobian))
-    {
-      infoStreamPrint(LOG_STDOUT, 0, "Jacobian or SparsePattern is not generated or failed to initialize! Switch back to normal.");
-      idaData->jacobianMethod = INTERNALNUMJAC;
-      if (idaData->linearSolverMethod == IDA_LS_KLU)
-      {
-        idaData->linearSolverMethod = IDA_LS_DENSE;
-        warningStreamPrint(LOG_STDOUT, 0, "IDA linear solver method also switched back to %s", IDA_LS_METHOD_DESC[idaData->linearSolverMethod]);
-      }
-    } else {
-      ANALYTIC_JACOBIAN* jac = &data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A];
-      infoStreamPrint(LOG_SIMULATION, 1, "Initialized colored Jacobian:");
-      infoStreamPrint(LOG_SIMULATION, 0, "columns: %d rows: %d", jac->sizeCols, jac->sizeRows);
-      infoStreamPrint(LOG_SIMULATION, 0, "NNZ:  %d colors: %d", jac->sparsePattern->numberOfNonZeros, jac->sparsePattern->maxColors);
-      messageClose(LOG_SIMULATION);
-    }
+  ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A]);
+  data->callback->initialAnalyticJacobianA(data, threadData, jacobian);
+  if(jacobian->availability == JACOBIAN_AVAILABLE || jacobian->availability == JACOBIAN_ONLY_SPARSITY) {
+    infoStreamPrint(LOG_SIMULATION, 1, "Initialized Jacobian:");
+    infoStreamPrint(LOG_SIMULATION, 0, "columns: %d rows: %d", jacobian->sizeCols, jacobian->sizeRows);
+    infoStreamPrint(LOG_SIMULATION, 0, "NNZ:  %d colors: %d", jacobian->sparsePattern->numberOfNonZeros, jacobian->sparsePattern->maxColors);
+    messageClose(LOG_SIMULATION);
+  }
+
+  // Compare user flag to availabe Jacobian methods
+  const char* flagValue;
+  if(omc_flag[FLAG_JACOBIAN]){
+    flagValue = omc_flagValue[FLAG_JACOBIAN];
+  } else {
+    flagValue = NULL;
+  }
+  idaData->jacobianMethod = setJacobianMethod(threadData, jacobian->availability, flagValue);
+
+  // change IDA specific jacobian method
+  if(idaData->jacobianMethod == SYMJAC) {
+    warningStreamPrint(LOG_STDOUT, 0, "Symbolic Jacobians without coloring are currently not supported by IDA."
+                                      " Colored symbolical Jacobian will be used.");
+    idaData->jacobianMethod = COLOREDSYMJAC;
+  }else if(idaData->jacobianMethod == NUMJAC) {
+    warningStreamPrint(LOG_STDOUT, 0, "Numerical Jacobians without coloring are currently not supported by IDA."
+                                      " Colored numerical Jacobian will be used.");
+    idaData->jacobianMethod = COLOREDNUMJAC;
+  }else if(idaData->jacobianMethod == INTERNALNUMJAC && idaData->linearSolverMethod == IDA_LS_KLU) {
+    warningStreamPrint(LOG_STDOUT, 0, "Internal Numerical Jacobians without coloring are currently not supported by IDA with KLU."
+                                      " Colored numerical Jacobian will be used.");
+    idaData->jacobianMethod = COLOREDNUMJAC;
   }
 
   /* Set NNZ */

OMCompiler/SimulationRuntime/c/simulation_data.h

@@ -119,7 +119,21 @@ typedef struct CALL_STATISTICS
   long functionAlgebraics;
 } CALL_STATISTICS;
 
-typedef enum {ERROR_AT_TIME,NO_PROGRESS_START_POINT,NO_PROGRESS_FACTOR,IMPROPER_INPUT} EQUATION_SYSTEM_ERROR;
+typedef enum
+{
+  ERROR_AT_TIME,
+  NO_PROGRESS_START_POINT,
+  NO_PROGRESS_FACTOR,
+  IMPROPER_INPUT
+} EQUATION_SYSTEM_ERROR;
+
+typedef enum
+{
+  JACOBIAN_UNKNOWN = 0,       /* availability of jacobian unknown (not initialized) */
+  JACOBIAN_NOT_AVAILABLE,     /* no symbolic jacobian and no sparsity pattern available */
+  JACOBIAN_ONLY_SPARSITY,     /* only sparsity pattern available */
+  JACOBIAN_AVAILABLE          /* symbolic jacobian and sparsity pattern available */
+} JACOBIAN_AVAILABILITY;
 
 /**
  * @brief Sparse pattern for Jacobian matrix.
@@ -143,14 +157,15 @@ typedef struct SPARSE_PATTERN
  */
 typedef struct ANALYTIC_JACOBIAN
 {
-  unsigned int sizeCols;          /* Number of columns of Jacobian */
-  unsigned int sizeRows;          /* Number of rows of Jacobian */
-  unsigned int sizeTmpVars;       /* Length of vector tmpVars */
-  SPARSE_PATTERN* sparsePattern;  /* Contain sparse pattern including coloring */
-  modelica_real* seedVars;        /* Seed vector for specifying which columns to evaluate */
+  JACOBIAN_AVAILABILITY availability;  /* Availability status */
+  unsigned int sizeCols;              /* Number of columns of Jacobian */
+  unsigned int sizeRows;              /* Number of rows of Jacobian */
+  unsigned int sizeTmpVars;           /* Length of vector tmpVars */
+  SPARSE_PATTERN* sparsePattern;      /* Contain sparse pattern including coloring */
+  modelica_real* seedVars;            /* Seed vector for specifying which columns to evaluate */
   modelica_real* tmpVars;
-  modelica_real* resultVars;      /* Result column for given seed vector */
-  modelica_real dae_cj;           /* Is the scalar in the system Jacobian, proportional to the inverse of the step size. From User Documentation for ida v5.4.0 equation (2.5). */
+  modelica_real* resultVars;          /* Result column for given seed vector */
+  modelica_real dae_cj;               /* Is the scalar in the system Jacobian, proportional to the inverse of the step size. From User Documentation for ida v5.4.0 equation (2.5). */
   int (*constantEqns)(void* data, threadData_t *threadData, void* thisJacobian, void* parentJacobian);  /* Constant equations independed of seed vector */
 } ANALYTIC_JACOBIAN;
 
@@ -253,7 +268,6 @@ typedef int (*analyticalJacobianColumn_func_ptr)(DATA* data, threadData_t* threa
 /**
  * @brief User data provided to residual functions.
  *
- * Created by (non-)linear solver before evaluating a residual.
  */
 typedef struct RESIDUAL_USERDATA {
   DATA* data;

OMCompiler/SimulationRuntime/c/util/jacobian_util.c

@@ -32,6 +32,7 @@
  */
 
 #include "jacobian_util.h"
+#include "../simulation/options.h"
 
 /**
  * @brief Initialize analytic jacobian.
@@ -55,6 +56,7 @@ void initAnalyticJacobian(ANALYTIC_JACOBIAN* jacobian, unsigned int sizeCols, un
   jacobian->tmpVars = (modelica_real*) calloc(sizeTmpVars, sizeof(modelica_real));
   jacobian->constantEqns = constantEqns;
   jacobian->sparsePattern = sparsePattern;
+  jacobian->availability = JACOBIAN_UNKNOWN;
   jacobian->dae_cj = 0;
 }
 
@@ -129,3 +131,90 @@ void freeSparsePattern(SPARSE_PATTERN *spp) {
     free(spp->leadindex); spp->leadindex = NULL;
   }
 }
+
+/**
+ * @brief Set Jacobian method from user flag and available Jacobian.
+ *
+ * @param threadData              Used for error handling.
+ * @param availability            Is the Jacobian available, only the sparsity pattern available or nothing available.
+ * @param flagValue               Flag value of FLAG_JACOBIAN. Can be NULL.
+ * @return enum JACOBIAN_METHOD   Returns jacobian method that is availble.
+ */
+enum JACOBIAN_METHOD setJacobianMethod(threadData_t* threadData, JACOBIAN_AVAILABILITY availability, const char* flagValue){
+  enum JACOBIAN_METHOD jacobianMethod = JAC_UNKNOWN;
+  assertStreamPrint(threadData, availability != JACOBIAN_UNKNOWN, "Jacobian availablity status is unknown.");
+
+  /* if FLAG_JACOBIAN is set, choose jacobian calculation method */
+  if (flagValue) {
+    for (int method=1; method < JAC_MAX; method++) {
+      if (!strcmp(flagValue, JACOBIAN_METHOD[method])) {
+        jacobianMethod = (enum JACOBIAN_METHOD) method;
+        break;
+      }
+    }
+    // Error case
+    if(jacobianMethod == JAC_UNKNOWN){
+      errorStreamPrint(LOG_STDOUT, 0, "Unknown value `%s` for flag `-jacobian`", flagValue);
+      infoStreamPrint(LOG_STDOUT, 1, "Available options are");
+      for (int method=1; method < JAC_MAX; method++) {
+        infoStreamPrint(LOG_STDOUT, 0, "%s", JACOBIAN_METHOD[method]);
+      }
+      messageClose(LOG_STDOUT);
+      omc_throw(threadData);
+    }
+  }
+
+  /* Check if method is available */
+  switch (availability)
+  {
+  case JACOBIAN_NOT_AVAILABLE:
+    if (jacobianMethod != INTERNALNUMJAC && jacobianMethod != JAC_UNKNOWN) {
+      warningStreamPrint(LOG_STDOUT, 0, "Jacobian not available, switching to internal numerical Jacobian.");
+    }
+    jacobianMethod = INTERNALNUMJAC;
+    break;
+  case JACOBIAN_ONLY_SPARSITY:
+    if (jacobianMethod == COLOREDSYMJAC) {
+      warningStreamPrint(LOG_STDOUT, 0, "Symbolic Jacobian not available, only sparsity pattern. Switching to colored numerical Jacobian.");
+      jacobianMethod = COLOREDNUMJAC;
+    } else if(jacobianMethod == SYMJAC) {
+      warningStreamPrint(LOG_STDOUT, 0, "Symbolic Jacobian not available, only sparsity pattern. Switching to numerical Jacobian.");
+      jacobianMethod = NUMJAC;
+    } else if(jacobianMethod == JAC_UNKNOWN) {
+      jacobianMethod = COLOREDNUMJAC;
+    }
+    break;
+  case JACOBIAN_AVAILABLE:
+    if (jacobianMethod == JAC_UNKNOWN) {
+      jacobianMethod = COLOREDSYMJAC;
+    }
+    break;
+  default:
+    throwStreamPrint(threadData, "Unhandled case in setJacobianMethod");
+    break;
+  }
+
+  /* Log Jacobian method */
+  switch (jacobianMethod)
+  {
+  case INTERNALNUMJAC:
+    infoStreamPrint(LOG_JAC, 0, "Using Jacobian method: Internal numerical Jacobian.");
+    break;
+  case NUMJAC:
+    infoStreamPrint(LOG_JAC, 0, "Using Jacobian method: Numerical Jacobian.");
+    break;
+  case COLOREDNUMJAC:
+    infoStreamPrint(LOG_JAC, 0, "Using Jacobian method: Colored numerical Jacobian.");
+    break;
+  case SYMJAC:
+    infoStreamPrint(LOG_JAC, 0, "Using Jacobian method: Symbolical Jacobian.");
+    break;
+  case COLOREDSYMJAC:
+    infoStreamPrint(LOG_JAC, 0, "Using Jacobian method: Colored symbolical Jacobian.");
+    break;
+  default:
+    throwStreamPrint(threadData, "Unhandled case in setJacobianMethod");
+    break;
+  }
+  return jacobianMethod;
+}

OMCompiler/SimulationRuntime/c/util/jacobian_util.h

@@ -46,6 +46,8 @@ void freeAnalyticJacobian(ANALYTIC_JACOBIAN* jac);
 
 SPARSE_PATTERN* allocSparsePattern(unsigned int n_leadIndex, unsigned int numberOfNonZeros, unsigned int maxColors);
 void freeSparsePattern(SPARSE_PATTERN *spp);
+enum JACOBIAN_METHOD setJacobianMethod(threadData_t* threadData, JACOBIAN_AVAILABILITY availability, const char* flagValue);
+
 
 #ifdef __cplusplus
 }

testsuite/openmodelica/cruntime/optimization/basic/ocpWithInputs.mos

@@ -42,7 +42,6 @@ val(cost,1.0);
 // --------------------------------------------------------
 // number of nonlinear constraints: 0
 // ========================================================
-// stdout            | info    | Jacobian or SparsePattern is not generated or failed to initialize! Switch back to normal.
 //
 // ******************************************************************************
 // This program contains Ipopt, a library for large-scale nonlinear optimization.

testsuite/openmodelica/cruntime/optimization/basic/testAlgLoop1.mos

@@ -48,7 +48,6 @@ getErrorString();
 // --------------------------------------------------------
 // number of nonlinear constraints: 1
 // ========================================================
-// stdout            | info    | Jacobian or SparsePattern is not generated or failed to initialize! Switch back to normal.
 //
 // ******************************************************************************
 // This program contains Ipopt, a library for large-scale nonlinear optimization.