Merge branch 'adjoint2' into dev

Analytical adjoint

blast/src/DCoupledAdjoint.h

@@ -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

@@ -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;
 }
 
 /**