From 9c8baff8a72eb06584587ec9f2c8ebb7afdb3d0d Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Wed, 13 Mar 2024 11:32:44 +0100
Subject: [PATCH] (feat) Adjoint working
Fixed node loop in DartAdjointInterface.
---
blast/interfaces/DSolversInterface.py | 22 +-
blast/interfaces/dart/DartAdjointInterface.py | 56 ++---
blast/interfaces/dart/DartInterface.py | 8 +-
blast/src/DCoupledAdjoint.cpp | 202 +++++++++---------
blast/src/DCoupledAdjoint.h | 1 +
blast/src/DDriver.cpp | 2 +-
blast/tests/dart/adjointVII.py | 43 +++-
7 files changed, 188 insertions(+), 146 deletions(-)
diff --git a/blast/interfaces/DSolversInterface.py b/blast/interfaces/DSolversInterface.py
index 5722230..fcf09e9 100644
--- a/blast/interfaces/DSolversInterface.py
+++ b/blast/interfaces/DSolversInterface.py
@@ -48,7 +48,7 @@ class SolversInterface:
for iSec in range(self.cfg['nSections'])]
#self.vtOld = [[np.zeros(0) for iReg in range(3)] for iSec in range(self.cfg['nSections'])]
- def setViscousSolver(self, couplIter):
+ def setViscousSolver(self, couplIter, _adj=0):
"""Interpolate solution to viscous mesh and sets the inviscid boundary in the viscous solver.
"""
self.interpolator.inviscidToViscous(self.iBnd, self.vBnd, couplIter)
@@ -118,12 +118,13 @@ class SolversInterface:
raise RuntimeError('Tried to initialize a viscous\
struture but name was', reg.name)
## External variables
- for iRegv in range(3):
- self.vSolver.setxxExt(iSec, iRegv, self.xxExt[iSec][iRegv])
- self.vSolver.setDeltaStarExt(iSec, iRegv,
- self.deltaStarExt[iSec][iRegv])
+ if _adj == 0:
+ for iRegv in range(3):
+ self.vSolver.setxxExt(iSec, iRegv, self.xxExt[iSec][iRegv])
+ self.vSolver.setDeltaStarExt(iSec, iRegv,
+ self.deltaStarExt[iSec][iRegv])
#self.vSolver.setVtOld(iSec, iRegv, self.vtOld[iSec][iRegv])
- def getBlowingBoundary(self):
+ def getBlowingBoundary(self, _adj=0):
for iSec in range(len(self.vBnd)):
for reg in self.vBnd[iSec]:
if 'vWake' in reg.name:
@@ -131,10 +132,11 @@ class SolversInterface:
elif 'vAirfoil' in reg.name:
reg.blowingVel = np.concatenate((np.flip(np.asarray(self.vSolver.extractBlowingVel(iSec, 0))),
np.asarray(self.vSolver.extractBlowingVel(iSec, 1))))
- for iRegv in range(3):
- self.xxExt[iSec][iRegv] = self.vSolver.extractxx(iSec, iRegv)
- self.deltaStarExt[iSec][iRegv] = self.vSolver.extractDeltaStar(iSec, iRegv)
- #self.vtOld[iSec][iRegv] = self.vSolver.extractUe(iSec, iRegv)
+ if _adj == 0:
+ for iRegv in range(3):
+ self.xxExt[iSec][iRegv] = self.vSolver.extractxx(iSec, iRegv)
+ self.deltaStarExt[iSec][iRegv] = self.vSolver.extractDeltaStar(iSec, iRegv)
+ #self.vtOld[iSec][iRegv] = self.vSolver.extractUe(iSec, iRegv)
self.interpolator.viscousToInviscid(self.iBnd, self.vBnd)
diff --git a/blast/interfaces/dart/DartAdjointInterface.py b/blast/interfaces/dart/DartAdjointInterface.py
index f8ece62..c1aa816 100644
--- a/blast/interfaces/dart/DartAdjointInterface.py
+++ b/blast/interfaces/dart/DartAdjointInterface.py
@@ -1,4 +1,5 @@
import numpy as np
+from matplotlib import pyplot as plt
class DartAdjointInterface():
def __init__(self, _coupledAdjointSolver, _iSolverAPI, _vSolver):
@@ -10,16 +11,18 @@ class DartAdjointInterface():
nDim = self.iSolverAPI.solver.pbl.nDim
nNd = self.iSolverAPI.solver.pbl.msh.nodes.size()
nEs = len(self.iSolverAPI.iBnd[0].elemsCoord[:,0])
- nNs = len(self.iSolverAPI.iBnd[0].nodesCoord[:,0])
- print(f"{nDim=}")
- print(f"{nNd=}")
- print(f"{nEs=}")
- print(f"{nNs=}")
+ nNs = len(self.iSolverAPI.blw[0].nodes)
dRblw_M = np.zeros((nEs, nNs))
dRblw_rho = np.zeros((nEs, nNs))
dRblw_v = np.zeros((nEs, nDim*nNs))
- delta = 1e-4
+ #self.__runViscous()
+
+ # self.deltaStarExt[0][iRegv]
+ # for iRegv in range(2):
+ # self.deltaStarExt[0][iRegv] = self.vSolver.extractDeltaStar(0, iRegv)
+
+ delta = 1e-5
saveBlw = np.zeros(nEs)
for i in range(nEs):
@@ -28,48 +31,51 @@ class DartAdjointInterface():
saveMach = np.zeros(nNs)
saverho = np.zeros(nNs)
savev = np.zeros(nDim*nNs)
- for i in range(nNs):
- saveMach[i] = self.iSolverAPI.solver.M[i]
- saverho[i] = self.iSolverAPI.solver.rho[i]
- for j in range(nDim):
- savev[i*nDim + j] = self.iSolverAPI.solver.U[i][j]
+ for iNo in range(nNs):
+ irow = self.iSolverAPI.blw[0].nodes[iNo].row
+ saveMach[iNo] = self.iSolverAPI.solver.M[irow]
+ saverho[iNo] = self.iSolverAPI.solver.rho[irow]
+ for jj in range(nDim):
+ savev[iNo*nDim + jj] = self.iSolverAPI.solver.U[irow][jj]
# Mach
for iNo in range(nNs):
- self.iSolverAPI.solver.M[i] += delta
+ irow = self.iSolverAPI.blw[0].nodes[iNo].row
+ self.iSolverAPI.solver.M[irow] += delta
blowing = self.__runViscous()
dRblw_M[:,iNo] = (blowing-saveBlw)/delta
- self.iSolverAPI.solver.M[i] = saveMach[i]
+ self.iSolverAPI.solver.M[irow] = saveMach[iNo]
# Rho
for iNo in range(nNs):
- self.iSolverAPI.solver.rho[i] += delta
+ irow = self.iSolverAPI.blw[0].nodes[iNo].row
+ self.iSolverAPI.solver.rho[irow] += delta
blowing = self.__runViscous()
dRblw_rho[:,iNo] = (blowing-saveBlw)/delta
- self.iSolverAPI.solver.rho[i] = saverho[i]
-
- # V
+ self.iSolverAPI.solver.rho[irow] = saverho[iNo]
+
+ # # V
for iNo in range(nNs):
+ irow = self.iSolverAPI.blw[0].nodes[iNo].row
for iDim in range(nDim):
- self.iSolverAPI.solver.U[iNo][iDim] += delta
+ self.iSolverAPI.solver.U[irow][iDim] += delta
blowing = self.__runViscous()
dRblw_v[:,iNo*nDim + iDim] = (blowing-saveBlw)/delta
- self.iSolverAPI.solver.U[iNo][iDim] = savev[iNo*nDim + iDim]
-
- # dRblw_M[abs(dRblw_M) < 1e-10] = 0.
- # dRblw_rho[abs(dRblw_rho) < 1e-10] = 0.
- # dRblw_v[abs(dRblw_v) < 1e-10] = 0.
+ self.iSolverAPI.solver.U[irow][iDim] = savev[iNo*nDim + iDim]
self.coupledAdjointSolver.setdRblw_M(dRblw_M)
self.coupledAdjointSolver.setdRblw_rho(dRblw_rho)
self.coupledAdjointSolver.setdRblw_v(dRblw_v)
def __runViscous(self):
- self.iSolverAPI.imposeInvBC(1)
+ self.iSolverAPI.imposeInvBC(1, adj=1)
# Run viscous solver.
exitCode = self.vSolver.run(1)
+ if exitCode != 0:
+ print('SOLVER NOT CONVERGED')
+ quit()
# Impose blowing boundary condition in the inviscid solver.
- self.iSolverAPI.imposeBlwVel()
+ self.iSolverAPI.imposeBlwVel(adj=1)
blowing = np.zeros(len(self.iSolverAPI.iBnd[0].elemsCoord[:,0]))
for i in range(len(self.iSolverAPI.iBnd[0].elemsCoord[:,0])):
diff --git a/blast/interfaces/dart/DartInterface.py b/blast/interfaces/dart/DartInterface.py
index b6fecc9..e2cc4f7 100644
--- a/blast/interfaces/dart/DartInterface.py
+++ b/blast/interfaces/dart/DartInterface.py
@@ -122,7 +122,7 @@ class DartInterface(SolversInterface):
def reset(self):
self.solver.reset()
- def imposeInvBC(self, couplIter):
+ def imposeInvBC(self, couplIter, adj=0):
""" Extract the inviscid paramaters at the edge of the boundary layer
and use it as a boundary condition in the viscous solver.
@@ -140,13 +140,13 @@ class DartInterface(SolversInterface):
self.iBnd[n].Rho[iRow] = self.solver.rho[row]
if self.cfg['nDim']==3:
self.iBnd[n].V[:,[1,2]] = self.iBnd[n].V[:,[2,1]]
- self.setViscousSolver(couplIter)
+ self.setViscousSolver(couplIter, adj)
- def imposeBlwVel(self):
+ def imposeBlwVel(self, adj=0):
""" Extract the solution of the viscous calculation (blowing velocity)
and use it as a boundary condition in the inviscid solver.
"""
- self.getBlowingBoundary()
+ self.getBlowingBoundary(_adj=adj)
# Impose blowing velocities.
for n in range(len(self.iBnd)):
if self.cfg['nDim'] == 2:
diff --git a/blast/src/DCoupledAdjoint.cpp b/blast/src/DCoupledAdjoint.cpp
index 47ed765..96f1a49 100644
--- a/blast/src/DCoupledAdjoint.cpp
+++ b/blast/src/DCoupledAdjoint.cpp
@@ -30,9 +30,11 @@
#include <Eigen/Sparse>
#include<Eigen/SparseLU>
#include <Eigen/Dense>
+#include<Eigen/SparseCholesky>
#include <fstream>
#include <deque>
#include <algorithm>
+#include <iomanip>
using namespace blast;
@@ -112,24 +114,24 @@ void CoupledAdjoint::buildAdjointMatrix(std::vector<std::vector<Eigen::SparseMat
std::vector<Eigen::Triplet<double>> triplets;
// Iterate over the rows and columns of the vector
-int rowOffset = 0;
-for (const auto& row : matrices) {
- int colOffset = 0;
- for (const auto& matrix : row) {
- // Add the triplets of the matrix to the list of triplets for the larger matrix
- for (int k = 0; k < matrix->outerSize(); ++k) {
- for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(*matrix, k); it; ++it) {
- triplets.push_back(Eigen::Triplet<double>(it.row() + rowOffset, it.col() + colOffset, it.value()));
+ int rowOffset = 0;
+ for (const auto& row : matrices) {
+ int colOffset = 0;
+ for (const auto& matrix : row) {
+ // Add the triplets of the matrix to the list of triplets for the larger matrix
+ for (int k = 0; k < matrix->outerSize(); ++k) {
+ for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(*matrix, k); it; ++it) {
+ triplets.push_back(Eigen::Triplet<double>(it.row() + rowOffset, it.col() + colOffset, it.value()));
+ }
}
+ colOffset += matrix->cols();
}
- colOffset += matrix->cols();
+ rowOffset += row[0]->rows();
}
- rowOffset += row[0]->rows();
-}
// Create a new matrix from the deque of triplets
matrixAdjoint.setFromTriplets(triplets.begin(), triplets.end());
matrixAdjoint.prune(0.);
- matrixAdjoint.makeCompressed();
+ //matrixAdjoint.makeCompressed();
}
void CoupledAdjoint::run()
@@ -159,64 +161,70 @@ void CoupledAdjoint::run()
/******** CL ********/
// Build right hand side for Cl
std::vector<double> rhsCl(A_adjoint.cols(), 0.0);
- std::copy(dCl_phi.data(), dCl_phi.data() + dCl_phi.size(), rhsCl.begin());
+ for (auto i = 0; i<dCl_phi.size(); ++i)
+ rhsCl[i] = dCl_phi[i];
+ // std::copy(dCl_phi.data(), dCl_phi.data() + dCl_phi.size(), rhsCl.begin());
+ Eigen::VectorXd rhsCl_ = Eigen::Map<const Eigen::VectorXd>(rhsCl.data(), rhsCl.size());
// Solution vector for lambdasCl
- Eigen::VectorXd lambdasCl = Eigen::VectorXd::Zero(A_adjoint.cols());
- Eigen::Map<Eigen::VectorXd> lambdasCl_(lambdasCl.data(), lambdasCl.size());
+ Eigen::VectorXd lambdaCl(A_adjoint.cols());
+ Eigen::Map<Eigen::VectorXd> lambdaCl_(lambdaCl.data(), lambdaCl.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 << "A_adjoint " << A_adjoint.rows() << " " << A_adjoint.cols() << " " << std::setprecision(12) << A_adjoint.norm() << std::endl;
+ // std::cout << "rhsCl " << std::setprecision(12) << Eigen::Map<Eigen::VectorXd>(rhsCl.data(), rhsCl.size()).norm() << std::endl;
+ // std::cout << "lambdasCl_ " << std::setprecision(12) << 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 " << std::setprecision(20) << A_adjoint.norm() << std::endl;
+ // // std::cout << "A_adjoint " << std::setprecision(3) << A_adjoint << std::endl;
+ std::cout << "rhsCl " << std::setprecision(20) << Eigen::Map<Eigen::VectorXd>(rhsCl.data(), rhsCl.size()).norm() << std::endl;
// for (auto i=0; i<rhsCl.size(); i++)
// {
- // std::cout << i << " " << rhsCl[i] << "\n";
+ // std::cout << i << " " << rhsCl[i] << "\n";
// }
- alinsol->compute(A_adjoint, Eigen::Map<Eigen::VectorXd>(rhsCl.data(), rhsCl.size()), lambdasCl_);
+ // solverTest.compute(A_adjoint);
+ // lambdaCl_ = solverTest.solve(Eigen::Map<Eigen::VectorXd>(rhsCl_.data(), rhsCl_.size()));
+
+ //std::cout << *alinsol << std::endl;
+ alinsol->compute(A_adjoint, Eigen::Map<Eigen::VectorXd>(rhsCl_.data(), rhsCl_.size()), lambdaCl_);
// 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 << "nIter " << alinsol->getIterations() << std::endl;
+ // std::cout << "error " << alinsol->getError() << std::endl;
- std::cout << "lambdasCl_ " << std::setprecision(10) << lambdasCl.norm() << std::endl;
+ std::cout << "lambdasCl_ " << std::setprecision(20) << lambdaCl.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;
+ // std::cout << "complete " << std::setprecision(12) << "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;
+ // std::cout << "lambdasCl " << std::setprecision(12) << lambdasCl.norm() << std::endl;
- Eigen::VectorXd lambdaClPhi = lambdasCl.block(0, 0, isol->pbl->msh->nodes.size(), 1);
+ Eigen::VectorXd lambdaClPhi = lambdaCl.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(12) << lambdaClPhi.norm() << 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;
+ // std::cout << "dCl_AoA (partial) " << std::setprecision(12) << dCl_AoA << std::endl;
+ // std::cout << "dRphi_AoA " << std::setprecision(12) << dRphi_AoA.transpose().norm() << std::endl;
+ // std::cout << "pdCl_AoA " << std::setprecision(12) << dCl_AoA << std::endl;
+ // std::cout << "dRphi_AoA " << std::setprecision(12) << dRphi_AoA.transpose().norm() << std::endl;
+ // std::cout << "lambdaClPhi " << std::setprecision(12) << lambdaClPhi.norm() << std::endl;
tdCl_AoA = dCl_AoA - dRphi_AoA.transpose()*lambdaClPhi;
- std::cout << "tdCl_AoA " << std::setprecision(10) << tdCl_AoA << std::endl;
- throw;
+ std::cout << "tdCl_AoA " << std::setprecision(12) << tdCl_AoA << std::endl;
/******** CD ********/
// Build right hand side for Cd
@@ -272,10 +280,10 @@ void CoupledAdjoint::gradientswrtInviscidFlow()
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));
+ 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, (isol->U[srfNode->row](i) - dv(i))/deltaPhi));
+ tripletsdv.push_back(Eigen::Triplet<double>(jj*pbl->nDim + i, n->row, (dv(i) - isol->U[srfNode->row](i))/deltaPhi));
jj++;
}
// Reset phi
@@ -442,35 +450,35 @@ void CoupledAdjoint::transposeMatrices(){
// // (Cout all norms)
- 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;
+ std::cout << "norm(dRphi_phi) " << std::setprecision(12) << dRphi_phi.norm() << std::endl;
+ std::cout << "norm(dRphi_M) " << std::setprecision(12) << dRphi_M.norm() << std::endl;
+ std::cout << "norm(dRphi_rho) " << std::setprecision(12) << dRphi_rho.norm() << std::endl;
+ std::cout << "norm(dRphi_v) " << std::setprecision(12) << dRphi_v.norm() << std::endl;
+ std::cout << "norm(dRphi_blw) " << std::setprecision(12) << dRphi_blw.norm() << std::endl;
+
+ std::cout << "norm(dRM_phi) " << std::setprecision(12) << dRM_phi.norm() << std::endl;
+ std::cout << "norm(dRM_M) " << std::setprecision(12) << dRM_M.norm() << std::endl;
+ std::cout << "norm(dRM_rho) " << std::setprecision(12) << dRM_rho.norm() << std::endl;
+ std::cout << "norm(dRM_v) " << std::setprecision(12) << dRM_v.norm() << std::endl;
+ std::cout << "norm(dRM_blw) " << std::setprecision(12) << dRM_blw.norm() << std::endl;
+
+ std::cout << "norm(dRrho_phi) " << std::setprecision(12) << dRrho_phi.norm() << std::endl;
+ std::cout << "norm(dRrho_M) " << std::setprecision(12) << dRrho_M.norm() << std::endl;
+ std::cout << "norm(dRrho_rho) " << std::setprecision(12) << dRrho_rho.norm() << std::endl;
+ std::cout << "norm(dRrho_v) " << std::setprecision(12) << dRrho_v.norm() << std::endl;
+ std::cout << "norm(dRrho_blw) " << std::setprecision(12) << dRrho_blw.norm() << std::endl;
+
+ std::cout << "norm(dRv_phi) " << std::setprecision(12) << dRv_phi.norm() << std::endl;
+ std::cout << "norm(dRv_M) " << std::setprecision(12) << dRv_M.norm() << std::endl;
+ std::cout << "norm(dRv_rho) " << std::setprecision(12) << dRv_rho.norm() << std::endl;
+ std::cout << "norm(dRv_v) " << std::setprecision(12) << dRv_v.norm() << std::endl;
+ std::cout << "norm(dRv_blw) " << std::setprecision(12) << dRv_blw.norm() << std::endl;
+
+ std::cout << "norm(dRblw_phi) " << std::setprecision(12) << dRblw_phi.norm() << std::endl;
+ std::cout << "norm(dRblw_M) " << std::setprecision(20) << dRblw_M.norm() << std::endl;
+ std::cout << "norm(dRblw_rho) " << std::setprecision(20) << dRblw_rho.norm() << std::endl;
+ std::cout << "norm(dRblw_v) " << std::setprecision(20) << dRblw_v.norm() << std::endl;
+ std::cout << "norm(dRblw_blw) " << std::setprecision(12) << dRblw_blw.norm() << std::endl;
// // Print all matrices after transpose
@@ -504,31 +512,31 @@ void CoupledAdjoint::transposeMatrices(){
// 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;
+ // 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) {
@@ -536,7 +544,7 @@ void CoupledAdjoint::writeMatrixToFile(const Eigen::MatrixXd& matrix, const std:
if (file.is_open()) {
for (int i = 0; i < matrix.rows(); ++i) {
for (int j = 0; j < matrix.cols(); ++j) {
- file << matrix(i, j);
+ file << std::setprecision(20) << matrix(i, j);
if (j != matrix.cols() - 1) {
file << " ";
}
diff --git a/blast/src/DCoupledAdjoint.h b/blast/src/DCoupledAdjoint.h
index 8b2a61e..5f69e3a 100644
--- a/blast/src/DCoupledAdjoint.h
+++ b/blast/src/DCoupledAdjoint.h
@@ -23,6 +23,7 @@
#include <Eigen/Sparse>
#include <iostream>
+#include <memory>
namespace blast
{
diff --git a/blast/src/DDriver.cpp b/blast/src/DDriver.cpp
index d41f579..0ccb893 100644
--- a/blast/src/DDriver.cpp
+++ b/blast/src/DDriver.cpp
@@ -306,7 +306,7 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, int forced)
double avLam = 1. - avTurb;
// Impose boundary condition.
- double Cteq_trans;
+ double Cteq_trans = 0.;
bl->closSolver->turbulentClosures(Cteq_trans, bl->u[(iPoint - 1) * nVar + 0], bl->u[(iPoint - 1) * nVar + 1], 0.03, bl->Ret[iPoint - 1], bl->blEdge->getMe(iPoint - 1), bl->name);
bl->ctEq[iPoint - 1] = Cteq_trans;
bl->u[(iPoint - 1) * nVar + 4] = 0.7 * bl->ctEq[iPoint - 1];
diff --git a/blast/tests/dart/adjointVII.py b/blast/tests/dart/adjointVII.py
index 8f1c6d3..5d10de6 100644
--- a/blast/tests/dart/adjointVII.py
+++ b/blast/tests/dart/adjointVII.py
@@ -66,7 +66,7 @@ def cfgInviscid(nthrds, verb):
'Wake' : 'wake', # LIST of names of physical group containing the wake
'WakeTip' : 'wakeTip', # LIST of names of physical group containing the edge of the wake
'Te' : 'te', # LIST of names of physical group containing the trailing edge
- 'dbc' : True,
+ 'dbc' : False,
'Upstream' : 'upstream',
# Freestream
'M_inf' : 0.2, # freestream Mach number
@@ -78,7 +78,7 @@ def cfgInviscid(nthrds, verb):
'y_ref' : 0.0, # reference point for moment computation (y)
'z_ref' : 0.0, # reference point for moment computation (z)
# Numerical
- 'LSolver' : 'GMRES', # inner solver (Pardiso, MUMPS or GMRES)
+ 'LSolver' : 'PARDISO', # inner solver (Pardiso, MUMPS or GMRES)
'G_fill' : 1, # fill-in factor for GMRES preconditioner
'G_tol' : 1e-5, # tolerance for GMRES
'G_restart' : 50, # restart for GMRES
@@ -93,7 +93,7 @@ def cfgBlast(verb):
'Minf' : 0.2, # Freestream Mach number (used for the computation of the time step only)
'CFL0' : 1, # Inital CFL number of the calculation
'Verb': verb, # Verbosity level of the solver
- 'couplIter': 25, # Maximum number of iterations
+ 'couplIter': 75, # Maximum number of iterations
'couplTol' : 1e-4, # Tolerance of the VII methodology
'iterPrint': 5, # int, number of iterations between outputs
'resetInv' : True, # bool, flag to reset the inviscid calculation at every iteration.
@@ -135,11 +135,11 @@ def main():
tms['solver'].start()
aeroCoeffs = coupler.run()
tms['solver'].stop()
- #iSolverAPI.run()
- print('Dart adjoint')
+ tms['TotalAdjoint'].start()
+ tms['dartAdjoint'].start()
iSolverAPI.adjointSolver.run()
- tms['CoupledAdjoint'].stop()
+ tms['dartAdjoint'].stop()
dCl_dAlpha_DartAdjoint = iSolverAPI.adjointSolver.dClAoa
dCd_dAlpha_DartAdjoint = iSolverAPI.adjointSolver.dCdAoa
print('Compute derivative')
@@ -147,10 +147,28 @@ def main():
pyAdjoint.computeDerivatives()
tms['FDDerivatives'].stop()
tms['CoupledAdjoint'].start()
- print('Coupled adjoint')
- tms['CouledAdjoint'].stop()
coupledAdjointSolver.run()
- print(coupledAdjointSolver.tdCl_AoA)
+ tms['CoupledAdjoint'].stop()
+ tms['TotalAdjoint'].stop()
+
+ # norm1 = []
+ # norm2 = []
+ # norm3 = []
+
+ # for i in range(1):
+ # pyAdjoint.computeDerivatives()
+ # quit()
+
+
+ # clv = []
+ # for i in range(50):
+ # coupledAdjointSolver.reset()
+ # coupledAdjointSolver.run()
+ # clv.append(coupledAdjointSolver.tdCl_AoA)
+ # from matplotlib import pyplot as plt
+ # print(clv)
+ # print(coupledAdjointSolver.tdCl_AoA)
+ # quit()
# quit()
@@ -307,12 +325,16 @@ def main():
# print(i)
# print("sucess", coupledAdjointSolver)
# quit()
+
+ #icfg['dbc'] = True
dalpha = 1*np.pi/180
aoaV = [alpha-dalpha, alpha+dalpha]
clV = []
cdV = []
+ tms['coupledFD'].start()
+
for aoa in aoaV:
coupler, iSolverAPI, vSolver = viscUtils.initBlast(icfg, vcfg)
iSolverAPI.argOut['pbl'].update(aoa)
@@ -321,10 +343,12 @@ def main():
cdV.append(iSolverAPI.getCd())
dCl_dalpha_FD_coupled = (clV[-1] - clV[0]) / (2*dalpha)
dCd_dalpha_FD_coupled = (cdV[-1] - cdV[0]) / (2*dalpha)
+ tms['coupledFD'].stop()
aoaV = [alpha-dalpha, alpha+dalpha]
clV = []
cdV = []
+ tms['dartFD'].start()
for aoa in aoaV:
coupler, iSolverAPI, vSolver = viscUtils.initBlast(icfg, vcfg)
iSolverAPI.argOut['pbl'].update(aoa)
@@ -333,6 +357,7 @@ def main():
cdV.append(iSolverAPI.getCd())
dCl_dalpha_FD_dart = (clV[-1] - clV[0]) / (2*dalpha)
dCd_dalpha_FD_dart = (cdV[-1] - cdV[0]) / (2*dalpha)
+ tms['dartFD'].stop()
print(' ------------------- RESULTS -------------------')
print('********************** CL **********************')
--
GitLab