epfl-ecps / channelflow

Channelflow is a software system for numerical analysis of the incompressible fluid flow in channel geometries, written in C++ and MPI-parallelized.
GNU General Public License v2.0
67 stars 29 forks source link

LaminarTest fails for small Vsuck #20

Closed sajjadazimi closed 4 years ago

sajjadazimi commented 5 years ago

Describe the bug The test `laminarTest' implemented in #9 fails for small values of suction velocity. The test is to calculate and check the laminar profile for 1000 randomly generated physical variables. When suction is small and dPdx is not zero the test fails. For example for Vs = 1e-7 and dPdx=0.01 the error is higher than the defined threshold and the test fails. The reason seems to be related to multiplication and division of the big and small numbers in the calculation of the laminar profile. As a result of the failed test, CodeCov shows wrong test coverage of the code in pull requests.

Expected Result

Passing the laminarTest for every value of Vs.

Actual Result

For small Vs with nonzero dPdx, laminarTest fails.

Steps to reproduce the issue Two ways: 1 - Run the test for enough times until it fails (~ 10 times)

2 - In the code of the test set manually:

Nsteps = 2;
dPdx = 0.01;
Vsuck = 1e-7;
Ubulk = 1;
a=-1;
b=1;
ua=-1;
ub=1;
Reynolds =1000;
nu=0.001;
johnfgibson commented 5 years ago

Interesting. I had expected the results to be single-precision accurate for |Vsuck H/nu| \approx 1e-8, near the crossover between the Vsuck == 0 and the |Vsuck H / nu| nonzero cases. The test maxerr = 1e-10 is too tight for that. I probably should rewrite the laminarProfile code for three cases instead of two:

Treat Vsuck == 0 with the simple quadratic formula. Treat 0 < |Vsuck H / nu| < eps with a Taylor expansion of the Vsuck nonzero formula Treat eps < |Vsuck H / nu| with the the Vsuck nonzero formula

and figure out, given this, what the best value for eps is to keep the test error within a few orders of magnitude of machine epsilon.

johnfgibson commented 4 years ago

Ok, I've finally worked this one out and should have it cleaned up for a PR relatively soon. The code is working for VsuckH/nu -> 0 but needs cleaning up and documentation.

It turns out the straightforward exponential formulae start to lose digits of precision for Vsuck H / nu ~ 0.1 (even using the C expm1(x) function instead of exp(x) -1) . So I did 6th-order Taylor expansions for the pressure-ASBL formula and the dP/dx / Ubulk relation and set the crossover at VsuckH/nu = 0.01. That keeps the error roughly 1e-12 or better for Vsuck H/nu ranging from 10 to 1e-14. I modified laminarText.cpp to do more thorough checking across orders of magnitude of bulk-flow constants, scaling dPdx H^2/nu and Ubulk over four orders of magnitude, and Vsuck H/nu from 1e-12 to 10. It all checks out.