lua_api.md

# GMSH DG solver Lua API

The solver employs the Lua programming language for its **configuration**. Lua was chosen first of all for simplicity: it is extremely lightweight and carries almost no dependencies, and this allow to keep the solver small and compact. Secondly, Lua was chosen to *deliberately* limit the possibilities of what the user can do in the configuration files. If configurations become full-fledged programs, it means that something is missing in the solver core or that the solver is being used in the wrong way.

This file documents the Lua API available in the solver. API has a *general* part and a *problem-specific* part.

The *general* part has to do with configuration not related to a specific problem (i.e. timestep, geometry file), whereas the *problem-specific* part configures all the parameter that make sense only on a given problem (i.e. materials and sources). This separation is reflected also in how the configuration is handled internally.

## API version
This document describes the API available on the version v0.1 of the solver.

## General interface

### Simulation variables

- `sim.name` (string): name of the simulation. Used also as the oname of the output directory.
- `sim.dt` (real): timestep duration.
- `sim.timesteps` (integer): number of timesteps to do.
- `sim.gmsh_model` (string): name of the file containing the GMSH model.
- `sim.use_gpu` (0/1): enable/disable GPU usage.
- `sim.approx_order` (integer): method approximation order.
- `sim.geom_order` (integer): geometry order. Only order 1 supported for now.
- `sim.time_integrator` (string): `"rk4"`, `"leapfrog"` or `"euler"`. Euler is only for test purporses, as it is numerically instable. Leapfrog works on CPU for now, use it if you have dissipative materials.

### Postprocessing general variables

- `postpro.silo_output_rate` (integer):
- `postpro.cycle_print_rate` (integer):

### Mesh variables
- `mesh.scalefactor` (real): apply a scaling factor to the mesh.

### Callable functions
- `enable_boundary_sources()`: takes a boolean parameter that specifies if boundary sources should be enabled or not
- `enable_interface_sources()`: takes a boolean parameter that specifies if the sources applied on internal interfaces should be enabled or not
- `enable_volume_sources()`: takes a boolean parameter that specifies if volumetric sources should be enabled or not

### Callbacks
None.

## Maxwell solver interface

### Materials
Materials are specified populating the `materials` table exposed by the solver. The material parameters are the relative electric permittivity `epsilon`, the relative magnetic permeability `mu` and the conducibility `sigma`. The empty space parameters are stored in `const.eps0` and `const.mu0`.

It is possible either to provide a function describing the material properties in the whole domain or specify the materials subdomain by subdomain. In the first case, appropriate callbacks have to be installed in the material table, whereas in the second case the appropriate variables must be set.

If, for example, the model has two subdomains with tags `1` and `2`, the electric permeability can be specified as

```
function epsilon_callback(tag, x, y, z)
    if tag == 1 then
        return 1
    return 2
end

materials.epsilon = epsilon_callback
```
or by setting the subdomain-specific variables as
```
materials[1] = {}
materials[1].epsilon = 1
materials[2] = {}
materials[2].epsilon = 2
```
Note that in this second case the subdomain-specific tables must be initialized with `{}` before setting any variable.

The two approaches can be mixed and the subdomain-specific variables take priority on the whole-domain function.

The solver expects that the callbacks for materials have 4 parameters:
- `tag`: the current subdomain tag
- `x`, `y` and `z`: the current coordinates at which the material needs to be evaluated. Note that currently the materials are assumed to be piecewise constant and that the material is evaluated in the barycenter of the element.

### Volumetric sources
Current sources in volumes can be specified by providing a function that evaluates the current density field in a given subdomain. For example the following code defines a volumetric source in the subdomain with tag `4`:
```
function J_source(tag, x, y, z, t)
    local Ex = 0
    local Ey = 0
    local Ez = math.sin(2*math.pi*freq*t)
    return Ex, Ey, Ez
end

sources[4] = J_source
```

### Boundary conditions
Boundary conditions are specified for each surface with tag `tag` via the table `bndconds`. If a surface does not have an entry in `bndconds` it is assumed to be PEC.

For example, to specify an impedance boundary condition on the surface with tag `4`, one proceeds as following:
```
bndconds[4] = {}
bndconds[4].kind = "impedance"
```

There are also nonhomogeneous boundary conditions, as the plane wave condition. A plane wave described in terms of electric field is configured as
```
function pw_callback(tag, x, y, z, t)
    --- Computation of the field components
    return Ex, Ey, Ez
end

bndconds[4] = {}
bndconds[4].kind = "plane_wave_E"
bndconds[4].source = pw_callback
```
The available homogeneous conditions are:
- `"pec"`: Perfect Electric Conductor.
- `"pmc"`: Perfect Magnetic Conductor.
- `"impedance"`: Impedance condition, actual impedance value is computed using the material parameters of the adjacent element.

The available nonhomogeneous conditions are:
- `"E_field"`: electric field. Not implemented yet.
- `"H_field"`: magnetic field. Not implemented yet.
- `"surface_current"`: surface current.
- `"plane_wave_E"`: plane wave specified by the electric field. Normal incidence is assumed.
- `"plane_wave_H"`: plane wave specified by the magnetic field. Not implemented yet.

### Interface conditions
The interface conditions follow the exact same logic as boundary conditions, except that they are applied on internal interfaces. The table in this case is `ifaceconds`.

Currently, only a surface current condition is available:
- `"surface_current"`: surface current.
### Callable functions
None.

### Callbacks
- `on_timestep(ts)`: if defined, the solver calls this function at every timestep. The solver passes the current timestep number in the parameter `ts`. The function is called **before** the solver state is advanced.
- `electric_initial_condition(x, y, z)`: if defined, the solver calls this function at the beginning of the simulation in order to initialize the electric field. The parameters `x`, `y` and `z` represent the point at which the field needs to be evaluated. The function has to return a triple specifying the three field components, i.e. `return Ex, Ey, Ez`.
- `magnetic_initial_condition(x, y, z)`: if defined, the solver calls this function at the beginning of the simulation in order to initialize the magnetic field. The parameters `x`, `y` and `z` represent the point at which the field needs to be evaluated. The function has to return a triple specifying the three field components, i.e. `return Hx, Hy, Hz`.

### Postprocessing
The export of the fields E, H and J to Silo is controlled using the appropriate sections in the `postpro` table.
- `postpro["E"].silo_mode` controls the export of the electric field E. Possible values are `"none"` (no export), `"nodal"` (export as a Silo nodal variable) and `"zonal"` (export as a Silo zonal variable) 
- `postpro["H"].silo_mode` controls the export of the magnetic field H. Possible values are `"none"` (no export), `"nodal"` (export as a Silo nodal variable) and `"zonal"` (export as a Silo zonal variable) 
- `postpro["J"].silo_mode` controls the export of the current density field J. Possible values are `"none"` (no export), `"nodal"` (export as a Silo nodal variable) and `"zonal"` (export as a Silo zonal variable)

Nodal variables represent the field at the mesh nodes. The nodal value is reconstructed by taking the average of the solution evaluated at the node in all the elements adjacent to that node.

Zonal variables represent the field as a piecewise-constant field evaluated in the barycenter of the element.

Default mode is `"nodal"`.

### Considerations on the evaluation of sources and the computation on GPU
Currently sources for the timestep `t+1` are evaluated asynchronously on the CPU while the GPU is computing the timestep `t`. The sources at `t+1` are then uploaded asynchronously from the CPU to the GPU. This usually works relatively well, especially on big domains with high approximation order, however on small domains at low approximation order you can notice slowdowns. For this reason, if the sources do not need to be evaluated along the whole simulation, it is suggested to turn them off at the appropriate timestep by using the `on_timestep()` callback and the appropriate functions described below. Despite of that, I am not sure that implementing repetitive boundary sources will be necessary.