bcs.c

来自「基于格子Boltzmann方法开源可视化软件的源代码 具有很高的实用价值。对学」· C语言 代码 · 共 2,058 行 · 第 1/5 页

  //  -- Velocity boundary on east side using outflow velocity condition.
  if( lattice->param.velocity_e_out[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_out;
    if( lattice->param.bc_poisseuille)
    {
      i = 0;
    }
    else
    {
      u = lattice->param.ux_out;
    }
    while( ftemp < ftemp_end)
    {
      // East, 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.ux_out)
                 /( .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++;

    } /* 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_out[subs]) */

  // V E L O C I T Y   W E S T   O U T F L O W   B C 
  //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  // velocity west outflow
  //  -- Velocity boundary on west side using outflow velocity condition.
  if( lattice->param.velocity_w_out[subs])
  {
//printf("bcs() %s %d >> velocity_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;
    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_out;
    if( lattice->param.bc_poisseuille)
    {
      i = 0;
    }
    else
    {
      u = lattice->param.ux_out;
    }
    while( ftemp < ftemp_end)
    {
      // West, 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.ux_out)
                 /( .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

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

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

#if INAMURO_SIGMA_COMPONENT

  // C O N S T A N T   C O N C   N O R T H   I N F L O W   B C
  //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  // constant conc north inflow
  //  -- Constant concentration boundary on north side using inflow value.
  if( lattice->param.constcon_n_in )
  {
    //printf("%s %d -- constcon_n_in\n", __FILE__, __LINE__);
    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->time >= lattice->param.sigma_start)
      {
        c = lattice->param.rho_sigma_in;
      }
      else
      {
        c = 0.;
      }
    
      if( lattice->param.simple_diffusion)
      {
        rho = 6.*( c - ( ftemp[0] + ftemp[1] + ftemp[3] + ftemp[2]));
        ftemp[4] = (1./8.)*rho;
      }
      else
      {
        rho = 6.*( c - ( ftemp[0] + ftemp[1] + ftemp[3]
                       + ftemp[2] + ftemp[5] + ftemp[6]));
        ftemp[4] = (1./9.)*rho;
        ftemp[7] = (1./36.)*rho;
        ftemp[8] = (1./36.)*rho;
      }

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

    } /* while( ftemp < ftemp_end) */

  } /* if( lattice->param.constcon_n_in ) */

  // C O N S T A N T   C O N C   S O U T H   I N F L O W   B C 
  //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  // constant conc south inflow
  //  -- Constant concentration boundary on south side using inflow value.
  if( lattice->param.constcon_s_in )
  {
    //printf("bcs() %s %d >> constcon_s_in\n", __FILE__, __LINE__);
    ftemp = lattice->pdf[subs][0].ftemp;
    ftemp_end = lattice->pdf[subs][lattice->param.LX].ftemp;
    while( ftemp < ftemp_end)
    {
      // South, Inflow

      if( lattice->time >= lattice->param.sigma_start)
      {
        c = lattice->param.rho_sigma_in;
      }
      else
      {
        c = 0.;
      }

      if( lattice->param.simple_diffusion)
      {
        rho = 8.*( c - ( ftemp[0] + ftemp[1] + ftemp[3] + ftemp[4]));
        ftemp[2] = (1./8.)*rho;
      }
      else
      {
        rho = 6.*( c - ( ftemp[0] + ftemp[1] + ftemp[3]
                       + ftemp[7] + ftemp[4] + ftemp[8]));
        ftemp[2] = (1./9.)*rho;
        ftemp[5] = (1./36.)*rho;
        ftemp[6] = (1./36.)*rho;
      }

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

    } /* while( ftemp < ftemp_end) */

  } /* if( lattice->param.constcon_s_in ) */

  // C O N S T A N T   C O N C   N O R T H   O U T F L O W   B C 
  //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  // constant conc north outflow
  //  -- Constant concentration boundary on north side using outflow value.
  if( lattice->param.constcon_n_out)
  {
    //printf("%s %d -- constcon_n_out\n", __FILE__, __LINE__);
    ftemp = lattice->pdf[subs][lattice->NumNodes-lattice->param.LX].ftemp;
    ftemp_end = lattice->pdf[subs][lattice->NumNodes].ftemp;
    while( ftemp < ftemp_end)
    {
      // North, Outflow
      if( lattice->time >= lattice->param.sigma_start)
      {
        c = lattice->param.rho_sigma_out;
      }
      else
      {
        c = 0.;
      }

      if( lattice->param.simple_diffusion)
      {
        rho = 8.*( c - ( ftemp[0] + ftemp[1] + ftemp[3] + ftemp[2]));
        ftemp[4] = (1./8.)*rho;
      }
      else
      {
//printf("%s %d >> constcon_n_out\n",__FILE__,__LINE__);
        rho = 6.*( c - ( ftemp[0] + ftemp[1] + ftemp[3]
                       + ftemp[2] + ftemp[5] + ftemp[6]));
        ftemp[4] = (1./9.)*rho;
        ftemp[7] = (1./36.)*rho;
        ftemp[8] = (1./36.)*rho;
      }

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

    } /* while( ftemp < ftemp_end) */

  } /* if( lattice->param.constcon_n_out) */

  // C O N S T A N T   C O N C   S O U T H   O U T F L O W   B C 
  //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  // constant conc south outflow
  //  -- Constant concentration boundary on south side with outflow value.
  if( lattice->param.constcon_s_out)
  {
    //printf("bcs() %s %d >> constcon_s_out\n", __FILE__, __LINE__);
    ftemp = lattice->pdf[subs][0].ftemp;
    ftemp_end = lattice->pdf[subs][lattice->param.LX].ftemp;
    while( ftemp < ftemp_end)
    {
      // South, Outflow

      if( lattice->time >= lattice->param.sigma_start)
      {
        c = lattice->param.rho_sigma_out;
      }
      else
      {
        c = 0.;
      }

      if( lattice->param.simple_diffusion)
      {
        rho = 8.*( c - ( ftemp[0] + ftemp[1] + ftemp[3] + ftemp[4]));
        ftemp[2] = (1./8.)*rho;
      }
      else
      {
        rho = 6.*( c - ( ftemp[0] + ftemp[1] + ftemp[3]
                       + ftemp[7] + ftemp[4] + ftemp[8]));
        ftemp[2] = (1./9.)*rho;
        ftemp[5] = (1./36.)*rho;
        ftemp[6] = (1./36.)*rho;
      }

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

    } /* while( ftemp < ftemp_end) */

  } /* if( constcon_s_out) */


  // C O N S T A N T   F L U X   N O R T H   I N F L O W   B C
  //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  // constant flux north inflow
  //  -- Constant flux boundary on north side with inflow value.
  if( lattice->param.constflx_n_in )
  {
    //printf("%s %d -- constflx_n_in\n", __FILE__, __LINE__);
    ftemp = lattice->pdf[subs][lattice->NumN