diff options
Diffstat (limited to 'src/spheres_poly/box.cpp')
| -rw-r--r-- | src/spheres_poly/box.cpp | 1157 |
1 files changed, 0 insertions, 1157 deletions
diff --git a/src/spheres_poly/box.cpp b/src/spheres_poly/box.cpp deleted file mode 100644 index 716090d..0000000 --- a/src/spheres_poly/box.cpp +++ /dev/null @@ -1,1157 +0,0 @@ -/* - 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]; - } - -} |