bcs.c

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

//##############################################################################
//
// Copyright (C), 2005, Michael Sukop and Danny Thorne
//
// bcs.c
//
//  - Boundary conditions.
//
//  - void bcs( lattice_ptr lattice);
//
//  - void process_bcs( char *filename, int **bcs);
//

#define RHO0_TEST 0

#if 0
void bcs( lattice_ptr lattice)
{
  int    n;
  int    subs;
  int    *bc;
  double *ftemp;
  double u_x, 
         u_y, 
         rho;

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

  bc = &( lattice->bc[subs][0].bc_type);
  ftemp = lattice->pdf[subs][0].ftemp;

  for( n=0; n<lattice->NumNodes; n++)
  {
    if(    *bc != BC_FLUID_NODE 
        && *bc != ( BC_FLUID_NODE | BC_FILM_NODE)
        && *bc != BC_SOLID_NODE)
    {
      if( *bc & BC_PRESSURE_N_IN ) 
      {
        // North, Inflow

//printf("bcs() -- North, Inflow at n = %d, ( %d, %d)\n", 
//  n, lattice->node[n].i, lattice->node[n].j);
//printf("  Before>> %f %f %f %f [%f] %f %f [%f] [%f]\n",
//    ftemp[0], ftemp[1], ftemp[2],
//    ftemp[3], ftemp[4], ftemp[5],
//    ftemp[6], ftemp[7], ftemp[8]  );

        u_y = -1.
            + ( ftemp[0] + ftemp[1] + ftemp[3]
              + 2.*( ftemp[2] + ftemp[5] + ftemp[6])) 
            / lattice->param.rho_in;
    
        ftemp[4] = ftemp[2] - (2./3.)*lattice->param.rho_in*u_y;
    
        ftemp[7] = ftemp[5] + (1./2.)*( ftemp[1] - ftemp[3])
                            - (1./6.)*lattice->param.rho_in*u_y;
    
        ftemp[8] = ftemp[6] + (1./2.)*( ftemp[3] - ftemp[1])
                            - (1./6.)*lattice->param.rho_in*u_y; 

//printf("  Before>> %f %f %f %f [%f] %f %f [%f] [%f]\n",
//    ftemp[0], ftemp[1], ftemp[2],
//    ftemp[3], ftemp[4], ftemp[5],
//    ftemp[6], ftemp[7], ftemp[8]  );

      } /* if( *bc & BC_PRESSURE_N_IN ) */

      if( *bc & BC_PRESSURE_S_IN ) 
      { 
        printf("bcs() -- ERROR: Support for South, Inflow pressure "
            "boundaries is pending. (Exiting!)");
        exit(1);
      }
      if( *bc & BC_PRESSURE_E_IN )
      { 
        printf("bcs() -- ERROR: Support for East, Inflow pressure "
            "boundaries is pending. (Exiting!)");
        exit(1);
      }
      if( *bc & BC_PRESSURE_W_IN )
      { 
        printf("bcs() -- ERROR: Support for West, Inflow pressure "
            "boundaries is pending. (Exiting!)");
        exit(1);
      }
      if( *bc & BC_PRESSURE_N_OUT)
      { 
        printf("bcs() -- ERROR: Support for North, Outflow pressure "
            "boundaries is pending. (Exiting!)");
        exit(1);
      }

      if( *bc & BC_PRESSURE_S_OUT) 
      { 
        // South, Outflow

//printf("bcs() -- South, Outflow at n = %d, ( %d, %d)\n",
//  n, lattice->node[n].i, lattice->node[n].j);
//printf("  Before>> %f %f [%f] %f %f [%f] [%f] %f %f\n",
//    ftemp[0], ftemp[1], ftemp[2],
//    ftemp[3], ftemp[4], ftemp[5],
//    ftemp[6], ftemp[7], ftemp[8]  );

        u_y = 1.
            - ( ftemp[0] + ftemp[1] + ftemp[3]
              + 2.*( ftemp[4] + ftemp[7] + ftemp[8])) 
            / lattice->param.rho_out;
    
        ftemp[2] = ftemp[4] + (2./3.)*lattice->param.rho_out*u_y;
    
        ftemp[5] = ftemp[7] + (1./2.)*( ftemp[3] - ftemp[1])
                            + (1./6.)*lattice->param.rho_out*u_y;
    
        ftemp[6] = ftemp[8] + (1./2.)*( ftemp[1] - ftemp[3])
                            + (1./6.)*lattice->param.rho_out*u_y; 

//printf("  Before>> %f %f [%f] %f %f [%f] [%f] %f %f\n",
//    ftemp[0], ftemp[1], ftemp[2],
//    ftemp[3], ftemp[4], ftemp[5],
//    ftemp[6], ftemp[7], ftemp[8]  );

      } /* if( *bc & BC_PRESSURE_S_OUT) */

      if( *bc & BC_PRESSURE_E_OUT)
      { 
        printf("bcs() -- ERROR: Support for East, Outflow pressure "
            "boundaries is pending. (Exiting!)");
        exit(1);
      }
      if( *bc & BC_PRESSURE_W_OUT)
      { 
        printf("bcs() -- ERROR: Support for West, Outflow pressure "
            "boundaries is pending. (Exiting!)");
        exit(1);
      }

    } /* if( *bc != 0 && *bc != BC_SOLID_NODE) */

    ftemp+=27;

    bc++;

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

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

} /* void bcs( lattice_ptr lattice) */
#else
#if 0
void bcs( lattice_ptr lattice)
{
  int    i, j, n;
  int    subs;
  int    *bc;
  double *ftemp, *ftemp_end, *ftemp_mid;
  double u_x,
         u_y;
  double u_in[2][2], 
         u_out[2][2], 
         u,
         rho;

#if 0
 u_in[0][0]  =  lattice->param.uy_in;
 u_in[1][0]  =  lattice->param.uy_in;
 u_in[0][1]  = -lattice->param.uy_in;
 u_in[1][1]  = -lattice->param.uy_in;

 u_out[0][0] =  lattice->param.uy_out;
 u_out[1][0] =  lattice->param.uy_out;
 u_out[0][1] = -lattice->param.uy_out;
 u_out[1][1] = -lattice->param.uy_out;
#else
 u_in[0][0]  =  0.;
 u_in[1][0]  =  lattice->param.uy_in;
 u_in[0][1]  = -lattice->param.uy_in;
 u_in[1][1]  =  0.;

 u_out[0][0] =  lattice->param.uy_out;
 u_out[1][0] =  0.;
 u_out[0][1] =  0.;
 u_out[1][1] = -lattice->param.uy_out;
#endif

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

#if 0
  bc = &( lattice->bc[subs][0].bc_type);
  ftemp = lattice->pdf[subs][0].ftemp;

  for( n=0; n<lattice->NumNodes; n++)
  {
    if(    *bc != BC_FLUID_NODE 
        && *bc != ( BC_FLUID_NODE | BC_FILM_NODE)
        && *bc != BC_SOLID_NODE)
    {
      if( *bc & BC_PRESSURE_N_IN ) 
      {
        // North, Inflow
        u_y = -1.
            + ( ftemp[0] + ftemp[1] + ftemp[3]
              + 2.*( ftemp[2] + ftemp[5] + ftemp[6])) 
            / lattice->param.rho_in;
    
        ftemp[4] = ftemp[2] - (2./3.)*lattice->param.rho_in*u_y;
    
        ftemp[7] = ftemp[5] + (1./2.)*( ftemp[1] - ftemp[3])
                            - (1./6.)*lattice->param.rho_in*u_y;
    
        ftemp[8] = ftemp[6] + (1./2.)*( ftemp[3] - ftemp[1])
                            - (1./6.)*lattice->param.rho_in*u_y; 

      } /* if( *bc & BC_PRESSURE_N_IN ) */

      if( *bc & BC_PRESSURE_S_IN ) 
      { 
        printf("bcs() -- ERROR: Support for South, Inflow pressure "
            "boundaries is pending. (Exiting!)");
        exit(1);
      }
      if( *bc & BC_PRESSURE_E_IN )
      { 
        printf("bcs() -- ERROR: Support for East, Inflow pressure "
            "boundaries is pending. (Exiting!)");
        exit(1);
      }
      if( *bc & BC_PRESSURE_W_IN )
      { 
        printf("bcs() -- ERROR: Support for West, Inflow pressure "
            "boundaries is pending. (Exiting!)");
        exit(1);
      }
      if( *bc & BC_PRESSURE_N_OUT)
      { 
        printf("bcs() -- ERROR: Support for North, Outflow pressure "
            "boundaries is pending. (Exiting!)");
        exit(1);
      }

      if( *bc & BC_PRESSURE_S_OUT) 
      { 
        // South, Outflow
        u_y = 1.
            - ( ftemp[0] + ftemp[1] + ftemp[3]
              + 2.*( ftemp[4] + ftemp[7] + ftemp[8])) 
            / lattice->param.rho_out;
    
        ftemp[2] = ftemp[4] + (2./3.)*lattice->param.rho_out*u_y;
    
        ftemp[5] = ftemp[7] + (1./2.)*( ftemp[3] - ftemp[1])
                            + (1./6.)*lattice->param.rho_out*u_y;
    
        ftemp[6] = ftemp[8] + (1./2.)*( ftemp[1] - ftemp[3])
                            + (1./6.)*lattice->param.rho_out*u_y; 

      } /* if( *bc & BC_PRESSURE_S_OUT) */

      if( *bc & BC_PRESSURE_E_OUT)
      { 
        printf("bcs() -- ERROR: Support for East, Outflow pressure "
            "boundaries is pending. (Exiting!)");
        exit(1);
      }
      if( *bc & BC_PRESSURE_W_OUT)
      { 
        printf("bcs() -- ERROR: Support for West, Outflow pressure "
            "boundaries is pending. (Exiting!)");
        exit(1);
      }

    } /* if( *bc != 0 && *bc != BC_SOLID_NODE) */

    ftemp+=27;

    bc++;

  } /* for( n=0; n<lattice->NumNodes; n++) */
#else
  // 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!
  if( lattice->param.pressure_n_in[subs] )
  {
//printf("bcs() -- pressure_n_in[%d]\n", 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
      u_y = -1.
          + ( ftemp[0] + ftemp[1] + ftemp[3]
            + 2.*( ftemp[2] + ftemp[5] + ftemp[6])) 
          / lattice->param.rho_in;
    
      ftemp[4] = ftemp[2] - (2./3.)*lattice->param.rho_in*u_y;
    
      ftemp[7] = ftemp[5] + (1./2.)*( ftemp[1] - ftemp[3])
                          - (1./6.)*lattice->param.rho_in*u_y;
    
      ftemp[8] = ftemp[6] + (1./2.)*( ftemp[3] - ftemp[1])
                          - (1./6.)*lattice->param.rho_in*u_y; 

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

    } /* while( ftemp < ftemp_end) */
  } /* if( lattice->param.pressure_n_in[subs] ) */

  if( lattice->param.pressure_s_in[subs] )
  {
  }
  if( lattice->param.pressure_n_out[subs])
  {
  }
  if( lattice->param.pressure_s_out[subs])
  {
//printf("bcs() -- pressure_s_out[%d]\n", subs);
    ftemp = lattice->pdf[subs][0].ftemp;
    ftemp_end = lattice->pdf[subs][lattice->param.LX].ftemp;
    while( ftemp < ftemp_end)
    {
      u_y = 1.
          - ( ftemp[0] + ftemp[1] + ftemp[3]
            + 2.*( ftemp[4] + ftemp[7] + ftemp[8])) 
          / lattice->param.rho_out;
    
      ftemp[2] = ftemp[4] + (2./3.)*lattice->param.rho_out*u_y;
    
      ftemp[5] = ftemp[7] + (1./2.)*( ftemp[3] - ftemp[1])
                          + (1./6.)*lattice->param.rho_out*u_y;
    
      ftemp[6] = ftemp[8] + (1./2.)*( ftemp[1] - ftemp[3])
                          + (1./6.)*lattice->param.rho_out*u_y; 

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

    } /* while( ftemp < ftemp_end) */

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

  if( lattice->param.velocity_n_in[subs] )
  {
//printf("bcs() -- velocity_n_in[%d]\n", subs);
    ftemp = lattice->pdf[subs][lattice->NumNodes-lattice->param.LX].ftemp;
ftemp_mid = lattice->pdf[subs][lattice->NumNodes-7*lattice->param.LX/8].ftemp;
    ftemp_end = lattice->pdf[subs][lattice->NumNodes].ftemp;
    //u = ((subs==1)?(-1):(1))*lattice->param.uy_in;
    while( ftemp < ftemp_mid)
    {
      // 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) */
ftemp_mid = lattice->pdf[subs][lattice->NumNodes-3*lattice->param.LX/4].ftemp;
    while( ftemp < ftemp_mid)
    {
      // North, Inflow
      //u = u_in[((double)rand()/(double)RAND_MAX<.5)?(0):(1)][subs];
      u = u_in[0][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) */
ftemp_mid = lattice->pdf[subs][lattice->NumNodes-lattice->param.LX/2].ftemp;
    while( ftemp < ftemp_mid)
    {
      // 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) */
ftemp_mid = lattice->pdf[subs][lattice->NumNodes-3*lattice->param.LX/8].ftemp;