compute.c

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

  bc_ptr bc;
  int    subs;

#if NON_LOCAL_FORCES

#if ZHANG_AND_CHEN_ENERGY_TRANSPORT

#define BIG_U_X( u_) \
        (u_) \
        + lattice->param.tau[subs]  \
          * lattice->param.gforce[subs][0]

#define BIG_U_Y( u_) \
        (u_) \
        + lattice->param.tau[subs]  \
          * lattice->param.gforce[subs][1]

#define BIG_U_X_BUOY( u_, rho1_, rho2_) 1.
#define BIG_U_Y_BUOY( u_, rho1_, rho2_) 1.

#else /* !( ZHANG_AND_CHEN_ENERGY_TRANSPORT) */

#define BIG_U_X( u_, rho_) \
        (u_) \
        + lattice->param.tau[subs]  \
          * lattice->force[subs][n].force[0]/(rho_) \
        + lattice->param.tau[subs]  \
          * lattice->force[subs][n].sforce[0]/(rho_) \
        + lattice->param.tau[subs]  \
          * lattice->param.gforce[subs][0]

#define BIG_U_Y( u_, rho_) \
        (u_) \
        + lattice->param.tau[subs]  \
          * lattice->force[subs][n].force[1]/(rho_) \
        + lattice->param.tau[subs]  \
          * lattice->force[subs][n].sforce[1]/(rho_) \
        + lattice->param.tau[subs]  \
          * lattice->param.gforce[subs][1]

#endif /* ZHANG_AND_CHEN_ENERGY_TRANSPORT */

  if( NUM_FLUID_COMPONENTS == 2)
  {
#if ZHANG_AND_CHEN_ENERGY_TRANSPORT
    compute_phase_force( lattice, 0);
#else /* !( ZHANG_AND_CHEN_ENERGY_TRANSPORT) */
    if( lattice->param.G != 0.)
    {
      compute_fluid_fluid_force( lattice);
    }
#endif /* ZHANG_AND_CHEN_ENERGY_TRANSPORT */
  }
  else if( NUM_FLUID_COMPONENTS == 1)
  {
    if( lattice->param.G != 0.)
    {
      compute_phase_force( lattice, 0);
    }
    if(    lattice->param.Gads[0] != 0.)
    {
      compute_single_fluid_solid_force( lattice, 0);
    }
  }
  else
  {
    printf(
      "compute_feq() -- "
      "Unhandled case NUM_FLUID_COMPONENTS = %d .  "
      "Exiting!\n", NUM_FLUID_COMPONENTS);
    exit(1);
  }
  //dump_forces( lattice);
  //force2bmp( lattice);
#else /* !( NON_LOCAL_FORCES) */

#if INAMURO_SIGMA_COMPONENT

#define BIG_U_X_BAK( u_, rho1_, rho2_) \
        (u_) \
      + lattice->param.tau[subs]  \
        * lattice->param.gforce[subs][0]*((rho1_)+(rho2_))\
                                        /lattice->param.rho_A[0]

#define BIG_U_X( u_) \
        (u_) \
      + lattice->param.tau[subs]  \
        * lattice->param.gforce[subs][0]\
      + .00000*( (double)rand()/(double)RAND_MAX - .5)

#define BIG_U_X_BUOY( u_, rho1_, rho2_) \
        (u_) \
      + lattice->param.tau[subs]  \
        * lattice->param.gforce[subs][0] \
        *((rho1_)+(get_buoyancy_sign(lattice))*(rho2_)) / (rho1_)

#define BIG_U_Y_BAK( u_, rho1_, rho2_) \
        (u_) \
      + lattice->param.tau[subs]  \
        * lattice->param.gforce[subs][1]*((rho1_)+(rho2_))\
                                        /lattice->param.rho_A[0]

#define BIG_U_Y( u_) \
        (u_) \
      + lattice->param.tau[subs]  \
        * lattice->param.gforce[subs][1]\
      + .00000*( (double)rand()/(double)RAND_MAX - .5)

#define BIG_U_Y_BUOY( u_, rho1_, rho2_) \
        (u_) \
      + lattice->param.tau[subs]  \
        * lattice->param.gforce[subs][1] \
        *((rho1_)+(get_buoyancy_sign(lattice))*(rho2_)) / (rho1_)

#else /* !( INAMURO_SIGMA_COMPONENT) */

#define BIG_U_X( u_, rho_) \
        (u_) \
      + lattice->param.tau[subs]  \
        * lattice->param.gforce[subs][0]

#define BIG_U_Y( u_, rho_) \
        (u_) \
      + lattice->param.tau[subs]  \
        * lattice->param.gforce[subs][1]

#endif /* INAMURO_SIGMA_COMPONENT */

#endif /* NON_LOCAL_FORCES */

//printf("%s %d >> rand() = %f\n", __FILE__, __LINE__, .00001*(double)rand()/(double)RAND_MAX);

 f1 = 3.;
 f2 = 9./2.;
 f3 = 3./2.;

 for( subs=0; subs<(NUM_FLUID_COMPONENTS)-(INAMURO_SIGMA_COMPONENT); subs++)
 {

  macro_var = &( lattice->macro_vars[subs][0].rho);
#if INAMURO_SIGMA_COMPONENT
  u0         =    lattice->macro_vars[0][0].u;
  u1         =    lattice->macro_vars[1][0].u;
#endif /* INAMURO_SIGMA_COMPONENT */
  feq       =    lattice->pdf[subs][0].feq;
#if STORE_UEQ
  ueq       =    lattice->ueq[0].u;
#endif /* STORE_UEQ */
  bc        =    lattice->bc[subs];

  for( n=0; n<lattice->NumNodes; n++)
  {
    if( 1 || !( bc[n].bc_type & BC_SOLID_NODE))
    {
      rt0 = *macro_var++; // Preserve raw density until after BIG_U.

     if( 1 || !( bc[n].bc_type & BC_SOLID_NODE))
     {
#if STORE_UEQ
      ux = BIG_U_X( *ueq, ((rt0!=0)?(rt0):(1))); *ueq++;// = ux;
      uy = BIG_U_Y( *ueq, ((rt0!=0)?(rt0):(1))); *ueq++;// = uy;
      macro_var+=2;
#else /* !( STORE_UEQ) */

#if INAMURO_SIGMA_COMPONENT
      if( get_buoyancy_flag( lattice))
      {
        ux = BIG_U_X_BUOY( *macro_var, rt0, lattice->macro_vars[1][n].rho); 
      }
      else
      {
//#if ZHANG_AND_CHEN_ENERGY_TRANSPORT
//        ux = BIG_U_X( *macro_var, rt0);
//#else /* !( ZHANG_AND_CHEN_ENERGY_TRANSPORT) */
        ux = BIG_U_X( *macro_var); 
//#endif /* ZHANG_AND_CHEN_ENERGY_TRANSPORT */
      }
      macro_var++;
      *u1++ = *u0++;
#else /* !( INAMURO_SIGMA_COMPONENT) */
      ux = BIG_U_X( *macro_var, rt0); macro_var++;
#endif /* INAMURO_SIGMA_COMPONENT */
#if INAMURO_SIGMA_COMPONENT
      if( get_buoyancy_flag( lattice))
      {
        uy = BIG_U_Y_BUOY( *macro_var, rt0, lattice->macro_vars[1][n].rho); 
      }
      else
      {
//#if ZHANG_AND_CHEN_ENERGY_TRANSPORT
//        uy = BIG_U_Y( *macro_var, rt0); macro_var++;
//#else /* !( ZHANG_AND_CHEN_ENERGY_TRANSPORT) */
        uy = BIG_U_Y( *macro_var); 
//#endif /* ZHANG_AND_CHEN_ENERGY_TRANSPORT */
      }
      macro_var++;
      *u1++ = *u0++;
#else /* !( INAMURO_SIGMA_COMPONENT) */
      uy = BIG_U_Y( *macro_var, rt0); macro_var++;
#endif /* INAMURO_SIGMA_COMPONENT */

#endif /* STORE_UEQ */
     }
     else
     {
#if STORE_UEQ
      ux = *ueq++;
      uy = *ueq++;
      macro_var+=2;
#else /* !( STORE_UEQ) */

      ux = *macro_var++;
#if INAMURO_SIGMA_COMPONENT
      *u1++ = *u0++;
#endif /* INAMURO_SIGMA_COMPONENT */
      uy = *macro_var++;
#if INAMURO_SIGMA_COMPONENT
      *u1++ = *u0++;
#endif /* INAMURO_SIGMA_COMPONENT */

#endif /* STORE_UEQ */
     }

  if( lattice->param.incompressible)
  {
    c  = rt0; // rt0 == rho
    rt1 = 1./9.; 
    rt2 = 1./36.;
    rt0 = 4./9.; // Overwrite rt0.
  }
  else // compressible
  {
    c  = 1.;
    rt1 = rt0/(9.); // rt0 == rho
    rt2 = rt0/(36.);// rt0 == rho
    rt0 *= (4./9.); // Update rt0 now that raw density is no longer needed.
  }


      uxsq = ux*ux;
      uysq = uy*uy;

      usq = uxsq + uysq;

      *feq++ = rt0 * (c - f3*usq);
      *feq++ = rt1 * (c + f1*ux + f2*uxsq - f3*usq);
      *feq++ = rt1 * (c + f1*uy + f2*uysq - f3*usq);
      *feq++ = rt1 * (c - f1*ux + f2*uxsq - f3*usq);
      *feq++ = rt1 * (c - f1*uy + f2*uysq - f3*usq);
      *feq++ = rt2*(c+f1*( ux+uy)+f2*( ux+uy)*( ux+uy)-f3*usq);
      *feq++ = rt2*(c+f1*(-ux+uy)+f2*(-ux+uy)*(-ux+uy)-f3*usq);
      *feq++ = rt2*(c+f1*(-ux-uy)+f2*(-ux-uy)*(-ux-uy)-f3*usq);
      *feq++ = rt2*(c+f1*( ux-uy)+f2*( ux-uy)*( ux-uy)-f3*usq);

      feq+=18;
#if INAMURO_SIGMA_COMPONENT
      u0++;
      u1++;
#endif /* INAMURO_SIGMA_COMPONENT */
    }
    else
    {
#if 1
      *feq++ = 0.;
      *feq++ = 0.;
      *feq++ = 0.;
      *feq++ = 0.;
      *feq++ = 0.;
      *feq++ = 0.;
      *feq++ = 0.;
      *feq++ = 0.;
      *feq++ = 0.;
      feq+=18;
#else
      feq+=27;
#endif
      macro_var+=3;
#if STORE_UEQ
      ueq+=2;
#endif /* STORE_UEQ */
    }

  } /* for( n=0; n<lattice->NumNodes; n++) */

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

#if INAMURO_SIGMA_COMPONENT
  subs=1;

  c  = 1.;

#if ZHANG_AND_CHEN_ENERGY_TRANSPORT
 f1 = 3.;
 f2 = 9./2.;
 f3 = 3./2.;
#else /* !( ZHANG_AND_CHEN_ENERGY_TRANSPORT) */
  f1 = 3.;
  f2 = 0.;
  f3 = 0.;
#endif /* ZHANG_AND_CHEN_ENERGY_TRANSPORT */

  rho       = &( lattice->macro_vars[subs][0].rho);
  u         =    lattice->macro_vars[subs][0].u;
  feq       =    lattice->pdf[subs][0].feq;
  bc        =    lattice->bc[subs];

  for( n=0; n<lattice->NumNodes; n++)
  {
    if( 1 || !( bc[n].bc_type & BC_SOLID_NODE))
    {
      rt0 = *rho; // Preserve raw density until after BIG_U.
      rho+=3; 

      ux = *u++;
      uy = *u++;

      if( lattice->param.simple_diffusion)
      {
        rt1 = rt0/(8.); // rt0 == rho
        rt0 *= (1./2.); // Update rt0 now that raw density is no longer needed.

        *feq++ = rt0 * ( 1.             );

        *feq++ = rt1 * ( c + f1*ux      );
        *feq++ = rt1 * ( c + f1*uy      );
        *feq++ = rt1 * ( c - f1*ux      );
        *feq++ = rt1 * ( c - f1*uy      );

        *feq++ = 0.; //rt2 * ( c + f1*( ux+uy));
        *feq++ = 0.; //rt2 * ( c + f1*(-ux+uy));
        *feq++ = 0.; //rt2 * ( c + f1*(-ux-uy));
        *feq++ = 0.; //rt2 * ( c + f1*( ux-uy));
      }
      else
      {
        rt1 = rt0/(9.); // rt0 == rho
        rt2 = rt0/(36.);// rt0 == rho
        rt0 *= (4./9.); // Update rt0 now that raw density is no longer needed.

#if ZHANG_AND_CHEN_ENERGY_TRANSPORT
      uxsq = ux*ux;
      uysq = uy*uy;

      usq = uxsq + uysq;

      *feq++ = rt0 * (c - f3*usq);
      *feq++ = rt1 * (c + f1*ux + f2*uxsq - f3*usq);
      *feq++ = rt1 * (c + f1*uy + f2*uysq - f3*usq);
      *feq++ = rt1 * (c - f1*ux + f2*uxsq - f3*usq);
      *feq++ = rt1 * (c - f1*uy + f2*uysq - f3*usq);
      *feq++ = rt2*(c+f1*( ux+uy)+f2*( ux+uy)*( ux+uy)-f3*usq);
      *feq++ = rt2*(c+f1*(-ux+uy)+f2*(-ux+uy)*(-ux+uy)-f3*usq);
      *feq++ = rt2*(c+f1*(-ux-uy)+f2*(-ux-uy)*(-ux-uy)-f3*usq);
      *feq++ = rt2*(c+f1*( ux-uy)+f2*( ux-uy)*( ux-uy)-f3*usq);
#else /* !( ZHANG_AND_CHEN_ENERGY_TRANSPORT) */
        *feq++ = rt0 * ( 1.             );

        *feq++ = rt1 * ( c + f1*ux      );
        *feq++ = rt1 * ( c + f1*uy      );
        *feq++ = rt1 * ( c - f1*ux      );
        *feq++ = rt1 * ( c - f1*uy      );

        *feq++ = rt2 * ( c + f1*( ux+uy));
        *feq++ = rt2 * ( c + f1*(-ux+uy));
        *feq++ = rt2 * ( c + f1*(-ux-uy));
        *feq++ = rt2 * ( c + f1*( ux-uy));
#endif /* ZHANG_AND_CHEN_ENERGY_TRANSPORT */
      }

      feq+=18;
      u++;
    }
    else
    {
#if 1
      *feq++ = 0.;
      *feq++ = 0.;
      *feq++ = 0.;
      *feq++ = 0.;
      *feq++ = 0.;
      *feq++ = 0.;
      *feq++ = 0.;
      *feq++ = 0.;
      *feq++ = 0.;
      feq+=18;
#else
      feq+=27;
#endif
      rho+=3;
      u+=3;
    }

  } /* for( n=0; n<lattice->NumNodes; n++) */

#endif /* INAMURO_SIGMA_COMPONENT */

 //dump_pdf( lattice, lattice->time);

} /* void compute_feq( struct lattice_struct *lattice) */

#if NON_LOCAL_FORCES
// void compute_fluid_fluid_force( lattice_ptr lattice)
//##############################################################################
//
// C O M P U T E   F L U I D   F L U I D   F O R C E
//
#if 1
void compute_fluid_fluid_force( lattice_ptr lattice)
{
  double ***psi; //psi[ NUM_FLUID_COMPONENTS][LX][LY];
  double psi_temp;

  double *rho;

  double *force[2];

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

  int    n, ni, nj;

  int    sj, ej;

  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( ni*sizeof(double));
    }

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

  // Initialize psi.
  for( subs=0; subs<NUM_FLUID_COMPONENTS; subs++)
  {
    for( j=0; j<nj; j++)
    {
      for( i=0; i<ni; i++, n++)
      {
        psi[subs][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);
    }
    if( !lattice->periodic_x[subs])
    {
      printf("%s %d >> Need to be periodic in x-direction to use "
          "NON_LOCAL_FORCES for now.  Exiting!\n",
          __FILE__, __LINE__);
      exit(1);
    }

    for( j=sj; j<=ej; 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;
        }

        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( subs=0; subs<NUM_FLUID_COMPONENTS; subs++)
  {
    if( lattice->periodic_y[subs])
    {
      sj = 0;
      ej = nj-1;
      force[subs] = lattice->force[subs][0].force;
    }