tIncreased size of pressure grid, solutions seem wrong - 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 fb42d1dbafff951cfe8ca1a4dc6798e31b59b62b
(DIR) parent eb664c2b1bdfae1ad3348fae95a81c0c5df27b2b
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Mon, 3 Mar 2014 18:49:37 +0100
Increased size of pressure grid, solutions seem wrong
Diffstat:
M ns2dfd.py | 10 +++++-----
M src/Makefile | 2 +-
M src/file_io.c | 8 ++++----
M src/main.c | 23 ++++++++++++-----------
M src/utility.c | 8 ++++----
5 files changed, 26 insertions(+), 25 deletions(-)
---
(DIR) diff --git a/ns2dfd.py b/ns2dfd.py
t@@ -42,7 +42,7 @@ class fluid:
self.ny = numpy.asarray(ny)
self.dx = numpy.asarray(dx)
self.dy = numpy.asarray(dy)
- self.P = numpy.zeros((nx+1, ny+1))
+ self.P = numpy.zeros((nx+2, ny+2))
self.U = numpy.zeros((nx+2, ny+2))
self.V = numpy.zeros((nx+2, ny+2))
t@@ -186,8 +186,8 @@ class fluid:
self.init_grid(dx = self.dx, dy = self.dy,\
nx = self.nx, ny = self.ny)
- for i in range(self.nx+1):
- for j in range(self.ny+1):
+ for i in range(self.nx+2):
+ for j in range(self.ny+2):
self.P[i,j] = \
numpy.fromfile(fh, dtype=numpy.float64, count=1)
t@@ -241,8 +241,8 @@ class fluid:
fh.write(self.nx.astype(numpy.int32))
fh.write(self.ny.astype(numpy.int32))
- for i in range(self.nx+1):
- for j in range(self.ny+1):
+ for i in range(self.nx+2):
+ for j in range(self.ny+2):
fh.write(self.P[i,j].astype(numpy.float64))
for i in range(self.nx+2):
(DIR) diff --git a/src/Makefile b/src/Makefile
t@@ -1,4 +1,4 @@
-CFLAGS=-Wall -pedantic -g
+CFLAGS=-Wall -pedantic -g -O2
LDLIBS=-lm
BIN=../ns2dfd
(DIR) diff --git a/src/file_io.c b/src/file_io.c
t@@ -45,8 +45,8 @@ int read_file(char *path,
allocate_memory(P, U, V, *nx, *ny);
- for (i=0; i<*nx+1; i++) {
- for (j=0; j<*ny+1; j++) {
+ for (i=0; i<*nx+2; i++) {
+ for (j=0; j<*ny+2; j++) {
fread(&tmp, sizeof(double), 1, fp);
(*P)[i][j] = tmp;
}
t@@ -113,8 +113,8 @@ int write_file(char *path,
fwrite(nx, sizeof(int), 1, fp);
fwrite(ny, sizeof(int), 1, fp);
- for (i=0; i<*nx+1; i++) {
- for (j=0; j<*ny+1; j++) {
+ for (i=0; i<*nx+2; i++) {
+ for (j=0; j<*ny+2; j++) {
tmp = (*P)[i][j];
fwrite(&tmp, sizeof(double), 1, fp);
}
(DIR) diff --git a/src/main.c b/src/main.c
t@@ -22,9 +22,10 @@ int main(int argc, char** argv)
double **P, **U, **V;
double **F, **G, **RHS, **RES;
- double dt, res;
+ double dt = 0.0;
+ double res;
long n = 0;
- int it;
+ int it = 0;
int nfile = 0;
double t_file_elapsed = 0.0;
char filename[50];
t@@ -75,10 +76,10 @@ int main(int argc, char** argv)
dot[0] = '\0';
printf("%s: simulation id: `%s`\n", argv[0], simulation_id);
- allocate_matrix(&F, nx+1, ny+1);
- allocate_matrix(&G, nx+1, ny+1);
- allocate_matrix(&RHS, nx+1, ny+1);
- allocate_matrix(&RES, nx+1, ny+1);
+ allocate_matrix(&F, nx+2, ny+2);
+ allocate_matrix(&G, nx+2, ny+2);
+ allocate_matrix(&RHS, nx+2, ny+2);
+ allocate_matrix(&RES, nx+2, ny+2);
while (t <= t_end+dt) {
t@@ -105,7 +106,7 @@ int main(int argc, char** argv)
compute_rhs(RHS, F, G, dt, dx, dy, nx, ny);
it = 0;
- res = 0.0;
+ res = 1.0e9;
while ((it <= itermax) && (res >= epsilon)) {
sor_cycle(P, RHS, omega, nx, ny, dx, dy);
t@@ -124,9 +125,9 @@ 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);
- free_matrix(&RES, nx);
+ free_matrix(&F, nx+2);
+ free_matrix(&G, nx+2);
+ free_matrix(&RHS, nx+2);
+ free_matrix(&RES, nx+2);
return 0;
}
(DIR) diff --git a/src/utility.c b/src/utility.c
t@@ -27,7 +27,7 @@ void free_matrix(double ***M, int nx)
int allocate_memory(double ***P, double ***U, double ***V, int nx, int ny)
{
- if (allocate_matrix(P, nx+1, ny+1) != 0) {
+ if (allocate_matrix(P, nx+2, ny+2) != 0) {
fprintf(stderr, "allocate_memory: Could not allocate memory for P\n");
return 1;
}
t@@ -44,9 +44,9 @@ int allocate_memory(double ***P, double ***U, double ***V, int nx, int ny)
void free_memory(double **P, double **U, double **V, int nx)
{
- free_matrix(&P, nx);
- free_matrix(&U, nx);
- free_matrix(&V, nx);
+ free_matrix(&P, nx+2);
+ free_matrix(&U, nx+2);
+ free_matrix(&V, nx+2);
}
double max_abs_value(double **M, int nx, int ny)