Skip to content
Snippets Groups Projects
Verified Commit edccc177 authored by Paul Dechamps's avatar Paul Dechamps :speech_balloon:
Browse files

(feat) Analytical gradients for adjoint 

Change potential and mesh derivatives to analytical. Also changed the residuals computation of the viscous solver to the analytical expression instead of solving the linear system. This commit reduces the computational cost of the coupled adjoint methodology by a factor 4, at least.
parent 3b935a0e
No related branches found
No related tags found
1 merge request!1BLASTER v1.0
Pipeline #48289 passed
Stage: build
    Stage: test
      Stage: validation
        Expand all Hide whitespace changes
        Inline Side-by-side
        blast/src/DCoupledAdjoint.cpp
        + 289
        − 316
        View file @ edccc177
        This diff is collapsed.
        blast/src/DCoupledAdjoint.h
        + 2
        − 11
        View file @ edccc177
        ......@@ -116,7 +116,8 @@ private:
        Eigen::VectorXd dCdp_x; ///< Partial derivative of the pressure drag coefficient with respect to the mesh coordinates
        Eigen::VectorXd dCdf_x; ///< Partial derivative of the friction drag coefficient with respect to the mesh coordinates
        fwk::Timers tms; //< internal timers
        fwk::Timers tms; //< internal timers
        fwk::Timers tms2; //< internal precise timers
        public:
        CoupledAdjoint(std::shared_ptr<dart::Adjoint> iAdjoint, std::shared_ptr<blast::Driver> vSolver);
        ......@@ -125,16 +126,6 @@ public:
        void run();
        void reset();
        void setdRdStar_M(std::vector<std::vector<double>> _dRdStar_M);
        void setdRdStar_v(std::vector<std::vector<double>> _dRdStar_v);
        void setdRdStar_dStar(std::vector<std::vector<double>> _dRdStar_dStar);
        void setdRdStar_x(std::vector<std::vector<double>> _dRdStar_x);
        void setdRblw_rho(std::vector<std::vector<double>> _dRblw_rho);
        void setdRblw_v(std::vector<std::vector<double>> _dRblw_v);
        void setdRblw_dStar(std::vector<std::vector<double>> _dRblw_ddStar);
        void setdRblw_x(std::vector<std::vector<double>> _dRblw_dx);
        private:
        void buildAdjointMatrix(std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor>*>> &matrices, Eigen::SparseMatrix<double, Eigen::RowMajor> &matrixAdjoint);
        ......
        blast/src/DFluxes.cpp
        + 47
        − 27
        View file @ edccc177
        ......@@ -61,11 +61,8 @@ VectorXd Fluxes::blLaws(size_t iPoint, BoundaryLayer *bl, std::vector<double> u)
        {
        size_t nVar = bl->getnVar();
        MatrixXd timeMatrix = MatrixXd::Zero(5, 5);
        VectorXd spaceVector = VectorXd::Zero(5);
        // Dissipation ratio.
        double dissipRatio;
        double dissipRatio = 0.;
        if (bl->name == "wake")
        dissipRatio = 0.9;
        else
        ......@@ -146,39 +143,62 @@ VectorXd Fluxes::blLaws(size_t iPoint, BoundaryLayer *bl, std::vector<double> u)
        double ctEqa = bl->ctEq[iPoint];
        double usa = bl->us[iPoint];
        //--------------------------------------------------------------------------------//
        // Integral boundary layer equations. //
        // Space part.
        spaceVector(0) = dt_dx + (2 + u[1] - Mea * Mea) * (u[0] / u[3]) * due_dx - cfa / 2;
        spaceVector(1) = u[0] * dHstar_dx + (2 * Hstar2a + Hstara * (1 - u[1])) * u[0] / u[3] * due_dx - 2 * cda + Hstara * cfa / 2;
        spaceVector(2) = dN_dx - ax;
        spaceVector(3) = due_dx - c * (u[1] * dt_dx + u[0] * dH_dx) - dueExt_dx + cExt * ddeltaStarExt_dx;
        // spaceVector(0) = dt_dx + (2 + u[1] - Mea * Mea) * (u[0] / u[3]) * due_dx - cfa / 2;
        // spaceVector(1) = u[0] * dHstar_dx + (2 * Hstar2a + Hstara * (1 - u[1])) * u[0] / u[3] * due_dx - 2 * cda + Hstara * cfa / 2;
        // spaceVector(2) = dN_dx - ax;
        // spaceVector(3) = due_dx - c * (u[1] * dt_dx + u[0] * dH_dx) - dueExt_dx + cExt * ddeltaStarExt_dx;
        // if (bl->regime[iPoint] == 1)
        // spaceVector(4) = (2 * deltaa / u[4]) * dct_dx - 5.6 * (ctEqa - u[4] * dissipRatio) - 2 * deltaa * (4 / (3 * u[1] * u[0]) * (cfa / 2 - (((Hka - 1) / (6.7 * Hka * dissipRatio)) * ((Hka - 1) / (6.7 * Hka * dissipRatio)))) - 1 / u[3] * due_dx);
        if (bl->regime[iPoint] == 1)
        spaceVector(4) = (2 * deltaa / u[4]) * dct_dx - 5.6 * (ctEqa - u[4] * dissipRatio) - 2 * deltaa * (4 / (3 * u[1] * u[0]) * (cfa / 2 - (((Hka - 1) / (6.7 * Hka * dissipRatio)) * ((Hka - 1) / (6.7 * Hka * dissipRatio)))) - 1 / u[3] * due_dx);
        // // Time part.
        // timeMatrix(0, 0) = u[1] / u[3];
        // timeMatrix(0, 1) = u[0] / u[3];
        // timeMatrix(0, 3) = u[0] * u[1] / (u[3] * u[3]);
        // Time part.
        timeMatrix(0, 0) = u[1] / u[3];
        timeMatrix(0, 1) = u[0] / u[3];
        timeMatrix(0, 3) = u[0] * u[1] / (u[3] * u[3]);
        // timeMatrix(1, 0) = (1 + u[1] * (1 - Hstara)) / u[3];
        // timeMatrix(1, 1) = (1 - Hstara) * u[0] / u[3];
        // timeMatrix(1, 3) = (2 - Hstara * u[1]) * u[0] / (u[3] * u[3]);
        timeMatrix(1, 0) = (1 + u[1] * (1 - Hstara)) / u[3];
        timeMatrix(1, 1) = (1 - Hstara) * u[0] / u[3];
        timeMatrix(1, 3) = (2 - Hstara * u[1]) * u[0] / (u[3] * u[3]);
        // timeMatrix(2, 2) = 1.;
        timeMatrix(2, 2) = 1.;
        // timeMatrix(3, 0) = -c * u[1];
        // timeMatrix(3, 1) = -c * u[0];
        // timeMatrix(3, 3) = 1.;
        timeMatrix(3, 0) = -c * u[1];
        timeMatrix(3, 1) = -c * u[0];
        timeMatrix(3, 3) = 1.;
        // if (bl->regime[iPoint] == 1)
        // {
        // timeMatrix(4, 3) = 2 * deltaa / (usa * u[3] * u[3]);
        // timeMatrix(4, 4) = 2 * deltaa / (u[3] * u[4] * usa);
        // }
        // else
        // timeMatrix(4, 4) = 1.0;
        //--------------------------------------------------------------------------------//
        if (bl->regime[iPoint] == 1)
        VectorXd result = VectorXd::Zero(5);
        if (bl->regime[iPoint] == 0)
        {
        timeMatrix(4, 3) = 2 * deltaa / (usa * u[3] * u[3]);
        timeMatrix(4, 4) = 2 * deltaa / (u[3] * u[4] * usa);
        result[0] = 0.5*(-2.0*u[0]*(1.0*u[1] - 2.0)*(c*(dH_dx*u[0] + dt_dx*u[1]) - cExt*ddeltaStarExt_dx + dueExt_dx - due_dx) + 1.0*(c*u[0]*u[1] + u[3])*(2.0*due_dx*u[0]*(2.*Hstar2a - Hstara*(u[1] - 1)) + 1.0*u[3]*(Hstara*cfa - 4*cda + 2.*dHstar_dx*u[0])) + 1.0*(2.*due_dx*u[0]*(-Mea*Mea + u[1] + 2) + u[3]*(-cfa + 2.*dt_dx))*(1.0*Hstara*c*u[0]*u[1] + 1.0*Hstara*u[3] - 2.0*c*u[0] - 1.0*u[3]))/(c*u[0]*u[1] + u[3]);
        result[1] = 0.5*(2.0*u[0]*u[1]*(u[1] - 1)*(c*(dH_dx*u[0] + dt_dx*u[1]) - cExt*ddeltaStarExt_dx + dueExt_dx - due_dx) - 1.0*u[1]*(c*u[0]*u[1] + u[3])*(2.*due_dx*u[0]*(2.*Hstar2a - Hstara*(u[1] - 1)) + u[3]*(Hstara*cfa - 4*cda + 2.*dHstar_dx*u[0])) + 1.0*(2.*due_dx*u[0]*(-Mea*Mea + u[1] + 2) + u[3]*(-cfa + 2.*dt_dx))*(-1.0*Hstara*c*u[0]*u[1]*u[1] - 1.0*Hstara*u[1]*u[3] + 2.0*c*u[0]*u[1] + 1.0*u[1]*u[3] + 1.0*u[3]))/(u[0]*(c*u[0]*u[1] + u[3]));
        result[2] = -1.0*ax + 1.0*dN_dx;
        result[3] = 1.0*u[3]*(-1.0*c*(dH_dx*u[0] + dt_dx*u[1]) + 0.5*c*(2.*due_dx*u[0]*(-Mea*Mea + u[1] + 2) + u[3]*(-cfa + 2.*dt_dx)) + 1.0*cExt*ddeltaStarExt_dx - 1.0*dueExt_dx + 1.0*due_dx)/(c*u[0]*u[1] + u[3]);
        result[4] = 0.;
        }
        else if (bl->regime[iPoint] == 1)
        {
        result[0] = (-u[0]*(u[1] - 2)*(c*(dH_dx*u[0] + dt_dx*u[1]) - cExt*ddeltaStarExt_dx + dueExt_dx - due_dx) + (c*u[0]*u[1] + u[3])*(2.*due_dx*u[0]*(2.*Hstar2a - Hstara*(u[1] - 1)) + u[3]*(Hstara*cfa - 4*cda + 2.*dHstar_dx*u[0]))/2. + (2.*due_dx*u[0]*(-Mea*Mea + u[1] + 2) + u[3]*(-cfa + 2.*dt_dx))*(Hstara*c*u[0]*u[1] + Hstara*u[3] - 2.*c*u[0] - u[3])/2.)/(c*u[0]*u[1] + u[3]);
        result[1] = (2.*u[0]*u[1]*(u[1] - 1)*(c*(dH_dx*u[0] + dt_dx*u[1]) - cExt*ddeltaStarExt_dx + dueExt_dx - due_dx) - u[1]*(c*u[0]*u[1] + u[3])*(2.*due_dx*u[0]*(2.*Hstar2a - Hstara*(u[1] - 1)) + u[3]*(Hstara*cfa - 4*cda + 2.*dHstar_dx*u[0])) + (2.*due_dx*u[0]*(-Mea*Mea + u[1] + 2) + u[3]*(-cfa + 2.*dt_dx))*(-Hstara*c*u[0]*u[1]*u[1] - Hstara*u[1]*u[3] + 2.*c*u[0]*u[1] + u[1]*u[3] + u[3]))/(2.*u[0]*(c*u[0]*u[1] + u[3]));
        result[2] = -ax + dN_dx;
        result[3] = u[3]*(-2.*c*(dH_dx*u[0] + dt_dx*u[1]) + c*(2.*due_dx*u[0]*(-Mea*Mea + u[1] + 2) + u[3]*(-cfa + 2.*dt_dx)) + 2.*cExt*ddeltaStarExt_dx - 2.*dueExt_dx + 2.*due_dx)/(2.*(c*u[0]*u[1] + u[3]));
        result[4] = (-Hka*Hka*c*deltaa*dissipRatio*dissipRatio*u[0]*u[1]*u[4]*(2.*due_dx*u[0]*(-Mea*Mea + u[1] + 2) + u[3]*(-cfa + 2.*dt_dx))/2. + Hka*Hka*deltaa*dissipRatio*dissipRatio*u[0]*u[1]*u[4]*(c*(dH_dx*u[0] + dt_dx*u[1]) - cExt*ddeltaStarExt_dx + dueExt_dx - due_dx) + usa*(c*u[0]*u[1] + u[3])*(6*Hka*Hka*dct_dx*deltaa*dissipRatio*dissipRatio*u[0]*u[1]*u[3] + 16.8*Hka*Hka*dissipRatio*dissipRatio*u[0]*u[1]*u[3]*u[4]*(-ctEqa + dissipRatio*u[4]) + 2.*deltaa*u[4]*(3*Hka*Hka*dissipRatio*dissipRatio*due_dx*u[0]*u[1] + u[3]*(-2.*Hka*Hka*cfa*dissipRatio*dissipRatio + 0.0891067052795723*(Hka - 1)*(Hka - 1))))/6)/(Hka*Hka*deltaa*dissipRatio*dissipRatio*u[0]*u[1]*(c*u[0]*u[1] + u[3]));
        }
        else
        timeMatrix(4, 4) = 1.0;
        return timeMatrix.partialPivLu().solve(spaceVector);
        return result;
        }
        /**
        ......
        0% or .
        You are about to add 0 people to the discussion. Proceed with caution.
        Finish editing this message first!
        Please register or to comment