#+SETUPFILE: ../packages/org-themes/src/readtheorg_inline/readtheorg_inline.theme
#+STARTUP: org-pretty-entities entitiespretty
#+PROPERTY: header-args :comments no
#+OPTIONS: timestamp:nil
#+OPTIONS: h:3 num:t toc:3
#+TITLE: The pugs user manual
#+OPTIONS: author:nil date:nil
#+OPTIONS: tex:t
#+LANGUAGE: en
#+ATTR_LATEX: :width 4cm
#+HTML_HEAD_EXTRA: <style> pre.src-pugs:before { content: 'pugs'; } </style>
#+HTML_HEAD_EXTRA: <style> pre.src-pugs-error:before { content: 'invalid pugs'; } </style>
#+HTML_HEAD_EXTRA: <style> .footref{ color: #2980b9; font-size: 100%; } </style>
#+LATEX_CLASS_OPTIONS: [10pt]
#+LATEX_HEADER: \usepackage[hmargin=2.5cm,vmargin=1.5cm]{geometry}
#+LATEX_HEADER: \usepackage{ae,lmodern}
#+LATEX_HEADER: \usepackage[OT1]{fontenc}
#+LATEX_HEADER: \usepackage{booktabs}
#+LATEX_HEADER: \usepackage{siunitx}
#+LATEX_HEADER: \usepackage[table]{xcolor}
#+LATEX_HEADER: \usepackage{colortbl}
#+LATEX_HEADER: \hypersetup{linktoc = all, colorlinks = true, linkcolor = blue}
#+LATEX_COMPILER: pdflatex --shell-escape
#+LATEX_HEADER_EXTRA: \usepackage{amsmath}
#+LATEX_HEADER_EXTRA: \usepackage{amsthm}
#+LATEX_HEADER_EXTRA: \usepackage{amssymb}
#+LATEX_HEADER_EXTRA: \usepackage{xcolor}
#+LATEX_HEADER_EXTRA: \usepackage{mdframed}
#+LATEX_HEADER_EXTRA: \usepackage{mathpazo}
#+LATEX_HEADER_EXTRA: \usepackage{inconsolata}
#+LATEX_HEADER_EXTRA: %Patch accsupp to avoid copying line numbers when copying from listing
#+LATEX_HEADER_EXTRA: \usepackage{accsupp}
#+LATEX_HEADER_EXTRA: \newcommand\emptyaccsupp[1]{\BeginAccSupp{ActualText={}}#1\EndAccSupp{}}
#+LATEX_HEADER_EXTRA: \let\theHFancyVerbLine\theFancyVerbLine
#+LATEX_HEADER_EXTRA: \def\theFancyVerbLine{\rmfamily\tiny\emptyaccsupp{\arabic{FancyVerbLine}}}
#+LATEX_HEADER_EXTRA: \BeforeBeginEnvironment{tabular}{\rowcolors[]{2}{orange!5}{orange!10}}
#+LATEX_HEADER_EXTRA: \BeforeBeginEnvironment{minted}{\begin{mdframed}[linecolor=blue,backgroundcolor=blue!10]}
#+LATEX_HEADER_EXTRA: \AfterEndEnvironment{minted}{\end{mdframed}}
#+LATEX_HEADER_EXTRA: \BeforeBeginEnvironment{verbatim}{\begin{mdframed}[linecolor=gray,backgroundcolor=gray!5]}
#+LATEX_HEADER_EXTRA: \AfterEndEnvironment{verbatim}{\end{mdframed}}
#+LATEX_HEADER_EXTRA: \newtheorem{note}{Note}
#+LATEX_HEADER_EXTRA: \BeforeBeginEnvironment{note}{\begin{mdframed}[linecolor=green,backgroundcolor=green!10]}
#+LATEX_HEADER_EXTRA: \AfterEndEnvironment{note}{\end{mdframed}}
#+LATEX_HEADER_EXTRA: \newtheorem{warning}{Warning}
#+LATEX_HEADER_EXTRA: \BeforeBeginEnvironment{warning}{\begin{mdframed}[linecolor=red,backgroundcolor=red!10]}
#+LATEX_HEADER_EXTRA: \AfterEndEnvironment{warning}{\end{mdframed}}
#+LATEX_HEADER_EXTRA: \usemintedstyle{perldoc}
#+LATEX_HEADER_EXTRA: \renewcommand{\MintedPygmentize}{${PUGS_SOURCE_DIR}/tools/pgs-pygments.sh}
#+begin_src latex :exports results
\let\OldTexttt\texttt
\renewcommand{\texttt}[1]{\fcolorbox{gray!50}{gray!5}{\OldTexttt{#1}}}
#+end_src
* Introduction
~pugs~[fn:pugs-def] is a general purpose solver collection built to
approximate solutions of partial differential equations. It is mainly
(but not only) designed to deal with finite-volume methods.
~pugs~ is a parallel software that uses ~MPI~[fn:MPI-def] and multi-threading
for parallelism. Multi-threading is achieved through an encapsulation
of some [[https:github.com/kokkos/kokkos][Kokkos]] mechanisms.
The philosophy of ~pugs~ is to provide "simple" numerical tools that are
assembled together through a high-level language (a DSL[fn:DSL-def]
close to the mathematics) to build more complex solvers. This approach
is inspired by the success of [[http://freefem.org][FreeFEM]], which uses a similar approach.
Before detailing the leading concepts and choices that we have made to
develop ~pugs~, we give a simple example.
** An example
For instance, the following code builds a uniform Cartesian grid of
$]-1,1[^2$ made of $20\times20$ cells and transforms it by displacing the
mesh nodes according to a user-defined vector field $T: \mathbb{R}^2
\to \mathbb{R}^2$.
#+NAME: introduction-example
#+BEGIN_SRC pugs :exports both :results output
import mesh;
import math;
import writer;
let pi:R, pi = acos(-1);
let theta:R^2 -> R, x -> 0.5*pi*(x[0]*x[0]-1)*(x[1]*x[1]-1);
let M:R^2 -> R^2x2, x -> [[cos(theta(x)), -sin(theta(x))],
[sin(theta(x)), cos(theta(x))]];
let T: R^2 -> R^2, x -> x + M(x)*x;
let m:mesh, m = cartesianMesh([-1,-1], [1,1], (20,20));
m = transform(m, T);
write_mesh(gnuplot_writer("transformed"), m);
#+END_SRC
The example is quite easy to read.
- First, some *modules* are loaded: the ~mesh~ module, which contains some
mesh manipulation functions. The ~math~ module provides a set of
classical mathematical functions ($\sin$, $\cos$, ...). The ~writer~
module is used to save meshes or discrete functions to files using
various formats.
- The second block of data defines variables of different types
- ~pi~ is a real value that is initialized by an approximation of $\pi$.
- ~theta~ is the real value function $\theta$ defined by
\begin{equation*}
\theta: \mathbb{R}^2 \to \mathbb{R},\quad\mathbf{x} \mapsto
\frac{\pi}{2} (x_0^2-1)(x_1^2-1)
\end{equation*}
- ~M~ is the $\mathbb{R}^{2\times2}$ matrix field $M$ defined by
\begin{equation*}
M: \mathbb{R}^2 \to \mathbb{R}^{2\times2},\quad\mathbf{x} \mapsto
\begin{pmatrix}
\cos(\theta(\mathbf{x})) & -\sin(\theta(\mathbf{x})) \\
\sin(\theta(\mathbf{x})) & \cos(\theta(\mathbf{x}))
\end{pmatrix}
\end{equation*}
- ~T~ is the vector field $T$ defined by
\begin{equation*}
T: \mathbb{R}^2 \to \mathbb{R}^2, \quad\mathbf{x} \mapsto (I+M(\mathbf{x}))\mathbf{x}
\end{equation*}
- Finally ~m~ is defined as the uniform Cartesian mesh grid of
$]-1,1[^2$. The last argument: ~(20,20)~ sets the number of cells in
each direction: $20$.
- The third block is the application of the transformation $T$ to the
mesh. Observe that if the resulting mesh is stored in the same
variable ~m~, the old one was not modified in the process. It is
important to already have in mind that the ~pugs~ language *does not
allow* the modifications of values of variables of *non-basic*
types. This is discussed in the section [[high-level-types]].
- Finally, the last block consists in saving the obtained mesh in a
~gnuplot~ file. The result is shown in Figure [[fig:intro-example]].
#+NAME: intro-transform-mesh-img
#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/intro-transform-mesh.png")
reset
unset grid
unset border
unset key
unset xtics
unset ytics
set terminal png truecolor enhanced
set size square
plot '<(sed "" $PUGS_SOURCE_DIR/doc/transformed.gnu)' w l
#+END_SRC
#+CAPTION: Obtained transformed mesh
#+NAME: fig:intro-example
#+ATTR_LATEX: :width 0.38\textwidth
#+ATTR_HTML: :width 300px;
#+RESULTS: intro-transform-mesh-img
Even if this first example is very simple, some key aspects can
already be discussed.
- There is no predefined constant in ~pugs~. Here a value is provided
for ~pi~.
- There are two kinds of variable in ~pugs~: variables of basic types
and variables of high-level types. These two kinds of variable behave
almost the same but one must know their differences to understand
better the underlying mechanisms and choices that we made. See
sections [[basic-types]] and [[high-level-types]] for details.
- Also, there are two types of functions: *user-defined* functions and
*builtin* functions. In this example, ~theta~, ~M~ and ~T~ are user-defined
functions. All other functions (~cos~, ~cartesianMesh~,...) are
builtin functions and are generally defined when importing a
module. These functions behave similarly, one should refer to
[[functions]] for details.
** Concepts and design
As it was already stated, ~pugs~ can be viewed as a collection of
numerical methods or utilities that can be assembled together, using a
user friendly language, to build simulation scenarios.
Utilities are tools that are often used in numerical
simulations. Examples of such utilities are mesh manipulations,
definition of initial conditions, post-processing, error
calculations,...
*** A C++ toolbox driven by a user friendly language
Numerical simulation packages are softwares of a particular
kind. Generally, in order to run a calculation, one has to define a
set of data and parameters. This can simply be the definition of a
discretization parameter such as the mesh size. One can also specify
boundary conditions, equations of state, source terms for a specific
model. Choosing a numerical method or even more, setting the model
itself, is common in large codes.
In ~pugs~, all these "parameters" are set through a DSL. Thus, when ~pugs~
is launched, it actually executes a provided script. A ~C++~ function is
associated to each instruction of the script. The ~C++~ components of
~pugs~ are completely unaware one of the others. ~pugs~ interpreter is
responsible for the data flow between the components: it manages the
data transfer between those ~C++~ components and ensures that the
workflow is properly defined.
**** Why?
In this section we motivate the choice of a language and not of a more
standard approach.
***** Data files are evil
There are lots of reasons not to use data files. By data file, we
refer to a set of options that describe physical models, numerical
methods or their settings.
- Data files are not flexible. This implies on the one hand that
application scenarios must be known somehow precisely to reflect
possible option combinations and on the other hand even defining a
specific initial data may require the creation of a new option and
its associated code (in ~C++~ for instance). \\
Usually, the last point is addressed by adding a local interpreter
to evaluate user functions.
- Data files may contain irrelevant information. Actually, it is quite
common to allow to define options that are valid but irrelevant to a
given scenario. This policy can be changed but it is generally not
an easy task and requires more work from the user (which can be a
good thing).
- Generally data files become rapidly obsolete. An option was not the
right one, or its type changed to allow other contexts... This puts
pressure on the user.
- Even worst, option meanings can depend on other
options. Unfortunately, this happens commonly. For instance, a
global option can change implicitly the treatment associated to
another one. This is dangerous since writing or reading the data
file requires an important and up to date knowledge of the code's
internals.
- Another major difficulty when dealing with data files is to check
the compatibility of provided options.
***** Embedded "data files" are not a solution
Using directly the general purpose language (~C~, ~C++~, ~Fortran~,...) used
to write the code can be tempting. It has the advantage that no
particular treatment is necessary to build a parser (to read data
files or scripts), but it presents several drawbacks.
- The first one is probably that it allows too much freedom. While
defining the model and numerical options, the user has generally
access to the whole code and can change almost everything, even
things that should not be changed.
- Again, one can easily have access to irrelevant options and it
requires a great knowledge of the code to find important ones.
- With this in mind, defining a simulation properly can be a difficult
task. For instance, in the early developments of ~pugs~ (when it was
just a raw ~C++~ code) it was tricky to change boundary conditions for
multiphysics problems.
- Another difficulty is related to the fact that code's internal API
is likely to change permanently in a research code. Thus valid
constructions or settings may become rapidly obsolete. In other
words keeping up to date embedded "data files" might be difficult.
- Finally it requires recompilation of pieces of code (which can be
large in some cases) even if one is just changing a simple
parameter.
***** Benefits of a DSL
Actually, an honest analysis cannot conclude that a DSL is the
solution to all problems. However, it offers some advantages.
- First, it allows a fine control on what the user can or cannot
perform. In some sense, it offers a chosen level of flexibility.
- It allows to structure the code in the sense that new developments
have to be designed not only focusing on the result but also on the
way it should be used (its interactions with the scripting language).
- In the same vein, it provides a framework that drives the desired
principle of "do simple things and do them well".
- There are no hidden dependencies between numerical options: the DSL
code is easier to read (than data files) and is less likely to
become obsolete (this is not that true in early developments since
the language itself and some concepts are still likely to change).
- The simulation scenario is *defined* by the script, it is the
responsibility of the user to check its meaning (not to the charge
of the code).
***** ~pugs~ language purpose
~pugs~ language is used to assemble ~C++~ tools which are well
tested. These tools should ideally be small pieces of ~C++~ code that
do one single thing and do it well.
Another purpose of the language is to allow to perform high-level
calculations. In other words, the language defines a data flow and
checks that each ~C++~ piece of code is used correctly. Since each piece
of code acts as a pure function (arguments are unchanged by calls),
the calling context is quite easy to check.
Finally it aims at simplifying the definition of new methods since
common utilities are available directly in scripts.
***** The framework: divide and conquer
~pugs~ is a research oriented software, thus generally the user is also
a developer. If this paragraph is indeed more dedicated to the
developer, by reading it, the user will have a better understanding of
the development choices and the underlying policy that the code
follows.
Actually the development framework imposed by the DSL is a guideline
for writing of new methods.
- In the process of writing a *new numerical method*, one must create
*new functions in the language*. Moreover, if a new method is close to
an existing one, it is generally *better* to use completely new
underlying ~C++~ code than to patch existing methods. Starting from a
*copy* of the existing ~C++~ code is *encouraged* for developments. This
may sound weird since classical development guidelines encourage
inheritance or early redesign. Actually, this policy is the result
of the following discussion.
- Using this approach, one *separates* clearly the *design* of a
numerical method and the *design* of the code.
- From the computer science point of view, early design for new
numerical methods is generally wrong: usually one cannot
anticipate precisely enough possible problems or method
corrections.
- It is much more difficult to introduce bugs in existing methods,
since previously validated methods are unchanged!
- For the same reason, existing method performances are naturally
unchanged by new developments.
- Also, when comparing methods, it is better to compare to the
original existing code.
- If the new method is abandoned or simply not kept, the existing
code is not polluted.
- Finally when a method is validated and ready to integrate the
mainline sources of the code, it is easier to see differences with
existing ones and *at this time* one can redesign the ~C++~ code
checking that results are unchanged and performances not
deteriorated. At this time, it is likely that the numerical method
design is finished, thus (re)designing the source code makes more
sense.
- Another consequence is that utilities are not developed again and
again.
- This implies an /a priori/ safer code: utilities are well tested and
validated.
- It saves development time obviously.
- The code of numerical methods is not polluted by environment
instructions (data initialization, error calculation,
post-processing,...)
- The counterpart is somehow classical. On the one hand, the
knowledge of existing utilities is required, this document tries
to address a part of it. On the other hand, if the developer
requires a new utility, a good practice is to discuss with the
other developers to check if it could benefit to them. Then one
can determine if it should integrate rapidly or not the main
development branch.
***** Why not python or any other scripting language?
As it was already pointed out above, general purpose languages offer
too much freedom: it is not easy to protect data. For instance in the
~pugs~ DSL, non basic variables are constant (see paragraph
[[high-level-types]]). It is important since it prevents the user from
modifying data in inconsistent ways. Also, one must keep in mind that
constraining the expressiveness is actually a strength. As said
before, one can warranty consistency of the data, perform calculations
without paying attention to the parallelism aspects,... Observe that
it is not a limitation: if the DSL's field of application needs to be
extended, it is always possible. But these extensions should never
break the rule that a DSL must not become a general purpose language.
Providing a language close to the application reduces the gap between
the application and its expression (code): Domain Specific Languages
are made for that.
Finally, python is ugly.
*** A high-level language
Following the previous discussion, the reader should now understand
the motivations that drove the design choices that conducted to build
~pugs~ as a ~C++~ toolbox driven by a user friendly language.
#+begin_verse
Keep it simple and relevant!
#+end_verse
In the development process of an application, the easy step is
generally the implementation itself. It is even more true when the
developer feels that changes to the code are natural and that the
modifications themselves look easy. Obviously, any experienced
programmer knows that writing the code is only the first step of a
much longer process which requires rigorous tests and validations (the
enrichment of a non-regression database). Generally one also desires
to ensure that the development is used within the correct bounds which
requires to implement data verification. In a perfect world, an
up-to-date documentation of the functionality and its domain of
validity should be provided.
This is even more true when defining a language (or a DSL). Enriching
a language syntax (or grammar) is not something that must be done to
answer a specific need. It must not be done /because it is possible to
do it/!
#+begin_verse
When designing a language, the difficulty is not to offer new functionalities,
it is generally to decide not to offer them.
--- Bjarne Stroustrup, C++ conference 2021.
#+end_verse
If the grammar of ~pugs~ is still likely to be extended, it should *never*
integrate low-level instructions. Low-level instructions give too much
freedom and thus are a source of errors. Several things are already
done to forbid this kind of evolution. The immutability of high-level
data is a good illustration. For instance, meshes or discrete
functions *cannot* be modified. This is not only a security to protect
the user from doing "dangerous" manipulations, but it also permits to
define high-level optimizations.
- Since meshes are constant objects, one can for instance compute
geometric data on demand. These data are kept into memory as long as
the mesh lives.
- Forbidding the modification of values of a discrete function ensures
that parallel communication instructions should never appear in a
~pugs~ script.
Another benefit of not providing low-level instructions is that the
scripts are easier to write and read, and it is more difficult to
write errors.
* Language
** Variables
In order to simplify the presentation, before going further, we
introduce the syntax that allows to print data to the terminal. It
follows ~C++~ streams construction for convenience. For instance
#+NAME: cout-preamble-example
#+BEGIN_SRC pugs :exports both :results output
cout << "2+3 = " << 2+3 << "\n";
#+END_SRC
produces the following output
#+results: cout-preamble-example
The code is quite obvious for ~C++~ users, note that ~"\n"~ is the
linefeed string (there is no character type in ~pugs~, just strings).
Actually, ~cout~ is itself a variable, we will come to this later.
~pugs~ is a strongly typed language. It means that the type of a
variable *cannot* change in its lifetime.
*** Declaration and affectation syntax
**** Declaration of simple variables
To declare a variable ~v~ of a given type ~V~, one simply writes
#+BEGIN_SRC pugs :exports source
let v:V;
#+END_SRC
This instruction reads as
#+begin_verse
Let $v\in V$.
#+end_verse
Actually,
- ~let~ is the declaration keyword,
- ~v~ is the variable name,
- ~:~ is a separation token (it can be read as "/in/" in this context),
- ~V~ is the identifier of the type, and
- ~;~ marks the end of the instruction.
For instance to declare a real variable ~x~, one writes
#+NAME: declare-R
#+BEGIN_SRC pugs :exports both :results none
let x:R;
#+END_SRC
The instructions that follow the declaration of a variable can use it
while defined in the same scope, but it cannot be used before. Also,
after its declaration, one cannot declare another variable with the
same name in the same scope (see [[blocks-and-life-time]]).
#+NAME: redeclare-variable
#+BEGIN_SRC pugs-error :exports both :results output
let x:R; // first declaration
let x:R; // second declaration
#+END_SRC
produces the following error
#+results: redeclare-variable
**** Affectation of simple variables
To affect the value of an expression ~expression~ to variable ~v~ one
simply uses the ~=~ operator. Thus assuming that a variable ~v~ has
already been /declared/, one writes simply
#+BEGIN_SRC pugs :exports source
v = expression;
#+END_SRC
There is not too much to comment, reading is quite natural
- ~v~ is the variable name,
- ~=~ is the affectation operator,
- ~expression~ is some code that provides a value of the same type as ~v~
(it can be another variable, an arithmetic expression, the result of
a function,...), and
- ~;~ marks the end of the instruction.
For instance, one can write
#+NAME: simple-affectations-example
#+BEGIN_SRC pugs :exports both :results output
import math; // to load the sin function
let a:N;
let b:N;
let x:R;
a=3;
b=2+a;
x=sin(b)+a;
cout << "a = " << a << " b = " << b << " x = " << x << "\n";
#+END_SRC
In this example, we import the ~math~ module which provides the ~sin~
function and we use the ~N~ data type of natural integers
($\mathbb{N}\equiv\mathbb{Z}_{\ge0}$).
Running the example gives the following result.
#+results: simple-affectations-example
#+BEGIN_warning
Actually, *separating* the *declaration* from the *initialization* of a
variable is quite *dangerous*. This is prone to errors and can lead to
the use of undefined values. During the compilation of scripts, ~pugs~
detects uninitialized values. For instance,
#+NAME: uninitialized-variable
#+BEGIN_SRC pugs-error :exports both :results output
let x:R;
let y:R;
y=x+1;
#+END_SRC
produces the following compilation error
#+results: uninitialized-variable
For more complex constructions, it can be very difficult to detect it
at /compile time/, this is why it is *encouraged* to use variable
definition (see [[definition-simple-variables]]).
Observe nonetheless that ~pugs~ checks at /run time/ that used variables
are correctly defined. If not, a /runtime/ error is produced.
#+END_warning
**** Definition of simple variables<<definition-simple-variables>>
The best way to define variables is the following. To define a
variable ~v~ of a given type ~V~, from an expression one writes
#+BEGIN_SRC pugs :exports source
let v:V, v = expression;
#+END_SRC
- ~let~ is the declaration keyword,
- ~v~ is the variable name,
- ~:~ is a separation token (can be read as "/in/" in this context),
- ~V~ is the identifier of the type,
- ~,~ is the separator that indicates the beginning of the affectation,
- ~expression~ is some code that provides a value of type ~V~ (it can be
another variable, an expression, the result of a function,...), and
- ~;~ marks the end of the instruction (here the definition of ~v~).
A practical example is
#+NAME: simple-definition-example
#+BEGIN_SRC pugs :exports both :results output
import math; // to load the sin function
let a:N, a=3;
let b:N, b=2+a;
let x:R, x=sin(b)+a;
cout << "a = " << a << " b = " << b << " x = " << x << "\n";
#+END_SRC
which produces the result
#+results: simple-definition-example
*** Blocks and lifetime of variables<<blocks-and-life-time>>
In pugs scripts, variables have a precise lifetime. They are defined
within scopes. The main scope is implicitly defined and sub-scopes are
enclosed between curly bracket pairs: ~{~ and ~}~. Following ~C++~, a
variable exists (can be used) as soon as it has been declared and
until the end of the scope (where it has been declared). At this point
we give a few examples.
**** A variable cannot be used before its declaration
#+NAME: undeclare-variable
#+BEGIN_SRC pugs-error :exports both :results output
n = 3;
let n:N;
#+END_SRC
#+results: undeclare-variable
**** A variable cannot be used after its definition scope
#+NAME: out-of-scope-variable-use
#+BEGIN_SRC pugs-error :exports both :results output
{
let n:N, n = 3;
n = n+2;
}
cout << n << "\n";
#+END_SRC
#+results: out-of-scope-variable-use
**** A variable name *cannot* be reused in an enclosed scope
#+NAME: nested-scope-variable-example
#+BEGIN_SRC pugs-error :exports both :results output
{
let n:N, n = 1;
}
{
let n:N, n = 2;
}
let n:N, n = 3;
{
let n:N, n = 4;
}
#+END_SRC
#+results: nested-scope-variable-example
This example is self explanatory. The same variable name can be used
in unrelated blocks (this is useful in loops for instance). However,
it cannot be reused if a variable has already been declared with the
same one in an embedding scope. This is a difference with ~C++~, the
reason is that this kind of construction is dangerous and difficult to
read.
*** Basic types<<basic-types>>
Basic types in ~pugs~ are boolean (~B~), natural integers (~N~), integers
(~Z~), real numbers (~R~), small vectors (~R^1~, ~R^2~ and ~R^3~), small
matrices (~R^1x1~, ~R^2x2~ and ~R^3x3~) and strings (~string~).
#+BEGIN_note
Observe that while mathematically, obviously $\mathbb{R} = \mathbb{R}^1
= \mathbb{R}^{1\times1}$, the data types ~R~, ~R^1~ and ~R^1x1~ are different
in ~pugs~ and are *not implicitly* convertible from one to the other!
This may sound strange but there are few reasons for that.
- First, these are the reflect of internal ~pugs~ ~C++~-data types that
are used to write algorithms. In its core design pugs aims at
writing numerical methods generically with regard to the
dimension. One of the ingredients to achieve this purpose is to use
dimension $1$ vectors and matrices when some algorithms reduce to
dimension $1$ instead of ~double~ values. To avoid ambiguity that may
arise in some situations (this can lead to very tricky code), we
decided to forbid automatic conversions of these types with
~double~. When designing the language, we adopted the same rule to
avoid ambiguity.
- A second reason is connected to the first one. Since ~pugs~ aims at
providing numerical methods for problems in dimension $1$, $2$ or
$3$, this allows to distinguish the nature of the underlying objects.
- It is natural to consider that the coordinates of the vertices
defining a mesh in dimension $d$ are elements of $\mathbb{R}^d$,
- or that a velocity or a displacement are also defined as
$\mathbb{R}^d$ values.
Thus using ~R^1~ in dimension $1$ for this kind of data makes precise
their nature in some sense .
#+END_note
**** Expression types
The language associates statically some types to several special
expressions.
- Special boolean (type ~B~) expressions. Two *keywords* allow to define
boolean values: ~true~ and ~false~.
- Integers (type ~Z~) are defined as a contiguous list of digits.
- For instance, the code ~0123~ is interpreted as the integer $123$.
- However the sequence ~123 456~ is *not interpreted* as $123456$ but as
the two integers $123$ and $456$.
- Real values (type ~R~) use the same syntax as in ~C++~. For instance,
the following expressions are accepted to define the number $1.23$.
#+BEGIN_SRC pugs :exports both
1.23;
0.123E1;
.123e+1;
123e-2;
12.3E-1;
#+END_SRC
- Small vectors values (types ~R^1~, ~R^2~ and ~R^3~) are written through
square brackets. Each component of vector values must be a scalar
expression.
#+BEGIN_SRC pugs :exports both
[1]; // R^1 value
[1,2]; // R^2 value
[1,2,3]; // R^3 value
#+END_SRC
- Small matrices values (types ~R^1x1~, ~R^2x2~ and ~R^3x3~) are written by
lines through square brackets. Each line is enclosed between square
brackets.
#+BEGIN_SRC pugs :exports both
[[1]]; // R^1x1 value
[[1,2],[3,4]]; // R^2x2 value
[[1,2,3],[4,5,6],[7,8,9]]; // R^3x3 value
#+END_SRC
- ~string~ values are defined as the set of characters enclosed between
two double quotes ( ~"~ ). The string /Hello world!/ would be simply
written as ~"Hello world!"~. Strings support the following escape
sequences (similarly to ~C++~):
#+LATEX: \definecolor{contiYellow}{RGB}{255,165,0}
#+LATEX: \rowcolors[]{2}{contiYellow!5}{contiYellow!20}
| escape | meaning |
|--------+-----------------|
| ~\'~ | single quote |
| ~\"~ | double quote |
| ~\?~ | question mark |
| ~\\~ | backslash |
| ~\a~ | audible bell |
| ~\b~ | backspace |
| ~\f~ | new page |
| ~\n~ | new line |
| ~\r~ | carriage return |
| ~\t~ | horizontal tab |
| ~\v~ | vertical tab |
These special characters are not interpreted by ~pugs~ itself but
interpreted by the system when an output is created. They are just
allowed in the definition of a ~string~.
**** Variables of basic types are stored by value
This means that each variable of basic type allocates its own memory
to store its data. This is the natural behavior for variables. It is
illustrated in the following example
#+NAME: basic-type-value-storage
#+BEGIN_SRC pugs :exports both :results output
let a:N, a = 3;
let b:N, b = a;
a = 1;
cout << "a = " << a << " b = " << b << "\n";
#+END_SRC
which produces the expected result
#+results: basic-type-value-storage
When defining ~b~, the *value* contained in ~a~ is copied to set the value
of ~b~. Thus changing ~a~'s value does not impact the variable ~b~.
**** Variables of basic types are mutable
In ~pugs~ the only variables that are mutable (their value can be
*modified*) are of basic types. Executing the following code
#+NAME: basic-type-mutable-value
#+BEGIN_SRC pugs :exports both :results output
let a:N, a = 2;
a += 3;
cout << "a = " << a << "\n";
#+END_SRC
gives
#+results: basic-type-mutable-value
which is not a surprise. However, the use of the ~+=~ operator results
in the modification of the stored value. There is no copy.
Actually, this is not really important from the user point of
view. One just has to keep in mind that, as it will be depicted below,
high-level variables *are not mutable*: their values can be *replaced* by
new ones but *cannot be modified*.
*** Implicit type conversions<<implicit-conversion>>
In order to avoid ambiguities, in ~pugs~, there is *no implicit*
conversion in general.
#+BEGIN_note
Actually, there are only three situations for which implicit type
conversion can occur: when the value is
- given as a parameter of a function,
- used as the returned value of a function, or
- used to define a tuple
This will be illustrated in section [[functions]] and [[tuples]]
#+END_note
This means that all affectations, unary operators and binary operators
are defined explicitly for supported types.
Here is a table of implicit type conversions *when allowed*.
| expected type | convertible types |
|---------------+------------------------------------------------|
| ~N~ | ~B~, ~Z~ |
| ~Z~ | ~B~, ~N~ |
| ~R~ | ~B~, ~N~, ~Z~ |
| ~R^1~ | ~B~, ~N~, ~Z~, ~R~ |
| ~R^2~ | ~0~ (special value) |
| ~R^3~ | ~0~ (special value) |
| ~R^1x1~ | ~B~, ~N~, ~Z~, ~R~ |
| ~R^2x2~ | ~0~ (special value) |
| ~R^3x3~ | ~0~ (special value) |
| ~string~ | ~B~, ~N~, ~Z~, ~R~, ~R^1~, ~R^2~, ~R^3~, ~R^1x1~, ~R^2x2~, ~R^3x3~ |
*** Operators
**** Affectation operators
In the ~pugs~ language, the affectation operators are the following.
| operator | description |
|----------+----------------------|
| ~=~ | affectation operator |
|----------+----------------------|
| ~+=~ | increment operator |
| ~-=~ | decrement operator |
| ~*=~ | multiply operator |
| ~/=~ | divide operator |
#+BEGIN_note
It is important to note that in ~pugs~ language, affectation operators
have *no* return value. This is a notable difference with ~C~ or
~C++~. Again, this is done to avoid common mistakes that can be
difficult to address. For instance, the following ~C++~ code is valid
but does not produce the expected result.
#+BEGIN_SRC C++ :exports source
bool b = true;
// do things ...
if (b=false) {
// do other things ...
}
#+END_SRC
Obviously the mistake is that the test should have been ~(b==false)~,
since otherwise, the conditional block is never executed (in ~C++~, the
result of an affectation is the value affected to the variable, which
is always false in that case.
This cannot happen with ~pugs~. A similar example
#+NAME: no-affectation-result
#+BEGIN_SRC pugs-error :exports both :results output
let b:B, b=true;
// do things
if (b=false) {
// do other things
}
#+END_SRC
produces the following error
#+results: no-affectation-result
Actually, affectations are /expressions/ in ~C++~. In ~pugs~, affectations
are /instructions/.
#+END_note
***** List of defined operator ~=~ for basic types.
As already mentioned, operator ~=~ is defined for *all* types in ~pugs~ if
the right hand side expression has the *same* type as the left hand side
variable. This is true for basic types as well as for high-level
types.
We now give the complete list of supported ~=~ affectations. The lists
are sorted by type of left hand side variable.
- ~B~: boolean left hand side variable. One is only allowed to affect boolean
values.
| ~B =~ allowed expression type |
|-----------------------------|
| ~B~ |
- ~N~: natural integer ($\mathbb{N}$ or $\mathbb{Z}_{\ge0}$) left hand side variable.
| ~N =~ allowed expression type |
|-----------------------------|
| ~B~ |
| ~N~ |
| ~Z~ (for convenience) |
- ~Z~: integer ($\mathbb{Z}$) left hand side variable.
| ~Z =~ allowed expression type |
|-----------------------------|
| ~B~ |
| ~N~ |
| ~Z~ |
- ~R~: real ($\mathbb{R}$) left hand side variable.
| ~R =~ allowed expression type |
|-----------------------------|
| ~B~ |
| ~N~ |
| ~Z~ |
| ~R~ |
- ~R^1~: vector of dimension 1 ($\mathbb{R}^1$) left hand side variable.
| ~R^1 =~ allowed expression type |
|-------------------------------|
| ~B~ |
| ~N~ |
| ~Z~ |
| ~R~ |
| ~R^1~ |
- ~R^2~: vector of dimension 2 ($\mathbb{R}^2$) left hand side variable.
| ~R^2 =~ allowed expression type |
|----------------------------------------------|
| ~R^2~ |
| ~0~ (special value) |
| list of 2 scalars (~B~, ~N~, ~Z~ or ~R~) expressions |
An example of initialization using an $\mathbb{R}^2$ value or the special value ~0~ is
#+NAME: affectation-to-R2-from-list
#+BEGIN_SRC pugs :exports both :results output
let u:R^2, u = [-3, 2.5];
let z:R^2, z = 0;
cout << "u = " << u << "\n";
cout << "z = " << z << "\n";
#+END_SRC
which produces
#+RESULTS: affectation-to-R2-from-list
Observe that ~0~ is a special value and *treated as a keyword*. There is no conversion from
integer values. For instance:
#+NAME: R2-invalid-integer-affectation
#+BEGIN_SRC pugs-error :exports both :results output
let z:R^2, z = 1-1;
#+END_SRC
produces the compilation error
#+results: R2-invalid-integer-affectation
- ~R^3~: vector of dimension 3 ($\mathbb{R}^3$) left hand side variable.
| ~R^3 =~ allowed expression type |
|---------------------------------------------|
| ~R^3~ |
| ~0~ (special value) |
| list of 3 scalar (~B~, ~N~, ~Z~ or ~R~) expressions |
An example of initialization is
#+NAME: affectation-to-R3-from-list
#+BEGIN_SRC pugs :exports both :results output
let u:R^3, u = [-3, 2.5, 1E-2];
let z:R^3, z = 0;
cout << "u = " << u << "\n";
cout << "z = " << z << "\n";
#+END_SRC
the output is
#+RESULTS: affectation-to-R3-from-list
- ~R^1x1~: matrix of dimensions $1\times1$ ($\mathbb{R}^{1\times1}$) left hand side variable.
| ~R^1x1 =~ allowed expression type |
|---------------------------------|
| ~B~ |
| ~N~ |
| ~Z~ |
| ~R~ |
| ~R^1x1~ |
- ~R^2x2~: matrix of dimension $2\times2$ ($\mathbb{R}^{2\times2}$) left hand side variable.
| ~R^2x2 =~ allowed expression type |
|---------------------------------------------|
| ~R^2x2~ |
| ~0~ (special value) |
| list of 4 scalar (~B~, ~N~, ~Z~ or ~R~) expressions |
An example of initialization using an $\mathbb{R}^{2\times2}$ value or
the special value ~0~ is
#+NAME: affectation-to-R2x2-from-list
#+BEGIN_SRC pugs :exports both :results output
let u:R^2x2, u = [[-3, 2.5],
[ 4, 1.2]];
let z:R^2x2, z = 0;
cout << "u = " << u << "\n";
cout << "z = " << z << "\n";
#+END_SRC
which produces
#+RESULTS: affectation-to-R2x2-from-list
- ~R^3x3~: matrix of dimension $3\times3$ ($\mathbb{R}^{3\times3}$) left hand side variable.
| ~R^3x3 =~ allowed expression type |
|---------------------------------------------|
| ~R^3x3~ |
| ~0~ (special value) |
| list of 9 scalar (~B~, ~N~, ~Z~ or ~R~) expressions |
An example of initialization is
#+NAME: affectation-to-R3x3-from-list
#+BEGIN_SRC pugs :exports both :results output
let u:R^3x3, u = [[ -3, 2.5, 1E-2],
[ 2, 1.7, -2],
[1.2, 4, 2.3]];
let z:R^3x3, z = 0;
cout << "u = " << u << "\n";
cout << "z = " << z << "\n";
#+END_SRC
the output is
#+RESULTS: affectation-to-R3x3-from-list
- ~string~ left hand side variable. Expressions of any basic type can be
used as the right hand side.
| ~string =~ allowed expression type |
|----------------------------------|
| ~B~ |
| ~N~ |
| ~Z~ |
| ~R~ |
| ~R^1~ |
| ~R^2~ |
| ~R^3~ |
| ~R^1x1~ |
| ~R^2x2~ |
| ~R^3x3~ |
| ~string~ |
The stored value is the same as the output value described
above. Here is an example.
#+NAME: affectation-to-string-example
#+BEGIN_SRC pugs :exports both :results output
let s_from_B:string,
s_from_B = 2>1;
let s_from_R2:string,
s_from_R2 = [ -3.5, 1.3];
let s_from_R3x3:string,
s_from_R3x3 = [[ -3, 2.5, 1E-2],
[ 2, 1.7, -2],
[1.2, 4, 2.3]];
cout << "s_from_B = " << s_from_B << "\n";
cout << "s_from_R2 = " << s_from_R2 << "\n";
cout << "s_from_R3x3 = " << s_from_R3x3 << "\n";
#+END_SRC
the output is
#+RESULTS: affectation-to-string-example
***** List of defined operator ~+=~ for basic types.
- ~B~: the ~+=~ operator is not defined for left hand side boolean
variables.
- ~N~: natural integer ($\mathbb{N}$ or $\mathbb{Z}_{\ge0}$) left hand side variable.
| ~N +=~ allowed expression type |
|------------------------------|
| ~B~ |
| ~N~ |
| ~Z~ (for convenience) |
- ~Z~: integer ($\mathbb{Z}$) left hand side variable.
| ~Z +=~ allowed expression type |
|------------------------------|
| ~B~ |
| ~N~ |
| ~Z~ |
- ~R~: real ($\mathbb{R}$) left hand side variable.
| ~R +=~ allowed expression type |
|------------------------------|
| ~B~ |
| ~N~ |
| ~Z~ |
| ~R~ |
- ~R^1~: vector of dimension 1 ($\mathbb{R}^1$) left hand side variable.
| ~R^1 +=~ allowed expression type |
|--------------------------------|
| ~R^1~ |
- ~R^2~: vector of dimension 2 ($\mathbb{R}^2$) left hand side variable.
| ~R^2 +=~ allowed expression type |
|--------------------------------|
| ~R^2~ |
- ~R^3~: vector of dimension 3 ($\mathbb{R}^3$) left hand side variable.
| ~R^3 +=~ allowed expression type |
|--------------------------------|
| ~R^3~ |
- ~R^1x1~: matrix of dimensions $1\times1$ ($\mathbb{R}^{1\times1}$) left hand side variable.
| ~R^1x1 +=~ allowed expression type |
|----------------------------------|
| ~R^1x1~ |
- ~R^2x2~: matrix of dimension $2\times2$ ($\mathbb{R}^{2\times2}$) left hand side variable.
| ~R^2x2 +=~ allowed expression type |
|----------------------------------|
| ~R^2x2~ |
- ~R^3x3~: matrix of dimension $3\times3$ ($\mathbb{R}^{3\times3}$) left hand side variable.
| ~R^3x3 +=~ allowed expression type |
|----------------------------------|
| ~R^3x3~ |
- ~string~ left hand side variable. Expressions of any basic type can be
used as the right hand side.
| ~string +=~ allowed expression type |
|-----------------------------------|
| ~B~ |
| ~N~ |
| ~Z~ |
| ~R~ |
| ~R^1~ |
| ~R^2~ |
| ~R^3~ |
| ~R^1x1~ |
| ~R^2x2~ |
| ~R^3x3~ |
| ~string~ |
The concatenated value is the same as the output value described
above, for instance
#+NAME: concatenate-to-string-example
#+BEGIN_SRC pugs :exports both :results output
let s:string, s = "foo";
s += [ -3.5, 1.3, -2];
s += "_";
s += 1>2;
s += "_";
s += [[ -3, 2.5],
[ 2, 1.7]];
cout << "s = " << s << "\n";
#+END_SRC
the output is
#+RESULTS: concatenate-to-string-example
***** List of defined operator ~-=~ for basic types.
- ~B~: the ~-=~ operator is not defined for left hand side boolean variables.
- ~N~: natural integer ($\mathbb{N}$ or $\mathbb{Z}_{\ge0}$) left hand side variable.
| ~N -=~ allowed expression type |
|------------------------------|
| ~B~ |
| ~N~ |
| ~Z~ (for convenience) |
- ~Z~: integer ($\mathbb{Z}$) left hand side variable.
| ~Z -=~ allowed expression type |
|------------------------------|
| ~B~ |
| ~N~ |
| ~Z~ |
- ~R~: real ($\mathbb{R}$) left hand side variable.
| ~R -=~ allowed expression type |
|------------------------------|
| ~B~ |
| ~N~ |
| ~Z~ |
| ~R~ |
- ~R^1~: vector of dimension 1 ($\mathbb{R}^1$) left hand side variable.
| ~R^1 -=~ allowed expression type |
|--------------------------------|
| ~R^1~ |
- ~R^2~: vector of dimension 2 ($\mathbb{R}^2$) left hand side variable.
| ~R^2 -=~ allowed expression type |
|--------------------------------|
| ~R^2~ |
- ~R^3~: vector of dimension 3 ($\mathbb{R}^3$) left hand side variable.
| ~R^3 -=~ allowed expression type |
|--------------------------------|
| ~R^3~ |
- ~R^1x1~: matrix of dimensions $1\times1$ ($\mathbb{R}^{1\times1}$) left hand side variable.
| ~R^1x1 -=~ allowed expression type |
|----------------------------------|
| ~R^1x1~ |
- ~R^2x2~: matrix of dimension $2\times2$ ($\mathbb{R}^{2\times2}$) left hand side variable.
| ~R^2x2 -=~ allowed expression type |
|----------------------------------|
| ~R^2x2~ |
- ~R^3x3~: matrix of dimension $3\times3$ ($\mathbb{R}^{3\times3}$) left hand side variable.
| ~R^3x3 -=~ allowed expression type |
|----------------------------------|
| ~R^3x3~ |
- ~string~: the ~-=~ operator is not defined for left hand side string variables.
***** List of defined operator ~*=~ for basic types.
- ~B~: the ~*=~ operator is not defined for left hand side boolean variables.
- ~N~: natural integer ($\mathbb{N}$ or $\mathbb{Z}_{\ge0}$) left hand side variable.
| ~N *=~ allowed expression type |
|------------------------------|
| ~B~ |
| ~N~ |
| ~Z~ (for convenience) |
- ~Z~: integer ($\mathbb{Z}$) left hand side variable.
| ~Z *=~ allowed expression type |
|------------------------------|
| ~B~ |
| ~N~ |
| ~Z~ |
- ~R~: real ($\mathbb{R}$) left hand side variable.
| ~R *=~ allowed expression type |
|------------------------------|
| ~B~ |
| ~N~ |
| ~Z~ |
| ~R~ |
- ~R^1~: vector of dimension 1 ($\mathbb{R}^1$) left hand side variable.
| ~R^1 *=~ allowed expression type |
|--------------------------------|
| ~B~ |
| ~N~ |
| ~Z~ |
| ~R~ |
- ~R^2~: vector of dimension 2 ($\mathbb{R}^2$) left hand side variable.
| ~R^2 *=~ allowed expression type |
|--------------------------------|
| ~B~ |
| ~N~ |
| ~Z~ |
| ~R~ |
- ~R^3~: vector of dimension 3 ($\mathbb{R}^3$) left hand side variable.
| ~R^3 *=~ allowed expression type |
|--------------------------------|
| ~B~ |
| ~N~ |
| ~Z~ |
| ~R~ |
- ~R^1x1~: matrix of dimensions $1\times1$ ($\mathbb{R}^{1\times1}$) left hand side variable.
| ~R^1x1 *=~ allowed expression type |
|----------------------------------|
| ~B~ |
| ~N~ |
| ~Z~ |
| ~R~ |
- ~R^2x2~: matrix of dimension $2\times2$ ($\mathbb{R}^{2\times2}$) left hand side variable.
| ~R^2x2 *=~ allowed expression type |
|----------------------------------|
| ~B~ |
| ~N~ |
| ~Z~ |
| ~R~ |
- ~R^3x3~: matrix of dimension $3\times3$ ($\mathbb{R}^{3\times3}$) left hand side variable.
| ~R^3x3 *=~ allowed expression type |
|----------------------------------|
| ~B~ |
| ~N~ |
| ~Z~ |
| ~R~ |
#+BEGIN_note
Observe that for these small matrix types ($\mathbb{R}^{d\times d}$) the
construction ~A *= B;~ where ~B~ is a matrix of the same type as ~A~ is not
allowed. The main reason for that is that for $d>1$ this operation has
no interest since it requires a temporary. One will see below that it
#+END_note
- ~string~: the ~*=~ operator is not defined for left hand side string variables.
***** List of defined operator ~/=~ for basic types.
- ~B~: the ~/=~ operator is not defined for left hand side boolean variables.
- ~N~: natural integer ($\mathbb{N}$ or $\mathbb{Z}_{\ge0}$) left hand side variable.
| ~N /=~ allowed expression type |
|------------------------------|
| ~N~ |
| ~Z~ (for convenience) |
- ~Z~: integer ($\mathbb{Z}$) left hand side variable.
| ~Z /=~ allowed expression type |
|------------------------------|
| ~N~ |
| ~Z~ |
- ~R~: real ($\mathbb{R}$) left hand side variable.
| ~R /=~ allowed expression type |
|------------------------------|
| ~N~ |
| ~Z~ |
| ~R~ |
- ~R^d~: the ~/=~ operator is not defined for left hand side vector (of
dimension $d\in\{1,2,3\}$) variables.
- ~R^dxd~: the ~/=~ operator is not defined for left hand side matrix (of
dimension $d\times d$ with $d\in\{1,2,3\}$) variables.
**** Unary operators
The ~pugs~ language allows the following tokens as unary operators
| operator | description |
|----------+----------------------|
| ~not~ | not operator |
| ~+~ | plus unary operator |
| ~-~ | minus unary operator |
|----------+----------------------|
| ~++~ | increment operator |
| ~--~ | decrement operator |
|----------+----------------------|
| ~[]~ | access operator |
The ~not~, ~+~ and ~-~ operators apply to the *expression* on their right. ~++~
and ~--~ operators apply only to a *variable* that can be positioned
before (pre increment/decrement) or after the token (post
increment/decrement). These operators are also inspired from their ~C++~
counterparts for convenience.
The ~+~ unary operator is a convenient operator that is *elided* when
parsing the script.
For basic types, when the operators ~not~, ~+~, ~-~, ~++~ or ~--~ are defined,
they return a value of the same type as the argument (except for the
operator ~-~ if the argument is a ~B~ or a ~N~, then the result is a
~Z~). These operators can be defined for high-level types.
- The ~not~ operator is only defined for boolean values (~B~).
- The ~-~ unary operator is defined for numeric basic types: ~B~,
~N~, ~Z~, ~R~, ~R^1~, ~R^2~, ~R^3~, ~R^1x1~, ~R^2x2~ and ~R^3x3~. It is not defined
for ~string~ variables.
- Pre and post increment operators, ~--~ and ~++~, are defined for integer
types: ~N~ and ~Z~. They are not defined for ~B~, ~R~, ~R^1~, ~R^2~, ~R^3~, ~R^1x1~,
~R^2x2~, ~R^3x3~ and ~string~ variables.
Note that the pre increment/decrement operators behave slightly
differently than their ~C++~ counterparts since they are not allowed to
be chained. In ~C++~, the following code is allowed
#+BEGIN_SRC C++ :exports source
int i = 0;
int j = ++ ++i;
#+END_SRC
In ~pugs~, it is forbidden:
#+NAME: double-pre-incr-result
#+BEGIN_SRC pugs-error :exports both :results output
let i:N, i=0;
let j:N, j = ++ ++i;
#+END_SRC
produces the compilation error
#+results: double-pre-incr-result
Again, this is done to simplify the syntax and to avoid weird
constructions.
- Access operators are only defined for small vectors ~R^d~ and small
matrices ~R^dxd~. To avoid use of uninitialized variables (or
partially uninitialized variables), these are ~read-only~ access
operators. Their syntax is the following.
#+NAME: Rd-Rdxd-access-operator
#+BEGIN_SRC pugs :exports both :results output
let x:R^2, x = [1,2];
let A:R^3x3, A = [[1,2,3],[4,5,6],[7,8,9]];
cout << "x[0] = " << x[0] << "\nx[1] = " << x[1] << "\n";
cout << "A[0,0] = " << A[0,0] << "\nA[2,1] = " << A[2,1] << "\n";
#+END_SRC
This code produces
#+results: Rd-Rdxd-access-operator
**** Binary operators
Syntax for binary operators follows again a classical structure: if
~exp1~ and ~exp2~ denote two expressions and if ~op~ denotes a binary
operator, one simply writes ~exp1 op exp2~.
Here is the list of binary operators
| keyword | operator |
|---------+--------------------------|
| ~and~ | logic and |
| ~or~ | logic or |
| ~xor~ | logic exclusive or |
|---------+--------------------------|
| ~==~ | equality |
| ~!=~ | non-equality |
| ~<~ | lower than |
| ~<=~ | lower than or equal to |
| ~>~ | greater than |
| ~>=~ | greater than or equal to |
|---------+--------------------------|
| ~<<~ | shift left |
| ~>>~ | shift right |
|---------+--------------------------|
| ~+~ | sum |
| ~-~ | difference |
| ~*~ | product |
| ~/~ | division |
Binary operators can be defined for high-level types. For basic types,
they follow a few rules.
- Logical operators ~and~, ~or~ and ~xor~ are defined for boolean operands
(type is ~B~) only. The result of the expression is a boolean.
#+begin_src latex :results drawer :exports results
\begin{equation*}
\left|
\begin{array}{rl}
\mathtt{and}:&\quad \mathbb{B} \times \mathbb{B} \to \mathbb{B}\\
\mathtt{or}:& \quad\mathbb{B} \times \mathbb{B} \to \mathbb{B}\\
\mathtt{xor}:& \quad \mathbb{B} \times \mathbb{B} \to \mathbb{B}
\end{array}
\right.
\end{equation*}
#+end_src
- Comparison operators ~==~, ~!=~, ~<~, ~<=~, ~>~ and ~>=~ are defined for all
basic scalar types and return a boolean value.
#+begin_src latex :results drawer :exports results
\begin{equation*}
\forall \mathbb{S}_1, \mathbb{S}_2 \in \{\mathbb{B},\mathbb{N},\mathbb{Z},\mathbb{R}\},
\quad
\left|
\begin{array}{rl}
\mathtt{==}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{B}\\
\mathtt{!=}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{B}\\
\mathtt{<}: & \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{B}\\
\mathtt{<=}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{B}\\
\mathtt{>}: & \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{B}\\
\mathtt{>=}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{B}
\end{array}
\right.
\end{equation*}
#+end_src
When comparing a boolean value (type ~B~) with another scalar value
type (~N~, ~Z~ or ~R~), the value ~true~ is interpreted as $1$ and the value
~false~ as $0$.
\\
For vector and matrix basic types, the only allowed operators are ~==~
and ~!=~.
#+begin_src latex :results drawer :exports results
\begin{equation*}
\forall d \in \{1,2,3\},\quad
\left|
\begin{array}{rl}
\mathtt{==}:& \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{B}\\
\mathtt{==}:& \mathbb{R}^{d \times d} \times \mathbb{R}^{d \times d} \to \mathbb{B}\\
\mathtt{!=}:& \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{B}\\
\mathtt{!=}:& \mathbb{R}^{d \times d} \times \mathbb{R}^{d \times d} \to \mathbb{B}
\end{array}
\right.
\end{equation*}
#+end_src
This is also the case for ~string~ values: only allowed operators are
~==~ and ~!=~.
#+begin_src latex :results drawer :exports results
\begin{equation*}
\left|
\begin{array}{rl}
\mathtt{==}:& \mathtt{string} \times \mathtt{string} \to \mathbb{B}\\
\mathtt{!=}:& \mathtt{string} \times \mathtt{string} \to \mathbb{B}
\end{array}
\right.
\end{equation*}
#+end_src
- Shift operators ~<<~ and ~>>~ are not used to define binary operators
between two basic types.
- Arithmetic operators (~+~, ~-~, ~*~ and ~/~) are defined for a set of
combinations of basic types. We classify them by their returned
types.
- There are no arithmetic operations that return a boolean ~B~.
- Operators that return a natural integer ~N~.
#+begin_src latex :results drawer :exports results
\begin{equation*}
\forall \mathbb{S}_1, \mathbb{S}_2 \in \{\mathbb{B},\mathbb{N}\},
\quad
\left|
\begin{array}{rl}
\mathtt{+}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{N}\\
\mathtt{*}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{N}
\end{array}
\right.
\end{equation*}
#+end_src
Boolean values (type ~B~) are not allowed as the right operand of a
division.
#+begin_src latex :results drawer :exports results
\begin{equation*}
\forall \mathbb{S} \in \{\mathbb{B},\mathbb{N}\},
\mathtt{/}: \mathbb{S} \times \mathbb{N} \to \mathbb{N}
\end{equation*}
#+end_src
Observe that ~-~ is *not* in the list.
- Operators that return an integer ~Z~.
#+begin_src latex :results drawer :exports results
\begin{equation*}
\forall \mathbb{S}_1, \mathbb{S}_2 \in \{\mathbb{B},\mathbb{N},\mathbb{Z}\},
\mbox{ such that }\mathbb{S}_1 = \mathbb{Z}\mbox{ or }\mathbb{S}_2 = \mathbb{Z},
\quad
\left|
\begin{array}{rl}
\mathtt{+}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{Z}\\
\mathtt{-}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{Z}\\
\mathtt{*}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{Z}\\
\mathtt{/}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{Z}
\end{array}
\right.
\end{equation*}
#+end_src
The divide operator does not allow boolean values as a right
operand.
#+begin_src latex :results drawer :exports results
\begin{align*}
\forall \mathbb{S}_1 \in \{\mathbb{B},\mathbb{N},\mathbb{Z}\}&\mbox{ and } \mathbb{S}_2 \in \{\mathbb{N},\mathbb{Z}\},
\mbox{ such that }\mathbb{S}_1 = \mathbb{Z}\mbox{ or }\mathbb{S}_2 = \mathbb{Z},\\
&\mathtt{/}: \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{Z}.
\end{align*}
#+end_src
Finally the following operator is also defined
#+begin_src latex :results drawer :exports results
\begin{equation*}
\mathtt{-}: \mathbb{N} \times \mathbb{N} \to \mathbb{Z}
\end{equation*}
#+end_src
The difference of two natural integers returns a /signed/ integer.
- Operators that return a real ~R~.
#+begin_src latex :results drawer :exports results
\begin{equation*}
\forall \mathbb{S}_1, \mathbb{S}_2 \in \{\mathbb{B},\mathbb{N},\mathbb{Z},\mathbb{R}\},
\mbox{ such that }\mathbb{S}_1 = \mathbb{R}\mbox{ or }\mathbb{S}_2 = \mathbb{R},
\quad
\left|
\begin{array}{rl}
\mathtt{+}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{R}\\
\mathtt{-}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{R}\\
\mathtt{*}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{R}
\end{array}
\right.
\end{equation*}
#+end_src
The divide operator does not allow boolean values as a right
operand.
#+begin_src latex :results drawer :exports results
\begin{align*}
\forall \mathbb{S}_1 \in \{\mathbb{B},\mathbb{N},\mathbb{Z},\mathbb{R}\}&\mbox{ and } \mathbb{S}_2 \in \{\mathbb{N},\mathbb{Z},\mathbb{R}\},
\mbox{ such that }\mathbb{S}_1 = \mathbb{R}\mbox{ or }\mathbb{S}_2 = \mathbb{R},\\
&\mathtt{/}: \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{Z}.
\end{align*}
#+end_src
- Operators that return a small vector ~R^d~.
#+begin_src latex :results drawer :exports results
\begin{equation*}
\forall d \in\{1,2,3\},
\quad
\left|
\begin{array}{rl}
\mathtt{+}:& \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}^d\\
\mathtt{-}:& \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}^d\\
\mathtt{*}:& \mathbb{B} \times \mathbb{R}^d \to \mathbb{R}^d\\
\mathtt{*}:& \mathbb{N} \times \mathbb{R}^d \to \mathbb{R}^d\\
\mathtt{*}:& \mathbb{Z} \times \mathbb{R}^d \to \mathbb{R}^d\\
\mathtt{*}:& \mathbb{R} \times \mathbb{R}^d \to \mathbb{R}^d\\
\mathtt{*}:& \mathbb{R}^{d\times d} \times \mathbb{R}^d \to \mathbb{R}^d
\end{array}
\right.
\end{equation*}
#+end_src
- Operators that return a small matrix ~R^dxd~.
#+begin_src latex :results drawer :exports results
\begin{equation*}
\forall d \in\{1,2,3\},
\quad
\left|
\begin{array}{rl}
\mathtt{+}:& \mathbb{R}^{d\times d} \times \mathbb{R}^{d\times d} \to \mathbb{R}^{d\times d}\\
\mathtt{-}:& \mathbb{R}^{d\times d} \times \mathbb{R}^{d\times d} \to \mathbb{R}^{d\times d}\\
\mathtt{*}:& \mathbb{B} \times \mathbb{R}^{d\times d} \to \mathbb{R}^{d\times d}\\
\mathtt{*}:& \mathbb{N} \times \mathbb{R}^{d\times d} \to \mathbb{R}^{d\times d}\\
\mathtt{*}:& \mathbb{Z} \times \mathbb{R}^{d\times d} \to \mathbb{R}^{d\times d}\\
\mathtt{*}:& \mathbb{R} \times \mathbb{R}^{d\times d} \to \mathbb{R}^{d\times d}\\
\mathtt{*}:& \mathbb{R}^{d\times d} \times \mathbb{R}^{d\times d} \to \mathbb{R}^{d\times d}
\end{array}
\right.
\end{equation*}
#+end_src
- Operators that return a ~string~.
#+begin_src latex :results drawer :exports results
\begin{equation*}
\forall \mathbb{S} \in \{\mathbb{B},\mathbb{N},\mathbb{Z},\mathbb{R}, \mathbb{R}^1, \mathbb{R}^2, \mathbb{R}^3, \mathbb{R}^{1\times1}, \mathbb{R}^{2\times2}, \mathbb{R}^{3\times3}\},
\quad
\left|
\begin{array}{rl}
\mathtt{+}:& \mathbb{S} \times \mathtt{string} \to \mathtt{string}\\
\mathtt{+}:& \mathtt{string} \times \mathbb{S} \to \mathtt{string}\\
\mathtt{+}:& \mathtt{string} \times \mathtt{string} \to \mathtt{string}
\end{array}
\right.
\end{equation*}
#+end_src
**** Operator precedence and associativity
To avoid confusions, the operators precedence in ~pugs~ language follows
the same rules as in ~C++~.
This is summarized in the following table, where ~a~ and ~b~ denote two
expressions.
| Precedence | Operator |
|------------+----------|
| 1 | ~a++~ |
| | ~a--~ |
| | ~a[]~ |
|------------+----------|
| 2 | ~++a~ |
| | ~--a~ |
| | ~not a~ |
|------------+----------|
| 3 | ~a*b~ |
| | ~a/b~ |
|------------+----------|
| 4 | ~a+b~ |
| | ~a-b~ |
|------------+----------|
| 5 | ~a<<b~ |
| | ~a>>b~ |
|------------+----------|
| 6 | ~a<b~ |
| | ~a<=b~ |
| | ~a>b~ |
| | ~a>=b~ |
|------------+----------|
| 7 | ~a==b~ |
| | ~a!=b~ |
|------------+----------|
| 8 | ~a=b~ |
| | ~a+=b~ |
| | ~a-=b~ |
| | ~a*=b~ |
| | ~a/=b~ |
As already said, we forbid some constructions to try to avoid spurious
constructions. By construction associativity of many operators makes
no sense in ~pugs~ language (affectations and increment/decrement
operators for instance).
For all other operators, the associativity rule is *left to right*. Thus
the following code
#+BEGIN_SRC pugs :exports source
- 1 + 3 - 4 + 2;
#+END_SRC
is equivalent to
#+BEGIN_SRC pugs :exports source
(((- 1) + 3) - 4) + 2;
#+END_SRC
Obviously ~pugs~ allows the use of parenthesis to write
expressions. It is enough to give a simple example.
#+NAME: parenthesis-precedence
#+BEGIN_SRC pugs :exports both :results output
cout << " 2 + 3 * 4 = " << 2 + 3 * 4 << "\n";
cout << "(2 + 3) * 4 = " << (2 + 3) * 4 << "\n";
#+END_SRC
the output is
#+RESULTS: parenthesis-precedence
*** High-level types<<high-level-types>>
Aside from the basic types described in the previous section, ~pugs~
also deals with "high-level" types. This term is better understood as
"non-basic types". The ~pugs~ language is not object oriented to keep it
simple.
To fix ideas, let us give examples of a few high-level data
types. These can be meshes (the ~mesh~ type), output streams (the
~ostream~ type), boundary descriptors (~boundary~), quadrature formula
(~quadrature~), discrete functions (~Vh~),...
One can already see that the complexity of these types may vary a lot.
The main difference between these types and basic types is that,
high-level types are not available directly in the language but they
are loaded on demand (using the ~import~ keyword in the preamble of the
script).
The second difference is that data of these types are *constant*. More
precisely, the content of high-level variables can be replaced by a
new one *but* it cannot be modified. For this reason, the following
operators can never be applied to variables of these kinds
| forbidden operators | description |
|---------------------+-----------------------------|
| ~++~ | pre/post increment operator |
| ~--~ | pre/post decrement operator |
| ~+=~ | assignment by sum |
| ~-=~ | assignment by difference |
| ~*=~ | assignment by product |
| ~/=~ | assignment by quotient |
We conclude by stating that if access operator ~[]~ can possibly be
defined for high-level types, it should be done with care. It is not
recommended.
Finally, the last difference lies in the fact that high-level types
use shallow copies and not value copies as it is the case for basic
types. This is transparent to the user and provides the intuitive
(similar) behavior since data of high-level variables are constant. To
illustrate this, let us consider the following example.
#+BEGIN_SRC pugs :exports both
import mesh;
let m1:mesh, m1 = cartesianMesh(0, [1,1], (10,10));
let m2:mesh, m2 = m1;
let m3:mesh, m3 = cartesianMesh(0, [1,1], (10,10));
m3 = m1;
#+END_SRC
In this example, we are dealing with 3 ~mesh~ variables.
- First, ~m1~ is defined as a uniform cartesian mesh in dimension 2 of
$]0,1[^2$ made of $10\times10$ identical square cells. Let us call this
resident mesh $\mathcal{M}_a$ for clarity.
- Next, the variable ~m2~ is defined as a copy of the variable ~m1~. It is
a copy of the variable, but its content is *the same*. The variable ~m2~
is also referring to the mesh $\mathcal{M}_a$: there is no
duplication of the mesh in memory.
- Then, the variable ~m3~ is defined to refer to a *new* mesh (called
$\mathcal{M}_b$. It is /similar/ to $\mathcal{M}_a$, but it is *not* the
same one! When ~m3~ is defined, two meshes $\mathcal{M}_a$ and
$\mathcal{M}_b$ (and two distinct connectivities) are resident in
memory.
- Finally, the last instruction (~m3 = m1;~) sets ~m3~ to also refer to
$\mathcal{M}_a$. Since no other variable refers to $\mathcal{M}_b$,
this mesh is destroyed (memory is freed). At the end of the program,
all the variables ~m1~, ~m2~ and ~m3~ are referring to $\mathcal{M}_a$
which is the only mesh that resides in memory.
*** Compound types
The ~pugs~ language allows to deal with compound types. The idea is to
define a list of variables as a member of a product space (each
variable belongs to one of them, each variable has a simple type:
basic or high-level).
**** Compound declaration
Let us provide an example to fix ideas.
#+NAME: compound-declaration
#+BEGIN_SRC pugs :exports both :results output
let (A,x,n): R^2x2*R^3*N;
A = [[1,2],[3,4]];
x = [2,4,6];
n = 2;
cout << "A = " << A << "\nx = " << x << "\nn = " << n << "\n";
#+END_SRC
the output is
#+RESULTS: compound-declaration
This is completely equivalent to declaring the variables one after the
other.
**** Compound definition
One can also use the following definition instruction
#+NAME: compound-definition
#+BEGIN_SRC pugs :exports both :results output
let (A,x,n): R^2x2*R^3*N, (A,x,n) = ([[1,2],[3,4]], [2,4,6], 2);
cout << "A = " << A << "\nx = " << x << "\nn = " << n << "\n";
#+END_SRC
which gives the same result
#+RESULTS: compound-definition
A potential mistake with this construction, is that variables
must appear in the same order in the declaration part and in the
affectation part of a definition. For instance, the following code is
invalid.
#+NAME: invalid-compound-definition
#+BEGIN_SRC pugs-error :exports both :results output
let (x,y):R*R, (y,x) = (0,1);
#+END_SRC
It produces the following error
#+results: invalid-compound-definition
which is easy to fix.
Another potential mistake is that all variables of the list are marked
as defined at the same time. Thus one cannot use construction like
this:
#+NAME: undeclared-compound-definition
#+BEGIN_SRC pugs-error :exports both :results output
let (x,y):R*R, (x,y) = (0,2+x);
#+END_SRC
It produces the following error
#+results: undeclared-compound-definition
While the variable ~x~ is defined *before* ~y~, this kind of construction is
forbidden. From a technical point of view, this behavior would be easy
to change (allow to use the fresh value of ~x~ in the definition of ~y~),
but this makes the code unclear and this is not the purpose of
compound types.
#+BEGIN_note
Observe that there is no implicit conversion when dealing with
compound types. The ~=~ operators are used sequentially to set the
different variables.
#+END_note
**** Compound affectation
The last way to use compound types is illustrated by the following
example.
#+NAME: compound-affectation
#+BEGIN_SRC pugs :exports both :results output
let A: R^2x2;
let x: R^3;
let n: N;
(A,x,n) = ([[1,2],[3,4]], [2,4,6], 2);
cout << "A = " << A << "\nx = " << x << "\nn = " << n << "\n";
#+END_SRC
It produces again
#+results: compound-affectation
#+BEGIN_note
Observe that the only operator allowed for this kind of construction
is the operator ~=~.
#+END_note
**** Use of compound types
Actually using compound types the way it is presented in this
paragraph is not recommend. The purpose of compound types and compound
affectations is related to functions. As one will see in section
[[functions]], functions can return compound values, thus compound
affectations (or definitions) are needed to get returned values in
that case.
*** Tuple types<<tuples>>
The last special type construction is the ability to deal with tuples
in the ~pugs~ language. The tuples we are describing here are lists of
data of a *unique* and *simple* type (one cannot use compound types to
define tuples). Tuples in the ~pugs~ language have *no relationship* with
~C++~'s ~std::tuple~ these are two completely different notions.
The list of values given to the tuple must be *implicitly* convertible
to the type of tuple elements (see the conversion table in section
[[implicit-conversion]]). There is no ambiguity, the implicit conversions
follow the rules of operator ~=~.
**** Tuple declaration and affectation
For instance, one may write
#+NAME: tuple-declaration-affectation
#+BEGIN_SRC pugs :exports both :results output
let x:(R);
x = (1,2,3.4,2);
cout << x << "\n";
#+END_SRC
Executing this code, one gets
#+results: tuple-declaration-affectation
**** Tuple definition
The definition syntax is also possible.
#+NAME: tuple-definition
#+BEGIN_SRC pugs :exports both :results output
let x:(R^2), x = ([1,2],[3.4,2], 0, [2,3]);
cout << x << "\n";
#+END_SRC
This code gives
#+results: tuple-definition
Observe that the special value ~0~ is used there.
**** Tuple purpose
Tuple variables are just lists of data of the same type. In the ~pugs~
language, one cannot access to a specific value of the list nor alter
one of them. This is not something that should ever change. Tuples are
not arrays! The ~pugs~ language is not meant to allow low-level
instructions.
The use case of tuples is to provide lists of data to the ~C++~
underlying methods. A classical example is to provide a set of
boundary conditions to a method.
**** Tuple and other variables
If one cannot access to specific values of tuple variables, one can
however dispatch their values to compound variable lists or to simple
variable. To do so the tuple content must be convertible to each of
the compound variable using the operator ~=~ rules (see the conversion
table in section [[implicit-conversion]]). The one to one matching is
checked at runtime.
#+NAME: tuple-to-compound
#+BEGIN_SRC pugs :exports both :results output
let t:(R), t = (1, 2, 3.4, -5);
let (x, y, u, A):R*R*R^1*R^1x1, (x, y, u, A) = t;
cout << "t = " << t << "\n";
cout << "x = " << x << "\n";
cout << "y = " << y << "\n";
cout << "u = " << u << "\n";
cout << "A = " << A << "\n";
#+END_SRC
This code gives
#+results: tuple-to-compound
#+BEGIN_note
One could have separate the declarations of ~x~, ~y~, ~u~ and ~A~ from the
affectation. Writing simply at some point
#+BEGIN_SRC pugs :exports code
(x, y, u, A) = t;
#+END_SRC
#+END_note
If the tuple contains only one entry, one can assign its value to a
simple variable.
#+NAME: tuple-to-variable
#+BEGIN_SRC pugs :exports both :results output
let t:(R^3), t = [1, 2, 3];
let x:R^3, x = t;
cout << "t = " << t << "\n";
cout << "x = " << x << "\n";
#+END_SRC
One gets
#+results: tuple-to-variable
If the size of the tuple does not match the number of space defining
the compound type, one gets a runtime error. For instance
#+NAME: tuple-to-compound-wrong-size
#+BEGIN_SRC pugs-error :exports both :results output
let t:(Z), t = (1, 2);
let (a, b, c):R*R*R, (a, b, c) = t;
#+END_SRC
produces
#+results: tuple-to-compound-wrong-size
or
#+NAME: tuple-to-variable-wrong-size
#+BEGIN_SRC pugs-error :exports both :results output
let t:(Z), t = (1, 2);
let a:R, a = t;
#+END_SRC
which gives
#+results: tuple-to-variable-wrong-size
** Statements
The ~pugs~ language supports classical statements to control the data
flow. For simplicity, these statements syntax follow their ~C++~
counterpart. The only statement that is not implemented in ~pugs~ is the
~switch...case~. This may change but on the one hand, up to now it has
never been necessary (up to now, we did not encountered the need to
chain ~if...else~ statements), and on the other hand, keeping the
language as simple as possible remains the policy in ~pugs~ development.
*** ~if...else~ statement
The simplest statement is the conditional statement ~if...else~. If a
condition is satisfied, the ~truestatement~ (a single instruction or
block of instructions) is executed. In the other case, if the /optional/
~else~ keyword is used, the ~falsestatement~ (a single instruction or
block of instructions) is executed.
#+BEGIN_SRC pugs :exports code
if (condition) truestatement
if (condition) truestatement else falsestatement
#+END_SRC
The condition itself *must be* a boolean value (of type ~B~).
#+NAME: if-instruction
#+BEGIN_SRC pugs :exports both :results output
if (2>1) cout << "yes: 2 > 1\n";
if (not (2>1)) cout << "hmm... 2 <= 1 ?\n";
#+END_SRC
generates the following expected output. The second output is
hopefully not executed.
#+results: if-instruction
It is probably better to use the ~else~ keyword to obtain the same
result.
An example of the use of the ~else~ keyword follows.
#+NAME: if-else-instruction
#+BEGIN_SRC pugs :exports both :results output
if (2<=1)
cout << "hmm... 2 <= 1 ?\n";
else
cout << "ok: 2 > 1\n";
#+END_SRC
#+results: if-else-instruction
However, it is recommended even in the case of single instructions in
the conditional statements to use instruction blocks:
#+NAME: if-else-block-of-one
#+BEGIN_SRC pugs :exports both :results output
import math;
let sin100_positive:B;
if (sin(100) > 0){
sin100_positive = true;
} else {
sin100_positive = false;
}
cout << "sin(100)>0: " << sin100_positive << "\n";
#+END_SRC
#+results: if-else-block-of-one
We give a final illustration
#+NAME: if-block
#+BEGIN_SRC pugs :exports both :results output
let m:R, m = 2.5;
let M:R, M = 2;
let has_swapped:B, has_swapped = false;
if (m>M){
let tmp:R, tmp = m;
m = M;
M = tmp;
has_swapped = true;
}
cout << "min = " << m << " max = " << M << " has_swapped = " << has_swapped << "\n";
#+END_SRC
#+results: if-block
*** ~for~ loops
~pugs~ allows to write ~for~ loops. It follows the ~C++~ syntax
#+BEGIN_SRC pugs :exports code
for (declarationinstruction ; condition ; postinstruction) statement
#+END_SRC
Following ~C++~, ~declarationinstruction~, ~condition~ and ~postinstruction~
are optional. The ~condition~ argument, if it is present, *must* be a
boolean value (type ~B~). If it is absent, the default value ~true~ is
used.
- The ~declarationinstruction~ is executed only *once* /before/ the beginning
of the loop. The lifetime of the variable declared here is defined
by the ~for~ instruction itself.
- The ~condition~ is evaluated /before/ each loop.
- The ~postinstruction~ is executed /after/ each loop.
- The ~statement~ is either a single instruction or a block of
instructions. The ~statement~ is executed if the ~condition~ has the
value ~true~.
For instance, one can write
#+NAME: for-block
#+BEGIN_SRC pugs :exports both :results output
let n:N, n = 10;
let sum:N, sum = 0;
for (let i:N, i=1; i<=n; ++i) {
sum += i;
}
cout << "sum = " << sum << "\n";
#+END_SRC
which gives as expected
#+results: for-block
The lifetime of the declared variable (in the ~declarationinstruction~
statement) is illustrated by the following example
#+NAME: for-scope-error
#+BEGIN_SRC pugs-error :exports both :results output
for (let i:N, i=0; i<2; ++i) {
cout << "i = " << i << "\n";
}
cout << "i = " << i << "\n";
#+END_SRC
Running this example produces the following error
#+results: for-scope-error
To fix the previous code, one can write
#+NAME: for-no-decl
#+BEGIN_SRC pugs :exports both :results output
let i:N, i=0;
for (; i<2; ++i) {
cout << "i = " << i << "\n";
}
cout << "i = " << i << "\n";
#+END_SRC
One gets
#+results: for-no-decl
*** ~do...while~ loops
The second kind of loops that is supported is the ~do...while~
construction which executes /at least/ one time the enclosed ~statement~.
#+BEGIN_SRC pugs :exports code
do statement while (condition);
#+END_SRC
The ~statement~ is either a single instruction or a block of
instructions. The ~condition~ is an expression of boolean value (type
~B~).
#+NAME: do-while-block
#+BEGIN_SRC pugs :exports both :results output
let sum:N, sum = 0;
let i:N, i = 0;
do {
sum += i;
++i;
} while (i<=10);
cout << "sum = " << sum << "\n";
#+END_SRC
It gives
#+results: do-while-block
*** ~while~ loops
The last kind of loops that is allowed in ~pugs~ language is the ~while~
loop.
#+BEGIN_SRC pugs :exports code
while (condition) statement
#+END_SRC
The ~statement~ is either a single instruction or a block of
instructions. The ~condition~ is an expression of boolean value (type
~B~).
This time if the ~condition~ is not satisfied (~false~ when reaching the
~while~ instruction), the ~statment~ is never executed.
An example of the ~while~ loop is the following.
#+NAME: while-block
#+BEGIN_SRC pugs :exports both :results output
let sum:N, sum = 0;
let i:N, i = 1;
while (i<=10) {
sum += i;
++i;
}
cout << "sum = " << sum << "\n";
#+END_SRC
The result is
#+results: while-block
*** ~break~ and ~continue~ keywords.
These *keywords* are used to control loops behavior from enclosed loop
statements. They follow their ~C++~ counterparts.
- The ~continue~ keyword is used to skip the instructions that follow in
the loop statement and continue the loop.
- The ~break~ keyword leaves the current loop.
An example of use of the ~continue~ keyword is
#+NAME: nested-continue
#+BEGIN_SRC pugs :exports both :results output
for (let i:N, i = 0; i < 5; ++i) {
cout << i << ": ";
for (let j:N, j=0; j < 5; ++j) {
if (j<i) continue;
cout << j << " ";
}
cout << "\n";
}
#+END_SRC
The result is
#+results: nested-continue
Replacing the ~continue~ keyword by a ~break~ (and changing the test)
#+NAME: nested-break
#+BEGIN_SRC pugs :exports both :results output
for (let i:N, i = 0; i < 5; ++i) {
cout << i << ": ";
for (let j:N, j=0; j < 5; ++j) {
if (j>i) break;
cout << j << " ";
}
cout << "\n";
}
#+END_SRC
changes the output to
#+results: nested-break
Obviously the behavior is the same using ~do...while~ or ~while~ loops.
** Functions<<functions>>
The ~pugs~ language allows the manipulation and definition of
functions. These are *mathematical functions* and are not functions as
functions in ~C++~: functions in the ~pugs~ language are *not subroutines*.
To be more precise, a function $f$ follows the following structure
\begin{align*}
\mbox{let }f:& X_1\times \cdots \times X_n \to Y_1\times \cdots \times Y_m,\\
& (x_1,\ldots,x_n)\mapsto (y_1,\ldots,y_m)=f(x_1,\ldots,x_n),
\end{align*}
where $n,m\in\mathbb{N}_{\ge1}$, and where ${(X_i)}_{1\le i\le n}$ and ${(Y_i)}_{1\le
i\le m}$ are /simple/ types. Actually $X_1\times \cdots \times X_n$ and
$Y_1\times \cdots \times Y_m$ are /compound/ types.
Thus assuming that the function $f$ is defined in the language as ~f~,
one can use the following syntax to evaluate it. This is a
pseudo-code, real examples will follow. Assuming that ~(x1...,xn)~ has
been properly defined in ~X1*...*Xn~ one can write
#+BEGIN_SRC pugs :exports code
let (y1,...,ym):Y1*...*Ym, (y1,...,ym) = f(x1,...,xn);
#+END_SRC
or if ~(y1,...,ym)~ has already been declared in ~Y1*...*Ym~
#+BEGIN_SRC pugs :exports code
(y1,...,ym) = f(x1,...,xn);
#+END_SRC
The ~pugs~ language handles two kinds of functions. User-defined
functions (see [[user-defined-functions]]) and builtin functions (see
[[builtin-functions]]). They behave essentially the same and are both
constant. In this context, it means that one cannot just declare or
modify a function.
*** Pure functions
In the ~pugs~ language, functions are *pure functions* in the sense that
arguments given to a function are *never* modified by the function. They
act as operators.
#+BEGIN_note
Actually these functions are not strictly /pure functions/ in the
computer science sense. The reason is that they may have side
effects. As an example, it is possible to modify the random seed used
by the code. In that case, the modified value is not a variable of the
language itself but the internal random seed.
#+END_note
*** Implicit type conversion for parameters and returned values
As already said, in general, there is no implicit type conversion in
the ~pugs~ language. The only exceptions are when initializing tuples
(see [[tuples]]), passing arguments to functions and dealing with returned
values. See the conversion table in section [[implicit-conversion]].
*** User-defined functions<<user-defined-functions>>
To define user functions, the syntax mimics mathematics. The
definition of the function itself is a *single* expression. This means
that one cannot use ~if...else~, loops and there is not such keyword as
~return~ in the ~C++~.
The syntax of the definition of functions is
#+BEGIN_SRC pugs :exports code
let f1 : X -> Y, x -> e;
let f2 : X -> Y1*...*Ym, x -> (e1,...,em);
let f3 : X1*...*Xn -> Y, (x1,...,xn) -> e;
let f4 : X1*...*Xn -> Y1*...*Ym, (x1,...,xn) -> (e1,...,em);
#+END_SRC
~X~, ~Y~, ~X1~, ..., ~Xn~ and ~Y1~, ..., ~Ym~ are /simple/ data types (not
compound). The type of ~x~ is ~X~ and for lists, each argument ~xi~ belongs
to ~Xi~. The function themselves are defined by the expressions ~e~ (or
~e1~, ..., ~em~) which types are set by ~Y~ (or ~Y1~, ..., ~Ym~).
#+BEGIN_warning
- ~X~, ~X1~, ..., ~Xn~ cannot be tuples spaces.
- ~Y~, ~Y1~, ..., ~Ym~ can be tuple spaces.
Also, ~X~ can be the empty space $\varnothing$, but none of ~X1~, ..., ~Xn~
(for $n>1$) and none of ~Y~, ~Y1~, ..., ~Ym~ can be $\varnothing$. In other
word one can define functions such as
\begin{align*}
\mbox{let }g:& \varnothing \to Y_1\times \cdots \times Y_m,\\
& \varnothing\mapsto (y_1,\ldots,y_m)=g(\varnothing).
\end{align*}
The *keyword* associated to $\varnothing$ is ~void~. Keep in mind that ~void~
is just a keyword in ~pugs~, it is not a type in the computer science
point of view (one cannot declare variables of type ~void~).
#+END_warning
Let us give a few examples.
#+NAME: R-to-R-function
#+BEGIN_SRC pugs :exports both :results output
let f: R -> R, x -> -x;
cout << "f(2.3) = " << f(2.3)
<< "\nf(1) = " << f(1)
<< "\nf(true) = " << f(true)
<< "\n";
#+END_SRC
One observes the implicit conversion of the values passed by argument.
#+results: R-to-R-function
The following example illustrates the implicit conversion to the
returned type
#+NAME: R-to-R1-function
#+BEGIN_SRC pugs :exports both :results output
let f: R -> R^1, x -> 2*x;
cout << "f(3.2) = " << f(3.2) << "\n";
#+END_SRC
#+results: R-to-R1-function
But one cannot use tuples to define the domain
#+NAME: no-tuple-in-domain
#+BEGIN_SRC pugs-error :exports both :results output
let f: (N)->N, n -> 3;
#+END_SRC
produces the following compilation time error
#+results: no-tuple-in-domain
Using compound types as input and output, one can write
#+NAME: R22-R-string-to-R-string-function
#+BEGIN_SRC pugs :exports both :results output
let f : R^2x2*R*string -> R*string*R^3x3,
(A,x,s) -> (x*A[0,0]*A[1,1]-A[1,0],
A+","+x+","+s,
[[A[0,0], A[0,1], 0],
[A[1,0], A[1,1], 0],
[ 0, 0, x]]);
let x : R, x=0;
let s : string, s= "";
let A : R^3x3, A = 0;
(x,s,A) = f([[1,2],[3,4]], 2.3, "foo");
cout << "x = " << x
<< "\ns = " << s
<< "\nA = " << A << "\n";
let (y,t,A2):R*string*R^3x3, (y,t,A2) = f([[3,1],[4,2]], -5.2, "bar");
cout << "y = " << y
<< "\nt = " << t
<< "\nA2 = " << A2 << "\n";
#+END_SRC
This meaningless example produces the following result.
#+results: R22-R-string-to-R-string-function
The following example shows how to use tuples as codomain.
#+NAME: R-to-tuple-R-function
#+BEGIN_SRC pugs :exports both :results output
let f: R -> (R), x -> (2*x, 2-x, 7*x-2, 2, x/3);
cout << "f(3.2) = " << f(3.2) << "\n";
#+END_SRC
#+results: R-to-tuple-R-function
Finally, we give an example of an empty function. Observe the special
role of the ~void~ keyword
#+NAME: void-to-tuple-R-function
#+BEGIN_SRC pugs :exports both :results output
let f: void -> (R^2), void -> ([2, 3.5], [6, 7]);
cout << "f() = " << f() << "\n";
#+END_SRC
#+results: void-to-tuple-R-function
**** Lifetime of function arguments
The arguments used to define a function are *local* variables that exist
only during the evaluation of the function.
Let us consider the following example
#+NAME: lifetime-of-function-args
#+BEGIN_SRC pugs :exports both :results output
let a:R, a = 1.4;
let plus: R -> R, a -> (a>0)*a;
cout << "plus(1) = " << plus(1) << "\n";
cout << "plus(-3.2) = " << plus(-3.2) << "\n";
cout << "a = " << a << "\n";
#+END_SRC
This gives the expected result: the value of the variable ~a~ is
unchanged.
#+results: lifetime-of-function-args
**** Non-argument variables in function expressions
Here we discuss rapidly of using variables (which are not arguments)
in function expressions.
#+NAME: non-arg-variables-in-functions
#+BEGIN_SRC pugs :exports both :results output
let a:R, a = 1.4;
let plus: R -> R, x -> (x>0)*a;
cout << "a = " << a << "\n";
cout << "plus(2.2) = " << plus(2.2) << "\n";
cout << "plus(-3.2) = " << plus(-3.2) << "\n";
a = 2;
cout << "a = " << a << "\n";
cout << "plus(2.2) = " << plus(2.2) << "\n";
cout << "plus(-3.2) = " << plus(-3.2) << "\n";
#+END_SRC
Running the example, one gets
#+results: non-arg-variables-in-functions
While the function itself is a constant object, one sees that since
the value of ~a~ is changed, the function value is implicitly
modified. /This is a dangerous feature and should be avoided!/
Since functions themselves are variables one can use functions in
function expressions.
#+NAME: functions-in-functions
#+BEGIN_SRC pugs :exports both :results output
let plus: R -> R, x -> (x>0)*x;
let minus: R -> R, x -> -(x<0)*x;
let pm : R -> R*R, x -> (plus(x), minus(x));
let toR2: R*R -> R^2, (x,y) -> [x,y];
cout << "pm(2) = " << toR2(pm(2)) << " pm(-3) = " << toR2(pm(-3)) << "\n";
#+END_SRC
One observes the utility function ~toR2~ that is used to perform the
output since ~cout~ does not handle compound types output. One gets
#+results: functions-in-functions
**** Lifetime of user-defined functions
Since functions are somehow variables, the lifetime of functions
follows similar rules.
Let us give an example
#+NAME: functions-lifetime
#+BEGIN_SRC pugs :exports both :results output
{
{
let f: R -> R, x -> 3.25*x-0.3;
cout << "block 1.1: f(2) = " << f(2) << "\n";
}
let f: R -> R, x -> 1.25*x+3;
cout << "block 1: f(2) = " << f(2) << "\n";
}
for (let i:N, i = 1; i<=2; ++i) {
let f: R -> R, x -> i*x+2;
cout << "for(i=" << i << "): f(2) = " << f(2) << "\n";
}
#+END_SRC
Running this example produces
#+results: functions-lifetime
**** Recursion
Since functions cannot be simply declared but must be defined, double
recursion is not possible. Thus, there is no equivalent construction
in ~pugs~ of the following ~C++~ code
#+BEGIN_SRC C++ :exports source
int f(int i);
int g(int i)
{
if (i < 0)
return 0;
else
return f(i-1);
}
int f(int i) { return g(i); }
#+END_SRC
Moreover simple recursion is forbidden for similar reasons, since
functions are made of a single expression (one cannot use statements),
there would be no way to end the recursion. Thus, the code
#+NAME: no-recursion
#+BEGIN_SRC pugs-error :exports both :results output
let f: N->N, n -> f(2*n);
#+END_SRC
produces the following compilation time error
#+results: no-recursion
*** Builtin functions<<builtin-functions>>
In ~pugs~ language, builtin functions are ~C++~ pieces of code that can be
called in scripts. Their usage is very similar to user-defined
functions. They differ from user-defined functions in three points.
- Builtin functions may have no returned value.
- Builtin functions can use tuples as arguments.
- Builtin functions are polymorphic. More precisely, this means that
the signature of a builtin function is also defined by its expected
argument types.
- Builtin functions can take user-defined functions as parameters.
- user-defined functions cannot take functions as parameters
- builtin functions cannot take builtin functions as parameters
(actually, this is not a limitation since it is trivial to embed a
builtin function into a user-defined one).
Here is a simple example of builtin function embedded in a user
function.
#+NAME: builtin-function-embedding
#+BEGIN_SRC pugs :exports both :results output
import math;
let cosinus: R -> R, x -> cos(x);
cout << "cosinus(2) = " << cosinus(2) << "\n";
#+END_SRC
Running this example produces
#+results: builtin-function-embedding
Builtin functions are defined when modules are imported, see [[modules]].
** Modules<<modules>>
Modules are sets of *builtin functions* and *high-level types* that are
not provided by default by ~pugs~. To load a module, one has to use an
~import~ instruction in the preamble of the script.
For instance to load the ~math~ module which contains builtin
mathematical functions, one writes in the preamble of the script
#+BEGIN_SRC pugs :exports source
import math;
#+END_SRC
#+BEGIN_warning
A work in progress
- At the time of writing this documentation, one should note that
module inter-dependencies are still not implemented.
- Also, (and especially with regard to the ~scheme~ module), module
contents are likely to change and to be reorganized.
- Finally it is almost sure that modules will be equipped with a
/namespace/-like functionality to avoid conflicts. This actually
should make the scripts cleaner since, the naming of functions would
be more natural.
#+END_warning
One can access the list of available modules inside the language.
#+NAME: get-available-modules
#+BEGIN_SRC pugs :exports both :results output
cout << getAvailableModules() << "\n";
#+END_SRC
The output lists all available modules
#+RESULTS: get-available-modules
Let us comment a bit this output. One notices that there are two kinds
of modules. Modules that are automatically imported (tagged with a ~*~)
and the other ones.
In this section we will not describe exhaustively the whole module
content but will give the basic information that should allow the user
to find his way. To do so, it is important to examine carefully the
content of the ~core~ module, since it contains some helper functions.
*** The ~core~ module
As already said, the ~core~ module cannot be imported manually since it
is imported automatically.
**** ~core~ provided types
***** ~ostream~
The ~core~ module provides an important type: ~ostream~. This type has
already been encountered in this documentation. This is the type of
some specific predefined objects: ~cout~, ~cerr~ and ~clog~. These objects
are very similar to their ~C++~ counterparts, respectively ~std::cout~,
~std::cerr~ and ~std::clog~.
Objects of type ~ostream~ can be used as the left operand for the binary
operator ~<<~ for basic types and their associated tuples.
| ~ostream <<~ allowed expression type |
|------------------------------------|
| ~B~ |
| ~N~ |
| ~Z~ |
| ~R~ |
| ~R^1~ |
| ~R^2~ |
| ~R^3~ |
| ~R^1x1~ |
| ~R^2x2~ |
| ~R^3x3~ |
| ~string~ |
As one should have noticed, when the left operand of the binary
operator ~<<~ is an ~ostream~, the result of the operation is also an
~ostream~, which allows to chain output.
One can overload the ~ostream <<~ construction for high-level types.
Other variables of the type ~ostream~ can be created (in order to write
to files for instance) as we will see below.
**** ~core~ provided functions
***** execution control functions
Here are functions that allow to control the execution of the
script. These can stop the execution if some conditions are met or
create checkpoint that may be used to stop and then resume the
execution.
****** ~exit: Z -> void~
This function interrupts the execution of the script. The integer (~Z~)
value is the code that is returned when ~pugs~ exits.
#+NAME: exit-function
#+BEGIN_SRC pugs :exports both :results output
for (let i:N, i = 1; i < 10; ++i) {
cout << "i = " << i << "\n";
if (i==3) {
exit(0);
}
}
#+END_SRC
The output shows that the execution is interrupted (in a clean way) if
the condition (here ~i == 3~) is met.
#+RESULTS: exit-function
****** ~stop: void -> B~
This is an interactive function that is used to notify the code to
stop its execution. It returns ~true~ if one of the following stopping
condition is met.
- If a ~stop~ file is created *after* the beginning of the execution in
the *execution directory*. This file can be created using the ~touch~
command for instance
#+BEGIN_SRC shell :exports source
touch stop
#+END_SRC
- If ~pugs~ was compiled using ~Slurm~ support (this can be checked using
the command ~./pugs -v~) then during batch execution the ~stop~ function
will also return ~true~ as soon as the remaining execution time is
less than $150s$.
In the following example the function ~stop~ returns always ~false~ since
none of the above conditions are satisfied.
#+NAME: stop-function
#+BEGIN_SRC pugs :exports both :results output
for (let i:N, i = 1; i <= 5; ++i) {
cout << "i = " << i << "\n";
if (stop()) {
exit(0);
}
}
#+END_SRC
The output shows that the execution reaches ~i == 5~.
#+RESULTS: stop-function
****** ~checkpoint: void -> void~
This function creates a checkpoint that can be used as a starting
point for another execution. This function can be placed anywhere in
the script and possibly multiple times. The checkpoint storage file is
~checkpoint.h5~ and it is created in the current execution directory.
*This functionality requires ~pugs~ to be compiled with ~HDF5~ support*.
The following example creates two checkpoints.
#+NAME: checkpoint-function
#+BEGIN_SRC pugs :exports both :results output
let n:N, n = 4;
cout << "n = " << n << "\n";
checkpoint();
let x:R^3, x = [1, 2, 3];
cout << "x = " << x << "\n";
checkpoint();
#+END_SRC
The output displays the creation of two checkpoints in the file ~checkpoint.h5~.
#+RESULTS: checkpoint-function
The command
#+BEGIN_SRC shell :exports code
pugs_checkpoint --info checkpoint.h5
#+END_SRC
displays a simple execution state for each checkpoint stored in the
file in order to help user to choose the appropriate resuming
checkpoint.
#+NAME: checkpoint-info
#+BEGIN_SRC shell :exports results :results output
${PUGS_CHECKPOINT} --no-color --info checkpoint.h5
#+END_SRC
#+RESULTS: checkpoint-info
One can notice that the ~x~ variable does not appear in the first
checkpoint since it was not created at this point. Also two pointers
are defined:
- ~last_checkpoint~ that displays the more recent checkpoint,
and
- ~resuming_checkpoint~ that indicates which checkpoint will be used when
resuming.
Proceeding with this example, one can change the resuming checkpoint
with the command
#+BEGIN_SRC shell :exports code
pugs_checkpoint --resume-from 0 checkpoint.h5
#+END_SRC
which displays
#+NAME: checkpoint-resume-from
#+BEGIN_SRC shell :exports results :results output
${PUGS_CHECKPOINT} --no-color --resume-from 0 checkpoint.h5
#+END_SRC
#+RESULTS: checkpoint-resume-from
Running again
#+BEGIN_SRC shell :exports code
pugs_checkpoint --info checkpoint.h5
#+END_SRC
now prints
#+NAME: checkpoint-info0
#+BEGIN_SRC shell :exports results :results output
${PUGS_CHECKPOINT} --no-color --info checkpoint.h5
#+END_SRC
#+RESULTS: checkpoint-info0
One notices that ~resuming_checkpoint~ now points to ~checkpoint_0~.
It is now possible to resume the execution using the command
#+BEGIN_SRC shell :exports code
pugs --resume checkpoint.h5
#+END_SRC
The output is now
#+NAME: checkpoint-resume0
#+BEGIN_SRC shell :exports results :results output
${PUGS} --no-exec-stat --no-preamble --no-color --threads=1 --resume checkpoint.h5
#+END_SRC
#+RESULTS: checkpoint-resume0
Observe that the random seed is reset to the value that was stored at
checkpoint 0.
Also each checkpoint whose number is greater than the resuming
checkpoint (here 0) are removed when the next checkpoint (here 1) is
written.
#+BEGIN_note
One may have noticed that the script file that is used for resuming is
actually stored in the ~checkpoint.h5~ file, the script is not a command
line argument.
By now, it is not possible to *simply* modify the script while
resuming. This is done on purpose since it remains unclear if such a
dangerous functionality should be make easy. Indeed allowing script
modifications may lead to undefined behaviors.
#+END_note
#+BEGIN_warning
The number of ~MPI~ ranks *cannot* be changed when resuming.
This is not likely to change since it is a very specific
functionality. However, a post-processing tool rewriting checkpoints
may be developed to achieve it, but it is not a priority.
#+END_warning
****** ~checkpoint_and_exit: void -> void~
This is an advanced version of the ~checkpoint~ function that cause a
clean ~exit~ of the code after writing a checkpoint.
The use case for this function is batch execution and is generally
triggered by the ~stop~ function.
We give here a simple example.
#+NAME: checkpoint-exit-function
#+BEGIN_SRC pugs :exports both :results output
for (let i:N, i = 1; i<=10; ++i) {
cout << "i = " << i << "\n";
if (i==5) {
checkpoint_and_exit();
}
}
#+END_SRC
The output is
#+RESULTS: checkpoint-exit-function
Then one can resume the execution as previously by
#+BEGIN_SRC shell :exports code
pugs --resume checkpoint.h5
#+END_SRC
It gives the output
The output is now
#+NAME: checkpoint-resume1
#+BEGIN_SRC shell :exports results :results output
${PUGS} --no-exec-stat --no-preamble --no-color --threads=1 --resume checkpoint.h5
#+END_SRC
#+RESULTS: checkpoint-resume1
# Clean-up
#+BEGIN_SRC shell :exports results :results none
/bin/rm -f checkpoint.h5
#+END_SRC
***** ~getAvailableModules: void -> string~
This function that is used in the preamble of this section is a
function that returns a ~string~ that contains the list of modules that
are available in the current version of ~pugs~.
***** ~getModuleInfo: string -> string~
This is a *very* helpful function. It lists *all* the available *builtin
functions* of a module and *all* the provided *types*. It is important to
remark that the name of the module is given as a ~string~ and not as is.
#+BEGIN_note
This function provides a very simple documentation of the
modules. Hopefully in future developments it would be really useful to
obtain a complete documentation of each function to offer a precise
online documentation.
#+END_note
In the following of the documentation we will run this function on all
modules.
#+NAME: get-module-info-core
#+BEGIN_SRC pugs :exports both :results output
cout << getModuleInfo("core") << "\n";
#+END_SRC
For the ~core~ module, it gives
#+RESULTS: get-module-info-core
This output is quite rustic, but it contains the main information: the
name of the function and its input and output sets.
#+BEGIN_warning
Observe that this function does not provide the list of operators that
are defined in the module (possibly associated to the defined types).
#+END_warning
***** ~getPugsBuildInfo: void -> string~
Sets the building information of the executable into a ~string~.
#+NAME: get-pugs-build-info
#+BEGIN_SRC pugs :exports both :results output
cout << getPugsBuildInfo() << "\n";
#+END_SRC
It gives for instance
#+RESULTS: get-pugs-build-info
***** ~getPugsVersion: void -> string~
Sets the ~pugs~ version of the executable into a ~string~.
#+NAME: get-pugs-version
#+BEGIN_SRC pugs :exports both :results output
cout << getPugsVersion() << "\n";
#+END_SRC
The output contains also ~git~ information
#+RESULTS: get-pugs-version
***** ~resetRandomSeed: void -> void~
At the start of ~pugs~, a *common* random seed is defined for all ~MPI~
processes. The seed is automatically written in the preamble of the
listing (in the documentation the ~--no-preamble~ option is used this is
why it is not displayed in examples).
The ~resetRandomSeed~ creates a *new shared* random and displays its value.
#+NAME: reset-random-seed
#+BEGIN_SRC pugs :exports both :results output
resetRandomSeed();
#+END_SRC
The output is
#+RESULTS: reset-random-seed
***** ~setRandomSeed: Z -> void~ <<set-random-seed>>
In order to reproduce exactly the same results of a calculation, it
can be interesting to set the random seed to some given value.
#+NAME: set-random-seed
#+BEGIN_SRC pugs :exports both :results output
setRandomSeed(123456789);
#+END_SRC
Running this example gives
#+RESULTS: set-random-seed
***** ~ofstream: string -> ostream~
This function is used to create an ~ostream~ that actually writes to the
file whose name is given by the ~string~ argument. One should notice
that the file is *created* at the function call. If a file with the same
name existed, it is *erased*.
#+NAME: ofstream-example
#+BEGIN_SRC pugs :exports both :results output
let fout:ostream, fout = ofstream("filename.txt");
fout << [1,2] << " is a vector of R^2\n";
#+END_SRC
Running this example produces no output but a file is created (in the
execution directory), with the name ~"filename.txt"~. Its content is
#+NAME: cat-filename-txt
#+BEGIN_SRC shell :exports results :results output
cat filename.txt
#+END_SRC
#+RESULTS: cat-filename.txt
Following ~C++~, the file is closed when the variable is destroyed.
*** The ~dev_utils~ module
This module contains utilities for ~pugs~ developers.
**** ~dev_utils~ provided functions
#+NAME: get-module-info-dev-utils
#+BEGIN_SRC pugs :exports both :results output
cout << getModuleInfo("dev_utils") << "\n";
#+END_SRC
#+RESULTS: get-module-info-dev-utils
***** ~getAST: void -> string~
Gets the abstract syntax tree (AST) that will be executed into a ~string~.
***** ~getFunctionAST: function -> string~
Gets the abstract syntax tree associated to a user function into a
~string~.
***** ~saveASTDot: string -> void~
Saves the AST of the script into the file (whose name is given as
argument) using the dot format.
***** Parallel checker utilities
Developing parallel applications, it is generally useful to be able to
check that results are reproducible bit-to-bit, whichever the number
of processes or threads are. Such a mechanism is provided by ~pugs~.
It works as follows. One first runs the code to build a reference (by
now the reference cannot be build using ~MPI~ parallelism), then a
second run will compare the obtained results and reference ones. If
ghost value differs, a warning is printed, stating that the data is
not synchronized (this may not be a mistake) and if non-ghost value
differs the execution stops indicating an error.
Parallel checking points can be placed directly in the source code
(see developer documentation), but for convenience, one can also check
parallelism directly in the scripting language.
The provided functions are
- ~parallel_check: Vh*string -> void~
- ~parallel_check: item_array*string -> void~
- ~parallel_check: item_value*string -> void~
- ~parallel_check: sub_item_array*string -> void~
- ~parallel_check: sub_item_value*string -> void~
They allow to check parallelism for ~Vh~ functions, ~item_value~,
~item_array~, ~sub_item_value~ or ~sub_item_arrays~ variables. These types
are discussed below when describing the ~mesh~ and ~schemes~ module (see
sections [[mesh]] and [[scheme]]). The second argument is just a string that is just
used as a tag to ease the reading of the output.
To create the reference, one launches the code as
#+BEGIN_SRC shell :exports source
./pugs example.pgs --parallel-checker-file="/tmp/my-ref.h5"
#+END_SRC
Observe that the option ~--parallel-checker-file="/tmp/my-ref.h5"~ is
not mandatory. If not specified, the reference is created in the
current directory in the file ~parallel_checker.h5~. Also observe that
this command runs the code in sequential (in the sense that there is
no message passing). In that case parallel checker runs in write mode
automatically. To force the read mode (for comparison), one can use
the ~--parallel-checker-mode=read~ option
#+BEGIN_SRC shell :exports source
./pugs example.pgs \
--parallel-checker-file="/tmp/my-ref.h5" \
--parallel-checker-mode=read
#+END_SRC
Running the code in parallel automatically triggers the
read/comparison mode.
#+BEGIN_SRC shell :exports source
mpirun -n 3 pugs example.pgs --parallel-checker-file="/tmp/my-ref.h5"
#+END_SRC
*** The ~math~ module
The ~math~ module is a small utility module that provides a set of
standard mathematical functions.
#+NAME: get-module-info-math
#+BEGIN_SRC pugs :exports both :results output
cout << getModuleInfo("math") << "\n";
#+END_SRC
For the ~math~ module, it gives
#+RESULTS: get-module-info-math
There is not much to say. One can see that this module does not
provide new types of data. All the classical mathematical functions
follow their ~C++~ counterparts. One can see that there are integer
versions (type ~Z~) of ~abs~, ~min~ and ~max~ functions.
The ~dot~ function family provides the dot product for vectors of
$\mathbb{R}^1$, $\mathbb{R}^2$ and $\mathbb{R}^3$.
#+NAME: dot-examples
#+BEGIN_SRC pugs :exports both :results output
import math;
cout << "([1],[2]) = " << dot([1],[2]) << "\n";
cout << "([1,2],[3,4]) = " << dot([1,2],[3,4]) << "\n";
cout << "([1,2,3],[4,5,6]) = " << dot([1,2,3],[4,5,6]) << "\n";
#+END_SRC
The output is
#+RESULTS: dot-examples
A set of ~det~ functions is defined to get the determinant of matrices
of $\mathbb{R}^{1\times1}$, $\mathbb{R}^{2\times2}$ and
$\mathbb{R}^{3\times3}$.
#+NAME: det-examples
#+BEGIN_SRC pugs :exports both :results output
import math;
cout << "det([[1.2]]) = " << det([[1.2]]) << "\n";
cout << "det([[1,2],[3,4]]) = " << det([[1,2],[3,4]]) << "\n";
cout << "det([[1,2,3],[4,5,6],[7,8,9]]) = "
<< det([[1,2,3],[4,5,6],[7,8,9]]) << "\n";
#+END_SRC
The output is
#+RESULTS: det-examples
The ~trace~ functions compute the trace of matrices of
$\mathbb{R}^{1\times1}$, $\mathbb{R}^{2\times2}$ and $\mathbb{R}^{3\times3}$.
#+NAME: trace-examples
#+BEGIN_SRC pugs :exports both :results output
import math;
cout << "trace([[1.2]]) = " << trace([[1.2]]) << "\n";
cout << "trace([[1,2],[3,4]]) = " << trace([[1,2],[3,4]]) << "\n";
cout << "trace([[1,2,3],[4,5,6],[7,8,9]]) = "
<< trace([[1,2,3],[4,5,6],[7,8,9]]) << "\n";
#+END_SRC
The output is
#+RESULTS: trace-examples
The ~doubleDot~ functions compute the double-dot ($A:B = \mathrm{tr}(A^T
B)$) of two matrices of $\mathbb{R}^{1\times1}$, $\mathbb{R}^{2\times2}$ and
$\mathbb{R}^{3\times3}$.
#+NAME: double-dot-examples
#+BEGIN_SRC pugs :exports both :results output
import math;
cout << "doubleDot([[1.2]],[[2.3]]) = " << doubleDot([[1.2]],[[2.3]]) << "\n";
cout << "doubleDot([[1,2],[3,4]],[[3,6],[2,5]]) = " << doubleDot([[1,2],[3,4]],[[3,6],[2,5]]) << "\n";
cout << "doubleDot([[1,2,3],[4,5,6],[7,8,9]], [[5,2,1],[4,2,8],[2,6,2]]) = "
<< doubleDot([[1,2,3],[4,5,6],[7,8,9]], [[5,2,1],[4,2,8],[2,6,2]]) << "\n";
#+END_SRC
The output is
#+RESULTS: double-dot-examples
Also, one can compute inverses of $\mathbb{R}^{1\times1}$,
$\mathbb{R}^{2\times2}$ and $\mathbb{R}^{3\times3}$ matrices using the
~inverse~ function set.
#+NAME: inverse-examples
#+BEGIN_SRC pugs :exports both :results output
import math;
cout << "inverse([[1.2]]) = " << inverse([[1.2]]) << "\n";
cout << "inverse([[1,2],[3,4]]) = " << inverse([[1,2],[3,4]]) << "\n";
cout << "inverse([[3,2,1],[5,6,4],[7,8,9]])\n = "
<< inverse([[3,2,1],[5,6,4],[7,8,9]]) << "\n";
#+END_SRC
The output is
#+RESULTS: inverse-examples
Transpose of matrices of $\mathbb{R}^{1\times1}$,
$\mathbb{R}^{2\times2}$ and $\mathbb{R}^{3\times3}$ are obtained using the
~transpose~ functions.
#+NAME: transpose-examples
#+BEGIN_SRC pugs :exports both :results output
import math;
cout << "transpose([[1.2]]) = " << transpose([[1.2]]) << "\n";
cout << "transpose([[1,2],[3,4]]) = " << transpose([[1,2],[3,4]]) << "\n";
cout << "transpose([[3,2,1],[5,6,4],[7,8,9]]) = "
<< transpose([[3,2,1],[5,6,4],[7,8,9]]) << "\n";
#+END_SRC
The output is
#+RESULTS: transpose-examples
#+BEGIN_note
Observe that the use of a proper rounding or truncation function is
the right way to convert a real value to an integer one. Available
rounding or truncation functions are ~ceil~, ~floor~, ~round~ and ~trunc~. See
their ~C++~ man pages for details.
#+END_note
#+BEGIN_note
Let us comment on the use of the ~pow~ function. Actually one can wonder
why we did not use a syntax like ~x^y~? The reason is that if
mathematically ${x^y}^z = x^{(y^z)}$, many softwares treat it (by mistake)
as ${(x^y)}^z$. Thus, using the ~pow~ function avoids any confusion.
#+END_note
*** The ~mesh~ module<<mesh>>
This is an important module. It provides mesh utility tools.
#+NAME: get-module-info-mesh
#+BEGIN_SRC pugs :exports both :results output
cout << getModuleInfo("mesh") << "\n";
#+END_SRC
#+RESULTS: get-module-info-mesh
**** ~mesh~ provided types
***** ~boundary~
The ~boundary~ type is a boundary descriptor: it refers to a boundary of
a ~mesh~ that is either designated by an integer or by a ~string~.
A ~boundary~ can designate a set of nodes, edges or faces. The ~boundary~
(descriptor) itself is not related to any ~mesh~, thus the nature of the
~boundary~ is made precise when it is used with a particular ~mesh~.
#+BEGIN_warning
A ~boundary~ *cannot* be used to refer to an interface (/ie/ an inner set of
items).
#+END_warning
***** ~zone~
Following the same idea, a ~zone~ is a descriptor of a set of cells. It
can be either defined by an integer or by a ~string~. Its meaning is
made precise when it is associated with a ~mesh~.
***** ~mesh~
The type ~mesh~ is an *abstract* type that is used to store meshes. A
variable of that type can refer for instance to unstructured meshes of
dimension 1, 2 or 3.
The following binary operator is provided.
| ~ostream <<~ allowed expression type |
|------------------------------------|
| ~mesh~ |
It enables to write mesh information to an ~ostream~, here is a simple
example.
#+NAME: ostream-mesh-example
#+BEGIN_SRC pugs :exports both :results output
import mesh;
let m:mesh, m = cartesianMesh(0,[1,1,2], (2,3,2));
cout << m << "\n";
#+END_SRC
Running this script generates the following output.
#+RESULTS: ostream-mesh-example
The information produced concerns
- the dimension of the connectivity,
- the number of cells and their references,
- the number of faces and their references,
- the number of edges and their references,
- and the number of nodes and their references.
***** ~item_type~
This type is used to designate kinds of items (cell, face, edge or node).
***** ~item_value~ and ~item_array~
The types ~item_value~ and ~item_array~ are abstract types use to
designate values or arrays of values defined on each entities of a
~mesh~. Actually, these set of values (or arrays) are not associated to
the mesh itself but to the *connectivity*. The values on the entities
can be of type ~B~, ~N~, ~Z~, ~R~, ~R^1~, ~R^2~, ~R^3~, ~R^1x1~, ~R^2x2~ and
~R^3x3~. Entities themselves can be cells, faces, edges or nodes.
These variables are used to pass data from one function to another.
#+BEGIN_warning
By now, no mathematical operation is defined on ~item_value~ or
~item_array~ variables.
#+END_warning
Moreover, ~item_value~ variables can be post processed. Observe that
edge or face values cannot be post processed since neither ~VTK~ nor
~Gnuplot~ can handle these data. Also ~item_raray~ can only be post
processed if array per item contain scalar data ( ~B~, ~N~, ~Z~, ~R~).
***** ~sub_item_value~ and ~sub_item_arrays~
These abstract type handles values or arrays defined on the sub items
of the items of a ~mesh~. Similarly these sets of values or arrays are
attached to a *connectivity* and not to a mesh. The values associated to
the sub items can be of type ~B~, ~N~, ~Z~, ~R~, ~R^1~, ~R^2~, ~R^3~, ~R^1x1~, ~R^2x2~
or ~R^3x3~. An example of ~sub_item_value~ is the $\mathbf{C}_{jr}$ vectors
which are defined at each node of each cell.
These variables are used to pass data from one function to
another. They cannot be post processed.
**** ~mesh~ provided functions
***** Boundary descriptor functions
****** ~boundaryName: string -> boundary~
Sets a boundary descriptor to a boundary name
#+BEGIN_SRC pugs :exports both :results none
import mesh;
let b:boundary, b = boundaryName("boundary_name");
#+END_SRC
****** ~boundaryTag: Z -> boundary~
Creates a boundary descriptor to a boundary reference
#+BEGIN_SRC pugs :exports both :results none
import mesh;
let b:boundary, b = boundaryTag(12);
#+END_SRC
***** Zone descriptor functions
****** ~zoneName: string -> zone~
Creates a zone descriptor from a ~string~ name
#+BEGIN_SRC pugs :exports both :results none
import mesh;
let z:zone, z = zoneName("zone_name");
#+END_SRC
****** ~zoneTag: Z -> zone~
Associates a zone descriptor from a zone tag
#+BEGIN_SRC pugs :exports both :results none
import mesh;
let z:zone, z = zoneTag(5);
#+END_SRC
***** ~cartesianMesh: Rˆd*Rˆd*(N) -> mesh~ (with $d\in\{1, 2, 3\}$)
Creates a cartesian mesh of dimension $d\in\{1, 2, 3\}$. The produced
cartesian grid is aligned with the axis and made of identical cells.
The first two arguments are two opposite corners of the box (or of the
segment in 1d) and the list of natural integers (type ~(N)~) sets the
number of *cells* in each direction. Thus the size of the list ~N~ is $d$.
For instance one can write:
#+BEGIN_SRC pugs :exports both :results none
import mesh;
let m1d:mesh, m1d = cartesianMesh([0], [1], 10);
let m2d:mesh, m2d = cartesianMesh([0,0], [1,1], (3,6));
let m3d:mesh, m3d = cartesianMesh([-1,-1,-1], [1,1,1], (3,2,4));
#+END_SRC
- The ~m1d~ variable contains a uniform grid of $]0,1[$ made of 10 cells.
- The ~m2d~ variable refers to a uniform grid of $]0,1[^2$, made of 3
cells in the $x$ direction, and of 6 cells along the $y$ axis.
- The ~m3d~ variable designates a mesh of $]-1,1[^3$ made of $3\times2\times4$
cells.
***** ~readGmsh: string -> mesh~
Reads a ~mesh~ from a file.
#+BEGIN_warning
The file must conform to the mesh format ~msh2~.
#+END_warning
#+BEGIN_SRC shell :exports results :results none
cat << EOF > hybrid-2d.geo
//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {1, 0, 0, 1.0};
//+
Point(3) = {1, 1, 0, 1.0};
//+
Point(4) = {0, 1, 0, 1.0};
//+
Point(5) = {2, 0, 0, 1.0};
//+
Point(6) = {2, 1, 0, 1.0};
//+
Line(1) = {4, 1};
//+
Line(2) = {2, 3};
//+
Line(3) = {5, 6};
//+
Line(4) = {3, 4};
//+
Line(5) = {6, 3};
//+
Line(6) = {1, 2};
//+
Line(7) = {2, 5};
//+
Curve Loop(1) = {1, 6, 2, 4};
//+
Surface(1) = {1};
//+
Curve Loop(2) = {-2, 7, 3, 5};
//+
Surface(2) = {2};
//+
Characteristic Length {4, 3, 6, 5, 2, 1} = 0.3;
//+
Recombine Surface {1};
//+
Physical Curve("XMIN") = {1};
//+
Physical Curve("XMAX") = {3};
//+
Physical Curve("YMAX") = {4, 5};
//+
Physical Curve("YMIN") = {6, 7};
//+
Physical Curve("INTERFACE") = {2};
//+
Physical Surface("LEFT") = {1};
//+
Physical Surface("RIGHT") = {2};
//+
Physical Point("XMINYMIN") = {1};
//+
Physical Point("XMINYMAX") = {4};
//+
Physical Point("XMAXYMIN") = {5};
//+
Physical Point("XMAXYMAX") = {6};
#+END_SRC
Here is an example of ~gmsh~ use that produces a mesh at format ~msh2~.
#+BEGIN_SRC shell :exports both :results none
gmsh -2 hybrid-2d.geo -format msh2
#+END_SRC
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import writer;
let m:mesh, m = readGmsh("hybrid-2d.msh");
write_mesh(gnuplot_writer("hybrid-2d"),m);
#+END_SRC
The ~mesh~ is represented in Figure [[fig:gmsh-hybrid-2d]].
#+NAME: gmsh-hybrid-2d-img
#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/gmsh-hybrid-2d.png")
reset
unset grid
unset border
unset key
unset xtics
unset ytics
set terminal png truecolor enhanced size 960,480
plot '<(sed "" $PUGS_SOURCE_DIR/doc/hybrid-2d.gnu)' w l
#+END_SRC
#+CAPTION: The mesh that was read from the file ~hydrid-2d.msh~ and then saved to the ~gnuplot~ format
#+NAME: fig:gmsh-hybrid-2d
#+ATTR_LATEX: :width 0.38\textwidth
#+ATTR_HTML: :width 300px;
#+RESULTS: gmsh-hybrid-2d-img
***** ~diamondDual: mesh -> mesh~
This function creates the diamond dual ~mesh~ of a primal ~mesh~.
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import writer;
let primal_mesh:mesh, primal_mesh = readGmsh("hybrid-2d.msh");
let dual_mesh:mesh, dual_mesh = diamondDual(primal_mesh);
write_mesh(gnuplot_writer("diamond"), dual_mesh);
#+END_SRC
The diamond dual mesh is defined by joining the nodes of the faces to
the center of the adjacent cells of the primal mesh.
The ~mesh~ is represented in Figure [[fig:gmsh-hybrid-2d]].
#+NAME: diamond-dual-img
#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/diamond-dual.png")
reset
unset grid
unset border
unset key
unset xtics
unset ytics
set terminal png truecolor enhanced size 960,480
plot '<(sed "" $PUGS_SOURCE_DIR/doc/hybrid-2d.gnu)' lt rgb "green" w l, '<(sed "" $PUGS_SOURCE_DIR/doc/diamond.gnu)' lt rgb "black" w l
#+END_SRC
#+CAPTION: The primal mesh in green and the diamond dual mesh in black
#+NAME: fig:diamond-dual
#+ATTR_LATEX: :width 0.38\textwidth
#+ATTR_HTML: :width 300px;
#+RESULTS: diamond-dual-img
#+BEGIN_note
The mesh storage mechanism in ~pugs~ is such that the diamond dual mesh
is built only once. This means that if one writes for instance
#+BEGIN_SRC pugs :exports both :results none
import mesh;
let primal_mesh:mesh, primal_mesh = readGmsh("hybrid-2d.msh");
let dual_mesh_1:mesh, dual_mesh_1 = diamondDual(primal_mesh);
let dual_mesh_2:mesh, dual_mesh_2 = diamondDual(primal_mesh);
#+END_SRC
then the variables ~dual_mesh_1~ and ~dual_mesh_2~ will refer to the same
~mesh~ in memory.
#+END_note
#+BEGIN_warning
The diamond dual mesh construction is not yet implemented in parallel
#+END_warning
***** ~medianDual: mesh -> mesh~
This function creates the median dual ~mesh~ of a primal ~mesh~.
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import writer;
let primal_mesh:mesh, primal_mesh = readGmsh("hybrid-2d.msh");
let dual_mesh:mesh, dual_mesh = medianDual(primal_mesh);
write_mesh(gnuplot_writer("median"), dual_mesh);
#+END_SRC
The median dual mesh is defined by joining the centers of the faces to
the centers of the adjacent cells of the primal mesh.
The ~mesh~ is represented in Figure [[fig:gmsh-hybrid-2d]].
#+NAME: median-dual-img
#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/median-dual.png")
reset
unset grid
unset border
unset key
unset xtics
unset ytics
set terminal png truecolor enhanced size 960,480
plot '<(sed "" $PUGS_SOURCE_DIR/doc/hybrid-2d.gnu)' lt rgb "green" w l, '<(sed "" $PUGS_SOURCE_DIR/doc/median.gnu)' lt rgb "black" w l
#+END_SRC
#+CAPTION: The primal mesh in green and the median dual mesh in black
#+NAME: fig:median-dual
#+ATTR_LATEX: :width 0.38\textwidth
#+ATTR_HTML: :width 300px;
#+RESULTS: median-dual-img
#+BEGIN_note
In ~pugs~, the storage mechanisms of median dual meshes follow the same
rules as the diamond dual meshes. As long as the primary mesh lives
and as long as the median dual mesh is referred, it is kept in memory,
thus constructed only once.
#+END_note
#+BEGIN_warning
The median dual mesh construction is not yet implemented in 3d and not
available in parallel
#+END_warning
***** Item types
The following functions are used to designate a specific ~item_type~
- ~cell: void -> item_type~
- ~face: void -> item_type~
- ~edge: void -> item_type~
- ~node: void -> item_type~
***** ~interpolate: mesh*item_type*function -> item_value~
This function is used to create an ~item_value~ by interpolating a user
function to a given ~item_type~ using a mesh. The interpolation takes
place at the center of the items.
#+BEGIN_SRC pugs :exports both :results none
import mesh;
let m:mesh, m = cartesianMesh([0,0], [1,1], (10,10));
let f:R^2 -> R, x -> 2*x[0]-x[1];
let node_values:item_value, node_values = interpolate(m, node(), f);
let face_values:item_value, face_values = interpolate(m, face(), f);
let cell_values:item_value, cell_values = interpolate(m, cell(), f);
#+END_SRC
In this example, we set three arrays defined at all nodes, all the
faces and all the cells of the underlying connectivity of the mesh ~m~.
- ~node_values~ interpolates ~f~ at the nodes of ~m~,
- ~face_values~ interpolates ~f~ at the centers of the faces of ~m~,
- ~cell_values~ interpolates ~f~ at the centers of the cells of ~m~.
***** ~load_balance:mesh -> mesh~
Creates a new partition of a mesh in parallel (~MPI~).
#+BEGIN_SRC pugs :exports both :results none
import mesh;
let m1:mesh, m1 = cartesianMesh([0,0], [1,1], (10,10));
let m2:mesh, m2 = load_balance(m1);
#+END_SRC
In this example the mesh ~m2~ is a new parallel partition of the mesh
~m1~. This implies that the two meshes, while being identical from the
mathematical point of view, are now different (and are not defined on
a same connectivity: the connectivity is actually balanced).
Discrete functions defined on ~m1~ are *not* implicitly transferred on the
mesh ~m2~! The function ~load_balance:(Vh)->(Vh)~ must be used to perform
mesh and data load balancing.
#+BEGIN_note
To simplify development of large calculations. A sequential call of
this function has the same effect: the mesh and its connectivity are
copied!
#+END_note
***** ~transform: mesh*function -> mesh~
This function allows to compute a new mesh as the transformation of a
given mesh by displacing its nodes through a user defined function.
For a mesh of dimension $d$, the transformation must be a function
$\mathbb{R}^d \to\mathbb{R}^d$.
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import math;
import writer;
let m0:mesh, m0 = cartesianMesh([0,0], [1,1], (10,10));
let pi_over_3: R, pi_over_3 = acos(-1)/3;
let t: R^2 -> R^2, x -> [(0.2+0.8*x[0])*cos(pi_over_3*x[1]), (0.2+0.8*x[0])*sin(pi_over_3*x[1])];
let m1: mesh, m1 = transform(m0, t);
write_mesh(gnuplot_writer("transformed"), m1);
#+END_SRC
#+BEGIN_note
One should keep in mind that the mesh produced by the ~transform~
function *shares* the same connectivity as the original mesh. This means
that in ~pugs~ internals, there is only one connectivity object for
these two meshes.
#+END_note
The result of the previous script is given in Figure
[[fig:transformed]]. They all share the same connectivity in memory.
#+NAME: transformed-img
#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/transformed.png")
reset
unset grid
unset border
unset key
unset xtics
unset ytics
set square
set terminal png truecolor enhanced size 640,640
plot '<(sed "" $PUGS_SOURCE_DIR/doc/transformed.gnu)' w l
#+END_SRC
#+CAPTION: The transformed unit square
#+NAME: fig:transformed
#+ATTR_LATEX: :width 0.38\textwidth
#+ATTR_HTML: :width 300px;
#+RESULTS: transformed-img
***** ~relax: mesh*mesh*R -> mesh~ <<relax-function>>
This function is a simple utility that computes a ~mesh~ as the /mean/ of
two other meshes that share the same connectivity. The coordinates of
the vertices of the relaxed mesh $\mathcal{M}_2$, are given by
\begin{equation*}
\forall r\in\mathcal{R},\quad\mathbf{x}_r^{\mathcal{M}_2} = (1-\theta) \mathbf{x}_r^{\mathcal{M}_0} + \theta \mathbf{x}_r^{\mathcal{M}_1}.
\end{equation*}
#+BEGIN_warning
The connectivity must be the *same in memory*, this means that
constructing two identical meshes with /equivalent/ connectivity is not
allowed.
#+END_warning
Thus for instance, the following code
#+NAME: relax-with-similar-connecticities
#+BEGIN_SRC pugs-error :exports both :results output
import mesh;
let m1:mesh, m1 = cartesianMesh(0, [1,1], (10,10));
let m2:mesh, m2 = cartesianMesh(0, [1,1], (10,10));
let m3:mesh, m3 = relax(m1,m2,0.3);
#+END_SRC
produces the runtime error
#+results: relax-with-similar-connecticities
A proper use is
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import math;
import writer;
let m0:mesh, m0 = cartesianMesh([0,0], [1,1], (10,10));
write_mesh(gnuplot_writer("relax_example_m0"), m0);
let t: R^2 -> R^2, x -> [x[0]+0.25*x[1]*x[1], x[1]+0.3*sin(x[0])];
let m1: mesh, m1 = transform(m0, t);
write_mesh(gnuplot_writer("relax_example_m1"), m1);
let m2: mesh, m2 = relax(m0, m1, 0.3);
write_mesh(gnuplot_writer("relax_example_m2"), m2);
#+END_SRC
In this example, the relaxation parameter is set to $\theta=0.3$. The
different meshes produced in this example are displayed in Figure
[[fig:relax]].
#+NAME: relax-img
#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/relax.png")
reset
unset grid
unset border
unset key
unset xtics
unset ytics
set square
set terminal png truecolor enhanced size 640,640
plot '<(sed "" $PUGS_SOURCE_DIR/doc/relax_example_m0.gnu)' lt rgb "black" w l, '<(sed "" $PUGS_SOURCE_DIR/doc/relax_example_m1.gnu)' lt rgb "blue" w l, '<(sed "" $PUGS_SOURCE_DIR/doc/relax_example_m2.gnu)' lw 2 lt rgb "red" w l
#+END_SRC
#+CAPTION: Example of meshes relaxation. The relaxed mesh $\mathcal{M}_2$ (red) and the original meshes in black ($\mathcal{M}_0$) and blue ($\mathcal{M}_1$).
#+NAME: fig:relax
#+ATTR_LATEX: :width 0.38\textwidth
#+ATTR_HTML: :width 300px;
#+RESULTS: relax-img
This function is mainly useful to reduce the displacement of nodes
when using the ~randomizeMesh~ functions (see section [[scheme]]) for
instance.
*** The ~scheme~ module<<scheme>>
This module provides a lot of numerical tools.
#+NAME: get-module-info-scheme
#+BEGIN_SRC pugs :exports both :results output
cout << getModuleInfo("scheme") << "\n";
#+END_SRC
#+RESULTS: get-module-info-scheme
#+BEGIN_warning
This module is very large and will be split in smaller ones in the future
#+END_warning
This module provides various types and functions. It also provides a
set of operators overloading.
**** ~scheme~ provided types
***** ~Vh~
The ~Vh~ type is the /abstract/ type of variables that refer to *discrete
functions*.
****** $\mathbb{P}_0$ scalar functions
These functions are defined on a ~mesh~ and have a constant value in
each cell.
The type of values in each cells for a $\mathbb{P}_0$ function can be
$\mathbb{R}$, $\mathbb{R}^1$, $\mathbb{R}^2$, $\mathbb{R}^3$,
$\mathbb{R}^{1\times1}$, $\mathbb{R}^{2\times2}$ or $\mathbb{R}^{3\times3}$.
For simplicity in this document, we denote these specific
$\mathbb{P}_0$ types as $\mathbb{P}_0(\mathbb{R})$,
$\mathbb{P}_0(\mathbb{R}^1)$, $\mathbb{P}_0(\mathbb{R}^2)$,
$\mathbb{P}_0(\mathbb{R}^3)$, $\mathbb{P}_0(\mathbb{R}^{1\times1})$,
$\mathbb{P}_0(\mathbb{R}^{2\times2})$ or
$\mathbb{P}_0(\mathbb{R}^{3\times3})$.
****** $\mathbb{P}_0$ vector functions
Additionally to scalar values per cell, one can define vectors of real
values per cell. The size of the vectors is the same for all
cells. This kind of variables is useful to define mass fractions for
instance.
Again for convenience, these types are denoted as
$\vec{\mathbb{P}}_0(\mathbb{R})$.
***** ~Vh_type~
The ~Vh_type~ allows to describe a type of discretization. The available
types of discretization are
- ~P0~ for $\mathbb{P}_0(\mathbb{R})$, $\mathbb{P}_0(\mathbb{R}^1)$,
$\mathbb{P}_0(\mathbb{R}^2)$, $\mathbb{P}_0(\mathbb{R}^3)$,
$\mathbb{P}_0(\mathbb{R}^{1\times1})$, $\mathbb{P}_0(\mathbb{R}^{2\times2})$
or $\mathbb{P}_0(\mathbb{R}^{3\times3})$.
- ~P0Vector~ for $\vec{\mathbb{P}}_0(\mathbb{R})$
***** ~boundary_condition~
This type is used to describe boundary conditions.
***** ~quadrature~
This type is used to describe quadrature types.
**** ~scheme~ provided functions
The ~scheme~ module provides a lot of functions, we categorize their
description.
***** Mathematical functions
****** ~Vh -> Vh~
The majority of these functions are defined for
$\mathbb{P}_0(\mathbb{R})$ data and the return value is also a
$\mathbb{P}_0(\mathbb{R})$ function. The value is simply the
application of the function to the cell values.
Here is the list of these functions
- ~abs: Vh -> Vh~
- ~acos: Vh -> Vh~
- ~acosh: Vh -> Vh~
- ~asin: Vh -> Vh~
- ~asinh: Vh -> Vh~
- ~atan: Vh -> Vh~
- ~atanh: Vh -> Vh~
- ~cos: Vh -> Vh~
- ~cosh: Vh -> Vh~
- ~exp: Vh -> Vh~
- ~log: Vh -> Vh~
- ~sin: Vh -> Vh~
- ~sinh: Vh -> Vh~
- ~sqrt: Vh -> Vh~
- ~tan: Vh -> Vh~
- ~tanh: Vh -> Vh~
A few functions are defined for $\mathbb{P}_0(\mathbb{R^{}}^{1\times1})$,
$\mathbb{P}_0(\mathbb{R^{}}^{2\times2})$ and
$\mathbb{P}_0(\mathbb{R^{}}^{3\times3})$ data and the return value is a
$\mathbb{P}_0(\mathbb{R})$ function. The value is simply the
application of the function to the cell values.
- ~det: Vh -> Vh~
- ~trace: Vh -> Vh~
Also, functions are defined for $\mathbb{P}_0(\mathbb{R}^{d\times d})$
data and the return value is a $\mathbb{P}_0(\mathbb{R}^{d\times d})$
function, with $d\in\{1,2,3\}$. The value is simply the application of
the function to the cell values.
- ~inverse: Vh -> Vh~
- ~transpose: Vh -> Vh~
****** ~Vh*Vh -> Vh~
These functions are defined for $\mathbb{P}_0(\mathbb{R})$ data and the
return value is also a $\mathbb{P}_0(\mathbb{R})$ function. These
functions require that the two arguments are defined on the *same
mesh*. The result is obtained by applying the function cell by cell.
Here is the function list
- ~atan2: Vh*Vh -> Vh~
- ~max: Vh*Vh -> Vh~
- ~min: Vh*Vh -> Vh~
- ~pow: Vh*Vh -> Vh~
Let us mention another function that applies to
$\mathbb{P}_0(\mathbb{R}^1)$, $\mathbb{P}_0(\mathbb{R}^2)$,
$\mathbb{P}_0(\mathbb{R}^3)$ and to $\vec{\mathbb{P}}_0(\mathbb{R})$
vector functions.
- ~dot: Vh*Vh -> Vh~
This function requires that both arguments are defined on the same
mesh and have the same data type. The result is a
$\mathbb{P}_0(\mathbb{R})$ function.
Finally the equivalent exists for discrete functions of matrices and
applies to $\mathbb{P}_0(\mathbb{R}^1x1)$, $\mathbb{P}_0(\mathbb{R}^2x2)$,
$\mathbb{P}_0(\mathbb{R}^3x3)$
- ~doubleDot: Vh*Vh -> Vh~
Both arguments must be defined on the same mesh and have the same data
type. The result is a $\mathbb{P}_0(\mathbb{R})$ function.
****** ~R*Vh -> Vh~ and ~Vh*R -> Vh~
These functions are defined for $\mathbb{P}_0(\mathbb{R})$ data and the
return value is also a $\mathbb{P}_0(\mathbb{R})$ function.
The following functions can be applied using a scalar ~R~ and a ~Vh~
operand.
- ~atan2: Vh*R -> Vh~
- ~atan2: R*Vh -> Vh~
- ~max: Vh*R -> Vh~
- ~max: R*Vh -> Vh~
- ~min: Vh*R -> Vh~
- ~min: R*Vh -> Vh~
- ~pow: Vh*R -> Vh~
- ~pow: R*Vh -> Vh~
****** ~R^1*Vh -> Vh~ and ~Vh*R^1 -> Vh~
These functions are defined for $\mathbb{P}_0(\mathbb{R}^1)$ data and the
return value is a $\mathbb{P}_0(\mathbb{R})$ function.
The following functions
- ~dot: Rˆ1*Vh -> Vh~
- ~dot: Vh*Rˆ1 -> Vh~
****** ~R^2*Vh -> Vh~ and ~Vh*R^2 -> Vh~
These functions are defined for $\mathbb{P}_0(\mathbb{R}^2)$ data and the
return value is a $\mathbb{P}_0(\mathbb{R})$ function.
The following functions
- ~dot: Rˆ2*Vh -> Vh~
- ~dot: Vh*Rˆ2 -> Vh~
****** ~R^3*Vh -> Vh~ and ~Vh*R^3 -> Vh~
These functions are defined for $\mathbb{P}_0(\mathbb{R}^3)$ data and the
return value is a $\mathbb{P}_0(\mathbb{R})$ function.
The following functions
- ~dot: Rˆ3*Vh -> Vh~
- ~dot: Vh*Rˆ3 -> Vh~
****** ~R^1x1*Vh -> Vh~ and ~Vh*R^1x1 -> Vh~
These functions are defined for $\mathbb{P}_0(\mathbb{R}^1x1)$ data and the
return value is a $\mathbb{P}_0(\mathbb{R})$ function.
The following functions
- ~doubleDot: Rˆ1x1*Vh -> Vh~
- ~doubleDot: Vh*Rˆ1x1 -> Vh~
****** ~R^2x2*Vh -> Vh~ and ~Vh*R^2x2 -> Vh~
These functions are defined for $\mathbb{P}_0(\mathbb{R}^2x2)$ data and the
return value is a $\mathbb{P}_0(\mathbb{R})$ function.
The following functions
- ~doubleDot: Rˆ2x2*Vh -> Vh~
- ~doubleDot: Vh*Rˆ2x2 -> Vh~
****** ~R^3x3*Vh -> Vh~ and ~Vh*R^3x3 -> Vh~
These functions are defined for $\mathbb{P}_0(\mathbb{R}^3x3)$ data and the
return value is a $\mathbb{P}_0(\mathbb{R})$ function.
The following functions
- ~doubleDot: Rˆ3x3*Vh -> Vh~
- ~doubleDot: Vh*Rˆ3x3 -> Vh~
****** ~Vh -> R~
These functions are defined for $\mathbb{P}_0(\mathbb{R})$ data and the
return value is real ($\mathbb{R}$).
The following functions
- ~min: Vh -> R~\\
returns the minimal value of all the cell values
- ~max: Vh -> R~\\
returns the maximal value of all the cell values
- ~integral_of_R: Vh -> R~\\
computes the integral of a $\mathbb{P}_0(\mathbb{R})$ discrete
function $f$. If $f_j \in\mathbb{R}$ denotes the value of $f$ in the
cell $j$, the return value is $$\sum_{j\in\mathcal{J}}V_j f_j,$$ where $V_j$
denotes the volume of the cell $j$.
- ~sum_of_R: Vh -> R~\\
computes the sum of the cell values of a $\mathbb{P}_0(\mathbb{R})$
discrete function $f$. If $f_j \in\mathbb{R}$ denotes the value of
$f$ in the cell $j$, the return value is $$\sum_{j\in\mathcal{J}} f_j.$$
****** ~Vh -> R^1~
These functions are defined for $\mathbb{P}_0(\mathbb{R}^1)$ data and
the return value is a $\mathbb{R}^1$ vector.
- ~integral_of_R1: Vh -> R^1~\\
computes the integral of a $\mathbb{P}_0(\mathbb{R}^1)$ discrete
function $\mathbf{u}$. If $\mathbf{u}_j \in\mathbb{R}^1$ denotes the
value of $\mathbf{u}$ in the cell $j$, the return value is
$$\sum_{j\in\mathcal{J}}V_j \mathbf{u}_j,$$ where $V_j$ denotes the volume of
the cell $j$.
- ~sum_of_R1: Vh -> R^1~\\
computes the sum of the cell values of a $\mathbb{P}_0(\mathbb{R}^1)$
discrete function $\mathbf{u}$. If $\mathbf{u}_j \in\mathbb{R}^1$
denotes the value of $\mathbf{u}$ in the cell $j$, the return value
is $$\sum_{j\in\mathcal{J}} \mathbf{u}_j.$$
****** ~Vh -> R^2~
These functions are defined for $\mathbb{P}_0(\mathbb{R}^2)$ data and
the return value is a $\mathbb{R}^2$ vector.
- ~integral_of_R2: Vh -> R^2~\\
computes the integral of a $\mathbb{P}_0(\mathbb{R}^2)$ discrete
function $\mathbf{u}$. If $\mathbf{u}_j \in\mathbb{R}^2$ denotes the
value of $\mathbf{u}$ in the cell $j$, the return value is
$$\sum_{j\in\mathcal{J}}V_j \mathbf{u}_j,$$ where $V_j$ denotes the volume of
the cell $j$.
- ~sum_of_R2: Vh -> R^2~\\
computes the sum of the cell values of a $\mathbb{P}_0(\mathbb{R}^2)$
discrete function $\mathbf{u}$. If $\mathbf{u}_j \in\mathbb{R}^2$
denotes the value of $\mathbf{u}$ in the cell $j$, the return value
is $$\sum_{j\in\mathcal{J}} \mathbf{u}_j.$$
****** ~Vh -> R^3~
These functions are defined for $\mathbb{P}_0(\mathbb{R}^3)$ data and
the return value is a $\mathbb{R}^3$ vector.
- ~integral_of_R3: Vh -> R^3~\\
computes the integral of a $\mathbb{P}_0(\mathbb{R}^3)$ discrete
function $\mathbf{u}$. If $\mathbf{u}_j \in\mathbb{R}^3$ denotes the
value of $\mathbf{u}$ in the cell $j$, the return value is
$$\sum_{j\in\mathcal{J}}V_j \mathbf{u}_j,$$ where $V_j$ denotes the volume of
the cell $j$.
- ~sum_of_R3: Vh -> R^3~\\
computes the sum of the cell values of a $\mathbb{P}_0(\mathbb{R}^3)$
discrete function $\mathbf{u}$. If $\mathbf{u}_j \in\mathbb{R}^3$
denotes the value of $\mathbf{u}$ in the cell $j$, the return value
is $$\sum_{j\in\mathcal{J}} \mathbf{u}_j.$$
****** ~Vh -> R^1x1~
These functions are defined for $\mathbb{P}_0(\mathbb{R}^{1\times1})$ data
and the return value is an $\mathbb{R}^{1\times1}$ matrix.
- ~integral_of_R1x1: Vh -> R^1x1~\\
computes the integral of a $\mathbb{P}_0(\mathbb{R}^{1\times1})$
discrete function $A$. If $A_j \in\mathbb{R}^{1\times1}$ denotes the
value of $A$ in the cell $j$, the return value is
$$\sum_{j\in\mathcal{J}}V_j A_j,$$ where $V_j$ denotes the volume of the cell
$j$.
- ~sum_of_R1x1: Vh -> R^1x1~\\
computes the sum of the cell values of a
$\mathbb{P}_0(\mathbb{R}^{1\times1})$ discrete function $A$. If $A_j
\in\mathbb{R}^{1\times1}$ denotes the value of $A$ in the cell $j$, the
return value is $$\sum_{j\in\mathcal{J}} A_j.$$
****** ~Vh -> R^2x2~
These functions are defined for $\mathbb{P}_0(\mathbb{R}^{2\times2})$ data
and the return value is an $\mathbb{R}^{2\times2}$ matrix.
- ~integral_of_R2x2: Vh -> R^2x2~\\
computes the integral of a $\mathbb{P}_0(\mathbb{R}^{2\times2})$
discrete function $A$. If $A_j \in\mathbb{R}^{2\times2}$ denotes the
value of $A$ in the cell $j$, the return value is
$$\sum_{j\in\mathcal{J}}V_j A_j,$$ where $V_j$ denotes the volume of the cell
$j$.
- ~sum_of_R2x2: Vh -> R^2x2~\\
computes the sum of the cell values of a
$\mathbb{P}_0(\mathbb{R}^{2\times2})$ discrete function $A$. If $A_j
\in\mathbb{R}^{2\times2}$ denotes the value of $A$ in the cell $j$, the
return value is $$\sum_{j\in\mathcal{J}} A_j.$$
****** ~Vh -> R^3x3~
These functions are defined for $\mathbb{P}_0(\mathbb{R}^{3\times3})$ data
and the return value is an $\mathbb{R}^{3\times3}$ matrix.
- ~integral_of_R3x3: Vh -> R^3x3~\\
computes the integral of a $\mathbb{P}_0(\mathbb{R}^{3\times3})$
discrete function $A$. If $A_j \in\mathbb{R}^{3\times3}$ denotes the
value of $A$ in the cell $j$, the return value is
$$\sum_{j\in\mathcal{J}}V_j A_j,$$ where $V_j$ denotes the volume of the cell
$j$.
- ~sum_of_R3x3: Vh -> R^3x3~\\
computes the sum of the cell values of a
$\mathbb{P}_0(\mathbb{R}^{3\times3})$ discrete function $A$. If $A_j
\in\mathbb{R}^{3\times3}$ denotes the value of $A$ in the cell $j$, the
return value is $$\sum_{j\in\mathcal{J}} A_j.$$
***** Utilities
****** ~sum_of_Vh: Vh -> Vh~
This function it computes the sum of all the components of
$\vec{\mathbb{P}}_0(\mathbb{R})$ function and stores the result into a
$\mathbb{P}_0(\mathbb{R})$.
****** ~vectorize: (Vh) -> Vh~
This function creates a $\vec{\mathbb{P}}_0(\mathbb{R})$ function using
a tuple of $\mathbb{P}_0(\mathbb{R})$.
***** Interpolation and integration functions
These functions are ways to define discrete functions from analytic
data.
****** ~interpolate: mesh*Vh_type*(function) -> Vh~
This function takes a ~mesh~, a type of discrete function (~Vh_type~) and
a list of user functions as arguments. It returns a $\mathbb{P}_0$
function defined at the mesh.
All the user functions of the list must be defined on $\mathbb{R}^d$
for a mesh of dimension $d$.
There are several situations according to the ~Vh_type~.
******* ~P0~
In that case the list of user functions *must* reduce to a *single* user
function.
The codomain (or image space) of the user function defines the type of
the returned discrete function.
Let us give an example.
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
let m:mesh, m = cartesianMesh([0,0], [1,1], (10,10));
let u:R^2 -> R^3, x -> [x[0], 2*x[1], x[1]-x[0]];
let uh:Vh, uh = interpolate(m, P0(), u);
#+END_SRC
In this example the discrete function ~uh~ is a
$\mathbb{P}_0(\mathbb{R}^3)$ function defined on a 2d ~mesh~. The ~P0()~
function returns the type of interpolation ($\mathbb{P}_0$).
#+BEGIN_note
In the case of $\mathbb{P}_0$ interpolation, the function is evaluated
at the mass center of the mesh cells $$\forall j\in\mathcal{J},\quad
\mathbf{x}_j=\frac{1}{V_j}\int_j \mathbf{x}.$$
#+END_note
******* ~P0Vector~
In that case the codomain of each user function in the list must be a
real function (values in $\mathbb{R}$). The instruction will define a
$\vec{\mathbb{P}}_0(\mathbb{R})$.
A first example is
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
let m:mesh, m = cartesianMesh([0,0], [1,1], (10,10));
let f:R^2 -> R, x -> 2*x[0]-x[1];
let fh:Vh, fh = interpolate(m, P0Vector(), f);
#+END_SRC
The obtained ~fh~ is a vector $\vec{\mathbb{P}}_0(\mathbb{R})$ where each
vector (in each cell) is of dimension 1.
The number of scalar user functions sets the size of the discrete
vector function.
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
import math;
let m:mesh, m = cartesianMesh([0,0], [1,1], (10,10));
let f0:R^2 -> R, x -> 2*x[0]-x[1];
let f1:R^2 -> R, x -> x[0]*x[1]-2;
let f2:R^2 -> R, x -> 2*dot(x,x);
let fh:Vh, fh = interpolate(m, P0Vector(), (f0,f1,f2));
#+END_SRC
Here we defined a $\vec{\mathbb{P}}_0(\mathbb{R})$ discrete function of
dimension 3.
****** ~interpolate: mesh*(zone)*Vh_type*(function) -> Vh~
This function is similar to the previous function. The additional
parameter, the ~zone~ list is used to define the cells where the user
function (or the user function list) is interpolated. For cells that
are not in the ~zone~ list, the discrete function is set to the value
$0$.
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
let m:mesh, m = readGmsh("hybrid-2d.msh");
let f:R^2 -> R, x -> 2*x[0]-x[1];
let fh:Vh, fh = interpolate(m, (zoneName("LEFT"), zoneName("RIGHT")), P0(), f);
let fh_l:Vh, fh_l = interpolate(m, zoneName("LEFT"), P0(), f);
#+END_SRC
In this example, we define two discrete functions. ~fh~ is defined as
the interpolation of the function ~f~ in *all* cells of the mesh since the
mesh is partitioned into two zones labeled ~LEFT~ and ~RIGHT~. The
discrete function ~fh_l~ has exactly the same values in the ~LEFT~ region,
but is $0$ in the cells that belong to ~RIGHT~.
For completeness, we give an example in the case of ~P0Vector~.
The number of scalar user functions sets the size of the discrete
vector function.
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
import math;
let m:mesh, m = readGmsh("hybrid-2d.msh");
let f0:R^2 -> R, x -> 2*x[0]-x[1];
let f1:R^2 -> R, x -> x[0]*x[1]-2;
let f2:R^2 -> R, x -> 2*dot(x,x);
let fh:Vh, fh = interpolate(m, zoneName("RIGHT"), P0Vector(), (f0,f1,f2));
#+END_SRC
Here we defined a $\vec{\mathbb{P}}_0(\mathbb{R})$ discrete function of
dimension 3.
****** ~integrate: mesh*quadrature*function -> Vh~ <<integrate-classic>>
This function integrates the given user function in each cell of a
~mesh~ using a prescribed ~quadrature~. The result is of type
$\mathbb{P}_0(\mathbb{R})$, $\mathbb{P}_0(\mathbb{R}^1)$,
$\mathbb{P}_0(\mathbb{R}^2)$, $\mathbb{P}_0(\mathbb{R}^3)$,
$\mathbb{P}_0(\mathbb{R}^{1\times1})$, $\mathbb{P}_0(\mathbb{R}^{2\times2})$
or $\mathbb{P}_0(\mathbb{R}^{3\times3})$, according to the user function
codomain.
Let us consider the following example
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
import math;
let m:mesh, m = readGmsh("hybrid-2d.msh");
let u:R^2 -> R^2, x -> [cos(x[0]), sin(x[1])];
let U:Vh, U = integrate(m, Gauss(5), u);
#+END_SRC
Here, for each cell $j$, the value of the discrete function
$\mathbf{U}_j$ is computed using a Gauss quadrature formula that is
exact for polynomials of degree $5$, $\mathbf{U}_j \approx\int_j
\mathbf{u}$. More details about quadrature formula will be given
below.
Thus if one wants to project $\mathbf{u}$ on
$\mathbb{P}_0(\mathbb{R}^2)$ to the sixth order, one can modify the
previous script by writing
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
import math;
let m:mesh, m = readGmsh("hybrid-2d.msh");
let u:R^2 -> R^2, x -> [cos(x[0]), sin(x[1])];
let uh:Vh, uh = (1/cell_volume(m)) * integrate(m, Gauss(5), u);
#+END_SRC
The function ~cell_volume: mesh -> Vh~ creates a
$\mathbb{P}_0(\mathbb{R})$ function whose values are the volume of the
cells.
****** ~integrate: mesh*quadrature*Vh_type*(function) -> Vh~ <<integrate-P1-vector>>
This function behaves similarly, the user function list size defines
the dimension of the vector value of the produced
$\vec{\mathbb{P}}_0(\mathbb{R})$ discrete function. Actually the
~Vh_type~ parameter is there to allow the construction of
$\vec{\mathbb{P}}_0(\mathbb{R})$ of dimension 1 (since passing a single
function as a list of user function, the previous function
[[integrate-classic]], would be used).
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
import math;
let m:mesh, m = readGmsh("hybrid-2d.msh");
let u:R^2 -> R, x -> cos(x[0]*x[1]);
let v:R^2 -> R, x -> x[0]*x[1];
let w:R^2 -> R, x -> x[0]+x[1];
let u1h:Vh, u1h = integrate(m, Gauss(5), P0Vector(), u);
let u2h:Vh, u2h = integrate(m, Gauss(5), u);
let u3h:Vh, u3h = integrate(m, Gauss(5), P0Vector(), (u,v,w));
#+END_SRC
In this example one observes the difference between ~u1h~ and ~u2h~. The
first one is a $\vec{\mathbb{P}}_0(\mathbb{R})$ where the vector size
is 1, and ~u2h~ is a $\mathbb{P}_0(\mathbb{R})$.
****** ~integrate: mesh*(zone)*quadrature*function -> Vh~
This function is an enhancement of the function defined in
[[integrate-classic]]. It allows to specify a ~zone~ list which defines the
set of cells where the integration is operated.
#+BEGIN_SRC shell :exports results :results none
cat << EOF > zones-1d.geo
//+
Point(1) = {-1, 0, 0, 0.01};
//+
Point(2) = {-0.3, 0, 0, 0.01};
//+
Point(3) = {0.3, 0, 0, 0.01};
//+
Point(4) = {1, 0, 0, 0.01};
//+
Line(1) = {1, 2};
//+
Line(2) = {2, 3};
//+
Line(3) = {3, 4};
//+
Physical Point("XMIN") = {1};
//+
Physical Point("XMAX") = {4};
//+
Physical Point("INTERFACE1") = {2};
//+
Physical Point("INTERFACE2") = {3};
//+
Physical Curve("LEFT") = {1};
//+
Physical Curve("MIDDLE") = {2};
//+
Physical Curve("RIGHT") = {3};
#+END_SRC
#+BEGIN_SRC shell :exports results :results none
gmsh -1 zones-1d.geo -format msh2
#+END_SRC
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
import math;
import writer;
let pi:R, pi = acos(-1);
let m:mesh, m = readGmsh("zones-1d.msh");
let f:R^1 -> R, x -> sin(2*pi*x[0]);
let fh:Vh,
fh = (1/cell_volume(m))
,* integrate(m, (zoneName("LEFT"), zoneName("RIGHT")), Gauss(5), f);
write(gnuplot_1d_writer("zone_integrate"), name_output(fh, "f"));
#+END_SRC
In this example, the ~mesh~ provided in the file ~zones-1d.msh~ is a 1d
~mesh~ of $]-1,1[$ made of $200$ cells that is partitioned into 3
connected subdomains. The zones corresponding to these 3 subdomains
are named ~LEFT~ for $]-1,-0.3[$, ~MIDDLE~ for $]-0.3, 0.3[$ and ~RIGHT~ for
$]0.3,1[$. The result is displayed in Figure [[fig:zone-integrate-1d]]. In
the ~MIDDLE~ region, cell values are set to 0.
#+NAME: zone-integrate-1d-img
#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/zone-integrate-1d.png")
reset
set grid
set border
unset key
set xtics
set ytics
set terminal png truecolor enhanced size 960,480
plot '<(sed "" $PUGS_SOURCE_DIR/doc/zone_integrate.gnu)' lw 2 w l
#+END_SRC
#+CAPTION: $L^1$ projection of order 6 of the function $\sin(2\pi x)$ on the zones ~LEFT~ and ~RIGHT~. The values in the zone ~MIDDLE~ are set to $0$.
#+NAME: fig:zone-integrate-1d
#+ATTR_LATEX: :width 0.38\textwidth
#+ATTR_HTML: :width 300px;
#+RESULTS: zone-integrate-1d-img
****** ~integrate: mesh*(zone)*quadrature*Vh_type*(function) -> Vh~
This function behaves essentially as the function described in
paragraph [[integrate-P1-vector]], it also adds the possibility to choose
sets of cells where to integrate the list of user functions.
***** Random mesh generators
For numerical tests purpose, it is often useful to create meshes with
random vertices positions. This is the aim of the functions that are
described in this section. These functions share some properties.
- The generated mesh is always suitable for calculations in the sense
that cell volumes are warrantied to be positive.
- Generated cells can be non-convex.
- One has to specify boundary conditions to drive the mesh
displacement on boundaries.
- The obtained mesh does not depend on parallelism: it is exactly the
same whichever the number of ~MPI~ processes is. It only depends on
the random seed (see in paragraph [[set-random-seed]] how to set the
random seed to obtain the exact same mesh through different runs).
****** ~randomizeMesh: mesh*(boundary_condition) -> mesh~
This function creates a random mesh by moving the nodes of a given
~mesh~ and a list of ~boundary_condition~.
The supported boundary conditions are the following:
- ~fixed~: the vertices of the ~boundary~ cannot be displaced
- ~axis~: the vertices are allowed to be displaced in the direction of
the ~boundary~
- ~symmetry~: the vertices are displaced in the plane formed by the
~boundary~
One should refer to the section [[boundary-condition-descriptor]] for a
documentation of the boundary condition descriptors.
#+BEGIN_note
Let us make precise these boundary conditions behavior
- In dimension 1, ~fixed~, ~axis~ and ~symmetry~ boundary conditions have
the same effect.
- In dimension 2, ~axis~ and ~symmetry~ behave similarly. Thus, boundaries
supporting this kind of boundary conditions *must* form *straight*
lines.
- In dimension 3, boundaries describing ~axis~ conditions *must* be
*straight* lines, and boundaries describing ~symmetry~ conditions *must*
be *planar*.
If a boundary does not satisfy geometrical requirements, ~pugs~ produces
a runtime error.
#+END_note
Let us consider a simple example
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
import writer;
setRandomSeed(123456789); // not required
let m:mesh, m = cartesianMesh([-1,-1],[1,1], (20,20));
let bc_list:(boundary_condition),
bc_list = ((fixed(boundaryName("XMINYMIN"))),
(fixed(boundaryName("XMINYMAX"))),
(fixed(boundaryName("XMAXYMIN"))),
(fixed(boundaryName("XMAXYMAX"))),
(symmetry(boundaryName("XMIN"))),
(symmetry(boundaryName("XMAX"))),
(symmetry(boundaryName("YMIN"))),
(symmetry(boundaryName("YMAX"))));
m = randomizeMesh(m, bc_list);
write_mesh(gnuplot_writer("random-mesh"), m);
#+END_SRC
Running this script one gets the ~mesh~ displayed in Figure
[[fig:random-mesh]]. To reduce the vertices displacement, one can use the
~relax~ function, see section [[relax-function]].
#+NAME: random-mesh-img
#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/random-mesh.png")
reset
unset grid
unset border
unset key
unset xtics
unset ytics
set square
set terminal png truecolor enhanced size 640,640
plot '<(sed "" $PUGS_SOURCE_DIR/doc/random-mesh.gnu)' w l
#+END_SRC
#+CAPTION: The obtained random mesh
#+NAME: fig:random-mesh
#+ATTR_LATEX: :width 0.38\textwidth
#+ATTR_HTML: :width 300px;
#+RESULTS: random-mesh-img
****** ~randomizeMesh: mesh*(boundary_condition)*function -> mesh~
This function is a variation of the previous one. It allows
to additionally provide a characteristic function that designates the
vertices that can be displaced.
The characteristic function *must* be a function of $\mathbb{R}^d
\to\mathbb{B}$ for a ~mesh~ of dimension $d$.
Here is a modification of the previous example, where the random
displacement is allowed for $x<2y$.
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
import writer;
import math;
setRandomSeed(123456789); // not required
let m:mesh, m = cartesianMesh([-1,-1],[1,1], (20,20));
let bc_list:(boundary_condition),
bc_list = ((fixed(boundaryName("XMINYMIN"))),
(fixed(boundaryName("XMINYMAX"))),
(fixed(boundaryName("XMAXYMIN"))),
(fixed(boundaryName("XMAXYMAX"))),
(symmetry(boundaryName("XMIN"))),
(symmetry(boundaryName("XMAX"))),
(symmetry(boundaryName("YMIN"))),
(symmetry(boundaryName("YMAX"))));
let chi:R^2 -> B, x -> x[0]<2*x[1];
m = randomizeMesh(m, bc_list, chi);
write_mesh(gnuplot_writer("random-mesh-chi"), m);
#+END_SRC
Running this script one gets the ~mesh~ displayed in Figure
[[fig:random-mesh-chi]].
#+NAME: random-mesh-chi-img
#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/random-mesh-chi.png")
reset
unset grid
unset border
unset key
unset xtics
unset ytics
set square
set terminal png truecolor enhanced size 640,640
plot '<(sed "" $PUGS_SOURCE_DIR/doc/random-mesh-chi.gnu)' w l
#+END_SRC
#+CAPTION: Random mesh with characteristic function. The original mesh is unchanged for $x \ge 2 y$.
#+NAME: fig:random-mesh-chi
#+ATTR_LATEX: :width 0.38\textwidth
#+ATTR_HTML: :width 300px;
#+RESULTS: random-mesh-chi-img
#+BEGIN_note
Since we set the random seed to the same value in both cases (with or
without using the characteristic function $\chi$), node displacements are
the same $\forall r\in\mathcal{R}$ such that $\chi(\mathbf{x}_r)$, see Figures
[[fig:random-mesh]] and [[fig:random-mesh-chi]].
This allows for instance to build a random mesh step-by-step.
#+END_note
***** Boundary condition descriptors <<boundary-condition-descriptor>>
In this paragraph, we discuss the set of boundary condition
descriptors that are provided by the ~scheme~ module.
#+BEGIN_note
These functions provide *descriptors*, these are not related to a
particular implementation. The way they are used in different
functions may vary.
#+END_note
#+BEGIN_note
There are three kinds of boundaries supported by ~pugs~, boundaries made
- of sets of nodes,
- of sets of edges, or
- of sets of faces.
In dimension 1, nodes, edges and faces denote the same entities. In
dimension 2, edges and faces refer the same entities.
#+END_note
#+BEGIN_note
~pugs~ integrates some automatic mechanisms to translate boundary types
in other ones.
For instance, if an algorithm or a method requires a set of nodes to
set some numerical treatment, it can be deduced from a set of faces.
Obviously, these conversions can be meaningless, for instance, if one
expects a *line* in 3d, it cannot be defined by a set of faces. ~pugs~
will forbid this kind of conversion at runtime.
#+END_note
#+BEGIN_note
When a method expects a *straight* line or a *planar* surface, ~pugs~ checks
that the given boundary is actually *straight* or *planar*.
#+END_note
#+BEGIN_note
~pugs~ checks that boundaries do not contain /inner/ items.
#+END_note
We regroup the boundary condition descriptor functions according to
their arguments
****** ~boundary -> boundary_condition~
- ~axis: boundary -> boundary_condition~\\
The boundary denotes a *straight* line of the mesh
- ~fixed: boundary -> boundary_condition~
- ~symmetry: boundary -> boundary_condition~\\
The boundary denotes a *planar* surface of the mesh
****** ~boundary*function -> boundary_condition~
The provided user function type depends on the numerical method of
function that utilizes the boundary condition.
- ~pressure: boundary*function -> boundary_condition~
- ~velocity: boundary*function -> boundary_condition~
***** Quadrature formulas descriptors
~pugs~ provides a quite complete set of quadrature formulas that can be
chosen inside the script.
#+BEGIN_warning
While ~pugs~ is written to deal with general polygonal and polyhedral
meshes, quadrature formulas are not yet implemented on the general
elements. Supported elements are triangles, quadrangles, tetrahedra,
pyramids, prisms and hexahedra.
For completeness, let us state that they are also available for
diamond cells that can be decomposed into two tetrahedra or two
pyramids.
#+END_warning
****** ~Gauss: N -> quadrature~
This function returns the quadrature descriptor associated to Gauss
formulas for the given ~N~.
In the following table, we summarize the *maximal degree* quadrature
(exact for polynomials of a given degree) that are available in ~pugs~
for various elements.
| element type | max. degree |
|------------------------+-------------|
| segment | 23 |
|------------------------+-------------|
| triangle | 20 |
| square | 21 |
|------------------------+-------------|
| tetrahedron | 20 |
| pyramid (square basis) | 20 |
| prism (triangle basis) | 20 |
| cube | 21 |
****** ~GaussLegendre: N -> quadrature~
Gets the Gauss-Legendre quadrature descriptor exact for polynomials of
degree given in argument.
The maximum allowed degree is 23.
For 2 or 3-dimensional elements, Gauss-Legendre formulas are defined
by tensorization. Conform transformations are used to map the cube
$]-1,1[^d$ to supported elements.
****** ~GaussLobatto: N -> quadrature~
Gets the Gauss-Lobatto quadrature descriptor exact for polynomials of
degree given in argument.
The maximum allowed degree is "only" 13.
For 2 or 3-dimensional elements, Gauss-Lobatto formulas are defined by
tensorization. Conform transformations are used to map the cube
$]-1,1[^d$ to supported elements.
***** ~lagrangian: mesh*Vh -> Vh~
This function is a special function whose purpose is to transport
lagrangian quantities from one mesh to the other. Obviously, this
function makes a lot of sense in the case of lagrangian calculations.
This is a zero cost function, since discrete functions are *constant*
in ~pugs~, it consists in associating the data of the given discrete
function to the provided ~mesh~. In other words, the underlying array of
values is shared by the two discrete functions, which are associated
to different meshes.
A good example of the use of this kind of function is mass fractions.
***** ~load_balance: (Vh) -> (Vh)~
This function performs a parallel load balancing of a list of ~Vh~
variables. All the input variables must be defined on a same mesh. The
return list consists is the balance variables using the input
ordering. The new (balanced) mesh is carried out by the new variables
and can be retrieved using the ~get_mesh: Vh -> mesh~ function.
#+BEGIN_warning
Using this function in a sequential calculation has the same behavior:
the discrete functions, the mesh and its connectivity are copied. This
is intended to simplify development of large calculations on simpler
tests, and for consistency.
#+END_warning
***** Numerical methods
We describe rapidly two functions that embed numerical methods. These
are in some sense /models/ that are used to test evolution of ~pugs~
itself and can be used as examples.
#+BEGIN_warning
These functions will become obsolete (soon?), since another interface
to numerical methods is in preparation.
#+END_warning
The functions are
- ~eucclhyd_solver: Vh*Vh*Vh*Vh*Vh*(boundary_condition)*R -> mesh*Vh*Vh*Vh~
- ~glace_solver: Vh*Vh*Vh*Vh*Vh*(boundary_condition)*R -> mesh*Vh*Vh*Vh~
Both functions share the same interface. The arguments are provided in
this order:
- the mass density $\rho$ of type $\mathbb{P}_0(\mathbb{R})$,
- the velocity $\mathbf{u}$ of type $\mathbb{P}_0(\mathbb{R}^d)$ in
dimension $d$,
- the total energy density $E$ of type $\mathbb{P}_0(\mathbb{R})$,
- the sound speed $c$ of type $\mathbb{P}_0(\mathbb{R})$,
- the pressure $p$ of type $\mathbb{P}_0(\mathbb{R})$,
- a list of boundary conditions,
- and a time step of type $\mathbb{R}$.
Observe that ~pugs~ checks the types of the parameters and that all
discrete functions are defined on the same mesh.
The functions return a compound type made of
- the new ~mesh~ $\mathcal{M}$,
- the new mass density $\rho$ of type $\mathbb{P}_0(\mathbb{R})$ defined
on $\mathcal{M}$,
- the new velocity $\mathbf{u}$ of type $\mathbb{P}_0(\mathbb{R}^d)$ in
dimension $d$, defined on the mesh $\mathcal{M}$,
- and the new total energy density $E$ of type
$\mathbb{P}_0(\mathbb{R})$, defined on the mesh $\mathcal{M}$.
The time step can be calculated through the ~acoustic_dt: Vh -> R~
function.
**** Operators overloading for ~Vh~ <<Vh-operators>>
The ~scheme~ module provides overload of binary operators. Since ~Vh~ is
an abstract type, some operators may not be defined (or allowed) for
concrete types.
***** Unary operators
The only supported unary operators for ~Vh~ are given in the following
table
| operator | description |
|----------+----------------------|
| ~+~ | plus unary operator |
| ~-~ | minus unary operator |
We recall that the unary ~+~ operator is a convenience operator that has
no effect.
***** Binary operators
The supported binary operators for ~Vh~ data types are arithmetic
operators.
#+begin_src latex :results drawer :exports results
\begin{equation*}
\left|
\begin{array}{rl}
\mathtt{+}:& \mathbb{V}_h \times \mathbb{V}_h \to \mathbb{V}_h\\
\mathtt{-}:& \mathbb{V}_h \times \mathbb{V}_h \to \mathbb{V}_h\\
\mathtt{*}:& \mathbb{V}_h \times \mathbb{V}_h \to \mathbb{V}_h\\
\mathtt{/}:& \mathbb{V}_h \times \mathbb{V}_h \to \mathbb{V}_h
\end{array}
\right.
\end{equation*}
#+end_src
Observe that in the case of $\vec{\mathbb{P}}_0(\mathbb{R})$, the only
available operators are ~+~ and ~-~.
#+BEGIN_note
In the case of $\mathbb{P}_0$ functions, an operator is defined as soon
as it is defined for the value type.
#+END_note
For instance, one can multiply a $\mathbb{P}_0(\mathbb{R}^{2\times2})$
discrete function by a $\mathbb{P}_0(\mathbb{R}^2)$. The result is then
a $\mathbb{P}_0(\mathbb{R}^2)$ discrete function.
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
let m:mesh, m = cartesianMesh(0,[1,1],(10,10));
let A:R^2 -> R^2x2, x -> [[2*x[0], x[1]],[-x[0], 3*x[1]]];
let u:R^2 -> R^2, x -> 2*[x[0], x[1]*x[0]];
let Ah:Vh, Ah = interpolate(m, P0(), A);
let uh:Vh, uh = interpolate(m, P0(), u);
let Auh:Vh, Auh = Ah*uh;
#+END_SRC
Another illustration is: trying to add $\mathbb{P}_0(\mathbb{R})$ and
$\mathbb{P}_0(\mathbb{R}^1)$
#+NAME: invalid-Vh-sum-type
#+BEGIN_SRC pugs-error :exports both :results output
import mesh;
import scheme;
let m:mesh, m = cartesianMesh([0],[1],10);
let f:R^1 -> R, x -> 2*x[0];
let f1:R^1 -> R^1, x -> 2*x;
let fh:Vh, fh = interpolate(m, P0(), f);
let f1h:Vh, f1h = interpolate(m, P0(), f1);
fh+f1h;
#+END_SRC
produces the following error
#+results: invalid-Vh-sum-type
****** Additional ~+~ and ~-~ operators
#+begin_src latex :results drawer :exports results
\begin{equation*}
\forall \mathbb{S}, \mathbb{S}_2 \in \{\mathbb{B},\mathbb{N},\mathbb{Z},\mathbb{R},\mathbb{R}^1,\mathbb{R}^2,\mathbb{R}^3,\mathbb{R}^{1\times1},\mathbb{R}^{2\times2},\mathbb{R}^{3\times3}\},
\quad
\left|
\begin{array}{rl}
\mathtt{+}:& \mathbb{S} \times \mathbb{V}_h \to \mathbb{V}_h\\
\mathtt{-}:& \mathbb{S} \times \mathbb{V}_h \to \mathbb{V}_h\\
\mathtt{+}:& \mathbb{V}_h \times \mathbb{S} \to \mathbb{V}_h\\
\mathtt{-}:& \mathbb{V}_h \times \mathbb{S} \to \mathbb{V}_h
\end{array}
\right.
\end{equation*}
#+end_src
Let us consider the following example
#+NAME: substract-mean-value-to-Vh
#+BEGIN_SRC pugs :exports both :results output
import mesh;
import scheme;
import math;
let m:mesh, m = cartesianMesh(0,[0.3,1.1],(5,15));
let mesh_volume:R, mesh_volume = sum_of_R(cell_volume(m));
let u:R^2 -> R, x -> 2*dot(x,x);
let uh:Vh, uh = interpolate(m, P0(), u);
let uh0:Vh, uh0 = uh - (integral_of_R(uh) / mesh_volume);
cout << "integral(uh) = " << integral_of_R(uh) << "\n";
cout << "integral(uh0) = " << integral_of_R(uh0) << "\n";
#+END_SRC
Here we subtract the mean value of a discrete function.
#+results: substract-mean-value-to-Vh
****** Additional ~*~ operators
The following constructions are allowed for ~*~ operator.
#+begin_src latex :results drawer :exports results
\begin{equation*}
\forall \mathbb{S} \in \{\mathbb{B},\mathbb{N},\mathbb{Z},\mathbb{R},\mathbb{R}^{1\times1},\mathbb{R}^{2\times2},\mathbb{R}^{3\times3}\},
\quad
\mathtt{*}: \mathbb{S} \times \mathbb{V}_h \to \mathbb{V}_h.
\end{equation*}
#+end_src
Obviously, if $\mathbb{S}=\mathbb{R}^{d\times d}$, for $d\in\{1,2,3\}$,
the right operand must have a compatible type, for instance
$\mathbb{P}_0(\mathbb{R}^{d\times d})$ or $\mathbb{P}_0(\mathbb{R}^d)$.
Additionally these operators are defined
#+begin_src latex :results drawer :exports results
\begin{equation*}
\forall \mathbb{S} \in \{\mathbb{B},\mathbb{N},\mathbb{Z},\mathbb{R},\mathbb{R}^1,\mathbb{R}^2,\mathbb{R}^3,\mathbb{R}^{1\times1},\mathbb{R}^{2\times2},\mathbb{R}^{3\times3}\},
\quad
\mathtt{*}: \mathbb{V}_h \times \mathbb{S}\to \mathbb{V}_h.
\end{equation*}
#+end_src
Following the logic, if for instance the right operand is an ~R^2~
expression, the left operand ~Vh~ must be of type
$\mathbb{P}_0(\mathbb{R})$ or $\mathbb{P}_0(\mathbb{R}^{2\times 2})$.
****** Additional ~/~ operators
The following operators are defined
#+begin_src latex :results drawer :exports results
\begin{equation*}
\forall \mathbb{S} \in \{\mathbb{N},\mathbb{Z},\mathbb{R}\},
\quad
\mathtt{/}: \mathbb{S} \times \mathbb{V}_h \to \mathbb{V}_h.
\end{equation*}
#+end_src
*** The ~linear_solver~ module
This module provides the following functions
#+NAME: get-module-info-linear-solver
#+BEGIN_SRC pugs :exports both :results output
cout << getModuleInfo("linear_solver") << "\n";
#+END_SRC
#+RESULTS: get-module-info-linear-solver
**** ~linear_solver~ provided functions
The set of provided functions is used to define the global behavior of
~pugs~ linear system solver.
***** Utility functions
An important function is
****** ~getLSAvailable: void -> string~
This shows the available options or libraries. It depends on the
compilation options of the code.
#+NAME: get-ls-available-example
#+BEGIN_SRC pugs :exports both :results output
import linear_solver;
cout << getLSAvailable() << "\n";
#+END_SRC
#+results: get-ls-available-example
****** ~getLSOptions: void -> string~
This function shows the current tuning of the global linear solver
#+NAME: get-ls-options-example
#+BEGIN_SRC pugs :exports both :results output
import linear_solver;
cout << getLSOptions() << "\n";
#+END_SRC
#+results: get-ls-options-example
***** Tuning functions
****** ~setLSEpsilon: R -> void~
Sets the relative epsilon criterion for iterative methods.
****** ~setLSLibrary: string -> void~
Selects the library to use.
****** ~setLSMaxIter: N -> void~
Sets the maximum number of iterations.
****** ~setLSMethod: string -> void~
Sets the method to solve linear systems.
****** ~setLSPrecond: string -> void~
Sets the preconditioner type.
****** ~setLSVerbosity: B -> void~
Sets verbosity mode to ~true~ or ~false~ for linear solvers.
*** The ~writer~ module
This module provides functionalities to write numerical results to
files for post processing.
It provides the following functions and types.
#+NAME: get-module-info-writer
#+BEGIN_SRC pugs :exports both :results output
cout << getModuleInfo("writer") << "\n";
#+END_SRC
#+RESULTS: get-module-info-writer
**** ~writer~ provided types
***** ~output~
This type is used to describe output of discrete functions (of type
~Vh~). Mainly, it associates a name to the discrete function.
#+BEGIN_note
At writing, ~pugs~ checks that the names in an ~output~ list are all
different.
#+END_note
***** ~writer~
Variables of this type manage outputs: which format is used and
possibly the writing policy. This policy sets for instance the time
period for time-dependent post processing.
**** ~writer~ provided functions
***** ~name_output: Vh*string -> output~
This function gives a name to a discrete function.
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
import writer;
let m:mesh, m = cartesianMesh([0], [1], 100);
let identity:R^1 -> R, x -> x[0];
let xh:Vh, xh = interpolate(m, P0(), identity);
let x2h:Vh, x2h = xh*xh;
write(gnuplot_writer("writer-example1"), name_output(x2h, "x_square"));
#+END_SRC
See section [[gnuplot-writers]] for details about these writers family.
#+NAME: writer-example1-img
#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/writer-example1.png")
reset
set grid
set border
unset key
set xtics
set ytics
set square
set terminal png truecolor enhanced size 640,640
plot '<(sed "" $PUGS_SOURCE_DIR/doc/writer-example1.gnu)' lw 2 w l
#+END_SRC
#+CAPTION: Simple illustration of the output.
#+NAME: fig:writer-example1
#+ATTR_LATEX: :width 0.38\textwidth
#+ATTR_HTML: :width 300px;
#+RESULTS: writer-example1-img
#+BEGIN_warning
It can be convenient to define variables or lists of ~output~. However,
one has to pay attention of the following fact. When declaring an
output variable, a reference *to the data* of the discrete function is
done. It is not a reference to the language variable (or
expression).
#+END_warning
Let us illustrate it by an important second example.
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
import writer;
import math;
let m:mesh, m = cartesianMesh([0], [1], 100);
let identity:R^1 -> R, x -> x[0];
let fh:Vh, fh = interpolate(m, P0(), identity);
let gh:Vh, gh = fh*fh;
let pi:R, pi = acos(-1);
let output_list:(output),
output_list = (name_output(fh, "f"),
name_output(gh, "g"),
name_output(sin(pi*fh), "sin"));
let abs_cos:R^1 -> R, x -> abs(cos(pi*x[0]));
fh = interpolate(m, P0(), abs_cos); // fh now refers to another function
write(gnuplot_writer("writer-example2"), output_list);
output_list = (name_output(fh, "f"),
name_output(gh, "g"));
write(gnuplot_writer("writer-example2-2"), output_list);
#+END_SRC
Running this code produces the gnuplot file displayed in Figure
[[fig:writer-example2]]. One sees that ~f~ is the $\mathbb{P}_0(\mathbb{R})$
function corresponding to the function $x \to x$ and not to the function
$x \to |\cos(\pi x)|$. This later function is plotted in Figure
[[fig:writer-example2-2]] since ~output_list~ is set with the updated value
of ~fh~.
#+NAME: writer-example2-img
#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/writer-example2.png")
reset
set grid
set border
set key
set xtics
set ytics
set square
set terminal png truecolor enhanced size 640,640
plot '<(sed "" $PUGS_SOURCE_DIR/doc/writer-example2.gnu)' u 1:2 lw 2 t "f" w l, '<(sed "" $PUGS_SOURCE_DIR/doc/writer-example2.gnu)' u 1:3 lw 2 t "g" w l, '<(sed "" $PUGS_SOURCE_DIR/doc/writer-example2.gnu)' u 1:4 lw 2 t "sin(pi*f)" w l
#+END_SRC
#+CAPTION: Illustration of the life time of variables. The output for ~fh~ corresponds to its value when ~output_list~ is created: the interpolation of $x \to x$.
#+NAME: fig:writer-example2
#+ATTR_LATEX: :width 0.38\textwidth
#+ATTR_HTML: :width 300px;
#+RESULTS: writer-example2-img
#+NAME: writer-example2-2-img
#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/writer-example2-2.png")
reset
set grid
set border
set key
set xtics
set ytics
set square
set terminal png truecolor enhanced size 640,640
plot '<(sed "" $PUGS_SOURCE_DIR/doc/writer-example2-2.gnu)' u 1:2 lw 2 t "f" w l, '<(sed "" $PUGS_SOURCE_DIR/doc/writer-example2-2.gnu)' u 1:3 lw 2 t "g" w l
#+END_SRC
#+CAPTION: Illustration of the life time of variables. ~output_list~ is updated after ~fh~ has been updated.
#+NAME: fig:writer-example2-2
#+ATTR_LATEX: :width 0.38\textwidth
#+ATTR_HTML: :width 300px;
#+RESULTS: writer-example2-2-img
***** ~name_output: item_value*string -> output~
This function gives a name to an ~item_value~.
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import writer;
let m:mesh, m = cartesianMesh([0], [1], 100);
let x2:R^1 -> R, x -> x[0]*x[0];
let cell_values:item_value, cell_values = interpolate(m, node(), x2);
write(gnuplot_writer("writer-item-value-example1"), m, name_output(cell_values, "x^2"));
#+END_SRC
#+NAME: writer-item-value-example1-img
#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/writer-item-value-example1.png")
reset
set grid
set border
unset key
set xtics
set ytics
set square
set terminal png truecolor enhanced size 640,640
plot '<(sed "" $PUGS_SOURCE_DIR/doc/writer-item-value-example1.gnu)' lw 2 w l
#+END_SRC
#+CAPTION: Simple illustration of the output.
#+NAME: fig:writer-item-value-example1
#+ATTR_LATEX: :width 0.38\textwidth
#+ATTR_HTML: :width 300px;
#+RESULTS: writer-item-value-example1-img
#+BEGIN_warning
One observes that we supplied a mesh to the writer function. The
reason for that is that ~item_value~ refers to a connectivity. It requires
a ~mesh~ with the same connectivity to be written.
However, if a discrete function (of type ~Vh~) build on the same
connectivity is provided, there is no need to specify the ~mesh~.
#+END_warning
Let us give an example.
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
import writer;
import math;
let m:mesh, m = cartesianMesh([0], [1], 100);
let f:R^1 -> R, x -> 2*x[0]*x[0];
let g:R^1 -> R, x -> pow(sin(5*x[0]*x[0]),2);
let fh:Vh, fh = interpolate(m, P0(), f);
let cell_value:item_value, cell_value = interpolate(m, cell(), g);
let output_list:(output),
output_list = (name_output(fh, "f"),
name_output(cell_value, "g"));
write(gnuplot_1d_writer("writer-item-value-example2"), output_list);
#+END_SRC
Running this code produces the gnuplot file displayed in Figure
[[fig:writer-item-value-example2]].
#+NAME: writer-item-value-example2-img
#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/writer-item-value-example2.png")
reset
set grid
set border
set key
set xtics
set ytics
set square
set terminal png truecolor enhanced size 640,640
plot '<(sed "" $PUGS_SOURCE_DIR/doc/writer-item-value-example2.gnu)' u 1:2 lw 2 t "f" w l, '<(sed "" $PUGS_SOURCE_DIR/doc/writer-item-value-example2.gnu)' u 1:3 lw 2 t "g" w l
#+END_SRC
#+CAPTION: The ~item_value~ is plotted on the same mesh as ~fh~.
#+NAME: fig:writer-item-value-example2
#+ATTR_LATEX: :width 0.38\textwidth
#+ATTR_HTML: :width 300px;
#+RESULTS: writer-item-value-example2-img
***** ~gnuplot~ writers <<gnuplot-writers>>
There are three ~gnuplot~ writers. One is dedicated to output of
dimension 1 results (~gnuplot_1d_writer~), the second allows post
processing in dimension 1 and 2 (~gnuplot_writer~) and the third can be
used to write arbitrary data in the ~gnuplot~ format, whichever the
space dimension is (~gnuplot_raw_writer~). The last one is useful to
create scatter plots for instance.
All of these writers can be used for single output or time series
outputs. In the case of single output, the filename is completed by
adding the extension ~.gnu~, in the case of time series, the filename is
extending by adding ~.abcd.gnu~, where ~abcd~ is the number in the output
series.
#+BEGIN_note
The ~gnuplot~ writers are implemented in parallel.
The ~gnuplot~ post processing of produced files does not depend on the
number of processors (as soon as the saved data is also the same,
which is ensured by ~pugs~ for explicit methods).
#+END_note
For an obvious practical reason, each ~gnuplot~ file starts with a
preamble that indicates the running information.
This information contains, the production date of the file, the ~pugs~
version (including the hash key of the last commit) and a summary of
the output data to ease its use. For ~gnuplot~ files, this is the
location where the provided name (to the discrete function) is used.
In the case of time series the "physical" time of output is also
stored as a comment before the variable list.
Here is an example of preamble of a produced ~gnuplot~ file.
#+NAME: head-gnuplot-txt
#+BEGIN_SRC shell :exports results :results output
head -n 8 writer-example1.gnu
#+END_SRC
#+RESULTS: head-gnuplot-txt
****** ~gnuplot_1d_writer~ functions
This writer family makes sense only in 1d.
#+BEGIN_note
In parallel, as soon as the saved data themselves are the same, the
~gnuplot_1d_writer~ generates *exactly* the same files (whichever the
number of processors) since the coordinates of the post processed data
are sorted according to growing abscissa.
#+END_note
******* ~gnuplot_1d_writer: string -> writer~
Defines a ~gnuplot~ writer for 1d data. This is a single file writer in
the sense that it does not take care of time series.
Let us consider an example.
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
import writer;
import math;
let m:mesh, m = cartesianMesh([0], [1], 5);
let identity:R^1 -> R, x -> x[0];
let xh:Vh, xh = interpolate(m, P0(), identity);
let gp:writer, gp = gnuplot_1d_writer("gp_1d_example");
write(gp, (name_output(xh*xh, "x^2"), name_output(3*xh+2, "3x+2")));
#+END_SRC
The produced file (~"gp_1d_example.gnu"~) contains the following
#+NAME: gp-1d-example-txt
#+BEGIN_SRC shell :exports results :results output
cat gp_1d_example.gnu
#+END_SRC
#+RESULTS: gp-1d-example-txt
As one can see the output for $\mathbb{P}_0$ data consists in saving
the values at cell centers in the order of growing abscissa. This
allows to plot "smoother" curves than the ~gnuplot_writer~ (see
[[gnuplot-writer-function]]).
A typical use of this writer is the following.
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
import writer;
import math;
let m:mesh, m = cartesianMesh([0], [2], 20);
let pi:R, pi = acos(-1);
let sin_pi_x:R^1 -> R, x -> sin(pi*x[0]);
let fh:Vh, fh = interpolate(m, P0(), sin_pi_x);
write(gnuplot_1d_writer("gp_1d_sin"), name_output(fh, "f"));
#+END_SRC
#+NAME: writer-gp-1d-sin-img
#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/gp-1d-sin.png")
reset
set grid
set border
unset key
set xtics
set ytics
set square
set terminal png truecolor enhanced size 628,400
plot '<(sed "" $PUGS_SOURCE_DIR/doc/gp_1d_sin.gnu)' lw 2 w lp
#+END_SRC
#+CAPTION: Example of produced gnuplot results from the ~gnuplot_1d_writer~. Since data are stored by growing abscissa, one can use the ~w lp~ plotting option in ~gnuplot~.
#+NAME: fig:writer-gp-1d-sin
#+ATTR_LATEX: :width 0.38\textwidth
#+ATTR_HTML: :width 300px;
#+RESULTS: writer-gp-1d-sin-img
******* ~gnuplot_1d_writer: string*R -> writer~ <<gnuplot-1d-series>>
This writer differs from the previous one by handling output
series. The real value argument defines the period between two
outputs. It can be viewed as helper to outputs.
Let us give an example to fix ideas.
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
import writer;
import math;
let m:mesh, m = cartesianMesh([0], [1], 20);
let pi:R, pi = acos(-1);
let t:R, t = 0;
let tmax:R, tmax = 1;
let dt:R, dt = 3.7e-2;
let period:R, period = 0.15; // sets the period output
let gp:writer, gp = gnuplot_1d_writer("gp_1d_exp_sin", period);
let f: R^1 -> R, x -> exp(-t)*sin(2*pi*x[0]);
let fh:Vh, fh = interpolate(m, P0(), f);
write(gp, name_output(fh, "f"), 0);
do {
if (dt > tmax-t) {
dt = tmax-t;
}
t +=dt;
write(gp, name_output(fh, "f"), t);
} while(t<tmax);
#+END_SRC
Running this example produces the following files
#+NAME: ls-produced-gp-1d-series
#+BEGIN_SRC shell :exports results :results output
ls gp_1d_exp_sin.*.gnu
#+END_SRC
#+RESULTS: ls-produced-gp-1d-series
Each of these files contains the numerical solution at following
saving times:
#+NAME: times-in-gp-1d-series
#+BEGIN_SRC shell :exports results :results output
grep -n "# time = " gp_1d_exp_sin.*.gnu | cut -d '=' -f 2
#+END_SRC
#+RESULTS: times-in-gp-1d-series
****** ~gnuplot_writer~ functions <<gnuplot-writer-function>>
These writers differ from the previous ones since it draws the cells
and affects the cell value to its nodes. This produces larger files
but allows for 2d representations. Also, if the saved data is exactly
the same in parallel, the order of the cells is generally different
since they are written processor by processor.
Additionally, this writer allows to write 2d meshes, see paragraph
[[write-mesh]].
******* ~gnuplot_writer: string -> writer~
Here is an illustrating example in dimension 1, see the result on
Figure [[fig:writer-gp-sin]].
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
import writer;
import math;
let m:mesh, m = cartesianMesh([0], [2], 20);
let pi:R, pi = acos(-1);
let sin_pi_x:R^1 -> R, x -> sin(pi*x[0]);
let fh:Vh, fh = interpolate(m, P0(), sin_pi_x);
write(gnuplot_writer("gp_sin"), name_output(fh, "f"));
#+END_SRC
#+NAME: writer-gp-sin-img
#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/gp-sin.png")
reset
set grid
set border
unset key
set xtics
set ytics
set square
set terminal png truecolor enhanced size 628,400
plot '<(sed "" $PUGS_SOURCE_DIR/doc/gp_sin.gnu)' lw 2 w lp
#+END_SRC
#+CAPTION: Example of produced gnuplot results from the ~gnuplot_writer~. One can compare the produced result to the one of the ~gnuplot_1d_writer~ given in Figure [[fig:writer-gp-1d-sin]]
#+NAME: fig:writer-gp-sin
#+ATTR_LATEX: :width 0.38\textwidth
#+ATTR_HTML: :width 300px;
#+RESULTS: writer-gp-sin-img
Let use give a 2d example.
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
import writer;
import math;
let m:mesh, m = cartesianMesh(-[1,1], [1,1], (40,40));
let pi:R, pi = acos(-1);
let f:R^2 -> R, x -> cos(pi*x[0])*sin(pi*x[1]);
write(gnuplot_writer("gp_2d_cos_sin"),
name_output(interpolate(m,P0(), f), "f"));
#+END_SRC
The gnuplot result is displayed on Figure [[fig:writer-gp-2d-cos-sin]].
#+NAME: writer-gp-2d-cos-sin-img
#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/gp-2d-cos-sin.png")
reset
set grid
set border
unset key
set xtics
set ytics
set square
unset colorbox
set palette rgb 33,13,10
set terminal png truecolor enhanced size 300,300
plot '<(sed "" $PUGS_SOURCE_DIR/doc/gp_2d_cos_sin.gnu)' w filledcurves closed palette, '<(sed "" $PUGS_SOURCE_DIR/doc/gp_2d_cos_sin.gnu)' lc rgb 0 w l
#+END_SRC
#+CAPTION: Example of 2d plot from the ~gnuplot_writer~
#+NAME: fig:writer-gp-2d-cos-sin
#+ATTR_LATEX: :width 0.38\textwidth
#+ATTR_HTML: :width 300px;
#+RESULTS: writer-gp-2d-cos-sin-img
******* ~gnuplot_writer: string*R -> writer~
This is the time series function in the case of the ~gnuplot_writer~. It
behaves similarly to [[gnuplot-1d-series]].
****** ~gnuplot_raw_writer~ functions <<gnuplot-raw-writer-function>>
These writers are used to store raw data to the ~gnuplot~ format. Doing
so, it is the user responsibility to store coordinates in the file if
needed. Therefore it can be used in any space dimension.
******* ~gnuplot_raw_writer: string -> writer~
Here is an illustrating example in dimension 3, see the result on
Figure [[fig:writer-raw-sin]].
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
import writer;
import math;
let m:mesh, m = cartesianMesh([0,0,0], [2,2,2], (10,10,10));
let pi:R, pi = acos(-1);
let sin_pi_x:R^3 -> R, x -> sin(pi*dot(x,x));
let fh:Vh, fh = interpolate(m, P0(), sin_pi_x);
let r:Vh, r = sqrt(dot(cell_center(m),cell_center(m)));
write(gnuplot_raw_writer("raw_sin"), (name_output(r,"r"), name_output(fh, "f")));
#+END_SRC
#+NAME: writer-raw-sin-img
#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/raw-sin.png")
reset
set grid
set border
unset key
set xtics
set ytics
set square
set terminal png truecolor enhanced size 628,400
plot '<(sed "" $PUGS_SOURCE_DIR/doc/raw_sin.gnu)' lw 2 w lp
#+END_SRC
#+CAPTION: Example of produced gnuplot results from the ~gnuplot_raw_writer~.
#+NAME: fig:writer-raw-sin
#+ATTR_LATEX: :width 0.38\textwidth
#+ATTR_HTML: :width 300px;
#+RESULTS: writer-raw-sin-img
******* ~gnuplot_raw_writer: string*R -> writer~
This is the time series function in the case of the ~gnuplot_raw_writer~. It
behaves similarly to [[gnuplot-1d-series]].
***** ~vtk~ writers
For more complex post processing (including 3d), ~pugs~ can generate ~vtk~
outputs.
The used format is one file in the ~vtu~ format for each parallel domain
(and each time). The output is produced using binary data for
performance reasons. For each time step a ~pvtu~ file is generated to
handle parallelism. And for a complete time series, a ~pvd~ file is
produced. This is the file that should be loaded.
Observe that each of these files (~vtu~, ~pvtu~ and ~pvd~) contains a
comment that stores the creation date and the version of ~pugs~ that was
run.
The use is exactly the same as for ~gnuplot~ writers so we do not
provide full examples.
~vtk~ writers are compatible with the ~write_mesh~ function, see paragraph
[[write-mesh]].
****** ~vtk_writer: string -> writer~
One should use this writer for single output (no time series). The
produced ~pvd~ file is built by adding ~.pvd~ to the provided ~string~.
****** ~vtk_writer: string*R -> writer~
This function follows the same rule. One just specifies the output
period. The generated ~pvd~ file is built similarly, one adds ~.pvd~ to
the provided ~string~.
***** ~write~, ~force_write~ and ~write_mesh~ functions
Once a mesh writer has been defined, these functions are called to
effectively generate the post processing files.
****** ~write: writer*(output) -> void~
As parameters, it takes a ~writer~ (that is not a time series writer)
and a list of ~output~ (which are named discrete functions that are
defined on the *same* mesh).
We have already shown a lot of examples of use of the ~write~
function. So let us focus on improper calls.
Trying to use a time series ~writer~
#+NAME: cannot-use-time-series-writer
#+BEGIN_SRC pugs-error :exports both :results output
import mesh;
import scheme;
import writer;
import math;
let m:mesh, m = cartesianMesh([0], [2], 20);
let pi:R, pi = acos(-1);
let sin_pi_x:R^1 -> R, x -> sin(pi*x[0]);
let fh:Vh, fh = interpolate(m, P0(), sin_pi_x);
let vtk:writer, vtk = vtk_writer("vtk_time_series", 0);
write(vtk, name_output(fh, "f"));
#+END_SRC
produces the following runtime error
#+results: cannot-use-time-series-writer
Trying to use variables defined on different meshes
#+NAME: cannot-use-diffrent-meshes-in-writer
#+BEGIN_SRC pugs-error :exports both :results output
import mesh;
import scheme;
import writer;
import math;
let m0:mesh, m0 = cartesianMesh([0], [2], 20);
let m1:mesh, m1 = cartesianMesh([0], [2], 20);
let pi:R, pi = acos(-1);
let sin_pi_x:R^1 -> R, x -> sin(pi*x[0]);
let fh:Vh, fh = interpolate(m0, P0(), sin_pi_x);
let gh:Vh, gh = interpolate(m1, P0(), sin_pi_x);
let vtk:writer, vtk = vtk_writer("vtk_different_meshes");
write(vtk, (name_output(fh, "f"), name_output(gh, "g")));
#+END_SRC
gives the runtime error
#+results: cannot-use-diffrent-meshes-in-writer
****** ~write: writer*mesh*(output) -> void~
This variation permits to specify a ~mesh~. Its main utility is to allow
to write data in the case where ~output~ list contains only ~item_value~
data which are not related to a mesh but to a *connectivity*. One checks
that the connectivity of the ~output~ data are the same as the
connectivity of the ~mesh~.
****** ~write: writer*(output)*R -> void~
This ~write~ function is used to save time series data. The real
parameter (of type ~R~) is the current time.
One cannot use a ~writer~ that does not support time series.
Trying to use a non time series ~writer~
#+NAME: cannot-use-non-time-series-writer
#+BEGIN_SRC pugs-error :exports both :results output
import mesh;
import scheme;
import writer;
import math;
let m:mesh, m = cartesianMesh([0], [2], 20);
let pi:R, pi = acos(-1);
let sin_pi_x:R^1 -> R, x -> sin(pi*x[0]);
let fh:Vh, fh = interpolate(m, P0(), sin_pi_x);
let vtk:writer, vtk = vtk_writer("vtk_time_series");
write(vtk, name_output(fh, "f"), 1.2);
#+END_SRC
gives the runtime error
#+results: cannot-use-non-time-series-writer
****** ~write: writer*mesh*(output)*R -> void~
This is a variation of the previous function where a mesh is
specified. The function checks the validity of the ~output~ list:
functions ~Vh~ must defined on the given ~mesh~ and ~item_value~ data mush
share the same connectivity with the ~mesh~.
****** ~force_write: writer*(output)*R -> void~
One probably noticed that using the ~write~ function with a time series
~writer~, last time of the calculation may not be written (see section
[[gnuplot-1d-series]]). The ~force_write~ function does not check that the
saving time has been reached. It only checks that the current time has
not already been saved.
Let us improve slightly the example given in section
[[gnuplot-1d-series]].
#+BEGIN_SRC pugs :exports both :results none
import mesh;
import scheme;
import writer;
import math;
let m:mesh, m = cartesianMesh([0], [1], 20);
let pi:R, pi = acos(-1);
let t:R, t = 0;
let tmax:R, tmax = 1;
let dt:R, dt = 3.7e-2;
let period:R, period = 0.15; // sets the period output
let gp:writer, gp = gnuplot_1d_writer("gp_1d_exp_sin_force", period);
let f: R^1 -> R, x -> exp(-t)*sin(2*pi*x[0]);
let fh:Vh, fh = interpolate(m, P0(), f);
write(gp, name_output(fh, "f"), 0);
do {
if (dt > tmax-t) {
dt = tmax-t;
}
t +=dt;
fh = interpolate(m, P0(), f);
if (t<tmax) {
write(gp, name_output(fh, "f"), t);
} else {
force_write(gp, name_output(fh, "f"), t);
}
} while(t<tmax);
#+END_SRC
Running this example produces the following files
#+NAME: ls-produced-gp-1d-series-force
#+BEGIN_SRC shell :exports results :results output
ls gp_1d_exp_sin_force.*.gnu
#+END_SRC
#+RESULTS: ls-produced-gp-1d-series-force
One can see the additional file.
Each of these files contains the numerical solution at following
saving times:
#+NAME: times-in-gp-1d-series-force
#+BEGIN_SRC shell :exports results :results output
grep -n "# time = " gp_1d_exp_sin_force.*.gnu | cut -d '=' -f 2
#+END_SRC
#+RESULTS: times-in-gp-1d-series-force
The last post processing time is now 1.
****** ~force_write: writer*mesh*(output)*R -> void~
This variation permits to specify a ~mesh~. Compatibility of output data
with the ~mesh~ is checked.
****** ~write_mesh: writer*mesh -> void~ <<write-mesh>>.
This function just saves a ~mesh~ using a ~writer~. The ~gnuplot_1d_writer~
cannot write mesh. And the ~writer~ argument must be a non time series
~writer~.
*** The ~socket~ module
This module is more dedicated to developers that intend to couple ~pugs~
to other codes using sockets.
We describe here language utilities that can be useful, but effective
coupling generally requires ~C++~ coding.
Here is the list of the functions and types that the module provides.
#+NAME: get-module-info-socket
#+BEGIN_SRC pugs :exports both :results output
cout << getModuleInfo("socket") << "\n";
#+END_SRC
#+RESULTS: get-module-info-socket
**** ~soket~ provided types
***** ~socket~
This is the ~socket~ descriptor type.
This module adds the following binary operator
| ~ostream <<~ allowed expression type |
|------------------------------------|
| ~socket~ |
For instance running the following code
#+NAME: print-socket
#+BEGIN_SRC pugs :exports both :results output
import socket;
let s:socket, s = createSocketServer(0);
cout << s << "\n";
#+END_SRC
gives
#+RESULTS: print-socket
The output uses the classical format: ~machine_name:port_number~.
**** ~socket~ provided functions
***** socket management functions
There are three socket management functions.
****** ~createSocketServer: N -> socket~
This function creates a socket server on the current machine. The
integer parameter (of type ~N~) is the port number. If one uses the
magical value $0$ a free port is assigned (this is the standard
behavior for server socket creation). The return value is the ~socket~
itself.
****** ~connectSocketServer: string*N -> socket~
Connects ~pugs~ to the server described by a name and a port number. It
returns a socket to the server.
****** ~acceptSocketClient: socket -> socket~
Accepts connection from a client. The argument is the server socket
and the return socket to the client.
***** writing and reading basic types through sockets
****** Writing functions
These polymorcphic overload send a basic value through the given ~socket~
- ~write: socket*B -> void~
- ~write: socket*N -> void~
- ~write: socket*Z -> void~
- ~write: socket*R -> void~
- ~write: socket*Rˆ1 -> void~
- ~write: socket*Rˆ2 -> void~
- ~write: socket*Rˆ3 -> void~
- ~write: socket*Rˆ1x1 -> void~
- ~write: socket*Rˆ2x2 -> void~
- ~write: socket*Rˆ3x3 -> void~
- ~write: socket*string -> void~
****** Reading functions
These are the reading counterparts of the previous functions. The
argument is the ~socket~ on which values are read.
- ~read_B: socket -> B~
- ~read_N: socket -> N~
- ~read_Z: socket -> Z~
- ~read_R: socket -> R~
- ~read_R1: socket -> Rˆ1~
- ~read_R2: socket -> Rˆ2~
- ~read_R3: socket -> Rˆ3~
- ~read_R1x1: socket -> Rˆ1x1~
- ~read_R2x2: socket -> Rˆ2x2~
- ~read_R3x3: socket -> Rˆ3x3~
- ~read_string: socket -> string~
[fn:pugs-def] ~pugs~: Parallel Unstructured Grid Solvers
[fn:MPI-def] ~MPI~: Message Passing Interface
[fn:DSL-def] ~DSL~: Domain Specific Language