# Ini-File Framework

Contents

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.