In the last three episodes about FEniCS we introduced the software and installed it. We also derived the weak form for the Helmholtz equation and developed some code to solve it. This is all nice and good, but writing scripts from scratch every time can be tedious and error prone. We can help ourselves by writing reusable code packages to serve as tools. For example, we could write a software package to solve the Helmholtz equation with FEniCS. Or any other package we might find useful. The idea is that we can create an entire toolbox to deal with common tasks when creating simulations. We started the development of some tools under the acoupy project. This series of articles will serve as a sort of “development” diary, illustrating how the packages are being developed and why. This will allow us to illustrate how to develop scientific software.

Why?

Computational science and engineering were born well before the era of electronic computers. Originally, computers were humans that would manually solve numerical problems. For example, the main idea behind FEM originates in the early 1900s. In this day and age instead we use electronic computers to carry out the computations. Hence, making numerical simulations involves writing code.

There exists a lot of ready to use code that one can use for simulations. In the previous articles we focused on Elmer, which is a great example. However, it is important to be able to develop, own and maintain code as well. This because we might be interested in pre-processing data to feed to our simulators. Or post-process the output of our simulations. Or create simulation code ourselves when the existing packages do not quite do what we need. An example would be Elmer acoustic solvers lacking Perfectly Matched Layers (PML) or Adiabatic Absorber boundary conditions for reflection-less boundaries. If our simulation needs such boundary condition then we need to code it ourselves. With good coding skills we can develop entire ecosystems of packages that can serve our purposes and ease the development of our simulations.

Unfortunately this skill is rarely taught to people that do not come from a software development background. As a result, code developed by engineers and scientists often results in being poor under one or more aspects. For example, the code might work but not being very readable or maintainable. Developing good software is very hard. But the basics are simple to pick up. In this series we will not focus on the very basics of programming: there are better sources for that. Rather, we will use the development of some tools in the acoupy project as an example to illustrate how scientists and engineers can attempt being better developers. Bear in mind though that this article is not written by software developers. If you are a software developer, and see some bad advice in here, please let us know!

These articles will be written while the packages we refer to are also being written. Hence, at some point in the future the software and the snippets in the article might differ. We will try to keep things as consistent as possible.

We will use Python as a language. This because it is the most straightforward language to use with FEniCS. But also because it offers many features that make it very suitable to write scientific code.

Writing Good Software

Software doesn’t simply need to run and provide the correct answer. Humans will also need to be able to understand it. Most often, the human that needs to understanding the most will be the developer(s). Understanding what you written, and why, can become quite a challenge even after a few weeks. If understanding the code is hard, then fixing it or extending it becomes much harder. In other words:

If it doesn’t look good it doesn’t work.

Code that is well written looks good. You might have found yourself in situations in which you typed your code fast, to get to a result. It is messy, but it works. Except that it doesn’t work. As soon as you have a bug, or need to change something (maybe update to reflect a change in a dependency…) then it is as good as broken code. You struggle to read it, and hence to fix it. Writing code which is bug free and completely future proof is impossible. Hence, you need to write it in a way that allows you to deal with problems the easiest way. Only code that looks good has a chance at being that. All the rest doesn’t.

Good Looking Code

There a few general rules that define good looking code. However, these are rules, not laws. We will list them below, with the understanding that one day you might find compelling reasons to not apply some of those in certain cases.

Style Guidelines

Refer to style guidelines for your language of choice. For Python this will be PEP 8. Style guidelines exist for every language and you should read them as they have specific reasons to exist. Keeping them all in mind, and using them, might seem like a huge task. But a good IDE will make things much easier. The IDE can highlight where your style is not up to standards and suggest fixes. A good IDE for Python that will check your style is PyCharm.

Use Descriptive Variable Names

You wrote some acoustics simulation code last year. Now it has a bug. You dig into the code and read:

T = 0.161 * V / (S * a)

What is T? What is 0.161? Does it have some physical dimension? What are all the other variable names? You will have to go through your code to the place of definition of each one of the variables to figure it out. Which is time consuming and can easily result in confusion.

This would be better:

# Refer to EQ. 12.3.4 of Fundamentals of Acoustics, Fourth Edition, page 336.
reverberation_time = (
    constants.get_metric_sabine_constant() * room_volume /
    (room_surface * average_sabine_absorptivity)
    )

Now we know what of each piece is just by reading the line. The magic number is returned by a function in a module, so that we know what it is. We also know it is metric, so all the variables should be SI units. Being the value provided by a function we change the function definition if we want to change accuracy, or having an algorithm delivering the value. This will be reflected everywhere the function is used. The comment tells where to grab the equation, so we can look the theory up. All of this makes the line clear, and we are ready to debug why it is giving the wrong result. For example, are we using imperial units for volume and surface?

On the downside the expression got longer. That is not quite a downside, as we gain way more in clarity. It took marginally more time to write but takes about the same time to read. Also, we used ( ) to expand it on multiple lines in a readable way: numerator and denominator are clearly separated.

Comments

Use good comments. Good comments explain why things are done. To learn effective use of comments you can read this excellent tutorial. Even though that is C++ the ideas apply just fine to Pyhton (or any other language). Just follow your language style guidelines when typing them.

Leverage Packaging Infrastructure

If your language, like Python, supports packages then use them. You will be able to create your own packages and install them anywhere. For an overview about crating Python packages and hosting them on Git repositories see Hosting Python Packages on Git Repositories. By crating packages you crate a set of tools that you can use anywhere you need!

Leverage Documentation Infrastructure

Your language of choice might offer quite a lot of infrastructure when it comes to write inline documentation. Not only that, it will most likely have tools to generate online documentation too. For Python one can use many different formats for inline documentation and different HTML documentation generators. Hosting Python Packages on Git Repositories goes over some of the basics.

Proper documentation will help you using your code once you written it. If your code is disseminated it will help others to make use of it.

Use Type Hints

In Python we can use type hints. It is a very good idea to use them. They mainly serve to document your code and make functions and classes clearer. However, they do serve also quite practical purposes. For example, PyCharm will flag potential issues and auto-complete much easier when they are present. We could write a function for Sabine reverberation formula as follows:

def sabine_reverberation_time(
    room_volume: float,
    room_surface: float,
    average_sabine_absorptivity: float
) -> float:
    return (
        constants.get_metric_sabine_constant() * room_volume /
        (room_surface * average_sabine_absorptivity)
    )

Now we know what kind of objects the function expects, and what object it returns. This will also appear in inline documentation.

The acoupy Project

acoupy aims to be a small collection of Python utilities to create acoustics simulations. Each utility aims to be minimal and serving a single purpose. The project doesn’t aim to be a full fledged toolbox, with proper development timeline and support. However, it can serve as an example on how to develop scientific packages.

At the moment of writing there are only two packages:

• acoupy_meshutil;
• acoupy_helmholtz;

Both of these have started development very recently. Their purposes are outlined below.

acoupy_meshutil

Creating meshes with FEniCS itself is really only suitable for simple problems. For most problems you will want to create meshes by using something like Gmsh or Salome Platform. The acoupy_meshutil package is meant to read a mesh and convert it to FEniCS format. At the time of writing it only supports MED files (Salome). The code is able to recognise entities (surfaces, bodies, etc…) and allows their uses within FEniCS. Plans for the future are:

• Add MSH file support (Gmsh);

acoupy_helmholtz

This package implements the model outlined in the Intro to FEniCS - Part 2 episode. At the time of writing it only supports uniform boundary conditions. Plans for the future are:

• Add Adiabatic Absorbers support;
• Add source terms support;
• Add nonuniform boundary conditions support;

Using These Packages

These packages can be easily installed with pip. By using these packages we can setup a Pulsating Sphere simulation as follows:

from acoupy_meshutil import readers
from acoupy_helmholtz import solver
import pathlib
import fenics
import numpy as np


reader = readers.Reader(
    file_format=readers.Reader.Format.MED,
    dimension=readers.Reader.Dimension.D_3D
)

success, exception = reader.read(pathlib.Path(__name__).parent.joinpath('mid_spherecut_m.med'))
if not success:
    raise exception
    
mesh, markers, sub_domains = reader.get_mesh()
cell_tags = reader.get_cell_tags()

dx_volumes = fenics.Measure('dx', domain=mesh, subdomain_data=markers)
ds_surfaces = fenics.Measure('ds', domain=mesh, subdomain_data=sub_domains[0])

frequency = 1000.0
medium = solver.MediumSettings()
outer_radius = 0.343
source_velocity = (1.0 + 1.0j) / np.sqrt(2.0)

wave_number = 2.0 * np.pi * frequency / medium.speed_of_sound
outer_impedance = (
    medium.density * medium.speed_of_sound * (1.0j * wave_number * outer_radius) / 
    (1.0 + 1.0j * wave_number * outer_radius)
)

solution = solver.solve(
    mesh=mesh,
    simulation_settings=solver.SimulationSettings(
        frequency=frequency, 
        medium_settings=medium
    ),
    main_domain_settings=solver.MainDomainSettings(measure=dx_volumes),
    element_settings=solver.ElementSettings(degree=2),
    uniform_neumann_bcs=solver.UniformNeumannBCS(
        [
            solver.UniformNeumannBC(
                velocity=source_velocity,
                measure=ds_surfaces(1)
            )
        ]
    ),
    uniform_robin_bcs=solver.UniformRobinBCS(
        [
            solver.UniformRobinBC(
                impedance=outer_impedance,
                velocity=0.0,
                measure=ds_surfaces(0)
            )
        ]
    ),
)

Where mid_spherecut_m.med is a MED mesh file we prepared beforehand. There is a lot to unpack here, so we will cover this (and more) in the next episodes.

Conclusion

Developing good simulations is not just a matter of knowing how to use software. It is also a matter of knowing how to write good software. In this episode we explained why. We also provided some tips and links to useful information. We will use the acoupy project as a mean to illustrate these concepts. In the next episodes we will explain why these packages are written the way they are. We will also show them in action by using the Pulsating Sphere benchmark problem.

Featured Picture

This article featured picture is from The Python Logo.

License Information

This work is licensed under a Creative Commons Attribution 4.0 International License.

• Python