Better zero crossing estimation for fifthwave

Had some difficulties getting decent zero crossing approximations for fifthwave. The main issue is that the data generated from the Starccm+ plots is quite messy in that it is rather zig-zaggy having some points very close to others, but not coincident or even the direction you might expect.

This has caused a number of problems with the code in the past (particularly when estimating wavelength and celerity which relies on good estimation of the zero crossings). The original code had a basic linear interpolation routine checking for a switch of sign and them some error checking after that to make sure it wasn’t an inflection. Unfortunately the messiness of the data made this process extremely unreliable.

Moving onto examining all the waves in the numerical wave tank (NWT) rather than just a region of it required much better estimation of the zero crossings in order to strongly identify which wave was which after a time step and thus, monitor its progress. I needed to smooth the data before it could be really useful and then best solution to this would seem to be splines.

Spline fitting is a pretty heavy bit of mathematics, but I am lucky that the authors of Scipy (the scientific Python library) created an interface to FITPACK (a FORTRAN curve fitting library) called scipy.interpolate. Thus generating splines was not too difficult, it’s just deciding what to do with them. This takes me on a slight tangent where I must express the discomfort caused by using something you don’t understand and the risks associated with that. I have a mathematics degree, but not enough time to learn splines again so I just merrily let the smoothing routine do it’s business and harvested the neat roots() method which gives the roots of the spline (or the zero crossings. Only until the results continued to be a bit rubbish did I actually investigate what was going on and this is what I found.

Basically the default spline is absolutely nowhere near the wave form. The smoothing parameter ‘s’ was set way to high (in fact, if it’s not specified its the length of the data i.e O(1000)). Ideally the smoothing in this first case should be O(0.01) which generates a results like this:

Close up though this still isn’t perfect and it may be that it won’t calculate this level of smoothing which means that I want to reduce the amount of data that the spline is fitted to, idealy to go from peak to trough to peak. To get a good estimate of the locations of the peaks and troughs, I create another graph of the 1st derivative of the above spline and fit another spline (with even lower s this time) and find it’s roots. This gives me the turning points of the original data. I do need to be careful here, however, as I want to get peak trough, so if there is another peak in there, I skip it till I find a trough. Finally I trim the data to a region between a peak and trough and fit a new spline (default smoothing is fine now), and voila, a good approximation to the zero crossing can be found. Here’s one I made earlier.

Finally, I’ll just put a quick summery of the code for getting these plots, etc, so that I can do it again if it all goes a bit airy.

# Import the fifthwave module
import fifthwave as fw
from numpy import linspace

# Set up the wave form data from a csv file
test_data = fw.TwoD_Wave_Form()
test_data.import_csv('mycsv.csv')

# Fit the first spline to that data and plot it
test_data.fit_spline(smoothing=0.01)
test_data.spline_plot('first_spline.png')

# Now build another wave form based on the gradient of the first spline
new_x = linspace(test_data.x[0], test_data.x[-1], len(test_data.x))
new_eta = test_data._spline(newx, 1)

work_data = fw.TwoD_Wave_Form(eta_list=new_eta, x_list=new_x)

# Get the spline approximation and the roots.
work_data.fit_spline(smoothing=0.001)
work_data.spline_plot('gradient_spline.png')

inflections = work_data._spline.roots()

# You can now filter the inflections to make sure they run peak to trough 
# Finally for each inflection move the work data and find the root, i.e there
# may be a crossing between the start of the data and the first inflection
fromindex = 0
toindex = test.closest_x(inflections[0]) + 1

work_data.x = test_data.x[fromindex:toindex]
work_data.eta = test_data.eta[fromindex:toindex]

work_data.fit_spline()
work_data.spline_plot('final_spline.png')

print 'Zero crossing is at ', work_data._spline.roots()[0]

# Note that the current smoothing factor is not easy to recover from the 
# spline, but in fact it is given by
test_data_spline_smoothing_factor = test_data._spline._data[6]

I'm a DOCTOR!

Well, not quite. I still need to pay off my dues and print and bind 7 copies(!) of the thing. Still, what a great feeling to have it finished. Viva and corrections felt a bit like going back to the bad old days of a one task life, but now it’s done it feels fabulous and I’m greatly looking forward to graduation, which will involve a KILT, no less.

Quite a bit has gone on actually since the last time I blogged but I’ve not really felt like I had the time to post. However, I did a time management course last Wednesday and I’ve worked out that Monday morning is a light year or two from my “Prime Time”, so I thought this might be a good opportunity to write something. Also, a friend just posted a blog and it made me feel guilty about never updating this, so twice the motivation.

So what’s the news? Well, fifthwave is nearly working and version 0.2 was released in time for my IES presentation, which explains a little about it and other things I do on the side. In fact, I might have a link to a pdf version of the presentation just click HERE to get the presentation. There is also a natty intro video to the presentation which is on youtube: I’ll embed it below.

Even with fifthwave coming on-line, there is still a lot of work to do. The primary job of a researcher is the production of papers, and I don’t have any to my name as yet. Thus, all attention is now focused on that. Hopefully fifthwave can give some backbone to one of the papers. The other is starting with a EWTEC submission which is based around the V&V report I did at the start of the job and then I’ll give it a canny twist of a hypothesis, to make it into a proper paper, hopefully. More on that to follow!

Right, lots to do! Need to plan plan plan! Oh, and check this out:

A youtube vid to be added…

Fixing Fifthwave

Ok, so the last two weeks of panic have been fixing something called fifthwave. To understand what fifthwave is, I need to introduce the project it is related to, that is actually supposed to be the most important part of my job. Strangly, this is only the second project I have taken on, but it has spawned a lot of work. Ironically, it is not even my part of the SuperGen project, and it’s doing something I never really wanted to do! Nonetheless, it has now occupied more than 2 years of my time.

What is this mysterious ball and chain, you may well ask? Simple, its modelling the impact of a current on a wave energy device. Now, this is easy to say, but the physics (in particular the numerical modelling challenge of doing this) is quite complicated. It is made exceptionally complicated by the fact that the current will be applied perpendicular to the wave, which means that it must be modelled in 3D. Argh! Add to this the fact that we have no in house code that could cope with this or the resources to produce one and you should be starting to see the challenges involved.

These issues are what have brought about ‘fifthwave’ as I was forced to use the commercial CFD package which is used to teach the students. Now to some extent this is reasonable package for our purposes – it generates waves and currents, it handles floating bodies with 6 degrees of freedom. Unfortunately, this was developed for preliminary investigations into ship dynamics and therefore lacks a great number of diagnostic tools for solving wave energy problems, in particular, testing the quality of the waves generated.

And that is what fifthwave is for – it is Python code that takes a waveform (in 2D) and records various metrics in graphs and reports about that wave form. Those metrics can then be used to do what is know as a convergence study that will inform me about the quality of the grid I am using. Very useful. What went wrong is that fifthwave originally just produced errors against an analytical wave form, but the convergence tests only work on the original data, and that was harder to extract than it might have been. After a two week slog, fifthwave can now output raw data and errors which makes more sense , say, if a wave with no analytical solution is being examined.

Fifthwave can do more too! In the course of developing fifthwave, a delivery system for sending CFD jobs to our compute cluster has also been added. Fifthwave is turning into an interesting little bit of code that might be useful to others (I have plans). If you want to keep an eye on how the development is going, I (used to) publish the state of the issue tracker (ditz, if you’re interested) here. Now I’d better fix the rest of the bugs and start using it for something.