tImplemented adaptive time step - 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 08a48407f92cd5c0c603c8e5207ab92b8a5df5a5
 (DIR) parent 712b7d25f92ae5ecd752f39c2d1578a39d973457
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Sun,  2 Mar 2014 20:51:49 +0100
       
       Implemented adaptive time step
       
       Diffstat:
         M src/Makefile                        |       3 ++-
         M src/main.c                          |      23 ++++++++++++++++-------
         M src/utility.c                       |      30 ++++++++++++++++++++++++++++++
         M src/utility.h                       |       3 +++
       
       4 files changed, 51 insertions(+), 8 deletions(-)
       ---
 (DIR) diff --git a/src/Makefile b/src/Makefile
       t@@ -1,9 +1,10 @@
        CFLAGS=-Wall -pedantic -g
       +LDLIBS=-lm
        
        BIN=../ns2dfd
        
        $(BIN): main.o file_io.o utility.o
       -        $(CC) $(CFLAGS) $^ -o $@
       +        $(CC) $(CFLAGS) $(LDLIBS) $^ -o $@
        
        main.o: file_io.h utility.h
        file_io.o: utility.h
 (DIR) diff --git a/src/main.c b/src/main.c
       t@@ -18,6 +18,9 @@ int main(int argc, char** argv)
            int nx, ny;
            double **P, **U, **V;
        
       +    double dt;
       +    long n = 0;
       +
            int c;
            while ((c = getopt(argc, argv, "hv")) != -1)
                switch (c)
       t@@ -30,7 +33,7 @@ int main(int argc, char** argv)
                    return 0;
                    break;
                case 'v':
       -            printf("%s: Navier Stokes solver in 2D using finite differences, "
       +            printf("%s: A Navier Stokes solver in 2D using finite differences, "
                            "version %.1f\n"
                            "Written by Anders Damsgaard, "
                            "https://github.com/anders-dc/ns2dfd\n", argv[0], VERSION);
       t@@ -57,14 +60,20 @@ int main(int argc, char** argv)
                return 1;
            }
        
       -    printf("omega = %f\n", omega);
       -    printf("P[0][0] = %f\n", P[0][0]);
       -    printf("V[%d][%d] = %f\n", nx, ny, V[nx][ny]);
       -
       -    write_file("unnamed2.dat", &t, &t_end, &tau, &itermax,
       +    /*write_file("unnamed2.dat", &t, &t_end, &tau, &itermax,
                    &epsilon, &omega, &gamma, 
                    &gx, &gy, &re, &w_left, &w_right, &w_top, &w_bottom,
       -            &dx, &dy, &nx, &ny, &P, &U, &V);
       +            &dx, &dy, &nx, &ny, &P, &U, &V);*/
       +
       +    /*while (t < t_end) {*/
       +
       +        dt = select_time_step(tau, re, dx, dy, nx, ny, U, V);
       +        printf("dt = %f\n", dt);
       +
       +
       +        /*t += dt;
       +        n++;
       +    }*/
        
            free_memory(P, U, V, nx);
            return 0;
 (DIR) diff --git a/src/utility.c b/src/utility.c
       t@@ -1,5 +1,6 @@
        #include <stdio.h>
        #include <stdlib.h>
       +#include <math.h>
        
        int allocate_matrix(double ***M, int nx, int ny)
        {
       t@@ -47,3 +48,32 @@ void free_memory(double **P, double **U, double **V, int nx)
            free_matrix(&U, nx);
            free_matrix(&V, nx);
        }
       +
       +double max_abs_value(double **M, int nx, int ny)
       +{
       +    int i, j;
       +    double tmp;
       +    double max = -1.0e19;
       +    for (j=1; j<=ny; j++) {
       +        for (i=1; i<=nx; i++) {
       +            tmp = fabs(M[i][j]);
       +            if (tmp > max) max = tmp;
       +        }
       +    }
       +    return max;
       +}
       +
       +double select_time_step(double tau, double re, double dx, double dy,
       +        int nx, int ny, double **U, double **V)
       +{
       +    double t_diff, t_adv_u, t_adv_v;
       +    double u_max, v_max;
       +    u_max = max_abs_value(U, nx, ny);
       +    v_max = max_abs_value(V, nx, ny);
       +
       +    t_diff = re/2.0*(1.0/(1.0/(dx*dx) + (1.0/(dy*dy))));
       +    t_adv_u = dx/(u_max + 1.0e-12);
       +    t_adv_v = dx/(v_max + 1.0e-12);
       +
       +    return fmin(t_diff, fmin(t_adv_u, t_adv_v));
       +}
 (DIR) diff --git a/src/utility.h b/src/utility.h
       t@@ -4,4 +4,7 @@
        int allocate_memory(double ***P, double ***U, double ***V, int nx, int ny);
        void free_memory(double **P, double **U, double **V, int nx);
        
       +double select_time_step(double tau, double re, double dx, double dy,
       +        int nx, int ny, double **U, double **V);
       +
        #endif