Prevent use of FlowReactor sensitivity analysis

Some implementation issues are fixed here, but the sensitivities currently do not match finite difference comparisons.
Cantera · Jun 14, 2023 · 1f5a8c7 · 1f5a8c7
1 parent 4cb4e18
commit 1f5a8c7
Show file tree

Hide file tree

Showing 5 changed files with 17 additions and 17 deletions.
diff --git a/include/cantera/zeroD/FlowReactor.h b/include/cantera/zeroD/FlowReactor.h
@@ -38,16 +38,7 @@ class FlowReactor : public IdealGasReactor
         throw NotImplementedError("FlowReactor::eval");
     }
 
-    /*!
-     * Evaluate the reactor governing equations. Called by ReactorNet::eval.
-     * @param[in] t time.
-     * @param[in] y solution vector, length neq()
-     * @param[in] ydot rate of change of solution vector, length neq()
-     * @param[in] params sensitivity parameter vector, length ReactorNet::nparams()
-     * @param[out] residual residuals vector, length neq()
-     */
-    void evalDae(double t, double* y, double* ydot, double* params,
-                 double* residual) override;
+    void evalDae(double t, double* y, double* ydot, double* residual) override;
 
     //! Given a vector of length neq(), mark which variables should be
     //! considered algebraic constraints

diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h
@@ -125,12 +125,9 @@ class Reactor : public ReactorBase
      * @param[in] t time.
      * @param[in] y solution vector, length neq()
      * @param[in] ydot rate of change of solution vector, length neq()
-     * @param[in] params sensitivity parameter vector, length ReactorNet::nparams()
      * @param[out] residual residuals vector, length neq()
      */
-    virtual void evalDae(double t, double* y, double* ydot, double* params,
-                         double* residual)
-    {
+    virtual void evalDae(double t, double* y, double* ydot, double* residual) {
         throw NotImplementedError("Reactor::evalDae");
     }
 

diff --git a/src/numerics/IdasIntegrator.cpp b/src/numerics/IdasIntegrator.cpp
@@ -307,6 +307,17 @@ void IdasIntegrator::initialize(double t0, FuncEval& func)
     if (flag != IDA_SUCCESS) {
         throw CanteraError("IdasIntegrator::initialize", "IDASetUserData failed.");
     }
+    if (func.nparams() > 0) {
+        throw CanteraError("IdasIntegrator::initialize", "Sensitivity analysis "
+                           "for DAE systems is not fully implemented");
+        sensInit(t0, func);
+        flag = IDASetSensParams(m_ida_mem, func.m_sens_params.data(),
+                                func.m_paramScales.data(), NULL);
+        if (flag != IDA_SUCCESS) {
+            throw CanteraError("IdasIntegrator::initialize",
+                               "IDASetSensParams failed.");
+        }
+    }
     applyOptions();
 }
 

diff --git a/src/zeroD/FlowReactor.cpp b/src/zeroD/FlowReactor.cpp
@@ -237,8 +237,7 @@ void FlowReactor::updateSurfaceState(double* y)
     }
 }
 
-void FlowReactor::evalDae(double time, double* y, double* ydot, double* params,
-                          double* residual)
+void FlowReactor::evalDae(double time, double* y, double* ydot, double* residual)
 {
     m_thermo->restoreState(m_state);
 

diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp
@@ -327,7 +327,9 @@ void ReactorNet::evalDae(double t, double* y, double* ydot, double* p, double* r
     m_time = t;
     updateState(y);
     for (size_t n = 0; n < m_reactors.size(); n++) {
-        m_reactors[n]->evalDae(t, y, ydot, p, residual);
+        m_reactors[n]->applySensitivity(p);
+        m_reactors[n]->evalDae(t, y, ydot, residual);
+        m_reactors[n]->resetSensitivity(p);
     }
     checkFinite("ydot", ydot, m_nv);
 }