initial commit

6 files changed, 286 insertions, 0 deletions

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..e1a8f01
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,4 @@
+*.dat
+*.o
+*.swp
+kun
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..becbc85
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,31 @@
+PRG=kun
+
+CC=gcc
+CFLAGS=-Wall -Wextra -pedantic -O0
+LFLAGS=
+SRC=io.c main.c
+LIB=-lm
+
+OBJ=$(SRC:.c=.o)
+
+all:$(PRG)
+
+$(PRG): $(OBJ)
+	$(CC) $(LFLAGS) $(CFLAGS) $(OBJ) -o $@ $(LIB)
+
+%.o: %.c
+	$(CC) $(LFLAGS) $(CFLAGS) -c $< -o $@
+
+run: $(PRG)
+	./$(PRG)
+
+clean:
+	rm -f $(PRG) $(OBJ)
+
+cleanall:
+	make clean
+	rm -f *.dat; rm -f *.png
+
+
+.PHONY: all clean cleanall
+
diff --git a/io.c b/io.c
new file mode 100644
index 0000000..522b3dc
--- /dev/null
+++ b/io.c
@@ -0,0 +1,28 @@
+#include <stdio.h>
+#include "io.h"
+
+/*
+ * write out an array
+ */
+void iowrite(char* filename, double* q, int n)
+{
+  FILE *f = fopen(filename, "w");
+  for (int i = 0; i < n; i ++) {
+    fprintf(f, "%f, %f, %f, %f\n", q[4*i], q[4*i+1], q[4*i+2], q[4*i+3]);
+  }
+  fprintf(f, "\n");
+  fclose(f);
+}
+
+/*
+ * write out an array
+ */
+void iocoord(char* filename, double* x, int n)
+{
+  FILE *f = fopen(filename, "w");
+  for (int i = 0; i < n; i ++) {
+    fprintf(f, "%f\n", x[i]);
+  }
+  fprintf(f, "\n");
+  fclose(f);
+}
diff --git a/io.h b/io.h
new file mode 100644
index 0000000..5f20c03
--- /dev/null
+++ b/io.h
@@ -0,0 +1,8 @@
+#ifndef __IOXXXX_INCLUDED__
+#define __IOXXXX_INCLUDED__
+
+void iowrite(char* filename, double* q, int n);
+
+void iocoord(char* filename, double* x, int n);
+
+#endif
diff --git a/main.c b/main.c
new file mode 100644
index 0000000..7411fc0
--- /dev/null
+++ b/main.c
@@ -0,0 +1,165 @@
+#include <stdio.h>
+#include <stdlib.h>
+
+#include <math.h>
+#include <string.h>
+
+#include "io.h"
+
+/*
+ * main.c : acoustic equations
+ */
+int
+main (void)
+{
+
+    int n = 500;
+    double xl = 0.0;
+    double xr = 1.0;
+    double dx = (xr - xl) / n;
+
+    double dt = 5.0e-4;
+
+    double* x = (double*) calloc (n, sizeof (double));
+    double* u = (double*) calloc (4 * n, sizeof (double));
+    double* b = (double*) calloc (n, sizeof (double));
+
+    double* u12 = (double*) calloc (4 * (n + 1), sizeof (double));
+    double* f12 = (double*) calloc (4 * (n + 1), sizeof (double));
+
+    double* unew = (double*) calloc (4 * n, sizeof (double));
+
+    // initial conditions
+
+    for (int i = 0; i < n; i ++) {
+        x[i] = 0.5 * dx + dx * i;
+    }
+
+    iocoord("coord.dat", x, n);
+
+    double hbar    = 0.2;
+    double xbar    = 0.075;
+    double x0      = 0.4;
+
+    for (int i = 0; i < n; i ++) {
+        if (fabs (x[i] - x0) < xbar) {
+            u[4*i] = hbar * sqrt (1.0 - pow ((x[i] - x0) / xbar, 2));
+        }
+    }
+
+    for (int i = 0; i < n; i ++) {
+        if (x[i] < 0.6) {
+            u[4*i+2] = 1.0;
+        } else {
+            u[4*i+2] = 4.0;
+        }
+        u[4*i+3] = 1.0;
+    }
+
+    iowrite("q00.dat", u, n);
+
+    //
+    //
+    //
+    // -- initialized --
+    //
+    //
+    //
+
+    int m = 0;
+    while (1) {
+
+        printf (">>> step %i to be taken. press any key to continue.", m);
+	printf (">>> time passed is %f s.", m * dt);
+        getchar ();
+
+        u12[0] = u[0];
+        u12[1] = u[1];
+        u12[2] = u[2];
+        u12[3] = u[3];
+        u12[4*n] = u[4*(n-1)];
+        u12[4*n+1] = u[4*(n-1)+1];
+        u12[4*n+2] = u[4*(n-1)+2];
+        u12[4*n+3] = u[4*(n-1)+3];
+
+        for (int i = 1; i < n; i ++) {
+            double pl = u[4*(i-1)];
+            double ul = u[4*(i-1)+1];
+            double rl = u[4*(i-1)+2];
+            double kl = u[4*(i-1)+3];
+
+            double pr = u[4*i];
+            double ur = u[4*i+1];
+            double rr = u[4*i+2];
+            double kr = u[4*i+3];
+
+            double hat_a = 0.5 * (1.0 / rl + 1.0 / rr);
+            double hat_k = 0.5 * (kl + kr);
+
+            double dq1 = pr - pl;
+            double dq2 = ur - ul;
+
+            f12[4*i]   = 0.5 * hat_k * dq2;
+            f12[4*i+1] = 0.5 * hat_a * dq1;
+            f12[4*i+2] = 0.0;
+            f12[4*i+3] = 0.0;
+
+            for (int j = 0; j < 4; j ++) {
+                u12[4*i+j] = 0.5 * (u[4*(i-1)+j] + u[4*i+j]) - dt * f12[4*i+j] / dx;
+            }
+        }
+
+        memset (f12, 0.0, (4*(n+1)) * sizeof (double));
+
+	iowrite("q12.dat", u12, n+1);
+
+        for (int i = 1; i < n + 1; i ++) {
+            double pl = u12[4*(i-1)];
+            double ul = u12[4*(i-1)+1];
+            double rl = u12[4*(i-1)+2];
+            double kl = u12[4*(i-1)+3];
+
+            double pr = u12[4*i];
+            double ur = u12[4*i+1];
+            double rr = u12[4*i+2];
+            double kr = u12[4*i+3];
+
+            double star_a = 0.5 * (1.0 / rl + 1.0 / rr);
+            double star_k = 0.5 * (kr + kl);
+
+            double dq1 = pr - pl;
+            double dq2 = ur - ul;
+
+            f12[4*i]   = 0.5 * star_k * dq2;
+            f12[4*i+1] = 0.5 * star_a * dq1;
+            f12[4*i+2] = 0.0;
+            f12[4*i+3] = 0.0;
+
+            for (int j = 0; j < 4; j ++) {
+	      unew[4*(i-1)+j] = 0.5 * (u12[4*(i-1)+j] + u12[4*i+j]) - dt * f12[4*i+j] / dx;
+	      if (j == 2)
+		unew[4*(i-1)+j] = u[4*(i-1)+j];
+            }
+        }
+
+	iowrite("qne.dat", unew, n);
+
+        memcpy (u,    unew, 4*n * sizeof (double));
+        memset (u12,  0.0, (4*(n+1)) * sizeof (double));
+        memset (f12,  0.0, (4*(n+1)) * sizeof (double));
+        memset (unew, 0.0,  4*n * sizeof (double));
+
+        m = m + 1;
+
+    }
+
+
+    free (unew);
+    free (f12);
+    free (u12);
+    free (b);
+    free (u);
+    free (x);
+
+    return 0;
+}
diff --git a/plot.py b/plot.py
new file mode 100755
index 0000000..80acee3
--- /dev/null
+++ b/plot.py
@@ -0,0 +1,50 @@
+#!/usr/bin/env python
+
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.style as style
+
+style.use('classic')
+
+x00 = np.genfromtxt('coord.dat', delimiter=',')
+dx  = x00[1] - x00[0]
+x12 = x00 - 0.5 * dx
+x12 = np.append(x12, x12[len(x12)-1] + dx)
+
+q00 = np.genfromtxt('q00.dat', delimiter=',')
+q12 = np.genfromtxt('q12.dat', delimiter=',')
+
+qne = np.genfromtxt('qne.dat', delimiter=',')
+
+fig, (ax1, ax2, ax3, ax4) = plt.subplots (4,1,sharex=True)
+
+ax1.set_title('pressure')
+ax1.plot (x00, q00[:,0], 'ko-', fillstyle='none', label='init')
+#ax1.plot (x12, q12[:,0], 'ro-', fillstyle='none', label='t+1/2')
+ax1.plot (x00, qne[:,0], 'g^-', fillstyle='none', label='updated')
+
+ax2.set_title('velocity')
+ax2.plot (x00, q00[:,1], 'ko-', fillstyle='none', label='init')
+#ax2.plot (x12, q12[:,1], 'ro-', fillstyle='none', label='t+1/2')
+ax2.plot (x00, qne[:,1], 'g^-', fillstyle='none', label='updated')
+
+ax3.set_title('density')
+ax3.plot (x00, q00[:,2], 'ko-', fillstyle='none', label='init')
+#ax3.plot (x12, q12[:,2], 'ro-', fillstyle='none', label='t+1/2')
+ax3.plot (x00, qne[:,2], 'g^-', fillstyle='none', label='updated')
+
+ax4.set_title('bulk modulus')
+ax4.plot (x00, q00[:,3], 'ko-', fillstyle='none', label='init')
+#ax4.plot (x12, q12[:,3], 'ro-', fillstyle='none', label='t+1/2')
+ax4.plot (x00, qne[:,3], 'g^-', fillstyle='none', label='updated')
+
+ax1.legend()
+#ax2.legend()
+#ax3.legend()
+#ax4.legend()
+
+ax4.set_xlabel('x (m)')
+
+plt.tight_layout()
+plt.savefig ('res104.png', dpi=200, format='png')
+