Google Summer of Code 2021 - Community Bonding

It’s time for a GSoC project update! May 18 to June 8 in the Google Summer of Code timeline was the community bonding period, where the student spends about a month knowing more about the organization and the community. In the previous post, I spoke about the general application process, converging on a topic for the summer of code. In this short post, I will be talking about what goes on during the bonding period. Note that this can vary between the various organizations in GSoC, for example, differing in the tech that they use for communication etc.

I will also try to briefly discuss the technical details of the work that I’ve been working on so far. I don’t intend to explain the entire thing in this single blog post, as it would make it long and hard to follow. I will however try to make a separate page documenting the routines, if possible.

Overview

To get started, I needed to know how to contribute to a project. This step is almost uniform throughout the open-source community. There are several articles online that describe the various terminologies I use throughout this post. An excellent place to start with is the official Github’s documentation page

• Fork a repo.

Gridap uses Github for hosting their projects. The first step was to create a fork of the official repository. A fork is a copy of the remote repository in your account, which can then be cloned to make a local copy on your computer. Any contributions/changes that you want to make to the code will be done locally and pushed to the fork. The next step is to submit a pull request to the original repository from your fork once you commit the local changes. The admins then review the pull request and merge the changes if everything looks good.

Task Lists

Aside from communicating through Slack, my mentors advised me to create a task list in the issue tracker of the fork, similar to the one here. This list is to keep track of where I am with the implementation. Github flavored markdown offers support to create a task list.

• Tasklist.

Evaluation of functions at arbitrary points

The first ingredient towards building the fast finite element interpolator in Gridap is to evaluate the function at arbitrary points in the domain. The algorithm works by finding the nearest neighbour of the point in a KD-Tree, built using the mesh vertices. Subsequently, we construct an inverse map to get the corresponding cell/element location and then evaluate the function using Gridap’s evaluate method. To make things simple, it is easy to think of points as vertices, and the vertices make up the cell/element. The goal is to make evaluate work on any point inside the cell.

This algorithm was written by Erik Schnetter for Gridap a couple of months ago. I will go through the algorithm briefly here. Feel free to skip to the conclusions if you are not familiar with FEM. The first step is to build a mesh of an arbitrary dimension D. The code snippet below works for any value of D, subject to computational power of your computer of course.

D = 3
domain = repeat([0, 1], D) # domain = [0,1]^D
partition = repeat([3], D) # 3 cells
model = CartesianDiscreteModel(domain, partition)

The next step is to create a finite element space on the mesh.

reffe = ReferenceFE(lagrangian, Float64, 2)
V = FESpace(model, reffe)

This creates a Lagrangian Finite Element Space of order 2. The reference finite element space is where most of the computation occurs and this makes the implementation efficient. Now let us define an affine function on the finite element space. We chose this for the test, since affine functions are represented exactly on a second order FESpace.

coeff0 = rand(Float64)
coeffs = rand(SVector{D,Float64})
f(x) = coeffs ⋅ SVector(Tuple(x)) + coeff0
fh = interpolate_everywhere(f, V)

The interpolate_everywhere method in Gridap interpolates the function on the degrees of freedom of the finite element space. The new evaluate method then works if we specify an arbitrary point/array of points in the domain, say something like

evaluate(fh, Point(√2-1, √3-1))
evaluate(fh, [Point(√2-1, √3-1), Point(0.121,0.123)])

The signature of the new method would be

evaluate!(f::CellField, xs::AbstractVector{<:Point}) ## xs can be a vector of points

By default, the evaluate method works for inputs of the form

evaluate(get_data(fh), Point(√2-1, √3-1))

but returns the piecewise finite element function evaluated on all the cells/elements. A k-dimensional tree is then used to perform an efficient search of the cell where the point is located in. This is done mainly using three new methods present in the Gridap.CellData module

_point_to_cell_cache(trian::Triangulation)
_point_to_cell!(cache, x::Point)
make_inverse_table(i2j::AbstractVector{<:Integer},nj::Int)

The _point_to_cell_cache method takes in the triangulation (of Gridap type Triangulation) and returns the KDTree of the mesh vertices along with other caches required to compute the location of the point in the triangulation. The method _point_to_cell! takes in the cache returned by _point_to_cell_cache along with the point itself, and returns the cell closest to the point. The method make_inverse_table then computes the global and local positions of the point in the cell, with local being the position of the point in the input array (remember evaluate now works for an array of points). This is useful if there are multiple points inside the same cell. Once the appropriate cell-positions cell_to_points are found, the old evaluate function can then be used,

cell_to_xs = lazy_map(Broadcasting(Reindex(xs)),cell_to_points)
evaluate(get_data(fh), cell_to_xs)

where cell_to_xs is the table containing the final positions of the input points xs in the cell array. This returns the function value on the queried points xs.

One of the key functions in Gridap is the lazy_map, which is described in detail in the developer tutorials

• Developer Tutorials

In short, lazy_map is one of, if not the most important function in Gridap. It returns a data structure known as LazyArray which are arrays that are conceptually defined at the indices. This means that they are not actually computed as a whole, but rather contain the information on how it is built. The actual value of the array is computed only when the indices are called, and hence the term Lazy. This results in more performant code as opposed to greedy evaluations which result in memory intensive and slower code.

Conclusion

This work was started a couple of months ago in a pull request by Erik Schnetter, who wrote most of the algorithm. I then volunteered to finish up Erik’s implementation, mostly cleaning up the code, extending the method to CellPoint structures, and completing the unit tests. The pull request details are below if you would like to have a look at them.

• Pull Request by Erik Schnetter
• My Pull Requests to Erik’s branch, #1 and #2

With this implementation, the geometrical part is complete, and my new task list is here. The objective of the bonding period was to take time and get introduced to the codebase, read the contributing guidelines, submit and review pull requests etc. Once all the tests have passed, the last step is to add the list of new features in the changelog.

• See this: keep a changelog.

Over the past month, I learned a lot about how contributing to open-source organization works, and also had a chance to interact with the developers of Gridap. I also thank my mentors for guiding me towards submitting my first ever successful pull request! The official coding period starts on June 8. Super stoked for GSoC this year!