/*
This file is part of EZRide.
Copyright 2007 Dwight Barkley,
Copyright 2010 Vadim Biktashev & Andrew Foulkes.
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published
by the Free Software Foundation, either version 3 of the License,
or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see .
*/
#include
#include
#include
#include
#include
#include
#include
#include
#include "ezride.h"
/*
* Global variables used throughout the EZride Code (see ezride.h)
* -------------------------------------------------------------------- */
Real *fields, *u, *v, *u1, *v1, *sigma_u, *sigma_v, *tip_x, *tip_y;
Real bet, gam, a_param, b_param, eps, one_o_eps, delta, grid_h, dt, Utip, Vtip,
dt_o_eps, dt_o_2eps, one_o_a, b_o_a, dt_o_wh2, dtDv_o_wh2, dteps, dt_o_2h,
LX2, LY2, X0, Y0, TH0, TH1, stab_cx, stab_cy, stab_om, u_root, v_root, one_o_h,
X1, Y1, Au, Av, a1, ResU, ResV, Om, b_o_a2, a2, one_o_a2, dt_h, x_tip, y_tip,
cx, cy, omega;
Real ulo, uhi, umin, umax, vmin, vmax; /* visualization benchmarks, model-dependent */
int nx, ny, field_size, nsteps, istep, write_tip, write_quot, ride_on,
graphics, show_tip, autostart, verbose, ntips, NX2, NY2, XX, YY,
plot_step, image_step, x_inc, y_inc, NBC, perturb_on, cxdiff, cydiff, omdiff;
/* Global variables for this file only
* ----------------------------------- */
static Real length_x, length_y, Dv, ts;
static int hist_step, hist_x, hist_y, binary_write, ini_x, ini_y;
static FILE *history_file, *tip_file, *quot_file, *recon_file;
/* Private functions
* ----------------- */
static void Initialize (void);
static void Allocate_memory (void);
static void Finish_up (void);
static void Recon_tip (void);
static void Write_history (int wrt_step);
static void Write_fc (void);
static void Write_fields (Real *ufiled, Real *vfield, char *filename);
static void Read_ic (void);
static void Read_fields (Real *ufield, Real *vfield, FILE *fp);
static void Generate_ic (void);
#if FHN
static Real Cube (Real g);
static void Root (Real x, Real y);
#endif
#define NEXT_LINE(fp) while(getc(fp)!='\n'); /* macro used to skip to end
of input line */
/* ========================================================================= */
int main(void)
{
Initialize();
time_t start,end;
double time_diff,secs;
int minute;
time(&start);
for(istep=0; istepNX ||
hist_y<1 || hist_y>NY ) ) {
fprintf(stderr,"***** history point out of range ******\n");
exit(1);
}
if(ts > 1.0 ) {
fprintf(stderr,"***** ts > 1 (the diffusion stability limit) ******\n");
exit(1);
}
if (plot_step == 0) {
fprintf(stderr,"***** plot_step cannot be 0 ******\n");
plot_step = nsteps+1;
}
/* ------------
* Final things
* ------------ */
Allocate_memory();
Read_ic();
Step_ini();
if(hist_step) history_file = fopen("history.dat", "w");
if(write_tip) tip_file = fopen("tip.dat", "w");
if(write_quot) quot_file = fopen("quot.dat", "w");
recon_file = fopen("recon.dat","w");
fprintf(recon_file, "# %3s %15s %15s %15s %15s %15s %15s %15s %15s %15s %15s\n",
"t", "X0", "Y0", "TH0", "X1", "Y1", "TH1", "X2", "Y2", "TH2", "ride_on");
#if GRAPHICS
if (graphics) Draw_ini(initial_field);
#endif
Tip_ini();
stab_cx = 0.5*grid_h*grid_h/(dt);
stab_cy = 0.5*grid_h*grid_h/(dt);
stab_om = 1/(NX*dt);
#if FHN
Root(bet,gam);
ulo = -1.0;
uhi = 1.0;
umin = -sqrt(3.0);
vmin = -2.0/3.0;
umax = sqrt(3.0);
vmax = 2.0/3.0;
#else
u_root = 0;
v_root = 0;
ulo = 0.1;
uhi = 0.9;
umin = 0.0;
vmin = 0.0;
umax = 1.0;
vmax = 0.8;
#endif
printf("stab_cx = %g, stab_cy = %g, stab_om = %g,\n", stab_cx, stab_cy, stab_om);
if (x_inc==0 && y_inc==0){
printf("\n\nWARNING: Second pinning point has been chosen to be at (0,0)\n");
printf("which will result in division by zero if advection is switched on!\n");
}
}
/* ========================================================================= */
void Allocate_memory (void)
{
/* -----------------------------------------------------------------------
* There are NX and NY actual grid points in each direction.
* field_size = (NX+2)*(NY+2) because of the way boundaries are
* treated. See ezride.h for definition of NX, NY.
*
* I allocate one big block of memory for everything (u, v, sigma_u,
* and if needed sigma_v) and then set appropriate pointers. The
* handling of memory can be changed with appropriate changes to the
* macros INDEX etc in ezride.h.
* ----------------------------------------------------------------------- */
int num_fields, next_field;
field_size = FIELD_SIZE;
num_fields=6;
#if V_DIFF_ON
num_fields+=2;
#endif
fields = (Real *) calloc(num_fields*field_size, sizeof(Real));
if (fields != NULL) {
/* ---------------------------------------------------
* We have one big block of memory. Set some pointers
* --------------------------------------------------- */
next_field=0;
u = fields + field_size*(next_field++);
u1 = fields + field_size*(next_field++);
v = fields + field_size*(next_field++);
v1 = fields + field_size*(next_field++);
sigma_u = fields + field_size*(next_field++); next_field++;
#if V_DIFF_ON
sigma_v = fields + field_size*(next_field++); next_field++;
#endif
if (next_field != num_fields) {
fprintf(stderr,"fields miscounted: %d vs %d\n",next_field,num_fields);
exit(1);
}
if(verbose>2) printf("Memory allocated in one block.\n");
}
else {
/* -------------------
* Did not get memory.
* ------------------- */
fprintf(stderr, "\n ****Error: Not enough memory available\n");
exit(-1);
}
}
/* ========================================================================= */
void Finish_up (void)
{
#if GRAPHICS
if (graphics) QuitX();
#endif
Write_fc();
if(hist_step) fclose(history_file);
if(write_tip) fclose(tip_file);
if(write_quot) fclose(quot_file);
}
/* ========================================================================= */
void Write_tip_data (float x, float y, float th)
{
/* This routine for outputting the tip data is called from Find_tips() in
* eztip.c where the tip is computed. */
fprintf(tip_file, "%.5f %.5f %.5f %.5f\n", (float)(dt*istep), x, y, th);
}
/* ========================================================================= */
void Write_quot (float cx, float cy, float omega)
{
if (verbose)
printf("t=%.5f cx=%.10f cy=%.10f om=%.10f\n", (float)(dt*istep), (float)cx, (float)cy, (float)omega);
if (quot_file)
fprintf(quot_file, "%.5f %.10f %.10f %.10f\n", (float)(dt*istep), (float)cx, (float)cy, (float)omega);
fflush(quot_file);
}
/* ========================================================================= */
void Recon_tip (void) {
Real cosTH1=cos(TH1);
Real sinTH1=sin(TH1);
Real X2 = X1 + cosTH1*X0 - sinTH1*Y0;
Real Y2 = Y1 + sinTH1*X0 + cosTH1*Y0;
Real TH2 = TH1 + TH0;
fprintf(recon_file, "%.5f %15.10f %15.10f %15.10f %15.10f %15.10f %15.10f %15.10f %15.10f %15.10f %1u\n",
(float)(dt*istep), X0, Y0, TH0, X1, Y1, TH1, X2, Y2, TH2, ride_on);
}
/* ========================================================================= */
#if FHN
Real Cube (Real g)
{
Real g1,g2,g3,g4,cube;
if (g==0) return 0.0;
g1 = log(abs(g));
g2 = 1./3.;
g3 = g1*g2;
g4 = abs(g)/g;
cube = exp(g3)*g4;
return(cube);
}
/* ========================================================================= */
void Root (Real x, Real y)
{
Real a,b,c,d;
if (y==0){ u_root = -1;
v_root = -2./3.;}
else{
a = (3*(1-y)/y);
b = (3*x/y);
c = (-(108*b)+(12*sqrt((12*a*a*a)+(81*b*b))));
if (c==0){ u_root = -1;
v_root = -2./3.;}
else{d = Cube(c);
if (d==0){ u_root = -1;
v_root = -2./3.;}
else { u_root = ((d/6)-(2*a/d));
v_root = (u_root-((u_root*u_root*u_root)/3));}
}
}
}
#endif
/* ========================================================================= */
void Write_history (int wrt_step)
{
fprintf(history_file, "%.5f %.5f %.5f\n", dt*(Real)wrt_step,
U(hist_x,hist_y), V(hist_x,hist_y) );
if(verbose>2) {
printf("istep = %d", wrt_step);
printf("; (u,v) = %.5f, %.5f",
U(hist_x,hist_y), V(hist_x,hist_y) );
printf("\n");
}
}
/* ========================================================================= */
void Write_fc (void) {
Write_fields(u,v,"fc.dat");
}
/* ========================================================================= */
void Write_fields (Real *ufield, Real *vfield, char *filename) {
/* Write final condition file */
FILE *fp;
time_t tt1;
int i,j;
time(&tt1);
fp = fopen(filename, "w");
#if FHN
fprintf(fp,"FHN Model Parameters: beta, gamma, eps, Lx, Ly, Dv = ");
fprintf(fp,"%g, %g, %g, %g, %g, %g\n", bet, gam, eps, length_x, length_y, Dv);
#else
fprintf(fp,"Barkley Model Parameters: a, b, 1/eps, Lx, Ly, Dv = ");
fprintf(fp,"%g, %g, %g, %g, %g, %g\n", a_param, b_param, one_o_eps, length_x, length_y, Dv);
#endif
fprintf(fp,"Numerical Parameters: NX, NY, ts, delta = ");
fprintf(fp,"%d, %d, %g, %g \n", NX, NY, ts, delta);
fprintf(fp,"Riding Coordinates: X1, Y1, TH1 = ");
fprintf(fp,"%15.10f, %15.10f, %15.10f\n", X1, Y1, TH1);
fprintf(fp,"File written: %s",ctime(&tt1));
fprintf(fp,"Comments: 2D \n");
fprintf(fp,"\n");
fprintf(fp,"\n");
fprintf(fp,"\n");
fprintf(fp,"\n");
if(binary_write) {
/* Write binary data */
fprintf(fp,"Binary values of u and v follow\n");
for(j=1;j<=NY;j++) {
fwrite(&(ufield[INDEX(1,j)]), sizeof(Real), NX, fp);
}
for(j=1;j<=NY;j++) {
fwrite(&(vfield[INDEX(1,j)]), sizeof(Real), NX, fp);
}
}
else {
/* Write ascii data */
fprintf(fp,"Ascii values of u and v follow\n");
for(j=1;j<=NY;j++) {
for(i=1;i<=NX;i++) {
fprintf(fp, "%g %g\n", (float)ufield[INDEX(i,j)], (float)vfield[INDEX(i,j)]);
}
}
}
fclose(fp);
}
/* ========================================================================= */
void Read_ic (void) {
/* Reads initial condition files given by std filenames if they exist,
* otherwise calls Generate_*ic() to generate new initial condition. */
FILE *fp;
if ((fp=fopen("ic.dat","r"))!=NULL) {
Read_fields(u,v,fp);
} else {
Generate_ic();
X1=Y1=TH1=0.0;
}
}
/* ========================================================================= */
void Read_fields (Real *ufield, Real *vfield, FILE *fp) {
double u_in, v_in;
Real *u_tmp, *v_tmp, rx, ry;
int nx_ic, ny_ic, npts_ic2, index, i_tmp, j_tmp, i, j;
char f_type, dummy;
/* Read nx_ic etc following = sign on second line of file */
NEXT_LINE(fp);
while( (dummy=getc(fp)) != '=');
fscanf(fp, "%d, %d", &nx_ic, &ny_ic);
if (ufield == u) {
/* Read riding coords following = sign on third line of file */
NEXT_LINE(fp);
while( ((dummy=getc(fp)) != '=') && (dummy != '\n'));
if (dummy=='=')
fscanf(fp, "%g, %g, %g", &X1, &Y1, &TH1);
else
X1=Y1=TH1=0.0;
if(verbose) printf("\nRiding is %s, initial coordinates %g, %g, %g\n",
ride_on?"on":"off", X1, Y1, TH1);
/* Skip to 10th line */
for(i=0;i<7;i++) NEXT_LINE(fp);
} else {
/* Skip to 10th line */
for(i=0;i<8;i++) NEXT_LINE(fp);
}
/* and read first character to determine type B(inary) or A(scii) */
f_type = getc(fp); NEXT_LINE(fp);
if ( (f_type !='B') && (f_type !='A') ) {
if(verbose)
printf("\n ic.dat exists but of unrecognized type Binary or Ascii \n");
exit(1);
}
if(verbose) printf("\nReading ic.dat with nx, ny = %d, %d... \n\n",
nx_ic, ny_ic);
npts_ic2 = nx_ic * ny_ic;
/* Allocate temporary memory and read from file */
u_tmp =(Real *) malloc((unsigned)(npts_ic2)*sizeof(Real));
v_tmp =(Real *) malloc((unsigned)(npts_ic2)*sizeof(Real));
if(f_type =='B') {
/* Binary data file */
fread(u_tmp, sizeof(Real), npts_ic2, fp);
fread(v_tmp, sizeof(Real), npts_ic2, fp);
}
else {
/* Ascii data file */
for(index=0;index(NY/2+ini_y)) U(i,j) = d;
}
}
}
/* ========================================================================= */
/***********************************************************************************************/
/* Save to a file a png image via ppm format generated independently on X graphics capability. */
/* Relies on pnmtopng utility from netpbm library to generate the png format. */
/* Uses forks and pipes so only usable under unix. */
/***************************************************/
#define MAXCHAR (255)
static unsigned char Byte(double x) {
int probyte=x*MAXCHAR;
return probyte<0?0:probyte>MAXCHAR?MAXCHAR:probyte;
}
void Make_images (void) {
static int num = 0;
Make_image(num,"",u,v,UMIN,UMAX,VMIN,VMAX);
num++;
}
void Make_image (int num, char *stem, Real *ufield, Real *vfield, Real umin, Real umax, Real vmin, Real vmax) {
char *fifoname=tempnam(NULL,NULL);
pid_t child_pid;
static char cmd[4096];
FILE *fifo;
int i,j;
double scaled_u, scaled_v;
unsigned char *red, *green, *blue;
int status;
red =(unsigned char *)calloc(NX*NY,1);
green=(unsigned char *)calloc(NX*NY,1);
blue =(unsigned char *)calloc(NX*NY,1);
#define cindex(i,j) (((i)-1)*NX+(j)-1)
if (-1==mkfifo(fifoname,0600)) {perror("Make_image could not make a fifo"); return;}
switch (child_pid=fork()) {
/******/
case -1: /* Fork was unsuccessful */
perror("pipeto could not fork");
return;
/******/
case 0: /* This is the child process. */
/* It spawns the pipe that converts fifo output to the image file */
sprintf(cmd,"cat %s | pnmscale -xsize 400 | pnmflip -tb | pnmtopng > %s%06d.png",fifoname,stem,num);
system(cmd);
_exit(0); /* child process does not have any other function so it terminates */
/******/
default: /* this is the parent process */
/* Open the fifo for writing */
if (NULL==(fifo=fopen(fifoname,"w"))) {perror("Make_image could not write to fifo"); return;}
/* Generate the picture pixels */
for(j=1;j