summaryrefslogtreecommitdiffstats
path: root/hydro1d.cu
blob: 357af9eeb3f909d974c112d26a5836780034d11c (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
#include <cuda.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define NSTEPS 128              // number of time steps
#define NX 12                   // number of cells in computational domain (EXCLUDING BOUNDARY) (0..11 / 0..NX-1)
#define NTHREADS_PER_BLOCK 4    // number of threads per block (0..3 / 0..NTHREADS_PER_BLOCK-1)
#define NBOUNDARY 1             // boundary cells added on each side of the domain





__device__ void syncBoundaries() {}





// all arrays/pointers passed are nthreads (*float) wide
template <int nthreads>
__device__ void advectStep(float *ui,float *qold, float *qnew,
                           int itx, int dx, int dt)
{
    __shared__ flux_s[nthreads];
    
    // Compute Flux
    if(ui[itx]>0.0) {
        flux[itx] = q[itx] * ui[itx];
    } else {
        flux[itx] = q[itx+1] * ui[itx];
    }
    __syncthreads();

    // Update Me
    qnew[itx] = qold[itx] - dt / dx * ( flux[itx] - flux[itx-1] );
}

template <int nthreads, int ntotal>
__global__ void hydroStep(float *rho_old_d, float *rhou_old_d, 
                          float *rho_new_d, float *rhou_new_d,
                          int dt, int dx, int itx)
{
    // global and local indices
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    int itx = threadIdx.x;
    
    // interface  velocity
    __shared__ float ui[nthreads]
    ui[itx] = 0.5 * ( rhou[itx]/rho[itx] - rhou[itx+1]/rho[itx+1] );

    // first hydro half-steps
    advectStep <nthreads> (ui, rho_old_s, rhou_new_s, itx, dx, dt);
    advectStep <nthreads> (ui, rhou_old_s, rhou_new_s, itx, dx, dt);

    // update old to new before we do the next
    ...

    // new do the pressure half step
    
    










        // update ghost cells (periodic BCs)
        rho_new_d[0] = rho_new[ntotal-1]
        rho_new_d[ntotal] = rho_new[1]
        rhou_new_d[0] = rhoi_new[ntotal-1]
        rhou_new_d[ntotal] = rhoi_new[0]
        __syncthreads();
    }

    // update old/new
    rho_old[idx] = rho_new[idx]
    rhou_old[idx] = rhou_new[idx]

    // compute pressure
    p[idx] = (gamma - 1.0) * rho_new_d * eint

    // hydro, second half-step
    if(idx>0 && idx<ntotal=1) {
        // source terms
        
        p[idx] = gamma - 1.0 ) * rho_new
        
        
        

        //printf("%i %.4f %.4f\n", idx, umid, rho_old_d[idx]);
    }
}

// rho_d and rhou_d are ncells wide
template <nthreads>
__global__ void hydroLoop(float *rho_d, float *rhou_d,
                          float cfl, int dx)
{
    // global and local indices
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    int itx = threadIdx.x;

    // prepare shared memory
    __shared__ float rho_new_s[nthreads];
    __shared__ float rho_old_s[nthreads];
    __shared__ float rhou_new_s[nthreads];
    __shared__ float rhou_old_s[nthreads];

    rho_old_s[itx] = rho_d[idx]; rhou_old_s[itx] = rho_d[idx];
    rho_new_s[itx] = 0.f; rhou_new_s[itx] = 0.f;
    
    // other variables
    int dt;
    int dx;
    
    // main loop, be liberal with barriers
    while(t < tmax) {
        dt = 1.0;
        t += dt;
        hydroStep <nthreads> (rho_old_s, rhou_old_s, rho_new_s, rhou_new_s, dt, dx, itx);
        __syncthreads();
        updateBoundaries();
        __syncthreads();
        copyNewToold();
        __syncthreads();
    }

    // quite kernel, now we can write back into main memory, write outputs, etc
    // when that's done, we can terminate (if we have reached tend)
    // or jump back on the hydro loop here (if we have not reached the final time)

}

__global__ void dumpMem(float *x_d, float *rho_old_d, float *rhou_old_d, float *rho_new_d, float *rhou_new_d)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    //printf("%i %.4f %.4f %.4f %.4f %.4f\n", idx, x_d[idx], rho_old_d[idx], rhou_old_d[idx], rho_new_d[idx], rhou_new_d[idx]);
    printf("%i %.4f %.4f %.4f\n", idx, x_d[idx], rho_old_d[idx], rhou_old_d[idx]);
}

int main()
{
    float *rho_h, *rhou_h;          // density and momentum on the host
    float *rho_old_d, *rhou_old_d;  // density and momentum on the device
    float *rho_new_d, *rhou_new_d;  // density and momentum on the device

    float *x_h, *x_d;     // cell center grid

    float xlo = -5.5;   // lower grid limit
    float xhi =  5.5;   // upper grid limit
    float dx;           // grid spacing

    // index transformations
    const int NTOTAL = NX + 2 * NBOUNDARY;

    // preamble
    cudaSetDevice(0);

    // set up grid
    dx = (xhi - xlo) / float(NX-1);
    printf("%.4f %.4f\n\n", dx, float(NX));
    x_h = (float *) malloc(sizeof(float)*NTOTAL);
    cudaMalloc((void**) &x_d, sizeof(float)*NTOTAL);

    // Computational Domain    
    printf("X-Grid:\n");
    for(int ii=NBOUNDARY; ii<NTOTAL; ii++) {
        x_h[ii] = xlo + (ii-NBOUNDARY) * dx;
    }
    // Boundary Cells (Only Works For One Boundary Cell!)
    x_h[0] = xlo - dx;
    x_h[NTOTAL-1] = xhi + dx;
    // Debug Print
    for(int ii=0; ii<NTOTAL; ii++) {
        printf("%.4f ", x_h[ii]);
    }
    printf("\n\n");
    cudaMemcpy(x_d, x_h, sizeof(float)*NTOTAL, cudaMemcpyHostToDevice);

    // initial conditions
    rho_h = (float*) malloc(sizeof(float)*NTOTAL);
    rhou_h = (float*) malloc(sizeof(float)*NTOTAL);
    printf("ICs\n");
    for(int ii=0; ii<NTOTAL; ii++) {
        rhou_h[ii] = 0.f;
        rho_h[ii] = 1.f + exp(-(x_h[ii]*x_h[ii]) / (1.0*1.0));
        printf("%.4f ", rho_h[ii]);
    }
    printf("\n\n");
    cudaMalloc((void**) &rho_old_d, sizeof(float)*NTOTAL);
    cudaMalloc((void**) &rhou_old_d, sizeof(float)*NTOTAL);
    cudaMemcpy(rho_old_d, rho_h, sizeof(float)*NTOTAL, cudaMemcpyHostToDevice);
    cudaMemcpy(rhou_old_d, rhou_h, sizeof(float)*NTOTAL, cudaMemcpyHostToDevice);

    // prepare new step
    cudaMemset(rho_new_d, 0, sizeof(float)*NTOTAL);
    cudaMemset(rhou_new_d, 0, sizeof(float)*NTOTAL);

    // call cuda kernel to check on memory contents here!!!!
    dumpMem <<< 1, 14 >>> (x_d, rho_old_d, rhou_old_d, rho_new_d, rhou_new_d);
    cudaDeviceSynchronize();
    printf("\n");

    // compute gpu parameters
    int nblocks = NX / NTHREADS_PER_BLOCK + NX % NTHREADS_PER_BLOCK;

    // hydro step
    hydroStep < 32, NTOTAL > <<< 1, 14 >>> (rho_old_d, rhou_old_d, rho_new_d, rhou_new_d);
    cudaDeviceSynchronize();
    printf("\n");

    // dump cpu memory
    // dumpMem <<< 1, 14 >>> (x_d, rho_old_d, rhou_old_d, rho_new_d, rhou_new_d); 
    
    // for(int istep=0; istep<NSTEPS; istep++) {
        // run our physics
        // <32> [32x32] is the biggest block size we can have!
        // alt: diffuse <32> <<< (1,1,1), (32,32,1) >>> (...)
        // diffuse <32,1> <<< nBlocks, nThreads  >>> (uold_d, unew_d, Ncells);
        // set uold to unew before the next step
        // updateMem <<< nBlocks, nThreads >>> (uold_d, unew_d, Ncells);
    // }
        
    // copy back into host memory
    // cudaMemcpy(u_h, unew_d, sizeof(float)*Ncells*Ncells, cudaMemcpyDeviceToHost);
  
    // output array
    // for(int i=0; i<Ncells; i++) {
    //    for(int j=0; j<Ncells; j++) {
    //        printf("%i %i %.16f \n", i, j, u_h[j*Ncells+i]);
    //    }
    // }

    // free memory
    free(x_h); free(rho_h); free(rhou_h);
    cudaFree(x_d);
    cudaFree(rho_old_d); cudaFree(rhou_old_d);
    cudaFree(rho_new_d); cudaFree(rhou_new_d);

    return 0;
}