# 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.

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 hedge.tools 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.

## 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:

neu_bc = join_fields(u, -v) rad_bc = join_fields( 0.5*(u - self.sign*numpy.dot(normal, v)), 0.5*normal*(numpy.dot(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) 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.

from hedge.tools import join_fields return ( - join_fields( -self.c*numpy.dot(nabla, v), -self.c*(nabla*u) ) + InverseMassOperator() * ( flux_op*w + 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

flux_op*w

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.

## 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:

from hedge.tools import join_fields flux_weak = join_fields( numpy.dot(v.avg, normal), u.avg * normal) if self.flux_type == "central": pass elif self.flux_type == "upwind": flux_weak -= self.sign*join_fields( 0.5*(u.int-u.ext), 0.5*(normal * numpy.dot(normal, v.int-v.ext))) else: 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`.