adjointkOmegaSST.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) 2014-2022 PCOpt/NTUA
9 Copyright (C) 2014-2022 FOSS GP
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "adjointkOmegaSST.H"
30#include "wallFvPatch.H"
33#include "linear.H"
34#include "reverseLinear.H"
38
39// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
40
41namespace Foam
42{
43namespace incompressibleAdjoint
44{
45namespace adjointRASModels
46{
47
48// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
49
52
53
54// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
55
57{
59 (
60 min
61 (
62 max
63 (
64 (scalar(1)/betaStar_)*sqrt(k())/(omega()*y_),
65 scalar(500)*nu()/(sqr(y_)*omega())
66 ),
68 ),
69 scalar(10)
70 );
71
72 return tanh(pow4(arg1));
73}
74
75
77{
79 (
80 max
81 (
82 (scalar(2)/betaStar_)*sqrt(k())/(omega()*y_),
83 scalar(500)*nu()/(sqr(y_)*omega())
84 ),
85 scalar(100)
86 );
87
88 return tanh(sqr(arg2));
89}
90
91
93(
94 const volScalarField& GbyNu0,
95 const volScalarField& F2,
96 const volScalarField& S2
97) const
98{
99 return min
100 (
101 GbyNu0,
103 *max(a1_*omega(), b1_*F2*sqrt(S2))
104 );
105}
106
107
109(
110 const volScalarField::Internal& GbyNu0,
113) const
114{
115 return min
116 (
117 GbyNu0,
118 (c1_/a1_)*betaStar_*omega()()
119 *max(a1_*omega()(), b1_*F2*sqrt(S2))
120 );
121}
122
123
125{
126 auto tzeroFirstCell =
128 (
130 (
131 "zeroFirstCell",
133 mesh_,
136 ),
137 mesh_,
139 );
140 volScalarField& zeroFirstCell = tzeroFirstCell.ref();
141
143 label counter(0);
144 for (const fvPatchScalarField& omegab : omega().boundaryField())
145 {
146 const fvPatch& patch = omegab.patch();
147 if (isA<omegaWallFunctionFvPatchScalarField>(omegab))
148 {
149 const label patchi = patch.index();
150 const labelList& faceCells = patch.faceCells();
151 fvPatchScalarField& bf = zeroFirstCell.boundaryFieldRef()[patchi];
152 forAll(faceCells, faceI)
153 {
154 const label cellI = faceCells[faceI];
155
156 zeroFirstCell[cellI] = 0.;
157 bf[faceI] = 0.;
158 firstCellIDs_[counter++] = cellI;
159 }
160 }
161 }
162 firstCellIDs_.setSize(counter);
163
164 zeroFirstCell.correctBoundaryConditions();
165
166 return tzeroFirstCell;
167}
168
169
171{
172 const volVectorField& U = primalVars_.U();
173 const volVectorField& Ua = adjointVars_.UaInst();
174 word scheme("coeffsDiff");
175 tmp<volScalarField> tdR_dnut
176 (
177 dNutdbMult(U, Ua, nutRef(), scheme)
178 + dNutdbMult(k(), ka(), alphaK_, nutRef(), scheme)
181 );
182 volScalarField& dRdnut = tdR_dnut.ref();
183
184 forAll(mesh_.boundary(), pI)
185 {
186 const fvPatch& patch = mesh_.boundary()[pI];
187 const fvPatchScalarField& nutb = nutRef().boundaryField()[pI];
188 const labelList& faceCells = patch.faceCells();
189 if (isA<nutkWallFunctionFvPatchScalarField>(nutb))
190 {
191 fvPatchScalarField& bf = dRdnut.boundaryFieldRef()[pI];
192 const scalarField kSnGrad(k().boundaryField()[pI].snGrad());
193 const scalarField omegaSnGrad
194 (
195 omega().boundaryField()[pI].snGrad()
196 );
197 const fvPatchScalarField& akb = alphaK_.boundaryField()[pI];
198 const fvPatchScalarField& aOmegab = alphaOmega_.boundaryField()[pI];
199 const vectorField USnGrad(U.boundaryField()[pI].snGrad());
200 const fvPatchTensorField& gradUb = gradU_.boundaryField()[pI];
201 const vectorField nf(mesh_.boundary()[pI].nf());
202 forAll(faceCells, fI)
203 {
204 const label cI(faceCells[fI]);
205 bf[fI] =
206 - wa()[cI]*zeroFirstCell_[cI]*aOmegab[fI]*omegaSnGrad[fI]
207 - ka()[cI]*akb[fI]*kSnGrad[fI]
208 - (Ua[cI] & (USnGrad[fI] + (dev2(gradUb[fI]) & nf[fI])));
209 }
210 }
211 }
212
213 return tdR_dnut;
214}
215
216
218{
219 return
220 (
221 - case_1_nut_*k()/sqr(omega())
222 - a1_*k()/(b1_*S_*F2_*F2_)*dF2_domega()
223 );
224}
225
226
228{
229 return
230 (
231 a1_/max(a1_*omega(), b1_*F2_*S_)
232 - a1_*k()/(b1_*S_*F2_*F2_)*dF2_dk()
233 );
234}
235
236
238{
240 (
241 max
242 (
243 (scalar(2)/betaStar_)*sqrt(k())/(omega()*y_),
244 scalar(500)*nu()/(sqr(y_)*omega())
245 ),
246 scalar(100)
247 );
248
249 return
250 - scalar(2)*arg2*(1 - F2_*F2_)*
251 (
252 case_2_nut_*scalar(2)*sqrt(k())/(betaStar_*sqr(omega())*y_)
253 + case_3_nut_*scalar(500)*nu()/sqr(omega()*y_)
254 );
255}
256
257
259{
261 (
262 max
263 (
264 (scalar(2)/betaStar_)*sqrt(k())/(omega()*y_),
265 scalar(500)*nu()/(sqr(y_)*omega())
266 ),
267 scalar(100)
268 );
269
270 return
271 case_2_nut_*scalar(2)*arg2*(1 - F2_*F2_)/(betaStar_*omega()*y_
272 *sqrt(k()));
273}
274
275
277{
278 return
279 (
281 *(
282 max(a1_*omega(), b1_*F2_*S_)
284 + (scalar(1) - case_1_nut_)*omega()*b1_*S_*dF2_domega()
285 )
286 );
287}
288
289
291{
293}
294
295
297{
298 const volVectorField& U = primalVars_.U();
302 word scheme("coeffsDiff");
303
305 (
306 nutRef()
307 *(
308 coeffsDifferentiation(k(), ka(), scheme)*(alphaK1_ - alphaK2_)
311 )
312 );
313 volScalarField& dRdF1 = tdRdF1.ref();
314
315 typedef fixedValueFvPatchScalarField fixedValue;
316 typedef zeroGradientFvPatchScalarField zeroGrad;
318 forAll(mesh_.boundary(), pI)
319 {
320 const fvPatchScalarField& kb = k().boundaryField()[pI];
321 const fvPatchScalarField& omegab = omega().boundaryField()[pI];
322 fvPatchScalarField& dRdF1b = dRdF1.boundaryFieldRef()[pI];
323 if
324 (
325 isA<zeroGrad>(kb)
326 && (isA<zeroGrad>(omegab) || isA<omegaWF>(omegab))
327 )
328 {
329 dRdF1b = dRdF1b.patchInternalField();
330 }
331 else if (isA<fixedValue>(kb) && isA<fixedValue>(omegab)
332 )
333 {
334 // Note: might need revisiting
335 dRdF1b = dRdF1b.patchInternalField();
336 }
337 }
338
339 volScalarField dR_dF1_coeffs
340 (
342 *(
343 (gamma1_ - gamma2_)*((2.0/3.0)*omega()*tdivU - tGPrime)
344 + omega()*omega()*(beta1_ - beta2_) + CDkOmega_
345 )
346 );
347
348 forAll(mesh_.boundary(), pI)
349 {
350 const fvPatchScalarField& kb = k().boundaryField()[pI];
351 const fvPatchScalarField& omegab = omega().boundaryField()[pI];
352 fvPatchScalarField& dRdF1b = dR_dF1_coeffs.boundaryFieldRef()[pI];
353 if
354 (
355 isA<zeroGrad>(kb)
356 && (isA<zeroGrad>(omegab) || isA<omegaWF>(omegab))
357 )
358 {
359 dRdF1b = Zero;
360 }
361 else if (isA<fixedValue>(kb) && isA<fixedValue>(omegab))
362 {
363 dRdF1b = dRdF1b.patchInternalField();
364 }
365 }
366
367 dRdF1 += dR_dF1_coeffs;
368 return tdRdF1;
369}
370
371
373(
374 const volScalarField& arg1
375) const
376{
377 return
378 (
379 4*pow3(arg1)*(scalar(1) - F1_*F1_)
380 *(
382 - case_2_F1_*scalar(500)*nu()/sqr(omega()*y_)
385 )
386 );
387}
388
389
391(
392 const volScalarField& arg1
393) const
394{
395 return
396 (
397 - case_3_F1_*scalar(4)*pow3(arg1)*(scalar(1) - F1_*F1_)
398 *scalar(8)*k()*sqr(alphaOmega2_/(CDkOmegaPlus_*y_))/omega()*gradK_
399 );
400}
401
402
404{
406 (
407 min
408 (
409 max
410 (
411 (scalar(1)/betaStar_)*sqrt(k())/(omega()*y_),
412 scalar(500)*nu()/(sqr(y_)*omega())
413 ),
415 ),
416 scalar(10)
417 );
418
422
423 surfaceScalarField dR_dGradOmega
424 (
425 interpolationScheme<vector>("div(dR_dGradOmega)")().
427 );
428
429 return
430 (
432 - fvc::div(dR_dGradOmega)
433 );
434}
435
436
438{
439 tmp<volVectorField> tsource
440 (
442 );
443 volVectorField& source = tsource.ref();
444
445 forAll(mesh_.boundary(), pI)
446 {
447 const fvPatchScalarField& omegab = omega().boundaryField()[pI];
448 fvPatchVectorField& sourceb = source.boundaryFieldRef()[pI];
449 if
450 (
451 isA<zeroGradientFvPatchScalarField>(omegab)
452 || isA<omegaWallFunctionFvPatchScalarField>(omegab)
453 )
454 {
455 sourceb = Zero;
456 }
457 else if (isA<fixedValueFvPatchScalarField>(omegab))
458 {
459 sourceb = sourceb.patchInternalField();
460 }
461 }
462
463 surfaceScalarField sourceFaces
464 (
465 interpolationScheme<vector>("sourceAdjontkOmegaSST")().
466 interpolate(source) & mesh_.Sf()
467 );
468
469 return
470 (
471 // Differentiation of omega in CDkOmega
473 // Differentiation of grad(omega) in CDkOmega
474 + fvc::div(sourceFaces)
475 );
476}
477
478
480{
481 tmp<volVectorField> tsource
482 (
484 );
485 volVectorField& source = tsource.ref();
486
487 forAll(mesh_.boundary(), pI)
488 {
489 const fvPatchScalarField& kb = k().boundaryField()[pI];
490 fvPatchVectorField& sourceb = source.boundaryFieldRef()[pI];
491 if (isA<zeroGradientFvPatchScalarField>(kb))
492 {
493 sourceb = Zero;
494 }
495 else if (isA<fixedValueFvPatchScalarField>(kb))
496 {
497 sourceb = sourceb.patchInternalField();
498 }
499 }
500
501 return
503 (
504 interpolationScheme<vector>("sourceAdjontkOmegaSST")().
505 interpolate(source) & mesh_.Sf()
506 );
507}
508
509
511(
512 const volScalarField& arg1
513) const
514{
515 return
516 (
517 4*pow3(arg1)*(scalar(1) - F1_*F1_)
518 *(
519 case_1_F1_*0.5/(betaStar_*omega()*sqrt(k())*y_)
521 )
522 );
523}
524
525
527(
528 const volScalarField& arg1
529) const
530{
531 return
532 (
533 - case_3_F1_*scalar(4)*pow3(arg1)*(scalar(1) - F1_*F1_)
534 *scalar(8)*k()*sqr(alphaOmega2_/(CDkOmegaPlus_*y_))/omega()
536 );
537}
538
539
541{
543 (
544 min
545 (
546 max
547 (
548 (scalar(1)/betaStar_)*sqrt(k())/(omega()*y_),
549 scalar(500)*nu()/(sqr(y_)*omega())
550 ),
552 ),
553 scalar(10)
554 );
555
557 volScalarField dF1_dk(this->dF1_dk(arg1));
559
560 surfaceScalarField dR_dGradK
561 (
562 interpolationScheme<vector>("div(dR_dGradK)")().
564 );
565
566 return
567 (
569 - fvc::div(dR_dGradK)
570 );
571}
572
573
575(
576 const volScalarField& primalField,
577 const volScalarField& adjointField,
578 const word& schemeName
579) const
580{
582 (interpolationScheme<scalar>(schemeName));
583 surfaceScalarField snGradPrimal(fvc::snGrad(primalField)*mesh_.magSf());
584 surfaceScalarField flux(scheme().interpolate(adjointField)*snGradPrimal);
585
586 forAll(mesh_.boundary(), pI)
587 {
588 const fvPatchScalarField& primalbf = primalField.boundaryField()[pI];
589 if (isA<zeroGradientFvPatchScalarField>(primalbf))
590 {
591 // Needless, but just to be safe
592 snGradPrimal.boundaryFieldRef()[pI] = Zero;
593 flux.boundaryFieldRef()[pI] = Zero;
594 }
595 else if (isA<fixedValueFvPatchScalarField>(primalbf))
596 {
597 // Note: to be potentially revisited
598 snGradPrimal.boundaryFieldRef()[pI] = Zero;
599 flux.boundaryFieldRef()[pI] = Zero;
600 }
601 }
602
603 return (fvc::div(flux) - adjointField*fvc::div(snGradPrimal));
604}
605
606
608(
609 const volScalarField& primalField,
610 const volScalarField& adjointField,
611 const volScalarField& coeffField,
612 const volScalarField& bcField,
613 const word& schemeName
614) const
615{
617 (interpolationScheme<scalar>(schemeName));
618 surfaceScalarField snGradPrimal(fvc::snGrad(primalField)*mesh_.magSf());
619 surfaceScalarField flux(scheme().interpolate(adjointField)*snGradPrimal);
620
621 forAll(mesh_.boundary(), pI)
622 {
623 const fvPatchScalarField& bc = bcField.boundaryField()[pI];
624 if (isA<zeroGradientFvPatchScalarField>(bc))
625 {
626 const fvPatchScalarField& coeffFieldb =
627 coeffField.boundaryField()[pI];
628 snGradPrimal.boundaryFieldRef()[pI] *=
629 coeffFieldb/coeffFieldb.patchInternalField();
630 flux.boundaryFieldRef()[pI] = Zero;
631 }
632 else if (isA<fixedValueFvPatchScalarField>(bc))
633 {
634 snGradPrimal.boundaryFieldRef()[pI] = Zero;
635 flux.boundaryFieldRef()[pI] = Zero;
636 }
637 }
638
639 return coeffField*(fvc::div(flux) - adjointField*fvc::div(snGradPrimal));
640}
641
642
644(
645 const volVectorField& U,
646 const volVectorField& Ua,
647 const volScalarField& nut,
648 const word& schemeName
649) const
650{
652 (interpolationScheme<vector>(schemeName));
654 surfaceScalarField flux(scheme().interpolate(Ua) & snGradU);
655
656 // Terms form the Laplacian part of the momentum stresses
657 forAll(mesh_.boundary(), pI)
658 {
659 const fvPatchScalarField& bc = nut.boundaryField()[pI];
660 if (isA<zeroGradientFvPatchScalarField>(bc))
661 {
662 flux.boundaryFieldRef()[pI] = Zero;
663 }
664 else if (isA<fixedValueFvPatchScalarField>(bc))
665 {
666 snGradU.boundaryFieldRef()[pI] = Zero;
667 flux.boundaryFieldRef()[pI] = Zero;
668 }
669 }
670
671 // Terms form the tranpose gradient in the momentum stresses
672 const surfaceVectorField& Sf = mesh_.Sf();
673 surfaceTensorField fluxTranspose
674 (
676 );
677 forAll(mesh_.boundary(), pI)
678 {
679 const fvPatchVectorField& Ub = U.boundaryField()[pI];
680 if (!isA<coupledFvPatchVectorField>(Ub))
681 {
682 const vectorField Uai(Ua.boundaryField()[pI].patchInternalField());
683 const vectorField& Sfb = Sf.boundaryField()[pI];
684 fluxTranspose.boundaryFieldRef()[pI] = Uai*Sfb;
685 }
686 }
687 volScalarField M(dev2(gradU_) && fvc::div(fluxTranspose));
689 forAll(mesh_.boundary(), pI)
690 {
691 const fvPatchScalarField& bc = nut.boundaryField()[pI];
692 if (isA<zeroGradientFvPatchScalarField>(bc))
693 {
694 const vectorField Uai(Ua.boundaryField()[pI].patchInternalField());
695 const tensorField dev2GradU
696 (
697 dev2(gradU_.boundaryField()[pI].patchInternalField())
698 );
699 const vectorField& Sfb = Sf.boundaryField()[pI];
700 const labelList& faceCells = mesh_.boundary()[pI].faceCells();
701 forAll(faceCells, fI)
702 {
703 const label celli = faceCells[fI];
704 M[celli] -= ((Uai[fI] & dev2GradU[fI]) & Sfb[fI])/V[celli];
705 }
706 }
707 }
708 M.correctBoundaryConditions();
709
710 //surfaceScalarField fluxTranspose =
711 // fvc::interpolate(UaGradU, schemeName) & mesh_.Sf();
712 //forAll(mesh_.boundary(), pI)
713 //{
714 // const fvPatchScalarField& bc = nut.boundaryField()[pI];
715 // const vectorField& Sf = mesh_.boundary()[pI].Sf();
716 // if (isA<zeroGradientFvPatchScalarField>(bc))
717 // {
718 // fluxTranspose.boundaryFieldRef()[pI] =
719 // (
720 // UaGradU.boundaryField()[pI].patchInternalField()
721 // - (
722 // Ua.boundaryField()[pI].patchInternalField()
723 // & gradU_.boundaryField()[pI]
724 // )
725 // ) & Sf;
726 // }
727 // else if (isA<fixedValueFvPatchScalarField>(bc))
728 // {
729 // fluxTranspose.boundaryFieldRef()[pI] =
730 // UaGradU.boundaryField()[pI].patchInternalField() & Sf;
731 // }
732 //}
733 return
734 fvc::div(flux) - (Ua & fvc::div(snGradU)) + M;
735 //fvc::div(flux + fluxTranspose) - (Ua & fvc::div(snGradU));
736}
737
738
740(
741 const volScalarField& primalField,
742 const volScalarField& adjointField
743) const
744{
745 // Grab the interpolation scheme from the primal convection term
746 tmp<surfaceInterpolationScheme<scalar>> primalInterpolationScheme
747 (
748 convectionScheme(primalField.name())
749 );
750
751 surfaceVectorField interpolatedPrimal
752 (
753 primalInterpolationScheme().interpolate(primalField)*mesh_.Sf()
754 );
756 (
757 //reverseLinear<scalar>(mesh_).interpolate(adjointField)
758 linear<scalar>(mesh_).interpolate(adjointField)
759 *interpolatedPrimal
760 );
761
762 const volVectorField& U = primalVars_.U();
763 forAll(mesh_.boundary(), pI)
764 {
765 const fvPatchVectorField& bc = U.boundaryField()[pI];
766 if (isA<zeroGradientFvPatchVectorField>(bc))
767 {
768 flux.boundaryFieldRef()[pI] = Zero;
769 }
770 else if (isA<fixedValueFvPatchVectorField>(bc))
771 {
772 interpolatedPrimal.boundaryFieldRef()[pI] = Zero;
773 flux.boundaryFieldRef()[pI] = Zero;
774 }
775 }
776
777 return (-fvc::div(flux) + adjointField*fvc::div(interpolatedPrimal));
778}
779
780
782(
783 tmp<volSymmTensorField>& GbyNuMult
784) const
785{
787 (
789 );
790
791 const volVectorField& U = primalVars_.U();
792 forAll(mesh_.boundary(), pI)
793 {
794 const fvPatchVectorField& bc = U.boundaryField()[pI];
795 if (isA<zeroGradientFvPatchVectorField>(bc))
796 {
797 flux.boundaryFieldRef()[pI] = Zero;
798 }
799 else if (isA<fixedValueFvPatchVectorField>(bc))
800 {
801 flux.boundaryFieldRef()[pI] =
802 mesh_.boundary()[pI].Sf()
803 & GbyNuMult().boundaryField()[pI].patchInternalField();
804 }
805 }
806
807 return fvc::div(flux);
808}
809
810
812(
813 tmp<volScalarField>& divUMult
814) const
815{
817 (
819 );
820
821 const volVectorField& U = primalVars_.U();
822 forAll(mesh_.boundary(), pI)
823 {
824 const fvPatchVectorField& bc = U.boundaryField()[pI];
825 if (isA<zeroGradientFvPatchVectorField>(bc))
826 {
827 flux.boundaryFieldRef()[pI] = Zero;
828 }
829 else if (isA<fixedValueFvPatchVectorField>(bc))
830 {
831 flux.boundaryFieldRef()[pI] =
832 mesh_.boundary()[pI].Sf()
833 *divUMult().boundaryField()[pI].patchInternalField();
834 }
835 }
836
837 divUMult.clear();
838
839 return -fvc::div(flux);
840}
841
842
844(
845 const volScalarField& primalField,
846 const volScalarField& adjointField,
847 const volScalarField& coeffField
848) const
849{
850 // Note: we could grab the snGrad scheme from the diffusion term directly
851 surfaceScalarField snGradPrimal(fvc::snGrad(primalField)*mesh_.magSf());
853 (
854 reverseLinear<scalar>(mesh_).interpolate(adjointField)*snGradPrimal
855 );
856
857 const volVectorField& U = primalVars_.U();
858 forAll(mesh_.boundary(), pI)
859 {
860 const fvPatchVectorField& bc = U.boundaryField()[pI];
861 if (!isA<coupledFvPatchVectorField>(bc))
862 {
863 flux.boundaryFieldRef()[pI] = Zero;
864 snGradPrimal.boundaryFieldRef()[pI] = Zero;
865 }
866 }
867 return (fvc::div(flux) - adjointField*fvc::div(snGradPrimal))*coeffField;
868}
869
870
872(
874) const
875{
877 (
878 mult*a1_*k()*(1 - case_1_nut_)/(b1_*F2_*S_*S2_)*twoSymm(gradU_)
879 );
880 M.correctBoundaryConditions();
881 mult.clear();
882
883 surfaceVectorField returnFieldFlux
884 (
886 );
887
888 const volVectorField& U = primalVars_.U();
889 forAll(mesh_.boundary(), pI)
890 {
891 const fvPatchVectorField& bc = U.boundaryField()[pI];
892 if (isA<zeroGradientFvPatchVectorField>(bc))
893 {
894 returnFieldFlux.boundaryFieldRef()[pI] = Zero;
895 }
896 else if (isA<fixedValueFvPatchVectorField>(bc))
897 {
898 returnFieldFlux.boundaryFieldRef()[pI] =
899 mesh_.boundary()[pI].Sf()
900 & M.boundaryField()[pI].patchInternalField();
901 }
902 }
903
904 // Note: potentially missing contributions form patches with a calculated
905 // nut bc
906
907 return fvc::div(returnFieldFlux);
908}
909
910
912(
913 fvScalarMatrix& kaEqn,
914 const volScalarField& dR_dnut
915)
916{
917 // Add source term from the differentiation of nutWallFunction
918 scalarField& source = kaEqn.source();
920 const volScalarField& ka = this->ka();
921 const volScalarField& wa = this->wa();
922 const volScalarField& k = this->k();
923 const volScalarField& omega = this->omega();
924 forAll(nutRef().boundaryFieldRef(), patchi)
925 {
926 fvPatchScalarField& nutWall = nutRef().boundaryFieldRef()[patchi];
927 if (isA<nutkWallFunctionFvPatchScalarField>(nutWall))
928 {
929 const fvPatch& patch = mesh_.boundary()[patchi];
930 const scalarField& magSf = patch.magSf();
931
934
935 const scalarField& y = turbModel().y()[patchi];
936 const tmp<scalarField> tnuw = turbModel().nu(patchi);
937 const scalarField& nuw = tnuw();
938
940 refCast<nutWallFunctionFvPatchScalarField>(nutWall);
941 const wallFunctionCoefficients& wallCoeffs = nutWF.wallCoeffs();
942 const scalar Cmu = wallCoeffs.Cmu();
943 const scalar kappa = wallCoeffs.kappa();
944 const scalar E = wallCoeffs.E();
945 const scalar yPlusLam = wallCoeffs.yPlusLam();
946
947 const scalar Cmu25 = pow025(Cmu);
948
949 const labelList& faceCells = patch.faceCells();
950 const fvPatchScalarField& dR_dnutw =
951 dR_dnut.boundaryField()[patchi];
952 const fvPatchScalarField& omegaw = omega.boundaryField()[patchi];
953 bool addTermsFromOmegaWallFuction
954 (isA<omegaWallFunctionFvPatchScalarField>(omegaw));
955
956 const fvPatchVectorField& Uw =
957 primalVars_.U().boundaryField()[patchi];
958 const scalarField magGradUw(mag(Uw.snGrad()));
959 forAll(nuw, facei)
960 {
961 const label celli = faceCells[facei];
962
963 const scalar sqrtkCell(sqrt(k[celli]));
964 const scalar yPlus = Cmu25*y[facei]*sqrtkCell/nuw[facei];
965 const scalar logEyPlus = log(E*yPlus);
966 const scalar dnut_dyPlus =
967 nuw[facei]*kappa*(logEyPlus - 1)/sqr(logEyPlus);
968 const scalar dyPlus_dk =
969 Cmu25*y[facei]/(2*nuw[facei]*sqrtkCell);
970 const scalar dnut_dk = dnut_dyPlus*dyPlus_dk;
971
972 if (yPlusLam < yPlus)
973 {
974 // Terms from nutkWallFunction
975 source[celli] -= dR_dnutw[facei]*dnut_dk*magSf[facei];
976 }
977 if (addTermsFromOmegaWallFuction)
978 {
979 const scalar denom(Cmu25*kappa*y[facei]);
980 const scalar omegaLog(sqrtkCell/denom);
981 // Terms from omegaLog in omegaWallFunction
982 source[celli] +=
983 wa[celli]*omegaLog/omega[celli]
984 /(2*sqrtkCell*denom);
985
986 // Terms from G at the first cell off the wall
987 source[celli] +=
988 case_1_Pk_[celli]*ka[celli]*V[celli]
989 *(
990 (nutWall[facei] + nuw[facei])
991 *magGradUw[facei]
992 *Cmu25
993 /(2.0*sqrtkCell*kappa*y[facei])
994 );
995
996 if (yPlusLam < yPlus)
997 {
998 source[celli] +=
999 case_1_Pk_[celli]*ka[celli]*V[celli]
1000 *dnut_dk
1001 *magGradUw[facei]
1002 *Cmu25*sqrtkCell
1003 /(kappa*y[facei]);
1004 }
1005 }
1006 }
1007 }
1008 }
1009}
1010
1011
1013{
1015 {
1016 Info<< "Updating primal-based fields of the adjoint turbulence "
1017 << "model ..." << endl;
1018
1019 const volVectorField& U = primalVars_.U();
1020 gradU_ = fvc::grad(U);
1022 gradK_ = fvc::grad(k());
1023
1024 S2_ = 2*magSqr(symm(gradU_))
1026 S_ = sqrt(S2_);
1028
1029 // Instead of computing G directly here, delegate to RASModelVariables
1030 // to make sure G is computed consistently when the primal fields are
1031 // averaged. The local value computed here is just used to update
1032 // the switch fields
1033 /*
1034 volScalarField G(primalVars_.turbulence()->GName(), nutRef()*GbyNu0_);
1035 omega().correctBoundaryConditions();
1036 */
1038 (
1039 IOobject
1040 (
1041 type() + ":G",
1042 mesh_.time().timeName(),
1043 mesh_,
1046 ),
1047 mesh_,
1049 );
1050 G.ref() = primalVars_.RASModelVariables()->G()();
1051
1052 CDkOmega_ =
1055 (
1056 CDkOmega_,
1057 dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)
1058 );
1059 F1_ = F1();
1060 F2_ = F2();
1061 case_1_Pk_ = neg(G - c1_*betaStar_*k()*omega());
1062 case_2_Pk_ = pos0(G - c1_*betaStar_*k()*omega());
1063 case_3_Pk_ = neg(a1_*omega() - b1_*F2_*S_);
1064
1065
1066 alphaK_ = alphaK(F1_);
1068 beta_ = beta(F1_);
1069 gamma_ = gamma(F1_);
1070
1071 // Switch fields for F1
1072 {
1073 volScalarField arg1_C3
1074 (
1075 (scalar(1)/betaStar_)*sqrt(k())/(omega()*y_)
1076 - scalar(500)*nu()/(sqr(y_)*omega())
1077 );
1078 volScalarField arg1_C2
1079 (
1080 max
1081 (
1082 (scalar(1)/betaStar_)*sqrt(k())/(omega()*y_),
1083 scalar(500)*nu()/(sqr(y_)*omega())
1084 )
1085 - (4*alphaOmega2_)*k()/(CDkOmegaPlus_*sqr(y_))
1086 );
1087 volScalarField arg1_C1
1088 (
1089 min
1090 (
1091 max
1092 (
1093 (scalar(1)/betaStar_)*sqrt(k())/(omega()*y_),
1094 scalar(500)*nu()/(sqr(y_)*omega())
1095 ),
1097 ) - scalar(10)
1098 );
1099 volScalarField CDkOmegaPlus_C1
1100 (
1101 CDkOmega_
1102 - dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)
1103 );
1104
1105 case_1_F1_ = pos(arg1_C3)*neg(arg1_C2)*neg(arg1_C1);
1106 case_2_F1_ = neg0(arg1_C3)*neg(arg1_C2)*neg(arg1_C1);
1107 case_3_F1_ = pos0(arg1_C2)*neg(arg1_C1)*pos(CDkOmegaPlus_C1);
1108 case_4_F1_ = pos0(arg1_C2)*neg(arg1_C1);
1109 }
1110
1111 // Switch fields for nut
1112 {
1113 volScalarField nut_C1(a1_*omega() - b1_*F2_*S_);
1114 volScalarField arg2_C2
1115 (
1116 (scalar(2)/betaStar_)*sqrt(k())/(omega()*y_)
1117 - scalar(500)*nu()/(sqr(y_)*omega())
1118 );
1119 volScalarField arg2_C1
1120 (
1121 max
1122 (
1123 (scalar(2)/betaStar_)*sqrt(k())/(omega()*y_),
1124 scalar(500)*nu()/(sqr(y_)*omega())
1125 ) - scalar(100)
1126 );
1127
1128 case_1_nut_ = pos(nut_C1);
1129 case_2_nut_ = neg0(nut_C1)*pos(arg2_C2)*neg(arg2_C1);
1130 case_3_nut_ = neg0(nut_C1)*neg0(arg2_C2)*neg(arg2_C1);
1131 }
1132
1133 {
1134 volScalarField GPrime_C1
1135 (
1136 GbyNu0_
1138 );
1139 case_1_GPrime_ = neg(GPrime_C1);
1140 case_2_GPrime_ = pos0(GPrime_C1);
1141 }
1142
1144 dnut_dk_ = dnut_dk();
1146 DkEff_ = DkEff(F1_);
1147
1148 changedPrimalSolution_ = false;
1149 }
1150}
1151
1152
1154(
1155 const word& varName
1156) const
1157{
1159 const surfaceScalarField& phiInst = primalVars_.phiInst();
1160 word divEntry("div(" + phiInst.name() + ',' + varName +')');
1161 ITstream& divScheme = mesh_.divScheme(divEntry);
1162 // Skip the first entry which might be 'bounded' or 'Gauss'.
1163 // If it is 'bounded', skip the second entry as well
1164 word discarded(divScheme);
1165 if (discarded == "bounded")
1166 {
1167 discarded = word(divScheme);
1168 }
1170}
1171
1172
1173// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1174
1176(
1177 incompressibleVars& primalVars,
1179 objectiveManager& objManager,
1180 const word& adjointTurbulenceModelName,
1181 const word& modelName
1182)
1183:
1185 (
1186 modelName,
1187 primalVars,
1188 adjointVars,
1189 objManager,
1190 adjointTurbulenceModelName
1191 ),
1192
1193 kappa_
1194 (
1195 dimensioned<scalar>::getOrAddToDict
1196 (
1197 "kappa",
1198 coeffDict_,
1199 0.41
1200 )
1201 ),
1202 alphaK1_
1203 (
1204 dimensioned<scalar>::getOrAddToDict
1205 (
1206 "alphaK1",
1207 this->coeffDict_,
1208 0.85
1209 )
1210 ),
1211 alphaK2_
1212 (
1213 dimensioned<scalar>::getOrAddToDict
1214 (
1215 "alphaK2",
1216 this->coeffDict_,
1217 1.0
1218 )
1219 ),
1220 alphaOmega1_
1221 (
1222 dimensioned<scalar>::getOrAddToDict
1223 (
1224 "alphaOmega1",
1225 this->coeffDict_,
1226 0.5
1227 )
1228 ),
1229 alphaOmega2_
1230 (
1231 dimensioned<scalar>::getOrAddToDict
1232 (
1233 "alphaOmega2",
1234 this->coeffDict_,
1235 0.856
1236 )
1237 ),
1238 gamma1_
1239 (
1240 dimensioned<scalar>::getOrAddToDict
1241 (
1242 "gamma1",
1243 this->coeffDict_,
1244 5.0/9.0
1245 )
1246 ),
1247 gamma2_
1248 (
1249 dimensioned<scalar>::getOrAddToDict
1250 (
1251 "gamma2",
1252 this->coeffDict_,
1253 0.44
1254 )
1255 ),
1256 beta1_
1257 (
1258 dimensioned<scalar>::getOrAddToDict
1259 (
1260 "beta1",
1261 this->coeffDict_,
1262 0.075
1263 )
1264 ),
1265 beta2_
1266 (
1267 dimensioned<scalar>::getOrAddToDict
1268 (
1269 "beta2",
1270 this->coeffDict_,
1271 0.0828
1272 )
1273 ),
1274 betaStar_
1275 (
1276 dimensioned<scalar>::getOrAddToDict
1277 (
1278 "betaStar",
1279 this->coeffDict_,
1280 0.09
1281 )
1282 ),
1283 a1_
1284 (
1285 dimensioned<scalar>::lookupOrAddToDict
1286 (
1287 "a1",
1288 this->coeffDict_,
1289 0.31
1290 )
1291 ),
1292 b1_
1293 (
1294 dimensioned<scalar>::lookupOrAddToDict
1295 (
1296 "b1",
1297 this->coeffDict_,
1298 1.0
1299 )
1300 ),
1301 c1_
1302 (
1303 dimensioned<scalar>::lookupOrAddToDict
1304 (
1305 "c1",
1306 this->coeffDict_,
1307 10.0
1308 )
1309 ),
1310 F3_
1311 (
1312 Switch::lookupOrAddToDict
1313 (
1314 "F3",
1315 this->coeffDict_,
1316 false
1317 )
1318 ),
1319
1320 y_(primalVars_.RASModelVariables()().d()),
1321
1322 //Primal Gradient Fields
1323 gradU_
1324 (
1325 IOobject
1326 (
1327 "rasModel::gradU",
1328 runTime_.timeName(),
1329 mesh_,
1330 IOobject::NO_READ,
1331 IOobject::NO_WRITE
1332 ),
1333 mesh_,
1335 ),
1336 gradOmega_
1337 (
1338 IOobject
1339 (
1340 "rasModel::gradOmega",
1341 runTime_.timeName(),
1342 mesh_,
1343 IOobject::NO_READ,
1344 IOobject::NO_WRITE
1345 ),
1346 mesh_,
1348 ),
1349 gradK_
1350 (
1351 IOobject
1352 (
1353 "rasModel::gradK",
1354 runTime_.timeName(),
1355 mesh_,
1356 IOobject::NO_READ,
1357 IOobject::NO_WRITE
1358 ),
1359 mesh_,
1361 ),
1362
1363 S2_
1364 (
1365 IOobject
1366 (
1367 "S2",
1368 runTime_.timeName(),
1369 mesh_,
1370 IOobject::NO_READ,
1371 IOobject::NO_WRITE
1372 ),
1373 mesh_,
1375 ),
1376 S_
1377 (
1378 IOobject
1379 (
1380 "kOmegaSST_S",
1381 runTime_.timeName(),
1382 mesh_,
1383 IOobject::NO_READ,
1384 IOobject::NO_WRITE
1385 ),
1386 mesh_,
1388 ),
1389 GbyNu0_
1390 (
1391 IOobject
1392 (
1393 "adjointRASModel::GbyNu0",
1394 runTime_.timeName(),
1395 mesh_,
1396 IOobject::NO_READ,
1397 IOobject::NO_WRITE
1398 ),
1399 mesh_,
1401 ),
1402 CDkOmega_
1403 (
1404 IOobject
1405 (
1406 "CDkOmega_",
1407 runTime_.timeName(),
1408 mesh_,
1409 IOobject::NO_READ,
1410 IOobject::NO_WRITE
1411 ),
1412 mesh_,
1414 ),
1415 CDkOmegaPlus_
1416 (
1417 IOobject
1418 (
1419 "CDkOmegaPlus",
1420 runTime_.timeName(),
1421 mesh_,
1422 IOobject::NO_READ,
1423 IOobject::NO_WRITE
1424 ),
1425 mesh_,
1427 ),
1428 F1_
1429 (
1430 IOobject
1431 (
1432 "F1",
1433 runTime_.timeName(),
1434 mesh_,
1435 IOobject::NO_READ,
1436 IOobject::NO_WRITE
1437 ),
1438 mesh_,
1440 ),
1441 F2_
1442 (
1443 IOobject
1444 (
1445 "F2",
1446 runTime_.timeName(),
1447 mesh_,
1448 IOobject::NO_READ,
1449 IOobject::NO_WRITE
1450 ),
1451 mesh_,
1453 ),
1454 // Model Field coefficients
1455 alphaK_
1456 (
1457 IOobject
1458 (
1459 "alphaK",
1460 runTime_.timeName(),
1461 mesh_,
1462 IOobject::NO_READ,
1463 IOobject::NO_WRITE
1464 ),
1465 mesh_,
1467 ),
1468 alphaOmega_
1469 (
1470 IOobject
1471 (
1472 "alphaOmega",
1473 runTime_.timeName(),
1474 mesh_,
1475 IOobject::NO_READ,
1476 IOobject::NO_WRITE
1477 ),
1478 mesh_,
1480 ),
1481 beta_
1482 (
1483 IOobject
1484 (
1485 "beta",
1486 runTime_.timeName(),
1487 mesh_,
1488 IOobject::NO_READ,
1489 IOobject::NO_WRITE
1490 ),
1491 mesh_,
1493 ),
1494 gamma_
1495 (
1496 IOobject
1497 (
1498 "gamma",
1499 runTime_.timeName(),
1500 mesh_,
1501 IOobject::NO_READ,
1502 IOobject::NO_WRITE
1503 ),
1504 mesh_,
1506 ),
1507
1508 case_1_F1_
1509 (
1510 IOobject
1511 (
1512 "case_1_F1",
1513 runTime_.timeName(),
1514 mesh_,
1515 IOobject::NO_READ,
1516 IOobject::NO_WRITE
1517 ),
1518 mesh_,
1520 ),
1521 case_2_F1_
1522 (
1523 IOobject
1524 (
1525 "case_2_F1",
1526 runTime_.timeName(),
1527 mesh_,
1528 IOobject::NO_READ,
1529 IOobject::NO_WRITE
1530 ),
1531 mesh_,
1533 ),
1534 case_3_F1_
1535 (
1536 IOobject
1537 (
1538 "case_3_F1",
1539 runTime_.timeName(),
1540 mesh_,
1541 IOobject::NO_READ,
1542 IOobject::NO_WRITE
1543 ),
1544 mesh_,
1546 ),
1547 case_4_F1_
1548 (
1549 IOobject
1550 (
1551 "case_4_F1",
1552 runTime_.timeName(),
1553 mesh_,
1554 IOobject::NO_READ,
1555 IOobject::NO_WRITE
1556 ),
1557 mesh_,
1559 ),
1560 case_1_Pk_
1561 (
1562 IOobject
1563 (
1564 "case_1_Pk",
1565 runTime_.timeName(),
1566 mesh_,
1567 IOobject::NO_READ,
1568 IOobject::NO_WRITE
1569 ),
1570 mesh_,
1572 ),
1573 case_2_Pk_
1574 (
1575 IOobject
1576 (
1577 "case_2_Pk",
1578 runTime_.timeName(),
1579 mesh_,
1580 IOobject::NO_READ,
1581 IOobject::NO_WRITE
1582 ),
1583 mesh_,
1585 ),
1586 case_3_Pk_
1587 (
1588 IOobject
1589 (
1590 "case_3_Pk",
1591 runTime_.timeName(),
1592 mesh_,
1593 IOobject::NO_READ,
1594 IOobject::NO_WRITE
1595 ),
1596 mesh_,
1598 ),
1599
1600 case_1_nut_
1601 (
1602 IOobject
1603 (
1604 "case_1_nut",
1605 runTime_.timeName(),
1606 mesh_,
1607 IOobject::NO_READ,
1608 IOobject::NO_WRITE
1609 ),
1610 mesh_,
1612 ),
1613 case_2_nut_
1614 (
1615 IOobject
1616 (
1617 "case_2_nut",
1618 runTime_.timeName(),
1619 mesh_,
1620 IOobject::NO_READ,
1621 IOobject::NO_WRITE
1622 ),
1623 mesh_,
1625 ),
1626 case_3_nut_
1627 (
1628 IOobject
1629 (
1630 "case_3_nut",
1631 runTime_.timeName(),
1632 mesh_,
1633 IOobject::NO_READ,
1634 IOobject::NO_WRITE
1635 ),
1636 mesh_,
1638 ),
1639 case_1_GPrime_
1640 (
1641 IOobject
1642 (
1643 "case_1_GPrime",
1644 runTime_.timeName(),
1645 mesh_,
1646 IOobject::NO_READ,
1647 IOobject::NO_WRITE
1648 ),
1649 mesh_,
1651 ),
1652 case_2_GPrime_
1653 (
1654 IOobject
1655 (
1656 "case_2_GPrime",
1657 runTime_.timeName(),
1658 mesh_,
1659 IOobject::NO_READ,
1660 IOobject::NO_WRITE
1661 ),
1662 mesh_,
1664 ),
1665
1666 // Zero 1rst cell field
1667 firstCellIDs_(0),
1668 zeroFirstCell_(zeroFirstCell()),
1669
1670 // Turbulence model multipliers
1671 dnut_domega_
1672 (
1673 IOobject
1674 (
1675 "dnut_domega",
1676 runTime_.timeName(),
1677 mesh_,
1678 IOobject::NO_READ,
1679 IOobject::NO_WRITE
1680 ),
1681 mesh_,
1683 ),
1684 dnut_dk_
1685 (
1686 IOobject
1687 (
1688 "dnut_dk",
1689 runTime_.timeName(),
1690 mesh_,
1691 IOobject::NO_READ,
1692 IOobject::NO_WRITE
1693 ),
1694 mesh_,
1696 ),
1697 DOmegaEff_
1698 (
1699 IOobject
1700 (
1701 "DomegaEff",
1702 runTime_.timeName(),
1703 mesh_,
1704 IOobject::NO_READ,
1705 IOobject::NO_WRITE
1706 ),
1707 mesh_,
1708 dimensionedScalar(nutRef().dimensions(), Zero)
1709 ),
1710 DkEff_
1711 (
1712 IOobject
1713 (
1714 "DkEff",
1715 runTime_.timeName(),
1716 mesh_,
1717 IOobject::NO_READ,
1718 IOobject::NO_WRITE
1719 ),
1720 mesh_,
1721 dimensionedScalar(nutRef().dimensions(), Zero)
1722 )
1723{
1727 // Read in adjoint fields
1729 (
1731 mesh_,
1732 "ka",
1733 adjointVars.solverName(),
1734 adjointVars.useSolverNameForFields()
1735 );
1737 (
1739 mesh_,
1740 "wa",
1741 adjointVars.solverName(),
1742 adjointVars.useSolverNameForFields()
1743 );
1744
1745 setMeanFields();
1746
1747 // No sensitivity contributions from the adjoint to the eikonal equation
1748 // for the moment
1749 includeDistance_ = false;
1750
1751 // Update the primal related fields here so that functions computing
1752 // sensitivities have the updated fields in case of continuation
1754}
1755
1756
1757// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1758
1760{
1761 const volVectorField& Ua = adjointVars_.UaInst();
1762 return devReff(Ua);
1763}
1764
1765
1767(
1768 const volVectorField& U
1769) const
1770{
1771 return
1773 (
1774 IOobject
1775 (
1776 "devRhoReff",
1778 mesh_,
1781 ),
1783 );
1784}
1785
1786
1788{
1789 tmp<volScalarField> tnuEff = nuEff();
1790 const volScalarField& nuEff = tnuEff();
1791
1792 return
1793 (
1794 - fvm::laplacian(nuEff, Ua)
1795 - fvc::div(nuEff*dev(fvc::grad(Ua)().T()))
1796 );
1797
1798 /* WIP
1799 const volVectorField& U = primalVars_.U();
1800 const surfaceVectorField& Sf = mesh_.Sf();
1801 tmp<surfaceTensorField> tflux =
1802 reverseLinear<vector>(mesh_).interpolate(Ua)*Sf;
1803 surfaceTensorField& flux = tflux.ref();
1804 forAll(mesh_.boundary(), pI)
1805 {
1806 const fvPatchVectorField& Ub = U.boundaryField()[pI];
1807 if (!isA<coupledFvPatchVectorField>(Ub))
1808 {
1809 const vectorField Uai = Ua.boundaryField()[pI].patchInternalField();
1810 const vectorField& Sfb = Sf.boundaryField()[pI];
1811 flux.boundaryFieldRef()[pI] = Uai*Sfb;
1812 }
1813 }
1814 volTensorField M(nuEff*dev2(fvc::div(flux)));
1815 const DimensionedField<scalar, volMesh>& V = mesh_.V();
1816
1817 forAll(mesh_.boundary(), pI)
1818 {
1819 const fvPatchVectorField& Ub = U.boundaryField()[pI];
1820 if (!isA<coupledFvPatchVectorField>(Ub))
1821 {
1822 const fvPatchScalarField& nuEffb = nuEff.boundaryField()[pI];
1823 const vectorField nf = mesh_.boundary()[pI].nf();
1824 const vectorField Uai = Ua.boundaryField()[pI].patchInternalField();
1825 const labelList& faceCells = mesh_.boundary()[pI].faceCells();
1826 const vectorField& Sfb = Sf.boundaryField()[pI];
1827
1828 forAll(faceCells, fI)
1829 {
1830 const label celli = faceCells[fI];
1831 const tensor t(dev2(Uai[fI]*Sfb[fI]));
1832 M[celli] -= nuEffb[fI]*(t - nf[fI]*(nf[fI] & t))/V[celli];
1833 }
1834 }
1835 }
1836 M.correctBoundaryConditions();
1837
1838 surfaceVectorField returnFlux =
1839 - (Sf & reverseLinear<tensor>(mesh_).interpolate(M));
1840 forAll(mesh_.boundary(), pI)
1841 {
1842 const fvPatchVectorField& Ub = U.boundaryField()[pI];
1843 if (isA<zeroGradientFvPatchVectorField>(Ub))
1844 {
1845 returnFlux.boundaryFieldRef()[pI] = Zero;
1846 }
1847 else if (isA<fixedValueFvPatchVectorField>(Ub))
1848 {
1849 const scalarField& deltaCoeffs = mesh_.boundary()[pI].deltaCoeffs();
1850 const fvPatchScalarField& nuEffb = nuEff.boundaryField()[pI];
1851 const vectorField Uai = Ua.boundaryField()[pI].patchInternalField();
1852 const vectorField nf = mesh_.boundary()[pI].nf();
1853 const vectorField& Sfb = Sf.boundaryField()[pI];
1854
1855 returnFlux.boundaryFieldRef()[pI] =
1856 - (Sfb & M.boundaryField()[pI].patchInternalField())
1857 + nuEffb*deltaCoeffs*(nf & dev2(Uai*Sfb));
1858 }
1859 }
1860
1861 return
1862 (
1863 - fvm::laplacian(nuEff, Ua)
1864 + fvc::div(returnFlux)
1865 );
1866 */
1867}
1868
1869
1871{
1872 return (ka()*gradK_ + wa()*gradOmega_);
1873}
1874
1875
1877{
1878 tmp<volVectorField> tmeanFlowSource
1879 (
1881 (
1882 IOobject
1883 (
1884 "adjointMeanFlowSource" + type(),
1885 mesh_.time().timeName(),
1886 mesh_,
1889 ),
1890 mesh_,
1892 )
1893 );
1894 volVectorField& meanFlowSource = tmeanFlowSource.ref();
1895
1896 // Contributions from the convection terms of the turbulence model
1897 meanFlowSource +=
1900
1901 // Contributions from GbyNu, including gradU
1902 tmp<volSymmTensorField> twoSymmGradU(twoSymm(gradU_));
1903 tmp<volSymmTensorField> GbyNuMult
1904 (
1905 // First part of GPrime and G from Pk
1906 2.*dev(twoSymmGradU())*zeroFirstCell_
1908 // Second part of GPrime
1909 + twoSymmGradU()*wa()*zeroFirstCell_*gamma_*case_2_GPrime_
1911 );
1912 twoSymmGradU.clear();
1913 meanFlowSource += GMeanFlowSource(GbyNuMult);
1914 GbyNuMult.clear();
1915
1916 // Contributions from divU
1917 tmp<volScalarField> divUMult
1918 (
1919 (2.0/3.0)*(zeroFirstCell_*wa()*omega()*gamma_ + ka()*k())
1920 );
1921 meanFlowSource += divUMeanFlowSource(divUMult);
1922
1923 // Contributions from S2, existing in nut
1924 const volVectorField& U = primalVars_.U();
1925 const volVectorField& Ua = adjointVars_.UaInst();
1926 tmp<volScalarField> nutMeanFlowSourceMult
1927 (
1928 // nut in the diffusion coefficients
1931 + dNutdbMult(U, Ua, nutRef(), "coeffsDiff")
1932 // nut in G
1934 );
1935 meanFlowSource += nutMeanFlowSource(nutMeanFlowSourceMult);
1936
1937 // G at the first cell includes mag(U.snGrad())
1938 // Add term here
1939 forAll(omega().boundaryFieldRef(), patchi)
1940 {
1941 fvPatchScalarField& omegaWall = omega().boundaryFieldRef()[patchi];
1942 if (isA<omegaWallFunctionFvPatchScalarField>(omegaWall))
1943 {
1944 const fvPatch& patch = mesh_.boundary()[patchi];
1945
1948
1949 const scalarField& y = turbModel().y()[patchi];
1950 const tmp<scalarField> tnuw = turbModel().nu(patchi);
1951 const scalarField& nuw = tnuw();
1952
1954 refCast<nutWallFunctionFvPatchScalarField>
1955 (nutRef().boundaryFieldRef()[patchi]);
1956 const wallFunctionCoefficients& wallCoeffs = nutw.wallCoeffs();
1957 const scalar Cmu = wallCoeffs.Cmu();
1958 const scalar kappa = wallCoeffs.kappa();
1959 const scalar Cmu25 = pow025(Cmu);
1960
1961 const labelList& faceCells = patch.faceCells();
1962
1963 const fvPatchVectorField& Uw =
1964 primalVars_.U().boundaryField()[patchi];
1965 vectorField snGradUw(Uw.snGrad());
1966 const scalarField& deltaCoeffs = patch.deltaCoeffs();
1967 forAll(faceCells, facei)
1968 {
1969 const label celli = faceCells[facei];
1970 // Volume will be added when meanFlowSource is added to UaEqn
1971 meanFlowSource[celli] +=
1972 ka()[celli]*case_1_Pk_[celli]
1973 *(nutw[facei] + nuw[facei])
1974 *snGradUw[facei].normalise()
1975 *Cmu25*sqrt(k()[celli])
1976 *deltaCoeffs[facei]
1977 /(kappa*y[facei]);
1978 }
1979 }
1980 }
1981 return tmeanFlowSource;
1982}
1983
1984
1986{
1987 return dnut_dk_;
1988}
1989
1990
1992{
1993 return dnut_domega_;
1994}
1995
1996
1998{
1999 return
2000 (
2001 alphaK_.boundaryField()[patchI]*nutRef().boundaryField()[patchI]
2002 + nu()().boundaryField()[patchI]
2003 );
2004}
2005
2006
2008{
2009 return
2010 (
2011 alphaOmega_.boundaryField()[patchI]*nutRef().boundaryField()[patchI]
2012 + nu()().boundaryField()[patchI]
2013 );
2014}
2015
2016
2018{
2020
2021 if (!adjointTurbulence_)
2022 {
2023 return;
2024 }
2025
2027
2028 // Primal and adjoint fields
2029 const volVectorField& U = primalVars_.U();
2033
2035 (
2036 fvm::div(-phi, wa())
2042 + fvm::Sp(zeroFirstCell_*scalar(2.)*beta_*omega(), wa())
2043 + fvm::SuSp
2044 (
2045 zeroFirstCell_()*gamma_*((2.0/3.0)*divU - dGPrime_domega().ref()()),
2046 wa()
2047 )
2048 - (case_2_Pk_*c1_ - scalar(1))*betaStar_*k()*ka()
2049 );
2050
2051 // Boundary manipulate changes the diagonal component, so relax has to
2052 // come after that
2053 waEqn.ref().boundaryManipulate(wa().boundaryFieldRef());
2054 waEqn.ref().relax();
2055
2056 // Sources from the objective should be added after the boundary
2057 // manipulation
2059 waEqn.ref().solve();
2060
2061 // Adjoint Turbulent kinetic energy equation
2063 (
2064 fvm::div(-phi, ka())
2065 + fvm::SuSp(fvc::div(phi), ka())
2067 + fvm::Sp(betaStar_*omega(), ka())
2068 - case_2_Pk_()*c1_*betaStar_*omega()()*ka()()
2069 + fvm::SuSp(scalar(2.0/3.0)*divU, ka())
2072 + dR_dnut()*dnut_dk_()
2074 );
2075
2076 kaEqn.ref().relax();
2077 kaEqn.ref().boundaryManipulate(ka().boundaryFieldRef());
2079 // Add sources from the objective functions
2081
2082 kaEqn.ref().solve();
2083
2085 {
2086 dimensionedScalar maxwa = max(mag(wa()));
2087 dimensionedScalar maxka = max(mag(ka()));
2088 Info<< "Max mag of adjoint dissipation = " << maxwa.value() << endl;
2089 Info<< "Max mag of adjoint kinetic energy = " << maxka.value() << endl;
2090 }
2091}
2092
2093
2095{
2096 return adjMomentumBCSourcePtr_();
2097}
2098
2099
2101{
2104
2105 forAll(mesh_.boundary(), patchi)
2106 {
2107 vectorField nf(mesh_.boundary()[patchi].nf());
2108 wallShapeSens[patchi] = nf & FITerm.boundaryField()[patchi];
2109 }
2110 return wallShapeSens;
2111}
2112
2113
2115{
2117}
2118
2119
2121{
2123 (
2124 IOobject
2125 (
2126 "adjointEikonalSource" + type(),
2128 mesh_,
2131 ),
2132 mesh_,
2134 );
2135}
2136
2137
2139{
2140 const volVectorField& U = primalVars_.U();
2141 const volScalarField& kInst =
2142 primalVars_.RASModelVariables()->TMVar1Inst();
2143 const volScalarField& omegaInst =
2144 primalVars_.RASModelVariables()->TMVar2Inst();
2145
2147 (
2148 min
2149 (
2150 max
2151 (
2152 (scalar(1)/betaStar_)*sqrt(k())/(omega()*y_),
2153 scalar(500)*nu()/(sqr(y_)*omega())
2154 ),
2156 ),
2157 scalar(10)
2158 );
2159
2160 // Interpolation schemes used by the primal convection terms
2161 auto kScheme(convectionScheme(kInst.name()));
2162 auto omegaScheme(convectionScheme(omegaInst.name()));
2163 const surfaceVectorField& Sf = mesh_.Sf();
2164 tmp<volTensorField> tFISens
2165 (
2167 (
2168 IOobject
2169 (
2170 type() + "FISensTerm",
2171 mesh_.time().timeName(),
2172 mesh_,
2175 ),
2176 mesh_,
2178 zeroGradientFvPatchTensorField::typeName
2179 )
2180 );
2181 volTensorField& FISens = tFISens.ref();
2182 FISens =
2183 // k convection
2184 - ka()*fvc::div
2185 (
2186 kScheme().interpolate(k())
2188 )
2189 // k diffusion
2190 + ka()*T(fvc::grad(DkEff_*gradK_))
2191 - DkEff_*(fvc::grad(ka())*gradK_)
2192 // omega convection
2194 (
2195 omegaScheme().interpolate(omega())
2197 )
2198 // omega diffusion
2201 // terms including GbyNu0
2202 + (
2204 + case_1_Pk_*ka()*nutRef()
2206 // S2 (includes contribution from nut in UEqn as well)
2207 + (
2208 dR_dnut()*a1_*k()/(b1_*S_*S_*S_*F2_)
2211 )*T(gradU_ & twoSymm(gradU_))*(1. - case_1_nut_)
2212 // CDkOmega in omegaEqn
2213 + 2.*wa()*(1. - F1_)*alphaOmega2_/omega()*zeroFirstCell_
2215 // F1
2216 - dR_dF1()
2217 *(dF1_dGradK(arg1)*gradK_ + dF1_dGradOmega(arg1)*gradOmega_);
2218
2220
2221 return tFISens;
2222}
2223
2224
2226(
2227 const word& designVarsName
2228) const
2229{
2230 // Missing proper terms - return zero for now
2232}
2233
2234
2236{
2239}
2240
2241
2243{
2245 {
2256 a1_.readIfPresent(this->coeffDict());
2257 b1_.readIfPresent(this->coeffDict());
2258 c1_.readIfPresent(this->coeffDict());
2259 F3_.readIfPresent("F3", this->coeffDict());
2260
2261 return true;
2262 }
2263 else
2264 {
2265 return false;
2266 }
2267}
2268
2269
2270// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
2271
2272} // End namespace adjointRASModels
2273} // End namespace incompressible
2274} // End namespace Foam
2275
2276// ************************************************************************* //
#define F2(B, C, D)
Definition: SHA1.C:151
#define M(I)
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
surfaceScalarField & phi
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
Generic GeometricBoundaryField class.
Boundary & boundaryFieldRef(const bool updateAccessTime=true)
Return a reference to the boundary field.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
void correctBoundaryConditions()
Correct boundary field.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
const word & name() const noexcept
Return the object name.
Definition: IOobjectI.H:65
An input stream of tokens.
Definition: ITstream.H:56
void setSize(const label n)
Alias for resize()
Definition: List.H:218
void resize(const label len)
Adjust allocated size of list.
Definition: ListI.H:139
A simple wrapper around bool so that it can be read as a word: true/false, on/off,...
Definition: Switch.H:78
bool readIfPresent(const word &key, const dictionary &dict)
Update the value of the Switch if it is found in the dictionary.
Definition: Switch.C:335
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
static word timeName(const scalar t, const int precision=precision_)
Definition: Time.C:780
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
const Type & value() const
Return const reference to value.
bool readIfPresent(const dictionary &dict)
Abstract base class for finite area calculus convection schemes.
Smooth ATC in cells next to a set of patches supplied by type.
Definition: faceCells.H:59
UPtrList< const labelUList > faceCells() const
Return a list of faceCells for each patch.
A special matrix type and solver, designed for finite volume solutions of scalar equations....
Definition: fvMatrix.H:121
Field< Type > & source() noexcept
Definition: fvMatrix.H:458
const DimensionedField< scalar, volMesh > & V() const
Return cell volumes.
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:712
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:290
const surfaceVectorField & Sf() const
Return cell face area vectors.
const surfaceScalarField & magSf() const
Return cell face area magnitudes.
virtual tmp< Field< Type > > patchInternalField() const
Return internal field next to patch as patch field.
Definition: fvPatchField.C:237
virtual tmp< Field< Type > > snGrad() const
Return patch-normal gradient.
Definition: fvPatchField.C:229
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
Definition: fvPatch.H:71
Manages the adjoint mean flow fields and their mean values.
const volVectorField & UaInst() const
Return const reference to velocity.
const solverControl & getSolverControl() const
Return const reference to solverControl.
Abstract base class for incompressible turbulence models.
autoPtr< boundaryVectorField > wallShapeSensitivitiesPtr_
Wall sensitivity term for shape optimisation.
virtual void correct()
Solve the adjoint turbulence equations.
Switch adjointTurbulence_
Turbulence on/off flag.
bool changedPrimalSolution_
Has the primal solution changed?
autoPtr< boundaryVectorField > adjMomentumBCSourcePtr_
const nearWallDist & y() const
Return the near wall distances.
autoPtr< volScalarField > adjointTMVariable1Ptr_
Adjoint turbulence model variable 1.
autoPtr< volScalarField > adjointTMVariable2Ptr_
Adjoint turbulence model variable 2.
objectiveManager & objectiveManager_
Reference to the objectiveManager.
wordList adjointTMVariablesBaseNames_
Base names of the adjoint fields.
autoPtr< boundaryVectorField > wallFloCoSensitivitiesPtr_
Wall sensitivity term for flow control optimisation.
const dictionary & coeffDict() const
Const access to the coefficients dictionary.
virtual bool read()
Read adjointRASProperties dictionary.
Continuous adjoint to the kOmegaSST turbulence model for incompressible flows.
virtual tmp< fvVectorMatrix > divDevReff(volVectorField &U) const
Return the transpose part of the adjoint momentum stresses.
virtual tmp< scalarField > diffusionCoeffVar1(label patchI) const
Diffusion coeff at the boundary for k.
virtual tmp< volScalarField > nutJacobianTMVar2() const
Jacobian of nut wrt to omega.
tmp< volScalarField > dGPrime_domega() const
GbyNu Jacobian wrt omega.
tmp< volVectorField > GMeanFlowSource(tmp< volSymmTensorField > &GbyNuMult) const
Contributions from the G.
virtual tmp< volTensorField > FISensitivityTerm()
Sensitivity derivative contributions when using the FI approach.
tmp< volVectorField > dF1_dGradK(const volScalarField &arg1) const
F1 Jacobian wrt grad(k)
tmp< volScalarField > dR_dF1() const
Derivative of the primal equations wrt F1.
tmp< volScalarField > dF1_dk(const volScalarField &arg1) const
F1 Jacobian wrt k (no contributions from grad(k))
tmp< volScalarField > dnut_dk() const
Nut Jacobian wrt k.
tmp< volScalarField > dGPrime_dk() const
GbyNu Jacobian wrt k.
tmp< volVectorField > convectionMeanFlowSource(const volScalarField &primalField, const volScalarField &adjointField) const
Contributions from the turbulence model convection terms.
virtual void correct()
Solve the adjoint turbulence equations.
volScalarField DkEff_
Diffusivity of the k equation.
volScalarField S2_
Primal cached fields involved in the solution of the.
tmp< volScalarField > dNutdbMult(const volScalarField &primalField, const volScalarField &adjointField, const volScalarField &coeffField, const volScalarField &bcField, const word &schemeName) const
Term multiplying dnut/db, coming from the turbulence model.
tmp< fvScalarMatrix > waEqnSourceFromCDkOmega() const
Source to waEqn from the differentiation of CDkOmega.
tmp< volScalarField > kaEqnSourceFromF1() const
Source to kaEqn from the differentiation of F1.
tmp< volScalarField > dR_dnut()
Derivative of the primal equations wrt nut.
void updatePrimalRelatedFields()
Update of the primal cached fields.
tmp< volScalarField > dF2_dk() const
F2 Jacobian wrt k.
virtual const boundaryVectorField & adjointMomentumBCSource() const
tmp< surfaceInterpolationScheme< scalar > > convectionScheme(const word &varName) const
virtual tmp< volSymmTensorField > devReff() const
Return the effective stress tensor including the laminar stress.
virtual void nullify()
Nullify all adjoint turbulence model fields and their old times.
virtual tmp< volScalarField > nutJacobianTMVar1() const
Jacobian of nut wrt to k.
tmp< volVectorField > nutMeanFlowSource(tmp< volScalarField > &mult) const
Contributions from nut(U)
volTensorField gradU_
Cached primal gradient fields.
volScalarField DOmegaEff_
Diffusivity of the omega equation.
virtual tmp< volScalarField > GbyNu(const volScalarField &GbyNu0, const volScalarField &F2, const volScalarField &S2) const
tmp< volScalarField > alphaK(const volScalarField &F1) const
tmp< volScalarField > diffusionNutMeanFlowMult(const volScalarField &primalField, const volScalarField &adjointField, const volScalarField &coeffField) const
virtual const boundaryVectorField & wallShapeSensitivities()
Sensitivity derivative contributions when using the (E)SI approach.
tmp< volScalarField > DomegaEff(const volScalarField &F1) const
Return the effective diffusivity for omega.
tmp< volScalarField > alphaOmega(const volScalarField &F1) const
void addWallFunctionTerms(fvScalarMatrix &kaEqn, const volScalarField &dR_dnut)
tmp< volScalarField > waEqnSourceFromF1() const
Source to waEqn from the differentiation of F1.
volScalarField case_1_Pk_
Switch fields for the production in the k Eqn.
tmp< volScalarField::Internal > beta(const volScalarField::Internal &F1) const
virtual const boundaryVectorField & wallFloCoSensitivities()
Sensitivity terms for flow control, emerging from the.
virtual tmp< volVectorField > nonConservativeMomentumSource() const
Non-conservative part of the terms added to the mean flow equations.
tmp< volScalarField > kaEqnSourceFromCDkOmega() const
Source to kaEqn from the differentiation of CDkOmega.
tmp< volScalarField > coeffsDifferentiation(const volScalarField &primalField, const volScalarField &adjointField, const word &schemeName) const
Differentiation of the turbulence model diffusion coefficients.
virtual tmp< volScalarField > distanceSensitivities()
Contributions to the adjoint eikonal equation (zero for now)
virtual tmp< scalarField > diffusionCoeffVar2(label patchI) const
Diffusion coeff at the boundary for omega.
tmp< volVectorField > divUMeanFlowSource(tmp< volScalarField > &divUMult) const
Contributions from the divU.
tmp< volVectorField > dF1_dGradOmega(const volScalarField &arg1) const
F1 Jacobian wrt grad(omega)
virtual tmp< scalarField > topologySensitivities(const word &designVarsName) const
tmp< volScalarField > dnut_domega() const
Nut Jacobian wrt omega.
tmp< volScalarField > DkEff(const volScalarField &F1) const
Return the effective diffusivity for k.
tmp< volScalarField > dF1_domega(const volScalarField &arg1) const
F1 Jacobian wrt omega (no contributions from grad(omega))
virtual bool read()
Read adjointRASProperties dictionary.
tmp< volScalarField > dF2_domega() const
F2 Jacobian wrt omega.
virtual tmp< volScalarField > nuEff() const
Return the effective viscosity.
tmp< volScalarField > nu() const
Return the laminar viscosity.
virtual const volScalarField & nut() const
Return the turbulence viscosity.
Base class for solution control classes.
const volVectorField & U() const
Return const reference to velocity.
const surfaceScalarField & phi() const
Return const reference to volume flux.
const surfaceScalarField & phiInst() const
Return const reference to volume flux.
const autoPtr< incompressible::turbulenceModel > & turbulence() const
Return const reference to the turbulence model.
const autoPtr< incompressible::RASModelVariables > & RASModelVariables() const
Return const reference to the turbulence model variables.
Central-differencing interpolation scheme class.
Definition: linear.H:58
The class nutWallFunction is an abstract base class that hosts calculation methods and common functi...
const wallFunctionCoefficients & wallCoeffs() const noexcept
Return wallFunctionCoefficients.
class for managing incompressible objective functions.
virtual void addTMEqn2Source(fvScalarMatrix &adjTMEqn2)=0
Add contribution to second adjoint turbulence model PDE.
virtual void addTMEqn1Source(fvScalarMatrix &adjTMEqn1)=0
Add contribution to first adjoint turbulence model PDE.
This boundary condition provides a wall function for the specific dissipation rate (i....
A class representing the concept of 1 (one) that can be used to avoid manipulating objects known to b...
Definition: one.H:62
label nCells() const noexcept
Number of mesh cells.
Inversed weight central-differencing interpolation scheme class.
Definition: reverseLinear.H:60
ITstream & divScheme(const word &name) const
Get div scheme for given name, or default.
bool printMaxMags() const
Print max mags of solver fields.
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &, const tmp< surfaceScalarField > &, const tmp< surfaceScalarField > &)
Return the face-interpolate of the given cell field.
A class for managing temporary objects.
Definition: tmp.H:65
void clear() const noexcept
Definition: tmpI.H:287
T & ref() const
Definition: tmpI.H:227
static void setField(autoPtr< GeometricField< Type, fvPatchField, volMesh > > &fieldPtr, const fvMesh &mesh, const word &baseName, const word &solverName, const bool useSolverNameForFields)
Read vol fields.
const word & solverName() const
Return solver name.
Definition: variablesSet.C:84
bool useSolverNameForFields() const
Append solver name to fields?
Definition: variablesSet.C:90
static void nullifyField(GeometricField< Type, PatchField, GeoMesh > &fieldPtr)
Nullify field and old times, if present.
Class to host the wall-function coefficients being used in the wall function boundary conditions.
scalar kappa() const noexcept
Return the object: kappa.
scalar E() const noexcept
Return the object: E.
scalar yPlusLam() const noexcept
Return the object: yPlusLam.
scalar Cmu() const noexcept
Return the object: Cmu.
A class for handling words, derived from Foam::string.
Definition: word.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
U
Definition: pEqn.H:72
const volScalarField & T
rDeltaT ref()
const scalar gamma
Definition: EEqn.H:9
scalar yPlus
scalar nut
zeroField divU
Definition: alphaSuSp.H:3
word timeName
Definition: getTimeIndex.H:3
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh > > grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:54
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:47
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:49
tmp< surfaceScalarField > absolute(const tmp< surfaceScalarField > &tphi, const volVectorField &U)
Return the given relative flux in absolute form.
Definition: fvcMeshPhi.C:190
tmp< fvMatrix< Type > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmLaplacian.C:48
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:48
zeroField SuSp(const Foam::zero, const GeometricField< Type, fvPatchField, volMesh > &)
A no-op source.
zeroField Sp(const Foam::zero, const GeometricField< Type, fvPatchField, volMesh > &)
A no-op source.
Namespace for OpenFOAM.
dimensionedSymmTensor dev2(const dimensionedSymmTensor &dt)
dimensionedSymmTensor dev(const dimensionedSymmTensor &dt)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
dimensionedScalar pos(const dimensionedScalar &ds)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
dimensionedSymmTensor symm(const dimensionedSymmTensor &dt)
const dimensionSet dimless
Dimensionless.
dimensionedScalar pos0(const dimensionedScalar &ds)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:52
dimensionedScalar pow3(const dimensionedScalar &ds)
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:53
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:59
dimensionedSymmTensor twoSymm(const dimensionedSymmTensor &dt)
const dimensionSet dimVelocity
dimensionedScalar tanh(const dimensionedScalar &ds)
messageStream Info
Information stream (stdout output on master, null elsewhere)
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
Definition: MSwindows.C:598
dimensionedScalar log(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
dimensionedScalar pow4(const dimensionedScalar &ds)
dimensionedScalar neg(const dimensionedScalar &ds)
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
bool interpolate(const vector &p1, const vector &p2, const vector &o, vector &n, scalar l)
Definition: curveTools.C:75
dimensionedScalar neg0(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
dimensionedScalar pow025(const dimensionedScalar &ds)
scalar kb
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333