Let Newton solver exploit symbolic Jacobians

OpenModelica · Apr 23, 2016 · 15ff300 · 15ff300
1 parent 6221cae
commit 15ff300
Show file tree

Hide file tree

Showing 2 changed files with 28 additions and 2 deletions.
diff --git a/Compiler/Template/CodegenCpp.tpl b/Compiler/Template/CodegenCpp.tpl
@@ -6247,17 +6247,34 @@ case SIMCODE(modelInfo = MODELINFO(__)) then
    let &preExp= buffer ""
 
   match eq
+  case SES_NONLINEAR(nlSystem = nls as NONLINEARSYSTEM(jacobianMatrix = SOME(jac))) then
+    let jacIndex = match jac
+      case (_, _, _, _, _, _, index) then index else "getSystemMatrix:ERROR"
+  <<
+
+  const matrix_t& <%modelname%>Algloop<%nls.index%>::getSystemMatrix()
+  {
+    return static_cast<<%modelname%>Mixed*>(_system)->getJacobian(<%jacIndex%>);
+  }
+
+  const sparsematrix_t& <%modelname%>Algloop<%nls.index%>::getSystemSparseMatrix()
+  {
+    throw ModelicaSimulationError(MATH_FUNCTION, "Sparse symbolic Jacobians not suported yet");
+  }
+  >>
   case SES_NONLINEAR(nlSystem = nls as NONLINEARSYSTEM(__)) then
   <<
 
   const matrix_t& <%modelname%>Algloop<%nls.index%>::getSystemMatrix()
   {
-    throw ModelicaSimulationError(MATH_FUNCTION,"Symbolic jacobians is not suported yet");
+    // return empty matrix to indicate that no symbolic Jacobian is available
+    static matrix_t empty(0, 0);
+    return empty;
   }
 
   const sparsematrix_t& <%modelname%>Algloop<%nls.index%>::getSystemSparseMatrix()
   {
-    throw ModelicaSimulationError(MATH_FUNCTION,"Symbolic jacobians is not suported yet");
+    throw ModelicaSimulationError(MATH_FUNCTION, "Sparse symbolic Jacobians not suported yet");
   }
   >>
  case SES_LINEAR(lSystem = ls as LINEARSYSTEM(__)) then

diff --git a/SimulationRuntime/cpp/Solver/Newton/Newton.cpp b/SimulationRuntime/cpp/Solver/Newton/Newton.cpp
@@ -347,6 +347,15 @@ void Newton::stepCompleted(double time)
 
 void Newton::calcJacobian()
 {
+  // Use analytic Jacobian if available
+  matrix_t A = _algLoop->getSystemMatrix();
+  if (A.size1() == _dimSys && A.size2() == _dimSys) {
+    const double* jac = A.data().begin();
+    std::copy(jac, jac + _dimSys*_dimSys, _jac);
+    return;
+  }
+
+  // Alternatively apply finite differences
   for (int j = 0; j < _dimSys; ++j) {
     // Reset variables for every column
     std::copy(_y, _y + _dimSys, _yHelp);