exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Recent cumulative port fixes related to nextRelease. Author: Vuko Vukcevic. Merge: Hrvoje Jasak.

Browse source
This commit is contained in:
Hrvoje Jasak 2018-09-05 12:18:24 +01:00
commit d7d07ea75b
10 changed files with 128 additions and 113 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

92
src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C
View file

@ -116,11 +116,11 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
forAll (owner, faceI)
{
label own = owner[faceI];
label nei = neighbour[faceI];
const label own = owner[faceI];
const label nei = neighbour[faceI];
vector d = C[nei] - C[own];
symmTensor wdd = (1.0/magSqr(d))*sqr(d);
const vector d = C[nei] - C[own];
const symmTensor wdd = (1.0/magSqr(d))*sqr(d);
dd[own] += wdd;
dd[nei] += wdd;
@ -134,21 +134,10 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
// Better version of d-vectors: Zeljko Tukovic, 25/Apr/2010
const vectorField pd = p.delta();
if (p.coupled())
forAll (pd, pFaceI)
{
forAll (pd, pFaceI)
{
const vector& d = pd[pFaceI];
dd[fc[pFaceI]] += (1.0/magSqr(d))*sqr(d);
}
}
else
{
forAll (pd, pFaceI)
{
const vector& d = pd[pFaceI];
dd[fc[pFaceI]] += (1.0/magSqr(d))*sqr(d);
}
const vector& d = pd[pFaceI];
dd[fc[pFaceI]] += (1.0/magSqr(d))*sqr(d);
}
}
@ -186,16 +175,20 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
// Least squares vectors on internal faces
forAll (owner, faceI)
{
label own = owner[faceI];
label nei = neighbour[faceI];
const label own = owner[faceI];
const label nei = neighbour[faceI];
vector d = C[nei] - C[own];
scalar magSfByMagSqrd = 1.0/magSqr(d);
const vector d = C[nei] - C[own];
const scalar magSfByMagSqrd = 1.0/magSqr(d);
lsPIn[faceI] = magSfByMagSqrd*(invDd[own] & d);
lsNIn[faceI] = -magSfByMagSqrd*(invDd[nei] & d);
}
// Get inverse dd boundary field to grab the neighbouring side eventually
const volSymmTensorField::GeometricBoundaryField& volInvDdb =
volInvDd.boundaryField();
// Least squares vectors on boundary faces
forAll (lsP.boundaryField(), patchI)
{
@ -204,13 +197,14 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
const fvPatch& p = mesh().boundary()[patchI];
const unallocLabelList& fc = p.faceCells();
const fvPatchSymmTensorField& volInvDdp = volInvDdb[patchI];
// Better version of d-vectors: Zeljko Tukovic, 25/Apr/2010
const vectorField pd = p.delta();
if (p.coupled())
if (volInvDdp.coupled())
{
const symmTensorField invDdNei =
volInvDd.boundaryField()[patchI].patchNeighbourField();
const symmTensorField invDdNei = volInvDdp.patchNeighbourField();
forAll (pd, pFaceI)
{
@ -347,8 +341,8 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
// Internal faces
forAll (owner, faceI)
{
label own = owner[faceI];
label nei = neighbour[faceI];
const label own = owner[faceI];
const label nei = neighbour[faceI];
if (uggIn[own] > SMALL)
{
@ -363,6 +357,10 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
}
}
// Get boundary field for cell volumes to fetch the other side eventually
const volScalarField::GeometricBoundaryField& cellVb =
cellV.boundaryField();
// Boundary faces
forAll (lsP.boundaryField(), patchI)
{
@ -372,16 +370,16 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
const unallocLabelList& fc = p.faceCells();
const fvsPatchScalarField& pw = w.boundaryField()[patchI];
const fvPatchScalarField& cellVp = cellVb[patchI];
const vectorField& pSf = Sf.boundaryField()[patchI];
// Get indicator for local side
const fvPatchScalarField& ugg = useGaussGrad.boundaryField()[patchI];
const scalarField pUgg = ugg.patchInternalField();
if (p.coupled())
if (cellVp.coupled())
{
const scalarField cellVInNei =
cellV.boundaryField()[patchI].patchNeighbourField();
const scalarField cellVInNei = cellVp.patchNeighbourField();
// Get indicator for neighbour side
const scalarField nUgg = ugg.patchNeighbourField();
@ -417,40 +415,6 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
}
}
// Manual check of least squares vectors
vectorField sumLsq(mesh().nCells(), vector::zero);
vectorField sumSf(mesh().nCells(), vector::zero);
const vectorField& sfIn = mesh().Sf().internalField();
// Least squares vectors on internal faces
forAll (owner, faceI)
{
sumLsq[owner[faceI]] += lsPIn[faceI];
sumLsq[neighbour[faceI]] += lsNIn[faceI];
sumSf[owner[faceI]] += sfIn[faceI];
sumSf[neighbour[faceI]] -= sfIn[faceI];
}
// Least squares vectors on boundary faces
forAll (lsP.boundaryField(), patchI)
{
const vectorField& patchLsP = lsP.boundaryField()[patchI];
const vectorField& patchSf = mesh().Sf().boundaryField()[patchI];
const unallocLabelList& fc = mesh().boundary()[patchI].faceCells();
forAll (fc, pFaceI)
{
//sumLsq[fc[pFaceI]] += 0.5*patchLsP[pFaceI]; // works!!!
sumLsq[fc[pFaceI]] += patchLsP[pFaceI];
sumSf[fc[pFaceI]] += patchSf[pFaceI];
}
}
if (debug)
{
Info<< "leastSquaresVectors::makeLeastSquaresVectors() :"

66
src/foam/interpolations/GGIInterpolation/GGIInterpolationTemplate.H
View file

@ -242,7 +242,7 @@ class GGIInterpolation
mutable scalarField* slaveFaceUncoveredFractionsPtr_;
// Private static data
// Private static data
//- Facet area error tolerance
static const debug::tolerancesSwitch areaErrorTol_;
@ -253,38 +253,35 @@ class GGIInterpolation
//- Facet bound box extension factor
static const debug::tolerancesSwitch faceBoundBoxExtendSpanFraction_;
//- Tolerance for uncovered face areas
static const debug::tolerancesSwitch uncoveredFaceAreaTol_;
//- The next 3 attributes are parameters controlling the creation
// of an octree search engine for the GGI facets neighbourhood
// determination.
//
// Tests using GGI patches of up to ~100K facets have shown that
// the following values gave the best compromise between the
// "octree creation time" vs "octree search time":
//
// octreeSearchMinNLevel_ : 3
// octreeSearchMaxLeafRatio_ : 3
// octreeSearchMaxShapeRatio_ : 1
//
// For GGI patches larger than ~100K facets, your mileage may vary.
// So these 3 control parameters are adjustable using the following
// global optimization switches:
//
// GGIOctreeSearchMinNLevel
// GGIOctreeSearchMaxLeafRatio
// GGIOctreeSearchMaxShapeRatio
//
//- The next 3 attributes are parameters controlling the creation
// of an octree search engine for the GGI facets neighbourhood
// determination.
//
// Tests using GGI patches of up to ~100K facets have shown that
// the following values gave the best compromise between the
// "octree creation time" vs "octree search time":
//
// octreeSearchMinNLevel_ : 3
// octreeSearchMaxLeafRatio_ : 3
// octreeSearchMaxShapeRatio_ : 1
//
// For GGI patches larger than ~100K facets, your mileage may vary.
// So these 3 control parameters are adjustable using the following
// global optimization switches:
//
// GGIOctreeSearchMinNLevel
// GGIOctreeSearchMaxLeafRatio
// GGIOctreeSearchMaxShapeRatio
//
//- Octree search: minNlevel parameter for octree constructor
static const debug::optimisationSwitch octreeSearchMinNLevel_;
//- Octree search: minNlevel parameter for octree constructor
static const debug::optimisationSwitch octreeSearchMinNLevel_;
//- Octree search: maxLeafRatio parameter for octree constructor
static const debug::optimisationSwitch octreeSearchMaxLeafRatio_;
//- Octree search: maxLeafRatio parameter for octree constructor
static const debug::optimisationSwitch octreeSearchMaxLeafRatio_;
//- Octree search: maxShapeRatio parameter for octree constructor
static const debug::optimisationSwitch octreeSearchMaxShapeRatio_;
//- Octree search: maxShapeRatio parameter for octree constructor
static const debug::optimisationSwitch octreeSearchMaxShapeRatio_;
// Private Member Functions
@ -307,6 +304,7 @@ class GGIInterpolation
//- Calculate point weights
void calcSlavePointWeights() const;
// Helper functions for parallel search
//- Is parallel cutting active?
@ -582,6 +580,14 @@ public:
~GGIInterpolation();
// Public static data
//- Tolerance for uncovered face areas. Public since the
// MixingPlaneInterpolation needs to be able to change it. Seems more
// appropriate than having friends.
static debug::tolerancesSwitch uncoveredFaceAreaTol_;
// Member Functions
// Access

2
src/foam/interpolations/GGIInterpolation/GGIInterpolationWeights.C
View file

@ -69,7 +69,7 @@ GGIInterpolation<MasterPatch, SlavePatch>::featureCosTol_
template<class MasterPatch, class SlavePatch>
const Foam::debug::tolerancesSwitch
Foam::debug::tolerancesSwitch
GGIInterpolation<MasterPatch, SlavePatch>::uncoveredFaceAreaTol_
(
"GGIUncoveredFaceAreaTol",

34
src/foam/interpolations/MixingPlaneInterpolation/MixingPlaneInterpolationAddressing.C
View file

@ -125,6 +125,30 @@ void MixingPlaneInterpolation<MasterPatch, SlavePatch>::calcAddressing() const
// introduced by the possibly different
// master and slave patch mesh resolution.
// Note: The weights in the GGI interpolation must be scaled when using
// interpolation from/to ribbon patches. Now, the GGI interpolation can also
// take care of partially overlapping faces where we must not scale the
// weights (e.g. bridgeOverlap option for partial symmetry plane
// treatment in uncovered faces). In order to avoid having partially
// overlapping faces in the mixing plane, we will manually set the
// corresponding tolerance here and reset it after we calculate the
// addressing. VV, 31/Aug/2018.
// Save old tolerance
const scalar oldTolerance =
GGIInterpolation
<
standAlonePatch,
standAlonePatch
>::uncoveredFaceAreaTol_();
// Set tolerance to zero such that no faces qualify as partially
// overlapping faces and all the weights are properly scaled
GGIInterpolation
<
standAlonePatch,
standAlonePatch
>::uncoveredFaceAreaTol_ = 0.0;
GGIInterpolation<standAlonePatch, standAlonePatch>
masterCircumAvgPatchToPatch
(
@ -202,6 +226,16 @@ void MixingPlaneInterpolation<MasterPatch, SlavePatch>::calcAddressing() const
slaveProfileToPatchWeightsPtr_ =
new scalarListList(slaveCircumAvgPatchToPatch.slaveWeights());
// Now that we have finished calculating the addressing and the weights,
// let's reset the old tolerances in case some other GGIInterpolation with
// stand alone patches uses the partially covered face machinery.
GGIInterpolation
<
standAlonePatch,
standAlonePatch
>::uncoveredFaceAreaTol_ = oldTolerance;
}

12
src/foam/interpolations/MixingPlaneInterpolation/MixingPlaneInterpolationTemplate.H
View file

@ -53,11 +53,11 @@ Contributor
Giving credit where credit is due:
1: Hakan Nilsson from the Chalmers University of Technology came up with
the initial idea of using a combination of 2 GGI interfaces sharing a
common single-face 360 degree ribbons patch in order to compute the
circumferential average of fields. This implementation of the
mixingPlane interpolation algorithm is an exploration of this simple
but rather powerful idea.
the initial idea of using a combination of 2 GGI interfaces sharing a
common single-face 360 degree ribbons patch in order to compute the
circumferential average of fields. This implementation of the
mixingPlane interpolation algorithm is an exploration of this simple
but rather powerful idea.
2: Maryse Page from Hydro-Quebec provided many test cases and many
simulation runs for testing this interpolation algorithm. Testing is
@ -240,7 +240,7 @@ class MixingPlaneInterpolation
//- Interpolation profile
// This list of points defines the profile generating the 'n'
// circular bands we are going to use for
// the circumferential averaging algorithm. For for 'n'
// the circumferential averaging algorithm. For 'n'
// circular bands, we need 'n + 1' points.
mutable pointField interpolationProfile_;

4
src/foam/interpolations/MixingPlaneInterpolation/MixingPlaneProfile.C
View file

@ -313,13 +313,13 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::calcProfile() const
if (debug > 1)
{
// Write histograms
forAllIter (profileHistogram, masterHistogram, zHi)
forAllConstIter (profileHistogram, masterHistogram, zHi)
{
Info<< "master histo (z, n): (" << zHi->first << " "
<< static_cast<int>(zHi->second.size()) << ")" << endl;
}
forAllIter (profileHistogram, slaveHistogram, zHi)
forAllConstIter (profileHistogram, slaveHistogram, zHi)
{
Info<< "slave histo (z, n): (" << zHi->first << " "
<< static_cast<int>(zHi->second.size()) << ")" << endl;

3
src/foam/meshes/polyMesh/polyMesh.C
View file

@ -1232,10 +1232,11 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
geometricD_ = Vector<label>::zero;
solutionD_ = Vector<label>::zero;
// Update all function objects
// Update all mesh objects
// Moved from fvMesh.C in 1.6.x merge. HJ, 29/Aug/2010
meshObjectBase::allMovePoints<polyMesh>(*this);
// Update all function objects
const_cast<Time&>(time()).functionObjects().movePoints(allPoints_);
return sweptVols;

5
src/foam/meshes/polyMesh/polyMeshUpdate.C
View file

@ -75,9 +75,12 @@ void Foam::polyMesh::updateMesh(const mapPolyMesh& mpm)
geometricD_ = Vector<label>::zero;
solutionD_ = Vector<label>::zero;
// Update all function objects
// Update all mesh objects
// Moved from fvMesh.C in 1.6.x merge. HJ, 29/Aug/2010
meshObjectBase::allUpdateTopology<polyMesh>(*this, mpm);
// Update all function objects
const_cast<Time&>(time()).functionObjects().updateMesh(mpm);
}

18
src/overset/oversetMesh/include/writeOversetMasks.H
View file

@ -18,7 +18,8 @@
IOobject::NO_WRITE
),
mesh,
dimless
dimless,
"zeroGradient"
);
scalarField& regionIndexIn = regionIndex.internalField();
@ -32,13 +33,18 @@
}
}
// Update boundary values
forAll (regionIndex.boundaryField(), patchI)
// Update boundary values, making sure that we skip the overset patch
volScalarField::GeometricBoundaryField& regionIndexb =
regionIndex.boundaryField();
forAll(regionIndexb, patchI)
{
if (!regionIndex.boundaryField()[patchI].coupled())
// Get the patch field
fvPatchScalarField& ripf = regionIndexb[patchI];
if (!isA<oversetFvPatchScalarField>(ripf))
{
regionIndex.boundaryField()[patchI] =
regionIndex.boundaryField()[patchI].patchInternalField();
ripf = ripf.patchInternalField();
}
}

5
tutorials/overset/parallelTests/NACA4412-ThreeLayerOverlapFringeAssemby-100Processors/Allrun
View file

@ -44,5 +44,6 @@ transformPoints -scale "(0.901 1 0.905512)" > log.transformPoints
runApplication decomposePar -cellDist
runApplication decomposeSets -writeEmptySets
runParallel calcOverset 100
runParallel oversetRegionIndex 100
# Run with -oversubscribe flag to circumvent mpi config issues on small machines
mpirun -oversubscribe -np 100 calcOverset -parallel > log.calcOverset