Commit b8632543 authored by Vaggelis Papoutsis's avatar Vaggelis Papoutsis Committed by Andrew Heather
Browse files

ENH: New adjont shape optimisation functionality

The adjoint library is enhanced with new functionality enabling
automated shape optimisation loops.  A parameterisation scheme based on
volumetric B-Splines is introduced, the control points of which act as
the design variables in the optimisation loop [1, 2].  The control
points of the volumetric B-Splines boxes can be defined in either
Cartesian or cylindrical coordinates.

The entire loop (solution of the flow and adjoint equations, computation
of sensitivity derivatives, update of the design variables and mesh) is
run within adjointOptimisationFoam. A number of methods to update the
design variables are implemented, including popular Quasi-Newton methods
like BFGS and methods capable of handling constraints like loop using
the SQP or constraint projection.

The software was developed by PCOpt/NTUA and FOSS GP, with contributions from

Dr. Evangelos Papoutsis-Kiachagias,
Konstantinos Gkaragounis,
Professor Kyriakos Giannakoglou,
Andy Heather

[1] E.M. Papoutsis-Kiachagias, N. Magoulas, J. Mueller, C. Othmer,
K.C.  Giannakoglou: 'Noise Reduction in Car Aerodynamics using a
Surrogate Objective Function and the Continuous  Adjoint Method with
Wall Functions', Computers & Fluids, 122:223-232, 2015

[2] E. M. Papoutsis-Kiachagias, V. G. Asouti, K. C. Giannakoglou,
K.  Gkagkas, S. Shimokawa, E. Itakura: ‘Multi-point aerodynamic shape
optimization of cars based on continuous adjoint’, Structural and
Multidisciplinary Optimization, 59(2):675–694, 2019
parent eb4fec37
......@@ -55,10 +55,6 @@ int main(int argc, char *argv[])
for (om++; !om.end(); om++)
{
Info<< "* * * * * * * * * * * * * * * * * * *" << endl;
Info<< "Time = " << runTime.timeName() << endl;
Info<< "* * * * * * * * * * * * * * * * * * *" << endl;
if (om.update())
{
// Update design variables and solve all primal equations
......
cumulativeDisplacement.C
EXE = $(FOAM_APPBIN)/cumulativeDisplacement
EXE_INC = \
-I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-ldynamicFvMesh \
-ldynamicMesh \
-lmeshTools \
-lfiniteVolume
pointField points0
(
pointIOField
(
IOobject
(
"points",
mesh.time().constant(),
polyMesh::meshSubDir,
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
)
)
);
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2007-2019 PCOpt/NTUA
Copyright (C) 2013-2019 FOSS GP
Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
cumulativeDisplacement
Description
Computes and writes the difference between the mesh points in each time
instance and the ones in the constant folder. Additionally, the projection
of this difference to the normal point vectors of the initial mesh is also
written
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "emptyFvPatch.H"
#include "coupledFvPatch.H"
#include "pointMesh.H"
#include "pointPatchField.H"
#include "pointPatchFieldsFwd.H"
#include "syncTools.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
timeSelector::addOptions();
#include "addRegionOption.H"
#include "setRootCase.H"
#include "createTime.H"
instantList timeDirs = timeSelector::select0(runTime, args);
#include "createNamedMesh.H"
#include "createFields.H"
// polyPatch::pointNormals will give the wrong result for points
// belonging to multiple patches or patch-processorPatch intersections.
// Keeping a mesh-wide field to allow easy reduction using syncTools.
// A bit expensive? Better way?
vectorField pointNormals(mesh.nPoints(), vector::zero);
for (const fvPatch& patch : mesh.boundary())
{
// Each local patch point belongs to these local patch faces.
// Local numbering
const labelListList& patchPointFaces = patch.patch().pointFaces();
if (!isA<coupledFvPatch>(patch) && !isA<emptyFvPatch>(patch))
{
const labelList& meshPoints = patch.patch().meshPoints();
const vectorField nf(patch.nf());
forAll(meshPoints, ppI)
{
const labelList& pointFaces = patchPointFaces[ppI];
forAll(pointFaces, pfI)
{
const label& localFaceIndex = pointFaces[pfI];
pointNormals[meshPoints[ppI]] += nf[localFaceIndex];
}
}
}
}
// Sum from all processors
syncTools::syncPointList
(
mesh, pointNormals, plusEqOp<vector>(), vector::zero
);
pointNormals /= mag(pointNormals) + VSMALL;
forAll(timeDirs, timeI)
{
runTime.setTime(timeDirs[timeI], timeI);
Info<< "Time = " << runTime.timeName() << endl;
mesh.readUpdate();
const pointMesh& pMesh(pointMesh::New(mesh));
// Point displacement projected to the
// unit point normal of the initial geometry
pointScalarField normalDisplacement
(
IOobject
(
"normalDisplacement",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
pMesh,
dimensionedScalar(dimless, Zero)
);
// Point displacement as a vector
pointVectorField displacement
(
IOobject
(
"displacement",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
pMesh,
dimensionedVector(dimless, Zero)
);
Info<< "Calculating cumulative mesh movement for time "
<< runTime.timeName() << endl;
// Normal displacement
const pointField& meshPoints = mesh.points();
forAll(mesh.boundary(), pI)
{
const polyPatch& patch = mesh.boundaryMesh()[pI];
const labelList& localPoints = patch.meshPoints();
forAll(localPoints, ppI)
{
label pointI = localPoints[ppI];
normalDisplacement[pointI] =
(meshPoints[pointI] - points0[pointI])
& pointNormals[pointI];
}
}
normalDisplacement.write();
// Vectorial displacement
displacement.primitiveFieldRef() = meshPoints - points0;
displacement.write();
}
// Print execution time
mesh.time().printExecutionTime(Info);
Info<< "End\n" << endl;
return(0);
}
// ************************************************************************* //
writeActiveDesignVariables.C
EXE = $(FOAM_APPBIN)/writeActiveDesignVariables
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/finiteVolume/cfdTools \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/optimisation/adjointOptimisation/adjoint/lnInclude
EXE_LIBS = \
-lfiniteVolume \
-lmeshTools \
-lfvOptions \
-ldynamicMesh \
-lfvMotionSolvers \
-ladjointOptimisation
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2007-2019 PCOpt/NTUA
Copyright (C) 2013-2019 FOSS GP
Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
writeActiveDesignVariables
Description
Writes the active design variables based on the selected parameterisation
scheme, as read from the meshMovement part of optimisationDict.
Keeps a back-up of the original optimisationDict in
system/optimisationDict.org, as comments in the dict will be lost.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "optMeshMovement.H"
#include "updateMethod.H"
#include "OFstream.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
IOdictionary optDict
(
IOobject
(
"optimisationDict",
mesh.time().caseSystem(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
// Back-up old optimisationDict, to maintain potential comments in it
if (Pstream::master())
{
cp(optDict.objectPath(), optDict.objectPath() + ".org");
}
// Construct mesh movement object and grab active design variables
// Will exit with error if not implemented for this type
const dictionary& movementDict =
optDict.subDict("optimisation").subDict("meshMovement");
// Empty patch list will do
labelList patchIDs(0);
autoPtr<optMeshMovement> movementPtr
(optMeshMovement::New(mesh, movementDict, patchIDs));
const labelList activeDesignVariables =
movementPtr().getActiveDesignVariables();
// Construct update method to grab the type
dictionary& updateMethodDict =
optDict.subDict("optimisation").subDict("updateMethod");
autoPtr<updateMethod> updMethod(updateMethod::New(mesh, updateMethodDict));
// Add to appropriate dictionary (creates it if it does not exist)
dictionary& coeffsDict = updateMethodDict.subDictOrAdd(updMethod().type());
coeffsDict.add<labelList>
(
"activeDesignVariables",
activeDesignVariables,
true
);
// Write modified dictionary
optDict.regIOobject::writeObject
(
IOstream::ASCII,
IOstream::currentVersion,
IOstream::UNCOMPRESSED,
true
);
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //
writeMorpherCPs.C
EXE = $(FOAM_APPBIN)/writeMorpherCPs
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/finiteVolume/cfdTools \
-I$(LIB_SRC)/optimisation/adjointOptimisation/adjoint/lnInclude
EXE_LIBS = \
-lfiniteVolume \
-lmeshTools \
-lfvOptions \
-ldynamicMesh \
-lfvMotionSolvers \
-ladjointOptimisation
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2007-2019 PCOpt/NTUA
Copyright (C) 2013-2019 FOSS GP
Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
writeMorpherCPs
Description
Writes the NURBS3DVolume control points corresponding to dynamicMeshDict
to csv files. Parametric coordinates are not computed.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "NURBS3DVolume.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
IOdictionary dict
(
IOobject
(
"dynamicMeshDict",
mesh.time().constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
const dictionary& coeffDict =
dict.subDict("volumetricBSplinesMotionSolverCoeffs");
wordList controlBoxes(coeffDict.get<wordList>("controlBoxes"));
forAll(controlBoxes, iNURB)
{
// Creating an object writes the control points in the
// constructor
NURBS3DVolume::New
(
coeffDict.subDict(controlBoxes[iNURB]),
mesh,
false // do not compute parametric coordinates
);
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //
......@@ -52,6 +52,9 @@ pointCells::pointCells
:
zeroATCcells(mesh, dict)
{
boolList isZeroed(mesh_.nCells(), false);
labelList zeroedIDs(mesh_.nCells(), -1);
label i(0);
forAll(mesh_.boundary(), patchI)
{
const fvPatch& patch = mesh_.boundary()[patchI];
......@@ -65,7 +68,14 @@ pointCells::pointCells
for (const label pointI : meshPoints)
{
const labelList& pointCells = mesh_.pointCells()[pointI];
zeroATCcells_.append(pointCells);
for (const label cellI : pointCells)
{
if (!isZeroed[cellI])
{
zeroedIDs[i++] = cellI;
isZeroed[cellI] = true;
}
}
}
}
}
......@@ -76,14 +86,22 @@ pointCells::pointCells
if (zoneID != -1)
{
const labelList& zoneCells = mesh_.cellZones()[zoneID];
zeroATCcells_.append(zoneCells);
for (const label cellI : zoneCells)
{
if (!isZeroed[cellI])
{
zeroedIDs[i++] = cellI;
isZeroed[cellI] = true;
}
}
}
}
/*
zeroedIDs.setSize(i);
zeroATCcells_ = zeroedIDs;
label size = zeroATCcells_.size();
reduce(size, sumOp<label>());
Info<< "Zeroing ATC on "<< size << " cells" << nl << endl;
*/
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -16,6 +16,7 @@ solvers/variablesSet/incompressibleAdjoint/incompressibleAdjointVars.C
solvers/solverControl/solverControl/solverControl.C
solvers/solverControl/SIMPLEControl/SIMPLEControl/SIMPLEControl.C
solvers/solverControl/SIMPLEControl/singleRun/SIMPLEControlSingleRun.C
solvers/solverControl/SIMPLEControl/optimisation/SIMPLEControlOpt.C
/* SOLVERS */
solvers/solver/solver.C
......@@ -52,6 +53,8 @@ objectives/incompressible/objectiveIncompressible/objectiveIncompressible.C
objectives/incompressible/objectiveForce/objectiveForce.C
objectives/incompressible/objectiveMoment/objectiveMoment.C
objectives/incompressible/objectivePtLosses/objectivePtLosses.C
objectives/incompressible/objectiveForceTarget/objectiveForceTarget.C
objectives/incompressible/objectivePartialVolume/objectivePartialVolume.C
/* OBJECTIVE MANAGER*/
objectiveManager/objectiveManager/objectiveManager.C
......@@ -87,8 +90,39 @@ adjointBoundaryConditions/adjointOutletVelocityFlux/adjointOutletVelocityFluxFvP
/* DELTA BOUNDARY */
deltaBoundary/deltaBoundary.C
/* NURBS */
NURBS=parameterization/NURBS
$(NURBS)/NURBSbasis/NURBSbasis.C
$(NURBS)/NURBS3DCurve/NURBS3DCurve.C
$(NURBS)/NURBS3DSurface/NURBS3DSurface.C
$(NURBS)/NURBS3DVolume/NURBS3DVolume/NURBS3DVolume.C
$(NURBS)/NURBS3DVolume/cartesian/NURBS3DVolumeCartesian.C
$(NURBS)/NURBS3DVolume/cylindrical/NURBS3DVolumeCylindrical.C
$(NURBS)/NURBS3DVolume/volBSplinesBase/volBSplinesBase.C
/* BEZIER */
parameterization/Bezier/Bezier.C
/* MOTION SOLVERS */
dynamicMesh/motionSolver/volumetricBSplinesMotionSolver/volumetricBSplinesMotionSolver.C
dynamicMesh/motionSolver/elasticityMotionSolver/elasticityMotionSolver.C
dynamicMesh/motionSolver/laplacianMotionSolver/laplacianMotionSolver.C
/* DISPLACEMENT METHOD */
displacementMethod/displacementMethod/displacementMethod.C
displacementMethod/displacementMethodvolumetricBSplinesMotionSolver/displacementMethodvolumetricBSplinesMotionSolver.C
displacementMethod/displacementMethoddisplacementLaplacian/displacementMethoddisplacementLaplacian.C
displacementMethod/displacementMethodvelocityLaplacian/displacementMethodvelocityLaplacian.C
displacementMethod/displacementMethodelasticityMotionSolver/displacementMethodelasticityMotionSolver.C
displacementMethod/displacementMethodlaplacianMotionSolver/displacementMethodlaplacianMotionSolver.C
/* INTERPOLATION SCHEMES */
interpolation/pointVolInterpolation/pointVolInterpolation.C
/* ADJOINT SENSITIVITY */
optimisation/adjointSensitivity/sensitivity/sensitivity.C
optimisation/adjointSensitivity/shapeSensitivitiesBase/shapeSensitivitiesBase.C
incoSens=optimisation/adjointSensitivity/incompressible