Try computing and using the exciter Jacobian.

Author: nkoukpaizan

GridKit/CommonMath.hpp

@@ -27,8 +27,9 @@ namespace GridKit
     __attribute__((always_inline)) inline ScalarT sigmoid(const ScalarT x)
     {
       using RealT               = typename GridKit::ScalarTraits<ScalarT>::RealT;
-      static constexpr RealT MU = 240.0;
-      return ONE<RealT> / (ONE<RealT> + std::exp(-MU * x));
+      static constexpr RealT MU = 40000.0;
+      return (HALF<RealT> * MU * x) / (ONE<RealT> + std::abs(MU * x)) + HALF<RealT>;
+      //return ONE<RealT> / (ONE<RealT> + std::exp(-MU * x));
     }
 
     /**

GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp

@@ -100,11 +100,6 @@ namespace GridKit
         int evaluateResidual() override;
         int evaluateJacobian() override;
 
-        bool hasJacobian() override
-        {
-          return false;
-        }
-
         /// Get the `ComponentSignals` from this `Ieeet1`
         auto getSignals()
             -> ComponentSignals<ScalarT,
@@ -150,7 +145,6 @@ namespace GridKit
         ScalarT vUEL_{0};
         ScalarT vOEL_{0};
         ScalarT vS_{0};
-        ScalarT Ec_{0}; // "Compensated" terminal measurment
 
         /// Component signal extension
         ComponentSignals<ScalarT, IdxT, Ieeet1InternalVariables, Ieeet1ExternalVariables> signals_;

GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp

@@ -27,9 +27,6 @@ namespace GridKit
         // std::cout << "Evaluate Jacobian for Ieeet1..." << std::endl;
         // std::cout << "Jacobian evaluation is experimental!" << std::endl;
 
-// The following will fail to build without a smooth approximation for
-// the saturation
-#if 0
         J_.zeroMatrix(); 
 
         GridKit::Enzyme::Sparse::InternalJacobianWithSignal<GridKit::PhasorDynamics::Exciter::Ieeet1<ScalarT, IdxT>,
@@ -57,7 +54,7 @@ namespace GridKit
                                                               bus_->getVariableIndices(),
                                                               y_.data(),
                                                               yp_.data(),
-                                                              wb_.data(),
+                                                              (bus_->y()).data(),
                                                               ws_.data(),
                                                               J_);
 
@@ -74,7 +71,6 @@ namespace GridKit
                                                               wb_.data(),
                                                               ws_.data(),
                                                               J_);
-#endif
 
         return 0;
       }

GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp

@@ -185,18 +185,18 @@ namespace GridKit
         // Terminal Voltage
         ScalarT vreal = bus_->Vr();
         ScalarT vimag = bus_->Vi();
-        Ec_           = std::sqrt(vreal * vreal + vimag * vimag);
+        ScalarT Ec    = std::sqrt(vreal * vreal + vimag * vimag);
 
         // Derived from External initial values
         ScalarT vr  = Ke_ * efd0;
         ScalarT vfx = Kf_ / Tf_ * efd0;
         ScalarT vtr = Ke_ / Ka_ * efd0;
 
         // Vref (setpoint = terminal + error)
-        vref_ = Ec_ + vtr;
+        vref_ = Ec + vtr;
 
         // IVP for Internal Variables
-        y_[0] = Ec_;  // y0 - vts  - Sensed term volt
+        y_[0] = Ec;  // y0 - vts  - Sensed term volt
         y_[1] = vr;   // y1 - vr   - Voltage reg
         y_[2] = efd0; // y2 - efdp - Efd pre mult
         y_[3] = vfx;  // y3 - vfx  - Exciter feedback
@@ -259,7 +259,7 @@ namespace GridKit
         // Read E comp (terminal voltage, unless compensation impedance)
         ScalarT vreal = wb[0];
         ScalarT vimag = wb[1];
-        Ec_           = std::sqrt(vreal * vreal + vimag * vimag);
+        ScalarT Ec    = std::sqrt(vreal * vreal + vimag * vimag);
 
         // Read Internal Variables
         ScalarT vts  = y[0]; // y0 - Sensed term volt
@@ -286,13 +286,13 @@ namespace GridKit
         ScalarT vr_ind = Math::indicator(Vrmin_, Vrmax_, vr, func);
 
         // Internal Differential Equations
-        f[0] = -vts_dot + (Ec_ - vts) / Tr_;
+        f[0] = -vts_dot + (Ec - vts) / Tr_;
         f[1] = -vr_dot + vr_ind * func / Ta_;
         f[2] = -efdp_dot + (vr - ve - Ke_ * efdp) / Te_;
         f[3] = -vfx_dot + vf / Tf_;
 
         // Internal Algebraic Equations
-        f[4] = -vtr + vref_ + vUEL_ + vOEL_ + vS_ - vf - vts;
+        f[4] = -vts + vref_ + vUEL_ + vOEL_ + vS_ - vtr - vf;
         f[5] = -vf + (efdp * Kf_) / Tf_ - vfx;
         f[6] = -ve + ksat * efdp;
         f[7] = -efd + efdp + omega * efdp * Ispdlim_;