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) 1991-2009 OpenCFD Ltd.
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "MULES.H"
  28 #include "upwind.H"
  29 #include "uncorrectedSnGrad.H"
  30 #include "gaussConvectionScheme.H"
  31 #include "gaussLaplacianScheme.H"
  32 #include "uncorrectedSnGrad.H"
  33 #include "surfaceInterpolate.H"
  34 #include "fvcSurfaceIntegrate.H"
  35 #include "slicedSurfaceFields.H"
  36 #include "syncTools.H"
  37
  38 #include "fvCFD.H"
  39
  40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  41
  42 template<class RhoType, class SpType, class SuType>
  43 void Foam::MULES::explicitSolve
  44 (
  45     const RhoType& rho,
  46     volScalarField& psi,
  47     const surfaceScalarField& phi,
  48     surfaceScalarField& phiPsi,
  49     const SpType& Sp,
  50     const SuType& Su,
  51     const scalar psiMax,
  52     const scalar psiMin
  53 )
  54 {
  55     Info<< "MULES: Solving for " << psi.name() << endl;
  56
  57     const fvMesh& mesh = psi.mesh();
  58     psi.correctBoundaryConditions();
  59
  60     surfaceScalarField phiBD = upwind<scalar>(psi.mesh(), phi).flux(psi);
  61
  62     surfaceScalarField& phiCorr = phiPsi;
  63     phiCorr -= phiBD;
  64
  65     scalarField allLambda(mesh.nFaces(), 1.0);
  66
  67     slicedSurfaceScalarField lambda
  68     (
  69         IOobject
  70         (
  71             "lambda",
  72             mesh.time().timeName(),
  73             mesh,
  74             IOobject::NO_READ,
  75             IOobject::NO_WRITE,
  76             false
  77         ),
  78         mesh,
  79         dimless,
  80         allLambda,
  81         false   // Use slices for the couples
  82     );
  83
  84     limiter
  85     (
  86         allLambda,
  87         rho,
  88         psi,
  89         phiBD,
  90         phiCorr,
  91         Sp.field(),
  92         Su.field(),
  93         psiMax,
  94         psiMin,
  95         3
  96     );
  97
  98     phiPsi = phiBD + lambda*phiCorr;
  99
 100     scalarField& psiIf = psi;
 101     const scalarField& psi0 = psi.oldTime();
 102     const scalar deltaT = mesh.time().deltaT().value();
 103
 104     psiIf = 0.0;
 105     fvc::surfaceIntegrate(psiIf, phiPsi);
 106
 107     if (mesh.moving())
 108     {
 109         psiIf =
 110         (
 111             mesh.V0()*rho.oldTime()*psi0/(deltaT*mesh.V())
 112           + Su.field()
 113           - psiIf
 114         )/(rho/deltaT - Sp.field());
 115     }
 116     else
 117     {
 118         psiIf =
 119         (
 120             rho.oldTime()*psi0/deltaT
 121           + Su.field()
 122           - psiIf
 123         )/(rho/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         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.field(),
 248             Su.field(),
 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 unallocLabelList& owner = mesh.owner();
 328     const unallocLabelList& neighb = mesh.neighbour();
 329     const scalarField& V = mesh.V();
 330     const scalar deltaT = mesh.time().deltaT().value();
 331
 332     const scalarField& phiBDIf = phiBD;
 333     const surfaceScalarField::GeometricBoundaryField& phiBDBf =
 334         phiBD.boundaryField();
 335
 336     const scalarField& phiCorrIf = phiCorr;
 337     const surfaceScalarField::GeometricBoundaryField& phiCorrBf =
 338         phiCorr.boundaryField();
 339
 340     slicedSurfaceScalarField lambda
 341     (
 342         IOobject
 343         (
 344             "lambda",
 345             mesh.time().timeName(),
 346             mesh,
 347             IOobject::NO_READ,
 348             IOobject::NO_WRITE,
 349             false
 350         ),
 351         mesh,
 352         dimless,
 353         allLambda,
 354         false   // Use slices for the couples
 355     );
 356
 357     scalarField& lambdaIf = lambda;
 358     surfaceScalarField::GeometricBoundaryField& lambdaBf =
 359         lambda.boundaryField();
 360
 361     scalarField psiMaxn(psiIf.size(), psiMin);
 362     scalarField psiMinn(psiIf.size(), psiMax);
 363
 364     scalarField sumPhiBD(psiIf.size(), 0.0);
 365
 366     scalarField sumPhip(psiIf.size(), VSMALL);
 367     scalarField mSumPhim(psiIf.size(), VSMALL);
 368
 369     forAll(phiCorrIf, facei)
 370     {
 371         label own = owner[facei];
 372         label nei = neighb[facei];
 373
 374         psiMaxn[own] = max(psiMaxn[own], psiIf[nei]);
 375         psiMinn[own] = min(psiMinn[own], psiIf[nei]);
 376
 377         psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]);
 378         psiMinn[nei] = min(psiMinn[nei], psiIf[own]);
 379
 380         sumPhiBD[own] += phiBDIf[facei];
 381         sumPhiBD[nei] -= phiBDIf[facei];
 382
 383         scalar phiCorrf = phiCorrIf[facei];
 384
 385         if (phiCorrf > 0.0)
 386         {
 387             sumPhip[own] += phiCorrf;
 388             mSumPhim[nei] += phiCorrf;
 389         }
 390         else
 391         {
 392             mSumPhim[own] -= phiCorrf;
 393             sumPhip[nei] -= phiCorrf;
 394         }
 395     }
 396
 397     forAll(phiCorrBf, patchi)
 398     {
 399         const fvPatchScalarField& psiPf = psiBf[patchi];
 400         const scalarField& phiBDPf = phiBDBf[patchi];
 401         const scalarField& phiCorrPf = phiCorrBf[patchi];
 402
 403         const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
 404
 405         if (psiPf.coupled())
 406         {
 407             scalarField psiPNf = psiPf.patchNeighbourField();
 408
 409             forAll(phiCorrPf, pFacei)
 410             {
 411                 label pfCelli = pFaceCells[pFacei];
 412
 413                 psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]);
 414                 psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]);
 415             }
 416         }
 417         else
 418         {
 419             forAll(phiCorrPf, pFacei)
 420             {
 421                 label pfCelli = pFaceCells[pFacei];
 422
 423                 psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]);
 424                 psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]);
 425             }
 426         }
 427
 428         forAll(phiCorrPf, pFacei)
 429         {
 430             label pfCelli = pFaceCells[pFacei];
 431
 432             sumPhiBD[pfCelli] += phiBDPf[pFacei];
 433
 434             scalar phiCorrf = phiCorrPf[pFacei];
 435
 436             if (phiCorrf > 0.0)
 437             {
 438                 sumPhip[pfCelli] += phiCorrf;
 439             }
 440             else
 441             {
 442                 mSumPhim[pfCelli] -= phiCorrf;
 443             }
 444         }
 445     }
 446
 447     psiMaxn = min(psiMaxn, psiMax);
 448     psiMinn = max(psiMinn, psiMin);
 449
 450     //scalar smooth = 0.5;
 451     //psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax);
 452     //psiMinn = max((1.0 - smooth)*psiIf + smooth*psiMinn, psiMin);
 453
 454     if (mesh.moving())
 455     {
 456         psiMaxn =
 457             V*((rho/deltaT - Sp)*psiMaxn - Su)
 458           - (mesh.V0()/deltaT)*rho.oldTime()*psi0
 459           + sumPhiBD;
 460
 461         psiMinn =
 462             V*(Su - (rho/deltaT - Sp)*psiMinn)
 463           + (mesh.V0()/deltaT)*rho.oldTime()*psi0
 464           - sumPhiBD;
 465     }
 466     else
 467     {
 468         psiMaxn =
 469             V*((rho/deltaT - Sp)*psiMaxn - (rho.oldTime()/deltaT)*psi0 - Su)
 470           + sumPhiBD;
 471
 472         psiMinn =
 473             V*((rho/deltaT)*psi0 - (rho.oldTime()/deltaT - Sp)*psiMinn + Su)
 474           - sumPhiBD;
 475     }
 476
 477     scalarField sumlPhip(psiIf.size());
 478     scalarField mSumlPhim(psiIf.size());
 479
 480     for(int j=0; j<nLimiterIter; j++)
 481     {
 482         sumlPhip = 0.0;
 483         mSumlPhim = 0.0;
 484
 485         forAll(lambdaIf, facei)
 486         {
 487             label own = owner[facei];
 488             label nei = neighb[facei];
 489
 490             scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
 491
 492             if (lambdaPhiCorrf > 0.0)
 493             {
 494                 sumlPhip[own] += lambdaPhiCorrf;
 495                 mSumlPhim[nei] += lambdaPhiCorrf;
 496             }
 497             else
 498             {
 499                 mSumlPhim[own] -= lambdaPhiCorrf;
 500                 sumlPhip[nei] -= lambdaPhiCorrf;
 501             }
 502         }
 503
 504         forAll(lambdaBf, patchi)
 505         {
 506             scalarField& lambdaPf = lambdaBf[patchi];
 507             const scalarField& phiCorrfPf = phiCorrBf[patchi];
 508
 509             const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
 510
 511             forAll(lambdaPf, pFacei)
 512             {
 513                 label pfCelli = pFaceCells[pFacei];
 514
 515                 scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei];
 516
 517                 if (lambdaPhiCorrf > 0.0)
 518                 {
 519                     sumlPhip[pfCelli] += lambdaPhiCorrf;
 520                 }
 521                 else
 522                 {
 523                     mSumlPhim[pfCelli] -= lambdaPhiCorrf;
 524                 }
 525             }
 526         }
 527
 528         forAll (sumlPhip, celli)
 529         {
 530             sumlPhip[celli] =
 531                 max(min
 532                 (
 533                     (sumlPhip[celli] + psiMaxn[celli])/mSumPhim[celli],
 534                     1.0), 0.0
 535                 );
 536
 537             mSumlPhim[celli] =
 538                 max(min
 539                 (
 540                     (mSumlPhim[celli] + psiMinn[celli])/sumPhip[celli],
 541                     1.0), 0.0
 542                 );
 543         }
 544
 545         const scalarField& lambdam = sumlPhip;
 546         const scalarField& lambdap = mSumlPhim;
 547
 548         forAll(lambdaIf, facei)
 549         {
 550             if (phiCorrIf[facei] > 0.0)
 551             {
 552                 lambdaIf[facei] = min
 553                 (
 554                     lambdaIf[facei],
 555                     min(lambdap[owner[facei]], lambdam[neighb[facei]])
 556                 );
 557             }
 558             else
 559             {
 560                 lambdaIf[facei] = min
 561                 (
 562                     lambdaIf[facei],
 563                     min(lambdam[owner[facei]], lambdap[neighb[facei]])
 564                 );
 565             }
 566         }
 567
 568
 569         forAll(lambdaBf, patchi)
 570         {
 571             fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
 572             const scalarField& phiCorrfPf = phiCorrBf[patchi];
 573
 574             const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
 575
 576             forAll(lambdaPf, pFacei)
 577             {
 578                 label pfCelli = pFaceCells[pFacei];
 579
 580                 if (phiCorrfPf[pFacei] > 0.0)
 581                 {
 582                     lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdap[pfCelli]);
 583                 }
 584                 else
 585                 {
 586                     lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdam[pfCelli]);
 587                 }
 588             }
 589         }
 590
 591         syncTools::syncFaceList(mesh, allLambda, minEqOp<scalar>(), false);
 592     }
 593 }
 594
 595
 596 // ************************************************************************* //