14
$\begingroup$

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 *)
$\endgroup$

2 Answers 2

13
$\begingroup$

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.

$\endgroup$
15
  • 1
    $\begingroup$ Side note: NDSolve`FEM`GetInactivePDE is built-in in v12.2. (There's no separate document page for it though. ) $\endgroup$
    – xzczd
    Apr 29, 2021 at 6:21
  • 1
    $\begingroup$ I take the liberty to add a line making GetInactivePDE in this post handle NDSolve`StateData, feel free to roll back if you don't like it :) . $\endgroup$
    – xzczd
    Apr 30, 2021 at 13:25
  • $\begingroup$ @user21 Looking at the inactive form of the pde: The Div-part defines the NeumannValue (-{{x, 0}, {0, 1}} . Grad [U[x, y], {x, y}]) . Is this conclusion correct? Thanks! $\endgroup$ May 16 at 14:08
  • $\begingroup$ @UlrichNeumann, yes, whatever is in the Div part defines NeumannValue. $\endgroup$
    – user21
    May 17 at 11:22
  • $\begingroup$ @user21 Thank you, in this question Understanding NeumannValue I 'm trying, in vain, to verify this result. Please read my question. $\endgroup$ May 18 at 9:47
8
$\begingroup$

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 :

enter image description here

The DampingCoefficients and MassCoefficients are explained beyond.

$\endgroup$
2
  • $\begingroup$ Oh, I didn't notice the Details section of InitializePDECoefficients… but I think it'll be good to have a function that can generate the formal PDE automatically, checking the table in the document is still a bit painful. $\endgroup$
    – xzczd
    Jul 11, 2020 at 8:09
  • 1
    $\begingroup$ It took me a while to notice this too (I don't rememder exactly, but it was after having read the tutorials). $\endgroup$
    – andre314
    Jul 11, 2020 at 8:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.