From d7b894f738547e3e385cf06e10b8347c747b0415 Mon Sep 17 00:00:00 2001
From: Volker Hoffmann <volker@cheleb.net>
Date: Thu, 19 Jun 2014 21:25:22 +0200
Subject: upwind advection

---
 advection.cu | 81 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 81 insertions(+)
 create mode 100644 advection.cu

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);
+}
+
+        
-- 
cgit v1.1