Writing PDE Solvers

The Finite Element Method Interface

To Solve PDEs

Needs["NDSolve`FEM`"]
$HistoryLength=0;

Problem setup

{x,y,..}∈Ω

PDE

Γ

Solve and Visualize

Solving with NDSolve

Problem Setup

{x,y,..}∈Ω

Ω=Rectangle[{0,0},{5,1}];

PDE

pde=-Laplacian[u[x,y],{x,y}]-1/2

Γ

Γ={DirichletCondition[u[x,y]==0,x==0],
DirichletCondition[u[x,y]==1,x==5]};

Solving with NDSolve

Solve

uif=NDSolveValue[{pde==0,Γ},u,{x,y}∈Ω]

Visualize

Plot3D[uif[x,y],{x,y}∈Ω]

You can stop here - You can do a lot with that

Why an Interface?

Understandable what NDSolve does

Extensibility - for you and for me

Interceptable

NDSolve`

ProcessEquations

{x,y,..}∈Ω

PDE

Γ

Iterate

Solves

ProcessSolution

InterpolatingFunction

NDSolve` steps

{ndstate}=NDSolve`ProcessEquations[{pde==0,Γ},u,{x,y}∈Ω]

NDSolve`Iterate[ndstate]

NDSolve`ProcessSolutions[ndstate]

NDSolve`

ProcessEquations

{x,y,..}∈Ω

PDE

Γ

Iterate

Solves

ProcessSolution

InterpolatingFunction

Inside ProcessEquations

{ndstate}=NDSolve`ProcessEquations[{pde==0,Γ},u,{x,y}∈Ω]

feData=ndstate["FiniteElementData"]

Inside ProcessEquations

Problem setup

{x,y,..}∈Ω

Variable Data

Solution Data

PDE

InitializePDECoefficients

Γ

InitializeBoundaryConditions

{x,y,..}∈Ω

InitializePDEMethodData

Inside ProcessEquations

{x,y,..}∈Ω

vd=NDSolve`VariableData[{"DependentVariables","Space"}->{{u},{x,y}}];

nr=ToNumericalRegion[Ω];
sd=NDSolve`SolutionData["Space"->nr];

PDE

InitializePDECoefficients[vd,sd,"DiffusionCoefficients"->{{-IdentityMatrix[2]}},"LoadCoefficients"->{{-1/2}}]

Γ

InitializeBoundaryConditions[vd,sd,{Γ}]

{x,y,..}∈Ω

InitializePDEMethodData[vd,sd]

Inside ProcessEquations

{ndstate}=NDSolve`ProcessEquations[{pde==0,Γ},u,{x,y}∈Ω]

feData=ndstate["FiniteElementData"]

feData["Properties"]

{pdeData,bcData,methodData}={feData["PDECoefficientData"],feData["BoundaryConditionData"],feData["FEMMethodData"]};

Inside ProcessEquations

NDSolve`ProcessEquations is an equation pre-processor

NDSolve`

ProcessEquations

{x,y,..}∈Ω

PDE

Γ

Iterate

Solves

ProcessSolution

InterpolatingFunction

Inside Iterate

dPDE=DiscretizePDE[pdeData,methodData,sd]

dBC=DiscretizeBoundaryConditions[bcData,methodData,sd]

Inside Iterate

{l,s,d,m}=dPDE["SystemMatrices"]

DeployBoundaryConditions[{l,s,d,m},dBC]

Inside Iterate

result=LinearSolve[s,l];

NDSolve`

ProcessEquations

{x,y,..}∈Ω

PDE

Γ

Iterate

Solves

ProcessSolution

InterpolatingFunction

Visualize

uif=ElementMeshInterpolation[sd,result]

Plot3D[uif[x,y],{x,y}∈Ω]

Summary

PDEtoMatrix[{pde_,Γ___},u_,r__]:=Module[{ndstate,feData,sd,bcData,methodData,pdeData},
{ndstate}=
NDSolve`ProcessEquations[Flatten[{pde,Γ}],u,Sequence@@{r}];
sd=ndstate["SolutionData"][[1]];
feData=ndstate["FiniteElementData"];
pdeData=feData["PDECoefficientData"];
bcData=feData["BoundaryConditionData"];
methodData=feData["FEMMethodData"];
{DiscretizePDE[pdeData,methodData,sd],DiscretizeBoundaryConditions[bcData,methodData,sd],sd,methodData}
]

PDEtoMatrix[{pde==0,Γ},u,{x,y}∈Ω];
dPDE

Let’s do something with it....

NDEigensystem

{dPDE,dBC,sd,md}=PDEtoMatrix[{D[u[t,x,y],t]==Laplacian[u[t,x,y],{x,y}],u[0,x,y]==0},u,{t,0,1},{x,y}∈Ω];

s=dPDE["StiffnessMatrix"];
d=dPDE["DampingMatrix"];

First[Reverse/@Eigensystem[{s,d}, -4]]

First[NDEigensystem[{-(1/2)-Derivative[0,2][u][x,y]-Derivative[2,0][u][x,y]},u,{x,y}∈Rectangle[{0,0},{5,1}],4,Method->"Direct"]]

NDSingularValueDecomposition
- Why not?

SingularValueDecomposition[s,4];

All of linear algebra....

NDLeastSquares

Matrix based minimization

NDxyDecomposition

....

NIntegrate

You want to guess...

NIntegrate[1,{x,y}∈Disk[],Method->"FiniteElement"]

New Boundary Condition

uif=NDSolveValue[(Subscript[∇, {x,y}]^2)u[x,y]==NeumannValue[Sin[x y],True],u,{x,y}∈Disk[]]

Plot3D[uif[x, y],{x, y}∈Disk[]]

NIntegrate[uif[x,y],{x,y}∈uif["ElementMesh"]]

Integral Boundary Conditions

{dPDE,dBC,sd,md}=PDEtoMatrix[{(Subscript[∇, {x,y}]^2)u[x,y]==NeumannValue[Sin[x y],True]},u,{x,y}∈Disk[]];

l=dPDE["LoadVector"];
s=dPDE["StiffnessMatrix"];
DeployBoundaryConditions[{l,s},dBC]

Integral Boundary Conditions

mesh=md["ElementMesh"];
dof=md["DegreesOfFreedom"];

dc=mkIntegralConstraint[43,mesh,dof]

DeployBoundaryConditions[{l,s},dc,"ConstraintMethod"->"Append" ];
sol=LinearSolve[s,l];

if2=ElementMeshInterpolation[{mesh},Take[sol,dof]];
Plot3D[if2[x,y],{x,y}∈Disk[]]

NIntegrate[if2[x,y],{x,y}∈mesh]

Transient PDEs

bmesh=;
Show[
bmesh,
ToBoundaryMesh[bmesh]["Wireframe"["MeshElement"->"BoundaryElements","MeshElementStyle"->{Directive[EdgeForm[],FaceForm[]],Directive[EdgeForm[Black],FaceForm[Red]],
Directive[EdgeForm[Black],FaceForm[Blue]]}]]]

mesh=ToElementMesh[bmesh,"MeshOrder"->1];

Length[mesh["Coordinates"]]

{dpde,dbcs,sd,md}=PDEtoMatrix[{D[u[t,x,y,z],t]==Laplacian[u[t,x,y,z],{x,y,z}],u[0,x,y,z]==0,DirichletCondition[u[t,x,y,z]==1,ElementMarker==1],DirichletCondition[u[t,x,y,z]==0,ElementMarker==2]},u[t,x,y,z],{t,0,1},{x,y,z}∈mesh];

Transient PDEs II

Extract all system matrices and deploy the boundary conditions.

{l,s,d,m}=dpde["SystemMatrices"];

DeployBoundaryConditions[{l,s,d},dbcs]

Set up initial conditions based on the boundary conditions.

init=Table[{0.},{md["DegreesOfFreedom"]}];
init[[dbcs["DirichletRows"]]]=dbcs["DirichletValues"];

Dynamic["time: "<>ToString[CForm[currentTime]]]
AbsoluteTiming[
tif=NDSolveValue[{
d.u'[t]+s.u[t]==l
,u[0]==init},u,{t,0,1}
,Method->{"EquationSimplification"->"Residual"}
,Jacobian->{Automatic,Sparse->d["PatternArray"]}
,EvaluationMonitor:>(currentTime=t;)
]]

Transient PDEs III

Visualize

allGraphics=
ElementMeshSurfacePlot3D[tif[#],mesh,Boxed->False,Background->Gray,ColorFunctionScaling->False]&/@{0.01,0.05,0.1,0.2,0.4,0.6,0.8,1.};

GraphicsGrid[Partition[allGraphics,4],ImageSize->Large]

Model order reduction

mor[{load_,stiffness_,damping_},size_Integer]:=Module[{lsf,tmp,vt,v},
lsf=LinearSolve[-stiffness];
vt=Orthogonalize[Transpose[-lsf[damping.load]]];
Do[
tmp=lsf[damping.vt[[-1]]];
vt=Orthogonalize[Join[vt,{tmp}]];
,{i,2,size}];
v=Transpose[vt];
{{vt.load,vt.stiffness.v,vt.damping.v},v}
]

Find an approximate system of equations that has only (Subscript[n, r])=50 degrees of freedom.

AbsoluteTiming[{sys,v}=mor[{l,s,d},50];]

Extract the reduced system matrices and inspect the dimensions.

{loadR,stiffnessR,dampingR}=sys;
Dimensions/@{loadR, stiffnessR, dampingR}

MOR II

Efficiently time integrate the reduced system.

initR=Transpose[v].init;
dinitR=LinearSolve[dampingR,-(stiffnessR.initR-loadR)];
Dynamic["time: "<>ToString[CForm[currentTime]]]
AbsoluteTiming[
{tifR} = u /. NDSolve[{
dampingR.u'[t] +stiffnessR.u[t] == loadR
,u[0] ==initR
,u'[0] ==dinitR
}, {u}, {t,0,1}
,Method->{"EquationSimplification"->"Residual"}
,EvaluationMonitor:>(currentTime=t;)
 ];
]

Compare the (Subscript[l, ∞]) norm of the full system and the reduced system.

nrms={};
errorGraphics= (
vals = tif[#];
valsR = v.tifR[#];
nrms = {nrms,Norm[vals- valsR,Infinity]};
)& /@ Range[0,1,0.1];
ListLinePlot[Flatten[nrms],PlotRange->All]

Documentation

All low level FEM functions are documented!

Initialization
InitializePDECoefficients — initialize partial differential equation coefficients
InitializeBoundaryConditions — initialize boundary conditions
InitializePDEMethodData — initialize partial differential equation method data

PDECoefficientData ▪ BoundaryConditionData ▪ FEMMethodData

Discretization
DiscretizePDE — discretize initialized partial differential equations
DiscretizeBoundaryConditions — discretize initialized boundary conditions

DiscretizedPDEData ▪ DiscretizedBoundaryConditionData

Solution
DeployBoundaryConditions — deploy discretized boundary conditions into discretized partial differential equations
LinearSolve — solve linear systems of equations

Tutorial

Finite Element Programming

Guide page

Finite Element Method Guide

Thanks - NDSolve rocks!

ruebenko@wolfram.com

Non-linear PDEs

Linearize

Interpolate dependent and re-instert in PDE

Navier-Stokes

Linear part

μ=10^-3;ρ=1;
l=2.2;h=0.41;

Ω=RegionDifference[Rectangle[{0,0},{l,h}],ImplicitRegion[(x-1/5)^2+(y-1/5)^2<(1/20)^2,{x,y}]];

RegionPlot[Ω,AspectRatio->Automatic]

Γ={
DirichletCondition[p[x,y]==0.,x==l],
DirichletCondition[{u[x,y]==4*0.3*y*(h-y)/h^2,v[x,y]==0},x==0],
DirichletCondition[{u[x,y]==0.,v[x,y]==0.},y==0||y==h||(x-1/5)^2+(y-1/5)^2≤(1/20)^2]};

stokes={
D[u[x,y],x]+D[v[x,y],y],
Div[{{-μ,0},{0,-μ}}.Grad[u[x,y],{x,y}],{x,y}]+D[p[x,y],x],
Div[{{-μ,0},{0,-μ}}.Grad[v[x,y],{x,y}],{x,y}]+D[p[x,y],y]
};

{dPDE,dBC,sd,md}=PDEtoMatrix[{stokes=={0,0,0},Γ},{p,u,v},{x,y}∈Ω,Method->{"FiniteElement","InterpolationOrder"->{p->1,u->2,v->2},"MeshOptions"->{"ImproveBoundaryPosition"->False}}];

linearLoad=dPDE["LoadVector"];
linearStiffness=dPDE["StiffnessMatrix"];

vd=md["VariableData"];
offsets=md["IncidentOffsets"];

Navier-Stokes

(*uOld=stationarySolution;*)
uOld=ConstantArray[{0.},md["DegreesOfFreedom"]];
mesh2=md["ElementMesh"];
mesh1=MeshOrderAlteration[mesh2,1];

ClearAll[rhs]
rhs[t_?NumericQ,ut_]:=Module[{uOld},
uOld=ut;
Do[
ClearAll[u0,v0,p0];

p0=ElementMeshInterpolation[{mesh1},uOld[[offsets[[1]]+1;;offsets[[2]]]]];
u0=ElementMeshInterpolation[{mesh2},uOld[[offsets[[2]]+1;;offsets[[3]]]]];
v0=ElementMeshInterpolation[{mesh2},uOld[[offsets[[3]]+1;;offsets[[4]]]]];

nlPdeCoeff=InitializePDECoefficients[vd,sd,
"LoadCoefficients"->{(* F *)
{-(D[u0[x,y],x]+D[v0[x,y],y])},
{-ρ(u0[x,y]*D[u0[x,y],x]+v0[x,y]*D[u0[x,y],y])-D[p0[x,y],x]},
{-ρ(u0[x,y]*D[v0[x,y],x]+v0[x,y]*D[v0[x,y],y])-D[p0[x,y],y]}
},
"LoadDerivativeCoefficients"->-{(* gamma *)
{{0,0}},
{{μ D[u0[x,y],x],μ D[u0[x,y],y]}},
{{μ D[v0[x,y],x],μ D[v0[x,y],y]}}
},
"ConvectionCoefficients"->{(*beta*)
{{{0,0}}, {{0,0}},{{0,0}}},
{{{0,0}}, {{ρ u0[x,y],ρ v0[x,y]}},{{0,0}}},
{{{0,0}}, {{0,0}},{{ρ u0[x,y],ρ v0[x,y]}}}
},
"ReactionCoefficients"->{(* a *)
{0,0,0},
{0,ρ D[u0[x,y],x],ρ D[u0[x,y],y]},
{0,ρ D[v0[x,y],x],ρ D[v0[x,y],y]}
}
];

nlsys=DiscretizePDE[nlPdeCoeff,md,sd];
nlLoad=nlsys["LoadVector"];
nlStiffness=nlsys["StiffnessMatrix"];

ns=nlStiffness+linearStiffness;
nl=nlLoad+linearLoad;

DeployBoundaryConditions[{nl,ns},dBC];
diriPos=dBC["DirichletRows"];
nl[[ diriPos ]]=nl[[ diriPos ]]-uOld[[diriPos]];

dU=LinearSolve[ns,nl];
Print[ i, " Residual: ", Norm[nl, Infinity], "  Correction: ", Norm[ dU, Infinity ]];
uOld=uOld+dU;

(*If[Norm[ dU, Infinity ]<10^-6,Break[]];*)

,{i,8}
];
uOld
]

uNew=rhs[0,uOld];

p0=ElementMeshInterpolation[{mesh1},uNew[[offsets[[1]]+1;;offsets[[2]]]]];
u0=ElementMeshInterpolation[{mesh2},uNew[[offsets[[2]]+1;;offsets[[3]]]]];
v0=ElementMeshInterpolation[{mesh2},uNew[[offsets[[3]]+1;;offsets[[4]]]]];

ContourPlot[u0[x,y],{x,y}∈mesh2,AspectRatio->Automatic, PlotRange-> All, ColorFunction-> ColorData["TemperatureMap"],Contours->10,ImageSize->Large]

ContourPlot[v0[x,y],{x,y}∈mesh2,AspectRatio->Automatic,PlotRange-> All, ColorFunction-> ColorData["TemperatureMap"], Contours->10,ImageSize->Large]

ContourPlot[p0[x,y],{x,y}∈mesh1,AspectRatio->Automatic,PlotRange-> All, ColorFunction-> ColorData["TemperatureMap"], Contours->10,ImageSize->Large]

Init

Needs["NDSolve`FEM`"]
$HistoryLength=0;

mkIntegralConstraint[level_?NumberQ, mesh_,dof_Integer]:=
Module[{s1,s2,d1,diriRows,diriValMat},
s1=SparseArray[{},{dof,1}];
s2=SparseArray[{},{dof,dof}];
d1=SparseArray[{ConstantArray[1., {dof}]}];
diriRows={1};
(* level the pure neumann to a specific value *)
diriValMat= SparseArray[{{N[level*dof/Total[mesh["MeshElementMeasure"],2]]}}];
DiscretizedBoundaryConditionData[{s1,s2,d1,diriRows,diriValMat,{dof,0,Length[diriRows]}},1]
]