latman.c

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

        // u(x) = K(a^2-x^2)
        x = (i-.5) - a;
        *macro_var_ptr++ = K * ( a*a - x*x);

#if 0
        ( 1.5*( lattice->param.uy_in)
             /( .25*(lattice->param.LX-2)*(lattice->param.LX-2)) )
             *( 
                .25*( lattice->param.LX-2)*( lattice->param.LX-2) 
              - 
                (i-.5*( lattice->param.LX-2)-.5)
               *(i-.5*( lattice->param.LX-2)-.5)
              )
          ;
#endif
//printf("%s (%d) -- %d %f\n", __FILE__, __LINE__, i, *(macro_var_ptr-1));

      } /* if( lattice->param.ic_poisseuille) */

      else // !lattice->param.ic_poisseuille
      {
#if INITIALIZE_WITH_UY_IN
        *macro_var_ptr++ = lattice->param.uy_in;
#else /* !( INITIALIZE_WITH_UY_IN) */
        *macro_var_ptr++ = 0.;
#endif /* INITIALIZE_WITH_UY_IN */

      } /* if( lattice->param.ic_poisseuille) else */

#else
      *macro_var_ptr++ = ( lattice->param.ux_in + lattice->param.ux_out)/2.;
      *macro_var_ptr++ = ( lattice->param.uy_in + lattice->param.uy_out)/2.;
#endif
    }
    else
    {
      *macro_var_ptr++ = 0.;
      *macro_var_ptr++ = 0.;
    }

    bc++;

#if STORE_UEQ
    *ueq++ = 0.;
    *ueq++ = 0.;
#endif /* STORE_UEQ */

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

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

  // Compute initial feq.
  compute_feq( lattice);

 for( subs=0; subs<NUM_FLUID_COMPONENTS; subs++)
 {
  for( n=0; n<lattice->NumNodes; n++)
  {
    f = lattice->pdf[subs][n].f;
    feq = lattice->pdf[subs][n].feq;
    ftemp = lattice->pdf[subs][n].ftemp;
    if( 1 || !( bc[n].bc_type & BC_SOLID_NODE))
    {
#if 1
      // Copy feq to f.
      f[0]= feq[0];
      f[1]= feq[1];
      f[2]= feq[2];
      f[3]= feq[3];
      f[4]= feq[4];
      f[5]= feq[5];
      f[6]= feq[6];
      f[7]= feq[7];
      f[8]= feq[8];

      // Initialize ftemp.
      ftemp[0]= 0.; //*(f-9);
      ftemp[1]= 0.; //*(f-9);
      ftemp[2]= 0.; //*(f-9);
      ftemp[3]= 0.; //*(f-9);
      ftemp[4]= 0.; //*(f-9);
      ftemp[5]= 0.; //*(f-9);
      ftemp[6]= 0.; //*(f-9);
      ftemp[7]= 0.; //*(f-9);
      ftemp[8]= 0.; //*(f-9);

      f+=9;
#else
      // Debug info.  To track during the streaming step.

      // f
      f[0] = (double)n + 0./10.;// *(f-9);
      f[1] = (double)n + 1./10.;// *(f-9);
      f[2] = (double)n + 2./10.;// *(f-9);
      f[3] = (double)n + 3./10.;// *(f-9);
      f[4] = (double)n + 4./10.;// *(f-9);
      f[5] = (double)n + 5./10.;// *(f-9);
      f[6] = (double)n + 6./10.;// *(f-9);
      f[7] = (double)n + 7./10.;// *(f-9);
      f[8] = (double)n + 8./10.;// *(f-9);

      // ftemp
      ftemp[0] = (double)n;// *(f-9);
      ftemp[1] = (double)n;// *(f-9);
      ftemp[2] = (double)n;// *(f-9);
      ftemp[3] = (double)n;// *(f-9);
      ftemp[4] = (double)n;// *(f-9);
      ftemp[5] = (double)n;// *(f-9);
      ftemp[6] = (double)n;// *(f-9);
      ftemp[7] = (double)n;// *(f-9);
      ftemp[8] = (double)n;// *(f-9);

      f+=9;
#endif
    }
    else
    {
      // f
      f[0] = 0.;
      f[1] = 0.;
      f[2] = 0.;
      f[3] = 0.;
      f[4] = 0.;
      f[5] = 0.;
      f[6] = 0.;
      f[7] = 0.;
      f[8] = 0.;
          
      // ftemp
      ftemp[0] = 0.;
      ftemp[1] = 0.;
      ftemp[2] = 0.;
      ftemp[3] = 0.;
      ftemp[4] = 0.;
      ftemp[5] = 0.;
      ftemp[6] = 0.;
      ftemp[7] = 0.;
      ftemp[8] = 0.;
   
      f+=9;
    }

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

#if NON_LOCAL_FORCES
  force = lattice->force[subs][0].force;
  for( n=0; n<lattice->NumNodes; n++)
  {
    *force++ = 0.;
    *force++ = 0.;
    *force++ = 0.;
    *force++ = 0.;
  }
#endif /* NON_LOCAL_FORCES */

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

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

  //dump_pdf( lattice, /*time=*/ 0);

  //compute_macro_vars( lattice);

  //dump_macro_vars( lattice, /*time=*/ 0);

#if VERBOSITY_LEVEL > 0
  printf("init_problem() -- Problem initialized.\n");
#endif /* VERBOSITY_LEVEL > 0 */

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

// void destruct_lattice( struct lattice_struct *lattice)
//##############################################################################
//
// D E S T R U C T   L A T T I C E 
//
//  - Destruct lattice.
//
void destruct_lattice( struct lattice_struct *lattice)
{
  int subs;

  assert( lattice!=NULL);

 for( subs=0; subs<NUM_FLUID_COMPONENTS; subs++)
 {
  assert( lattice->pdf[subs]!=NULL);
  free(   lattice->pdf[subs]);

  assert( lattice->macro_vars[subs]!=NULL);
  free(   lattice->macro_vars[subs]);

  assert( lattice->bc[subs]!=NULL);
  free(   lattice->bc[subs]);
#if NON_LOCAL_FORCES
  assert( lattice->force[subs]!=NULL);
  free(   lattice->force[subs]);
#endif /* NON_LOCAL_FORCES */
 }

#if STORE_UEQ
  free(   lattice->ueq);
#endif /* STORE_UEQ */

  free(   lattice);

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

#if INAMURO_SIGMA_COMPONENT && STORE_BTC
void sigma_stuff( lattice_ptr lattice)
{
  int    btc_time;
  int    i, j, n;
  double btc_val01, btc_val02, btc_val03, btc_val04;
  int    width01,   width02,   width03,   width04;

  // Turn off concentration boundaries when time is up.
  if(    lattice->param.sigma_stop >= 0 
      && lattice->time > lattice->param.sigma_stop)
  {
    lattice->param.rho_sigma_in  = 0.;
    lattice->param.rho_sigma_out = 0.;

  } /* if( lattice->param.sigma_stop >= 0 && ... */

  // Accumulate break through curve.
  if( !( lattice->param.sigma_btc_rate <= 0 || lattice->FlowDir==0))
  {
    if(    lattice->param.sigma_start <= lattice->time
       &&  ( 
              lattice->param.sigma_stop  < 0
            ||
              lattice->param.sigma_stop  >= lattice->time
           )
       && !( lattice->time % lattice->param.sigma_btc_rate))
    {
      btc_time = ( lattice->time 
                 - ((lattice->param.sigma_start>0)
                   ?(lattice->param.sigma_start)
                   :(0))                           ) 
                 / lattice->param.sigma_btc_rate;

//printf("%s (%d) >> btc_time = %d (%d)\n",
//    __FILE__,__LINE__,btc_time,lattice->SizeBTC);

      btc_val01 = 0.;
      btc_val02 = 0.;
      btc_val03 = 0.;
      btc_val04 = 0.;

      if( lattice->FlowDir == /*Horizontal*/1)
      {
        i = ( (lattice->param.sigma_btc_spot >= 0)
            ? (lattice->param.sigma_btc_spot)
            : (lattice->param.LX-2));
        width01 = 0;
        width02 = 0;
        width03 = 0;
        width04 = 0;
        for( j=0; j<lattice->param.LY; j++)
        {
          //
          // Concentration at sigma_spot-1
          //
          n = i-1 + j*lattice->param.LX;
          if( !( lattice->bc[0][n].bc_type & BC_SOLID_NODE))
          {
            btc_val01 += lattice->macro_vars[1][n].rho;
            width01++;

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

          //
          // Concentration at sigma_spot
          //
          n = i   + j*lattice->param.LX;
          if( !( lattice->bc[0][n].bc_type & BC_SOLID_NODE))
          {
            btc_val02 += lattice->macro_vars[1][n].rho;
            width02++;

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

          //
          // Concentration at sigma_spot+1
          //
          n = i+1 + j*lattice->param.LX;
          if( !( lattice->bc[0][n].bc_type & BC_SOLID_NODE))
          {
            btc_val03 += lattice->macro_vars[1][n].rho;
            width03++;

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

          //
          // Velocity at sigma_spot
          //
          n = i + j*lattice->param.LX;
          if( !( lattice->bc[0][n].bc_type & BC_SOLID_NODE))
          {
            btc_val04 += lattice->macro_vars[1][n].u[0];
            width04++;

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


        } /* for( j=0; j<lattice->param.LY; j++) */

        btc_val01 /= width01;
        btc_val02 /= width02;
        btc_val03 /= width03;
        btc_val04 /= width04;

      } /* if( lattice->FlowDir == 1) */

      else if( lattice->FlowDir == /*Vertical*/2)
      {
        j = ( (lattice->param.sigma_btc_spot >= 0)
            ? (lattice->param.sigma_btc_spot)
            : (lattice->param.LY-2));
        width01 = 0;
        width02 = 0;
        width03 = 0;
        width04 = 0;
        for( i=0; i<lattice->param.LX; i++)
        {
          //
          // Concentration at sigma_spot-1
          //
          n = i + (j-1)*lattice->param.LX;
          if( !( lattice->bc[0][n].bc_type & BC_SOLID_NODE))
          {
            btc_val01 += lattice->macro_vars[1][n].rho;
            width01++;

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

          //
          // Concentration at sigma_spot
          //
          n = i + (j  )*lattice->param.LX;
          if( !( lattice->bc[0][n].bc_type & BC_SOLID_NODE))
          {
            btc_val02 += lattice->macro_vars[1][n].rho;
            width02++;

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

          //
          // Concentration at sigma_spot+1
          //
          n = i + (j+1)*lattice->param.LX;
          if( !( lattice->bc[0][n].bc_type & BC_SOLID_NODE))
          {
            btc_val03 += lattice->macro_vars[1][n].rho;
            width03++;

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

          //
          // Velocity at sigma_spot
          //
          n = i + j*lattice->param.LX;
          if( !( lattice->bc[0][n].bc_type & BC_SOLID_NODE))
          {
            btc_val04 += lattice->macro_vars[1][n].u[1];
            width04++;

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


        } /* for( i=0; i<lattice->param.LX; i++) */

        btc_val01 /= width01;
        btc_val02 /= width02;
        btc_val03 /= width03;
        btc_val04 /= width04;

      } /* else if( lattice->FlowDir == 2) */

      else // Flow direction is undetermined.
      {
        // Unhandled case.
        // TODO: Warning message?

      } /* if( lattice->FlowDir == ?) else */

//printf("%s (%d) >> btc_val01 = %f\n",__FILE__,__LINE__,btc_val01);
      lattice->param.sigma_btc[5*btc_time+0] = lattice->time;
      lattice->param.sigma_btc[5*btc_time+1] = btc_val01;
      lattice->param.sigma_btc[5*btc_time+2] = btc_val02;
      lattice->param.sigma_btc[5*btc_time+3] = btc_val03;
      lattice->param.sigma_btc[5*btc_time+4] = btc_val04;

    } /* if( !( lattice->time % lattice->param.sigma_btc_rate)) */

  } /* if( lattice->param.sigma_start <= lattice->time) */

} /* void sigma_stuff( lattice_ptr lattice) */
#endif /* INAMURO_SIGMA_COMPONENT && STORE_BTC */

// int get_sizeof_lattice_structure( lattice_ptr lattice)
//##############################################################################
//
// G E T _ S I Z E O F _ L A T T I C E _ S T R U C T U R E
//
//  - Return size of struct lattice_struct in bytes.
//
int get_sizeof_lattice_structure( lattice_ptr lattice)
{
  return sizeof( struct lattice_struct);

} /* int get_sizeof_lattice_structure( lattice_ptr lattice) */

// int get_sizeof_lattice( lattice_ptr lattice)
//##############################################################################
//
// G E T _ S I Z E O F _ L A T T I C E
//
//  - Return size of lattice in bytes.
//
int get_sizeof_lattice( lattice_ptr lattice)
{
  return
      sizeof(int)
    + sizeof(int)
    + sizeof(int)
    + sizeof(int)*NUM_FLUID_COMPONENTS
    + sizeof(int)*NUM_FLUID_COMPONENTS

    + sizeof(struct param_struct)

    + lattice->NumNodes
    * ( 
        NUM_FLUID_COMPONENTS*sizeof(struct pdf_struct)
      + NUM_FLUID_COMPONENTS*sizeof(struct macro_vars_struct)
      + NUM_FLUID_COMPONENTS*sizeof(struct bc_struct) 
#if NON_LOCAL_FORCES
      + NUM_FLUID_COMPONENTS*sizeof(struct force_struct) 
#endif /* NON_LOCAL_FORCES */
#if STORE_UEQ
      + sizeof(struct ueq_struct) 
#endif /* STORE_UEQ */
#if POROUS_MEDIA
      + sizeof(struct ns_struct) 
#endif /* POROUS_MEDIA */
      );

} /* int get_sizeof_lattice( lattice_ptr lattice) */

// int get_num_active_nodes( lattice_ptr lattice)
//##############################################################################
//
// G E T _ N U M B E R _ A C T I V E _ N O D E S
//
//  - Return number of active nodes.
//
int get_num_active_nodes( lattice_ptr lattice)
{
  int n, k;

  k = 0;

  for( n=0; n<lattice->NumNodes; n++)
  {
    if( lattice->bc[0][n].bc_type != INACTIVE_NODE)
    {
      k++;
    }
  }

  return k;
  
} /* int get_num_active_nodes( lattice_ptr lattice) */