FreeFEM Documentation on GitHub

stars - forks

Time dependent schema optimization for heat equations

First, it is possible to define variational forms, and use this forms to build matrix and vector to make very fast script (4 times faster here).

For example solve the ThermalConduction problem, we must solve the temperature equation in \(\Omega\) in a time interval (0,T).

\[\begin{split}\begin{array}{rcll} \partial_t u -\nabla\cdot(\kappa\nabla u) &=& 0 &\hbox{ in } \Omega\times(0,T)\\ u(x,y,0) &=& u_0 + x u_1\\ u &=& 30 &\hbox{ on } \Gamma_{24}\times(0,T)\\ \kappa\frac{\partial u}{\partial n} + \alpha(u-u_e) &=& 0 &\hbox{ on } \Gamma\times(0,T) \end{array}\end{split}\]

The variational formulation is in \(L^2(0,T;H^1(\Omega))\); we shall seek \(u^n\) satisfying:

\[\forall w \in V_{0};\ \int_\Omega \frac{u^n-u^{n-1}}{\delta t} w + \kappa\nabla u^n\nabla w) +\int_\Gamma\alpha(u^n-u_{ue})w=0\]

where \(V_0 = \{w\in H^1(\Omega)/ w_{|\Gamma_{24}}=0\}\).

So, to code the method with the matrices \(A=(A_{ij})\), \(M=(M_{ij})\), and the vectors \(u^n, b^n, b',b", b_{cl}\) (notation if \(w\) is a vector then \(w_i\) is a component of the vector).

\[\begin{split}u^n = A^{-1} b^n, \quad b' = b_0 + M u^{n-1}, \quad b"= \frac{1}{\varepsilon} \; b_{cl} , \quad b^n_i = \left\{ \begin{array}{cl} b''_i & \mbox{if }\ i \in \Gamma_{24} \\ b'_i & \mbox{else } \end{array}\right. \label{eq tgv}\end{split}\]

Where with \(\frac{1}{\varepsilon} = \mathtt{tgv} = 10^{30}\):

\[\begin{split}\begin{array}{rcl} A_{ij} &=& \left\{\begin{array}{cl} \frac{1}{\varepsilon} & \mbox{if } i \in \Gamma_{24}, \mbox{and}\quad j=i\\ \displaystyle{\int_{\Omega} w_j w_i / dt + k (\nabla w_j. \nabla w_i ) + \int_{\Gamma_{13}} \alpha w_j w_i} & \mbox{else} \end{array}\right.\\ M_{ij} &=& \left\{\begin{array}{cl} \frac{1}{\varepsilon} & \mbox{if } i \in \Gamma_{24}, \mbox{and}\quad j=i \\ \displaystyle n{\int_{\Omega} w_j w_i / dt} & \mbox{else} \end{array}\right. \\ b_{0,i} &=& n{\int_{\Gamma_{13}} \alpha u_{ue} w_i } \\ b_{cl} &=& u^{0} \quad \mbox{the initial data} \end{array}\end{split}\]

The Fast version script:

1...
2Vh u0=fu0, u=u0;

Create three variational formulation, and build the matrices \(A\),\(M\).

 1varf vthermic (u, v)
 2    = int2d(Th)(
 3          u*v/dt
 4        + k*(dx(u)*dx(v) + dy(u)*dy(v))
 5    )
 6    + int1d(Th, 1, 3)(
 7          alpha*u*v
 8    )
 9    + on(2,4,u=1)
10    ;
11
12varf vthermic0 (u, v)
13    = int1d(Th, 1, 3)(
14          alpha*ue*v
15    )
16    ;
17varf vMass (u,v)
18    = int2d(Th)(
19          u*v/dt
20    )
21    + on(2, 4, u=1)
22    ;
23
24real tgv = 1e30;
25matrix A = vthermic(Vh, Vh, tgv=tgv, solver=CG);
26matrix M = vMass(Vh, Vh);

Now, to build the right hand size; we need 4 vectors.

 1real[int] b0 = vthermic0(0,Vh); //constant part of RHS
 2real[int] bcn = vthermic(0,Vh); //tgv on Dirichlet part
 3real[int] bcl = tgv*u0[];   //the Dirichlet B.C. part
 4
 5// The fast loop
 6for(real t = 0; t < T; t += dt){
 7    real[int] b = b0;   //the RHS
 8    b += M*u[]; //add the the time dependent part
 9    b = bcn ? bcl : b; //do $\forall i$: b[i] = bcn[i] ? bcl[i] : b[i];
10    u[] = A^-1*b; //solve linear problem
11    plot(u);
12}
Table of content