summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorVolker Hoffmann <volker@cheleb.net>2014-06-19 21:25:22 +0200
committerVolker Hoffmann <volker@cheleb.net>2014-06-19 21:25:22 +0200
commitd7b894f738547e3e385cf06e10b8347c747b0415 (patch)
tree0e70c76b506fb77c6701e74f1d55405b46e038a1
parente0e74374c61a60c24001a7593ffe2a68feb2da80 (diff)
upwind advection
-rw-r--r--advection.cu81
1 files changed, 81 insertions, 0 deletions
diff --git a/advection.cu b/advection.cu
new file mode 100644
index 0000000..7851ecd
--- /dev/null
+++ b/advection.cu
@@ -0,0 +1,81 @@
+// advection.cu
+#include <stdio.h>
+#include <assert.h>
+#include <cuda.h>
+
+__global__ void advect(float *uold, float *unew, float cs, float dx, float dt, int N)
+{
+ int idx = blockIdx.x * blockDim.x + threadIdx.x;
+ if (idx==0) {
+ unew[idx] = uold[idx] - dt * cs * ( uold[idx] - uold[N-1] ) / dx;
+ } else {
+ unew[idx] = uold[idx] - dt * cs * ( uold[idx] - uold[idx-1] ) / dx;
+ }
+ printf("(%i,%.3f,%.3f)\n", idx, uold[idx], unew[idx]);
+}
+
+int main(void)
+{
+ float *u_h; // pointer to host memory
+ float *uold_d, *unew_d; // pointers to device memory
+
+ float cs = 1.0; // bulk velocity
+ float dt = 0.1; // timestep
+ float dx = 0.2; // space step
+
+ int N = 32; // cells
+
+ int i; // counters
+
+ // check CFL condition
+ assert(cs*dt < dx);
+
+ // allocate host memory
+ u_h = (float *)malloc(sizeof(float)*N);
+
+ // allocate device memory
+ cudaMalloc((void **) &uold_d, sizeof(float)*N);
+ cudaMalloc((void **) &unew_d, sizeof(float)*N);
+
+ // fill initial array on host
+ for (i=0; i<N; i++) {
+ u_h[i] = 0.f;
+ }
+
+ // put a square in
+ u_h[N/2-2] = 1.f;
+ u_h[N/2-1] = 1.f;
+ u_h[N/2] = 1.f;
+ u_h[N/2+1] = 1.f;
+ u_h[N/2+2] = 1.f;
+
+ // copy to device
+ cudaMemcpy(uold_d, u_h, sizeof(float)*N, cudaMemcpyHostToDevice);
+
+ // launch configuration
+ int blockSize = N;
+ int nBlocks = 1;
+
+ // run simulation
+ for(i=0; i<10; i++) {
+ printf("*** %i\n", i);
+ advect <<< nBlocks, blockSize >>> (uold_d, unew_d, cs, dx, dt, N);
+ cudaDeviceSynchronize();
+ printf("\n");
+ cudaMemcpy(uold_d, unew_d, sizeof(float)*N, cudaMemcpyDeviceToDevice);
+ }
+ printf("\n");
+
+ // copy back
+ cudaMemcpy(u_h, unew_d, sizeof(float)*N, cudaMemcpyDeviceToHost);
+
+ // dump output
+ for (i=0; i<N; i++) {
+ printf("%.3f\n", u_h[i]);
+ }
+
+ // cleanup
+ free(u_h); cudaFree(unew_d); cudaFree(uold_d);
+}
+
+