tadded NaN check in residual, added verbose array debugging - 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 e3c43ff247f84337b66a8782d923b7e293ae03d5
 (DIR) parent 8dde47c85ef90b41df4d605fb96f55bce13628be
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Tue,  4 Mar 2014 10:55:06 +0100
       
       added NaN check in residual, added verbose array debugging
       
       Diffstat:
         M src/main.c                          |      13 ++++++++++++-
         M src/solution.c                      |       6 ++++++
         M src/utility.c                       |      12 ++++++++++++
         M src/utility.h                       |       2 ++
       
       4 files changed, 32 insertions(+), 1 deletion(-)
       ---
 (DIR) diff --git a/src/main.c b/src/main.c
       t@@ -81,6 +81,7 @@ int main(int argc, char** argv)
            allocate_matrix(&RHS, nx+2, ny+2);
            allocate_matrix(&RES, nx+2, ny+2);
        
       +    print_matrix("P", P, nx+2, ny+2);
            while (t <= t_end+dt) {
        
                dt = select_time_step(tau, re, dx, dy, nx, ny, U, V);
       t@@ -95,15 +96,19 @@ int main(int argc, char** argv)
                    nfile++;
                }
        
       -        printf("\rt = %f/%.2f s, dt = %f s, it = %d, last output: %d   ",
       +        /*printf("\rt = %f/%.2f s, dt = %f s, it = %d, last output: %d   ",*/
       +        printf("t = %f/%.2f s, dt = %f s, it = %d, last output: %d\n",
                        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);
       +        print_matrix("F", F, nx+2, ny+2);
       +        print_matrix("G", G, nx+2, ny+2);
        
                compute_rhs(RHS, F, G, dt, dx, dy, nx, ny);
       +        print_matrix("RHS", RHS, nx+2, ny+2);
        
                it = 0;
                res = 1.0e9;
       t@@ -113,10 +118,16 @@ int main(int argc, char** argv)
        
                    res = max_residual_norm(RES, P, RHS, nx, ny, dx, dy);
        
       +            printf("\nit = %d\n", it);
       +            print_matrix("P", P, nx+2, ny+2);
       +            print_matrix("RES", RES, nx+2, ny+2);
       +
                    it++;
                }
        
                correct_velocities(U, V, F, G, P, nx, ny, dt, dx, dy);
       +        print_matrix("U", U, nx+2, ny+2);
       +        print_matrix("V", V, nx+2, ny+2);
        
                t += dt;
                n++;
 (DIR) diff --git a/src/solution.c b/src/solution.c
       t@@ -177,6 +177,12 @@ void calculate_residuals(double **RES, double **P, double **RHS,
                        /(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];
       +            if (RES[i][j] != RES[i][j]) {
       +                fprintf(stderr, "calculate_residuals: "
       +                        "Error, residual is %f in cell i=%d, j=%d\n",
       +                        RES[i][j], i, j);
       +                exit(1);
       +            }
                }
            }
        }
 (DIR) diff --git a/src/utility.c b/src/utility.c
       t@@ -88,3 +88,15 @@ double select_time_step(double tau, double re, double dx, double dy,
            }
            return dt;
        }
       +
       +void print_matrix(char* description, double **M, int nx, int ny)
       +{
       +    int i, j;
       +    printf("%s:\n", description);
       +    for (j=0; j<ny; j++) {
       +        for (i=0; i<nx; i++) {
       +            printf("%.2f\t", M[i][j]);
       +        }
       +        puts("\n");
       +    }
       +}
 (DIR) diff --git a/src/utility.h b/src/utility.h
       t@@ -10,4 +10,6 @@ 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);
        
       +void print_matrix(char* description, double **M, int nx, int ny);
       +
        #endif