Skip to content

Instantly share code, notes, and snippets.

@dcyang
Last active February 28, 2017 15:40
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save dcyang/fb838408bae1e0ce0ba5c5c4744a168a to your computer and use it in GitHub Desktop.
Save dcyang/fb838408bae1e0ce0ba5c5c4744a168a to your computer and use it in GitHub Desktop.
populating and multiplying scalar field grids
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
const double x0 = -10.0;
const double y0 = -10.0;
const double z0 = -10.0;
const int xRes = 400;
const int yRes = 400;
const int zRes = 400;
double dx, dV;
//int num_points;
void fill_grid_1(double ***gridR);
double R_5_0(double r);
int main() {
int i, j, k;
double ***gridR, s;
dx = (fabs(x0)+fabs(x0))/(double)xRes;
dV = pow((fabs(x0)+fabs(x0))/(double)xRes, 3);
//num_points = (xRes+1)*(yRes+1)*(zRes+1);
gridR = (double ***)malloc(sizeof(double **)*(zRes+1));
for (k = 0; k <= zRes; ++k) {
gridR[k] = (double **)malloc(sizeof(double *)*(yRes+1));
for (j = 0; j <= yRes; ++j)
gridR[k][j] = (double *)malloc(sizeof(double)*(xRes+1));
}
fill_grid_1(gridR);
s = 0.0;
for (k = 0; k <= zRes; ++k)
for (j = 0; j <= yRes; ++j)
for (i = 0; i <= xRes; ++i)
s += gridR[i][j][k]*gridR[i][j][k];
printf("integral = %.10lf\n", sqrt(s*dV));
for (k = 0; k <= zRes; ++k) {
for (j = 0; j <= yRes; ++j) free(gridR[k][j]);
free(gridR[k]);
}
free(gridR);
return 0;
}
void fill_grid_1(double ***gridR) {
int i, j, k;
double rVec[3];
for (k = 0; k <= zRes; ++k) {
for (j = 0; j <= yRes; ++j) {
for (i = 0; i <= xRes; ++i) {
rVec[0] = x0 + dx*(double)i;
rVec[1] = y0 + dx*(double)j;
rVec[2] = z0 + dx*(double)k;
gridR[i][j][k] = R_5_0(sqrt(rVec[0]*rVec[0] + rVec[1]*rVec[1] + rVec[2]*rVec[2]));
}
}
}
}
double R_5_0(double r) {
double rho5 = (r+r)/5.0;
return exp(-0.5*rho5) * (120.0 - 240.0*rho5 + 120.0*rho5*rho5 - 20.0*rho5*rho5*rho5 + rho5*rho5*rho5*rho5) / (300.0 * sqrt(5.0));
}
/* Result ...
% time ./test-sfield_c
integral = 0.8816812671
./test-sfield_c 6.58s user 0.09s system 99% cpu 6.670 total
*/
#http://en.citizendium.org/wiki/Hydrogen-like_atom#List_of_radial_functions
const x0 = y0 = z0 = -10.0::Float64
const xRes = yRes = zRes = 400::Int
const dx = 2*abs(x0) / xRes
const dV = dx*dx*dx
#const num_points = (xRes+1)*(yRes+1)*(zRes+1)
function main()
gridR = zeros(Float64, (xRes+1,yRes+1,zRes+1))
#fill_grid_3!(gridR) # warm up in case needed
@timev fill_grid_1!(gridR)
@time println(((gridRgridR)*dV))
println()
@timev fill_grid_2!(gridR)
@time println(((gridRgridR)*dV))
println()
@timev fill_grid_3!(gridR)
@time println(((gridRgridR)*dV))
#@timev println(norm(gridR[:])*√dV)
end
# radial wave function (5,0) for hydrogenic atom
function R_5_0(r::Float64)
ρ5 = 2r/5
return exp(-ρ5/2) * (120 - 240*ρ5 + 120*ρ5*ρ5 - 20*ρ5*ρ5*ρ5 + ρ5*ρ5*ρ5*ρ5) / (300 * 5)
end
function fill_grid_1!(gridR::Array{Float64,3})
rVec = zeros(Float64, 3)
s = 1
for k in 0:zRes
for j in 0:yRes
for i in 0:xRes
rVec[1] = x0 + i*dx
rVec[2] = y0 + j*dx
rVec[3] = z0 + k*dx
gridR[s] = R_5_0(norm(rVec))
#gridR[i+1,j+1,k+1] = R_5_0(norm(rVec))
s += 1
end
end
end
nothing
end
function fill_grid_2!(gridR::Array{Float64,3})
s = 1
for k in 0:zRes
for j in 0:yRes
for i in 0:xRes
rVec = [x0 + i*dx, y0 + j*dx, z0 + k*dx]
gridR[s] = R_5_0(norm(rVec))
#gridR[i+1,j+1,k+1] = R_5_0(norm(rVec))
s += 1
end
end
end
nothing
end
function fill_grid_3!(gridR::Array{Float64,3})
gridR[:] = [R_5_0(norm([x0 + i*dx, y0 + j*dx, z0 + k*dx])) for k in 0:zRes for j in 0:yRes for i in 0:xRes]
nothing
end
main()
#=
Result ...
32.136163 seconds (387.05 M allocations: 5.770 GiB, 0.42% gc time)
elapsed time (ns): 32136162747
gc time (ns): 135494323
bytes allocated: 6195387024
pool allocs: 387049114
GC pauses: 270
0.8816812670740048
0.076864 seconds (19.94 k allocations: 1.021 MiB)
14.145461 seconds (516.01 M allocations: 13.457 GiB, 3.03% gc time)
elapsed time (ns): 14145460586
gc time (ns): 428191150
bytes allocated: 14448972944
pool allocs: 516011403
GC pauses: 630
0.8816812670740048
0.028958 seconds (6 allocations: 128 bytes)
88.498985 seconds (1.55 G allocations: 53.879 GiB, 4.99% gc time)
elapsed time (ns): 88498985076
gc time (ns): 4417725888
bytes allocated: 57851827744
pool allocs: 1549317246
non-pool GC allocs:9
malloc() calls: 1
realloc() calls: 12
GC pauses: 2522
full collections: 7
0.8816812670740048
0.029477 seconds (6 allocations: 128 bytes)
=#
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment