Domain Decomposition for Poisson equation in 1D

It’s not really Poisson equation, just a second order ODE in that sense. But it serves a good example to demonstrate the idea of domain decomposition.

Using Dirichlet-Neumann method, one of the two subdomains gives out Neumann boundary condition and receives Dirichlet boundary data at the interface, and the other one does the opposite. The Mathematica modul would be somthing like this:

DNDecomp:   Main module;
DBCLaplace: linear FEM module for Dirichlet data ua at xa, and ub at xb;
MBCLaplace: linear FEM module for Mixed data ua at xa, and flux db at xb;
(xi,ui):    interface position and data;
theta:      relaxation coefficient to update ui in each iteration.

DNDecomp[xa_, ua_, xb_, ub_, xi_, ui0_, n1_, n2_, iter_, \[Theta]_] :=
Module[{ui, di, x1, u1, x2, u2, h1, h2, sol},
  f[t_] := 2;
  ui = ui0;
  h1 = (xi - xa)/n1; h2 = (xb - xi)/n2;
  sol = {};
  Do[{
    {x2, u2} = DBCLaplace[xi, ui, xb, ub, n2];
    di = (u2[[2]] - u2[[1]])/h2;
    {x1, u1} = MBCLaplace[xa, ua, xi, di, n1];
    ui = \[Theta] u1[[-1]] + (1 - \[Theta]) u2[[1]];
    AppendTo[sol, {{x1, u1}\[Transpose], {x2, u2}\[Transpose]}]},
   {iter}
   ];
  sol
  ]

In above example force function gives a parabolic exact solution. The effect of the relaxation coefficient is obvious in plots below.

Above: theta=0.7, below: theta=0.5

Below is another example with exact solution .

Iteration = 6, theta = 0.5

Three faces of FEM

Finite Element Method (FEM) acquires several interpretations through its developments over decades. Assuming we are solving a PDE of form

where is partial-differential operator, is unknown function. FEM could be interpreted in the following ways (let me know if you have extra answers):

• Rayleigh-Ritz
• Galerkin method
• Method of Weighted Residuals (MWR)

The first one, Rayleigh-Ritz method, is actually how structure engineers were inspired to devise FEM. The connection of this method and FEM is the two way street of variational formulation of PDE, which gives the Euler equation of certain functional over . The derivation is done through searching for the minimum of the functional (usually some form of energy ). Rayleigh-Ritz method addresses the PDE problem by assuming the form of as

then bringing this formulation into PDE’s functional formulation, in which ‘s are acquired by solving a minimization problem. Original Rayleigh-Ritz method adopts functions with global support as . A textbook example is the boundary value problem of a beam’s modes, in which we can assume as sinusoidal functions. are called trial functions.

In order to have a PDE’s variational formulation in a general way, and in the same time remove some regularity requirement of the solution, Galerkin method is used. The so called test function are the ones with expected regularity (usually infinitive). By the procedure of Galerkin method, test functions are applied to and the integral formulation is acquired. If is linear, this integral formulation is the same as functional formulation in minimization problem by Rayleigh-Ritz method. Here we have introduced two sets of functions, trial functions which define the finite dimensional approximation of , and test functions which defines how close (accurate) the PDE is solved. The name of “test” can be understood in the following way. How close a variable is close to zero could be measured by test it with other variables:

where is known variable used for testing, and is proper inner product. Conceptually, the greater the variable domain consisting of , the closer to zero is. And if can make every satisfy above equation, it should be zero itself by definition.

So we can see, the space of test function defines the extension, and in turn, accuracy, the PDE stands by having

A subcategory of Galerkin method is Bubnov-Galerkin method, where test functions and trial functions are the same. Otherwise, it’s called Petrov-Galerkin method. Historically, trial function is also commonly referred as shape function.

Another abstraction of FEM is by MWR. Since all the numerical methods for solving PDE relies on moving from infinite dimensional solution spaces to finite ones, errors are introduced in all the methods. So instead of having equal to exact zero, ‘s numerical approximation introduces residual :

If numerical solution is in the form of trial function’s expansion, the coefficients (coordinates) could be acquired by restricting . In the light of above philosophy, could be restricted as to be zero in some sense, specifically,

for some and . This can be looked as putting weight on some locations defined by : the greater is on certain location, the more restricted is there. By this ‘s are also called weighted functions.

Though mathematically Galerkin method and MWR are in the same look, they are actually on different emphasis. In the former, we are specifically looking for integral form, with shifting regularity to test functions in mind, while in MWR we are focusing on minimizing residuals, and following integral form is just the math “trick” played on certain weight functions. In fact, by using different weight function in MWR, we can recover other numerical discretizations. Besides Bubnov/Petrov Galerkin methods, gives collocation method, within certain cell and otherwise gives FVM, and gives least squares method.

Log law at a smooth plate

Log-law is THE most important and elegant analytical result in turbulence, considering many others are dauntingly complicated,  though justified. What log law says is  the fact that for a turbulent flow over flat plate, there is this section of flow away from the wall, within the turbulent boundary layer,  in which the relation between the distance and mean flow velocity is simply of log function. This conclusion could be drawn by an asymptotic analysis matching outer flow and inner flow.   Any CFD trying tempted to address turbulence models should recover the “law of the wall“, among which log law is a part of, before moving forward to more advanced problems.

The delicacy of the modeling in this case is the boundary condition (shocking!) at the wall, because log law is only valid when y+ is approximately above 30 ( the other cap depends on Reynolds number ), and the boundary conditions there for unknown variables in turbulence models are generally not known a priori, where is all kinds of wall function come into play.

Channel Flow

Plot of CFD result against log law

Linear shape function, mass matrix, and element volume

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 , where is shape function. When shape function is linear, there is a better way to demonstrate 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 ,  and stands for the shape function whose unit value is reached at A, B, C correspondingly.

If we interpret local mass entry over element ABC shown above as the mass of a piece of wood splinter in the shape of tetrahedron A1-A-B-C, where is the thickness function on the element, while is the density function, then it’s not difficult to see that , 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 () 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 as the line parameter from C to line A-B, i.e. on AB, at C. Then, for mass of prism and tetrahedron, we both have

mass

where gives the area for cross section parallel to A1-A-B-B1 at . is the thickness of ABC, while is the density distribution. We have

,
and for prism, for tetrahedron. Then we have

which proves

Since is nothing but the sum of the mass of three tetrahedrons, which is the mass of the prism. This sum is then . Using , finally we have

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:

What’s the difference?

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.