applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Application
  25     ajointShapeOptimizationFoam
  26
  27 Description
  28     Steady-state solver for incompressible, turbulent flow of non-Newtonian
  29     fluids with optimisation of duct shape by applying "blockage" in regions
  30     causing pressure loss as estimated using an adjoint formulation.
  31
  32     References:
  33     \verbatim
  34         "Implementation of a continuous adjoint for topology optimization of
  35          ducted flows"
  36         C. Othmer,
  37         E. de Villiers,
  38         H.G. Weller
  39         AIAA-2007-3947
  40         http://pdf.aiaa.org/preview/CDReadyMCFD07_1379/PV2007_3947.pdf
  41     \endverbatim
  42
  43     Note that this solver optimises for total pressure loss whereas the
  44     above paper describes the method for optimising power-loss.
  45
  46 \*---------------------------------------------------------------------------*/
  47
  48 #include "fvCFD.H"
  49 #include "singlePhaseTransportModel.H"
  50 #include "RASModel.H"
  51 #include "simpleControl.H"
  52
  53 template<class Type>
  54 void zeroCells
  55 (
  56     GeometricField<Type, fvPatchField, volMesh>& vf,
  57     const labelList& cells
  58 )
  59 {
  60     forAll(cells, i)
  61     {
  62         vf[cells[i]] = pTraits<Type>::zero;
  63     }
  64 }
  65
  66
  67 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  68
  69 int main(int argc, char *argv[])
  70 {
  71     #include "setRootCase.H"
  72
  73     #include "createTime.H"
  74     #include "createMesh.H"
  75     #include "createFields.H"
  76     #include "initContinuityErrs.H"
  77     #include "initAdjointContinuityErrs.H"
  78
  79     simpleControl simple(mesh);
  80
  81     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  82
  83     Info<< "\nStarting time loop\n" << endl;
  84
  85     while (simple.loop())
  86     {
  87         Info<< "Time = " << runTime.timeName() << nl << endl;
  88
  89         p.storePrevIter();
  90
  91         laminarTransport.lookup("lambda") >> lambda;
  92
  93         //alpha +=
  94         //    mesh.relaxationFactor("alpha")
  95         //   *(lambda*max(Ua & U, zeroSensitivity) - alpha);
  96         alpha +=
  97             mesh.relaxationFactor("alpha")
  98            *(min(max(alpha + lambda*(Ua & U), zeroAlpha), alphaMax) - alpha);
  99
 100         zeroCells(alpha, inletCells);
 101         //zeroCells(alpha, outletCells);
 102
 103         // Pressure-velocity SIMPLE corrector
 104         {
 105             // Momentum predictor
 106
 107             tmp<fvVectorMatrix> UEqn
 108             (
 109                 fvm::div(phi, U)
 110               + turbulence->divDevReff(U)
 111               + fvm::Sp(alpha, U)
 112             );
 113
 114             UEqn().relax();
 115
 116             solve(UEqn() == -fvc::grad(p));
 117
 118             p.boundaryField().updateCoeffs();
 119             volScalarField rAU(1.0/UEqn().A());
 120             U = rAU*UEqn().H();
 121             UEqn.clear();
 122             phi = fvc::interpolate(U) & mesh.Sf();
 123             adjustPhi(phi, U, p);
 124
 125             // Non-orthogonal pressure corrector loop
 126             for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
 127             {
 128                 fvScalarMatrix pEqn
 129                 (
 130                     fvm::laplacian(rAU, p) == fvc::div(phi)
 131                 );
 132
 133                 pEqn.setReference(pRefCell, pRefValue);
 134                 pEqn.solve();
 135
 136                 if (nonOrth == simple.nNonOrthCorr())
 137                 {
 138                     phi -= pEqn.flux();
 139                 }
 140             }
 141
 142             #include "continuityErrs.H"
 143
 144             // Explicitly relax pressure for momentum corrector
 145             p.relax();
 146
 147             // Momentum corrector
 148             U -= rAU*fvc::grad(p);
 149             U.correctBoundaryConditions();
 150         }
 151
 152         pa.storePrevIter();
 153
 154         // Adjoint Pressure-velocity SIMPLE corrector
 155         {
 156             // Adjoint Momentum predictor
 157
 158             volVectorField adjointTransposeConvection((fvc::grad(Ua) & U));
 159             //volVectorField adjointTransposeConvection
 160             //(
 161             //    fvc::reconstruct
 162             //    (
 163             //        mesh.magSf()*(fvc::snGrad(Ua) & fvc::interpolate(U))
 164             //    )
 165             //);
 166
 167             zeroCells(adjointTransposeConvection, inletCells);
 168
 169             tmp<fvVectorMatrix> UaEqn
 170             (
 171                 fvm::div(-phi, Ua)
 172               - adjointTransposeConvection
 173               + turbulence->divDevReff(Ua)
 174               + fvm::Sp(alpha, Ua)
 175             );
 176
 177             UaEqn().relax();
 178
 179             solve(UaEqn() == -fvc::grad(pa));
 180
 181             pa.boundaryField().updateCoeffs();
 182             volScalarField rAUa(1.0/UaEqn().A());
 183             Ua = rAUa*UaEqn().H();
 184             UaEqn.clear();
 185             phia = fvc::interpolate(Ua) & mesh.Sf();
 186             adjustPhi(phia, Ua, pa);
 187
 188             // Non-orthogonal pressure corrector loop
 189             for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
 190             {
 191                 fvScalarMatrix paEqn
 192                 (
 193                     fvm::laplacian(rAUa, pa) == fvc::div(phia)
 194                 );
 195
 196                 paEqn.setReference(paRefCell, paRefValue);
 197                 paEqn.solve();
 198
 199                 if (nonOrth == simple.nNonOrthCorr())
 200                 {
 201                     phia -= paEqn.flux();
 202                 }
 203             }
 204
 205             #include "adjointContinuityErrs.H"
 206
 207             // Explicitly relax pressure for adjoint momentum corrector
 208             pa.relax();
 209
 210             // Adjoint momentum corrector
 211             Ua -= rAUa*fvc::grad(pa);
 212             Ua.correctBoundaryConditions();
 213         }
 214
 215         turbulence->correct();
 216
 217         runTime.write();
 218
 219         Info<< "ExecutionTime = "
 220             << runTime.elapsedCpuTime()
 221             << " s\n\n" << endl;
 222     }
 223
 224     Info<< "End\n" << endl;
 225
 226     return(0);
 227 }
 228
 229
 230 // ************************************************************************* //