latman.c

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

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

      fclose( in);

    } /* if( in = fopen( filename, "r+")) */

    else /* !( in = fopen( filename, "r+")) */
    {
      // Can't read ns.bmp file, so use default values.
      printf("%s %d >> WARNING: Can't read \"%s\". "
          "Using default ns values.\n",__FILE__,__LINE__,filename);
    } /* if( in = fopen( filename, "r+")) else */


      break;
    }
    case 2:
    {
    // Allocate space for ns values.
    (*lattice)->ns = 
      (struct ns_struct*)malloc( (*lattice)->NumNodes*sizeof(struct ns_struct));
    if( (*lattice)->ns==NULL)
    {
      printf(
        "construct_lattice() -- ERROR:  "
        "Attempt to allocate %d struct ns_struct types failed.  "
        "Exiting!\n",
        (*lattice)->NumNodes
        );
      exit(1);
    }

    // Try to read ns<LX>x<LY>.bmp file.
    sprintf( filename, "./in/ns%dx%d.bmp", 
           (*lattice)->param.LX, 
           (*lattice)->param.LY);
    if( in = fopen( filename, "r+"))
    {
      printf("%s %d >> Reading file \"%s\".\n",__FILE__,__LINE__,filename);
      bmp_read_header( in, &bmih);
      for( n=0; n<(*lattice)->NumNodes; n++)
      {
        bmp_read_entry( in, bmih, &r, &g, &b);

        // Verify grayscale.
        if(    (double)r != (double)g 
            || (double)g != (double)b 
            || (double)r != (double)b)
        {
          printf(
            "%s %d >> latman.c: construct_lattice() -- "
            "n=%d:  [ r g b] = [ %3u %3u %3u]\n",__FILE__,__LINE__,
            n, (unsigned int)r%256, (unsigned int)g%256, (unsigned int)b%256);
          printf(
            "%s %d >> latman.c: construct_lattice() -- "
            "ERROR: File %s needs to be grayscale. "
            "Exiting!\n",__FILE__,__LINE__, filename);
          exit(1);
        }

        if( ((unsigned int)r%256) != 0 && ((unsigned int)r%256) != 255 )
        {
          printf(
            "%s %d >> latman.c: construct_lattice() -- "
            "ERROR: File %s needs to be black and white. "
            "Exiting!\n",__FILE__,__LINE__, filename);
          exit(1);
        }

        // Assign ns value.
        if( ((unsigned int)r%256) == 0)
        {
          (*lattice)->ns[n].ns = (*lattice)->param.ns;
        }
        else
        {
          (*lattice)->ns[n].ns = 0.;
        }
#if 0 && VERBOSITY_LEVEL>0
        printf("%s %d >> n=%d, ns=%f\n", 
            __FILE__, __LINE__, n, (*lattice)->ns[n].ns);
#endif /* 1 && VERBOSITY_LEVEL>0 */

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

      fclose( in);

    } /* if( in = fopen( filename, "r+")) */

    else /* !( in = fopen( filename, "r+")) */
    {
      // Can't read ns.bmp file, so use default values.
      printf("%s %d >> WARNING: Can't read \"%s\". "
          "Using default ns values.\n",__FILE__,__LINE__,filename);
    } /* if( in = fopen( filename, "r+")) else */
      break;
    }
    default:
    {
      printf("%s %d >> construct_lattice() -- Unhandled case: "
        "ns_flag = %d . (Exiting!)\n",
        __FILE__,__LINE__,(*lattice)->param.ns_flag < 0.);
      exit(1);
      break;
    }

  } /* switch( (*lattice)->param.ns_flag) */
#endif /* POROUS_MEDIA */

#if INAMURO_SIGMA_COMPONENT && DETERMINE_FLOW_DIRECTION
  //
  // Try to determine the direction of flow.  
  //
  // NOTE: This determination informs the breakthrough curve mechanism which
  // should be used in a simple situation with either pressure/velocity
  // boundaries driving the flow in one direction or gravity driving the flow
  // in one direction.  If the direction of flow cannot be determined, FlowDir
  // will be set to indeterminate (=0) and a BTC will not be stored.
  //
  // NOTE: This determination also informs the sigma slip boundary which
  // should only be used in the simple situation of flow through a channel
  // where geometry is trivial and the direction of flow is obvious.
  //
  if( // Pressure/Velocity boundaries moving the flow vertically.
 (
    ( (*lattice)->param.pressure_n_in[0] && (*lattice)->param.pressure_s_out[0])
  ||
    ( (*lattice)->param.pressure_n_out[0] && (*lattice)->param.pressure_s_in[0])
  ||
    ( (*lattice)->param.velocity_n_in[0] && (*lattice)->param.velocity_s_out[0])
  ||
    ( (*lattice)->param.velocity_n_out[0] && (*lattice)->param.velocity_s_in[0])
 )
&&
!(
    ( (*lattice)->param.pressure_e_in[0] && (*lattice)->param.pressure_w_out[0])
  ||
    ( (*lattice)->param.pressure_e_out[0] && (*lattice)->param.pressure_w_in[0])
  ||
    ( (*lattice)->param.velocity_e_in[0] && (*lattice)->param.velocity_w_out[0])
  ||
    ( (*lattice)->param.velocity_e_out[0] && (*lattice)->param.velocity_w_in[0])
 )
    )
  {
    (*lattice)->FlowDir = /*Vertical*/2;
  }
  else if( // Pressure/Velocity boundaries moving the flow horizontally.
!( 
    ( (*lattice)->param.pressure_n_in[0] && (*lattice)->param.pressure_s_out[0])
  ||
    ( (*lattice)->param.pressure_n_out[0] && (*lattice)->param.pressure_s_in[0])
  ||
    ( (*lattice)->param.velocity_n_in[0] && (*lattice)->param.velocity_s_out[0])
  ||
    ( (*lattice)->param.velocity_n_out[0] && (*lattice)->param.velocity_s_in[0])
 )
&&
 (
    ( (*lattice)->param.pressure_e_in[0] && (*lattice)->param.pressure_w_out[0])
  ||
    ( (*lattice)->param.pressure_e_out[0] && (*lattice)->param.pressure_w_in[0])
  ||
    ( (*lattice)->param.velocity_e_in[0] && (*lattice)->param.velocity_w_out[0])
  ||
    ( (*lattice)->param.velocity_e_out[0] && (*lattice)->param.velocity_w_in[0])
 )
      )
  {
    (*lattice)->FlowDir = /*Horizontal*/1;
  }
  else
  {
    if( // Gravity driving flow vertically.
        (*lattice)->param.gforce[0][1] != 0.
     && (*lattice)->param.gforce[0][0] == 0.
      )
    {
      (*lattice)->FlowDir = /*Vertical*/2;
    }
    else if( // Gravity driving flow horizontally.
             (*lattice)->param.gforce[0][0] != 0.
          && (*lattice)->param.gforce[0][1] == 0.
           )
    {
      (*lattice)->FlowDir = /*Horizontal*/1;
    }
    else // Cannot determine direction of flow.
    {
      (*lattice)->FlowDir = /*Indeterminate*/0;
    }
  }
#if INITIALIZE_WITH_UX_IN
  (*lattice)->FlowDir = /*Horizontal*/1;
#endif /* INITIALIZE_WITH_UX_IN */
#if INITIALIZE_WITH_UY_IN
  (*lattice)->FlowDir = /*Vertical*/2;
#endif /* INITIALIZE_WITH_UY_IN */
#endif /* INAMURO_SIGMA_COMPONENT && DETERMINE_FLOW_DIRECTION */

#if INAMURO_SIGMA_COMPONENT && STORE_BTC
  // Allocate space for a break through curve if necessary.
  if( (*lattice)->param.sigma_btc_rate > 0 && (*lattice)->FlowDir!=0)
  {
    // Compute "size" of break through curve: number of readings
    // to store.
    (*lattice)->SizeBTC = 
      (int)ceil((double)(
                         (
        (*lattice)->NumTimeSteps
      -
        (((*lattice)->param.sigma_start>0)
        ?((*lattice)->param.sigma_start)
        :(0))                            
                         ) 
       / (*lattice)->param.sigma_btc_rate))+1;

    //
    // Allocate 4*SizeBTC elements.
    //
    // Readings will come in groups of four (r1,r2,r3,r4):
    //
    //   r0: Timestep.
    //
    //   r1: Concentration at sigma_spot-1
    //
    //   r2: Concentration at sigma_spot-0
    //
    //   r3: Concentration at sigma_spot+1
    //
    //   r4: Velocity in direction of flow.
    //
    // Then
    //
    //   Cf = ( C*v - D*dCdx)/v = ( (r2+r3)/(2*r4) - D*(r3-r2)) / r4
    //
    (*lattice)->param.sigma_btc =
      ( double*)malloc( 5*(*lattice)->SizeBTC*sizeof(double));

  } /* if( sigma_btc_rate > 0) */
#endif /* INAMURO_SIGMA_COMPONENT && STORE_BTC */

  dump_params( *lattice);

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

// void init_problem( struct lattice_struct *lattice)
//##############################################################################
//
// I N I T   P R O B L E M 
//
//  - Initialize the problem on the lattice.
//
//    - Set the initial density and velocity.
//
//    - Compute the initial feq.
//
void init_problem( struct lattice_struct *lattice)
{
#if WRITE_CHEN_DAT_FILES
  FILE   *o;
  char   filename[1024];
#endif /* WRITE_CHEN_DAT_FILES */
  int    n, i, j;
  double a, x, u_max, K, drho, m;
  double *macro_var_ptr;
  double *f, *feq, *ftemp;
#if STORE_UEQ
  double *ueq;
#endif /* STORE_UEQ */
#if NON_LOCAL_FORCES
  double *force;
#endif /* NON_LOCAL_FORCES */
  bc_ptr bc;
  int    subs;
  double kappa;
  double ti;
  double y;
  if( lattice->param.initial_condition == IC_WOLF_GLADROW_DIFFUSION)
  {
    kappa = 1.*( lattice->param.tau[1] - .5);
    ti = 15./kappa;
    printf("IC_WOLF_GLADROW_DIFFUSION\n");
    printf("  ti = 15./kappa = 15./%f = %f\n", kappa, ti);
    printf("  tf = 75./kappa = 75./%f = %f\n", kappa, 75./kappa);
    printf("  tf-ti = (75.-15.)/kappa = 60./%f = %f\n", kappa, 60./kappa);
    printf("\n");
  }

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

#if WRITE_CHEN_DAT_FILES
  //
  // Create empty chen_*.dat files.
  //
  sprintf( filename, "%s", "./out/chen_xyrho.dat");
  if( !( o = fopen( filename,"w+")))
  {
    printf("Error creating \"%s\".  Exiting!\n", filename);
    exit(1);
  }
  fclose( o);

  sprintf( filename, "%s", "./out/chen_xy_ux_uy.dat");
  if( !( o = fopen( filename,"w+")))
  {
    printf("Error creating \"%s\".  Exiting!\n", filename);
    exit(1);
  }
  fclose( o);

  sprintf( filename, "%s", "./out/chen_time.dat");
  if( !( o = fopen( filename,"w+")))
  {
    printf("Error creating \"%s\".  Exiting!\n", filename);
    exit(1);
  }
  fclose( o);
 
#endif /* WRITE_CHEN_DAT_FILES */

#if 0 // Want to allow mix of velocity and pressure boundaries.
 if( lattice->param.ic_poisseuille)
 {
   if( lattice->param.uy_in != lattice->param.uy_out)
   {
     printf("\n");
     printf("\n");
     printf("%s (%d) -- ERROR: "
       "Need uy_in == uy_out to initialize with poisseuille profile.  "
       "Exiting.\n",
       __FILE__,__LINE__);
     printf("\n");
     printf("\n");
     exit(1);

   } /* if( lattice->param.uy_in != lattice->param.uy_out) */

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

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

  macro_var_ptr = &( lattice->macro_vars[subs][0].rho);
  bc            = lattice->bc[subs];
#if STORE_UEQ
  ueq           = lattice->ueq[0].u;
#endif /* STORE_UEQ */

  for( n=0; n<lattice->NumNodes; n++)
  {
    i = n%lattice->param.LX;
    j = n/lattice->param.LX;

    // Set initial density.
    if( ( 1 || !( bc->bc_type & BC_SOLID_NODE)) )
    {
      switch( lattice->param.initial_condition)
      {
        case IC_UNIFORM_RHO_A:
        {
#if INAMURO_SIGMA_COMPONENT
          if( subs==0)
          {
      if( 1 && lattice->param.ic_poisseuille)
      {
        if( !lattice->periodic_x[subs])
        {
        // Density gradient corresponding to the desired velocity.
        // Based on formula for poisseuille velocity profile
        //
        //   u(x) = K*( a^2 - x^2)
        //
        // where, in terms of the density/pressure gradient,
        //
        //   K = dP/(2 L rho0 nu) ==> drho = 6 L rho0 nu K
        //
        // using the equation of state rho = 3 P.
        //

        // u_max = (3/2)ux_in
        u_max = 1.5*(lattice->param.ux_in);

        // a = .5*(LY-2) .  
        a     = .5*(lattice->param.LY-2);

        //     u(y) = K(a^2-y^2) 
        // ==> u_max = Ka^2 
        // ==> K = u_max/a^2
        K     = u_max/(a*a);

        drho = 6.*lattice->param.LX
                 *lattice->param.rho_A[subs]
                 *((1./3.)*(lattice->param.tau[subs]-.5))
                 *K;
        if( lattice->param.incompressible)
        {
          drho = drho/lattice->param.rho_A[subs];
        }
        //drho = .0188121666666667;
        //drho = .62680456666658;
        //drho = 6.2680456666658;
        //printf("%s (%d) -- drho = %f\n",__FILE__,__LINE__,drho);

        m = drho/(lattice->param.LX-1);

        *macro_var_ptr++ = 
          ( lattice->param.rho_A[subs] + drho/2.) - m*i;

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

        }
        else if( !lattice->periodic_y[subs])
        {
        // Density gradient corresponding to the desired velocity.
        // Based on formula for poisseuille velocity profile
        //
        //   u(x) = K*( a^2 - x^2)
        //
        // where, in terms of the density/pressure gradient,
        //
        //   K = dP/(2 L rho0 nu) ==> drho = 6 L rho0 nu K
        //
        // using the equation of state rho = 3 P.
        //

        // u_max = (3/2)uy_in
        u_max = 1.5*(lattice->param.uy_in);

        // a = .5*(LX-2) .  
        a     = .5*(lattice->param.LX-2);

        //     u(x) = K(a^2-x^2) 
        // ==> u_max = Ka^2 
        // ==> K = u_max/a^2
        K     = u_max/(a*a);

        drho = 6.*lattice->param.LY
                 *lattice->param.rho_A[subs]
                 *((1./3.)*(lattice->param.tau[subs]-.5))
                 *K;
        if( lattice->param.incompressible)
        {
          drho = drho/lattice->param.rho_A[subs];
        }
        //drho = .0188121666666667;
        //drho = .62680456666658;
        //drho = 6.2680456666658;
        //printf("%s (%d) -- drho = %f\n",__FILE__,__LINE__,drho);

        m = drho/(lattice->param.LY-1);

        *macro_var_ptr++ = 
          ( lattice->param.rho_A[subs] + drho/2.) - m*j;

#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));


        }
        else
        {
        }

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

      else // !lattice->param.ic_poisseuille
      {
            *macro_var_ptr++ = lattice->param.rho_A[subs];

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

          }
          else // subs==1
          {
            *macro_var_ptr++ = lattice->param.rho_sigma;