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 =
 146        mesh.solverDict(psi.name()).subDict("MULESImplicit");
 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             mesh.time().deltaT()*mesh.surfaceInterpolation::deltaCoeffs()
 173            *mag(phi)/mesh.magSf();
 174
 175         slicedSurfaceScalarField CoLambda
 176         (
 177             IOobject
 178             (
 179                 "CoLambda",
 180                 mesh.time().timeName(),
 181                 mesh,
 182                 IOobject::NO_READ,
 183                 IOobject::NO_WRITE,
 184                 false
 185             ),
 186             mesh,
 187             dimless,
 188             allCoLambda,
 189             false   // Use slices for the couples
 190         );
 191
 192         CoLambda == 1.0/max(CoCoeff*Cof, scalar(1));
 193     }
 194
 195     scalarField allLambda(allCoLambda);
 196     //scalarField allLambda(mesh.nFaces(), 1.0);
 197
 198     slicedSurfaceScalarField lambda
 199     (
 200         IOobject
 201         (
 202             "lambda",
 203             mesh.time().timeName(),
 204             mesh,
 205             IOobject::NO_READ,
 206             IOobject::NO_WRITE,
 207             false
 208         ),
 209         mesh,
 210         dimless,
 211         allLambda,
 212         false   // Use slices for the couples
 213     );
 214
 215     linear<scalar> CDs(mesh);
 216     upwind<scalar> UDs(mesh, phi);
 217     //fv::uncorrectedSnGrad<scalar> snGrads(mesh);
 218
 219     fvScalarMatrix psiConvectionDiffusion
 220     (
 221         fvm::ddt(rho, psi)
 222       + fv::gaussConvectionScheme<scalar>(mesh, phi, UDs).fvmDiv(phi, psi)
 223         //- fv::gaussLaplacianScheme<scalar, scalar>(mesh, CDs, snGrads)
 224         //.fvmLaplacian(Dpsif, psi)
 225       - fvm::Sp(Sp, psi)
 226       - Su
 227     );
 228
 229     surfaceScalarField phiBD = psiConvectionDiffusion.flux();
 230
 231     surfaceScalarField& phiCorr = phiPsi;
 232     phiCorr -= phiBD;
 233
 234     for (label i=0; i<maxIter; i++)
 235     {
 236         if (i != 0 && i < 4)
 237         {
 238             allLambda = allCoLambda;
 239         }
 240
 241         limiter
 242         (
 243             allLambda,
 244             rho,
 245             psi,
 246             phiBD,
 247             phiCorr,
 248             Sp.field(),
 249             Su.field(),
 250             psiMax,
 251             psiMin,
 252             nLimiterIter
 253         );
 254
 255         solve
 256         (
 257             psiConvectionDiffusion + fvc::div(lambda*phiCorr),
 258             MULEScontrols.lookup("solver")
 259         );
 260
 261         scalar maxPsiM1 = gMax(psi.internalField()) - 1.0;
 262         scalar minPsi = gMin(psi.internalField());
 263
 264         scalar unboundedness = max(max(maxPsiM1, 0.0), -min(minPsi, 0.0));
 265
 266         if (unboundedness < maxUnboundedness)
 267         {
 268             break;
 269         }
 270         else
 271         {
 272             Info<< "MULES: max(" << psi.name() << " - 1) = " << maxPsiM1
 273                 << " min(" << psi.name() << ") = " << minPsi << endl;
 274
 275             phiBD = psiConvectionDiffusion.flux();
 276
 277             /*
 278             word gammaScheme("div(phi,gamma)");
 279             word gammarScheme("div(phirb,gamma)");
 280
 281             const surfaceScalarField& phir =
 282                 mesh.lookupObject<surfaceScalarField>("phir");
 283
 284             phiCorr =
 285                 fvc::flux
 286                 (
 287                     phi,
 288                     psi,
 289                     gammaScheme
 290                 )
 291               + fvc::flux
 292                 (
 293                     -fvc::flux(-phir, scalar(1) - psi, gammarScheme),
 294                     psi,
 295                     gammarScheme
 296                 )
 297                 - phiBD;
 298             */
 299         }
 300     }
 301
 302     phiPsi = psiConvectionDiffusion.flux() + lambda*phiCorr;
 303 }
 304
 305
 306 template<class RhoType, class SpType, class SuType>
 307 void Foam::MULES::limiter
 308 (
 309     scalarField& allLambda,
 310     const RhoType& rho,
 311     const volScalarField& psi,
 312     const surfaceScalarField& phiBD,
 313     const surfaceScalarField& phiCorr,
 314     const SpType& Sp,
 315     const SuType& Su,
 316     const scalar psiMax,
 317     const scalar psiMin,
 318     const label nLimiterIter
 319 )
 320 {
 321     const scalarField& psiIf = psi;
 322     const volScalarField::GeometricBoundaryField& psiBf = psi.boundaryField();
 323
 324     const scalarField& psi0 = psi.oldTime();
 325
 326     const fvMesh& mesh = psi.mesh();
 327
 328     const unallocLabelList& owner = mesh.owner();
 329     const unallocLabelList& neighb = mesh.neighbour();
 330     const scalarField& V = mesh.V();
 331     const scalar deltaT = mesh.time().deltaT().value();
 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             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         psiMaxn =
 458             V*((rho/deltaT - Sp)*psiMaxn - Su)
 459           - (mesh.V0()/deltaT)*rho.oldTime()*psi0
 460           + sumPhiBD;
 461
 462         psiMinn =
 463             V*(Su - (rho/deltaT - Sp)*psiMinn)
 464           + (mesh.V0()/deltaT)*rho.oldTime()*psi0
 465           - sumPhiBD;
 466     }
 467     else
 468     {
 469         psiMaxn =
 470             V*((rho/deltaT - Sp)*psiMaxn - (rho.oldTime()/deltaT)*psi0 - Su)
 471           + sumPhiBD;
 472
 473         psiMinn =
 474             V*((rho/deltaT)*psi0 - (rho.oldTime()/deltaT - Sp)*psiMinn + Su)
 475           - sumPhiBD;
 476     }
 477
 478     scalarField sumlPhip(psiIf.size());
 479     scalarField mSumlPhim(psiIf.size());
 480
 481     for(int j=0; j<nLimiterIter; j++)
 482     {
 483         sumlPhip = 0.0;
 484         mSumlPhim = 0.0;
 485
 486         forAll(lambdaIf, facei)
 487         {
 488             label own = owner[facei];
 489             label nei = neighb[facei];
 490
 491             scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
 492
 493             if (lambdaPhiCorrf > 0.0)
 494             {
 495                 sumlPhip[own] += lambdaPhiCorrf;
 496                 mSumlPhim[nei] += lambdaPhiCorrf;
 497             }
 498             else
 499             {
 500                 mSumlPhim[own] -= lambdaPhiCorrf;
 501                 sumlPhip[nei] -= lambdaPhiCorrf;
 502             }
 503         }
 504
 505         forAll(lambdaBf, patchi)
 506         {
 507             scalarField& lambdaPf = lambdaBf[patchi];
 508             const scalarField& phiCorrfPf = phiCorrBf[patchi];
 509
 510             const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
 511
 512             forAll(lambdaPf, pFacei)
 513             {
 514                 label pfCelli = pFaceCells[pFacei];
 515
 516                 scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei];
 517
 518                 if (lambdaPhiCorrf > 0.0)
 519                 {
 520                     sumlPhip[pfCelli] += lambdaPhiCorrf;
 521                 }
 522                 else
 523                 {
 524                     mSumlPhim[pfCelli] -= lambdaPhiCorrf;
 525                 }
 526             }
 527         }
 528
 529         forAll (sumlPhip, celli)
 530         {
 531             sumlPhip[celli] =
 532                 max(min
 533                 (
 534                     (sumlPhip[celli] + psiMaxn[celli])/mSumPhim[celli],
 535                     1.0), 0.0
 536                 );
 537
 538             mSumlPhim[celli] =
 539                 max(min
 540                 (
 541                     (mSumlPhim[celli] + psiMinn[celli])/sumPhip[celli],
 542                     1.0), 0.0
 543                 );
 544         }
 545
 546         const scalarField& lambdam = sumlPhip;
 547         const scalarField& lambdap = mSumlPhim;
 548
 549         forAll(lambdaIf, facei)
 550         {
 551             if (phiCorrIf[facei] > 0.0)
 552             {
 553                 lambdaIf[facei] = min
 554                 (
 555                     lambdaIf[facei],
 556                     min(lambdap[owner[facei]], lambdam[neighb[facei]])
 557                 );
 558             }
 559             else
 560             {
 561                 lambdaIf[facei] = min
 562                 (
 563                     lambdaIf[facei],
 564                     min(lambdam[owner[facei]], lambdap[neighb[facei]])
 565                 );
 566             }
 567         }
 568
 569
 570         forAll(lambdaBf, patchi)
 571         {
 572             fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
 573             const scalarField& phiCorrfPf = phiCorrBf[patchi];
 574
 575             const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
 576
 577             forAll(lambdaPf, pFacei)
 578             {
 579                 label pfCelli = pFaceCells[pFacei];
 580
 581                 if (phiCorrfPf[pFacei] > 0.0)
 582                 {
 583                     lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdap[pfCelli]);
 584                 }
 585                 else
 586                 {
 587                     lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdam[pfCelli]);
 588                 }
 589             }
 590         }
 591
 592         syncTools::syncFaceList(mesh, allLambda, minEqOp<scalar>(), false);
 593     }
 594 }
 595
 596
 597 // ************************************************************************* //