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.

__init__

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__`.

bind

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:

op_template

flux

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.

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.

op_template

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:

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)
InverseMassOperator()

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.

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.

flux

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:

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.