init.c

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

  for(i=1;i<=imax;i++)
     for(j=1;j<=jmax;j++)
        FLAG[i][j] = C_F;

                   /* problem dependent obstacle cells in the interior */
                   /*--------------------------------------------------*/
  if(strcmp(problem,"fluidtrap")==0)
    {
     for(i=9*imax/22+1;i<=13*imax/22;i++)
       {
	for(j=1;j<=4*jmax/11;j++)
           FLAG[i][j] = C_B;
        for(j=8*jmax/11+1;j<=jmax;j++)
           FLAG[i][j] = C_B;          
       }
    }

  if(strcmp(problem,"plate")==0)
    {
                   /* flow past an inclined plate */
                   /*---------------------------------------*/
     low = 2*jmax/5;          /* lower and upper bound of the plate */
     up  = 3*jmax/5;
     FLAG[low][low]     = C_B;
     FLAG[low][low+1]   = C_B;
     FLAG[up][up-1]     = C_B;
     FLAG[up][up]       = C_B;
     for (i=low+1;i<=up-1;i++)
        for (j=i-1;j<=i+1;j++)
           FLAG[i][j] = C_B;
    }

  if(strcmp(problem,"backstep")==0 || strcmp(problem,"wave")==0)
    {
                    /* flow past a backward facing step */
                    /*----------------------------------*/
     for (i=1;i<=jmax;i++)
        for (j=1;j<=jmax/2;j++)
           FLAG[i][j] = C_B;
    }

  if(strcmp(problem,"circle")==0)
    {
                   /* flow past a cylinder/circle */
                   /*-----------------------------*/  
     mx = 20.0/41.0*jmax*dely;
     my = mx;
     rad1 = 5.0/41.0*jmax*dely;
     for (i=1;i<=imax;i++)
        for (j=1;j<=jmax;j++)   
          {
           x = (i-0.5)*delx;
           y = (j-0.5)*dely;
           if ((x-mx)*(x-mx)+(y-my)*(y-my) <= rad1*rad1)
              FLAG[i][j] = C_B;
	  }
    }

  if(strcmp(problem,"molding")==0)
    {
                   /* circular obstacle */
                   /*-------------------*/  
     mx = jmax*dely/2;
     my = jmax*dely/2;
     rad1 = jmax*dely/6;
     for (i=1;i<=imax;i++)
        for (j=1;j<=jmax;j++)   
          {
           x = (i-0.5)*delx;
           y = (j-0.5)*dely;
           if ((x-mx)*(x-mx)+(y-my)*(y-my) <= rad1*rad1)
              FLAG[i][j] = C_B;
	  }
    }

                     /* Printing the geometry of the fluid domain */
                     /*-------------------------------------------*/
  printf ("\nGeometry of the fluid domain:\n\n");
  for(j=jmax+1;j>=0;j--)
    {
     for(i=0;i<=imax+1;i++)
        if (!(FLAG[i][j] & C_F))
           printf("**");
        else      
           printf("  ");                                    
     printf ("\n");
    }
  printf ("\n");
  printf ("\n");
                    /* FLAGs for boundary cells */
                    /*--------------------------*/
  (*ibound) = 0;
  for(i=1;i<=imax;i++)
     for(j=1;j<=jmax;j++){
        if (!(FLAG[i][j] & C_F))
           (*ibound)++;
        FLAG[i][j] += ((FLAG[i-1][j] & C_F)*B_W + (FLAG[i+1][j] & C_F)*B_O +
                      (FLAG[i][j-1] & C_F)*B_S + (FLAG[i][j+1] & C_F)*B_N)/C_F;
        switch (FLAG[i][j]){
           case 0x0003:
           case 0x0007:
           case 0x000b:
           case 0x000c:
           case 0x000d:
           case 0x000e:
           case 0x000f:{           
                     printf("Illegal obstacle cell [%d][%d]\n",i,j);
                     exit(0);
                    }  
	 }
      }
}


/*-------------------------------------------------------------------------*/
/* Writing the values of U,V,P,TEMP,FLAG into a file for subsequent        */ 
/* calculations                                                            */
/*-------------------------------------------------------------------------*/
void WRITE_bin(REAL **U,REAL **V,REAL **P,REAL **TEMP,int **FLAG,
		      int imax,int jmax,char* file)
{
 int i;
 FILE *fp;

 fp = fopen(file, "wb"); 

 fwrite(&imax, sizeof(int), 1, fp);
 fwrite(&jmax, sizeof(int), 1, fp);

 for(i=0;i<=imax+1;i+=1)
   fwrite(U[i], sizeof(REAL), jmax+2, fp);
 for(i=0;i<=imax+1;i+=1)
   fwrite(V[i], sizeof(REAL), jmax+2, fp);
 for(i=0;i<=imax+1;i+=1)
   fwrite(P[i], sizeof(REAL), jmax+2, fp);
 for(i=0;i<=imax+1;i+=1)
   fwrite(TEMP[i], sizeof(REAL), jmax+2, fp);
 for(i=0;i<=imax+1;i+=1)
   fwrite(FLAG[i], sizeof(int), jmax+2, fp);
 fclose(fp);
}

/*-------------------------------------------------------------------------*/
/* Reading initial values from a file                                      */
/*-------------------------------------------------------------------------*/
int READ_bin(REAL **U,REAL **V,REAL **P,REAL **TEMP,int **FLAG,
		 int imax,int jmax,char* file)
{
 int i,j;
 FILE *fp;

 if(strcmp(file, "none") == 0) return(-1);

 if( (fp = fopen(file,"rb")) == NULL){
   printf("Error in READ_bin: File %s cannot be opened!\n", file);
   return(1);
 }

 fread(&i, sizeof(int), 1, fp);
 fread(&j, sizeof(int), 1, fp);

 if(i!=imax || j!=jmax){
    printf("ATTENTION: imax and jmax have wrong values in %s\n",file);
    printf("imax = %d  jmax = %d\n", i, j);
    return(2);
 }

 for(i=0;i<=imax+1;i+=1)
   fread(U[i], sizeof(REAL), jmax+2, fp);
 for(i=0;i<=imax+1;i+=1)
   fread(V[i], sizeof(REAL), jmax+2, fp);
 for(i=0;i<=imax+1;i+=1)
   fread(P[i], sizeof(REAL), jmax+2, fp);
 for(i=0;i<=imax+1;i+=1)
   fread(TEMP[i], sizeof(REAL), jmax+2, fp);
 for(i=0;i<=imax+1;i+=1)
   fread(FLAG[i], sizeof(int), jmax+2, fp);
 fclose(fp);

 return(0);
}


/*-------------------------------------------------------------------------*/
/* RMATRIX allocates memory for a [nrl,nrh]x[ncl,nch]-array of REAL-type   */
/*-------------------------------------------------------------------------*/
REAL **RMATRIX(int nrl,int nrh,int ncl,int nch)
{
  int i;
  REAL **m;
  if((m = (REAL**) malloc((unsigned) (nrh-nrl+1)*sizeof(double*))) == NULL)
     {	
      printf("no memory\n");
      exit(0);
     }
  m -= nrl;
  for(i=nrl;i<=nrh;i++)
    {
     if((m[i] = (REAL*) malloc((unsigned) (nch-ncl+1)*sizeof(double)))==NULL)
       {	
        printf("no memory\n");
        exit(0);
       }
     m[i] -= ncl;
    }
  return m;
} 

/*-------------------------------------------------------------------------*/
/* FREE_RMATRIX frees the memory of an array allocated with RMATRIX        */
/*-------------------------------------------------------------------------*/
void FREE_RMATRIX(REAL** m,int nrl,int nrh,int ncl,int nch)
{
  int i;
  for (i=nrh;i>=nrl;i--) free((void*) (m[i]+ncl));
  free((char*) (m+nrl));
}

/*-------------------------------------------------------------------------*/
/* IMATRIX allocates memory for a [nrl,nrh]x[ncl,nch]-array of integer-type*/
/*-------------------------------------------------------------------------*/
int **IMATRIX(int nrl,int nrh,int ncl,int nch)
{
  int i;
  int **m;
  if((m = (int**) malloc((unsigned) (nrh-nrl+1)*sizeof(int*))) == NULL)
     {	
      printf("no memory\n");
      exit(0);
     }
  m -= nrl;
  for(i=nrl;i<=nrh;i++)
    {
     if((m[i] = (int*) malloc((unsigned) (nch-ncl+1)*sizeof(int)))==NULL)
       {	
        printf("no memory\n");
        exit(0);
       }
     m[i] -= ncl;
    }
  return m;
} 

/*-------------------------------------------------------------------------*/
/* FREE_IMATRIX frees the memory of an array allocated with IMATRIX        */
/*-------------------------------------------------------------------------*/
void FREE_IMATRIX(int** m,int nrl,int nrh,int ncl,int nch)
{
  int i;
  for (i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
  free((char*) (m+nrl));
}