particle simulation using fast multiple method

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

Moderators: tonybarry, MSimon

happyjack27
Posts: 1439
Joined: Wed Jul 14, 2010 5:27 pm

Re: particle simulation using fast multiple method

Post by happyjack27 »

runnable .jar file available for download at: http://kbtrader.info/plasma_sim.php

happyjack27
Posts: 1439
Joined: Wed Jul 14, 2010 5:27 pm

Re: particle simulation using fast multiple method

Post by happyjack27 »

* generic slider panel and slider frame classes are complete and tested.
now adding slider controls is a breeze! as in 3 lines of code!:

Code: Select all

globalSliders = new SliderFrame();
globalSliders.setProducer(elements, "colors");
globalSliders.show();
sample3png.png
sample3png.png (53.13 KiB) Viewed 8001 times

happyjack27
Posts: 1439
Joined: Wed Jul 14, 2010 5:27 pm

Re: particle simulation using fast multiple method

Post by happyjack27 »

more sliders...
sample4.png
sample4.png (177.77 KiB) Viewed 7986 times

happyjack27
Posts: 1439
Joined: Wed Jul 14, 2010 5:27 pm

Re: particle simulation using fast multiple method

Post by happyjack27 »

Code: Select all

coils : [
    {x:0,y:0,z:-1.0, nx:0,ny:0,nz:1,element_number:3,current:299792458,radius:1,segments:32},
    {x:0,y:0,z:-1.0, nx:0,ny:0,nz:1,element_number:4,current:0,radius:1,segments:32},
],
x,y,z, = center of coil
nx,ny,nz = normal to plane on which coil should be created.

will create a coil with the given elements with the given current (speed of travel around the coil), with the given number of segments.

much easier to e.g. create a polywell this way.

i might add more stuff like this, such as planes, etc.

i'm working on testing the magnetic fields now. not nearly as easy to verify as electric fields...

then i think i'm going to move on to making element "groups", that can be controlled together,
and then probably do some deep optimization, such as turning it into a sparse octree, with its depth automatically increasing or decreasing so that the bottom ("leaves") expand to another depth level if and only if they have more than x particles in them (i'm thinking 8 )

looks pretty good. i added "trails" that show the particle velocity vector - or actually, current, because the length is a function of the charge..

prestonbarrows
Posts: 78
Joined: Sat Aug 03, 2013 4:41 pm

Re: particle simulation using fast multiple method

Post by prestonbarrows »

If you wanted to do some trivial magnetic geometries, I could benchmark your results against some established software.

happyjack27
Posts: 1439
Joined: Wed Jul 14, 2010 5:27 pm

Re: particle simulation using fast multiple method

Post by happyjack27 »

...go on.... ....reveal to me these magical geometries...

prestonbarrows
Posts: 78
Joined: Sat Aug 03, 2013 4:41 pm

Re: particle simulation using fast multiple method

Post by prestonbarrows »

happyjack27 wrote:i'm working on testing the magnetic fields now. not nearly as easy to verify as electric fields...
I mean that if you get your magnetic fields working, run something simple like a cylindrical solenoid. Let us know the dimensions etc. and I'll run the same problem on established modeling software so you can compare it to your results for verification.

happyjack27
Posts: 1439
Joined: Wed Jul 14, 2010 5:27 pm

Re: particle simulation using fast multiple method

Post by happyjack27 »

They're working, but i'm just seeing particles moving around erratically. it's hard to tell by the particle movement if i have the equations and all that correct.

not sure how we'd "compare"...

but here's my math if anyone wants to double-check it...

Code: Select all

      float deltax = x-p.x;
      float deltay = y-p.y;
      float deltaz = z-p.z;
      if( deltax == 0 && deltay == 0 && deltaz == 0) return new float[]{0,0,0};
      float rds = (float)1.0/(deltax*deltax + deltay*deltay + deltaz*deltaz); //reciprocal distance squared
      float reciprocal_norm = (float)Math.sqrt(rds);
      
      float electrostatic_force = 
    		  -1*COLOUMB_MULTIPLIER*
    		  engine.elements[element_number].net_charge*quantity*
    		  engine.elements[p.element_number].net_charge*p.quantity*
    		  rds;

      //electric force vector
      float ex = electrostatic_force*deltax*reciprocal_norm;
      float ey = electrostatic_force*deltay*reciprocal_norm;
      float ez = electrostatic_force*deltaz*reciprocal_norm;
      
      float mx = 0;
      float my = 0;
      float mz = 0;
      if( !Octet.is2d) {
          //magnetic force vector
          float deltadx = dx-p.dx;
          float deltady = dy-p.dy;
          float deltadz = dz-p.dz;
          
          float bx = RC2*(deltady*ez - deltadz*ey);
          float by = RC2*(deltadz*ex - deltadx*ez);
          float bz = RC2*(deltadx*ey - deltady*ex);
          
          mx = -(deltady*bz - deltadz*by);
          my = -(deltadz*bx - deltadx*bz);
          mz = -(deltadx*by - deltady*bx);
      }
      if( !DO_MAGNETIC) {
    	  mx = 0;
    	  my = 0;
    	  mz = 0;
      }
      
      float[] ff = new float[]{
            ex+mx,
            ey+my,
            Octet.is2d ? 0 : ez+mz,
            };
      if( ff[0] != ff[0]) //if the x coordinate is "NaN" (out of range), don't do anything.
    	  return new float[]{0,0,0};
      return ff;
   }



prestonbarrows
Posts: 78
Joined: Sat Aug 03, 2013 4:41 pm

Re: particle simulation using fast multiple method

Post by prestonbarrows »

Here are the most basic tests for confirming your B-field is working correctly with particles.

Define a uniform B-field and send off a single particle with a known velocity. If everything is correct, you will get a circular orbit.
Check for

-Constant radius
-Plane of orbit perpendicular to both field and velocity
-Orbit is constant over many cycles (i.e. energy is conserved at a constant value because B-fields do no work)
-Radius and cycle period matches with theory for a number of field strengths and velocities

You can find the gyrofrequency and gyroperiod listed on page 28 in the NRL plasma formulary if you are unfamiliar.
http://wwwppd.nrl.navy.mil/nrlformulary ... ARY_07.pdf

You can repeat this with a constant E-field. In that case, you should get constant acceleration in direction of the field. This acceleration should be a=(q/m)*E.

Then, test a B-field parallel to an E-field. You should see a helical orbit with the same radius as the magnetic-only test and the center of motion should drift with the same acceleration as the electric-only test.

I am not very clear on how your algorithms work, but you will need timesteps much less than a gyroperiod so that the full orbits are resolved. This is true for any particle-based code unless there is some trickery I am not aware of.

happyjack27
Posts: 1439
Joined: Wed Jul 14, 2010 5:27 pm

Re: particle simulation using fast multiple method

Post by happyjack27 »

thanks much for the test geometries, preston.

for now i made a simple MHD thruster - a pair of charged lates and a pair of solenoids. the thrust seems to be going the wrong way. i would think it should be going out the sides - current X mag field - current is between the plates, mag field is between the solenoids (the solenoids are just a current loop) but in my test they are going out the solenoids (the top). as shown below.
sampl1.png
sampl1.png (7.21 KiB) Viewed 7923 times
so i'm guessing there's something wrong with my cross products on the mag field, in the code i posted above.

btw, as you can see, in addition to coils, i now added plates. i also added wires (not shown).

happyjack27
Posts: 1439
Joined: Wed Jul 14, 2010 5:27 pm

Re: particle simulation using fast multiple method

Post by happyjack27 »

prestonbarrows wrote: I am not very clear on how your algorithms work, but you will need timesteps much less than a gyroperiod so that the full orbits are resolved. This is true for any particle-based code unless there is some trickery I am not aware of.
yeah, it's lagrangian. so yeah, timestep thing. though i added an adjustment that seems to work well (oddx = ddx of last time step. initially the ddx's are actualy force vectors, so they're divided by mass to convert to accelleration vectors. the "adjustment" i'm tlaking about is the "adjustLastTimeStep" part. it makes it so that, instead of accelleration changing suddenly at each sampled time step, position and speed update is calculated as if accelleration changed linearly. this results, among other things, in not having to use as small of a time-step to accurately simulate orbits.):

Code: Select all

     //do initial setup calculations
      float dtt = dt*dt/2;
      float dttt = dt*dtt/3;
      float rdt = (float)1.0/dt;
      
      float mass = engine.elements[element_number].net_mass*quantity;
      //System.out.println("integrate "+ddx+" "+ddy+" "+ddz);
      
      //handle NaN's
      if( ddx != ddx) ddx = 0;
      if( ddy != ddy) ddy = 0;
      if( ddz != ddz) ddz = 0;

      ddx /= mass;
      ddy /= mass;
      ddz /= mass;
      
      float dddx = (ddx-oddx)*rdt;
      float dddy = (ddy-oddy)*rdt;
      float dddz = (ddz-oddz)*rdt;
      
      if( USE_ADJUST && !RELATIVISTIC) {
          //adjustLastTimeStep
          x += dddx*dttt;
          y += dddy*dttt;
          z += dddz*dttt;
            
          dx += dddx*dtt;
          dy += dddy*dtt;
          dz += dddz*dtt;
      }


as to "magical tricks", i'm not aware of anything either, but "vorticity" comes to mind. :-) ( http://en.wikipedia.org/wiki/Vorticity )

happyjack27
Posts: 1439
Joined: Wed Jul 14, 2010 5:27 pm

Re: particle simulation using fast multiple method

Post by happyjack27 »

i made it so it now shows the electric and mag field vectors on the particles (as colored lines coming off of the particles), and that helps immensely. now i see there appears to be an issue with the electric field vectors in the presence of current loops. i'm just using a circulating charge to create the magnetic field, and cancelling out its electric field with an inverse non-circulating charge. i'm guessing that might be messing it up, since it means really large charges. so i'm going separate the charge used for calculating current from the charge used for calculating electric field, so that i can have electromagnets with no electric field to begin with, instead of having to cancel it out. then i'll see how it goes from there. i might eventually add visualization to the octets so i can make sure they're updating correctly.

happyjack27
Posts: 1439
Joined: Wed Jul 14, 2010 5:27 pm

Re: particle simulation using fast multiple method

Post by happyjack27 »

Success!

That was the issue - the huge charges on the magnetic coils were leading to precision loss. now that they are chargeless, it works just as expected.

simple MHD thruster:
sampl2.png
sampl2.png (16.43 KiB) Viewed 7916 times
white = columnb force vector
red = b-field vector
blue = lorentz force vector
orange = "trail" (velocity vector, pointing backwards)

charged particles traveling perpendicular to charge and mag field, as expected.

happyjack27
Posts: 1439
Joined: Wed Jul 14, 2010 5:27 pm

Re: particle simulation using fast multiple method

Post by happyjack27 »

prestonbarrows wrote:Here are the most basic tests for confirming your B-field is working correctly with particles.

Define a uniform B-field and send off a single particle with a known velocity. If everything is correct, you will get a circular orbit.
Check for

-Constant radius - CHECK
-Plane of orbit perpendicular to both field and velocity - CHECK
-Orbit is constant over many cycles (i.e. energy is conserved at a constant value because B-fields do no work) - CHECK
-Radius and cycle period matches with theory for a number of field strengths and velocities - ehh, for that i'd have to calculate b field from a bunch of moving charged particles... i'm not real concerned, though, this would just mean its off by a constant factor, and it's trivial to just multiply by that factor to make it match.
e field checks out fine, as well. e field + b field work in combination per the MHD thruster test. only concern is I might be off by a constant factor on e-field, b-field, or both. which is fine -- i can tune that later if need be just by multiplying by a constant. (i used physical constants to known precision so presumably they'll be exact.) i've validated what i needed to validate. now i need to focus on optimizing.

.

happyjack27
Posts: 1439
Joined: Wed Jul 14, 2010 5:27 pm

Re: particle simulation using fast multiple method

Post by happyjack27 »

as part of my optimization, i'm moving lorentz force calculation to the integrateTimeStep function, so it's no longer per particle pair, now it's only per particle. the mag field is accumulated per particle pair, and then that is used in the integrateTimeStep function.

so now in the integrateTimeStep function i can make sure that energy is conserved by the magnetic field. (currently it slowly adds energy).
i'm taking the force vector from the lorentz force, and treating that as a centripetal force, and then using that to calculate a gyro-radius, and a radians-per-second.

now there's two things: a position update and a speed update. i might need to do some calculus... but i was thinking, well the position update will just be the speed update, integrated, which is just the force vector integrated. so basically i'd calculate the new position as if there were no other forces on the particle, and use that as the delta. and then do the same for speed. i add the new velocity vector and subtract the old one, to get my velocity vector delta. and then i add those to the delta from the electric force, and viola! a mag field that perfectly conserves energy! and in the absence of electric fields, has exact gyroradii to boot!

there's room to get more ambitious if i feel like it, with some caclulus...


bddx is mag field.
dx is particle velocity.
mddx is instantaneous lorentz force.

Code: Select all

      float b2 = bddx[0]*bddx[0]+bddx[1]*bddx[1]+bddx[2]*bddx[2];
      float rnb = (float)Math.sqrt(1.0/b2);

      float v2 = dx*dx+dy*dy+dz*dz;
      float v = (float)Math.sqrt(v2);
      float cpf2 = mddx[0]*mddx[0]+mddx[1]*mddx[1]+mddx[2]*mddx[2]; //centripetal force squared
      float gyro_radius = (float)(mass * v2 / Math.sqrt(cpf2));
      float radians_per_second = v/gyro_radius;
      float cpfrn = (float)Math.sqrt(1.0/cpf2);
      float[] gyro_center = new float[]{
    		  gyro_radius*mddx[0]*cpfrn+x,
    		  gyro_radius*mddx[1]*cpfrn+y,
    		  gyro_radius*mddx[2]*cpfrn+z,
      };
      float[] gyro_plane_norm = new float[]{
    		  bddx[0]*rnb,
    		  bddx[1]*rnb,
    		  bddx[2]*rnb,
      };
I am not very clear on how your algorithms work, but you will need timesteps much less than a gyroperiod so that the full orbits are resolved. This is true for any particle-based code unless there is some trickery I am not aware of.
I just invented said trickery! 8)

EDIT: changed "/rnb" to "*rnb".
Last edited by happyjack27 on Thu Jan 02, 2014 4:35 pm, edited 1 time in total.

Post Reply