[NB] find homotopy calls in algebraic loops (#13754)

kabdelhak · web-flow · commit 6dd49e61d281 · 2025-03-28T10:10:18.000+01:00
* [NB] find homotopy calls in algebraic loops

* [NSim] collect algebraic loops and jacobians from lambda_0 systems
diff --git a/OMCompiler/Compiler/NBackEnd/Classes/NBStrongComponent.mo b/OMCompiler/Compiler/NBackEnd/Classes/NBStrongComponent.mo
@@ -56,6 +56,7 @@ protected
   import Causalize = NBCausalize;
   import BVariable = NBVariable;
   import NBEquation.{Equation, EquationPointer, EquationPointers, EquationAttributes, Iterator};
+  import Initialization = NBInitialization;
   import NBJacobian.JacobianType;
   import Matching = NBMatching;
   import Resizable = NBResizable;
@@ -153,6 +154,7 @@ public
     Option<Tearing> casual;
     Boolean linear              "true if the loop is linear";
     Boolean mixed               "true for systems that have discrete variables";
+    Boolean homotopy            "true if contains homotopy()";
     Solve.Status status;
   end ALGEBRAIC_LOOP;
 
@@ -210,7 +212,7 @@ public
       then str;
 
       case ALGEBRAIC_LOOP() algorithm
-        str := StringUtil.headline_3("BLOCK" + indexStr + ": Algebraic Loop (Linear = " + boolString(comp.linear) + ", Mixed = " + boolString(comp.mixed) + ")");
+        str := StringUtil.headline_3("BLOCK" + indexStr + ": Algebraic Loop (Linear = " + boolString(comp.linear) + ", Mixed = " + boolString(comp.mixed) + ", Homotopy = " + boolString(comp.homotopy) + ")");
         str := str + Tearing.toString(comp.strict, "Strict Tearing Set");
         if isSome(comp.casual) then
           str := str + Tearing.toString(Util.getOption(comp.casual), "Casual Tearing Set");
@@ -771,7 +773,7 @@ public
       case RESIZABLE_COMPONENT()then Equation.isDiscrete(Slice.getT(comp.eqn));
       case ENTWINED_COMPONENT() then List.all(list(isDiscrete(c) for c in comp.entwined_slices), bool_ident);
       case GENERIC_COMPONENT()  then Equation.isDiscrete(Slice.getT(comp.eqn));
-      case ALGEBRAIC_LOOP()     then not comp.mixed;
+      case ALGEBRAIC_LOOP()     then comp.mixed;
       case ALIAS()              then isDiscrete(comp.original);
       else algorithm
         Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " failed because of wrong component: " + toString(comp)});
@@ -818,6 +820,7 @@ public
         Tearing tearingSet;
         Slice<VariablePointer> var_slice;
         Slice<EquationPointer> eqn_slice;
+        Pointer<Boolean> homotopy = Pointer.create(false);
 
       // Size 1 strong component
       // - case 1: sliced equation because of for-equation
@@ -880,12 +883,16 @@ public
               residual_eqns   = comp_eqns,
               innerEquations  = listArray({}),
               jac             = NONE());
+            for eqn in comp_eqns loop
+              Equation.map(Pointer.access(Slice.getT(eqn)), function Initialization.containsHomotopyCall(b = homotopy));
+            end for;
           then ALGEBRAIC_LOOP(
             idx     = -1,
             strict  = tearingSet,
             casual  = NONE(),
             linear  = false,
             mixed   = false,
+            homotopy = Pointer.access(homotopy),
             status  = NBSolve.Status.IMPLICIT);
         end match;
       then comp;
diff --git a/OMCompiler/Compiler/NBackEnd/Modules/1_Main/NBInitialization.mo b/OMCompiler/Compiler/NBackEnd/Modules/1_Main/NBInitialization.mo
@@ -632,24 +632,29 @@ public
     input Boolean init "if init then replace with simplified, else replace with actual";
     input Pointer<Boolean> hasHom   "output, determines if partition contains homotopy()";
   algorithm
-    exp := match exp
-      local
-        Expression e;
-        String name;
-        Call call;
-      case Expression.CALL(call = call as Call.TYPED_CALL()) algorithm
-        name := AbsynUtil.pathString(Function.nameConsiderBuiltin(call.fn));
-        e := match name
-          case "homotopy" algorithm
-            Pointer.update(hasHom, true);
-          then listGet(Call.arguments(exp.call), if init then 2 else 1);
-          else exp;
-        end match;
-      then e;
-      else exp;
-    end match;
+    if Expression.isCallNamed(exp, "homotopy") then
+      exp := match exp
+        local
+          Call call;
+
+        case Expression.CALL(call = call as Call.TYPED_CALL()) algorithm
+          Pointer.update(hasHom, true);
+        then listGet(Call.arguments(exp.call), if init then 2 else 1);
+
+        else exp;
+      end match;
+    end if;
   end cleanupHomotopy;
 
+  function containsHomotopyCall
+    input output Expression exp;
+    input Pointer<Boolean> b;
+  algorithm
+    if not Pointer.access(b) then
+      Pointer.update(b, Expression.isCallNamed(exp, "homotopy"));
+    end if;
+  end containsHomotopyCall;
+
   function removeWhenEquation
     "this function checks if an equation has to be removed before initialization.
     true for: when branch without condition initial()"
diff --git a/OMCompiler/Compiler/NBackEnd/Modules/3_Post/NBTearing.mo b/OMCompiler/Compiler/NBackEnd/Modules/3_Post/NBTearing.mo
@@ -55,10 +55,11 @@ protected
   // Backend imports
   import Adjacency = NBAdjacency;
   import NBAdjacency.Solvability;
+  import Causalize = NBCausalize;
   import BEquation = NBEquation;
+  import Initialization = NBInitialization;
   import BJacobian = NBJacobian;
   import BVariable = NBVariable;
-  import Causalize = NBCausalize;
   import Differentiate = NBDifferentiate;
   import Inline = NBInline;
   import Jacobian = NBackendDAE.BackendDAE;
@@ -163,36 +164,43 @@ public
     // dummy adjacency matrix, don't need it for implicit
     Adjacency.Matrix dummy = Adjacency.EMPTY(NBAdjacency.MatrixStrictness.FULL);
     StrongComponent new_comp;
+    Pointer<Boolean> homotopy = Pointer.create(false);
   algorithm
     (comp, dummy, funcTree, index) := match comp
       // create implicit equations
       case StrongComponent.SINGLE_COMPONENT() algorithm
+        Equation.map(Pointer.access(comp.eqn), function Initialization.containsHomotopyCall(b = homotopy));
         new_comp := StrongComponent.ALGEBRAIC_LOOP(
           idx     = index,
           strict  = singleImplicit(comp.var, comp.eqn),
           casual  = NONE(),
           linear  = false,
           mixed   = false,
+          homotopy = Pointer.access(homotopy),
           status  = NBSolve.Status.IMPLICIT);
       then finalize(new_comp, dummy, funcTree, index, VariablePointers.empty(), EquationPointers.empty(), Pointer.create(0), kind);
 
       case StrongComponent.MULTI_COMPONENT() algorithm
+        Equation.map(Pointer.access(Slice.getT(comp.eqn)), function Initialization.containsHomotopyCall(b = homotopy));
         new_comp := StrongComponent.ALGEBRAIC_LOOP(
           idx     = index,
           strict  = singleImplicit(Slice.getT(listHead(comp.vars)), Slice.getT(comp.eqn)), // this is wrong! need to take all vars
           casual  = NONE(),
           linear  = false,
           mixed   = false,
+          homotopy = Pointer.access(homotopy),
           status  = NBSolve.Status.IMPLICIT);
       then finalize(new_comp, dummy, funcTree, index, VariablePointers.empty(), EquationPointers.empty(), Pointer.create(0), kind);
 
       case StrongComponent.RESIZABLE_COMPONENT() algorithm
+        Equation.map(Pointer.access(Slice.getT(comp.eqn)), function Initialization.containsHomotopyCall(b = homotopy));
         new_comp := StrongComponent.ALGEBRAIC_LOOP(
           idx     = index,
           strict  = singleImplicit(Slice.getT(comp.var), Slice.getT(comp.eqn)),
           casual  = NONE(),
           linear  = false,
           mixed   = false,
+          homotopy = Pointer.access(homotopy),
           status  = NBSolve.Status.IMPLICIT);
       then finalize(new_comp, dummy, funcTree, index, VariablePointers.empty(), EquationPointers.empty(), Pointer.create(0), kind);
 
diff --git a/OMCompiler/Compiler/NSimCode/NSimCode.mo b/OMCompiler/Compiler/NSimCode/NSimCode.mo
@@ -409,7 +409,8 @@ public
               //(daeModeData, modelInfo, jacA, simcode_map, simCodeIndices) := DaeModeData.createSparsityJacobian(daeModeData, modelInfo, Util.getOption(bdae.dae), simcode_map, simCodeIndices);
             //else
 
-            (linearLoops, nonlinearLoops, jacobians, simCodeIndices) := collectAlgebraicLoops(init, ode, algebraic, daeModeData, simCodeIndices, simcode_map);
+            (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)));
@@ -577,6 +578,7 @@ public
       "Collects algebraic loops from all systems (ode, init, init_0, dae, ...).
       ToDo: Add other systems once implemented!"
       input list<SimStrongComponent.Block> init;
+      input list<SimStrongComponent.Block> init_0;
       input list<list<SimStrongComponent.Block>> ode;
       input list<list<SimStrongComponent.Block>> algebraic;
       input Option<DaeModeData> daeModeData;
@@ -588,7 +590,7 @@ public
     protected
       list<list<SimStrongComponent.Block>> dae_mode_blcks;
     algorithm
-      (linearLoops, nonlinearLoops, jacobians, simCodeIndices) := SimStrongComponent.Block.collectAlgebraicLoops({init}, linearLoops, nonlinearLoops, jacobians, simCodeIndices, simcode_map);
+      (linearLoops, nonlinearLoops, jacobians, simCodeIndices) := SimStrongComponent.Block.collectAlgebraicLoops({init, init_0}, linearLoops, nonlinearLoops, jacobians, simCodeIndices, simcode_map);
       (linearLoops, nonlinearLoops, jacobians, simCodeIndices) := SimStrongComponent.Block.collectAlgebraicLoops(ode, linearLoops, nonlinearLoops, jacobians, simCodeIndices, simcode_map);
       (linearLoops, nonlinearLoops, jacobians, simCodeIndices) := SimStrongComponent.Block.collectAlgebraicLoops(algebraic, linearLoops, nonlinearLoops, jacobians, simCodeIndices, simcode_map);
       if isSome(daeModeData) then
diff --git a/OMCompiler/Compiler/NSimCode/NSimStrongComponent.mo b/OMCompiler/Compiler/NSimCode/NSimStrongComponent.mo
@@ -728,7 +728,7 @@ public
             indexSystem   = simCodeIndices.nonlinearSystemIndex,
             size          = listLength(crefs),
             jacobian      = Pointer.create(jacobian),
-            homotopy      = false,
+            homotopy      = comp.homotopy,
             mixed         = comp.mixed,
             torn          = true
           );
diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/nonlinearSystem.c b/OMCompiler/SimulationRuntime/c/simulation/solver/nonlinearSystem.c
@@ -403,7 +403,7 @@ void initializeNonlinearSystemData(DATA *data, threadData_t *threadData, NONLINE
   nonlinsys->numberOfIterations = 0;
 
   /* check if residual function pointer are valid */
-  assertStreamPrint(threadData, ((0 != nonlinsys->residualFunc)) || ((nonlinsys->strictTearingFunctionCall != NULL) ? (0 != nonlinsys->strictTearingFunctionCall) : 0), "residual function pointer is invalid" );
+  assertStreamPrint(threadData, (nonlinsys->residualFunc != NULL) || (nonlinsys->strictTearingFunctionCall != NULL), "residual function pointer is invalid");
 
   /* check if analytical jacobian is created */
   if(nonlinsys->jacobianIndex != -1)
