latman.c

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

//##############################################################################
//
// Copyright (C), 2005, Michael Sukop and Danny Thorne
//
// all_nodes_latman.c
//
//  - Lattice Manager.
//
//  - Routines for managing a lattice:
//
//    - construct
//    - init
//    - destruct
//
//  - This file, with prefix "all_nodes_", is for the version of the code
//    that stores all nodes of the domain even if they are interior solid
//    nodes that are not involved in any computations.  This is more
//    efficient, in spite of storing unused nodes, if the ratio of
//    interior solid nodes to total nodes is sufficiently low.  How
//    low is sufficient? is a difficult question.  If storage is the
//    only consideration, then the two approaches balance at somewhere
//    around .75 .  But more observations need to be made to characterize
//    the trade-off in terms of computational efficiency.
//

// void process_matrix( struct lattice_struct *lattice, int **matrix)
//##############################################################################
//
// P R O C E S S   M A T R I X 
//
//  - Process the soil matrix.
//
//  - Convert the full matrix representation into sparse lattice form.
//
//  - This routine is useful when the domain is read from a BMP file
//    representing the full lattice.  If more general mechanisms
//    are developed (e.g., reading parameterized fracture information),
//    then determining the sparse lattice data structure will require a
//    corresponding mechanism (and storing the full lattice, even for
//    pre-/post-processing, may become undesirable if not impossible).
//
void process_matrix( struct lattice_struct *lattice, int **matrix, int subs)
{
  // Variable declarations.
  int i,  j;
  int in, jn;
  int ip, jp;
  int ei, ej;
  int n;
  int NumActiveNodes;

  // Ending indices.
  ei = lattice->param.LX-1;
  ej = lattice->param.LY-1;

  // Mark solid nodes that have fluid nodes as neighbors and
  // count the total number of nodes requiring storage.
  NumActiveNodes = 0;
  if( lattice->periodic_y[subs])
  {
    for( j=0; j<=ej; j++)
    {
      jp = ( j<ej) ? ( j+1) : ( 0 );
      jn = ( j>0 ) ? ( j-1) : ( ej);

      for( i=0; i<=ei; i++)
      {
        ip = ( i<ei) ? ( i+1) : ( 0 );
        in = ( i>0 ) ? ( i-1) : ( ei);

        if( matrix[j ][i ] == 0) 
        { 
          NumActiveNodes++;

          if( matrix[j ][ip] == 1) { matrix[j ][ip] = 2; NumActiveNodes++;}
          if( matrix[jp][i ] == 1) { matrix[jp][i ] = 2; NumActiveNodes++;}
          if( matrix[j ][in] == 1) { matrix[j ][in] = 2; NumActiveNodes++;}
          if( matrix[jn][i ] == 1) { matrix[jn][i ] = 2; NumActiveNodes++;}
          if( matrix[jp][ip] == 1) { matrix[jp][ip] = 2; NumActiveNodes++;}
          if( matrix[jp][in] == 1) { matrix[jp][in] = 2; NumActiveNodes++;}
          if( matrix[jn][in] == 1) { matrix[jn][in] = 2; NumActiveNodes++;}
          if( matrix[jn][ip] == 1) { matrix[jn][ip] = 2; NumActiveNodes++;}

        }

      } /* for( i=0; i<=ei; i++) */
    } /* for( j=0; j<=ej; j++) */
  } /* if( lattice->periodic_y[subs]) */

  else /* !lattice->periodic_y[subs] */
  {
    j  = 0;
    jp = 1;
    for( i=0; i<=ei; i++)
    {
      ip = ( i<ei) ? ( i+1) : ( 0 );
      in = ( i>0 ) ? ( i-1) : ( ei);

      if( matrix[j ][i ] == 0) 
      { 
        NumActiveNodes++;

        if( matrix[j ][ip] == 1) { matrix[j ][ip] = 2; NumActiveNodes++;}
        if( matrix[jp][i ] == 1) { matrix[jp][i ] = 2; NumActiveNodes++;}
        if( matrix[j ][in] == 1) { matrix[j ][in] = 2; NumActiveNodes++;}
        //if( matrix[jn][i ] == 1) { matrix[jn][i ] = 2; NumActiveNodes++;}
        if( matrix[jp][ip] == 1) { matrix[jp][ip] = 2; NumActiveNodes++;}
        if( matrix[jp][in] == 1) { matrix[jp][in] = 2; NumActiveNodes++;}
        //if( matrix[jn][in] == 1) { matrix[jn][in] = 2; NumActiveNodes++;}
        //if( matrix[jn][ip] == 1) { matrix[jn][ip] = 2; NumActiveNodes++;}

      }

    } /* for( i=0; i<=ei; i++) */


    for( j=1; j<ej; j++)
    {
      jp = j+1;
      jn = j-1;

      for( i=0; i<=ei; i++)
      {
        ip = ( i<ei) ? ( i+1) : ( 0 );
        in = ( i>0 ) ? ( i-1) : ( ei);

        if( matrix[j ][i ] == 0) 
        { 
          NumActiveNodes++;

          if( matrix[j ][ip] == 1) { matrix[j ][ip] = 2; NumActiveNodes++;}
          if( matrix[jp][i ] == 1) { matrix[jp][i ] = 2; NumActiveNodes++;}
          if( matrix[j ][in] == 1) { matrix[j ][in] = 2; NumActiveNodes++;}
          if( matrix[jn][i ] == 1) { matrix[jn][i ] = 2; NumActiveNodes++;}
          if( matrix[jp][ip] == 1) { matrix[jp][ip] = 2; NumActiveNodes++;}
          if( matrix[jp][in] == 1) { matrix[jp][in] = 2; NumActiveNodes++;}
          if( matrix[jn][in] == 1) { matrix[jn][in] = 2; NumActiveNodes++;}
          if( matrix[jn][ip] == 1) { matrix[jn][ip] = 2; NumActiveNodes++;}

        }

      } /* for( i=0; i<=ei; i++) */
    } /* for( j=0; j<=ej; j++) */

    j  = ej;
    jn = ej-1;
    for( i=0; i<=ei; i++)
    {
      ip = ( i<ei) ? ( i+1) : ( 0 );
      in = ( i>0 ) ? ( i-1) : ( ei);

      if( matrix[j ][i ] == 0) 
      { 
        NumActiveNodes++;

        if( matrix[j ][ip] == 1) { matrix[j ][ip] = 2; NumActiveNodes++;}
        //if( matrix[jp][i ] == 1) { matrix[jp][i ] = 2; NumActiveNodes++;}
        if( matrix[j ][in] == 1) { matrix[j ][in] = 2; NumActiveNodes++;}
        if( matrix[jn][i ] == 1) { matrix[jn][i ] = 2; NumActiveNodes++;}
        //if( matrix[jp][ip] == 1) { matrix[jp][ip] = 2; NumActiveNodes++;}
        //if( matrix[jp][in] == 1) { matrix[jp][in] = 2; NumActiveNodes++;}
        if( matrix[jn][in] == 1) { matrix[jn][in] = 2; NumActiveNodes++;}
        if( matrix[jn][ip] == 1) { matrix[jn][ip] = 2; NumActiveNodes++;}

      }

    } /* for( i=0; i<=ei; i++) */

  } /* if( lattice->periodic_y[subs]) else */

#if VERBOSITY_LEVEL > 0
  printf( "[%s,%d] process_matrix() -- NumActiveNodes = %d\n", 
    __FILE__, __LINE__, NumActiveNodes);
#endif /* VERBOSITY_LEVEL > 0 */


#if 0 // Dump the matrix contents to the screen.
  for( j=0; j<=ej; j++)
  {
    for( i=0; i<=ei; i++)
    {
      printf(" %d", matrix[j][i]);
    }
    printf("\n");
  }
  //exit(1);
#endif


  // Set lattice->NumNodes in the lattice.
  lattice->NumNodes = lattice->param.LX * lattice->param.LY;

#if VERBOSITY_LEVEL > 0
  printf("[%s,%d] process_matrix() -- NumNodes = %d\n", 
    __FILE__,__LINE__, lattice->NumNodes);
#endif /* VERBOSITY_LEVEL > 0 */

  // Allocate memory for lattice->NumNodes boundary conditions.
  lattice->bc[subs]= 
    ( struct bc_struct*)malloc( lattice->NumNodes*sizeof( struct bc_struct));
  assert( lattice->bc[subs]!=NULL);

  // Set coordinates and index of lattice nodes.
  // Use matrix entries to store pointers to associated nodes to
  // facilitate finding neighbors on a following traversal.
  for( j=0; j<=ej; j++)
  {
    n = j*lattice->param.LX;

    for( i=0; i<=ei; i++, n++)
    {
      switch( matrix[j][i])
      {
        case 0: lattice->bc[subs][n].bc_type = 0;             break;
        case 1: lattice->bc[subs][n].bc_type = INACTIVE_NODE; break;
        case 2: lattice->bc[subs][n].bc_type = BC_SOLID_NODE; break;
        default:
          printf("all_nodes_latman.c: process_matrix() -- "
              "Unhandled case matrix[%d][%d]=%d. Exiting!\n",
              j, i, matrix[j][i] );
          exit(1);
          break;
      }
    } /* for( i=0; i<=ei; i++) */
  } /* for( j=0; j<=ej; j++) */


#if 0 // Dump the matrix contents to the screen.
  for( j=0; j<=ej; j++)
  {
    for( i=0; i<=ei; i++)
    {
      printf(" %d", matrix[j][i]);
    }
    printf("\n");
  }
#endif

#if 0 // Dump BCs to screen.
  for( n=0; n<lattice->NumNodes; n++)
  {
    printf("%d (%d,%d), %d\n", 
      n, 
      n%lattice->param.LX, 
      n/lattice->param.LX, 
      lattice->bc[subs][n].bc_type);
  }
  printf("\n");
  for( n=0, j=0; j<=ej; j++)
  {
    for( i=0; i<=ei; i++, n++)
    {
      printf(" %d", lattice->bc[subs][n].bc_type);
    }
    printf("\n");
  }
#endif

} /* void process_matrix( struct lattice_struct *lattice, int **matrix) */

// void construct_lattice( struct lattice_struct *lattice)
//##############################################################################
//
// C O N S T R U C T   L A T T I C E 
//
//  - Construct lattice.
//
void construct_lattice( lattice_ptr *lattice, int argc, char **argv)
{
  // Variable declarations
  int    **matrix;
  int    i, 
         j;
  int    n;
  int    subs;
  int    width, 
         height;
  char   filename[1024];
#if POROUS_MEDIA
  FILE   *in;
  char   r, g, b;
  struct bitmap_info_header bmih;
#endif /* POROUS_MEDIA */

  // Allocate the lattice structure.
  *lattice = ( struct lattice_struct*)malloc( sizeof(struct lattice_struct));

  assert(*lattice!=NULL);

  if( argc == 2)
  {
    printf("argv = \"%s\"\n", argv[1]);
    strcpy( filename, argv[1]);
    printf("filename = \"%s\"\n", filename);
  }
  else if( argc == 1)
  {
    sprintf(filename, "./in/%s", "params.in");
    printf("filename = \"%s\"\n", filename);
  }
  else
  {
    printf("\n\nusage: ./lb2d [infile]\n\n\n");
    exit(1);
  }

  // Read problem parameters
  read_params( *lattice, filename);

  // Read soil matrix from bmp file.
  matrix = (int**)malloc( (*lattice)->param.LY*sizeof(int*));
  for( j=0; j<(*lattice)->param.LY; j++)
  {
    matrix[j] = (int*)malloc( (*lattice)->param.LX*sizeof(int));
  }

  // Initialize matrix[][].
  for( j=0; j<(*lattice)->param.LY; j++)
  {
    for( i=0; i<(*lattice)->param.LX; i++)
    {
      matrix[j][i] = 0;
    }
  }

 for( subs=0; subs<NUM_FLUID_COMPONENTS; subs++)
 {
  // Default periodicity. This will be adjusted in read_bcs() 
  // or process_bcs() depending on flow boundaries.
  (*lattice)->periodic_x[subs] = 1;
  (*lattice)->periodic_y[subs] = 1;

  // Get solids.
  sprintf( filename, "./in/%dx%d.bmp", (*lattice)->param.LX, (*lattice)->param.LY);
  spy_bmp( filename, *lattice, matrix);

  // Determine active nodes.
  process_matrix( *lattice, matrix, subs);
  assert( (*lattice)->bc[subs]!=NULL);

#if 0
  // Read boundary conditions from BMP files.
  // Eventually this mechanism will be very general to handle 
  // (somewhat?) arbitrary arrangements of boundary conditions.
  read_bcs( *lattice, matrix);
#else
  // Process boundary conditions based on flags read from params.in .
  // This will only support standard inflow and outflow boundaries along
  // entire sides of the lattice.
  process_bcs( *lattice, subs);
#endif

 } /* for( subs=0; subs<NUM_FLUID_COMPONENTS; subs++) */

  // Deallocate memory used for storing the full matrix.
#if VERBOSITY_LEVEL > 0
  printf("latman.c: contruct_lattice() -- Free the matrix.\n");
#endif /* VERBOSITY_LEVEL > 0 */
  for( n=0; n<(*lattice)->param.LY; n++)
  {
    free( matrix[n]);
  }
  free( matrix);
#if VERBOSITY_LEVEL > 0
  printf("latman.c: contruct_lattice() -- Matrix is free.\n");
#endif /* VERBOSITY_LEVEL > 0 */

 for( subs=0; subs<NUM_FLUID_COMPONENTS; subs++)
 {
  // Allocate NumNodes particle distribution functions.
  (*lattice)->pdf[subs] = 
    ( struct pdf_struct*)malloc( 
        (*lattice)->NumNodes*sizeof( struct pdf_struct));
  if( (*lattice)->pdf[subs] == NULL)
  {
    printf(
      "construct_lattice() -- ERROR:  "
      "Attempt to allocate %d struct pdf_struct types failed.  "
      "Exiting!\n",
      (*lattice)->NumNodes
      );
    exit(1);
  }

  // Allocate NumNodes macroscopic variables.
  (*lattice)->macro_vars[subs] = 
    ( struct macro_vars_struct*)malloc( 
        (*lattice)->NumNodes*sizeof( struct macro_vars_struct));
  if( (*lattice)->macro_vars[subs]==NULL)
  {
    printf(
      "construct_lattice() -- ERROR:  "
      "Attempt to allocate %d struct macro_vars_struct types failed.  "
      "Exiting!\n",
      (*lattice)->NumNodes
      );
    exit(1);
  }

#if NON_LOCAL_FORCES
  // Allocate NumNodes elements for force.
  (*lattice)->force[subs] = 
    ( struct force_struct*)malloc( 
        (*lattice)->NumNodes*sizeof( struct force_struct));
  if( (*lattice)->force[subs]==NULL)
  {
    printf(
      "construct_lattice() -- ERROR:  "
      "Attempt to allocate %d struct force_struct types failed.  "
      "Exiting!\n",
      (*lattice)->NumNodes
      );
    exit(1);
  }
#endif /* NON_LOCAL_FORCES */
 } /* for( subs=0; subs<NUM_FLUID_COMPONENTS; subs++) */

#if STORE_UEQ
  // Allocate NumNodes elements for ueq.
  (*lattice)->ueq = 
    ( struct ueq_struct*)malloc( 
        (*lattice)->NumNodes*sizeof( struct ueq_struct));
  if( (*lattice)->ueq==NULL)
  {
    printf(
      "construct_lattice() -- ERROR:  "
      "Attempt to allocate %d struct ueq_struct types failed.  "
      "Exiting!\n",
      (*lattice)->NumNodes
      );
    exit(1);
  }
#endif /* STORE_UEQ */

#if POROUS_MEDIA
  switch( (*lattice)->param.ns_flag)
  {
    case 0:
    {
      if( (*lattice)->param.ns > 1.)
      {
        printf(
          "latman.c: construct_lattice() -- "
          "ERROR: ns = %f. "
          "Should have 0 <= ns <=1. "
          "Exiting!\n", (*lattice)->param.ns
          );
        exit(1);
      }
      break;
    }
    case 1:
    {
    // Allocate space for ns values.
    (*lattice)->ns = 
      (struct ns_struct*)malloc( (*lattice)->NumNodes*sizeof(struct ns_struct));
    if( (*lattice)->ns==NULL)
    {
      printf(
        "construct_lattice() -- ERROR:  "
        "Attempt to allocate %d struct ns_struct types failed.  "
        "Exiting!\n",
        (*lattice)->NumNodes
        );
      exit(1);
    }

    // Try to read ns<LX>x<LY>.bmp file.
    sprintf( filename, "./in/ns%dx%d.bmp", 
           (*lattice)->param.LX, 
           (*lattice)->param.LY);
    if( in = fopen( filename, "r+"))
    {
      printf("%s %d >> Reading file \"%s\".\n",__FILE__,__LINE__,filename);
      bmp_read_header( in, &bmih);
      for( n=0; n<(*lattice)->NumNodes; n++)
      {
        bmp_read_entry( in, bmih, &r, &g, &b);

        // Verify grayscale.
        if(    (double)r != (double)g 
            || (double)g != (double)b 
            || (double)r != (double)b)
        {
          printf(
            "%s %d >> latman.c: construct_lattice() -- "
            "n=%d:  [ r g b] = [ %3u %3u %3u]\n",__FILE__,__LINE__,
            n, (unsigned int)r%256, (unsigned int)g%256, (unsigned int)b%256);
          printf(
            "%s %d >> latman.c: construct_lattice() -- "
            "ERROR: File %s needs to be grayscale. "
            "Exiting!\n",__FILE__,__LINE__, filename);
          exit(1);
        }

        // Assign ns value.
        (*lattice)->ns[n].ns = ((double)((unsigned int)r%256))/255.;
#if 0 && VERBOSITY_LEVEL>0
        printf("%s %d >> n=%d, ns=%f\n", 
            __FILE__, __LINE__, n, (*lattice)->ns[n].ns);
#endif /* 1 && VERBOSITY_LEVEL>0 */