From edccc1776122bfa5eb62e79f7eab66037609db44 Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Mon, 7 Oct 2024 10:46:08 +0200
Subject: [PATCH] (feat) Analytical gradients for adjoint
Change potential and mesh derivatives to analytical. Also changed the residuals computation of the viscous solver to the analytical expression instead of solving the linear system. This commit reduces the computational cost of the coupled adjoint methodology by a factor 4, at least.
---
blast/src/DCoupledAdjoint.cpp | 605 ++++++++++++++++------------------
blast/src/DCoupledAdjoint.h | 13 +-
blast/src/DFluxes.cpp | 74 +++--
3 files changed, 338 insertions(+), 354 deletions(-)
diff --git a/blast/src/DCoupledAdjoint.cpp b/blast/src/DCoupledAdjoint.cpp
index 8ea9318..4cb1449 100644
--- a/blast/src/DCoupledAdjoint.cpp
+++ b/blast/src/DCoupledAdjoint.cpp
@@ -151,7 +151,7 @@ void CoupledAdjoint::reset()
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::vector<Eigen::Triplet<double>> triplets;
+ std::deque<Eigen::Triplet<double>> triplets;
// Sanity check
for (size_t irow = 0; irow < matrices.size(); ++irow)
@@ -192,36 +192,37 @@ void CoupledAdjoint::buildAdjointMatrix(std::vector<std::vector<Eigen::SparseMat
void CoupledAdjoint::run()
{
- tms["0-adj_total"].start();
- tms["1-adj_inviscid"].start();
+ tms["0-Total"].start();
+ tms["1-Derivatives"].start();
+ tms2["1-adj_inviscid"].start();
gradientswrtInviscidFlow();
- tms["1-adj_inviscid"].stop();
- tms["2-adj_mach"].start();
+ tms2["1-adj_inviscid"].stop();
+ tms2["2-adj_mach"].start();
gradientswrtMachNumber();
- tms["2-adj_mach"].stop();
- tms["3-adj_density"].start();
+ tms2["2-adj_mach"].stop();
+ tms2["3-adj_density"].start();
gradientswrtDensity();
- tms["3-adj_density"].stop();
- tms["4-adj_velocity"].start();
+ tms2["3-adj_density"].stop();
+ tms2["4-adj_velocity"].start();
gradientswrtVelocity();
- tms["4-adj_velocity"].stop();
- tms["5-adj_dStar"].start();
+ tms2["4-adj_velocity"].stop();
+ tms2["5-adj_dStar"].start();
gradientswrtDeltaStar();
- tms["5-adj_dStar"].stop();
- tms["6-adj_blw"].start();
+ tms2["5-adj_dStar"].stop();
+ tms2["6-adj_blw"].start();
gradientswrtBlowingVelocity();
- tms["6-adj_blw"].stop();
- tms["7-adj_aoa"].start();
+ tms2["6-adj_blw"].stop();
+ tms2["7-adj_aoa"].start();
gradientwrtAoA();
- tms["7-adj_aoa"].stop();
- tms["8-adj_mesh"].start();
+ tms2["7-adj_aoa"].stop();
+ tms2["8-adj_mesh"].start();
gradientwrtMesh();
- tms["8-adj_mesh"].stop();
+ tms2["8-adj_mesh"].stop();
+ tms["1-Derivatives"].stop();
transposeMatrices();
- tms["9-build"].start();
- // Store pointers to the matrices in a 2D vector
+ 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},
@@ -236,7 +237,7 @@ void CoupledAdjoint::run()
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);
- tms["9-build"].stop();
+ tms2["9-build"].stop();
size_t phiIdx = 0;
size_t machIdx = dRphi_phi.cols();
@@ -245,10 +246,12 @@ void CoupledAdjoint::run()
size_t dStarIdx = vIdx + dRv_v.cols();
size_t blwIdx = dStarIdx + dRdStar_dStar.cols();
- //***** Gradient wrt angle of attack *****//
- //****************************************//
+ //**************************************************************************//
+ // Gradients wrt angle of attack //
+ //**************************************************************************//
// CL
- tms["10-solveCxAoA"].start();
+ 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
@@ -297,15 +300,18 @@ void CoupledAdjoint::run()
// Total gradient
tdCd_AoA = dCd_AoA - dRphi_AoA.transpose()*lambdaCd.segment(phiIdx, dRphi_phi.cols()); // Total gradient
- tms["10-solveCxAoA"].stop();
-
- //***** Gradient 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"
-
- tms["12-solveCxMesh"].start();
+ 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
@@ -339,8 +345,32 @@ void CoupledAdjoint::run()
tdCd_x[n->row](m) = lambdaCd_x[isol->pbl->nDim * (n->row) + m];
}
}
- tms["12-solveCxMesh"].stop();
- //std::cout << tms << std::endl;
+ tms2["12-solveCxMesh"].stop();
+ tms["2-Solve"].stop();
+
+ // Print output
+ // aoa gradients
+ std::cout << "-------------- Adjoint solution --------------" << 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 << "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;
+
+ std::cout << "--------------- Adjoint timers ---------------" << std::endl;
+ std::cout << tms << "----------------------------------------------" << std::endl;
+ if (vsol->verbose > 2)
+ std::cout << tms2 << "----------------------------------------------" << std::endl;
}
/**
@@ -350,119 +380,67 @@ void CoupledAdjoint::run()
*/
void CoupledAdjoint::gradientswrtInviscidFlow()
{
- //**** dRphi_phi ****//
+ //**************************************************************************//
+ // dRphi_phi //
+ //--------------------------------------------------------------------------//
+ // Analytical //
+ //**************************************************************************//
dRphi_phi = iadjoint->getRu_U();
- //**** dRM_phi, dRrho_phi, dRv_phi ****//
+ //**************************************************************************//
+ // dRM_phi, dRrho_phi, dRv_phi //
+ //--------------------------------------------------------------------------//
+ // Analytical //
+ //**************************************************************************//
auto pbl = isol->pbl;
- // Finite difference
- std::vector<double> Phi_save = isol->phi; // Differential of the independent variable (phi)
- double deltaPhi = 1.e-8; // perturbation
-
- std::vector<Eigen::Triplet<double>> tripletsdM;
- std::vector<Eigen::Triplet<double>> tripletsdrho;
- std::vector<Eigen::Triplet<double>> tripletsdv;
-
- for (auto n : pbl->msh->nodes){
- Phi_save[n->row] = isol->phi[n->row] + deltaPhi;
-
- // Recompute the quantities on surface nodes.
- size_t jj = 0;
- for (auto sec: vsol->sections[0]){
- int iReg = 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 dM = 0.; double drho = 0.;
- Eigen::Vector3d dv = Eigen::Vector3d::Zero();
- // Average of the quantity on the elements adjacent to the considered node.
- for (auto e : pbl->fluid->neMap[srfNode]){
- dM += pbl->fluid->mach->eval(*e, Phi_save, 0);
- drho += pbl->fluid->rho->eval(*e, Phi_save, 0);
- Eigen::VectorXd grad = e->computeGradient(Phi_save, 0);
- for (int i = 0; i < grad.size(); ++i)
- dv(i) += grad(i);
+
+ 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));
}
- dM /= static_cast<double>(pbl->fluid->neMap[srfNode].size());
- drho /= static_cast<double>(pbl->fluid->neMap[srfNode].size());
- dv /= static_cast<double>(pbl->fluid->neMap[srfNode].size());
-
- // Form triplets
- // Minus because M = f(phi) and therefore RM = M - f(phi) = 0
- tripletsdM.push_back(Eigen::Triplet<double>(jj, n->row, -(dM - isol->M[srfNode->row])/deltaPhi));
- tripletsdrho.push_back(Eigen::Triplet<double>(jj, n->row, -(drho - isol->rho[srfNode->row])/deltaPhi));
- for (int i = 0; i < pbl->nDim; ++i)
- tripletsdv.push_back(Eigen::Triplet<double>(jj*pbl->nDim + i, n->row, -(dv(i) - isol->U[srfNode->row](i))/deltaPhi));
- jj++;
}
+ ++jj;
}
- // Reset phi
- Phi_save[n->row] = isol->phi[n->row];
}
- // Fill matrices with gradients
- dRM_phi.setFromTriplets(tripletsdM.begin(), tripletsdM.end());
- dRrho_phi.setFromTriplets(tripletsdrho.begin(), tripletsdrho.end());
- dRv_phi.setFromTriplets(tripletsdv.begin(), tripletsdv.end());
- // Remove zeros
+ 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.);
- // // Analytical
- // std::vector<Eigen::Triplet<double>> tripletsdM_an;
-
- // 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]];
- // for (auto e: pbl->fluid->neMap[srfNode])
- // {
- // // Get pre-computed values
- // tbox::Cache &cache = e->getVCache();
- // tbox::Gauss &gauss = cache.getVGauss();
- // Eigen::VectorXd grad = Eigen::VectorXd::Zero(e->nodes.size());
- // for (size_t k = 0; k < gauss.getN(); ++k)
- // {
- // // Gauss point and determinant
- // double wdj = gauss.getW(k) * e->getDetJ(k);
- // // Shape functions and solution gradient
- // Eigen::MatrixXd const &dSf = e->getJinv(k) * cache.getDsf(k);
- // Eigen::VectorXd dPhi = e->computeGradient(isol->phi, k);
- // grad += isol->pbl->fluid->mach->evalGrad(*e, isol->phi, k) * dSf.transpose() * dPhi;
- // }
- // for (size_t k = 0; k < e->nodes.size(); ++k)
- // {
- // tripletsdM_an.push_back(Eigen::Triplet<double>(jj, e->nodes[k]->row, -grad(k)/e->nodes.size()));
- // }
- // }
- // ++jj;
- // }
- // }
- // Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_phi_an(dRM_phi.rows(), dRM_phi.cols());
- // dRM_phi_an.setFromTriplets(tripletsdM_an.begin(), tripletsdM_an.end());
- // dRM_phi_an.prune(0.);
- // writeMatrixToFile(dRM_phi_an, "dRM_phi_an.txt");
- // writeMatrixToFile(dRM_phi, "dRM_phi.txt");
- // writeMatrixToFile((dRM_phi_an - dRM_phi), "dif.txt");
- // std::cout << "written" << std::endl;
- // std::cout << dRM_phi_an - dRM_phi << std::endl;
-
- // throw;
-
-
- //**** partials dCl_phi, dCd_phi ****//
+ //**************************************************************************//
+ // 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());
@@ -475,10 +453,16 @@ void CoupledAdjoint::gradientswrtInviscidFlow()
*/
void CoupledAdjoint::gradientswrtMachNumber()
{
- /**** dRM_M ****/
+ //**************************************************************************//
+ // dRM_M //
+ //**************************************************************************//
dRM_M.setIdentity();
- //**** dRdStar_M ****//
+ //**************************************************************************//
+ // dRdStar_M //
+ //--------------------------------------------------------------------------//
+ // Central finite-difference //
+ //**************************************************************************//
std::deque<Eigen::Triplet<double>> T;
double delta = 1e-6;
size_t i = 0;
@@ -524,10 +508,16 @@ void CoupledAdjoint::gradientswrtDensity()
for (auto bBC : isol->pbl->bBCs)
nEs += bBC->uE.size(); // Number of blowing elements
- //**** dRrho_rho ****//
+ //**************************************************************************//
+ // dRrho_rho //
+ //**************************************************************************//
dRrho_rho.setIdentity();
- //**** dRblw_rho ****//
+ //**************************************************************************//
+ // dRblw_rho //
+ //--------------------------------------------------------------------------//
+ // Analytical //
+ //**************************************************************************//
std::deque<Eigen::Triplet<double>> T_blw;
int offSetElms = 0;
int offSetNodes = 0;
@@ -560,7 +550,9 @@ void CoupledAdjoint::gradientswrtVelocity()
for (auto bBC : isol->pbl->bBCs)
nEs += bBC->uE.size(); // Number of blowing elements
- /**** dRv_v ****/
+ //**************************************************************************//
+ // dRv_v //
+ //**************************************************************************//
dRv_v.setIdentity();
//**** Get velocity****/
@@ -583,7 +575,11 @@ void CoupledAdjoint::gradientswrtVelocity()
}
}
- //**** dvt_v (analytical) ****//
+ //**************************************************************************//
+ // dvt_v //
+ //--------------------------------------------------------------------------//
+ // Analytical //
+ //**************************************************************************//
Eigen::SparseMatrix<double, Eigen::RowMajor> dVt_v(nNs, nNs * nDim);
std::deque<Eigen::Triplet<double>> T_vt;
@@ -599,7 +595,11 @@ void CoupledAdjoint::gradientswrtVelocity()
dVt_v.setFromTriplets(T_vt.begin(), T_vt.end());
dVt_v.prune(0.);
- //**** dRdStar_vt (finite-difference) ****//
+ //**************************************************************************//
+ // dRdStar_vt //
+ //--------------------------------------------------------------------------//
+ // Central finite-difference //
+ //**************************************************************************//
std::deque<Eigen::Triplet<double>> T;
double delta = 1e-6;
i = 0;
@@ -630,14 +630,22 @@ void CoupledAdjoint::gradientswrtVelocity()
dRdStar_vt.setFromTriplets(T.begin(), T.end());
dRdStar_vt.prune(0.);
- //**** dRdStar_v ****//
+ //**************************************************************************//
+ // dRdStar_v //
+ //**************************************************************************//
dRdStar_v = dRdStar_vt * dVt_v;
dRdStar_v.prune(0.);
- //**** dRCdf_v ****//
+ //**************************************************************************//
+ // dCdf_v //
+ //**************************************************************************//
dCd_v = dCd_vt.transpose() * dVt_v; // [1xnNs] x [nNs x nNs*nDim] = [1 x nNs*nDim]
- //**** dRblw_vt (analytical) ****//
+ //**************************************************************************//
+ // dRblw_vt //
+ //--------------------------------------------------------------------------//
+ // Analytical //
+ //**************************************************************************//
std::deque<Eigen::Triplet<double>> T_blw;
int offSetElms = 0;
int offSetNodes = 0;
@@ -655,7 +663,9 @@ void CoupledAdjoint::gradientswrtVelocity()
dRblw_vt.setFromTriplets(T_blw.begin(), T_blw.end());
dRblw_vt.prune(0.);
- //**** dRblw_v ****//
+ //**************************************************************************//
+ // dRblw_v //
+ //**************************************************************************//
dRblw_v = dRblw_vt * dVt_v;
dRblw_v.prune(0.);
}
@@ -667,7 +677,11 @@ void CoupledAdjoint::gradientswrtVelocity()
*/
void CoupledAdjoint::gradientswrtDeltaStar()
{
- //**** dRdStar_dStar (finite-difference) ****//
+ //**************************************************************************//
+ // dRdStar_dStar //
+ //--------------------------------------------------------------------------//
+ // Central finite-difference //
+ //**************************************************************************//
std::deque<Eigen::Triplet<double>> T;
double delta = 1e-6;
double saveDStar = 0.;
@@ -701,7 +715,11 @@ void CoupledAdjoint::gradientswrtDeltaStar()
dRdStar_dStar -= dRdStar_dStar_dum;
dRdStar_dStar.prune(0.);
- //**** dRblw_dStar (analytical) ****//
+ //**************************************************************************//
+ // dRblw_dStar //
+ //--------------------------------------------------------------------------//
+ // Analytical //
+ //**************************************************************************//
std::deque<Eigen::Triplet<double>> T_blw;
int offSetElms = 0;
int offSetNodes = 0;
@@ -725,7 +743,11 @@ void CoupledAdjoint::gradientswrtDeltaStar()
*/
void CoupledAdjoint::gradientswrtBlowingVelocity()
{
- /**** dRphi_blw ****/
+ //**************************************************************************//
+ // dRphi_blw //
+ //--------------------------------------------------------------------------//
+ // Analytical //
+ //**************************************************************************//
std::deque<Eigen::Triplet<double>> T;
size_t jj = 0;
@@ -762,9 +784,19 @@ void CoupledAdjoint::gradientswrtBlowingVelocity()
* 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();
}
@@ -775,22 +807,43 @@ void CoupledAdjoint::gradientwrtAoA(){
*/
void CoupledAdjoint::gradientwrtMesh(){
- /**** dCk_dx ****/
+ 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 ****/
+ //**************************************************************************//
+ // dRphi_x //
+ //--------------------------------------------------------------------------//
+ // Analytical //
+ //**************************************************************************//
dRphi_x = iadjoint->getRu_X();
- /**** dRx_x ****/
+ //**************************************************************************//
+ // dRx_x //
+ //--------------------------------------------------------------------------//
+ // Analytical //
+ //**************************************************************************//
dRx_x = iadjoint->getRx_X();
- /**** dRdStar_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;
- size_t i = 0;
int iReg = 0;
- double saveX = 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])
@@ -805,7 +858,7 @@ void CoupledAdjoint::gradientwrtMesh(){
{
int rowj = isol->pbl->bBCs[iReg]->nodes[sec->rows[j]]->row;
// x
- saveX = sec->x[j];
+ double saveX = sec->x[j];
sec->x[j] = saveX + delta;
sec->computeLocalCoordinates();
@@ -819,36 +872,43 @@ void CoupledAdjoint::gradientwrtMesh(){
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
- saveX = sec->y[j];
+ double saveY = sec->y[j];
- sec->y[j] = saveX + delta;
+ sec->y[j] = saveY + delta;
sec->computeLocalCoordinates();
sec->xxExt = sec->loc;
ddStar[0] = this->runViscous();
dCdf[0] = vsol->Cdf;
- sec->y[j] = saveX - delta;
+ sec->y[j] = saveY - delta;
sec->computeLocalCoordinates();
sec->xxExt = sec->loc;
ddStar[1] = this->runViscous();
dCdf[1] = vsol->Cdf;
- sec->y[j] = saveX;
+ 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);
- ++i;
}
}
dRdStar_x.setFromTriplets(T_dStar_x.begin(), T_dStar_x.end());
dRdStar_x.prune(0.);
- /**** dRblw_x ****/
+ //**************************************************************************//
+ // dRblw_x //
+ //--------------------------------------------------------------------------//
+ // Analytical //
+ //**************************************************************************//
std::deque<Eigen::Triplet<double>> T_blw_x;
int offSetElms = 0;
for (const auto& sec: vsol->sections[0])
@@ -873,167 +933,69 @@ void CoupledAdjoint::gradientwrtMesh(){
dRblw_x.setFromTriplets(T_blw_x.begin(), T_blw_x.end());
dRblw_x.prune(0.);
- /**** dRM_x, dRrho_x, dRv_x ****/
- double delta_x = 1e-8;
+ //**************************************************************************//
+ // dRM_x, dRrho_x, dRv_x //
+ //--------------------------------------------------------------------------//
+ // Analytical //
+ //**************************************************************************//
auto pbl = isol->pbl;
- std::vector<Eigen::Triplet<double>> tripletsdM;
- std::vector<Eigen::Triplet<double>> tripletsdrho;
- std::vector<Eigen::Triplet<double>> tripletsdv;
-
- for (auto n : pbl->msh->nodes){
- for (int idim = 0; idim < isol->pbl->nDim; ++idim){
- // Modify node position
- double save = n->pos[idim];
- n->pos[idim] = save + delta_x;
- for (auto e: pbl->msh->elems)
- if (e->type() == tbox::ElType::LINE2 || e->type() == tbox::ElType::TRI3) e->update();
-
- // Recompute the quantities on surface nodes.
- size_t jj = 0;
- for (auto sec: vsol->sections[0]){
- int iReg = 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 dM = 0.; double drho = 0.;
- Eigen::Vector3d dv = Eigen::Vector3d::Zero();
- // Average of the quantity on the elements adjacent to the considered node.
- for (auto e : pbl->fluid->neMap[srfNode]){
- dM += pbl->fluid->mach->eval(*e, isol->phi, 0);
- drho += pbl->fluid->rho->eval(*e, isol->phi, 0);
- Eigen::VectorXd grad = e->computeGradient(isol->phi, 0);
- for (int i = 0; i < grad.size(); ++i)
- dv(i) += grad(i);
- }
- dM /= static_cast<double>(pbl->fluid->neMap[srfNode].size());
- drho /= static_cast<double>(pbl->fluid->neMap[srfNode].size());
- dv /= static_cast<double>(pbl->fluid->neMap[srfNode].size());
-
- // Form triplets
- // Minus because M = f(phi) and therefore RM = M - f(phi) = 0
- tripletsdM.push_back(Eigen::Triplet<double>(jj, n->row*pbl->nDim + idim, -(dM - isol->M[srfNode->row])/delta_x));
- tripletsdrho.push_back(Eigen::Triplet<double>(jj, n->row*pbl->nDim + idim, -(drho - isol->rho[srfNode->row])/delta_x));
- for (int i = 0; i < pbl->nDim; ++i)
- tripletsdv.push_back(Eigen::Triplet<double>(jj*pbl->nDim + i, n->row*pbl->nDim + idim, -(dv(i) - isol->U[srfNode->row](i))/delta_x));
- jj++;
- }
- }
- // Reset node pos
- n->pos[idim] = save;
- }
- }
- // Fill matrices with gradients
- dRM_x.setFromTriplets(tripletsdM.begin(), tripletsdM.end());
- dRrho_x.setFromTriplets(tripletsdrho.begin(), tripletsdrho.end());
- dRv_x.setFromTriplets(tripletsdv.begin(), tripletsdv.end());
- // Remove zeros
- dRM_x.prune(0.); dRrho_x.prune(0.); dRv_x.prune(0.);
-}
-
-void CoupledAdjoint::setdRdStar_M(std::vector<std::vector<double>> _dRdStar_dM){
- std::vector<Eigen::Triplet<double>> triplets;
- for (size_t i = 0; i < _dRdStar_dM.size(); i++){
- for (size_t j = 0; j < _dRdStar_dM[i].size(); j++){
- triplets.push_back(Eigen::Triplet<double>(i, j, _dRdStar_dM[i][j]));
- }
- }
- // Set matrix
- dRdStar_M.setFromTriplets(triplets.begin(), triplets.end());
- dRdStar_M.prune(0.);
-}
-
-void CoupledAdjoint::setdRdStar_v(std::vector<std::vector<double>> _dRdStar_dv){
- std::vector<Eigen::Triplet<double>> triplets;
- for (size_t i = 0; i < _dRdStar_dv.size(); i++){
- for (size_t j = 0; j < _dRdStar_dv[i].size(); j++){
- triplets.push_back(Eigen::Triplet<double>(i, j, _dRdStar_dv[i][j]));
- }
- }
- // Set matrix
- dRdStar_v.setFromTriplets(triplets.begin(), triplets.end());
- dRdStar_v.prune(0.);
-}
-void CoupledAdjoint::setdRdStar_dStar(std::vector<std::vector<double>> _dRdStar_dStar){
- std::vector<Eigen::Triplet<double>> triplets;
- for (size_t i = 0; i < _dRdStar_dStar.size(); i++){
- for (size_t j = 0; j < _dRdStar_dStar[i].size(); j++){
- triplets.push_back(Eigen::Triplet<double>(i, j, _dRdStar_dStar[i][j]));
- }
- }
- // Set matrix
- dRdStar_dStar.setFromTriplets(triplets.begin(), triplets.end());
- dRdStar_dStar.prune(0.);
-}
+ std::deque<Eigen::Triplet<double>> T_dM_x;
+ std::deque<Eigen::Triplet<double>> T_drho_x;
+ std::deque<Eigen::Triplet<double>> T_dv_x;
-void CoupledAdjoint::setdRdStar_x(std::vector<std::vector<double>> _dRdStar_x){
- std::vector<Eigen::Triplet<double>> triplets;
- for (size_t i = 0; i < _dRdStar_x.size(); i++){
- for (size_t j = 0; j < _dRdStar_x[i].size(); j++){
- triplets.push_back(Eigen::Triplet<double>(i, j, _dRdStar_x[i][j]));
- }
- }
- // Set matrix
- dRdStar_x.setFromTriplets(triplets.begin(), triplets.end());
- dRdStar_x.prune(0.);
-}
-
-void CoupledAdjoint::setdRblw_rho(std::vector<std::vector<double>> _dRblw_rho){
- std::vector<Eigen::Triplet<double>> triplets;
- for (size_t i = 0; i < _dRblw_rho.size(); i++){
- for (size_t j = 0; j < _dRblw_rho[i].size(); j++){
- triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_rho[i][j]));
- }
- }
- // Set matrix
- dRblw_rho.setFromTriplets(triplets.begin(), triplets.end());
- dRblw_rho.prune(0.);
-}
-
-void CoupledAdjoint::setdRblw_v(std::vector<std::vector<double>> _dRblw_v){
- std::vector<Eigen::Triplet<double>> triplets;
- for (size_t i = 0; i < _dRblw_v.size(); i++){
- for (size_t j = 0; j < _dRblw_v[i].size(); j++){
- triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_v[i][j]));
- }
- }
- // Set matrix
- dRblw_v.setFromTriplets(triplets.begin(), triplets.end());
- dRblw_v.prune(0.);
-}
-
-void CoupledAdjoint::setdRblw_dStar(std::vector<std::vector<double>> _dRblw_dStar){
- std::vector<Eigen::Triplet<double>> triplets;
- for (size_t i = 0; i < _dRblw_dStar.size(); i++){
- for (size_t j = 0; j < _dRblw_dStar[i].size(); j++){
- triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_dStar[i][j]));
- }
- }
- // Set matrix
- dRblw_dStar.setFromTriplets(triplets.begin(), triplets.end());
- dRblw_dStar.prune(0.);
-}
+ int jj = 0;
+ iReg = 0;
-void CoupledAdjoint::setdRblw_x(std::vector<std::vector<double>> _dRblw_x){
- std::vector<Eigen::Triplet<double>> triplets;
- for (size_t i = 0; i < _dRblw_x.size(); i++){
- for (size_t j = 0; j < _dRblw_x[i].size(); j++){
- triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_x[i][j]));
+ 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 (int 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;
}
}
- // Set matrix
- dRblw_x.setFromTriplets(triplets.begin(), triplets.end());
- dRblw_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.);
}
+/**
+ * @brief Run the viscous solver and extract the deltaStar
+ * on all the stations
+ *
+ * @return Eigen::VectorXd deltaStar
+ */
Eigen::VectorXd CoupledAdjoint::runViscous()
{
int exitCode = vsol->run();
+ if (exitCode != 0)
+ std::cout << "vsol terminated with exit code " << exitCode << std::endl;
// Extract deltaStar
std::vector<Eigen::VectorXd> dStar_p;
@@ -1108,6 +1070,17 @@ void CoupledAdjoint::transposeMatrices(){
}
+/**
+ * @brief Write the matrices to files
+ * @param matrix Eigen matrix to write
+ * @param filename Name of the file to write
+ * @return void
+ *
+ * The matrices are written in a file with the following format:
+ * - The first line contains the column headers
+ * - The following lines contain the matrix data
+ * - 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()) {
diff --git a/blast/src/DCoupledAdjoint.h b/blast/src/DCoupledAdjoint.h
index 086a70d..2cfa4ab 100644
--- a/blast/src/DCoupledAdjoint.h
+++ b/blast/src/DCoupledAdjoint.h
@@ -116,7 +116,8 @@ private:
Eigen::VectorXd dCdp_x; ///< Partial derivative of the pressure drag coefficient with respect to the mesh coordinates
Eigen::VectorXd dCdf_x; ///< Partial derivative of the friction drag coefficient with respect to the mesh coordinates
- fwk::Timers tms; //< internal timers
+ fwk::Timers tms; //< internal timers
+ fwk::Timers tms2; //< internal precise timers
public:
CoupledAdjoint(std::shared_ptr<dart::Adjoint> iAdjoint, std::shared_ptr<blast::Driver> vSolver);
@@ -125,16 +126,6 @@ public:
void run();
void reset();
- void setdRdStar_M(std::vector<std::vector<double>> _dRdStar_M);
- void setdRdStar_v(std::vector<std::vector<double>> _dRdStar_v);
- void setdRdStar_dStar(std::vector<std::vector<double>> _dRdStar_dStar);
- void setdRdStar_x(std::vector<std::vector<double>> _dRdStar_x);
-
- void setdRblw_rho(std::vector<std::vector<double>> _dRblw_rho);
- void setdRblw_v(std::vector<std::vector<double>> _dRblw_v);
- void setdRblw_dStar(std::vector<std::vector<double>> _dRblw_ddStar);
- void setdRblw_x(std::vector<std::vector<double>> _dRblw_dx);
-
private:
void buildAdjointMatrix(std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor>*>> &matrices, Eigen::SparseMatrix<double, Eigen::RowMajor> &matrixAdjoint);
diff --git a/blast/src/DFluxes.cpp b/blast/src/DFluxes.cpp
index 3cf81cc..168b60e 100644
--- a/blast/src/DFluxes.cpp
+++ b/blast/src/DFluxes.cpp
@@ -61,11 +61,8 @@ VectorXd Fluxes::blLaws(size_t iPoint, BoundaryLayer *bl, std::vector<double> u)
{
size_t nVar = bl->getnVar();
- MatrixXd timeMatrix = MatrixXd::Zero(5, 5);
- VectorXd spaceVector = VectorXd::Zero(5);
-
// Dissipation ratio.
- double dissipRatio;
+ double dissipRatio = 0.;
if (bl->name == "wake")
dissipRatio = 0.9;
else
@@ -146,39 +143,62 @@ VectorXd Fluxes::blLaws(size_t iPoint, BoundaryLayer *bl, std::vector<double> u)
double ctEqa = bl->ctEq[iPoint];
double usa = bl->us[iPoint];
+ //--------------------------------------------------------------------------------//
+ // Integral boundary layer equations. //
+
// Space part.
- spaceVector(0) = dt_dx + (2 + u[1] - Mea * Mea) * (u[0] / u[3]) * due_dx - cfa / 2;
- spaceVector(1) = u[0] * dHstar_dx + (2 * Hstar2a + Hstara * (1 - u[1])) * u[0] / u[3] * due_dx - 2 * cda + Hstara * cfa / 2;
- spaceVector(2) = dN_dx - ax;
- spaceVector(3) = due_dx - c * (u[1] * dt_dx + u[0] * dH_dx) - dueExt_dx + cExt * ddeltaStarExt_dx;
+ // spaceVector(0) = dt_dx + (2 + u[1] - Mea * Mea) * (u[0] / u[3]) * due_dx - cfa / 2;
+ // spaceVector(1) = u[0] * dHstar_dx + (2 * Hstar2a + Hstara * (1 - u[1])) * u[0] / u[3] * due_dx - 2 * cda + Hstara * cfa / 2;
+ // spaceVector(2) = dN_dx - ax;
+ // spaceVector(3) = due_dx - c * (u[1] * dt_dx + u[0] * dH_dx) - dueExt_dx + cExt * ddeltaStarExt_dx;
+
+ // if (bl->regime[iPoint] == 1)
+ // spaceVector(4) = (2 * deltaa / u[4]) * dct_dx - 5.6 * (ctEqa - u[4] * dissipRatio) - 2 * deltaa * (4 / (3 * u[1] * u[0]) * (cfa / 2 - (((Hka - 1) / (6.7 * Hka * dissipRatio)) * ((Hka - 1) / (6.7 * Hka * dissipRatio)))) - 1 / u[3] * due_dx);
- if (bl->regime[iPoint] == 1)
- spaceVector(4) = (2 * deltaa / u[4]) * dct_dx - 5.6 * (ctEqa - u[4] * dissipRatio) - 2 * deltaa * (4 / (3 * u[1] * u[0]) * (cfa / 2 - (((Hka - 1) / (6.7 * Hka * dissipRatio)) * ((Hka - 1) / (6.7 * Hka * dissipRatio)))) - 1 / u[3] * due_dx);
+ // // Time part.
+ // timeMatrix(0, 0) = u[1] / u[3];
+ // timeMatrix(0, 1) = u[0] / u[3];
+ // timeMatrix(0, 3) = u[0] * u[1] / (u[3] * u[3]);
- // Time part.
- timeMatrix(0, 0) = u[1] / u[3];
- timeMatrix(0, 1) = u[0] / u[3];
- timeMatrix(0, 3) = u[0] * u[1] / (u[3] * u[3]);
+ // timeMatrix(1, 0) = (1 + u[1] * (1 - Hstara)) / u[3];
+ // timeMatrix(1, 1) = (1 - Hstara) * u[0] / u[3];
+ // timeMatrix(1, 3) = (2 - Hstara * u[1]) * u[0] / (u[3] * u[3]);
- timeMatrix(1, 0) = (1 + u[1] * (1 - Hstara)) / u[3];
- timeMatrix(1, 1) = (1 - Hstara) * u[0] / u[3];
- timeMatrix(1, 3) = (2 - Hstara * u[1]) * u[0] / (u[3] * u[3]);
+ // timeMatrix(2, 2) = 1.;
- timeMatrix(2, 2) = 1.;
+ // timeMatrix(3, 0) = -c * u[1];
+ // timeMatrix(3, 1) = -c * u[0];
+ // timeMatrix(3, 3) = 1.;
- timeMatrix(3, 0) = -c * u[1];
- timeMatrix(3, 1) = -c * u[0];
- timeMatrix(3, 3) = 1.;
+ // if (bl->regime[iPoint] == 1)
+ // {
+ // timeMatrix(4, 3) = 2 * deltaa / (usa * u[3] * u[3]);
+ // timeMatrix(4, 4) = 2 * deltaa / (u[3] * u[4] * usa);
+ // }
+ // else
+ // timeMatrix(4, 4) = 1.0;
+ //--------------------------------------------------------------------------------//
- if (bl->regime[iPoint] == 1)
+ VectorXd result = VectorXd::Zero(5);
+
+ if (bl->regime[iPoint] == 0)
{
- timeMatrix(4, 3) = 2 * deltaa / (usa * u[3] * u[3]);
- timeMatrix(4, 4) = 2 * deltaa / (u[3] * u[4] * usa);
+ result[0] = 0.5*(-2.0*u[0]*(1.0*u[1] - 2.0)*(c*(dH_dx*u[0] + dt_dx*u[1]) - cExt*ddeltaStarExt_dx + dueExt_dx - due_dx) + 1.0*(c*u[0]*u[1] + u[3])*(2.0*due_dx*u[0]*(2.*Hstar2a - Hstara*(u[1] - 1)) + 1.0*u[3]*(Hstara*cfa - 4*cda + 2.*dHstar_dx*u[0])) + 1.0*(2.*due_dx*u[0]*(-Mea*Mea + u[1] + 2) + u[3]*(-cfa + 2.*dt_dx))*(1.0*Hstara*c*u[0]*u[1] + 1.0*Hstara*u[3] - 2.0*c*u[0] - 1.0*u[3]))/(c*u[0]*u[1] + u[3]);
+ result[1] = 0.5*(2.0*u[0]*u[1]*(u[1] - 1)*(c*(dH_dx*u[0] + dt_dx*u[1]) - cExt*ddeltaStarExt_dx + dueExt_dx - due_dx) - 1.0*u[1]*(c*u[0]*u[1] + u[3])*(2.*due_dx*u[0]*(2.*Hstar2a - Hstara*(u[1] - 1)) + u[3]*(Hstara*cfa - 4*cda + 2.*dHstar_dx*u[0])) + 1.0*(2.*due_dx*u[0]*(-Mea*Mea + u[1] + 2) + u[3]*(-cfa + 2.*dt_dx))*(-1.0*Hstara*c*u[0]*u[1]*u[1] - 1.0*Hstara*u[1]*u[3] + 2.0*c*u[0]*u[1] + 1.0*u[1]*u[3] + 1.0*u[3]))/(u[0]*(c*u[0]*u[1] + u[3]));
+ result[2] = -1.0*ax + 1.0*dN_dx;
+ result[3] = 1.0*u[3]*(-1.0*c*(dH_dx*u[0] + dt_dx*u[1]) + 0.5*c*(2.*due_dx*u[0]*(-Mea*Mea + u[1] + 2) + u[3]*(-cfa + 2.*dt_dx)) + 1.0*cExt*ddeltaStarExt_dx - 1.0*dueExt_dx + 1.0*due_dx)/(c*u[0]*u[1] + u[3]);
+ result[4] = 0.;
+ }
+ else if (bl->regime[iPoint] == 1)
+ {
+ result[0] = (-u[0]*(u[1] - 2)*(c*(dH_dx*u[0] + dt_dx*u[1]) - cExt*ddeltaStarExt_dx + dueExt_dx - due_dx) + (c*u[0]*u[1] + u[3])*(2.*due_dx*u[0]*(2.*Hstar2a - Hstara*(u[1] - 1)) + u[3]*(Hstara*cfa - 4*cda + 2.*dHstar_dx*u[0]))/2. + (2.*due_dx*u[0]*(-Mea*Mea + u[1] + 2) + u[3]*(-cfa + 2.*dt_dx))*(Hstara*c*u[0]*u[1] + Hstara*u[3] - 2.*c*u[0] - u[3])/2.)/(c*u[0]*u[1] + u[3]);
+ result[1] = (2.*u[0]*u[1]*(u[1] - 1)*(c*(dH_dx*u[0] + dt_dx*u[1]) - cExt*ddeltaStarExt_dx + dueExt_dx - due_dx) - u[1]*(c*u[0]*u[1] + u[3])*(2.*due_dx*u[0]*(2.*Hstar2a - Hstara*(u[1] - 1)) + u[3]*(Hstara*cfa - 4*cda + 2.*dHstar_dx*u[0])) + (2.*due_dx*u[0]*(-Mea*Mea + u[1] + 2) + u[3]*(-cfa + 2.*dt_dx))*(-Hstara*c*u[0]*u[1]*u[1] - Hstara*u[1]*u[3] + 2.*c*u[0]*u[1] + u[1]*u[3] + u[3]))/(2.*u[0]*(c*u[0]*u[1] + u[3]));
+ result[2] = -ax + dN_dx;
+ result[3] = u[3]*(-2.*c*(dH_dx*u[0] + dt_dx*u[1]) + c*(2.*due_dx*u[0]*(-Mea*Mea + u[1] + 2) + u[3]*(-cfa + 2.*dt_dx)) + 2.*cExt*ddeltaStarExt_dx - 2.*dueExt_dx + 2.*due_dx)/(2.*(c*u[0]*u[1] + u[3]));
+ result[4] = (-Hka*Hka*c*deltaa*dissipRatio*dissipRatio*u[0]*u[1]*u[4]*(2.*due_dx*u[0]*(-Mea*Mea + u[1] + 2) + u[3]*(-cfa + 2.*dt_dx))/2. + Hka*Hka*deltaa*dissipRatio*dissipRatio*u[0]*u[1]*u[4]*(c*(dH_dx*u[0] + dt_dx*u[1]) - cExt*ddeltaStarExt_dx + dueExt_dx - due_dx) + usa*(c*u[0]*u[1] + u[3])*(6*Hka*Hka*dct_dx*deltaa*dissipRatio*dissipRatio*u[0]*u[1]*u[3] + 16.8*Hka*Hka*dissipRatio*dissipRatio*u[0]*u[1]*u[3]*u[4]*(-ctEqa + dissipRatio*u[4]) + 2.*deltaa*u[4]*(3*Hka*Hka*dissipRatio*dissipRatio*due_dx*u[0]*u[1] + u[3]*(-2.*Hka*Hka*cfa*dissipRatio*dissipRatio + 0.0891067052795723*(Hka - 1)*(Hka - 1))))/6)/(Hka*Hka*deltaa*dissipRatio*dissipRatio*u[0]*u[1]*(c*u[0]*u[1] + u[3]));
}
- else
- timeMatrix(4, 4) = 1.0;
- return timeMatrix.partialPivLu().solve(spaceVector);
+ return result;
}
/**
--
GitLab