tWorking on Chorin's projection method: horizontal projection done - ns2dfd - 2D finite difference Navier Stokes solver for fluid dynamics
(HTM) git clone git://src.adamsgaard.dk/ns2dfd
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
(DIR) commit ccc17b2737056c0104d110fa098e7f330be5c919
(DIR) parent dd1cbc3067f419e3cc5263c2c892d4db0dd63a6c
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Mon, 3 Mar 2014 14:37:58 +0100
Working on Chorin's projection method: horizontal projection done
Diffstat:
M src/Makefile | 4 ++--
M src/boundary.c | 78 ++++++++++++++++++++++++++-----
M src/main.c | 13 +++++++++++++
M src/utility.h | 3 +++
4 files changed, 85 insertions(+), 13 deletions(-)
---
(DIR) diff --git a/src/Makefile b/src/Makefile
t@@ -3,10 +3,10 @@ LDLIBS=-lm
BIN=../ns2dfd
-$(BIN): main.o file_io.o utility.o boundary.o
+$(BIN): main.o file_io.o utility.o boundary.o solution.o
$(CC) $(CFLAGS) $(LDLIBS) $^ -o $@
-main.o: file_io.h utility.h boundary.h
+main.o: file_io.h utility.h boundary.h solution.h
file_io.o: utility.h
clean:
(DIR) diff --git a/src/boundary.c b/src/boundary.c
t@@ -1,29 +1,85 @@
-
-#define BC_NO_SLIP 1
-#define BC_FREE_SLIP 2
+#define BC_FREE_SLIP 1
+#define BC_NO_SLIP 2
#define BC_OUTFLOW 3
#define BC_PERIODIC 4
/* Set boundary values of u, v and p according to the boundary conditions.
* The BC flags are as follows:
- * 1: No-slip condition
- * 2: Free-slip condition
- * 3: Outflow condition
- * 4: Periodic condition
+ * 1: Free-slip condition
+ * 2: No-slip condition
+ * 3: Outflow condition
+ * 4: Periodic condition
*/
void set_boundary_conditions(int w_left, int w_right, int w_top, int w_bottom,
double **P, double **U, double **V, int nx, int ny)
{
int i, j;
- if (w_left == BC_NO_SLIP) {
- i = 0;
- for (j=1; j<ny+1; j++) {
- U[i][j] = 0.0;
+ /* vertical boundaries (left and right) */
+ for (j=1; j<ny+1; j++) {
+ if (w_left == BC_FREE_SLIP) {
+ U[0][j] = 0.0;
+ V[0][j] = V[1][j];
+ } else if (w_left == BC_NO_SLIP) {
+ U[0][j] = 0.0;
+ V[0][j] = -1.0*V[1][j];
+ } else if (w_left == BC_OUTFLOW) {
+ U[0][j] = U[1][j];
+ V[0][j] = V[1][j];
+ } else if (w_left == BC_PERIODIC) {
+ U[0][j] = U[nx-1][j];
+ V[0][j] = V[nx-1][j];
+ V[1][j] = V[nx][j];
+ P[1][j] = P[nx][j];
+ }
+ if (w_right == BC_FREE_SLIP) {
+ U[nx][j] = 0.0;
+ V[nx+1][j] = V[nx][j];
+ } else if (w_right == BC_NO_SLIP) {
+ U[nx][j] = 0.0;
+ V[nx+1][j] = -1.0*V[nx][j];
+ } else if (w_right == BC_OUTFLOW) {
+ U[nx][j] = U[nx-1][j];
+ V[nx+1][j] = V[nx][j];
+ } else if (w_right == BC_PERIODIC) {
+ U[nx][j] = U[1][j];
+ V[nx+1][j] = V[2][j];
}
}
+ /* horizontal boundaries (bottom and top) */
+ for (i=1; i<nx+1; i++) {
+
+ if (w_bottom == BC_FREE_SLIP) {
+ U[i][0] = U[i][1];
+ V[i][0] = 0.0;
+ } else if (w_bottom == BC_NO_SLIP) {
+ U[i][0] = -1.0*U[i][1];
+ V[i][0] = 0.0;
+ } else if (w_bottom == BC_OUTFLOW) {
+ U[i][0] = U[i][1];
+ V[i][0] = V[i][1];
+ } else if (w_bottom == BC_PERIODIC) {
+ U[i][0] = U[i][ny-1];
+ U[i][1] = U[i][nx-1];
+ V[i][0] = V[i][ny-1];
+ P[i][1] = P[i][ny];
+ }
+ if (w_top == BC_FREE_SLIP) {
+ U[i][ny+1] = U[i][ny];
+ V[i][ny] = 0.0;
+ } else if (w_top == BC_NO_SLIP) {
+ U[i][ny+1] = -1.0*U[i][ny];
+ V[i][ny] = 0.0;
+ } else if (w_top == BC_OUTFLOW) {
+ U[i][ny+1] = U[i][ny];
+ V[i][ny] = V[i][ny-1];
+ } else if (w_top == BC_PERIODIC) {
+ U[i][ny+1] = U[i][2];
+ V[i][ny] = V[i][1];
+ }
+ }
}
(DIR) diff --git a/src/main.c b/src/main.c
t@@ -6,6 +6,7 @@
#include "file_io.h"
#include "utility.h"
#include "boundary.h"
+#include "solution.h"
#define VERSION 0.1
t@@ -19,6 +20,7 @@ int main(int argc, char** argv)
double dx, dy;
int nx, ny;
double **P, **U, **V;
+ double **F, **G, **RHS;
double dt;
long n = 0;
t@@ -72,6 +74,10 @@ int main(int argc, char** argv)
dot[0] = '\0';
printf("%s: simulation id: `%s`\n", argv[0], simulation_id);
+ allocate_matrix(&F, nx, ny);
+ allocate_matrix(&G, nx, ny);
+ allocate_matrix(&RHS, nx, ny);
+
while (t <= t_end+dt) {
dt = select_time_step(tau, re, dx, dy, nx, ny, U, V);
t@@ -92,6 +98,10 @@ int main(int argc, char** argv)
set_boundary_conditions(w_left, w_right, w_top, w_bottom, P, U, V,
nx, ny);
+ compute_f_g(F, G, U, V, P, re, nx, ny, dx, dy, gx, gy, dt, gamma);
+
+
+
t += dt;
n++;
t_file_elapsed += dt;
t@@ -99,5 +109,8 @@ int main(int argc, char** argv)
puts("\n");
free_memory(P, U, V, nx);
+ free_matrix(&F, nx);
+ free_matrix(&G, nx);
+ free_matrix(&RHS, nx);
return 0;
}
(DIR) diff --git a/src/utility.h b/src/utility.h
t@@ -1,6 +1,9 @@
#ifndef UTILITY_H_
#define UTILITY_H_
+int allocate_matrix(double ***M, int nx, int ny);
+void free_matrix(double ***M, int nx);
+
int allocate_memory(double ***P, double ***U, double ***V, int nx, int ny);
void free_memory(double **P, double **U, double **V, int nx);