-
-
Save dcyang/fb838408bae1e0ce0ba5c5c4744a168a to your computer and use it in GitHub Desktop.
populating and multiplying scalar field grids
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#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 | |
*/ |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#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(√((gridR⋅gridR)*dV)) | |
println() | |
@timev fill_grid_2!(gridR) | |
@time println(√((gridR⋅gridR)*dV)) | |
println() | |
@timev fill_grid_3!(gridR) | |
@time println(√((gridR⋅gridR)*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