Main function

This page has not been updated. However almost all explanations are fully relevant.

track.cpp

The main source is in the file track.cpp.
The simulation is performed by calling an instance of the class Csimulator. This is the last function called in the main. This simulation (which uses templates) takes different arguments that are instances of Multivac classes (e.g. SpeedType, MeshType...). These instances (by the way, all the classes use templates) are to be initialized, and that is precisely what is done in the first part of the main. We have for example to choose the type of speed function we want to use, the size of the mesh, which initial curve we want, whether we want to save every curve or just the last one and whether we want to use the regular Level set method, the narrow band Level set method (this will be explained further) or the Fast marching method...

The main is split in two parts: So each new set of values, written in track.cpp, will trigger a new and different simulation (depending on these new arguments).

From now on, the variables that can be modified to create a new simulation will be followed, in the source file, by:
@***** simulation variable ******@.

The easiest way to proceed is to take the main and explain it in a linear way, which is the natural "program order".

Note: A general explanation about the different types of changes that have to be made by the user to initialize the different class instances used by Multivac is given in DOMAIN MESH(which is used as an example). This procedure, explained with DOMAIN MESH is quite the same for SPEED FUNCTION, INITIAL CURVE, LEVEL SET FUNCTION, INITIALIZER, UPDATER and SAVER.

Includes and options
Time
Domain mesh
Speed function
Initial curve
Level set function
Initializer
Updater
Saver
Simulation


INCLUDES AND OPTIONS

00001 /************************
00002  * INCLUDES AND OPTIONS *
00003  ************************/
00004 
00005 // Seldon library and Multivac project provide error management
00006 // using exception handling.  Exceptions that may be raised
00007 // are selected according to debug levels.
00008 #define SELDON_DEBUG_LEVEL_1
00009 #define MULTIVAC_DEBUG_LEVEL_1
00010 
00011 // Define MULTIVAC_REPORT if you want to clock the time needed for
00012 // updates and initializations (results will be displayed on screen).
00013 #define MULTIVAC_REPORT_NO
00014 
00015 #define ERR(x) cout << "Hermes - " #x << endl
00016 #define DISPLAY(x) cout << #x ": " << x << endl
00017 #define DISP(x) cout << #x ": " << x << endl
00018 
00019 // Seldon library: matrices, vectors and lists.
00020 #include "seldon.hxx"
00021 
00022 // Multivac includes.
00023 #include "includes.hxx"
00024 
00025 using namespace Multivac;
	

This is the header of the file track.cpp, before the definition of the main function.
One notices the use of the namespace Multivac and the two inclusions seldon.hxx and includes.hxx. These two header files include functions and classes that will be used in the main.
As for the namespace, Multivac uses two namespaces Multivac and Seldon.
Actually, Seldon is another "little library" that defines different useful tools such as a class vector, a class matrix, a class list... which are used by Multivac.
There are also three macros defined: ERR(x), DISPLAY(x) and DISP(x).
Finally, there are three other items "defined":
#define MULTIVAC_DEBUG_LEVEL_1
#define SELDON_DEBUG_LEVEL_1
#define MULTIVAC_REPORT_NO

#define MULTIVAC_DEBUG_LEVEL_1 & #define SELDON_DEBUG_LEVEL_1

(MULTIVAC_DEBUG_LEVEL_1 implies SELDON_DEBUG_LEVEL_1 so here, it is redundant)
Multivac and Seldon Libraries are provided with different levels of error management (four for each):
MULTIVAC_DEBUG_LEVEL_1, MULTIVAC_DEBUG_LEVEL_2, MULTIVAC_DEBUG_LEVEL_3, MULTIVAC_DEBUG_LEVEL_4
SELDON_DEBUG_LEVEL_1, SELDON_DEBUG_LEVEL_2, SELDON_DEBUG_LEVEL_3, SELDON_DEBUG_LEVEL_4.

The user is then able to use 4 levels of errors handling (more or less strict).
The level 1 is the least "strict" and the level 4 in the most "complete" one. Each level has the same features than the previous levels plus new ones. (level 2 is level 1 plus new functions). In this example, for this track.cpp, the debug level chosen is the number 1.

Let k be an integer between 1 and 4.
if MULTIVAC_DEBUG_LEVEL_k is defined, it automatically defines SELDON_DEBUG_LEVEL_k.
MULTIVAC_DEBUG_LEVEL_k and SELDON_DEBUG_LEVEL_k handle the same types of errors (for example, they may both take care of memory management) but respectively for the MULTIVAC part and the SELDON part of the program.

What are the features specific to each level?

Level 1
The exception handling mode is "off". No "throws" can be used. All exception specifications are disabled.

Level 2
Memory management & Input/Output operations are checked. (example: creating a new variable, reading/writing a file...). If an error of one of these types occurs, then the program catches the exception and processes a return (an error message is also sent).

Level 3
Memory management & Input/Output operations are still handled (cf Level 2)
Indices are also checked. Multivac & Seldon use matrices and vectors and this Level enable the program to check and report any problem of indices (negative or greater than the size).

Level 4
This level checks everything
But for the moment, it does exactly the same thing as Level 3.

#define MULTIVAC_REPORT_NO, as mentioned in the comments is to assess the duration of the updates and initializations and display it on the screen.
#define MULTIVAC_REPORT could have been another possible "define" here, it is used to monitor the whole process: (initMesh, initphiandF, saveatcurrentiteration...). It records any operation of Multivac and prints the corresponding times on the screen.

Note about a programming technique: The definition of the content of each level or of MULTIVAC_REPORT can be found in the files includes.hxx and seldon/seldon.hxx.
Then in the other files, with a #ifdef ... #endif or a #ifndef ... #endif structure, we can decide to tell the compiler to include that or that action or to disregard it depending on wether or not an item has been "defined". That is why defining MULTIVAC_REPORT_NO or MULTIVAC_REPORT for example, will not give the same results. This use of the tools #define, #ifdef, #ifndef ... which are commands treated by the compiler (because of the "#") is a C++ programming tip. When the compiler sees a #ifdef NAME ... #endif for example, it will take into account the define situation,i.e wether or not NAME has been previously defined by the command #define NAME, if so, it will include the code between #idef NAME and #endif in the compilation, ortherwise it won't. The program can then have different features just by changing the tools defined.


TIME

00028  /*****************
00029   * MAIN FUNCTION *
00030   *****************/
00031 
00032 int main()
00033 {
00034 
00035   clock_t timer = clock();
00036 
00037   // To catch exceptions.
00038   TRY
00039 
00040     // real type: double, float...
00041     typedef double real;
00042 
00043 
00045   // TIME //
00047   
00048   real FinalTime = 0.1;	@***** simulation variable ******@	
00049   real Delta_t = 0.0001;	@***** simulation variable ******@
00050   //  FinalTime += Delta_t;
00051   // Number of iterations.
00052   int NbIterations = int (FinalTime / Delta_t);
	

We can first notice the timer which will enable us to clock the program, the typedef (double will be referred to as real) and the TRY to catch the exceptions. Note: In C++ the try syntax is the following: try{...} catch(), but a lighter notation: TRY ... END has been defined in Multivac in the file includes.hxx.
Then FinalTime, which is the time of the simulation and Delta_t the time step have to be set to a new value by the user depending on the simulation he wishes to do. Delta_t is the step used by the numerical schemes. The program will compute the values of our functions (here, it is the function phi, our "Level set function" which has the value 0 on the front, and which is evaluated at each time step). The simulation stops when the "FinalTime" is reached.
So the number of iterations is known when we have the time step and the final time. The user doesn't have to change it.


DOMAIN MESH

00056   // DOMAIN & MESH //
00058 
00059   // Choose the type of the mesh:
00060   //   1) COrthogonalMesh<real>
00061   //   2) --
00062   typedef COrthogonalMesh<real> MeshType;
00063 
00064   // Domain bounds (the domain is a rectangle).
00065   real Xmin = 0.0;	@***** simulation variable ******@	
00066   real Xmax = 3.0;	@***** simulation variable ******@	
00067   real Ymin = 0.0;	@***** simulation variable ******@	
00068   real Ymax = 3.0;	@***** simulation variable ******@	
00069 
00070   // For an orthogonal mesh, Nx and Ny are the number of
00071   // grid points along (x'x) and (y'y) (respectively).
00072   int Nx = 601;		@***** simulation variable ******@	
00073   int Ny = Nx;		@***** simulation variable ******@	
00074 
00075   // Constructor for an orthogonal mesh.
00076   MeshType Mesh(Xmin, Xmax, Ymin, Ymax, Nx, Ny);
	

This first part of the initialization aims at creating the mesh. For the moment, only one type of mesh can be used: an orthogonal one. Xmin, Xmax, Ymin, Ymax which defines the domain and Nx, Ny (number of grid points) have to be defined by the user.

The structure "typedef COrthogonalMesh<real> MeshType;" is used here so that MeshType will represent, in the rest of the file, COrthogonalMesh<real>. If another mesh type called, let say, CTriangleMesh<real> had been available for example, we would have written "typedef CTriangleMesh<real> MeshType;" to actually use the triangle mesh instead of "typedef COrthogonalMesh<real> MeshType;". The use of this typedef is really convenient because in the rest of the document, the mesh type class is referred to as MeshType (and stands for COrthogonalMesh or CTriangleMesh in our example) and we only have to change one line of the code to select which one we need and change the mesh type.
Then we create a variable Mesh of this class Meshtype (MeshType Mesh(Xmin, Xmax, Ymin, Ymax, Nx, Ny)) with a set of arguments. Indeed, the mesh is stored and saved in an instance, an object of the class MeshType.

General Note: This technique will be used for each class we will study further in this document. Each time a typedef is used to associate a specific class to a generic name, name which is used in the rest of the file.
Then we update a set of variables which are given as an argument to a contructor of the class, to create an instance with the properties needed.
In Multivac, most properties and datas are set in specific "Multivac class instances". For example, the mesh is stored in an instance of MeshType (we have seen previously that MeshType refers to one of the mesh classes because there could be different meshes possible, these different mesh classes are in fact sub-classes of a mother mesh class: C++ heritage). The files containing the definitions of these classes and their sub-classes can be found in the directory includes of Multivac and have explicit names to enable the user to find his way in the program easily.
Then for each mother class, let say for example the speed function class (cf next section), the user just has to select which version (i.e which sub-class) he wants to use, in this precise case CConstantSpeed, CSpeedFunction1 or CFireModel and write it in the typedef (typedef sub-class Speedtype). Then he has to initialize the different variables given as an argument to the instance of the class created. These arguments can be different(most of the time, their number is different) from one sub-class to another.
Finally, depending on which sub-class is used, the user has to give the constructor of the class the corresponding set of arguments. In our example of the speed function, we have to choose between:
1) F(SpeedRate) if the speed type is CConstantSpeed
2) F(SpeedRate, SpeedRate0, Limit) if the speed type is CSpeedFunction1
3) F(U, a, b, epsilon0, epsilon1, c1, mu, CFL, coeff) if the speed type is CFireModel
If the third one is chosen, "SpeedType F(U, a, b, epsilon0, epsilon1, c1, mu, CFL, coeff)" creates an instance named "F" of SpeedType (which stands for CFireModel here) with all the paramaters it needs (9 here).
If the first one, a constant speed had been chosen we have had:
typedef CConstantSpeed SpeedType;
initialization of SpeedRate
SpeedType F(SpeedRate);
Of course, in this case we just need one argument SpeedRate so we don't have to change or write the others (unused). Each time, the user just has to initialize the parameters needed (however in this version of track.cpp, all parameters for every cases are written to help the user, we then just need to modify the parameters we want). Moreover, the different classes and sub-classes available are written in the source code as comments. Note: We will discuss the nature of these "speed function" parameters in the next section.


SPEED FUNCTION

00080   // SPEED FUNCTION //
00082   
00083   // Choose the speed function:
00084   //   1) CConstantSpeed<real>
00085   //   2) CSpeedFunction1<real>
00086   //   3) CFireModel<real>
00087   typedef CFireModel<real> SpeedType;
00088 
00089   // Speed rate for a constant speed function.
00090   real SpeedRate = 0.5;       @***** simulation variable ******@	
00091 
00092   // Second speed rate.
00093   real SpeedRate0 = 0.2;	    @***** simulation variable ******@	
00094 
00095   // Limit (for CSpeedFunction1).
00096   real Limit = 0.9;	@***** simulation variable ******@	
00097 
00098   // Fire model parameters:
00099   // v_f(x) = epsilon0 + c1 * sqrt(x).  <-- head
00100   // beta(x) = epsilon0 + a * x * exp(- b * x).  <-- flank
00101   // epsilon(x) = epsilon0 * exp(- epsilon1 * x).  <-- rear
00102   real U = 100.0;				@***** simulation variable ******@	
00103   real a = 0.1;				@***** simulation variable ******@	
00104   real b = 1.0;				@***** simulation variable ******@	
00105   real epsilon0 = 0.02;			@***** simulation variable ******@	
00106   real epsilon1 = 0.003;			@***** simulation variable ******@	
00107   real c1 = 0.5;				@***** simulation variable ******@	
00108   real mu = 2.0;				@***** simulation variable ******@	
00109 
00110   real CFL = FinalTime / NbIterations
00111   / min( (Xmax - Xmin) / (Nx-1), (Ymax - Ymin) / (Ny-1));
00112   real coeff = 2.0;			@***** simulation variable ******@	
00113   real MaxGrad = 0.5 * sqrt( pow(Xmax - Xmin, 2.0)
00114                              + pow(Ymax - Ymin, 2.0) );
00115 
00116   // Choose the right constructor:
00117   //   1) F(SpeedRate)
00118   //   2) F(SpeedRate, SpeedRate0, Limit)
00119   //   3) F(U, a, b, epsilon0, epsilon1, c1, mu, CFL, coeff)
00120   CFL = 0.9;					@***** simulation variable ******@	
00121   SpeedType F(U, a, b, epsilon0, epsilon1, c1, mu, CFL, coeff);
	

This part deals with the speed function, (the speed is known for every point of the curve).
Programming note: One more time, the speed values (on the grid) are saved in an instance of the speed class. But before being saved, these values have to be assessed depending on which speed model we want to use. This job is done thanks to heritage, virtual functions and sub-classes of the speed class.
There are, as mentionned in the previous section, 3 possible sub-classes defined for the moment:
Note: There are in fact 4 fire models defined, plus the possibility to rotate the wind speed along with the time (the wind speed is a vector, with a norm equal to the parameter U, which has otherwise a constant direction and is then rotated of an angle of 10 x current time). As seen in the previous section, the use of #ifdef enable us here to define several models and select just one during the compilation using the key word #define NAME where NAME belongs to:
MULTIVAC_FIREMODEL4, MULTIVAC_FIREMODEL3, MULTIVAC_FIREMODEL2 or MULTIVAC_FIREMODEL1
The same thing can be applied to MULTIVAC_ROTATING_WIND0 (constant wind direction), MULTIVAC_ROTATING_WIND (rotating wind).
But this time, in this version of the code, #define MULTIVAC_FIREMODEL4 and #define MULTIVAC_ROTATING_WIND0 are called in the head of the file /userpath/includes/speedfunctions/firemodel.cxx. So if we want to try to change it, we have to modify that in this file. These different models are defined in the end of the file /userpath/includes/speedfunctions/firemodel.cxx.

More information can be found on this subject on the website of Vivien Mallet:
http://spacetown.free.fr/fronts/
For any point M, let theta be the angle between the direction of the wind and the normal of the curve at this point M. Let cos2 = cos(theta)*cos(theta) and sin2 = sin(theta)*sin(theta).

Speed function : F

MULTIVAC_FIREMODEL1
F(theta) = epsilon0 * sin2 + a * U * sin2 Exp(-b * U * sin2) + G(theta)
where G(theta) is defined as:
If |theta|<pi/2 then G(theta) = epsilon0 * cos2 + c1 * srqt(U) * costheta^(5/2)
else G(theta) = epsilon0 * cos2 * Exp(-epsilon1 * U * cos2)

MULTIVAC_FIREMODEL2
F(theta) = epsilon0 * sin2 + a * U * sin2 Exp(-b * U * sin2) + G(theta)
where G(theta) is defined as:
If |theta|<pi/2 then G(theta) = epsilon0 * cos2 + c1 * srqt(U) * costheta^(3/2)
else G(theta) = epsilon0 * cos2 * Exp(-epsilon1 * U * cos2)

MULTIVAC_FIREMODEL3
F(theta) = epsilon0 * sin2 + a * U * sin2 Exp(-b * U * sin2) + G(theta)
where G(theta) is defined as:
If |theta|<pi/2 then G(theta) = epsilon0 * cos2 + c1 * srqt(U) * costheta^(5/2)
else G(theta) = epsilon0 * cos2 * Exp(-epsilon1 * U * cos2)

MULTIVAC_FIREMODEL4
Model used here.
if (costheta > 0)
F(theta) = epsilon0 + c1 * sqrt(U) * pow(costheta, 1.5) + a * U * sin2 * exp(-b * U * sin2) );
else
F(theta) = epsilon0 * sin2 + a * U * sin2 * exp(-b * U * sin2)+ epsilon0 * cos2 * exp(-epsilon1 * U * cos2) );


These models calculate the speed knowing the wind speed, a few parameters and the normal of the curve.
With: Note: In the previous versions of Multivac, the speed rate was defined as followed:
If |theta| < pi / 2
F(U, theta) = vm_f(U cos(theta)) + beta(U sin^mu(theta))
else
F(U, theta) = beta(U sin^mu(theta)) + epsilon(U cos^2(theta)).

This definition is written in the function "Value_Normal" in the file "userpath/includes/speedfunction/firemodel.cxx". But this function is not called anymore in the new versions, instead, we use the 4 models previously defined (which are written in the function: "Model" in the same file).
That's why in the new version, the parameter mu is not used anymore (it is included indirectly in the different models).


INITIAL CURVE

00125   // INITIAL CURVE //
00127   
00128   // Choose the initial curve:
00129   //   1) CCircle<real>
00130   //   2) CTwoCircles<real>
00131   //   3) CThreeCircles<real>
00132   //   4) CIsland<real>
00133   //   5) CIsland0<real>
00134   //   6) CSetOfPoints<real>
00135   typedef CCircle<real> InitialCurveType;
00136 
00137   // For a circle, set center coordinates and radius.
00138   real CircleCenterX = 1.0;		@***** simulation variable ******@	
00139   real CircleCenterY = 1.5;		@***** simulation variable ******@	
00140   real CircleRadius = 0.5;		@***** simulation variable ******@	
00141 
00142   // For a second circle, set center coordinates and radius.
00143   real CircleCenterX0 = 1.65;		@***** simulation variable ******@	
00144   real CircleCenterY0 = 1.6;		@***** simulation variable ******@	
00145   real CircleRadius0 = 0.3;		@***** simulation variable ******@	
00146 
00147   // For a third circle, set center coordinates and radius.
00148   real CircleCenterX1 = 0.50;		@***** simulation variable ******@	
00149   real CircleCenterY1 = 1.0;		@***** simulation variable ******@	
00150   real CircleRadius1 = 0.25;		@***** simulation variable ******@	
00151 
00152   // File containing a front defined by a set of points.
00153   // Available shapes in F:/Fronts/multivac/includes/initialcurves/:
00154   // M.pts
00155   string InitialFrontFile = "F:/Fronts/multivac/includes/initialcurves/M.pts";
@***** simulation variable ******@	
00156 
00157   // Choose the right constructor:
00158   //   1) InitialCurve(CircleCenterX, CircleCenterY, CircleRadius)
00159   //   2) InitialCurve(CircleCenterX, CircleCenterY, CircleRadius,
00160   //                   CircleCenterX0, CircleCenterY0, CircleRadius0)
00161   //   3) InitialCurve(CircleCenterX, CircleCenterY, CircleRadius,
00162   //                   CircleCenterX0, CircleCenterY0, CircleRadius0,
00163   //                   CircleCenterX1, CircleCenterY1, CircleRadius1)
00164   //   4) InitialCurve(CircleCenterX, CircleCenterY, CircleRadius,
00165   //                   CircleCenterX0, CircleCenterY0, CircleRadius0)
00166   //   6) InitialCurve(InitialFrontFile, 0)
00167   InitialCurveType InitialCurve(CircleCenterX, CircleCenterY, CircleRadius);
      

The notations are explicit. This class creates the first curve (from a set of datas).
Example: The file M.pts creates the shape of the letter M.

M.pts

0.51 0.2
0.5 0.21
0.5 2.79
0.51 2.8
0.9 2.8
1.49 2.0
1.51 2.0
2.1 2.8
2.49 2.8
2.5 2.79
2.5 0.21
2.49 0.2
2.11 0.2
2.1 0.21
2.1 1.6
2.09 1.61
1.51 1.2
1.49 1.2
0.91 1.61
0.9 1.6
0.9 0.21
0.89 0.2
	


Initial M curve:

M.pts

LEVEL SET FUNCTION

00171   // LEVEL SET FUNCTION //
00173   
00174   // Choose the level set function type:
00175   //   1) COrthogonalLevelSet<real>
00176   //   2) --
00177   typedef COrthogonalLevelSet<real> LevelSetType;
00178 
00179   // Constructor for a level set function defined
00180   // on an orthogonal mesh.
00181   LevelSetType Phi;
	

We only have one choice, we have to create phi on an orthogonal mesh because we have seen previously that it is the only type of mesh for the moment. By the way, creating a 3D simulation or using a mesh using triangles could be future improvements.
So, there is nothing to be done with this class.


INITIALIZER

00185   // INITIALIZER //
00187   
00188   // Choose an initializer:
00189   //   1) CNeverInit<real> *
00190   //   2) CNarrowBandNeverInit<real> **
00191   //   2) CNarrowBandInitDomain<real> **
00192   //   4) CFastMarchingNeverInit<real> ***
00193   //         * Level set method
00194   //        ** Narrow band level set method
00195   //       *** Fast marching method
00196   typedef CNarrowBandNeverInit<real> InitializerType;
00197 
00198   InitializerType Initializer;
	

From now on, a new feature is now introduced: the method used to realize the simulation.
In the previous initializations, it had no inluence.
The Level Set Method and the Fast Marching Method are both implemented under Multivac.
By the way, the Level Set Method is implemented "twice". In fact, a more efficient version of the Level Set Method has been implemented as well: this is the Narrow Band Level Set Method. Phi is computed within a narrow band instead of being assessed on the whole grid. That is why this method, developped for programming, is better than the regular algorithm. In Multivac, the term NarrowBand is used in the name to refer to this method. We also refer to the Fast Marching Method by using FastMarching in the names as well.
More information about these different techniques can be found in the book of J.A Sethian Level Set Methods and Fast Marching Methods Cambridge University press.

In the initializer, if we want to use the regular Level Set Method or the Fast Marching Method we have to use respectively the first (CNeverInit) and the fourth class (CFastMarchingNeverInit).
These initializers only provide first initializations. It doesn't provide any reinitialization or update.
The mesh is initialized. The level set is initialized according to the initial curve position. Then, the speed function is initialized according to the initial level set function. The reinitialization only calls the speed function updater.

But if we use the narrow band version, then we have two choices: "never init" and "init domain".
This initializer can be used to perform initializations and reinitializations of the level set function and the tube involved in the narrow band level set method.
First, the mesh is initialized. The level set is initialized according to the initial curve position. Then, the speed function is initialized according to the initial level set function. Reinitializations build the tube and init the corresponding level set function. As for the speed function, reinitializations only call the speed function updater.
In the narrow band method, the tube has to be recalculated any time a point of the front enters a limit zone defined in the updater (cf next section). In fact, with this method, the front evolves in a certain tube and sometimes tends to go out of it. In such situations, a new tube is recomputed with its associated phi (this is handled by the initializer).
The only difference between "never init" and "init domain" is that "init domain", reinitializes a new domain (newXmin, newXmax, newYmin, newYmax) from the initial domain (Xmin, Xmax, Ymin, Ymax) anytime to fit the current situation. Indeed, we only work with a narrow band, so shrinking the domain to a smaller one (of course still containing the curve and the band) and using the same number of grid points (Nx x Ny) is more accurate since the new space steps (newDelta_x & newDelta_y) are smaller.
For "never init", the initial domain and the initial grid pattern are kept.


UPDATER

00202   // UPDATER //
00204 
00205   // Choose an initializer:
00206   //   1) CFirstOrderEngquistOsher<real> *
00207   //   2) CFirstOrderLaxFriedrichs<real> *
00208   //   3) CFirstOrderGodunov<real> *
00209   //   4) CNarrowBandFirstOrderEngquistOsher<real> **
00210   //   5) CNarrowBandFirstOrderLaxFriedrichs<real> **
00211   //   6) CNarrowBandFirstOrderGodunov<real> **
00212   //   7) CFastMarchingFirstOrderEngquistOsher<real> ***
00213   //         * Level set method
00214   //        ** Narrow band level set method
00215   //       *** Fast marching method
00216   typedef CNarrowBandFirstOrderEngquistOsher<real> UpdaterType;
00217 
00218   // Choose the right constructor and choose its parameters:
00219   //   1, 2, 3) Updater
00220   //              | Default constructor is enough for 1).
00221   //   4, 5, 6) Updater(TubeSemiWidth, BarrierWidth, OutSpaceWidth)
00222   //              | TubeSemiWidth: number of grid points on each side of the front.
00223   //              | BarrierWidth: number of grid points (on each side of the front)
00224   //                        that imply tube reconstruction when reached.
00225   //              | OutSpaceWidth: number of grid points (on each side of the front)
00226   //                         that must not be reached.
00227   //              | Example: Updater(6, 3, 1) or Updater(12, 5, 1).
00228   //   7) Updater(TMax)
00229   //        | TMax: time greater than all arrival times that will be computed.
00230   UpdaterType Updater(6, 3, 1);		@***** simulation variable ******@	
	


The updater is the class that handles the numerical scheme and "updates" phi (i.e: evaluates phi at the time t+dt knowing phi at the time t). It also handles the tube in the Narrow Band version (checks wether the tube's "limit barrier" is reached or not) and the binary tree in the Fast Marching Method.
We have to choose our numerical scheme between EngquistOsher, LaxFriedrichs and Godunov. These are all first order schemes. Ones again, if we use the regular Level Set Method we have to choose among the three first: 1) 2) or 3), if we use the Narrow Band version we can use 4) 5) or 6) and we only have one possibilty for the Fast Marching Method 7).
As for the arguments, the regular Level Set Method takes none and the Fast Marching Method only needs "a time greater than all arrival times that will be computed".
But for the Narrow Band Algorithm, three arguments are required: TubeSemiWidth, BarrierWidth, OutSpaceWidth.
As mentionned as a comment in the code:



Here is a scheme of the situation: (the bold line is the front)

narrow band


The user has to choose these values carefully. If the tube is too large, the computation time increases. But if it is too small, then it has to be reinitialized a lot, because the curve reaches the second zone too often, and the cost of such operation is quite high.
In this reason the choice of these parameters is an important issue, as far as computational time is concerned. A reconstruction of the tube every 6 to 10 iterations seems to be a fine compromise and gives good results. That's why the parameters should be chosen (knowing the grid and having an idea of the speed) so that this range (6-10 iterations) could be possible. It would give the best results.


SAVER

00234   // SAVER //
00236 
00237   // Choose the saver type:
00238   //   1) CNeverSave<real>
00239   //        | Nothing is saved.
00240   //   2) CCurvesSaver<real>
00241   //        | The front is constructed and saved at each time step.
00242   //        | Not relevant for the fast marching method.
00243   //   3) CSaveLastCurve<real>
00244   //        | Saves the last curve, after all calculations.
00245   //        | Not relevant for the fast marching method.
00246   //   4) CSaveAtTheEnd<real>
00247   //        | Saves the level set, after all calculations.
00248   //        | Relevant for the fast marching method.
00249   //   5) CTrack<real>
00250   //        | Relevant for optimization.
00251   typedef CTrack<real> SaverType;
00252 
00253   // Number of curves that will be saved.
00254   // Set NbCurves to 0 in order to save all curves.
00255   int NbCurves = 10;		@***** simulation variable ******@
00256 
00257   // flag identifies data type (double: 'f', ... -- cf fwrite from stdio.h --)
00258   char flag = 'f';
00259   // 
00260   int width = 10;	@***** simulation variable ******@
00261   // 
00262   int precision = 10;	@***** simulation variable ******@
00263 
00264   // Directory where files will be created.
00265   string Directory = "/home/vivien/Computations/Fronts/"; @***** simulation variable ******@
00266 
00267   // Stores time at each time step.
00268   string TimeFile = Directory + "Time";
00269 
00270   // Stores front points.
00271   string CurvesFile = Directory + "Curves";
00272   // Stores number of points on fronts at each time step.
00273   string CurveLengthsFile = Directory + "CurveLengths";
00274   // Stores the level set function.
00275   string PhiFile = Directory + "Phi";
00276   // Stores the speed function.
00277   string FFile = Directory + "F";
00278 
00279   // Stores points coordinates along the (x'x) axe.
00280   string XFile = Directory + "X";
00281   // Stores points coordinates along the (y'y) axe.
00282   string YFile = Directory + "Y";
00283   // Stores mesh points.
00284   string PointsFile = Directory + "Points";
00285   // Stores mesh edges.
00286   string EdgesFile = Directory + "Edges";
00287   // Stores mesh triangles.
00288   string TrianglesFile = Directory + "Triangles";
00289 
00290   // Indicates the number of iterations between saves.
00291   int Period = NbIterations / (NbCurves==0?1:NbCurves);
00292 
00293   SaverType Saver(flag, width, precision, TimeFile, CurvesFile,
00294                   CurveLengthsFile, PhiFile, FFile, XFile, YFile,
00295                   PointsFile, EdgesFile, TrianglesFile, Period);
	

As we work with the level set function phi, the front is not always reconstructed and if we want to save the front we sometimes have to build it first. That is what the Saver class does. We can first notice that there is a variable called NbCurves that enable the user to choose how many curves he wants to save. A value 0 means that all the curves will be saved. Otherwise, if we choose to save let say 4 curves, it won't be the first 4 curves, but a curve will be chosen every period (the period is defined as the total number of iterations divided by the number of curves to save). string Directory = "/home/vivien/Computations/Fronts/" has also to be changed and set by the user to the output directory he wants. Example: In the "organisation tip section", the output advised is:
string Directory = "userpath/runs/runk/output/".

The user just has to set the path in which he wants the output files to be stored. These output files have explicit names and contain the different values saved (curves, phi, mesh...).



SIMULATION

To spawn a simulation, an instance called Simulator of the class CSimulator, taking as parameters the different instances of every Multivac class previously initialized (MeshType, SpeedType, Phi...), is created.
Then, Simulator is initialized using the member function Init: "Simulator.Init()" and the simulation is lauched with the member function run: "Simulator.Run()".

Note about the optimization:
Multivac has also an optimization class that we can use instead of the simulation.
The optimization part is still being improved. More informations about the algorithm chosen to compute the optimization can be found on the following website: http://spacetown.free.fr/fronts/.
Both optimization and simulation are called with the same set of parameters. Everything we have seen in this document is unchanged when we use the optimization function.
Coptimizer2nd is the name of the class chosen if we want an optimization.
So, the user just has to select wether he wants to do a simulation or an optimization by creating the instance Simulator from the right class.
We have only two choices.
Note: the structure is the following: the name of the class (e.g. Csimulator), a list of templates used by the class (between <...>), the name of the instance created (Simulator) and the list of arguments.