Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073 #include <postCalc/calc.H>
00074 #include <finiteVolume/fvc.H>
00075
00076
00077
00078 namespace Foam
00079 {
00080 tmp<volScalarField> Co(const surfaceScalarField& Cof)
00081 {
00082 const fvMesh& mesh = Cof.mesh();
00083
00084 tmp<volScalarField> tCo
00085 (
00086 new volScalarField
00087 (
00088 IOobject
00089 (
00090 "Co",
00091 mesh.time().timeName(),
00092 mesh
00093 ),
00094 mesh,
00095 dimensionedScalar("0", Cof.dimensions(), 0)
00096 )
00097 );
00098
00099 volScalarField& Co = tCo();
00100
00101
00102 const unallocLabelList& owner = mesh.owner();
00103 const unallocLabelList& neighbour = mesh.neighbour();
00104
00105 forAll(owner, facei)
00106 {
00107 label own = owner[facei];
00108 label nei = neighbour[facei];
00109
00110 Co[own] = max(Co[own], Cof[facei]);
00111 Co[nei] = max(Co[nei], Cof[facei]);
00112 }
00113
00114 forAll(Co.boundaryField(), patchi)
00115 {
00116 Co.boundaryField()[patchi] = Cof.boundaryField()[patchi];
00117 }
00118
00119 return tCo;
00120 }
00121 }
00122
00123
00124 void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
00125 {
00126 bool writeResults = !args.optionFound("noWrite");
00127
00128 IOobject phiHeader
00129 (
00130 "phi",
00131 runTime.timeName(),
00132 mesh,
00133 IOobject::MUST_READ
00134 );
00135
00136 if (phiHeader.headerOk())
00137 {
00138 autoPtr<surfaceScalarField> CoPtr;
00139
00140 Info<< " Reading phi" << endl;
00141 surfaceScalarField phi(phiHeader, mesh);
00142 Info<< " Calculating Co" << endl;
00143
00144 if (phi.dimensions() == dimensionSet(1, 0, -1, 0, 0))
00145 {
00146
00147 volScalarField rho
00148 (
00149 IOobject
00150 (
00151 "rho",
00152 runTime.timeName(),
00153 mesh,
00154 IOobject::MUST_READ
00155 ),
00156 mesh
00157 );
00158
00159 CoPtr.set
00160 (
00161 new surfaceScalarField
00162 (
00163 IOobject
00164 (
00165 "Cof",
00166 runTime.timeName(),
00167 mesh,
00168 IOobject::NO_READ
00169 ),
00170 (
00171 mesh.surfaceInterpolation::deltaCoeffs()
00172 * (mag(phi)/(fvc::interpolate(rho)*mesh.magSf()))
00173 * runTime.deltaT()
00174 )
00175 )
00176 );
00177 }
00178 else if (phi.dimensions() == dimensionSet(0, 3, -1, 0, 0))
00179 {
00180
00181 CoPtr.set
00182 (
00183 new surfaceScalarField
00184 (
00185 IOobject
00186 (
00187 "Cof",
00188 runTime.timeName(),
00189 mesh,
00190 IOobject::NO_READ
00191 ),
00192 (
00193 mesh.surfaceInterpolation::deltaCoeffs()
00194 * (mag(phi)/mesh.magSf())
00195 * runTime.deltaT()
00196 )
00197 )
00198 );
00199 }
00200 else
00201 {
00202 FatalErrorIn(args.executable())
00203 << "Incorrect dimensions of phi: " << phi.dimensions()
00204 << abort(FatalError);
00205 }
00206
00207 Info<< "Co max : " << max(CoPtr()).value() << endl;
00208
00209 if (writeResults)
00210 {
00211 CoPtr().write();
00212 Co(CoPtr())().write();
00213 }
00214 }
00215 else
00216 {
00217 Info<< " No phi" << endl;
00218 }
00219
00220 Info<< "\nEnd\n" << endl;
00221 }
00222
00223