An exploration to visualize finite element data with a DSL

TL;DR

This paper explores extending the Diderot language to visualize finite element data, demonstrating simple sampling and volume rendering techniques to enhance scientific visualization capabilities.

Contribution

It introduces an extension of the Diderot language specifically for finite element data visualization, enabling new visualization methods for scientific computations.

Findings

Successful visualization of FE data using Diderot

Demonstrated simple sampling and volume rendering techniques

Extended Diderot to support FE data visualization

Abstract

The scientific community use PDEs to model a range of problems. The people in this domain are interested in visualizing their results, but existing mechanisms for visualization can not handle the full richness of computations in the domain. We did an exploration to see how Diderot, a domain specific language for scientific visualization and image analysis, could be used to solve this problem. We demonstrate our first and modest approach of visualizing FE data with Diderot and provide examples. Using Diderot, we do a simple sampling and a volume rendering of a FE field. These examples showcase Diderot's ability to provide a visualization result for Firedrake. This paper describes the extension of the Diderot language to include FE data.

Figures

Equations

$$-\nabla {}^{2}u+u=f$$

$$\nabla u\cdot \vec{n}=0\text{ on }\Gamma .$$

$${\int }_{\Omega }\nabla u\cdot \nabla v+uvdx={\int }_{\Omega }vfdx$$

$$f=(1.0+8.0{\pi }^2)cos(2\pi x)cos(2\pi y)$$

Full text

We demonstrate our first and modest approach of visualizing FE data with Diderot and provide examples. Using Diderot, we do a simple sampling and a volume rendering of a FE field. These examples showcase Diderot’s ability to provide a visualization result for Firedrake. This paper describes the extension of the Diderot language to include FE data.

1 Introduction

Computational scientists compute solutions to systems of partial differential equations (PDEs) on large finite meshes using numerical techniques, such as the finite element method (FEM). These PDEs can be used to describe complex phenomena like turbulent fluid flow. The solution to PDEs or the output to software that solves PDEs are sometimes referred to as “FE fields”. The Diderot language does not know how to solve PDEs, or represent FE fields, but it does have the computational model to visualize fields. The work in this paper takes a step towards using Diderot to visualize FE fields. With this work we hope to use Diderot to help debug visualizations of FE fields and enable more interesting visualizations.

Solving PDEs and visualizing PDEs require two different techniques and entirely different code. To simplify the transition from solving a PDE to visualizing its solution, scientists may turn to standard visualization tools to analyze their data. The problem with this approach is that there is not a universal solution to accurately visualize every PDE. For example, visualizing finite-element data created with higher order elements and a small number of cells can lead to images that do not accurately represent the original solution. We believe that Diderot can be useful.

We want the user to be able to use visualization programs enabled by Diderot on fields created by FEM. Expecting users to transition from visualization toolkits to writing in a new programming languages and developing an expertise in scientific visualization is a big ask. Our goal is to be able to augment any existing Diderot program (written for discrete data) and apply it to FE data with minimal changes to the program. That way Diderot could compile programs to extract interesting visualization features from FE data [13].

Our work demonstrates a modest step towards visualizing FE data with Diderot. This paper is organized as follows. Section 2 offers some background about this problem area and provides a motivating example. Section 3 describes the implementation details to our approach. Section 4 demonstrates an application of our approach by providing one example of the Helmholtz equation and interpolating a function. We end with a discussion in Section 5.

2 Motivation

The Finite Element Method [4] gives a general framework for computing solutions to differential equations. In Section 2.1, we more closely describe FE fields and the software used to create them. In Section 2.2 we describe the state of using visualization toolkits to visualize a certain domain of FE fields. We provide an example of a problem that can not be visualized with the current state of visualization.

2.1 Background

Diderot

Diderot is a language designed for scientific visualization and image analysis. Algorithms in this domain depend on a high level of math. Diderot eases the transition from visualization algorithms to code by providing the user with a rich mathematical notation. We recently provided examples of more mature visualizations [13] that can created with Diderot due to developments with our intermediate representation EIN [5,6].The computational core to visualization programs involve tensor and field operators applied to fields. In Diderot, we have previously defined fields as the convolution between discrete image data and a kernel. Our goal is to extend our definition of a field to include other types of data.

Diderot has approached evaluating the correctness in the Diderot compiler in two ways; an automated testing infrastructure [7] and presenting properties of the rewriting system [9].There has yet to be an evaluation of an image, created by an outside source. In this paper we evaluate the results based on visualization principles [14].

FE fields

FE fields are created from a solving a PDE on a finite element mesh, which involves discretizing the domain into small finite mesh elements and using a set of basis functions (derived from the mesh) to span the domain space. These fields approximate numeric solutions to PDEs.

The Python code builds on the description of the problem. There are a variety of different meshes that are built-in or could be created by outside tools. Often, we use a unit square mesh.

⬇

m = UnitSquareMesh(2, 2)

The code creates a 2x2 mesh of a square. Each smaller square is divided into two triangles for a total of eight elements. The mesh is one of the arguments when defining a function space.

⬇

V = FunctionSpace(m, $``P"$, K)

$``$P” refers to a family of finite element spaces (other families include $``$DP”, $``$RT”, and $``$BDM”). The basis functions are linear when K=1 and cubic when K=3.

A FE field can be created from solving a PDE, but it can also be generated by interpolating an analytically defined expression.

⬇

f = Function(V).interpolate(Expression(exp))

In this chapter, when we wish to create simple examples we choose to generate fields from expressions.

Software

Computer scientists build software to solve PDEs that represent a wide range of problems. On the surface the software (or programming language) ideally represents a high-level math notation that is easy to understand but under the hood it is more complicated. The translation between the notation and computer code involves several steps and pieces of software. Solving a PDE involves the discretization of differential equations and uses the finite element method to provide an approximate solution. Optimizing the translation from PDE equation to approximate solution is pursed by many groups including the FEniCS Project [10] and Firedrake [19].

The FEniCS [10,17] project is an automated system to find solutions for partial differential equations using the finite element method. It enables users to employ a wide range of discretization to a variety of PDEs. On the surface it uses the Unified Form Language (UFL) [2],a domain-specific language to represent weak formulations of partial differential equations. UFL supports tensor algebra, high-level expressions with domain-driven abstraction. UFL does not provide the problem solving environment, instead it creates an abstract representation that is used by form compilers, such as the FEniCS Form compiler (FFC) [3],to generate low-level code. FFC can implement tensor reduction for finite element assembly [16] and aims to accept input from any multilinear variational form and any finite element to generate efficient code.

Firedrake [19] is a similar system that is also used to solve PDEs. In addition to UFL, it uses a modified version of FFC [3], FIAT [15], a PyOP2 interface [20], and COFFEE [18] 111We describe the Firedrake project at the time of our work but Firedrake has developed new parts of their infrastructure: FINAT augments FIAT and their new form compiler TSFC replaces FFC.. FFC is the FEniCS form compiler for generation of low-level C kernels from UFL forms. FIAT is the finite element automatic tabulator. It presents an abstract description of elements and has a wide range of finite element families. PyOP2 provides a framework for carrying out parallel computations on unstructured meshes. The COFFEE compiler optimizes the abstract syntax trees generated by FFC.

2.2 Creating and Visualizing FE data

There are a number of approaches to supporting scientific visualization. A common way is to use a toolkit such as the Visualization Toolkit (VTK) [21],ParaView [1], and the Insight Toolkit (ITK) [12] In this domain it is important for the visualization tool to understand how the data is represented and that restraint limits the options available to the FEM user. Toolkits are commonly used to visualize the solutions created by FEM software.

Existing practices to visualize FEM are insufficient. The strategy to solve PDEs is very different from the algorithms used to visualize the result. Firedrake uses a VTK file format for its visualization output. The format only supports linear and quadratic data. Firedrake takes the output and writes the output to a linear file format222 The images in this paper are created when Firedrake did an L2 projection, but now Firedrake uses interpolation to generate linear output.. Paraview might then accurately assume linear basis functions to represent the Firedrake output even though the original solution was created with higher-order elements. As a result the image may not accurately represent the PDE solution.

In this section, we present the problem with the current state of visualizing FE data by providing a simple example of a how a bug can occur. We create a field using higher-order data with a small number of cells. We present the results of using Paraview and Diderot side by side. The implementation process that enables the communication between Diderot and Firedrake is introduced in more detail in Section 3.

Creating FE data

In the following example, we use Firedrake to create FE data. Instead of solving a PDE, we define a field by interpolating the expression given the function space. We show the Python code in Figure 1.

We define a polynomial expression $x^{2}(1-x)$. Given a higher-order polynomial expression we can also assume that linear functions will not correctly be able to represent it. On the other hand, cubic functions should be able to offer a reasonable approximation for this problem.

Visualizing FE data

We chose to visualize the expression several different ways in order to provide a means of comparison. Figure 2 and Figure 3 present these visualizations. The images in Figure 2 offer the ground truth for this example. Figure 3 uses Firedrake data created with linear and cubic elements333As a note subdividing the field would create a more accurate solution but the context of the problem required creating a field with higher-order elements and a small number of cells.. Regardless of the visualization strategy we expect the images (left) created with linear elements to be an inaccurate representation of the solution. On the other hand, we expect the images created with cubic elements to be correct. The image (bottom right) created with Diderot is correct while the image (top right) created with Paraview is not.

When the results do not represent the field it can be difficult to understand and use visualizations to debug. From the user’s perspective the issue could be with the user’s UFL code, Firedrake’s evaluation, or the visualization program. Our goal is to provide Diderot as an alternative tool that can be used in these instances.

3 Our Approach

It is not a long term solution but we created our current prototype to establish communication between Firedrake and Diderot. Firedrake creates FE data by solving a PDE or interpolating an expression over a function space. Diderot visualizes the results and provides a few programs that can be initialized. Firedrake provides a way to evaluate a field at a position and provides Diderot with a call back. In the Diderot program a field is defined for FE data instead of using field reconstruction and convolution. Everything else in the Diderot program is the same. Lastly, when possible we try to evaluate the visualization based on established visualization principles. We worked in collaboration with the Firedrake team [19] at the Imperial College of London.

FE data

A field is created using FE data with a Firedrake program. The field can be the result of solving a PDE or interpolating an expression over a function space. The Diderot syntax offers a limited way to represent tensor fields created by FE data.

Diderot program

The core of visualization programs is independent of the source of data, but Diderot could only represent regular imaging grids. We addressed this limitation by allowing the Diderot compiler to declare a field that is defined by FEM. Instead of creating a field by using field reconstruction and convolution the user can use the following Diderot line to indicate the new kind of field data.

⬇

**input** **fem**#0(2)[] g;

**field**#0(2)[] F = toField(g)

As a result, a visualization program written for regular imaging grids could be easily changed so it is applied to the result of a PDE solution. We have created a few Diderot programs that have been compiled to a library.

Firedrake program

The Diderot program provides the framework to sample fields and do volume rendering of 3-D fields. A Firedrake program makes a call to the Diderot library. A field reference is used to initialize a call to the Diderot library. The following is a sample of Python code that would be written in a Firedrake program.

⬇

res = 200 # resolution to MIP program

step = 0.01 # step size to ray tracer

diderot.mip(file_name, f, res, res, step)

It initializes visualization parameters (such as resolution and step size) and provides a pointer to the field.

Point Evaluation

Firedrake can evaluate fields (represented as a FE field) at arbitrary physical points [11] It determines which cell to look at. Diderot does not know how to evaluate a FE field at a given position. For each inside test, probe, and gradient operation the field has to be evaluated at a position. The generated code will use Firedrake’s point evaluation functions for a given field and position.

4 Demonstration

We provide two examples to demonstrate results of this work. In the first example, a field is created by interpolating an expression over an indicated function space and visualized with a volume rendering program. In the second example, we provide a classic example of a PDE and simply sample the result.

4.1 Communication between Diderot and Firedrake

Python code

We create the field defined by interpolating the expression given the function space using Firedrake. We show the Python code under the “Firedrake Field” box in Figure 6. We make a call to the Diderot library by passing the field and initializing the resolution and step size (for the Diderot program).

Diderot code

Diderot visualizes the result by using an augmented MIP program. MIP or maximum intensity projection is a minimal volume visualization tool for 3-d scalar images. The Diderot code for the Diderot program is presented in Figure 4.

Results

We visualize the FE field using a Diderot implementation of a MIP program. We vary the expression created by Firedrake, where the camera in the Diderot program is pointed, and how the data is included in the visualization program. We try these different variations to mimic how a user might use an existing Diderot program (or template) to visualize their data.

We provide the results in Figure 5. We set the camera to look at the center of the data (center and right column) or off-center (left column). We also created fields with the Firedrake expression centered (two top rows) and off-centered (third row).

The red box indicates the setting where the Diderot camera and firedrake expression are both centered and result is as expected (symmetric sphere). The other images in the figure are subject to bugs, because the function space defined by Firedrake does not match the image space probed by Diderot (or where the camera is pointed). This example illustrates the care that is needed when setting up the Firedrake Diderot pipeline.

Evaluate Visualization

The data representation, either by FE field, or discrete field, should not change the visual representation (according to the Principle of Representation Invariance). We compare the output for the field created by Firedrake and the field created by Diderot in Figure 6. As can be expected, the image created by Firedrake is essentially the same as that created by Diderot.

4.2 PDE Example

The Helmholtz problem is a symmetric problem and is a classic example of a PDE. Consider the Helmholtz equation on a unit square $\Omega$ with boundary $\Gamma$.

[TABLE]

[TABLE]

The solution to the equation is some function u $\in V$ for some suitable function space V that satisfies both equations. After transforming the equation into weak form, applying a test function $V$, and integrating over the domain we get the following variational problem444The approach to solve PDEs with FEM is described in more detail in by Brenner and Scott [4]..

[TABLE]

We choose function f555The example and code are provided by Firedrake

[TABLE]

The Python code (shown in Figure 7) solves the PDE and connects the solution to Diderot. The original code uses linear elements (k=1), but we choose to use linear and cubic elements.

We illustrate the results using Paraview and Diderot. The results are shown in Figure 8. The image on the bottom right is a difference image and illustrates the difference between the Diderot results. The data with higher-order elements and visualized with Diderot is the more clear and refined. There is a smoothness captured with the higher-order data that is not in the linear data.

5 Conclusion

The Finite-element community has difficulty visualizing fields that use higher order elements. The communication pipeline allowed us to demonstrate Diderot’s ability to visualize FEM data by sampling the data and doing volume renderings. We have shown that Diderot can be used to correctly visualize fields created with higher order elements and a small number of cells. We have also shown a way to evaluate the results for simple test cases by comparing the result with an external data representation.

It would be beneficial if a Firedrake user could use Diderot with little hassle. We proposed that we could ease this transition by providing a Diderot template (or existing Diderot programs) that can be initialized by Firedrake code, however, as we have shown in the example in Section 4.1 this approach can still lead to errors. It takes care, by the user to ensure that a Diderot program is set-up to correctly visualize the Firedrake data.

We have successfully established communication between Firedrake and Diderot but it is limited. We suggest further ideas for future work that help develop Diderot and benefit the community.

Debugging FE data with Diderot

The work on visualizing FE data opens the door to more useful applications. Diderot could potentially be used to visually debug and validate fields created by FEM and possibly reveal hidden details in data created with higher-order elements. We have shown an example of using Diderot to correctly visualize a field created by Firedrake, but have not explored its full potential.

Mature visualization for FE data

We have not visualized FE data with Diderot programs that use higher-order code. Partly, because the Diderot compiler cannot yet represent FE data in more intimate level. In the future, we look forward to being able to apply more complicated visualization programs to FE fields.

The callback that Firedrake offers, to evaluate a field at a point, is also limited. As far as we know, Firedrake does not offer a call back to take third or fourth derivatives of fields (which is necessary in some visualization programs as shown in [13].Additionally, the call back can only evaluate some fields. There is no callback supported for some niche problems such as extruded meshes, manifolds, and mixed element meshes. These type of structures could possibly require more syntax consideration on the Diderot side as well.

Better communication between Diderot and FEM

A call is made to Firedrake every time a field needs to be evaluated, which makes the communication process expensive. Each function call creates multiple tensor operations in order to do the right transformations and find the right cell. The change in coordinates from a reference element to one being computed involves the calculation of the Jacobian matrix, its determinant, and its inverse.

These operations can be similar to previous calls, leading to redundant and expensive computations. Diderot does not know how to evaluate a FE field at a point or have the syntax to find the right cell, but it can represent tensor computations. If we are able to move some of the computations into Diderot then Diderot can catch these redundant computations and the entire process is less expensive.