Vlasov Solver [work in progress]

Discuss how polywell fusion works; share theoretical questions and answers.

Moderators: tonybarry, MSimon

olivier
Posts: 155
Joined: Thu Feb 14, 2008 5:21 pm
Location: Cherbourg, France

Post by olivier »

Tom Ligon wrote:That's a constant irratation in the working world ... MatLab and its $^*^% licenses!
You could try an open source equivalent (at home if not within your company!) such as Scilab which provides a good level of compatibility.
Scilab has a broad community, is developed by INRIA, a large research institute, is backed by industry and is a member of the Numerical Mathematics Consortium.
I have also heard of GNU Octave but have never tried.

kcdodd
Posts: 722
Joined: Tue Jun 03, 2008 3:36 am
Location: Austin, TX

Post by kcdodd »

As an update I have been doing debugging and continued making memory and speed optimization and I think its starting to get into a workable range.

I need to come up with a file format to get the geometry and configuration into the simulation, and then output the results. I'm thinking of using vtk format for outputting the data it looks simple, but am unfamiliar with the program itself. I have implemented particle sources, but I am still unsure how to do the sinks, since they will be just about every surface.

I'm excited! lol
Carter

drmike
Posts: 825
Joined: Sat Jul 14, 2007 11:54 pm
Contact:

Post by drmike »

I'm going to turn my data into *.obj format. That way anybody with a 3D viewer can rotate the view on their own computer and get a very good idea of what is going on.

I think picking a format that makes it possible for many people to view is a tough choice. But it sure is fun getting things work, so good luck!!

kcdodd
Posts: 722
Joined: Tue Jun 03, 2008 3:36 am
Location: Austin, TX

Post by kcdodd »

I thought obj was only a simple mesh formal
Carter

drmike
Posts: 825
Joined: Sat Jul 14, 2007 11:54 pm
Contact:

Post by drmike »

It has color, mesh, texture and grouping. I will just use the color and mesh.
I might use transparency, but I think just putting enough space between the vectors will work ok. It will be fun to experiment with.

kcdodd
Posts: 722
Joined: Tue Jun 03, 2008 3:36 am
Location: Austin, TX

Post by kcdodd »

I am starting to think maybe my method will not be efficient enough. Considering the scales one which the interesting stuff will happen and the refinement needed to see them. But I am continuing to work on it. I have finished the electric field solver and import xml for meshes. Here is a screenshot of the "projection" of the 6D phase mesh onto three dimensions. Projection in quotes because the points are not 1-1, there are many points in the 6D mesh which overlap one point in 3D. In the screenshot the mesh is maximally refined around the (0,0,0,0,0,0) phase coordinate, with ~500k divisions.

Image
Image
Carter

drmike
Posts: 825
Joined: Sat Jul 14, 2007 11:54 pm
Contact:

Post by drmike »

One thing you might try is Indrek's symmetry method. Only 1/48 of the volume actually needs to be modeled if you assume symmetry. There are 3 cut planes which give you 1/8th of a volume, then 3 cut planes through the corners that give you 1/6th of that. If you take the derivatives to be zero at the interfaces, you can fold the whole system around to get a final picture.

That factor of 48 can help your phase space tremendously. Same storage, more accuracy.

kcdodd
Posts: 722
Joined: Tue Jun 03, 2008 3:36 am
Location: Austin, TX

Post by kcdodd »

Ah, maybe a good idea. The 1/8 symmetry seems fairly easy to change since its still cubic. I'll have to think about 1/48 with what I have already done, but it already makes a bit of a difference! I think I will be using VisIt in the future for visualization, so hopefully it will look a bit prettier too. haha

Image

Image
Carter

drmike
Posts: 825
Joined: Sat Jul 14, 2007 11:54 pm
Contact:

Post by drmike »

Yup, already a huge savings in data - great increase in accuracy!

kcdodd
Posts: 722
Joined: Tue Jun 03, 2008 3:36 am
Location: Austin, TX

Post by kcdodd »

drmike, or anyone having experience with this, I am having an issue with QNAN apearing in my data. I'm not quite sure why and I have't been able to find any information on what operations actually produce this. Which makes it hard to trace the code producing it. Everything I have tried to reproduce it makes either IND or INF, but not QNAN. I'm using g++ on windows 2k with mingw.
Carter

drmike
Posts: 825
Joined: Sat Jul 14, 2007 11:54 pm
Contact:

Post by drmike »

Nope, sorry. A NAN usually comes from underflow or double overflow (INF/0 for example). It's a pain to find usually, I have to reduce the problem to specific zones and then print a hell of a lot to debug that kind of thing.

MSimon
Posts: 14335
Joined: Mon Jul 16, 2007 7:37 pm
Location: Rockford, Illinois
Contact:

Post by MSimon »

drmike wrote:Nope, sorry. A NAN usually comes from underflow or double overflow (INF/0 for example). It's a pain to find usually, I have to reduce the problem to specific zones and then print a hell of a lot to debug that kind of thing.
Interrupt on exception and go into the debug routine.

Ooops. Wrong domain.

Of course identifying the immediate source of a problem is no guarantee of finding the ultimate source.
Engineering is the art of making what you want from what you can get at a profit.

Indrek
Posts: 113
Joined: Wed Jul 04, 2007 1:51 pm
Location: Estonia
Contact:

Post by Indrek »

kcdodd wrote:drmike, or anyone having experience with this, I am having an issue with QNAN apearing in my data. I'm not quite sure why and I have't been able to find any information on what operations actually produce this. Which makes it hard to trace the code producing it. Everything I have tried to reproduce it makes either IND or INF, but not QNAN. I'm using g++ on windows 2k with mingw.
Not that I have any experience but google gave me this:

Code: Select all

indrek@home:/opt/soft/fpexc> cat test.cc

#include <fenv.h>
#include <math.h>

int main ()
{
  double a, b, c;
  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW );
  a = 0;
  b = 0;
  c = a / b;
  return 0;
}

indrek@home:/opt/soft/fpexc> g++ -g -Wall test.cc -o test -lm
indrek@home:/opt/soft/fpexc> gdb ./test
GNU gdb 6.5
Copyright (C) 2006 Free Software Foundation, Inc.
GDB is free software, covered by the GNU General Public License, and you are
welcome to change it and/or distribute copies of it under certain conditions.
Type "show copying" to see the conditions.
There is absolutely no warranty for GDB.  Type "show warranty" for details.
This GDB was configured as "x86_64-suse-linux"...Using host libthread_db library "/lib64/libthread_db.so.1".

(gdb) r
Starting program: /opt/soft/fpexc/test

Program received signal SIGFPE, Arithmetic exception.
0x0000000000400631 in main () at test.cc:11
11        c = a / b;
(gdb)
Not sure mingw has feenableexcept() though.

- Indrek

edit: feenableexcept(FE_INVALID); allows division by zero but gives you an exception on operation that produces a nan. So should be useful for filtering.

kcdodd
Posts: 722
Joined: Tue Jun 03, 2008 3:36 am
Location: Austin, TX

Post by kcdodd »

Well, I am having lots of problems. I'm sorry I can't post more positive thing, heh. My discretization algorithm seems to hang up at certain configurations, and the integration of density from phase space to position space doesn't seem to be working correctly.

In this picture the density is supposed to all be in the center. But when it integrates it, well, appears in off center in a place it shouldn't be. Plus if you notice there are negative values (ugh, to say the least).

http://www.andromedaspace.com/files/wrongdensity.png


I feel somewhat stupid right now on the whole thing since I can't find the problem after many hours of staring at the code. Here is an animation of the simulation I ran for 10 steps with an electron source at x=2. Of course you'll notice it shows up at x=4, but I already know that part is bugged of something. There is a polywell not rendered with the coils at 1.5 on each axis, at least its supposed to be there not sure what effect its having on the sim.

http://www.andromedaspace.com/files/firstsim.mpeg
Carter

tonybarry
Posts: 219
Joined: Sun Jul 08, 2007 4:32 am
Location: Sydney, Australia
Contact:

Post by tonybarry »

Hello Carter,
Do you need somebody to look at your Matlab code?

Regards,
Tony Barry

Post Reply