compute.c

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

    else
    {
      sj = 1;
      ej = nj-2;
      force[subs] = 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);

//printf("compute_fluid_fluid_force() -- "
//    "subs %d, ( i, j) = ( %2d, %2d), | f - f0| = %d\n", 
//    subs, i, j,
//    force[subs] - lattice->force[subs][0].force );

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

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

        force[subs] += ( sizeof( struct force_struct)/8);

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

  force[0] = lattice->force[0][0].force;
  force[1] = lattice->force[1][0].force;
  for( j=0; j<nj; j++)
  {
    for( i=0; i<ni; i++)
    {
      if( !( lattice->bc[0][ j*ni+i].bc_type & BC_SOLID_NODE))
      {
        psi_temp = *( force[1]  );
        *( force[1]  ) = -lattice->param.G*psi[1][j][i]*( *(force[0]  ));
        *( force[0]  ) = -lattice->param.G*psi[0][j][i]*( psi_temp     );

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

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

      force[0] += ( sizeof( struct force_struct)/8);
      force[1] += ( sizeof( struct force_struct)/8);

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

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

} /* void compute_fluid_fluid_force( lattice_ptr lattice) */
#else
void compute_fluid_fluid_force( lattice_ptr lattice)
{
  double ***psi; //psi[ NUM_FLUID_COMPONENTS][LX][LY];
  double psi_x[2];
  double psi_y[2];

  double *rho;

  double *force[2];

  int    i,  j, 
         in, jn, 
         ip, jp;

  int    n, ni, nj;

  int    subs;

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

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

  // Initialize psi.
  for( subs=0; subs<NUM_FLUID_COMPONENTS; subs++)
  {
    rho = &( lattice->macro_vars[subs][0].rho);

    if( !lattice->periodic_x[subs] || !lattice->periodic_y[subs])
    {
      printf("%s %d >> Need to be fully periodic to use "
          "NON_LOCAL_FORCES for now.  Exiting!\n"
          __FILE__, __LINE__);
      exit(1);
    }

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

        rho+=3;

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

//printf("sizeof( struct force_struct) = %d\n", sizeof( struct force_struct));
//printf("NumNodes*sizeof( struct force_struct) = %d\n", 
//    lattice->NumNodes*sizeof( struct force_struct));


    for( j=0; j<nj; 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);

        if( !( lattice->bc[0][ j*ni+i].bc_type & BC_SOLID_NODE))
        {

  for( subs=0; subs<NUM_FLUID_COMPONENTS; subs++)
  {
    force[subs] = lattice->force[subs][j*ni+i].force;

        psi_x[subs] = 0.;
        psi_y[subs] = 0.;

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

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

lattice->force[0][j*ni+i].force[0] 
  = -lattice->param.G*psi[0][j][i]*( psi_x[1]);
lattice->force[0][j*ni+i].force[1] 
  = -lattice->param.G*psi[0][j][i]*( psi_y[1]);
lattice->force[1][j*ni+i].force[0] 
  = -lattice->param.G*psi[1][j][i]*( psi_x[0]);
lattice->force[1][j*ni+i].force[1] 
  = -lattice->param.G*psi[1][j][i]*( psi_y[0]);
      }
      else
      {
lattice->force[0][j*ni+i].force[0] = 0.;
lattice->force[0][j*ni+i].force[1] = 0.;
lattice->force[1][j*ni+i].force[0] = 0.;
lattice->force[1][j*ni+i].force[1] = 0.;
      }

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

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

} /* void compute_fluid_fluid_force( lattice_ptr lattice) */
#endif

//void compute_phase_force( lattice_ptr lattice)
//##############################################################################
//
// C O M P U T E   P H A S E   F O R C E
//
#if ZHANG_AND_CHEN_ENERGY_TRANSPORT
//
// Based on Zhang & Chen, PRE 67, 0066711 (2003)
//
// Supposed to give thermodynamic consistency unlike old Shan & Chen method.
// And supports general equation of state P = P(rho,T).
// Utilizes the Inamuro component for evolution of the energy transport
// equation.  Employs modified compute_phase_force routine to compute
// body force term representing non-local interaction potential U among
// particles.
//
void compute_phase_force( lattice_ptr lattice, int subs)
{
  double **U;

  double *rho;
  double *T;

  double *force;

  double a, b, R;

  int    i,  j, 
         in, jn, 
         ip, jp;

  int    n, ni, nj;

  int    sj, ej;

  double y;

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

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

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

  if( !lattice->periodic_x[0] || !lattice->periodic_y[0])
  {
    printf("%s %d >> Need to be fully periodic to use "
        "ZHANG_AND_CHEN_ENERGY_TRANSPORT for now.  Exiting!\n"
        __FILE__, __LINE__);
    exit(1);

  } /* if( !lattice->periodic_x[0] || !lattice->periodic_y[0]) */

  sj = 0;
  ej = nj-1;

  //############################################################################
  //
  // Set U = P(rho,T) - rho*T0
  //
  //       = R*T*( rho/(1-rho*b)) - a*rho^2
  //
  // Maxwell ==> rho_l = 10.8657, rho_v = 4.98648
  //
  a  = 3.592;
  b  = 0.04267;
  R  = 0.082057;
  rho = &( lattice->macro_vars[/*subs*/0][0].rho);
  T   = &( lattice->macro_vars[/*subs*/1][0].rho);
  for( j=sj; j<=ej; j++)
  {
    n = j*ni;

    for( i=0; i<ni; i++, n++)
    {
      if( !( lattice->bc[0][n].bc_type & BC_SOLID_NODE))
      {
#if 0
        // Van der Walls
        U[j][i] = R*(*T)*( (*rho)/( 1. - (*rho)*b))
                - a*(*rho)*(*rho);
#else
        // Carnahan-Starling
        y = (*rho)*b/4.;
        U[j][i] = R*(*T)*(*rho)*( ( 1. + y + y*y - y*y*y) 
                / ( (1.-y)*(1.-y)*(1.-y)))
                - a*(*rho)*(*rho);
//printf("%s (%d) >> U[%d][%d](rho=%f,T=%f) = %20.17f\n", 
//    __FILE__, __LINE__, j, i, *rho, *T, U[j][i]);
#endif

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

      rho+=3;
      T  +=3;

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

//printf("sizeof( struct force_struct) = %d\n", sizeof( struct force_struct));
//printf("NumNodes*sizeof( struct force_struct) = %d\n", 
//    lattice->NumNodes*sizeof( struct force_struct));

  //############################################################################
  //
  // Compute F(x,t) = -\sum_i ( D/(b*c_i^2)) e_i U( x+e_i,t)
  //

  sj = 0;
  ej = nj-1;
  force = ( lattice->force[/*subs*/0][0].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[0][ j*ni+i].bc_type & BC_SOLID_NODE))
      {
/* 1 */ if( !( lattice->bc[0][ j *ni+ip].bc_type & BC_SOLID_NODE)) {
          *( force  ) += 1.*vx[1]*U[j ][ip];
          *( force+1) += 1.*vy[1]*U[j ][ip]; }
/* 2 */ if( !( lattice->bc[0][ jp*ni+i ].bc_type & BC_SOLID_NODE)) {
          *( force  ) += 1.*vx[2]*U[jp][i ];
          *( force+1) += 1.*vy[2]*U[jp][i ]; }
/* 3 */ if( !( lattice->bc[0][ j *ni+in].bc_type & BC_SOLID_NODE)) {
          *( force  ) += 1.*vx[3]*U[j ][in];
          *( force+1) += 1.*vy[3]*U[j ][in]; }
/* 4 */ if( !( lattice->bc[0][ jn*ni+i ].bc_type & BC_SOLID_NODE)) {
          *( force  ) += 1.*vx[4]*U[jn][i ];
          *( force+1) += 1.*vy[4]*U[jn][i ]; }
/* 5 */ if( !( lattice->bc[0][ jp*ni+ip].bc_type & BC_SOLID_NODE)) {
          *( force  ) +=    vx[5]*U[jp][ip];
          *( force+1) +=    vy[5]*U[jp][ip]; }
/* 6 */ if( !( lattice->bc[0][ jp*ni+in].bc_type & BC_SOLID_NODE)) {
          *( force  ) +=    vx[6]*U[jp][in];
          *( force+1) +=    vy[6]*U[jp][in]; }
/* 7 */ if( !( lattice->bc[0][ jn*ni+in].bc_type & BC_SOLID_NODE)) {
          *( force  ) +=    vx[7]*U[jn][in];
          *( force+1) +=    vy[7]*U[jn][in]; }
/* 8 */ if( !( lattice->bc[0][ jn*ni+ip].bc_type & BC_SOLID_NODE)) {
          *( force  ) +=    vx[8]*U[jn][ip];
          *( force+1) +=    vy[8]*U[jn][ip]; }

          *( force  ) = -(1./8.)*( *(force  ));
          *( force+1) = -(1./8.)*( *(force+1));
//printf("%s (%d) >> Fx = %20.17f\n", __FILE__, __LINE__, *(force  ));
//printf("%s (%d) >> Fy = %20.17f\n", __FILE__, __LINE__, *(force+1));

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

      else
      {
          *( force  ) = 0.;
          *( force+1) = 0.;
//printf("%s (%d) >> Fx = ZERO\n", __FILE__, __LINE__);
//printf("%s (%d) >> Fy = ZERO\n", __FILE__, __LINE__);
      }

      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( U[j]);
  }
  free( U);

} /* void compute_phase_force( lattice_ptr lattice, int subs) */

#else /* !( ZHANG_AND_CHEN_ENERGY_TRANSPORT) */
void compute_phase_force( lattice_ptr lattice, int subs)
{
  double **psi; //psi[ NUM_FLUID_COMPONENTS][LX][LY];
  double psi_temp;

  double *rho;

  double *force;

  int    i,  j, 
         in, jn, 
         ip, jp;

  int    n, ni, nj;

  int    sj, ej;

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

  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++, n++)
    {
      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));