title: Using the full toolchain

Navigate: [&larr; Set up your first simulation with Musubi](tut_1_mus_config.html) 
| [Overview](index.html)
| [Tracking quantities during the simualtion run &rarr;](tut_3_tracking.html) 

In this tutorial, you will learn how to use the complete tool-chain of the APES suite to deploy
simulations. This includes the generation of the mesh with Seeder, 
Running the actual simulation with Musubi and finally post-process the results with mus_harvesting
which generates files to be visualized with standard tools like Gnuplot and Paraview.
Therefore we use a channel as an example.

## Preparing the Simulation Folder ##

Generate the folders

- mesh
- tracking
- restart
- harvest
- output


## Generate the Mesh ##

Make sure that Seeder is compiled and you know how to call it.

At first, we define the common header of the file containing general information like the name and
the refinement level of the mesh. 
```lua
--! [common header]

level =  8                -- general level of mesh refinement
refinementLevel  = 0      -- special levels for local refinement
refinementLevel2 = 0
label   = 'channel'       
idLabel = 'channel'
simName = 'channel'
outputname = 'channel_cylinder'
outputpreview = true
folder = 'mesh/' 

--! [common header]
```
In the next part we have some variables with boolian values. We need them to activate or deactivate 
some parts of the script. With this script we can have a look at the effect of qvalues, periodic boundaries
and the use of different STL files (3D model).  
```lua
--! [geometry definition]

qValues = true            -- activates qvalues
usePeriodic = true        -- activates periodic boundaries
periodicX = false         -- 
useObstacle = false       -- activates obstacle in the channel
case2d = false            -- creates a 2D channel
walls = true              -- walls are used
testIntersection = false  -- ?
verbose = true            -- ?
fixHeight = true          -- the channel's height is fixed and is used to compute dx

--! [geometry definition]
```
After that, we can specify some basic relations of the channel's geometry.
```lua
--! [geometry relations]

height = 1.               -- will be adapted to be exactly resolved by the discretization (physical)
ratio  = 8.               -- will be adapted to be exactly resolved by the discretization 
length = ratio*height     -- absolute length of the channel (physical)
radius = 0.1*height       -- radius of the channel if a tube is used (physical)

--! [geometry relations]
```
You might remember the `fixHeight` value from above. If it is set to `true` the `dx` is computed by
the specified height of the channel. Else `dx` is computed defining the length of the channel first.
```lua
--! [dx computed by height]

if fixHeight then

  nElemsHeight = 2^(level-4)                              -- number of elements in the height (lattice)   
  dx = height/nElemsHeight                                
  nElemsLength = nElemsHeight*ratio + 2                   -- need inlet and outlet element
  level = math.ceil( math.log(nElemsLength)/math.log(2))  -- compute the required level
  lengthSeeder = 2^level*dx                               -- length of a bounding cube element


--! [dx computed by height]
```
```lua
--! [dx computed by length]

else
  lengthSeeder = length/(1.-2./2.^level)                 -- length of a bounding cube element 
  dx  = lengthSeeder/(2^level)                           
  nElemsHeight = 2.*math.floor(length/(dx*2.*ratio))     -- number of elements in the height (lattice)
  height = nElemsHeight*dx                               -- absolute height of the channel (physical)
end 

--! [dx computed by length]
```
Here, we have some reference values that have effect on the geometry of the channel. 
You might change them and compare the effects later on in Paraview.
```lua
--! [reference values]

minlevel=level
maxLevel = level+math.max(refinementLevel, refinementLevel2) 
nElemsMax = 2^maxLevel
dxMin  = lengthSeeder/(2^maxLevel)
dxDash = 0.5*dxMin
lPhys = length  -- m
dxPhys = dx
heightLB = nElemsHeight
heightPhys = heightLB*dx

ebug = {debugMode = true, debugMesh = './prototree/'}
size_x = 0.1*length
start_x = -0.5*size_x
if useObstacle then
  size_x = 0.68*length
  start_x = -0.45*length
  size_y = 2.*height/5.-5.*dxDash
  if testIntersection then
    start_x = -0.40*length
  end
else
  size_y = height/5-5.*dxDash
end
start_y = -size_y/2
size_z = size_y
start_z = start_y

start2_x = -0.32*length
size2_x  = 0.37*size_x
size2_y = size_y*0.48
start2_y = -size2_y/2
size2_z = size2_y
start2_z = start2_y

--! [reference values]
```
Next, we define the bounding cube. 
The bounding cube, i.e. the box that includes the complete simulation domain.
It must be large enough, otherwise it just truncates all geometries that go beyond this cube.
```lua
--! [bounding cube]

bounding_cube = {origin = {-lengthSeeder*0.5, -lengthSeeder*.5, -lengthSeeder*0.5},
               length = lengthSeeder}

--! [bounding cube]
```
This part is followed by the `spatial_object` table. Here we define the seed, the channel itself
and obstacles. As you see, we have only the seed defined here. The missing parts depend on the use of 
periodic boundaries, obstacles and more that are activated and deactivated above.
```lua
--! [spatial objects]

spatial_object = {
  { attribute = { kind = 'seed', },
    geometry = { kind = 'canoND',
                 object = { origin = { 0.0, 0.0, dx*0.5 },
               }
    } -- geometry
  } -- seed
}

--! [spatial objects]
```
In this script we make use of an area with special refined mesh. Therefore we defined the levels above.
You can enter a lua table from any place after its definition. So we are able to deactivate
the following refinement boxes if we think that they are not necessary. To do this, you just make use of `--`
in front of each line.
```lua
--! [refinement box]

table.insert( spatial_object,
{
    attribute = {
      kind = 'refinement',
      level = refinementLevel+level,
      label='box1'
    },
    geometry = {
      kind = 'canoND',
      object = {
        origin = {start_x, start_y, start_z },
        vec = {{size_x, 0., 0.},
              {0.,size_y, 0.},
              {0., 0., size_z}}
      }
    }
  })

table.insert( spatial_object,
  {
    attribute = {
      kind = 'refinement',
      level = refinementLevel2+level,
      label='box2'
    },
    geometry = {
      kind = 'canoND',
      object = {
        origin = {start2_x, start2_y, start2_z },
        vec = {{size2_x, 0., 0.},
              {0.,size2_y, 0.},
              {0., 0., size2_z}}
      }
    }
  })

--! [refinement box]
```
Here we can define the inlet and the outlet of the channel. They are boundaries and belong to 
the `spatial_object` table. 
```lua
--! [inlet and outlet boundary conditions]
  
table.insert(spatial_object,  { 
    attribute = {
      kind = 'boundary',
      label='east' -- outlet
    },
    geometry = {
      kind = 'canoND',
      object = {
        origin = {length*0.5+dxDash, -length*0.5, -length*0.5},
        vec = {{0.0, length, 0.},
              {0.,0.0, length}}
      }
    }
    }) 

table.insert(spatial_object,  { 
    attribute = {
      kind = 'boundary',
      label='west' -- inlet
    },
    geometry = {
      kind = 'canoND',
      object = {
        origin = {-length*0.5-dxDash, -length*0.5, -length*0.5},
        vec = {{0.0, length, 0.},
              {0.,0.0, length}}
      }
    }
  })
--! [inlet and outlet boundary conditions]
```
It is possible for the rectangular channel to make use of periodic boundaries. Therefore you need to define 
two planes close to each other that are repeated along the channel length (`usePeriodic=true`) or the height
(`periodicX=true`). There exist these scripts:
```lua
--! [use periodicX boundaries]

if periodicX == true then
  table.insert(spatial_object, { 
    attribute = { 
      kind = 'periodic', 
      label = 'periodicX'
    },
    geometry = {
      kind = 'periodic',
      object = {
        plane1 = {
          origin = { -length/2+1.3*dx, -length/2, -length*0.5},
          vec = { 
            { 0.0, 0.0, length },
            { 0.0, length, 0.0 },
            }
        }, -- plane 1
        plane2 = {
          origin = { length/2-1.3*dx,  -length/2, -length*0.5},
          vec = { 
            { 0.0, length, 0.0 },
            { 0.0, 0.0, length },
            }
        }, -- plane 2
      } -- object
    } -- geometry

  }) 

end

--! [use periodicX boundaries]
```
```lua
--! [use periodic boundaries]

if usePeriodic == true then
  table.insert(spatial_object, { 
    attribute = { 
      kind = 'periodic', 
      label = 'periodic'
    },
    geometry = {
      kind = 'periodic',
      object = {
        plane1 = {
          origin = { -length/2, -length/2, dx+dxDash},
          vec = { { length, 0.0, 0.0},
                  { 0.0, length, 0.0},}
        }, -- plane 1
        plane2 = {
          origin = { -length/2,  -length/2, -dxDash},
          vec = { { 0., length, 0.0, 0.0},
                  { length, 0., 0.0},}
        }, -- plane 2
      } -- object
    } -- geometry

  }) 

--! [use periodic boundaries]
```
Of course you do not have to use periodic boundaries. Instead you could use a script that defines specific
boundaries for `top` and `bottom` of the channel. 
```lua
--! [non-periodic]

else
  if case2d then
    depth = 4*dx
  else
    depth=height
  end
  table.insert(spatial_object, { 
    attribute = { 
      kind = 'boundary', 
      label = 'top'
    },
    geometry = {
      kind = 'canoND',
      object = {
          origin = { -length/2, -length/2, depth/2+dxDash},
          vec = { { length, 0.0, 0.0},
                  { 0.0, length, 0.0},} 
        } -- object
      } -- geometry 
    } -- spatial object    
  )

  table.insert(spatial_object, {
    attribute = {
      kind = 'boundary',
      label = 'bottom'
    },
    geometry = {
      kind = 'canoND',
      object = {
          origin = { -length/2, -length/2, -depth*0.5-dxDash},
          vec = { { length, 0.0, 0.0},
                  { 0.0, length, 0.0},}
        },-- object
      } -- geometry
    } -- spatial object
  ) 
end

--! [non-periodic]
```
In order to simulate a flow around an obstacle you might activate `useObstacle` (`=true`). Here you can choose
between two stl files for a sphere (2D) or a cylinder (3D). You have to make sure that the files are located in the 
defined directory. Now it makes sense to activate qvalues. This allows to refine the mesh at i.e. circular
boundaries. This is done with `calc_dist = true` in the attribute table of a spatial object.
```lua
--! [use of obstacles]

if useObstacle ==true then

  if case2d == false then
    stlfile='stl/cylinder.stl'
    stlLabel='cylinder'
  else
    stlfile='stl/sphere.stl'
    stlLabel='sphere'
  end 

  if qValues == true then
    if stlLabel == 'cylinder' then
      table.insert(spatial_object,  { 
        attribute = { 
          kind = 'boundary', 
          level = maxLevel,
          label = stlLabel, 
          calc_dist = qValues,
        },
        geometry = {
          kind = 'stl', -- was: sphere
          object = { filename = stlfile
  --          { origin = {-length*0.3,-0.01*height,0.},
  --            radius = radius }
          }
        },
       transformation = {
          deformation =  2.0,
          translation =  {-2., 0., 0. }
        }
        }) 
    elseif stlLabel == 'sphere' then
      table.insert(spatial_object,  { 
        attribute = { 
          kind = 'boundary', 
          level = maxLevel,
          label = stlLabel, 
          calc_dist = qValues,
        },
        geometry = {
          kind = 'stl', -- was: sphere
          object = { filename = stlfile
  --          { origin = {-length*0.3,-0.01*height,0.},
  --            radius = radius }
          }
        },
       transformation = {
          deformation =  2.0,
          translation =  {-2., 0., 0. }
        }
        }) 
      -- elseif stlLabel == 'sphere' then
    end -- stlLabel
  else
--! [use of obstacles]
```
If you do not want to use qvalues you need this script:
```lua
--! [cylinderObj]
    table.insert(spatial_object,  { 
      attribute = { 
        kind = 'boundary', 
        level = maxLevel,
        label = stlLabel, 
        calc_dist = false,
      },
      geometry = {
        kind = stlLabel, -- was: sphere
        object = { filename = stlfile,
         -- { 
            --origin = {-length*0.40,-0.01*height,-lengthSeeder*0.5},
            --vec = {0.0,.0,lengthSeeder}, -- length and axis of the cylinder
            --radius = radius }
        }
      }
      }) 
--! [cylinderObj]
  end -- qValue == true
end -- useObstacle
```
Now you can start Seeder. If it successfully finishes, the simulation mesh was generated.

## Control the generated Mesh ##

Harvester (integrated in Seeder) can process the mesh and generate a VTK file, which can then be visualized by Paraview.
Let's prepare the harvester.lua for visualisation.
This is the simplest form of mesh visualisation. For more possibilities have a look at the
seeder documentation itself.

```lua
simulation_name = 'channel'

logging = {level=3}

mesh = 'mesh/'
-- define the input
restart = {
--          read = 'debug/final100_restart_ProtoData_header_100.lua'
--          read = 'debug/0_restart_ProtoData_header_0.lua'
        }

-- define the output_folder and output table if tracking table is not defined
output_folder = 'harvest/'
output = {  -- some general information
            format = 'vtk',   -- Output format 
         }
```
Now run `~/apes/seeder/build/sdr_harvesting harvester.lua`.

Visualize the mesh with paraview. 
If you set up `usePeriodic = false` in seeder.lua 
the VTU file should look like this:

![periodic_false](use_periodic_false.png "3D channel")

Otherwise you can build a two-dimensional channel like this:

![periodic_false](use_periodic_true.png "2D channel")

## Run Musubi ##

Let's have a look at how to configure Musubi. 

@note Many things were mentioned in the last tutorial

@note Describe here the most important parameters for running the simulation.

In order to repeat the most important things, you can make use of a checklist. 
For the simulation you need to define the 

  1. Basic information
  2. `mesh` The location of the directory is needed.
  3. `identify = { layout = "..", kind = "..", relaxation = ".."}`
     There exist default values, which are explained in the [last tutorial](tut_1_mus_config.html)
  4. Physical and lattice reference values
  5. time settings that describe the time interval, the time maximum and an abort criteria
  6. `fluid` table with information about physics and LBM units
  7. `initial_condition` table with the initial pressure and velocity
  8. `boundary_condition` This table aims to tell Musubi how to handle the boundaries that were set up
      in seeder.lua. It reads the boundaries' labels and converts them into a wall, an inlet or an outlet.
  9. `tracking` table for information about the target output. It exist a description module
      for tracking that you find [here](../treelm/page/features/tracking.html).
      
### 1. Basic information

Set up the simulation name that is used for the labels:
```lua
require "seeder" -- Tells musubi.lua to use defined variables from seeder.lua
simulation_name = simName
```

Specify the scaling (diffusive or acoustic) and interpolation method [[mus_interpolate_module]]:
```lua
scaling='acoustic'
interpolation_method = {
  method = 'linear', -- or 'quadratic'
}
```
We define local variables that we need to turn off or on a few features:
```lua
initChannel     = false
testIntp        = true
control_routine = true
track           = true
```

### 2. Location of the mesh

```lua
mesh = './mesh/'
```

### 3. Identify table

The corresponding elements are defined in the [[mus_scheme_header_module]] in [[mus_scheme_header_type]].
Here we define the label, layout [[tem_stencil_module]], 
the used model and the relaxation.

```lua
model       = 'lbm_incomp'
relaxation  = 'bgk'
scaling     = 'acoustic'

identify = {
  layout      = stencil, 
  kind        = model,
  relaxation  = relaxation
}
```

### 4. Physical and lattice reference values

```lua
--! [reference physical values]
u0Phys    = 1.
viscPhys  = u0Phys*heightPhys/Re
csPhys    = 300    -- m/s
rho0Phys  = 1.
--! [reference physical values]
--! [reference LB values]
csLB    = 1/math.sqrt(3)
cs2     = 1./3.
u0LB    = 0.05
rho0LB  = 1.
rho0LB0 = 0.
omega0  = 1.7 

-- Reynolds number
Re = heightPhys*u0Phys/viscPhys
ReLB = heightLB*uLB/viscLB
if verbose then
  print('           omega   = '..omega)
  print('          viscLB   = '..viscLB)
  print('        heightLB   = '..heightLB)
  print('             uLB   = '..uLB)
  print('Reynolds number Re = '..ReLB)
  print('          viscPh   = '..viscPhys)
  print('        heightPh   = '..heightPhys)
  print('        lengthPh   = '..length)
  print('             uPh   = '..u0Phys)
  print('Reynolds number Re = '..Re)
end  
```
Moreover we need these:
```lua
amplitude = u_in
--! [reference pressure]
if model == 'lbm_incomp' then
  p0 = 0.
else
  p0 = rho0LB*cs2*dx*dx/dt/dt
end 
--! [reference pressure]

originX = -length*0.25
originY = height*0.
originZ = 00.0
halfwidth = length*0.015
amplitude = 0.00
background = p0
function ic_2Dgauss_pulse(x, y, z)
  return background+amplitude*math.exp(-0.5/(halfwidth^2)*(( x - originX )^2+( y - originY )^2))
end
-- Reference values for the flow state in the 2d stationary channel at laminar
-- flow state
--! [velocity function]
function uX(x,y,z)
    uPhys = u0Phys*(1-(2.*y/heightPhys)^2)
    return uPhys
  end
--! [velocity function]
--! [pressure function]
function pressureRef(x,y,z)
  dp = u0Phys*8.*viscPhys*rho0Phys/heightPhys^2*length
  return p0 + dp*0.5 - dp/length*x
  end
--! [pressure function]
--! [boundary pressure]
pressureIn=pressureRef(-length/2,0.,0.)
pressureOut=pressureRef(length/2,0.,0.)
--! [boundary pressure]
function Sxx(x,y,z)
  return 0.
  end
function Syy(x,y,z)
  return 0.
  end
--! [shear stress function]
function Sxy(x,y,z)
  tauxy= -viscPhys*rho0Phys*8./heightPhys^2*u0Phys*y
  S_xy = tauxy/viscPhys/rho0Phys
  return S_xy
  end
--! [shear stress function]
function stressRef(x,y,z)
    return Sxy(x,y,z)*rho0Phys*viscPhys
  end
  
-- Consistent initial conditions for the channel
if initChannel then
  function ic_uX(x,y,z)
    return uX(x,y,z)
  end
  function ic_pressure(x,y,z)
    if( model=='lbm') then
      return pressureRef(x,y,z) + rho0Phys*cs2*dx/dt
    else
      return pressureRef(x,y,z)
    end 
    end
  function ic_Sxx(x,y,z)
    return Sxx(x,y,z)
    end
  function ic_Syy(x,y,z)
    return Syy(x,y,z)
    end
  function ic_Sxy(x,y,z)
    return Sxy(x,y,z)
    end
else
  function ic_uX(x,y,z)
    return 0.
  end
  function ic_pressure(x,y,z)
    return p0
  end
  function ic_Sxx(x,y,z)
    return 0.
  end
  function ic_Syy(x,y,z)
    return 0.
  end
  function ic_Sxy(x,y,z)
    return 0.
  end
end 
```

### 5. Simulation control

To define the time steps we need to disinguish between `diffusive` and `acoustic` scaling.
```lua
tEnd = 10.0
if scaling == 'acoustic' then
  uLB     = u0LB
  dt      = uLB/u0Phys*dx
  viscLB  = viscPhys*dt/dx/dx
  omega   = 1./(3.*viscLB + 0.5)
else
  omega   = omega0
  viscLB  = 1./6.*(2./omega - 1.)
  dt      = viscLB/viscPhys*dx*dx  
  uLB     = u0Phys*dt/dx
end 
```

We define the `sim_control` table [[tem_simControl_module]] containing information about `time_control` [[tem_timeControl_module]] and `abort_criteria` [[tem_abortCriteria_module]].

```lua
-- Time step settigs
interval = tEnd/10. 
tRamping = tEnd/10.
-- simulation will end when simulation wall clock time reaches 1 hr
-- or simulation time reaches tEnd
sim_control = {
  time_control = { 
    max = tEnd, 
    interval = interval, 
    clock = 3600 --s
  },
  abort_criteria = { 
    stop_file    = 'stop',
    steady_state = true,
    convergence = {
      variable = {'pressure_phy'}, 
      shape = {
        kind = 'canoND', 
        object = {{origin = {0.,0.,0.} }}
      },
      time_control = {
        interval = interval/20.,
        min = tEnd/2, 
        max = tEnd, 
      },
      format = 'convergence',
      reduction = {'average'},
      norm = 'average', 
      nvals = 100, 
      absolute = false,
      condition = { 
        threshold = 1.e-5, 
        operator = '<=' 
      }
    }
  }
}
```

### 6. Fluid table and physics 

We specify the [[mus_fluid_module]] information in this table.
At least `omega` or `kine_viscosity` is needed here.
```lua
fluid = { 
       omega = omega 
  }
```
The physics table is needed to convert between lattice and physical values.
See [[mus_load_physics]].
```lua
physics = { dt = dt, rho0 = rho0Phys }
```

### 7. Initial condition

Here we define the initial conditions for the simulation.
```lua
initial_condition = { pressure  = ic_pressure,
                      velocityX = ic_uX,
                      velocityY = 0.0,
                      velocityZ = 0.0,
                      Sxx = 0.,    
                      Syy = 0.,    
                      Sxy = ic_Sxy 
                      }
```

### 8. Boundary Conditions

Next, we define the boundary conditions to handle the boundaries we have set up in seeder.lua.
```lua
boundary_condition = {  
{  label = 'north', 
   kind = 'wall' },
{  label = 'south', 
   kind = 'wall' },
{  label = 'west', 
   kind = 'outlet_expol',
   pressure= 'pressureIn'},
{  label = 'east', 
   kind = 'outlet_expol',
   pressure= 'pressureOut'}}

if usePeriodic ==false then
  table.insert( boundary_condition, 
{ label = 'top', 
   kind = 'wall' } )
  table.insert( boundary_condition, 
{ label = 'bottom', 
   kind = 'wall' } )
end
if useObstacle ==true then
  table.insert( boundary_condition, 
{ label = 'sphere', 
   kind = 'wall_linearInterpolation' } )
 end
```

### 9. Tracking table

During the simulation we track variables to have a look at e.g. the pressure over time.
Therefore we use the [[tem_tracking_module]]. If we want to make use of different than the
default values we define the `variable` table that is explained [here](../treelm/page/features/variables/index.html).
```lua
variable = {
  { name='pressureRef', ncomponents=1, vartype = 'st_fun', st_fun = pressureRef },
}

if track then
-- Tracking
tracking = {
  {
    label = 'probePressure', 
    variable = {'pressure_phy','pressureRef'},
    shape = {kind = 'canoND', object = {origin ={0.0,0.,0.} } },
    time_control = {min = {iter=1}, interval = {iter=1}},
    output = {format = 'ascii'}, folder = './tracking/'                 
  },
 -- {
 --   ...
 -- },
  {
  label = 'hvs', 
  variable = { 
    'pressure_phy',
    'velocity_phy',
    'shear_stress_phy',
  },
  shape = 'global', 
  time_control = {min = 0., max = tEnd, interval = interval},
  output = {format = 'vtk'},
  folder = tracking_folder
 },
}
else
  tracking = { }
end
```

Now run the simulation with:
`~/apes/musubi/build/musubi musubi.lua`

## Post-process the Results ##

After the simulation run you can have a look at the vtk files with Paraview.
Here, you can see a possible output from Paraview:

![pressure_paraview2D](cylinder_pressure_screenshot_2D.png "2D channel with tracked pressure")

Now you are free to go on with the next chapter that deals with Tracking.

Next chapter: [Tracking quantities during the simulation run &rarr;](tut_3_tracking.html)
