Initial commit with uploaded files.

Compute.cu

@@ -0,0 +1,140 @@
+#include "globals.h"
+#include "timing.h"
+#include "Compute.h"
+#include <sstream>
+
+void write_kspace_data(const char*, complex<float>*);
+__global__ void d_multiplyComplex(cufftComplex*, cufftComplex*,
+    cufftComplex*, int);
+__global__ void d_prepareDensity(int, float*, cufftComplex*, int);
+
+Compute::Compute() {};
+Compute::~Compute() {};
+
+void Compute::Initialize(string command) {
+    this->input_command = command;
+    this->compute_id = n_computes - 1;
+
+    istringstream iss(command);
+    string word;
+
+    // Put the first word, which should be "compute"
+    iss >> word;
+
+    // Detect compute style
+    iss >> word;
+    if (word == "avg_sk") {
+
+        // For avg_sk, will calculate the average S(k) for one type of particle.
+        // iparams[0] = particle type
+        // iparams[1] will track the number of calls to doCompute()
+
+        this->style = "avg_sk";
+
+
+        iss >> word;
+        this->iparams[0] = stoi(word) - 1;    // Store relevant type, shifting by one to zero index.
+        if (this->iparams[0] >= ntypes)
+            die("Invalid type to calculate avg_sk");
+
+        this->iparams[1] = 0;
+
+
+        // Set defaults for optional arguments
+        this->compute_wait = 0;
+        this->compute_freq = 100;
+
+
+        // Check for optional arguments
+        while (!iss.eof()) {
+            iss >> word;
+            if (word == "freq") {
+                iss >> word;
+                this->compute_freq = stoi(word);
+            }
+            else if (word == "wait") {
+                iss >> word;
+                this->compute_wait = stoi(word);
+            }
+        }
+        cout << "  Calculating <S(k)> for component " << iparams[0] + 1 << " every " << this->compute_freq <<
+            " steps after " << this->compute_wait << " steps have passed." << endl;
+
+    }// word == avg_sk
+
+    else {
+        die("Undefined compute style!");
+    }
+}
+
+
+void Compute::allocStorage() {
+
+    if ( this->style == "avg_sk" ) {
+
+        this->cpx.resize(M);
+        for (int i = 0; i < M; i++) 
+            this->cpx.at(i) = 0.0f;
+
+        cout << " this->cpx has initial size " << this->cpx.size() << " and capacity " << this->cpx.capacity() << endl;
+    }
+
+}
+void Compute::doCompute() {
+    compute_t_in = int(time(0));
+
+    if (this->style == "avg_sk") {
+
+        // Extract the density of the relevant type
+        d_prepareDensity<<<M_Grid, M_Block>>>(this->iparams[0], d_all_rho, d_cpx1, M);
+        check_cudaError("Compute->doCompute.style = avg_sk prepare density");
+
+        // fourier from d_cpx1 to d_cpx2 forward
+        cufftExecC2C(fftplan, d_cpx1, d_cpx2, CUFFT_FORWARD);
+        check_cudaError("Compute->doCompute.style = avg_sk cufftExec");
+
+
+        // Multiply by the complex conjugate and scale by 1/M
+        // Store it in d_cpx1 as the values inside are not needed at this point
+        d_multiplyComplex<<<M_Grid, M_Block>>> (d_cpx2, d_cpx2, d_cpx1, M);
+        check_cudaError("Compute->doCompute.style = avg_sk multiplyComplex");
+
+
+        // Copy data to host and store
+        // NOTE: this should probably be stored on the device and only 
+        // communicated when writing, but may be OK for now.
+        cudaMemcpy(cpx1, d_cpx1, M * sizeof(cufftComplex), cudaMemcpyDeviceToHost);
+
+        for (int i = 0; i < M; i++)
+            this->cpx.at(i) += cpx1[i].x + I * cpx1[i].y;
+
+        this->iparams[1] += 1;
+    }
+
+    compute_t_out = int(time(0));
+    compute_tot_time += compute_t_out - compute_t_in;
+
+}
+
+
+// Results are going to be written at log_freq intervals
+void Compute::writeResults(int compute_id) {
+
+
+    if (this->style == "avg_sk") {
+        for (int i = 0; i < M; i++) {
+            if (this->iparams[1] > 0)
+                k_tmp[i] = this->cpx[i] / float(this->iparams[1]);
+            else
+                k_tmp[i] = 0.0f;
+        }
+        char nm[50];
+
+        // compute_id is used in the name instead of "type" in case multiple 
+        // computes operate on the same type
+        sprintf(nm, "avg_sk_%d.dat", compute_id);
+        write_kspace_data(nm, k_tmp);
+    }
+}
+
+

Compute.h

@@ -0,0 +1,43 @@
+//////////////////////////////////////////////
+// Rob Riggleman                8/18/2021   //
+// Defining compute group to calculate      //
+// various quantities at regular intervals. //
+// Examples: Average S(k) for a group, avg  //
+// density field, etc.                      //
+//////////////////////////////////////////////
+
+#ifndef _COMPUTE
+#define _COMPUTE
+
+class Compute {
+private:
+    int GRID;               // Size of thread grid for cell operations
+    int BLOCK;              // Num of blocks per thread
+    string group_name;      // Name of the group in the n list
+    int group_index;        // Index of the group
+    string out_name;        // Name of the output file
+    string style;           // Compute style (e.g., avg_sk, avg_rho)
+    int iparams[5];         // Integer parameters
+    int compute_id;
+    float fparams[5];       // Float parameters
+    float* fstore1;         // Variable to store data
+    float* fstore2;         // Variable to store data
+    vector<complex<float>> cpx;   // Variable to store data
+
+
+public:
+
+
+    int compute_freq;       // Frequency for calling doCompute (timesteps), default 100
+    int compute_wait;       // Time to wait before calling doCompute (timesteps), default 0
+    string input_command;   // Full text of input command.
+    void Initialize(string);
+    void allocStorage(void);
+    void doCompute(void);
+    void writeResults(int);
+    Compute();
+    ~Compute();
+
+};
+
+#endif

array_utils.cu

@@ -0,0 +1,94 @@
+#include "globals.h"
+
+void unstack(int, int*);
+
+void get_r(int id, float* r) {
+    int i, * n = new int[Dim];
+
+    unstack(id, n);
+
+    for (i = 0; i < Dim; i++)
+        r[i] = dx[i] * float(n[i]);
+}
+
+
+
+
+float get_k(int id, float* k, int Dim) {
+    // id between 0 and M-1 (i value in loop), float k kx,ky,ky (that vecotor), Dim is dimensionality
+    // declare a vector for this
+    float kmag = 0.0f;
+    int i, *n;
+    n = new int[Dim];
+
+    unstack(id, n);
+
+    for (i = 0; i < Dim; i++) {
+        if (float(n[i]) < float(Nx[i]) / 2.)
+            k[i] = PI2 * float(n[i]) / L[i];
+
+        else
+            k[i] = PI2 * float(n[i] - Nx[i]) / L[i];
+
+        kmag += k[i] * k[i];
+    }
+    delete n;
+    return kmag;
+
+}
+
+
+// Receives index id in [0 , M ) and makes array
+// nn[Dim] in [ 0 , Nx[Dim] )
+void unstack(int id, int* nn) {
+
+    if (Dim == 1) {
+        nn[0] = id;
+        return;
+    }
+    else if (Dim == 2) {
+        nn[1] = id / Nx[0];
+        nn[0] = (id - nn[1] * Nx[0]);
+        return;
+    }
+    else if (Dim == 3) {
+        nn[2] = id / Nx[1] / Nx[0];
+        nn[1] = id / Nx[0] - nn[2] * Nx[1];
+        nn[0] = id - (nn[1] + nn[2] * Nx[1]) * Nx[0];
+        return;
+    }
+    else {
+        cout << "Dim is goofy!" << endl;
+        return;
+    }
+}
+
+
+void unstack_like_device(int id, int* nn) {
+
+    if (Dim == 1) {
+        nn[0] = id;
+        return;
+    }
+    else if (Dim == 2) {
+        nn[1] = id / Nx[0];
+        nn[0] = (id - nn[1] * Nx[0]);
+        return;
+    }
+    else if (Dim == 3) {
+        int idx = id;
+        nn[2] = idx / (Nx[1] * Nx[0]);
+        idx -= (nn[2] * Nx[0] * Nx[1]);
+        nn[1] = idx / Nx[0];
+        nn[0] = idx % Nx[0];
+        return;
+    }
+    else {
+        cout << "Dim is goofy!" << endl;
+        return;
+    }
+}
+
+
+
+

bonds.cu

@@ -0,0 +1,52 @@
+#include "globals.h"
+
+void host_bonds() {
+
+
+	int i, j, k, id1, id2, btype;
+	float mdr2, mdr, dr[3], delr, mf;
+
+	Ubond = 0.0f;
+	for (i = 0; i < n_P_comps; i++)
+		bondVir[i] = 0.f;
+
+	for (i = 0; i < ns; i++) {
+		id1 = i;
+
+		for (j = 0; j < n_bonds[id1]; j++) {
+			id2 = bonded_to[id1][j];
+
+			if (id2 < id1)
+				continue;
+
+			btype = bond_type[id1][j];
+
+			mdr2 = pbc_mdr2(x[id1], x[id2], dr);
+			mdr = sqrt(mdr2);
+			delr = mdr - bond_req[btype];
+
+			Ubond += delr * delr * bond_k[btype] ;
+
+			mf = 2.0f * bond_k[btype] * delr / mdr;
+			for (k = 0; k < Dim; k++) {
+				f[id1][k] -= mf * dr[k];
+				f[id2][k] += mf * dr[k];
+			}
+
+			bondVir[0] += -mf * dr[0] * dr[0];
+			bondVir[1] += -mf * dr[1] * dr[1];
+			if ( Dim == 2 )
+				bondVir[2] += -mf * dr[0] * dr[1];
+			else if (Dim == 3)
+			{
+				bondVir[2] += -mf * dr[2] * dr[2];
+				bondVir[3] += -mf * dr[0] * dr[1];
+				bondVir[4] += -mf * dr[0] * dr[2];
+				bondVir[5] += -mf * dr[1] * dr[2];
+			}
+		}//for ( j=0 ; j<n_bonds ;
+	}// for ( i=0 ; i<ns_loc
+
+	for (i = 0; i < n_P_comps; i++)
+		bondVir[i] *= 1.0 / V / float(Dim);
+}