src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "MULES.H"
  27 #include "upwind.H"
  28 #include "uncorrectedSnGrad.H"
  29 #include "gaussConvectionScheme.H"
  30 #include "gaussLaplacianScheme.H"
  31 #include "uncorrectedSnGrad.H"
  32 #include "surfaceInterpolate.H"
  33 #include "fvcSurfaceIntegrate.H"
  34 #include "slicedSurfaceFields.H"
  35 #include "syncTools.H"
  36
  37 #include "fvm.H"
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40
  41 template<class RhoType, class SpType, class SuType>
  42 void Foam::MULES::explicitSolve
  43 (
  44     const RhoType& rho,
  45     volScalarField& psi,
  46     const surfaceScalarField& phi,
  47     surfaceScalarField& phiPsi,
  48     const SpType& Sp,
  49     const SuType& Su,
  50     const scalar psiMax,
  51     const scalar psiMin
  52 )
  53 {
  54     Info<< "MULES: Solving for " << psi.name() << endl;
  55
  56     const fvMesh& mesh = psi.mesh();
  57     psi.correctBoundaryConditions();
  58
  59     surfaceScalarField phiBD(upwind<scalar>(psi.mesh(), phi).flux(psi));
  60
  61     surfaceScalarField& phiCorr = phiPsi;
  62     phiCorr -= phiBD;
  63
  64     scalarField allLambda(mesh.nFaces(), 1.0);
  65
  66     slicedSurfaceScalarField lambda
  67     (
  68         IOobject
  69         (
  70             "lambda",
  71             mesh.time().timeName(),
  72             mesh,
  73             IOobject::NO_READ,
  74             IOobject::NO_WRITE,
  75             false
  76         ),
  77         mesh,
  78         dimless,
  79         allLambda,
  80         false   // Use slices for the couples
  81     );
  82
  83     limiter
  84     (
  85         allLambda,
  86         rho,
  87         psi,
  88         phiBD,
  89         phiCorr,
  90         Sp,
  91         Su,
  92         psiMax,
  93         psiMin,
  94         3
  95     );
  96
  97     phiPsi = phiBD + lambda*phiCorr;
  98
  99     scalarField& psiIf = psi;
 100     const scalarField& psi0 = psi.oldTime();
 101     const scalar deltaT = mesh.time().deltaTValue();
 102
 103     psiIf = 0.0;
 104     fvc::surfaceIntegrate(psiIf, phiPsi);
 105
 106     if (mesh.moving())
 107     {
 108         psiIf =
 109         (
 110             mesh.Vsc0()().field()*rho.oldTime().field()
 111            *psi0/(deltaT*mesh.Vsc()().field())
 112           + Su.field()
 113           - psiIf
 114         )/(rho.field()/deltaT - Sp.field());
 115     }
 116     else
 117     {
 118         psiIf =
 119         (
 120             rho.oldTime().field()*psi0/deltaT
 121           + Su.field()
 122           - psiIf
 123         )/(rho.field()/deltaT - Sp.field());
 124     }
 125
 126     psi.correctBoundaryConditions();
 127 }
 128
 129
 130 template<class RhoType, class SpType, class SuType>
 131 void Foam::MULES::implicitSolve
 132 (
 133     const RhoType& rho,
 134     volScalarField& psi,
 135     const surfaceScalarField& phi,
 136     surfaceScalarField& phiPsi,
 137     const SpType& Sp,
 138     const SuType& Su,
 139     const scalar psiMax,
 140     const scalar psiMin
 141 )
 142 {
 143     const fvMesh& mesh = psi.mesh();
 144
 145     const dictionary& MULEScontrols = mesh.solverDict(psi.name());
 146
 147     label maxIter
 148     (
 149         readLabel(MULEScontrols.lookup("maxIter"))
 150     );
 151
 152     label nLimiterIter
 153     (
 154         readLabel(MULEScontrols.lookup("nLimiterIter"))
 155     );
 156
 157     scalar maxUnboundedness
 158     (
 159         readScalar(MULEScontrols.lookup("maxUnboundedness"))
 160     );
 161
 162     scalar CoCoeff
 163     (
 164         readScalar(MULEScontrols.lookup("CoCoeff"))
 165     );
 166
 167     scalarField allCoLambda(mesh.nFaces());
 168
 169     {
 170         tmp<surfaceScalarField> Cof =
 171             mesh.time().deltaT()*mesh.surfaceInterpolation::deltaCoeffs()
 172            *mag(phi)/mesh.magSf();
 173
 174         slicedSurfaceScalarField CoLambda
 175         (
 176             IOobject
 177             (
 178                 "CoLambda",
 179                 mesh.time().timeName(),
 180                 mesh,
 181                 IOobject::NO_READ,
 182                 IOobject::NO_WRITE,
 183                 false
 184             ),
 185             mesh,
 186             dimless,
 187             allCoLambda,
 188             false   // Use slices for the couples
 189         );
 190
 191         CoLambda == 1.0/max(CoCoeff*Cof, scalar(1));
 192     }
 193
 194     scalarField allLambda(allCoLambda);
 195     //scalarField allLambda(mesh.nFaces(), 1.0);
 196
 197     slicedSurfaceScalarField lambda
 198     (
 199         IOobject
 200         (
 201             "lambda",
 202             mesh.time().timeName(),
 203             mesh,
 204             IOobject::NO_READ,
 205             IOobject::NO_WRITE,
 206             false
 207         ),
 208         mesh,
 209         dimless,
 210         allLambda,
 211         false   // Use slices for the couples
 212     );
 213
 214     linear<scalar> CDs(mesh);
 215     upwind<scalar> UDs(mesh, phi);
 216     //fv::uncorrectedSnGrad<scalar> snGrads(mesh);
 217
 218     fvScalarMatrix psiConvectionDiffusion
 219     (
 220         fvm::ddt(rho, psi)
 221       + fv::gaussConvectionScheme<scalar>(mesh, phi, UDs).fvmDiv(phi, psi)
 222         //- fv::gaussLaplacianScheme<scalar, scalar>(mesh, CDs, snGrads)
 223         //.fvmLaplacian(Dpsif, psi)
 224       - fvm::Sp(Sp, psi)
 225       - Su
 226     );
 227
 228     surfaceScalarField phiBD(psiConvectionDiffusion.flux());
 229
 230     surfaceScalarField& phiCorr = phiPsi;
 231     phiCorr -= phiBD;
 232
 233     for (label i=0; i<maxIter; i++)
 234     {
 235         if (i != 0 && i < 4)
 236         {
 237             allLambda = allCoLambda;
 238         }
 239
 240         limiter
 241         (
 242             allLambda,
 243             rho,
 244             psi,
 245             phiBD,
 246             phiCorr,
 247             Sp,
 248             Su,
 249             psiMax,
 250             psiMin,
 251             nLimiterIter
 252         );
 253
 254         solve
 255         (
 256             psiConvectionDiffusion + fvc::div(lambda*phiCorr),
 257             MULEScontrols
 258         );
 259
 260         scalar maxPsiM1 = gMax(psi.internalField()) - 1.0;
 261         scalar minPsi = gMin(psi.internalField());
 262
 263         scalar unboundedness = max(max(maxPsiM1, 0.0), -min(minPsi, 0.0));
 264
 265         if (unboundedness < maxUnboundedness)
 266         {
 267             break;
 268         }
 269         else
 270         {
 271             Info<< "MULES: max(" << psi.name() << " - 1) = " << maxPsiM1
 272                 << " min(" << psi.name() << ") = " << minPsi << endl;
 273
 274             phiBD = psiConvectionDiffusion.flux();
 275
 276             /*
 277             word gammaScheme("div(phi,gamma)");
 278             word gammarScheme("div(phirb,gamma)");
 279
 280             const surfaceScalarField& phir =
 281                 mesh.lookupObject<surfaceScalarField>("phir");
 282
 283             phiCorr =
 284                 fvc::flux
 285                 (
 286                     phi,
 287                     psi,
 288                     gammaScheme
 289                 )
 290               + fvc::flux
 291                 (
 292                     -fvc::flux(-phir, scalar(1) - psi, gammarScheme),
 293                     psi,
 294                     gammarScheme
 295                 )
 296                 - phiBD;
 297             */
 298         }
 299     }
 300
 301     phiPsi = psiConvectionDiffusion.flux() + lambda*phiCorr;
 302 }
 303
 304
 305 template<class RhoType, class SpType, class SuType>
 306 void Foam::MULES::limiter
 307 (
 308     scalarField& allLambda,
 309     const RhoType& rho,
 310     const volScalarField& psi,
 311     const surfaceScalarField& phiBD,
 312     const surfaceScalarField& phiCorr,
 313     const SpType& Sp,
 314     const SuType& Su,
 315     const scalar psiMax,
 316     const scalar psiMin,
 317     const label nLimiterIter
 318 )
 319 {
 320     const scalarField& psiIf = psi;
 321     const volScalarField::GeometricBoundaryField& psiBf = psi.boundaryField();
 322
 323     const scalarField& psi0 = psi.oldTime();
 324
 325     const fvMesh& mesh = psi.mesh();
 326
 327     const labelUList& owner = mesh.owner();
 328     const labelUList& neighb = mesh.neighbour();
 329     tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
 330     const scalarField& V = tVsc();
 331     const scalar deltaT = mesh.time().deltaTValue();
 332
 333     const scalarField& phiBDIf = phiBD;
 334     const surfaceScalarField::GeometricBoundaryField& phiBDBf =
 335         phiBD.boundaryField();
 336
 337     const scalarField& phiCorrIf = phiCorr;
 338     const surfaceScalarField::GeometricBoundaryField& phiCorrBf =
 339         phiCorr.boundaryField();
 340
 341     slicedSurfaceScalarField lambda
 342     (
 343         IOobject
 344         (
 345             "lambda",
 346             mesh.time().timeName(),
 347             mesh,
 348             IOobject::NO_READ,
 349             IOobject::NO_WRITE,
 350             false
 351         ),
 352         mesh,
 353         dimless,
 354         allLambda,
 355         false   // Use slices for the couples
 356     );
 357
 358     scalarField& lambdaIf = lambda;
 359     surfaceScalarField::GeometricBoundaryField& lambdaBf =
 360         lambda.boundaryField();
 361
 362     scalarField psiMaxn(psiIf.size(), psiMin);
 363     scalarField psiMinn(psiIf.size(), psiMax);
 364
 365     scalarField sumPhiBD(psiIf.size(), 0.0);
 366
 367     scalarField sumPhip(psiIf.size(), VSMALL);
 368     scalarField mSumPhim(psiIf.size(), VSMALL);
 369
 370     forAll(phiCorrIf, facei)
 371     {
 372         label own = owner[facei];
 373         label nei = neighb[facei];
 374
 375         psiMaxn[own] = max(psiMaxn[own], psiIf[nei]);
 376         psiMinn[own] = min(psiMinn[own], psiIf[nei]);
 377
 378         psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]);
 379         psiMinn[nei] = min(psiMinn[nei], psiIf[own]);
 380
 381         sumPhiBD[own] += phiBDIf[facei];
 382         sumPhiBD[nei] -= phiBDIf[facei];
 383
 384         scalar phiCorrf = phiCorrIf[facei];
 385
 386         if (phiCorrf > 0.0)
 387         {
 388             sumPhip[own] += phiCorrf;
 389             mSumPhim[nei] += phiCorrf;
 390         }
 391         else
 392         {
 393             mSumPhim[own] -= phiCorrf;
 394             sumPhip[nei] -= phiCorrf;
 395         }
 396     }
 397
 398     forAll(phiCorrBf, patchi)
 399     {
 400         const fvPatchScalarField& psiPf = psiBf[patchi];
 401         const scalarField& phiBDPf = phiBDBf[patchi];
 402         const scalarField& phiCorrPf = phiCorrBf[patchi];
 403
 404         const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
 405
 406         if (psiPf.coupled())
 407         {
 408             const scalarField psiPNf(psiPf.patchNeighbourField());
 409
 410             forAll(phiCorrPf, pFacei)
 411             {
 412                 label pfCelli = pFaceCells[pFacei];
 413
 414                 psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]);
 415                 psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]);
 416             }
 417         }
 418         else
 419         {
 420             forAll(phiCorrPf, pFacei)
 421             {
 422                 label pfCelli = pFaceCells[pFacei];
 423
 424                 psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]);
 425                 psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]);
 426             }
 427         }
 428
 429         forAll(phiCorrPf, pFacei)
 430         {
 431             label pfCelli = pFaceCells[pFacei];
 432
 433             sumPhiBD[pfCelli] += phiBDPf[pFacei];
 434
 435             scalar phiCorrf = phiCorrPf[pFacei];
 436
 437             if (phiCorrf > 0.0)
 438             {
 439                 sumPhip[pfCelli] += phiCorrf;
 440             }
 441             else
 442             {
 443                 mSumPhim[pfCelli] -= phiCorrf;
 444             }
 445         }
 446     }
 447
 448     psiMaxn = min(psiMaxn, psiMax);
 449     psiMinn = max(psiMinn, psiMin);
 450
 451     //scalar smooth = 0.5;
 452     //psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax);
 453     //psiMinn = max((1.0 - smooth)*psiIf + smooth*psiMinn, psiMin);
 454
 455     if (mesh.moving())
 456     {
 457         tmp<volScalarField::DimensionedInternalField> V0 = mesh.Vsc0();
 458
 459         psiMaxn =
 460             V
 461            *(
 462                (rho.field()/deltaT - Sp.field())*psiMaxn
 463              - Su.field()
 464             )
 465           - (V0().field()/deltaT)*rho.oldTime().field()*psi0
 466           + sumPhiBD;
 467
 468         psiMinn =
 469             V
 470            *(
 471                Su.field()
 472              - (rho.field()/deltaT - Sp.field())*psiMinn
 473             )
 474           + (V0().field()/deltaT)*rho.oldTime().field()*psi0
 475           - sumPhiBD;
 476     }
 477     else
 478     {
 479         psiMaxn =
 480             V
 481            *(
 482                (rho.field()/deltaT - Sp.field())*psiMaxn
 483              - Su.field()
 484              - (rho.oldTime().field()/deltaT)*psi0
 485             )
 486           + sumPhiBD;
 487
 488         psiMinn =
 489             V
 490            *(
 491                Su.field()
 492              - (rho.field()/deltaT - Sp.field())*psiMinn
 493              + (rho.oldTime().field()/deltaT)*psi0
 494             )
 495           - sumPhiBD;
 496     }
 497
 498     scalarField sumlPhip(psiIf.size());
 499     scalarField mSumlPhim(psiIf.size());
 500
 501     for (int j=0; j<nLimiterIter; j++)
 502     {
 503         sumlPhip = 0.0;
 504         mSumlPhim = 0.0;
 505
 506         forAll(lambdaIf, facei)
 507         {
 508             label own = owner[facei];
 509             label nei = neighb[facei];
 510
 511             scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
 512
 513             if (lambdaPhiCorrf > 0.0)
 514             {
 515                 sumlPhip[own] += lambdaPhiCorrf;
 516                 mSumlPhim[nei] += lambdaPhiCorrf;
 517             }
 518             else
 519             {
 520                 mSumlPhim[own] -= lambdaPhiCorrf;
 521                 sumlPhip[nei] -= lambdaPhiCorrf;
 522             }
 523         }
 524
 525         forAll(lambdaBf, patchi)
 526         {
 527             scalarField& lambdaPf = lambdaBf[patchi];
 528             const scalarField& phiCorrfPf = phiCorrBf[patchi];
 529
 530             const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
 531
 532             forAll(lambdaPf, pFacei)
 533             {
 534                 label pfCelli = pFaceCells[pFacei];
 535
 536                 scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei];
 537
 538                 if (lambdaPhiCorrf > 0.0)
 539                 {
 540                     sumlPhip[pfCelli] += lambdaPhiCorrf;
 541                 }
 542                 else
 543                 {
 544                     mSumlPhim[pfCelli] -= lambdaPhiCorrf;
 545                 }
 546             }
 547         }
 548
 549         forAll(sumlPhip, celli)
 550         {
 551             sumlPhip[celli] =
 552                 max(min
 553                 (
 554                     (sumlPhip[celli] + psiMaxn[celli])/mSumPhim[celli],
 555                     1.0), 0.0
 556                 );
 557
 558             mSumlPhim[celli] =
 559                 max(min
 560                 (
 561                     (mSumlPhim[celli] + psiMinn[celli])/sumPhip[celli],
 562                     1.0), 0.0
 563                 );
 564         }
 565
 566         const scalarField& lambdam = sumlPhip;
 567         const scalarField& lambdap = mSumlPhim;
 568
 569         forAll(lambdaIf, facei)
 570         {
 571             if (phiCorrIf[facei] > 0.0)
 572             {
 573                 lambdaIf[facei] = min
 574                 (
 575                     lambdaIf[facei],
 576                     min(lambdap[owner[facei]], lambdam[neighb[facei]])
 577                 );
 578             }
 579             else
 580             {
 581                 lambdaIf[facei] = min
 582                 (
 583                     lambdaIf[facei],
 584                     min(lambdam[owner[facei]], lambdap[neighb[facei]])
 585                 );
 586             }
 587         }
 588
 589
 590         forAll(lambdaBf, patchi)
 591         {
 592             fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
 593             const scalarField& phiCorrfPf = phiCorrBf[patchi];
 594
 595             const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
 596
 597             forAll(lambdaPf, pFacei)
 598             {
 599                 label pfCelli = pFaceCells[pFacei];
 600
 601                 if (phiCorrfPf[pFacei] > 0.0)
 602                 {
 603                     lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdap[pfCelli]);
 604                 }
 605                 else
 606                 {
 607                     lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdam[pfCelli]);
 608                 }
 609             }
 610         }
 611
 612         syncTools::syncFaceList(mesh, allLambda, minEqOp<scalar>());
 613     }
 614 }
 615
 616
 617 // ************************************************************************* //