Arduino Plasma Simulation

/*
 * Example hybrid electrostatic Particle-in-Cell (ES-PIC) simulation running on an Arduino
 * More information at https://www.particleincell.com/2020/arduino-plasma-simulation/
 * The code simulates ions moving from left boundary to right where potential oscillates
 * Written by Lubos Brieda, April 2020
 * Tested with MKR Vidor 400
 */
 
// declare simulation parameters
const float x0 = 0.0;        // mesh origin
const float xm = 0.1;        // mesh extend
const int ni = 21;           // number of nodes
const float dx = xm/(ni-1);  // cell spacing
const float dt = 1e-6;       // time step
 
const float n0 = 1e9;        // reference number density
const float phi0 = 0;        // reference potential
const float kTe0 = 1;        // reference electron temperature
 
// physical constants
const float EPS_0 = 8.854187e-12;  // permittivity of free space
const float QE = 1.602e-19;        // elementary charge
const float M = 131*1.660539e-27;  // ion mass, 131 AMU
 
const int np_max = 1000;           // maximum number of simulation particles
const float w_mp = 1e4;            // macroparticle weight
 
float* ndi;             // ion number density array
float* phi;             // potential array
float* ef;              // electric field array
float* xp;              // particle positions
float* up;              // particle velocities
 
unsigned int np = 0;    // actual number of particles
unsigned long ts = 0;   // current time step
 
// random number generator, returns [0,1)
float rnd() {
  unsigned int rnd_max = 2147483647;
  float r;
  do { r=random(rnd_max)/(float)rnd_max; }
  while (r>=1.0);  // sample again if we get 1.0
  return r;
}
 
// this is the function that runs when the Arduino is powered up
void setup() {
  Serial.begin(57600);    // open serial port communication
  while(!Serial) {};      // wait for the initialization to complete  
 
  pinMode(LED_BUILTIN, OUTPUT);  // enable writing to the built in LED
  digitalWrite(LED_BUILTIN,0);   // turn off the LED
 
  ndi = (float*)malloc(ni*sizeof(float));  // allocate memory for ion denstity
  ef = (float*)malloc(ni*sizeof(float));   // allocate memory for electric field
  phi = (float*)malloc(ni*sizeof(float));  // allocate memory for potential
 
  for (int i=0;i<ni;i++) { 
    ndi[i] = 0;   phi[i] = 0;  ef[i] = 0;  // initialize to zero
  }
 
  xp = (float*)malloc(np_max*sizeof(float)); // allocate memory for ion positions
  up = (float*)malloc(np_max*sizeof(float)); // allocate memory for ion velocities
 
  for (int i=0;i<10;i++) {                   // blink the LED 10x to update status
    digitalWrite(LED_BUILTIN,0); delay(200);
    digitalWrite(LED_BUILTIN,1); delay(200);
  }
 
  // optionally, wait to receive 'go' from the serial port before proceeding
  if (0) {
    Serial.println("Done initializing, type 'go' to start");
 
    while(1) {
      if (Serial.available()>=2) {         // if we have at least two characters
        String str = Serial.readString();  // get the entire string
        Serial.println("read: "+str);      // write out what we received
        if (str[0]=='g' && str[1]=='o') break; // if we get 'g' and 'o', continue
      }
    }
  }
 
  Serial.println("Starting");   // write to serial port
  digitalWrite(LED_BUILTIN,0);  // turn off LED
}
 
// this is the code that runs continuosly once setup complets
// it is effective the PIC main loop
void loop() {
  int np_inject = 10;     // number of particles to inject at
  if (np+np_inject>np_max) // check for particle overflow
    np_inject = np_max-np; // adjust particle count if too many
 
  // inject np_inject particles
  for (int i=0;i<np_inject;i++) {
    xp[np] = x0;          
    do {
      up[np] = 10+50*rnd();  // pick random velocity, shoud sample from Maxwellian here!
                             // should also rewind for Leapfrog
    } while(up[np]<0);       // pick again if moving backwards
 
    np++;                    // update number of particles
  }
 
  // clear number density array
  for (int i=0;i<ni;i++) {ndi[i] = 0;} 
 
  // scatter particles to the grid to compute number density
  for (int p=0;p<np;p++) {
    float li = (xp[p]-x0)/dx;   // get logical coordinate
    int i = (int)li;            // truncate to get cell index
    float di = li-i;            // fraction distance in the cell
    ndi[i]+= w_mp*(1-di);       // scatter macroparticle weight to the nodes
    ndi[i+1]+=w_mp*(di);       
  }
 
  // divide by cell volume to get number density
  for (int i=0;i<ni;i++) ndi[i]/=dx;
  ndi[0]*=2.0; ndi[ni-1]*=2.0;     // only half-volume on boundaries
 
  // update potential on the right boundary
  const int freq_ts = 200;
  float f = (ts%freq_ts)/(float)freq_ts;    
  phi[0] = 0;
  phi[ni-1] = -10*cos(2*PI*f);
 
  // compute nonlinear potential, GS for now for simplicity  
  for (int k=0;k<100;k++) {   // perform 100 iterations, no convergence check
    for (int i=1;i<ni-1;i++) {
      float rho = QE*(ndi[i] - n0*exp((phi[i]-phi0)/(kTe0))); // charge density
      phi[i] = ((rho/EPS_0)*(dx*dx) + phi[i-1] + phi[i+1])/2; // update potential on node i
    }
  }
 
  // compute electric field
  ef[0] = (phi[0]-phi[1])/dx;          // one sided difference on left
  ef[ni-1] = (phi[ni-2]-phi[ni-1])/dx; // one sided difference on right
  for (int i=1;i<ni-1;i++) 
    ef[i] = (phi[i-1]-phi[i+1])/(2*dx); // central difference everywhere else
 
  // move particles
  for (int p=0;p<np;p++) {
    float li = (xp[p]-x0)/dx;      // logical coordinate
    int i = (int)li;               // cell index
    float di = li-i;               // fractional distance
    float ef_p = ef[i]*(1-di)+ef[i+1]*(di);  // electric field seen by the particle
    up[p] += QE*ef_p/M*dt;         // update velocity 
    xp[p] += up[p]*dt;             // upate position
 
    if (xp[p]<x0 || xp[p]>=xm) {   // check for particles leaving the domain
      xp[p] = xp[np-1];            // remove particle by replacing with the last one
      up[p] = up[np-1];            
      np--;                        // decrease number of particles
      p--;                         // decrement p to recheck this position again
    }
  }
 
  ts++;  // update time step
 
  // send field data over the serial port every 5 time steps
  // data sent as ts,999,np,999,phi,999,999,999,...,ndi,999,999,999,...
  if (ts%5==0) {         
     String str;
     str="ts,"+String(ts);          // time step
     str+=",np,"+String(np);        // number of particles
     str+=",phi";                   // potential
     for (int i=0;i<ni;i++) str+=","+String(phi[i]);
     str+=",ndi";                   // number density
     for (int i=0;i<ni;i++) str+=","+String(ndi[i]);
     Serial.println(str);           // write out
  }  
 
  /*
  if (ts%100==0) {
    Serial.print(ts);
    Serial.print(":");
    Serial.println(millis()/1000.0);
  }
  */
}

#include <cstdlib>
#include <iostream>
#include <chrono>
#include <random>
 
using namespace std;
 
/* same parameters as in Arduino version */
const float PI = 3.1415926;
 
// holds starting time
chrono::time_point<chrono::high_resolution_clock> time_start;
 
// object for sampling random numbers
class Rnd {
public:
  //constructor: set initial random seed and distribution limits
  Rnd(): mt_gen{/*std::random_device()()*/0}, rnd_dist{0,1.0} {}
  void reseed(int seed) {mt_gen.seed(seed);}
  double operator() () {return rnd_dist(mt_gen);}
 
protected:
  std::mt19937 mt_gen;      //random number generator
  std::uniform_real_distribution<double> rnd_dist;  //uniform distribution
};
 
Rnd rnd;                    // make an rnd object
 
// this is the function that runs when the Arduino is powered up
void setup() {
  /* same as Arduino, minus Serial communication and LED blinking */
 
  //set time start
  time_start =  chrono::high_resolution_clock::now();	//time at simulation start
}
 
// this is the code that runs continuosly once setup complets
// it is effective the PIC main loop
void loop() {
  /* same as Arduino */
  if (ts%10000==0) {
    auto time_now = chrono::high_resolution_clock::now();
    chrono::duration<double> time_delta = time_now-time_start;
    cout<<ts<<":"<<time_delta.count()<<endl;
  }
}
 
// main adds the microcontroller loop functionality
int main() {
  setup();
  while(1) {loop();}
  return 0;
}