#include
#include
#include
#include
#include
#include “alloc.h”
#include “boundary.h”
#include “datadef.h”
#include “init.h”
#include “simulation.h”
void write_bin(float **u, float **v, float **p, char **flag,
int imax, int jmax, float xlength, float ylength, char *file);
int read_bin(float **u, float **v, float **p, char **flag,
int imax, int jmax, float xlength, float ylength, char *file);
static void print_usage(void);
static void print_version(void);
static void print_help(void);
static char *progname;
int proc = 0; /* Rank of the current process */
int nprocs = 0; /* Number of processes in communicator */
int *ileft, *iright; /* Array bounds for each processor */
#define PACKAGE “karman”
#define VERSION “1.0”
/* Command line options */
static struct option long_opts[] = {
{ “del-t”, 1, NULL, ‘d’ },
{ “help”, 0, NULL, ‘h’ },
{ “imax”, 1, NULL, ‘x’ },
{ “infile”, 1, NULL, ‘i’ },
{ “jmax”, 1, NULL, ‘y’ },
{ “outfile”, 1, NULL, ‘o’ },
{ “t-end”, 1, NULL, ‘t’ },
{ “verbose”, 1, NULL, ‘v’ },
{ “version”, 1, NULL, ‘V’ },
{ 0, 0, 0, 0 }
};
#define GETOPTS “d:hi:o:t:v:Vx:y:”
int main(int argc, char *argv[])
{
int verbose = 1; /* Verbosity level */
float xlength = 22.0; /* Width of simulated domain */
float ylength = 4.1; /* Height of simulated domain */
int imax = 660; /* Number of cells horizontally */
int jmax = 120; /* Number of cells vertically */
char *infile; /* Input raw initial conditions */
char *outfile; /* Output raw simulation results */
float t_end = 2.1; /* Simulation runtime */
float del_t = 0.003; /* Duration of each timestep */
float tau = 0.5; /* Safety factor for timestep control */
int itermax = 100; /* Maximum number of iterations in SOR */
float eps = 0.001; /* Stopping error threshold for SOR */
float omega = 1.7; /* Relaxation parameter for SOR */
float gamma = 0.9; /* Upwind differencing factor in PDE
discretisation */
float Re = 150.0; /* Reynolds number */
float ui = 1.0; /* Initial X velocity */
float vi = 0.0; /* Initial Y velocity */
float t, delx, dely;
int i, j, itersor = 0, ifluid = 0, ibound = 0;
float res;
float **u, **v, **p, **rhs, **f, **g;
char **flag;
int init_case, iters = 0;
int show_help = 0, show_usage = 0, show_version = 0;
progname = argv[0];
infile = strdup(“karman.bin”);
outfile = strdup(“karman.bin”);
int optc;
while ((optc = getopt_long(argc, argv, GETOPTS, long_opts, NULL)) != -1) {
switch (optc) {
case ‘h’:
show_help = 1;
break;
case ‘V’:
show_version = 1;
break;
case ‘v’:
verbose = atoi(optarg);
break;
case ‘x’:
imax = atoi(optarg);
break;
case ‘y’:
jmax = atoi(optarg);
break;
case ‘i’:
free(infile);
infile = strdup(optarg);
break;
case ‘o’:
free(outfile);
outfile = strdup(optarg);
break;
case ‘d’:
del_t = atof(optarg);
break;
case ‘t’:
t_end = atof(optarg);
break;
default:
show_usage = 1;
}
}
if (show_usage || optind < argc) {
print_usage();
return 1;
}
if (show_version) {
print_version();
if (!show_help) {
return 0;
}
}
if (show_help) {
print_help();
return 0;
}
delx = xlength/imax;
dely = ylength/jmax;
/* Allocate arrays */
u = alloc_floatmatrix(imax+2, jmax+2);
v = alloc_floatmatrix(imax+2, jmax+2);
f = alloc_floatmatrix(imax+2, jmax+2);
g = alloc_floatmatrix(imax+2, jmax+2);
p = alloc_floatmatrix(imax+2, jmax+2);
rhs = alloc_floatmatrix(imax+2, jmax+2);
flag = alloc_charmatrix(imax+2, jmax+2);
if (!u || !v || !f || !g || !p || !rhs || !flag) {
fprintf(stderr, "Couldn't allocate memory for matrices.\n");
return 1;
}
/* Read in initial values from a file if it exists */
init_case = read_bin(u, v, p, flag, imax, jmax, xlength, ylength, infile);
if (init_case > 0) {
/* Error while reading file */
return 1;
}
if (init_case < 0) {
/* Set initial values if file doesn't exist */
for (i=0;i<=imax+1;i++) {
for (j=0;j<=jmax+1;j++) {
u[i][j] = ui;
v[i][j] = vi;
p[i][j] = 0.0;
}
}
init_flag(flag, imax, jmax, delx, dely, &ibound);
apply_boundary_conditions(u, v, flag, imax, jmax, ui, vi);
}
/* Main loop */
for (t = 0.0; t < t_end; t += del_t, iters++) {
set_timestep_interval(&del_t, imax, jmax, delx, dely, u, v, Re, tau);
ifluid = (imax * jmax) - ibound;
compute_tentative_velocity(u, v, f, g, flag, imax, jmax,
del_t, delx, dely, gamma, Re);
compute_rhs(f, g, rhs, flag, imax, jmax, del_t, delx, dely);
if (ifluid > 0) {
itersor = poisson(p, rhs, flag, imax, jmax, delx, dely,
eps, itermax, omega, &res, ifluid);
} else {
itersor = 0;
}
if (proc == 0 && verbose > 1) {
printf(“%d t:%g, del_t:%g, SOR iters:%3d, res:%e, bcells:%d\n”,
iters, t+del_t, del_t, itersor, res, ibound);
}
update_velocity(u, v, f, g, p, flag, imax, jmax, del_t, delx, dely);
apply_boundary_conditions(u, v, flag, imax, jmax, ui, vi);
} /* End of main loop */
if (outfile != NULL && strcmp(outfile, “”) != 0 && proc == 0) {
write_bin(u, v, p, flag, imax, jmax, xlength, ylength, outfile);
}
free_matrix(u);
free_matrix(v);
free_matrix(f);
free_matrix(g);
free_matrix(p);
free_matrix(rhs);
free_matrix(flag);
return 0;
}
/* Save the simulation state to a file */
void write_bin(float **u, float **v, float **p, char **flag,
int imax, int jmax, float xlength, float ylength, char* file)
{
int i;
FILE *fp;
fp = fopen(file, “wb”);
if (fp == NULL) {
fprintf(stderr, “Could not open file ‘%s’: %s\n”, file,
strerror(errno));
return;
}
fwrite(&imax, sizeof(int), 1, fp);
fwrite(&jmax, sizeof(int), 1, fp);
fwrite(&xlength, sizeof(float), 1, fp);
fwrite(&ylength, sizeof(float), 1, fp);
for (i=0;i