applications/solvers/multiphase/interFoam/LTSInterFoam/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 "fvCFD.H"
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40
  41 template<class RhoType, class SpType, class SuType>
  42 void Foam::MULES::explicitLTSSolve
  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.field(),
  91         Su.field(),
  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
 102     const volScalarField& rDeltaT =
 103         mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
 104
 105     psiIf = 0.0;
 106     fvc::surfaceIntegrate(psiIf, phiPsi);
 107
 108     if (mesh.moving())
 109     {
 110         psiIf =
 111         (
 112             mesh.Vsc0()*rho.oldTime()*psi0*rDeltaT/mesh.Vsc()
 113           + Su.field()
 114           - psiIf
 115         )/(rho*rDeltaT - Sp.field());
 116     }
 117     else
 118     {
 119         psiIf =
 120         (
 121             rho.oldTime()*psi0*rDeltaT
 122           + Su.field()
 123           - psiIf
 124         )/(rho*rDeltaT - Sp.field());
 125     }
 126
 127     psi.correctBoundaryConditions();
 128 }
 129
 130
 131 template<class RhoType, class SpType, class SuType>
 132 void Foam::MULES::implicitSolve
 133 (
 134     const RhoType& rho,
 135     volScalarField& psi,
 136     const surfaceScalarField& phi,
 137     surfaceScalarField& phiPsi,
 138     const SpType& Sp,
 139     const SuType& Su,
 140     const scalar psiMax,
 141     const scalar psiMin
 142 )
 143 {
 144     const fvMesh& mesh = psi.mesh();
 145
 146     const dictionary& MULEScontrols = mesh.solverDict(psi.name());
 147
 148     label maxIter
 149     (
 150         readLabel(MULEScontrols.lookup("maxIter"))
 151     );
 152
 153     label nLimiterIter
 154     (
 155         readLabel(MULEScontrols.lookup("nLimiterIter"))
 156     );
 157
 158     scalar maxUnboundedness
 159     (
 160         readScalar(MULEScontrols.lookup("maxUnboundedness"))
 161     );
 162
 163     scalar CoCoeff
 164     (
 165         readScalar(MULEScontrols.lookup("CoCoeff"))
 166     );
 167
 168     scalarField allCoLambda(mesh.nFaces());
 169
 170     {
 171         surfaceScalarField Cof
 172         (
 173             mesh.time().deltaT()*mesh.surfaceInterpolation::deltaCoeffs()
 174            *mag(phi)/mesh.magSf()
 175         );
 176
 177         slicedSurfaceScalarField CoLambda
 178         (
 179             IOobject
 180             (
 181                 "CoLambda",
 182                 mesh.time().timeName(),
 183                 mesh,
 184                 IOobject::NO_READ,
 185                 IOobject::NO_WRITE,
 186                 false
 187             ),
 188             mesh,
 189             dimless,
 190             allCoLambda,
 191             false   // Use slices for the couples
 192         );
 193
 194         CoLambda == 1.0/max(CoCoeff*Cof, scalar(1));
 195     }
 196
 197     scalarField allLambda(allCoLambda);
 198     //scalarField allLambda(mesh.nFaces(), 1.0);
 199
 200     slicedSurfaceScalarField lambda
 201     (
 202         IOobject
 203         (
 204             "lambda",
 205             mesh.time().timeName(),
 206             mesh,
 207             IOobject::NO_READ,
 208             IOobject::NO_WRITE,
 209             false
 210         ),
 211         mesh,
 212         dimless,
 213         allLambda,
 214         false   // Use slices for the couples
 215     );
 216
 217     linear<scalar> CDs(mesh);
 218     upwind<scalar> UDs(mesh, phi);
 219     //fv::uncorrectedSnGrad<scalar> snGrads(mesh);
 220
 221     fvScalarMatrix psiConvectionDiffusion
 222     (
 223         fvm::ddt(rho, psi)
 224       + fv::gaussConvectionScheme<scalar>(mesh, phi, UDs).fvmDiv(phi, psi)
 225         //- fv::gaussLaplacianScheme<scalar, scalar>(mesh, CDs, snGrads)
 226         //.fvmLaplacian(Dpsif, psi)
 227       - fvm::Sp(Sp, psi)
 228       - Su
 229     );
 230
 231     surfaceScalarField phiBD(psiConvectionDiffusion.flux());
 232
 233     surfaceScalarField& phiCorr = phiPsi;
 234     phiCorr -= phiBD;
 235
 236     for (label i=0; i<maxIter; i++)
 237     {
 238         if (i != 0 && i < 4)
 239         {
 240             allLambda = allCoLambda;
 241         }
 242
 243         limiter
 244         (
 245             allLambda,
 246             rho,
 247             psi,
 248             phiBD,
 249             phiCorr,
 250             Sp.field(),
 251             Su.field(),
 252             psiMax,
 253             psiMin,
 254             nLimiterIter
 255         );
 256
 257         solve
 258         (
 259             psiConvectionDiffusion + fvc::div(lambda*phiCorr),
 260             MULEScontrols
 261         );
 262
 263         scalar maxPsiM1 = gMax(psi.internalField()) - 1.0;
 264         scalar minPsi = gMin(psi.internalField());
 265
 266         scalar unboundedness = max(max(maxPsiM1, 0.0), -min(minPsi, 0.0));
 267
 268         if (unboundedness < maxUnboundedness)
 269         {
 270             break;
 271         }
 272         else
 273         {
 274             Info<< "MULES: max(" << psi.name() << " - 1) = " << maxPsiM1
 275                 << " min(" << psi.name() << ") = " << minPsi << endl;
 276
 277             phiBD = psiConvectionDiffusion.flux();
 278
 279             /*
 280             word gammaScheme("div(phi,gamma)");
 281             word gammarScheme("div(phirb,gamma)");
 282
 283             const surfaceScalarField& phir =
 284                 mesh.lookupObject<surfaceScalarField>("phir");
 285
 286             phiCorr =
 287                 fvc::flux
 288                 (
 289                     phi,
 290                     psi,
 291                     gammaScheme
 292                 )
 293               + fvc::flux
 294                 (
 295                     -fvc::flux(-phir, scalar(1) - psi, gammarScheme),
 296                     psi,
 297                     gammarScheme
 298                 )
 299                 - phiBD;
 300             */
 301         }
 302     }
 303
 304     phiPsi = psiConvectionDiffusion.flux() + lambda*phiCorr;
 305 }
 306
 307
 308 template<class RhoType, class SpType, class SuType>
 309 void Foam::MULES::limiter
 310 (
 311     scalarField& allLambda,
 312     const RhoType& rho,
 313     const volScalarField& psi,
 314     const surfaceScalarField& phiBD,
 315     const surfaceScalarField& phiCorr,
 316     const SpType& Sp,
 317     const SuType& Su,
 318     const scalar psiMax,
 319     const scalar psiMin,
 320     const label nLimiterIter
 321 )
 322 {
 323     const scalarField& psiIf = psi;
 324     const volScalarField::GeometricBoundaryField& psiBf = psi.boundaryField();
 325
 326     const scalarField& psi0 = psi.oldTime();
 327
 328     const fvMesh& mesh = psi.mesh();
 329
 330     const labelUList& owner = mesh.owner();
 331     const labelUList& neighb = mesh.neighbour();
 332     tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
 333     const scalarField& V = tVsc();
 334
 335     const volScalarField& rDeltaT =
 336         mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
 337
 338     const scalarField& phiBDIf = phiBD;
 339     const surfaceScalarField::GeometricBoundaryField& phiBDBf =
 340         phiBD.boundaryField();
 341
 342     const scalarField& phiCorrIf = phiCorr;
 343     const surfaceScalarField::GeometricBoundaryField& phiCorrBf =
 344         phiCorr.boundaryField();
 345
 346     slicedSurfaceScalarField lambda
 347     (
 348         IOobject
 349         (
 350             "lambda",
 351             mesh.time().timeName(),
 352             mesh,
 353             IOobject::NO_READ,
 354             IOobject::NO_WRITE,
 355             false
 356         ),
 357         mesh,
 358         dimless,
 359         allLambda,
 360         false   // Use slices for the couples
 361     );
 362
 363     scalarField& lambdaIf = lambda;
 364     surfaceScalarField::GeometricBoundaryField& lambdaBf =
 365         lambda.boundaryField();
 366
 367     scalarField psiMaxn(psiIf.size(), psiMin);
 368     scalarField psiMinn(psiIf.size(), psiMax);
 369
 370     scalarField sumPhiBD(psiIf.size(), 0.0);
 371
 372     scalarField sumPhip(psiIf.size(), VSMALL);
 373     scalarField mSumPhim(psiIf.size(), VSMALL);
 374
 375     forAll(phiCorrIf, facei)
 376     {
 377         label own = owner[facei];
 378         label nei = neighb[facei];
 379
 380         psiMaxn[own] = max(psiMaxn[own], psiIf[nei]);
 381         psiMinn[own] = min(psiMinn[own], psiIf[nei]);
 382
 383         psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]);
 384         psiMinn[nei] = min(psiMinn[nei], psiIf[own]);
 385
 386         sumPhiBD[own] += phiBDIf[facei];
 387         sumPhiBD[nei] -= phiBDIf[facei];
 388
 389         scalar phiCorrf = phiCorrIf[facei];
 390
 391         if (phiCorrf > 0.0)
 392         {
 393             sumPhip[own] += phiCorrf;
 394             mSumPhim[nei] += phiCorrf;
 395         }
 396         else
 397         {
 398             mSumPhim[own] -= phiCorrf;
 399             sumPhip[nei] -= phiCorrf;
 400         }
 401     }
 402
 403     forAll(phiCorrBf, patchi)
 404     {
 405         const fvPatchScalarField& psiPf = psiBf[patchi];
 406         const scalarField& phiBDPf = phiBDBf[patchi];
 407         const scalarField& phiCorrPf = phiCorrBf[patchi];
 408
 409         const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
 410
 411         if (psiPf.coupled())
 412         {
 413             scalarField psiPNf(psiPf.patchNeighbourField());
 414
 415             forAll(phiCorrPf, pFacei)
 416             {
 417                 label pfCelli = pFaceCells[pFacei];
 418
 419                 psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]);
 420                 psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]);
 421             }
 422         }
 423         else
 424         {
 425             forAll(phiCorrPf, pFacei)
 426             {
 427                 label pfCelli = pFaceCells[pFacei];
 428
 429                 psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]);
 430                 psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]);
 431             }
 432         }
 433
 434         forAll(phiCorrPf, pFacei)
 435         {
 436             label pfCelli = pFaceCells[pFacei];
 437
 438             sumPhiBD[pfCelli] += phiBDPf[pFacei];
 439
 440             scalar phiCorrf = phiCorrPf[pFacei];
 441
 442             if (phiCorrf > 0.0)
 443             {
 444                 sumPhip[pfCelli] += phiCorrf;
 445             }
 446             else
 447             {
 448                 mSumPhim[pfCelli] -= phiCorrf;
 449             }
 450         }
 451     }
 452
 453     psiMaxn = min(psiMaxn, psiMax);
 454     psiMinn = max(psiMinn, psiMin);
 455
 456     //scalar smooth = 0.5;
 457     //psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax);
 458     //psiMinn = max((1.0 - smooth)*psiIf + smooth*psiMinn, psiMin);
 459
 460     if (mesh.moving())
 461     {
 462         tmp<volScalarField::DimensionedInternalField> V0 = mesh.Vsc0();
 463
 464         psiMaxn =
 465             V*((rho*rDeltaT - Sp)*psiMaxn - Su)
 466           - (V0()*rDeltaT)*rho.oldTime()*psi0
 467           + sumPhiBD;
 468
 469         psiMinn =
 470             V*(Su - (rho*rDeltaT - Sp)*psiMinn)
 471           + (V0*rDeltaT)*rho.oldTime()*psi0
 472           - sumPhiBD;
 473     }
 474     else
 475     {
 476         psiMaxn =
 477             V*((rho*rDeltaT - Sp)*psiMaxn - (rho.oldTime()*rDeltaT)*psi0 - Su)
 478           + sumPhiBD;
 479
 480         psiMinn =
 481             V*((rho*rDeltaT)*psi0 - (rho.oldTime()*rDeltaT - Sp)*psiMinn + Su)
 482           - sumPhiBD;
 483     }
 484
 485     scalarField sumlPhip(psiIf.size());
 486     scalarField mSumlPhim(psiIf.size());
 487
 488     for(int j=0; j<nLimiterIter; j++)
 489     {
 490         sumlPhip = 0.0;
 491         mSumlPhim = 0.0;
 492
 493         forAll(lambdaIf, facei)
 494         {
 495             label own = owner[facei];
 496             label nei = neighb[facei];
 497
 498             scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
 499
 500             if (lambdaPhiCorrf > 0.0)
 501             {
 502                 sumlPhip[own] += lambdaPhiCorrf;
 503                 mSumlPhim[nei] += lambdaPhiCorrf;
 504             }
 505             else
 506             {
 507                 mSumlPhim[own] -= lambdaPhiCorrf;
 508                 sumlPhip[nei] -= lambdaPhiCorrf;
 509             }
 510         }
 511
 512         forAll(lambdaBf, patchi)
 513         {
 514             scalarField& lambdaPf = lambdaBf[patchi];
 515             const scalarField& phiCorrfPf = phiCorrBf[patchi];
 516
 517             const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
 518
 519             forAll(lambdaPf, pFacei)
 520             {
 521                 label pfCelli = pFaceCells[pFacei];
 522
 523                 scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei];
 524
 525                 if (lambdaPhiCorrf > 0.0)
 526                 {
 527                     sumlPhip[pfCelli] += lambdaPhiCorrf;
 528                 }
 529                 else
 530                 {
 531                     mSumlPhim[pfCelli] -= lambdaPhiCorrf;
 532                 }
 533             }
 534         }
 535
 536         forAll (sumlPhip, celli)
 537         {
 538             sumlPhip[celli] =
 539                 max(min
 540                 (
 541                     (sumlPhip[celli] + psiMaxn[celli])/mSumPhim[celli],
 542                     1.0), 0.0
 543                 );
 544
 545             mSumlPhim[celli] =
 546                 max(min
 547                 (
 548                     (mSumlPhim[celli] + psiMinn[celli])/sumPhip[celli],
 549                     1.0), 0.0
 550                 );
 551         }
 552
 553         const scalarField& lambdam = sumlPhip;
 554         const scalarField& lambdap = mSumlPhim;
 555
 556         forAll(lambdaIf, facei)
 557         {
 558             if (phiCorrIf[facei] > 0.0)
 559             {
 560                 lambdaIf[facei] = min
 561                 (
 562                     lambdaIf[facei],
 563                     min(lambdap[owner[facei]], lambdam[neighb[facei]])
 564                 );
 565             }
 566             else
 567             {
 568                 lambdaIf[facei] = min
 569                 (
 570                     lambdaIf[facei],
 571                     min(lambdam[owner[facei]], lambdap[neighb[facei]])
 572                 );
 573             }
 574         }
 575
 576
 577         forAll(lambdaBf, patchi)
 578         {
 579             fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
 580             const scalarField& phiCorrfPf = phiCorrBf[patchi];
 581
 582             const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
 583
 584             forAll(lambdaPf, pFacei)
 585             {
 586                 label pfCelli = pFaceCells[pFacei];
 587
 588                 if (phiCorrfPf[pFacei] > 0.0)
 589                 {
 590                     lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdap[pfCelli]);
 591                 }
 592                 else
 593                 {
 594                     lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdam[pfCelli]);
 595                 }
 596             }
 597         }
 598
 599         syncTools::syncFaceList(mesh, allLambda, minEqOp<scalar>());
 600     }
 601 }
 602
 603
 604 // ************************************************************************* //