latman.c

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

      {
        ip = ( i<ei) ? ( i+1) : ( 0 );
        in = ( i>0 ) ? ( i-1) : ( ei);
         
        if( matrix[j ][i ] != 1)
        {

          if( matrix[j ][ip] != 1) 
          { lattice->node[n].nn[0] = ((node_ptr)matrix[j ][ip])->n;}
          else
          { lattice->node[n].nn[0] = -1;}

          if( matrix[jp][i ] != 1) 
          { lattice->node[n].nn[1] = ((node_ptr)matrix[jp][i ])->n;}
          else
          { lattice->node[n].nn[1] = -1;}

          if( matrix[j ][in] != 1) 
          { lattice->node[n].nn[2] = ((node_ptr)matrix[j ][in])->n;}
          else
          { lattice->node[n].nn[2] = -1;}

          if( matrix[jn][i ] != 1) 
          { lattice->node[n].nn[3] = ((node_ptr)matrix[jn][i ])->n;}
          else
          { lattice->node[n].nn[3] = -1;}

          if( matrix[jp][ip] != 1) 
          { lattice->node[n].nn[4] = ((node_ptr)matrix[jp][ip])->n;}
          else
          { lattice->node[n].nn[4] = -1;}

          if( matrix[jp][in] != 1) 
          { lattice->node[n].nn[5] = ((node_ptr)matrix[jp][in])->n;}
          else
          { lattice->node[n].nn[5] = -1;}

          if( matrix[jn][in] != 1) 
          { lattice->node[n].nn[6] = ((node_ptr)matrix[jn][in])->n;}
          else
          { lattice->node[n].nn[6] = -1;}

          if( matrix[jn][ip] != 1) 
          { lattice->node[n].nn[7] = ((node_ptr)matrix[jn][ip])->n;}
          else
          { lattice->node[n].nn[7] = -1;}

          n++;

        } /* if( matrix[j][i] != 1) */
      } /* 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 ] != 1)
      {

        if( matrix[j ][ip] != 1) 
        { lattice->node[n].nn[0] = ((node_ptr)matrix[j ][ip])->n;}
        else
        { lattice->node[n].nn[0] = -1;}

        //if( matrix[jp][i ] != 1) 
        //{ lattice->node[n].nn[1] = ((node_ptr)matrix[jp][i ])->n;}
        //else
        //{ lattice->node[n].nn[1] = -1;}
        lattice->node[n].nn[1] = -1;

        if( matrix[j ][in] != 1) 
        { lattice->node[n].nn[2] = ((node_ptr)matrix[j ][in])->n;}
        else
        { lattice->node[n].nn[2] = -1;}

        if( matrix[jn][i ] != 1) 
        { lattice->node[n].nn[3] = ((node_ptr)matrix[jn][i ])->n;}
        else
        { lattice->node[n].nn[3] = -1;}

        //if( matrix[jp][ip] != 1) 
        //{ lattice->node[n].nn[4] = ((node_ptr)matrix[jp][ip])->n;}
        //else
        //{ lattice->node[n].nn[4] = -1;}
        lattice->node[n].nn[4] = -1;

        //if( matrix[jp][in] != 1) 
        //{ lattice->node[n].nn[5] = ((node_ptr)matrix[jp][in])->n;}
        //else
        //{ lattice->node[n].nn[5] = -1;}
        lattice->node[n].nn[5] = -1;

        if( matrix[jn][in] != 1) 
        { lattice->node[n].nn[6] = ((node_ptr)matrix[jn][in])->n;}
        else
        { lattice->node[n].nn[6] = -1;}

        if( matrix[jn][ip] != 1) 
        { lattice->node[n].nn[7] = ((node_ptr)matrix[jn][ip])->n;}
        else
        { lattice->node[n].nn[7] = -1;}

        n++;

      } /* if( matrix[j][i] != 1) */
    } /* for( i=0; i<=ei; i++) */

  } /* if( lattice->periodic_y) else */

} /* 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( struct lattice_struct **lattice)
{
  // Variable declarations
  int **matrix;
  int i, j;
  int n;
  int width, height;
  char filename[1024];

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

  assert(*lattice!=NULL);

  // Read problem parameters
  read_params( *lattice);

  // Default periodicity. This will be adjusted in read_bcs() 
  // depending on flow boundaries.
  (*lattice)->periodic_x = 1;
  (*lattice)->periodic_y = 1;

  // 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));
  }

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

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

  // Get boundary conditions.
  sprintf( filename, "%dx%d.bmp", (*lattice)->param.LX, (*lattice)->param.LY);
  read_bcs( filename, *lattice, matrix);

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

  // Allocate NumNodes particle distribution functions.
  (*lattice)->pdf = 
    ( struct pdf_struct*)malloc( 
        (*lattice)->NumNodes*sizeof( struct pdf_struct));
  assert( (*lattice)->pdf!=NULL);

  // Allocate NumNodes macroscopic variables.
  (*lattice)->macro_vars = 
    ( struct macro_vars_struct*)malloc( 
        (*lattice)->NumNodes*sizeof( struct macro_vars_struct));
  assert( (*lattice)->macro_vars!=NULL);

  // Write info to files
  dump_lattice_info( *lattice);

  dump_node_info( *lattice);

  dump_params( *lattice);

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

// void init_problem( struct lattice_struct *lattice)
//##############################################################################
//
// I N I T   P R O B L E M 
//
//  - Initialize the problem on the lattice.
//
//    - Set the initial density and velocity.
//
//    - Compute the initial feq.
//
void init_problem( struct lattice_struct *lattice)
{
  int n;
  double *macro_var_ptr;
  double *f;
  bc_ptr bc;

#if VERBOSITY_LEVEL > 0
  printf("init_problem() -- Initilizing problem...\n");
#endif /* VERBOSITY_LEVEL > 0 */

  macro_var_ptr = lattice->macro_vars[0].rho;
  bc            = lattice->bc;

  for( n=0; n<lattice->NumNodes; n++)
  {
    // Set initial density.
    if( 1 || !( bc++->bc_type & BC_SOLID_NODE))
    {
      *macro_var_ptr++ = 1.;
    }
    else
    {
      *macro_var_ptr++ = 0.;
    }

    // Set initial velocty.
    *macro_var_ptr++ = 0.;
    *macro_var_ptr++ = 0.;

  } /* for( n=0; n<lattice->NumNodes; n++) */

  // Compute initial feq.
  compute_feq( lattice);

  f = lattice->pdf[0].f;

  for( n=0; n<lattice->NumNodes; n++)
  {
    if( 1 || !( bc[n].bc_type & BC_SOLID_NODE))
    {
#if 1
      *f++ = *(f-9);
      *f++ = *(f-9);
      *f++ = *(f-9);
      *f++ = *(f-9);
      *f++ = *(f-9);
      *f++ = *(f-9);
      *f++ = *(f-9);
      *f++ = *(f-9);
      *f++ = *(f-9);

      *f++ = *(f-9);
      *f++ = *(f-9);
      *f++ = *(f-9);
      *f++ = *(f-9);
      *f++ = *(f-9);
      *f++ = *(f-9);
      *f++ = *(f-9);
      *f++ = *(f-9);
      *f++ = *(f-9);

      f+=9;
#else
      // Debug info.  To track during the streaming step.
      *f++ = (double)n + 0./10.;// *(f-9);
      *f++ = (double)n + 1./10.;// *(f-9);
      *f++ = (double)n + 2./10.;// *(f-9);
      *f++ = (double)n + 3./10.;// *(f-9);
      *f++ = (double)n + 4./10.;// *(f-9);
      *f++ = (double)n + 5./10.;// *(f-9);
      *f++ = (double)n + 6./10.;// *(f-9);
      *f++ = (double)n + 7./10.;// *(f-9);
      *f++ = (double)n + 8./10.;// *(f-9);

      *f++ = (double)n;// *(f-9);
      *f++ = (double)n;// *(f-9);
      *f++ = (double)n;// *(f-9);
      *f++ = (double)n;// *(f-9);
      *f++ = (double)n;// *(f-9);
      *f++ = (double)n;// *(f-9);
      *f++ = (double)n;// *(f-9);
      *f++ = (double)n;// *(f-9);
      *f++ = (double)n;// *(f-9);

      f+=9;
#endif
    }
    else
    {
      *f++ = 0.;
      *f++ = 0.;
      *f++ = 0.;
      *f++ = 0.;
      *f++ = 0.;
      *f++ = 0.;
      *f++ = 0.;
      *f++ = 0.;
      *f++ = 0.;
          
      *f++ = 0.;
      *f++ = 0.;
      *f++ = 0.;
      *f++ = 0.;
      *f++ = 0.;
      *f++ = 0.;
      *f++ = 0.;
      *f++ = 0.;
      *f++ = 0.;
   
      f+=9;
    }

  } /* for( n=0; n<lattice->NumNodes; n++) */

  //dump_pdf( lattice, /*time=*/ 0);

  compute_macro_vars( lattice);

  dump_macro_vars( lattice, /*time=*/ 0);

#if VERBOSITY_LEVEL > 0
  printf("init_problem() -- Problem initialized.\n");
#endif /* VERBOSITY_LEVEL > 0 */

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

// void destruct_lattice( struct lattice_struct *lattice)
//##############################################################################
//
// D E S T R U C T   L A T T I C E 
//
//  - Destruct lattice.
//
void destruct_lattice( struct lattice_struct *lattice)
{
  assert( lattice!=NULL);

  assert( lattice->node!=NULL);
  assert( lattice->pdf!=NULL);
  assert( lattice->macro_vars!=NULL);
  assert( lattice->bc!=NULL);

  free( lattice->node);
  free( lattice->pdf);
  free( lattice->macro_vars);
  free( lattice->bc);

  free( lattice);

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