Reveal the formal PDE of FiniteElement

Question

When FiniteElement method is used, the differential equations will first be transformed to certain standard form (named as formal PDE in recent FEM document), and it turns out to be critical to check what the standard form is when analyzing various issues related to FEM. Here are some examples:

Position of discontinuous coefficient influences the solution of PDE

How to input Robin boundary conditions for nonstandard Laplace equation?

Sign of conservative convection coefficient in a formal (Inactive) PDE

Stress analysis in axisymmetric bodies

The coefficient of formal PDE is available from PDECoefficientData, but its output is just hard to read. For example, with

{state} = NDSolve`ProcessEquations[
            With[{u = u[x, y]}, {-2 D[u, y, y] - 3 D[u, x, x] == 1, 
                                 DirichletCondition[u == 0, True]}], 
            u, {x, 0, 1}, {y, 0, 1}];

data = state["FiniteElementData"]["PDECoefficientData"];

data["All"]
(* {{{{1}}, {{{{0}, {0}}}}}, {{{{{3, 0}, {0, 2}}}}, {{{{0}, {0}}}}, {{{{0, 
  0}}}}, {{0}}}, {{{0}}}, {{{0}}}} *)

at hand, can you tell what's what? Can you label $d$, $c$, $\alpha$, etc. in the formal PDE

$$d\frac{\partial }{\partial t}u+\nabla \cdot (-c \nabla u-\alpha u+\gamma ) +\beta \cdot \nabla u+ a u -f=0$$

with corresponding values, without doubt?

Can we have a function that shows the formal PDE of FiniteElement in a easy-to-read way? A possible (but not necessary of course) input-output in my mind:

showFormalPDE@With[{u = u[x, y]}, -2 D[u, y, y] - 3 D[u, x, x] == 1]
(* -1 + Inactive[Div][(-{{3, 0}, {0, 2}}.Inactive[Grad][u[x, y], {x, y}]), {x, y}] == 0 *)

Answer 1

Coincidence had it that I needed a code to reconstruct the inactive PDE that has been parsed for a customer a few weeks back. I have then added this function to the kernel and it will be available in 12.2.

The details of the operators and their specification can be found in the documentation and @andre already added links to that documentation.

In versions after 12.2 you can use:

NDSolve`FEM`GetInactivePDE 

Else, here is the code to get the inactive PDE from the NDSolve state data:

Needs["NDSolve`FEM`"]
zeroCoefficientQ[c_] := Union[N[Flatten[c]]] === {0.}
ClearAll[GetInactivePDE]
GetInactivePDE[pdec_PDECoefficientData, vd_] := 
 Module[{lif, sif, dif, mif, hasTimeQ, tvar, vars, depVars, neqn, 
   nspace, dep, load, dload, diff, cconv, conv, react, 
   pde},
  {lif, sif, dif, mif} = pdec["All"];      

  tvar = NDSolve`SolutionDataComponent[vd, "Time"];
  If[tvar === None || tvar === {}, hasTimeQ = False;
   tvar = Sequence[];, hasTimeQ = True;];
  
  vars = NDSolve`SolutionDataComponent[vd, "Space"];
  depVars = NDSolve`SolutionDataComponent[vd, "DependentVariables"];
  neqn = Length[depVars];
  nspace = Length[vars];
  dep = (# @@ Join[{tvar}, vars]) & /@ depVars;
  
  {load, dload} = lif;
  {diff, cconv, conv, react} = sif;
  
  load = load[[All, 1]];
  dload = dload[[All, 1, All, 1]];
  conv = conv[[All, All, 1, All]];
  cconv = cconv[[All, All, All, 1]];      

  pde = If[hasTimeQ, 
    mif[[1]].D[dep, {tvar, 2}] + dif[[1]].D[dep, tvar], 
    ConstantArray[0, {Length[dep]}]];
  
  If[! zeroCoefficientQ[diff], 
   pde += (Plus @@@ 
       Table[Inactive[
          Div][-diff[[r, c]].Inactive[Grad][dep[[c]], vars], 
         vars], {r, neqn}, {c, neqn}]);];
  
  If[! zeroCoefficientQ[cconv], 
   pde += (Plus @@@ 
       Table[Inactive[Div][-cconv[[r, c]]*dep[[c]], vars], {r, 
         neqn}, {c, neqn}]);];
  
  If[! zeroCoefficientQ[dload], 
   pde += (Inactive[Div][#, vars] & /@ dload);];
  
  If[! zeroCoefficientQ[conv], 
   pde += (Plus @@@ 
       Table[conv[[r, c]].Inactive[Grad][dep[[c]], vars], {r, 
         neqn}, {c, neqn}]);];
  
  pde += react.dep;
  
  pde -= load;
  
  pde
  ]

GetInactivePDE[state_NDSolve`StateData] := Module[{femd = state["FiniteElementData"],
   vd = state["VariableData"], pdec},
  pdec = femd["PDECoefficientData"];
  GetInactivePDE[pdec, vd]]

Here is an example of it's usage:

op = -x D[u[x, y], {x, 2}] - D[u[x, y], {y, 2}] - 1;
{state} = 
  NDSolve`ProcessEquations[{op == 0, 
    DirichletCondition[u[x, y] == 0, True]}, 
   u, {x, y} ∈ Disk[]
   ];

pde = GetInactivePDE[state];
pde // InputForm

{-1 + {1, 0} . Inactive[Grad][u[x, y], {x, y}] + 
  Inactive[Div][-{{x, 0}, {0, 1}} . Inactive[Grad][u[x, y], {x, y}], {x, y}]}

Note, how the x in front of the D got pulled into the Div - Grad and how that is compensated by a convection component. See for example FEMDocumentation/tutorial/FiniteElementBestPractice#588198981 that explains this behavior.

Answer 2

I don't know if you are aware that this is documented in details.

The problem is that the informations are dispatched over the documentation of PDECoefficentData and InitializePDECoefficients.

your code :

{state} = 
  NDSolve`ProcessEquations[
   With[{u = u[x, y]}, {-2 D[u, y, y] - 3 D[u, x, x] == 1, 
     DirichletCondition[u == 0, True]}], u, {x, 0, 1}, {y, 0, 1}];

data = state["FiniteElementData"]["PDECoefficientData"];

data["All"]
(*{{{{1}},{{{{0},{0}}}}},{{{{{3,0},{0,2}}}},{{{{0},{0}}}},{{{{0,0}}}},\
{{0}}},{{{0}}},{{{0}}}}*)  

The PDECoefficentData documentation explains this :

data["ConvectionCoefficients"]
data["DampingCoefficients"]
data["MassCoefficients"]
data["LoadCoefficients"]
(* etc ... *) 

{{{{0, 0}}}}

{{0}}

{{0}}

{{1}}

InitializePDECoefficients documentation :

The DampingCoefficients and MassCoefficients are explained beyond.

Answer 3

The problem has been resolved by user21. I'd like to add a wrapper for NDSolve`FEM`GetInactivePDE to make it easier to use:

ClearAll[checkFormalPDE]

toeq[sys_] := 
  Replace[Flatten@{sys}, 
   expr : Except[_Equal | _DirichletCondition | _PeriodicBoundaryCondition] :> 
    expr == 0, {1}];

SetAttributes[checkFormalPDE, HoldAll]; 

With[{pe = NDSolve`ProcessEquations, gip = NDSolve`FEM`GetInactivePDE},
 checkFormalPDE[(head : NDSolve | NDSolveValue | NDEigenvalues | NDEigensystem)[sys_,
     depend_, domain__, n___Integer, rest : OptionsPattern[]]] :=
   Module[{state = Sow@First@pe[sys // toeq, depend, domain, rest]}, state // gip]];

Usage

Simply wrap it on NDSolve/NDSolveValue/NDEigenvalues/NDEigensystem:

checkFormalPDE@
 NDSolve[With[{u = u[x, y]}, {-2 D[u, y, y] - 3 D[u, x, x] == 1, 
    DirichletCondition[u == 0, True]}], u, {x, 0, 1}, {y, 0, 1}]
(* {-1 + Inactive[Div][(-{{3, 0}, {0, 2}} . Inactive[Grad][u[x, y], {x, y}]), {x, y}]} *)