Particles.c

Writing a particle system with fields per particle. The example is a simulation of colliding galaxies.

 
#include <hdf5.h>
#include <math.h>
 
#include <F5/F5particles.h>
#include <F5/F5coordinates.h>
#include <F5/F5X.h>
 
#if 1
#define N       511
#define SHELLS  19
#define T_END   800
#elif 0
#define N       1002
#define SHELLS  30
#define T_END   800
#else
#define N       6
#define SHELLS  1
#define T_END   1
#endif
 
#ifndef M_PI
#define M_PI 3.141592
#endif
 
 
int main(int argc, char*argv[])
{
        /* Needed Variables */
F5_vec3_point_t P[N];
F5_vec3_float_t V[N]; 
F5_vec3_float_t A[N];
double          scalarValue[N];
float           sinPhi[N],
                mass[N],
                weirdness[N];
int             i;
hid_t           FileID;
double          t = 0, dt = .1, t_end = T_END; 
int             iteration = 0, out_every = 10; 
        
const double    G = 1; /* Gravitational Constant */
 
        P[0].x = P[0].y =  0;   P[0].z = 10;
        V[0].x = V[0].y = -0.02; V[0].z = -.05;
        mass[0] = .01;
        weirdness[0] = 1;
        scalarValue[0] = sinPhi[0] = 0;
 
        P[1].x = P[1].y = P[1].z = 0;
        V[1].x = V[1].y = V[1].z = 0;
        mass[1] = .04;
        weirdness[1] = -1;
        scalarValue[1] = sinPhi[1] = 1;
 
 
        for(i=2;i<N;i++) /* Init the Field */
        {
        int     ir = (i-2) / SHELLS,
                ip = (i-2) % SHELLS;
        double  r = 4.0    * (ir+1) / ( (N-2)/SHELLS),
                p = 2*M_PI * ip / SHELLS,
                sp = sin(p),
                cp = cos(p);
 
                P[i].x = r*sp; P[i].y = r*cp; P[i].z = 0;
                /*
                   v^2 / r = M / r^2
                   -->
                   |v| = M/r  Kepler rotation
                */
                V[i].x = cp * sqrt(G*mass[1]/ r); V[i].y = -sp * sqrt(G*mass[1]/ r); V[i].z = 0; 
 
                // some artificial formula to set up radially increasing mass 
                // this does not make sense physically but is for demonstration
                mass[i] = .0001*(r+1.0);
 
                weirdness[i] = r + .1*p;
 
                sinPhi[i] =
                scalarValue[i] = sp;
        }
 
 
        /* Save Data to F5 File - this is the real magic */
        /*-----------------------------------------------*/
        /* Create the File */
        FileID = H5Fcreate("Particles.f5",    /* The FileName */
                           H5F_ACC_TRUNC,    /* Creation Flags */
                           H5P_DEFAULT,      /* File creation property list identifier */
                           H5P_DEFAULT);     /* File access property list identifier */
 
 
        /* Time iteration (euler integration) */
        for(t = 0; t<t_end; t+= dt, iteration++)
        {
                for(i=0;i<N;i++) /* Outer Loop: for all particles */
                {
                int     j;
        #if 1
                int     Ngrav = 2;
        #else
                int     Ngrav = N;
        #endif
                        /* Inner loop: compute acceleration */
                F5_vec3_float_t a;
                        a.x = a.y = a.z = 0;
 
                        for(j=0; j<Ngrav; j++)
                        {
                        F5_vec3_float_t Rij;
                        double  R2, R3;
 
                                if (i==j) continue;
                        
                                Rij.x = P[i].x - P[j].x;
                                Rij.y = P[i].y - P[j].y;
                                Rij.z = P[i].z - P[j].z;
 
                                R2 = Rij.x*Rij.x + Rij.y*Rij.y + Rij.z*Rij.z;
                                R3 = R2 * sqrt(R2);
 
                                a.x += -G*mass[j] * Rij.x / R3;
                                a.y += -G*mass[j] * Rij.y / R3;
                                a.z += -G*mass[j] * Rij.z / R3;
                        }
 
                        /* Update velocities */
                        V[i].x += dt * a.x;
                        V[i].y += dt * a.y;
                        V[i].z += dt * a.z; 
 
                        /* Set acceleration */
                        A[i].x = a.x;
                        A[i].y = a.y;
                        A[i].z = a.z;                   
                }
                for(i=0;i<N;i++) /* Update positions */
                {
                        P[i].x += dt * V[i].x;
                        P[i].y += dt * V[i].y;
                        P[i].z += dt * V[i].z;
                } 
 
                if (iteration % out_every == 0)
                {
                        printf("T=%03.3lg      \r", t); 
 
                        /* Write the Data */ 
                        F5write_particle_cartesian3Dv(FileID, t*0.01, "GalaxyStars", P, N, 0,
                                                      "Velocity", F5T_VEC3_FLOAT, V,
                                                      "Acceleration", F5T_VEC3_FLOAT, A,
                                                      "mass", H5T_NATIVE_FLOAT, mass,
                                                      "weirdness", H5T_NATIVE_FLOAT, weirdness,
                                                      "sinusPhi0", H5T_NATIVE_DOUBLE, scalarValue,
                                                      "sinusPhi0f", H5T_NATIVE_FLOAT, sinPhi,
                                                      0);
                }               
        }
 
        /* Close the File */
        F5Xclose(FileID);
 
        return 0;
}