49 alphaPorosityPtr_(nullptr),
50 laminarTransportPtr_(nullptr),
51 turbulencePtr_(nullptr),
52 daTurbulenceModelPtr_(nullptr),
53 daFvSourcePtr_(nullptr),
54 fvSourcePtr_(nullptr),
65 const fvSolution& myFvSolution =
meshPtr_->thisDb().lookupObject<fvSolution>(
"fvSolution");
66 if (myFvSolution.found(
"relaxationFactors"))
68 if (myFvSolution.subDict(
"relaxationFactors").found(
"equations"))
70 if (myFvSolution.subDict(
"relaxationFactors").subDict(
"equations").found(
"U"))
72 relaxUEqn_ = myFvSolution.subDict(
"relaxationFactors").subDict(
"equations").getScalar(
"U");
76 solverDictU_ = myFvSolution.subDict(
"solvers").subDict(
"U");
77 solverDictP_ = myFvSolution.subDict(
"solvers").subDict(
"p");
88 Info <<
"Initializing fields for DASimpleFoam" << endl;
95 const word turbModelName(
98 "turbulenceProperties",
99 mesh.time().constant(),
105 .lookup(
"RASModel"));
112 if (
allOptions.subDict(
"fvSource").toc().size() != 0)
115 Info <<
"Initializing DASource" << endl;
116 word sourceName =
allOptions.subDict(
"fvSource").toc()[0];
117 word fvSourceType =
allOptions.subDict(
"fvSource").subDict(sourceName).getWord(
"type");
131 #include "createFvOptions.H"
136 Info <<
"\nStarting time loop\n"
144 Info <<
"Time = " <<
runTime.timeName() << nl << endl;
214 label printInterval =
daOptionPtr_->getOption<label>(
"printInterval");
216 word adjEqnSolMethod =
daOptionPtr_->getOption<word>(
"adjEqnSolMethod");
218 if (adjEqnSolMethod ==
"fixedPoint")
220 Info <<
"Solving the adjoint using consistent fixed-point iteration method..."
221 <<
" Execution Time: " <<
meshPtr_->time().elapsedCpuTime() <<
" s" << endl;
223 label fpMaxIters =
daOptionPtr_->getSubDictOption<label>(
"adjEqnOption",
"fpMaxIters");
224 scalar fpRelTol =
daOptionPtr_->getSubDictOption<scalar>(
"adjEqnOption",
"fpRelTol");
225 scalar fpMinResTolDiff =
daOptionPtr_->getSubDictOption<scalar>(
"adjEqnOption",
"fpMinResTolDiff");
227 const objectRegistry& db =
meshPtr_->thisDb();
228 volVectorField&
U =
const_cast<volVectorField&
>(db.lookupObject<volVectorField>(
"U"));
229 volScalarField&
p =
const_cast<volScalarField&
>(db.lookupObject<volScalarField>(
"p"));
230 surfaceScalarField&
phi =
const_cast<surfaceScalarField&
>(db.lookupObject<surfaceScalarField>(
"phi"));
231 volScalarField& nuTilda =
const_cast<volScalarField&
>(db.lookupObject<volScalarField>(
"nuTilda"));
245 volVectorField URes(
"URes", 0.0 *
U);
246 volScalarField pRes(
"pRes", 0.0 *
p);
247 surfaceScalarField phiRes(
"phiRes", 0.0 *
phi);
248 volScalarField nuTildaRes(
"nuTildaRes", 0.0 * nuTilda);
264 volVectorField UPsi(
"UPsi", 0.0 *
U);
265 volScalarField pPsi(
"pPsi", 0.0 *
p);
266 surfaceScalarField phiPsi(
"phiPsi", 0.0 *
phi);
267 volScalarField nuTildaPsi(
"nuTildaPsi", 0.0 * nuTilda);
303 volVectorField adjURes(
"adjURes", 0.0 *
U);
304 volScalarField adjPRes(
"adjPRes", 0.0 *
p);
305 surfaceScalarField adjPhiRes(
"adjPhiRes", 0.0 *
phi);
306 volScalarField adjNuTildaRes(
"adjNuTildaRes", 0.0 * nuTilda);
309 vector initNormAdjURes = vector::zero;
310 scalar initNormAdjPRes = 0.0;
311 scalar initNormAdjPhiRes = 0.0;
312 scalar initNormAdjNuTildaRes = 0.0;
340 volScalarField
p1(
"p1",
p);
341 volScalarField
p2(
"p2",
p);
342 this->
meshPtr_->setFluxRequired(p2.name());
343 surfaceScalarField
phi1(
"phi1",
phi);
344 volVectorField
U1(
"U1",
U);
347 volVectorField gradP(
"gradP", fvc::grad(
p));
348 surfaceScalarField pEqnFlux(
"pEqnFlux", 0.0 *
phi);
349 volScalarField divPhi(
"divPhi", fvc::div(
phi));
350 surfaceScalarField fluxU(
"fluxU", 0.0 *
phi);
365 while (cnt < fpMaxIters)
367 if (cnt % printInterval == 0)
369 Info <<
"Step = " << cnt <<
" Execution Time: " <<
meshPtr_->time().elapsedCpuTime() <<
" s" << endl;
378 this->
calcAdjointResidual(URes, pRes, phiRes, nuTildaRes, dFdW, UPsi, pPsi, phiPsi, nuTildaPsi, adjURes, adjPRes, adjPhiRes, adjNuTildaRes, cnt);
417 vector normAdjURes =
L2norm(adjURes);
418 scalar normAdjPRes =
L2norm(adjPRes);
419 scalar normAdjPhiRes =
L2norm(adjPhiRes);
420 scalar normAdjNuTildaRes =
L2norm(adjNuTildaRes);
424 initNormAdjURes = normAdjURes;
425 initNormAdjPRes = normAdjPRes;
426 initNormAdjPhiRes = normAdjPhiRes;
427 initNormAdjNuTildaRes = normAdjNuTildaRes;
431 normAdjURes[0] /= initNormAdjURes[0];
432 normAdjURes[1] /= initNormAdjURes[1];
433 normAdjURes[2] /= initNormAdjURes[2];
434 normAdjPRes /= initNormAdjPRes;
435 normAdjPhiRes /= initNormAdjPhiRes;
436 normAdjNuTildaRes /= initNormAdjNuTildaRes;
438 if (cnt % printInterval == 0)
440 Info <<
"Normalized L2 norm of adjURes: " << normAdjURes[0] <<
" " << normAdjURes[1] <<
" " << normAdjURes[2] << endl;
441 Info <<
"Normalized L2 norm of adjPRes: " << normAdjPRes << endl;
442 Info <<
"Normalized L2 norm of adjPhiRes: " << normAdjPhiRes << endl;
443 Info <<
"Normalized L2 norm of adjNuTildaRes: " << normAdjNuTildaRes << endl;
447 if (normAdjURes[0] < fpRelTol && normAdjURes[1] < fpRelTol && normAdjURes[2] < fpRelTol && normAdjPRes < fpRelTol && normAdjPhiRes < fpRelTol && normAdjNuTildaRes < fpRelTol)
449 Info <<
"Residual drop of " << fpRelTol <<
" has been achieved after " << cnt <<
" steps! Execution Time: " <<
meshPtr_->time().elapsedCpuTime() <<
" s" << endl;
450 Info <<
"Adjoint is converged and succussful.." << endl;
451 Info <<
"Normalized L2 norm of adjURes: " << normAdjURes[0] <<
" " << normAdjURes[1] <<
" " << normAdjURes[2] << endl;
452 Info <<
"Normalized L2 norm of adjPRes: " << normAdjPRes << endl;
453 Info <<
"Normalized L2 norm of adjPhiRes: " << normAdjPhiRes << endl;
454 Info <<
"Normalized L2 norm of adjNuTildaRes: " << normAdjNuTildaRes << endl;
468 adjURes[cellI] /=
UDiag[cellI];
479 volScalarField& step_2 = adjPRes;
487 this->globalADTape_.setActive();
492 this->globalADTape_.registerInput(
p1[cellI]);
496 gradPStart_ = this->globalADTape_.getPosition();
499 p1.correctBoundaryConditions();
500 U.correctBoundaryConditions();
503 gradP = fvc::grad(
p1);
506 gradP.primitiveFieldRef()[cellI] *=
U.mesh().V()[cellI];
510 gradPEnd_ = this->globalADTape_.getPosition();
513 this->globalADTape_.setPassive();
520 gradP[cellI][0].setGradient(adjURes[cellI][0].getValue());
521 gradP[cellI][1].setGradient(adjURes[cellI][1].getValue());
522 gradP[cellI][2].setGradient(adjURes[cellI][2].getValue());
526 this->globalADTape_.evaluate(gradPEnd_, gradPStart_);
529 step_2[cellI] +=
p1[cellI].getGradient();
533 this->globalADTape_.clearAdjoints();
537 scalar alphaP =
p.mesh().fieldRelaxationFactor(
p.name());
549 fvScalarMatrix pEqn1(
550 fvm::laplacian(
rAU,
p2));
554 this->globalADTape_.setActive();
559 this->globalADTape_.registerInput(
p2[cellI]);
563 pEqnFluxStart_ = this->globalADTape_.getPosition();
566 p2.correctBoundaryConditions();
567 U.correctBoundaryConditions();
570 pEqnFlux = pEqn1.flux();
573 pEqnFluxEnd_ = this->globalADTape_.getPosition();
576 this->globalADTape_.setPassive();
582 forAll(pEqnFlux.primitiveFieldRef(), faceI)
584 pEqnFlux.primitiveFieldRef()[faceI].setGradient(adjPhiRes.primitiveFieldRef()[faceI].getValue());
587 forAll(pEqnFlux.boundaryFieldRef(), patchI)
589 forAll(pEqnFlux.boundaryFieldRef()[patchI], faceI)
591 pEqnFlux.boundaryFieldRef()[patchI][faceI].setGradient(adjPhiRes.boundaryFieldRef()[patchI][faceI].getValue());
596 this->globalADTape_.evaluate(pEqnFluxEnd_, pEqnFluxStart_);
599 step_2[cellI] +=
p2[cellI].getGradient();
603 this->globalADTape_.clearAdjoints();
621 volScalarField& pSource = step_2;
626 volScalarField& step_3 =
pseudoP;
646 surfaceScalarField& phiTemp = adjPhiRes;
652 this->globalADTape_.setActive();
659 this->globalADTape_.registerInput(
phi1.primitiveFieldRef()[faceI]);
664 forAll(
phi1.boundaryFieldRef()[patchI], faceI)
666 this->globalADTape_.registerInput(
phi1.boundaryFieldRef()[patchI][faceI]);
671 divPhiStart_ = this->globalADTape_.getPosition();
679 divPhi = fvc::div(
phi1);
682 divPhi.primitiveFieldRef()[cellI] *=
p.mesh().V()[cellI];
686 divPhiEnd_ = this->globalADTape_.getPosition();
689 this->globalADTape_.setPassive();
696 divPhi[cellI].setGradient(step_3[cellI].getValue());
700 this->globalADTape_.evaluate(divPhiEnd_, divPhiStart_);
703 phiTemp.primitiveFieldRef()[faceI] +=
phi1.primitiveFieldRef()[faceI].getGradient();
707 forAll(
phi1.boundaryFieldRef()[patchI], faceI)
709 phiTemp.boundaryFieldRef()[patchI][faceI] +=
phi1.boundaryFieldRef()[patchI][faceI].getGradient();
714 this->globalADTape_.clearAdjoints();
731 UPsi[cellI] -= adjURes[cellI];
737 volVectorField& step_4 = adjURes;
742 this->globalADTape_.setActive();
747 this->globalADTape_.registerInput(
U1[cellI][0]);
748 this->globalADTape_.registerInput(
U1[cellI][1]);
749 this->globalADTape_.registerInput(
U1[cellI][2]);
753 fluxUStart_ = this->globalADTape_.getPosition();
756 U1.correctBoundaryConditions();
757 p.correctBoundaryConditions();
760 fluxU = fvc::flux(
U1);
763 fluxUEnd_ = this->globalADTape_.getPosition();
766 this->globalADTape_.setPassive();
772 forAll(fluxU.primitiveFieldRef(), faceI)
774 fluxU.primitiveFieldRef()[faceI].setGradient(phiTemp.primitiveFieldRef()[faceI].getValue());
777 forAll(fluxU.boundaryFieldRef(), patchI)
779 forAll(fluxU.boundaryFieldRef()[patchI], faceI)
781 fluxU.boundaryFieldRef()[patchI][faceI].setGradient(phiTemp.boundaryFieldRef()[patchI][faceI].getValue());
786 this->globalADTape_.evaluate(fluxUEnd_, fluxUStart_);
789 for (label cmpt = 0; cmpt < 3; cmpt++)
791 step_4[cellI][cmpt] -=
U1[cellI][cmpt].getGradient() /
UDiag[cellI];
796 this->globalADTape_.clearAdjoints();
805 step_4.primitiveFieldRef() = -
pseudoUEqn.lduMatrix::H(step_4);
819 volVectorField& USource = step_4;
824 volVectorField& step_5 =
pseudoU;
830 UPsi[cellI] += step_5[cellI];
845 volScalarField& nuTildaSource = adjNuTildaRes;
863 Info <<
"Prescribed adjoint convergence of " << fpRelTol <<
" has not been met after " << cnt <<
" steps! Execution Time: " <<
meshPtr_->time().elapsedCpuTime() <<
" s" << endl;
864 Info <<
"Now checking if the relaxed convergence is met..." << endl;
865 this->
calcAdjointResidual(URes, pRes, phiRes, nuTildaRes, dFdW, UPsi, pPsi, phiPsi, nuTildaPsi, adjURes, adjPRes, adjPhiRes, adjNuTildaRes, cnt);
867 vector normAdjURes =
L2norm(adjURes);
868 scalar normAdjPRes =
L2norm(adjPRes);
869 scalar normAdjPhiRes =
L2norm(adjPhiRes);
870 scalar normAdjNuTildaRes =
L2norm(adjNuTildaRes);
873 normAdjURes[0] /= initNormAdjURes[0];
874 normAdjURes[1] /= initNormAdjURes[1];
875 normAdjURes[2] /= initNormAdjURes[2];
876 normAdjPRes /= initNormAdjPRes;
877 normAdjPhiRes /= initNormAdjPhiRes;
878 normAdjNuTildaRes /= initNormAdjNuTildaRes;
880 scalar relaxedFpRelTol = fpRelTol * fpMinResTolDiff;
881 if (normAdjURes[0] < relaxedFpRelTol && normAdjURes[1] < relaxedFpRelTol && normAdjURes[2] < relaxedFpRelTol && normAdjPRes < relaxedFpRelTol && normAdjPhiRes < relaxedFpRelTol && normAdjNuTildaRes < relaxedFpRelTol)
883 Info <<
"Relaxed residual drop of " << relaxedFpRelTol <<
" has been achieved! Execution Time: " <<
meshPtr_->time().elapsedCpuTime() <<
" s" << endl;
884 Info <<
"Adjoint is still considered succussful.." << endl;
889 Info <<
"Relaxed residual drop of " << relaxedFpRelTol <<
" has not been met! Execution Time: " <<
meshPtr_->time().elapsedCpuTime() <<
" s" << endl;
890 Info <<
"Adjoint has failed to converge.." << endl;
892 Info <<
"Normalized L2 norm of adjURes: " << normAdjURes[0] <<
" " << normAdjURes[1] <<
" " << normAdjURes[2] << endl;
893 Info <<
"Normalized L2 norm of adjPRes: " << normAdjPRes << endl;
894 Info <<
"Normalized L2 norm of adjPhiRes: " << normAdjPhiRes << endl;
895 Info <<
"Normalized L2 norm of adjNuTildaRes: " << normAdjNuTildaRes << endl;
899 this->
vec2Fields(
"field2Vec", psi, UPsi, pPsi, phiPsi, nuTildaPsi);
903 FatalErrorIn(
"adjEqnSolMethod not valid") << exit(FatalError);
913 volVectorField& UField,
914 volScalarField& pField,
915 surfaceScalarField& phiField,
916 volScalarField& nuTildaField)
919 PetscScalar* cVecArray;
920 if (
mode ==
"vec2Field")
922 VecGetArray(cVec, &cVecArray);
927 for (label comp = 0; comp < 3; comp++)
929 label adjLocalIdx =
daIndexPtr_->getLocalAdjointStateIndex(
"U", cellI, comp);
930 UField[cellI][comp] = cVecArray[adjLocalIdx];
936 label adjLocalIdx =
daIndexPtr_->getLocalAdjointStateIndex(
"p", cellI);
937 pField[cellI] = cVecArray[adjLocalIdx];
942 label adjLocalIdx =
daIndexPtr_->getLocalAdjointStateIndex(
"phi", faceI);
944 if (faceI < daIndexPtr_->nLocalInternalFaces)
946 phiField[faceI] = cVecArray[adjLocalIdx];
950 label relIdx = faceI -
daIndexPtr_->nLocalInternalFaces;
953 phiField.boundaryFieldRef()[patchIdx][faceIdx] = cVecArray[adjLocalIdx];
959 label adjLocalIdx =
daIndexPtr_->getLocalAdjointStateIndex(
"nuTilda", cellI);
960 nuTildaField[cellI] = cVecArray[adjLocalIdx];
963 VecRestoreArray(cVec, &cVecArray);
965 else if (
mode ==
"field2Vec")
967 VecGetArray(cVec, &cVecArray);
972 for (label comp = 0; comp < 3; comp++)
974 label adjLocalIdx =
daIndexPtr_->getLocalAdjointStateIndex(
"U", cellI, comp);
975 cVecArray[adjLocalIdx] = UField[cellI][comp].value();
981 label adjLocalIdx =
daIndexPtr_->getLocalAdjointStateIndex(
"p", cellI);
982 cVecArray[adjLocalIdx] = pField[cellI].value();
987 label adjLocalIdx =
daIndexPtr_->getLocalAdjointStateIndex(
"phi", faceI);
989 if (faceI < daIndexPtr_->nLocalInternalFaces)
991 cVecArray[adjLocalIdx] = phiField[faceI].value();
995 label relIdx = faceI -
daIndexPtr_->nLocalInternalFaces;
998 cVecArray[adjLocalIdx] = phiField.boundaryFieldRef()[patchIdx][faceIdx].value();
1004 label adjLocalIdx =
daIndexPtr_->getLocalAdjointStateIndex(
"nuTilda", cellI);
1005 cVecArray[adjLocalIdx] = nuTildaField[cellI].value();
1008 VecRestoreArray(cVec, &cVecArray);
1012 FatalErrorIn(
"mode not valid") << exit(FatalError);
1018 const volVectorField& mySource,
1074 const volScalarField& mySource,
1141 volVectorField& URes,
1142 volScalarField& pRes,
1143 surfaceScalarField& phiRes)
1145 const objectRegistry& db =
meshPtr_->thisDb();
1146 const volVectorField&
U = db.lookupObject<volVectorField>(
"U");
1147 const volScalarField&
p = db.lookupObject<volScalarField>(
"p");
1148 const surfaceScalarField&
phi = db.lookupObject<surfaceScalarField>(
"phi");
1151 fvVectorMatrix
UEqn(
1153 - fvm::laplacian(nuEff,
U)
1154 - fvc::div(nuEff * dev2(
T(fvc::grad(
U)))));
1156 List<vector>& USource =
UEqn.source();
1162 volVectorField gradP(fvc::grad(
p));
1165 for (label i = 0; i <
U.size(); i++)
1167 URes[i] =
UDiag[i] *
U[i] - USource[i] +
U.mesh().V()[i] * gradP[i];
1169 URes.primitiveFieldRef() -=
UEqn.lduMatrix::H(
U);
1172 forAll(
U.boundaryField(), patchI)
1174 const fvPatch& pp =
U.boundaryField()[patchI].patch();
1179 label cellI = pp.faceCells()[faceI];
1182 for (label cmpt = 0; cmpt < 3; cmpt++)
1184 URes[cellI][cmpt] +=
UEqn.internalCoeffs()[patchI][faceI][cmpt] *
U[cellI][cmpt];
1187 URes[cellI] -=
UEqn.boundaryCoeffs()[patchI][faceI];
1192 URes.correctBoundaryConditions();
1196 volScalarField
rAU(1.0 /
UEqn.A());
1197 volVectorField
HbyA(
"HbyA",
U);
1199 surfaceScalarField
phiHbyA(
"phiHbyA", fvc::flux(
HbyA));
1201 fvScalarMatrix pEqn(
1204 List<scalar>& pSource = pEqn.source();
1205 List<scalar>& pDiag = pEqn.diag();
1208 for (label i = 0; i <
p.size(); i++)
1210 pRes[i] = pDiag[i] *
p[i] - pSource[i];
1212 pRes.primitiveFieldRef() -= pEqn.lduMatrix::H(
p);
1215 forAll(
p.boundaryField(), patchI)
1217 const fvPatch& pp =
p.boundaryField()[patchI].patch();
1222 label cellI = pp.faceCells()[faceI];
1226 pRes[cellI] += pEqn.internalCoeffs()[patchI][faceI] *
p[cellI];
1227 pRes[cellI] -= pEqn.boundaryCoeffs()[patchI][faceI];
1232 pRes.correctBoundaryConditions();
1240 volVectorField& URes,
1241 volScalarField& pRes,
1242 surfaceScalarField& phiRes,
1243 volScalarField& nuTildaRes,
1245 volVectorField& UPsi,
1246 volScalarField& pPsi,
1247 surfaceScalarField& phiPsi,
1248 volScalarField& nuTildaPsi,
1249 volVectorField& adjURes,
1250 volScalarField& adjPRes,
1251 surfaceScalarField& adjPhiRes,
1252 volScalarField& adjNuTildaRes,
1256 volVectorField&
U =
const_cast<volVectorField&
>(
meshPtr_->thisDb().lookupObject<volVectorField>(
"U"));
1257 volScalarField&
p =
const_cast<volScalarField&
>(
meshPtr_->thisDb().lookupObject<volScalarField>(
"p"));
1258 volScalarField& nuTilda =
const_cast<volScalarField&
>(
meshPtr_->thisDb().lookupObject<volScalarField>(
"nuTilda"));
1259 surfaceScalarField&
phi =
const_cast<surfaceScalarField&
>(
meshPtr_->thisDb().lookupObject<surfaceScalarField>(
"phi"));
1262 this->
vec2Fields(
"vec2Field", dFdW, adjURes, adjPRes, adjPhiRes, adjNuTildaRes);
1265 adjPhiRes = -adjPhiRes;
1266 adjNuTildaRes = -adjNuTildaRes;
1270 this->globalADTape_.reset();
1271 this->globalADTape_.setActive();
1277 this->globalADTape_.registerInput(
U[cellI][0]);
1278 this->globalADTape_.registerInput(
U[cellI][1]);
1279 this->globalADTape_.registerInput(
U[cellI][2]);
1284 this->globalADTape_.registerInput(
p[cellI]);
1291 this->globalADTape_.registerInput(
phi.primitiveFieldRef()[faceI]);
1296 forAll(
phi.boundaryFieldRef()[patchI], faceI)
1298 this->globalADTape_.registerInput(
phi.boundaryFieldRef()[patchI][faceI]);
1304 this->globalADTape_.registerInput(nuTilda[cellI]);
1308 adjResStart_ = this->globalADTape_.getPosition();
1311 U.correctBoundaryConditions();
1312 p.correctBoundaryConditions();
1313 nuTilda.correctBoundaryConditions();
1354 adjResEnd_ = this->globalADTape_.getPosition();
1357 this->globalADTape_.setPassive();
1364 URes[cellI][0].setGradient(UPsi[cellI][0].getValue());
1365 URes[cellI][1].setGradient(UPsi[cellI][1].getValue());
1366 URes[cellI][2].setGradient(UPsi[cellI][2].getValue());
1371 pRes[cellI].setGradient(pPsi[cellI].getValue());
1375 forAll(phiRes.primitiveFieldRef(), faceI)
1377 phiRes.primitiveFieldRef()[faceI].setGradient(phiPsi.primitiveFieldRef()[faceI].getValue());
1380 forAll(phiRes.boundaryFieldRef(), patchI)
1382 forAll(phiRes.boundaryFieldRef()[patchI], faceI)
1384 phiRes.boundaryFieldRef()[patchI][faceI].setGradient(phiPsi.boundaryFieldRef()[patchI][faceI].getValue());
1388 forAll(nuTildaRes, cellI)
1390 nuTildaRes[cellI].setGradient(nuTildaPsi[cellI].getValue());
1394 this->globalADTape_.evaluate(adjResEnd_, adjResStart_);
1397 adjURes[cellI][0] +=
U[cellI][0].getGradient();
1398 adjURes[cellI][1] +=
U[cellI][1].getGradient();
1399 adjURes[cellI][2] +=
U[cellI][2].getGradient();
1403 adjPRes[cellI] +=
p[cellI].getGradient();
1407 adjPhiRes.primitiveFieldRef()[faceI] +=
phi.primitiveFieldRef()[faceI].getGradient();
1411 forAll(
phi.boundaryFieldRef()[patchI], faceI)
1413 adjPhiRes.boundaryFieldRef()[patchI][faceI] +=
phi.boundaryFieldRef()[patchI][faceI].getGradient();
1418 adjNuTildaRes[cellI] += nuTilda[cellI].getGradient();
1422 this->globalADTape_.clearAdjoints();
1430 scalar L2normV = 0.0;
1434 L2normV += sqr(v[cellI] /
meshPtr_->V()[cellI]);
1436 L2normV = sqrt(L2normV);
1444 vector L2normU = vector::zero;
1448 L2normU[0] += sqr(
U[cellI][0] /
meshPtr_->V()[cellI]);
1449 L2normU[1] += sqr(
U[cellI][1] /
meshPtr_->V()[cellI]);
1450 L2normU[2] += sqr(
U[cellI][2] /
meshPtr_->V()[cellI]);
1452 L2normU[0] = sqrt(L2normU[0]);
1453 L2normU[1] = sqrt(L2normU[1]);
1454 L2normU[2] = sqrt(L2normU[2]);
1462 scalar L2normPhi = 0.0;
1464 forAll(Phi.primitiveField(), faceI)
1466 L2normPhi += sqr(Phi.primitiveField()[faceI]);
1468 forAll(Phi.boundaryField(), patchI)
1470 forAll(Phi.boundaryField()[patchI], faceI)
1472 L2normPhi += sqr(Phi.boundaryField()[patchI][faceI]);
1475 L2normPhi = sqrt(L2normPhi);
1483 List<scalar> temp = a;