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?
//

// 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, *oo;
  int    *node_ptr;
  int    n;
  double *macro_vars_ptr;
  int    frame;
#if WRITE_RHO_AND_U_TO_TXT
  int    i, j;
#endif /* WRITE_RHO_AND_U_TO_TXT */
  double max_u[2], ave_u[2];
  double max_rho, ave_rho;
  double rho_ratio, u_x_ratio, u_y_ratio;

  // A P P E N D   T I M E S T E P 
  //
  //  - Append the current timestep to the frames dat file.
  //
  sprintf( filename, "frames%dx%d.dat", 
    lattice->param.LX, 
    lattice->param.LY);
  // On the first timestep, make sure we start with a new file.
  if( 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   max_u[0]   max_u[1]   ave_u[0]   ave_u[1] "
          "max/ave[0] max/ave[1]    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_max_u( lattice, max_u);
  compute_ave_u( lattice, ave_u);
  compute_max_rho( lattice, &max_rho);
  compute_ave_rho( lattice, &ave_rho);
  rho_ratio = ( ave_rho  != 0.) ? ( max_rho /ave_rho ):( 1.);
  u_x_ratio = ( ave_u[0] != 0.) ? ( max_u[0]/ave_u[0]):( 1.);
  u_y_ratio = ( ave_u[1] != 0.) ? ( max_u[1]/ave_u[1]):( 1.);
  fprintf( o, 
    "%10d %10.7f %10.7f %10.7f %10.7f %10.7f %10.7f %10.7f %10.7f %10.7f\n", 
    time, 
    max_u[0], max_u[1], 
    ave_u[0], ave_u[1], 
    u_x_ratio, 
    u_y_ratio, 
    max_rho,
    ave_rho,
    rho_ratio );
  fclose(o);
#if VERBOSITY_LEVEL > 0
  printf("dump_macro_vars() -- Wrote file \"%s\"\n", filename);
#endif /* VERBOSITY_LEVEL > 0 */

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

  // W R I T E   C O O R D I N A T E S   O F   A C T I V E   N O D E S 
  //
  //  - Write the coordinates of the active nodes to the coords dat file.
  //
  if( frame==0)
  {
    sprintf( filename, "coords%04d.dat", frame);
    if( !( o = fopen( filename, "w+")))
    {
      printf("ERROR: fopen( \"%s\", \"w+\") = NULL.  Bye, bye!\n", filename);
      exit(1);
    }

    node_ptr = &(lattice->node[0].i);
    for( n=0; n<lattice->NumNodes; n++)
    {
      fprintf( o, "%d %d\n", *node_ptr++, *node_ptr++);
      node_ptr+=9;
    }
    fclose(o);
#if VERBOSITY_LEVEL > 0
    printf("dump_macro_vars() -- Wrote file \"%s\"\n", filename);
#endif /* VERBOSITY_LEVEL > 0 */
  } /* if( frame==0) */

  // 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, "rho%04d.dat", frame);
  if( !( o = fopen( filename, "w+")))
  {
    printf("ERROR: fopen( \"%s\", \"w+\") = NULL.  Bye, bye!\n", filename);
    exit(1);
  }

  sprintf( filename, "u%04d.dat", frame);
  if( !( oo = fopen( filename, "w+")))
  {
    printf("ERROR: fopen( \"%s\", \"w+\") = NULL.  Bye, bye!\n", filename);
    exit(1);
  }

  macro_vars_ptr = lattice->macro_vars[0].rho;

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

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

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

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

  }

  fclose(oo);
#if VERBOSITY_LEVEL > 0
  sprintf( filename, "u%04d.dat", frame);
#endif /* VERBOSITY_LEVEL > 0 */
  printf("dump_macro_vars() -- Wrote file \"%s\"\n", filename);
  fclose(o);
#if VERBOSITY_LEVEL > 0
  sprintf( filename, "rho%04d.dat", frame);
#endif /* VERBOSITY_LEVEL > 0 */
  printf("dump_macro_vars() -- Wrote file \"%s\"\n", filename);

#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, "rho%04d.txt", frame);
  if( !( o = fopen( filename, "w+")))
  {
    printf("ERROR: fopen( \"%s\", \"w+\") = NULL.  Bye, bye!\n", filename);
    exit(1);
  }
  sprintf( filename, "u%04d.txt", frame);
  if( !( oo = fopen( filename, "w+")))
  {
    printf("ERROR: fopen( \"%s\", \"w+\") = NULL.  Bye, bye!\n", filename);
    exit(1);
  }
  for( j=0; j<lattice->param.LY; j++)
  {
    for( i=0; i<lattice->param.LX; i++)
    {
      for( n=0; n<lattice->NumNodes; n++)
      {
        if(    lattice->node[n].i == i 
            && lattice->node[n].j == j  )
        {
          fprintf( o, "%10.7f ", lattice->macro_vars[n].rho[0]);
          fprintf( oo, "( %10.7f, %10.7f) ", 
            lattice->macro_vars[n].u[0],
            lattice->macro_vars[n].u[1] );
          break;
        }
      }
      if( n==lattice->NumNodes)
      {
        fprintf( o, "%10.7f ", 0.);
        fprintf( oo, "( %10.7f, %10.7f) ", 0., 0.);
      }
    }
    fprintf( o, "\n");
    fprintf( oo, "\n");
  }

  fclose(oo);
#if VERBOSITY_LEVEL > 0
  sprintf( filename, "u%04d.txt", frame);
#endif /* VERBOSITY_LEVEL > 0 */
  printf("dump_macro_vars() -- Wrote file \"%s\"\n", filename);
  fclose(o);
#if VERBOSITY_LEVEL > 0
  sprintf( filename, "rho%04d.txt", frame);
#endif /* VERBOSITY_LEVEL > 0 */
  printf("dump_macro_vars() -- Wrote file \"%s\"\n", filename);

#endif /* WRITE_RHO_AND_U_TO_TXT */

} /* void dump_macro_vars( struct lattice_struct *lattice, int time) */

#if 1
// void read_macro_vars( struct lattice_struct *lattice)
//##############################################################################
// 
// R E A D   M A C R O S C O P I C 
//
//  - Read the macro_vars variables from files.
//
void read_macro_vars( struct lattice_struct *lattice, int time)
{
  char   filename[1024];
  FILE   *in, *rho_in, *u_in;
  int    *node_ptr;
  int    n;
  double *macro_vars_ptr;
  int    frame;
  double max_u[2], ave_u[2];
  double max_rho, ave_rho;
  double rho_ratio, u_x_ratio, u_y_ratio;

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

  // R E A D   R H O   A N D   U 
  //
  //  - Read the density and velocity values at the active nodes to
  //    the rho and u dat files.
  //
  sprintf( filename, "rho%04d.dat", frame);
  if( !( rho_in = fopen( filename, "r+")))
  {
    printf("ERROR: fopen( \"%s\", \"r+\") = NULL.  Bye, bye!\n", filename);
    exit(1);
  }

  sprintf( filename, "u%04d.dat", frame);
  if( !( u_in = fopen( filename, "r+")))
  {
    printf("ERROR: fopen( \"%s\", \"r+\") = NULL.  Bye, bye!\n", filename);
    exit(1);
  }

  macro_vars_ptr = lattice->macro_vars[0].rho;

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

    fscanf( rho_in, "%lf\n", macro_vars_ptr++);

    fscanf( u_in,   "%lf ",  macro_vars_ptr++);

    fscanf( u_in,   "%lf\n", macro_vars_ptr++);

  }

  fclose(u_in);
#if VERBOSITY_LEVEL > 0
  sprintf( filename, "u%04d.dat", frame);
#endif /* VERBOSITY_LEVEL > 0 */
  printf("read_macro_vars() -- Read file \"%s\"\n", filename);
  fclose(rho_in);
#if VERBOSITY_LEVEL > 0
  sprintf( filename, "rho%04d.dat", frame);
#endif /* VERBOSITY_LEVEL > 0 */
  printf("read_macro_vars() -- Read file \"%s\"\n", filename);

} /* void read_macro_vars( struct lattice_struct *lattice, int time) */

#endif
// void dump_pdf( struct lattice_struct *lattice, int time)
//##############################################################################
//
// D U M P   P D F 
//
//  - Output the particle distribution functions to a text file.
//
//  - This is useful mainly for debugging with small problems.
//
void dump_pdf( struct lattice_struct *lattice, int time)
{
  char   filename[1024];
  FILE   *o_feq, 
         *o_f, 
         *o_ftemp;
  double *fptr, 
         *end_ptr;
  bc_ptr bc;
  int    frame;

  frame = time/lattice->param.FrameRate;

  sprintf( filename, "feq%04d.dat", frame);
  if( !( o_feq = fopen( filename, "w+")))
  {
    printf("ERROR: fopen( \"%s\", \"w+\") = NULL.  Bye, bye!\n", filename);
    exit(1);
  }

  sprintf( filename, "f%04d.dat", frame);
  if( !( o_f = fopen( filename, "w+")))
  {
    printf("ERROR: fopen( \"%s\", \"w+\") = NULL.  Bye, bye!\n", filename);
    exit(1);
  }

  sprintf( filename, "ftemp%04d.dat", frame);
  if( !( o_ftemp = fopen( filename, "w+")))
  {
    printf("ERROR: fopen( \"%s\", \"w+\") = NULL.  Bye, bye!\n", filename);
    exit(1);
  }

  bc      = lattice->bc;
  fptr    = lattice->pdf[0].feq;
  end_ptr = &(lattice->pdf[ lattice->NumNodes-1].ftemp[8]) + 1;
  while( fptr!=end_ptr)
  {
    if( 1 || !( bc++->bc_type & BC_SOLID_NODE))
    {
      fprintf( o_feq  , "%10.7f ", *fptr++);
      fprintf( o_feq  , "%10.7f ", *fptr++);
      fprintf( o_feq  , "%10.7f ", *fptr++);
      fprintf( o_feq  , "%10.7f ", *fptr++);
      fprintf( o_feq  , "%10.7f ", *fptr++);
      fprintf( o_feq  , "%10.7f ", *fptr++);
      fprintf( o_feq  , "%10.7f ", *fptr++);
      fprintf( o_feq  , "%10.7f ", *fptr++);
      fprintf( o_feq  , "%10.7f ", *fptr++);
    }
    else
    {
      fprintf( o_feq  , "%10.7f ", 0.);
      fprintf( o_feq  , "%10.7f ", 0.);
      fprintf( o_feq  , "%10.7f ", 0.);
      fprintf( o_feq  , "%10.7f ", 0.);
      fprintf( o_feq  , "%10.7f ", 0.);
      fprintf( o_feq  , "%10.7f ", 0.);
      fprintf( o_feq  , "%10.7f ", 0.);
      fprintf( o_feq  , "%10.7f ", 0.);
      fprintf( o_feq  , "%10.7f ", 0.);
      fptr+=9;
    }
    fprintf( o_feq  , "\n");

    fprintf( o_f    , "%10.7f ", *fptr++);
    fprintf( o_f    , "%10.7f ", *fptr++);
    fprintf( o_f    , "%10.7f ", *fptr++);
    fprintf( o_f    , "%10.7f ", *fptr++);
    fprintf( o_f    , "%10.7f ", *fptr++);
    fprintf( o_f    , "%10.7f ", *fptr++);
    fprintf( o_f    , "%10.7f ", *fptr++);
    fprintf( o_f    , "%10.7f ", *fptr++);
    fprintf( o_f    , "%10.7f ", *fptr++);
    fprintf( o_f    , "\n");

    fprintf( o_ftemp, "%10.7f ", *fptr++);
    fprintf( o_ftemp, "%10.7f ", *fptr++);
    fprintf( o_ftemp, "%10.7f ", *fptr++);
    fprintf( o_ftemp, "%10.7f ", *fptr++);
    fprintf( o_ftemp, "%10.7f ", *fptr++);
    fprintf( o_ftemp, "%10.7f ", *fptr++);
    fprintf( o_ftemp, "%10.7f ", *fptr++);
    fprintf( o_ftemp, "%10.7f ", *fptr++);
    fprintf( o_ftemp, "%10.7f ", *fptr++);
    fprintf( o_ftemp, "\n");

  } /* while( fptr!=end_ptr) */

  fclose( o_feq);
  fclose( o_f);
  fclose( o_ftemp);

} /* void dump_pdf( struct lattice_struct *lattice, int time) */

// void dump_lattice_info( struct lattice_struct *lattice)
//##############################################################################
//
// D U M P   L A T T I C E   I N F O 
//
//  - Output information about the lattice.
//
void dump_lattice_info( struct lattice_struct *lattice)
{
  FILE *o;

#if VERBOSITY_LEVEL >= 1
  printf("----------------------------------------");
  printf("----------------------------------------\n");

  printf("   lattice              = %x ", lattice);
  printf("(size = %3d bytes)\n", sizeof( struct lattice_struct));
  printf("&( lattice->NumNodes)   = %x ",  &( (lattice)->NumNodes));
  printf("(size = %3d bytes)\n",sizeof(int));
  printf("&( lattice->param)      = %x ",  &( (lattice)->param));
  printf("(size = %3d bytes)\n",sizeof(struct param_struct));
  printf("&( lattice->node)       = %x ",  &( (lattice)->node));
  printf("(size = %3d B/node * ",sizeof(struct node_struct));
  printf("%d nodes = %d bytes)\n", 
      (lattice)->NumNodes,
      (lattice)->NumNodes*sizeof(struct node_struct));
  printf("&( lattice->pdf)        = %x ",  &( (lattice)->pdf));
  printf("(size = %3d B/node * ", sizeof(struct pdf_struct));
  printf("%d nodes = %d bytes)\n", 
      (lattice)->NumNodes,
      (lattice)->NumNodes*sizeof(struct pdf_struct));
  printf("&( lattice->macro_vars) = %x ",  &( (lattice)->macro_vars));
  printf("(size = %3d B/node * ", sizeof(struct macro_vars_struct));
  printf("%d nodes = %d bytes)\n", 
      (lattice)->NumNodes, 
      (lattice)->NumNodes*sizeof(struct macro_vars_struct));
  printf("&( lattice->bc)         = %x ",  &( (lattice)->bc));
  printf("(size = %3d B/node * ", sizeof(struct bc_struct));
  printf("%d nodes = %d bytes)\n", 
      (lattice)->NumNodes, 
      (lattice)->NumNodes*sizeof(struct bc_struct));

  printf("----------------------------------------");
  printf("----------------------------------------\n");

  printf("Size of \"sparse\" lattice:  %d bytes = %f MB\n",

     sizeof(int)
   + sizeof(struct param_struct)
   + (lattice)->NumNodes*sizeof(struct node_struct)
   + (lattice)->NumNodes*sizeof(struct pdf_struct)
   + (lattice)->NumNodes*sizeof(struct macro_vars_struct)
   + (lattice)->NumNodes*sizeof(struct bc_struct),

   ( (double)sizeof(int)
   + (double)sizeof(struct param_struct)
   + (double)(lattice)->NumNodes*sizeof(struct node_struct)
   + (double)(lattice)->NumNodes*sizeof(struct pdf_struct)
   + (double)(lattice)->NumNodes*sizeof(struct macro_vars_struct)