MathematiceToFEniCS

Introduction
Installation
Examples
- 1D Laplace
- Oldroyd-B in Taylor-Couette flow

Introduction

The aim of this library is utilize symbolic power of Wolfram Mathematica and numeric power of FEniCS. The ideal work-flow is that you define symbolic form of your PDE in a weak form in Mathematica and then you call single function which generates FEniCS script which solves this PDE numerically.

Installation

Just simply download MathematicaToFEniCS with git by command git clone https://github.com/lecopivo/MathematicaToFEniCS

TODO Mention files which needs to be added to paths for Mathematica and python

Examples

1D Laplace

For a start we will solve simple system of two coupled Laplace equations

\begin{equation} \Delta u = 1-v \qquad \Delta v = u-1 \end{equation}

on the interval $[0,1]$ with boundary conditions $u(0)=u(1)=0$ and $v(0)=0,v(1)=0.5$.

Weak form of these equations is

\begin{equation} 0= \int_0^1 \nabla u \cdot \nabla \tilde u + (1-v) \tilde u \, dx \qquad 0= \int_0^1 \nabla v \cdot \nabla \tilde v + (u-1) \tilde v \, dx \end{equation}

where $\tilde u, \tilde v$ are test functions.

Note: We solve system of two equations because MathematicaToFEniCS for now does not support PDEs with only one unknown function.

Mathematica code

Open the Mathematica notebook in /tutorial/Laplace1D.nb. We will now go line by line of this code and explain what it does.

It starts with the line

<< (FileNameJoin[{NotebookDirectory[], "../MathematicaToFEniCS.m"}])

which loads the MathematicaToFEniCS library.

Next line says on which mesh the calculation should be preformed

mesh := "meshLine";

It has to be a name of a python script which defines the mesh object which has the same name i.e. meshLine in this example and also it has to define object boundary_parts which is facet function on this mesh marking different boundary parts. Example of such a file is

# file: meshLine.py
from dolfin import *

meshLine = IntervalMesh(100, 0, 1)

# Create boundary markers
boundary_parts = FacetFunction('size_t', meshLine)
left   = AutoSubDomain(lambda x: near(x[0], lineMin))
right  = AutoSubDomain(lambda x: near(x[0], lineMax))
left  .mark(boundary_parts, 1)
right .mark(boundary_parts, 2)

Next two lines defines variables names, unknown functions and their corresponding test functions.

vars := {x}; 
funs := { u, v};
testFuns := {tu, tv};

Now we define boundary conditions.

bcs := {
{{0, 1}, {0, 2}},
{{0, 1}, {0.5, 2}}
}

Symbol bcs has to be a list whose n-th element describes boundary conditions of n-th function and it has to be a list of pairs in the form {value,boundary-id}. So the first function $u$ has zero boundary condition on boundaries marked with 1 and 2, but the second function $v$ is zero on the boundary marked with 1 and 0.5 on the boundary marked with 2.

Next we specify which finite elements will be used. femSpace is list of pairs where n-th pair specifies what type of finite element should be used for n-th function

femSpaces = {{"Lagrange", 1}, {"Lagrange", 2}};

In this example we use Lagrange elements of order one for $u$ and of order 2 for $v$. For list of available finite element see FEniCS web.

The most important thing is to define weak forms of PDEs.

weak1 := Grad[u[x], vars].Grad[tu[x], vars] + (1 - v[x]) tu[x];
weak2 := Grad[v[x], vars].Grad[tv[x], vars] + (u[x] - 1) tv[x];
weakForm = {weak1, weak2};

The last line of code generates python script called Laplace1D.py

GenerateCode["Laplace1D", mesh, vars, funs, bcs, testFuns, femSpaces, weakForm];

So if you execute notebook Laplace1D.nb a python script Laplace1D.py should appear in the same directory and when executed you should see two graphs

TODO Mention Laplace1DExtra.py

Oldroyd-B in Taylor-Couette flow

Next example is Taylor-Couette flow of Oldroyd-B model.

The equation for stationary flow of Oldroyd-B are

\begin{align} (v \cdot \nabla) v &= \text{div}\, \mathbb{T}\\ \nabla \cdot v &= 0 \\ \mathbb{T} &= - p \mathbb{I} + 2 \mu \mathbb{D} + (\mathbb{A} - \mathbb{I}) \\ (v\cdot \nabla)\mathbb{A} - \mathbb{LA} - \mathbb{AL}^T &= - (\mathbb{A} - \mathbb{I}) \end{align}

Their weak form is

\begin{align} 0 &= \int (v \cdot \nabla) v \cdot \tilde v + \mathbb{T} : \nabla \tilde v \, dx \\ 0 &= \int v \cdot \nabla \tilde p \, dx \\ 0 &= \int \left( (v\cdot \nabla)\mathbb{A} - \mathbb{LA} - \mathbb{AL}^T + \mathbb{A} - \mathbb{I} \right) : \tilde A \, dx \end{align}

where $\tilde v,\tilde p, \tilde A$ are test functions.

Mathematica code

Code for Taylor-Couette flow of Oldroyd-B is /tutorial/TaylorCouetteFlowOldroydB.nb.

We first say that we want to calculate in cylindrical coordinates and the used variables are $\rho,\theta,z$, also we define jacobian for cylindrical coordinates.

coordType = "Cylindrical";
vars = {ρ, θ, z};
Jac = ρ;

Furthermore we define velocity V, pressure P, conformation tensor 𝔸 and specify their special forms in Taylor-Couette flow e.g. the only non-zero component of velocity is $θ$-component which can depend only on $\rho$.

V    = {0, vt[ρ] , 0};
TV   = {0, tvt[ρ], 0};
P    = p[ρ];
TP   = tp[ρ]; 
𝔸 =  {
       {Arr[ρ], Art[ρ], 0},
       {Art[ρ], Att[ρ], 0},
       {     0,      0, Azz[ρ]}
       };
t𝔸 = {
       {tArr[ρ], tArt[ρ], 0},
       {tArt[ρ], tAtt[ρ], 0},
       {      0,       0, tAzz[ρ]}
       };

Next we define other useful quantities such as velocity gradient or Cauchy stress tensor

𝕃 = Grad[V, vars, coordType];
𝔻 = 1/2 (𝕃+ Transpose@𝕃);
𝕀 = IdentityMatrix[3];
μ = 1;
𝕋 = -P 𝕀+  2 μ 𝔻+ (𝔸 - 𝕀);

Definitions of weak forms, weak1 is momentum equation, weak2 is continuity equation and weak3 is evolution equation for $𝔸$. Since continuity equation is automatically satisfied we do not include it to the final list of weak forms weakForm

weak1 = (𝕃 .V.TV + Tr[𝕋.Grad[TV, vars,coordType]]) Jac;
weak2 =  V.Grad[TP, vars,coordType] Jac;
weak3 = Tr[(Grad[𝔸, vars, coordType].V -𝕃. 𝔸 -𝔸.Transpose @𝕃 + (𝔸 - 𝕀)).t𝔸 ] Jac;
weakForm = {weak1, weak3};

Now we specify mesh used for computation, variable in which we do the numerical computation, unknown functions, boundary conditions - Dirichlet boundary conditions for velocity and no boundary conditions for components of 𝔸, test functions and as last we specify used finite elements.

mesh := "meshLine";
vars = {ρ};
funs = {vt, Arr, Art, Att, Azz};
bcs := {
	{{1, 1}, {0, 2}},
	{},
	{},
	{},
	{}
	};
testFuns := {tvt, tArr, tArt, tAtt, tAzz}; 
femSpaces = {{"Lagrange", 1}, {"Lagrange", 1}, {"Lagrange", 1}, {"Lagrange", 1},{"Lagrange", 1}};

The last line just generates the python script

GenerateCode["TaylorCouetteFlowOldroydB", mesh, vars, funs, bcs, testFuns, femSpaces, weakForm];

The generated python script should display these plots .

MathematiceToFEniCS

Table of Contents

Introduction

Installation

Examples

1D Laplace

Mathematica code

Oldroyd-B in Taylor-Couette flow

Mathematica code