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
An -conforming variational formulation of (1) reads
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