0
Comment:

← Revision 11 as of 20130116 23:23:13 ⇥
11280

Deletions are marked like this.  Additions are marked like this. 
Line 1:  Line 1: 
## page was renamed from SciComp/Hedge/HowTo/IniFileFramework ## page was renamed from SciComp/HedgeHowTo/IniFileFramework = IniFile Framework = <<TableOfContents>> 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 (RungeKutta, AdamsBashforth, 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 disklike 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: <<Anchor(simpmesh)>> {{{ 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. {{attachment:../rectang_mesh_small.png}} == Discretization of the Mesh == As nodal DGFramework 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. <<Anchor(discr)>> {{{ 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 [[#simpmeshmesh]] by: {{{ from hedge.backends.jit import Discretization discr = Discretization(mesh, order=3) }}} will give you a fine mesh with 20 nodes. {{attachment:../rectang_3O_fine_mesh_small.png}} === Example (VerticesCoordinates) === After [[#discrdiscretizing]] the [[#simpmeshmesh]] 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 vectorfield 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 (SourceFunction) === The wavemin.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 [[#simpmeshmesh]]. {{{ [ 8.15987835e04 2.76077257e06 2.76077257e06 1.60381089e28 2.10422364e17 2.10422364e17 1.60381089e28 2.10422364e17 2.10422364e17 1.60381089e28 8.15987835e04 2.76077257e06 2.76077257e06 1.60381089e28 2.10422364e17 2.10422364e17 1.60381089e28 2.10422364e17 2.10422364e17 1.60381089e28] }}} === Example (BoundaryFunction) === 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 xcoordinate 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: * SpaceDimensions (1D, 2D, 3D) * SourceFunctions * BoundaryConditions (Functions) * FluxType (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 `wavemin.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 [[#OperatorBuildBuilding 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 vectorfield of zeros with the length of the volumevector describing the discretization. === Example (2D WaveEquation) === 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 zerovectorfields 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 RungeKutte scheme * AdamsBashfort scheme * MultiRateAdamsBashforth 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 AdamsBashfort 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. 
IniFile 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 (RungeKutta, AdamsBashforth, 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 disklike 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 DGFramework 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 (VerticesCoordinates)
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 vectorfield 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 (SourceFunction)
The wavemin.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.15987835e04 2.76077257e06 2.76077257e06 1.60381089e28 2.10422364e17 2.10422364e17 1.60381089e28 2.10422364e17 2.10422364e17 1.60381089e28 8.15987835e04 2.76077257e06 2.76077257e06 1.60381089e28 2.10422364e17 2.10422364e17 1.60381089e28 2.10422364e17 2.10422364e17 1.60381089e28]
Example (BoundaryFunction)
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 xcoordinate 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:
 SpaceDimensions (1D, 2D, 3D)
 SourceFunctions
 BoundaryConditions (Functions)
 FluxType (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 wavemin.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 vectorfield of zeros with the length of the volumevector describing the discretization.
Example (2D WaveEquation)
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 zerovectorfields 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 RungeKutte scheme
 AdamsBashfort scheme
 MultiRateAdamsBashforth 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 AdamsBashfort 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.