bjodah / batemaneq

The Bateman equation in C++ and a Python wrapper thereof.
BSD 2-Clause "Simplified" License
4 stars 2 forks source link

U-238 Chain #11

Open Rolleroo opened 4 years ago

Rolleroo commented 4 years ago

Hi Björn,

I've been working with my friend Tobi on a decay calculator website using your module. The code is very rough atm, as you can see it is my first proper project! Its public in my repository although not currently working.

As a trial I put the U-238 decay chain into your excellent batemaneq module and ended up with daughter products which were a little odd. Values were higher than was possible i.e. disintegrations per second were higher for daugher products than the parent. I believe greater uncertainty is to be expected with bateman as the Thalf approach each other (U-238, U234 & Th230 Thalf = 4.9e9, 2.4e5 & 7.4e4y). The errors approach 4 significant figures on the 4th daugher product which I'd like to reduce.

I noticed in your C++ proof of the code you increased the precision of the calculations. Could a similar thing be done in python by using a 64 or 128 bit float? Example here?

I appreciate your help and time up to now and any further assistance is gratefully received. If you want me to send the code/results through let me know.

All the best, Kenny

bjodah commented 4 years ago

Hi Kenny,

Not sure if you're experiencing catastrophic cancellation or not, but you can use the pure Python implementation: https://github.com/bjodah/batemaneq/blob/master/batemaneq/bateman.py#L6

which can be used together with SymPy's arbitrary precision floats. That should give you a clue.

Unfortunately numpy's float128 is just a wrapper around 80-bit float on x86, which means that you'll only gains ~2 digits or so (epsilon for 64-bit is ~1e-16 and 80-bit has an epsilon of 1e-19 or so).

Using SymPy is something I've described in a blog post a few years ago: https://bjodah.github.io/blog/posts/bateman-equation.html

I actually only use this package for analytic solutions for decay chains for comparison with numerical solutions from ODE solvers (basically I use it in test suites), so using an ODE solver would be another option for you.

bjodah commented 4 years ago

Also, if you haven't already, I think you also should give PyNE a look. Might have what you're looking for already.

Rolleroo commented 4 years ago

Björn,

No catastrophic failure even in Thalf values which are quite close, just larger movement from correct values.

Thanks for the pointer to your blog post. ODE solvers is not something I have any experience with. I'll chat with my maths guru (AKA my dad) and see if this meets our criteria and if there is a simple implementation!

I did glance at Pyne but the sheer scale of it's reach did discourage me. Since you suggested I had a further look. Are you referring to the "transmute" module? Might this work even with a zero neutron flux? Moving a long way out of my comfort zone....

Thanks for the time and information again. All the best, Kenny

bjodah commented 4 years ago

Gotcha. I haven't actually used PyNE, but I did look into it at some point a few years back and I thought I remembered at least seeing some ENDF half-lives in there somewhere, not sure if they offer any convenience methods for calculating populations though.

All right, let me know if you can verify it's a bug with this package or not (I've used it to verify correctness of some ODE solvers, so in if there's a bug in here somewhere it'll probably be subtle since it hasn't manifested itself in those tests).

Rolleroo commented 4 years ago

Hi Bjorn,

Just thought I would let you know that I am going to be benchmarking your bateman code against Pyne and a third code, https://github.com/alexmalins/radioactivedecay along with the author.

So far the comparison has been very close we we are looking at doing something relatively rigorous using the ICRP 107 decay data https://github.com/pyne/pyne/blob/develop/examples/open_origen_data.ipynb

If you are interested, I'll let you know how it goes.

All the best, Kenny

P.S. Your code is now live in an alpha website. http://ec2-18-132-47-249.eu-west-2.compute.amazonaws.com/

bjodah commented 4 years ago

Hi Kenny,

Yes, that sounds interesting, so yes, please do let me know how it goes.

That web app of yours really do look good, great work!

Best regards, Björn