Fix for #21 and #24

include/evolution.h

@@ -56,4 +56,16 @@ void evolve(Grid &par,
             unsigned int gstate,
             std::string buffer);
 
+void apply_gauge(Grid &par, double2 *wfc, double2 *Ax, double2 *Ay,
+                 double2 *Az, double renorm_factor_x,
+                 double renorm_factor_y, double renorm_factor_z, bool flip,
+                 cufftHandle plan_1d, cufftHandle plan_dim2,
+                 cufftHandle plan_dim3, double dx, double dy, double dz,
+                 double time, int yDim, int size);
+
+void apply_gauge(Grid &par, double2 *wfc, double2 *Ax, double2 *Ay,
+                 double renorm_factor_x, double renorm_factor_y, bool flip,
+                 cufftHandle plan_1d, cufftHandle plan_dim2, double dx,
+                 double dy, double dz, double time);
+
 #endif

include/kernels.h

@@ -17,6 +17,32 @@
 #define KERNELS_H
 #include<stdio.h>
 
+/**
+* @brief        derivative of data
+* @param        input data
+* @param        output data
+* @param        stride of derivative, for (xDim, yDim, zDim) derivative,
+                use stride (1, xDim, xDim*yDim)
+* @param        grid size for simulation
+* @param        dx value for derivative
+* @ingroup      gpu
+*/
+__global__ void derive(double *data, double *out, int stride, int gsize,
+                       double dx);
+
+/**
+* @brief        derivative of data
+* @param        input data
+* @param        output data
+* @param        stride of derivative, for (xDim, yDim, zDim) derivative,
+                use stride (1, xDim, xDim*yDim)
+* @param        grid size for simulation
+* @param        dx value for derivative
+* @ingroup      gpu
+*/
+__global__ void derive(double2 *data, double2 *out, int stride, int gsize,
+                       double dx);
+
 /**
 * @brief	subtraction operation for 2 double2 values
 * @ingroup	gpu
@@ -120,6 +146,15 @@ __host__ __device__ double2 complexMultiply(double2 in1, double2 in2);
 */
 __device__ double2 make_complex(double in, int evolution_type);
 
+/**
+* @brief        copies a double2 value
+* @ingroup      gpu
+* @param        complex input
+* @return       complex output
+*/
+
+__global__ void copy(double2 *in, double2 *out);
+
 /**
 * @brief        Sums the absolute value of two complex arrays
 * @ingroup      gpu
@@ -129,6 +164,44 @@ __device__ double2 make_complex(double in, int evolution_type);
 */
 __global__ void complexAbsSum(double2 *in1, double2 *in2, double *out);
 
+/**
+* @brief        Sums double2* and double2* energies
+* @ingroup      gpu
+* @param        Array 1
+* @param        Array 2
+* @param        Output
+*/
+__global__ void energy_sum(double2 *in1, double2 *in2, double *out);
+
+/**
+* @brief        Sums double* and double2* energies for angular momentum
+* @ingroup      gpu
+* @param        Array 1
+* @param        Array 2
+* @param        Output
+*/
+__global__ void energy_lsum(double *in1, double2 *in2, double *out);
+
+/**
+* @brief        Sums double2* and double2* to an output double2
+* @ingroup      gpu
+* @param        Array 1
+* @param        Array 2
+* @param        Output
+*/
+__global__ void sum(double2 *in1, double2 *in2, double2 *out);
+
+/**
+* @brief        Sums the absolute value of two complex arrays
+* @ingroup      gpu
+* @param        Array 1
+* @param        Array 2
+* @param        Array 3
+* @param        Output
+*/
+__global__ void complexAbsSum(double2 *in1, double2 *in2, double2 *in3,
+                              double *out);
+
 /**
 * @brief        Complex magnitude of a double2 array
 * @ingroup      gpu
@@ -203,11 +276,10 @@ __global__ void cMultPhi(double2* in1, double* in2, double2* out);
 * @param	in2 Evolution operator input
 * @param	out Pass by reference output for multiplication result
 * @param	dt Timestep for evolution
-* @param	mass Atomic species mass
 * @param	gState If performing real (1) or imaginary (0) time evolution
 * @param	gDenConst a constant for evolution
 */
-__global__ void cMultDensity(double2* in1, double2* in2, double2* out, double dt, double mass, int gstate, double gDenConst);
+__global__ void cMultDensity(double2* in1, double2* in2, double2* out, double dt, int gstate, double gDenConst);
 
 /**
 * @brief        Kernel for complex multiplication with nonlinear density term
@@ -221,14 +293,13 @@ __global__ void cMultDensity(double2* in1, double2* in2, double2* out, double dt
 * @param        time
 * @param        element number in AST
 * @param        dt Timestep for evolution
-* @param        mass Atomic species mass
 * @param        gState If performing real (1) or imaginary (0) time evolution
 * @param        gDenConst a constant for evolution
 */
 
 __global__ void cMultDensity_ast(EqnNode_gpu *eqn, double2* in, double2* out,
                                  double dx, double dy, double dz, double time,
-                                 int e_num, double dt, double mass, int gstate,
+                                 int e_num, double dt, int gstate,
                                  double gDenConst);
 
 
@@ -251,6 +322,33 @@ __global__ void pinVortex(double2* in1, double2* in2, double2* out);
 */
 __global__ void vecMult(double2 *in, double *factor, double2 *out);
 
+/**
+* @brief        Complex field Summation
+* @ingroup      gpu
+* @param        in Complex field to be scaled (divided, not multiplied)
+* @param        factor Scaling vector to be used
+* @param        out Pass by reference output for result
+*/
+__global__ void vecSum(double2 *in, double *factor, double2 *out);
+
+/**
+* @brief        field scaling
+* @ingroup      gpu
+* @param        in field to be scaled (divided, not multiplied)
+* @param        factor Scaling vector to be used
+* @param        out Pass by reference output for result
+*/
+__global__ void vecMult(double *in, double *factor, double *out);
+
+/**
+* @brief        field Summation
+* @ingroup      gpu
+* @param        in field to be scaled (divided, not multiplied)
+* @param        factor Scaling vector to be used
+* @param        out Pass by reference output for result
+*/
+__global__ void vecSum(double *in, double *factor, double *out);
+
 /**
 * @brief        performs the l2 normalization of the provided terms
 * @ingroup      gpu
@@ -297,11 +395,30 @@ __global__ void scalarDiv(double* in, double factor, double* out);
 * @brief        Complex field scaling and renormalisation. Used mainly post-FFT.
 * @ingroup      gpu
 * @param        in Complex field to be scaled (multiplied, not divided)
-* @param        factor Scaling factor to be used
+* @param        scaling factor to be used
 * @param        out Pass by reference output for result
 */
 __global__ void scalarMult(double2* in, double factor, double2* out);
 
+/**
+* @brief        field scaling and renormalisation. Used mainly post-FFT.
+* @ingroup      gpu
+* @param        in field to be scaled (multiplied, not divided)
+* @param        scalaing factor to be used
+* @param        out Pass by reference output for result
+*/
+__global__ void scalarMult(double* in, double factor, double* out);
+
+/**
+* @brief        Complex field scaling and renormalisation. Used mainly post-FFT.
+* @ingroup      gpu
+* @param        in Complex field to be scaled (multiplied, not divided)
+* @param        complex scaling factor to be used
+* @param        out Pass by reference output for result
+*/
+__global__ void scalarMult(double2* in, double2 factor, double2* out);
+
+
 /**
 * @brief        Complex field raised to a power
 * @ingroup      gpu
@@ -428,7 +545,11 @@ __device__ double2 im_ast(double val, double dt);
 * @brief        Sets boolean array to 0
 * @ingroup      gpu
 */
-__global__ void zeros(bool *in, bool *out);
+__global__ void zeros(bool *out);
+
+__global__ void zeros(double *out);
+
+__global__ void zeros(double2 *out);
 
 /**
 * @brief        Sets in2 to be equal to in1

include/operators.h

@@ -21,6 +21,15 @@
 #include <unordered_map>
 //#include <boost/math/special_functions.hpp>
 
+// Laplacian functions
+void laplacian(Grid &par, double2 *data, double2* out, int xDim, int yDim,
+              int zDim, double dx, double dy, double dz);
+
+void laplacian(Grid &par, double2 *data, double2* out, int xDim, int yDim,
+              double dx, double dy);
+
+void laplacian(Grid &par, double2 *data, double2* out, int xDim, double dx);
+
 // function to find Bz, the curl of the gauge field
  /**
  * @brief       Finds Bz, the curl of the gauge field
@@ -245,4 +254,5 @@ __global__ void aux_fields(double *V, double *K, double gdt, double dt,
                            double2* EpAx, double2* EpAy, double2* EpAz);
 // Function to generate grid and treads
 void generate_grid(Grid &par);
+
 #endif

julia/calc_energy.jl

@@ -0,0 +1,133 @@
+using DelimitedFiles
+using FFTW
+
+function derive(data, dim, dx)
+    out = Array{Complex,2}(undef,size(data,1), size(data,2))
+    dim2 = dim+1
+    if dim2 == 3
+        dim2 = 1
+    end
+    for i = 1:size(data,dim)
+        for j = 1:size(data,dim2)
+            #println(j,'\t',i)
+            if (dim == 1)
+                if (i == size(data,dim))
+                    out[i,j] = (data[1, j] - data[i, j])/dx
+                else
+                    out[i,j] = (data[i + 1, j] - data[i, j])/dx
+                end
+            else
+                if (i == size(data,dim))
+                    out[j,i] = (data[j,1] - data[j,i])/dx
+                else 
+                    out[j,i] = (data[j,i + 1] - data[j,i])/dx
+                end
+            end
+        end
+    end
+    return out
+end
+
+# We are calculating the energy to check <Psi|H|Psi>
+function calculate_energy(wfc, H_k, H_r, Ax, Ay, g, xDim, yDim, dx, dy)
+    hbar = 1.05457148e-34
+    omega = 0.6
+    omegaX = 6.283
+
+    # Creating momentum and conjugate wavefunctions
+    wfc_k = fft(wfc)
+    wfc_c = conj(wfc)
+
+    # Finding the momentum and real-space energy terms
+    energy_k = wfc_c.*ifft((H_k) .* wfc_k)
+    energy_r = wfc_c.* (H_r).* wfc
+    energy_i = wfc_c.* (0.5*g*abs2.(wfc)).* wfc
+
+    energy_l = wfc_c.*(im*hbar*(Ax.*derive(wfc,1,dx) + Ay.*derive(wfc,2,dy)))
+
+    # Integrating over all space
+    energy_if = 0
+    energy_lf = 0
+    energy_kf = 0
+    energy_rf = 0
+    for i = 1:xDim*yDim
+        energy_if += real(energy_i[i])*dx*dy
+        energy_rf += real(energy_r[i])*dx*dy
+        energy_lf += real(energy_l[i])*dx*dy
+        energy_kf += real(energy_k[i])*dx*dy
+    end 
+
+    println("Kinetic energy:", "\t\t\t", energy_kf)
+    println("Potential energy:", "\t\t", energy_rf)
+    println("Internal energy:", "\t\t", energy_if)
+    println("Angular Momentum energy:", '\t', energy_lf)
+    println("Total energy:", "\t\t\t", energy_kf+energy_if+energy_rf+energy_lf)
+
+end
+
+# Read in param.cfg file
+function calculate(param_file::String, data_dir::String)
+    parameters = Dict()
+
+    for line in readlines(data_dir*param_file)
+        if line != "[Params]"
+            tmp = split(line,"=")
+            parameters[tmp[1]] = tmp[2]
+        end
+    end
+
+    xDim = parse(Int64, parameters["xDim"])
+    yDim = parse(Int64, parameters["yDim"])
+    dx = parse(Float64, parameters["dx"])
+    dy = parse(Float64, parameters["dy"])
+
+    omega = parse(Float64, parameters["omega"])
+    omegaX = parse(Float64, parameters["omegaX"])
+
+    g = parse(Float64, parameters["gDenConst"])
+
+    Ax = readdlm(data_dir*"Ax_0")
+    Ay = readdlm(data_dir*"Ay_0")
+    K = readdlm(data_dir*"K_0")
+    V = readdlm(data_dir*"V_0")
+
+    Ax = reshape(Ax, xDim, yDim)
+    Ay = reshape(Ay, xDim, yDim)
+    K = reshape(K, xDim, yDim)
+    V = reshape(V, xDim, yDim)
+
+    start = 0
+    ende = parse(Int64, parameters["esteps"])
+    endg = parse(Int64, parameters["gsteps"])
+    incr = parse(Int64, parameters["printSteps"])
+
+    # Ground State Evolution
+    println("Starting imaginary time energy calculation")
+    for i = start:incr:endg
+        wfc = readdlm(data_dir*"wfc_0_const_"*string(i)) +
+              readdlm(data_dir*"wfc_0_consti_"*string(i))*im
+        println(data_dir*"wfc_0_const_"*string(i), '\t',
+                data_dir*"wfc_0_consti_"*string(i))
+        wfc = reshape(wfc, xDim, yDim)
+        calculate_energy(wfc, K, V, Ax, Ay, g, xDim, yDim, dx, dy)
+        println()
+    end
+
+    println()
+
+    # Ground State Evolution
+    println("Starting real time energy calculation")
+    for i = start:incr:ende
+        wfc = readdlm(data_dir*"wfc_ev_"*string(i)) +
+              readdlm(data_dir*"wfc_evi_"*string(i))*im
+        println(data_dir*"wfc_0_const_"*string(i), '\t', 
+                data_dir*"wfc_0_consti_"*string(i))
+        wfc = reshape(wfc, xDim, yDim)
+        calculate_energy(wfc, K, V, Ax, Ay, g, xDim, yDim, dx, dy)
+        println()
+    end
+
+
+end
+
+calculate("Params.dat", "../data/")

src/ds.cu

@@ -46,7 +46,7 @@ void generate_plan_other3d(cufftHandle *plan_fft1d, Grid &par, int axis){
 
     cufftResult result;
 
-    // Along first dimension (z)
+    // Along first dimension (x)
     if (axis == 0){
         result = cufftPlan1d(plan_fft1d, xDim, CUFFT_Z2Z, yDim*zDim);
     }
@@ -71,7 +71,7 @@ void generate_plan_other3d(cufftHandle *plan_fft1d, Grid &par, int axis){
 
     }
 
-    // Along third dimension (x)
+    // Along third dimension (z)
     else if (axis == 2){
 
         int batch = xDim*yDim;