First pass is pretty "wow". This version of electron_fluid.c never went NAN or +/- INF on me, so that was one reason I decided to start there. The qd results were off by two orders of magnitude from the original electron_fluid.c results, but then I realized something else: I hadn't added the Intel rounding hack required! So that probably explains part of it. I also found one constant that should really have been configured as a quad double because it was fractional. Both of those together mean I need to move the constant builds to some point after I call fpu_fix_start().
I ran the previous incarnation of qd_ef with "time" for the following results on my 3-year-old i386 desktop machine:
real 380m21.669s
user 294m57.375s
sys 0m18.993s
So it took about 6:20. I'm going to re-run the new version with the above corrections and see if the results agree better. Here's the first few lines from the original electron_fluid.c:
Brent wrote:Kind of a newbie here. I wrote a few thoughts about the Polywell concept on my newly created weblog. I spoke some about modeling. In reality, this part is basically a summary of what I have been learning in an advanced thermodynamics course.
Thanks for the comment at IEC Fusion Tech.
Engineering is the art of making what you want from what you can get at a profit.
Just for giggles, I changed the output to write in straight ASCII in case I was missing something in the electron_fluid.c version. Still writes the same values.
Having re-run it again, it looks like there are a few significant changes, but they're all past digit 10 or so. The overall message is that the whole calculation took a jump to the left, and then a step to the right, and, well, you can forget about the Time Warp, but it sure did a results warp. I haven't looked too closely at the underlying numbers, but I have a bad feeling that all those divisions make this code very unstable. After reading those PDFs that MSimon linked to (IIRC in this thread) I have no idea where the problem really lies, and if it's necessary to have a PhD in math to figure it out, I'm hopelessly outclassed.
In any event -- I would like to take a crack at drmike's code that was causing underflows. I don't think you ever published it, did you, drmike? Shoot me a copy someplace and I'll see what I can do with it.
Finding ways to make particle codes stable is definitly PhD thesis work, if not life long carrer stuff. There are lots of good reasons controlled fusion hasn't been done yet. It's hard!!
Incidentally, I tried running the qd library on J-M Muller's recurrence from the "Mindless" paper MSimon pointed out above. Eventually, it, too, came up with the wrong convergence (100), albeit it took longer to get there than conventional double precision (started to blow up around step 45 rather than step 12 for regular double precision). But without a detailed analysis, it's not clear what we're calculating, or how deep the error noise goes, etc.
Well, after about an hour's worth of banging away on it, I clipped the wings of the electron_fluid_magnetic.c code so that the test section at the bottom is removed. I figure if we don't get NaN's by the end of the first section I can move on to run that, and if I do, we're hosed anyway. I also took the liberty of moving certain constants that were being recalculated multiple times in the integration process (pi/4 in particular, but I should probably have also included cos(pi/4)) as global quasi-constants. Anyway, I'm going to leave it run overnight and see what turns up in the morning. It may not be done even by then! That's a lot more kerchunking going on there than in the original version of electron_fluid.c.
Actually, I think I see a couple things I may have missed in the qd version of electron_fluid.c that may have caused some weirdness, a couple constants I needed to convert to qd_real's that are involved in divisions. I'll worry about that in the morning.
It's certainly a great experiment. I'm starting to write up the description of the brute force phase space math, it will be interesting to see if we need the quad precision for that as well.
the complexity of plasma is really mind boggling. following electrons in a vacuum is almost doable, and it should be interesting. but once you add ionization and scattering - yuck. "Challenging" becomes an understatement.
HOWEVER, I should add that I also clipped out the initialization files for magnetic and electric potential (which I did not have), so those values were all set to zero.
The potentials are computed in separate files. You'll probably have to make those quad precision too. That will change the loading from double size to quad size, but otherwise should be the same.
Well, huh. Looks like this may take a while longer as the potential function uses the GSL quadrature package (I'm starting to pick up all this lingo -- fun stuff!). I may have to retool the GSL package to use C++ and the qd libs.
Wait -- are you saying the magnetic_potential.c itself underflows?
I noticed there's a couple constants called "TOOSMALL" or something like that. Are those there to prevent underflow? 'Cause I'm wondering if I should just yank those and let the qd_real variables manage it themselves.
