lbio.c

基于格子Boltzmann方法开源可视化软件的源代码 具有很高的实用价值。对学习LBM方法及其软件开发非常游泳

//##############################################################################
//
// lbio.c
//
//  - Lattice Boltzmann I/O routines.
//
//  - Mainly, dump the data to files that can be read by Matlab.
//
//  - Also, output routines for facilitating debugging.
//
//  - Should have a routine that dumps a matlab script?
//

// Some compilers, e.g., VC++, don't have the usual round() function
// in their math library.  Alternatively, ROUND can be defined as
// ceil or floor or some other rounding function.  It is used in
// the below routines for converting the real number valued of 
// quantities at a lattice node into integer RGB values for writing 
// to BMP files.
#define ROUND floor

#if SWAP_BYTE_ORDER
// Swap byte order.
#define ENDIAN2(w) ((((w)&0x00ff)<<8)|(((w)&0xff00)>>8))
#define ENDIAN4(w) ((((w)&0x000000ff)<<24)|(((w)&0xff000000)>>24)|(((w)&0x0000ff00)<<8)|(((w)&0x00ff0000)>>8))
#else /* !( SWAP_BYTE_ORDER) */
#define ENDIAN2(w) (w)
#define ENDIAN4(w) (w)
#endif /* SWAP_BYTE_ORDER */

//void output_frame( lattice_ptr lattice)
//##############################################################################
//
// O U T P U T   F R A M E
//
void output_frame( lattice_ptr lattice)
{
  double s, u[2];
  double nu;
  double L;
#if VERBOSITY_LEVEL > 0
      printf("\n");
      printf( "========================================"
              "========================================\n");
      printf("Begin file I/O at time = %d, frame = %d.\n", 
          lattice->time, lattice->time/lattice->param.FrameRate);
      printf("\n");
#endif /* VERBOSITY_LEVEL > 0 */
      dump_frame_summary( lattice);
#if WRITE_MACRO_VAR_DAT_FILES
      dump_macro_vars( lattice, lattice->time);
#endif /* WRITE_MACRO_VAR_DAT_FILES */
#if WRITE_PDF_DAT_FILES
      dump_pdf( lattice, lattice->time);
#endif /* WRITE_PDF_DAT_FILES */
      if( lattice->param.dump_rho) { rho2bmp( lattice, lattice->time);}
      if( lattice->param.dump_u  ) { u2bmp( lattice, lattice->time);}
      if( lattice->param.dump_vor) { vor2bmp( lattice, lattice->time);}
#if NON_LOCAL_FORCES
      if( lattice->param.G != 0.) 
      { 
        if( lattice->param.dump_force) { force2bmp( lattice);}
      }
      if(    lattice->param.Gads[0] != 0.
          || lattice->param.Gads[0] != 0.)
      {
        if( lattice->param.dump_force) { sforce2bmp( lattice);}
      }
#endif /* NON_LOCAL_FORCES */
      slice( lattice);
#if WRITE_CHEN_DAT_FILES
      chen_output( lattice);
#endif /* WRITE_CHEN_DAT_FILES */
#if VERBOSITY_LEVEL > 0
      printf("\n");
      printf("File I/O done.\n");
      printf("--\n");
#endif /* VERBOSITY_LEVEL > 0 */

  nu = (1./3.)*(lattice->param.tau[0] - .5);
  L = lattice->param.length_scale;

  compute_ave_u( lattice, u, 0);
  s = sqrt( u[0]*u[0] + u[1]*u[1]);
  printf("subs 0: Re = ux_ave*L/nu = %f * %f / %f = %f\n",
    u[0], L, nu, u[0]*L/nu );
  printf("subs 0: Re = uy_ave*L/nu = %f * %f / %f = %f\n",
    u[1], L, nu, u[1]*L/nu );
  printf("subs 0: Re = u_ave*L/nu  = %f * %f / %f = %f\n",
    s, L, nu, s*L/nu );

#if NUM_FLUID_COMPONENTS == 2
  compute_ave_u( lattice, u, 1);
  s = sqrt( u[0]*u[0] + u[1]*u[1]);
  printf("subs 1: Re = ux_ave*L/nu = %f * %f / %f = %f\n",
    u[0], L, nu, u[0]*L/nu );
  printf("subs 1: Re = uy_ave*L/nu = %f * %f / %f = %f\n",
    u[1], L, nu, u[1]*L/nu );
  printf("subs 1: Re = u_ave*L/nu  = %f * %f / %f = %f\n",
    s, L, nu, s*L/nu );
#endif /* NUM_FLUID_COMPONENTS == 2 */

#if STORE_UEQ
  compute_ave_ueq( lattice, u);
  s = sqrt( u[0]*u[0] + u[1]*u[1]);
  printf("eq:     Re = ux_ave*L/nu = %f * %f / %f = %f\n",
    u[0], L, nu, u[0]*L/nu );
  printf("eq:     Re = uy_ave*L/nu = %f * %f / %f = %f\n",
    u[1], L, nu, u[1]*L/nu );
  printf("eq:     Re = u_ave*L/nu  = %f * %f / %f = %f\n",
    s, L, nu, s*L/nu );
#endif /* STORE_UEQ */

} /* void output_frame( lattice_ptr lattice) */

// void dump_frame_info( struct lattice_struct *lattice)
//##############################################################################
// 
// D U M P   F R A M E   I N F O
//
void dump_frame_summary( struct lattice_struct *lattice)
{
  char   filename[1024];
  FILE   *o;
  double min_u[5], max_u[5], ave_u[5], flux[3];
  double min_rho,  max_rho,  ave_rho;
  double rho_ratio, u_x_ratio, u_y_ratio;
  int    subs;

 for( subs = 0; subs < NUM_FLUID_COMPONENTS; subs++)
 {
  sprintf( filename, "./out/frames%dx%d_subs%02d.dat", 
    lattice->param.LX, 
    lattice->param.LY,
    subs);
  // On the first timestep, make sure we start with a new file.
  if( lattice->time==0)
  {
      if( !( o = fopen(filename,"w+")))
      {
        printf("ERROR: fopen(\"%s\",\"w+\") = NULL.  Bye, bye!\n", filename);
        exit(1);
      }
      else
      {
        // Put a header on the file.
        fprintf( o, "\n");
        fprintf( o, "       time  "
                    "        |j|  "
                    "        j_x  "
                    "        j_y  "
                    "    ave |u|  "
                    "   ave |u_x| "
                    "   ave |u_y| "
                    "    ave u_x  "
                    "    ave u_y  "
                    "    min |u|  "
                    "   min |u_x| "
                    "   min |u_y| "
                    "    min u_x  "
                    "    min u_y  "
                    "    max |u|  "
                    "   max |u_x| "
                    "   max |u_y| "
                    "    max u_x  "
                    "    max u_y  "
                    "  max/ave_x  "
                    "  max/ave_y  "
                    "    min rho  "
                    "    max rho  "
                    "    ave rho  "
                    "    max/ave  "
                    "\n");
        fprintf( o, " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    " ------------"
                    "\n");
        fclose(o);
      }
  }
  if( !( o = fopen(filename,"a+")))
  {
    printf("ERROR: fopen(\"%s\",\"a+\") = NULL.  Bye, bye!\n", filename);
    exit(1);
  }
  compute_min_u_all( lattice, min_u, subs);
  compute_max_u_all( lattice, max_u, subs);
  compute_ave_u_all( lattice, ave_u, subs);
  compute_min_rho( lattice, &min_rho, subs);
  compute_max_rho( lattice, &max_rho, subs);
  compute_ave_rho( lattice, &ave_rho, subs);
  rho_ratio = ( ave_rho  != 0.) ? ( max_rho /ave_rho ):( 1.);
  u_x_ratio = ( ave_u[1] != 0.) ? ( max_u[1]/ave_u[0]):( 1.);
  u_y_ratio = ( ave_u[2] != 0.) ? ( max_u[2]/ave_u[1]):( 1.);

  compute_flux( lattice, flux, subs);

  fprintf( o, 
    "%12d "
    "%12.7f %12.7f %12.7f %12.7f %12.7f %12.7f %12.7f %12.7f %12.7f "
    "%12.7f %12.7f %12.7f %12.7f %12.7f %12.7f %12.7f %12.7f %12.7f "
    "%12.7f %12.7f %12.7f %12.7f %12.7f %12.7f\n", 
    lattice->time, 
    flux [0], flux [1], flux [2], 
    ave_u[0], ave_u[1], ave_u[2], ave_u[3], ave_u[4], 
    min_u[0], min_u[1], min_u[2], min_u[3], min_u[4], 
    max_u[0], max_u[1], max_u[2], max_u[3], max_u[4], 
    (u_x_ratio<=9999.)?(u_x_ratio):(9999.), 
    (u_y_ratio<=9999.)?(u_y_ratio):(9999.), 
    min_rho,
    max_rho,
    ave_rho,
    (rho_ratio<=9999.)?(rho_ratio):(9999.) );

  fclose(o);

#if VERBOSITY_LEVEL > 0
  printf("dump_frame_info() -- Wrote file \"%s\"\n", filename);
#endif /* VERBOSITY_LEVEL > 0 */

#if VERBOSITY_LEVEL > 0
  printf("dump_frame_info() -- frame = %d/%d = %d\n", 
      lattice->time,
      lattice->param.FrameRate,
      (int)((double)lattice->time/(double)lattice->param.FrameRate));
#endif /* VERBOSITY_LEVEL > 0 */
 }

} /* void dump_frame_info( struct lattice_struct *lattice) */

// void dump_macro_vars( struct lattice_struct *lattice)
//##############################################################################
// 
// D U M P   M A C R O S C O P I C 
//
//  - Output the macro_vars variables to files.
//
void dump_macro_vars( struct lattice_struct *lattice, int time)
{
  char   filename[1024];
  FILE   *o, *o_u, *o_rho, *o_ux, *o_uy, *o_ueq, *o_ueq_x, *o_ueq_y;
  int    *node_ptr;
  int    n;
  double *macro_vars_ptr;
  double *ueq;
  int    frame;
#if WRITE_RHO_AND_U_TO_TXT
  int    i, j;
#endif /* WRITE_RHO_AND_U_TO_TXT */
  double min_u[2], max_u[2],  ave_u[2];
  double min_rho, max_rho,   ave_rho;
  double rho_ratio, u_x_ratio, u_y_ratio;
  int    subs;

  frame = (int)((double)lattice->time/(double)lattice->param.FrameRate);

 for( subs = 0; subs < NUM_FLUID_COMPONENTS; subs++)
 {

  // W R I T E   R H O   A N D   U 
  //
  //  - Write the density and velocity values at the active nodes to
  //    the rho and u dat files.
  //
  sprintf( filename, "./out/rho_frame%04d_subs%02d.dat", frame, subs);
  if( !( o_rho = fopen( filename, "w+")))
  {
    printf("ERROR: fopen( \"%s\", \"w+\") = NULL.  Bye, bye!\n", filename);
    exit(1);
  }

  sprintf( filename, "./out/u_frame%04d_subs%02d.dat", frame, subs);
  if( !( o_u = fopen( filename, "w+")))
  {
    printf("ERROR: fopen( \"%s\", \"w+\") = NULL.  Bye, bye!\n", filename);
    exit(1);
  }

#if STORE_UEQ
  sprintf( filename, "./out/ueq_frame%04d_subs%02d.dat", frame, subs);
  if( !( o_ueq = fopen( filename, "w+")))
  {
    printf("ERROR: fopen( \"%s\", \"w+\") = NULL.  Bye, bye!\n", filename);
    exit(1);
  }

  ueq = lattice->ueq[0].u;

#endif /* STORE_UEQ */

  macro_vars_ptr = &( lattice->macro_vars[subs][0].rho);

  for( n=0; n<lattice->NumNodes; n++)
  {

    fprintf( o_rho, "%20.17f\n", *macro_vars_ptr++);

    fprintf( o_u, "%20.17f ", *macro_vars_ptr++);

    fprintf( o_u, "%20.17f\n", *macro_vars_ptr++);

#if STORE_UEQ
    fprintf( o_ueq, "%20.17f ", *ueq++);
    fprintf( o_ueq, "%20.17f\n", *ueq++);
#endif /* STORE_UEQ */

  }

  fclose(o_u);
#if VERBOSITY_LEVEL > 0
  sprintf( filename, "./out/u_frame%04d_subs%02d.dat", frame, subs);
#endif /* VERBOSITY_LEVEL > 0 */
  printf("dump_macro_vars() -- Wrote file \"%s\"\n", filename);
  fclose(o_rho);
#if VERBOSITY_LEVEL > 0
  sprintf( filename, "./out/rho_frame%04d_subs%02d.dat", frame, subs);
#endif /* VERBOSITY_LEVEL > 0 */
  printf("dump_macro_vars() -- Wrote file \"%s\"\n", filename);
#if STORE_UEQ
  fclose(o_ueq);
#if VERBOSITY_LEVEL > 0
  sprintf( filename, "./out/ueq_frame%04d_subs%02d.dat", frame, subs);
#endif /* VERBOSITY_LEVEL > 0 */
  printf("dump_macro_vars() -- Wrote file \"%s\"\n", filename);
#endif /* STORE_UEQ */

#if WRITE_RHO_AND_U_TO_TXT
  // NOTE: This is very inefficient.  But it's only intended
  // for debugging purposes on small problems.
  sprintf( filename, "./out/rho_frame%04d_subs%02d.txt", frame, subs);
  if( !( o_rho = fopen( filename, "w+")))
  {
    printf("ERROR: fopen( \"%s\", \"w+\") = NULL.  Bye, bye!\n", filename);
    exit(1);
  }
  sprintf( filename, "./out/ux_frame%04d_subs%02d.txt", frame, subs);
  if( !( o_ux = fopen( filename, "w+")))
  {
    printf("ERROR: fopen( \"%s\", \"w+\") = NULL.  Bye, bye!\n", filename);
    exit(1);
  }
  sprintf( filename, "./out/uy_frame%04d_subs%02d.txt", frame, subs);
  if( !( o_uy = fopen( filename, "w+")))
  {
    printf("ERROR: fopen( \"%s\", \"w+\") = NULL.  Bye, bye!\n", filename);
    exit(1);
  }
#if STORE_UEQ
  sprintf( filename, "./out/ueq_x_frame%04d_subs%02d.txt", frame, subs);
  if( !( o_ueq_x = fopen( filename, "w+")))
  {
    printf("ERROR: fopen( \"%s\", \"w+\") = NULL.  Bye, bye!\n", filename);
    exit(1);
  }
  sprintf( filename, "./out/ueq_y_frame%04d_subs%02d.txt", frame, subs);
  if( !( o_ueq_y = fopen( filename, "w+")))
  {
    printf("ERROR: fopen( \"%s\", \"w+\") = NULL.  Bye, bye!\n", filename);
    exit(1);
  }
#endif /* STORE_UEQ */
  for( j=lattice->param.LY-1; j>=0; j--)
  {
    n = j*lattice->param.LX;
    for( i=0; i<lattice->param.LX; i++, n++)
    {
      fprintf( o_rho, "%12.7f ", lattice->macro_vars[subs][n].rho);
      fprintf( o_ux,  "%12.7f ", lattice->macro_vars[subs][n].u[0]);
      fprintf( o_uy,  "%12.7f ", lattice->macro_vars[subs][n].u[1]);
#if STORE_UEQ
      fprintf( o_ueq_x,  "%12.7f ", lattice->ueq[n].u[0]);
      fprintf( o_ueq_y,  "%12.7f ", lattice->ueq[n].u[1]);
#endif /* STORE_UEQ */
      if( n==lattice->NumNodes)
      {
        fprintf( o_rho, "%12.7f ", 0.);
        fprintf( o_ux, "%12.7f ", 0.);
        fprintf( o_uy, "%12.7f ", 0.);
#if STORE_UEQ
        fprintf( o_ueq_x, "%12.7f ", 0.);
        fprintf( o_ueq_y, "%12.7f ", 0.);
#endif /* STORE_UEQ */
      }
    }
    fprintf( o_rho, "\n");
    fprintf( o_ux, "\n");
    fprintf( o_uy, "\n");
#if STORE_UEQ
    fprintf( o_ueq_x, "\n");
    fprintf( o_ueq_y, "\n");
#endif /* STORE_UEQ */
  }

  fclose(o_ux);
#if VERBOSITY_LEVEL > 0
  sprintf( filename, "./out/ux_frame%04d_subs%02d.txt", frame, subs);
#endif /* VERBOSITY_LEVEL > 0 */
  printf("dump_macro_vars() -- Wrote file \"%s\"\n", filename);
  fclose(o_uy);
#if VERBOSITY_LEVEL > 0
  sprintf( filename, "./out/uy_frame%04d_subs%02d.txt", frame, subs);
#endif /* VERBOSITY_LEVEL > 0 */
  printf("dump_macro_vars() -- Wrote file \"%s\"\n", filename);
  fclose(o_rho);
#if VERBOSITY_LEVEL > 0
  sprintf( filename, "./out/rho_frame%04d_subs%02d.txt", frame, subs);
#endif /* VERBOSITY_LEVEL > 0 */
  printf("dump_macro_vars() -- Wrote file \"%s\"\n", filename);
#if STORE_UEQ
  fclose(o_ueq_x);
#if VERBOSITY_LEVEL > 0
  sprintf( filename, "./out/ueq_x_frame%04d_subs%02d.txt", frame, subs);
#endif /* VERBOSITY_LEVEL > 0 */
  printf("dump_macro_vars() -- Wrote file \"%s\"\n", filename);
  fclose(o_ueq_y);
#if VERBOSITY_LEVEL > 0
  sprintf( filename, "./out/ueq_y_frame%04d_subs%02d.txt", frame, subs);
#endif /* VERBOSITY_LEVEL > 0 */
  printf("dump_macro_vars() -- Wrote file \"%s\"\n", filename);
#endif /* STORE_UEQ */

#endif /* WRITE_RHO_AND_U_TO_TXT */