From 3fc8c16b0d46addb865b15968a2cf2af4e5e0941 Mon Sep 17 00:00:00 2001 From: sergio Date: Tue, 22 May 2012 16:29:41 +0100 Subject: [PATCH 01/42] BUG: Fix to source in energy Eq for solid --- .../porousSolid/setPorousRegionSolidFields.H | 3 +++ .../chtMultiRegionSimpleFoam/porousSolid/solvePorousSolid.H | 2 +- .../porousSolid/setPorousRegionSolidFields.H | 3 +++ .../chtMultiRegionFoam/porousSolid/solvePorousSolid.H | 2 +- 4 files changed, 8 insertions(+), 2 deletions(-) diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/setPorousRegionSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/setPorousRegionSolidFields.H index bbbd079ef2..3ce6b03b9d 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/setPorousRegionSolidFields.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/setPorousRegionSolidFields.H @@ -15,6 +15,9 @@ const volScalarField& kappa = tkappa(); //const volSymmTensorField& K = tK(); + tmp trhoCp = cp*rho; + const volScalarField& rhoCp = trhoCp(); + volScalarField& T = thermo.T(); IObasicSourceList& sources = solidHeatSources[i]; diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/solvePorousSolid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/solvePorousSolid.H index 2dd09b974f..5d4bb47b22 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/solvePorousSolid.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/solvePorousSolid.H @@ -4,7 +4,7 @@ tmp TEqn ( - fvm::laplacian(betav*kappa, T, "laplacian(K,T)") - + sources(rho, T) + + sources(rhoCp, T) ); TEqn().relax(); diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/setPorousRegionSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/setPorousRegionSolidFields.H index 6f4834000d..2ade49d0ac 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/setPorousRegionSolidFields.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/setPorousRegionSolidFields.H @@ -15,6 +15,9 @@ const volScalarField& kappa = tkappa(); //const volSymmTensorField& K = tK(); + tmp trhoCp = cp*rho; + const volScalarField& rhoCp = trhoCp(); + volScalarField& T = thermo.T(); IObasicSourceList& sources = solidHeatSources[i]; diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/solvePorousSolid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/solvePorousSolid.H index e8941cf2ad..dfe0bff675 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/solvePorousSolid.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/solvePorousSolid.H @@ -10,7 +10,7 @@ if (finalIter) ( fvm::ddt(betav*rho*cp, T) - fvm::laplacian(betav*kappa, T, "laplacian(K,T)") - + sources(rho, T) + + sources(rhoCp, T) ); TEqn().relax(); From f22ba88d8ae2845a27edecce2f9ec1952627e519 Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 22 May 2012 17:10:26 +0100 Subject: [PATCH 02/42] ENH: extrudeMesh: allow flipNormals for extrude-from-patch --- .../extrude/extrudeMesh/extrudeMesh.C | 63 ++++++++++++++++--- .../extrude/extrudeMesh/extrudeMeshDict | 3 +- 2 files changed, 58 insertions(+), 8 deletions(-) diff --git a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C index 1c769a4346..6cd54bb461 100644 --- a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C +++ b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C @@ -251,10 +251,10 @@ int main(int argc, char *argv[]) if (mode == PATCH || mode == MESH) { - if (flipNormals) + if (flipNormals && mode == MESH) { FatalErrorIn(args.executable()) - << "Flipping normals not supported for extrusions from patch." + << "Flipping normals not supported for extrusions from mesh." << exit(FatalError); } @@ -297,6 +297,60 @@ int main(int argc, char *argv[]) addPatchCellLayer layerExtrude(mesh, (mode == MESH)); const labelList meshFaces(patchFaces(patches, sourcePatches)); + + if (mode == PATCH && flipNormals) + { + // Cheat. Flip patch faces in mesh. This invalidates the + // mesh (open cells) but does produce the correct extrusion. + polyTopoChange meshMod(mesh); + forAll(meshFaces, i) + { + label meshFaceI = meshFaces[i]; + + label patchI = patches.whichPatch(meshFaceI); + label own = mesh.faceOwner()[meshFaceI]; + label nei = -1; + if (patchI == -1) + { + nei = mesh.faceNeighbour()[meshFaceI]; + } + + label zoneI = mesh.faceZones().whichZone(meshFaceI); + bool zoneFlip = false; + if (zoneI != -1) + { + label index = mesh.faceZones()[zoneI].whichFace(meshFaceI); + zoneFlip = mesh.faceZones()[zoneI].flipMap()[index]; + } + + meshMod.modifyFace + ( + mesh.faces()[meshFaceI].reverseFace(), // modified face + meshFaceI, // label of face + own, // owner + nei, // neighbour + true, // face flip + patchI, // patch for face + zoneI, // zone for face + zoneFlip // face flip in zone + ); + } + + // Change the mesh. No inflation. + autoPtr map = meshMod.changeMesh(mesh, false); + + // Update fields + mesh.updateMesh(map); + + // Move mesh (since morphing does not do this) + if (map().hasMotionPoints()) + { + mesh.movePoints(map().preMotionPoints()); + } + } + + + indirectPrimitivePatch extrudePatch ( IndirectList @@ -471,11 +525,6 @@ int main(int argc, char *argv[]) displacement[pointI] = extrudePt - layer0Points[pointI]; } - if (flipNormals) - { - Info<< "Flipping faces." << nl << endl; - displacement = -displacement; - } // Check if wedge (has layer0 different from original patch points) // If so move the mesh to starting position. diff --git a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMeshDict b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMeshDict index f1985fe7f6..e075acbc4e 100644 --- a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMeshDict +++ b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMeshDict @@ -31,7 +31,8 @@ exposedPatchName movingWall; // If construct from surface: surface "movingWall.stl"; -// Flip surface normals before usage. +// Flip surface normals before usage. Valid only for extrude from surface or +// patch. flipNormals false; //- Linear extrusion in point-normal direction From 0136deb88e007e5034ff403f46897e36ff0f3bc5 Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 22 May 2012 17:34:03 +0100 Subject: [PATCH 03/42] ENH: extrudeMesh: make edge collape optional --- .../extrude/extrudeMesh/extrudeMesh.C | 22 ++++++++++++++++++- .../extrude/extrudeMesh/extrudeMeshDict | 4 ++++ .../system/extrudeMeshDict | 12 +++++++--- 3 files changed, 34 insertions(+), 4 deletions(-) diff --git a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C index 6cd54bb461..1f72f1d221 100644 --- a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C +++ b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C @@ -238,6 +238,25 @@ int main(int argc, char *argv[]) dict.lookup("constructFrom") ); + // Any merging of small edges + const scalar mergeTol(dict.lookupOrDefault("mergeTol", 1e-4)); + + Info<< "Extruding from " << ExtrudeModeNames[mode] + << " using model " << model().type() << endl; + if (flipNormals) + { + Info<< "Flipping normals before extruding" << endl; + } + if (mergeTol > 0) + { + Info<< "Collapsing edges < " << mergeTol << " of bounding box" << endl; + } + else + { + Info<< "Not collapsing any edges after extrusion" << endl; + } + Info<< endl; + // Generated mesh (one of either) autoPtr meshFromMesh; @@ -715,7 +734,7 @@ int main(int argc, char *argv[]) const boundBox& bb = mesh.bounds(); const vector span = bb.span(); - const scalar mergeDim = 1e-4 * bb.minDim(); + const scalar mergeDim = mergeTol * bb.minDim(); Info<< "Mesh bounding box : " << bb << nl << " with span : " << span << nl @@ -726,6 +745,7 @@ int main(int argc, char *argv[]) // Collapse edges // ~~~~~~~~~~~~~~ + if (mergeDim > 0) { Info<< "Collapsing edges < " << mergeDim << " ..." << nl << endl; diff --git a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMeshDict b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMeshDict index e075acbc4e..ecbc647797 100644 --- a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMeshDict +++ b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMeshDict @@ -88,4 +88,8 @@ sigmaRadialCoeffs // degree wedges. mergeFaces false; //true; +// Merge small edges. Fraction of bounding box. +mergeTol 0; + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_simpleFoam/system/extrudeMeshDict b/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_simpleFoam/system/extrudeMeshDict index b695ab7e5d..0bcac7b5ce 100644 --- a/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_simpleFoam/system/extrudeMeshDict +++ b/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_simpleFoam/system/extrudeMeshDict @@ -10,12 +10,14 @@ FoamFile version 2.0; format ascii; class dictionary; - object extrudeProperties; + object extrudeMeshDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // What to extrude: // patch : from patch of another case ('sourceCase') +// mesh : as above but with original case included +// surface : from externally read surface constructFrom patch; sourceCase "../wingMotion_snappyHexMesh"; @@ -24,7 +26,8 @@ sourcePatches (front); // If construct from patch: patch to use for back (can be same as sourcePatch) exposedPatchName back; -// Flip surface normals before usage. +// Flip surface normals before usage. Valid only for extrude from surface or +// patch. flipNormals false; //- Linear extrusion in point-normal direction @@ -41,7 +44,10 @@ linearNormalCoeffs // Do front and back need to be merged? Usually only makes sense for 360 // degree wedges. -mergeFaces false; +mergeFaces false; //true; + +// Merge small edges. Fraction of bounding box. +mergeTol 0; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // From 46764219625a07b70ea7f3f6c38123b3856f07ee Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 22 May 2012 17:34:37 +0100 Subject: [PATCH 04/42] ENH: foamFormatConvert: convert cvMesh's internalDelaunayVertices --- .../foamFormatConvert/foamFormatConvert.C | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C b/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C index 6f3f040a0e..8104f9e7c3 100644 --- a/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C +++ b/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -206,6 +206,14 @@ int main(int argc, char *argv[]) runTime ); + // cvMesh vertices + writeMeshObject + ( + "internalDelaunayVertices", + regionPrefix, + runTime + ); + if (runTime.writeFormat() == IOstream::ASCII) { // Only do zones when converting from binary to ascii From cb30c8389299a70575bcffcc95a0ddccda3f4437 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 23 May 2012 11:24:00 +0100 Subject: [PATCH 05/42] ENH: foamPack: enforce presence of .build --- bin/foamPack | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/bin/foamPack b/bin/foamPack index 9688cd039e..e3d090e2f2 100755 --- a/bin/foamPack +++ b/bin/foamPack @@ -3,7 +3,7 @@ # ========= | # \\ / F ield | OpenFOAM: The Open Source CFD Toolbox # \\ / O peration | -# \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation +# \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation # \\/ M anipulation | #------------------------------------------------------------------------------ # License @@ -107,15 +107,23 @@ then echo "pack manually" 1>&2 foamPackSource $packDir $packFile else - echo "pack with git-archive" 1>&2 - ( cd $packDir && git archive --format=tar --prefix=$packDir/ HEAD) > $packBase.tar + if [ ! -f $packDir/.build ] + then + echo "Error: $packDir/.build does not exists" 1>&2 + echo " Please update this by running e.g. 'wmakePrintBuild -update' in $packDir" 1>&2 + exit 2 + fi - echo "add in time-stamp and lnInclude directories" 1>&2 - tar cf $packBase.tar2 $packDir/.timeStamp $packDir/.build `find -H $packDir -type d -name lnInclude` - tar Af $packBase.tar $packBase.tar2 - echo "gzip tar file" 1>&2 - gzip -c9 $packBase.tar > $packFile + echo "pack with git-archive" 1>&2 && + ( cd $packDir && git archive --format=tar --prefix=$packDir/ HEAD) > $packBase.tar && + + echo "add in time-stamp and lnInclude directories" 1>&2 && + tar cf $packBase.tar2 $packDir/.timeStamp $packDir/.build `find -H $packDir -type d -name lnInclude` && + tar Af $packBase.tar $packBase.tar2 && + + echo "gzip tar file" 1>&2 && + gzip -c9 $packBase.tar > $packFile && rm -f $packBase.tar $packBase.tar2 2>/dev/null fi From 5a80198bef39013b51b600c4a0c836d04230992f Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 23 May 2012 15:11:55 +0100 Subject: [PATCH 06/42] ENH: Updates to rotor momentum source --- .../rotorDiskSource/bladeModel/bladeModel.C | 49 ++++--- .../profileModel/lookup/lookupProfile.C | 37 +++--- .../rotorDiskSource/rotorDiskSource.C | 120 ++++++++---------- .../rotorDiskSource/rotorDiskSource.H | 3 + .../trimModel/fixed/fixedTrim.C | 5 - .../trimModel/targetForce/targetForceTrim.C | 5 - 6 files changed, 100 insertions(+), 119 deletions(-) diff --git a/src/fieldSources/basicSource/rotorDiskSource/bladeModel/bladeModel.C b/src/fieldSources/basicSource/rotorDiskSource/bladeModel/bladeModel.C index a1f76d3e99..af12f6b0ce 100644 --- a/src/fieldSources/basicSource/rotorDiskSource/bladeModel/bladeModel.C +++ b/src/fieldSources/basicSource/rotorDiskSource/bladeModel/bladeModel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -47,33 +47,34 @@ void Foam::bladeModel::interpolateWeights scalar& ddx ) const { - scalar x = -GREAT; + i2 = 0; label nElem = values.size(); - i2 = 0; - while ((x < xIn) && (i2 < nElem)) + if (nElem == 1) { - x = values[i2]; - i2++; - } - - if (i2 == 0) - { - i1 = i2; - ddx = 0.0; - return; - } - else if (i2 == values.size()) - { - i2 = values.size() - 1; i1 = i2; ddx = 0.0; return; } else { - i1 = i2 - 1; - ddx = (xIn - values[i1])/(values[i2] - values[i1]); + while ((values[i2] < xIn) && (i2 < nElem)) + { + i2++; + } + + if (i2 == nElem) + { + i2 = nElem - 1; + i1 = i2; + ddx = 0.0; + return; + } + else + { + i1 = i2 - 1; + ddx = (xIn - values[i1])/(values[i2] - values[i1]); + } } } @@ -120,14 +121,8 @@ Foam::bladeModel::bladeModel(const dictionary& dict) } else { - FatalErrorIn - ( - "Foam::bladeModel::bladeModel" - "(" - "const dictionary&, " - "const word&" - ")" - ) << "No blade data specified" << exit(FatalError); + FatalErrorIn("Foam::bladeModel::bladeModel(const dictionary&)") + << "No blade data specified" << exit(FatalError); } } diff --git a/src/fieldSources/basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.C b/src/fieldSources/basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.C index 1dda6b403b..93e54bc8ef 100644 --- a/src/fieldSources/basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.C +++ b/src/fieldSources/basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.C @@ -49,33 +49,34 @@ void Foam::lookupProfile::interpolateWeights scalar& ddx ) const { - scalar x = -GREAT; + i2 = 0; label nElem = values.size(); - i2 = 0; - while ((x < xIn) && (i2 < nElem)) + if (nElem == 1) { - x = values[i2]; - i2++; - } - - if (i2 == 0) - { - i1 = i2; - ddx = 0.0; - return; - } - else if (i2 == nElem) - { - i2 = nElem - 1; i1 = i2; ddx = 0.0; return; } else { - i1 = i2 - 1; - ddx = (xIn - values[i1])/(values[i2] - values[i1]); + while ((values[i2] < xIn) && (i2 < nElem)) + { + i2++; + } + + if (i2 == nElem) + { + i2 = nElem - 1; + i1 = i2; + ddx = 0.0; + return; + } + else + { + i1 = i2 - 1; + ddx = (xIn - values[i1])/(values[i2] - values[i1]); + } } } diff --git a/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C b/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C index 20161147df..6014b623fb 100644 --- a/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C +++ b/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C @@ -67,6 +67,7 @@ namespace Foam void Foam::rotorDiskSource::checkData() { + // set inflow type switch (selectionMode()) { case smCellSet: @@ -129,11 +130,10 @@ void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct) const label nInternalFaces = mesh_.nInternalFaces(); const polyBoundaryMesh& pbm = mesh_.boundaryMesh(); - const vectorField& Sf = mesh_.Sf().internalField(); - const scalarField& magSf = mesh_.magSf().internalField(); + const vectorField& Sf = mesh_.Sf(); + const scalarField& magSf = mesh_.magSf(); vector n = vector::zero; - label nFace = 0; // calculate cell addressing for selected cells labelList cellAddr(mesh_.nCells(), -1); @@ -158,7 +158,6 @@ void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct) // correct for parallel running syncTools::swapBoundaryFaceList(mesh_, nbrFaceCellAddr); - // add internal field contributions for (label faceI = 0; faceI < nInternalFaces; faceI++) { @@ -173,7 +172,6 @@ void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct) { area_[own] += magSf[faceI]; n += Sf[faceI]; - nFace++; } } else if ((own == -1) && (nbr != -1)) @@ -184,7 +182,6 @@ void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct) { area_[nbr] += magSf[faceI]; n -= Sf[faceI]; - nFace++; } } } @@ -210,7 +207,6 @@ void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct) { area_[own] += magSfp[j]; n += Sfp[j]; - nFace++; } } } @@ -222,11 +218,10 @@ void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct) const label own = cellAddr[mesh_.faceOwner()[faceI]]; const vector nf = Sfp[j]/magSfp[j]; - if ((own != -1) && ((nf & axis) > 0.8)) + if ((own != -1) && ((nf & axis) > tol)) { area_[own] += magSfp[j]; n += Sfp[j]; - nFace++; } } } @@ -359,32 +354,30 @@ void Foam::rotorDiskSource::constructGeometry() forAll(cells_, i) { - const label cellI = cells_[i]; - - // position in (planar) rotor co-ordinate system - x_[i] = coordSys_.localPosition(C[cellI]); - - // cache max radius - rMax_ = max(rMax_, x_[i].x()); - - // swept angle relative to rDir axis [radians] in range 0 -> 2*pi - scalar psi = x_[i].y(); - if (psi < 0) + if (area_[i] > ROOTVSMALL) { - psi += mathematical::twoPi; + const label cellI = cells_[i]; + + // position in (planar) rotor co-ordinate system + x_[i] = coordSys_.localPosition(C[cellI]); + + // cache max radius + rMax_ = max(rMax_, x_[i].x()); + + // swept angle relative to rDir axis [radians] in range 0 -> 2*pi + scalar psi = x_[i].y(); + + // blade flap angle [radians] + scalar beta = + flap_.beta0 - flap_.beta1*cos(psi) - flap_.beta2*sin(psi); + + // determine rotation tensor to convert from planar system into the + // rotor cone system + scalar c = cos(beta); + scalar s = sin(beta); + R_[i] = tensor(c, 0, -s, 0, 1, 0, s, 0, c); + invR_[i] = R_[i].T(); } - - // blade flap angle [radians] - scalar beta = flap_.beta0 - flap_.beta1*cos(psi) - flap_.beta2*sin(psi); - - // determine rotation tensor to convert from planar system into the - // rotor cone system - scalar cNeg = cos(-beta); - scalar sNeg = sin(-beta); - R_[i] = tensor(cNeg, 0.0, -sNeg, 0.0, 1.0, 0.0, sNeg, 0.0, cNeg); - scalar cPos = cos(beta); - scalar sPos = sin(beta); - invR_[i] = tensor(cPos, 0.0, -sPos, 0.0, 1.0, 0.0, sPos, 0.0, cPos); } } @@ -440,6 +433,7 @@ Foam::rotorDiskSource::rotorDiskSource : basicSource(name, modelType, dict, mesh), rhoName_("none"), + rhoRef_(1.0), omega_(0.0), nBlades_(0), inletFlow_(ifLocal), @@ -477,13 +471,15 @@ void Foam::rotorDiskSource::calculate const bool output ) const { - tmp trho; - if (rhoName_ != "none") - { - trho = mesh_.lookupObject(rhoName_); - } - const scalarField& V = mesh_.V(); + const bool compressible = rhoName_ != "none"; + + tmp trho + ( + compressible + ? mesh_.lookupObject(rhoName_) + : volScalarField::null() + ); // logging info @@ -503,17 +499,17 @@ void Foam::rotorDiskSource::calculate // velocity in local cylindrical reference frame vector Uc = coordSys_.localVector(U[cellI]); - // apply correction in local system due to coning + // transform from rotor cylindrical into local coning system Uc = R_[i] & Uc; // set radial component of velocity to zero Uc.x() = 0.0; - // remove blade linear velocity from blade normal component - Uc.y() -= radius*omega_; + // set blade normal component of velocity + Uc.y() = radius*omega_ - Uc.y(); // determine blade data for this radius - // i1 = index of upper bound data point in blade list + // i2 = index of upper radius bound data point in blade list scalar twist = 0.0; scalar chord = 0.0; label i1 = -1; @@ -529,17 +525,7 @@ void Foam::rotorDiskSource::calculate } // effective angle of attack - scalar alphaEff = - mathematical::pi + atan2(Uc.z(), Uc.y()) - alphaGeom; - - if (alphaEff > mathematical::pi) - { - alphaEff -= mathematical::twoPi; - } - if (alphaEff < -mathematical::pi) - { - alphaEff += mathematical::twoPi; - } + scalar alphaEff = alphaGeom - atan2(-Uc.z(), Uc.y()); AOAmin = min(AOAmin, alphaEff); AOAmax = max(AOAmax, alphaEff); @@ -564,21 +550,21 @@ void Foam::rotorDiskSource::calculate // calculate forces perpendicular to blade scalar pDyn = 0.5*magSqr(Uc); - if (trho.valid()) + if (compressible) { pDyn *= trho()[cellI]; } - scalar f = pDyn*chord*nBlades_*area_[i]/mathematical::twoPi; - vector localForce = vector(0.0, f*Cd, tipFactor*f*Cl); + scalar f = pDyn*chord*nBlades_*area_[i]/radius/mathematical::twoPi; + vector localForce = vector(0.0, -f*Cd, tipFactor*f*Cl); + + // accumulate forces + dragEff += rhoRef_*localForce.y(); + liftEff += rhoRef_*localForce.z(); // convert force from local coning system into rotor cylindrical localForce = invR_[i] & localForce; - // accumulate forces - dragEff += localForce.y(); - liftEff += localForce.z(); - // convert force to global cartesian co-ordinate system force[cellI] = coordSys_.globalVector(localForce); @@ -616,6 +602,7 @@ void Foam::rotorDiskSource::addSup(fvMatrix& eqn, const label fieldI) } else { + coeffs_.lookup("rhoRef") >> rhoRef_; dims.reset(dimForce/dimVolume/dimDensity); } @@ -666,6 +653,8 @@ bool Foam::rotorDiskSource::read(const dictionary& dict) coeffs_.lookup("fieldNames") >> fieldNames_; applied_.setSize(fieldNames_.size(), false); + + // read co-ordinate system/geometry invariant properties scalar rpm(readScalar(coeffs_.lookup("rpm"))); omega_ = rpm/60.0*mathematical::twoPi; @@ -683,14 +672,17 @@ bool Foam::rotorDiskSource::read(const dictionary& dict) flap_.beta1 = degToRad(flap_.beta1); flap_.beta2 = degToRad(flap_.beta2); - trim_->read(coeffs_); - - checkData(); + // create co-ordinate system createCoordinateSystem(); + // read co-odinate system dependent properties + checkData(); + constructGeometry(); + trim_->read(coeffs_); + if (debug) { writeField("alphag", trim_->thetag()(), true); diff --git a/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H b/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H index 8e26425e10..a2f856988d 100644 --- a/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H +++ b/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H @@ -148,6 +148,9 @@ protected: //- Name of density field word rhoName_; + //- Reference density if rhoName = 'none' + scalar rhoRef_; + //- Rotational speed [rad/s] // Positive anti-clockwise when looking along -ve lift direction scalar omega_; diff --git a/src/fieldSources/basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.C b/src/fieldSources/basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.C index d5f4843d58..28b208c7c9 100644 --- a/src/fieldSources/basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.C +++ b/src/fieldSources/basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.C @@ -71,11 +71,6 @@ void Foam::fixedTrim::read(const dictionary& dict) forAll(thetag_, i) { scalar psi = x[i].y(); - if (psi < 0) - { - psi += mathematical::twoPi; - } - thetag_[i] = theta0 + theta1c*cos(psi) + theta1s*sin(psi); } } diff --git a/src/fieldSources/basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.C b/src/fieldSources/basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.C index 409015d885..f22226f84c 100644 --- a/src/fieldSources/basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.C +++ b/src/fieldSources/basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.C @@ -142,11 +142,6 @@ Foam::tmp Foam::targetForceTrim::thetag() const forAll(t, i) { scalar psi = x[i].y(); - if (psi < 0) - { - psi += mathematical::twoPi; - } - t[i] = theta_[0] + theta_[1]*cos(psi) + theta_[2]*sin(psi); } From 89bcfe0441861a6a39fd1035b63963fcf1a540bd Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 23 May 2012 18:12:39 +0100 Subject: [PATCH 07/42] BUG: ptScotchDecomp: handle zero-sized domains --- .../ptscotchDecomp/dummyPtscotchDecomp.C | 20 +- .../decompose/ptscotchDecomp/ptscotchDecomp.C | 494 +++++++++--------- .../decompose/ptscotchDecomp/ptscotchDecomp.H | 26 +- 3 files changed, 278 insertions(+), 262 deletions(-) diff --git a/src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.C b/src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.C index 701de4e23a..36d405d597 100644 --- a/src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.C +++ b/src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -58,7 +58,7 @@ void Foam::ptscotchDecomp::check(const int retVal, const char* str) {} -Foam::label Foam::ptscotchDecomp::decomposeZeroDomains +Foam::label Foam::ptscotchDecomp::decompose ( const fileName& meshPath, const List& initxadj, @@ -72,6 +72,7 @@ Foam::label Foam::ptscotchDecomp::decomposeZeroDomains ( "label ptscotchDecomp::decompose" "(" + "onst fileName&," "const List&, " "const List&, " "const scalarField&, " @@ -84,8 +85,10 @@ Foam::label Foam::ptscotchDecomp::decomposeZeroDomains Foam::label Foam::ptscotchDecomp::decompose ( const fileName& meshPath, - const List& adjncy, - const List& xadj, + const int adjncySize, + const int adjncy[], + const int xadjSize, + const int xadj[], const scalarField& cWeights, List& finalDecomp ) const @@ -94,9 +97,12 @@ Foam::label Foam::ptscotchDecomp::decompose ( "label ptscotchDecomp::decompose" "(" - "const List&, " - "const List&, " - "const scalarField&, " + "const fileName&," + "const int," + "const int," + "const int," + "const int," + "const scalarField&," "List&" ")" ) << notImplementedMessage << exit(FatalError); diff --git a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C index e7c4995a36..81d3851837 100644 --- a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C +++ b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -103,9 +103,9 @@ License - Note: writes out .dgr files for debugging. Run with e.g. + Note: writeGraph=true : writes out .dgr files for debugging. Run with e.g. - mpirun -np 4 dgpart 2 'processor%r.grf' + mpirun -np 4 dgpart 2 'region0_%r.dgr' - %r gets replaced by current processor rank - decompose into 2 domains @@ -167,192 +167,192 @@ void Foam::ptscotchDecomp::check(const int retVal, const char* str) } -//- Does prevention of 0 cell domains and calls ptscotch. -Foam::label Foam::ptscotchDecomp::decomposeZeroDomains -( - const fileName& meshPath, - const List& initadjncy, - const List& initxadj, - const scalarField& initcWeights, - - List& finalDecomp -) const -{ - globalIndex globalCells(initxadj.size()-1); - - bool hasZeroDomain = false; - for (label procI = 0; procI < Pstream::nProcs(); procI++) - { - if (globalCells.localSize(procI) == 0) - { - hasZeroDomain = true; - break; - } - } - - if (!hasZeroDomain) - { - return decompose - ( - meshPath, - initadjncy, - initxadj, - initcWeights, - finalDecomp - ); - } - - - if (debug) - { - Info<< "ptscotchDecomp : have graphs with locally 0 cells." - << " trickling down." << endl; - } - - // Make sure every domain has at least one cell - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - // (scotch does not like zero sized domains) - // Trickle cells from processors that have them up to those that - // don't. - - - // Number of cells to send to the next processor - // (is same as number of cells next processor has to receive) - List nSendCells(Pstream::nProcs(), 0); - - for (label procI = nSendCells.size()-1; procI >=1; procI--) - { - label nLocalCells = globalCells.localSize(procI); - if (nLocalCells-nSendCells[procI] < 1) - { - nSendCells[procI-1] = nSendCells[procI]-nLocalCells+1; - } - } - - // First receive (so increasing the sizes of all arrays) - - Field xadj(initxadj); - Field adjncy(initadjncy); - scalarField cWeights(initcWeights); - - if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0) - { - // Receive cells from previous processor - IPstream fromPrevProc(Pstream::blocking, Pstream::myProcNo()-1); - - Field prevXadj(fromPrevProc); - Field prevAdjncy(fromPrevProc); - scalarField prevCellWeights(fromPrevProc); - - if (prevXadj.size() != nSendCells[Pstream::myProcNo()-1]) - { - FatalErrorIn("ptscotchDecomp::decompose(..)") - << "Expected from processor " << Pstream::myProcNo()-1 - << " connectivity for " << nSendCells[Pstream::myProcNo()-1] - << " nCells but only received " << prevXadj.size() - << abort(FatalError); - } - - // Insert adjncy - prepend(prevAdjncy, adjncy); - // Adapt offsets and prepend xadj - xadj += prevAdjncy.size(); - prepend(prevXadj, xadj); - // Weights - prepend(prevCellWeights, cWeights); - } - - - // Send to my next processor - - if (nSendCells[Pstream::myProcNo()] > 0) - { - // Send cells to next processor - OPstream toNextProc(Pstream::blocking, Pstream::myProcNo()+1); - - int nCells = nSendCells[Pstream::myProcNo()]; - int startCell = xadj.size()-1 - nCells; - int startFace = xadj[startCell]; - int nFaces = adjncy.size()-startFace; - - // Send for all cell data: last nCells elements - // Send for all face data: last nFaces elements - toNextProc - << Field::subField(xadj, nCells, startCell)-startFace - << Field::subField(adjncy, nFaces, startFace) - << - ( - cWeights.size() - ? static_cast - ( - scalarField::subField(cWeights, nCells, startCell) - ) - : scalarField(0) - ); - - // Remove data that has been sent - if (cWeights.size()) - { - cWeights.setSize(cWeights.size()-nCells); - } - adjncy.setSize(adjncy.size()-nFaces); - xadj.setSize(xadj.size() - nCells); - } - - - // Do decomposition as normal. Sets finalDecomp. - label result = decompose(meshPath, adjncy, xadj, cWeights, finalDecomp); - - - if (debug) - { - Info<< "ptscotchDecomp : have graphs with locally 0 cells." - << " trickling up." << endl; - } - - - // If we sent cells across make sure we undo it - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - // Receive back from next processor if I sent something - if (nSendCells[Pstream::myProcNo()] > 0) - { - IPstream fromNextProc(Pstream::blocking, Pstream::myProcNo()+1); - - List nextFinalDecomp(fromNextProc); - - if (nextFinalDecomp.size() != nSendCells[Pstream::myProcNo()]) - { - FatalErrorIn("parMetisDecomp::decompose(..)") - << "Expected from processor " << Pstream::myProcNo()+1 - << " decomposition for " << nSendCells[Pstream::myProcNo()] - << " nCells but only received " << nextFinalDecomp.size() - << abort(FatalError); - } - - append(nextFinalDecomp, finalDecomp); - } - - // Send back to previous processor. - if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0) - { - OPstream toPrevProc(Pstream::blocking, Pstream::myProcNo()-1); - - int nToPrevious = nSendCells[Pstream::myProcNo()-1]; - - toPrevProc << - SubList - ( - finalDecomp, - nToPrevious, - finalDecomp.size()-nToPrevious - ); - - // Remove locally what has been sent - finalDecomp.setSize(finalDecomp.size()-nToPrevious); - } - return result; -} +////- Does prevention of 0 cell domains and calls ptscotch. +//Foam::label Foam::ptscotchDecomp::decomposeZeroDomains +//( +// const fileName& meshPath, +// const List& initadjncy, +// const List& initxadj, +// const scalarField& initcWeights, +// +// List& finalDecomp +//) const +//{ +// globalIndex globalCells(initxadj.size()-1); +// +// bool hasZeroDomain = false; +// for (label procI = 0; procI < Pstream::nProcs(); procI++) +// { +// if (globalCells.localSize(procI) == 0) +// { +// hasZeroDomain = true; +// break; +// } +// } +// +// if (!hasZeroDomain) +// { +// return decompose +// ( +// meshPath, +// initadjncy, +// initxadj, +// initcWeights, +// finalDecomp +// ); +// } +// +// +// if (debug) +// { +// Info<< "ptscotchDecomp : have graphs with locally 0 cells." +// << " trickling down." << endl; +// } +// +// // Make sure every domain has at least one cell +// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// // (scotch does not like zero sized domains) +// // Trickle cells from processors that have them up to those that +// // don't. +// +// +// // Number of cells to send to the next processor +// // (is same as number of cells next processor has to receive) +// List nSendCells(Pstream::nProcs(), 0); +// +// for (label procI = nSendCells.size()-1; procI >=1; procI--) +// { +// label nLocalCells = globalCells.localSize(procI); +// if (nLocalCells-nSendCells[procI] < 1) +// { +// nSendCells[procI-1] = nSendCells[procI]-nLocalCells+1; +// } +// } +// +// // First receive (so increasing the sizes of all arrays) +// +// Field xadj(initxadj); +// Field adjncy(initadjncy); +// scalarField cWeights(initcWeights); +// +// if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0) +// { +// // Receive cells from previous processor +// IPstream fromPrevProc(Pstream::blocking, Pstream::myProcNo()-1); +// +// Field prevXadj(fromPrevProc); +// Field prevAdjncy(fromPrevProc); +// scalarField prevCellWeights(fromPrevProc); +// +// if (prevXadj.size() != nSendCells[Pstream::myProcNo()-1]) +// { +// FatalErrorIn("ptscotchDecomp::decompose(..)") +// << "Expected from processor " << Pstream::myProcNo()-1 +// << " connectivity for " << nSendCells[Pstream::myProcNo()-1] +// << " nCells but only received " << prevXadj.size() +// << abort(FatalError); +// } +// +// // Insert adjncy +// prepend(prevAdjncy, adjncy); +// // Adapt offsets and prepend xadj +// xadj += prevAdjncy.size(); +// prepend(prevXadj, xadj); +// // Weights +// prepend(prevCellWeights, cWeights); +// } +// +// +// // Send to my next processor +// +// if (nSendCells[Pstream::myProcNo()] > 0) +// { +// // Send cells to next processor +// OPstream toNextProc(Pstream::blocking, Pstream::myProcNo()+1); +// +// int nCells = nSendCells[Pstream::myProcNo()]; +// int startCell = xadj.size()-1 - nCells; +// int startFace = xadj[startCell]; +// int nFaces = adjncy.size()-startFace; +// +// // Send for all cell data: last nCells elements +// // Send for all face data: last nFaces elements +// toNextProc +// << Field::subField(xadj, nCells, startCell)-startFace +// << Field::subField(adjncy, nFaces, startFace) +// << +// ( +// cWeights.size() +// ? static_cast +// ( +// scalarField::subField(cWeights, nCells, startCell) +// ) +// : scalarField(0) +// ); +// +// // Remove data that has been sent +// if (cWeights.size()) +// { +// cWeights.setSize(cWeights.size()-nCells); +// } +// adjncy.setSize(adjncy.size()-nFaces); +// xadj.setSize(xadj.size() - nCells); +// } +// +// +// // Do decomposition as normal. Sets finalDecomp. +// label result = decompose(meshPath, adjncy, xadj, cWeights, finalDecomp); +// +// +// if (debug) +// { +// Info<< "ptscotchDecomp : have graphs with locally 0 cells." +// << " trickling up." << endl; +// } +// +// +// // If we sent cells across make sure we undo it +// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// +// // Receive back from next processor if I sent something +// if (nSendCells[Pstream::myProcNo()] > 0) +// { +// IPstream fromNextProc(Pstream::blocking, Pstream::myProcNo()+1); +// +// List nextFinalDecomp(fromNextProc); +// +// if (nextFinalDecomp.size() != nSendCells[Pstream::myProcNo()]) +// { +// FatalErrorIn("parMetisDecomp::decompose(..)") +// << "Expected from processor " << Pstream::myProcNo()+1 +// << " decomposition for " << nSendCells[Pstream::myProcNo()] +// << " nCells but only received " << nextFinalDecomp.size() +// << abort(FatalError); +// } +// +// append(nextFinalDecomp, finalDecomp); +// } +// +// // Send back to previous processor. +// if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0) +// { +// OPstream toPrevProc(Pstream::blocking, Pstream::myProcNo()-1); +// +// int nToPrevious = nSendCells[Pstream::myProcNo()-1]; +// +// toPrevProc << +// SubList +// ( +// finalDecomp, +// nToPrevious, +// finalDecomp.size()-nToPrevious +// ); +// +// // Remove locally what has been sent +// finalDecomp.setSize(finalDecomp.size()-nToPrevious); +// } +// return result; +//} // Call scotch with options from dictionary. @@ -362,13 +362,42 @@ Foam::label Foam::ptscotchDecomp::decompose const List& adjncy, const List& xadj, const scalarField& cWeights, + List& finalDecomp +) const +{ + List dummyAdjncy(1); + List dummyXadj(1); + dummyXadj[0] = 0; + + return decompose + ( + meshPath, + adjncy.size(), + (adjncy.size() ? adjncy.begin() : dummyAdjncy.begin()), + xadj.size(), + (xadj.size() ? xadj.begin() : dummyXadj.begin()), + cWeights, + finalDecomp + ); +} + + +// Call scotch with options from dictionary. +Foam::label Foam::ptscotchDecomp::decompose +( + const fileName& meshPath, + const int adjncySize, + const int adjncy[], + const int xadjSize, + const int xadj[], + const scalarField& cWeights, List& finalDecomp ) const { if (debug) { - Pout<< "ptscotchDecomp : entering with xadj:" << xadj.size() << endl; + Pout<< "ptscotchDecomp : entering with xadj:" << xadjSize << endl; } // Dump graph @@ -387,7 +416,7 @@ Foam::label Foam::ptscotchDecomp::decompose Pout<< "Dumping Scotch graph file to " << str.name() << endl << "Use this in combination with dgpart." << endl; - globalIndex globalCells(xadj.size()-1); + globalIndex globalCells(xadjSize-1); // Distributed graph file (.grf) label version = 2; @@ -400,17 +429,17 @@ Foam::label Foam::ptscotchDecomp::decompose // Total number of vertices (vertglbnbr) str << globalCells.size(); // Total number of connections (edgeglbnbr) - str << ' ' << returnReduce(xadj[xadj.size()-1], sumOp