diff --git a/blast/api/mdaAPI.py b/blast/api/mdaAPI.py
new file mode 100644
index 0000000000000000000000000000000000000000..43fc618cf91f0024418d3b763a26f61f3c9d11db
--- /dev/null
+++ b/blast/api/mdaAPI.py
@@ -0,0 +1,294 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2024 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.
+
+
+## BLASTER api for openMDAO
+# Paul Dechamps
+
+import numpy as np
+import openmdao.api as om
+import blast
+import blast.utils as viscUtils
+
+# Solver openMDAO object
+class BlasterSolver(om.ExplicitComponent):
+ """OpenMDAO component for the Blaster solver
+
+ Parameters
+ ----------
+ - cfg {dict} -- Configuration dictionary for the VII coupling
+ - icfg {dict} -- Configuration dictionary for the inviscid solver interface
+ - iSolverName {str} -- Inviscid solver name (default: 'DART')
+
+ Inputs
+ ------
+ - aoa {float} -- Angle of attack
+ - x_aero {np.array} -- Mesh coordinates
+
+ Outputs
+ -------
+ - cl {float} -- Lift coefficient
+ - cd {float} -- Drag coefficient
+ """
+
+ def __init__(self, cfg, icfg, iSolverName='DART'):
+ super().__init__()
+
+ self.nCpt = 0
+
+ if iSolverName == 'DART':
+ from blast.interfaces.dart.DartInterface import DartInterface as interface
+ else:
+ raise RuntimeError('Solver '+iSolverName+' currently not implemented')
+
+ for i in range(len(cfg['xtrF'])):
+ if cfg['xtrF'][i] is None:
+ cfg['xtrF'][i] = -1
+ self.vsol = blast.Driver(cfg['Re'], cfg['Minf'], cfg['CFL0'], cfg['nSections'], xtrF_top=cfg['xtrF'][0], xtrF_bot=cfg['xtrF'][1], _span=cfg['span'], _verbose=cfg['Verb'])
+
+ self.isol = interface(icfg, self.vsol, cfg)
+
+ # Coupler
+ import blast.coupler as blastCoupler
+ self.coupler = blastCoupler.Coupler(self.isol, self.vsol, _maxCouplIter=cfg['couplIter'], _couplTol=cfg['couplTol'], _iterPrint=cfg['iterPrint'], _resetInv=cfg['resetInv'])
+
+ # Adjoint
+ self.adjointSolver = blast.CoupledAdjoint(self.isol.adjointSolver, self.vsol)
+
+ # Single coupler run to initialize viscous sections (which depend on size of upper/lower sides)
+ self.coupler.run(write=False)
+ self.coupler.reset()
+ print('object created')
+
+ def setup(self):
+ self.add_input('aoa', val=0.0)
+ self.add_input('x_aero', val=np.zeros(self.isol.iobj['bnd'].nodes.size() * self.isol.getnDim()))
+ self.add_output('cl', val=0.0)
+ self.add_output('cd', val=0.0)
+
+ def setup_partials(self):
+ self.declare_partials('cl', 'aoa', method='exact')
+ self.declare_partials('cd', 'aoa', method='exact')
+ self.declare_partials('cl', 'x_aero', method='exact')
+ self.declare_partials('cd', 'x_aero', method='exact')
+
+
+ def compute_partials(self, inputs, partials, discrete_inputs=None):
+ """Compute partial derivatives using the adjoint solver
+ """
+ # Run adjoint solver
+ self.isol.iobj['adj'].run()
+ self.adjointSolver.reset()
+ self.adjointSolver.run()
+
+ # AoA derivatives
+ partials['cl', 'aoa'] = self.adjointSolver.tdCl_AoA
+ partials['cd', 'aoa'] = self.adjointSolver.tdCd_AoA
+
+ # Mesh derivative
+ nDim = self.isol.getnDim()
+ # Get mesh matrix
+ for inod, n in enumerate(self.isol.iobj['bnd'].nodes):
+ for idim in range(nDim):
+ partials['cl', 'x_aero'][0][inod*nDim + idim] = self.adjointSolver.tdCl_x[n.row][idim]
+ partials['cd', 'x_aero'][0][inod*nDim + idim] = self.adjointSolver.tdCd_x[n.row][idim]
+
+ def compute(self, inputs, outputs):
+ """Update the mesh, aoa and run the solver
+ """
+ nDim = self.isol.getnDim()
+ # Modify surface and deform mesh
+ self.isol.iobj['mrf'].savePos()
+ for inod, n in enumerate(self.isol.iobj['bnd'].nodes):
+ for idim in range(nDim):
+ n.pos[idim] = inputs['x_aero'][inod*nDim + idim]
+ self.isol.iobj['mrf'].deform()
+
+ # Update viscous mesh
+ self.isol.updateMesh()
+
+ # Update angle of attack
+ aoa = inputs['aoa'][0]
+ self.isol.setAoA(aoa*np.pi/180) # aoa is given in deg
+
+ # Reset coupler
+ self.coupler.reset()
+
+ # Run forward problem
+ self.coupler.run(write=False)
+ outputs['cl'] = self.isol.getCl()
+ outputs['cd'] = self.isol.getCd() + self.vsol.Cdf
+
+ # Write results
+ self.isol.writeCp(sfx='_'+str(self.nCpt))
+ viscUtils.getSolution(self.isol.sec, write=True, toW=['x', 'y', 'cf'], sfx='_'+str(self.nCpt))
+ self.nCpt += 1
+
+ def getBoundaryMesh(self):
+ """Initial surface mesh coordinates
+ """
+ return BlasterMesh(nDim=self.isol.getnDim(), bnd=self.isol.iobj['bnd'])
+
+ def getBoundaryGeometry(self):
+ """Airfoil geometry
+ """
+ return BlasterGeometry(nDim=self.isol.getnDim(), nNodes=self.isol.iobj['bnd'].nodes.size(), nodeRows=self.isol.vBnd[0][0].connectListNodes)
+
+# Surface mesh
+class BlasterMesh(om.IndepVarComp):
+ """Initial 2D surface mesh coordinates
+
+ Parameters
+ ----------
+ - nDim {int} -- Number of dimensions
+ - bnd {object} -- Boundary object containing the mesh
+
+ Outputs
+ -------
+ - x_aero0 {np.array} -- Initial mesh coordinates
+ """
+
+ def initialize(self):
+ self.options.declare('nDim')
+ self.options.declare('bnd', recordable=False)
+
+ def setup(self):
+ bnd = self.options['bnd']
+ nDim = self.options['nDim']
+ x_aero0 = np.zeros(bnd.nodes.size() * nDim)
+ for inod, n in enumerate(bnd.nodes):
+ for idim in range(nDim):
+ x_aero0[inod*nDim + idim] = n.pos[idim]
+ self.add_output('x_aero0', val=x_aero0)
+
+class BlasterGeometry(om.ExplicitComponent):
+ """Airfoil parametrization as a NACA 4-digit airfoil.
+ Uses the NACA 4-digit airfoil parametrization to generate an airfoil from a given x
+ distribution.
+ Parameters become tau (thickness), eps (max camber) and p (max camber position).
+
+ Parameters
+ ----------
+ - nDim {int} -- Number of dimensions
+ - nNodes {int} -- Number of nodes
+ - nodeRows {np.array} -- Nodal indexes
+
+ Inputs
+ ------
+ - x_aero0 {np.array} -- Initial mesh
+ - tau {float} -- Airfoil thickness
+ - eps {float} -- Airfoil max camber
+ - p {float} -- Airfoil max camber position
+
+ Outputs
+ -------
+ - x_aero_out {np.array} -- Airfoil coordinates
+ """
+
+ def initialize(self):
+ self.options.declare('nDim')
+ self.options.declare('nNodes')
+ self.options.declare('nodeRows')
+
+ def setup(self):
+ nNodes = self.options['nNodes']
+ nDim = self.options['nDim']
+ self.add_input('x_aero0', np.zeros(nNodes*nDim)) # Initial mesh
+ self.add_input('tau', val = 0.12, desc='Airfoil thickness')
+ self.add_input('eps', val = 0., desc='Airfoil max camber')
+ self.add_input('p', val = 0., desc='Airfoil max camber position')
+ self.add_output('x_aero_out', val = np.zeros(nNodes*nDim))
+
+ def setup_partials(self):
+ self.declare_partials('x_aero_out', 'tau', method='fd')
+ self.declare_partials('x_aero_out', 'eps', method='fd')
+ self.declare_partials('x_aero_out', 'p', method='fd')
+
+ def compute(self, inputs, outputs):
+ """Compute airfoil coordinates
+ Points are sorted in seilig format to identify upper and lower surfaces.
+ Then the NACA 4-digit airfoil is generated and the points are sorted back
+ to the original mesh format."""
+
+ x_aero0 = inputs['x_aero0']
+ tau = inputs['tau'][0]
+ eps = inputs['eps'][0]
+ p = inputs['p'][0]
+ nDim = self.options['nDim']
+
+ # Retreive selig coordinates
+ x_seilig = np.zeros(len(self.options['nodeRows']))
+ for i, val in enumerate(self.options['nodeRows']):
+ x_seilig[i] = x_aero0[val*nDim]
+ # Generate airfoil
+ # Format tau eps and p with 2 digits after coma
+ print(f"tau: {tau:.10f}, eps: {eps:.10f}, p: {p:.10f}")
+ xb, yb = self.__generateAirfoil(x_seilig, tau, eps, p)
+ # Go back to normal mesh
+ x_aero_out = np.zeros(len(x_aero0))
+ for i, val in enumerate(self.options['nodeRows']):
+ x_aero_out[val*nDim] = xb[i]
+ x_aero_out[val*nDim+1] = yb[i]
+ outputs['x_aero_out'] = x_aero_out
+
+ def __generateAirfoil(self, x, tau, eps, p):
+ """Generate NACA 4-digit airfoil from x coordinates points
+
+ Arguments
+ ---------
+ - x {np.array} -- x coordinates sorted in seilig format
+ - tau {float} -- Airfoil thickness
+ - eps {float} -- Airfoil max camber
+ - p {float} -- Airfoil max camber position
+
+ Returns
+ -------
+ - tuple -- x and y coordinates
+ """
+
+ T = 10 * tau * (
+ 0.2969 * np.sqrt(x)
+ - 0.1260 * (x)
+ - 0.3537 * (x)**2
+ + 0.2843 * (x)**3
+ - 0.1015 * (x)**4) # Thickness
+
+ YBar = np.zeros(len(T))
+ dYBar_dx = np.zeros(len(T))
+ Y = np.zeros(len(T))
+ xB = np.zeros(len(T))
+ yB = np.zeros(len(T))
+
+ for i, xi in enumerate(x):
+ # Camber
+ if xi < p:
+ YBar[i] = ((eps * xi)/(p**2)) * (2 * p - xi)
+ dYBar_dx[i] = (2 * eps)/(p**2) * (p - xi)
+ elif xi >= p:
+ YBar[i] = (eps * (1 - xi)/((1-p)**2)) * (1 + xi - 2 * p)
+ dYBar_dx[i] = (2 * eps)/((1 - p)**2) * (p - xi)
+
+ # Final x & y coordinates
+ if i >= len(x)/2 - 1:
+ # Lower surface
+ yB[i] = YBar[i] - 0.5 * T[i] * np.cos(np.arctan(dYBar_dx[i]))
+ xB[i] = x[i] + 0.5 * T[i] * np.sin(np.arctan(dYBar_dx[i]))
+ else:
+ # Upper surface
+ yB[i] = YBar[i] + 0.5 * T[i] * np.cos(np.arctan(dYBar_dx[i]))
+ xB[i] = x[i] - 0.5 * T[i] * np.sin(np.arctan(dYBar_dx[i]))
+ return (xB, yB)
diff --git a/blast/coupler.py b/blast/coupler.py
index b7213a4c4c9ed54404316d8377338bdf45d600cd..2878d67f8d1bb29607c63671310a522894251a5a 100644
--- a/blast/coupler.py
+++ b/blast/coupler.py
@@ -85,6 +85,7 @@ class Coupler:
# Write inviscid Cp distribution.
if couplIter == 0:
self.isol.initializeViscousSolver()
+ self.isol.updateStagnation()
if write:
self.isol.writeCp(sfx='_inviscid'+self.filesfx)
@@ -142,20 +143,6 @@ class Coupler:
# 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:
@@ -166,3 +153,28 @@ class Coupler:
for i in range(len(side.blowingVel)):
side.blowingVel[i] = 0.
self.isol.setBlowingVelocity()
+
+ # Viscous solver
+ if self.isol.vinit:
+ for s in self.isol.sec:
+ for reg in s:
+ reg.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())
+ self.isol.vtExt[iSec][i] = np.zeros(reg.getnNodes())
+ reg.setxxExt(self.isol.xxExt[iSec][i])
+ reg.setVtExt(self.isol.vtExt[iSec][i])
+ reg.setDeltaStarExt(self.isol.deltaStarExt[iSec][i])
+
diff --git a/blast/interfaces/DAdjointInterface.py b/blast/interfaces/DAdjointInterface.py
deleted file mode 100644
index cf3f60ff13934985b9e8c65a44c291658d1f055e..0000000000000000000000000000000000000000
--- a/blast/interfaces/DAdjointInterface.py
+++ /dev/null
@@ -1,102 +0,0 @@
-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 4a63dc6302cf5a6839dbdbe6122b9a1c018b0f35..a5049f6ce9f6cef58c64e9b45f064ef73bc90264 100644
--- a/blast/interfaces/DDataStructure.py
+++ b/blast/interfaces/DDataStructure.py
@@ -69,6 +69,23 @@ class Group:
else:
raise RuntimeError('Invalid region', reg)
+ def getNodesRow(self, reg:str):
+ """Return the row of the nodes
+
+ args:
+ ----
+ reg: str
+ Region of the boundary layer (upper, lower, wake)
+ """
+ if reg == "wake":
+ return self.connectListNodes
+ elif reg == "upper":
+ return np.flip(self.connectListNodes[:self.stagPoint+1])
+ elif reg == "lower":
+ return self.connectListNodes[self.stagPoint:]
+ else:
+ raise RuntimeError('Invalid region', reg)
+
def getVelocity(self, reg:str, dim:int):
"""Return the velocity in the direction dim
@@ -140,6 +157,23 @@ class Group:
return self.Rho[self.stagPoint:]
else:
raise RuntimeError('Invalid region', reg)
+
+ def getConnectList(self, reg:str):
+ """Return the connectivity list of the elements
+
+ args:
+ ----
+ reg: str
+ Region of the boundary layer (upper, lower, wake)
+ """
+ if reg == "wake":
+ return self.connectListElems
+ elif reg == "upper":
+ return np.flip(self.connectListElems[:self.stagPoint])
+ elif reg == "lower":
+ return self.connectListElems[self.stagPoint:]
+ else:
+ raise RuntimeError('Invalid region', reg)
def setConnectList(self, _connectListNodes, _connectListElems):
"""Set connectivity lists for viscous structures.
diff --git a/blast/interfaces/DSolversInterface.py b/blast/interfaces/DSolversInterface.py
index dc85a470676297b4f60acb7850aed6040a3db6f9..b7542bc56ab9b64d9600921f416f5607e88d3053 100644
--- a/blast/interfaces/DSolversInterface.py
+++ b/blast/interfaces/DSolversInterface.py
@@ -1,7 +1,7 @@
import numpy as np
import blast
+import tbox
from blast.interfaces.DDataStructure import Group
-import time
class SolversInterface:
def __init__(self, _vSolver, _cfg):
@@ -16,23 +16,22 @@ class SolversInterface:
self.vBnd = [[Group('vAirfoil'+str(k) if iReg == 0 else 'vWake'+str(k), self.getnDim()) for iReg in range(2)] for k in range(self.cfg['nSections'])]
# Initialize interpolator and meshes
- initialTime = time.time()
- self.setMesh()
+ self.getWing()
if self.cfg['interpolator'] == 'Matching':
from blast.interfaces.interpolators.DMatchingInterpolator import MatchingInterpolator as interp
# Initialize viscous structures for matching computations
for iSec in range(len(self.vBnd)):
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
+ self.vBnd[iSec][iReg].nodesCoord = self.iBnd[iReg].nodesCoord.copy()
+ self.vBnd[iSec][iReg].elemsCoord = self.iBnd[iReg].elemsCoord.copy()
+ self.vBnd[iSec][iReg].connectListNodes = self.iBnd[iReg].connectListNodes.copy()
+ self.vBnd[iSec][iReg].connectListElems = self.iBnd[iReg].connectListElems.copy()
elif self.cfg['interpolator'] == 'Rbf':
from blast.interfaces.interpolators.DRbfInterpolator import RbfInterpolator as interp
- print('RBF interpolator initialized.')
- print('Loading viscous mesh... ', end='')
+ print('Initializing RBF interpolator.')
self.getVWing()
- print('done in {:.2f}s.'.format(time.time()-initialTime))
else:
raise RuntimeError('Incorrect interpolator specified in config.')
@@ -47,7 +46,7 @@ class SolversInterface:
for iSec in range(self.cfg['nSections'])]
self.vinit = False
- def setViscousSolver(self, adjoint=False):
+ def setViscousSolver(self):
"""Interpolate solution to viscous mesh and sets the inviscid boundary in the viscous solver.
"""
for iSec in range(self.cfg['nSections']):
@@ -63,11 +62,11 @@ class SolversInterface:
reg.setRhoe(self.vBnd[iSec][iReg].getRho(reg.name))
reg.setVt(self.vBnd[iSec][iReg].getVt(reg.name))
- if adjoint == False:
- reg.setDeltaStarExt(self.deltaStarExt[iSec][i])
- reg.setxxExt(self.xxExt[iSec][i])
+ reg.setDeltaStarExt(self.deltaStarExt[iSec][i])
+ reg.setxxExt(self.xxExt[iSec][i])
+ reg.setVtExt(self.vtExt[iSec][i])
- def getBlowingBoundary(self, couplIter, adjoint=False):
+ def getBlowingBoundary(self, couplIter):
"""Get blowing boundary condition from the viscous solver.
args:
----
@@ -82,11 +81,11 @@ 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)))
-
- 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)
+
+ 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.vtExt[iSec][i] = np.asarray(self.sec[iSec][i].vt)
def interpolate(self, dir):
if dir == 'i2v':
@@ -96,6 +95,102 @@ class SolversInterface:
else:
raise RuntimeError('Invalid direction', dir)
+ def updateMesh(self):
+ """Collect the mesh on blowing boundaries and update the mesh in the viscous solver.
+ This function does not update connectivity lists
+ """
+ if self.getnDim() != 2:
+ raise RuntimeError('SolverInterface::updateMesh - Only 2D meshes are supported.')
+ if self.cfg['interpolator'] != 'Matching':
+ raise RuntimeError('SolverInterface::updateMesh - Only Matching interpolator is supported.')
+
+ meshList = self.getMesh()
+
+ # Update the mesh in data structures
+ for ireg, reg in enumerate(self.iBnd):
+ for i, row in enumerate(reg.connectListNodes):
+ for idim in range(self.getnDim()):
+ reg.nodesCoord[i, idim] = meshList[ireg][row, idim]
+
+ # Recompute elements data
+ for ireg, reg in enumerate(self.iBnd):
+ for i in range(len(reg.nodesCoord[:,0])-1):
+ for idim in range(self.getnDim()):
+ reg.elemsCoord[i, idim] = (reg.nodesCoord[i, idim] + reg.nodesCoord[i+1, idim])/2
+
+ # Update the mesh in the data structures
+ for iSec in range(len(self.vBnd)):
+ for iReg, reg in enumerate(self.vBnd[iSec]):
+ self.vBnd[iSec][iReg].nodesCoord = self.iBnd[iReg].nodesCoord.copy()
+ self.vBnd[iSec][iReg].elemsCoord = self.iBnd[iReg].elemsCoord.copy()
+
+ # Set mesh in the viscous solver
+ for iSec in range(self.cfg['nSections']):
+ # Compute section chord
+ xMin = np.min(self.vBnd[iSec][0].nodesCoord[:,0])
+ xMax = np.max(self.vBnd[iSec][0].nodesCoord[:,0])
+ chord = xMax - xMin
+ for reg in self.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)
+
+ xIdx = 0
+ if self.getnDim() == 2:
+ yIdx = 1
+ zIdx = 2
+ elif self.getnDim() == 3:
+ yIdx = 2
+ zIdx = 1
+ else:
+ raise RuntimeError('Incorrect number of dimensions', self.getnDim())
+
+ nodeRows = tbox.std_vector_int()
+ if self.cfg['interpolator'] == 'Matching':
+ for val in self.vBnd[iSec][iReg].getNodesRow(reg.name):
+ nodeRows.push_back(int(val))
+ else:
+ for i in range(len(self.vBnd[iSec][iReg].getNodesCoord(reg.name, xIdx))):
+ nodeRows.push_back(int(0))
+
+ connectElems = tbox.std_vector_int()
+ if self.cfg['interpolator'] == 'Matching':
+ for val in self.vBnd[iSec][iReg].getConnectList(reg.name):
+ connectElems.push_back(int(val))
+ else:
+ for i in range(len(self.vBnd[iSec][iReg].getNodesCoord(reg.name, xIdx))-1):
+ connectElems.push_back(int(0))
+
+ # x, y, z
+ xv = tbox.std_vector_double()
+ for val in self.vBnd[iSec][iReg].getNodesCoord(reg.name, xIdx):
+ xv.push_back(val)
+ yv = tbox.std_vector_double()
+ for val in self.vBnd[iSec][iReg].getNodesCoord(reg.name, yIdx):
+ yv.push_back(val)
+ zv = tbox.std_vector_double()
+ for val in self.vBnd[iSec][iReg].getNodesCoord(reg.name, zIdx):
+ zv.push_back(val)
+
+ # Set mesh
+ reg.setMesh(xv,
+ yv,
+ zv,
+ chord,
+ xMin,
+ nodeRows,
+ connectElems)
+
+ def updateStagnation(self):
+ """Update the stagnation points in the viscous solver.
+ """
+ for iSec in range(self.cfg['nSections']):
+ for reg in self.vBnd[iSec]:
+ reg.computeStagPoint()
+
def initializeViscousSolver(self):
"""Initialize the viscous sections in the viscous solver.
"""
@@ -117,7 +212,6 @@ class SolversInterface:
name = "upper"
xtrF = self.cfg['xtrF'][j]
-
self.sec[i].append(blast.BoundaryLayer(xtrF, name))
self.vSolver.addSection(i, self.sec[i][j])
@@ -153,12 +247,41 @@ class SolversInterface:
else:
raise RuntimeError('Incorrect number of dimensions', self.getnDim())
- reg.setMesh(self.vBnd[iSec][iReg].getNodesCoord(reg.name, xIdx),
- self.vBnd[iSec][iReg].getNodesCoord(reg.name, yIdx),
- self.vBnd[iSec][iReg].getNodesCoord(reg.name, zIdx),
- chord,
- xMin)
+ nodeRows = tbox.std_vector_int()
+ if self.cfg['interpolator'] == 'Matching':
+ for val in self.vBnd[iSec][iReg].getNodesRow(reg.name):
+ nodeRows.push_back(int(val))
+ else:
+ for i in range(len(self.vBnd[iSec][iReg].getNodesCoord(reg.name, xIdx))):
+ nodeRows.push_back(int(0))
+ connectElems = tbox.std_vector_int()
+ if self.cfg['interpolator'] == 'Matching':
+ for val in self.vBnd[iSec][iReg].getConnectList(reg.name):
+ connectElems.push_back(int(val))
+ else:
+ for i in range(len(self.vBnd[iSec][iReg].getNodesCoord(reg.name, xIdx))-1):
+ connectElems.push_back(int(0))
+
+ # x, y, z
+ xv = tbox.std_vector_double()
+ for val in self.vBnd[iSec][iReg].getNodesCoord(reg.name, xIdx):
+ xv.push_back(val)
+ yv = tbox.std_vector_double()
+ for val in self.vBnd[iSec][iReg].getNodesCoord(reg.name, yIdx):
+ yv.push_back(val)
+ zv = tbox.std_vector_double()
+ for val in self.vBnd[iSec][iReg].getNodesCoord(reg.name, zIdx):
+ zv.push_back(val)
+
+ # Set mesh
+ reg.setMesh(xv,
+ yv,
+ zv,
+ chord,
+ xMin,
+ nodeRows,
+ connectElems)
# External variables
for i, reg in enumerate(self.sec[iSec]):
if reg.name == 'wake':
@@ -172,7 +295,9 @@ class SolversInterface:
loc[iPoint] = reg.loc[iPoint]
self.xxExt[iSec][i] = loc
self.deltaStarExt[iSec][i] = np.zeros(reg.getnNodes())
+ self.vtExt[iSec][i] = np.zeros(reg.getnNodes())
self.vinit = True
+ return 0
def getVWing(self):
"""Obtain the nodes' location and row and cg of all elements.
@@ -256,4 +381,82 @@ class SolversInterface:
self.vBnd[iy][iReg].initStructures(nodesSec.shape[0])
self.vBnd[iy][iReg].nodesCoord = nodesSec
- self.vBnd[iy][iReg].elemsCoord = cgSec
\ No newline at end of file
+ self.vBnd[iy][iReg].elemsCoord = cgSec
+
+ # Inviscid solver methods
+ def run()->int:
+ """Run inviscid solver
+ """
+ raise NotImplementedError('SolverInterface - run not implemented')
+
+ def reset()->None:
+ """Reset solution
+ """
+ raise NotImplementedError('SolverInterface - reset not implemented')
+
+ def save(self)->None:
+ """Save solution
+ """
+ raise NotImplementedError('SolverInterface - save not implemented')
+
+ def writeCp(self)->None:
+ """Write pressure coefficient
+ """
+ raise NotImplementedError('SolverInterface - writeCp not implemented')
+
+ # Inviscid solver properties
+ def getVerbose(self)->int:
+ """Get verbose
+ """
+ raise NotImplementedError('SolverInterface - getVerbose not implemented')
+
+ # Flow properties and aerodynamic coefficients
+ def getAoA(self)->float:
+ """Get angle of attack
+ """
+ raise NotImplementedError('SolverInterface - getAoA not implemented')
+
+ def getCl(self)->float:
+ """Get lift coefficient
+ """
+ raise NotImplementedError('SolverInterface - getCl not implemented')
+
+ def getCd(self)->float:
+ """Get inviscid drag coefficient
+ """
+ raise NotImplementedError('SolverInterface - getCd not implemented')
+
+ def getCm(self)->float:
+ """Get moment coefficient
+ """
+ raise NotImplementedError('SolverInterface - getCm not implemented')
+
+ def getnDim(self)->int:
+ """Get number of dimensions
+ """
+ raise NotImplementedError('SolverInterface - getnDim not implemented')
+
+ def getWing(self)->None:
+ """Get wing mesh
+ """
+ raise NotImplementedError('SolverInterface - getWing not implemented')
+
+ def getMesh(self)->list:
+ """Get mesh
+ """
+ raise NotImplementedError('SolverInterface - getMesh not implemented')
+
+ def getInviscidBC(self)->None:
+ """Get inviscid boundary condition
+ """
+ raise NotImplementedError('SolverInterface - getInviscidBC not implemented')
+
+ def setBlowingVelocity(self)->None:
+ """Set blowing velocity
+ """
+ raise NotImplementedError('SolverInterface - setBlowingVelocity not implemented')
+
+ def getnNodesDomain(self)->int:
+ """Get number of nodes in domain
+ """
+ raise NotImplementedError('SolverInterface - getnNodesDomain not implemented')
diff --git a/blast/interfaces/dart/DartInterface.py b/blast/interfaces/dart/DartInterface.py
index 09e06bead26a75a8cdc60298f0c38a097c67f004..7913ab549b1b745aa332d945590932b7f0633bca 100644
--- a/blast/interfaces/dart/DartInterface.py
+++ b/blast/interfaces/dart/DartInterface.py
@@ -28,11 +28,11 @@ class DartInterface(SolversInterface):
def __init__(self, dartCfg, vSolver, _cfg):
try:
from modules.dartflo.dart.api.core import init_dart
- self.inviscidObjects = init_dart(dartCfg, task=dartCfg['task'], viscous=True)
- self.solver = self.inviscidObjects.get('sol') # Solver
- self.blw = [self.inviscidObjects.get('blwb'), self.inviscidObjects.get('blww')]
+ self.iobj = init_dart(dartCfg, task=dartCfg['task'], viscous=True)
+ self.solver = self.iobj.get('sol') # Solver
+ self.blw = [self.iobj.get('blwb'), self.iobj.get('blww')]
if dartCfg['task'] == 'optimization':
- self.adjointSolver = self.inviscidObjects.get('adj')
+ self.adjointSolver = self.iobj.get('adj')
print('Dart successfully loaded for optimization.')
else:
print('Dart successfully loaded for analysis.')
@@ -76,25 +76,27 @@ class DartInterface(SolversInterface):
for j in range(size):
data[i,3+j] = vData[vElems[i%vElems.size()].nodes[0].row][j]
i += 1
- print('writing data file ' + 'Cp' +sfx + '.dat')
+ if self.getVerbose() > 1:
+ print('writing data file ' + 'Cp' +sfx + '.dat')
np.savetxt('Cp'+sfx+'.dat', data, fmt='%1.6e', delimiter=', ', header='x, y, z, Cp', comments='')
elif self.getnDim() == 3:
# Save surface cp
- data = np.zeros((self.inviscidObjects['bnd'].nodes.size(), 4))
- for i, n in enumerate(self.inviscidObjects['bnd'].nodes):
+ data = np.zeros((self.iobj['bnd'].nodes.size(), 4))
+ for i, n in enumerate(self.iobj['bnd'].nodes):
data[i,0] = n.pos[0]
data[i,1] = n.pos[1]
data[i,2] = n.pos[2]
data[i,3] = self.solver.Cp[n.row]
- np.savetxt(self.inviscidObjects['msh'].name+'_Cp_surface'+sfx+'.dat', data, fmt='%1.6e', delimiter=', ', header='x, y, z, Cp', comments='')
+ np.savetxt(self.iobj['msh'].name+'_Cp_surface'+sfx+'.dat', data, fmt='%1.6e', delimiter=', ', header='x, y, z, Cp', comments='')
import modules.dartflo.dart.utils as invUtils
if self.cfg['Format'] == 'vtk':
- print('Saving Cp files in vtk format. Msh {:s}, Tag {:.0f}'.format(self.inviscidObjects['msh'].name, self.cfg['saveTag']))
- invUtils.write_slices(self.inviscidObjects['msh'].name, self.cfg['writeSections'], self.cfg['saveTag'], sfx=sfx)
+ if self.getVerbose() > 1:
+ print('Saving Cp files in vtk format. Msh {:s}, Tag {:.0f}'.format(self.iobj['msh'].name, self.cfg['saveTag']))
+ invUtils.write_slices(self.iobj['msh'].name, self.cfg['writeSections'], self.cfg['saveTag'], sfx=sfx)
def save(self, sfx='inviscid'):
- self.solver.save(self.inviscidObjects['wrt'], sfx)
+ self.solver.save(self.iobj['wrt'], sfx)
def getAoA(self):
return self.solver.pbl.alpha
@@ -118,7 +120,7 @@ class DartInterface(SolversInterface):
return self.solver.verbose
def getnDim(self):
- return self.inviscidObjects['pbl'].nDim
+ return self.iobj['pbl'].nDim
def getnNodesDomain(self):
return self.solver.pbl.msh.nodes.size()
@@ -154,17 +156,27 @@ class DartInterface(SolversInterface):
for i, blw in enumerate(self.iBnd[n].blowingVel):
self.blw[n].setU(i, blw)
- def setMesh(self):
+ def getWing(self):
if self.getnDim() == 2:
- self.getAirfoil()
+ self.__getWing2D()
elif self.getnDim() == 3:
- self.getWing()
+ self.__getWing3D()
else:
raise RuntimeError('dartInterface::Could not resolve how to set\
the mesh. Either ''nDim'' is incorrect or ''msh'' was not set for\
3D computations')
+
+ def getMesh(self):
+ """Return the mesh on blowing boundaries
+ """
+ meshList = [np.zeros((b.nodes.size(), self.getnDim())) for b in self.blw]
+ for ib, b in enumerate(self.blw):
+ for inod, n in enumerate(b.nodes):
+ for idim in range(self.getnDim()):
+ meshList[ib][inod, idim] = n.pos[idim]
+ return meshList
- def getAirfoil(self):
+ def __getWing2D(self):
"""Create data structure for information transfer.
"""
self.iBnd[0].initStructures(self.blw[0].nodes.size())
@@ -266,9 +278,14 @@ class DartInterface(SolversInterface):
data[:,2] = data[connectListNodes,2]
data[:,3] = data[connectListNodes,3]
data[:,4] = data[connectListNodes,4]
+
+ connect2 = np.zeros((len(connectListElems)), dtype=int)
+ for i, e in enumerate(self.blw[0].tag.elems):
+ connect2[connectListElems[i]] = i
+
self.iBnd[0].nodesCoord = np.column_stack((data[:,1], data[:,2],\
data[:,3], data[:,4]))
- self.iBnd[0].setConnectList(connectListNodes, connectListElems)
+ self.iBnd[0].setConnectList(connectListNodes, connect2)
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])
@@ -277,7 +294,6 @@ class DartInterface(SolversInterface):
elemCoord[i,3] = i
self.iBnd[0].elemsCoord = elemCoord
-
# Wake
self.iBnd[1].initStructures(self.blw[1].nodes.size())
# Node number
@@ -314,7 +330,11 @@ class DartInterface(SolversInterface):
dataW[:,4] = dataW[connectListNodes,4]
self.iBnd[1].nodesCoord = np.column_stack((dataW[:,1], dataW[:,2], \
dataW[:,3], dataW[:,4]))
- self.iBnd[1].setConnectList(connectListNodes, connectListElems)
+
+ connect2_w = np.zeros((len(connectListElems)), dtype=int)
+ for i, e in enumerate(self.blw[1].tag.elems):
+ connect2_w[connectListElems[i]] = i
+ self.iBnd[1].setConnectList(connectListNodes, connect2_w)
elemCoordW = np.zeros((len(self.iBnd[1].nodesCoord)-1, 4))
for i in range(len(elemCoordW[:,0])):
@@ -324,7 +344,7 @@ class DartInterface(SolversInterface):
elemCoordW[i,3] = i
self.iBnd[1].elemsCoord = elemCoordW
- def getWing(self):
+ def __getWing3D(self):
"""Obtain the nodes' location and row and cg of all elements.
"""
diff --git a/blast/mdao/naca0012_mdao.py b/blast/mdao/naca0012_mdao.py
new file mode 100644
index 0000000000000000000000000000000000000000..ec34a6ead36e2ad5799ead7bb2584b3b4f2385ba
--- /dev/null
+++ b/blast/mdao/naca0012_mdao.py
@@ -0,0 +1,201 @@
+import blast.utils as viscUtils
+import numpy as np
+
+from fwk.wutils import parseargs
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+from blast.api.mdaAPI import BlasterSolver
+
+import openmdao.api as om
+
+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.dirname(os.path.abspath(__file__))) + '/models/dart/n0012.geo', # Input file containing the model
+ 'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'msF' : 30, 'msTe' : 0.01, 'msLe' : 0.0075}, # 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': 50, # Maximum number of iterations
+ 'couplTol' : 1e-3, # 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' : [0.05, 0.05], # Forced transition location
+ 'interpolator' : 'Matching',# Interpolator for the coupling
+ 'nSections' : 1, # Number of sections in the domain
+ 'span' : 1.0 # Wing span
+ }
+
+class groupMDA(om.Group):
+ def setup(self):
+ args = parseargs()
+ icfg = cfgInviscid(args.k, args.verb)
+ vcfg = cfgBlast(args.verb)
+
+ #cycle = self.add_subsystem('cycle', om.Group())
+
+ # Add the MDAO components
+ blaster = BlasterSolver(vcfg, icfg)
+
+ # Design variables box
+ self.add_subsystem('dvs', om.IndepVarComp(), promotes=['*'])
+ self.add_subsystem('mesh', blaster.getBoundaryMesh())
+ self.add_subsystem('geometry', blaster.getBoundaryGeometry())
+ self.add_subsystem('solver', blaster)
+
+ # Add objective function
+ self.add_subsystem('obj_cmp', om.ExecComp('obj = cd'), promotes=['obj', 'cd'])
+
+ # Connect meshes
+ self.connect('mesh.x_aero0', 'geometry.x_aero0')
+ self.connect('geometry.x_aero_out', 'solver.x_aero')
+
+ #nonlinear_solver = om.NewtonSolver(solve_subsystems=False)
+ self.linear_solver = om.DirectSolver()
+
+ def configure(self):
+ # Design vars
+ self.dvs.add_output('aoa', val=0.0)
+ self.connect('aoa', 'solver.aoa')
+ #self.add_design_var('aoa', lower=0.0, upper=3.0)
+
+ self.dvs.add_output('tau', val=0.12)
+ self.connect('tau', 'geometry.tau')
+ #self.add_design_var('tau', lower=0.10, upper=0.2)
+
+ self.dvs.add_output('eps', val=0.0)
+ self.connect('eps', 'geometry.eps')
+ self.add_design_var('eps', lower=0.0, upper=0.05)
+
+ self.dvs.add_output('p', val=0.4)
+ self.connect('p', 'geometry.p')
+ self.add_design_var('p', lower=0.35, upper=0.45)
+
+ self.add_constraint('solver.cl', equals=0.4)
+
+ #self.connect('solver.cl', 'cl')
+ self.connect('solver.cd', 'cd')
+ self.add_objective('obj')
+
+def main():
+ # Timer.
+ tms = fwk.Timers()
+ tms['total'].start()
+
+ prob = om.Problem()
+ prob.model = groupMDA()
+
+ prob.driver = om.ScipyOptimizeDriver()
+ prob.driver.options['optimizer'] = 'SLSQP'
+ prob.driver.options['tol'] = 1e-3
+
+ recorder = om.SqliteRecorder('reports/BlasterSolver_file.sql')
+
+ prob.driver.add_recorder(recorder)
+ prob.driver.recording_options['record_objectives'] = True
+ prob.driver.recording_options['record_desvars'] = True
+ prob.driver.recording_options['includes'] = ['aoa', 'solver.cl', 'solver.cd', 'mesh.x_aero0', 'geometry.x_aero_out']
+
+ prob.setup()
+ #prob.check_partials(includes=['cl', 'aoa'], compact_print=True)
+ prob.run_driver()
+
+ # Instantiate your CaseReader
+ cr = om.CaseReader(prob.get_reports_dir() / "../BlasterSolver_file.sql")
+ # Isolate "problem" as your source
+ driver_cases = cr.list_cases('driver', out_stream=None)
+
+ # Get evolution of f and x
+ aoa_evol = []
+ cl_evol = []
+ cd_evol = []
+ shape = []
+ shape0 = []
+
+ for c in driver_cases:
+ case = cr.get_case(c)
+ aoa_evol.append(case.get_val('aoa')[0]*180/np.pi)
+ cl_evol.append(case.get_val('solver.cl')[0])
+ cd_evol.append(case.get_val('solver.cd')[0])
+ shape0.append(case.get_val('mesh.x_aero0'))
+ shape.append(case.get_val('geometry.x_aero_out'))
+
+ print('success')
+ print('aoa:', prob['aoa']*np.pi/180)
+ print('cl:', prob['solver.cl'])
+ print('cd:', prob['solver.cd'])
+ print('obj:', prob['obj'])
+
+ # Recover selig format
+ shape0_selig = np.array(shape0[0])
+ shape1_selig = np.array(shape[-1])
+ for i, val in enumerate(prob.model.solver.isol.vBnd[0][0].connectListNodes):
+ for j in range(prob.model.solver.isol.getnDim()):
+ shape0_selig[i*prob.model.solver.isol.getnDim()+j] = shape0[0][val*prob.model.solver.isol.getnDim()+j]
+ shape1_selig[i*prob.model.solver.isol.getnDim()+j] = shape[-1][val*prob.model.solver.isol.getnDim()+j]
+ from matplotlib import pyplot as plt
+
+ plt.plot(shape0_selig[0::2], shape0_selig[1::2], '--', color='grey', label='Initial shape - cl = {:.2f} - cd = {:.5f}'.format(cl_evol[0], cd_evol[0]))
+ plt.plot(shape1_selig[0::2], shape1_selig[1::2], '-', color='blue', label='Final shape - cl = {:.2f} - cd = {:.5f}'.format(cl_evol[-1], cd_evol[-1]))
+ plt.legend()
+ plt.axis('equal')
+ plt.show()
+ quit()
+
+
+
+ plt.plot(shape0[0][::2], shape0[0][1::2], 'x', color='grey')
+ for xx in shape:
+ plt.plot(xx[::2], xx[1::2], 'x', color='blue')
+ plt.axis('equal')
+ plt.show()
+ quit()
+
+
+if __name__ == "__main__":
+ main()
diff --git a/blast/src/DBoundaryLayer.cpp b/blast/src/DBoundaryLayer.cpp
index ad18c9951b515213bb3107d8bd1db034b88aeaba..5f84a6996231bd1bd69a902c5118b11fde9606d0 100644
--- a/blast/src/DBoundaryLayer.cpp
+++ b/blast/src/DBoundaryLayer.cpp
@@ -25,18 +25,51 @@ BoundaryLayer::~BoundaryLayer()
delete closSolver;
}
-void BoundaryLayer::setMesh(std::vector<double> _x, std::vector<double> _y, std::vector<double> _z, double _chord, double xMin)
+void BoundaryLayer::setMesh(std::vector<double> _x, std::vector<double> _y, std::vector<double> _z, double _chord, double _xMin, std::vector<int> const _rows, std::vector<int> const _no)
{
if (_x.size() != _y.size() || _x.size() != _z.size())
throw std::runtime_error("Mesh coordinates are not consistent.");
nNodes = _x.size();
nElms = nNodes - 1;
+ if (_rows.size() == 0)
+ {
+ rows.resize(0);
+ for (size_t i = 0; i < nNodes; ++i)
+ rows.push_back(static_cast<int>(i));
+ }
+ else if (_rows.size() != nNodes)
+ throw std::runtime_error("Row vector is not consistent.");
+ else
+ {
+ rows = _rows;
+ }
+ // Elements
+ if (_no.size() == 0)
+ {
+ no.resize(0);
+ for (size_t i = 0; i < nElms; ++i)
+ no.push_back(static_cast<int>(i));
+ }
+ else if (_no.size() != nElms)
+ throw std::runtime_error("No vector is not consistent.");
+ else
+ {
+ no = _no;
+ }
+
x = _x;
y = _y;
z = _z;
chord = _chord;
+ xMin = _xMin;
// Compute the local coordinate.
+ this->computeLocalCoordinates();
+ this->reset();
+}
+
+void BoundaryLayer::computeLocalCoordinates()
+{
loc.resize(nNodes, 0.);
alpha.resize(nElms, 0.);
dx.resize(nElms, 0.);
@@ -55,8 +88,6 @@ void BoundaryLayer::setMesh(std::vector<double> _x, std::vector<double> _y, std:
xoc.resize(nNodes, 0.);
for (size_t iPoint = 0; iPoint < nNodes; ++iPoint)
xoc[iPoint] = (x[iPoint] - xMin) / chord;
-
- this->reset();
}
void BoundaryLayer::reset()
@@ -94,7 +125,7 @@ void BoundaryLayer::reset()
*/
void BoundaryLayer::setStagnationBC(double Re)
{
- double dv0 = (vt[1] - vt[0]) / (loc[1] - loc[0]);
+ double dv0 = std::abs((vt[1] - vt[0]) / (loc[1] - loc[0]));
if (dv0 > 0)
u[0] = sqrt(0.075 / (Re * dv0)); // Theta
else
@@ -173,11 +204,9 @@ void BoundaryLayer::setWakeBC(double Re, std::vector<double> UpperCond, std::vec
*/
void BoundaryLayer::computeBlowingVelocity()
{
- //deltaStar[0] = u[0] * u[1];
blowingVelocity.resize(nNodes - 1, 0.);
for (size_t i = 1; i < nNodes; ++i)
{
- //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]))
{
diff --git a/blast/src/DBoundaryLayer.h b/blast/src/DBoundaryLayer.h
index 030bb1d2d40e6015c6745b7bbc93a9a60b4aa3a7..5f7adb2b930c5ac2fdd6696da7ab1b975732c8b7 100644
--- a/blast/src/DBoundaryLayer.h
+++ b/blast/src/DBoundaryLayer.h
@@ -26,16 +26,19 @@ public:
std::string name; ///< Name of the region.
// Mesh
- size_t nNodes; ///< Number of points in the domain.
- size_t nElms; ///< Number of cells in the domain.
- std::vector<double> x; ///< x coordinate in the global frame of reference.
- std::vector<double> y; ///< y coordinate in the global frame of reference.
- std::vector<double> z; ///< z coordinate in the global frame of reference.
- std::vector<double> xoc; ///< x/c coordinate in the global frame of reference.
- std::vector<double> loc; ///< xi coordinate in the local frame of reference.
- std::vector<double> dx; ///< Cell size.
- std::vector<double> alpha; ///< Angle of the cell wrt the x axis of the global frame of reference.
- double chord; ///< Chord of the section.
+ size_t nNodes; ///< Number of points in the domain.
+ size_t nElms; ///< Number of cells in the domain.
+ std::vector<double> x; ///< x coordinate in the global frame of reference.
+ std::vector<double> y; ///< y coordinate in the global frame of reference.
+ std::vector<double> z; ///< z coordinate in the global frame of reference.
+ std::vector<double> xoc; ///< x/c coordinate in the global frame of reference.
+ std::vector<double> loc; ///< xi coordinate in the local frame of reference.
+ std::vector<int> rows; ///< Reference numbers of the nodes.
+ std::vector<int> no; ///< Reference numbers of the nodes.
+ std::vector<double> dx; ///< Cell size.
+ std::vector<double> alpha; ///< Angle of the cell wrt the x axis of the global frame of reference.
+ double chord; ///< Chord of the section.
+ double xMin; ///< Minimum x coordinate of the mesh (for local coordinates computation).
// Boundary layer variables.
std::vector<double> u; ///< Solution vector (theta, H, N, ue, Ct).
@@ -81,7 +84,8 @@ public:
double getMaxMach() const { return *std::max_element(Me.begin(), Me.end()); };
// Setters
- void setMesh(std::vector<double> const _x, std::vector<double> const y, std::vector<double> const z, double const _chord, double const xMin);
+ void setMesh(std::vector<double> const _x, std::vector<double> const y, std::vector<double> const z, double const _chord, double const _xMin, std::vector<int> const _rows=std::vector<int>(0), std::vector<int> const _no=std::vector<int>(0));
+ void computeLocalCoordinates();
void setMe(std::vector<double> const _Me) {if (_Me.size() != nNodes) throw std::runtime_error("Mach number vector is not consistent."); Me = _Me;};
void setVt(std::vector<double> const _vt) {if (_vt.size() != nNodes) throw std::runtime_error("Velocity vector is not consistent."); vt = _vt;};
void setRhoe(std::vector<double> const _rhoe) {if (_rhoe.size() != nNodes) throw std::runtime_error("Density vector is not consistent."); rhoe = _rhoe;};
diff --git a/blast/src/DClosures.cpp b/blast/src/DClosures.cpp
index c411d1c41c3de7937bfdf85b382e138898481999..2cf9f44bc79ed170eb5a9a32bd3a04bb6875c8e8 100644
--- a/blast/src/DClosures.cpp
+++ b/blast/src/DClosures.cpp
@@ -169,8 +169,8 @@ void Closures::turbulentClosures(std::vector<double> &closureVars, double theta,
if (us > 0.95)
us = 0.98;
n = 1.0;
- cf = (1 / Fc) * (0.3 * std::exp(arg) * std::pow(logRt / 2.3026, (-1.74 - 0.31 * Hk)) + 0.00011 * (std::tanh(4 - (Hk / 0.875)) - 1));
- Hkc = std::max(Hk - 1 - 18 / Ret, 0.01);
+ cf = (1 / Fc) * (0.3 * std::exp(arg) * std::pow(logRt / 2.3026, (-1.74 - 0.31 * Hk)) + 0.00011 * (std::tanh(4. - (Hk / 0.875)) - 1.));
+ Hkc = std::max(Hk - 1. - 18. / Ret, 0.01);
// Dissipation coefficient.
Hmin = 1 + 2.1 / std::log(Ret);
Fl = (Hk - 1) / (Hmin - 1);
diff --git a/blast/src/DCoupledAdjoint.cpp b/blast/src/DCoupledAdjoint.cpp
index 8139d32f12a89a518de264b9377b4f9cb1522243..8ea9318add45b1a8ec63c746970f578903e819fb 100644
--- a/blast/src/DCoupledAdjoint.cpp
+++ b/blast/src/DCoupledAdjoint.cpp
@@ -15,6 +15,7 @@
*/
#include "DCoupledAdjoint.h"
+#include "DDriver.h"
#include "wAdjoint.h"
#include "wNewton.h"
#include "wSolver.h"
@@ -23,9 +24,12 @@
#include "wFluid.h"
#include "wMshData.h"
#include "wNode.h"
+#include "wBlowingResidual.h"
#include "wLinearSolver.h"
#include "wGmres.h"
+#include "wCache.h"
+#include "wGauss.h"
#include <Eigen/Sparse>
#include<Eigen/SparseLU>
@@ -35,6 +39,7 @@
#include <deque>
#include <algorithm>
#include <iomanip>
+#include <chrono>
using namespace blast;
@@ -58,53 +63,74 @@ CoupledAdjoint::CoupledAdjoint(std::shared_ptr<dart::Adjoint> _iadjoint, std::sh
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
+ size_t nNs = 0;
+ for (auto bBC : isol->pbl->bBCs)
+ nNs += bBC->nodes.size(); // Number of surface nodes
+ nNs += 1; // Duplicate stagnation point
+ size_t nEs = 0;
+ for (auto bBC : isol->pbl->bBCs)
+ nEs += bBC->uE.size(); // Number of blowing elements
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);
+ dRdStar_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(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);
+ dRdStar_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(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);
+ dRdStar_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(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);
+ dRdStar_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim*nNs);
dRblw_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nDim*nNs);
+ dRphi_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nNs);
+ dRM_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
+ dRrho_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
+ dRv_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNs, nNs);
+ dRdStar_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
+ dRblw_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, 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);
+ dRdStar_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(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);
+ dRdStar_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(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();
+
+ // Gradient wrt objective function
+ dRphi_AoA = Eigen::VectorXd::Zero(nNd);
// Objective functions
+ // Cl = f(phi) & Cd = Cdp + Cdf = f(phi, M, v)
dCl_phi = Eigen::VectorXd::Zero(nNd);
dCd_phi = Eigen::VectorXd::Zero(nNd);
+ dCd_M = Eigen::VectorXd::Zero(nNs);
+ dCd_v = Eigen::VectorXd::Zero(nDim*nNs);
+ dCd_dStar = Eigen::VectorXd::Zero(nNs);
// Angle of attack
dCl_AoA = 0.0;
@@ -115,7 +141,8 @@ void CoupledAdjoint::reset()
// 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);
+ dCdp_x = Eigen::VectorXd::Zero(nDim*nNd);
+ dCdf_x = Eigen::VectorXd::Zero(nDim*nNd);
tdCl_x.resize(nNd, Eigen::Vector3d::Zero());
tdCd_x.resize(nNd, Eigen::Vector3d::Zero());
@@ -165,22 +192,43 @@ void CoupledAdjoint::buildAdjointMatrix(std::vector<std::vector<Eigen::SparseMat
void CoupledAdjoint::run()
{
+ tms["0-adj_total"].start();
+ tms["1-adj_inviscid"].start();
gradientswrtInviscidFlow();
+ tms["1-adj_inviscid"].stop();
+ tms["2-adj_mach"].start();
+ gradientswrtMachNumber();
+ tms["2-adj_mach"].stop();
+ tms["3-adj_density"].start();
+ gradientswrtDensity();
+ tms["3-adj_density"].stop();
+ tms["4-adj_velocity"].start();
+ gradientswrtVelocity();
+ tms["4-adj_velocity"].stop();
+ tms["5-adj_dStar"].start();
+ gradientswrtDeltaStar();
+ tms["5-adj_dStar"].stop();
+ tms["6-adj_blw"].start();
gradientswrtBlowingVelocity();
+ tms["6-adj_blw"].stop();
+ tms["7-adj_aoa"].start();
gradientwrtAoA();
+ tms["7-adj_aoa"].stop();
+ tms["8-adj_mesh"].start();
gradientwrtMesh();
+ tms["8-adj_mesh"].stop();
transposeMatrices();
- // Adjoint matrix
-
+ tms["9-build"].start();
// 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}
+ {&dRphi_phi, &dRM_phi, &dRrho_phi, &dRv_phi, &dRdStar_phi, &dRblw_phi},
+ {&dRphi_M, &dRM_M, &dRrho_M, &dRv_M, &dRdStar_M, &dRblw_M},
+ {&dRphi_rho, &dRM_rho, &dRrho_rho, &dRv_rho, &dRdStar_rho, &dRblw_rho},
+ {&dRphi_v, &dRM_v, &dRrho_v, &dRv_v, &dRdStar_v, &dRblw_v},
+ {&dRphi_dStar, &dRM_dStar, &dRrho_dStar, &dRv_dStar, &dRdStar_dStar, &dRblw_dStar},
+ {&dRphi_blw, &dRM_blw, &dRrho_blw, &dRv_blw, &dRdStar_blw, &dRblw_blw}
};
int rows = 0; int cols = 0;
@@ -188,60 +236,118 @@ void CoupledAdjoint::run()
for (const auto& mat1 : matrices[0]) cols += mat1->cols(); // All matrices in a column have the same number of columns
Eigen::SparseMatrix<double, Eigen::RowMajor> A_adjoint(rows, cols);
buildAdjointMatrix(matrices, A_adjoint);
+ tms["9-build"].stop();
+
+ size_t phiIdx = 0;
+ size_t machIdx = dRphi_phi.cols();
+ size_t rhoIdx = machIdx + dRM_M.cols();
+ size_t vIdx = rhoIdx + dRrho_rho.cols();
+ size_t dStarIdx = vIdx + dRv_v.cols();
+ size_t blwIdx = dStarIdx + dRdStar_dStar.cols();
//***** Gradient wrt angle of attack *****//
//****************************************//
// CL
+ tms["10-solveCxAoA"].start();
std::vector<double> rhsCl(A_adjoint.cols(), 0.0);
- for (auto i = 0; i<dCl_phi.size(); ++i)
- rhsCl[i] = dCl_phi[i];
+ for (auto i = phiIdx; i<phiIdx+dCl_phi.size(); ++i)
+ rhsCl[i] = dCl_phi[i]; // Only dCl_dphi is non-zero
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;
+ tdCl_AoA = dCl_AoA - dRphi_AoA.transpose()*lambdaCl.segment(phiIdx, dRphi_phi.cols());
// 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);
+ std::vector<double> rhsCd(A_adjoint.cols(), 0.0);
+ int jj = 0;
+ for (auto i = phiIdx; i<phiIdx+dCd_phi.size(); ++i)
+ {
+ rhsCd[i] = dCd_phi[jj];
+ jj++;
+ }
+ jj = 0;
+ for (auto i = machIdx; i<machIdx+dCd_M.size(); ++i)
+ {
+ rhsCd[i] = dCd_M[jj];
+ jj++;
+ }
+ jj = 0;
+ for (auto i = vIdx; i<vIdx+dCd_v.size(); ++i)
+ {
+ rhsCd[i] = dCd_v[jj];
+ jj++;
+ }
+ jj = 0;
+ for (auto i = dStarIdx; i<dStarIdx+dCd_dStar.size(); ++i)
+ {
+ rhsCd[i] = dCd_dStar[jj];
+ jj++;
+ }
+ 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());
- // tdCd_AoA = dCd_AoA - dRphi_AoA.transpose()*lambdaCdPhi; // Total gradient
+ alinsol->compute(A_adjoint, Eigen::Map<Eigen::VectorXd>(rhsCd_.data(), rhsCd_.size()), lambdaCd_);
+
+ // Total gradient
+ tdCd_AoA = dCd_AoA - dRphi_AoA.transpose()*lambdaCd.segment(phiIdx, dRphi_phi.cols()); // Total gradient
+ tms["10-solveCxAoA"].stop();
//***** Gradient wrt mesh coordinates *****//
//*****************************************//
- // CL
- // Solution of (from augmented Lagrangian, eqn d/dx )
- // "dCl_x + dRphi_x * lambdaCl_phi + dRblw_x * lambdaCl_blw + dRx_x * lambdaCl_x"
+ // Ck (cL, cd)
+ // Solution lambdaCl_x of (from augmented Lagrangian, eqn d/dx )
+ // "dCk_x + dRphi_x * lambdaCk_phi + dRdStar_x * lambdaCk_x + dRblw_x * lambdaCk_blw + dRx_x * lambdaCk_x = 0"
+
+ tms["12-solveCxMesh"].start();
+ Eigen::VectorXd rhsCl_x = dCl_x;
+ rhsCl_x -= dRphi_x.transpose() * lambdaCl.segment(phiIdx, dRphi_phi.cols()); // Potential contribution
+ rhsCl_x -= dRM_x.transpose() * lambdaCl.segment(machIdx, dRM_M.cols()); // Mach number contribution
+ rhsCl_x -= dRrho_x.transpose() * lambdaCl.segment(rhoIdx, dRrho_rho.cols()); // Density contribution
+ rhsCl_x -= dRv_x.transpose() * lambdaCl.segment(vIdx, dRv_v.cols()); // Velocity contribution
+ rhsCl_x -= dRdStar_x.transpose() * lambdaCl.segment(dStarIdx, dRdStar_dStar.cols()); // Displacement thickness contribution
+ rhsCl_x -= dRblw_x.transpose() * lambdaCl.segment(blwIdx, dRblw_blw.cols()); // Blowing contribution
+
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_);
+ Eigen::VectorXd rhsCd_x = dCdp_x; // Pressure drag coefficient contribution
+ rhsCd_x += dCdf_x; // Friction drag coefficient contribution
+ rhsCd_x -= dRphi_x.transpose() * lambdaCd.segment(phiIdx, dRphi_phi.cols()); // Potential contribution
+ rhsCd_x -= dRM_x.transpose() * lambdaCd.segment(machIdx, dRM_M.cols()); // Mach number contribution
+ rhsCd_x -= dRrho_x.transpose() * lambdaCd.segment(rhoIdx, dRrho_rho.cols()); // Density contribution
+ rhsCd_x -= dRv_x.transpose() * lambdaCd.segment(vIdx, dRv_v.cols()); // Velocity contribution
+ rhsCd_x -= dRdStar_x.transpose() * lambdaCd.segment(dStarIdx, dRdStar_dStar.cols()); // Displacement thickness contribution
+ rhsCd_x -= dRblw_x.transpose() * lambdaCd.segment(blwIdx, dRblw_blw.cols()); // Blowing contribution
+
+ Eigen::VectorXd lambdaCd_x = Eigen::VectorXd::Zero(isol->pbl->nDim * isol->pbl->msh->nodes.size());
+ Eigen::Map<Eigen::VectorXd> lambdaCd_x_(lambdaCd_x.data(), lambdaCd_x.size());
+ alinsol->compute(dRx_x.transpose(), Eigen::Map<Eigen::VectorXd>(rhsCd_x.data(), rhsCd_x.size()), lambdaCd_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];
+ tdCd_x[n->row](m) = lambdaCd_x[isol->pbl->nDim * (n->row) + m];
}
}
+ tms["12-solveCxMesh"].stop();
+ //std::cout << tms << std::endl;
}
+/**
+ * @brief Compute all the gradients wrt the inviscid flow (phi) in the domain
+ * Non-zero are: dRphi_phi, dRM_phi, dRrho_phi, dRv_phi
+ * Zeros are: dRdStar_phi, dRblw_phi, dRx_phi
+ */
void CoupledAdjoint::gradientswrtInviscidFlow()
{
//**** dRphi_phi ****//
@@ -251,7 +357,7 @@ void CoupledAdjoint::gradientswrtInviscidFlow()
auto pbl = isol->pbl;
// Finite difference
std::vector<double> Phi_save = isol->phi; // Differential of the independent variable (phi)
- double deltaPhi = 1e-6; // perturbation
+ double deltaPhi = 1.e-8; // perturbation
std::vector<Eigen::Triplet<double>> tripletsdM;
std::vector<Eigen::Triplet<double>> tripletsdrho;
@@ -261,30 +367,39 @@ void CoupledAdjoint::gradientswrtInviscidFlow()
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));
+ size_t jj = 0;
+ for (auto sec: vsol->sections[0]){
+ int iReg = 0;
+ if (sec->name == "wake")
+ iReg = 1;
+ else if (sec->name == "upper" || sec->name == "lower")
+ iReg = 0;
+ else
+ throw std::runtime_error("gradientswrtInviscidFlow: Unknown section name");
+ for (size_t iNode = 0; iNode < sec->nNodes; ++iNode){
+ tbox::Node* srfNode = pbl->bBCs[iReg]->nodes[sec->rows[iNode]];
+ double dM = 0.; double drho = 0.;
+ Eigen::Vector3d dv = Eigen::Vector3d::Zero();
+ // Average of the quantity on the elements adjacent to the considered node.
+ for (auto e : pbl->fluid->neMap[srfNode]){
+ dM += pbl->fluid->mach->eval(*e, Phi_save, 0);
+ drho += pbl->fluid->rho->eval(*e, Phi_save, 0);
+ Eigen::VectorXd grad = e->computeGradient(Phi_save, 0);
+ for (int i = 0; i < grad.size(); ++i)
+ dv(i) += grad(i);
+ }
+ 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());
- jj++;
+ // Form triplets
+ // Minus because M = f(phi) and therefore RM = M - f(phi) = 0
+ tripletsdM.push_back(Eigen::Triplet<double>(jj, n->row, -(dM - isol->M[srfNode->row])/deltaPhi));
+ tripletsdrho.push_back(Eigen::Triplet<double>(jj, n->row, -(drho - isol->rho[srfNode->row])/deltaPhi));
+ for (int i = 0; i < pbl->nDim; ++i)
+ tripletsdv.push_back(Eigen::Triplet<double>(jj*pbl->nDim + i, n->row, -(dv(i) - isol->U[srfNode->row](i))/deltaPhi));
+ jj++;
+ }
}
// Reset phi
Phi_save[n->row] = isol->phi[n->row];
@@ -296,111 +411,619 @@ void CoupledAdjoint::gradientswrtInviscidFlow()
// Remove zeros
dRM_phi.prune(0.); dRrho_phi.prune(0.); dRv_phi.prune(0.);
- //**** dCl_phi, dCd_phi ****//
+ // // Analytical
+ // std::vector<Eigen::Triplet<double>> tripletsdM_an;
+
+ // size_t jj = 0;
+ // int iReg = 0;
+ // for (auto sec : vsol->sections[0])
+ // {
+ // if (sec->name == "wake")
+ // iReg = 1;
+ // else if (sec->name == "upper" || sec->name == "lower")
+ // iReg = 0;
+ // else
+ // throw std::runtime_error("gradientswrtInviscidFlow: Unknown section name");
+ // for (size_t iNode = 0; iNode < sec->nNodes; ++iNode)
+ // {
+ // tbox::Node* srfNode = pbl->bBCs[iReg]->nodes[sec->rows[iNode]];
+ // for (auto e: pbl->fluid->neMap[srfNode])
+ // {
+ // // Get pre-computed values
+ // tbox::Cache &cache = e->getVCache();
+ // tbox::Gauss &gauss = cache.getVGauss();
+ // Eigen::VectorXd grad = Eigen::VectorXd::Zero(e->nodes.size());
+ // for (size_t k = 0; k < gauss.getN(); ++k)
+ // {
+ // // Gauss point and determinant
+ // double wdj = gauss.getW(k) * e->getDetJ(k);
+ // // Shape functions and solution gradient
+ // Eigen::MatrixXd const &dSf = e->getJinv(k) * cache.getDsf(k);
+ // Eigen::VectorXd dPhi = e->computeGradient(isol->phi, k);
+ // grad += isol->pbl->fluid->mach->evalGrad(*e, isol->phi, k) * dSf.transpose() * dPhi;
+ // }
+ // for (size_t k = 0; k < e->nodes.size(); ++k)
+ // {
+ // tripletsdM_an.push_back(Eigen::Triplet<double>(jj, e->nodes[k]->row, -grad(k)/e->nodes.size()));
+ // }
+ // }
+ // ++jj;
+ // }
+ // }
+ // Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_phi_an(dRM_phi.rows(), dRM_phi.cols());
+ // dRM_phi_an.setFromTriplets(tripletsdM_an.begin(), tripletsdM_an.end());
+ // dRM_phi_an.prune(0.);
+ // writeMatrixToFile(dRM_phi_an, "dRM_phi_an.txt");
+ // writeMatrixToFile(dRM_phi, "dRM_phi.txt");
+ // writeMatrixToFile((dRM_phi_an - dRM_phi), "dif.txt");
+ // std::cout << "written" << std::endl;
+ // std::cout << dRM_phi_an - dRM_phi << std::endl;
+
+ // throw;
+
+
+ //**** partials dCl_phi, dCd_phi ****//
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());
+/**
+ * @brief Compute all the gradients wrt the Mach number on the surface
+ * Non-zero are: dRM_M, dRdStar_M
+ * Zeros are: dRphi_M, dRrho_M, dRv_M, dRblw_M, dRx_M
+ */
+void CoupledAdjoint::gradientswrtMachNumber()
+{
+ /**** dRM_M ****/
+ dRM_M.setIdentity();
+
+ //**** dRdStar_M ****//
+ std::deque<Eigen::Triplet<double>> T;
+ double delta = 1e-6;
+ size_t i = 0;
+ double saveM = 0.;
+ std::vector<Eigen::VectorXd> ddStar(2, Eigen::VectorXd::Zero(dRdStar_M.cols()));
+ std::vector<double> dCdf(2, 0.);
+ for (auto sec: vsol->sections[0])
+ for (size_t j = 0; j < sec->nNodes; ++j)
+ {
+ saveM = sec->Me[j];
+
+ sec->Me[j] = saveM + delta;
+ ddStar[0] = this->runViscous();
+ dCdf[0] = vsol->Cdf;
+ sec->Me[j] = saveM - delta;
+ ddStar[1] = this->runViscous();
+ dCdf[1] = vsol->Cdf;
+
+ sec->Me[j] = saveM;
+
+ for (int k = 0; k < ddStar[0].size(); ++k)
+ T.push_back(Eigen::Triplet<double>(k, i, -(ddStar[0][k] - ddStar[1][k])/(2.*delta)));
+ // Collect Cd gradients
+ dCd_M(i) = (dCdf[0] - dCdf[1])/(2.*delta);
+ ++i;
+ }
+ dRdStar_M.setFromTriplets(T.begin(), T.end());
+ dRdStar_M.prune(0.);
}
-void CoupledAdjoint::gradientswrtDensity(){
- // dRphi_dRho = 0
- // dRM_dRho = 0
- // dRrho_dRho =
- // dRv_dRho = 0
- // dRblw_dRho =
+/**
+ * @brief Compute all the gradients wrt the density on the surface
+ * Non-zero are: dRrho_rho, dRblw_rho
+ * Zeros are: dRphi_rho, dRM_rho, dRv_rho, dRdStar_rho, dRx_rho
+ */
+void CoupledAdjoint::gradientswrtDensity()
+{
+ size_t nNs = 0;
+ for (auto bBC : isol->pbl->bBCs)
+ nNs += bBC->nodes.size(); // Number of surface nodes
+ nNs += 1; // Duplicate stagnation point
+ size_t nEs = 0;
+ for (auto bBC : isol->pbl->bBCs)
+ nEs += bBC->uE.size(); // Number of blowing elements
+
+ //**** dRrho_rho ****//
+ dRrho_rho.setIdentity();
+
+ //**** dRblw_rho ****//
+ std::deque<Eigen::Triplet<double>> T_blw;
+ int offSetElms = 0;
+ int offSetNodes = 0;
+ for (const auto& sec: vsol->sections[0])
+ {
+ std::vector<std::vector<double>> grad = sec->evalGradwrtRho();
+ for (size_t i = 0; i < grad.size(); ++i)
+ for (size_t j = 0; j < grad[i].size(); ++j)
+ T_blw.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes, -grad[i][j]));
+ offSetElms += sec->nElms;
+ offSetNodes += sec->nNodes;
+ }
+ dRblw_rho.setFromTriplets(T_blw.begin(), T_blw.end());
+ dRblw_rho.prune(0.);
+}
+
+/**
+ * @brief Compute all the gradients wrt the velocity on the surface
+ * Non-zeros are: dRv_v, dRdStar_v, dRblw_v
+ * Zeros are: dRphi_v, dRM_v, dRrho_v, dRx_v
+ */
+void CoupledAdjoint::gradientswrtVelocity()
+{
+ size_t nDim = isol->pbl->nDim; // Number of dimensions of the problem
+ size_t nNs = 0;
+ for (auto bBC : isol->pbl->bBCs)
+ nNs += bBC->nodes.size(); // Number of surface nodes
+ nNs += 1; // Duplicate stagnation point
+ size_t nEs = 0;
+ for (auto bBC : isol->pbl->bBCs)
+ nEs += bBC->uE.size(); // Number of blowing elements
+
+ /**** dRv_v ****/
+ dRv_v.setIdentity();
+
+ //**** Get velocity****/
+ Eigen::MatrixXd v = Eigen::MatrixXd::Zero(nNs, nDim);
+ size_t i = 0;
+ int iReg = -1;
+ for (const auto &sec: vsol->sections[0])
+ {
+ if (sec->name == "wake")
+ iReg = 1;
+ else if (sec->name == "upper" || sec->name == "lower")
+ iReg = 0;
+ else
+ throw std::runtime_error("Unknown section name");
+ for (size_t j = 0; j < sec->nNodes; ++j)
+ {
+ for (size_t idim = 0; idim < nDim; ++idim)
+ v(i, idim) = isol->U[isol->pbl->bBCs[iReg]->nodes[sec->rows[j]]->row](idim);
+ ++i;
+ }
+ }
+
+ //**** dvt_v (analytical) ****//
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dVt_v(nNs, nNs * nDim);
+ std::deque<Eigen::Triplet<double>> T_vt;
+
+ i = 0;
+ for (const auto& sec : vsol->sections[0]) {
+ for (size_t j = 0; j < sec->nNodes; ++j) {
+ for (size_t idim = 0; idim < nDim; ++idim) {
+ T_vt.push_back(Eigen::Triplet<double>(i, i * nDim + idim, 1 / std::sqrt(v(i, 0) * v(i, 0) + v(i, 1) * v(i, 1)) * v(i, idim)));
+ }
+ ++i;
+ }
+ }
+ dVt_v.setFromTriplets(T_vt.begin(), T_vt.end());
+ dVt_v.prune(0.);
+
+ //**** dRdStar_vt (finite-difference) ****//
+ std::deque<Eigen::Triplet<double>> T;
+ double delta = 1e-6;
+ i = 0;
+ double saveM = 0.;
+ std::vector<Eigen::VectorXd> dvt(2, Eigen::VectorXd::Zero(dRdStar_M.cols()));
+ std::vector<double> dCdf(2, 0.);
+ Eigen::VectorXd dCd_vt = Eigen::VectorXd::Zero(dRM_M.cols());
+ for (auto sec: vsol->sections[0])
+ for (size_t j = 0; j < sec->nNodes; ++j)
+ {
+ saveM = sec->vt[j];
+
+ sec->vt[j] = saveM + delta;
+ dvt[0] = this->runViscous();
+ dCdf[0] = vsol->Cdf;
+ sec->vt[j] = saveM - delta;
+ dvt[1] = this->runViscous();
+ dCdf[1] = vsol->Cdf;
- // dCl_dRho =
- // dCd_dRho =
+ sec->vt[j] = saveM;
+
+ for (int k = 0; k < dvt[0].size(); ++k)
+ T.push_back(Eigen::Triplet<double>(k, i, -(dvt[0][k] - dvt[1][k])/(2.*delta)));
+ dCd_vt(i) = (dCdf[0] - dCdf[1])/(2.*delta);
+ ++i;
+ }
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_vt(dRdStar_M.rows(), dRdStar_M.cols());
+ dRdStar_vt.setFromTriplets(T.begin(), T.end());
+ dRdStar_vt.prune(0.);
+
+ //**** dRdStar_v ****//
+ dRdStar_v = dRdStar_vt * dVt_v;
+ dRdStar_v.prune(0.);
+
+ //**** dRCdf_v ****//
+ dCd_v = dCd_vt.transpose() * dVt_v; // [1xnNs] x [nNs x nNs*nDim] = [1 x nNs*nDim]
+
+ //**** dRblw_vt (analytical) ****//
+ std::deque<Eigen::Triplet<double>> T_blw;
+ int offSetElms = 0;
+ int offSetNodes = 0;
+ for (const auto& sec: vsol->sections[0])
+ {
+ std::vector<std::vector<double>> grad = sec->evalGradwrtVt();
+ for (size_t i = 0; i < grad.size(); ++i)
+ for (size_t j = 0; j < grad[i].size(); ++j)
+ T_blw.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes, -grad[i][j]));
+ offSetElms += sec->nElms;
+ offSetNodes += sec->nNodes;
+ }
+
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_vt(nEs, nNs);
+ dRblw_vt.setFromTriplets(T_blw.begin(), T_blw.end());
+ dRblw_vt.prune(0.);
+
+ //**** dRblw_v ****//
+ dRblw_v = dRblw_vt * dVt_v;
+ dRblw_v.prune(0.);
}
-void CoupledAdjoint::gradientswrtVelocity(){
- // dRphi_dV = 0
- // dRM_dV = 0
- // dRrho_dV = 0
- // dRv_dV =
- // dRblw_dV =
+/**
+ * @brief Compute all the gradients wrt the displacement thickness on the surface
+ * Non-zero are: dRdStar_dStar, dRblw_dStar
+ * Zeros are: dRphi_dStar, dRM_dStar, dRrho_dStar, dRv_dStar, dRx_dStar
+ */
+void CoupledAdjoint::gradientswrtDeltaStar()
+{
+ //**** dRdStar_dStar (finite-difference) ****//
+ std::deque<Eigen::Triplet<double>> T;
+ double delta = 1e-6;
+ double saveDStar = 0.;
+ std::vector<Eigen::VectorXd> ddstar(2, Eigen::VectorXd::Zero(dRdStar_dStar.cols()));
+ std::vector<double> dCdf(2, 0.);
+ size_t i = 0;
+ for (auto sec: vsol->sections[0])
+ for (size_t j = 0; j < sec->nNodes; ++j)
+ {
+ saveDStar = sec->deltaStarExt[j];
+
+ sec->deltaStarExt[j] = saveDStar + delta;
+ ddstar[0] = this->runViscous();
+ dCdf[0] = vsol->Cdf;
+ sec->deltaStarExt[j] = saveDStar - delta;
+ ddstar[1] = this->runViscous();
+ dCdf[1] = vsol->Cdf;
+
+ sec->deltaStarExt[j] = saveDStar;
+
+ for (int k = 0; k < ddstar[0].size(); ++k)
+ T.push_back(Eigen::Triplet<double>(k, i, (ddstar[0][k] - ddstar[1][k])/(2.*delta)));
+ dCd_dStar(i) = (dCdf[0] - dCdf[1])/(2.*delta);
+ ++i;
+ }
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_dStar_dum(dRdStar_dStar.rows(), dRdStar_dStar.cols());
+ // RdStar = dStar - f(dStar) = 0
+ // dRdStar_dStar = 1 - df_dStar (minus is 3 lines below and not in the triplets)
+ dRdStar_dStar_dum.setFromTriplets(T.begin(), T.end());
+ dRdStar_dStar.setIdentity();
+ dRdStar_dStar -= dRdStar_dStar_dum;
+ dRdStar_dStar.prune(0.);
- // dCl_dV =
- // dCd_dV =
+ //**** dRblw_dStar (analytical) ****//
+ std::deque<Eigen::Triplet<double>> T_blw;
+ int offSetElms = 0;
+ int offSetNodes = 0;
+ for (const auto& sec: vsol->sections[0])
+ {
+ std::vector<std::vector<double>> grad = sec->evalGradwrtDeltaStar();
+ for (size_t i = 0; i < grad.size(); ++i)
+ for (size_t j = 0; j < grad[i].size(); ++j)
+ T_blw.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes, -grad[i][j]));
+ offSetElms += sec->nElms;
+ offSetNodes += sec->nNodes;
+ }
+ dRblw_dStar.setFromTriplets(T_blw.begin(), T_blw.end());
+ dRblw_dStar.prune(0.);
}
-void CoupledAdjoint::gradientswrtBlowingVelocity(){
- dRphi_blw = iadjoint->getRu_Blw();
+/**
+ * @brief Compute all the gradients wrt the blowing velocity
+ * Non-zero are: dRphi_blw, dRblw_blw
+ * Zeros are: dRM_blw, dRrho_blw, dRv_blw, dRdStar_blw, dRx_blw
+ */
+void CoupledAdjoint::gradientswrtBlowingVelocity()
+{
+ /**** dRphi_blw ****/
+ std::deque<Eigen::Triplet<double>> T;
+
+ size_t jj = 0;
+ int iReg = 0;
+ for (auto sec: vsol->sections[0]){
+ if (sec->name == "wake")
+ iReg = 1;
+ else if (sec->name == "upper" || sec->name == "lower")
+ iReg = 0;
+ else
+ throw std::runtime_error("Unknown section name");
+ for (size_t iElm = 0; iElm < sec->nElms; ++iElm){
+ auto blowingElement = isol->pbl->bBCs[iReg]->uE[sec->no[iElm]].first;
+ Eigen::VectorXd be = dart::BlowingResidual::buildGradientBlowing(*blowingElement);
+ for (size_t ii = 0; ii < blowingElement->nodes.size(); ii++)
+ {
+ tbox::Node *nodi = blowingElement->nodes[ii];
+ T.push_back(Eigen::Triplet<double>(isol->rows[nodi->row], jj, be(ii)));
+ }
+ ++jj;
+ }
+ }
+ dRphi_blw.setFromTriplets(T.begin(), T.end());
+ dRphi_blw.prune(0.);
+ dRphi_blw.makeCompressed();
- // All others are zero.
+ /**** dRblw_blw ****/
+ dRblw_blw.setIdentity();
}
+/**
+ * @brief Compute all the gradients wrt the angle of attack
+ * Non-zero are: dRphi_AoA
+ * Zeros are: dRM_AoA, dRrho_AoA, dRv_AoA, dRdStar_AoA, dRblw_AoA, dRx_AoA
+ */
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.
}
+/**
+ * @brief Compute all the gradients wrt the mesh coordinates
+ * Non-zero are: dRphi_x, dRM_x, dRrho_x, dRv_x, dRdStar_x, dRblw_x, dRx_x
+ * Zeros are: -
+ */
void CoupledAdjoint::gradientwrtMesh(){
+
+ /**** dCk_dx ****/
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());
+ dCdp_x = Eigen::Map<Eigen::VectorXd>(dCx_x[1].data(), dCx_x[1].size());
+ /**** dRphi_x ****/
dRphi_x = iadjoint->getRu_X();
+ /**** dRx_x ****/
dRx_x = iadjoint->getRx_X();
+
+ /**** dRdStar_x ****/
+ std::deque<Eigen::Triplet<double>> T_dStar_x;
+ double delta = 1e-8;
+ size_t i = 0;
+ int iReg = 0;
+ double saveX = 0.;
+ std::vector<Eigen::VectorXd> ddStar(2, Eigen::VectorXd::Zero(dRdStar_M.cols()));
+ std::vector<double> dCdf(2, 0.);
+ for (auto sec: vsol->sections[0])
+ {
+ if (sec->name == "wake")
+ iReg = 1;
+ else if (sec->name == "upper" || sec->name == "lower")
+ iReg = 0;
+ else
+ throw std::runtime_error("Unknown section name");
+ for (size_t j = 0; j < sec->nNodes; ++j)
+ {
+ int rowj = isol->pbl->bBCs[iReg]->nodes[sec->rows[j]]->row;
+ // x
+ saveX = sec->x[j];
+
+ sec->x[j] = saveX + delta;
+ sec->computeLocalCoordinates();
+ sec->xxExt = sec->loc;
+ ddStar[0] = this->runViscous();
+ dCdf[0] = vsol->Cdf;
+ sec->x[j] = saveX - delta;
+ sec->computeLocalCoordinates();
+ sec->xxExt = sec->loc;
+ ddStar[1] = this->runViscous();
+ dCdf[1] = vsol->Cdf;
+
+ sec->x[j] = saveX;
+
+ for (int k = 0; k < ddStar[0].size(); ++k)
+ T_dStar_x.push_back(Eigen::Triplet<double>(k, rowj*isol->pbl->nDim+0, -(ddStar[0][k] - ddStar[1][k])/(2.*delta)));
+ dCdf_x(rowj*isol->pbl->nDim+0) += (dCdf[0] - dCdf[1])/(2.*delta);
+ // y
+ saveX = sec->y[j];
+
+ sec->y[j] = saveX + delta;
+ sec->computeLocalCoordinates();
+ sec->xxExt = sec->loc;
+ ddStar[0] = this->runViscous();
+ dCdf[0] = vsol->Cdf;
+ sec->y[j] = saveX - delta;
+ sec->computeLocalCoordinates();
+ sec->xxExt = sec->loc;
+ ddStar[1] = this->runViscous();
+ dCdf[1] = vsol->Cdf;
+
+ sec->y[j] = saveX;
+
+ for (int k = 0; k < ddStar[0].size(); ++k)
+ T_dStar_x.push_back(Eigen::Triplet<double>(k, rowj*isol->pbl->nDim+1, -(ddStar[0][k] - ddStar[1][k])/(2.*delta)));
+ dCdf_x(rowj*isol->pbl->nDim+1) += (dCdf[0] - dCdf[1])/(2.*delta);
+ ++i;
+ }
+ }
+ dRdStar_x.setFromTriplets(T_dStar_x.begin(), T_dStar_x.end());
+ dRdStar_x.prune(0.);
+
+ /**** dRblw_x ****/
+ std::deque<Eigen::Triplet<double>> T_blw_x;
+ int offSetElms = 0;
+ for (const auto& sec: vsol->sections[0])
+ {
+ if (sec->name == "wake")
+ iReg = 1;
+ else if (sec->name == "upper" || sec->name == "lower")
+ iReg = 0;
+ else
+ throw std::runtime_error("Unknown section name");
+ std::vector<std::vector<double>> gradX = sec->evalGradwrtX();
+ std::vector<std::vector<double>> gradY = sec->evalGradwrtY();
+ for (size_t i = 0; i < gradX.size(); ++i)
+ for (size_t j = 0; j < gradX[i].size(); ++j)
+ {
+ int rowj = isol->pbl->bBCs[iReg]->nodes[sec->rows[j]]->row;
+ T_blw_x.push_back(Eigen::Triplet<double>(i + offSetElms, rowj*isol->pbl->nDim+0, -gradX[i][j]));
+ T_blw_x.push_back(Eigen::Triplet<double>(i + offSetElms, rowj*isol->pbl->nDim+1, -gradY[i][j]));
+ }
+ offSetElms += sec->nElms;
+ }
+ dRblw_x.setFromTriplets(T_blw_x.begin(), T_blw_x.end());
+ dRblw_x.prune(0.);
+
+ /**** dRM_x, dRrho_x, dRv_x ****/
+ double delta_x = 1e-8;
+ auto pbl = isol->pbl;
+ std::vector<Eigen::Triplet<double>> tripletsdM;
+ std::vector<Eigen::Triplet<double>> tripletsdrho;
+ std::vector<Eigen::Triplet<double>> tripletsdv;
+
+ for (auto n : pbl->msh->nodes){
+ for (int idim = 0; idim < isol->pbl->nDim; ++idim){
+ // Modify node position
+ double save = n->pos[idim];
+ n->pos[idim] = save + delta_x;
+ for (auto e: pbl->msh->elems)
+ if (e->type() == tbox::ElType::LINE2 || e->type() == tbox::ElType::TRI3) e->update();
+
+ // Recompute the quantities on surface nodes.
+ size_t jj = 0;
+ for (auto sec: vsol->sections[0]){
+ int iReg = 0;
+ if (sec->name == "wake")
+ iReg = 1;
+ else if (sec->name == "upper" || sec->name == "lower")
+ iReg = 0;
+ else
+ throw std::runtime_error("gradientswrtInviscidFlow: Unknown section name");
+ for (size_t iNode = 0; iNode < sec->nNodes; ++iNode){
+ tbox::Node* srfNode = pbl->bBCs[iReg]->nodes[sec->rows[iNode]];
+ double dM = 0.; double drho = 0.;
+ Eigen::Vector3d dv = Eigen::Vector3d::Zero();
+ // Average of the quantity on the elements adjacent to the considered node.
+ for (auto e : pbl->fluid->neMap[srfNode]){
+ dM += pbl->fluid->mach->eval(*e, isol->phi, 0);
+ drho += pbl->fluid->rho->eval(*e, isol->phi, 0);
+ Eigen::VectorXd grad = e->computeGradient(isol->phi, 0);
+ for (int i = 0; i < grad.size(); ++i)
+ dv(i) += grad(i);
+ }
+ dM /= static_cast<double>(pbl->fluid->neMap[srfNode].size());
+ drho /= static_cast<double>(pbl->fluid->neMap[srfNode].size());
+ dv /= static_cast<double>(pbl->fluid->neMap[srfNode].size());
+
+ // Form triplets
+ // Minus because M = f(phi) and therefore RM = M - f(phi) = 0
+ tripletsdM.push_back(Eigen::Triplet<double>(jj, n->row*pbl->nDim + idim, -(dM - isol->M[srfNode->row])/delta_x));
+ tripletsdrho.push_back(Eigen::Triplet<double>(jj, n->row*pbl->nDim + idim, -(drho - isol->rho[srfNode->row])/delta_x));
+ for (int i = 0; i < pbl->nDim; ++i)
+ tripletsdv.push_back(Eigen::Triplet<double>(jj*pbl->nDim + i, n->row*pbl->nDim + idim, -(dv(i) - isol->U[srfNode->row](i))/delta_x));
+ jj++;
+ }
+ }
+ // Reset node pos
+ n->pos[idim] = save;
+ }
+ }
+ // Fill matrices with gradients
+ dRM_x.setFromTriplets(tripletsdM.begin(), tripletsdM.end());
+ dRrho_x.setFromTriplets(tripletsdrho.begin(), tripletsdrho.end());
+ dRv_x.setFromTriplets(tripletsdv.begin(), tripletsdv.end());
+ // Remove zeros
+ dRM_x.prune(0.); dRrho_x.prune(0.); dRv_x.prune(0.);
}
-void CoupledAdjoint::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
+void CoupledAdjoint::setdRdStar_M(std::vector<std::vector<double>> _dRdStar_dM){
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]));
+ for (size_t i = 0; i < _dRdStar_dM.size(); i++){
+ for (size_t j = 0; j < _dRdStar_dM[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRdStar_dM[i][j]));
}
}
// Set matrix
- dRblw_M.setFromTriplets(triplets.begin(), triplets.end());
- dRblw_M.prune(0.);
+ dRdStar_M.setFromTriplets(triplets.begin(), triplets.end());
+ dRdStar_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
+
+void CoupledAdjoint::setdRdStar_v(std::vector<std::vector<double>> _dRdStar_dv){
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]));
+ for (size_t i = 0; i < _dRdStar_dv.size(); i++){
+ for (size_t j = 0; j < _dRdStar_dv[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRdStar_dv[i][j]));
}
}
// Set matrix
- dRblw_v.setFromTriplets(triplets.begin(), triplets.end());
- dRblw_v.prune(0.);
+ dRdStar_v.setFromTriplets(triplets.begin(), triplets.end());
+ dRdStar_v.prune(0.);
+}
+
+void CoupledAdjoint::setdRdStar_dStar(std::vector<std::vector<double>> _dRdStar_dStar){
+ std::vector<Eigen::Triplet<double>> triplets;
+ for (size_t i = 0; i < _dRdStar_dStar.size(); i++){
+ for (size_t j = 0; j < _dRdStar_dStar[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRdStar_dStar[i][j]));
+ }
+ }
+ // Set matrix
+ dRdStar_dStar.setFromTriplets(triplets.begin(), triplets.end());
+ dRdStar_dStar.prune(0.);
}
-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
+
+void CoupledAdjoint::setdRdStar_x(std::vector<std::vector<double>> _dRdStar_x){
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]));
+ for (size_t i = 0; i < _dRdStar_x.size(); i++){
+ for (size_t j = 0; j < _dRdStar_x[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRdStar_x[i][j]));
+ }
+ }
+ // Set matrix
+ dRdStar_x.setFromTriplets(triplets.begin(), triplets.end());
+ dRdStar_x.prune(0.);
+}
+
+void CoupledAdjoint::setdRblw_rho(std::vector<std::vector<double>> _dRblw_rho){
+ std::vector<Eigen::Triplet<double>> triplets;
+ for (size_t i = 0; i < _dRblw_rho.size(); i++){
+ for (size_t j = 0; j < _dRblw_rho[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_rho[i][j]));
}
}
// Set matrix
dRblw_rho.setFromTriplets(triplets.begin(), triplets.end());
dRblw_rho.prune(0.);
}
-void CoupledAdjoint::setdRblw_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
+
+void CoupledAdjoint::setdRblw_v(std::vector<std::vector<double>> _dRblw_v){
+ std::vector<Eigen::Triplet<double>> triplets;
+ for (size_t i = 0; i < _dRblw_v.size(); i++){
+ for (size_t j = 0; j < _dRblw_v[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_v[i][j]));
+ }
+ }
+ // Set matrix
+ dRblw_v.setFromTriplets(triplets.begin(), triplets.end());
+ dRblw_v.prune(0.);
+}
+
+void CoupledAdjoint::setdRblw_dStar(std::vector<std::vector<double>> _dRblw_dStar){
+ std::vector<Eigen::Triplet<double>> triplets;
+ for (size_t i = 0; i < _dRblw_dStar.size(); i++){
+ for (size_t j = 0; j < _dRblw_dStar[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_dStar[i][j]));
+ }
+ }
+ // Set matrix
+ dRblw_dStar.setFromTriplets(triplets.begin(), triplets.end());
+ dRblw_dStar.prune(0.);
+}
+
+void CoupledAdjoint::setdRblw_x(std::vector<std::vector<double>> _dRblw_x){
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]));
+ for (size_t i = 0; i < _dRblw_x.size(); i++){
+ for (size_t j = 0; j < _dRblw_x[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_x[i][j]));
}
}
// Set matrix
@@ -408,47 +1031,99 @@ void CoupledAdjoint::setdRblw_x(std::vector<std::vector<double>> _dRblw_dx){
dRblw_x.prune(0.);
}
+Eigen::VectorXd CoupledAdjoint::runViscous()
+{
+ int exitCode = vsol->run();
+
+ // Extract deltaStar
+ std::vector<Eigen::VectorXd> dStar_p;
+
+ // Extract parts in a loop
+ for (size_t i = 0; i < vsol->sections[0].size(); ++i) {
+ std::vector<double> deltaStar_vec = this->vsol->sections[0][i]->deltaStar;
+ Eigen::VectorXd deltaStar_eigen = Eigen::Map<Eigen::VectorXd>(deltaStar_vec.data(), deltaStar_vec.size());
+ dStar_p.push_back(deltaStar_eigen);
+ }
+ // Calculate the total size
+ int nNs = 0;
+ for (const auto& p : dStar_p)
+ nNs += p.size();
+
+ // Concatenate the vectors
+ Eigen::VectorXd deltaStar(nNs);
+ int idx = 0;
+ for (const auto& p : dStar_p) {
+ deltaStar.segment(idx, p.size()) = p;
+ idx += p.size();
+ }
+ return deltaStar;
+}
+
/**
* @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_dStar = dRphi_dStar.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_dStar = dRM_dStar.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_dStar = dRrho_dStar.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_dStar = dRv_dStar.transpose();
dRv_blw = dRv_blw.transpose();
+ dRdStar_phi = dRdStar_phi.transpose();
+ dRdStar_M = dRdStar_M.transpose();
+ dRdStar_rho = dRdStar_rho.transpose();
+ dRdStar_v = dRdStar_v.transpose();
+ dRdStar_dStar = dRdStar_dStar.transpose();
+ dRdStar_blw = dRdStar_blw.transpose();
+
dRblw_phi = dRblw_phi.transpose();
dRblw_M = dRblw_M.transpose();
dRblw_rho = dRblw_rho.transpose();
dRblw_v = dRblw_v.transpose();
+ dRblw_dStar = dRblw_dStar.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()) {
+ // Write the column headers
+ for (int j = 0; j < matrix.cols(); ++j) {
+ file << std::setw(15) << j;
+ if (j != matrix.cols() - 1) {
+ file << " ";
+ }
+ }
+ file << "\n";
+
+ // Write the matrix data
for (int i = 0; i < matrix.rows(); ++i) {
for (int j = 0; j < matrix.cols(); ++j) {
- file << std::setprecision(20) << matrix(i, j);
+ file << std::fixed << std::setprecision(12) << std::setw(15) << matrix(i, j);
if (j != matrix.cols() - 1) {
file << " ";
}
diff --git a/blast/src/DCoupledAdjoint.h b/blast/src/DCoupledAdjoint.h
index 3ae603da7eb523f8da19ef8243770e714620d478..086a70d31cd853d8508bf4d59e111bd4d8dac390 100644
--- a/blast/src/DCoupledAdjoint.h
+++ b/blast/src/DCoupledAdjoint.h
@@ -20,6 +20,7 @@
#include "blast.h"
#include "wObject.h"
#include "dart.h"
+#include "wTimers.h"
#include <Eigen/Sparse>
#include <iostream>
@@ -52,6 +53,7 @@ private:
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_dStar; ///< Derivative of the inviscid flow residual with respect to delta*
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
@@ -59,6 +61,7 @@ private:
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_dStar; ///< Derivative of the Mach number residual with respect to delta*
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
@@ -66,6 +69,7 @@ private:
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_dStar; ///< Derivative of the density residual with respect to delta*
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
@@ -73,29 +77,46 @@ private:
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_dStar; ///< Derivative of the velocity residual with respect to delta*
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> dRdStar_phi; ///< Derivative of delta* residual with respect to the inviscid flow
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_M; ///< Derivative of delta* residual with respect to the Mach number
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_v; ///< Derivative of delta* residual with respect to the velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_rho; ///< Derivative of delta* residual with respect to the density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_dStar; ///< Derivative of delta* residual with respect to delta*
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_blw; ///< Derivative of delta* residual with respect to the blowing velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_x; ///< Derivative of delta* 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_dStar; ///< Derivative of the blowing velocity residual with respect to delta*
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
+ Eigen::VectorXd dCd_M; ///< Derivative of the drag coefficient with respect to the Mach number
+ Eigen::VectorXd dCd_rho; ///< Derivative of the drag coefficient with respect to the density
+ Eigen::VectorXd dCd_v; ///< Derivative of the drag coefficient with respect to the velocity
+ Eigen::VectorXd dCd_dStar; ///< Derivative of the drag coefficient with respect to delta*
// 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
+ double dCd_AoA; ///< Partial derivative of the total 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
+ Eigen::VectorXd dCdp_x; ///< Partial derivative of the pressure drag coefficient with respect to the mesh coordinates
+ Eigen::VectorXd dCdf_x; ///< Partial derivative of the friction drag coefficient with respect to the mesh coordinates
+
+ fwk::Timers tms; //< internal timers
public:
CoupledAdjoint(std::shared_ptr<dart::Adjoint> iAdjoint, std::shared_ptr<blast::Driver> vSolver);
@@ -104,18 +125,26 @@ public:
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);
+ void setdRdStar_M(std::vector<std::vector<double>> _dRdStar_M);
+ void setdRdStar_v(std::vector<std::vector<double>> _dRdStar_v);
+ void setdRdStar_dStar(std::vector<std::vector<double>> _dRdStar_dStar);
+ void setdRdStar_x(std::vector<std::vector<double>> _dRdStar_x);
+
+ void setdRblw_rho(std::vector<std::vector<double>> _dRblw_rho);
+ void setdRblw_v(std::vector<std::vector<double>> _dRblw_v);
+ void setdRblw_dStar(std::vector<std::vector<double>> _dRblw_ddStar);
+ void setdRblw_x(std::vector<std::vector<double>> _dRblw_dx);
private:
void buildAdjointMatrix(std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor>*>> &matrices, Eigen::SparseMatrix<double, Eigen::RowMajor> &matrixAdjoint);
+ Eigen::VectorXd runViscous();
+
void gradientswrtInviscidFlow();
void gradientswrtMachNumber();
void gradientswrtDensity();
void gradientswrtVelocity();
+ void gradientswrtDeltaStar();
void gradientswrtBlowingVelocity();
void gradientwrtAoA();
diff --git a/blast/src/DDriver.cpp b/blast/src/DDriver.cpp
index 2fee439bf1b7c806d1418c4900c24f19a73fb287..8a4fb03d44511956c40cce0a98e8bfcc4cfca715 100644
--- a/blast/src/DDriver.cpp
+++ b/blast/src/DDriver.cpp
@@ -266,7 +266,7 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced)
double avTurb = (bl->x[iPoint] - bl->xtr) / (bl->x[iPoint] - bl->x[iPoint-1]);
double avLam = 1. - avTurb;
- //std::cout << "avTurb " << avTurb << " avLam " << avLam << std::endl;
+ // std::cout << "avTurb " << avTurb << " avLam " << avLam << std::endl;
if (avTurb < 0. || avTurb > 1. || avLam < 0. || avLam > 1.)
{
bl->printSolution(iPoint);
diff --git a/blast/tests/dart/adjointVII_2D.py b/blast/tests/dart/adjointVII_2D.py
index 7341e9836585a889184a840658e0ea09ebb30db0..ea6db8211ea1d332f2aaf6157edf445efe66570b 100644
--- a/blast/tests/dart/adjointVII_2D.py
+++ b/blast/tests/dart/adjointVII_2D.py
@@ -18,12 +18,11 @@
# @author Paul Dechamps
# @date 2024
-# Test the blast adjoint implementation.
+# Test the blaster adjoint implementation.
# Imports.
import blast.utils as viscUtils
-from blast.interfaces.DAdjointInterface import AdjointInterface
import blast
import numpy as np
@@ -37,8 +36,6 @@ from fwk.testing import *
from fwk.coloring import ccolors
import numpy as np
-import math
-
def cfgInviscid(nthrds, verb):
import os.path
# Parameters
@@ -50,7 +47,7 @@ def cfgInviscid(nthrds, verb):
'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
+ 'Pars' : {'xLgt' : 5, 'yLgt' : 5, 'msF' : 3, 'msTe' : 0.1, 'msLe' : 0.1}, # parameters for input file model
'Dim' : 2, # problem dimension
'Format' : 'gmsh', # save format (vtk or gmsh)
# Markers
@@ -72,7 +69,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' : 'SparseLu', # 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
@@ -87,8 +84,8 @@ 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': 200, # Maximum number of iterations
- 'couplTol' : 1e-8, # Tolerance of the VII methodology
+ 'couplIter': 100, # Maximum number of iterations
+ 'couplTol' : 1e-10, # 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
@@ -108,25 +105,12 @@ def main():
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
+ morpher = isol.iobj['mrf']
# Adjoint objects
cAdj = blast.CoupledAdjoint(isol.adjointSolver, vsol)
- pyAdjoint = AdjointInterface(cAdj, isol, vsol)
# Run coupled case
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
@@ -140,35 +124,103 @@ def main():
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()
+ # Get mesh matrix
+ dCl_x = np.zeros((isol.iobj['bnd'].nodes.size(), isol.getnDim()))
+ dCd_x = np.zeros((isol.iobj['bnd'].nodes.size(), isol.getnDim()))
+ for inod, n in enumerate(isol.blw[0].nodes):
+ for i in range(isol.getnDim()):
+ dCl_x[inod, i] = cAdj.tdCl_x[n.row][i]
+ dCd_x[inod, i] = cAdj.tdCd_x[n.row][i]
+
+
+ print(ccolors.ANSI_BLUE + 'PyFiniteDifference...' + ccolors.ANSI_RESET)
+ # Finite difference on AoA
+ tms['fd_blaster_aoa'].start()
dcl = []
+ dcd = []
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()
+ dcd.append(isol.getCd()+vsol.Cdf)
+ dCl_AoA_FD = (dcl[-1] - dcl[-2])/(2.*dalpha)
+ dCd_AoA_FD = (dcd[-1] - dcd[-2])/(2.*dalpha)
+ tms['fd_blaster_aoa'].stop()
+ isol.setAoA(icfg['AoA']*np.pi/180) # Reset AoA after FD
+
+ # Recover grad AoA from grad x
+ drot = np.array([[-np.sin(alpha), np.cos(alpha)], [-np.cos(alpha), -np.sin(alpha)]])
+
+ dCl_AoA_fromX = 0.
+ dCd_AoA_fromX = 0.
+ for inod, n in enumerate(isol.iobj['bnd'].nodes):
+ dx = drot.dot([n.pos[0] - 0.25, n.pos[1]])
+ dCl_AoA_fromX += np.array([cAdj.tdCl_x[n.row][0], cAdj.tdCl_x[n.row][1]]).dot(dx)
+ dCd_AoA_fromX += np.array([cAdj.tdCd_x[n.row][0], cAdj.tdCd_x[n.row][1]]).dot(dx)
+
+ ### BLASTER FD MESH ###
+ tms['fd_blaster_msh'].start()
+ saveNodes = np.zeros((len(isol.iobj['pbl'].msh.nodes), isol.getnDim()))
+ for n in isol.iobj['pbl'].msh.nodes:
+ for i in range(isol.getnDim()):
+ saveNodes[n.row, i] = n.pos[i]
+ dCl_x_blaster_fd = np.zeros((len(isol.blw[0].nodes), isol.getnDim()))
+ dCd_x_blaster_fd = np.zeros((len(isol.blw[0].nodes), isol.getnDim()))
+ delta = 1e-6
+ cmp = 0
+ import time
+ start_time = time.time()
+ for inod, n in enumerate(isol.blw[0].nodes):
+ # Safety check
+ if isol.iobj['bnd'].nodes[inod].row != n.row:
+ raise RuntimeError('Nodes do not match', n.row, isol.iobj['bnd'].nodes[inod].row)
+ for idim in range(isol.getnDim()):
+ cl = [0, 0]
+ cd = [0, 0]
+ for isgn, sign in enumerate([-1, 1]):
+ morpher.savePos()
+ n.pos[idim] = saveNodes[n.row, idim] + sign * delta
+ morpher.deform()
+ isol.updateMesh()
+
+ coupler.reset()
+ # prevent print
+ from contextlib import redirect_stdout
+ with redirect_stdout(None):
+ coupler.run(write=False)
+ cmp += 1
+ if cmp % 20 == 0:
+ print('F-D opers', cmp, '/', isol.blw[0].nodes.size() * isol.getnDim() * 2, '(node '+str(n.row)+')', 'time', str(round(time.time() - start_time, 3))+'s')
+ cl[isgn] = isol.getCl()
+ cd[isgn] = isol.getCd() + vsol.Cdf
+
+ # Reset mesh
+ for node in isol.iobj['sol'].pbl.msh.nodes:
+ for d in range(isol.getnDim()):
+ node.pos[d] = saveNodes[node.row, d]
+ dCl_x_blaster_fd[inod, idim] = (cl[1] - cl[0]) / (2. * delta)
+ dCd_x_blaster_fd[inod, idim] = (cd[1] - cd[0]) / (2. * delta)
+ tms['fd_blaster_msh'].stop()
+
+ print('Statistics')
print(tms)
+ print(ccolors.ANSI_BLUE + 'PyTesting' + ccolors.ANSI_RESET)
+ tests = CTests()
+ tests.add(CTest('dCl_dAoA', cAdj.tdCl_AoA, dCl_AoA_FD, 1e-5))
+ tests.add(CTest('dCd_dAoA', cAdj.tdCd_AoA, dCd_AoA_FD, 1e-4))
+ tests.add(CTest('norm dCl_dx (adj - fd)', np.linalg.norm(dCl_x), np.linalg.norm(dCl_x_blaster_fd), 1e-6))
+ tests.add(CTest('norm dCd_dx (adj - fd)', np.linalg.norm(dCd_x), np.linalg.norm(dCd_x_blaster_fd), 1e-4))
+ tests.add(CTest('dCl_dAoA (msh)', dCl_AoA_fromX, cAdj.tdCl_AoA, 1e-2)) # The value gets closer if the mesh if refined
+ tests.add(CTest('dCd_dAoA (msh)', dCd_AoA_fromX, cAdj.tdCd_AoA, 1e-1)) # The value gets closer if the mesh if refined
+ tests.run()
+
+ # eof
+ print('')
if __name__ == '__main__':
main()
\ No newline at end of file
diff --git a/blast/utils.py b/blast/utils.py
index 93efffcbdd4cc05225ece4500730d58ef26ca91a..db3d5236f2a249159a775f876f8c5728f0b78602 100644
--- a/blast/utils.py
+++ b/blast/utils.py
@@ -44,7 +44,7 @@ def initBL(Re, Minf, CFL0, nSections, xtrF = [None, None], span=0, verb=None):
if xtrF[i] is None:
xtrF[i] = -1
if xtrF[i] != -1 and not(0<= xtrF[i] <= 1):
- raise RuntimeError('Incorrect forced transition location')
+ raise RuntimeError('Incorrect forced transition location')
solver = blast.Driver(Re, Minf, CFL0, nSections, xtrF_top=xtrF[0], xtrF_bot=xtrF[1], _span=span, _verbose=verbose)
return solver
@@ -70,7 +70,7 @@ def initBlast(iconfig, vconfig, iSolver='DART'):
if 'sfx' not in vconfig:
vconfig['sfx'] = ''
# Viscous solver
- vSolver = initBL(vconfig['Re'], vconfig['Minf'], vconfig['CFL0'], vconfig['nSections'], xtrF=vconfig['xtrF'])
+ vSolver = initBL(vconfig['Re'], vconfig['Minf'], vconfig['CFL0'], vconfig['nSections'], xtrF=vconfig['xtrF'], verb=vconfig['Verb'])
# Inviscid solver
if iSolver == 'DART':
@@ -123,7 +123,7 @@ def getSolution(sections, write=False, toW=['x', 'xelm', 'theta', 'H', 'deltaSta
varKeys = ['theta', 'H', 'N', 'ue', 'Ct']
attrKeys = ['x', 'y', 'z', 'xoc', 'deltaStar', \
'cd', 'cf', 'Hstar', 'Hstar2', 'Hk', 'ctEq', 'us', 'delta',\
- 'vt', 'rhoe', 'Ret']
+ 'vt', 'vtExt', 'deltaStarExt', 'rhoe', 'Ret', 'regime']
elemsKeys = ['xelm', 'yelm', 'zelm']
for iSec, sec in enumerate(sections):
diff --git a/blast/validation/raeValidation.py b/blast/validation/raeValidation.py
index 2ee1444ccf710a15eeda3c39fde086f3b75477e8..ef462b89c1a17f7710be828104fba0cd4dd64599 100644
--- a/blast/validation/raeValidation.py
+++ b/blast/validation/raeValidation.py
@@ -132,12 +132,12 @@ def main():
# Test solution
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
- tests.add(CTest('Cl', isol.getCl(), 0.759, 5e-2))
- tests.add(CTest('Cd wake', vsol.Cdt, 0.0093, 1e-3, forceabs=True))
- tests.add(CTest('Cd integral', isol.getCd() + vsol.Cdf, 0.0134, 1e-3, forceabs=True))
+ tests.add(CTest('Cl', isol.getCl(), 0.730, 5e-2))
+ tests.add(CTest('Cd wake', vsol.Cdt, 0.0095, 1e-3, forceabs=True))
+ tests.add(CTest('Cd integral', isol.getCd() + vsol.Cdf, 0.0132, 1e-3, forceabs=True))
tests.add(CTest('Cdf', vsol.Cdf, 0.0068, 1e-3, forceabs=True))
if icfg['LSolver'] == 'SparseLu':
- tests.add(CTest('Iterations', len(aeroCoeffs['Cl']), 29, 0, forceabs=True))
+ tests.add(CTest('Iterations', len(aeroCoeffs['Cl']), 21, 0, forceabs=True))
tests.run()
expResults = np.loadtxt(os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + '/models/references/rae2822_AR138_case6.dat')
diff --git a/modules/dartflo b/modules/dartflo
index 2cdf73b03d98e1dae2a094fe9f1edb9f4f7e1eab..7bb5ce07ce84f4620798b99c68e2f94181870936 160000
--- a/modules/dartflo
+++ b/modules/dartflo
@@ -1 +1 @@
-Subproject commit 2cdf73b03d98e1dae2a094fe9f1edb9f4f7e1eab
+Subproject commit 7bb5ce07ce84f4620798b99c68e2f94181870936