2
$\begingroup$

I am using MMA 12.0.0 on a Mac running OS version 10.13.6.

Before attempting my own problem I am attempting to follow exactly the code in the FEM Documentation Tutorial, as an exercise:

https://reference.wolfram.com/language/FEMDocumentation/tutorial/SolvingPDEwithFEM.html#1405680271

My version (after its first two lines) is supposed to be identical to that in the above tutorial as follows:

Clear["Global`*"];
NotebookSave[];
Needs["NDSolve`FEM`"];
rules = {length -> 22/10, hight -> 41/100};
\[CapitalOmega] = 
  RegionDifference[Rectangle[{0, 0}, {length, hight}], 
    Disk[{1/5, 1/5}, 1/20]] /. rules;
region = RegionPlot[\[CapitalOmega], AspectRatio -> Automatic]
op = {\!\(
\*SubscriptBox[\(\[Del]\), \({x, 
        y}\)] . \(({{\(-\[Mu]\), 0}, {0, \(-\[Mu]\)}} . 
\*SubscriptBox[\(\[Del]\), \({x, y}\)]u[x, y])\)\) + \[Rho] {{u[x, y],
          v[x, y]}}.\!\(
\*SubscriptBox[\(\[Del]\), \({x, y}\)]\(u[x, y]\)\) + 
     Derivative[1, 0][p][x, y], \!\(
\*SubscriptBox[\(\[Del]\), \({x, 
        y}\)] . \(({{\(-\[Mu]\), 0}, {0, \(-\[Mu]\)}} . 
\*SubscriptBox[\(\[Del]\), \({x, y}\)]v[x, y])\)\) + \[Rho] {{u[x, y],
          v[x, y]}}.\!\(
\*SubscriptBox[\(\[Del]\), \({x, y}\)]\(v[x, y]\)\) + 
     Derivative[0, 1][p][x, y], 
    Derivative[1, 0][u][x, y] + 
     Derivative[0, 1][v][x, y]} /. {\[Mu] -> 10^-3, \[Rho] -> 1};
pde = op == {0, 0, 0};
bcs = {DirichletCondition[u[x, y] == 4*0.3*y (hight - y)/hight^2, 
     x == 0], DirichletCondition[v[x, y] == 0, x == 0], 
    DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, 
     x > 0 && x < length], 
    DirichletCondition[p[x, y] == 0., x == length]} /. rules;
refinementRegion = 
  ImplicitRegion[
    a^4 (23 (y - 1/5))^2 + 
      b^2 (5/2 x - 2.05)^3 (2 a + (5/2 x - 2.05)) < 
     0, {{x, 0, length}, {y, 0, hight}}] /. 
   Flatten[{rules, a -> 1, b -> 5/2}];
Show[RegionPlot[\[CapitalOmega]], RegionPlot[refinementRegion], 
 AspectRatio -> Automatic]
mrf = With[{rmf = RegionMember[refinementRegion]}, 
   Function[{vertices, area}, Block[{x, y}, {x, y} = Mean[vertices];
     If[rmf[{x, y}], area > 0.00025, area > 0.0025]]]];
{xVel, yVel, pressure} = 
  NDSolveValue[{pde, bcs}, {u, v, p}, 
   Element[{x, y}, \[CapitalOmega]], 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}, 
     "MeshOptions" -> {"IncludePoints" -> {{0.15, 0.2}, {0.25, 0.2}}, 
       AccuracyGoal -> 5, PrecisionGoal -> 5, 
       MeshRefinementFunction -> mrf}}];

After printing out the regions correctly, it returns the following warning/error messages:

NDSolveValue: The maximum derivative order of the nonlinear PDE coefficients for the Finite Element Method is larger than 1. It may help to rewrite the PDE in inactive form.

followed by a message that the two sides of my last command are not the same shape.

Perhaps I transcribed something wrong from the tutorial, but I can't see where. Otherwise, why does the tutorial not work in my version? I thought FEM was supposed to work in MMA 12?

Thanks.

$\endgroup$
  • $\begingroup$ When you evaluate the code in the tutorial does that work? $\endgroup$ – user21 Feb 14 at 11:49
  • $\begingroup$ @user21. Actually yes! So, a transcription error. Sorry, I should have just cut and pasted from the start. When I cut-and-paste, many of the differential operators turn up as Inactive in my notebook. But those commands are "invisible" in the tutorial (perhaps that could be fixed in the tutorial, to prevent some other user having the same issue as I had)? Thanks for your help. $\endgroup$ – Paul Harrison Feb 14 at 11:59
  • $\begingroup$ Sorry for the trouble, what do you mean by 'fixed' : Do you mean that the equations should be typeset Inactive[Div][...] instead of what we have now? I mean I understand your point but how would I make that explicit in all the places where there are Inactive PDEs. Writing Inactive[...] everywhere will not be approved. A bit like what follows in FEMDocumentation/tutorial/SolvingPDEwithFEM#156906496. Also note this section: FEMDocumentation/tutorial/SolvingPDEwithFEM#205761584. $\endgroup$ – user21 Feb 14 at 12:07
  • $\begingroup$ Thanks for helping! For some reason, we are not fully understanding each other. What I meant by "fixed" was to make it explicit somewhere in the tutorial that you have to set the Div operators Inactive. The cut-and-paste from the tutorial in my .nb, show-up as: "Inactive[Div]", which is, I guess, what you typed in your notebook. Why is the command "Inactive" suppressed in the code block in the tutorial? The typical user will not know that they are actually set Inactive. If you need to save space, you could just write below the code that all the Del operators need to be set Inactive. $\endgroup$ – Paul Harrison Feb 14 at 12:24
  • $\begingroup$ OK, let's try to improve this. If you paste {Inactive[ Div][{{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][ u[x, y], {x, y}], {x, y}] + Derivative[1, 0][p][x, y], Inactive[Div][{{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][ v[x, y], {x, y}], {x, y}] + Derivative[0, 1][p][x, y], Derivative[0, 1][v][x, y] + Derivative[1, 0][u][x, y]} into your notebook if will format as a proper equation. You can use that formatted equation as input to NDSolve. Are you saying that it is better in the notebook to use the Inactive InputForm that the formatted equation. $\endgroup$ – user21 Feb 14 at 12:38