zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

blockMesh::projectFace: Improved robustness of convergence check

Browse Source
This commit is contained in:
Henry Weller
2019-11-13 14:45:58 +00:00
parent e7128a0852
commit 80bb6e2b0f
3 changed files with 48 additions and 49 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.C
View File

@ -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)
{

26
src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C
View File

@ -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<pointConstraint> 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)

55
src/mesh/blockMesh/blockFaces/projectFace/projectFace.C
View File

@ -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;