compute.c

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

      }
      else
      {
        psi[j][i] = 0.;
      }

      rho+=3;

    } /* for( i=0; i<ni; i++) */
  } /* for( j=0; j<nj; j++) */

  if( lattice->periodic_y[subs])
  {
    sj = 0;
    ej = nj-1;
    force = lattice->force[subs][0].force;
  }
  else
  //{
  //  sj = 1;
  //  ej = nj-2;
  //  force = lattice->force[subs][ni].force;
  //}
  {
    j = 0;
    jp = j+1;
    force = lattice->force[subs][0].force;
    for( i=0; i<ni; i++)
    {
      ip = ( i<ni-1)?( i+1):( 0   );
      in = ( i>0   )?( i-1):( ni-1);

      *( force  ) = 0.;
      *( force+1) = 0.;

      if( !( lattice->bc[subs][ j*ni+i].bc_type & BC_SOLID_NODE))
      {
        if( !( lattice->bc[subs][ j *ni+ip].bc_type & BC_SOLID_NODE)) {
          *( force  ) += WM*vx[1]*psi[j ][ip];
          *( force+1) += WM*vy[1]*psi[j ][ip]; }
        if( !( lattice->bc[subs][ jp*ni+i ].bc_type & BC_SOLID_NODE)) {
          *( force  ) += WM*vx[2]*psi[jp][i ];
          *( force+1) += WM*vy[2]*psi[jp][i ]; }
        if( !( lattice->bc[subs][ j *ni+in].bc_type & BC_SOLID_NODE)) {
          *( force  ) += WM*vx[3]*psi[j ][in];
          *( force+1) += WM*vy[3]*psi[j ][in]; }
/**/    if( !( lattice->bc[subs][ j *ni+i ].bc_type & BC_SOLID_NODE)) {
/**/      *( force  ) += WM*vx[4]*psi[j ][i ];
/**/      *( force+1) += WM*vy[4]*psi[j ][i ]; }
        if( !( lattice->bc[subs][ jp*ni+ip].bc_type & BC_SOLID_NODE)) {
          *( force  ) += WD*vx[5]*psi[jp][ip];
          *( force+1) += WD*vy[5]*psi[jp][ip]; }
        if( !( lattice->bc[subs][ jp*ni+in].bc_type & BC_SOLID_NODE)) {
          *( force  ) += WD*vx[6]*psi[jp][in];
          *( force+1) += WD*vy[6]*psi[jp][in]; }
/**/    if( !( lattice->bc[subs][ j *ni+in].bc_type & BC_SOLID_NODE)) {
/**/      *( force  ) += WD*vx[7]*psi[j ][in];
/**/      *( force+1) += WD*vy[7]*psi[j ][in]; }
/**/    if( !( lattice->bc[subs][ j *ni+ip].bc_type & BC_SOLID_NODE)) {
/**/      *( force  ) += WD*vx[8]*psi[j ][ip];
/**/      *( force+1) += WD*vy[8]*psi[j ][ip]; }

          *( force  ) = -lattice->param.G*psi[j][i]*( *(force  ));
          *( force+1) = -lattice->param.G*psi[j][i]*( *(force+1));

      } /* if( !( lattice->bc[subs][n].bc_type & BC_SOLID_NODE)) */

      else
      {
          *( force  ) = 0.;
          *( force+1) = 0.;
      }

      force += ( sizeof( struct force_struct)/8);

    } /* for( i=0; i<ni; i++) */

    j = nj-1;
    jn = j-1;
    force = lattice->force[subs][j*ni].force;
    for( i=0; i<ni; i++)
    {
      ip = ( i<ni-1)?( i+1):( 0   );
      in = ( i>0   )?( i-1):( ni-1);

      *( force  ) = 0.;
      *( force+1) = 0.;

      if( !( lattice->bc[subs][ j*ni+i].bc_type & BC_SOLID_NODE))
      {
        if( !( lattice->bc[subs][ j *ni+ip].bc_type & BC_SOLID_NODE)) {
          *( force  ) += WM*vx[1]*psi[j ][ip];
          *( force+1) += WM*vy[1]*psi[j ][ip]; }
/**/    if( !( lattice->bc[subs][ j *ni+i ].bc_type & BC_SOLID_NODE)) {
/**/      *( force  ) += WM*vx[2]*psi[j ][i ];
/**/      *( force+1) += WM*vy[2]*psi[j ][i ]; }
        if( !( lattice->bc[subs][ j *ni+in].bc_type & BC_SOLID_NODE)) {
          *( force  ) += WM*vx[3]*psi[j ][in];
          *( force+1) += WM*vy[3]*psi[j ][in]; }
        if( !( lattice->bc[subs][ jn*ni+i ].bc_type & BC_SOLID_NODE)) {
          *( force  ) += WM*vx[4]*psi[jn][i ];
          *( force+1) += WM*vy[4]*psi[jn][i ]; }
/**/    if( !( lattice->bc[subs][ j *ni+ip].bc_type & BC_SOLID_NODE)) {
/**/      *( force  ) += WD*vx[5]*psi[j ][ip];
/**/      *( force+1) += WD*vy[5]*psi[j ][ip]; }
/**/    if( !( lattice->bc[subs][ j *ni+in].bc_type & BC_SOLID_NODE)) {
/**/      *( force  ) += WD*vx[6]*psi[j ][in];
/**/      *( force+1) += WD*vy[6]*psi[j ][in]; }
        if( !( lattice->bc[subs][ jn*ni+in].bc_type & BC_SOLID_NODE)) {
          *( force  ) += WD*vx[7]*psi[jn][in];
          *( force+1) += WD*vy[7]*psi[jn][in]; }
        if( !( lattice->bc[subs][ jn*ni+ip].bc_type & BC_SOLID_NODE)) {
          *( force  ) += WD*vx[8]*psi[jn][ip];
          *( force+1) += WD*vy[8]*psi[jn][ip]; }

          *( force  ) = -lattice->param.G*psi[j][i]*( *(force  ));
          *( force+1) = -lattice->param.G*psi[j][i]*( *(force+1));

      } /* if( !( lattice->bc[subs][n].bc_type & BC_SOLID_NODE)) */

      else
      {
          *( force  ) = 0.;
          *( force+1) = 0.;
      }

      force += ( sizeof( struct force_struct)/8);

    } /* for( i=0; i<ni; i++) */

    sj = 1;
    ej = nj-2;
    force = lattice->force[subs][ni].force;
  }
  for( j=sj; j<=ej; j++)
  {
    jp = ( j<nj-1)?( j+1):( 0   );
    jn = ( j>0   )?( j-1):( nj-1);

    for( i=0; i<ni; i++)
    {
      ip = ( i<ni-1)?( i+1):( 0   );
      in = ( i>0   )?( i-1):( ni-1);

      *( force  ) = 0.;
      *( force+1) = 0.;

      if( !( lattice->bc[subs][ j*ni+i].bc_type & BC_SOLID_NODE))
      {
/* 1 */ if( !( lattice->bc[subs][ j *ni+ip].bc_type & BC_SOLID_NODE)) {
          *( force  ) += WM*vx[1]*psi[j ][ip];
          *( force+1) += WM*vy[1]*psi[j ][ip]; }
/* 2 */ if( !( lattice->bc[subs][ jp*ni+i ].bc_type & BC_SOLID_NODE)) {
          *( force  ) += WM*vx[2]*psi[jp][i ];
          *( force+1) += WM*vy[2]*psi[jp][i ]; }
/* 3 */ if( !( lattice->bc[subs][ j *ni+in].bc_type & BC_SOLID_NODE)) {
          *( force  ) += WM*vx[3]*psi[j ][in];
          *( force+1) += WM*vy[3]*psi[j ][in]; }
/* 4 */ if( !( lattice->bc[subs][ jn*ni+i ].bc_type & BC_SOLID_NODE)) {
          *( force  ) += WM*vx[4]*psi[jn][i ];
          *( force+1) += WM*vy[4]*psi[jn][i ]; }
/* 5 */ if( !( lattice->bc[subs][ jp*ni+ip].bc_type & BC_SOLID_NODE)) {
          *( force  ) += WD*vx[5]*psi[jp][ip];
          *( force+1) += WD*vy[5]*psi[jp][ip]; }
/* 6 */ if( !( lattice->bc[subs][ jp*ni+in].bc_type & BC_SOLID_NODE)) {
          *( force  ) += WD*vx[6]*psi[jp][in];
          *( force+1) += WD*vy[6]*psi[jp][in]; }
/* 7 */ if( !( lattice->bc[subs][ jn*ni+in].bc_type & BC_SOLID_NODE)) {
          *( force  ) += WD*vx[7]*psi[jn][in];
          *( force+1) += WD*vy[7]*psi[jn][in]; }
/* 8 */ if( !( lattice->bc[subs][ jn*ni+ip].bc_type & BC_SOLID_NODE)) {
          *( force  ) += WD*vx[8]*psi[jn][ip];
          *( force+1) += WD*vy[8]*psi[jn][ip]; }

          *( force  ) = -lattice->param.G*psi[j][i]*( *(force  ));
          *( force+1) = -lattice->param.G*psi[j][i]*( *(force+1));

      } /* if( !( lattice->bc[subs][n].bc_type & BC_SOLID_NODE)) */

      else
      {
          *( force  ) = 0.;
          *( force+1) = 0.;
      }

      force += ( sizeof( struct force_struct)/8);

    } /* for( i=0; i<ni; i++) */
  } /* for( j=0; j<nj; j++) */

  for( j=0; j<nj; j++)
  {
    free( psi[j]);
  }
  free( psi);

} /* void compute_phase_force( lattice_ptr lattice) */
#endif /* ZHANG_AND_CHEN_ENERGY_TRANSPORT */

//void compute_double_fluid_solid_force( lattice_ptr lattice)
//##############################################################################
//
// C O M P U T E   F L U I D   S O L I D   F O R C E
//
void compute_double_fluid_solid_force( lattice_ptr lattice)
{
  // Eq. 20 of Martys and Chen, 1996

  // Declare local variables.
  double sum_x, 
         sum_y;
  int    x, y, 
         xn, yn, 
         xp, yp;

  int subs;

  int LX = lattice->param.LX;
  int LY = lattice->param.LY;

  for( y = 0; y < LY; y++)
  {
    yp = ( y<LY-1)?( y+1):( 0   );
    yn = ( y>0   )?( y-1):( LY-1);

    for( x = 0; x < LX; x++)
    {
      xp = ( x<LX-1)?( x+1):( 0   );
      xn = ( x>0   )?( x-1):( LX-1);

      if( !( lattice->bc[0][ y*LX + x].bc_type & BC_SOLID_NODE)) 
      {
        sum_x=0.;
        sum_y=0.;

        // neighbor 1
        //if( b[y][xp]) 
        if( lattice->bc[0][ y*LX + xp].bc_type & BC_SOLID_NODE)
        {
          sum_x = 2.*vx[1] ;
          sum_y = 2.*vy[1] ;
        }
        // neighbor 2
        //if( b[yp][x]) 
        if( lattice->bc[0][ yp*LX + x].bc_type & BC_SOLID_NODE)
        {
          sum_x = sum_x + 2.*vx[2] ;
          sum_y = sum_y + 2.*vy[2] ;
        }
        // neighbor 3
        //if( b[y][xn]) 
        if( lattice->bc[0][ y*LX + xn].bc_type & BC_SOLID_NODE)
        {
          sum_x = sum_x + 2.*vx[3] ;
          sum_y = sum_y + 2.*vy[3] ;
        }
        // neighbor 4
        //if( b[yn][x]) 
        if( lattice->bc[0][ yn*LX + x].bc_type & BC_SOLID_NODE)
        {
          sum_x = sum_x + 2.*vx[4] ;
          sum_y = sum_y + 2.*vy[4] ;
        }
        // neighbor 5
        //if( b[yp][xp]) 
        if( lattice->bc[0][ yp*LX + xp].bc_type & BC_SOLID_NODE)
        {
          sum_x = sum_x + vx[5] ;
          sum_y = sum_y + vy[5] ;
        }
        // neighbor 6
        //if( b[yp][xn]) 
        if( lattice->bc[0][ yp*LX + xn].bc_type & BC_SOLID_NODE)
        {
          sum_x = sum_x + vx[6] ;
          sum_y = sum_y + vy[6] ;
        }
        // neighbor 7
        //if( b[yn][xn]) 
        if( lattice->bc[0][ yn*LX + xn].bc_type & BC_SOLID_NODE)
        {
          sum_x = sum_x + vx[7] ;
          sum_y = sum_y + vy[7] ;
        }
        // neighbor 8
        //if( b[yn][xp]) 
        if( lattice->bc[0][ yn*LX + xp].bc_type & BC_SOLID_NODE)
        {
          sum_x = sum_x + vx[8] ;
          sum_y = sum_y + vy[8] ;
        }

        for( subs=0; subs<NUM_FLUID_COMPONENTS; subs++)
        {
          if( lattice->macro_vars[subs][y*LX+x].rho != 0)
          {
            lattice->force[ subs][ y*LX+x].sforce[0] 
              = -lattice->param.Gads[subs]*sum_x;
            lattice->force[ subs][ y*LX+x].sforce[1] 
              = -lattice->param.Gads[subs]*sum_y;
          }
          else
          {
            lattice->force[ subs][ y*LX+x].sforce[0] = 0.;
            lattice->force[ subs][ y*LX+x].sforce[1] = 0.;
          }
        }

      } /* if( !obst[y][x]) */
    } /* for( x = 1; x <= LX; x++) */
  } /* for( y = 1; y <= LY; y++) */

} /* void compute_double_fluid_solid_force( lattice_ptr lattice) */


//void compute_single_fluid_solid_force( lattice_ptr lattice)
//##############################################################################
//
// C O M P U T E   F L U I D   S O L I D   F O R C E
//
void compute_single_fluid_solid_force( lattice_ptr lattice, int subs)
{
  // Eq. 20 of Martys and Chen, 1996

  // Declare local variables.
  double sum_x, 
         sum_y;
  int    x, y, 
         xn, yn, 
         xp, yp;

  int LX = lattice->param.LX;
  int LY = lattice->param.LY;

  int ni = lattice->param.LX;
  int nj = lattice->param.LY;

  int i, j, si, sj, ei, ej;

  int n;

  double **psi;

  double *rho;

  psi = ( double**)malloc( nj*sizeof(double*));
  for( j=0; j<nj; j++)
  {
    psi[j] = ( double*)malloc( ni*sizeof(double));
  }

  // Initialize psi.
  for( j=0; j<nj; j++)
  {
    for( i=0; i<ni; i++)
    {
      psi[j][i] = 0.;
    }
  }

  if( lattice->periodic_y[subs])
  {
    sj = 0;
    ej = nj-1;
    rho = &( lattice->macro_vars[subs][0].rho);
  }
  else
  {
    sj = 1;
    ej = nj-2;
    rho = &( lattice->macro_vars[subs][ni].rho);
  }

  for( j=sj; j<=ej; j++)
  {
    n = j*ni;
    for( i=0; i<ni; i++, n++)
    {
      if( *rho != 0 && !( lattice->bc[subs][n].bc_type & BC_SOLID_NODE))
      {
        psi[j][i] = 4.*exp(-200./(*rho));
      }
      else
      {
        psi[j][i] = 0.;
      }

      rho+=3;

    } /* for( i=0; i<ni; i++) */
  } /* for( j=0; j<nj; j++) */

  //(MOSIX) write(6,*) "Sforce", LX,LY

  //rho = &( lattice->macro_vars[subs][0].rho);
  for( y = 0; y < LY; y++)
  {
    yp = ( y<LY-1)?( y+1):( 0   );
    yn = ( y>0   )?( y-1):( LY-1);

    for( x = 0; x < LX; x++)//, rho+=3)
    {
      xp = ( x<LX-1)?( x+1):( 0   );
      xn = ( x>0   )?( x-1):( LX-1);

      if( !( lattice->bc[0][ y*LX + x].bc_type & BC_SOLID_NODE)) 
      {
        sum_x=0.;
        sum_y=0.;

        // neighbor 1
        //if( b[y][xp]) 
        if( lattice->bc[0][ y*LX + xp].bc_type & BC_SOLID_NODE)
        {
          sum_x = WM*vx[1] ;
          sum_y = WM*vy[1] ;
        }
        // neighbor 2
        //if( b[yp][x]) 
        if( lattice->bc[0][ yp*LX + x].bc_type & BC_SOLID_NODE)
        {
          sum_x = sum_x + WM*vx[2] ;
          sum_y = sum_y + WM*vy[2] ;
        }
        // neighbor 3
        //if( b[y][xn]) 
        if( lattice->bc[0][ y*LX + xn].bc_type & BC_SOLID_NODE)
        {
          sum_x = sum_x + WM*vx[3] ;
          sum_y = sum_y + WM*vy[3] ;
        }
        // neighbor 4
        //if( b[yn][x]) 
        if( lattice->bc[0][ yn*LX + x].bc_type & BC_SOLID_NODE)
        {
          sum_x = sum_x + WM*vx[4] ;
          sum_y = sum_y + WM*vy[4] ;
        }
        // neighbor 5
        //if( b[yp][xp]) 
        if( lattice->bc[0][ yp*LX + xp].bc_type & BC_SOLID_NODE)
        {
          sum_x = sum_x + WD*vx[5] ;
          sum_y = sum_y + WD*vy[5] ;
        }
        // neighbor 6
        //if( b[yp][xn]) 
        if( lattice->bc[0][ yp*LX + xn].bc_type & BC_SOLID_NODE)
        {
          sum_x = sum_x + WD*vx[6] ;
          sum_y = sum_y + WD*vy[6] ;
        }
        // neighbor 7
        //if( b[yn][xn]) 
        if( lattice->bc[0][ yn*LX + xn].bc_type & BC_SOLID_NODE)
        {
          sum_x = sum_x + WD*vx[7] ;
          sum_y = sum_y + WD*vy[7] ;
        }
        // neighbor 8
        //if( b[yn][xp]) 
        if( lattice->bc[0][ yn*LX + xp].bc_type & BC_SOLID_NODE)
        {
          sum_x = sum_x + WD*vx[8] ;
          sum_y = sum_y + WD*vy[8] ;
        }

        if( lattice->macro_vars[subs][y*LX+x].rho != 0)
        {
          lattice->force[ subs][ y*LX+x].sforce[0] 
            = -lattice->param.Gads[subs]*psi[y][x]*sum_x;///(*rho);
          lattice->force[ subs][ y*LX+x].sforce[1] 
            = -lattice->param.Gads[subs]*psi[y][x]*sum_y;///(*rho);
        }
        else
        {
          lattice->force[ subs][ y*LX+x].sforce[0] = 0.;
          lattice->force[ subs][ y*LX+x].sforce[1] = 0.;
        }

      } /* if( !obst[y][x]) */
    } /* for( x = 1; x <= LX; x++) */
  } /* for( y = 1; y <= LY; y++) */

  for( j=0; j<nj; j++)
  {
    free( psi[j]);
  }
  free( psi);

} /* void compute_single_fluid_solid_force( lattice_ptr lattice) */

#endif /* NON_LOCAL_FORCES */

//void compute_max_rho( lattice_ptr lattice, double *max_rho, int subs)
void compute_max_rho( lattice_ptr lattice, double *max_rho, int subs)
{
  int n;
  *max_rho = 0.;
  for( n=0; n<lattice->NumNodes; n++)
  {
    if( !( lattice->bc[subs][n].bc_type & BC_SOLID_NODE))