From 234f6403ce7c369b41623b2168c20137a27d74aa Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 24 Aug 2011 13:40:10 +0100 Subject: [PATCH] Corrected delta coefficients and weighting factors on coupled patches. Henrik Rusche --- .../constraint/cyclic/cyclicFvPatch.C | 18 ++++++---- .../constraint/cyclicGgi/cyclicGgiFvPatch.C | 20 ++++++++--- .../fvPatches/constraint/ggi/ggiFvPatch.C | 16 +++++++-- .../constraint/overlapGgi/overlapGgiFvPatch.C | 16 +++++++-- .../constraint/processor/processorFvPatch.C | 35 ++++++++++++------- .../regionCouple/regionCoupleFvPatch.C | 16 +++++++-- 6 files changed, 89 insertions(+), 32 deletions(-) diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclic/cyclicFvPatch.C b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclic/cyclicFvPatch.C index 09ec2c3ba..a85ffa7c7 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclic/cyclicFvPatch.C +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclic/cyclicFvPatch.C @@ -46,7 +46,12 @@ void cyclicFvPatch::makeWeights(scalarField& w) const { const scalarField& magFa = magSf(); - scalarField deltas = nf() & fvPatch::delta(); + // Note: mag in the dot-product. + // For all valid meshes, the non-orthogonality will be less that + // 90 deg and the dot-product will be positive. For invalid + // meshes (d & s <= 0), this will stabilise the calculation + // but the result will be poor. HJ, 24/Aug/2011 + scalarField deltas = mag(nf() & fvPatch::delta()); label sizeby2 = deltas.size()/2; scalar maxMatchError = 0; @@ -100,15 +105,14 @@ void cyclicFvPatch::makeWeights(scalarField& w) const // Make patch face - neighbour cell distances void cyclicFvPatch::makeDeltaCoeffs(scalarField& dc) const { - scalarField deltas = nf() & fvPatch::delta(); - label sizeby2 = deltas.size()/2; + vectorField d = delta(); + vectorField n = nf(); + label sizeby2 = d.size()/2; for (label facei = 0; facei < sizeby2; facei++) { - scalar di = deltas[facei]; - scalar dni = deltas[facei + sizeby2]; - - dc[facei] = 1.0/(di + dni); + // Stabilised form for bad meshes. HJ, 24/Aug/2011 + dc[facei] = 1.0/mag(n[facei] & d[facei]), 0.05*mag(d[facei]); dc[facei + sizeby2] = dc[facei]; } } diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicGgi/cyclicGgiFvPatch.C b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicGgi/cyclicGgiFvPatch.C index 9dd888c32..b0766aa89 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicGgi/cyclicGgiFvPatch.C +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicGgi/cyclicGgiFvPatch.C @@ -54,10 +54,19 @@ void Foam::cyclicGgiFvPatch::makeWeights(scalarField& w) const if (cyclicGgiPolyPatch_.master()) { vectorField n = nf(); - scalarField nfc = - n & (cyclicGgiPolyPatch_.reconFaceCellCentres() - Cf()); - w = nfc/((n & (Cf() - Cn())) + nfc); + // Note: mag in the dot-product. + // For all valid meshes, the non-orthogonality will be less that + // 90 deg and the dot-product will be positive. For invalid + // meshes (d & s <= 0), this will stabilise the calculation + // but the result will be poor. HJ, 24/Aug/2011 + scalarField nfc = + mag + ( + n & (cyclicGgiPolyPatch_.reconFaceCellCentres() - Cf()) + ); + + w = nfc/(mag(n & (Cf() - Cn())) + nfc); if (bridgeOverlap()) { @@ -89,7 +98,10 @@ void Foam::cyclicGgiFvPatch::makeDeltaCoeffs(scalarField& dc) const { if (cyclicGgiPolyPatch_.master()) { - dc = 1.0/max(nf() & delta(), 0.05*mag(delta())); + // Stabilised form for bad meshes. HJ, 24/Aug/2011 + vectorField d = delta(); + + dc = 1.0/max(nf() & d, 0.05*mag(d)); if (bridgeOverlap()) { diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C b/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C index e3cbf0d89..5881daa64 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C @@ -64,9 +64,16 @@ void Foam::ggiFvPatch::makeWeights(scalarField& w) const if (ggiPolyPatch_.master()) { vectorField n = nf(); - scalarField nfc = n & (ggiPolyPatch_.reconFaceCellCentres() - Cf()); - w = nfc/((n & (Cf() - Cn())) + nfc); + // Note: mag in the dot-product. + // For all valid meshes, the non-orthogonality will be less that + // 90 deg and the dot-product will be positive. For invalid + // meshes (d & s <= 0), this will stabilise the calculation + // but the result will be poor. HJ, 24/Aug/2011 + scalarField nfc = + mag(n & (ggiPolyPatch_.reconFaceCellCentres() - Cf())); + + w = nfc/(mag(n & (Cf() - Cn())) + nfc); if (bridgeOverlap()) { @@ -100,7 +107,10 @@ void Foam::ggiFvPatch::makeDeltaCoeffs(scalarField& dc) const { if (ggiPolyPatch_.master()) { - dc = 1.0/max(nf() & delta(), 0.05*mag(delta())); + // Stabilised form for bad meshes. HJ, 24/Aug/2011 + vectorField d = delta(); + + dc = 1.0/max(nf() & d, 0.05*mag(d)); if (bridgeOverlap()) { diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/overlapGgi/overlapGgiFvPatch.C b/src/finiteVolume/fvMesh/fvPatches/constraint/overlapGgi/overlapGgiFvPatch.C index 7e42c8c7d..4220705a4 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/overlapGgi/overlapGgiFvPatch.C +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/overlapGgi/overlapGgiFvPatch.C @@ -51,9 +51,16 @@ void Foam::overlapGgiFvPatch::makeWeights(scalarField& w) const if (overlapGgiPolyPatch_.master()) { vectorField n = nf(); + + // Note: mag in the dot-product. + // For all valid meshes, the non-orthogonality will be less that + // 90 deg and the dot-product will be positive. For invalid + // meshes (d & s <= 0), this will stabilise the calculation + // but the result will be poor. HJ, 24/Aug/2011 scalarField nfc = - n & (overlapGgiPolyPatch_.reconFaceCellCentres() - Cf()); - w = nfc/((n & (Cf() - Cn())) + nfc); + mag(n & (overlapGgiPolyPatch_.reconFaceCellCentres() - Cf())); + + w = nfc/(mag(n & (Cf() - Cn())) + nfc); } else { @@ -71,7 +78,10 @@ void Foam::overlapGgiFvPatch::makeDeltaCoeffs(scalarField& dc) const { if (overlapGgiPolyPatch_.master()) { - dc = 1.0/max(nf() & delta(), 0.05*mag(delta())); + // Stabilised form for bad meshes. HJ, 24/Aug/2011 + vectorField d = delta(); + + dc = 1.0/max(nf() & d, 0.05*mag(d)); } else { diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/processor/processorFvPatch.C b/src/finiteVolume/fvMesh/fvPatches/constraint/processor/processorFvPatch.C index f3c3d029e..6a4c53e27 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/processor/processorFvPatch.C +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/processor/processorFvPatch.C @@ -46,18 +46,26 @@ void processorFvPatch::makeWeights(scalarField& w) const if (Pstream::parRun()) { // The face normals point in the opposite direction on the other side - scalarField neighbFaceCentresCn - ( - ( - procPolyPatch_.neighbFaceAreas() - /(mag(procPolyPatch_.neighbFaceAreas()) + VSMALL) - ) - & ( - procPolyPatch_.neighbFaceCentres() - - procPolyPatch_.neighbFaceCellCentres()) - ); - w = neighbFaceCentresCn/((nf() & fvPatch::delta()) + // Note: mag in the dot-product. + // For all valid meshes, the non-orthogonality will be less that + // 90 deg and the dot-product will be positive. For invalid + // meshes (d & s <= 0), this will stabilise the calculation + // but the result will be poor. HJ, 24/Aug/2011 + scalarField neighbFaceCentresCn = + mag + ( + ( + procPolyPatch_.neighbFaceAreas()/ + (mag(procPolyPatch_.neighbFaceAreas()) + VSMALL) + ) + & ( + procPolyPatch_.neighbFaceCentres() + - procPolyPatch_.neighbFaceCellCentres() + ) + ); + + w = neighbFaceCentresCn/(mag(nf() & fvPatch::delta()) + neighbFaceCentresCn); } else @@ -71,7 +79,10 @@ void processorFvPatch::makeDeltaCoeffs(scalarField& dc) const { if (Pstream::parRun()) { - dc = (1.0 - weights())/(nf() & fvPatch::delta()); + vectorField d = delta(); + + // Stabilised form for bad meshes + dc = 1.0/max((nf() & d), 0.05*mag(d)); } else { diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.C b/src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.C index bbcbc4c55..57e7df5b5 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.C +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.C @@ -53,9 +53,16 @@ void Foam::regionCoupleFvPatch::makeWeights(scalarField& w) const if (rcPolyPatch_.attached()) { vectorField n = nf(); - scalarField nfc = n & (rcPolyPatch_.reconFaceCellCentres() - Cf()); - w = nfc/((n & (Cf() - Cn())) + nfc); + // Note: mag in the dot-product. + // For all valid meshes, the non-orthogonality will be less that + // 90 deg and the dot-product will be positive. For invalid + // meshes (d & s <= 0), this will stabilise the calculation + // but the result will be poor. HJ, 24/Aug/2011 + scalarField nfc = + mag(n & (rcPolyPatch_.reconFaceCellCentres() - Cf())); + + w = nfc/(mag(n & (Cf() - Cn())) + nfc); } else { @@ -69,7 +76,10 @@ void Foam::regionCoupleFvPatch::makeDeltaCoeffs(scalarField& dc) const { if (rcPolyPatch_.attached()) { - dc = (1.0 - weights())/(nf() & fvPatch::delta()); + // Stabilised form for bad meshes. HJ, 24/Aug/2011 + vectorField d = delta(); + + dc = 1.0/max(nf() & d, 0.05*mag(d)); } else {