latman.c

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

          }
#else /* !( INAMURO_SIGMA_COMPONENT) */
          *macro_var_ptr++ = lattice->param.rho_A[subs];
#endif /* INAMURO_SIGMA_COMPONENT */
          break;
        }
        case IC_UNIFORM_RHO_B:
        {
#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_B[subs]
                 *((1./3.)*(lattice->param.tau[subs]-.5))
                 *K;
        if( lattice->param.incompressible)
        {
          drho = drho/lattice->param.rho_B[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_B[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_B[subs]
                 *((1./3.)*(lattice->param.tau[subs]-.5))
                 *K;

        if( lattice->param.incompressible)
        {
          drho = drho/lattice->param.rho_B[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_B[subs] + drho/2.) - m*j;

        if(    lattice->param.pressure_n_out[0] == 1
            && lattice->param.velocity_s_in[0]  == 1)
        {
          lattice->param.rho_out = 
              ( lattice->param.rho_B[subs] + drho/2.) 
          - m*( lattice->param.LY-1);
        }

#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_B[subs];

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

          }
          else // subs==1
          {
            *macro_var_ptr++ = 0.;
          }
#else /* !( INAMURO_SIGMA_COMPONENT) */
          *macro_var_ptr++ = lattice->param.rho_B[subs];
#endif /* INAMURO_SIGMA_COMPONENT */
          break;
        }
        case IC_UNIFORM_RHO_IN:
        {
#if INAMURO_SIGMA_COMPONENT
          if( subs==0)
          {
            *macro_var_ptr++ = lattice->param.rho_in;
          }
          else // subs==1
          {
            *macro_var_ptr++ = lattice->param.rho_in;
          }
#else /* !( INAMURO_SIGMA_COMPONENT) */
          *macro_var_ptr++ = lattice->param.rho_in;
#endif /* INAMURO_SIGMA_COMPONENT */
          break;
        }
        case IC_BUBBLE:
        {
#if INAMURO_SIGMA_COMPONENT
          if( subs==0)
          {
            *macro_var_ptr++ = lattice->param.rho_A[subs];
          }
          else
          {
            if( (i-lattice->param.x0)*(i-lattice->param.x0) 
              + (j-lattice->param.y0)*(j-lattice->param.y0) 
              < lattice->param.r0*lattice->param.r0)
            {
              *macro_var_ptr++ = lattice->param.rho_sigma;
            }
            else
            {
              *macro_var_ptr++ = 0.;
            }
          }
#else /* !( INAMURO_SIGMA_COMPONENT) */
          if( (i-lattice->param.x0)*(i-lattice->param.x0) 
            + (j-lattice->param.y0)*(j-lattice->param.y0) 
            < lattice->param.r0*lattice->param.r0)
          {
            *macro_var_ptr++ = lattice->param.rho_A[subs];
          }
          else
          {
            *macro_var_ptr++ = lattice->param.rho_B[subs];
          }
#endif /* INAMURO_SIGMA_COMPONENT */
          break;
        }
        case IC_YIN_YANG:
        {
          if( i < lattice->param.LX/2)
          {
            if( (i-lattice->param.x0/2.)*(i-lattice->param.x0/2.) 
              + (j-lattice->param.y0)*(j-lattice->param.y0) 
              < lattice->param.r0*lattice->param.r0/4.)
            {
              *macro_var_ptr++ = lattice->param.rho_A[subs];
            }
            else
            {
              *macro_var_ptr++ = lattice->param.rho_B[subs];
            }
          }
          else
          {
            if( (i-3.*lattice->param.x0/2.)*(i-3.*lattice->param.x0/2.) 
              + (j-lattice->param.y0)*(j-lattice->param.y0) 
              < lattice->param.r0*lattice->param.r0/4.)
            {
              *macro_var_ptr++ = lattice->param.rho_B[subs];
            }
            else
            {
              *macro_var_ptr++ = lattice->param.rho_A[subs];
            }
          }

          break;
        }
        case IC_DIAGONAL:
        {
          if( i+j < lattice->param.LX)
          {
            *macro_var_ptr++ = lattice->param.rho_A[subs];
          }
          else
          {
            *macro_var_ptr++ = lattice->param.rho_B[subs];
          }
          break;
        }
        case IC_2X2_CHECKERS:
        {
          if( ( ( i < lattice->param.LX/2) && ( j < lattice->param.LX/2))
           || ( ( i >= lattice->param.LX/2) && ( j >= lattice->param.LX/2)) )
          {
            *macro_var_ptr++ = lattice->param.rho_A[subs];
          }
          else
          {
            *macro_var_ptr++ = lattice->param.rho_B[subs];
          }
          break;
        }
        case IC_STATIC:
        {
          if( NUM_FLUID_COMPONENTS==2)
          {
            if( ((double)rand()/(double)RAND_MAX) < lattice->param.cut)
            {
              *macro_var_ptr++ = lattice->param.rho_A[subs];
            }
            else
            {
              *macro_var_ptr++ = lattice->param.rho_B[subs];
            }
          }
          else if( NUM_FLUID_COMPONENTS==1)
          {
            *macro_var_ptr++ = lattice->param.rho_in
                             + ((double)rand()/(double)RAND_MAX);
          }
          else
          {
            printf(
              "init_lattice() -- "
              "Unhandled case NUM_FLUID_COMPONENTS = %d .  "
              "Exiting!\n", NUM_FLUID_COMPONENTS);
            exit(1);
          }
          break;
        }
        case IC_RECTANGLE:
        {
#if INAMURO_SIGMA_COMPONENT
#if ZHANG_AND_CHEN_ENERGY_TRANSPORT
          if( subs==0)
          {
            if( ( i >= lattice->param.x1) && ( j >= lattice->param.y1)
             && ( i <= lattice->param.x2) && ( j <= lattice->param.y2))
            {
              *macro_var_ptr++ = lattice->param.rho_A[subs];
            }
            else
            {
              *macro_var_ptr++ = lattice->param.rho_B[subs];
            }
          }
          else // subs==1
          {
            if( ( i >= lattice->param.x1) && ( j >= lattice->param.y1)
             && ( i <= lattice->param.x2) && ( j <= lattice->param.y2))
            {
              *macro_var_ptr++ = lattice->param.rho_sigma;
            }
            else
            {
              *macro_var_ptr++ = lattice->param.rho_sigma;
            }
          }
#else /* !( ZHANG_AND_CHEN_ENERGY_TRANSPORT) */
          if( subs==0)
          {
            *macro_var_ptr++ = lattice->param.rho_A[subs];
          }
          else
          {
            if( ( i >= lattice->param.x1) && ( j >= lattice->param.y1)
             && ( i <= lattice->param.x2) && ( j <= lattice->param.y2))
            {
              *macro_var_ptr++ = lattice->param.rho_sigma;
            }
            else
            {
//            if( ( lattice->param.LX-i-1 >= lattice->param.x1) 
//             && ( j >= lattice->param.y1)
//             && ( lattice->param.LX-i-1 <= lattice->param.x2) 
//             && ( j <= lattice->param.y2))
//            {
//            *macro_var_ptr++ = lattice->param.rho_sigma;
//            }
//            else
//            {
              *macro_var_ptr++ = 0.;
//            }
            }
          }
#endif /* ZHANG_AND_CHEN_ENERGY_TRANSPORT */
#else /* !( INAMURO_SIGMA_COMPONENT) */
          if( ( i >= lattice->param.x1) && ( j >= lattice->param.y1)
           && ( i <= lattice->param.x2) && ( j <= lattice->param.y2))
          {
            *macro_var_ptr++ = lattice->param.rho_B[subs];
          }
          else
          {
            *macro_var_ptr++ = lattice->param.rho_A[subs];
          }
#endif /* INAMURO_SIGMA_COMPONENT */
          break;
        }
        case IC_DOT:
        {
          if( i == lattice->param.x0 && j == lattice->param.y0)
          {
            *macro_var_ptr++ = lattice->param.rho_A[subs];
          }
          else
          {
            *macro_var_ptr++ = lattice->param.rho_B[subs];
          }
          break;
        }
        case IC_WOLF_GLADROW_DIFFUSION:
        {
#if INAMURO_SIGMA_COMPONENT
          if( subs==0)
          {
            *macro_var_ptr++ = lattice->param.rho_A[subs];
          }
          else // subs==1
          {
            y = (j - ( lattice->param.LY-1)/2.)/((lattice->param.LY-1)/200.);
            *macro_var_ptr++ = (1./(2.*sqrt(PI*15.)))*exp(-y*y/(4.*15.));
          }
#else /* !( INAMURO_SIGMA_COMPONENT) */
          printf("%s (%d) >> Unhandled case. Exiting!\n", __FILE__, __LINE__);
          exit(1);
#endif /* INAMURO_SIGMA_COMPONENT */
          break;
        }
        case IC_HYDROSTATIC:
        {
#if INAMURO_SIGMA_COMPONENT
          if( subs==0)
          {
            *macro_var_ptr++ = 
              lattice->param.rho_A[subs]
            -
            j*( lattice->param.rho_A[subs]
              - lattice->param.rho_B[subs] )
              / lattice->param.LY;
          }
          else // subs==1
          {
            *macro_var_ptr++ = 0.;
          }
#else /* !( INAMURO_SIGMA_COMPONENT) */
          printf("%s (%d) >> Unhandled case. Exiting!\n", __FILE__, __LINE__);
          exit(1);
#endif /* INAMURO_SIGMA_COMPONENT */
          break;
        }
        default:
        {
          printf( 
            "%s (%d) >> init_problem() -- Unhandled case  "
            "lattice->param.initial_condition = %d.  "
            "Exiting!\n", __FILE__, __LINE__,
            lattice->param.initial_condition );
          exit(1);
          break;
        }
      } /* switch( lattice->param.initial_condition) */
    } /* if( ( 1 || !( bc->bc_type & BC_SOLID_NODE)) ) */
    else
    {
      if( bc->bc_type & BC_SOLID_NODE)
      {
        //*macro_var_ptr++ = lattice->param.rho_A[subs];
        *macro_var_ptr++ = lattice->param.rho_in;
      }
      else
      {
        *macro_var_ptr++ = 0.;
      }
    } /* if( ( 1 || !( bc->bc_type & BC_SOLID_NODE)) ) else */

    // Set initial velocty.
    if( ( 0 || !( bc->bc_type & BC_SOLID_NODE)) )
    {
#if 1
      // u_x
      if(    lattice->param.ic_poisseuille
          && !lattice->periodic_x[subs])
      {
        // Poisseuille flow profile in the x- direction.  Assuming 
        // one-lattice-unit walls on both sides.  

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

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

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

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

#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) 
              - 
                (j-.5*( lattice->param.LY-2)-.5)
               *(j-.5*( lattice->param.LY-2)-.5)
              )
          ;
#endif
//printf("%s (%d) -- %d %f\n", __FILE__, __LINE__, j, *(macro_var_ptr-1));

      }
      else
      {
#if INITIALIZE_WITH_UX_IN
        *macro_var_ptr++ = lattice->param.ux_in;
#else /* !( INITIALIZE_WITH_UX_IN) */
        *macro_var_ptr++ = 0.;
#endif /* INITIALIZE_WITH_UX_IN */
      }

      // u_y
      if(    lattice->param.ic_poisseuille
          && !lattice->periodic_y[subs])
      {
        // Poisseuille flow profile in the vertical/y- direction.  Assuming 
        // one-lattice-unit walls on both sides.  

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

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

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