fixing energy calculation. Test passes for 2D!

Author: leios

include/kernels.h

@@ -59,9 +59,12 @@ __global__ void is_eq(bool *a, bool *b, bool *ans);
 * @return	Magnitude of complex number
 */
 __device__ double complexMagnitude(double2 in);
+__global__ void complexMultiply(double2 *in1, double2 *in2, double2 *out);
+__host__ __device__ double2 complexMultiply(double2 in1, double2 in2);
 
 __device__ double2 make_complex(double in, int evolution_type);
 
+__global__ void complexAbsSum(double2 *in1, double2 *in2, double *out);
 __global__ void complexMagnitude(double2 *in, double *out);
 /**
 * @brief	Return the squared magnitude of a complex number. $|(a+\textrm{i}b)*(a-\textrm{i}b)|$
@@ -286,6 +289,8 @@ __global__ void zeros(bool *in, bool *out);
 
 __global__ void set_eq(double *in1, double *in2);
 
+__global__ void print_ds(double *vector);
+
 //##############################################################################
 /**
  * Non-implemented functions.

include/split_op.h

@@ -97,5 +97,6 @@ void optLatSetup(const std::shared_ptr<Vtx::Vortex> centre, const double* V,
 * @return	$\langle \Psi | H | \Psi \rangle$
 */
 double energy_angmom(double2 *gpuWfc, int gState, Grid &par);
+double energy_calc(Grid &par, double2* wfc);
 
 #endif

src/evolution.cu

@@ -191,8 +191,8 @@ void evolve(Grid &par,
             time_spent = (double) (end - begin) / CLOCKS_PER_SEC;
             printf("Time spent: %lf\n", time_spent);
             std::string fileName = "";
-            printf("ramp=%d        gstate=%d    rg=%d        \n",
-                   ramp, gstate, ramp | (gstate << 1));
+            //printf("ramp=%d        gstate=%d    rg=%d        \n",
+            //       ramp, gstate, ramp | (gstate << 1));
             switch (ramp | (gstate << 1)) {
                 case 0: //Groundstate solver, constant Omega value.
                 {
@@ -464,7 +464,6 @@ void evolve(Grid &par,
                         num_vortices[1] = num_vortices[0];
                         vortCoords->getVortices().swap(vortCoordsP->getVortices());
 		                vortCoords->getVortices().clear();
-			            std::cout << "I am here" << std::endl;
 
                     }
                     fileName = "wfc_ev";
@@ -488,7 +487,7 @@ void evolve(Grid &par,
             }
             //std::cout << "written" << '\n';
             if (par.bval("energy_calc")){
-                double energy = energy_angmom(gpuWfc,gstate, par);
+                double energy = energy_calc(par,gpuWfc);
                 // Now opening and closing file for writing.
                 std::ofstream energy_out;
                 std::string mode = "energyi.dat";

src/kernels.cu

@@ -159,6 +159,15 @@ __device__ double2 realCompMult(double scalar, double2 comp){
 __device__ double complexMagnitude(double2 in){
     return sqrt(in.x*in.x + in.y*in.y);
 }
+
+__global__ void complexAbsSum(double2 *in1, double2 *in2, double *out){
+    int gid = getGid3d3d();
+    double2 temp;
+    temp.x = in1[gid].x + in2[gid].x;
+    temp.y = in1[gid].y + in2[gid].y;
+    out[gid] = sqrt(temp.x*temp.x + temp.y*temp.y);
+}
+
 __global__ void complexMagnitude(double2 *in, double *out){
     int gid = getGid3d3d();
     out[gid] = sqrt(in[gid].x*in[gid].x + in[gid].y*in[gid].y);
@@ -186,6 +195,12 @@ __host__ __device__ double2 complexMultiply(double2 in1, double2 in2){
     return result;
 }
 
+__global__ void complexMultiply(double2 *in1, double2 *in2, double2 *out){
+    int gid = getGid3d3d();
+    out[gid] = complexMultiply(in1[gid], in2[gid]);
+}
+
+
 /*
 * Used to perform conj(in1)*in2; == < in1 | in2 >
 */
@@ -220,8 +235,8 @@ __global__ void cMultPhi(double2* in1, double* in2, double2* out){
 __global__ void vecMult(double2 *in, double *factor, double2 *out){
     double2 result;
     unsigned int gid = getGid3d3d();
-    result.x = (in[gid].x * factor[gid]);
-    result.y = (in[gid].y * factor[gid]);
+    result.x = in[gid].x * factor[gid];
+    result.y = in[gid].y * factor[gid];
     out[gid] = result;
 }
 
@@ -381,6 +396,7 @@ __global__ void scalarPow(double2* in, double param, double2* out){
 __global__ void vecConjugate(double2 *in, double2 *out){
     double2 result;
     unsigned int gid = getGid3d3d(); 
+    result.x = in[gid].x;
     result.y = -in[gid].y;
     out[gid] = result;
 }
@@ -465,7 +481,6 @@ __global__ void multipass(double* input, double* output){
     }
 }
 
-
 /*
 * Calculates all of the energy of the current state. sqrt_omegaz_mass = sqrt(omegaZ/mass), part of the nonlin interaction term
 */
@@ -650,6 +665,11 @@ __global__ void set_eq(double *in1, double *in2){
     in2[gid] = in1[gid];
 }
 
+__global__ void print_ds(double *vector){
+    int gid = getGid3d3d();
+    printf("%d\t%e\n",gid,vector[gid]);
+
+}
 
 //##############################################################################
 //##############################################################################

src/operators.cu

@@ -826,12 +826,10 @@ void generate_fields(Grid &par){
         cudaFree(V_gpu);
     }
     else{
-        par.store("V_gpu", V_gpu);
-        par.store("K_gpu", K_gpu);
+        par.store("V_gpu",V_gpu);
     }
 
     par.store("V",V);
-    //par.store("V_gpu",V_gpu);
     par.store("items", items);
     //par.store("items_gpu", items_gpu);
     par.store("wfc", wfc);

src/split_op.cu

@@ -241,6 +241,89 @@ void optLatSetup(std::shared_ptr<Vtx::Vortex> centre, const double* V,
     par.store("V_opt",v_opt);
 }
 
+double energy_calc(Grid &par, double2* wfc){
+    double* K = par.dsval("K_gpu");
+    double* V = par.dsval("V_gpu");
+
+    dim3 grid = par.grid;
+    dim3 threads = par.threads;
+
+    int xDim = par.ival("xDim");
+    int yDim = par.ival("yDim");
+    int zDim = par.ival("zDim");
+    int gsize = xDim*yDim*zDim;
+
+    double dx = par.dval("dx");
+    double dy = par.dval("dy");
+    double dz = par.dval("dz");
+    double dg = dx*dy*dz;
+
+    cufftHandle plan;
+
+    if (par.ival("dimnum") == 1){
+        plan = par.ival("plan_1d");
+    }
+    if (par.ival("dimnum") == 2){
+        plan = par.ival("plan_2d");
+    }
+    if (par.ival("dimnum") == 3){
+        plan = par.ival("plan_3d");
+    }
+
+    double renorm_factor = 1.0/pow(gsize,0.5);
+
+    double2 *wfc_c, *wfc_k;
+    double2 *energy_r, *energy_k;
+    double *energy;
+
+    cudaMalloc((void **) &wfc_c, sizeof(double2)*gsize);
+    cudaMalloc((void **) &wfc_k, sizeof(double2)*gsize);
+    cudaMalloc((void **) &energy_r, sizeof(double2)*gsize);
+    cudaMalloc((void **) &energy_k, sizeof(double2)*gsize);
+
+    cudaMalloc((void **) &energy, sizeof(double)*gsize);
+
+    // Finding conjugate
+    vecConjugate<<<grid, threads>>>(wfc, wfc_c);
+
+    // Momentum-space energy
+    cufftExecZ2Z(plan, wfc, wfc_k, CUFFT_FORWARD);
+    scalarMult<<<grid, threads>>>(wfc_k, renorm_factor, wfc_k);
+
+    vecMult<<<grid, threads>>>(wfc_k, K, energy_k);
+
+    cufftExecZ2Z(plan, energy_k, energy_k, CUFFT_INVERSE);
+    scalarMult<<<grid, threads>>>(energy_k, renorm_factor, energy_k);
+
+    cMult<<<grid, threads>>>(wfc_c, energy_k, energy_k);
+
+    // Position-space energy
+    vecMult<<<grid, threads>>>(wfc, V, energy_r);
+    cMult<<<grid, threads>>>(wfc_c, energy_r, energy_r);
+
+    complexAbsSum<<<grid, threads>>>(energy_r, energy_k, energy);
+
+    double *energy_cpu;
+    energy_cpu = (double *)malloc(sizeof(double)*gsize);
+
+    cudaMemcpy(energy_cpu, energy, sizeof(double)*gsize,
+               cudaMemcpyDeviceToHost);
+
+    double sum = 0;
+    for (int i = 0; i < gsize; ++i){
+        sum += energy_cpu[i]*dg;
+    }
+
+    free(energy_cpu);
+    cudaFree(energy_r);
+    cudaFree(energy_k);
+    cudaFree(energy);
+    cudaFree(wfc_c);
+    cudaFree(wfc_k);
+
+    return sum;
+}
+
 /**
 ** Calculates energy and angular momentum of current state.
 ** Implementation not fully finished.

src/unit_test.cu

@@ -388,9 +388,11 @@ void dynamic_test(){
 
     cudaMemcpy(array, array_gpu, sizeof(double)*n, cudaMemcpyDeviceToHost);
 
+/*
     for (int i = 0; i < n; ++i){
         std::cout << array[i] << '\n';
     }
+*/
 
     std::cout << "Dynamic tests passed" <<'\n';
 }
@@ -441,6 +443,8 @@ void bessel_test(){
 // These will be checked against 1d 
 void fft_test(){
 
+    std::cout << "Beginning cufft test.\n";
+
     // For these tests, we are assuming that the x, y and z dimensions are 
     // All the same (2x2x2)
     // Note that yDim needs to be singled out differently, but z/x need no loops
@@ -500,9 +504,11 @@ void fft_test(){
         exit(1);
     }
 
+/*
     for (int i = 0; i < gsize; i++){
         std::cout << array[i].x << '\t' << array[i].y << '\n';
     }
+*/
 
     // Now to try the inverse direction
 
@@ -521,9 +527,13 @@ void fft_test(){
         exit(1);
     }
 
+/*
     for (int i = 0; i < gsize; i++){
         std::cout << array[i].x << '\t' << array[i].y << '\n';
     }
+*/
+
+    std::cout << "cufft test passed!\n";
 
 
 
@@ -764,6 +774,8 @@ void grid_test3d(){
 // Test of the parSum function in 3d
 void parSum_test(){
 
+    std::cout << "Beginning test of parallel summation.\n";
+
     // Setting error
     cudaError_t err;
 
@@ -886,6 +898,8 @@ void parSum_test(){
         }
     }
 
+    std::cout << "Parallel summation test passed in 2 and 3D!\n";
+
 }
 
 // Test for the Grid structure with paramters in it
@@ -1141,7 +1155,7 @@ void evolve_test(){
     par.store("thresh_const", 1.0);
 
 
-    double thresh = 0.0001;
+    double thresh = 0.01;
     std::string buffer;
     int gsteps = 30001;
     int esteps = 30001;
@@ -1154,14 +1168,15 @@ void evolve_test(){
     par.store("omegaY", 1.0);
     par.store("esteps", esteps);
     par.store("gsteps", gsteps);
-    par.store("printSteps", 1000);
+    par.store("printSteps", 30000);
     par.store("write_file", false);
-    par.store("write_it", true);
+    par.store("write_it", false);
     par.store("energy_calc", true);
     par.store("box_size", 0.00007);
     par.store("yDim", 1);
     par.store("zDim", 1);
 
+
     // Running through all the dimensions to check the energy
     for (int i = 2; i <= 3; ++i){
         if (i == 2){
@@ -1173,6 +1188,10 @@ void evolve_test(){
         par.store("dimnum",i);
         init(par);
 
+        if (par.bval("write_file")){
+            FileIO::writeOutParam(buffer, par, "data/Params.dat");
+        }
+
         double omegaX = par.dval("omegaX");
         set_variables(par, 0);