compute.c

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

//##############################################################################
//
// compute.c
//
//  - Routines for computing on the lattice:
//
//    - compute_rho_and_u
//    - compute_feq
//    - compute_big_u
//    - compute_gforce
//    - compute_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.
//
void compute_macro_vars( struct lattice_struct *lattice)
{
  int n, a;

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

  double *rho, *u_x, *u_y;
  double *ftemp;
  bc_ptr bc;

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

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

    if( 1 || !( bc[n].bc_type & BC_SOLID_NODE))
    {
      for( a=0; a<9; a++)
      {
        (*rho) += (*ftemp);
//printf("RHO: n=%d, a=%d, *ftemp = %f, *rho = %f\n", n, a, *ftemp, *rho);
        (*u_x) += vx[a]*(*ftemp);
        (*u_y) += vy[a]*(*ftemp);
        ftemp++;

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

//printf("compute_macro_vars() -- ( i, j) = ( %d, %d), *rho = %f\n",
//  lattice->node[n].i,
//  lattice->node[n].j,
//  *rho
//  );

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

      if( fabs(*u_x) > 0)
      {
        *u_x = *u_x / *rho;
      }
      else
      {
        *u_x = 0.;
      }

      if( fabs(*u_y) > 0)
      {
        *u_y = *u_y / *rho;
      }
      else
      {
        *u_y = 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+=3;
    u_x+=3;
    u_y+=3;
    ftemp+=18;

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

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

// void compute_feq( struct lattice_struct *lattice)
//##############################################################################
//
// C O M P U T E   F E Q 
//
//  - Compute equilibriam 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;
  double *feq;
  bc_ptr bc;

//compute_force(  lattice);
//compute_gforce( lattice);
//compute_pforce( lattice);

#define BIG_U_X( u) \
        (u) \
      + lattice->param.gforce[0]

#define BIG_U_Y( u) \
        (u) \
      + lattice->param.gforce[1]

  c  = 1.;

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

  macro_var = lattice->macro_vars[0].rho;
  feq       = lattice->pdf[0].feq;
  bc        = lattice->bc;

  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( n==3) printf("compute_feq() -- n=%d: u_x = %f\n", n, *macro_var);
      ux = BIG_U_X( *macro_var++);
//if( n==3) printf("compute_feq() -- n=%d: ux = %f\n", n, ux);
//if( n==3) printf("compute_feq() -- n=%d: u_y = %f\n", n, *macro_var);
      uy = BIG_U_Y( *macro_var++);
//if( n==3) printf("compute_feq() -- n=%d: uy = %f\n", n, uy);
     }
     else
     {
      ux = *macro_var++;
      uy = *macro_var++;
     }

      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 * (1.0 - f3*usq);
//printf("FEQ: n=%d, *feq = %f\n", n, *(feq-1));

      *feq++ = rt1 * (c + f1*ux + f2*uxsq - f3*usq);
//printf("FEQ: n=%d, *feq = %f\n", n, *(feq-1));
      *feq++ = rt1 * (c + f1*uy + f2*uysq - f3*usq);
//printf("FEQ: n=%d, *feq = %f\n", n, *(feq-1));
      *feq++ = rt1 * (c - f1*ux + f2*uxsq - f3*usq);
//printf("FEQ: n=%d, *feq = %f\n", n, *(feq-1));
      *feq++ = rt1 * (c - f1*uy + f2*uysq - f3*usq);
//printf("FEQ: n=%d, *feq = %f\n", n, *(feq-1));

      *feq++ = rt2*(c+f1*( ux+uy)+f2*( ux+uy)*( ux+uy)-f3*usq);
//printf("FEQ: n=%d, *feq = %f\n", n, *(feq-1));
      *feq++ = rt2*(c+f1*(-ux+uy)+f2*(-ux+uy)*(-ux+uy)-f3*usq);
//printf("FEQ: n=%d, *feq = %f\n", n, *(feq-1));
      *feq++ = rt2*(c+f1*(-ux-uy)+f2*(-ux-uy)*(-ux-uy)-f3*usq);
//printf("FEQ: n=%d, *feq = %f\n", n, *(feq-1));
      *feq++ = rt2*(c+f1*( ux-uy)+f2*( ux-uy)*( ux-uy)-f3*usq);
//printf("FEQ: n=%d, *feq = %f\n", n, *(feq-1));

      feq+=18;
    }
    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;
    }

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

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

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

//void compute_min_rho( lattice_ptr lattice, double *min_rho)
void compute_min_rho( lattice_ptr lattice, double *min_rho)
{
  int n;
  compute_max_rho( lattice, min_rho);

  for( n=0; n<lattice->NumNodes; n++)
  {
    if( !( lattice->bc[n].bc_type & BC_SOLID_NODE))
    {
      if( fabs( lattice->macro_vars[n].rho[0]) < *min_rho)
      {
        *min_rho = fabs( lattice->macro_vars[n].rho[0]);
      }
    }
  }
} /* void compute_min_rho( lattice_ptr lattice, double *min_rho) */

//void compute_max_u( lattice_ptr lattice, double *max_u)
void compute_max_u( lattice_ptr lattice, double *max_u)
{
  int n;
  *max_u = 0.;
  *(max_u+1) = 0.;
  for( n=0; n<lattice->NumNodes; n++)
  {
    if( !( lattice->bc[n].bc_type & BC_SOLID_NODE))
    {
      if( fabs( lattice->macro_vars[n].u[0]) > *(max_u))
      {
        *max_u = fabs( lattice->macro_vars[n].u[0]);
      }
      if( fabs( lattice->macro_vars[n].u[1]) > *(max_u+1))
      {
        *(max_u+1) = fabs( lattice->macro_vars[n].u[1]);
      }
    }
  }
} /* void compute_max_u( lattice_ptr lattice, double *max_u) */

//void compute_min_u( lattice_ptr lattice, double *min_u)
void compute_min_u( lattice_ptr lattice, double *min_u)
{
  int n;
  compute_max_u( lattice, min_u);
  compute_max_u( lattice, min_u+1);
  for( n=0; n<lattice->NumNodes; n++)
  {
    if( !( lattice->bc[n].bc_type & BC_SOLID_NODE))
    {
      if( fabs( lattice->macro_vars[n].u[0]) < *(min_u))
      {
        *min_u = fabs( lattice->macro_vars[n].u[0]);
      }
      if( fabs( lattice->macro_vars[n].u[1]) < *(min_u+1))
      {
        *(min_u+1) = fabs( lattice->macro_vars[n].u[1]);
      }
    }
  }
} /* void compute_min_u( lattice_ptr lattice, double *min_u) */

//void compute_ave_rho( lattice_ptr lattice, double *ave_rho)
void compute_ave_rho( lattice_ptr lattice, double *ave_rho)
{
  int n, nn;
  *ave_rho = 0.;
  nn = 0;
  for( n=0; n<lattice->NumNodes; n++)
  {
    if( !( lattice->bc[n].bc_type & BC_SOLID_NODE))
    {
      *ave_rho += fabs( lattice->macro_vars[n].rho[0]);
      nn++;
    }
  }
  *ave_rho = (*ave_rho) / nn;

} /* void compute_ave_rho( lattice_ptr lattice, double *ave_rho) */

//void compute_ave_u( lattice_ptr lattice, double *ave_u)
void compute_ave_u( lattice_ptr lattice, double *ave_u)
{
  int n, nn;
  *ave_u = 0.;
  *(ave_u+1) = 0.;
  nn = 0;
  for( n=0; n<lattice->NumNodes; n++)
  {
    if( !( lattice->bc[n].bc_type & BC_SOLID_NODE))
    {
      *ave_u += fabs( lattice->macro_vars[n].u[0]);
      *(ave_u+1) += fabs( lattice->macro_vars[n].u[1]);
      nn++;
    }
  }
  *ave_u     = (*ave_u)/nn;
  *(ave_u+1) = (*(ave_u+1))/nn;

} /* void compute_ave_u( lattice_ptr lattice, double *ave_u) */

//void compute_vorticity( lattice_ptr lattice, int i, int j, int n)
void compute_vorticity( lattice_ptr lattice, int i, int j, int n, double *vor)
{
  double duyx, duxy;

  // NOTE: Assuming dx=1 .  TODO: Generalize dx?

  // Derivative of uy wrt x.
  if( lattice->bc[lattice->node[n].nn[1]].bc_type == 0)
  {
    // Forward difference.
    duyx = lattice->macro_vars[lattice->node[n].nn[1]].u[1]
        - lattice->macro_vars[n].u[1];
  }
  else if( lattice->bc[lattice->node[n].nn[3]].bc_type == 0)
  {
    // Backward difference.
    duyx = lattice->macro_vars[n].u[1]
        - lattice->macro_vars[lattice->node[n].nn[3]].u[1];
  }
  else
  {
    duyx = 0.;
  }

  // Derivative of ux wrt y.
  if( lattice->bc[lattice->node[n].nn[2]].bc_type == 0)
  {
    // Forward difference.
    duxy = lattice->macro_vars[lattice->node[n].nn[2]].u[0]
        - lattice->macro_vars[n].u[0];
  }
  else if( lattice->bc[lattice->node[n].nn[4]].bc_type == 0)
  {
    // Backward difference.
    duxy = lattice->macro_vars[n].u[0]
        - lattice->macro_vars[lattice->node[n].nn[4]].u[0];
  }
  else
  {
    duxy = 0.;
  }

  *vor = duyx - duxy;

} /* void compute_vorticity( lattice_ptr lattice, int i, int j, int n) */

//void compute_max_vor(lattice_ptr lattice,double *max_vor_p,double *max_vor_n)
void compute_max_vor( lattice_ptr lattice, double *max_vor_p, double *max_vor_n)
{
  int n;
  double vor;
  int nnz;

  *max_vor_p = 0.;
  *max_vor_n = 0.;
  nnz = 0;

  for( n=0; n<=lattice->NumNodes; n++)
  {
    if( lattice->bc[n].bc_type == 0)
    {
      compute_vorticity( lattice,
                         lattice->node[n].i,
                         lattice->node[n].j,
                         lattice->node[n].n,
                         &vor );
      if( vor != 0.) { nnz++;}

      if( vor > *max_vor_p)
      {
        *max_vor_p = vor;

      } /* if( vor > *max_vor_p) */

      else if( vor < *max_vor_n)
      {
        *max_vor_n = vor;

      } /* if( vor > *max_vor_p) */

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

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

  printf("compute_max_vor() -- nnz = %d.  nnz/NumNodes = %f\n", 
    nnz, (double)nnz/(double)lattice->NumNodes);

} /* void compute_max_vor( lattice_ptr lattice, double *max_vor) */

//void compute_ave_vor(lattice_ptr lattice,double *ave_vor_p,double *ave_vor_n)
void compute_ave_vor( lattice_ptr lattice, double *ave_vor_p, double *ave_vor_n)
{
  int n;
  double vor;
  int nnz;
  int num_p;
  int num_n;

  *ave_vor_p = 0.;
  *ave_vor_n = 0.;
  nnz = 0;
  num_p = 0;
  num_n = 0;

  for( n=0; n<=lattice->NumNodes; n++)
  {
    if( lattice->bc[n].bc_type == 0)
    {
      compute_vorticity( lattice,
                         lattice->node[n].i,
                         lattice->node[n].j,
                         lattice->node[n].n,
                         &vor );
      if( vor != 0.) { nnz++;}

      if( vor > *ave_vor_p)
      {
        *ave_vor_p += vor;
        num_p++;

      } /* if( vor > *ave_vor_p) */

      else if( vor < *ave_vor_n)
      {
        *ave_vor_n += vor;
        num_n++;

      } /* if( vor > *ave_vor_p) */

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

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

  *ave_vor_p /= num_p;
  *ave_vor_n /= num_n;

  printf("compute_ave_vor() -- nnz = %d.  nnz/NumNodes = %f\n", 
    nnz, (double)nnz/(double)lattice->NumNodes);

} /* void compute_ave_vor( lattice_ptr lattice, double *ave_vor) */