ElmerSolver Manual Ftp://nic.funet.fi/pub/sci/physics/elmer/doc/Elmer Solver Elmer

User Manual: ftp://nic.funet.fi/pub/sci/physics/elmer/doc/ElmerSolverManual

Open the PDF directly: View PDF PDF.
Page Count: 158 [warning: Documents this large are best viewed by clicking the View PDF Link!]

ElmerSolver Manual
Juha Ruokolainen, Mika Malinen, Peter Råback,
Thomas Zwinger, Antti Pursula and Mikko Byckling
CSC – IT Center for Science
May 22, 2017
ElmerSolver Manual
About this document
The ElmerSolver Manual is part of the documentation of Elmer finite element software. ElmerSolver Man-
ual describes the Elmer Solver options common for all specific equation solvers. The different equations
solver options are described separately in Elmer Models Manual. The ElmerSolver Manual is best used as a
reference manual rather than a concise introduction to the matter.
The present manual corresponds to Elmer software version 8.3 for Windows NT and Unix platforms.
The latest documentation and program versions of Elmer are available (or links are provided) at
http://www.csc.fi/elmer.
Copyright information
The original copyright of this document belongs to CSC – IT Center for Science, Finland, 1995–2015. This
document is licensed under the Creative Commons Attribution-No Derivative Works 3.0 License. To view a
copy of this license, visit http://creativecommons.org/licenses/by-nd/3.0/.
Elmer program is free software; you can redistribute it and/or modify it under the terms of the GNU
General Public License as published by the Free Software Foundation; either version 2 of the License, or (at
your option) any later version. Elmer software is distributed in the hope that it will be useful, but without
any warranty. See the GNU General Public License for more details.
Elmer includes a number of libraries licensed also under free licensing schemes compatible with the
GPL license. For their details see the copyright notices in the source files.
All information and specifications given in this document have been carefully prepared by the best ef-
forts of CSC, and are believed to be true and accurate as of time writing. CSC assumes no responsibility or
liability on any errors or inaccuracies in Elmer software or documentation. CSC reserves the right to modify
Elmer software and documentation without notice.
2
Contents
Table of Contents 3
1 Solving a multiphysics problem with the solver of Elmer: Fundamentals 6
1.1 Basic concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2 Handling interactions of multiple physical phenomena . . . . . . . . . . . . . . . . . . . . 7
1.3 The key abilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2 Structure of the Solver Input File 10
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2 The sections of solver input file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3 Keyword syntax . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.4 Running several sequences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3 Restart from existing solutions 18
3.1 Restart file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.2 Initialization of dependent variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4 Solution methods for linear systems 21
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.2 Direct methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.3 Preconditioned Krylov methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
4.4 Multilevel methods ....................................... 23
4.5 Preconditioning via inner iterations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.6 Keywords related to linear system solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.7 Implementation issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
5 Nonlinear System Options 33
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
5.2 Keywords related to solution of nonlinear systems . . . . . . . . . . . . . . . . . . . . . . . 33
6 Integration of time-dependent systems 36
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
6.2 Time discretization strategies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
6.3 Keywords related to time discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
6.4 On the treatment of time derivatives in Elmer Solver code . . . . . . . . . . . . . . . . . . . 39
6.5 Predictor Corrector Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
7 Solving eigenvalue problems 41
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
7.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
7.3 Keywords related to eigenvalue problems . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
7.4 Constructing matrices M and D in Solver code . . . . . . . . . . . . . . . . . . . . . . . . . 44
3
CONTENTS
4
8 Dirichlet boundary conditions and limiters 45
8.1 Dirichlet constraints ....................................... 45
8.2 Soft limiters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
9 Enforcing shared conditions between boundaries 48
9.1 Periodic conditions ....................................... 48
9.2 Mortar conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
10 Linear systems with constraints 54
10.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
10.2 Keywords related to treatment of linear constraints . . . . . . . . . . . . . . . . . . . . . . 55
11 Solving constraint modes 56
11.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
11.2 Keywords related to constraint modes analysis . . . . . . . . . . . . . . . . . . . . . . . . . 56
12 Working with nodal loads 57
12.1 Setting nodal loads ....................................... 57
12.2 Computing nodal loads . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
12.3 Computing nodal weights . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
12.4 Energy norm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
13 Generic solver utilities 59
13.1 Solver activation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
13.2 Variable names . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
13.3 Exported variables ....................................... 60
13.4 Active and passive elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
13.5 Timing the solvers ....................................... 61
13.6 Consistency testing of solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
14 Meshing Utilities 62
14.1 Coordinate transformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
14.2 Mesh multiplication ....................................... 63
14.3 Mesh extrusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
14.4 Discontinuous interfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
15 Adaptive Solution 66
15.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
15.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
15.3 Keywords related to the adaptive solution . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
15.4 Implementing own error estimators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
16 Parallel runs 71
16.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
16.2 Preprocessing of Parallel Runs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
16.3 Parallel Computations in Elmer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
16.4 Post-processing of Parallel Runs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
17 Compilation and Linking 82
17.1 Compiling with CMake . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
17.2 Compiling a user-defined subroutine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
CONTENTS
5
18 Basic Programming 84
18.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
18.2 Basic Elmer Functions and Structures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
18.3 Writing a User Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
18.4 Writing a Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
18.5 Compilation and Linking of User Defined Routines/Functions . . . . . . . . . . . . . . . . 116
A Format of mesh files 117
A.1 The format of header file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
A.2 The format of node file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
A.3 The format of element file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
A.4 The format of boundary element file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
A.5 The format of shared nodes file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
A.6 Exceptions on parallel mesh format . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
B Format of result output files 120
B.1 Format versions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
B.2 General structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
B.3 The positions file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
C Format of ElmerPost Input File 124
D Basic element types 126
E Higher-order finite elements 129
E.1 Higher-order elements in Elmer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
E.2 ElmerSolver services for higher-order elements . . . . . . . . . . . . . . . . . . . . . . . . 132
E.3 Basis functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
F Face and edge elements 144
F.1 The construction of face element interpolation . . . . . . . . . . . . . . . . . . . . . . . . . 144
F.2 The construction of edge element interpolation . . . . . . . . . . . . . . . . . . . . . . . . 148
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
G Advanced finite element definitions 154
G.1 Specification of the degrees of freedom . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154
Index 155
Chapter 1
Solving a multiphysics problem with the
solver of Elmer: Fundamentals
Elmer software has been developed multiphysics simulations in mind. Thus, in addition to offering ways
to produce computational solutions to single-physics models (the available collection of which is described
in the Elmer Models Manual), Elmer provides procedures for creating computational models which describe
interactions of multiple physical phenomena. Our intention here is to give an overview how this functionality
is built into the solver of Elmer associated with the independent program executable ElmerSolver.
1.1 Basic concepts
The models handled by Elmer may generally be stationary or evolutionary, with nonlinearities possible
in both the cases. Starting from the weak formulation of the problem, finite element approximation and
advancing in time with implicit methods are typically applied in order to obtain the computational version
of the model. In the simplest case of single-physics models we are then lead to solving equations
F(u) = 0,(1.1)
where urepresents either the vector of coefficients in the finite element expansion of the stationary solution
or the coefficient vector to describe the evolutionary finite element solution at a given time t=tk. Thus, in
the case of evolution, the problems of the type (1.1) are solved repeatedly when advancing in time.
For linear models the problem (1.1) reduces to solving a linear system via defining
F(u) = b−Ku
where the coefficient matrix Kis referred to as the stiffness matrix and bcorresponds to the right-hand side
vector in the linear system. Otherwise Fis a nonlinear mapping and an iteration is needed to handle the
solution of the problem (1.1). Available nonlinear iteration methods generally depend on the model, as the
definition of the linearization strategy is a part of the computational description of each physical model.
We note that many single-physics models offer the possibility of using the Newton iteration where the
current nonlinear iterate u(m)to approximate uis updated at each iteration step as
DF (u(m))[δ(m)] = −F(u(m)),
u(m+1) =u(m)+δ(m),(1.2)
where DF (u(m))is the derivative of Fat u(m). Thus, performing the nonlinear solution update again
entails the solution of the linear system at each iteration step. As an alternate to the Newton method, lin-
earization strategies based on lagged-value approximations are also often available. In addition, relaxation
is conventionally offered as a way to enable convergence in cases where the basic nonlinear iteration fails
CSC – IT Center for Science
1. Solving a multiphysics problem with the solver of Elmer: Fundamentals 7
to produce convergence. Given the current nonlinear iterate u(m)and a computed correction δu(m)to the
approximation, the new nonlinear iterate is then defined by
u(m+1) =u(m)+λ(m)δ(m),
where λ(m)is an adjustable parameter referred to as the relaxation parameter.
1.2 Handling interactions of multiple physical phenomena
Having considered the basic concepts in the context of single-physics models, we now proceed to describe
how the modularity employed in the design of Elmer allows us to create models which represent interactions
of multiple physical phenomena. To this end, we assume that the complete model describes the interaction
of Nconstituent models, the computational versions of which are primarily associated with the coefficient
vectors ui,i= 1,2,...,N. As before, the coefficients contained in uiare associated with the finite element
expansion of either the stationary solution or the evolutionary solution at a time level t=tk.
The fully discrete version of the coupled model leads to handling a problem of the form
F1(u1, u2,...,uN) = 0,
F2(u1, u2,...,uN) = 0,
··· (1.3)
FN(u1, u2,...,uN) = 0.
If all the constituent models are linear, the problem (1.3) corresponds to solving the linear system where the
coefficient matrix is a N×Nblock matrix. Otherwise (1.3) describes a nonlinear problem. Although the
solution of (1.3) could in principle be done in the same way as explained in the context of single-physics
models in Section 1.1, i.e. by performing either a coupled linear solve or Newton iteration, the coupled
problems have usually been handled differently in order to enable the reuse of solvers for single-physics
models and the easy extendability of the code to handle new applications.
To this end, the nonlinear Gauss-Seidel iteration is usually applied, so that the coupling of the models is
resolved via generating new coupled system iterates u(j)= (u(j)
1, u(j)
2,...,u(j)
N)as
F1(u(j)
1, u(j−1)
2, u(j−1)
3,...,u(j−1)
N) = 0,
F2(u(j)
1, u(j)
2, u(j−1)
3,...,u(j−1)
N) = 0,
··· (1.4)
FN(u(j)
1, u(j)
2,...,u(j)
N) = 0.
It is noted that the kth discrete model description in (1.4) depends implicitly only on the coupled system
iterate to its primary variable uk, while the dependencies on the other constituent model variables are treated
explicitly. This brings us to solving a nonlinear single-field problem
F(u(j)
k) = Fk(v1,...,vk−1, u(j)
k, vk+1,...,vN) = 0,with all vjgiven,(1.5)
which is handled by using the methods already described in Section 1.1. We also note that if all the con-
stituent models are linear the nonlinear Gauss-Seidel iteration (1.4) reduces to the block Gauss-Seidel itera-
tion for linear systems. Relaxation may again be applied as an attempt to improve the convergence behaviour
of the basic iteration (1.4).
It is good to pause here to stress that the main advantage of the adopted nonlinear Gauss-Seidel scheme is
its support for the modular software design. Also, it brings us to handling coupled problems via solving linear
systems which are smaller than those which would result from treating all constraints in (1.3) simultaneously.
Despite these merits, the suitability of the loosely coupled iteration (1.4) generally is case-dependent as
convergence problems may occur in cases where a strong physical coupling is involved. Such problems
are often best handled by methods which treat all the constituent models in (1.3) simultaneously. Certain
physical models available in Elmer indeed employ this alternate tightly coupled solution strategy. However,
CSC – IT Center for Science
1. Solving a multiphysics problem with the solver of Elmer: Fundamentals 8
these models have initially been developed independently, as common Elmer utilities for creating tightly
coupled iteration methods in a general manner are less developed.
To summarize, the following pseudo-code presentation describes the basic loosely coupled iteration
scheme employed by the solver of Elmer. This rough description may be helpful in summarizing what
needs to be controlled overall to create a working computational solution procedure for a coupled problem.
! The time integration loop
for k= 1 : M
Generate an initial guess u(0) = (u(0)
1, u(0)
2,...,u(0)
N)for the coupled system solution at t=tk
! The nonlinear Gauss-Seidel iteration
for j= 1,2,...
! Generate the next coupled system iterate u(j)by performing single-field updates
for i= 1 : N
Set vl=u(j)
lfor l= 1,2,...,i−1
Set vl=u(j−1)
lfor l=i+ 1 : N
Perform nonlinear iteration to solve Fi(v1,...,vi−1, u(j)
i, vi+1,...,vN) = 0
Apply relaxation to set u(j)
i:= u(j−1)
i+αi(u(j)
i−u(j−1)
i)
end
end
end
Here the descriptions of the termination criteria for the iterations have been omitted. It is also noted that,
obviously, the time integration loop is not needed in the case of a stationary problem. On the other hand,
in the case of stationary simulation it is possible to replace the time integration loop by a pseudo-version of
time stepping to enable performing multiple solutions for a range of model parameter values.
1.3 The key abilities
In the following, we give some additional information on the key abilities of the solver of Elmer to create
computational solution procedures.
1.3.1 Extendability by modular design
A module of the Elmer software which enables the creation of the discrete model description of the type (1.5)
and its solution with respect to the primary variable is generally called a solver. The solvers of Elmer are
truly modular in this manner and have a standard interface. Thus, each solver usually contains an implemen-
tation of the nonlinear iteration, instructions to assemble the corresponding linear systems from elementwise
contributions, and standard subroutine invocations to actually solve the linear systems assembled.
It follows that enabling an interaction with another field, designated by vjin (1.5), is simply a matter of
solver-level implementation. Therefore, interactions which have not been implemented yet can be enabled
by making modifications which are localized to the solvers. In addition, a completely new physical model
may be added by introducing a new solver which comprises a separate software module and which can be
developed independently with respect to the main program. As a result, a loosely coupled solution procedure
for a coupled problem based on the new physical model may again be achieved by making only solver-level
modifications.
1.3.2 Model-specific finite element approximation
In the most basic setting all constituent model variables ui,i= 1,...,N, of a coupled problem are approx-
imated by using the same piecewise polynomial basis functions defined over a single mesh. In addition to
this, the solver of Elmer offers a built-in functionality to perform a coupled problem simulation by using
solver-specific finite element meshes. The solver description is then augmented by the specification of the
CSC – IT Center for Science
1. Solving a multiphysics problem with the solver of Elmer: Fundamentals 9
independent mesh which the solver uses. To make this functional in connection with the solution of cou-
pled problems, Elmer has the capability of performing the solution data transfer, which is needed between
the solvers in the loosely coupled solution procedure, even when the meshes are non-matching. It must be
understood, however, that the loss of high-resolution details is unavoidable when the high-resolution field is
represented by using a coarser finite element mesh.
1.3.3 Approximation by various finite element formulations
Elmer has traditionally employed the Galerkin finite element approximation of weak formulation based on
the Lagrange interpolation basis functions. In this connection, the piecewise polynomial approximation of
degree 1≤p≤3is possible for 2-D models, while 3-D models may be discretized by using the elements
of degree 1≤p≤2. The isoparametric mapping to describe curved element shapes is also supported with
these traditional elements.
Discrete models based on more recent versions of the Galerkin finite element approximation are also
possible. As an alternate to using the standard Lagrange interpolation basis functions, the Galerkin approxi-
mation based on using high-degree piecewise polynomials can be employed. In this connection, the degree
of polynomial approximation can also be defined elementwise, with Elmer providing an in-built mechanism
to guarantee the continuity of any solution candidate, so that in effect the use of the hp-version of the finite
element method is enabled. However, generic ways to describe curved body surfaces accurately in connec-
tion with the high-degree finite elements have not been implemented yet which may limit the utility of these
elements.
The way to define the high-degree approximation is based on the idea that a background mesh for rep-
resenting the standard lowest-degree continuous finite element expansion is first provided so that a specific
element type definition in relation to elements present in the background mesh may then be given to enhance
the approximation. The same idea has been adapted to create other alternate finite element formulations.
For example, finite element formulations which enhance the approximation defined on the background mesh
by a subscale approximation spanned by elementwise bubble basis functions can be obtained in this way.
We note that this strategy is widely used in Elmer to stabilize otherwise unstable formulations. Another
example of the use of the user-supplied element definition relates to creating approximations based on the
discontinuous Galerkin method. As a final example we mention that enhancing the approximation on the
background mesh by either face or edge degrees of freedom and then omitting the original nodal degrees of
freedom is also possible. This leads to a suitable set of unknowns for creating discretizations based on the
Raviart-Thomas or edge element interpolation.
1.3.4 Parallel computing
A strength of the solver of Elmer is that it supports the use of parallel computing. This opportunity sig-
nificantly widens the range of problems which can be considered. Additional details on utilizing parallel
computers are found elsewhere in this manual.
1.3.5 Linear algebra abilities
This exposition should have made it clear that having the ability to solve large linear systems efficiently is a
central aspect of the simulation process with Elmer. As explained, in the basic setting a linear solve is needed
to obtain the solution update at each step of the nonlinear iteration. In practice linear solves are usually done
iteratively, revealing one unexposed iteration level in relation to the pseudo-code presentation in the end of
Section 1.2.
The solver of Elmer offers a large selection of strategies to construct linear solvers. The majority of
them are directly implemented into Elmer software, but interfaces to exploit external linear algebra libraries
are also available. Typically the most demanding aspect in the use of linear solvers relates to identifying
an effective preconditioning strategy for the problem at hand. Traditionally Elmer has employed generic
preconditioning strategies based on the fully algebraic approach, but recently alternate block preconditioning
strategies which typically try to utilize physics-based splittings have also been developed.
CSC – IT Center for Science
Chapter 2
Structure of the Solver Input File
2.1 Introduction
Solving partial differential equation (PDE) models with the solver of Elmer requires that a precise description
of the problem is given using the so-called solver input file, briefly referred to as the sif file. This file contains
user-prepared input data which specify the location of mesh files and control the selection of physical models,
material parameters, boundary conditions, initial conditions, stopping tolerances for iterative solvers, etc. In
this chapter, the general structure of the file is described. We explain how the input data is organized into
different sections and describe the general keyword syntax which is used in these sections to define the
values of various model parameters and to control the solution procedures.
In the case of simple problem setups the solver input file may be written automatically by the prepro-
cessor of Elmer software, so that knowing the solver input file format may be unnecessary. Creating a more
complicated setup, or using keywords introduced by the user, however, requires the knowledge of the file
format and keyword syntax.
In the following the general structure of the input file is first illustrated by using simple examples, without
trying to explain all possibilities in an exhaustive manner. We then describe the keyword syntax in more
detail, showing also how model parameters whose values depend on solution fields can be created. The later
chapters of this manual, and Elmer Models Manual, which focuses on describing the PDE models Elmer can
handle, provide more detailed material on specific issues. Elmer Tutorials also gives complete examples of
solver input files.
2.2 The sections of solver input file
The material of the solver input file is organized into different sections. Each section is generally started
with a row containing the name of the section, followed by a number of keyword commands, and ended with
a row containing the word End. The names for starting new sections are
•Header
•Simulation
•Constants
•Body n
•Material n
•Body Force n
•Equation n
CSC – IT Center for Science
2. Structure of the Solver Input File 11
•Solver n
•Boundary Condition n
•Initial Condition n
•Component n
Here nassociated with the section name represents an integer identifier needed for distinguishing between
sections of the same type. A basic keyword command included in a section is nothing more than a statement
which sets the value of a keyword with an equal sign.
In the following we describe how the sections are basically arranged without trying to explain all possi-
bilities in an exhaustive manner. The later chapters of this manual and Elmer Models Manual provide more
detailed material on specific issues. Elmer Tutorials also gives complete examples of solver input files.
Header section
The location of mesh files is usually given in the header section. Often this is also the only declaration given
in the header section. If the Elmer mesh files (see Appendix A) are located in the directory ./mymesh, the
header section may simply be
Header
Mesh DB "." "mymesh"
End
Note that separate equations can nevertheless be discretized using different meshes if the location of mesh
files is given in the solver section described below.
Simulation section
The simulation section is used for giving general information that is not specific to a particular PDE model
involved in the simulation. This information describes the coordinate system used, indicates whether the
problem is stationary or evolutionary, defines the file names for outputting, etc. Without trying to describe
many possibilities and the details of commands, we only give the following simple example:
Simulation
Coordinate System = "Cartesian 2D"
Coordinate Mapping(3) = 1 2 3
Coordinate Scaling = 0.001
Simulation Type = Steady State
Steady State Max Iterations = 1
Output Intervals(1) = 1
Post File = "case.ep"
Output File = "case.dat"
End
Constants section
The constants section is used for defining certain physical constants. For example the gravity vector and the
Stefan-Boltzmann constant may be defined using the commands
Constants
Gravity(4) = 0 -1 0 9.82
Stefan Boltzmann = 5.67e-08
End
If the constants are not actually needed in the simulation, this section can also be left empty.
CSC – IT Center for Science
2. Structure of the Solver Input File 12
Body, material, body force and initial condition sections
The Elmer mesh files contain information on how the computational region is divided into parts referred
to as bodies. A body section associates each body with an equation set, material properties, body forces,
and initial conditions by referring to definitions given in a specified equation section, material section, body
force section, and initial condition section. To manage to do this, the different sections of the same type
are distinguished by integer identifiers that are parts of the section names. Note that the integer in the body
section name is an identifier for the body itself.
For example, one may define
Body 1
Material = 1
Body Force = 1
Equation = 1
Initial Condition = 2
End
Material properties, body forces, an equation set, and initial conditions are then defined in the material
section
Material 1
...
End
the body force section
Body Force 1
...
End
the equation section
Equation 1
...
End
and the initial condition section
Initial Condition 2
...
End
What material properties and body forces need to be specified depends on the mathematical models involved
in the simulation, and the initial condition section used for giving initial values is only relevant in the so-
lution of evolutionary problems. We here omit the discussion of these very model-dependent issues; after
reading this introductory chapter the reader should be able to understand the related documentation given in
Elmer Models Manual, which focuses on describing the different mathematical models, while the contents
of equation section will be described next.
Equation and solver sections
Equation section provides us a way to associate each body with a set of solvers, where each solver is typically
associated with the ability to solve a certain physical model; cf. the definition of solver given in the beginning
of Section 1.3.1. That is, if the set defined consists of more than one solver, several physical phenomena may
be considered to occur simultaneously over the same region of space. The actual definitions of the solvers
are given in solver sections, the contents of an equation section being basically a list of integer identifiers
for finding the solver sections that define the solvers. The keyword commands given in the solver sections
then control the selection of physical models, linearization procedures of nonlinear models, the selection of
solution methods for resulting linear equations, convergence tolerances, etc.
For example, if only two solvers are needed, one may simply define a list of two solver identifiers
CSC – IT Center for Science
2. Structure of the Solver Input File 13
Equation 1
Active Solvers(2) = 1 2
End
Then the solver definitions are read from the solver sections
Solver 1
...
End
and
Solver 2
...
End
Finally, we give an example of solver definitions, without trying to explain the commands at this point:
Solver 1
Equation = "Poisson"
Variable = "Potential"
Variable DOFs = 1
Procedure = "Poisson" "PoissonSolver"
Linear System Solver = "Direct"
Steady State Convergence Tolerance = 1e-06
End
Boundary condition section
Boundary condition sections define the boundary conditions for the different models. The Elmer mesh files
contain information on how the boundaries of the bodies are divided into parts distinguished by their own
boundary numbers. The keyword Target Boundaries is used to list the boundary numbers that form
the domain for imposing the boundary condition. For example the declaration
Boundary Condition 1
Target Boundaries(2) = 1 2
...
End
means that the boundary condition definitions that follow concern the two parts having the boundary numbers
1 and 2.
Component section
Component is a physical entity that is not associated with the mesh. It could be a representation of a 0D
object such as electrical component, equation of state, chemical reaction, rigid object etc. The component
could be independent, or it can be target for model reduction in the higher dimensional parts.
The following is an example of a component which is used as a reduction target for two bodies.
Component 1
Name = "MyComponent"
Master Bodies(2) = 1 2
...
End
This section that has been added much later than the other sections and it is still quite little used.
CSC – IT Center for Science
2. Structure of the Solver Input File 14
Text outside sections
We finally note that some commands, such as comments started with the symbol ! and MATC expres-
sions described below, may also be placed outside section definitions. An exception of this type is also the
command
Check Keywords "Warn"
usually placed in the beginning of the input file. When this command is given, the solver outputs warning
messages if the input file contains keywords that are not listed in the file of known keywords. If new
keywords are introduced, misleading warning messages can be avoided by adding the new keywords to the
keyword file SOLVER.KEYWORDS, located in the directory of the shared library files of ElmerSolver. The
other options include ignore, abort, silent.
There is also the commands echo on and echo off that may be used to control the output of the
parser. This is mainly intended for debugging purposes. The default is off.
2.3 Keyword syntax
As already illustrated, a basic keyword command used in the solver input file is a statement which sets the
value of a solution parameter with the equal sign. Such a command in its full form also contains the data
type declaration; for example
Density = Real 1000.0
Valid data types generally are
•Real
•Integer
•Logical
•String
•File
If the keyword is listed in the keyword file SOLVER.KEYWORDS, the data type declaration may be omitted.
Therefore, in the case of our example, we may also define
Density = 1000.0
The value of a keyword may also be an array of elements of specified data type, with the array size
definition associated with the keyword. We recall our previous examples of the equation and boundary
condition sections, where we defined two lists of integers using the commands
Active Solvers(2) = 1 2
and
Target Boundaries(2) = 1 2
Two-dimensional arrays are also possible and may be defined as
My Parameter Array(3,3) = Real 1 2 3 \
4 5 6 \
7 8 9
CSC – IT Center for Science
2. Structure of the Solver Input File 15
Defining parameters depending on field variables
Most solver parameters may depend on time, or on the field variables defined in the current simulation
run. Such dependencies can generally be created by means of tabular data, MATC functions, or Fortran
functions. MATC has the benefit of being an interpreted language, making an additional compilation step
with a compiler unnecessary.
Simple interpolating functions can be created by means of tabular data. The following example defines
the parameter Density the value of which depends on the variable Temperature:
Density = Variable Temperature
Real
0 900
273 1000
300 1020
400 1000
End
This means that the value of Density is 900 when Temperature is 0, and the following lines are
interpreted similarly. Elmer then uses linear interpolation to approximate the parameter for argument values
not given in the table. If the value of the independent variable is outside the data set, the first or the last
interpolating function which can be created from the tabulated values is used to extrapolate the value of the
parameter.
If the field variable has several independent components, such as the components of displacement vector,
the independent components may be used as arguments in a function definition. For example, if a three-
component field variable is defined in a solver section using the commands
Variable = "Displ"
Variable DOFs = 3
then the solver of Elmer knows, in addition to the three-component vector Displ, three scalar fields Displ
1,Displ 2 and Displ 3. These scalar fields may be used as independent variables in parameter defini-
tions and used in the definitions of initial and boundary conditions, etc.
More complicated functions can be defined using MATC language. Here the basic usage of MATC in
connection with the solver input file is illustrated; for an additional documentation, see a separate manual
for MATC. For example, one may define
Density = Variable Temperature
MATC "1000*(1-1.0e-4*(tx-273))"
This means that the parameter Density depends on the value of Temperature as
ρ=ρ0(1 −β(T−T0)),(2.1)
with ρ0= 1000, β = 10−4and T0= 273. Note that the value of the independent variable is known as tx
in the MATC expression.
If the independent variable has more than one component, the variable tx will contain all the compo-
nents in values tx(0),tx(1),...,tx(n-1), where nis the number of the components of the independent
variable. A MATC expression may also take several scalar arguments; one may define, for example,
My Parameter = Variable Time, Displ 1
Real MATC "..."
The values of the scalar fields Time and Displ 1 can then be referred in the associated MATC expression
by the names tx(0) and tx(1), respectively.
In addition to using MATC functions, Fortran 90 functions may also be used to create parameter defini-
tions, etc. In the same manner as MATC functions are used, we may define
Density = Variable Temperature
Procedure "filename" "proc"
CSC – IT Center for Science
2. Structure of the Solver Input File 16
In this case the file "filename" should contain a shareable .so (Unix) or .dll (Windows) code for the user
function whose name is "proc". The call interface for the Fortran function is as follows
FUNCTION proc( Model, n, T ) RESULT(dens)
USE DefUtils)
IMPLICIT None
TYPE(Model_t) :: Model
INTEGER :: n
REAL(KIND=dp) :: T, dens
dens = 1000*(1-1.0d-4(T-273.0d0))
END FUNCTION proc
The Model structure contains pointers to all information about the model, and may be used to obtain field
variable values, node coordinates, etc. The argument n is the index of the node to be processed, and T is the
value of the independent variable at the node. The function should finally return the value of the dependent
variable.
The independent variable can also be composed of several independent components. We may thus define
Density = Variable Coordinate
Procedure "filename" "proc"
Now the argument T in the Fortran function interface should be a real array of three values, which give the
x,y and z coordinates of the current node.
Parameterized keyword commands
The solver input file also offers possibilities for creating parameterized commands that utilize MATC. In the
solver input file an expression following the symbol $ is generally interpreted to be in MATC language. If
the solver input file contains the lines
$solvertype = "Iterative"
$tol = 1.0e-6
then one may define, e.g.,
Solver 1
...
Linear System Solver = $solvertype
Linear System Convergence Tolerance = $tol
...
End
Solver 2
...
Linear System Solver = $solvertype
Linear System Convergence Tolerance = $100*tol
...
End
Alternative keyword syntax
Some keyword commands may be given by using an alternative syntax that may sometimes be needed. The
size of an integer or real number arrays may be given in parenthesis in connection with the keyword, but also
with the Size declaration. Therefore the following are exactly the same
Timestep Intervals(3) = 1 10 100
Timestep Intervals = Size 3; 1 10 100
CSC – IT Center for Science
2. Structure of the Solver Input File 17
This feature is useful when giving vectors and matrices in ElmerGUI since there the keyword body is fixed
and cannot include any size declaration. Note that in the above the semicolon is used as an alternative
character for newline.
Another convention is to use two colons to make in-lined definitions in the sif files. The following two
expressions are equal
Body Force 1
Heat Source = 1.0
End
and
Body Force 1 :: Heat Source = 1.0
2.4 Running several sequences
Execution within command file
When reading the string RUN in the command file, the solver stops the reading and performs the computation
with the instructions so far obtained. After a successful execution the solver continues to interpret the
command file. Using this functionality it is therefore possible to create scripts where some parameter value
is changed and the problem is recomputed. For example, adding the same sequence to the end of the sif file
could be used to test the solution with another linear solver
RUN
Solver 1::Linear System Iterative Method = BiCgstabl
RUN
It should be noted that not quite all features support this procedure. For example, some preconditioners
create static structures that will not be recreated.
CSC – IT Center for Science
Chapter 3
Restart from existing solutions
Often, the user wants to restart a run. This may be needed simply to continue an interrupted simulation
— no matter what caused the interruption — but also to read in values needed either as initial conditions or
as boundary conditions. These features are controlled by some keywords in the Simulation section that
will be presented in the following sections.
3.1 Restart file
Any output file obtained by using the keyword command
Output File = String ...
can be used to define a restart point for a new simulation. By convention the suffix of the file has been
.result but basically it may be chosen by the user.
By default all distributed fields will be saved to the restart file. Additionally the user may save global
variables with the following keyword
Output Global Variables = Logical
If one wants to select the fields to be saved they may be given individually by
Output Variable i = String
where i= 1,2,3,....
The format of the output is either ascii or binary. Ascii is the default, and binary format is enforced by
setting
Binary Output = True
The limitation of the restart functionality is that the mesh on which the previous case has been run must
be identical to that on which the new run is performed. In parallel runs, additionally, also the partitions of
the mesh have to coincide. The permuation vectors may vary, and also the field variables do not have to be
the same.
The commands for restarting are then given in the Simulation section by declaring the restart file
name as well as a restart position. For example we may specify
Simulation
Restart File = "previousrun.result"
Restart Position = 101
End
This would perform the current simulation by restarting from the time/iteration level 101 of the previously
stored result file previousrun.result.
If one wants to select the fields to be restarted they may be given individually by
CSC – IT Center for Science
3. Restart from existing solutions 19
Restart Variable i = String
where i= 1,2,3,....
Upon running the new simulation, a message similar to the following example should be seen in the
standard output of Elmer:
LoadRestartFile:
LoadRestartFile: --------------------------------------------
LoadRestartFile: Reading data from file: previousrun.result
LoadRestartFile: BINARY 3.L
LoadRestartFile:
LoadRestartFile: Total number of dofs to load: 13
LoadRestartFile: Reading time sequence: 2.000E-03
LoadRestartFile: Reading timestep: 10
LoadRestartFile: Time spent for restart (s): 7.9990E-03
LoadRestartFile: All done
LoadRestartFile: --------------------------------------------
If the number of stored time/iteration levels is not known a priori, the user can insert the command
Restart Position = 0
in order to make sure that the results for the lastly stored time/iteration level are reloaded.
Result files arising from steady state simulations often contain results for multiple iteration steps (with
the result for the last step containing the converged solution). Nevertheless, these instances of solutions are
— if reloaded — interpreted to describe the solution at different time levels. In this case the user might want
to redefine the value of time variable for the restarted simulation, especially if continuing with transient runs.
The keyword command
Restart Time = Real ...
may be given in order to manually set the time to correspond the zeroth time level of the new simulation.
3.2 Initialization of dependent variables
The initialization of variables and their boundary conditions is done by default before reading in the previous
results. That has two main implications:
1. Values set in the section Initial Condition are overwritten by the corresponding values of the
variable loaded afterwards from the restart file
2. Variable values given as initial or boundary conditions and specified to depend on other variables are
not initiated with those values from the restart file by default.
The latter can be influenced with two keywords, Restart Before Initial Conditions (by de-
fault False) and Initialize Dirichlet Condition (by default True).
By setting
Restart Before Initial Conditions = Logical True
Elmer would first load the variables from the restart file and then apply initial conditions to those variables
that have not been set by the earlier solution. This is necessary if one of the initial conditions depends on
the earlier solution. By default, first the initial conditions from the solver input file are set and thereafter the
restart file (if existing) is read.
The value of the keyword Initialize Dirichlet Condition is by default set to be true, which
means that Dirichlet conditions are set before the simulation and thus also before the particular solver for
handling that variable is executed. If a boundary condition for a certain variable now depends on the value
of another, the first-time Dirichlet condition is set using the initial values of variables — either set or read in
from the restart file. If this is not wanted, the user can switch to using the option
CSC – IT Center for Science
3. Restart from existing solutions 20
Initialize Dirichlet Condition = False
which will set the Dirichlet condition on the fly during the execution of the solver associated with the
variable.
CSC – IT Center for Science
Chapter 4
Solution methods for linear systems
4.1 Introduction
Discretization and linearization of a system of partial differential equations generally leads to solving linear
systems
Ax =b, (4.1)
where Aand bare of orders n×nand n×1, respectively. The coefficient matrix Aresulting from the finite
element discretization has a specific feature that the matrix is sparse, i.e. only a few of the matrix entries in
each row differ from zero. In many applications the system can also have a very large order n, so that the
chief part of the computation time in performing the simulation is typically spent by solvers for the linear
systems.
Solution methods for linear systems fall into two large categories: direct methods and iterative methods.
Direct methods determine the solution of the linear system exactly up to a machine precision. They per-
form in a robust manner leading to the solution after a predetermined number of floating-point operations.
Nevertheless, the drawback of direct methods is that they are expensive in computation time and computer
memory requirements and therefore cannot be applied to linear systems of very large order. The efficient
solution of large systems generally requires the use of iterative methods which work by generating sequences
of (hopefully) improving approximate solutions.
ElmerSolver provides access to both direct and iterative methods. The iterative methods available fall
into two main categories: preconditioned Krylov subspace methods and multilevel methods. Iteration meth-
ods that combine the ideas of these two approaches may also be constructed. Such methods may be very
efficient leading to a solution after a nearly optimal number of operation counts.
The development of efficient solution methods for linear systems is still an active area of research, the
amount of literature on the topic being nowadays vast. The aim of the following discussion is to provide
the user the basic knowledge of the solution methods available in ElmerSolver. The detailed description of
methods is omitted. For a more comprehensive treatment the reader is referred to references mentioned.
4.2 Direct methods
A linear system may be solved in a robust way by using direct methods. ElmerSolver offers two main options
for using direct methods. The default method utilizes the well-known LAPACK collection of subroutines for
band matrices. In practice, this solution method can only be used for the solution of small linear systems as
the operation count for this method is of order n3.
The other direct solver employs the Umfpack routines to solve sparse linear systems [1]. Umfpack uses
the Unsymmetric MultiFrontal method. In practice it may be the most efficient method for solving 2D
problems as long as there is enough memory available.
It should be noted that the success of the direct solvers depends very much on the bandwidth of the sparse
matrix. In 3D these routines therefore usually fail miserably.
CSC – IT Center for Science
4. Solution methods for linear systems 22
Elmer may be also compiled with Mumps,SuperLU, and Pardiso. The licensing scheme of these
software do not allow the distribution of precompiled binaries but everybody may themselves compile a
version that includes these solvers. Many times the best linear solver for a particular problem may be found
among these.
4.3 Preconditioned Krylov methods
ElmerSolver contains a set of Krylov subspace methods for the iterative solution of linear systems. These
methods may be applied to large linear systems but rapid convergence generally requires the use of precon-
ditioning.
4.3.1 Krylov subspace methods
The Krylov subspace methods available in ElmerSolver are
•Conjugate Gradient (CG)
•Conjugate Gradient Squared (CGS)
•Biconjugate Gradient Stabilized (BiCGStab)
•BiCGStab(ℓ)
•Transpose-Free Quasi-Minimal Residual (TFQMR)
•Generalized Minimal Residual (GMRES)
•Generalized Conjugate Residual (GCR)
Both real and complex systems can be solved using these algorithms. For the detailed description of some
of these methods see [3] and [4].
A definite answer to the question of what is the best iteration method for a particular case cannot be
given. In the following only some remarks on the applicability of the methods are made.
The CG method is an ideal solution algorithm for cases where the coefficient matrix Ais symmetric and
positive definite. The other methods may also be applied to cases where Ais non-symmetric. It is noted that
the convergence of the CGS method may be irregular. The BiCGStab and TFQMR methods are expected
to give smoother convergence. In cases where BiCGStab does not work well it may be advantageous to use
the BiCGStab(ℓ) method, with ℓ≥2a parameter. Faster convergence in terms of iteration counts may be
expected for increasing values of the parameter ℓ. However, since more work is required to obtain the iterate
as ℓincreases, optimal performance in terms of computational work may be realized for quite a small value
of ℓ. Starting with the value ℓ= 2 is recommended.
The GMRES and GCR methods generate gradually improving iterates that satisfy an optimality condi-
tion. The optimality may however come with a significant cost since the computational work and computer
memory requirements of these methods increase as the number of iterations grows. In practice these meth-
ods may be restarted after msolution updates have been performed in order to avoid the increasing work
and storage requirements. The resulting methods are referred to as the GMRES(m) and GCR(m) meth-
ods. Here the choice of mhas to be controlled by the user. It should be noted that the convergence of the
restarted algorithms may be considerably slower than that of full versions. Unfortunately, general guidelines
for determining a reasonable value for mcannot be given as this value is case-dependent.
4.3.2 Basic preconditioning strategies
The performance of iterative Krylov methods depends greatly on the spectrum of the coefficient matrix
A. The rate at which an iteration method converges can often be improved by transforming the original
system into an equivalent one that has more favorable spectral properties. This transformation is called
preconditioning and a matrix which determines the transformation is called a preconditioner.
CSC – IT Center for Science
4. Solution methods for linear systems 23
In ElmerSolver preconditioning is usually done by transforming (4.1) into the system
AM−1z=b, (4.2)
where the preconditioner Mis an approximation to Aand zis related to the solution xby z=M x. In
practice, the explicit construction of the inverse M−1is not needed, since only a subroutine that for a given
vreturns a solution uto the system
Mu =v(4.3)
is required.
The solver of Elmer provides several ways to obtain preconditioners. The basic strategies described in
this section may be used to form the Jacobi preconditioner and incomplete factorization preconditioners.
In addition to these preconditioners which remain stationary during the linear system iteration, the solver
of Elmer allows the use of variable preconditioning where the action of the preconditioner is defined in
terms of applying another iteration method to auxiliary preconditioning systems. Such inner-outer iterations
are typically applied in connection with block preconditioners and preconditioning by multilevel methods.
These advanced preconditioning strategies are considered in Section 4.5 below.
The Jacobi preconditioner is simply based on taking Mto be the diagonal of A. More sophisticated pre-
conditioners may be created by computing incomplete LU factorizations of A. The resulting preconditioners
are referred to as the ILU preconditioners. This approach gives the preconditioner matrix Min the form
M=LU where Land Uare lower and upper triangular with certain elements that arise in the factorization
process ignored.
There are several ways to choose a set of matrix positions that are allowed to be filled with nonzero
elements. ILU preconditioners of fill level Nreferred to as the ILU(N) preconditioners are built so that
ILU(0) accepts nonzero elements in the positions in which Ahas nonzero elements. ILU(1) allows nonzero
elements in the positions that are filled if the first step of Gaussian elimination is performed for A. ILU(2)
accepts fill in positions that are needed if the next step of Gaussian elimination is performed with ILU(1)
factorization, etc.
Another approximate factorization strategy is based on numerical tolerances. The resulting precondi-
tioner is referred to as the ILUT preconditioner. In the creation of this preconditioner Gaussian elimination
is performed so that elements of a given row of the LU factorization are obtained but only elements whose
absolute value (scaled by the norm of all values of the row) is over a given threshold value are accepted in
the preconditioner matrix.
Obviously, the additional computation time that is spent in creating the preconditioner matrix and solving
systems of the type (4.3) should be compensated by faster convergence. Finding an optimal ILU precondi-
tioner for a particular case may require the use of trial and error. Start with ILU(0) and try to increase the
fill level N. As N increases, more and more elements in the incomplete LU factorization of the coefficient
matrix are computed, so the preconditioner should in principle be better and the number of iterations needed
to obtain a solution should decrease. At the same time the memory usage grows rapidly and so does the time
spent in building the preconditioner matrix and in applying the preconditioner during iterations. The same
applies to the ILUT preconditioner with decreasing threshold values.
4.4 Multilevel methods
A class of iterative methods referred to as multilevel methods provides an efficient way to solve large linear
systems. For certain class of problems they perform nearly optimally, the operation count needed to obtain a
solution being nearly of order n. Two different multilevel-method approaches are available in ElmerSolver,
namely the geometric multigrid (GMG) and algebraic multigrid (AMG).
4.4.1 Geometric multigrid
Given a mesh T1for the finite element discretization of problem the geometric multigrid method utilizes
a set of coarser meshes Tk,k= 2, ..., N, to solve the linear system arising from the discretization. One
of the fundamental ideas underlying the method is based on the idea of coarse grid correction. That is, a
CSC – IT Center for Science
4. Solution methods for linear systems 24
coarser grid is utilized to obtain an approximation to the error in the current approximate solution of the
linear system. The recursive application of this strategy leads us to multigrid methods.
To utilize different meshes multigrid methods require the development of methods for transferring vec-
tors between fine and coarse meshes. Projection operators are used to transfer vectors from a fine mesh Tk
to a coarse mesh Tk+1 and will be denoted by Ik+1
k, while interpolation operators Ik
k+1 transfer vectors from
a coarse mesh to a fine mesh.
The multigrid method is defined by the following recursive algorithm: Given A,band an initial guess y
for the solution of the system Ax =bset i= 1 and do the following steps:
1. If i=N, then solve the system Ax =bby using the direct method and return.
2. Do pre-smoothing by applying some iterative algorithm for a given number of times to obtain a new
approximate solution y.
3. Perform coarse grid correction by starting a new application of this algorithm with A=Ii+1
iAIi
i+1,
b=Ii+1
i(Ay −b),i=i+ 1 and the initial guess e= 0.
4. Compute a new approximate solution by setting y=y+Ii
i+1e
5. Do post-smoothing by applying some iterative algorithm for a given number of times to obtain a new
approximate solution y.
6. If the solution has not yet converged, go to step 2.
In ElmerSolver one may choose the Jacobi, CG or BiCGStab algorithm as the method for smoothing itera-
tions.
The full success of multigrid methods is based on the favorable combination of the properties of ba-
sic iteration methods and methods for transferring vectors between meshes. The smoothing iterations give
rapid convergence for oscillatory solution components while coarse grid correction entails an efficient solu-
tion method for smooth solution components. For a comprehensive introduction to the geometric multigrid
method the reader is referred to [2].
4.4.2 Algebraic multigrid
In many cases the geometric multigrid may not be applied because we do not have the luxury of having
a set of appropriate hierarchical meshes. The alternative is the algebraic multigrid (AMG) method which
uses only the matrix Ato construct the projectors and the coarse level equations. AMG is best suited for
symmetric and positive semidefinite problems. For other types of problems the standard algorithm may fail.
For more information on AMG see reference [5].
The AMG method has two main phases. The set-up phase includes the recursive selection of the coarser
levels and definition of the transfer and coarse-grid operators. The solution phase uses the resulting compo-
nents to perform a normal multigrid cycling until a desired accuracy is reached. The solution phase is similar
to that of the GMG.
Note that the AMG solvers in ElmerSolver are not fully mature. They may provide good solutions for
some problems while desperately failing for others.
Classical Ruge-Stüben algorithm
The coarsening is performed using a standard Ruge-Stüben coarsening algorithm. The possible connections
are defined by the entries in the matrix A. The variable iis strongly coupled to another variable jif
aij <−c−max |aik|or aij > c+max |aik|,(4.4)
where 0< c−<1and 0< c+<1are parameters. Typically c−≈0.2and c+≈0.5. Once the negative
(P−) and positive (P+) strong couplings have been determined the variables are divided into coarse (C) and
fine (F) variables using the standard coarsening scheme.
CSC – IT Center for Science
4. Solution methods for linear systems 25
The interpolation matrix may be constructed using the C/F -splitting and the strong couplings of the
matrix. The interpolation of coarse nodes is simple as they remain unchanged. The interpolation of fine
nodes starts from the fact the smooth error emust roughly satisfy the condition Ae = 0 or
aiiei+X
j6=i
aij ej= 0.(4.5)
By manipulation
aiiei+αiX
j∈C∩P−
i
aij ej+βiX
j∈C∩P+
i
aij ej= 0,(4.6)
where
αi=Pj∈P−
iaij
Pj∈C∩P−
iaij
and βi=Pj∈P+
iaij
Pj∈C∩P+
iaij
.(4.7)
The interpolation thus becomes
ei=X
j∈C∩Pi
wij ejwith wij =−αiaij /aii, j ∈P−
i,
−βiaij /aii, j ∈P+
i.(4.8)
This is known as direct interpolation. It may be modified by using also the strong F-nodes in the
interpolation. This means that in formula (4.5) the following elimination is made for each j∈F∩Pi
ej→ − X
k∈C∩Pj
ajkek/ajj .(4.9)
This is known as standard interpolation. In practice it means that the number of nodes used in the interpo-
lation is increased. This may be important to the quality of the interpolation particularly if the number of
direct C-neighbors is small.
After the interpolation weights have been computed the smallest coefficients may be truncated if they
are small, i.e.,wj< cwmax |wk|, where cw≈0.2. The other values must accordingly be increased so that
the sum of weights remains constant. The truncation is essential in preventing the filling of the coarse level
matrices.
Cluster multigrid
There is also an implementation of the agglomeration or cluster multigrid method. It is a variant of the
algebraic multigrid method. In this method the components are grouped and the coarse-level matrices are
created simply by summing up the corresponding rows and columns. In other words, the projection matrix
includes just ones and zeros.
The cluster multigrid method should be more robust for problems where it is difficult to generate an
optimal projection matrix. However, for simple problems it is usually beaten by the standard Ruge-Stüben
method.
4.5 Preconditioning via inner iterations
In the following, we return to the subject of preconditioning and consider the construction of more advanced
preconditioners that utilize the idea of variable preconditioning. We recall that the procedural description
of preconditioning allows us to consider schemes where the action of M−1in (4.2) is replaced by perform-
ing an inexact solution of the preconditioning system of the type (4.3) with other iterative methods. In this
connection, the iteration steps of the method applied to the preconditioning system are referred to as inner
iterations, while the iteration steps of the method which is being preconditioned are referred to as outer iter-
ations. In addition, the concept of variable preconditioning is motivated by the fact that the preconditioner
can change during the outer iteration process. It should be noted that all the Krylov methods do not neces-
sarily admit using variable preconditioning. Anyhow, GCR algorithm can always be employed as the outer
iterative method, and we generally recommend using it in connection with inner-outer iterations.
CSC – IT Center for Science
4. Solution methods for linear systems 26
4.5.1 Preconditioning by multilevel methods
Multilevel methods are iteration methods on their own but they can also be applied to define preconditioners
for the Krylov subspace methods. This preconditioning approach corresponds to taking M=Ain (4.3) and
performing an inaccurate solution of the resulting system using multilevel methods to obtain u. A rather
mild stopping criterion may be used in this connection. Preconditioning by multilevel methods may lead to
very efficient solution methods for large linear systems.
4.5.2 Block preconditioning
Performing preconditioning by solving auxiliary systems iteratively also arises naturally in connection with
a strategy referred to as the block preconditioning. Tacit in the concept of block preconditioning is that
the unknowns of the linear system are grouped into sets of variables so that a natural block partitioning of
the linear system is induced. Often such partitioning is naturally induced by physics. Block factorizations,
block matrix splittings, or even other iteration schemes based on updating different types of variables inde-
pendently can then be utilized to derive preconditioners. The bottom line is that such block preconditioners
typically lead to solving subsidiary problems which can be handled by employing inner iterations.
The designs of efficient block preconditioners tend to involve problem-specific features, especially, if
block factorizations or other segregated iteration methods are used to derive preconditioners. Therefore the
model solvers in Elmer that offer the possibility of using block preconditioning have originally been devel-
oped independently, with the block preconditioning ability being a part of the solver-level implementation.
Recently, routines which are suited for constructing block preconditioners on a more general level have also
been implemented.
4.6 Keywords related to linear system solvers
The following keywords may be given in Solver section of the solver input file (.sif file).
Linear System Solver String
Using this keyword the type of linear system solver is selected. This keyword may take the following
values:
•Direct
•Iterative
•Multigrid
Here Iterative and Multigrid refer to the Krylov and multilevel methods, respectively.
Linear System Direct Method String
If the value of the Linear System Solver keyword is set to be Direct, one may choose a band
matrix solver with the value Banded or a sparse matrix solver with the value umfpack,mumps,
Pardiso or superlu, . The default is Banded.
Linear System Iterative Method String
If the value of the Linear System Solver keyword is set to be Iterative, one should choose
a Krylov method by setting the value of this keyword to be one of the following alternatives:
•CG
•CGS
•BiCGStab
•BiCGStabl
•TFQMR
•GMRES
•GCR
CSC – IT Center for Science
4. Solution methods for linear systems 27
See also the MG Smoother keyword.
Linear System GMRES Restart Integer [10]
The restart parameter mfor the GMRES(m) method may be given using this keyword.
Linear System GCR Restart Integer
The restart parameter mfor the GCR(m) method may be given using this keyword. The default option
is that restarting is not performed, i.e. the full GCR is used.
BiCGstabl polynomial degree Integer
The parameter ℓfor the BiCGStab(ℓ) method may be given. By default the minimal applicable value
ℓ= 2 is used.
Linear System Preconditioning String
A preconditioner for the Krylov methods may be declared by setting the value of this keyword to be
one of the following alternatives:
•None
•Diagonal
•ILUn, where the literal nmay take values 0,1,...,9.
•ILUT
•Multigrid
See also the MG Preconditioning keyword.
Linear System ILUT Tolerance Real [0.0]
This keyword is used to define the value of the numerical tolerance for the ILUT preconditioner.
Linear System Convergence Tolerance Real [0.0]
This keyword is used to define a stopping criterion for the Krylov methods. The approximate solution
is considered to be accurate enough if the iterate satisfies
||Ax −b||
||b|| ≤ǫ
where ǫis the value of this keyword. See also MG Tolerance.
Linear System Max Iterations Integer [0]
This keyword is used to define the maximum number of the iterations the Krylov methods are permit-
ted to perform. If this limit is reached and the approximate solution does not satisfy the stopping crite-
rion, ElmerSolver either continues the run using the current approximate solution as the solution of the
system or aborts the run depending on the value of Linear System Abort Not Converged
keyword. See also MG Max Iterations keyword.
Linear System Abort Not Converged Logical [True]
If the value of this keyword is set to be True, ElmerSolver aborts the run when the maximum number
of iterations the algorithm is permitted to perform is reached and the approximate solution does not
satisfy the stopping criterion. Otherwise the run will be continued using the current approximate
solution as the solution of the system (this may lead to troubles at later steps of computation).
Linear System Residual Mode Logical
Usually the solvers of Elmer are written so that the linear system is solved for the full amplitude of the
solution. However, often it is numerically to solve the residual equation where a change in the system
of solved such that the residual would be minimized. This has the advantage that the linear system
does not need to be solved as accurately. In practive we take the linear equation Ax =band make a
change of variables such that A dx =r, where r=b−Ax0and dx =x1−x0. After solution of the
modified equation we revert back to the original system by setting x1=x0+dx.
CSC – IT Center for Science
4. Solution methods for linear systems 28
Linear System Residual Output Integer [1]
By default the iterative algorithms display the value of the (scaled) residual norm after each iteration
step. Giving a value n > 1for this keyword may be used to display the residual norm only after every
n iterations. If the value 0 is given, the residual output is disabled.
Linear System Precondition Recompute Integer [1]
By default the ElmerSolver computes the preconditioner when a new application of iterative algorithm
is started. If the value of this keyword is set to be n, the preconditioner is computed only after
n successive subroutine calls for linear systems arising from same source. This may speed up the
solution procedure especially in cases where the coefficient matrix does not change much between
successive subroutine calls. On the other hand if the coefficient matrix has changed significantly, the
preconditioner may not be efficient anymore.
Optimize Bandwidth Logical [True]
If the value of this keyword is set to be True, the Cuthill-McKee bandwidth optimization scheme is
used to order the unknowns in such a way that band matrices can be handled efficiently. The bandwidth
optimization is recommended when the direct solver or incomplete factorization preconditioners are
used.
The keywords beginning with MG are activated only if either the Linear System Solver or Linear
System Preconditioning keyword has the value Multigrid. If a multigrid method is used as the
linear system solver, some keywords starting with MG may be replaced by corresponding keywords starting
with phrase Linear System. It should be noted that in the case of a multigrid solver there are some
limitations to what values the keywords starting with the phrase Linear System may take, see below.
MG Levels Integer [1]
This keyword is used to define the number of levels for the multigrid method.
MG Equal Split Logical [False]
A hierarchy of meshes utilized by the multigrid method may be generated automatically by setting the
value of this keyword to be True. The coarsest mesh must be supplied by the user and is declared in
the usual way in the Header section of the solver input file. The other meshes are obtained using an
equal division of the coarse mesh. The solution of the problem will be sought for the finest mesh.
MG Mesh Name File
A hierarchy of meshes utilized by the multigrid method may be supplied by the user. A base name of
the mesh directories is declared using this keyword. The names of mesh directories must be composed
of the base name appended with a level number such that if the base name is mgridmesh, the mesh
directories should have names mgridmesh1,mgridmesh2, etc. The meshes are numbered starting
from the coarsest mesh. In addition, the finest mesh must be declared in the Header section of the
solver input file. It should be noted that the MG Equal Split keyword must be set to be False
to enable the use of user-supplied meshes.
MG Max Iterations Integer [0]
If a multigrid method is used as a preconditioner for the Krylov methods, the value of this keyword
defines the maximum number of iterations the multigrid method is allowed to perform to solve the
preconditioning equation (4.3). Usually one or two iterations are sufficient. If a multigrid method
is the linear system solver, the use of this keyword is similar to that of the Linear System Max
Iterations keyword.
MG Convergence Tolerance Real [0.0]
If a multigrid method is used as a preconditioner for the Krylov methods, this keyword defines the
solution accuracy for the preconditioning equation (4.3). This keyword is not usually needed if the MG
Max Iterations keyword has a small value. If a multigrid method is the linear system solver,
the use of this keyword is similar to that of the Linear System Convergence Tolerance
keyword.
CSC – IT Center for Science
4. Solution methods for linear systems 29
MG Smoother String
This keyword defines the algorithm for pre- and post-smoothing. It may take one of the following
values:
•Jacobi
•CG
•BiCGStab
If the linear system solver is a multigrid method, the Linear System Iterative Method
keyword may be used instead of this keyword, but only the three algorithms mentioned here can be
applied.
MG Pre Smoothing Iterations Integer [0]
This keyword defines the number of pre-smoothing iterations.
MG Post Smoothing Iterations Integer [0]
This keyword defines the number of post-smoothing iterations.
MG Preconditioning String
This keyword declares the preconditioner for the algorithm which is used in smoothing iterations. It
may take one of the following values:
•None
•ILUn, where the literal nmay take values 0,1,...,9.
•ILUT
Note that this keyword is not related to using multigrid method as a preconditioner. It is also noted
that preconditioning the smoothing algorithms does not seem to work well if a multigrid method is
used as a preconditioner for Krylov methods.
MG ILUT Tolearance Real [0.0]
This keyword defines the numerical tolerance for the ILUT preconditioner in connection with smooth-
ing iterations.
The keywords for the algebraic multigrid solver are in a large part the same as for the geometric multigrid.
There are however some keywords that are related only to AMG.
MG Lowest Linear Solver Limit Integer
This value gives a lower limit for the set of coarse nodes after which the recursive multilevel routine
is terminated. A proper value might be around 100.
MG Recompute Projector Logical
This flag may be used to enforce recomputation of the projector each time the algebraic multigrid
solver is called. The default is False as usually the same projector is appropriate for all computations.
MG Eliminate Dirichlet Logical
At the highest level the fixed nodes may all be set to be coarse since their value is not affected by the
lower levels. The default is True
MG Eliminate Dirichlet Limit Real
Gives the maximum fraction of non-diagonal entries for a Dirichlet node.
MG Smoother String
In addition to the selection for the GMG option sor (symmetric over relaxation) is possible.
MG SOR Relax String
The relaxation factor for the SOR method. The default is 1.
CSC – IT Center for Science
4. Solution methods for linear systems 30
MG Strong Connection Limit Real
The coefficient c−in the coarsening scheme. Default is 0.25.
MG Positive Connection Limit Real
The coefficient c+in the coarsening scheme. Default is 1.0.
MG Projection Limit Real
The coefficient cwin the truncation of the small weights. The default is 0.1.
MG Direct Interpolate Logical
Chooses between direct and standard interpolation. The default is False.
MG Direct Interpolate Limit Integer
The standard interpolation may also be applied to nodes with only a small number of coarse connec-
tion. This gives the smallest number of nodes for which direct interpolation is used.
Finally, there are also some keywords related only to the clustering multigrid.
MG Cluster Size Integer
The desired choice of the cluster. Possible choices are 2,3,4,5,. . . and zero which corresponds to the
maximum cluster.
MG Cluster Alpha Real
In the clustering algorithm the coarse level matrix is not optimal for getting the correct convergence.
Tuning this value between 1 and 2 may give better performance.
MG Strong Connection Limit Real
This is used similarly as in the AMG method except it is related to positive and negative connections
alike.
MG Strong Connection Minimum Integer
If the number of strong connections with the given limit is smaller than this number then increase the
set of strong connection if available connections exist.
4.7 Implementation issues
4.7.1 The sparse matrix storage
To be efficient, iteration methods require that a matrix-vector product for sparse matrices is efficiently im-
plemented. A special storage scheme called the Compressed Row Storage (CRS) [3] is used in ElmerSolver
to store only those matrix coefficients that differ from zero.
The matrix structure is defined in module Types as:
TYPE Matrix_t
...
INTEGER :: NumberOfRows
REAL(KIND=dp), POINTER :: Values(:)
INTEGER, POINTER :: Rows(:), Cols(:), Diag(:)
...
END TYPE Matrix_t
The matrix type has several additional fields, but the basic storage scheme can be implemented using the
fields shown. The array Values is used to store the nonzero elements of the coefficient matrix. The array
Cols contains the column numbers for the elements stored in the array Values, while the array Rows
contains indices to elements that start new rows. In addition, Row(n+1) gives the number of nonzero
matrix elements + 1. The array Diag is used to store the indices of the diagonal elements.
For example, to go through the matrix row by row the following loop may be used
CSC – IT Center for Science
4. Solution methods for linear systems 31
USE Types
TYPE(Matrix_t), POINTER :: A
REAL(KIND=dp):: val
INTEGER :: i, j, row, col
DO i=1, A % NumberOfRows
PRINT *, ’Diagonal element for row ’, i, ’ is ’, A % Values( A % Diag(i) )
DO j=A % Rows(i), A % Rows(i+1)-1
row = i
col = A % Cols(j)
val = A % Values(j)
PRINT *, ’Matrix element at position: ’, row,col, ’ is ’, val
END DO
END DO
4.7.2 Subroutine calls
Most of the functionality of the sparse linear system solver of the ElmerSolver is available by using the
function call
Norm = DefaultSolve().
The return value Norm is a norm of the solution vector.
Sometimes it may be convenient to modify the linear system before solving it. A Fortran function which
performs this modification can be written by the user with the name of the function being declared in the
solver input file. For example, this technique may be used to define a user-supplied linear system solver.
If the name of the user-supplied Fortran function is proc and it can be found in the file having the name
Filename, the declaration
Before Linsolve File Filename proc
in the solver input file has the effect that the function will be called just before the default call of linear
system solver. The arguments the function can take are fixed and are declared as
FUNCTION proc( Model, Solver, A, b, x, n, DOFs, Norm ) RESULT(stat)
USE SolverUtils
TYPE(Model_t) :: Model
TYPE(Solver_t) :: Solver
TYPE(Matrix_t), POINTER :: A
REAL(KIND=dp) :: b(:), x(:), Norm
INTEGER :: n, DOFs, stat
...
END FUNCTION proc
Here the Model structure contains the whole definition of the elmer run. The Solver structure contains
information for the equation solver from which this linear system originates. The coefficient matrix Ais in
CRS format, bis the right-hand side vector, and xcontains the previous solution. The argument nis the
number of unknowns, and DOFs is the number of unknowns at a single node.
If the return value from this function is zero, the (possibly) modified linear system is solved after the
return. If the return value is 1, the linear system is assumed to be already solved and the vector xshould
contain the result. It is noted that the user-supplied Fortran function may also call the default linear equation
solver within the function, i.e. one may write the subroutine call
CALL SolveLinearSystem( A, b, x, Norm, DOFs, Solver )
Here Aand bmay be modified so that the linear system which is solved need not be same as the input system.
In a similar way the user may also supply a user-defined Fortran function which will be called just after
the solution of linear system. This is done using the declaration
CSC – IT Center for Science
BIBLIOGRAPHY 32
After Linsolve File Filename proc
in the solver input file. The arguments of this function are the same as for a function in the context of
Before Linsolve keyword.
Bibliography
[1] Umfpack home page. http://www.cise.ufl.edu/research/sparse/umfpack/.
[2] W.L. Briggs. A Multigrid Tutorial. SIAM, 1987.
[3] Richard Barrett et al. Templates for the Solution of Linear Systems: Building Blocks for Iterative Meth-
ods. SIAM, 1993.
[4] R.W. Freund. A transpose-free quasi-minimal residual algorithm for non-hermitian linear systems. SIAM
J. Sci. Comput., 14:470–482, 1993.
[5] K. Stüben. Algebraic Multigrid (AMG): An introduction with applications. GMD – Forschungszentrum
Informationstechnik GmbH, 1999.
CSC – IT Center for Science
Chapter 5
Nonlinear System Options
5.1 Introduction
Numerical methods in linear algebra are usually intended for the solution of linear problems. However,
there are many problems which are not linear in nature. The nonlinearity may a intrinsic characteristics
of the equation, such as is the case with intertial forces in the Navier-Stokes equation. The nonlinerity
might also a result of nonlinear material parameters that depend on the solution. What ever the reason for
nonlinearity the equations in Elmer are always first linearized to the form
A(ui−1)ui=b(ui−1),(5.1)
where irefers to the iteration cycle.
How the equations are linearized varies from solver toanother. For example, in the Navier-Stokes solver
there are tow different methods – the Picard linearization and the Newton linearization that may be used.
Also a hybrid scheme where the Picard type of scheme is switched to the Newton kind of scheme when
certain criteria are met is available. Therefore this section will not deal with the particular linearization
technique of different solver but tries to give some light to the generic keywords that are available. Some
keywords may also be defined in the Models Manual related to particular solvers.
In multiphysical simulations there are also a number of keywords related to the solution of coupled sys-
tems. Basically one may may define how many times a system of equations is solved repeatedly at maximum
and how what are the convergence criteria of the individual solvers that must be met simulataneously.
5.2 Keywords related to solution of nonlinear systems
These keywords are located in the Solver section of each solver, if requited at all.
Nonlinear System Convergence Measure String
The change of solution between two consecutive iterations may be estimated by a number of different
measures which are envoked by values norm,solution and residual. The default way of
checking for convergence is to test the change of norm
δ= 2 ∗ ||ui| − |ui−1||/(|ui|+|ui−1|).(5.2)
This measure is rather liberal since the norm of two solutions may be the same even though the
solutions would not. Therefore it is often desirable to look at the norm of change,
δ= 2 ∗ |ui−ui−1|/(|ui|+|ui−1|).(5.3)
The third choice is to use a backward norm of the residual where the old solution is used with the new
matrix.
δ=|Axi−1−b|/|b|.(5.4)
CSC – IT Center for Science
5. Nonlinear System Options 34
In the current implementation this norm lags one step behind and therefore always performs one extra
iteration.
Nonlinear System Norm Degree Integer
The choice of norms used in the evaluation of the convergence measures is not self evident. The
default is the L2norm. This keyword may be used to replace this by Ln norm where value n= 0
corresponds to the infinity (i.e. maximum) norm.
Nonlinear System Norm Dofs Integer
For vector valued field variables by default all components are used in the computation of the norm.
However, sometimes it may be desirable only to use some of them. This keyword may be used to give
the number of components used in the evaluation. For example, in the Navier-Stokes equations the
norm is only taken in respect to the velocity components while pressure is omitted.
Nonlinear System Convergence Absolute Logical
This keyword may be used to enforce absolute convergence measures rather than relative. The default
is False.
Nonlinear System Convergence Tolerance Real
This keyword gives a criterion to terminate the nonlinear iteration after the relative change of the norm
of the field variable between two consecutive iterations is small enough δ < ǫ, where ǫis the value
given with this keyword.
Nonlinear System Max Iterations Integer
The maxmimum number of nonlinear iterations the solver is allowed to do.
Nonlinear System Newton After Iterations Integer
Change the nonlinear solver type to Newton iteration after a number of Picard iterations have been
performed. If a given convergence tolerance between two iterations is met before the iteration count
is met, it will switch the iteration type instead. This applies only to some few solvers (as the Navier-
Stokes) where different linearization strategies are available.
Nonlinear System Newton After Tolerance Real
Change the nonlinear solver type to Newton iteration, if the relative change of the norm of the field
variable meets a tolerance criterion:
δ < ǫ,
where ǫis the value given with this keyword.
Nonlinear System Relaxation Factor Real
Giving this keyword triggers the use of relaxation in the nonlinear equation solver. Using a factor
below unity is sometimes required to achive convergence of the nonlinear system. Typical values
range between 0.3 and unity. If one must use smaller values for the relaxation factor some other
methods to boost up the convergence might be needed to improve the convergence. A factor above
unity might rarely speed up the convergence. Relaxed variable is defined as follows:
u′
i=λui+ (1 −λ)ui−1,
where λis the factor given with this keyword. The default value for the relaxation factor is unity.
Many of the keywords used to control the Nonlinear System have a corresponding keyword for the
Steady State. Basically the operation is similar except the reference value for the current solution uiis the
last converged value of the nonlinear system before starting a new loosely coupled iteration cycle. Otherwise
the explanations given above are valid.
Steady State Convergence Measure String
Steady State Norm Degree Integer
Steady State Norm Dofs Integer
CSC – IT Center for Science
5. Nonlinear System Options 35
Steady State Convergence Tolerance Real
Steady State Relaxation Factor Real
Additionally these keywords are located in the Simulation section of the command file.
Steady State Max Iterations Integer
The maximum number of coupled system iterations. For steady state analysis this means it litelarly,
for transient analysis this is the maximum number of iterations within each timestep.
Steady State Min Iterations Integer
Sometimes the coupling is such that nontrivial solutions are obtained only after some basic cycle
is repeated. Therefore the user may sometimes need to set also the minimum number of iterations.
Sometimes the steady state loop is also used in a dirty way to do some systematic procedures – for
example computing the capacitance matrix, or lumped elastic springs. Then this value may be set to
an a priori known constant value.
CSC – IT Center for Science
Chapter 6
Integration of time-dependent systems
6.1 Introduction
Solving time-dependent systems is becoming more and more common in various branches of computational
science, as the computer resources grow steadily. ElmerSolver may be adapted to solve such systems. The
first order time derivatives may be discretized by using the following methods:
•the Crank-Nicolson method
•the Backward Differences Formulae (BDF) of several orders
In the case of the first order BDF scheme adaptive time-stepping strategy may also be used.
The second order time derivatives are approximated by either using the Bossak method or reformulating
the second order equations as equivalent systems of first order equations.
In addition, a predictor-corrector scheme is available to automatically control the time steps according to
the local truncation error. The analysis and implementation details can be seen in [1].
6.2 Time discretization strategies
Consider the numerical solution of the evolutionary field equation
∂φ
∂t +Kφ=f, (6.1)
where the differential operator Kdoes not involve differentiation with respect to time tand fis a given
function of spatial coordinates and time. The spatial discretization of (6.1) leads to the algebraic equations
M∂Φ
∂t +KΦ = F, (6.2)
where M,Kand Fresult from the discretization of the identity operator, the operator Kand f, respectively.
The vector Φcontains the values of the unknown field φat nodes.
The applications of the first three BDF methods to discretizate the time derivative term in (6.2) yield the
following systems: 1
∆tM+KΦi+1 =Fi+1 +1
∆tMΦi,(6.3)
1
∆tM+2
3KΦi+1 =2
3Fi+1 +1
∆tM4
3Φi−1
3Φi−1,(6.4)
1
∆tM+6
11KΦi+1 =6
11Fi+1 +1
∆tM18
11Φi−9
11Φi−1+2
11Φi−2,(6.5)
CSC – IT Center for Science
6. Integration of time-dependent systems 37
where ∆tis the time step and Φiis the solution at time step i. Similarly, Fiis the value of Fat time step i.
All the BDF methods are implicit in time and stable. The accuracies of the methods increase along with
the increasing order. The starting values for the BDF schemes of order n > 1are computed using the BDF
schemes of order 1, ..., n −1as starting procedures. It should be noted that the BDF discretizations of order
n > 3do not allow the use of variable time-step size. Adaptive time-stepping strategy may also be used in
the case of the first order BDF scheme.
The adaptive time-stepping is accomplished by first solving the system using a trial time step and then
using two time steps the lengths of which equal to the half of that of the trial time step and comparing the
results. If the difference between the results is found to be sufficiently small, the use of the trial time step is
accepted. Otherwise a new trial time step is defined by dividing the previous trial time step into two steps
of equal length and then the procedure is repeated. One may define one’s own criterion for determining
whether the use of the current time step is accepted. The default criterion is that the norms of the solutions
to each system of field equations do not differ more than the given threshold value.
The time discretization of the second order equation
M∂2Φ
∂t2+B∂Φ
∂t +KΦ = F(6.6)
using the Bossak method leads to the system
1−α
β(∆t)2M+γ
β∆tB+KΦi+1 =Fi+1 +M1−α
β(∆t)2Φi+γ
β∆tVi+(1 −α)
2βAi+
Bγ
β∆tΦi+γ
β−1Vi+1−γ
2β∆tAi,
(6.7)
where
Vi+1 =Vi+ ∆t(1 −γ)Ai+γAi+1,
Ai+1 =1
β(∆t)2(Φi+1 −Φi)−1
β∆tVi+1−1
2βAi,
β=1
4(1 −α)2, γ =1
2−α.
(6.8)
In the following the matrices Mand Bare referred to as the mass and damping matrix, respectively.
6.3 Keywords related to time discretization
All the keywords related to the time discretization may be given in Simulation section of the solver input file
(.sif file). A number of keywords may also be given in Solver section, so that each system of field equations
may be discretized using independently chosen time-stepping method. If keywords are not given in the
Solver section, the values of the keywords are taken to be those given in the Simulation section. It should
be noted that certain keywords such as those controlling the number of time steps, time-step sizes etc. may
only be given in the Simulation section.
Timestepping Method String
This keyword is used to declare the time discretization strategy for the first order equations. The value
of this keyword may be set to be either ”BDF” or ”Crank-Nicolson” and may be given in either
Simulation section or Solver section of the solver input file.
BDF Order Integer
This keyword is used to define the order of the BDF method and may take values 1,...,5. This keyword
may be given in either Simulation section or Solver section of the solver input file.
Time Derivative Order Integer
If a second order equation is discretized, this keyword must be given the value 2 in the Solver section
of the solver input file. It should be noted that the second order time derivatives are always discretized
using the Bossak method.
CSC – IT Center for Science
6. Integration of time-dependent systems 38
Bossak Alpha Real [-0.05]
This keyword is used to define the value for αin the Bossak method used in the time discretization of
second order equations. This keyword may be given in either Simulation section or Solver section of
the solver input file.
Timestep Intervals Integer array
This keyword is used to define the number of time steps. It may be array-valued so that different
time-step lengths may be used for different time intervals of the entire simulation. For example, if one
wishes to take first 50 time steps and then to use a different time-step length for the following 100 time
steps, one may define
Timestep Intervals(2) = 50 100
and use the Timestep Sizes keyword to define time-step lengths for the two sets of time steps.
Timestep Sizes Real array
This keyword is used to define the length of time step. If the value of the Timestep Intervals
keyword is array-valued, the value of this keyword must also be an array of the same size. For example,
if one has defined
Timestep Intervals(2) = 50 100
the declaration
Timestep Sizes(2) = 0.1 1.0
sets the time-step length for the first 50 time steps to be 0.1 and for the remaining 100 time steps 1.0.
Timestep Size Real
Instead of using the Timestep Sizes keyword the length of time step may be defined by using
this keyword (mind the singular). The value of this keyword is evaluated at the beginning of each time
step. A variable time-step length may conveniently be defined using a MATC or Fortran function.
This replaces the old keyword Timestep Function.
Output Intervals Integer array
This keyword is used to define the time-step interval for writing the results on disk. As in the case of
the Timestep Sizes keyword the size of the value of this keyword must be compatible with that
of the Timestep Intervals keyword. The value at a step mis saved if for the corresponding
output interval omod(m-1,o)==0. An exception is output interval equal to zero for which output is
not saved at all. However, the last step of the simulation is always saved.
Lumped Mass Matrix Logical [false]
The use of a lumped mass matrix may be activated by setting the value of this keyword to be True in
the Solver section of solver input file. The default lumping is defined by
M′
ii =Mii PiPjMij
PiMii
.(6.9)
The keywords related to the adaptive time-stepping may only be given in the Simulation section of the
solver input file. When the adaptive time-stepping strategy is used, a set of trial time steps is defined using
the keywords introduced above. The adaptive procedure is executed for each of these trial steps. Note that
the adaptive time-stepping is possible only in the case of the first order BDF scheme.
Adaptive Timestepping Logical [false]
The value of this keyword must be set to be True if the adaptive time integration is to be used.
CSC – IT Center for Science
6. Integration of time-dependent systems 39
Adaptive Time Error Real
This keyword is used to define the threshold value for the criterion for determining whether the use of
the current time step is accepted.
Adaptive Error Measure Real
Using this keyword one may define one’s own measure for evaluating the difference between the
computed results. This measure and the threshold value, which is given using the Adaptive Time
Error keyword, may be used to define a user-defined criterion for determining whether the use of the
current time step is accepted. The value of the Adaptive Error Measure keyword is evaluated
twice for each trial time step. For the first time the value of the keyword is evaluated after the system
is solved using the trial time step. The second time is after the system is solved using two time steps
the lengths of which equal to the half of that of the trial time step. The absolute value of the relative
difference between these two values is compared to the threshold value given by the Adaptive
Time Error keyword to determine whether the use of the current time step is accepted. If several
systems of field equations are solved, all the solutions must satisfy the similar criterion. If this keyword
is not used, the default criterion is based on comparing the norms of the solution fields.
Adaptive Min Timestep Real
Using this keyword one can limit the subsequent division of the trial time steps by giving the minimum
time-step length which is allowed.
Adaptive Keep Smallest Integer [1]
By default the adaptive scheme tries to double the length of the time step after the acceptable time
step is found. If a value n > 1is given for this keyword, the adaptive scheme tries to increase the step
length after taking n steps which are at most as long as the step length accepted.
6.4 On the treatment of time derivatives in Elmer Solver code
In the following a number of issues that may be useful if one is writing a code to solve one’s own application
are explained.
By default Elmer Solver does not generate or use global mass or damping matrices in the solution of
time-dependent systems. Mass and damping matrices need to be computed only element-wise, as the linear
system resulting from the time discretization, such as (6.3), is first formed element-wise and this local
contribution is later assembled to the global system. In the case of the first order equation (6.2) the local
linear system may be formed by using the subroutine call
CALL Default1stOrderTime( M, K, F ),
where Mis the element mass matrix, Kis the element stiffness matrix and Fis the element force vector. In
a similar manner, in the case of the second order equation (6.6) one may use the subroutine call
CALL Default2ndOrderTime( M, B, K, F ),
where Bis the element damping matrix.
Note that these subroutines must also be called for the local matrices and vectors that result from the
discretization of neumann and newton boundary conditions. If the boundary conditions do not contain any
time derivatives, the Mand Bmatrices should be set to be zero before calling the above subroutines.
If the global mass matrix is required, it may be generated by using the subroutine call
CALL DefaultUpdateMass( M )
Similarly, the global damping matrix may be generated by using the subroutine call
CALL DefaultUpdateDamp( B ).
Global mass (and possibly damping) matrices are required, for example, in the solution of eigenvalue prob-
lems. One may also implement one’s own time-stepping scheme at the global level using these matrices.
CSC – IT Center for Science
BIBLIOGRAPHY 40
6.5 Predictor Corrector Scheme
The keyword Timestep Intervals is used to set the number of steps and the keyword Timestep
Sizes is used to set the initial step size only. The step size will be changed automatically by the predictor-
corrector scheme based on an error estimator and PI controller. The solver used for predictor-corrector error
estimator is set by having
Exec Solver = "Predictor-Corrector"
There can be more than one solver with predictor-corrector scheme, but only the one with the smallest index
will be used for the error estimate and the time step will be changed accordingly. In order to use a predictor-
corrector scheme, one keyword should be given in the Simulation section and several keywords may also be
given in the Solver section.
Predictor-Corrector Control Logical
The predictor-corrector scheme may be activated by setting the value of this keyword to be True in
the Simulation section.
Skip First Timestep Logical
This keyword has to be set to True in the Solver section for most cases to get the solutions from the
other solvers which are used as the initial condition for the predictor-corrector solver. If this keyword
is false, the predictor-corrector solver will use the values from the Initial Condition section directly.
Predictor method String
This keyword is used to define the time scheme for the predictor in the Solver section. Currently, the
value may only be Adams-Bashforth. But, there will be more options in the future release.
Corrector method String
This keyword is used to define the time scheme for the corrector in the Solver section. Currently, the
value may only be Adams-Moulton. But, there will be more options in the future release.
Predictor-Corrector Scheme Order Integer[2]
This keyword is used to define the order of the predictor-corrector scheme. The value can be set to 1
or 2which corresponds the first or second order semi-implicit schemes mentioned in [1].
Predictor-Corrector Control Tolerance Real[1e-6]
This keyword is used to define the tolerance ǫfor the PI control process in [1]. The default value is
1e-6. A smaller value puts more restriction on the accuracy which may reduce the step sizes. A
larger value of the tolerance allows larger step sizes, but the solution may not be stable in some cases.
Predictor-Corrector Save Error Logical[false]
This keywords is used to save the time steps and local truncation errors for each step.
The following two keywords are optional. They are related to the PI controller used for the time step
control. The default values are given as in [1] for first and second order scheme.
Predictor-Corrector Control Beta 1 Real
Predictor-Corrector Control Beta 2 Real
Bibliography
[1] Gong Cheng, Per Lötstedt, and Lina von Sydow. Accurate and stable time stepping in ice sheet modeling.
Journal of Computational Physics, 329:29–47, 2017.
CSC – IT Center for Science
Chapter 7
Solving eigenvalue problems
7.1 Introduction
Eigenvalue problems form an important class of numerical problems, especially in the field of structural
analysis. Also some other application fields may have eigenvalue problems, such as those in density func-
tional theory. This manual, however, introduces eigenvalue computation in Elmer using terminology from
elasticity.
Several different eigenvalue problems can be formulated in elasticity. These include the “standard”
generalized eigenvalue problems, problems with geometric stiffness or with damping, as well as stability
(buckling) analysis. All of the aforementioned problems can be solved with Elmer. The eigenproblems can
be solved using direct, iterative or multigrid solution methods.
7.2 Theory
The steady-state equation for elastic deformation of solids may be written as
− ∇ · τ=~
f, (7.1)
where τis the stress tensor. When considering eigen frequency analysis, the force term ~
fis replaced by the
inertia term,
− ∇ · τ=−ρ∂2~
d
∂t2,(7.2)
where ρis the density.
The displacement can now be assumed to oscillate harmonically with the eigen frequency ωin a form
defined by the eigenvector ~
d. Inserting this into the above equation yields
− ∇ · τ(~
d) = ω2ρ~
d, (7.3)
or in discretized form
Ku =ω2Mu, (7.4)
where Kis the stiffness matrix, Mis the mass matrix, and uis a vector containing the values of ~
dat
discretization points. The equation (7.4) is called the generalized eigenproblem.
Including the effects of pre-stresses into the eigenproblem is quite straightforward. Let us assume that
there is a given tension field σin the solid. The tension is included by an extra term in the steady-state
equation
− ∇ · τ− ∇ · (σ∇u) = ~
f. (7.5)
The pre-stress term includes a component KGto the stiffness matrix of the problem and thus the eigenvalue
equation including pre-stresses is
(K+KG)u=ω2Mu. (7.6)
CSC – IT Center for Science
7. Solving eigenvalue problems 42
The pre-stress in Elmer may be a known pre-tension, due to external loading or due to thermal stress,
for example. The stress tensor containing the pre-stresses σis first computed by a steady-state analysis and
after that the eigenvalue problem is solved. It should be noted though that the eigenvalue problem in a pre-
stressed state is solved using first order linearization, which means that the eigenvalues are solved about the
non-displaced state. If the pre-loading influences large deformations the eigenvalues are not accurate.
The eigenvalue problem with pre-stresses may be used to study the stability of the system. Some initial
loading is defined and a pre-stress tensor σis computed. This tensor is then multiplied by a test scalar λ. The
critical load for stability, or buckling, is found by setting the force on the right hand side of the equation (7.5)
equal to zero. The problem then is to solve λfrom
Ku =−λKGu, (7.7)
which again is formally an eigenvalue problem for the test parameter. The critical loading is found by
multiplying the given test load with the value of λ. In other words, if λ > 1the loading is unstable.
7.2.1 Damped eigenvalue problem
Finally, let us consider the damped eigenproblem, also called quadratic eigenvalue problem. In this case
there is a force component proportional to the first time derivative of the displacement in addition to the
inertia term
− ∇ · τ=−δ∂~
d
∂t −ρ∂2~
d
∂t2,(7.8)
where δis a damping coefficient. The problem is transformed into a more suitable form for numerical
solution by using a new variable ~v′defined as ~v′=∂~
d
∂t . This yields
− ∇ · τ=−δ~v′−ρ∂~v′
∂t .(7.9)
Working out the time derivatives and moving into the matrix form, the equation may be written as
Ku =−Dv −iωM v, (7.10)
or,
iω I0
0Mu
v=0I
−K−Du
v,(7.11)
where iis the imaginary unit, Dis the damping matrix, and va vector containing the values of ~v′at the
discretization points. Now the damped eigenproblem is transformed into a generalized eigenproblem for
complex eigenvalues.
Finally, to improve the numerical behavior of the damped eigenproblem, a scaling constant sis intro-
duced into the definition s~v′=s∂~
d
∂t . In the matrix equation (7.11) this influences only the identity matrix
blocks Ito be replaced by sI. Good results for numerical calculations are found when
s=||M||∞= max |Mi,j |.(7.12)
7.3 Keywords related to eigenvalue problems
An eigenvalue analysis in Elmer is set up just as the corresponding steady-state elasticity analysis. An
eigenvalue analysis is then defined by a few additional keywords in the Solver section of the sif file. The
solver in question can be linear elasticity solver called Stress Analysis, linear plate elasticity solver, or even
nonlinear elasticity solver, though the eigen analysis is, of course, linear.
Many of the standard equation solver keywords affect also the eigen analysis, e.g. the values given for
Linear System Solver and Linear System Iterative Method in case of an iterative solver. More information
about these settings is given in this Manual under the chapter concerning linear system solvers. The specific
keywords for eigen analysis are listed below
CSC – IT Center for Science
7. Solving eigenvalue problems 43
Eigen Analysis Logical
Instructs Elmer to use eigensystem solvers. Must be set to True in all eigenvalue problems.
Eigen System Values Integer
Determines the number of eigen values and eigen vectors computed.
Eigen System Select String
This keyword allows the user to select, which eigenvalues are computed. The allowable choices are
•Smallest Magnitude
•Largest Magnitude
•Smallest Real Part
•Largest Real Part
•Smallest Imag Part
•Largest Imag Part
Smallest magnitude is the default.
Eigen System Shift Real
Gives an offset in the eigenvalues. May be useful if the user is interested in eigenvalues in a specific
pre-defined frequency-regime, for example. In practice the shift ω2
0is achieved by defining K:=
K−ω2
0Mand solving the eigenvalues of this modified system.
Eigen System Convergence Tolerance Real
The convergence tolerance for iterative eigensystem solver. The default is 100 times Linear System
Convergence Tolerance.
Eigen System Max Iterations Integer
The number of iterations for iterative eigensystem solver. The default is 300.
Eigen System Complex Logical
Should be given value True if the eigen system is complex, i.e. the system matrices are complex. Not
to be given in damped eigen value analysis.
Geometric Stiffness Logical
Defines geometric stiffness (pre-stress) to be taken into account in eigen analysis. This feature is only
available with linear bulk elasticity solver.
Stability Analysis Logical
Defines stability analysis. This feature is only available with linear bulk elasticity solver.
Eigen System Damped Logical
Defines a damped eigen analysis. Damped eigen analysis is available only when using iterative solver.
Eigen System Use Identity Logical
If True defines the relation displacement and its derivative to be s~v′=s∂~
d
∂t . The other possibility is to
use Mv =iωM u. The default is True.
Eigen System Lanczos Vectors Integer
Number of Lanczos vectors that are used to compute the eigenvalues. The default is 3Ne+ 1, where
Neis the number of eigenvalues.
Eigen System Compute Residuals Logical
Computes the residuals of the eigen value system.
CSC – IT Center for Science
7. Solving eigenvalue problems 44
7.4 Constructing matrices M and D in Solver code
In eigen analysis the mass matrix Mand the damping matrix Dhave to be separately constructed. Usually
in Elmer the different matrices are summed into a single matrix structure, since the final linear equation is
of the form Ax =b, and there is no need for separate values of the mass matrix and the stiffness matrix.
The matrix is represented in Elmer using compressed row storage (CRS) format, as explained in chapter
about Linear system solvers. The matrix structure holds also vectors for the values of the mass and damping
matrices
TYPE Matrix_t
...
REAL(KIND=dp), POINTER :: MassValues(:), DampValues(:)
...
END TYPE Matrix_t
These arrays use the same Rows and Cols tables than the normal Values array.
The mass and damping matrices are constructed elementwise in a similar manner as the stiffness matrix.
After each element the local contributions are updated to the equation matrices by the following subroutine
calls
CALL DefaultUpdateEquations( STIFF, FORCE )
IF ( Solver % NOFEigenValues > 0 ) THEN
CALL DefaultUpdateMass( MASS )
CALL DefaultUpdateDamp( DAMP )
END IF
In this segment of code the variables STIFF,MASS,DAMP and FORCE store the local values of the
stiffness matrix, the mass matrix, the damping matrix, and the right hand side of the equation, respectively.
The integer NOFEigenValues in the Solver data structure gives the number of eigenvalues requested.
Here it is used as an indicator of whether the mass and damping matrices need to be constructed.
The eigenvalues and eigenvectors are stored in the arrays Solver % Variable % EigenValues
and Solver % Variable % EigenVectors,
TYPE Variable_t
...
COMPLEX(KIND=dp), POINTER :: EigenValues(:)
COMPLEX(KIND=dp), POINTER :: EigenVectors(:,:)
...
END TYPE Matrix_t
and Solver % Variable % EigenVectors(i,:) gives the eigenvector corresponding to the eigen-
value i.
CSC – IT Center for Science
Chapter 8
Dirichlet boundary conditions and
limiters
Except for special cases, basically any partial differential equation (PDE) model needs boundary con-
ditions. Often their definition is solver-specific in Elmer. This chapter focuses on constraints which are
referred to as the Dirichlet conditions.
8.1 Dirichlet constraints
In the context of variational formulations of PDE models different kinds of boundary conditions are met. In
connection with the second-order PDE models, boundary conditions of Dirichlet type refer to cases where
the value of the field is set on the boundary. In engineering literature these are also knows as the essential
boundary conditions. The other types of boundary conditions involve derivatives of the solution and often
correspond to prescribing a flux or force on the boundary. They are referred to as the Neumann bound-
ary conditions or natural boundary conditions. More general boundary conditions of Robin type are also
possible. For problem-specific examples of them the reader is directed to the Models Manual.
Technically the Dirichlet conditions in ElmerSolver are set through manipulating only the values of the
linear system entries rather than the structure of the linear system. To be more specific, in setting the degree
of freedom with index i, the ith row of the matrix is set to be zero, except for the diagonal entry which is set
to be unity. When also the corresponding entry of the right-hand side (RHS) vector is set to have the desired
value, the solution will satisfy the Dirichlet condition.
A handicap of the procedure described is that the symmetry of the original matrix is not maintained.
This may affect the applicability of linear system solvers. There are two remedies to ensure that the matrix
remains symmetric. The column may also be zeroed except for the diagonal entry and the contribution of
the known values may be subtracted from the RHS vector. The second option is to eliminate all the rows
and columns related to the known values. This reduces the size of the matrix but has an additional cost as a
secondary matrix is created and the values are copied into it.
The standard Dirichlet conditions are given by referring to objects which have a lower dimension than
the leading dimension in the geometry definition, i.e. for 3D problems values are usually fixed only at 2D
faces. The Dirichlet conditions may hence be set for the degrees of freedom of existing boundary elements.
Additionally Dirichlet-like conditions may be set for a set of nodes that is created on-the-fly. It is possible
also to set Dirichlet-like conditions for bodies also. This may be particularly useful when the constraint is
only conditional.
Sometimes a Dirichlet condition should depend on other variables in a way which defines whether or not
to apply the condition at all. For example, the temperature at a boundary might only be defined in the case of
inflow, while for outflow giving the temperature would not be physically justified. For this kind of purposes
the user may define an additional real-valued variable that depends on the variable itself. Then, the Dirichlet
condition can be applied conditionally in cases where the additional variable has a positive value.
Simulation
CSC – IT Center for Science
8. Dirichlet boundary conditions and limiters 46
Set Dirichlet BCs by BC Numbering Logical
By default the Dirichlet conditions are set in the order that they just happen to appear in the
boundary element list. This is fine when there are no conflicting boundary conditions. However,
when such is the case the user may want to control the order of the conditions – the last Dirichlet
condition always prevails. This flag enforces the BCs to be set in the order of boundary condition
numbering in the sif file.
Boundary Condition bc id
Target Boundaries(n) Integer
The set of boundaries on which the Dirichlet conditions will be applied.
Target Nodes(n) Integer
Sets pointwise conditions on-the-fly. These points are referred via the absolute indexing of the
nodes.
Target Coordinates(n,DIM) Real
Sets pointwise conditions via transforming the coordinate values into nodal indeces correspond-
ing to the nearest nodes at the time of the first call. Target groups defined by Target Boundaries,
Target Nodes, and Target Coordinates should not reside in the same boundary con-
dition definition.
Varname Real
The value of each variable which has an active equation associated may be set by giving its value
in a boundary condition section. If the variable name is not listed in the keyword listing, the user
should also define the value type which is always Real.
Varname i Real
For multicomponent fields the Dirichlet condition may be set for each field component sepa-
rately.
Varname Condition Real
The Dirichlet condition related to the variable is set to be active only when the condition variable
is positive.
Varname {e} i Real
This keyword is used to define the vector gso that its tangential trace g×nis approximated by
uh×n, with uhthe edge element interpolating function corresponding to the solver variable.
The value of this keyword defines the components gi,i∈ {1,2,3}, with respect to the global
Cartesian coordinate system. For the edge finite elements see the appendix F.
Varname {e} Real
If the edge finite element has only degrees of freedom which are associated with element edges,
this keyword can be used to define a scalar field ˆgτso that the values of the prescribed degrees
of freedom are computed as (cf. the appendix F)
D=Z
Eij
[g(x)·τ(x) + ˆgτ]φ(x)ds,
with Eij an element edge. This keyword should not be used if additional degrees of freedom
associated with the element faces are present (that is, this feature should not be used for edge
elements of the second degree or edge elements from the optimal family).
The Dirichlet-like conditions for bodies. It is also possible to set the values of an exported variable with
exactly the same logic.
Body Force body force id
Varname Real
The setting of Dirichlet-like conditions for the whole body follows the same logic as for the
boundaries. When the body force definition is assigned to a body, the variable values will be
fixed as defined.
CSC – IT Center for Science
8. Dirichlet boundary conditions and limiters 47
Varname Condition Real
The conditional Dirichlet-like condition may also be given for bodies.
The following may be used in conjunction with the Dirichlet conditions to affect the setup of the linear
system equation.
Solver solver id
Linear System Symmetric Logical True
Make the matrix symmetric by eliminating the known values from the RHS and by zeroing the
appropriate matrix entries.
Before Linsolve "EliminateDirichlet" "EliminateDirichlet"
Creates the secondary matrix of a reduced size by eliminating Dirichlet conditions and passing
this to the linear system solver.
8.2 Soft limiters
The user may set soft lower and upper limits to the values of the field variable. For example, concentration
can never be negative and hence a zero lower limit could be applied for it. The limits are applied in an
iterative manner defining a contact set where Dirichlet-like conditions are applied. A node is included in the
contact set when its value passes the limits, and a node is released from the contact set when the nodal load
related to it is negative (or positive). The limiters may be applied to both boundary conditions and bodies.
Solver solver id
Apply Limiter Logical True
Activates the search for limiters and also activates the computation of nodal loads.
Limiter Value Tolerance Real
Defines the tolerance of the field value by which a node is added to the contact set.
Limiter Load Tolerance Real
Defines the tolerance of the nodal load to determine whether a node is removed from the contact
set.
Limiter Load Sign Negative Logical
The default sign used for determining the contact set is derived from the default discretization of
the Poisson equation. If the equation would be multiplied by −1, the load would also change its
sign. With that in mind, a possibility to influence the rules is given with this keyword.
Boundary Condition bc id
Varname Lower Limit Real
The lower limit of field variable.
Varname Upper Limit Real
The upper limit of field variable.
Body Force body force id
Varname Lower Limit Real
The lower limit of field variable.
Varname Upper Limit Real
The upper limit of field variable.
CSC – IT Center for Science
Chapter 9
Enforcing shared conditions between
boundaries
Special classes of boundary conditions arise from treating constraints that involve two boundaries si-
multaneously. For example enforcing the periodicity of the solution is considered here. We also consider
conditions relating to the application of the mortar finite element method, which offers a general strategy to
approximate the solution continuity across a shared boundary when the meshes of two independently dis-
cretizated bodies are not conforming on the shared boundary (that is, the locations of the mesh nodes do not
coincide on the shared boundary).
Historically the implementation of periodic boundary conditions in Elmer has been based on projec-
tors which are derived directly using the interpolation scheme. However, over the time the development
has evolved towards mortar finite element strategies which lead to projectors based on weak formulations.
Currently both approaches are still used but the mortar finite element approach should be employed when
striving for the best possible accuracy over nonconforming meshes.
9.1 Periodic conditions
Periodic BCs may be considered to be a special case of Dirichlet conditions where the fixed value is given as
a linear combination of other unknown values. The periodic boundary conditions in Elmer are very flexible.
In fact they may even be antiperiodic.
Boundary Condition bc id
Periodic BC Integer
Periodic boundaries generally come in pairs. This command defines a pointer to another bound-
ary condition definition which is considered to be the counterpart.
Periodic BC Explicit Logical
Sometimes the implicit periodic BCs (the default) leads to convergence problems. It also com-
plicates the matrix structure by adding additional connections. Then the explicit type of periodic
conditions may be a good alternative. Note that it requires iteration even for a linear operator.
Periodic BC Translate(3) Real
The correspondence of the periodic boundaries should be obtained by performing three different
operations: translation, rotation and scaling. This generality is not usually needed and therefore
default values are usually appropriate. For translation the default translation vector is obtained
via moving in the normal direction of the primary boundary until the target boundary is hit. If
this is not desired, the user may give another translation vector using this keyword.
Periodic BC Rotate(3) Real
By default no rotation is performed prior to the construction of the constraints. This keyword
may be used to give the angles of rotation. There is no automatic rotation, but there is some
detection of the correctness of the rotation.
CSC – IT Center for Science
9. Enforcing shared conditions between boundaries 49
Periodic BC Scale(3) Real
This keyword may be used to give a scaling vector if this is desired. If this is not given, an
isotropic scaling will be used to map the bounding boxes of the two boundaries to the same size.
Periodic BC VarName Logical
The user should define the variables that are to be periodic in nature. This is done by attaching
their names into logical expressions following the string Periodic BC.
Anti Periodic BC VarName Logical
The variable may be also antiperiodic, i.e. the absolute value is the same but the sign is different.
Then this should be used instead.
Periodic BC Use Lagrange Coefficients Logical
By default the periodic BCs are handled by elimination. The other alternative is to use a Lagrange
multiplier. This strategy is enforced by setting this keyword to be True.
9.2 Mortar conditions
The motivation for implementing mortar finite element framework into Elmer came from the simulation of
rotating electrical machines. There the continuity of solution between rotor and stator must be approxi-
mated. Both parts usually have their own fixed meshes that move with respect to each other. The principal
complication associated with the moving mesh approach is that the finite element meshes are then generally
non-matching across the model boundary which represents the interaction surface Γshared by the rotating
and surrounding bodies. The usual continuity requirements to obtain a conforming finite element approxi-
mation over the entire domain cannot thus be satisfied. To this aim a monolithic discretization strategy based
on the mortar finite element method has been implemented.
The implementation of the mortar conditions in Elmer is related to the periodic boundary conditions in
that some common routines and techniques are shared. The mortar conditions may be dynamic, whereas the
periodic conditions are assumed to be fixed.
The mortar finite element formulations approximate the solution continuity via a weak formulation and,
in the case of the scalar unknown, lead to evaluating integrals of the type
Z(u1−u2)φSdΓS= 0 (9.1)
where the integration is over the element(s) on the slave (non-mortar) side, u1the value of the matched
quantity at the slave side and u2the quantity at the master (mortar) side. In addition, φSis a test function.
The integration is accomplished in practice by searching the locations of the integration points within the
mortar elements. If the set of the test functions coincides with the trace space of admissible solutions on the
non-mortar side, the degrees of freedom at each side are related via the linear equations
Mu1=Nu2,(9.2)
with
Mij =Rφj
Sφi
SdΓS,
Nij =Rφj
Mφi
SdΓS.(9.3)
All the integration points of the non-mortar side must be locatable in the set of mortar elements. There-
fore the mortar boundary may be larger than the non-mortar boundary but never vice versa. To enforce the
constraint in the original system a constraint matrix block
C= [M−N](9.4)
is constructed. This constraint block is then added to the original system by using the Lagrange multipliers
as A CT
C0u
λ=f
0.(9.5)
CSC – IT Center for Science
9. Enforcing shared conditions between boundaries 50
Re-parametrization of the shared boundary
When we are studying problems where the potential contact surface (or line) is known, a re-parametrization
of the surface meshes may be useful. For example, in the case of a cylindrical surface the element faces
may be mapped from (x, y, z)coordinates to (θ, z)coordinates since the r-coordinate should be the same
by construction. Finding the non-mortar points from the mortar side becomes then somewhat easier as
we avoid the problem of two mismatching polygons. Similarly, (r, z)coordinates may be used instead of
(x, y, z)coordinates when the solutions on two boundaries which are planes in the (r, z)coordinates are
matched. The similar strategy works also for constraints that need to be expressed at certain coordinate
level. For example, we may use the (x, y)-plane to express the constraints instead of the full 3D coordinates.
When using the reduction to (θ, z)plane, we may also consider periodicity of the projector. Assuming
that the mortar nodes are in interval [θmin, θmax]such that 2π/(θmax −θmin)is an integer we may always
map the angle φto the existing interval from anywhere in [0,2π]. Similar repeated projector could basically
also be taking place in translation between two planes if their structure would be periodic.
Integration of mortar projector
When performing the integration, there are several ways how the integration and search is carried out. Which
one is optimal depends greatly on how the mesh is configured.
If the mesh is basically conforming and the need for the projector is caused by the internally created
discontinuity (using keyword Discontinuous Boundary), we have the luxury that there is an exact
correspondence between a master and a slave element. Then an exact projector may be computed with a
minimal effort. By default the projector is of mass matrix type but also a diagonal projector may be created
by mass lumping. This branch of projectors is invoked if the corresponding Mortar BC is given the value
zero. This approach works also in parallel, at least if the boundary halo is saved during the partitioning
process.
In some cases of rotating machines created by extrusion we have the luxury that the nodes are only at
specific coordinate levels. This makes it easier to have an accurate integration over the intersection of the
elements. In the current implementation of Level Projector it is assumed that the surface elements
are quadrilaterals and two of the element sides remain parallel. The integration is performed by checking
whether any of master and slave elements cross each other and creating a temporal triangulation for their
intersection. This is a loop of complexity O(N2)but for typical cases it is still not a bottle-neck. This
approach may work in parallel when using a tailored partitioning algorithm that distributes the elements
such that they by construction will be in the same partitioning as their counterparts.
A generic approach is used if we have no simplifications that could be utilized when integrating over
the non-mortar elements. The integration utilizes oct-tree based search algorithm. Unfortunately the inte-
gration is not accurate (a temporal triangulation for improved integration accuracy is not generated in this
case) which may result to accuracy and convergence problems when treating the resulting constrained linear
system. Currently there is no way to perform parallel communication when performing the integration.
Biorthogonal basis for mortar projectors
There is a possibility to use biorthogonal basis functions when creating the mortar projector with some
integrators. The biorthogonal basis is advantageous since either the non-mortar or mortar degrees of freedom
may then be eliminated from the linear system, thereby significantly simplifying the solution.
Contact mechanics
ElmerSolver library includes some generic features for contact mechanics. The functionality relies heavily
on contact mortar projectors which are used to ensure the match on the contact surface. The contact pressure,
on the other hand, is used to determine the active contact set quite similarly as for the case of soft limiters.
The contact features are applicable to all elasticity equations having the displacement degrees of freedom
at each node. The contact mechanics is still under development and the documentation is therefore not
complete. There are severe limitations when it comes to parallel operation and conflicting contact conditions
on different boundaries.
CSC – IT Center for Science
9. Enforcing shared conditions between boundaries 51
Keywords for mortar projector
Solver solver id
Apply Mortar BCs Logical
Enforces the application of mortar BCs. This is not by default set to be True partly because at
the time of writing this machinery is still under development. If this flag is set to be true, then
the mortar BCs are applied to the variable of this solver. Otherwise the mortar projectors are not
applied to the solver variable.
Edge Basis Logical
Activates the use of edge element basis functions for the mortar conditions.
Boundary Condition bc id
Mortar BC Integer
This is set for the non-mortar (slave) boundary to refer to the boundary condition index of the
corresponding mortar boundary. The mortar boundaries always come in pairs, but the settings
are only given in the context of the non-mortar boundary.
Rotational Projector Logical
The interface lying in 3D space may be flattened to 2D plane by several ways. Using the flattened
geometry results to better interpolation accuracy since faceting is not needed. This keyword
activates the use of the (θ, z)coordinate system when computing the projector. The angle may
have a discontinuity around 180 degrees. The Level Projector option offers a special way
to deal with cases that make a full circle.
Cylindrical Projector Logical
Similar to Rotational Projector except there is automatic fitting to get a cylinder surface
such that the cylinder axis is determined via computation. Also, if the angle is not full 2π, it will
not be enforced.
Anti Rotational Projector Logical
As the previous one but if the module or the periodic cycle is odd use a negative coefficient in
enforcing the periodicity. This way the mortar projector becomes antiperiodic.
Sliding Projector Logical
Perform dimensional reduction so that only two of the coordinates are accounted for. By default
this means mapping to (x, y)plane.
Anti Sliding Projector Logical
As the sliding projector except if the mapping of the sliding projector cycle is odd use a negative
coefficient in the mortar projector.
Radial Projector Logical
Use the (r, z)coordinate system when computing the projector.
Mortar BC Scaling Real
The coefficient with which the values are scaled when enforcing continuity over the mortar
interface. The default is one.
Level Projector Logical
Level projector assumes that the re-parametrized surfaces consist of quadrilaterals with two par-
allel edges at constant z-levels. This projector may use either strong (interpolation based) or
weak strategies to find the constraints for edge and nodal degrees of freedom. There are many
additional keywords that apply to the level projector.
Plane Projector Logical
This is similar to Level Projector except the plane over which the projection will be made
is fitted automatically such that the best fit of a plane associated with either slave or master
boundary will be used.
Level Projector Nodes Strong Logical
This uses strong (interpolation based) projector for the nodal degrees of freedom. The default is
False, i.e. a weak projector is used.
CSC – IT Center for Science
9. Enforcing shared conditions between boundaries 52
Level Projector Plane Edges Strong Logical
For extruded geometry some edges will always be aligned with the plane. For these edges it
seems natural to use only the other edges on the same constant z-plane to create the projector.
Then this keyword should be set True. The default is False.
Level Projector Extruded Edges Strong Logical
The extruded edges are the ones that are created when the 2D mesh is extruded into the third di-
rection. So they are aligned with the z-axis. This keyword controls how the edge element degrees
of freedom are treated in the z-direction. Using the default value False is recommended.
Projector Skip Nodes Logical
The level projector may skip to create the projector for the nodal degrees of freedom. This would
also be the choice when there are no nodal degrees of freedom for the solver. Also there could
be a physical reason that would make them obsolete.
Projector Skip Edges Logical
Also the edge element degrees of freedom may be skipped. They will be treated by default if
there are degrees of freedom associated with the edges.
Save Projector Logical
Save the projector values to an external file for mainly debugging purposes. Default is False.
The following keywords are related to the use of biorthogonal basis.
Biorthogonal Basis Logical [False]
This will enable the elimination of slave degrees of freedom by making the constraint matrix
block corresponding to the slave entries diagonal.
Biorthogonal Dual Master Logical [True]
This option also makes the constraint matrix block corresponding to the slave entries diagonal.
The value of this keyword defines whether the nodal basis or the biorthogonal basis should be
used on the master side. The default is the latter.
Biorthogonal Dual LCoeff Logical [False]
Here nodal basis functions are used for expressing the traces of the fields to be matched, while
the Lagrange multiplier is treated with the biorthogonal basis functions.
The following keywords are related to the use of contact mechanics.
Solver solver id
Apply Contact BCs Logical
Enforces the solver to apply contact mechanics conditions for the boundaries where the flag is
given. When doing contact mechanics, it may be useful to set Linear System Residual
Mode = Logical True so that only the difference is solved making the solution process
more robust.
Boundary Condition bc id
Contact Depth Offset Initial Real
Sometimes the contact set is initially empty which may result to rigid body motion making the
system indefinite. Therefore it may advantageous to introduce some initial offset to the contact
set detection such that the initial contact set is increased for the first iteration.
Contact Depth Offset Real
This has the same effect as the previous one but is used for each iteration. This could be used to
model a tight fit of cylinders, for example.
Contact Active Condition Real
This keyword may be used to enforce active contact somewhere. For positive values the contact
is active, otherwise the flag has no effect.
CSC – IT Center for Science
9. Enforcing shared conditions between boundaries 53
Contact Passive Condition Real
This keyword may be used to enforce no-contact somewhere. For positive values the contact is
passive, otherwise the flag has no effect.
Stick Active Condition Real
This keyword may be used to enforce active stick condition somewhere. For positive values the
contact is active, otherwise the flag has no effect.
Stick Passive Condition Real
This keyword may be used to enforce slip contact somewhere. For positive values the slip con-
tact, otherwise the flag has no effect.
Static Friction Coefficient Real
This keyword is used to give the static friction coefficient on a surface. The magnitude of the
tangential force will then be the friction coefficient times the normal contact force.
Dynamic Friction Coefficient Real
This keyword is used to give the dynamic friction coefficient on a surface.
Stick Contact Coefficient Real
This keyword is used to give the stick coefficient that allows a small tangential movement even
if stick model is used.
Contact Velocity(dim) Real
For steady state contact the velocity cannot be retrieved from history. Then the user may give the
relative velocity between the slave and master surfaces. The velocity magnitude is not important
for linear friction models but the direction of the velocity is used to set the direction of the friction
forces.
Normal-Tangential VarName Logical
If the normal of the surface is not suitably aligned with some coordinate direction, it is rec-
ommended to use a normal-tangential coordinate system except for the tie contact. Here the
Varname would typically be Displacement.
Contact Active Set Minimum Integer
If the contact set is empty, there may be numerical problems. Therefore a minimum size for the
contact set may be defined. Then the requested number of nodes with the shortest distance are
added to the contact set.
Tie Contact Logical
Contact type where the contact set is active on the whole surface. Hence there is no testing
performed on the direction of the contact pressure.
Stick Contact Logical
Contact type where the tangential contact set is the same as the normal contact set.
Slide Contact Logical
Contact type where the tangential contact set is empty.
Friction Contact Logical
Contact type where the tangential contact set is tuned after some conditions. This flag is not fully
operational yet.
Contact Type String [stick, tie, friction, slide]
Instead of using the above flags to set the contact type, the user may also give a string that
chooses between the different options.
CSC – IT Center for Science
Chapter 10
Linear systems with constraints
10.1 Introduction
There are a number of cases where a set of linear equations under additional constraints has to be solved.
For example, in conjunction with the mortar finite element methods the application of the mortar projection
introduces additional constraints that are dealt with introducing Lagrange multipliers.
There are two different ways to treat the linear constraints. The first one is to use a large monolithic
system where new constraint equations are added. This is generally applicable but may result in numerically
difficult problems. The other approach uses elimination to reduce the size of the system. It is numerically
favorable method but has some limitations related to its generality.
10.1.1 Solving linear systems with constraints
The Lagrange multiplier procedure leads to a matrix that is indefinite. Any Krylov method capable of
handling indefinite matrices should be applicable to the constrained system. In practice the GCR and
BiCGStab(L) methods in Elmer has been used succesfully in applications. The constraint equations can
be effectively preconditioned using ILU (for example).
In the case of (9.5) the current implementation in Elmer uses a preconditioning matrix
P=P0CT
C0(10.1)
For the ungauged edge element discretization of the A-V formulation we select
P0=I(10.2)
for magnetostatic cases, and
P0=I0
0AV(10.3)
for dynamic cases, with AVthe diagonal block corresponding to the scalar potential of the A-V formulation.
As we are not using any gauging for the curl-curl systems currently, the direct solvers are not available for
solution of these systems.
Currently the mortar projector assumes that each node on the master boundary can be found in the slave
boundary of the same partition. This limits the current parallelization strategy such that the ownership of
the mortar boundaries must be carefully considered when making the partitioning. If the ownership issues is
honoured, standard iterative parallel linear solvers may be used for solving the constrained system.
CSC – IT Center for Science
10. Linear systems with constraints 55
10.1.2 Eliminating linear constraints
Sometimes the constraint block Callows for applying an elimination procedure. The elimination may be
highly beneficial since it removes the saddle-point nature of the problem and makes it much easier to solve
with standard linear solvers.
10.2 Keywords related to treatment of linear constraints
The following keywords should be put in the solver section of the solver in question.
Apply Mortar BCs Logical
This flag introduces mortar BCs that results in linear constraints. More information on mortar BCs are
provided in Section 9.2.
Apply Contact BCs Logical
This flag introduces contact BCs that results in linear constraints.
Export Lagrange Multiplier Logical
If linear constraints are involved, this defines whether the related Lagrange multipliers should be saved
as a field. Often this field has a useful physical interpretation.
Lagrange Multiplier Name String
Name of the Lagrange multiplier to be exported.
Use Transpose values Logical
Should we use the constraint matrix in a transposed form.
Nonlinear System Convergence Without Constraints Logical
Taking into account the linear constraints (i.e. the Lagrange multipliers) in the nonlinear system norm
may have a large impact on the value. Often it is preferable just to study the norm of the solution and
neglect the additional dofs associated with the constraints. This is enforced by setting this keyword to
True. The same effect is obtained also by setting NonLinear System Consistent Norm to
True.
Eliminate Linear Constraints Logical
This specifies whether the Lagrange multipliers are eliminated. A sufficient condition for the elimina-
tion is that the slave part of the constraint matrix is diagonal. This may be obtained by either using a
biorthogonal basis with the accurate integration or in conjunction with a discontinuous boundary.
Eliminate Slave Logical
This enforces also the elimination of slave degrees of freedom using the constraint equations. This is
not possible if the constraint matrix is not diagonal (for example when using Biorthogonal Dual
LCoeff = True), or there are diagonal entries in the constraint equation rows.
Eliminate From Master Logical
Instead of eliminating the Lagrange multiplier using the slave part use the master part. This could be
applicable if we use the keywords Discontinuous Boundary and Mortar BC Diag and the
constraint equation involves the diagonal. In this way both the Lagrange multiplier and the slave dofs
can be eliminated.
CSC – IT Center for Science
Chapter 11
Solving constraint modes
11.1 Introduction
Elmer includes a possibility to compute constraint modes associated with a selected part of the boundary
referred to as the attachment boundary. These are solution modes that are obtained by successively applying
a unit field value to each degree of freedom associated with the attachment boundary. The results may be
used for model lumping, for example. Together with solving eigenproblems where the degrees of freedom
associated with the attachment boundary are fixed, the constraint modes may be employed to perform the
Craig-Bampton component mode synthesis.
11.2 Keywords related to constraint modes analysis
Solver solver id
Constraint Modes Analysis Logical
When true, this instructs Elmer to compute constrained modes using the field equation of the
solver.
Constraint Modes Lumped Logical
This defines whether the constraint modes are lumped such that all the dofs associated with the
attachment boundary are given the fixed value, or whether all the dofs have a separate permu-
tation. The default is that the system is not lumped and each dof has a separate permutation.
Setting this flag to True enforces lumping.
Constraint Modes Fluxes Logical
Computes fluxes between the constraint modes. The fluxes are used to determine a conductivity
matrix between the constraint modes. For example, if the equation in question is electrostatic
equation, the fluxes will be surface charges, and the resulting lumped matrix is the capacitance
matrix.
Constraint Modes Fluxes Symmetric Logical
This is used to assume that fluxes between constraint modes are symmetric in order to enforce
the symmetry of the conductivity matrix.
Boundary Condition bc id
Constrained Modes Varname String
This keyword can be used to specify the attachment boundary and to specify that the field variable
Varname is used for obtaining the constraint modes.
CSC – IT Center for Science
Chapter 12
Working with nodal loads
Usually source data for equations are specified as spatial fields. Elmer then internally transforms these
definitions to entries of the right-hand side vector by following the weak formulation of the problem. This
is convenient as the specification of data is independent of the underlying mesh.
However, it may sometimes be attractive to define discrete loads directly as entries of the right-hand
side vector. Similarly, after the solution of the discrete system has been accomplished, it may be useful
to consider discrete reactions which occur at nodes to maintain the balance. Both options are supported in
Elmer. In addition, integrals of the basis functions can be computed to generate weights which can then be
employed to reproduce a spatial data field approximately corresponding to the entries of the given right-hand
side vector.
12.1 Setting nodal loads
Similarly to the Dirichlet constraints one may also set nodal loads, i.e. the entries of the right-hand side
vector in the linear system. This has to be done with care as proper values generally depend on the mesh.
There are anyhow some uses for this option. For example, in multiphysical couplings it may sometimes be
a good choice to transfer the forces directly in the nodal form as this is the consistent way to compute the
forces resulting from the solution of the discrete system.
Body Force bf id
Varname Load Real
This can be used to specify directly the entries of the right-hand side vector in the linear system.
The nodal loads are given exactly as the Dirichlet conditions except that a string Load is attached
to the name of the variable occuring in the place of Varname.
12.2 Computing nodal loads
It is possible to evaluate the nodal loads after the solution has been computed. This however requires that
the original matrix A0is maintained without eliminating the rows and columns that correspond to the known
Dirichlet conditions. Then the nodal forces may be computed from
f=A0x−b. (12.1)
It should be noted that the nodal values depend on the mesh. If SI units are employed, for heat equation they
are expressed in Watts and for electrostatic equation in Coulombs, for example.
Solver solver id
Calculate Loads Logical True
This keyword activates the computation of nodal loads. The resulting values will be saved to
the variable the name of which is derived from the primary variable name by adding the suffix
Loads to it.
CSC – IT Center for Science
12. Working with nodal loads 58
12.3 Computing nodal weights
While the computation of the nodal loads often gives the best approximation of boundary fluxes, for example,
they are however in cumbersome units as the nodal loads depend on the mesh. It may thus be more ideal to
transform the nodal loads into a distributed field. To this aim there is a possibility to compute the values of
the integrals
wi=ZΩ
φidΩ(12.2)
where φiis the basis function associated with the mesh node i. When the basis forms a partition of unity, the
sum of all weights wigives the volume (or area) of the domain occupied by the elements. An approximate
reproduction of the data field ¯gcorresponding to the nodal load f may be obtained as ¯g=fi/wi.
Solver solver id
Calculate Weights Logical True
This keyword activates the computation of the nodal weights wi.
12.4 Energy norm
When the initial matrix is known, the energy norm of xmay be computed as
E=xTA0x. (12.3)
Solver solver id
Calculate Energy Norm Logical True
This activates the computation of the energy norm. The result will be written to the Simulation
block with name res: VarName Energy Norm so that the SaveScalars subroutine may be
used to save the result. The energy norm may only be computed when the loads are also com-
puted.
CSC – IT Center for Science
Chapter 13
Generic solver utilities
When the solvers use the default procedure for solving the partial differential equations, there are a
number of generic features that may be used with any model. This chapter describes these features.
13.1 Solver activation
There is a large number of different ways how solvers need to be activated and deactivated. Mostly these
needs are related to different kinds of multiphysical coupling schemes. There is also a large number of
auxiliary solvers that are needed, for example, only when results are saved or inspected. In the Solver
section one may give the following keywords.
Exec Solver String
The value options are here never, always, before all, before timestep, after
timestep, before saving, after saving, after all. If nothing else is specified,
the solver is called every time in its order of appearance. The saving instance refers to the one defined
by Output Intervals which controls the saving of the results.
Exec Intervals Integer
This keyword gives an interval at which the solver is active. At other intervals the solver is not used.
This should have the same length as Timestep Intervals in the Simulation section.
Start Time Real
The start time after which the execution of the solver is started.
Stop Time Real
The stop time after which the execution of the solver is halted.
Exec Condition Real
Given this keyword, activate the solver only upon a positive value for the keyword.
13.2 Variable names
The variable name is presented in the Solver section by keyword Variable, for example
Variable = Varname
This name is used when setting Dirichlet conditions and initial conditions. Also the name is used as a basis
for other features appending it with suffixes such as Load,Condition and Passive, described later in
this chapter.
Sometimes one may want to rename the components of the primary variable. This may be done by
defining the component names in the brackets. For example the command
CSC – IT Center for Science
13. Generic solver utilities 60
Variable = Flow[Velo:2 Pres:1]
decleares that variable Flow consists of Velo with two components and Pres with one component. If the
number of components is 2 or 3, the variable will be treated as a vector in the ElmerPost files.
If one does not require output for a certain variable, one may prevent it with the -nooutput option,
e.g.
Variable = -nooutput Dummy
One may also use the -dofs option to define the number of components in a variable, e.g.
Variable = -dofs 3 Flow
By default variables are assumed to be field variables. However, they may also be scalars which have
globally the same value. These variables may be introduced with the -global option, e.g.
Variable = -global Frequency
After defining a global variable, it may be used similarly to time in giving dependencies.
These different options should not be fully mixed.
13.3 Exported variables
The intent of exported variables is to enable automatic allocation and treatment of additional data that may
usually be derived from the primary fields. Often this is done within a solver and many times the machinery
is used transparently from the user.
Each solver that has a primary variable (defined by the Variable keyword) may also have exported
variables. When this is made active, the values of exported variables may be set either in Body Force
or in Boundary Condition Section. The exported variable inherits the permutation of the primary
variable. The values of the exported variable can hence be updated at the same nodes as the primary variable,
and the variable may have the same optional arguments.
Upon request the exported variables may also be defined by the user in the Body Force and Boundary
Condition sections. The operation is set behind keywords to circumvent unwanted definitions.
Solver solver id
Exported Variable i Varname
A name for an additional variable created by the solver, i=1,2,3,....
Update Exported Variables Logical
Update the exported variable values after leaving the iterative solution for the next solver, i.e. in
the steady-state iteration level.
Nonlinear Update Exported Variables Logical
Update the exported variables within the nonlinear solution to get the current iterate.
13.4 Active and passive elements
In Elmer it is possible to define certain areas of the modeled geometry to be passive during the solution. This
feature also allows for deactivating and reactivating of the elements. Being a passive element means that
its contribution is not included in the global matrix assembly. One could, for example, model two separate
bodies heated with different heating power and connect them with a third structure only after suitable time
has elapsed. This all could be modeled within a single simulation.
The geometry of the whole system is meshed as usual, and the passive elements are only omitted in the
assembly of equations. The passive definition is done solverwise and elementwise. The former means that
for example the temperature may be passive and the displacements active at the same element. The passive
property of elements is defined in the Body Force section using a real valued parameter with the name
constructed from the name of the variable followed by Passive. When the parameter obtains a value
greater than zero, the element is passive.
CSC – IT Center for Science
13. Generic solver utilities 61
Body Force body force id
Varname Passive Real
If this variable obtains a positive value, the element is set passive. Note that it is not possible to
control components of vector-valued variables separately.
13.5 Timing the solvers
Often it is of interest to time the performance of the different steps in the solution sequence. For that purpose
there are some keywords for activating timers for the solvers.
Solver solver id
Linear System Timing Logical True
This keyword activates the timing for the linear system.
Linear System Timing Cumulative Logical True
This keyword sums up the cumulative time used in solving the linear systems.
Solver Timing Logical True
This keyword activates the timing for the whole solver including iteration over nonlinear and
linear systems etc. The time used for assembly may roughly be estimated as the difference
between the linear solution time and the total time used in the solver.
Solver Timing Cumulative Logical True
As the previous one but sums up the cumulative time used in the solver.
13.6 Consistency testing of solvers
Each solver may be tested for consistency. The way the testing is performed is that it is assumed that the
correct solution will have a predefined norm. After the simulation has been carried out the computed norm
and predefined norm are compared. If they are the same with a given tolerance the test is assumed to be
passed.
The result of a test is outputted on the screen. Additionally a file named TEST.PASSED is saved in
the working directory. If all the tests of the case are passed then the number 1 is written to the file in ascii
format. If any of the tests fails then the number 0 is written to the file, correspondingly.
Reference Norm Real
The target norm for the solver. The computed norm is tested against this one.
Reference Norm Tolerance Real
The tolerance for testing. The default is 1.0e-6. A sloppier tolerance could be used if the case is known
to be sensitive againts numerical fluctuations.
CSC – IT Center for Science
Chapter 14
Meshing Utilities
ElmerSolver includes some internal possibilities to affect the mesh. These include the ability to split
the mesh repeatedly on a partitioned level (mesh multiplication), the possibility to extrude meshes within
the code, and the possibility to create discontinuous meshes on-the-fly. These features make it possible to
eliminate some meshing bottle-necks and to implement more versatile physical boundary conditions, for
example.
14.1 Coordinate transformation
By default Elmer read in coordinates ~r = (x, y, z)as they are in the mesh files. The user may perform
scaling such that ~r := C ~r. Here Cmay be either a scalar or a diagonal matrix that maps each coordinate
differently. It could be used to scale the units of the mesh back to metres if they were initially given in mm.
The coordinates (x, y, z)may be mapped to any order by a suitable permutation which is any ordered
list of numbers 1, 2 and 3. The mapping could be useful, for example, if the initial coordinate system has the
rotational axis at x-axis whereas Elmer always assumes that the rotational axis is the y-axis.
Coordinate transformations can be used to map the cartesian coordinates to some other coordinate system
and vice versa. Currently supported transformations are from cartesian to cylindrical,
r=r0+px2+y2(14.1)
φ= arctan(y, x)(14.2)
z=z(14.3)
and cylindrical to cartesian,
x= cos(φ)(r+r0)(14.4)
y= sin(φ)(r+r0)(14.5)
z=z. (14.6)
The following keywords are related to coordinate tranformations
Simulation
Coordinate Scaling Real
The coefficient used to scale the coordinate.
Coordinate Mapping(dim) Integer
The permutation of the coordinates.
Coordinate Transformation String
Options are "cartesian to cylindrical", "cylindrical to cartesian" and "previous".
Coordinate Transformation Radius Real
The offset r0in radius.
CSC – IT Center for Science
14. Meshing Utilities 63
Note that the Coordinate Transformation may also be given in the Solver section which can then
be used to invoke different coordinate systems for different solvers.
14.2 Mesh multiplication
Mesh multiplication is the process where each mesh edge is splitted into two resulting to an increased number
of elements. In 2D each mesh multiplication increases the number of elements by factor of 4, and in 3D by
factor of 8 correspondingly. For hexahedral elements this results to eight elements that are of original shape.
For tetrahedral elements the innermost octahedron may be split into four tetrahedra in three ways of which
only one preserves the minimum angle condition (Delaunay splitting) and this is naturally implemented in
Elmer.
The mesh multiplication may be used recursively. This makes it possible to create meshes with a very
large number of elements. The functionality of mesh multiplication may also be used in conjunction with the
geometric multigrid of Elmer. Saving the meshes of different density gives an optimal hierarchy of meshes.
The main handicap of mesh multiplication is that the geometric description will not be improved during the
mesh multiplication. Therefore curved boundaries will be as faceted as they were in the coarsest mesh level.
Elmer includes also a feature for maintaining mesh grading. For simple geometries the mesh grading
may be conserved with the following scheme:
1. Solve at the coarsest level with Galerkin method the nodal mesh parameter from the equation h=
V1/m
e, where his the element size factor, Veis the element volume and mis the dimension of grading.
2. When splitting an edge for the new mesh, inherit the size factor and coordinate from h′=√h1h2and
~r′= (h1~r2+h2~r1)/(h1+h2).
The procedure may be applied recursively. It provides ideal grading for structured boundary layers with
m= 1 and homogeneous Delaunay meshes with m=dim.
The functionality of mesh multiplication is implemented also in parallel and is very efficient since it
requires only quite little communication.
The following keywords are related to the mesh multiplication.
Simulation
Mesh Levels Integer
The number of mesh levels when using the equal split feature.
Mesh Keep Integer
The user may choose to optionally keep more than one level. This could be needed for coarser
postprocessing, for example.
Mesh Keep Grading Logical
When creating meshes using equal split, the elements are by default split to equally sized pieces.
When this flag is on, the solver tries to maintain the mesh grading.
Mesh Grading Power Real
The mesh grading is evaluted from the element sizes of the intial mesh. The size is a scalar
function and cannot therefore handle very complicated meshes. For boundary layer type of
meshes the optimal value is one, while for meshes of Delaunay type the optimal value is the
space dimension.
14.3 Mesh extrusion
The feature of mesh extrusion was initially implemented for the needs of computational glaciology where
the meshing had become a serious bottle-neck. The current implementation of mesh multiplication takes
a 1D or 2D mesh and extrudes it to the desired direction. In the extrusion process the dimension of the
elements also grow. Line segments are extruded into quadrilaterals, triangles into wedges, and quadrilaterals
into hexahedrons, correspondingly.
CSC – IT Center for Science
14. Meshing Utilities 64
When the mesh is extruded, the boundary condition will also be extruded. New boundary conditions
may be set to the bottom and top of the mesh. These boundaries will get the first free boundary condition
index and may thereby be referenced in the sif file.
The mesh density of the extrusion is defined only once and this same division is copied for all nodes. So
each node gives rise to exactly similar 1D mesh stride in the desired direction.
Because only one division needs to be defined it can be more costly than would be otherwise feasible.
Therefore also an iterative process for setting the division has been implemented. There the user can give
an arbitrary mesh parameter function h(x)and a 1D mesh is created so that the mesh position as accurately
obeys the desired distribution. The algorithm is an iterative one based on a Gauss-Seidel type of iteration.
First initialize the distribution:
do i= 0, n
wi=i/n
end
Then iteratate the distribution until it has converged:
do while( max |dwi|> ǫ )
do i= 1, n
hi=f((wi+wi−1)/2)
end
do i= 1, n −1
wi= (wi−1hi+1 +wi+1hi)/(hi+hi+1)
end do
do i=n−1,1,−1
wi= (wi−1hi+1 +wi+1hi)/(hi+hi+1)
end do
end do
The following keywords are used to control the internal extrusion
Simulation
Extruded Mesh Levels Integer
The number of node layers in the extruded mesh. This keyword not only defines the number but
also activates the whole feature. The number of node layers is the same over the whole extruded
mesh.
Extruded Coordinate Index Integer
The index (1,2,3) of the coordinate to be extruded. The default direction is dim+1. So for 2D
mesh the extrusion is by default carried out to the z-direction.
Extruded Mesh Ratio Real
The ratio of the first and last element when extruding the mesh. The element size is assumed to
follow a geometric division. The same division is applied for all strides resulting from extruding
a node.
Extruded Mesh Density Real
This keyword may be used to define an arbitrary functional dependence that defines the local
mesh parameter h. The function will be automatically called at the interval [0,1]. The same
division is applied for the strides.
Extruded Min Coordinate Real
The mesh is by default extruded to interval [0,1]. This keyword may be used to set a new
minimum coordinate for the extrusion.
Extruded Max Coordinate Real
The mesh is by default extruded to interval [0,1]. This keyword may be used to set a new
maximum coordinate for the extrusion.
Preserve Baseline Integer
Sometimes the user wants to give a boundary condition for the initial baseline of the 2D mesh.
CSC – IT Center for Science
14. Meshing Utilities 65
This keyword may be used to preserve it for setting the boundary condition. Otherwise only the
extruded version of the baseline, i.e. higher dimensional creature will be available.
14.4 Discontinuous interfaces
Discontinuous interfaces are sometimes needed for internal boundary conditions. For example, in the heat
equation a poor thermal contact between objects may lead to a temperature jump at the interface. Such a
jump may be difficult to model with standard elements since they cannot reproduce a discontinuity in the
mesh. For this purpose a second set of nodes at the interface may be created. All the elements that share
nodes at the interface must then be carefully checked to decide on which side of the interface the elements
are located, and set the new element indeces accordingly.
For some nodes on the boundary there is some room for interpretation when creating a discontinuous
boundary. By default all nodes that have even just one hit on the boundary are doubled, i.e. a greedy method
is used. The other alternative is to not include nodes which are attached also to other than boundary nodes
and their immediate parents. This could be applied mainly to structured meshes with some cut through
the mesh. For triangular and tetrahedral meshes the greedy method should be used since in these meshes
elements will always have some orphan nodes on the boundary.
After a discontinuity has been created into the mesh the solution will by default be totally discontinuous,
i.e. no flux will pass through the boundary. To allow again continuous solutions mortar and periodic projec-
tors may be used to ensure continuity. Some field variables may be continuous and some discontinuous over
the interface. The projectors with related jump and offset terms will allow this.
The following keywords are used to control the creation of discontinuous boundaries.
Simulation
Discontinuous Boundary Greedy Logical [True]
Whether to use the greedy method when defining the discontinuous boundary. The user may
choose to set the flag to False freeing some corner nodes that are located in the junction of
continuous and discontinuous interface from being treated as discontinuous.
Discontinuous Bulk Greedy Logical [True]
Similar as the one for boundary element but now the nodes are treated as continuous if they
can be associated with some non-parental bulk element. This probably works correctly only
for structured meshes. In unstructured meshes there are usually many nodes at the boundary
associated with non-parent elements.
Boundary Condition bc_id
Discontinuous Boundary Logical
Create discontinuous boundary on-the-fly. The nodes on the boundary are doubled and heuristics
is used to determine for each element which of the doubled nodes it should use. Discontinuous
boundaries may only be created for internal boundaries. Currently there can be only one discon-
tinuous boundary at a time.
Discontinuous Target Bodies Integer
This keyword may be used to control the owner of the discontinuous boundary (i.e. the "left"
parent element). If not specified the ownership will be given to the parent element with a higher
body index. There may be several discontinuous target bodies. Then the size of the array must
also be specified.
Discontinuous BC Integer
This keyword may be used to activate the creation of the other side of the discontinuous bound-
ary. The other side will be needed if the two boundaries, for example, move with respect to each
other. The value of this keyword relates to the boundary condition index of the new boundary.
Note that also Mortar BC and Periodic BC keywords will have the same effect but only
when they are nonzero! If mortar or periodic boundary condition is requested with a zero ar-
gument, the discontinuous boundary will not be doubled. Instead a cheaper way to create the
mapping between the boundary nodes is used.
CSC – IT Center for Science
Chapter 15
Adaptive Solution
15.1 Introduction
A posteriori error analysis and adaptive mesh refinement are nowadays standard tools in finite element anal-
ysis when cracks, boundary layers, corner singularities, shock waves, and other irregularities are present.
A posteriori error indicators can be used to reveal flaws in finite element discretizations and well-designed
adaptive mesh refinements can reduce the computational costs drastically.
15.2 Theory
Let us consider equilibrium equations of the form
− ∇ · q=fin Ω,(15.1)
q·n=gon Γ,(15.2)
where qis either a flux vector or a second order stress tensor, Ωis a computational domain, Γis a boundary
part, fis an external source or body force, gis an external flux or traction and nis the unit outward normal
to the boundary.
Most symmetric steady state problems described in the model manual of Elmer [] fit in the above
framework of equilibrium equations. To fix ideas, suppose that qis the heat flux satisfying Fourier’s law
q=−k∇T, where Tis the temperature and kis the heat conductivity of the material. We could also think
of qas the stress tensor of linear elasticity. In this case Hooke’s law states that q=E:ε, where Eis
the fourth order tensor of elastic coefficients, ε=symm(∇u)is the linearized strain tensor and uis the
displacement vector.
15.2.1 A posteriori estimate
Let us denote the finite element approximation of qby qhand measure the error q−qhas
ERROR =sZΩ|q−qh|2dΩ(15.3)
Our primary goal is to ensure the accuracy of the solution by imposing the condition
ERROR ≤T OLERANCE (15.4)
where T OLERANCE > 0is an error tolerance prescribed by the user.
In practise, the goal must be replaced by a stronger condition
EST IM AT E ≤T OLERAN CE (15.5)
CSC – IT Center for Science
15. Adaptive Solution 67
where EST IMAT E is a computable functional (of all available data) satisfying
ERROR ≤EST IMAT E (15.6)
Then, if (15.5) holds, (15.4) is satisfied and the quality of the numerical solution is guaranteed.
In Elmer the a posteriori estimate (15.5) is computed from local residuals of the finite element solution
as a weighted sum over the elements,
EST IM AT E =sX
E
η2
E,(15.7)
where ηEis the local error indicator for an individual element E:
η2
E=αEh2
EZE∇ · qh+f2dΩ
+βEX
e in Ω
heZe[[qh·ne]]e2dΓ(15.8)
+γEX
e on Γ
heZeqh·ne−g2dΓ
Here αE,βE, and γE, are local positive constants. The values of these constants depend, among other
things, on the problem to be solved, and must be estimated carefully case by case [].
The first sum in (15.8) is taken over all edges eof Einside the computational domain, the second sum
is taken over all edges on the boundary part Γ,[[·]]eis the jump in (·)across e, and neis a unit normal to the
edge. hEis the size of the element and heis the size of the edge.
The first term on the right-hand-side of (15.8) measures the local residual of the finite element solution
with respect to the equilibrium equation (15.1). The second term measures the discontinuity in the numerical
flux inside Ωand the third term the residual with respect to the boundray condition (15.2).
15.2.2 Adaptivity
The secondary goal of our numerical computations is to find a solution satisfying (15.4) as efficiently as
possible. A nearly optimal solution strategy is obtained by utilizing the property (here we need to impose
some minor restrictions on fand g, see [])
LOCAL ERROR ≥ηE(15.9)
where
LOCAL ERROR =sZE|q−qh|2dΩ(15.10)
The estimate suggests that the error in the numerical solution should be reduced efficiently if the mesh is
refined locally where the indicators ηEare large. Naturally, we can think of coarsening the mesh where the
values of the indicators are small.
The adaptive mesh refinement strategy of Elmer is based on the local estimate (15.9) and on the following
additional assumptions and heuristic optimality conditions:
•The local error behaves as
ηE=CEhpE
E(15.11)
for some constants CEand pE.
•In the optimal mesh the error is uniformly distributed over the elements:
ηE=T OLERAN CE/Nelements (15.12)
CSC – IT Center for Science
15. Adaptive Solution 68
The constants CEand pEin (15.11) can be solved locally for each element if the local errors and the
local mesh sizes are known from at least two different solutions. The second rule (15.12) can then be applied
to extrapolate a new nearly optimal mesh density for the subsequent calculations.
The mesh refinements can be performed eiher by splitting the existing elements into smaller using the
so called RGB-refinement strategy described in [], or by permorming a complete remeshing of the computa-
tional domain using the built-in unstructured mesh generators that produce high quality Delaunay triangula-
tions. In the latter alternative not only mesh refinement is possible, but also local adaptive coarsening.
15.3 Keywords related to the adaptive solution
The adaptive solver of Elmer is activated and controlled by the following keywords in the Solver block of
the solver-input-file.
Adaptive Mesh Refinement Logical
If set to true, then after the solution of the linear system the program computes residual error indicators
for all active elements, estimates the global error, computes a new mesh density and refines the mesh
accordingly.
Adaptive Remesh Logical
If set to true, then a complete remeshing is performed after error estimation using the Mesh2D or
Mesh3D generators. The new mesh density is written in file “bgmesh”. If set to false, then the RGB-
splitting strategy for triangles is applied to perform the refinements.
Adaptive Save Mesh Logical
If set to true, the subsequent meshes are stored in directories RefinedMeshN, where Nis the number
of the adaptive iterate.
Adaptive Error Limit Real
Error tolerance for the adaptive solution.
Adaptive Min H Real
Imposes a restriction on the mesh size. Defualt is zero.
Adaptive Max H Real
Imposes a restriction on the mesh size. Default is infinite.
Adaptive Max Change Real
Controls the change in local mesh density between two subsequent adaptive iterates. Using this key-
word the user can restrict the refinement/coarsening to stabilize the adaptive solution process.
15.4 Implementing own error estimators
Suppose that we are given a subroutine called MySolver for solving the Poisson equation, and we would
like to enhance the code by implementing an a posteriori error indicator for adaptive mesh refinement. The
first thing to do is to take the module Adaptive in use, an define the local error indicators as functions in
an interface block. The beginning of the subroutine should look like the following:
SUBROUTINE MySolver( Model,Solver,Timestep,TransientSimulation )
USE DefUtils
USE Adaptive
INTERFACE
FUNCTION InsideResidual( Model, Element, Mesh, &
Quant, Perm, Fnorm ) RESULT( Indicator )
USE Types
CSC – IT Center for Science
15. Adaptive Solution 69
TYPE(Element_t), POINTER :: Element
TYPE(Model_t) :: Model
TYPE(Mesh_t), POINTER :: Mesh
REAL(KIND=dp) :: Quant(:), Indicator, Fnorm
INTEGER :: Perm(:)
END FUNCTION InsideResidual
FUNCTION EdgeResidual( Model, Edge, Mesh, &
Quant, Perm ) RESULT( Indicator )
USE Types
TYPE(Element_t), POINTER :: Edge
TYPE(Model_t) :: Model
TYPE(Mesh_t), POINTER :: Mesh
REAL(KIND=dp) :: Quant(:), Indicator
INTEGER :: Perm(:)
END FUNCTION EdgeResidual
FUNCTION BoundaryResidual( Model, Edge, Mesh, &
Quant, Perm, Gnorm ) RESULT( Indicator )
USE Types
TYPE(Element_t), POINTER :: Edge
TYPE(Model_t) :: Model
TYPE(Mesh_t), POINTER :: Mesh
REAL(KIND=dp) :: Quant(:), Indicator, Gnorm
INTEGER :: Perm(:)
END FUNCTION BoundaryResidual
END INTERFACE
After these fixed declarations we may proceed normally by defining the local variables, allocate memory
for local tables, integrate the stiffness matrix, set boundary conditions, and solve the problem. Error esti-
mation and adaptive mesh refinements are then performed by calling the subroutine RefineMesh, which
should appear in the code right after the function DefaultSolve.
Norm = DefaultSolve()
IF ( ListGetLogical( Solver % Values, ’Adaptive Mesh Refinement’ ) ) &
CALL RefineMesh( Model, Solver, Potential, Permutation, &
InsideResidual, EdgeResidual, BoundaryResidual )
The functions InsideResidual,EdgeResidual and BoundaryResidual defined in the inter-
face block should finally be contained in MySolve, and return the values of the error indicators described
in the previous section.
As an example, suppose that we are using linear triangles or tetrahedra for solving the Poisson equation.
In this case it holds ∇ · qh= 0 on each element E, and the contribution of the firtst term in (7.1) is simply
InsideResidual =hEsZE|f|2dΩ(15.13)
The function that computes the value of the inside redisual could be written as follows.
FUNCTION InsideResidual( Model, Element, Mesh, &
Quant, Perm, Fnorm ) RESULT( Indicator )
IMPLICIT NONE
TYPE(Model_t) :: Model
CSC – IT Center for Science
15. Adaptive Solution 70
INTEGER :: Perm(:)
REAL(KIND=dp) :: Quant(:), Indicator, Fnorm
TYPE( Mesh_t ), POINTER :: Mesh
TYPE( Element_t ), POINTER :: Element
TYPE(GaussIntegrationPoints_t), TARGET :: IP
TYPE(ValueList_t), POINTER :: BodyForce
REAL(KIND=dp) :: f, hK, detJ, Basis(MAX_NODES), &
dBasisdx(MAX_NODES,3), ddBasisddx(MAX_NODES,3,3), &
Source(MAX_NODES)
LOGICAL :: stat
INTEGER :: n
Indicator = 0.0d0
Fnorm = 0.0d0
hK = element % hK
BodyForce => GetBodyForce( Element )
Source = GetReal( Element, ’Source’ )
IP = GaussPoints( Element )
DO n = 1, IP % n
stat = ElementInfo( Element, Nodes, IP % u(n), IP % v(n), &
IP % w(n), detJ, Basis, dBasisdx, ddBasisddx, .FALSE. )
f = SUM( Source *Basis )
Fnorm = Fnorm + f**2*detJ % IP % s(n)
Indicator = Indicator + f**2*detJ *IP % s(n)
END DO
Fnorm = SQRT( Fnorm )
Indicator = hK *SQRT( Indicator )
END FUNCTION Inside Residual
For the boundary and edge residuals refer to the example Poisson.f90 in the tutorial manual of Elmer.
CSC – IT Center for Science
Chapter 16
Parallel runs
16.1 Introduction
In times of even simple desktop PCs containing multiple CPUs or at least multiple cores, parallel computing
is a necessity to exploit the complete potential of those architectures. As on multi-core architectures multi-
threading (e.g., OpenMP) would be a feasible concept, Elmer utilizes the well established Message Passing
Interface standard for inter-process communication. This approach makes it possible to run Elmer on both,
multi-core as well as multi processor environments.
16.1.1 Parallel Implementation in Elmer
The general concept of a parallel run within Elmer is displayed in Fig. 16.1. Elmer uses domain decompo-
unpartitioned
mesh
combined result
partitioned mesh
parallel
solution
domain
decomposition
unification of
result
Figure 16.1: The principle steps to be taken for a parallel run of Elmer
CSC – IT Center for Science
16. Parallel runs 72
sition for distributing the load to multiple processes that are being run on either different cores or CPUs. To
that end, the initial mesh has to be split into parts that – with respect to the applied models – lead to similar
loads of the processors1. This will be discussed in section 16.2.
The solver stage mostly will demand from serial runs differing numerical techniques, as solution strate-
gies have to take care of the by the domain boundaries limited possibilities of memory access. In general,
convergence of linear systems are more difficult to achieve compared to serial runs. These issues will be
addressed in section 16.3.1.
Finally, as also the output of the parallel runs is split into domains, the post-processing usually demands
an additional step of unifying the split results. Alternatively, new visualization software is capable to do that
either on the fly or to also visualize the results using multiple processes. Especially the latter method in-
evitably will gain importance with the increasing size of models that cannot be handled on a single CPU/core
due to memory and computational constraints. Concepts of post-processing parallel results are discussed in
section 16.4.
16.2 Preprocessing of Parallel Runs
In order to utilize the decomposition, the mesh has to be split into the same amount of partitions, N, as there
are different processes within the parallel computation. The plain and easy way is to start from a mesh for a
serial run. The typical structure of a single domain mesh database is the following:
meshdirectoryname|
|-mesh.header
|-mesh.nodes
|-mesh.elements
|-mesh.boundary
The mesh consists of a header file, containing the basic information (e.g., numbers of nodes and elements),
a file containing all nodes and two further files defining the bulk- and boundary-elements.
The parallel mesh consisting of 2 partitions the is written under the same directory within the sub-
directory partitioning.2:
meshdirectoryname|
|-mesh.header
|-mesh.nodes
|-mesh.elments
|-mesh.boundary
|-partitioning.2|
|-part.1.header
|-part.1.nodes
|-part.1.elements
|-part.1.boundary
|-part.1.shared
|-part.2.header
|-part.2.nodes
|-part.2.elements
|-part.2.boundary
|-part.2.shared
These files basically reflect the structure of a single domain mesh on the partition level. Additionally, a
file names part.N.shared (with Nbeing the partition number) is introduced. It contains – as the name
suggests – information on between domains shared nodes.
1Currently Elmer is not able to perform internal load balancing.
CSC – IT Center for Science
16. Parallel runs 73
16.2.1 Partitioning with ElmerGrid
Provided, a single domain mesh exists, the corresponding ElmerGrid command to create a with respect to
the x-direction split mesh (in our case 2×1×1 = 2 partitions) would read as
ElmerGrid 2 2 meshdirectoryname -partition 2 1 1 0
There are different methods of partitioning built into ElmerGrid. they are summarized in table 16.1
option purpose parameters
-partition NxNyNzFpartition with respect to di-
rections
Nx/y/z . . . number of par-
titions in x/y/z-direction,
F= 0 . . . element-wise par-
titioning, 1 . . . node-wise
partitioning
-partorder nxnynz(optional in additional to
previous) direction of or-
dering
nx/y/z . . . components of
normal vector of ordering
-metis N M uses metis library for parti-
tioning
N. . . number of partitions,
M. . . method
M=0 . . . PartMeshNodal
M=1 . . . PartMeshDual
M=2 . . . PartGraphRecursive
M=3 . . . PartGraphKway
M=4 . . . PartGraphVKway
Table 16.1: Partition methods for ElmerGrid
Depending on what partitioning method was used, additional parameters may be used for adaption of
mesh specifications Those parameters and their purpose are listed in 16.2
CSC – IT Center for Science
16. Parallel runs 74
option purpose parameters
-halo create halo for the parti-
tioning
-indirect create indirect connections
in the partitioning
-periodic FxFyFzdeclare the periodic coor-
dinate directions for paral-
lel meshes and sets peri-
odic points into same par-
titions
Fx.y,z = 1 . . . periodic, 0
. . . not periodic
-partoptim apply aggressive optimiza-
tion to node sharing
-partbw minimize the bandwidth
of partition-partition cou-
plings
-parthypre hypre type numbering
(number the nodes contin-
uously partition-wise)
Table 16.2: Additional mesh generation options for ElmerGrid
Figure 16.2 shows the different distribution of partitions obtained with two different methods. In gen-
Figure 16.2: Distribution of four partitions using the options -partition 2 2 1 (left) and -metis 4
1(right). It comes clear that the partitioning to the left contains more partition-partition boundaries and
consequently will perform worse in a parallel run
eral, the user should utilize the metis options, if more complex geometries (like in Fig. 16.2) are to be
decomposed. This ensures that the number of shared nodes and consequently also the amount of inter-
process communication during the parallel computation is minimized. More simple objects, especially those
using regular meshes, can be split according to the axis using the partition option without having to
compromise on communication speed.
CSC – IT Center for Science
16. Parallel runs 75
Halo Elements
One of the additional options displayed in Tab. 16.2 are so called halo elements. As displayed in Fig. 16.3,
halo-elements are additional elements that do not belong to the partition (i.e., they are not contributing in
the domain-wise solution procedure), but rather are replicas of the neighbor elements of adjoining partitions.
Thus, in a parallel run, the values of variables as well as the geometry of the whole element are at disposition
withing the domain partition. These may be needed by a specific FE method, such as the Discontinuous
Figure 16.3: The concept of halo-elements. Each partition contains information on the neighbor elements
along the domain boundary (red) of the adjoining partitions thus leading to a redundant stripe of elements
(light-gray) that is shared between the domains
Galerkin method or by specific solvers/functions that need additional geometric information from the other
domain (e.g., element-averaged surface normals).
16.3 Parallel Computations in Elmer
As mentioned before, Elmer utilizes Message Passing Interface (MPI) for inter-process communication
while doing a parallel computation. To that end, a special parallel executable that is linked to a MPI li-
brary (the minimum requirement). The compilation process for the MPI version is shortly explained in
chapter 17 of this guide. The executable file of the parallel version of Elmer has a to the serial call different
name, ElmerSolver_mpi. Typically it is executed as an argument to an additional call that is specific
to the parallel (MPI) environment of the platform. For instance, in a typical MPI installation (OpenMPI,
MPICH) the command
mpirun -np 4 ElmerSolver_mpi
will run a four-process parallel Elmer run. The typical screen-output upon launching ElmerSolver_mpi
indicating the number of processes is
ELMER SOLVER (v 5.5.0) STARTED AT: 2009/09/09 10:28:28
ELMER SOLVER (v 5.5.0) STARTED AT: 2009/09/09 10:28:28
ELMER SOLVER (v 5.5.0) STARTED AT: 2009/09/09 10:28:28
ELMER SOLVER (v 5.5.0) STARTED AT: 2009/09/09 10:28:28
ParCommInit: Initialize #PEs: 4
CSC – IT Center for Science
16. Parallel runs 76
MAIN:
MAIN: ==========================================
MAIN: E L M E R S O L V E R S T A R T I N G
MAIN: Library version: 5.5.0 (Rev: 4195)
MAIN: Running in parallel using 4 tasks.
MAIN: HYPRE library linked in.
MAIN: MUMPS library linked in.
MAIN: ==========================================
MAIN:
MAIN:
It is – unlike in the serial version of Elmer – not possible to explicitly add the Solver Input File (SIF, suffix
*.sif) as an argument to the calling command. The user rather has to provide a file called ELMERSOLVER_STARTINFO
containing the file name. Additionally, the in the SIF declared mesh-directory has to contain a mesh with –
in this specific case – four partitions.
16.3.1 Numerical Strategies in Parallel
The concept of domain decomposition means that ElmerSolver is run on N > 1separate parts of a domain
that are interlinked at the boundaries. If no special solvers (see later in this section) are utilized, this in-
herently means that iterative methods have to be used in order to achieve convergence for the linear(ized)
system solution procedure. The selection of available iterative methods, which all fall within Krylov sub-
space methods, is to be found in section 4.3. These methods in general have similar convergence compared
to a single process run. The biggest difference introduced by domain decomposition is, that preconditioning
strategies are altered. To give an example: As only applied to the local matrix, the LU factorization of a
parallel run in comparison to a serial drops the gray zones indicated in Fig. 16.4. This not necessarily will,
Domain 1
Domain 2
Domain 3
Figure 16.4: Difference of ILU factorization between serial and domain decomposition runs. If the factor-
ization is applied only locally within the domain, contributions from the light-gray zones are not accounted
for in the latter
but can negatively affect the convergence of the iterative method.
CSC – IT Center for Science
16. Parallel runs 77
Hypre
Hypre is a library for solving large, sparse linear systems of equations on massively parallel computers.
Hypre was developed at the Center for Applied Scientific Computing (CASC) at Lawrence Livermore Na-
tional Laboratory. Hypre is distributed under the GNU Lesser General Public License and thus not included
in the Elmer distribution. It rather has to be downloaded, compiled and linked together with the Elmer
sources.
The principle keyword for utilizing Hypre is given in the solver section
Linear System Use Hypre Logical
if set to True, Hypre will be used to solve the linear system.
In Elmer, the only Krylov sub-space method being implemented is the Conjugate Gradient Stabilized (BiCGStab)
method, which is taken into use by
Linear System Solver = "Iterative"
Linear System Iterative Method = "BiCGStab"
In combination with the BiCGStab method, the following preconditioner can be taken into use If ILUn
method for preconditioning, the following settings have to be set in the solver section (here with ILU fill-in
level 1):
Linear System Preconditioning String "ILUN"
with Nbeing the fill-in level (just in the built-in Elmer preconditioner). The only significant difference
to Elmer’s built-in ILU preconditioner is, that in case of Hypre, the missing parts (illustrated in Fig.
16.4) are now being passed from one domain to the other. In other words, the precodnitioner should
behave exactly as if it would be applied in a serial, single domain run. This can improve convergence,
but comes at the cost of increased inter-processor communication.
Linear System Preconditioning String "ParaSails"
Parasails is a sparse approximate inverse preconditioner it is preconditioner for sparse matrix systems.
It has the following additional parameters
ParaSails Threshold is a Real value that determines the typical value for matrix entries
being dropped. Its suggested range for direct input (positive sign) is from 0.0 to 0.1, with lower
values demanding higher accuracy and consequently computing time/memory. Alternatively, if
negative values are entered, they are interpreted as the fraction of nonzero elements that are being
dropped (e.g., -0.9 leads to 90/
ParaSails Filter is a Real value that determines the typical value for matrix entries
in the in the computed approximate inverse that are dropped. Its suggested range for direct
input (positive sign) is from 0.0 to 0.05. Alternatively, if negative values are entered, they are
interpreted as the fraction of nonzero elements that are being dropped (see earlier item).
ParaSails Maxlevel is an Integer value that determines the accuracy of the precondi-
tioner. Usually a value of 0 or 1 is within the applicable frame
ParaSails Symmetry is an Integer value that determines the nature of the original ma-
trix. the following settings are to be found from the Hypre manual:
0 non-symmetric and/or indefinite problem, non-symmetric preconditioner
1 Semi Positive Definite (SPD) problem, SPD (factored) preconditioner
2 non-symmetric definite problem, SPD (factored) preconditioner
A typical section for the Navier-Stokes solver being solved with BiCGStab and ParaSails could look
as follows
Solver 1
Equation = "Navier-Stokes"
Optimize Bandwidth = Logical True
Linear System Solver = "Iterative"
CSC – IT Center for Science
16. Parallel runs 78
Linear System Iterative Method = "BiCGStab"
Linear System Max Iterations = 500
Linear System Convergence Tolerance = 1.0E-06
Linear System Abort Not Converged = True
Linear System Preconditioning = "ILU1"
Linear System Residual Output = 1
Linear System Use Hypre = Logical True
Linear System Preconditioning = "ParaSails"
ParaSails Threshold = Real -0.95
ParaSails Filter = Real -0.95
ParaSails Maxlevel = Integer 1
ParaSails Symmetry = Integer 0
Stabilization Method = Stabilized
Nonlinear System Convergence Tolerance = 1.0E-04
Nonlinear System Max Iterations = 30
Nonlinear System Newton After Iterations = 1
Nonlinear System Newton After Tolerance = 1.0E-04
Steady State Convergence Tolerance = 1.0E-03
End
Linear System Preconditioning String "BoomerAMG"
The multi-level procedure given by BoomerAMG can also be utilized as a preconditioner. See the
following part for the particular parameters that can be passed to BoomerAMG
Additionally, the Algebraic MultiGrid (AMG) solver, BoomerAMG, can be used directly to solve the system
by
Linear System Solver = "Iterative"
Linear System Iterative Method = "BoomerAMG"
The following parameters for BoomerAMG can be set
BoomerAMG Relax Type Integer
Defines the smoother on the fine grid and the up and the down cycle. Possible choices
0Jacobi
1Gauss-Seidel, sequential
2Gauss-Seidel, interior points in parallel, boundary sequential
3hybrid Gauss-Seidel or SOR, forward solve (default value)
4hybrid Gauss-Seidel or SOR, backward solve
5hybrid chaotic Gauss-Seidel (does not work with MPI!)
6hybrid symmetric Gauss-Seidel or SSOR
9Gaussian elimination (only on coarsest level)
The solver on the coarsest level is set to Gaussian elimination
BoomerAMG Coarsen Type Integer
Sets the parallel coarsening algorithm to be used. Possible options are
CSC – IT Center for Science
16. Parallel runs 79
0CLJP-coarsening (a parallel coarsening algorithm using independent
sets (default)
1classical Ruge-Stüben coarsening on each processor, no boundary
treatme
3classical Ruge-Stüben coarsening on each processor, followed by a third
pass, which adds coarse points on the boundaries
6Falgout coarsening (uses 1 first, followed by CLJP using the interior
coarse points generated by 1 as its first independent set)
7CLJP-coarsening (using a fixed random vector, for debugging purposes
only)
8PMIS-coarsening (a parallel coarsening algorithm using independent
sets, generating lower complexities than CLJP, might also lead to slower
convergence)
9PMIS-coarsening (using a fixed random vector, for debugging purposes
only)
10 HMIS-coarsening (uses one pass Ruge-Stüben on each processor inde-
pendently, followed by PMIS using the interior C-points generated as
its first independent set)
11 one-pass Ruge-Stüben coarsening on each processor, no boundary treat-
ment
BoomerAMG Num Sweeps Integer
sets the number of sweeps on the finest level (default value = 1)
Boomeramg Max Levels Integer
sets maximum number of MG levels (default value = 25)
BoomerAMG Interpolation Type Integer
Sets parallel interpolation operator. Possible options are
0classical modified interpolation (default)
1LS interpolation (for use with GSMG)
2classical modified interpolation for hyperbolic PDEs
3direct interpolation (with separation of weights)
4multipass interpolation
5multipass interpolation (with separation of weights)
6extended classical modified interpolation
7extended (if no common C neighbor) classical modified interpolation
8standard interpolation
9standard interpolation (with separation of weights)
10 classical block interpolation (for use with nodal systems version only)
11 classical block interpolation (for use with nodal systems version only)
with diagonalized diagonal blocks
12 FF interpolation
13 FF1 interpolation
BoomerAMG Smooth Type Integer
For the use of more complex smoothers. possible options are
6Schwarz smoothers (default and recommended)
7Pilut
8ParaSails
9Euclid
BoomerAMG Cycle Type Integer
For a V-cycle (default) give the value 1for a W-cycle 2
CSC – IT Center for Science
16. Parallel runs 80
BoomerAMG Num Functions Integer
has to be equal to the value given in Variable DOFs
Usually, the default values deliver a good performance and should hence be used as a reference constellation.
Mind also, that BoomerAMG would have more (partly obsolete) option that have not directly been made
accessible through its Elmer interface.
MUMPS
The only implementation of a direct solver in the parallel version of Elmer is MUMPS (http://mumps.enseeiht.fr/,
a sparse matrix multi-frontal parallel direct solver. MUMPS is not licensed under the GNU license terms
and hence is not able to be distributed together with the Elmer sources or binaries, but has to be downloaded
by the user. How MUMPS is integrated into Elmer is explained in chapter 17.
There are certain limits what comes to memory consumption, as for the time being the analysis phase is
done on a master process. As a rule of thumb, about 100k Elements distributed on a quad core machine with
2 Gb of memory per process still works, but does not work fast.
MUMPS is invoked by choosing Linear System Solver = Direct in combination with Linear
System Direct Method = Mumps. A typical call of the Navier-Stokes solver reads like:
Solver 1
Equation = "Navier-Stokes"
Optimize Bandwidth = Logical True
Linear System Solver = Direct
Linear System Direct Method = Mumps
Linear System Convergence Tolerance = 1.0E-06
Steady State Convergence Tolerance = 1.0E-03
Stabilization Method = Stabilized
Nonlinear System Convergence Tolerance = 1.0E-03
Nonlinear System Max Iterations = 1
Nonlinear System Newton After Iterations = 1
Nonlinear System Newton After Tolerance = 1.0E-03
End
Mind, that MUMPS will not work in serial runs. There the UMFPack library should be applied in order
to utilize the same method.
16.4 Post-processing of Parallel Runs
During a parallel run of Elmer the results are also written in parallel. That means, a run with N > 1
partitions/processes produces N > 1output files. If the base name of the output file is parallelrun,
given in the SIF as
Post File = "parallelrun.ep"
the results of a N= 4 parallel run will be written into the mesh-directory as
parallelrun.ep.0
parallelrun.ep.1
parallelrun.ep.2
parallelrun.ep.3
These files contain the results of each domain (starting with zero) for ElmerPost. Similar, if the result file
base name is given as
Post File = "parallelrun.result"
the results of a N= 4 parallel run will be written into the mesh-directory as
CSC – IT Center for Science
16. Parallel runs 81
parallelrun.result.0
parallelrun.result.1
parallelrun.result.2
parallelrun.result.3
From these files, a new Elmer run (on the same partitioned mesh) can be restarted.
16.4.1 Combination of Results with ElmerGrid
The traditional way of post-processing Elmer results is to view them using ElmerPost. It is possible to
directly load one of the domain ElmerPost output files into ElmerPost, viewing only the part of the mesh
corresponding to the part of the domain decompoistion. Usually, the user wants to combine those splitted
results into a single. With the results as displayed above, the correct ElmerGrid command to achieve that
results is
ElmerGrid 15 3 parallelrun
This will go through all timesteps (transient) or steady state levels (steady state) of all partitions and add
them to the file parallelrun.ep (default output name can be changed using the -out option). If lesser
timesteps/steady-state levels have been needed for convergence of the run, ElmerGrid still would try to add
the maximium given number, thus filling in zeros into the combined file. This can be avoided by giving the
additional option -saveinterval first last intervall. For instance, if the user knows that
the run has converged at steady-state level number 8, and is just interested to combine the converged result
the command
ElmerGrid 15 3 parallelrun -saveinterval 8 8 1
executed in the mesh-directory would deliver the result. ElmerGrid asumes that all partitions in the directory
have to be unified into the single output file. If for instance a run of 4 partitions has been run with the
same output name after a run with 6 partitions, the files parallelrun.ep.4 and parallelrun.ep.5
would still reside in the mesh directory, but should not be added. This can be achieved by the additional
option -partjoin 4. This is also usefull if only some parts of the full result should be included.
16.4.2 User Defined Functions and Subroutines in Parallel
In principle, the user does not have to distinguish between a user defined function/solver in parallel and serial
runs. Communication between processes is being taken care of within the built-in Elmer routines, such that
the single solver does not really "see" anything from the parallel environment it is running on. In special
cases a few lines of code might be necessary to deal with differences imposed by parallel runs. Those are
explained in section 18.4.3.
CSC – IT Center for Science
Chapter 17
Compilation and Linking
The main Elmer distribution channel is currently GitHub. The source code may hence be found at
https://github.com/ElmerCSC/elmerfem
and a local copy of the repository can be made by giving the command
git clone git://www.github.com/ElmerCSC/elmerfem
or
git clone https://www.github.com/ElmerCSC/elmerfem
with the files placed under the working directory.
There are different branches in the repository. The most important branch devel contains the latest
updates of Elmer and should preferably be used for compilation. In order to see what is the current branch
use the command
git status
If the branch is not currently devel, one may switch to it by using the command
git checkout devel
The following simple instructions are intended mainly for Unix/Linux users. Additional compilation
instructions may be found via the www-site of Elmer and Wiki pages. For example compiling a parallel
version is not covered in the following.
17.1 Compiling with CMake
To compile, one should first create a separate directory for building. Now, if the working directory is the
place where git clone command was given, one may create the building directory and proceed to it as
mkdir build
cd build
CMake (version 2.8.9 or newer is needed) may then be called to generate makefiles. The following command
may already produce what is needed:
cmake -DWITH_ELMERGUI:BOOL=FALSE -DWITH_MPI:BOOL=FALSE
-DCMAKE_INSTALL_PREFIX=../install ../elmerfem
The command line may include additional options. For example the place of BLAS and LAPACK libraries
can be defined by using options -DBLAS_LIBRARIES and -DLAPACK_LIBRARIES.
If no errors are encountered in connection with running the cmake command, one may proceed to
compilation:
CSC – IT Center for Science
17. Compilation and Linking 83
make install
Finally, one may test the fresh version of Elmer as
ctest
or just perform a specific test (cf. the contents of ../elmerfem/fem/tests for existing tests) for ex-
ample as
ctest -R elasticity
In order to be able to run Elmer from anywhere, one may define the following environment variables:
export ELMER_HOME=$HOME/elmer/install
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$ELMER_HOME/lib
export PATH=$PATH:$ELMER_HOME/bin
17.2 Compiling a user-defined subroutine
The elmerf90 command is provided to help compiling your own solvers. It is a wrapper script to the
compiler that was used to produce the Elmer installation and may be used as
elmerf90 -o MySolver MySolver.f90
In the MinGW system in Windows the suffix .dll should preferably be used
elmerf90 -o MySolver.dll MySolver.f90
After successful compilation, the file mysolver.dll should be found in the local directory. In the file-
name declaration of the Procedure-keyword given in the solver input file, the suffix .dll can be omitted,
so one should be able to define
Solver 1
Procedure = "MySolver" "SubroutineName"
...
End
CSC – IT Center for Science
Chapter 18
Basic Programming
18.1 Introduction
The Elmer distribution contains a large set of different solvers and also the possibility to declare dependence
of material properties or boundary conditions on certain variables (e.g., using the MATC language). Never-
theless, there always may occur certain physical problems that are too difficult to be handled via the solver
input file. Such problems can be coped by introducing new user functions or even complete new solvers.
Elmer is very flexible if certain models have to be added by the user providing her/his own code.
This chapter shall provide a brief introduction to basic programming of user functions as well as solvers
that can be added to Elmer. This will be done by mainly using certain examples and explaining the program-
ming steps occurring in these.
The Elmer Solver source is written in the programming language Fortran 90. Since the Elmer Solver
binaries are compiled as shared objects, it is sufficient to just newly compile the code contributed by the user
as an executable of a shared object (.so in UNIX and .dll in Windows) that dynamically is linked to the
rest of the distribution. In order to provide Elmer with the needed information to be able to load an external
function or solver, the following entry in the solver input file (suffix .sif) has to be given:
Procedure "filename" "procedure"
Where the file filename is the above mentioned executable that should contain the Fortran 90 subroutine
or function procedure. The file filename should reside in the same directory where the solver input
file is. Else, the relative or absolute path to that file should be added in front of the entry filename.
18.2 Basic Elmer Functions and Structures
In order to provide a better understanding for the following mainly example-based explanation some of the
most often needed functions and routines provided by Elmer shall be discussed in this section. Most of these
routines and functions are defined in the Fortran 90 module DefUtils. It has to be included in the code
provided by the user giving the keyword
USE DefUtils
It is important to notice that – due to the nature of the Finite Element Method – the basic data structure
in the Elmer Solver is the single element, rather than single points. That simplifies data manipulation in
solver subroutines, but makes things a little bit more difficult if dealing with the coding of pointwise defined
boundary and initial condition as well as body forces and parameter functions. In the Elmer Solver the type
Element_t contains information on elements.
CSC – IT Center for Science
18. Basic Programming 85
18.2.1 How to Read Values from the Input File
In the Elmer Solver the entries of each section of the solver input file – such as material, body force and
initial condition – are accessed via pointer of the defined data type ValueList_t, further referred to as
“list”. A list provides access to all the information that has been passed to the Elmer Solver from the solver
input files, related to the specific section.
The principal connection between the solver input file and the access from a user function is depicted in
Fig. 18.1
Variable = "MyVar"
Fortran 90 code (*.f90):
SUBROUTINE MySolver( Model,Solver,dt,Transient )
LOGICAL :: Transient
IMPLICIT NONE
TYPE(Solver_t) :: Solver
Type(Model_t) :: Model
REAL(KIND=dp) :: dt
!local variables
TYPE(ValueList_t), POINTER :: listsol
CHARACTER(LEN=MAX_NAME_LEN) :: varname
varname = GetString(listsol,’Variable’,GotIt)
!get keyword Variable
listsol => GetSolverParams()
!get list on Solver Section
LOGICAL :: GotIt
IF (.NOT.GotIt) THEN
CALL FATAL(’MySolver’,’Variable not found’)
END IF
END SUBROUTINE MySolver
...
Solver 1
Equation = "Poisson"
Linear System Solver = "Direct"
Steady State Convergence Tolerance = 1e−06
End
solver input file (*.sif):
Header
Mesh DB "." "mymesh"
End
Simulation
Coordinate System = "Cartesian 2D"
Coordinate Main (3) = 1 2 3
End
Body 1
Equation = 1
Body Force =1
End
Body Force 1
Source = Real 1
End
Equation 1
Active Solvers (1) = 1
End
Solver 1
Equation = "MyEquation"
Output Intervals(1) = 1
Simulation Type = Steady State
Steady State Max Iterations = 1
Post File = "myresult.ep"
Output File = "myresult.dat"
END
Steady State Convergence Tolerance = 1E−06
Procedure = "Poisson" "PoissonSolver"
Variable DOFs = 1
End
MyVar=0
Boundary Condition 1
Linear System Solver = "Direct"
Procedure = "File" "MySolver"
Variable DOFs = 1
Variable = "MyVar"
Target Boundaries (4) = 1 2 3 4
Variable = "MyVar"
Figure 18.1: Scheme of the access of structures in the solver input file from a Fortran 90 user subroutine.
The example shows, how a string is read in from the Solver section.
How to Access Different Sections
The following table shows the definition of the functions defined in the module DefUtils to provide the
correct list for parameters and constants concerning the simulation and solvers
function corresponding section
GetSimulation() Simulation
GetConstants() Constants
GetSolverParams() Solver 1,...
For instance, the following source code lines provide access to the entries in the simulation section of the
solver input file
CSC – IT Center for Science
18. Basic Programming 86
! variable declaration for pointer on list
TYPE(ValueList_t), POINTER :: Simulation
...
! assign pointer to list
Simulation => GetSimulation()
...
Lists that provide information connected to a certain element are
function corresponding section
GetMaterial( Element, Found ) Material 1,...
GetBodyForce( Element, Found ) Bodyforce 1,...
GetEquation( Element, Found ) Equation 1,...
GetBC( UElement ) Boundary Condition 1,...
In the first three of these functions shown above the optional variable Found of type LOGICAL is set to
.TRUE. upon successful search in the solver input file. Hence, it can be used for error handling. The argu-
ments Element and UElement are of type Element_t. If writing a solver, the current element is known
and hence can directly be passed to the functions listed above. Else, this argument may also be omitted.
However, Elmer Solver needs to have the information upon the element in order to inquire the number of the
material/bodyforce/equation/boundary condition section from the solver input file. Hence, if this function
argument is dropped, Elmer Solver falls back to the structure CurrentModel % CurrentElement,
which by the active solver has to be assigned to the address of the current element (see section 18.4).
The functions for input of different values from the solver input file need the assigned pointer to the
corresponding to the specific section.
Reading Constants from the Solver Input File
The following value types are defined for the solver input file:
Value in Input File Variable in Elmer Solver
Real Real(KIND=dp)
Integer INTEGER
Logical LOGICAL
String CHARACTER(LEN=MAX_NAME_LEN)
File CHARACTER(LEN=*)
The defined interface of such a function is
FUNCTION FunctionName( List, Name, Found ) Result(x)
TYPE(ValueList_t), POINTER :: List
CHARACTER(LEN=*) :: Name
LOGICAL, OPTIONAL :: Found
The arguments have the following purpose
List List from which the value has to be read. This pointer has to be obtained by one
of the previously introduced functions
Name The keyword in the particular section for the value
Found Optional boolean variable that contains the value .TRUE. upon successful read
in
The type of the returned of value, x, is depending on the function. The following functions are declared in
the DefUtils module:
•A value of type REAL is read in using the function
CSC – IT Center for Science
18. Basic Programming 87
REAL(KIND=dp) :: r
...
r = GetConstReal( List, Name, Found )
•A variable of type INTEGER is read in using the function
INTEGER :: i
...
i = GetInteger( List, Name, Found )
•A string is read into a user function or solver by the following code line
CHARACTER(LEN=MAX_NAME_LEN) :: str
...
str = GetString( List, Name, Found )
It is important to note that these routines are only meant for reading in constant values. Consequently, these
values must not be dependent on other variables.
Reading Mesh-values from the Solver Input File
The previously introduced function GetConstReal is defined for reading in a constant value of type
REAL(KIND=dp). In the case if values have to be obtained for nodes of an element defined on the mesh
(e.g., an initial condition, a boundary condition or a material parameter), the following function has to be
used
FUNCTION GetReal( List, Name, Found, UElement ) RESULT(x)
TYPE(ValueList_t), POINTER :: List
CHARACTER(LEN=*) :: Name
LOGICAL, OPTIONAL :: Found
TYPE(Element_t), OPTIONAL, TARGET :: UElement
REAL(KIND=dp) :: x(CurrentModel % CurrentElement % Type % NumberOfNodes)
The returned value, x, is a one-dimensional array of type REAL(KIND=dp) with entries for every node of
the either given element UElement or alternatively the default structure CurrentModel % CurrentElement.
For instance, reading in the material parameter Viscosity from an already assigned pointer of type
ValueList_t for a given element, Element, is done by the following code lines
REAL(KIND=dp), ALLOCATABLE :: viscosity(:)
INTEGER :: NoNodes
TYPE(ValueList_t), POINTER :: Material
TYPE(Element_t), POINTER :: Element
LOGICAL :: Found
...
allocate viscosity, set pointers Material and Element
...
NoNodes = GetElementNOFNodes(Element)
...
viscosity(1:NoNodes) = GetReal(Material, ’Viscosity’, Found, Element)
The user has to make sure that the array that later contains the nodal values is of sufficient size. This, for
instance, can be guaranteed by allocating it to the maximal occurring number of nodes for an element in the
model
ALLOCATE(viscosity(CurrentModel % MaxElementNodes))
CSC – IT Center for Science
18. Basic Programming 88
Physical Time as Argument of User Function
If a user function needs physical time as an input, it can be passed as an argument. For instance, if a boundary
condition for the normal component of the velocity would have the physical time as the input variable, the
function call in the solver input file then would look as follows (see section 18.3 for more details on user
functions)
Boundary Condition BCNo
Name = "time_dependent_outlet"
Target Boundaries = BoundaryNo
Normal-Tangential Velocity = True
Velocity 2 = 0.0
Velocity 1
Variable Time
Real Procedure "executable" "timeOutletCondition"
End
End
Here the entries BCNo and BoundaryNo have to be replaced by the correct boundary condition and bound-
ary target number. The file executable should contain the compiled user function timeOutletCondition.
18.2.2 How to Communicate with Structures Inside Elmer Solver
Often it is necessary to get information from inside the Elmer Solver, such as mesh coordinates or field
variables associated to another solver procedure. If writing a solver subroutine, all information of that kind
is accessible via the type TYPE(Solver_t) :: Solver. In the case of a user function (boundary
condition, initial condition, material parameter), the default structure CurrentModel % Solver has to
be used.
Inquiring Information on the Element
As mentioned earlier, most of the pre-defined functions and subroutines inside Elmer Solver apply on the
whole element rather than on single nodes. Information on elements can be accessed via the pre-defined
type Element_t. We list the functions/subroutines for the mostly needed purposes:
•Setting the active element (bulk):
TYPE(Element_t), POINTER :: Element
Type(Solver_t), Target :: Solver
INTEGER :: ElementNumber
...
Element => GetActiveElement(ElementNumber)
The argument Solver is optional. If it is not given, CurrentModel % Solver is used. This
function also automatically sets the pointer CurrentModel % CurrentElement to the element
with the given element number ElementNumber. This is important if sub-sequentially called func-
tions rely on this default value to be set.
The total number of active bulk elements for a specific solver is to be inquired using the value Solver
% NumberOfActiveElements.
•Setting the active element (boundary):
TYPE(Element_t), POINTER :: BoundaryElement
INTEGER :: ElementNumber
...
BoundaryElement => GetBoundaryElement(ElementNumber)
CSC – IT Center for Science
18. Basic Programming 89
This routine also sets the structure CurrentModel % CurrentElement to the boundary ele-
ment.
In contrary to the domain (i.e., bulk) it is a priory not known which boundary element is part of a
boundary condition of a specific solver. This information may be obtained using the function
Type(Element_) :: BoundaryElement
LOGICAL :: IsActiveBoundary
...
IsActiveBoundary = BoundaryElement(BoundaryElement,Solver)
where both arguments are optional. If they are omitted, Elmer Solver takes the values CurrentModel
% CurrentElement and CurrentModel % Solver, respectively. The boundary element num-
ber, ElementNumber may vary between 1 and the maximum value
Solver % Mesh % NumberOfBoundaryElements
•Inquire number of nodes in an element:
INTEGER :: N
TYPE(Element_t) :: Element
...
N = GetElementNOFNodes( Element )
The argument Element is optional. The default value is CurrentModel % CurrentElement
•Get nodal coordinates for element:
TYPE(Nodes_t) :: ElementNodes
TYPE(Element_t) :: Element
TYPE(Solver_t) :: Solver
...
CALL GetElementNodes( ElementNodes, Element, Solver )
The arguments Element and Solver are optional. The default values are CurrentModel %
CurrentElement and CurrentModel % Solver, respectively. The argument ElementNodes
is of the pre-defined type Nodes_t. The different components of the coordinates for the i-th node
can be accessed by
REAL(KIND=dp) :: Xcoord, Ycoord, Zcoord
...
Xcoord = ElementNodes % x(i)
Ycoord = ElementNodes % y(i)
Zcoord = ElementNodes % z(i)
They correspond to the axes of the defined coordinate system in the solver input file.
•Get local coordinates of the i-th node for assigned element:
REAL(KIND=dp) :: U, V, W
TYPE(Element_t), POINTER :: Element
INTEGER :: i
...
U = Element % Type % NodeU(i)
V = Element % Type % NodeV(i)
W = Element % Type % NodeW(i)
CSC – IT Center for Science
18. Basic Programming 90
Local coordinates are corresponding to the position inside the prototype element that is used inside
the Elmer Solver. They are important if parameter values have to be obtained by summation over the
basis functions.
•Get normal vector at the i-th node of the assigned boundary element:
REAL(KIND=dp) :: U, V, Normal(3)
TYPE(Element_t), POINTER :: BoundaryElement
LOGICAL :: CheckDirection
...
U = BoundaryElement % Type % NodeU(i)
V = BoundaryElement % Type % NodeV(i)
Normal = NormalVector( BoundaryElement, Nodes, U, V, CheckDirection )
The function needs the boundary element as well as the local coordinates of the point, where the sur-
face (edge) normal shall be evaluated. The optional last argument, CheckDirection, is a boolean
variable. If set to .TRUE., the direction of the normal is set correctly to the rules given in section
18.3.2. The surface normal is returned in model coordinates and is of unity length.
Inquiring Nodal Values of Field Variables
Nodal values for an element of a scalar variables are read by the subroutine
SUBROUTINE GetScalarLocalSolution( x,name,UElement,USolver )
REAL(KIND=dp) :: x(:)
CHARACTER(LEN=*), OPTIONAL :: name
TYPE(Solver_t) , OPTIONAL, TARGET :: USolver
TYPE(Element_t), OPTIONAL, TARGET :: UElement
The returned value is an array containing the nodal values of the variable name. If this variable name is not
provided, it is assumed that the corresponding solver USolver has only one variable with a single degree
of freedom. If the optional parameters USolver and UElement are not provided, then the default values
CurrentModel % Solver and CurrentModel % CurrentElement, respectively, are used.
For instance, the following code lines read in the nodal element values for the variable Temperature
(from the heat solver)
REAL(KIND=dp), ALLOCATABLE :: localTemp(:)
ALLOCATE(localTemp(CurrentModel % MaxElementNodes))
...
CALL GetScalarLocalSolution(localTemp, ’Temperature’)
In the case of a vector field variable, the analog function
GetVectorLocalSolution has to be used. For instance, if the user wants to read in the local velocity
of an deforming mesh (from the MeshSolver), the following syntax has to be applied
REAL(KIND=dp) , ALLOCATABLE :: localMeshVelocity(:,:)
ALLOCATE(localMeshVelocity(3,Solver % Mesh % MaxElementNodes)
...
CALL GetVectorLocalSolution( MeshVelocity, ’Mesh Velocity’)
Inquiring Values of Field Variables for the Whole Mesh
Sometimes, the user also would like to have values for a field variable of the complete domain accessible.
This is done by assigning a pointer to the variable using the function VariableGet
VariablePointer => VariableGet(Solver % Mesh % Variables, ’Variable’ )
The argument Variable has to be replaced by the variable name. The returned pointer is of type Variable_t.
This type contains the following components
CSC – IT Center for Science
18. Basic Programming 91
component purpose
INTEGER, POINTER :: Perm(:) Contains permutations for the variable. Since
Elmer Solver tries to optimize the matrix struc-
ture, the numbering of the nodal values of the
variable usually does not coincide with the num-
bering of the mesh-nodes. The is used to identify
the mesh-node for a variable-entry. Hence, the
field VariablePointer % Perm(i) con-
tains the nodal number of the i-th value of the
field variable Variable.
INTEGER :: DOFs Contains the amount of vector components of the
variable. DOFs is 1 in case of a scalar, 2 or 3 in
case of a two- or three-dimensional vector field.
REAL(KIND=dp), POINTER ::
Values(:)
contains the values of the field variable
For instance, in order to get access to the temperature field (similar as in the example above), the following
code lines may be used
TYPE(Variable_t), POINTER :: TempVar
INTEGER, POINTER :: TempPerm(:)
REAL(KIND=dp), POINTER :: Temperature(:)
INTEGER :: ElmentNo, N
REAL(KIND=dp), ALLOCATABLE :: localTemp(:)
ALLOCATE(localTemp(CurrentModel % MaxElementNodes))
TYPE(Element_t), POINTER :: Element
...
TempVar => VariableGet( Solver % Mesh % Variables, ’Temperature’ )
IF ( ASSOCIATED( TempVar) ) THEN
TempPerm => TempVar % Perm
Temperature => TempVar % Values
!!!! stop if temperature field has not been found !!!!
ELSE IF
CALL Fatal(’MyOwnSolver’, ’No variable Temperature found’)
END IF
...
DO ElementNo = 1,Solver % NumberOfActiveElements
Element => GetActiveElement(ElementNo)
N = GetElementNOFNodes()
NodeIndexes => Element % NodeIndexes
localTemp(1:N) = Temperature(TempPerm(Element % NodeIndexes))
END DO
It is recommended to check whether the pointer to the variable has been assigned correctly. In our little
example the call of the routine Fatal would stop the simulation run if the assessment would lead to a
negative result.
Inquiring the Current Time
In certain situations in transient simulations the physical time might be needed in a user function. In Elmer
Solver the physical time is treated as a variable and hence has to be read in using the type Variable_t
TYPE(Variable_t), POINTER :: TimeVar
Real(KIND=dp) :: Time
...
TimeVar => VariableGet( Solver % Mesh % Variables, ’Time’ )
Time = TimeVar % Values(1)
CSC – IT Center for Science
18. Basic Programming 92
How to Post Messages
Including PRINT or WRITE statements to stdio in numeric-codes can lead excessive output (large number
of iterations and/or large model sizes) and consequently to a reduction in performance. It is recommended
to use stdio-output routines provided by Elmer Solver, for which to a certain extend the amount of output
can be controlled via the solver input file. The three pre-defined subroutines for output of messages are:
•Info is used for displaying messages (e.g., upon convergence) on the screen. The syntax is
CALL Info(’FunctionName’,’The displayed message’, level=levelnumber)
The first string shall indicate which function the displayed message comes from. The second entry is
a string that contains the message itself. The integer value levelnumber indicates the importance
of the message, starting from 1 (most important). The maximum level of messages being displayed
can be determined in the simulation section of the solver input file
max output level = 3
•Warn is used for displaying warnings. Warnings are always displayed and should be used if conditions
in the code occur that might lead to possible errors
CALL Warn(’FunctionName’,’The displayed warning’)
•Fatal is used to terminate the simulation displaying a message. Consequently, it is used in conditions
in the code where a continuation of the computation is impossible
CALL Fatal(’FunctionName’,’The displayed error message’)
Of course the strings do not have to be hard-coded but can be composed during run-time, using the WRITE
command. The string variable Message is provided for that purpose
WRITE(Message,formatstring)list, of, variables
CALL Info(’FunctionName’,Message, level=3)
The format-string has to be set according to the list of variables.
18.3 Writing a User Function
User functions can be used to provide a pointwise (not element wise!) value for a certain property. They are
used for setting values of material properties inside the domain and/or to set values for certain variables on
boundary conditions at the domain boundary.
The defined interface for a user function looks as follows
FUNCTION myfunction( model, n, var ) RESULT(result)
USE DefUtils
IMPLICIT None
TYPE(Model_t) :: model
INTEGER :: n
REAL(KIND=dp) :: var, result
... contents of the function
leading to a value for variable result ...
END FUNCTION myfunction
CSC – IT Center for Science
18. Basic Programming 93
The function returns the value result. If this is supposed to be a nodal value of a material property or a
variable condition, the variable type in Elmer has to be real double precision, hence, REAL(KIND=dp).
The first argument passed to the function, model, is of the declared type Model_t. It is basically the
interface to all information accessible at runtime, such as variable values, mesh coordinates and boundary
conditions.
The second argument, n, is the integer number of the node for which the function - i.e., the value result -
is is evaluated. Through the last argument, var, a value of an input variable on which the result depends is
passed to the function. In the solver input file this variable is indicated using the Variable-keyword. For
instance (see examples later in this chapter), if the user wants to provide the function above with the result
being the density which – in this particular case – depends on the temperature, the syntax looks as follows
Density = Variable Temperature
Procedure "filename" "myfunction"
Mind that there always has to be an input variable to be given using this keyword. In the case that there is
no need for an input, this may occur as a dummy argument in the function.
18.3.1 User Functions for Domain Properties
In the following we will give an outline of the main issues concerning the preparation of a user function for
a domain property. This might be of use if a material parameter (from material section of the solver input
file), an initial condition (from solver input file initial condition section) or a body force (from solver input
file body force section) of somewhat higher complexity has to be defined for the domain.
Some basic aspects concerning the syntax of such functions shall be explained by the following exam-
ples:
Example: Viscosity as a Function of Temperature
This example demonstrates the following topics:
•definition of a material property dependent on one input variable
•how to read in a value from the material section of the solver input file
We want to implement the following relation for the dynamic viscosity, µ, of a fluid
µ(T) = µ293Kexp C·293
T−1 (18.1)
where Tis the temperature. The parameters µ293K(i.e., the reference viscosity at 293 Kelvin) and Chave to
be read into our function from the material section of the solver input file. Thus, the material section looks
as follows:
! Material section (ethanol)
! --------------------------
Material 1
...
Viscosity = Variable Temperature
Procedure "fluidproperties" "getViscosity"
...
Reference Viscosity = Real 1.2E-03
Parameter C = Real 5.72
...
End
The values µ293K= 1.2·10−3and C= 5.72 are the numerical values of the parameter occurring in
(18.1) for pure ethanol. The executable containing the function named getViscosity will be called
fluidproperties. The header – including the variable declaration – of the function then reads as
follows:
CSC – IT Center for Science
18. Basic Programming 94
!-----------------------------------------------------
! material property function for ELMER:
! dynamic fluid viscosity as a function of temperature
!-----------------------------------------------------
FUNCTION getViscosity( model, n, temp ) RESULT(visc)
! modules needed
USE DefUtils
IMPLICIT None
! variables in function header
TYPE(Model_t) :: model
INTEGER :: n
REAL(KIND=dp) :: temp, visc
! variables needed inside function
REAL(KIND=dp) :: refVisc, C
Logical :: GotIt
TYPE(ValueList_t), POINTER :: material
In order to get the pointer to the correct material-list, we use the function GetMaterial
! get pointer on list for material
material => GetMaterial()
IF (.NOT. ASSOCIATED(material)) THEN
CALL Fatal(’getViscosity’, ’No material found’)
END IF
We immediately check, whether the pointer assignment was successful. In the case of the NULL-pointer
being returned, the pre-defined procedure Fatal will stop the simulation and stdio will display the the
message: (getViscosity): No material-id found
The next step is to read in the numerical values for the parameter:
! read in reference viscosity
refvisc = GetConstReal( material, ’Reference Viscosity’, GotIt)
IF(.NOT. GotIt) THEN
CALL Fatal(’getViscosity’, ’Reference Viscosity not found’)
END IF
! read in parameter C
C = GetConstReal( material, ’Parameter C’, GotIt)
IF(.NOT. GotIt) THEN
CALL Fatal(’getViscosity’, ’Parameter C not found’)
END IF
The variable GotIt contains the state of success of the input. In the case of unsuccessful read-in (i.e., the
variable GotIt shows the value .FALSE.), the simulation will be stopped by the routine Fatal.
Finally, after a check upon the absolute temperature being larger than zero, the viscosity is computed ac-
cording to equation (18.1)
! compute viscosity
IF (temp <= 0.0D00) THEN ! check for physical reasonable temperature
CALL Warn(’getViscosity’, ’Negative absolute temperature.’)
CALL Warn(’getViscosity’, ’Using viscosity reference value’)
visc = refVisc(1)
ELSE
CSC – IT Center for Science
18. Basic Programming 95
visc = refVisc *EXP(C *(2.93D02/temp - 1.0D00))
END IF
END FUNCTION getViscosity
In the case of negative absolute temperature, the reference value will be returned. The pre-defined routine
Warn will display a specific warning on stdio.
Example: Body Force as a Function of Space
For the use as body force for the solver presented in 18.4 (i.e. heat source distribution for heat conduction
equation), we want to write a function that provides a scalar in the domain as a function of space. This
example demonstrates the following topics:
•definition of a scalar in the domain as function of space in two dimensions
•how to inquire the dimension of the problem
We want to implement the following two-dimensional spatial distribution for the scalar h:
h(x, y) = −cos(2 π x)·sin(2 π y), x, y ∈[0,1] (18.2)
where xcorresponds to Coordinate 1 and yto Coordinate 2 of the solver input file.
Since the function given in (18.2) is dependent on two input variables, the single argument that is able to
be passed to the function is not sufficient. Hence it will just be used as dummy argument. Consequently, the
user can provide any (existing!) variable as argument in the solver input file. The header reads as follows
!-----------------------------------------------------
! body force function for ELMER:
! scalar load as function of coordinates x and y
! -cos(2*pi*x)*sin(2*pi*y)
!-----------------------------------------------------
FUNCTION getLoad( model, n, dummyArgument ) RESULT(load)
! modules needed
USE DefUtils
IMPLICIT None
! variables in function header
TYPE(Model_t) :: model
INTEGER :: n
REAL(KIND=dp) :: dummyArgument, load
! variables needed inside function
INTEGER :: DIM
REAL(KIND=dp) :: x, y
Logical :: FirstVisited = .TRUE.
! remember these variables
SAVE DIM, FirstVisited
Further we want to check whether the problem is two-dimensional and hence suitable for our function.
This is done only the first time the function is called, since - hopefully - the dimension of the prob-
lem does not change during all the following calls. The function returning the problem dimension is
CoordinateSystemDimension()
! things done only the first time the function is visited
IF (FirstVisited) THEN
CSC – IT Center for Science
18. Basic Programming 96
! check whether we are dealing with a 2d model
DIM = CoordinateSystemDimension()
IF (DIM /= 2) THEN
CALL FATAL(’getLoad’, ’Dimension of model has to be 2d’)
END IF
FirstVisited = .FALSE.
END IF
The next step to inquire the coordinates of the current point the function is evaluated for. The structure
model contains that information
! get the local coordinates of the point
x = model % Nodes % x(n)
y = model % Nodes % y(n)
Finally, the result is computed
! compute load
load = -COS(2.0D00*PI*x) *SIN(2.0D00*PI*y)
END FUNCTION getLoad
Figure 18.2 shows the result of a simulation using the solver defined in section 18.4 together with the function
getLoad as body force. The entry in the body force section of the solver input file reads as follows:
Figure 18.2: Result obtained with the routine getLoad as body force input for the heat conduction solver
presented in 18.4. The z-coordinate is set proportional to the result obtained in the x-yplane.
Body Force 1
Heat Source
Variable Temp !just a dummy argument
CSC – IT Center for Science
18. Basic Programming 97
Real Procedure "load" "getLoad"
End
where the shared object file has been given the name load. All boundaries are set to adiabatic boundary
condition, i.e., ∇T·~n = 0. This is possible if – and only if – the integral over the load vector over the whole
domain balances to zero, like in our case. Since no Dirichtlet condition is set, the result is not unique and
contains an undetermined offset T0.
18.3.2 User Functions for Boundary Conditions
As for the user functions for bulk properties presented in section 18.3.1, the function for a boundary property
is evaluated pointwise. Hence, the identical function interface is used. The major difference between bulk
elements and elements on boundaries are, that the latter can be associated to more than one body. That is
the case on boundaries between bodies. This is important in such cases where the boundary condition is
dependent on properties inside a specific body to which the boundary is attached. Possible configurations of
boundaries are depicted in Fig. 18.3.
bulk bulk
b) body−body boundary (interface)
a) outside boundary
corner point
other body
outside
Figure 18.3: Possible constellations for boundaries and the normal vectors, ~n at element nodes. Mind, that
for a body-body interface (case b) the default orientation of the surface normal may vary from element to
element. At boundary points that have a discontinuous first derivative of the surface (i.e. at corners and
edges), the evaluation of the normal at the same point for two adjacent boundary elements leads to different
surface normals. Parent elements of boundary elements for the specific body are marked as filled.
If the keyword for checking the directions in the function NormalVector (see section 18.2.2) is set to
.TRUE., the following rules apply:
•In the case of an outside boundary the surface normal, ~n, is always pointing outwards of the body.
•By default on a body-body boundary, the orientation is such that the normal always is pointing towards
the body with lower density, ̺, which is evaluated by comparison of the values given for the keyword
Density in the corresponding material sections of the adjacent bodies.
In certain cases, if densities on both sides are either equal or varying functions of other variables, this may
lead to a varying orientation of the surface. This effect and the effect of different directions of surface
normal for the same point at surface corners and edges is depicted in Fig. 18.3. Whereas the latter effect can
only be dealt with by either producing smooth surfaces or averaging surface normals, the first problem can
CSC – IT Center for Science
18. Basic Programming 98
be overcome by providing the keyword Normal Target Body in the solver input file. This keyword
defines the direction of the surface normals for the pre-defined subroutine NormalVector. For instance,
the following sequence fixes the surface normal at boundary condition number 1 to point into body number
2
Boundary Condition 1
Target Boundaries = 1
Normal Target Body = Integer 2
...
End
Example: Idealized Radiation Law for External Heat Transfer
As an illustrative example we want to show how to implement a very idealized radiation boundary condition
for heat transfer. This example explains:
•how to identify the type of boundary according to Fig. 18.3
•how to get material parameters from domain parent elements of the boundary element
•how to identify the local node number in an elment
•how to inquire boundary normals
The net flux into the body at such a boundary shall be approximated by
qn=ǫ qext −ǫσ ·T4−T4
ext,(18.3a)
where Text is the external temperature, σstands for the Stefan-Boltzmann constant and ǫis the emissivity.
The external heat flux shall be defined as
qext =(−I~s ·~n, ~s ·~n < 0,
0,else,(18.3b)
where Iis the intensity and ~s the direction of the incoming radiation vector. The surface normal, ~n, is
assumed to point outwards the body surface.
Since we are planning to use this boundary condition in connection with the solver presented in section
18.4.2, we have to provide the load vector l=qn/(cp̺)occurring in the force vector of (18.7). This means
that we have to inquire the material parameters cpand ̺for the parent element from the material section of
the adjacent body.
The header of the function reads as
!-----------------------------------------------------
! boundary condition function for ELMER:
! simplified radiation BC
!-----------------------------------------------------
FUNCTION simpleRadiation( model, n, temp ) RESULT(load)
! modules needed
USE DefUtils
IMPLICIT None
! variables in function header
TYPE(Model_t) :: model
INTEGER :: n
REAL(KIND=dp) :: temp, load
! variables needed inside function
CSC – IT Center for Science
18. Basic Programming 99
REAL(KIND=dp), ALLOCATABLE :: Density(:), HeatCapacity(:), ExtTemp(:)
REAL(KIND=dp), POINTER :: Work(:,:)
REAL(KIND=dp) :: radvec(3), Normal(3), NormDir, U, V,&
Emissivity, normalabsorbtion, emission, StefanBoltzmann
INTEGER :: DIM, other_body_id, nboundary, nparent,&
BoundaryElementNode, ParentElementNode, istat
Logical :: GotIt, FirstTime=.TRUE., Absorption = .TRUE.
TYPE(ValueList_t), POINTER :: ParentMaterial, BC
TYPE(Element_t), POINTER :: BoundaryElement, ParentElement
TYPE(Nodes_t) :: ElementNodes
SAVE FirstTime, Density, HeatCapacity, ExtTemp
!-----------------------------------------------------------------------
The boundary element and the pointer to the list of the corresponding Boundary Condition-entry in the solver
input file are assigned
! -----------------------
! get element information
! -----------------------
BoundaryElement => CurrentModel % CurrentElement
IF ( .NOT. ASSOCIATED(BoundaryElement) ) THEN
CALL FATAL(’simpleRadiation’,’No boundary element found’)
END IF
BC => GetBC()
IF ( .NOT. ASSOCIATED(BC) ) THEN
CALL FATAL(’simpleRadiation’,’No boundary condition found’)
END IF
Thereafter, a case distinction between the two possible constellations in Fig. 18.3 (i.e.,outside or body-body
boundary). A outside boundary is indicated by the value -1 of BoundaryElement % BoundaryInfo
% outbody. In the case of a boundary-boundary interface the surface normal is supposed to point out-
wards. I.e., the body hosting the parent element is taken the one for which ParentElement % BodyId
does not coincide with BoundaryElement % BoundaryInfo % outbody
other_body_id = BoundaryElement % BoundaryInfo % outbody
! only one body in simulation
! ---------------------------
IF (other_body_id < 1) THEN
ParentElement => BoundaryElement % BoundaryInfo % Right
IF ( .NOT. ASSOCIATED(ParentElement) )&
ParentElement => BoundaryElement % BoundaryInfo % Left
! we are dealing with a body-body boundary
! and assume that the normal is pointing outwards
! -----------------------------------------------
ELSE
ParentElement => BoundaryElement % BoundaryInfo % Right
IF (ParentElement % BodyId == other_body_id)&
ParentElement => BoundaryElement % BoundaryInfo % Left
END IF
! just to be on the save side, check again
! ----------------------------------------
IF ( .NOT. ASSOCIATED(ParentElement) ) THEN
WRITE(Message, *)&
CSC – IT Center for Science
18. Basic Programming 100
’No parent element found for boundary element no. ’, n
CALL FATAL(’simpleRadiation’,Message)
END IF
After that procedure the pointer ParentElement is set on the adjacent element of the boundary element
inside the body for which the radiation boundary condition is evaluated.
We further need the total number of nodes in the boundary element and the parent element given by
nboundary and nparent, respectively. Also the corresponding number of the boundary node number, n,
in the parent element, ParentElementNode, as well as in the boundary element itself, BoundaryElementNode,
is evaluated. At the end of this code sequence, the routine GetElementNodes sets the information on the
nodes of the boundary element
! get the corresponding node in the elements
! ------------------------------------------
nboundary = GetElementNOFNodes(BoundaryElement)
DO BoundaryElementNode=1,nboundary
IF ( n == BoundaryElement % NodeIndexes(BoundaryElementNode) ) EXIT
END DO
nparent = GetElementNOFNodes(ParentElement)
DO ParentElementNode=1,nboundary
IF ( n == ParentElement % NodeIndexes(ParentElementNode) ) EXIT
END DO
! get element nodes
! -----------------
CALL GetElementNodes(ElementNodes, BoundaryElement)
The needed space for reading material parameter fro the parent element as well as boundary condition
parameters for the boundary element is allocated. In the case of the function being re-visited, we first do a
deallocation, since the values of nboundary or nparent may change from element to element (hybrid
mesh)
! -----------
! Allocations
! -----------
IF (.NOT.FirstTime) THEN
DEALLOCATE(Density, HeatCapacity, ExtTemp)
ELSE
FirstTime = .FALSE.
END IF
ALLOCATE(Density( nparent ), HeatCapacity( nparent ),&
ExtTemp(nboundary), stat=istat)
IF (istat /= 0) CALL FATAL(’simpleRadiation’, ’Allocations failed’)
The following code lines read the values for the parameters associated with the boundary element and the
Stefan-Boltzmann constant from the solver input file
! --------------------------------------
! get parameters from constants section
! and boundary condition section
! --------------------------------------
DIM = CoordinateSystemDimension()
StefanBoltzmann = ListGetConstReal( Model % Constants, &
’Stefan Boltzmann’,GotIt)
IF (.NOT. GotIt) THEN ! default value in SI units
StefanBoltzmann = 5.6704D-08
END IF
CSC – IT Center for Science
18. Basic Programming 101
Emissivity = GetConstReal( BC,’Emissivity’,GotIt)
IF ((.NOT. GotIt) .OR. &
((Emissivity < 0.0d0) .OR. (Emissivity > 1.0d0))) THEN
load = 0.0d00
CALL WARN(’simpleRadiation’,’No Emissivity found.’)
RETURN ! no flux without or with unphysical emissivity
END IF
ExtTemp(1:nboundary) = GetReal( BC,’External Temperature’,GotIt)
IF (.NOT.GotIt) THEN
WRITE (Message,*) ’No external temperature defined at point no. ’, n
CALL Warn(’simpleRadiation’, Message)
ExtTemp(1::nboundary)= temp
END IF
Work => ListGetConstRealArray( BC,’Radiation Vector’,GotIt)
IF ( GotIt ) THEN
radvec = 0.0D00
NormDir = SQRT(SUM(Work(1:DIM,1) *Work(1:DIM,1)))
IF (NormDir /= 0.0d00) THEN
radvec(1:DIM) = Work(1:DIM,1)*Work(4,1)/NormDir
ELSE ! no radiation for weird radiation vector
Absorption = .FALSE.
END IF
ELSE ! no absorption without radiation vector
Absorption = .FALSE.
END IF
If absorption of an incoming radiation vector has to be computed, the surface normal has to be inquired
! ---------------------------------
! get surface normals ( if needed )
! ---------------------------------
IF (Absorption) THEN
U = BoundaryElement % Type % NodeU(BoundaryElementNode)
V = BoundaryElement % Type % NodeV(BoundaryElementNode)
Normal = NormalVector(BoundaryElement, ElementNodes, U, V, .TRUE.)
END IF
Thereafter, the needed material parameters are read from the material section of the solver input file that
associated with the body for which the radiation boundary condition is computed
! -------------------------------------------
! get material parameters from parent element
! -------------------------------------------
ParentMaterial => GetMaterial(ParentElement)
! next line is needed, if the consequently read
! parameters are themselves user functions
! ---------------------------------------------
CurrentModel % CurrentElement => ParentElement
Density(1:nparent) = GetReal(ParentMaterial, ’Density’, GotIt)
IF (.NOT.GotIt) Density(1:nparent) = 1.0d00
HeatCapacity(1:nparent) = GetReal(ParentMaterial, ’Heat Capacity’, GotIt)
IF (.NOT.GotIt) HeatCapacity(1:nparent) = 1.0d00
! put default pointer back to where it belongs
CSC – IT Center for Science
18. Basic Programming 102
! --------------------------------------------
CurrentModel % CurrentElement => BoundaryElement
Since these material parameters themselves may be given in form of user functions, the default pointer
CurrentModel % CurrentElement has to be set to the parent element upon call of the function
GetReal.
Finally the two parts of the total normal heat flux are evaluated. The external load is obtained by dividing
this Laue by the inquired values for Density and HeatCapacity.
!------------------------------------
! compute flux and subsequently load
!-----------------------------------
IF (Absorption) THEN
normalabsorbtion = -1.0D00 & ! since normal pointing outwards
*Emissivity *DOT_PRODUCT(Normal,radvec)
IF (normalabsorbtion < 0.0d0) &
normalabsorbtion = 0.0d00
ELSE
normalabsorbtion = 0.0d00
END IF
emission = StefanBoltzmann *Emissivity *&
( temp**4 - ExtTemp(BoundaryElementNode)**4) &
/ (HeatCapacity(ParentElementNode) *Density(ParentElementNode))
load = normalabsorbtion + emission
END FUNCTION simpleRadiation
Figure 18.4 shows the result of the heat conduction solver presented in section 18.4 in combination with
the function simpleRadiation as boundary condition on two oppositely directed boundaries (boundary
condition no. 1). Since the radiation vector is aligned with the x-direction and hence perpendicular with
respect to these two boundaries, the absorption part vanishes for one of these boundaries. For the sake of
simplicity, the material parameters ̺,cpand khave been set to unity. The corresponding entries of the solver
input file for the boundary conditions of the case shown in Fig. 18.4 are:
Boundary Condition 1
Target Boundaries(2) = 1 2
Boundary Type = String "Heat Flux"
External Load
Variable Temp
Real Procedure "radiation_flux.exe" "simpleRadiation"
External Temperature = Real -20.0E00
Radiation Vector(4) = Real -1.0E00 0.0E00 0.0E00 1.0E01
Emissivity = Real 0.5
End
Boundary Condition 2
Target Boundaries = 3
Boundary Type = String "Given Temperature"
Temp = Real 0
End
Boundary Condition 3
Target Boundaries = 4
Boundary Type = String "Adiabatic"
End
CSC – IT Center for Science
18. Basic Programming 103
Figure 18.4: Result of the heat conduction solver applying the simplified radiation boundary condition
described in this section. The coordinate directions as well as the number of the corresponding boundary
condition section are shown. The latter can be compared to the solver input file entries shown in this section.
CSC – IT Center for Science
18. Basic Programming 104
18.4 Writing a Solver
In Elmer an additional solver may be provided by dynamically linked subroutines. The pre-defined interface
of a solver subroutine is
SUBROUTINE mysolver( Model,Solver,dt,TransientSimulation )
TYPE(Model_t) :: Model
TYPE(Solver_t) :: Solver
REAL(KIND=dp) :: dt
LOGICAL :: TransientSimulation
The first argument, Model, is the same structure also passed to a user function (see section 18.3). The
second argument, Solver, is of type Solver_t and contains all information upon options set for this
particular solver. The timestep size, dt, and a boolean variable, TransientSimulation, indicating
whether the solver is to be run as a part of a transient (value .TRUE.) or steady state (value .FALSE.)
simulation are further arguments of the solver subroutine.
18.4.1 Structure of a Solver
The well known structure of a linearized partial differential equation (PDE) for the scalar Tin the variational
formulation is
Mij
∂Tj
∂t +Aij Tj=Fi,(18.4)
with the mass matrix, Mij , the stiffness matrix, Aij and the force vector, Fi.
Time-stepping and the coupled solver iteration – i.e., the steady state or time level iteration of several
solvers of the simulation – is taken care of by the main part of the Elmer Solver. The task of the user
supplying a user defined solver subroutine is to linearize an eventually nonlinear PDE and to supply the
Elmer Solver with the element-wise components for the matrices as well as the force vector.
This leads to a principle structure of a user defined solver subroutine as it is depicted in Fig. 18.5. We
further will explain the construction of a user solver subroutine by discussing an example case.
18.4.2 Example: Heat Conduction Equation
In order to provide a simple example, we want to explain how the solution of the purely convective heat
transport equation is implemented in Elmer. This example explains:
•reading solver parameters from the solver input file
•assembly of the local element matrix components for the domain and a Neumann condition including
explanation of the most important Elmer Solver routines needed for linear system matrix assembly
and solution
•how to deal with non-linear and transient PDE’s
•how to implement Dirichlet boundary conditions
For constant density, ̺, and heat capacity, cpthis equation may be written as
∂T
∂t − ∇ · (k
cp̺∇T) = h
cp̺,(18.5)
where Tstands for the temperature, kfor the heat conductivity and his the heat source.
The variational formulation of (18.5) reads after partial integration of the conduction term and application
of Green’s theorem
Z
V
∂T
∂t γidV +Z
V
k
cp̺∇T· ∇γidV =
Z
V
h
cp̺γidV +I
∂V
1
cp̺(k∇T)·~n
|{z }
=qn
γidA, (18.6)
CSC – IT Center for Science
18. Basic Programming 105
subroutine inside
the solver routine
subroutine inside
the solver routine
often provided as
often provided as
ElmerSolver Main
until last timestep
relative change of norms < Steady State Tolerance
Timestepping loop
User Subroutine
Steady state iteration (coupled system)
Initialization
Nonlinear iteration loop
relative change of norms < Nonlinear Tolerance
Domain element loop
Matrix assembly for domain element
until last bulk element
Boundary element loop
Matrix assembly for von Neumann and
Newton conditions at boundary element
until last boundary element
set Dirichlet boundary conditions
solve the system
or
nonlinear max. iterations exceeded
Figure 18.5: Flowchart of a user defined solver subroutine within Elmer Solver. The grey underlayed area
indicates the tasks of a simulation that are provided by Elmer, whereas the white area contains the flowchart
for the source code of the user subroutine. Arrows pointing into this field indicate needed subroutine/function
calls to Elmer Solver.
CSC – IT Center for Science
18. Basic Programming 106
where γiis the basis-function, Vand ∂V is the element volume and its enclosing surface, respectively. The
surface normal of ∂V is denoted by ~n. According to the Galerkin method, the variable is given as T=Tjγj,
with the sum taken over the index j. Comparing with (18.4) leads to the matrices and vectors
Mij =Z
V
γjγidV,
Aij =Z
V
k
cp̺∇γj· ∇γidV,
Fi=Z
V
h
cp̺∇γidV +I
∂V
qn
cp̺
|{z}
=l
γidA.
(18.7)
Although the external heat flux perpendicular to the surface normal, qn(T), in general is a function of
the temperature we want to keep it formal as being prescribed. Hence, only a contribution in the force
vector occurs in our formulation. Mind, that a linear or linearized expression of qn(T)directly could give a
contribution to the stiffness matrix at the entries corresponding to boundary nodes. In our approach, even in
the case of a linear dependency qn(T)∝Twe have to iterate the solution because of our treatment of the
boundary condition.
The header contains the declaration needed variables – we tried to give them self explaining identi-
fiers. Furthermore, allocations of the needed field arrays are done for the first time the subroutine is visited
(checked by the boolean variable AllocationsDone). The variable names of these arrays then have to
be included in the SAVE-statement at the end of the variable declaration block.
SUBROUTINE MyHeatSolver( Model,Solver,dt,TransientSimulation )
USE DefUtils
IMPLICIT NONE
!----------------------------------------------------------------
TYPE(Solver_t) :: Solver
TYPE(Model_t) :: Model
REAL(KIND=dp) :: dt
LOGICAL :: TransientSimulation
!----------------------------------------------------------------
! Local variables
!----------------------------------------------------------------
TYPE(Element_t),POINTER :: Element
LOGICAL :: AllocationsDone = .FALSE., Found, Converged
INTEGER :: n, t, istat, other_body_id, iter, NonlinearIter
REAL(KIND=dp) :: Norm, PrevNorm=0.0d00, NonlinearTol, RelativeChange
TYPE(ValueList_t), POINTER :: BodyForce, Material, BC, SolverParams
REAL(KIND=dp), ALLOCATABLE :: MASS(:,:), STIFF(:,:), LOAD(:), FORCE(:)
REAL(KIND=dp), ALLOCATABLE :: HeatCapacity(:), HeatConductivity(:),&
Density(:), ExternalTemp(:), TransferCoeff(:), HeatFlux(:)
CHARACTER(LEN=MAX_NAME_LEN) :: BoundaryType
SAVE MASS, STIFF, LOAD, FORCE,&
HeatCapacity, HeatConductivity,&
Density, ExternalTemp, TransferCoeff, HeatFlux, &
CSC – IT Center for Science
18. Basic Programming 107
AllocationsDone, PrevNorm
!----------------------------------------------------------------
!Allocate some permanent storage, this is done first time only:
!--------------------------------------------------------------
IF ( .NOT. AllocationsDone ) THEN
N = Solver % Mesh % MaxElementNodes !big enough for elemental arrays
ALLOCATE( FORCE(N), LOAD(N), MASS(N,N), STIFF(N,N), STAT=istat )
IF ( istat /= 0 ) THEN
CALL Fatal( ’MyHeatSolve’,&
’Memory allocation error for matrix/vectors.’ )
END IF
ALLOCATE( HeatCapacity(N), HeatConductivity(N), Density(N),&
ExternalTemp(N), TransferCoeff(N), HeatFlux(N), STAT=istat )
IF ( istat /= 0 ) THEN
CALL Fatal( ’MyHeatSolve’,&
’Memory allocation error for parameter arrays.’ )
END IF
AllocationsDone = .TRUE.
END IF
In the next step, information on the nonlinear iteration is being read from the solver section of the solver
input file
!Read in solver parameters
!-------------------------
SolverParams => GetSolverParams()
IF (.NOT. ASSOCIATED(SolverParams))&
CALL FATAL(’MyHeatSolve’,’No Solver section found’)
NonlinearIter = GetInteger(SolverParams, &
’Nonlinear System Max Iterations’, Found)
IF ( .NOT.Found ) NonlinearIter = 1
NonlinearTol = GetConstReal( SolverParams, &
’Nonlinear System Convergence Tolerance’, Found )
IF ( .NOT.Found ) NonlinearTol = 1.0D-03
Therafter, the nonlinear iteration loop (outermost loop of the white underlayed area in Fig. 18.5) is started
and the linear system solver is initialized (routine DefaultInitialize)
!----------------------------------------------------------------
! Nonlinear iteration loop
!----------------------------------------------------------------
DO iter=1,NonlinearIter
Converged = .FALSE.
WRITE(Message,’(A,I5,A,I5)’) ’Nonlinear iteration no.’,iter,&
’ of max. ’, NonlinearIter
CALL INFO(’MyHeatSolve’,Message,level=1)
!Initialize the system and do the assembly:
!------------------------------------------
CALL DefaultInitialize()
The next loop is over all elements in the simulation domain our solver has been assigned to (Solver %
NumberOfActiveElements). The function GetActiveElement inquires the element associated
with the element number t. This call at the same time also sets the default pointer CurrentModel %
CurrentElement to that particular element, which is important if subsequently called functions rely on
CSC – IT Center for Science
18. Basic Programming 108
this pointer to be set correctly (see section 18.3). After inquiring the number of nodes the nodal material pa-
rameter values cp→HeatCapacity(1:n),k→HeatConductivity(1:n) and ̺→Density(1:n)
are read in. If one of these parameters is not found (i.e., Found == .FALSE.), a default value of 1 will
be set in order to avoid division by zero.
!----------------------------------------------------------------
! Assembly for the domain
!----------------------------------------------------------------
DO t=1,Solver % NumberOfActiveElements
! get element info
!-----------------
Element => GetActiveElement(t)
n = GetElementNOFNodes()
! get material parameters
! ----------------------
Material => GetMaterial()
IF (.NOT. ASSOCIATED(Material)) THEN
WRITE(Message,’(A,I5,A)’) &
’No material for bulk element no. ’,t,’ found.’
CALL FATAL(’MyHeatSolve’,Message)
END IF
HeatCapacity(1:n) = GetReal(Material, ’Heat Capacity’, Found )
IF (.NOT. Found) HeatCapacity(1:n) = 1.0D00
HeatConductivity(1:n) = &
GetReal(Material, ’Heat Conductivity’, Found )
IF (.NOT. Found) HeatCapacity(1:n) = 1.0D00
Density(1:n) = GetReal(Material, ’Density’, Found )
IF (.NOT. Found) Density(1:n) = 1.0D00
In order to call the subroutine taking care of the composition of the element matrices and force vector
(subroutine LocalMatrix), the load vector – in our case the heat source, h→LOAD(1:n) – has to
be read from the body section of the solver input file. In the case of a transient simulation (indicated by
TransientSimulation == .TRUE.) the first order time discretization is accounted for using the
subroutine Default1stOrderTime. Mind, that also higher order time discretization routines would be
at hand. The local matrix is added to the global coefficient matrix of Elmer Solver calling the subroutine
DefaultUpdateEquations
!Get load for force vector
!-------------------------
LOAD = 0.0d0
BodyForce => GetBodyForce()
IF ( ASSOCIATED(BodyForce) ) &
LOAD(1:n) = GetReal( BodyForce, ’Heat Source’, Found )
!Get element local matrix and rhs vector:
!----------------------------------------
CALL LocalMatrix( MASS, STIFF, FORCE, LOAD, Element, n,&
HeatCapacity, HeatConductivity, Density, TransientSimulation)
!Update global matrix and rhs vector from local matrix & vector:
!---------------------------------------------------------------
IF ( TransientSimulation ) THEN
CALL Default1stOrderTime( MASS,STIFF,FORCE )
END IF
CSC – IT Center for Science
18. Basic Programming 109
CALL DefaultUpdateEquations( STIFF, FORCE )
!----------------------------------------------------------------
END DO ! end Assembly for the domain
!----------------------------------------------------------------
After the bulk elements, the contribution to the coefficient matrix and the force vector from a Neumann type
of boundary condition has to be taken into account. To this end, we are looping over all boundary elements.
Their total number is given by Solver % Mesh % NumberOfBoundaryElements. The routine
ActiveBoundaryElement checks whether the previously inquired element is part of a boundary condi-
tion that has been assigned to our solver. If the value 1 is returned from the function GetElementFamily
– i.e. we are dealing with boundary element given at a point element – the element also will be skipped,
since Neumann conditions cannot be set on such elements. Finally, the list-pointer to the associated boundary
condition section (GetBC) is set and the local matrices and vectors are initiated to zero.
!----------------------------------------------------------------
! assembly of Neumann boundary conditions
!----------------------------------------------------------------
DO t=1, Solver % Mesh % NumberOfBoundaryElements
! get element and BC info
! -----------------------
Element => GetBoundaryElement(t)
IF ( .NOT.ActiveBoundaryElement() ) CYCLE
n = GetElementNOFNodes()
! no evaluation of Neumann BC’s on points
IF ( GetElementFamily() == 1 ) CYCLE
BC => GetBC()
FORCE = 0.0d00
MASS = 0.0d00
STIFF = 0.0d00
Since we have to define between different types of boundary conditions, we inquire the contents of a keyword
Boundary Type from the solver input file. If this string is equal to ’heat flux’, the variable with the
name ’External Load’ will be read in from the boundary condition list BC. Thereafter, the contribution
to the force vector will be computed by the internal subroutine BoundaryCondition (see later in this
code). Mind, that the external load is not the given heat flux, qn, but its value divided by the heat capacity
and the density, l=qn/(cp̺). This has to be taken care of if a numerical value or even a user function
is provided in the solver input file (see section 18.3.2). In the case of no boundary type being found or
an unknown string being detected, the natural boundary condition, zero flux perpendicular to the surface,
will be assumed. This is equivalent to the ’adiabatic’ boundary condition. In the case of ’given
temperature’ the natural boundary condition will be altered by the matrix manipulation arising from
the Dirichlet boundary condition (see later in this code).
! check type of boundary and set BC accordingly
!----------------------------------------------
BoundaryType = GetString(BC,’Boundary Type’,Found)
IF (.NOT. Found) CYCLE
! natural boundary condition
IF ((BoundaryType == ’adiabatic’)&
.OR. (BoundaryType == ’given temperature’)) THEN
CYCLE
! external heat flux
ELSE IF(BoundaryType == ’heat flux’) THEN
! get external load; mind that this is the heat flux
CSC – IT Center for Science
18. Basic Programming 110
! divided by the density and heat capacity
LOAD(1:n) = LOAD(1:n) + GetReal(BC,’External Load’, Found)
! do the assembly of the force vector
CALL BoundaryCondition(LOAD, FORCE, Element, n)
ELSE
WRITE(Message,’(A,I3,A)’)&
’No boundary condition given for BC no. ’,GetBCId(),&
’. Setting it to adiabatic.’
CALL WARN(’MyHeatSolve’,Message)
CYCLE
END IF
The boundary element loop is closed after the component system matrix and vector are updated for the
current boundary element.
IF ( TransientSimulation ) THEN
MASS = 0.d0
CALL Default1stOrderTime( MASS, STIFF, FORCE )
END IF
CALL DefaultUpdateEquations( STIFF, FORCE )
!----------------------------------------------------------------
END DO ! end of assembly of Neumann boundary conditions
!----------------------------------------------------------------
Before setting the Dirichlet conditions (i.e., given boundary temperature T) using the subroutine DefaultDirichletBCs()
it is important to finish the element-wise assembly of the Elmer Solver system matrix calling DefaultFinishAssembly
CALL DefaultFinishAssembly()
! call Elmer Solver routine for Dirichlet BCs
! ------------------------------------------
CALL DefaultDirichletBCs()
The system is solved by the function call DefaultSolve, which returns the norm, Nnof the solution
vector Tjfor the n-th nonlinear iteration step. This is needed in order to inquire the change of the solution
between two steps. If the relative norm
R= 2 |Nn−1−Nn|
Nn−1+Nn,
is smaller than the given tolerance Nonlinear System Tolerance of the solver section, then the
nonlinear iteration is taken to be converged.
! Solve the system
! ----------------
Norm = DefaultSolve()
! compute relative change of norm
! -------------------------------
IF ( PrevNorm + Norm /= 0.0d0 ) THEN
RelativeChange = 2.0d0 *ABS( PrevNorm-Norm ) / (PrevNorm + Norm)
ELSE
RelativeChange = 0.0d0
END IF
WRITE( Message, *) ’Result Norm : ’,Norm
CSC – IT Center for Science
18. Basic Programming 111
CALL Info( ’MyHeatSolve’, Message, Level=4 )
WRITE( Message, *) ’Relative Change : ’,RelativeChange
CALL Info( ’MyHeatSolve’, Message, Level=4 )
! do we have to do another round?
! -------------------------------
IF ( RelativeChange < NonlinearTol ) THEN ! NO
Converged = .TRUE.
EXIT
ELSE ! YES
PrevNorm = Norm
END IF
!----------------------------------------------------------------
END DO ! of the nonlinear iteration
!----------------------------------------------------------------
After leaving the nonlinear iteration loop the status of convergence shall be displayed on stdio
! has non-linear solution converged?
! ----------------------------------
IF ((.NOT.Converged) .AND. (NonlinearIter > 1)) THEN
WRITE( Message, *) ’Nonlinear solution has not converged’,&
’Relative Change=’,RelativeChange,’>’,NonlinearTol
CALL Warn(’MyHeatSolve’, Message)
ELSE
WRITE( Message, *) ’Nonlinear solution has converged after ’,&
iter,’ steps.’
CALL Info(’MyHeatSolve’,Message,Level=1)
END IF
In the code lines given above, the user could exchange the routine Warn by Fatal if the simulation should
stop upon failed nonlinear iteration convergence.
Further we have to include the needed local subroutines using the Fortran 90 command
!----------------------------------------------------------------
! internal subroutines of MyHeatSolver
!----------------------------------------------------------------
CONTAINS
The subroutine LocalMatrix composes the local matrices and vectors for a bulk element. The header
with the variable declaration reads as follows
!----------------------------------------------------------------
SUBROUTINE LocalMatrix(MASS, STIFF, FORCE, LOAD, Element, n, &
HeatCapacity, HeatConductivity, Density, TransientSimulation)
IMPLICIT NONE
!----------------------------------------------------------------
REAL(KIND=dp), DIMENSION(:,:) :: MASS, STIFF
REAL(KIND=dp), DIMENSION(:) :: FORCE, LOAD, &
HeatCapacity, HeatConductivity, Density
INTEGER :: n
TYPE(Element_t), POINTER :: Element
LOGICAL :: TransientSimulation
!----------------------------------------------------------------
REAL(KIND=dp) :: Basis(n),dBasisdx(n,3),ddBasisddx(n,3,3)
REAL(KIND=dp) :: detJ, LoadAtIP,&
CSC – IT Center for Science
18. Basic Programming 112
LocalHeatCapacity, LocalHeatConductivity, LocalDensity
LOGICAL :: Stat
INTEGER :: t,i,j,DIM
TYPE(GaussIntegrationPoints_t) :: IP
TYPE(Nodes_t) :: Nodes
SAVE Nodes
For the sake of simplicity we use the same identifiers as in the solver subroutine for the variables in the
argument list.
The next step is to inquire the dimension of the coordinate system. Thereafter, we get the nodes of the
element using the already introduced function GetElementNodes. Since the values in CurrentModel
% CurrentElement and CurrentModel % Solver have been set, no additional arguments to the
variable Nodes have to be set. After we have initialized the local matrix and vector components to zero, the
information upon the Gauss-points needed for integration is inquired by the function GaussPoints. They
returned variable IP is of type GaussIntegrationPoints_t.
DIM = CoordinateSystemDimension()
CALL GetElementNodes( Nodes )
STIFF = 0.0d0
FORCE = 0.0d0
MASS = 0.0d0
!Numerical integration:
!----------------------
IP = GaussPoints( Element )
The integration over the element is done by summing over all Gauss-points (from 1 to NIP →IP % n.
The square root of the determinant of the element coordinate system metric tensor pdet(JT·J)→DetJ
as well as the local basis functions, γi→Basis, their derivatives, ∇γi→dBasisdx, are evaluated
for every Gauss-point using the function ElementInfo. The variable ddBasisddx is just passed as a
dummy argument, since the last argument, getSecondDerivatives is set to .FALSE.. The pointer to
the element, Element, and its nodes, Nodes and the local variables of the Gauss-points IP % U(t),IP
% V(t) and IP % W(t), are needed as input.
!----------------------------------------------------------------
! Loop over Gauss-points (element Integration)
!----------------------------------------------------------------
DO t=1,IP % n
!Basis function values & derivatives at the integration point:
!-------------------------------------------------------------
getSecondDerivatives = .FALSE.
stat = ElementInfo( Element, Nodes, IP % U(t), IP % V(t), &
IP % W(t), detJ, Basis, dBasisdx, ddBasisddx, &
getSecondDerivatives)
Thereafter, the material parameters at the Gauss-points are evaluated, using the basis function. For instance,
the local density, ̺|IP →LocalDensity at the Gauss-point is evaluated as follows:
̺|IP =̺iγi|IP,
with the sum taken over the nodal index i. The load vector h/(̺ cp)|IP →LoadAtIP is evaluated in a
similar way.
CSC – IT Center for Science
18. Basic Programming 113
!Material parameters at integration point:
!----------------------------------------
LocalHeatCapacity = SUM( HeatCapacity(1:n) *Basis(1:n) )
LocalHeatConductivity = SUM( HeatConductivity(1:n) *Basis(1:n) )
LocalDensity = SUM( Density(1:n) *Basis(1:n) )
!The source term at the integration point:
!-----------------------------------------
LoadAtIP = SUM( Basis(1:n) *LOAD(1:n) ) &
/(LocalHeatCapacity *LocalDensity)
The force vector is obtained by the integral over the element, which is approximated by the sum over all
Gauss-point contributions
Fj=Z
V
h
̺ cp
γjDV ≈
NIP
X
t=1 √ds2qdet(JT·J)h
̺ cp
γj|IP.
The model coordinate system metric √ds2→IP % s(t) as well as the previously inquired element
coordinate system metric pdet(JT·J)→DetJ have to be taken into account.
The matrix components are evaluated in complete analogy to the force vector
Mij =Z
V
γjγidV ≈
NIP
X
t=1 √ds2qdet(JT·J)γjγi|IP,
Aij =Z
V
k
cp̺∇γj· ∇γidV ≈
NIP
X
t=1 √ds2qdet(JT·J)k
cp̺(∇γj· ∇γi)|IP,
where the dot product of the first derivatives of the basis function is implemented using the expression
∇γj· ∇γi→SUM(dBasisdx(i,1:DIM) *dBasisdx(j,1:DIM))
!----------------------------------------------------------------
! Loop over j-components of matrices and over force vector
!----------------------------------------------------------------
DO j=1,n
FORCE(j) = FORCE(j) + IP % s(t) *DetJ *LoadAtIP *Basis(j)
!----------------------------------------------------------------
! Loop over i-components of matrices
!----------------------------------------------------------------
DO i=1,n
!The mass matrix, if needed
!--------------------------
IF (TransientSimulation) THEN
MASS(i,j) = MASS(i,j)+ IP % s(t) *DetJ *&
Basis(i)*Basis(j)
END IF
!Finally, the stiffness matrix:
!------------------------------
STIFF(i,j) = STIFF(i,j) &
+ IP % s(t) *DetJ *LocalHeatConductivity &
*SUM(dBasisdx(i,1:DIM) *dBasisdx(j,1:DIM))&
/(LocalDensity *LocalHeatCapacity)
!----------------------------------------------------------------
END DO ! end Loop over i-components of matrices
!----------------------------------------------------------------
CSC – IT Center for Science
18. Basic Programming 114
END DO ! end Loop over j-components of matrices and vector
!----------------------------------------------------------------
END DO ! end Loop over Gauss-points (element Integration)
!----------------------------------------------------------------
END SUBROUTINE LocalMatrix
!----------------------------------------------------------------
The last two statements in the code sequence given above close the loop over the Gauss-points and provide
the end statement of the local subroutine LocalMatrix.
The subroutine BoundaryCondition evaluates the contribution to the force vector at the boundary ele-
ments with given external load l=qn/(cp̺)→LOAD
Fj=I
∂V
l dV ≈
NIP
X
t=1 √ds2qdet(JT·J)l|IP.
Since this is implemented in complete analogy to the assembly of the force vector in the previously discussed
subroutine LocalMatrix, a detailed explanation can be omitted
!----------------------------------------------------------------
SUBROUTINE BoundaryCondition(LOAD, FORCE, Element, n)
IMPLICIT NONE
!----------------------------------------------------------------
REAL(KIND=dp), DIMENSION(:) :: FORCE, LOAD
INTEGER :: n
TYPE(Element_t), POINTER :: Element
!----------------------------------------------------------------
REAL(KIND=dp) :: Basis(n),dBasisdx(n,3),ddBasisddx(n,3,3)
REAL(KIND=dp) :: detJ, LoadAtIP,&
LocalHeatCapacity, LocalDensity
LOGICAL :: stat, getSecondDerivatives
INTEGER :: t,j
TYPE(GaussIntegrationPoints_t) :: IP
TYPE(Nodes_t) :: Nodes
SAVE Nodes
!----------------------------------------------------------------
CALL GetElementNodes( Nodes )
FORCE = 0.0d0
!Numerical integration:
!----------------------
IP = GaussPoints( Element )
!-----------------------------------------------------------------
! Loop over Gauss-points (boundary element Integration)
!-----------------------------------------------------------------
DO t=1,IP % n
!Basis function values & derivatives at the integration point:
!-------------------------------------------------------------
getSecondDerivatives = .FALSE.
stat = ElementInfo( Element, Nodes, IP % U(t), IP % V(t), &
IP % W(t), detJ, Basis, dBasisdx, ddBasisddx, &
getSecondDerivatives)
!The source term at the integration point:
CSC – IT Center for Science
18. Basic Programming 115
!-----------------------------------------
LoadAtIP = SUM( Basis(1:n) *LOAD(1:n) )
!-----------------------------------------------------------------
! Loop over j-components of matrices and over force vector
!-----------------------------------------------------------------
DO j=1,n
FORCE(j) = FORCE(j) + IP % s(t) *DetJ *LoadAtIP *Basis(j)
END DO
END DO
!----------------------------------------------------------------
END SUBROUTINE BoundaryCondition
!----------------------------------------------------------------
The Fortran 90 file is completed providing the end statement for the user solver subroutine MyHeatSolver.
!----------------------------------------------------------------
END SUBROUTINE MyHeatSolver
!----------------------------------------------------------------
18.4.3 User Defined Functions and Subroutines in Parallel
User defined functions (boundary and initial conditions as well as body forces) and routines (solvers) do not
have to be especially adopted to be run within a ElmerSolver_mpi call. Nevertheless, there might be
special situation that demand additional code-lines to avoid issues caused by parallel execution. To that end,
one has to be able to retrieve certain information on the parallel environment. Therefore, a special type exist
within Elmer with the following entries:
TYPE ParEnv_t
INTEGER :: PEs
INTEGER :: MyPE
LOGICAL :: Initialized
INTEGER :: ActiveComm
LOGICAL, DIMENSION(:), POINTER :: Active
LOGICAL, DIMENSION(:), POINTER :: IsNeighbour
LOGICAL, DIMENSION(:), POINTER :: SendingNB
INTEGER :: NumOfNeighbours
END TYPE ParEnv_t
The following list shall give an overview of the most important operations in this context:
How many parallel processes/domains are there in the current run? This can be accessed by
the construct ParEnv % PEs. For instance, the following line can be used to introduce a branch
between instructions for a serial and a parallel Elmer run:
! serial run
IF ( ParEnv % PEs <= 1 ) THEN
ELSE ! parallel run
ENDIF
What parallel process/domain is the current task being run on? This number can be inquired via
the construct ParEnv % myPe
What parallel process/domain my element belongs to? If the pointer CurrentElement is point-
ing to the actual element of interest, the host process/domain this CurrentElement % partIndex.
For instance, in combination with the instruction in the previous point, this can be used to inquire
whether an element actually is being processed within the domain (is part of the domain’s matrix) or
is just an halo-element:
CSC – IT Center for Science
18. Basic Programming 116
! skip rest of routine, because element is not
! part of the domain (halo-element)
IF (ParEnv % myPe .NE. Element % partIndex) CYCLE
Is the parallel process/domain with index ia neighbour of my current process/domain? This is
inquired by the logical value in ParEnv % IsNeighbour(i)
How many neighbour processes/domains does my current process/domain have? This is inquired
by the logical value in ParEnv % NumOfNeighbours
++++++++++++++++++++++++++++++++
18.5 Compilation and Linking of User Defined Routines/Functions
Elmer provides a script elmerf90 (residing in $ELMER_HOME/bin) for compilation of external routines.
This makes sure that the same compiler the Elmer installation has been compiled with is used for compi-
lation of the additional routines. This is necessary, as usually modules of different Fortran compilers are
incompatible. Additionally, the script contains the necessary compiler options in order to take care that the
essential libraries are being linked to the resulting shared object.
Compilation of ones own code is quite straight forward. It is essential, that the wrapper elmerf90 is in
the path of the system (or alternatively called with its path preceded). If now the routines/functions are in the
file mySolver.f90 then the dynamically linked object is obtained under UNIX/Linux with the command
elmerf90 mySolver.f90 -o mySolver
Under Windows XP, compilation can be made via the ElmerGUI or manually using the command window
(Windows XP: Start →Run, enter cmd). Preferably, the source file should be in the tree of the Elmer (e.g.,
Figure 18.6: Manual compilation of additional modules in Windows XP
C:\Elmer6.0) installation under the directory C:\Elmer6.0\share\elmersolver\lib. From
within this directory, the command
elmerf90 myOwnRoutines.f90
should create the shared object mySolver.dll that should be found automatically in every Elmer run.
CSC – IT Center for Science
Appendix A
Format of mesh files
In this appendix the format of ElmerSolver mesh files is desribed. The mesh data are arranged into four
separate files: mesh.header,mesh.nodes,mesh.elements, and mesh.boundary.
In parallel computations the names of the files are changed so that mesh is replaced by string part.n
where nrefers to the number of the partition. Additionally in parallel there is the part.n.shared file.
In the mesh files numeric codes are used for distinguishing different element types. For the element type
codes and the node numbering order of the elements see also appendix D.
A.1 The format of header file
The header file mesh.header tells how many nodes and elements are present in the mesh. The lines of
this file are organized as
nodes elements boundary-elements
nof_types
type-code nof_elements
type-code nof_elements
...
In the first line the numbers of nodes, elements, and boundary elements are given, while the count in the
second line is the number of different element types used in the mesh. The lines which follow give the
numbers of elements as sorted into different element types.
For example, the following header file
300 261 76
2
404 261
202 76
tells us that the mesh is composed of 300 nodes, 261 elements, and 76 boundary elements. Two different
element types are used in the mesh: there are 261 elements of type code 404 (bilinear quadrilateral) and 76
elements of type code 202 (linear line element).
A.2 The format of node file
The file mesh.nodes contains node data so that each line defines one node. Each line starts with two
integers followed by three real numbers:
n1 p x y z
n2 p x y z
...
nn p x y z
CSC – IT Center for Science
A. Format of mesh files 118
The first integer is the identification number for the node. The second integer is a partition index for parallel
execution and is not usually referenced by the solver in the case of sequential runs. If the partition index is
not of particular use, it may be set to be -1 (or 1). The real numbers are the spatial coordinates of the node.
Three coordinates should always be given, even if the simulation was done in 1D or 2D. It is also noted that
the nodes may be listed in any order.
A.3 The format of element file
The mesh.elements file contains element data arranged as
e1 body type n1 ... nn
e2 body type n1 ... nn
...
en body type n1 ... nn
Each line starts with an integer which is used for identifying the element. The integer body defines the
material body which this element is part of. Then the element type code and element nodes are listed. For
example, the element file might start with the following lines:
1 1 404 1 2 32 31
2 1 404 2 3 33 32
3 1 404 3 4 34 33
4 1 404 4 5 35 34
...
A.4 The format of boundary element file
The elements that form the boundary are listed in the file mesh.boundary. This file is similar to the
mesh.elements file and is organized as
e1 bndry p1 p2 type n1 ... nn
e2 bndry p1 p2 type n1 ... nn
...
en bndry p1 p2 type n1 ... nn
The first integer is again the identification number of the boundary element. Next the identification number
of the part of the boundary where this element is located is given. Whether the boundary element can be
represented as the side of a parent element defined in the file mesh.elements is indicated using the two
parent element numbers p1 and p2. If the boundary element is located on an outer boundary of the body,
it has only one parent element and either of these two integers may be set to be zero. It is also possible that
both parent element numbers are zeros. Finally the element type code and element nodes are listed.
A.5 The format of shared nodes file
The nodes that are shared between different partitions are written in the part.n.shared file. The file is
organized as
n1 np p1 p2 ... pm
n2 np p1 p2 ... pm
...
nn np p1 p2 ... pm
Here the first index on the row refers to a node that is shared by different partitions, the partition under study
being one of the owners. Integer np gives the number of partitions sharing the node and this is followed
by the number of the sharing partitions. By convention the owner partition is the first one of the partitions.
Most shared nodes are shared just by two partitions but the can be also more owners on some nodes.
CSC – IT Center for Science
A. Format of mesh files 119
A.6 Exceptions on parallel mesh format
When the mesh is saved in parallel there is an exception that must be dealt with. Assume that you have a
boundary condition for which you need to assign Dirichlet conditions. If it would happen to be so that some
partition shares such boundary node without being the owner of the boundary element itself, the boundary
nodes must still somehow be available for setting boundary conditions. In practice this is done so that new
single-node elements of type 101 are added to the mesh.n.boundary list. These must be given the same
bndry index as the initial boundary that it was deduced from.
CSC – IT Center for Science
Appendix B
Format of result output files
B.1 Format versions
Result files can be written as either ASCII text or in binary. This is controlled by the parameter
Binary output = logical true|false
in the ’Simulation’ section of the .sif file. Default is false.
The format of the file is recorded on it’s first line;1it’s either
BINARY v.e
or
ASCII v
The vat denotes the version number of the format, and the ein the binary format denotes an endianess-
marker; either Lfor little endian or Bfor big endian.
ElmerSolver can read files of older format versions for restarting, but all new files are written in the
newest formats. The current formats documented here are ASCII version 1 and BINARY version 2.
B.2 General structure
Both binary and ASCII files have the following general structure. In the binary files, the header is separated
from the rest by a null byte. The ASCII format has no such separator.
[File format version line]
[Header]
[<null byte> (binary only)]
[timestep 1]
[timestep 2]
[timestep 3]
.
.
.
[timestep n]
1except for old ASCII files, that lack the format version line, and start with ”!File started at:”
CSC – IT Center for Science
B. Format of result output files 121
B.2.1 The header
The header looks the same for both binary and ASCII (ans is written in ASCII also for binary files):
!File started at: date time
Degrees of freedom:
variable 1 n1 :fs
variable 2 n2 :fs
variable 3 n3 :fs
.
.
.
Total DOFs: nTotDOFs
Number Of Nodes: nNodes
Note that in the list of variables, vectors appear both as vectors (DOF > 1) and separate components (DOF =
1).
B.2.2 Timesteps
For each time step, the time and the values of the variables are stored. For vector variables, the components
are stored as separate variables, typically named varname 1,varname 2, etc.
If the parameter
Omit unchanged variables in output = logical true|false
in the ’Simulation’ section of the .sif file is set to true, only variables whose values has changes since last
time are saved. Default is false
For the binary format, the following type specifiers are used in this document:
<s(str)> Null terminated string of characters.
<i(i4)> 4 byte integer.
<j(i8)> 8 byte integer.
<x(dbl)> 8 byte floating point number.
For this part of the file, the general structure of binary and ASCII files are essantially the same, with just
some minor differences:
CSC – IT Center for Science
B. Format of result output files 122
ASCII Binary
Time: ns nt t<Time:(str)><ns(i4)><nt(i4)><t(dbl)>
Variablename_1 <Variablename_1(str)>
[Permutation table p1] [Permutation table p1]
.
.
..
.
.Variable 1’s values
∀is.t. p1(i)>0.
Var1(p1(i)) <Var 1(p1(i))(dbl)>
.
.
..
.
.
Variablename_2 <Variablename_2(str)>
[Permutation table p2] [Permutation table p2]
.
.
..
.
.Variable 2’s values
∀is.t. p2(i)>0.
Var2(p2(i)) <Var 2(p2(i))(dbl)>
.
.
..
.
.
Variablename_3 <Variablename_3(str)>
[Permutation table p3] [Permutation table p3]
.
.
..
.
.
.
.
..
.
.Variable 3’s values
∀is.t. p3(i)>0.
Var3(p3(i)) <Var 3(p3(i))(dbl)>
.
.
..
.
.
nt= time step number, ns= saved time steps number, t= time.
The permutation tables
The permutation tables can be stored in three different ways:
1. As an explicit table:
ASCII Binary
Perm: size np <size(i4)><np(i4)>
.
.
..
.
.Permutation indexes i
and values p(i)
∀is.t. p(i)>0.
i p(i) <i(i4)><p(i)(i4)>
.
.
..
.
.
size = total size of the permutation table (> 0), and np = number of positive values in the table.
2. If the permutation table is the same as for the previous variable, there’s no need to store it again. This
case is written as
ASCII Binary
Perm: use previous <−1(i4)><Pos(i8)>
Pos in the binary format is the position in bytes of the previous table.
3. No permutation table; corresponds to the case
size =np =nNodes,and p(i) = i∀i.
This case is stored as
ASCII Binary
Perm: NULL <0(i4)>
CSC – IT Center for Science
B. Format of result output files 123
B.3 The positions file
For binary format, a positions file named ’outputfile.pos’ will be created. It contains the positions (in bytes)
of the timesteps and variables in the result file, expressed as 8 byte integers. It has the following format
(nVar = number of variables):
<Endianess-marker (char)>
<nVar(i4)>
<varname 1(str)>
<varname 2(str)>
.
.
.
<varname nVar(str)>
<Pos. for Timestep 1(i8)>
<Pos. for variable 1(i8)>
<Pos. for variable 2(i8)>
.
.
.
<Pos. for variable nVar(i8)>
<Pos. for Timestep 2(i8)>
<Pos. for variable 1(i8)>
<Pos. for variable 2(i8)>
.
.
.
<Pos. for variable nVar(i8)>
<Pos. for Timestep 3(i8)>
<Pos. for variable 1(i8)>
<Pos. for variable 2(i8)>
.
.
.
<Pos. for variable nVar(i8)>
.
.
.
Note: Positions are saved for all variables for every time step; even if an unchanged variable isn’t saved
to the result file for a time step, it will still have a position in the .pos file (i.e. the position of where it was
saved last time). Because of this all timesteps has the same size of (nVar +1)×8bytes. Finding the position
of the n:th variable of the m:th time step is therefore easy; it’s found at the
(size-of-header + ((nVar + 1) ×(m−1) + n)×8) : th
byte in the positions file.
CSC – IT Center for Science
Appendix C
Format of ElmerPost Input File
The lines of ElmerPost input file are organized as
nn ne nf nt scalar: name vector: name ...
x0 y0 z0
... ! nn rows of node coordinates (x,y,z)
xn yn zn
group-name element-type i0 ... in
... ! group data and element decriptions
group-name element-type i0 ... in
#time 1 timestep1 time1
vx vy vz p ...
... ! nn rows
vx vy vz p ...
#time 2 timestep2 time2
vx vy vz p
... ! nn rows
vx vy vz p ...
....
#time n timestepn timen
vx vy vz p ...
... ! nn rows
vx vy vz p ...
The header
The file starts with the header line which contains the following information:
•nn: the total number of nodes
•ne: the total number of the elements including boundary elements
•nf: the total number of degrees of freedom, i.e. the total number of scalar unknowns in the model
•nt: the number of time steps for which solution data is stored
•scalar: name vector: name ... : the list which pairs variable names with their types.
The mesh and group data
After the header the node coordinates are given, each coordinate triplet on its own row. Three coordinates
shoud be given even if the model was two-dimensional.
CSC – IT Center for Science
C. Format of ElmerPost Input File 125
Group data consist of the following information:
•group-name: the name of the element group (having the same material, body etc.)
•element-type: the numeric code giving the element type; see also Appendix D.
•The numbers i0 ... in are the indeces of the element nodes. The nodes are referenced using the
row indeces of the node coordinate array at the beginning of the file The first node in the array has the
index zero.
It is noted that there is also another level of element grouping that can be employed as follows
#group group1
element-definitions
...
#group group2
element-definitions
...
#endgroup group2
element-definitions
...
#endgroup group1
The number of element groups is not restricted in any way.
The solution data
For each timestep the following solution data is written:
•#time n t time:nis the order number of the solution data set, tis the simulation timestep
number, and time is the current simulation time.
•The next nn rows give the node values of the degrees of freedom. The values are listed in the same
order as given in the header with the keywords scalar: and vector:
An example file
Here a very simple example file is given. There is only one element, three nodes, one variable, and the
solution data are given for a single timestep:
3 1 1 1 scalar: Temperature
0 0 0
1 0 0
0 1 0
#group all
body1 303 0 1 2
#endgroup all
#time 1 1 0
1
2
3
CSC – IT Center for Science
Appendix D
Basic element types
The basic nodal element types which ElmerSolver can handle are the linear and quadratic elements in one,
two, and three dimensions. In addition, cubic elements in one and two dimensions are also available. Each
element type is associated with a unique code which consists of the element family (the count of hundreds)
and the number of nodes.
•nodal element (101)
•linear (element type code 202), quadratic (203) and cubic (204) line segments
•linear (303), quadratic (306) and cubic (310) triangles, see Figure D.1
u
v
1 2
3
u
v
1 2
3
4
56
Figure D.1: The linear (303) and quadratic (306) triangular elements.
•bilinear (404), quadratic (408,409) and cubic (412,416) quadrilaterals, see Figure D.2
•linear (504) and quadratic (510) tetrahedrons, see Figure D.3
•linear (605) and quadratic (613) pyramids
•linear (706) and quadratic (715) wedges
•trilinear (808) and quadratic (820,827) hexahedrons (or bricks), see Figure D.4.
The element types are defined in the file elements.def and new element types may be added by adding
their definitions to this file.
CSC – IT Center for Science
D. Basic element types 127
u
v
(0,0)
1 2
34
u
v
(0,0)
(−1,0) (1,0)
(0,−1)
(0,1)
1 2
34
5
6
7
8
Figure D.2: The four-node (404) and eight-node (408) quadrilateral elements.
u
w
v
1 2
4
3
u
w
v
1 2
4
3
5
98
7 6
10
Figure D.3: The linear (504) and quadratic (510) tetrahedron elements.
CSC – IT Center for Science
D. Basic element types 128
u
w
v
(0,0,0)
1
2
6
5
4
3
7
8
u
w
v
1
2
6
5
4
3
7
8
9
14
17
13
12
10
18
20
11
15
19
16
Figure D.4: The 8-node (808) and 20-node (820) hexahedron elements.
CSC – IT Center for Science
Appendix E
Higher-order finite elements
Higher-order finite elements are elements for which the degree of basis functions is higher than 1. If better
solution accuracy is sought by increasing the approximation degree, higher-order elements are also called p-
elements. They differ from the usual Lagrange-type elements in that, in addition to nodal basis functions of
the lowest degree, there exists additional basis functions which are associated with edges, faces and interiors
of elements.
•Side modes get their values along some element edge. They vanish towards other edges and all nodal
points of the element. Side modes are defined for all 2d and 3d elements.
•Face modes get their values along some element face. They vanish towards other faces and all edges
and nodal points of the element. Face modes are only defined for 3d elements.
•Internal modes get their values inside element and vanish towards element faces, edges and nodal
points. They are defined for all 1d, 2d and 3d elements.
Thus, for p-elements the usual element interpolant, now denoted as uh,p, is the sum of nodal, edge, face
and bubble interpolants as
uh,p =uv
h,p +ue
h,p +uf
h,p +ub
h,p (E.1)
where uv
h,p is the lowest-order nodal interpolant while ue
h,p,uf
h,p and ub
h,p are edge, face and bubble inter-
polants, respectively. Contribution of one p-element to the global system is similar to that of h-element.
Naturally for higher-order elements the size of the local stiffness matrix is larger, because of the larger
number of basis functions.
Properties for good p-elements include computational efficiency, at least partial orthogonality and hierar-
chy of basis functions. Hierarchy means that if basis for some element of a given degree pis Bpit holds for
p+1 that Bp⊂ Bp+1. Orthogonality properties of basis functions ensure that condition number of the global
stiffness matrix does not increase as dramatically as for the nodal (Lagrange) elements of higher order. This
ensures good numerical stability. Some good references to higher-order finite elements in literature are [3]
by Szabo and Babuska and [4] by Solin et al.
E.1 Higher-order elements in Elmer
Elements implemented in Elmer follow the ones presented in [3]. Now let us define some orthogonal poly-
nomials based on Legendre polynomials Pi(x), i ≥0. The so-called Lobatto shape functions φkare defined
as
φk(ξ) = s1
2(2k−1)(Pk(ξ)−Pk−2(ξ)), k = 2,3,... (E.2)
where Pkare Legendre polynomials. Function φk(ξ)has two of its roots at ξ=±1, so define now another
function ϕkas
CSC – IT Center for Science
E. Higher-order finite elements 130
ϕk(ξ) = 4φk(ξ)
1−ξ2, k = 2,...,p (E.3)
Functions φkand ϕkare used to define higher order elements. Different element shapes and their basis
functions are defined in Appendix E.3. The pyramidal element used in Elmer is based loosely to Devloo’s
representation in [2].
In Elmer elements of varying polynomial degree pmay be used in the same mesh. The mesh may also
consist of different element types, as the side and face modes for different element types are compatible with
one another. Pyramidal and wedge elements of high order to connect tetrahedral and brick elements are also
supported. To achieve the best possible converge the use of pyramidal elements in a mesh should be kept to
a minimum.
The global continuity of the higher-order finite element approximation is enforced by the solver when
the function ElementInfo is used for obtaining basis functions values for elements. To combine elements
of varying degree a maximum rule is used. Thus, if two or more elements share an edge and have different
polynomial degrees, the maximum of among the values is chosen as the degree of the global edge.
The approximation degree can be specified either in the solver section or in the equation section of the
solver input file. The solver section may for example contain the command
Element = "p:4"
so that the spatial discretization is based on polynomials of degree 4. Basically the same command may also
be used in the equation section if the intention is to apply the same basis functions for all active solvers,
while a more advanced element definition command is needed if the user wants to specify the approximation
degree solver-wise. One may define for example
Equation 1
Active Solvers(2) = 1 2
Element{1} = String "p:4"
Element{2} = String "p:1"
End
so that the fourth-order approximation is applied in connection with the solver having the integer identifier
1, while the first-order approximation is used for the solver with the identifier 2.
A body-wise changing approximation degree can also be defined by a combined use of the body and
equation sections. If we have for example split the entire computational domain into two bodies and we
want to use basis functions of degree p= 1 in the body 1 and the polynomial degree p= 3 in the body 2,
we may make the desired element choices by including the following definitions in the sif-file:
Body 1
Equation = 1
Material = 1
End
Body 2
Equation = 2
Material = 1
End
Equation 1
Active Solvers(2) = 1 2
Element{1} = String "p:1"
Element{2} = String "p:1"
End
CSC – IT Center for Science
E. Higher-order finite elements 131
Equation 2
Active Solvers(2) = 1 2
Element{1} = String "p:3"
Element{2} = String "p:3"
End
It is worth noting, however, that all available solvers of Elmer may not support p-approximation without
modifying the code. Usually this only consists of making local stiffness matrix and force vector larger.
Nevertheless several solvers already support the approximation with high-order polynomials.
The actual number of degrees of freedom associated with the element edges, faces or interiors depends
on the element type and the local polynomial degree p. Each degree of freedom is naturally associated with
some basis function. The following table gives the number of degrees of freedom for elements used in Elmer:
Element Nodes Edges Faces Bubbles
Line 2- - p−1
Quadrilateral 4 4(p−1) -(p−2)(p−3)
2
Triangle 3 3(p−1) -(p−1)(p−2)
2
Brick 8 12(p−1) 3(p−2)(p−3) (p−3)(p−4)(p−5)
6
Tetrahedron 4 6(p−1) 2(p−1)(p−2) (p−1)(p−2)(p−3)
6
Wedge 6 9(p−1) -(p−2)(p−3)(p−4)
6
(quad. face) - - 3(p−2)(p−3)
2-
(triang. face) - - (p−1)(p−2) -
Pyramid 5 8(p−1) -(p−1)(p−2)(p−3)
6
(quad. face) - - (p−2)(p−3)
2-
(triang. face) - - 2(p−1)(p−2) -
Ap-element passed to the Elmer Gaussian point generator GaussPoints contained in the module
Integration returns enough integration points to integrate the worst-case product of two element basis
functions. Here the worst case is the integration for two basis functions of degree pm= max{pe, pf, pb}.
As the Gaussian quadrature is accurate up to the degree p= 2n−1, where nis the number of points used,
the number of points for each element is calculated from
n=2pm+ 1
2(E.4)
and rounded up to the nearest integer. To get the final number of points for multiple integrals, nis raised to
the power of element dimension. If integral includes a non-constant factor, i.e RKαφiφjdK where αis a
function of degree k, the default quadrature is not accurate and the number of integration points needs to be
set manually. Now the minimum number of gaussian points to integrate accurately becomes
n=min {2pm+k, 3pm}+ 1
2(E.5)
which may again be rounded up to the nearest integer and raised to the power of element dimension to get
the actual number of integration points.
E.1.1 Boundary conditions
From the point of view of programming, nonhomogeneous boundary conditions of Dirichlet type, e.g. u=g
on a boundary, are more difficult to handle for higher order elements. Even though the nodal values can be
set easily, finding coefficients of additional basis functions such that the data is represented (approximately)
as a linear combination of the basis functions over the element boundary is not as straightforward. The
subroutine DefaultDirichletBCs solves the unknown coefficients of higher-order basis functions by
seeking the best approximation of the data over the boundary. Thus, from the viewpoint of the user, Dirichlet
boundary conditions need no extra actions as compared to the use of the conventional elements. Neumann
CSC – IT Center for Science
E. Higher-order finite elements 132
and Robin boundary conditions can also be treated in the same way as done in the case of the standard
Lagrange elements.
E.1.2 Some practical aspects
Generally using p-elements yields a better approximation than using normal linear elements. In fact, conver-
gence for p-elements is exponential when there are no singularities inside or on the boundary of the solution
domain. Typical singular points in the solution are points where boundary conditions or material parameters
change abruptly or vertex-type singularities (such as the inner corner of a L-shaped beam or a crack tip).
When there are singular points inside the domain, convergence is algebraic. If the singular point is a nodal
point, convergence rate is twice that of h-method, otherwise it is equal to the h-method.
Therefore, if possible, the mesh should be designed such that near nodal singularities small low-degree
(p= 1) elements were used. In other parts of the solution domain, where the solution is smoother, large
elements with high polynomial degree are adviced. As Elmer is not hp-adaptive so that the element degree
is not modified by the solver, mesh design issues must be taken into account for computational efficiency.
It is well known that it is possible to reduce the size of the linear system by the elimination of bubble
functions. This procedure is often called the static condensation. In Elmer the static condensation for the
local stiffness matrix may be done (and is adviced to be done) when the bubble part need not be retrieved
afterwards (for example time stepping algorithms need to know the complete solution). Whether the bub-
ble degrees of freedom are assembled into the global system is controlled by the keyword Bubbles in
Global System.
E.2 ElmerSolver services for higher-order elements
This section describes some of the services related to p-elements, which are included in different parts of the
Solver.
E.2.1 Properties of p-elements
There are several utilities for determining p-element properties. First of all, it is possible to check if some
element is a p-element by checking the isPElement flag of the element. If the flag is set to be true,
element is a p-element. Functions
isPTriangle( Element )
isPTetra( Element )
isPPyramid( Element )
isPWedge( Element )
check if the given element is a p-type triangle, tetrahedron, pyramid or wedge. They are implemented
because the reference p-element for triangles, tetrahedra, pyramids and wedges are different from those
defined for the Lagrange elements. The routines
getEdgeP( Element, Mesh )
getFaceP( Element, Mesh )
return the maximum polynomial degrees associated with the element edges or faces when Element and
finite element mesh Mesh are given.
E.2.2 Fields related to p-elements
In the module Types, the type Element_t has the following p-element related fields
INTEGER :: TetraType
LOGICAL :: isPElement
LOGICAL :: isEdge
CSC – IT Center for Science
E. Higher-order finite elements 133
INTEGER :: localNumber
INTEGER :: GaussPoints
Tetratype defines the type of the tetrahedral p-element: for non-tetrahedral elements Tetratype=0,
for tetrahedral elements Tetratype={1,2}.
isPElement defines if an element is of higher order: isPElement=.TRUE. for p-elements, other-
wise .FALSE..
isEdge defines if the element is an edge for some higher entity, i.e. an edge or a face of a 2d or 3d
element. If isEdge=.TRUE., the element is an edge, .FALSE. otherwise.
For boundary elements localNumber defines which local edge or face number the boundary element
has with respect to its parent element.
GaussPoints defines the number nof the Gauss points for the element. The value is calculated from
n= (2pm+1
2)d, where dis the element dimension and pmthe maximum polynomial degree. nis rounded
up to the nearest integer. The value GaussPoints corresponds to the requirement that a polynomial of
degree 2pmis integrated accurately.
When modifying a solver to support higher-order elements, the maximum number of degrees of freedom
associated with an element may be obtained from the mesh variable MaxElementDOFs. This variable is
set by the read-in process of the mesh.
E.2.3 Higher-order basis and element mappings
Basis functions for higher-order elements are defined in the module PElementBase. The module contains
also definitions for the functions φkand ϕkand Legendre polynomials. These definitions have been gener-
ated into explicit form with the symbolic program Maple [1] up to pmax = 20. This means that no recursion
is needed for the generation of the values of Legendre polynomials or other lower-level components based
on them, if p < pmax.
Generally higher-order basis functions take as their arguments the point in which to calculate the function
value and indexing i,m(i, j)or m(i, j, k)depending on the function type. All edge functions take in addition
to these parameters a special optional flag, namely invertEdge, which defines if the direction of the edge
basis function needs to be inverted. In Elmer the parametrization of edges depends on the global node
numbering. To describe this, let Aand Bbe the global node numbers of an edge. The varying parameter
of edge function then varies between [−1,1] as A→B. Inversion is then used for enforcing the global
continuity of edge basis functions which are not properly aligned. The edge rule is presented in figure E.2.3.
BA
Figure E.1: The global direction of edge for global node indexes A < B
Most of the face functions take an optional orientation argument in order to specify which face functions
are formed. This is formed according to the global numbers of face nodes. There are rules for triangular
and square faces. Let A, B, C be the global nodes of a triangular face. Globally the face is aligned so that
A < B < C. For square faces let A= min{vi}where viare the global node indeces of square face
and B, C correspond to the nodes next to the node Aon the face. The square face is aligned by the rule
A < B < C for these nodes. These rules are presented in figure E.2.3.
The tetrahedral element uses an exception to the above interface rules, as the edge and face functions
of tetrahedral elements take the type of the tetrahedral element as their optional argument. This is due to
fact that it is possible to reduce any tetrahedral element to one of the two reference tetrahedral elements for
which all edges and faces are defined so that their local orientation matches global orientation. This means,
that for tetrahedral elements, global continuity does not need to be enforced, if proper reduction to one of
the two reference elements has been made.
Mappings from element nodal numbers to describe different p-element edges or faces are defined in
the module PElementMaps. These mappings generally define which nodes of the element belong to
a certain local edge or face of the element. Mappings to elements edges, faces and from faces to local
CSC – IT Center for Science
E. Higher-order finite elements 134
B
C
AA B
C D
Figure E.2: The global direction of triangular and quadrilateral faces for global node indexes A < B < C;
Ahas the lowest index among the indeces of the face.
edge numbers may be obtained from the routines GetElementEdgeMap,GetElementFaceMap and
GetElementFaceEdgeMap. The mappings may also be accessed via methods getTePeMap, where Te
is the element name and Pe={Edge,Face}is the part of the element for which the map is requested. The
routine getElementBoundaryMap returns mappings for element boundaries depending on the element
type.
For example, to get the global nodes for the brick face number 4, one would use the following Fortran90
code
map(1:4) = getBrickFaceMap(4)
nodes(1:4) = Element % NodeIndexes(map)
E.3 Basis functions
Let λ1, λ2, λ3∈ {±ξ, ±η, ±ζ}and additionally Tiλi=φ.
E.3.1 Line
ξ
−1 1
21
Figure E.3: Line element
Nodal basis
L1=1−ξ
2
L2=1 + ξ
2
Bubble basis
L(0)
i=φi(ξ), i = 2,...,p
CSC – IT Center for Science
E. Higher-order finite elements 135
E.3.2 Quadrilateral
−1
−1
1
1
v2
v3
v4
v1
ξ
η
Figure E.4: Quadrilateral element
Nodal basis
N1=1
4(1 −ξ)(1 −η)
N2=1
4(1 + ξ)(1 −η)
N3=1
4(1 + ξ)(1 + η)
N4=1
4(1 −ξ)(1 + η)
Edge basis
N(1,2)
i=1
2(1 −η)φi(ξ), i = 2,...,p
N(2,3)
i=1
2(1 + ξ)φi(η), i = 2,...,p
N(4,3)
i=1
2(1 + η)φi(ξ), i = 2,...,p
N(1,4)
i=1
2(1 −ξ)φi(η), i = 2,...,p
Bubble basis
N(0)
m(i,j)=φi(ξ)φj(η)
where i, j ≥2, i +j= 4, . . . , p
CSC – IT Center for Science
E. Higher-order finite elements 136
−1 1
0
ξ
v3
η√3
v2
v1
Figure E.5: Triangle element
E.3.3 Triangle
Nodal basis
L1=1
2(1 −ξ−1
√3η)
L2=1
2(1 + ξ−1
√3η)
L3=η
√3
Edge basis
N(1,2)
i=L1L2ϕi(L2−L1), i = 2,...,p
N(2,3)
i=L2L3ϕi(L3−L2), i = 2,...,p
N(3,1)
i=L3L1ϕi(L1−L3), i = 2,...,p
Bubble basis
N(0)
m(j,n)=L1L2L3P1(L2−L1)jP1(2L3−1)n
where j, n = 0,...,i−3,j+n=i−3, i = 3,...,p
CSC – IT Center for Science
E. Higher-order finite elements 137
v2
v3
v7
v8
v4
v6
v1
v5
ξ
ζ
η
Figure E.6: Brick element
E.3.4 Brick
Nodal basis
N1=1
8(1 −ξ)(1 −η)(1 −ζ)
N2=1
8(1 + ξ)(1 −η)(1 −ζ)
N3=1
8(1 + ξ)(1 + η)(1 −ζ)
N4=1
8(1 −ξ)(1 + η)(1 −ζ)
N5=1
8(1 −ξ)(1 −η)(1 + ζ)
N6=1
8(1 + ξ)(1 −η)(1 + ζ)
N7=1
8(1 + ξ)(1 + η)(1 + ζ)
N8=1
8(1 −ξ)(1 + η)(1 + ζ)
CSC – IT Center for Science
E. Higher-order finite elements 138
Edge basis
N1,2
i−1=1
4φi(ξ)(1 −η)(1 −ζ)
N2,3
i−1=1
4φi(η)(1 + ξ)(1 −ζ)
N4,3
i−1=1
4φi(ξ)(1 + η)(1 −ζ)
N1,4
i−1=1
4φi(η)(1 −ξ)(1 −ζ)
N1,5
i−1=1
4φi(ζ)(1 −ξ)(1 −η)
N2,6
i−1=1
4φi(ζ)(1 + ξ)(1 −η)
N3,7
i−1=1
4φi(ζ)(1 + ξ)(1 + η)
N4,8
i−1=1
4φi(ζ)(1 −ξ)(1 + η)
N5,6
i−1=1
4φi(ξ)(1 −η)(1 + ζ)
N6,7
i−1=1
4φi(η)(1 + ξ)(1 + ζ)
N8,7
i−1=1
4φi(ξ)(1 + η)(1 + ζ)
N5,8
i−1=1
4φi(η)(1 −ξ)(1 + ζ)
Face basis
N(1,2,5,6)
m(i,j)=1
2(1 −η)φi(ξ)φj(ζ)
N(1,2,4,3)
m(i,j)=1
2(1 −ζ)φi(ξ)φj(η)
N(1,4,5,8)
m(i,j)=1
2(1 −ξ)φi(η)φj(ζ)
N(4,3,8,7)
m(i,j)=1
2(1 + η)φi(ξ)φj(ζ)
N(5,6,8,7)
m(i,j)=1
2(1 + ζ)φi(ξ)φj(η)
N(2,3,6,7)
m(i,j)=1
2(1 + ξ)φi(η)φj(ζ)
where i, j = 2,3,...,p−2,i+j= 4,5,...,p
Bubble basis
N(0)
m(i,j,k)=φi(ξ)φj(η)φk(ζ)
where i, j, k = 2,3,...,p−4,i+j+k= 6,7,...,p
CSC – IT Center for Science
E. Higher-order finite elements 139
v1
v4
v3
v2
ζ
ξ
η
ζ
ξ
η
v1
v4
v3
v2
Figure E.7: Tetrahedral elements of types 1 and 2
E.3.5 Tetrahedron
Nodal basis
L1=1
2(1 −ξ−1
√3η−1
√6ζ)
L2=1
2(1 + ξ−1
√3η−1
√6ζ)
L3=√3
3(η−1
√8ζ)
L4=r3
8ζ
Edge basis
Type 1
N(1,2)
i−1=L1L2ϕi(L2−L1), i = 2,...,p
N(1,3)
i−1=L1L3ϕi(L3−L1), i = 2,...,p
N(1,4)
i−1=L1L4ϕi(L4−L1), i = 2,...,p
N(2,3)
i−1=L2L3ϕi(L3−L2), i = 2,...,p
N(2,4)
i−1=L2L4ϕi(L4−L2), i = 2,...,p
N(3,4)
i−1=L3L4ϕi(L4−L3), i = 2,...,p
Type 2
N(3,2)
i−1=L3L2ϕi(L2−L3), i = 2,...,p
Edges (1,2),(1,3),(1,4),(2,4) and (3,4) according to type 1.
Face basis
Type 1
CSC – IT Center for Science
E. Higher-order finite elements 140
N(1,2,3)
m(i,j)=L1L2L3Pi(L2−L1)Pj(2L3−1)
N(1,2,4)
m(i,j)=L1L2L4Pi(L2−L1)Pj(2L4−1)
N(1,3,4)
m(i,j)=L1L4L3Pi(L3−L1)Pj(2L4−1)
N(2,3,4)
m(i,j)=L2L3L4Pi(L3−L2)Pj(2L4−1)
Type 2
N(1,3,2)
m(i,j)=L1L3L2Pi(L3−L1)Pj(2L2−1)
N(3,2,4)
m(i,j)=L3L2L4Pi(L2−L3)Pj(2L4−1)
where i, j = 0,1,2,...,p−3,i+j= 0,1,...,p−3. Faces (1,2,4) and (1,3,4) defined according to type
1.
Bubble basis
N(0)
m(i,j,k)=L1L2L3L4Pi(L2−L1)Pj(2L3−1)Pk(2L4−1)
where i, j, k = 0,1,...,p−4,i+j+k= 0,1,...,p−4
E.3.6 Pyramid
ξ
η
ζ
v1v2
v5
v3
v4
Figure E.8: Pyramidal element
Nodal basis
T0(c, t) = (1 −t
√2)−c
2(1 −t
√2)
T1(c, t) = (1 −t
√2) + c
2(1 −t
√2)
CSC – IT Center for Science
E. Higher-order finite elements 141
P1=T0(ξ, ζ)T0(η, ζ)(1 −ζ
√2)
P2=T1(ξ, ζ)T0(η, ζ)(1 −ζ
√2)
P3=T1(ξ, ζ)T1(η, ζ)(1 −ζ
√2)
P4=T0(ξ, ζ)T1(η, ζ)(1 −ζ
√2)
P5=1
√2ζ
Edge basis
P(1,2)
i−1=P1(ξ, η, ζ)P2(ξ, η, ζ)ϕi(ξ)
P(2,3)
i−1=P2(ξ, η, ζ)P3(ξ, η, ζ)ϕi(η)
P(4,3)
i−1=P4(ξ, η, ζ)P3(ξ, η, ζ)ϕi(ξ)
P(1,4)
i−1=P1(ξ, η, ζ)P4(ξ, η, ζ)ϕi(η)
P(1,5)
i−1=P1(ξ, η, ζ)P5(ξ, η, ζ)ϕi(ξ
2+η
2+ζ
√2)
P(2,5)
i−1=P2(ξ, η, ζ)P5(ξ, η, ζ)ϕi(−ξ
2+η
2+ζ
√2)
P(3,5)
i−1=P3(ξ, η, ζ)P5(ξ, η, ζ)ϕi(−ξ
2−η
2+ζ
√2)
P(4,5)
i−1=P4(ξ, η, ζ)P5(ξ, η, ζ)ϕi(ξ
2−η
2+ζ
√2)
Face basis
Square face
P(1,2,3,4)
m(i,j)=P1(ξ, η, ζ)P3(ξ, η, ζ)ϕi(ξ)ϕj(η)
where i, j = 2,...,p−2,i+j= 4,...,p.
Triangular faces
P(1,2,5)
m(i,j)=P1(ξ, η, ζ)P2(ξ, η, ζ)P5(ξ, η, ζ)Pi(P2(ξ, η, ζ)−P1(ξ, η, ζ))Pj(2P5(ξ, η, ζ)−1)
P(2,3,5)
m(i,j)=P2(ξ, η, ζ)P3(ξ, η, ζ)P5(ξ, η, ζ)Pi(P3(ξ, η, ζ)−P2(ξ, η, ζ))Pj(2P5(ξ, η, ζ)−1)
P(3,4,5)
m(i,j)=P3(ξ, η, ζ)P4(ξ, η, ζ)P5(ξ, η, ζ)Pi(P4(ξ, η, ζ)−P3(ξ, η, ζ))Pj(2P5(ξ, η, ζ)−1)
P(4,1,5)
m(i,j)=P4(ξ, η, ζ)P1(ξ, η, ζ)P5(ξ, η, ζ)Pi(P1(ξ, η, ζ)−P4(ξ, η, ζ))Pj(2P5(ξ, η, ζ)−1)
where i, j = 0,...,p−3,i+j= 0,...,p−3and Pi, PjLegendre polynomials.
CSC – IT Center for Science
E. Higher-order finite elements 142
Bubble basis
P(0)
m(i,j,k)=P1(ξ, η, ζ)P3(ξ, η, ζ)P5(ξ, η, ζ)Pi(ξ
1−ζ
√2
)Pj(η
1−ζ
√2
)Pk(ζ
√2)
where i, j, k = 0,...,p−4,i+j+k= 0....,p−4and Pi, Pj, PkLegendre polynomials
E.3.7 Wedge
ζ
ξ
η
v2
v5
v1
v4
v6
v3
Figure E.9: Wedge element
Nodal basis
L1=1
2(1 −ξ−1
√3η)
L2=1
2(1 + ξ−1
√3η)
L3=√3
3η
H1=1
2L1(1 −ζ)
H2=1
2L2(1 −ζ)
H3=1
2L3(1 −ζ)
H4=1
2L1(1 + ζ)
H5=1
2L2(1 + ζ)
H6=1
2L3(1 + ζ)
CSC – IT Center for Science
BIBLIOGRAPHY 143
Edge basis
H(1,2)
i−1=1
2L1L2ϕi(L2−L1)(1 −ζ)
H(2,3)
i−1=1
2L2L3ϕi(L3−L2)(1 −ζ)
H(3,1)
i−1=1
2L3L1ϕi(L1−L3)(1 −ζ)
H(4,5)
i−1=1
2L4L5ϕi(L5−L4)(1 + ζ)
H(5,6)
i−1=1
2L5L6ϕi(L6−L5)(1 + ζ)
H(6,4)
i−1=1
2L6L4ϕi(L4−L6)(1 + ζ)
H(1,4)
i−1=L1φi(ζ)
H(2,5)
i−1=L2φi(ζ)
H(3,6)
i−1=L3φi(ζ)
Face basis
Triangular faces
H(1,2,3)
m(i,j)=1
2(1 −ζ)Pi(L2−L1)Pj(2L3−1)L1L2L3
H(4,5,6)
m(i,j)=1
2(1 + ζ)Pi(L2−L1)Pj(2L3−1)L1L2L3
where i, j = 0,1,...,p−3,i+j= 0,1,...,p−3and Pi, PjLegendre polynomials.
Square faces
H(1,2,5,4)
m(i,j)=ϕi(L2−L1)φj(ζ)L1L2
H(2,3,6,5)
m(i,j)=ϕi(L3−L2)φj(ζ)L2L3
H(3,1,4,6)
m(i,j)=ϕi(L1−L3)φj(ζ)L3L1
where i, j = 2,...,p−2,i+j= 4,...,p.
Bubble basis
H(0)
m(i,j,k)=φk(ζ)L1L2L3Pi(L2−L1)Pj(2L3−1)
where i, j = 0,...,p−5,k= 2,...,p−3,i+j+k= 2,...,p−3.
Bibliography
[1] Maplesoft home page. http://www.maplesoft.com/.
[2] P.R.B. Devloo. On the definition of high order shape functions for finite elements. Available online:
http://www.fec.unicamp.br/ phil/downloads/shape.zip.
[3] B. Szabo and I. Babuska. Finite Element Analysis. John Wiley & Sons Ltd., 1991.
[4] P. Šolin et al. Higher-Order Finite Element Methods. Chapman & Hall / CRC, 2004.
CSC – IT Center for Science
Appendix F
Face and edge elements
In Elmer the way to define a high-degree finite element discretization is based on the idea that a background
mesh for representing the standard lowest-degree continuous finite element expansion is first provided and
then additional degrees of freedom are introduced in order to enhance the approximation defined originally
over the background mesh. The same idea has also been adapted to create other alternate finite element
formulations. For example, it is possible to enhance the approximation on the background mesh by intro-
ducing additional degrees of freedom associated with either faces or edges of the elements and omitting then
the original nodal degrees of freedom. This leads to a suitable set of unknowns for creating discretizations
based on the face or edge element interpolation. If L2(Ω) is used to denote the set of square-integrable scalar
functions and Ω⊂Rd, we are thus led to bases for approximating vector fields in finite dimensional versions
of the spaces
H(div,Ω) = {v∈L2(Ω)d|div v∈L2(Ω)}
or
H(curl,Ω) = {v∈L2(Ω)d|curl v∈L2(Ω)d}.
A physical motivation for using these spaces is that fields with only either normal or tangential continuity
on boundaries can then be approximated in a natural manner. They may also be used to approximate various
mixed variational formulations. In this appendix we shall detail the construction of available bases for
defining face and edge interpolation elements. Currently the selection of edge interpolation elements is
more extensive.
F.1 The construction of face element interpolation
Here we focus on the construction of finite element spaces Xh(Ω) ⊂H(div,Ω). To this end, we use local
spaces X(K)which are defined over elements Kin the finite element mesh Mhassociated with the region
Ω. Then H(div,Ω)-conforming approximations are obtained as members of the space
Xh={vh∈H(div,Ω) |vh|K∈X(K)∀K∈ Mh}.
We note that the functions belonging to Xhmust have continuous normal components across the interior
element faces of the mesh Mh.
Any real element K⊂Ωis obtained as an image of a reference element ˆ
Kunder the element mapping
fKso that K=fK(ˆ
K). The element space X(K)is defined via introducing a local space X(ˆ
K)on the
reference element and showing how its members are transformed in connection with the element mapping
to get the basis for spanning X(K). To this end, we use the Piola transformation so that the connection
between the two element spaces is given by
X(K) = {vh|vh(x) = 1
det F(f−1
K(x)) F(f−1
K(x))(ˆ
vh◦f−1
K)(x),with ˆ
vh∈X(ˆ
K)}
where F=∇fKis the gradient of fKwith respect to ˆ
x∈ˆ
K. It then follows that the basis functions for
X(K)are obtained by applying the Piola transformation to the basis of X(ˆ
K).
CSC – IT Center for Science
F. Face and edge elements 145
Currently the Elmer solver contains the implementation of face element bases for triangular, quadrilateral
and tetrahedral elements. In the case of triangular and quadrilateral elements the reference element is chosen
to be either an equilateral triangle or a square, while tetrahedral elements use a regular tetrahedron as the
reference element. These choices are convenient as using reference elements which have the faces of an
equal size simplifies implementation. These reference elements have also been employed in connection with
the high-order elements. The element mapping fK:ˆ
K→R3is always defined in terms of the lowest-order
Lagrange interpolation basis functions ˆ
Ljassociated with the reference element as
x=fK(ˆ
x) = X
j
ˆ
Lj(ˆ
x)xj,
where ˆ
xis a point on the reference element and the points xjgive the vertices of the real element K.
F.1.1 The Raviart–Thomas space and its extension to three dimensions
In two dimensions (d= 2) the simplest face element is the triangular Raviart–Thomas element of order
k= 0. The reference element space is then chosen to be
X(ˆ
K) = RT 0(ˆ
K) = [P0(ˆ
K)]d⊕P0(ˆ
K)ˆx1
ˆx2(F.1)
where Pk(ˆ
K)is the space of polynomials of at most degree kon ˆ
K. The dimension of X(ˆ
K)is thus 3.
To define a basis for RT 0(ˆ
K)each element face ˆ
Fij connecting the reference element vertices ˆ
xiand
ˆ
xjis associated with one element basis function ˆ
ψij . We choose
ˆ
ψ12 =−1
20
1+1
2√3ˆx1
ˆx2,
ˆ
ψ23 =1
2√31
0+1
2√3ˆx1
ˆx2,
ˆ
ψ31 =−1
2√31
0+1
2√3ˆx1
ˆx2.
(F.2)
A finite element expansion
ˆ
vh(ˆ
x) = ˆ
d12 ˆ
ψ12(ˆ
x) + ˆ
d23 ˆ
ψ23(ˆ
x) + ˆ
d31 ˆ
ψ31(ˆ
x)
is referred to as the Raviart-Thomas interpolation function. Here the degrees of freedom ˆ
dij describe the
flux across the element faces as
ˆ
dij =Z
ˆ
Fij
ˆ
vh(ˆ
x)·ˆ
n(ˆ
x)dˆs,
where ˆ
ndenotes the outward unit normal associated with the element face.
The tetrahedral element of the Nédélec face element family of the first kind generalizes the use of the
space (F.1) to three dimensions (d= 3). If we introduce the Whitney form ˆ
φijk by
ˆ
φijk =ˆ
λi∇ˆ
λj×∇ˆ
λk−ˆ
λj∇ˆ
λi×∇ˆ
λk+ˆ
λk∇ˆ
λi×∇ˆ
λj,
the reference element basis functions ˆ
ψijk associated with the element faces ˆ
Fijk may be defined as
ˆ
ψ213 =−2ˆ
φ123,ˆ
ψ124 = 2ˆ
φ124,ˆ
ψ234 = 2ˆ
φ234,ˆ
ψ314 =−2ˆ
φ134,
so that
ˆ
ψ213 =
0
−√6/12
−1/√3
+√2
4
ˆx1
ˆx2
ˆx3
,ˆ
ψ124 =
0
−√6/4
0
+√2
4
ˆx1
ˆx2
ˆx3
,
ˆ
ψ234 =
√2/4
0
0
+√2
4
ˆx1
ˆx2
ˆx3
,ˆ
ψ314 =
−√2/4
0
0
+√2
4
ˆx1
ˆx2
ˆx3
.
(F.3)
CSC – IT Center for Science
F. Face and edge elements 146
The degrees of freedom for the corresponding finite element expansion ˆ
vh(ˆ
x)expressed in terms of these
basis functions are then characterized by
ˆ
dijk =Z
ˆ
Fijk
ˆ
vh(ˆ
x)·ˆ
n(ˆ
x)dˆ
S.
It can be verified that under the Piola transformation the basis functions of the reference element space
X(ˆ
K) = RT 0(ˆ
K)transform to a basis for X(K) = RT 0(K). In addition, the degrees of freedom
associated with the resulting basis functions ψij or ψijk are then characterized by
dij =Z
Fij
vh(x)·n(x)ds or dijk =Z
Fijk
vh(x)·n(x)dS.
In order to define a unique set of degrees of freedom for assembling global finite element equations, each
element face in the mesh Mhis associated with a uniquely defined unit normal m, so that either n(x)·
m(x) = 1 or n(x)·m(x) = −1for all xon the same element face. The degrees of freedom for assembling
the global equations from elementwise contributions are then selected to be
Dij =Z
Fij
vh(x)·m(x)ds or Dijk =Z
Fijk
vh(x)·m(x)dS.
The restriction of the global finite element expansion vhon the element can then be written in terms of the
local basis functions and the global degrees of freedom as
vh|K=Xdij ψij =XDij (n·m)x∈Fij ψij
or
vh|K=Xdijk ψijk =XDijk (n·m)x∈Fijk ψijk .
In practice this is implemented by applying sign reversions to the local basis functions when n·m=−1.
It should be noted that the normal continuity of the resulting face element expansion is then in-built.
F.1.2 The Brezzi-Douglas-Marini space and its extension to three dimensions
In the following we consider choosing the reference element space as X(ˆ
K) = [P1(ˆ
K)]d. In two dimensions
this choice leads to the construction of the Brezzi-Douglas-Marini space of order 1 (BDM1), while in three
dimensions it brings us to the definition of the Nédélec face elements of the second kind.
It seems that ways to obtain computational bases for face element discretizations are far less standardized
in comparison with implementing the basic element types. In Elmer we have followed the idea employed
in [1] to obtain a basis for the reference element space. Each basis function may then be associated with a
certain numerical integration point located on one of the element faces. The normal component of the basis
function is normalized to yield the unity at the integration point to which the basis function is associated,
while the normal component of this basis function is required to vanish at other integration points located
on the element faces. It follows from this strategy that convenient integral descriptions of the degrees of
freedom are also obtained.
For triangular elements the dimension of X(ˆ
K)is 6 and hence two basis functions ˆ
ψk
ij ,k= 1,2, are
associated with each face ˆ
Fij . A possible choice would be to use ˆ
λirot ˆ
λjand ˆ
λjrot ˆ
λias these basis
functions, but the basis functions implemented into Elmer are chosen to be their linear combinations such
that the normal components of the resulting alternate basis functions satisfy the version of the Kronecker
delta property explained above. We note that if the element face is parametrized as ˆ
Fij ={ˆ
x=ˆ
s(ξ)|−1≤
ξ≤1}, the two basis functions are associated with the Gaussian points ξ=−1/√3and ξ= 1/√3. The
degrees of freedom associated with the face element interpolation function ˆ
vh(ˆ
x)can then be characterized
by
ˆ
dk
ij =Z
ˆ
Fij
ˆ
vh(ˆ
x)·ˆ
n(ˆ
x)ˆ
φk(ˆ
x)dˆs,
CSC – IT Center for Science
F. Face and edge elements 147
where ˆ
φk(ˆ
x)are the linear Lagrange interpolation functions which have the Gaussian points of the element
face as the interpolation points.
The construction of the basis for the reference element space in three dimensions is analogous to the
two-dimensional case. The dimension of X(ˆ
K)is then 12 and each face ˆ
Fijk of the reference tetrahedron
is associated with three basis functions. The basis functions ˆ
ψl
ijk,l= 1,2,3, associated with each face ˆ
Fijk
are determined as linear combinations of the alternate basis functions (see, for example, [4]) ˆ
λk∇ˆ
λi×∇ˆ
λj,
ˆ
λj∇ˆ
λi×∇ˆ
λkand ˆ
λi∇ˆ
λj×∇ˆ
λk. Parametrizing the element face in the coordinates ξand ηas shown in
Figure E.5, the three basis functions for the element face are associated with the numerical integration points
(−1/2,√3/6),(1/2,√3/6), and (0,2/√3). The degrees of freedom for the face element interpolation
function ˆ
vh(ˆ
x)in terms of this basis may again be characterized as
ˆ
dl
ijk =Z
ˆ
Fijk
ˆ
vh(ˆ
x)·ˆ
n(ˆ
x)ˆ
φl(ˆ
x)dˆ
S,
where ˆ
φl(ˆ
x)are the linear Lagrange interpolation functions having the numerical integration points of the
element face as the interpolation points.
The Piola transformation is again used to transform the basis functions of the reference element space
X(ˆ
K) = [P1(ˆ
K)]dto a basis for X(K) = [P1(K)]d. Each element face in the mesh Mhis also associated
with a uniquely defined unit normal m, so that global finite element equations can be assembled using the
unique set of degrees of freedom as explained in connection with the Raviart-Thomas element. Additionally
an in-built mechanism is also applied to permute the order of basis functions so that the continuity of the
normal components is maintained when elementwise contributions are assembled to the global linear system
using the standard assembly subroutines of the Elmer solver.
F.1.3 The quadrilateral Arnold-Boffi-Falk element
The quadrilateral Arnold-Boffi-Falk edge element [3] is designed to provide the optimal accuracy even when
the physical elements are not obtained as affine images of the reference element. In this case the reference
element space is taken to be
X(ˆ
K) = ABF 0(ˆ
K) = Q2,0(ˆ
K)×Q0,2(ˆ
K)(F.4)
where Qm,n(ˆ
K)denotes the space of polynomial functions on ˆ
Kof degree at most min ˆx1and at most nin
ˆx2. The dimension of this space is 6. Four degrees of freedom are associated with the faces of the element
and the remaining two are local.
F.1.4 Using integration over the reference element
As with standard finite elements, it is convenient to convert integrals over elements Kto integrals over
the reference element Kin order to perform the numerical integration. Here the effect of using the Piola
transformation needs to be taken into account and hence the implementation of elementwise integration is
somewhat different as compared with using the standard transformation laws for the basic finite elements.
In the following u:K→Rdand v:K→Rdare the Piola transformations of ˆ
u:ˆ
K→Rdand
ˆ
v:ˆ
K→Rd. A scalar field φ:K→Ris obtained from a field ˆ
φ:ˆ
K→Rvia φ(x) = ( ˆ
φ◦f−1
K)(x).
In addition, grad and div are the gradient and divergence operators with respect to x, while ∇and Div are
used to denote the gradient and divergence operators with respect to ˆ
x. The following transformation laws
CSC – IT Center for Science
F. Face and edge elements 148
can then be used to convert integrals over elements Kto integrals over the reference element K:
Z
∂K
u(x)·n(x)φ(x)dS =Z
∂ˆ
K
ˆ
u(ˆ
x)·ˆ
n(ˆ
x)ˆ
φ(ˆ
x)dˆ
S,
Z
K
u(x)·v(x)dΩ = Z
ˆ
K
Fˆ
u(ˆ
x)·Fˆ
v(ˆ
x)1
det Fdˆ
Ω,
Z
K
divu(x)φ(x)dΩ = Z
ˆ
K
Div ˆ
u(ˆ
x)ˆ
φ(ˆ
x)dˆ
Ω,
Z
K
u(x)·grad φ(x)dΩ = Z
ˆ
K
ˆ
u(ˆ
x)·∇ˆ
φ(ˆ
x)dˆ
Ω.
(F.5)
F.1.5 Examples
The classic Poisson equation
−div grad u=fon Ω
can be written as a mixed problem such that the problem consists of finding simultaneously uin L2(Ω) and
the flux vector q=grad uin H(div,Ω). The files in the directories
•../trunk/fem/tests/RaviartThomas2D
•../trunk/fem/tests/BDM2D
•../trunk/fem/tests/RaviartThomas3D
•../trunk/fem/tests/BDM3D
contained in the source code repository of Elmer exemplify how the approximation of this mixed formulation
can be implemented using the face elements which have been considered in this appendix.
F.2 The construction of edge element interpolation
The definition of finite element spaces Xh(Ω) ⊂H(curl,Ω) is also given by using local spaces X(K)so
that H(curl,Ω)-conforming approximations are obtained as members of the space
Xh={vh∈H(curl,Ω) |vh|K∈X(K)∀K∈ Mh}.
We note that the functions belonging to Xhmust now have continuous tangential components across the
interior element faces of the finite element mesh Mhassociated with the region Ω.
The element space X(K)is again constructed via introducing a local space X(ˆ
K)on the reference
element and showing how its members are transformed in connection with the element mapping to obtain the
basis for spanning X(K). To this end, we use a transformation which is similar to the Piola transformation
used in the construction of the face element interpolation. Here the connection between the two element
spaces is given by
X(K) = {vh|vh(x) = [F(f−1
K(x))]−Tˆ
vh(f−1
K(x)),with ˆ
vh∈X(ˆ
K)}.
It then follows that the basis functions for X(K)are obtained by applying the H(curl)-version of the Piola
transformation to the basis of X(ˆ
K).
CSC – IT Center for Science
F. Face and edge elements 149
F.2.1 The lowest-order triangular and tetrahedral edge elements of Nédélec’s first
and second family
In the case of triangular and tetrahedral elements the Nédélec edge elements of the first kind are based on
choosing the reference element space as
X(ˆ
K) = [Pk−1(ˆ
K)]d+Sk(ˆ
K)(F.6)
where
Sk(ˆ
K) = {q∈[˜
Pk]d|q·x= 0 }(F.7)
with ˜
Pkthe set of homogeneous polynomials of degree k, i.e. the span of monomials of degree k. On the
other hand, the Nédélec edge elements of the second kind are obtained by selecting
X(ˆ
K) = [Pk(ˆ
K)]d.(F.8)
In the case of Nédélec’s first family of degree k= 1 each reference element edge ˆ
Eij connecting the
reference element vertices ˆ
xiand ˆ
xjis associated with one basis function ˆ
ψij . Given a finite element
expansion ˆ
vh(ˆ
x)in terms of these basis functions, the associated degrees of freedom ˆ
dij satisfy
ˆ
dij =Z
ˆ
Eij
ˆ
vh(ˆ
x)·ˆ
t(ˆ
x)dˆs,
where ˆ
tdenotes the unit tangent vector associated with the element edge.
The H(curl,Ω)-version of the Piola transformation is applied to transform the basis functions of the
reference element space X(ˆ
K)to a basis for X(K). The degrees of freedom associated with the resulting
basis functions ψij are then characterized by
dij =Z
Eij
vh(x)·t(x)ds,
with tthe unit tangent vector that is oriented similarly with ˆ
t. In order to define a unique set of degrees of
freedom for assembling global finite element equations, each element edge in the mesh Mhis associated
with a uniquely defined unit tangent τ, so that either t(x)·τ(x) = 1 or t(x)·τ(x) = −1for all xon
the same element edge. The degrees of freedom for assembling the global equations from elementwise
contributions are then selected to be
Dij =Z
Eij
vh(x)·τ(x)ds.
The restriction of the global finite element expansion vhon the element can then be written in terms of the
local basis functions and the global degrees of freedom as
vh|K=Xdij ψij =XDij (t·τ)x∈Eij ψij .
In practice this is implemented by applying sign reversions to the local basis functions when t·τ=−1. It
should be noted that the tangential continuity of the resulting edge element expansion is then in-built.
The strategy to construct basis functions for the Nédélec edge elements of the second kind is analogous to
that described in Section F.1.2 for the face elements. That is, each basis function is associated with a certain
numerical integration point located on one of the element edges. The tangential component of the basis
function is normalized to yield the unity at the integration point to which the basis function is associated,
while the tangential component of this basis function is required to vanish at other integration points located
on the element edges. Considering the lowest-degree case, in both two and three dimensions each element
edge is associated with two basis functions (as the dimension of X(ˆ
K)is 6 for the triangular element and 12
CSC – IT Center for Science
F. Face and edge elements 150
for the tetrahedral element). The degrees of freedom associated with the edge element interpolation function
ˆ
vh(ˆ
x)can then be characterized by
ˆ
dk
ij =Z
ˆ
Eij
ˆ
vh(ˆ
x)·ˆ
t(ˆ
x)ˆ
φk(ˆ
x)dˆs,
with ˆ
φk(ˆ
x)the linear Lagrange interpolation functions which have the Gaussian points of the element edge
as the interpolation points.
The basis of the element space X(K)is again obtained via applying the adopted version of the Piola
transformation. The global degrees of freedom are also defined in the same way as in the case of Nédélec’s
first family. In addition, an in-built mechanism is now applied to permute the order of basis functions so that
the continuity of the tangential components is maintained when elementwise contributions are assembled to
the global linear system using the standard assembly subroutines of the Elmer solver.
F.2.2 The optimal edge element family of the lowest degree
The simplest edge elements which have other shapes than a triangle or a tetrahedron may not provide op-
timal accuracy when the physical elements are not obtained as affine images of the reference element. In
this section, we describe a family of edge elements designed to provide optimal accuracy even for non-affine
element shapes. The triangular and tetrahedral elements of this family correspond to the Nédélec edge ele-
ment family of the first kind; cf. the previous section. The other edge element types including quadrilateral,
pyramidic, prismatic and hexahedral elements are described in the following.
The curl-conforming version of the quadrilateral Arnold-Boffi-Falk element
In this case the reference element space is taken to be
X(ˆ
K) = Q0,2(ˆ
K)×Q2,0(ˆ
K).(F.9)
The dimension of this space is 6. Four degrees of freedom are associated with the edges of the element and
the remaining two are local.
A 10-DOF pyramidic element
The construction of the pyramidic element is based on [2]. Using the reference element of the high-order
family, the reference element space is given by
X(ˆ
K) = [P0(ˆ
K)]3⊕span{ˆ
ϕj|j= 4,...,10 }(F.10)
where
ˆ
ϕ4=1
(1 −ˆx3/a)2
ˆx2(1 −ˆx3/a)
ˆx1(1 −ˆx3/a)
ˆx1ˆx2/a
,ˆ
ϕ5=
1−ˆx3/a
0
ˆx1/a
,ˆ
ϕ6=
0
1−ˆx3/a
ˆx2/a
,
ˆ
ϕ7=ˆx2
1−ˆx3/a
1−ˆx3/a
0
ˆx1/a
,ˆ
ϕ8=ˆx1
1−ˆx3/a
0
1−ˆx3/a
ˆx2/a
,
ˆ
ϕ9=ˆx2
2
(1 −ˆx3/a)2
1−ˆx3/a
0
ˆx1/a
,ˆ
ϕ10 =ˆx2
1
(1 −ˆx3/a)2
0
1−ˆx3/a
ˆx2/a
,
(F.11)
with a=√2.
CSC – IT Center for Science
F. Face and edge elements 151
A 15-DOF prismatic element
The construction of the prismatic element is also given in [2]. If the prismatic reference element is repre-
sented as ˆ
K=ˆ
C× { ˆx3| − 1≤ˆx3≤1}, with ˆ
Ca triangle on a plane ˆx3=const, the reference element
space is defined by
X(ˆ
K) = [N d1(ˆ
C)⊗P2({ˆx3| − 1≤ˆx3≤1})] ×P2(ˆ
C)(F.12)
where N d1(ˆ
C)denotes the Nédélec edge element space of the first kind and P2(D)is the space of the
second-order polynomials defined over the domain D.
A 27-DOF hexahedral element
The optimal hexahedral edge element has been developed in [5]. In this case the reference element space is
taken to be
X(ˆ
K) = Q0,2,2(ˆ
K)×Q2,0,2(ˆ
K)×Q2,2,0(ˆ
K).(F.13)
F.2.3 Second-order edge elements
Second-order edge finite elements are also available. They are members of the conventional Nédélec’s first
family. The background mesh can then consist of either the lowest-order or second-order nodal elements.
For example, a tetrahedral backgound mesh may consist of elements of type 504 or 510.
In the case of triangles and tetrahedra the reference element space is given by (F.6) with k= 2 and thus
has the dimension 8 or 20, respectively. In the case of quadrilateral elements, the reference element space is
X(ˆ
K) = Q1,2(ˆ
K)×Q2,1(ˆ
K),(F.14)
having the dimension 12.
If the prismatic reference element is represented as ˆ
K=ˆ
C׈
I, with ˆ
Ca triangle on a plane ˆx3=const
and ˆ
I={ˆx3| − 1≤ˆx3≤1}, the reference element space for the second-order prism is defined by
X(ˆ
K) = [N d2(ˆ
C)⊗P2(ˆ
I)] ×[P2(ˆ
C)⊗P1(ˆ
I)] (F.15)
where N d2(ˆ
C)denotes the second-order Nédélec edge element space of the first kind and Pk(ˆ
I)is the
space of the kth order polynomials defined over the domain ˆ
I. This space has the dimension 36. On the
other hand, the reference element space for bricks is
X(ˆ
K) = Q1,2,2(ˆ
K)×Q2,1,2(ˆ
K)×Q2,2,1(ˆ
K),(F.16)
having the dimension 54. Finally, the construction of the pyramidic element of degree 2 follows [2] by
omitting some of the basis functions associated with the quadrilateral face to obtain the conformity with the
elements of the Nédélec’s first family. The resulting element has 31 DOFs.
F.2.4 Using integration over the reference element
The effect of the applied version of the Piola transformation needs to be taken into account when integrals
over elements Kare converted to integrals over the reference element K. Therefore the implementation
of elementwise integration is again somewhat different as compared with using the standard transformation
laws for the basic finite elements.
In the following u:K→Rdand v:K→Rdare the transformations of ˆ
u:ˆ
K→Rdand ˆ
v:ˆ
K→Rd
under the H(curl)-version of the Piola transformation. A scalar field φ:K→Ris obtained from a field
ˆ
φ:ˆ
K→Rvia φ(x) = ( ˆ
φ◦f−1
K)(x). In addition, the curl operator with respect to xis denoted by curl,
while Curl is used to denote the curl operator with respect to ˆ
x. The following transformation laws can
then be used to convert integrals over elements Kto integrals over the reference element K:
Z
∂K
u(x)·t(x)φ(x)dS =Z
∂ˆ
K
ˆ
u(ˆ
x)·ˆ
t(ˆ
x)ˆ
φ(ˆ
x)dˆ
S,
Z
K
u(x)·v(x)dΩ = Z
ˆ
K
[F(ˆ
x)]−Tˆ
u(ˆ
x)·[F(ˆ
x)]−Tˆ
v(ˆ
x) det F(ˆ
x)dˆ
Ω.
(F.17)
CSC – IT Center for Science
BIBLIOGRAPHY 152
In addition, the referential description of the field curl u(x)is given by
(curl u)(fK(ˆ
x)) = 1
det F(ˆ
x)F(ˆ
x)Curl ˆ
u(ˆ
x).(F.18)
In two dimensions the only nontrivial component of the curl field transforms in particular as
(curl u·e3)(fK(ˆ
x)) = 1
det F(ˆ
x)Curl ˆ
u(ˆ
x)·ˆ
e3.(F.19)
In view of (F.18), we thus have
Z
K
curl u(x)·curl v(x)dΩ = Z
ˆ
K
F(ˆ
x)Curl ˆ
u(ˆ
x)·F(ˆ
x)Curl ˆ
v(ˆ
x)1
detF(ˆ
x)dˆ
Ω.(F.20)
F.2.5 Setting Dirichlet boundary conditions
In connection with edge element interpolation Dirichlet constraints can be specified to prescribe the tangen-
tial components of the unknown field on the boundary. In the case of the simplest edge elements the values
of the prescribed degrees of freedom can be generated by computing integrals of the type
D=Z
Eij
uh(x)·τ(x)φ(x)ds
over element edges Eij . In the case of elements of the second degree or elements from the optimal family
some degrees of freedom may be associated with element faces and a more general approach must be fol-
lowed. In that case the values of the degrees of freedom associated with the element faces are obtained by
seeking the best approximation of the data in L2.
Specific syntax must be used in the solver input file to define boundary conditions for the edge element
approximation. The following commands can be used in the Boundary Condition section to define
the boundary values of the solver variable Uwhich corresponds to a continuous field uand is approximated
with the edge element basis functions:
U {e} j Real
This keyword is used to define the vector gso that its tangential trace g×nis approximated by uh×n,
with uhthe finite element interpolating function. The value of this keyword defines the components
gj,j∈ {1,2,3}, with respect to the global Cartesian coordinate system.
U {e} Real
If the edge element has only degrees of freedom which are associated with edges, this keyword can be
used to define a scalar field ˆgτso that the values of the prescribed degrees of freedom are computed as
D=Z
Eij
[g(x)·τ(x) + ˆgτ]φ(x)ds.
This keyword should not be used if additional degrees of freedom associated with the element faces
are present (that is, this feature should not be used for elements of the second degree or elements from
the optimal family).
Bibliography
[1] Ervin V. J. Computational bases for RTkand BDMkon triangles. Computers & Mathematics with
Applications, 64:2765–2774, 2012.
[2] Bergot M. and Duruflé M. High-order optimal edge elements for pyramids, prisms and hexahedra.
Journal of Computational Physics, 232:189–213, 2013.
CSC – IT Center for Science
BIBLIOGRAPHY 153
[3] Arnold D. N., Boffi D., and Falk R. S. Quadrilateral H(div) finite elements. SIAM Journal on Numerical
Analysis, 42:2429–2451, 2005.
[4] Arnold D. N., Falk R. S., and Winther R. Geometric decompositions and local bases for spaces of finite
element differential forms. Computer Methods in Applied Mechanics and Engineering, 198:1660–1672,
2009.
[5] Falk R. S., Gatto P., and Monk P. Hexahedral H(div) and H(curl) finite elements. ESAIM: Mathematical
Modelling and Numerical Analysis, 45:115–143, 2011.
CSC – IT Center for Science
Appendix G
Advanced finite element definitions
Advanced finite element formulations may require that several fields are approximated simultaneously using
different types of interpolation. For example curl-conforming (edge finite element) approximation may need
to be used in combination with the standard Lagrange interpolation of another field. This requires advanced
flexibility in creating suitable degrees of freedom for finite elements.
Elmer offers general support for creating advanced finite element formulations when a background mesh
for representing the standard lowest-degree continuous finite element expansion is available. Additional
degrees of freedom may then be added by associating them with the edges, faces or interiors of the elements
of the background mesh. This enhancement strategy is used internally for example in the creation of the
second-order curl-conforming finite elements. In general, the writer of a solver must nevertheless decide
how the additional degrees of freedom are utilized in the finite element approximation.
G.1 Specification of the degrees of freedom
Each solver specification in the sif file may contain an additional command of the type
Element = "n:N... "
where the string between the quotation marks controls what degrees of freedom (DOFs) are allocated. The
specification string can contain the following parts:
•"n:N"to create N≥0DOFs associated with each node of the background mesh.
•"e:E"to create E≥0DOFs associated with each edge of the background mesh.
•"f:F"to create F≥0DOFs associated with each face contained in the background mesh.
•As a more general alternate for the previous definition, "-quad_face b:FQ-tri_face b:FT"
associates FQDOFs with every quadrilateral face and FTDOFs with every triangular face of the back-
ground mesh.
•"b:B"to create B≥0DOFs associated with each element interior
•As a more general alternate for the previous definition, use
"-tetra b:BT-pyramid b:Bp-prism b:BW-brick b:BB"
controls the number of the local DOFs associated with the element interiors in a manner that depends
on the element shape (tetrahedron, pyramid, prism or brick).
CSC – IT Center for Science
Index
-dofs, 58
-global, 58
-nooutput, 58
Adaptive Error Limit, 66
Adaptive Error Measure, 37
Adaptive Keep Smallest, 37
Adaptive Max Change, 66
Adaptive Max H, 66
Adaptive Mesh Refinement, 66
Adaptive Min H, 66
Adaptive Min Timestep, 37
Adaptive Remesh, 66
Adaptive Save Mesh, 66
Adaptive Time Error, 37
Adaptive Timestepping, 36
After Linsolve, 30
algebraic multigrid, 21
Anti Periodic BC VarName, 47
Anti Rotational Projector, 49
Anti Sliding Projector, 49
Apply Contact BCs, 50,53
Apply Limiter, 45
Apply Mortar BCs, 49,53
Backward Differences Formulae, 34
bandwidth optimization, 26
BDF, 34
BDF Order, 35
Before Linsolve, 29,45
BiCGstab(ℓ), 24
BiCGstab, biconjugate gradient stabilized, 24
BiCGstabl polynomial degree, 25
Binary output, 118
Biorthogonal Basis, 50
Biorthogonal Dual LCoeff, 50
Biorthogonal Dual Master, 50
Body, 8
Body Force, 8,44,45,55,59
BoomerAMG, 76
BoomerAMG Coarsen Type, 76
BoomerAMG Cycle Type, 77
BoomerAMG Interpolation Type, 77
Boomeramg Max Levels, 77
BoomerAMG Num Functions, 77
BoomerAMG Num Sweeps, 77
BoomerAMG Relax Type, 76
BoomerAMG Smooth Type, 77
Bossak Alpha, 36
Boundary Condition, 9,44–46,49,50,54
Boundary Condition bc_id, 63
Calculate Energy Norm, 56
Calculate Loads, 55
Calculate Weights, 56
CG, conjugate gradient, 24
CGS, conjugate gradient squared, 24
Check Keywords, 12
compilation instructions, 80
Component, 9
component name, 57
Compressed Row Storage, 28
Constants, 8
Constrained Modes Varname, 54
constraint modes, 54
Constraint Modes Analysis, 54
Constraint Modes Fluxes, 54
Constraint Modes Fluxes Symmetric, 54
Constraint Modes Lumped, 54
Contact Active Condition, 50
Contact Active Set Minimum, 51
Contact Depth Offset, 50
Contact Depth Offset Initial, 50
Contact Passive Condition, 51
Contact Type, 51
Contact Velocity(dim), 51
Convergence Tolerance, 32
Coordinate Mapping(dim), 60
Coordinate Scaling, 60
Coordinate Transformation , 60
Coordinate Transformation Radius, 60
Corrector method, 38
Crank-Nicolson method, 34
Cuthill-McKee algorithm, 26
Cylindrical Projector, 49
Direct methods, 19
direct solver, 24
Dirichlet constraints, 43
Discontinuous BC, 63
Discontinuous Boundary, 63
Discontinuous Boundary Greedy, 63
CSC – IT Center for Science
INDEX 156
Discontinuous Bulk Greedy, 63
Discontinuous Galerkin method, 73
Discontinuous Target Bodies, 63
domain decomposition, 69,74
Dynamic Friction Coefficient, 51
Edge Basis, 49
Eigen Analysis, 41
Eigen System Complex, 41
Eigen System Compute Residuals, 41
Eigen System Convergence Tolerance, 41
Eigen System Damped, 41
Eigen System Lanczos Vectors, 41
Eigen System Max Iterations, 41
Eigen System Select, 41
Eigen System Shift, 41
Eigen System Use Identity, 41
Eigen System Values, 41
Eigenvalue problems, 39
element definition command, 128
Eliminate From Master, 53
Eliminate Linear Constraints, 53
Eliminate Slave, 53
ElmerPost input file, 122
ElmerSolver mesh file, 115
Equation, 8
error indicator, 64
Exec Condition, 57
Exec Intervals, 57
Exec Solver, 57
Export Lagrange Multiplier, 53
Exported Variable i, 58
Exported variables, 58
Extruded Coordinate Index, 62
Extruded Max Coordinate, 62
Extruded Mesh Density, 62
Extruded Mesh Levels, 62
Extruded Mesh Ratio, 62
Extruded Min Coordinate, 62
Friction Contact, 51
GCR, generalized conjugate residual, 24
geometric multigrid, 21
Geometric Stiffness, 41
GMRES, generalized minimal residual, 24
Header, 8
ILU preconditioners, 21
ILU, incomplete LU-decomposition, 25
Initial Condition, 9
Initialization of dependent variables, 17
Initialize Dirichlet Condition, 17
iterative solver, 24
keyword syntax, 12
Krylov methods, 24
Krylov subspace methods, 20
Lagrange Multiplier Name, 53
LAPACK, 19
Level Projector, 49
Level Projector Extruded Edges Strong, 50
Level Projector Nodes Strong, 49
Level Projector Plane Edges Strong, 50
limiter, 45
Limiter Load Sign Negative, 45
Limiter Load Tolerance, 45
Limiter Value Tolerance, 45
linear algebra, 19
Linear System Abort Not Converged, 25
Linear System Convergence Tolerance, 25
Linear System Direct Method, 24
Linear System GCR Restart, 25
Linear System GMRES Restart, 25
Linear System ILUT Tolerance, 25
Linear System Iterative Method, 24
Linear System Max Iterations, 25
Linear System Precondition Recompute, 26
Linear System Preconditioning, 25,75,76
Linear System Residual Mode, 25
Linear System Residual Output, 26
Linear System Solver, 24
Linear System Symmetric, 45
Linear System Timing, 59
Linear System Timing Cumulative, 59
Linear System Use Hypre, 75
linear systems, 19
lower limit, 45
Lumped Mass Matrix, 36
MATC, 12,13
parameterized keyword commands, 14
Material, 8
mesh database, 70
Mesh DB, 9
Mesh Grading Power, 61
Mesh Keep, 61
Mesh Keep Grading, 61
Mesh Levels, 61
mesh refinement, 64
MG Cluster Alpha, 28
MG Cluster Size, 28
MG Convergence Tolerance, 26
MG Direct Interpolate, 28
MG Direct Interpolate Limit, 28
MG Eliminate Dirichlet, 27
MG Eliminate Dirichlet Limit, 27
MG Equal Split, 26
MG ILUT Tolearance, 27
CSC – IT Center for Science
INDEX 157
MG Levels, 26
MG Lowest Linear Solver Limit, 27
MG Max Iterations, 26
MG Mesh Name, 26
MG Positive Connection Limit, 28
MG Post Smoothing Iterations, 27
MG Pre Smoothing Iterations, 27
MG Preconditioning, 27
MG Projection Limit, 28
MG Recompute Projector, 27
MG Smoother, 27
MG SOR Relax, 27
MG Strong Connection Limit, 28
MG Strong Connection Minimum, 28
Mortar BC, 49
Mortar BC Scaling, 49
multigrid solver, 24
Mumps, 20
nodal loads, 55
Nonlinear System Convergence Absolute, 32
Nonlinear System Convergence Measure, 31
Nonlinear System Convergence Tolerance, 32
Nonlinear System Convergence Without Constraints,
53
Nonlinear System Max Iterations, 32
Nonlinear System Newton After Iterations, 32
Nonlinear System Newton After Tolerance, 32
Nonlinear System Norm Degree, 32
Nonlinear System Norm Dofs, 32
Nonlinear System Relaxation Factor, 32
Nonlinear Update Exported Variables, 58
nonlinearity, 31
Normal-Tangential VarName, 51
Optimize Bandwidth, 26
Output Intervals, 36
Parasails, 75
Pardiso, 20
Periodic BC, 46
Periodic BC Explicit, 46
Periodic BC Rotate(3), 46
Periodic BC Scale(3), 47
Periodic BC Translate(3), 46
Periodic BC Use Lagrange Coefficients, 47
Periodic BC VarName, 47
Periodic conditions, 46
Plane Projector, 49
preconditioning, 21
Predictor method, 38
Predictor-Corrector Control, 38
Predictor-Corrector Control Beta 1, 38
Predictor-Corrector Control Beta 2, 38
Predictor-Corrector Control Tolerance, 38
Predictor-Corrector Save Error, 38
Predictor-Corrector Scheme Order, 38
Preserve Baseline, 62
Projector Skip Edges, 50
Projector Skip Nodes, 50
Radial Projector, 49
Reference Norm, 59
Reference Norm Tolerance, 59
Relaxation Factor, 32
Restart Before Initial Conditions, 17
Restart file, 16
Restart Position, 16
Rotational Projector, 49
Save Projector, 50
Set Dirichlet BCs by BC Numbering, 44
sif file, 8
Simulation, 8,43,60–63
Skip First Timestep, 38
Slide Contact, 51
Sliding Projector, 49
Solver, 9,45,49,50,54–56,58,59
Solver activation, 57
solver input file, 8
section names, 8
Solver Timing, 59
Solver Timing Cumulative, 59
SOLVER.KEYWORDS, 12
Stability Analysis, 41
Start Time, 57
Static Friction Coefficient, 51
Steady State Convergence Measure, 32
Steady State Convergence Tolerance, 33
Steady State Max Iterations, 33
Steady State Min Iterations, 33
Steady State Norm Degree, 32
Steady State Norm Dofs, 32
Steady State Relaxation Factor, 33
Stick Active Condition, 51
Stick Contact, 51
Stick Contact Coefficient, 51
Stick Passive Condition, 51
Stop Time, 57
SuperLU, 20
Target Boundaries(n), 44
Target Coordinates(n,DIM), 44
Target Nodes(n), 44
Tie Contact, 51
Time Derivative Order, 35
time-dependent, 34
timeintegration, 34
Timestep Intervals, 36
Timestep Size, 36
CSC – IT Center for Science
INDEX 158
Timestep Sizes, 36
Timestepping Method, 35
U {e}, 150
U {e} j, 150
Umfpack, 19
Update Exported Variables, 58
upper limit, 45
Use Transpose values, 53
Variable name, 57
Varname, 44
Varname {e}, 44
Varname {e} i, 44
Varname Condition, 44,45
Varname i, 44
Varname Load , 55
Varname Lower Limit, 45
Varname Passive, 59
Varname Upper Limit, 45
CSC – IT Center for Science

Navigation menu