PDEs : automatic method choice : TensorProductGrid or FiniteElement?

Question

A problem that arises when we solve PDEs with NDSolve[] and Method->Automatic is how to know which method has been chosen : TensorProductGrid or FiniteElement ?

The question is important because :

• NDSolve often use TensorProductGrid (wrongly I think) instead of FiniteElement. I guess the reason is historical : TensorProductGrid is maybe 15 years older than FiniteElement. For compatibility reasons all forms of syntax that was leading to TensorProductGrid should continue to lead to this method.

• The wrong choice may lead to problems that are very, very advanced (here for example)

• If we don't know which method is used, we don't know which documentation to read (the one for TensorProductGrid or the one for FiniteElement) and these documentations are very consequent (and it's crucial to understand error messages).
• Even advanced user can't deduce the choice of the method from the syntax. For example, the use of the expression NeumannValue impose FiniteElement, but Laplacian not, though historically they appeared at the same time.

Here are some ways to find the answer to this question :

• some error messages may give the information
• after having done sol=NDSolveValue[...], then if sol["ElementMesh"] gives a object ElementMesh[ blabla ], FiniteElement has been used
• one can do also NDSolve`ProcessEquations[...] instead of NDSolve[...] and examine the resulting object.

Is there something more friendly, above all for beginners who will use preferently the automatic choice ?

Answer 1

Update:

There is a new section in the FEM documentation What Triggers the Use of the Finite Element Method that has a detailed list of what triggers the usage of the finite element method.

Previous:

Ideally, the method you choose does not matter because either finds a solution. However, currently (V12.0) the TensorProductGrid (TPG) can not deal with arbitrarily shaped regions.

I understand your question that you are looking for a simple rule of thumb. The simplest rules that trigger FEM I can think of are the following:

• If you specify a region ({x,...} ∈ reg)
• If you use DirichletCondition, NeumannValue or PeriodicBoundaryCondition
• If the PDE is elliptic and embedding dimension is larger than 1
• If the PDE contains Inactive components
• If you write a wave equation but do not specify (enough) initial conditions. The equation will be treated as stationary PDE and solved with the FEM

Specifying NeumannValue will lead to FEM since the correspondence to Derivative is not one to one (See the section The Relation between NeumannValue and Boundary Derivatives).

Most reliable is to check if the resulting interpolation function contains an ElementMesh.

if = NDSolveValue[{-Laplacian[u[x, y], {x, y}] == 1, u[0, y] == 0, 
    u[1, y] == 0}, u, {x, 0, 1}, {y, 0, 1}];
if["ElementMesh"]

(* NDSolve`FEM`ElementMesh[{{0., 1.}, {0., 1.}}, 
{NDSolve`FEM`QuadElement["<" 400 ">"]}]*)

In case the TPG was used None will be returned.

Another hint is that all FEM error message tags start with fem. So if you see one of those you know NDSolve is trying to using FEM.

I'll continue to work to make the FEM more versatile in future releases.