From 36afd89a97061df6d67276a7957ce5951d18db08 Mon Sep 17 00:00:00 2001 From: Paul Dechamps <paul.dechamps@uliege.be> Date: Mon, 11 Mar 2024 09:34:24 +0100 Subject: [PATCH] (feat) Adjoint: Debug features Linear solver not converging --- blast/src/DCoupledAdjoint.cpp | 295 ++++++++++++++++++----------- blast/tests/dart/adjointVII.py | 337 +++++++++++++++++---------------- 2 files changed, 351 insertions(+), 281 deletions(-) diff --git a/blast/src/DCoupledAdjoint.cpp b/blast/src/DCoupledAdjoint.cpp index b9bac11..47ed765 100644 --- a/blast/src/DCoupledAdjoint.cpp +++ b/blast/src/DCoupledAdjoint.cpp @@ -29,6 +29,7 @@ #include <Eigen/Sparse> #include<Eigen/SparseLU> +#include <Eigen/Dense> #include <fstream> #include <deque> #include <algorithm> @@ -91,11 +92,18 @@ void CoupledAdjoint::reset() dRv_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNs, nEs); dRblw_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nEs); + dRM_M.setIdentity(); + dRrho_rho.setIdentity(); + dRv_v.setIdentity(); + dRblw_blw.setIdentity(); + dCl_phi = Eigen::VectorXd::Zero(nNd); dCd_phi = Eigen::VectorXd::Zero(nNd); dCl_AoA = 0.0; dCd_AoA = 0.0; + tdCl_AoA = 0.0; + tdCd_AoA = 0.0; } void CoupledAdjoint::buildAdjointMatrix(std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor>*>> &matrices, Eigen::SparseMatrix<double, Eigen::RowMajor> &matrixAdjoint) @@ -158,23 +166,57 @@ void CoupledAdjoint::run() Eigen::Map<Eigen::VectorXd> lambdasCl_(lambdasCl.data(), lambdasCl.size()); // Solve coupled system + // std::cout << "A_adjoint " << A_adjoint.rows() << " " << A_adjoint.cols() << " " << std::setprecision(10) << A_adjoint.norm() << std::endl; + // std::cout << "rhsCl " << std::setprecision(10) << Eigen::Map<Eigen::VectorXd>(rhsCl.data(), rhsCl.size()).norm() << std::endl; + // std::cout << "lambdasCl_ " << std::setprecision(10) << lambdasCl_.norm() << std::endl; + + std::cout << "Supposed size " << isol->pbl->msh->nodes.size() + isol->pbl->bBCs[0]->nodes.size() + 2*isol->pbl->bBCs[0]->nodes.size() + isol->pbl->bBCs[0]->nodes.size() + isol->pbl->bBCs[0]->uE.size() << std::endl; std::cout << "A_adjoint " << A_adjoint.rows() << " " << A_adjoint.cols() << " " << std::setprecision(10) << A_adjoint.norm() << std::endl; + // std::cout << "A_adjoint " << std::setprecision(3) << A_adjoint << std::endl; std::cout << "rhsCl " << std::setprecision(10) << Eigen::Map<Eigen::VectorXd>(rhsCl.data(), rhsCl.size()).norm() << std::endl; - std::cout << "lambdasCl_ " << std::setprecision(10) << lambdasCl_.norm() << std::endl; + // for (auto i=0; i<rhsCl.size(); i++) + // { + // std::cout << i << " " << rhsCl[i] << "\n"; + // } - std::cout << "Before solver. NIter " << alinsol->getIterations() << ". Error " << alinsol->getError() << std::endl; - alinsol->mit = 500; alinsol->compute(A_adjoint, Eigen::Map<Eigen::VectorXd>(rhsCl.data(), rhsCl.size()), lambdasCl_); - std::cout << "After solver. NIter " << alinsol->getIterations() << ". Error " << alinsol->getError() << std::endl; - std::cout << "lambdasCl " << std::setprecision(10) << lambdasCl.norm() << std::endl; + // Eigen::MatrixXd denseMat = Eigen::MatrixXd(A_adjoint); + // std::cout << "det " << denseMat.determinant() << std::endl; + // std::cout << "norm inv " << denseMat.inverse().norm() << std::endl; + // Eigen::VectorXd denseVec = Eigen::Map<Eigen::VectorXd>(rhsCl.data(), rhsCl.size()); + // lambdasCl_ = denseMat.partialPivLu().solve(denseVec); + + std::cout << "nIter " << alinsol->getIterations() << std::endl; + std::cout << "error " << alinsol->getError() << std::endl; + + std::cout << "lambdasCl_ " << std::setprecision(10) << lambdasCl.norm() << std::endl; + + // Eigen::VectorXd lambdasClTest = Eigen::VectorXd::Zero(dRphi_phi.cols()); + // Eigen::Map<Eigen::VectorXd> lambdasClTest_(lambdasClTest.data(), lambdasClTest.size()); + // alinsol->compute(dRphi_phi, Eigen::Map<Eigen::VectorXd>(dCl_phi.data(), dCl_phi.size()), lambdasClTest_); + + // std::cout << "complete " << std::setprecision(10) << "lambdaClTest " << lambdasClTest.norm() << std::endl; + + // tdCl_AoA = dCl_AoA - dRphi_AoA.transpose()*lambdasClTest; + + // std::cout << "After solver. NIter " << alinsol->getIterations() << ". Error " << alinsol->getError() << std::endl; + // writeMatrixToFile(lambdasCl_, "lambdasCl4.txt"); + // std::cout << "lambdasCl " << std::setprecision(10) << lambdasCl.norm() << std::endl; + Eigen::VectorXd lambdaClPhi = lambdasCl.block(0, 0, isol->pbl->msh->nodes.size(), 1); + // writeMatrixToFile(lambdaClPhi, "lambdasClPhi4.txt"); - std::cout << "lambdasClphi " << std::setprecision(10) << lambdaClPhi.norm() << std::endl; + // std::cout << "lambdasClphi " << std::setprecision(10) << lambdaClPhi.norm() << std::endl; - // Compute total derivative dCl_dAoA - std::cout << "dCl_AoA (partial) " << std::setprecision(10) << dCl_AoA << std::endl; + // // Compute total derivative dCl_dAoA + // std::cout << "dCl_AoA (partial) " << std::setprecision(10) << dCl_AoA << std::endl; + // std::cout << "dRphi_AoA " << std::setprecision(10) << dRphi_AoA.transpose().norm() << std::endl; + std::cout << "pdCl_AoA " << std::setprecision(10) << dCl_AoA << std::endl; std::cout << "dRphi_AoA " << std::setprecision(10) << dRphi_AoA.transpose().norm() << std::endl; + std::cout << "lambdaClPhi " << std::setprecision(10) << lambdaClPhi.norm() << std::endl; tdCl_AoA = dCl_AoA - dRphi_AoA.transpose()*lambdaClPhi; + std::cout << "tdCl_AoA " << std::setprecision(10) << tdCl_AoA << std::endl; + throw; /******** CD ********/ // Build right hand side for Cd @@ -199,52 +241,52 @@ void CoupledAdjoint::gradientswrtInviscidFlow() dRphi_phi = iadjoint->getRu_U(); //**** dRM_phi, dRrho_phi, dRv_phi ****// - // auto pbl = isol->pbl; - // // Finite difference - // std::vector<double> Phi_save = isol->phi; // Differential of the independant variable (phi) - // double deltaPhi = 1e-6; // 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 srfNode : pbl->bBCs[0]->nodes){ - - // 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); - // } - // 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 - // 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 (size_t 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++; - // } - // // 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.prune(0.); dRrho_phi.prune(0.); dRv_phi.prune(0.); + auto pbl = isol->pbl; + // Finite difference + std::vector<double> Phi_save = isol->phi; // Differential of the independant variable (phi) + double deltaPhi = 1e-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 srfNode : pbl->bBCs[0]->nodes){ + + 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); + } + 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 + tripletsdM.push_back(Eigen::Triplet<double>(jj, n->row, (isol->M[srfNode->row] - dM)/deltaPhi)); + tripletsdrho.push_back(Eigen::Triplet<double>(jj, n->row, (isol->rho[srfNode->row] - drho)/deltaPhi)); + for (size_t i = 0; i < pbl->nDim; ++i) + tripletsdv.push_back(Eigen::Triplet<double>(jj*pbl->nDim + i, n->row, (isol->U[srfNode->row](i) - dv(i))/deltaPhi)); + 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.prune(0.); dRrho_phi.prune(0.); dRv_phi.prune(0.); //**** dCl_phi, dCd_phi ****// std::vector<std::vector<double>> dCx_phi = iadjoint->computeGradientCoefficientsFlow(); @@ -282,7 +324,7 @@ void CoupledAdjoint::gradientswrtVelocity(){ } void CoupledAdjoint::gradientswrtBlowingVelocity(){ - // dRphi_blw = iadjoint->getRu_Blw(); + dRphi_blw = iadjoint->getRu_Blw(); // All others are zero. } @@ -400,68 +442,93 @@ void CoupledAdjoint::transposeMatrices(){ // // (Cout all norms) - std::cout << "norm(dRphi_phi) " << dRphi_phi.norm() << std::endl; - std::cout << "norm(dRphi_M) " << dRphi_M.norm() << std::endl; - std::cout << "norm(dRphi_rho) " << dRphi_rho.norm() << std::endl; - std::cout << "norm(dRphi_v) " << dRphi_v.norm() << std::endl; - std::cout << "norm(dRphi_blw) " << dRphi_blw.norm() << std::endl; - - std::cout << "norm(dRM_phi) " << dRM_phi.norm() << std::endl; - std::cout << "norm(dRM_M) " << dRM_M.norm() << std::endl; - std::cout << "norm(dRM_rho) " << dRM_rho.norm() << std::endl; - std::cout << "norm(dRM_v) " << dRM_v.norm() << std::endl; - std::cout << "norm(dRM_blw) " << dRM_blw.norm() << std::endl; - - std::cout << "norm(dRrho_phi) " << dRrho_phi.norm() << std::endl; - std::cout << "norm(dRrho_M) " << dRrho_M.norm() << std::endl; - std::cout << "norm(dRrho_rho) " << dRrho_rho.norm() << std::endl; - std::cout << "norm(dRrho_v) " << dRrho_v.norm() << std::endl; - std::cout << "norm(dRrho_blw) " << dRrho_blw.norm() << std::endl; - - std::cout << "norm(dRv_phi) " << dRv_phi.norm() << std::endl; - std::cout << "norm(dRv_M) " << dRv_M.norm() << std::endl; - std::cout << "norm(dRv_rho) " << dRv_rho.norm() << std::endl; - std::cout << "norm(dRv_v) " << dRv_v.norm() << std::endl; - std::cout << "norm(dRv_blw) " << dRv_blw.norm() << std::endl; - - std::cout << "norm(dRblw_phi) " << dRblw_phi.norm() << std::endl; - std::cout << "norm(dRblw_M) " << dRblw_M.norm() << std::endl; - std::cout << "norm(dRblw_rho) " << dRblw_rho.norm() << std::endl; - std::cout << "norm(dRblw_v) " << dRblw_v.norm() << std::endl; - std::cout << "norm(dRblw_blw) " << dRblw_blw.norm() << std::endl; + std::cout << "norm(dRphi_phi) " << std::setprecision(8) << dRphi_phi.norm() << std::endl; + std::cout << "norm(dRphi_M) " << std::setprecision(8) << dRphi_M.norm() << std::endl; + std::cout << "norm(dRphi_rho) " << std::setprecision(8) << dRphi_rho.norm() << std::endl; + std::cout << "norm(dRphi_v) " << std::setprecision(8) << dRphi_v.norm() << std::endl; + std::cout << "norm(dRphi_blw) " << std::setprecision(8) << dRphi_blw.norm() << std::endl; + + std::cout << "norm(dRM_phi) " << std::setprecision(8) << dRM_phi.norm() << std::endl; + std::cout << "norm(dRM_M) " << std::setprecision(8) << dRM_M.norm() << std::endl; + std::cout << "norm(dRM_rho) " << std::setprecision(8) << dRM_rho.norm() << std::endl; + std::cout << "norm(dRM_v) " << std::setprecision(8) << dRM_v.norm() << std::endl; + std::cout << "norm(dRM_blw) " << std::setprecision(8) << dRM_blw.norm() << std::endl; + + std::cout << "norm(dRrho_phi) " << std::setprecision(8) << dRrho_phi.norm() << std::endl; + std::cout << "norm(dRrho_M) " << std::setprecision(8) << dRrho_M.norm() << std::endl; + std::cout << "norm(dRrho_rho) " << std::setprecision(8) << dRrho_rho.norm() << std::endl; + std::cout << "norm(dRrho_v) " << std::setprecision(8) << dRrho_v.norm() << std::endl; + std::cout << "norm(dRrho_blw) " << std::setprecision(8) << dRrho_blw.norm() << std::endl; + + std::cout << "norm(dRv_phi) " << std::setprecision(8) << dRv_phi.norm() << std::endl; + std::cout << "norm(dRv_M) " << std::setprecision(8) << dRv_M.norm() << std::endl; + std::cout << "norm(dRv_rho) " << std::setprecision(8) << dRv_rho.norm() << std::endl; + std::cout << "norm(dRv_v) " << std::setprecision(8) << dRv_v.norm() << std::endl; + std::cout << "norm(dRv_blw) " << std::setprecision(8) << dRv_blw.norm() << std::endl; + + std::cout << "norm(dRblw_phi) " << std::setprecision(8) << dRblw_phi.norm() << std::endl; + std::cout << "norm(dRblw_M) " << std::setprecision(8) << dRblw_M.norm() << std::endl; + std::cout << "norm(dRblw_rho) " << std::setprecision(8) << dRblw_rho.norm() << std::endl; + std::cout << "norm(dRblw_v) " << std::setprecision(8) << dRblw_v.norm() << std::endl; + std::cout << "norm(dRblw_blw) " << std::setprecision(8) << dRblw_blw.norm() << std::endl; // // Print all matrices after transpose - // std::cout << "dRphi_phi " << dRphi_phi.rows() << " " << dRphi_phi.cols() << "\n" << dRphi_phi << std::endl; - // std::cout << "dRphi_M " << dRphi_M.rows() << " " << dRphi_M.cols() << "\n" << dRphi_M << std::endl; - // std::cout << "dRphi_rho " << dRphi_rho.rows() << " " << dRphi_rho.cols() << "\n" << dRphi_rho << std::endl; - // std::cout << "dRphi_v " << dRphi_v.rows() << " " << dRphi_v.cols() << "\n" << dRphi_v << std::endl; - // std::cout << "dRphi_blw " << dRphi_blw.rows() << " " << dRphi_blw.cols() << "\n" << dRphi_blw << std::endl; - - // std::cout << "dRM_phi " << dRM_phi.rows() << " " << dRM_phi.cols() << "\n" << dRM_phi << std::endl; - // std::cout << "dRM_M " << dRM_M.rows() << " " << dRM_M.cols() << "\n" << dRM_M << std::endl; - // std::cout << "dRM_rho " << dRM_rho.rows() << " " << dRM_rho.cols() << "\n" << dRM_rho << std::endl; - // std::cout << "dRM_v " << dRM_v.rows() << " " << dRM_v.cols() << "\n" << dRM_v << std::endl; - // std::cout << "dRM_blw " << dRM_blw.rows() << " " << dRM_blw.cols() << "\n" << dRM_blw << std::endl; - - // std::cout << "dRrho_phi " << dRrho_phi.rows() << " " << dRrho_phi.cols() << "\n" << dRrho_phi << std::endl; - // std::cout << "dRrho_M " << dRrho_M.rows() << " " << dRrho_M.cols() << "\n" << dRrho_M << std::endl; - // std::cout << "dRrho_rho " << dRrho_rho.rows() << " " << dRrho_rho.cols() << "\n" << dRrho_rho << std::endl; - // std::cout << "dRrho_v " << dRrho_v.rows() << " " << dRrho_v.cols() << "\n" << dRrho_v << std::endl; - // std::cout << "dRrho_blw " << dRrho_blw.rows() << " " << dRrho_blw.cols() << "\n" << dRrho_blw << std::endl; - - // std::cout << "dRv_phi " << dRv_phi.rows() << " " << dRv_phi.cols() << "\n" << dRv_phi << std::endl; - // std::cout << "dRv_M " << dRv_M.rows() << " " << dRv_M.cols() << "\n" << dRv_M << std::endl; - // std::cout << "dRv_rho " << dRv_rho.rows() << " " << dRv_rho.cols() << "\n" << dRv_rho << std::endl; - // std::cout << "dRv_v " << dRv_v.rows() << " " << dRv_v.cols() << "\n" << dRv_v << std::endl; - // std::cout << "dRv_blw " << dRv_blw.rows() << " " << dRv_blw.cols() << "\n" << dRv_blw << std::endl; - - // std::cout << "dRblw_phi " << dRblw_phi.rows() << " " << dRblw_phi.cols() << "\n" << dRblw_phi << std::endl; - // std::cout << "dRblw_M " << dRblw_M.rows() << " " << dRblw_M.cols() << "\n" << dRblw_M << std::endl; - // std::cout << "dRblw_rho " << dRblw_rho.rows() << " " << dRblw_rho.cols() << "\n" << dRblw_rho << std::endl; - // std::cout << "dRblw_v " << dRblw_v.rows() << " " << dRblw_v.cols() << "\n" << dRblw_v << std::endl; - // std::cout << "dRblw_blw " << dRblw_blw.rows() << " " << dRblw_blw.cols() << "\n" << dRblw_blw << std::endl; - + // std::cout << "dRphi_phi " << dRphi_phi.rows() << " " << dRphi_phi.cols() << "\n" << dRphi_phi << std::endl; + // std::cout << "dRphi_M " << dRphi_M.rows() << " " << dRphi_M.cols() << "\n" << dRphi_M << std::endl; + // std::cout << "dRphi_rho " << dRphi_rho.rows() << " " << dRphi_rho.cols() << "\n" << dRphi_rho << std::endl; + // std::cout << "dRphi_v " << dRphi_v.rows() << " " << dRphi_v.cols() << "\n" << dRphi_v << std::endl; + // std::cout << "dRphi_blw " << dRphi_blw.rows() << " " << dRphi_blw.cols() << "\n" << dRphi_blw << std::endl; + + // std::cout << "dRM_phi " << dRM_phi.rows() << " " << dRM_phi.cols() << "\n" << dRM_phi << std::endl; + // std::cout << "dRM_M " << dRM_M.rows() << " " << dRM_M.cols() << "\n" << dRM_M << std::endl; + // std::cout << "dRM_rho " << dRM_rho.rows() << " " << dRM_rho.cols() << "\n" << dRM_rho << std::endl; + // std::cout << "dRM_v " << dRM_v.rows() << " " << dRM_v.cols() << "\n" << dRM_v << std::endl; + // std::cout << "dRM_blw " << dRM_blw.rows() << " " << dRM_blw.cols() << "\n" << dRM_blw << std::endl; + + // std::cout << "dRrho_phi " << dRrho_phi.rows() << " " << dRrho_phi.cols() << "\n" << dRrho_phi << std::endl; + // std::cout << "dRrho_M " << dRrho_M.rows() << " " << dRrho_M.cols() << "\n" << dRrho_M << std::endl; + // std::cout << "dRrho_rho " << dRrho_rho.rows() << " " << dRrho_rho.cols() << "\n" << dRrho_rho << std::endl; + // std::cout << "dRrho_v " << dRrho_v.rows() << " " << dRrho_v.cols() << "\n" << dRrho_v << std::endl; + // std::cout << "dRrho_blw " << dRrho_blw.rows() << " " << dRrho_blw.cols() << "\n" << dRrho_blw << std::endl; + + // std::cout << "dRv_phi " << dRv_phi.rows() << " " << dRv_phi.cols() << "\n" << dRv_phi << std::endl; + // std::cout << "dRv_M " << dRv_M.rows() << " " << dRv_M.cols() << "\n" << dRv_M << std::endl; + // std::cout << "dRv_rho " << dRv_rho.rows() << " " << dRv_rho.cols() << "\n" << dRv_rho << std::endl; + // std::cout << "dRv_v " << dRv_v.rows() << " " << dRv_v.cols() << "\n" << dRv_v << std::endl; + // std::cout << "dRv_blw " << dRv_blw.rows() << " " << dRv_blw.cols() << "\n" << dRv_blw << std::endl; + + // std::cout << "dRblw_phi " << dRblw_phi.rows() << " " << dRblw_phi.cols() << "\n" << dRblw_phi << std::endl; + // std::cout << "dRblw_M " << dRblw_M.rows() << " " << dRblw_M.cols() << "\n" << dRblw_M << std::endl; + // std::cout << "dRblw_rho " << dRblw_rho.rows() << " " << dRblw_rho.cols() << "\n" << dRblw_rho << std::endl; + // std::cout << "dRblw_v " << dRblw_v.rows() << " " << dRblw_v.cols() << "\n" << dRblw_v << std::endl; + // std::cout << "dRblw_blw " << dRblw_blw.rows() << " " << dRblw_blw.cols() << "\n" << dRblw_blw << std::endl; + + std::cout << "dRphi_phi " << dRphi_phi.rows() << " " << dRphi_phi.cols() << std::endl; + std::cout << "dRphi_M " << dRphi_M.rows() << " " << dRphi_M.cols() << std::endl; + std::cout << "dRphi_rho " << dRphi_rho.rows() << " " << dRphi_rho.cols() << std::endl; + std::cout << "dRphi_v " << dRphi_v.rows() << " " << dRphi_v.cols() << std::endl; + std::cout << "dRphi_blw " << dRphi_blw.rows() << " " << dRphi_blw.cols() << std::endl; + std::cout << "dRM_phi " << dRM_phi.rows() << " " << dRM_phi.cols() << std::endl; + std::cout << "dRM_M " << dRM_M.rows() << " " << dRM_M.cols() << std::endl; + std::cout << "dRM_rho " << dRM_rho.rows() << " " << dRM_rho.cols() << std::endl; + std::cout << "dRM_v " << dRM_v.rows() << " " << dRM_v.cols() << std::endl; + std::cout << "dRM_blw " << dRM_blw.rows() << " " << dRM_blw.cols() << std::endl; + std::cout << "dRrho_phi " << dRrho_phi.rows() << " " << dRrho_phi.cols() << std::endl; + std::cout << "dRrho_M " << dRrho_M.rows() << " " << dRrho_M.cols() << std::endl; + std::cout << "dRrho_rho " << dRrho_rho.rows() << " " << dRrho_rho.cols() << std::endl; + std::cout << "dRrho_v " << dRrho_v.rows() << " " << dRrho_v.cols() << std::endl; + std::cout << "dRrho_blw " << dRrho_blw.rows() << " " << dRrho_blw.cols() << std::endl; + std::cout << "dRv_phi " << dRv_phi.rows() << " " << dRv_phi.cols() << std::endl; + std::cout << "dRv_M " << dRv_M.rows() << " " << dRv_M.cols() << std::endl; + std::cout << "dRv_rho " << dRv_rho.rows() << " " << dRv_rho.cols() << std::endl; + std::cout << "dRv_v " << dRv_v.rows() << " " << dRv_v.cols() << std::endl; + std::cout << "dRv_blw " << dRv_blw.rows() << " " << dRv_blw.cols() << std::endl; + std::cout << "dRblw_phi " << dRblw_phi.rows() << " " << dRblw_phi.cols() << std::endl; + std::cout << "dRblw_M " << dRblw_M.rows() << " " << dRblw_M.cols() << std::endl; + std::cout << "dRblw_rho " << dRblw_rho.rows() << " " << dRblw_rho.cols() << std::endl; + std::cout << "dRblw_v " << dRblw_v.rows() << " " << dRblw_v.cols() << std::endl; + std::cout << "dRblw_blw " << dRblw_blw.rows() << " " << dRblw_blw.cols() << std::endl; } void CoupledAdjoint::writeMatrixToFile(const Eigen::MatrixXd& matrix, const std::string& filename) { diff --git a/blast/tests/dart/adjointVII.py b/blast/tests/dart/adjointVII.py index 5a9cf09..8f1c6d3 100644 --- a/blast/tests/dart/adjointVII.py +++ b/blast/tests/dart/adjointVII.py @@ -51,12 +51,12 @@ def cfgInviscid(nthrds, verb): return { # Options 'scenario' : 'aerodynamic', - 'task' : 'analysis', + 'task' : 'optimization', 'Threads' : nthrds, # number of threads 'Verb' : verb, # verbosity # Model (geometry or mesh) 'File' : os.path.dirname(os.path.abspath(__file__)) + '/../../models/dart/n0012.geo', # Input file containing the model - 'Pars' : {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}, # parameters for input file model + 'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'msF' : 20, 'msTe' : 0.01, 'msLe' : 0.0075}, # parameters for input file model 'Dim' : 2, # problem dimension 'Format' : 'gmsh', # save format (vtk or gmsh) # Markers @@ -79,14 +79,12 @@ def cfgInviscid(nthrds, verb): 'z_ref' : 0.0, # reference point for moment computation (z) # Numerical 'LSolver' : 'GMRES', # inner solver (Pardiso, MUMPS or GMRES) - 'G_fill' : 2, # fill-in factor for GMRES preconditioner + 'G_fill' : 1, # fill-in factor for GMRES preconditioner 'G_tol' : 1e-5, # tolerance for GMRES 'G_restart' : 50, # restart for GMRES 'Rel_tol' : 1e-6, # relative tolerance on solver residual 'Abs_tol' : 1e-8, # absolute tolerance on solver residual 'Max_it' : 20, # solver maximum number of iterations - - 'task': 'optimization' } def cfgBlast(verb): @@ -100,7 +98,7 @@ def cfgBlast(verb): 'iterPrint': 5, # int, number of iterations between outputs 'resetInv' : True, # bool, flag to reset the inviscid calculation at every iteration. 'sections' : [0], - 'xtrF' : [None, 0.4],# Forced transition location + 'xtrF' : [0.2, 0.2],# Forced transition location 'interpolator' : 'Matching', 'nDim' : 2 } @@ -135,176 +133,182 @@ def main(): print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET) tms['solver'].start() - #aeroCoeffs = coupler.run() + aeroCoeffs = coupler.run() tms['solver'].stop() - iSolverAPI.run() + #iSolverAPI.run() - tms['FDDerivatives'].start() - tms['FDDerivatives'].stop() - tms['CoupledAdjoint'].start() + print('Dart adjoint') iSolverAPI.adjointSolver.run() tms['CoupledAdjoint'].stop() dCl_dAlpha_DartAdjoint = iSolverAPI.adjointSolver.dClAoa dCd_dAlpha_DartAdjoint = iSolverAPI.adjointSolver.dCdAoa - #pyAdjoint.computeDerivatives() + print('Compute derivative') + tms['FDDerivatives'].start() + pyAdjoint.computeDerivatives() + tms['FDDerivatives'].stop() + tms['CoupledAdjoint'].start() + print('Coupled adjoint') + tms['CouledAdjoint'].stop() coupledAdjointSolver.run() print(coupledAdjointSolver.tdCl_AoA) - quit() - - print('loop') - dRphi_phi = [] - dRphi_M = [] - dRphi_rho = [] - dRphi_v = [] - dRphi_blw = [] - dRM_phi = [] - dRM_M = [] - dRM_rho = [] - dRM_v = [] - dRM_blw = [] - dRrho_phi = [] - dRrho_M = [] - dRrho_rho = [] - dRrho_v = [] - dRrho_blw = [] - dRv_phi = [] - dRv_M = [] - dRv_rho = [] - dRv_v = [] - dRv_blw = [] - dRblw_phi = [] - dRblw_M = [] - dRblw_rho = [] - dRblw_v = [] - dRblw_blw = [] - A_adjoint = [] - rhsCl = [] - for i in range(10): - coupledAdjointSolver.reset() - #pyAdjoint.computeDerivatives() - coupledAdjointSolver.run() - dRphi_phi.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRphi_phi.txt')) - dRphi_M.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRphi_M.txt')) - dRphi_rho.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRphi_rho.txt')) - dRphi_v.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRphi_v.txt')) - dRphi_blw.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRphi_blw.txt')) - dRM_phi.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRM_phi.txt')) - dRM_M.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRM_M.txt')) - dRM_rho.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRM_rho.txt')) - dRM_v.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRM_v.txt')) - dRM_blw.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRM_blw.txt')) - dRrho_phi.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRrho_phi.txt')) - dRrho_M.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRrho_M.txt')) - dRrho_rho.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRrho_rho.txt')) - dRrho_v.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRrho_v.txt')) - dRrho_blw.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRrho_blw.txt')) - dRv_phi.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRv_phi.txt')) - dRv_M.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRv_M.txt')) - dRv_rho.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRv_rho.txt')) - dRv_v.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRv_v.txt')) - dRv_blw.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRv_blw.txt')) - dRblw_phi.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRblw_phi.txt')) - dRblw_M.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRblw_M.txt')) - dRblw_rho.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRblw_rho.txt')) - dRblw_v.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRblw_v.txt')) - dRblw_blw.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRblw_blw.txt')) - A_adjoint.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/A_adjoint.txt')) - rhsCl.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/rhsCl.txt')) - print(i, coupledAdjointSolver.tdCl_AoA) - np.set_printoptions(threshold=np.inf) - - #print(dRblw_M[0] - dRblw_M[1]) - - for i in range(len(dRphi_phi)-1): - if np.any(dRphi_phi[i+1] - dRphi_phi[i] != 0): - print('dRphi_phi') - break - if np.any(dRphi_M[i+1] - dRphi_M[i] != 0): - print('dRphi_M') - break - if np.any(dRphi_rho[i+1] - dRphi_rho[i] != 0): - print('dRphi_rho') - break - if np.any(dRphi_v[i+1] - dRphi_v[i] != 0): - print('dRphi_v') - break - if np.any(dRphi_blw[i+1] - dRphi_blw[i] != 0): - print('dRphi_blw') - break - if np.any(dRM_phi[i+1] - dRM_phi[i] != 0): - print('dRM_phi') - break - if np.any(dRM_M[i+1] - dRM_M[i] != 0): - print('dRM_M') - break - if np.any(dRM_rho[i+1] - dRM_rho[i] != 0): - print('dRM_rho') - break - if np.any(dRM_v[i+1] - dRM_v[i] != 0): - print('dRM_v') - break - if np.any(dRM_blw[i+1] - dRM_blw[i] != 0): - print('dRM_blw') - break - if np.any(dRrho_phi[i+1] - dRrho_phi[i] != 0): - print('dRrho_phi') - break - if np.any(dRrho_M[i+1] - dRrho_M[i] != 0): - print('dRrho_M') - break - if np.any(dRrho_rho[i+1] - dRrho_rho[i] != 0): - print('dRrho_rho') - break - if np.any(dRrho_v[i+1] - dRrho_v[i] != 0): - print('dRrho_v') - break - if np.any(dRrho_blw[i+1] - dRrho_blw[i] != 0): - print('dRrho_blw') - break - if np.any(dRv_phi[i+1] - dRv_phi[i] != 0): - print('dRv_phi') - break - if np.any(dRv_M[i+1] - dRv_M[i] != 0): - print('dRv_M') - break - if np.any(dRv_rho[i+1] - dRv_rho[i] != 0): - print('dRv_rho') - break - if np.any(dRv_v[i+1] - dRv_v[i] != 0): - print('dRv_v') - break - if np.any(dRv_blw[i+1] - dRv_blw[i] != 0): - print('dRv_blw') - break - if np.any(dRblw_phi[i+1] - dRblw_phi[i] != 0): - print('dRblw_phi') - break - if np.any(dRblw_M[i+1] - dRblw_M[i] >= 1e-8): - print('dRblw_M') - print(np.any(dRblw_M[i+1] - dRblw_M[i] >= 1e-8)) - #break - if np.any(dRblw_rho[i+1] - dRblw_rho[i] >= 1e-8): - print('dRblw_rho') - #break - if np.any(dRblw_v[i+1] - dRblw_v[i] >= 1e-8): - print('dRblw_v') - #break - if np.any(dRblw_blw[i+1] - dRblw_blw[i] >= 1e-8): - print('dRblw_blw') - break - if np.any(A_adjoint[i+1] - A_adjoint[i] != 0): - print('A_adjoint') - break - if np.any(rhsCl[i+1] - rhsCl[i] != 0): - print('rhsCl') - break - for i in rhsCl: - print(i) - print("sucess", coupledAdjointSolver) - quit() + # quit() + + # print('loop') + # dRphi_phi = [] + # dRphi_M = [] + # dRphi_rho = [] + # dRphi_v = [] + # dRphi_blw = [] + # dRM_phi = [] + # dRM_M = [] + # dRM_rho = [] + # dRM_v = [] + # dRM_blw = [] + # dRrho_phi = [] + # dRrho_M = [] + # dRrho_rho = [] + # dRrho_v = [] + # dRrho_blw = [] + # dRv_phi = [] + # dRv_M = [] + # dRv_rho = [] + # dRv_v = [] + # dRv_blw = [] + # dRblw_phi = [] + # dRblw_M = [] + # dRblw_rho = [] + # dRblw_v = [] + # dRblw_blw = [] + # A_adjoint = [] + # rhsCl = [] + # for i in range(10): + # coupledAdjointSolver.reset() + # #pyAdjoint.computeDerivatives() + # coupledAdjointSolver.run() + # # dRphi_phi.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRphi_phi.txt')) + # # dRphi_M.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRphi_M.txt')) + # # dRphi_rho.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRphi_rho.txt')) + # # dRphi_v.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRphi_v.txt')) + # # dRphi_blw.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRphi_blw.txt')) + # # dRM_phi.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRM_phi.txt')) + # # dRM_M.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRM_M.txt')) + # # dRM_rho.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRM_rho.txt')) + # # dRM_v.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRM_v.txt')) + # # dRM_blw.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRM_blw.txt')) + # # dRrho_phi.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRrho_phi.txt')) + # # dRrho_M.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRrho_M.txt')) + # # dRrho_rho.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRrho_rho.txt')) + # # dRrho_v.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRrho_v.txt')) + # # dRrho_blw.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRrho_blw.txt')) + # # dRv_phi.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRv_phi.txt')) + # # dRv_M.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRv_M.txt')) + # # dRv_rho.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRv_rho.txt')) + # # dRv_v.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRv_v.txt')) + # # dRv_blw.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRv_blw.txt')) + # # dRblw_phi.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRblw_phi.txt')) + # # dRblw_M.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRblw_M.txt')) + # # dRblw_rho.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRblw_rho.txt')) + # # dRblw_v.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRblw_v.txt')) + # # dRblw_blw.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRblw_blw.txt')) + # # A_adjoint.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/A_adjoint.txt')) + # # rhsCl.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/rhsCl.txt')) + # print(i, coupledAdjointSolver.tdCl_AoA) + # if abs(coupledAdjointSolver.tdCl_AoA - 5.499) >= 1e-3: + # print('MINCE ', coupledAdjointSolver.tdCl_AoA) + # np.set_printoptions(threshold=np.inf) + + # #print(dRblw_M[0] - dRblw_M[1]) + + # for i in range(len(dRphi_phi)-1): + # if np.any(dRphi_phi[i+1] - dRphi_phi[i] != 0): + # print('dRphi_phi') + # break + # if np.any(dRphi_M[i+1] - dRphi_M[i] != 0): + # print('dRphi_M') + # break + # if np.any(dRphi_rho[i+1] - dRphi_rho[i] != 0): + # print('dRphi_rho') + # break + # if np.any(dRphi_v[i+1] - dRphi_v[i] != 0): + # print('dRphi_v') + # break + # if np.any(dRphi_blw[i+1] - dRphi_blw[i] != 0): + # print('dRphi_blw') + # break + # if np.any(dRM_phi[i+1] - dRM_phi[i] != 0): + # print('dRM_phi') + # break + # if np.any(dRM_M[i+1] - dRM_M[i] != 0): + # print('dRM_M') + # break + # if np.any(dRM_rho[i+1] - dRM_rho[i] != 0): + # print('dRM_rho') + # break + # if np.any(dRM_v[i+1] - dRM_v[i] != 0): + # print('dRM_v') + # break + # if np.any(dRM_blw[i+1] - dRM_blw[i] != 0): + # print('dRM_blw') + # break + # if np.any(dRrho_phi[i+1] - dRrho_phi[i] != 0): + # print('dRrho_phi') + # break + # if np.any(dRrho_M[i+1] - dRrho_M[i] != 0): + # print('dRrho_M') + # break + # if np.any(dRrho_rho[i+1] - dRrho_rho[i] != 0): + # print('dRrho_rho') + # break + # if np.any(dRrho_v[i+1] - dRrho_v[i] != 0): + # print('dRrho_v') + # break + # if np.any(dRrho_blw[i+1] - dRrho_blw[i] != 0): + # print('dRrho_blw') + # break + # if np.any(dRv_phi[i+1] - dRv_phi[i] != 0): + # print('dRv_phi') + # break + # if np.any(dRv_M[i+1] - dRv_M[i] != 0): + # print('dRv_M') + # break + # if np.any(dRv_rho[i+1] - dRv_rho[i] != 0): + # print('dRv_rho') + # break + # if np.any(dRv_v[i+1] - dRv_v[i] != 0): + # print('dRv_v') + # break + # if np.any(dRv_blw[i+1] - dRv_blw[i] != 0): + # print('dRv_blw') + # break + # if np.any(dRblw_phi[i+1] - dRblw_phi[i] != 0): + # print('dRblw_phi') + # break + # if np.any(dRblw_M[i+1] - dRblw_M[i] >= 1e-8): + # print('dRblw_M') + # print(np.any(dRblw_M[i+1] - dRblw_M[i] >= 1e-8)) + # #break + # if np.any(dRblw_rho[i+1] - dRblw_rho[i] >= 1e-8): + # print('dRblw_rho') + # #break + # if np.any(dRblw_v[i+1] - dRblw_v[i] >= 1e-8): + # print('dRblw_v') + # #break + # if np.any(dRblw_blw[i+1] - dRblw_blw[i] >= 1e-8): + # print('dRblw_blw') + # break + # if np.any(A_adjoint[i+1] - A_adjoint[i] != 0): + # print('A_adjoint') + # break + # if np.any(rhsCl[i+1] - rhsCl[i] != 0): + # print('rhsCl') + # break + # for i in rhsCl: + # print(i) + # print("sucess", coupledAdjointSolver) + # quit() + dalpha = 1*np.pi/180 - dalpha = 1e-4 aoaV = [alpha-dalpha, alpha+dalpha] clV = [] cdV = [] @@ -318,7 +322,6 @@ def main(): dCl_dalpha_FD_coupled = (clV[-1] - clV[0]) / (2*dalpha) dCd_dalpha_FD_coupled = (cdV[-1] - cdV[0]) / (2*dalpha) - dalpha = 1e-4 aoaV = [alpha-dalpha, alpha+dalpha] clV = [] cdV = [] -- GitLab