numerical « Mesh Mess

Linear shape function, mass matrix, and element volume

Posted in FEM, Flows by Yi Zhang on 04/22/2010

In unsteady CFD problem, under certain time integration scheme, L2 inner product of shape functions gives what so called mass matrix, i.e. the matrix with its entries in the form $M_{ij}=\int_{\Omega}{N_i N_j}\,d{\Omega}$ , where $N_i$ is shape function. When shape function is linear, there is a better way to demonstrate $M_{ij}$ is proportional to local element volume (area) than trivially plugging in and calculating or using Gaussian integral. Assume in 2D triangular element ABC whose area is $A^{ABC}$ , and $N_i(x,y), i=1,2,3$ stands for the shape function whose unit value is reached at A, B, C correspondingly.

If we interpret local mass entry $M_{13}^{ABC}=(N_1, N_3)$ over element ABC shown above as the mass of a piece of wood splinter in the shape of tetrahedron A1-A-B-C, where $N_1(x,y)$ is the thickness function on the element, while $N_3(x,y)$ is the density function, then it’s not difficult to see that $M_{13}=(N_1,N_3)=(N_2,N_3)=M_{23}$ , because they are both the mass of tetrahedron with the same volume and density distribution from bottom (A1-A-B or B1-A-B) to top (C). The next step is to show that the mass of tetrahedron C1-A-B-C ( $M_{33}$ ) is the half of the total mass of the prism A1-B1-C1-A-B-C with the same density distribution. The argument is based on linearity, which allows us to assign a single parameter for above integrals. Set this single coordinate $s$ as the line parameter from C to line A-B, i.e. $s=0$ on AB, $s=1$ at C. Then, for mass of prism and tetrahedron, we both have

mass $=\int{L(s) H(s) D(s)}\,ds$

where $L\times H$ gives the area for cross section parallel to A1-A-B-B1 at $s$ . $H$ is the thickness of ABC, while $D$ is the density distribution. We have

$L\propto (1-s),\quad D\propto s$ ,
and $H=1$ for prism, $H\propto s$ for tetrahedron. Then we have

$M_p\propto \int{(1-s) s}\,ds$

$M_{33}\propto \int{(1-s) s^2}\,ds$

which proves

$2 M_{33}=M_p=\frac{1}{3} A_{ABC}$

Since $M_{33}+M_{13}+M_{23}$ is nothing but the sum of the mass of three tetrahedrons, which is the mass of the prism. This sum is then $\frac{1}{3}A_{ABC}$ . Using $M_{13}=M_{23}$ , finally we have

$M_{11}=M_{22}=M_{33}=\frac{1}{6}A_{ABC}$

$M_{12}=M_{13}=M_{23}=\frac{1}{12}A_{ABC}$

Note: for 3D tetrahedron element, linear shape function gives similar relation, in this case, above ratio over the volume of a tetrahedron would be 1/10 and 1/20:

$M_{11}=M_{22}=M_{33}=M_{44}=\frac{1}{10}V_{ABCD}$

$M_{12}=M_{13}=M_{14}=M_{23}=M_{24}=M_{34}=\frac{1}{20}V_{ABCD}$

Tagged with: FEM, mathematica, numerical

leave a comment

What’s the difference?

Posted in FEM by Yi Zhang on 01/04/2010

Today I was asked for a one-sentence comment of the difference between finite difference method and finite element method. What I came up was:

In finite difference method we approximate the differential expression, or the partial derivative, while in finite element method we approximate the space in which the solution lies.

Hopefully this shed some light on the issue.

Tagged with: FEM, Finite difference, numerical, PDE, Sobolev space

leave a comment

Saddle point problem in Navier-Stokes equation

Posted in Flows by Yi Zhang on 12/21/2009

It is well known that the the optimization problem

$\min J(x)= \frac{1}{2}(Ax, x)-(f, x)$

is equal to solving a linear operator equation

$Ax=f$

provided that $A$ is symmetric and positive definite. Here $x\in X$ which is a Hilbert space, and $f\in X'$ while $A$ is a bilinear operator.

This is a very general problem considered widely in many areas. However, in the real world, some further condition are to be applied on the space where $x$ lies on. If this is described as constraint equation

$Bx=g$

then we are dealing with constraint optimization problem ,to which Lagrange multiplier is usually the first one can come up. We then try to solve the stationary point of Lagrangian

$L(x,\lambda)=\frac{1}{2}(Ax, x)-(f, x)+\lambda (Bx-g)$

Now supposing we set $X=R^n$ , using stationary point condition on above equation, we arrive at

$\left[\begin{array}{cc} A & B^T\\ B & O \end{array}\right]\left[\begin{array}{c} x\\ y\end{array}\right]=\left[\begin{array}{c} f\\ g\end{array}\right]$

where $\lambda$ is replaced by $y$ indicating the parallel unknown as $x$ in the problem. The stationary point nature of this solution gives the name of saddle point problem, whose general case is

$L(x,v)\le L(x,y) \le L(u,y), \forall u,v \Leftrightarrow \min_u\max_v L(u,v)=L(x,y)=\max_v\min_u L(u,v)$

So, putting a constraint to a minimization problem, one try to solve a saddle problem instead of a simple linear system with SPD matrix.

For a linear system, one can always decompose it in to the form of

$\left[\begin{array}{cc} A & B_1^T\\ B_2 & -C \end{array}\right]\left[\begin{array}{c} x\\ y\end{array}\right]=\left[\begin{array}{c} f\\ g\end{array}\right]$

and this is the most general form of saddle problem.

For incompressible Navier-Stokes equation, we have (without mentioning the boundary conditions)

$\frac{\partial\mathbf{u}}{\partial t}+\mathbf{u}\cdot\nabla\mathbf{u}-\nu\Delta \mathbf{u}+\nabla p/\rho=\mathbf{f}$

$\nabla\cdot\mathbf{u}=0$

Following $\theta$ -scheme by Glowinski, one has the time discretized plan as:

$(\mathbf{u^{n+\theta}}-\mathbf{u}^n)/(\theta\Delta t)-\alpha\nu\Delta\mathbf{u}^{n+\theta}+\nabla p^{n+\theta}=\mathbf{f}^{n+\theta}+\beta\nu\Delta \mathbf{u}^n-(\mathbf{u}^n\cdot\Delta)\mathbf{u}^n$

$\nabla\cdot\mathbf{u}^{n+\theta}=0$

then

$(\mathbf{u^{n+1-\theta}}-\mathbf{u}^{n+\theta})/((1-2\theta)\Delta t)-\beta\nu\Delta\mathbf{u}^{n+1-\theta}+(\mathbf{u}^{n+1-\theta}\cdot\Delta)\mathbf{u}^{n+1-\theta}=\mathbf{f}^{n+\theta}+\alpha\nu\Delta \mathbf{u}^{n+\theta}-\nabla p^{n+\theta}$

then

$(\mathbf{u}^{n+\theta}-\mathbf{u}^{n+1-\theta})/(\theta\Delta t)-\alpha\nu\Delta\mathbf{u}^{n+1}+\nabla p^{n+1}=\mathbf{f}^{n+1}+\beta\nu\Delta \mathbf{u}^{n+1-\theta}-(\mathbf{u}^{n+1-\theta}\cdot\Delta)\mathbf{u}^{n+1-\theta}$

$\nabla\cdot\mathbf{u}^{n+1}=0$

The first and last step in above scheme are to solve Stokes problem like

$-\Delta\mathbf{u}+\nabla p=\mathbf{f}$

$\nabla\cdot \mathbf{u}=0$

Multiply $\mathbf{u}$ on both sides of the first equation, by using Green’s formula, we can see it is just the saddle point problem, with $A$ as $-\Delta$ and $B$ as $\nabla$ :

$\left[\begin{array}{cc} -\Delta & \nabla\\ \nabla & O \end{array}\right]\left[\begin{array}{c} \mathbf{u}\\ p\end{array}\right]=\left[\begin{array}{c} \mathbf{f}\\ 0\end{array}\right]$

Here pressure acts as Lagrange multiplier, and the Lagrangian is

$L(\mathbf{u}, p)=\frac{1}{2}|\Delta \mathbf{u}|_{\Omega}^2-(\mathbf{f}, \mathbf{u})+p(\nabla\cdot\mathbf{u}) \Longrightarrow$

$L(\mathbf{u})=\frac{1}{2}|\Delta \mathbf{u}|_{\Omega}^2-(\mathbf{f}, \mathbf{u})$

where $|\cdot|_{\Omega}$ is the seminorm on $\Omega$ and $\mathbf{u}$ is in $H_0^1(\Omega)$ .

Tagged with: CFD, numerical, PDE

leave a comment

Jacobi method demo

Posted in Numerical by Yi Zhang on 11/02/2009

Code below is to solve the static version of the equation in previous FDM code demo, with B.C. changed as homogeneous. The Jacobi method is used to replace the LAPACK solver which gives inexact method. The exact solution is
$u(x)=-x^2/2+x/2$
Using initial guess as
$u_0(x)=\sin(2\pi x)+0.5\sin(10\pi x)+0.3\sin(20\pi x)+0.1\sin(40\pi x)$
one can see that most of the computing cost is spent on smoothing lowest mode.


#include <stdio.h>
#include <math.h>
#include <iostream>
#include <fstream>
//#include <vector>
//#include <sunperf.h>
using namespace std;
//typedef std::vector<float> fvector_s;

//Maximum of an array.
float max(float a[], int n) {
 float v = a[0];
 for (int i=1; i<n; i++) {
 if (a[i] > v) {
 v = a[i];
 }
 }
 return v;
}

int main ()
{
 int n=500;
 float tmp, b[n], a_d[n], a_u[n-1], a_l[n-1],    \
 tol=1.e-7, x0[n], x1[n], x2[n], error=1.;
 float h=1./(n+1), pi=3.1416;
 int count=0;
 
 //initializing a
 for(int i=0; i<n; ++i){
 a_d[i]=2.0;
 }
 for(int i=0; i<n-1; i++){
 a_l[i]=-1.;
 a_u[i]=-1.;
 }
  
 //initializaing b[]
 for(int i=0; i<n; ++i){
 b[i]=h*h;
 }

 //initial guess
 for(int i=0; i<n; ++i){
 x0[i]=sin(2*pi*h*(i+1))+0.5*sin(10*pi*h*(i+1))+\
       0.3*sin(20*pi*h*(i+1))+0.1*sin(40*pi*h*(i+1));
 x1[i]=x0[i];
 }
 
 ofstream outfile;
 outfile.open("jac.out");
 do{
 count+=1;
 for(int i=0; i<n; ++i){
 if(i==0)
 x2[i]=(b[i]-a_u[i]*x1[i+1])/a_d[i];
 else if(i==n-1)
 x2[i]=(b[i]-a_l[i-1]*x1[i-1])/a_d[i];
 else
 x2[i]=(b[i]-a_u[i]*x1[i+1]-a_l[i-1]*x1[i-1])/a_d[i];
 }
 for(int i=0; i<n; i++){
 x1[i]=x2[i]-x1[i];
 }
 error=max(x1, n)/max(x0, n);
 for(int i=0; i<n; i++){
 x1[i]=fabs(x2[i]-(-(i*h)*(i*h)/2+i*h/2));
 }
 outfile << count << "\t" << error << endl;
 for(int i=0; i<n; ++i){
 x1[i]=x2[i];
 }
 }while(error>tol);
 cout << "Number of iterations:" << count << endl;
 for(int i=0; i<n; i++){
 outfile << count << "\t" << x2[i] << endl;
 }
 
 outfile.close();
 return 0;
}

Result of Jacobi method at different iteration

Tagged with: algorithm, numerical

leave a comment

A sample of toy FDM code

Posted in Numerical by Yi Zhang on 10/22/2009

Solving equation

$u_t-u''=f(x,u)$

$f(x,u)=-2+2*(\sin(u)-\sin(x^2))$

Time integration is full implicit, and during each time step Newton’s method is employed.

PROGRAM MAIN
!--------------------------------------------------
!Implicit time evolution of equation
!u_t-u''=f(x,u)
!f(x,u)=-2+r*(sin(u)-sin(x^2)),r=2
!--------------------------------------------------
! t=time points
! dt=time step
! h=neighboring nodes distance
! u=unknown
! x=node location
! n=No. of unknowns

!USE SUNPERF

IMPLICIT NONE
INTEGER(4) :: n
REAL(4) :: elapse,time(2),etime
REAL(8) :: t=0,dt=0.1
REAL(8) :: h
REAL(8),ALLOCATABLE :: u(:),x(:)
INTEGER(4) :: count,i

n=50
h=1/(DBLE(n)+1.)
ALLOCATE (x(n),u(n))
x=(/(i*h,i=1,n)/)
u=(/((i*h)**0.5,i=1,n)/)

OPEN(UNIT=100,FILE='data',STATUS='UNKNOWN',ACTION='&
     &WRITE')
DO i=1,n
   WRITE(100,*) t,x(i),u(i)
END DO

DO WHILE (t<1.0)    
CALL newton(n,h,dt,u,x)    
t=t+dt    
DO i=1,n       
WRITE(100,*) t,x(i),u(i)    
END DO END DO CLOSE(UNIT=100) 
elapse=etime(time) 
WRITE(*,*) 'elapsed:',elapse 
END PROGRAM MAIN 
!-------------------------------------------------- 
SUBROUTINE newton(n,h,dt,u,x)   
IMPLICIT NONE   
INTEGER(4),INTENT(IN) :: n   
REAL(8),INTENT(IN) :: x(n),h,dt   
REAL(8),INTENT(INOUT) :: u(n)   
REAL(8) :: r=2,c(n),error,b(n,1),b0(n,1),a(2,n)   
REAL(8) :: alpha=1.0,beta=-1.0   
INTEGER(4) :: info,i,iter   
c=u*h*h/dt   
iter=0   
error=1   
DO WHILE (error>1.e-13)
     iter=iter+1
     a(1,:)=-1
     a(2,:)=2+h*h/dt
     b(:,1)=(-2+r*(SIN(u)-SIN(x*x)))*h*h+c
     b(n,1)=b(n,1)+1
     CALL DSBMV('u',n,1,alpha,a,2,u,1,beta,b,1)
     b=-b
     a(2,:)=a(2,:)-r*COS(u)*h*h
     CALL DGTSV(n,1,a(1,2:),a(2,:),a(1,2:),b,n,info)
     error=MAXVAL(ABS(b))
     u=u+b(:,1)
  END DO
  WRITE(*,*) 'INFO',info,'Number of iterations:',iter
END SUBROUTINE newton
!--------------------------------------------------

Output would be like:

% f95 -dalign implicit_newton.f95 -xlic_lib=sunperf
% a.out
 INFO 0 Number of iterations: 4
 INFO 0 Number of iterations: 4
 INFO 0 Number of iterations: 4
 INFO 0 Number of iterations: 4
 INFO 0 Number of iterations: 4
 INFO 0 Number of iterations: 3
 INFO 0 Number of iterations: 3
 INFO 0 Number of iterations: 3
 INFO 0 Number of iterations: 3
 INFO 0 Number of iterations: 3
 elapsed: 0.061407

Of course the exact solution is $u(x)=x^2$ . Here is time evolution plot of the numerical solution.

Tagged with: Finite difference, numerical, time integration

1 comment

Older Posts »

M	T	W	T	F	S	S
« Mar
	1	2	3	4	5	6
7	8	9	10	11	12	13
14	15	16	17	18	19	20
21	22	23	24	25	26	27
28	29	30	31

Mesh Mess

Linear shape function, mass matrix, and element volume

What’s the difference?

Saddle point problem in Navier-Stokes equation

Jacobi method demo

A sample of toy FDM code

Recent Posts

Tags

Meta

Blogroll

Follow “Mesh Mess”