compute.c

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

//##############################################################################
//
// Copyright (C), 2005, Michael Sukop and Danny Thorne
//
// compute.c
//
//  - Routines for computing on the lattice:
//
//    - compute_rho_and_u
//    - compute_feq
//    - compute_big_u
//    - compute_gforce
//    - compute_fluid_fluid_force
//    - etc...
//

// void compute_macro_vars( struct lattice_struct *lattice)
//##############################################################################
//
// C O M P U T E   M A C R O   V A R S 
//
//  - Compute macroscopic variables.
//
#if INAMURO_SIGMA_COMPONENT
void compute_macro_vars( struct lattice_struct *lattice)
{
  int n, a;

  double *rho[ NUM_FLUID_COMPONENTS], 
         *u_x[ NUM_FLUID_COMPONENTS], 
         *u_y[ NUM_FLUID_COMPONENTS];

  double *ueq;

  double *ftemp;

  bc_ptr bc;

  int    subs;

  double c;

  //############################################################################
  //
  // For first component, compute rho and u.
  //
  subs=0;

  rho[subs]   = &( lattice->macro_vars[subs][0].rho);
  u_x[subs]   =    lattice->macro_vars[subs][0].u;
  u_y[subs]   =    lattice->macro_vars[subs][0].u + 1;
  ftemp       =    lattice->pdf[subs][0].ftemp;
  bc          =    lattice->bc[subs];

  for( n=0; n<lattice->NumNodes; n++)
  {
    *rho[subs] = 0.;
    *u_x[subs] = 0.;
    *u_y[subs] = 0.;

    if( 1 || !( bc[n].bc_type & BC_SOLID_NODE))
    {
      for( a=0; a<9; a++)
      {
        (*rho[subs]) += (*ftemp);
        (*u_x[subs]) += vx[a]*(*ftemp);
        (*u_y[subs]) += vy[a]*(*ftemp);
        ftemp++;

      } /* for( a=0; a<9; a++) */

#if PUKE_NEGATIVE_DENSITIES
      if( *rho[subs] < 0.)
      {
        printf("\n");
        printf(
          "compute_macro_vars() -- "
          "Node %d (%d,%d) has negative density %20.17f "
          "at timestep %d. Exiting!\n", 
          n, n%lattice->param.LX,
             n/lattice->param.LX, *rho[subs], 
             lattice->time             );
        printf("\n");
        exit(1);
      }
#endif /* PUKE_NEGATIVE_DENSITIES */

      if( *rho[subs] == 0)
      {
        printf("\n");
        printf("\n");
        printf("%s (%d) -- "
               "ERROR:  rho[subs=%d][j=%d][i=%d] = 0. "
               "at timestep %d.  "
               "Exiting!\n",
               __FILE__,__LINE__,
               subs,
               n/lattice->param.LX, 
               n%lattice->param.LX,
               lattice->time        );
        printf("\n");
        printf("\n");
        exit(1);
      }

    } /* if( !( bc++->bc_type & BC_SOLID_NODE)) */

    else // bc++->bc_type & BC_SOLID_NODE
    {
//printf("RHO: n=%d, Solid node.\n", n);
      ftemp+=9;

    } /* if( !( bc++->bc_type & BC_SOLID_NODE)) else */

    rho[subs]+=3;
    u_x[subs]+=3;
    u_y[subs]+=3;
    ftemp+=18;

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

  //############################################################################
  //
  // For second component, compute just rho.
  //
  subs=1;

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

  for( n=0; n<lattice->NumNodes; n++)
  {
    *rho[subs] = 0.;

    if( 1 || !( bc[n].bc_type & BC_SOLID_NODE))
    {
      for( a=0; a<9; a++)
      {
        (*rho[subs]) += (*ftemp);
        ftemp++;

      } /* for( a=0; a<9; a++) */

#if PUKE_NEGATIVE_CONCENTRATIONS
      if( *rho[subs] < 0.)
      {
        printf("\n");
        printf(
          "compute_macro_vars() -- "
          "Node %d (%d,%d) has negative density %20.17f "
          "at timestep %d. Exiting!\n", 
          n, n%lattice->param.LX,
             n/lattice->param.LX, *rho[subs], 
             lattice->time             );
        printf("\n");
        exit(1);
      }
#endif /* PUKE_NEGATIVE_CONCENTRATIONS */
      //assert( *rho[subs] != 0);

    } /* if( !( bc++->bc_type & BC_SOLID_NODE)) */

    else // bc++->bc_type & BC_SOLID_NODE
    {
//printf("RHO: n=%d, Solid node.\n", n);
      ftemp+=9;

    } /* if( !( bc++->bc_type & BC_SOLID_NODE)) else */

    rho[subs]+=3;
    ftemp+=18;

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


  rho[0]   = &( lattice->macro_vars[0][0].rho);
  u_x[0]   =    lattice->macro_vars[0][0].u;
  u_x[1]   =    lattice->macro_vars[1][0].u;
  u_y[0]   =    lattice->macro_vars[0][0].u + 1;
  u_y[1]   =    lattice->macro_vars[1][0].u + 1;

  for( n=0; n<lattice->NumNodes; n++)
  {
    if( 1 || !( bc[n].bc_type & BC_SOLID_NODE))
    {
      if( *rho[0] == 0)
      {
        printf("\n");
        printf("\n");
        printf("%s (%d) -- "
               "ERROR:  rho[subs=%d][j=%d][i=%d] = 0. "
               "at timestep %d.  "
               "Exiting!\n",
               __FILE__,__LINE__,
               0,
               n/lattice->param.LX, 
               n%lattice->param.LX,
               lattice->time        );
        printf("\n");
        printf("\n");
        exit(1);
      }


      if( lattice->param.incompressible)
      {
        c = 1.;
      }
      else
      {
        c = *rho[0];
      }

      *u_x[0] = *u_x[0] / c;
      *u_x[1] = *u_x[0];// / c;
      *u_y[0] = *u_y[0] / c;
      *u_y[1] = *u_y[0];// / c;

    } /* if( 1 || !( bc[n].bc_type & BC_SOLID_NODE)) */
#if STORE_UEQ
    else
    {
      ueq+=2;
    } /* if( 1 || !( bc[n].bc_type & BC_SOLID_NODE)) else */
#endif /* STORE_UEQ */

    rho[0]+=3;
    u_x[0]+=3;
    u_x[1]+=3;
    u_y[0]+=3;
    u_y[1]+=3;

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

} /* void compute_macro_vars( struct lattice_struct *lattice) */
#else /* !( INAMURO_SIGMA_COMPONENT) */
void compute_macro_vars( struct lattice_struct *lattice)
{
  int n, a;

  double *rho[ NUM_FLUID_COMPONENTS], 
         *u_x[ NUM_FLUID_COMPONENTS], 
         *u_y[ NUM_FLUID_COMPONENTS];

  double ux_sum, uy_sum;

  double *ueq;

  double *ftemp;

  bc_ptr bc;

  int    subs;

  double tau0,
         tau1;

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

  rho[subs]   = &( lattice->macro_vars[subs][0].rho);
  u_x[subs]   =    lattice->macro_vars[subs][0].u;
  u_y[subs]   =    lattice->macro_vars[subs][0].u + 1;
  ftemp       =    lattice->pdf[subs][0].ftemp;
  bc          =    lattice->bc[subs];

  for( n=0; n<lattice->NumNodes; n++)
  {
    *rho[subs] = 0.;
    *u_x[subs] = 0.;
    *u_y[subs] = 0.;

    if( 0 || !( bc[n].bc_type & BC_SOLID_NODE))
    {
      for( a=0; a<9; a++)
      {
        (*rho[subs]) += (*ftemp);
        (*u_x[subs]) += vx[a]*(*ftemp);
        (*u_y[subs]) += vy[a]*(*ftemp);
        ftemp++;

      } /* for( a=0; a<9; a++) */

#if PUKE_NEGATIVE_DENSITIES
      if( *rho[subs] < 0.)
      {
        printf("\n");
        printf(
          "compute_macro_vars() -- "
          "Node %d (%d,%d) has negative density %20.17f "
          "at timestep %d. Exiting!\n", 
          n, n%lattice->param.LX,
             n/lattice->param.LX, *rho[subs], 
             lattice->time             );
        printf("\n");
        exit(1);
      }
#endif /* PUKE_NEGATIVE_DENSITIES */
      //assert( *rho[subs] != 0);

    } /* if( !( bc++->bc_type & BC_SOLID_NODE)) */

    else // bc++->bc_type & BC_SOLID_NODE
    {
//printf("RHO: n=%d, Solid node.\n", n);
      ftemp+=9;

    } /* if( !( bc++->bc_type & BC_SOLID_NODE)) else */

    rho[subs]+=3;
    u_x[subs]+=3;
    u_y[subs]+=3;
    ftemp+=18;

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

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

 for( subs=0; subs<NUM_FLUID_COMPONENTS; subs++)
 {
  rho[subs]   = &( lattice->macro_vars[subs][0].rho);
  u_x[subs]   =    lattice->macro_vars[subs][0].u;
  u_y[subs]   =    lattice->macro_vars[subs][0].u + 1;

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

#if STORE_UEQ
  ueq      =    lattice->ueq[0].u;
#endif /* STORE_UEQ */

 if( NUM_FLUID_COMPONENTS==2)
 {
  tau0 = lattice->param.tau[0];
  tau1 = lattice->param.tau[1];
  for( n=0; n<lattice->NumNodes; n++)
  {
    if( 0 || !( bc[n].bc_type & BC_SOLID_NODE))
    {
      ux_sum =  *u_x[0]/tau0 + *u_x[1]/tau1;
      uy_sum =  *u_y[0]/tau0 + *u_y[1]/tau1;

#if STORE_UEQ
      //assert( *rho[0] + *rho[1] != 0.);
      if( *rho[0] + *rho[1] != 0)
      {
        *ueq++ = ( ux_sum) / ( *rho[0]/tau0 + *rho[1]/tau1);
        *ueq++ = ( uy_sum) / ( *rho[0]/tau0 + *rho[1]/tau1);
      }
      else
      {
        *ueq++ = 0.;
        *ueq++ = 0.;
      }
#endif /* STORE_UEQ */

      //assert( *rho[0] != 0.);
      //assert( *rho[1] != 0.);

      if( ux_sum != 0.) 
      { 
        if( *rho[0] != 0)
        {
          *u_x[0] = ux_sum / *rho[0];
        }
        else
        {
          *u_x[0] = 0.;
        }

        if( *rho[1] != 0)
        {
          *u_x[1] = ux_sum / *rho[1];
        }
        else
        {
          *u_x[1] = 0.;
        }
      }
      else
      {
        *u_x[0] = 0.;
        *u_x[1] = 0.;
      }

      if( uy_sum != 0.) 
      { 
        if( *rho[0] != 0)
        {
          *u_y[0] = uy_sum / *rho[0];
        }
        else
        {
          *u_y[0] = 0.;
        }

        if( *rho[1] != 0)
        {
          *u_y[1] = uy_sum / *rho[1];
        }
        else
        {
          *u_y[1] = 0.;
        }
      }
      else
      {
        *u_y[0] = 0.;
        *u_y[1] = 0.;
      }

    } /* if( 1 || !( bc[n].bc_type & BC_SOLID_NODE)) */
#if STORE_UEQ
    else
    {
      ueq+=2;
    } /* if( 1 || !( bc[n].bc_type & BC_SOLID_NODE)) else */
#endif /* STORE_UEQ */

    rho[0]+=3;
    u_x[0]+=3;
    u_y[0]+=3;

    rho[1]+=3;
    u_x[1]+=3;
    u_y[1]+=3;

  } /* for( n=0; n<lattice->NumNodes; n++) */
 }
 else if( NUM_FLUID_COMPONENTS == 1)
 {
  for( n=0; n<lattice->NumNodes; n++)
  {
    if( 0 || !( bc[n].bc_type & BC_SOLID_NODE))
    {
      //assert( *rho[0] != 0.);

      if( *rho[0] != 0 && *u_x[0] != 0.) 
      { 
        *u_x[0] = *u_x[0] / *rho[0];
      }
      else
      {
        *u_x[0] = 0.;
      }

      if( *rho[0] != 0 && *u_y[0] != 0.) 
      { 
        *u_y[0] = *u_y[0] / *rho[0];
      }
      else
      {
        *u_y[0] = 0.;
      }

    } /* if( 1 || !( bc[n].bc_type & BC_SOLID_NODE)) */

    rho[0]+=3;
    u_x[0]+=3;
    u_y[0]+=3;

  } /* for( n=0; n<lattice->NumNodes; n++) */
 }
 else
 {
  printf(
    "compute_macro_vars() -- "
    "Unhandled case "
    "NUM_FLUID_COMPONENTS = %d . "
    "Exiting!\n",
    NUM_FLUID_COMPONENTS);
  exit(1);
 }

} /* void compute_macro_vars( struct lattice_struct *lattice) */
#endif /* INAMURO_SIGMA_COMPONENT */

// void compute_feq( struct lattice_struct *lattice)
//##############################################################################
//
// C O M P U T E   F E Q 
//
//  - Compute equilibrium distribution function, feq.
//
void compute_feq( struct lattice_struct *lattice)
{
  int n, a;

  double rt0,  rt1,  rt2;
  double f1,   f2,   f3;
  double ux,   uy, 
         uxsq, uysq, usq;
  double c;

  double *macro_var;
#if INAMURO_SIGMA_COMPONENT
  double *rho, *u, *u0, *u1;
#endif /* INAMURO_SIGMA_COMPONENT */
  double *feq;
#if STORE_UEQ
  double *ueq;
#endif /* STORE_UEQ */