tfinished solver, still to be tested - 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 eb664c2b1bdfae1ad3348fae95a81c0c5df27b2b
(DIR) parent 108a9b876d69391731de42da1ab178384198f345
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Mon, 3 Mar 2014 18:37:05 +0100
finished solver, still to be tested
Diffstat:
M src/main.c | 30 +++++++++++++++++++++++-------
M src/solution.c | 129 ++++++++++++++++++++++++++++++-
M src/solution.h | 11 +++++++++++
3 files changed, 161 insertions(+), 9 deletions(-)
---
(DIR) diff --git a/src/main.c b/src/main.c
t@@ -20,10 +20,11 @@ int main(int argc, char** argv)
double dx, dy;
int nx, ny;
double **P, **U, **V;
- double **F, **G, **RHS;
+ double **F, **G, **RHS, **RES;
- double dt;
+ double dt, res;
long n = 0;
+ int it;
int nfile = 0;
double t_file_elapsed = 0.0;
char filename[50];
t@@ -74,9 +75,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);
+ 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);
while (t <= t_end+dt) {
t@@ -92,15 +94,28 @@ int main(int argc, char** argv)
nfile++;
}
- printf("\rt = %f/%.2f s, dt = %f s, last output: %d ",
- t, t_end, dt, nfile-1);
+ printf("\rt = %f/%.2f s, dt = %f s, it = %d, last output: %d ",
+ t, t_end, dt, it, nfile-1);
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);
+ compute_rhs(RHS, F, G, dt, dx, dy, nx, ny);
+ it = 0;
+ res = 0.0;
+ while ((it <= itermax) && (res >= epsilon)) {
+
+ sor_cycle(P, RHS, omega, nx, ny, dx, dy);
+
+ res = max_residual_norm(RES, P, RHS, nx, ny, dx, dy);
+
+ it++;
+ }
+
+ correct_velocities(U, V, F, G, P, nx, ny, dt, dx, dy);
t += dt;
n++;
t@@ -112,5 +127,6 @@ int main(int argc, char** argv)
free_matrix(&F, nx);
free_matrix(&G, nx);
free_matrix(&RHS, nx);
+ free_matrix(&RES, nx);
return 0;
}
(DIR) diff --git a/src/solution.c b/src/solution.c
t@@ -1,3 +1,5 @@
+#include <stdio.h>
+#include <stdlib.h>
#include <math.h>
/* Second order central difference (d^2 u)/(dx^2) */
t@@ -71,19 +73,142 @@ void compute_f_g(double **F, double **G,
double gx, double gy, double dt, double gamma)
{
int i, j;
- for (i=1; i<=nx; i++) {
- for (j=1; j<=ny; j++) {
+ for (i=1; i<nx; i++)
+ for (j=1; j<=ny; j++)
F[i][j] = U[i][j]
+ dt*(1.0/re*(ddu_dxx(U, i, j, dx) + ddu_dyy(U, i, j, dy))
- duu_dx(U, i, j, dx, gamma)
- duv_dy(U, V, i, j, dy, gamma)
+ gx);
+ for (i=1; i<=nx; i++)
+ for (j=1; j<ny; j++)
G[i][j] = V[i][j]
+ dt*(1.0/re*(ddv_dxx(V, i, j, dx) + ddv_dyy(V, i, j, dy))
- duv_dx(U, V, i, j, dx, gamma)
- dvv_dy(V, i, j, dy, gamma)
+ gy);
+}
+
+/* Find the right hand side value of the Poisson equation */
+void compute_rhs(double **RHS, double **F, double **G,
+ double dt, double dx, double dy, int nx, int ny)
+{
+ int i, j;
+ for (i=1; i<=nx; i++) {
+ for (j=1; j<=ny; j++) {
+ RHS[i][j] = 1.0/dt*(
+ (F[i][j] - F[i-1][j])/dx + (G[i][j] - G[i][j-1])/dy);
}
}
}
+
+void e_parameters(double *e_left, double *e_right,
+ double *e_bottom, double *e_top,
+ int i, int j, int nx, int ny)
+{
+ if (i == 1) {
+ *e_left = 0.0;
+ } else if (i > 1) {
+ *e_left = 1.0;
+ } else {
+ fprintf(stderr, "e_parameters: i value %d not understood\n", i);
+ exit(1);
+ }
+
+ if (i < nx) {
+ *e_right = 1.0;
+ } else if (i == nx) {
+ *e_right = 0.0;
+ } else {
+ fprintf(stderr, "e_parameters: i value %d not understood\n", i);
+ exit(1);
+ }
+
+ if (j == 1) {
+ *e_bottom = 0.0;
+ } else if (j > 1) {
+ *e_bottom = 1.0;
+ } else {
+ fprintf(stderr, "e_parameters: j value %d not understood\n", j);
+ exit(1);
+ }
+
+ if (j < ny) {
+ *e_top = 1.0;
+ } else if (j == ny) {
+ *e_top = 0.0;
+ } else {
+ fprintf(stderr, "e_parameters: j value %d not understood\n", j);
+ exit(1);
+ }
+}
+
+/* Compute a successive overrelaxation (SOR) cycle */
+void sor_cycle(double **P, double **RHS, double omega,
+ int nx, int ny, double dx, double dy)
+{
+ int i, j;
+ double e_left, e_right, e_top, e_bottom;
+
+ for (i=1; i<=nx; i++) {
+ for (j=1; j<=ny; j++) {
+ e_parameters(&e_left, &e_right, &e_bottom, &e_top, i, j, nx, ny);
+ P[i][j] = (1.0 - omega)*P[i][j] +
+ omega/((e_right + e_left)/(dx*dx) + (e_top + e_bottom)/(dy*dy))
+ *( (e_right*P[i+1][j] + e_left*P[i-1][j])/(dx*dx) +
+ (e_top*P[i][j+1] + e_bottom*P[i][j-1])/(dy*dy) -
+ RHS[i][j]);
+ }
+ }
+}
+
+void calculate_residuals(double **RES, double **P, double **RHS,
+ int nx, int ny, double dx, double dy)
+{
+ int i, j;
+ double e_left, e_right, e_top, e_bottom;
+
+ for (i=1; i<=nx; i++) {
+ for (j=1; j<=ny; j++) {
+ e_parameters(&e_left, &e_right, &e_bottom, &e_top, i, j, nx, ny);
+ RES[i][j] =
+ (e_right*(P[i+1][j] - P[i][j]) - e_left*(P[i][j] - P[i-1][j]))
+ /(dx*dx) +
+ (e_top*(P[i][j+1] - P[i][j]) - e_bottom*(P[i][j] - P[i][j-1]))
+ /(dy*dy) - RHS[i][j];
+ }
+ }
+}
+
+double max_residual_norm(double **RES, double **P, double **RHS,
+ int nx, int ny, double dx, double dy)
+{
+ int i, j;
+ double val;
+ double max = -1.0e9;
+
+ calculate_residuals(RES, P, RHS, nx, ny, dx, dy);
+
+ for (i=1; i<=nx; i++) {
+ for (j=1; j<=ny; j++) {
+ val = fabs(RES[i][j]);
+ if (val > max)
+ max = val;
+ }
+ }
+ return max;
+}
+
+void correct_velocities(double **U, double **V, double **F, double **G,
+ double **P, int nx, int ny, double dt, double dx, double dy)
+{
+ int i, j;
+ for (i=1; i<nx; i++)
+ for (j=1; j<=ny; j++)
+ U[i][j] = F[i][j] - dt/dx*(P[i+1][j] - P[i][j]);
+
+ for (i=1; i<=nx; i++)
+ for (j=1; j<ny; j++)
+ V[i][j] = G[i][j] - dt/dy*(P[i][j+1] - P[i][j]);
+}
(DIR) diff --git a/src/solution.h b/src/solution.h
t@@ -6,4 +6,15 @@ void compute_f_g(double **F, double **G,
int nx, int ny, double dx, double dy,
double gx, double gy, double dt, double gamma);
+void compute_rhs(double **RHS, double **F, double **G,
+ double dt, double dx, double dy, int nx, int ny);
+
+void sor_cycle(double **P, double **RHS, double omega,
+ int nx, int ny, double dx, double dy);
+
+double max_residual_norm(double **RES, double **P, double **RHS,
+ int nx, int ny, double dx, double dy);
+
+void correct_velocities(double **U, double **V, double **F, double **G,
+ double **P, int nx, int ny, double dt, double dx, double dy);
#endif