From 3bf01e74ece0b90b2fabc9c1e4beee0819c73fdc Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Tue, 7 May 2024 17:15:24 +0200
Subject: [PATCH] (feat) Added adjoint
Adjoint is not working because the delta star box is missing
---
blast/_src/blastw.h | 3 +-
blast/_src/blastw.i | 10 +-
blast/coupler.py | 56 ++-
blast/interfaces/DAdjointInterface.py | 102 ++++++
blast/interfaces/DDataStructure.py | 1 -
blast/interfaces/DSolversInterface.py | 29 +-
blast/interfaces/dart/DartInterface.py | 31 +-
blast/src/CMakeLists.txt | 7 +-
blast/src/DBoundaryLayer.cpp | 168 ++++++++-
blast/src/DBoundaryLayer.h | 6 +
blast/src/DCoupledAdjoint.cpp | 460 +++++++++++++++++++++++++
blast/src/DCoupledAdjoint.h | 128 +++++++
blast/src/DDriver.cpp | 25 +-
blast/src/DDriver.h | 3 +
blast/src/blast.h | 3 +
blast/tests/dart/adjointVII_2D.py | 174 ++++++++++
blast/tests/dart/naca0012_2D.py | 22 +-
17 files changed, 1168 insertions(+), 60 deletions(-)
create mode 100644 blast/interfaces/DAdjointInterface.py
create mode 100644 blast/src/DCoupledAdjoint.cpp
create mode 100644 blast/src/DCoupledAdjoint.h
create mode 100644 blast/tests/dart/adjointVII_2D.py
diff --git a/blast/_src/blastw.h b/blast/_src/blastw.h
index f0b51a2..dd883e5 100644
--- a/blast/_src/blastw.h
+++ b/blast/_src/blastw.h
@@ -15,4 +15,5 @@
*/
#include "DDriver.h"
-#include "DBoundaryLayer.h"
\ No newline at end of file
+#include "DBoundaryLayer.h"
+#include "DCoupledAdjoint.h"
\ No newline at end of file
diff --git a/blast/_src/blastw.i b/blast/_src/blastw.i
index a28e208..3528a3b 100644
--- a/blast/_src/blastw.i
+++ b/blast/_src/blastw.i
@@ -43,11 +43,10 @@ threads="1"
// ----------- blast CLASSES ----------------
%include "blast.h"
-%shared_ptr(blast::Driver);
%shared_ptr(blast::BoundaryLayer);
%include "DBoundaryLayer.h"
-
+%shared_ptr(blast::Driver);
%immutable blast::Driver::Re; // read-only variable (no setter)
%immutable blast::Driver::Cdt;
%immutable blast::Driver::Cdt_sec;
@@ -57,4 +56,9 @@ threads="1"
%immutable blast::Driver::Cdp_sec;
%immutable blast::Driver::CFL0;
%immutable blast::Driver::verbose;
-%include "DDriver.h"
\ No newline at end of file
+%include "DDriver.h"
+
+%shared_ptr(blast::CoupledAdjoint);
+%immutable blast::CoupledAdjoint::tdCl_AoA;
+%immutable blast::CoupledAdjoint::tdCd_AoA;
+%include "DCoupledAdjoint.h"
diff --git a/blast/coupler.py b/blast/coupler.py
index 91ecd93..b7213a4 100644
--- a/blast/coupler.py
+++ b/blast/coupler.py
@@ -33,8 +33,6 @@ class Coupler:
self.iterPrint = _iterPrint if self.isol.getVerbose() == 0 and self.vsol.verbose == 0 else 1
self.tms = fwk.Timers()
self.filesfx = sfx
-
- def run(self):
print('')
print(' ######################################################## ')
print('| ___ __ ___ _________________ |')
@@ -54,13 +52,14 @@ class Coupler:
print('Interpolator:', self.isol.cfg['interpolator'])
print('')
print('-- Coupler parameters --')
- print('MaxIter:{:>3.0f}'.format(self.maxIter))
- print('Tolerance:{:>6.2f}'.format(np.log10(self.tol)))
- print('IterPrint:{:>2.0f}'.format(self.iterPrint))
- print('ResetInv:{:>5s}'.format(str(self.resetInviscid)))
- print('Verb:{:>2.0f}'.format(self.isol.getVerbose()))
+ print('MaxIter: {:>3.0f}'.format(self.maxIter))
+ print('Tolerance: {:>6.2f}'.format(np.log10(self.tol)))
+ print('IterPrint: {:>2.0f}'.format(self.iterPrint))
+ print('ResetInv: {:>5s}'.format(str(self.resetInviscid)))
+ print('Verb: {:>2.0f}'.format(self.isol.getVerbose()))
print('')
+ def run(self, write=True):
# Aerodynamic coefficients.
aeroCoeffs = {'Cl':[], 'Cd': [], 'Cdwake': []}
@@ -72,6 +71,7 @@ class Coupler:
# Impose blowing boundary condition in the inviscid solver.
self.tms['processing'].start()
self.isol.getBlowingBoundary(couplIter)
+ self.isol.interpolate('v2i')
self.isol.setBlowingVelocity()
self.tms['processing'].stop()
@@ -85,11 +85,13 @@ class Coupler:
# Write inviscid Cp distribution.
if couplIter == 0:
self.isol.initializeViscousSolver()
- self.isol.writeCp(sfx='_inviscid'+self.filesfx)
+ if write:
+ self.isol.writeCp(sfx='_inviscid'+self.filesfx)
# Impose inviscid boundary in the viscous solver.
self.tms['processing'].start()
self.isol.getInviscidBC()
+ self.isol.interpolate('i2v')
self.isol.setViscousSolver()
self.tms['processing'].stop()
@@ -109,7 +111,8 @@ class Coupler:
if error <= self.tol:
print(ccolors.ANSI_GREEN, '{:>4.0f}| {:>7.5f} {:>7.5f} {:>7.5f} | {:>6.4f} {:>6.4f} | {:>6.3f}\n'.format(couplIter, self.isol.getCl(), self.isol.getCd()+self.vsol.Cdf, self.vsol.Cdt, self.vsol.getAverageTransition(0), self.vsol.getAverageTransition(1), np.log10(error)), ccolors.ANSI_RESET)
- self.isol.writeCp(sfx='_viscous'+self.filesfx)
+ if write:
+ self.isol.writeCp(sfx='_viscous'+self.filesfx)
return aeroCoeffs
cdPrev = cd
@@ -128,5 +131,38 @@ class Coupler:
print('{:>5s}| {:>7s} {:>7s} {:>7s} | {:>6s} {:>8s} | {:>6s}'.format('It', 'Cl', 'Cd', 'Cdwake', 'xtrT', 'xtrB', 'Error'))
print(' ----------------------------------------------------------')
print(ccolors.ANSI_RED, '{:>4.0f}| {:>7.5f} {:>7.5f} {:>7.5f} | {:>6.4f} {:>7.4f} | {:>6.3f}\n'.format(couplIter-1, self.isol.getCl(), self.isol.getCd()+self.vsol.Cdf, self.vsol.Cdt, self.vsol.getAverageTransition(0), self.vsol.getAverageTransition(1), np.log10(error)), ccolors.ANSI_RESET)
- self.isol.writeCp(sfx='_viscous'+self.filesfx)
+ if write:
+ self.isol.writeCp(sfx='_viscous'+self.filesfx)
return aeroCoeffs
+
+ def reset(self):
+ """Reset the entire VII coupling; inviscid solver, transfer quantities
+ and blowing velocity. Viscous section objects are not deleted
+ """
+ # Inviscid solver
+
+ self.isol.reset()
+
+ for iSec in range(self.isol.cfg['nSections']):
+ for i, reg in enumerate(self.isol.sec[iSec]):
+ if reg.name == 'wake':
+ iReg = 1
+ elif reg.name == 'lower' or reg.name == 'upper':
+ iReg = 0
+ else:
+ raise RuntimeError('Invalid region', reg.name)
+ loc = np.zeros(reg.getnNodes())
+ for iPoint in range(reg.getnNodes()):
+ loc[iPoint] = reg.loc[iPoint]
+ self.isol.xxExt[iSec][i] = loc
+ self.isol.deltaStarExt[iSec][i] = np.zeros(reg.getnNodes())
+
+ # Blowing velocity
+ for b in self.isol.iBnd:
+ for i in range(len(b.blowingVel)):
+ b.blowingVel[i] = 0.
+ for sec in self.isol.vBnd:
+ for side in sec:
+ for i in range(len(side.blowingVel)):
+ side.blowingVel[i] = 0.
+ self.isol.setBlowingVelocity()
diff --git a/blast/interfaces/DAdjointInterface.py b/blast/interfaces/DAdjointInterface.py
new file mode 100644
index 0000000..cf3f60f
--- /dev/null
+++ b/blast/interfaces/DAdjointInterface.py
@@ -0,0 +1,102 @@
+import numpy as np
+from matplotlib import pyplot as plt
+
+class AdjointInterface():
+ def __init__(self, _coupledAdjointSolver, _iSolverAPI, _vSolver):
+ self.coupledAdjointSolver = _coupledAdjointSolver
+ self.isol = _iSolverAPI
+ self.vsol = _vSolver
+
+ def build(self):
+ nDim = self.isol.solver.pbl.nDim
+ nNd = self.isol.solver.pbl.msh.nodes.size()
+ nEs = len(self.isol.iBnd[0].elemsCoord[:,0])
+ nNs = len(self.isol.blw[0].nodes)
+ dRblw_M = np.zeros((nEs, nNs))
+ dRblw_rho = np.zeros((nEs, nNs))
+ dRblw_v = np.zeros((nEs, nDim*nNs))
+ dRblw_x = np.zeros((nEs, nDim*nNd))
+
+ delta = 1e-6
+
+ saveBlw = np.zeros(nEs)
+ for i in range(nEs):
+ saveBlw[i] = self.isol.blw[0].getU(i)
+
+ # Mach
+ print('Mach', end='...')
+ for i, n in enumerate(self.isol.blw[0].nodes):
+ savemach = self.isol.solver.M[n.row]
+ self.isol.solver.M[n.row] = savemach + delta
+ blowing1 = self.__runViscous()
+
+ self.isol.solver.M[n.row] = savemach - delta
+ blowing2 = self.__runViscous()
+ self.isol.solver.M[n.row] = savemach
+ dRblw_M[:,i] = -(blowing1 - blowing2)/(2*delta)
+ print('done')
+
+ # Rho
+ print('Density', end='...')
+ for i, n in enumerate(self.isol.blw[0].nodes):
+ saverho = self.isol.solver.rho[n.row]
+ self.isol.solver.rho[n.row] = saverho + delta
+ blowing1 = self.__runViscous()
+
+ self.isol.solver.rho[n.row] = saverho - delta
+ blowing2 = self.__runViscous()
+ self.isol.solver.rho[n.row] = saverho
+ dRblw_rho[:,i] = -(blowing1 - blowing2)/(2*delta)
+ print('done')
+
+ # # V
+ print('Velocity', end='...')
+ for i, n in enumerate(self.isol.blw[0].nodes):
+ for jj in range(nDim):
+ savev = self.isol.solver.U[n.row][jj]
+ self.isol.solver.U[n.row][jj] = savev + delta
+ blowing1 = self.__runViscous()
+
+ self.isol.solver.U[n.row][jj] = savev - delta
+ blowing2 = self.__runViscous()
+ self.isol.solver.U[n.row][jj] = savev
+ dRblw_v[:,i*nDim + jj] = -(blowing1 - blowing2)/(2*delta)
+ print('done')
+
+ """# # Mesh coordinates
+ print('Mesh')
+ for i, n in enumerate(self.iSolverAPI.blw[0].nodes):
+ vidx = np.where(self.iSolverAPI.iBnd[0].nodesCoord[:,3] == n.row)[0][0]
+ for jj in range(nDim):
+ if self.iSolverAPI.iBnd[0].nodesCoord[vidx,jj] != n.pos[jj]:
+ print('WRONG NODE MESH')
+ quit()
+ savemesh = n.pos[jj]
+ self.iSolverAPI.iBnd[0].nodesCoord[vidx, jj] = savemesh + delta
+ blowing1 = self.__runViscous()
+
+ self.iSolverAPI.iBnd[0].nodesCoord[vidx, jj] = savemesh - delta
+ blowing2 = self.__runViscous()
+ self.iSolverAPI.iBnd[0].nodesCoord[vidx, jj] = savemesh
+ dRblw_x[:, n.row*nDim+jj] = -(blowing1 - blowing2)/(2*delta)"""
+
+ self.coupledAdjointSolver.setdRblw_M(dRblw_M)
+ self.coupledAdjointSolver.setdRblw_rho(dRblw_rho)
+ self.coupledAdjointSolver.setdRblw_v(dRblw_v)
+ self.coupledAdjointSolver.setdRblw_x(dRblw_x)
+
+ def __runViscous(self):
+ self.isol.getInviscidBC()
+ self.isol.interpolate('i2v')
+ self.isol.setViscousSolver(adjoint=True)
+ ext = self.vsol.run()
+ if ext != 0:
+ raise ValueError('Viscous solver did not converge')
+ self.isol.getBlowingBoundary(1, adjoint=True)
+ self.isol.interpolate('v2i')
+ self.isol.setBlowingVelocity()
+
+ blowing = np.zeros(len(self.isol.iBnd[0].elemsCoord[:,0]))
+ for i in range(len(self.isol.iBnd[0].elemsCoord[:,0])):
+ blowing[i] = self.isol.blw[0].getU(i)
+ return blowing
diff --git a/blast/interfaces/DDataStructure.py b/blast/interfaces/DDataStructure.py
index 8fa91c0..4a63dc6 100644
--- a/blast/interfaces/DDataStructure.py
+++ b/blast/interfaces/DDataStructure.py
@@ -25,7 +25,6 @@ class Group:
self.nPoints = nPoints
self.nElems = nPoints - 1
self.nodesCoord = np.zeros((self.nPoints,4))
- self.elemsCoord = np.zeros((self.nPoints,4))
self.V = np.zeros((self.nPoints,3))
self.M = np.zeros(self.nPoints)
diff --git a/blast/interfaces/DSolversInterface.py b/blast/interfaces/DSolversInterface.py
index c1f9ed1..dc85a47 100644
--- a/blast/interfaces/DSolversInterface.py
+++ b/blast/interfaces/DSolversInterface.py
@@ -25,6 +25,7 @@ class SolversInterface:
for iReg, reg in enumerate(self.vBnd[iSec]):
reg.initStructures(self.iBnd[iReg].nPoints)
self.vBnd[iSec][iReg].nodesCoord = self.iBnd[iReg].nodesCoord
+ self.vBnd[iSec][iReg].elemsCoord = self.iBnd[iReg].elemsCoord
elif self.cfg['interpolator'] == 'Rbf':
from blast.interfaces.interpolators.DRbfInterpolator import RbfInterpolator as interp
@@ -44,12 +45,11 @@ class SolversInterface:
for iSec in range(self.cfg['nSections'])]
self.vtExt = [[np.zeros(0) for iReg in range(3)]\
for iSec in range(self.cfg['nSections'])]
+ self.vinit = False
def setViscousSolver(self, adjoint=False):
"""Interpolate solution to viscous mesh and sets the inviscid boundary in the viscous solver.
"""
- self.interpolator.inviscidToViscous(self.iBnd, self.vBnd)
-
for iSec in range(self.cfg['nSections']):
for i, reg in enumerate(self.sec[iSec]):
if reg.name == 'wake':
@@ -67,8 +67,8 @@ class SolversInterface:
reg.setDeltaStarExt(self.deltaStarExt[iSec][i])
reg.setxxExt(self.xxExt[iSec][i])
- def getBlowingBoundary(self, couplIter):
- """ Get blowing boundary condition from the viscous solver.
+ def getBlowingBoundary(self, couplIter, adjoint=False):
+ """Get blowing boundary condition from the viscous solver.
args:
----
couplIter: int
@@ -82,15 +82,26 @@ class SolversInterface:
elif 'vAirfoil' in reg.name:
reg.blowingVel = np.concatenate((np.flip(np.asarray(self.sec[iSec][0].blowingVelocity)),
np.asarray(self.sec[iSec][1].blowingVelocity)))
- for i in range(len(self.sec[iSec])):
- self.xxExt[iSec][i] = np.asarray(self.sec[iSec][i].loc)
- self.deltaStarExt[iSec][i] = np.asarray(self.sec[iSec][i].deltaStar)
- self.interpolator.viscousToInviscid(self.iBnd, self.vBnd)
+ if adjoint == False:
+ for i in range(len(self.sec[iSec])):
+ self.xxExt[iSec][i] = np.asarray(self.sec[iSec][i].loc)
+ self.deltaStarExt[iSec][i] = np.asarray(self.sec[iSec][i].deltaStar)
+
+ def interpolate(self, dir):
+ if dir == 'i2v':
+ self.interpolator.inviscidToViscous(self.iBnd, self.vBnd)
+ elif dir == 'v2i':
+ self.interpolator.viscousToInviscid(self.iBnd, self.vBnd)
+ else:
+ raise RuntimeError('Invalid direction', dir)
def initializeViscousSolver(self):
"""Initialize the viscous sections in the viscous solver.
"""
+ if self.vinit == True:
+ return None
+
self.sec = [[] for _ in range(self.cfg['nSections'])]
self.vSolver.reset()
for i in range(self.cfg['nSections']):
@@ -161,6 +172,7 @@ class SolversInterface:
loc[iPoint] = reg.loc[iPoint]
self.xxExt[iSec][i] = loc
self.deltaStarExt[iSec][i] = np.zeros(reg.getnNodes())
+ self.vinit = True
def getVWing(self):
"""Obtain the nodes' location and row and cg of all elements.
@@ -169,7 +181,6 @@ class SolversInterface:
If not, the nodes are pseudo-duplicated (without changing the row). They therefore
exist only for the interpolation but not for the solution.
"""
-
data = [np.zeros((0,4)) for _ in range(2)]
lwRows = []
for e in self.cfg['vMsh'].elems:
diff --git a/blast/interfaces/dart/DartInterface.py b/blast/interfaces/dart/DartInterface.py
index 9c6f040..c505ec0 100644
--- a/blast/interfaces/dart/DartInterface.py
+++ b/blast/interfaces/dart/DartInterface.py
@@ -28,10 +28,14 @@ class DartInterface(SolversInterface):
def __init__(self, dartCfg, vSolver, _cfg):
try:
from modules.dartflo.dart.api.core import initDart
- self.inviscidObjects = initDart(dartCfg, viscous=True)
+ self.inviscidObjects = initDart(dartCfg, task=dartCfg['task'], viscous=True)
self.solver = self.inviscidObjects.get('sol') # Solver
self.blw = [self.inviscidObjects.get('blwb'), self.inviscidObjects.get('blww')]
- print('Dart successfully loaded.')
+ if dartCfg['task'] == 'optimization':
+ self.adjointSolver = self.inviscidObjects.get('adj')
+ print('Dart successfully loaded for optimization.')
+ else:
+ print('Dart successfully loaded for analysis.')
except:
print(ccolors.ANSI_RED, 'Could not load DART. Make sure it is installed.', ccolors.ANSI_RESET)
raise RuntimeError('Import error')
@@ -95,6 +99,9 @@ class DartInterface(SolversInterface):
def getAoA(self):
return self.solver.pbl.alpha
+ def setAoA(self, aoa):
+ self.solver.pbl.update(aoa)
+
def getMinf(self):
return self.solver.pbl.M_inf
@@ -113,6 +120,9 @@ class DartInterface(SolversInterface):
def getnDim(self):
return self.inviscidObjects['pbl'].nDim
+ def getnNodesDomain(self):
+ return self.solver.pbl.msh.nodes.size()
+
def reset(self):
self.solver.reset()
@@ -137,7 +147,6 @@ class DartInterface(SolversInterface):
""" Extract the solution of the viscous calculation (blowing velocity)
and use it as a boundary condition in the inviscid solver.
"""
-
# Impose blowing velocities.
for n in range(len(self.iBnd)):
if self.getnDim() == 2:
@@ -260,6 +269,14 @@ class DartInterface(SolversInterface):
self.iBnd[0].nodesCoord = np.column_stack((data[:,1], data[:,2],\
data[:,3], data[:,4]))
self.iBnd[0].setConnectList(connectListNodes, connectListElems)
+ elemCoord = np.zeros((len(self.iBnd[0].nodesCoord)-1, 4))
+ for i in range(len(elemCoord[:,0])):
+ elemCoord[i,0] = 0.5 * (self.iBnd[0].nodesCoord[i,0] + self.iBnd[0].nodesCoord[i+1,0])
+ elemCoord[i,1] = 0.5 * (self.iBnd[0].nodesCoord[i,1] + self.iBnd[0].nodesCoord[i+1,1])
+ elemCoord[i,2] = 0.5 * (self.iBnd[0].nodesCoord[i,2] + self.iBnd[0].nodesCoord[i+1,2])
+ elemCoord[i,3] = i
+ self.iBnd[0].elemsCoord = elemCoord
+
# Wake
self.iBnd[1].initStructures(self.blw[1].nodes.size())
@@ -299,6 +316,14 @@ class DartInterface(SolversInterface):
dataW[:,3], dataW[:,4]))
self.iBnd[1].setConnectList(connectListNodes, connectListElems)
+ elemCoordW = np.zeros((len(self.iBnd[1].nodesCoord)-1, 4))
+ for i in range(len(elemCoordW[:,0])):
+ elemCoordW[i,0] = 0.5 * (self.iBnd[1].nodesCoord[i,0] + self.iBnd[1].nodesCoord[i+1,0])
+ elemCoordW[i,1] = 0.5 * (self.iBnd[1].nodesCoord[i,1] + self.iBnd[1].nodesCoord[i+1,1])
+ elemCoordW[i,2] = 0.5 * (self.iBnd[1].nodesCoord[i,2] + self.iBnd[1].nodesCoord[i+1,2])
+ elemCoordW[i,3] = i
+ self.iBnd[1].elemsCoord = elemCoordW
+
def getWing(self):
"""Obtain the nodes' location and row and cg of all elements.
"""
diff --git a/blast/src/CMakeLists.txt b/blast/src/CMakeLists.txt
index 96caacb..cd3f09b 100644
--- a/blast/src/CMakeLists.txt
+++ b/blast/src/CMakeLists.txt
@@ -18,8 +18,11 @@ FILE(GLOB SRCS *.h *.cpp *.inl *.hpp)
ADD_LIBRARY(blast SHARED ${SRCS})
MACRO_DebugPostfix(blast)
-TARGET_INCLUDE_DIRECTORIES(blast PUBLIC ${PROJECT_SOURCE_DIR}/blast/src)
-
+TARGET_INCLUDE_DIRECTORIES(blast PUBLIC ${PROJECT_SOURCE_DIR}/blast/src
+ ${PROJECT_SOURCE_DIR}/modules/dartflo/ext/amfe/tbox/src
+ ${PROJECT_SOURCE_DIR}/modules/dartflo/ext/amfe/fwk/src
+ ${PROJECT_SOURCE_DIR}/modules/dartflo/dart/src
+ ${PYTHON_INCLUDE_PATH})
TARGET_LINK_LIBRARIES(blast dart fwk tbox)
INSTALL(TARGETS blast DESTINATION ${CMAKE_INSTALL_PREFIX})
diff --git a/blast/src/DBoundaryLayer.cpp b/blast/src/DBoundaryLayer.cpp
index 86fc294..ad18c99 100644
--- a/blast/src/DBoundaryLayer.cpp
+++ b/blast/src/DBoundaryLayer.cpp
@@ -173,27 +173,169 @@ void BoundaryLayer::setWakeBC(double Re, std::vector<double> UpperCond, std::vec
*/
void BoundaryLayer::computeBlowingVelocity()
{
- deltaStar[0] = u[0] * u[1];
+ //deltaStar[0] = u[0] * u[1];
blowingVelocity.resize(nNodes - 1, 0.);
- for (size_t iPoint = 1; iPoint < nNodes; ++iPoint)
+ for (size_t i = 1; i < nNodes; ++i)
{
- deltaStar[iPoint] = u[iPoint * nVar + 0] * u[iPoint * nVar + 1];
- blowingVelocity[iPoint - 1] = (rhoe[iPoint] * vt[iPoint] * deltaStar[iPoint] - rhoe[iPoint-1] * vt[iPoint-1] * deltaStar[iPoint - 1]) / (rhoe[iPoint] * (loc[iPoint] - loc[iPoint-1]));
- if (std::isnan(blowingVelocity[iPoint - 1]))
+ //deltaStar[i] = u[i * nVar + 0] * u[i * nVar + 1];
+ blowingVelocity[i - 1] = (rhoe[i] * vt[i] * deltaStar[i] - rhoe[i-1] * vt[i-1] * deltaStar[i - 1]) / (rhoe[i] * (loc[i] - loc[i-1]));
+ if (std::isnan(blowingVelocity[i - 1]))
{
- if (rhoe[iPoint] == 0.)
- std::cout << "Density is zero at point " << iPoint << std::endl;
- if ((loc[iPoint] - loc[iPoint-1]) == 0.)
- std::cout << "Points " << iPoint -1 << " and " << iPoint << " are at the same position " << loc[iPoint] << ", resulting in an infinite gradient." << std::endl;
- if (std::isnan(deltaStar[iPoint]))
- std::cout << "Displacement thickness is nan at position " << iPoint << std::endl;
- if (std::isnan(deltaStar[iPoint - 1]))
- std::cout << "Displacement thickness is nan at position " << iPoint - 1 << std::endl;
+ if (rhoe[i] == 0.)
+ std::cout << "Density is zero at point " << i << std::endl;
+ if ((loc[i] - loc[i-1]) == 0.)
+ std::cout << "Points " << i -1 << " and " << i << " are at the same position " << loc[i] << ", resulting in an infinite gradient." << std::endl;
+ if (std::isnan(deltaStar[i]))
+ std::cout << "Displacement thickness is nan at position " << i << std::endl;
+ if (std::isnan(deltaStar[i - 1]))
+ std::cout << "Displacement thickness is nan at position " << i - 1 << std::endl;
throw std::runtime_error("NaN detected in blowing velocity computation (see reason above)");
}
}
}
+std::vector<std::vector<double>> BoundaryLayer::evalGradwrtDeltaStar()
+{
+ std::vector<std::vector<double>> gradVe = std::vector<std::vector<double>>(nElms);
+ for (auto &vec : gradVe)
+ vec.resize(nNodes, 0.);
+
+ for (size_t i = 1; i < nNodes; ++i)
+ {
+ gradVe[i-1][i] = vt[i] / (loc[i] - loc[i-1]);
+ gradVe[i-1][i-1] = -1/rhoe[i] * rhoe[i-1] * vt[i-1] / (loc[i] - loc[i-1]);
+ }
+ return gradVe;
+}
+
+std::vector<std::vector<double>> BoundaryLayer::evalGradwrtVt()
+{
+ std::vector<std::vector<double>> gradVe = std::vector<std::vector<double>>(nElms);
+ for (auto &vec : gradVe)
+ vec.resize(nNodes, 0.);
+
+ for (size_t i = 1; i < nNodes; ++i)
+ {
+ gradVe[i-1][i] = deltaStar[i] / (loc[i] - loc[i-1]);
+ gradVe[i-1][i - 1] = - rhoe[i-1] / rhoe[i] * deltaStar[i-1] / (loc[i] - loc[i-1]);
+ }
+ return gradVe;
+}
+
+std::vector<std::vector<double>> BoundaryLayer::evalGradwrtRho()
+{
+ std::vector<std::vector<double>> gradVe = std::vector<std::vector<double>>(nElms);
+ for (auto &vec : gradVe)
+ vec.resize(nNodes, 0.);
+
+ for (size_t i = 1; i < nNodes; ++i)
+ {
+ gradVe[i-1][i] = -1/(rhoe[i]*rhoe[i]) * ((rhoe[i] * vt[i] * deltaStar[i] - rhoe[i-1] * vt[i-1] * deltaStar[i-1])/(loc[i] - loc[i-1]))
+ + 1/rhoe[i] * vt[i] * deltaStar[i] / (loc[i] - loc[i-1]);
+ gradVe[i-1][i-1] = -vt[i-1] * deltaStar[i-1] / (rhoe[i] * (loc[i] - loc[i-1]));
+ }
+ return gradVe;
+}
+
+/**
+ * @brief Compute the gradient of the velocity with respect to the x coordinate.
+ *
+ * dVe_dx = dVe_dloc * dloc_dx
+ *
+ * @return std::vector<double> Gradient of the velocity with respect to the x coordinate.
+ */
+std::vector<std::vector<double>> BoundaryLayer::evalGradwrtX()
+{
+ std::vector<std::vector<double>> gradVe = std::vector<std::vector<double>>(nElms);
+ for (auto &vec : gradVe)
+ vec.resize(nNodes, 0.);
+
+ // Compute dVe/dloc
+ std::vector<std::vector<double>> gradVe_loc = std::vector<std::vector<double>>(nElms);
+ for (auto &vec : gradVe_loc)
+ vec.resize(nNodes, 0.);
+
+ for (size_t i = 1; i < nNodes; ++i)
+ {
+ double val = 1/rhoe[i] * (rhoe[i] * vt[i] * deltaStar[i] - rhoe[i-1] * vt[i-1] * deltaStar[i-1]) / ((loc[i] - loc[i-1]) * (loc[i] - loc[i-1]));
+ gradVe_loc[i-1][i] = -val;
+ gradVe_loc[i-1][i-1] = val;
+ }
+
+ // Compute dloc/dx
+ // Derivative wrt x_j = (x_j - x_j-1) / sqrt((x_j - x_j-1)^2 + (y_j - y_j-1)^2)
+ // Derivative wrt x_j-1 = - (x_j - x_j-1) / sqrt((x_j - x_j-1)^2 + (y_j - y_j-1)^2) = - dloc/dx_j
+ // For loc_i, there are two contributions of x_j except if j = i.
+ // loc_n = sqrt((xn-xn-1)^2) + sqrt((xn-1-xn-2)^2) + ... + sqrt((x1-x0)^2)
+ std::vector<std::vector<double>> gradloc_x = std::vector<std::vector<double>>(nNodes);
+ for (auto &vec : gradloc_x)
+ vec.resize(nNodes, 0.);
+
+ for (size_t j = 1; j < nNodes; ++j)
+ {
+ gradloc_x[j][j] = (x[j] - x[j-1]) / std::sqrt((x[j] - x[j-1]) * (x[j] - x[j-1]) + (y[j] - y[j-1]) * (y[j] - y[j-1]));
+ for (size_t i = 0; i < nNodes - j; ++i)
+ gradloc_x[j+i][j-1] = gradloc_x[j-1][j-1] - gradloc_x[j][j];
+ }
+
+ // Matrix product of gradVe_loc and gradloc_x
+ for (size_t i = 0; i < nElms; ++i)
+ {
+ for (size_t j = 0; j < nNodes; ++j)
+ {
+ for (size_t k = 0; k < nNodes; ++k)
+ gradVe[i][j] += gradVe_loc[i][k] * gradloc_x[k][j];
+ }
+ }
+ return gradVe;
+}
+
+std::vector<std::vector<double>> BoundaryLayer::evalGradwrtY()
+{
+ std::vector<std::vector<double>> gradVe = std::vector<std::vector<double>>(nElms);
+ for (auto &vec : gradVe)
+ vec.resize(nNodes, 0.);
+
+ // Compute dVe/dloc
+ std::vector<std::vector<double>> gradVe_loc = std::vector<std::vector<double>>(nElms);
+ for (auto &vec : gradVe_loc)
+ vec.resize(nNodes, 0.);
+
+ for (size_t i = 1; i < nNodes; ++i)
+ {
+ double val = 1/rhoe[i] * (rhoe[i] * vt[i] * deltaStar[i] - rhoe[i-1] * vt[i-1] * deltaStar[i-1]) / ((loc[i] - loc[i-1]) * (loc[i] - loc[i-1]));
+ gradVe_loc[i-1][i] = -val;
+ gradVe_loc[i-1][i-1] = val;
+ }
+
+ // Compute dloc/dy
+ // Derivative wrt y_j = (y_j - y_j-1) / sqrt((y_j - y_j-1)^2 + (y_j - y_j-1)^2)
+ // Derivative wrt y_j-1 = - (y_j - y_j-1) / sqrt((y_j - y_j-1)^2 + (y_j - y_j-1)^2) = - dloc/dy_j
+ // For loc_i, there are two contributions of y_j except if j = i.
+ // loc_n = sqrt((yn-yn-1)^2) + sqrt((yn-1-yn-2)^2) + ... + sqrt((y1-y0)^2)
+ std::vector<std::vector<double>> gradloc_y = std::vector<std::vector<double>>(nNodes);
+ for (auto &vec : gradloc_y)
+ vec.resize(nNodes, 0.);
+
+ for (size_t j = 1; j < nNodes; ++j)
+ {
+ gradloc_y[j][j] = (y[j] - y[j-1]) / std::sqrt((x[j] - x[j-1]) * (x[j] - x[j-1]) + (y[j] - y[j-1]) * (y[j] - y[j-1]));
+ for (size_t i = 0; i < nNodes - j; ++i)
+ gradloc_y[j+i][j-1] = gradloc_y[j-1][j-1] - gradloc_y[j][j];
+ }
+
+ // Matrix product of gradVe_loc and gradloc_x
+ for (size_t i = 0; i < nElms; ++i)
+ {
+ for (size_t j = 0; j < nNodes; ++j)
+ {
+ for (size_t k = 0; k < nNodes; ++k)
+ gradVe[i][j] += gradVe_loc[i][k] * gradloc_y[k][j];
+ }
+ }
+ return gradVe;
+}
+
/**
* @brief Print solution at a given station.
*
diff --git a/blast/src/DBoundaryLayer.h b/blast/src/DBoundaryLayer.h
index ec04200..030bb1d 100644
--- a/blast/src/DBoundaryLayer.h
+++ b/blast/src/DBoundaryLayer.h
@@ -88,6 +88,7 @@ public:
void setxxExt(std::vector<double> const _xxExt) {if (_xxExt.size() != nNodes) throw std::runtime_error("External mesh vector is not consistent."); xxExt = _xxExt;};
void setDeltaStarExt(std::vector<double> const _deltaStarExt) {if (_deltaStarExt.size() != nNodes) throw std::runtime_error("Displacement thickness vector is not consistent."); deltaStarExt = _deltaStarExt;};
void setVtExt(std::vector<double> const _vtExt) {if (_vtExt.size() != nNodes) throw std::runtime_error("Velocity vector is not consistent."); vtExt = _vtExt;};
+ void setxtrF(double const _xtrF) {xtrF = _xtrF;};
// Boundary conditions.
void setStagnationBC(double Re);
@@ -101,6 +102,11 @@ public:
// Others
void printSolution(size_t iPoint) const;
void computeBlowingVelocity();
+ std::vector<std::vector<double>> evalGradwrtRho();
+ std::vector<std::vector<double>> evalGradwrtVt();
+ std::vector<std::vector<double>> evalGradwrtDeltaStar();
+ std::vector<std::vector<double>> evalGradwrtX();
+ std::vector<std::vector<double>> evalGradwrtY();
};
} // namespace blast
#endif // DBOUNDARYLAYER_H
diff --git a/blast/src/DCoupledAdjoint.cpp b/blast/src/DCoupledAdjoint.cpp
new file mode 100644
index 0000000..8139d32
--- /dev/null
+++ b/blast/src/DCoupledAdjoint.cpp
@@ -0,0 +1,460 @@
+/*
+ * Copyright 2020 University of Liège
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "DCoupledAdjoint.h"
+#include "wAdjoint.h"
+#include "wNewton.h"
+#include "wSolver.h"
+#include "wProblem.h"
+#include "wBlowing.h"
+#include "wFluid.h"
+#include "wMshData.h"
+#include "wNode.h"
+
+#include "wLinearSolver.h"
+#include "wGmres.h"
+
+#include <Eigen/Sparse>
+#include<Eigen/SparseLU>
+#include <Eigen/Dense>
+#include<Eigen/SparseCholesky>
+#include <fstream>
+#include <deque>
+#include <algorithm>
+#include <iomanip>
+
+using namespace blast;
+
+CoupledAdjoint::CoupledAdjoint(std::shared_ptr<dart::Adjoint> _iadjoint, std::shared_ptr<blast::Driver> vSolver) : iadjoint(_iadjoint), vsol(vSolver)
+{
+ isol = iadjoint->sol;
+ // Linear solvers (if GMRES, use the same, but with a thighter tolerance)
+ if (isol->linsol->type() == tbox::LSolType::GMRES_ILUT)
+ {
+ std::shared_ptr<tbox::Gmres> gmres = std::make_shared<tbox::Gmres>(*dynamic_cast<tbox::Gmres *>(isol->linsol.get()));
+ gmres->setTolerance(1e-8);
+ alinsol = gmres;
+ }
+ else
+ alinsol = isol->linsol;
+
+ // Initalize all derivatives
+ this->reset();
+}
+
+void CoupledAdjoint::reset()
+{
+ size_t nDim = isol->pbl->nDim; // Number of dimensions of the problem
+ size_t nNs = isol->pbl->bBCs[0]->nodes.size(); // Number of surface nodes
+ size_t nNd = isol->pbl->msh->nodes.size(); // Number of domain nodes
+ size_t nEs = isol->pbl->bBCs[0]->uE.size(); // Number of elements on the surface
+
+ dRphi_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nNd);
+ dRM_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNd);
+ dRrho_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNd);
+ dRv_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNs, nNd);
+ dRblw_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nNd);
+
+ dRphi_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nNs);
+ dRM_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
+ dRrho_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
+ dRv_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNs, nNs);
+ dRblw_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nNs);
+
+ dRphi_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nNs);
+ dRM_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
+ dRrho_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
+ dRv_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNs, nNs);
+ dRblw_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nNs);
+
+ dRphi_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nDim*nNs);
+ dRM_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim*nNs);
+ dRrho_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim*nNs);
+ dRv_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNs, nDim*nNs);
+ dRblw_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nDim*nNs);
+
+ dRphi_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nEs);
+ dRM_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nEs);
+ dRrho_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nEs);
+ dRv_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNs, nEs);
+ dRblw_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nEs);
+
+ dRphi_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nDim*nNd);
+ dRM_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim*nNd);
+ dRrho_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim*nNd);
+ dRv_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNs, nDim*nNd);
+ dRblw_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nDim*nNd);
+ dRM_M.setIdentity();
+ dRrho_rho.setIdentity();
+ dRv_v.setIdentity();
+ dRblw_blw.setIdentity();
+
+ // Objective functions
+ dCl_phi = Eigen::VectorXd::Zero(nNd);
+ dCd_phi = Eigen::VectorXd::Zero(nNd);
+
+ // Angle of attack
+ dCl_AoA = 0.0;
+ dCd_AoA = 0.0;
+ tdCl_AoA = 0.0;
+ tdCd_AoA = 0.0;
+
+ // Mesh
+ dRx_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNd, nDim*nNd);
+ dCl_x = Eigen::VectorXd::Zero(nDim*nNd);
+ dCd_x = Eigen::VectorXd::Zero(nDim*nNd);
+
+ tdCl_x.resize(nNd, Eigen::Vector3d::Zero());
+ tdCd_x.resize(nNd, Eigen::Vector3d::Zero());
+}
+
+void CoupledAdjoint::buildAdjointMatrix(std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor>*>> &matrices, Eigen::SparseMatrix<double, Eigen::RowMajor> &matrixAdjoint)
+{
+ // Create a list of triplets for the larger matrix
+ std::vector<Eigen::Triplet<double>> triplets;
+
+ // Sanity check
+ for (size_t irow = 0; irow < matrices.size(); ++irow)
+ {
+ if (irow > 0)
+ if (matrices[irow].size() != matrices[irow-1].size()){
+ std::cout << "WRONG ROW SIZE " << std::endl;
+ throw;
+ }
+ for (size_t icol = 0; icol < matrices[irow].size(); ++icol)
+ if (icol > 0)
+ if (matrices[irow][icol]->rows() != matrices[irow][icol-1]->rows()){
+ std::cout << "WRONG COL SIZE " << std::endl;
+ throw;
+ }
+ }
+
+ // Iterate over the rows and columns of the vector
+ int rowOffset = 0;
+ for (const auto& row : matrices) {
+ int colOffset = 0;
+ for (const auto& matrix : row) {
+ // Add the triplets of the matrix to the list of triplets for the larger matrix
+ for (int k = 0; k < matrix->outerSize(); ++k) {
+ for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(*matrix, k); it; ++it) {
+ triplets.push_back(Eigen::Triplet<double>(it.row() + rowOffset, it.col() + colOffset, it.value()));
+ }
+ }
+ colOffset += matrix->cols();
+ }
+ rowOffset += row[0]->rows();
+ }
+ // Create a new matrix from the deque of triplets
+ matrixAdjoint.setFromTriplets(triplets.begin(), triplets.end());
+ matrixAdjoint.prune(0.);
+ matrixAdjoint.makeCompressed();
+}
+
+void CoupledAdjoint::run()
+{
+ gradientswrtInviscidFlow();
+ gradientswrtBlowingVelocity();
+ gradientwrtAoA();
+ gradientwrtMesh();
+
+ transposeMatrices();
+
+ // Adjoint matrix
+
+ // Store pointers to the matrices in a 2D vector
+ std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor>*>> matrices = {
+ {&dRphi_phi, &dRM_phi, &dRv_phi, &dRrho_phi, &dRblw_phi},
+ {&dRphi_M, &dRM_M, &dRv_M, &dRrho_M, &dRblw_M},
+ {&dRphi_v, &dRM_v, &dRv_v, &dRrho_v, &dRblw_v},
+ {&dRphi_rho, &dRM_rho, &dRv_rho, &dRrho_rho, &dRblw_rho},
+ {&dRphi_blw, &dRM_blw, &dRv_blw, &dRrho_blw, &dRblw_blw}
+ };
+
+ int rows = 0; int cols = 0;
+ for (const auto& row : matrices) rows += row[0]->rows(); // All matrices in a row have the same number of rows
+ for (const auto& mat1 : matrices[0]) cols += mat1->cols(); // All matrices in a column have the same number of columns
+ Eigen::SparseMatrix<double, Eigen::RowMajor> A_adjoint(rows, cols);
+ buildAdjointMatrix(matrices, A_adjoint);
+
+ //***** Gradient wrt angle of attack *****//
+ //****************************************//
+ // CL
+ std::vector<double> rhsCl(A_adjoint.cols(), 0.0);
+ for (auto i = 0; i<dCl_phi.size(); ++i)
+ rhsCl[i] = dCl_phi[i];
+ Eigen::VectorXd rhsCl_ = Eigen::Map<const Eigen::VectorXd>(rhsCl.data(), rhsCl.size());
+
+ Eigen::VectorXd lambdaCl(A_adjoint.cols()); // Solution vector for lambdasCl
+ Eigen::Map<Eigen::VectorXd> lambdaCl_(lambdaCl.data(), lambdaCl.size());
+ alinsol->compute(A_adjoint, Eigen::Map<Eigen::VectorXd>(rhsCl_.data(), rhsCl_.size()), lambdaCl_);
+ Eigen::VectorXd lambdaClPhi = lambdaCl.segment(0, isol->pbl->msh->nodes.size());
+ Eigen::VectorXd lambdaClBlw = lambdaCl.segment(lambdaCl.size() - isol->pbl->bBCs[0]->uE.size(), isol->pbl->bBCs[0]->uE.size());
+
+ // Total gradient
+ tdCl_AoA = dCl_AoA - dRphi_AoA.transpose()*lambdaClPhi;
+
+ // CD
+ // std::vector<double> rhsCd(A_adjoint.cols(), 0.0);
+ // for (auto i = 0; i<dCd_phi.size(); ++i)
+ // rhsCd[i] = dCd_phi[i];
+ // Eigen::VectorXd rhsCd_ = Eigen::Map<const Eigen::VectorXd>(rhsCd.data(), rhsCd.size());
+
+ // Eigen::VectorXd lambdaCd(A_adjoint.cols()); // Solution vector for lambdasCd
+ // Eigen::Map<Eigen::VectorXd> lambdaCd_(lambdaCd.data(), lambdaCd.size());
+ // alinsol->compute(A_adjoint, Eigen::Map<Eigen::VectorXd>(rhsCd_.data(), rhsCd_.size()), lambdaCd_);
+ // Eigen::VectorXd lambdaCdPhi = lambdaCd.block(0, 0, isol->pbl->msh->nodes.size(), 1);
+
+ // tdCd_AoA = dCd_AoA - dRphi_AoA.transpose()*lambdaCdPhi; // Total gradient
+
+ //***** Gradient wrt mesh coordinates *****//
+ //*****************************************//
+ // CL
+ // Solution of (from augmented Lagrangian, eqn d/dx )
+ // "dCl_x + dRphi_x * lambdaCl_phi + dRblw_x * lambdaCl_blw + dRx_x * lambdaCl_x"
+ Eigen::VectorXd lambdaCl_x = Eigen::VectorXd::Zero(isol->pbl->nDim * isol->pbl->msh->nodes.size());
+ Eigen::Map<Eigen::VectorXd> lambdaCl_x_(lambdaCl_x.data(), lambdaCl_x.size());
+
+ Eigen::VectorXd rhsCl_x = dCl_x - dRphi_x.transpose() * lambdaClPhi - dRblw_x.transpose() * lambdaClBlw;
+ dRx_x.setIdentity();
+ alinsol->compute(dRx_x.transpose(), Eigen::Map<Eigen::VectorXd>(rhsCl_x.data(), rhsCl_x.size()), lambdaCl_x_);
+
+ for (auto n : isol->pbl->msh->nodes)
+ {
+ for (int m = 0; m < isol->pbl->nDim; m++)
+ {
+ tdCl_x[n->row](m) = lambdaCl_x[isol->pbl->nDim * (n->row) + m];
+ // TODO: CHANGE FOR CD
+ tdCd_x[n->row](m) = lambdaCl_x[isol->pbl->nDim * (n->row) + m];
+ }
+ }
+}
+
+void CoupledAdjoint::gradientswrtInviscidFlow()
+{
+ //**** dRphi_phi ****//
+ 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 independent 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.
+ int 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 (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++;
+ }
+ // 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();
+ 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());
+}
+
+void CoupledAdjoint::gradientswrtMachNumber(){
+ // dRphi_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(isol->pbl->msh->nodes.size(), isol->pbl->bBCs[0]->nodes.size());
+ // dRrho_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(isol->pbl->bBCs[0]->nodes.size(), isol->pbl->bBCs[0]->nodes.size());
+ // dRv_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(isol->pbl->nDim*isol->pbl->bBCs[0]->nodes.size(), isol->pbl->bBCs[0]->nodes.size());
+ // dRblw_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(isol->pbl->bBCs[0]->uE.size(), isol->pbl->bBCs[0]->nodes.size());
+}
+
+void CoupledAdjoint::gradientswrtDensity(){
+ // dRphi_dRho = 0
+ // dRM_dRho = 0
+ // dRrho_dRho =
+ // dRv_dRho = 0
+ // dRblw_dRho =
+
+ // dCl_dRho =
+ // dCd_dRho =
+}
+
+void CoupledAdjoint::gradientswrtVelocity(){
+ // dRphi_dV = 0
+ // dRM_dV = 0
+ // dRrho_dV = 0
+ // dRv_dV =
+ // dRblw_dV =
+
+ // dCl_dV =
+ // dCd_dV =
+}
+
+void CoupledAdjoint::gradientswrtBlowingVelocity(){
+ dRphi_blw = iadjoint->getRu_Blw();
+
+ // All others are zero.
+}
+
+void CoupledAdjoint::gradientwrtAoA(){
+ std::vector<double> dCx_AoA = iadjoint->computeGradientCoefficientsAoa();
+ dCl_AoA = dCx_AoA[0]; dCd_AoA = dCx_AoA[1];
+
+ dRphi_AoA = iadjoint->getRu_A();
+
+ // All others are zero.
+}
+
+void CoupledAdjoint::gradientwrtMesh(){
+ std::vector<std::vector<double>> dCx_x = iadjoint->computeGradientCoefficientsMesh();
+ dCl_x = Eigen::Map<Eigen::VectorXd>(dCx_x[0].data(), dCx_x[0].size());
+ dCd_x = Eigen::Map<Eigen::VectorXd>(dCx_x[1].data(), dCx_x[1].size());
+
+ dRphi_x = iadjoint->getRu_X();
+ dRx_x = iadjoint->getRx_X();
+}
+
+void CoupledAdjoint::setdRblw_M(std::vector<std::vector<double>> _dRblw_dM){
+ //dRblw_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(_dRblw_dM.size(), _dRblw_dM[0].size());
+ // Convert std::vector<double> to Eigen::Triplets
+ std::vector<Eigen::Triplet<double>> triplets;
+ for (size_t i = 0; i < _dRblw_dM.size(); i++){
+ for (size_t j = 0; j < _dRblw_dM[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_dM[i][j]));
+ }
+ }
+ // Set matrix
+ dRblw_M.setFromTriplets(triplets.begin(), triplets.end());
+ dRblw_M.prune(0.);
+}
+void CoupledAdjoint::setdRblw_v(std::vector<std::vector<double>> _dRblw_dv){
+ //dRblw_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(_dRblw_dv.size(), _dRblw_dv[0].size());
+ // Convert std::vector<double> to Eigen::Triplets
+ std::vector<Eigen::Triplet<double>> triplets;
+ for (size_t i = 0; i < _dRblw_dv.size(); i++){
+ for (size_t j = 0; j < _dRblw_dv[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_dv[i][j]));
+ }
+ }
+ // Set matrix
+ dRblw_v.setFromTriplets(triplets.begin(), triplets.end());
+ dRblw_v.prune(0.);
+}
+void CoupledAdjoint::setdRblw_rho(std::vector<std::vector<double>> _dRblw_drho){
+ //dRblw_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(_dRblw_drho.size(), _dRblw_drho[0].size());
+ // Convert std::vector<double> to Eigen::Triplets
+ std::vector<Eigen::Triplet<double>> triplets;
+ for (size_t i = 0; i < _dRblw_drho.size(); i++){
+ for (size_t j = 0; j < _dRblw_drho[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_drho[i][j]));
+ }
+ }
+ // Set matrix
+ dRblw_rho.setFromTriplets(triplets.begin(), triplets.end());
+ dRblw_rho.prune(0.);
+}
+void CoupledAdjoint::setdRblw_x(std::vector<std::vector<double>> _dRblw_dx){
+ //dRblw_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(_dRblw_dx.size(), _dRblw_dx[0].size());
+ // Convert std::vector<double> to Eigen::Triplets
+ std::vector<Eigen::Triplet<double>> triplets;
+ for (size_t i = 0; i < _dRblw_dx.size(); i++){
+ for (size_t j = 0; j < _dRblw_dx[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_dx[i][j]));
+ }
+ }
+ // Set matrix
+ dRblw_x.setFromTriplets(triplets.begin(), triplets.end());
+ dRblw_x.prune(0.);
+}
+
+/**
+ * @brief Transpose all matrices of the adjoint problem included in the system. Not the others.
+ */
+void CoupledAdjoint::transposeMatrices(){
+ dRphi_phi = dRphi_phi.transpose();
+ dRphi_M = dRphi_M.transpose();
+ dRphi_rho = dRphi_rho.transpose();
+ dRphi_v = dRphi_v.transpose();
+ dRphi_blw = dRphi_blw.transpose();
+
+ dRM_phi = dRM_phi.transpose();
+ dRM_M = dRM_M.transpose();
+ dRM_rho = dRM_rho.transpose();
+ dRM_v = dRM_v.transpose();
+ dRM_blw = dRM_blw.transpose();
+
+ dRrho_phi = dRrho_phi.transpose();
+ dRrho_M = dRrho_M.transpose();
+ dRrho_rho = dRrho_rho.transpose();
+ dRrho_v = dRrho_v.transpose();
+ dRrho_blw = dRrho_blw.transpose();
+
+ dRv_phi = dRv_phi.transpose();
+ dRv_M = dRv_M.transpose();
+ dRv_rho = dRv_rho.transpose();
+ dRv_v = dRv_v.transpose();
+ dRv_blw = dRv_blw.transpose();
+
+ dRblw_phi = dRblw_phi.transpose();
+ dRblw_M = dRblw_M.transpose();
+ dRblw_rho = dRblw_rho.transpose();
+ dRblw_v = dRblw_v.transpose();
+ dRblw_blw = dRblw_blw.transpose();
+}
+
+void CoupledAdjoint::writeMatrixToFile(const Eigen::MatrixXd& matrix, const std::string& filename) {
+ std::ofstream file(filename);
+ if (file.is_open()) {
+ for (int i = 0; i < matrix.rows(); ++i) {
+ for (int j = 0; j < matrix.cols(); ++j) {
+ file << std::setprecision(20) << matrix(i, j);
+ if (j != matrix.cols() - 1) {
+ file << " ";
+ }
+ }
+ file << "\n";
+ }
+ file.close();
+ }
+}
diff --git a/blast/src/DCoupledAdjoint.h b/blast/src/DCoupledAdjoint.h
new file mode 100644
index 0000000..3ae603d
--- /dev/null
+++ b/blast/src/DCoupledAdjoint.h
@@ -0,0 +1,128 @@
+/*
+ * Copyright 2020 University of Liège
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef DCOUPLEDADJOINT_H
+#define DCOUPLEDADJOINT_H
+
+#include "blast.h"
+#include "wObject.h"
+#include "dart.h"
+#include <Eigen/Sparse>
+
+#include <iostream>
+#include <memory>
+
+namespace blast
+{
+
+/**
+ * @brief Adjoint of the coupled system.
+ * @author Paul Dechamps
+ */
+class BLAST_API CoupledAdjoint : public fwk::wSharedObject
+{
+
+public:
+ std::shared_ptr<dart::Adjoint> iadjoint; ///< Adjoint solver
+ std::shared_ptr<dart::Newton> isol; ///< Inviscid Newton solver
+ std::shared_ptr<blast::Driver> vsol; ///< Viscous solver
+ std::shared_ptr<tbox::LinearSolver> alinsol; ///< linear solver (adjoint aero)
+
+ double tdCl_AoA; ///< Total derivative of the lift coefficient with respect to the angle of attack
+ double tdCd_AoA; ///< Total erivative of the drag coefficient with respect to the angle of attack
+
+ std::vector<Eigen::Vector3d> tdCl_x; ///< Total derivative of the lift coefficient with respect to the mesh coordinates
+ std::vector<Eigen::Vector3d> tdCd_x; ///< Total derivative of the drag coefficient with respect to the mesh coordinates
+
+private:
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_phi; ///< Inviscid flow Jacobian
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_M; ///< Derivative of the inviscid flow residual with respect to the Mach number
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_v; ///< Derivative of the inviscid flow residual with respect to the velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_rho; ///< Derivative of the inviscid flow residual with respect to the density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_blw; ///< Derivative of the inviscid flow residual with respect to the blowing velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_x; ///< Derivative of the potential residual with respect to the mesh coordinates
+
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_phi; ///< Derivative of the Mach number residual with respect to the inviscid flow
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_M; ///< Derivative of the Mach number residual with respect to the Mach number
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_v; ///< Derivative of the Mach number residual with respect to the velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_rho; ///< Derivative of the Mach number residual with respect to the density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_blw; ///< Derivative of the Mach number residual with respect to the blowing velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_x; ///< Derivative of the Mach number residual with respect to the mesh coordinates
+
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_phi; ///< Derivative of the density residual with respect to the inviscid flow
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_M; ///< Derivative of the density residual with respect to the Mach number
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_v; ///< Derivative of the density residual with respect to the velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_rho; ///< Derivative of the density residual with respect to the density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_blw; ///< Derivative of the density residual with respect to the blowing velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_x; ///< Derivative of the density residual with respect to the mesh coordinates
+
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_phi; ///< Derivative of the velocity residual with respect to the inviscid flow
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_M; ///< Derivative of the velocity residual with respect to the Mach number
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_v; ///< Derivative of the velocity residual with respect to the velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_rho; ///< Derivative of the velocity residual with respect to the density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_blw; ///< Derivative of the velocity residual with respect to the blowing velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_x; ///< Derivative of the velocity residual with respect to the mesh coordinates
+
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_phi; ///< Derivative of the blowing velocity residual with respect to the inviscid flow
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_M; ///< Derivative of the blowing velocity residual with respect to the Mach number
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_v; ///< Derivative of the blowing velocity residual with respect to the velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_rho; ///< Derivative of the blowing velocity residual with respect to the density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_blw; ///< Derivative of the blowing velocity residual with respect to the blowing velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_x; ///< Derivative of the blowing velocity residual with respect to the mesh coordinates
+
+ // Objective functions
+ Eigen::VectorXd dCl_phi; ///< Derivative of the lift coefficient with respect to the inviscid flow
+ Eigen::VectorXd dCd_phi; ///< Derivative of the drag coefficient with respect to the inviscid flow
+
+ // Angle of attack
+ Eigen::VectorXd dRphi_AoA; ///< Derivative of the inviscid flow residual with respect to the angle of attack
+ double dCl_AoA; ///< Partial derivative of the lift coefficient with respect to the angle of attack
+ double dCd_AoA; ///< Partial erivative of the drag coefficient with respect to the angle of attack
+
+ // Mesh
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRx_x; ///< Derivative of the mesh residual with respect to the mesh coordinates
+ Eigen::VectorXd dCl_x; ///< Partial derivative of the lift coefficient with respect to the mesh coordinates
+ Eigen::VectorXd dCd_x; ///< Partial derivative of the drag coefficient with respect to the mesh coordinates
+
+public:
+ CoupledAdjoint(std::shared_ptr<dart::Adjoint> iAdjoint, std::shared_ptr<blast::Driver> vSolver);
+ virtual ~CoupledAdjoint () {std::cout << "~blast::CoupledAdjoint()\n";};
+
+ void run();
+ void reset();
+
+ void setdRblw_M(std::vector<std::vector<double>> dRblw_dM);
+ void setdRblw_v(std::vector<std::vector<double>> dRblw_dv);
+ void setdRblw_rho(std::vector<std::vector<double>> dRblw_drho);
+ 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);
+
+ void gradientswrtInviscidFlow();
+ void gradientswrtMachNumber();
+ void gradientswrtDensity();
+ void gradientswrtVelocity();
+ void gradientswrtBlowingVelocity();
+
+ void gradientwrtAoA();
+ void gradientwrtMesh();
+
+ void transposeMatrices();
+ void writeMatrixToFile(const Eigen::MatrixXd& matrix, const std::string& filename);
+};
+} // namespace blast
+#endif // DDARTADJOINT
diff --git a/blast/src/DDriver.cpp b/blast/src/DDriver.cpp
index ecaf966..2fee439 100644
--- a/blast/src/DDriver.cpp
+++ b/blast/src/DDriver.cpp
@@ -104,7 +104,7 @@ int Driver::run()
// Reset transition
if (reg->name == "wake")
reg->xtr = 0.;
- else if (reg-> name == "upper" || reg->name == "lower")
+ else if (reg->name == "upper" || reg->name == "lower")
reg->xtr = 1.;
else
{
@@ -222,15 +222,15 @@ int Driver::run()
void Driver::checkWaves(size_t iPoint, BoundaryLayer *bl)
{
- double logRet_crit = 2.492 * std::pow(1 / (bl->Hk[iPoint] - 1), 0.43) + 0.7 * (std::tanh(14 * (1 / (bl->Hk[iPoint] - 1)) - 9.24) + 1.0);
+ double logRet_crit = 2.492 * std::pow(1. / (bl->Hk[iPoint] - 1), 0.43) + 0.7 * (std::tanh(14 * (1. / (bl->Hk[iPoint] - 1)) - 9.24) + 1.0);
if (bl->Ret[iPoint] > 0) // Avoid log of negative number.
{
if (std::log10(bl->Ret[iPoint]) <= logRet_crit - 0.08)
- bl->u[iPoint * bl->getnVar() + 2] = 0; // Set N to 0 at that point if waves do not grow.
+ bl->u[iPoint * bl->getnVar() + 2] = 0.; // Set N to 0 at that point if waves do not grow.
}
else
- bl->u[iPoint * bl->getnVar() + 2] = 0; // Set N to 0 at that point if waves do not grow.
+ bl->u[iPoint * bl->getnVar() + 2] = 0.; // Set N to 0 at that point if waves do not grow.
}
/**
@@ -340,7 +340,6 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced)
void Driver::computeSectionalDrag(std::vector<BoundaryLayer *> const &sec, size_t i)
{
size_t nVar = sec[0]->getnVar();
- size_t lastWkPt = (sec[2]->nNodes - 1) * nVar;
// Normalize friction and dissipation coefficients.
std::vector<double> frictioncoeff_upper(sec[0]->nNodes, 0.0);
@@ -361,8 +360,20 @@ void Driver::computeSectionalDrag(std::vector<BoundaryLayer *> const &sec, size_
// Friction drag coefficient.
Cdf_sec[i] = Cdf_upper + Cdf_lower; // No friction in the wake
- // Total drag coefficient.
- Cdt_sec[i] = 2 * sec[2]->u[lastWkPt + 0] * pow(sec[2]->u[lastWkPt + 3], (sec[2]->u[lastWkPt + 1] + 5) / 2);
+
+ // Total drag coefficient. Compute from upper and lower TE if there is no wake.
+ if (sec.size() > 2 && sec[2]->getName() == "wake")
+ {
+ size_t lastWkPt = (sec[2]->nNodes - 1) * nVar;
+ Cdt_sec[i] = 2 * sec[2]->u[lastWkPt + 0] * pow(sec[2]->u[lastWkPt + 3], (sec[2]->u[lastWkPt + 1] + 5) / 2);
+ }
+ else
+ {
+ size_t lastUpPt = (sec[0]->nNodes - 1) * nVar;
+ size_t lastLwPt = (sec[1]->nNodes - 1) * nVar;
+ Cdt_sec[i] = sec[0]->u[lastUpPt + 0] * pow(sec[0]->u[lastUpPt + 3], (sec[0]->u[lastUpPt + 1] + 5) / 2) +
+ sec[1]->u[lastLwPt + 0] * pow(sec[1]->u[lastLwPt + 3], (sec[1]->u[lastLwPt + 1] + 5) / 2);
+ }
// Pressure drag coefficient.
Cdp_sec[i] = Cdt_sec[i] - Cdf_sec[i];
}
diff --git a/blast/src/DDriver.h b/blast/src/DDriver.h
index ed598aa..c04b1e8 100644
--- a/blast/src/DDriver.h
+++ b/blast/src/DDriver.h
@@ -47,6 +47,9 @@ public:
double getRe() const { return Re; };
std::vector<std::vector<double>> getSolution(size_t iSec);
+ // Setters.
+ void setVerbose(unsigned int _verbose) { verbose = _verbose; };
+
void addSection(size_t iSec, BoundaryLayer * &bl){sections[iSec].push_back(bl);};
private:
diff --git a/blast/src/blast.h b/blast/src/blast.h
index 382e040..0fcf795 100644
--- a/blast/src/blast.h
+++ b/blast/src/blast.h
@@ -41,6 +41,9 @@ namespace blast
class Driver;
class Solver;
+// Adjoint
+class CoupledAdjoint;
+
// Data structures
class BoundaryLayer;
diff --git a/blast/tests/dart/adjointVII_2D.py b/blast/tests/dart/adjointVII_2D.py
new file mode 100644
index 0000000..7341e98
--- /dev/null
+++ b/blast/tests/dart/adjointVII_2D.py
@@ -0,0 +1,174 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2022 University of Liège
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+
+# @author Paul Dechamps
+# @date 2024
+# Test the blast adjoint implementation.
+
+# Imports.
+
+import blast.utils as viscUtils
+from blast.interfaces.DAdjointInterface import AdjointInterface
+import blast
+import numpy as np
+
+from fwk.wutils import parseargs
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+import numpy as np
+
+import math
+
+def cfgInviscid(nthrds, verb):
+ import os.path
+ # Parameters
+ return {
+ # Options
+ 'scenario' : 'aerodynamic',
+ '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' : 50, 'yLgt' : 50, 'msF' : 20, 'msTe' : 0.01, 'msLe' : 0.01}, # parameters for input file model
+ 'Dim' : 2, # problem dimension
+ 'Format' : 'gmsh', # save format (vtk or gmsh)
+ # Markers
+ 'Fluid' : 'field', # name of physical group containing the fluid
+ 'Farfield' : ['upstream', 'farfield', 'downstream'], # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
+ 'Wing' : 'airfoil', # LIST of names of physical groups containing the lifting surface boundary
+ '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' : False,
+ 'Upstream' : 'upstream',
+ # Freestream
+ 'M_inf' : 0.2, # freestream Mach number
+ 'AoA' : 2., # freestream angle of attack
+ # Geometry
+ 'S_ref' : 1., # reference surface length
+ 'c_ref' : 1., # reference chord length
+ 'x_ref' : .25, # reference point for moment computation (x)
+ '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)
+ 'G_fill' : 2, # 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
+ }
+
+def cfgBlast(verb):
+ return {
+ 'Re' : 1e7, # Freestream Reynolds number
+ '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': 200, # Maximum number of iterations
+ 'couplTol' : 1e-8, # Tolerance of the VII methodology
+ 'iterPrint': 20, # int, number of iterations between outputs
+ 'resetInv' : True, # bool, flag to reset the inviscid calculation at every iteration.
+ 'sections' : [0], # List of sections for boundary layer calculation
+ 'xtrF' : [0.2, 0.2], # Forced transition location
+ 'interpolator' : 'Matching', # Interpolator for the coupling
+ }
+
+def main():
+ # Timer.
+ tms = fwk.Timers()
+ tms['total'].start()
+
+ args = parseargs()
+ icfg = cfgInviscid(args.k, args.verb)
+ vcfg = cfgBlast(args.verb)
+
+ alpha = icfg['AoA']*np.pi/180
+ dalpha = 1e-6 * np.pi/180
+ daoa = [alpha - dalpha, alpha + dalpha]
+ dcl = []
+
+ coupler, isol, vsol = viscUtils.initBlast(icfg, vcfg)
+
+ print(ccolors.ANSI_BLUE + 'PyInviscid...' + ccolors.ANSI_RESET)
+ isol.run()
+ # Inviscid finite difference
+ for aa in daoa:
+ isol.setAoA(aa)
+ isol.reset()
+ isol.run()
+ dcl.append(isol.getCl())
+ dcl_daoa_dart_fd = (dcl[-1] - dcl[-2])/(2*dalpha)
+ isol.adjointSolver.run()
+ dcl_daoa_dart_ad = isol.adjointSolver.dClAoa
+
+ # Adjoint objects
+ cAdj = blast.CoupledAdjoint(isol.adjointSolver, vsol)
+ pyAdjoint = AdjointInterface(cAdj, isol, vsol)
+
+ # Run coupled case
+ print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+ tms['forward'].start()
+ aeroCoeffs = coupler.run()
+ tms['forward'].stop()
+
+ print(ccolors.ANSI_BLUE + 'PyDartAdjoint...' + ccolors.ANSI_RESET)
+ tms['adj_dart'].start()
+ isol.adjointSolver.run()
+ tms['adj_dart'].stop()
+
+ print(ccolors.ANSI_BLUE + 'PyBlasterAdjoint...' + ccolors.ANSI_RESET)
+ tms['adj_py'].start()
+ pyAdjoint.build()
+ tms['adj_py'].stop()
+ tms['adj_blast'].start()
+ cAdj.run()
+ tms['adj_blast'].stop()
+
+ # Coupled finite difference
+ tms['fd_blast'].start()
+ dcl = []
+ for aa in daoa:
+ isol.setAoA(aa)
+ coupler.reset()
+ aeroCoeffs = coupler.run(write=False)
+ dcl.append(isol.getCl())
+ dcl_daoa_blast_fd = (dcl[-1] - dcl[-2])/(2*dalpha)
+ tms['fd_blast'].stop()
+
+ print(ccolors.ANSI_BLUE + 'PyResults' + ccolors.ANSI_RESET)
+ print('--------- Dart ---------')
+ print('AD: ' + str(dcl_daoa_dart_ad))
+ print('FD: ' + str(dcl_daoa_dart_fd))
+ print('Error: ' + str(abs(dcl_daoa_dart_ad - dcl_daoa_dart_fd)/dcl_daoa_dart_ad))
+ print('--------- Blaster ---------')
+ print('AD: ' + str(cAdj.tdCl_AoA))
+ print('FD: ' + str(dcl_daoa_blast_fd))
+ print('Error: ' + str(abs(cAdj.tdCl_AoA - dcl_daoa_blast_fd)/cAdj.tdCl_AoA))
+ tms['total'].stop()
+ print(tms)
+
+if __name__ == '__main__':
+ main()
\ No newline at end of file
diff --git a/blast/tests/dart/naca0012_2D.py b/blast/tests/dart/naca0012_2D.py
index 361b760..04ab790 100644
--- a/blast/tests/dart/naca0012_2D.py
+++ b/blast/tests/dart/naca0012_2D.py
@@ -88,17 +88,17 @@ def cfgInviscid(nthrds, verb):
def cfgBlast(verb):
return {
- 'Re' : 1e7, # Freestream Reynolds number
- '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': 50, # 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.
- 'sections' : [0],
- 'xtrF' : [None, 0.4],# Forced transition location
- 'interpolator' : 'Matching'
+ 'Re' : 1e7, # Freestream Reynolds number
+ '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': 50, # 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.
+ 'sections' : [0], # List of sections for boundary layer calculation
+ 'xtrF' : [None, 0.4], # Forced transition location
+ 'interpolator' : 'Matching' # Interpolator for the coupling
}
def main():
--
GitLab