Download User Manual Release 1.1

User Manual
Release 1.1
S. Turek, Chr. Becker
University of Heidelberg, Institute for Applied Mathematics
Im Neuenheimer Feld 294, D–69120 Heidelberg, Germany
Heidelberg, 02/01/1998
Contents
Short description
2
1 Mathematical Background
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . .
1.2 Discretization and solution techniques in featflow
1.2.1 Discretization schemes . . . . . . . . . . . . .
1.2.2 Nonlinear solvers . . . . . . . . . . . . . . . .
1.2.3 Linear solvers . . . . . . . . . . . . . . . . . .
1.2.4 Discrete projection method . . . . . . . . . .
1.3 The Navier–Stokes solver packages in featflow . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
5
5
6
10
11
11
12
14
2 Structure of featflow
2.1 The programming structure of featflow . . . . . . . . . . . . . . . . .
2.2 The input data of featflow . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 General comments for user–input . . . . . . . . . . . . . . . . . .
2.2.2 indat2d.f/indat3d.f/parq2d.f/parq3d.f – data files for user
2.2.3 trigen2d.dat – parameter file for trigen2d . . . . . . . . . . .
2.2.4 tr2to3.dat – parameter file for tr2to3 . . . . . . . . . . . . . .
2.2.5 trigen3d.dat – parameter file for trigen3d . . . . . . . . . . .
2.2.6 pp2d.dat/pp3d.dat – parameter files for pp2d and pp3d . . . .
2.2.7 cc2d.dat/cc3d.dat – parameter files for cc2d and cc3d . . . .
2.3 The file structure of featflow . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Subdirectory source – source code for featflow . . . . . . . .
2.3.2 Subdirectory manual – manuals for featflow . . . . . . . . . .
2.3.3 Subdirectory object – system software for featflow . . . . . .
2.3.4 Subdirectory application – applications under featflow . . .
2.3.5 Subdirectory graphic – graphic tools for featflow . . . . . . .
2.3.6 Subdirectory utility – utilities for featflow . . . . . . . . . .
2.4 Installation of featflow . . . . . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
15
15
18
18
19
29
32
33
34
39
43
43
43
43
44
45
45
45
3 Examples for the use of featflow
3.1 The installation example . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 The 2D example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 The 3D example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
48
48
50
60
Bibliography
67
A Appendix: Troubleshooting with featflow
69
B Appendix: The featflow group
71
C Appendix: Future projects in featflow
72
1
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Short description
The progam package featflow is both a user oriented as well as a general purpose
subroutine system for the numerical solution of the incompressible Navier–Stokes equations
in two and three space dimensions. featflow is part of the feast project which has
the aim to develop software which realizes our new mathematical and algorithmic ideas
in combination with high performance computational techniques. For more information
about this project ask the authors or look at the Internet–URL:
http://www.iwr.uni-heidelberg.de/˜featflow
which contains most of our research activities, including the ‘Virtual Album of Fluid
Motion’. featflow is designed for the following three classes of applications:
Education:
Students coming from Mathematics, Physics, Engineering sciences or Computer sciences
learn to handle modern numerical software. This helps for future work in industry or for
doing Diploma or Ph.D. works since the basics of today‘s numerical software engineering, concerning software and hardware, are provided. Actually, we realize this education
program as ‘Software–Praktikum’ at the University of Heidelberg.
Research:
This software (and hence the underlying Mathematics) is running at several universities in
Europe and North–America. As a basic tool it helps to verify own codes by comparisons
with reference configurations, and it can be used for adding new components to solve
more complex problems. For instance, in our group it is the basic solver which is used,
after slight modifications, for exploring new algorithms for compressible media, turbulence
modelling, multi–phase and non–newtonian flows or shape optimization.
Industry:
Since this software is purely based on a mathematically optimized approach it should be
able (if everything is implemented in a correct way) to be much more efficient than most
actual software tools used by industry. Therefore, we try not only to solve ”mathematical”
test problems, but realistic applications, too. We hope to demonstrate, that this approach
is not only faster and more efficient than most other programs, but also applicable for the
same problems. Therefore, one of our main goals is to demonstrate the ”practical” flexibility of our software. A first ”proof” for these statements could be found by performing the
DFG-benchmarks of the project ”Flow Simulation with High Performance Computers”.
2
CONTENTS
3
Mainly active members of the feast-project are the ‘numeric group’ in Heidelberg around
Rannacher and Turek. A list containing (almost) all participants (helping practically and
theoretically) can be found in the Appendix.
featflow 1.1 is based on the FORTRAN77 finite element packages feat2d and feat3d
which are not user oriented systems. They only provide subroutines for several main
steps in a finite element program. The user should be familiar with the mathematical
formulation of the discrete problems. The data structure of feat is transparent so that
modifications or augmentations of the program package are very easy, for instance the
implementation of new elements or the application of error control. For details concerning tools and data structures of these finite element codes the reader is referred to the
manuals [3] and [9]. We consider to use other, more sophisticated finite element tools in
FORTRAN90 under NETSCAPE in the future.
This manual has the following structure: Chapter 1 shows the mathematical background
and explains the temporal and spatial discretizations. Additionally, subjects like nonlinear
solution schemes, multigrid solvers for linear subproblems, discrete projection schemes for
coupled systems, adaptive time step control and other optimization tools are treated explicitly. Inbetween, there is much more literature cited which explains all pointed subjects
more in detail. We describe the discretization and solution process in form of flow charts
since they can be directly translated into the programming structure of featflow.
This structure of programs and subroutines is explained in Chapter 2. We divide all
software into three classes: Preprocessing, solver and postprocessing.
The preprocessing part contains the programs omega2d (description of 2D domains,
graphical generation of coarse tringulations, developed by Matthies/Schieweck in Magdeburg, and which will be soon replaced by our own JAVA–based preprocessing tool), trigen2d (adaptively refined meshes in 2D, output of triangulations), tr2to3 (generation of
3D meshes out of 2D meshes) and trigen3d (analogous as trigen2d but in 3D). These
programs mainly provide (coarse) triangulations for the following solution and graphical
output routines.
The solution part contains the solver packages pp2d/pp3d (as purely time dependent
projection schemes) and cc2d/cc3d (solving stationary and nonstationary problems in
a fully coupled way). The use of both solvers, two- as well as three–dimensional, is
demonstrated and their input parameters are explained. The systems cp2d/cp3d are not
finished yet and will be added to the next release (in fact, you may find in the recent
version a test version of cp2d, together with some other ’new’ software tools!).
The postprocessing part, finally, contains files and shell–scripts for supported graphical packages. These are, in this version, gnuplot (for 1D pictures) and movie.byu or
cquel.byu, resp., avs for 2D- and 3D–graphics. Moreover, it has shown to be easy to
incoorporate output files for grape or some other appropriate graphics, and to add our
own particle tracing tool ptrac.
At the end of this chapter, we show the content of the other subdirectories of featflow
and explain the structure of featflow with its installation, backup, makefile and library
parts.
4
CONTENTS
In addition, most important for the user, sample programs for some standard applications
(installation, preprocessing and the solution of a 2D and 3D problem) are included in
Chapter 3. These can be used as starting programs which may be modified for the actual
application.
There are several main levels to be explained in detail:
–
Installation of featflow, editing of makefiles and other preparation for the use
–
Generation of coarse and refined triangulations for the solution process
–
Input data (data files and input parameter, as for instance boundary conditions and
right hand side)
–
Code execution and explanation of output data
–
Graphical presentation of the data
At the end of this manual, the Appendix contains a discussion of the ”most usual” errors
and problems with featflow and how to survive them, followed by a list with names and
contact addresses of involved persons, and a remark how to get featflow. Additionally,
there is a short section which shows projects being under development for future versions
of featflow
We hope that you enjoy this software and are very grateful for every comment concerning
bugs and critical aspects or subjects to be added.
In fact, this version featflow 1.1 contains only slight modifications compared to the
‘original’ version featflow 1.0! However, we may announce still for 1998 the new version featflow 2.0 which will contain the essentially improved solvers cp2d and some
tools for Boussinesq and/or nonnewtonian flows, including also the improved pre- and
processing tools together with (freely downloadable!!!) AVS EXPRESS modules for high–
end visualization. For a mathematical description of these software releases we recommend
the papers in our paper archive (see the following URL) or our book ‘Efficient solvers for
incompressible flow problems: An algorithmic approach in view of computational aspects’
which will appear by Springer Verlag. Further, we hope to present a first version of our
new feast package.
For more information about these project, take a look at the Internet–URL:
http://www.iwr.uni-heidelberg.de/˜featflow
which contains most of our research activities, including the ‘Virtual Album of Fluid
Motion’.
Heidelberg, August 1998
1. Mathematical Background
1.1. Introduction
The following mathematical description for the implemented software is part of our paper
archive (see the given URL) or of our book:
‘Efficient solvers for incompressible flow problems: An algorithmic approach in view of
computational aspects’
which will appear by Springer Verlag. There, the reader will find much more details.
We consider numerical solution techniques for the nonstationary (or stationary, without
the term ut ) incompressible Navier–Stokes equations,
ut − ν∆u + u · ∇u + ∇p = f
,
∇·u = 0 ,
in Ω × (0, T ] ,
(1.1)
for a given force f and viscosity ν, with any prescribed boundary values on the boundary
∂Ω (of Dirichlet- or Neumann–type, see [10]) and an initial condition at t = 0. Solving
this problem numerically still seems to be a considerable task in the case of long time
calculations and higher Reynolds numbers, particularly in 3D.
The corresponding discretized nonlinear systems of equations may be treated by a coupled
approach in u and p (see [23]) which promises the best stability behaviour, but also
entails the largest numerical effort. Another variant, known as projection method (see [7]),
decouples pressure and velocity, which reduces the problem to the solution of a sequence
of ”simple” (scalar) problems. However, at the same time, it leads to smaller time steps
due to the inherently more explicit character and often suffers from spurious pressure
boundary layers.
These different approaches lead to a large variety of schemes all of which are occuring
in practice since years (see [22] and [21]). Theoretical considerations could provide some
ideas, concerning stability of these schemes, convergence rates for subproblems, necessary
time step sizes, or qualitative behaviour for large Reynolds numbers, but a complete
analysis or quantitative prediction is not possible even today. Therefore, the only way to
make a judgement was to perform numerical tests, at least for some classes of problems
which seem to be representative. This was done in [21].
5
6
Mathematical Background
What we reached is a (finite element) discretization and a solution procedure such that:
1) The finite element spaces for u and p are stable, i.e., satisfy the LBB–condition ([8]).
2) A robust and efficient coupled solver is available.
3) A robust and efficient solver of projection type is available.
4) An efficient nonlinear solution strategy is available.
5) An efficient time step control is available.
The method which seems to satisfy all these requirements consists of discrete projection
schemes with nonconforming linear or rotated multilinear finite elements for u and piecewise constant approximations for p (see [15],[22],[21]). With this approach, we can develop
very efficient solution schemes of both coupled and projection type, with a special nonlinear or linearized treatment of the advection. The resulting solutions are coincident (as
soon as the time steps are small enough), and no spurious pressure oscillations occur. This
approach is the basis of our following theoretical and numerical investigations.
1.2. Discretization and solution techniques in featflow
We first discretize the time derivative in the Navier–Stokes equations (1.1) by one of the
usual time–stepping schemes, with prescribed boundary values for every time step.
Given un and the time step k = tn+1 − tn , then solve for u = un+1 and p = pn+1
u − un
+ θ[−ν∆u + u · ∇u] + ∇p = gn+1
k
,
∇·u = 0
,
in Ω ,
(1.2)
with right hand side
gn+1 := θf n+1 + (1 − θ)f n − (1 − θ)[−ν∆un + un · ∇un ] .
(1.3)
In the past, explicit time–stepping schemes have been commonly used in nonstationary
flow calculations, but because of the severe stability problems inherent in this approach,
the required small time steps prohibited the long time solution of really time–dependent
flows. Due to the high stiffness, one seems to be limited to implicit schemes in the choice
of time–stepping methods for solving this problem. Since implicit methods have become
feasible, thanks to more efficient linear solvers, the schemes most frequently used are either
the simple first–order Backward Euler–scheme (BE), with θ = 1, or the second–order
Crank–Nicolson–scheme (CN), with θ = 1/2. These two methods belong to the group of
one–step–θ–schemes. The CN–scheme occasionally suffers from unexpected instabilities
because of its only weak damping property (not strongly A–stable), while the BE–scheme
is of first order accuracy only. Another method which seems to have the potential to excel
in this competition is the Fractional–step–θ–scheme (FS). It uses three different values for
θ and for the time step k at each time level. For a realistic comparison we define a macro
time step with K = tn+1 − tn as a sequence of 3 time steps of possibly (variable) size k.
Then, in the case of the Backward Euler- or the Crank–Nicolson–scheme, we perform 3
substeps with the same θ as above and time step k = K/3.
1.2 Discretization and solution techniques in featflow
7
√
For the Fractional–step–θ–scheme we proceed as follows. Choosing θ = 1− 22 , θ′ = 1−2θ,
and α = 1−2θ
1−θ , β = 1 − α, the macro time step tn → tn+1 = tn + K is split into three
consecutive substeps (with θ˜ := αθK = βθ′ K):
˜ (un+θ )]un+θ + θK∇pn+θ = [I − βθKN (un )]un + θKf n
[I + θN
∇·un+θ = 0 ,
˜ (un+1−θ )]un+1−θ + θ′ K∇pn+1−θ = [I − αθ′ KN (un+θ )]un+θ + θ′ Kf n+1−θ
[I + θN
∇·un+1−θ = 0 ,
˜ (un+1 )]un+1 + θK∇pn+1 = [I − βθKN (un+1−θ )]un+1−θ + θKf n+1−θ
[I + θN
∇·un+1 = 0 .
Here and in the following, we use the more compact form for the diffusive and advective
part
N (u)u := −ν∆u + u · ∇u .
(1.4)
Being a strongly A–stable scheme, the FS–method possesses the full smoothing property
which is important in the case of rough initial or boundary values. Further, it contains
only very little numerical dissipation which is crucial in the computation of non–enforced
temporal oscillations in the flow. A rigorous theoretical analysis of the FS–scheme (see
[11],[12]) applied to the Navier–Stokes problem establishes second order accuracy for this
special choice of θ. Corresponding numerical tests are performed in [21]. Therefore, this
scheme can combine the advantages of both the classical CN–scheme (2nd order accuracy)
and the BE–scheme (strongly A–stable), but with the same numerical effort.
So, in each time step we have to solve nonlinear problems of the following type:
Given un , parameters k = k(tn+1 ) , θ = θ(tn+1 ) and θi = θi (tn+1 ), i = 1, . . . , 3, then solve
for u = un+1 and p = pn+1
[I + θkN (u)]u + k∇p = [I − θ1 kN (un )]un + θ2 kf n+1 + θ3 kf n
,
∇·u = 0 .
(1.5)
For spatial discretization, we choose a finite element approach. In setting up a finite
element model of the Navier–Stokes equations, one starts with a variational formulation.
On the finite mesh Th (triangles, quadrilaterals or their analogues in 3D) covering the
domain Ω with local element width h, one defines polynomial trial functions for velocity
and pressure. These spaces Hh and Lh should lead to numerically stable approximations,
as h → 0, i.e., they should satisfy the Babuska–Brezzi condition with a mesh–independent
constant γ (see [8]),
(ph , ∇·vh )
min max
≥ γ > 0.
(1.6)
ph ∈Lh vh ∈Hh kph k0 k∇vh k0
Many stable pairs of finite element spaces have been proposed in the literature. Our
favorite candidate is a quadrilateral element which, in 2D, uses piecewise rotated bilinear shape functions for the velocities, spanned by hx2 − y 2 , x, y, 1i, resp., rotated trilinear
shape functions in 3D, spanned by hx2 − y 2 , x2 − z 2 , x, y, z, 1i, and piecewise constant
pressure approximations (see Figure 1.1). The nodal values are the mean values of the
8
Mathematical Background
u,v
p
u,v
u,v
u,v
Figure 1.1: Nodal points of the nonconforming finite element pair in 2D
velocity vector over the element edges / faces (element type 30) or in the corresponding
midpoints (element type 31), and the mean values of the pressure over the elements rendering this approach ”nonconforming”. This element is the natural quadrilateral analogue
of the well–known triangular Stokes element of Crouzeix–Raviart (see [5]). A convergence
analysis of both the parametric (elements E030/E031) and the nonparametric versions
(elements EM30/EM31), which are preferrable on non–tensor product grids or meshes with
large aspect ratios, is given in [15] and very promising computational results are reported
in [19],[23],[21].
This element pair has several important features. It admits, beside the typical Galerkin–
streamline diffusion technique, simple upwind strategies which lead to matrices with certain M–matrix properties. Further, efficient multigrid solvers are available which work
satisfactorily over the whole range of relevant Reynolds numbers, 1 ≤ Re ≤ 105 , and also
on nonuniform meshes. In [22], we have even shown by a complexity analysis that this
pair of elements is most efficient compared to other finite element pairs, especially in the
case of highly nonstationary flows. In combination with the discrete projection methods,
see [22], it works very robust and efficient in a multigrid code also on highly stretched and
anisotropic grids.
Using the same symbols u and p also for the coefficient vectors in the nodal representation
for the functions u and p, the discrete version of problem (1.5) may be written as a
(nonlinear) algebraic system of the form:
Given un , a right hand side g and a time step k, then solve for u = un+1 and p = pn+1
Su + kBp = g
,
BT u = 0 ,
(1.7)
with matrix S and right hand side g such that
Su = [M + θkN (u)]u
,
g = [M − θ1 kN (un )]un + θ2 kf n+1 + θ3 kf n .
(1.8)
Here, M is the mass matrix and N (.) the advection matrix containing the diffusive and
convective parts corresponding to the nonlinear form in (1.4). For dominant transport
the advection part may include some stabilization, for instance, some upwind mechanism
(see [23]) or the streamline diffusion method ([26]). B is the gradient matrix, and −B T
1.2 Discretization and solution techniques in featflow
9
the transposed divergence matrix. With Ml we denote the lumped mass matrix which is
diagonal.
Two possible approaches for solving these discrete nonlinear problems are:
1) We first treat the nonlinearity by an outer nonlinear iteration of fixed point- or quasi–
Newton type or by a linearization technique through extrapolation in time, and we obtain
linear indefinite subproblems of Oseen type which can be solved by a coupled (versions
cc2d/ cc3d) or a splitting approach (versions cp2d/ cp3d).
2) We first split the coupled problem and obtain definite problems in u (Burgers–equations)
as well as in p (linear pressure–Poisson-, resp., quasi–Poisson problems). Then we treat
the nonlinear problems in u by an appropriate nonlinear iteration or a linearization technique (versions pp2d/ pp3d).
In our applications the nonlinear problems are solved by the adaptive fixed point defect
correction method (see [23]), while for the solution of the coupled problems, resp., for their
decoupling, the discrete projection method formalism is used (see [22]). To understand
completely the algorithm, we will shortly present other important tools of the solver,
namely the multilevel solvers for the corresponding linear subproblems (Oseen problems,
(convection–diffusion problems for u, quasi–Poisson problems for p), the adaptive time
step control and defect optimization procedures. In compact form the components of the
solvers are explained in the following subsections.
Discretization schemes
Linear solvers
FEATFLOW
Navier-Stokes
Discrete Proj. Scheme
tree
Nonlinear solvers
Figure 1.2: The solver structure of featflow
10
Mathematical Background
1.2.1. Discretization schemes
Spatial discretization
–
nonconforming nonparametric rotated bi / trilinear finite elements
with edge oriented d.o.f.’s for velocity
–
piecewise constant finite elements for pressure
–
a–priori adapted coarse meshes with exact boundary adaptation
during successive refinements
–
adaptive refinement leading to conforming triangulations possible
–
adaptive FE-upwinding for convective terms (Samarskij upwinding)
–
adaptive streamline–diffusion for convective terms
–
natural do nothing b.c.’s
–
flux- and pressure drop b.c.’s
Temporal discretization
–
One–step–θ–scheme (Implicit Euler, Crank–Nicolson scheme)
–
Fractional–step–θ–scheme as strongly A–stable implicit time stepping scheme of 2nd order
–
adaptive time stepping by estimating the local discretization error
(3 substeps with ∆t and 1 substep with 3∆t)
–
extrapolation of partial solutions for higher accuracy
1.2 Discretization and solution techniques in featflow
1.2.2. Nonlinear solvers
Nonlinear iteration
–
for stationary and nonstationary problems of Navier–Stokes- or
Burgers type
–
fixed point iteration (quasi–Newton) with adaptive step length control by nonlinear defect minimization
–
additional defect correction possible
Linearization
–
for nonstationary problems only
–
semi–implicit treatment of nonlinear convective terms by linear
extrapolation (in time of 2nd order)
1.2.3. Linear solvers
Multigrid for velocity and pressure simultanously (Oseen):
–
nonconforming mesh adapted makro-elementwise interpolation for
grid–transfer
–
adaptive step length control for correction step (with F–cycle)
–
Vanca–like block–Gauß-Seidel scheme as smoother and solver
Multigrid for velocity (Burgers) and pressure (quasi–Poisson):
–
nonconforming mesh adapted makro-elementwise interpolation for
grid–transfer
–
adaptive step length control for correction step (with F–cycle)
–
ILU/SOR–scheme as smoother with special renumbering strategies
11
12
Mathematical Background
1.2.4. Discrete projection method
For linear (or nonlinear) coupled problems:
Su + kBp = g
,
BT u = 0
Given p0 , solve for l = 1, . . . , L:
1
pl = pl−1 − αl [B T C −1 B]−1 (B T S −1 Bpl−1 − B T S −1 g)
k
Set p := pL , and determine u through
Su = g − kBp +
(⇒ B T u = 0)
k l
[α I − SC −1 ]B(pL − pL−1 )
αl
Classical Richardson–scheme for Schur–complement formulation with
preconditioner P −1 = [B T C −1 B]−1
C=S
C = Ml
,
L=1
for cc2d/ cc3d
, L≥1
for pp2d/ pp3d
If C = Ml ⇒ P corresponds to discrete Laplacian with minimal matrix
stencil (5 in 2D/7 in 3D) independent of the mesh
Further improvements:
–
Additional preconditioner possible (of diffusive type)
–
Elimination of pressure boundary layers
–
Multilevel Discrete Projection scheme (cp2d/ cp3d)
−→
−→
”fastest” solver/discretization for fully nonstationary problems
”robust” solver/discretization for quasi–stationary problems
All proposed discretization and solution steps can be represented by the following flow
chart (”the Navier-Stokes tree”), leading to the methods mentioned above. It is remarkable
that most of all existing Navier-Stokes solvers can be interpreted by this tree structure.
This is a very powerfull tool for the understanding and, then, optimization of existing
algorithms and software.
1.2 Discretization and solution techniques in featflow
NSE :
13
u t + N(u)u + ∇p = f , ∇⋅ u = 0 in 2D/3D + B.C .
BE, CN, FR
[ Ι + θkN(u)]u + k∇p = g , ∇⋅ u = 0
FEM, FV, FD
[M + θkN(u h)]u h + kΒph = g h , ΒT u = 0
Outer: N nonlinear steps
Inner: L DPM steps (Oseen)
C=S
L=1
C=Ml
L>=1
S(u h)u h + kΒph = g h
ΒT uh = 0
Outer: L DPM steps (C=Ml )
Inner: N nonlinear steps (Burgers)
L>1
L=1
Chorin (p 0 = 0 )
CC2D/CC3D
CP2D/CP3D
PP2D/PP3D
Van Kan (p 0 = p old )
N=1
"Galerkin schemes"
N>1
"Projection schemes"
Structure of (almost) all solvers
Figure 1.3: The Navier-Stokes tree
14
Mathematical Background
1.3. The Navier–Stokes solver packages in featflow
The different solution schemes for the Navier–Stokes equations we can apply are the following ones:
Coupled solver cc2d/ cc3d:
We use the adaptive fixed point defect correction method as outer iteration, while for the
linear coupled subproblems we choose C = S and α = 1. Hence, the linear equations in
the nonlinear process are solved in ”one iteration step”, meaning that L = 1. This solver
is a multilevel based approach (see [19]) with a block Gauß-Seidel–scheme as smoothing
operation (Vanca–smoother). There may be problems on very anisotropic meshes, and
there is no gain in efficiency for fully nonstationary problems if ∆t → 0. This solver is,
actually, preferrable for the case of low Reynolds numbers (stationary or quasi–stationary
problems).
Mixed solver cp2d/ cp3d:
We use again the adaptive fixed point defect correction method as outer iteration and
select the preconditioner C = Ml . We will obtain the same solutions as by cc2d/ cc3d,
however in a more robust way. First tests show that this solver works in a very efficient
way, independent of the mesh and Re number, being very efficient for fully nonstationary
problems. The essential tool is a so called ”multilevel discrete projection” algorithm for the
linear coupled subproblems. This solver will be added to featflow in the next release.
Projection solver pp2d/ pp3d:
First, we apply a decoupling step for u and p as outer iteration (”nonlinear discrete
projection scheme”). We choose again C = Ml , but perform only L = 1 iteration in each
time step. In the fully nonlinear case, we use the same fixed point iteration, as explained
above, for the transport–diffusion–step with nonlinear operator S (”Burgers equations”).
In an analogous way, we treat (by extrapolation) the linearized schemes.
All versions of cc2d/ cc3d and cp2d/ cp3d (and even pp2d/ pp3d for L large enough)
lead to the ”same” solutions if the numbers of nonlinear steps N and discrete projection
steps L are large enough. The complexity analysis in [22] shows that in 2D one single
iteration of the coupled scheme (C = S) can be expected to cost at least 10 (2D) until
20 (3D) times more than one iteration of the operator splitting method (C = Ml ). At
the same time, we have to use more nonlinear sweeps or smaller time steps for C = Ml ,
due to the more explicit character of the scheme. However, extensive numerical tests in
[21] show that the projection solvers are superior in most cases compared to cc2d/ cc3d,
especially for highly nonstationary flows.
In the next chapter one will see that featflow follows exactly this flow chart (”the
Navier-Stokes tree”) and that the programming structure is a direct translation of this
tree structure.
2. Structure of featflow
2.1. The programming structure of featflow
As pointed in the last section, the programming structure of featflow follows directly
the tree structure in diagram 1.3.
The corresponding programs cc2d/ cc3d and pp2d/ pp3d (cp2d/ cp3d are not yet
finished) are the main programs which form the body of the code. Their task is to
initialize all data, to build triangulations, linear matrices, pointer structures and right
hand sides on all grid levels (all in the file init1.f) and to generate all vectors needed
(also in init1.f). After selecting all parameters for the chosen time stepping scheme
and corresponding adaptive time step size, the discretization process is finished and the
time stepping loop may begin. That means, everything is prepared for solving the pointed
nonlinear coupled generalized stationary Navier–Stokes problems corresponding to each
time level. After doing this in the subroutine prostp.f, resp., mgstp.f, the main program
cc2d.f (or cc3d.f, pp2d.f, pp3d.f) continues with the time step control, selecting a
new time step size for the next (macro) step or repeating the last one if necessary, then
evaluates the actual solution vector with some extrapolation in time if needed, and enters
the output–subroutines fpost.f or error.f. These main programs provide most of the
output for monitoring all actions and give the needed workspace and computer times for
the performed application.
The subroutine on the next level is prostp.f, resp., mgstp.f, solving (or approximating
only in the case of the pure projection schemes pp2d/ pp3d) a number of generalized
stationary Navier–Stokes equations corresponding to the actual macro time level. The
number of these problems and their parameters depend on the scheme and the time step
control chosen. Furthermore, the corresponding boundary conditions are set in this file
(by bndry.f).
As shown in the previous tree diagram, at this point there are the differences between
the codes pp2d/ pp3d and cc2d/ cc3d. While the fully coupled versions cc2d/ cc3d
perform an outer linearization by applying first the adaptive fixed point control as nonlinear solver, pp2d/ pp3d first decouple velocity and pressure, and then solve a nonlinear
Burgers–equation and a linear quasi–Poisson problem.
In detail, pp2d/ pp3d do the following for ”solving” one nonlinear coupled equation: With
given pressure as right hand side, a nonlinear Burgers–equation is solved in nsdef.f, with
adaptive defect optimization (optcnl.f) and multigrid solvers for the linear convection–diffusion problems (m011.f). With the divergence of this intermediate velocity as right
hand side a corresponding update problem of quasi–Poisson type for the pressure is solved,
15
16
Structure of featflow
by a multigrid scheme (again m011.f), and the auxiliary solution is added to the old
pressure to obtain the new one. Finally, the diffusive preconditioner may be additionally
obtained to improve the new pressure. This approach corresponds to the special version
L = 1 of the nonlinear discrete projection scheme. We perform this decoupling process
only once, and, hence, we obtain an approximate solution in each time step only. However,
that is very similar to the idea of the (classical) projection schemes (of Van Kan [27], resp.,
Chorin [4]).
In contrast, cc2d/ cc3d (the upcoming cp2d/ cp3d analogously) perform one nonlinear
coupled solution step in a different way: The file nsdef.f handles the solution process
being a solver of the nonlinear coupled equation, with performing Oseen equations (linearized Navier–Stokes equations) as preconditioner in each nonlinear step (m011.f, with
a special block Gauss-Seidel smoother), and following defect minimization (optcnl.f).
The versions cp2d/ cp3d will differ in solving the linearized coupled equations by a special multilevel scheme involving the linear discrete projection method (”multilevel discrete
projection method”).
This approach guarantees the fully coupled solution and, hence, a larger stability and improved accuracy compared to pp2d/ pp3d. However, the numerical effort increases and,
in this actual release, cc2d/ cc3d is preferrable in the case of low Reynolds numbers and
on moderate meshes only.
Following these remarks the programming structure of cc2d/ cc3d and pp2d/ pp3d can
be represented by the following diagrams.
cc2d.f/cc3d.f
mgstp.f
init1.f
fpost.f
error.f
gdat.f
xmrout.f
bbuild.f
nsdef.f
dfkt.f
gupwd.f
supwg.f
bndry.f
m011.f
optcnl.f
vanca.f
mgrout.f
Figure 2.1: The cc2d/ cc3d structure
2.1 The programming structure of featflow
17
pp2d.f/pp3d.f
prostp.f
init1.f
fpost.f
error.f
gdat.f
xmrout.f
bbuild.f
projma.f
nsdef.f
dfkt.f
gupwd.f
supwg.f
bndry.f
m011.f
optcnl.f
m011.f
dfkt.f
mgrout.f
mgrout.f
Figure 2.2: The pp2d/ pp3d structure
Beside the explained solver tools there are some other programs belonging to featflow
which are essential for the preprocessing (description of domain Ω, coarsest triangulation, pre–adapted meshes, output for graphics) and for the preparation of an application
(generation and interpolation of a start solution). These are in detail:
omega2d:
omega2d is a graphical preprocessor which provides tools for describing and defining the
boundary components of a 2D domain Ω and for generating a first (and very coarse) triangulation of this domain. This process is fully interactive and mouse–oriented, however,
actually, running under DOS and WINDOWS 3.1 only. It has been created by Matthies/
Schieweck in Magdeburg, and further development is under progress. In the next chapter
we demonstrate how to use it, and the complete manual of Version 1.6 can be found in
the manual directory of featflow.
trigen2d:
trigen2d is a preprocessing tool for designing ”better” 2D coarse triangulations and
to write the corresponding data in some special formats onto hard disc. It reads in a
parametrization (in standard feat- or in omega2d–format) and a corresponding mesh
(in feat–format), and refines it locally in a successive manner. The user has different
possibilities for mesh refinement, and in this release the elements to be refined have to
be given in a list. We hope to provide a graphical interface for this marking process in
the next version, and, additionally, to use this tool in a fully adaptive mesh refinement
18
Structure of featflow
process. Furthermore, a boundary check on consistency for complex domains is included.
The formats for data–output are: formatted or unformatted feat–style, byu–style and
avs-ucd–format.
tr2to3:
tr2to3 is a tool providing 3D mesh generation from 2D meshes, all in feat–format.
This tool is designed for applications when the 3D mesh is simply constructed by popping
2D mesh layers in a sandwich–like technique. For instance, problems in ducts around
cylinders or other ”prism”-like bodies are belonging to this class (see the example in the
next chapter).
trigen3d:
trigen3d is doing the same job as trigen2d, but in 3 dimensions. However, up to now,
only feat–format is supported (that means, the boundary of Ω must be defined by the
standard feat parametrization file). We hope to add some appropriate CAD–tools in
future, for being able to model much more complex domains. Additionally, the adaptive
grid refinement procedure is not finished yet, and will be part of the next version. However,
at least the data–output is performing analogously as in 2D.
intpol2d and intpol3d:
These programs will be tools for interpolating a given solution vector on a given mesh to
an arbitrary different triangulation, in 2D as well as in 3D. Both data and grid have to be
given in feat–format. These tools will be available in the next release.
2.2. The input data of featflow
In this section, we explain all input data which are needed to start an application with
featflow. Let’s start with the preprocessing tools, and with some general comments. In
the next chapter, we demonstrate explicitely the meaning of all input data.
2.2.1. General comments for user–input
First of all, files containing pure input parameter for user applications have to be in the
directory #data. The file name is the name of the corresponding program, ended with
.dat, for instance #data/pp2d.dat. These files will be explained in detail in the following
subsections. Examples for different ”styles” of parameter parameter files, concerning robustness, efficiency and small stoarge amount, can be found in the directory application
/data example.
Furthermore, the only source–files which have to be modified (by editing or copying) are
the parametrization files parq2d.f, resp., parq3d.f, and the data files indat2d.f, resp.,
2.2 The input data of featflow
19
indat3d.f. They contain all information about boundary values, right hand sides, exact solutions and their derivatives and output constants like lift and drag or integral or
pointwise values in some special coordinates. The last file which has to be edited is a
corresponding .inc–file, for instance pp2d.inc, containing the needed storage size for the
applied program (for instance, NNWORK=1.000.000 means 1 million double precision elements needed, or, resp., a storage amount of 8 Mbyte). These files are usually located in the
directory input files. The best way is to collect all data files together in a subdirectory,
corresponding to the actual application. This is usually done by featflow, for instance
in the directories application/user start, application/example and application/
comp. Furthermore, the corresponding makefiles are also located in input files, such
that all data files are concentrated at one place.
2.2.2. indat2d.f/indat3d.f/parq2d.f/parq3d.f – data files for user
indat2d.f:
************************************************************************
DOUBLE PRECISION FUNCTION FDATIN(ITYP,IBLOC,X,Y,TIMENS,RE)
*
*
Prescribed data for files coeff.f and bndry.f
************************************************************************
IMPLICIT DOUBLE PRECISION(A,C-H,O-U,W-Z),LOGICAL(B)
PARAMETER (PI=3.1415926535897931D0)
C
FDATIN=0D0
C
C=======================================================================
C *** Case 1: Velocity boundary values and/or exact solution
C=======================================================================
C
IF (ITYP.EQ.1) THEN
C
IF (IBLOC.EQ.1) THEN
IF (X.EQ.0D0) FDATIN=4D0*0.3D0/0.1681D0*Y*(0.41D0-Y)
*
*SIN(0.5D0*PI*MIN(TIMENS,1D0)/1D0)
ENDIF
C
IF (IBLOC.EQ.2) THEN
IF (X.EQ.0.0D0) FDATIN=0D0
ENDIF
C
ENDIF
C
C=======================================================================
C *** Case 2: Velocity x-derivative of exact solution
C=======================================================================
C
IF (ITYP.EQ.2) THEN
C
IF (IBLOC.EQ.1) THEN
IF (X.EQ.0.0D0) FDATIN=0D0
ENDIF
C
IF (IBLOC.EQ.2) THEN
IF (X.EQ.0.0D0) FDATIN=0D0
ENDIF
C
ENDIF
C
C=======================================================================
C *** Case 3: Velocity y-derivative of exact solution
20
Structure of featflow
C=======================================================================
C
IF (ITYP.EQ.3) THEN
C
IF (IBLOC.EQ.1) THEN
FDATIN=0D0
ENDIF
C
IF (IBLOC.EQ.2) THEN
FDATIN=0D0
ENDIF
C
ENDIF
C
C=======================================================================
C *** Case 4: Exact pressure solution
C=======================================================================
C
IF (ITYP.EQ.4) THEN
C
FDATIN=0D0
C
ENDIF
C
C=======================================================================
C *** Case 5: Right hand side for momentum equation
C=======================================================================
C
IF (ITYP.EQ.5) THEN
C
IF (IBLOC.EQ.1) THEN
FDATIN=0D0
ENDIF
C
IF (IBLOC.EQ.2) THEN
FDATIN=0D0
ENDIF
C
ENDIF
C
C=======================================================================
C *** Case 6: Right hand side for continuity equation
C=======================================================================
C
IF (ITYP.EQ.6) THEN
C
FDATIN=0D0
C
ENDIF
C
C=======================================================================
C *** Case 7: Mean pressure values
C=======================================================================
C
IF (ITYP.EQ.7) THEN
DPAR=X
INPR=IBLOC
C
IF ((DPAR.GT.1D0).AND.(DPAR.LT.2D0).AND.(INPR.EQ.1)) THEN
FDATIN=0D0
ENDIF
C
IF ((DPAR.GT.3D0).AND.(DPAR.LT.4D0).AND.(INPR.EQ.1)) THEN
FDATIN=0D0
ENDIF
C
ENDIF
C
99999 END
C
2.2 The input data of featflow
************************************************************************
SUBROUTINE NEUDAT(INPART,INPRN,DPARN1,DPARN2,TIMENS)
*
*
Neumann-boundary part
************************************************************************
IMPLICIT DOUBLE PRECISION(A,C-H,O-U,W-Z),LOGICAL(B)
C
C=======================================================================
C *** Case 0: Set number of Neumann-boundary parts
C=======================================================================
C
IF (INPART.EQ.0) THEN
INPART=1
ENDIF
C
C=======================================================================
C *** Case <>0: Specify Neumann-boundary parts
C=======================================================================
C
IF (INPART.GT.0) THEN
C
IF (INPART.EQ.1) THEN
DPARN1=1D0
DPARN2=2D0
INPRN =1
ENDIF
C
ENDIF
C
99999 END
C
************************************************************************
SUBROUTINE PTSDAT(TIMENS,DNU)
*
*
Data for Point-output (for fpost and bdpres and bdforc)
************************************************************************
IMPLICIT DOUBLE PRECISION(A,C-H,O-U,W-Z),LOGICAL(B)
COMMON /NSPTS/ KPU(2),KPP(4),KPX(4),KPI(2),DPI(2,2),DPF(2)
SAVE
C
EXTERNAL UE
C
C=======================================================================
C *** Points for velocity, pressure and flux-tracing
C=======================================================================
C
KPU(1)=22
KPU(2)=2
C
KPP(1)=42
KPP(2)=46
KPP(3)=2
KPP(4)=22
C
KPX(1)=41
KPX(2)=10
C
C=======================================================================
C *** Parameters for 2 integral pressures in bdpres.f
C=======================================================================
C
KPI(1) =2
DPI(1,1)=0D0
DPI(2,1)=1D0
C
KPI(2) =2
DPI(1,2)=0.00D0
DPI(2,2)=0.25D0
C
C=======================================================================
21
22
Structure of featflow
C *** Parameters for lift (DFW) and drag (DAW) in bdforc.f (INPR=2)
C ***
C *** dfw=2 int_s [dpf(1) dut/dn n_y - p n_x] ds / dpf(2)
C *** daw=2 int_s [dpf(1) dut/dn n_x + p n_y] ds / dpf(2)
C ***
C=======================================================================
C
RHO =1.0D0
DIST =0.1D0
UMEAN=0.2D0
C
DPF(1)=RHO*DNU
DPF(2)=RHO*DIST*UMEAN**2
C
99999 END
Most command lines have not to be explained more in detail (even being written in
FORTRAN77 they explain themselves); only a few comments are necessary:
In the case IF (ITYP.EQ.7) THEN ... in FDATIN the mean pressure values on boundary
parts corresponding to natural b.c.’s can be described. This is done by prescribing the
starting parameter value (DPAR1) and the ending parameter value (DPAR2) corresponding to
the parametrization which is used for describing the geometrical boundary. Additionally,
the number of the boundary component has to be defined (INPR).
In NEUDAT the boundary parts containing natural boundary conditions are defined. This is
done by prescribing the number of such boundary parts (by setting INPART), and then again
by setting the starting parameter value (DPARN1), the ending parameter value (DPARN2)
and the corresponding boundary component (INPRN).
In PTSDAT certain mesh points are defined in which some velocity components (KPU(1),
KPU(2)), some pressure values (KPP(1),KPP(2),KPP(3),KPP(4)) and a flux value (defined
as difference of streamfunction values KPX(1) and KPX(2)) are printed for the runtime
protocol. Furthermore, 2 integral pressure values on some fixed boundary parts (parameter
values DPI(1,1),DPI(2,1) and boundary component KPI(1), resp., DPI(1,2),DPI(2,2)
and KPI(2)) are defined appropriately. And finally, some constants for the calculation of
lift and drag on boundary component 2 (!) are prescribed. If not desired, set DPF(1),
DPF(2) to 0, or KPI(1), KPI(2) to 0.
2.2 The input data of featflow
indat3d.f:
************************************************************************
DOUBLE PRECISION FUNCTION FDATIN(ITYP,IBLOC,X,Y,Z,TIMENS,RE)
*
*
Prescribed data for files coeff.f and bndry.f
************************************************************************
IMPLICIT DOUBLE PRECISION(A,C-H,O-U,W-Z),LOGICAL(B)
PARAMETER (PI=3.1415926535897931D0)
C
FDATIN=0D0
C
C=======================================================================
C *** Case 1: Velocity boundary values and/or exact solution
C=======================================================================
C
IF (ITYP.EQ.1) THEN
C
IF (IBLOC.EQ.1) THEN
IF (X.EQ.0D0)
*
FDATIN=16D0*0.45D0/(0.41D0**4)*Y*(0.41D0-Y)*Z*(0.41D0-Z)
ENDIF
C
IF (IBLOC.EQ.2) THEN
IF (X.EQ.0.0D0) FDATIN=0D0
ENDIF
C
IF (IBLOC.EQ.3) THEN
IF (X.EQ.0.0D0) FDATIN=0D0
ENDIF
C
ENDIF
C
C=======================================================================
C *** Case 2: Velocity x-derivative of exact solution
C=======================================================================
C
IF (ITYP.EQ.2) THEN
C
IF (IBLOC.EQ.1) THEN
IF (X.EQ.0.0D0) FDATIN=0D0
ENDIF
C
IF (IBLOC.EQ.2) THEN
IF (X.EQ.0.0D0) FDATIN=0D0
ENDIF
C
IF (IBLOC.EQ.3) THEN
IF (X.EQ.0.0D0) FDATIN=0D0
ENDIF
ENDIF
C
C=======================================================================
C *** Case 3: Velocity y-derivative of exact solution
C=======================================================================
C
IF (ITYP.EQ.3) THEN
C
IF (IBLOC.EQ.1) THEN
FDATIN=0D0
ENDIF
C
IF (IBLOC.EQ.2) THEN
FDATIN=0D0
ENDIF
C
IF (IBLOC.EQ.3) THEN
FDATIN=0D0
ENDIF
23
24
Structure of featflow
C
ENDIF
C
C=======================================================================
C *** Case 4: Velocity z-derivative of exact solution
C=======================================================================
C
IF (ITYP.EQ.4) THEN
C
IF (IBLOC.EQ.1) THEN
FDATIN=0D0
ENDIF
C
IF (IBLOC.EQ.2) THEN
FDATIN=0D0
ENDIF
C
IF (IBLOC.EQ.3) THEN
FDATIN=0D0
ENDIF
C
ENDIF
C
C=======================================================================
C *** Case 5: Exact pressure solution
C=======================================================================
C
IF (ITYP.EQ.5) THEN
C
FDATIN=0D0
C
ENDIF
C
C=======================================================================
C *** Case 6: Right hand side for momentum equation
C=======================================================================
C
IF (ITYP.EQ.6) THEN
C
IF (IBLOC.EQ.1) THEN
FDATIN=0D0
ENDIF
C
IF (IBLOC.EQ.2) THEN
FDATIN=0D0
ENDIF
C
IF (IBLOC.EQ.3) THEN
FDATIN=0D0
ENDIF
C
ENDIF
C
C=======================================================================
C *** Case 7: Right hand side for continuity equation
C=======================================================================
C
IF (ITYP.EQ.7) THEN
C
FDATIN=0D0
C
ENDIF
C
C=======================================================================
C *** Case 8: Mean pressure values
C=======================================================================
C
IF (ITYP.EQ.8) THEN
C
FDATIN=0D0
2.2 The input data of featflow
C
ENDIF
C
99999 END
C
************************************************************************
SUBROUTINE NEUDAT(IEL,INPR,PX,PY,PZ,TIMENS,IFLAG)
*
*
Neumann-boundary part
************************************************************************
IMPLICIT DOUBLE PRECISION(A,C-H,O-U,W-Z),LOGICAL(B)
C
C=======================================================================
C *** Set Neumann-boundary parts
C=======================================================================
C
IF (PX.EQ.2.5D0) THEN
IFLAG=1
ENDIF
C
99999 END
C
************************************************************************
SUBROUTINE BDPDAT(IEL,INPR,PX,PY,PZ,TIMENS,IFLAG1,IFLAG2)
*
*
Pressure integral boundary part
************************************************************************
IMPLICIT DOUBLE PRECISION(A,C-H,O-U,W-Z),LOGICAL(B)
C
C=======================================================================
C *** Set pressure integral boundary parts
C=======================================================================
C
IF (((PZ.GT.0.00D0).AND.(PZ.LT.0.41D0)).AND.
*
((PX.GE.0.45D0).AND.(PX.LE.0.55D0)).AND.
*
((PY.GE.0.15D0).AND.(PY.LE.0.25D0))) THEN
IFLAG1=1
ENDIF
C
IF (((PZ.GT.0.00D0).AND.(PZ.LT.0.41D0)).AND.
*
(PX.EQ.0.45D0)
.AND.
*
((PY.GE.0.15D0).AND.(PY.LE.0.25D0))) THEN
IFLAG2=1
ENDIF
C
99999 END
C
************************************************************************
SUBROUTINE BDFDAT(IEL,INPR,PX,PY,PZ,TIMENS,DNU,IFLAG,DPF1,DPF2)
*
*
lift and drag data
************************************************************************
IMPLICIT DOUBLE PRECISION(A,C-H,O-U,W-Z),LOGICAL(B)
C
C=======================================================================
C *** Parameters for lift (DFW) and drag (DAW) in bdforc.f (INPR=2)
C ***
C *** dfw=2 int_s [dpf(1) dut/dn n_y - p n_x] ds / dpf(2)
C *** daw=2 int_s [dpf(1) dut/dn n_x + p n_y] ds / dpf(2)
C ***
C=======================================================================
C
RHO =1.0D0
DIST =0.041D0
UMEAN=0.2D0
C
DPF1=RHO*DNU
DPF2=RHO*DIST*UMEAN**2
C
C=======================================================================
25
26
Structure of featflow
C
*
*
IF (((PZ.GT.0.00D0).AND.(PZ.LT.0.41D0)).AND.
((PX.GE.0.45D0).AND.(PX.LE.0.55D0)).AND.
((PY.GE.0.15D0).AND.(PY.LE.0.25D0))) THEN
IFLAG=1
ENDIF
C
99999 END
C
************************************************************************
SUBROUTINE PTSDAT(TIMENS,DNU)
*
*
Data for Point-output (for fpost)
************************************************************************
IMPLICIT DOUBLE PRECISION(A,C-H,O-U,W-Z),LOGICAL(B)
COMMON /NSPTS/ KPU(2),KPP(4)
SAVE
C
C=======================================================================
C *** Points for velocity and pressure
C=======================================================================
C
KPU(1)=7263
KPU(2)=3444
C
KPP(1)=3491
KPP(2)=3557
KPP(3)=3575
KPP(4)=3444
C
99999 END
In NEUDAT the boundary parts containing natural boundary conditions are defined. This is
done by prescribing the cartesian coordinates for this manifold. Analogously, in BDPDAT the
coordinates for the pressure integral boundary parts are set, and in BDFDAT the coordinates
for calculating lift and drag. If not desired, set the corresponding IFLAG parameters to 0.
Again, in PTSDAT certain mesh points are defined in which some velocity components
(KPU(1), KPU(2)) and some pressure values (KPP(1),KPP(2),KPP(3),KPP(4)) are printed
for the runtime protocol.
2.2 The input data of featflow
27
parq2d.f:
This file is either a standard feat parametrization file which is created by your own (see
the feat2d manual) or the special omega2d parametrization file. In this case, the file
parpre.f has to be copied onto parq2d.f.
parq3d.f:
************************************************************************
DOUBLE PRECISION FUNCTION PARX(T1,T2,T3,IBCT)
IMPLICIT REAL*8 (A-H,O-Z)
PARX=T1
99999 END
C
DOUBLE PRECISION FUNCTION PARY(T1,T2,T3,IBCT)
IMPLICIT REAL*8 (A-H,O-Z)
PARY=T2
99999 END
C
DOUBLE PRECISION FUNCTION PARZ(T1,T2,T3,IBCT)
IMPLICIT REAL*8 (A-H,O-Z)
PARZ=T3
GOTO 99999
99999 END
C
************************************************************************
SUBROUTINE
TRPARV (DCORVG,KNPR,KVEL,NVT,NVEL)
IMPLICIT DOUBLE PRECISION (A,C-H,O-U,W-Z),LOGICAL(B)
DIMENSION DCORVG(3,*),KNPR(*),KVEL(NVEL,*)
C
C----------------------------------------------------------------------C
DO 10 IVT=1,NVT
INPR=KNPR(IVT)
IF (INPR.EQ.0) GOTO 10
C
IEL=KVEL(4,IVT)
PX=DCORVG(1,IVT)
PY=DCORVG(2,IVT)
PZ=DCORVG(3,IVT)
C
IF ((ABS(PZ-0.0D0).GT.1D-8).AND.(ABS(PZ-0.41D0).GT.1D-8).AND.
*
(ABS(PX-0.0D0).GT.1D-8).AND.(ABS(PX-2.50D0).GT.1D-8).AND.
*
(ABS(PY-0.0D0).GT.1D-8).AND.(ABS(PY-0.41D0).GT.1D-8)) THEN
PXM=0.50D0
PYM=0.20D0
RAD=0.05D0
DL=SQRT((PX-PXM)**2+(PY-PYM)**2)
DCORVG(1,IVT)=PXM+RAD/DL*(PX-PXM)
DCORVG(2,IVT)=PYM+RAD/DL*(PY-PYM)
GOTO 10
ENDIF
C
C
IF ((ABS(PZ-0.0D0).LE.1D-8).OR.(ABS(PZ-0.41D0).LE.1D-8)) THEN
PXM=0.50D0
PYM=0.20D0
RAD=0.05D0
RADH=0.05D0+1D-8
IF ((ABS(PX-PXM).LE.RADH).AND.(ABS(PY-PYM).LE.RADH).AND.
*
(IEL.EQ.0)) THEN
DL=SQRT((PX-PXM)**2+(PY-PYM)**2)
DCORVG(1,IVT)=PXM+RAD/DL*(PX-PXM)
DCORVG(2,IVT)=PYM+RAD/DL*(PY-PYM)
GOTO 10
ENDIF
ENDIF
28
C
10
C
Structure of featflow
CONTINUE
END
The functions parx, pary, parz simply prescribe the cartesian coordinates, while the
subroutine trparv provides transformations to more complex domains. In this case, a
channel with an innner cylinder around (0.5, 0.2) and radius r = 0.05 is prescribed. This
file corresponds to the example in chapter 3.
2.2 The input data of featflow
29
2.2.3. trigen2d.dat – parameter file for trigen2d
M
MT
ICHECK
IMESH
CPARM
IBDCHK
IBYU
IAVS
NLEV
IFMT
CFILEI
CFILEO
ITYPEL
feat parameter for output
= 0: no output
> 0: output, see feat2d manual
feat parameter for terminal output
= 0: no output
> 0: output, see feat2d manual
feat parameter for subroutine tracing
= 0: no tracing
> 0: tracing, see feat2d manual
parameter for type of parametrization
= 0: feat parametrization
= 1: omega2d parametrization
name of omega2d parametrization file
parameter for boundary checking for further refinements
= 0: no checking
> 0: = number of finer boundary points for check of boundary consistency for further refinement
level for byu output
level for avs output
number of refined levels-1
parameter for level of output
= i: formatted output of mesh data until level i
= −i: unformatted output of mesh data until level i
name of input coarse mesh
name of created output coarse mesh
parameter for type of element refinement
= 0: no adaptive refinement strategy
= 1: element is partitioned into 9 finer elements
for: all elements allowed
DELMA: distance of element boundary to first refinement line (in percentage)
DELMB: = no use
usual: DELMA = 0.3333333 ∼ equidistant
= 2: element is partitioned (tangential to domain boundary) into 3 finer stripes
for: only elements belonging to 1 (!) boundary component allowed
DELMA: distance of domain boundary to first refinement line (in percentage)
DELMB: distance of domain boundary to second refinement line (in percentage)
usual: DELMA = 0.3333333, DELMA = 0.6666667 ∼ equidistant
= 3: element is partitioned (tangential to domain boundary) into 2 finer stripes
for: only elements belonging to 1 (!) boundary component allowed
DELMA: distance of domain boundary to first refinement line (in percentage)
DELMB: no use
usual: DELMA = 0.5 ∼ equidistant
= 4: element is partitioned (normal to domain boundary) into 3 finer stripes
for: only elements belonging to 1 (!) boundary component allowed
DELMA: distance of element boundary to first refinement line (in percentage)
DELMB: no use
usual: DELMA = 0.3333333 ∼ equidistant
= 5: element is partitioned into 5 finer elements
for: all elements allowed
DELMA: distance of element boundary to first refinement line (in percentage)
DELMB: no use
usual: DELMA = 0.3333333 ∼ equidistant
= 6: element is partitioned into 3 finer elements
for: only elements belonging to 2 (!) boundary component allowed
30
Structure of featflow
NELMOD
DELMA
DELMB
for: ”corner element”
DELMA: distance of element boundary to new inner vertex (in percentage)
DELMB: no use
usual: DELMA = 0.5 ∼ equidistant
name of elements to be adaptively refined
list of elements to be refined follows after DELMB
one element per line, separated by ”return”
parameter for distance in refinement process
= 0: no adaptive refinement strategy
= −i: unformatted output of mesh data until level i
parameter for distance in refinement process
= 0: no adaptive refinement strategy
= −i: unformatted output of mesh data until level i
}
DELMA
Figure 2.3: ITYPEL = 1: Refinement into 9 finer elements
BOUNDARY
}
INT.
INT.
}
DELMA
DELMB
INTERIOR
Figure 2.4: ITYPEL = 2: Refinement into 3 finer elements
BOUNDARY
}
B.
DELMA
B.
INT.
INT.
INTERIOR
Not
allowed
Figure 2.5: ITYPEL = 3: Refinement into 2 finer elements
!!!
2.2 The input data of featflow
31
BOUNDARY
B.
INT.
INT.
}
}
INTERIOR
DELMA
B.
Not
DELMA
allowed
!!!
Figure 2.6: ITYPEL = 4: Refinement into 3 finer elements
}
DELMA
Figure 2.7: ITYPEL = 5: Refinement into 5 finer elements
INTERIOR
B.
INT.
B.
}
DELMA
B.
B.
BOUNDARY
BOUNDARY
Not
allowed
!!!
Figure 2.8: ITYPEL = 6: Refinement into 3 finer elements
32
Structure of featflow
2.2.4. tr2to3.dat – parameter file for tr2to3
M
MT
ICHECK
IMESH
CPARM
IBDCHK
IBYU
IAVS
CFILEI
CFILEO
NPZ
PZMIN
PZMAX
feat parameter for output
= 0: no output
> 0: output, see feat2d manual
feat parameter for terminal output
= 0: no output
> 0: output, see feat2d manual
feat parameter for subroutine tracing
= 0: no tracing
> 0: tracing, see feat2d manual
parameter for type of parametrization
= 0: feat parametrization
= 1: omega2d parametrization
name of omega2d parametrization file
parameter for boundary checking for further refinements
= 0: no checking
> 0: = number of finer boundary points for check of boundary consistency for further refinement
level for byu output
level for avs output
name of input coarse mesh
name of created output coarse mesh
parameter for number of z-slides
list of interior z-coordinates follows after PZMAX
one element per line, separated by ”return”
parameter for minimum z-coordinate
parameter for maximum z-coordinate
2.2 The input data of featflow
2.2.5. trigen3d.dat – parameter file for trigen3d
M
MT
ICHECK
IBYU
IAVS
NLEV
IFMT
CFILEI
CFILEO
feat parameter for output
= 0: no output
> 0: output, see feat2d manual
feat parameter for terminal output
= 0: no output
> 0: output, see feat2d manual
feat parameter for subroutine tracing
= 0: no tracing
> 0: tracing, see feat2d manual
level for byu output
level for avs output
number of refined levels-1
parameter for level of output
= i: formatted output of mesh data until level i
= −i: unformatted output of mesh data until level i
name of input coarse mesh
name of created output coarse mesh
33
34
Structure of featflow
2.2.6. pp2d.dat/pp3d.dat – parameter files for pp2d and pp3d
The input parameter for the 2D- and the 3D version are almost identical. The parameter
files pp2d.dat/pp3d.dat and cc2d.dat/cc3d.dat differ with respect to a few parameters
only.
Specials in pp2d.dat/pp3d.dat:
IMASSL, ICUBML, ISORTU, ICYCU, ILMINU, ILMAXU, IINTU, ISMU, ISLU, NSMU,
NSLU, NSMUFA, ISORTP, ICYCP, ILMINP, ILMAXP, IINTP, ISMP, ISLP, NSMP,
NSLP, NSMPFA, EPSUR, EPSUD, DMPUD, DMPUMG, DMPUSL, RLXSMU, RLXSLU, AMINU,
AMAXU, EPSP, DMPPMG, DMPPSL, RLXSMP, AMINP, AMAXP, IPROJ, PRDIF1, PRDIF2
IMESH
IRMESH
CPARM
CMESH
CFILE
ISTART
CSTART
ISOL
CSOL
M
MT
ICHECK
MSHOW
NLMIN
NLMAX
IELT
parameter for type of parametrization
= 0: feat parametrization
= 1: omega2d parametrization (active only in pp2d!!!)
input of mesh data
= 0: create mesh
> 0: read mesh (created by trigen2d)
> 1: formatted mesh data (= 1 unformatted)
name of omega2d parametrization file
name of coarse mesh (not longer than 15 characters)
name of protocol file
input of start vector
= 0: start with homogeneous vector (only b.c.’s)
= 1: read (unformatted) start vector from same level (= −1: formatted)
= 2: read (unformatted) start vector from level-1 (= −2: formatted)
name of start vector file (not longer than 15 characters)
output of solution vector
= 0: no output
= 1: unformatted output on finest level
> 1: formatted output on finest level
name of solution vector file (not longer than 15 characters)
feat parameter for output
= 0: no output
> 0: output, see feat2d manual
feat parameter for terminal output
= 0: no output
> 0: output, see feat2d manual
feat parameter for subroutine tracing
= 0: no tracing
> 0: tracing, see feat2d manual
parameter for protocol
= 0: reduced protocol in file CFILE
= 1: expanded protocol in file CFILE
= 2: expanded protocol in file CFILE and terminal
= 3: full protocol in file CFILE
= 4: full protocol in file CFILE and terminal
parameter for smallest level number
parameter for highest level number
parameter for element type
2.2 The input data of featflow
ISTOK
IRHS
IBDR
IERANA
IMASS
IMASSL
IUPW
IPRECA
IPRECB
ICUBML
ICUBM
ICUBA
ICUBN
ICUBB
ICUBF
INLMIN
INLMAX
= 0: E031 parametric
= 1: E030 parametric
= 2: E031 nonparametric
= 3: E030 nonparametric
parameter for Stokes calculation
= 1: Stokes calculation
<> 1: Navier–Stokes calculation
parameter for right hand side
= 0: homogeneous r.h.s
= 1: steady inhomogeneous r.h.s
= 2: nonsteady inhomogeneous r.h.s
parameter for boundary
= 0: Dirichlet + Neumann boundary values
= 1: Dirichlet + Neumann + pressure drop boundary values
= 2: = 1 + time dependent Dirichlet/Neumann conditions
parameter for error analysis
= 0: no error analysis
> 0: error analysis with quadrature formula IERANA
parameter for mass matrix type
= 0: lumped mass matrix
= 1: real mass matrix
parameter for element type of lumped mass matrix
= 0: E031 parametric
= 1: E030 parametric
= 2: E031 nonparametric
= 3: E030 nonparametric
parameter for convective terms
= 0: streamline diffusion
= 1: upwinding
parameter for diffusive and reactive matrices
= 0: single precision in RAM
= 1: double precision in RAM
= 2: single precision from hard disc
= 3: double precision from hard disc
= 4: double precision built up every time
parameter for gradient and divergence matrices
= 0: single precision with usual quadrature, matrix in RAM
= 1: double precision with usual quadrature, matrix in RAM
= 2: double precision, exact matrix entries, elementwise application
= 3: single precision with exact evaluation, matrix in RAM
= 4: double precision with exact evaluation, matrix in RAM
quadrature formula for lumped mass matrix
> 0: usual lumping
< 0: diagonal lumping
quadrature formula for real mass matrix
quadrature formula for Laplacian matrix
quadrature formula for convective matrix
quadrature formula for gradient matrix
quadrature formula for right hand side
parameter for smallest number of nonlinear steps matrix
= 1: linear extrapolation in time (if INLMAX=1)
= −1: constant extrapolation in time (if INLMAX=-1)
parameter for largest number of nonlinear steps matrix
= 1: linear extrapolation in time (if INLMIN=1)
= −1: constant extrapolation in time (if INLMIN=-1)
35
36
ISORTU
ICYCU
ILMINU
ILMAXU
IINTU
ISMU
ISLU
NSMU
NSLU
NSMUFA
ISORTP
ICYCP
ILMINP
ILMAXP
IINTP
ISMP
ISLP
NSMP
NSLP
Structure of featflow
parameter for renumbering of edges
= 1: with respect to x–coordinate
= 2: with respect to y–coordinate
= 3: with Cuthill-McKee algorithm
parameter for mg-cycle for velocity
= 0: F–cycle
= 1: V–cycle
= 2: W–cycle
parameter for smallest number of mg-steps for velocity
parameter for largest number of mg-steps for velocity
parameter for mg-interpolation for velocity
= 1: E031 (parametric or nonparametric)
= 2: E030 (parametric or nonparametric)
parameter for mg-smoother for velocity
= 1: Jacobi
= 2: SOR
= 3: SSOR
= 4: ILU
parameter for mg-solver for velocity
= 1: SOR
= 2: BiCGSTAB
= 3: ILU
= 4: BiCGSTAB + ILU
parameter for number of mg-smoothing steps for velocity
parameter for number of mg-solving steps for velocity
parameter for change of number of mg-smoothing steps on coarser levels for velocity
= n: NSMU ∗ nNLMAX−ILEV smoothing steps if on level ILEV
parameter for renumbering of elements
= 1: with respect to x–coordinate
= 2: with respect to y–coordinate
= 3: with Cuthill-McKee algorithm
parameter for mg-cycle for pressure
= 0: F–cycle
= 1: V–cycle
= 2: W–cycle
parameter for smallest number of mg-steps for pressure
parameter for largest number of mg-steps for pressure
parameter for mg-interpolation for pressure
= 1: constant interpolation
= 2: rotated trilinear interpolation
= 3: trilinear interpolation
= 3: modified + optimized rotated trilinear interpolation
parameter for mg-smoother for pressure
= 1: Jacobi
= 2: SOR
= 3: SSOR
= 4: ILU
parameter for mg-solver for pressure
= 1: SOR
= 2: CG
= 3: ILU
= 4: CG + ILU
parameter for number of mg-smoothing steps for pressure
parameter for number of mg-solving steps for pressure
2.2 The input data of featflow
37
NSMPFA
parameter for change of number of mg-smoothing steps on coarser levels for pressure
= n: NSMP ∗ nNLMAX−ILEV smoothing steps if on level ILEV
RE
UPSAM
parameter for viscosity 1/NU
parameter for convectice discretization
> 0: upwind: parameter for Samarskij upwinding (usual: UPSAM=1)
> 0: SD: leading coefficient with spatial adaptivity (usual: UPSAM=1)
< 0: upwind: simple first order upwinding
< 0: SD: leading coefficient without spatial adaptivity (usual: UPSAM=-1)
parameter for lower limit for optimal relaxation in nonlinear iteration
> 0: relative changes are calculated
< 0: no relative changes are calculated (if OMGMIN = OMGMAX)
parameter for upper limit for optimal relaxation in nonlinear iteration
> 0: relative changes are calculated
< 0: no relative changes are calculated (if OMGMIN = OMGMAX)
parameter for start value for calculation of optimal relaxation in nonlinear iteration
stopping criterion for relative changes in velocity in nonlinear iteration
stopping criterion for defect in velocity in nonlinear iteration
stopping criterion for defect improvement in velocity in nonlinear iteration
stopping criterion for defect improvement in velocity in linear mg-iteration
stopping criterion for defect improvement for solver in velocity in linear mg-iteration
relaxation parameter for mg-smoother in velocity
relaxation parameter for mg-solver in velocity
lower limit for optimal correction for mg-solver in velocity
upper limit for optimal correction for mg-solver in velocity
stopping criterion for divergence of velocity in pressure equation
stopping criterion for defect improvement in pressure in linear mg-iteration
stopping criterion for defect improvement for solver in pressure in linear mg-iteration
relaxation parameter for mg-smoother in pressure
relaxation parameter for mg-solver in pressure
lower limit for optimal correction for mg-solver in pressure
upper limit for optimal correction for mg-solver in pressure
OMGMIN
OMGMAX
OMGINI
EPSUR
EPSUD
DMPUD
DMPUMG
DMPUSL
RLXSMU
RLXSLU
AMINU
AMAXU
EPSP
DMPPMG
DMPPSL
RLXSMP
RLXSLP
AMINP
AMAXP
IPROJ
NITNS
EPSNS
TIMENS
THETA
TSTEP
IFRSTP
INSAV
INSAVN
DTFILM
DTAVS
DTBYU
IFUSAV
parameter for type of projection scheme
= 0: first order (Chorin)
> 0: second order (Van Kan)
< 0: ABS(IPROJ) steps of first order, then second order
maximum number of macro time steps
stopping criterion for time derivative
parameter for absolute start time
parameter for time stepping value
= 1: Implicit Euler (first order), if IFRSTP=0
= 0.5: Crank–Nicolson (second order), if IFRSTP=0
starting time step
parameter for time stepping scheme
= 0: one step schemes
= 1: fractional step scheme (second order)
parameter for unformatted saving of solution vector
= 0: no saving
> 0: macro step size for saving procedure (in #ns)
modulo number of files for unformatted saving of solution vector (maximum = 10)
time difference for unformatted film output
time difference for avs output
time difference for byu output
level for unformatted velocity film output
38
IFPSAV
IFXSAV
IAVS
IBYU
IFINIT
IADTIM
TIMEMX
DTMIN
DTMAX
DTFAC
TIMEIN
EPSADI
EPSADL
EPSADU
IEPSAD
IADIN
IREPIT
PRDIF1
PRDIF2
Structure of featflow
level for unformatted pressure film output
level for unformatted streamline film output
active only in pp2d!!!
level for avs output
level for byu output
start file number for film output
parameter for adaptive time step control
= 0: no control, fixed time step TSTEP is used
= 1: prediction without repetition
= 2: prediction with repetition, if nonlinear stopping criteria too large
= 3: prediction with repetition, if time error or nonlinear stopping criteria too large
< 0: the same as ABS(IADTIM), but with extrapolation in time
maximum absolute time for stopping
parameter for smallest time step during adaptive control
parameter for largest time step during adaptive control
factor for largest possible time step changes
absolute time for start procedure
parameter for time error limit in start phase
parameter for time error limit after start phase
upper limit for acceptance for ABS(IADTIM)=3
parameter for type of error control
= 1: control of u(L2)
= 2: control of u(MAX)
= 3: control of p(L2)
= 4: control of p(MAX)
= 5: control of MAX(u(L2),p(L2))
= 6: control of MAX(u(MAX),p(MAX))
= 7: control of MAX(u(L2),p(L2),u(MAX),p(MAX))
= 8: control of MIN(u(L2),p(L2),u(MAX),p(MAX))
parameter for error control in start phase
= 0: EPSADL=EPSADI for T < TIMEIN
= 1: EPSADL=linear combination(EPSADI,EPSADL) for T < TIMEIN
= 2: EPSADL=logarithmic combination(EPSADI,EPSADL) for T < TIMEIN
maximim number of repetitions for ABS(IADTIM)=3
parameter for reactive preconditioner (usual: PRDIF1=1)
parameter for diffusive preconditioner (usual: PRDIF2=0 or 1)
2.2 The input data of featflow
39
2.2.7. cc2d.dat/cc3d.dat – parameter files for cc2d and cc3d
The input parameter for the 2D- and the 3D version are almost identical. The parameter
files cc2d.dat/cc3d.dat and pp2d.dat/pp3d.dat differ with respect to a few parameters
only.
Specials in cc2d.dat/cc3d.dat:
IMASSL, ICYCLE, ILMIN, ILMAX, IINT, ISM, ISL, NSM, NSL, NSMFAC,
EPSD, EPSDIV, EPSUR, EPSPR, DMPD, DMPMG, EPSMG, DMPSL, EPSSL,
RLXSM, RLXSL, AMINMG, AMAXMG, ISTAT
IMESH
IRMESH
CPARM
CMESH
CFILE
ISTART
CSTART
ISOL
CSOL
M
MT
ICHECK
MSHOW
NLMIN
NLMAX
IELT
parameter for type of parametrization
= 0: feat parametrization
= 1: omega2d parametrization (active only in cc2d!!!)
input of mesh data
= 0: create mesh
> 0: read mesh (created by trigen2d)
> 1: formatted mesh data (= 1 unformatted)
name of omega2d parametrization file
name of coarse mesh (not longer than 15 characters)
name of protocol file
input of start vector
= 0: start with homogeneous vector (only b.c.’s)
= 1: read (unformatted) start vector from same level (= −1: formatted)
= 2: read (unformatted) start vector from level-1 (= −2: formatted)
name of start vector file (not longer than 15 characters)
output of solution vector
= 0: no output
= 1: unformatted output on finest level
> 1: formatted output on finest level
name of solution vector file (not longer than 15 characters)
feat parameter for output
= 0: no output
> 0: output, see feat2d manual
feat parameter for terminal output
= 0: no output
> 0: output, see feat2d manual
feat parameter for subroutine tracing
= 0: no tracing
> 0: tracing, see feat2d manual
parameter for protocol
= 0: reduced protocol in file CFILE
= 1: expanded protocol in file CFILE
= 2: expanded protocol in file CFILE and terminal
= 3: full protocol in file CFILE
= 4: full protocol in file CFILE and terminal
parameter
parameter
parameter
= 0: E031
= 1: E030
for smallest level number
for highest level number
for element type
parametric
parametric
40
ISTOK
IRHS
IBDR
IERANA
IMASS
IMASSL
IUPW
IPRECA
IPRECB
ICUBM
ICUBA
ICUBN
ICUBB
ICUBF
INLMIN
INLMAX
ICYCLE
ILMIN
ILMAX
IINT
Structure of featflow
= 2: E031 nonparametric
= 3: E030 nonparametric
parameter for Stokes calculation
= 1: Stokes calculation
<> 1: Navier–Stokes calculation
parameter for right hand side
= 0: homogeneous r.h.s
= 1: steady inhomogeneous r.h.s
= 2: nonsteady inhomogeneous r.h.s
parameter for boundary
= 0: Dirichlet + Neumann boundary values
= 1: Dirichlet + Neumann + pressure drop boundary values
= 2: = 1 + time dependent Dirichlet/Neumann conditions
parameter for error analysis
= 0: no error analysis
> 0: error analysis with quadrature formula IERANA
parameter for mass matrix type
= 0: lumped mass matrix
= 1: real mass matrix
parameter for mass matrix lumping
= 0: usual lumping
= 10: diagonal lumping
parameter for convective terms
= 0: streamline diffusion
= 1: upwinding
parameter for diffusive and reactive matrices
= 0: single precision in RAM
= 1: double precision in RAM
= 2: single precision from hard disc
= 3: double precision from hard disc (not yet)
= 4: double precision built up every time
parameter for gradient and divergence matrices
= 0: single precision with usual quadrature, matrix in RAM
= 1: double precision with usual quadrature, matrix in RAM (not yet)
= 2: double precision, exact matrix entries, elementwise application (not yet)
= 3: single precision with exact evaluation, matrix in RAM
= 4: double precision with exact evaluation, matrix in RAM (not yet)
quadrature formula for real mass matrix
quadrature formula for Laplacian matrix
quadrature formula for convective matrix
quadrature formula for gradient matrix
quadrature formula for right hand side
parameter for smallest number of nonlinear steps matrix
= 1: linear extrapolation in time (if INLMAX=1)
= −1: constant extrapolation in time (if INLMAX=-1)
parameter for largest number of nonlinear steps matrix
= 1: linear extrapolation in time (if INLMIN=1)
= −1: constant extrapolation in time (if INLMIN=-1)
parameter for mg-cycle
= 0: F–cycle
= 1: V–cycle
= 2: W–cycle
parameter for smallest number of mg-steps
parameter for largest number of mg-steps
parameter for mg-interpolation
2.2 The input data of featflow
ISM
ISL
NSM
NSL
NSMFAC
RE
UPSAM
OMGMIN
OMGMAX
OMGINI
EPSD
EPSDIV
EPSUR
EPSPR
DMPD
DMPMG
EPSMG
DMPSL
EPSSL
RLXSM
RLXSL
AMINMG
AMAXMG
ISTAT
NITNS
EPSNS
TIMENS
THETA
TSTEP
IFRSTP
INSAV
INSAVN
41
= 1: rotated trilinear for velocity + constant for pressure
= 2: rotated trilinear for velocity + pressure
parameter for mg-smoother
= 1: Vanca
parameter for mg-solver
= 1: Vanca
parameter for number of mg-smoothing steps
parameter for number of mg-solving steps
parameter for change of number of mg-smoothing steps on coarser levels
= n: NSM ∗ nNLMAX−ILEV smoothing steps if on level ILEV
parameter for viscosity 1/NU
parameter for convectice discretization
> 0: upwind: parameter for Samarskij upwinding (usual: UPSAM=1)
> 0: SD: leading coefficient with spatial adaptivity (usual: UPSAM=1)
< 0: upwind: simple first order upwinding
< 0: SD: leading coefficient without spatial adaptivity (usual: UPSAM=-1)
parameter for lower limit for optimal relaxation in nonlinear iteration
> 0: relative changes are calculated
< 0: no relative changes are calculated (if OMGMIN = OMGMAX)
parameter for upper limit for optimal relaxation in nonlinear iteration
> 0: relative changes are calculated
< 0: no relative changes are calculated (if OMGMIN = OMGMAX)
parameter for start value for calculation of optimal relaxation in nonlinear iteration
stopping criterion for defect in velocity in nonlinear iteration
stopping criterion for defect in divergence in nonlinear iteration
stopping criterion for relative changes in velocity in nonlinear iteration
stopping criterion for relative changes in pressure in nonlinear iteration
stopping criterion for defect improvement in nonlinear iteration
stopping criterion for defect improvement in linear mg-iteration
stopping criterion for defect limit in linear mg-iteration
stopping criterion for defect improvement for solver in linear mg-iteration
stopping criterion for defect limit for solver in linear mg-iteration
relaxation parameter for mg-smoother
relaxation parameter for mg-solver
lower limit for optimal correction for mg-solver
upper limit for optimal correction for mg-solver
parameter for type of problem
= 0: steady Navier-Stokes problem
= 1: nonsteady Navier-Stokes problem
maximum number of macro time steps
stopping criterion for time derivative
parameter for absolute start time
parameter for time stepping value
= 1: Implicit Euler (first order), if IFRSTP=0
= 0.5: Crank–Nicolson (second order), if IFRSTP=0
starting time step
parameter for time stepping scheme
= 0: one step schemes
= 1: fractional step scheme (second order)
parameter for unformatted saving of solution vector
= 0: no saving
> 0: macro step size for saving procedure (in #ns)
modulo number of files for unformatted saving of solution vector (maximum = 10)
42
DTFILM
DTAVS
DTBYU
IFUSAV
IFPSAV
IFXSAV
IAVS
IBYU
IFINIT
IADTIM
TIMEMX
DTMIN
DTMAX
DTFAC
TIMEIN
EPSADI
EPSADL
EPSADU
IEPSAD
IADIN
IREPIT
Structure of featflow
time difference for unformatted film output
time difference for avs output
time difference for byu output
level for unformatted velocity film output
level for unformatted pressure film output
level for unformatted streamline film output
active only in pp2d!!!
level for avs output
level for byu output
start file number for film output
parameter for adaptive time step control
= 0: no control, fixed time step TSTEP is used
= 1: prediction without repetition
= 2: prediction with repetition, if nonlinear stopping criteria too large
= 3: prediction with repetition, if time error or nonlinear stopping criteria too large
< 0: the same as ABS(IADTIM), but with extrapolation in time
maximum absolute time for stopping
parameter for smallest time step during adaptive control
parameter for largest time step during adaptive control
factor for largest possible time step changes
absolute time for start procedure
parameter for time error limit in start phase
parameter for time error limit after start phase
upper limit for acceptance for ABS(IADTIM)=3
parameter for type of error control
= 1: control of u(L2)
= 2: control of u(MAX)
= 3: control of p(L2)
= 4: control of p(MAX)
= 5: control of MAX(u(L2),p(L2))
= 6: control of MAX(u(MAX),p(MAX))
= 7: control of MAX(u(L2),p(L2),u(MAX),p(MAX))
= 8: control of MIN(u(L2),p(L2),u(MAX),p(MAX))
parameter for error control in start phase
= 0: EPSADL=EPSADI for T < TIMEIN
= 1: EPSADL=linear combination(EPSADI,EPSADL) for T < TIMEIN
= 2: EPSADL=logarithmic combination(EPSADI,EPSADL) for T < TIMEIN
maximim number of repetitions for ABS(IADTIM)=3
2.3 The file structure of featflow
43
2.3. The file structure of featflow
In reversed order we arrive at the point to explain the internal structure of featflow.
featflow consists of 6 directories: application, graphic, manual, object, source and
utility. We want to explain their tasks in detail, not following the literal order.
2.3.1. Subdirectory source – source code for featflow
The directory source contains (almost) all source code for the preprocessing and solver
tools. These are omega2d, trigen2d, trigen3d, tr2to3, intpol2d, intpol3d, cc2d,
cc3d, pp2d and pp3d. Each of them is split into two, resp., three directories, namely src
and dev, resp., mg.
The idea is that src contains all information to build up the corresponding system libraries
(see object and the instructions for installation). These are needed for running an user–
application. mg is similar and is needed for building the corresponding libraries containing
the multigrid components for solving (if necessary). They are only needed by the solvers
cc2d, cc3d, pp2d and pp3d.
The directory dev is a corresponding directory containing all source- and makefiles, and
it is thought to be a ”developer directory” which is needed by the experienced user to
modify the featflow software.
Furthermore, the source directory contains the code for blas, feat2d, feat3d and
omega2d. In an analogous way, src contains all files to build up the corresponding system
libraries during installation. Additionally, there are two testdir directories containing a
small test program to get familiar with feat. Additionally, we think about adding the
GNU FORTRAN77 compiler to provide a FORTRAN compiler for everybody.
2.3.2. Subdirectory manual – manuals for featflow
The directory manual contains (almost) all LATEX–source files (in src) and postscript–
files (in ps) for needed manuals. These are omega2d, feat2d, feat3d and this featflow manual.
2.3.3. Subdirectory object – system software for featflow
The directory object contains the system software of featflow which is needed for
installation, for user–applications, and for generating a (compressed) binary data–file containing the featflow code. There are three subdirectories, extract, makefiles and
libraries.
extract is a directory having shell scripts for generating a (compressed) tarfile containing
the featflow package. There are several levels of featflow data files:
44
Level
Level
Level
Level
Level
Structure of featflow
00
01
02
10
20
complete featflow (even with all libraries)
complete featflow sources (but without compiled libraries)
partial featflow sources (without dev directories)
partial featflow sources + libraries (without multigrid sources)
featflow libraries (without any sources)
The subdirectory makefiles contains all makefiles for building up the system libraries and
for user–applications. They are available for a class of machines (see the file MACHINES) and
have to be modified for computers or compilers not belonging to this list. Each of these
machine–dependent directories contains some shell scripts, beginning with make , which
have to be executed for installation (see the section on installation). In the next versions,
there will be some more computers added to our list, containing optimized machine–
dependent compiler options.
Finally, the directory libraries contains the compiled system libraries which have to be
linked to each user–application. The actual featflow libraries created by the installation (depending on the makefile used) is in the subdirectory libgen, while in libspec
machine- and compiler–dependent system libraries are offered in correspondance to the
special makefiles (see above).
2.3.4. Subdirectory application – applications under featflow
application is a directory thought to be the place for user–applications. The idea is to
use a separate subdirectory for each application. New subdirectories should be generated
by linking or copying, at least the input files and the #data directory.
featflow is installed with 5 subdirectories: user start, comp, dat example, example
and workspace.
user start is a typical working directory for the user, containing all needed input parameter files in #data and a special subdirectory input files. In this directory, the user
can find all makefiles, copied during installation and adapted to the machine used, and
the input files parq2d.f, parq3d.f indat2d.f and indat3d.f (see above). These files
have to be edited and modified according to the actual application, and the corresponding
makefile has to be executed. Additionally, before compiling, the corresponding storage has
to be defined, by defining the corresponding .inc file (see above and the later subsection
on the workspace subdirectory).
After compiling the tool (by executing the makefile) the compiled program should be
copied or moved to the parent directory (= user start), and the application may be
started there, after editing the corresponding parameter file in #data. All other applications should follow these instructions to make life easier.
In a similar way, comp is a pre–installed application directory, performing test calculcations
(similar to the DFG–benchmark configuration of 1995, see [16]). These applications are
thought to provide reference results and CPU times for verifying the installed featflow
version and to obtain benchmark results for different computer types for comparisons.
These results can be found in results.
2.4 Installation of featflow
45
example is a pre–installed application directory, containing the example programs of Chapter 3.
Finally, there is the directory workspace containing some shell scripts. Their only use is to
select the corresponding include–files containing the size of storage amount. For instance, if
the directory user start is in progress, a link has to be set to the corresponding .inc files
in the directory user start. This is done by executing the appropriate links.user start
file (or other files analogously). This process has to be done before compiling.
2.3.5. Subdirectory graphic – graphic tools for featflow
This directory provides graphic features which are helpful for featflow. In the present
version there are the subdirectories avs, byu and gnuplot. These features will be massively
expanded in future versions.
avs contains a .avsrc file (as link to avsrc) and a directory appl containing special
applications. This .avsrc file has to be edited and adapted to the existing configuration.
Up to now, avs has to be executed over network (if available), and the manual has to be
taken from there, too.
The subdirectories byu and gnuplot contain some helpful files for handling the packages
movie.byu or cquel.byu (over network, if available) and for working with gnuplot,
which may be used for 1D graphics.
2.3.6. Subdirectory utility – utilities for featflow
This directory contains software which can be helpful for the use of featflow. In the
present version there are only programs which provide an a priori–estimate for the needed
storage amount. They simply have to be compiled (f77 file.f -o file) and to be
moved to the directory used for the application, for instance to application/user start.
Then, it reads the corresponding parameter file and gives an estimate for the NNWORK value
needed.
2.4. Installation of featflow
The usual installation process is the following:
Step 1: Unzip and tar
Create a directory called featflow or similar, which will contain all featflow data.
Usually, the user has obtained a binary file featflow.tar.gz or levi file.tar.gz (”i”
stands for level, see above). This has to be moved to the created directory in which the
featflow tree structure will be installed. It has to be decompressed
gunzip featflow.tar.gz
46
Structure of featflow
and then, a tar has to be started
tar -xvf featflow.tar.
This generates the complete featflow tree structure, and the file featflow.tar may be
removed.
In most cases, all makefiles are already prepared if you work on a ‘normal’ platform (SUN,
IBM, SGI, DEC,, PC+LINUX, etc.) such that the following installation step is very easy:
Step 2a: Automatic installation
Read the README file and execute the installation script:
install feat
After selecting the supported platform, the installation procedure starts! You may continue with Step 3
If you work on a ‘nonstandard’ platform (whatsoever), you must manually perform the
following two steps:
Step 2b: Editing of makefiles and installation scripts
In the next step go to featflow/object/makefiles. There are different subdirectories
corresponding to various computer types and compiler options. If not, copy one to create
your owen, or modify the makefiles in example. Let us assume you are in the directory
featflow/object/makefiles/example. There wait two tasks for you:
1) Edit the shell scripts make copy and make lib, and change the shell variable FEATFLOW,
containing the right location of featflow, and the variable MAKEFILES where you are.
2) Edit the makefiles and modify, if necessary, the compiler options. These are the
lines beginning with COMOPT =. This can be done by hand, or by using the shell script
make change. Analogously, the shell variable FEATFLOW has to be defined correctly as
above, using make change.
Step 2c: Installation
Assuming that you are in featflow/object/makefiles/examples (or any other directory
containing your edited makefiles and shell scripts of step 2). Then, you have to do the
following:
0.1) Be sure that you use a C-shell /bin/csh.
0.2) Be sure that featflow/object/libraries/libgen exists.
0.3) Be sure that the correct ztime.f file in featflow/source/feat2d/src is taken.
1.0) Execute the shell script make copy.
2.0) Execute the shell script make lib.
2.4 Installation of featflow
47
This last step takes between 10 and 60 minutes. All libraries will be generated, and
all makefiles needed for applications are copied to featflow/application/user start,
featflow/application/example and featflow/application/comp.
Step 3: Test and application
As described in subsection 2.3.4 go to featflow/application/user start, to featflow
/application/example or to featflow/application/comp, and start an application, or
create a new directory for your application as described above.
3. Examples for the use of featflow
This chapter demonstrates explicitely the use of featflow, performing first the installation, and then applying featflow for solving a 2D and 3D problem which represents the
typical application of a flow in a channel around a cylinder. The values to be computed are
lift and drag and pressure differences on the cylinder surface, in a stationary as well as nonstationary configuration. These examples are almost identical with the DFG–benchmark
1995, see [16]), only the 2D case is a little bit different (longer channel!!!).
We show (in short form) how featflow is manually installed, followed by designing a
first coarse mesh with omega2d. This triangulation is refined in a pre–adaptive way by
trigen2d, generating the full sequence of nested meshes which are neccessary for the
multigrid solvers. We perform a stationary (by cc2d) and a nonstationary calculation
(by pp2d), and demonstrate how to interpret the data output and how to use avs for
visualization.
These 2D calculations are followed by a 3D application. We first use tr2to3 to construct
an adequate 3D coarse mesh for this channel configuration. We demonstrate the use
of trigen3d to generate the 3D multigrid structure and, finally, we perform the same
calculations with cc3d and pp3d as in the 2D case.
3.1. The installation example
We assume that you got the binary file featflow.tar.gz and that you created the directory /home/people/example/featflow in which you should move this binary file. Being
there, you have to perform the following commands:
gunzip featflow.tar.gz
tar -xvf featflow.tar
rm featflow.tar.
Now, the complete source code of featflow is installed.
In most cases, all makefiles are already prepared if you work on a ‘normal’ platform (SUN,
IBM, SGI, DEC, HP, PC+LINUX, etc.) such that the following installation step is very
easy:
Read the README file and execute the installation script:
install feat
48
3.1 The installation example
49
After selecting the supported platform, the installation procedure starts! You may continue with the 2D and 3D examples!
If you work on a ‘nonstandard’ platform (whatsoever), you must manually perform the
following two steps:
For the following, we assume you have a SUN Sparc 10 (or compatible) with the SUN
FORTRAN77 compiler, version 2.0. Next you go to /home/people/example/featflow/
object/makefiles, and do the following copy:
cp -rp sun sparc10 2.0 my compiler
cd my compiler.
Edit, for instance, the file pp2d.m to figure out how the shell–variables FEATFLOW and
COMOPT are defined. FEATFLOW should be set to:
FEATFLOW=/home/people/example/featflow,
while COMOPT has to be modified for other machines or compilers only. These changes
could be done by hand, or by using the shell–script make change which performs (after
appropriate replacements) an sed command on all files. It is necessary that you use
a C-shell /bin/csh. Additionally, the variable MAKEFILES in make lib should be set
appropriate. From now on, we use FEATFLOW as abbreviation for /home/people/example
/featflow.
Before we continue be sure that FEATFLOW/object/libraries/libgen exists (if not, create them by the mkdir command), and that the ”right” subroutine for time measurements
are taken. That means, for SUN SPARC 10 (and many others),
cd FEATFLOW/source/feat2d/src
cp ztime.sun sparc ztime.f
cd FEATFLOW/object/makefiles/my compiler.
Now, do the following:
Execute the shell script make copy by typing: make copy.
Execute the shell script make lib by typing: time make lib.
This last step takes between 10 and 60 minutes (the time command is not necessary).
All featflow libraries will be generated, and all makefiles needed for applications are
copied to FEATFLOW/application/example (and to FEATFLOW/application/user start
and FEATFLOW/application/comp).
Next, go to FEATFLOW/application/workspace and run the shell script links.example
which activates the applications in FEATFLOW/application/example (more precise: the
.inc files in FEATFLOW/application/example/input files will be used for defining the
storage amount).
Now, we are ready to start the 2D and 3D examples.
50
Examples for the use of featflow
3.2. The 2D example
We start with describing our 2D domain which shall look like:
U=V=0
(0,H)
0.16m
0.45m
U=V=0
outlet
0.1m
y
inlet
0.15m
U=V=0
x
(0,0)
2.5m
Figure 3.1: 2D domain: channel with cylinder
This can be easily done with omega2d, requiring a computer with PC emulation (at least
an AT 286 or higher) and WINDOWS 3.1 (or higher).
We start in omega2d with describing the two boundary components: first, we prescribe
the outer rectangle, by setting the first point at the origin, and then proceeding in counterclockwise sense, until this boundary ”curve” is closed. Second, we draw a circle, with starting point at (0.45, 0.20) and center point (0.50, 0.20). The end point is again at (0.45, 0.20),
and it is very important to perform a ”negative” circle, that means the parametrization
is proceeding in a clockwise sense. For more details, look at the omega2d manual.
Next we define mesh points at the boundary components and in the interior. It is remarkable that boundary points cannot be set arbitrarily, but next to another (existing)
boundary point, and then they can be moved along the boundary curve. Furthermore,
you are always able to ”undo” your last actions or to remove points (and even makros and
boundary components when needed).
In parallel to defining mesh points, ”makros” (that means quadrilateral elements) can be
defined, by clicking at four mesh points to form an element. Finally, this session has to be
saved (let us assume as c2d.geb, c2d.prm and c2d.tri), and the files c2d.geb, c2d.prm
and c2d.tri have to be transferred (by ftp or whatever) to the subdirectory FEATFLOW/
application/example/#pre, and additionally the coarse mesh c2d.tri as c2d0.tri to
FEATFLOW/application/example/#adc. Now, we are ready to proceed with featflow
on our UNIX workstation.
Next, we want to ”improve” our coarse mesh c2d0.tri to obtain a ”better” triangulation.
Therefore, we go toFEATFLOW/application/example/input files, execute the makefiles (ending with .m) and move the compiled objects to the parent directory, namely to
FEATFLOW/application/example. More in detail, for performing the 2D example:
3.2 The 2D example
51
Figure 3.2: Coarsest mesh c2d0.tri: 30 elements
cd FEATFLOW/application/example/input files
make -f trigen2d.m ; mv trigen2d ..
cd FEATFLOW/application/example.
Next, we want to refine our mesh c2d0.tri around and before and behind the circle. This
is done by trigen2d. If we perform
cp #data/trigen2d.dat 0 #data/trigen2d.dat
time trigen2d
we obtain, in #adc, a new coarse mesh c2d1.tri which is refined in an appropriate way
(with ITYPEL=1, that means the marked elements are refined into 9 elements). The next
steps generate the ”final” coarse mesh c2d2.tri which is additionally refined around the
circle only (with ITYPEL=3).
cp #data/trigen2d.dat 1 #data/trigen2d.dat
time trigen2d.
Figure 3.3: Refined coarse mesh c2d1.tri: 130 elements
Figure 3.4: Further refined coarse mesh c2d2.tri: 148 elements
52
Examples for the use of featflow
Now, we are ready to start a Navier–Stokes calculation. Let us begin with a stationary example, for Reynolds number 20. Doing that, we use the data file input files/
indat2d.f stat, and the following procedure has to be performed
cd FEATFLOW/application/example/input files
cp indat2d.f stat indat2d.f ; make -f cc2d.m ; mv cc2d ..
cd FEATFLOW/application/example.
This data file input files/indat2d.f stat contains the following definitions:
–
The inflow velocity profile at x = 0.0 is parabolic, with a maximum value umax = 0.3.
–
There is one boundary part containing natural boundary conditions (outflow !). This
is the boundary segment with x = 2.5, and the corresponding parameter values range
from DPARN1=1D0 until DPARN1=2D0.
–
We show the velocity values at the mesh points KPU(1)=264 (corresponds to the
coordinates (0.65, 0.20)) and KPU(2)=17 (= (0.85, 0.20)), the pressure values at
KPP(1)=27 (= (0.45, 0.20)), KPP(2)=30 (= (0.55, 0.20)), KPP(3)=316 (= (0.50, 0.25))
and KPP(4)=17 (= (0.85, 0.20)), and the flux between the mesh points KPX(1)=31 and
KPX(2)=6 (under the circle). The difference of the pressure values for KPP(1)=27 and
KPP(2)=30 determines a typical pressure difference. Additionally, our runtime protocol will give us the mean pressure over the whole circle (by setting DPI(1,1)=0D0
and DPI(2,1)=0D0), and over half of the circle (by DPI(2,2)=0.5D0). Finally, we define parameters for the calculation of drag and lift, where RHO is a density parameter,
DIST a typical length scale, and UMEAN is a ”mean velocity”.
Now, we can execute cc2d with given parameter file #data/cc2d.dat. A solution is calculated on level NLMAX=4, and the corresponding solution vector is saved as (unformatted)
file #data/#DX4 stat. The chosen discretization scheme for the convective parts is the
streamline–diffusion ansatz, and the stopping criterions are 10−3 for the maximum of relative changes. As result we obtain the protocol file #data/cc2d.stat which has to be
explained in more detail.
Protocol file #data/cc2d.stat:
-------------------------------------------------------------------------------INPUT DATA
-------------------------------------------------------------------------------Parametrization file =
#pre/c2d.prm
Coarse grid file
=
#adc/c2d2.tri
Integer parameters of /IPARM/,etc. :
-------------------------------------------------------------------------------minimum mg-level: NLMIN =
1
maximum mg-level: NLMAX =
4
element type
=
3
Stokes calculation: ISTOK =
0
RHS
generation
=
0
Boundary generation =
0
Error evaluation
=
0
mass evaluation
=
0
3.2 The 2D example
lumped mass eval.
=
0
convective part
=
0
Accuracy for ST
=
4
Accuracy for B
=
3
ICUB mass matrix
=
3
ICUB diff. matrix
=
4
ICUB conv. matrix
=
4
ICUB matrices B1,B2
=
8
ICUB right hand side
=
1
minimum of nonlinear iterations: INLMIN =
1
maximum of nonlinear iterations: INLMAX =
10
type of mg-cycle: ICYCLE =
0
minimum of linear mg steps : ILMIN =
1
maximum of linear mg steps : ILMAX =
5
type of interpolation: IINT =
2
type of smoother : ISM =
1
type of solver : ISL =
1
number of smoothing steps : NSM =
4
number of solver steps : NSL =
100
factor sm. steps on coarser lev.:NSMFAC= 2
KPRSM,KPOSM ON LEVEL:
1 32 32
KPRSM,KPOSM ON LEVEL:
2 16 16
KPRSM,KPOSM ON LEVEL:
3 8 8
KPRSM,KPOSM ON LEVEL:
4 4 4
KPRSM,KPOSM ON LEVEL:
5 4 4
KPRSM,KPOSM ON LEVEL:
6 4 4
KPRSM,KPOSM ON LEVEL:
7 4 4
KPRSM,KPOSM ON LEVEL:
8 4 4
KPRSM,KPOSM ON LEVEL:
9 4 4
Real parameters of /RPARM/,etc. :
-------------------------------------------------------------------------------Viscosity parameter: 1/NU =
1000.00000000000
parameter for Samarskij-upwind: UPSAM =
0.50000000000000
lower limit for optimal OMEGA: OMGMIN =
0.
upper limit for optimal OMEGA: OMGMAX =
2.0000000000000
start value for optimal OMEGA: OMGINI =
1.0000000000000
limit for U-defects
:
EPSD
=
1.0000000000000D-05
limit for DIV-defects :
EPSDIV =
1.0000000000000D-08
limit for U-changes :
EPSUR =
1.0000000000000D-03
limit for P-changes :
EPSPR =
1.0000000000000D-03
defect improvement :
DMPD
=
1.0000000000000D-01
damping of MG residuals
: DMPMG =
1.0000000000000D-01
limit for MG residuals
: EPSMG =
1.0000000000000D-01
damping of residuals for solving: DMPSL =
1.0000000000000D-01
limit of changes for solving:
EPSSL =
1.0000000000000D-01
relaxation for the U-smoother: RLXSM =
1.0000000000000
relaxation for the U-solver : RLXSL =
1.0000000000000
lower limit optimal MG-ALPHA: AMINMG =
-10.0000000000000
upper limit optimal MG-ALPHA: AMAXMG =
10.0000000000000
Parameters of /NS.../ :
-------------------------------------------------------------------------------Time dependency
: ISTAT =
0
Number of time steps
: NITNS =
9
limit for time derivative: EPSNS =
1.0000000000000D-05
Total time
: TIMENS =
0.
Theta
: THETA =
1.0000000000000
Time step
: TSTEP =
1.0000000000000D-02
Fractional step
: IFRSTP =
1
Stepsize for nonsteady savings: INSAV =
0
Number of files
: INSAVN =
0
Time step for Film
: DTFILM =
0.
Time step for AVS
: DTAVS =
1.0000000000000
Time step for BYU
: DTBYU =
1.0000000000000
Level for velocity
: IFUSAV =
0
Level for pressure
: IFPSAV =
0
Level for streamlines
: IFXSAV =
0
Level for AVS
: IAVS =
4
Level for BYU
: IBYU =
4
Start file
: IFINIT =
1
Type of adaptivity
: IADTIM =
-3
53
54
Examples for the use of featflow
Max. Time
: TIMEMX =
100.000000000000
Min. Timestep
: DTMIN =
1.0000000000000D-06
Max. Timestep
: DTMAX =
1.0000000010000
Max. Timestep change
: DTFACT =
9.0000000010000
Time for start procedure
: TIMEIN =
0.50000000000000
EPS for start procedure
: EPSADI =
0.12500000000000
EPS for acceptance
: EPSADL =
1.2500000000000D-03
EPS for not acceptance
: EPSADU =
0.50000000000000
Acceptance criterion
: IEPSAD =
1
Start procedure
: IADIN =
2
Max.numbers of repetitions
: IREPIT =
3
-------------------------------------------------------------------------------ILEV,NVT,NMT,NEL,NVBD:
1 165 313 148 34
ILEV,NVT,NMT,NEL,NVBD:
2 626 1218 592 68
ILEV,NVT,NMT,NEL,NVBD:
3 2436 4804 2368 136
ILEV,NVT,NMT,NEL,NVBD:
4 9608 19080 9472 272
time for grid generation :
ILEV,NU,NA,NB:
ILEV,NU,NA,NB:
ILEV,NU,NA,NB:
ILEV,NU,NA,NB:
1
2
3
4
3.6833333298564
313 2089 592
1218 8322 2368
4804 33220 9472
19080 132744 37888
time for initialization of linear operators :
2.3499999046326
-------------------------------------------------------------------------------IT RELU
RELP
DEF-U
DEF-DIV DEF-TOT RHONL
OMEGNL
RHOMG
-------------------------------------------------------------------------------0
0.58D-02 0.26D-01 0.26D-01
-------------------------------------------------------------------------------1 0.10D+01 0.10D+01 0.34D-02 0.34D-04 0.34D-02 0.13D+00 0.99D+00 0.20D+00
2 0.41D+00 0.49D+00 0.25D-03 0.85D-05 0.25D-03 0.98D-01 0.10D+01 0.87D-01
3 0.21D+00 0.14D+00 0.87D-04 0.82D-06 0.87D-04 0.15D+00 0.11D+01 0.24D+00
4 0.64D-01 0.52D-01 0.20D-04 0.27D-06 0.20D-04 0.17D+00 0.11D+01 0.31D+00
5 0.13D-01 0.45D-02 0.58D-05 0.44D-07 0.58D-05 0.19D+00 0.11D+01 0.33D+00
6 0.33D-02 0.19D-02 0.14D-05 0.16D-07 0.14D-05 0.19D+00 0.11D+01 0.29D+00
7 0.62D-03 0.39D-03 0.33D-06 0.38D-08 0.33D-06 0.20D+00 0.11D+01 0.31D+00
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
1(
1) TIME= 0.000D+00 NORM(U)= 0.2092401D+00 NORM(P)= 0.4859169D-01
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
P(VELO)
P(PRES)
I(PRES)
I(FORCE)
P(FLUX)
0.20913D-01
0.12008D+00
0.47304D-01
0.53363D+01
0.39419D-01
0.14274D-02 0.16732D+00 0.14309D-02
0.14282D-01-0.81432D-02 0.28207D-01
0.29039D-01
0.66688D-02
STATISTICS :
NWORK :
1000000
IWORKG:
116168
IWMAXG:
119323
IWORKI:
685234
IWMAXI:
685234
IWORK :
685234
IWMAX :
704450
-------------------------------------------------------------------------------total time :
212.94999999553
appr. time :
212.83334153146
grid time :
3.6833333298564
post time :
7.5666656494141
lin. time :
49.299993515015
-> mavec time :
6.1333212852478
-> konv. time :
39.716608047485
-> bdry time :
4.9992084503174D-02
-> LC
time :
3.4000720977783
mg
time :
152.28334903717
#substeps :
1
3.2 The 2D example
55
#nonlinear :
7
#mg
:
14
-------------------------------------------------------------------------------MULTIGRID COMPONENTS [in percent]:
smoothing
:
74.324147510282
solver
:
6.0085270812218
defect calc. :
8.7993897069504
prolongation :
9.7515820090053
restriction
:
1.0944625011425
--------------------------------------------------------------------------------
Most statement have not to be explained, only a few ones:
–
ILEV,NVT,NMT,NEL,NVBD means number of vertices/midpoints/elements/boundary
points on level ILEV
–
ILEV,NU,NA,NB means number of midpoints and nonzero matrix entries for the
Laplacian/gradient matrix on level ILEV
–
RHONL convergence rate of nonlinear iteration
–
OMEGNL optimally chosen relaxation parameter for nonlinear iteration
–
RHOMG multigrid convergence rate for Oseen equation
–
P(VELO) contains the u- and v–velocity for both grid points defined in input files
/indat2d.f.
–
P(PRES) contains the 4 pressure values for the grid points defined before.
–
I(PRES) contains the 2 integral mean pressure values defined before.
–
I(FORCE) contains the drag and lift values defined before.
–
P(FLUX) contains the flux value defined before.
–
IWMAX contains the value for NNWORK needed. NWORK shows the actually defined
parameter.
Additionally, we obtain files for graphical output, namely in #avs the file u1.inp, and in
#byu the byu files u1.vec (velocity), p1.scl (pressure) and x1.vec (streamfunction). A
typical example for a avs isoline plot of the pressure can be found in the following picture.
Figure 3.5: Pressure isolines for the stationary 2D calculation
56
Examples for the use of featflow
Next, we perform a nonstationary calculation, for Reynolds number about 100. We use the
data file input files/indat2d.f non, and the following procedure has to be performed
cd FEATFLOW/application/example/input files
cp indat2d.f non indat2d.f ; make -f pp2d.m ; mv pp2d ..
cd FEATFLOW/application/example.
This data file input files/indat2d.f non is almost identical to the preceeding one:
–
The inflow velocity profile at x = 0.0 is again parabolic, but with a maximum value
umax = 1.5.
–
UMEAN is changed for defining drag and lift.
Now, we can execute pp2d with given parameter file #data/pp2d.dat. A solution is
calculated on level NLMAX=4, with (stationary) start solution file #data/#DX4 stat. Now,
the chosen discretization scheme for the convective parts is the upwinding ansatz, and
we perform our time stepping with the fractional step scheme and fully adaptive time
step control until TIMEMX=4D0. Files for graphical output, namely for #avs and #byu, are
written out all 1 ”second”. As result we obtain the protocol file #data/pp2d.non which is
very similar to #data/cc2d.stat, but containing information for all time steps performed.
We explain this protocol file for only one macro time step.
Protocol file #data/pp2d.non:
..
.
ILEV,NVT,NMT,NEL,NVBD:
ILEV,NVT,NMT,NEL,NVBD:
ILEV,NVT,NMT,NEL,NVBD:
ILEV,NVT,NMT,NEL,NVBD:
1
2
3
4
time for grid generation :
ILEV,NU,NA,NB:
ILEV,NU,NA,NB:
ILEV,NU,NA,NB:
ILEV,NU,NA,NB:
1
2
3
4
ILEV,NP,NC:
ILEV,NP,NC:
ILEV,NP,NC:
ILEV,NP,NC:
148 706
592 2892
2368 11704
9472 47088
1
2
3
4
165 313 148 34
626 1218 592 68
2436 4804 2368 136
9608 19080 9472 272
3.6333331540227
313 2089 592
1218 8322 2368
4804 33220 9472
19080 132744 37888
..
.
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
MACRO STEP
68 AT TIME = 0.294D+01 WITH 1 STEP : DT1 = 0.288D-01
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
-------------------------------------------------------------------------------IT REL-U1
REL-U2
DEF-U1
DEF-U2
RHONL
OMEGNL
RHOMG1
RHOMG2
3.2 The 2D example
-------------------------------------------------------------------------------0
0.17D-03 0.23D-03
-------------------------------------------------------------------------------1 0.20D+00 0.26D+00 0.14D-04 0.77D-05 0.62D-01 0.10D+01 0.18D-01 0.13D-01
2 0.45D-01 0.28D-01 0.11D-05 0.14D-05 0.78D-01 0.10D+01 0.37D-01 0.57D-01
-------------------------------------------------------------------------------IT DIV-L2
RHOMGP
-------------------------------------------------------------------------------6 0.11D-12 0.37D-01
-------------------------------------------------------------------------------++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
1( 68) TIME= 0.297D+01 REL2(P)= 0.705D+00 RELM(P)= 0.221D+01
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
MACRO STEP
68 AT TIME = 0.294D+01 WITH 3 STEPS: DT3 = 0.960D-02
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
-------------------------------------------------------------------------------IT REL-U1
REL-U2
DEF-U1
DEF-U2
RHONL
OMEGNL
RHOMG1
RHOMG2
-------------------------------------------------------------------------------0
0.49D-04 0.67D-04
-------------------------------------------------------------------------------1 0.67D-01 0.78D-01 0.15D-05 0.73D-06 0.23D-01 0.10D+01 0.51D-02 0.40D-02
-------------------------------------------------------------------------------IT DIV-L2
RHOMGP
-------------------------------------------------------------------------------5 0.39D-12 0.35D-01
-------------------------------------------------------------------------------++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
1( 68) TIME= 0.295D+01 REL2(P)= 0.151D+01 RELM(P)= 0.581D+01
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-------------------------------------------------------------------------------IT REL-U1
REL-U2
DEF-U1
DEF-U2
RHONL
OMEGNL
RHOMG1
RHOMG2
-------------------------------------------------------------------------------0
0.69D-04 0.94D-04
-------------------------------------------------------------------------------1 0.93D-01 0.12D+00 0.21D-05 0.10D-05 0.23D-01 0.10D+01 0.52D-02 0.41D-02
2 0.72D-02 0.62D-02 0.66D-07 0.80D-07 0.29D-01 0.10D+01 0.13D-01 0.26D-01
-------------------------------------------------------------------------------IT DIV-L2
RHOMGP
-------------------------------------------------------------------------------5 0.31D-12 0.33D-01
-------------------------------------------------------------------------------++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
2( 68) TIME= 0.296D+01 REL2(P)= 0.153D+01 RELM(P)= 0.537D+01
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-------------------------------------------------------------------------------IT REL-U1
REL-U2
DEF-U1
DEF-U2
RHONL
OMEGNL
RHOMG1
RHOMG2
-------------------------------------------------------------------------------0
0.50D-04 0.66D-04
-------------------------------------------------------------------------------1 0.68D-01 0.83D-01 0.15D-05 0.68D-06 0.22D-01 0.10D+01 0.58D-02 0.45D-02
-------------------------------------------------------------------------------IT DIV-L2
RHOMGP
-------------------------------------------------------------------------------5 0.31D-12 0.35D-01
-------------------------------------------------------------------------------++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
3( 68) TIME= 0.297D+01 REL2(P)= 0.279D+01 RELM(P)= 0.961D+01
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
OLD DT= 0.96D-02 U(L2)= 0.13D-02 U(MX)= 0.51D-02 P(L2)= 0.45D-02 P(MX)= 0.95D-02
CHOICE 0( 68)
---- NEW DT = 0.954D-02 -- OLD DT = 0.960D-02
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
P(VELO) -0.12490D-01 0.45097D+00 0.10542D+01-0.67446D+00
P(PRES) 0.17141D+01-0.41863D+00-0.64658D+00 0.17289D+00
57
58
Examples for the use of featflow
I(PRES) 0.30836D+00-0.41448D-02
I(FORCE) 0.29159D+01-0.54744D+00
P(FLUX) 0.20401D+00
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
# 68 ( 272) TIME= 0.297D+01 RELU(L2)= 0.25D+01 RELP(L2)= 0.20D+01 REL= 0.25D+01
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
total
mavec
konv.
bdry
LC
ILU
U-mg
P-mg
time
time
time
time
time
time
time
time
:
:
:
:
:
:
:
:
3364.2498697899
309.36813378334
406.31425285339
2.0668334960938
87.866734504700
323.13556003571
1063.4651870728
1141.2314567566
..
.
STATISTICS :
NWORK :
1300000
IWORKG:
116168
IWMAXG:
119323
IWORKI:
967903
IWMAXI:
967903
IWORK :
1039351
IWMAX :
1203452
-------------------------------------------------------------------------------total time :
5327.6336588524
appr. time :
5320.4342214540
grid time :
3.6333331540227
post time :
38.716011047363
lin. time :
1797.2014658451
-> mavec time :
487.26168847084
-> konv. time :
647.08451652527
-> bdry time :
3.0346069335938
-> LC
time :
139.76541614532
-> ILU
time :
520.05523777008
U-mg time :
1709.3506851196
P-mg time :
1771.5327262878
#substeps :
428
#mg P
:
2326
#nonlinear :
604
#mg U
:
1220
-------------------------------------------------------------------------------U-MULTIGRID COMPONENTS [in percent]:
smoothing
:
42.218687909096
solver
:
2.2564278443300
defect calc. :
32.417565448349
prolongation :
20.300764952555
restriction
:
2.6163307795102
P-MULTIGRID COMPONENTS [in percent]:
smoothing
:
50.378172496605
solver
:
4.9000111811739
defect calc. :
16.941610100891
prolongation :
19.692759428133
restriction
:
7.8014891041156
--------------------------------------------------------------------------------
3.2 The 2D example
59
Some explanations:
–
ILEV,NVT,NMT,NEL,NVBD, ILEV,NU,NA,NB as before
–
ILEV,NP,NC means number of pressure unknowns and nonzero matrix entries for
pressure matrix
–
We perform macro time step 68, at absolute time T = 2.94s, and the actually chosen
time step is t = 0.0096. As predictor step we perform a calculation with a three
times larger time step DT1=0.288D-01, and then three substeps with the actual
time step size, DT3=0.960D-02. Anyway, we perform first the nonlinear iteration for
the Burgers–step, which multigrid convergence rates RHOMG1 and RHOMG2 for each
convection–diffusion equation. Then, the corresponding pressure equation is solved,
with multigrid convergence rate RHOMGP.
–
The estimated time error in different norms (L2 and L∞ for velocity and pressure)
are written out, followed by the new adaptively chosen time step value.
–
RELU(L2)=0.25D+01 and RELP(L2)=0.20D+01 are measures for the time derivative
of the velocity, resp., the pressure.
Additionally, the point values for P(VELO), . . ., P(FLUX) are printed separately into corresponding files in the subdirectory #points and can be visualized by GNUPLOT, for
instance. Furthermore, we obtain other files for graphical output, namely in #avs and in
#byu, written with a ”1 second” delay. In the following picture, the corresponding pressure
plot for T = 4 (file #avs/u91.inp) is shown.
1
ca-value
0.5
0
-0.5
-1
0
0.5
1
1.5
2
time
2.5
3
3.5
4
Figure 3.6: Lift values for the nonstationary 2D calculation
Figure 3.7: Pressure isolines for the nonstationary 2D calculation for T = 4
60
Examples for the use of featflow
3.3. The 3D example
In this example we want to perform similar calculations in the following 3D domain:
Outflow plane
U=V=W=0
2.5m
U=V=W=0
0.16m
(0,H,0)
1.95m
0.1m
0.1m
0.41m
y
0.45m
0.15m
x
U=V=W=0
0.41m
z
Inflow plane
(0,0,0)
(0,0,H)
Figure 3.8: 3D domain: channel with cylinder
A corresponding coarse mesh, which is simply a 3D extension of our c2d2.tri grid, can
be obtained by using tr2to3. First of all, do the following:
cd FEATFLOW/application/example/input files
make -f tr2to3.m ; mv tr2to3 ..
make -f trigen3d.m ; mv trigen3d ..
cd FEATFLOW/application/example.
Next, the corresponding 3D coarse mesh #adc/c3d2.tri (see #data/tr2to3.dat) is generated by executing tr2to3. Here, in z–direction, 8 non–equidistant layers are defined.
The following avs picture can be generated by performing trigen3d.
Now, we are ready to start our 3D Navier–Stokes calculation. Again, let us begin with the
stationary example, for Reynolds number 20. In this case, we use the data file input files
/indat3d.f stat, and the following, already well known procedure has to be performed
cd FEATFLOW/application/example/input files
cp indat3d.f stat indat3d.f ; make -f cc3d.m ; mv cc3d ..
cd FEATFLOW/application/example.
3.3 The 3D example
61
Figure 3.9: Coarse mesh c3d2.tri
This data file input files/indat3d.f stat contains the following definitions:
–
The inflow velocity profile at x = 0.0 is parabolic, with a maximum value umax =
0.45, but leading to the same mean velocity as in the 2D case.
–
There is one boundary part containing natural boundary conditions (outflow !). This
is the boundary segment with x = 2.5, where the corresponding parameter IFLAG is
activated in subroutine NEUDAT.
–
The integral mean pressure is calculated over the full and half of the circular cylinder.
–
We define similar parameters for the calculation of drag and lift as in 2D, now over
the full cylinder surface. RHO is a density parameter, DIST a typical length scale,
and UMEAN is a ”mean velocity”.
–
We show the velocity values at the mesh points KPU(1)=3421 (corresponds to the
coordinates (0.65, 0.20, 0.205)) and KPU(2)=677 (= (0.85, 0.20, 0.205)), and the pressure values at KPP(1)=687 (= (0.45, 0.20, 0.205)), KPP(2)=690 (= (0.55, 0.20, 0.205)),
KPP(3)=3498 (= (0.50, 0.25, 0.205)) and KPP(4)=677 (= (0.85, 0.20, 0.205)). The difference of the pressure values for KPP(1)=687 and KPP(2)=690 can be used again for
determing a typical pressure difference on the cylinder.
Now, we can execute cc3d with given parameter file #data/cc3d.dat. A solution is
calculated on level NLMAX=3, and the corresponding solution vector is saved as (unformatted) file #data/#DX3 stat. The chosen discretization scheme for the convective parts
is the streamline–diffusion ansatz, and the stopping criterions are 1 · 10−2 for the maximum of relative changes. It is remarkable, that the number of smoothing steps NSM=64 is
surprisingly large, due to the chosen anisotropic mesh in combination with the streamline–
diffusion method. As result of cc3d we obtain the protocol file #data/cc3d.stat which
is almost identical to the 2D–version #data/cc2d.stat.
62
Examples for the use of featflow
Protocol file #data/cc3d.stat:
-------------------------------------------------------------------------------INPUT DATA
-------------------------------------------------------------------------------Parametrization file =
#pre/c3d.prm
Coarse grid file
=
#adc/c3d2.tri
Integer parameters of /IPARM/,etc. :
-------------------------------------------------------------------------------minimum mg-level: NLMIN =
1
maximum mg-level: NLMAX =
3
element type
=
3
Stokes calculation: ISTOK =
0
RHS
generation
=
0
Boundary generation =
0
Error evaluation
=
0
mass evaluation
=
0
lumped mass eval.
=
0
convective part
=
0
Accuracy for ST
=
4
Accuracy for B
=
3
ICUB mass matrix
=
7
ICUB diff. matrix
=
7
ICUB conv. matrix
=
7
ICUB matrices B1,B2
=
7
ICUB right hand side
=
1
minimum of nonlinear iterations: INLMIN =
1
maximum of nonlinear iterations: INLMAX =
20
type of mg-cycle: ICYCLE =
0
minimum of linear mg steps : ILMIN =
1
maximum of linear mg steps : ILMAX =
5
type of interpolation: IINT =
2
type of smoother : ISM =
1
type of solver : ISL =
1
number of smoothing steps : NSM =
64
number of solver steps : NSL =
500
factor sm. steps on coarser lev.:NSMFAC= 8
KPRSM,KPOSM ON LEVEL:
1 4096 4096
KPRSM,KPOSM ON LEVEL:
2 512 512
KPRSM,KPOSM ON LEVEL:
3 64 64
KPRSM,KPOSM ON LEVEL:
4 64 64
KPRSM,KPOSM ON LEVEL:
5 64 64
KPRSM,KPOSM ON LEVEL:
6 64 64
KPRSM,KPOSM ON LEVEL:
7 64 64
KPRSM,KPOSM ON LEVEL:
8 64 64
KPRSM,KPOSM ON LEVEL:
9 64 64
Real parameters of /RPARM/,etc. :
-------------------------------------------------------------------------------Viscosity parameter: 1/NU =
1000.00000000000
parameter for Samarskij-upwind: UPSAM =
1.0000000000000
lower limit for optimal OMEGA: OMGMIN =
0.
upper limit for optimal OMEGA: OMGMAX =
1.0000000000000
start value for optimal OMEGA: OMGINI =
1.0000000000000
limit for U-defects
:
EPSD
=
1.0000000000000D-05
limit for DIV-defects :
EPSDIV =
1.0000000000000D-05
limit for U-changes :
EPSUR =
5.0000000000000D-02
limit for P-changes :
EPSPR =
5.0000000000000D-02
defect improvement :
DMPD
=
1.0000000000000D-01
damping of MG residuals
: DMPMG =
0.50000000000000
limit for MG residuals
: EPSMG =
0.50000000000000
damping of residuals for solving: DMPSL =
1.0000000000000D-01
limit of changes for solving:
EPSSL =
1.0000000000000D-01
relaxation for the U-smoother: RLXSM =
1.0000000000000
relaxation for the U-solver : RLXSL =
0.80000000000000
lower limit optimal MG-ALPHA: AMINMG =
-10.0000000000000
upper limit optimal MG-ALPHA: AMAXMG =
10.0000000000000
Parameters of /NS.../ :
3.3 The 3D example
63
-------------------------------------------------------------------------------Time dependency
: ISTAT =
0
Number of time steps
: NITNS =
50
limit for time derivative: EPSNS =
1.0000000000000D-05
Total time
: TIMENS =
0.
Theta
: THETA =
1.0000000000000
Time step
: TSTEP =
1.0000000000000D-02
Fractional step
: IFRSTP =
1
Stepsize for nonsteady savings: INSAV =
0
Number of files
: INSAVN =
0
Time step for Film
: DTFILM =
0.
Time step for AVS
: DTAVS =
1.0000000000000
Time step for BYU
: DTBYU =
1.0000000000000
Level for velocity
: IFUSAV =
0
Level for pressure
: IFPSAV =
0
Level for streamlines
: IFXSAV =
0
Level for AVS
: IAVS =
3
Level for BYU
: IBYU =
2
Start file
: IFINIT =
1
Type of adaptivity
: IADTIM =
-1
Max. Time
: TIMEMX =
5.0000000000000
Min. Timestep
: DTMIN =
1.0000000000000D-06
Max. Timestep
: DTMAX =
9.0000000010000
Max. Timestep change
: DTFACT =
9.0000000010000
Time for start procedure
: TIMEIN =
1.0000000000000
EPS for start procedure
: EPSADI =
0.12500000000000
EPS for acceptance
: EPSADL =
1.2500000000000D-03
EPS for not acceptance
: EPSADU =
0.50000000000000
Acceptance criterion
: IEPSAD =
3
Start procedure
: IADIN =
2
Max.numbers of repetitions
: IREPIT =
1
-------------------------------------------------------------------------------ILEV,NVT,NAT,NEL,NEL0,NEL1,NEL2: 1 1485 3836 1184 1120 0 64
ILEV,NVT,NAT,NEL,NEL0,NEL1,NEL2: 2 10642 29552 9472 8864 96 512
ILEV,NVT,NAT,NEL,NEL0,NEL1,NEL2: 3 80388 231872 75776 70512 1168 4096
time for grid initialization :
ILEV,NU,NA,NB:
ILEV,NU,NA,NB:
ILEV,NU,NA,NB:
1
2
3
36.883334092796
3836 39356 7104
29552 313712 56832
231872 2505152 454656
time for initialization of linear operators :
62.316661834717
-------------------------------------------------------------------------------IT RELU
RELP
DEF-U
DEF-DIV DEF-TOT RHONL
OMEGNL
RHOMG
-------------------------------------------------------------------------------0
0.26D-03 0.68D-03 0.72D-03
-------------------------------------------------------------------------------1 0.10D+01 0.10D+01 0.93D-04 0.48D-05 0.93D-04 0.13D+00 0.10D+01 0.12D+00
2 0.51D+00 0.42D+00 0.14D-03 0.84D-06 0.14D-03 0.43D+00 0.10D+01 0.18D+00
3 0.26D+00 0.26D+00 0.28D-04 0.11D-06 0.28D-04 0.34D+00 0.10D+01 0.77D-01
4 0.10D+00 0.39D-01 0.34D-05 0.20D-06 0.34D-05 0.26D+00 0.99D+00 0.22D+00
5 0.40D-01 0.31D-01 0.29D-05 0.79D-07 0.29D-05 0.33D+00 0.99D+00 0.69D+00
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
1(
1) TIME= 0.000D+00 NORM(U)= 0.2160847D+00 NORM(P)= 0.7659621D-01
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
P(VELO)
P(PRES)
I(PRES)
I(FORCE)
0.32444D-01
0.19791D+00
0.61808D-01
0.63585D+01
STATISTICS :
NWORK :
IWORKG:
IWMAXG:
IWORKI:
0.12893D-02 0.21419D-04 0.20393D+00 0.16592D-02-0.28822D-04
0.36067D-01 0.19174D-02 0.45314D-01
0.95865D-01
0.38097D-02
10000000
1339516
3450964
9732825
64
Examples for the use of featflow
IWMAXI:
9732825
IWORK :
9732825
IWMAX :
9893601
-------------------------------------------------------------------------------total time :
21627.316276040
appr. time :
21626.999564104
grid time :
36.883334092796
post time :
51.515625000000
lin. time :
1290.7828254700
-> mavec time :
109.70226669312
-> konv. time :
1106.0506057739
-> bdry time :
0.75054168701172
-> LC
time :
74.279411315918
mg
time :
20247.817779541
#substeps :
1
#nonlinear :
5
#mg
:
6
-------------------------------------------------------------------------------MULTIGRID COMPONENTS [in percent]:
smoothing
:
98.178536576809
solver
:
0.46786725951263
defect calc. :
0.70117728389642
prolongation :
0.59852068627434
restriction
:
5.3898193507251D-02
--------------------------------------------------------------------------------
Only a few differences with respect to the 2D example have to be explained:
–
ILEV,NEL0,NEL1,NEL2 appears in combination with streamline–diffusion only. It
presents the number of element requiring fully trilinear/linear/axiparallel transformations on the reference element.
–
P(VELO), . . ., I(IFORCE) are analogously defined as in the 2D case.
As before, we obtain files for graphical output, namely in #avs the file u1.inp, and in
#byu the byu files u1.vec (velocity) and p1.scl (pressure). A typical example for a avs
velocity plot in the midplane is shown in the following picture.
Figure 3.10: Velocity plot in the midplane z = 0.205 for the stationary 3D calculation
3.3 The 3D example
65
Next, we perform a nonstationary calculation, again for Reynolds number 100. We use
the data file input files/indat3d.f non, and the same ”well known” procedure has to
be performed:
cd FEATFLOW/application/example/input files
cp indat3d.f non indat3d.f ; make -f pp3d.m ; mv pp3d ..
cd FEATFLOW/application/example.
This data file input files/indat3d.f non is almost identical to the stationary one:
–
The inflow velocity profile at x = 0.0 is again parabolic, but with a maximum value
umax = 2.25.
–
UMEAN is changed for defining drag and lift.
Now, we can execute pp3d with given parameter file #data/pp3d.dat. A solution is calculated on level NLMAX=3, with (stationary) start solution file #data/#DX3 stat. Now, the
chosen discretization scheme for the convective parts is done by upwinding, and we perform our time stepping with the fractional step scheme and fully adaptive time step control
until TIMEMX=9D0. Files for graphical output, namely for #avs and #byu, are written out
all 1 ”second”. We obtain the protocol file #data/pp3d.non which is almost identical to
the preceeding protocol files. Therefore, we renounce a more precise explanation.
Again, the point values for P(VELO), . . ., I(FORCE) are printed separately into corresponding files in the subdirectory #points and be visualized by GNUPLOT, for instance (see
Figure 3.6). Furthermore, we obtain other files for graphical output, namely in #avs and
in #byu, written with a ”1 second” delay. In the following Figure 3.12, the corresponding
velocity plot for T = 9 (file #avs/u110.inp) is shown.
0.08
0.06
ca-value
0.04
0.02
0
-0.02
-0.04
-0.06
-0.08
0
1
2
3
4
5
6
7
8
9
time
Figure 3.11: Lift values for the nonstationary 3D calculation
We hope, that the presented examples are helpful in understanding the features of featflow, and give a first feeling how to use it. For more questions, the author is hopefully
prepared to your problems.
66
Examples for the use of featflow
Figure 3.12: Velocity plot in the midplane z = 0.205 for T = 9
Bibliography
[1]
Axelsson, O., Barker, V.A.: Finite Element Solution of Boundary Value Problems,
Academic Press, 1984
[2]
Becker, R., Rannacher, R.: Finite element discretization of the Stokes and Navier–
Stokes equations on anisotropic grids, Proc. 10th GAMM-Seminar, Kiel, January
14–16, 1994 (G. Wittum, W. Hackbusch, eds.), Vieweg
[3]
Blum, H., Harig, J., M¨
uller, S., Turek, S.: FEAT2D . Finite element analysis
tools. User Manual. Release 1.3, Technical report, University Heidelberg, 1992
[4]
Chorin, A.J.: Numerical solution of the Navier–Stokes equations, Math. Comp.,
22, 745–762 (1968)
[5]
Crouzeix, M., Raviart, P.A.: Conforming and non–conforming finite element methods for solving the stationary Stokes equations, R.A.I.R.O. R–3, 77–104 (1973)
[6]
Cuvelier, C., Segal, A., Steenhoven, A.: Finite element methods and Navier Stokes
equations, D. Reidel Publishing Company, Dordrecht 1986
[7]
Gresho, P.M.: On the theory of semi–implicit projection methods for viscous incompressible flow and its implementation via a finite element method that also
introduces a nearly consistent mass matrix, Part 1: Theory, Int. J. Numer. Meth.
Fluids, 11, 587–620 (1990). Part 2: Implementation, Int. J. Numer. Meth. Fluids,
11, 621–659 (1990)
[8]
Girault, V., Raviart, P.A.: Finite Element Methods for Navier–Stokes Equations,
Springer, Berlin–Heidelberg 1986
[9]
Harig, J., Schreiber, P., Turek, S.: FEAT3D. Finite element analysis tools in 3
dimensions. User Manual. Release 1.1, Technical report, University Heidelberg,
1993
[10]
Heywood, J., Rannacher, R., Turek, S.: Artificial boundaries and flux and pressure
conditions for the incompressible Navier–Stokes equations, to appear in: Int. J.
Numer. Meth. Fluids
[11]
Kloucek, P., Rys, F.S.: On the stability of the fractional step–θ–scheme for the
Navier–Stokes equations, SIAM J. Numer. Anal., 31, 1312–1335 (1994)
[12]
M¨
uller, S., Prohl, A., Rannacher, R., Turek, S.: Implicit time–discretization of
the nonstationary incompressible Navier–Stokes equations, Proc. 10th GAMMSeminar, Kiel, January 14–16, 1994 (G. Wittum, W. Hackbusch, eds.), Vieweg
67
68
BIBLIOGRAPHY
[13]
Rannacher, R.: Numerical analysis of the Navier–Stokes equations, Appl. Math.,
38, 361–380 (1993)
[14]
Rannacher, R.: On Chorin’s projection method for the incompressible Navier–
Stokes equations, in ”Navier–Stokes Equations: Theory and Numerical Methods” (R. Rautmann, et al., eds.), Proc. Oberwolfach Conf., August 19–23, 1991,
Springer, 1992
[15]
Rannacher, R., Turek, S.: A simple nonconforming quadrilateral Stokes element,
Numer. Meth. Part. Diff. Equ., 8, 97–111 (1992)
[16]
Sch¨afer, M., Turek, S. (with support by F. Durst, E. Krause, R. Rannacher):
Benchmark computations of laminar flow around cylinder, in E.H. Hirschel (editor), Flow Simulation with High-Performance Computers II, Volume 52 of Notes
on Numerical Fluid Mechanics, 547–566, Vieweg, 1996.
[17]
Schieweck, F.: A parallel multigrid algorithm for solving the Navier–Stokes equation, Impact Comp. Sci. Engnrg., 5, 345–378 (1993)
[18]
Schreiber, P.: A new finite element solver for the nonstationary incompressible
Navier–Stokes equations in three dimensions, Thesis, University of Heidelberg,
1995
[19]
Schreiber, P., Turek, S.: An efficient finite element solver for the nonstationary
incompressible Navier–Stokes equations in two and three dimensions, Proc. Workshop ”Numerical Methods for the Navier–Stokes Equations”, Heidelberg, Oct. 25–
28, 1993, Vieweg
[20]
Thomasset, F.: Implementation of Finite Element Methods for Navier–Stokes
Equations, Springer, New York 1981
[21]
Turek, S.: A comparative study of time stepping techniques for the incompressible
Navier–Stokes equations: From fully implicit nonlinear schemes to semi–implicit
projection methods, Int. J. Numer. Meth. Fluids, 22, 987 – 1011 (1996)
[22]
Turek, S.: On discrete projection methods for the incompressible Navier–Stokes
equations: An algorithmical approach, Comput. Methods Appl. Mech. Engrg., 143,
271 – 288 (1997)
[23]
Turek, S.: Tools for simulating nonstationary incompressible flow via discretely
divergence–free finite element models, Int. J. Numer. Meth. Fluids, 18, 71–105
(1994)
[24]
Turek, S.: Multilevel Pressure Schur Complement techniques for the numerical
solution of the incompressible Navier–Stokes equations, Habilitation Thesis, University of Heidelberg, 1997
[25]
Turek, S.: Multigrid techniques for a divergence–freefinite element discretization,
East-West J. Numer. Math., Vol. 2, No. 3, 229–255 (1994)
[26]
Turek, S.: Efficient solvers for incompressible flow problems: An algorithmic approach in view of computational aspects, Springer, 1998
[27]
Van Kan, J.: A second–order accurate pressure–correction scheme for viscous incompressible flow, SIAM J. Sci. Stat. Comp., 7, 870–891 (1986)
A. Appendix: Troubleshooting with featflow
In the following we list some problems which may apparently occur (as far as the author
knows) during the installation and execution of featflow. The following list cannot contain everything, and the author will be very grateful for showing him even more problems
(but, hopefully, with solution strategies).
1) Known problems and errors during installation:
Q.:
A.:
I have no gunzip to decompress my featflow binary file?
The author may tell you how to get gzip/gunzip or can send you the featflow
data in another format!
Q.:
A.:
I have no FORTRAN77 compiler?
The author may tell you how to get an ”older”, but ”free” FORTRAN77 or the
GNU FORTRAN77 compiler!
Q.:
A.:
I get a tar error after decompression?
Your gzip/gunzip command may not work well, or you did forget to activate the
binary–mode during ftp transfer!
I have problems with the execution of the make shell scripts while installing the
system libraries?
A.: Be sure that you start the installation in a /bin/csh C–shell!
Q.:
Q.:
A.:
The installation stops directly after calling make lib?
Be sure that you the directories FEATFLOW/object/libraries/libgen exist!
Q.:
A.:
I get errors during the compilation process with make lib?
Check your compiler options used, perhaps together with your system administrator!
Q.:
A.:
My application stops with an error immediately after starting?
Check your ztime.f subroutine in FEATFLOW/source/feat2d/src!
Q.:
If I start the makefiles for cc2d or cc3d in my application, I get the error message
that the VANCA routines are not linked?
A.: This happens on some computers. Copy the file FEATFLOW/source/cc2d/src/
vanca.f to FEATFLOW/source/cc2d/mg. Then, start there again the shell–script
cc2d mg.m. Do the analogous procedure for the 3D case!
69
70
Appendix: Troubleshooting with featflow
2) Known problems and errors during preprocessing:
Q.:
A.:
Q.:
A.:
Q.:
A.:
I have problems with omega2d?
Check your PC–emulation and that you have (at least) WINDOWS 3.1. If so, ask
the author for more advice!
I have problems with the german manual of omega2d?
Sorry, but in summer 1998 our first version of the DEVISOR will be finished...
I have problems with trigen2d?
Check in your manual of the actual release, that you marked an admissible set of
elements for your desired adaptive refinement strategy!
3) Known problems and errors during solution process:
Q.:
A.:
Q.:
A.:
I have problems while seemingly the triangulations on every level are created?
The most usual error is that your parametrization or coarse mesh is wrong. For
instance, both parametrization curves follow the same direction (error!!!), or the
starting point of your parametrization is not captured by a mesh point. Additionally,
be sure that all mesh points, which were created while using omega2d, do really
belong to an element. In most cases, one of these errors causes your trouble! Another
reason might be, if you have more than 1 boundary component, that you must have
positioned at least 3 mesh points on every boundary component.
I have problems with reading an unformatted solution file?
This happens on some computers. Use formatted output!
Q.:
The code starts correctly, but before finishing the first iteration step, it stops with
an error?
A.: Be sure that the corresponding provided memory size NNWORK is large enough!
Q.:
A.:
The code starts correctly, but the multigrid rates become (almost) identical 1?
Perhaps your mesh is so bad! But in most cases, you prescribed an inflow profile
only, and your definition of the boundary part containing natural b.c.’s is wrong
(so, your flow cannot get incompressible!!!). Check this in the data file indat2d.f
(analogously in 3D)!
Q.:
My solution schemes (linear muligrid, nonlinear schemes, time stepping schemes)
are crashing?
A.: Perhaps your problem is so bad! If you perform a stationary calculation, try the
same with the nonstationary code (with adaptively chosen time step size). If you still
have problems, then your mesh might be too hard (too large aspect ratios!!!). Check
your triangulation. Take the parameter files belonging to the most robust version.
However, in most cases it is sufficient to check again the data and the parameter file!
4) Known problems and errors during postprocessing:
Q.:
A.:
I have none of the proposed graphic tools?
That’s a really hard problem! Try to find another one which has similar features, or
ask your local computer center, or ask the author . . . for the (hopefully soon) freely
available AVS EXPRESS modules!
B. Appendix: The featflow group
We are indebted to the following persons who were involved, with theoretical and practical
help or suggestions, to develop featflow. We hope they won’t be too angry about arising
questions because of some misleading comments in this manual. However, this is still the
first version, and source code and man pages of featflow will grow hopefully.
Whoever is interested in getting featflow, or whoever has questions or trouble, please,
send an email to S.Turek/Chr.Becker or to
[email protected]
List of involved persons:
Christian Becker
Roland Becker
Heribert Blum
Phil Gresho
Joachim Harig
Jaroslav Hron
John Heywood
Susanne Kilian
Steffen M¨
uller–Urbaniak
Hubertus Oswald
Rolf Rannacher
Ludmilla Rivkind
Friedhelm Schieweck
University of Heidelberg
University of Heidelberg
University of Dortmund
LLNL Livermore
University of Heidelberg
Charles University of Prague
UBC Vancouver
University of Heidelberg
FORD K¨
oln
University of Heidelberg
University of Heidelberg
University of Heidelberg
University of Magdeburg
Rainer Schmachtel
Peter Schreiber
Stefan Turek
John Wallis
Owen Walsh
University of Heidelberg
University of Heidelberg
University of Heidelberg
University of Heidelberg
UBC Vancouver
71
[email protected]
[email protected]
[email protected]
[email protected]
[email protected]
[email protected]
[email protected]
[email protected]
[email protected]
[email protected]
Friedhelm.Schieweck@Mathematik.
Uni-Magdeburg.de
[email protected]
[email protected]
[email protected]
[email protected]
[email protected]
C. Appendix: Future projects in featflow
In the following we list the current projects which are under development, and which shall
be added to one of the next releases of featflow.
Concerning the improvement of the discretization and solution schemes we are just (or
will soon begin with) developing software containing:
–
The ”mixed coupled” solvers cp2d and cp3d
–
Periodic and moving boundaries
–
A Galerkin method with adaptive error control in space
–
A space–time Galerkin method with adaptive error control
–
A parallel version of all codes (for shared and distributed memory machines)
–
An improved feat version (feast !)
Software for a improved preprocessing:
–
omega2d as CAD–tool under JAVA (see the DEVISOR !)
–
omega3d as CAD–tool under JAVA (see the DEVISOR !)
–
trigen2d with mixing of refinement types
–
trigen3d with adaptive refinement strategies
–
intpol2d and intpol3d for interpolation between arbitrary meshes
Software for a improved postprocessing:
–
”Free” modules of AVS EXPRESS
–
Cooperation with other graphic packages
72
Appendix: Future projects in featflow
Software for more complex models:
–
Adding and testing of Boussinesq- and more complex turbulence models
–
Weakly compressible flows
–
Non–newtonian and multiphase flows
–
Shape optimization (with respect to lift and drag, for instance)
And finally, for the ”users”:
–
Other FORTRAN compilers (GNU, for instance)
–
A FORTRAN90 version (see feast !)
–
A C or C++ version (?)
–
Industrial and commercial applications !!!
If you have more suggestions, please, let them be known the authors.
73