Commit 026598df authored by henry's avatar henry
Browse files
parents bcfa40c4 febecdcc
Short overview of the changes to have cyclics split into two halves.
Cyclics
-------
The two cyclic halves are now split like processor patches. There should be no
difference in running.
Advantages:
- decomposed cyclics can now be handled properly. It just needs to preserve
the cyclic patch it originates from.
- We can now construct a table of global transformations and handle
points/edges/cells with transformations.
- face ordering after topological changes becomes much easier since we
now preserve what half the face comes from.
- cyclic handling becomes more consistent with processor handling and can
quite often be handled in the same condition.
- transformation tensors now become single entry.
The disadvantages:
- a patch-wise loop now might need to store data to go to the neighbour half
since it is no longer handled in a single patch.
- decomposed cyclics now require overlapping communications so will
only work in non-blocking mode. Hence the underlying message passing library
will require overlapping communications with message tags.
- it is quite a code-change and there might be some oversights.
- once converted (see foamUpgradeCyclics below) cases are not backwards
compatible with previous versions.
blockMesh
---------
blockMeshDict now allows patch definition using the construct-from-dictionary
constructor. This helps defining patches that require additional input e.g.
directMapped and now cyclic:
boundary
(
sides2_half0
{
type cyclic;
neighbourPatch sides2_half1;
faces ((2 4 5 3));
}
The syntax is - like the polyMesh/boundary file - a list of dictionaries with
one additional entry 'faces' for the block faces. Above shows the new
required entry 'neighbourPatch' for cyclic.
blockMesh still reads the old format. For a cyclic it will automatically
introduce two patches for the halves, with names xxx_half0 and xxx_half1.
foamUpgradeCyclics
------------------
This is a tool which reads the polyMesh/boundary file and any vol/surface/point
fields and converts them.
It will check if anything needs to be converted, backup the current file to .old
and split any cyclic patchFields into two entries.
decomposePar
------------
Decomposes cyclics into processorCyclic:
procBoundary0to1throughsides1_half0
{
type processorCyclic;
nFaces 1000;
startFace 91350;
myProcNo 0;
neighbProcNo 1;
referPatch sides1_half0;
}
They have an additional 'referPatch' entry which gives the (cyclic) patch
to use for any transformation.
Details
-------
- the cyclic patch dictionary has an entry neighbourPatch. The
patch has new member functions:
//- Get neighbouring patchID
label neighbPatchID() const
//- Get neighbouring patch
const cyclicPolyPatch& neighbPatch()
//- Am I the owner half
bool owner()
The cyclic still has forward() and reverse() transformations (with
the reverse() equal to the neighbPatch().forward()).
There is no transformLocalFace anymore - the ordering is the same for
both halves.
- 'pure' processor patches now are always coincident - they (should) have no
transformation. As said above cyclics are decomposed into a derived
type 'processorCyclic'.
- processor patches use overlapping communication using a different message
tag. This maps straight through into the MPI message tag.
- when constructing a GeometricField from a dictionary it will explicitly
check for non-existing entries for cyclic patches and exit with an error message
warning to run foamUpgradeCyclics.
- allocate/free tags.
Not tested.
OK - test blockMesh with cyclics
unitTestCases/singleCyclic/
OK - test cyclics sequential running.
unitTestCases/singleCyclic/
OK - test decomposePar
tested channel395
OK - FaceCellWave
unitTestCases/twoCavityCyclicForWallDistance/
OK - non-parallel finite volume : channelFoam
unitTestCases/channel395-splitCyclic vs. channel395-dev
OK - parallel finite volume with processorCyclic: channelFoam
unitTestCases/channel395-splitCyclic vs. channel395-dev
OK - preProcessing/foamUpgradeCyclics
OK - gamg - sequential.
Tested on channel395-splitCyclic with GAMG.
OK - gamg parallel.
Tested on channel395-splitCyclic with GAMG.
- initTransfer in GAMGprocessorInterfaces using nonblocking+tags
untested.
OK - cyclic baffles.
Tested on t-junction-with-fan
OK. - jumpCyclics/fanFvPatchField. All usages of jump() now need to account
for being owner() or not.
Tested on t-junction-with-fan.
OK - regionSplit
tested on singleCyclic
OK - pointFields on cyclics. volPointInterpolation.
tested on channel395-splitCyclic
OK - fvMeshSubset
tested on singleCyclic
OK - pointEdgeWave (maybe test through inversePointDistanceDiffusivity?)
tested on twoCavityCyclicForWallDistance
OK - scotchDecomp
tested with testCalcCSR on twoCavityCyclicForWallDistance
NOT WORKING - fvMeshDistribute to split cyclic patches into ones
with different separation.
tested on singleCyclic
OK - test createPatch pointSync
note: only works if face-centre position of 0th faces is ok since uses
this for separation. Should in fact make cyclic planar using patch centre and
normal?
test on twoCavityCyclicForWallDistance with point (0 1 0) set to (0 1.001 0)
NO PROBLEM - renumberMesh
It doesn't do renumbering through cyclics.
OK - rotational cyclics.
Tested on movingCone-with-cyclics
OK - LUscalarMatrix::convert still expects interfaces to be cyclic
tested on channel395 with 'directSolveCoarsest true;'
OK - grep for size()/2
- all tutorials with cyclics:
OK - incompressible/DNS/dnsFoam/boxTurb16
OK - incompressible/channelFoam/channel395
slight differences due to divergence. combustion/XiFoam/les/pitzDaily3D
OK - no cyclics. combustion/fireFoam/les/smallPoolFire2D
discreteMethods/dsmcFoam/freeSpacePeriodic
discreteMethods/dsmcFoam/wedge15Ma5
discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCubeArgon
OK - incompressible/boundaryFoam/boundaryLaunderSharma
OK - incompressible/boundaryFoam/boundaryWallFunctions
OK - incompressible/boundaryFoam/boundaryWallFunctionsProfile
OK - needs createBaffles. incompressible/pimpleFoam/t-junction-with-fan
OK - incompressible/simpleSRFFoam/mixer
OK - needs createBaffles. lagrangian/porousExplicitSourceReactingParcelFoam/filter
needs special coupledbcs. lagrangian/reactingParcelFilmFoam/multipleBoxes
OK - createBaffles
- have foamUpgradeCyclics split 'value' field
- activeBaffleVelocity
- kivaToFoam/readKivaGrid.H sorts cyclics (but in incorrect order?)
- isoSurface.C
- referredCellList.C
- work out scheduled communication?
OK - add neighbourPatch checking to 16x.
Info<< "Reading transportProperties\n" << endl;
Info<< "\nReading transportProperties\n" << endl;
IOdictionary transportProperties
(
......@@ -31,7 +31,7 @@
rhoInfValue
);
Info<< "\nReading field U\n" << endl;
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
......@@ -127,7 +127,7 @@
if (HdotGradHheader.headerOk())
{
Info<< "\nReading field HdotGradH\n" << endl;
Info<< "Reading field HdotGradH" << endl;
HdotGradHPtr_.reset
(
......
......@@ -295,8 +295,7 @@ label mergePatchFaces
const faceZone& fZone = mesh.faceZones()[zoneID];
zoneFlip = fZone.flipMap()[fZone.whichFace(newMasterI)];
}
label patchI = mesh.boundaryMesh().whichPatch(newMasterI);
label patchID = mesh.boundaryMesh().whichPatch(newMasterI);
Pout<< "Restoring new master face " << newMasterI
<< " to vertices " << setFaceVerts[0] << endl;
......@@ -311,7 +310,7 @@ label mergePatchFaces
own, // owner
-1, // neighbour
false, // face flip
patchI, // patch for face
patchID, // patch for face
false, // remove from zone
zoneID, // zone for face
zoneFlip // face flip in zone
......@@ -336,7 +335,7 @@ label mergePatchFaces
-1, // masterEdgeID,
newMasterI, // masterFaceID,
false, // flipFaceFlux,
patchI, // patchID,
patchID, // patchID,
zoneID, // zoneID,
zoneFlip // zoneFlip
)
......
......@@ -206,25 +206,8 @@ int main(int argc, char *argv[])
Info<< nl << "Creating polyMesh from blockMesh" << endl;
wordList patchNames = blocks.patchNames();
wordList patchTypes = blocks.patchTypes();
word defaultFacesName = "defaultFaces";
word defaultFacesType = emptyPolyPatch::typeName;
wordList patchPhysicalTypes = blocks.patchPhysicalTypes();
preservePatchTypes
(
runTime,
runTime.constant(),
polyMeshDir,
patchNames,
patchTypes,
defaultFacesName,
defaultFacesType,
patchPhysicalTypes
);
polyMesh mesh
(
IOobject
......@@ -236,11 +219,9 @@ int main(int argc, char *argv[])
xferCopy(blocks.points()), // could we re-use space?
blocks.cells(),
blocks.patches(),
patchNames,
patchTypes,
blocks.patchDicts(),
defaultFacesName,
defaultFacesType,
patchPhysicalTypes
defaultFacesType
);
......
......@@ -158,7 +158,6 @@ void Foam::createShellMesh::calcPointRegions
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::createShellMesh::createShellMesh
......@@ -184,7 +183,6 @@ Foam::createShellMesh::createShellMesh
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::createShellMesh::setRefinement
(
const pointField& thickness,
......
......@@ -130,6 +130,7 @@ Usage
#include "volFields.H"
#include "surfaceFields.H"
#include "cyclicPolyPatch.H"
#include "syncTools.H"
using namespace Foam;
......@@ -595,6 +596,237 @@ void createDummyFvMeshFiles(const polyMesh& mesh, const word& regionName)
}
//XXXXXXXXX
label findUncoveredPatchFace
(
const fvMesh& mesh,
const UIndirectList<label>& extrudeMeshFaces,// mesh faces that are extruded
const label meshEdgeI // mesh edge
)
{
// Make set of extruded faces.
labelHashSet extrudeFaceSet(extrudeMeshFaces.size());
forAll(extrudeMeshFaces, i)
{
extrudeFaceSet.insert(extrudeMeshFaces[i]);
}
label patchI = -1;
const labelList& eFaces = mesh.edgeFaces()[meshEdgeI];
forAll(eFaces, i)
{
label faceI = eFaces[i];
if (!mesh.isInternalFace(faceI) && !extrudeFaceSet.found(faceI))
{
patchI = mesh.boundaryMesh().whichPatch(faceI);
break;
}
}
return patchI;
}
// Count the number of faces in patches that need to be created
void countExtrudePatches
(
const fvMesh& mesh,
const primitiveFacePatch& extrudePatch,
const label nZones,
const labelList& zoneID,
const labelList& extrudeMeshFaces,
const labelList& extrudeMeshEdges,
labelList& zoneSidePatch,
labelList& zoneZonePatch
)
{
const labelListList& edgeFaces = extrudePatch.edgeFaces();
forAll(edgeFaces, edgeI)
{
const labelList& eFaces = edgeFaces[edgeI];
if (eFaces.size() == 2)
{
label zone0 = zoneID[eFaces[0]];
label zone1 = zoneID[eFaces[1]];
if (zone0 != zone1)
{
label minZone = min(zone0,zone1);
label maxZone = max(zone0,zone1);
zoneZonePatch[minZone*nZones+maxZone]++;
}
}
else
{
// Check whether we are on a mesh edge with external patches. If
// so choose any uncovered one. If none found put face in
// undetermined zone 'side' patch
label patchI = findUncoveredPatchFace
(
mesh,
UIndirectList<label>(extrudeMeshFaces, eFaces),
extrudeMeshEdges[edgeI]
);
if (patchI == -1)
{
// Determine the min zone of all connected zones.
label minZone = zoneID[eFaces[0]];
for (label i = 1; i < eFaces.size(); i++)
{
minZone = min(minZone, zoneID[eFaces[i]]);
}
zoneSidePatch[minZone]++;
}
}
}
Pstream::listCombineGather(zoneSidePatch, plusEqOp<label>());
Pstream::listCombineScatter(zoneSidePatch);
Pstream::listCombineGather(zoneZonePatch, plusEqOp<label>());
Pstream::listCombineScatter(zoneZonePatch);
}
bool lessThan(const point& x, const point& y)
{
for (direction dir = 0; dir < point::nComponents; dir++)
{
if (x[dir] < y[dir]) return true;
if (x[dir] > y[dir]) return false;
}
return false;
}
class minEqVectorListOp
{
public:
void operator()(List<vector>& x, const List<vector>& y) const
{
if (y.size())
{
if (x.size())
{
forAll(x, i)
{
if (lessThan(y[i], x[i]))
{
x[i] = y[i];
}
}
}
else
{
x = y;
}
}
}
};
// Constrain&sync normals on points that are on coupled patches.
void constrainCoupledNormals
(
const fvMesh& mesh,
const primitiveFacePatch& extrudePatch,
const labelList& regionToPoint,
vectorField& regionNormals
)
{
// Invert regionToPoint to create pointToRegions.
labelListList pointToRegions
(
invertOneToMany
(
extrudePatch.nPoints(),
regionToPoint
)
);
// Sort acc. to region so (hopefully) coupled points will do the same.
forAll(pointToRegions, pointI)
{
sort(pointToRegions[pointI]);
}
const polyBoundaryMesh& patches = mesh.boundaryMesh();
// Constrain displacement on cyclic patches
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Note: bit contentious to always do this on cyclic - should user use
// different patch type, e.g. 'cyclicSlip' instead?
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (isA<cyclicPolyPatch>(pp))
{
forAll(pp.meshPoints(), pointI)
{
Map<label>::const_iterator fnd =
extrudePatch.meshPointMap().find
(
pp.meshPoints()[pointI]
);
if (fnd != extrudePatch.meshPointMap().end())
{
// fnd() is a point on this cyclic.
const vector& cycNormal = pp.pointNormals()[pointI];
const labelList& pRegions = pointToRegions[fnd()];
forAll(pRegions, i)
{
// Remove cyclic normal component.
vector& regionNormal = regionNormals[pRegions[i]];
regionNormal -= (regionNormal&cycNormal)*cycNormal;
}
}
}
}
}
// Synchronise regionNormals
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Re-work regionNormals into multiple normals per point
List<List<vector> > pointNormals(mesh.nPoints());
forAll(pointToRegions, pointI)
{
const labelList& pRegions = pointToRegions[pointI];
label meshPointI = extrudePatch.meshPoints()[pointI];
List<vector>& pNormals = pointNormals[meshPointI];
pNormals.setSize(pRegions.size());
forAll(pRegions, i)
{
pNormals[i] = regionNormals[pRegions[i]];
}
}
// Synchronise
syncTools::syncPointList
(
mesh,
pointNormals,
minEqVectorListOp(),
List<vector>(), // nullValue
false // applySeparation
);
// Re-work back into regionNormals
forAll(pointToRegions, pointI)
{
const labelList& pRegions = pointToRegions[pointI];
label meshPointI = extrudePatch.meshPoints()[pointI];
const List<vector>& pNormals = pointNormals[meshPointI];
forAll(pRegions, i)
{
regionNormals[pRegions[i]] = pNormals[i];
}
}
}
//XXXXXXXXX
tmp<pointField> calcOffset
(
const primitiveFacePatch& extrudePatch,
......@@ -739,6 +971,16 @@ int main(int argc, char *argv[])
<< endl;
// Determine corresponding mesh edges
const labelList extrudeMeshEdges
(
extrudePatch.meshEdges
(
mesh.edges(),
mesh.pointEdges()
)
);
// Check whether the zone is internal or external faces to determine
......@@ -772,7 +1014,7 @@ int main(int argc, char *argv[])
Info<< "FaceZone " << fz.name() << " has boundary faces" << endl;
}
}
Info<< endl;
......@@ -866,36 +1108,18 @@ int main(int argc, char *argv[])
labelList zoneSidePatch(faceZones.size(), 0);
labelList zoneZonePatch(faceZones.size()*faceZones.size(), 0);
forAll(edgeFaces, edgeI)
{
const labelList& eFaces = edgeFaces[edgeI];
if (eFaces.size() == 2)
{
label zone0 = zoneID[eFaces[0]];
label zone1 = zoneID[eFaces[1]];
countExtrudePatches
(
mesh,
extrudePatch,
faceZones.size(),
zoneID,
extrudeMeshFaces,