Domain Specific Languages

In order to overcome the Programming and Portability walls, I studied and developed different tools to assist the user when writing his code, in a high level language such as Python, without sacrificing the performances. Moreover, the user can get a low level code that he/she may integrate to an existing code without having to go through the trouble of developing everything whenever they are new changes in the studied model.

In order to make this possible, we created two important libraries:

  • First, we developed a Python-Fortran converter, which is in fact, a static compiler based on Fortran.
  • A mathematical language that capture the semantic of mathematical expressions. Its design was inspired by notions from Category Theory and Homotopy Type Theory.

As an example, we consider the Poisson equation

(1)   \begin{align*}-\nabla^2 u = f \quad \text{in~$\Omega$}, \quad \quad  u = 0            \quad \text{on~$\Gamma_D$}, \quad \quad  \partial_n u = g \quad \text{on~$\Gamma_N$}.  \end{align*}

An H^1-conforming variational formulation of (1) reads

(2)   \begin{align*} \text{find $u \in V$ such that} \quad a(u,v) = l(v) \quad \forall v \in V,  \end{align*}

where V \subset H^1(\Omega),
a(u,v) := \int_{\Omega} \nabla u \cdot \nabla v ~ d\Omega, and
l(v) := \int_{\Omega} f v ~ d\Omega + \int_{\Gamma_N} g v ~ d\Gamma.

The associated formal model is given in the following Python code, using the Python library Sympde

# define a topological domain as a square
domain = Square()
# Neumann boundary 
B_neumann = domain.get_boundary('Gamma_1')

# Dirichlet boundary 
B_dirichlet = domain.boundary.complement(B_neumann)

# define a scalar function space
V = ScalarFunctionSpace('V', domain)

u,v = element_of_space(V*V, ['u', 'v'])

a  = BilinearForm((u,v), dot(grad(v), grad(u)))
l0 = LinearForm  ( v,    f*v)

# boundary term
l_B_neumann = LinearForm(v, v*trace_1(grad(solution), B_neumann))

l = LinearForm(v, l0(v) + l_B_neumann(v))

bc = EssentialBC(u, 0, B_dirichlet)

# declare the variational problem as an equation
equation = find(u, forall=v, lhs=a(u,v), rhs=l(v), bc=bc)

The next step now, is to have a discrete version of our formal model. This is provided by the Python library Psydac. The user can call the magic function discretize to convert abstract notions to discrete concepts, by adding some parameters as in the following Python code:

# create the computational domain from a topological domain
domain_h = discretize(domain, ncells=ncells, comm=comm)

# discrete spaces
Vh = discretize(V, domain_h, degree=degree)

# dsicretize the equation using Dirichlet bc
equation_h = discretize(equation, domain_h, [Vh, Vh])

# discretize norms
l2norm_h = discretize(l2norm, domain_h, Vh)
h1norm_h = discretize(h1norm, domain_h, Vh)

# solve the discrete equation
x = equation_h.solve()

# store the solution in a Field
phi = FemField( Vh, x )

# compute norms
l2_error = l2norm_h.assemble(F=phi)
h1_error = h1norm_h.assemble(F=phi)

The previous Python code will generate automatically a Fortran code, and can also run in parallel using MPI. By doing so, we can easily perform numerical simulations and study the behavior of the linear system.

The following diagram shows how one can go from the formal language to a discrete representation

The spectral analysis is performed the same way