Boolean 7 - velocity (FHP-III) and interactions Copyright (C) 1997 D.H. Rothman and S. Zaleski. Covered by the GNU General Public License, see the file ../README. This model corresponds to the main program Fhp7_msite.c and a number of subroutines. You are invited to report any bugs or changes required to get this to run on your system to John Lapeyre lapeyre@physics.arizona.edu. 0. How to Use To compile the code type 'make clean' then 'make'. You may want to edit the Makefile first. You can set parameters for experiments using preprocessor #define statements as described below. But it may be easier to use command line arguments (command line arguments will override parameters set with the preprocessor.) For example: Fhp7_msite -x 64 -y 16 -t 100000 -p 2000 -f 0.00001 -r 1 -s 10001 -C -P -S \ --fnprofile=./output/profile.q.1 --solidflag=1 --separate_malloc=1 To make using the arguments easier, a perl script 'dolg' is provided that sets the command line arguments for you. You can edit dolg quickly to change parameters and don't need to recompile the program. The parameters and arguments are described in comments in the executable file 'dolg'. If you want set parameters via the preprocessor, you should first edit the lines in the included file "parameters_7.h", before you compile the program. #define NX 64 #define NY 64 #define TPRINT 1000 #define TMAX 3000 #define NOSEED 137 /* Density is reduced density (nr per channel) */ #define DENSITY 0.5 #define SOLID_SITES #define TSTART_AVER 1000 I. Implementation details, non-interacting code This code uses boolean operations to perform the collisions. The presence of a particle of a given velocity is coded in a single bit. A 32-bit word encodes velocities on a line segment of 32 lattice points oriented in the x-direction. Then propagating the particles to the left, for instance, consists of left shifting the 32 bit word and carrying one bit to the next word to the left. As a consequence, the user is required to choose the number of sites in the x-direction, nx, to be an integral multiple of the word length (usually 32). The code implements the 7 velocity (FHP-III) model with periodic boundary conditions along each axis where it meets and edge. To make propagation faster, sites are stored using a diagonal pattern: site at (x,y) = (i + j/2, \sqrt 3 / 2 j ) are stored in location "pointer" + i + NX*j . In other words, i,j correspond to the non-cartesian axes. Indices start at zero, e.g. 0 <= i < NX . / j / / / / / / ----------------------->i This implies that the periodicity of the lattice is a bit unusual. In fact the following pattern is reproduced periodically -------- -------- ... / / / / / / / / / / / / / * / * / ---------------- .... / / / / / / / / / / / / / * / * / ----------------- ..... Any structure in this pattern (a wall, a vortex) such as the point marked * is periodically located in a monohedral, not rectangular, pattern. If nx=ny, the pattern is hexagonal. This does not affect most applications. For instance channel flow is easily simulated using the Solid.c program to initialize the walls. ---------------------------------------- FORCING_RATE The type of forcing used in Fhp7_msite is a bit more sophisticated than the forcing in Fhp6_simp. First we subdivide the NX x NY lattice in blocks of size word_length x block_length bits, which are fixed parameters. There are then two possibilities, a high forcing rate and a low forcing rate. A high forcing rate occurs when FORCING_RATE*word_len*block_len is larger than 1. Then we want at least one forcing in each block in each time step. We then take the integer part force_push = (integer part of) FORCING_RATE*word_len*block_len and we perform force_push flips of an A particle into a D particle. To perform the flips, we select at random a site in a given block. If the flip is possible, that is there is an A and there is no D, we perform the flip. Otherwise, we try the next site in the block, possibly wrapping around in the block. When the forcing rate is low, then for each block we draw a random number to decide whether we push a particle in the block. Once it is decided that a particle should be pushed, the site is selected at random as above. In both procedures, whenever we fail to find a site to flip, a message is printed. This error message often occurs when the forcing rate is too high. A forcing rate too high results in all particles packed in one direction, with none left to flip. COUNTING and PARABOLAS parameters If you uncomment the lines #define COUNTING #define PARABOLAS the code will perform experiments setting a Poiseuille flow in the channel created by Solid. The code then automatically creates a directory with names related to the intensity of the forcing. In this form the code is rather slow. Because it needs to record accurately the momentum at each time step, it performs a lot of additional calculations. SOLID_SITES If you do not wish solid sites, it may accelerate the code then comment out the line: #define SOLID_SITES II. Interactions Two possibilities existing for interacting liquid gas models. These have not been tested with the new command line arguments and other changes. #define INTERACTIONS then one may define #define INTERACTIONS_ISO or leave the code to use #define INTERACTIONS_CRAS Details are given in the file parameters_7.h III. More details of implementation. The random number generator used is by Knuth. There are two files taken from his site that I renamed rng.c.old.knuth and rng.c.new.knuth. The new file reduces correlations in a very particular circumstance that will not occur with our code. I am currently using his older routine. Knuth prohibits distribution of his code if it has been modified. In order to make my changes available, I distribute his code and a patch that can be applied with the patch program. The main Makefile should handle this by calling the Makefile in the subdirectory rng. You can contact me ( John Lapeyre lapeyre@physics.arizona.edu) if you have some problem with this. Use mkdist to make a tar file of this distribution. Read notes in mkdist.