visual.c

The 2D CFD Program NaSt2D The program is a 2D solver for the incompressible, transient Navier-Sto

#include <stdio.h>
#include <stdlib.h>
#include "datadef.h"
#include "visual.h"

/*---------------------------------------------------------------------------*/
/* Writing U,V,P,PSI, and ZETA into "vecfile" for visualization              */
/*---------------------------------------------------------------------------*/
void OUTPUTVEC_bin(REAL **U,REAL **V,REAL **P,REAL **TEMP,
		   REAL **PSI,REAL **ZETA,REAL **HEAT,int **FLAG,
                   REAL xlength,REAL ylength,int imax,int jmax,
                   char* vecfile)
{
 int i,j;
 float temp;
 FILE *fp;
 fp = fopen(vecfile, "wb");

 temp = xlength;
 fwrite(&temp, sizeof(float), 1, fp);
 temp = ylength;
 fwrite(&temp, sizeof(float), 1, fp);
 temp = imax;
 fwrite(&temp, sizeof(float), 1, fp);
 temp = jmax;
 fwrite(&temp, sizeof(float), 1, fp);

 for(j=1;j<=jmax;j+=1)
  for(i=1;i<=imax;i+=1){
   if( (FLAG[i][j] & C_F) && (FLAG[i][j] < C_E) ) 
      temp = (U[i][j]+U[i-1][j])/2.0;
   else			 temp = 0.0;
   fwrite(&temp, sizeof(float), 1, fp);
  }
 for(j=1;j<=jmax;j+=1)
  for(i=1;i<=imax;i+=1){
   if( (FLAG[i][j] & C_F) && (FLAG[i][j] < C_E) ) 
      temp = (V[i][j]+V[i][j-1])/2.0;
   else			 temp = 0.0;
   fwrite(&temp, sizeof(float), 1, fp);
  }
 for(j=1;j<=jmax;j+=1)
  for(i=1;i<=imax;i+=1){
   if( (FLAG[i][j] & C_F) && (FLAG[i][j] < C_E) ) 
      temp = P[i][j];
   else			 temp = 0.0;
   fwrite(&temp, sizeof(float), 1, fp);
  }
 for(j=1;j<=jmax;j+=1)
  for(i=1;i<=imax;i+=1){
   if( (FLAG[i][j] & C_F) && (FLAG[i][j] < C_E) ) 
        temp = TEMP[i][j];
   else			 temp = -0.5;

   temp = TEMP[i][j];
   fwrite(&temp, sizeof(float), 1, fp);
  }
 for(j=1;j<=jmax-1;j+=1)
  for(i=1;i<=imax-1;i+=1){
    temp = ZETA[i][j];
    fwrite(&temp, sizeof(float), 1, fp);
  }
 for(j=0;j<=jmax;j+=1)
  for(i=0;i<=imax;i+=1){
    temp = PSI[i][j];
    fwrite(&temp, sizeof(float), 1, fp);
  }
 for(j=0;j<=jmax;j+=1)
  for(i=0;i<=imax;i+=1){
    temp = HEAT[i][j];
    fwrite(&temp, sizeof(float), 1, fp);
  }

 fclose(fp);
}

/*-----------------------------------------------------*/
/* Computation of stream function and vorticity        */
/*-----------------------------------------------------*/
void COMPPSIZETA(REAL **U,REAL **V,REAL **PSI,REAL **ZETA,int **FLAG,
                int imax,int jmax,REAL delx,REAL dely)
{
 int i,j;

 /* Computation of the vorticity zeta at the upper right corner     */
 /* of cell (i,j) (only if the corner is surrounded by fluid cells) */
 /*-----------------------------------------------------------------*/
 for (i=1;i<=imax-1;i++)
    for (j=1;j<=jmax-1;j++)
        if( ((FLAG[i][j] & C_F) && (FLAG[i][j] < C_E))  &&
            ((FLAG[i+1][j] & C_F) && (FLAG[i+1][j] < C_E))  &&
            ((FLAG[i][j+1] & C_F) && (FLAG[i][j+1] < C_E))  &&
            ((FLAG[i+1][j+1] & C_F) && (FLAG[i+1][j+1] < C_E)) )
          ZETA[i][j] = (U[i][j+1]-U[i][j])/dely - (V[i+1][j]-V[i][j])/delx;
       else
          ZETA[i][j] = 0.0;

 /* Computation of the stream function at the upper right corner    */
 /* of cell (i,j) (only if bother lower cells are fluid cells)      */
 /*-----------------------------------------------------------------*/
 for (i=0;i<=imax;i++)
   {
    PSI[i][0] = 0.0;
    for(j=1;j<=jmax;j++)
        if( ((FLAG[i][j] & C_F) && (FLAG[i][j] < C_E))  ||
            ((FLAG[i+1][j] & C_F) && (FLAG[i+1][j] < C_E)) )
          PSI[i][j] = PSI[i][j-1] + U[i][j]*dely;
       else
          PSI[i][j] = PSI[i][j-1];
  } 
}

/*-----------------------------------------------------*/
/* Computation of the heat function                    */
/*-----------------------------------------------------*/
void COMP_HEAT(REAL **U,REAL **V,REAL **TEMP,REAL **HEAT,int **FLAG,
               REAL Re,REAL Pr,int imax,int jmax,REAL delx,REAL dely)
{
 int i,j;

 /* Computation at the upper right corner of cell (i,j) */
 /*-----------------------------------------------------*/
 for (i=0;i<=imax;i++)
   {
    HEAT[i][0] = 0.0;
    for(j=1;j<=jmax;j++)
       if( ((FLAG[i][j] & C_F) && (FLAG[i][j] < C_E))  ||
           ((FLAG[i+1][j] & C_F) && (FLAG[i+1][j] < C_E)) )
          HEAT[i][j] = HEAT[i][j-1] + 
                    dely*(U[i][j]*0.5*(1.0+TEMP[i+1][j]+TEMP[i][j])*Re*Pr-
                          (TEMP[i+1][j]-TEMP[i][j])/delx );
   }
}

/*-----------------------------------------------------------------*/
/* Allocate memory for a particle and set the coordinates to (x,y) */
/*-----------------------------------------------------------------*/
struct particle *PARTALLOC(REAL x, REAL y){
  struct particle *part;

  if((part=(struct particle *)malloc(sizeof(struct particle))) == NULL)
    {	
     printf("no memory\n");
     exit(0);
    }
  part->x = x; part->y = y;
  part->next = NULL;
  return( part );
}

/*------------------------------------------------------------------*/
/* Allocate memory for a particleline and set coordinates where     */
/* particles are injected.                                          */
/*------------------------------------------------------------------*/
struct particleline *SET_PARTICLES(int N,REAL pos1x,REAL pos1y,
				   	 REAL pos2x,REAL pos2y){
   int i;
   REAL hx,hy;
   struct particleline *Partlines;

   if((Partlines=(struct particleline *)
                 malloc((unsigned)(N) * sizeof(struct particleline))) == NULL)
     {	
      printf("no memory\n");
      exit(0);
     }
   Partlines -= 1;  

   if (N>=2)
     {
      hx  = (pos2x-pos1x)/(N-1);
      hy  = (pos2y-pos1y)/(N-1);
      for(i=1; i<=N; i++){
         Partlines[i].length = 0;
         Partlines[i].Particles = 
             PARTALLOC(pos1x+hx*(i-1),pos1y+hy*(i-1));
         Partlines[i].Particles->next = 
             PARTALLOC(pos1x+hx*(i-1),pos1y+hy*(i-1));
         Partlines[i].length++;
       }
     }
   return(Partlines);
}/*End SET_PARTICLES*/

/*-------------------------------------------------------------*/
/* Moving particles                                            */
/*-------------------------------------------------------------*/
void ADVANCE_PARTICLES(int imax,int jmax,REAL delx,REAL dely,REAL delt,
                       REAL **U,REAL **V,int **FLAG,
                       int N,struct particleline *Partlines)
{
  int i, j, k;
  REAL x, y, x1, y1, x2, y2, u, v; 
  struct particle *part,*help;

  for(k=1;k<=N;k++){
     for(part=Partlines[k].Particles; part->next != NULL; part=part->next){
        /* first element is only a dummy element           */
        /*-------------------------------------------------*/
        x = part->next->x; y = part->next->y;

        /* Computation of new x-coordinates by discretizing dx/dt=u */
        /*----------------------------------------------------------*/
        i = (int)(x/delx)+1;  j = (int)((y+0.5*dely)/dely)+1;

        x1 = (i-1)*delx;    y1 = ((j-1)-0.5)*dely;
        x2 = i*delx;        y2 = (j-0.5)*dely;

        /* bilinear interpolation */
        /*------------------------*/
        u = ((x2-x)*(y2-y)*U[i-1][j-1] +
  	     (x-x1)*(y2-y)*U[i][j-1]   +
	     (x2-x)*(y-y1)*U[i-1][j]   +
	     (x-x1)*(y-y1)*U[i][j])/delx/dely;

        /* Computation of new y-coordinates by discretizing dy/dt=v */
        /*----------------------------------------------------------*/
        i = (int)((x+0.5*delx)/delx)+1; j = (int)(y/dely)+1;

        x1 = ((i-1)-0.5)*delx;    y1 = (j-1)*dely;
        x2 = (i-0.5)*delx;        y2 = j*dely;

        /* bilinear interpolation */
        /*------------------------*/
        v = ((x2-x)*(y2-y)*V[i-1][j-1] +
	     (x-x1)*(y2-y)*V[i][j-1]   +
	     (x2-x)*(y-y1)*V[i-1][j]   +
	     (x-x1)*(y-y1)*V[i][j])/delx/dely;

        x += delt*u;   y += delt*v; 

        /*-------------------------------------*/
        /* determine new cell for the particle */
        /*-------------------------------------*/
        i = (int)(x/delx)+1;   j = (int)(y/dely)+1;

        /* if particle left the fluid domain, delete it */
        /*----------------------------------------------*/
        if (x>=imax*delx || y>=jmax*dely || x<=0 || y<=0){
          help = part->next->next;
          free(part->next);
          part->next = help;
          Partlines[k].length--;
	  if (help == NULL)
	    break; 
        }
        else{
          /*-------------------------------------*/
          /* special treatment if particle would */
	  /* be in an inner obstacle cell        */