ex1: solve

exercise1/wave_1d.c

@@ -43,44 +43,71 @@ void domain_save(int_t step) {
 // Set up our three buffers, fill two with an initial cosine wave,
 // and set the time step.
 void domain_initialize(void) {
-  // BEGIN: T1
-  ;
-  // END: T1
+  // allocate memory
+  for (int i = 0; i < 3; i++)
+    buffers[i] = malloc((N + 2) * sizeof(real_t));
+
+  // apply initial condition
+  real_t x = 0.0f;
+  for (int i = 0; i < N; i++) {
+    U_prv(i) = cos(2 * M_PI * x);
+    U(i) = U_prv(i);
+    U_nxt(i) = 0.0f;
+    x += (real_t)dx / (real_t)N;
+  }
+
+  // set the time-step dt according to CFL
+  dt = dx / c;
 }
 
 // TASK T2:
 // Return the memory to the OS.
-// BEGIN: T2
 void domain_finalize(void) {
-  ;
+  for (int i = 0; i < 3; i++) {
+    free(buffers[i]);
+    buffers[i] = NULL;
+  }
 }
-// END: T2
 
 // TASK: T3
 // Rotate the time step buffers.
-// BEGIN: T3
-;
-// END: T3
+void rotate_buffers(void) {
+  real_t *temp = buffers[0];
+  buffers[0] = buffers[1];
+  buffers[1] = buffers[2];
+  buffers[2] = temp;
+}
 
 // TASK: T4
 // Derive step t+1 from steps t and t-1.
-// BEGIN: T4
-;
-// END: T4
+void time_step(void) {
+  for (int i = 0; i < N; i++) {
+    U_nxt(i) = -U_prv(i) + 2 * U(i) + dt * dt * c * c / (dx * dx) * (U(i - 1) + U(i + 1) - 2 * U(i));
+  }
+}
 
 // TASK: T5
 // Neumann (reflective) boundary condition.
-// BEGIN: T5
-;
-// END: T5
+void boundary(void) {
+  U_prv(-1) = U_prv(1);
+  U(-1) = U(1);
+  U_nxt(-1) = U_nxt(1);
+  U_prv(N) = U_prv(N - 2);
+  U(N) = U(N - 2);
+  U_nxt(N) = U_nxt(N - 2);
+}
 
 // TASK: T6
 // Main time integration.
 void simulate(void) {
-  // BEGIN: T6
   int_t iteration = 0;
-  domain_save(iteration / snapshot_freq);
-  // END: T6
+  while (iteration < max_iteration) {
+    boundary();
+    time_step();
+    rotate_buffers();
+    domain_save(iteration / snapshot_freq);
+    iteration++;
+  }
 }
 
 int main(void) {