adjointEikonalSolverIncompressible.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 "wallFvPatch.H"
32 #include "patchDistMethod.H"
33 
34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38 
39 namespace incompressible
40 {
41 
42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
43 
45 
46 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
47 
49 {
50  wordList daTypes
51  (
52  mesh_.boundary().size(),
53  fixedValueFvPatchScalarField::typeName
54  );
55 
56  for (const label patchi : wallPatchIDs_)
57  {
58  daTypes[patchi] = zeroGradientFvPatchScalarField::typeName;
59  }
60 
61  return daTypes;
62 }
63 
64 
66 {
67  nEikonalIters_ = dict_.getOrDefault<label>("iters", 1000);
68  tolerance_ = dict_.getOrDefault<scalar>("tolerance", 1e-6);
69  epsilon_ = dict_.getOrDefault<scalar>("epsilon", 0.1);
70 }
71 
72 
74 {
75  // Primal distance field
76  const volScalarField& d = RASModelVars_().d();
77 
79  (
80  IOobject
81  (
82  "ny",
83  mesh_.time().timeName(),
84  mesh_,
87  false
88  ),
89  mesh_,
91  patchDistMethod::patchTypes<vector>(mesh_, wallPatchIDs_)
92  );
93 
94  const fvPatchList& patches = mesh_.boundary();
95  volVectorField::Boundary& nybf = ny.boundaryFieldRef();
96 
97  for (const label patchi : wallPatchIDs_)
98  {
99  nybf[patchi] == -patches[patchi].nf();
100  }
101 
102  ny = fvc::grad(d);
103 
105 
106  return tmp<surfaceScalarField>::New("yPhi", mesh_.Sf() & nf);
107 }
108 
109 
110 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
111 
112 adjointEikonalSolver::adjointEikonalSolver
113 (
114  const fvMesh& mesh,
115  const dictionary& dict,
116  const autoPtr<incompressible::RASModelVariables>& RASModelVars,
117  incompressibleAdjointVars& adjointVars,
118  const labelHashSet& sensitivityPatchIDs
119 )
120 :
121  mesh_(mesh),
122  dict_(dict.subOrEmptyDict("adjointEikonalSolver")),
123  RASModelVars_(RASModelVars),
124  adjointTurbulence_(adjointVars.adjointTurbulence()),
125  sensitivityPatchIDs_(sensitivityPatchIDs),
126  nEikonalIters_(-1),
127  tolerance_(-1),
128  epsilon_(Zero),
129  wallPatchIDs_(mesh_.boundaryMesh().findPatchIDs<wallPolyPatch>()),
130  da_
131  (
132  IOobject
133  (
134  word
135  (
136  adjointVars.useSolverNameForFields() ?
137  "da" + adjointTurbulence_().adjointSolverName() :
138  "da"
139  ),
140  mesh_.time().timeName(),
141  mesh_,
144  ),
145  mesh_,
147  patchTypes()
148  ),
149  source_
150  (
151  IOobject
152  (
153  "sourceEikonal",
154  mesh_.time().timeName(),
155  mesh_,
158  ),
159  mesh_,
161  ),
162  distanceSensPtr_(createZeroBoundaryPtr<vector>(mesh_))
163 {
164  read();
165 }
166 
167 
168 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
169 
171 {
172  dict_ = dict.subOrEmptyDict("adjointEikonalSolver");
173 
174  return true;
175 }
176 
177 
179 {
180  // Accumulate integrand from the current time step
182 }
183 
184 
186 {
187  read();
188 
189  // Primal distance field
190  const volScalarField& d = RASModelVars_().d();
191 
192  // Convecting flux
194  const surfaceScalarField& yPhi = tyPhi();
195 
196  // Iterate the adjoint to the eikonal equation
197  for (label iter = 0; iter < nEikonalIters_; ++iter)
198  {
199  read();
200 
201  Info<< "Adjoint Eikonal Iteration : " << iter << endl;
202 
203  fvScalarMatrix daEqn
204  (
205  2*fvm::div(-yPhi, da_)
208  + source_
209  );
210 
211  daEqn.relax();
212  scalar residual = daEqn.solve().initialResidual();
213  Info<< "Max da " << gMax(mag(da_)()) << endl;
214 
216 
217  // Check convergence
218  if (residual < tolerance_)
219  {
220  Info<< "\n***Reached adjoint eikonal convergence limit, iteration "
221  << iter << "***\n\n";
222  break;
223  }
224  }
225  if (debug)
226  {
227  da_.write();
228  }
229 }
230 
231 
233 {
234  source_ == dimensionedScalar(source_.dimensions(), Zero);
236 }
237 
238 
240 {
241  Info<< "Calculating distance sensitivities " << endl;
242 
243  boundaryVectorField& distanceSens = distanceSensPtr_();
244 
245  const volScalarField& d = RASModelVars_().d();
246  for (const label patchi : sensitivityPatchIDs_)
247  {
248  vectorField nf(mesh_.boundary()[patchi].nf());
249 
250  // No surface area included. Will be done by the actual sensitivity tool
251  distanceSens[patchi] =
252  -2.*da_.boundaryField()[patchi]
253  *d.boundaryField()[patchi].snGrad()
254  *d.boundaryField()[patchi].snGrad()*nf;
255  }
256  return distanceSens;
257 }
258 
259 
261 {
262  Info<< "Calculating distance sensitivities " << endl;
263 
264  const volScalarField& d = RASModelVars_().d();
265  const volVectorField gradD(fvc::grad(d));
266 
267  volVectorField gradDDa
268  (
269  IOobject
270  (
271  "gradDDa",
272  mesh_.time().timeName(),
273  mesh_,
276  ),
277  mesh_,
278  dimensionedVector(d.dimensions()*da_.dimensions()/dimLength, Zero),
279  patchDistMethod::patchTypes<vector>(mesh_, wallPatchIDs_)
280  );
281  gradDDa = fvc::grad(d*da_);
282 
283  tmp<volTensorField> tdistanceSens
284  (
285  new volTensorField
286  (
287  IOobject
288  (
289  "distanceSensFI",
290  mesh_.time().timeName(),
291  mesh_,
294  ),
295  mesh_,
296  dimensionedTensor(da_.dimensions(), Zero)
297  )
298  );
299  volTensorField& distanceSens = tdistanceSens.ref();
300 
301  distanceSens =
302  - 2.*da_*gradD*gradD
303  - epsilon_*gradD*gradDDa
304  + epsilon_*da_*d*fvc::grad(gradD);
305 
306  return tdistanceSens;
307 }
308 
309 
311 {
312  return da_;
313 }
314 
315 
317 {
318  const volScalarField& d = RASModelVars_().d();
319  volVectorField gradD(fvc::grad(d));
320  return tmp<volVectorField>::New("gradEikonal", 2*gradD & fvc::grad(gradD));
321 }
322 
323 
324 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
325 
326 } // End namespace incompressible
327 } // End namespace Foam
328 
329 // ************************************************************************* //
Foam::variablesSet::useSolverNameForFields
bool useSolverNameForFields() const
Append solver name to fields?
Definition: variablesSet.C:90
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::IOobject::NO_WRITE
Definition: IOobject.H:195
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:169
Foam::incompressible::adjointEikonalSolver::nEikonalIters_
label nEikonalIters_
Definition: adjointEikonalSolverIncompressible.H:174
Foam::incompressible::adjointEikonalSolver::reset
void reset()
Reset source term.
Definition: adjointEikonalSolverIncompressible.C:232
Foam::IOobject::AUTO_WRITE
Definition: IOobject.H:194
Foam::incompressible::adjointEikonalSolver::getFISensitivityTerm
tmp< volTensorField > getFISensitivityTerm() const
Return the volume-based sensitivity term depending on da.
Definition: adjointEikonalSolverIncompressible.C:260
Foam::incompressible::adjointEikonalSolver::accumulateIntegrand
void accumulateIntegrand(const scalar dt)
Accumulate source term.
Definition: adjointEikonalSolverIncompressible.C:178
Foam::incompressible::adjointEikonalSolver::distanceSensitivities
boundaryVectorField & distanceSensitivities()
Return the sensitivity term depending on da.
Definition: adjointEikonalSolverIncompressible.C:239
Foam::fvc::grad
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:54
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:65
Foam::incompressible::defineTypeNameAndDebug
defineTypeNameAndDebug(adjointEikonalSolver, 0)
Foam::dimLength
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:52
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::incompressible::adjointEikonalSolver::wallPatchIDs_
labelHashSet wallPatchIDs_
Definition: adjointEikonalSolverIncompressible.H:180
Foam::read
bool read(const char *buf, int32_t &val)
Same as readInt32.
Definition: int32.H:108
Foam::incompressible::adjointEikonalSolver::solve
void solve()
Calculate the adjoint distance field.
Definition: adjointEikonalSolverIncompressible.C:185
Foam::Time::timeName
static word timeName(const scalar t, const int precision=precision_)
Definition: Time.C:780
wallFvPatch.H
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::incompressible::adjointEikonalSolver::da_
volScalarField da_
Definition: adjointEikonalSolverIncompressible.H:182
Foam::HashSet< label, Hash< label > >
Foam::incompressible::adjointEikonalSolver::da
const volScalarField & da()
Return the adjoint distance field.
Definition: adjointEikonalSolverIncompressible.C:310
Foam::incompressibleAdjointVars
Class including all adjoint fields for incompressible flows.
Definition: incompressibleAdjointVars.H:52
patchDistMethod.H
Foam::incompressible::adjointEikonalSolver::sensitivityPatchIDs_
const labelHashSet & sensitivityPatchIDs_
Definition: adjointEikonalSolverIncompressible.H:172
patchTypes
wordList patchTypes(nPatches)
Foam::dimensionedVector
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Definition: dimensionedVector.H:50
Foam::incompressible::adjointEikonalSolver
Solver of the adjoint to the eikonal PDE.
Definition: adjointEikonalSolverIncompressible.H:146
Foam::fvc::laplacian
tmp< GeometricField< Type, fvPatchField, volMesh > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcLaplacian.C:47
Foam::Time::printExecutionTime
Ostream & printExecutionTime(OSstream &os) const
Print the elapsed ExecutionTime (cpu-time), ClockTime.
Definition: TimeIO.C:618
Foam::dimTime
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:53
Foam::fvm::SuSp
tmp< fvMatrix< Type > > SuSp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
Foam::incompressible::adjointEikonalSolver::patchTypes
wordList patchTypes() const
Return the boundary condition types for da.
Definition: adjointEikonalSolverIncompressible.C:48
Foam::Field< vector >
Foam::pow3
dimensionedScalar pow3(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:89
Foam::incompressible::adjointEikonalSolver::dict_
dictionary dict_
Definition: adjointEikonalSolverIncompressible.H:165
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::fvm::laplacian
tmp< fvMatrix< Type > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmLaplacian.C:48
Foam::fvMatrix::relax
void relax(const scalar alpha)
Relax matrix (for steady-state solution).
Definition: fvMatrix.C:1183
Foam::incompressible::adjointEikonalSolver::tolerance_
scalar tolerance_
Definition: adjointEikonalSolverIncompressible.H:176
Foam::IOobject::READ_IF_PRESENT
Definition: IOobject.H:187
Foam::incompressible::adjointEikonalSolver::d
const volScalarField & d()
Return the distance field.
Foam::incompressible::adjointEikonalSolver::read
void read()
Read options each time a new solution is found.
Definition: adjointEikonalSolverIncompressible.C:65
Foam::dimensionedScalar
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Definition: dimensionedScalarFwd.H:42
Foam::PtrList
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: List.H:59
Foam::incompressible::adjointEikonalSolver::computeYPhi
tmp< surfaceScalarField > computeYPhi()
Compute convecting velocity.
Definition: adjointEikonalSolverIncompressible.C:73
Foam::incompressible::adjointEikonalSolver::source_
volScalarField source_
Definition: adjointEikonalSolverIncompressible.H:184
Foam::incompressible::adjointEikonalSolver::gradEikonal
tmp< volVectorField > gradEikonal()
Return the gradient of the eikonal equation.
Definition: adjointEikonalSolverIncompressible.C:316
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::incompressibleAdjoint::adjointRASModel::distanceSensitivities
virtual tmp< volScalarField > distanceSensitivities()=0
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::wallPolyPatch
Foam::wallPolyPatch.
Definition: wallPolyPatch.H:50
Foam::incompressible::adjointEikonalSolver::adjointTurbulence_
autoPtr< Foam::incompressibleAdjoint::adjointRASModel > & adjointTurbulence_
Definition: adjointEikonalSolverIncompressible.H:170
Foam::fvMesh
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:85
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::fvc::interpolate
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
Foam::incompressible::adjointEikonalSolver::readDict
virtual bool readDict(const dictionary &dict)
Read dict if changed.
Definition: adjointEikonalSolverIncompressible.C:170
Foam::fvMesh::boundary
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:685
Foam::incompressible::adjointEikonalSolver::epsilon_
scalar epsilon_
Definition: adjointEikonalSolverIncompressible.H:178
Foam::autoPtr
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: HashPtrTable.H:53
Foam::incompressible::adjointEikonalSolver::mesh_
const fvMesh & mesh_
Definition: adjointEikonalSolverIncompressible.H:163
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::GeometricField::Boundary
The boundary fields.
Definition: GeometricField.H:115
adjointEikonalSolverIncompressible.H
Foam::GeometricField::ref
Internal & ref(const bool updateAccessTime=true)
Return a reference to the dimensioned internal field.
Definition: GeometricField.C:749
Foam::incompressibleAdjointVars::adjointTurbulence
const autoPtr< incompressibleAdjoint::adjointRASModel > & adjointTurbulence() const
Return const reference to the adjointRASModel.
Definition: incompressibleAdjointVars.C:70
Foam::List< word >
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::tmp::New
static tmp< T > New(Args &&... args)
Construct tmp of T with forwarding arguments.
Foam::incompressible::adjointEikonalSolver::RASModelVars_
const autoPtr< incompressible::RASModelVariables > & RASModelVars_
Definition: adjointEikonalSolverIncompressible.H:167
patches
const polyBoundaryMesh & patches
Definition: convertProcessorPatches.H:65
Foam::fvMatrix
A special matrix type and solver, designed for finite volume solutions of scalar equations....
Definition: fvPatchField.H:68
Foam::fvMatrix::solve
SolverPerformance< Type > solve(const dictionary &)
Solve returning the solution statistics.
Definition: fvMatrixSolve.C:319
Foam::fvMesh::time
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:280
Foam::VectorSpace< Vector< scalar >, scalar, 3 >::zero
static const Vector< scalar > zero
Definition: VectorSpace.H:115
Foam::incompressible::adjointEikonalSolver::distanceSensPtr_
autoPtr< boundaryVectorField > distanceSensPtr_
Wall face sens w.r.t. (x,y.z)
Definition: adjointEikonalSolverIncompressible.H:187
Foam::dictionary::getOrDefault
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:148
Foam::GeometricField< scalar, fvPatchField, volMesh >
Foam::IOobject::NO_READ
Definition: IOobject.H:188
Foam::fvm::div
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:48
Foam::dimensionedTensor
dimensioned< tensor > dimensionedTensor
Dimensioned tensor obtained from generic dimensioned type.
Definition: dimensionedTensor.H:52
Foam::gMax
Type gMax(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:592
Foam::GeometricField::boundaryField
const Boundary & boundaryField() const
Return const-reference to the boundary field.
Definition: GeometricFieldI.H:62
Foam::dimless
const dimensionSet dimless
Dimensionless.
Definition: dimensionSets.C:189
Foam::fvMesh::Sf
const surfaceVectorField & Sf() const
Return cell face area vectors.
Definition: fvMeshGeometry.C:319