applications/utilities/errorEstimation/estimateScalarError/estimateScalarError.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 Application
  25     estimateScalarError
  26
  27 Description
  28     Estimates the error in the solution for a scalar transport equation in the
  29     standard form
  30
  31 \*---------------------------------------------------------------------------*/
  32
  33 #include "fvCFD.H"
  34 #include "errorEstimate.H"
  35 #include "resError.H"
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38
  39 int main(int argc, char *argv[])
  40 {
  41     timeSelector::addOptions();
  42
  43 #   include "setRootCase.H"
  44 #   include "createTime.H"
  45
  46     instantList timeDirs = timeSelector::select0(runTime, args);
  47
  48 #   include "createMesh.H"
  49
  50     Info<< "\nEstimating error in scalar transport equation\n"
  51         << "Reading transportProperties\n" << endl;
  52
  53     IOdictionary transportProperties
  54     (
  55         IOobject
  56         (
  57             "transportProperties",
  58             runTime.constant(),
  59             mesh,
  60             IOobject::MUST_READ,
  61             IOobject::NO_WRITE
  62         )
  63     );
  64
  65
  66     Info<< "Reading diffusivity DT\n" << endl;
  67
  68     dimensionedScalar DT
  69     (
  70         transportProperties.lookup("DT")
  71     );
  72
  73
  74     forAll(timeDirs, timeI)
  75     {
  76         runTime.setTime(timeDirs[timeI], timeI);
  77
  78         Info<< "Time = " << runTime.timeName() << endl;
  79
  80         mesh.readUpdate();
  81
  82         IOobject THeader
  83         (
  84             "T",
  85             runTime.timeName(),
  86             mesh,
  87             IOobject::MUST_READ
  88         );
  89
  90         IOobject Uheader
  91         (
  92             "U",
  93             runTime.timeName(),
  94             mesh,
  95             IOobject::MUST_READ
  96         );
  97
  98         if (THeader.headerOk() && Uheader.headerOk())
  99         {
 100             Info << "Reading T" << endl;
 101             volScalarField T(THeader, mesh);
 102
 103             Info << "Reading U" << endl;
 104             volVectorField U(Uheader, mesh);
 105
 106 #           include "createPhi.H"
 107
 108             errorEstimate<scalar> ee
 109             (
 110                 resError::div(phi, T)
 111               - resError::laplacian(DT, T)
 112             );
 113
 114             ee.residual()().write();
 115             volScalarField e = ee.error();
 116             e.write();
 117             mag(e)().write();
 118         }
 119         else
 120         {
 121             Info<< "    No T or U" << endl;
 122         }
 123
 124         Info<< endl;
 125     }
 126
 127     Info<< "End\n" << endl;
 128
 129     return 0;
 130 }
 131
 132
 133 // ************************************************************************* //