tutility.c - ns2dfd - 2D finite difference Navier Stokes solver for fluid dynamics
  HTML git clone git://src.adamsgaard.dk/ns2dfd
   DIR Log
   DIR Files
   DIR Refs
   DIR LICENSE
       ---
       tutility.c (2516B)
       ---
            1 #include <stdio.h>
            2 #include <stdlib.h>
            3 #include <math.h>
            4 
            5 int allocate_matrix(double ***M, int nx, int ny)
            6 {
            7     int i;
            8     (*M) = (double**) malloc(sizeof(double*)*nx);
            9     if (*M == NULL)
           10         return 1;
           11 
           12     for (i=0; i<nx; i++) {
           13         (*M)[i] = (double*) malloc(sizeof(double)*ny);
           14         if ((*M)[i] == NULL)
           15             return 1;
           16     }
           17     return 0;
           18 }
           19 
           20 void free_matrix(double ***M, int nx)
           21 {
           22     int i;
           23     for (i=0; i<nx; i++)
           24         free((*M)[i]);
           25     free(*M);
           26 }
           27 
           28 int allocate_memory(double ***P, double ***U, double ***V, int nx, int ny)
           29 {
           30     if (allocate_matrix(P, nx+2, ny+2) != 0) {
           31         fprintf(stderr, "allocate_memory: Could not allocate memory for P\n");
           32         return 1;
           33     }
           34     if (allocate_matrix(U, nx+2, ny+2) != 0) {
           35         fprintf(stderr, "allocate_memory: Could not allocate memory for U\n");
           36         return 1;
           37     }
           38     if (allocate_matrix(V, nx+2, ny+2) != 0) {
           39         fprintf(stderr, "allocate_memory: Could not allocate memory for V\n");
           40         return 1;
           41     }
           42     return 0;
           43 }
           44 
           45 void free_memory(double **P, double **U, double **V, int nx)
           46 {
           47     free_matrix(&P, nx+2);
           48     free_matrix(&U, nx+2);
           49     free_matrix(&V, nx+2);
           50 }
           51 
           52 double max_abs_value(double **M, int nx, int ny)
           53 {
           54     int i, j;
           55     double tmp;
           56     double max = -1.0e19;
           57     for (j=1; j<=ny; j++) {
           58         for (i=1; i<=nx; i++) {
           59             tmp = fabs(M[i][j]);
           60             if (tmp > max) max = tmp;
           61         }
           62     }
           63     return max;
           64 }
           65 
           66 double select_time_step(double tau, double re, double dx, double dy,
           67         int nx, int ny, double **U, double **V)
           68 {
           69     double t_diff, t_adv_u, t_adv_v;
           70     double u_max, v_max;
           71     double dt;
           72     u_max = max_abs_value(U, nx, ny);
           73     v_max = max_abs_value(V, nx, ny);
           74 
           75     t_diff = re/2.0*(1.0/(1.0/(dx*dx) + (1.0/(dy*dy))));
           76     t_adv_u = dx/(u_max + 1.0e-12);
           77     t_adv_v = dx/(v_max + 1.0e-12);
           78 
           79     dt = fmin(t_diff, fmin(t_adv_u, t_adv_v));
           80 
           81     if (dt < 1.0e-12) {
           82         fprintf(stderr, "select_time_step: Error, the time step is too small: "
           83                 "%f s.\n"
           84                 "tau = %f, Re = %f, u_max = %f, v_max = %f\n"
           85                 "t_diff = %f s, t_adv_u = %f s, t_adv_v = %f s\n",
           86                 dt, tau, re, u_max, v_max, t_diff, t_adv_u, t_adv_v);
           87         exit(1);
           88     }
           89     return dt;
           90 }
           91 
           92 void print_matrix(char* description, double **M, int nx, int ny)
           93 {
           94     int i, j;
           95     printf("%s:\n", description);
           96     for (j=0; j<ny; j++) {
           97         for (i=0; i<nx; i++) {
           98             printf("%.2f\t", M[i][j]);
           99         }
          100         puts("\n");
          101     }
          102 }