Ini-File Framework

To run hedge you need to make some fundamental setting like:

• Mesh (1D, 2D, 3D, circle, cylinder, square, cube, etc.)
• Order of discretization
• Boundary Conditions (Dirichlet, Neuman, Inflow, Outflow, etc.)
• Timestepping method (Runge-Kutta, Adams-Bashforth, etc.)
• Total Time
• etc.

All these initial settings will always be done in the same framework. The following sections will give short description of every step in the setup of a computation.

Building the Mesh

The mesh in hedge is provided by the module hedge.mesh with a set of several simple mesh geometries for all kinds of applications in one, two and three dimensions. You can find them in hedge/src/python/mesh.py in the source folder of hedge. To create a 2D disk you have to type:

from hedge.mesh import make_disk_mesh
mesh = make_disk_mesh(r=0.5,faces=100,max_area=0.008)

r is the disk's radius, faces is the number of Faces used to shape the disk and max_area is the maximum area of one element. By changing the number of faces you can decide how disk-like the disk shall be. Three faces will approximate the disk by a triangle. Four faces will create a rectangle. More faces will make the shape of the disk look more like a circle of course. With the max_area you can trigger how many elements will be approximating the disk and how detailed your results will be. A small max_area number will create many elements.

The object mesh now provides all geometrical information about the specified disk. It serves as a base for further discretization.

Example (simple rectangular mesh)

A very basic mesh which will allow us to show some features provided by hedge is the following 2D rectangle:

from hedge.mesh import make_rect_mesh
mesh = make_rect_mesh(a=(-0.5,-0.5),b=(0.5,0.5),max_area=0.8)

You will get a rectangle approximated by two triangles.

Discretization of the Mesh

As nodal DG-Framework needs a different amount of interpolation points - nodes - per element w.r.t. the spatial order, the mesh needs to be discretized. This can be done for a certain order by the Discretization module.

from hedge.backends.jit import Discretization
discr = Discretization(mesh, order=4)

The object discr provides a wide range of methods to gain information about the mesh with the specific discretization. It allows to build functions w.r.t. to the coordinates of the mesh. The following method

discr.interpolate_volume_function(f_u)

will interpolate the function f_u(x,element) - with x[0:dim] - on the volume vector (x,y) of all nodes in the mesh. The result will be a vector (-field) with the values of f_u calculated from every node in the mesh.

Another important feature provided by the object discr is the possibility to evaluate a function f_bc only on boundary nodes of the mesh. The method

discr.interpolate_boundary_function(f_bc)

will give a vector of values for each node on the boundary of the mesh.

Besides the two mentioned methods discr provides a lot of other methods which are important for building the operator or to implement the boundary conditions. The module where all methods are defined can be found in discretization.py.

Example (Discretizing to 3rd order)

Discretizing the before generated mesh by:

from hedge.backends.jit import Discretization
discr = Discretization(mesh, order=3)

will give you a fine mesh with 20 nodes.

Example (Vertices-Coordinates)

After discretizing the mesh you might want to look at the coordinates of the nodes. This can be done by

discr.nodes

which will give you a numpy array of the coordinates.

As a slower alternative you could also use the interpolate_volume_function to create a vector field of the coordinates:

from hedge.tools import join_fields
coords = join_fields(discr.interpolate_volume_function(lambda x,el : x[0]),
                     discr.interpolate_volume_function(lambda x,el : x[1]))

You will receive a vector-field of two vectors showing each the x and y component of a node.

[ [-0.16666667 -0.2236068   0.2236068   0.5   0.2236068  -0.2236068  -0.5  -0.5        -0.5        -0.5   0.16666667  0.2236068  -0.2236068  -0.5   -0.2236068   0.2236068   0.5    0.5         0.5         0.5       ]
  [-0.16666667  0.2236068  -0.2236068  -0.5  -0.5        -0.5        -0.5   0.2236068  -0.2236068   0.5   0.16666667 -0.2236068   0.2236068   0.5    0.5         0.5         0.5   -0.2236068   0.2236068  -0.5       ]]

As the dimensions might vary you might write the statement in a more general way using a for loop:

coords = join_fields([discr.interpolate_volume_function(lambda x,el : x[i]) for i in range(discr.dimensions)])

The result will be the same.

Example (Source-Function)

The wave-min.py example uses source function source_u to provide a state u.

def source_u(x, el):
    return exp(-numpy.dot(x, x)*128)

To get the state from the source function at every point of the mesh you have to type:

source_u_vec = discr.interpolate_volume_function(source_u)

You will receive a vector of values showing the result of the source_u function having been evaluated at every node of the mesh.

[  8.15987835e-04   2.76077257e-06   2.76077257e-06   1.60381089e-28
   2.10422364e-17   2.10422364e-17   1.60381089e-28   2.10422364e-17
   2.10422364e-17   1.60381089e-28   8.15987835e-04   2.76077257e-06
   2.76077257e-06   1.60381089e-28   2.10422364e-17   2.10422364e-17
   1.60381089e-28   2.10422364e-17   2.10422364e-17   1.60381089e-28]

Example (Boundary-Function)

You might want to define a certain state on the boundary nodes of the mesh defined by a function

def bound_u(x,el):
    if x[0] > 0:
        return 1
    else:
        return 0

which sets the state u=0 if the the x-coordinate of the node is <= 0 and u=1 if x > 0. To evaluate this function only on the boundary nodes you can use the interpolate_boundary_function method of the discr object.

bound_u_vec = discr.interpolate_boundary_function(bound_u)

You will receive a vector showing the result only for the boundary nodes

[ 1.  1.  1.  1.  0.  0.  1.  1.  0.  0.  0.  0.  1.  1.  0.  0.]

which are less than all nodes of the mesh. In this example there are 16 boundary nodes. For meshes consisting of more than two elements the difference of the magnitude between boundary nodes and all nodes will be much bigger of course.

Visualization Output Setting

In order to provide an output of the results you can choose between '*.silo' or '*.vtk' output file formats:

from hedge.visualization import VtkVisualizer
vis = VtkVisualizer(discr, None, "fld")

or

from hedge.visualization import SiloVisualizer
vis = SiloVisualizer(discr, None)

Operator Setup

At this point you should know which kind of PDE you want to solve. The operator provides the RHS of the equation which goes into the timestepper. A selection of different operators can be found in pde.py in the source folder. Depending on the problem you might have several inputs for the operator. At least some of the most common inputs might be:

• Space-Dimensions (1D, 2D, 3D)
• Source-Functions
• Boundary-Conditions (Functions)
• Flux-Type (central, upwind, etc.)

The input of an operator depends on the problem and should be looked up in `pde.py'.

Example (Initializing an Operator)

The wave-min.py example uses the StrongWaveOperator. You can define an object op as instance of the class StrongWaveOperator.

    from hedge.pde import StrongWaveOperator
    from hedge.mesh import TAG_ALL, TAG_NONE
    op = StrongWaveOperator(1, discr.dimensions,
            source_vec_getter,
            dirichlet_tag=TAG_NONE,
            dirichlet_bc_f=bound_vec_getter,
            neumann_tag=TAG_NONE,
            radiation_tag=TAG_ALL,
            flux_type="upwind")

To get the RHS operator the method bind has to be used:

rhs = op.bind(discr)

As the bind method actually builds the entire operator an a lot of things a more detailed explanation to this part will be given in section Building Operators

Initial State

To start a computation you will have to define an initial state of the solution. Using the feature

discr.volume_zeros()

of the discr object you can create a vector-field of zeros with the length of the volume-vector describing the discretization.

Example (2D Wave-Equation)

The 2D wave equation has three different arguments for each interpolation point (u, v_x, v_y). To build an initial state you will have to combine three zero-vector-fields created by discr.volume_zeros() with the join_field method from the hedge.tools module. The code would look like this:

from hedge.tools import join_fields
fields = join_fields(discr.volume_zeros(),
                     discr.volume_zeros(),
                     discr.volume_zeros())

You will receive the following vector field:

[[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]]

Timestep setup

After dealing with setup of the spatial part of the computation some temporal settings have to be made.

A method to find the right timestep ensures that the computation will be stable

dt = discr.dt_factor(op.max_eigenvalue()) 

The timestep.py source file provides a module including different stepping methods:

• 4th Order Runge-Kutte scheme
• Adams-Bashfort scheme
• Multi-Rate-Adams-Bashforth scheme

The timestepper can be set by:

from hedge.timestep import RK4TimeStepper
stepper = RK4TimeStepper()

Example (Timestepper)

As a variation to the 4th order RK scheme you can choose the Adams-Bashfort scheme with a certain order by typing:

from hedge.timestep import AdamsBashforthTimeStepper
stepper = AdamsBashforthTimeStepper(3,None)

Timestep Loop

Finally all setup has been done and the computational loop can start. A typical way how this works is shown below:

nsteps = int(700/dt)
for step in xrange(nsteps):
    logmgr.tick()

    t = step*dt
    if step % 5 == 0:
        visf = vis.make_file("fld-%04d" % step)
        vis.add_data(visf, [ ("u", u), ], 
                    time=t, 
                    step=step
                    )
        visf.close()

    u = stepper(u, t, dt, rhs)

The loop includes the output to the screen and to the output files. The field data u gets updated every time and the rhs provides the spatial evolution of the solution.