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.
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
):
Matrix<T> Values
: a
matrix that stores speed rates at mesh points.
bool dependence_*
: booleans
that describe the speed-function dependencies.
Init
: is called before the first
iteration to initialize the speed function. Usually it
allocates Values
, if needed.
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.
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.
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.
operator() (int i, int j)
:
returns the speed rate at mesh point (i, j).
Get*Derivatives
: not used in
currently distributed version of Multivac. Discard them for
the moment.
CSpeedFunction
and should not be redefined in
your derived class.
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.
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...
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, ...).
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.
You can now edit track.cpp
(main example provided
with Multivac). Add #include "newvelocity.cxx"
after #include "multivac.hxx"
. Replace
"typedef CSimplifiedFireModel
" with
"typedef CNewVelocity
". Replace
"SpeedType F(U, m, c1, epsilon0, alpha);
" with
"SpeedType F;
" (call to the default constructor).
Compile and run track
. The result should be:
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.