fix numerical overflow in pressure solve

[johnfgibson]

This PR fixes issue https://github.com/epfl-ecps/channelflow/issues/4 , a numerical overflow issue for a boundary value problem in PressureSolver. It only gets triggered for y domains [a,b] with large a or b and very fine discretizations (large values of kx/Lx or kz/Lz).

The prior formulation of the general solution of the BVP was g(y) = c exp(mu y) + d exp(-mu y), where mu = 2 pi sqrt( (kx/Lx)^2 + (kz/Lz)^2). The resultant formulae for c and d in terms of the boundary conditions have terms with exp(mu a) and exp(mu b). When I ran an ASBL simulation with [a,b] = [0,4], Lx=4, Lz=2, Nx = 128, Nz=162, this produced kxmax=41, kzmax=53 and and mu = 180, approximately. Then exp(mu b) = exp(720) = 5e312, which is overflow in double precision!

The solution is to reformulate the general solution as g(y) = c exp(mu(y-a)) + d exp(-mu(y-b)). Note that exp(mu(y-a)) and exp(-mu(y-b)) are bounded between 0 and 1 on the interval y in [a,b] for mu>0, and thus the formulae for the constants c and d are much better behaved --they involve exponentials only of the form exp(-mu H) where H = b-a, and it's fine if those underflow to zero.

The reformulation passes pressureTest and my triggering ASBL case now computes fine.

[johnfgibson]

Is there a problem travis builds? The revision compiles fine and all tests pass on my local machine.

[sajjadazimi]

style format test failed

[johnfgibson]

Aha. Can you remind me how I can run a style test locally?

EDIT: nevermind, I see the share/qa/style_test.sh file. now to install llvm-6.0....

[johnfgibson]

One LLVM 6.0.1 build and a few clang-format iterations later, I believe this PR should pass the style test. I got channelflow/poissonsolver.cpp: 0 changes needed to conform to style guide from share/qa/style_test.sh.

However I think there are some overly-restrictive style requirements that interfere with code readability. For example I had to change the more readable

Complex alpha = nu_ * vkyy.eval_a() - pky.eval_a();
Complex beta  = nu_ * vkyy.eval_b() - pky.eval_b();

to the less readable

Complex alpha = nu_ * vkyy.eval_a() - pky.eval_a();
Complex beta = nu_ * vkyy.eval_b() - pky.eval_b();

to satisfy an exactly-one-space-between symbols requirement.

[codecov-io]

Codecov Report

Merging #5 into master will increase coverage by <.01%. The diff coverage is 100%.

@@            Coverage Diff             @@
##           master       #5      +/-   ##
==========================================
+ Coverage    27.1%   27.11%   +<.01%     
==========================================
  Files          63       63              
  Lines       15320    15321       +1     
  Branches     7256     7256              
==========================================
+ Hits         4153     4154       +1     
  Misses      11167    11167

Impacted Files	Coverage Δ
channelflow/poissonsolver.cpp	65.23% <100%> (+0.13%)	:arrow_up:

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update f1488d3...fd88a21. Read the comment docs.

[johnfgibson]

My first PR to channelflow! :-) :-) :-)!!!