diff options
Diffstat (limited to 'src/spheres_poly')
-rw-r--r-- | src/spheres_poly/box.cpp | 1157 | ||||
-rw-r--r-- | src/spheres_poly/box.h | 169 | ||||
-rw-r--r-- | src/spheres_poly/event.cpp | 65 | ||||
-rw-r--r-- | src/spheres_poly/event.h | 58 | ||||
-rw-r--r-- | src/spheres_poly/grid_field.h | 148 | ||||
-rw-r--r-- | src/spheres_poly/heap.cpp | 186 | ||||
-rw-r--r-- | src/spheres_poly/heap.h | 42 | ||||
-rw-r--r-- | src/spheres_poly/neighbor.cpp | 40 | ||||
-rw-r--r-- | src/spheres_poly/sphere.cpp | 66 | ||||
-rw-r--r-- | src/spheres_poly/sphere.h | 44 | ||||
-rw-r--r-- | src/spheres_poly/spheres.cpp | 64 | ||||
-rw-r--r-- | src/spheres_poly/spheres.h | 4 | ||||
-rw-r--r-- | src/spheres_poly/vector.h | 471 |
13 files changed, 2514 insertions, 0 deletions
diff --git a/src/spheres_poly/box.cpp b/src/spheres_poly/box.cpp new file mode 100644 index 0000000..716090d --- /dev/null +++ b/src/spheres_poly/box.cpp @@ -0,0 +1,1157 @@ +/* + Updated July 24, 2009 to include hardwall boundary condition option, as + well as polydisperse spheres. + + Packing of hard spheres via molecular dynamics + Developed by Monica Skoge, 2006, Princeton University + Contact: Monica Skoge (mskoge.princeton.edu) with questions + This code may be used, modified and distributed freely. + Please cite: + + "Packing Hyperspheres in High-Dimensional Euclidean Spaces" + M. Skoge, A. Donev, F. H. Stillinger and S. Torquato, 2006 + + if you use these codes. +*/ + +#include "box.h" +#include <iostream> +#include <math.h> +#include <stdlib.h> +#include <time.h> +#include <iomanip> + +//============================================================== +//============================================================== +// Class Box: Fills box with hardspheres to given packing fraction +// and evolves spheres using molecular dynamics! +//============================================================== +//============================================================== + + +//============================================================== +// Constructor +//============================================================== +box::box(int N_i, double r_i, double growthrate_i, double maxpf_i, + double bidispersityratio_i, double bidispersityfraction_i, + double massratio_i, int hardwallBC_i): + r(r_i), + N(N_i), + growthrate(growthrate_i), + h(N_i+1), + maxpf(maxpf_i), + bidispersityratio(bidispersityratio_i), + bidispersityfraction(bidispersityfraction_i), + massratio(massratio_i), + hardwallBC(hardwallBC_i) +{ + ngrids = Optimalngrids2(r); + cells.set_size(ngrids); + + s = new sphere[N]; + binlist = new int[N]; + x = new vector<DIM>[N]; + h.s = s; + + gtime = 0.; + rtime = 0.; + ncollisions = 0; + ntransfers = 0; + nchecks = 0; + neventstot = 0; + ncycles = 0; + xmomentum = 0.; + pressure = 0.; + collisionrate = 0.; + + cells.set_size(ngrids); + cells.initialize(-1); // initialize cells to -1 + srandom(::time(0)); // initialize the random number generator + for (int i=0; i<N; i++) // initialize binlist to -1 + binlist[i] = -1; + + time(&start); +} + + +//============================================================== +// Destructor +//============================================================== +box::~box() +{ + delete[] s; + delete[] binlist; + delete[] x; +} + + +//============================================================== +// ReadFile +//============================================================== +void box::ReadPositions(const char* filename) +{ + // open file to read in arrays + std::ifstream infile(filename); + + infile.ignore(256, '\n'); // ignore the dim line + infile.ignore(256, '\n'); // ignore the #sphere 1 line + infile.ignore(256, '\n'); // ignore the #sphere line + infile.ignore(256, '\n'); // ignore the diameter line + infile.ignore(1000, '\n'); // ignore the 100 010 001 line + infile.ignore(256, '\n'); // ignore the T T T line + + for (int i=0; i<N; i++) + { + infile >> s[i].r; // read in radius + infile >> s[i].gr; // read in growth rate + infile >> s[i].m; // read in mass + for (int k=0; k<DIM; k++) + infile >> s[i].x[k]; // read in position + } + infile.close(); +} + + +//============================================================== +// Recreates all N spheres at random positions +//============================================================== +void box::RecreateSpheres(const char* filename, double temp) +{ + ReadPositions(filename); // reads in positions of spheres + VelocityGiver(temp); // gives spheres initial velocities + AssignCells(); // assigns spheres to cells + SetInitialEvents(); +} + + +//============================================================== +// Creates all N spheres at random positions +//============================================================== +void box::CreateSpheres(double temp) +{ + int Ncurrent = 0; + for(int i=0; i<N; i++) + { + CreateSphere(Ncurrent); + Ncurrent++; + } + if (Ncurrent != N) + std::cout << "problem! only made " << Ncurrent << " out of " << N << " desired spheres" << std::endl; + + VelocityGiver(temp); + SetInitialEvents(); +} + + +//============================================================== +// Creates a sphere of random radius r at a random unoccupied position +//============================================================== +void box::CreateSphere(int Ncurrent) +{ + int keeper; // boolean variable: 1 means ok, 0 means sphere already there + int counter = 0; // counts how many times sphere already exists + vector<DIM> xrand; // random new position vector + double d = 0.; + double radius; + double growth_rate; + double mass; + int species; + + if (Ncurrent < bidispersityfraction*N) + { + radius = r; + growth_rate = growthrate; + mass = 1.; + species = 1; + } + else + { + radius = r*bidispersityratio; + growth_rate = growthrate*bidispersityratio; + mass = massratio; + species = 2; + } + + while (counter<1000) + { + keeper = 1; + + for(int k=0; k<DIM; k++) + xrand[k] = ((double)random()/(double)RAND_MAX)*SIZE; + + for (int i=0; i<Ncurrent; i++) // check if overlapping other spheres + { + d=0.; + if (hardwallBC) + { + for (int k=0; k<DIM; k++) + { + d += (xrand[k] - s[i].x[k])*(xrand[k] - s[i].x[k]); + } + } + else + { + for (int k=0; k<DIM; k++) + { + if ((xrand[k] - s[i].x[k])*(xrand[k] - s[i].x[k]) > SIZE*SIZE/4.) + { + if (xrand[k] > SIZE/2.) // look at right image + d += (xrand[k] - (s[i].x[k]+SIZE))* + (xrand[k] - (s[i].x[k]+SIZE)); + else // look at left image + d += (xrand[k] - (s[i].x[k]-SIZE))* + (xrand[k] - (s[i].x[k]-SIZE)); + } + else + d += (xrand[k] - s[i].x[k])*(xrand[k] - s[i].x[k]); + } + } + if (d <= (radius + s[i].r)*(radius + s[i].r)) // overlapping! + { + keeper = 0; + counter++; + break; + } + } + + if (hardwallBC) + { + for (int k=0; k<DIM; k++) // check if overlapping wall + { + if ((xrand[k] <= radius) || (SIZE - xrand[k] <= radius)) // touching wall + { + keeper = 0; + counter++; + break; + } + } + } + + if (keeper == 1) + break; + + if (counter >= 1000) + { + std::cout << "counter >= 1000" << std::endl; + exit(-1); + } + } + + // now convert xrand into index vector for cells + vector<DIM,int> cell; + cell = vector<DIM>::integer(xrand*((double)(ngrids))/SIZE); + + s[Ncurrent] = sphere(Ncurrent, xrand, cell, gtime, radius, growth_rate, + mass, species); + + //first check to see if entry at cell + if (cells.get(cell) == -1) //if yes, add Ncurrent to cells gridfield + cells.get(cell) = Ncurrent; + + else // if no, add i to right place in binlist + { + int iterater = cells.get(cell); // now iterate through to end and add Ncurrent + int pointer = iterater; + while (iterater != -1) + { + pointer = iterater; + iterater = binlist[iterater]; + } + binlist[pointer] = Ncurrent; + } + +} + + +//============================================================== +// Assign cells to spheres read in from existing configuration +//============================================================== +void box::AssignCells() +{ + for (int i=0; i<N; i++) + { + // now convert x into index vector for cells + vector<DIM,int> cell; + cell = vector<DIM>::integer(s[i].x*((double)(ngrids))/SIZE); + s[i].cell = cell; + + //first check to see if entry at cell + if (cells.get(cell) == -1) //if yes, add Ncurrent to cells gridfield + cells.get(cell) = i; + + else // if no, add i to right place in binlist + { + int iterater = cells.get(cell); // now iterate through to end and add Ncurrent + int pointer = iterater; + while (iterater != -1) + { + pointer = iterater; + iterater = binlist[iterater]; + } + binlist[pointer] = i; + } + } +} + + +//============================================================== +// Velocity Giver, assigns initial velocities from Max/Boltz dist. +//============================================================== +void box::VelocityGiver(double T) +{ + for (int i=0; i<N; i++) + { + for (int k=0; k<DIM; k++) + { + if (T==0.) + s[i].v[k] = 0.; + else + s[i].v[k] = Velocity(T); + } + } +} + + +//============================================================== +// Velocity, gives a single velocity from Max/Boltz dist. +//============================================================== +double box::Velocity(double T) +{ + double rand; // random number between -0.5 and 0.5 + double sigmasquared = T; // Assumes M = mass of sphere = 1 + double sigma = sqrt(sigmasquared); // variance of Gaussian + double stepsize = 1000.; // stepsize for discretization of integral + double vel = 0.0; // velocity + double dv=sigma/stepsize; + double p=0.0; + + rand = (double)random() / (double)RAND_MAX - 0.5; + if(rand < 0) + { + rand = -rand; + dv = -dv; + } + + while(fabs(p) < rand) // integrate until the integral equals rand + { + p += dv * 0.39894228 * exp(-vel*vel/(2.*sigmasquared))/sigma; + vel += dv; + } + return vel; +} + + +//============================================================== +// Finds next events for all spheres..do this once at beginning +//============================================================== +void box::SetInitialEvents() +{ + for (int i=0; i<N; i++) // set all events to checks + { + event e(gtime, i, INF); + s[i].nextevent = e; + h.insert(i); + } +} + + +//============================================================== +// Finds next event for sphere i +//============================================================== +event box::FindNextEvent(int i) +{ + double outgrowtime; + outgrowtime = (.5/ngrids - (s[i].r+s[i].gr*gtime))/s[i].gr + gtime; + + event t = FindNextTransfer(i); + event c = FindNextCollision(i); + + if ((outgrowtime < c.time)&&(outgrowtime < t.time)&&(ngrids>1)) + { + event o = event(outgrowtime,i,INF-1); + return o; + } + + if ((c.time < t.time)&&(c.j == INF)) // next event is check at DBL infinity + return c; + else if (c.time < t.time) // next event is collision! + { + CollisionChecker(c); + return c; + } + else // next event is transfer! + return t; +} + + +//============================================================== +// Checks events of predicted collision partner to keep collisions +// symmetric +//============================================================== +void box::CollisionChecker(event c) +{ + int i = c.i; + int j = c.j; + event cj(c.time,j,i,c.v*(-1)); + + // j should have NO event before collision with i! + if (!(c.time < s[j].nextevent.time)) + std::cout << i << " " << j << " error collchecker, s[j].nextevent.time= " << s[j].nextevent.time << " " << s[j].nextevent.j << ", c.time= " << c.time << std::endl; + + int k = s[j].nextevent.j; + if ((k < N) && (k!=i)) // j's next event was collision so give k a check + s[k].nextevent.j = INF; + + // give collision cj to j + s[j].nextevent = cj; + h.upheap(h.index[j]); +} + + +//============================================================== +// Find next transfer for sphere i +//============================================================== +event box::FindNextTransfer(int i) +{ + double ttime = dblINF; + int wallindex = INF; // -(k+1) left wall, (k+1) right wall + + vector<DIM> xi = s[i].x + s[i].v*(gtime - s[i].lutime); + vector<DIM> vi = s[i].v; + + for (int k=0; k<DIM; k++) + { + double newtime; + if (vi[k]==0.) + newtime= dblINF; + else if (vi[k]>0) // will hit right wall, need to check if last wall + { + if ((hardwallBC)&&(s[i].cell[k] == ngrids - 1)) + newtime = ((double)(s[i].cell[k]+1)*SIZE/((double)(ngrids)) + - (xi[k]+s[i].r+s[i].gr*gtime))/(vi[k]+s[i].gr); + else + newtime = ((double)(s[i].cell[k]+1)*SIZE/((double)(ngrids)) + - xi[k])/(vi[k]); + + if (newtime<ttime) + { + wallindex = k+1; + ttime = newtime; + } + } + else if (vi[k]<0) // will hit left wall + { + if ((hardwallBC)&&(s[i].cell[k] == 0)) + newtime = ((double)(s[i].cell[k])*SIZE/((double)(ngrids)) + - (xi[k]-(s[i].r+s[i].gr*gtime)))/(vi[k]-s[i].gr); + else + newtime = ((double)(s[i].cell[k])*SIZE/((double)(ngrids)) + - xi[k])/(vi[k]); + + if (newtime<ttime) + { + wallindex = -(k+1); + ttime = newtime; + } + } + } + + + if (ttime < 0) + ttime = 0; + // make the event and return it + event e = event(ttime+gtime,i,wallindex+DIM+N+1); + return e; +} + + +//============================================================== +// Check all nearest neighbor cells for collision partners +//============================================================== +void box::ForAllNeighbors(int i, vector<DIM,int> vl, vector<DIM,int> vr, + neighbor& operation) +{ + vector<DIM,int> cell = s[i].cell; + + // now iterate through nearest neighbors + vector<DIM, int> offset; // nonnegative neighbor offset + vector<DIM, int> pboffset; // nearest image offset + + vector<DIM,int> grid; + + int ii ; + + grid=vl; + while(1) + { + //if (vr[0] > 1) + //std::cout << grid << "..." << cell+grid << "\n"; + for(int k=0; k<DIM; k++) + { + offset[k]=grid[k]+ngrids; // do this so no negatives + if (cell[k]+grid[k]<0) //out of bounds to left + pboffset[k] = -1; + else if (cell[k]+grid[k]>=ngrids) // out of bounds to right + pboffset[k] = 1; + else + pboffset[k] = 0; + } + int j = cells.get((cell+offset)%ngrids); + while(j!=-1) + { + operation.Operation(j,pboffset); + j = binlist[j]; + } + + // A. Donev: + // This code makes this loop dimension-independent + // It is basically a flattened-out loop nest of depth DIM + for(ii=0;ii<DIM;ii++) + { + grid[ii] += 1; + if(grid[ii]<=vr[ii]) break; + grid[ii]=vl[ii]; + } + if(ii>=DIM) break; + } +} + + +//============================================================== +// PredictCollision +//============================================================== +void box::PredictCollision(int i, int j, vector<DIM, int> pboffset, + double& ctime, int& cpartner, + vector<DIM, int>& cpartnerpboffset) +{ + double ctimej; + + if (i!=j) + { + ctimej = CalculateCollision(i,j,pboffset.Double())+gtime; + + if (ctimej < gtime) + std::cout << "error in find collision ctimej < 0" << std::endl; + + if ((ctimej < ctime)&&(ctimej < s[j].nextevent.time)) + { + ctime = ctimej; + cpartner = j; + cpartnerpboffset = pboffset; + } + } +} + + +//============================================================== +// Find next collision +//============================================================== +event box::FindNextCollision(int i) +{ + collision cc(i, this); + + vector<DIM, int> vl, vr; + + for (int k=0; k<DIM; k++) // check all nearest neighbors + { + vl[k] = -1; + vr[k] = 1; + } + + ForAllNeighbors(i,vl,vr,cc); + + event e; + if (cc.cpartner == i) // found no collisions in neighboring cells + { + if (cc.ctime != dblINF) + std::cout << "ctime != dblINF" << std::endl; + e = event(dblINF,i,INF); // give check at double INF + } + else + e = event(cc.ctime,i,cc.cpartner,cc.cpartnerpboffset); + + return e; +} + + +//============================================================== +// Calculates collision time between i and image of j using quadratic formula +//============================================================== +double box::CalculateCollision(int i, int j, vector<DIM> pboffset) +{ + if ((hardwallBC)&&(pboffset.norm_squared() > 1E-12)) + { + return dblINF; + } + else + { + // calculate updated position and velocity of i and j + vector<DIM> xi = s[i].x + s[i].v*(gtime - s[i].lutime); + vector<DIM> vi = s[i].v; + vector<DIM> xj = s[j].x + pboffset*SIZE + s[j].v*(gtime - s[j].lutime); + vector<DIM> vj = s[j].v; + + double r_now_i = s[i].r + gtime*s[i].gr; + double r_now_j = s[j].r + gtime*s[j].gr; + + double A,B,C; + A = vector<DIM>::norm_squared(vi - vj) - (s[i].gr+s[j].gr)*(s[i].gr+s[j].gr); + B = vector<DIM>::dot(xi - xj, vi - vj) - (r_now_i+r_now_j)*(s[i].gr+s[j].gr); + C = vector<DIM>::norm_squared(xi - xj) - (r_now_i+r_now_j)*(r_now_i+r_now_j); + + if (C < -1E-12*(r_now_i+r_now_j)) + { + std::cout << "error, " << i << " and " << j << + " are overlapping at time "<< gtime << std::cout; + std::cout << "A, B, C = " << A << " " << " " << B << + " " << " " << C << std::endl; + if (CheckSphereDiameters()>0) + std::cout << "a sphere has grown greater than unit cell" << + std::endl; + else + std::cout << "unknown error" << std::cout; + exit(-1); + } + + return QuadraticFormula(A, B, C); + } +} + + +//============================================================== +// Quadratic Formula ax^2 + bx + c = 0 +//============================================================== + double box::QuadraticFormula(double a, double b, double c) +{ + double x = dblINF; + double xpos; + double xneg; + double det = b*b - a*c; + + if (c <= 0.) + { + if(b < 0.) // spheres already overlapping and approaching + { + //std::cout << "spheres overlapping and approaching" << std::endl; + //std::cout << "# events= " << neventstot << std::endl; + x = 0.; + } + } + else if (det > -10.*DBL_EPSILON) + { + if (det < 0.) // determinant can be very small for double roots + det = 0.; + if (b < 0.) + x = c/(-b + sqrt(det)); + else if ((a < 0.)&&(b > 0.)) + x = -(b + sqrt(det))/a; + else + x = dblINF; + } + return x; +} + + +//============================================================== +// Returns first event +//============================================================== +void box::ProcessEvent() +{ + neventstot++; + // Extract first event from heap + int i = h.extractmax(); + event e = s[i].nextevent; // current event + event f; // replacement event + + if ((e.j>=0)&&(e.j<N)) // collision! + { + ncollisions++; + //std::cout << "collision between " << e.i << " and " << e.j << " at time " << e.time << std::endl; + Collision(e); + f = FindNextEvent(i); + s[i].nextevent = f; + h.downheap(1); + if (f.time < e.time) + { + std::cout << "error, replacing event with < time" << std::endl; + exit(-1); + } + + // make sure collision was symmetric and give j a check + if ((s[e.j].nextevent.j != i)||(s[e.j].nextevent.time != gtime)) + { + std::cout << "error collisions not symmetric" << std::endl; + std::cout << "collision between " << e.i << " and " << e.j << " at time " << e.time << std::endl; + std::cout << "but " << e.j << " thinks it has " << s[e.j].nextevent.j<< " " << s[e.j].nextevent.time << std::endl; + exit(-1); + } + else // give j a check + s[e.j].nextevent.j = INF; + } + else if (e.j==INF) // check! + { + nchecks++; + //std::cout << "check for " << e.i << " at time " << e.time << std::endl; + f = FindNextEvent(i); + s[i].nextevent = f; + h.downheap(1); + } + else if (e.j==INF-1) // sphere outgrowing unit cell, decrease ngrids! + { + gtime = e.time; + Synchronize(false); + ngrids = ngrids - 1; +// std::cout << "need to reduce ngrids to " << ngrids << std::endl; + ChangeNgrids(ngrids); + h.downheap(1); + } + else // transfer! + { + ntransfers++; + //std::cout << "transfer for " << e.i << " at time " << e.time << std::endl; + Transfer(e); + f = FindNextEvent(i); + s[i].nextevent = f; + h.downheap(1); + //r = FindNextEvent(i, e.j-N-DIM-1); + //if (f.time <= e.time) + if (f.time < e.time) + { + std::cout << "error after transfer, replacing new event with < time" << " " << std::endl; + std::cout << std::setprecision(16) << "e.time= " << e.time << ", f.time= " << f.time << ", f.i= " << f.i << ", f.j= " << f.j << "e.i= " << e.i << ", e.j= " << e.j << std::endl; + std::cout << std::setprecision(16) << "difference= " << e.time - f.time << std::endl; + exit(-1); + } + } +} + + +//============================================================== +// Processes a collision +//============================================================= +void box::Collision(event e) +{ + double ctime = e.time; + int i = e.i; + int j = e.j; + vector<DIM,int> v = e.v; // virtual image + gtime = ctime; + + // Update positions and cells of i and j to ctime + s[i].x += s[i].v*(gtime-s[i].lutime); + s[j].x += s[j].v*(gtime-s[j].lutime); + + // Check to see if a diameter apart + double r_sum = s[i].r + s[j].r + (s[i].gr+s[j].gr)*gtime; + double distance = vector<DIM>::norm_squared(s[i].x - s[j].x- v.Double()*SIZE) - r_sum*r_sum; + if (distance*distance > 10.*DBL_EPSILON) + std::cout << "overlap " << distance << std::endl; + + s[i].lutime = gtime; + s[j].lutime = gtime; + + vector<DIM,double> vipar; // parallel comp. vi + vector<DIM,double> vjpar; // parallel comp. vj + vector<DIM,double> viperp; // perpendicular comp. vi + vector<DIM,double> vjperp; // perpendicular comp. vj + + // make unit vector out of displacement vector + vector<DIM,double> dhat; + dhat = s[i].x - s[j].x - v.Double()*SIZE; // using image of j!! + double dhatmagnitude = sqrt(dhat.norm_squared()); + dhat /= dhatmagnitude; + + vipar = dhat*vector<DIM>::dot(s[i].v, dhat); + vjpar = dhat*vector<DIM>::dot(s[j].v, dhat); + viperp = s[i].v - vipar; + vjperp = s[j].v - vjpar; + + s[i].v = vjpar + dhat*(s[i].gr+s[j].gr)*2 + viperp; + s[j].v = vipar - dhat*(s[i].gr+s[j].gr)*2 + vjperp; + + // momentum exchange + double xvelocity; // exchanged velocity + xvelocity = vector<DIM>::dot(s[i].v - s[j].v, dhat) - (s[i].gr+s[j].gr); + xmomentum += xvelocity*dhatmagnitude*s[i].m*s[j].m*2/(s[i].m+s[j].m); +} + + +//============================================================== +// Transfer, takes care of boundary events too +//============================================================= +void box::Transfer(event e) +{ + gtime = e.time; + int i = e.i; + int j = e.j; + int k=0; // dimension perpendicular to wall it crosses + + // update position and lutime (velocity doesn't change) + s[i].x += s[i].v*(gtime-s[i].lutime); + s[i].lutime = gtime; + + vector<DIM,int> celli; // new cell for i + celli = s[i].cell; // this is not redundant + + // update cell + if (j>N+DIM+1) // right wall + { + k = j-N-DIM-2; + celli[k] = s[i].cell[k] + 1; + + if (hardwallBC) + { + // if in right-most cell, reflect back + if (s[i].cell[k] == ngrids - 1) + s[i].v[k] = -s[i].v[k] - 4.*s[i].gr; + else + { + if (ngrids>1) + UpdateCell(i, celli); + } + } + else + { + // if in right-most cell, translate x and cell + if (s[i].cell[k] == ngrids - 1) + { + s[i].x[k] -= SIZE; + celli[k] -= ngrids; + } + if (ngrids>1) + UpdateCell(i, celli); + } + } + else if (j<N+DIM+1) // left wall + { + k = -j+N+DIM; + celli[k] = s[i].cell[k] - 1; + + if (hardwallBC) + { + // if in left-most cell, reflect back + if (s[i].cell[k] == 0) + s[i].v[k] = -s[i].v[k] + 4.*s[i].gr; + else + { + if (ngrids>1) + UpdateCell(i, celli); + } + } + else + { + // if in left-most cell, translate x and cell + if (s[i].cell[k] == 0) + { + s[i].x[k] += SIZE; + celli[k] += ngrids; + } + if (ngrids>1) + UpdateCell(i, celli); + } + } + else + std::cout << "error in Transfer" << std::endl; + +} + + +//============================================================== +// Updates cell of a sphere to time +//============================================================= +void box::UpdateCell(int i, vector<DIM,int>& celli) +{ + if (celli == s[i].cell) + std::cout << "error in update cell..shouldn't be the same" << std::endl; + + // delete i from cell array at cell + + if (cells.get(s[i].cell) == i) + { + if (binlist[i] == -1) + cells.get(s[i].cell) = -1; + else + { + cells.get(s[i].cell) = binlist[i]; + binlist[i] = -1; + } + } + + else if (cells.get(s[i].cell) == -1) + { + std::cout << "error " << i << " not in claimed cell UpdateCell" << std::endl; + OutputCells(); + } + + else // if no, find i in binlist + { + int iterater = cells.get(s[i].cell); + int pointer = iterater; + while ((iterater != i)&&(iterater != -1)) + { + pointer = iterater; + iterater = binlist[iterater]; + } + if (iterater == -1) // got to end of list without finding i + { + std::cout << "problem " << i << " wasn't in claimed, cell iterater = -1" << std::endl; + OutputCells(); + } + else // we found i! + { + binlist[pointer] = binlist[i]; + binlist[i] = -1; + } + } + + // now add i to cell array at celli + s[i].cell = celli; + + //first check to see if entry at celli + if (cells.get(celli) == -1) //if yes, add i to cells gridfield + cells.get(celli) = i; + else // if no, add i to right place in binlist + { + int iterater = cells.get(celli); // now iterate through to end and add i + int pointer = iterater; + while (iterater != -1) // find the end of the list + { + pointer = iterater; + iterater = binlist[iterater]; + } + binlist[pointer] = i; + binlist[i] = -1; // redundant + } +} + + +//============================================================== +// Output event heap...purely used for debugging +//============================================================== +void box::OutputEvents() +{ + h.print(); +} + + +//============================================================== +// Output positions of spheres and their cells...purely used for debugging +//============================================================== +void box::OutputCells() +{ + for (int i=0; i<N; i++) + std::cout << i << " " << s[i].x << " " << s[i].v << " " << s[i].cell << std::endl; +} + + +//============================================================== +// Update positions...purely for graphical display +//============================================================== +void box::TrackPositions() +{ + for (int i=0; i<N; i++) + x[i] = s[i].x + s[i].v*(gtime-s[i].lutime); +} + + +//============================================================== +// Computes the total energy +//============================================================== +double box::Energy() +{ + double E=0; + for (int i=0; i<N; i++) + E += 0.5*s[i].m*s[i].v.norm_squared(); + + return E/N; +} + + +//============================================================== +// Calculates the packing fraction +//============================================================== +double box::PackingFraction() +{ + double rfactor = 0.; + for (int i=0; i<N; i++) + { + rfactor += pow(s[i].r + gtime*s[i].gr, DIM); + } + double v = (rfactor*pow(sqrt(PI), DIM))/(exp(lgamma(1.+((double)(DIM))/2.))); + return v/(pow(SIZE, DIM)); +} + + +//============================================================== +// Checks to make sure all sphere diameters are less than dimension +// of unit cell +//============================================================== +int box::CheckSphereDiameters() +{ + int offender = 0; + for (int i=0; i<N; i++) + { + if (s[i].r*2 > 1./ngrids){ + offender = i; + break; + } + } + return offender; +} + + +//============================================================== +// Change ngrids +//============================================================== +void box::ChangeNgrids(int newngrids) +{ + cells.set_size(newngrids); + cells.initialize(-1); // initialize cells to -1 + for (int i=0; i<N; i++) // initialize binlist to -1 + binlist[i] = -1; + AssignCells(); + for (int i=0; i<N; i++) + s[i].nextevent = event(0., i, INF); + Process(N); +} + + +//============================================================== +// Calculates the optimal ngrids for the initial configuration +// and assumes that ngrids gets updated (reduced) as the packing +// proceeds +//============================================================== +int box::Optimalngrids2(double currentradius) +{ + return (int)(1./(2.*currentradius)); +} + + +//============================================================== +// Calculates the optimal ngrids, assuming ngrids is not updated +// automatically and is very conservative +//============================================================== +int box::Optimalngrids() +{ + double maxr; + + maxr = pow(exp(lgamma(1.+((double)(DIM))/2.))*maxpf/ + (N*(bidispersityfraction + (1.-bidispersityfraction)* + pow(bidispersityratio, DIM))), + 1./DIM)/sqrt(PI); + + return (int)(1./(2.*maxr)); +} + + +//============================================================== +// Processes n events +//============================================================== +void box::Process(int n) +{ + double deltat = gtime; + for (int i=0; i<n; i++) + { + ProcessEvent(); + } + pf = PackingFraction(); // packing fraction + deltat = gtime - deltat; + double oldenergy = energy; + energy = Energy(); // kinetic energy + + energychange = ((oldenergy - energy)/oldenergy)*100; // percent change in energy + + if (deltat != 0.) + { + pressure = 1+xmomentum/(2.*energy*N*deltat); + collisionrate = ((double)(ncollisions))/deltat; + } + + // reset to 0 + ncollisions = 0; + ntransfers = 0; + nchecks = 0; + xmomentum = 0.; + ncycles++; +} + + +//============================================================== +// Prints statistics for n events +//============================================================== +void box::PrintStatistics() +{ + std::cout << "packing fraction = " << pf << std::endl; + std::cout << "gtime = " << gtime << std::endl; + std::cout << "total time = " << rtime+gtime << std::endl; + std::cout << "kinetic energy = " << energy << std::endl; + std::cout << "total # events = " << neventstot << std::endl; + std::cout << "# events = " << ncollisions+ntransfers+nchecks << ", # collisions = " << ncollisions << ", # transfers = " << ntransfers << ", # checks =" << nchecks << std::endl; + std::cout << "growthrate = " << growthrate << std::endl; + std::cout << "collisionrate = " << collisionrate << std::endl; + std::cout << "reduced pressure = " << pressure << std::endl; + std::cout << "-----------------" << std::endl; +} + + +//============================================================== +// Updates spheres to gtime, synchronizes, and can change growth rate +//============================================================== +void box::Synchronize(bool rescale) +{ + double vavg = sqrt(2.*M*energy); + + for (int i=0; i<N; i++) + { + s[i].x = s[i].x + s[i].v*(gtime-s[i].lutime); + s[i].nextevent.time -= gtime; + + if (s[i].nextevent.time < 0.) + std::cout << "error, event times negative after synchronization" << std::endl; + if (rescale == true) // give everyone checks + { + s[i].nextevent = event(0., i, INF); + s[i].v /= vavg; + } + + s[i].lutime = 0.; + s[i].r += gtime*s[i].gr; + } + + //r += gtime*growthrate; // r defined at gtime = 0 + rtime += gtime; + gtime = 0.; + + if (rescale == true) + Process(N); +} + + +//============================================================== +// Run time +//============================================================== +void box::RunTime() +{ + time(&end); + std::cout << "run time = " << difftime(end, start) << std::endl; +} + + +//============================================================== +// Write configuration +//============================================================== +void box::WriteConfiguration(double *data) +{ + int count1; + + if (gtime != 0.) // synchronize spheres if not currently synchronized + Synchronize(false); + + + for (int i=0; i<N; i++) // output radius and positions + { + for (int k=0; k<DIM; k++) + data[DIM * i + k] = s[i].x[k]; + } + +} diff --git a/src/spheres_poly/box.h b/src/spheres_poly/box.h new file mode 100644 index 0000000..69418f4 --- /dev/null +++ b/src/spheres_poly/box.h @@ -0,0 +1,169 @@ +/* + Packing of hard spheres via molecular dynamics + Developed by Monica Skoge, 2006, Princeton University + Contact: Aleksandar Donev (adonev@math.princeton.edu) with questions + This code may be used, modified and distributed freely. + Please cite: + + "Packing Hyperspheres in High-Dimensional Euclidean Spaces" + M. Skoge, A. Donev, F. H. Stillinger and S. Torquato, 2006 + + if you use these codes. +*/ + +//----------------------------------------------------------------------------- +// Box maker +//--------------------------------------------------------------------------- + +#ifndef BOX_H +#define BOX_H + + +#include <vector> +#include <math.h> + +#include "vector.h" +#include "grid_field.h" +#include "event.h" +#include "sphere.h" +#include "heap.h" + + +#define PI 3.141592653589793238462643 +#define SIZE 1.0 // size of box +#define VOLUMESPHERE pow(PI,((double)(DIM))/2.)/exp(lgamma(1+((double)(DIM))/2.)) // volume prefactor for sphere +#define DBL_EPSILON 2.2204460492503131e-016 // smallest # such that 1.0+DBL_EPSILON!=1.0 +#define M 1.0 + +//--------------------------------------------------------------------------- +// Class neighbor +//--------------------------------------------------------------------------- +class neighbor +{ + public: + int i; + + neighbor(int i_i); + + public: + virtual void Operation(int j, vector<DIM, int>& pboffset) = 0; +}; + + +class box { + + public: + + // constructor and destructor + box(int N_i, double r_i, double growthrate_i, double maxpf_i, + double bidispersityratio, double bidispersityfraction, + double massratio, int hardwallBC); + ~box(); + + // Creating configurations + int Optimalngrids(); + int Optimalngrids2(double maxr); + void CreateSpheres(double temp); + void CreateSphere(int Ncurrent); + double Velocity(double temp); + void VelocityGiver(double temp); + void SetInitialEvents(); + void RecreateSpheres(const char* filename, double temp); + void ReadPositions(const char* filename); + void AssignCells(); + + // Predicting next event + event FindNextEvent(int i); + void CollisionChecker(event c); + event FindNextTransfer(int i); + event FindNextCollision(int i); + void ForAllNeighbors(int, vector<DIM, int>, vector<DIM,int>, neighbor&); + void PredictCollision(int i, int j, vector<DIM, int> pboffset, + double& ctime, int& cpartner, + vector<DIM, int>& cpartnerpboffset); + double CalculateCollision(int i, int j, vector<DIM> pboffset); + double QuadraticFormula(double a, double b, double c); + + // Processing an event + void Process(int n); + void ProcessEvent(); + void Collision(event e); + void Transfer(event e); + void UpdateCell(int i, vector<DIM,int>& celli); + void Synchronize(bool rescale); + void ChangeNgrids(int newngrids); + + // Debugging + void TrackPositions(); + void OutputEvents(); + void OutputCells(); + void GetInfo(); + int CheckSphereDiameters(); + + // Statistics + double Energy(); + double PackingFraction(); + void PrintStatistics(); + void RunTime(); + void WriteConfiguration(double* data); + + + //variables + + const int N; // number of spheres + + int ngrids; // number of cells in one direction + double maxpf; + double growthrate; // growth rate of the spheres + double r; // radius, defined at gtime = 0 + double gtime; // this is global clock + double rtime; // reset time, total time = rtime + gtime + double collisionrate; // average rate of collision between spheres + double bidispersityratio; // ratio of sphere radii + double bidispersityfraction; // fraction of smaller spheres + double massratio; // ratio of sphere masses + int hardwallBC; // =0 for periodic BC, =1 for hard wall + + + // statistics + double pressure; // pressure + double xmomentum; // exchanged momentum + double pf; // packing fraction + double energy; // kinetic energy + double energychange; + int ncollisions; // number of collisions + int ntransfers; // number of transfers + int nchecks; // number of checks + int neventstot; // total number of events + int ncycles; // counts # cycles for output + + time_t start, error, end; // run time of program + + // arrays + sphere *s; // array of spheres + grid_field<DIM, int> cells; // array that keeps track of spheres in each cell + int *binlist; // linked-list for cells array + heap h; // event heap + vector<DIM> *x; // positions of spheres.used for graphics +}; + + +//--------------------------------------------------------------------------- +// Predicts collisions, inherits neighbor operation +//--------------------------------------------------------------------------- +class collision : public neighbor +{ + public: + + box *b; + double ctime; + int cpartner; + vector<DIM,int> cpartnerpboffset; + + public: + collision(int i_i, box *b); + + virtual void Operation(int j, vector<DIM, int>& pboffset); +}; + +#endif diff --git a/src/spheres_poly/event.cpp b/src/spheres_poly/event.cpp new file mode 100644 index 0000000..c8c104d --- /dev/null +++ b/src/spheres_poly/event.cpp @@ -0,0 +1,65 @@ +#include "event.h" + +//============================================================== +//============================================================== +// Class Event +//============================================================== +//============================================================== + + +//============================================================== +// Constructor +//============================================================== +event::event(double time_i, int i_i, int j_i, vector<DIM,int> v_i): + time(time_i), + i(i_i), + j(j_i), + v(v_i) +{ +} + +event::event(double time_i, int i_i, int j_i): + time(time_i), + i(i_i), + j(j_i) +{ +} + +event::event(const event& e) +{ + time = e.time; + i = e.i; + j = e.j; + v = e.v; +} + +event::event() +{ +} + + +//============================================================== +// Destructor +//============================================================== +event::~event() +{ +} + +void event::erase() +{ + time = dblINF; + i = 0; + j = 0; +} + +bool event::operator<(const event& e) const +{ + return e.time < time; +} + +bool event::operator>(const event& e) const +{ + return e.time > time; +} + + diff --git a/src/spheres_poly/event.h b/src/spheres_poly/event.h new file mode 100644 index 0000000..55dd1fc --- /dev/null +++ b/src/spheres_poly/event.h @@ -0,0 +1,58 @@ +#ifndef EVENT_H +#define EVENT_H + +#include "vector.h" +#define INF 100000000 +#define dblINF 100000000. + +class event { + + public: + + // constructor and destructor + event(double time_i, int i_i, int j_i, vector<DIM,int> v_i); + event(double time_i, int i_i, int j_i); + event(const event& e); + event(); + + ~event(); + + bool operator<(const event&) const; + bool operator>(const event&) const; + void erase(); + + //variables + + + double time; // time of next collision + int i; // collision partner with lower number + int j; // collision partner with higher number + vector<DIM,int> v; // virtual image + + /* 0<=j<=N binary collision between i and j + j=N+DIM+1+x transfer where x=-(k+1) for left wall + and x=k+1 for right wall + j=INF both check after event that did not altered motion of i and check after event that altered motion of i, i.e rescaling of velocities. I currently don't see need t o separate the two + + j=-1 check after collision + + Virtual identifiers as scalars...I think bad idea, but here's my work + there will be easily fixed problems if DIM>=10 + -x<=v<=x x=k+1, image in k direction + v=xy x,y + =-xy -x,y + =-yx x,-y + =yx -x,-y + v=xyz x,y,z + =-xyz -x,y,z + =-yxz x,-y,z + =-zxy x,y,-z + =zyx -x,-y,z + =xzy x,-y,-z + =yzx -x,y,-z + =-zyx -x,-y,-z + */ + +}; + +#endif diff --git a/src/spheres_poly/grid_field.h b/src/spheres_poly/grid_field.h new file mode 100644 index 0000000..10430ed --- /dev/null +++ b/src/spheres_poly/grid_field.h @@ -0,0 +1,148 @@ +/* + Packing of hard spheres via molecular dynamics + Developed by Monica Skoge, 2006, Princeton University + Contact: Aleksandar Donev (adonev@math.princeton.edu) with questions + This code may be used, modified and distributed freely. + Please cite: + + "Packing Hyperspheres in High-Dimensional Euclidean Spaces" + M. Skoge, A. Donev, F. H. Stillinger and S. Torquato, 2006 + + if you use these codes. +*/ + + +#ifndef GRID_FIELD_H +#define GRID_FIELD_H + +#include "vector.h" + +// ====================================================================== +// grid_field +// ====================================================================== + +// A field of V-vectors on a D dimensional manifold + +template<int D, class T> +class grid_field { + + public: + int elements; + + private: + T* f; + vector<D, int> size; // number of grid points for each dimension + vector<D, int> offset; + + public: + + grid_field(); + grid_field(const vector<D, int>&); + ~grid_field(); + + T& get(const vector<D, int>&); + + vector<D, int> get_size() const; + void set_size(const vector<D, int>&); + void set_size(const int); + void initialize(const int i); +}; + + +// grid_field +// ~~~~~~~~~~~~ +template<int D, class T> +grid_field<D, T>::grid_field() + : f(0), elements(0) +{ +} + + +// grid_field +// ~~~~~~~~~~~~ +template<int D, class T> +grid_field<D, T>::grid_field(const vector<D, int>& s) + : f(0) +{ + set_size(s); +} + + +// ~grid_field +// ~~~~~~~~~~~~~ +template <int D, class T> +grid_field<D, T>::~grid_field() +{ + if(f != 0) + delete[] f; +} + + +// get_size +// ~~~~~~~~ +template<int D, class T> +inline vector<D, int> grid_field<D, T>::get_size() const +{ + return size; +} + + +// set_size +// ~~~~~~~~ +template<int D, class T> +void grid_field<D, T>::set_size(const vector<D, int>& s) +{ + if(f != 0) + delete[] f; + + size = s; + + elements = 1; + for(int i=0; i<D; i++) { + offset[i] = elements; + elements *= size.x[i]; + } + + f = new T[elements]; +} + + +// set_size +// ~~~~~~~~ +template<int D, class T> +void grid_field<D, T>::set_size(const int s) +{ + vector<D, int> square; + + for(int k=0; k<D; k++) + square[k] = s; + + set_size(square); +} + + +// get +// ~~~ +template<int D, class T> +inline T& grid_field<D, T>::get(const vector<D, int>& pos) +{ + int p=0; + for(int i=0; i<D; i++) + p += pos.x[i]*offset[i]; + + return f[p]; +} + + +// initialize +// ~~~ +template<int D, class T> +void grid_field<D, T>::initialize(const int value) +{ + for(int i=0; i<elements; i++) + f[i] = value; +} + + + +#endif diff --git a/src/spheres_poly/heap.cpp b/src/spheres_poly/heap.cpp new file mode 100644 index 0000000..c83865c --- /dev/null +++ b/src/spheres_poly/heap.cpp @@ -0,0 +1,186 @@ +#include "heap.h" +#include "event.h" +#include <iostream> +#include "box.h" + + +//============================================================== +//============================================================== +// Class heap: Event heap used in hard spheres calculation +//============================================================== +//============================================================== + + +//============================================================== +// Constructor +//============================================================== +heap::heap(int maxsize_i): + maxsize(maxsize_i) +{ + a = new int[maxsize]; + index = new int[maxsize]; + + N = 0; // current number of events in heap +} + + +//============================================================== +// Constructor +//============================================================== +heap::heap(const heap &h) +{ + maxsize = h.maxsize; + a = h.a; + index = h.index; + N = h.N; // current number of events in heap + s = h.s; +} + + +//============================================================== +// Destructor +//============================================================== +heap::~heap() +{ + delete[] a; + delete[] index; +} + +//============================================================== +// Upheap +//============================================================== +void heap::upheap(int k) +{ + int i = a[k]; + + while ((k>1)&&(s[a[k/2]].nextevent.time > s[i].nextevent.time)) + { + a[k] = a[k/2]; + index[a[k/2]] = k; + k = k/2; + } + a[k] = i; + index[i] = k; +} + +//============================================================== +// Downheap +//============================================================== +void heap::downheap(int k) +{ + int j; + int i = a[k]; + + while(k <= N/2) + { + j = k+k; + if ((j < N)&&(s[a[j]].nextevent.time > s[a[j+1]].nextevent.time)) + j++; + if (s[i].nextevent.time <= s[a[j]].nextevent.time) + break; + a[k] = a[j]; + index[a[j]] = k; + k = j; + } + a[k] = i; + index[i] = k; +} + +//============================================================== +// Insert +//============================================================== +void heap::insert(int i) +{ + if (N >= maxsize) + std::cout << "error, N >= maxsize, cannot insert another event" << std::endl; + else + { + N++; + a[N] = i; + index[i] = N; + upheap(N); + } +} + + +//============================================================== +// Extract max +//============================================================== +int heap::extractmax() +{ + return a[1]; +} + +/* +//============================================================== +// Replace +//============================================================== +void heap::replace(int i) +{ + a[1] = i; + s[i].nextevent = t; + + if (!(e.time > s[a[1]].nextevent.time)) + { + if (!(s[a[1]].nextevent.j == INF))// must be check i'm changing to coll. at same time + std::cout << "error replaced non check with earlier time" << std::endl; + a[1] = i; + index[i] = 1; + } + else + { + a[0] = i; + downheap(0); + } +} +*/ + +//============================================================== +// Search +//============================================================== +int heap::search(int j) +{ + for (int k=1; k<=N; k++) + { + if (a[k] == j) + return k; + } + return -1; +} + +/* +//============================================================== +// Change +//============================================================== +void heap::change(int i) +{ + int iindex = index[i]; + + if (s[i].nextevent.time == s[a[iindex]].nextevent.time) + std::cout << "error changing an event to an equal time" << std::endl; + else if (s[i].nextevent.time > s[a[iindex]].nextevent.time) + std::cout << "error changing an event to a greater time" << std::endl; + + a[iindex] = i; + upheap(iindex); +} +*/ + +//============================================================== +// Print +//============================================================== +void heap::print() +{ + for (int k=1; k<=N; k++) + std::cout << k << " " << a[k] << " " << s[a[k]].nextevent.j << " " << s[a[k]].nextevent.time << std::endl; +} + + +//============================================================== +// Check index array +//============================================================== +void heap::checkindex() +{ + for (int k=1; k<=N; k++) + std::cout << k << " " << a[k] << " " << index[a[k]] << std::endl; +} diff --git a/src/spheres_poly/heap.h b/src/spheres_poly/heap.h new file mode 100644 index 0000000..91180e3 --- /dev/null +++ b/src/spheres_poly/heap.h @@ -0,0 +1,42 @@ +//--------------------------------------------------------------------------- +// Event heap maker +//--------------------------------------------------------------------------- + +#ifndef HEAP_H +#define HEAP_H + +#include "event.h" +#include "sphere.h" + +class heap { + + public: + + // constructor and destructor + heap(int maxsize); + heap(const heap &h); + ~heap(); + + // variables + int maxsize; // max allowed number of events + int N; // current number of events + int *a; + sphere *s; + int *index; // array of indices for each sphere + //event minevent; + + + // functions which operate on a binary heap + + void upheap(int k); + void downheap(int k); + void insert(int i); + void replace(int i); + int search(int j); + void change(int i); + int extractmax(); + void print(); + void checkindex(); + +}; +#endif diff --git a/src/spheres_poly/neighbor.cpp b/src/spheres_poly/neighbor.cpp new file mode 100644 index 0000000..920c099 --- /dev/null +++ b/src/spheres_poly/neighbor.cpp @@ -0,0 +1,40 @@ +#include "vector.h" +#include <iostream> +#include "box.h" + + +//============================================================== +//============================================================== +// Class neighbor +//============================================================== +//============================================================== + +neighbor::neighbor(int i_i) + : i(i_i) +{ +} + + +//============================================================== +//============================================================== +// Class collision +//============================================================== +//============================================================== + +collision::collision(int i_i, box *b_i) + : neighbor(i_i), + b(b_i) +{ + ctime = dblINF; + cpartner = i; +} + + +//============================================================== +// Operation is finding the next collision from a given cell +//============================================================== +void collision::Operation(int j, vector<DIM, int>& pboffset) +{ + b->PredictCollision(i, j, pboffset, ctime, cpartner, cpartnerpboffset); +} + diff --git a/src/spheres_poly/sphere.cpp b/src/spheres_poly/sphere.cpp new file mode 100644 index 0000000..0fbbd68 --- /dev/null +++ b/src/spheres_poly/sphere.cpp @@ -0,0 +1,66 @@ +#include "box.h" +#include <iostream> +#include <math.h> +#include <stdlib.h> +#include <time.h> + +#include "vector.h" + +//============================================================== +//============================================================== +// Class Sphere: +//============================================================== +//============================================================== + + +//============================================================== +// Constructor +//============================================================== +sphere::sphere() +{ +} + + +//============================================================== +// Constructor +//============================================================== +sphere::sphere(const sphere& s) +{ + i = s.i; + x = s.x; + v = s.v; + cell = s.cell; + lutime = s.lutime; + nextevent = s.nextevent; + nextcollision = s.nextcollision; + r = s.r; + gr = s.gr; + m = s.m; + species=s.species; +} + +//============================================================== +// Constructor +//============================================================== +sphere::sphere(int i_i, vector<DIM> x_i, vector<DIM, int> cell_i, + double lutime_i, double r_i, double gr_i, double m_i, int species_i): + i(i_i), + x(x_i), + cell(cell_i), + lutime(lutime_i), + r(r_i), + gr(gr_i), + m(m_i), + species(species_i) +{ +} + +//============================================================== +// Destructor +//============================================================== +sphere::~sphere() +{ + +} + + diff --git a/src/spheres_poly/sphere.h b/src/spheres_poly/sphere.h new file mode 100644 index 0000000..28c64b4 --- /dev/null +++ b/src/spheres_poly/sphere.h @@ -0,0 +1,44 @@ +#ifndef SPHERE_H +#define SPHERE_H + + +#include "vector.h" + + +class sphere { + + public: + + // constructor and destructor + + sphere(); + sphere(const sphere& s); + sphere(int i_i, vector<DIM> x, vector<DIM, int> cell_i, double lutime_i, + double r_i, double gr_i, double m_i, int species_i); + ~sphere(); + + //variables + + int i; // sphere ID + + // impending event + event nextevent; // next event...can be collision or transfer + event nextcollision; // next collision if next event is transfer + // maybe nextnext event + + // past information + double lutime; // last update time + vector<DIM, int> cell; // cell that it belongs to + vector<DIM, double> x; // position + vector<DIM, double> v; // velocity + double r; // sphere radius + double gr; // sphere growth rate + double m; // sphere mass + int species; // species number (not used during the MD) + // make sure efficent in memory + + + +}; + +#endif diff --git a/src/spheres_poly/spheres.cpp b/src/spheres_poly/spheres.cpp new file mode 100644 index 0000000..9b14b01 --- /dev/null +++ b/src/spheres_poly/spheres.cpp @@ -0,0 +1,64 @@ +//=========================================================== +//=========================================================== +//=========================================================== +// +// Molecular dynamics simulation of hardspheres +// +//=========================================================== +//=========================================================== +//=========================================================== + +#include <iostream> +#include <math.h> +#include <fstream> +#include <vector> +#include <time.h> +#include <string.h> +#include <stdlib.h> +#include <cstdlib> + +#include "box.h" +#include "sphere.h" +#include "event.h" +#include "heap.h" + + +double *spherespp(int N, double bidispersityratio, double bidispersityfraction, double maxpf, double maxpressure, double maxcollisionrate) +{ + int eventspercycle = 20; // # events per sphere per cycle + double initialpf = 0.01; // initial packing fraction + double temp = 0.2; // initial temp., use 0 for zero velocities + double growthrate = 0.001; // growth rate + double massratio = 1.; // ratio of sphere masses + int hardwallBC = 0; // 0 for periodic, 1 for hard wall BC + + double d, r; // initial diameter and radius of spheres + + r = pow(initialpf*pow(SIZE, DIM)/(N*VOLUMESPHERE), 1.0/((double)(DIM))); + + box b(N, r, growthrate, maxpf, bidispersityratio, + bidispersityfraction, massratio, hardwallBC); + + b.CreateSpheres(temp); + + + while ((b.collisionrate < maxcollisionrate) && (b.pf < maxpf) && (b.pressure < maxpressure)) + { + b.Process(eventspercycle*N); + b.Synchronize(true); + } + + double *data = (double *)malloc(DIM * N * sizeof(double)); + + b.WriteConfiguration(data); + + return data; +} + +extern "C" { + double *spheres(int N, double bidispersityratio, double bidispersityfraction, double maxpf, double maxpressure, double maxcollisionrate) { + return spherespp(N, bidispersityratio, bidispersityfraction, maxpf, maxpressure, maxcollisionrate); + } +} + + diff --git a/src/spheres_poly/spheres.h b/src/spheres_poly/spheres.h new file mode 100644 index 0000000..44d8052 --- /dev/null +++ b/src/spheres_poly/spheres.h @@ -0,0 +1,4 @@ + +#pragma once + +double *spheres(int N, double bidispersityratio, double bidispersityfraction, double maxpf, double maxpressure, double maxcollisionrate); diff --git a/src/spheres_poly/vector.h b/src/spheres_poly/vector.h new file mode 100644 index 0000000..9ee481f --- /dev/null +++ b/src/spheres_poly/vector.h @@ -0,0 +1,471 @@ +#ifndef VECTOR_H +#define VECTOR_H + +#include <iostream> +#include <fstream> + +#define DIM 2 + +template <int D, typename T=double> +class vector { + + public: + T x[D]; + + public: + + vector(); + vector(const T[D]); + vector(const vector&); + ~vector(); + + vector<D, T>& operator+=(const vector<D, T>&); + vector<D, T>& operator-=(const vector<D, T>&); + vector<D, T>& operator*=(const T); + vector<D, T>& operator/=(const T); + vector<D, T> operator+(const vector<D, T>&) const; + vector<D, T> operator-(const vector<D, T>&) const; + vector<D, T> operator*(const T) const; + vector<D, T> operator/(const T) const; + vector<D, T> operator%(const T) const; + bool operator==(const vector<D, T> &a) const; + + vector<D, int> integer() const; + vector<D, double> Double() const; + static vector<D, int> integer(const vector<D, T>&); + static vector<D, double> Double(const vector<D, T>&); + + T& operator[](const unsigned int); + + double dot(const vector<D, T>&) const; + static double dot(const vector<D, T>&, const vector<D, T>&); + + double norm_squared() const; + static double norm_squared(const vector<D, T>&); + + void read(std::ifstream&); + void write(std::ofstream&) const; +}; + + +template <int D, typename T> +std::ostream& operator<<(std::ostream&, const vector<D, T>&); + + +// constructor +// ~~~~~~~~~~~ +template <int D, typename T> +vector<D, T>::vector() +{ + for(int k=0; k<D; k++) + x[k] = 0; +} + +template <int D, typename T> +vector<D, T>::vector(const T x_i[D]) +{ + for(int k=0; k<D; k++) + x[k] = x_i[k]; +} + +template <int D, typename T> +vector<D, T>::vector(const vector<D, T> &v) +{ + for(int k=0; k<D; k++) + x[k] = v.x[k]; +} + + +// destructor +// ~~~~~~~~~~ +template <int D, typename T> +vector<D, T>::~vector() +{ +} + + +// += +// ~~ +template <int D, typename T> +inline vector<D, T>& vector<D, T>::operator+=(const vector<D, T> &v) +{ + for(int k=0; k<D; k++) + x[k] += v.x[k]; + + return *this; +} + + +// -= +// ~~ +template <int D, typename T> +inline vector<D, T>& vector<D, T>::operator-=(const vector<D, T> &v) +{ + for(int k=0; k<D; k++) + x[k] -= v.x[k]; + + return *this; +} + + +// *= +// ~~ +template <int D, typename T> +inline vector<D, T>& vector<D, T>::operator*=(const T s) +{ + for(int k=0; k<D; k++) + x[k] *= s; + + return *this; +} + +// /= +// ~~ +template <int D, typename T> +inline vector<D, T>& vector<D, T>::operator/=(const T s) +{ + for(int k=0; k<D; k++) + x[k] /= s; + + return *this; +} + + +// + +// ~ +template <int D, typename T> +inline vector<D, T> vector<D, T>::operator+(const vector<D, T> &a) const +{ + vector<D, T> c; + + for(int k=0; k<D; k++) + c.x[k] = x[k] + a.x[k]; + + return c; +} + + +// - +// ~ +template <int D, typename T> +inline vector<D, T> vector<D, T>::operator-(const vector<D, T> &a) const +{ + vector<D, T> c; + + for(int k=0; k<D; k++) + c[k] = x[k] - a.x[k]; + + return c; +} + + +// * +// ~ +template <int D, typename T> +inline vector<D, T> vector<D, T>::operator*(const T s) const +{ + vector<D, T> c; + + for(int k=0; k<D; k++) + c[k] = x[k] * s; + + return c; +} + + +// / +// ~ +template <int D, typename T> +inline vector<D, T> vector<D, T>::operator/(const T s) const +{ + vector<D, T> c; + + for(int k=0; k<D; k++) + c[k] = x[k] / s; + + return c; +} + + +// == +// ~ +template <int D, typename T> +inline bool vector<D, T>::operator==(const vector<D, T> &a) const +{ + for(int k=0; k<D; k++) + { + if (!(x[k]==a.x[k])) + return false; + } + return true; +} + + +// % +// ~ +template <int D, typename T> +inline vector<D, T> vector<D, T>::operator%(const T s) const +{ + vector<D, T> c; + + for(int k=0; k<D; k++) + c[k] = x[k] % s; + + return c; +} + + +// integer +// ~~~~~~~ +template <int D, typename T> +inline vector<D, int> vector<D, T>::integer() const +{ + vector<D, int> c; + + for(int k=0; k<D; k++) + c[k] = (int)x[k]; + + return c; +} + +template <int D, typename T> +inline vector<D, int> vector<D, T>::integer(const vector<D, T>& v) +{ + return v.integer(); +} + + +// double +// ~~~~~~~ +template <int D, typename T> +inline vector<D, double> vector<D, T>::Double() const +{ + vector<D, double> c; + + for(int k=0; k<D; k++) + c[k] = (double)x[k]; + + return c; +} + +template <int D, typename T> +inline vector<D, double> vector<D, T>::Double(const vector<D, T>& v) +{ + return v.Double(); +} + + + +// [] +// ~~ +template <int D, typename T> +inline T& vector<D, T>::operator[](const unsigned int i) +{ + return x[i]; +} + + +// Dot +// ~~~ +template <int D, typename T> +inline double vector<D, T>::dot(const vector<D, T> &a) const +{ + double d=0; + + for(int k=0; k<D; k++) + d += x[k] * a.x[k]; + + return d; +} + +template <int D, typename T> +inline double vector<D, T>::dot(const vector<D, T> &a, const vector<D, T> &b) +{ + return a.dot(b); +} + + +// NormSquared +// ~~~~~~~~~~~ +template <int D, typename T> +inline double vector<D, T>::norm_squared() const +{ + return dot(*this, *this); +} + +template <int D, typename T> +inline double vector<D, T>::norm_squared(const vector<D, T>& v) +{ + return v.norm_squared(); +} + + +// read +// ~~~~ +template <int D, typename T> +void vector<D, T>::read(std::ifstream& in) +{ + in.read((char*)x, sizeof(T)*D); +} + +// write +// ~~~~~ +template <int D, typename T> +void vector<D, T>::write(std::ofstream& out) const +{ + out.write((const char*)x, sizeof(T)*D); +} + + + +// Insertion +// ~~~~~~~~~ +template <int D, typename T> +std::ostream& operator<<(std::ostream& os, const vector<D, T>& v) +{ + os << "("; + + for(int k=0; k<D-1; k++) + os << v.x[k] << ", "; + + os << v.x[D-1] << ")"; + + return os; +} + + +// ====================================================================== +// vector_field +// ====================================================================== + +// A field of V-vectors on a D dimensional manifold + +template<int V, int D, typename T=double> +class vector_field { + + public: + int elements; + + private: + vector<V, T>* f; + vector<D, int> size; // number of grid points for each dimension + vector<D, int> offset; + + public: + + vector_field(); + vector_field(const vector<D, int>&); + ~vector_field(); + + vector<D, int> get_size() const; + void set_size(const vector<D, int>&); + + vector<V, T>& get(const vector<D, int>&); + + void read(std::ifstream&); + void write(std::ofstream&) const; + + static void swap(vector_field<V, D, T>&, vector_field<V, D, T>&); +}; + + +// vector_field +// ~~~~~~~~~~~~ +template<int V, int D, typename T> +vector_field<V, D, T>::vector_field() + : f(0), elements(0) +{ +} + + +// vector_field +// ~~~~~~~~~~~~ +template<int V, int D, typename T> +vector_field<V, D, T>::vector_field(const vector<D, int>& s) + : f(0) +{ + set_size(s); +} + +// ~vector_field +// ~~~~~~~~~~~~~ +template <int V, int D, typename T> +vector_field<V, D, T>::~vector_field() +{ + if(f != 0) + delete[] f; +} + +// get_size +// ~~~~~~~~ +template<int V, int D, typename T> +inline vector<D, int> vector_field<V, D, T>::get_size() const +{ + return size; +} + +// set_size +// ~~~~~~~~ +template<int V, int D, typename T> +void vector_field<V, D, T>::set_size(const vector<D, int>& s) +{ + if(f != 0) + delete[] f; + + size = s; + + elements = 1; + for(int i=0; i<D; i++) { + offset[i] = elements; + elements *= size.x[i]; + } + + f = new vector<V, T>[elements]; +} + +// get +// ~~~ +template<int V, int D, typename T> +inline vector<V, T>& vector_field<V, D, T>::get(const vector<D, int>& pos) +{ + int p=0; + for(int i=0; i<D; i++) + p += pos.x[i]*offset[i]; + + return f[p]; +} + + +// read +// ~~~~ +template<int V, int D, typename T> +void vector_field<V, D, T>::read(std::ifstream& in) +{ + in.read((char*)f, elements*sizeof(T)*V); +} + + +// write +// ~~~~~ +template<int V, int D, typename T> +void vector_field<V, D, T>::write(std::ofstream& out) const +{ + out.write((const char*)f, elements*sizeof(T)*V); +} + +// swap +// ~~~~ +template<int V, int D, typename T> +void vector_field<V, D, T>::swap(vector_field<V, D, T>& v1, + vector_field<V, D, T>& v2) +{ + vector<V, T>* f; + + f = v1.f; + v1.f = v2.f; + v2.f = f; +} + + + +#endif |