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 #include "cellLimitedGrad.H"
00027 #include <finiteVolume/gaussGrad.H>
00028 #include <finiteVolume/fvMesh.H>
00029 #include <finiteVolume/volMesh.H>
00030 #include <finiteVolume/surfaceMesh.H>
00031 #include <finiteVolume/volFields.H>
00032 #include <finiteVolume/fixedValueFvPatchFields.H>
00033
00034
00035
00036 namespace Foam
00037 {
00038
00039
00040
00041 namespace fv
00042 {
00043
00044
00045
00046 makeFvGradScheme(cellLimitedGrad)
00047
00048
00049
00050 template<>
00051 inline void cellLimitedGrad<scalar>::limitFace
00052 (
00053 scalar& limiter,
00054 const scalar& maxDelta,
00055 const scalar& minDelta,
00056 const scalar& extrapolate
00057 )
00058 {
00059 if (extrapolate > maxDelta + VSMALL)
00060 {
00061 limiter = min(limiter, maxDelta/extrapolate);
00062 }
00063 else if (extrapolate < minDelta - VSMALL)
00064 {
00065 limiter = min(limiter, minDelta/extrapolate);
00066 }
00067 }
00068
00069 template<class Type>
00070 inline void cellLimitedGrad<Type>::limitFace
00071 (
00072 Type& limiter,
00073 const Type& maxDelta,
00074 const Type& minDelta,
00075 const Type& extrapolate
00076 )
00077 {
00078 for(direction cmpt=0; cmpt<Type::nComponents; cmpt++)
00079 {
00080 cellLimitedGrad<scalar>::limitFace
00081 (
00082 limiter.component(cmpt),
00083 maxDelta.component(cmpt),
00084 minDelta.component(cmpt),
00085 extrapolate.component(cmpt)
00086 );
00087 }
00088 }
00089
00090
00091
00092
00093 template<>
00094 tmp<volVectorField> cellLimitedGrad<scalar>::grad
00095 (
00096 const volScalarField& vsf
00097 ) const
00098 {
00099 const fvMesh& mesh = vsf.mesh();
00100
00101 tmp<volVectorField> tGrad = basicGradScheme_().grad(vsf);
00102
00103 if (k_ < SMALL)
00104 {
00105 return tGrad;
00106 }
00107
00108 volVectorField& g = tGrad();
00109
00110 const unallocLabelList& owner = mesh.owner();
00111 const unallocLabelList& neighbour = mesh.neighbour();
00112
00113 const volVectorField& C = mesh.C();
00114 const surfaceVectorField& Cf = mesh.Cf();
00115
00116 scalarField maxVsf(vsf.internalField());
00117 scalarField minVsf(vsf.internalField());
00118
00119 forAll(owner, facei)
00120 {
00121 label own = owner[facei];
00122 label nei = neighbour[facei];
00123
00124 scalar vsfOwn = vsf[own];
00125 scalar vsfNei = vsf[nei];
00126
00127 maxVsf[own] = max(maxVsf[own], vsfNei);
00128 minVsf[own] = min(minVsf[own], vsfNei);
00129
00130 maxVsf[nei] = max(maxVsf[nei], vsfOwn);
00131 minVsf[nei] = min(minVsf[nei], vsfOwn);
00132 }
00133
00134
00135 const volScalarField::GeometricBoundaryField& bsf = vsf.boundaryField();
00136
00137 forAll(bsf, patchi)
00138 {
00139 const fvPatchScalarField& psf = bsf[patchi];
00140
00141 const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
00142
00143 if (psf.coupled())
00144 {
00145 scalarField psfNei = psf.patchNeighbourField();
00146
00147 forAll(pOwner, pFacei)
00148 {
00149 label own = pOwner[pFacei];
00150 scalar vsfNei = psfNei[pFacei];
00151
00152 maxVsf[own] = max(maxVsf[own], vsfNei);
00153 minVsf[own] = min(minVsf[own], vsfNei);
00154 }
00155 }
00156 else
00157 {
00158 forAll(pOwner, pFacei)
00159 {
00160 label own = pOwner[pFacei];
00161 scalar vsfNei = psf[pFacei];
00162
00163 maxVsf[own] = max(maxVsf[own], vsfNei);
00164 minVsf[own] = min(minVsf[own], vsfNei);
00165 }
00166 }
00167 }
00168
00169 maxVsf -= vsf;
00170 minVsf -= vsf;
00171
00172 if (k_ < 1.0)
00173 {
00174 scalarField maxMinVsf = (1.0/k_ - 1.0)*(maxVsf - minVsf);
00175 maxVsf += maxMinVsf;
00176 minVsf -= maxMinVsf;
00177
00178
00179
00180 }
00181
00182
00183
00184 scalarField limiter(vsf.internalField().size(), 1.0);
00185
00186 forAll(owner, facei)
00187 {
00188 label own = owner[facei];
00189 label nei = neighbour[facei];
00190
00191
00192 limitFace
00193 (
00194 limiter[own],
00195 maxVsf[own],
00196 minVsf[own],
00197 (Cf[facei] - C[own]) & g[own]
00198 );
00199
00200
00201 limitFace
00202 (
00203 limiter[nei],
00204 maxVsf[nei],
00205 minVsf[nei],
00206 (Cf[facei] - C[nei]) & g[nei]
00207 );
00208 }
00209
00210 forAll(bsf, patchi)
00211 {
00212 const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
00213 const vectorField& pCf = Cf.boundaryField()[patchi];
00214
00215 forAll(pOwner, pFacei)
00216 {
00217 label own = pOwner[pFacei];
00218
00219 limitFace
00220 (
00221 limiter[own],
00222 maxVsf[own],
00223 minVsf[own],
00224 (pCf[pFacei] - C[own]) & g[own]
00225 );
00226 }
00227 }
00228
00229 if (fv::debug)
00230 {
00231 Info<< "gradient limiter for: " << vsf.name()
00232 << " max = " << gMax(limiter)
00233 << " min = " << gMin(limiter)
00234 << " average: " << gAverage(limiter) << endl;
00235 }
00236
00237 g.internalField() *= limiter;
00238 g.correctBoundaryConditions();
00239 gaussGrad<scalar>::correctBoundaryConditions(vsf, g);
00240
00241 return tGrad;
00242 }
00243
00244
00245 template<>
00246 tmp<volTensorField> cellLimitedGrad<vector>::grad
00247 (
00248 const volVectorField& vsf
00249 ) const
00250 {
00251 const fvMesh& mesh = vsf.mesh();
00252
00253 tmp<volTensorField> tGrad = basicGradScheme_().grad(vsf);
00254
00255 if (k_ < SMALL)
00256 {
00257 return tGrad;
00258 }
00259
00260 volTensorField& g = tGrad();
00261
00262 const unallocLabelList& owner = mesh.owner();
00263 const unallocLabelList& neighbour = mesh.neighbour();
00264
00265 const volVectorField& C = mesh.C();
00266 const surfaceVectorField& Cf = mesh.Cf();
00267
00268 vectorField maxVsf(vsf.internalField());
00269 vectorField minVsf(vsf.internalField());
00270
00271 forAll(owner, facei)
00272 {
00273 label own = owner[facei];
00274 label nei = neighbour[facei];
00275
00276 const vector& vsfOwn = vsf[own];
00277 const vector& vsfNei = vsf[nei];
00278
00279 maxVsf[own] = max(maxVsf[own], vsfNei);
00280 minVsf[own] = min(minVsf[own], vsfNei);
00281
00282 maxVsf[nei] = max(maxVsf[nei], vsfOwn);
00283 minVsf[nei] = min(minVsf[nei], vsfOwn);
00284 }
00285
00286
00287 const volVectorField::GeometricBoundaryField& bsf = vsf.boundaryField();
00288
00289 forAll(bsf, patchi)
00290 {
00291 const fvPatchVectorField& psf = bsf[patchi];
00292 const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
00293
00294 if (psf.coupled())
00295 {
00296 vectorField psfNei = psf.patchNeighbourField();
00297
00298 forAll(pOwner, pFacei)
00299 {
00300 label own = pOwner[pFacei];
00301 const vector& vsfNei = psfNei[pFacei];
00302
00303 maxVsf[own] = max(maxVsf[own], vsfNei);
00304 minVsf[own] = min(minVsf[own], vsfNei);
00305 }
00306 }
00307 else
00308 {
00309 forAll(pOwner, pFacei)
00310 {
00311 label own = pOwner[pFacei];
00312 const vector& vsfNei = psf[pFacei];
00313
00314 maxVsf[own] = max(maxVsf[own], vsfNei);
00315 minVsf[own] = min(minVsf[own], vsfNei);
00316 }
00317 }
00318 }
00319
00320 maxVsf -= vsf;
00321 minVsf -= vsf;
00322
00323 if (k_ < 1.0)
00324 {
00325 vectorField maxMinVsf = (1.0/k_ - 1.0)*(maxVsf - minVsf);
00326 maxVsf += maxMinVsf;
00327 minVsf -= maxMinVsf;
00328
00329
00330
00331 }
00332
00333
00334
00335 vectorField limiter(vsf.internalField().size(), vector::one);
00336
00337 forAll(owner, facei)
00338 {
00339 label own = owner[facei];
00340 label nei = neighbour[facei];
00341
00342
00343 limitFace
00344 (
00345 limiter[own],
00346 maxVsf[own],
00347 minVsf[own],
00348 (Cf[facei] - C[own]) & g[own]
00349 );
00350
00351
00352 limitFace
00353 (
00354 limiter[nei],
00355 maxVsf[nei],
00356 minVsf[nei],
00357 (Cf[facei] - C[nei]) & g[nei]
00358 );
00359 }
00360
00361 forAll(bsf, patchi)
00362 {
00363 const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
00364 const vectorField& pCf = Cf.boundaryField()[patchi];
00365
00366 forAll(pOwner, pFacei)
00367 {
00368 label own = pOwner[pFacei];
00369
00370 limitFace
00371 (
00372 limiter[own],
00373 maxVsf[own],
00374 minVsf[own],
00375 ((pCf[pFacei] - C[own]) & g[own])
00376 );
00377 }
00378 }
00379
00380 if (fv::debug)
00381 {
00382 Info<< "gradient limiter for: " << vsf.name()
00383 << " max = " << gMax(limiter)
00384 << " min = " << gMin(limiter)
00385 << " average: " << gAverage(limiter) << endl;
00386 }
00387
00388 tensorField& gIf = g.internalField();
00389
00390 forAll(gIf, celli)
00391 {
00392 gIf[celli] = tensor
00393 (
00394 cmptMultiply(limiter[celli], gIf[celli].x()),
00395 cmptMultiply(limiter[celli], gIf[celli].y()),
00396 cmptMultiply(limiter[celli], gIf[celli].z())
00397 );
00398 }
00399
00400 g.correctBoundaryConditions();
00401 gaussGrad<vector>::correctBoundaryConditions(vsf, g);
00402
00403 return tGrad;
00404 }
00405
00406
00407
00408
00409 }
00410
00411
00412
00413 }
00414
00415