Adding a speed function

Introduction

Multivac is designed to handle a wide range of front propagation problems. One may want to provide his own initial front and his own speed function. Many initial fronts may be defined thanks to the class CSetOfPoints. Defining a new velocity is easy (at least with the explanations of this section) but it requires to code a new class.

Class design

All speed functions have to be implemented according to a given interface. This interface is defined by a base class, see includes/speedfunctions/baseclass.hxx. This base class, called CSpeedFunction is then derived in a speed-function class. To define your speed function, you have to create such a derived class.

Here is a description of members of the class CSpeedFunction (open includes/speedfunctions/baseclass.hxx):

• Attribute Matrix<T> Values: a matrix that stores speed rates at mesh points.
• Attributes bool dependence_*: booleans that describe the speed-function dependencies.
• Method Init: is called before the first iteration to initialize the speed function. Usually it allocates Values, if needed.
• Method operator() (T x, T y, T time): is called by the fast marching algorithms to compute the velocity at mesh points and at a given time.
• Method operator() (T x, T y, T time, T nx, T ny, T curvature): is called by the narrow band algorithms to compute the velocity at mesh points, at a given time, for a given normal to the front and a given curvature.
• Methods GetMaxF*(T Xmin, T Xmax, T Ymin, T Ymax, T norm2): returns an estimation of a given quantity that is necessary in the Lax-Friedrichs scheme.
• Method operator() (int i, int j): returns the speed rate at mesh point (i, j).
• Methods Get*Derivatives: not used in currently distributed version of Multivac. Discard them for the moment.
• Other methods are defined in the base class CSpeedFunction and should not be redefined in your derived class.

Defining a new speed function

The point is to define a few methods among those previously described. The best way is to duplicate another speed function and to modify it. Let's try to implement the following speed function:

F_example(x, y, time, nx, ny, curvature) = x^2 + y^2

where 'x' and 'y' are the coordinates in space, 'time' is the current date (the initial time is so that 'time'=0), 'nx' and 'ny' are the two components of the normal to the front and 'curvature' is the curvature of the front.

Files

In includes/speedfunctions/, copy firemodel.hxx into newvelocity.hxx and firemodel.cxx into newvelocity.cxx. Open newvelocity.hxx and newvelocity.cxx.

Put relevant include guards: rename FILE_SPEEDFUNCTIONS_FIREMODEL_CXX into FILE_SPEEDFUNCTIONS_NEWVELOCITY_CXX and FILE_SPEEDFUNCTIONS_FIREMODEL_HXX into FILE_SPEEDFUNCTIONS_NEWVELOCITY_HXX. Replace #include "firemodel.hxx" with #include "newvelocity.hxx". Finally replace CFireModel with CNewVelocity in both files. So, you have a new valid class, but it is still the fire model...

Initialization

First, modify the constructors. You can keep only one constructor. If the default copy-constructor is not enough, you must provide a copy constructor because Multivac needs it. But you should put heavy initializations (in terms of memory or computational requirements) in the method Init. Since Multivac uses the copy constructor (and therefore duplicates the object), the method Init may be a good replacement for the constructor. Modify the destructor too. For our speed function F_example, keep only the default constructor (no argument) where dependence_position should be set to true and other dependencies to false (refer to the expression of F_example, above). Don't forget to update newvelocity.hxx.

Method Init is called by Multivac once, at the very beginning of the simulation (after the initialization of the mesh and of the initial front). Keep at least the reallocation of Values. As for our new velocity F_example, keep Init as it is. No initialization is needed there. That would be different if we had to read data (e.g. topography for fire propagation, an image, ...).

Speed rate

operator() (T x, T y, T time) is needed by the fast marching method and operator() (T x, T y, T time, T nx, T ny, T curvature) is called by the narrow band method. Define the one you need. You should define it carefully because the main computational costs could come from this part of the code.

As for F_example, we could use the fast marching method. However this method cannot be applied to lots of problems. Assuming that you want to use the narrow band method, define operator() (T x, T y, T time, T nx, T ny, T curvature). Just replace "return Model(nx, ny, time);" with "return x*x + y*y;". By the way, you can remove the method Model.

If you want to use the Lax-Friedrichs scheme, you must define GetMaxF1 and GetMaxF2. F_example doesn't depend on the normal to the front, so replace c_1 * m * sqrt(U) with 0.

Update newvelocity.hxx according to your changes and your speed function is ready.

Test with F_example

You can now edit track.cpp (main example provided with Multivac). Add #include "newvelocity.cxx" after #include "multivac.hxx". Replace "typedef CSimplifiedFireModel<real> SpeedType;" with "typedef CNewVelocity<real> SpeedType;". Replace "SpeedType F(U, m, c1, epsilon0, alpha);" with "SpeedType F;" (call to the default constructor).

Compile and run track. The result should be:

Remarks

Of course, it is perfectly possible to call another code from Multivac. For instance, you could call a Fortran code that computes the velocity (in operator() (T x, T y, T time, T nx, T ny, T curvature)).

Multivac uses the C++ library Seldon which provides vectors and matrices. Seldon is provided with Multivac in the directory includes/seldon/. A user's guide for Seldon is available at http://spacetown.free.fr/lib/seldon/. While implementing your speed function, you could need vectors and matrices and using Seldon would make sense.

If your speed function is intricate, you can compile it with strict options: -Wall -ansi -pedantic. This way, g++ issues warnings that may be useful. Multivac doesn't generate any warning with those options.