bcs.c

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

      if( lattice->param.incompressible)
      {
        rho = -u + (      ftemp[0] + ftemp[1] + ftemp[3]
                   + 2.*( ftemp[2] + ftemp[5] + ftemp[6]));
        c = u;
      }
      else // compressible
      {
        rho = (      ftemp[0] + ftemp[1] + ftemp[3]
              + 2.*( ftemp[2] + ftemp[5] + ftemp[6])) 
              / ( 1. + u);
        c = rho*u;
      }
    
      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;

#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);
//------------------------------------------------------------------[ TEST ]----
#endif /* RHO0_TEST */

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

     bc_type++;

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


  }

  // V E L O C I T Y   S O U T H   I N F L O W   B C 
  //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  // velocity south inflow
  //  -- Velocity boundary on south side using inflow velocity condition.
  if( lattice->param.velocity_s_in[subs] )
  {
//printf("bcs() %s %d >> velocity_s_in[%d]\n", __FILE__, __LINE__, subs);
    ftemp = lattice->pdf[subs][0].ftemp;
    ftemp_end = lattice->pdf[subs][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.uy_in;
    if( lattice->param.bc_poisseuille)
    {
      i = 0;
    }
    else
    {
      u = lattice->param.uy_in;
    }
    while( ftemp < ftemp_end)
    {
      // South, 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)
                  )
              ;
        i++;
      }
      if( lattice->param.incompressible)
      {
        rho = u + (      ftemp[0] + ftemp[1] + ftemp[3]
                  + 2.*( ftemp[4] + ftemp[7] + ftemp[8]));
        c = u;
      }
      else
      {
        rho = ( ftemp[0] + ftemp[1] + ftemp[3]
              + 2.*( ftemp[4] + ftemp[7] + ftemp[8])) 
              / ( 1. - u);
        c = rho*u;
      }
    
      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;
#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);
//------------------------------------------------------------------[ TEST ]----
#endif /* RHO0_TEST */

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

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

     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   N O R T H   O U T F L O W   B C 
  //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  // velocity north outflow
  //  -- Velocity boundary on north side using outflow velocity condition.
  if( lattice->param.velocity_n_out[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_out;
    if( lattice->param.bc_poisseuille)
    {
      i = 0;
    }
    else
    {
      u = lattice->param.uy_out;
    }
    while( ftemp < ftemp_end)
    {
      // North, Inflow
     if( 1 || !( *bc_type & BC_SOLID_NODE))
     {
      //u = u_out[((double)rand()/(double)RAND_MAX<.5)?(0):(1)][subs];
      //u = 0.;
      if( lattice->param.bc_poisseuille)
      {
        u = ( 1.5*( lattice->param.uy_out)
                 /( .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++;
      }
      if( lattice->param.incompressible)
      {
        rho = -u + (      ftemp[0] + ftemp[1] + ftemp[3]
                   + 2.*( ftemp[2] + ftemp[5] + ftemp[6]));
        c = u;
      }
      else // compressible
      {
        rho = (      ftemp[0] + ftemp[1] + ftemp[3]
              + 2.*( ftemp[2] + ftemp[5] + ftemp[6])) 
              / ( 1. + u);
        c = rho*u;
      }
    
      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;

#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);
//------------------------------------------------------------------[ TEST ]----
#endif /* RHO0_TEST */

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

     bc_type++;

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



  }

  // V E L O C I T Y   S O U T H   O U T F L O W   B C 
  //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  // velocity south outflow
  //  -- Velocity boundary on south side using outflow velocity condition.
  if( lattice->param.velocity_s_out[subs])
  {
//printf("bcs() %s %d >> velocity_s_out[%d]\n", __FILE__, __LINE__, subs);
    ftemp = lattice->pdf[subs][0].ftemp;
    ftemp_end = lattice->pdf[subs][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.uy_out;
    if( lattice->param.bc_poisseuille)
    {
      i = 0;
    }
    else
    {
      u = lattice->param.uy_out;
    }
    while( ftemp < ftemp_end)
    {
      // South, Outflow
     if( 1 || !( *bc_type & BC_SOLID_NODE))
     {
      //u = u_out[((double)rand()/(double)RAND_MAX<.5)?(0):(1)][subs];
      //u = 0.;
      if( lattice->param.bc_poisseuille)
      {
        u = ( 1.5*( lattice->param.uy_out)
                 /( .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)
                  )
              ;
        i++;
      }
      if( lattice->param.incompressible)
      {
        rho = u + (      ftemp[0] + ftemp[1] + ftemp[3]
                  + 2.*( ftemp[4] + ftemp[7] + ftemp[8]));
        c = u;
      }
      else
      {
        rho = ( ftemp[0] + ftemp[1] + ftemp[3]
              + 2.*( ftemp[4] + ftemp[7] + ftemp[8])) 
              / ( 1. - u);
        c = rho*u;
      }
    
      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;
#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);
//------------------------------------------------------------------[ TEST ]----
#endif /* RHO0_TEST */

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

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

     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;