Dealing with Operators

At some point the examples might not serve your purposes any more. Then you need to modify or to create an entire new operator for your problem. If this happens the following section will help you to understand what happens by calling the different methods of the operator and how the final operator which serves as the RHS of the equation will be assembled.


Here all internal used variables of the class should be defined. By calling the specific operator class they will be provided. Read more about `__init__`.


Putting everything together. The bind method returns the method rhs as an operator which can be used on the discretized domain. It serves as an assembly method bringing all parts of the operator together.

The different parts are:

Stepwise through the method bind:

1. Going into the op_template method to create an operator template which only needs to be called with the specific input data (Source, BC, etc.).

compiled_op_template = discr.compile(self.op_template())

2. Defining the vector-field for the state and the boundary conditions.

3. Building the flux operator with the flux method.

 flux_op = get_flux_operator(self.flux())

4. Building all other requested operators (Nabla, Mass-Matrix, Stiffness-Matrix, Inverse Mass-Matrix, etc.).

5. Putting all separate operators together.

6. Returning the compiled operator template.

7. Building a rhs method by feeding the compiled operator with source functions and boundary values.

def rhs(t, w):
    return rhs

The entire building process works with place holders or empty vector-fields and not with the actual vector-fields containing the data of the entire mesh.

Step 1. - 6. will be described in more detail in section op_template. Step 7 will be described in the next section.

rhs - Building the Right Hand Side

After the operator template has been assembled the link between the place holders and the actual functions behind them needs to be provided. This happens in the part where the method rhs is defined:

def rhs(t, w):
    from import join_fields, ptwise_dot
    rhs = compiled_op_template(w=w, dir_bc_u=self.dirichlet_bc_f(t))

    if self.source_f is not None:
        rhs[0] += self.source_f(t)

        return rhs

Here the keywords w and dir_bc_u which served as place holders in the operator template got linked to the input variables w and the boundary function dirichlet_bc_f(t) which has been passed to the operator during initialisation. As a special feature the boundary function uses the time t as input. Depending on how the function looks like it also could have been w or parts of it - w[0] = u (state) oder w[1:] = v (velocity). However this is the crucial part of the implementation of an operator. All new functions used in the operator have to be linked with the field-vectors at this point.


The operator template will not use any vector-fields containing data of the entire mesh but place holders which have to be linked with the actual vector-field.

As an example the StrongWaveOperator will serve.

First the structure of the state needs to be defined:

w = make_vector_field("w", dimensions+1)
u = w[0]
v = w[1:]

The structure of the Operator is a vector-field. Later u and v will be linked to the actual vector-fields containing the value for each node. At this stage the vector-field appears as

[w[0] w[1] w[2]]

when printing it to the screen.

In the next step the vector-fields containing the information about the boundary conditions will be defined. Again only the structure gets defined here and later the actual values get linked to it. The vector-field for the Dirichlet BC's is defined as:

dir_bc = join_fields(Field("dir_bc_u"), v)

In the case of the StrongWaveOperator only the state u will be defined as a BC place holder but not the velocity v. Later dir_bc_u will be linked to an external function calculating the state u on the boundary node. The BC-vectro-field appears as:

[dir_bc_u w[1] w[2]]

at this stage.

As the StrongWaveOperator has different BC's included another two possibilities can be found:

neu_bc = join_fields(u, -v)
rad_bc = join_fields(
         0.5*(u - self.sign*, v)),
         0.5*normal*(, v) - self.sign*u)

neu_bc describes the Neumann-BC's.

SUMMARY: The vector-field for the operator and the vector-field for the BC's have been defined. Both have the same structure [u, v] and [u_bc, v_bc].

Now that we have the framework for the state vector-fields we need to define the flux operator.

flux_op = get_flux_operator(self.flux())

How the flux operator will be built is described in the section flux

The next step is to build the DG specific operators (Mass-Matrix, Stiffness-Matrix, etc.). In the case of the StrongWaveOperator we only need the Nabla-Operator and the inverse Massmatrix

nabla = make_nabla(d)

As both the nabla and the InverseMassOperator are external features from the optemplate module they have to be imported from this module at first:

from hedge.optemplate import make_nabla, InverseMassOperator

SUMMARY: The operator template now has all important parts - vector-fields, flux-operator, DG-operators (Mass-Matrix, Stiffness-Matrix, etc.) to assemble the final expression.

The last step is the assembly of the operator template which gets passed back to the bind method. In case of the StrongWaveOperator the returned expression is a sum of vector fields for the field and the flux.

from import join_fields
    return (
            - join_fields(
                -self.c*, v), 
            InverseMassOperator() * (
                + flux_op * pair_with_boundary(w, dir_bc, self.dirichlet_tag)
                + flux_op * pair_with_boundary(w, neu_bc, self.neumann_tag)
                + flux_op * pair_with_boundary(w, rad_bc, self.radiation_tag)

Important to mention is the flux formulation w.r.t. the BC's. Actually the flux is the sum of the internal flux between the elements


and the fluxes at the boundaries

+ flux_op * pair_with_boundary(w, dir_bc, self.dirichlet_tag)
+ flux_op * pair_with_boundary(w, neu_bc, self.neumann_tag)
+ flux_op * pair_with_boundary(w, rad_bc, self.radiation_tag)

The pair_with_boundary method hast three inputs. w describes the volume vector of the internal part of the mesh and dir_bc describes a boundary vector. The third argument is the boundary tag. If a face has a self.dirichlet_tag as BC then the flux can be calculated by using w as the internal part and dir_bc as the external part. As the StrongWaveOperator has three different BC's all three possibilities are added to the flux. If one BC is not activated this part of the flux will be zero and has no influence.


To assemble the flux operator place holders for the different parts of the state vector will be used. The place holders will later be linked to the values of the vector-fields.

w = FluxVectorPlaceholder(1+dim)
u = w[0]
v = w[1:]

With the place holders the different types of fluxes can be described:

from import join_fields
flux_weak = join_fields(, normal),
     u.avg * normal)

if self.flux_type == "central":
elif self.flux_type == "upwind":
    flux_weak -= self.sign*join_fields(
    0.5*(normal *,
    raise ValueError, "invalid flux type '%s'" % self.flux_type

Here again the join_fields method is used to assemble the different components of the flux vector. Finally the flux operator gets returned to the op_template and there it will be used to assemble the entire WaveOpterator.