Gascoigne3D Manual
Gascoigne3D_manual
Gascoigne3D_manual
User Manual:
Open the PDF directly: View PDF 
.
Page Count: 42
| Download | |
| Open PDF In Browser | View PDF |
Introduction to
Gascoigne 3D
High Performance Adaptive Finite Element Toolkit
Tutorial and manuscript by Thomas Richter
thomas.richter@iwr.uni-heidelberg.de
March 9, 2014
2
Contents
1 Introduction
1.1 Installation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 A minimal example solving a partial differential equation . . . . . . . . . . .
5
5
7
2 The parameter file
13
3 Description of the problem
3.1 The right hand side . . . . . . . . . . . . . . . . . . . . . .
3.2 Boundary Data . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Dirichlet boundary data . . . . . . . . . . . . . . . .
3.2.2 Neumann conditions . . . . . . . . . . . . . . . . . .
3.2.3 Robin conditions . . . . . . . . . . . . . . . . . . . .
3.3 Definition of the partial differential equations . . . . . . . .
3.3.1 Nonlinear equations . . . . . . . . . . . . . . . . . .
3.3.2 Equations on the boundary & Robin boundary data
3.4 Exact Solution and Evaluation of Errors . . . . . . . . . . .
3.5 Systems of partial differential equations . . . . . . . . . . .
3.5.1 Modification of the right hand side . . . . . . . . . .
3.5.2 Modification of Dirichlet data . . . . . . . . . . . . .
3.5.3 Modification for the equation, form and matrix . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
15
15
16
17
18
19
19
21
22
23
24
24
25
25
4 Time dependent problems
4.1 Time discretization of the heat equation . . . . . . . . . . . . . . . . . . . . .
4.2 Time dependent problem data . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Non standard time-discretizations . . . . . . . . . . . . . . . . . . . . . . . . .
27
27
30
30
5 Mesh handling
5.1 Definition of coarse meshes . . . . .
5.2 Curved boundaries . . . . . . . . . .
5.3 Adaptive refinement of meshes . . .
5.3.1 Estimating the Energy-Error
5.3.2 Picking Nodes for Refinement
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
33
33
35
37
37
38
6 Flow problems and stabilization
6.1 The function point . . . . . . . .
6.2 Stabilization of the Stokes system
6.3 Stabilization of convective terms
6.4 LPS Stabilization . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
39
39
39
40
41
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
3
Contents
4
1 Introduction
When referring to files, the following directories are used:
HOME Your home directory
GAS The base-directory of Gascoigne. For the workshop:
GAS/src The source- and include-files of Gascoigne.
Gascoigne is a software-library written in C++. To solve a partial differential equation with
Gascoigne, a small program has to be written using the routines provided by Gascoigne. All
functionality is gathered in the library-file File: libGascoigneStd . This library is linked
to every application.
1.1 Installation
In this section we explain the usage of Gascoigne as a library for own applications.
Gascoigne requires a configuration file:
File: HOME/.gascoignerc
Here, some setting like “where is Gascoigne”, “where are the libraries” are given. The
minimal file looks like:
1
2
#############################################
SET(GASCOIGNE STD ” ˜/ S o f t w a r e / s r c / G a s c o i g ne ” )
3
4
5
6
7
#############################################
# Gascoig n e L i b r a r i e s & B i n a r i e s
SET(GASCOIGNE LIBRARY OUTPUT PATH ” ˜/ S o f t w a r e / x 8 6 6 4 / l i b ” )
SET(GASCOIGNE EXECUTABLE OUTPUT PATH ” . ” )
You need a file like this in the home directory. To build an application using Gascoigne at
least two files and directories are required:
File: TEST/bin/
File: TEST/src/main.cc
File: TEST/src/CMakeLists.txt
5
1 Introduction
The first directory will be used to store the executable file for this application. The other
two files are the source code of the application and a file to control the compilation of the
application. This could look like:
File: TEST/src/CMakeLists.txt
1
INCLUDE($ENV{HOME} / . g a s c o i g n e r c )
2
3
4
INCLUDE( ${GASCOIGNE STD}/ CMakeGlobal . t x t )
LINK LIBRARIES ( $ {GASCOIGNE LIBS} )
5
6
ADD EXECUTABLE( ” Test ”
main . c c )
In the first line, the user’s config file is read. The following two lines set global parameters for Gascoigne. The variable GASCOIGNE STD is provided in the config-file, the variable
GASCOIGNE LIBS is given in the file CMakeGlobal.txt. The last line indicates the name of
the application: Test, and lists all code-files to be compiled into the application.
A minimal example for a code-file (it does not really calculate something) looks like this:
File: TEST/src/main.cc
1
2
#include
#include ” s t d l o o p . h”
3
4
5
6
7
8
int main ( int argc , char ∗∗ argv )
{
s t d : : c o u t << ” C r e a t e a G a s c oi g n e Loop ! ” << s t d : : e n d l ;
Gascoi g n e : : StdLoop Loop ;
}
For building the application we use the program cmake to generate Makefiles. In the directory
TEST/bin/ type:
1
2
cmake . . / s r c
make
The first line creates the File: TEST/bin/Makefile . It reads the File: CMakeLists.txt
to put everything together. The second command starts the compilation and linking of the
application and if everything works out, File: TEST/bin/Test should be the executable
file which then can be started.
The command cmake is only needed to create new Makefiles. This is necessary, if new
code-files are added to the application. If only source code is changed, e.g. the file File:
TEST/src/main.cc , you only need to call make in the directory TEST/bin.
6
1.2 A minimal example solving a partial differential equation
1.2 A minimal example solving a partial differential equation
Assume there exists a directory FIRSTEXAMPLE with sub-dirs FIRSTEXAMPLE/bin and
FIRSTEXAMPLE/src. In the src-directory you need the files File: main.cc and File:
CMakeLists.txt . in the bin-directory you call cmake ../src for configuration and then
make for compilation.
File: FIRSTEXAMPLE/src/main.cc
1
2
3
4
5
#include
#include
#include
#include
#include
” s t d l o o p . h”
” p r o b l e m d e s c r i p t o r b a s e . h”
” c o n s t a n t r i g h t h a n d s i d e . h”
” z e r o d i r i c h l e t d a t a . h”
” l a p l a c e 2 d . h”
6
7
using namespace G a s c oi g n e ;
8
9
10
11
12
13
14
15
16
17
c l a s s Problem : public P r o b l e m D e s c r i p t o r B a s e
{
public :
s t d : : s t r i n g GetName ( ) const { return ” F i r s t Problem ” ; }
void B a s i c I n i t ( const ParamFile ∗ p a r a m f i l e )
{
G e t E q u a t i o n P o i n t e r ( ) = new L ap l a ce 2 d ;
GetRightHandSidePointer ( ) = new OneRightHandSideData ( 1 ) ;
G e t D i r i c h l e t D a t a P o i n t e r ( ) = new Z e r o D i r i c h l e t D a t a ;
18
ProblemDescriptorBase : : B a s i c I n i t ( paramfile ) ;
19
}
20
21
};
22
23
24
25
26
int main ( int argc , char ∗∗ argv )
{
ParamFile p a r a m f i l e ( ” f i r s t . param” ) ;
i f ( argc >1) p a r a m f i l e . SetName ( argv [ 1 ] ) ;
27
Problem PD;
PD. B a s i c I n i t (& p a r a m f i l e ) ;
28
29
30
ProblemContainer PC;
PC. AddProblem ( ” l a p l a c e ” , &PD) ;
31
32
33
StdLoop l o o p ;
l o o p . B a s i c I n i t (& p a r a m f i l e , &PC ) ;
l o o p . run ( ” l a p l a c e ” ) ;
34
35
36
37
}
7
1 Introduction
The Problem
Here, the problem to be solved is defined. The Problem gathers all information on the
problem to be solved. This is: which pde, what kind of boundary data, what is the right
hand side, and so on.
By PD.BasicInit() we initialize the problem (will be described later on). This problem is
then added to another class, the ProblemContainer. Here, we can gather different problems,
each together with a keyword, which is laplace in this case. There are applications where
it is necessary to solve different problems in the same program.
The problem set pointers to a large variety of classes. Here, we define the equation to be of
type Laplace2d, an implementation of the two dimensional Laplace equation. For right hand
side and Dirichlet data we define self written classes (will be explained later). In Gascoigne
a very large set of parameters can be set for a problem, all in form of pointers. However these
three classes, Equation, RightHandSide and DirichletData are necessary for every instance
of Gascoigne. The class Problem is derived from the base-class ProblemDescriptorBase:
File: GAS/src/Problem/problemdescriptorbase.h
The equation Laplace2d is part of the Gascoigne-Library:
File: GAS/src/Problem/laplace2d.h
This is explained in detail in Section 3.
The Loop
Then, a Loop is specified, the loop is the class controlling what is happening at run-time.
and what Gascoigne is doing: this is usually a sequence like: initialize, solve a pde, write
the solution to the disk, evaluate functionals, compute errors, refine the mesh and start over.
In BasicInit() we initialize the most basic structure and read data from parameter files.
Here we also pass the ProblemContainer to the different instances.
Finally, the equation is solved by calling loop.run("laplace"), where the keyword laplace
indicates what problem to solve. Let us now have a closer look at this function which is
defined in
File: FIRSTEXAMPLE/src/loop.cc
First, in run(), we define two vectors u,f for storing the solution and the right hand side.
Then, in the for-loop we initialize the problem, solve the problem, write out the solution and
refine the mesh. This is very typical for every program run, since for a reliable simulation it
is necessary to observe the convergence of the solution on a sequence of meshes. Next, we
have a closer look at every step of the inner loop:
8
1.2 A minimal example solving a partial differential equation
G e t M u l t i L e v e l S o l v e r ()−> R e I n i t ( p r o b l e m l a b e l ) ;
G e t M u l t i L e v e l S o l v e r ()−> R e I n i t V e c t o r ( u ) ;
G e t M u l t i L e v e l S o l v e r ()−> R e I n i t V e c t o r ( f ) ;
1
2
3
4
G e t S o l v e r I n f o s ()−>GetNLInfo ( ) . c o n t r o l ( ) . m a t r i x m u s t b e b u i l d ( ) = 1 ;
5
The first line initializes the basic ingredients of Gascoigne. In particular, we pass the equation to be solved and all the necessary data to different parts like “Solver”, “Discretization”.
Then, by GetMultiLevelSolver()->ReInitVector(u) we initialize the vector on the current mesh. This basically means that we resize the vector to match the number of degrees
of freedom (which changes with the mesh-size).
Finally, in the last line we tell Gascoigne to assemble the system matrix. This is necessary
whenever the system matrix changes, which is the case if the mesh changes, or if the matrix
depends on the solution u itself.
File: Template1/src/loop.cc
S o l v e ( u , f , ” R e s u l t s /u” ) ;
GetMeshAgent()−> g l o b a l r e f i n e ( ) ;
1
2
By calling Solve(u,f) we solve the PDE. The subsequent line refines the mesh.
File: Template1/src/loop.cc
1
2
3
4
5
s t r i n g Loop : : S o l v e ( V e c t o r I n t e r f a c e& u , V e c t o r I n t e r f a c e& f , s t r i n g
name )
{
G e t M u l t i L e v e l S o l v e r ()−> G e t S o l v e r ()−>Zero ( f ) ;
G e t M u l t i L e v e l S o l v e r ()−> G e t S o l v e r ()−>Rhs ( f ) ;
6
G e t M u l t i L e v e l S o l v e r ()−> G e t S o l v e r ()−>SetBoundaryVector ( f ) ;
G e t M u l t i L e v e l S o l v e r ()−> G e t S o l v e r ()−>SetBoundaryVector ( u ) ;
7
8
9
s t r i n g s t a t u s = G e t M u l t i L e v e l S o l v e r ()−>
S o l v e ( u , f , G e t S o l v e r I n f o s ()−>GetNLInfo ( ) ) ;
10
11
12
Output ( u , name ) ;
13
14
}
Here, we first assemble the right hand side f of our problem. Then, we need to initialize
Dirichlet boundary values. Finally, we tell the multigrid solver to solve the equation.
In the last line, we write the solution into a file for visualization.
9
1 Introduction
Output of Gascoigne
After starting the code, Gascoigne first prints out some information on the current mesh,
solver, Discretization and problem data.
Then, for every iteration of the inner loop the output looks like:
...
================== 2 ================ [l,n,c] 5 1089 1024
0: 9.77e-04
M 1: 7.73e-09 [0.00 0.00]
- 6.26e-08 [0.012] {3}
[u.00002.vtk]
...
In the first line the current iteration count (2) , the number of mesh-refinemenet levels (5),
the number of mesh nodes (1089) and the number of mesh elements (quads) (1024) is given.
Then, Gascoigne prints control and stastic data for the solution of the algebraic systems.
Gascoigne solves every equation with a Newton-method, even linear ones. This is for reasons
of simplicity as well as having a defect correction method for a better treatment of roundoff
errors. On the left side of the output the convergence history of the newton method is
printed. Here we only need one Newton step. Before this step, the residual of the equation
was 9.77e − 4, after this one step the residual was reduced to 7.73e − 9. The letter M indicates
that Gascoigne assembled a new system matrix in this step. The folling two numbers [0.00
0.00] indicate the convergance rates of the newton iteration by printing the factor, by which
the residual was reduced. We only print the first 2 digits, hence Gascoigne gives a 0. The
first number is the reduction rate in the current step, the second number is the average
reduction rate over all Newton steps.
On the right side, the convergence history of the linear multigrid solver is printed. Here,
the values indicate that in newton step 1 we needed {3} iteration of the multigrid solver. In
every of these three steps the error was reduced by a factor of [0.012] and after these three
steps the residual had the value 6.26e-08.
Finally Gascoigne prints the name of the output file.
visusimple u.00002.vtk in the terminal window.
The parameter file
In this directory you can also specify the parameter-file
File: FIRSTEXAMPLE/src/first.param
// B l o c k Loop
niter 5
1
2
3
// B l o c k Mesh
d im e n s i o n 2
4
5
10
You can visualize it by calling
1.2 A minimal example solving a partial differential equation
6
7
gridname s q u a r e . i n p
prerefine 3
8
9
10
// B l o c k BoundaryManager
dirichlet 4 1 2 3 4
11
12
13
14
15
dirichletcomp
dirichletcomp
dirichletcomp
dirichletcomp
1
2
3
4
1
1
1
1
0
0
0
0
16
17
// B l o c k Nix
This minimal example tells Gascoigne the dimension of the problem, which domain to use,
where to apply boundary data.
When starting the program, Gascoigne reads in a parameter file. Here, different parameters
controlling the behaviour of the code can be supplied. For instance in //Block Loop of this
file you can change the values of niter, which tells, how many times Gascoigne solves the
equation and refine the mesh.
In //Block Mesh we supply the mesh file as well as the number of refinement steps realized
before solving the first problem. You can change this parameter in prerefine. Thoe //Block
BoundaryManager will be explained in detail in Section 3.2
11
1 Introduction
12
2 The parameter file
Gascoigne can process a parameter file during runtime. Here, we can set some parameters
and options that can be changed at run-time without need for recompilation of the program.
For the file to be processed, go to the directory containing the parameter file and call
1
. . / b i n / Test t e s t . param
The data in the param-file is organized in blocks which are read at different points of the
code. A sample parameter file could look like: File: TEST/src/test.param
1
2
3
4
// B l o c k S e t t i n g s
month
10
year
2009
name
FirstTestOfGascoigne
5
6
7
// B l o c k S o l v e r
...
8
9
// B l o c k Nix
Each block has to be initialized by the line
1
// B l o c k
Note that the slashes in front of the blocks are not treated as a comment here. The parameter
file has to end with a dummy-block, e.g. //Block Nix.
In Gascoigne, the parameter file is stored in the class ParamFile. Reading values from
the parameter-file is done using the classes FileScanner and DataFormatHandler. These
classes are declared in the header File: filescanner.h . Here is a modified File: main.cc
reading values from the parameter file:
File: TEST/src/main.cc
1
2
3
4
#include
#include
#include
#include
” s t d l o o p . h”
” f i l e s c a n n e r . h”
” p a r a m f i l e . h”
5
6
7
using namespace G a s c oi g n e ;
using namespace s t d ;
8
13
2 The parameter file
9
10
11
12
13
14
15
int main ( int argc , char ∗∗ argv )
{
i f ( argc <2)
{ c e r r << ” c a l l : Test p a r a m f i l e ! ! ” << e n d l ;
abort ( ) ; }
ParamFile p a r a m f i l e ;
p a r a m f i l e . SetName ( argv [ 1 ] ) ;
16
int month , y e a r ;
s t r i n g name ;
DataFormatHandler DFH;
DFH. i n s e r t ( ”month” ,&month ) ;
DFH. i n s e r t ( ” y e a r ” , &year , 2009 ) ;
DFH. i n s e r t ( ”name” ,&name , ” t e s t ” ) ;
F i l e S c a n n e r FS (DFH, &p a r a m f i l e , ” S e t t i n g s ” ) ;
c o u t << month << ” / ” << y e a r << ” \ t ” << name << e n d l ;
17
18
19
20
21
22
23
24
25
}
To read in one block, we first initialize a DataFormatHandler (line 19). Next, we tell the
DataFormatHandler which keywords to look for using the function insert. The first argument of the insert function is the keyword we are searching, the second argument the
variable that its value will be written to. Optionally, we can specify a default value as a
third argument which is used if the keyword is not present in the parameter file. Finally, we
initialize an object FileScanner passing the DataFormatHandler, a parameter file and the
block to be read. Note that only this block is read! To read from a different block within
one file, we have to define a second DataFormatHandler.
14
3 Description of the problem
File: GAS/src/Problem/problemdescriptorbase.h
The class ProblemDescriptor provides all information necessary to set up an application.
This includes all problem data like: equation, boundary data, right hand side, initial condition, etc. The most important function is ProblemDescriptor::BasicInit(), where the
problem is defined by setting pointers to the describing classes. A minimal example might
be:
1
2
3
4
5
void L o c a l P r o b l e m D e s c r i p t o r : :
{
GetEquationPointer ( )
GetRightHandSidePointer ( )
GetDirichletDataPointer ()
B a s i c I n i t ( const ParamFile ∗ p f )
= new L o c a l E q u a t i o n ;
= new LocalRightHandSide ;
= new L o c a l D i r i c h l e t D a t a ;
6
ProblemDescriptorBase : : B a s i c I n i t ( pf ) ;
7
8
}
For an problem to be well defined, at least the equation and the Dirichlet data have to be
specified. All possible pointers to be set are listed in
File: GAS/src/Problem/problemdescriptorbase.h
For the most basic problems Gascoigne provides predefined classes, e.g.
1
2
3
G e t E q u a t i o n P o i n t e r ( ) = new L ap l a ce 2 d ;
GetRightHandSidePointer ( ) = new OneRightHanSideData ( 1 ) ;
G e t D i r i c h l e t D a t a P o i n t e r ( ) = new Z e r o D i r i c h l e t D a t a ;
which defines the Laplace problem with homogeneous Dirichlet data and constant right-hand
side f = 1.
In the following we describe the classes needed to specify the most common ingredients of a
problem.
3.1 The right hand side
File: GAS/src/Interface/domainrighthandside.h
The interface class to specify a function f : Ω → R as right hand side for the problem is
DomainRightHandSide. Let f (x, y) = sin(x) cos(y) be the right hand side:
15
3 Description of the problem
c l a s s LocalRightHandSide : public DomainRightHandSide
{
public :
int GetNcomp ( ) const { return 1 ; }
1
2
3
4
5
s t d : : s t r i n g GetName ( ) const { return ” L o c a l Right Hand S i d e ” ; }
6
7
double operator ( ) ( int c , const Vertex2d& v ) const
{
return s i n ( v . x ( ) ) ∗ c o s ( v . y ( ) ) ;
}
8
9
10
11
}
12
The most important function of this class is double operator()(int c, const Vertex2d&
v) const, here, the value of the right hand side function f at point v ∈ Ω is returned. The
class Vertex2d is explained in File: GAS/src/Common/vertex.h . For three dimensional
problems the operator double operator()(int c, const Vertex3d& v) const has to be
written.
The second parameter int c as well as the function int GetNcomp() const is only needed
if we solve a system of equations. This is explained in detail in Section 3.5.
A function called GetName() is necessary in every class used in the ProblemDescriptor to
provide a label.
3.2 Boundary Data
In the geometry file, a certain color is assigned to every boundary part of the mesh, see
Section 5. You can visualize meshes (inp-files) with the program visusimple. Here you can
also look at the boundary data, use Select->Scalar. In the parameter file we specify all
boundary colors where Dirichlet boundary conditions are used.
To assign Dirichlet boundary condition to a certain boundary color, the parameter file has
to contain:
1
2
3
4
// B l o c k BoundaryManager
dirichlet
2 4 8
dirichletcomp 4 1 0
dirichletcomp 8 1 0
This definition lets the boundaries with colors 4 and 8 have Dirichlet condition. The parameter dirichlet is given as a vector: the first value is the number of Dirichlet colors, then
the colors follow. If only the color 4 is picked for Dirichlet it would look like:
1
2
3
// B l o c k BoundaryManager
dirichlet
1 4
dirichletcomp 4 1 0
16
3.2 Boundary Data
We have to specify a second parameter dirichletcomp for every boundary color used with
Dirichlet values. Here, we choose the components of the solution where Dirichlet conditions
shall be used. These values are again given as vectors: First the color is specified, then the
number of solution components that have Dirichlet values on this color. Finally the solution
components are listed. For scalar equations with only one component, the last two values
are always 1 0. See Section 3.5 for using Dirichlet colors for systems of partial differential
equations with more than one solution component.
Every boundary color, that is not listed as Dirichlet boundary uses natural boundary conditions given due to integration by parts. For the Laplace Equation
(−∆u, φ)Ω + h∂n u, φiΓ = (∇u, ∇φ)Ω ,
this is the homogeneous Neumann condition ∂n u = 0. How to use non-homogeneous Neumann is explained in Section 3.2.2 and details on Robin-boundary conditions are given in
Section 3.2.3.
3.2.1 Dirichlet boundary data
File: src/Problem/diricheltdata.h
Assume, that the colors 4 and 8 are picked as Dirichlet boundary colors in the parameter
file.
The Dirichlet values to be used for the boundaries with colors specified in the parameter file
are described in the class DirichletData, see File: GAS/src/Interface/dirichletdata.h
Here we show an example how to define a class LocalDirichletData that sets Dirichlet
values
1
2
#include ” d i r i c h l e t d a t a . h”
3
4
5
6
7
8
9
10
11
12
13
14
c l a s s L o c a l D i r i c h l e t D a t a : public D i r i c h l e t D a t a
{
public :
s t d : : s t r i n g GetName ( ) const { return ” L o c a l D i r i c h l e t Data ; ” }
void operator ( ) ( DoubleVector& b , const Vertex2d& v , int c o l ) const
{
b . zero ( ) ;
i f ( c o l ==4) b [ 0 ] = 0 . 0 ;
i f ( c o l ==8) b [ 0 ] = v . x ( ) + v . y ( ) ;
}
}
This example sets the boundary value to zero on the boundary with color 4 and to the
function g(x, y) = x+y on color 8. The values are not returned as in the case of the right hand
side, but written in the vector b. The value b[i] defines the Dirichlet value of the i-th solution
17
3 Description of the problem
component. For scalar pde’s we only set b[0]. The parameter Vertex2d& v gives the
coordinate and the parameter col gives the color of the node. In three dimensions, the same
operator exists using a Vertex3d& v to indicate the coordinate. In the ProblemDescriptor
we need to set a pointer to this new class: File: TEST/src/problem.h
c l a s s P r o b l e m D e s c r i p t o r : public P r o b l e m D e s c r i p t o r B a s e
{
// . . .
void B a s i c I n i t ( const ParamFile ∗ p a r a m f i l e )
{
// . . .
GetDirichletDataPointer () =
new L o c a l D i r i c h l e t D a t a ;
}
// . . .
};
1
2
3
4
5
6
7
8
9
10
11
3.2.2 Neumann conditions
For non-homogenous Neumann conditions of the type
h∂n u, φiΓN = hg N , φiΓN ,
we have to add the term (g N , φ) to the right hand side. In Gascoigne, this Neumann
right hand side is derived from the class BoundaryRightHandSide specified in the File:
GAS/src/Interface/boundaryrighthandside.h . An example for a 2d scalar problem is
given by
#i n c l u d e ” b o u n d a r y r i g h t h a n d s i d e . h”
1
2
c l a s s LocalBoundaryRightHandSide : public BoundaryRightHandSide
{
public :
int GetNcomp ( ) const { return 1 ; }
s t r i n g GetName ( ) const { return ” L o c a l B−RHS” ; }
3
4
5
6
7
8
double operator ( ) ( int c , const Vertex2d& v , const Vertex2d& n ,
int c o l o r ) const
{
i f ( c o l o r ==0) return 1 . 0 ;
i f ( c o l o r ==1) return v . x ( ) ∗ n . x ( ) + v . y ( ) ∗ n . y ( ) ;
}
9
10
11
12
13
14
15
};
16
Here, we define Neumann data on two different parts of the boundary with colors 0 and 1:
g0 (x) = 1,
18
g1 (x) = n(x) · x,
3.3 Definition of the partial differential equations
where n(x) is the outward unit-normal vector in the point x.
Boundary right hand sides need to be specified in the parameter file in //Block BoundaryManager.
Otherwise these boundary terms are not taken into account: File: TEST/src/test.param
1
2
3
4
// B l o c k BoundaryManager
dirichlet
1 4
dirichletcomp 4 1 0
righthandside 2 0 1
The new parameter righthandside specifies a vector. Here it tells Gascoigne that 2 boundary colors have Neumann data. These are the colors 0 and 1. In the ProblemDescriptor
we need to set a pointer to this new class: File: TEST/src/problem.h
1
2
3
4
5
6
7
8
9
10
11
c l a s s P r o b l e m D e s c r i p t o r : public P r o b l e m D e s c r i p t o r B a s e
{
// . . .
void B a s i c I n i t ( const ParamFile ∗ p a r a m f i l e )
{
// . . .
GetBoundaryRightHandSidePointer ( ) =
new LocalBoundaryRightHandSide ;
}
// . . .
};
3.2.3 Robin conditions
Robin boundary data includes very general conditions to be fulfilled on the boundary of the
domain. We can have:
hG(u), φiΓR = 0,
where G(·) can be some operator, e.g. G(u) = ∂n u + u2 . Hence, Robin boundary data
means that we have an additional equation that is valid on the boundary. We will explain
this concept in detail after explaining how to specify equations. See Section 3.3.2.
3.3 Definition of the partial differential equations
File: GAS/src/Interface/equation.h
Let the pde to be solved be given in the weak formulation
a(u)(φ) = (f, φ)
∀φ,
19
3 Description of the problem
where a(·)(·) is a semi-linear form, linear in the second argument. Gascoigne solves every
problem with a Newton method (also linear problems). With an initial guess u0 updates wk
are defined by
a0 (uk )(wk , φ) = (f, φ) − a(uk )(φ) ∀φ,
uk+1 := uk + wk .
The Jacobi matrix is the matrix of the directional derivatives of a(·)(·) in the point uk and
defined by
d
a0 (u)(w, φ) := a(u + sw)(φ)
.
ds
s=0
For linear problems, it holds
a(u)(w, φ) = a(w)(φ).
To solve, Gascoigne needs to know about the form a(uk )(φ) and its derivative, the matrix
a0 (uk )(wk , φ). Form and matrix are given in the class Equation:
1
2
3
c l a s s L o c a l E q u a t i o n : public Equation
{
public :
4
int GetNcomp ( ) const ;
s t r i n g GetName ( ) const ;
5
6
7
void p o i n t ( double h , const Vertex2d& v ) const ;
void p o i n t ( double h , const Vertex3d& v ) const ;
8
9
10
void Form ( V e c t o r I t e r a t o r b , const FemFunction& U,
T e s t F u n c t i o n& N) const ;
void Matrix ( EntryMatrix& A, const FemFunction& U,
const T e s t F u n c t i o n& M, const T e s t F u n c t i o n& N) const ;
11
12
13
14
15
16
};
The first function GetNcomp() returns the number of solution components for systems of
partial differential equations. For scalar equations, this function returns 1. GetName() is
a label for the equation. The most important functions are Form(b,U,N), which defines
a(u)(φ) and Matrix(A,U,M,N) which gives the derivative a0 (u)(w, φ). The parameters b and
A are the return values, U is the last approximation u, N is the test-function φ and M stands
for the direction w. The function point() is meant to set parameters depending on the
mesh size h or on the coordinate v. point is called before each call of Form() or Matrix()
and can be used to set local variables.
For the Laplace equation, the implementation is given by
1
2
3
void Form ( V e c t o r I t e r a t o r b , const FemFunction& U,
const , T e s t F u n c t i o n& N) const
{
20
3.3 Definition of the partial differential equations
b [ 0 ] += U [ 0 ] . x ( ) ∗ N. x ( ) + U [ 0 ] . y ( ) ∗ N. y ( ) ;
4
5
}
6
7
8
9
10
11
void Matrix ( EntryMatrix& A, const FemFunction& U,
const T e s t F u n c t i o n& M, const T e s t F u n c t i o n& N) const
{
A( 0 , 0 ) += M. x ( ) ∗ N. x ( ) + M. y ( ) ∗ N. y ( ) ;
}
The class TestFunction describes the values and derivatives of a discrete function in a
certain point. By N.m() the value is accessed, by N.x(), N.y() and N.z() the directional
derivatives. The class FemFunction is a vector of TestFunctions and used for the solution
function u. For systems of partial differential equations, the index gives the number of the
solution component. For scalar equations, it is always U[0].m().
3.3.1 Nonlinear equations
As example we now consider the nonlinear partial differential equation given by
X
a(u)(φ) = (∇u, ∇φ) + (h∇u, ∇ui, φ), hx, yi :=
xi yi .
The form is given by
1
2
3
4
5
6
void Form ( V e c t o r I t e r a t o r b , const FemFunction& U,
const , T e s t F u n c t i o n& N) const
{
b [ 0 ] += U [ 0 ] . x ( ) ∗ N. x ( ) + U [ 0 ] . y ( ) ∗ N. y ( )
+ (U [ 0 ] . x ( ) ∗ U [ 0 ] . x ( ) + U [ 0 ] . y ( ) ∗ U [ 0 ] . y ( ) ) ∗ N.m( ) ;
}
To define the matrix we first build the derivative:
d
a(u + sw)(φ)|s=0
ds
d
d
=
∇(u + sw), ∇φ +
h∇(u + sw), ∇(u + sw)i, φ |s=0
ds
ds
= (∇w, ∇φ) + (h∇(u + sw), ∇wi, φ) + (h∇w, ∇(u + sw)i, φ) |s=0
a0 (u)(w, φ) =
= (∇w, ∇φ) + 2(h∇u, ∇wi, φ).
The matrix function is then given by
1
2
3
4
5
6
void Matrix ( EntryMatrix& A, const FemFunction& U,
const T e s t F u n c t i o n& M, const T e s t F u n c t i o n& N) const
{
A( 0 , 0 ) += M. x ( ) ∗ N. x ( ) + M. y ( ) ∗ N. y ( ) ;
A( 0 , 0 ) += 2 . 0 ∗ (U [ 0 ] . x ( ) ∗ M. x ( ) + U [ 0 ] . y ( ) ∗ M. y ) ∗ N.m( ) ;
}
21
3 Description of the problem
3.3.2 Equations on the boundary & Robin boundary data
Robin boundary data is realized by setting an equation on the boundary. First, we need to
tell Gascoigne to include boundary equations. Herefore, section //Block BoundaryManager
of the parameter file needs to know on what boundary colors we set an equation: File:
TEST/src/test.param
// B l o c k BoundaryManager
dirichlet
1 4
dirichletcomp 4 1 0
equation
1 1
1
2
3
4
Here we tell Gascoigne to use the equation on 1 color, the color 1. The actual equation is
of type BoundaryEquation defined in the File: GAS/src/Interface/boundaryequation.h
The boundary equation is declared as the equation: you need to specify the Form and the
Matrix. Major difference is that these two function get the outward unit normal vector
and the color of the boundary line in addition. Boundary equations are then defined in the
ProblemDescriptor by setting: File: TEST/src/problem.h
c l a s s P r o b l e m D e s c r i p t o r : public P r o b l e m D e s c r i p t o r B a s e
{
// . . .
void B a s i c I n i t ( const ParamFile ∗ p a r a m f i l e )
{
// . . .
GetBoundaryEquationPointer ( ) =
new LocalBoundaryEquation ;
}
// . . .
};
1
2
3
4
5
6
7
8
9
10
11
We next give an example of a boundary equation which will realize the Robin data ∂n u+u2 =
0, see File: TEST/src/problem.h
1
#include ” b o u n d a r y e q u a t i o n . h”
2
3
4
5
6
c l a s s LocalBoundaryEquation : public BoundaryEquation
{
public :
7
int GetNcomp ( ) const { return 1 ; }
s t r i n g GetName ( ) const { return ” L o c a l B−EQ” ; }
8
9
10
11
12
13
void pointboundary ( double h , const Vertex2d& v , const Vertex2d& n )
const
{}
14
22
3.4 Exact Solution and Evaluation of Errors
void Form ( V e c t o r I t e r a t o r b , const FemFunction& U,
T e s t F u n c t i o n& N, int c o l ) const
{
i f ( c o l ==1) b [ 0 ] += U [ 0 ] .m( ) ∗ U [ 0 ] .m( ) ∗ N.m( ) ;
}
15
16
17
18
19
20
21
22
23
24
25
26
void Matrix ( EntryMatrix& A, const FemFunction& U,
const T e s t F u n c t i o n& M, const T e s t F u n c t i o n& N, int c o l )
const
{
i f ( c o l ==1) A( 0 , 0 ) += 2 . 0 ∗ U [ 0 ] .m( ) ∗ M.m( ) ∗ N.m( ) ;
}
27
28
};
3.4 Exact Solution and Evaluation of Errors
If an analytic solution is known it can be added to the ProblemDescriptor and then be
used for evaluating the error in the L2 , H 1 and L∞ norms. The operator in the class
ExactSolution is used to define the solution function. The parameter int c specifies the
solution component for systems of pdes.
File: ../src/problem.h
1
#i n c l u d e ” e x a c t s o l u t i o n . h”
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
c l a s s MyExactSolution : public E x a c t S o l u t i o n
{
public :
s t d : : s t r i n g GetName ( ) const { return ”My e x a c t s o l u t i o n ” ; }
double operator ( ) ( int c , const Vertex2d& v ) const
{
return v . x ( ) ∗ v . y ( ) ;
}
};
c l a s s Problem : public P r o b l e m D e s c r i p t o r B a s e
{
...
public :
void B a s i c I n i t ( const ParamFile ∗ p f )
{
...
G e t E x a c t S o l u t i o n P o i n t e r ( ) = new MyExactSolution ;
}
};
23
3 Description of the problem
The function ComputeGlobalErrors which calculates the L2 , L∞ and H 1 norm error has to
be called in the Loop after solving the equation with: File: ../src/myloop.cc
1
...
Solve (u , f ) ;
ComputeGlobalErrors ( u ) ;
2
3
4
...
3.5 Systems of partial differential equations
One of the main features of Gascoigne is the efficient solution of large systems of partial
differential equations. We call the number of equations in a system (and thus the number
of solution variables) the number of components, denoted by ncomp in Gascoigne. Nearly all
classes have the number of components as member. As example we discuss the two-equation
system in Ω:
Γ8
(∇u0 , ∇φ0 ) + (u0 u1 , φ0 ) = (f0 , φ0 ),
(∇u1 , ∇φ1 ) + (u0 sin(u1 ), φ1 ) = (f1 , φ1 ),
u0 = 0
on Γ4 ∪ Γ8 ,
u1 = 1
on Γ4 ,
∂n u1 = 0
on Γ8 .
Ω
Γ4
Γ4
Γ8
Several modifications are necessary:
3.5.1 Modification of the right hand side
See 3.1 for comparison. The right hand side has to know the number of equations, in this
case:
1
int GetNcomp ( ) const { return 2 ; }
The operator for specifying the functions f0 and f1 gets the current component as first
argument. As an example, let us set f0 (x, y) = 5 and f1 (x, y) = xy:
1
2
3
4
5
6
double operator ( ) ( int c , const Vertex2d& v ) const
{
i f ( c==0) return 5 ;
i f ( c==1) return v . x ( ) ∗ v . y ( ) ;
return 0 . 0 ;
}
24
3.5 Systems of partial differential equations
3.5.2 Modification of Dirichlet data
For Dirichlet data two separate issues have to be considered: first the declaration in the
parameter file has to be modified. Second, the function specifying the Dirichlet data has to
be altered. In this example, where the component u0 has Dirichlet data on both and u1 just
on one boundary color, we set:
1
2
3
4
// B l o c k BoundaryManager
dirichlet 2
4 8
dirichletcomp 4
2
0
dirichletcomp 8
1
0
1
This is to be read as: there are 2 dirichlet colors, the colors 4 and 8. Dirichlet color 4 is
applied to 2 components, components 0 and 1. And so on. Gascoigne uses 0 as index for the
first component. The implementation of the boundary data is as follows:
1
2
3
4
5
6
7
8
9
10
11
12
c l a s s L o c a l D i r i c h l e t D a t a : public D i r i c h l e t D a t a
{
public :
s t d : : s t r i n g GetName ( ) const { return ” L o c a l D i r i c h l e t Data ; ” }
void operator ( ) ( DoubleVector& b , const Vertex2d& v , int c o l ) const
{
b . zero ( ) ;
i f ( c o l ==4) b [ 0 ] = 0 . 0 ;
i f ( c o l ==4) b [ 1 ] = 1 . 0 ;
i f ( c o l ==8) b [ 0 ] = 0 . 0 ;
}
};
Here, the values are sorted into a vector whose components correspond to the components
of the solution u.
3.5.3 Modification for the equation, form and matrix
Like the right hand side, the equation needs to know about the number of solution components given in the header file as
1
int GetNcomp ( ) const { return 2 ; }
In both functions, Form and Matrix the last approximate solution U is given as a parameter.
This parameter is a vector and U[c] specifies component c of the solution. The modifications
for the Form are straight-forward, the values for i-th equation are written to the component
b[i]. The test function N is the same for all equations. We give an example:
1
2
3
void Form ( V e c t o r I t e r a t o r& b , const FemFunction& U,
const T e s t F u n c t i o n& N) const
{
25
3 Description of the problem
b [ 0 ] += U [ 0 ] . x ( ) ∗ N. x ( ) + U [ 0 ] . y ( ) ∗ N. y ( ) ;
b [ 0 ] += U [ 0 ] .m( ) ∗ U [ 1 ] .m( ) ∗ N.m( ) ;
4
5
6
b [ 1 ] += U [ 1 ] . x ( ) ∗ N. x ( ) + U [ 1 ] . y ( ) ∗ N. y ( ) ;
b [ 1 ] += U [ 0 ] .m( ) ∗ s i n (U [ 1 ] .m( ) ) ∗ N.m( ) ;
7
8
9
}
The Matrix works in a simular way. We assemble a local matrix of size ncomp times ncomp,
where A(i,c) specifies the derivative of equation i with respect to the component U[c]. In
this example:
d
(∇(u0 + sw), ∇φ0 ) + ((u0 + sw)u1 , φ0 )
ds
d
= (∇u0 , ∇φ0 ) + (u0 (u1 + sw), φ0 ) s=0
ds
d
= (∇u1 , ∇φ1 ) + ((u0 + sw) sin(u1 ), φ1 )
ds
d
= (∇u1 , ∇φ1 ) + (u0 sin(u1 + sw), φ1 )
ds
A00 =
A01
A10
A11
s=0
= (∇w, ∇φ0 ) + (wu1 , φ0 ).
= (u0 w, φ0 )
= (w sin(u1 ), φ1 ),
= (∇w, ∇φ1 ) + (u0 cos(u1 )w, φ1 ),
and the implementation:
1
2
3
4
5
6
void Matrix ( EntryMatrix& A, const FemFunction& U,
const T e s t F u n c t i o n& M, const T e s t F u n c t i o n& N) const
{
A( 0 , 0 ) += M. x ( ) ∗ N. x ( ) + M. y ( ) ∗ N. y ( ) ;
A( 0 , 0 ) += M.m( ) ∗ U [ 1 ] .m( ) ∗ N.m( ) ;
A( 0 , 1 ) += U [ 0 ] .m( ) ∗ M.m( ) ∗ N.m( ) ;
7
A( 1 , 0 ) += M.m( ) ∗ s i n (U [ 0 ] .m( ) ) ∗ N.m( ) ;
A( 1 , 1 ) += M. x ( ) ∗ N. x ( ) + M. y ( ) ∗ N. y ( ) ;
A( 1 , 1 ) += U [ 0 ] .m( ) ∗ c o s (U [ 1 ] .m( ) ) ∗ M.m( ) ∗ N.m( ) ;
8
9
10
11
}
26
4 Time dependent problems
Gascoigne provides a simple technique for the modeling of time-dependent partial differential
equations like
∂t u + A(u) = f in [0, T ] × Ω,
with a spatial differential operator A(u). Time-discretization is accomplished with Rothe’s
method by first discretizing in time with simple time-stepping methods. For time-dependent
problems, we must provide some additional information:
• We must supply an initial value u(x, 0) = u0 (x).
• All problem data like the right hand side, boundary values and differential operator
can depend on the time t ∈ [0, T ].
• We must describe the time-stepping scheme and choose a time-step size.
4.1 Time discretization of the heat equation
Here, we examplarily demonstrate the time discretization of the heat equation
∂t u − ∆u = f in [0, T ] × Ω,
with Dirichlet boundary
u = g on [0, T ] × Ω,
and the initial value
u(·, 0) = u0 (·) on Ω.
In Gascoigne, time discretization is based on the general θ-scheme
1 m
1−θ
1
1−θ
(u , φ) + a(um , φ) = (f (tm ), φ) +
(f (tm−1 ), φ) + (um−1 , φ) −
a(um−1 , φ),
θk
θ
θk
θ
where θ ∈ (0, 1]. For θ = 21 , this scheme corresponds to the Crank-Nicolson method, for
θ = 1 to the implicit Euler method.
The class StdTimeLoop includes a simple loop for time-stepping problems and replaces the
StdLoop:
27
4 Time dependent problems
1
2
3
4
5
6
7
8
void StdTimeLoop : : run ( const s t d : : s t r i n g& p r o b l e m l a b e l )
{
[...]
for ( i t e r =1; i t e r <= n i t e r ;
i t e r ++)
{
[...]
t i m e i n f o . SpecifyScheme ( i t e r ) ;
TimeInfoBroadcast ( ) ;
9
// RHS
G e t M u l t i L e v e l S o l v e r ()−> G e t S o l v e r ()−>GetGV( f ) . z e r o ( ) ;
G e t M u l t i L e v e l S o l v e r ()−> G e t S o l v e r ()−>TimeRhsOperator ( f , u ) ;
G e t M u l t i L e v e l S o l v e r ()−> G e t S o l v e r ()−>TimeRhs ( 1 , f ) ;
10
11
12
13
14
timeinfo . iteration ( i t e r );
TimeInfoBroadcast ( ) ;
15
16
17
G e t M u l t i L e v e l S o l v e r ()−> G e t S o l v e r ()−>TimeRhs ( 2 , f ) ;
18
19
// SOLVE
SolveTimePrimal ( u , f ) ;
20
21
22
[...]
23
}
24
25
}
Every time-step is treated as a stationary problem. The differential operator evaluated at
the old time step is known in advance and thus, Gascoigne automatically includes it when
assembling the right-hand side, if θ < 1. The function TimeInfoBroadcast() announces the
current time t and the time-step k.
The class StdTimeSolver is used instead of StdSolver. It has to be invoked by the option solver instat in block //Block MultiLevelSolver of the parameter file. Gascoigne
assembles the time-dependent matrix
1
M + A,
θk
where M is the mass matrix and A the stiffness matrix which has to be specified as usual
in the class Equation.
1
2
3
void StdTimeSolver : : AssembleMatrix ( const V e c t o r I n t e r f a c e& gu , double d )
{
S t d S o l v e r : : AssembleMatrix ( gu , d ) ;
4
double s c a l e = d / ( d t ∗ t h e t a ) ;
GetMatrix()−> A d d M a s s W i t h D i f f e r e n t S t e n c i l ( GetMassMatrix ( ) ,
GetTimePattern ( ) , s c a l e ) ;
5
6
7
28
4.1 Time discretization of the heat equation
8
}
First, the main part regarding a(·, ·) is assembled. Then, the part regarding the timederivative is added with the proper scaling depending on the time-step k and parameter θ.
The parameters θ and the time-step k are provided in the parameter file within the //Block
Loop
1
2
3
4
5
// B l o c k
scheme
theta
dt
neuler
Loop
Theta
0.5
0.01
0
By the option neuler we can tell Gascoigne to start with a number of neuler initial steps
using the implicit Euler scheme. In this way irregular initial data may be smoothened which
can be necessary for stability reasons. By setting neuler 0, no Euler steps are used and the
loop directly starts with the θ-time-stepping scheme.
The function GetTimePattern() tells Gascoigne, on which solution components the massmatrix will act. This TimePattern is specified in the Equation by a further function:
1
2
3
4
void HeatEquation : : SetTimePattern ( TimePattern& TP) const
{
TP( 0 , 0 ) = 1 . ;
}
TimePattern TP is a matrix of size ncomp times ncomp. Typically, this matrix is the diagonal
unit-matrix. See Section 4.3 for special choices of the TimePattern.
Time dependent problems usually are initial-boundary value problems and require some
initial data, like
u(x, 0) = u0 (x),
at time t = 0. The initial data is specified via the function GetInitialConditionPointer()
in the class ProblemDescriptor. The initial condition can be set in the same way as the
right-hand side:
1
2
3
4
5
c l a s s M y I n i t i a l : public DomainRightHandSide
{
public :
s t d : : s t r i n g GetName ( ) const { return ” I n i t i a l C o n d i t i o n ” ; }
int GetNcomp ( ) const { return 1 ; }
6
double operator ( ) ( int c , const Vertex2d& v ) const
{
return v . x ( ) ∗ v . y ( ) ;
}
7
8
9
10
11
};
12
29
4 Time dependent problems
13
[...]
14
15
16
17
18
19
20
21
22
23
24
25
c l a s s P r o b l e m D e s c r i p t o r : public P r o b l e m D e s c r i p t o r B a s e
{
public :
s t d : : s t r i n g GetName ( ) const { return ”Time−Dependent ” ; }
void B a s i c I n i t ( const ParamFile ∗ p f )
{
[...]
G e t I n i t i a l C o n d i t i o n P o i n t e r ( ) = new M y I n i t i a l ;
ProblemDescriptorBase : : B a s i c I n i t ( pf ) ;
}
};
4.2 Time dependent problem data
The problem data like boundary data and right hand side can depend on time. Let us for
example consider the right hand side
f (x, t) = sin(πt)(1 − x2 )(1 − y 2 ).
All data classes like DomainRightHandSide, DirichletData or Equation are derived from
the class Application. We can access the current time and the time step via the functions
double GetTime() const and double GetTimeStep() const.
1
2
3
4
5
c l a s s RHS : public DomainRightHandSide
{
public :
int GetNcomp ( ) const { return 1 ; }
s t r i n g GetName ( ) const { return ”RHS” ; }
6
double operator ( ) ( int c , const Vertex2d& v ) const
{
double t = GetTime ( ) ;
return s i n ( M PI ∗ t ) ∗ ( 1 . 0 − v . x ( ) ∗ v . x ( ) ) ∗ (1.0 − v . y ( ) ∗ v . y ( ) ) ;
}
7
8
9
10
11
12
};
4.3 Non standard time-discretizations
As example for a non-standard time-depending problem we consider the following system of
partial differential equations
∂t u1 + ∂t u2 − ∆u1 = f1
∂t u2 + u1 − ∆u2 = f2
30
4.3 Non standard time-discretizations
For the implementation of this system, Matrix, Form and SetTimePattern must be given
as:
1
2
3
void EQ : : Form ( . . . )
{
b [ 0 ] += U [ 0 ] . x ( ) ∗ N. x ( ) + U [ 0 ] . y ( ) ∗ N. y ( ) ;
4
b [ 1 ] += U [ 0 ] .m( ) ∗ N.m( ) ;
b [ 1 ] += U [ 1 ] . x ( ) ∗ N. x ( ) + U [ 1 ] . y ( ) ∗ N. y ( ) ;
5
6
7
}
8
9
10
11
void EQ : : Matrix ( . . . )
{
A( 0 , 0 ) += M. x ( ) ∗ N. x ( ) + M. y ( ) ∗ N. y ( ) ;
12
A( 1 , 0 ) += M.m( ) ∗ N.m( ) ;
A( 1 , 1 ) += M. x ( ) ∗ N. x ( ) + M. y ( ) ∗ N. y ( ) ;
13
14
15
}
16
17
18
19
20
21
22
void EQ : : SetTimePattern ( TimePattern& TP) const
{
TP( 0 , 0 ) = 1 . ;
TP( 0 , 1 ) = 1 . ;
TP( 1 , 1 ) = 1 . ;
}
31
4 Time dependent problems
32
5 Mesh handling
Mesh handling includes definition of the domain Ω, refinement of meshes, managing of curved
boundaries, input and output of meshes.
The mesh to be used in the calculation is indicated in the parameter file. In the Section
//Block Mesh the mesh to be read is specified:
1
// B l o c k Mesh
2
3
4
5
d im en si o n
gridname
prerefine
2
mesh . i n p
3
Gascoigne distinguishes two types of meshes: inp-files are coarse meshes describing the domain Ω. They are user-specified and the file format is explained in Section 5.1. These meshes
are usually as coarse as possible. Refinement is applied during run-time. The parameter
prerefine declares global refinements to be executed before the first solution cycle. It is
also possible to write out and read in locally refines meshes. These meshes are given in the
gup-format:
1
// B l o c k Mesh
2
3
4
5
d im en si o n
gridname
prerefine
2
mesh . gup
0
During run-time all mesh handling like I/O and refinement is done by the class MeshAgent.
The MeshAgent is a member of the class Loop and is created in the Loop::BasicInit. The
MeshAgent takes care of the mesh refinement hierarchy HierarchicalMesh, it creates the
sequence of multigrid meshes GascoigneMultigridMesh and it finally provides the meshes
on every level used for the computation, the GascoigneMeshes. See the files with the same
names in File: GAS/src/Mesh/* .
If domains with curved boundaries are used, the shape definition is also done in the class
MeshAgent.
5.1 Definition of coarse meshes
Coarse meshes are specified in the inp-format. We indicate the nodes of the mesh, the
elements (that are quadrilaterals in two and hexahedrals in three dimensions) and finally
33
5 Mesh handling
the boundary elements. In two dimensions, these are the lines of quadrilaterals that touch
the boundary, in three dimensions we lists the quadrilaterals that are boundaries of the
domain and of elements. We never list interior lines (or quadrilaterals in three dimensions),
as these are not needed for the calculation. The following example is the coarse mesh for
the L-shaped domain:
8
0
1
2
3
4
5
6
7
0
1
2
0
1
2
3
4
5
6
7
11 0 0
0 -1 0
1 -1 0
-1 0 0
0 0 0
1 0 0
-1 1 0
0 1 0
1 1 0
0 quad
0 quad
0 quad
4 line
4 line
4 line
4 line
2 line
2 line
8 line
8 line
0
0
2
3
0
1
4
0
2
2
5
6
1 4 3
3 6 5
4 7 6
1
4
7
3
3
5
6
7
1
5
6
7
0
2
3
4
0
1
−1
−1
0
1
n7
n3
n6
n2
n2
n3
n4
n0
n1
n0
n5
n1
The file format is very simple: in the row, first the number of nodes in the coarse mesh is
given, followed by the sum of mesh elements plus elements on the boundary. The last three
values are always zero. This example has 8 nodes, 3 quads and 8 lines on the boundary.
After the header line all nodes of the mesh are declared in the format:
n x y z
where n is the index of the node (starting with zero) and x y z are the coordinates. We
always have to give all three coordinates (also in two dimensions).
The nodes are followed by the mesh elements. In two dimensions the format is:
n 0 quad n0 n1 n2 n3
where n is a consecutive index starting with 0. The second parameter is not used, thus
always 0. The mesh element, always quad is indicated by the third entry. The last four
values give the nodes describing the quadrilateral in counter-clockwise sense.
In three spatial dimension, the format is:
34
5.2 Curved boundaries
n 0 hex n0 n1 n2 n3 n4 n5 n6 n7
Here, first the nodes n0 to n3 in the front are given in counter-clock wise sense followed by
the nodes in the back.
Finally, all the boundary elements are given. In two dimensional meshes, these are all lines,
where both nodes are on the boundary. In three dimensions the boundary elements are all
quads with four nodes on the boundary. The format is:
n c line n0 n1
in two dimensions and
n c quad n0 n1 n2 n3
in three dimensions. n is again an iteration index. This index starts with 0. The value c
defines a color for the boundary element. These colors are necessary to distinguish between
different types of boundary conditions, see Sections 3.2.1 and 3.2.2. In two dimensions, the
order of the two nodes n0 n1 is arbitrary. In three dimensions however the four nodes have to
be given in a counter clock-wise sense seen from the outside of the domain! This is necessary
to get a unique definition of the normal vectors on the boundary. Thus, for the hex shown
in Picture above, the inp-file would be:
8
0
1
2
3
4
5
6
7
7
0
1
1
0
0
1
1
0
0
0
0
1
1
0
0
1
1
0
0
1
2
3
4
5
0 0
0
0
0
0
1
1
1
1
0
4
4
4
4
4
4
hex 0 1 2 3 4 5 6 7
quad 0 1 2 3
quad 1 5 6 2
quad 3 2 6 7
quad 4 0 3 7
quad 5 4 7 6
quad 4 5 1 0
5.2 Curved boundaries
A part of the boundary can be assigned a curved boundary function. The idea is easy: Let
the boundary with color col be a curved boundary. We specify a function Rd → R which
has the desired curved boundary as zero-iso-contour. For instance if the boundary Γ of the
domain should have the shape of a circle with radius r and midpoint (mx , my ), one possible
function fcircle would be
fcircle (x, y) = (x − mx )2 + (y − my )2 − r2 .
Whenever an element with boundary nodes on a line with color col is refined, the newly generated nodes on the boundary line are then pulled onto the zero-iso-contour of the boundary
function. This is done by solving for the root of fcircle with Newton’s method.
35
5 Mesh handling
Refine
Adjust Curved Boundary
In Gascoigne boundary functions are described by the class BoundaryFunction in
File: GAS/src/Mesh/BoundaryFunction . The circle is given by:
1
2
3
4
5
6
c l a s s C i r c l e : public BoundaryFunction <2>
{
double
r , mx , my ;
public :
C i r c l e ( double r , double mx, double my) :
{ }
r (r),
mx (mx) ,
my (my)
7
double operator ( ) ( const Vertex2d& v ) const
{
return ( v . x()− mx ) ∗ ( v . x()− mx )+(v . y()− my ) ∗ ( v . y()− my)−
}
8
9
10
11
12
}
The boundary function object has to be passed to the MeshAgent. For this, a derived mesh
agent class has to be specified and created in the Loop. The following minimal steps have to
be taken:
1
2
3
4
5
6
7
8
9
c l a s s LocalMeshAgent : public MeshAgent
{
Circle c i r c l e ;
public :
LocalMeshAgent ( ) : MeshAgent ( )
{
AddShape ( 4 , & c i r c l e ) ;
}
};
10
11
12
13
c l a s s LocalLoop : public StdLoop
{
public :
14
void B a s i c I n i t ( const ParamFile ∗ p a r a m f i l e ,
const ProblemContainer ∗ PC,
const F u n c t i o n a l C o n t a i n e r ∗ FC=NULL)
{
15
16
17
18
36
r∗
r;
5.3 Adaptive refinement of meshes
GetMeshAgentPointer ( ) = new LocalMeshAgent ( ) ;
StdLoop : : B a s i c I n i t ( p a r a m f i l e , PC, FC ) ;
19
20
}
21
22
};
In this example, the boundary with color 4 is used as curved boundary.
5.3 Adaptive refinement of meshes
For local refinement, the MeshAgent has the function refine nodes(nodes). The nvector
nodes contains indices of mesh nodes to be refined. If a node is to be refined, all elements
having this node as a corner node will be refined.
File: ../src/myloop.cc
1
2
3
4
5
6
7
8
9
10
for ( i t e r =1; i t e r < n i t e r ;++ i t e r )
{
...
Solve (u , f ) ;
...
// E s t i m a t e Error
n v e c t o r r e f ;
// p i c k nodes f o r r e f i n e m e n t i n r e f
GetMeshAgent()−> r e f i n e n o d e s ( r e f ) ;
}
5.3.1 Estimating the Energy-Error
The class DwrQ1Q2 contains a function for estimation of an energy-error. It can be used in
the following way:
File: ../src/myloop.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
#include ” dwrq1q2 . h”
...
for ( i t e r =1; i t e r < n i t e r ;++ i t e r )
{
...
Solve (u , f ) ;
...
// E s t i m a t e Error
DoubleVector e t a ;
DwrQ1Q2 dwr ( ∗ G e t M u l t i L e v e l S o l v e r ()−> G e t S o l v e r ( ) ) ;
dwr . EstimatorEnergy ( eta , f , u ) ;
...
}
37
5 Mesh handling
Now, for each node of the mesh the vector η contains a measure for the local contribution
to the total energy error.
5.3.2 Picking Nodes for Refinement
There are several methods for picking nodes for refinement. Common to all is to choose nodes
n with large estimator values eta[n]. The class MalteAdaptor uses a global optimization
scheme for picking nodes:
File: ../src/myloop.cc
1
2
3
4
5
6
7
8
9
10
11
#include ” m a l t e a d a p t o r . h ’ ’
...
for ( i t e r =1; i t e r < n i t e r ;++ i t e r )
{
...
Solve (u , f ) ;
...
// E s t i m a t e Error
DoubleVector e t a ;
DwrQ1Q2 dwr ( G e t M u l t i L e v e l S o l v e r ()−> G e t S o l v e r ( ) ) ;
dwr . EstimatorEnergy ( eta , f , u ) ;
12
// Pick Elements
n v e c t o r r e f n o d e s ;
MalteAdaptor A( p a r a m f i l e , e t a ) ;
A. r e f i n e ( r e f n o d e s ) ;
13
14
15
16
17
// R e f i n e
GetMeshAgent()−> r e f i n e n o d e s ( r e f n o d e s ) ;
18
19
}
20
38
6 Flow problems and stabilization
In order to discretize the Stokes or Navier-Stokes systems one can either use inf-sup- stable
finite element pairs for pressure and velocity which fulfill an inf-sup condition or include
stabilization terms in the equations that ensure uniqueness of the solution and stability.
Gascoigne always uses equal-order elements which are not inf-sup stable, in combination
with different pressure stabilization techniques. If in the Navier-Stokes case the convective
term is dominating (high Reynolds number) further stabilization is required.
6.1 The function point
To set element-wise constants in Gascoigne, the function point is needed:
1
2
3
void Equation : :
p o i n t ( double h , const FemFunction& U, const Vertex2d& v ) const
{
4
5
}
This function is called before every call of Matrix and Form and here the local mesh-size
h and the coordinates of the node v are given as a parameter. That way one can specify
certain criteria for the equation based on location or local mesh width.
To stabilize a given equation or system of equations (like Stokes and Navier-Stokes) one
usually modifies it by adding certain stabilizing terms. These terms often depend on the
local mesh width and have therefore to be defined in this function.
Since the function point is defined as const no class members can be altered. Thus, every
variable which shall be modified here has to be defined as mutable in the header file:
1
mutable double v a r i a b l e ;
6.2 Stabilization of the Stokes system
The Stokes equations are given by
∇ · v = 0,
−ν∆v + ∇p = f in Ω.
39
6 Flow problems and stabilization
The simplest stabilization technique is the addition of an artificial diffusion. Stability is
ensured by modifying the divergence equation to
X
(∇ · v, ξ) +
αK (∇p, ∇ξ)K ,
K∈Ωh
where αK is an element-wise constant that can be chosen as
αK = α0
h2K
.
6ν
To get the local cell size hK of an element K, we use the function point which is called
before each call of Form and Matrix. In point we can calculate the value αK and save it to
a local variable:
void S t o k e s : :
p o i n t ( double h , const FemFunction& U, const Vertex2d& v ) const
{
alphaK =
a l p h a 0 ∗ h∗h / 6 . 0 / n u ;
}
1
2
3
4
5
The parameter alpha0 can be read from the parameter file. The variable
be defined as mutable in the header file:
alphaK ;
mutable double
1
alphaK has to
6.3 Stabilization of convective terms
The Navier-Stokes system
∇ · v = 0,
v · ∇v − ν∆v + ∇p = f in Ω,
needs additional stabilization of the convective term if the viscosity is low. One popular
approach consists of adding diffusion in streamline direction (Streamline Upwind PetrovGalerkin (SUPG) method). In combination with the artificial diffusion ansatz for pressure
stabilization, the full stabilized system reads
X
(∇ · v, ξ) +
(αK ∇p, ∇ξ)K = 0,
K
X
(v · ∇v, φ) + ν(∇v, ∇φ) − (p, ·∇φ) +
(δK v · ∇v, v · ∇φ)K = (f, φ).
K
The parameters αK and δK are set to
αK = α0
40
kvkK,∞
6ν
+
2
hK
hK
−1
,
δK = δ0
kvkK,∞
6ν
+
2
hK
hK
−1
,
6.4 LPS Stabilization
where kvkK,∞ is the norm of the local velocity vector. The parameters might be chosen as
1
α0 = δ 0 ≈ .
2
For programming the Matrix it is usually not necessary to implement all derivatives. It
should be sufficient to have:
1
2
void Equation : : Matrix ( . . . )
{
3
4
5
6
7
8
9
10
11
...
// d e r i v a t i v e o f t h e c o n v e c t i o n s t a b i l i z a t i o n terms
A( 1 , 1 ) +=
d e l t a ∗ (U [ 1 ] .m( ) ∗M. x ( ) + U [ 2 ] .m( ) ∗ M. y ( ) )
∗ (U [ 1 ] .m( ) ∗N. x ( ) + U [ 2 ] .m( ) ∗ N. y ( ) ) ;
A( 2 , 2 ) +=
d e l t a ∗ (U [ 1 ] .m( ) ∗M. x ( ) + U [ 2 ] .m( ) ∗ M. y ( ) )
∗ (U [ 1 ] .m( ) ∗N. x ( ) + U [ 2 ] .m( ) ∗ N. y ( ) ) ;
...
}
6.4 LPS Stabilization
Another way to stabilize these equations is the so called Local Projection Stabilization (LPS).
The idea is to suppress local fluctuations by adding the difference between the function itself
and its projection onto a coarser grid.
We use the fluctuation operator
π = id − i2h
where i2h is a projection to the coarser grid. The LPS stabilized formulation for the Stokes
system reads
X
(∇ · v, ξ) +
αK (∇πp, ∇πξ)K = 0
K∈Ωh
and for the Navier-Stokes systems
(∇ · v, ξ) +
X
(αK ∇πp, ∇πξ)K = 0,
K
(v · ∇v, φ) + ν(∇v, ∇φ) − (p, ·∇φ) +
X
(δK v · ∇πv, v · ∇πφ)K = (f, φ).
K
αK and δK can be taken from above.
To implement the projection in Gascoigne, one has to use the class LpsEquation instead of
Equation. It provides the additional functions
41
6 Flow problems and stabilization
1
2
3
4
5
6
7
8
void l p s p o i n t ( double h , const FemFunction& U,
const Vertex2d& v ) const { . . . }
void StabForm ( V e c t o r I t e r a t o r b , const FemFunction& U,
const FemFunction& UP, const T e s t F u n c t i o n& Np)
const { . . . }
void StabMatrix ( EntryMatrix& A, const FemFunction& U,
const T e s t F u n c t i o n& Np , const T e s t F u n c t i o n& Mp)
const { . . . }
which have to be used for the stabilizing terms.
The variables UP, Np and Mp stand for the fluctuation operator π applied to N and M. Gascoigne provides these fluctuation operators if the discretization specified by the keyword
discname in the parameterfile is chosen Q1 or Q2Lps.
Similar as in the case of the SUPG method, it should be sufficient to implement the following
terms in the StabMatrix
1
2
void LpsEquation : : StabMatrix ( . . . )
{
3
4
5
6
7
8
9
10
11
...
// d e r i v a t i v e o f t h e c o n v e c t i o n s t a b i l i z a t i o n terms
A( 1 , 1 ) +=
d e l t a ∗ (U [ 1 ] .m( ) ∗Mp. x ( ) + U [ 2 ] .m( ) ∗ Mp. y ( ) )
∗ (U [ 1 ] .m( ) ∗Np . x ( ) + U [ 2 ] .m( ) ∗ Np . y ( ) ) ;
A( 2 , 2 ) +=
d e l t a ∗ (U [ 1 ] .m( ) ∗Mp. x ( ) + U [ 2 ] .m( ) ∗ Mp. y ( ) )
∗ (U [ 1 ] .m( ) ∗Np . x ( ) + U [ 2 ] .m( ) ∗ Np . y ( ) ) ;
...
}
42
Source Exif Data:
File Type : PDF File Type Extension : pdf MIME Type : application/pdf PDF Version : 1.5 Linearized : Yes Create Date : 2014:03:09 05:39:42+01:00 Creator : TeX Modify Date : 2014:03:10 07:48:06-05:00 PTEX Fullbanner : This is pdfTeX, Version 3.1415926-2.4-1.40.13 (TeX Live 2012/Debian) kpathsea version 6.1.0 XMP Toolkit : Adobe XMP Core 5.2-c001 63.139439, 2010/09/27-13:37:26 Creator Tool : TeX Metadata Date : 2014:03:10 07:48:06-05:00 PTEX Fullbanner : This is pdfTeX, Version 3.1415926-2.4-1.40.13 (TeX Live 2012/Debian) kpathsea version 6.1.0 Producer : pdfTeX-1.40.13 Trapped : False Document ID : uuid:799f9eeb-9dad-4594-8906-b0adab493472 Instance ID : uuid:2bbc756d-d559-405b-9017-0e4545fd460a Format : application/pdf Page Count : 42EXIF Metadata provided by EXIF.tools