From 43080d70428b0973c64794591921ff3100e21670 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Thu, 13 Jun 2019 11:07:18 +0200 Subject: [PATCH 1/4] Fixed typos in comments --- .../decomposeTools/finiteVolume/decomposeMesh.C | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/decomposeMesh.C b/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/decomposeMesh.C index 57fe4ec5b..db73122ac 100644 --- a/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/decomposeMesh.C +++ b/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/decomposeMesh.C @@ -459,7 +459,7 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches) // cell with multiple faces, we need to collect neighbour cell // indices on the other side and possibly swap the order of // adding faces. The same "swapping" of insertion order needs to - // happend on the slave side, but now we sort on the remote + // happen on the slave side, but now we sort on the remote // (master) data and not on the local (slave) data as we do for // master processor. VV, 16/Feb/2019. @@ -588,7 +588,7 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches) // above, owner and neighbour proc are the same. In // order to avoid using cellToProc list for slave // processor (where ownCellI is actually found on the - // other side), use curNbrPtc with patchFace indes + // other side), use curNbrPtc with patchFace index const label& ownerProc = curNbrPtc[patchFaceI]; // Add the patch face into the list in the correct From ce6f4f6473aeeb2b2e8ab800678fa4042d802aef Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Thu, 13 Jun 2019 11:10:32 +0200 Subject: [PATCH 2/4] Update in fvcReconstruct to enable hinv Reconstruction equation reformulated in a way such that the inverse is always calculated for a dimensionless dyadic tensor (Sf*Sf/magSf^2) instead of one that scales with magSf (Sf*Sf/magSf). This proved to be problematic when we had extremely small cells which triggered the calculation of inverse using Householder's method. --- .../finiteVolume/fvc/fvcReconstruct.C | 76 +++++++++++++------ 1 file changed, 54 insertions(+), 22 deletions(-) diff --git a/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C b/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C index ed2e8fcab..06982ba1e 100644 --- a/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C +++ b/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C @@ -76,34 +76,66 @@ reconstruct GeometricField& reconField = treconField(); - // Note: - // 1) Reconstruction is only available in cell centres: there is no need + // Notes regarding boundary: + // 1. Reconstruction is only available in cell centres: there is no need // to invert the tensor on the boundary - // 2) For boundaries, the only reconstructed data is the flux times - // normal. Based on this guess, boundary conditions can adjust - // patch values - // HJ, 12/Aug/2011 + // 2. For boundaries, the only reconstructed data is the flux times + // normal. Based on this guess, boundary conditions can adjust + // patch values. HJ, 12/Aug/2011 - GeometricField fluxTimesNormal = - surfaceSum((mesh.Sf()/mesh.magSf())*ssf); + // Notes regarding procedure: + // 1. hinv inverse must be used to stabilise the inverse on bad meshes + // but it gives strange failures because it unnecessarily avoids + // performing ordinary inverse for meshes with reasonably sized + // determinant (e.g. if SfSf/magSf is small). HJ, 19/Aug/2015 + // 2. hinv has been stabilised now. HJ, 22/Mar/2019 + // 3. But we still need to make sure that the determinant is not extremely + // small, which may happen for extremely small meshes. We avoid this + // issue by dividing the reconstruction equation with magSf^2 (instead of + // magSf), which basically makes the dyadic tensor that we need to invert + // dimensionless. VV, 13/Jun/2019 - // Note: hinv inverse must be used to stabilise the inverse on bad meshes - // but it gives strange failures. Fixed hinv. HJ, 22/Mar/2019 - // HJ, 19/Aug/2015 - // Temporarily reverting hinv: Vuko Vukcevic, 6/Jun/2019 + // Here's a short derivation in a Latex--like notation, where: + // - Sf is the surface area vector + // - uf is the face velocity factor (or field to be reconstructed) + // - Ff is the face flux + // - magSf is the magnitude of the surface area vector + // - uP is the velocity field in the cell centre + // - G is the surface area dyadic tensor + // - Fn is the vector representing surface sum of directional fluxes + // - \dprod is a dot product + + // 1. Sf \dprod uf = Ff + // Multiply Eq (1) with Sf/magSf^2 + // 2. \frac{Sf Sf}{magSf^2} \dprod uf = \frac{Sf Ff}{magSf^2} + // Sum Eq (2) over all the faces + // 3. \sum_f(\frac{Sf Sf}{magSf^2} \dprod uf) = \sum_f(\frac{Sf Ff}{magSf^2}) + // Assume first order extrapolation of uf, e.g. uP = uf + // 4. \sum_f(\frac{Sf Sf}{magSf^2}) \dprod uP) = \sum_f(\frac{Sf Ff}{magSf^2}) + // Use shorthand notation + // 5. G \dprod uP = Fn + // 6. uP = G^-1 \dprod Fn + + // Fn -> fluxTimesNormal + // G -> G + + // Calculate sum of the directional fluxes + const surfaceScalarField magSfSqr = sqr(mesh.magSf()); + const GeometricField fluxTimesNormal = + surfaceSum((mesh.Sf()/magSfSqr)*ssf); + + // Calculate the G tensor + const volSymmTensorField G = surfaceSum(sqr(mesh.Sf())/magSfSqr); + + // Finally calculate the reconstructed field using hinv for stabilisation on + // really bad fvMesh bits (uses ordinary inverse most of the time, see + // tensor.C) reconField.internalField() = - ( - // hinv - inv - ( - surfaceSum(sqr(mesh.Sf())/mesh.magSf())().internalField() - ) - & fluxTimesNormal.internalField() - ); + hinv(G.internalField()) & fluxTimesNormal.internalField(); + // Boundary value update reconField.boundaryField() = fluxTimesNormal.boundaryField(); - - treconField().correctBoundaryConditions(); + reconField.correctBoundaryConditions(); return treconField; } From 07c16f14545bfd8fec77dc8b7360358feb08f4a4 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Thu, 13 Jun 2019 11:14:49 +0200 Subject: [PATCH 3/4] Updates to oversetMesh include files --- .../oversetMesh/include/oversetAlphaCourantNo.H | 17 ++++++++++++----- .../oversetMesh/include/oversetCourantNo.H | 6 ++++-- 2 files changed, 16 insertions(+), 7 deletions(-) diff --git a/src/overset/oversetMesh/include/oversetAlphaCourantNo.H b/src/overset/oversetMesh/include/oversetAlphaCourantNo.H index aefb835fe..dd04e33dd 100644 --- a/src/overset/oversetMesh/include/oversetAlphaCourantNo.H +++ b/src/overset/oversetMesh/include/oversetAlphaCourantNo.H @@ -38,21 +38,28 @@ scalar maxAlphaCo scalar alphaCoNum = 0.0; scalar meanAlphaCoNum = 0.0; -surfaceScalarField alpha1f = +const surfaceScalarField alpha1f = fvc::interpolate(min(max(alpha1, scalar(0)), scalar(1))); -const dimensionedScalar alphaOffset("alphaOffset", dimless, dAlpha); +const dimensionedScalar alphaOffset +( + "alphaOffset", + dimless, + runTime.controlDict().lookupOrDefault("dAlpha", 0.01) +); if (mesh.nInternalFaces()) { - surfaceScalarField magAlphaPhi + const oversetMesh& om = oversetMesh::New(mesh); + + const surfaceScalarField magAlphaPhi ( pos(alpha1f - alphaOffset)* pos(scalar(1) - alphaOffset - alpha1f)* - mag(faceOversetMask*phi) + mag(om.sGamma()*phi) ); - surfaceScalarField SfUfbyDelta = + const surfaceScalarField SfUfbyDelta = mesh.surfaceInterpolation::deltaCoeffs()*magAlphaPhi; const scalar deltaT = runTime.deltaT().value(); diff --git a/src/overset/oversetMesh/include/oversetCourantNo.H b/src/overset/oversetMesh/include/oversetCourantNo.H index 261ef8b77..008b0a727 100644 --- a/src/overset/oversetMesh/include/oversetCourantNo.H +++ b/src/overset/oversetMesh/include/oversetCourantNo.H @@ -36,9 +36,11 @@ scalar velMag = 0.0; if (mesh.nInternalFaces()) { - surfaceScalarField magPhi = mag(faceOversetMask*phi); + const oversetMesh& om = oversetMesh::New(mesh); - surfaceScalarField SfUfbyDelta = + const surfaceScalarField magPhi = mag(om.sGamma()*phi); + + const surfaceScalarField SfUfbyDelta = mesh.surfaceInterpolation::deltaCoeffs()*magPhi; CoNum = max(SfUfbyDelta/mesh.magSf()) From a5a326a83ba4525fcde6ad4e1ad54765bf11c7f8 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Thu, 13 Jun 2019 15:14:49 +0200 Subject: [PATCH 4/4] Update to extrudeMesh.C Use default mergeTolerance of 1e-4 to have backward compatibility with the tutorials. --- .../utilities/mesh/generation/extrudeMesh/extrudeMesh.C | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeMesh.C b/applications/utilities/mesh/generation/extrudeMesh/extrudeMesh.C index 42aa2f5b5..d2170d0a9 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeMesh.C +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeMesh.C @@ -25,7 +25,8 @@ Description Extrude mesh from existing patch (by default outwards facing normals; optional flips faces) or from patch read from file. - Note: Merges close points so be careful. + Note: Merges close points so be careful. This can be controlled with the + optional mergeTolerance parameter (1e-4) by default. Type of extrusion prescribed by run-time selectable model. @@ -181,7 +182,8 @@ int main(int argc, char *argv[]) const vector span = bb.span(); // Read merge tolerance - const scalar mergeTolerance = readScalar(dict.lookup("mergeTolerance")); + const scalar mergeTolerance = + dict.lookupOrDefault("mergeTolerance", 1e-4); const scalar mergeDim = mergeTolerance*bb.minDim(); Info<< "Mesh bounding box : " << bb << nl