4
$\begingroup$

I want to solve a PDE, but I'm struggling a lot getting it into Mathematica.

This is what I currently have (thanks to @user21). However, there are still a couple of issues which I explain below.

Needs["NDSolve`FEM`"]
{dpde, dbc, vd, sd, md} = 
 ProcessPDEEquations[{D[C1[x, y]*u[x, y], x] + D[C2[x, y]*u[x, y], y] +
 D[C3[x, y]*u[x, y], x, x] + D[C4[x, y]*u[x, y], x, y] + 
D[C5[x, y]*u[x, y], y, y] == 
    NeumannValue[0, x <= 0] + NeumannValue[0, y <= 0] + 
     NeumannValue[0, x + y >= 1], 
   DirichletCondition[u[x, y] == 0, x == 1 && y == 0]}, 
  u, {x, y} \[Element] region, 
  Method -> {"PDEDiscretization" -> {"FiniteElement", 
      "MeshOptions" -> {"MaxCellMeasure" -> 0.000125}}}]
l = dpde["LoadVector"];
s = dpde["StiffnessMatrix"];
DeployBoundaryConditions[{l, s}, dbc];
mesh = md["ElementMesh"];
dof = md["DegreesOfFreedom"];

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]]

dc = mkIntegralConstraint[1, mesh, dof]
DeployBoundaryConditions[{l, s}, dc, "ConstraintMethod" -> "Append"];

sol = LinearSolve[s, l];
if2 = ElementMeshInterpolation[{mesh}, Take[sol, dof]];
Plot3D[if2[x, y], {x, y} \[Element] region, PlotRange -> All]

However, I believe for these parameter values, the peak of the pdf should not be occurring at the edge of the domain like this.

$\endgroup$
3
  • 2
    $\begingroup$ You've said the derivative is zero on the boundary and there is no flux, but for this particular equation isn't $u(x,y)=0$ a valid solution as your source condition u[x,1] == 1 is a single point on the top left corner of the triangle? There are no source terms in your equation, and no non-zero boundary conditions. $\endgroup$
    – flinty
    Jun 24 '20 at 13:43
  • $\begingroup$ @flinty Ah yeah, thanks for reminding me. I left out another condition from the problem. I also need the condition that u integrates to 1 over the domain! $\endgroup$
    – user112495
    Jun 24 '20 at 14:00
  • $\begingroup$ @flinty Also, I made a mistake above, the absorbing state is when $x=1$, not $y=1$. (So $u[1,y]$, which on my domain is only the point $u[1,0]$). $\endgroup$
    – user112495
    Jun 24 '20 at 14:40
4
$\begingroup$

Here I am sketching out a solution path. My feeling is that the equations are not quite right, or something with the setup. More on this later.

This setup solves the equation given. Make sure it is correct as the examples shown in the question differ slightly.

uval = NDSolveValue[{D[C1[x, y]*u[x, y], x] + 
     D[C2[x, y]*u[x, y], y] + (1/(2 NN))*
      D[C3[x, y]*u[x, y], x, x] - (1/NN)*
      D[C4[x, y]*u[x, y], x, y] + (1/(2*NN))*
      D[C5[x, y]*u[x, y], y, y] == 0,
   DirichletCondition[u[x, y] == 1, x == 1 && y == 0]}, 
  u, {x, y} \[Element] region]

This gives a message:

enter image description here

If the equation and setup are indeed correct then there are two ways to keep the convection dominance of this at bay. One is to add artificial diffusion which will change the PDE and the other is to mesh it better. I have gone with the second approach; and because I did not want to think about the PDE too much I have chosen the mesh-it-do-death approach.

uval = NDSolveValue[{D[C1[x, y]*u[x, y], x] + 
      D[C2[x, y]*u[x, y], y] + (1/(2 NN))*
       D[C3[x, y]*u[x, y], x, x] - (1/NN)*
       D[C4[x, y]*u[x, y], x, y] + (1/(2*NN))*
       D[C5[x, y]*u[x, y], y, y] == 0,
    DirichletCondition[u[x, y] == 1, x == 1 && y == 0]}, 
   u, {x, y} \[Element] region, 
   Method -> {"PDEDiscretization" -> {"FiniteElement", 
       "MeshOptions" -> {"MaxCellMeasure" -> 0.0001}}}];

Plot3D[uval[x, y], {x, y} \[Element] region, PlotRange -> All]

enter image description here

When we check the integral constraint we see that:

NIntegrate[uval[x, y], {x, y} \[Element] region]
NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.
(* 0.000345862 *)

which is not the requested value. How could it be, we did not improve a constraint. Now, currently NDSolve can not handle integral constraints but many moons back I gave a presentation where I showed how to do it. I am going to reproduce the relevant part here.

Parse and discretize the PDE:

Needs["NDSolve`FEM`"]
{dpde, dbc, vd, sd, md} = 
  ProcessPDEEquations[{D[C1[x, y]*u[x, y], x] + 
      D[C2[x, y]*u[x, y], y] + (1/(2 NN))*
       D[C3[x, y]*u[x, y], x, x] - (1/NN)*
       D[C4[x, y]*u[x, y], x, y] + (1/(2*NN))*
       D[C5[x, y]*u[x, y], y, y] == 0,
    DirichletCondition[u[x, y] == 1, x == 1 && y == 0]}, 
   u, {x, y} \[Element] region, 
   Method -> {"PDEDiscretization" -> {"FiniteElement", 
       "MeshOptions" -> {"MaxCellMeasure" -> 0.0000125}}}];

Extract the system matrices and deploy the given boundary conditions:

l = dpde["LoadVector"];
s = dpde["StiffnessMatrix"];
DeployBoundaryConditions[{l, s}, dbc];

Extract the mesh and the degrees of freedom:

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

Next, we write a little helper function to create a boundary condition to satisfy boundary integral constrains:

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]
  ]

This function works by creating a weighted value for every node in the mesh in such a way that once these are summed up they satisfy the boundary integral value.

Make use of that:

dc = mkIntegralConstraint[1, mesh, dof]

Now, we deploy that boundary condition in such that the values are appended to the system matrices as Lagrangian multipliers.

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

Solve the equations and post process the equations:

sol = LinearSolve[s, l];
if2 = ElementMeshInterpolation[{mesh}, Take[sol, dof]];
Plot3D[if2[x, y], {x, y} \[Element] region, PlotRange -> All]

enter image description here

Check the integral:

NIntegrate[if2[x, y], {x, y} \[Element] region, 
 Method -> "FiniteElement"]

enter image description here

(* 0.990778 *)

This is not perfect but that's do to the nature of the PDE used that is convection dominated and unstable. In that light, I think, the result is OK.

$\endgroup$
9
  • $\begingroup$ Hi, Thanks you very much for your response. I'll have a proper read through your answer, and compare it to an approximation I have. Thanks! Edit: Actually, in the example you gave, you've used Dirichlet boundary conditions. The majority of my boundary conditions I want to be neumann, as in the second exmaple in my post (at least, that's what I've tried to implement). Would you be able to show me how to put in the neumann conditions I require? (Is this what you meant about the equation not being quite right?) $\endgroup$
    – user112495
    Jun 26 '20 at 10:02
  • $\begingroup$ I've tried running your code as is. It produces the first plot, as it did for you. However, when I get to the point of running the last section, I get the following error messages: 1) Take::seqs: Sequence specification (+n, -n, {+n}, {-n}, {m, n}, or {m, n, s}) expected at position 2 in Take[LinearSolve[dpde[StiffnessMatrix],dpde[LoadVector]],md[DegreesOfFreedom]]. 2) Plot3D::idomdim: {x,y}[Element]region does not have a valid dimension as a plotting domain. $\endgroup$
    – user112495
    Jun 26 '20 at 11:08
  • $\begingroup$ Also, thinking about it a bit more, I might want the absorbing boundary condition to be u[1,0]=0,rather than u[1,0]=1. $\endgroup$
    – user112495
    Jun 26 '20 at 12:36
  • $\begingroup$ Do you have any thoughts on the issues I'm currently having? $\endgroup$
    – user112495
    Jul 3 '20 at 13:16
  • 1
    $\begingroup$ @user112495, make sure that you re-add definitions for NN and region. $\endgroup$
    – user21
    Jul 16 '20 at 10:18
1
$\begingroup$

Infered from the NeumannValue settings the region is

RegionPlot[
 x + y >= 1 \[Or] x <= 0 \[Or] y <= 0, {x, -.1, 1.1}, {y, -.1, 1.1}]

region of the given question

This can be enhanced using PlotPoints option.

As Region it looks:

Region[MeshRegion[{{0, 0}, {1, 0}, {0, 1}}, Line[{{1, 2, 3, 1}}]]]

region in the given case

Error message of the original question is than: enter image description here

Better work with

polygon = Region[Polygon[{{0, 0}, {1, 0}, {0, 1}}]]

region of the question

From the fundamentals of mathemtics follows that without restriction all C_i functions can be set equal one.

Then the questions looks like this:

{dpde, dbc, vd, sd, md} = 
 ProcessPDEEquations[{D[u[x, y], x] + D[u[x, y], y] + 
     D[u[x, y], x, x] + D[u[x, y], x, y] + D[u[x, y], y, y] == 
    NeumannValue[0, x <= 0] + NeumannValue[0, y <= 0] + 
     NeumannValue[0, x + y >= 1], 
   DirichletCondition[u[x, y] == 0, x == 1 && y == 0]}, 
  u, {x, y} \[Element] polygon, 
  Method -> {"PDEDiscretization" -> {"FiniteElement", 
      "MeshOptions" -> {"MaxCellMeasure" -> 0.000125}}}]
l = dpde["LoadVector"];
s = dpde["StiffnessMatrix"];
DeployBoundaryConditions[{l, s}, dbc];
mesh = md["ElementMesh"];
dof = md["DegreesOfFreedom"];

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]]

dc = mkIntegralConstraint[1, mesh, dof]
DeployBoundaryConditions[{l, s}, dc, "ConstraintMethod" -> "Append"];

sol = LinearSolve[s, l];
if2 = ElementMeshInterpolation[{mesh}, Take[sol, dof]];
Plot3D[if2[x, y], {x, y} \[Element] region, PlotRange -> All]

output

Plot3D graphics

This confirms the expection of the question in all manners:

MaxValue[if2[x, y], {x, y} \[Element] polygon]
if2[0, 1]

2.14897

2.14897

The minimum is in a corner too:

if2[1, 0]

-2.74605*10^-7

These confirm that the given boundary conditions hold only for the Dirichlet condition!

From the documentation of DirichletCondition

Typically, at least one Dirichlet-type boundary condition needs to be specified to make the differential equation uniquely solvable. Dirichlet conditions are also called essential boundary conditions.

Something different from the coordinate function is definitely wrong after this changes.

$\endgroup$

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.