adjointBoundaryCondition.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2007-2020 PCOpt/NTUA
9  Copyright (C) 2013-2020 FOSS GP
10  Copyright (C) 2019-2020 OpenCFD Ltd.
11 -------------------------------------------------------------------------------
12 License
13  This file is part of OpenFOAM.
14 
15  OpenFOAM is free software: you can redistribute it and/or modify it
16  under the terms of the GNU General Public License as published by
17  the Free Software Foundation, either version 3 of the License, or
18  (at your option) any later version.
19 
20  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
21  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
22  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
23  for more details.
24 
25  You should have received a copy of the GNU General Public License
26  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
27 
28 \*---------------------------------------------------------------------------*/
29 
31 #include "emptyFvPatch.H"
32 #include "adjointSolverManager.H"
33 #include "HashTable.H"
35 #include "ATCUaGradU.H"
36 
37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38 
39 namespace Foam
40 {
41 
42 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
43 
44 template<class Type>
45 template<class Type2>
46 tmp
47 <
49 >
51 {
52  // Return field
53  typedef typename outerProduct<vector, Type2>::type GradType;
54  auto tresGrad = tmp<Field<GradType>>::New(patch_.size(), Zero);
55  auto& resGrad = tresGrad.ref();
56 
57  const labelList& faceCells = patch_.faceCells();
58  const fvMesh& mesh = patch_.boundaryMesh().mesh();
59  const cellList& cells = mesh.cells();
60 
61  // Go through the surfaceInterpolation scheme defined in gradSchemes for
62  // consistency
63  const GeometricField<Type2, fvPatchField, volMesh>& field =
65 
66  // Gives problems when grad(AdjointVar) is computed using a limited scheme,
67  // since it is not possible to know a priori how many words to expect in the
68  // stream.
69  // Interpolation scheme is now read through interpolation schemes.
70  /*
71  word gradSchemeName("grad(" + name + ')');
72  ITstream& is = mesh.gradScheme(gradSchemeName);
73  word schemeData(is);
74  */
75 
76  tmp<surfaceInterpolationScheme<Type2>> tinterpScheme
77  (
79  (
80  mesh,
81  mesh.interpolationScheme("interpolate(" + name + ")")
82  )
83  );
84 
85  GeometricField<Type2, fvsPatchField, surfaceMesh> surfField
86  (
87  tinterpScheme().interpolate(field)
88  );
89 
90  // Auxiliary fields
91  const surfaceVectorField& Sf = mesh.Sf();
92  tmp<vectorField> tnf = patch_.nf();
93  const vectorField& nf = tnf();
94  const scalarField& V = mesh.V();
95  const labelUList& owner = mesh.owner();
96 
97  // Compute grad value of cell adjacent to the boundary
98  forAll(faceCells, fI)
99  {
100  const label cI = faceCells[fI];
101  const cell& cellI = cells[cI];
102  for (const label faceI : cellI) // global face numbering
103  {
104  label patchID = mesh.boundaryMesh().whichPatch(faceI);
105  if (patchID == -1) //face is internal
106  {
107  const label own = owner[faceI];
108  tensor flux = Sf[faceI]*surfField[faceI];
109  if (cI == own)
110  {
111  resGrad[fI] += flux;
112  }
113  else
114  {
115  resGrad[fI] -= flux;
116  }
117  }
118  else // Face is boundary. Covers coupled patches as well
119  {
120  if (!isA<emptyFvPatch>(mesh.boundary()[patchID]))
121  {
122  const fvPatch& patchForFlux = mesh.boundary()[patchID];
123  const label boundaryFaceI = faceI - patchForFlux.start();
124  const vectorField& Sfb = Sf.boundaryField()[patchID];
125  resGrad[fI] +=
126  Sfb[boundaryFaceI]
127  *surfField.boundaryField()[patchID][boundaryFaceI];
128  }
129  }
130  }
131  resGrad[fI] /= V[cI];
132  }
133 
134  // This has concluded the computation of the grad at the cell next to the
135  // boundary. We now need to compute the grad at the boundary face
136  const fvPatchField<Type2>& bField = field.boundaryField()[patch_.index()];
137  resGrad = nf*bField.snGrad() + (resGrad - nf*(nf & resGrad));
138 
139  return tresGrad;
140 }
141 
142 
143 template<class Type>
145 {
146  if (!addATCUaGradUTerm_)
147  {
148  addATCUaGradUTerm_.reset(new bool(isA<ATCUaGradU>(getATC())));
149  }
150  return addATCUaGradUTerm_();
151 }
152 
153 
154 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
155 
156 template<class Type>
158 (
159  const adjointBoundaryCondition<Type>& adjointBC
160 )
161 :
162  patch_(adjointBC.patch_),
163  managerName_(adjointBC.managerName_),
164  adjointSolverName_(adjointBC.adjointSolverName_),
165  simulationType_(adjointBC.simulationType_),
166  boundaryContrPtr_
167  (
169  (
170  adjointBC.managerName_,
171  adjointBC.adjointSolverName_,
172  adjointBC.simulationType_,
173  adjointBC.patch_
174  )
175  ),
176  addATCUaGradUTerm_(adjointBC.addATCUaGradUTerm_)
177 {}
178 
179 
180 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
181 
182 template<class Type>
184 (
185  const fvPatch& p,
187  const word& solverName
188 )
189 :
190  patch_(p),
191  managerName_("objectiveManager" + solverName),
192  adjointSolverName_(solverName),
193  simulationType_("incompressible"),
194  boundaryContrPtr_(nullptr),
195  addATCUaGradUTerm_(nullptr)
196 {
197  // Set the boundaryContribution pointer
198  setBoundaryContributionPtr();
199 }
200 
201 
202 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
203 
204 template<class Type>
206 {
207  return managerName_;
208 }
209 
210 
211 template<class Type>
213 {
214  return adjointSolverName_;
215 }
216 
217 
218 template<class Type>
220 {
221  return simulationType_;
222 }
223 
224 
225 template<class Type>
227 {
228  // Note:
229  // Check whether there is an objectiveFunctionManager object in the registry
230  // Necessary for decomposePar if the libadjoint is loaded
231  // through controlDict. A nicer way should be found
232  const fvMesh& meshRef = patch_.boundaryMesh().mesh();
233  if (meshRef.foundObject<regIOobject>(managerName_))
234  {
235  boundaryContrPtr_.reset
236  (
238  (
239  managerName_,
240  adjointSolverName_,
241  simulationType_,
242  patch_
243  ).ptr()
244  );
245  }
246  else
247  {
249  << "No objectiveManager " << managerName_ << " available." << nl
250  << "Setting boundaryAdjointContributionPtr to nullptr. " << nl
251  << "OK for decomposePar."
252  << endl;
253  }
254 }
255 
256 
257 template<class Type>
260 {
261  return boundaryContrPtr_();
262 }
263 
264 
265 template<class Type>
267 {
268  return
269  patch_.boundaryMesh().mesh().template
270  lookupObject<ATCModel>("ATCModel" + adjointSolverName_);
271 }
272 
273 
274 template<class Type>
275 tmp
276 <
278 >
280 {
281  return
283  (
284  patch_.size(),
285  Zero
286  );
287 }
288 
289 
290 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
291 
292 } // End namespace Foam
293 
294 // ************************************************************************* //
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:67
Foam::adjointBoundaryCondition::getBoundaryAdjContribution
boundaryAdjointContribution & getBoundaryAdjContribution()
Get boundaryContribution.
Definition: adjointBoundaryCondition.C:259
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:52
Foam::surfaceInterpolationScheme::New
static tmp< surfaceInterpolationScheme< Type > > New(const fvMesh &mesh, Istream &schemeData)
Return new tmp interpolation scheme.
Definition: surfaceInterpolationScheme.C:40
p
volScalarField & p
Definition: createFieldRefs.H:8
adjointBoundaryCondition.H
Foam::fvc::flux
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
Foam::adjointBoundaryCondition::adjointBoundaryCondition
adjointBoundaryCondition(const fvPatch &p, const DimensionedField< Type, volMesh > &iF, const word &solverName)
Construct from field and base name.
Definition: adjointBoundaryCondition.C:184
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:65
HashTable.H
Foam::adjointBoundaryCondition::patch_
const fvPatch & patch_
Reference to patch.
Definition: adjointBoundaryCondition.H:59
Foam::adjointBoundaryCondition::simulationType
const word & simulationType() const
Return the simulationType.
Definition: adjointBoundaryCondition.C:219
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:61
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::outerProduct::type
typeOfRank< typename pTraits< arg1 >::cmptType, direction(pTraits< arg1 >::rank)+direction(pTraits< arg2 >::rank) >::type type
Definition: products.H:114
Foam::primitiveMesh::cells
const cellList & cells() const
Definition: primitiveMeshCells.C:138
Foam::adjointBoundaryCondition
Base class for solution control classes.
Definition: adjointBoundaryCondition.H:52
Foam::polyMesh::boundaryMesh
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:444
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::ATCModel
Base class for selecting the adjoint transpose convection model. Inherits from regIOobject to add loo...
Definition: ATCModel.H:60
Foam::objectRegistry::foundObject
bool foundObject(const word &name, const bool recursive=false) const
Is the named Type found?
Definition: objectRegistryTemplates.C:379
Foam::adjointBoundaryCondition::getATC
const ATCModel & getATC() const
ATC type might be useful for a number of BCs. Return here.
Definition: adjointBoundaryCondition.C:266
adjointSolverManager.H
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::vectorField
Field< vector > vectorField
Specialisation of Field<T> for vector.
Definition: primitiveFieldsFwd.H:54
Foam::adjointBoundaryCondition::setBoundaryContributionPtr
void setBoundaryContributionPtr()
Set the ptr to the correct boundaryAdjointContribution.
Definition: adjointBoundaryCondition.C:226
Foam::interpolate
bool interpolate(const vector &p1, const vector &p2, const vector &o, vector &n, scalar l)
Definition: curveTools.C:75
Foam::polyBoundaryMesh::mesh
const polyMesh & mesh() const noexcept
Return the mesh reference.
Definition: polyBoundaryMesh.H:152
ATCUaGradU.H
Foam::Field
Generic templated field type.
Definition: Field.H:63
Foam::schemesLookup::interpolationScheme
ITstream & interpolationScheme(const word &name) const
Get interpolation scheme for given name, or default.
Definition: schemesLookup.C:199
Foam::adjointBoundaryCondition::adjointSolverName_
word adjointSolverName_
adjointSolver name corresponding to field
Definition: adjointBoundaryCondition.H:65
patchID
label patchID
Definition: boundaryProcessorFaPatchPoints.H:5
Foam::fvPatch
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
Definition: fvPatch.H:65
Foam::cellList
List< cell > cellList
A List of cells.
Definition: cellListFwd.H:47
Foam::polyBoundaryMesh::whichPatch
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
Definition: polyBoundaryMesh.C:812
field
rDeltaTY field()
Foam::objectRegistry::lookupObject
const Type & lookupObject(const word &name, const bool recursive=false) const
Definition: objectRegistryTemplates.C:434
Foam::adjointBoundaryCondition::simulationType_
word simulationType_
simulationType corresponding to field.
Definition: adjointBoundaryCondition.H:69
Foam::adjointBoundaryCondition::addATCUaGradUTerm
bool addATCUaGradUTerm()
Whether to add the extra term from the UaGradU formulation.
Definition: adjointBoundaryCondition.C:144
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::adjointBoundaryCondition::computePatchGrad
tmp< Field< typename Foam::outerProduct< Foam::vector, Type2 >::type > > computePatchGrad(word name)
Get gradient of field on a specific boundary.
Foam::adjointBoundaryCondition::managerName_
word managerName_
objectiveManager name corresponding to field
Definition: adjointBoundaryCondition.H:62
Foam::fvMesh
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:85
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
emptyFvPatch.H
Foam::fvMesh::owner
const labelUList & owner() const
Internal face owner. Note bypassing virtual mechanism so.
Definition: fvMesh.H:407
Foam::volVectorField
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:62
Foam::New
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
Global function forwards to reuseTmpDimensionedField::New.
Definition: DimensionedFieldReuseFunctions.H:105
Foam::fvMesh::boundary
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:685
Foam::regIOobject
regIOobject is an abstract class derived from IOobject to handle automatic object registration with t...
Definition: regIOobject.H:73
Foam::boundaryAdjointContribution
Abstract base class for computing contributions of the objective functions to the adjoint boundary co...
Definition: boundaryAdjointContribution.H:58
Foam::adjointBoundaryCondition::dxdbMult
virtual tmp< Field< typename Foam::outerProduct< Foam::vector, Type >::type > > dxdbMult() const
Return contribution to sensitivity derivatives.
Foam::nl
constexpr char nl
Definition: Ostream.H:404
Foam::adjointBoundaryCondition::objectiveManagerName
const word & objectiveManagerName() const
Return objectiveManager name.
Definition: adjointBoundaryCondition.C:205
Foam::primitiveMesh::reset
void reset(const label nPoints, const label nInternalFaces, const label nFaces, const label nCells)
Reset this primitiveMesh given the primitive array sizes.
Definition: primitiveMesh.C:207
Foam::roots::type
type
Types of root.
Definition: Roots.H:54
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
cells
const cellShapeList & cells
Definition: gmvOutputHeader.H:3
Foam::adjointBoundaryCondition::adjointSolverName
const word & adjointSolverName() const
Return adjointSolverName.
Definition: adjointBoundaryCondition.C:212
Foam::labelUList
UList< label > labelUList
A UList of labels.
Definition: UList.H:85
Foam::surfaceVectorField
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
Definition: surfaceFieldsFwd.H:59
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328
Foam::tensor
Tensor< scalar > tensor
Tensor of scalars, i.e. Tensor<scalar>.
Definition: symmTensor.H:61
surfaceInterpolationScheme.H
Foam::boundaryAdjointContribution::New
static autoPtr< boundaryAdjointContribution > New(const word &managerName, const word &adjointSolverName, const word &simulationType, const fvPatch &patch)
Return a reference to the selected turbulence model.
Definition: boundaryAdjointContribution.C:59
Foam::DimensionedField
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
Definition: DimensionedField.H:54
Foam::fvMesh::V
const DimensionedField< scalar, volMesh > & V() const
Return cell volumes.
Definition: fvMeshGeometry.C:179
Foam::adjointBoundaryCondition::addATCUaGradUTerm_
autoPtr< bool > addATCUaGradUTerm_
Whether to add the extra term from the UaGradU formulation.
Definition: adjointBoundaryCondition.H:79
Foam::fvMesh::Sf
const surfaceVectorField & Sf() const
Return cell face area vectors.
Definition: fvMeshGeometry.C:319