@@ -606,6 +606,153 @@ function generateAlgebraicLoopSolvers(id: string, dae: ModelicaDAE, result: FmuR
606606 return lines . join ( "\n" ) ;
607607}
608608
609+ /**
610+ * Generate _initializeSolve() — Newton-Raphson solver for initial equations
611+ * using exact AD Jacobians from StaticTapeBuilder.
612+ */
613+ function generateInitializeSolve ( id : string , dae : ModelicaDAE , result : FmuResult ) : string [ ] {
614+ const lines : string [ ] = [ ] ;
615+
616+ if ( dae . initialEquations . length === 0 ) {
617+ lines . push ( `void ${ id } ${ id } ) ;
618+ return lines ;
619+ }
620+
621+ // Collect initial equations that are simple equations (LHS = RHS)
622+ const initEqs : { lhs : ModelicaExpression ; rhs : ModelicaExpression } [ ] = [ ] ;
623+ for ( const eq of dae . initialEquations ) {
624+ if ( "expression1" in eq && "expression2" in eq ) {
625+ const se = eq as { expression1 : ModelicaExpression ; expression2 : ModelicaExpression } ;
626+ initEqs . push ( { lhs : se . expression1 , rhs : se . expression2 } ) ;
627+ }
628+ }
629+
630+ if ( initEqs . length === 0 ) {
631+ lines . push ( `void ${ id } ${ id } ) ;
632+ return lines ;
633+ }
634+
635+ const N = initEqs . length ;
636+
637+ // Identify unknowns: variables referenced in initial equations that are not parameters
638+ const paramNames = new Set < string > ( ) ;
639+ for ( const v of dae . variables ) {
640+ if ( v . variability === ModelicaVariability . PARAMETER || v . variability === ModelicaVariability . CONSTANT ) {
641+ paramNames . add ( v . name ) ;
642+ }
643+ }
644+
645+ const referencedNames = new Set < string > ( ) ;
646+ for ( const eq of initEqs ) {
647+ collectReferencedNames ( eq . lhs , referencedNames ) ;
648+ collectReferencedNames ( eq . rhs , referencedNames ) ;
649+ }
650+
651+ const unknowns : string [ ] = [ ] ;
652+ for ( const name of referencedNames ) {
653+ if ( ! paramNames . has ( name ) && name !== "time" ) {
654+ const sv = result . scalarVariables . find ( ( v ) => v . name === name ) ;
655+ if ( sv ) unknowns . push ( name ) ;
656+ }
657+ }
658+
659+ const nUnknowns = Math . min ( unknowns . length , N ) ;
660+ if ( nUnknowns === 0 ) {
661+ lines . push ( `void ${ id } ${ id } ) ;
662+ return lines ;
663+ }
664+
665+ // Build AD tapes for each residual: R_i = LHS_i - RHS_i
666+ const tapes : StaticTapeBuilder [ ] = [ ] ;
667+ const residualOutputIndices : number [ ] = [ ] ;
668+ for ( const eq of initEqs ) {
669+ const tape = new StaticTapeBuilder ( ) ;
670+ const lhsIdx = tape . walk ( eq . lhs ) ;
671+ const rhsIdx = tape . walk ( eq . rhs ) ;
672+ residualOutputIndices . push ( tape . pushOp ( { type : "sub" , a : lhsIdx , b : rhsIdx } ) ) ;
673+ tapes . push ( tape ) ;
674+ }
675+
676+ // Variable name -> inst->vars[VR] resolver
677+ const varResolver = ( name : string ) : string => {
678+ if ( name === "time" ) return "inst->time" ;
679+ const sv = result . scalarVariables . find ( ( v ) => v . name === name ) ;
680+ return sv ? `inst->vars[${ sv . valueReference } : `0.0 /* ${ name } ;
681+ } ;
682+
683+ lines . push ( `/* Initial Equation Solver with Exact AD Jacobian */` ) ;
684+ lines . push ( `void ${ id } ${ id } ) ;
685+ lines . push ( ` double R[${ N } ${ N * nUnknowns } ${ nUnknowns } ) ;
686+ lines . push ( ` int iter;` ) ;
687+ lines . push ( ` for (iter = 0; iter < 50; iter++) {` ) ;
688+
689+ // Forward pass: compute residuals R[i]
690+ for ( let i = 0 ; i < N ; i ++ ) {
691+ const tape = tapes [ i ] ;
692+ if ( ! tape ) continue ;
693+ const outIdx = residualOutputIndices [ i ] ;
694+ if ( outIdx === undefined ) continue ;
695+
696+ lines . push ( ` { /* Residual ${ i } ) ;
697+ const fwdCode = tape . emitForwardC ( varResolver ) ;
698+ lines . push ( ...fwdCode . map ( ( c ) => " " + c ) ) ;
699+ lines . push ( ` R[${ i } ${ outIdx } ) ;
700+ lines . push ( ` }` ) ;
701+ }
702+
703+ // Convergence check
704+ lines . push ( ` double err = 0.0;` ) ;
705+ lines . push ( ` for (int i = 0; i < ${ N } ) ;
706+ lines . push ( ` if (err < 1e-10) break;` ) ;
707+
708+ // Reverse pass: compute Jacobian J[row * nUnknowns + col]
709+ lines . push ( ` /* Exact Jacobian via Reverse-Mode AD */` ) ;
710+ for ( let row = 0 ; row < N ; row ++ ) {
711+ const tape = tapes [ row ] ;
712+ if ( ! tape ) continue ;
713+ const outIdx = residualOutputIndices [ row ] ;
714+ if ( outIdx === undefined ) continue ;
715+
716+ lines . push ( ` { /* J row ${ row } ) ;
717+ const fwdCode = tape . emitForwardC ( varResolver ) ;
718+ lines . push ( ...fwdCode . map ( ( c ) => " " + c ) ) ;
719+ const { code : revCode , gradients } = tape . emitReverseC ( outIdx ) ;
720+ lines . push ( ...revCode . map ( ( c ) => " " + c ) ) ;
721+
722+ for ( let col = 0 ; col < nUnknowns ; col ++ ) {
723+ const varName = unknowns [ col ] ;
724+ if ( ! varName ) continue ;
725+ const gradIdx = gradients . get ( varName ) ;
726+ if ( gradIdx !== undefined ) {
727+ lines . push ( ` J[${ row * nUnknowns + col } ${ gradIdx } ) ;
728+ } else {
729+ lines . push ( ` J[${ row * nUnknowns + col } ) ;
730+ }
731+ }
732+ lines . push ( ` }` ) ;
733+ }
734+
735+ // Newton update: solve J * dx = -R
736+ lines . push ( ` for (int i = 0; i < ${ N } ) ;
737+ lines . push ( ` solve_linear_sys(inst, ${ nUnknowns } ) ;
738+
739+ // Apply dx to unknowns
740+ for ( let j = 0 ; j < nUnknowns ; j ++ ) {
741+ const varName = unknowns [ j ] ;
742+ if ( ! varName ) continue ;
743+ const sv = result . scalarVariables . find ( ( v ) => v . name === varName ) ;
744+ if ( sv ) {
745+ lines . push ( ` inst->vars[${ sv . valueReference } ${ j } ${ varName } ) ;
746+ }
747+ }
748+
749+ lines . push ( ` }` ) ;
750+ lines . push ( `}` ) ;
751+ lines . push ( `` ) ;
752+
753+ return lines ;
754+ }
755+
609756function generateModelC ( id : string , dae : ModelicaDAE , result : FmuResult ) : string {
610757 // Reset buffer counters for this codegen invocation
611758 delayBufferCounter = 0 ;
@@ -762,9 +909,16 @@ function generateModelC(id: string, dae: ModelicaDAE, result: FmuResult): string
762909 }
763910 }
764911 }
912+ // Call the initial equation solver (generated below) after start values
913+ if ( dae . initialEquations . length > 0 ) {
914+ lines . push ( ` ${ id } ) ;
915+ }
765916 lines . push ( "}" ) ;
766917 lines . push ( "" ) ;
767918
919+ // ── Initial equation solver via Newton-Raphson with exact AD Jacobian ──
920+ lines . push ( ...generateInitializeSolve ( id , dae , result ) ) ;
921+
768922 // ── getDerivatives function ──
769923 lines . push ( `void ${ id } ${ id } ) ;
770924 lines . push ( ` ${ id } ) ;
0 commit comments