User:Egm6936.s11/FENICS

 Could you write a doc for "dummies" for installing and using fenics and comment on the knowledge requirements (e.g., FEM, C++ or python) ? Like you said, the knowledge of C++ or python is required, in addition to FEM. Keep in mind that students would know matlab, but not necessarily C++ or python; so the minimum knowledge of computer language should be spelled out. Egm6321.f11 14:33, 23 June 2011 (UTC)

Here is a documentation of how to install fenics and build up from the Poisson example problem to an incompressible Navier-Stokes code.

It is assumed the reader has some knowledge of C++, linux, and variational calculus (to obtain the weak form of a problem).

 Write as follows: Here is a documentation of how to install and use FENiCS. Egm6321.f11 21:07, 28 June 2011 (UTC)


 * fenics C++ Demos: Good for working examples
 * fenics-book: Perhaps the most useful reference.
 * Computational Turbulent Incompressible Flow: Reference for CFD

The work that follows was performed on Ubuntu 11.04 (64 bit) within a VirtualBox machine. The environment was chosen for developing code with the fenics project as Ubuntu packages are maintained directly by the fenics project team.

Installation
Installation notes for installing from Ubuntu's debian package manager, and from source can be found on Fenics Installation

Compiling Examples
The installation contains several demos for the fenics project in python and C++. The demos are installed into the path /usr/share/dolfin/demo

To compile and run the C++ demos


 * Copy the demo out of the usr/share/dolfin/demo directory to working directory such as home.

For example, to copy the Poisson demo into your home directory


 * Compile using cmake to generate the make file, then calling make itself.

The demos use cmake to configure the C++ demos, for more on cmake see a tutorial such as: CMake tutorial

CMakeBuilder
There is an eclipse plugin for working with CMake projects.

http://www.cmakebuilder.com/

Poisson's Equation Tutorial in C++
Poisson's Equation in C++ is a great first coding attempt as the full C++ code exists in the demo folder for comparison and the test case is fully documented in the fenics book linked at the top (though in Python). This first project is for understanding the details of the code and process for developing a fenics C++ code.

First things first, the CMakeLists.txt
The core of the project is the CMakeLists.txt file. This file defines the "project" in cmake terms. The file describes which libraries to include, the pathing, the output executable and which files are its source files. The particulars of using fenics (DOLPHIN) are handled with the demo's CMakeLists.txt files, so we will start with copying one of their files.

In this file there are three commands/lines which need to be changed for a project. The first, and perhaps most obvious is the project name near the top. The other two are the last two commands listed below.

The second to last line defines the executable name and its source code (main.cpp). The last line links the binary to the libraries. In these last two lines you will need to add any additional source files and change the output binary name as desired.

Generating an IDE project
CMake can generate project files for a number of IDEs (integrated development environments). In this work we use Eclipse CDT4, but the choice is a matter of preference.

To generate a project for Eclipse CDT4,

This is a particularly convienient way for obtaining a project for a demo. To open the project in Eclipse CDT4, start Eclipse, file, import, general > Existing Projects into Workspace. Browse to the path where cmake was called and finish to import the project.

Note that I had to manually update the includes to get the indexer up to date.

To build a project from command line and not within an IDE,

The binary can then be executed from within the directory. To clean a directory with make,

The uniform format language file
For a 'static' language to achieve the dynamic ability code generation is necessary. This generation is accomplished in the .ufl (uniform format language) file extensions. These files generate C++ header files to be included by the C++ code file. The generated header file is synonymous with the ufl file. It's just been compiled (interpreted). Note that python, being a dynamic scripting language, voids this step as it is performed at runtime.

The .ufl file from the demo for the Poisson problem

 Explain $$\displaystyle dx$$ and $$\displaystyle ds$$ that these are fixed keywords, i.e., they are part of the fenics language. Egm6321.f11 14:35, 30 June 2011 (UTC)

This defines an element with triangular elements using Lagrange interpolation of first order. If you were to run the problem in one-dimension you could use an "interval" instead of "triangle".

The file then defines the trial and test functions, u and v respectively, for that element along with coefficients f and g (not these are not constant in space).

The next step is to generate Poisson.h from Poisson.ufl, with the following command

We now include the header file into a basic C++ program. Note that the details of Poisson.h can be blackboxed for now.

Poisson C++ code
Now for the interesting part of the code, the main C++ file. The first step is to include the necessary header files and include the dolfin namespace for convenience.

Now we also need to define our source terms, normal derivatives at the boundary, and the Dirichlet boundary conditions. These are defined as classes that inherit from dolfin classes and specify the functional behavior of these terms (e.g., as of space). In C++ these classes must be define before they are used, so we place the three before the main function.

The main function
The entry point for the program is the function main, which returns zero at the end as is standard. We will look at main in two parts, setting up the problem, and solving/plotting the problem. The code is discussed through comments in code

Running the code
Compile the code, and run it to view the solution which pops up in 3D

Poisson's problem in one-dimension
Most, if not all of the examples are two-dimensions. Running in one-dimension, the same problem as listed above involves:


 * Using "interval" instead of "triangle" in the ufl definition of the element
 * Updating the Source class to only reference x[0], and not x[1], the second dimension which no longer exists.
 * Defining a one-dimensional mesh with something like: Interval mesh(20, -1, 1);

Other than that, no change is necessary. The resulting plot is automatically reduce to a 2-d scatter.

Newton's law of cooling on a bi-unit square
This problem is borrowed from EML5526 Spring 2011 at the University of Florida as taught by Prof. Loc Vu-Quoc, with wiki presence here: Eml5526.s11

The problem has several solutions from students in the course:


 * Team 1
 * Team 2
 * Team 3
 * Team 4
 * Team 5
 * Team 6
 * Team 7

Problem statement


Consider Newton's Law of Cooling. The PDE for the heat problem is given by:

Solution
We note that the time derivative and source term are zero. Further $$k=1$$, so that this is a specific case of the Poission equation, the Laplace equation.

The additional challenge is then to define the direchlet, neuman and mixed (Robin) boundary conditions.

Setting up the variational problem
The first step is to account for the additional term from the mixed boundary. Note from one of the solutions the weak form from team 4 is


 * {| style="width:100%" border="0"

$$  \displaystyle \begin{align} \left(w H\left(u - u_\infty\right) \right|_{x=-1}      +       \left( w h \right|_{y=1} - \int_\Omega\left( \frac{\partial w}{\partial x}\frac{\partial u}{\partial x}       +       \frac{\partial w}{\partial y}\frac{\partial u}{\partial y} \right) = 0  \end{align} $$     (7.4.12)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

Notation aside, there are two things that may come to attention. First is the varying location of the boundaries $$x=-1$$, and $$y=1$$. The second thing to note is in the mixed boundary, where we have a term that is a function of $$u$$ and another that is a constant $$u_\infty$$.

To handle these two problems we update the ufl file as follows

where we have 'marked' the boundaries with numerical indexes with $$ds(0), ds(1)$$. Note that the mixed boundary has been split between a and L. This splitting is the key step to getting a mixed boundary condition to work.

The main.cpp file is mostly the same with the addition of handling the boundary condition parts and the new constants. These changes are perhaps best seen with a code listing,

The solution
The result which is plotted to screen is shown below. Note pressing 'h' within the plot window will list the shortcuts to console. Note the solution is the same as those of the teams.



Solution
This problem introduces non-linearity in the convective term. With the non-linear term the solution method becomes iterative and we must linearize the problem in some manner. There are several approaches, four of which are discussed in the fenics book. We will linearize the PDE directly.

To linearize, we first approximate $$u$$ as $$\displaystyle u^{k+1}=u^k+\zeta$$, and substitute into our PDE. Then we multiply by a test function, and perform our partial integration to obtain the variational form. The resulting variational formulation is

The variation problem is then to find $$\displaystyle \zeta$$ such that $$\displaystyle a(\zeta, v) = L(v)$$, $$\displaystyle \forall v$$.

Note the separation in the left hand and right hand sides by

 In mechanics, often the notation $$\displaystyle \delta u$$ is used to denote a virtual quantity (displacement or velocity), which is equivalent to the weighting (or test) function, and which would be eliminated upon approximation to yield the algebraic equation. Here, you are using $$\displaystyle \delta u$$ to denote a perturbation or increment from the known solution at a previous time step. It is possible to use $$\displaystyle \Delta u$$, but $$\displaystyle \Delta$$ may be misunderstood with the Laplacian. Maybe $$\displaystyle w$$ can be used to designate the unknown function (increment) to be found. Egm6321.f11 13:45, 30 June 2011 (UTC) I chose $$\displaystyle \zeta$$ to avoid confusion Rodbourn 21:35, 13 July 2011 (UTC)

 The complete statement of the weak form (to avoid confusion with the notation, e.g., $$\displaystyle \delta u$$ misunderstood as virtual velocity), you can write as follows: Find $$\displaystyle \delta u$$ such that $$\displaystyle a(\delta u, v) = L(v)$$, $$\displaystyle \forall v$$. Egm6321.f11 13:45, 30 June 2011 (UTC)

Key in this is the assumption that $$u^k$$ is a coefficient, known at the time of solving. Then we continually update our approximation by iterating through our k index, correcting our solution by $$\zeta$$.

First our UFL file

Then our main C++ program

Convergence history, with an initial condition of u(x) = 0.0



Burger's Equation
Here we work with Burger's equation, without any terms removed. Building off of the previous code we have to handle the time derivative and integrate in time with some discritization.

Solution
Building upon the previous solution, we need to discritize in time as well as in space. For a simple discritization we shall discritize with

Here, $$u^{n-1}$$ is known and is a coefficient in the terms of fenics. The 'trick' is that this step is before the linearization of the PDE performed above, so that $$u^n$$ is further approximated by $$\scriptstyle u_{k+1}^n = u_k^n+\delta u^n$$, where k is the sub-iteration index used to converge to the non-linear solution (quadratic convergence).

The solution procedure is then to step in time and then iterate at each step in time to converge for that point in time.



$$\theta$$-implicit integration
We now update the time discritization to a $$\theta$$-implicit scheme.

Now, with $$\theta=0$$ we have an explicit integration scheme and with $$\theta=1$$ we have a fully implicit scheme. Values of $$\theta$$ in between are mixed methods.

Comments
One should not the stability of the $$\theta$$-implicit scheme. Fullly explicit integration requires 'small' timesteps while larger time steps may be taken with an implicit integration.

For example a time step size of 0.1 and $$\theta$$ value of 0.1 results in aliasing. The solution at $$t=1$$ s is shown below with the aliasing. Simply changing $$\theta$$ to a value of 0.6 results in convergence to the steady state solution.



One-Dimensional Incompressible Navier-Stokes
At this point we can solve the non-linear Burger's equation. An additional level of complexity can be added by solving the incompressible Navier-Stokes equations.


 * {| style="width:100%" border="0"

$$  \begin{align} \nabla \cdot {} & = 0 \\ u_t+u\nabla u {} & = -\nabla p + \frac{1}{Re}\nabla^2 u  \end{align} $$
 * style="width:95%;" |
 * style="width:95%;" |
 * 
 * }

The trick is to split the time derivative operator in a process known as operator splitting.


 * {| style="width:100%" border="0"

$$  \begin{align} \frac{u^*-u^n}{dt}+u\nabla u {} & = \frac{1}{Re}\nabla^2 u \\ \frac{u^{n+1}-u^*}{dt} {} & =  -\nabla p    \end{align} $$
 * style="width:95%;" |
 * style="width:95%;" |
 * 
 * }

Where $$\scriptstyle u^*$$ is not divergence free. Take the divergence of the pressure equation and use continuity to find the Poisson equation for $$u^*$$. Then we use the pressure equation to recover $$u^{n+1}$$


 * {| style="width:100%" border="0"

$$  \begin{align} \frac{u^*}{dt} {} & = - u\nabla u + \frac{1}{Re}\nabla^2 u + \frac{u^n}{dt} \\ \nabla \cdot u^* {} & = dt \nabla^2 p \\ u^{n+1} {} & =  - dt \nabla p + u^* \end{align} $$
 * style="width:95%;" |
 * style="width:95%;" |
 * 
 * }

The first equation is Burger's equation for $$\scriptstyle u^*$$, the second is a pressure-Poisson equation as we solved early on, and the last equation is a pressure correction step where we obtain a divergence free velocity field.

So our task is to add a pressure field and solve the two additional equations. The two additional equations, however, are linear.

Solution
We are limited to one variational problem per UFL file, so we now have three UFL files.

Our updated program

The solution is not all to interesting as the only divergence one dimensional flow is a constant velocity. It is, however, a nice transition into the 2D incompressible NS code.

Observations
One thing, which I didn't see document anywhere, is that one BilinearForm and one LinearForm is allowed per ufl file. If there are two terms, say a diffusive and convective term, sum them in the ufl file, don't try and create two BilinearForms and add them later.