Running 1D-PIC simulation on Arduino-Nano
Sample output of 1D Electrostatic Particle-In-Cell code run on Arduino-Nano Every board (price ~10.33 USD)
/*
* 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
* Modified and Adapted for Nano-Every by Rupak Mukherjee, 2020
*/
// 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 = 1e5; // 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 = 100; // 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 = 2147;
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(9600); // 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
}
Serial.print(ts);
Serial.print(":");
Serial.println(millis()/50.0);
}