LimitedSchemes: Improved handling of division by very small number

This commit is contained in:
Henry
2012-05-30 15:25:17 +01:00
parent 3e03e3ba07
commit 53bd62e50d
2 changed files with 34 additions and 18 deletions

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -81,11 +81,15 @@ public:
gradcf = d & gradcN; gradcf = d & gradcN;
} }
// Stabilise for division if (mag(gradf) >= 1000*mag(gradcf))
gradcf = stabilise(gradcf, VSMALL); {
return 1 - 0.5*1000*sign(gradcf)*sign(gradf);
}
else
{
return 1 - 0.5*gradf/gradcf; return 1 - 0.5*gradf/gradcf;
} }
}
scalar r scalar r
@ -111,11 +115,15 @@ public:
gradcf = d & gradcN; gradcf = d & gradcN;
} }
// Stabilise for division if (mag(gradcf) >= 1000*mag(gradf))
gradf = stabilise(gradf, VSMALL); {
return 2*1000*sign(gradcf)*sign(gradf) - 1;
}
else
{
return 2*(gradcf/gradf) - 1; return 2*(gradcf/gradf) - 1;
} }
}
}; };

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -81,11 +81,15 @@ public:
gradcf = gradfV & (d & gradcN); gradcf = gradfV & (d & gradcN);
} }
// Stabilise for division if (mag(gradf) >= 1000*mag(gradcf))
gradcf = stabilise(gradcf, VSMALL); {
return 1 - 0.5*1000*sign(gradcf)*sign(gradf);
}
else
{
return 1 - 0.5*gradf/gradcf; return 1 - 0.5*gradf/gradcf;
} }
}
scalar r scalar r
@ -112,11 +116,15 @@ public:
gradcf = gradfV & (d & gradcN); gradcf = gradfV & (d & gradcN);
} }
// Stabilise for division if (mag(gradcf) >= 1000*mag(gradf))
gradf = stabilise(gradf, VSMALL); {
return 2*1000*sign(gradcf)*sign(gradf) - 1;
}
else
{
return 2*(gradcf/gradf) - 1; return 2*(gradcf/gradf) - 1;
} }
}
}; };