[NB] update DAE mode jacobians (#14137)

* [NB] update DAE mode jacobians

  - update pipeline to only do required DAE mode jacobians and not ODE jacobians
  - update pipeline to do modules in DAE mode and remove them from the main DAE module
  - remove empty partitions from DAE mode (full discrete)
  - when differentiating algebraic loops for jacobians they are known to be linear
  - [SimCode] use correct DAE mode jacobian instead of ODE jacobian
  - [NB] add missing DAEMode residual variable checks
  - works towards #13980

Author: kabdelhak

OMCompiler/Compiler/NBackEnd/Classes/NBackendDAE.mo

@@ -267,6 +267,7 @@ public
     list<String> followEquations = Flags.getConfigStringList(Flags.DEBUG_FOLLOW_EQUATIONS);
     Option<UnorderedSet<String>> eq_filter_opt;
     list<DAE.InlineType> inline_types = {DAE.NORM_INLINE(), DAE.BUILTIN_EARLY_INLINE(), DAE.EARLY_INLINE(), DAE.DEFAULT_INLINE()};
+    NBPartition.Kind kind;
   algorithm
     // if we filter dump for equations
     if listEmpty(followEquations) then
@@ -292,10 +293,13 @@ public
 
     if Flags.getConfigBool(Flags.DAE_MODE) then
       mainModules := {(DAEMode.main, "DAE-Mode")};
+      kind := NBPartition.Kind.DAE;
     else
       mainModules := {};
+      kind := NBPartition.Kind.ODE;
     end if;
 
+    // all main modules are always done in ODE mode
     mainModules := listAppend({
       (function Partitioning.main(kind = NBPartition.Kind.ODE),             "Partitioning"),
       (function Causalize.main(kind = NBPartition.Kind.ODE),                "Causalize"),
@@ -307,10 +311,10 @@ public
     // (do not change order SOLVE -> JACOBIAN)
     postOptModules := {
       (Evaluation.removeDummies,    "Remove Dummies"),
-      (function Tearing.main(kind = NBPartition.Kind.ODE),    "Tearing"),
-      (Partitioning.categorize,                               "Categorize"),
-      (Solve.main,                                            "Solve"),
-      (function Jacobian.main(kind = NBPartition.Kind.ODE),   "Jacobian")
+      (function Tearing.main(kind = kind),    "Tearing"),
+      (Partitioning.categorize,               "Categorize"),
+      (Solve.main,                            "Solve"),
+      (function Jacobian.main(kind = kind),   "Jacobian")
     };
 
     (bdae, preOptClocks)  := applyModules(bdae, preOptModules, eq_filter_opt, ClockIndexes.RT_CLOCK_NEW_BACKEND_MODULE);

OMCompiler/Compiler/NBackEnd/Modules/1_Main/NBDAEMode.mo

@@ -75,10 +75,8 @@ public
         then fail();
       end match;
 
-      // Modules
+      // fake causalize system to fulfill pipeline requirements
       bdae := Causalize.main(bdae, NBPartition.Kind.DAE);
-      bdae := Tearing.main(bdae, NBPartition.Kind.DAE);
-      bdae := Jacobian.main(bdae, NBPartition.Kind.DAE);
     else
       Error.addMessage(Error.INTERNAL_ERROR, {getInstanceName() + " failed."});
     end try;
@@ -130,7 +128,7 @@ protected
           // remove strong components
           part.strongComponents := NONE();
           // accumulate new partitions
-          then part :: new_partitions;
+          then if Partition.Partition.isEmpty(part) then new_partitions else part :: new_partitions;
         else new_partitions;
       end match;
     end for;

OMCompiler/Compiler/NBackEnd/Modules/3_Post/NBJacobian.mo

@@ -137,6 +137,7 @@ public
           () := match kind
             case NBPartition.Kind.ODE algorithm bdae.ode := newPartitions; then ();
             case NBPartition.Kind.DAE algorithm bdae.dae := SOME(newPartitions); then ();
+            else ();
           end match;
           bdae.ode_event := newEvents;
           bdae.funcTree := funcTree;
@@ -421,7 +422,7 @@ public
           // create row-wise sparsity pattern
           for cref in listReverse(partial_vars) loop
             // only create rows for derivatives
-            if jacType == JacobianType.NLS or BVariable.checkCref(cref, BVariable.isStateDerivative) then
+            if jacType == JacobianType.NLS or BVariable.checkCref(cref, BVariable.isStateDerivative) or BVariable.checkCref(cref, BVariable.isResidual) then
               if UnorderedMap.contains(cref, map) then
                 tmp := UnorderedSet.unique_list(UnorderedMap.getOrFail(cref, map), ComponentRef.hash, ComponentRef.isEqual);
                 rows := (cref, tmp) :: rows;
@@ -539,7 +540,8 @@ public
       if jacType == JacobianType.NLS then
         partials := listArray(sparsityPattern.partial_vars);
       else
-        partials := listArray(list(cref for cref guard(BVariable.checkCref(cref, BVariable.isStateDerivative)) in sparsityPattern.partial_vars));
+        partials := listArray(list(cref for cref guard(BVariable.checkCref(cref, BVariable.isStateDerivative) or
+          BVariable.checkCref(cref, BVariable.isResidual)) in sparsityPattern.partial_vars));
       end if;
 
       // create cref -> index maps
@@ -674,17 +676,23 @@ protected
     input String name                                     "Context name for jacobian";
     input Module.jacobianInterface func;
   protected
+    JacobianType jacType;
+    VariablePointers unknowns;
     list<Pointer<Variable>> derivative_vars, state_vars;
     VariablePointers seedCandidates, partialCandidates;
     Option<Jacobian> jacobian                             "Resulting jacobian";
     Partition.Kind kind = Partition.Partition.getKind(part);
   algorithm
     partialCandidates := part.unknowns;
-    derivative_vars := list(var for var guard(BVariable.isStateDerivative(var)) in VariablePointers.toList(part.unknowns));
+    unknowns  := if Partition.Partition.getKind(part) == NBPartition.Kind.DAE then Util.getOption(part.daeUnknowns) else part.unknowns;
+    jacType   := if Partition.Partition.getKind(part) == NBPartition.Kind.DAE then JacobianType.DAE else JacobianType.ODE;
+
+    derivative_vars := list(var for var guard(BVariable.isStateDerivative(var)) in VariablePointers.toList(unknowns));
     state_vars := list(Util.getOption(BVariable.getVarState(var)) for var in derivative_vars);
     seedCandidates := VariablePointers.fromList(state_vars, partialCandidates.scalarized);
 
-    (jacobian, funcTree) := func(name, JacobianType.ODE, seedCandidates, partialCandidates, part.equations, knowns, part.strongComponents, funcTree, kind ==  NBPartition.Kind.INI);
+    (jacobian, funcTree) := func(name, jacType, seedCandidates, partialCandidates, part.equations, knowns, part.strongComponents, funcTree, kind ==  NBPartition.Kind.INI);
+
     part.association := Partition.Association.CONTINUOUS(kind, jacobian);
     if Flags.isSet(Flags.JAC_DUMP) then
       print(Partition.Partition.toString(part, 2));
@@ -821,14 +829,14 @@ protected
   end jacobianNone;
 
   function getTmpFilterFunction
-    " - ODE/DAE filter by state derivative / algebraic
-      - LS/NLS filter by residual / inner"
+    " - ODE filter by state derivative / algebraic
+      - LS/NLS/DAE filter by residual / inner"
     input JacobianType jacType;
     output BVariable.checkVar func;
   algorithm
     func := match jacType
       case JacobianType.ODE then BVariable.isStateDerivative;
-      case JacobianType.DAE then BVariable.isStateDerivative;
+      case JacobianType.DAE then BVariable.isResidual;
       case JacobianType.LS  then BVariable.isResidual;
       case JacobianType.NLS then BVariable.isResidual;
       else algorithm

OMCompiler/Compiler/NBackEnd/Modules/3_Post/NBTearing.mo

@@ -149,7 +149,9 @@ public
           (partitions, funcTree) := tearingTraverser(partitions, funcs, funcTree, eq_index, kind);
           bdae.dae := SOME(partitions);
           bdae.funcTree := funcTree;
-      then bdae;
+          // recursively call this function to also apply to the ODE section (used for events)
+          // ToDo: only create event partitions, disregard rest
+      then main(bdae, NBPartition.Kind.ODE);
 
     // ToDo: all the other cases: e.g. Jacobian, Hessian
     end match;

OMCompiler/Compiler/NBackEnd/Modules/NBModule.mo

@@ -272,7 +272,7 @@ public
     [!] This function can not only be used as an optimization module but also for
     nonlinear partitions, state sets, linearization and dynamic optimization."
     input String name                                     "Name of jacobian";
-    input JacobianType jacType                            "Type of jacobian (sim/nonlin)";
+    input JacobianType jacType                            "Type of jacobian (ode/dae/lin/nonlin)";
     input VariablePointers seedCandidates                 "differentiate by these";
     input VariablePointers partialCandidates              "solve the equations for these";
     input EquationPointers equations                      "Equations array";

OMCompiler/Compiler/NBackEnd/Util/NBDifferentiate.mo

@@ -182,6 +182,7 @@ public
         DifferentiationArguments diffArguments;
         Tearing strict;
         Option<Tearing> casual;
+        Boolean linear;
 
       case StrongComponent.SINGLE_COMPONENT() algorithm
         new_var := differentiateVariablePointer(comp.var, diffArguments_ptr);
@@ -221,7 +222,9 @@ public
       case StrongComponent.ALGEBRAIC_LOOP() algorithm
         strict := differentiateTearing(comp.strict, diffArguments_ptr, idx, context, name);
         casual := Util.applyOption(comp.casual, function differentiateTearing(diffArguments_ptr=diffArguments_ptr, idx=idx, context=context, name=name));
-      then StrongComponent.ALGEBRAIC_LOOP(-1, strict, casual, comp.linear, false, comp.homotopy, comp.status);
+        // if we differentiate for jacobian, the algebraic loops will always be linear
+        linear := match Pointer.access(diffArguments_ptr) case DIFFERENTIATION_ARGUMENTS(diffType = NBDifferentiate.DifferentiationType.JACOBIAN) then true; else comp.linear; end match;
+      then StrongComponent.ALGEBRAIC_LOOP(-1, strict, casual, linear, false, comp.homotopy, comp.status);
 
       case StrongComponent.ENTWINED_COMPONENT() algorithm
         Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " not implemented for entwined equation:\n" + StrongComponent.toString(comp)});

OMCompiler/Compiler/NSimCode/NSimCode.mo

@@ -406,15 +406,11 @@ public
 
             (linearLoops, nonlinearLoops, jacobians, simCodeIndices) := collectAlgebraicLoops(init, init_0, ode, algebraic, daeModeData, simCodeIndices, simcode_map);
 
-            for jac in jacobians loop
-              if Util.isSome(jac.jac_map) then
-                vars := SimVars.addSeedAndJacobianVars(vars, UnorderedMap.toList(Util.getOption(jac.jac_map)));
-              end if;
-            end for;
-
-            (jacA, simCodeIndices) := SimJacobian.createSimulationJacobian(bdae.ode, bdae.ode_event, simCodeIndices, simcode_map);
             if isSome(daeModeData) then
+              (jacA, simCodeIndices) := SimJacobian.createSimulationJacobian(Util.getOption(bdae.dae), simCodeIndices, simcode_map);
               daeModeData := DaeModeData.addJacobian(daeModeData, jacA);
+            else
+              (jacA, simCodeIndices) := SimJacobian.createSimulationJacobian(listAppend(bdae.ode, bdae.ode_event), simCodeIndices, simcode_map);
             end if;
 
             (jacB, simCodeIndices) := SimJacobian.empty("B", simCodeIndices);
@@ -424,6 +420,13 @@ public
             (jacH, simCodeIndices) := SimJacobian.empty("H", simCodeIndices);
             //jacobians := jacA :: jacB :: jacC :: jacD :: jacF :: jacobians;
             jacobians := listReverse(jacH :: jacF :: jacD :: jacC :: jacB :: jacA :: jacobians);
+
+            for jac in jacobians loop
+              if Util.isSome(jac.jac_map) then
+                vars := SimVars.addSeedAndJacobianVars(vars, UnorderedMap.toList(Util.getOption(jac.jac_map)));
+              end if;
+            end for;
+
             // jacobian blocks only from simulation jacobians
             jac_blocks := SimJacobian.getJacobiansBlocks({jacA, jacB, jacC, jacD, jacF, jacH});
             (jac_blocks, simCodeIndices) := SimStrongComponent.Block.fixIndices(jac_blocks, {}, simCodeIndices);

OMCompiler/Compiler/NSimCode/NSimJacobian.mo

@@ -276,6 +276,16 @@ public
                   UnorderedMap.add(cref, var.index, idx_map);
                 end if;
               end for;
+
+              // also add residuals if its DAE Mode
+              if jacobian.jacType == NBJacobian.JacobianType.DAE then
+                for var in resVars loop
+                  cref := SimVar.getName(var);
+                  UnorderedMap.add(cref, var.index, idx_map);
+                  //cref := BVariable.getPartnerCref(cref, BVariable.getVarPDer);
+                  //UnorderedMap.add(cref, var.index, idx_map);
+                end for;
+              end if;
             else
               for var in seedVars loop
                 cref := SimVar.getName(var);
@@ -325,13 +335,11 @@ public
     end create;
 
     function createSimulationJacobian
-      input list<Partition.Partition> ode;
-      input list<Partition.Partition> ode_event;
+      input list<Partition.Partition> partitions;
       output SimJacobian simJac;
       input output SimCode.SimCodeIndices simCodeIndices;
       input UnorderedMap<ComponentRef, SimVar> simcode_map;
     protected
-      list<Partition.Partition> partitions = listAppend(ode, ode_event);
       list<BackendDAE> jacobians = {};
       BackendDAE simJacobian;
       Option<SimJacobian> simJac_opt;

OMCompiler/Compiler/SimCode/SerializeModelInfo.mo

@@ -1267,6 +1267,10 @@ algorithm
       equation
         File.write(file,"iteration variable for solving an algebraic loop");
       then ();
+    case BackendDAE.DAE_RESIDUAL_VAR()
+      equation
+        File.write(file,"residual variable for dae mode");
+      then ();
     else
       equation
         Error.addMessage(Error.INTERNAL_ERROR, {"serializeVarKind failed for " + SimCodeUtil.simVarString(var)});