bcs.c

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

    while( ftemp < ftemp_mid)
    {
      // North, Inflow
      //u = u_in[((double)rand()/(double)RAND_MAX<.5)?(0):(1)][subs];
      u = u_in[1][subs];
      rho = ( ftemp[0] + ftemp[1] + ftemp[3]
             + 2.*( ftemp[2] + ftemp[5] + ftemp[6])) 
           / ( 1. + u);
    
      ftemp[4] = ftemp[2] - (2./3.)*rho*u;
    
      ftemp[7] = ftemp[5] + (1./2.)*( ftemp[1] - ftemp[3])
                          - (1./6.)*rho*u;
    
      ftemp[8] = ftemp[6] + (1./2.)*( ftemp[3] - ftemp[1])
                          - (1./6.)*rho*u;

      ftemp += 27;//( sizeof(struct pdf_struct)/8);

    } /* while( ftemp < ftemp_end) */
    while( ftemp < ftemp_end)
    {
      // North, Inflow
      //u = u_in[((double)rand()/(double)RAND_MAX<.5)?(0):(1)][subs];
      u = 0.;
      rho = ( ftemp[0] + ftemp[1] + ftemp[3]
             + 2.*( ftemp[2] + ftemp[5] + ftemp[6])) 
           / ( 1. + u);
    
      ftemp[4] = ftemp[2] - (2./3.)*rho*u;
    
      ftemp[7] = ftemp[5] + (1./2.)*( ftemp[1] - ftemp[3])
                          - (1./6.)*rho*u;
    
      ftemp[8] = ftemp[6] + (1./2.)*( ftemp[3] - ftemp[1])
                          - (1./6.)*rho*u;

      ftemp += 27;//( sizeof(struct pdf_struct)/8);

    } /* while( ftemp < ftemp_end) */
  }
  if( lattice->param.velocity_s_in[subs] )
  {
  }
  if( lattice->param.velocity_n_out[subs])
  {
  }
  if( lattice->param.velocity_s_out[subs])
  {
//printf("bcs() -- velocity_s_out[%d]\n", subs);
    ftemp = lattice->pdf[subs][0].ftemp;
    ftemp_mid = lattice->pdf[subs][lattice->param.LX/8].ftemp;
    ftemp_end = lattice->pdf[subs][lattice->param.LX].ftemp;
    //u = ((subs==1)?(-1):(1))*lattice->param.uy_out;
    while( ftemp < ftemp_mid)
    {
      //u = u_out[((double)rand()/(double)RAND_MAX<.5)?(0):(1)][subs];
      u = 0.;
      rho = ( ftemp[0] + ftemp[1] + ftemp[3]
             + 2.*( ftemp[4] + ftemp[7] + ftemp[8])) 
           / ( 1. - u);
    
      ftemp[2] = ftemp[4] + (2./3.)*rho*u;
    
      ftemp[5] = ftemp[7] + (1./2.)*( ftemp[3] - ftemp[1])
                          + (1./6.)*rho*u;
    
      ftemp[6] = ftemp[8] + (1./2.)*( ftemp[1] - ftemp[3])
                          + (1./6.)*rho*u;

      ftemp += 27;//( sizeof(struct pdf_struct)/8);

    } /* while( ftemp < ftemp_end) */
    ftemp_mid = lattice->pdf[subs][lattice->param.LX/4].ftemp;
    while( ftemp < ftemp_mid)
    {
      //u = u_out[((double)rand()/(double)RAND_MAX<.5)?(0):(1)][subs];
      u = u_out[0][subs];
      rho = ( ftemp[0] + ftemp[1] + ftemp[3]
             + 2.*( ftemp[4] + ftemp[7] + ftemp[8])) 
           / ( 1. - u);
    
      ftemp[2] = ftemp[4] + (2./3.)*rho*u;
    
      ftemp[5] = ftemp[7] + (1./2.)*( ftemp[3] - ftemp[1])
                          + (1./6.)*rho*u;
    
      ftemp[6] = ftemp[8] + (1./2.)*( ftemp[1] - ftemp[3])
                          + (1./6.)*rho*u;

      ftemp += 27;//( sizeof(struct pdf_struct)/8);

    } /* while( ftemp < ftemp_end) */
    ftemp_mid = lattice->pdf[subs][lattice->param.LX/2].ftemp;
    while( ftemp < ftemp_mid)
    {
      //u = u_out[((double)rand()/(double)RAND_MAX<.5)?(0):(1)][subs];
      u = 0.;
      rho = ( ftemp[0] + ftemp[1] + ftemp[3]
             + 2.*( ftemp[4] + ftemp[7] + ftemp[8])) 
           / ( 1. - u);
    
      ftemp[2] = ftemp[4] + (2./3.)*rho*u;
    
      ftemp[5] = ftemp[7] + (1./2.)*( ftemp[3] - ftemp[1])
                          + (1./6.)*rho*u;
    
      ftemp[6] = ftemp[8] + (1./2.)*( ftemp[1] - ftemp[3])
                          + (1./6.)*rho*u;

      ftemp += 27;//( sizeof(struct pdf_struct)/8);

    } /* while( ftemp < ftemp_end) */
    ftemp_mid = lattice->pdf[subs][5*lattice->param.LX/8].ftemp;
    while( ftemp < ftemp_mid)
    {
      //u = u_out[((double)rand()/(double)RAND_MAX<.5)?(0):(1)][subs];
      u = u_out[1][subs];
      rho = ( ftemp[0] + ftemp[1] + ftemp[3]
             + 2.*( ftemp[4] + ftemp[7] + ftemp[8])) 
           / ( 1. - u);
    
      ftemp[2] = ftemp[4] + (2./3.)*rho*u;
    
      ftemp[5] = ftemp[7] + (1./2.)*( ftemp[3] - ftemp[1])
                          + (1./6.)*rho*u;
    
      ftemp[6] = ftemp[8] + (1./2.)*( ftemp[1] - ftemp[3])
                          + (1./6.)*rho*u;

      ftemp += 27;//( sizeof(struct pdf_struct)/8);

    } /* while( ftemp < ftemp_end) */
    while( ftemp < ftemp_end)
    {
      //u = u_out[((double)rand()/(double)RAND_MAX<.5)?(0):(1)][subs];
      u = 0.;
      rho = ( ftemp[0] + ftemp[1] + ftemp[3]
             + 2.*( ftemp[4] + ftemp[7] + ftemp[8])) 
           / ( 1. - u);
    
      ftemp[2] = ftemp[4] + (2./3.)*rho*u;
    
      ftemp[5] = ftemp[7] + (1./2.)*( ftemp[3] - ftemp[1])
                          + (1./6.)*rho*u;
    
      ftemp[6] = ftemp[8] + (1./2.)*( ftemp[1] - ftemp[3])
                          + (1./6.)*rho*u;

      ftemp += 27;//( sizeof(struct pdf_struct)/8);

    } /* while( ftemp < ftemp_end) */

  } /* if( lattice->param.velocity_s_out[subs]) */
#endif

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

} /* void bcs( lattice_ptr lattice) */
#else
// void bcs( lattice_ptr lattice)
//##############################################################################
//
// B C S
//
//  - Apply boundary conditions.
//
void bcs( lattice_ptr lattice)
{
  int    i, j, n;
  int    subs;
  int    *bc_type;
  double *ftemp, *ftemp_end, *ftemp_mid;
  double  v;
  double  c0;
  double  D;
  double  c2;

#if RHO0_TEST
//------------------------------------------------------------------[ TEST ]----
  double *rho0;
//------------------------------------------------------------------[ TEST ]----
#endif /* RHO0_TEST */
  double u_x,
         u_y;
  double u_in[2][2], 
         u_out[2][2], 
         u,
         rho;
  double c;

  // NOTE: Should previously (in initialization stage) have checked to 
  // insure no solid nodes on inflow/outflow boundaries. Do not do it here 
  // inside the loop!

 for( subs=0; subs<(NUM_FLUID_COMPONENTS)-(INAMURO_SIGMA_COMPONENT); subs++)
 {
  // P R E S S U R E   N O R T H   I N F L O W   B C 
  //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  // pressure north inflow
  //  -- Pressure boundary on north side using inflow pressure condition.
  if( lattice->param.pressure_n_in[subs] )
  {
//printf("bcs() %s %d >> pressure_n_in[%d]\n", __FILE__, __LINE__, subs);
    ftemp = lattice->pdf[subs][lattice->NumNodes-lattice->param.LX].ftemp;
    ftemp_end = lattice->pdf[subs][lattice->NumNodes].ftemp;
    while( ftemp < ftemp_end)
    {
      // North, Inflow
      if( lattice->param.incompressible)
      {
        u_y = -lattice->param.rho_in
            + ( ftemp[0] + ftemp[1] + ftemp[3]
              + 2.*( ftemp[2] + ftemp[5] + ftemp[6]));
        c = u_y;
      }
      else // compressible
      {
        u_y = -1.
            + ( ftemp[0] + ftemp[1] + ftemp[3]
              + 2.*( ftemp[2] + ftemp[5] + ftemp[6])) 
            / lattice->param.rho_in;
        c = u_y*lattice->param.rho_in;
      }
    
      ftemp[4] = ftemp[2] - (2./3.)*c;
    
      ftemp[7] = ftemp[5] + (1./2.)*( ftemp[1] - ftemp[3])
                          - (1./6.)*c;
    
      ftemp[8] = ftemp[6] + (1./2.)*( ftemp[3] - ftemp[1])
                          - (1./6.)*c; 

      ftemp += ( sizeof(struct pdf_struct)/8);

    } /* while( ftemp < ftemp_end) */

  } /* if( lattice->param.pressure_n_in[subs] ) */

  // P R E S S U R E   S O U T H   I N F L O W   B C 
  //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  // pressure south inflow
  //  -- Pressure boundary on south side using inflow pressure condition.
  if( lattice->param.pressure_s_in[subs] )
  {
//printf("bcs() %s %d >> pressure_s_in[%d]\n", __FILE__, __LINE__, subs);
    ftemp = lattice->pdf[subs][0].ftemp;
    ftemp_end = lattice->pdf[subs][lattice->param.LX].ftemp;
    while( ftemp < ftemp_end)
    {
      // South, Inflow
      if( lattice->param.incompressible)
      {
        u_y = lattice->param.rho_in
            - ( ftemp[0] + ftemp[1] + ftemp[3]
              + 2.*( ftemp[4] + ftemp[7] + ftemp[8]));
        c = u_y;
      }
      else // compressible
      {
        u_y = 1.
            - ( ftemp[0] + ftemp[1] + ftemp[3]
              + 2.*( ftemp[4] + ftemp[7] + ftemp[8])) 
            / lattice->param.rho_in;
        c = u_y*lattice->param.rho_in;
      }
    
      ftemp[2] = ftemp[4] + (2./3.)*c;
    
      ftemp[5] = ftemp[7] + (1./2.)*( ftemp[3] - ftemp[1])
                          + (1./6.)*c;
    
      ftemp[6] = ftemp[8] + (1./2.)*( ftemp[1] - ftemp[3])
                          + (1./6.)*c; 

      ftemp += ( sizeof(struct pdf_struct)/8);

    } /* while( ftemp < ftemp_end) */

  } /* if( lattice->param.pressure_s_in[subs] ) */

  // P R E S S U R E   N O R T H   O U T F L O W   B C 
  //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  // pressure north outflow
  //  -- Pressure boundary on north side using outflow pressure condition.
  if( lattice->param.pressure_n_out[subs])
  {
//printf("bcs() %s %d >> pressure_n_out[%d]\n", __FILE__, __LINE__, subs);
    ftemp = lattice->pdf[subs][lattice->NumNodes-lattice->param.LX].ftemp;
    ftemp_end = lattice->pdf[subs][lattice->NumNodes].ftemp;
    while( ftemp < ftemp_end)
    {
      // North, Inflow
      if( lattice->param.incompressible)
      {
        u_y = -lattice->param.rho_out
            + ( ftemp[0] + ftemp[1] + ftemp[3]
              + 2.*( ftemp[2] + ftemp[5] + ftemp[6]));
        c = u_y;
      }
      else // compressible
      {
        u_y = -1.
            + ( ftemp[0] + ftemp[1] + ftemp[3]
              + 2.*( ftemp[2] + ftemp[5] + ftemp[6])) 
            / lattice->param.rho_out;
        c = u_y*lattice->param.rho_out;
      }
    
      ftemp[4] = ftemp[2] - (2./3.)*c;
    
      ftemp[7] = ftemp[5] + (1./2.)*( ftemp[1] - ftemp[3])
                          - (1./6.)*c;
    
      ftemp[8] = ftemp[6] + (1./2.)*( ftemp[3] - ftemp[1])
                          - (1./6.)*c; 

      ftemp += ( sizeof(struct pdf_struct)/8);

    } /* while( ftemp < ftemp_end) */

  } /* if( lattice->param.pressure_n_out[subs]) */

  // P R E S S U R E   S O U T H   O U T F L O W   B C 
  //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  // pressure south outflow
  //  -- Pressure boundary on south side using outflow pressure condition.
  if( lattice->param.pressure_s_out[subs])
  {
//printf("bcs() %s %d >> pressure_s_out[%d]\n", __FILE__, __LINE__, subs);
    ftemp = lattice->pdf[subs][0].ftemp;
    ftemp_end = lattice->pdf[subs][lattice->param.LX].ftemp;
    while( ftemp < ftemp_end)
    {
      // South, Outflow
      if( lattice->param.incompressible)
      {
        u_y = lattice->param.rho_out
            - ( ftemp[0] + ftemp[1] + ftemp[3]
              + 2.*( ftemp[4] + ftemp[7] + ftemp[8]));
        c = u_y;
      }
      else // compressible
      {
        u_y = 1.
            - ( ftemp[0] + ftemp[1] + ftemp[3]
              + 2.*( ftemp[4] + ftemp[7] + ftemp[8])) 
            / lattice->param.rho_out;
        c = u_y*lattice->param.rho_out;
      }
    
      ftemp[2] = ftemp[4] + (2./3.)*c;
    
      ftemp[5] = ftemp[7] + (1./2.)*( ftemp[3] - ftemp[1])
                          + (1./6.)*c;
    
      ftemp[6] = ftemp[8] + (1./2.)*( ftemp[1] - ftemp[3])
                          + (1./6.)*c; 

      ftemp += ( sizeof(struct pdf_struct)/8);

    } /* while( ftemp < ftemp_end) */

  } /* if( pressure_s_out[subs]) */


  // V E L O C I T Y   N O R T H   I N F L O W   B C 
  //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  // velocity north inflow
  //  -- Velocity boundary on north side using inflow velocity condition.
  if( lattice->param.velocity_n_in[subs])
  {
//printf("bcs() %s %d >> velocity_n_in[%d]\n", __FILE__, __LINE__, subs);
    ftemp = lattice->pdf[subs][lattice->NumNodes-lattice->param.LX].ftemp;
    ftemp_end = lattice->pdf[subs][lattice->NumNodes].ftemp;
    bc_type = &( lattice->bc[subs][lattice->NumNodes-lattice->param.LX].bc_type);
#if RHO0_TEST
//------------------------------------------------------------------[ TEST ]----
    rho0 = 
      &( lattice->macro_vars[subs][lattice->NumNodes-lattice->param.LX].rho);
//------------------------------------------------------------------[ TEST ]----
#endif /* RHO0_TEST */
    //u = ((subs==1)?(-1):(1))*lattice->param.uy_in;
    if( lattice->param.bc_poisseuille)
    {
      i = 0;
    }
    else
    {
      u = lattice->param.uy_in;
    }
    while( ftemp < ftemp_end)
    {
      // North, Inflow
     if( 1 || !( *bc_type & BC_SOLID_NODE))
     {
      //u = u_in[((double)rand()/(double)RAND_MAX<.5)?(0):(1)][subs];
      //u = 0.;
      if( lattice->param.bc_poisseuille)
      {
        u = ( 1.5*( lattice->param.uy_in)
                 /( .25*(lattice->param.LX-2)*(lattice->param.LX-2)) )
                 *( 
                    .25*( lattice->param.LX-2)*( lattice->param.LX-2) 
                  - 
                    (i-.5*( lattice->param.LX-2)-.5)
                   *(i-.5*( lattice->param.LX-2)-.5)
                  )
              ;
//printf("%s (%d) -- %d %f\n", __FILE__, __LINE__, i, u);
        i++;
      }