(revert) Back to finite-diff for jacobian computation

For now the tape recording does not work for adjoint and jacobian at the same time

blast/src/blFluxes.cpp

@@ -26,7 +26,7 @@ blVectorXd Fluxes::computeResiduals(size_t inod, std::shared_ptr<BoundaryLayer>
 }
 
 /**
- * @brief Compute Jacobian of the integral boundary layer equations using 
+ * @brief Compute Jacobian of the integral boundary layer equations using
  * automatic differentiation or finite differences.
  *
  * @param inod Marker id.
@@ -36,7 +36,7 @@ blVectorXd Fluxes::computeResiduals(size_t inod, std::shared_ptr<BoundaryLayer>
  * @return blMatrixXd
  */
 blMatrixXd Fluxes::computeJacobian(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
-                                 blVectorXd const &sysRes, bldouble eps)
+                                   blVectorXd const &sysRes, bldouble eps)
 {
     size_t nVar = bl->getnVar();
     blMatrixXd JacMatrix = blMatrixXd::Zero(nVar, nVar);
@@ -44,44 +44,53 @@ blMatrixXd Fluxes::computeJacobian(size_t inod, std::shared_ptr<BoundaryLayer> &
     for (size_t k = 0; k < nVar; ++k)
         sol[k] = bl->u[inod * nVar + k];
 
-    #ifdef USE_CODI
-        // Compute the Jacobian using automatic differentiation.
-        bltape &tape = bldouble::getTape();
-        tape.reset();
-        tape.setActive();
-        // Register inputs
-        for (auto &ui : sol)
-            tape.registerInput(ui);
-        // Compute residuals
-        blVectorXd res = blLaws(inod, bl, sol);
-        // Register outputs
-        for (auto& ri : res)
-            tape.registerOutput(ri);
-        // Desactivate tape
-        tape.setPassive();
-
-        // Form the Jacobian
-        for (size_t i = 0; i < 5; ++i)
-        {
-            res[i].gradient() = 1.0;
+    // #ifdef USE_CODI
+    //     // Compute the Jacobian using automatic differentiation.
+    //     bltape &tape = bldouble::getTape();
 
-            tape.evaluate();
-            for(size_t j = 0; j < 5; ++j) {
-                JacMatrix(i, j) = sol[j].getGradient();
-            }
-            tape.clearAdjoints();
-        }
-    #else
-        // Compute the Jacobian using finite differences.
-        bldouble varSave = 0.;
-        for (size_t iVar = 0; iVar < nVar; ++iVar)
-        {
-            varSave = sol[iVar];
-            sol[iVar] += eps;
-            JacMatrix.col(iVar) = (blLaws(inod, bl, sol) - sysRes) / eps;
-            sol[iVar] = varSave;
-        }
-    #endif
+    //     // The tape could be active for adjoint calculation
+    //     bool wasActive = true;
+    //     if (!tape.isActive())
+    //     {
+    //         tape.reset();
+    //         tape.setActive();
+    //         wasActive = false;
+    //     }
+    //     // Register inputs
+    //     for (auto &ui : sol)
+    //         tape.registerInput(ui);
+    //     // Compute residuals
+    //     blVectorXd res = blLaws(inod, bl, sol);
+    //     // Register outputs
+    //     for (auto &ri : res)
+    //         tape.registerOutput(ri);
+    //     // Desactivate tape
+    //     tape.setPassive();
+
+    //     // Form the Jacobian
+    //     for (size_t i = 0; i < 5; ++i)
+    //     {
+    //         res[i].gradient() = 1.0;
+
+    //         tape.evaluate();
+    //         for (size_t j = 0; j < 5; ++j)
+    //             JacMatrix(i, j) = sol[j].getGradient();
+    //         tape.clearAdjoints();
+    //     }
+    //     if (wasActive)
+    //         tape.setActive();
+    //     // tape.printStatistics();
+    // #else
+    // Compute the Jacobian using finite differences.
+    bldouble varSave = 0.;
+    for (size_t iVar = 0; iVar < nVar; ++iVar)
+    {
+        varSave = sol[iVar];
+        sol[iVar] += eps;
+        JacMatrix.col(iVar) = (blLaws(inod, bl, sol) - sysRes) / eps;
+        sol[iVar] = varSave;
+    }
+    // #endif
     return JacMatrix;
 }
 
@@ -94,7 +103,7 @@ blMatrixXd Fluxes::computeJacobian(size_t inod, std::shared_ptr<BoundaryLayer> &
  * @return blVectorXd
  */
 blVectorXd Fluxes::blLaws(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
-                        std::vector<bldouble> &u)
+                          std::vector<bldouble> &u)
 {
     size_t nVar = bl->getnVar();
     auto &currNode = bl->nodes[inod];
@@ -363,7 +372,7 @@ bldouble Fluxes::amplificationRoutine(bldouble Hk, bldouble Ret, bldouble theta)
 
         // Normal envelope amplification rate
         bldouble m = -0.05 + 2.7 * Hmi - 5.5 * Hmi * Hmi + 3. * Hmi * Hmi * Hmi +
-                   0.1 * std::exp(-20. * Hmi);
+                     0.1 * std::exp(-20. * Hmi);
         bldouble arg = 3.87 * Hmi - 2.52;
         bldouble dndRet = 0.028 * (Hk - 1.) - 0.0345 * std::exp(-(arg * arg));
         ax = (m * dndRet / theta) * Rfac;

blast/src/blFluxes.h

@@ -23,11 +23,11 @@ public:
     ~Fluxes() {}
     blVectorXd computeResiduals(size_t inod, std::shared_ptr<BoundaryLayer> &bl);
     blMatrixXd computeJacobian(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
-                                    blVectorXd const &sysRes, bldouble eps);
+                               blVectorXd const &sysRes, bldouble eps);
 
 private:
     blVectorXd blLaws(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
-                           std::vector<bldouble> &u);
+                      std::vector<bldouble> &u);
     bldouble amplificationRoutine(bldouble Hk, bldouble Ret, bldouble theta) const;
 };
 } // namespace blast