diff --git a/blast/src/DBoundaryLayer.cpp b/blast/src/DBoundaryLayer.cpp index 148ed29c5834bfe3e8ba4b65b5e17afddb6e066c..f93d1b3f58faeecc9e58152b9d089485a2381280 100644 --- a/blast/src/DBoundaryLayer.cpp +++ b/blast/src/DBoundaryLayer.cpp @@ -11,12 +11,13 @@ using namespace blast; * * @param _xtrF Forced transition location (default=-1; free transition). */ -BoundaryLayer::BoundaryLayer(double _xtrF, std::string _name) { - xtrF = _xtrF; - name = _name; - xtr = 1.0; - closSolver = new Closures(); - transitionNode = 0; +BoundaryLayer::BoundaryLayer(double _xtrF, std::string _name) +{ + xtrF = _xtrF; + name = _name; + xtr = 1.0; + closSolver = new Closures(); + transitionNode = 0; } BoundaryLayer::~BoundaryLayer() { delete closSolver; } @@ -24,90 +25,100 @@ BoundaryLayer::~BoundaryLayer() { delete closSolver; } void BoundaryLayer::setMesh(std::vector<double> _x, std::vector<double> _y, std::vector<double> _z, double _chord, double _xMin, std::vector<int> const _rows, - std::vector<int> const _no) { - if (_x.size() != _y.size() || _x.size() != _z.size()) - throw std::runtime_error("Mesh coordinates are not consistent."); - nNodes = _x.size(); - nElms = nNodes - 1; - if (_rows.size() == 0) { - rows.resize(0); - for (size_t i = 0; i < nNodes; ++i) - rows.push_back(static_cast<int>(i)); - } else if (_rows.size() != nNodes) - throw std::runtime_error("Row vector is not consistent."); - else { - rows = _rows; - } - // Elements - if (_no.size() == 0) { - no.resize(0); - for (size_t i = 0; i < nElms; ++i) - no.push_back(static_cast<int>(i)); - } else if (_no.size() != nElms) - throw std::runtime_error("No vector is not consistent."); - else { - no = _no; - } - - x = _x; - y = _y; - z = _z; - chord = _chord; - xMin = _xMin; - - // Compute the local coordinate. - this->computeLocalCoordinates(); - this->reset(); + std::vector<int> const _no) +{ + if (_x.size() != _y.size() || _x.size() != _z.size()) + throw std::runtime_error("Mesh coordinates are not consistent."); + nNodes = _x.size(); + nElms = nNodes - 1; + if (_rows.size() == 0) + { + rows.resize(0); + for (size_t i = 0; i < nNodes; ++i) + rows.push_back(static_cast<int>(i)); + } + else if (_rows.size() != nNodes) + throw std::runtime_error("Row vector is not consistent."); + else + { + rows = _rows; + } + // Elements + if (_no.size() == 0) + { + no.resize(0); + for (size_t i = 0; i < nElms; ++i) + no.push_back(static_cast<int>(i)); + } + else if (_no.size() != nElms) + throw std::runtime_error("No vector is not consistent."); + else + { + no = _no; + } + + x = _x; + y = _y; + z = _z; + chord = _chord; + xMin = _xMin; + + // Compute the local coordinate. + this->computeLocalCoordinates(); + this->reset(); } -void BoundaryLayer::computeLocalCoordinates() { - loc.resize(nNodes, 0.); - alpha.resize(nElms, 0.); - dx.resize(nElms, 0.); - - for (size_t iPoint = 0; iPoint < nElms; ++iPoint) { - alpha[iPoint] = - std::atan2(y[iPoint + 1] - y[iPoint], x[iPoint + 1] - x[iPoint]); - // dx = sqrt(delta x ^2 + delta y ^2). - dx[iPoint] = - std::sqrt((x[iPoint + 1] - x[iPoint]) * (x[iPoint + 1] - x[iPoint]) + - (y[iPoint + 1] - y[iPoint]) * (y[iPoint + 1] - y[iPoint])); - loc[iPoint + 1] = loc[iPoint] + dx[iPoint]; - } - - // Compute the x/c coordinate. - if (chord <= 0.) - throw std::runtime_error("Chord length is zero or negative."); - xoc.resize(nNodes, 0.); - for (size_t iPoint = 0; iPoint < nNodes; ++iPoint) - xoc[iPoint] = (x[iPoint] - xMin) / chord; +void BoundaryLayer::computeLocalCoordinates() +{ + loc.resize(nNodes, 0.); + alpha.resize(nElms, 0.); + dx.resize(nElms, 0.); + + for (size_t iPoint = 0; iPoint < nElms; ++iPoint) + { + alpha[iPoint] = + std::atan2(y[iPoint + 1] - y[iPoint], x[iPoint + 1] - x[iPoint]); + // dx = sqrt(delta x ^2 + delta y ^2). + dx[iPoint] = + std::sqrt((x[iPoint + 1] - x[iPoint]) * (x[iPoint + 1] - x[iPoint]) + + (y[iPoint + 1] - y[iPoint]) * (y[iPoint + 1] - y[iPoint])); + loc[iPoint + 1] = loc[iPoint] + dx[iPoint]; + } + + // Compute the x/c coordinate. + if (chord <= 0.) + throw std::runtime_error("Chord length is zero or negative."); + xoc.resize(nNodes, 0.); + for (size_t iPoint = 0; iPoint < nNodes; ++iPoint) + xoc[iPoint] = (x[iPoint] - xMin) / chord; } -void BoundaryLayer::reset() { - u.resize(nNodes * getnVar(), 0.); - Ret.resize(nNodes, 0.); - cd.resize(nNodes, 0.); - cf.resize(nNodes, 0.); - Hstar.resize(nNodes, 0.); - Hstar2.resize(nNodes, 0.); - Hk.resize(nNodes, 0.); - ctEq.resize(nNodes, 0.); - us.resize(nNodes, 0.); - deltaStar.resize(nNodes, 0.); - delta.resize(nNodes, 0.); - deltaStar.resize(nNodes, 0.); - regime.resize(nNodes, 0); - blowingVelocity.resize(nElms, 0.); - - Me.resize(nNodes, 0.); - vt.resize(nNodes, 0.); - rhoe.resize(nNodes, 0.); - deltaStarExt.resize(nNodes, 0.); - vtExt.resize(nNodes, 0.); - xxExt.resize(nNodes, 0.); - - transitionNode = 0; - xtr = 1.0; +void BoundaryLayer::reset() +{ + u.resize(nNodes * getnVar(), 0.); + Ret.resize(nNodes, 0.); + cd.resize(nNodes, 0.); + cf.resize(nNodes, 0.); + Hstar.resize(nNodes, 0.); + Hstar2.resize(nNodes, 0.); + Hk.resize(nNodes, 0.); + ctEq.resize(nNodes, 0.); + us.resize(nNodes, 0.); + deltaStar.resize(nNodes, 0.); + delta.resize(nNodes, 0.); + deltaStar.resize(nNodes, 0.); + regime.resize(nNodes, 0); + blowingVelocity.resize(nElms, 0.); + + Me.resize(nNodes, 0.); + vt.resize(nNodes, 0.); + rhoe.resize(nNodes, 0.); + deltaStarExt.resize(nNodes, 0.); + vtExt.resize(nNodes, 0.); + xxExt.resize(nNodes, 0.); + + transitionNode = 0; + xtr = 1.0; } /** @@ -115,35 +126,37 @@ void BoundaryLayer::reset() { * * @param Re Freestream Reynolds number. */ -void BoundaryLayer::setStagnationBC(double Re) { - double dv0 = std::abs((vt[1] - vt[0]) / (loc[1] - loc[0])); - if (dv0 > 0) - u[0] = sqrt(0.075 / (Re * dv0)); // Theta - else - u[0] = 1e-6; // Theta - u[1] = 2.23; // H - u[2] = 0.; // N - u[3] = vt[0]; // ue - u[4] = 0.; // Ctau - - Ret[0] = u[3] * u[0] * Re; - if (Ret[0] < 1.) { - Ret[0] = 1.; - u[3] = Ret[0] / (u[0] * Re); - } - deltaStar[0] = u[0] * u[1]; - - std::vector<double> lamParam(8, 0.); - closSolver->laminarClosures(lamParam, u[0], u[1], Ret[0], Me[0], name); - - Hstar[0] = lamParam[0]; - Hstar2[0] = lamParam[1]; - Hk[0] = lamParam[2]; - cf[0] = lamParam[3]; - cd[0] = lamParam[4]; - delta[0] = lamParam[5]; - ctEq[0] = lamParam[6]; - us[0] = lamParam[7]; +void BoundaryLayer::setStagnationBC(double Re) +{ + double dv0 = std::abs((vt[1] - vt[0]) / (loc[1] - loc[0])); + if (dv0 > 0) + u[0] = sqrt(0.075 / (Re * dv0)); // Theta + else + u[0] = 1e-6; // Theta + u[1] = 2.23; // H + u[2] = 0.; // N + u[3] = vt[0]; // ue + u[4] = 0.; // Ctau + + Ret[0] = u[3] * u[0] * Re; + if (Ret[0] < 1.) + { + Ret[0] = 1.; + u[3] = Ret[0] / (u[0] * Re); + } + deltaStar[0] = u[0] * u[1]; + + std::vector<double> lamParam(8, 0.); + closSolver->laminarClosures(lamParam, u[0], u[1], Ret[0], Me[0], name); + + Hstar[0] = lamParam[0]; + Hstar2[0] = lamParam[1]; + Hk[0] = lamParam[2]; + cf[0] = lamParam[3]; + cd[0] = lamParam[4]; + delta[0] = lamParam[5]; + ctEq[0] = lamParam[6]; + us[0] = lamParam[7]; } /** @@ -154,116 +167,127 @@ void BoundaryLayer::setStagnationBC(double Re) { * @param LowerCond Variables at the last point on the lower side. */ void BoundaryLayer::setWakeBC(double Re, std::vector<double> UpperCond, - std::vector<double> LowerCond) { - if (name != "wake") - std::cout << "Warning: Imposing wake boundary condition on " << name - << std::endl; - u[0] = UpperCond[0] + LowerCond[0]; // (Theta_up+Theta_low). - u[1] = (UpperCond[0] * UpperCond[1] + LowerCond[0] * LowerCond[1]) / - u[0]; // ((Theta*H)_up+(Theta*H)_low)/(Theta_up+Theta_low). - u[2] = 9.; // Amplification factor. - u[3] = vt[0]; // Inviscid velocity. - u[4] = (UpperCond[0] * UpperCond[4] + LowerCond[0] * LowerCond[4]) / - u[0]; // ((Ct*Theta)_up+(Theta)_low)/(Ct*Theta_up+Theta_low). - - Ret[0] = u[3] * u[0] * Re; // Reynolds number based on the momentum thickness. - - if (Ret[0] < 1.) { - Ret[0] = 1.; - u[3] = Ret[0] / (u[0] * Re); - } - deltaStar[0] = u[0] * u[1]; - - // Laminar closures. - std::vector<double> lamParam(8, 0.); - closSolver->laminarClosures(lamParam, u[0], u[1], Ret[0], Me[0], name); - Hstar[0] = lamParam[0]; - Hstar2[0] = lamParam[1]; - Hk[0] = lamParam[2]; - cf[0] = lamParam[3]; - cd[0] = lamParam[4]; - delta[0] = lamParam[5]; - ctEq[0] = lamParam[6]; - us[0] = lamParam[7]; - - for (size_t k = 0; k < nNodes; ++k) - regime[k] = 1; + std::vector<double> LowerCond) +{ + if (name != "wake") + std::cout << "Warning: Imposing wake boundary condition on " << name + << std::endl; + u[0] = UpperCond[0] + LowerCond[0]; // (Theta_up+Theta_low). + u[1] = (UpperCond[0] * UpperCond[1] + LowerCond[0] * LowerCond[1]) / + u[0]; // ((Theta*H)_up+(Theta*H)_low)/(Theta_up+Theta_low). + u[2] = 9.; // Amplification factor. + u[3] = vt[0]; // Inviscid velocity. + u[4] = (UpperCond[0] * UpperCond[4] + LowerCond[0] * LowerCond[4]) / + u[0]; // ((Ct*Theta)_up+(Theta)_low)/(Ct*Theta_up+Theta_low). + + Ret[0] = u[3] * u[0] * Re; // Reynolds number based on the momentum thickness. + + if (Ret[0] < 1.) + { + Ret[0] = 1.; + u[3] = Ret[0] / (u[0] * Re); + } + deltaStar[0] = u[0] * u[1]; + + // Laminar closures. + std::vector<double> lamParam(8, 0.); + closSolver->laminarClosures(lamParam, u[0], u[1], Ret[0], Me[0], name); + Hstar[0] = lamParam[0]; + Hstar2[0] = lamParam[1]; + Hk[0] = lamParam[2]; + cf[0] = lamParam[3]; + cd[0] = lamParam[4]; + delta[0] = lamParam[5]; + ctEq[0] = lamParam[6]; + us[0] = lamParam[7]; + + for (size_t k = 0; k < nNodes; ++k) + regime[k] = 1; } /** * @brief Compute the blowing velocity in each station of the region. * */ -void BoundaryLayer::computeBlowingVelocity() { - blowingVelocity.resize(nNodes - 1, 0.); - for (size_t i = 1; i < nNodes; ++i) { - blowingVelocity[i - 1] = (rhoe[i] * vt[i] * deltaStar[i] - - rhoe[i - 1] * vt[i - 1] * deltaStar[i - 1]) / - (rhoe[i] * (loc[i] - loc[i - 1])); - if (std::isnan(blowingVelocity[i - 1])) { - if (rhoe[i] == 0.) - std::cout << "Density is zero at point " << i << std::endl; - if ((loc[i] - loc[i - 1]) == 0.) - std::cout << "Points " << i - 1 << " and " << i - << " are at the same position " << loc[i] - << ", resulting in an infinite gradient." << std::endl; - if (std::isnan(deltaStar[i])) - std::cout << "Displacement thickness is nan at position " << i - << std::endl; - if (std::isnan(deltaStar[i - 1])) - std::cout << "Displacement thickness is nan at position " << i - 1 - << std::endl; - throw std::runtime_error( - "NaN detected in blowing velocity computation (see reason above)"); +void BoundaryLayer::computeBlowingVelocity() +{ + blowingVelocity.resize(nNodes - 1, 0.); + for (size_t i = 1; i < nNodes; ++i) + { + blowingVelocity[i - 1] = (rhoe[i] * vt[i] * deltaStar[i] - + rhoe[i - 1] * vt[i - 1] * deltaStar[i - 1]) / + (rhoe[i] * (loc[i] - loc[i - 1])); + if (std::isnan(blowingVelocity[i - 1])) + { + if (rhoe[i] == 0.) + std::cout << "Density is zero at point " << i << std::endl; + if ((loc[i] - loc[i - 1]) == 0.) + std::cout << "Points " << i - 1 << " and " << i + << " are at the same position " << loc[i] + << ", resulting in an infinite gradient." << std::endl; + if (std::isnan(deltaStar[i])) + std::cout << "Displacement thickness is nan at position " << i + << std::endl; + if (std::isnan(deltaStar[i - 1])) + std::cout << "Displacement thickness is nan at position " << i - 1 + << std::endl; + throw std::runtime_error( + "NaN detected in blowing velocity computation (see reason above)"); + } } - } } -std::vector<std::vector<double>> BoundaryLayer::evalGradwrtDeltaStar() { - std::vector<std::vector<double>> gradVe = - std::vector<std::vector<double>>(nElms); - for (auto &vec : gradVe) - vec.resize(nNodes, 0.); - - for (size_t i = 1; i < nNodes; ++i) { - gradVe[i - 1][i] = vt[i] / (loc[i] - loc[i - 1]); - gradVe[i - 1][i - 1] = - -1 / rhoe[i] * rhoe[i - 1] * vt[i - 1] / (loc[i] - loc[i - 1]); - } - return gradVe; +std::vector<std::vector<double>> BoundaryLayer::evalGradwrtDeltaStar() +{ + std::vector<std::vector<double>> gradVe = + std::vector<std::vector<double>>(nElms); + for (auto &vec : gradVe) + vec.resize(nNodes, 0.); + + for (size_t i = 1; i < nNodes; ++i) + { + gradVe[i - 1][i] = vt[i] / (loc[i] - loc[i - 1]); + gradVe[i - 1][i - 1] = + -1 / rhoe[i] * rhoe[i - 1] * vt[i - 1] / (loc[i] - loc[i - 1]); + } + return gradVe; } -std::vector<std::vector<double>> BoundaryLayer::evalGradwrtVt() { - std::vector<std::vector<double>> gradVe = - std::vector<std::vector<double>>(nElms); - for (auto &vec : gradVe) - vec.resize(nNodes, 0.); - - for (size_t i = 1; i < nNodes; ++i) { - gradVe[i - 1][i] = deltaStar[i] / (loc[i] - loc[i - 1]); - gradVe[i - 1][i - 1] = - -rhoe[i - 1] / rhoe[i] * deltaStar[i - 1] / (loc[i] - loc[i - 1]); - } - return gradVe; +std::vector<std::vector<double>> BoundaryLayer::evalGradwrtVt() +{ + std::vector<std::vector<double>> gradVe = + std::vector<std::vector<double>>(nElms); + for (auto &vec : gradVe) + vec.resize(nNodes, 0.); + + for (size_t i = 1; i < nNodes; ++i) + { + gradVe[i - 1][i] = deltaStar[i] / (loc[i] - loc[i - 1]); + gradVe[i - 1][i - 1] = + -rhoe[i - 1] / rhoe[i] * deltaStar[i - 1] / (loc[i] - loc[i - 1]); + } + return gradVe; } -std::vector<std::vector<double>> BoundaryLayer::evalGradwrtRho() { - std::vector<std::vector<double>> gradVe = - std::vector<std::vector<double>>(nElms); - for (auto &vec : gradVe) - vec.resize(nNodes, 0.); - - for (size_t i = 1; i < nNodes; ++i) { - gradVe[i - 1][i] = - -1 / (rhoe[i] * rhoe[i]) * - ((rhoe[i] * vt[i] * deltaStar[i] - - rhoe[i - 1] * vt[i - 1] * deltaStar[i - 1]) / - (loc[i] - loc[i - 1])) + - 1 / rhoe[i] * vt[i] * deltaStar[i] / (loc[i] - loc[i - 1]); - gradVe[i - 1][i - 1] = - -vt[i - 1] * deltaStar[i - 1] / (rhoe[i] * (loc[i] - loc[i - 1])); - } - return gradVe; +std::vector<std::vector<double>> BoundaryLayer::evalGradwrtRho() +{ + std::vector<std::vector<double>> gradVe = + std::vector<std::vector<double>>(nElms); + for (auto &vec : gradVe) + vec.resize(nNodes, 0.); + + for (size_t i = 1; i < nNodes; ++i) + { + gradVe[i - 1][i] = + -1 / (rhoe[i] * rhoe[i]) * + ((rhoe[i] * vt[i] * deltaStar[i] - + rhoe[i - 1] * vt[i - 1] * deltaStar[i - 1]) / + (loc[i] - loc[i - 1])) + + 1 / rhoe[i] * vt[i] * deltaStar[i] / (loc[i] - loc[i - 1]); + gradVe[i - 1][i - 1] = + -vt[i - 1] * deltaStar[i - 1] / (rhoe[i] * (loc[i] - loc[i - 1])); + } + return gradVe; } /** @@ -274,104 +298,114 @@ std::vector<std::vector<double>> BoundaryLayer::evalGradwrtRho() { * @return std::vector<double> Gradient of the velocity with respect to the x * coordinate. */ -std::vector<std::vector<double>> BoundaryLayer::evalGradwrtX() { - std::vector<std::vector<double>> gradVe = - std::vector<std::vector<double>>(nElms); - for (auto &vec : gradVe) - vec.resize(nNodes, 0.); - - // Compute dVe/dloc - std::vector<std::vector<double>> gradVe_loc = - std::vector<std::vector<double>>(nElms); - for (auto &vec : gradVe_loc) - vec.resize(nNodes, 0.); - - for (size_t i = 1; i < nNodes; ++i) { - double val = 1 / rhoe[i] * - (rhoe[i] * vt[i] * deltaStar[i] - - rhoe[i - 1] * vt[i - 1] * deltaStar[i - 1]) / - ((loc[i] - loc[i - 1]) * (loc[i] - loc[i - 1])); - gradVe_loc[i - 1][i] = -val; - gradVe_loc[i - 1][i - 1] = val; - } - - // Compute dloc/dx - // Derivative wrt x_j = (x_j - x_j-1) / sqrt((x_j - x_j-1)^2 + (y_j - - // y_j-1)^2) Derivative wrt x_j-1 = - (x_j - x_j-1) / sqrt((x_j - x_j-1)^2 + - // (y_j - y_j-1)^2) = - dloc/dx_j For loc_i, there are two contributions of - // x_j except if j = i. loc_n = sqrt((xn-xn-1)^2) + sqrt((xn-1-xn-2)^2) + ... - // + sqrt((x1-x0)^2) - std::vector<std::vector<double>> gradloc_x = - std::vector<std::vector<double>>(nNodes); - for (auto &vec : gradloc_x) - vec.resize(nNodes, 0.); - - for (size_t j = 1; j < nNodes; ++j) { - gradloc_x[j][j] = - (x[j] - x[j - 1]) / std::sqrt((x[j] - x[j - 1]) * (x[j] - x[j - 1]) + - (y[j] - y[j - 1]) * (y[j] - y[j - 1])); - for (size_t i = 0; i < nNodes - j; ++i) - gradloc_x[j + i][j - 1] = gradloc_x[j - 1][j - 1] - gradloc_x[j][j]; - } - - // Matrix product of gradVe_loc and gradloc_x - for (size_t i = 0; i < nElms; ++i) { - for (size_t j = 0; j < nNodes; ++j) { - for (size_t k = 0; k < nNodes; ++k) - gradVe[i][j] += gradVe_loc[i][k] * gradloc_x[k][j]; +std::vector<std::vector<double>> BoundaryLayer::evalGradwrtX() +{ + std::vector<std::vector<double>> gradVe = + std::vector<std::vector<double>>(nElms); + for (auto &vec : gradVe) + vec.resize(nNodes, 0.); + + // Compute dVe/dloc + std::vector<std::vector<double>> gradVe_loc = + std::vector<std::vector<double>>(nElms); + for (auto &vec : gradVe_loc) + vec.resize(nNodes, 0.); + + for (size_t i = 1; i < nNodes; ++i) + { + double val = 1 / rhoe[i] * + (rhoe[i] * vt[i] * deltaStar[i] - + rhoe[i - 1] * vt[i - 1] * deltaStar[i - 1]) / + ((loc[i] - loc[i - 1]) * (loc[i] - loc[i - 1])); + gradVe_loc[i - 1][i] = -val; + gradVe_loc[i - 1][i - 1] = val; + } + + // Compute dloc/dx + // Derivative wrt x_j = (x_j - x_j-1) / sqrt((x_j - x_j-1)^2 + (y_j - + // y_j-1)^2) Derivative wrt x_j-1 = - (x_j - x_j-1) / sqrt((x_j - x_j-1)^2 + + // (y_j - y_j-1)^2) = - dloc/dx_j For loc_i, there are two contributions of + // x_j except if j = i. loc_n = sqrt((xn-xn-1)^2) + sqrt((xn-1-xn-2)^2) + ... + // + sqrt((x1-x0)^2) + std::vector<std::vector<double>> gradloc_x = + std::vector<std::vector<double>>(nNodes); + for (auto &vec : gradloc_x) + vec.resize(nNodes, 0.); + + for (size_t j = 1; j < nNodes; ++j) + { + gradloc_x[j][j] = + (x[j] - x[j - 1]) / std::sqrt((x[j] - x[j - 1]) * (x[j] - x[j - 1]) + + (y[j] - y[j - 1]) * (y[j] - y[j - 1])); + for (size_t i = 0; i < nNodes - j; ++i) + gradloc_x[j + i][j - 1] = gradloc_x[j - 1][j - 1] - gradloc_x[j][j]; + } + + // Matrix product of gradVe_loc and gradloc_x + for (size_t i = 0; i < nElms; ++i) + { + for (size_t j = 0; j < nNodes; ++j) + { + for (size_t k = 0; k < nNodes; ++k) + gradVe[i][j] += gradVe_loc[i][k] * gradloc_x[k][j]; + } } - } - return gradVe; + return gradVe; } -std::vector<std::vector<double>> BoundaryLayer::evalGradwrtY() { - std::vector<std::vector<double>> gradVe = - std::vector<std::vector<double>>(nElms); - for (auto &vec : gradVe) - vec.resize(nNodes, 0.); - - // Compute dVe/dloc - std::vector<std::vector<double>> gradVe_loc = - std::vector<std::vector<double>>(nElms); - for (auto &vec : gradVe_loc) - vec.resize(nNodes, 0.); - - for (size_t i = 1; i < nNodes; ++i) { - double val = 1 / rhoe[i] * - (rhoe[i] * vt[i] * deltaStar[i] - - rhoe[i - 1] * vt[i - 1] * deltaStar[i - 1]) / - ((loc[i] - loc[i - 1]) * (loc[i] - loc[i - 1])); - gradVe_loc[i - 1][i] = -val; - gradVe_loc[i - 1][i - 1] = val; - } - - // Compute dloc/dy - // Derivative wrt y_j = (y_j - y_j-1) / sqrt((y_j - y_j-1)^2 + (y_j - - // y_j-1)^2) Derivative wrt y_j-1 = - (y_j - y_j-1) / sqrt((y_j - y_j-1)^2 + - // (y_j - y_j-1)^2) = - dloc/dy_j For loc_i, there are two contributions of - // y_j except if j = i. loc_n = sqrt((yn-yn-1)^2) + sqrt((yn-1-yn-2)^2) + ... - // + sqrt((y1-y0)^2) - std::vector<std::vector<double>> gradloc_y = - std::vector<std::vector<double>>(nNodes); - for (auto &vec : gradloc_y) - vec.resize(nNodes, 0.); - - for (size_t j = 1; j < nNodes; ++j) { - gradloc_y[j][j] = - (y[j] - y[j - 1]) / std::sqrt((x[j] - x[j - 1]) * (x[j] - x[j - 1]) + - (y[j] - y[j - 1]) * (y[j] - y[j - 1])); - for (size_t i = 0; i < nNodes - j; ++i) - gradloc_y[j + i][j - 1] = gradloc_y[j - 1][j - 1] - gradloc_y[j][j]; - } - - // Matrix product of gradVe_loc and gradloc_x - for (size_t i = 0; i < nElms; ++i) { - for (size_t j = 0; j < nNodes; ++j) { - for (size_t k = 0; k < nNodes; ++k) - gradVe[i][j] += gradVe_loc[i][k] * gradloc_y[k][j]; +std::vector<std::vector<double>> BoundaryLayer::evalGradwrtY() +{ + std::vector<std::vector<double>> gradVe = + std::vector<std::vector<double>>(nElms); + for (auto &vec : gradVe) + vec.resize(nNodes, 0.); + + // Compute dVe/dloc + std::vector<std::vector<double>> gradVe_loc = + std::vector<std::vector<double>>(nElms); + for (auto &vec : gradVe_loc) + vec.resize(nNodes, 0.); + + for (size_t i = 1; i < nNodes; ++i) + { + double val = 1 / rhoe[i] * + (rhoe[i] * vt[i] * deltaStar[i] - + rhoe[i - 1] * vt[i - 1] * deltaStar[i - 1]) / + ((loc[i] - loc[i - 1]) * (loc[i] - loc[i - 1])); + gradVe_loc[i - 1][i] = -val; + gradVe_loc[i - 1][i - 1] = val; + } + + // Compute dloc/dy + // Derivative wrt y_j = (y_j - y_j-1) / sqrt((y_j - y_j-1)^2 + (y_j - + // y_j-1)^2) Derivative wrt y_j-1 = - (y_j - y_j-1) / sqrt((y_j - y_j-1)^2 + + // (y_j - y_j-1)^2) = - dloc/dy_j For loc_i, there are two contributions of + // y_j except if j = i. loc_n = sqrt((yn-yn-1)^2) + sqrt((yn-1-yn-2)^2) + ... + // + sqrt((y1-y0)^2) + std::vector<std::vector<double>> gradloc_y = + std::vector<std::vector<double>>(nNodes); + for (auto &vec : gradloc_y) + vec.resize(nNodes, 0.); + + for (size_t j = 1; j < nNodes; ++j) + { + gradloc_y[j][j] = + (y[j] - y[j - 1]) / std::sqrt((x[j] - x[j - 1]) * (x[j] - x[j - 1]) + + (y[j] - y[j - 1]) * (y[j] - y[j - 1])); + for (size_t i = 0; i < nNodes - j; ++i) + gradloc_y[j + i][j - 1] = gradloc_y[j - 1][j - 1] - gradloc_y[j][j]; } - } - return gradVe; + + // Matrix product of gradVe_loc and gradloc_x + for (size_t i = 0; i < nElms; ++i) + { + for (size_t j = 0; j < nNodes; ++j) + { + for (size_t k = 0; k < nNodes; ++k) + gradVe[i][j] += gradVe_loc[i][k] * gradloc_y[k][j]; + } + } + return gradVe; } /** @@ -380,23 +414,24 @@ std::vector<std::vector<double>> BoundaryLayer::evalGradwrtY() { * @param iPoint Node id. * @note Function used for debugging. */ -void BoundaryLayer::printSolution(size_t iPoint) const { - if (iPoint < 0 || iPoint > nNodes - 1) - throw std::runtime_error("Tried to access element outside of region size."); - std::cout << "Pt " << iPoint << "Reg " << name << " Turb " << regime[iPoint] - << std::endl; - std::cout << std::setprecision(5) << "x " << x[iPoint] << ", xoc " - << xoc[iPoint] << ", xx " << loc[iPoint] << "xxExt " - << xxExt[iPoint] << std::endl; - std::cout << std::scientific << std::setprecision(15); - std::cout << std::setw(10) << "T " << u[iPoint * nVar + 0] << "\n" - << std::setw(10) << "H " << u[iPoint * nVar + 1] << "\n" - << std::setw(10) << "N " << u[iPoint * nVar + 2] << "\n" - << std::setw(10) << "ue " << u[iPoint * nVar + 3] << "\n" - << std::setw(10) << "Ct " << u[iPoint * nVar + 4] << "\n" - << std::setw(10) << "dStar (BL) " << deltaStar[iPoint] << "\n" - << std::setw(10) << "dStar (ext) " << deltaStarExt[iPoint] << "\n" - << std::setw(10) << "vt " << vt[iPoint] << "\n" - << std::setw(10) << "Me " << Me[iPoint] << "\n" - << std::setw(10) << "rhoe " << rhoe[iPoint] << std::endl; +void BoundaryLayer::printSolution(size_t iPoint) const +{ + if (iPoint < 0 || iPoint > nNodes - 1) + throw std::runtime_error("Tried to access element outside of region size."); + std::cout << "Pt " << iPoint << "Reg " << name << " Turb " << regime[iPoint] + << std::endl; + std::cout << std::setprecision(5) << "x " << x[iPoint] << ", xoc " + << xoc[iPoint] << ", xx " << loc[iPoint] << "xxExt " + << xxExt[iPoint] << std::endl; + std::cout << std::scientific << std::setprecision(15); + std::cout << std::setw(10) << "T " << u[iPoint * nVar + 0] << "\n" + << std::setw(10) << "H " << u[iPoint * nVar + 1] << "\n" + << std::setw(10) << "N " << u[iPoint * nVar + 2] << "\n" + << std::setw(10) << "ue " << u[iPoint * nVar + 3] << "\n" + << std::setw(10) << "Ct " << u[iPoint * nVar + 4] << "\n" + << std::setw(10) << "dStar (BL) " << deltaStar[iPoint] << "\n" + << std::setw(10) << "dStar (ext) " << deltaStarExt[iPoint] << "\n" + << std::setw(10) << "vt " << vt[iPoint] << "\n" + << std::setw(10) << "Me " << Me[iPoint] << "\n" + << std::setw(10) << "rhoe " << rhoe[iPoint] << std::endl; } \ No newline at end of file diff --git a/blast/src/DBoundaryLayer.h b/blast/src/DBoundaryLayer.h index 87102059ad31837068fbb5396c675b511ff1115a..a4be6203544795e1ed2037c7107a1d217783fa45 100644 --- a/blast/src/DBoundaryLayer.h +++ b/blast/src/DBoundaryLayer.h @@ -7,153 +7,162 @@ #include <algorithm> -namespace blast { +namespace blast +{ /** * @brief Boundary layer region upper/lower side or wake. * @author Paul Dechamps */ -class BLAST_API BoundaryLayer : public fwk::wSharedObject { +class BLAST_API BoundaryLayer : public fwk::wSharedObject +{ private: - unsigned int const nVar = - 5; ///< Number of variables of the partial differential problem. - double nCrit = 9.0; ///< Critical amplification factor. + unsigned int const nVar = + 5; ///< Number of variables of the partial differential problem. + double nCrit = 9.0; ///< Critical amplification factor. public: - Closures *closSolver; ///< Closure relations class. - std::string name; ///< Name of the region. - - // Mesh - size_t nNodes; ///< Number of points in the domain. - size_t nElms; ///< Number of cells in the domain. - std::vector<double> x; ///< x coordinate in the global frame of reference. - std::vector<double> y; ///< y coordinate in the global frame of reference. - std::vector<double> z; ///< z coordinate in the global frame of reference. - std::vector<double> xoc; ///< x/c coordinate in the global frame of reference. - std::vector<double> loc; ///< xi coordinate in the local frame of reference. - std::vector<int> rows; ///< Reference numbers of the nodes. - std::vector<int> no; ///< Reference numbers of the nodes. - std::vector<double> dx; ///< Cell size. - std::vector<double> alpha; ///< Angle of the cell wrt the x axis of the global - ///< frame of reference. - double chord; ///< Chord of the section. - double xMin; ///< Minimum x coordinate of the mesh (for local coordinates - ///< computation). - - // Boundary layer variables. - std::vector<double> u; ///< Solution vector (theta, H, N, ue, Ct). - std::vector<double> - Ret; ///< Reynolds number based on the momentum thickness (theta). - std::vector<double> cd; ///< Local dissipation coefficient. - std::vector<double> cf; ///< Local friction coefficient. - std::vector<double> - Hstar; ///< Kinetic energy shape parameter (thetaStar/theta). - std::vector<double> - Hstar2; ///< Density shape parameter (deltaStarStar/theta). - std::vector<double> Hk; ///< Kinematic shape parameter (int(1-u/ue d_eta)). - std::vector<double> - ctEq; ///< Equilibrium shear stress coefficient (turbulent BL). - std::vector<double> us; ///< Equivalent normalized wall slip velocity. - std::vector<double> delta; ///< Boundary layer thickness. - std::vector<double> - deltaStar; ///< Displacement thickness (int(1-rho*u/rhoe*ue d_eta)). - - std::vector<double> Me; ///< Mach number. - std::vector<double> vt; ///< Velocity. - std::vector<double> rhoe; ///< Density. - std::vector<double> xxExt; ///< xi coordinate in the local frame of reference, - ///< at the edge of the boundary layer. - std::vector<double> deltaStarExt; ///< Displacement thickness. - std::vector<double> vtExt; ///< Previous velocity. - - std::vector<double> blowingVelocity; ///< Blowing velocity. - - // Transition related variables. - std::vector<int> regime; ///< Laminar (0) or turbulent (1) regime. - size_t transitionNode; ///< Node id of the transition location. - double xtrF; ///< Forced transition location in the local frame of reference - ///< (x/c). - double xtr; ///< Transition location in the local frame of reference (x/c). - - BoundaryLayer(double _xtrF = -1.0, std::string _name = "None"); - ~BoundaryLayer(); - - void reset(); - - // Getters - std::string getName() const { return name; }; - size_t getnNodes() const { return nNodes; }; - size_t getnElms() const { return nElms; }; - - std::vector<double> getDeltaStar() const { return deltaStar; }; - std::vector<double> getUe() const { - std::vector<double> ue(nNodes, 0.); - for (size_t i = 0; i < nNodes; ++i) - ue[i] = u[i * nVar + 3]; - return ue; - }; - std::vector<double> getBlowing() const { return blowingVelocity; }; - double getMaxMach() const { return *std::max_element(Me.begin(), Me.end()); }; - - // Setters - void setMesh(std::vector<double> const _x, std::vector<double> const y, - std::vector<double> const z, double const _chord, - double const _xMin, - std::vector<int> const _rows = std::vector<int>(0), - std::vector<int> const _no = std::vector<int>(0)); - void computeLocalCoordinates(); - void setMe(std::vector<double> const _Me) { - if (_Me.size() != nNodes) - throw std::runtime_error("Mach number vector is not consistent."); - Me = _Me; - }; - void setVt(std::vector<double> const _vt) { - if (_vt.size() != nNodes) - throw std::runtime_error("Velocity vector is not consistent."); - vt = _vt; - }; - void setRhoe(std::vector<double> const _rhoe) { - if (_rhoe.size() != nNodes) - throw std::runtime_error("Density vector is not consistent."); - rhoe = _rhoe; - }; - void setxxExt(std::vector<double> const _xxExt) { - if (_xxExt.size() != nNodes) - throw std::runtime_error("External mesh vector is not consistent."); - xxExt = _xxExt; - }; - void setDeltaStarExt(std::vector<double> const _deltaStarExt) { - if (_deltaStarExt.size() != nNodes) - throw std::runtime_error( - "Displacement thickness vector is not consistent."); - deltaStarExt = _deltaStarExt; - }; - void setVtExt(std::vector<double> const _vtExt) { - if (_vtExt.size() != nNodes) - throw std::runtime_error("Velocity vector is not consistent."); - vtExt = _vtExt; - }; - void setxtrF(double const _xtrF) { xtrF = _xtrF; }; - - // Boundary conditions. - void setStagnationBC(double Re); - void setWakeBC(double Re, std::vector<double> upperConditions, - std::vector<double> lowerConditions); - - // Getters. - size_t getnVar() const { return nVar; }; - double getnCrit() const { return nCrit; }; - void setnCrit(double const _nCrit) { nCrit = _nCrit; }; - - // Others - void printSolution(size_t iPoint) const; - void computeBlowingVelocity(); - std::vector<std::vector<double>> evalGradwrtRho(); - std::vector<std::vector<double>> evalGradwrtVt(); - std::vector<std::vector<double>> evalGradwrtDeltaStar(); - std::vector<std::vector<double>> evalGradwrtX(); - std::vector<std::vector<double>> evalGradwrtY(); + Closures *closSolver; ///< Closure relations class. + std::string name; ///< Name of the region. + + // Mesh + size_t nNodes; ///< Number of points in the domain. + size_t nElms; ///< Number of cells in the domain. + std::vector<double> x; ///< x coordinate in the global frame of reference. + std::vector<double> y; ///< y coordinate in the global frame of reference. + std::vector<double> z; ///< z coordinate in the global frame of reference. + std::vector<double> xoc; ///< x/c coordinate in the global frame of reference. + std::vector<double> loc; ///< xi coordinate in the local frame of reference. + std::vector<int> rows; ///< Reference numbers of the nodes. + std::vector<int> no; ///< Reference numbers of the nodes. + std::vector<double> dx; ///< Cell size. + std::vector<double> alpha; ///< Angle of the cell wrt the x axis of the global + ///< frame of reference. + double chord; ///< Chord of the section. + double xMin; ///< Minimum x coordinate of the mesh (for local coordinates + ///< computation). + + // Boundary layer variables. + std::vector<double> u; ///< Solution vector (theta, H, N, ue, Ct). + std::vector<double> + Ret; ///< Reynolds number based on the momentum thickness (theta). + std::vector<double> cd; ///< Local dissipation coefficient. + std::vector<double> cf; ///< Local friction coefficient. + std::vector<double> + Hstar; ///< Kinetic energy shape parameter (thetaStar/theta). + std::vector<double> + Hstar2; ///< Density shape parameter (deltaStarStar/theta). + std::vector<double> Hk; ///< Kinematic shape parameter (int(1-u/ue d_eta)). + std::vector<double> + ctEq; ///< Equilibrium shear stress coefficient (turbulent BL). + std::vector<double> us; ///< Equivalent normalized wall slip velocity. + std::vector<double> delta; ///< Boundary layer thickness. + std::vector<double> + deltaStar; ///< Displacement thickness (int(1-rho*u/rhoe*ue d_eta)). + + std::vector<double> Me; ///< Mach number. + std::vector<double> vt; ///< Velocity. + std::vector<double> rhoe; ///< Density. + std::vector<double> xxExt; ///< xi coordinate in the local frame of reference, + ///< at the edge of the boundary layer. + std::vector<double> deltaStarExt; ///< Displacement thickness. + std::vector<double> vtExt; ///< Previous velocity. + + std::vector<double> blowingVelocity; ///< Blowing velocity. + + // Transition related variables. + std::vector<int> regime; ///< Laminar (0) or turbulent (1) regime. + size_t transitionNode; ///< Node id of the transition location. + double xtrF; ///< Forced transition location in the local frame of reference + ///< (x/c). + double xtr; ///< Transition location in the local frame of reference (x/c). + + BoundaryLayer(double _xtrF = -1.0, std::string _name = "None"); + ~BoundaryLayer(); + + void reset(); + + // Getters + std::string getName() const { return name; }; + size_t getnNodes() const { return nNodes; }; + size_t getnElms() const { return nElms; }; + + std::vector<double> getDeltaStar() const { return deltaStar; }; + std::vector<double> getUe() const + { + std::vector<double> ue(nNodes, 0.); + for (size_t i = 0; i < nNodes; ++i) + ue[i] = u[i * nVar + 3]; + return ue; + }; + std::vector<double> getBlowing() const { return blowingVelocity; }; + double getMaxMach() const { return *std::max_element(Me.begin(), Me.end()); }; + + // Setters + void setMesh(std::vector<double> const _x, std::vector<double> const y, + std::vector<double> const z, double const _chord, + double const _xMin, + std::vector<int> const _rows = std::vector<int>(0), + std::vector<int> const _no = std::vector<int>(0)); + void computeLocalCoordinates(); + void setMe(std::vector<double> const _Me) + { + if (_Me.size() != nNodes) + throw std::runtime_error("Mach number vector is not consistent."); + Me = _Me; + }; + void setVt(std::vector<double> const _vt) + { + if (_vt.size() != nNodes) + throw std::runtime_error("Velocity vector is not consistent."); + vt = _vt; + }; + void setRhoe(std::vector<double> const _rhoe) + { + if (_rhoe.size() != nNodes) + throw std::runtime_error("Density vector is not consistent."); + rhoe = _rhoe; + }; + void setxxExt(std::vector<double> const _xxExt) + { + if (_xxExt.size() != nNodes) + throw std::runtime_error("External mesh vector is not consistent."); + xxExt = _xxExt; + }; + void setDeltaStarExt(std::vector<double> const _deltaStarExt) + { + if (_deltaStarExt.size() != nNodes) + throw std::runtime_error( + "Displacement thickness vector is not consistent."); + deltaStarExt = _deltaStarExt; + }; + void setVtExt(std::vector<double> const _vtExt) + { + if (_vtExt.size() != nNodes) + throw std::runtime_error("Velocity vector is not consistent."); + vtExt = _vtExt; + }; + void setxtrF(double const _xtrF) { xtrF = _xtrF; }; + + // Boundary conditions. + void setStagnationBC(double Re); + void setWakeBC(double Re, std::vector<double> upperConditions, + std::vector<double> lowerConditions); + + // Getters. + size_t getnVar() const { return nVar; }; + double getnCrit() const { return nCrit; }; + void setnCrit(double const _nCrit) { nCrit = _nCrit; }; + + // Others + void printSolution(size_t iPoint) const; + void computeBlowingVelocity(); + std::vector<std::vector<double>> evalGradwrtRho(); + std::vector<std::vector<double>> evalGradwrtVt(); + std::vector<std::vector<double>> evalGradwrtDeltaStar(); + std::vector<std::vector<double>> evalGradwrtX(); + std::vector<std::vector<double>> evalGradwrtY(); }; } // namespace blast #endif // DBOUNDARYLAYER_H diff --git a/blast/src/DClosures.cpp b/blast/src/DClosures.cpp index 1b1fd249193dc183c789854cd0052b2120e9a973..6e608589fa3ce0a3d31899bf6b87eee92818ad70 100644 --- a/blast/src/DClosures.cpp +++ b/blast/src/DClosures.cpp @@ -16,66 +16,72 @@ using namespace blast; */ void Closures::laminarClosures(std::vector<double> &closureVars, double theta, double H, double Ret, const double Me, - const std::string name) { - H = std::max(H, 1.00005); - - // Kinematic shape factor H*. - double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me); - if (name == "wake") - Hk = std::max(Hk, 1.00005); - else - Hk = std::max(Hk, 1.05000); - - // Density shape parameter. - double Hstar2 = (0.064 / (Hk - 0.8) + 0.251) * Me * Me; - - // Boundary layer thickness. - double delta = std::min((3.15 + H + (1.72 / (Hk - 1))) * theta, 12 * theta); - - double Hstar = 0.; - double Hks = Hk - 4.35; - if (Hk <= 4.35) - Hstar = 0.0111 * Hks * Hks / (Hk + 1) - - 0.0278 * Hks * Hks * Hks / (Hk + 1) + 1.528 - - 0.0002 * (Hks * Hk) * (Hks * Hk); - else - Hstar = 1.528 + 0.015 * Hks * Hks / Hk; - - // Whitfield's minor additional compressibility correction. - Hstar = (Hstar + 0.028 * Me * Me) / (1 + 0.014 * Me * Me); - - // Friction coefficient. - double tmp = 0.; - double cf = 0.; - if (Hk < 5.5) { - tmp = (5.5 - Hk) * (5.5 - Hk) * (5.5 - Hk) / (Hk + 1.0); - cf = (-0.07 + 0.0727 * tmp) / Ret; - } else if (Hk >= 5.5) { - tmp = 1.0 - 1.0 / (Hk - 4.5); - cf = (-0.07 + 0.015 * tmp * tmp) / Ret; - } - - // Dissipation coefficient. - double Cd = 0.; - if (Hk < 4) - Cd = (0.00205 * std::pow(4.0 - Hk, 5.5) + 0.207) * (Hstar / (2 * Ret)); - else { - double HkCd = (Hk - 4) * (Hk - 4); - double denCd = 1 + 0.02 * HkCd; - Cd = (-0.0016 * HkCd / denCd + 0.207) * (Hstar / (2 * Ret)); - } - - // Wake relations. - if (name == "wake") { - Cd = 2 * (1.10 * (1.0 - 1.0 / Hk) * (1.0 - 1.0 / Hk) / Hk) * - (Hstar / (2 * Ret)); - cf = 0.; - } - - double us = 0.; - double ctEq = 0.; - - closureVars = {Hstar, Hstar2, Hk, cf, Cd, delta, ctEq, us}; + const std::string name) +{ + H = std::max(H, 1.00005); + + // Kinematic shape factor H*. + double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me); + if (name == "wake") + Hk = std::max(Hk, 1.00005); + else + Hk = std::max(Hk, 1.05000); + + // Density shape parameter. + double Hstar2 = (0.064 / (Hk - 0.8) + 0.251) * Me * Me; + + // Boundary layer thickness. + double delta = std::min((3.15 + H + (1.72 / (Hk - 1))) * theta, 12 * theta); + + double Hstar = 0.; + double Hks = Hk - 4.35; + if (Hk <= 4.35) + Hstar = 0.0111 * Hks * Hks / (Hk + 1) - + 0.0278 * Hks * Hks * Hks / (Hk + 1) + 1.528 - + 0.0002 * (Hks * Hk) * (Hks * Hk); + else + Hstar = 1.528 + 0.015 * Hks * Hks / Hk; + + // Whitfield's minor additional compressibility correction. + Hstar = (Hstar + 0.028 * Me * Me) / (1 + 0.014 * Me * Me); + + // Friction coefficient. + double tmp = 0.; + double cf = 0.; + if (Hk < 5.5) + { + tmp = (5.5 - Hk) * (5.5 - Hk) * (5.5 - Hk) / (Hk + 1.0); + cf = (-0.07 + 0.0727 * tmp) / Ret; + } + else if (Hk >= 5.5) + { + tmp = 1.0 - 1.0 / (Hk - 4.5); + cf = (-0.07 + 0.015 * tmp * tmp) / Ret; + } + + // Dissipation coefficient. + double Cd = 0.; + if (Hk < 4) + Cd = (0.00205 * std::pow(4.0 - Hk, 5.5) + 0.207) * (Hstar / (2 * Ret)); + else + { + double HkCd = (Hk - 4) * (Hk - 4); + double denCd = 1 + 0.02 * HkCd; + Cd = (-0.0016 * HkCd / denCd + 0.207) * (Hstar / (2 * Ret)); + } + + // Wake relations. + if (name == "wake") + { + Cd = 2 * (1.10 * (1.0 - 1.0 / Hk) * (1.0 - 1.0 / Hk) / Hk) * + (Hstar / (2 * Ret)); + cf = 0.; + } + + double us = 0.; + double ctEq = 0.; + + closureVars = {Hstar, Hstar2, Hk, cf, Cd, delta, ctEq, us}; } /** @@ -91,113 +97,117 @@ void Closures::laminarClosures(std::vector<double> &closureVars, double theta, */ void Closures::turbulentClosures(std::vector<double> &closureVars, double theta, double H, double Ct, double Ret, - const double Me, const std::string name) { - H = std::max(H, 1.00005); - Ct = std::min(Ct, 0.3); - // Ct = std::max(std::min(Ct, 0.30), 0.0000001); - - double Hk = (H - 0.29 * Me * Me) / (1. + 0.113 * Me * Me); - if (name == "wake") - Hk = std::max(Hk, 1.00005); - else - Hk = std::max(Hk, 1.05000); - - double Hstar2 = ((0.064 / (Hk - 0.8)) + 0.251) * Me * Me; - double gamma = 1.4 - 1.; - double Fc = std::sqrt(1. + 0.5 * gamma * Me * Me); - - double H0 = 0.; - if (Ret <= 400.) - H0 = 4.0; - else - H0 = 3. + (400. / Ret); - if (Ret <= 200.) - Ret = 200.; - - double Hstar = 0.; - if (Hk <= H0) - Hstar = - ((0.5 - 4. / Ret) * ((H0 - Hk) / (H0 - 1.)) * ((H0 - Hk) / (H0 - 1.))) * - (1.5 / (Hk + 0.5)) + - 1.5 + 4. / Ret; - else - Hstar = (Hk - H0) * (Hk - H0) * - (0.007 * std::log(Ret) / - ((Hk - H0 + 4. / std::log(Ret)) * - (Hk - H0 + 4. / std::log(Ret))) + - 0.015 / Hk) + + const double Me, const std::string name) +{ + H = std::max(H, 1.00005); + Ct = std::min(Ct, 0.3); + // Ct = std::max(std::min(Ct, 0.30), 0.0000001); + + double Hk = (H - 0.29 * Me * Me) / (1. + 0.113 * Me * Me); + if (name == "wake") + Hk = std::max(Hk, 1.00005); + else + Hk = std::max(Hk, 1.05000); + + double Hstar2 = ((0.064 / (Hk - 0.8)) + 0.251) * Me * Me; + double gamma = 1.4 - 1.; + double Fc = std::sqrt(1. + 0.5 * gamma * Me * Me); + + double H0 = 0.; + if (Ret <= 400.) + H0 = 4.0; + else + H0 = 3. + (400. / Ret); + if (Ret <= 200.) + Ret = 200.; + + double Hstar = 0.; + if (Hk <= H0) + Hstar = + ((0.5 - 4. / Ret) * ((H0 - Hk) / (H0 - 1.)) * ((H0 - Hk) / (H0 - 1.))) * + (1.5 / (Hk + 0.5)) + 1.5 + 4. / Ret; + else + Hstar = (Hk - H0) * (Hk - H0) * + (0.007 * std::log(Ret) / + ((Hk - H0 + 4. / std::log(Ret)) * + (Hk - H0 + 4. / std::log(Ret))) + + 0.015 / Hk) + + 1.5 + 4. / Ret; + + // Whitfield's minor additional compressibility correction. + Hstar = (Hstar + 0.028 * Me * Me) / (1. + 0.014 * Me * Me); + + double logRt = std::log(Ret / Fc); + logRt = std::max(logRt, 3.0); + double arg = std::max(-1.33 * Hk, -20.0); + + // Equivalent normalized wall slip velocity. + double us = (Hstar / 2.) * (1. - 4. * (Hk - 1.) / (3. * H)); + + // Boundary layer thickness. + double delta = std::min((3.15 + H + (1.72 / (Hk - 1.))) * theta, 12. * theta); + + double Ctcon = 0.5 / (6.7 * 6.7 * 0.75); + double Hkc = 0.; + double Cdw = 0.; + double Cdd = 0.; + double Cdl = 0.; + + double Hmin = 0.; + double Fl = 0.; + double Dfac = 0.; + + double cf = 0.; + double Cd = 0.; + + double n = 1.0; + + double ctEq = 0.; + + if (name == "wake") + { + if (us > 0.99995) + us = 0.99995; + n = 2.0; // two wake halves + cf = 0.0; // no friction inside the wake + Hkc = Hk - 1.; + Cdw = 0.0; // Wall contribution. + Cdd = (0.995 - us) * Ct * Ct; // Turbulent outer layer contribution. + Cdl = 0.15 * (0.995 - us) * (0.995 - us) / + Ret; // Laminar stress contribution to outer layer. + ctEq = std::sqrt(4. * Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) / + ((1. - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5. + } + else + { + if (us > 0.95) + us = 0.98; + n = 1.0; + cf = (1 / Fc) * + (0.3 * std::exp(arg) * std::pow(logRt / 2.3026, (-1.74 - 0.31 * Hk)) + + 0.00011 * (std::tanh(4. - (Hk / 0.875)) - 1.)); + Hkc = std::max(Hk - 1. - 18. / Ret, 0.01); + // Dissipation coefficient. + Hmin = 1. + 2.1 / std::log(Ret); + Fl = (Hk - 1.) / (Hmin - 1); + Dfac = 0.5 + 0.5 * std::tanh(Fl); + Cdw = 0.5 * (cf * us) * Dfac; + Cdd = (0.995 - us) * Ct * Ct; + Cdl = 0.15 * (0.995 - us) * (0.995 - us) / Ret; + ctEq = std::sqrt(Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) / + ((1. - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5 + // ctEq = std::sqrt(Hstar * 0.015/(1-us) * (Hk-1)*(Hk-1)*(Hk-1)/(Hk*Hk*H)); + // // Drela 1987 + } + if (n * Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) / + ((1. - us) * (Hk * Hk) * H) < + 0.) + std::cout << "Negative sqrt encountered " << std::endl; - // Whitfield's minor additional compressibility correction. - Hstar = (Hstar + 0.028 * Me * Me) / (1. + 0.014 * Me * Me); - - double logRt = std::log(Ret / Fc); - logRt = std::max(logRt, 3.0); - double arg = std::max(-1.33 * Hk, -20.0); - - // Equivalent normalized wall slip velocity. - double us = (Hstar / 2.) * (1. - 4. * (Hk - 1.) / (3. * H)); - - // Boundary layer thickness. - double delta = std::min((3.15 + H + (1.72 / (Hk - 1.))) * theta, 12. * theta); - - double Ctcon = 0.5 / (6.7 * 6.7 * 0.75); - double Hkc = 0.; - double Cdw = 0.; - double Cdd = 0.; - double Cdl = 0.; - - double Hmin = 0.; - double Fl = 0.; - double Dfac = 0.; - - double cf = 0.; - double Cd = 0.; - - double n = 1.0; - - double ctEq = 0.; - - if (name == "wake") { - if (us > 0.99995) - us = 0.99995; - n = 2.0; // two wake halves - cf = 0.0; // no friction inside the wake - Hkc = Hk - 1.; - Cdw = 0.0; // Wall contribution. - Cdd = (0.995 - us) * Ct * Ct; // Turbulent outer layer contribution. - Cdl = 0.15 * (0.995 - us) * (0.995 - us) / - Ret; // Laminar stress contribution to outer layer. - ctEq = std::sqrt(4. * Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) / - ((1. - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5. - } else { - if (us > 0.95) - us = 0.98; - n = 1.0; - cf = (1 / Fc) * - (0.3 * std::exp(arg) * std::pow(logRt / 2.3026, (-1.74 - 0.31 * Hk)) + - 0.00011 * (std::tanh(4. - (Hk / 0.875)) - 1.)); - Hkc = std::max(Hk - 1. - 18. / Ret, 0.01); // Dissipation coefficient. - Hmin = 1. + 2.1 / std::log(Ret); - Fl = (Hk - 1.) / (Hmin - 1); - Dfac = 0.5 + 0.5 * std::tanh(Fl); - Cdw = 0.5 * (cf * us) * Dfac; - Cdd = (0.995 - us) * Ct * Ct; - Cdl = 0.15 * (0.995 - us) * (0.995 - us) / Ret; - ctEq = std::sqrt(Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) / - ((1. - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5 - // ctEq = std::sqrt(Hstar * 0.015/(1-us) * (Hk-1)*(Hk-1)*(Hk-1)/(Hk*Hk*H)); - // // Drela 1987 - } - if (n * Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) / - ((1. - us) * (Hk * Hk) * H) < - 0.) - std::cout << "Negative sqrt encountered " << std::endl; - - // Dissipation coefficient. - Cd = n * (Cdw + Cdd + Cdl); - closureVars = {Hstar, Hstar2, Hk, cf, Cd, delta, ctEq, us}; + Cd = n * (Cdw + Cdd + Cdl); + closureVars = {Hstar, Hstar2, Hk, cf, Cd, delta, ctEq, us}; } /** @@ -213,65 +223,69 @@ void Closures::turbulentClosures(std::vector<double> &closureVars, double theta, */ void Closures::turbulentClosures(double &closureVars, double theta, double H, double Ct, double Ret, const double Me, - const std::string name) { - H = std::max(H, 1.00005); - Ct = std::min(Ct, 0.3); - - double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me); - if (name == "wake") - Hk = std::max(Hk, 1.00005); - else - Hk = std::max(Hk, 1.05000); - - double H0 = 0.; - if (Ret <= 400.) - H0 = 4.0; - else - H0 = 3. + (400. / Ret); - if (Ret <= 200.) - Ret = 200.; - - double Hstar = 0.; - if (Hk <= H0) - Hstar = - ((0.5 - 4. / Ret) * ((H0 - Hk) / (H0 - 1.)) * ((H0 - Hk) / (H0 - 1.))) * - (1.5 / (Hk + 0.5)) + - 1.5 + 4. / Ret; - else - Hstar = (Hk - H0) * (Hk - H0) * - (0.007 * std::log(Ret) / - ((Hk - H0 + 4. / std::log(Ret)) * - (Hk - H0 + 4. / std::log(Ret))) + - 0.015 / Hk) + + const std::string name) +{ + H = std::max(H, 1.00005); + Ct = std::min(Ct, 0.3); + + double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me); + if (name == "wake") + Hk = std::max(Hk, 1.00005); + else + Hk = std::max(Hk, 1.05000); + + double H0 = 0.; + if (Ret <= 400.) + H0 = 4.0; + else + H0 = 3. + (400. / Ret); + if (Ret <= 200.) + Ret = 200.; + + double Hstar = 0.; + if (Hk <= H0) + Hstar = + ((0.5 - 4. / Ret) * ((H0 - Hk) / (H0 - 1.)) * ((H0 - Hk) / (H0 - 1.))) * + (1.5 / (Hk + 0.5)) + 1.5 + 4. / Ret; - - // Whitfield's minor additional compressibility correction. - Hstar = (Hstar + 0.028 * Me * Me) / (1. + 0.014 * Me * Me); - - // Equivalent normalized wall slip velocity. - double us = (Hstar / 2.) * (1 - 4 * (Hk - 1.) / (3. * H)); - - double Ctcon = 0.5 / (6.7 * 6.7 * 0.75); - - double Hkc = 0.; - double ctEq = 0.; - - if (name == "wake") { - if (us > 0.99995) - us = 0.99995; - Hkc = Hk - 1.; - ctEq = std::sqrt(4. * Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) / - ((1. - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5 - } else { - if (us > 0.95) - us = 0.98; - Hkc = std::max(Hk - 1. - 18. / Ret, 0.01); - ctEq = std::sqrt(Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) / - ((1. - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5 - } - if (Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) / ((1. - us) * (Hk * Hk) * H) < - 0.) - std::cout << "Negative sqrt encountered " << std::endl; - - closureVars = ctEq; + else + Hstar = (Hk - H0) * (Hk - H0) * + (0.007 * std::log(Ret) / + ((Hk - H0 + 4. / std::log(Ret)) * + (Hk - H0 + 4. / std::log(Ret))) + + 0.015 / Hk) + + 1.5 + 4. / Ret; + + // Whitfield's minor additional compressibility correction. + Hstar = (Hstar + 0.028 * Me * Me) / (1. + 0.014 * Me * Me); + + // Equivalent normalized wall slip velocity. + double us = (Hstar / 2.) * (1 - 4 * (Hk - 1.) / (3. * H)); + + double Ctcon = 0.5 / (6.7 * 6.7 * 0.75); + + double Hkc = 0.; + double ctEq = 0.; + + if (name == "wake") + { + if (us > 0.99995) + us = 0.99995; + Hkc = Hk - 1.; + ctEq = std::sqrt(4. * Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) / + ((1. - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5 + } + else + { + if (us > 0.95) + us = 0.98; + Hkc = std::max(Hk - 1. - 18. / Ret, 0.01); + ctEq = std::sqrt(Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) / + ((1. - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5 + } + if (Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) / ((1. - us) * (Hk * Hk) * H) < + 0.) + std::cout << "Negative sqrt encountered " << std::endl; + + closureVars = ctEq; } \ No newline at end of file diff --git a/blast/src/DClosures.h b/blast/src/DClosures.h index 1f3841f326f3cd10c07c575f70ab98bff8922b74..f5ffe557dd4b3eae504646fd58b355f6fe2eefa9 100644 --- a/blast/src/DClosures.h +++ b/blast/src/DClosures.h @@ -5,24 +5,26 @@ #include <string> #include <vector> -namespace blast { +namespace blast +{ /** * @brief Boundary layer closure relations. * @author Paul Dechamps */ -class BLAST_API Closures { +class BLAST_API Closures +{ public: - static void laminarClosures(std::vector<double> &closureVars, double theta, - double H, double Ret, const double Me, - const std::string name); - static void turbulentClosures(std::vector<double> &closureVars, double theta, - double H, double Ct, double Ret, double Me, - const std::string name); - static void turbulentClosures(double &closureVars, double theta, double H, - double Ct, double Ret, const double Me, + static void laminarClosures(std::vector<double> &closureVars, double theta, + double H, double Ret, const double Me, const std::string name); + static void turbulentClosures(std::vector<double> &closureVars, double theta, + double H, double Ct, double Ret, double Me, + const std::string name); + static void turbulentClosures(double &closureVars, double theta, double H, + double Ct, double Ret, const double Me, + const std::string name); }; } // namespace blast #endif // DCLOSURES_H diff --git a/blast/src/DCoupledAdjoint.cpp b/blast/src/DCoupledAdjoint.cpp index 1911ebe15d2c15424711e1657e08625c10ed3c7c..bdf2f5715fd020a195323f50e3c439d61cf7f004 100644 --- a/blast/src/DCoupledAdjoint.cpp +++ b/blast/src/DCoupledAdjoint.cpp @@ -45,408 +45,427 @@ using namespace blast; CoupledAdjoint::CoupledAdjoint(std::shared_ptr<dart::Adjoint> _iadjoint, std::shared_ptr<blast::Driver> vSolver) - : iadjoint(_iadjoint), vsol(vSolver) { - isol = iadjoint->sol; - // Linear solvers (if GMRES, use the same, but with a thighter tolerance) - if (isol->linsol->type() == tbox::LSolType::GMRES_ILUT) { - std::shared_ptr<tbox::Gmres> gmres = std::make_shared<tbox::Gmres>( - *dynamic_cast<tbox::Gmres *>(isol->linsol.get())); - gmres->setTolerance(1e-8); - alinsol = gmres; - } else - alinsol = isol->linsol; - - // Initalize all derivatives - this->reset(); + : iadjoint(_iadjoint), vsol(vSolver) +{ + isol = iadjoint->sol; + // Linear solvers (if GMRES, use the same, but with a thighter tolerance) + if (isol->linsol->type() == tbox::LSolType::GMRES_ILUT) + { + std::shared_ptr<tbox::Gmres> gmres = std::make_shared<tbox::Gmres>( + *dynamic_cast<tbox::Gmres *>(isol->linsol.get())); + gmres->setTolerance(1e-8); + alinsol = gmres; + } + else + alinsol = isol->linsol; + + // Initalize all derivatives + this->reset(); } -void CoupledAdjoint::reset() { - size_t nDim = isol->pbl->nDim; // Number of dimensions of the problem - size_t nNd = isol->pbl->msh->nodes.size(); // Number of domain nodes - size_t nNs = 0; - for (auto bBC : isol->pbl->bBCs) - nNs += bBC->nodes.size(); // Number of surface nodes - nNs += 1; // Duplicate stagnation point - size_t nEs = 0; - for (auto bBC : isol->pbl->bBCs) - nEs += bBC->uE.size(); // Number of blowing elements - - dRphi_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nNd); - dRM_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNd); - dRrho_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNd); - dRv_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim * nNs, nNd); - dRdStar_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNd); - dRblw_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nNd); - - dRphi_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nNs); - dRM_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs); - dRrho_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs); - dRv_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim * nNs, nNs); - dRdStar_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs); - dRblw_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nNs); - - dRphi_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nNs); - dRM_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs); - dRrho_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs); - dRv_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim * nNs, nNs); - dRdStar_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs); - dRblw_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nNs); - - dRphi_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nDim * nNs); - dRM_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim * nNs); - dRrho_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim * nNs); - dRv_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim * nNs, nDim * nNs); - dRdStar_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim * nNs); - dRblw_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nDim * nNs); - - dRphi_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nNs); - dRM_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs); - dRrho_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs); - dRv_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim * nNs, nNs); - dRdStar_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs); - dRblw_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nNs); - - dRphi_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nEs); - dRM_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nEs); - dRrho_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nEs); - dRv_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim * nNs, nEs); - dRdStar_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nEs); - dRblw_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nEs); - - dRphi_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nDim * nNd); - dRM_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim * nNd); - dRrho_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim * nNd); - dRv_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim * nNs, nDim * nNd); - dRdStar_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim * nNd); - dRblw_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nDim * nNd); - - // Gradient wrt objective function - dRphi_AoA = Eigen::VectorXd::Zero(nNd); - - // Objective functions - // Cl = f(phi) & Cd = Cdp + Cdf = f(phi, M, v) - dCl_phi = Eigen::VectorXd::Zero(nNd); - dCd_phi = Eigen::VectorXd::Zero(nNd); - dCd_M = Eigen::VectorXd::Zero(nNs); - dCd_v = Eigen::VectorXd::Zero(nDim * nNs); - dCd_dStar = Eigen::VectorXd::Zero(nNs); - - // Angle of attack - dCl_AoA = 0.0; - dCd_AoA = 0.0; - tdCl_AoA = 0.0; - tdCd_AoA = 0.0; - - // Mesh - dRx_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim * nNd, nDim * nNd); - dCl_x = Eigen::VectorXd::Zero(nDim * nNd); - dCdp_x = Eigen::VectorXd::Zero(nDim * nNd); - dCdf_x = Eigen::VectorXd::Zero(nDim * nNd); - - tdCl_x.resize(nNd, Eigen::Vector3d::Zero()); - tdCd_x.resize(nNd, Eigen::Vector3d::Zero()); +void CoupledAdjoint::reset() +{ + size_t nDim = isol->pbl->nDim; // Number of dimensions of the problem + size_t nNd = isol->pbl->msh->nodes.size(); // Number of domain nodes + size_t nNs = 0; + for (auto bBC : isol->pbl->bBCs) + nNs += bBC->nodes.size(); // Number of surface nodes + nNs += 1; // Duplicate stagnation point + size_t nEs = 0; + for (auto bBC : isol->pbl->bBCs) + nEs += bBC->uE.size(); // Number of blowing elements + + dRphi_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nNd); + dRM_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNd); + dRrho_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNd); + dRv_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim * nNs, nNd); + dRdStar_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNd); + dRblw_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nNd); + + dRphi_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nNs); + dRM_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs); + dRrho_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs); + dRv_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim * nNs, nNs); + dRdStar_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs); + dRblw_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nNs); + + dRphi_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nNs); + dRM_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs); + dRrho_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs); + dRv_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim * nNs, nNs); + dRdStar_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs); + dRblw_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nNs); + + dRphi_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nDim * nNs); + dRM_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim * nNs); + dRrho_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim * nNs); + dRv_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim * nNs, nDim * nNs); + dRdStar_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim * nNs); + dRblw_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nDim * nNs); + + dRphi_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nNs); + dRM_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs); + dRrho_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs); + dRv_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim * nNs, nNs); + dRdStar_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs); + dRblw_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nNs); + + dRphi_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nEs); + dRM_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nEs); + dRrho_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nEs); + dRv_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim * nNs, nEs); + dRdStar_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nEs); + dRblw_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nEs); + + dRphi_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nDim * nNd); + dRM_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim * nNd); + dRrho_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim * nNd); + dRv_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim * nNs, nDim * nNd); + dRdStar_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim * nNd); + dRblw_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nDim * nNd); + + // Gradient wrt objective function + dRphi_AoA = Eigen::VectorXd::Zero(nNd); + + // Objective functions + // Cl = f(phi) & Cd = Cdp + Cdf = f(phi, M, v) + dCl_phi = Eigen::VectorXd::Zero(nNd); + dCd_phi = Eigen::VectorXd::Zero(nNd); + dCd_M = Eigen::VectorXd::Zero(nNs); + dCd_v = Eigen::VectorXd::Zero(nDim * nNs); + dCd_dStar = Eigen::VectorXd::Zero(nNs); + + // Angle of attack + dCl_AoA = 0.0; + dCd_AoA = 0.0; + tdCl_AoA = 0.0; + tdCd_AoA = 0.0; + + // Mesh + dRx_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim * nNd, nDim * nNd); + dCl_x = Eigen::VectorXd::Zero(nDim * nNd); + dCdp_x = Eigen::VectorXd::Zero(nDim * nNd); + dCdf_x = Eigen::VectorXd::Zero(nDim * nNd); + + tdCl_x.resize(nNd, Eigen::Vector3d::Zero()); + tdCd_x.resize(nNd, Eigen::Vector3d::Zero()); } void CoupledAdjoint::buildAdjointMatrix( std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor> *>> &matrices, - Eigen::SparseMatrix<double, Eigen::RowMajor> &matrixAdjoint) { - // Create a list of triplets for the larger matrix - std::deque<Eigen::Triplet<double>> triplets; - - // Sanity check - for (size_t irow = 0; irow < matrices.size(); ++irow) { - if (irow > 0) - if (matrices[irow].size() != matrices[irow - 1].size()) { - std::cout << "WRONG ROW SIZE " << std::endl; - throw; - } - for (size_t icol = 0; icol < matrices[irow].size(); ++icol) - if (icol > 0) - if (matrices[irow][icol]->rows() != matrices[irow][icol - 1]->rows()) { - std::cout << "WRONG COL SIZE " << std::endl; - throw; - } - } - - // Iterate over the rows and columns of the vector - int rowOffset = 0; - for (const auto &row : matrices) { - int colOffset = 0; - for (const auto &matrix : row) { - // Add the triplets of the matrix to the list of triplets for the larger - // matrix - for (int k = 0; k < matrix->outerSize(); ++k) { - for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it( - *matrix, k); - it; ++it) { - triplets.push_back(Eigen::Triplet<double>( - it.row() + rowOffset, it.col() + colOffset, it.value())); + Eigen::SparseMatrix<double, Eigen::RowMajor> &matrixAdjoint) +{ + // Create a list of triplets for the larger matrix + std::deque<Eigen::Triplet<double>> triplets; + + // Sanity check + for (size_t irow = 0; irow < matrices.size(); ++irow) + { + if (irow > 0) + if (matrices[irow].size() != matrices[irow - 1].size()) + { + std::cout << "WRONG ROW SIZE " << std::endl; + throw; + } + for (size_t icol = 0; icol < matrices[irow].size(); ++icol) + if (icol > 0) + if (matrices[irow][icol]->rows() != matrices[irow][icol - 1]->rows()) + { + std::cout << "WRONG COL SIZE " << std::endl; + throw; + } + } + + // Iterate over the rows and columns of the vector + int rowOffset = 0; + for (const auto &row : matrices) + { + int colOffset = 0; + for (const auto &matrix : row) + { + // Add the triplets of the matrix to the list of triplets for the larger + // matrix + for (int k = 0; k < matrix->outerSize(); ++k) + { + for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it( + *matrix, k); + it; ++it) + { + triplets.push_back(Eigen::Triplet<double>( + it.row() + rowOffset, it.col() + colOffset, it.value())); + } + } + colOffset += matrix->cols(); } - } - colOffset += matrix->cols(); + rowOffset += row[0]->rows(); } - rowOffset += row[0]->rows(); - } - // Create a new matrix from the deque of triplets - matrixAdjoint.setFromTriplets(triplets.begin(), triplets.end()); - matrixAdjoint.prune(0.); - matrixAdjoint.makeCompressed(); + // Create a new matrix from the deque of triplets + matrixAdjoint.setFromTriplets(triplets.begin(), triplets.end()); + matrixAdjoint.prune(0.); + matrixAdjoint.makeCompressed(); } -void CoupledAdjoint::run() { - tms["0-Total"].start(); - tms["1-Derivatives"].start(); - tms2["1-adj_inviscid"].start(); - gradientswrtInviscidFlow(); - tms2["1-adj_inviscid"].stop(); - tms2["2-adj_mach"].start(); - gradientswrtMachNumber(); - tms2["2-adj_mach"].stop(); - tms2["3-adj_density"].start(); - gradientswrtDensity(); - tms2["3-adj_density"].stop(); - tms2["4-adj_velocity"].start(); - gradientswrtVelocity(); - tms2["4-adj_velocity"].stop(); - tms2["5-adj_dStar"].start(); - gradientswrtDeltaStar(); - tms2["5-adj_dStar"].stop(); - tms2["6-adj_blw"].start(); - gradientswrtBlowingVelocity(); - tms2["6-adj_blw"].stop(); - tms2["7-adj_aoa"].start(); - gradientwrtAoA(); - tms2["7-adj_aoa"].stop(); - tms2["8-adj_mesh"].start(); - gradientwrtMesh(); - tms2["8-adj_mesh"].stop(); - tms["1-Derivatives"].stop(); - - transposeMatrices(); - - tms2["9-build"].start(); - std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor> *>> - matrices = {{&dRphi_phi, &dRM_phi, &dRrho_phi, &dRv_phi, &dRdStar_phi, - &dRblw_phi}, - {&dRphi_M, &dRM_M, &dRrho_M, &dRv_M, &dRdStar_M, &dRblw_M}, - {&dRphi_rho, &dRM_rho, &dRrho_rho, &dRv_rho, &dRdStar_rho, - &dRblw_rho}, - {&dRphi_v, &dRM_v, &dRrho_v, &dRv_v, &dRdStar_v, &dRblw_v}, - {&dRphi_dStar, &dRM_dStar, &dRrho_dStar, &dRv_dStar, - &dRdStar_dStar, &dRblw_dStar}, - {&dRphi_blw, &dRM_blw, &dRrho_blw, &dRv_blw, &dRdStar_blw, - &dRblw_blw}}; - - int rows = 0; - int cols = 0; - for (const auto &row : matrices) - rows += - row[0]->rows(); // All matrices in a row have the same number of rows - for (const auto &mat1 : matrices[0]) - cols += mat1->cols(); // All matrices in a column have the same number of - // columns - Eigen::SparseMatrix<double, Eigen::RowMajor> A_adjoint(rows, cols); - buildAdjointMatrix(matrices, A_adjoint); - tms2["9-build"].stop(); - - size_t phiIdx = 0; - size_t machIdx = dRphi_phi.cols(); - size_t rhoIdx = machIdx + dRM_M.cols(); - size_t vIdx = rhoIdx + dRrho_rho.cols(); - size_t dStarIdx = vIdx + dRv_v.cols(); - size_t blwIdx = dStarIdx + dRdStar_dStar.cols(); - - //**************************************************************************// - // Gradients wrt angle of attack // - //**************************************************************************// - // CL - tms["2-Solve"].start(); - tms2["10-solveCxAoA"].start(); - std::vector<double> rhsCl(A_adjoint.cols(), 0.0); - for (auto i = phiIdx; i < phiIdx + dCl_phi.size(); ++i) - rhsCl[i] = dCl_phi[i]; // Only dCl_dphi is non-zero - Eigen::VectorXd rhsCl_ = - Eigen::Map<const Eigen::VectorXd>(rhsCl.data(), rhsCl.size()); - - Eigen::VectorXd lambdaCl(A_adjoint.cols()); // Solution vector for lambdasCl - Eigen::Map<Eigen::VectorXd> lambdaCl_(lambdaCl.data(), lambdaCl.size()); - - alinsol->compute(A_adjoint, - Eigen::Map<Eigen::VectorXd>(rhsCl_.data(), rhsCl_.size()), - lambdaCl_); - - // Total gradient - tdCl_AoA = dCl_AoA - - dRphi_AoA.transpose() * lambdaCl.segment(phiIdx, dRphi_phi.cols()); - - // CD - std::vector<double> rhsCd(A_adjoint.cols(), 0.0); - int jj = 0; - for (auto i = phiIdx; i < phiIdx + dCd_phi.size(); ++i) { - rhsCd[i] = dCd_phi[jj]; - jj++; - } - jj = 0; - for (auto i = machIdx; i < machIdx + dCd_M.size(); ++i) { - rhsCd[i] = dCd_M[jj]; - jj++; - } - jj = 0; - for (auto i = vIdx; i < vIdx + dCd_v.size(); ++i) { - rhsCd[i] = dCd_v[jj]; - jj++; - } - jj = 0; - for (auto i = dStarIdx; i < dStarIdx + dCd_dStar.size(); ++i) { - rhsCd[i] = dCd_dStar[jj]; - jj++; - } - Eigen::VectorXd rhsCd_ = - Eigen::Map<const Eigen::VectorXd>(rhsCd.data(), rhsCd.size()); - - Eigen::VectorXd lambdaCd(A_adjoint.cols()); // Solution vector for lambdasCd - Eigen::Map<Eigen::VectorXd> lambdaCd_(lambdaCd.data(), lambdaCd.size()); - - alinsol->compute(A_adjoint, - Eigen::Map<Eigen::VectorXd>(rhsCd_.data(), rhsCd_.size()), - lambdaCd_); - - // Total gradient - tdCd_AoA = dCd_AoA - - dRphi_AoA.transpose() * - lambdaCd.segment(phiIdx, dRphi_phi.cols()); // Total gradient - tms2["10-solveCxAoA"].stop(); - - //**************************************************************************// - // Gradients wrt mesh coordinates // - //--------------------------------------------------------------------------// - // Ck (cL, cd) // - // Solution lambdaCl_x of (from augmented Lagrangian, eqn d/dx ) // - // "dCk_x + dRphi_x * lambdaCk_phi + dRdStar_x * lambdaCk_x // - // + dRblw_x * lambdaCk_blw + dRx_x * lambdaCk_x = 0" // - //**************************************************************************// - - tms2["12-solveCxMesh"].start(); - Eigen::VectorXd rhsCl_x = dCl_x; - rhsCl_x -= - dRphi_x.transpose() * - lambdaCl.segment(phiIdx, dRphi_phi.cols()); // Potential contribution - rhsCl_x -= - dRM_x.transpose() * - lambdaCl.segment(machIdx, dRM_M.cols()); // Mach number contribution - rhsCl_x -= dRrho_x.transpose() * - lambdaCl.segment(rhoIdx, dRrho_rho.cols()); // Density contribution - rhsCl_x -= dRv_x.transpose() * - lambdaCl.segment(vIdx, dRv_v.cols()); // Velocity contribution - rhsCl_x -= dRdStar_x.transpose() * - lambdaCl.segment( - dStarIdx, - dRdStar_dStar.cols()); // Displacement thickness contribution - rhsCl_x -= dRblw_x.transpose() * - lambdaCl.segment(blwIdx, dRblw_blw.cols()); // Blowing contribution - - Eigen::VectorXd lambdaCl_x = - Eigen::VectorXd::Zero(isol->pbl->nDim * isol->pbl->msh->nodes.size()); - Eigen::Map<Eigen::VectorXd> lambdaCl_x_(lambdaCl_x.data(), lambdaCl_x.size()); - alinsol->compute(dRx_x.transpose(), - Eigen::Map<Eigen::VectorXd>(rhsCl_x.data(), rhsCl_x.size()), - lambdaCl_x_); - - Eigen::VectorXd rhsCd_x = dCdp_x; // Pressure drag coefficient contribution - rhsCd_x += dCdf_x; // Friction drag coefficient contribution - rhsCd_x -= - dRphi_x.transpose() * - lambdaCd.segment(phiIdx, dRphi_phi.cols()); // Potential contribution - rhsCd_x -= - dRM_x.transpose() * - lambdaCd.segment(machIdx, dRM_M.cols()); // Mach number contribution - rhsCd_x -= dRrho_x.transpose() * - lambdaCd.segment(rhoIdx, dRrho_rho.cols()); // Density contribution - rhsCd_x -= dRv_x.transpose() * - lambdaCd.segment(vIdx, dRv_v.cols()); // Velocity contribution - rhsCd_x -= dRdStar_x.transpose() * - lambdaCd.segment( - dStarIdx, - dRdStar_dStar.cols()); // Displacement thickness contribution - rhsCd_x -= dRblw_x.transpose() * - lambdaCd.segment(blwIdx, dRblw_blw.cols()); // Blowing contribution - - Eigen::VectorXd lambdaCd_x = - Eigen::VectorXd::Zero(isol->pbl->nDim * isol->pbl->msh->nodes.size()); - Eigen::Map<Eigen::VectorXd> lambdaCd_x_(lambdaCd_x.data(), lambdaCd_x.size()); - alinsol->compute(dRx_x.transpose(), - Eigen::Map<Eigen::VectorXd>(rhsCd_x.data(), rhsCd_x.size()), - lambdaCd_x_); - - for (auto n : isol->pbl->msh->nodes) { - for (int m = 0; m < isol->pbl->nDim; m++) { - tdCl_x[n->row](m) = lambdaCl_x[isol->pbl->nDim * (n->row) + m]; - tdCd_x[n->row](m) = lambdaCd_x[isol->pbl->nDim * (n->row) + m]; +void CoupledAdjoint::run() +{ + tms["0-Total"].start(); + tms["1-Derivatives"].start(); + tms2["1-adj_inviscid"].start(); + gradientswrtInviscidFlow(); + tms2["1-adj_inviscid"].stop(); + tms2["2-adj_mach"].start(); + gradientswrtMachNumber(); + tms2["2-adj_mach"].stop(); + tms2["3-adj_density"].start(); + gradientswrtDensity(); + tms2["3-adj_density"].stop(); + tms2["4-adj_velocity"].start(); + gradientswrtVelocity(); + tms2["4-adj_velocity"].stop(); + tms2["5-adj_dStar"].start(); + gradientswrtDeltaStar(); + tms2["5-adj_dStar"].stop(); + tms2["6-adj_blw"].start(); + gradientswrtBlowingVelocity(); + tms2["6-adj_blw"].stop(); + tms2["7-adj_aoa"].start(); + gradientwrtAoA(); + tms2["7-adj_aoa"].stop(); + tms2["8-adj_mesh"].start(); + gradientwrtMesh(); + tms2["8-adj_mesh"].stop(); + tms["1-Derivatives"].stop(); + + transposeMatrices(); + + tms2["9-build"].start(); + std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor> *>> + matrices = {{&dRphi_phi, &dRM_phi, &dRrho_phi, &dRv_phi, &dRdStar_phi, + &dRblw_phi}, + {&dRphi_M, &dRM_M, &dRrho_M, &dRv_M, &dRdStar_M, &dRblw_M}, + {&dRphi_rho, &dRM_rho, &dRrho_rho, &dRv_rho, &dRdStar_rho, + &dRblw_rho}, + {&dRphi_v, &dRM_v, &dRrho_v, &dRv_v, &dRdStar_v, &dRblw_v}, + {&dRphi_dStar, &dRM_dStar, &dRrho_dStar, &dRv_dStar, + &dRdStar_dStar, &dRblw_dStar}, + {&dRphi_blw, &dRM_blw, &dRrho_blw, &dRv_blw, &dRdStar_blw, + &dRblw_blw}}; + + int rows = 0; + int cols = 0; + for (const auto &row : matrices) + rows += + row[0]->rows(); // All matrices in a row have the same number of rows + for (const auto &mat1 : matrices[0]) + cols += mat1->cols(); // All matrices in a column have the same number of + // columns + Eigen::SparseMatrix<double, Eigen::RowMajor> A_adjoint(rows, cols); + buildAdjointMatrix(matrices, A_adjoint); + tms2["9-build"].stop(); + + size_t phiIdx = 0; + size_t machIdx = dRphi_phi.cols(); + size_t rhoIdx = machIdx + dRM_M.cols(); + size_t vIdx = rhoIdx + dRrho_rho.cols(); + size_t dStarIdx = vIdx + dRv_v.cols(); + size_t blwIdx = dStarIdx + dRdStar_dStar.cols(); + + //**************************************************************************// + // Gradients wrt angle of attack // + //**************************************************************************// + // CL + tms["2-Solve"].start(); + tms2["10-solveCxAoA"].start(); + std::vector<double> rhsCl(A_adjoint.cols(), 0.0); + for (auto i = phiIdx; i < phiIdx + dCl_phi.size(); ++i) + rhsCl[i] = dCl_phi[i]; // Only dCl_dphi is non-zero + Eigen::VectorXd rhsCl_ = + Eigen::Map<const Eigen::VectorXd>(rhsCl.data(), rhsCl.size()); + + Eigen::VectorXd lambdaCl(A_adjoint.cols()); // Solution vector for lambdasCl + Eigen::Map<Eigen::VectorXd> lambdaCl_(lambdaCl.data(), lambdaCl.size()); + + alinsol->compute(A_adjoint, + Eigen::Map<Eigen::VectorXd>(rhsCl_.data(), rhsCl_.size()), + lambdaCl_); + + // Total gradient + tdCl_AoA = dCl_AoA - + dRphi_AoA.transpose() * lambdaCl.segment(phiIdx, dRphi_phi.cols()); + + // CD + std::vector<double> rhsCd(A_adjoint.cols(), 0.0); + int jj = 0; + for (auto i = phiIdx; i < phiIdx + dCd_phi.size(); ++i) + { + rhsCd[i] = dCd_phi[jj]; + jj++; + } + jj = 0; + for (auto i = machIdx; i < machIdx + dCd_M.size(); ++i) + { + rhsCd[i] = dCd_M[jj]; + jj++; + } + jj = 0; + for (auto i = vIdx; i < vIdx + dCd_v.size(); ++i) + { + rhsCd[i] = dCd_v[jj]; + jj++; + } + jj = 0; + for (auto i = dStarIdx; i < dStarIdx + dCd_dStar.size(); ++i) + { + rhsCd[i] = dCd_dStar[jj]; + jj++; } - } - tms2["12-solveCxMesh"].stop(); - tms["2-Solve"].stop(); - - // Print output - std::cout << "-------------- Adjoint solution --------------" << std::endl; - // aoa gradients - std::cout << std::setw(15) << std::left << "Adjoint AoA" << std::setw(15) - << std::right << "Res[Cl_AoA]" << std::setw(15) << std::right - << "Res[Cd_AoA]" << std::endl; - std::cout << std::fixed << std::setprecision(5); - std::cout << std::setw(15) << std::left << "" << std::setw(15) << std::right - << log10((A_adjoint * Eigen::Map<Eigen::VectorXd>(lambdaCl.data(), - lambdaCl.size()) - - Eigen::Map<Eigen::VectorXd>(rhsCl.data(), rhsCl.size())) - .norm()) - << std::setw(15) << std::right - << log10((A_adjoint * Eigen::Map<Eigen::VectorXd>(lambdaCd.data(), - lambdaCd.size()) - - Eigen::Map<Eigen::VectorXd>(rhsCd.data(), rhsCd.size())) - .norm()) - << std::endl; - std::cout << std::setw(15) << std::left << "AoA gradients" << std::setw(15) - << std::right << "dCl_dAoA" << std::setw(15) << std::right - << "dCd_dAoA" << std::endl; - std::cout << std::fixed << std::setprecision(5); - std::cout << std::setw(15) << std::left << "" << std::setw(15) << std::right - << tdCl_AoA << std::setw(15) << std::right << tdCd_AoA << std::endl; - // mesh gradients - std::cout << std::setw(15) << std::left << "Adjoint mesh" << std::setw(15) - << std::right << "Res[Cl_x]" << std::setw(15) << std::right - << "Res[Cd_x]" << std::endl; - std::cout << std::fixed << std::setprecision(5); - std::cout - << std::setw(15) << std::left << "" << std::setw(15) << std::right - << log10((dRx_x.transpose() * Eigen::Map<Eigen::VectorXd>( - lambdaCl_x.data(), lambdaCl_x.size()) - - Eigen::Map<Eigen::VectorXd>(rhsCl_x.data(), rhsCl_x.size())) - .norm()) - << std::setw(15) << std::right - << log10((dRx_x.transpose() * Eigen::Map<Eigen::VectorXd>( - lambdaCd_x.data(), lambdaCd_x.size()) - - Eigen::Map<Eigen::VectorXd>(rhsCd_x.data(), rhsCd_x.size())) - .norm()) - << std::endl; - std::cout << std::setw(15) << std::left << "Mesh gradients" << std::setw(15) - << std::right << "Norm[dCl_dX]" << std::setw(15) << std::right - << "Norm[dCd_dX]" << std::endl; - std::cout << std::fixed << std::setprecision(5); - std::cout << std::setw(15) << std::left << "" << std::setw(15) << std::right - << Eigen::Map<Eigen::VectorXd>(lambdaCl_x.data(), lambdaCl_x.size()) - .norm() - << std::setw(15) << std::right - << Eigen::Map<Eigen::VectorXd>(lambdaCd_x.data(), lambdaCd_x.size()) - .norm() - << std::endl; - tms["0-Total"].stop(); - std::cout << "--------------- Adjoint timers ---------------" << std::endl; - std::cout << tms << "----------------------------------------------" - << std::endl; - if (vsol->verbose > 2) - std::cout << tms2 << "----------------------------------------------" + Eigen::VectorXd rhsCd_ = + Eigen::Map<const Eigen::VectorXd>(rhsCd.data(), rhsCd.size()); + + Eigen::VectorXd lambdaCd(A_adjoint.cols()); // Solution vector for lambdasCd + Eigen::Map<Eigen::VectorXd> lambdaCd_(lambdaCd.data(), lambdaCd.size()); + + alinsol->compute(A_adjoint, + Eigen::Map<Eigen::VectorXd>(rhsCd_.data(), rhsCd_.size()), + lambdaCd_); + + // Total gradient + tdCd_AoA = dCd_AoA - + dRphi_AoA.transpose() * + lambdaCd.segment(phiIdx, dRphi_phi.cols()); // Total gradient + tms2["10-solveCxAoA"].stop(); + + //**************************************************************************// + // Gradients wrt mesh coordinates // + //--------------------------------------------------------------------------// + // Ck (cL, cd) // + // Solution lambdaCl_x of (from augmented Lagrangian, eqn d/dx ) // + // "dCk_x + dRphi_x * lambdaCk_phi + dRdStar_x * lambdaCk_x // + // + dRblw_x * lambdaCk_blw + dRx_x * lambdaCk_x = 0" // + //**************************************************************************// + + tms2["12-solveCxMesh"].start(); + Eigen::VectorXd rhsCl_x = dCl_x; + rhsCl_x -= + dRphi_x.transpose() * + lambdaCl.segment(phiIdx, dRphi_phi.cols()); // Potential contribution + rhsCl_x -= + dRM_x.transpose() * + lambdaCl.segment(machIdx, dRM_M.cols()); // Mach number contribution + rhsCl_x -= dRrho_x.transpose() * + lambdaCl.segment(rhoIdx, dRrho_rho.cols()); // Density contribution + rhsCl_x -= dRv_x.transpose() * + lambdaCl.segment(vIdx, dRv_v.cols()); // Velocity contribution + rhsCl_x -= dRdStar_x.transpose() * + lambdaCl.segment( + dStarIdx, + dRdStar_dStar.cols()); // Displacement thickness contribution + rhsCl_x -= dRblw_x.transpose() * + lambdaCl.segment(blwIdx, dRblw_blw.cols()); // Blowing contribution + + Eigen::VectorXd lambdaCl_x = + Eigen::VectorXd::Zero(isol->pbl->nDim * isol->pbl->msh->nodes.size()); + Eigen::Map<Eigen::VectorXd> lambdaCl_x_(lambdaCl_x.data(), lambdaCl_x.size()); + alinsol->compute(dRx_x.transpose(), + Eigen::Map<Eigen::VectorXd>(rhsCl_x.data(), rhsCl_x.size()), + lambdaCl_x_); + + Eigen::VectorXd rhsCd_x = dCdp_x; // Pressure drag coefficient contribution + rhsCd_x += dCdf_x; // Friction drag coefficient contribution + rhsCd_x -= + dRphi_x.transpose() * + lambdaCd.segment(phiIdx, dRphi_phi.cols()); // Potential contribution + rhsCd_x -= + dRM_x.transpose() * + lambdaCd.segment(machIdx, dRM_M.cols()); // Mach number contribution + rhsCd_x -= dRrho_x.transpose() * + lambdaCd.segment(rhoIdx, dRrho_rho.cols()); // Density contribution + rhsCd_x -= dRv_x.transpose() * + lambdaCd.segment(vIdx, dRv_v.cols()); // Velocity contribution + rhsCd_x -= dRdStar_x.transpose() * + lambdaCd.segment( + dStarIdx, + dRdStar_dStar.cols()); // Displacement thickness contribution + rhsCd_x -= dRblw_x.transpose() * + lambdaCd.segment(blwIdx, dRblw_blw.cols()); // Blowing contribution + + Eigen::VectorXd lambdaCd_x = + Eigen::VectorXd::Zero(isol->pbl->nDim * isol->pbl->msh->nodes.size()); + Eigen::Map<Eigen::VectorXd> lambdaCd_x_(lambdaCd_x.data(), lambdaCd_x.size()); + alinsol->compute(dRx_x.transpose(), + Eigen::Map<Eigen::VectorXd>(rhsCd_x.data(), rhsCd_x.size()), + lambdaCd_x_); + + for (auto n : isol->pbl->msh->nodes) + { + for (int m = 0; m < isol->pbl->nDim; m++) + { + tdCl_x[n->row](m) = lambdaCl_x[isol->pbl->nDim * (n->row) + m]; + tdCd_x[n->row](m) = lambdaCd_x[isol->pbl->nDim * (n->row) + m]; + } + } + tms2["12-solveCxMesh"].stop(); + tms["2-Solve"].stop(); + + // Print output + std::cout << "-------------- Adjoint solution --------------" << std::endl; + // aoa gradients + std::cout << std::setw(15) << std::left << "Adjoint AoA" << std::setw(15) + << std::right << "Res[Cl_AoA]" << std::setw(15) << std::right + << "Res[Cd_AoA]" << std::endl; + std::cout << std::fixed << std::setprecision(5); + std::cout << std::setw(15) << std::left << "" << std::setw(15) << std::right + << log10((A_adjoint * Eigen::Map<Eigen::VectorXd>(lambdaCl.data(), + lambdaCl.size()) - + Eigen::Map<Eigen::VectorXd>(rhsCl.data(), rhsCl.size())) + .norm()) + << std::setw(15) << std::right + << log10((A_adjoint * Eigen::Map<Eigen::VectorXd>(lambdaCd.data(), + lambdaCd.size()) - + Eigen::Map<Eigen::VectorXd>(rhsCd.data(), rhsCd.size())) + .norm()) + << std::endl; + std::cout << std::setw(15) << std::left << "AoA gradients" << std::setw(15) + << std::right << "dCl_dAoA" << std::setw(15) << std::right + << "dCd_dAoA" << std::endl; + std::cout << std::fixed << std::setprecision(5); + std::cout << std::setw(15) << std::left << "" << std::setw(15) << std::right + << tdCl_AoA << std::setw(15) << std::right << tdCd_AoA << std::endl; + // mesh gradients + std::cout << std::setw(15) << std::left << "Adjoint mesh" << std::setw(15) + << std::right << "Res[Cl_x]" << std::setw(15) << std::right + << "Res[Cd_x]" << std::endl; + std::cout << std::fixed << std::setprecision(5); + std::cout + << std::setw(15) << std::left << "" << std::setw(15) << std::right + << log10((dRx_x.transpose() * Eigen::Map<Eigen::VectorXd>( + lambdaCl_x.data(), lambdaCl_x.size()) - + Eigen::Map<Eigen::VectorXd>(rhsCl_x.data(), rhsCl_x.size())) + .norm()) + << std::setw(15) << std::right + << log10((dRx_x.transpose() * Eigen::Map<Eigen::VectorXd>( + lambdaCd_x.data(), lambdaCd_x.size()) - + Eigen::Map<Eigen::VectorXd>(rhsCd_x.data(), rhsCd_x.size())) + .norm()) + << std::endl; + std::cout << std::setw(15) << std::left << "Mesh gradients" << std::setw(15) + << std::right << "Norm[dCl_dX]" << std::setw(15) << std::right + << "Norm[dCd_dX]" << std::endl; + std::cout << std::fixed << std::setprecision(5); + std::cout << std::setw(15) << std::left << "" << std::setw(15) << std::right + << Eigen::Map<Eigen::VectorXd>(lambdaCl_x.data(), lambdaCl_x.size()) + .norm() + << std::setw(15) << std::right + << Eigen::Map<Eigen::VectorXd>(lambdaCd_x.data(), lambdaCd_x.size()) + .norm() << std::endl; + tms["0-Total"].stop(); + std::cout << "--------------- Adjoint timers ---------------" << std::endl; + std::cout << tms << "----------------------------------------------" + << std::endl; + if (vsol->verbose > 2) + std::cout << tms2 << "----------------------------------------------" + << std::endl; } /** @@ -454,80 +473,85 @@ void CoupledAdjoint::run() { * Non-zero are: dRphi_phi, dRM_phi, dRrho_phi, dRv_phi * Zeros are: dRdStar_phi, dRblw_phi, dRx_phi */ -void CoupledAdjoint::gradientswrtInviscidFlow() { - //**************************************************************************// - // dRphi_phi // - //--------------------------------------------------------------------------// - // Analytical // - //**************************************************************************// - dRphi_phi = iadjoint->getRu_U(); - - //**************************************************************************// - // dRM_phi, dRrho_phi, dRv_phi // - //--------------------------------------------------------------------------// - // Analytical // - //**************************************************************************// - auto pbl = isol->pbl; - - std::deque<Eigen::Triplet<double>> T_dM; - std::deque<Eigen::Triplet<double>> T_drho; - std::deque<Eigen::Triplet<double>> T_dv; - - size_t jj = 0; - int iReg = 0; - for (auto sec : vsol->sections[0]) { - if (sec->name == "wake") - iReg = 1; - else if (sec->name == "upper" || sec->name == "lower") - iReg = 0; - else - throw std::runtime_error( - "gradientswrtInviscidFlow: Unknown section name"); - for (size_t iNode = 0; iNode < sec->nNodes; ++iNode) { - tbox::Node *srfNode = pbl->bBCs[iReg]->nodes[sec->rows[iNode]]; - double nAdjElms = static_cast<double>(pbl->fluid->neMap[srfNode].size()); - for (auto e : pbl->fluid->neMap[srfNode]) { - // dv/dphi = dv/dgradPhi * dgradPhi/dphi = dgradPhi/dphi = d(gradN_k - // phi_k) / dphi_j dM/dphi = dM/dgradPhi * dgradPhi/dphi = dM/dgradPhi * - // dv/dphi drho/dphi = drho/dgradPhi * dgradPhi/dphi = drho/dgradPhi * - // dv/dphi - Eigen::MatrixXd gradV = e->getJinv(0) * e->getVCache().getDsf(0); - Eigen::VectorXd gradM = - isol->pbl->fluid->mach->evalGrad(*e, isol->phi, 0) * - e->computeGradient(isol->phi, 0).transpose() * gradV; - Eigen::VectorXd gradrho = - isol->pbl->fluid->rho->evalGrad(*e, isol->phi, 0) * - e->computeGradient(isol->phi, 0).transpose() * gradV; - for (size_t k = 0; k < e->nodes.size(); ++k) { - T_dM.push_back(Eigen::Triplet<double>(jj, e->nodes[k]->row, - -gradM(k) / nAdjElms)); - T_drho.push_back(Eigen::Triplet<double>(jj, e->nodes[k]->row, - -gradrho(k) / nAdjElms)); - for (int idim = 0; idim < pbl->nDim; ++idim) - T_dv.push_back(Eigen::Triplet<double>( - jj * pbl->nDim + idim, e->nodes[k]->row, - -gradV(k * pbl->nDim + idim) / nAdjElms)); +void CoupledAdjoint::gradientswrtInviscidFlow() +{ + //**************************************************************************// + // dRphi_phi // + //--------------------------------------------------------------------------// + // Analytical // + //**************************************************************************// + dRphi_phi = iadjoint->getRu_U(); + + //**************************************************************************// + // dRM_phi, dRrho_phi, dRv_phi // + //--------------------------------------------------------------------------// + // Analytical // + //**************************************************************************// + auto pbl = isol->pbl; + + std::deque<Eigen::Triplet<double>> T_dM; + std::deque<Eigen::Triplet<double>> T_drho; + std::deque<Eigen::Triplet<double>> T_dv; + + size_t jj = 0; + int iReg = 0; + for (auto sec : vsol->sections[0]) + { + if (sec->name == "wake") + iReg = 1; + else if (sec->name == "upper" || sec->name == "lower") + iReg = 0; + else + throw std::runtime_error( + "gradientswrtInviscidFlow: Unknown section name"); + for (size_t iNode = 0; iNode < sec->nNodes; ++iNode) + { + tbox::Node *srfNode = pbl->bBCs[iReg]->nodes[sec->rows[iNode]]; + double nAdjElms = static_cast<double>(pbl->fluid->neMap[srfNode].size()); + for (auto e : pbl->fluid->neMap[srfNode]) + { + // dv/dphi = dv/dgradPhi * dgradPhi/dphi = dgradPhi/dphi = d(gradN_k + // phi_k) / dphi_j dM/dphi = dM/dgradPhi * dgradPhi/dphi = dM/dgradPhi * + // dv/dphi drho/dphi = drho/dgradPhi * dgradPhi/dphi = drho/dgradPhi * + // dv/dphi + Eigen::MatrixXd gradV = e->getJinv(0) * e->getVCache().getDsf(0); + Eigen::VectorXd gradM = + isol->pbl->fluid->mach->evalGrad(*e, isol->phi, 0) * + e->computeGradient(isol->phi, 0).transpose() * gradV; + Eigen::VectorXd gradrho = + isol->pbl->fluid->rho->evalGrad(*e, isol->phi, 0) * + e->computeGradient(isol->phi, 0).transpose() * gradV; + for (size_t k = 0; k < e->nodes.size(); ++k) + { + T_dM.push_back(Eigen::Triplet<double>(jj, e->nodes[k]->row, + -gradM(k) / nAdjElms)); + T_drho.push_back(Eigen::Triplet<double>(jj, e->nodes[k]->row, + -gradrho(k) / nAdjElms)); + for (int idim = 0; idim < pbl->nDim; ++idim) + T_dv.push_back(Eigen::Triplet<double>( + jj * pbl->nDim + idim, e->nodes[k]->row, + -gradV(k * pbl->nDim + idim) / nAdjElms)); + } + } + ++jj; } - } - ++jj; } - } - dRM_phi.setFromTriplets(T_dM.begin(), T_dM.end()); - dRrho_phi.setFromTriplets(T_drho.begin(), T_drho.end()); - dRv_phi.setFromTriplets(T_dv.begin(), T_dv.end()); - dRM_phi.prune(0.); - dRrho_phi.prune(0.); - dRv_phi.prune(0.); - - //**************************************************************************// - // dCl_phi, dCd_phi // - //--------------------------------------------------------------------------// - // Analytical // - //**************************************************************************// - std::vector<std::vector<double>> dCx_phi = - iadjoint->computeGradientCoefficientsFlow(); - dCl_phi = Eigen::VectorXd::Map(dCx_phi[0].data(), dCx_phi[0].size()); - dCd_phi = Eigen::VectorXd::Map(dCx_phi[1].data(), dCx_phi[1].size()); + dRM_phi.setFromTriplets(T_dM.begin(), T_dM.end()); + dRrho_phi.setFromTriplets(T_drho.begin(), T_drho.end()); + dRv_phi.setFromTriplets(T_dv.begin(), T_dv.end()); + dRM_phi.prune(0.); + dRrho_phi.prune(0.); + dRv_phi.prune(0.); + + //**************************************************************************// + // dCl_phi, dCd_phi // + //--------------------------------------------------------------------------// + // Analytical // + //**************************************************************************// + std::vector<std::vector<double>> dCx_phi = + iadjoint->computeGradientCoefficientsFlow(); + dCl_phi = Eigen::VectorXd::Map(dCx_phi[0].data(), dCx_phi[0].size()); + dCd_phi = Eigen::VectorXd::Map(dCx_phi[1].data(), dCx_phi[1].size()); } /** @@ -535,46 +559,48 @@ void CoupledAdjoint::gradientswrtInviscidFlow() { * Non-zero are: dRM_M, dRdStar_M * Zeros are: dRphi_M, dRrho_M, dRv_M, dRblw_M, dRx_M */ -void CoupledAdjoint::gradientswrtMachNumber() { - //**************************************************************************// - // dRM_M // - //**************************************************************************// - dRM_M.setIdentity(); - - //**************************************************************************// - // dRdStar_M // - //--------------------------------------------------------------------------// - // Central finite-difference // - //**************************************************************************// - std::deque<Eigen::Triplet<double>> T; - double delta = 1e-6; - size_t i = 0; - double saveM = 0.; - std::vector<Eigen::VectorXd> ddStar(2, - Eigen::VectorXd::Zero(dRdStar_M.cols())); - std::vector<double> dCdf(2, 0.); - for (auto sec : vsol->sections[0]) - for (size_t j = 0; j < sec->nNodes; ++j) { - saveM = sec->Me[j]; - - sec->Me[j] = saveM + delta; - ddStar[0] = this->runViscous(); - dCdf[0] = vsol->Cdf; - sec->Me[j] = saveM - delta; - ddStar[1] = this->runViscous(); - dCdf[1] = vsol->Cdf; - - sec->Me[j] = saveM; - - for (int k = 0; k < ddStar[0].size(); ++k) - T.push_back(Eigen::Triplet<double>( - k, i, -(ddStar[0][k] - ddStar[1][k]) / (2. * delta))); - // Collect Cd gradients - dCd_M(i) = (dCdf[0] - dCdf[1]) / (2. * delta); - ++i; - } - dRdStar_M.setFromTriplets(T.begin(), T.end()); - dRdStar_M.prune(0.); +void CoupledAdjoint::gradientswrtMachNumber() +{ + //**************************************************************************// + // dRM_M // + //**************************************************************************// + dRM_M.setIdentity(); + + //**************************************************************************// + // dRdStar_M // + //--------------------------------------------------------------------------// + // Central finite-difference // + //**************************************************************************// + std::deque<Eigen::Triplet<double>> T; + double delta = 1e-6; + size_t i = 0; + double saveM = 0.; + std::vector<Eigen::VectorXd> ddStar(2, + Eigen::VectorXd::Zero(dRdStar_M.cols())); + std::vector<double> dCdf(2, 0.); + for (auto sec : vsol->sections[0]) + for (size_t j = 0; j < sec->nNodes; ++j) + { + saveM = sec->Me[j]; + + sec->Me[j] = saveM + delta; + ddStar[0] = this->runViscous(); + dCdf[0] = vsol->Cdf; + sec->Me[j] = saveM - delta; + ddStar[1] = this->runViscous(); + dCdf[1] = vsol->Cdf; + + sec->Me[j] = saveM; + + for (int k = 0; k < ddStar[0].size(); ++k) + T.push_back(Eigen::Triplet<double>( + k, i, -(ddStar[0][k] - ddStar[1][k]) / (2. * delta))); + // Collect Cd gradients + dCd_M(i) = (dCdf[0] - dCdf[1]) / (2. * delta); + ++i; + } + dRdStar_M.setFromTriplets(T.begin(), T.end()); + dRdStar_M.prune(0.); } /** @@ -582,39 +608,41 @@ void CoupledAdjoint::gradientswrtMachNumber() { * Non-zero are: dRrho_rho, dRblw_rho * Zeros are: dRphi_rho, dRM_rho, dRv_rho, dRdStar_rho, dRx_rho */ -void CoupledAdjoint::gradientswrtDensity() { - size_t nNs = 0; - for (auto bBC : isol->pbl->bBCs) - nNs += bBC->nodes.size(); // Number of surface nodes - nNs += 1; // Duplicate stagnation point - size_t nEs = 0; - for (auto bBC : isol->pbl->bBCs) - nEs += bBC->uE.size(); // Number of blowing elements - - //**************************************************************************// - // dRrho_rho // - //**************************************************************************// - dRrho_rho.setIdentity(); - - //**************************************************************************// - // dRblw_rho // - //--------------------------------------------------------------------------// - // Analytical // - //**************************************************************************// - std::deque<Eigen::Triplet<double>> T_blw; - int offSetElms = 0; - int offSetNodes = 0; - for (const auto &sec : vsol->sections[0]) { - std::vector<std::vector<double>> grad = sec->evalGradwrtRho(); - for (size_t i = 0; i < grad.size(); ++i) - for (size_t j = 0; j < grad[i].size(); ++j) - T_blw.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes, - -grad[i][j])); - offSetElms += sec->nElms; - offSetNodes += sec->nNodes; - } - dRblw_rho.setFromTriplets(T_blw.begin(), T_blw.end()); - dRblw_rho.prune(0.); +void CoupledAdjoint::gradientswrtDensity() +{ + size_t nNs = 0; + for (auto bBC : isol->pbl->bBCs) + nNs += bBC->nodes.size(); // Number of surface nodes + nNs += 1; // Duplicate stagnation point + size_t nEs = 0; + for (auto bBC : isol->pbl->bBCs) + nEs += bBC->uE.size(); // Number of blowing elements + + //**************************************************************************// + // dRrho_rho // + //**************************************************************************// + dRrho_rho.setIdentity(); + + //**************************************************************************// + // dRblw_rho // + //--------------------------------------------------------------------------// + // Analytical // + //**************************************************************************// + std::deque<Eigen::Triplet<double>> T_blw; + int offSetElms = 0; + int offSetNodes = 0; + for (const auto &sec : vsol->sections[0]) + { + std::vector<std::vector<double>> grad = sec->evalGradwrtRho(); + for (size_t i = 0; i < grad.size(); ++i) + for (size_t j = 0; j < grad[i].size(); ++j) + T_blw.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes, + -grad[i][j])); + offSetElms += sec->nElms; + offSetNodes += sec->nNodes; + } + dRblw_rho.setFromTriplets(T_blw.begin(), T_blw.end()); + dRblw_rho.prune(0.); } /** @@ -622,137 +650,145 @@ void CoupledAdjoint::gradientswrtDensity() { * Non-zeros are: dRv_v, dRdStar_v, dRblw_v * Zeros are: dRphi_v, dRM_v, dRrho_v, dRx_v */ -void CoupledAdjoint::gradientswrtVelocity() { - size_t nDim = isol->pbl->nDim; // Number of dimensions of the problem - size_t nNs = 0; - for (auto bBC : isol->pbl->bBCs) - nNs += bBC->nodes.size(); // Number of surface nodes - nNs += 1; // Duplicate stagnation point - size_t nEs = 0; - for (auto bBC : isol->pbl->bBCs) - nEs += bBC->uE.size(); // Number of blowing elements - - //**************************************************************************// - // dRv_v // - //**************************************************************************// - dRv_v.setIdentity(); - - //**** Get velocity****/ - Eigen::MatrixXd v = Eigen::MatrixXd::Zero(nNs, nDim); - size_t i = 0; - int iReg = -1; - for (const auto &sec : vsol->sections[0]) { - if (sec->name == "wake") - iReg = 1; - else if (sec->name == "upper" || sec->name == "lower") - iReg = 0; - else - throw std::runtime_error("Unknown section name"); - for (size_t j = 0; j < sec->nNodes; ++j) { - for (size_t idim = 0; idim < nDim; ++idim) - v(i, idim) = - isol->U[isol->pbl->bBCs[iReg]->nodes[sec->rows[j]]->row](idim); - ++i; +void CoupledAdjoint::gradientswrtVelocity() +{ + size_t nDim = isol->pbl->nDim; // Number of dimensions of the problem + size_t nNs = 0; + for (auto bBC : isol->pbl->bBCs) + nNs += bBC->nodes.size(); // Number of surface nodes + nNs += 1; // Duplicate stagnation point + size_t nEs = 0; + for (auto bBC : isol->pbl->bBCs) + nEs += bBC->uE.size(); // Number of blowing elements + + //**************************************************************************// + // dRv_v // + //**************************************************************************// + dRv_v.setIdentity(); + + //**** Get velocity****/ + Eigen::MatrixXd v = Eigen::MatrixXd::Zero(nNs, nDim); + size_t i = 0; + int iReg = -1; + for (const auto &sec : vsol->sections[0]) + { + if (sec->name == "wake") + iReg = 1; + else if (sec->name == "upper" || sec->name == "lower") + iReg = 0; + else + throw std::runtime_error("Unknown section name"); + for (size_t j = 0; j < sec->nNodes; ++j) + { + for (size_t idim = 0; idim < nDim; ++idim) + v(i, idim) = + isol->U[isol->pbl->bBCs[iReg]->nodes[sec->rows[j]]->row](idim); + ++i; + } } - } - - //**************************************************************************// - // dvt_v // - //--------------------------------------------------------------------------// - // Analytical // - //**************************************************************************// - Eigen::SparseMatrix<double, Eigen::RowMajor> dVt_v(nNs, nNs * nDim); - std::deque<Eigen::Triplet<double>> T_vt; - - i = 0; - for (const auto &sec : vsol->sections[0]) { - for (size_t j = 0; j < sec->nNodes; ++j) { - for (size_t idim = 0; idim < nDim; ++idim) { - T_vt.push_back(Eigen::Triplet<double>( - i, i * nDim + idim, - 1 / std::sqrt(v(i, 0) * v(i, 0) + v(i, 1) * v(i, 1)) * v(i, idim))); - } - ++i; + + //**************************************************************************// + // dvt_v // + //--------------------------------------------------------------------------// + // Analytical // + //**************************************************************************// + Eigen::SparseMatrix<double, Eigen::RowMajor> dVt_v(nNs, nNs * nDim); + std::deque<Eigen::Triplet<double>> T_vt; + + i = 0; + for (const auto &sec : vsol->sections[0]) + { + for (size_t j = 0; j < sec->nNodes; ++j) + { + for (size_t idim = 0; idim < nDim; ++idim) + { + T_vt.push_back(Eigen::Triplet<double>( + i, i * nDim + idim, + 1 / std::sqrt(v(i, 0) * v(i, 0) + v(i, 1) * v(i, 1)) * v(i, idim))); + } + ++i; + } } - } - dVt_v.setFromTriplets(T_vt.begin(), T_vt.end()); - dVt_v.prune(0.); - - //**************************************************************************// - // dRdStar_vt // - //--------------------------------------------------------------------------// - // Central finite-difference // - //**************************************************************************// - std::deque<Eigen::Triplet<double>> T; - double delta = 1e-6; - i = 0; - double saveM = 0.; - std::vector<Eigen::VectorXd> dvt(2, Eigen::VectorXd::Zero(dRdStar_M.cols())); - std::vector<double> dCdf(2, 0.); - Eigen::VectorXd dCd_vt = Eigen::VectorXd::Zero(dRM_M.cols()); - for (auto sec : vsol->sections[0]) - for (size_t j = 0; j < sec->nNodes; ++j) { - saveM = sec->vt[j]; - - sec->vt[j] = saveM + delta; - dvt[0] = this->runViscous(); - dCdf[0] = vsol->Cdf; - sec->vt[j] = saveM - delta; - dvt[1] = this->runViscous(); - dCdf[1] = vsol->Cdf; - - sec->vt[j] = saveM; - - for (int k = 0; k < dvt[0].size(); ++k) - T.push_back(Eigen::Triplet<double>( - k, i, -(dvt[0][k] - dvt[1][k]) / (2. * delta))); - dCd_vt(i) = (dCdf[0] - dCdf[1]) / (2. * delta); - ++i; + dVt_v.setFromTriplets(T_vt.begin(), T_vt.end()); + dVt_v.prune(0.); + + //**************************************************************************// + // dRdStar_vt // + //--------------------------------------------------------------------------// + // Central finite-difference // + //**************************************************************************// + std::deque<Eigen::Triplet<double>> T; + double delta = 1e-6; + i = 0; + double saveM = 0.; + std::vector<Eigen::VectorXd> dvt(2, Eigen::VectorXd::Zero(dRdStar_M.cols())); + std::vector<double> dCdf(2, 0.); + Eigen::VectorXd dCd_vt = Eigen::VectorXd::Zero(dRM_M.cols()); + for (auto sec : vsol->sections[0]) + for (size_t j = 0; j < sec->nNodes; ++j) + { + saveM = sec->vt[j]; + + sec->vt[j] = saveM + delta; + dvt[0] = this->runViscous(); + dCdf[0] = vsol->Cdf; + sec->vt[j] = saveM - delta; + dvt[1] = this->runViscous(); + dCdf[1] = vsol->Cdf; + + sec->vt[j] = saveM; + + for (int k = 0; k < dvt[0].size(); ++k) + T.push_back(Eigen::Triplet<double>( + k, i, -(dvt[0][k] - dvt[1][k]) / (2. * delta))); + dCd_vt(i) = (dCdf[0] - dCdf[1]) / (2. * delta); + ++i; + } + Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_vt(dRdStar_M.rows(), + dRdStar_M.cols()); + dRdStar_vt.setFromTriplets(T.begin(), T.end()); + dRdStar_vt.prune(0.); + + //**************************************************************************// + // dRdStar_v // + //**************************************************************************// + dRdStar_v = dRdStar_vt * dVt_v; + dRdStar_v.prune(0.); + + //**************************************************************************// + // dCdf_v // + //**************************************************************************// + dCd_v = + dCd_vt.transpose() * dVt_v; // [1xnNs] x [nNs x nNs*nDim] = [1 x nNs*nDim] + + //**************************************************************************// + // dRblw_vt // + //--------------------------------------------------------------------------// + // Analytical // + //**************************************************************************// + std::deque<Eigen::Triplet<double>> T_blw; + int offSetElms = 0; + int offSetNodes = 0; + for (const auto &sec : vsol->sections[0]) + { + std::vector<std::vector<double>> grad = sec->evalGradwrtVt(); + for (size_t i = 0; i < grad.size(); ++i) + for (size_t j = 0; j < grad[i].size(); ++j) + T_blw.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes, + -grad[i][j])); + offSetElms += sec->nElms; + offSetNodes += sec->nNodes; } - Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_vt(dRdStar_M.rows(), - dRdStar_M.cols()); - dRdStar_vt.setFromTriplets(T.begin(), T.end()); - dRdStar_vt.prune(0.); - - //**************************************************************************// - // dRdStar_v // - //**************************************************************************// - dRdStar_v = dRdStar_vt * dVt_v; - dRdStar_v.prune(0.); - - //**************************************************************************// - // dCdf_v // - //**************************************************************************// - dCd_v = - dCd_vt.transpose() * dVt_v; // [1xnNs] x [nNs x nNs*nDim] = [1 x nNs*nDim] - - //**************************************************************************// - // dRblw_vt // - //--------------------------------------------------------------------------// - // Analytical // - //**************************************************************************// - std::deque<Eigen::Triplet<double>> T_blw; - int offSetElms = 0; - int offSetNodes = 0; - for (const auto &sec : vsol->sections[0]) { - std::vector<std::vector<double>> grad = sec->evalGradwrtVt(); - for (size_t i = 0; i < grad.size(); ++i) - for (size_t j = 0; j < grad[i].size(); ++j) - T_blw.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes, - -grad[i][j])); - offSetElms += sec->nElms; - offSetNodes += sec->nNodes; - } - - Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_vt(nEs, nNs); - dRblw_vt.setFromTriplets(T_blw.begin(), T_blw.end()); - dRblw_vt.prune(0.); - - //**************************************************************************// - // dRblw_v // - //**************************************************************************// - dRblw_v = dRblw_vt * dVt_v; - dRblw_v.prune(0.); + + Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_vt(nEs, nNs); + dRblw_vt.setFromTriplets(T_blw.begin(), T_blw.end()); + dRblw_vt.prune(0.); + + //**************************************************************************// + // dRblw_v // + //**************************************************************************// + dRblw_v = dRblw_vt * dVt_v; + dRblw_v.prune(0.); } /** @@ -760,67 +796,70 @@ void CoupledAdjoint::gradientswrtVelocity() { * surface Non-zero are: dRdStar_dStar, dRblw_dStar Zeros are: dRphi_dStar, * dRM_dStar, dRrho_dStar, dRv_dStar, dRx_dStar */ -void CoupledAdjoint::gradientswrtDeltaStar() { - //**************************************************************************// - // dRdStar_dStar // - //--------------------------------------------------------------------------// - // Central finite-difference // - //**************************************************************************// - std::deque<Eigen::Triplet<double>> T; - double delta = 1e-6; - double saveDStar = 0.; - std::vector<Eigen::VectorXd> ddstar( - 2, Eigen::VectorXd::Zero(dRdStar_dStar.cols())); - std::vector<double> dCdf(2, 0.); - size_t i = 0; - for (auto sec : vsol->sections[0]) - for (size_t j = 0; j < sec->nNodes; ++j) { - saveDStar = sec->deltaStarExt[j]; - - sec->deltaStarExt[j] = saveDStar + delta; - ddstar[0] = this->runViscous(); - dCdf[0] = vsol->Cdf; - sec->deltaStarExt[j] = saveDStar - delta; - ddstar[1] = this->runViscous(); - dCdf[1] = vsol->Cdf; - - sec->deltaStarExt[j] = saveDStar; - - for (int k = 0; k < ddstar[0].size(); ++k) - T.push_back(Eigen::Triplet<double>( - k, i, (ddstar[0][k] - ddstar[1][k]) / (2. * delta))); - dCd_dStar(i) = (dCdf[0] - dCdf[1]) / (2. * delta); - ++i; +void CoupledAdjoint::gradientswrtDeltaStar() +{ + //**************************************************************************// + // dRdStar_dStar // + //--------------------------------------------------------------------------// + // Central finite-difference // + //**************************************************************************// + std::deque<Eigen::Triplet<double>> T; + double delta = 1e-6; + double saveDStar = 0.; + std::vector<Eigen::VectorXd> ddstar( + 2, Eigen::VectorXd::Zero(dRdStar_dStar.cols())); + std::vector<double> dCdf(2, 0.); + size_t i = 0; + for (auto sec : vsol->sections[0]) + for (size_t j = 0; j < sec->nNodes; ++j) + { + saveDStar = sec->deltaStarExt[j]; + + sec->deltaStarExt[j] = saveDStar + delta; + ddstar[0] = this->runViscous(); + dCdf[0] = vsol->Cdf; + sec->deltaStarExt[j] = saveDStar - delta; + ddstar[1] = this->runViscous(); + dCdf[1] = vsol->Cdf; + + sec->deltaStarExt[j] = saveDStar; + + for (int k = 0; k < ddstar[0].size(); ++k) + T.push_back(Eigen::Triplet<double>( + k, i, (ddstar[0][k] - ddstar[1][k]) / (2. * delta))); + dCd_dStar(i) = (dCdf[0] - dCdf[1]) / (2. * delta); + ++i; + } + Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_dStar_dum( + dRdStar_dStar.rows(), dRdStar_dStar.cols()); + // RdStar = dStar - f(dStar) = 0 + // dRdStar_dStar = 1 - df_dStar (minus is 3 lines below and not in the + // triplets) + dRdStar_dStar_dum.setFromTriplets(T.begin(), T.end()); + dRdStar_dStar.setIdentity(); + dRdStar_dStar -= dRdStar_dStar_dum; + dRdStar_dStar.prune(0.); + + //**************************************************************************// + // dRblw_dStar // + //--------------------------------------------------------------------------// + // Analytical // + //**************************************************************************// + std::deque<Eigen::Triplet<double>> T_blw; + int offSetElms = 0; + int offSetNodes = 0; + for (const auto &sec : vsol->sections[0]) + { + std::vector<std::vector<double>> grad = sec->evalGradwrtDeltaStar(); + for (size_t i = 0; i < grad.size(); ++i) + for (size_t j = 0; j < grad[i].size(); ++j) + T_blw.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes, + -grad[i][j])); + offSetElms += sec->nElms; + offSetNodes += sec->nNodes; } - Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_dStar_dum( - dRdStar_dStar.rows(), dRdStar_dStar.cols()); - // RdStar = dStar - f(dStar) = 0 - // dRdStar_dStar = 1 - df_dStar (minus is 3 lines below and not in the - // triplets) - dRdStar_dStar_dum.setFromTriplets(T.begin(), T.end()); - dRdStar_dStar.setIdentity(); - dRdStar_dStar -= dRdStar_dStar_dum; - dRdStar_dStar.prune(0.); - - //**************************************************************************// - // dRblw_dStar // - //--------------------------------------------------------------------------// - // Analytical // - //**************************************************************************// - std::deque<Eigen::Triplet<double>> T_blw; - int offSetElms = 0; - int offSetNodes = 0; - for (const auto &sec : vsol->sections[0]) { - std::vector<std::vector<double>> grad = sec->evalGradwrtDeltaStar(); - for (size_t i = 0; i < grad.size(); ++i) - for (size_t j = 0; j < grad[i].size(); ++j) - T_blw.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes, - -grad[i][j])); - offSetElms += sec->nElms; - offSetNodes += sec->nNodes; - } - dRblw_dStar.setFromTriplets(T_blw.begin(), T_blw.end()); - dRblw_dStar.prune(0.); + dRblw_dStar.setFromTriplets(T_blw.begin(), T_blw.end()); + dRblw_dStar.prune(0.); } /** @@ -828,41 +867,45 @@ void CoupledAdjoint::gradientswrtDeltaStar() { * Non-zero are: dRphi_blw, dRblw_blw * Zeros are: dRM_blw, dRrho_blw, dRv_blw, dRdStar_blw, dRx_blw */ -void CoupledAdjoint::gradientswrtBlowingVelocity() { - //**************************************************************************// - // dRphi_blw // - //--------------------------------------------------------------------------// - // Analytical // - //**************************************************************************// - std::deque<Eigen::Triplet<double>> T; - - size_t jj = 0; - int iReg = 0; - for (auto sec : vsol->sections[0]) { - if (sec->name == "wake") - iReg = 1; - else if (sec->name == "upper" || sec->name == "lower") - iReg = 0; - else - throw std::runtime_error("Unknown section name"); - for (size_t iElm = 0; iElm < sec->nElms; ++iElm) { - auto blowingElement = isol->pbl->bBCs[iReg]->uE[sec->no[iElm]].first; - Eigen::VectorXd be = - dart::BlowingResidual::buildGradientBlowing(*blowingElement); - for (size_t ii = 0; ii < blowingElement->nodes.size(); ii++) { - tbox::Node *nodi = blowingElement->nodes[ii]; - T.push_back( - Eigen::Triplet<double>(isol->getRow(nodi->row), jj, be(ii))); - } - ++jj; +void CoupledAdjoint::gradientswrtBlowingVelocity() +{ + //**************************************************************************// + // dRphi_blw // + //--------------------------------------------------------------------------// + // Analytical // + //**************************************************************************// + std::deque<Eigen::Triplet<double>> T; + + size_t jj = 0; + int iReg = 0; + for (auto sec : vsol->sections[0]) + { + if (sec->name == "wake") + iReg = 1; + else if (sec->name == "upper" || sec->name == "lower") + iReg = 0; + else + throw std::runtime_error("Unknown section name"); + for (size_t iElm = 0; iElm < sec->nElms; ++iElm) + { + auto blowingElement = isol->pbl->bBCs[iReg]->uE[sec->no[iElm]].first; + Eigen::VectorXd be = + dart::BlowingResidual::buildGradientBlowing(*blowingElement); + for (size_t ii = 0; ii < blowingElement->nodes.size(); ii++) + { + tbox::Node *nodi = blowingElement->nodes[ii]; + T.push_back( + Eigen::Triplet<double>(isol->getRow(nodi->row), jj, be(ii))); + } + ++jj; + } } - } - dRphi_blw.setFromTriplets(T.begin(), T.end()); - dRphi_blw.prune(0.); - dRphi_blw.makeCompressed(); + dRphi_blw.setFromTriplets(T.begin(), T.end()); + dRphi_blw.prune(0.); + dRphi_blw.makeCompressed(); - /**** dRblw_blw ****/ - dRblw_blw.setIdentity(); + /**** dRblw_blw ****/ + dRblw_blw.setIdentity(); } /** @@ -870,22 +913,23 @@ void CoupledAdjoint::gradientswrtBlowingVelocity() { * Non-zero are: dRphi_AoA * Zeros are: dRM_AoA, dRrho_AoA, dRv_AoA, dRdStar_AoA, dRblw_AoA, dRx_AoA */ -void CoupledAdjoint::gradientwrtAoA() { - //**************************************************************************// - // dCl_Aoa dCd_Aoa // - //--------------------------------------------------------------------------// - // Analytical // - //**************************************************************************// - std::vector<double> dCx_AoA = iadjoint->computeGradientCoefficientsAoa(); - dCl_AoA = dCx_AoA[0]; - dCd_AoA = dCx_AoA[1]; - - //**************************************************************************// - // dRphi_Aoa // - //--------------------------------------------------------------------------// - // Analytical // - //**************************************************************************// - dRphi_AoA = iadjoint->getRu_A(); +void CoupledAdjoint::gradientwrtAoA() +{ + //**************************************************************************// + // dCl_Aoa dCd_Aoa // + //--------------------------------------------------------------------------// + // Analytical // + //**************************************************************************// + std::vector<double> dCx_AoA = iadjoint->computeGradientCoefficientsAoa(); + dCl_AoA = dCx_AoA[0]; + dCd_AoA = dCx_AoA[1]; + + //**************************************************************************// + // dRphi_Aoa // + //--------------------------------------------------------------------------// + // Analytical // + //**************************************************************************// + dRphi_AoA = iadjoint->getRu_A(); } /** @@ -893,198 +937,207 @@ void CoupledAdjoint::gradientwrtAoA() { * Non-zero are: dRphi_x, dRM_x, dRrho_x, dRv_x, dRdStar_x, dRblw_x, dRx_x * Zeros are: - */ -void CoupledAdjoint::gradientwrtMesh() { - - int nDim = isol->pbl->nDim; - - //**************************************************************************// - // dCl_x, dCd_x // - //--------------------------------------------------------------------------// - // Analytical // - //**************************************************************************// - std::vector<std::vector<double>> dCx_x = - iadjoint->computeGradientCoefficientsMesh(); - dCl_x = Eigen::Map<Eigen::VectorXd>(dCx_x[0].data(), dCx_x[0].size()); - dCdp_x = Eigen::Map<Eigen::VectorXd>(dCx_x[1].data(), dCx_x[1].size()); - - //**************************************************************************// - // dRphi_x // - //--------------------------------------------------------------------------// - // Analytical // - //**************************************************************************// - dRphi_x = iadjoint->getRu_X(); - //**************************************************************************// - // dRx_x // - //--------------------------------------------------------------------------// - // Analytical // - //**************************************************************************// - dRx_x = iadjoint->getRx_X(); - - //**************************************************************************// - // dRdStar_x, dCdf_dx // - //--------------------------------------------------------------------------// - // Central finite-difference // - //--------------------------------------------------------------------------// - // We don't do dRdStar_x = dRdStar_loc * dloc_dx because Cdf = f(xi(x), x) // - // and we would have to compute dCdf/dx and dCdf/dloc. Same for deltaStar. // - // It is more computationally expansive to compute // - // ddeltaStar/dx, ddeltaStar/dloc than ddeltaStar_dx, ddeltaStar_dy // - //**************************************************************************// - std::deque<Eigen::Triplet<double>> T_dStar_x; - double delta = 1e-8; - int iReg = 0; - std::vector<Eigen::VectorXd> ddStar(2, - Eigen::VectorXd::Zero(dRdStar_M.cols())); - std::vector<double> dCdf(2, 0.); - for (auto sec : vsol->sections[0]) { - if (sec->name == "wake") - iReg = 1; - else if (sec->name == "upper" || sec->name == "lower") - iReg = 0; - else - throw std::runtime_error("Unknown section name"); - for (size_t j = 0; j < sec->nNodes; ++j) { - int rowj = isol->pbl->bBCs[iReg]->nodes[sec->rows[j]]->row; - // x - double saveX = sec->x[j]; - - sec->x[j] = saveX + delta; - sec->computeLocalCoordinates(); - sec->xxExt = sec->loc; - ddStar[0] = this->runViscous(); - dCdf[0] = vsol->Cdf; - sec->x[j] = saveX - delta; - sec->computeLocalCoordinates(); - sec->xxExt = sec->loc; - ddStar[1] = this->runViscous(); - dCdf[1] = vsol->Cdf; - - sec->x[j] = saveX; - sec->xxExt = sec->loc; - sec->computeLocalCoordinates(); - - for (int k = 0; k < ddStar[0].size(); ++k) - T_dStar_x.push_back(Eigen::Triplet<double>( - k, rowj * isol->pbl->nDim + 0, - -(ddStar[0][k] - ddStar[1][k]) / (2. * delta))); - dCdf_x(rowj * isol->pbl->nDim + 0) += (dCdf[0] - dCdf[1]) / (2. * delta); - // y - double saveY = sec->y[j]; - - sec->y[j] = saveY + delta; - sec->computeLocalCoordinates(); - sec->xxExt = sec->loc; - ddStar[0] = this->runViscous(); - dCdf[0] = vsol->Cdf; - sec->y[j] = saveY - delta; - sec->computeLocalCoordinates(); - sec->xxExt = sec->loc; - ddStar[1] = this->runViscous(); - dCdf[1] = vsol->Cdf; - - sec->y[j] = saveY; - sec->xxExt = sec->loc; - sec->computeLocalCoordinates(); - - for (int k = 0; k < ddStar[0].size(); ++k) - T_dStar_x.push_back(Eigen::Triplet<double>( - k, rowj * isol->pbl->nDim + 1, - -(ddStar[0][k] - ddStar[1][k]) / (2. * delta))); - dCdf_x(rowj * isol->pbl->nDim + 1) += (dCdf[0] - dCdf[1]) / (2. * delta); +void CoupledAdjoint::gradientwrtMesh() +{ + + int nDim = isol->pbl->nDim; + + //**************************************************************************// + // dCl_x, dCd_x // + //--------------------------------------------------------------------------// + // Analytical // + //**************************************************************************// + std::vector<std::vector<double>> dCx_x = + iadjoint->computeGradientCoefficientsMesh(); + dCl_x = Eigen::Map<Eigen::VectorXd>(dCx_x[0].data(), dCx_x[0].size()); + dCdp_x = Eigen::Map<Eigen::VectorXd>(dCx_x[1].data(), dCx_x[1].size()); + + //**************************************************************************// + // dRphi_x // + //--------------------------------------------------------------------------// + // Analytical // + //**************************************************************************// + dRphi_x = iadjoint->getRu_X(); + //**************************************************************************// + // dRx_x // + //--------------------------------------------------------------------------// + // Analytical // + //**************************************************************************// + dRx_x = iadjoint->getRx_X(); + + //**************************************************************************// + // dRdStar_x, dCdf_dx // + //--------------------------------------------------------------------------// + // Central finite-difference // + //--------------------------------------------------------------------------// + // We don't do dRdStar_x = dRdStar_loc * dloc_dx because Cdf = f(xi(x), x) // + // and we would have to compute dCdf/dx and dCdf/dloc. Same for deltaStar. // + // It is more computationally expansive to compute // + // ddeltaStar/dx, ddeltaStar/dloc than ddeltaStar_dx, ddeltaStar_dy // + //**************************************************************************// + std::deque<Eigen::Triplet<double>> T_dStar_x; + double delta = 1e-8; + int iReg = 0; + std::vector<Eigen::VectorXd> ddStar(2, + Eigen::VectorXd::Zero(dRdStar_M.cols())); + std::vector<double> dCdf(2, 0.); + for (auto sec : vsol->sections[0]) + { + if (sec->name == "wake") + iReg = 1; + else if (sec->name == "upper" || sec->name == "lower") + iReg = 0; + else + throw std::runtime_error("Unknown section name"); + for (size_t j = 0; j < sec->nNodes; ++j) + { + int rowj = isol->pbl->bBCs[iReg]->nodes[sec->rows[j]]->row; + // x + double saveX = sec->x[j]; + + sec->x[j] = saveX + delta; + sec->computeLocalCoordinates(); + sec->xxExt = sec->loc; + ddStar[0] = this->runViscous(); + dCdf[0] = vsol->Cdf; + sec->x[j] = saveX - delta; + sec->computeLocalCoordinates(); + sec->xxExt = sec->loc; + ddStar[1] = this->runViscous(); + dCdf[1] = vsol->Cdf; + + sec->x[j] = saveX; + sec->xxExt = sec->loc; + sec->computeLocalCoordinates(); + + for (int k = 0; k < ddStar[0].size(); ++k) + T_dStar_x.push_back(Eigen::Triplet<double>( + k, rowj * isol->pbl->nDim + 0, + -(ddStar[0][k] - ddStar[1][k]) / (2. * delta))); + dCdf_x(rowj * isol->pbl->nDim + 0) += (dCdf[0] - dCdf[1]) / (2. * delta); + // y + double saveY = sec->y[j]; + + sec->y[j] = saveY + delta; + sec->computeLocalCoordinates(); + sec->xxExt = sec->loc; + ddStar[0] = this->runViscous(); + dCdf[0] = vsol->Cdf; + sec->y[j] = saveY - delta; + sec->computeLocalCoordinates(); + sec->xxExt = sec->loc; + ddStar[1] = this->runViscous(); + dCdf[1] = vsol->Cdf; + + sec->y[j] = saveY; + sec->xxExt = sec->loc; + sec->computeLocalCoordinates(); + + for (int k = 0; k < ddStar[0].size(); ++k) + T_dStar_x.push_back(Eigen::Triplet<double>( + k, rowj * isol->pbl->nDim + 1, + -(ddStar[0][k] - ddStar[1][k]) / (2. * delta))); + dCdf_x(rowj * isol->pbl->nDim + 1) += (dCdf[0] - dCdf[1]) / (2. * delta); + } } - } - dRdStar_x.setFromTriplets(T_dStar_x.begin(), T_dStar_x.end()); - dRdStar_x.prune(0.); - - //**************************************************************************// - // dRblw_x // - //--------------------------------------------------------------------------// - // Analytical // - //**************************************************************************// - std::deque<Eigen::Triplet<double>> T_blw_x; - int offSetElms = 0; - for (const auto &sec : vsol->sections[0]) { - if (sec->name == "wake") - iReg = 1; - else if (sec->name == "upper" || sec->name == "lower") - iReg = 0; - else - throw std::runtime_error("Unknown section name"); - std::vector<std::vector<double>> gradX = sec->evalGradwrtX(); - std::vector<std::vector<double>> gradY = sec->evalGradwrtY(); - for (size_t i = 0; i < gradX.size(); ++i) - for (size_t j = 0; j < gradX[i].size(); ++j) { - int rowj = isol->pbl->bBCs[iReg]->nodes[sec->rows[j]]->row; - T_blw_x.push_back(Eigen::Triplet<double>( - i + offSetElms, rowj * isol->pbl->nDim + 0, -gradX[i][j])); - T_blw_x.push_back(Eigen::Triplet<double>( - i + offSetElms, rowj * isol->pbl->nDim + 1, -gradY[i][j])); - } - offSetElms += sec->nElms; - } - dRblw_x.setFromTriplets(T_blw_x.begin(), T_blw_x.end()); - dRblw_x.prune(0.); - - //**************************************************************************// - // dRM_x, dRrho_x, dRv_x // - //--------------------------------------------------------------------------// - // Analytical // - //**************************************************************************// - auto pbl = isol->pbl; - - std::deque<Eigen::Triplet<double>> T_dM_x; - std::deque<Eigen::Triplet<double>> T_drho_x; - std::deque<Eigen::Triplet<double>> T_dv_x; - - int jj = 0; - iReg = 0; - - for (auto sec : vsol->sections[0]) { - if (sec->name == "wake") - iReg = 1; - else if (sec->name == "upper" || sec->name == "lower") - iReg = 0; - else - throw std::runtime_error( - "gradientswrtInviscidFlow: Unknown section name"); - for (size_t iNode = 0; iNode < sec->nNodes; ++iNode) { - tbox::Node *srfNode = pbl->bBCs[iReg]->nodes[sec->rows[iNode]]; - double nAdjElms = static_cast<double>(pbl->fluid->neMap[srfNode].size()); - for (auto e : pbl->fluid->neMap[srfNode]) { - // dv/dx = dv/dgradPhi * dgradPhi/dx = dgradPhi/dx = (-J^(-1) * - // dJ/dx_j)*grad_xPhi dM/dx = dM/dgradPhi * dgradPhi/dx = dM/dgradPhi * - // gradPhi * dv/dx drho/dx = drho/dgradPhi * dgradPhi/dx = drho/dgradPhi - // * gradPhi * dv/dx - for (size_t inodOnElm = 0; inodOnElm < e->nodes.size(); inodOnElm++) - for (int idim = 0; idim < nDim; ++idim) { - Eigen::VectorXd gradV = (-e->getJinv(0)) * - e->getGradJ(0)[inodOnElm * nDim + idim] * - e->computeGradient(isol->phi, 0); - double gradM = pbl->fluid->mach->evalGrad(*e, isol->phi, 0) * - e->computeGradient(isol->phi, 0).transpose() * gradV; - double gradrho = pbl->fluid->rho->evalGrad(*e, isol->phi, 0) * - e->computeGradient(isol->phi, 0).transpose() * - gradV; - T_dM_x.push_back(Eigen::Triplet<double>( - jj, e->nodes[inodOnElm]->row * nDim + idim, -gradM / nAdjElms)); - T_drho_x.push_back(Eigen::Triplet<double>( - jj, e->nodes[inodOnElm]->row * nDim + idim, - -gradrho / nAdjElms)); - for (int idim2 = 0; idim2 < nDim; idim2++) - T_dv_x.push_back(Eigen::Triplet<double>( - jj * nDim + idim2, e->nodes[inodOnElm]->row * nDim + idim, - -gradV(idim2) / nAdjElms)); - } - } - ++jj; + dRdStar_x.setFromTriplets(T_dStar_x.begin(), T_dStar_x.end()); + dRdStar_x.prune(0.); + + //**************************************************************************// + // dRblw_x // + //--------------------------------------------------------------------------// + // Analytical // + //**************************************************************************// + std::deque<Eigen::Triplet<double>> T_blw_x; + int offSetElms = 0; + for (const auto &sec : vsol->sections[0]) + { + if (sec->name == "wake") + iReg = 1; + else if (sec->name == "upper" || sec->name == "lower") + iReg = 0; + else + throw std::runtime_error("Unknown section name"); + std::vector<std::vector<double>> gradX = sec->evalGradwrtX(); + std::vector<std::vector<double>> gradY = sec->evalGradwrtY(); + for (size_t i = 0; i < gradX.size(); ++i) + for (size_t j = 0; j < gradX[i].size(); ++j) + { + int rowj = isol->pbl->bBCs[iReg]->nodes[sec->rows[j]]->row; + T_blw_x.push_back(Eigen::Triplet<double>( + i + offSetElms, rowj * isol->pbl->nDim + 0, -gradX[i][j])); + T_blw_x.push_back(Eigen::Triplet<double>( + i + offSetElms, rowj * isol->pbl->nDim + 1, -gradY[i][j])); + } + offSetElms += sec->nElms; + } + dRblw_x.setFromTriplets(T_blw_x.begin(), T_blw_x.end()); + dRblw_x.prune(0.); + + //**************************************************************************// + // dRM_x, dRrho_x, dRv_x // + //--------------------------------------------------------------------------// + // Analytical // + //**************************************************************************// + auto pbl = isol->pbl; + + std::deque<Eigen::Triplet<double>> T_dM_x; + std::deque<Eigen::Triplet<double>> T_drho_x; + std::deque<Eigen::Triplet<double>> T_dv_x; + + int jj = 0; + iReg = 0; + + for (auto sec : vsol->sections[0]) + { + if (sec->name == "wake") + iReg = 1; + else if (sec->name == "upper" || sec->name == "lower") + iReg = 0; + else + throw std::runtime_error( + "gradientswrtInviscidFlow: Unknown section name"); + for (size_t iNode = 0; iNode < sec->nNodes; ++iNode) + { + tbox::Node *srfNode = pbl->bBCs[iReg]->nodes[sec->rows[iNode]]; + double nAdjElms = static_cast<double>(pbl->fluid->neMap[srfNode].size()); + for (auto e : pbl->fluid->neMap[srfNode]) + { + // dv/dx = dv/dgradPhi * dgradPhi/dx = dgradPhi/dx = (-J^(-1) * + // dJ/dx_j)*grad_xPhi dM/dx = dM/dgradPhi * dgradPhi/dx = dM/dgradPhi * + // gradPhi * dv/dx drho/dx = drho/dgradPhi * dgradPhi/dx = drho/dgradPhi + // * gradPhi * dv/dx + for (size_t inodOnElm = 0; inodOnElm < e->nodes.size(); inodOnElm++) + for (int idim = 0; idim < nDim; ++idim) + { + Eigen::VectorXd gradV = (-e->getJinv(0)) * + e->getGradJ(0)[inodOnElm * nDim + idim] * + e->computeGradient(isol->phi, 0); + double gradM = pbl->fluid->mach->evalGrad(*e, isol->phi, 0) * + e->computeGradient(isol->phi, 0).transpose() * gradV; + double gradrho = pbl->fluid->rho->evalGrad(*e, isol->phi, 0) * + e->computeGradient(isol->phi, 0).transpose() * + gradV; + T_dM_x.push_back(Eigen::Triplet<double>( + jj, e->nodes[inodOnElm]->row * nDim + idim, -gradM / nAdjElms)); + T_drho_x.push_back(Eigen::Triplet<double>( + jj, e->nodes[inodOnElm]->row * nDim + idim, + -gradrho / nAdjElms)); + for (int idim2 = 0; idim2 < nDim; idim2++) + T_dv_x.push_back(Eigen::Triplet<double>( + jj * nDim + idim2, e->nodes[inodOnElm]->row * nDim + idim, + -gradV(idim2) / nAdjElms)); + } + } + ++jj; + } } - } - dRM_x.setFromTriplets(T_dM_x.begin(), T_dM_x.end()); - dRrho_x.setFromTriplets(T_drho_x.begin(), T_drho_x.end()); - dRv_x.setFromTriplets(T_dv_x.begin(), T_dv_x.end()); - dRM_x.prune(0.); - dRrho_x.prune(0.); - dRv_x.prune(0.); + dRM_x.setFromTriplets(T_dM_x.begin(), T_dM_x.end()); + dRrho_x.setFromTriplets(T_drho_x.begin(), T_drho_x.end()); + dRv_x.setFromTriplets(T_dv_x.begin(), T_dv_x.end()); + dRM_x.prune(0.); + dRrho_x.prune(0.); + dRv_x.prune(0.); } /** @@ -1093,83 +1146,87 @@ void CoupledAdjoint::gradientwrtMesh() { * * @return Eigen::VectorXd deltaStar */ -Eigen::VectorXd CoupledAdjoint::runViscous() { - int exitCode = vsol->run(); - if (exitCode != 0 && vsol->verbose > 0) - std::cout << "vsol terminated with exit code " << exitCode << std::endl; - - // Extract deltaStar - std::vector<Eigen::VectorXd> dStar_p; - - // Extract parts in a loop - for (size_t i = 0; i < vsol->sections[0].size(); ++i) { - std::vector<double> deltaStar_vec = this->vsol->sections[0][i]->deltaStar; - Eigen::VectorXd deltaStar_eigen = - Eigen::Map<Eigen::VectorXd>(deltaStar_vec.data(), deltaStar_vec.size()); - dStar_p.push_back(deltaStar_eigen); - } - // Calculate the total size - int nNs = 0; - for (const auto &p : dStar_p) - nNs += p.size(); - - // Concatenate the vectors - Eigen::VectorXd deltaStar(nNs); - int idx = 0; - for (const auto &p : dStar_p) { - deltaStar.segment(idx, p.size()) = p; - idx += p.size(); - } - return deltaStar; +Eigen::VectorXd CoupledAdjoint::runViscous() +{ + int exitCode = vsol->run(); + if (exitCode != 0 && vsol->verbose > 0) + std::cout << "vsol terminated with exit code " << exitCode << std::endl; + + // Extract deltaStar + std::vector<Eigen::VectorXd> dStar_p; + + // Extract parts in a loop + for (size_t i = 0; i < vsol->sections[0].size(); ++i) + { + std::vector<double> deltaStar_vec = this->vsol->sections[0][i]->deltaStar; + Eigen::VectorXd deltaStar_eigen = + Eigen::Map<Eigen::VectorXd>(deltaStar_vec.data(), deltaStar_vec.size()); + dStar_p.push_back(deltaStar_eigen); + } + // Calculate the total size + int nNs = 0; + for (const auto &p : dStar_p) + nNs += p.size(); + + // Concatenate the vectors + Eigen::VectorXd deltaStar(nNs); + int idx = 0; + for (const auto &p : dStar_p) + { + deltaStar.segment(idx, p.size()) = p; + idx += p.size(); + } + return deltaStar; } /** * @brief Transpose all matrices of the adjoint problem included in the system. * Not the others. */ -void CoupledAdjoint::transposeMatrices() { - - dRphi_phi = dRphi_phi.transpose(); - dRphi_M = dRphi_M.transpose(); - dRphi_rho = dRphi_rho.transpose(); - dRphi_v = dRphi_v.transpose(); - dRphi_dStar = dRphi_dStar.transpose(); - dRphi_blw = dRphi_blw.transpose(); - - dRM_phi = dRM_phi.transpose(); - dRM_M = dRM_M.transpose(); - dRM_rho = dRM_rho.transpose(); - dRM_v = dRM_v.transpose(); - dRM_dStar = dRM_dStar.transpose(); - dRM_blw = dRM_blw.transpose(); - - dRrho_phi = dRrho_phi.transpose(); - dRrho_M = dRrho_M.transpose(); - dRrho_rho = dRrho_rho.transpose(); - dRrho_v = dRrho_v.transpose(); - dRrho_dStar = dRrho_dStar.transpose(); - dRrho_blw = dRrho_blw.transpose(); - - dRv_phi = dRv_phi.transpose(); - dRv_M = dRv_M.transpose(); - dRv_rho = dRv_rho.transpose(); - dRv_v = dRv_v.transpose(); - dRv_dStar = dRv_dStar.transpose(); - dRv_blw = dRv_blw.transpose(); - - dRdStar_phi = dRdStar_phi.transpose(); - dRdStar_M = dRdStar_M.transpose(); - dRdStar_rho = dRdStar_rho.transpose(); - dRdStar_v = dRdStar_v.transpose(); - dRdStar_dStar = dRdStar_dStar.transpose(); - dRdStar_blw = dRdStar_blw.transpose(); - - dRblw_phi = dRblw_phi.transpose(); - dRblw_M = dRblw_M.transpose(); - dRblw_rho = dRblw_rho.transpose(); - dRblw_v = dRblw_v.transpose(); - dRblw_dStar = dRblw_dStar.transpose(); - dRblw_blw = dRblw_blw.transpose(); +void CoupledAdjoint::transposeMatrices() +{ + + dRphi_phi = dRphi_phi.transpose(); + dRphi_M = dRphi_M.transpose(); + dRphi_rho = dRphi_rho.transpose(); + dRphi_v = dRphi_v.transpose(); + dRphi_dStar = dRphi_dStar.transpose(); + dRphi_blw = dRphi_blw.transpose(); + + dRM_phi = dRM_phi.transpose(); + dRM_M = dRM_M.transpose(); + dRM_rho = dRM_rho.transpose(); + dRM_v = dRM_v.transpose(); + dRM_dStar = dRM_dStar.transpose(); + dRM_blw = dRM_blw.transpose(); + + dRrho_phi = dRrho_phi.transpose(); + dRrho_M = dRrho_M.transpose(); + dRrho_rho = dRrho_rho.transpose(); + dRrho_v = dRrho_v.transpose(); + dRrho_dStar = dRrho_dStar.transpose(); + dRrho_blw = dRrho_blw.transpose(); + + dRv_phi = dRv_phi.transpose(); + dRv_M = dRv_M.transpose(); + dRv_rho = dRv_rho.transpose(); + dRv_v = dRv_v.transpose(); + dRv_dStar = dRv_dStar.transpose(); + dRv_blw = dRv_blw.transpose(); + + dRdStar_phi = dRdStar_phi.transpose(); + dRdStar_M = dRdStar_M.transpose(); + dRdStar_rho = dRdStar_rho.transpose(); + dRdStar_v = dRdStar_v.transpose(); + dRdStar_dStar = dRdStar_dStar.transpose(); + dRdStar_blw = dRdStar_blw.transpose(); + + dRblw_phi = dRblw_phi.transpose(); + dRblw_M = dRblw_M.transpose(); + dRblw_rho = dRblw_rho.transpose(); + dRblw_v = dRblw_v.transpose(); + dRblw_dStar = dRblw_dStar.transpose(); + dRblw_blw = dRblw_blw.transpose(); } /** @@ -1184,29 +1241,36 @@ void CoupledAdjoint::transposeMatrices() { * - The matrix data is written with a precision of 12 digits */ void CoupledAdjoint::writeMatrixToFile(const Eigen::MatrixXd &matrix, - const std::string &filename) { - std::ofstream file(filename); - if (file.is_open()) { - // Write the column headers - for (int j = 0; j < matrix.cols(); ++j) { - file << std::setw(15) << j; - if (j != matrix.cols() - 1) { - file << " "; - } - } - file << "\n"; - - // Write the matrix data - for (int i = 0; i < matrix.rows(); ++i) { - for (int j = 0; j < matrix.cols(); ++j) { - file << std::fixed << std::setprecision(12) << std::setw(15) - << matrix(i, j); - if (j != matrix.cols() - 1) { - file << " "; + const std::string &filename) +{ + std::ofstream file(filename); + if (file.is_open()) + { + // Write the column headers + for (int j = 0; j < matrix.cols(); ++j) + { + file << std::setw(15) << j; + if (j != matrix.cols() - 1) + { + file << " "; + } + } + file << "\n"; + + // Write the matrix data + for (int i = 0; i < matrix.rows(); ++i) + { + for (int j = 0; j < matrix.cols(); ++j) + { + file << std::fixed << std::setprecision(12) << std::setw(15) + << matrix(i, j); + if (j != matrix.cols() - 1) + { + file << " "; + } + } + file << "\n"; } - } - file << "\n"; + file.close(); } - file.close(); - } } diff --git a/blast/src/DCoupledAdjoint.h b/blast/src/DCoupledAdjoint.h index f9a006d8be2cd7f872af773feacbcf24ef049ea8..3087e12d881179ffeb3c50ef6235e7fc99689106 100644 --- a/blast/src/DCoupledAdjoint.h +++ b/blast/src/DCoupledAdjoint.h @@ -26,226 +26,228 @@ #include <iostream> #include <memory> -namespace blast { +namespace blast +{ /** * @brief Adjoint of the coupled system. * @author Paul Dechamps */ -class BLAST_API CoupledAdjoint : public fwk::wSharedObject { +class BLAST_API CoupledAdjoint : public fwk::wSharedObject +{ public: - std::shared_ptr<dart::Adjoint> iadjoint; ///< Adjoint solver - std::shared_ptr<dart::Newton> isol; ///< Inviscid Newton solver - std::shared_ptr<blast::Driver> vsol; ///< Viscous solver - std::shared_ptr<tbox::LinearSolver> alinsol; ///< linear solver (adjoint aero) - - double tdCl_AoA; ///< Total derivative of the lift coefficient with respect to - ///< the angle of attack - double tdCd_AoA; ///< Total erivative of the drag coefficient with respect to - ///< the angle of attack - - std::vector<Eigen::Vector3d> - tdCl_x; ///< Total derivative of the lift coefficient with respect to the - ///< mesh coordinates - std::vector<Eigen::Vector3d> - tdCd_x; ///< Total derivative of the drag coefficient with respect to the - ///< mesh coordinates + std::shared_ptr<dart::Adjoint> iadjoint; ///< Adjoint solver + std::shared_ptr<dart::Newton> isol; ///< Inviscid Newton solver + std::shared_ptr<blast::Driver> vsol; ///< Viscous solver + std::shared_ptr<tbox::LinearSolver> alinsol; ///< linear solver (adjoint aero) + + double tdCl_AoA; ///< Total derivative of the lift coefficient with respect to + ///< the angle of attack + double tdCd_AoA; ///< Total erivative of the drag coefficient with respect to + ///< the angle of attack + + std::vector<Eigen::Vector3d> + tdCl_x; ///< Total derivative of the lift coefficient with respect to the + ///< mesh coordinates + std::vector<Eigen::Vector3d> + tdCd_x; ///< Total derivative of the drag coefficient with respect to the + ///< mesh coordinates private: - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRphi_phi; ///< Inviscid flow Jacobian - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRphi_M; ///< Derivative of the inviscid flow residual with respect to the - ///< Mach number - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRphi_v; ///< Derivative of the inviscid flow residual with respect to the - ///< velocity - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRphi_rho; ///< Derivative of the inviscid flow residual with respect to - ///< the density - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRphi_dStar; ///< Derivative of the inviscid flow residual with respect to - ///< delta* - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRphi_blw; ///< Derivative of the inviscid flow residual with respect to - ///< the blowing velocity - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRphi_x; ///< Derivative of the potential residual with respect to the - ///< mesh coordinates - - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRM_phi; ///< Derivative of the Mach number residual with respect to the - ///< inviscid flow - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRM_M; ///< Derivative of the Mach number residual with respect to the - ///< Mach number - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRM_v; ///< Derivative of the Mach number residual with respect to the - ///< velocity - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRM_rho; ///< Derivative of the Mach number residual with respect to the - ///< density - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRM_dStar; ///< Derivative of the Mach number residual with respect to - ///< delta* - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRM_blw; ///< Derivative of the Mach number residual with respect to the - ///< blowing velocity - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRM_x; ///< Derivative of the Mach number residual with respect to the - ///< mesh coordinates - - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRrho_phi; ///< Derivative of the density residual with respect to the + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRphi_phi; ///< Inviscid flow Jacobian + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRphi_M; ///< Derivative of the inviscid flow residual with respect to the + ///< Mach number + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRphi_v; ///< Derivative of the inviscid flow residual with respect to the + ///< velocity + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRphi_rho; ///< Derivative of the inviscid flow residual with respect to + ///< the density + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRphi_dStar; ///< Derivative of the inviscid flow residual with respect to + ///< delta* + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRphi_blw; ///< Derivative of the inviscid flow residual with respect to + ///< the blowing velocity + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRphi_x; ///< Derivative of the potential residual with respect to the + ///< mesh coordinates + + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRM_phi; ///< Derivative of the Mach number residual with respect to the ///< inviscid flow - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRrho_M; ///< Derivative of the density residual with respect to the Mach - ///< number - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRrho_v; ///< Derivative of the density residual with respect to the + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRM_M; ///< Derivative of the Mach number residual with respect to the + ///< Mach number + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRM_v; ///< Derivative of the Mach number residual with respect to the ///< velocity - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRrho_rho; ///< Derivative of the density residual with respect to the + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRM_rho; ///< Derivative of the Mach number residual with respect to the ///< density - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRrho_dStar; ///< Derivative of the density residual with respect to + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRM_dStar; ///< Derivative of the Mach number residual with respect to ///< delta* - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRrho_blw; ///< Derivative of the density residual with respect to the + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRM_blw; ///< Derivative of the Mach number residual with respect to the ///< blowing velocity - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRrho_x; ///< Derivative of the density residual with respect to the mesh - ///< coordinates + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRM_x; ///< Derivative of the Mach number residual with respect to the + ///< mesh coordinates - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRv_phi; ///< Derivative of the velocity residual with respect to the - ///< inviscid flow - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRv_M; ///< Derivative of the velocity residual with respect to the Mach - ///< number - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRv_v; ///< Derivative of the velocity residual with respect to the - ///< velocity - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRv_rho; ///< Derivative of the velocity residual with respect to the - ///< density - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRv_dStar; ///< Derivative of the velocity residual with respect to delta* - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRv_blw; ///< Derivative of the velocity residual with respect to the - ///< blowing velocity - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRv_x; ///< Derivative of the velocity residual with respect to the mesh - ///< coordinates - - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRdStar_phi; ///< Derivative of delta* residual with respect to the + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRrho_phi; ///< Derivative of the density residual with respect to the ///< inviscid flow - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRdStar_M; ///< Derivative of delta* residual with respect to the Mach + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRrho_M; ///< Derivative of the density residual with respect to the Mach ///< number - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRdStar_v; ///< Derivative of delta* residual with respect to the velocity - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRdStar_rho; ///< Derivative of delta* residual with respect to the + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRrho_v; ///< Derivative of the density residual with respect to the + ///< velocity + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRrho_rho; ///< Derivative of the density residual with respect to the ///< density - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRdStar_dStar; ///< Derivative of delta* residual with respect to delta* - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRdStar_blw; ///< Derivative of delta* residual with respect to the + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRrho_dStar; ///< Derivative of the density residual with respect to + ///< delta* + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRrho_blw; ///< Derivative of the density residual with respect to the ///< blowing velocity - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRdStar_x; ///< Derivative of delta* residual with respect to the mesh + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRrho_x; ///< Derivative of the density residual with respect to the mesh ///< coordinates - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRblw_phi; ///< Derivative of the blowing velocity residual with respect - ///< to the inviscid flow - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRblw_M; ///< Derivative of the blowing velocity residual with respect to - ///< the Mach number - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRblw_v; ///< Derivative of the blowing velocity residual with respect to - ///< the velocity - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRblw_rho; ///< Derivative of the blowing velocity residual with respect - ///< to the density - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRblw_dStar; ///< Derivative of the blowing velocity residual with respect - ///< to delta* - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRblw_blw; ///< Derivative of the blowing velocity residual with respect - ///< to the blowing velocity - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRblw_x; ///< Derivative of the blowing velocity residual with respect to - ///< the mesh coordinates - - // Objective functions - Eigen::VectorXd dCl_phi; ///< Derivative of the lift coefficient with respect - ///< to the inviscid flow - Eigen::VectorXd dCd_phi; ///< Derivative of the drag coefficient with respect - ///< to the inviscid flow - Eigen::VectorXd dCd_M; ///< Derivative of the drag coefficient with respect to - ///< the Mach number - Eigen::VectorXd dCd_rho; ///< Derivative of the drag coefficient with respect - ///< to the density - Eigen::VectorXd dCd_v; ///< Derivative of the drag coefficient with respect to - ///< the velocity - Eigen::VectorXd - dCd_dStar; ///< Derivative of the drag coefficient with respect to delta* - - // Angle of attack - Eigen::VectorXd dRphi_AoA; ///< Derivative of the inviscid flow residual with - ///< respect to the angle of attack - double dCl_AoA; ///< Partial derivative of the lift coefficient with respect - ///< to the angle of attack - double dCd_AoA; ///< Partial derivative of the total drag coefficient with - ///< respect to the angle of attack - - // Mesh - Eigen::SparseMatrix<double, Eigen::RowMajor> - dRx_x; ///< Derivative of the mesh residual with respect to the mesh - ///< coordinates - Eigen::VectorXd dCl_x; ///< Partial derivative of the lift coefficient with - ///< respect to the mesh coordinates - 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 tms2; //< internal precise timers + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRv_phi; ///< Derivative of the velocity residual with respect to the + ///< inviscid flow + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRv_M; ///< Derivative of the velocity residual with respect to the Mach + ///< number + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRv_v; ///< Derivative of the velocity residual with respect to the + ///< velocity + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRv_rho; ///< Derivative of the velocity residual with respect to the + ///< density + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRv_dStar; ///< Derivative of the velocity residual with respect to delta* + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRv_blw; ///< Derivative of the velocity residual with respect to the + ///< blowing velocity + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRv_x; ///< Derivative of the velocity residual with respect to the mesh + ///< coordinates + + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRdStar_phi; ///< Derivative of delta* residual with respect to the + ///< inviscid flow + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRdStar_M; ///< Derivative of delta* residual with respect to the Mach + ///< number + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRdStar_v; ///< Derivative of delta* residual with respect to the velocity + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRdStar_rho; ///< Derivative of delta* residual with respect to the + ///< density + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRdStar_dStar; ///< Derivative of delta* residual with respect to delta* + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRdStar_blw; ///< Derivative of delta* residual with respect to the + ///< blowing velocity + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRdStar_x; ///< Derivative of delta* residual with respect to the mesh + ///< coordinates + + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRblw_phi; ///< Derivative of the blowing velocity residual with respect + ///< to the inviscid flow + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRblw_M; ///< Derivative of the blowing velocity residual with respect to + ///< the Mach number + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRblw_v; ///< Derivative of the blowing velocity residual with respect to + ///< the velocity + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRblw_rho; ///< Derivative of the blowing velocity residual with respect + ///< to the density + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRblw_dStar; ///< Derivative of the blowing velocity residual with respect + ///< to delta* + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRblw_blw; ///< Derivative of the blowing velocity residual with respect + ///< to the blowing velocity + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRblw_x; ///< Derivative of the blowing velocity residual with respect to + ///< the mesh coordinates + + // Objective functions + Eigen::VectorXd dCl_phi; ///< Derivative of the lift coefficient with respect + ///< to the inviscid flow + Eigen::VectorXd dCd_phi; ///< Derivative of the drag coefficient with respect + ///< to the inviscid flow + Eigen::VectorXd dCd_M; ///< Derivative of the drag coefficient with respect to + ///< the Mach number + Eigen::VectorXd dCd_rho; ///< Derivative of the drag coefficient with respect + ///< to the density + Eigen::VectorXd dCd_v; ///< Derivative of the drag coefficient with respect to + ///< the velocity + Eigen::VectorXd + dCd_dStar; ///< Derivative of the drag coefficient with respect to delta* + + // Angle of attack + Eigen::VectorXd dRphi_AoA; ///< Derivative of the inviscid flow residual with + ///< respect to the angle of attack + double dCl_AoA; ///< Partial derivative of the lift coefficient with respect + ///< to the angle of attack + double dCd_AoA; ///< Partial derivative of the total drag coefficient with + ///< respect to the angle of attack + + // Mesh + Eigen::SparseMatrix<double, Eigen::RowMajor> + dRx_x; ///< Derivative of the mesh residual with respect to the mesh + ///< coordinates + Eigen::VectorXd dCl_x; ///< Partial derivative of the lift coefficient with + ///< respect to the mesh coordinates + 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 tms2; //< internal precise timers public: - CoupledAdjoint(std::shared_ptr<dart::Adjoint> iAdjoint, - std::shared_ptr<blast::Driver> vSolver); - virtual ~CoupledAdjoint() { std::cout << "~blast::CoupledAdjoint()\n"; }; + CoupledAdjoint(std::shared_ptr<dart::Adjoint> iAdjoint, + std::shared_ptr<blast::Driver> vSolver); + virtual ~CoupledAdjoint() { std::cout << "~blast::CoupledAdjoint()\n"; }; - void run(); - void reset(); + void run(); + void reset(); private: - void buildAdjointMatrix( - std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor> *>> - &matrices, - Eigen::SparseMatrix<double, Eigen::RowMajor> &matrixAdjoint); - - Eigen::VectorXd runViscous(); - - void gradientswrtInviscidFlow(); - void gradientswrtMachNumber(); - void gradientswrtDensity(); - void gradientswrtVelocity(); - void gradientswrtDeltaStar(); - void gradientswrtBlowingVelocity(); - - void gradientwrtAoA(); - void gradientwrtMesh(); - - void transposeMatrices(); - void writeMatrixToFile(const Eigen::MatrixXd &matrix, - const std::string &filename); + void buildAdjointMatrix( + std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor> *>> + &matrices, + Eigen::SparseMatrix<double, Eigen::RowMajor> &matrixAdjoint); + + Eigen::VectorXd runViscous(); + + void gradientswrtInviscidFlow(); + void gradientswrtMachNumber(); + void gradientswrtDensity(); + void gradientswrtVelocity(); + void gradientswrtDeltaStar(); + void gradientswrtBlowingVelocity(); + + void gradientwrtAoA(); + void gradientwrtMesh(); + + void transposeMatrices(); + void writeMatrixToFile(const Eigen::MatrixXd &matrix, + const std::string &filename); }; } // namespace blast #endif // DDARTADJOINT diff --git a/blast/src/DDriver.cpp b/blast/src/DDriver.cpp index 575682b2ab8b837a4a96295993f853ffdf85e7c4..8a0787c4411c9b480d2e5416ae11eb4d35f147f0 100644 --- a/blast/src/DDriver.cpp +++ b/blast/src/DDriver.cpp @@ -22,54 +22,57 @@ using namespace blast; */ Driver::Driver(double _Re, double _Minf, double _CFL0, size_t nSections, double xtrF_top, double xtrF_bot, double _span, - unsigned int _verbose) { - // Freestream parameters. - Re = _Re; - CFL0 = _CFL0; - verbose = _verbose; - - std::vector<double> xtrF = {xtrF_top, xtrF_bot, 0.}; - - // Initialzing boundary layer - sections.resize(nSections); - solverExitCode = 0; - convergenceStatus.resize(nSections); - for (size_t iSec = 0; iSec < nSections; ++iSec) - convergenceStatus[iSec].resize(3); - - Cdt_sec.resize(sections.size(), 0.); - Cdf_sec.resize(sections.size(), 0.); - Cdp_sec.resize(sections.size(), 0.); - Cdt = 0.; - Cdf = 0.; - Cdp = 0.; - span = _span; - - // Time solver initialization. - tSolver = new Solver(_CFL0, _Minf, Re); - - // Numerical parameters. - std::cout << "--- Viscous solver setup ---\n" - << "Reynolds number: " << Re << " \n" - << "Freestream Mach number: " << _Minf << " \n" - << "Initial CFL number: " << CFL0 << " \n" - << std::endl; - std::cout << "Boundary layer initialized." << std::endl; + unsigned int _verbose) +{ + // Freestream parameters. + Re = _Re; + CFL0 = _CFL0; + verbose = _verbose; + + std::vector<double> xtrF = {xtrF_top, xtrF_bot, 0.}; + + // Initialzing boundary layer + sections.resize(nSections); + solverExitCode = 0; + convergenceStatus.resize(nSections); + for (size_t iSec = 0; iSec < nSections; ++iSec) + convergenceStatus[iSec].resize(3); + + Cdt_sec.resize(sections.size(), 0.); + Cdf_sec.resize(sections.size(), 0.); + Cdp_sec.resize(sections.size(), 0.); + Cdt = 0.; + Cdf = 0.; + Cdp = 0.; + span = _span; + + // Time solver initialization. + tSolver = new Solver(_CFL0, _Minf, Re); + + // Numerical parameters. + std::cout << "--- Viscous solver setup ---\n" + << "Reynolds number: " << Re << " \n" + << "Freestream Mach number: " << _Minf << " \n" + << "Initial CFL number: " << CFL0 << " \n" + << std::endl; + std::cout << "Boundary layer initialized." << std::endl; } /** * @brief Driver and subsequent dependencies destruction. */ -Driver::~Driver() { - delete tSolver; - std::cout << "~blast::Driver()\n"; +Driver::~Driver() +{ + delete tSolver; + std::cout << "~blast::Driver()\n"; } // end ~Driver -void Driver::reset() { - size_t nSections = sections.size(); - for (size_t iSec = 0; iSec < sections.size(); ++iSec) - sections[iSec].resize(0); - sections.resize(nSections); +void Driver::reset() +{ + size_t nSections = sections.size(); + for (size_t iSec = 0; iSec < sections.size(); ++iSec) + sections[iSec].resize(0); + sections.resize(nSections); } /** @@ -79,129 +82,143 @@ void Driver::reset() { * * @return int Solver exit code. */ -int Driver::run() { - bool lockTrans = false; // Flag activation of transition routines. - int pointExitCode; // Output of pseudo time integration (0: converged; -1: - // unsuccessful, failed; >0: unsuccessful nb iter). - - double maxMach = 0.; - - for (size_t iSec = 0; iSec < sections.size(); ++iSec) { - if (verbose > 1) - std::cout << "Sec " << iSec << " "; - - // Loop over the regions (upper, lower and wake). - int iRegion = 0; - for (auto reg : sections[iSec]) { - convergenceStatus[iSec][iRegion].resize(0); - - if (reg->getMaxMach() > maxMach) - maxMach = reg->getMaxMach(); - - // Reset transition - if (reg->name == "wake") - reg->xtr = 0.; - else if (reg->name == "upper" || reg->name == "lower") - reg->xtr = 1.; - else { - std::cout << "reg->name " << reg->name << std::endl; - throw std::runtime_error("Wrong region name\n"); - } - - // Reset regime. - lockTrans = false; - if (reg->name == "upper" || reg->name == "lower") { - for (size_t k = 0; k < reg->nNodes; k++) - reg->regime[k] = 0; - reg->setStagnationBC(Re); - } - - else if (reg->name == "wake") // Wake - { - for (size_t k = 0; k < reg->nNodes; k++) - reg->regime[k] = 1; - - // Impose wake boundary condition. - std::vector<double> upperConditions(reg->getnVar(), 0.); - std::vector<double> lowerConditions(reg->getnVar(), 0.); - for (size_t k = 0; k < upperConditions.size(); ++k) { - upperConditions[k] = - sections[iSec][0]->u[(sections[iSec][0]->nNodes - 1) * - sections[iSec][0]->getnVar() + - k]; - lowerConditions[k] = - sections[iSec][1]->u[(sections[iSec][1]->nNodes - 1) * - sections[iSec][1]->getnVar() + - k]; - } - reg->setWakeBC(Re, upperConditions, lowerConditions); - lockTrans = true; - } else - throw std::runtime_error("Wrong region name\n"); - - // Loop over points - for (size_t iPoint = 1; iPoint < reg->nNodes; ++iPoint) { - // Initialize solution at point - tSolver->initialCondition(iPoint, reg); - // Solve equations - tms["1-Solver"].start(); - pointExitCode = tSolver->integration(iPoint, reg); - tms["1-Solver"].stop(); - - if (pointExitCode != 0) - convergenceStatus[iSec][iRegion].push_back(iPoint); - - // Transition - tms["2-Transition"].start(); - if (!lockTrans) { - // Forced transition - if (reg->xtrF != -1 && reg->xoc[iPoint] >= reg->xtrF) { - reg->u[iPoint * reg->getnVar() + 2] = reg->getnCrit(); - averageTransition(iPoint, reg, true); - if (verbose > 1) { - std::cout << std::fixed; - std::cout << std::setprecision(2); - std::cout << reg->xtr << " (" << reg->transitionNode << ") "; +int Driver::run() +{ + bool lockTrans = false; // Flag activation of transition routines. + int pointExitCode; // Output of pseudo time integration (0: converged; -1: + // unsuccessful, failed; >0: unsuccessful nb iter). + + double maxMach = 0.; + + for (size_t iSec = 0; iSec < sections.size(); ++iSec) + { + if (verbose > 1) + std::cout << "Sec " << iSec << " "; + + // Loop over the regions (upper, lower and wake). + int iRegion = 0; + for (auto reg : sections[iSec]) + { + convergenceStatus[iSec][iRegion].resize(0); + + if (reg->getMaxMach() > maxMach) + maxMach = reg->getMaxMach(); + + // Reset transition + if (reg->name == "wake") + reg->xtr = 0.; + else if (reg->name == "upper" || reg->name == "lower") + reg->xtr = 1.; + else + { + std::cout << "reg->name " << reg->name << std::endl; + throw std::runtime_error("Wrong region name\n"); + } + + // Reset regime. + lockTrans = false; + if (reg->name == "upper" || reg->name == "lower") + { + for (size_t k = 0; k < reg->nNodes; k++) + reg->regime[k] = 0; + reg->setStagnationBC(Re); } - lockTrans = true; - } - // Free transtion - else if (reg->xtrF == -1) { - // Check if perturbation waves are growing or not. - checkWaves(iPoint, reg); - - // Amplification factor is compared to critical amplification - // factor. - if (reg->u[iPoint * reg->getnVar() + 2] >= reg->getnCrit()) { - averageTransition(iPoint, reg, false); - if (verbose > 1) { - std::cout << std::fixed; - std::cout << std::setprecision(2); - std::cout << reg->xtr << " (" << reg->transitionNode << ") "; - } - lockTrans = true; + + else if (reg->name == "wake") // Wake + { + for (size_t k = 0; k < reg->nNodes; k++) + reg->regime[k] = 1; + + // Impose wake boundary condition. + std::vector<double> upperConditions(reg->getnVar(), 0.); + std::vector<double> lowerConditions(reg->getnVar(), 0.); + for (size_t k = 0; k < upperConditions.size(); ++k) + { + upperConditions[k] = + sections[iSec][0]->u[(sections[iSec][0]->nNodes - 1) * + sections[iSec][0]->getnVar() + + k]; + lowerConditions[k] = + sections[iSec][1]->u[(sections[iSec][1]->nNodes - 1) * + sections[iSec][1]->getnVar() + + k]; + } + reg->setWakeBC(Re, upperConditions, lowerConditions); + lockTrans = true; } - } + else + throw std::runtime_error("Wrong region name\n"); + + // Loop over points + for (size_t iPoint = 1; iPoint < reg->nNodes; ++iPoint) + { + // Initialize solution at point + tSolver->initialCondition(iPoint, reg); + // Solve equations + tms["1-Solver"].start(); + pointExitCode = tSolver->integration(iPoint, reg); + tms["1-Solver"].stop(); + + if (pointExitCode != 0) + convergenceStatus[iSec][iRegion].push_back(iPoint); + + // Transition + tms["2-Transition"].start(); + if (!lockTrans) + { + // Forced transition + if (reg->xtrF != -1 && reg->xoc[iPoint] >= reg->xtrF) + { + reg->u[iPoint * reg->getnVar() + 2] = reg->getnCrit(); + averageTransition(iPoint, reg, true); + if (verbose > 1) + { + std::cout << std::fixed; + std::cout << std::setprecision(2); + std::cout << reg->xtr << " (" << reg->transitionNode << ") "; + } + lockTrans = true; + } + // Free transtion + else if (reg->xtrF == -1) + { + // Check if perturbation waves are growing or not. + checkWaves(iPoint, reg); + + // Amplification factor is compared to critical amplification + // factor. + if (reg->u[iPoint * reg->getnVar() + 2] >= reg->getnCrit()) + { + averageTransition(iPoint, reg, false); + if (verbose > 1) + { + std::cout << std::fixed; + std::cout << std::setprecision(2); + std::cout << reg->xtr << " (" << reg->transitionNode << ") "; + } + lockTrans = true; + } + } + } + tms["2-Transition"].stop(); + } + tms["3-Blowing"].start(); + reg->computeBlowingVelocity(); + tms["3-Blowing"].stop(); + iRegion++; } - tms["2-Transition"].stop(); - } - tms["3-Blowing"].start(); - reg->computeBlowingVelocity(); - tms["3-Blowing"].stop(); - iRegion++; + if (verbose > 1) + std::cout << std::endl; } - if (verbose > 1) - std::cout << std::endl; - } - for (size_t i = 0; i < sections.size(); ++i) - computeSectionalDrag(sections[i], i); - computeTotalDrag(); + for (size_t i = 0; i < sections.size(); ++i) + computeSectionalDrag(sections[i], i); + computeTotalDrag(); - // Output information. - solverExitCode = outputStatus(maxMach); + // Output information. + solverExitCode = outputStatus(maxMach); - return solverExitCode; // exit code + return solverExitCode; // exit code } /** @@ -211,20 +228,22 @@ int Driver::run() { * @param iPoint Marker id. * @param bl BoundaryLayer region. */ -void Driver::checkWaves(size_t iPoint, BoundaryLayer *bl) { - - double logRet_crit = - 2.492 * std::pow(1. / (bl->Hk[iPoint] - 1), 0.43) + - 0.7 * (std::tanh(14 * (1. / (bl->Hk[iPoint] - 1)) - 9.24) + 1.0); - - if (bl->Ret[iPoint] > 0) // Avoid log of negative number. - { - if (std::log10(bl->Ret[iPoint]) <= logRet_crit - 0.08) - bl->u[iPoint * bl->getnVar() + 2] = - 0.; // Set N to 0 at that point if waves do not grow. - } else - bl->u[iPoint * bl->getnVar() + 2] = - 0.; // Set N to 0 at that point if waves do not grow. +void Driver::checkWaves(size_t iPoint, BoundaryLayer *bl) +{ + + double logRet_crit = + 2.492 * std::pow(1. / (bl->Hk[iPoint] - 1), 0.43) + + 0.7 * (std::tanh(14 * (1. / (bl->Hk[iPoint] - 1)) - 9.24) + 1.0); + + if (bl->Ret[iPoint] > 0) // Avoid log of negative number. + { + if (std::log10(bl->Ret[iPoint]) <= logRet_crit - 0.08) + bl->u[iPoint * bl->getnVar() + 2] = + 0.; // Set N to 0 at that point if waves do not grow. + } + else + bl->u[iPoint * bl->getnVar() + 2] = + 0.; // Set N to 0 at that point if waves do not grow. } /** @@ -235,114 +254,116 @@ void Driver::checkWaves(size_t iPoint, BoundaryLayer *bl) { * @param bl BoundaryLayer region. * @param forced Flag for forced transition. */ -void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced) { - // Averages solution on transition marker. - size_t nVar = bl->getnVar(); - - // Save laminar solution. - std::vector<double> lamSol(nVar, 0.0); - for (size_t k = 0; k < lamSol.size(); ++k) - lamSol[k] = bl->u[iPoint * nVar + k]; - - // Set up turbulent portion boundary condition. - bl->transitionNode = iPoint; // Save transition marker. - - // Loop over remaining points in the region. - for (size_t k = iPoint; k < bl->nNodes; ++k) - bl->regime[k] = 1; // Set Turbulent regime. - - bl->xtr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * - ((bl->xoc[iPoint] - bl->xoc[iPoint - 1]) / - (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + - bl->xoc[iPoint - 1]; - - double avTurb = - (bl->xoc[iPoint] - bl->xtr) / (bl->xoc[iPoint] - bl->xoc[iPoint - 1]); - double avLam = 1. - avTurb; - - // std::cout << "avTurb " << avTurb << " avLam " << avLam << std::endl; - if (avTurb < 0. || avTurb > 1. || avLam < 0. || avLam > 1.) { - bl->printSolution(iPoint); - bl->printSolution(iPoint - 1); - std::cout << "dx " << std::abs(bl->xoc[iPoint] - bl->xoc[iPoint - 1]) - << std::endl; - std::cout << " (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2]) " - << (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2]) - << std::endl; - std::cout << " bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2] " - << bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2] << std::endl; - std::cout << "xtr " << bl->xtr << std::endl; - std::cout << "avTurb " << avTurb << " avLam " << avLam << std::endl; - throw std::runtime_error("Transition location out of bounds."); - } - - // Impose boundary condition. - double Cteq_trans; - bl->closSolver->turbulentClosures(Cteq_trans, bl->u[(iPoint - 1) * nVar + 0], - bl->u[(iPoint - 1) * nVar + 1], 0.03, - bl->Ret[iPoint - 1], bl->Me[iPoint - 1], - bl->name); - bl->ctEq[iPoint - 1] = Cteq_trans; - bl->u[(iPoint - 1) * nVar + 4] = 0.7 * bl->ctEq[iPoint - 1]; - - // Avoid starting with ill conditioned IC for turbulent computation. (Regime - // was switched above). These initial guess values do not influence the - // converged solution. - bl->u[iPoint * nVar + 4] = - bl->u[(iPoint - 1) * nVar + 4]; // IC of transition point = turbulent BC - // (imposed @ previous point). - bl->u[iPoint * nVar + 1] = - 1.515; // Because we expect a significant drop of the shape factor. - - // Solve point in turbulent configuration. - int exitCode = tSolver->integration(iPoint, bl); - - if (exitCode != 0) - std::cout << "Warning: Transition point turbulent computation terminated " - "with exit code " - << exitCode << std::endl; - - // Average both solutions (Now solution @ iPoint is the turbulent solution). - bl->u[iPoint * nVar + 0] = - avLam * lamSol[0] + avTurb * bl->u[(iPoint)*nVar + 0]; // Theta. - bl->u[iPoint * nVar + 1] = - avLam * lamSol[1] + avTurb * bl->u[(iPoint)*nVar + 1]; // H. - bl->u[iPoint * nVar + 2] = bl->getnCrit(); // N. - bl->u[iPoint * nVar + 3] = - avLam * lamSol[3] + avTurb * bl->u[(iPoint)*nVar + 3]; // ue. - bl->u[iPoint * nVar + 4] = avTurb * bl->u[(iPoint)*nVar + 4]; // Ct. - - bl->deltaStar[iPoint - 1] = - bl->u[(iPoint - 1) * nVar + 0] * bl->u[(iPoint - 1) * nVar + 1]; - bl->deltaStar[iPoint] = bl->u[iPoint * nVar + 0] * bl->u[iPoint * nVar + 1]; - - // Recompute closures. (The turbulent BC @ iPoint - 1 does not influence - // laminar closure relations). - std::vector<double> lamParam(8, 0.); - bl->closSolver->laminarClosures( - lamParam, bl->u[(iPoint - 1) * nVar + 0], bl->u[(iPoint - 1) * nVar + 1], - bl->Ret[iPoint - 1], bl->Me[iPoint - 1], bl->name); - bl->Hstar[iPoint - 1] = lamParam[0]; - bl->Hstar2[iPoint - 1] = lamParam[1]; - bl->Hk[iPoint - 1] = lamParam[2]; - bl->cf[iPoint - 1] = lamParam[3]; - bl->cd[iPoint - 1] = lamParam[4]; - bl->delta[iPoint - 1] = lamParam[5]; - bl->ctEq[iPoint - 1] = lamParam[6]; - bl->us[iPoint - 1] = lamParam[7]; - - std::vector<double> turbParam(8, 0.); - bl->closSolver->turbulentClosures( - turbParam, bl->u[(iPoint)*nVar + 0], bl->u[(iPoint)*nVar + 1], - bl->u[(iPoint)*nVar + 4], bl->Ret[iPoint], bl->Me[iPoint], bl->name); - bl->Hstar[iPoint] = turbParam[0]; - bl->Hstar2[iPoint] = turbParam[1]; - bl->Hk[iPoint] = turbParam[2]; - bl->cf[iPoint] = turbParam[3]; - bl->cd[iPoint] = turbParam[4]; - bl->delta[iPoint] = turbParam[5]; - bl->ctEq[iPoint] = turbParam[6]; - bl->us[iPoint] = turbParam[7]; +void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced) +{ + // Averages solution on transition marker. + size_t nVar = bl->getnVar(); + + // Save laminar solution. + std::vector<double> lamSol(nVar, 0.0); + for (size_t k = 0; k < lamSol.size(); ++k) + lamSol[k] = bl->u[iPoint * nVar + k]; + + // Set up turbulent portion boundary condition. + bl->transitionNode = iPoint; // Save transition marker. + + // Loop over remaining points in the region. + for (size_t k = iPoint; k < bl->nNodes; ++k) + bl->regime[k] = 1; // Set Turbulent regime. + + bl->xtr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * + ((bl->xoc[iPoint] - bl->xoc[iPoint - 1]) / + (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + + bl->xoc[iPoint - 1]; + + double avTurb = + (bl->xoc[iPoint] - bl->xtr) / (bl->xoc[iPoint] - bl->xoc[iPoint - 1]); + double avLam = 1. - avTurb; + + // std::cout << "avTurb " << avTurb << " avLam " << avLam << std::endl; + if (avTurb < 0. || avTurb > 1. || avLam < 0. || avLam > 1.) + { + bl->printSolution(iPoint); + bl->printSolution(iPoint - 1); + std::cout << "dx " << std::abs(bl->xoc[iPoint] - bl->xoc[iPoint - 1]) + << std::endl; + std::cout << " (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2]) " + << (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2]) + << std::endl; + std::cout << " bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2] " + << bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2] << std::endl; + std::cout << "xtr " << bl->xtr << std::endl; + std::cout << "avTurb " << avTurb << " avLam " << avLam << std::endl; + throw std::runtime_error("Transition location out of bounds."); + } + + // Impose boundary condition. + double Cteq_trans; + bl->closSolver->turbulentClosures(Cteq_trans, bl->u[(iPoint - 1) * nVar + 0], + bl->u[(iPoint - 1) * nVar + 1], 0.03, + bl->Ret[iPoint - 1], bl->Me[iPoint - 1], + bl->name); + bl->ctEq[iPoint - 1] = Cteq_trans; + bl->u[(iPoint - 1) * nVar + 4] = 0.7 * bl->ctEq[iPoint - 1]; + + // Avoid starting with ill conditioned IC for turbulent computation. (Regime + // was switched above). These initial guess values do not influence the + // converged solution. + bl->u[iPoint * nVar + 4] = + bl->u[(iPoint - 1) * nVar + 4]; // IC of transition point = turbulent BC + // (imposed @ previous point). + bl->u[iPoint * nVar + 1] = + 1.515; // Because we expect a significant drop of the shape factor. + + // Solve point in turbulent configuration. + int exitCode = tSolver->integration(iPoint, bl); + + if (exitCode != 0) + std::cout << "Warning: Transition point turbulent computation terminated " + "with exit code " + << exitCode << std::endl; + + // Average both solutions (Now solution @ iPoint is the turbulent solution). + bl->u[iPoint * nVar + 0] = + avLam * lamSol[0] + avTurb * bl->u[(iPoint)*nVar + 0]; // Theta. + bl->u[iPoint * nVar + 1] = + avLam * lamSol[1] + avTurb * bl->u[(iPoint)*nVar + 1]; // H. + bl->u[iPoint * nVar + 2] = bl->getnCrit(); // N. + bl->u[iPoint * nVar + 3] = + avLam * lamSol[3] + avTurb * bl->u[(iPoint)*nVar + 3]; // ue. + bl->u[iPoint * nVar + 4] = avTurb * bl->u[(iPoint)*nVar + 4]; // Ct. + + bl->deltaStar[iPoint - 1] = + bl->u[(iPoint - 1) * nVar + 0] * bl->u[(iPoint - 1) * nVar + 1]; + bl->deltaStar[iPoint] = bl->u[iPoint * nVar + 0] * bl->u[iPoint * nVar + 1]; + + // Recompute closures. (The turbulent BC @ iPoint - 1 does not influence + // laminar closure relations). + std::vector<double> lamParam(8, 0.); + bl->closSolver->laminarClosures( + lamParam, bl->u[(iPoint - 1) * nVar + 0], bl->u[(iPoint - 1) * nVar + 1], + bl->Ret[iPoint - 1], bl->Me[iPoint - 1], bl->name); + bl->Hstar[iPoint - 1] = lamParam[0]; + bl->Hstar2[iPoint - 1] = lamParam[1]; + bl->Hk[iPoint - 1] = lamParam[2]; + bl->cf[iPoint - 1] = lamParam[3]; + bl->cd[iPoint - 1] = lamParam[4]; + bl->delta[iPoint - 1] = lamParam[5]; + bl->ctEq[iPoint - 1] = lamParam[6]; + bl->us[iPoint - 1] = lamParam[7]; + + std::vector<double> turbParam(8, 0.); + bl->closSolver->turbulentClosures( + turbParam, bl->u[(iPoint)*nVar + 0], bl->u[(iPoint)*nVar + 1], + bl->u[(iPoint)*nVar + 4], bl->Ret[iPoint], bl->Me[iPoint], bl->name); + bl->Hstar[iPoint] = turbParam[0]; + bl->Hstar2[iPoint] = turbParam[1]; + bl->Hk[iPoint] = turbParam[2]; + bl->cf[iPoint] = turbParam[3]; + bl->cd[iPoint] = turbParam[4]; + bl->delta[iPoint] = turbParam[5]; + bl->ctEq[iPoint] = turbParam[6]; + bl->us[iPoint] = turbParam[7]; } /** @@ -357,55 +378,60 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced) { * @param i Considered section. */ void Driver::computeSectionalDrag(std::vector<BoundaryLayer *> const &sec, - size_t i) { - size_t nVar = sec[0]->getnVar(); - - // Normalize friction and dissipation coefficients. - std::vector<double> frictioncoeff_upper(sec[0]->nNodes, 0.0); - for (size_t k = 0; k < frictioncoeff_upper.size(); ++k) - frictioncoeff_upper[k] = sec[0]->vt[k] * sec[0]->vt[k] * sec[0]->cf[k]; - - std::vector<double> frictioncoeff_lower(sec[1]->nNodes, 0.0); - for (size_t k = 0; k < frictioncoeff_lower.size(); ++k) - frictioncoeff_lower[k] = sec[1]->vt[k] * sec[1]->vt[k] * sec[1]->cf[k]; - - // Compute friction drag coefficient. - double Cdf_upper = 0.0; - for (size_t k = 1; k < sec[0]->nNodes; ++k) - Cdf_upper += (frictioncoeff_upper[k] + frictioncoeff_upper[k - 1]) * - (sec[0]->x[k] - sec[0]->x[k - 1]) / 2; - double Cdf_lower = 0.0; - for (size_t k = 1; k < sec[1]->nNodes; ++k) - Cdf_lower += (frictioncoeff_lower[k] + frictioncoeff_lower[k - 1]) * - (sec[1]->x[k] - sec[1]->x[k - 1]) / 2; - Cdf_sec[i] = Cdf_upper + Cdf_lower; // No friction in the wake - - // Total drag coefficient. Computed from upper and lower TE if there is no - // wake. - if (sec.size() > 2 && sec[2]->getName() == "wake") { - size_t lastWkPt = (sec[2]->nNodes - 1) * nVar; - Cdt_sec[i] = - 2 * sec[2]->u[lastWkPt + 0] * - pow(sec[2]->u[lastWkPt + 3], (sec[2]->u[lastWkPt + 1] + 5) / 2); - } else { - size_t lastUpPt = (sec[0]->nNodes - 1) * nVar; - size_t lastLwPt = (sec[1]->nNodes - 1) * nVar; - Cdt_sec[i] = - sec[0]->u[lastUpPt + 0] * - pow(sec[0]->u[lastUpPt + 3], (sec[0]->u[lastUpPt + 1] + 5) / 2) + - sec[1]->u[lastLwPt + 0] * - pow(sec[1]->u[lastLwPt + 3], (sec[1]->u[lastLwPt + 1] + 5) / 2); - } - // Pressure drag coefficient. - Cdp_sec[i] = Cdt_sec[i] - Cdf_sec[i]; + size_t i) +{ + size_t nVar = sec[0]->getnVar(); + + // Normalize friction and dissipation coefficients. + std::vector<double> frictioncoeff_upper(sec[0]->nNodes, 0.0); + for (size_t k = 0; k < frictioncoeff_upper.size(); ++k) + frictioncoeff_upper[k] = sec[0]->vt[k] * sec[0]->vt[k] * sec[0]->cf[k]; + + std::vector<double> frictioncoeff_lower(sec[1]->nNodes, 0.0); + for (size_t k = 0; k < frictioncoeff_lower.size(); ++k) + frictioncoeff_lower[k] = sec[1]->vt[k] * sec[1]->vt[k] * sec[1]->cf[k]; + + // Compute friction drag coefficient. + double Cdf_upper = 0.0; + for (size_t k = 1; k < sec[0]->nNodes; ++k) + Cdf_upper += (frictioncoeff_upper[k] + frictioncoeff_upper[k - 1]) * + (sec[0]->x[k] - sec[0]->x[k - 1]) / 2; + double Cdf_lower = 0.0; + for (size_t k = 1; k < sec[1]->nNodes; ++k) + Cdf_lower += (frictioncoeff_lower[k] + frictioncoeff_lower[k - 1]) * + (sec[1]->x[k] - sec[1]->x[k - 1]) / 2; + Cdf_sec[i] = Cdf_upper + Cdf_lower; // No friction in the wake + + // Total drag coefficient. Computed from upper and lower TE if there is no + // wake. + if (sec.size() > 2 && sec[2]->getName() == "wake") + { + size_t lastWkPt = (sec[2]->nNodes - 1) * nVar; + Cdt_sec[i] = + 2 * sec[2]->u[lastWkPt + 0] * + pow(sec[2]->u[lastWkPt + 3], (sec[2]->u[lastWkPt + 1] + 5) / 2); + } + else + { + size_t lastUpPt = (sec[0]->nNodes - 1) * nVar; + size_t lastLwPt = (sec[1]->nNodes - 1) * nVar; + Cdt_sec[i] = + sec[0]->u[lastUpPt + 0] * + pow(sec[0]->u[lastUpPt + 3], (sec[0]->u[lastUpPt + 1] + 5) / 2) + + sec[1]->u[lastLwPt + 0] * + pow(sec[1]->u[lastLwPt + 3], (sec[1]->u[lastLwPt + 1] + 5) / 2); + } + // Pressure drag coefficient. + Cdp_sec[i] = Cdt_sec[i] - Cdf_sec[i]; } -double Driver::getAverageTransition(size_t const iRegion) const { - double xtr = 0.0; - for (size_t iSec = 0; iSec < sections.size(); ++iSec) - xtr += sections[iSec][iRegion]->xtr; - xtr /= sections.size(); - return xtr; +double Driver::getAverageTransition(size_t const iRegion) const +{ + double xtr = 0.0; + for (size_t iSec = 0; iSec < sections.size(); ++iSec) + xtr += sections[iSec][iRegion]->xtr; + xtr /= sections.size(); + return xtr; } /** @@ -413,30 +439,35 @@ double Driver::getAverageTransition(size_t const iRegion) const { * Performs the sectional drag integration for 3D computations. * */ -void Driver::computeTotalDrag() { - Cdt = 0.; - Cdf = 0.; - Cdp = 0.; - - if (sections.size() == 1 || span == 0.) { - Cdt = Cdt_sec[0]; - Cdf = Cdf_sec[0]; - Cdp = Cdp_sec[0]; - } else { - Cdt += (sections[0][0]->z[0] - 0) * Cdt_sec[0]; - Cdp += (sections[0][0]->z[0] - 0) * Cdp_sec[0]; - Cdf += (sections[0][0]->z[0] - 0) * Cdf_sec[0]; - double dz = 0.; - for (size_t k = 1; k < sections.size(); ++k) { - dz = (sections[k][0]->z[0] - sections[k - 1][0]->z[0]); - Cdt += dz * Cdt_sec[k]; - Cdp += dz * Cdp_sec[k]; - Cdf += dz * Cdf_sec[k]; +void Driver::computeTotalDrag() +{ + Cdt = 0.; + Cdf = 0.; + Cdp = 0.; + + if (sections.size() == 1 || span == 0.) + { + Cdt = Cdt_sec[0]; + Cdf = Cdf_sec[0]; + Cdp = Cdp_sec[0]; + } + else + { + Cdt += (sections[0][0]->z[0] - 0) * Cdt_sec[0]; + Cdp += (sections[0][0]->z[0] - 0) * Cdp_sec[0]; + Cdf += (sections[0][0]->z[0] - 0) * Cdf_sec[0]; + double dz = 0.; + for (size_t k = 1; k < sections.size(); ++k) + { + dz = (sections[k][0]->z[0] - sections[k - 1][0]->z[0]); + Cdt += dz * Cdt_sec[k]; + Cdp += dz * Cdp_sec[k]; + Cdf += dz * Cdf_sec[k]; + } + Cdt /= span; + Cdp /= span; + Cdf /= span; } - Cdt /= span; - Cdp /= span; - Cdf /= span; - } } /** @@ -445,78 +476,88 @@ void Driver::computeTotalDrag() { * @param maxMach Maximum Mach number identified on the configuration. * @return int Solver exit code (0 converged, !=0 not converged). */ -int Driver::outputStatus(double maxMach) const { - int solverExitCode = 0; - for (size_t iSec = 0; iSec < sections.size(); ++iSec) - for (size_t iReg = 0; iReg < sections[iSec].size(); ++iReg) - if (convergenceStatus[iSec][iReg].size() != 0) { - solverExitCode = -1; - break; - } - - if (verbose > 2 && solverExitCode != 0) { - // Output unconverged points. - std::cout << "Unconverged points" << std::endl; - std::cout << std::setw(10) << "Section" << std::setw(12) << "Region" - << std::setw(16) << "Markers" << std::endl; - for (size_t iSec = 0; iSec < sections.size(); ++iSec) { - unsigned int count = 0; - for (size_t iReg = 0; iReg < sections[iSec].size(); ++iReg) { - std::cout << std::setw(10) << iSec; - std::cout << std::setw(12) << iReg << " "; - if (convergenceStatus[iSec][iReg].size() != 0) - for (size_t iPoint = 0; iPoint < convergenceStatus[iSec][iReg].size(); - ++iPoint) { - std::cout << convergenceStatus[iSec][iReg][iPoint] << " (" - << sections[iSec][iReg] - ->x[convergenceStatus[iSec][iReg][iPoint]] - << "), "; - ++count; - } - std::cout << "\n"; - } - std::cout << "total: " << count << "\n"; - } - std::cout << "" << std::endl; - } - - if (verbose > 0) { - // Compute mean transition locations. - std::vector<double> meanTransitions(2, 0.); +int Driver::outputStatus(double maxMach) const +{ + int solverExitCode = 0; for (size_t iSec = 0; iSec < sections.size(); ++iSec) - for (size_t iReg = 0; iReg < 2; ++iReg) - meanTransitions[iReg] += sections[iSec][iReg]->xtr; - for (size_t iReg = 0; iReg < meanTransitions.size(); ++iReg) - meanTransitions[iReg] /= sections.size(); - - std::cout << std::fixed << std::setprecision(2) << "Viscous solver for Re " - << Re / 1e6 << "e6, Mach max " << maxMach << std::endl; - - std::cout << std::setw(10) << "xTr Upper" << std::setw(12) << "xTr Lower" - << std::setw(12) << "Cd" << std::endl; - std::cout << std::fixed << std::setprecision(2); - std::cout << std::setw(10) << meanTransitions[0] << std::setw(12) - << meanTransitions[1]; - std::cout << std::fixed << std::setprecision(5); - std::cout << std::setw(12) << Cdt << "\n"; - - if (verbose > 0) { - if (solverExitCode == 0) - std::cout << ANSI_COLOR_GREEN << "Viscous solver converged." - << ANSI_COLOR_RESET << std::endl; - else - std::cout << ANSI_COLOR_YELLOW << "Viscous solver not converged." - << ANSI_COLOR_RESET << std::endl; + for (size_t iReg = 0; iReg < sections[iSec].size(); ++iReg) + if (convergenceStatus[iSec][iReg].size() != 0) + { + solverExitCode = -1; + break; + } + + if (verbose > 2 && solverExitCode != 0) + { + // Output unconverged points. + std::cout << "Unconverged points" << std::endl; + std::cout << std::setw(10) << "Section" << std::setw(12) << "Region" + << std::setw(16) << "Markers" << std::endl; + for (size_t iSec = 0; iSec < sections.size(); ++iSec) + { + unsigned int count = 0; + for (size_t iReg = 0; iReg < sections[iSec].size(); ++iReg) + { + std::cout << std::setw(10) << iSec; + std::cout << std::setw(12) << iReg << " "; + if (convergenceStatus[iSec][iReg].size() != 0) + for (size_t iPoint = 0; iPoint < convergenceStatus[iSec][iReg].size(); + ++iPoint) + { + std::cout << convergenceStatus[iSec][iReg][iPoint] << " (" + << sections[iSec][iReg] + ->x[convergenceStatus[iSec][iReg][iPoint]] + << "), "; + ++count; + } + std::cout << "\n"; + } + std::cout << "total: " << count << "\n"; + } + std::cout << "" << std::endl; } - if (verbose > 2) { - std::cout << "Viscous solver CPU" << std::endl << tms; + if (verbose > 0) + { + // Compute mean transition locations. + std::vector<double> meanTransitions(2, 0.); + for (size_t iSec = 0; iSec < sections.size(); ++iSec) + for (size_t iReg = 0; iReg < 2; ++iReg) + meanTransitions[iReg] += sections[iSec][iReg]->xtr; + for (size_t iReg = 0; iReg < meanTransitions.size(); ++iReg) + meanTransitions[iReg] /= sections.size(); + + std::cout << std::fixed << std::setprecision(2) << "Viscous solver for Re " + << Re / 1e6 << "e6, Mach max " << maxMach << std::endl; + + std::cout << std::setw(10) << "xTr Upper" << std::setw(12) << "xTr Lower" + << std::setw(12) << "Cd" << std::endl; + std::cout << std::fixed << std::setprecision(2); + std::cout << std::setw(10) << meanTransitions[0] << std::setw(12) + << meanTransitions[1]; + std::cout << std::fixed << std::setprecision(5); + std::cout << std::setw(12) << Cdt << "\n"; + + if (verbose > 0) + { + if (solverExitCode == 0) + std::cout << ANSI_COLOR_GREEN << "Viscous solver converged." + << ANSI_COLOR_RESET << std::endl; + else + std::cout << ANSI_COLOR_YELLOW << "Viscous solver not converged." + << ANSI_COLOR_RESET << std::endl; + } + + if (verbose > 2) + { + std::cout << "Viscous solver CPU" << std::endl + << tms; + } } - } - if (verbose > 0) - std::cout << "" << std::endl; - return solverExitCode; + if (verbose > 0) + std::cout << "" << std::endl; + return solverExitCode; } /** @@ -526,102 +567,106 @@ int Driver::outputStatus(double maxMach) const { * @param iSec id of the considered section. * @return std::vector<std::vector<double>> */ -std::vector<std::vector<double>> Driver::getSolution(size_t iSec) { - size_t nMarkersUp = sections[iSec][0]->nNodes; - size_t nMarkersLw = sections[iSec][1]->nNodes; // No stagnation point doublon. - size_t nVar = sections[iSec][0]->getnVar(); - - std::vector<std::vector<double>> Sln(20); - for (size_t iVector = 0; iVector < Sln.size(); ++iVector) - Sln[iVector].resize(nMarkersUp + nMarkersLw + sections[iSec][2]->nNodes, - 0.); - // Upper side. - size_t revPoint = nMarkersUp - 1; - for (size_t iPoint = 0; iPoint < nMarkersUp; ++iPoint) { - Sln[0][iPoint] = sections[iSec][0]->u[revPoint * nVar + 0]; - Sln[1][iPoint] = sections[iSec][0]->u[revPoint * nVar + 1]; - Sln[2][iPoint] = sections[iSec][0]->u[revPoint * nVar + 2]; - Sln[3][iPoint] = sections[iSec][0]->u[revPoint * nVar + 3]; - Sln[4][iPoint] = sections[iSec][0]->u[revPoint * nVar + 4]; - Sln[5][iPoint] = sections[iSec][0]->deltaStar[revPoint]; - - Sln[6][iPoint] = sections[iSec][0]->Ret[revPoint]; - Sln[7][iPoint] = sections[iSec][0]->cd[revPoint]; - Sln[8][iPoint] = sections[iSec][0]->cf[revPoint]; - Sln[9][iPoint] = sections[iSec][0]->Hstar[revPoint]; - Sln[10][iPoint] = sections[iSec][0]->Hstar2[revPoint]; - Sln[11][iPoint] = sections[iSec][0]->Hk[revPoint]; - Sln[12][iPoint] = sections[iSec][0]->ctEq[revPoint]; - Sln[13][iPoint] = sections[iSec][0]->us[revPoint]; - Sln[14][iPoint] = sections[iSec][0]->delta[revPoint]; - - Sln[15][iPoint] = sections[iSec][0]->x[revPoint]; - Sln[16][iPoint] = sections[iSec][0]->xoc[revPoint]; - Sln[17][iPoint] = sections[iSec][0]->y[revPoint]; - Sln[18][iPoint] = sections[iSec][0]->z[revPoint]; - Sln[19][iPoint] = sections[iSec][0]->vt[revPoint]; - --revPoint; - } - // Lower side. - for (size_t iPoint = 0; iPoint < sections[iSec][1]->nNodes; ++iPoint) { - Sln[0][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 0]; - Sln[1][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 1]; - Sln[2][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 2]; - Sln[3][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 3]; - Sln[4][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 4]; - Sln[5][nMarkersUp + iPoint] = sections[iSec][1]->deltaStar[iPoint]; - - Sln[6][nMarkersUp + iPoint] = sections[iSec][1]->Ret[iPoint]; - Sln[7][nMarkersUp + iPoint] = sections[iSec][1]->cd[iPoint]; - Sln[8][nMarkersUp + iPoint] = sections[iSec][1]->cf[iPoint]; - Sln[9][nMarkersUp + iPoint] = sections[iSec][1]->Hstar[iPoint]; - Sln[10][nMarkersUp + iPoint] = sections[iSec][1]->Hstar2[iPoint]; - Sln[11][nMarkersUp + iPoint] = sections[iSec][1]->Hk[iPoint]; - Sln[12][nMarkersUp + iPoint] = sections[iSec][1]->ctEq[iPoint]; - Sln[13][nMarkersUp + iPoint] = sections[iSec][1]->us[iPoint]; - Sln[14][nMarkersUp + iPoint] = sections[iSec][1]->delta[iPoint]; - - Sln[15][nMarkersUp + iPoint] = sections[iSec][1]->x[iPoint]; - Sln[16][nMarkersUp + iPoint] = sections[iSec][1]->xoc[iPoint]; - Sln[17][nMarkersUp + iPoint] = sections[iSec][1]->y[iPoint]; - Sln[18][nMarkersUp + iPoint] = sections[iSec][1]->z[iPoint]; - Sln[19][nMarkersUp + iPoint] = sections[iSec][1]->vt[iPoint]; - } - // Wake. - for (size_t iPoint = 0; iPoint < sections[iSec][2]->nNodes; ++iPoint) { - Sln[0][nMarkersUp + nMarkersLw + iPoint] = - sections[iSec][2]->u[iPoint * nVar + 0]; - Sln[1][nMarkersUp + nMarkersLw + iPoint] = - sections[iSec][2]->u[iPoint * nVar + 1]; - Sln[2][nMarkersUp + nMarkersLw + iPoint] = - sections[iSec][2]->u[iPoint * nVar + 2]; - Sln[3][nMarkersUp + nMarkersLw + iPoint] = - sections[iSec][2]->u[iPoint * nVar + 3]; - Sln[4][nMarkersUp + nMarkersLw + iPoint] = - sections[iSec][2]->u[iPoint * nVar + 4]; - Sln[5][nMarkersUp + nMarkersLw + iPoint] = - sections[iSec][2]->deltaStar[iPoint]; - - Sln[6][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->Ret[iPoint]; - Sln[7][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->cd[iPoint]; - Sln[8][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->cf[iPoint]; - Sln[9][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->Hstar[iPoint]; - Sln[10][nMarkersUp + nMarkersLw + iPoint] = - sections[iSec][2]->Hstar2[iPoint]; - Sln[11][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->Hk[iPoint]; - Sln[12][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->ctEq[iPoint]; - Sln[13][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->us[iPoint]; - Sln[14][nMarkersUp + nMarkersLw + iPoint] = - sections[iSec][2]->delta[iPoint]; - - Sln[15][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->x[iPoint]; - Sln[16][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->xoc[iPoint]; - Sln[17][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->y[iPoint]; - Sln[18][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->z[iPoint]; - Sln[19][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->vt[iPoint]; - } - - if (verbose > 2) - std::cout << "Solution structure sent." << std::endl; - return Sln; +std::vector<std::vector<double>> Driver::getSolution(size_t iSec) +{ + size_t nMarkersUp = sections[iSec][0]->nNodes; + size_t nMarkersLw = sections[iSec][1]->nNodes; // No stagnation point doublon. + size_t nVar = sections[iSec][0]->getnVar(); + + std::vector<std::vector<double>> Sln(20); + for (size_t iVector = 0; iVector < Sln.size(); ++iVector) + Sln[iVector].resize(nMarkersUp + nMarkersLw + sections[iSec][2]->nNodes, + 0.); + // Upper side. + size_t revPoint = nMarkersUp - 1; + for (size_t iPoint = 0; iPoint < nMarkersUp; ++iPoint) + { + Sln[0][iPoint] = sections[iSec][0]->u[revPoint * nVar + 0]; + Sln[1][iPoint] = sections[iSec][0]->u[revPoint * nVar + 1]; + Sln[2][iPoint] = sections[iSec][0]->u[revPoint * nVar + 2]; + Sln[3][iPoint] = sections[iSec][0]->u[revPoint * nVar + 3]; + Sln[4][iPoint] = sections[iSec][0]->u[revPoint * nVar + 4]; + Sln[5][iPoint] = sections[iSec][0]->deltaStar[revPoint]; + + Sln[6][iPoint] = sections[iSec][0]->Ret[revPoint]; + Sln[7][iPoint] = sections[iSec][0]->cd[revPoint]; + Sln[8][iPoint] = sections[iSec][0]->cf[revPoint]; + Sln[9][iPoint] = sections[iSec][0]->Hstar[revPoint]; + Sln[10][iPoint] = sections[iSec][0]->Hstar2[revPoint]; + Sln[11][iPoint] = sections[iSec][0]->Hk[revPoint]; + Sln[12][iPoint] = sections[iSec][0]->ctEq[revPoint]; + Sln[13][iPoint] = sections[iSec][0]->us[revPoint]; + Sln[14][iPoint] = sections[iSec][0]->delta[revPoint]; + + Sln[15][iPoint] = sections[iSec][0]->x[revPoint]; + Sln[16][iPoint] = sections[iSec][0]->xoc[revPoint]; + Sln[17][iPoint] = sections[iSec][0]->y[revPoint]; + Sln[18][iPoint] = sections[iSec][0]->z[revPoint]; + Sln[19][iPoint] = sections[iSec][0]->vt[revPoint]; + --revPoint; + } + // Lower side. + for (size_t iPoint = 0; iPoint < sections[iSec][1]->nNodes; ++iPoint) + { + Sln[0][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 0]; + Sln[1][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 1]; + Sln[2][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 2]; + Sln[3][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 3]; + Sln[4][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 4]; + Sln[5][nMarkersUp + iPoint] = sections[iSec][1]->deltaStar[iPoint]; + + Sln[6][nMarkersUp + iPoint] = sections[iSec][1]->Ret[iPoint]; + Sln[7][nMarkersUp + iPoint] = sections[iSec][1]->cd[iPoint]; + Sln[8][nMarkersUp + iPoint] = sections[iSec][1]->cf[iPoint]; + Sln[9][nMarkersUp + iPoint] = sections[iSec][1]->Hstar[iPoint]; + Sln[10][nMarkersUp + iPoint] = sections[iSec][1]->Hstar2[iPoint]; + Sln[11][nMarkersUp + iPoint] = sections[iSec][1]->Hk[iPoint]; + Sln[12][nMarkersUp + iPoint] = sections[iSec][1]->ctEq[iPoint]; + Sln[13][nMarkersUp + iPoint] = sections[iSec][1]->us[iPoint]; + Sln[14][nMarkersUp + iPoint] = sections[iSec][1]->delta[iPoint]; + + Sln[15][nMarkersUp + iPoint] = sections[iSec][1]->x[iPoint]; + Sln[16][nMarkersUp + iPoint] = sections[iSec][1]->xoc[iPoint]; + Sln[17][nMarkersUp + iPoint] = sections[iSec][1]->y[iPoint]; + Sln[18][nMarkersUp + iPoint] = sections[iSec][1]->z[iPoint]; + Sln[19][nMarkersUp + iPoint] = sections[iSec][1]->vt[iPoint]; + } + // Wake. + for (size_t iPoint = 0; iPoint < sections[iSec][2]->nNodes; ++iPoint) + { + Sln[0][nMarkersUp + nMarkersLw + iPoint] = + sections[iSec][2]->u[iPoint * nVar + 0]; + Sln[1][nMarkersUp + nMarkersLw + iPoint] = + sections[iSec][2]->u[iPoint * nVar + 1]; + Sln[2][nMarkersUp + nMarkersLw + iPoint] = + sections[iSec][2]->u[iPoint * nVar + 2]; + Sln[3][nMarkersUp + nMarkersLw + iPoint] = + sections[iSec][2]->u[iPoint * nVar + 3]; + Sln[4][nMarkersUp + nMarkersLw + iPoint] = + sections[iSec][2]->u[iPoint * nVar + 4]; + Sln[5][nMarkersUp + nMarkersLw + iPoint] = + sections[iSec][2]->deltaStar[iPoint]; + + Sln[6][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->Ret[iPoint]; + Sln[7][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->cd[iPoint]; + Sln[8][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->cf[iPoint]; + Sln[9][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->Hstar[iPoint]; + Sln[10][nMarkersUp + nMarkersLw + iPoint] = + sections[iSec][2]->Hstar2[iPoint]; + Sln[11][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->Hk[iPoint]; + Sln[12][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->ctEq[iPoint]; + Sln[13][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->us[iPoint]; + Sln[14][nMarkersUp + nMarkersLw + iPoint] = + sections[iSec][2]->delta[iPoint]; + + Sln[15][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->x[iPoint]; + Sln[16][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->xoc[iPoint]; + Sln[17][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->y[iPoint]; + Sln[18][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->z[iPoint]; + Sln[19][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->vt[iPoint]; + } + + if (verbose > 2) + std::cout << "Solution structure sent." << std::endl; + return Sln; } \ No newline at end of file diff --git a/blast/src/DDriver.h b/blast/src/DDriver.h index 1ff4ba258f00701696a17b595ce0340638a7c41c..07287d4587d0e8cdf6454d68ad09a5b40ae10a7d 100644 --- a/blast/src/DDriver.h +++ b/blast/src/DDriver.h @@ -5,65 +5,68 @@ #include "blast.h" #include "wObject.h" #include "wTimers.h" -namespace blast { +namespace blast +{ /** * @brief Boundary layer solver driver. Performs the space marching and set up * time marching for each point. * @authors Paul Dechamps */ -class BLAST_API Driver : public fwk::wSharedObject { +class BLAST_API Driver : public fwk::wSharedObject +{ public: - double Cdt; ///< Total drag coefficient. - double Cdf; ///< Total friction drag coefficient. - double Cdp; ///< Total pressure drag coefficient. - std::vector<double> Cdt_sec; ///< Sectional total drag coefficient. - std::vector<double> Cdf_sec; ///< Sectional friction drag coefficient. - std::vector<double> Cdp_sec; ///< Sectional pressure drag coefficient. - std::vector<std::vector<BoundaryLayer *>> - sections; ///< Boundary layer regions sorted by section. - unsigned int verbose; ///< Verbosity level of the solver 0<=verbose<=1. + double Cdt; ///< Total drag coefficient. + double Cdf; ///< Total friction drag coefficient. + double Cdp; ///< Total pressure drag coefficient. + std::vector<double> Cdt_sec; ///< Sectional total drag coefficient. + std::vector<double> Cdf_sec; ///< Sectional friction drag coefficient. + std::vector<double> Cdp_sec; ///< Sectional pressure drag coefficient. + std::vector<std::vector<BoundaryLayer *>> + sections; ///< Boundary layer regions sorted by section. + unsigned int verbose; ///< Verbosity level of the solver 0<=verbose<=1. private: - double Re; ///< Freestream Reynolds number. - double CFL0; ///< Initial CFL number for pseudo time integration. - double span; ///< Wing Span (Used only for drag computation, not used if 2D - ///< case). - Solver *tSolver; ///< Pseudo-time solver. - int solverExitCode; ///< Exit code of viscous calculations. - std::vector<std::vector<std::vector<size_t>>> - convergenceStatus; ///< Vector containing points that did not converged. + double Re; ///< Freestream Reynolds number. + double CFL0; ///< Initial CFL number for pseudo time integration. + double span; ///< Wing Span (Used only for drag computation, not used if 2D + ///< case). + Solver *tSolver; ///< Pseudo-time solver. + int solverExitCode; ///< Exit code of viscous calculations. + std::vector<std::vector<std::vector<size_t>>> + convergenceStatus; ///< Vector containing points that did not converged. protected: - fwk::Timers tms; ///< Internal timers + fwk::Timers tms; ///< Internal timers public: - Driver(double _Re, double _Minf, double _CFL0, size_t nSections, - double xtrF_top = -1., double xtrF_bot = -1., double _span = 0., - unsigned int _verbose = 1); - ~Driver(); - int run(); - void reset(); + Driver(double _Re, double _Minf, double _CFL0, size_t nSections, + double xtrF_top = -1., double xtrF_bot = -1., double _span = 0., + unsigned int _verbose = 1); + ~Driver(); + int run(); + void reset(); - // Getters. - double getAverageTransition(size_t const iRegion) const; - double getRe() const { return Re; }; - std::vector<std::vector<double>> getSolution(size_t iSec); + // Getters. + double getAverageTransition(size_t const iRegion) const; + double getRe() const { return Re; }; + std::vector<std::vector<double>> getSolution(size_t iSec); - // Setters. - void setVerbose(unsigned int _verbose) { verbose = _verbose; }; + // Setters. + void setVerbose(unsigned int _verbose) { verbose = _verbose; }; - void addSection(size_t iSec, BoundaryLayer *&bl) { - sections[iSec].push_back(bl); - }; + void addSection(size_t iSec, BoundaryLayer *&bl) + { + sections[iSec].push_back(bl); + }; private: - void checkWaves(size_t iPoint, BoundaryLayer *bl); - void averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced); - void computeSectionalDrag(std::vector<BoundaryLayer *> const &bl, size_t i); - void computeTotalDrag(); - void computeBlwVel(); - int outputStatus(double maxMach) const; + void checkWaves(size_t iPoint, BoundaryLayer *bl); + void averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced); + void computeSectionalDrag(std::vector<BoundaryLayer *> const &bl, size_t i); + void computeTotalDrag(); + void computeBlwVel(); + int outputStatus(double maxMach) const; }; } // namespace blast #endif // DDRIVER_H diff --git a/blast/src/DFluxes.cpp b/blast/src/DFluxes.cpp index cb9e69d5942fb63c4b57c18f48b9e7099d741f2f..73e89db8542c20ecbe4376ed4ab103ee3cf88b78 100644 --- a/blast/src/DFluxes.cpp +++ b/blast/src/DFluxes.cpp @@ -12,12 +12,13 @@ using namespace Eigen; * @param bl Boundary layer region. * @return VectorXd */ -VectorXd Fluxes::computeResiduals(size_t iPoint, BoundaryLayer *bl) { - size_t nVar = bl->getnVar(); - std::vector<double> u(nVar, 0.); - for (size_t k = 0; k < nVar; ++k) - u[k] = bl->u[iPoint * nVar + k]; - return blLaws(iPoint, bl, u); +VectorXd Fluxes::computeResiduals(size_t iPoint, BoundaryLayer *bl) +{ + size_t nVar = bl->getnVar(); + std::vector<double> u(nVar, 0.); + for (size_t k = 0; k < nVar; ++k) + u[k] = bl->u[iPoint * nVar + k]; + return blLaws(iPoint, bl, u); } /** @@ -30,21 +31,23 @@ VectorXd Fluxes::computeResiduals(size_t iPoint, BoundaryLayer *bl) { * @return MatrixXd */ MatrixXd Fluxes::computeJacobian(size_t iPoint, BoundaryLayer *bl, - VectorXd const &sysRes, double eps) { - size_t nVar = bl->getnVar(); - MatrixXd JacMatrix = MatrixXd::Zero(5, 5); - std::vector<double> uUp(nVar, 0.); - for (size_t k = 0; k < nVar; ++k) - uUp[k] = bl->u[iPoint * nVar + k]; + VectorXd const &sysRes, double eps) +{ + size_t nVar = bl->getnVar(); + MatrixXd JacMatrix = MatrixXd::Zero(5, 5); + std::vector<double> uUp(nVar, 0.); + for (size_t k = 0; k < nVar; ++k) + uUp[k] = bl->u[iPoint * nVar + k]; - double varSave = 0.; - for (size_t iVar = 0; iVar < nVar; ++iVar) { - varSave = uUp[iVar]; - uUp[iVar] += eps; - JacMatrix.col(iVar) = (blLaws(iPoint, bl, uUp) - sysRes) / eps; - uUp[iVar] = varSave; - } - return JacMatrix; + double varSave = 0.; + for (size_t iVar = 0; iVar < nVar; ++iVar) + { + varSave = uUp[iVar]; + uUp[iVar] += eps; + JacMatrix.col(iVar) = (blLaws(iPoint, bl, uUp) - sysRes) / eps; + uUp[iVar] = varSave; + } + return JacMatrix; } /** @@ -56,229 +59,239 @@ MatrixXd Fluxes::computeJacobian(size_t iPoint, BoundaryLayer *bl, * @return VectorXd */ VectorXd Fluxes::blLaws(size_t iPoint, BoundaryLayer *bl, - std::vector<double> u) { - size_t nVar = bl->getnVar(); + std::vector<double> u) +{ + size_t nVar = bl->getnVar(); - // Dissipation ratio. - double dissipRatio = 0.; - if (bl->name == "wake") - dissipRatio = 0.9; - else - dissipRatio = 1.; + // Dissipation ratio. + double dissipRatio = 0.; + if (bl->name == "wake") + dissipRatio = 0.9; + else + dissipRatio = 1.; - bl->Ret[iPoint] = - std::max(u[3] * u[0] * Re, 1.0); // Reynolds number based on theta. - bl->deltaStar[iPoint] = u[0] * u[1]; // Displacement thickness. + bl->Ret[iPoint] = + std::max(u[3] * u[0] * Re, 1.0); // Reynolds number based on theta. + bl->deltaStar[iPoint] = u[0] * u[1]; // Displacement thickness. - // Boundary layer closure relations. - if (bl->regime[iPoint] == 0) { - // Laminar closure relations. - std::vector<double> lamParam(8, 0.); - bl->closSolver->laminarClosures(lamParam, u[0], u[1], bl->Ret[iPoint], - bl->Me[iPoint], bl->name); - bl->Hstar[iPoint] = lamParam[0]; - bl->Hstar2[iPoint] = lamParam[1]; - bl->Hk[iPoint] = lamParam[2]; - bl->cf[iPoint] = lamParam[3]; - bl->cd[iPoint] = lamParam[4]; - bl->delta[iPoint] = lamParam[5]; - bl->ctEq[iPoint] = lamParam[6]; - bl->us[iPoint] = lamParam[7]; - } + // Boundary layer closure relations. + if (bl->regime[iPoint] == 0) + { + // Laminar closure relations. + std::vector<double> lamParam(8, 0.); + bl->closSolver->laminarClosures(lamParam, u[0], u[1], bl->Ret[iPoint], + bl->Me[iPoint], bl->name); + bl->Hstar[iPoint] = lamParam[0]; + bl->Hstar2[iPoint] = lamParam[1]; + bl->Hk[iPoint] = lamParam[2]; + bl->cf[iPoint] = lamParam[3]; + bl->cd[iPoint] = lamParam[4]; + bl->delta[iPoint] = lamParam[5]; + bl->ctEq[iPoint] = lamParam[6]; + bl->us[iPoint] = lamParam[7]; + } - else if (bl->regime[iPoint] == 1) { - // Turbulent closure relations. - std::vector<double> turbParam(8, 0.); - bl->closSolver->turbulentClosures( - turbParam, u[0], u[1], u[4], bl->Ret[iPoint], bl->Me[iPoint], bl->name); - bl->Hstar[iPoint] = turbParam[0]; - bl->Hstar2[iPoint] = turbParam[1]; - bl->Hk[iPoint] = turbParam[2]; - bl->cf[iPoint] = turbParam[3]; - bl->cd[iPoint] = turbParam[4]; - bl->delta[iPoint] = turbParam[5]; - bl->ctEq[iPoint] = turbParam[6]; - bl->us[iPoint] = turbParam[7]; - } else { - std::cout << "Wrong regime\n" << std::endl; - throw; - } + else if (bl->regime[iPoint] == 1) + { + // Turbulent closure relations. + std::vector<double> turbParam(8, 0.); + bl->closSolver->turbulentClosures( + turbParam, u[0], u[1], u[4], bl->Ret[iPoint], bl->Me[iPoint], bl->name); + bl->Hstar[iPoint] = turbParam[0]; + bl->Hstar2[iPoint] = turbParam[1]; + bl->Hk[iPoint] = turbParam[2]; + bl->cf[iPoint] = turbParam[3]; + bl->cd[iPoint] = turbParam[4]; + bl->delta[iPoint] = turbParam[5]; + bl->ctEq[iPoint] = turbParam[6]; + bl->us[iPoint] = turbParam[7]; + } + else + { + std::cout << "Wrong regime\n" + << std::endl; + throw; + } - // Gradients computation. - double dx = (bl->loc[iPoint] - bl->loc[iPoint - 1]); - double dxExt = (bl->xxExt[iPoint] - bl->xxExt[iPoint - 1]); - double dt_dx = (u[0] - bl->u[(iPoint - 1) * nVar + 0]) / dx; - double dH_dx = (u[1] - bl->u[(iPoint - 1) * nVar + 1]) / dx; - double due_dx = (u[3] - bl->u[(iPoint - 1) * nVar + 3]) / dx; - double dct_dx = (u[4] - bl->u[(iPoint - 1) * nVar + 4]) / dx; - double dHstar_dx = (bl->Hstar[iPoint] - bl->Hstar[iPoint - 1]) / dx; + // Gradients computation. + double dx = (bl->loc[iPoint] - bl->loc[iPoint - 1]); + double dxExt = (bl->xxExt[iPoint] - bl->xxExt[iPoint - 1]); + double dt_dx = (u[0] - bl->u[(iPoint - 1) * nVar + 0]) / dx; + double dH_dx = (u[1] - bl->u[(iPoint - 1) * nVar + 1]) / dx; + double due_dx = (u[3] - bl->u[(iPoint - 1) * nVar + 3]) / dx; + double dct_dx = (u[4] - bl->u[(iPoint - 1) * nVar + 4]) / dx; + double dHstar_dx = (bl->Hstar[iPoint] - bl->Hstar[iPoint - 1]) / dx; - double dueExt_dx = (bl->vt[iPoint] - bl->vt[iPoint - 1]) / dxExt; - double ddeltaStarExt_dx = - (bl->deltaStarExt[iPoint] - bl->deltaStarExt[iPoint - 1]) / dxExt; + double dueExt_dx = (bl->vt[iPoint] - bl->vt[iPoint - 1]) / dxExt; + double ddeltaStarExt_dx = + (bl->deltaStarExt[iPoint] - bl->deltaStarExt[iPoint - 1]) / dxExt; - double c = 4. / (M_PI * dx); - double cExt = 4. / (M_PI * dxExt); + double c = 4. / (M_PI * dx); + double cExt = 4. / (M_PI * dxExt); - // Amplification routine. - double dN_dx = 0.0; - double ax = 0.0; - if (bl->xtrF == -1.0 && bl->regime[iPoint] == 0) { - // No forced transition and laminar regime. - ax = amplificationRoutine(bl->Hk[iPoint], bl->Ret[iPoint], u[0]); - dN_dx = (u[2] - bl->u[(iPoint - 1) * nVar + 2]) / dx; - } + // Amplification routine. + double dN_dx = 0.0; + double ax = 0.0; + if (bl->xtrF == -1.0 && bl->regime[iPoint] == 0) + { + // No forced transition and laminar regime. + ax = amplificationRoutine(bl->Hk[iPoint], bl->Ret[iPoint], u[0]); + dN_dx = (u[2] - bl->u[(iPoint - 1) * nVar + 2]) / dx; + } - double Mea = bl->Me[iPoint]; - double Hstara = bl->Hstar[iPoint]; - double Hstar2a = bl->Hstar2[iPoint]; - double Hka = bl->Hk[iPoint]; - double cfa = bl->cf[iPoint]; - double cda = bl->cd[iPoint]; - double deltaa = bl->delta[iPoint]; - double ctEqa = bl->ctEq[iPoint]; - double usa = bl->us[iPoint]; + double Mea = bl->Me[iPoint]; + double Hstara = bl->Hstar[iPoint]; + double Hstar2a = bl->Hstar2[iPoint]; + double Hka = bl->Hk[iPoint]; + double cfa = bl->cf[iPoint]; + double cda = bl->cd[iPoint]; + double deltaa = bl->delta[iPoint]; + double ctEqa = bl->ctEq[iPoint]; + double usa = bl->us[iPoint]; - //--------------------------------------------------------------------------------// - // Integral boundary layer equations. // + //--------------------------------------------------------------------------------// + // 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; + // 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; - // 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) + // { + // 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; + //--------------------------------------------------------------------------------// - VectorXd result = VectorXd::Zero(5); + VectorXd result = VectorXd::Zero(5); - if (bl->regime[iPoint] == 0) { - 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 * + if (bl->regime[iPoint] == 0) + { + 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)) * - (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 * + (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)) * - (-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])); - } + (-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])); + } - return result; + return result; } /** @@ -289,56 +302,59 @@ VectorXd Fluxes::blLaws(size_t iPoint, BoundaryLayer *bl, * @param theta Momentum thickness. * @return double */ -double Fluxes::amplificationRoutine(double Hk, double Ret, double theta) const { - double dgr = 0.08; - double Hk1 = 3.5; - double Hk2 = 4.; - double Hmi = (1 / (Hk - 1)); - double logRet_crit = - 2.492 * std::pow(1. / (Hk - 1.), 0.43) + - 0.7 * (std::tanh(14. * Hmi - 9.24) + - 1.0); // value at which the amplification starts to grow +double Fluxes::amplificationRoutine(double Hk, double Ret, double theta) const +{ + double dgr = 0.08; + double Hk1 = 3.5; + double Hk2 = 4.; + double Hmi = (1 / (Hk - 1)); + double logRet_crit = + 2.492 * std::pow(1. / (Hk - 1.), 0.43) + + 0.7 * (std::tanh(14. * Hmi - 9.24) + + 1.0); // value at which the amplification starts to grow - double ax = 0.; - if (Ret <= 0.) - Ret = 1.; - if (std::log10(Ret) < logRet_crit - dgr) - ax = 0.; - else { - double rNorm = (std::log10(Ret) - (logRet_crit - dgr)) / (2. * dgr); - double Rfac = 0.; - if (rNorm <= 0) - Rfac = 0.; - if (rNorm >= 1) - Rfac = 1.; + double ax = 0.; + if (Ret <= 0.) + Ret = 1.; + if (std::log10(Ret) < logRet_crit - dgr) + ax = 0.; else - Rfac = 3. * rNorm * rNorm - 2. * rNorm * rNorm * rNorm; + { + double rNorm = (std::log10(Ret) - (logRet_crit - dgr)) / (2. * dgr); + double Rfac = 0.; + if (rNorm <= 0) + Rfac = 0.; + if (rNorm >= 1) + Rfac = 1.; + else + Rfac = 3. * rNorm * rNorm - 2. * rNorm * rNorm * rNorm; - // Normal envelope amplification rate - double m = -0.05 + 2.7 * Hmi - 5.5 * Hmi * Hmi + 3. * Hmi * Hmi * Hmi + - 0.1 * std::exp(-20. * Hmi); - double arg = 3.87 * Hmi - 2.52; - double dndRet = 0.028 * (Hk - 1.) - 0.0345 * std::exp(-(arg * arg)); - ax = (m * dndRet / theta) * Rfac; + // Normal envelope amplification rate + double m = -0.05 + 2.7 * Hmi - 5.5 * Hmi * Hmi + 3. * Hmi * Hmi * Hmi + + 0.1 * std::exp(-20. * Hmi); + double arg = 3.87 * Hmi - 2.52; + double dndRet = 0.028 * (Hk - 1.) - 0.0345 * std::exp(-(arg * arg)); + ax = (m * dndRet / theta) * Rfac; - // Set correction for separated profiles - if (Hk > Hk1) { - double Hnorm = (Hk - Hk1) / (Hk2 - Hk1); - double Hfac = 0.; - if (Hnorm >= 1.) - Hfac = 1.; - else - Hfac = 3. * Hnorm * Hnorm - 2. * Hnorm * Hnorm * Hnorm; - double ax1 = ax; - double gro = 0.3 + 0.35 * std::exp(-0.15 * (Hk - 5.)); - double tnr = std::tanh(1.2 * (std::log10(Ret) - gro)); - double ax2 = (0.086 * tnr - 0.25 / std::pow((Hk - 1.), 1.5)) / theta; - if (ax2 < 0.) - ax2 = 0.; - ax = Hfac * ax2 + (1. - Hfac) * ax1; - if (ax < 0.) - ax = 0.; + // Set correction for separated profiles + if (Hk > Hk1) + { + double Hnorm = (Hk - Hk1) / (Hk2 - Hk1); + double Hfac = 0.; + if (Hnorm >= 1.) + Hfac = 1.; + else + Hfac = 3. * Hnorm * Hnorm - 2. * Hnorm * Hnorm * Hnorm; + double ax1 = ax; + double gro = 0.3 + 0.35 * std::exp(-0.15 * (Hk - 5.)); + double tnr = std::tanh(1.2 * (std::log10(Ret) - gro)); + double ax2 = (0.086 * tnr - 0.25 / std::pow((Hk - 1.), 1.5)) / theta; + if (ax2 < 0.) + ax2 = 0.; + ax = Hfac * ax2 + (1. - Hfac) * ax1; + if (ax < 0.) + ax = 0.; + } } - } - return ax; + return ax; } \ No newline at end of file diff --git a/blast/src/DFluxes.h b/blast/src/DFluxes.h index c3f1fad643e08e5a728efb94a929780fe2eb0f05..ef2f15f157bd7b024cb5138b595a9c736bf2fa57 100644 --- a/blast/src/DFluxes.h +++ b/blast/src/DFluxes.h @@ -5,28 +5,30 @@ #include "blast.h" #include <Eigen/Dense> -namespace blast { +namespace blast +{ /** * @brief Residual and Jacobian computation of the unsteady integral boundary * layer equations. * @author Paul Dechamps */ -class BLAST_API Fluxes { +class BLAST_API Fluxes +{ public: - double Re; ///< Freestream Reynolds number. + double Re; ///< Freestream Reynolds number. public: - Fluxes(double _Re) : Re(_Re){}; - ~Fluxes(){}; - Eigen::VectorXd computeResiduals(size_t iPoint, BoundaryLayer *bl); - Eigen::MatrixXd computeJacobian(size_t iPoint, BoundaryLayer *bl, - Eigen::VectorXd const &sysRes, double eps); + Fluxes(double _Re) : Re(_Re){}; + ~Fluxes(){}; + Eigen::VectorXd computeResiduals(size_t iPoint, BoundaryLayer *bl); + Eigen::MatrixXd computeJacobian(size_t iPoint, BoundaryLayer *bl, + Eigen::VectorXd const &sysRes, double eps); private: - Eigen::VectorXd blLaws(size_t iPoint, BoundaryLayer *bl, - std::vector<double> u); - double amplificationRoutine(double Hk, double Ret, double theta) const; + Eigen::VectorXd blLaws(size_t iPoint, BoundaryLayer *bl, + std::vector<double> u); + double amplificationRoutine(double Hk, double Ret, double theta) const; }; } // namespace blast #endif // DFLUXES_H diff --git a/blast/src/DSolver.cpp b/blast/src/DSolver.cpp index 38f46c110ae7acc26979c01e91cf0d2f7a34e644..75fe55c4d55e2fe8de01120c7a1aa2cfe79ad2ce 100644 --- a/blast/src/DSolver.cpp +++ b/blast/src/DSolver.cpp @@ -20,22 +20,24 @@ using namespace blast; * frozen. */ Solver::Solver(double _CFL0, double _Minf, double Re, unsigned int _maxIt, - double _tol, unsigned int _itFrozenJac) { - CFL0 = _CFL0; - maxIt = _maxIt; - tol = _tol; - itFrozenJac0 = _itFrozenJac; - Minf = std::max(_Minf, 0.1); - sysMatrix = new Fluxes(Re); + double _tol, unsigned int _itFrozenJac) +{ + CFL0 = _CFL0; + maxIt = _maxIt; + tol = _tol; + itFrozenJac0 = _itFrozenJac; + Minf = std::max(_Minf, 0.1); + sysMatrix = new Fluxes(Re); } /** * @brief Destroy the Time Solver:: Time Solver object. * */ -Solver::~Solver() { - delete sysMatrix; - std::cout << "~blast::Solver()\n"; +Solver::~Solver() +{ + delete sysMatrix; + std::cout << "~blast::Solver()\n"; } /** @@ -44,24 +46,25 @@ Solver::~Solver() { * @param iPoint Marker id. * @param bl Boundary layer region. */ -void Solver::initialCondition(size_t iPoint, BoundaryLayer *bl) { - size_t nVar = bl->getnVar(); - for (size_t k = 0; k < nVar; ++k) - bl->u[iPoint * nVar + k] = bl->u[(iPoint - 1) * nVar + k]; - bl->u[iPoint * nVar + 3] = bl->vt[iPoint]; - - bl->Ret[iPoint] = bl->Ret[iPoint - 1]; - bl->cd[iPoint] = bl->cd[iPoint - 1]; - bl->cf[iPoint] = bl->cf[iPoint - 1]; - bl->Hstar[iPoint] = bl->Hstar[iPoint - 1]; - bl->Hstar2[iPoint] = bl->Hstar2[iPoint - 1]; - bl->Hk[iPoint] = bl->Hk[iPoint - 1]; - bl->ctEq[iPoint] = bl->ctEq[iPoint - 1]; - bl->us[iPoint] = bl->us[iPoint - 1]; - bl->delta[iPoint] = bl->delta[iPoint - 1]; - bl->deltaStar[iPoint] = bl->deltaStar[iPoint - 1]; - if (bl->regime[iPoint] == 1 && bl->u[iPoint * nVar + 4] <= 0.) - bl->u[iPoint * nVar + 4] = 0.03; +void Solver::initialCondition(size_t iPoint, BoundaryLayer *bl) +{ + size_t nVar = bl->getnVar(); + for (size_t k = 0; k < nVar; ++k) + bl->u[iPoint * nVar + k] = bl->u[(iPoint - 1) * nVar + k]; + bl->u[iPoint * nVar + 3] = bl->vt[iPoint]; + + bl->Ret[iPoint] = bl->Ret[iPoint - 1]; + bl->cd[iPoint] = bl->cd[iPoint - 1]; + bl->cf[iPoint] = bl->cf[iPoint - 1]; + bl->Hstar[iPoint] = bl->Hstar[iPoint - 1]; + bl->Hstar2[iPoint] = bl->Hstar2[iPoint - 1]; + bl->Hk[iPoint] = bl->Hk[iPoint - 1]; + bl->ctEq[iPoint] = bl->ctEq[iPoint - 1]; + bl->us[iPoint] = bl->us[iPoint - 1]; + bl->delta[iPoint] = bl->delta[iPoint - 1]; + bl->deltaStar[iPoint] = bl->deltaStar[iPoint - 1]; + if (bl->regime[iPoint] == 1 && bl->u[iPoint * nVar + 4] <= 0.) + bl->u[iPoint * nVar + 4] = 0.03; } /** @@ -71,67 +74,70 @@ void Solver::initialCondition(size_t iPoint, BoundaryLayer *bl) { * @param bl Boundary layer region. * @return int */ -int Solver::integration(size_t iPoint, BoundaryLayer *bl) { - size_t nVar = bl->getnVar(); - - // Save initial condition. - std::vector<double> sln0(nVar, 0.); - for (size_t i = 0; i < nVar; ++i) - sln0[i] = bl->u[iPoint * nVar + i]; - - // Initialize time integration variables. - double dx = bl->loc[iPoint] - bl->loc[iPoint - 1]; - - // Initial time step. - double CFL = CFL0; - unsigned int itFrozenJac = itFrozenJac0; - double timeStep = setTimeStep(CFL, dx, bl->vt[iPoint]); - MatrixXd IMatrix(5, 5); - IMatrix = MatrixXd::Identity(5, 5) / timeStep; +int Solver::integration(size_t iPoint, BoundaryLayer *bl) +{ + size_t nVar = bl->getnVar(); + + // Save initial condition. + std::vector<double> sln0(nVar, 0.); + for (size_t i = 0; i < nVar; ++i) + sln0[i] = bl->u[iPoint * nVar + i]; + + // Initialize time integration variables. + double dx = bl->loc[iPoint] - bl->loc[iPoint - 1]; + + // Initial time step. + double CFL = CFL0; + unsigned int itFrozenJac = itFrozenJac0; + double timeStep = setTimeStep(CFL, dx, bl->vt[iPoint]); + MatrixXd IMatrix(5, 5); + IMatrix = MatrixXd::Identity(5, 5) / timeStep; - // Initial system residuals. - VectorXd sysRes = sysMatrix->computeResiduals(iPoint, bl); - double normSysRes0 = sysRes.norm(); - double normSysRes = normSysRes0; + // Initial system residuals. + VectorXd sysRes = sysMatrix->computeResiduals(iPoint, bl); + double normSysRes0 = sysRes.norm(); + double normSysRes = normSysRes0; - MatrixXd JacMatrix = MatrixXd::Zero(5, 5); - VectorXd slnIncr = VectorXd::Zero(5); + MatrixXd JacMatrix = MatrixXd::Zero(5, 5); + VectorXd slnIncr = VectorXd::Zero(5); - unsigned int innerIters = 0; // Inner (non-linear) iterations + unsigned int innerIters = 0; // Inner (non-linear) iterations - while (innerIters < maxIt) { - // Jacobian computation. - if (innerIters % itFrozenJac == 0) - JacMatrix = sysMatrix->computeJacobian(iPoint, bl, sysRes, 1e-8); + while (innerIters < maxIt) + { + // Jacobian computation. + if (innerIters % itFrozenJac == 0) + JacMatrix = sysMatrix->computeJacobian(iPoint, bl, sysRes, 1e-8); - // Linear system. - slnIncr = (JacMatrix + IMatrix).partialPivLu().solve(-sysRes); + // Linear system. + slnIncr = (JacMatrix + IMatrix).partialPivLu().solve(-sysRes); - // Solution increment. - for (size_t k = 0; k < nVar; ++k) - bl->u[iPoint * nVar + k] += slnIncr(k); - bl->u[iPoint * nVar + 0] = std::max(bl->u[iPoint * nVar + 0], 1e-6); - bl->u[iPoint * nVar + 1] = std::max(bl->u[iPoint * nVar + 1], 1.0005); + // Solution increment. + for (size_t k = 0; k < nVar; ++k) + bl->u[iPoint * nVar + k] += slnIncr(k); + bl->u[iPoint * nVar + 0] = std::max(bl->u[iPoint * nVar + 0], 1e-6); + bl->u[iPoint * nVar + 1] = std::max(bl->u[iPoint * nVar + 1], 1.0005); - sysRes = sysMatrix->computeResiduals(iPoint, bl); - normSysRes = (sysRes + slnIncr / timeStep).norm(); + sysRes = sysMatrix->computeResiduals(iPoint, bl); + normSysRes = (sysRes + slnIncr / timeStep).norm(); - if (normSysRes <= tol * normSysRes0) - return 0; // Successfull exit. + if (normSysRes <= tol * normSysRes0) + return 0; // Successfull exit. - // CFL adaptation. - CFL = std::max(CFL0 * pow(normSysRes0 / normSysRes, 0.7), 0.1); - timeStep = setTimeStep(CFL, dx, bl->vt[iPoint]); - IMatrix = MatrixXd::Identity(5, 5) / timeStep; + // CFL adaptation. + CFL = std::max(CFL0 * pow(normSysRes0 / normSysRes, 0.7), 0.1); + timeStep = setTimeStep(CFL, dx, bl->vt[iPoint]); + IMatrix = MatrixXd::Identity(5, 5) / timeStep; - innerIters++; - } + innerIters++; + } - if (std::isnan(normSysRes) || normSysRes / normSysRes0 > 1e-3) { - resetSolution(iPoint, bl, sln0); - return -1; - } - return innerIters; + if (std::isnan(normSysRes) || normSysRes / normSysRes0 > 1e-3) + { + resetSolution(iPoint, bl, sln0); + return -1; + } + return innerIters; } /** @@ -142,16 +148,17 @@ int Solver::integration(size_t iPoint, BoundaryLayer *bl) { * @param invVel Inviscid velocity. * @return double */ -double Solver::setTimeStep(double CFL, double dx, double invVel) const { - // Speed of sound. - double gradU2 = (invVel * invVel); - double soundSpeed = computeSoundSpeed(gradU2); +double Solver::setTimeStep(double CFL, double dx, double invVel) const +{ + // Speed of sound. + double gradU2 = (invVel * invVel); + double soundSpeed = computeSoundSpeed(gradU2); - // Velocity of the fastest traveling wave. - double advVel = soundSpeed + invVel; + // Velocity of the fastest traveling wave. + double advVel = soundSpeed + invVel; - // Time step computation. - return CFL * dx / advVel; + // Time step computation. + return CFL * dx / advVel; } /** @@ -160,8 +167,9 @@ double Solver::setTimeStep(double CFL, double dx, double invVel) const { * @param gradU2 Inviscid velocity squared in the cell. * @return double */ -double Solver::computeSoundSpeed(double gradU2) const { - return sqrt(1. / (Minf * Minf) + 0.5 * (gamma - 1.) * (1. - gradU2)); +double Solver::computeSoundSpeed(double gradU2) const +{ + return sqrt(1. / (Minf * Minf) + 0.5 * (gamma - 1.) * (1. - gradU2)); } /** @@ -174,56 +182,61 @@ double Solver::computeSoundSpeed(double gradU2) const { * @param usePrevPoint if 1, use the solution at previous point instead of sln0. */ void Solver::resetSolution(size_t iPoint, BoundaryLayer *bl, - std::vector<double> sln0, bool usePrevPoint) { - size_t nVar = bl->getnVar(); - - // Reset solution. - if (usePrevPoint) - for (size_t k = 0; k < nVar; ++k) - bl->u[iPoint * nVar + k] = bl->u[(iPoint - 1) * nVar + k]; - else - for (size_t k = 0; k < nVar; ++k) - bl->u[iPoint * nVar + k] = sln0[k]; - - // Reset closures. - bl->Ret[iPoint] = std::max(bl->u[iPoint * nVar + 3] * - bl->u[iPoint * nVar + 0] * sysMatrix->Re, - 1.0); // Reynolds number based on theta. - bl->deltaStar[iPoint] = bl->u[iPoint * nVar + 0] * - bl->u[iPoint * nVar + 1]; // Displacement thickness. - - if (bl->regime[iPoint] == 0) { - std::vector<double> lamParam(8, 0.); - bl->closSolver->laminarClosures(lamParam, bl->u[iPoint * nVar + 0], - bl->u[iPoint * nVar + 1], bl->Ret[iPoint], - bl->Me[iPoint], bl->name); - bl->Hstar[iPoint] = lamParam[0]; - bl->Hstar2[iPoint] = lamParam[1]; - bl->Hk[iPoint] = lamParam[2]; - bl->cf[iPoint] = lamParam[3]; - bl->cd[iPoint] = lamParam[4]; - bl->delta[iPoint] = lamParam[5]; - bl->ctEq[iPoint] = lamParam[6]; - bl->us[iPoint] = lamParam[7]; - } else if (bl->regime[iPoint] == 1) { - std::vector<double> turbParam(8, 0.); - bl->closSolver->turbulentClosures( - turbParam, bl->u[iPoint * nVar + 0], bl->u[iPoint * nVar + 1], - bl->u[iPoint * nVar + 4], bl->Ret[iPoint], bl->Me[iPoint], bl->name); - bl->Hstar[iPoint] = turbParam[0]; - bl->Hstar2[iPoint] = turbParam[1]; - bl->Hk[iPoint] = turbParam[2]; - bl->cf[iPoint] = turbParam[3]; - bl->cd[iPoint] = turbParam[4]; - bl->delta[iPoint] = turbParam[5]; - bl->ctEq[iPoint] = turbParam[6]; - bl->us[iPoint] = turbParam[7]; - } + std::vector<double> sln0, bool usePrevPoint) +{ + size_t nVar = bl->getnVar(); + + // Reset solution. + if (usePrevPoint) + for (size_t k = 0; k < nVar; ++k) + bl->u[iPoint * nVar + k] = bl->u[(iPoint - 1) * nVar + k]; + else + for (size_t k = 0; k < nVar; ++k) + bl->u[iPoint * nVar + k] = sln0[k]; + + // Reset closures. + bl->Ret[iPoint] = std::max(bl->u[iPoint * nVar + 3] * + bl->u[iPoint * nVar + 0] * sysMatrix->Re, + 1.0); // Reynolds number based on theta. + bl->deltaStar[iPoint] = bl->u[iPoint * nVar + 0] * + bl->u[iPoint * nVar + 1]; // Displacement thickness. + + if (bl->regime[iPoint] == 0) + { + std::vector<double> lamParam(8, 0.); + bl->closSolver->laminarClosures(lamParam, bl->u[iPoint * nVar + 0], + bl->u[iPoint * nVar + 1], bl->Ret[iPoint], + bl->Me[iPoint], bl->name); + bl->Hstar[iPoint] = lamParam[0]; + bl->Hstar2[iPoint] = lamParam[1]; + bl->Hk[iPoint] = lamParam[2]; + bl->cf[iPoint] = lamParam[3]; + bl->cd[iPoint] = lamParam[4]; + bl->delta[iPoint] = lamParam[5]; + bl->ctEq[iPoint] = lamParam[6]; + bl->us[iPoint] = lamParam[7]; + } + else if (bl->regime[iPoint] == 1) + { + std::vector<double> turbParam(8, 0.); + bl->closSolver->turbulentClosures( + turbParam, bl->u[iPoint * nVar + 0], bl->u[iPoint * nVar + 1], + bl->u[iPoint * nVar + 4], bl->Ret[iPoint], bl->Me[iPoint], bl->name); + bl->Hstar[iPoint] = turbParam[0]; + bl->Hstar2[iPoint] = turbParam[1]; + bl->Hk[iPoint] = turbParam[2]; + bl->cf[iPoint] = turbParam[3]; + bl->cd[iPoint] = turbParam[4]; + bl->delta[iPoint] = turbParam[5]; + bl->ctEq[iPoint] = turbParam[6]; + bl->us[iPoint] = turbParam[7]; + } } -void Solver::setCFL0(double _CFL0) { - if (_CFL0 > 0) - CFL0 = _CFL0; - else - throw std::runtime_error("Negative CFL numbers are not allowed.\n"); +void Solver::setCFL0(double _CFL0) +{ + if (_CFL0 > 0) + CFL0 = _CFL0; + else + throw std::runtime_error("Negative CFL numbers are not allowed.\n"); } \ No newline at end of file diff --git a/blast/src/DSolver.h b/blast/src/DSolver.h index 604c0061858ee5bcf97404ab4e8f47218ee2781b..89c4b4175fb1ddff24e7e7886dc9c20dd132ae63 100644 --- a/blast/src/DSolver.h +++ b/blast/src/DSolver.h @@ -4,49 +4,52 @@ #include "DFluxes.h" #include "blast.h" -namespace blast { +namespace blast +{ /** * @brief Pseudo-time integration. * @author Paul Dechamps */ -class BLAST_API Solver { +class BLAST_API Solver +{ private: - double CFL0; ///< Initial CFL number. - unsigned int maxIt; ///< Maximum number of iterations. - double tol; ///< Convergence tolerance. - unsigned int itFrozenJac0; ///< Number of iterations between which the - ///< Jacobian is frozen. - double const gamma = 1.4; ///< Air heat capacity ratio. - double Minf; ///< Freestream Mach number. + double CFL0; ///< Initial CFL number. + unsigned int maxIt; ///< Maximum number of iterations. + double tol; ///< Convergence tolerance. + unsigned int itFrozenJac0; ///< Number of iterations between which the + ///< Jacobian is frozen. + double const gamma = 1.4; ///< Air heat capacity ratio. + double Minf; ///< Freestream Mach number. - Fluxes *sysMatrix; + Fluxes *sysMatrix; public: - Solver(double _CFL0, double _Minf, double Re, unsigned int _maxIt = 100, - double _tol = 1e-8, unsigned int _itFrozenJac = 1); - ~Solver(); - void initialCondition(size_t iPoint, BoundaryLayer *bl); - int integration(size_t iPoint, BoundaryLayer *bl); - - // Getters. - double getCFL0() const { return CFL0; }; - unsigned int getmaxIt() const { return maxIt; }; - unsigned int getitFrozenJac() const { return itFrozenJac0; }; - - // Setters. - void setCFL0(double _CFL0); - void setmaxIt(unsigned int _maxIt) { maxIt = _maxIt; }; - void setitFrozenJac(unsigned int _itFrozenJac) { - itFrozenJac0 = _itFrozenJac; - }; + Solver(double _CFL0, double _Minf, double Re, unsigned int _maxIt = 100, + double _tol = 1e-8, unsigned int _itFrozenJac = 1); + ~Solver(); + void initialCondition(size_t iPoint, BoundaryLayer *bl); + int integration(size_t iPoint, BoundaryLayer *bl); + + // Getters. + double getCFL0() const { return CFL0; }; + unsigned int getmaxIt() const { return maxIt; }; + unsigned int getitFrozenJac() const { return itFrozenJac0; }; + + // Setters. + void setCFL0(double _CFL0); + void setmaxIt(unsigned int _maxIt) { maxIt = _maxIt; }; + void setitFrozenJac(unsigned int _itFrozenJac) + { + itFrozenJac0 = _itFrozenJac; + }; private: - double setTimeStep(double CFL, double dx, double invVel) const; - double computeSoundSpeed(double gradU2) const; - void resetSolution(size_t iPoint, BoundaryLayer *bl, std::vector<double> Sln0, - bool usePrevPoint = false); + double setTimeStep(double CFL, double dx, double invVel) const; + double computeSoundSpeed(double gradU2) const; + void resetSolution(size_t iPoint, BoundaryLayer *bl, std::vector<double> Sln0, + bool usePrevPoint = false); }; } // namespace blast #endif // DSOLVER_H diff --git a/blast/src/blast.h b/blast/src/blast.h index 9b526a72d2ff83cca967adc14a6346f2e42452ae..0fcf7951412b4df50b17eec5806697d7cd61474e 100644 --- a/blast/src/blast.h +++ b/blast/src/blast.h @@ -34,7 +34,8 @@ /** * @brief Namespace for blast module */ -namespace blast { +namespace blast +{ // Solvers class Driver;