compute.c

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

    {
      if(     ( lattice->macro_vars[subs][n].rho) > *max_rho)
      {
        *max_rho =     ( lattice->macro_vars[subs][n].rho);
      }
    }
  }
} /* void compute_max_rho( lattice_ptr lattice, double *max_rho, int subs) */

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

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

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


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

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

//void compute_max_u_all( lattice_ptr lattice, double *max_u, int subs)
void compute_max_u_all( lattice_ptr lattice, double *max_u, int subs)
{
  int n;
  double rho, u, u_x, u_y;
  *(max_u+0) = 0.;
  *(max_u+1) = 0.;
  *(max_u+2) = 0.;
  *(max_u+3) = 0.;
  *(max_u+4) = 0.;
  for( n=0; n<lattice->NumNodes; n++)
  {
    if( !( lattice->bc[subs][n].bc_type & BC_SOLID_NODE))
    {
      rho = lattice->macro_vars[subs][n].rho;
      u_x = lattice->macro_vars[subs][n].u[0];
      u_y = lattice->macro_vars[subs][n].u[1];

      u = sqrt(u_x*u_x+u_y*u_y);
      if( u > *(max_u+0)) { *(max_u+0) = u; }

      if( fabs( u_x) > *(max_u+1)) { *(max_u+1) = fabs( u_x); }
      if( fabs( u_y) > *(max_u+2)) { *(max_u+2) = fabs( u_y); }

      if(     ( u_x) > *(max_u+3)) { *(max_u+3) =     ( u_x); }
      if(     ( u_y) > *(max_u+4)) { *(max_u+4) =     ( u_y); }
    }
  }
} /* void compute_max_u_all( lattice_ptr lattice, double *max_u, int subs) */

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

//void compute_min_u_all( lattice_ptr lattice, double *min_u, int subs)
void compute_min_u_all( lattice_ptr lattice, double *min_u, int subs)
{
  int n;
  double rho, u, u_x, u_y;

  compute_max_u_all( lattice, min_u,   subs);

  for( n=0; n<lattice->NumNodes; n++)
  {
    if( !( lattice->bc[subs][n].bc_type & BC_SOLID_NODE))
    {
      rho = lattice->macro_vars[subs][n].rho;
      u_x = lattice->macro_vars[subs][n].u[0];
      u_y = lattice->macro_vars[subs][n].u[1];

      u = sqrt(u_x*u_x+u_y*u_y);
      if( u < *(min_u+0)) { *(min_u+0) = u; }

      if( fabs( u_x) < *(min_u+1)) { *(min_u+1) = fabs( u_x); }
      if( fabs( u_y) < *(min_u+2)) { *(min_u+2) = fabs( u_y); }

      if(     ( u_x) < *(min_u+3)) { *(min_u+3) =     ( u_x); }
      if(     ( u_y) < *(min_u+4)) { *(min_u+4) =     ( u_y); }

    }
  }
} /* void compute_min_u_all( lattice_ptr lattice, double *min_u) */

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

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

//void compute_ave_u_all( lattice_ptr lattice, double *ave_u, int subs)
void compute_ave_u_all( lattice_ptr lattice, double *ave_u, int subs)
{
  int n, nn;
  double rho, u_x, u_y;
  *(ave_u+0) = 0.;
  *(ave_u+1) = 0.;
  *(ave_u+2) = 0.;
  *(ave_u+3) = 0.;
  *(ave_u+4) = 0.;

  nn = 0;
  for( n=0; n<lattice->NumNodes; n++)
  {
    if( !( lattice->bc[subs][n].bc_type & BC_SOLID_NODE))
    {
      rho = lattice->macro_vars[subs][n].rho;
      u_x = lattice->macro_vars[subs][n].u[0];
      u_y = lattice->macro_vars[subs][n].u[1];

      *(ave_u+0) += sqrt(u_x*u_x+u_y*u_y);

      *(ave_u+1) += fabs(u_x);
      *(ave_u+2) += fabs(u_y);

      *(ave_u+3) += u_x;
      *(ave_u+4) += u_y;

      nn++;
    }
  }
  if( nn != 0) 
  { 
    *(ave_u+0) = (*(ave_u+0))/nn;
    *(ave_u+1) = (*(ave_u+1))/nn;
    *(ave_u+2) = (*(ave_u+2))/nn;
    *(ave_u+3) = (*(ave_u+3))/nn;
    *(ave_u+4) = (*(ave_u+4))/nn;
  }

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


// void compute_flux( lattice_ptr lattice, double *flux, int subs)
void compute_flux( lattice_ptr lattice, double *flux, int subs)
{
  int n, nn;
  double rho, u_x, u_y;
  *(flux+0) = 0.;
  *(flux+1) = 0.;
  *(flux+2) = 0.;
  nn = 0;
  for( n=0; n<lattice->NumNodes; n++)
  {
    rho = lattice->macro_vars[subs][n].rho;
    u_x = lattice->macro_vars[subs][n].u[0];
    u_y = lattice->macro_vars[subs][n].u[1];

    if( !( lattice->bc[subs][n].bc_type & BC_SOLID_NODE))
    {
      *(flux+0) += rho*sqrt(u_x*u_x+u_y*u_y);
      *(flux+1) += rho*u_x;
      *(flux+2) += rho*u_y;
      nn++;
    }
  }
  if( nn != 0) 
  { 
    *(flux+0) = (*(flux+0))/nn;
    *(flux+1) = (*(flux+1))/nn;
    *(flux+2) = (*(flux+2))/nn;
  }

} /* void compute_flux( lattice_ptr lattice, double *flux, int subs) */


#if STORE_UEQ
//void compute_max_ueq( lattice_ptr lattice, double *max_u)
void compute_max_ueq( 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[0][n].bc_type & BC_SOLID_NODE))
    {
      if( fabs( lattice->ueq[n].u[0]) > *(max_u))
      {
        *max_u = fabs( lattice->ueq[n].u[0]);
      }
      if( fabs( lattice->ueq[n].u[1]) > *(max_u+1))
      {
        *(max_u+1) = fabs( lattice->ueq[n].u[1]);
      }
    }
  }
} /* void compute_max_ueq( lattice_ptr lattice, double *max_u) */

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

//void compute_ave_ueq( lattice_ptr lattice, double *ave_u)
void compute_ave_ueq( 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[0][n].bc_type & BC_SOLID_NODE))
    {
      *ave_u += fabs( lattice->ueq[n].u[0]);
      *(ave_u+1) += fabs( lattice->ueq[n].u[1]);
      nn++;
    }
  }
  if( nn != 0)
  {
    *ave_u     = (*ave_u)/nn;
    *(ave_u+1) = (*(ave_u+1))/nn;
  }

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

#endif /* STORE_UEQ */

//void compute_vorticity( 
//       lattice_ptr lattice, int i, int j, int n, double *vor, int subs)
void compute_vorticity( 
       lattice_ptr lattice, int i, int j, int n, double *vor, int subs)
{
  double duyx, duxy;
  int    nn[5]; // Indices of neighbors;

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

  int    ip, in,
         jp, jn;

  ip = ( i<ni-1)?(i+1):(0   );
  in = ( i>0   )?(i-1):(ni-1);
  jp = ( j<nj-1)?(j+1):(0   );
  jn = ( j>0   )?(j-1):(nj-1);

  nn[1] = j *ni + ip;
  nn[2] = jp*ni + i ;
  nn[3] = j *ni + in;
  nn[4] = jn*ni + i ;

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

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

//printf("compute_vorticity() -- "
//    "n=%d, (i,j)=(%d,%d), duyx=%f, "
//    "nn1 = %d, nn3 = %d,"
//    "bc1 = %d, bc3 = %d"
//    "\n", 
//    n, i, j, duyx, 
//    nn[1], nn[3],
//    lattice->bc[subs][nn[1]].bc_type, 
//    lattice->bc[subs][nn[3]].bc_type);

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

  if( duxy*duyx != 0.)
  {
    *vor = duyx - duxy;
  }
  else
  {
    *vor = 0.;
  }

} /* 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, int subs)
void compute_max_vor( 
       lattice_ptr lattice, double *max_vor_p, double *max_vor_n, int subs)
{
  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[subs][n].bc_type == BC_FLUID_NODE)
    {
      compute_vorticity( lattice,
                         n%lattice->param.LX,
                         n/lattice->param.LX,
                         n,
                         &vor,
                         subs  );
      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[subs][n].bc_type == 0) */

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

#if 0 && VERBOSITY_LEVEL > 0
  printf("compute_max_vor() -- nnz = %d.  nnz/NumNodes = %f\n", 
    nnz, (double)nnz/(double)lattice->NumNodes);
#endif /* 0 && VERBOSITY_LEVEL > 0 */

} /* void compute_max_vor( lattice_ptr lattice, double *max_vor_p, ... */

//void compute_ave_vor(
//       lattice_ptr lattice,double *ave_vor_p,double *ave_vor_n, int subs)
void compute_ave_vor( 
       lattice_ptr lattice, double *ave_vor_p, double *ave_vor_n, int subs)
{
  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[subs][n].bc_type == BC_FLUID_NODE)
    {
      compute_vorticity( lattice,
                         n%lattice->param.LX,
                         n/lattice->param.LX,
                         n,
                         &vor,
                         subs  );
      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[subs][n].bc_type == 0) */