[FreeFem++] Time dependent problem

(Revised: 2022/05/29)

In this webpage, I will introduce how time-dependent differential equations are solved in FreeFEM. All information are obtained in an official FreeFEM document. 

Our contents are as follows; 

1) Linear & time-dependent differential equations 
2) Nonlinear & time-dependent differential equations 

---

1. Time-dependent linear differential equations 

As an example of this problem, we will solve thermal transport. Differential equations, boundary conditions, and initial conditions for thermal transport are described as follows;

1) Differential equations 

\begin{equation} 
\partial_{t} Q = \nabla \cdot (\kappa \nabla (Q)) 
\end{equation}

※ Note that this differential equation is time dependent "parabolic" problem.
※ Note that In our case, we will assume that $\kappa$ is constant. 

---

2) Boundary conditions & initial states 

a. Boundary condition 
\begin{equation}
\kappa \vec{n} \cdot \nabla Q + \alpha (Q - Q_e) = 0 
\end{equation}
b. Initial states  
\begin{equation}
Q(x,y; t = 0) = Q_0 + x Q_1, (Q_0 = 10, Q_1 = 50)
\end{equation}

Based on this the weak form of our problem become as follows; 

\begin{equation}
\begin{split}
\int \frac{\partial}{\partial t} Q \psi - \kappa \nabla ^2 Q \psi dxdy &=& \int \frac{\partial}{\partial t} Q \psi + \kappa (\nabla Q) \cdot (\nabla \psi)  dxdy - \int (\vec{n} \cdot \nabla Q) \psi dl = 0 
\end{split}
\end{equation}

where $\psi$ is a test function. In order to deal with time difference, explict Euler finite difference approximation in time is used; 

\begin{equation}
\begin{split}
\int \frac{Q_{n} - Q_{n-1}}{\delta t} \psi + \kappa (\nabla Q_{n}) \cdot (\nabla \psi)  dxdy + \int \alpha (Q^n - Q_e) \psi dl 
 = 0 
\end{split}
\end{equation}

---

3) FreeFEM Code 

a. System Geometry 
For simplicity, we assume that our system has a simple geometry; 

The codes corresponding to this geometry are as follows; 

mesh Th = square(60, 30, [10 * x, 3 * y]);

 

b. Codes for weak formalism in FreeFem++ 

real dt        = 0.1;
real rdt       = 1. / dt;
real kappa     = 1;
real alpha     = 2.0;
real Qout      = 20;
func InitialCondition = 10.0 + x * 50;
macro grad(u) [dx(u), dy(u)] //

// Declaration for finite Element functions
// Applying Initial Conditkions
Vh Qn = InitialCondition, psi;
Vh Qold;  // Qnp: (n-1) th Heat distribution

problem HeatEquation(Qn, psi)
= int2d(Th)(  rdt * Qn * psi + kappa * grad(Qn) '* grad(psi))
+ int1d(Th)( alpha * (Qn) * psi)
- int1d(Th)( alpha * Qout * psi)
- int2d(Th)( rdt * Qold * psi)
+ on(2,4,  Qn = 10);

 

 

Appendix 1 

For comparison, I will show codes with explict matrix for weak form. First of all, let's construct the explict matrix form as follows; For discretized unknown function $Q$ and test functions, 

\begin{equation}
\begin{split}
Q \approx \sum_i N_i Q_i \\
\psi \approx \sum_i N_i \\
\end{split}
\end{equation} 

If we plug these functions into weak form, then 

\begin{equation}
\begin{split}
\sum_{i,j} \int  \frac{Q^{n}_{i} - Q^{n-1}_{i}}{\Delta t} N_i N_j + \kappa Q^{n}_{i}  (\nabla N_i) \cdot (\nabla N_j)  dxdy + \alpha (Q^{n}_{i})\int  N_{i} N_{j} dl  = 0
\end{split}
\end{equation}

\begin{equation}
\begin{split}
\sum_{i,j} Q^{n}_{i} \left(\int  \frac{N_i N_j}{\Delta t} + \kappa  (\nabla N_i) \cdot (\nabla N_j)  dxdy + \alpha \int  N_{i} N_{j} dl\right) - Q^{n-1}_{i} \int \frac{N_{i}N_{j}}{\Delta t} -\alpha Q_{out} \int N_{j} dl &=& 0 
\end{split}
\end{equation}

The above equations can be written in the matrix form;

\begin{eqnarray}
A \vec{Q}^{n} + M \vec{Q}^{n-1} - \vec{b} = 0 \\ 
\rightarrow \vec{Q}^{n} = A^{n-1}( - M \vec{Q}^{n-1} + \vec{b})
\end{eqnarray}

The $A$, $M$, and $\vec{b}$ are as belows; 

\begin{eqnarray}
A_{ij} &=&\int  \frac{N_i N_j}{\Delta t} + \kappa  (\nabla N_i) \cdot (\nabla N_j)  dxdy  + \alpha \int  N_{i} N_{j} dl \\ 
M_{ij} &=& \int \frac{N_{i}N_{j}}{\Delta t} dxdy \\ 
b_j &=& \alpha \int N_{j} dl
\end{eqnarray}

Now, based on the above equations, we can write our matrix in FreeFEM++ code. Basically, we need to write $A$, $M$, and $\vec{b}$. Let's write the corresponding ones; 

// Declaration for finite Element functions
// Applying Initial Condition
Vh Qn = InitialCondition;
Vh Qold = Qn;  // Qnp: (n-1) th Heat distribution
Vh psi;

// Declaration for explicit matrix form of the weak form
varf MatrixA(Qn, psi)
  = int2d(Th)(rdt * Qn * psi + kappa * grad(Qn) '* grad(psi))
  + int1d(Th,1,3)( alpha * Qn * psi)
  + on(2,4,  Qn = 10);
varf MatrixM(Qn, psi)
  = int2d(Th)( rdt * Qn * psi)
  + on(2,4,  Qn = 10);
varf RightPart(unused, psi)
  =  int1d(Th,1,3)( alpha * Qout * psi);
varf DirichletBC(unused, psi)
  = on(2,4, unused = 10);

// assignments for the matrix 
real tgv = tgv;
matrix A = MatrixA(Vh, Vh);
matrix M = MatrixM(Vh, Vh);
real[int] b0  = RightPart(0, Vh);
real[int] bcn = DirichletBC(0, Vh);
real[int] bcl =  bcn;


for (int t = 1; t <=100; t++){
    real[int] bnew = b0;
    bnew += M * Qn[];
    bnew = (bcn ? bcl : bnew);
    Qn[] = A^-1 * bnew;
    plot(Qn);
}

- According to the official documents, FreeFem++ uses an algorithm for the boundary penalty method in order to deal with Dirichlet boundary conditions. If you are interested in the algorithm, the below references can be helpful. 
  [1] https://www.sciencedirect.com/science/article/pii/0045782582900573   
  [2] https://www.jstor.org/stable/2005611?seq=1      

---

Nonlinear & Time-dependent diffential equations 

As an example, let's consider stoke's equations; 

 

 

Another example: Nonlinear Schrodinger equation

1. Reference: http://www.georges-sadaka.fr/FreeFem++/Non_linear_Schrodinger.edp  
- You can find the entire code from the above link. I will write a few comments about the code for my understanding. 

2.