published or not. The documents may come

The roles of code in computational science

It's a safe bet that most readers of CiSE have written some code during their professional life. Many of us write code regularly as part of our scientific activity, perhaps even as a full-time job. But even though we write, and use, more and more code, we rarely think about the roles that this code will have in our research, in our publications, and ultimately in the scientific record.

In this article, I will outline some frequent roles of code in computational science. These roles are not exclusive, and in fact it is common for a piece of code to have several roles, at the same time or as an evolution over time. Thinking about these roles, ideally before starting to write the code, is a good habit to develop.

The most obvious role of code is its use as a tool for doing computations. In fact, it is this role that sets code apart from other types of information. In its tool role, code is evaluated using criteria such as suitability for a task, robustness, flexibility, performance, ease of use, etc. Since this role should be familiar to CiSE readers, I won't say much more about it and concentrate on the other, less obvious roles.

Code as a scientific notation

Another important family of roles can be summarized as executable expressions of scientific information. In these roles, code can be considered a scientific notation. Criteria such as clarity of expression and exposition are therefore important, just like for any other form of scientific communication. However, executability matters as well, otherwise we would stick to more traditional notations such as plain language or mathematical formulas. Being able to execute the code is a proof of completeness. In an article, you can forget to define a quantity without anybody noticing for a long time. In a computer program, such an omission yields an error message. Executability is also a proof of absence of certain types of mistakes: the mistakes that a compiler can identify, and the mistakes that cause a program to fail during its execution. I call this a family of roles because there are very different kinds of scientific information that can be expressed by code. Each kind requires its specific style, both for the code itself and for the documentation that accompanies it. In this respect, code is not different from traditional scientific writing: a journal article and a textbook for students have different goals and are therefore written in very different styles. This analogy between code and scientific writing turns out to be useful for identifying different types of executable scientific knowledge.

A paper in a scientific journal reports on a scientific study, explaining the context and motivation, the exact system being studied, and the methods that were applied. Results are then presented and discussed. The computational equivalent is a detailed record of all computations that were performed, including the software, the input data, and the results.

Note that this detailed record is not just code, but code in its scientific context. This is true in general for the role of code as scientific notation. In fact, it is true as well for other scientific notations. A mathematical formula, or even a sentence in plain English, does not convey much meaning on its own. It acquires meaning only in the context of other formulas and sentences. Similarly, the code that was run in a computational research project makes sense only in the context of an overall description of that project. This insight has led to the development of computational notebooks, pioneered by Mathematica (http://www.wolfram.com/mathematica/) and more recently implemented for several programming languages by the Jupyter project (http://jupyter.org/). Notebooks have become popular because they make it possible to combine the computational protocol used in a study with a scientific narrative into a single document. The reader of such a notebook can not only read the narrative, but also run the code, inspect intermediate results, and change parameters to see their impact on the results. This ability to work with the code makes it easier to understand the scientific narrative and vice versa.

If journal articles were the only type of scientific document, we would have drowned in the ocean of millions of individual scientific studies long ago. What makes it possible to keep an overview of progress in a specific field is condensed summaries, which take the form of review articles. They provide a snapshot of the state of the art in a field that still undergoes rapid development. In the world of code, the equivalent of a review article is a library or program package that implements state-of-the-art algorithms with the goal of making them easily accessible to other researchers. Like the readers of a review article, the users of such a library are mainly experienced practitioners who can be expected to know the ideas behind the algorithms. However, precise and up-to-date technical documentation is essential: users must be able to figure out easily what exactly the code does, and in particular what its limitations are.

Beyond review articles, there is a continuous consolidation process of scientific information whose final stage is textbook knowledge. Knowledge that is presented in textbooks is information that is well known and well understood by research practitioners from many specialty fields. The code equivalent of textbook knowledge is a library or program package that implements well-known and widely used algorithms in such a way that they can safely be used as black-box ingredients in research code. This requires particular attention to robustness. The code should behave reasonably even when used in an unforeseen way. Documentation must be written with non-expert users in mind, and should address both technical conditions and appropriate scientific contexts for the use of the code.

Unfortunately, we have not yet found a good way to present libraries and program packages to human readers in their scientific context. Ideally, readers should see the code embedded into a discussion of the methods that are implemented. This discussion should be illustrated by example applications, and test cases should point out assumptions about the inputs, the handling of edge cases, and other subtleties. Literate programming [Knuth1984] goes a long way towards this goal. Its basic premise is that software should be written as an essay for human readers, with markup that permits extracting pure code files for execution.

Note that literate programming is more complex than the computational notebooks I mentioned above. A notebook describes a linear sequence of computational steps with is inputs and outputs alongside a scientific narrative. In contrast, algorithms and computational methods are not linear, and are designed to work with a wide range of inputs. That is perhaps one reason why literate programming hasn't been adopted yet by computational scientists. Another reason is the unsatisfactory state of support tools, which are incomplete and incompatible with each other. But the main reason is that the scientific community has only recently acknowledged the importance of the scientific notation role played by code.

The notation roles that I have discussed up to here deal with scientific models and methods, which have occupied the central stage of research for a few centuries. Today, data plays an ever more important role, and datasets are increasingly recognized as scientific publications in their own right. Most published datasets are raw or processed observational data, but it can also be of interest to publish simulation output if re-generating it is particularly costly or difficult. Whatever the nature of the data, every electronic dataset has been created with the help of software, and that software's source code contains valuable information for interpreting the data. Ideally, the syntax and semantics of datasets should be well-defined and documented (a topic I have written about earlier [Hinsen2012]), making it less important to have access to the code that generated the published files. However, just like descriptions of algorithms in plain language tend to be incomplete, descriptions of data models are not always sufficient in practice. In the case of processed data, a detailed record of the processing steps in the form of the code that was actually run is also very helpful for interpretation. This role of code overlaps with the role of documenting a scientific study. The main difference is a focus on the process in one case compared to a focus on the results in the other case. Another difference is that code documenting datasets should ideally be attached to the datasets as metadata and not require fetching another document from a different place.

Code as the shared asset of a scientific community

A final role of code that deserves some discussion is its role in the structuring of scientific communities. Major program packages or libraries, but also more widely used infrastructure tools such a programming languages, have often been the nucleus for the formation of sizable user and developer communities which, in the case of programming languages and infrastructure libraries, transcend the traditional scientific disciplines. The social impact of code has been much larger than for any other technical aspect of doing research. I see the main reason for this in the enormous complexity of software. Keep in mind that when you run a 10-line Python script on your computer, you are really deploying a stack of software items that includes the Python interpreter, the libraries installed along with it, and an operating system such as Linux. Moreover, getting this software stack working on your computer has required other software that you may no longer see, such as a C compiler. The development, maintenance, and deployment of this enormous amount of code cannot be handled by any individual, nor even a team of people working in close contact.

Consider a computational scientist working on the development of domainspecific code. This task requires a good understanding of the level below in the stack -a programming language, various libraries, and various tools such as compilers. Acquiring a useful level of competence with all this is an important effort, so few people can become proficient in more than one or two such environments. The programming environment thus becomes an important aspect of one's work, and discussing this work is much easier with people who share the same background. Collaborating on code development imposes even stricter compatibility: all developers must be proficient using almost identical environments. In comparison, traditional scientific instruments are simple. A scientist familiar with the principles of electron microscopy can rapidly learn to use any particular electron microscope, and has no difficulty to talk to or collanorate with colleagues who use different models. Likewise, a theoretician can easily read articles written by scientists from a different school that favors different mathematical notations or uses different definitions.

In practice this means that the choice of a program package or a programming language implies the choice of a community one will interact with, and maybe become a part of. But every community also has shared values and attitudes that are not directly related to its shared technologies. One community may be conservative, valuing long-term compatibility over progress, whereas another one may like its code on the bleeding edge. Moving even further away from technical aspects, communities can be more or less welcoming to newcomers, and more or less open to new ideas. At the extreme, communities can have sexist, racist, or elitist attitudes that can be serious obstacles to productive research work.

The social roles of code are likely to gain in importance as code itself is more and more recognized as important in research. At this time, publishing code has become common though not yet universal. But unlike journal articles, code is not reviewed in the publication process. One obstacle to independent code reviewing is that a potential reviewer must have both the right scientific domain knowledge and solid technical competence in the development technologies employed by the code authors. If reviewing code becomes common, this will likely create social pressure towards standardization of such technologies inside each scientific domain.

Tensions between the roles

The different roles that a piece of code has are not independent. On the contrary, they can easily be in conflict, requiring scientific software developers to find a workable compromise. This is why I recommend to consider the roles of one's code before starting to write it. Compromises tend to be better when they are result of conscious decisions rather than accidents of development.

An often discussed tension is the one between performance in the tool role and clarity in the notation role. Performance optimization tends to make code less understandable. In the extreme, for example when parallelizing code at a low level using MPI, it can become very hard to figure out what the optimized code does exactly. It is possible that better notations, i.e. new programming languages designed from the start to act like scientific notations, can alleviate this tension, but for the moment this is just wishful thinking.

Another important conflict exists between the generality that is desirable for the tool and social roles on one hand, and the requirement for simplicity in the notation role on the other hand. Tools applicable to a wide range of problems are advantageous for their users, who can get away without learning too many tools. They are also advantageous for their developer communities, because these communities can recruit members from a wider base and thus share the development work among more people. However, general tools lead to larger and more complex code bases that are usually more difficult to understand than simple code written specifically to solve one problem.

A related tension stems from the reliance on other code, in particular libraries. It has been discussed in detail by Matt Turk in this department [Turk2015]. From the point of view of the tool role, using existing library code, assuming it is of good quality, is usually an advantage. In the social role, code with library dependencies creates asymmetric links between two communities, which has good and bad sides. From the notation point of view, dependencies on libraries implementing well-known techniques is an advantage, as readers are likely to be already familiar with them. But dependencies on new or little known libraries make understanding the code more difficult for the average reader because there is more code to read -in fact, this is a special case of the generality dilemma.

Interactions between different categories of code

Another aspect to consider in writing code that fulfills its intended roles well is that no piece of scientific code exists in isolation. The computational protocol of a specific study relies on program packages and libraries of varying degrees of maturity, and each program package and library adds dependencies of its own. From the tool perspective, the computational protocol has almost no value, because its work is done. It may serve as an inspiration for the development of other tools, but it is unlikely to be reused in its entirety. As a consequence, there is no point in doing any maintenance. In contrast, the domain-specific or general libraries it depends on will be used by larger communities for a longer time. They will be maintained and extended, but that also means that they undergo permanent change.

This poses a problem from the notation perspective, which requires that a published computational protocol should remain reproducible for as long as the study itself remains of scientific interest, which can easily be several decades. This is the core issue of the reproducibility problem, which has come up repeatedly in CiSE, in particular in two theme issues (January 2009 and July 2012). In the past, it is the tool perspective that has dominated scientific software, with the result that individual computational protocols are effectively built on quicksand foundations consisting of tools that change through continuous improvement. The increasing emphasis on reproducibility, and the ongoing tendency to publish scientific software, are thus likely to create pressure towards stability at least in widely used libraries and program packages.

A final aspect to consider is the evolution of roles in time. From the notation perspective, it is desirable to write new code as the status of a computational method advanced from innovative through established to well-known. The requirements for these respective roles are too different to be assumed by a single code base. The analogy with scientific writing is useful again: it would be ludicrous to compose a review article from pasted-together sections of original research articles. But scientific code is still rarely written with its notation role in mind. From the tool point of view, reusing working code may well seem preferable, and that is also the way of least immediate effort.

Readers waiting for a final coherent message from this article will be disappointed: I cannot offer any simple conclusion. The only advice I can offer is to take all the issues I have outlined into account when starting a coding project. Consider why you are writing the code right now, and for how long you expect to use it. Consider who else might use the code in the foreseeable future, and who might collaborate with you on its future development. Do you plan to publish it? This might actually not be a choice for much longer, so better prepare yourself for answering "yes".

For any code meant to be published, the notation roles should be explored carefully. As with any form of scientific communication, you have to write for your audience. For a code base, that audience consists of the users of the code, but also of the audience that its users are addressing with their own work. If you want me to use your code, you have to make it easy for me to explain to my readers what your code does.

The social role of code is particularly important to consider when building a community around a software package is your explicit goal. But even if you write small and specialized software that is unlikely to ever become widely used, it is worth thinking about how your code will be judged by the user community of your development and infrastructure tools, if only because in the not too remote future this judgment might become part of the scientific reviewing process.

Adopting the best practices of this community is the equivalent of making an effort to write your papers in good English -something that I haven't heard anyone argue against yet.


published or not. The documents may come

Domain-specific languages in scientific computing

If you have been following developments in software engineering over the recent years, you have probably noticed that the term DSL (domain-specific language) has become a minor buzzword in that field. You may have concluded that this is a hot new idea that is certainly not ready for application in real life. What I want to show in this article is that computational scientists (and others) have been using DSLs for decades and will continue to do so. What is new is not DSLs per se, but the name and the attention given to them.

As the name says, a DSL is first of all a language. In the context of computing, this means a formal language, i.e. a language that can be processed mechanically by algorithms. The definition of a formal language consists of two parts: syntax defines which byte sequences are valid elements in the language, and semantics assigns a meaning to each valid element. For example, the syntax of HTML (HyperText Markup Language) defines start tags (e.g. <h1>) and end tags (e.g. </h1>), and prescribes that they must occur as matching pairs, whereas HTML semantics says that the text between <h1> and </h1> is a top-level heading.

The qualifier "domain-specific" indicates that a language was designed specifically for a restricted application domain, as opposed to general-purpose or broad-spectrum languages that aim at being useful in many application domains. HTML is a DSL, because it was designed to represent Web documents and nothing else. Python is an example of a general-purpose language, because you can encode almost any kind of digital information in Python.

Given that there are general-purpose languages like Python, what's the point of more limited languages like HTML? Why don't we just store Web documents as Python scripts, using the ElementTree module for example? There are several good reasons. To start with, the same Web document would be longer and more complex when written in Python than it is in HTML. Web documents would thus be more difficult to read and write for human authors. Moreover, a Web browser would have to know the entire Python language, and come up with some way to detect if the Python script presented to it really is a Web document. It can't just run the script and see what happens, because the script might be malicious and delete all the files on the computer. From a security point of view, a limited language that doesn't allow malicious operations is a clear advantage. And for a language intended to be processed by many application programs, simplicity makes everybody's life easier.

As this example shows, a DSL is (1) more than, (2) less than, and (3) different from a general-purpose language. It is more in providing domain-specific functionality that a general-purpose language would have to acquire via a library. It is less in not providing superfluous functionality that makes the information more difficult to process and at worst could become dangerous. It is different from a general-purpose language in being optimized for a narrow range of information types right from the start, allowing a more concise representation.

To end the theoretical part of my discussion of DSLs, let me add that the distinction between "domain-specific" and "general-purpose" is more one of intended use than of technical characteristics. As an example, consider the TeX language that many scientists use for writing their publications. Its intended purpose is being a domain-specific language for typesetting, but technically it has all the required properties for a general-purpose language, in particular Turing-completeness. TeX can therefore be (ab-)used for other applications, such as the simulating a Mars rover (http://sdh33b.blogspot.fr/2008/07/icfpcontest-2008.html). Exploring the limits of a language has its merits beyond the intellectual challenge of winning a contest, in particular when it comes to detecting security problems. However, I still consider TeX a DSL rather than a general-purpose language, in spite of the possibility of using it for Mars rover simulations.

The two examples I have given until now, HTML and TeX, are both related to electronic documents, and that application domain has many more DSLs to offer. HTML is complemented by other Web-standard DSLs such as MathML for formulas or SVG for vector graphics. More recently, lightweight markup languages such as Markdown or reStructuredText have become very popular. And at the printing-press end of the publishing chain, we find DSLs such as PostScript and PDF.

At this point I can almost hear protesting murmurs. PDF is not a language, it's a file format! But then, isn't HTML a file format as well? Yes, if you prefer. There is in fact a significant overlap between file formats and DSLs. The two terms represent different points of view rather than distinct concepts. The term "file format" emphasizes the technical aspect of data management: how to store data in a file such that a program can read it. The term "language" emphasizes the human-computer interfacing aspect: representing information in such a way that is makes sense to both computers and their users. Purely machine-oriented file formats, and in particular binary formats, would not be considered DSLs, and therefore PDF is indeed a borderline case (it's a container file format inside which page contents are defined in PostScript, which is a DSL by any reasonable definition). But every file format meant to be read or written by humans, be they software developers or users, qualifies as a DSL.

In scientific computing, we find two dominant approaches to the humancomputer interface, not considering graphical users interfaces because they are used for a very different type of interaction. The traditional one is to write application programs, usually in a compiled language such as Fortran or C++, that have a broad but fixed functionality. Such programs tend to require a large number of parameters and input files, which they read from a file typically referred to as an "input", "control" or "configuration file". Since these files are written by humans, their formats should be treated as DSLs. A more recent alternative approach is to create a collection of interoperating libraries, and let users write scripts that combine items from these libraries for doing a specific computation. The scientific Python ecosystem is probably the best-known example for this approach. The problem-specific information is contained in a script, i.e. it is expressed in a general-purpose language rather than a DSL. I will illustrate these two approaches by a toy example: the computation of the electrostatic potential of a set of point charges. A traditional program would require an input file looking like this:

0. 0. 0. 1. 0. 1. 0. -1. 1. 0. 0. -1. 1. 1. 0. 1.
The program's user manual would explain that each line describes one point charge by four floating-point numbers separated by spaces, the first three numbers being the Cartesian coordinates in nanometers, and the last number being the charge in elementary charge units. Without this explanation, the file would be impossible to interpret for a human reader. Even for a reader who knows that the file describes four point charges, it is not obvious if they are stored by line or by column, nor if the charge comes first or last, nor what the units are.

Next, let's look at a Python script using a hypothetical library called ElectrostaticPotential: from ElectrostaticPotential import ChargeSystem import numpy as np charges = [(np.array([0., 0., 0.]), 1.), (np.array([0., 1., 0.]), -1.), (np.array([1., 0., 0.]), -1.), (np.array([1., 1., 0.]), 1.)] energy = ChargeSystem(charges).potential_energy() print("The energy is", energy)

In this case it would be the library's user manual that explains how to store the positions and charges in Python data structures, and in which units the values must be given. For the reader of the file, the Python code provides much useful context. Contrary to the input file containing nothing but 16 numbers, the script explicitly refers to charges and potentials. A reader with a basic knowledge of electrostatics and the Python language could easily identify the position and charge values from this context, but would remain in the dark concerning their units.

The sixteen-number program input file clearly suffers from a lack of explicit context (the type of information, including the units) and from a lack of structure (which values are the charges?) that would permit basic error checking. The Python script adds the structure and some of the context, but still lacks unit information. On the other hand, the Python script also contains too much information for the purpose of describing a point charge system: it refers to a specific library (ElectrostaticPotential) and to a specific computation (potential en-ergy). Moreover, it also contains scientifically irrelevant technical details, such as the use of three different sequence data structures (lists, tuples, arrays). If I want to compute something else, say the total charge of my system, using a different library, I have to write another script from scratch, probably using somewhat different data structures. I end up having two copies of my point charge system definition that I may have to keep in sync for a while as my research project evolves.

So what would a proper DSL for point charge systems look like? A good starting point is to consider how the information would be presented in a scientific article. My choice would be a table: This table contains exactly the information we want to store, neither more nor less. This could be translated straightforwardly into a machine-readable file:

| position [nm] | charge [e] | |---------------|------------| | 0. 0. 0. | 1. | | 0. 1. 0. | -1. | | 1. 0. 0. | -1. | | 1. 1. 0. | 1. |
To make sure that this table syntax is a valid formal language, I would have to show that it can be parsed unambiguously by a computer program. The proof would consist of actually writing a parser, but I can spare myself that effort because others have already done it: I have borrowed the table syntax from an extension to the lightweight markup language Markdown, which is parsed by the popular file format converter pandoc (http://www.pandoc.org/), among others.

Compared to the original input file, this point-charge-systems DSL has the advantage of being immediately comprehensible to a human reader, because all information is explicit. Note that the presence of a unit indication does not imply that every program reading such a file must be able to handle arbitrary length and charge units. It might well report an error if the units are different from what it expects. In fact, the definition of the DSL (which I have not given) might well prescribe that the input table must have exactly two columns with exactly the headings given in the example, making the first line completely redundant. It's up to the DSL designer to choose the right compromise between generality and simplicity. The main difference between the DSL and file format point of view is that the former can lead to requiring redundant information for the benefit of human readers, and require software to verify it as part of error checking.

The main obstacle to human-friendly DSLs is the effort required to implement them. Computational scientists are not specialists in parser development, so it's not something they are particularly eager to do. There are, however, various libraries and tools, called parser generators, that make the task easier. But even if you prefer to avoid writing parsers altogether (like myself, as I will happily admit), there are useful compromises that are much easier to implement while retaining many of the advantages of human-optimized DSLs. The general idea is to pick a syntax for generic data structures and express your DSL in terms of them. You can then use an off-the-shelf parser for the data structure syntax, and don't have to write your own. As a bonus, you may find that your favorite text editor already supports such data structure syntax.

Let me show you some concrete examples. A popular generic data structure syntax is YAML (http://www.yaml.org/), which calls itself a data serialization language. One way to encode our point charge system in YAML is columns:

-position: nm -charge: e records:

- (point-charges ((position nm) (charge e)) ((0. 0. 0.) 1.) ((0. 1. 0.) -1.) ((1. 0. 0.) -1.)

One of the big ideas of Lisp is to use the flexible data structure notation of s-expressions not only for data, but also for code. This feels weird at first to anyone coming to Lisp from other languages, but to those who persist until they get used to it, the syntactical uniformity starts to show its advantages. One of them is the ease of adding small DSLs to one's programs using Lisp macros. If macros were invented today, they would perhaps be called DSL compilers, because that's exactly what they are. Note however that these are so-called embedded DSLs, because they are used inside a program, not in separate files. Python programmers can try this out for themselves using Hy (http://docs.hylang.org/), which is an alternative s-expression syntax for Python. Hy can also be considered a Pythonic dialect of Lisp, because it supports macros. Using Hy macros, you can develop DSLs that are sub-or supersets of Python, or both, with relatively little effort, as long as your DSL syntax fits into the overall s-expression syntax.

In case you end up addicted to DSL development, your ultimate drug is the Racket language (http://www.racket-lang.org/), which claims to be "the world's first ecosystem for developing and deploying new languages." Racket also has a Lisp heritage, and favors s-expressions for everything, but it also lets you plug in your own parser for your own syntax if you prefer. With such a languagedevelopment toolkit, designing and implementing DSLs requires no bigger effort than designing and implementing graphical user interfaces (GUIs). Note that I am not saying that either one is trivial, but both are within the reach of motivated software developers. DSLs are a better match than GUIs for many scientific computing tasks, so I hope scientific software developers will explore this option more intensively than they did in the past. Some examples from the Racket ecosystem worth looking at to appreciate the potential of DSLs are Scribble (https://docs.racket-lang.org/scribble/), a DSL for electronic documents, Slideshow (https://docs.racket-lang.org/slideshow/index.html?q=slideshow), a DSL for slide-based presentations, or Video ((http://lang.video/), a DSL for video editing.

The one guiding principle to keep in mind when designing DSLs is to focus on the human user's perspective, starting with your own one. Ask yourself how you would ideally write down your scientific data or computation for use by a computer. Assume first that there are no software-related constraints, as if you were taking notes for sharing with a colleague. Then try to express the same information using generic data structures, for example in YAML or in s-expressions. Finally, if you believe that moving to a nicer syntax is worth the additional effort, write your own parser.

Konrad Hinsen is a researcher at the Centre de Biophysique Moléculaire in Orléans and at the Synchrotron SOLEIL in Saint Aubin. His research interests include protein structure and dynamics and scientific computing. Hinsen has a PhD in theoretical physics from RWTH Aachen University. Contact him at konrad.hinsen@cnrs.fr.

  Parsing this with the PyYAML library (https://github.com/yaml/pyyaml) returns a simple Python data structure that any Python programmer should know how to process: {'columns': [{'position': 'nm'}, {'charge': 'e'}], 'records: [[[0.0, 0.0, 0.0], 1.0], [[0.0, 1.0, 0.0], -1.0], [[1.0, 0.0, 0.0], -1.0], [[1.0, 1.0, 0.0], 1.0]]}Another generic data structure syntax is known as s-expressions. It is probably the oldest one, having been developed originally for the Lisp programming language in the 1950s. A possible encoding of our point charge system in terms of s-expressions looks like this:


Isabelle Raspo

Stéphane Viazzo

Anthony Randriamampianina

Gabriel Meletti

Uwe Harlander

Torsten Seelig

Andreas Krebs

Low Mach

Low Mach Number Modeling of Stratorotational Instability in a water-filled Taylor-Couette cavity

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

The Stratorotational Instability (SRI) is a purely hydrodynamic instability with distinctive local features occuring in centrifugally stable flows. It is suspected to contribute to the outward angular momentum transport in accretion discs, which is a central problem of planet formation. Only turbulence can achieve such a large angular momentum transport. SRI can be studied from specifically designed laboratory experiments and numerical simulations in an axially-stratified Taylor-Couette setup [START_REF] Gellert | Stratorotational instability in Taylor-Couette flow heated from above[END_REF][START_REF] Bars | Experimental analysis of the stratorotational instability in a cylindrical Couette flow[END_REF]. Recent studies [START_REF] Ibanez | Observations of the stratorotational instability in rotating concentric cylinders[END_REF][START_REF] Rüdiger | The stratorotational instability of Taylor-Couette flows with moderate Reynolds numbers[END_REF] mentioned that, at fixed values of the Froude number, F r = Ω i /N (with N the buoyancy frequency and Ω i the angular velocity of the inner cylinder), and of the angular velocity ratio µ = Ω o /Ω i (with Ω o the angular velocity of the outer cylinder), there exists an upper limit of the Reynolds number Re = Ω i R i (R o -R i )/ν for the occurrence of SRI. These authors also showed that the value of Re where the flow becomes stable again increases as F r decreases. Therefore, at high Reynolds number values relevant to transition to turbulence regimes, SRI can only be obtained at small Froude number values and, consequently, for large temperature difference ∆T between the top and bottom endplates. When liquid water is used in simulations, the validity of the Boussinesq approximation may be questioned under such conditions. Indeed, Shalybkov and Rüdiger [START_REF] Shalybkov | Stability of density-stratified viscous Taylor-Couette flows[END_REF] suspected non-Boussinesq effects to be responsible for large discrepancies between numerical and experimental results for strong stratifications, e.g small values of F r. Medale and Haddad [START_REF] Medale | A 3D Low Mach Number model for high performance computations in natural or mixed convection Newtonian liquid flows[END_REF] reported non-Boussinesq effects in a mixed convection problem using liquid water. The Boussinesq approximation can theoretically be applied with liquid water at 298.15K for temperature differences ∆T less than 5.2K [START_REF] Gray | The validity of the Boussinesq approximation for liquids and gases[END_REF]. We investigated the influence of non-Boussinesq effects on the characteristics of SRI by performing 3D simulations with both Low Mach number and Boussinesq approximations based on spectral methods in the configuration of experiments at the BTU Cottbus-Senftenberg, Germany [START_REF] Rüdiger | The stratorotational instability of Taylor-Couette flows with moderate Reynolds numbers[END_REF].

For the Low Mach number modeling, the specific equation of state proposed for water by Pátek et al. [START_REF] Pátek | Reference correlations for thermophysical properties of liquid water at 0.1MPa[END_REF] is implemented, which allows for an accurate representation of the density and other thermodynamic properties of water, with a mean deviation of 0.01% for the density compared with data provided by NIST. The simulations show that non-Boussinesq effects modify the value of the Froude number for the threshold of occurrence of SRI. Significant changes are also observed for the azimuthal wavenumber, as illustrated in figure 1, and for the fundamental frequency of the instability. The present studies will serve to guide future simulations based on High Performance Computing in order to reach higher values of governing parameters relevant to turbulence regimes [START_REF] Abide | High-Order Compact scheme for High-Performance Computing of stratified rotating flows[END_REF].

Figure 1 :

 1 Figure 1: Iso-contours of radial velocity at mid-height of the cavity for Re = 1000, F r = 0.58 (∆T = 24K) and µ = 0.35. Left: Low Mach number approximation (azimuthal wavenumber m = 1, frequency ω/Ω i =0.54). Right: Boussinesq approximation (m=2, ω/Ω i =1.12).

  Low Mach Number Modeling of Stratorotational Instability in a water-filled Taylor-Couette cavity Isabelle Raspo 1 , Stéphane Viazzo 1 , Anthony Randriamampianina 1 , Gabriel Meletti 2 , Uwe Harlander 2 , Torsten Seelig 2 , Andreas Krebs 2 1 Aix-Marseille Univ, CNRS, Centrale Marseille, M2P2, Marseille, France

2

Department of Aerodynamics and Fluid Mechanics, Brandenburg University of Technology, Cottbus-Senftenberg, Germany

Acknowledgements

This work was granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR-10-EQPX-29-01) of the program "Investissements d'Avenir" supervised by the "Agence Nationale de la Recherche". The authors acknowledge financial support from "Ministère des Affaires étrangères", France, in the framework of a PROCOPE project (PHC 35299PM).


Rodrigue Imad

Guillaume Sicot

Sebastien Houcke

Blind frame synchronization for error correcting codes having a sparse parity check matrix

Keywords: Blind frame synchronization, channel coding, MAP, LDPC codes, convolutional codes

We present in this letter a blind frame synchronization method based on a Maximum A Posteriori probability (MAP) approach. Applied to coded communication systems, this method is based on the calculation of the Log-likelihood Ratios (LLR) of the syndrome obtained according to the parity check matrix of the code. After presenting the proposed method, we compare it to an existing blind synchronization method previously introduced. Simulation results show that the proposed method outperforms the existing one and presents good performance when applied to codes having a sparse parity check matrix such as Low Density Parity Check (LDPC) codes and convolutional codes.

I. INTRODUCTION

C HANNEL coding is considered to be an important step in digital communication systems. At the reception and before decoding the received symbols, one should be able to find the beginning of each codeword so that the decoder functions properly. This step is called Frame Synchronization. Conventional frame synchronization is achieved by adding to the transmitted symbols a synchronization sequence known by the receiver. This sequence is detected at the reception by running the received symbols through a correlator [START_REF] Barker | Group synchronization of binary digital systems[END_REF], [START_REF] Massey | Optimum frame synchronization[END_REF], [START_REF] Lui | Frame synchronization for gaussian channels[END_REF]. Actual codes are very powerful and are able to decode heavily corrupted sequences. In order to maintain good performance in such cases, conventional synchronization methods should increase the length of the inserted synchronization word, which reduces the spectral efficiency of the system. Another solution to improve the system performance is to take advantage of the code structure and consider frame synchronization jointly with the decoding as in [START_REF] Robertson | A generalized frame synchronizer[END_REF], [START_REF] Cassaro | Frame synchronization for coded systems over AWGN channels[END_REF]. Rather than placing the synchronization bits in a separate header, the authors of [START_REF] Howlader | Decoder-assisted frame synchronization for packet transmission[END_REF] propose to place them in a midamble. In [START_REF] Huh | A unified approach to optimum frame synchronization[END_REF], a unified approach to frame synchronization has been described where the observed symbols are modeled as the output of a Markov chain corrupted by additive white Gaussian noise. In this letter, we present a blind method of frame synchronization wherein no additional sequence is added to the coded one. This synchronization method is based on a MAP approach in the sense of minimizing the probability of false frame synchronization of a coded system, given the received signal. It has been first introduced for a Binary Symmetric Channel (BSC) [START_REF] Houcke | Blind frame synchronization for block code[END_REF] and then generalized for a Gaussian Channel [START_REF] Imad | Blind frame synchronization on gaussian channel[END_REF]. This method has been also applied in an Interleaved Division Multiple Access (IDMA) context [START_REF] Houcke | Blind detection for block interleaved division multiple access[END_REF] and it gave promising results. In the literature, we found only one method of blind frame synchronization introduced by Sun and Valenti [START_REF] Sun | Synchronization for capacity-approaching coded communication systems[END_REF], [START_REF] Sun | Optimum frame synchronization for preamble-less packet transmission of turbo codes[END_REF]. The authors claim to have developed an optimal MAP frame synchronizer. However, they have missed an important detail in the expression of the a posteriori probability, which leads to an erroneous MAP estimator. In the remaining of this letter, Sun and Valenti's method will be referred to as the "existing" one and our method as the "proposed" one. This letter is organized as follows. In Section II, we introduce the proposed synchronization method and briefly present the existing one. Simulation results are shown in Section III where the two synchronization methods are applied to LDPC and convolutional codes. Finally, Section IV concludes the work.

II. PROPOSED METHOD OF BLIND FRAME

SYNCHRONIZATION

In this letter, we consider that the transmitter is sending a binary sequence of codewords and is using a Binary Phase Shift Keying (BPSK) modulation. The propagation channel is corrupted by a quasi-static additive white Gaussian noise. For simplicity reasons, we assume that perfect bit synchronization is achieved. Indeed, our proposed synchronization method works well even on low signal to noise ratios and therefore, it is able to deal with small errors in the timing recovery since such errors just increase the bit error rate of the system. For a given code of rate ρ = nc-nr nc , it is often possible to find the corresponding parity check matrix H of size n r ×n c , where n c represents the length of a codeword and n r the number of parity relations. Let b(i) = ±1 be the i th coded and modulated bit to be transmitted. At the reception, the i th received sample is given by r(i) = b(it 0 ) + w(i), where t 0 is an integer representing the shift of the transmitted symbols due to the delay introduced by the propagation channel and w(i) is a white Gaussian noise. The received sequence of N samples can then be written as r = [r(1), . . . , r(N )]. The main target of frame synchronization is to find the position of a codeword in the received sequence. In other words, we have to estimate the delay t 0 that we assume (without loss of generality) to be lower than n c . At a position t on the received sequence (t ∈ [0, n c -1]), we apply a sliding window of length Kn c samples, where K is an integer greater or equal to one. The study of all the blocks contained in the n c possible positions of this sliding window allows us to determine whether a position corresponds to the correct synchronization moment or not. In the remaining of this letter we assume that K = 1. This means that only one block is contained in the sliding window. Our proposed method is based on a MAP approach in the sense of maximizing the probability that a position t corresponds to the correct synchronization moment, given the received samples. In other words, it maximizes the following a posteriori probability:

P r[t/r], t ∈ [0, n c -1]. (1)
In order to confirm that a position t is the correct synchronization position, we should verify that the block transmitted at this position is a valid codeword and that all the other blocks in the remaining (n c -1) tested positions are not.

The easiest way to check whether a block corresponds to a valid codeword or not is by calculating its syndrome which is obtained according to the parity check matrix H of the code. Each element of the syndrome is calculated using one parity check equation defined by one row of H. At a position t of the sliding window, the syndrome can be written as

S t = [S t (1), S t (2), . . . , S t (n r )].
At the synchronized position, the probability of having a verified parity check equation is greater than the one at a nonsynchronized position. This is due to the fact that at the synchronized position, the sliding window contains a block corresponding to a valid codeword.

In the logarithmic domain, these probabilities can be expressed in Log-Likelihood Ratios (LLR). Let ψ(t) be the LLR of the syndrome calculated on position t of the sliding window.

ψ(t) = log P r [S t (1)
, . . . , S t (n r )] = 0 P r [S t (1), . . . , S t (n r )] = 0 = log 1 -P t P t .

(2) For codes having a sparse parity check matrix, the syndrome elements can be assumed independent and therefore:

P t = P r [S t (1), . . . , S t (n r )] = 0 = nr k=1 P r[S t (k) = 0]. (3) Let L(S t (k)) = log Pr[St(k)=1] P r[St(k)=0] be the LLR of the k th syndrome element, P r[S t (k) = 0] = 1 1 + exp L(S t (k) . (4)
Using ( 2),( 3), ( 4) and having

nr k=1 1 + exp L S t (k)
much larger than 1, the LLR of the syndrome can be written as:

ψ(t) = nr k=1 log 1 + exp L S t (k) (5)
and the frame synchronization position is estimated by the one that minimizes the LLR of the syndrome:

t0 = argmin t=0,...,nc-1 {ψ(t)}. (6)
For simplification reasons, (5) can be approximated by:

φ(t) = nr k=1 L S t (k) . (7)
Equation ( 7) is equivalent to define the LLR of the syndrome by:

φ(t) = log P r [S t (1) = 1, . . . , S t (n r ) = 1] P r [S t (1), . . . , S t (n r )] = 0 , (8)
which can be interpreted as the worst case scenario: either all parity check equations are verified or none. Finally, the proposed blind frame synchronizer estimates the frame synchronization position by: t0 = argmin t=0,...,nc-1

{φ(t)}. (9)
Note that the LLR of a syndrome element is given by:

L(S t (k)) = (-1) u k +1 atanh u k j=1 tanh r(t + k j ) 2 , (10)
where r(i) = 2 σ 2 r(i) is the LLR of the i th received sample and σ 2 is the variance of the noise. u k and k j represent the number of ones in the k th row of the parity check matrix of the code and the position of the j th non zero element in this k th row, respectively. According to [START_REF] Hagenauer | Iterative decoding of binary block and convolutional codes[END_REF], an approximation of ( 10) is given by:

L S t (k) = (-1) u k +1 u k j=1 sign r(t+k j ) min j=1,...,u k r(t+k j ) .
(11) Having r(i) and r(i) proportional and by assuming the variance of the noise constant during a transmission, r(t + k j ) can be replaced by r(t + k j ) in [START_REF] Sun | Synchronization for capacity-approaching coded communication systems[END_REF]. In this way, no a priori information about the AWGN channel is required [START_REF] Chen | Near optimum universal belief propagation based decoding of low-density parity check codes[END_REF]. L S t (k) and φ(t) are then proportional to:

L S t (k) = (-1) u k +1 u k j=1 sign r(t + k j ) . min j=1,...,u k r(t + k j ) (12)
and

φ(t) = nr k=1 L S t (k) , (13)
and the frame synchronization position is estimated by:

t0 = argmin t=0,...,nc-1 { φ(t)}. ( 14
)

Comparison with an existing method of blind frame synchronization

Applied in our context, the existing method introduced in [START_REF] Sun | Synchronization for capacity-approaching coded communication systems[END_REF], [START_REF] Sun | Optimum frame synchronization for preamble-less packet transmission of turbo codes[END_REF] aims to maximize the log-likelihood function:

L(t) = - (φ(t)-nrMt 0 ) 2 2κnrMt 0
, where n r M t0 is the mean of φ(t 0 ) and κ is a coefficient introduced to describe the dependence of the syndrome elements. The concept of this method is very close to ours but with one main drawback. In order to estimate the synchronization position, the authors examine only whether φ at a position t follows the same law as the one at the synchronized position. However, they do not verify that the remaining φ(i) i =t corresponding to the other (n c -1) tested positions follow the law of φ at a nonsynchronized position. In other words, comparing the variables φ to the mean of φ(t 0 ) is not enough to guarantee an optimum estimation of t 0 . At a nonsynchronized position, we might have a value of φ closer to the mean of φ(t 0 ) than a correct value corresponding to the synchronized one. Hence, it is obvious that the existing method does not afford an optimum estimation of the correct synchronization position.

III. SIMULATION RESULTS

In order to study the performance of our synchronization method, we estimated the probability of false synchronization. The evaluation of this probability is realized by Monte Carlo simulation. For each configuration, the noise, information bits and the delay of the channel were randomly chosen. Fig. 1 shows a comparison between the proposed synchronization method and the existing one once applied to an LDPC code of length 512 bits, having 5 non zero elements in each row of its parity check matrix. First of all, we can see that using [START_REF] Cassaro | Frame synchronization for coded systems over AWGN channels[END_REF] or [START_REF] Huh | A unified approach to optimum frame synchronization[END_REF] to estimate the frame synchronization position gives almost the same results. Therefore, in the remaining of this letter and for complexity reasons, our proposed method will be based on [START_REF] Huh | A unified approach to optimum frame synchronization[END_REF]. From Fig. 1, it is clear that the proposed method outperforms the existing one. A gain of around 0.5 dB is achieved at a probability of false synchronization equal to 10 -3 . Fig. 2 shows the same comparison when applying these methods to convolutional codes. Two systematic and recursive convolutional codes of generator polynomials (1, 5/7) and (1, 23/35) were tested. These two codes have the same length (512 bits), the same rate (1/2) and their constraint lengths equal to 3 and 5, respectively. Note that the number of non zero elements in the parity check matrix increases with the constraint length of the code. Consequently, a degradation in performance is observed [START_REF] Imad | Blind frame synchronization on gaussian channel[END_REF]. Compared to the existing method, our proposed blind frame synchronization method yields always better results. We also compared our results to the ones obtained with Massey's synchronization approach, which locates a synchronization word periodically embedded in the binary data [START_REF] Massey | Optimum frame synchronization[END_REF]. This comparison is done in Fig. 3 for the convolutional code (1, 5/7). Three different lengths of a synchronization word were tested, L = 15, 25 and 40 bits. We can clearly see that our method converges quicker than the one introduced by Massey. Even for a synchronization word of length L = 40 bits (around 8% of the total length of the code), our method presents better performance in the zone were the convolutional code is generally used (E b /N 0 > 3.5 dB). In order to evaluate the robustness of the proposed synchronization method, we plotted in Fig. 4 the Frame Error Rate (FER) curves obtained by decoding the previous convolutional codes using a MAP decoder, which is applied after achieving the blind frame synchronization step. As we can see, the FER curves for code (1, 5/7) obtained after applying the proposed synchronization method and the ones obtained in the perfect synchronization case are almost the same. A small degradation in performance is observed for code (1, 23/35). However, for E b /N 0 > 3.5 dB, the two curves of FER become one and no degradation in performance is detected.

IV. CONCLUSION

In this letter, we have proposed an optimum blind frame synchronization method based on a MAP approach. Applied to codes having a sparse parity check matrix such as LDPC codes, the proposed synchronization method presents good performance and clearly outperforms an existing blind frame synchronization method. Furthermore, once the proposed synchronization method is applied to convolutional codes, the frame error rates of the system are almost the same as the ones obtained in the perfect synchronization case.

Fig. 1 .

 1 Fig. 1. Comparison between the proposed and existing methods applied to an LDPC code of rate R = 1/2 and length nc = 512 bits.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Comparison between the proposed and existing methods applied to convolutional codes of rate R = 1/2 and length nc = 512 bits.

Fig. 4 .

 4 Fig. 4. Frame error rate of convolutional codes with code rate R = 1/2 and length nc = 512 bits.