@@ -29,6 +29,7 @@ import {
2929 ModelicaWhenEquation ,
3030} from "./dae.js" ;
3131import type { FmuOptions , FmuResult } from "./fmi.js" ;
32+ import { differentiateExpr , simplifyExpr } from "./symbolic-diff.js" ;
3233import { ModelicaBinaryOperator , ModelicaUnaryOperator , ModelicaVariability } from "./syntax.js" ;
3334
3435/** Generated C source files. */
@@ -963,9 +964,91 @@ function generateFmi2FunctionsC(
963964 lines . push ( " *state = (fmi2FMUstate)copy;" ) ;
964965 lines . push ( " return fmi2OK;" ) ;
965966 lines . push ( "}" ) ;
967+ // ── fmi2GetDirectionalDerivative ──
968+ // Compute Δż = J · Δz where J is the Jacobian ∂f/∂x
969+ // J[i][j] = ∂(derivative[i])/∂(state[j]) — precomputed symbolically
966970 lines . push (
967- "fmi2Status fmi2GetDirectionalDerivative(fmi2Component c, const fmi2ValueReference unknown[], size_t nUnknown, const fmi2ValueReference known[], size_t nKnown, const fmi2Real dvKnown[], fmi2Real dvUnknown[]) { (void)c; (void)unknown; (void)nUnknown; (void)known; (void)nKnown; (void)dvKnown; (void)dvUnknown; return fmi2Error; } " ,
971+ "fmi2Status fmi2GetDirectionalDerivative(fmi2Component c, const fmi2ValueReference unknown[], size_t nUnknown, const fmi2ValueReference known[], size_t nKnown, const fmi2Real dvKnown[], fmi2Real dvUnknown[]) {" ,
968972 ) ;
973+ lines . push ( " FMUInstance* inst = (FMUInstance*)c;" ) ;
974+
975+ // Extract derivative equations: der(x) = f(x, y, t)
976+ const derEquations : { stateName : string ; rhs : ModelicaExpression } [ ] = [ ] ;
977+ for ( const eq of dae . equations ) {
978+ if ( ! ( "expression1" in eq && "expression2" in eq ) ) continue ;
979+ const simpleEq = eq as { expression1 : ModelicaExpression ; expression2 : ModelicaExpression } ;
980+ const lhsDer = extractDerName ( simpleEq . expression1 ) ;
981+ const rhsDer = extractDerName ( simpleEq . expression2 ) ;
982+ if ( lhsDer ) {
983+ derEquations . push ( { stateName : lhsDer , rhs : simpleEq . expression2 } ) ;
984+ } else if ( rhsDer ) {
985+ derEquations . push ( { stateName : rhsDer , rhs : simpleEq . expression1 } ) ;
986+ }
987+ }
988+
989+ // Build VR→derivative index and VR→state variable name mappings
990+ const derVRs = result . modelStructure . derivatives ; // VRs of derivative variables
991+ const jacStateVRs : number [ ] = [ ] ; // VRs of state variables
992+ const stateNames : string [ ] = [ ] ; // Names of state variables
993+ for ( const derEq of derEquations ) {
994+ const sv = result . scalarVariables . find ( ( v ) => v . name === derEq . stateName ) ;
995+ if ( sv ) {
996+ jacStateVRs . push ( sv . valueReference ) ;
997+ stateNames . push ( derEq . stateName ) ;
998+ }
999+ }
1000+
1001+ if ( derEquations . length > 0 && jacStateVRs . length > 0 ) {
1002+ // Emit local variable aliases
1003+ lines . push ( ` double time = inst->model.time;` ) ;
1004+ for ( const sv of result . scalarVariables ) {
1005+ if ( sv . causality === "independent" ) continue ;
1006+ const cName = varToC ( sv . name ) ;
1007+ lines . push ( ` double ${ cName } ${ sv . valueReference } ) ;
1008+ }
1009+ lines . push ( "" ) ;
1010+ lines . push ( " /* Zero output */" ) ;
1011+ lines . push ( " for (size_t i = 0; i < nUnknown; i++) dvUnknown[i] = 0.0;" ) ;
1012+ lines . push ( "" ) ;
1013+ lines . push ( " /* Accumulate Jacobian-vector product: dvUnknown[i] += J[i][j] * dvKnown[j] */" ) ;
1014+ lines . push ( " for (size_t j = 0; j < nKnown; j++) {" ) ;
1015+ lines . push ( " for (size_t i = 0; i < nUnknown; i++) {" ) ;
1016+
1017+ // For each (derivative, state) pair, emit the symbolic Jacobian entry
1018+ // Use switch on unknown VR, then switch on known VR
1019+ lines . push ( " switch (unknown[i]) {" ) ;
1020+ for ( let di = 0 ; di < derEquations . length ; di ++ ) {
1021+ const derVR = derVRs [ di ] ;
1022+ const derEq = derEquations [ di ] ;
1023+ if ( derVR === undefined || ! derEq ) continue ;
1024+ lines . push ( ` case ${ derVR } ${ derEq . stateName } ) ;
1025+ lines . push ( " switch (known[j]) {" ) ;
1026+ for ( let si = 0 ; si < stateNames . length ; si ++ ) {
1027+ const jacVR = jacStateVRs [ si ] ;
1028+ const stateName = stateNames [ si ] ;
1029+ if ( jacVR === undefined || ! stateName ) continue ;
1030+ // Symbolically differentiate rhs w.r.t. state variable
1031+ const jacobianEntry = simplifyExpr ( differentiateExpr ( derEq . rhs , stateName ) ) ;
1032+ const jacobianC = exprToC ( jacobianEntry ) ;
1033+ lines . push (
1034+ ` case ${ jacVR } ${ jacobianC } ${ stateName } ,
1035+ ) ;
1036+ }
1037+ lines . push ( " default: break;" ) ;
1038+ lines . push ( " }" ) ;
1039+ lines . push ( " break;" ) ;
1040+ }
1041+ lines . push ( " default: break;" ) ;
1042+ lines . push ( " }" ) ;
1043+ lines . push ( " }" ) ;
1044+ lines . push ( " }" ) ;
1045+ } else {
1046+ lines . push ( " (void)inst; (void)unknown; (void)nUnknown; (void)known; (void)nKnown; (void)dvKnown;" ) ;
1047+ lines . push ( " for (size_t i = 0; i < nUnknown; i++) dvUnknown[i] = 0.0;" ) ;
1048+ }
1049+
1050+ lines . push ( " return fmi2OK;" ) ;
1051+ lines . push ( "}" ) ;
9691052 lines . push ( "" ) ;
9701053
9711054 return lines . join ( "\n" ) ;
0 commit comments