lattice.h

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

//##############################################################################
//
// Copyright (C), 2005, Michael Sukop and Danny Thorne
//
// lattice.h
//
//  - Lattice structure and supporting structures for lb2d_prime.
//

// struct pdf_struct
//
//  - Structure to hold the particle distribution functions.
//
struct pdf_struct
{
  double feq[   /*NUM_DIRS*/ 9];
  double f[     /*NUM_DIRS*/ 9];
  double ftemp[ /*NUM_DIRS*/ 9];

};
typedef struct pdf_struct *pdf_ptr;

// struct macro_vars_struct
//
//  - Structure to hold the macroscopic variables.
//
struct macro_vars_struct
{
  double rho;
  double u[   /*NUM_DIMS*/ 2];

};

#if STORE_UEQ
// struct ueq_struct
//
//  - Structure to hold the composite velocity used to compute feq.
//
struct ueq_struct
{
  double u[   /*NUM_DIMS*/ 2];

};
#endif /* STORE_UEQ */

// struct bc_struct
//
//  - Structure to hold boundary condition information.
//
//  - Currently it is just an integer flag. I'm putting it inside
//    a struct in case we want a more complicated boundary 
//    conditions mechanism later.
//
struct bc_struct
{
  int bc_type;

};
typedef struct bc_struct *bc_ptr;

#if DO_NOT_STORE_SOLIDS
// struct node_struct
//
//  - Structure to hold node information.
//
struct node_struct
{
  int i, j; // Lattice coordinates of node.

  int n; // Node index.

  int nn[8]; // Indices of neighboring nodes.

};
typedef struct node_struct *node_ptr;
#endif /* DO_NOT_STORE_SOLIDS */

#if NON_LOCAL_FORCES
struct force_struct
{
  double force[  /*NUM_DIMS*/2];
  double sforce[ /*NUM_DIMS*/2];
};
#endif /* NON_LOCAL_FORCES */
#if POROUS_MEDIA
struct ns_struct
{
  double ns;
};
#endif /* POROUS_MEDIA */

// struct param_struct
//
// Structure to hold parameters that specify the configuration of the problem
// to solve.
//
// Values for these parameters are read from file "./in/params.in". (Note that
// there are also some compile-time configuration options in file
// "./in/flags.h".)
//
// Each parameter is preceeded by documentation consisting of three kinds of
// fields: "Param", "Type", and "Comments". There may be multiple "Param"
// fields for sets of related parameters, e.g. LX and LY. There is only one
// occurance of the "Type" field and "Comments" field for each (set of)
// parameter(s). The "Param" field(s) and the "Type" field are each one line.
// The "Comments" field can be multiple lines.
//
// Distinct fluids are referred to both as components and as substances.
// Normally in conversation we talk about fluid components, but variables are
// indexed by "subs". The term "component" can be used to refer to a fluid
// component or a solute component (c.f. Inamuro) as well as the more broadly
// conventional mathematical usage of a vector component. The reader should be
// wary of context, therefore.
//
// In the "Param" field(s) of the documentation preceeding each parameter, the
// parameters that are arrays are stated with typical and/or suggestive index
// names between square brackets. For example, tau[subs] is written to suggest
// the fact that tau contains an array of values, one for each substance or
// fluid component, where "subs" is a typical index name for indexing fluid
// components. The commentary section ("Comments" field) for each parameter is
// explicit about the meaning of the indices.
//
// A common mistake made when configuring a problem is to specify in the file
// "params.in" an integer number for a type double parameter or vice versa.
// This can cause mystifying symptoms. Be careful to include a decimal point in
// all type double parameter values in "params.in" regardless of whether there
// is a fractional part to the value. Likewise, do not include decimal points
// in type integer parameter values. Remember! If things have been going
// smoothly and then suddenly go bad in weird ways (sorry I can't recollect any
// specific examples) be sure to check this issue in "params.in". Note that
// file "params.dat" is output by the code mostly as a sanity check: if values
// in "params.dat" do not match values in "params.in", the double versus
// integer parameter type issue may be the reason.
//
// Here are some useful Unix/Linux command lines for processing this 
// documentation (the "params.txt" file):
//
//  >> grep Param params.txt
//
//  >> grep Comment params.txt
//
//  >> TODO
//
//  >> TODO
//
// Here is the command line that generates the "params.txt" file from the
// "lattice.h" file:
//
//  >> cat ./src/lattice.h | \
//     sed -n '/^\/\/ struct param_struct/,/\/\* struct param_struct/p' | \
//     grep "\/\/" | \
//     sed 's/[ ]*\/\/ //; s/[ ]*\/\/$//' > ./doc/params.txt
//
// If you are developing under Windows, I recommend Cygwin
// 
//   http://www.cygwin.com
//
// or something comparable. The Unix/Linux command line is indispensable. See
//
//   http://www.cs.usfca.edu/~parrt/course/601/lectures/unix.util.html
//   http://www.student.northpark.edu/pemente/sed/sed1line.txt
//
// for example.
//
// -----------------------------------------------------------------------------
//
struct param_struct
{
  // Param: LX
  // Param: LY
  // Type: int
  // Comments: LX is the number of lattice nodes in x-direction. LY is the
  // number of lattice nodes in y-direction. These are be used to construct
  // the file name of the BMP file "<LX>x<LY>.bmp" containing the domain. For
  // example, if LX=160 and LY=120, then the file name is "160x120.bmp".
  //
  int    LX,
         LY;

  // Param: length_scale
  // Type: int 
  // Comments: length_scale is a characteristic length. This is used for
  // computing Re.
  //
  int    length_scale;

  // Param: NumFrames
  // Type: int
  // Comments: NumFrames is the number of frames to output.
  //
  int    NumFrames;

  // Param: FrameRate
  // Type: int
  // Comments: FrameRate is the number of timesteps per frame.
  //
  int    FrameRate;

  // Param: NumTimeSteps
  // Type: int
  // Comments: NumTimeSteps is the number of time steps computed as
  // NumFrames*FrameRate.
  //
  int    NumTimeSteps;

  // Param: tau[subs]
  // Type: double*
  // Comments: tau is the relaxation time(s). If NUM_FLUID_COMPONENTS==1, then
  // tau[0] is the relaxation time. If NUM_FLUID_COMPONENTS==2, then tau[0] is
  // the relaxation time for component 0 and tau[1] is the relaxation time for
  // component 1.
  //
  double tau[ NUM_FLUID_COMPONENTS];

  // Param: gforce[subs][dir]
  // Type: double**
  // Comments: gforce[subs][dir] is the gravitational (body) force(s). If
  // NUM_FLUID_COMPONENTS==1, then gforce[0][dir] holds the gravitational force.
  // If NUM_FLUID_COMPONENTS==2, then gforce[0][dir] holds the gravitational
  // force for component 0 and gforce[1][dir] holds the gravitational force for
  // component 1. In either case, gforce[subs][0] is the gravitational force
  // in the x-direction and gforce[subs][1] is the gravitational force in the
  // y-direction.
  //
  double gforce[   NUM_FLUID_COMPONENTS][2];

  // Param: end_grav[subs]
  // Type: int
  // Comments: end_grav[subs] is the timestep at which to end the gravitational
  // force. If NUM_FLUID_COMPONENTS==1, then end_grav[0] holds the ending
  // timestep. If NUM_FLUID_COMPONENTS==2, then end_grav[0] holds the ending
  // timestep for fluid compenent 0 and end_grav[1] holds the ending timestep
  // for fluid component 1. This mechanism is only enforced when the
  // MANAGE_BODY_FORCE flag is turned on in "flags.h". It is implemented as a
  // call to the "manage_body_force()" function at each time step, so it is a
  // special feature that should generally be turned off. In the future this
  // mechanism should mature.
  //
  int    end_grav[ NUM_FLUID_COMPONENTS];

  // Param: buoyancy
  // Type: int
  // Comments: The "buoyancy" flag toggles buoyancy effects due to the density
  // of solute. This flag is only relevant when INAMURO_SIGMA_COMPONENT is
  // turned on in "flags.h". More recently, the buoyancy flag can be +1 or -1
  // to switch between solute induced buoyancy effects and temperature induced
  // buoyancy effects. Use the accessor functions get_buoyancy_flag(lattice),
  // get_buoyancy_sign(lattice) and get_buoyancy(lattice).
  //
  int    buoyancy;

  // Param: incompressible
  // Type: int
  // Comments: The "incompressible" flag toggles an incompressible lattice
  // Boltzmann method (c.f. TODO: cite). The usual method that we experiment
  // with is the common "weakly" compressible BGK method. For some things
  // (e.g., TODO) an incompressible model is desired.
  //
  int    incompressible;

  // Param: simple_diffusion
  // Type: int Comments: The simple_diffusion flag toggles a simplified
  // diffusion computation. It simplifies the computation of feq for the solute
  // component. It is based on the observation that the process of diffusion
  // requires less symmetry than flow (c.f. TODO: cite).
  //
  int    simple_diffusion;

  // Param: rho_A[subs]
  // Param: rho_B[subs]
  // Type: double*
  // Comments: rho_A[subs] and rho_B[subs] are the initial density values. The
  // usage of these values varies with respect to the chosen initial condition
  // (see the initial_condition parameter below). If NUM_FLUID_COMPONENTS==1,
  // then rho_A[0] and/or rho_B[0] are used to initialize the fluid density.
  // If NUM_FLUID_COMPONENTS==2, then rho_A[0] and/or rho_B[0] are used to
  // initialize density of fluid component 0 and rho_A[1] and/or rho_B[1] are
  // used to initialize density of fluid component 1. All the values are not
  // necessarily required for all the initial conditions. See the documentation
  // related to the initial_condition parameter for specifics.
  //
  double rho_A[    NUM_FLUID_COMPONENTS];
  double rho_B[    NUM_FLUID_COMPONENTS];

#if INAMURO_SIGMA_COMPONENT
  // Initial/boundary condition values for solute component.
  double rho_sigma;
  double rho_sigma_in;
  double rho_sigma_out;
  double u_sigma;
  double u_sigma_in;
  double u_sigma_out;

  // Param: sigma_start;
  // Type: int    sigma_start;
  // Comments: Timestep at which solute is activated. If sigma_start == 0, solute will be activated from the beginning. The reason for this option is to make it easy to allow the fluid to equilibrate before introducing solute into the flow.
  //
  int    sigma_start;

  // Param: sigma_stop;
  // Type: int    sigma_stop;
  // Comments: Timestep at which solute is deactivated. If sigma_stop < 0 then the solute will never stop. WARNING: if 0 < sigma_stop < sigma_start, solute will never start. This option is useful, for instance, if it is desired to observe how a solute washes out of the domain.
  //
  int    sigma_stop;

  // Param: sigma_btc_rate;
  // Type: int    sigma_btc_rate;
  // Comments: Rate at which to accumulate temporal breakthrough curve. If sigma_btc_rate <= 0, no btc is stored.
  //
  int    sigma_btc_rate;

  // Param: sigma_btc_spot;
  // Type: int    sigma_btc_spot;
  // Comments: The spot in the domain (as distance from inflow) for measuring btc. If sigma_btc_spot is greater than the length of the domain or <0 then the btc will be measured at the outflow.
  //