Eval jacobian (#113)

* Evaluate Jacobian of non-linear system.

src/NonLinearSystemNeuralNetworkFMU.jl

@@ -41,6 +41,8 @@ if !isdefined(Base, :get_extension)
 end
 
 include("types.jl")
+export fmiEvaluateRes
+export fmiEvaluateEq
 export EqInfo
 export MinMaxBoundaryValues
 export ProfilingInfo

src/fmiExtensions.jl

@@ -72,7 +72,7 @@ for given equation number.
   - `x::Array{Float64}` Values of iteration variables at which to evaluate the residuals.
 
 # Returns
-  - Status of Libdl.ccall for `:myfmi2EvaluateEq`.
+  - Status of Libdl.ccall for `:myfmi2EvaluateRes`.
   - Array of residual values.
 """
 function fmiEvaluateRes(fmu::FMIImport.FMU2, eqNumber::Integer, x::Array{Float64})::Tuple{fmi2Status, Array{Float64}}
@@ -96,3 +96,41 @@ function fmiEvaluateRes(comp::FMICore.FMU2Component, eq::Integer, x::Array{Float
 
   return status, res
 end
+
+
+"""
+    fmiEvaluateJacobian(fmu, eqNumber, x)
+
+Evaluate Jacobian Matrix of given equation system 'eq' at a vector of iteration
+variables 'x'.
+
+# Arguments
+  - `fmu::FMIImport.FMU2`: FMU object containing C void pointer to FMU component.
+  - `eqNumber::Int`: Equation index specifying equation to evaluate jacobian of.
+  - `x::Array{Float64}` Values of iteration variables at which to evaluate the
+    jacobian.
+
+# Returns
+  - Status of Libdl.ccall for `:myfmi2EvaluateJacobian`.
+  - jacobian of eq at x.
+"""
+function fmiEvaluateJacobian(fmu::FMIImport.FMU2, eqNumber::Integer, x::Array{Float64})::Tuple{fmi2Status, Array{Float64}}
+  return fmiEvaluateJacobian(fmu.components[1], eqNumber, x)
+end
+
+function fmiEvaluateJacobian(comp::FMICore.FMU2Component, eqNumber::Integer, x::Array{Float64})::Tuple{fmi2Status, Array{Float64}}
+  @assert eqNumber>=0 "Residual index has to be non-negative!"
+
+  # this is a pointer to Jacobian matrix in row-major-format or NULL in error case.
+  fmiEvaluateJacobian = Libdl.dlsym(comp.fmu.libHandle, :myfmi2EvaluateJacobian)
+
+  jac = Array{Float64}(undef, length(x)*length(x))
+
+  eqCtype = Csize_t(eqNumber)
+
+  status = ccall(fmiEvaluateJacobian,
+                 Cuint,
+                 (Ptr{Nothing}, Csize_t, Ptr{Cdouble}, Ptr{Cdouble}),
+                 comp.compAddr, eqCtype, x, jac)
+  return status, jac
+end

src/templates/special_interface.tpl.c

@@ -149,27 +149,21 @@ fmi2Status myfmi2EvaluateRes(fmi2Component c, const size_t eqNumber, double* x,
 }
 
 /**
- * @brief Get Jacobian of non-linear equation system.
+ * @brief Get pointer to Jacobian of non-linear equation system.
  *
  * @param c             Pointer to FMU component.
  * @param sysNumber     Number of non-linear system.
- * @return jac          Return pointer to Jacobian matrix in row-major-format or NULL in error case.
+ * @return jac          Return pointer to Jacobian matrix in row-major-format or
+ *                      NULL in error case.
  */
 double* getJac(DATA* data, const size_t sysNumber) {
   double* jac;
   NONLINEAR_SYSTEM_DATA* nlsSystem = &(data->simulationInfo->nonlinearSystemData[sysNumber]);
 
   switch (nlsSystem->nlsMethod)
   {
-  //case NLS_HYBRID:
-  //  DATA_HYBRD* solverData = (DATA_HYBRD*) nlsSystem->solverData;
-  //  jac = solverData->fjac;
-  //  break;
-  //case NLS_NEWTON:
-  //  DATA_NEWTON* solverData = (DATA_NEWTON*) nlsSystem->solverData;
-  //  jac = solverData->fjac;
-  //  break;
   case NLS_HOMOTOPY:
+    // Get pointer from DATA_HOMOTOPY struct
     return getHomotopyJacobian(nlsSystem);
   default:
     printf("Unknown NLS method  %d in myfmi2GetJac\n", (int)nlsSystem->nlsMethod);
@@ -204,3 +198,37 @@ int scaleResidual(double* jac, double* res, size_t n) {
 
   return isRegular;
 }
+
+/**
+ * @brief Evaluate Jacobian for a given nonlinear system.
+ *
+ * This function evaluates the Jacobian matrix for a given nonlinear system.
+ * The Jacobian matrix is computed based on the method specified in the system's
+ * configuration. Currently, only the homotopy method is supported.
+ *
+ * @param c             Pointer to FMU component.
+ * @param sysNumber     Number of non-linear system.
+ * @param x             Iteration variables of non-linear system.
+ * @param jac           Pointer to allocated memory of size length(x)^2.
+ *                      On exit values of Jacobian matrix in row-major-format
+ * @return fmi2Status
+ */
+fmi2Status myfmi2EvaluateJacobian(fmi2Component c, const size_t sysNumber, double* x, double* jac)
+{
+  ModelInstance *comp = (ModelInstance *)c;
+  DATA* data = comp->fmuData;
+  threadData_t *threadData = comp->threadData;
+  NONLINEAR_SYSTEM_DATA* nlsSystem = &(data->simulationInfo->nonlinearSystemData[sysNumber]);
+
+  switch(nlsSystem->nlsMethod)
+  {
+    case NLS_HOMOTOPY:
+      DATA_HOMOTOPY* solverData = (DATA_HOMOTOPY*) nlsSystem->solverData;
+      int status = getAnalyticalJacobianHomotopy(solverData, jac);
+      break;
+    default:
+      printf("Unknown NLS method  %d in myfmi2GetJac\n", (int)nlsSystem->nlsMethod);
+      return fmi2Fatal;
+  }
+  return fmi2OK;
+}

src/templates/special_interface.tpl.h

@@ -28,6 +28,7 @@ extern "C" {
 
 FMI2_Export fmi2Status myfmi2EvaluateEq(fmi2Component c, const size_t eqNumber);
 fmi2Status myfmi2EvaluateRes(fmi2Component c, const size_t eqNumber, double* x, double* res);
+fmi2Status myfmi2EvaluateJacobian(fmi2Component c, const size_t eqNumber, double* x, double* res);
 double* getJac(DATA* data, const size_t sysNumber);
 int scaleResidual(double* jac, double* res, size_t n);