[C-Runtime] Parallelise DASSL and IDA Jacobian evaluation

Add support for parallel Jacobian evaluation in DASSL and IDA for symbolical jacobian of right hand side of DAE. For this linear systems can now be solved in parallel with all linear solvers, since a Jacobian column can contain linear loops. - Parallel Jacobian evaluation will be encapsulated by USE_PARJAC defines and compiler directives. - To compile omc with parallel jacobians pass `--enable-parjac` to `OMCompiler/configure`. This will check for OpenMP support and add corresponding Flags to simulation makefile and c runtime makefile. - To simulate a model with parallel jacobian evaluations use a) Pass `-jacobianThreads=<numberOfThreads>` as simulationflag to use desired ammount of threads for parallel jacobian evaluation. b) environment variable `OMP_NUM_THREADS=N` This order also gives the order of the precedence, i.e. if both are specified value of -jacobianThreads is taken over `OMP_NUM_THREADS`. If nothing is specified by the user `omp_get_max_thread()` is used to set the number of threads to use for parallel jacobian evaluation. Co-authored-by: mflehmig <martin.schroschk@tu-dresden.de>
OpenModelica · Nov 5, 2019 · aeb08ca · aeb08ca
1 parent 549fe86
commit aeb08ca
Show file tree

Hide file tree

Showing 34 changed files with 976 additions and 276 deletions.
diff --git a/OMCompiler/Compiler/Template/CodegenC.tpl b/OMCompiler/Compiler/Template/CodegenC.tpl
@@ -117,6 +117,7 @@ end translateModel;
     #include "simulation/simulation_info_json.h"
     #include "simulation/simulation_runtime.h"
     #include "util/omc_error.h"
+    #include "util/parallel_helper.h"
     #include "simulation/solver/model_help.h"
     #include "simulation/solver/delay.h"
     #include "simulation/solver/linearSystem.h"
@@ -1957,7 +1958,6 @@ template functionInitialLinearSystemsTemp(list<SimEqSystem> linearSystems, Strin
            linearSystemData[<%ls.indexLinearSystem%>].setA = setLinearMatrixA<%ls.index%>;
            linearSystemData[<%ls.indexLinearSystem%>].setb = setLinearVectorb<%ls.index%>;
            linearSystemData[<%ls.indexLinearSystem%>].initializeStaticLSData = initializeStaticLSData<%ls.index%>;
-           linearSystemData[<%ls.indexLinearSystem%>].parentJacobian = NULL;
            >>
          case SOME(__) then
            let size = listLength(ls.vars)
@@ -1976,10 +1976,9 @@ template functionInitialLinearSystemsTemp(list<SimEqSystem> linearSystems, Strin
            linearSystemData[<%ls.indexLinearSystem%>].analyticalJacobianColumn = <%generatedJac%>;
            linearSystemData[<%ls.indexLinearSystem%>].initialAnalyticalJacobian = <%initialJac%>;
            linearSystemData[<%ls.indexLinearSystem%>].jacobianIndex = <%jacIndex%>;
-           linearSystemData[<%ls.indexLinearSystem%>].setA = NULL;//setLinearMatrixA<%ls.index%>;
-           linearSystemData[<%ls.indexLinearSystem%>].setb = NULL; //setLinearVectorb<%ls.index%>;
+           linearSystemData[<%ls.indexLinearSystem%>].setA = NULL;  //setLinearMatrixA<%ls.index%>;
+           linearSystemData[<%ls.indexLinearSystem%>].setb = NULL;  //setLinearVectorb<%ls.index%>;
            linearSystemData[<%ls.indexLinearSystem%>].initializeStaticLSData = initializeStaticLSData<%ls.index%>;
-           linearSystemData[<%ls.indexLinearSystem%>].parentJacobian = NULL;
            >>
          else
          error(sourceInfo(), ' No jacobian create for linear system <%ls.index%>.')
@@ -2006,7 +2005,6 @@ template functionInitialLinearSystemsTemp(list<SimEqSystem> linearSystems, Strin
            linearSystemData[<%ls.indexLinearSystem%>].setA = setLinearMatrixA<%ls.index%>;
            linearSystemData[<%ls.indexLinearSystem%>].setb = setLinearVectorb<%ls.index%>;
            linearSystemData[<%ls.indexLinearSystem%>].initializeStaticLSData = initializeStaticLSData<%ls.index%>;
-           linearSystemData[<%ls.indexLinearSystem%>].parentJacobian = NULL;
 
            assertStreamPrint(NULL, nLinearSystems > <%at.indexLinearSystem%>, "Internal Error: nLinearSystems mismatch!");
            linearSystemData[<%at.indexLinearSystem%>].equationIndex = <%at.index%>;
@@ -2017,7 +2015,6 @@ template functionInitialLinearSystemsTemp(list<SimEqSystem> linearSystems, Strin
            linearSystemData[<%at.indexLinearSystem%>].setA = setLinearMatrixA<%at.index%>;
            linearSystemData[<%at.indexLinearSystem%>].setb = setLinearVectorb<%at.index%>;
            linearSystemData[<%at.indexLinearSystem%>].initializeStaticLSData = initializeStaticLSData<%at.index%>;
-           linearSystemData[<%ls.indexLinearSystem%>].parentJacobian = NULL;
            linearSystemData[<%at.indexLinearSystem%>].checkConstraints = checkConstraints<%at.index%>;
            >>
          case SOME(__) then
@@ -2043,9 +2040,8 @@ template functionInitialLinearSystemsTemp(list<SimEqSystem> linearSystems, Strin
            linearSystemData[<%ls.indexLinearSystem%>].analyticalJacobianColumn = <%generatedJac%>;
            linearSystemData[<%ls.indexLinearSystem%>].initialAnalyticalJacobian = <%initialJac%>;
            linearSystemData[<%ls.indexLinearSystem%>].jacobianIndex = <%jacIndex%>;
-           linearSystemData[<%ls.indexLinearSystem%>].parentJacobian = NULL;
-           linearSystemData[<%ls.indexLinearSystem%>].setA = NULL;//setLinearMatrixA<%ls.index%>;
-           linearSystemData[<%ls.indexLinearSystem%>].setb = NULL; //setLinearVectorb<%ls.index%>;
+           linearSystemData[<%ls.indexLinearSystem%>].setA = NULL;  //setLinearMatrixA<%ls.index%>;
+           linearSystemData[<%ls.indexLinearSystem%>].setb = NULL;  //setLinearVectorb<%ls.index%>;
            linearSystemData[<%ls.indexLinearSystem%>].initializeStaticLSData = initializeStaticLSData<%ls.index%>;
 
            assertStreamPrint(NULL, nLinearSystems > <%at.indexLinearSystem%>, "Internal Error: indexlinearSystem mismatch!");
@@ -2058,9 +2054,8 @@ template functionInitialLinearSystemsTemp(list<SimEqSystem> linearSystems, Strin
            linearSystemData[<%at.indexLinearSystem%>].analyticalJacobianColumn = <%generatedJac2%>;
            linearSystemData[<%at.indexLinearSystem%>].initialAnalyticalJacobian = <%initialJac2%>;
            linearSystemData[<%at.indexLinearSystem%>].jacobianIndex = <%jacIndex2%>;
-           linearSystemData[<%ls.indexLinearSystem%>].parentJacobian = NULL;
-           linearSystemData[<%at.indexLinearSystem%>].setA = NULL;//setLinearMatrixA<%at.index%>;
-           linearSystemData[<%at.indexLinearSystem%>].setb = NULL; //setLinearVectorb<%at.index%>;
+           linearSystemData[<%at.indexLinearSystem%>].setA = NULL;  //setLinearMatrixA<%at.index%>;
+           linearSystemData[<%at.indexLinearSystem%>].setb = NULL;  //setLinearVectorb<%at.index%>;
            linearSystemData[<%at.indexLinearSystem%>].initializeStaticLSData = initializeStaticLSData<%at.index%>;
            linearSystemData[<%at.indexLinearSystem%>].checkConstraints = checkConstraints<%at.index%>;
            >>
@@ -2131,7 +2126,7 @@ template functionSetupLinearSystemsTemp(list<SimEqSystem> linearSystems, String
          threadData_t *threadData = (threadData_t*) ((void**)dataIn[1]);
          const int equationIndexes[2] = {1,<%ls.index%>};
          <% if ls.partOfJac then
-           'ANALYTIC_JACOBIAN* parentJacobian = data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].parentJacobian;'
+           'ANALYTIC_JACOBIAN* parentJacobian = data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].parDynamicData[omc_get_thread_num()].parentJacobian;'
          %>
          ANALYTIC_JACOBIAN* jacobian = NULL;
          <%varDeclsRes%>
@@ -2182,7 +2177,7 @@ template functionSetupLinearSystemsTemp(list<SimEqSystem> linearSystems, String
          const int equationIndexes[2] = {1,<%ls.index%>};
          DATA* data = (DATA*) inData;
          LINEAR_SYSTEM_DATA* linearSystemData = (LINEAR_SYSTEM_DATA*) systemData;
-         <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parentJacobian;'%>
+         <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parDynamicData[omc_get_thread_num()].parentJacobian;'%>
          <%varDecls%>
          <%MatrixA%>
        }
@@ -2192,7 +2187,7 @@ template functionSetupLinearSystemsTemp(list<SimEqSystem> linearSystems, String
          const int equationIndexes[2] = {1,<%ls.index%>};
          DATA* data = (DATA*) inData;
          LINEAR_SYSTEM_DATA* linearSystemData = (LINEAR_SYSTEM_DATA*) systemData;
-         <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parentJacobian;'%>
+         <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parDynamicData[omc_get_thread_num()].parentJacobian;'%>
          <%varDecls2%>
          <%vectorb%>
        }
@@ -2270,7 +2265,9 @@ template functionSetupLinearSystemsTemp(list<SimEqSystem> linearSystems, String
          DATA *data = (DATA*) ((void**)dataIn[0]);
          threadData_t *threadData = (threadData_t*) ((void**)dataIn[1]);
          const int equationIndexes[2] = {1,<%ls.index%>};
-         <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].parentJacobian;'%>
+         <% if ls.partOfJac then
+           'ANALYTIC_JACOBIAN* parentJacobian = data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].parDynamicData[omc_get_thread_num()].parentJacobian;'
+         %>
          ANALYTIC_JACOBIAN* jacobian = NULL;
          <%varDeclsRes%>
          <% if profileAll() then 'SIM_PROF_TICK_EQ(<%ls.index%>);' %>
@@ -2299,7 +2296,9 @@ template functionSetupLinearSystemsTemp(list<SimEqSystem> linearSystems, String
          DATA *data = (DATA*) ((void**)dataIn[0]);
          threadData_t *threadData = (threadData_t*) ((void**)dataIn[1]);
          const int equationIndexes[2] = {1,<%at.index%>};
-         <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].parentJacobian;'%>
+         <% if ls.partOfJac then
+           'ANALYTIC_JACOBIAN* parentJacobian = data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].parDynamicData[omc_get_thread_num()].parentJacobian;'
+         %>
          ANALYTIC_JACOBIAN* jacobian = NULL;
          <%varDeclsRes2%>
          <% if profileAll() then 'SIM_PROF_TICK_EQ(<%at.index%>);' %>
@@ -2373,7 +2372,7 @@ template functionSetupLinearSystemsTemp(list<SimEqSystem> linearSystems, String
          const int equationIndexes[2] = {1,<%ls.index%>};
          DATA* data = (DATA*) inData;
          LINEAR_SYSTEM_DATA* linearSystemData = (LINEAR_SYSTEM_DATA*) systemData;
-         <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parentJacobian;'%>
+         <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parDynamicData[omc_get_thread_num()].parentJacobian;'%>
          <%varDecls%>
          <%MatrixA%>
        }
@@ -2383,7 +2382,7 @@ template functionSetupLinearSystemsTemp(list<SimEqSystem> linearSystems, String
          const int equationIndexes[2] = {1,<%ls.index%>};
          DATA* data = (DATA*) inData;
          LINEAR_SYSTEM_DATA* linearSystemData = (LINEAR_SYSTEM_DATA*) systemData;
-         <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parentJacobian;'%>
+         <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parDynamicData[omc_get_thread_num()].parentJacobian;'%>
          <%varDecls2%>
          <%vectorb%>
        }
@@ -2403,7 +2402,7 @@ template functionSetupLinearSystemsTemp(list<SimEqSystem> linearSystems, String
          const int equationIndexes[2] = {1,<%at.index%>};
          DATA* data = (DATA*) inData;
          LINEAR_SYSTEM_DATA* linearSystemData = (LINEAR_SYSTEM_DATA*) systemData;
-         <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parentJacobian;'%>
+         <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parDynamicData[omc_get_thread_num()].parentJacobian;'%>
          <%varDecls3%>
          <%MatrixA2%>
        }
@@ -2413,7 +2412,7 @@ template functionSetupLinearSystemsTemp(list<SimEqSystem> linearSystems, String
          const int equationIndexes[2] = {1,<%at.index%>};
          DATA* data = (DATA*) inData;
          LINEAR_SYSTEM_DATA* linearSystemData = (LINEAR_SYSTEM_DATA*) systemData;
-         <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parentJacobian;'%>
+         <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parDynamicData[omc_get_thread_num()].parentJacobian;'%>
          <%varDecls4%>
          <%vectorb2%>
        }
@@ -5585,7 +5584,7 @@ case e as SES_LINEAR(lSystem=ls as LINEARSYSTEM(__), alternativeTearing = at) th
   }
   <% if profileSome() then 'SIM_PROF_TICK_EQ(modelInfoGetEquation(&data->modelData->modelDataXml,<%ls.index%>).profileBlockIndex);' %>
   <% if ls.partOfJac then
-     'data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].parentJacobian = jacobian;'
+     'data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].parDynamicData[omc_get_thread_num()].parentJacobian = jacobian;'
   %>
 
   retValue = solve_linear_system(data, threadData, <%ls.indexLinearSystem%>, &aux_x[0]);
@@ -5599,6 +5598,7 @@ case e as SES_LINEAR(lSystem=ls as LINEARSYSTEM(__), alternativeTearing = at) th
   /* write solution */
   <%ls.vars |> SIMVAR(__) hasindex i0 => '<%contextCref(name, context, auxFunctions)%> = aux_x[<%i0%>];' ;separator="\n"%>
   <% if profileSome() then 'SIM_PROF_ACC_EQ(modelInfoGetEquation(&data->modelData->modelDataXml,<%ls.index%>).profileBlockIndex);' %>
+
   <%returnval%>
   >>
 end equationLinear;

diff --git a/OMCompiler/SimulationRuntime/c/Makefile.common b/OMCompiler/SimulationRuntime/c/Makefile.common
@@ -103,6 +103,7 @@ RUNTIMEUTIL_HEADERS = \
 ./util/generic_array.h \
 ./util/index_spec.h \
 ./util/integer_array.h \
+./util/jacobian_util.h \
 ./util/java_interface.h \
 ./util/modelica.h \
 ./util/modelica_string.h \
@@ -111,6 +112,7 @@ RUNTIMEUTIL_HEADERS = \
 ./util/omc_mmap.h \
 ./util/omc_msvc.h \
 ./util/omc_spinlock.h \
+./util/parallel_helper.h \
 ./util/read_matlab4.c \
 ./util/read_matlab4.h \
 ./util/read_csv.c \

diff --git a/OMCompiler/SimulationRuntime/c/Makefile.objs b/OMCompiler/SimulationRuntime/c/Makefile.objs
@@ -28,7 +28,7 @@ else
 UTIL_OBJS_NO_FMI=
 endif
 
-UTIL_OBJS_MINIMAL=base_array$(OBJ_EXT) boolean_array$(OBJ_EXT) omc_error$(OBJ_EXT) omc_file$(OBJ_EXT) division$(OBJ_EXT) generic_array$(OBJ_EXT) index_spec$(OBJ_EXT) integer_array$(OBJ_EXT) list$(OBJ_EXT) modelica_string$(OBJ_EXT) real_array$(OBJ_EXT) ringbuffer$(OBJ_EXT) string_array$(OBJ_EXT) utility$(OBJ_EXT) varinfo$(OBJ_EXT) ModelicaUtilities$(OBJ_EXT) omc_msvc$(OBJ_EXT) simulation_options$(OBJ_EXT) rational$(OBJ_EXT) modelica_string_lit$(OBJ_EXT) omc_init$(OBJ_EXT) omc_mmap$(OBJ_EXT) $(UTIL_OBJS_NO_FMI)
+UTIL_OBJS_MINIMAL=base_array$(OBJ_EXT) boolean_array$(OBJ_EXT) omc_error$(OBJ_EXT) omc_file$(OBJ_EXT) division$(OBJ_EXT) generic_array$(OBJ_EXT) index_spec$(OBJ_EXT) integer_array$(OBJ_EXT) list$(OBJ_EXT) modelica_string$(OBJ_EXT) real_array$(OBJ_EXT) ringbuffer$(OBJ_EXT) string_array$(OBJ_EXT) utility$(OBJ_EXT) varinfo$(OBJ_EXT) ModelicaUtilities$(OBJ_EXT) omc_msvc$(OBJ_EXT) parallel_helper$(OBJ_EXT) simulation_options$(OBJ_EXT) rational$(OBJ_EXT) modelica_string_lit$(OBJ_EXT) omc_init$(OBJ_EXT) omc_mmap$(OBJ_EXT) jacobian_util$(OBJ_EXT) $(UTIL_OBJS_NO_FMI)
 UTIL_HFILES_MINIMAL=base_array.h boolean_array.h division.h generic_array.h omc_error.h omc_file.h index_spec.h integer_array.h list.h modelica.h modelica_string.h read_write.h real_array.h ringbuffer.h rtclock.h string_array.h utility.h varinfo.h simulation_options.h omc_mmap.h modelica_string_lit.h omc_init.h
 
 ifeq ($(OMC_MINIMAL_RUNTIME),)

diff --git a/OMCompiler/SimulationRuntime/c/simulation/simulation_runtime.cpp b/OMCompiler/SimulationRuntime/c/simulation/simulation_runtime.cpp
@@ -86,6 +86,7 @@
 #include "simulation/solver/initialization/initialization.h"
 #include "simulation/solver/dae_mode.h"
 #include "dataReconciliation/dataReconciliation.h"
+#include "util/parallel_helper.h"
 
 #ifdef _OMC_QSS_LIB
   #include "solver_qss/solver_qss.h"
@@ -940,6 +941,36 @@ int initRuntimeAndSimulation(int argc, char**argv, DATA *data, threadData_t *thr
   data->simulationInfo->minStepSize = 4.0 * DBL_EPSILON * fmax(fabs(data->simulationInfo->startTime),fabs(data->simulationInfo->stopTime));
   rt_accumulate(SIM_TIMER_INIT_XML);
 
+  /* Set the maximum number of threads prior to any allocation w.r.t.
+   * linear systems and Jacobians in order to avoid memory leaks.
+   */
+#ifdef USE_PARJAC
+  int num_threads = omc_get_max_threads();
+  if (omc_flag[FLAG_JACOBIAN_THREADS]) {
+    int num_threads_tmp = atoi(omc_flagValue[FLAG_JACOBIAN_THREADS]);
+    infoStreamPrint(LOG_STDOUT, 0,
+         "Number of threads passed via -jacobianThreads: %d",
+         num_threads_tmp);
+    if (0 >= num_threads_tmp) {
+      warningStreamPrint(LOG_STDOUT, 0,
+          "Number of desired OpenMP threads for parallel Jacobian evaluation is <= 0.");
+      warningStreamPrint(LOG_STDOUT, 0, "Use omp_get_max_threads().");
+    } else {
+      num_threads = num_threads_tmp;
+    }
+  }
+  omp_set_num_threads(num_threads);
+
+  infoStreamPrint(LOG_STDOUT, 0,
+      "Number of OpenMP threads for parallel Jacobian evaluation: %d",
+      omc_get_max_threads());
+#else
+  if (omc_flag[FLAG_JACOBIAN_THREADS]) {
+      warningStreamPrint(LOG_STDOUT, 0,
+          "Simulation flag jacobianThreads not available. Make sure you have configured omc with \"--enable-parjac\" and build with a compiler supporting OpenMP.");
+  }
+#endif
+
   /* initialize static data of mixed/linear/non-linear system solvers */
   initializeMixedSystems(data, threadData);
   initializeLinearSystems(data, threadData);