visual.c

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

          /*-------------------------------------*/
          if (FLAG[i][j] < C_F)
             ADVANCE_AT_BOUND(i,j,&x,&y,u,v,U,V,FLAG,delx,dely,delt);
            
          part->next->x = x; part->next->y = y;
         }
      }
   }
}/*End ADVANCE_PARTICLES*/


/*-------------------------------------------------------------------------*/
/* Computation of new particle location of a particle near a no-slip wall, */
/* guaranteeing, that the new position is not in the obstacle cell.        */
/* Here a modified interpolation algorithm is applied, using the fact that */
/* at no-skip walls, the velocity is not only given at the midpoint of the */
/* edge but on the whole edge                                              */
/*-------------------------------------------------------------------------*/
void ADVANCE_AT_BOUND(int i,int j,REAL *x,REAL *y,REAL u,REAL v,
                      REAL **U,REAL **V,int **FLAG,
                      REAL delx,REAL dely,REAL delt)
{
 int ialt,jalt;
 REAL xalt,yalt;
 REAL ul,ur,vo,vu;     
 REAL x1,x2,y1,y2;

 /* get old particle position */
 xalt = (*x)-delt*u;          yalt = (*y)-delt*v;
 ialt = (int)(xalt/delx)+1;   jalt = (int)(yalt/dely)+1; 

 if (i != ialt){      /* compute new x */
   if (FLAG[ialt+1][jalt] < C_F)
      ur = 0.0;      else{
      if (yalt>= (jalt-0.5)*dely)
         if (FLAG[ialt+1][jalt+1] < C_F){
            y2 = jalt *dely;
            ur = U[ialt][jalt]*(y2-yalt)*2.0/dely;
           }
         else{  
            y1 = (jalt-0.5)*dely;
            y2 = (jalt+0.5)*dely;
            ur = (U[ialt][jalt]*(y2-yalt)+U[ialt][jalt+1]*(yalt-y1))/dely;
	   }
      else      
         if (FLAG[ialt+1][jalt-1] < C_F){
            y1 = (jalt-1.0)*dely;
            ur = U[ialt][jalt]*(yalt-y1)*2.0/dely;
	   }
         else{ 
            y1 = (jalt-1.5)*dely;
            y2 = (jalt-0.5)*dely;
            ur = (U[ialt][jalt-1]*(y2-yalt)+U[ialt][jalt]*(yalt-y1))/dely;  
	   }
     }
   if (FLAG[ialt-1][jalt] < C_F)
      ul = 0.0;  
   else{
      if (yalt>= (jalt-0.5)*dely) 
         if (FLAG[ialt-1][jalt+1] < C_F){
            y2 = jalt *dely;
            ul = U[ialt-1][jalt]*(y2-yalt)*2.0/dely;
           }
         else{   
            y1 = (jalt-0.5)*dely;
            y2 = (jalt+0.5)*dely;
            ul = (U[ialt-1][jalt]*(y2-yalt)+U[ialt-1][jalt+1]*(yalt-y1))/dely;
	   }
      else       
         if (FLAG[ialt-1][jalt-1] < C_F){
            y1 = (jalt-1.0)*dely;
            ul = U[ialt-1][jalt]*(yalt-y1)*2.0/dely;
	   }
         else{ 
            y1 = (jalt-1.5)*dely;
            y2 = (jalt-0.5)*dely;
            ul = (U[ialt-1][jalt-1]*(y2-yalt)+U[ialt-1][jalt]*(yalt-y1))/dely;
	   }
     }
   u = (ul*(ialt*delx-xalt)+ur*(xalt-(ialt-1)*delx))/delx;
   (*x) = xalt+u*delt;
  }      /* end new x */

 if (j != jalt){        /* copute new y */
   if (FLAG[ialt][jalt+1] < C_F)
      vo = 0.0;   
   else{
      if (xalt>= (ialt-0.5)*delx) 
         if (FLAG[ialt+1][jalt+1] < C_F){
            x2 = ialt*delx;
            vo = V[ialt][jalt]*(x2-xalt)*2.0/delx;
           }
         else{  
            x1 = (ialt-0.5)*delx;
            x2 = (ialt+0.5)*delx;
            vo = (V[ialt][jalt]*(x2-xalt)+V[ialt+1][jalt]*(xalt-x1))/delx;
	   }
      else      
         if (FLAG[ialt-1][jalt+1] < C_F){
            x1 = (ialt-1.0)*delx;
            vo = V[ialt][jalt]*(xalt-x1)*2.0/delx;
	   }
         else{   
            x1 = (ialt-1.5)*delx;
            x2 = (ialt-0.5)*delx;
            vo = (V[ialt-1][jalt]*(x2-xalt)+V[ialt][jalt]*(xalt-x1))/delx;
	   }
     }
   if (FLAG[ialt][jalt-1] < C_F)
      vu = 0.0;  
   else{
      if (xalt>= (ialt-0.5)*delx) 
         if (FLAG[ialt+1][jalt-1] < C_F){
            x2 = ialt*delx;
            vu = V[ialt][jalt-1]*(x2-xalt)*2.0/delx;
           }
         else{   
            x1 = (ialt-0.5)*delx;
            x2 = (ialt+0.5)*delx;
            vu = (V[ialt][jalt-1]*(x2-xalt)+V[ialt+1][jalt-1]*(xalt-x1))/delx;
	   }
      else       
         if (FLAG[ialt-1][jalt-1] < C_F){
            x1 = (ialt-1.0)*delx;
            vu = V[ialt][jalt-1]*(xalt-x1)*2.0/delx;
	   }
         else{ 
            x1 = (ialt-1.5)*delx;
            x2 = (ialt-0.5)*delx;
            vu = (V[ialt-1][jalt-1]*(x2-xalt)+V[ialt][jalt-1]*(xalt-x1))/delx;
	   }
     }
   v = (vu*(jalt*dely-yalt)+vo*(yalt-(jalt-1)*dely))/dely;
   (*y) = yalt+v*delt;
  }    /* end new y */
}/* End ADVANCE_AT_BOUND */


/*--------------------------------------------------------------------*/
/* Injection of new particles for streaklines                         */
/*--------------------------------------------------------------------*/
void INJECT_PARTICLES(int N, struct particleline *Partlines){
  int i;
  struct particle *part;

  for(i=1; i<=N; i++){
    part = PARTALLOC(Partlines[i].Particles->x,Partlines[i].Particles->y);
    part->next = Partlines[i].Particles->next;
    Partlines[i].Particles->next = part;
    Partlines[i].length++;       
  }
}/*End INJECT_PARTICLES*/

/*--------------------------------------------------------------*/
/* Append particle positions to file "partfile" in ascii format */
/*--------------------------------------------------------------*/
void WRITE_PARTICLES(char *partfile, int N, struct particleline *Partlines){
  int i;
  FILE *fp;
  struct particle *part;

  fp = fopen(partfile,"ab");

  for(i=1; i<=N; i++){
    fprintf(fp,"%d\n",Partlines[i].length);
    for(part=Partlines[i].Particles; part->next != NULL; part=part->next)
      fprintf(fp,"%3.3f %3.3f\n", part->next->x, part->next->y);
  }

  fclose(fp);
}/*End  WRITE_PARTICLES*/


/*---------------------------------------------------------------------*/
/* Append particle positions to file "partfile" in binary format       */
/*---------------------------------------------------------------------*/
void WRITE_PARTICLES_bin(char *partfile, int N, struct particleline *Partlines)
{
 int i;
 FILE *fp;
 float temp, temp2[2];
 struct particle *part;

 fp = fopen(partfile, "ab");
 for(i=1; i<=N; i++){
    temp=Partlines[i].length;
    fwrite(&temp, sizeof(float), 1, fp);
    part=Partlines[i].Particles;
    for(; part->next != NULL; part=part->next){
      temp2[0]=part->next->x;
      temp2[1]=part->next->y;
      fwrite(temp2, sizeof(float), 2, fp);
    }
 }
 fclose(fp);
}

/*----------------------------------------------------------------------*/
/* Move particle positions and append them to a file if wanted          */
/*----------------------------------------------------------------------*/
void PARTICLE_TRACING(char* outputfile,REAL t,int imax,int jmax,
		      REAL delx,REAL dely,REAL delt,
                      REAL **U,REAL **V,int **FLAG,
                      int N, struct particleline *Partlines, int write)
{
  FILE *fp;

  if(t == 0){
    fp = fopen(outputfile, "wb");
    fprintf(fp,"%d\n%d\n%f\n%f\n%d\n", imax, jmax, delx, dely, N);
    fclose(fp);
    WRITE_PARTICLES(outputfile,N,Partlines);
  }

  ADVANCE_PARTICLES(imax,jmax,delx,dely,delt,U,V,FLAG,N,Partlines);

  if(write & 1)
    WRITE_PARTICLES(outputfile,N,Partlines);
    
}/*End PARTICLE_TRACING*/

/*---------------------------------------------------------------------*/
/* Move particles for streaklines, inject and write particle positions */
/*---------------------------------------------------------------------*/
void STREAKLINES(char* streakfile,int write,
                 int imax, int jmax, REAL delx, REAL dely,REAL delt,REAL t, 
                 REAL **U,REAL **V,int **FLAG,
                 int N, struct particleline *Partlines)
{
  FILE *fp;

  if(t==0){
    fp = fopen(streakfile, "wb");
    fprintf(fp, "%d\n", imax);
    fprintf(fp, "%d\n", jmax);
    fprintf(fp, "%g\n", delx);
    fprintf(fp, "%g\n", dely);
    fprintf(fp, "%d\n", N);
    fclose(fp);
    WRITE_PARTICLES_bin(streakfile,N,Partlines);
   }

  ADVANCE_PARTICLES(imax,jmax,delx,dely,delt,U,V,FLAG,N,Partlines);

  if(write & 2)  INJECT_PARTICLES(N,Partlines);
  if(write & 4)  WRITE_PARTICLES_bin(streakfile,N,Partlines);

}/*End STREAKLINES*/