Skip to content

Commit

Permalink
Fix for #21 and #24 (#26)
Browse files Browse the repository at this point in the history
* refactoring code to make it easier to compute energy. Test required for 3D vortex states.

* rough draft of energy_calc() function

* modifying energy calculation to use slightly less energy, we should modify for AST's and recommend AST's for use of this calculation in the future.

* adding changes to gauge application functions.

* adding derive function and implementing it for energy calculation

* added laplacian functions

* fixing laplacian functions

* updating unit test to include angular momentum calculation

* adding unit test for derive() function

* adding julia energy script

* adding new energy calculation methods

* adding new energy calc script with derivatives.

* working on energy calculation

* adding *working* energy calculations in gpue

* fixing bug with phase imprinting process. Still largescale error in energy calc

* setting default mask size to 0 to prevent spontaneous seg faults

* modifying julia calculation script to read in Params.dat file

* adding meory check for energy calculation

* fixing memory errors in unit tests
  • Loading branch information
leios authored Apr 5, 2019
1 parent c095e68 commit 07afaaa
Show file tree
Hide file tree
Showing 12 changed files with 883 additions and 289 deletions.
12 changes: 12 additions & 0 deletions include/evolution.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
133 changes: 127 additions & 6 deletions include/kernels.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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);


Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
10 changes: 10 additions & 0 deletions include/operators.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
133 changes: 133 additions & 0 deletions julia/calc_energy.jl
Original file line number Diff line number Diff line change
@@ -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/")
4 changes: 2 additions & 2 deletions src/ds.cu
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
Expand All @@ -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;
Expand Down
Loading

0 comments on commit 07afaaa

Please sign in to comment.