src/thermophysicalModels/radiation/radiationModel/viewFactor/viewFactor.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "viewFactor.H"
  27 #include "surfaceFields.H"
  28 #include "addToRunTimeSelectionTable.H"
  29 #include "mathematicalConstants.H"
  30 #include "radiationConstants.H"
  31 #include "greyDiffusiveViewFactorFixedValueFvPatchScalarField.H"
  32 #include "typeInfo.H"
  33
  34
  35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  36
  37 namespace Foam
  38 {
  39     namespace radiation
  40     {
  41         defineTypeNameAndDebug(viewFactor, 0);
  42
  43         addToRunTimeSelectionTable
  44         (
  45             radiationModel,
  46             viewFactor,
  47             dictionary
  48         );
  49     }
  50 }
  51
  52
  53 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  54
  55 void Foam::radiation::viewFactor::initialise()
  56 {
  57     const polyBoundaryMesh& coarsePatches = coarseMesh_.boundaryMesh();
  58     const volScalarField::GeometricBoundaryField& Qrp = Qr_.boundaryField();
  59
  60     label count = 0;
  61     forAll(Qrp, patchI)
  62     {
  63         //const polyPatch& pp = mesh_.boundaryMesh()[patchI];
  64         const fvPatchScalarField& QrPatchI = Qrp[patchI];
  65
  66         if ((isA<fixedValueFvPatchScalarField>(QrPatchI)))
  67         {
  68             selectedPatches_[count] = QrPatchI.patch().index();
  69             nLocalCoarseFaces_ += coarsePatches[patchI].size();
  70             count++;
  71         }
  72     }
  73
  74     selectedPatches_.resize(count--);
  75
  76     if (debug)
  77     {
  78         Pout<< "Selected patches:" << selectedPatches_ << endl;
  79         Pout<< "Number of coarse faces:" << nLocalCoarseFaces_ << endl;
  80     }
  81
  82     totalNCoarseFaces_ = nLocalCoarseFaces_;
  83     reduce(totalNCoarseFaces_, sumOp<label>());
  84
  85     if (Pstream::master())
  86     {
  87         Info<< "Total number of clusters : " << totalNCoarseFaces_ << endl;
  88     }
  89
  90     labelListIOList subMap
  91     (
  92         IOobject
  93         (
  94             "subMap",
  95             mesh_.facesInstance(),
  96             mesh_,
  97             IOobject::MUST_READ,
  98             IOobject::NO_WRITE,
  99             false
 100         )
 101     );
 102
 103     labelListIOList constructMap
 104     (
 105         IOobject
 106         (
 107             "constructMap",
 108             mesh_.facesInstance(),
 109             mesh_,
 110             IOobject::MUST_READ,
 111             IOobject::NO_WRITE,
 112             false
 113         )
 114     );
 115
 116     IOList<label> consMapDim
 117     (
 118         IOobject
 119         (
 120             "constructMapDim",
 121             mesh_.facesInstance(),
 122             mesh_,
 123             IOobject::MUST_READ,
 124             IOobject::NO_WRITE,
 125             false
 126         )
 127     );
 128
 129     map_.reset
 130     (
 131         new mapDistribute
 132         (
 133             consMapDim[0],
 134             Xfer<labelListList>(subMap),
 135             Xfer<labelListList>(constructMap)
 136         )
 137     );
 138
 139     scalarListIOList FmyProc
 140     (
 141         IOobject
 142         (
 143             "F",
 144             mesh_.facesInstance(),
 145             mesh_,
 146             IOobject::MUST_READ,
 147             IOobject::NO_WRITE,
 148             false
 149         )
 150     );
 151
 152     labelListIOList globalFaceFaces
 153     (
 154         IOobject
 155         (
 156             "globalFaceFaces",
 157             mesh_.facesInstance(),
 158             mesh_,
 159             IOobject::MUST_READ,
 160             IOobject::NO_WRITE,
 161             false
 162         )
 163     );
 164
 165     List<labelListList> globalFaceFacesProc(Pstream::nProcs());
 166     globalFaceFacesProc[Pstream::myProcNo()] = globalFaceFaces;
 167     Pstream::gatherList(globalFaceFacesProc);
 168
 169     List<scalarListList> F(Pstream::nProcs());
 170     F[Pstream::myProcNo()] = FmyProc;
 171     Pstream::gatherList(F);
 172
 173     globalIndex globalNumbering(nLocalCoarseFaces_);
 174
 175     if (Pstream::master())
 176     {
 177         Fmatrix_.reset
 178         (
 179             new scalarSquareMatrix(totalNCoarseFaces_, 0.0)
 180         );
 181
 182         Info<< "Insert elements in the matrix..." << endl;
 183
 184         for (label procI = 0; procI < Pstream::nProcs(); procI++)
 185         {
 186             insertMatrixElements
 187             (
 188                 globalNumbering,
 189                 procI,
 190                 globalFaceFacesProc[procI],
 191                 F[procI],
 192                 Fmatrix_()
 193             );
 194         }
 195
 196
 197         bool smoothing = readBool(coeffs_.lookup("smoothing"));
 198         if (smoothing)
 199         {
 200             Info<< "Smoothing the matrix..." << endl;
 201
 202             for (label i=0; i<totalNCoarseFaces_; i++)
 203             {
 204                 scalar sumF = 0.0;
 205                 for (label j=0; j<totalNCoarseFaces_; j++)
 206                 {
 207                     sumF += Fmatrix_()[i][j];
 208                 }
 209                 scalar delta = 1.0 - sumF;
 210                 for (label j=0; j<totalNCoarseFaces_; j++)
 211                 {
 212                     Fmatrix_()[i][j] *= (1.0 - delta/(sumF + 0.001));
 213                 }
 214             }
 215         }
 216
 217         constEmissivity_ = readBool(coeffs_.lookup("constantEmissivity"));
 218         if (constEmissivity_)
 219         {
 220             CLU_.reset
 221             (
 222                 new scalarSquareMatrix
 223                 (
 224                     totalNCoarseFaces_,
 225                     0.0
 226                 )
 227             );
 228
 229             pivotIndices_.setSize(CLU_().n());
 230         }
 231     }
 232 }
 233
 234
 235 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 236
 237 Foam::radiation::viewFactor::viewFactor(const volScalarField& T)
 238 :
 239     radiationModel(typeName, T),
 240     finalAgglom_
 241     (
 242         IOobject
 243         (
 244             "finalAgglom",
 245             mesh_.facesInstance(),
 246             mesh_,
 247             IOobject::MUST_READ,
 248             IOobject::NO_WRITE,
 249             false
 250         )
 251     ),
 252     map_(),
 253     coarseMesh_
 254     (
 255         IOobject
 256         (
 257             mesh_.name(),
 258             mesh_.polyMesh::instance(),
 259             mesh_.time(),
 260             IOobject::NO_READ,
 261             IOobject::NO_WRITE
 262         ),
 263         mesh_,
 264         finalAgglom_
 265     ),
 266     Qr_
 267     (
 268         IOobject
 269         (
 270             "Qr",
 271             mesh_.time().timeName(),
 272             mesh_,
 273             IOobject::MUST_READ,
 274             IOobject::AUTO_WRITE
 275         ),
 276         mesh_
 277     ),
 278     Fmatrix_(),
 279     CLU_(),
 280     selectedPatches_(mesh_.boundary().size(), -1),
 281     totalNCoarseFaces_(0),
 282     nLocalCoarseFaces_(0),
 283     constEmissivity_(false),
 284     iterCounter_(0),
 285     pivotIndices_(0)
 286 {
 287     initialise();
 288 }
 289
 290
 291 Foam::radiation::viewFactor::viewFactor
 292 (
 293     const dictionary& dict,
 294     const volScalarField& T
 295 )
 296 :
 297     radiationModel(typeName, T),
 298     finalAgglom_
 299     (
 300         IOobject
 301         (
 302             "finalAgglom",
 303             mesh_.facesInstance(),
 304             mesh_,
 305             IOobject::MUST_READ,
 306             IOobject::NO_WRITE,
 307             false
 308         )
 309     ),
 310     map_(),
 311     coarseMesh_
 312     (
 313         IOobject
 314         (
 315             mesh_.name(),
 316             mesh_.polyMesh::instance(),
 317             mesh_.time(),
 318             IOobject::NO_READ,
 319             IOobject::NO_WRITE
 320         ),
 321         mesh_,
 322         finalAgglom_
 323     ),
 324     Qr_
 325     (
 326         IOobject
 327         (
 328             "Qr",
 329             mesh_.time().timeName(),
 330             mesh_,
 331             IOobject::MUST_READ,
 332             IOobject::AUTO_WRITE
 333         ),
 334         mesh_
 335     ),
 336     Fmatrix_(),
 337     CLU_(),
 338     selectedPatches_(mesh_.boundary().size(), -1),
 339     totalNCoarseFaces_(0),
 340     nLocalCoarseFaces_(0),
 341     constEmissivity_(false),
 342     iterCounter_(0),
 343     pivotIndices_(0)
 344 {
 345     initialise();
 346 }
 347
 348
 349 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 350
 351 Foam::radiation::viewFactor::~viewFactor()
 352 {}
 353
 354
 355 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 356
 357 bool Foam::radiation::viewFactor::read()
 358 {
 359     if (radiationModel::read())
 360     {
 361         return true;
 362     }
 363     else
 364     {
 365         return false;
 366     }
 367 }
 368
 369
 370 void Foam::radiation::viewFactor::insertMatrixElements
 371 (
 372     const globalIndex& globalNumbering,
 373     const label procI,
 374     const labelListList& globalFaceFaces,
 375     const scalarListList& viewFactors,
 376     scalarSquareMatrix& Fmatrix
 377 )
 378 {
 379     forAll(viewFactors, faceI)
 380     {
 381         const scalarList& vf = viewFactors[faceI];
 382         const labelList& globalFaces = globalFaceFaces[faceI];
 383
 384         label globalI = globalNumbering.toGlobal(procI, faceI);
 385         forAll(globalFaces, i)
 386         {
 387             Fmatrix[globalI][globalFaces[i]] = vf[i];
 388         }
 389     }
 390 }
 391
 392
 393 void Foam::radiation::viewFactor::calculate()
 394 {
 395     // Store previous iteration
 396     Qr_.storePrevIter();
 397
 398     scalarField compactCoarseT(map_->constructSize(), 0.0);
 399     scalarField compactCoarseE(map_->constructSize(), 0.0);
 400     scalarField compactCoarseHo(map_->constructSize(), 0.0);
 401
 402     globalIndex globalNumbering(nLocalCoarseFaces_);
 403
 404     // Fill local averaged(T), emissivity(E) and external heatFlux(Ho)
 405     DynamicList<scalar> localCoarseTave(nLocalCoarseFaces_);
 406     DynamicList<scalar> localCoarseEave(nLocalCoarseFaces_);
 407     DynamicList<scalar> localCoarseHoave(nLocalCoarseFaces_);
 408
 409     forAll(selectedPatches_, i)
 410     {
 411         label patchID = selectedPatches_[i];
 412
 413         const scalarField& Tp = T_.boundaryField()[patchID];
 414         const scalarField& sf = mesh_.magSf().boundaryField()[patchID];
 415
 416         fvPatchScalarField& QrPatch = Qr_.boundaryField()[patchID];
 417
 418         greyDiffusiveViewFactorFixedValueFvPatchScalarField& Qrp =
 419             refCast
 420             <
 421                 greyDiffusiveViewFactorFixedValueFvPatchScalarField
 422             >(QrPatch);
 423
 424         const scalarList eb = Qrp.emissivity();
 425
 426         const scalarList& Hoi = Qrp.Qro();
 427
 428         const polyPatch& pp = coarseMesh_.boundaryMesh()[patchID];
 429         const labelList& coarsePatchFace = coarseMesh_.patchFaceMap()[patchID];
 430
 431         scalarList Tave(pp.size(), 0.0);
 432         scalarList Eave(Tave.size(), 0.0);
 433         scalarList Hoiave(Tave.size(), 0.0);
 434
 435         if (pp.size() > 0)
 436         {
 437             const labelList& agglom = finalAgglom_[patchID];
 438             label nAgglom = max(agglom) + 1;
 439
 440             labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
 441
 442             forAll(coarseToFine, coarseI)
 443             {
 444                 const label coarseFaceID = coarsePatchFace[coarseI];
 445                 const labelList& fineFaces = coarseToFine[coarseFaceID];
 446                 UIndirectList<scalar> fineSf
 447                 (
 448                     sf,
 449                     fineFaces
 450                 );
 451                 scalar area = sum(fineSf());
 452                 // Temperature, emissivity and external flux area weighting
 453                 forAll(fineFaces, j)
 454                 {
 455                     label faceI = fineFaces[j];
 456                     Tave[coarseI] += (Tp[faceI]*sf[faceI])/area;
 457                     Eave[coarseI] += (eb[faceI]*sf[faceI])/area;
 458                     Hoiave[coarseI] += (Hoi[faceI]*sf[faceI])/area;
 459                 }
 460             }
 461         }
 462
 463         localCoarseTave.append(Tave);
 464         localCoarseEave.append(Eave);
 465         localCoarseHoave.append(Hoiave);
 466     }
 467
 468     // Fill the local values to distribute
 469     SubList<scalar>(compactCoarseT,nLocalCoarseFaces_).assign(localCoarseTave);
 470     SubList<scalar>(compactCoarseE,nLocalCoarseFaces_).assign(localCoarseEave);
 471     SubList<scalar>
 472         (compactCoarseHo,nLocalCoarseFaces_).assign(localCoarseHoave);
 473
 474     // Distribute data
 475     map_->distribute(compactCoarseT);
 476     map_->distribute(compactCoarseE);
 477     map_->distribute(compactCoarseHo);
 478
 479     // Distribute local global ID
 480     labelList compactGlobalIds(map_->constructSize(), 0.0);
 481
 482     labelList localGlobalIds(nLocalCoarseFaces_);
 483
 484     for(label k = 0; k < nLocalCoarseFaces_; k++)
 485     {
 486         localGlobalIds[k] = globalNumbering.toGlobal(Pstream::myProcNo(), k);
 487     }
 488
 489     SubList<label>
 490     (
 491         compactGlobalIds,
 492         nLocalCoarseFaces_
 493     ).assign(localGlobalIds);
 494
 495     map_->distribute(compactGlobalIds);
 496
 497     // Create global size vectors
 498     scalarField T(totalNCoarseFaces_, 0.0);
 499     scalarField E(totalNCoarseFaces_, 0.0);
 500     scalarField QrExt(totalNCoarseFaces_, 0.0);
 501
 502     // Fill lists from compact to global indexes.
 503     forAll(compactCoarseT, i)
 504     {
 505         T[compactGlobalIds[i]] = compactCoarseT[i];
 506         E[compactGlobalIds[i]] = compactCoarseE[i];
 507         QrExt[compactGlobalIds[i]] = compactCoarseHo[i];
 508     }
 509
 510     Pstream::listCombineGather(T, maxEqOp<scalar>());
 511     Pstream::listCombineGather(E, maxEqOp<scalar>());
 512     Pstream::listCombineGather(QrExt, maxEqOp<scalar>());
 513
 514     Pstream::listCombineScatter(T);
 515     Pstream::listCombineScatter(E);
 516     Pstream::listCombineScatter(QrExt);
 517
 518     // Net radiation
 519     scalarField q(totalNCoarseFaces_, 0.0);
 520
 521     if (Pstream::master())
 522     {
 523         // Variable emissivity
 524         if (!constEmissivity_)
 525         {
 526             scalarSquareMatrix C(totalNCoarseFaces_, 0.0);
 527
 528             for (label i=0; i<totalNCoarseFaces_; i++)
 529             {
 530                 for (label j=0; j<totalNCoarseFaces_; j++)
 531                 {
 532                     scalar invEj = 1.0/E[j];
 533                     scalar sigmaT4 =
 534                         radiation::sigmaSB.value()*pow(T[j], 4.0);
 535
 536                     if (i==j)
 537                     {
 538                         C[i][j] = invEj - (invEj - 1.0)*Fmatrix_()[i][j];
 539                         q[i] += (Fmatrix_()[i][j] - 1.0)*sigmaT4 - QrExt[j];
 540                     }
 541                     else
 542                     {
 543                         C[i][j] = (1.0 - invEj)*Fmatrix_()[i][j];
 544                         q[i] += Fmatrix_()[i][j]*sigmaT4 - QrExt[j];
 545                     }
 546
 547                 }
 548             }
 549
 550             Info<< "\nSolving view factor equations..." << endl;
 551             // Negative coming into the fluid
 552             scalarSquareMatrix::LUsolve(C, q);
 553         }
 554         else //Constant emissivity
 555         {
 556             // Initial iter calculates CLU and chaches it
 557             if (iterCounter_ == 0)
 558             {
 559                 for (label i=0; i<totalNCoarseFaces_; i++)
 560                 {
 561                     for (label j=0; j<totalNCoarseFaces_; j++)
 562                     {
 563                         scalar invEj = 1.0/E[j];
 564                         if (i==j)
 565                         {
 566                             CLU_()[i][j] = invEj-(invEj-1.0)*Fmatrix_()[i][j];
 567                         }
 568                         else
 569                         {
 570                             CLU_()[i][j] = (1.0 - invEj)*Fmatrix_()[i][j];
 571                         }
 572                     }
 573                 }
 574                 Info<< "\nDecomposing C matrix..." << endl;
 575                 scalarSquareMatrix::LUDecompose(CLU_(), pivotIndices_);
 576             }
 577
 578             for (label i=0; i<totalNCoarseFaces_; i++)
 579             {
 580                 for (label j=0; j<totalNCoarseFaces_; j++)
 581                 {
 582                     scalar sigmaT4 =
 583                         radiation::sigmaSB.value()
 584                        *pow(T[j], 4.0);
 585
 586                     if (i==j)
 587                     {
 588                         q[i] += (Fmatrix_()[i][j] - 1.0)*sigmaT4 - QrExt[j];
 589                     }
 590                     else
 591                     {
 592                         q[i] += Fmatrix_()[i][j]*sigmaT4 - QrExt[j];
 593                     }
 594                 }
 595             }
 596
 597             Info<< "\nLU Back substitute C matrix.." << endl;
 598             scalarSquareMatrix::LUBacksubstitute(CLU_(), pivotIndices_, q);
 599             iterCounter_ ++;
 600         }
 601     }
 602
 603     // Scatter q and fill Qr
 604     Pstream::listCombineScatter(q);
 605     Pstream::listCombineGather(q, maxEqOp<scalar>());
 606
 607
 608     label globCoarseId = 0;
 609     forAll(selectedPatches_, i)
 610     {
 611         const label patchID = selectedPatches_[i];
 612         const polyPatch& pp = mesh_.boundaryMesh()[patchID];
 613         if (pp.size() > 0)
 614         {
 615             scalarField& Qrp = Qr_.boundaryField()[patchID];
 616             const scalarField& sf = mesh_.magSf().boundaryField()[patchID];
 617             const labelList& agglom = finalAgglom_[patchID];
 618             label nAgglom = max(agglom)+1;
 619
 620             labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
 621
 622             const labelList& coarsePatchFace =
 623                 coarseMesh_.patchFaceMap()[patchID];
 624
 625             scalar heatFlux = 0.0;
 626             forAll(coarseToFine, coarseI)
 627             {
 628                 label globalCoarse =
 629                     globalNumbering.toGlobal(Pstream::myProcNo(), globCoarseId);
 630                 const label coarseFaceID = coarsePatchFace[coarseI];
 631                 const labelList& fineFaces = coarseToFine[coarseFaceID];
 632                 forAll(fineFaces, k)
 633                 {
 634                     label faceI = fineFaces[k];
 635
 636                     Qrp[faceI] = q[globalCoarse];
 637                     heatFlux += Qrp[faceI]*sf[faceI];
 638                 }
 639                 globCoarseId ++;
 640             }
 641         }
 642     }
 643
 644     if (debug)
 645     {
 646         forAll(Qr_.boundaryField(), patchID)
 647         {
 648             const scalarField& Qrp = Qr_.boundaryField()[patchID];
 649             const scalarField& magSf = mesh_.magSf().boundaryField()[patchID];
 650             scalar heatFlux = gSum(Qrp*magSf);
 651             Info<< "Total heat flux at patch: "
 652                 << patchID << " "
 653                 << heatFlux << " [W]" << endl;
 654         }
 655     }
 656
 657     // Relax Qr if necessary
 658     Qr_.relax();
 659 }
 660
 661
 662 Foam::tmp<Foam::volScalarField> Foam::radiation::viewFactor::Rp() const
 663 {
 664     return tmp<volScalarField>
 665     (
 666         new volScalarField
 667         (
 668             IOobject
 669             (
 670                 "Rp",
 671                 mesh_.time().timeName(),
 672                 mesh_,
 673                 IOobject::NO_READ,
 674                 IOobject::NO_WRITE,
 675                 false
 676             ),
 677             mesh_,
 678             dimensionedScalar
 679             (
 680                 "zero",
 681                 dimMass/pow3(dimTime)/dimLength/pow4(dimTemperature),
 682                 0.0
 683             )
 684         )
 685     );
 686 }
 687
 688
 689 Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
 690 Foam::radiation::viewFactor::Ru() const
 691 {
 692     return tmp<DimensionedField<scalar, volMesh> >
 693     (
 694         new DimensionedField<scalar, volMesh>
 695         (
 696             IOobject
 697             (
 698                 "Ru",
 699                 mesh_.time().timeName(),
 700                 mesh_,
 701                 IOobject::NO_READ,
 702                 IOobject::NO_WRITE,
 703                 false
 704             ),
 705             mesh_,
 706             dimensionedScalar("zero", dimMass/dimLength/pow3(dimTime), 0.0)
 707         )
 708     );
 709 }
 710
 711 // ************************************************************************* //