diff --git a/blast/src/DCoupledAdjoint.cpp b/blast/src/DCoupledAdjoint.cpp
index 8ea9318add45b1a8ec63c746970f578903e819fb..4cb1449ab6b6acacdaaa12f225c187019582ed4d 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 086a70d31cd853d8508bf4d59e111bd4d8dac390..2cfa4ab515e890956d72b19bc3cbf31f85fbbc5f 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 3cf81cc903ffe961ded40dfeff47b0c1a17bba1d..168b60e94e1cd9263e7a802da3d9e17fe7afcebd 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;
}
/**