update test with model call including Qprev, removed step-0, wrapped conductionT and added test against analytical version

[steo85it]

updates test with #20

[steo85it]

As mentioned in the commit, I set-up a wrapper for conductionT (leveraging thermal.pyx) and added a test against a port of the analytical version in Planetary-Code-Collection/master/Common/Test/test_Tprofile.m. The results are close but not equal: again, I'm not familiar with these models, so I might have introduced a stupid error while porting them.

[nschorgh]

The results are close but not equal:

How are they different? Roundoff or deviations at the bottom of the domain would be okay.

[sampotter]

Yes, re: what @nschorgh asked, it would be helpful if you could include the numbers here for quick reference

[steo85it]

yes, sorry, busy with weekend stuff ^^ @sampotter absolutely, I appreciate the opportunity to improve my coding! I'll fix this during the day, in between my plans+in the evening. @nschorgh here is the difference (the output you'd get printing the validation array from the test) after 120 steps over P = 670. * 88775.244 (seconds?), for each z layer

[ 1.8598130e-01  1.7000260e-01  1.5386400e-01  1.3744120e-01
  1.2030090e-01  1.0232100e-01  8.3528100e-02  6.3942700e-02
  4.3547300e-02  2.2308100e-02  1.9510000e-04 -2.2815200e-02
 -4.6742600e-02 -7.1606300e-02 -9.7424800e-02 -1.2421500e-01
 -1.5199120e-01 -1.8076420e-01 -2.1054170e-01 -2.4132770e-01
 -2.7312130e-01 -3.0591660e-01 -3.3970160e-01 -3.7445800e-01
 -4.1015960e-01 -4.4677160e-01 -4.8424990e-01 -5.2253970e-01
 -5.6157450e-01 -6.0127460e-01 -6.4154620e-01 -6.8227980e-01
 -7.2334930e-01 -7.6461010e-01 -8.0589860e-01 -8.4703080e-01
 -8.8780210e-01 -9.2798660e-01 -9.6733800e-01 -1.0055910e+00
 -1.0424638e+00 -1.0776626e+00 -1.1108883e+00 -1.1418456e+00
 -1.1702561e+00 -1.1958754e+00 -1.2185154e+00 -1.2380733e+00
 -1.2545664e+00 -1.2681757e+00 -1.2792951e+00 -1.2885880e+00
 -1.2970451e+00 -1.3060378e+00 -1.3173525e+00 -1.3331847e+00
 -1.3560534e+00 -1.3885744e+00 -1.4329965e+00 -1.4903547e+00]

Not so marginal, I'd say (but again, have a look at whether I used/wrapped/translated the code correctly).

[sampotter]

Overall, we shouldn't expect more than a bit of rounding error on the order of machine epsilon (~2.2*10^-16 for double precision), and even then with a bit of luck we should produce exactly the same floating point numbers regardless of whether we call Norbert's code from Fortran or Cython originally. So there is definitely some bug that needs to be fixed.

(And, as always... the bug may be in the code being tested or the test code itself. ;-) )

[steo85it]

ok, one thing I didn't have the time to do was to compare the, so to say, numerical and analytical results in Fortran/C, to get an idea of the differences we should expect. Will do (that should also help me see if I implemented something wrong).

[nschorgh]

P = 670. * 88775.244 is a Mars year in seconds. The deviations in the temperature profile increase with depth, but Fgeo is the same, so I can't guess why this deviates. I managed to write Python versions of conductionT/Q, which I'll push soon.

[nschorgh]

I added conductionT.py and conductionQ.py, which give the same results as the Fortran version. They are 150 times slower though. See https://github.com/nschorgh/Planetary-Code-Collection/tree/master/Common/Python

[sampotter]

In your test are you able to call solve_banded once and then reuse it for each triangle on the surface? If it gets recreated for each facet it will be pretty slow I think.

[steo85it]

Makes sense, I didn't make the connection. Anyway, whatever period would have been fine for testing. That's great, thanks: the "issue" when translating code to python is usually that to make it fast, one needs to avoid any for loop and vectorize as much as possible (applying 1 operation to all elements at once, instead of N times for N elements). Then one can usually get as fast as compiled languages (but yes, it's tricky sometimes with pre-existing code, since the logic is different).

[nschorgh]

In your test are you able to call solve_banded once and then reuse it for each triangle on the surface? If it gets recreated for each facet it will be pretty slow I think.

No. b[1] changes due to the nonlinear upper-boundary condition. Also, the user might want to use time-variable thermal properties, e.g., radiative heat transfer is temperature dependent. (I wrote conductionT2 and conductionQ2 for time-constant thermal properties, but they use private variables in a Fortran-specific way.)

[sampotter]

Ah, OK, didn't realize this. In that case I don't think the Python versions of conductionQ or conductionT are likely to be useful. The 150x slowdown isn't surprising to me. (The systems being solved simply aren't large enough to amortize Python's slowness.)

[nschorgh]

I'm sure someone with better Python skills than me can improve the performance of these functions. The slowness might also be due to the outermost 50000 time step loop. The function conductionQ.py might not be that slow by itself.

[steo85it]

I'm not sure what we'll be useful or not now: we are even too well equipped, having both the wrapper and the python code (it might be quite convenient to use the python code to validate the "faster" wrapper in the unittest, for example). :) In any case, I pushed the modifications @sampotter suggested (it's much cleaner now, especially the step functions).

[nschorgh]

Whichever option works for you. The new Python functions are just an alternative.

[sampotter]

Thanks @steo85it, the code looks good to me, now.

It's nice to have the extra Python versions to use as a groundtruth to validate against, but the best will be to get the Cython wrapper debugged and the tests passing.

Also, right now we have duplicated versions of C files from @nschorgh's PlanetaryCodeCollection in this repository, but eventually the ideal setup will be as described in #22. (Made an issue for this since it isn't a very high priority ATM but worth thinking about as we go.)

[steo85it]

I think I found the issue in the initialization of some parameters (starting from z[1:] instead of from the surface). I'll double check tomorrow and push the changes, but I'd say that now both the python code and the wrapper give the same results.

[steo85it]

As anticipated, I pushed my corrections to the conductionQ wrapper AND updates to the tests, which now run over the same 50k steps as the pcc/python version Norbert provided last week, and compare with success to his results for both Tprofile and Tsurface (conductionQ only, I didn't re-check conductionT and just commented its assertion). LMK, but for me this should be good to go (and to use for our applications), now. :)

The only thing that I'm not sure about, is whether the various commits (setting up the test, then modifying it) are maybe a bit messy. Still, they are separated, so that should be acceptable.

[steo85it]

Thanks for the comments, they should be ready to go (with the usual mess in the commits ... I'm still figuring out how to apply changes and amend commits without force pushes).

[nschorgh]

It seems dangerous to return z[1:] rather than just z[:] from setgrid (l.188 in thermal.pyx), because in the calling routine it will make z[0]>0. In conductionQ.c, z[0] is ignored and implicitly assumed to be zero.