research/transmission_spectroscopy/universesim/src/update.c

   1 /*
   2  */
   3 #include <common.h>
   4 #include <stdio.h>
   5 #include <video.h>      // For the thread functions
   6
   7 // === CONSTANTS ===
   8 #define G_CONSTANT      200     // Gravitational Constant
   9 #define E_CONSTANT      400     // Electric Constant
  10 #define N_CONSTANT      50      // Nuclear Constant
  11 #define WN_CONSTANT     100     // Weak Nuclear Constant
  12
  13 #define MAX_THREADS     32
  14
  15 #define SHOW_FORCES     1
  16
  17 #define DEBUG   0
  18 #define DEBUG_SCAN      (DEBUG&&)
  19 #define DEBUG_CHECKS    (DEBUG&&0)
  20 #define DEBUG_INNER     (DEBUG&&0)
  21 #define DEBUG_STATE     (DEBUG&&0)
  22
  23 // === PROTOTYPES ===
  24 void    Update_DoTimestep(void);
  25 void    Update_ParticleForces(tParticle *Part);
  26 void    Update_ParticleVelocity(tParticle *Part);
  27 void    Update_ParticlePosition(tParticle *Part);
  28
  29 // === GLOBALS ===
  30  int    giNumberOfJobs = 16;
  31 volatile int    gaThreadStatuses[MAX_THREADS];  // -1 = Die, 0 = Completed, 1 = Working
  32
  33 // === CODE ===
  34 int Update_ProcessThread(void *Offset)
  35 {
  36          int    ofs = (intptr_t)Offset; // Also Thread ID
  37          int    i;
  38
  39         while( gaThreadStatuses[ofs] != -1 )
  40         {
  41                 while( gaThreadStatuses[ofs] == 0 )
  42                         Thread_Yield();
  43                 if( gaThreadStatuses[ofs] == -1 )
  44                         break;
  45                 for( i = ofs; i < giNumParticles; i += giNumberOfJobs )
  46                 {
  47                         Update_ParticleForces( &gpParticles[i] );
  48                 }
  49                 for( i = ofs; i < giNumParticles; i += giNumberOfJobs )
  50                 {
  51                         Update_ParticleVelocity( &gpParticles[i] );
  52                 }
  53                 gaThreadStatuses[ofs] = 0;
  54         }
  55         return 0;
  56 }
  57
  58 /**
  59  * Start update threads
  60  */
  61 void Update_Initialise(void)
  62 {
  63          int    i;
  64         for( i = 1; i < giNumberOfJobs; i++ ) {
  65                 Thread_Create(Update_ProcessThread, (void*)(intptr_t)i);
  66         }
  67 }
  68
  69 void Update_Cleanup(void)
  70 {
  71          int    i;
  72         for( i = 1; i < giNumberOfJobs; i++ ) {
  73                 gaThreadStatuses[i] = -1;
  74         }
  75 }
  76
  77 void Update_DoTimestep(void)
  78 {
  79          int    i;
  80
  81         // Start jobs again
  82         for( i = 1; i < giNumberOfJobs; i++ ) {
  83                 gaThreadStatuses[i] = 1;        // You have work to do
  84         }
  85
  86         for( i = 0; i < giNumParticles; i += giNumberOfJobs )
  87         {
  88                 #if DEBUG | SHOW_FORCES
  89                 printf(" Forces for #%i          \r", i);
  90                 fflush(stdout);
  91                 #endif
  92                 Update_ParticleForces( &gpParticles[i] );
  93         }
  94         #if SHOW_FORCES
  95         printf("                         \r");
  96         #endif
  97
  98         #if DEBUG
  99         printf(" Forces Done             \n");
 100         printf(" Velocities\n");
 101         #endif
 102         for( i = 0; i < giNumParticles; i += giNumberOfJobs )
 103         {
 104                 #if DEBUG
 105                 printf(" Accel & Velocity for %i          %c", i, DEBUG_SCAN?'\r':'\n');
 106                 #endif
 107                 Update_ParticleVelocity( &gpParticles[i] );
 108         }
 109
 110         // Wait for jobs to terminate
 111         for( i = 1; i < giNumberOfJobs; i++ ) {
 112                 while( gaThreadStatuses[i] == 1 )
 113                         Thread_Yield();
 114         }
 115
 116         #if DEBUG
 117         printf(" Positions                        \n");
 118         #endif
 119         for( i = 0; i < giNumParticles; i++ )
 120         {
 121                 Update_ParticlePosition( &gpParticles[i] );
 122         }
 123 }
 124
 125 void Update_ParticleForces(tParticle *Part)
 126 {
 127          int    i;
 128         tLargeInt       dist, dist2, dist3, dist4;
 129         tLargeInt       force, tmpQuad;
 130         tVector line;
 131         tLargeInt       *partMass = &Part->Mass;
 132
 133         Vector_Zero(&Part->Forces);
 134
 135         for( i = 0; i < giNumParticles; i++ )
 136         {
 137                 if( &gpParticles[i] == Part )   continue;
 138                 #if DEBUG_CHECKS
 139                 printf("  Checking against %4i ", i);
 140                 # if DEBUG_INNER
 141                 printf("\n");
 142                 # else
 143                 printf("\r");
 144                 # endif
 145                 #endif
 146
 147                 Vector_Set( &line, &gpParticles[i].Location );
 148                 Vector_Sub( &line, &Part->Location );
 149                 Vector_Magnitude( &dist, &line );
 150                 LargeInt_PowNative( &dist2, &dist, 2 );
 151                 LargeInt_PowNative( &dist3, &dist, 3 );
 152                 LargeInt_PowNative( &dist4, &dist, 4 );
 153
 154                 #if DEBUG_INNER
 155                 printf("   Calculating Gravity\n");
 156                 #endif
 157                 // Gravitational Forces
 158                 // -(mM)*G/(dist^2)
 159                 LargeInt_Mul( &tmpQuad, NULL, partMass, &gpParticles[i].Mass );
 160                 //printf("   tmpQuad = %s\n", LargeInt_DumpEx(&tmpQuad, 0) );
 161                 LargeInt_MulNative( &tmpQuad, NULL, &tmpQuad, G_CONSTANT );
 162                 LargeInt_Div( &force, NULL, &tmpQuad, &dist2 );
 163                 LargeInt_Neg( &force, &force );
 164
 165                 #if DEBUG_INNER
 166                 printf("   Calculating Electric\n");
 167                 #endif
 168                 // Electric Forces
 169                 // - Opposites Attract, Like Repels, so no need to negate
 170                 // (e1e2)*E/(dist^2)
 171                 //printf("   Charges %s and %s\n",
 172                 //      LargeInt_DumpEx(&Part->Charge, 0),
 173                 //      LargeInt_DumpEx(&gpParticles[i].Charge, 1)
 174                 //      );
 175                 LargeInt_Mul( &tmpQuad, NULL, &Part->Charge, &gpParticles[i].Charge );
 176                 //printf("   tmpQuad = %s\n", LargeInt_DumpEx(&tmpQuad, 0) );
 177                 LargeInt_MulNative( &tmpQuad, NULL, &tmpQuad, E_CONSTANT );
 178                 LargeInt_Div( &tmpQuad, NULL, &tmpQuad, &dist2 );
 179                 LargeInt_Add( &force, &force, &tmpQuad );
 180
 181                 #if DEBUG_INNER
 182                 printf("   Calculating Nuclear\n");
 183                 #endif
 184                 // Nuclear Forces
 185                 // > Strong Nuclear (Attractive, 1/d^3)
 186                 // -(mM)*N/(dist^3)
 187                 LargeInt_Mul( &tmpQuad, NULL, &Part->Mass, &gpParticles[i].Mass );
 188                 LargeInt_MulNative( &tmpQuad, NULL, &tmpQuad, N_CONSTANT );
 189                 LargeInt_Div( &tmpQuad, NULL, &tmpQuad, &dist3 );
 190                 LargeInt_Neg( &tmpQuad, &tmpQuad );
 191                 LargeInt_Add( &force, &force, &tmpQuad );
 192                 // > Weak Nuclear (Repulsive, d^4)
 193                 // (mM)*N/(dist^4)
 194                 LargeInt_Mul( &tmpQuad, NULL, &Part->Mass, &gpParticles[i].Mass );
 195                 LargeInt_MulNative( &tmpQuad, NULL, &tmpQuad, WN_CONSTANT );
 196                 LargeInt_Div( &tmpQuad, NULL, &tmpQuad, &dist4 );
 197                 LargeInt_Add( &force, &force, &tmpQuad );
 198
 199                 #if DEBUG_INNER
 200                 printf("   Calculating final force\n");
 201                 #endif
 202                 #if DEBUG_INNER >= 2
 203                 printf("   Mul (%s,%s) *= %s\n",
 204                         LargeInt_DumpEx(&line.C[0], 0),
 205                         LargeInt_DumpEx(&line.C[1], 1),
 206                         LargeInt_DumpEx(&force, 2)
 207                         );
 208                 #endif
 209                 Vector_MulScalar( &line, &force );
 210
 211                 Vector_DivScalar( &line, &dist );
 212                 Vector_Add(&Part->Forces, &line);
 213         }
 214
 215         #if 0
 216         printf("   Part->Forces = (%s,%s)\n",
 217                 LargeInt_DumpEx(&Part->Forces.C[0], 0),
 218                 LargeInt_DumpEx(&Part->Forces.C[1], 1)
 219                 );
 220         #endif
 221
 222         #if DEBUG_CHECKS
 223         # if DEBUG_INNER
 224         printf("\n");
 225         # endif
 226         #endif
 227 }
 228
 229 void Update_ParticleVelocity(tParticle *Part)
 230 {
 231         tVector accel;
 232         #if DEBUG_STATE
 233         printf("   Set\n");
 234         #endif
 235         Vector_Set( &accel, &Part->Forces );
 236         #if DEBUG_STATE
 237         printf("   Div (%s,%s) /= %s\n",
 238                 LargeInt_DumpEx(&accel.C[0], 0),
 239                 LargeInt_DumpEx(&accel.C[1], 1),
 240                 LargeInt_DumpEx(&Part->Mass, 2)
 241                 );
 242         #endif
 243         Vector_DivScalar( &accel, &Part->Mass );
 244         #if DEBUG_STATE
 245         printf("   = (%s,%s)\n",
 246                 LargeInt_DumpEx(&accel.C[0], 0),
 247                 LargeInt_DumpEx(&accel.C[1], 1)
 248                 );
 249         printf("   Add\n");
 250         #endif
 251         Vector_Add( &Part->Velocity, &accel );
 252 }
 253
 254 void Update_ParticlePosition(tParticle *Part)
 255 {
 256         Vector_Add( &Part->Location, &Part->Velocity );
 257 }