0
Comment:

← Revision 9 as of 20130116 23:25:44 ⇥
8654

Deletions are marked like this.  Additions are marked like this. 
Line 1:  Line 1: 
## page was renamed from SciComp/Hedge/HowTo/DealingwithOperators ## page was renamed from SciComp/HedgeHowTo/DealingwithOperators ## page was renamed from SciComp/HedgeHowTo/BuildingOperators <<Anchor(OperatorBuild)>> = Dealing with Operators = <<TableOfContents>> 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 [[http://docs.python.org/reference/datamodel.html?highlight=__init__#object.__init__`__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 vectorfield 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, MassMatrix, StiffnessMatrix, Inverse MassMatrix, 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 vectorfields and not with the actual vectorfields containing the data of the entire mesh. Step 1.  6. will be described in more detail in section [[#OpTemplateop_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 fieldvectors at this point. <<Anchor(OpTemplate)>> == op_template == The operator template will not use any vectorfields containing data of the entire mesh but place holders which have to be linked with the actual vectorfield. 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 vectorfield. Later `u` and `v` will be linked to the actual vectorfields containing the value for each node. At this stage the vectorfield appears as {{{ [w[0] w[1] w[2]] }}} when printing it to the screen. In the next step the vectorfields 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 vectorfield 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 BCvectrofield 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 NeumannBC's. SUMMARY: The vectorfield for the operator and the vectorfield 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 vectorfields 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 [[#BuildFluxflux]] The next step is to build the DG specific operators (MassMatrix, StiffnessMatrix, etc.). In the case of the `StrongWaveOperator` we only need the NablaOperator 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  vectorfields, fluxoperator, DGoperators (MassMatrix, StiffnessMatrix, 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. <<Anchor(BuildFlux)>> == 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 vectorfields. {{{ 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.intu.ext), 0.5*(normal * numpy.dot(normal, v.intv.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`. 
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 vectorfield 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, MassMatrix, StiffnessMatrix, Inverse MassMatrix, 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 vectorfields and not with the actual vectorfields 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 fieldvectors at this point.
op_template
The operator template will not use any vectorfields containing data of the entire mesh but place holders which have to be linked with the actual vectorfield.
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 vectorfield. Later u and v will be linked to the actual vectorfields containing the value for each node. At this stage the vectorfield appears as
[w[0] w[1] w[2]]
when printing it to the screen.
In the next step the vectorfields 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 vectorfield 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 BCvectrofield 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 NeumannBC's.
SUMMARY: The vectorfield for the operator and the vectorfield 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 vectorfields 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 (MassMatrix, StiffnessMatrix, etc.). In the case of the StrongWaveOperator we only need the NablaOperator 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  vectorfields, fluxoperator, DGoperators (MassMatrix, StiffnessMatrix, 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 vectorfields.
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.intu.ext), 0.5*(normal * numpy.dot(normal, v.intv.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.