bcs.c

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

      ftemp[6] = ( ftemp[-27+6] + ftemp[27*lattice->param.LX + 6]) /2;

    }
#endif

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

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

      ftemp += ( sizeof(struct pdf_struct)/8)*lattice->param.LX;

    } /* while( ftemp < ftemp_end) */

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

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

      ftemp += ( sizeof(struct pdf_struct)/8)*lattice->param.LX;

    } /* while( ftemp < ftemp_end) */

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

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

      ftemp += ( sizeof(struct pdf_struct)/8)*lattice->param.LX;

    } /* while( ftemp < ftemp_end) */

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

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

      ftemp += ( sizeof(struct pdf_struct)/8)*lattice->param.LX;

    } /* while( ftemp < ftemp_end) */

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

  // V E L O C I T Y   E A S T   I N F L O W   B C 
  //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  // velocity east inflow
  //  -- Velocity boundary on east side using inflow velocity condition.
  if( lattice->param.velocity_e_in[subs])
  {
//printf("bcs() %s %d >> velocity_e_in[%d]\n", __FILE__, __LINE__, subs);
    ftemp = lattice->pdf[subs][lattice->param.LX-1].ftemp;
    ftemp_end = lattice->pdf[subs][lattice->NumNodes].ftemp;
    bc_type = &( lattice->bc[subs][lattice->param.LX-1].bc_type);
#if RHO0_TEST
//------------------------------------------------------------------[ TEST ]----
    rho0 = 
      &( lattice->macro_vars[subs][lattice->param.LX-1].rho);
//------------------------------------------------------------------[ TEST ]----
#endif /* RHO0_TEST */
    //u = ((subs==1)?(-1):(1))*lattice->param.ux_in;
    if( lattice->param.bc_poisseuille)
    {
      i = 0;
    }
    else
    {
      u = lattice->param.ux_in;
    }
    while( ftemp < ftemp_end)
    {
      // East, 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.ux_in)
                 /( .25*(lattice->param.LY-2)*(lattice->param.LY-2)) )
                 *( 
                    .25*( lattice->param.LY-2)*( lattice->param.LY-2) 
                  - 
                    (i-.5*( lattice->param.LY-2)-.5)
                   *(i-.5*( lattice->param.LY-2)-.5)
                  )
              ;
//printf("%s (%d) -- %d %f\n", __FILE__, __LINE__, i, u);
        i++;
      }
      if( lattice->param.incompressible)
      {
        rho = -u + (      ftemp[0] + ftemp[2] + ftemp[4]
                   + 2.*( ftemp[1] + ftemp[5] + ftemp[8]));
        c = u;
      }
      else // compressible
      {
        rho = (      ftemp[0] + ftemp[2] + ftemp[4]
              + 2.*( ftemp[1] + ftemp[5] + ftemp[8])) 
              / ( 1. + u);
        c = rho*u;
      }
    
      ftemp[3] = ftemp[1] - (2./3.)*c;
    
      ftemp[7] = ftemp[5] + (1./2.)*( ftemp[2] - ftemp[4])
                          - (1./6.)*c;
    
      ftemp[6] = ftemp[8] + (1./2.)*( ftemp[4] - ftemp[2])
                          - (1./6.)*c;

#if RHO0_TEST
//------------------------------------------------------------------[ TEST ]----
      ftemp[0] = *rho0 - ( ftemp[1]
                         + ftemp[2]
                         + ftemp[3]
                         + ftemp[4]
                         + ftemp[5]
                         + ftemp[6]
                         + ftemp[7]
                         + ftemp[8]);
      rho0 += ( sizeof(struct macro_vars_struct)/8)*lattice->param.LX;
//------------------------------------------------------------------[ TEST ]----
#endif /* RHO0_TEST */

      ftemp += ( sizeof(struct pdf_struct)/8)*lattice->param.LX;
      
     }
     else
     {
        if( lattice->param.bc_poisseuille) { i++;}
#if RHO0_TEST
//------------------------------------------------------------------[ TEST ]----
        rho0 += ( sizeof(struct macro_vars_struct)/8)*lattice->param.LX;
//------------------------------------------------------------------[ TEST ]----
#endif /* RHO0_TEST */
        ftemp += ( sizeof(struct pdf_struct)/8)*lattice->param.LX;
     }

     bc_type+=lattice->param.LX;

    } /* while( ftemp < ftemp_end) */

#if 0
    if( lattice->param.bc_poisseuille)
    {
      // Fix corners
      ftemp = lattice->pdf[subs][lattice->NumNodes-lattice->param.LX].ftemp;
      ftemp[4] = ( ftemp[27+4] + ftemp[-27*lattice->param.LX + 4]) /2;
      ftemp[8] = ( ftemp[27+8] + ftemp[-27*lattice->param.LX + 8]) /2;

      ftemp = lattice->pdf[subs][lattice->NumNodes-1].ftemp;
      ftemp[4] = ( ftemp[-27+4] + ftemp[-27*lattice->param.LX + 4]) /2;
      ftemp[7] = ( ftemp[-27+7] + ftemp[-27*lattice->param.LX + 7]) /2;

    }
#endif

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

  // V E L O C I T Y   W E S T   I N F L O W   B C 
  //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  // velocity west inflow
  //  -- Velocity boundary on west side using inflow velocity condition.
  if( lattice->param.velocity_w_in[subs] )
  {
//printf("bcs() %s %d >> velocity_w_in[%d]\n", __FILE__, __LINE__, subs);
    ftemp     = lattice->pdf[subs][0].ftemp;
    ftemp_end = 
      lattice->pdf[subs][(lattice->param.LY+1)*lattice->param.LX].ftemp;
    bc_type = &( lattice->bc[subs][0].bc_type);
#if RHO0_TEST
//------------------------------------------------------------------[ TEST ]----
    rho0 = &( lattice->macro_vars[subs][0].rho);
//------------------------------------------------------------------[ TEST ]----
#endif /* RHO0_TEST */
    //u = ((subs==1)?(-1):(1))*lattice->param.ux_in;
    if( lattice->param.bc_poisseuille)
    {
      i = 0;
    }
    else
    {
      u = lattice->param.ux_in;
    }
    while( ftemp < ftemp_end)
    {
      // West, 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.ux_in)
                 /( .25*( lattice->param.LY-2)*( lattice->param.LY-2)) )
                 *( 
                    .25*( lattice->param.LY-2)*( lattice->param.LY-2) 
                  - 
                    (i-.5*( lattice->param.LY-2)-.5)
                   *(i-.5*( lattice->param.LY-2)-.5)
                  )
              ;
        i++;
      }
      if( lattice->param.incompressible)
      {
        rho = u + (      ftemp[0] + ftemp[2] + ftemp[4]
                  + 2.*( ftemp[3] + ftemp[7] + ftemp[6]));
        c = u;
      }
      else
      {
        rho = (      ftemp[0] + ftemp[2] + ftemp[4]
              + 2.*( ftemp[3] + ftemp[7] + ftemp[6])) 
              / ( 1. - u);
        c = rho*u;
      }
    
      ftemp[1] = ftemp[3] + (2./3.)*c;
    
      ftemp[5] = ftemp[7] + (1./2.)*( ftemp[4] - ftemp[2])
                          + (1./6.)*c;
    
      ftemp[8] = ftemp[6] + (1./2.)*( ftemp[2] - ftemp[4])
                          + (1./6.)*c;
#if RHO0_TEST
//------------------------------------------------------------------[ TEST ]----
      ftemp[0] = *rho0 - ( ftemp[1]
                         + ftemp[2]
                         + ftemp[3]
                         + ftemp[4]
                         + ftemp[5]
                         + ftemp[6]
                         + ftemp[7]
                         + ftemp[8]);
      rho0 += ( sizeof(struct macro_vars_struct)/8)*lattice->param.LX;
//------------------------------------------------------------------[ TEST ]----
#endif /* RHO0_TEST */

      ftemp += ( sizeof(struct pdf_struct)/8)*lattice->param.LX;

     }
     else
     {
        if( lattice->param.bc_poisseuille) { i++;}
#if RHO0_TEST
//------------------------------------------------------------------[ TEST ]----
        rho0 += ( sizeof(struct macro_vars_struct)/8)*lattice->param.LX;
//------------------------------------------------------------------[ TEST ]----
#endif /* RHO0_TEST */
        ftemp += ( sizeof(struct pdf_struct)/8)*lattice->param.LX;
     }

     bc_type++;

    } /* while( ftemp < ftemp_end) */

#if 0
    if( lattice->param.bc_poisseuille)
    {
      // Fix corners
      ftemp = lattice->pdf[subs][0].ftemp;
      ftemp[2] = ( ftemp[27+2] + ftemp[27*lattice->param.LX + 2]) /2;
      ftemp[5] = ( ftemp[27+5] + ftemp[27*lattice->param.LX + 5]) /2;

      ftemp = lattice->pdf[subs][lattice->param.LX-1].ftemp;
      ftemp[2] = ( ftemp[-27+2] + ftemp[27*lattice->param.LX + 2]) /2;
      ftemp[6] = ( ftemp[-27+6] + ftemp[27*lattice->param.LX + 6]) /2;

    }
#endif

  }

  // V E L O C I T Y   E A S T   O U T F L O W   B C 
  //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  // velocity east outflow