tmain.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
---
tmain.c (4384B)
---
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <ctype.h>
4 #include <unistd.h>
5 #include <string.h>
6 #include "file_io.h"
7 #include "utility.h"
8 #include "boundary.h"
9 #include "solution.h"
10
11 #define VERSION 0.1
12
13 /*#define DEBUG*/
14
15 int main(int argc, char** argv)
16 {
17 double t, t_file, t_end, tau;
18 int itermax;
19 double epsilon, omega, gamma;
20 double gx, gy, re;
21 int w_left, w_right, w_top, w_bottom;
22 double dx, dy;
23 int nx, ny;
24 double **P, **U, **V;
25 double **F, **G, **RHS, **RES;
26
27 double dt = 0.0;
28 double res;
29 long n = 0;
30 int it = 0;
31 int nfile = 0;
32 double t_file_elapsed = 0.0;
33 char filename[50];
34 char *simulation_id;
35 char *dot;
36
37 int c;
38 while ((c = getopt(argc, argv, "hv")) != -1)
39 switch (c)
40 {
41 case 'h':
42 printf("usage: %s [OPTIONS] FILENAME\n"
43 "where FILENAME is a valid input file.\n"
44 "-h\tDisplay help\n"
45 "-v\tDisplay version information\n", argv[0]);
46 return 0;
47 break;
48 case 'v':
49 printf("%s: A Navier Stokes solver in 2D using finite differences, "
50 "version %.1f\n"
51 "Written by Anders Damsgaard, "
52 "https://github.com/anders-dc/ns2dfd\n", argv[0], VERSION);
53 return 0;
54 break;
55 case '?':
56 if (isprint(optopt))
57 fprintf(stderr, "Unknown option `-%c`.\n", optopt);
58 else
59 fprintf(stderr, "Unknown option character `\\x%x`.\n", optopt);
60 return 1;
61
62 default:
63 abort();
64 }
65
66 if (optind == argc - 1) {
67 read_file(argv[optind], &t, &t_end, &t_file, &tau, &itermax,
68 &epsilon, &omega, &gamma,
69 &gx, &gy, &re, &w_left, &w_right, &w_top, &w_bottom,
70 &dx, &dy, &nx, &ny, &P, &U, &V);
71 } else {
72 fprintf(stderr, "%s: no input file specified.\n", argv[0]);
73 return 1;
74 }
75
76 simulation_id = argv[optind];
77 dot = strchr(simulation_id, '.');
78 dot[0] = '\0';
79 printf("%s: simulation id: `%s`\n", argv[0], simulation_id);
80
81 allocate_matrix(&F, nx+2, ny+2);
82 allocate_matrix(&G, nx+2, ny+2);
83 allocate_matrix(&RHS, nx+2, ny+2);
84 allocate_matrix(&RES, nx+2, ny+2);
85
86 #ifdef DEBUG
87 print_matrix("V", P, nx+2, ny+2);
88 print_matrix("U", U, nx+2, ny+2);
89 print_matrix("V", V, nx+2, ny+2);
90 #endif
91
92 while (t <= t_end+dt) {
93
94 dt = select_time_step(tau, re, dx, dy, nx, ny, U, V);
95
96 if (t_file_elapsed >= t_file || n == 0) {
97 sprintf(filename, "%s%05d.dat", simulation_id, nfile);
98 write_file(filename, &t, &t_end, &t_file, &tau, &itermax,
99 &epsilon, &omega, &gamma,
100 &gx, &gy, &re, &w_left, &w_right, &w_top, &w_bottom,
101 &dx, &dy, &nx, &ny, &P, &U, &V);
102 t_file_elapsed = 0.0;
103 nfile++;
104 }
105
106 /*printf("\rt = %f/%.2f s, dt = %f s, it = %d, last output: %d ",*/
107 printf("t = %f/%.2f s, dt = %f s, it = %d, last output: %d\n",
108 t, t_end, dt, it, nfile-1);
109
110 set_boundary_conditions(w_left, w_right, w_top, w_bottom, P, U, V,
111 nx, ny);
112
113 compute_f_g(F, G, U, V, P, re, nx, ny, dx, dy, gx, gy, dt, gamma);
114
115 #ifdef DEBUG
116 print_matrix("F", F, nx+2, ny+2);
117 print_matrix("G", G, nx+2, ny+2);
118 #endif
119
120 compute_rhs(RHS, F, G, dt, dx, dy, nx, ny);
121
122 #ifdef DEBUG
123 print_matrix("RHS", RHS, nx+2, ny+2);
124 #endif
125
126 it = 0;
127 res = 1.0e9;
128 while ((it <= itermax) && (res >= epsilon)) {
129
130 sor_cycle(P, RHS, omega, nx, ny, dx, dy);
131
132 res = max_residual_norm(RES, P, RHS, nx, ny, dx, dy);
133
134 #ifdef DEBUG
135 printf("\nit = %d\n", it);
136 print_matrix("P", P, nx+2, ny+2);
137 print_matrix("RES", RES, nx+2, ny+2);
138 #endif
139
140 it++;
141 }
142
143 correct_velocities(U, V, F, G, P, nx, ny, dt, dx, dy);
144
145 #ifdef DEBUG
146 print_matrix("U", U, nx+2, ny+2);
147 print_matrix("V", V, nx+2, ny+2);
148 #endif
149
150 t += dt;
151 n++;
152 t_file_elapsed += dt;
153 }
154 puts("\n");
155
156 free_memory(P, U, V, nx);
157 free_matrix(&F, nx+2);
158 free_matrix(&G, nx+2);
159 free_matrix(&RHS, nx+2);
160 free_matrix(&RES, nx+2);
161 return 0;
162 }