From 53bd62e50da02d7c77a3d5a78c103e05eb1f97aa Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Wed, 30 May 2012 15:25:17 +0100
Subject: [PATCH] LimitedSchemes: Improved handling of division by very small
 number

---
 .../limitedSchemes/LimitedScheme/NVDTVD.H     | 26 ++++++++++++-------
 .../limitedSchemes/LimitedScheme/NVDVTVDV.H   | 26 ++++++++++++-------
 2 files changed, 34 insertions(+), 18 deletions(-)

diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/NVDTVD.H b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/NVDTVD.H
index f4ef12ff64e..688c26cd777 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/NVDTVD.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/NVDTVD.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,10 +81,14 @@ public:
                 gradcf = d & gradcN;
             }
 
-            // Stabilise for division
-            gradcf = stabilise(gradcf, VSMALL);
-
-            return 1 - 0.5*gradf/gradcf;
+            if (mag(gradf) >= 1000*mag(gradcf))
+            {
+                return 1 - 0.5*1000*sign(gradcf)*sign(gradf);
+            }
+            else
+            {
+                return 1 - 0.5*gradf/gradcf;
+            }
         }
 
 
@@ -111,10 +115,14 @@ public:
                 gradcf = d & gradcN;
             }
 
-            // Stabilise for division
-            gradf = stabilise(gradf, VSMALL);
-
-            return 2*(gradcf/gradf) - 1;
+            if (mag(gradcf) >= 1000*mag(gradf))
+            {
+                return 2*1000*sign(gradcf)*sign(gradf) - 1;
+            }
+            else
+            {
+                return 2*(gradcf/gradf) - 1;
+            }
         }
 };
 
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/NVDVTVDV.H b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/NVDVTVDV.H
index df6c110a659..ac15c38976f 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/NVDVTVDV.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/NVDVTVDV.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,10 +81,14 @@ public:
                 gradcf = gradfV & (d & gradcN);
             }
 
-            // Stabilise for division
-            gradcf = stabilise(gradcf, VSMALL);
-
-            return 1 - 0.5*gradf/gradcf;
+            if (mag(gradf) >= 1000*mag(gradcf))
+            {
+                return 1 - 0.5*1000*sign(gradcf)*sign(gradf);
+            }
+            else
+            {
+                return 1 - 0.5*gradf/gradcf;
+            }
         }
 
 
@@ -112,10 +116,14 @@ public:
                 gradcf = gradfV & (d & gradcN);
             }
 
-            // Stabilise for division
-            gradf = stabilise(gradf, VSMALL);
-
-            return 2*(gradcf/gradf) - 1;
+            if (mag(gradcf) >= 1000*mag(gradf))
+            {
+                return 2*1000*sign(gradcf)*sign(gradf) - 1;
+            }
+            else
+            {
+                return 2*(gradcf/gradf) - 1;
+            }
         }
 };
 
-- 
GitLab