From 549e72dff2cd01897b9a0db0cb4dfee432a2a8ff Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 26 May 2009 21:43:00 +0100 Subject: [PATCH 01/35] decomposing meshes in time directories --- .../parallelProcessing/decomposePar/decomposePar.C | 8 ++++---- .../decomposePar/domainDecomposition.C | 10 +++++----- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C index 65dd679ecd..a26179507a 100644 --- a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C +++ b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C @@ -566,7 +566,7 @@ int main(int argc, char *argv[]) IOobject ( "cellProcAddressing", - "constant", + procMesh.facesInstance(), procMesh.meshSubDir, procMesh, IOobject::MUST_READ, @@ -579,7 +579,7 @@ int main(int argc, char *argv[]) IOobject ( "boundaryProcAddressing", - "constant", + procMesh.facesInstance(), procMesh.meshSubDir, procMesh, IOobject::MUST_READ, @@ -603,7 +603,7 @@ int main(int argc, char *argv[]) IOobject ( "faceProcAddressing", - "constant", + procMesh.facesInstance(), procMesh.meshSubDir, procMesh, IOobject::MUST_READ, @@ -645,7 +645,7 @@ int main(int argc, char *argv[]) IOobject ( "pointProcAddressing", - "constant", + procMesh.facesInstance(), procMesh.meshSubDir, procMesh, IOobject::MUST_READ, diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C index 2940577d16..e310cc758f 100644 --- a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C +++ b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C @@ -269,7 +269,7 @@ bool domainDecomposition::writeDecomposition() IOobject ( this->polyMesh::name(), // region name of undecomposed mesh - "constant", + pointsInstance(), processorDb ), xferMove(procPoints), @@ -620,7 +620,7 @@ bool domainDecomposition::writeDecomposition() IOobject ( "pointProcAddressing", - "constant", + procMesh.facesInstance(), procMesh.meshSubDir, procMesh, IOobject::NO_READ, @@ -635,7 +635,7 @@ bool domainDecomposition::writeDecomposition() IOobject ( "faceProcAddressing", - "constant", + procMesh.facesInstance(), procMesh.meshSubDir, procMesh, IOobject::NO_READ, @@ -650,7 +650,7 @@ bool domainDecomposition::writeDecomposition() IOobject ( "cellProcAddressing", - "constant", + procMesh.facesInstance(), procMesh.meshSubDir, procMesh, IOobject::NO_READ, @@ -665,7 +665,7 @@ bool domainDecomposition::writeDecomposition() IOobject ( "boundaryProcAddressing", - "constant", + procMesh.facesInstance(), procMesh.meshSubDir, procMesh, IOobject::NO_READ, From 0e75278db97b89b02fec68ebbc168b88bc75ed5f Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 27 May 2009 18:29:49 +0100 Subject: [PATCH 02/35] bugfix: << operator --- .../ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C index bbf58e49a0..b836c1c412 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C +++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C @@ -276,15 +276,15 @@ Foam::Ostream& Foam::operator<< scalarField YSolidLoc = p.YSolid()*p.Y()[2]; if (os.format() == IOstream::ASCII) { - os << static_cast&>(p) + os << static_cast&>(p) << token::SPACE << YGasLoc << token::SPACE << YLiquidLoc << token::SPACE << YSolidLoc; } else { - os << static_cast&>(p); - os << YGasLoc << YLiquidLoc << YSolidLoc; + os << static_cast&>(p); + os << YGasLoc << YLiquidLoc << YSolidLoc; } // Check state of Ostream From 66cd36309275430aaf5eda9907d82fc98b4c3638 Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 27 May 2009 19:22:35 +0100 Subject: [PATCH 03/35] updated --- .../mutRoughWallFunctionFvPatchScalarField.C | 28 +++++++++---------- .../mutRoughWallFunctionFvPatchScalarField.H | 7 +---- 2 files changed, 15 insertions(+), 20 deletions(-) diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutRoughWallFunction/mutRoughWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutRoughWallFunction/mutRoughWallFunctionFvPatchScalarField.C index 38a275aeb5..a6b699d20d 100644 --- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutRoughWallFunction/mutRoughWallFunctionFvPatchScalarField.C +++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutRoughWallFunction/mutRoughWallFunctionFvPatchScalarField.C @@ -45,25 +45,23 @@ namespace RASModels scalar mutRoughWallFunctionFvPatchScalarField::fnRough ( const scalar KsPlus, - const scalar Cs, - const scalar kappa + const scalar Cs ) const { - // Set deltaB based on non-dimensional roughness height - scalar deltaB = 0.0; + // Return fn based on non-dimensional roughness height + if (KsPlus < 90.0) { - deltaB = - 1.0/kappa - *log((KsPlus - 2.25)/87.75 + Cs*KsPlus) - *sin(0.4258*(log(KsPlus) - 0.811)); + return pow + ( + (KsPlus - 2.25)/87.75 + Cs*KsPlus, + sin(0.4258*(log(KsPlus) - 0.811)) + ); } else { - deltaB = 1.0/kappa*log(1.0 + Cs*KsPlus); + return (1.0 + Cs*KsPlus); } - - return exp(min(deltaB*kappa, 50.0)); } @@ -216,8 +214,8 @@ void mutRoughWallFunctionFvPatchScalarField::updateCoeffs() scalar yPlusLamNew = yPlusLam; if (KsPlus > 2.25) { - Edash = E/fnRough(KsPlus, Cs_[faceI], kappa); - yPlusLam = rasModel.yPlusLam(kappa, Edash); + Edash = E/fnRough(KsPlus, Cs_[faceI]); + yPlusLamNew = rasModel.yPlusLam(kappa, Edash); } if (debug) @@ -231,7 +229,9 @@ void mutRoughWallFunctionFvPatchScalarField::updateCoeffs() if (yPlus > yPlusLamNew) { - mutw[faceI] = muw[faceI]*(yPlus*kappa/log(Edash*yPlus) - 1); + mutw[faceI] = + muw[faceI] + *(yPlus*kappa/log(max(Edash*yPlus, 1+1e-4)) - 1); } else { diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutRoughWallFunction/mutRoughWallFunctionFvPatchScalarField.H b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutRoughWallFunction/mutRoughWallFunctionFvPatchScalarField.H index 88209f4184..b12ad6bc57 100644 --- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutRoughWallFunction/mutRoughWallFunctionFvPatchScalarField.H +++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutRoughWallFunction/mutRoughWallFunctionFvPatchScalarField.H @@ -83,12 +83,7 @@ class mutRoughWallFunctionFvPatchScalarField // Private member functions //- Compute the roughness function - scalar fnRough - ( - const scalar KsPlus, - const scalar Cs, - const scalar kappa - ) const; + scalar fnRough(const scalar KsPlus, const scalar Cs) const; public: From fd9cf0c1f6690023043b2a10036359bd21cf0192 Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 27 May 2009 19:23:05 +0100 Subject: [PATCH 04/35] removed kappa from argument list - not used --- .../nutRoughWallFunctionFvPatchScalarField.C | 5 ++--- .../nutRoughWallFunctionFvPatchScalarField.H | 7 +------ 2 files changed, 3 insertions(+), 9 deletions(-) diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.C index 7b6e14493e..afb08a5d52 100644 --- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.C +++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.C @@ -45,8 +45,7 @@ namespace RASModels scalar nutRoughWallFunctionFvPatchScalarField::fnRough ( const scalar KsPlus, - const scalar Cs, - const scalar kappa + const scalar Cs ) const { // Return fn based on non-dimensional roughness height @@ -205,7 +204,7 @@ void nutRoughWallFunctionFvPatchScalarField::updateCoeffs() if (KsPlus > 2.25) { - Edash = E/fnRough(KsPlus, Cs_[faceI], kappa); + Edash = E/fnRough(KsPlus, Cs_[faceI]); } if (yPlus > yPlusLam) diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.H index 968c5c3f8f..3371cbaa6d 100644 --- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.H +++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.H @@ -80,12 +80,7 @@ class nutRoughWallFunctionFvPatchScalarField // Private member functions //- Compute the roughness function - scalar fnRough - ( - const scalar KsPlus, - const scalar Cs, - const scalar kappa - ) const; + scalar fnRough(const scalar KsPlus, const scalar Cs) const; public: From 0a976928574d9cc355babbde63ce6cd5de6bf0b4 Mon Sep 17 00:00:00 2001 From: andy Date: Thu, 28 May 2009 11:46:54 +0100 Subject: [PATCH 05/35] bugfix: changed order of logic to prevent hand on redude in parallel --- src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C index 5ad5dca8cd..5eae3a9f80 100644 --- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C +++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C @@ -480,7 +480,7 @@ void Foam::fvMatrix::setReference const bool forceReference ) { - if (celli >= 0 && (psi_.needReference() || forceReference)) + if ((forceReference || psi_.needReference()) && celli >= 0) { source()[celli] += diag()[celli]*value; diag()[celli] += diag()[celli]; From 55f5241fc826472a0b7acda87a225b63a1cb3a04 Mon Sep 17 00:00:00 2001 From: mattijs Date: Thu, 28 May 2009 18:14:38 +0100 Subject: [PATCH 06/35] output wallDistance --- .../postProcessing/wall/yPlusLES/Make/options | 1 + .../postProcessing/wall/yPlusLES/yPlusLES.C | 17 ++++++++++++++++- .../postProcessing/wall/yPlusRAS/Make/options | 1 + .../postProcessing/wall/yPlusRAS/yPlusRAS.C | 16 +++++++++++++++- 4 files changed, 33 insertions(+), 2 deletions(-) diff --git a/applications/utilities/postProcessing/wall/yPlusLES/Make/options b/applications/utilities/postProcessing/wall/yPlusLES/Make/options index 79fbbac91e..f6131ce41c 100644 --- a/applications/utilities/postProcessing/wall/yPlusLES/Make/options +++ b/applications/utilities/postProcessing/wall/yPlusLES/Make/options @@ -1,4 +1,5 @@ EXE_INC = \ + -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/turbulenceModels \ -I$(LIB_SRC)/turbulenceModels/incompressible/LES/LESModel \ -I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \ diff --git a/applications/utilities/postProcessing/wall/yPlusLES/yPlusLES.C b/applications/utilities/postProcessing/wall/yPlusLES/yPlusLES.C index 82ee27fa0a..67dc0df3f0 100644 --- a/applications/utilities/postProcessing/wall/yPlusLES/yPlusLES.C +++ b/applications/utilities/postProcessing/wall/yPlusLES/yPlusLES.C @@ -34,6 +34,7 @@ Description #include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H" #include "LESModel.H" #include "nearWallDist.H" +#include "wallDist.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -49,7 +50,18 @@ int main(int argc, char *argv[]) { runTime.setTime(timeDirs[timeI], timeI); Info<< "Time = " << runTime.timeName() << endl; - mesh.readUpdate(); + fvMesh::readUpdateState state = mesh.readUpdate(); + + // Wall distance + if (timeI == 0 || state != fvMesh::UNCHANGED) + { + Info<< "Calculating wall distance\n" << endl; + wallDist y(mesh, true); + Info<< "Writing wall distance to field " + << y.name() << nl << endl; + y.write(); + } + volScalarField yPlus ( @@ -116,6 +128,9 @@ int main(int argc, char *argv[]) } } + Info<< "Writing yPlus to field " + << yPlus.name() << nl << endl; + yPlus.write(); } diff --git a/applications/utilities/postProcessing/wall/yPlusRAS/Make/options b/applications/utilities/postProcessing/wall/yPlusRAS/Make/options index 89632547e0..9ec5fca52b 100644 --- a/applications/utilities/postProcessing/wall/yPlusRAS/Make/options +++ b/applications/utilities/postProcessing/wall/yPlusRAS/Make/options @@ -1,4 +1,5 @@ EXE_INC = \ + -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/transportModels \ -I$(LIB_SRC)/turbulenceModels \ -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \ diff --git a/applications/utilities/postProcessing/wall/yPlusRAS/yPlusRAS.C b/applications/utilities/postProcessing/wall/yPlusRAS/yPlusRAS.C index 7d1e4796f2..a93712263e 100644 --- a/applications/utilities/postProcessing/wall/yPlusRAS/yPlusRAS.C +++ b/applications/utilities/postProcessing/wall/yPlusRAS/yPlusRAS.C @@ -34,6 +34,7 @@ Description #include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H" #include "RASModel.H" #include "wallFvPatch.H" +#include "wallDist.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -49,7 +50,17 @@ int main(int argc, char *argv[]) { runTime.setTime(timeDirs[timeI], timeI); Info<< "Time = " << runTime.timeName() << endl; - mesh.readUpdate(); + fvMesh::readUpdateState state = mesh.readUpdate(); + + // Wall distance + if (timeI == 0 || state != fvMesh::UNCHANGED) + { + Info<< "Calculating wall distance\n" << endl; + wallDist y(mesh, true); + Info<< "Writing wall distance to field " + << y.name() << nl << endl; + y.write(); + } volScalarField yPlus ( @@ -106,6 +117,9 @@ int main(int argc, char *argv[]) } } + Info<< "Writing yPlus to field " + << yPlus.name() << nl << endl; + yPlus.write(); } From 6c2996cac456cd14495626fd2d521a1ad7117508 Mon Sep 17 00:00:00 2001 From: mattijs Date: Thu, 28 May 2009 18:14:53 +0100 Subject: [PATCH 07/35] correction from Hilary --- .../schemes/FitData/FitData.C | 23 ++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C index 18585bc645..7f432462b1 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C @@ -191,6 +191,13 @@ void Foam::FitData::calcFit ); } + // Additional weighting for constant and linear terms + for(label i = 0; i < B.n(); i++) + { + B[i][0] *= wts[0]; + B[i][1] *= wts[0]; + } + // Set the fit label stencilSize = C.size(); coeffsi.setSize(stencilSize); @@ -205,7 +212,7 @@ void Foam::FitData::calcFit for(label i=0; i maxCoeff) { maxCoeff = mag(coeffsi[i]); @@ -229,7 +236,8 @@ void Foam::FitData::calcFit // << " sing vals " << svd.S() << endl; // } - if (!goodFit) // (not good fit so increase weight in the centre) + if (!goodFit) // (not good fit so increase weight in the centre and weight + // for constant and linear terms) { // if (iIt == 7) // { @@ -237,7 +245,10 @@ void Foam::FitData::calcFit // ( // "FitData::calcFit" // "(const List& C, const label facei" - // ) << "Cannot fit face " << facei + // ) << "Cannot fit face " << facei << " iteration " << iIt + // << " with sum of weights " << sum(coeffsi) << nl + // << " Weights " << coeffsi << nl + // << " Linear weights " << wLin << " " << 1 - wLin << nl // << " sing vals " << svd.S() << endl; // } @@ -249,6 +260,12 @@ void Foam::FitData::calcFit B[0][j] *= 10; B[1][j] *= 10; } + + for(label i = 0; i < B.n(); i++) + { + B[i][0] *= 10; + B[i][1] *= 10; + } } } From d4d4f7282796bcc638626de7c6209d6b4bd743d0 Mon Sep 17 00:00:00 2001 From: mattijs Date: Thu, 28 May 2009 18:15:16 +0100 Subject: [PATCH 08/35] optional reading of writeGraph --- .../scotchDecomp/scotchDecomp.C | 53 ++++++++++--------- 1 file changed, 28 insertions(+), 25 deletions(-) diff --git a/src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.C b/src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.C index a33e49b828..f1e1cd2780 100644 --- a/src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.C +++ b/src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.C @@ -153,37 +153,40 @@ Foam::label Foam::scotchDecomp::decompose const dictionary& scotchCoeffs = decompositionDict_.subDict("scotchCoeffs"); - Switch writeGraph(scotchCoeffs.lookup("writeGraph")); - - if (writeGraph) + if (scotchCoeffs.found("writeGraph")) { - OFstream str(mesh_.time().path() / mesh_.name() + ".grf"); + Switch writeGraph(scotchCoeffs.lookup("writeGraph")); - Info<< "Dumping Scotch graph file to " << str.name() << endl - << "Use this in combination with gpart." << endl; - - label version = 0; - str << version << nl; - // Numer of vertices - str << xadj.size()-1 << ' ' << adjncy.size() << nl; - // Numbering starts from 0 - label baseval = 0; - // Has weights? - label hasEdgeWeights = 0; - label hasVertexWeights = 0; - label numericflag = 10*hasEdgeWeights+hasVertexWeights; - str << baseval << ' ' << numericflag << nl; - for (label cellI = 0; cellI < xadj.size()-1; cellI++) + if (writeGraph) { - label start = xadj[cellI]; - label end = xadj[cellI+1]; - str << end-start; + OFstream str(mesh_.time().path() / mesh_.name() + ".grf"); - for (label i = start; i < end; i++) + Info<< "Dumping Scotch graph file to " << str.name() << endl + << "Use this in combination with gpart." << endl; + + label version = 0; + str << version << nl; + // Numer of vertices + str << xadj.size()-1 << ' ' << adjncy.size() << nl; + // Numbering starts from 0 + label baseval = 0; + // Has weights? + label hasEdgeWeights = 0; + label hasVertexWeights = 0; + label numericflag = 10*hasEdgeWeights+hasVertexWeights; + str << baseval << ' ' << numericflag << nl; + for (label cellI = 0; cellI < xadj.size()-1; cellI++) { - str << ' ' << adjncy[i]; + label start = xadj[cellI]; + label end = xadj[cellI+1]; + str << end-start; + + for (label i = start; i < end; i++) + { + str << ' ' << adjncy[i]; + } + str << nl; } - str << nl; } } } From 1779687e3dda3cfad12b67003cdebe0a1040a44d Mon Sep 17 00:00:00 2001 From: mattijs Date: Thu, 28 May 2009 18:22:38 +0100 Subject: [PATCH 09/35] duplicate function removed --- .../fvMesh/extendedStencil/extendedStencil.C | 162 ------------------ .../fvMesh/extendedStencil/extendedStencil.H | 7 - 2 files changed, 169 deletions(-) diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C b/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C index 1ff00db7ba..38804ccbc1 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C +++ b/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C @@ -31,168 +31,6 @@ License // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -// Calculates per face a list of global cell/face indices. -void Foam::extendedStencil::calcFaceStencil -( - const labelListList& globalCellCells, - labelListList& faceStencil -) -{ - const polyBoundaryMesh& patches = mesh_.boundaryMesh(); - const label nBnd = mesh_.nFaces()-mesh_.nInternalFaces(); - const labelList& own = mesh_.faceOwner(); - const labelList& nei = mesh_.faceNeighbour(); - - - // Determine neighbouring global cell Cells - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - labelListList neiGlobalCellCells(nBnd); - forAll(patches, patchI) - { - const polyPatch& pp = patches[patchI]; - - if (pp.coupled()) - { - label faceI = pp.start(); - - forAll(pp, i) - { - neiGlobalCellCells[faceI-mesh_.nInternalFaces()] = - globalCellCells[own[faceI]]; - faceI++; - } - } - } - syncTools::swapBoundaryFaceList(mesh_, neiGlobalCellCells, false); - - - - // Construct stencil in global numbering - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - faceStencil.setSize(mesh_.nFaces()); - - labelHashSet faceStencilSet; - - for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++) - { - faceStencilSet.clear(); - - const labelList& ownCCells = globalCellCells[own[faceI]]; - label globalOwn = ownCCells[0]; - // Insert cellCells - forAll(ownCCells, i) - { - faceStencilSet.insert(ownCCells[i]); - } - - const labelList& neiCCells = globalCellCells[nei[faceI]]; - label globalNei = neiCCells[0]; - // Insert cellCells - forAll(neiCCells, i) - { - faceStencilSet.insert(neiCCells[i]); - } - - // Guarantee owner first, neighbour second. - faceStencil[faceI].setSize(faceStencilSet.size()); - label n = 0; - faceStencil[faceI][n++] = globalOwn; - faceStencil[faceI][n++] = globalNei; - forAllConstIter(labelHashSet, faceStencilSet, iter) - { - if (iter.key() != globalOwn && iter.key() != globalNei) - { - faceStencil[faceI][n++] = iter.key(); - } - } - //Pout<< "internalface:" << faceI << " toc:" << faceStencilSet.toc() - // << " faceStencil:" << faceStencil[faceI] << endl; - } - forAll(patches, patchI) - { - const polyPatch& pp = patches[patchI]; - label faceI = pp.start(); - - if (pp.coupled()) - { - forAll(pp, i) - { - faceStencilSet.clear(); - - const labelList& ownCCells = globalCellCells[own[faceI]]; - label globalOwn = ownCCells[0]; - forAll(ownCCells, i) - { - faceStencilSet.insert(ownCCells[i]); - } - - // And the neighbours of the coupled cell - const labelList& neiCCells = - neiGlobalCellCells[faceI-mesh_.nInternalFaces()]; - label globalNei = neiCCells[0]; - forAll(neiCCells, i) - { - faceStencilSet.insert(neiCCells[i]); - } - - // Guarantee owner first, neighbour second. - faceStencil[faceI].setSize(faceStencilSet.size()); - label n = 0; - faceStencil[faceI][n++] = globalOwn; - faceStencil[faceI][n++] = globalNei; - forAllConstIter(labelHashSet, faceStencilSet, iter) - { - if (iter.key() != globalOwn && iter.key() != globalNei) - { - faceStencil[faceI][n++] = iter.key(); - } - } - - //Pout<< "coupledface:" << faceI - // << " toc:" << faceStencilSet.toc() - // << " faceStencil:" << faceStencil[faceI] << endl; - - faceI++; - } - } - else if (!isA(pp)) - { - forAll(pp, i) - { - faceStencilSet.clear(); - - const labelList& ownCCells = globalCellCells[own[faceI]]; - label globalOwn = ownCCells[0]; - forAll(ownCCells, i) - { - faceStencilSet.insert(ownCCells[i]); - } - - // Guarantee owner first, neighbour second. - faceStencil[faceI].setSize(faceStencilSet.size()); - label n = 0; - faceStencil[faceI][n++] = globalOwn; - forAllConstIter(labelHashSet, faceStencilSet, iter) - { - if (iter.key() != globalOwn) - { - faceStencil[faceI][n++] = iter.key(); - } - } - - //Pout<< "boundaryface:" << faceI - // << " toc:" << faceStencilSet.toc() - // << " faceStencil:" << faceStencil[faceI] << endl; - - faceI++; - } - } - } -} - - Foam::autoPtr Foam::extendedStencil::calcDistributeMap ( const globalIndex& globalNumbering, diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H b/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H index 6faff8f3a6..a1c6d41474 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H +++ b/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H @@ -75,13 +75,6 @@ protected: // Protected Member Functions - //- Collect cell neighbours into extended stencil - void calcFaceStencil - ( - const labelListList& globalCellCells, - labelListList& faceStencil - ); - //- Calculate distribute map autoPtr calcDistributeMap ( From 1d14951bf57873b40fa121087800ee612cf6106c Mon Sep 17 00:00:00 2001 From: henry Date: Fri, 29 May 2009 09:32:32 +0100 Subject: [PATCH 10/35] Added turbulent diffusivity. --- applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H index c126025368..c27a3e752e 100644 --- a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H +++ b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H @@ -3,7 +3,7 @@ ( fvm::ddt(alpha1) + fvm::div(phi, alpha1) - - fvm::laplacian(Dab, alpha1) + - fvm::laplacian(turbulence->nut() + Dab, alpha1) ); alpha1Eqn.solve(); From 47beddc1c15e4d9173a611cda1fa416645fc52b5 Mon Sep 17 00:00:00 2001 From: henry Date: Fri, 29 May 2009 09:33:47 +0100 Subject: [PATCH 11/35] Added output of the patch area magnitude and the integral over the patch area magnitude. --- .../patch/patchIntegrate/patchIntegrate.C | 37 ++++++++++++++----- 1 file changed, 27 insertions(+), 10 deletions(-) diff --git a/applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C b/applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C index dc92aade57..54acae0025 100644 --- a/applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C +++ b/applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C @@ -78,7 +78,7 @@ int main(int argc, char *argv[]) // Give patch area if (isType(mesh.boundaryMesh()[patchi])) { - Info<< " Cyclic patch area: " << nl; + Info<< " Cyclic patch vector area: " << nl; label nFaces = mesh.boundaryMesh()[patchi].size(); vector sum1 = vector::zero; vector sum2 = vector::zero; @@ -92,12 +92,18 @@ int main(int argc, char *argv[]) Info<< " - half 1 = " << sum1 << ", " << mag(sum1) << nl << " - half 2 = " << sum2 << ", " << mag(sum2) << nl << " - total = " << (sum1 + sum2) << ", " - << mag(sum1 + sum2) << endl;; + << mag(sum1 + sum2) << endl; + Info<< " Cyclic patch area magnitude = " + << gSum(mesh.magSf().boundaryField()[patchi])/2.0 << endl; } else { - Info<< " Patch area = " + Info<< " Area vector of patch " + << patchName << '[' << patchi << ']' << " = " << gSum(mesh.Sf().boundaryField()[patchi]) << endl; + Info<< " Area magnitude of patch " + << patchName << '[' << patchi << ']' << " = " + << gSum(mesh.magSf().boundaryField()[patchi]) << endl; } // Read field and calc integral @@ -107,15 +113,26 @@ int main(int argc, char *argv[]) << fieldName << endl; volScalarField field(fieldHeader, mesh); - vector sumField = gSum - ( - mesh.Sf().boundaryField()[patchi] - *field.boundaryField()[patchi] - ); - Info<< " Integral of " << fieldName << " over patch " + Info<< " Integral of " << fieldName + << " over vector area of patch " << patchName << '[' << patchi << ']' << " = " - << sumField << nl; + << gSum + ( + mesh.Sf().boundaryField()[patchi] + *field.boundaryField()[patchi] + ) + << nl; + + Info<< " Integral of " << fieldName + << " over area magnitude of patch " + << patchName << '[' << patchi << ']' << " = " + << gSum + ( + mesh.magSf().boundaryField()[patchi] + *field.boundaryField()[patchi] + ) + << nl; } else if ( From 20dcad4e57cd6c6b071c0cc66bffe2cb73beeccb Mon Sep 17 00:00:00 2001 From: henry Date: Fri, 29 May 2009 11:03:20 +0100 Subject: [PATCH 12/35] Added the reciprocal of the turbulent Schmidt number. --- .../solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H | 6 +++++- .../solvers/multiphase/twoLiquidMixingFoam/createFields.H | 3 +++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H index c27a3e752e..d152baba0b 100644 --- a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H +++ b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H @@ -3,7 +3,11 @@ ( fvm::ddt(alpha1) + fvm::div(phi, alpha1) - - fvm::laplacian(turbulence->nut() + Dab, alpha1) + - fvm::laplacian + ( + Dab + alphatab*turbulence->nut(), alpha1, + "laplacian(Dab,alpha1)" + ) ); alpha1Eqn.solve(); diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H index a38135dc19..3b5e064372 100644 --- a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H +++ b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H @@ -50,6 +50,9 @@ dimensionedScalar Dab(twoPhaseProperties.lookup("Dab")); + // Read the reciprocal of the turbulent Schmidt number + dimensionedScalar alphatab(twoPhaseProperties.lookup("alphatab")); + // Need to store rho for ddt(rho, U) volScalarField rho("rho", alpha1*rho1 + (scalar(1) - alpha1)*rho2); rho.oldTime(); From 80175df2ddfe9f0703726f03c3edbb675c41cf2b Mon Sep 17 00:00:00 2001 From: henry Date: Fri, 29 May 2009 11:57:15 +0100 Subject: [PATCH 13/35] Changed the BC to switch between "pd" and static pressure according to the variable name. --- .../fixedFluxBuoyantPressureFvPatchScalarField.C | 12 +++++++++++- .../fixedFluxBuoyantPressureFvPatchScalarField.H | 5 ++++- 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxBuoyantPressure/fixedFluxBuoyantPressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxBuoyantPressure/fixedFluxBuoyantPressureFvPatchScalarField.C index c31f5acc5c..b6f48ff6ad 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxBuoyantPressure/fixedFluxBuoyantPressureFvPatchScalarField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxBuoyantPressure/fixedFluxBuoyantPressureFvPatchScalarField.C @@ -113,7 +113,17 @@ void fixedFluxBuoyantPressureFvPatchScalarField::updateCoeffs() const fvPatchField& rho = patch().lookupPatchField("rho"); - gradient() = -rho.snGrad()*(g.value() & patch().Cf()); + // If the variable name is "pd" assume it is p - rho*g.h + // and set the gradient appropriately. + // Otherwise assume the variable is the static pressure. + if (dimensionedInternalField().name() == "pd") + { + gradient() = -rho.snGrad()*(g.value() & patch().Cf()); + } + else + { + gradient() = rho*(g.value() & patch().nf()); + } fixedGradientFvPatchScalarField::updateCoeffs(); } diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxBuoyantPressure/fixedFluxBuoyantPressureFvPatchScalarField.H b/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxBuoyantPressure/fixedFluxBuoyantPressureFvPatchScalarField.H index 3ce584bf87..8e433cb8f1 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxBuoyantPressure/fixedFluxBuoyantPressureFvPatchScalarField.H +++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxBuoyantPressure/fixedFluxBuoyantPressureFvPatchScalarField.H @@ -26,7 +26,10 @@ Class Foam::fixedFluxBuoyantPressureFvPatchScalarField Description - Foam::fixedFluxBuoyantPressureFvPatchScalarField + Set the pressure gradient boundary condition appropriately for buoyant flow. + + If the variable name is "pd" assume it is p - rho*g.h and set the gradient + appropriately. Otherwise assume the variable is the static pressure. SourceFiles fixedFluxBuoyantPressureFvPatchScalarField.C From feb846bd435bdd9a8e2b1df328a7a3226bee5f43 Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 29 May 2009 15:27:55 +0100 Subject: [PATCH 14/35] work around scotch fpe bug --- .../decomposePar/decomposeParDict | 12 ++- .../scotchDecomp/scotchDecomp.C | 94 +++++++++++++++++-- 2 files changed, 95 insertions(+), 11 deletions(-) diff --git a/applications/utilities/parallelProcessing/decomposePar/decomposeParDict b/applications/utilities/parallelProcessing/decomposePar/decomposeParDict index 0c566d5081..aad15ee459 100644 --- a/applications/utilities/parallelProcessing/decomposePar/decomposeParDict +++ b/applications/utilities/parallelProcessing/decomposePar/decomposeParDict @@ -55,7 +55,17 @@ metisCoeffs } scotchCoeffs -{} +{ + //processorWeights + //( + // 1 + // 1 + // 1 + // 1 + //); + //writeGraph true; + //strategy "b"; +} manualCoeffs { diff --git a/src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.C b/src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.C index f1e1cd2780..994a824192 100644 --- a/src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.C +++ b/src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.C @@ -75,6 +75,24 @@ extern "C" } +// Hack: scotch generates floating point errors so need to switch of error +// trapping! +#if defined(linux) || defined(linuxAMD64) || defined(linuxIA64) +# define LINUX +#endif + +#if defined(LINUX) && defined(__GNUC__) +# define LINUX_GNUC +#endif + +#ifdef LINUX_GNUC +# ifndef __USE_GNU +# define __USE_GNU +# endif +# include +#endif + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam @@ -113,13 +131,30 @@ Foam::label Foam::scotchDecomp::decompose { // Strategy // ~~~~~~~~ + // Default. SCOTCH_Strat stradat; check(SCOTCH_stratInit(&stradat), "SCOTCH_stratInit"); - //SCOTCH_stratGraphMap(&stradat, &argv[i][2]); - //fprintf(stdout, "S\tStrat="); - //SCOTCH_stratSave(&stradat, stdout); - //fprintf(stdout, "\n"); + + if (decompositionDict_.found("scotchCoeffs")) + { + const dictionary& scotchCoeffs = + decompositionDict_.subDict("scotchCoeffs"); + + + string strategy; + if (scotchCoeffs.readIfPresent("strategy", strategy)) + { + if (debug) + { + Info<< "scotchDecomp : Using strategy " << strategy << endl; + } + SCOTCH_stratGraphMap(&stradat, strategy.c_str()); + //fprintf(stdout, "S\tStrat="); + //SCOTCH_stratSave(&stradat, stdout); + //fprintf(stdout, "\n"); + } + } // Graph @@ -198,12 +233,36 @@ Foam::label Foam::scotchDecomp::decompose SCOTCH_Arch archdat; check(SCOTCH_archInit(&archdat), "SCOTCH_archInit"); - check - ( - // SCOTCH_archCmpltw for weighted. - SCOTCH_archCmplt(&archdat, nProcessors_), - "SCOTCH_archCmplt" - ); + + List