From 80bb6e2b0f03c8ae09c2fe6ca2192f282e1f0adf Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Wed, 13 Nov 2019 14:45:58 +0000 Subject: [PATCH] blockMesh::projectFace: Improved robustness of convergence check --- .../projectCurveEdge/projectCurveEdge.C | 16 +++--- .../blockEdges/projectEdge/projectEdge.C | 26 +++++---- .../blockFaces/projectFace/projectFace.C | 55 ++++++++----------- 3 files changed, 48 insertions(+), 49 deletions(-) diff --git a/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.C b/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.C index 69f906c045..be58bcb71b 100644 --- a/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.C +++ b/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2016-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2016-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -147,11 +147,12 @@ Foam::projectCurveEdge::position(const scalarList& lambdas) const // Upper limit for number of iterations const label maxIter = 10; + // Residual tolerance const scalar relTol = 0.1; const scalar absTol = 1e-4; - scalar initialResidual = 0.0; + scalar initialResidual = 0; for (label iter = 0; iter < maxIter; iter++) { @@ -174,7 +175,7 @@ Foam::projectCurveEdge::position(const scalarList& lambdas) const { points[0] = startPt; } - if (lambdas.last() > 1.0-small) + if (lambdas.last() > 1 - small) { points.last() = endPt; } @@ -191,10 +192,11 @@ Foam::projectCurveEdge::position(const scalarList& lambdas) const // Calculate lambdas (normalised coordinate along edge) scalarField projLambdas(points.size()); { - projLambdas[0] = 0.0; + projLambdas[0] = 0; for (label i = 1; i < points.size(); i++) { - projLambdas[i] = projLambdas[i-1] + mag(points[i]-points[i-1]); + projLambdas[i] = + projLambdas[i-1] + mag(points[i] - points[i-1]); } projLambdas /= projLambdas.last(); } @@ -214,10 +216,10 @@ Foam::projectCurveEdge::position(const scalarList& lambdas) const { predicted += weights[indexi]*points[indices[indexi]]; } - residual[i] = predicted-points[i]; + residual[i] = predicted - points[i]; } - scalar scalarResidual = sum(mag(residual)); + const scalar scalarResidual = sum(mag(residual)); if (debug) { diff --git a/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C b/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C index 0b6e112e77..8f8bd86ecf 100644 --- a/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C +++ b/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2016-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2016-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -51,7 +51,7 @@ void Foam::projectEdge::findNearest { if (surfaces_.size()) { - const scalar distSqr = magSqr(points_[end_]-points_[start_]); + const scalar distSqr = magSqr(points_[end_] - points_[start_]); pointField boundaryNear(1); List boundaryConstraint(1); @@ -111,11 +111,14 @@ Foam::projectEdge::projectEdge Foam::point Foam::projectEdge::position(const scalar lambda) const { // Initial guess - const point start(points_[start_] + lambda*(points_[end_]-points_[start_])); + const point start + ( + points_[start_] + lambda*(points_[end_] - points_[start_]) + ); point near(start); - if (lambda >= small && lambda < 1.0-small) + if (lambda >= small && lambda < 1 - small) { pointConstraint constraint; findNearest(start, near, constraint); @@ -159,11 +162,12 @@ Foam::projectEdge::position(const scalarList& lambdas) const // Upper limit for number of iterations const label maxIter = 10; + // Residual tolerance const scalar relTol = 0.1; const scalar absTol = 1e-4; - scalar initialResidual = 0.0; + scalar initialResidual = 0; for (label iter = 0; iter < maxIter; iter++) { @@ -186,7 +190,7 @@ Foam::projectEdge::position(const scalarList& lambdas) const { points[0] = startPt; } - if (lambdas.last() > 1.0-small) + if (lambdas.last() > 1 - small) { points.last() = endPt; } @@ -203,10 +207,11 @@ Foam::projectEdge::position(const scalarList& lambdas) const // Calculate lambdas (normalised coordinate along edge) scalarField projLambdas(points.size()); { - projLambdas[0] = 0.0; + projLambdas[0] = 0; for (label i = 1; i < points.size(); i++) { - projLambdas[i] = projLambdas[i-1] + mag(points[i]-points[i-1]); + projLambdas[i] = + projLambdas[i-1] + mag(points[i] - points[i-1]); } projLambdas /= projLambdas.last(); } @@ -226,10 +231,10 @@ Foam::projectEdge::position(const scalarList& lambdas) const { predicted += weights[indexi]*points[indices[indexi]]; } - residual[i] = predicted-points[i]; + residual[i] = predicted - points[i]; } - scalar scalarResidual = sum(mag(residual)); + const scalar scalarResidual = sum(mag(residual)); if (debug) { @@ -250,7 +255,6 @@ Foam::projectEdge::position(const scalarList& lambdas) const break; } - if (debugStr.valid()) { forAll(points, i) diff --git a/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C index f692466722..5fa544600c 100644 --- a/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C +++ b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C @@ -96,10 +96,10 @@ void Foam::blockFaces::projectFace::calcLambdas { const label ij = index(n, labelPair(i, j)); const label iMin1j = index(n, labelPair(i-1, j)); - lambdaI[ij] = lambdaI[iMin1j] + mag(points[ij]-points[iMin1j]); + lambdaI[ij] = lambdaI[iMin1j] + mag(points[ij] - points[iMin1j]); const label ijMin1 = index(n, labelPair(i, j-1)); - lambdaJ[ij] = lambdaJ[ijMin1] + mag(points[ij]-points[ijMin1]); + lambdaJ[ij] = lambdaJ[ijMin1] + mag(points[ij] - points[ijMin1]); } } @@ -206,10 +206,9 @@ void Foam::blockFaces::projectFace::project // Residual tolerance const scalar relTol = 0.1; + const scalar absTol = 1e-4; scalar initialResidual = 0; - scalar iResidual = 0; - scalar jResidual = 0; for (label iter = 0; iter < maxIter; iter++) { @@ -237,26 +236,6 @@ void Foam::blockFaces::projectFace::project } } - if (debug) - { - Pout<< "Iter:" << iter << " initialResidual:" << initialResidual - << " iResidual + jResidual:" << iResidual + jResidual << endl; - } - - - if - ( - iter > 0 - && ( - initialResidual < small - || (iResidual + jResidual)/initialResidual < relTol - ) - ) - { - break; - } - - // Predict along i vectorField residual(points.size(), vector::zero); @@ -274,7 +253,7 @@ void Foam::blockFaces::projectFace::project const label iMin1j = index(n, labelPair(i-1, j)); projLambdas[i] = projLambdas[i-1] - +mag(points[ij]-points[iMin1j]); + + mag(points[ij] - points[iMin1j]); } projLambdas /= projLambdas.last(); @@ -293,7 +272,7 @@ void Foam::blockFaces::projectFace::project index(n, labelPair(indices[indexi], j)); predicted += weights[indexi]*points[ptIndex]; } - residual[ij] = predicted-points[ij]; + residual[ij] = predicted - points[ij]; } } @@ -306,7 +285,7 @@ void Foam::blockFaces::projectFace::project } } - iResidual = sum(mag(residual)); + scalar scalarResidual = sum(mag(residual)); // Update points before doing j. Note: is this needed? Complicates // residual checking. @@ -346,7 +325,7 @@ void Foam::blockFaces::projectFace::project index(n, labelPair(i, indices[indexi])); predicted += weights[indexi]*points[ptIndex]; } - residual[ij] = predicted-points[ij]; + residual[ij] = predicted - points[ij]; } } @@ -359,11 +338,25 @@ void Foam::blockFaces::projectFace::project } } - jResidual = sum(mag(residual)); + scalarResidual += sum(mag(residual)); - if (iter == 0) + if (debug) { - initialResidual = iResidual + jResidual; + Pout<< "Iter:" << iter << " initialResidual:" << initialResidual + << " scalarResidual:" << scalarResidual << endl; + } + + if (scalarResidual < absTol*lambdaI.size()) + { + break; + } + else if (iter == 0) + { + initialResidual = scalarResidual; + } + else if (scalarResidual/initialResidual < relTol) + { + break; } points += residual;