diff --git a/blast/src/blFluxes.cpp b/blast/src/blFluxes.cpp index 8b974ba90330a5bebac1b8a533a4b18c2bcef64d..0ada3cbfc75a978f5086b6d6515a0dc06865132f 100644 --- a/blast/src/blFluxes.cpp +++ b/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; diff --git a/blast/src/blFluxes.h b/blast/src/blFluxes.h index dee71ab5b46e1ad4101ed88136bf2cb1c33ca83b..cf7a54492a2f7ef89ba667370926dcf1ad32c46b 100644 --- a/blast/src/blFluxes.h +++ b/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