47namespace incompressible
73 createZeroBoundaryPointFieldPtr<vector>(
mesh_)
77 createZeroBoundaryPointFieldPtr<vector>(
mesh_)
81 for (
const label patchI : sensitivityPatchIDs_)
92 for (
auto& fun : functions)
94 dSdbMultiplierTot += fun.weight()*fun.dSdbMultiplier(patchI);
95 dndbMultiplierTot += fun.weight()*fun.dndbMultiplier(patchI);
99 const labelList& meshPoints = patch.patch().meshPoints();
104 const labelListList& patchPointFaces = patch.patch().pointFaces();
106 const label patchStartIndex = patch.start();
114 const labelList& pointFaces = patchPointFaces[ppI];
115 for (label localFaceIndex : pointFaces)
117 label globalFaceIndex = patchStartIndex + localFaceIndex;
118 const face& faceI = faces[globalFaceIndex];
124 if (faceI[facePointI] == meshPoints[ppI])
135 const tensor& deltaSf = deltaNormals[1];
137 dSdbMultiplierTot[localFaceIndex] & deltaSf;
140 const tensor& deltaNf = deltaNormals[2];
142 dndbMultiplierTot[localFaceIndex] & deltaNf;
150 for (
const label patchI : sensitivityPatchIDs_)
156 const label globaPointI = meshPoints[ppI];
157 dSdbGlobal[globaPointI] += pointSensdSdb()[patchI][ppI];
158 dndbGlobal[globaPointI] += pointSensdndb()[patchI][ppI];
171 for (
const label patchI : sensitivityPatchIDs_)
174 const labelList& meshPoints = patch.patch().meshPoints();
176 pointSensdSdb()[patchI].map(dSdbGlobal, meshPoints);
177 pointSensdndb()[patchI].map(dndbGlobal, meshPoints);
194 wallFaceSensVecPtr_()[patchI] += dSdbFace + dndbFace;
224 const label
iters(
dict().getOrDefault<label>(
"iters", 500));
225 const scalar tolerance(
dict().getOrDefault<scalar>(
"tolerance", 1.e-06));
246 Info<<
"Reading the already constructed faMesh" <<
endl;
266 Info<<
"Reading faMeshDefinition from system " <<
endl;
273 Info<<
"Constructing faMeshDefinition from sensitivity patches"
275 wordList polyMeshPatches(sensitivityPatchIDs_.size());
277 for (
const label
patchID : sensitivityPatchIDs_)
281 faMeshDefinition.
add<
wordList>(
"polyMeshPatches", polyMeshPatches);
292 const scalar Rphysical
295 <<
"Physical radius of the sensitivity smoothing "
296 << Rphysical <<
nl <<
endl;
333 for (label iter = 0; iter <
iters; ++iter)
335 Info<<
"Sensitivity smoothing iteration " << iter <<
endl;
347 const scalar residual(
mag(smoothEqn.
solve().initialResidual()));
350 <<
"Max smoothSens " <<
gMax(
mag(smoothedSens)()) <<
endl;
356 if (residual < tolerance)
358 Info<<
"\n***Reached smoothing equation convergence limit, "
359 "iteration " << iter <<
"***\n\n";
369 "smoothedSurfaceSens" + surfaceFieldSuffix_,
381 volSmoothedSens.
write();
392 if (geometricD[iDir] == -1)
394 averageArea /= bounds.
span()[iDir];
414 includeSurfaceArea_(false),
415 includePressureTerm_(false),
416 includeGradStressTerm_(false),
417 includeTransposeStresses_(false),
418 useSnGradInTranposeStresses_(false),
419 includeDivTerm_(false),
420 includeDistance_(false),
421 includeMeshMovement_(false),
422 includeObjective_(false),
423 writeGeometricInfo_(false),
424 smoothSensitivities_(false),
425 eikonalSolver_(nullptr),
426 meshMovementSolver_(nullptr),
428 nfOnPatchPtr_(nullptr),
429 SfOnPatchPtr_(nullptr),
430 CfOnPatchPtr_(nullptr)
436 wallFaceSensVecPtr_.reset(createZeroBoundaryPtr<vector>(
mesh_));
437 wallFaceSensNormalPtr_.reset(createZeroBoundaryPtr<scalar>(
mesh_));
438 wallFaceSensNormalVecPtr_.reset(createZeroBoundaryPtr<vector>(
mesh_));
554 sensitivityPatchIDs_,
564 if (sensitivity::readDict(
dict))
586 for (
const label patchI : sensitivityPatchIDs_)
618 adjointTurbulence->wallShapeSensitivities();
621 <<
" Calculating adjoint sensitivity. " <<
endl;
641 for (
const label patchI : sensitivityPatchIDs_)
646 wallFaceSensVecPtr_()[patchI] -=
667 if (isA<wallFvPatch>(patch))
670 gradUbf[patchI] = tnf*
U.boundaryField()[patchI].snGrad();
679 for (label idir = 0; idir < pTraits<vector>::nComponents; ++idir)
683 for (
const label patchI : sensitivityPatchIDs_)
688 wallFaceSensVecPtr_()[patchI] +=
703 tgradUa.
cref().boundaryField();
705 for (
const label patchI : sensitivityPatchIDs_)
714 : (gradUabf[patchI] & nf)
716 wallFaceSensVecPtr_()[patchI] -=
718 *(gradUaNf &
U.boundaryField()[patchI].snGrad())*nf;
722 for (
const label patchI : sensitivityPatchIDs_)
745 &
U.boundaryField()[patchI].snGrad()
757 &
U.boundaryField()[patchI].snGrad()
769 &
U.boundaryField()[patchI].snGrad()
783 functions[funcI].weight()
785 functions[funcI].dxdbDirectMultiplier(patchI)
791 wallFaceSensVecPtr_()[patchI] +=
795 + adjointTMsensitivities[patchI]
802 (wallFaceSensVecPtr_(), sensitivityPatchIDs_, dt);
815 for (
const label patchI : sensitivityPatchIDs_)
835 distanceSensPtr.
reset(createZeroBoundaryPtr<vector>(
mesh_));
838 for (
const label patchI : sensitivityPatchIDs_)
840 distanceSensPtr()[patchI] = sens[patchI];
848 meshMovementSensPtr.
reset(createZeroBoundaryPtr<vector>(
mesh_));
851 for (
const label patchI : sensitivityPatchIDs_)
853 meshMovementSensPtr()[patchI] = sens[patchI];
859 label nPassedFaces(0);
860 for (
const label patchI : sensitivityPatchIDs_)
869 wallFaceSensVecPtr_()[patchI] += distanceSensPtr()[patchI];
875 wallFaceSensVecPtr_()[patchI] += meshMovementSensPtr()[patchI];
880 wallFaceSensVecPtr_()[patchI] *= patch.magSf();
883 wallFaceSensNormalPtr_()[patchI] = wallFaceSensVecPtr_()[patchI] & nf;
884 wallFaceSensNormalVecPtr_()[patchI] =
885 wallFaceSensNormalPtr_()[patchI] * nf;
890 = wallFaceSensNormalPtr_()[patchI][fI];
892 nPassedFaces += patch.size();
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
const Field< Type > & field() const
Return field.
tmp< Field< cmptType > > component(const direction) const
Return a component field of the field.
Generic GeometricBoundaryField class.
Internal::FieldType & primitiveFieldRef(const bool updateAccessTime=true)
Return a reference to the internal field.
Boundary & boundaryFieldRef(const bool updateAccessTime=true)
Return a reference to the boundary field.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Defines the attributes of an object for which implicit objectRegistry management is supported,...
bool typeHeaderOk(const bool checkType=true, const bool search=true, const bool verbose=true)
Read header (uses typeFilePath to find file) and check its info.
void setSize(const label n)
Alias for resize()
Interpolation class within a primitive patch. Allows interpolation from points to faces and vice vers...
tmp< Field< Type > > pointToFaceInterpolate(const Field< Type > &pf) const
Interpolate from points to faces.
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
fileName caseSystem() const
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
word findInstance(const fileName &dir, const word &name=word::null, const IOobject::readOption rOpt=IOobject::MUST_READ, const word &stopInstance=word::null) const
static word timeName(const scalar t, const int precision=precision_)
Ostream & printExecutionTime(OSstream &os) const
Print the elapsed ExecutionTime (cpu-time), ClockTime.
void size(const label n)
Older name for setAddressableSize.
label size() const noexcept
The number of elements in the list.
Templated 3D Vector derived from VectorSpace adding construction from 3 components,...
Base class for adjoint solvers.
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
T & ref()
Return reference to the managed object without nullptr checking.
void reset(autoPtr< T > &&other) noexcept
Delete managed object and set to new given pointer.
A bounding box defined in terms of min/max extrema points.
vector span() const
The bounding box span (from minimum to maximum)
Differentiation of the mesh data structure.
vectorField makeFaceCentresAndAreas_d(const pointField &p, const pointField &p_d)
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
dictionary & subDictOrAdd(const word &keyword, enum keyType::option matchOpt=keyType::REGEX)
Find and return a sub-dictionary for manipulation.
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
entry * add(entry *entryPtr, bool mergeEntry=false)
Add a new entry.
A special matrix type and solver, designed for finite area solutions of scalar equations....
void relax(const scalar alpha)
Relax matrix (for steady-state solution).
SolverPerformance< Type > solve(const dictionary &)
Solve returning the solution statistics.
Finite area mesh (used for 2-D non-Euclidian finite area method) defined using a patch of faces on a ...
const DimensionedField< scalar, areaMesh > & S() const
Return face areas.
static word meshSubDir
The mesh sub-directory name (usually "faMesh")
A face is a list of labels corresponding to mesh vertices.
pointField points(const UList< point > &pts) const
Return the points corresponding to this face.
virtual bool write()
Write the output fields.
Mesh data needed to do the Finite Volume discretisation.
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
const Time & time() const
Return the top-level database.
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
const volScalarField & pa() const
Return const reference to pressure.
const volVectorField & Ua() const
Return const reference to velocity.
Base class for incompressibleAdjoint solvers.
virtual void additionalSensitivityMapTerms(boundaryVectorField &sensitivityMap, const labelHashSet &patchIDs, const scalar dt)
const autoPtr< incompressibleAdjoint::adjointRASModel > & adjointTurbulence() const
Return const reference to the adjointRASModel.
const volVectorField & U() const
Return const reference to velocity.
const volScalarField & p() const
Return const reference to pressure.
const autoPtr< incompressible::RASModelVariables > & RASModelVariables() const
Return const reference to the turbulence model variables.
Solver of the adjoint to the eikonal PDE.
Solver of the adjoint to the Laplace grid displacement equation.
Abstract base class for adjoint-based sensitivities in incompressible flows.
virtual void clearSensitivities()
Zero sensitivity fields and their constituents.
incompressibleAdjointSolver & adjointSolver_
objectiveManager & objectiveManager_
incompressibleAdjointVars & adjointVars_
const incompressibleVars & primalVars_
Calculation of adjoint based sensitivities at wall faces.
bool includeTransposeStresses_
Include the transpose part of the adjoint stresses.
virtual void clearSensitivities()
Zero sensitivity fields and their constituents.
autoPtr< adjointEikonalSolver > & getAdjointEikonalSolver()
Get adjoint eikonal solver.
bool includeDistance_
Include distance variation in sens computation.
autoPtr< volVectorField > CfOnPatchPtr_
bool useSnGradInTranposeStresses_
Use snGrad in the transpose part of the adjoint stresses.
autoPtr< volVectorField > nfOnPatchPtr_
autoPtr< adjointEikonalSolver > eikonalSolver_
virtual void assembleSensitivities()
Assemble sensitivities.
void computeDerivativesSize()
Compute the number of faces on sensitivityPatchIDs_.
void setSuffixName()
Set suffix name for sensitivity fields.
autoPtr< volVectorField > SfOnPatchPtr_
bool includePressureTerm_
Include the adjoint pressure term in sens computation.
bool includeGradStressTerm_
Include the term containing the grad of the stress at the boundary.
bool includeSurfaceArea_
Include surface area in sens computation.
bool writeGeometricInfo_
Write geometric info for use by external programs.
void read()
Read controls and update solver pointers if necessary.
bool includeMeshMovement_
Include mesh movement variation in sens computation.
bool includeDivTerm_
Include the term from the deviatoric part of the stresses.
scalar computeRadius(const faMesh &aMesh)
virtual void accumulateIntegrand(const scalar dt)
Accumulate sensitivity integrands.
autoPtr< adjointMeshMovementSolver > meshMovementSolver_
bool smoothSensitivities_
bool includeObjective_
Include terms directly emerging from the objective function.
void smoothSensitivities()
PtrList< objective > & getObjectiveFunctions()
Return reference to objective functions.
A class representing the concept of 1 (one) that can be used to avoid manipulating objects known to b...
label nGeometricD() const
Return the number of valid geometric dimensions in the mesh.
virtual const faceList & faces() const
Return raw faces.
virtual const fileName & dbDir() const
Override the objectRegistry dbDir for a single-region case.
const boundBox & bounds() const
Return mesh bounding box.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
virtual const pointField & points() const
Return raw points.
const Vector< label > & geometricD() const
Return the vector of geometric directions in mesh.
label nPoints() const noexcept
Number of mesh points.
virtual bool write(const bool valid=true) const
Write using setting from DB.
const dictionary & dict() const
Return the construction dictionary.
void clearSensitivities()
Zero sensitivity fields and their constituents.
void setSuffix(const word &suffix)
Set suffix.
void write()
Write sensitivity fields.
A class for managing temporary objects.
void clear() const noexcept
const word & solverName() const
Return solver name.
Volume to surface and surface to volume mapping.
tmp< Field< Type > > mapToSurface(const GeometricBoundaryField< Type, fvPatchField, volMesh > &df) const
Map volume boundary field to surface.
void mapToVolume(const GeometricField< Type, faPatchField, areaMesh > &af, GeometricBoundaryField< Type, fvPatchField, volMesh > &bf) const
Map surface field to volume boundary field.
A class for handling words, derived from Foam::string.
Author Zeljko Tukovic, FMENA Hrvoje Jasak, Wikki Ltd.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
volSurfaceMapping vsm(aMesh)
Calculate the matrix for the laplacian of the field.
Calculate the finiteArea matrix for implicit and explicit sources.
const labelList nFaces(UPstream::listGatherValues< label >(aMesh.nFaces()))
#define DebugInfo
Report an information message using Foam::Info.
tmp< faMatrix< Type > > laplacian(const GeometricField< Type, faPatchField, areaMesh > &vf)
zeroField Sp(const Foam::zero, const GeometricField< Type, faPatchField, areaMesh > &)
A no-op source.
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh > > grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
const dimensionSet dimless
Dimensionless.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const dimensionSet dimArea(sqr(dimLength))
dimensionedSymmTensor twoSymm(const dimensionedSymmTensor &dt)
const dimensionSet dimVelocity
messageStream Info
Information stream (stdout output on master, null elsewhere)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Type gAverage(const FieldField< Field, Type > &f)
static constexpr const zero Zero
Global zero (0)
void unzipRow(const FieldField< Field, SymmTensor< Cmpt > > &input, const direction idx, FieldField< Field, Vector< Cmpt > > &result)
Extract a symmTensor field field row (x,y,z) == (0,1,2)
Type gMax(const FieldField< Field, Type > &f)
constexpr char nl
The newline '\n' character (0x0a)
loopControl iters(runTime, aMesh.solutionDict(), "solution")
#define forAll(list, i)
Loop across all elements in list.
A non-counting (dummy) refCount.