zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: added adjoint-based topology optimisation (topO) capabilities

Browse Source 
Both porosity-based and level-set-based topO frameworks are included
through the topO and levelSet designVariables, respectively.

Both frameworks work by manipulating an underlying field of design
variables, defined in all cells of the computational domain. That field
is then regularised through a Helmholtz-like filter, before being
processed in a different way from the two topO frameworks (the
porosity-based topO sharpens/projects it while the level-set-based topO
computes signed distances around its zero iso-surface). The result of
this processing is then fed into functions that define source terms to
be added to the mean flow and turbulence model equations, to block
off/solidify parts of the mesh that are counterproductive with respect
to the objective function. These source terms are added through
fvOptions.

Since the designed walls are only simulated through source terms, the
outcome of topO should be re-analyzed on a body-fitted grid, to quantify
the actual gain in the objective function. Both topO frameworks output
the designed wall in STL format which can be used, for instance with
snappyHexMesh, to construct such a body fitted grid.
This commit is contained in:
Vaggelis Papoutsis
2023-07-19 11:28:22 +03:00
committed by Andrew Heather
parent 2b2c78309c
commit b23b09fdbd
81 changed files with 12098 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files

38
src/optimisation/adjointOptimisation/adjoint/Make/files
View File

@ -6,6 +6,13 @@ turbulenceModels/turbulenceModelVariables/RAS/kOmegaSST/kOmegaSST.C
turbulenceModels/turbulenceModelVariables/RAS/kEpsilon/kEpsilon.C
turbulenceModels/turbulenceModelVariables/RAS/LaunderSharmaKE/LaunderSharmaKE.C
/* FVOPTIONS */
fvOptions/sources/TopO/topOSource/topOSource.C
/* FUNCTION1 */
OpenFOAM/primitives/functions/Function1/stepRamp/stepRamp.C
OpenFOAM/primitives/functions/Function1/reverseRamp/reverseRamp.C
/* VARIABLES SET */
solvers/variablesSet/variablesSet/variablesSet.C
solvers/variablesSet/incompressible/incompressibleVars.C
@ -57,6 +64,8 @@ objectives/incompressible/objectiveUniformityPatch/objectiveUniformityPatch.C
objectives/incompressible/objectiveUniformityCellZone/objectiveUniformityCellZone.C
objectives/geometric/objectiveGeometric/objectiveGeometric.C
objectives/geometric/objectivePartialVolume/objectivePartialVolume.C
objectives/geometric/objectiveTopOVolume/objectiveTopOVolume.C
objectives/geometric/objectiveTopOSolidVolume/objectiveTopOSolidVolume.C
/* OBJECTIVE MANAGER*/
objectiveManager/objectiveManager.C
@ -152,6 +161,7 @@ optimisation/adjointSensitivity/adjointSensitivity/multiple/sensitivityMultiple.
/* LINE SEARCH */
optimisation/lineSearch/lineSearch/lineSearch.C
optimisation/lineSearch/ArmijoConditions/ArmijoConditions.C
optimisation/lineSearch/GCMMA/GCMMA.C
optimisation/lineSearch/stepUpdate/stepUpdate/stepUpdate.C
optimisation/lineSearch/stepUpdate/bisection/bisection.C
optimisation/lineSearch/stepUpdate/quadratic/quadratic.C
@ -167,6 +177,7 @@ $(updateMethod)/DBFGS/DBFGS.C
$(updateMethod)/LBFGS/LBFGS.C
$(updateMethod)/SR1/SR1.C
$(updateMethod)/conjugateGradient/conjugateGradient.C
$(updateMethod)/MMA/MMA.C
$(updateMethod)/constraintProjection/constraintProjection.C
$(updateMethod)/SQPBase/SQPBase.C
$(updateMethod)/SQP/SQP.C
@ -174,6 +185,33 @@ $(updateMethod)/ISQP/ISQP.C
/* DESIGN VARIABLES */
optimisation/designVariables/designVariables/designVariables.C
topoVars=optimisation/designVariables/topODesignVariables
$(topoVars)/topOVariablesBase/topOVariablesBase.C
$(topoVars)/topODesignVariables.C
$(topoVars)/dynamicTopODesignVariables/dynamicTopODesignVariables.C
$(topoVars)/betaMax/betaMax/betaMax.C
$(topoVars)/betaMax/value/betaMaxValue.C
$(topoVars)/betaMax/Darcy/betaMaxDarcy.C
$(topoVars)/betaMax/stepRamp/betaMaxStepRamp.C
$(topoVars)/topOZones/topOZones.C
$(topoVars)/regularisation/fieldRegularisation.C
$(topoVars)/regularisation/regularisationRadius/regularisationRadius/regularisationRadius.C
$(topoVars)/regularisation/regularisationRadius/isotropic/regularisationRadiusIsotropic.C
$(topoVars)/regularisation/regularisationPDE/regularisationPDE/regularisationPDE.C
$(topoVars)/regularisation/regularisationPDE/Helmoltz/Helmholtz.C
$(topoVars)/marchingCells/marchingCells.C
$(topoVars)/interpolationFunctions/interpolationFunction/topOInterpolationFunction.C
$(topoVars)/interpolationFunctions/BorrvallPetersson/BorrvallPeterssonInterpolation.C
$(topoVars)/interpolationFunctions/invBP/invBP.C
$(topoVars)/interpolationFunctions/linear/linearInterpolation.C
$(topoVars)/interpolationFunctions/SIMP/SIMPInterpolation.C
$(topoVars)/interpolationFunctions/sinh/sinhInterpolation.C
$(topoVars)/interpolationFunctions/tanh/tanhInterpolation.C
$(topoVars)/interpolationFunctions/exp/expInterpolation.C
levelSetVars=optimisation/designVariables/levelSet
$(levelSetVars)/levelSetDesignVariables.C
$(levelSetVars)/interpolationFunctions/sigmoidalHeaviside/sigmoidalHeaviside.C
$(levelSetVars)/interpolationFunctions/smoothHeaviside/smoothHeaviside.C
shapeVars=optimisation/designVariables/shape
$(shapeVars)/shapeDesignVariables/shapeDesignVariables.C
$(shapeVars)/volumetricBSplines/morphingBoxConstraints/morphingBoxConstaint/morphingBoxConstraint.C

60
src/optimisation/adjointOptimisation/adjoint/OpenFOAM/primitives/functions/Function1/reverseRamp/reverseRamp.C Normal file
View File

@ -0,0 +1,60 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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/>.
\*---------------------------------------------------------------------------*/
#include "reverseRamp.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace Function1Types
{
makeScalarFunction1(reverseRamp);
}
}
// * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * //
Foam::Function1Types::reverseRamp::reverseRamp
(
const word& entryName,
const dictionary& dict,
const objectRegistry* obrPtr
)
:
ramp(entryName, dict),
minValue_(dict.getOrDefault<scalar>("minValue", Zero)),
interval_(dict.get<scalar>("interval")),
steps_(duration_/interval_)
{}
// ************************************************************************* //

128
src/optimisation/adjointOptimisation/adjoint/OpenFOAM/primitives/functions/Function1/reverseRamp/reverseRamp.H Normal file
View File

@ -0,0 +1,128 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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/>.
Class
Foam::Function1Types::reverseRamp
Description
Reverse ramp function starting from 1 and decreasing stepRamply to a min.
user-specified value from \c start over the \c duration and remaining
constant thereafter.
See also
Foam::Function1Types::ramp
SourceFiles
reverseRamp.C
\*---------------------------------------------------------------------------*/
#ifndef reverseRamp_H
#define reverseRamp_H
#include "ramp.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace Function1Types
{
/*---------------------------------------------------------------------------*\
Class reverseRamp Declaration
\*---------------------------------------------------------------------------*/
class reverseRamp
:
public ramp
{
protected:
// Protected Data Members
//- The min value
scalar minValue_;
//- Interval for reducing the base value
scalar interval_;
//- Steps to reach the min. value
scalar steps_;
private:
// Private Data Members
//- No copy assignment
void operator=(const reverseRamp&) = delete;
public:
//- Runtime type information
TypeName("reverseRamp");
// Constructors
//- Construct from entry and dictionary
reverseRamp
(
const word& entryName,
const dictionary& dict,
const objectRegistry* obrPtr = nullptr
);
//- Destructor
virtual ~reverseRamp() = default;
// Member Functions
//- Return value for time t
inline virtual scalar value(const scalar t) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Function1Types
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "reverseRampI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

53
src/optimisation/adjointOptimisation/adjoint/OpenFOAM/primitives/functions/Function1/reverseRamp/reverseRampI.H Normal file
View File

@ -0,0 +1,53 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
inline Foam::scalar Foam::Function1Types::reverseRamp::value
(
const scalar t
) const
{
// Reverse linear ramp
//return min(max(scalar(1) - (t - start_)/duration_, scalar(0)), scalar(1));
// Reverse step ramp
return min
(
max
(
-(floor((t - start_)/interval_))/steps_
*(scalar(1) - minValue_) + scalar(1),
minValue_
),
scalar(1)
);
}
// ************************************************************************* //

56
src/optimisation/adjointOptimisation/adjoint/OpenFOAM/primitives/functions/Function1/stepRamp/stepRamp.C Normal file
View File

@ -0,0 +1,56 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017 OpenFOAM Foundation
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "stepRamp.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace Function1Types
{
makeScalarFunction1(stepRamp);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::Function1Types::stepRamp::stepRamp
(
const word& entryName,
const dictionary& dict,
const objectRegistry* obrPtr
)
:
ramp(entryName, dict, obrPtr),
interval_(dict.get<scalar>("interval")),
steps_(max(floor(duration_/interval_), scalar(1)))
{}
// ************************************************************************* //

118
src/optimisation/adjointOptimisation/adjoint/OpenFOAM/primitives/functions/Function1/stepRamp/stepRamp.H Normal file
View File

@ -0,0 +1,118 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017 OpenFOAM Foundation
-------------------------------------------------------------------------------
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/>.
Class
Foam::Function1Types::stepRamp
Description
Linear ramp function starting from 0 and increasing stepRamply to 1 from
\c start over the \c duration and remaining at 1 thereafter.
See also
Foam::Function1Types::ramp
SourceFiles
stepRamp.C
\*---------------------------------------------------------------------------*/
#ifndef stepRamp_H
#define stepRamp_H
#include "ramp.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace Function1Types
{
/*---------------------------------------------------------------------------*\
Class stepRamp Declaration
\*---------------------------------------------------------------------------*/
class stepRamp
:
public ramp
{
// Private data
//- Interval for increasing the base value
scalar interval_;
//- Steps to reach maximum value
scalar steps_;
// Private Member Functions
//- No copy assignment
void operator=(const stepRamp&) = delete;
public:
// Runtime type information
TypeName("stepRamp");
// Constructors
//- Construct from entry name and dictionary
stepRamp
(
const word& entryName,
const dictionary& dict,
const objectRegistry* obrPtr = nullptr
);
//- Destructor
virtual ~stepRamp() = default;
// Member Functions
//- Return value for time t
virtual inline scalar value(const scalar t) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Function1Types
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "stepRampI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

41
src/optimisation/adjointOptimisation/adjoint/OpenFOAM/primitives/functions/Function1/stepRamp/stepRampI.H Normal file
View File

@ -0,0 +1,41 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017 OpenFOAM Foundation
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "stepRamp.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline Foam::scalar Foam::Function1Types::stepRamp::value
(
const scalar t
) const
{
return max(min((floor((t -start_)/interval_) + 1)/steps_, 1), 0);
}
// ************************************************************************* //

229
src/optimisation/adjointOptimisation/adjoint/fvOptions/sources/TopO/topOSource/topOSource.C Normal file
View File

@ -0,0 +1,229 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "topOSource.H"
#include "fvMatrices.H"
#include "fvmSup.H"
#include "topOVariablesBase.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace fv
{
defineTypeNameAndDebug(topOSource, 1);
addToRunTimeSelectionTable
(
option,
topOSource,
dictionary
);
}
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh>>
Foam::fv::topOSource::getSource()
{
auto tinterpolant
(
tmp<DimensionedField<scalar, volMesh>>::New
(
IOobject
(
"source",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimless/dimTime,
scalarField(mesh_.nCells(), Zero)
)
);
DimensionedField<scalar, volMesh>& interpolant = tinterpolant.ref();
if (mesh_.foundObject<topOVariablesBase>("topoVars"))
{
const topOVariablesBase& vars =
mesh_.lookupObject<topOVariablesBase>("topoVars");
vars.sourceTerm
(interpolant, interpolation_(), betaMax_, interpolationFieldName_);
if (darcyFlow_)
{
interpolant.field() += betaMax_*Da_();
}
}
return tinterpolant;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::fv::topOSource::topOSource
(
const word& name,
const word& modelType,
const dictionary& dict,
const fvMesh& mesh
)
:
option(name, modelType, dict, mesh),
interpolation_(topOInterpolationFunction::New(mesh, dict)),
interpolationFieldName_(word::null),
betaMax_(0),
darcyFlow_(false),
Da_(nullptr)
{
read(dict);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::fv::topOSource::addSup
(
fvMatrix<vector>& eqn,
const label fieldi
)
{
DebugInfo
<< "Adding Brinkman source to " << eqn.psi().name() << endl;
eqn -= fvm::Sp(getSource(), eqn.psi());
}
void Foam::fv::topOSource::addSup
(
fvMatrix<scalar>& eqn,
const label fieldi
)
{
DebugInfo
<< "Adding Brinkman source to " << eqn.psi().name() << endl;
eqn -= fvm::Sp(getSource(), eqn.psi());
}
void Foam::fv::topOSource::addSup
(
const volScalarField& rho,
fvMatrix<vector>& eqn,
const label fieldi
)
{
DebugInfo
<< "Adding Brinkman source to " << eqn.psi().name() << endl;
eqn -= fvm::Sp(rho*getSource(), eqn.psi());
}
void Foam::fv::topOSource::addSup
(
const volScalarField& rho,
fvMatrix<scalar>& eqn,
const label fieldi
)
{
DebugInfo
<< "Adding Brinkman source to " << eqn.psi().name() << endl;
eqn -= fvm::Sp(rho*getSource(), eqn.psi());
}
void Foam::fv::topOSource::postProcessSens
(
scalarField& sens,
const word& fieldName,
const word& designVariablesName
)
{
const label fieldi = applyToField(fieldName);
if
(
fieldi != -1
&& mesh_.foundObject<topOVariablesBase>("topoVars")
)
{
DebugInfo
<< "Postprocessing Brinkman sensitivities for field "
<< fieldName << endl;
const topOVariablesBase& vars =
mesh_.lookupObject<topOVariablesBase>("topoVars");
vars.sourceTermSensitivities
(
sens,
interpolation_(),
betaMax_,
designVariablesName,
interpolationFieldName_
);
}
}
bool Foam::fv::topOSource::read(const dictionary& dict)
{
if (option::read(dict))
{
fieldNames_ = coeffs_.get<wordList>("names");
interpolationFieldName_ = coeffs_.get<word>("interpolationField");
applied_.setSize(fieldNames_.size(), false);
if (mesh_.foundObject<topOVariablesBase>("topoVars"))
{
const topOVariablesBase& vars =
mesh_.lookupObject<topOVariablesBase>("topoVars");
betaMax_ =
coeffs_.getOrDefault<scalar>("betaMax", vars.getBetaMax());
}
darcyFlow_ = coeffs_.getOrDefault<bool>("darcyFlow", false);
if (darcyFlow_)
{
Da_.reset(new scalar(coeffs_.getOrDefault<scalar>("Da", 1.e-5)));
}
return true;
}
return false;
}
// ************************************************************************* //

183
src/optimisation/adjointOptimisation/adjoint/fvOptions/sources/TopO/topOSource/topOSource.H Normal file
View File

@ -0,0 +1,183 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::fv::topOSource
Group
grpFvOptionsSources
Description
Impelements Brinkman penalisation terms for topology optimisation.
Looks up the indicator field (beta) from the registry, through
topOVariablesBase
SourceFiles
topOSource.C
\*---------------------------------------------------------------------------*/
#ifndef topOSource_H
#define topOSource_H
#include "cellSetOption.H"
#include "topOInterpolationFunction.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace fv
{
/*---------------------------------------------------------------------------*\
Class topOSource Declaration
\*---------------------------------------------------------------------------*/
class topOSource
:
public option
{
protected:
// Protected data
//- Interpolation function
autoPtr<topOInterpolationFunction> interpolation_;
//- Interpolation field name
word interpolationFieldName_;
//- Optional betaMax
// If not found, the one known by topOVariablesBase will be used.
scalar betaMax_;
//- Does this option apply to a Darcy flow model
bool darcyFlow_;
//- Dimensionless Darcy number
autoPtr<scalar> Da_;
// Protected Member Functions
//- Compute the source term based on the indicator field
virtual tmp<DimensionedField<scalar, volMesh>> getSource();
private:
// Private Member Functions
//- No copy construct
topOSource(const topOSource&) = delete;
//- No copy assignment
void operator=(const topOSource&) = delete;
public:
//- Runtime type information
TypeName("topOSource");
// Constructors
//- Construct from components
topOSource
(
const word& name,
const word& modelType,
const dictionary& dict,
const fvMesh& mesh
);
//- Destructor
virtual ~topOSource() = default;
// Member Functions
//- Add implicit contribution to momentum equation
virtual void addSup
(
fvMatrix<vector>& eqn,
const label fieldi
);
//- Add implicit contribution to scalar equations
//- (e.g. turbulence model)
virtual void addSup
(
fvMatrix<scalar>& eqn,
const label fieldi
);
//- Add implicit contribution to compressible momentum equation
virtual void addSup
(
const volScalarField& rho,
fvMatrix<vector>& eqn,
const label fieldi
);
//- Add implicit contribution to compressible scalar equation
virtual void addSup
(
const volScalarField& rho,
fvMatrix<scalar>& eqn,
const label fieldi
);
//- Multiply sensitivities with the derivative of the interpolation
//- function
virtual void postProcessSens
(
scalarField& sensField,
const word& fieldName = word::null,
const word& designValue = word::null
);
//- Read source dictionary
virtual bool read(const dictionary& dict);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace fv
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

132
src/optimisation/adjointOptimisation/adjoint/objectives/geometric/objectiveTopOSolidVolume/objectiveTopOSolidVolume.C Normal file
View File

@ -0,0 +1,132 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2007-2023 PCOpt/NTUA
Copyright (C) 2013-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "objectiveTopOSolidVolume.H"
#include "createZeroField.H"
#include "IOmanip.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace objectives
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(objectiveTopOSolidVolume, 1);
addToRunTimeSelectionTable
(
objectiveGeometric,
objectiveTopOSolidVolume,
dictionary
);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
objectiveTopOSolidVolume::objectiveTopOSolidVolume
(
const fvMesh& mesh,
const dictionary& dict,
const word& adjointSolverName,
const word& primalSolverName
)
:
objectiveGeometric(mesh, dict, adjointSolverName, primalSolverName),
targetPercentage_(Function1<scalar>::New("percentage", dict)),
percentInDenom_(dict.getOrDefault<bool>("percentInDenom", true))
{
// Allocate boundary field pointers
dJdbPtr_.reset(createZeroFieldPtr<scalar>(mesh_, "dJdb", dimless));
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
scalar objectiveTopOSolidVolume::J()
{
J_ = Zero;
if (mesh_.foundObject<volScalarField>("beta"))
{
const volScalarField& beta = mesh_.lookupObject<volScalarField>("beta");
const DimensionedField<scalar, volMesh>& V = mesh_.V();
const scalar time = mesh_.time().timeOutputValue();
J_ =
gSum(beta.primitiveField()*V)/gSum(V)
- targetPercentage_->value(time);
if (percentInDenom_)
{
J_ /= targetPercentage_->value(time);
}
}
else
{
WarningInFunction
<< "Beta field not yet registered in database. OK for start-up"
<< endl;
}
return J_;
}
void objectiveTopOSolidVolume::update_dJdb()
{
const scalar time = mesh_.time().timeOutputValue();
dJdbPtr_().primitiveFieldRef() = scalar(1)/gSum(mesh_.V());
if (percentInDenom_)
{
dJdbPtr_().primitiveFieldRef() /= targetPercentage_->value(time);
}
}
void objectiveTopOSolidVolume::addHeaderColumns() const
{
objFunctionFilePtr_()
<< setw(width_) << "TargetVolume" << " ";
}
void objectiveTopOSolidVolume::addColumnValues() const
{
const scalar time = mesh_.time().timeOutputValue();
objFunctionFilePtr_()
<< setw(width_) << targetPercentage_->value(time) << " ";
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace objectives
} // End namespace Foam
// ************************************************************************* //

118
src/optimisation/adjointOptimisation/adjoint/objectives/geometric/objectiveTopOSolidVolume/objectiveTopOSolidVolume.H Normal file
View File

@ -0,0 +1,118 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2007-2023 PCOpt/NTUA
Copyright (C) 2013-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::objectives::objectiveTopOSolidVolume
Description
Objective quantifying the difference between the volume occupied by solid
in topology optimisation and a target percentage; the latter can change
throughout the optimisation cycles through a Function1.
SourceFiles
objectiveTopOSolidVolume.C
\*---------------------------------------------------------------------------*/
#ifndef objectiveTopOSolidVolume_H
#define objectiveTopOSolidVolume_H
#include "objectiveGeometric.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace objectives
{
/*---------------------------------------------------------------------------*\
Class objectiveTopOSolidVolume Declaration
\*---------------------------------------------------------------------------*/
class objectiveTopOSolidVolume
:
public objectiveGeometric
{
// Private data
autoPtr<Function1<scalar>> targetPercentage_;
bool percentInDenom_;
public:
//- Runtime type information
TypeName("topOSolidVolume");
// Constructors
//- from components
objectiveTopOSolidVolume
(
const fvMesh& mesh,
const dictionary& dict,
const word& adjointSolverName,
const word& primalSolverName
);
//- Destructor
virtual ~objectiveTopOSolidVolume() = default;
// Member Functions
//- Return the objective function value
scalar J();
//- Contribution to field sensitivities
virtual void update_dJdb();
// Helper write functions
//- Write headers for additional columns
virtual void addHeaderColumns() const;
//- Write information to additional columns
virtual void addColumnValues() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace objectives
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

132
src/optimisation/adjointOptimisation/adjoint/objectives/geometric/objectiveTopOVolume/objectiveTopOVolume.C Normal file
View File

@ -0,0 +1,132 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2007-2023 PCOpt/NTUA
Copyright (C) 2013-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "objectiveTopOVolume.H"
#include "createZeroField.H"
#include "IOmanip.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace objectives
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(objectiveTopOVolume, 1);
addToRunTimeSelectionTable
(
objectiveGeometric,
objectiveTopOVolume,
dictionary
);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
objectiveTopOVolume::objectiveTopOVolume
(
const fvMesh& mesh,
const dictionary& dict,
const word& adjointSolverName,
const word& primalSolverName
)
:
objectiveGeometric(mesh, dict, adjointSolverName, primalSolverName),
targetPercentage_(Function1<scalar>::New("percentage", dict)),
percentInDenom_(dict.getOrDefault<bool>("percentInDenom", true))
{
// Allocate boundary field pointers
dJdbPtr_.reset(createZeroFieldPtr<scalar>(mesh_, "dJdb", dimless));
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
scalar objectiveTopOVolume::J()
{
J_ = Zero;
if (mesh_.foundObject<volScalarField>("beta"))
{
const volScalarField& beta = mesh_.lookupObject<volScalarField>("beta");
const DimensionedField<scalar, volMesh>& V = mesh_.V();
const scalar time = mesh_.time().timeOutputValue();
J_ =
scalar(1) - gSum(beta.primitiveField()*V)/gSum(V)
- targetPercentage_->value(time);
if (percentInDenom_)
{
J_ /= targetPercentage_->value(time);
}
}
else
{
WarningInFunction
<< "Beta field not yet registered in database. OK for start-up"
<< endl;
}
return J_;
}
void objectiveTopOVolume::update_dJdb()
{
const scalar time = mesh_.time().timeOutputValue();
dJdbPtr_().primitiveFieldRef() = -scalar(1)/gSum(mesh_.V());
if (percentInDenom_)
{
dJdbPtr_().primitiveFieldRef() /= targetPercentage_->value(time);
}
}
void objectiveTopOVolume::addHeaderColumns() const
{
objFunctionFilePtr_()
<< setw(width_) << "TargetVolume" << " ";
}
void objectiveTopOVolume::addColumnValues() const
{
const scalar time = mesh_.time().timeOutputValue();
objFunctionFilePtr_()
<< setw(width_) << targetPercentage_->value(time) << " ";
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace objectives
} // End namespace Foam
// ************************************************************************* //

119
src/optimisation/adjointOptimisation/adjoint/objectives/geometric/objectiveTopOVolume/objectiveTopOVolume.H Normal file
View File

@ -0,0 +1,119 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2007-2023 PCOpt/NTUA
Copyright (C) 2013-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::objectives::objectiveTopOVolume
Description
Objective quantifying the difference between the volume occupied by fluid
in topology optimisation and a target percentage; the latter can change
throughout the optimisation cycles through a Function1.
SourceFiles
objectiveTopOVolume.C
\*---------------------------------------------------------------------------*/
#ifndef objectiveTopOVolume_H
#define objectiveTopOVolume_H
#include "objectiveGeometric.H"
#include "Function1.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace objectives
{
/*---------------------------------------------------------------------------*\
Class objectiveTopOVolume Declaration
\*---------------------------------------------------------------------------*/
class objectiveTopOVolume
:
public objectiveGeometric
{
// Private data
autoPtr<Function1<scalar>> targetPercentage_;
bool percentInDenom_;
public:
//- Runtime type information
TypeName("topOVolume");
// Constructors
//- from components
objectiveTopOVolume
(
const fvMesh& mesh,
const dictionary& dict,
const word& adjointSolverName,
const word& primalSolverName
);
//- Destructor
virtual ~objectiveTopOVolume() = default;
// Member Functions
//- Return the objective function value
scalar J();
//- Contribution to field sensitivities
virtual void update_dJdb();
// Helper write functions
//- Write headers for additional columns
virtual void addHeaderColumns() const;
//- Write information to additional columns
virtual void addColumnValues() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace objectives
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

233
src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/adjointSensitivity/topO/sensitivityTopO.C Normal file
View File

@ -0,0 +1,233 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2007-2023 PCOpt/NTUA
Copyright (C) 2013-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "adjointSensitivity.H"
#include "sensitivityTopO.H"
#include "adjointSolver.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
defineTypeNameAndDebug(sensitivityTopO, 0);
addToRunTimeSelectionTable(adjointSensitivity, sensitivityTopO, dictionary);
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
void sensitivityTopO::zeroSensInFixedPorousZones(scalarField& sens)
{
const labelList& adjointPorousIDs = zones_.adjointPorousZoneIDs();
if (adjointPorousIDs.empty())
{
for (label cellZoneID : zones_.fixedPorousZoneIDs())
{
const labelList& zoneCells = mesh_.cellZones()[cellZoneID];
for (label cellI : zoneCells)
{
sens[cellI] = 0.;
}
}
for (label cellZoneID : zones_.fixedZeroPorousZoneIDs())
{
const labelList& zoneCells = mesh_.cellZones()[cellZoneID];
for (label cellI : zoneCells)
{
sens[cellI] = 0.;
}
}
for (label cellI : zones_.IOCells())
{
sens[cellI] = 0.;
}
}
else
{
// if adjointPorousZones are not empty, zero sensitivities in all cells
// not within them
scalarField mask(sens.size(), Zero);
for (label cellZoneID : adjointPorousIDs)
{
const labelList& zoneCells = mesh_.cellZones()[cellZoneID];
for (label cellI : zoneCells)
{
mask[cellI] = 1.;
}
}
sens *= mask;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
sensitivityTopO::sensitivityTopO
(
const fvMesh& mesh,
const dictionary& dict,
adjointSolver& adjointSolver
)
:
adjointSensitivity(mesh, dict, adjointSolver),
zones_(mesh, dict.parent()),
designVariablesName_("beta")
{
if (includeDistance_)
{
eikonalSolver_.reset
(
new adjointEikonalSolver
(
mesh_,
dict_,
adjointSolver,
labelHashSet(0)
)
);
}
// Allocate sensitivities field
fieldSensPtr_.reset
(
createZeroFieldPtr<scalar>
(
mesh_,
"topologySens" + adjointSolver.solverName(),
pow5(dimLength)/sqr(dimTime)
)
);
// Set return field size
derivatives_ = scalarField(mesh_.nCells(), Zero);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool sensitivityTopO::readDict(const dictionary& dict)
{
if (adjointSensitivity::readDict(dict))
{
if (includeDistance_)
{
if (eikonalSolver_)
{
eikonalSolver_->readDict(dict);
}
else
{
eikonalSolver_.reset
(
new adjointEikonalSolver
(
mesh_,
dict_,
adjointSolver_,
labelHashSet(0)
)
);
}
}
return true;
}
return false;
}
void sensitivityTopO::accumulateIntegrand(const scalar dt)
{
// Accumulate source for additional post-processing PDEs, if necessary
if (eikonalSolver_)
{
eikonalSolver_->accumulateIntegrand(dt);
}
adjointSolver_.topOSensMultiplier
(fieldSensPtr_().primitiveFieldRef(), designVariablesName_, dt);
}
void sensitivityTopO::assembleSensitivities
(
autoPtr<designVariables>& designVars
)
{
scalarField& sens = fieldSensPtr_().primitiveFieldRef();
if (eikonalSolver_)
{
eikonalSolver_->solve();
sens += eikonalSolver_->topologySensitivities(designVariablesName_);
}
zeroSensInFixedPorousZones(sens);
adjointSensitivity::assembleSensitivities(designVars);
}
void sensitivityTopO::postProcessSens
(
scalarField& sens,
scalarField& auxSens,
fv::options& fvOptions,
const word& fieldName,
const word& designVariablesName
)
{
if (fvOptions.appliesToField(fieldName))
{
DebugInfo
<< "Computing SD contributions from the interpolation of "
<< fieldName << endl;
fvOptions.postProcessSens(auxSens, fieldName, designVariablesName);
sens += auxSens;
}
}
void sensitivityTopO::postProcessSens
(
scalarField& sens,
scalarField& auxSens,
const word& fieldName
)
{
fv::options& fvOptions(fv::options::New(this->mesh_));
postProcessSens(sens, auxSens, fvOptions, fieldName, designVariablesName_);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

173
src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/adjointSensitivity/topO/sensitivityTopO.H Normal file
View File

@ -0,0 +1,173 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2007-2020 PCOpt/NTUA
Copyright (C) 2013-2020 FOSS GP
-------------------------------------------------------------------------------
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/>.
Description
Calculation of adjoint based sensitivities for topology optimisation.
This returns just the field part, with contributions from regularisation
and projection added by topODesignVariables.
Class
Foam::sensitivityTopO
SourceFiles
sensitivityTopO.C
\*---------------------------------------------------------------------------*/
#ifndef sensitivityTopO_H
#define sensitivityTopO_H
#include "adjointSensitivity.H"
#include "fvOptions.H"
#include "topOZones.H"
#include "adjointEikonalSolver.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class sensitivityTopO Declaration
\*---------------------------------------------------------------------------*/
class sensitivityTopO
:
public adjointSensitivity
{
protected:
// Protected data
//- Zones related to topology optimisation
topOZones zones_;
//- Name used as the argument for the post-processing of the
//- sensitivities through fvOptions
word designVariablesName_;
// Protected Member Functions
//- Zero sensitivities in fixed porous zones
void zeroSensInFixedPorousZones(scalarField& sens);
private:
// Private Member Functions
//- No copy construct
sensitivityTopO(const sensitivityTopO&);
//- No copy assignment
void operator=(const sensitivityTopO&);
public:
//- Runtime type information
TypeName("topO");
// Constructors
//- Construct from components
sensitivityTopO
(
const fvMesh& mesh,
const dictionary& dict,
adjointSolver& adjointSolver
);
//- Destructor
virtual ~sensitivityTopO() = default;
// Member Functions
//- Read dictionary if changed
virtual bool readDict(const dictionary& dict);
//- Accumulate sensitivity integrands
virtual void accumulateIntegrand(const scalar dt);
//- Assemble sensitivities
virtual void assembleSensitivities
(
autoPtr<designVariables>& designVars
);
//- Add part of the sensitivities coming from fvOptions
static void postProcessSens
(
scalarField& sens,
scalarField& auxSens,
fv::options& fvOptions,
const word& fieldName,
const word& designVariablesName
);
//- Add part of the sensitivities coming from fvOptions
void postProcessSens
(
scalarField& sens,
scalarField& auxSens,
const word& fieldName
);
//- Potentially manipulate auxSens within the fvOption before adding
//- the part related to the design variables
template<class Type, template<class> class PatchField, class GeoMesh>
static void postProcessAuxSens
(
const GeometricField<Type, PatchField, GeoMesh>& primal,
const GeometricField<Type, PatchField, GeoMesh>& adjoint,
fv::options& fvOptions,
const word& fieldName,
scalarField& auxSens
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "sensitivityTopOTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

51
src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/adjointSensitivity/topO/sensitivityTopOTemplates.C Normal file
View File

@ -0,0 +1,51 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2022 PCOpt/NTUA
Copyright (C) 2022 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type, template<class> class PatchField, class GeoMesh>
void Foam::sensitivityTopO::postProcessAuxSens
(
const GeometricField<Type, PatchField, GeoMesh>& primal,
const GeometricField<Type, PatchField, GeoMesh>& adjoint,
fv::options& fvOptions,
const word& fieldName,
scalarField& auxSens
)
{
if (fvOptions.appliesToField(fieldName))
{
DebugInfo
<< "Manupulation auxiliary SD contributions from field "
<< fieldName << endl;
fvOptions.postProcessAuxSens(primal, adjoint, auxSens, fieldName);
}
}
// ************************************************************************* //

134
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/levelSet/interpolationFunctions/sigmoidalHeaviside/sigmoidalHeaviside.C Normal file
View File

@ -0,0 +1,134 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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/>.
\*---------------------------------------------------------------------------*/
#include "sigmoidalHeaviside.H"
#include "mathematicalConstants.H"
#include "volFields.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * * //
defineTypeNameAndDebug(sigmoidalHeaviside, 1);
addToRunTimeSelectionTable
(
topOInterpolationFunction,
sigmoidalHeaviside,
dictionary
);
// * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
scalar sigmoidalHeaviside::computeNearBandWidth() const
{
scalar averageVol(gAverage(mesh_.V().field()));
const Vector<label>& solutionD = mesh_.solutionD();
const boundBox& bounds = mesh_.bounds();
forAll(solutionD, idir)
{
if (solutionD[idir] == -1)
{
averageVol /= bounds.span()[idir];
break;
}
}
scalar width = pow(averageVol, scalar(1)/scalar(mesh_.nGeometricD()));
scalar multMeanRadius =
dict_.getOrDefaultCompat<scalar>
(
"meanRadiusMult", {{"scale", 2306}}, 2
);
DebugInfo
<< "Computed near-band width :: " << width
<< " and multiplying with " << multMeanRadius << endl;
return multMeanRadius*width;
}
// * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * //
sigmoidalHeaviside::sigmoidalHeaviside
(
const fvMesh& mesh,
const dictionary& dict
)
:
topOInterpolationFunction(mesh, dict),
dNB_(dict.getOrDefault<scalar>("d", computeNearBandWidth()))
//dNB_(Function1<scalar>::New("d", dict))
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void sigmoidalHeaviside::interpolate
(
const scalarField& arg,
scalarField& res
) const
{
const scalar timeValue = mesh_.time().timeOutputValue();
const scalar t(timeValue == 0 ? 1. : timeValue);
const scalar pi = constant::mathematical::pi;
DebugInfo
<< type() << "::interpolate:: t, dNB " << t << ", " << dNB_ << endl;
res = 0.5*(scalar(1) + arg/dNB_ + sin(pi*arg/dNB_)/pi);
res = max(min(scalar(1), res), scalar(0));
}
tmp<scalarField> sigmoidalHeaviside::derivative(const scalarField& arg) const
{
tmp<scalarField> tderiv = tmp<scalarField>::New(arg.size(), Zero);
scalarField& deriv = tderiv.ref();
const scalar timeValue = mesh_.time().timeOutputValue();
const scalar t(timeValue == 0 ? 1. : timeValue);
const scalar pi = constant::mathematical::pi;
scalarField argLimited(max(min(dNB_, arg), -dNB_));
DebugInfo
<< type() << "::interpolate:: t, dNB " << t << ", " << dNB_ << endl;
deriv = 0.5*(scalar(1) + cos(pi*argLimited/dNB_))/dNB_;
return tderiv;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

123
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/levelSet/interpolationFunctions/sigmoidalHeaviside/sigmoidalHeaviside.H Normal file
View File

@ -0,0 +1,123 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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/>.
Class
Foam::sigmoidalHeaviside
Description
A smooth Heaviside function to project the signed distance field in
level set topology optimization.
\*---------------------------------------------------------------------------*/
#ifndef sigmoidalHeaviside_H
#define sigmoidalHeaviside_H
#include "topOInterpolationFunction.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class sigmoidalHeaviside Declaration
\*---------------------------------------------------------------------------*/
class sigmoidalHeaviside
:
public topOInterpolationFunction
{
protected:
// Protected Data Members
//- The near-band distance
scalar dNB_;
// Protected Member Functions
//- Compute the near-band width of the fluid-solid interface as
// the mean edge length of mesh cells
scalar computeNearBandWidth() const;
private:
// Private Member Functions
//- No copy construct
sigmoidalHeaviside(const sigmoidalHeaviside&) = delete;
//- No copy assignment
void operator=(const sigmoidalHeaviside&) = delete;
public:
//- Runtime type information
TypeName("sigmoidalHeaviside");
// Constructors
//- Construct from mesh and dictionary
sigmoidalHeaviside
(
const fvMesh& mesh,
const dictionary& dict
);
// Destructor
virtual ~sigmoidalHeaviside() = default;
// Member Functions
//- Interpolate argument to result
virtual void interpolate
(
const scalarField& arg,
scalarField& res
) const;
//- Return of function with respect to the argument field
virtual tmp<scalarField> derivative(const scalarField& arg) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

93
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/levelSet/interpolationFunctions/smoothHeaviside/smoothHeaviside.C Normal file
View File

@ -0,0 +1,93 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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/>.
\*---------------------------------------------------------------------------*/
#include "smoothHeaviside.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * * //
defineTypeNameAndDebug(smoothHeaviside, 0);
addToRunTimeSelectionTable
(
topOInterpolationFunction,
smoothHeaviside,
dictionary
);
// * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * //
smoothHeaviside::smoothHeaviside
(
const fvMesh& mesh,
const dictionary& dict
)
:
topOInterpolationFunction(mesh, dict),
b_(Function1<scalar>::New("b", dict))
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void smoothHeaviside::interpolate
(
const scalarField& arg,
scalarField& res
) const
{
const scalar timeValue = mesh_.time().timeOutputValue();
const scalar t(timeValue == 0 ? 1. : timeValue);
const scalar b(b_->value(t));
res = 0.5*(scalar(1) + tanh(b*arg));
}
tmp<scalarField> smoothHeaviside::derivative(const scalarField& arg) const
{
tmp<scalarField> tderiv = tmp<scalarField>::New(arg.size(), Zero);
scalarField& deriv = tderiv.ref();
const scalar timeValue = mesh_.time().timeOutputValue();
const scalar t(timeValue == 0 ? 1. : timeValue);
const scalar b(b_->value(t));
deriv = 0.5*b*(scalar(1) - sqr(tanh(b*arg)));
return tderiv;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

115
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/levelSet/interpolationFunctions/smoothHeaviside/smoothHeaviside.H Normal file
View File

@ -0,0 +1,115 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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/>.
Class
Foam::smoothHeaviside
Description
A smooth Heaviside function to project the signed distance field in
level set topology optimization.
\*---------------------------------------------------------------------------*/
#ifndef smoothHeaviside_H
#define smoothHeaviside_H
#include "topOInterpolationFunction.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class smoothHeaviside Declaration
\*---------------------------------------------------------------------------*/
class smoothHeaviside
:
public topOInterpolationFunction
{
protected:
// Protected Data Members
autoPtr<Function1<scalar>> b_;
private:
// Private Member Functions
//- No copy construct
smoothHeaviside(const smoothHeaviside&) = delete;
//- No copy assignment
void operator=(const smoothHeaviside&) = delete;
public:
//- Runtime type information
TypeName("smoothHeaviside");
// Constructors
//- Construct from mesh and dictionary
smoothHeaviside
(
const fvMesh& mesh,
const dictionary& dict
);
// Destructor
virtual ~smoothHeaviside() = default;
// Member Functions
//- Interpolate argument to result
virtual void interpolate
(
const scalarField& arg,
scalarField& res
) const;
//- Return of function with respect to the argument field
virtual tmp<scalarField> derivative(const scalarField& arg) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

581
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/levelSet/levelSetDesignVariables.C Normal file
View File

@ -0,0 +1,581 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2022-2023 PCOpt/NTUA
Copyright (C) 2022-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "HashSet.H"
#include "levelSetDesignVariables.H"
#include "wallDist.H"
#include "zeroGradientFvPatchField.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * * //
defineTypeNameAndDebug(levelSetDesignVariables, 1);
addToRunTimeSelectionTable
(
designVariables,
levelSetDesignVariables,
designVariables
);
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
void levelSetDesignVariables::readField()
{
scalarField& vars = *this;
if (localIOdictionary::found("alpha"))
{
vars = (scalarField("alpha", *this, vars.size()));
}
else
{
// Initialise as the distance from the wall patches of the initial
// domain
const labelHashSet wallPatchIDs =
mesh_.boundaryMesh().findPatchIDs<wallPolyPatch>();
volScalarField y
(
IOobject
(
"yLevelSet",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar(dimLength, Zero),
patchDistMethod::patchTypes<scalar>(mesh_, wallPatchIDs)
);
patchDistMethod::New
(
dict_.subDict("initialisation"),
mesh_,
wallPatchIDs
)->correct(y);
vars = y.primitiveField();
if (debug)
{
writeDesignVars();
}
}
}
void levelSetDesignVariables::applyFixedPorosityValues()
{
scalarField& betaIf = beta_.primitiveFieldRef();
// Safety - set beta equal to zero next to the IO cells
for (const label IOcell : zones_.IOCells())
{
betaIf[IOcell] = Zero;
}
const labelList& fixedZeroPorousZones = zones_.fixedZeroPorousZoneIDs();
const labelList& fixedPorousZones = zones_.fixedPorousZoneIDs();
const scalarList& fixedPorousValues = zones_.fixedPorousValues();
// Apply fixed porosity
for (const label zoneID : fixedZeroPorousZones)
{
const labelList& zone = mesh_.cellZones()[zoneID];
for (const label cellI : zone)
{
betaIf[cellI] = Zero;
}
}
// Apply fixed porosity
forAll(fixedPorousZones, zI)
{
const label zoneID = fixedPorousZones[zI];
const scalar value = fixedPorousValues[zI];
const labelList& zone = mesh_.cellZones()[zoneID];
for (const label cellI : zone)
{
betaIf[cellI] = value >= 0 ? 0 : 1;
}
}
beta_.correctBoundaryConditions();
}
void levelSetDesignVariables::setActiveDesignVariables(bool activeIO)
{
activeDesignVariables_.setSize(mesh_.nCells(), -1);
if (zones_.adjointPorousZoneIDs().empty())
{
boolList isActiveVar(mesh_.nCells(), true);
const labelList& fixedZeroPorousZones =
zones_.fixedZeroPorousZoneIDs();
for (const label zoneID : fixedZeroPorousZones)
{
const labelList& zone = mesh_.cellZones()[zoneID];
for (const label cellI : zone)
{
isActiveVar[cellI] = false;
}
}
const labelList& fixedPorousZones = zones_.fixedPorousZoneIDs();
for (const label zoneID : fixedPorousZones)
{
const labelList& zone = mesh_.cellZones()[zoneID];
for (const label cellI : zone)
{
isActiveVar[cellI] = false;
}
}
if (!activeIO)
{
for (label cellI : zones_.IOCells())
{
isActiveVar[cellI] = false;
}
}
label iVar(0);
forAll(isActiveVar, cI)
{
if (isActiveVar[cI])
{
activeDesignVariables_[iVar++] = cI;
}
}
activeDesignVariables_.setSize(iVar);
}
else
{
const labelList& adjointPorousZoneIDs = zones_.adjointPorousZoneIDs();
label iVar(0);
for (const label cellZoneID : adjointPorousZoneIDs)
{
const labelList& zone = mesh_.cellZones()[cellZoneID];
for (const label cellI : zone)
{
activeDesignVariables_[iVar] = cellI;
iVar++;
}
}
activeDesignVariables_.setSize(iVar);
}
}
void levelSetDesignVariables::updateBeta()
{
// Compute the beta field by passing the distance field through
// a Heaviside function
scalarField& beta = beta_.primitiveFieldRef();
interpolation_->interpolate(aTilda_.primitiveField(), beta);
beta = 1 - beta;
// Apply fixed values if necessary
applyFixedPorosityValues();
beta_.correctBoundaryConditions();
}
void Foam::levelSetDesignVariables::updateSignedDistances()
{
Info<< "Re-initilising the level-set distance field" << nl << endl;
volScalarField y
(
IOobject
(
"yLevelSet",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar(dimLength, Zero),
wordList
(
mesh_.boundary().size(),
zeroGradientFvPatchField<scalar>::typeName
)
);
y.primitiveFieldRef() = aTilda_.primitiveFieldRef();
y.correctBoundaryConditions();
labelList changedFaces(mesh_.nFaces(), -1);
List<wallPoint> changedFacesInfo(mesh_.nFaces());
writeFluidSolidInterface(aTilda_, 0, changedFaces, changedFacesInfo);
List<wallPoint> allFaceInfo(mesh_.nFaces());
List<wallPoint> allCellInfo(mesh_.nCells());
FaceCellWave<wallPoint> wave
(
mesh_,
changedFaces,
changedFacesInfo,
allFaceInfo,
allCellInfo,
mesh_.globalData().nTotalCells() + 1
);
// Transfer the distance from cellInfo to the alphaTilda field
forAll(allCellInfo, celli)
{
if (allCellInfo[celli].valid(wave.data()))
{
aTilda_[celli] =
sign(aTilda_[celli])*Foam::sqrt(allCellInfo[celli].distSqr());
}
}
aTilda_.correctBoundaryConditions();
}
// * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * //
levelSetDesignVariables::levelSetDesignVariables
(
fvMesh& mesh,
const dictionary& dict
)
:
topOVariablesBase(mesh, dict),
designVariables(mesh, dict, mesh.nCells()),
radius_
(regularisationRadius::New(mesh, dict.subDict("regularisation"), false)),
regularisation_
(regularisationPDE::New(mesh, dict.subDict("regularisation"), zones_)),
aTilda_
(
IOobject
(
"signedDistances",
mesh_.time().timeName(),
mesh_,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh_,
dimensionedScalar(dimless, Zero),
zeroGradientFvPatchField<scalar>::typeName
),
interpolation_
(
topOInterpolationFunction::New(mesh_, dict_.subDict("interpolation"))
),
beta_
(
IOobject
(
"beta",
mesh_.time().timeName(),
mesh_,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh_,
dimensionedScalar(dimless, Zero)
),
fixATildaValues_(dict.getOrDefault<bool>("fixATildaValues", true)),
writeAllDistanceFields_
(dict.getOrDefault<bool>("writeAllDistanceFields", false))
{
// Read the alpha field if present, or set it based on the distance field
readField();
// Read bounds of design variables if present.
// If not, use the maximum distance field in the fluid to set them.
scalar maxDist = gMax(*this);
scalar lowerBound =
localIOdictionary::getOrDefault<scalar>("lowerBound", -maxDist - SMALL);
scalar upperBound =
localIOdictionary::getOrDefault<scalar>("upperBound", maxDist + SMALL);
readBounds
(
autoPtr<scalar>(new scalar(lowerBound)),
autoPtr<scalar>(new scalar(upperBound))
);
DebugInfo
<< "Using lower/upper bounds "
<< lowerBounds_()[0] << "/" << upperBounds_()[0]
<< endl;
// Update the beta field based on the initial design vars
scalarField zeroUpdate(scalarField::size(), Zero);
update(zeroUpdate);
// Determine which design variables are active
setActiveDesignVariables();
}
// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
autoPtr<levelSetDesignVariables> levelSetDesignVariables::New
(
fvMesh& mesh,
const dictionary& dict
)
{
return autoPtr<levelSetDesignVariables>
(
new levelSetDesignVariables(mesh, dict)
);
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
const volScalarField& levelSetDesignVariables::beta() const
{
return beta_;
}
void levelSetDesignVariables::update(scalarField& correction)
{
scalarField::operator+=(correction);
// Compute the regularised design variables
regularisation_->regularise
(
aTilda_, *this, aTilda_.primitiveFieldRef(),
true, radius_(), upperBounds_()[0], fixATildaValues_
);
aTilda_.correctBoundaryConditions();
if (writeAllDistanceFields_)
{
writeDesignVars();
aTilda_.rename("alphaSmoothed");
aTilda_.write();
aTilda_.rename("signedDistances");
}
// Make aTilda a signed distance field
updateSignedDistances();
// Set beta based on aTilda
updateBeta();
if (writeAllDistanceFields_)
{
aTilda_.write();
beta_.write();
}
// Though the mesh is kept constant, the distance from wall may change
// due to fvOptions depending on beta. Trick wallDist into updating it
if (mesh_.foundObject<UpdateableMeshObject<fvMesh>>("wallDist"))
{
mesh_.lookupObjectRef<UpdateableMeshObject<fvMesh>>("wallDist").
movePoints();
}
}
scalar levelSetDesignVariables::computeEta(scalarField& correction)
{
// Back-up the old design variables and signed distances
scalarField& dvs = getVars();
scalarField oldDVs(getVars());
scalarField oldSignedDistances(aTilda_.primitiveField());
// Compute the smooth alpha field corresponding to the initial variables
// Can't use current aTilda_ values since they correspond to signed
// distances at this point
scalarField oldATilda(aTilda_.primitiveField());
regularisation_->regularise
(
aTilda_, dvs, oldATilda,
true, radius_(), upperBounds_()[0], fixATildaValues_
);
// Compute the smooth alpha field corresponding to the updated variables
dvs += correction;
regularisation_->regularise
(
aTilda_, dvs, aTilda_.primitiveFieldRef(),
true, radius_(), upperBounds_()[0], fixATildaValues_
);
aTilda_.correctBoundaryConditions();
// We want to locate the min value of aTilda_ and scale the correction
// appropriately such that this min value takes on the prescribed one.
// A bit tricky in parallel since we need not only the min value but
// its cellId/processor too
const label proci = Pstream::myProcNo();
scalarList minVs(Pstream::nProcs(), pTraits<scalar>::max);
labelList minCells(Pstream::nProcs(), Zero);
scalarField diff(aTilda_.primitiveField() - oldATilda);
label minId = findMin(diff);
if (minId != -1)
{
minVs[proci] = diff[minId];
minCells[proci] = minId;
}
// Collect info from all processors
Pstream::allGatherList(minVs);
Pstream::allGatherList(minCells);
minId = findMin(minVs);
scalar aTildaAtMinChange(Zero);
if (proci == minId)
{
const label cellId = minCells[minId];
aTildaAtMinChange = aTilda_.primitiveField()[cellId];
}
reduce(aTildaAtMinChange, sumOp<scalar>());
DebugInfo
<< "AlphaSmoothed at min(alphaSmoothedUpdate) with eta 1/"
<< "min desirable value "
<< minVs[minId] << '/' << maxInitChange_()
<< endl;
// Compute eta
const scalar eta((maxInitChange_() - aTildaAtMinChange)/minVs[minId] + 1);
Info<< "Setting eta value to " << eta << endl;
correction *= eta;
// Restore the dvs
dvs = oldDVs;
aTilda_.primitiveFieldRef() = oldSignedDistances;
return eta;
}
bool levelSetDesignVariables::globalSum() const
{
return true;
}
tmp<scalarField> levelSetDesignVariables::assembleSensitivities
(
adjointSensitivity& adjointSens
)
{
// Raw sensitivities field
const scalarField& fieldSens = adjointSens.fieldSensPtr()->primitiveField();
// Return field
auto tobjectiveSens(tmp<scalarField>::New(fieldSens));
scalarField& objectiveSens = tobjectiveSens.ref();
// Multiply with dBetadAtilda
objectiveSens *= -interpolation_->derivative(aTilda_.primitiveField());
// Solve the adjoint to the regularisation equation
regularisation_->
regularise(aTilda_, objectiveSens, objectiveSens, false, radius_());
objectiveSens *= mesh_.V();
if (writeAllDistanceFields_)
{
volScalarField sens
(
IOobject
(
"sens" + adjointSens.getAdjointSolver().solverName(),
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh_,
dimensionedScalar(dimless, Zero)
);
sens.primitiveFieldRef() = objectiveSens;
sens.write();
}
return tobjectiveSens;
}
void levelSetDesignVariables::writeDesignVars()
{
if (writeAllDistanceFields_ || mesh_.time().writeTime())
{
volScalarField alpha
(
IOobject
(
"alpha",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh_,
dimensionedScalar(dimLength, Zero)
);
alpha.primitiveFieldRef() = *this;
alpha.correctBoundaryConditions();
alpha.write();
}
}
bool levelSetDesignVariables::writeData(Ostream& os) const
{
// Lower and upper bound values might be computed based on the initial
// distance field. Write to file to enable continuation
os.writeEntry("lowerBound", lowerBounds_()[0]);
os.writeEntry("upperBound", upperBounds_()[0]);
scalarField::writeEntry("alpha", os);
return true;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

186
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/levelSet/levelSetDesignVariables.H Normal file
View File

@ -0,0 +1,186 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2022-2023 PCOpt/NTUA
Copyright (C) 2022-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::levelSetDesignVariables
Description
Signed distance field design variables for level-set based topology
optimization.
\*---------------------------------------------------------------------------*/
#ifndef levelSetDesignVariables_H
#define levelSetDesignVariables_H
#include "topOVariablesBase.H"
#include "designVariables.H"
#include "regularisationRadius.H"
#include "regularisationPDE.H"
#include "topOInterpolationFunction.H"
#include "adjointSensitivity.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class levelSetDesignVariables Declaration
\*---------------------------------------------------------------------------*/
class levelSetDesignVariables
:
public topOVariablesBase,
public designVariables
{
protected:
// Protected Data Members
//- The regularisation radius
autoPtr<regularisationRadius> radius_;
//- Regularisation mechanism
autoPtr<regularisationPDE> regularisation_;
//- The regularised field that is also transformed
//- into signed distances
volScalarField aTilda_;
//- Function to transorm signed distances to the indicator field beta_
autoPtr<topOInterpolationFunction> interpolation_;
//- The indicator field
volScalarField beta_;
//- Fix aTilda values in fixed{Zero}PorousZones and IOcells
bool fixATildaValues_;
//- Write all fields related to the distance calculation (debuging)
bool writeAllDistanceFields_;
// Protected Member Functions
//- Read the design variables field
void readField();
//- Apply fixed values in certain zones
void applyFixedPorosityValues();
//- Determine which design variables are active
void setActiveDesignVariables(bool activeIO = false);
//- Make aTilda a signed distance field based on the zero iso-surface
//- of the current aTilda field
void updateSignedDistances();
//- Update the beta field
void updateBeta();
private:
// Private Member Functions
//- No copy construct
levelSetDesignVariables(const levelSetDesignVariables&) = delete;
//- No copy assignment
void operator=(const levelSetDesignVariables&) = delete;
public:
//- Runtime type information
TypeName("levelSet");
// Constructors
//- Construct from dictionary
levelSetDesignVariables
(
fvMesh& mesh,
const dictionary& dict
);
// Selectors
//- Return a reference to the selected turbulence model
static autoPtr<levelSetDesignVariables> New
(
fvMesh& mesh,
const dictionary& dict
);
//- Destructor
virtual ~levelSetDesignVariables() = default;
// Member Functions
//- Const reference to the beta field
virtual const volScalarField& beta() const;
//- Update design variables based on a given correction
virtual void update(scalarField& correction);
//- Compute eta if not set in the first step
virtual scalar computeEta(scalarField& correction);
//- Whether to use global sum when computing matrix-vector products
//- in update methods
virtual bool globalSum() const;
//- Add contributions from the adjoint to the regularization PDE,
//- the derivative of the interpolation function and the mesh volume
//- to the final sensitivities
virtual tmp<scalarField> assembleSensitivities
(
adjointSensitivity& adjointSens
);
//- Write useful quantities to files
virtual void writeDesignVars();
//- The writeData function required by the regIOobject write operation
virtual bool writeData(Ostream&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

247
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/betaMax/Darcy/betaMaxDarcy.C Normal file
View File

@ -0,0 +1,247 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "betaMaxDarcy.H"
#include "EdgeMap.H"
#include "syncTools.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(betaMaxDarcy, 0);
addToRunTimeSelectionTable(betaMax, betaMaxDarcy, dictionary);
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
Foam::scalar Foam::betaMaxDarcy::computeLength(const dictionary& dict) const
{
scalar length = Zero;
// If length is not provided explicitly, loop over the provided patches
// and compute the hydraulic diamater
const dictionary& DarcyDict = dict.subDict(type() + "Coeffs");
if (!DarcyDict.readIfPresent("length", length))
{
const labelHashSet inletPatches =
mesh_.boundaryMesh().patchSet
(
DarcyDict.get<wordRes>("inletPatches")
);
// If 2D, use the inlet area divided by the depth in the empty direction
if (mesh_.nGeometricD() != label(3))
{
// Accumulate area
for (const label pI : inletPatches)
{
const fvPatch& patch = mesh_.boundary()[pI];
length += gSum(patch.magSf());
}
// Divide with the span in the empty direction
const Vector<label>& geometricD = mesh_.geometricD();
const boundBox& bounds = mesh_.bounds();
forAll(geometricD, iDir)
{
if (geometricD[iDir] == -1)
{
length /= bounds.span()[iDir];
}
}
}
// If 3D, use the inlet hydraulic diameter
else
{
scalar perimeter = Zero;
scalar area = Zero;
for (const label pI : inletPatches)
{
const polyPatch& patch = mesh_.boundaryMesh()[pI];
// Accumulate perimeter
const edgeList& edges = patch.edges();
const vectorField& points = patch.localPoints();
const label nInternalEdges = patch.nInternalEdges();
const label nEdges = patch.nEdges();
// Processor edges should not be accounted for
boolList isProcessorEdge = markProcessorEdges(patch);
label nActiveEdges(0);
forAll(isProcessorEdge, beI)
{
if (!isProcessorEdge[beI])
{
perimeter += edges[nInternalEdges + beI].mag(points);
nActiveEdges++;
}
}
if (debug > 1)
{
Info<< "patch:: " << patch.name() << nl
<< "Number of boundary edges "
<< returnReduce(nEdges - nInternalEdges, sumOp<label>())
<< ", Number of non-processor edges "
<< returnReduce(nActiveEdges, sumOp<label>()) << endl;
}
// Accumulate area
area += sum(patch.magFaceAreas());
}
reduce(perimeter, sumOp<scalar>());
reduce(area, sumOp<scalar>());
length = scalar(4)*area/perimeter;
}
}
return length;
}
Foam::boolList Foam::betaMaxDarcy::markProcessorEdges
(
const polyPatch& patch
) const
{
const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
const label nNonProcessor = pbm.nNonProcessor();
// Processor edges will artificially increase the perimeter of the inlet.
// We need to purge them.
// Mark all edges connected to a processor patch
label nProcEdges(0);
for (label procI = nNonProcessor; procI < pbm.size() ; ++procI)
{
const polyPatch& procPatch = pbm[procI];
nProcEdges += procPatch.nEdges() - procPatch.nInternalEdges();
}
EdgeMap<bool> isInletEdge(nProcEdges);
for (label procI = nNonProcessor; procI < pbm.size() ; ++procI)
{
const polyPatch& procPatch = pbm[procI];
const labelList& procMp = procPatch.meshPoints();
const label procInternalEdges = procPatch.nInternalEdges();
const label procEdges = procPatch.nEdges();
for (label edgeI = procInternalEdges; edgeI < procEdges; ++edgeI)
{
const edge& e = procPatch.edges()[edgeI];
const edge meshE = edge(procMp[e[0]], procMp[e[1]]);
isInletEdge.insert(meshE, false);
}
}
// Loop over inlet edges to mark common edges with processor patches
const label nInternalEdges = patch.nInternalEdges();
const label nEdges = patch.nEdges();
const edgeList& edges = patch.edges();
const labelList& mp = patch.meshPoints();
for (label edgeI = nInternalEdges; edgeI < nEdges; ++edgeI)
{
const edge& e = edges[edgeI];
const edge meshE = edge(mp[e[0]], mp[e[1]]);
auto iter = isInletEdge.find(meshE);
if (iter.found())
{
iter.val() = true;
}
}
// A special case is a processor patch intersecting the inlet patches on a
// (true) boundary edge of the latter. Use syncEdgeMap to make sure these
// edges don't make it to the final list
syncTools::syncEdgeMap
(
pbm.mesh(),
isInletEdge,
andEqOp<bool>()
);
// Now loop over the inlet faces once more and mark the common edges with
// processor boundaries
boolList isProcessorEdge(nEdges - nInternalEdges, false);
for (label edgeI = nInternalEdges; edgeI < nEdges; ++edgeI)
{
const edge& e = edges[edgeI];
const edge meshE = edge(mp[e[0]], mp[e[1]]);
const auto iter = isInletEdge.cfind(meshE);
if (iter.found() && iter.val())
{
isProcessorEdge[edgeI - nInternalEdges] = true;
}
}
return isProcessorEdge;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::betaMaxDarcy::betaMaxDarcy
(
const fvMesh& mesh,
const dictionary& dict
)
:
betaMax(mesh, dict),
DarcyNumber_
(
dict.subDict(type() + "Coeffs").getOrDefault<scalar>("Da", 1.e-5)
),
length_(computeLength(dict))
{
scalar nu =
IOdictionary
(
IOobject
(
"transportProperties",
mesh.time().constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
)
).get<dimensionedScalar>("nu").value();
value_ = nu/DarcyNumber_/length_/length_;
DebugInfo
<< "Computed a betaMax value of " << value_
<< " based on a length of " << length_ << nl << endl;
}
// ************************************************************************* //

127
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/betaMax/Darcy/betaMaxDarcy.H Normal file
View File

@ -0,0 +1,127 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::betaMaxDarcy
Description
Compute betaMax through the definition of the Darcy number, quantifying
the viscous-to-porous forces ratio
Da = nu/betaMax/L/L
where nu is the bulk viscosity and L is a characteristic length.
The latter is either supplied directly or computed as the
- 2D: area of the inlet patches divided by the span in the empty direction
- 3D: the hydraulic diameter of the inlet patches
SourceFiles
betaMaxDarcy.C
\*---------------------------------------------------------------------------*/
#ifndef betaMaxDarcy_H
#define betaMaxDarcy_H
#include "betaMax.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class betaMaxDarcy Declaration
\*---------------------------------------------------------------------------*/
class betaMaxDarcy
:
public betaMax
{
private:
// Private Member Functions
//- No copy construct
betaMaxDarcy(const betaMaxDarcy&) = delete;
//- No copy assignment
void operator=(const betaMaxDarcy&) = delete;
protected:
// Protected Data
//- The Darcy number expressing the ratio of viscous to porous forces
scalar DarcyNumber_;
//- Characteristic length of the case
// Either supplied directly or computed as the hydraulic diameter of
// the inlet
scalar length_;
// Protected Member Functions
//- Compute the characteristic length
scalar computeLength(const dictionary& dict) const;
//- Mark all common inlet - processor edges
boolList markProcessorEdges(const polyPatch& patch) const;
public:
//- Runtime type information
TypeName("Darcy");
// Constructors
//- Construct from components
betaMaxDarcy
(
const fvMesh& mesh,
const dictionary& dict
);
//- Destructor
virtual ~betaMaxDarcy() = default;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

90
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/betaMax/betaMax/betaMax.C Normal file
View File

@ -0,0 +1,90 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "betaMax.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(betaMax, 0);
defineRunTimeSelectionTable(betaMax, dictionary);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::betaMax::betaMax
(
const fvMesh& mesh,
const dictionary& dict
)
:
mesh_(mesh),
value_(dict.getOrDefault<scalar>("betaMax", Zero))
{}
// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
Foam::autoPtr<Foam::betaMax> Foam::betaMax::New
(
const fvMesh& mesh,
const dictionary& dict
)
{
const word modelType(dict.getOrDefault<word>("betaMaxType", "value"));
auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
Info<< "betaMax type " << modelType << endl;
if (!cstrIter.found())
{
FatalIOErrorInLookup
(
dict,
"betaMaxType",
modelType,
*dictionaryConstructorTablePtr_
) << exit(FatalIOError);
}
return autoPtr<betaMax>(cstrIter()(mesh, dict));
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::scalar Foam::betaMax::value() const
{
return value_;
}
// ************************************************************************* //

138
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/betaMax/betaMax/betaMax.H Normal file
View File

@ -0,0 +1,138 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::betaMax
Description
Base class for selecting the betaMax value, i.e. the value multiplying the
Brinkman penalisation term, for topology optimisation.
SourceFiles
betaMax.C
\*---------------------------------------------------------------------------*/
#ifndef betaMax_H
#define betaMax_H
#include "fvMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class betaMax Declaration
\*---------------------------------------------------------------------------*/
class betaMax
{
private:
// Private Member Functions
//- No copy construct
betaMax(const betaMax&) = delete;
//- No copy assignment
void operator=(const betaMax&) = delete;
protected:
// Protected Data
//- Reference to mesh
const fvMesh& mesh_;
//- betaMax value
scalar value_;
public:
//- Runtime type information
TypeName("betaMax");
// Declare run-time constructor selection table
declareRunTimeSelectionTable
(
autoPtr,
betaMax,
dictionary,
(
const fvMesh& mesh,
const dictionary& dict
),
(mesh, dict)
);
// Constructors
//- Construct from components
betaMax
(
const fvMesh& mesh,
const dictionary& dict
);
// Selectors
//- Construct and return the selected betaMax model
static autoPtr<betaMax> New
(
const fvMesh& mesh,
const dictionary& dict
);
//- Destructor
virtual ~betaMax() = default;
// Member Functions
//- Get value
virtual scalar value() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

79
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/betaMax/stepRamp/betaMaxStepRamp.C Normal file
View File

@ -0,0 +1,79 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "betaMaxStepRamp.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(betaMaxStepRamp, 0);
addToRunTimeSelectionTable
(
betaMax,
betaMaxStepRamp,
dictionary
);
}
// * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * * //
Foam::betaMaxStepRamp::betaMaxStepRamp
(
const fvMesh& mesh,
const dictionary& dict
)
:
betaMax(mesh, dict),
betaMin_(dict.getOrDefault<scalar>("betaMin", Zero)),
funcPtr_(nullptr)
{
funcPtr_.reset
(Function1<scalar>::New("betaMaxStepRamp", dict, "stepRamp"));
DebugInfo
<< "betaMaxStepRamp:: creating interpolation function of type "
<< funcPtr_().type() << endl;
}
// * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * * //
Foam::scalar Foam::betaMaxStepRamp::value() const
{
const scalar t = mesh_.time().timeOutputValue();
const scalar value = betaMin_ + (value_ - betaMin_)*funcPtr_().value(t);
DebugInfo
<< "stepRamp betaMax:: t, betaMax value "
<< t << ", " << value << endl;
return value;
}
// ************************************************************************* //

112
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/betaMax/stepRamp/betaMaxStepRamp.H Normal file
View File

@ -0,0 +1,112 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::betaMaxStepRamp
Description
Computes the betaMax using a stepRamp function which updates itself
based on the current time
\*---------------------------------------------------------------------------*/
#ifndef betaMaxStepRamp_H
#define betaMaxStepRamp_H
#include "betaMax.H"
#include "stepRamp.H"
#include "Scale.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class betaMaxStepRamp Declaration
\*---------------------------------------------------------------------------*/
class betaMaxStepRamp
:
public betaMax
{
private:
// Private Member Functions
//- No copy construct
betaMaxStepRamp(const betaMaxStepRamp&) = delete;
//- No copy assignment
void operator=(const betaMaxStepRamp&) = delete;
protected:
// Protected Data Members
//- Minimum betaMax value (to be used at the first opt cycles)
scalar betaMin_;
//- The step ramp function
autoPtr<Function1<scalar>> funcPtr_;
public:
//- Runtime type information
TypeName("stepRamp");
// Constructors
//- Construct from components
betaMaxStepRamp
(
const fvMesh& mesh,
const dictionary& dict
);
//- Destructor
virtual ~betaMaxStepRamp() = default;
// Member Functions
//- Get value
virtual scalar value() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

53
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/betaMax/value/betaMaxValue.C Normal file
View File

@ -0,0 +1,53 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "betaMaxValue.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(betaMaxValue, 0);
addToRunTimeSelectionTable(betaMax, betaMaxValue, dictionary);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::betaMaxValue::betaMaxValue
(
const fvMesh& mesh,
const dictionary& dict
)
:
betaMax(mesh, dict)
{}
// ************************************************************************* //

98
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/betaMax/value/betaMaxValue.H Normal file
View File

@ -0,0 +1,98 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::betaMaxValue
Description
Read the betaMax value directly.
SourceFiles
betaMaxValue.C
\*---------------------------------------------------------------------------*/
#ifndef betaMaxValue_H
#define betaMaxValue_H
#include "betaMax.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class betaMaxValue Declaration
\*---------------------------------------------------------------------------*/
class betaMaxValue
:
public betaMax
{
private:
// Private Member Functions
//- No copy construct
betaMaxValue(const betaMaxValue&) = delete;
//- No copy assignment
void operator=(const betaMaxValue&) = delete;
public:
//- Runtime type information
TypeName("value");
// Constructors
//- Construct from components
betaMaxValue
(
const fvMesh& mesh,
const dictionary& dict
);
//- Destructor
virtual ~betaMaxValue() = default;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

119
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/dynamicTopODesignVariables/dynamicTopODesignVariables.C Normal file
View File

@ -0,0 +1,119 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "dynamicTopODesignVariables.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(dynamicTopODesignVariables, 0);
addToRunTimeSelectionTable
(
designVariables,
dynamicTopODesignVariables,
designVariables
);
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
void Foam::dynamicTopODesignVariables::setActiveDesignVariables
(
const label fluidID,
const bool activeIO
)
{
cellZoneMesh& cellZones = mesh_.cellZones();
marchCells_.addFixedCells(cellZones, zones_.fixedPorousZoneIDs());
marchCells_.addFixedCells(cellZones, zones_.fixedZeroPorousZoneIDs());
if (!activeIO)
{
marchCells_.addFixedCells(zones_.IOCells());
}
marchCells_.update();
activeDesignVariables_ = marchCells_.getActiveCells();
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::dynamicTopODesignVariables::dynamicTopODesignVariables
(
fvMesh& mesh,
const dictionary& dict
)
:
dynamicTopODesignVariables(mesh, dict, mesh.nCells())
{}
Foam::dynamicTopODesignVariables::dynamicTopODesignVariables
(
fvMesh& mesh,
const dictionary& dict,
const label size
)
:
topODesignVariables(mesh, dict, size),
marchCells_(mesh, dict.subDict("marchingCoeffs"))
{
// Rest of the contrsuctor initialization
initialize();
}
// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
Foam::autoPtr<Foam::dynamicTopODesignVariables>
Foam::dynamicTopODesignVariables::New
(
fvMesh& mesh,
const dictionary& dict
)
{
return autoPtr<dynamicTopODesignVariables>
(
new dynamicTopODesignVariables(mesh, dict)
);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * //
void Foam::dynamicTopODesignVariables::evolveNumber()
{
marchCells_.update();
activeDesignVariables_ = marchCells_.getActiveCells();
}
// ************************************************************************* //

146
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/dynamicTopODesignVariables/dynamicTopODesignVariables.H Normal file
View File

@ -0,0 +1,146 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::dynamicTopODesignVariables
Description
Design variables for topology optimisation problems.
Active design variables are added gradually, to help the optimiser from
being trapped to a sub-optimal solution track in the first optimisation
cycles.
SourceFiles
dynamicTopODesignVariables.C
\*---------------------------------------------------------------------------*/
#ifndef dynamicTopODesignVariables_H
#define dynamicTopODesignVariables_H
#include "topODesignVariables.H"
#include "marchingCells.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class dynamicTopODesignVariables Declaration
\*---------------------------------------------------------------------------*/
class dynamicTopODesignVariables
:
public topODesignVariables
{
protected:
// Protected data
//- Mechanism for gradually activating design variables
marchingCells marchCells_;
// Protected Member Functions
//- Set active design variables
virtual void setActiveDesignVariables
(
const label fluidID = 0,
const bool activeIO = false
);
private:
// Private Member Functions
//- Disallow default bitwise copy construct
dynamicTopODesignVariables
(
const dynamicTopODesignVariables&
) = delete;
//- Disallow default bitwise assignment
void operator=(const dynamicTopODesignVariables&) = delete;
public:
//- Runtime type information
TypeName("dynamicTopO");
// Constructors
//- Construct from dictionary
dynamicTopODesignVariables
(
fvMesh& mesh,
const dictionary& dict
);
//- Construct from dictionary and size
dynamicTopODesignVariables
(
fvMesh& mesh,
const dictionary& dict,
const label size
);
// Selectors
//- Construct and return the selected design variables
static autoPtr<dynamicTopODesignVariables> New
(
fvMesh& mesh,
const dictionary& dict
);
//- Destructor
virtual ~dynamicTopODesignVariables() = default;
// Member Functions
//- Update the active design variables
virtual void evolveNumber();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

94
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/interpolationFunctions/BorrvallPetersson/BorrvallPeterssonInterpolation.C Normal file
View File

@ -0,0 +1,94 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "BorrvallPeterssonInterpolation.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(BorrvallPeterssonInterpolation, 0);
addToRunTimeSelectionTable
(
topOInterpolationFunction,
BorrvallPeterssonInterpolation,
dictionary
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::BorrvallPeterssonInterpolation::BorrvallPeterssonInterpolation
(
const fvMesh& mesh,
const dictionary& dict
)
:
topOInterpolationFunction(mesh, dict),
b_(Function1<scalar>::New("b", dict))
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::BorrvallPeterssonInterpolation::interpolate
(
const scalarField& arg,
scalarField& res
) const
{
const scalar time(mesh_.time().timeOutputValue());
const scalar t(time == 0 ? 1. : time);
const scalar b(b_->value(t));
DebugInfo
<< type() << "::interpolate:: t, b value " << t << " " << b << endl;
res = arg/(scalar(1) + b*(scalar(1) - arg));
}
Foam::tmp<Foam::scalarField>
Foam::BorrvallPeterssonInterpolation::derivative(const scalarField& arg) const
{
tmp<scalarField> tderiv(tmp<scalarField>::New(arg.size(), Zero));
scalarField& deriv = tderiv.ref();
const scalar t(mesh_.time().timeOutputValue());
const scalar b(b_->value(t));
DebugInfo
<< type() << "::derivative:: t, b " << t << " " << b << endl;
deriv = (b + scalar(1))/sqr(scalar(1) + b*(scalar(1) - arg));
return tderiv;
}
// ************************************************************************* //

130
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/interpolationFunctions/BorrvallPetersson/BorrvallPeterssonInterpolation.H Normal file
View File

@ -0,0 +1,130 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::BorrvallPeterssonInterpolation
Description
Computes the Brickmann penalization function for topO optimisation,
proposed by Borrvall and Petersson.
Reference:
\verbatim
Borrvall, T., & Petersson, J. (2002).
TopO optimization of fluids in Stokes flow.
International Journal for Numerical Methods in Fluids, 41(1), 77-107.
https://doi.org/10.1002/fld.426
\endverbatim
SourceFiles
BorrvallPeterssonInterpolation.C
\*---------------------------------------------------------------------------*/
#ifndef BorrvallPeterssonInterpolation_H
#define BorrvallPeterssonInterpolation_H
#include "topOInterpolationFunction.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class BorrvallPeterssonInterpolation Declaration
\*---------------------------------------------------------------------------*/
class BorrvallPeterssonInterpolation
:
public topOInterpolationFunction
{
private:
// Private Data
//- Steepness parameter
// High b values lead to more steep interpolations - higher contrasts.
// May vary in time
autoPtr<Function1<scalar>> b_;
// Private Member Functions
//- No copy construct
BorrvallPeterssonInterpolation
(
const BorrvallPeterssonInterpolation&
) = delete;
//- No copy assignment
void operator=(const BorrvallPeterssonInterpolation&) = delete;
public:
//- Runtime type information
TypeName("BorrvallPetersson");
// Constructors
//- Construct from components
BorrvallPeterssonInterpolation
(
const fvMesh& mesh,
const dictionary& dict
);
//- Destructor
virtual ~BorrvallPeterssonInterpolation() = default;
// Member Functions
//- Interpolate argument and write to result
virtual void interpolate
(
const scalarField& arg,
scalarField& res
) const;
//- Return of function with respect to the argument field
virtual tmp<scalarField> derivative(const scalarField& arg) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

91
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/interpolationFunctions/SIMP/SIMPInterpolation.C Normal file
View File

@ -0,0 +1,91 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Interpolation CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "SIMPInterpolation.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(SIMPInterpolation, 0);
addToRunTimeSelectionTable
(
topOInterpolationFunction,
SIMPInterpolation,
dictionary
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::SIMPInterpolation::SIMPInterpolation
(
const fvMesh& mesh,
const dictionary& dict
)
:
topOInterpolationFunction(mesh, dict),
b_(Function1<scalar>::New("b",dict))
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::SIMPInterpolation::interpolate
(
const scalarField& arg,
scalarField& res
) const
{
// Compute exponent
const scalar time(mesh_.time().timeOutputValue());
const scalar t(time == 0 ? 1. : time);
const scalar b(b_->value(t));
res = pow(arg, b);
}
Foam::tmp<Foam::scalarField>
Foam::SIMPInterpolation::derivative(const scalarField& arg) const
{
tmp<scalarField> tderiv(tmp<scalarField>::New(arg.size(), Zero));
scalarField& deriv = tderiv.ref();
const scalar t(mesh_.time().timeOutputValue());
const scalar b(b_->value(t));
deriv = b*pow(arg, b - scalar(1));
return tderiv;
}
// ************************************************************************* //

125
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/interpolationFunctions/SIMP/SIMPInterpolation.H Normal file
View File

@ -0,0 +1,125 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Interpolation CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::SIMPInterpolation
Description
Computes the Brickmann penalization function for topology optimisation,
using the SIMP approach proposed by Bendsøe
Reference:
\verbatim
Bendsøe, M. P. (1989)
Optimal shape design as a material distribution problem.
Structural optimization, 1(4), 193–202.
https://doi.org/10.1007/BF01650949
\endverbatim
InterpolationFiles
SIMPInterpolation.C
\*---------------------------------------------------------------------------*/
#ifndef SIMPInterpolation_H
#define SIMPInterpolation_H
#include "topOInterpolationFunction.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class SIMPInterpolation Declaration
\*---------------------------------------------------------------------------*/
class SIMPInterpolation
:
public topOInterpolationFunction
{
// Private Data
//- Steepness parameter
// High b values lead to more steep interpolations - higher contrasts.
// May vary in time
autoPtr<Function1<scalar>> b_;
private:
// Private Member Functions
//- No copy construct
SIMPInterpolation(const SIMPInterpolation&);
//- No copy assignment
void operator=(const SIMPInterpolation&);
public:
//- Runtime type information
TypeName("SIMP");
// Constructors
//- Construct from components
SIMPInterpolation
(
const fvMesh& mesh,
const dictionary& dict
);
//- Destructor
virtual ~SIMPInterpolation() = default;
// Member Functions
//- Interpolate argument and write to result
virtual void interpolate
(
const scalarField& arg,
scalarField& res
) const;
//- Return of function with respect to the argument field
virtual tmp<scalarField> derivative(const scalarField& arg) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

95
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/interpolationFunctions/exp/expInterpolation.C Normal file
View File

@ -0,0 +1,95 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "expInterpolation.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(expInterpolation, 1);
addToRunTimeSelectionTable
(
topOInterpolationFunction,
expInterpolation,
dictionary
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::expInterpolation::expInterpolation
(
const fvMesh& mesh,
const dictionary& dict
)
:
topOInterpolationFunction(mesh, dict),
b_(Function1<scalar>::New("b", dict))
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::expInterpolation::interpolate
(
const scalarField& arg,
scalarField& res
) const
{
const scalar time(mesh_.time().timeOutputValue());
const scalar t(time == 0 ? 1. : time);
const scalar b(b_->value(t));
DebugInfo
<< type() << "::interpolate:: t, b value " << t << " " << b << endl;
res = exp(-b*(1.-arg)) - (1. - arg)*exp(-b);
}
Foam::tmp<Foam::scalarField> Foam::expInterpolation::derivative
(
const scalarField& arg
) const
{
tmp<scalarField> tderiv(tmp<scalarField>::New(arg.size(), Zero));
scalarField& deriv = tderiv.ref();
const scalar t(mesh_.time().timeOutputValue());
const scalar b(b_->value(t));
DebugInfo
<< type() << "derivative:: t, b " << t << " " << b << endl;
deriv = b*exp(-b*(1. - arg)) + exp(-b);
return tderiv;
}
// ************************************************************************* //

129
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/interpolationFunctions/exp/expInterpolation.H Normal file
View File

@ -0,0 +1,129 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::expInterpolation
Description
Computes the Brickmann penalization function for topology optimisation,
relying on exp functions and proposed in
Reference:
\verbatim
Wang, F., Lazarov, B., Sigmund, O.(2011).
On projection methods, convergence and robust formulations in topology
optimization.
Structural and Multidisciplinary Optimization 43(6), 767–784.
https://doi:10.1007/s00158-010-0602-y
\endverbatim
SourceFiles
expInterpolation.C
\*---------------------------------------------------------------------------*/
#ifndef expInterpolation_H
#define expInterpolation_H
#include "topOInterpolationFunction.H"
#include "Function1.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class expInterpolation Declaration
\*---------------------------------------------------------------------------*/
class expInterpolation
:
public topOInterpolationFunction
{
private:
// Private Data
//- Steepness parameter
// High b values lead to more steep interpolations - higher contrasts.
// May vary in time
autoPtr<Function1<scalar>> b_;
// Private Member Functions
//- No copy construct
expInterpolation(const expInterpolation&);
//- No copy assignment
void operator=(const expInterpolation&);
public:
//- Runtime type information
TypeName("exp");
// Constructors
//- Construct from components
expInterpolation
(
const fvMesh& mesh,
const dictionary& dict
);
//- Destructor
virtual ~expInterpolation() = default;
// Member Functions
//- Interpolate argument and write to result
virtual void interpolate
(
const scalarField& arg,
scalarField& res
) const;
//- Return of function with respect to the argument field
virtual tmp<scalarField> derivative(const scalarField& arg) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

96
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/interpolationFunctions/interpolationFunction/topOInterpolationFunction.C Normal file
View File

@ -0,0 +1,96 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "topOInterpolationFunction.H"
#include "runTimeSelectionTables.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(topOInterpolationFunction, 1);
defineRunTimeSelectionTable(topOInterpolationFunction, dictionary);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::topOInterpolationFunction::topOInterpolationFunction
(
const fvMesh& mesh,
const dictionary& dict
)
:
mesh_(mesh),
dict_(dict)
{}
// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
Foam::autoPtr<Foam::topOInterpolationFunction>
Foam::topOInterpolationFunction::New
(
const fvMesh& mesh,
const dictionary& dict
)
{
const word type = dict.get<word>("function");
Info<< "topOInterpolationFunction type : " << type << endl;
auto cstrIter = dictionaryConstructorTablePtr_->cfind(type);
if (!cstrIter.found())
{
FatalIOErrorInLookup
(
dict,
"function",
type,
*dictionaryConstructorTablePtr_
) << exit(FatalIOError);
}
return autoPtr<topOInterpolationFunction>(cstrIter()(mesh, dict));
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::topOInterpolationFunction::setLengthScaleParam
(
const scalar lengthScale
)
{
WarningInFunction
<< "Function " << type() << " has no length scale parameter" << endl;
}
// ************************************************************************* //

149
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/interpolationFunctions/interpolationFunction/topOInterpolationFunction.H Normal file
View File

@ -0,0 +1,149 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
topOInterpolationFunction
Description
Abstract base class for interpolation functions used in topology
optimisation
SourceFiles
topOInterpolationFunction.C
\*---------------------------------------------------------------------------*/
#ifndef topOInterpolationFunction_H
#define topOInterpolationFunction_H
#include "fvMesh.H"
#include "Function1.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class topOInterpolationFunction Declaration
\*---------------------------------------------------------------------------*/
class topOInterpolationFunction
{
protected:
// Protected data
const fvMesh& mesh_;
dictionary dict_;
private:
// Private Member Functions
//- No copy construct
topOInterpolationFunction(const topOInterpolationFunction&);
//- No copy assignment
void operator=(const topOInterpolationFunction&);
public:
//- Runtime type information
TypeName("topOInterpolationFunction");
// Declare run-time constructor selection table
declareRunTimeSelectionTable
(
autoPtr,
topOInterpolationFunction,
dictionary,
(
const fvMesh& mesh,
const dictionary& dict
),
(mesh, dict)
);
// Constructors
//- Construct from components
topOInterpolationFunction
(
const fvMesh& mesh,
const dictionary& dict
);
// Selectors
//- Return an autoPtr to the selected interpolation type
static autoPtr<topOInterpolationFunction> New
(
const fvMesh& mesh,
const dictionary& dict
);
//- Destructor
virtual ~topOInterpolationFunction() = default;
// Member Functions
//- Interpolate argument to result
virtual void interpolate
(
const scalarField& arg,
scalarField& res
) const = 0;
//- Return of function with respect to the argument field
virtual tmp<scalarField> derivative(const scalarField& arg) const = 0;
//- Set the parameter determining length scale
// Warning if the selected interpolation function has no such param
virtual void setLengthScaleParam(const scalar lengthScale);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

94
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/interpolationFunctions/invBP/invBP.C Normal file
View File

@ -0,0 +1,94 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "invBP.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(invBP, 0);
addToRunTimeSelectionTable
(
topOInterpolationFunction,
invBP,
dictionary
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::invBP::invBP
(
const fvMesh& mesh,
const dictionary& dict
)
:
topOInterpolationFunction(mesh, dict),
b_(Function1<scalar>::New("b", dict))
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::invBP::interpolate
(
const scalarField& arg,
scalarField& res
) const
{
const scalar time(mesh_.time().timeOutputValue());
const scalar t(time == 0 ? 1. : time);
const scalar b(b_->value(t));
DebugInfo
<< type() << "::interpolate:: t, b value " << t << " " << b << endl;
res = arg*(b + 1)/( b*arg +1);
}
Foam::tmp<Foam::scalarField>
Foam::invBP::derivative(const scalarField& arg) const
{
tmp<scalarField> tderiv(tmp<scalarField>::New(arg.size(), Zero));
scalarField& deriv = tderiv.ref();
const scalar t(mesh_.time().timeOutputValue());
const scalar b(b_->value(t));
DebugInfo
<< type() << "::derivative:: t, b " << t << " " << b << endl;
deriv = (b + scalar(1))/sqr(scalar(1) + b*arg);
return tderiv;
}
// ************************************************************************* //

127
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/interpolationFunctions/invBP/invBP.H Normal file
View File

@ -0,0 +1,127 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::invBP
Description
Computes the Brickmann penalization function for topology optimisation,
based in the inverse of the function proposed by Borrvall and Petersson.
Reference:
\verbatim
Borrvall, T., & Petersson, J. (2002).
TopO optimization of fluids in Stokes flow.
International Journal for Numerical Methods in Fluids, 41(1), 77-107.
https://doi.org/10.1002/fld.426
\endverbatim
SourceFiles
invBP.C
\*---------------------------------------------------------------------------*/
#ifndef invBP_H
#define invBP_H
#include "topOInterpolationFunction.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class invBP Declaration
\*---------------------------------------------------------------------------*/
class invBP
:
public topOInterpolationFunction
{
private:
// Private Data
//- Steepness parameter
// High b values lead to more steep interpolations - higher contrasts.
// May vary in time
autoPtr<Function1<scalar>> b_;
// Private Member Functions
//- No copy construct
invBP(const invBP&);
//- No copy assignment
void operator=(const invBP&);
public:
//- Runtime type information
TypeName("invBP");
// Constructors
//- Construct from components
invBP
(
const fvMesh& mesh,
const dictionary& dict
);
//- Destructor
virtual ~invBP() = default;
// Member Functions
//- Interpolate argument and write to result
virtual void interpolate
(
const scalarField& arg,
scalarField& res
) const;
//- Return of function with respect to the argument field
virtual tmp<scalarField> derivative(const scalarField& arg) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

77
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/interpolationFunctions/linear/linearInterpolation.C Normal file
View File

@ -0,0 +1,77 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Interpolation CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "linearInterpolation.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(linearInterpolation, 0);
addToRunTimeSelectionTable
(
topOInterpolationFunction,
linearInterpolation,
dictionary
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::linearInterpolation::linearInterpolation
(
const fvMesh& mesh,
const dictionary& dict
)
:
topOInterpolationFunction(mesh, dict)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::linearInterpolation::interpolate
(
const scalarField& arg,
scalarField& res
) const
{
res = arg;
}
Foam::tmp<Foam::scalarField>
Foam::linearInterpolation::derivative(const scalarField& arg) const
{
return tmp<scalarField>::New(arg.size(), scalar(1));
}
// ************************************************************************* //

110
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/interpolationFunctions/linear/linearInterpolation.H Normal file
View File

@ -0,0 +1,110 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Interpolation CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::linearInterpolation
Description
Implements a linear (with respect to beta) term for topology optimisation
InterpolationFiles
linearInterpolation.C
\*---------------------------------------------------------------------------*/
#ifndef linearInterpolation_H
#define linearInterpolation_H
#include "topOInterpolationFunction.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class linearInterpolation Declaration
\*---------------------------------------------------------------------------*/
class linearInterpolation
:
public topOInterpolationFunction
{
private:
// Private Member Functions
//- No copy construct
linearInterpolation(const linearInterpolation&);
//- No copy assignment
void operator=(const linearInterpolation&);
public:
//- Runtime type information
TypeName("linear");
// Constructors
//- Construct from components
linearInterpolation
(
const fvMesh& mesh,
const dictionary& dict
);
//- Destructor
virtual ~linearInterpolation() = default;
// Member Functions
//- Interpolate argument and write to result
virtual void interpolate
(
const scalarField& arg,
scalarField& res
) const;
//- Return of function with respect to the argument field
virtual tmp<scalarField> derivative(const scalarField& arg) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

90
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/interpolationFunctions/sinh/sinhInterpolation.C Normal file
View File

@ -0,0 +1,90 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Interpolation CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "sinhInterpolation.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(sinhInterpolation, 0);
addToRunTimeSelectionTable
(
topOInterpolationFunction,
sinhInterpolation,
dictionary
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sinhInterpolation::sinhInterpolation
(
const fvMesh& mesh,
const dictionary& dict
)
:
topOInterpolationFunction(mesh, dict),
b_(Function1<scalar>::New("b",dict))
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::sinhInterpolation::interpolate
(
const scalarField& arg,
scalarField& res
) const
{
const scalar time(mesh_.time().timeOutputValue());
const scalar t(time == 0 ? 1. : time);
const scalar b(b_->value(t));
res = scalar(1) - sinh(b*(scalar(1) - arg))/sinh(b);
}
Foam::tmp<Foam::scalarField>
Foam::sinhInterpolation::derivative(const scalarField& arg) const
{
tmp<scalarField> tderiv(tmp<scalarField>::New(arg.size(), Zero));
scalarField& deriv = tderiv.ref();
const scalar t(mesh_.time().timeOutputValue());
const scalar b(b_->value(t));
deriv = b*cosh(b*(scalar(1) - arg))/sinh(b);
return tderiv;
}
// ************************************************************************* //

116
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/interpolationFunctions/sinh/sinhInterpolation.H Normal file
View File

@ -0,0 +1,116 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Interpolation CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::sinhInterpolation
Description
Updates the porosity values in topology optimization using a sinh function
InterpolationFiles
sinhInterpolation.C
\*---------------------------------------------------------------------------*/
#ifndef sinhInterpolation_H
#define sinhInterpolation_H
#include "topOInterpolationFunction.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class sinhInterpolation Declaration
\*---------------------------------------------------------------------------*/
class sinhInterpolation
:
public topOInterpolationFunction
{
private:
// Private Data
//- Steepness parameter
// High b values lead to more steep interpolations - higher contrasts.
// May vary in time
autoPtr<Function1<scalar>> b_;
// Private Member Functions
//- No copy construct
sinhInterpolation(const sinhInterpolation&);
//- No copy assignment
void operator=(const sinhInterpolation&);
public:
//- Runtime type information
TypeName("sinh");
// Constructors
//- Construct from components
sinhInterpolation
(
const fvMesh& mesh,
const dictionary& dict
);
//-Destructor
virtual ~sinhInterpolation() = default;
// Member Functions
//- Interpolate argument and write to result
virtual void interpolate
(
const scalarField& arg,
scalarField& res
) const;
//- Return of function with respect to the argument field
virtual tmp<scalarField> derivative(const scalarField& arg) const;
};
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

108
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/interpolationFunctions/tanh/tanhInterpolation.C Normal file
View File

@ -0,0 +1,108 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "tanhInterpolation.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(tanhInterpolation, 0);
addToRunTimeSelectionTable
(
topOInterpolationFunction,
tanhInterpolation,
dictionary
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::tanhInterpolation::tanhInterpolation
(
const fvMesh& mesh,
const dictionary& dict
)
:
topOInterpolationFunction(mesh, dict),
b_(Function1<scalar>::New("b", dict)),
eta_(dict.getOrDefault<scalar>("eta", scalar(0.5)))
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::tanhInterpolation::interpolate
(
const scalarField& arg,
scalarField& res
) const
{
const scalar time(mesh_.time().timeOutputValue());
const scalar t(time == 0 ? 1. : time);
const scalar b(b_->value(t));
if (debug > 1)
{
Info<< type() << "::interpolate:: t, b value " << t << " " << b
<< " eta " << eta_ << endl;
}
res = (tanh(b*eta_) + tanh(b*(arg - eta_)))
/(tanh(b*eta_) + tanh(b*(scalar(1) - eta_)));
}
Foam::tmp<Foam::scalarField> Foam::tanhInterpolation::derivative
(
const scalarField& arg
) const
{
tmp<scalarField> tderiv(tmp<scalarField>::New(arg.size(), Zero));
scalarField& deriv = tderiv.ref();
const scalar t(mesh_.time().timeOutputValue());
const scalar b(b_->value(t));
DebugInfo
<< type() << "::interpolate:: t, b value " << t << " " << b
<< " eta " << eta_ << endl;
deriv = b*(scalar(1) - sqr(tanh(b*(arg - eta_))))
/(tanh(b*eta_) + tanh(b*(scalar(1) - eta_)));
return tderiv;
}
void Foam::tanhInterpolation::setLengthScaleParam(const scalar lengthScale)
{
eta_ = lengthScale;
}
// ************************************************************************* //

136
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/interpolationFunctions/tanh/tanhInterpolation.H Normal file
View File

@ -0,0 +1,136 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::tanhInterpolation
Description
Computes the Brickmann penalization function for topology optimisation,
relying on tanh functions and proposed in
Reference:
\verbatim
Wang, F., Lazarov, B., Sigmund, O.(2011).
On projection methods, convergence and robust formulations in topO
optimization.
Structural and Multidisciplinary Optimization 43(6), 767–784.
https://doi:10.1007/s00158-010-0602-y
\endverbatim
SourceFiles
tanhInterpolation.C
\*---------------------------------------------------------------------------*/
#ifndef tanhInterpolation_H
#define tanhInterpolation_H
#include "topOInterpolationFunction.H"
#include "Function1.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class tanhInterpolation Declaration
\*---------------------------------------------------------------------------*/
class tanhInterpolation
:
public topOInterpolationFunction
{
private:
// Private Data
//- Steepness parameter
// High b values lead to more steep interpolations - higher contrasts.
// May vary in time
autoPtr<Function1<scalar>> b_;
//- Length scale parameter
// Default to 0.5
scalar eta_;
// Private Member Functions
//- No copy construct
tanhInterpolation(const tanhInterpolation&);
//- No copy assignment
void operator=(const tanhInterpolation&);
public:
//- Runtime type information
TypeName("tanh");
// Constructors
//- Construct from components
tanhInterpolation
(
const fvMesh& mesh,
const dictionary& dict
);
//- Destructor
virtual ~tanhInterpolation() = default;
// Member Functions
//- Interpolate argument and write to result
virtual void interpolate
(
const scalarField& arg,
scalarField& res
) const;
//- Return of function with respect to the argument field
virtual tmp<scalarField> derivative(const scalarField& arg) const;
//- Set the parameter determining length scale
virtual void setLengthScaleParam(const scalar lengthScale);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

262
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/marchingCells/marchingCells.C Normal file
View File

@ -0,0 +1,262 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "marchingCells.H"
#include "fvMeshSubset.H"
#include "cellSet.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(marchingCells, 1);
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
void Foam::marchingCells::initialise()
{
label nChangedFaces(0);
labelList changedFaces(mesh_.nFaces(), -1);
// Gather initial faces from the seeding patches
for (const label patchI : seedPatches_)
{
const fvPatch& patch = mesh_.boundary()[patchI];
const label start = patch.start();
forAll(patch, faceI)
{
changedFaces[nChangedFaces++] = start + faceI;
}
}
// Gather the boundary faces from the seeding cell zones
for (label cellZoneID : seedCellZoneIDs_)
{
const labelList& zoneCells = mesh_.cellZones()[cellZoneID];
// Create sub mesh from the given zone
fvMeshSubset subSetMesh(mesh_, zoneCells);
const fvMesh& subMesh = subSetMesh.subMesh();
// Get the exposed faces of the subMesh and use them as seeds to
// patchWave
const word patchName(fvMeshSubset::exposedPatchName);
const label patchID(subMesh.boundaryMesh().findPatchID(patchName));
const polyPatch& subMeshPatch = subMesh.boundaryMesh()[patchID];
// Map to faces of the original mesh
const labelList& faceMap = subSetMesh.faceMap();
const label start = subMeshPatch.start();
// Go through a primitivePatch field (faceCentres) since the type
// of the patch containing the exposed faces is empty and a zero size
// is returned for the fields accessed through fvPatch
forAll(subMeshPatch.faceCentres(), faceI)
{
changedFaces[nChangedFaces++] = faceMap[start + faceI];
}
}
// Gather faces from the seeding face zones
for (label faceZoneID : seedFaceZoneIDs_)
{
const labelList& zoneFaces = mesh_.faceZones()[faceZoneID];
for (label faceI : zoneFaces)
{
changedFaces[nChangedFaces++] = faceI;
}
}
changedFaces.setSize(nChangedFaces);
List<wallPointData<bool>> changedFacesInfo(nChangedFaces);
const vectorField& faceCentres = mesh_.faceCentres();
forAll(changedFaces, fI)
{
changedFacesInfo[fI] =
wallPointData<bool>(faceCentres[changedFaces[fI]], true, 0.0);
}
meshWave_.setFaceInfo(changedFaces, changedFacesInfo);
// Initialization has been completed
initialised_ = true;
}
void Foam::marchingCells::appendSeedCell(const label cellID)
{
if (!isFixedCell_[cellID])
{
isActiveCell_[cellID] = true;
activeCells_.append(cellID);
}
}
void Foam::marchingCells::march
(
label nVisited,
const label cI,
labelList& newlyAddedCells
)
{
++nVisited;
if (nVisited < marchingStep_ + 1)
{
const labelList& cellCells = mesh_.cellCells()[cI];
for (label neiCellI : cellCells)
{
if (!isFixedCell_[neiCellI])
{
if (!isActiveCell_[neiCellI])
{
isActiveCell_[neiCellI] = true;
newlyAddedCells.append(neiCellI);
}
march(nVisited, neiCellI, newlyAddedCells);
}
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::marchingCells::marchingCells
(
const fvMesh& mesh,
const dictionary& dict
)
:
mesh_(mesh),
seedPatches_
(
mesh_.boundaryMesh().patchSet
(
dict.getOrDefault<wordRes>("seedPatches", wordRes(0))
)
),
seedCellZoneIDs_
(
mesh_.cellZones().indices
(
dict.getOrDefault<wordRes>("seedCellZones", wordRes(0))
)
),
seedFaceZoneIDs_
(
mesh_.faceZones().indices
(
dict.getOrDefault<wordRes>("seedFaceZones", wordRes(0))
)
),
marchingStep_(dict.get<label>("marchingStep")),
isActiveCell_(mesh_.nCells(), false),
isFixedCell_(mesh_.nCells(), false),
activeCells_(0),
addedCells_(0),
initialised_(false),
nIters_(0),
allFaceInfo_(mesh_.nFaces()),
allCellInfo_(mesh.nCells()),
meshWave_(mesh_, allFaceInfo_, allCellInfo_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::marchingCells::update()
{
// If the first seeding has been performed, do it now
if (!initialised_)
{
initialise();
}
// Iterate the meshWave algorithm
meshWave_.iterate(marchingStep_);
// Grab the newly added cells
addedCells_ = labelList(mesh_.nCells(), -1);
label nAddedCells(0);
forAll(allCellInfo_, cI)
{
if (allCellInfo_[cI].data() && !isActiveCell_[cI] && !isFixedCell_[cI])
{
addedCells_[nAddedCells++] = cI;
isActiveCell_[cI] = true;
}
}
addedCells_.setSize(nAddedCells);
// Add newly found cells to the list of activeCells
activeCells_.append(addedCells_);
// Write cell set
if (debug)
{
cellSet activeCellList
(
mesh_,
"activeCells" + name(nIters_),
activeCells_.size()
);
for (label cellI : activeCells_)
{
activeCellList.insert(cellI);
}
activeCellList.write();
}
++nIters_;
}
void Foam::marchingCells::addFixedCells
(
const cellZoneMesh& cellZoneMesh,
const labelList& fixedCellZoneIDs
)
{
for (const label cellZoneID : fixedCellZoneIDs)
{
for (const label cI : cellZoneMesh[cellZoneID])
{
isFixedCell_[cI] = true;
isActiveCell_[cI] = false;
}
}
}
void Foam::marchingCells::addFixedCells(const labelList& fixedCellIDs)
{
for (const label cI : fixedCellIDs)
{
isFixedCell_[cI] = true;
isActiveCell_[cI] = false;
}
}
// ************************************************************************* //

207
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/marchingCells/marchingCells.H Normal file
View File

@ -0,0 +1,207 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
marchingCells
Description
Marches cells from given seed patches, faceZones or the boundary faces of
cellZones, towards the interior of the mesh, with a given step at a time.
Used in topology optimisation, to gradually activate the design variables,
in an attempt to prevent the optimiser from acting too agressively in the
first optimisation cycles. Use meshWave to do the marching.
SourceFiles
marchingCells.C
\*---------------------------------------------------------------------------*/
#ifndef marchingCells_H
#define marchingCells_H
#include "fvMesh.H"
#include "wallPointData.H"
#include "FaceCellWave.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class marchingCells Declaration
\*---------------------------------------------------------------------------*/
class marchingCells
{
protected:
// Protected data
const fvMesh& mesh_;
//- Patches used as seeds in the marching algorithm
labelHashSet seedPatches_;
//- Cell zones, the boundary faces of which are used as seeds in the
//- marching algorithm
labelList seedCellZoneIDs_;
//- Face zones used as seeds in the marching algorithm.
labelList seedFaceZoneIDs_;
//- Marching step
label marchingStep_;
//- Whether each cell is curently active or not
boolList isActiveCell_;
//- Should this cell remain incative
boolList isFixedCell_;
//- Which are the active cells
DynamicList<label> activeCells_;
//- Which are the added cells
labelList addedCells_;
//- Has the initial seeding been conducted
bool initialised_;
//- Iterations conducted thus far
label nIters_;
//- Information for all faces
List<wallPointData<bool>> allFaceInfo_;
//- Information for all cells
List<wallPointData<bool>> allCellInfo_;
//- Engine propagating the active cells.
// Uses infrastructure from FaceCellWave, to easily handle parallel
// interfaces
FaceCellWave<wallPointData<bool>> meshWave_;
// Protected Member Functions
//- Initialise the active cells from the seeding patches
void initialise();
//- Append cell to seed cells
void appendSeedCell(const label cellID);
// Updates the isActiveCells_ list and gathers newlyAddedCells
void march
(
label nVisited,
const label cI,
labelList& newlyAddedCells
);
private:
// Private Member Functions
//- No copy construct
marchingCells(const marchingCells&) = delete;
//- No copy assignment
void operator=(const marchingCells&) = delete;
public:
//- Runtime type information
TypeName("marchingCells");
// Constructors
//- Construct from components
marchingCells
(
const fvMesh& mesh,
const dictionary& dict
);
//- Destructor
virtual ~marchingCells() = default;
// Member Functions
// Access
//- Get the active cells
inline const labelList& getActiveCells() const
{
return activeCells_;
}
//- Number of newly added cells in the last iteration
inline label numberOfAddedCells() const
{
return addedCells_.size();
}
// Edit
//- Update active cells
void update();
//- Set marching step
inline void setMarchingStep(const label step)
{
marchingStep_ = step;
}
//- Add fixed cells through cellZone IDs
void addFixedCells
(
const cellZoneMesh& cellZoneMesh,
const labelList& fixedCellZoneIDs
);
//- Add fixed cells
void addFixedCells(const labelList& fixedCellIDs);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

175
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/regularisation/fieldRegularisation.C Normal file
View File

@ -0,0 +1,175 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2007-2023 PCOpt/NTUA
Copyright (C) 2013-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "fieldRegularisation.H"
#include "fixedValueFvPatchFields.H"
#include "zeroGradientFvPatchFields.H"
#include "wallFvPatch.H"
#include "fvm.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(fieldRegularisation, 1);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::fieldRegularisation::fieldRegularisation
(
fvMesh& mesh,
const scalarField& alpha,
const topOZones& zones,
const dictionary& dict
)
:
mesh_(mesh),
dict_(dict),
zones_(zones),
regularise_(dict.getOrDefault<bool>("regularise", false)),
project_(dict.getOrDefault<bool>("project", regularise_)),
radius_(regularisationRadius::New(mesh, dict, false)),
alpha_(alpha),
alphaTilda_
(
regularise_
? new volScalarField
(
IOobject
(
"alphaTilda",
mesh_.time().timeName(),
mesh_,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh_,
dimensionedScalar(dimless/dimTime, Zero),
fixedValueFvPatchScalarField::typeName
)
: nullptr
),
sharpenFunction_
(
project_ ?
topOInterpolationFunction::New(mesh, dict) :
nullptr
),
regularisationPDE_(regularisationPDE::New(mesh, dict, zones)),
betaArg_(regularise_ ? alphaTilda_().primitiveField() : alpha_),
growFromWalls_(dict.getOrDefault<bool>("growFromWalls", false)),
beta_
(
IOobject
(
"beta",
mesh_.time().timeName(),
mesh_,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh_,
dimensionedScalar(dimless, Zero),
zeroGradientFvPatchScalarField::typeName
)
{
DebugInfo
<< "Regularise " << Switch(regularise_) << endl;
DebugInfo
<< "Project " << Switch(project_) << endl;
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::fieldRegularisation::updateBeta()
{
if (regularise_)
{
regularise(alpha_, alphaTilda_(), true);
}
if (project_)
{
sharpenFunction_->interpolate(betaArg_, beta_.primitiveFieldRef());
}
else
{
beta_.primitiveFieldRef() = betaArg_;
}
beta_.correctBoundaryConditions();
}
void Foam::fieldRegularisation::regularise
(
const scalarField& source,
scalarField& result,
const bool isTopoField,
const regularisationRadius& radius
)
{
regularisationPDE_->
regularise(alphaTilda_(), source, result, isTopoField, radius);
}
void Foam::fieldRegularisation::regularise
(
const scalarField& source,
scalarField& result,
const bool isTopoField
)
{
regularisationPDE_->
regularise(alphaTilda_(), source, result, isTopoField, radius_());
}
void Foam::fieldRegularisation::postProcessSens(scalarField& sens)
{
// Add dBeta/dBetaArg
if (project_)
{
sens *= sharpenFunction_->derivative(betaArg_);
}
// Add part due to regulatisation
if (regularise_)
{
// Solve the adjoint to the regularisation equation
regularise(sens, sens, false);
}
// Add volume
sens *= mesh_.V();
}
// ************************************************************************* //

197
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/regularisation/fieldRegularisation.H Normal file
View File

@ -0,0 +1,197 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2007-2023 PCOpt/NTUA
Copyright (C) 2013-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
fieldRegularisation
Description
Regularises the given field using a PDE smoother and then sharpens
(projects) the result using a smoothed Heaviside function. Used in topology
optimisation to regularise the designed geometries and obtain 'grid
independent' solutions.
SourceFiles
fieldRegularisation.C
\*---------------------------------------------------------------------------*/
#ifndef fieldRegularisation_H
#define fieldRegularisation_H
#include "fvMesh.H"
#include "volFields.H"
#include "topOZones.H"
#include "regularisationPDE.H"
#include "regularisationRadius.H"
#include "topOInterpolationFunction.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class fieldRegularisation Declaration
\*---------------------------------------------------------------------------*/
class fieldRegularisation
{
protected:
// Protected data
fvMesh& mesh_;
dictionary dict_;
//- Cell zones related to topology optimisation
const topOZones& zones_;
//- Perform regulaisation on alpha before inputing it on beta?
bool regularise_;
//- Perform the projection (sharpening) step?
bool project_;
//- Smoothing radius
// May be isotropic or differ per spatial direction
autoPtr<regularisationRadius> radius_;
//- Alpha field (design variables of topology optimisation)
const scalarField& alpha_;
//- The regularised alpha field, if regulatisation is performed
autoPtr<volScalarField> alphaTilda_;
//- Function used to sharpen the field after regularisation
autoPtr<topOInterpolationFunction> sharpenFunction_;
//- PDE used for the regularisation
autoPtr<regularisationPDE> regularisationPDE_;
//- Argument of the beta field.
// Is either alpha_ (no regularisation)
// or alphaTilda_ (with regularisation)
const scalarField& betaArg_;
//- Whether to apply a fixedValue BC or zeroGradient one to alphaTilda,
//- when regularisation is performed
bool growFromWalls_;
//- Beta is the field used for all interpolations between fluid and
//- solid in topology optimisation
volScalarField beta_;
private:
// Private Member Functions
//- No copy construct
fieldRegularisation(const fieldRegularisation&);
//- No copy assignment
void operator=(const fieldRegularisation&);
public:
//- Runtime type information
TypeName("fieldRegularisation");
// Constructors
//- Construct from components
fieldRegularisation
(
fvMesh& mesh,
const scalarField& alpha,
const topOZones& zones,
const dictionary& dict
);
//- Destructor
virtual ~fieldRegularisation() = default;
// Member Functions
//- Return beta field
inline const volScalarField& beta() const
{
return beta_;
}
//- Grow beta from walls
inline bool growFromWalls() const
{
return growFromWalls_;
}
//- Should regularisation be executed
inline bool shouldRegularise() const
{
return regularise_;
}
//- Update the beta field.
// Performs regulaisation of alpha_, if necessary
virtual void updateBeta();
//- Regularise an externally provided radius
void regularise
(
const scalarField& source,
scalarField& result,
const bool isTopoField,
const regularisationRadius& radius
);
//- Regularise field (or sensitivities) using a Laplace PDE
void regularise
(
const scalarField& source,
scalarField& result,
const bool isTopoField
);
//- Update part of fieldRegularisation to the sensitivitiy derivatives
void postProcessSens(scalarField& sens);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

251
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/regularisation/regularisationPDE/Helmoltz/Helmholtz.C Normal file
View File

@ -0,0 +1,251 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 PCOpt/NTUA
Copyright (C) 2021 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "Helmholtz.H"
#include "fixedValueFvPatchFields.H"
#include "zeroGradientFvPatchFields.H"
#include "wallFvPatch.H"
#include "DynamicList.H"
#include "fvMeshSubset.H"
#include "fvm.H"
#include "bound.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(Helmholtz, 1);
addToRunTimeSelectionTable(regularisationPDE, Helmholtz, dictionary);
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
void Foam::Helmholtz::solveEqn
(
const volScalarField& aTilda,
const scalarField& source,
scalarField& result,
const bool isTopoField,
const regularisationRadius& radius,
const scalar minSetValue,
const bool fixATildaValues
)
{
const fvMesh& mesh = aTilda.internalField().mesh();
// Convergence criteria
const label iters = dict_.getOrDefault<label>("iters", 500);
const scalar tolerance = dict_.getOrDefault<scalar>("tolerance", 1.e-06);
dimensionedScalar one("1", dimless, 1.);
// Smoothed field
volScalarField bTilda
(
IOobject
(
"bTilda",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar(dimless, Zero),
(
growFromWalls_ ?
fixedValueFvPatchScalarField::typeName :
zeroGradientFvPatchScalarField::typeName
)
);
// If solution corresponds to the topology density field, modify boundary
// conditions accordingly
if (isTopoField && growFromWalls_)
{
// Apply a unit alphaTilda value next to all walls
forAll(mesh.boundary(), patchI)
{
const fvPatch& patch = mesh.boundary()[patchI];
if (isA<wallFvPatch>(patch))
{
bTilda.boundaryFieldRef()[patchI] == scalar(1);
}
}
}
// Source field
DimensionedField<scalar, volMesh> sourceField
(
IOobject
(
"source",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimless,
source
);
for (label iter = 0; iter < iters; ++iter)
{
fvScalarMatrix smoothEqn
(
fvm::Sp(one, bTilda)
==
sourceField
);
radius.addRegularisationTerm(smoothEqn, isTopoField);
// Set solution in given zones
if (fixATildaValues)
{
setValues(smoothEqn, isTopoField, minSetValue);
}
const scalar residual(mag(smoothEqn.solve().initialResidual()));
// if (isTopoField)
// {
// bound(bTilda, dimensionedScalar(bTilda.dimensions(), minSetValue));
// }
// Print execution time
mesh.time().printExecutionTime(Info);
// Check convergence
if (residual < tolerance)
{
Info<< "\n***Reached regularisation equation convergence limit, "
"iteration " << iter << "***\n\n";
break;
}
}
// Replace field with its regularised counterpart
result = bTilda.primitiveField();
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::Helmholtz::Helmholtz
(
const fvMesh& mesh,
const dictionary& dict,
const topOZones& zones
)
:
regularisationPDE(mesh, dict, zones),
solveOnActiveCells_(dict.getOrDefault<bool>("solveOnActiveCells", false))
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::Helmholtz::regularise
(
const volScalarField& aTilda,
const scalarField& source,
scalarField& result,
const bool isTopoField,
const regularisationRadius& radius,
const scalar minSetValue,
const bool fixATildaValues
)
{
// Set values for all constant cells here.
// If a subsetMesh is used, cells outside its domain will not be changed,
// potentially leading to fixed cells not getting their correct values
if (fixATildaValues)
{
DynamicList<label> cells(0);
DynamicList<scalar> values(0);
setValues(mesh_, cells, values, isTopoField);
result.rmap(values, cells);
}
/*
// Solve the regularisation equation on a mesh including only the active
// cells, if needed
if (solveOnActiveCells_)
{
const labelList& activeZones = zones_.adjointPorousZoneIDs();
if (!activeZones.empty())
{
Info<< "Solving regularisation equation on active cells only"
<< endl;
DynamicList<label> allActiveCells(0);
for (const label zI : activeZones)
{
allActiveCells.append(mesh_.cellZones()[zI]);
}
fvMeshSubset::exposedPatchType = wallPolyPatch::typeName;
fvMeshSubset subSetMesh(mesh_, allActiveCells);
fvMesh& subMesh = subSetMesh.subMesh();
schemesLookup& fvSchemes = static_cast<schemesLookup&>(subMesh);
fvSchemes.readOpt() = IOobject::MUST_READ_IF_MODIFIED;
fvSchemes.read();
fvSolution& solution = static_cast<fvSolution&>(subMesh);
solution.readOpt(IOobject::MUST_READ_IF_MODIFIED);
solution.read();
const labelList& cellMap = subSetMesh.cellMap();
// Map input fields to subSetMesh fields
volScalarField aTildaSub(subSetMesh.interpolate(aTilda));
scalarField sourceSub(source, cellMap);
scalarField resultSub(result, cellMap);
// Solve the regularisation equation
solveEqn
(
aTildaSub, sourceSub, resultSub,
isTopoField, radius, minSetValue, fixATildaValues
);
// Map result back to original field
result.rmap(resultSub, cellMap);
Info<< "min max " << gMin(result) << " " << gMax(result) << endl;
return;
}
}
*/
solveEqn
(
aTilda,
source,
result,
isTopoField,
radius,
minSetValue,
fixATildaValues
);
}
// ************************************************************************* //

141
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/regularisation/regularisationPDE/Helmoltz/Helmholtz.H Normal file
View File

@ -0,0 +1,141 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 PCOpt/NTUA
Copyright (C) 2021 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::Helmholtz
Description
A regulatisation PDE based on a Helmholtz filter
SourceFiles
Helmholtz.C
\*---------------------------------------------------------------------------*/
#ifndef Helmholtz_H
#define Helmholtz_H
#include "regularisationPDE.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class Helmholtz Declaration
\*---------------------------------------------------------------------------*/
class Helmholtz
:
public regularisationPDE
{
private:
// Private Member Functions
//- No copy construct
Helmholtz(const Helmholtz&) = delete;
//- No copy assignment
void operator=(const Helmholtz&) = delete;
protected:
// Protected Data
//- Solve the regularisationPDE only on a subset mesh made of the
//- active cell zones
bool solveOnActiveCells_;
// Protected Member Functions
//- Solve the regulatisation equation
// The mesh used is the one of aTilda
void solveEqn
(
const volScalarField& aTilda,
const scalarField& source,
scalarField& result,
const bool isTopoField,
const regularisationRadius& radius,
const scalar minSetValue = Zero,
const bool fixATildaValues = true
);
public:
//- Runtime type information
TypeName("Helmholtz");
// Constructors
//- Construct from components
Helmholtz
(
const fvMesh& mesh,
const dictionary& dict,
const topOZones& zones
);
//- Destructor
virtual ~Helmholtz() = default;
// Member Functions
//- Regularise field (or sensitivities) using a Laplace PDE
virtual void regularise
(
const volScalarField& aTilda,
const scalarField& source,
scalarField& result,
const bool isTopoField,
const regularisationRadius& radius,
const scalar minSetValue = Zero,
const bool fixATildaValues = true
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

163
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/regularisation/regularisationPDE/regularisationPDE/regularisationPDE.C Normal file
View File

@ -0,0 +1,163 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 PCOpt/NTUA
Copyright (C) 2021 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "regularisationPDE.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(regularisationPDE, 0);
defineRunTimeSelectionTable(regularisationPDE, dictionary);
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
void Foam::regularisationPDE::setValues
(
fvScalarMatrix& bTildaEqn,
const bool isTopoField,
const scalar minSetValue
)
{
const fvMesh& mesh = bTildaEqn.psi().internalField().mesh();
DynamicList<label> cells(mesh.nCells()/100);
DynamicList<scalar> values(mesh.nCells()/100);
setValues(mesh, cells, values, isTopoField, minSetValue);
bTildaEqn.setValues(cells, values);
}
void Foam::regularisationPDE::setValues
(
const fvMesh& mesh,
DynamicList<label>& cells,
DynamicList<scalar>& values,
bool isTopoField,
const scalar minSetValue
)
{
const labelList& IOcells = mesh.cellZones()[zones_.IOzoneID()];
cells.append(IOcells);
values.append(scalarField(IOcells.size(), minSetValue));
forAll(zones_.fixedPorousZoneIDs(), zI)
{
const label cellZoneID = zones_.fixedPorousZoneIDs()[zI];
const labelList& zoneCells = mesh.cellZones()[cellZoneID];
cells.append(zoneCells);
values.append
(
scalarField
(
zoneCells.size(),
isTopoField ?
zones_.fixedPorousValues()[zI] : scalar(0)
)
);
}
for (label cellZoneID : zones_.fixedZeroPorousZoneIDs())
{
const labelList& zoneCells = mesh.cellZones()[cellZoneID];
cells.append(zoneCells);
values.append(scalarField(zoneCells.size(), minSetValue));
}
}
Foam::scalar Foam::regularisationPDE::computeRadius()
{
scalar averageVol(gAverage(mesh_.V().field()));
const Vector<label>& geometricD = mesh_.geometricD();
const boundBox& bounds = mesh_.bounds();
forAll(geometricD, iDir)
{
if (geometricD[iDir] == -1)
{
averageVol /= bounds.span()[iDir];
}
}
scalar radius = pow(averageVol, scalar(1)/scalar(mesh_.nGeometricD()));
Info<< "Computed a mean radius of " << radius << endl;
return radius;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::regularisationPDE::regularisationPDE
(
const fvMesh& mesh,
const dictionary& dict,
const topOZones& zones
)
:
mesh_(mesh),
dict_(dict),
zones_(zones),
growFromWalls_(dict.getOrDefault<bool>("growFromWalls", false))
{}
// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
Foam::autoPtr<Foam::regularisationPDE> Foam::regularisationPDE::New
(
const fvMesh& mesh,
const dictionary& dict,
const topOZones& zones
)
{
const word modelType =
dict.getOrDefault<word>("regularisationPDE", "Helmholtz");
auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
Info<< "regularisationPDE type " << modelType << endl;
if (!cstrIter.found())
{
FatalIOErrorInLookup
(
dict,
"regularisationPDE",
modelType,
*dictionaryConstructorTablePtr_
) << exit(FatalIOError);
}
return autoPtr<regularisationPDE>
(
cstrIter()(mesh, dict, zones)
);
}
// ************************************************************************* //

184
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/regularisation/regularisationPDE/regularisationPDE/regularisationPDE.H Normal file
View File

@ -0,0 +1,184 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 PCOpt/NTUA
Copyright (C) 2021 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::regularisationPDE
Description
Base class for selecting the regulatisation PDE
SourceFiles
regularisationPDE.C
\*---------------------------------------------------------------------------*/
#ifndef regularisationPDE_H
#define regularisationPDE_H
#include "fvMesh.H"
#include "topOZones.H"
#include "fvMatrices.H"
#include "regularisationRadius.H"
#include "runTimeSelectionTables.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class regularisationPDE Declaration
\*---------------------------------------------------------------------------*/
class regularisationPDE
{
private:
// Private Member Functions
//- No copy construct
regularisationPDE(const regularisationPDE&) = delete;
//- No copy assignment
void operator=(const regularisationPDE&) = delete;
protected:
// Protected data
const fvMesh& mesh_;
const dictionary dict_;
//- Cell zones related to topology optimisation
const topOZones& zones_;
//- Whether to apply a fixedValue BC or zeroGradient one to alphaTilda,
//- when regularisation is performed
bool growFromWalls_;
// Protected Member Functions
//- Set fixed bTilda values
void setValues
(
fvScalarMatrix& bTildaEqn,
const bool isTopoField,
const scalar minSetValue = Zero
);
//- Gather the cells with constant values from all constrained zones
void setValues
(
const fvMesh& mesh,
DynamicList<label>& cells,
DynamicList<scalar>& values,
bool isTopoField,
const scalar minSetValue = Zero
);
//- Compute smoothing radius, if not directly given
scalar computeRadius();
public:
//- Runtime type information
TypeName("regularisationPDE");
// Declare run-time constructor selection table
declareRunTimeSelectionTable
(
autoPtr,
regularisationPDE,
dictionary,
(
const fvMesh& mesh,
const dictionary& dict,
const topOZones& zones
),
(mesh, dict, zones)
);
// Constructors
//- Construct from components
regularisationPDE
(
const fvMesh& mesh,
const dictionary& dict,
const topOZones& zones
);
// Selectors
//- Construct and return the selected regularisationPDE
static autoPtr<regularisationPDE> New
(
const fvMesh& mesh,
const dictionary& dict,
const topOZones& zones
);
//- Destructor
virtual ~regularisationPDE() = default;
// Member Functions
//- Regularise field (or sensitivities)
virtual void regularise
(
const volScalarField& aTilda,
const scalarField& source,
scalarField& result,
const bool isTopoField,
const regularisationRadius& radius,
const scalar minSetValue = Zero,
const bool fixATildaValues = true
) = 0;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

112
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/regularisation/regularisationRadius/isotropic/regularisationRadiusIsotropic.C Normal file
View File

@ -0,0 +1,112 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 PCOpt/NTUA
Copyright (C) 2021 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "regularisationRadiusIsotropic.H"
#include "fvm.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(isotropic, 1);
addToRunTimeSelectionTable(regularisationRadius, isotropic, dictionary);
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
Foam::scalar Foam::isotropic::computeRadius(const dictionary& dict)
{
scalar averageVol(gAverage(mesh_.V().field()));
const Vector<label>& geometricD = mesh_.geometricD();
const boundBox& bounds = mesh_.bounds();
forAll(geometricD, iDir)
{
if (geometricD[iDir] == -1)
{
averageVol /= bounds.span()[iDir];
}
}
scalar radius = pow(averageVol, scalar(1)/scalar(mesh_.nGeometricD()));
scalar multMeanRadius = dict.getOrDefault<scalar>("meanRadiusMult", 10);
Info<< "Computed a mean radius of " << radius
<< " and multiplying with " << multMeanRadius << endl;
return multMeanRadius*radius;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::isotropic::isotropic
(
const fvMesh& mesh,
const dictionary& dict,
bool adjustWallThickness
)
:
regularisationRadius(mesh, dict, adjustWallThickness),
radius_
(
"radius",
dimLength,
scalar
(
dict_.getOrDefault<scalar>("radius", computeRadius(dict))
/(2.*::sqrt(3.))
)
)
{
if (adjustWallThickness)
{
const scalar wallThicknessMult =
dict.getOrDefault<scalar>("wallThicknessMult", 0.75);
DebugInfo<<
"Adjusting wall thickness by " << wallThicknessMult << endl;
radius_ *= wallThicknessMult;
}
DebugInfo
<< "Using radius " << radius_ << endl;
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::isotropic::addRegularisationTerm
(
fvScalarMatrix& matrix,
bool isTopoField
) const
{
const volScalarField& field = matrix.psi();
matrix -= fvm::laplacian(sqr(radius_), field);
}
// ************************************************************************* //

128
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/regularisation/regularisationRadius/isotropic/regularisationRadiusIsotropic.H Normal file
View File

@ -0,0 +1,128 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 PCOpt/NTUA
Copyright (C) 2021 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::isotropic
Description
An isotropic regularisationRadius (same in all spatial directions)
SourceFiles
isotropic.C
\*---------------------------------------------------------------------------*/
#ifndef isotropic_H
#define isotropic_H
#include "regularisationRadius.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class isotropic Declaration
\*---------------------------------------------------------------------------*/
class isotropic
:
public regularisationRadius
{
private:
// Private Member Functions
//- No copy construct
isotropic(const isotropic&) = delete;
//- No copy assignment
void operator=(const isotropic&) = delete;
protected:
// Protected data
//- Multiplier of the wall thickness, used to obtain the radius for
//- the second regularisation in bi-fluid topO
scalar wallThicknessMult_;
//- Smoothing radius of the first regulatisation
dimensionedScalar radius_;
// Protected Member Functions
//- Compute smoothing radius, if not directly given
scalar computeRadius(const dictionary& dict);
public:
//- Runtime type information
TypeName("isotropic");
// Constructors
//- Construct from components
isotropic
(
const fvMesh& mesh,
const dictionary& dict,
bool adjustWallThickness
);
//- Destructor
virtual ~isotropic() = default;
// Member Functions
//- Add a Laplacian term with an isotropic diffusivity
virtual void addRegularisationTerm
(
fvScalarMatrix& matrix,
bool isTopoField
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

87
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/regularisation/regularisationRadius/regularisationRadius/regularisationRadius.C Normal file
View File

@ -0,0 +1,87 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 PCOpt/NTUA
Copyright (C) 2021 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "regularisationRadius.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(regularisationRadius, 0);
defineRunTimeSelectionTable(regularisationRadius, dictionary);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::regularisationRadius::regularisationRadius
(
const fvMesh& mesh,
const dictionary& dict,
bool adjustWallThickness
)
:
mesh_(mesh),
dict_(dict)
{}
// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
Foam::autoPtr<Foam::regularisationRadius> Foam::regularisationRadius::New
(
const fvMesh& mesh,
const dictionary& dict,
bool adjustWallThickness
)
{
const word modelType = dict.getOrDefault<word>("type", "isotropic");
auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
Info<< "regularisationRadius type " << modelType << endl;
if (!cstrIter.found())
{
FatalIOErrorInLookup
(
dict,
"regularisationRadius",
modelType,
*dictionaryConstructorTablePtr_
) << exit(FatalIOError);
}
return autoPtr<regularisationRadius>
(
cstrIter()(mesh, dict, adjustWallThickness)
);
}
// ************************************************************************* //

146
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/regularisation/regularisationRadius/regularisationRadius/regularisationRadius.H Normal file
View File

@ -0,0 +1,146 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 PCOpt/NTUA
Copyright (C) 2021 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::regularisationRadius
Description
Base class for selecting the regulatisation radius
SourceFiles
regularisationRadius.C
\*---------------------------------------------------------------------------*/
#ifndef regularisationRadius_H
#define regularisationRadius_H
#include "fvMatrices.H"
#include "runTimeSelectionTables.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class regularisationRadius Declaration
\*---------------------------------------------------------------------------*/
class regularisationRadius
{
private:
// Private Member Functions
//- No copy construct
regularisationRadius(const regularisationRadius&) = delete;
//- No copy assignment
void operator=(const regularisationRadius&) = delete;
protected:
// Protected data
const fvMesh& mesh_;
const dictionary dict_;
public:
//- Runtime type information
TypeName("regularisationRadius");
// Declare run-time constructor selection table
declareRunTimeSelectionTable
(
autoPtr,
regularisationRadius,
dictionary,
(
const fvMesh& mesh,
const dictionary& dict,
bool adjustWallThickness
),
(mesh, dict, adjustWallThickness)
);
// Constructors
//- Construct from components
regularisationRadius
(
const fvMesh& mesh,
const dictionary& dict,
bool adjustWallThickness
);
// Selectors
//- Construct and return the selected regularisationRadius
// Default to isotropic
static autoPtr<regularisationRadius> New
(
const fvMesh& mesh,
const dictionary& dict,
bool adjustWallThickness
);
//- Destructor
virtual ~regularisationRadius() = default;
// Member Functions
//- Term including the regulatisation of the field
virtual void addRegularisationTerm
(
fvScalarMatrix& matrix,
bool isTopoField
) const = 0;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

574
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/topODesignVariables.C Normal file
View File

@ -0,0 +1,574 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2007-2023 PCOpt/NTUA
Copyright (C) 2013-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "localIOdictionary.H"
#include "topODesignVariables.H"
#include "MeshObject.H"
#include "wallFvPatch.H"
#include "cutFaceIso.H"
#include "cutCellIso.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(topODesignVariables, 0);
addToRunTimeSelectionTable
(
designVariables,
topODesignVariables,
designVariables
);
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
void Foam::topODesignVariables::updateField
(
const scalarField& correction,
const label fluidID
)
{
DebugInfo
<< "Updating design variables for field " << fluidID << endl;
const label n = mesh_.nCells();
SubField<scalar> localCorrection(correction, n, fluidID*n);
SubField<scalar> field(*this, n, fluidID*n);
// Update porosity in adjoint porous cells
if (zones_.adjointPorousZoneIDs().empty())
{
forAll(field, cellI)
{
field[cellI] +=
min
(
max
(
field[cellI] + localCorrection[cellI],
scalar(0)
),
1.
)
- field[cellI];
}
}
else
{
for (label cellZoneID : zones_.adjointPorousZoneIDs())
{
const labelList& zoneCells = mesh_.cellZones()[cellZoneID];
for (label cellI : zoneCells)
{
field[cellI] +=
min
(
max
(
field[cellI] + localCorrection[cellI],
scalar(0)
),
1.
)
- field[cellI];
}
}
}
}
void Foam::topODesignVariables::applyFixedValues()
{
SubField<scalar> alpha(*this, mesh_.nCells());
// Zero porosity in the cells next to IO
for (label cellI : zones_.IOCells())
{
alpha[cellI] = 0.;
}
// Apply fixed porosity
forAll(zones_.fixedPorousZoneIDs(), zI)
{
const label cellZoneID = zones_.fixedPorousZoneIDs()[zI];
const labelList& zoneCells = mesh_.cellZones()[cellZoneID];
const scalar alphaValue(zones_.fixedPorousValues()[zI]);
for (label cellI : zoneCells)
{
alpha[cellI] = alphaValue;
}
}
// Apply fixed zero porosity
for (label cellZoneID : zones_.fixedZeroPorousZoneIDs())
{
const labelList& zoneCells = mesh_.cellZones()[cellZoneID];
for (label cellI : zoneCells)
{
alpha[cellI] = 0.;
}
}
}
Foam::scalar Foam::topODesignVariables::computeEta(scalarField& correction)
{
// Since the correction might want to increase values above 1 or
// reduce them below 0, using it directly to compute eta is a bad idea.
// Instead, do a pseudo-update and compute the max alpha change directly
// Back-up old design variables (dvs)
scalarField& dvs = getVars();
scalarField oldDVs(dvs);
// Do update to compute alpha field with eta 1
update(correction);
// Compute difference and through that, the max dv change
const scalar maxChange(gMax(mag(dvs - oldDVs)));
// Restore alpha back
dvs = oldDVs;
// Compute eta
Info<< "maxInitChange/maxChange \t"
<< maxInitChange_() << "/" << maxChange << endl;
const scalar eta(maxInitChange_() / maxChange);
Info<< "Setting eta value to " << eta << endl;
correction *= eta;
return eta;
}
void Foam::topODesignVariables::setActiveDesignVariables
(
const label fluidID,
const bool activeIO
)
{
const label offset(fluidID*mesh_.nCells());
label varI(activeDesignVariables_.size());
activeDesignVariables_.setSize(offset + mesh_.nCells(), -1);
// Set active design variables
// If specific porosity zones are perscribed, use them directly
if (!zones_.adjointPorousZoneIDs().empty())
{
for (label cellZoneID : zones_.adjointPorousZoneIDs())
{
for (const label var : mesh_.cellZones()[cellZoneID])
{
activeDesignVariables_[varI++] = var + offset;
}
}
}
// Else, pick up all cells in non-constant porosity zones
else
{
boolList isActiveDV(mesh_.nCells(), true);
// Exclude cells with fixed porosity
for (label cellZoneID : zones_.fixedPorousZoneIDs())
{
for (label cellI : mesh_.cellZones()[cellZoneID])
{
isActiveDV[cellI] = false;
}
}
for (label cellZoneID : zones_.fixedZeroPorousZoneIDs())
{
for (label cellI : mesh_.cellZones()[cellZoneID])
{
isActiveDV[cellI] = false;
}
}
if (!activeIO)
{
for (label cellI : zones_.IOCells())
{
isActiveDV[cellI] = false;
}
}
// Set active design variables
forAll(isActiveDV, cellI)
{
if (isActiveDV[cellI])
{
activeDesignVariables_[varI++] = offset + cellI;
}
}
}
activeDesignVariables_.setSize(varI);
}
void Foam::topODesignVariables::readField
(
const word& name,
const label fluidID,
const bool setIOValues
)
{
const label offset(fluidID*mesh_.nCells());
if (localIOdictionary::found(name))
{
SubField<scalar>(*this, mesh_.nCells(), offset) =
scalarField(name, *this, mesh_.nCells());
/*
// Set values next to IO boundaries if needed
if (setIOValues)
{
forAll(mesh_.boundary(), patchI)
{
const fvPatch& patch = mesh_.boundary()[patchI];
if (patch.type() == "patch")
{
const labelList& faceCells = patch.faceCells();
const scalarField& pf = volField.boundaryField()[patchI];
forAll(faceCells, fI)
{
scalarField::operator[](offset + faceCells[fI]) = pf[fI];
}
}
}
}
*/
}
}
void Foam::topODesignVariables::initialize()
{
// Set active design variables
setActiveDesignVariables();
// Read in values from file, if present
readField("alpha", 0, true);
if (regularisation_.growFromWalls())
{
scalarField& alpha = *this;
for (const fvPatch& patch : mesh_.boundary())
{
if (isA<wallFvPatch>(patch))
{
const labelList& faceCells = patch.faceCells();
forAll(faceCells, cI)
{
alpha[faceCells[cI]] = 1.;
}
}
}
}
// Make sure alpha has fixed values where it should
scalarField dummyCorrection(mesh_.nCells(), Zero);
update(dummyCorrection);
evolveNumber();
// Read bounds for design variables, if present
readBounds(autoPtr<scalar>::New(0), autoPtr<scalar>::New(1));
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::topODesignVariables::topODesignVariables
(
fvMesh& mesh,
const dictionary& dict
)
:
topODesignVariables(mesh, dict, mesh.nCells())
{}
Foam::topODesignVariables::topODesignVariables
(
fvMesh& mesh,
const dictionary& dict,
const label size
)
:
topOVariablesBase(mesh, dict),
designVariables(mesh, dict, size),
alpha_(SubField<scalar>(*this, mesh.nCells(), 0)),
regularisation_
(
mesh,
alpha_,
zones_,
dict_.subDict("regularisation")
),
writeAllFields_
(
dict.getOrDefaultCompat<bool>
(
"writeAllFields", {{"writeAllAlphaFields", 2306}}, false
)
),
addFvOptions_(dict.getOrDefault<bool>("addFvOptions", false))
{
// Rest of the contrsuctor initialisation
initialize();
}
// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
Foam::autoPtr<Foam::topODesignVariables> Foam::topODesignVariables::New
(
fvMesh& mesh,
const dictionary& dict
)
{
return autoPtr<topODesignVariables>
(
new topODesignVariables(mesh, dict)
);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * //
const Foam::volScalarField& Foam::topODesignVariables::beta() const
{
return regularisation_.beta();
}
const Foam::scalarField& Foam::topODesignVariables::interpolationField
(
const word& interpolationFieldName
) const
{
return beta().primitiveField();
}
void Foam::topODesignVariables::interpolate
(
volScalarField& field,
const topOInterpolationFunction& interpolationFunc,
const FieldField<Field, scalar>& fluidValues,
const scalarField& solidValues,
const label fieldi,
const word& interpolationFieldName
) const
{
const scalarField& indicator = interpolationField(interpolationFieldName);
scalarField interpolant(indicator.size(), Zero);
interpolationFunc.interpolate(indicator, interpolant);
// Interpolate field values
const scalar diff(solidValues[fieldi] - fluidValues[0][fieldi]);
field.primitiveFieldRef() = fluidValues[0][fieldi] + interpolant*diff;
field.correctBoundaryConditions();
}
void Foam::topODesignVariables::interpolationSensitivities
(
scalarField& sens,
const topOInterpolationFunction& interpolationFunc,
const FieldField<Field, scalar>& fluidValues,
const scalarField& solidValues,
const label fieldi,
const word& designVariablesName,
const word& interpolationFieldName
) const
{
const scalarField& indicator = interpolationField(interpolationFieldName);
sens *=
(solidValues[fieldi] - fluidValues[0][fieldi])
*interpolationFunc.derivative(indicator);
}
void Foam::topODesignVariables::nullifyInSolid
(
scalarField& field,
const topOInterpolationFunction& interpolationFunc
) const
{
const scalarField& beta = this->beta().primitiveField();
scalarField interpolant(beta.size(), Zero);
interpolationFunc.interpolate(beta, interpolant);
field *= scalar(1) - interpolant;
}
void Foam::topODesignVariables::nullifyInSolidSensitivities
(
scalarField& sens,
const topOInterpolationFunction& interpolationFunc,
const word& designVariablesName
) const
{
const scalarField& beta = this->beta().primitiveField();
sens *= - interpolationFunc.derivative(beta);
}
Foam::tmp<Foam::scalarField> Foam::topODesignVariables::penalty
(
const word& interpolationFieldName,
const topOInterpolationFunction& interpolationFunc
) const
{
const scalarField& beta = this->beta().primitiveField();
tmp<scalarField> tinterpolant(tmp<scalarField>::New(beta.size(), Zero));
interpolationFunc.interpolate(beta, tinterpolant.ref());
return tinterpolant;
}
Foam::tmp<Foam::scalarField> Foam::topODesignVariables::penaltySensitivities
(
const word& interpolationFieldName,
const topOInterpolationFunction& interpolationFunc
) const
{
const scalarField& beta = this->beta().primitiveField();
return interpolationFunc.derivative(beta);
}
void Foam::topODesignVariables::update(scalarField& correction)
{
// Update alpha values
updateField(correction);
// Fix alpha in zones
applyFixedValues();
// Update the beta field
regularisation_.updateBeta();
// Though the mesh is kept constant, the distance from wall may change
// if the method computing it includes fvOptions that depend on the
// indicator field.
// Trick wallDist into updating it
if (mesh_.foundObject<UpdateableMeshObject<fvMesh>>("wallDist"))
{
mesh_.lookupObjectRef<UpdateableMeshObject<fvMesh>>("wallDist").
movePoints();
}
// Write the 0.5 beta iso-line to files, as an indication of the
// fluid-solid interface
labelList changedFaces(mesh_.nFaces(), -1);
List<wallPoint> changedFacesInfo(mesh_.nFaces());
writeFluidSolidInterface(-beta(), -0.5, changedFaces, changedFacesInfo);
}
bool Foam::topODesignVariables::globalSum() const
{
return true;
}
Foam::tmp<Foam::scalarField> Foam::topODesignVariables::assembleSensitivities
(
adjointSensitivity& adjointSens
)
{
// Raw sensitivities field.
// Does not include the adjoint to the regularisation and projection steps
const scalarField& fieldSens = adjointSens.fieldSensPtr()->primitiveField();
// Return field (complete sensitivities)
auto tobjectiveSens(tmp<scalarField>::New(fieldSens));
scalarField& objectiveSens = tobjectiveSens.ref();
// Add part due to regularisation and projection
regularisation_.postProcessSens(objectiveSens);
return tobjectiveSens;
}
void Foam::topODesignVariables::addFvOptions
(
const PtrList<primalSolver>& primalSolvers,
const PtrList<adjointSolverManager>& adjointSolverManagers
)
{
// WIP
if (addFvOptions_)
{
for (const primalSolver& solver : primalSolvers)
{
solver.addTopOFvOptions();
}
for (const adjointSolverManager& manager : adjointSolverManagers)
{
const PtrList<adjointSolver>& adjointSolvers =
manager.adjointSolvers();
for (const adjointSolver& solver : adjointSolvers)
{
solver.addTopOFvOptions();
}
}
}
}
void Foam::topODesignVariables::writeDesignVars()
{
if (writeAllFields_ && mesh_.time().writeTime())
{
volScalarField alpha
(
IOobject
(
"alpha",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh_,
dimensionedScalar(dimless/dimTime, Zero)
);
alpha.primitiveFieldRef() = alpha_;
alpha.write();
}
}
bool Foam::topODesignVariables::writeData(Ostream& os) const
{
alpha_.writeEntry("alpha", os);
return true;
}
// ************************************************************************* //

278
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/topODesignVariables.H Normal file
View File

@ -0,0 +1,278 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2007-2023 PCOpt/NTUA
Copyright (C) 2013-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::topODesignVariables
Description
Design variables for density-based topology optimisation (topO) problems
SourceFiles
topODesignVariables.C
\*---------------------------------------------------------------------------*/
#ifndef topODesignVariables_H
#define topODesignVariables_H
#include "designVariables.H"
#include "topOVariablesBase.H"
#include "fieldRegularisation.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class topODesignVariables Declaration
\*---------------------------------------------------------------------------*/
class topODesignVariables
:
public topOVariablesBase,
public designVariables
{
protected:
// Protected data
//- A subfield of the design variables correponding to the porosity
//- field.
// Usually the same as *this
SubField<scalar> alpha_;
//- Mechanism to regularise the field of design variables
fieldRegularisation regularisation_;
//- Write all fields related to imposition of the Brinkman
//- penalisation (i.e. design variables, regularised design variables
//- and the indicator field)
bool writeAllFields_;
//- Add the fvOptions necessary for topO automatically
// WIP
bool addFvOptions_;
// Protected Member Functions
//- Update the design variables given their correction
virtual void updateField
(
const scalarField& correction,
const label fluidID = 0
);
//- Apply fixed values in certain zones
void applyFixedValues();
//- Set active design variables
virtual void setActiveDesignVariables
(
const label fluidID = 0,
const bool activeIO = false
);
//- Read field with (path of) the design variables and store input
//- in the design variables list with optional offset
void readField
(
const word& name,
const label fluidID = 0,
const bool setIOValues = false
);
//- Part of the constructor initialisation
virtual void initialize();
private:
// Private Member Functions
//- Disallow default bitwise copy construct
topODesignVariables(const topODesignVariables&) = delete;
//- Disallow default bitwise assignment
void operator=(const topODesignVariables&) = delete;
public:
//- Runtime type information
TypeName("topO");
// Constructors
//- Construct from dictionary
topODesignVariables
(
fvMesh& mesh,
const dictionary& dict
);
//- Construct from dictionary and size
topODesignVariables
(
fvMesh& mesh,
const dictionary& dict,
const label size
);
// Selectors
//- Construct and return the selected design variables
static autoPtr<topODesignVariables> New
(
fvMesh& mesh,
const dictionary& dict
);
//- Destructor
virtual ~topODesignVariables() = default;
// Member Functions
//- Get the indicator field
virtual const volScalarField& beta() const;
//- Return interpolant.
// Defaults to beta in mono-fluid topO optimization
virtual const scalarField& interpolationField
(
const word& interpolationFieldName = "beta"
) const;
//- Interpolate between the given field and solid values
virtual void interpolate
(
volScalarField& field,
const topOInterpolationFunction& interpolationFunc,
const FieldField<Field, scalar>& fluidValues,
const scalarField& solidValues,
const label fieldi,
const word& interpolationFieldName = "beta"
) const;
//- Post-processing sensitivities due to interpolations based on
//- the indicator fields
virtual void interpolationSensitivities
(
scalarField& sens,
const topOInterpolationFunction& interpolationFunc,
const FieldField<Field, scalar>& fluidValues,
const scalarField& solidValues,
const label fieldi,
const word& designVariablesName,
const word& interpolationFieldName = "beta"
) const;
//- Nullify given field in the solid domain
virtual void nullifyInSolid
(
scalarField& field,
const topOInterpolationFunction& interpolationFunc
) const;
//- Nullify given field in the solid domain
virtual void nullifyInSolidSensitivities
(
scalarField& sens,
const topOInterpolationFunction& interpolationFunc,
const word& designVariablesName
) const;
//- Return the Brinkman penalisation term
virtual tmp<scalarField> penalty
(
const word& interpolationFieldName,
const topOInterpolationFunction& interpolationFunc
) const;
//- Return the penalty term derivative
virtual tmp<scalarField> penaltySensitivities
(
const word& interpolationFieldName,
const topOInterpolationFunction& interpolationFunc
) const;
//- Update design variables based on a given correction
virtual void update(scalarField& correction);
//- Compute eta if not set in the first step
virtual scalar computeEta(scalarField& correction);
//- Whether to use global sum when computing matrix-vector products
//- in update methods
virtual bool globalSum() const;
//- Assemble sensitivity derivatives, by combining the part related
//- to the primal and adjoint solution with the part related to the
//- design variables.
// Takes into consideration the adjoint to the regularisation process
// and the derivative of the interpolationFunc
virtual tmp<scalarField> assembleSensitivities
(
adjointSensitivity& sens
);
//- Automatically add fvOptions depending on the design variables
//- to the primal and adjoint solvers
// WIP
virtual void addFvOptions
(
const PtrList<primalSolver>& primalSolver,
const PtrList<adjointSolverManager>& adjointSolverManagers
);
//- Get access to the regularisation object
inline fieldRegularisation& getRegularisation()
{
return regularisation_;
}
//- Write porosity field to file
virtual void writeDesignVars();
//- The writeData function required by the regIOobject write operation
virtual bool writeData(Ostream&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

510
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/topOVariablesBase/topOVariablesBase.C Normal file
View File

@ -0,0 +1,510 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2007-2023 PCOpt/NTUA
Copyright (C) 2013-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "topOVariablesBase.H"
#include "pointMesh.H"
#include "topOZones.H"
#include "volFields.H"
#include "volPointInterpolation.H"
#include "wallPolyPatch.H"
#include "coupledFvPatch.H"
#include "emptyFvPatch.H"
#include "MeshedSurfaceProxy.H"
#include "surfaceWriter.H"
#include "foamVtkSurfaceWriter.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(topOVariablesBase, 0);
}
// * * * * * * * * * * * Protected Members Functions * * * * * * * * * * * * //
bool Foam::topOVariablesBase::addCutBoundaryFaceToIsoline
(
const label facei,
const cutFaceIso& cutFace,
DynamicList<vector>& isoSurfPts,
DynamicList<face>& isoSurfFaces,
DynamicList<label>& zoneIDs,
label& nIsoSurfPts
) const
{
const fvMesh& mesh = zones_.mesh();
// If this is a boundary face being cut, append it to the iso-surface
if (!mesh.isInternalFace(facei))
{
const label cutPatchi = mesh.boundaryMesh().whichPatch(facei);
const fvPatch& cutPatch = mesh.boundary()[cutPatchi];
if
(
!isA<coupledFvPatch>(cutPatch)
&& !isA<emptyFvPatch>(cutPatch)
&& (cutFace.subFacePoints().size() > 2)
)
{
const DynamicList<point>& facePts = cutFace.subFacePoints();
face isoFace(identity(facePts.size(), nIsoSurfPts));
isoSurfPts.append(facePts);
isoSurfFaces.append(isoFace);
zoneIDs.append(cutPatchi);
nIsoSurfPts += facePts.size();
return true;
}
}
return false;
}
bool Foam::topOVariablesBase::addCuttingFaceToIsoline
(
const DynamicList<point>& facePoints,
DynamicList<vector>& isoSurfPts,
DynamicList<face>& isoSurfFaces,
DynamicList<label>& zoneIDs,
label& nIsoSurfPts,
const label nSerialPatches
) const
{
if (facePoints.size() > 2)
{
face isoFace(identity(facePoints.size(), nIsoSurfPts));
isoSurfPts.append(facePoints);
isoSurfFaces.append(isoFace);
zoneIDs.append(nSerialPatches);
nIsoSurfPts += facePoints.size();
return true;
}
return false;
}
void Foam::topOVariablesBase::addBoundaryFacesToIsoline
(
const pointScalarField& pointY,
const Map<label>& addedFaces,
const scalar isoValue,
DynamicList<vector>& isoSurfPts,
DynamicList<face>& isoSurfFaces,
DynamicList<label>& zoneIDs,
label& nIsoSurfPts
) const
{
const fvMesh& mesh = zones_.mesh();
const pointField& points = mesh.points();
const faceList& faces = mesh.faces();
forAll(mesh.boundary(), patchi)
{
const fvPatch& patch = mesh.boundary()[patchi];
if
(
!isA<emptyFvPatch>(patch)
&& !isA<coupledFvPatch>(patch)
)
{
const label start = patch.start();
forAll(patch, facei)
{
const label gFacei = start + facei;
// On rare occasions, a face value might be negative but all
// point values are positive. In such a case, the face won't be
// cut and is seen as a fluid face. Hence, the point values are
// used to determine the status of the boundary face.
bool isFluid(true);
for (const label pti : mesh.faces()[gFacei])
{
isFluid = isFluid && (pointY[pti] >= isoValue);
}
if (isFluid && !addedFaces.found(gFacei))
{
const pointField facePoints = faces[gFacei].points(points);
const face isoFace
(identity(facePoints.size(), nIsoSurfPts));
isoSurfPts.append(facePoints);
isoSurfFaces.append(isoFace);
zoneIDs.append(patchi);
nIsoSurfPts += facePoints.size();
}
}
}
}
}
void Foam::topOVariablesBase::writeSurfaceFiles
(
DynamicList<vector>& patchPoints,
DynamicList<face>& patchFaces,
DynamicList<label>& zoneIDs,
const label nSerialPatches
) const
{
const fvMesh& mesh = zones_.mesh();
pointField pts(std::move(patchPoints));
faceList faces(std::move(patchFaces));
labelList zoneIds(std::move(zoneIDs));
// Write vtp file
const word timeName = mesh.time().timeName();
fileName localName("topOIsoSurface" + timeName);
fileName fname(isoSurfFolder_/localName);
vtk::surfaceWriter vtkFile(pts, faces, fname);
vtkFile.beginFile();
vtkFile.writeGeometry();
vtkFile.beginCellData(1);
vtkFile.writeCellData("zoneIds", zoneIds);
vtkFile.close();
// It appears that there is no interface to write a multi-zone stl file (?)
// Follow a process similar to proxySurfaceWriter, but also use the zoneIds.
// Dummy faceIds.
// faceMap is built after merging the geometry from all processors, to
// be based on the global addressing
labelList faceIds;
autoPtr<meshedSurf> surf(nullptr);
if (Pstream::parRun())
{
surf.reset
(
autoPtr<meshedSurf>::NewFrom<mergedSurf>
(
pts, faces, zoneIds, faceIds, surfaceWriter::defaultMergeDim
)
);
}
else
{
surf.reset
(
autoPtr<meshedSurf>::NewFrom<meshedSurfRef>
(
pts, faces, zoneIds, faceIds
)
);
}
if (Pstream::master())
{
const faceList& surfFaces = surf().faces();
const labelList& surfZoneIds = surf().zoneIds();
// Size per zone
labelList zoneSizes(nSerialPatches + 1, 0);
forAll(surfFaces, fi)
{
++zoneSizes[surfZoneIds[fi]];
}
// Faces passed in previous zones
labelList cumulZoneSizes(surfZoneIds.size(), 0);
for (label zi = 1; zi < surfZoneIds.size() - 1; ++zi)
{
cumulZoneSizes[zi] = cumulZoneSizes[zi - 1] + zoneSizes[zi - 1];
}
// Construction of the faceMap
// ---------------------------
labelList faceMap(surfFaces.size(), -1);
// Faces visited so far per zone
labelList passedFacesPerZone(surfZoneIds.size(), 0);
forAll(surfFaces, fi)
{
const label zi = surfZoneIds[fi];
faceMap[cumulZoneSizes[zi] + passedFacesPerZone[zi]] = fi;
++passedFacesPerZone[zi];
}
// Construction of the surfZoneList
// --------------------------------
List<surfZone> surfZones(nSerialPatches + 1);
label zi = 0;
forAll(mesh.boundary(), pi)
{
const fvPatch& patch = mesh.boundary()[pi];
const word& name = patch.name();
if (!isA<coupledFvPatch>(patch) && !isA<emptyFvPatch>(patch))
{
surfZones[zi] =
surfZone(name, zoneSizes[zi], cumulZoneSizes[zi], zi);
zi++;
}
}
surfZones[zi] =
surfZone
(
"topOPatch",
zoneSizes[nSerialPatches],
cumulZoneSizes[nSerialPatches],
zi
);
surfZones.setSize(zi + 1);
MeshedSurfaceProxy<face>
(
surf().points(),
surfFaces,
surfZones,
faceMap
).write
(
fname + ".stl"
// BINARY seems to loose the patch names
//IOstreamOption
//(
// IOstreamOption::streamFormat::ASCII,
// IOstreamOption::compressionType::COMPRESSED
//)
);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::topOVariablesBase::topOVariablesBase
(
fvMesh& mesh,
const dictionary& dict
)
:
localIOdictionary
(
IOobject
(
"topoVars",
mesh.time().timeName(),
fileName("uniform"),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
word::null
),
zones_(mesh, dict),
isoSurfFolder_
(mesh.time().globalPath()/"optimisation"/"topOIsoSurfaces")
{
mkDir(isoSurfFolder_);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * //
void Foam::topOVariablesBase::sourceTerm
(
DimensionedField<scalar, volMesh>& field,
const topOInterpolationFunction& interpolationFunc,
const scalar betaMax,
const word& interpolationFieldName
) const
{
// Interpolate based on beta values
const scalarField& beta = this->beta().primitiveField();
interpolationFunc.interpolate(beta, field.field());
field *= betaMax;
}
void Foam::topOVariablesBase::sourceTermSensitivities
(
scalarField& sens,
const topOInterpolationFunction& interpolationFunc,
const scalar betaMax,
const word& designVariablesName,
const word& interpolationFieldName
) const
{
if (designVariablesName == "beta")
{
// Multiply with derivative of the interpolation function
const scalarField& beta = this->beta().primitiveField();
sens *= betaMax*interpolationFunc.derivative(beta);
}
}
void Foam::topOVariablesBase::writeFluidSolidInterface
(
const volScalarField& indicator,
const scalar isoValue,
labelList& changedFaces,
List<wallPoint>& changedFacesInfo
)
{
const fvMesh& mesh = zones_.mesh();
// Current number of interface faces
label nSurfFaces(0);
// Indicator at the mesh points, used to compute the iso-surface
pointScalarField pointY(volPointInterpolation(mesh).interpolate(indicator));
if (debug && mesh.time().writeTime())
{
pointY.rename("pointY");
pointY.write();
}
// Storage for the cutting face/cell mechanism
cutCellIso cutCell(mesh, pointY);
cutFaceIso cutFace(mesh, pointY);
DynamicList<face> isoSurfFaces(mesh.nFaces()/100);
DynamicList<point> isoSurfPts(mesh.nPoints()/100);
DynamicList<label> zoneIDs(mesh.nFaces()/100);
label nIsoSurfPts(0);
// Number of patches before processor ones
label nSerialPatches = mesh.boundaryMesh().nNonProcessor();
// Map between iso-surface and mesh faces (internal and boundary)
Map<label> addedFaces;
// Loop over cells and check whether they are cut by zero iso-surface
const cellList& cells = mesh.cells();
forAll(cells, celli)
{
if (!cutCell.calcSubCell(celli, isoValue))
{
const vector& cuttingCf = cutCell.faceCentre();
vector cuttingNf = cutCell.faceArea();
cuttingNf.normalise();
// Loop over faces of the cut cell and check whether they are cut
// too, setting their distance from the cutting face
for (const label facei : cells[celli])
{
if (!cutFace.calcSubFace(facei, isoValue))
{
const vector& Cf = mesh.Cf()[facei];
const scalar distSqr =
magSqr((Cf - cuttingCf) & cuttingNf);
const vector lsPoint =
Cf + ((cuttingCf - Cf) & cuttingNf)*cuttingNf;
if (!addedFaces.found(facei))
{
// If the face being cut is a boundary one, part of
// it belongs to the iso-surface
bool addedToIsoSurf = addCutBoundaryFaceToIsoline
(
facei,
cutFace,
isoSurfPts,
isoSurfFaces,
zoneIDs,
nIsoSurfPts
);
if (mesh.isInternalFace(facei) || addedToIsoSurf)
{
addedFaces.insert(facei, nSurfFaces);
}
changedFacesInfo[nSurfFaces] =
wallPoint(lsPoint, distSqr);
changedFaces[nSurfFaces++] = facei;
}
else
{
const label visitedFace = addedFaces.at(facei);
if (distSqr < changedFacesInfo[visitedFace].distSqr())
{
changedFacesInfo[visitedFace] =
wallPoint(lsPoint, distSqr);
}
}
}
}
addCuttingFaceToIsoline
(
cutCell.facePoints(),
isoSurfPts,
isoSurfFaces,
zoneIDs,
nIsoSurfPts,
nSerialPatches
);
}
}
// Insert wall faces if they belong to the outer boundary
labelHashSet wallPatchIDs =
mesh.boundaryMesh().findPatchIDs<wallPolyPatch>();
for (const label patchi : wallPatchIDs)
{
const fvPatch& patch = mesh.boundary()[patchi];
const labelList& faceCells = patch.faceCells();
const label start = patch.start();
const vectorField& Cf = patch.Cf();
forAll(Cf, facei)
{
if
(
indicator[faceCells[facei]] >= 0 &&
!addedFaces.found(start + facei)
)
{
changedFaces[nSurfFaces] = start + facei;
changedFacesInfo[nSurfFaces] = wallPoint(Cf[facei], 0);
++nSurfFaces;
}
}
}
// Add boundary faces on the initial domain with a positive values on all
// points to the zero iso-surface
addBoundaryFacesToIsoline
(
pointY,
addedFaces,
isoValue,
isoSurfPts,
isoSurfFaces,
zoneIDs,
nIsoSurfPts
);
changedFaces.setSize(nSurfFaces);
changedFacesInfo.setSize(nSurfFaces);
// vtp and multi-region stl files holding the current geometry
writeSurfaceFiles(isoSurfPts, isoSurfFaces, zoneIDs, nSerialPatches);
}
// ************************************************************************* //

229
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/topOVariablesBase/topOVariablesBase.H Normal file
View File

@ -0,0 +1,229 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2007-2023 PCOpt/NTUA
Copyright (C) 2013-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::topOVariablesBase
Description
Base class for all design variables related to topology optimisation (topO).
Provides the lookup functionality necessary for all fvOptions associated
with topO and access to topOZones.
SourceFiles
topOVariablesBase.C
\*---------------------------------------------------------------------------*/
#ifndef topOVariablesBase_H
#define topOVariablesBase_H
#include "localIOdictionary.H"
#include "topOZones.H"
#include "topOInterpolationFunction.H"
#include "FaceCellWave.H"
#include "wallPoint.H"
#include "cutFaceIso.H"
#include "cutCellIso.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class topOVariablesBase Declaration
\*---------------------------------------------------------------------------*/
class topOVariablesBase
:
public localIOdictionary
{
protected:
// Protected data
//- Cell zones usefull for defining the constant and changing parts
//- of the domain in topO
topOZones zones_;
//- Folder name holding the zero level-set iso-surface
fileName isoSurfFolder_;
// Protected Member Functions
//- Based on the current pseudo-distance field y, identity the zero
//- level-set iso-surface, write it to files for post-processing and
//- set the seeds for the distance calculation from this iso-surface
void setCutInterfaceFaces
(
const volScalarField& y,
labelList& changedFaces,
List<wallPoint>& changedFacesInfo
) const;
//- Check whether the cutFace intersects the boundary of the initial
//- domain and add fluid part of the intersected face to the isoline
bool addCutBoundaryFaceToIsoline
(
const label facei,
const cutFaceIso& cutFace,
DynamicList<vector>& isoSurfPts,
DynamicList<face>& isoSurfFaces,
DynamicList<label>& zoneIDs,
label& nIsoSurfPts
) const;
//- Add the cutting face to the zero level-set iso-surface
bool addCuttingFaceToIsoline
(
const DynamicList<point>& facePoints,
DynamicList<vector>& isoSurfPts,
DynamicList<face>& isoSurfFaces,
DynamicList<label>& zoneIDs,
label& nIsoSurfPts,
const label nSerialPatches
) const;
//- Check whether the boundary faces of the initial domain belong
//- to the fluid part and add them to the surface describing the
//- fluid domain.
// Used only to close the output surface, does not set distance seeds
void addBoundaryFacesToIsoline
(
const pointScalarField& pointY,
const Map<label>& addedFaces,
const scalar isoValue,
DynamicList<vector>& isoSurfPts,
DynamicList<face>& isoSurfFaces,
DynamicList<label>& zoneIDs,
label& nIsoSurfPts
) const;
//- Write the surface describing the fluid domain to stl and vtp files
void writeSurfaceFiles
(
DynamicList<vector>& patchPoints,
DynamicList<face>& patchFaces,
DynamicList<label>& zoneIDs,
const label nSerialPatches
) const;
private:
// Private Member Functions
//- Disallow default bitwise copy construct
topOVariablesBase(const topOVariablesBase&) = delete;
//- Disallow default bitwise assignment
void operator=(const topOVariablesBase&) = delete;
public:
//- Runtime type information
TypeName("topOVariablesBase");
// Constructors
//- Construct from dictionary
topOVariablesBase
(
fvMesh& mesh,
const dictionary& dict
);
//- Destructor
virtual ~topOVariablesBase() = default;
// Member Functions
//- Get topOZones
inline const topOZones& getTopOZones() const
{
return zones_;
}
//- Get betaMax value
inline scalar getBetaMax() const
{
return zones_.getBetaMax();
}
//- Get field used for physical interpolations
virtual const volScalarField& beta() const = 0;
//- Populate source terms for the flow equations
virtual void sourceTerm
(
DimensionedField<scalar, volMesh>& field,
const topOInterpolationFunction& interpolationFunc,
const scalar betaMax,
const word& interpolationFieldName = "beta"
) const;
//- Post-processing sensitivities due to interpolations based on
//- the indicator fields
virtual void sourceTermSensitivities
(
scalarField& sens,
const topOInterpolationFunction& interpolationFunc,
const scalar betaMax,
const word& designVariablesName,
const word& interpolationFieldName = "beta"
) const;
//- Write the fluid-solid interface to files.
// The interface is computed as an isoValue contour of the indicator
// field
// - 0 distance contour for levelSet or
// - 0.5 beta contour for density-based topO.
// For levelSet topO, the process of identitying the contour sets also
// the seeds for computing the distance field in the entire domain
void writeFluidSolidInterface
(
const volScalarField& indicator,
const scalar isoValue,
labelList& changedFaces,
List<wallPoint>& changedFacesInfo
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

118
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/topOZones/topOZones.C Normal file
View File

@ -0,0 +1,118 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2007-2023 PCOpt/NTUA
Copyright (C) 2013-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "topOZones.H"
#include "cellSet.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(topOZones, 0);
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
Foam::labelList Foam::topOZones::getZoneIDs
(
const dictionary& dict,
const word& zoneGroup
)
{
wordList zoneNames(dict.getOrDefault<wordList>(zoneGroup, wordList(0)));
labelList IDs(zoneNames.size(), -1);
forAll(zoneNames, zI)
{
IDs[zI] = mesh_.cellZones().findZoneID(zoneNames[zI]);
}
return IDs;
}
void Foam::topOZones::addIOcellsZone()
{
// Gather IO cells
DynamicList<label> IOcells(0);
for (const fvPatch& patch : mesh_.boundary())
{
if (patch.type() == "patch")
{
IOcells.append(patch.faceCells());
}
}
// Add zone to cellZoneMesh and populate it
cellZoneMesh& cellZones = const_cast<fvMesh&>(mesh_).cellZones();
cellZone& IOcellsZone = cellZones("IOcells", true);
IOcellsZone = IOcells;
IOcellsID_ = cellZones.size() - 1;
// Write as set
cellSet IOcellList(mesh_, "IOcellList", IOcells);
IOcellList.write();
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::topOZones::topOZones
(
const fvMesh& mesh,
const dictionary& dict
)
:
mesh_(mesh),
dict_(dict),
fixedPorousIDs_(getZoneIDs(dict, "fixedPorousZones")),
fixedPorousValues_
(
dict.getOrDefault<scalarList>
(
"fixedPorousValues", scalarList(fixedPorousIDs_.size(), 1.)
)
),
fixedZeroPorousIDs_(getZoneIDs(dict, "fixedZeroPorousZones")),
adjointPorousIDs_(getZoneIDs(dict, "adjointPorousZones")),
IOcellsID_(-1),
betaMaxPtr_(betaMax::New(mesh, dict))
{
addIOcellsZone();
if (fixedPorousValues_.size() != fixedPorousIDs_.size())
{
FatalErrorInFunction
<< "Number of fixedPorousValues and fixedPorousZones don't agree!"
<< "\nSize of fixedPorousIDs is " << fixedPorousIDs_.size()
<< " and size of fixedPorousValues is " << fixedPorousValues_.size()
<< endl << endl
<< exit(FatalError);
}
}
// ************************************************************************* //

190
src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/topOZones/topOZones.H Normal file
View File

@ -0,0 +1,190 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2007-2023 PCOpt/NTUA
Copyright (C) 2013-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
topOZones
Description
Class holding cells zones used in topology optimisation
SourceFiles
topOZones.C
\*---------------------------------------------------------------------------*/
#ifndef topOZones_H
#define topOZones_H
#include "fvMesh.H"
#include "betaMax.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class topOZones Declaration
\*---------------------------------------------------------------------------*/
class topOZones
{
protected:
// Protected data
//- Mesh reference
const fvMesh& mesh_;
//- TopO design variables dictionary
const dictionary dict_;
//- IDs of cellZones holding cells with constant alpha values
labelList fixedPorousIDs_;
//- The constant alpha values of fixedPorousIDs_
scalarList fixedPorousValues_;
//- IDs of cellZones holding cells with zero alpha values
labelList fixedZeroPorousIDs_;
//- IDs of cellZones with cells that can change their alpha value
//- throughout the optimisation
labelList adjointPorousIDs_;
//- ID of the cellZone holding the cells next to inlets & outlets
label IOcellsID_;
//- The multiplier of the Brinkman penalisation term
autoPtr<betaMax> betaMaxPtr_;
// Protected Member Functions
//- Get zone IDs corresponding to a wordList, read from a dict.
// Avoid going through ZoneMesh.indices() since this practically sorts
// the IDs from smallest to largest while we need to keep them in the
// same order as that perscribed in the wordList
labelList getZoneIDs(const dictionary& dict, const word& zoneGroup);
//- Add a cellZone containing the cells next to IO patches
void addIOcellsZone();
private:
// Private Member Functions
//- No copy construct
topOZones(const topOZones&);
//- No copy assignment
void operator=(const topOZones&);
public:
//- Runtime type information
TypeName("topOZones");
// Constructors
//- Construct from components
topOZones
(
const fvMesh& mesh,
const dictionary& dict
);
//- Destructor
virtual ~topOZones() = default;
// Member Functions
//- Const reference to mesh
inline const fvMesh& mesh() const
{
return mesh_;
}
//- Cell zone IDs with fixed porosity values
inline const labelList& fixedPorousZoneIDs() const
{
return fixedPorousIDs_;
}
//- Values of alpha for each fixed porous zone
inline const scalarList& fixedPorousValues() const
{
return fixedPorousValues_;
}
//- Cell zone IDs with fixed zero porosity values
inline const labelList& fixedZeroPorousZoneIDs() const
{
return fixedZeroPorousIDs_;
}
//- Cell zone IDs in which porosity is allowed to change
inline const labelList& adjointPorousZoneIDs() const
{
return adjointPorousIDs_;
}
//- Cells next to IO boundaries
inline const cellZone& IOCells() const
{
return mesh_.cellZones()[IOcellsID_];
}
//- ID of the cellZone holding the IOcells
inline label IOzoneID() const
{
return IOcellsID_;
}
//- Get betaMax
inline scalar getBetaMax() const
{
return betaMaxPtr_().value();
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

185
src/optimisation/adjointOptimisation/adjoint/optimisation/lineSearch/GCMMA/GCMMA.C Normal file
View File

@ -0,0 +1,185 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
#include "GCMMA.H"
#include "IOmanip.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(GCMMA, 1);
addToRunTimeSelectionTable
(
lineSearch,
GCMMA,
dictionary
);
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
void Foam::GCMMA::writeToFiles(const bool isConverged)
{
const PtrList<scalarField>& objValues = mma_.getValuesAndApproximations();
const scalarField& rhoValues = mma_.getRho();
const label m(rhoValues.size() - 1);
if (Pstream::master())
{
unsigned int width = IOstream::defaultPrecision() + 5;
if (writeHeader_)
{
GCMMAFile_
<< setw(width) << "#OuterIter" << " "
<< setw(width) << "InnerIter" << " "
<< setw(width) << "rhoObj" << " ";
costFile_
<< setw(width) << "#nCycle" << " "
<< setw(width) << "cumulativeCost" << " "
<< setw(width) << "Objective" << " ";
for (label i = 0; i < m; ++i)
{
GCMMAFile_
<< setw(width) << "rhoConst" << " ";
costFile_
<< setw(width) << "Constraint" << " ";
}
GCMMAFile_
<< setw(width) << "J" << " "
<< setw(width) << "JTilda" << " ";
for (label i = 0; i < m; ++i)
{
GCMMAFile_
<< setw(width) << "C" << " "
<< setw(width) << "CTilda" << " ";
}
GCMMAFile_<< endl;
costFile_<< endl;
writeHeader_ = false;
}
GCMMAFile_
<< setw(width) << iter_ + 2 << " "
<< setw(width) << innerIter_ << " ";
forAll(rhoValues, i)
{
GCMMAFile_
<< setw(width) << rhoValues[i] << " ";
}
forAll(rhoValues, i)
{
GCMMAFile_
<< setw(width) << objValues[0][i] << " "
<< setw(width) << objValues[1][i] << " ";
}
GCMMAFile_ << endl;
if (isConverged)
{
// The cost of this cycle is equal to the number of innerIters
// plus one for the adjoint solution
cost_ += innerIter_ + 1;
costFile_
<< setw(width) << iter_ + 2 << " "
<< setw(width) << cost_ << " ";
forAll(rhoValues, i)
{
costFile_
<< setw(width) << objValues[0][i] << " ";
}
costFile_<< endl;
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::GCMMA::GCMMA
(
const dictionary& dict,
const Time& time,
updateMethod& UpdateMethod
)
:
lineSearch(dict, time, UpdateMethod),
mma_(refCast<MMA>(UpdateMethod)),
GCMMAFile_
(time.globalPath()/"optimisation"/"objective"/time.timeName()/"GCMMA"),
costFile_
(
time.globalPath()/"optimisation"/"objective"/time.timeName()
/"GCMMACost"
),
cost_(2),
writeHeader_(true)
{
// Force rho constants to be updated in each optimisation cycle
mma_.setVariableRho();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * //
bool Foam::GCMMA::converged()
{
// Calls MMA::updateValuessAndApproximations() and needs to be called
// before writeToFiles
bool isConverged = mma_.converged();
writeToFiles(isConverged);
DebugInfo
<< "GCMMA converged ... " << Switch(isConverged) << endl;
return isConverged;
}
void Foam::GCMMA::updateStep()
{
DebugInfo
<< "GCMMA:: recomputing direction "<< endl;
// Update rho values
mma_.updateRho();
// Re-solve the subproblem
mma_.solveSubproblem();
}
void Foam::GCMMA::updateCorrection(scalarField& correction)
{
correction = mma_.returnCorrection();
}
// ************************************************************************* //

157
src/optimisation/adjointOptimisation/adjoint/optimisation/lineSearch/GCMMA/GCMMA.H Normal file
View File

@ -0,0 +1,157 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::GCMMA
Description
Class implementing a workaround for using the Globally converging wrapper
of the MMA update method. Implemented as a lineSearch, even though the
search direction changes within each inner iteration of GCMMA, due the
existing framework for re-solving the primal equations and checking the
convergence of a merit function.
Reference:
\verbatim
Svanberg, K. (2002)
A class of globally convergent optimization methods based on
conservative convex separable approximations.
SIAM Journal of Optimization, 12), 555-573.
https://doi.org/10.1137/S1052623499362822
\endverbatim
SourceFiles
GCMMA.C
\*---------------------------------------------------------------------------*/
#ifndef GCMMA_H
#define GCMMA_H
#include "lineSearch.H"
#include "MMA.H"
#include "OFstream.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class GCMMA Declaration
\*---------------------------------------------------------------------------*/
class GCMMA
:
public lineSearch
{
protected:
// Protected data
//- Cast of the update method to MMA
MMA& mma_;
//- File with rho values in each inner loop
OFstream GCMMAFile_;
//- File with CPU cost in each inner loop
OFstream costFile_;
//- Total cost, measured in Equivalent Flow Solutions, up to this cycle
label cost_;
//- Write the header of the output files
bool writeHeader_;
// Protected Member Functions
//- Write rho and objective/constraint values & approx to files
void writeToFiles(bool isConverged);
private:
// Private Member Functions
//- No copy construct
GCMMA(const GCMMA&) = delete;
//- No copy assignment
void operator=(const GCMMA&) = delete;
public:
//- Runtime type information
TypeName("GCMMA");
// Constructors
//- Construct from components
GCMMA
(
const dictionary& dict,
const Time& time,
updateMethod& UpdatheMethod
);
// Destructor
virtual ~GCMMA() = default;
// Member Functions
//- Check whether linear search has converged
virtual bool converged();
//- Actually computes a new direction entirely,
//- targeting the satisfaction of the GCMMA condition
virtual void updateStep();
//- Replace the correction with the one computed in updateStep
virtual void updateCorrection(scalarField& correction);
// Member operators
//- Increment iteration number and store merit value corresponding to
//- the previous optimisation cycle
//virtual GCMMA& operator++();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

1110
src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/MMA/MMA.C Normal file
View File

File diff suppressed because it is too large Load Diff

416
src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/MMA/MMA.H Normal file
View File

@ -0,0 +1,416 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020-2023 PCOpt/NTUA
Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
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/>.
Class
Foam::MMA
Description
Update design variables using the Method of Moving Assymptotes.
Can handle inequality constraints.
Reference:
\verbatim
Svanberg, K. (1987)
The method of moving asymptotes—a new method for structural
optimization
International Journal for Numerical Methods in Engineering, 24(2),
359-373.
https://doi.org/10.1002/nme.1620240207
\endverbatim
SourceFiles
MMA.C
\*---------------------------------------------------------------------------*/
#ifndef MMA_H
#define MMA_H
#include "constrainedOptimisationMethod.H"
#include "updateMethod.H"
#include "PtrList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class MMA Declaration
\*---------------------------------------------------------------------------*/
class MMA
:
public constrainedOptimisationMethod,
public updateMethod
{
protected:
// Protected data
//- The previous design variables. Used to compute new bounds
scalarField x0_;
//- The twice previous design variables. Used to compute new bounds
scalarField x00_;
//- The set of design variables being updated during the subproblem.
// Correction will end up being the difference of this and the
// old values of the design variables
scalarField xNew_;
//- Old objective value
// Needed by GCMMA
scalar oldObjectiveValue_;
//- Old constraint values
// Needed by GCMMA
scalarField oldCValues_;
//- Objective/Constraint values and approximations in the new point
// Needed by GCMMA
PtrList<scalarField> valsAndApproxs_;
//- Second term in the approximation function
scalar z_;
//- Term multiplying z in the objective function
scalar alpha0_;
//- Terms multyplying z in the constraint functions
scalarField alpha_;
//- y in the constraints functions
scalarField y_;
//- y multipliers in the objective function
scalarField c_;
//- y^2 multipliers in the objective function
scalarField d_;
//- Lower design variables asymptotes
scalarField lower_;
//- Upper design variables asymptotes
scalarField upper_;
//- Lower design variables constraints
scalarField a_;
//- Upper design variables constraints
scalarField b_;
//- Constants in the (p,q) functions.
// Size is cValues_.size() + 1 (the objective)
scalarField rho_;
//- Bound the rho value with an upper value?
bool boundRho_;
//- Correct the design variables
bool correctDVs_;
// Dual variables
//- Constraint Lagrange multipliers
scalarField lamda_;
//- Lagrange multipliers for the lower limits constraints
scalarField ksi_;
//- Lagrange multipliers for the upper limits constraints
scalarField Eta_;
//- Lagrange multipliers for the y constraints
scalarField mu_;
//- Lagrange multiplier for the z constraint
scalar zeta_;
//- Slack variables for the inequality constraints
scalarField s_;
// Fields holding updates of the design, Lagrange and slack variables
scalarField deltaLamda_;
scalar deltaZ_;
scalarField deltaX_;
scalarField deltaY_;
scalarField deltaS_;
scalar deltaZeta_;
scalarField deltaMu_;
scalarField deltaEta_;
scalarField deltaKsi_;
//- Infinitesimal quantity
// Updated during the inner iterations of the subproblem
scalar eps_;
//- Maxmimum number of Newton iterations for the subproblem
label maxNewtonIters_;
//- Maxmimum number of line search iterations for each iteration of the
//- subproblem
label maxLineSearchIters_;
//- Initialize every subproblem with x = 0.5*(a + b)
//- and reset Lagrange multipliers
// The solution of the subproblem can become relatively expensive
// in some cases, hence the GCMMA subproblems are usually
// initialized with the solution of the previous inner iteration,
// to save up some time. Occassionally, this can cause the solution
// of the subproblem to be trapped in a non-feasible point.
// In such cases, the solution of each subproblem could be
// initialized from scratch
bool initializeEverySubproblem_;
// Constant parameters
//- Lower assymprote reduction rate
scalar asymDecr_;
//- Upper assymprote increase rate
scalar asymIncr_;
//- Movement of the unatainable upper and lower bounds
scalar move_;
//- Constant in p, q functions
scalar raa0_;
// Parameters related to the update of rho in GCMMA
//- Multiplier of the mean derivative used for the initialisation
//- of rho
scalar maxInitRhoMult_;
//- Multiplier for the max. rho value during its update
scalar maxRhoMult_;
//- Change rho constants in each iteration?
bool variableRho_;
//- Change rho constants in each iteration?
bool dynamicRhoInitialisation_;
//- The rho values obtained by the dynamic rho initialisation
//- might be too conservative. Multiply with this number to relax
//- them
scalar dynamicRhoMult_;
// Protected Member Functions
//- Update sizes of fields related to the active design variables
void updateSizes();
//- Initialize rho constants in the (p, q) functions
// These control how aggressive or conservative the method will be
void initializeRho();
//- Update the bounds associated with the design variables
void updateBounds();
//- Allocate fields related to constraints
void initialize();
//- Store old objcective and constraint values
// Needed by GCMMA
void storeOldValues();
//- p and q functions, used to approximate the objective and contraints
// Computed based on the current set of design variables, not updated
// during the iterations of the subproblem
tmp<scalarField> p
(
const scalarField& derivs,
const scalar r,
const scalarField& x
);
tmp<scalarField> p(const scalarField& derivs, const scalar r);
tmp<scalarField> q
(
const scalarField& derivs,
const scalar r,
const scalarField& x
);
tmp<scalarField> q(const scalarField& derivs, const scalar r);
//- g of all constraint functions
// Computed using the current set of design variables, updated during
// the solution of the subproblem
tmp<scalarField> gConstr(const scalarField& vars);
//- The rhs of the contraint approximations.
// Always computed using the current (frozen) set of design variables
tmp<scalarField> b();
//- p and q functions, using the dual lamda
tmp<scalarField> pLamda();
tmp<scalarField> qLamda();
//- Zero the updates computed in the previous optimisation cycle
void zeroUpdates();
//- Computation of the approximate function
scalar fTilda
(
const scalarField& xInit,
const scalarField& x,
const scalar f,
const scalarField& dfdx,
const scalar rho
);
//- Read in scalarField from self (potentially in binary), or allocate
//- field with the size of the active design variables and given value
tmp<scalarField> getOrDefaultScalarField
(
const word& name,
const label size,
const scalar value = Zero
);
//- Read in scalarField from self (potentially in binary), or allocate
//- field with the size of the active design variables and given value
void setOrDefaultScalarField
(
scalarField& field,
const word& name,
const label size,
const scalarField& defaultField
);
// Functions related to the solution of the primal-dual subproblem
//- Compute direction for the Newton problem
void computeNewtonDirection();
//- Perform line search and return max residual corresponding to
//- the updated solution
scalar lineSearch();
//- Update the current solution using the known direction and the
//- given step length
void updateSolution(const scalar step);
//- Adjust step to satisfy cireteria
void adjustStep
(
scalar& step,
const scalar value,
const scalar update
);
//- Compute and return residuals based on the current solution
tmp<scalarField> computeResiduals();
private:
// Private Member Functions
//- No copy construct
MMA(const MMA&);
//- No copy assignment
void operator=(const MMA&);
public:
//- Runtime type information
TypeName("MMA");
// Constructors
//- Construct from components
MMA
(
const fvMesh& mesh,
const dictionary& dict,
autoPtr<designVariables>& designVars,
const label nConstraints,
const word& type
);
//- Destructor
virtual ~MMA() = default;
// Member Functions
//- Compute design variables correction
void computeCorrection();
//- Solve subproblem using a Newton optimiser
void solveSubproblem();
// Functions needed by GCMMA
//- Compute objective/constraint values and their approximations
// Needed to check convergence and update rho, if necessary
void updateValuesAndApproximations();
//- Get objective/constraint values and their approximations
const PtrList<scalarField>& getValuesAndApproximations() const;
//- Update rho value if needed
void updateRho();
//- Get rho value if needed
const scalarField& getRho() const;
//- Set variable rho
void setVariableRho(bool varRho = true);
//- Checks convergence of GCMMA
bool converged();
//- Write useful quantities to files
virtual bool writeData(Ostream& os) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

12
src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/adjointSolver/adjointSolver.H
View File

@ -303,6 +303,18 @@ public:
const scalar dt
)
{}
// Topology optimisation
//- Compute the multiplier of beta
virtual void topOSensMultiplier
(
scalarField& betaMult,
const word& designVariablesName,
const scalar dt
)
{}
};

62
src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/incompressibleAdjointSolver/incompressibleAdjointSolver.C
View File

@ -30,6 +30,7 @@ License
#include "incompressibleAdjointSolver.H"
#include "incompressiblePrimalSolver.H"
#include "wallFvPatch.H"
#include "sensitivityTopO.H"
#include "adjointBoundaryConditions.H"
#include "adjointBoundaryConditionsFwd.H"
#include "addToRunTimeSelectionTable.H"
@ -460,4 +461,65 @@ void Foam::incompressibleAdjointSolver::accumulateOptionsDxDbMultiplier
}
void Foam::incompressibleAdjointSolver::topOSensMultiplier
(
scalarField& betaMult,
const word& designVariablesName,
const scalar dt
)
{
const incompressibleAdjointVars& adjointVars = getAdjointVars();
const volVectorField& U = primalVars_.U();
const volVectorField& Ua = adjointVars.Ua();
const autoPtr<incompressibleAdjoint::adjointRASModel>& adjointTurbulence =
adjointVars.adjointTurbulence();
fv::options& fvOptions(fv::options::New(mesh_));
// Term from the momentum equations
scalarField momSens((U.primitiveField() & Ua.primitiveField())*dt);
Foam::sensitivityTopO::postProcessSens
(betaMult, momSens, fvOptions, U.name(), designVariablesName);
if (debug > 2)
{
volScalarField IvSens
(
IOobject
(
"IvSens" + solverName(),
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar(dimless, Zero)
);
IvSens.primitiveFieldRef() = momSens;
IvSens.write();
}
// Term from the turbulence model.
// Includes already the derivative of the interpolation function
betaMult +=
(adjointTurbulence->topologySensitivities(designVariablesName))*dt;
// Terms resulting directly from the objective function
PtrList<objective>& functions = objectiveManager_.getObjectiveFunctions();
for (objective& objI : functions)
{
const scalar weight(objI.weight());
if (objI.hasdJdb())
{
betaMult += weight*objI.dJdb()*dt;
}
if (objI.hasdJdbField())
{
SubField<scalar> betaSens(objI.dJdbField(), mesh_.nCells(), 0);
betaMult += weight*betaSens*dt;
}
}
}
// ************************************************************************* //

11
src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/incompressibleAdjointSolver/incompressibleAdjointSolver.H
View File

@ -252,6 +252,17 @@ public:
);
// Topology optimisation
//- Compute the multiplier of beta
virtual void topOSensMultiplier
(
scalarField& betaMult,
const word& designVariablesName,
const scalar dt
);
// IO
//- In case of multi-point runs with turbulent flows,

26
src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/null/adjointNull.C
View File

@ -235,4 +235,30 @@ void Foam::adjointNull::accumulateGeometryVariationsMultipliers
}
void Foam::adjointNull::topOSensMultiplier
(
scalarField& betaMult,
const word& designVariablesName,
const scalar dt
)
{
// Terms resulting directly from the objective function
PtrList<objective>& functions = objectiveManager_.getObjectiveFunctions();
for (objective& objI : functions)
{
const scalar weight = objI.weight();
if (objI.hasdJdb())
{
betaMult += weight*objI.dJdb()*dt;
}
if (objI.hasdJdbField())
{
SubField<scalar> betaSens(objI.dJdbField(), mesh_.nCells(), 0);
betaMult += weight*betaSens*dt;
}
}
}
// ************************************************************************* //

11
src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/null/adjointNull.H
View File

@ -171,6 +171,17 @@ public:
const labelHashSet& sensitivityPatchIDs,
const scalar dt
);
// Topology optimisation
//- Compute the multiplier of beta
virtual void topOSensMultiplier
(
scalarField& betaMult,
const word& designVariablesName,
const scalar dt
);
};

13
src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointLaminar/adjointLaminar.C
View File

@ -5,8 +5,8 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2007-2019 PCOpt/NTUA
Copyright (C) 2013-2019 FOSS GP
Copyright (C) 2007-2023 PCOpt/NTUA
Copyright (C) 2013-2023 FOSS GP
Copyright (C) 2019-2023 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@ -185,6 +185,15 @@ tmp<volTensorField> adjointLaminar::FISensitivityTerm()
}
tmp<scalarField> adjointLaminar::topologySensitivities
(
const word& designVarsName
) const
{
return tmp<scalarField>::New(mesh_.nCells(), Zero);
}
void adjointLaminar::nullify()
{
// Does nothing. No fields to nullify

10
src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointLaminar/adjointLaminar.H
View File

@ -5,8 +5,8 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2007-2019 PCOpt/NTUA
Copyright (C) 2013-2019 FOSS GP
Copyright (C) 2007-2023 PCOpt/NTUA
Copyright (C) 2013-2023 FOSS GP
Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@ -128,6 +128,12 @@ public:
//- Returns zero field
virtual tmp<volTensorField> FISensitivityTerm();
//- Returns zero field
virtual tmp<scalarField> topologySensitivities
(
const word& designVarsName
) const;
//- Nullify all adjoint turbulence model fields and their old times
virtual void nullify();

9
src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointRASModel/adjointRASModel.H
View File

@ -331,6 +331,15 @@ public:
// Misses grad(dxdb), to be added by the assembling the sensitivities
virtual tmp<volTensorField> FISensitivityTerm() = 0;
//- Term contributing to the computation of topology optimisation
//- sensitivities
// Misses betaMax, dBeta/dAlpha and the mesh volume, to be added
// during the assembly of the sensitivities
virtual tmp<scalarField> topologySensitivities
(
const word& designVarsName
) const = 0;
//- Solve the adjoint turbulence equations
virtual void correct();

25
src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointSpalartAllmaras/adjointSpalartAllmaras.C
View File

@ -36,6 +36,7 @@ License
#include "coupledFvPatch.H"
#include "ATCModel.H"
#include "fvOptions.H"
#include "sensitivityTopO.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -1035,6 +1036,24 @@ tmp<volTensorField> adjointSpalartAllmaras::FISensitivityTerm()
return tFISens;
}
tmp<scalarField> adjointSpalartAllmaras::topologySensitivities
(
const word& designVarsName
) const
{
auto tres(tmp<scalarField>::New(nuTilda().primitiveField().size(), Zero));
scalarField nuTildaSens
(nuTilda().primitiveField()*nuaTilda().primitiveField());
fv::options& fvOptions(fv::options::New(mesh_));
sensitivityTopO::postProcessSens
(
tres.ref(), nuTildaSens, fvOptions, nuTilda().name(), designVarsName
);
return tres;
}
@ -1070,7 +1089,7 @@ void adjointSpalartAllmaras::correct()
nuaTilda().storePrevIter();
fv::options& fvOptions(fv::options::New(this->mesh_));
fv::options& fvOptions(fv::options::New(mesh_));
tmp<fvScalarMatrix> nuaTildaEqn
(
@ -1082,9 +1101,9 @@ void adjointSpalartAllmaras::correct()
+ fvc::laplacian(2.0*Cb2_*oneOverSigmaNut*nuaTilda(), nuTilda())
+ gradNua*oneOverSigmaNut
==
// always a negative contribution to the lhs. No Sp used!
// Always a negative contribution to the lhs. No Sp used!
Cb1_*Stilda_*nuaTilda()
//always a positive contribution to the lhs. no need for SuSp
// Always a positive contribution to the lhs. no need for SuSp
- fvm::Sp(Cw1_*fw_*nuTilda()/sqr(y_), nuaTilda())
- Cdnut_*gradUaR
+ fvOptions(nuaTilda())

5
src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointSpalartAllmaras/adjointSpalartAllmaras.H
View File

@ -343,6 +343,11 @@ public:
virtual tmp<volTensorField> FISensitivityTerm();
virtual tmp<scalarField> topologySensitivities
(
const word& designVarsName
) const;
//- Nullify all adjoint turbulence model fields and their old times
virtual void nullify();

31
src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSST.C
View File

@ -34,6 +34,8 @@ License
#include "reverseLinear.H"
#include "nutkWallFunctionFvPatchScalarField.H"
#include "omegaWallFunctionFvPatchScalarField.H"
#include "fvOptions.H"
#include "sensitivityTopO.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -2149,6 +2151,8 @@ void adjointkOmegaSST::correct()
volScalarField dR_dnut(this->dR_dnut());
volScalarField::Internal divU(fvc::div(fvc::absolute(phi, U)));
fv::options& fvOptions(fv::options::New(mesh_));
tmp<fvScalarMatrix> waEqn
(
fvm::div(-phi, wa())
@ -2164,6 +2168,8 @@ void adjointkOmegaSST::correct()
wa()
)
- (case_2_Pk_*c1_ - scalar(1))*betaStar_*k()*ka()
==
fvOptions(wa())
);
// Boundary manipulate changes the diagonal component, so relax has to
@ -2174,6 +2180,7 @@ void adjointkOmegaSST::correct()
// Sources from the objective should be added after the boundary
// manipulation
objectiveManager_.addSource(waEqn.ref());
fvOptions.constrain(waEqn.ref());
waEqn.ref().solve();
// Adjoint Turbulent kinetic energy equation
@ -2189,11 +2196,14 @@ void adjointkOmegaSST::correct()
+ kaEqnSourceFromF1()
+ dR_dnut()*dnut_dk_()
- zeroFirstCell_()*gamma_*dGPrime_dk().ref()()*wa()
==
fvOptions(ka())
);
kaEqn.ref().relax();
kaEqn.ref().boundaryManipulate(ka().boundaryFieldRef());
addWallFunctionTerms(kaEqn.ref(), dR_dnut);
fvOptions.constrain(kaEqn.ref());
// Add sources from the objective functions
objectiveManager_.addSource(kaEqn.ref());
@ -2345,8 +2355,25 @@ tmp<scalarField> adjointkOmegaSST::topologySensitivities
const word& designVarsName
) const
{
// Missing proper terms - return zero for now
return tmp<scalarField>::New(mesh_.nCells(), Zero);
fv::options& fvOptions(fv::options::New(this->mesh_));
auto tres(tmp<scalarField>::New(mesh_.nCells(), Zero));
// Sensitivity from the source term in the k equation
scalarField auxSens
(k().primitiveField()*ka().primitiveField());
sensitivityTopO::postProcessSens
(
tres.ref(), auxSens, fvOptions, k().name(), designVarsName
);
// Sensitivity from the source term in the omega equation
auxSens = omega().primitiveField()*wa().primitiveField();
sensitivityTopO::postProcessSens
(
tres.ref(), auxSens, fvOptions, omega().name(), designVarsName
);
return tres;
}