FEM Mesh related errors

Question

Bug introduced in 12.0. Fixed in 12.1.0.

I try to solve the simple heat transfer equation but there are some errors at the stage of the pre-processing.

Needs["NDSolve`FEM`"]
c = 1380*0.6 + 4200*0.4;
ρ = 1250.*0.6 + 1000*0.4;
λ = (0.111*0.6 + 0.56*0.4 );
Q = 3000000;
σ = 5.675*10^-8;
ϵ = 0.33;
T0 = 300.;
R0 = 0.00007;  
rr = 0.0005;
h = 0.001;
MaxT=5;

Ω = ToElementMesh[Rectangle[{0., 0.}, {rr, h}], 
   MaxCellMeasure -> 0.5 10^-10, MeshQualityGoal -> 1];

The mesh is ElementMesh[{{0., 0.0005}, {0., 0.001}}, {QuadElement["<" 5000 ">"]}] with 5000 elements.

op = D[u[t, r, z],t] - λ/(c ρ)*(D[r^2 D[u[t, r, z], r], r]/r^2 + 
      D[u[t, r, z], {z, 2}]);
Γ = NeumannValue[
    Piecewise[{{Q - ϵ σ (u[t, r, z]^4 - T0^4), 
       0 <= r <= R0}, {-ϵ σ (u[t, r, z]^4 - T0^4), 
       R0 < r <= rr}}], z == 0] +
   NeumannValue[-ϵ σ (u[t, r, z]^4 - T0^4), z == h] +
   NeumannValue[-ϵ σ (u[t, r, z]^4 - T0^4), r == rr];

fun=NDSolveValue[{op == Γ, 
  u[0, r, z] == T0}, u, {t, 0, 
  MaxT}, {r, z} ∈ Ω, 
 EvaluationMonitor :> (tp = t), 
 Method -> {"PDEDiscretization" -> {"MethodOfLines", 
"SpatialDiscretization" -> {"FiniteElement","IntegrationOrder" -> 4}}}, MaxSteps -> 100000, 
 MaxStepFraction -> 0.01]

Running the NDSolveValue leads to generation of a lot of errors:

Part::partw: Part {10154,10155,30390} of {} does not exist.

Part::partw: Part {10155,10156,30392} of {} does not exist.

Part::partw: Part {10156,10157,30394} of {} does not exist.

General::stop: Further output of Part::partw will be suppressed during this calculation.

Dot::dotsh: Tensors {NDSolve`FEM`FEMBoundaryConditionsDump`IntegratedShapeFunctionLookup[142]} and {{}[[{10154,10155,30390}]],{}[[{10155,10156,30392}]],{}[[{10156,10157,30394}]],{}[[{10157,10158,30396}]],{}[[{10158,10159,30398}]],{}[[{10159,10160,30400}]],{}[[{10160,10161,30402}]],{}[[{10161,10162,30404}]],{}[[{10162,10163,30406}]],{}[[{10163,10164,30408}]],{}[[{10164,10165,30410}]],{}[[{10165,10166,30412}]],{}[[{10166,10167,30414}]],<<26>>,{}[[{10193,10194,30468}]],{}[[{10194,10195,30470}]],{}[[{10195,10196,30472}]],{}[[{10196,10197,30474}]],{}[[{10197,10198,30476}]],{}[[{10198,10199,30478}]],{}[[{10199,10200,30480}]],{}[[{10200,10201,30482}]],{}[[{10201,10202,30484}]],{}[[{10202,10203,30486}]],{}[[{10203,10204,30488}]],<<92>>} have incompatible shapes.

CompiledFunction::cfta: Argument {NDSolve`FEM`FEMBoundaryConditionsDump`IntegratedShapeFunctionLookup[142]} at position 1 should be a rank 2 tensor of machine-size real numbers.

Transpose::nmtx: The first two levels of {{NDSolve`FEM`FEMBoundaryConditionsDump`IntegratedShapeFunctionLookup[142]}.{{}[[{10154,10155,30390}]],{}[[{10155,10156,30392}]],{}[[{10156,10157,30394}]],{}[[{10157,10158,30396}]],{}[[{10158,10159,30398}]],{}[[{10159,10160,30400}]],{}[[{10160,10161,30402}]],{}[[{10161,10162,30404}]],{}[[{10162,10163,30406}]],{}[[{10163,10164,30408}]],{}[[{10164,10165,30410}]],<<29>>,{}[[{10194,10195,30470}]],{}[[{10195,10196,30472}]],{}[[{10196,10197,30474}]],{}[[{10197,10198,30476}]],{}[[{10198,10199,30478}]],{}[[{10199,10200,30480}]],{}[[{10200,10201,30482}]],{}[[{10201,10202,30484}]],{}[[{10202,10203,30486}]],{}[[{10203,10204,30488}]],<<92>>}} cannot be transposed.

CompiledFunction::cfta: Argument Transpose[{{NDSolve`FEM`FEMBoundaryConditionsDump`IntegratedShapeFunctionLookup[142]}.{{}[[{10154,10155,30390}]],{}[[{10155,10156,30392}]],{}[[{10156,10157,30394}]],{}[[{10157,10158,30396}]],{}[[{10158,10159,30398}]],{}[[{10159,10160,30400}]],{}[[{10160,10161,30402}]],{}[[{10161,10162,30404}]],{}[[{10162,10163,30406}]],{}[[{10163,10164,30408}]],{}[[{10164,10165,30410}]],<<29>>,{}[[{10194,10195,30470}]],{}[[{10195,10196,30472}]],{}[[{10196,10197,30474}]],{}[[{10197,10198,30476}]],{}[[{10198,10199,30478}]],{}[[{10199,10200,30480}]],{}[[{10200,10201,30482}]],{}[[{10201,10202,30484}]],{}[[{10202,10203,30486}]],{}[[{10203,10204,30488}]],<<92>>}}] at position 1 should be a rank 2 tensor of machine-size real numbers.

Dot::dotsh: Tensors {NDSolve`FEM`FEMBoundaryConditionsDump`IntegratedShapeFunctionLookup[142]} and {{}[[{10154,10155,30390}]],{}[[{10155,10156,30392}]],{}[[{10156,10157,30394}]],{}[[{10157,10158,30396}]],{}[[{10158,10159,30398}]],{}[[{10159,10160,30400}]],{}[[{10160,10161,30402}]],{}[[{10161,10162,30404}]],{}[[{10162,10163,30406}]],{}[[{10163,10164,30408}]],{}[[{10164,10165,30410}]],{}[[{10165,10166,30412}]],{}[[{10166,10167,30414}]],<<26>>,{}[[{10193,10194,30468}]],{}[[{10194,10195,30470}]],{}[[{10195,10196,30472}]],{}[[{10196,10197,30474}]],{}[[{10197,10198,30476}]],{}[[{10198,10199,30478}]],{}[[{10199,10200,30480}]],{}[[{10200,10201,30482}]],{}[[{10201,10202,30484}]],{}[[{10202,10203,30486}]],{}[[{10203,10204,30488}]],<<92>>} have incompatible shapes.

Dot::dotsh: Tensors {NDSolve`FEM`FEMBoundaryConditionsDump`IntegratedShapeFunctionLookup[142]} and {{}[[{10154,10155,30390}]],{}[[{10155,10156,30392}]],{}[[{10156,10157,30394}]],{}[[{10157,10158,30396}]],{}[[{10158,10159,30398}]],{}[[{10159,10160,30400}]],{}[[{10160,10161,30402}]],{}[[{10161,10162,30404}]],{}[[{10162,10163,30406}]],{}[[{10163,10164,30408}]],{}[[{10164,10165,30410}]],{}[[{10165,10166,30412}]],{}[[{10166,10167,30414}]],<<26>>,{}[[{10193,10194,30468}]],{}[[{10194,10195,30470}]],{}[[{10195,10196,30472}]],{}[[{10196,10197,30474}]],{}[[{10197,10198,30476}]],{}[[{10198,10199,30478}]],{}[[{10199,10200,30480}]],{}[[{10200,10201,30482}]],{}[[{10201,10202,30484}]],{}[[{10202,10203,30486}]],{}[[{10203,10204,30488}]],<<92>>} have incompatible shapes.

General::stop: Further output of Dot::dotsh will be suppressed during this calculation.

CompiledFunction::cfta: Argument {{{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}},{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}},{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}}},{{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}},{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}},{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}}},{{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}},{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}},{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}}},<<46>>,{{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}},{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}},{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}}},<<92>>} at position 1 should be a rank 3 tensor of machine-size real numbers.

General::stop: Further output of CompiledFunction::cfta will be suppressed during this calculation.

AssembleMatrix::badmat: {7.04225*10^-6 MapThreadDot[{{{0.687298},{-0.0872983},{0.4}},{{-0.0872983},{0.687298},{0.4}},{{0.},{0.},{1.}}},{{{-1.87275*10^-8 Plus[<<2>>]}},{{-1.87275*10^-8 Plus[<<2>>]}},{{-1.87275*10^-8 Plus[<<2>>]}}}],7.04225*10^-6 MapThreadDot[{{{0.687298},{-0.0872983},{0.4}},{{-0.0872983},{0.687298},{0.4}},{{0.},{0.},{1.}}},{{{-1.87275*10^-8 Plus[<<2>>]}},{{-1.87275*10^-8 Plus[<<2>>]}},{{-1.87275*10^-8 Plus[<<2>>]}}}],<<48>>,<<92>>} is not a valid element matrix to be assembled.

AssembleMatrix::badmat: {7.04225*10^-6 MapThreadDot[{{{0.687298},{-0.0872983},{0.4}},{{-0.0872983},{0.687298},{0.4}},{{0.},{0.},{1.}}},{MapThreadDot[{{{<<1>>}},{{<<1>>}},{{<<1>>}}},{{{<<3>>}},{{<<3>>}},{{<<3>>}}}],MapThreadDot[{{{<<1>>}},{{<<1>>}},{{<<1>>}}},{{{<<3>>}},{{<<3>>}},{{<<3>>}}}],MapThreadDot[{{{<<1>>}},{{<<1>>}},{{<<1>>}}},{{{<<3>>}},{{<<3>>}},{{<<3>>}}}],<<45>>,MapThreadDot[{{{<<1>>}},{{<<1>>}},{{<<1>>}}},{{{<<3>>}},{{<<3>>}},{{<<3>>}}}],MapThreadDot[{{{<<1>>}},{{<<1>>}},{{<<1>>}}},{{{<<3>>}},{{<<3>>}},{{<<3>>}}}],<<92>>}],7.04225*10^-6 <<1>>,<<47>>,7.04225*10^-6 <<1>>,<<92>>} is not a valid element matrix to be assembled.

AssembleMatrix::badmat: {7.04225*10^-6 MapThreadDot[{{{0.687298},{-0.0872983},{0.4}},{{-0.0872983},{0.687298},{0.4}},{{0.},{0.},{1.}}},{MapThreadDot[{{{<<1>>}},{{<<1>>}},{{<<1>>}}},{{{<<3>>}},{{<<3>>}},{{<<3>>}}}],MapThreadDot[{{{<<1>>}},{{<<1>>}},{{<<1>>}}},{{{<<3>>}},{{<<3>>}},{{<<3>>}}}],MapThreadDot[{{{<<1>>}},{{<<1>>}},{{<<1>>}}},{{{<<3>>}},{{<<3>>}},{{<<3>>}}}],<<45>>,MapThreadDot[{{{<<1>>}},{{<<1>>}},{{<<1>>}}},{{{<<3>>}},{{<<3>>}},{{<<3>>}}}],MapThreadDot[{{{<<1>>}},{{<<1>>}},{{<<1>>}}},{{{<<3>>}},{{<<3>>}},{{<<3>>}}}],<<92>>}],7.04225*10^-6 <<1>>,<<47>>,7.04225*10^-6 <<1>>,<<92>>} is not a valid element matrix to be assembled.

General::stop: Further output of AssembleMatrix::badmat will be suppressed during this calculation.

Transpose::nmtx: The first two levels of {{NDSolve`FEM`FEMBoundaryConditionsDump`IntegratedShapeFunctionLookup[71]}.{{}[[{286,143,10722}]],{}[[{429,286,11008}]],{}[[{572,429,11293}]],{}[[{715,572,11578}]],{}[[{858,715,11863}]],{}[[{1001,858,12148}]],{}[[{1144,1001,12433}]],{}[[{1287,1144,12718}]],{}[[{1430,1287,13003}]],{}[[{1573,1430,13288}]],{}[[{1716,1573,13573}]],{}[[{1859,1716,13858}]],{}[[{2002,1859,14143}]],<<26>>,{}[[{5863,5720,21838}]],{}[[{6006,5863,22123}]],{}[[{6149,6006,22408}]],{}[[{6292,6149,22693}]],{}[[{6435,6292,22978}]],{}[[{6578,6435,23263}]],{}[[{6721,6578,23548}]],{}[[{6864,6721,23833}]],{}[[{7007,6864,24118}]],{}[[{7150,7007,24403}]],{}[[{7293,7150,24688}]],<<21>>}} cannot be transposed.

Transpose::nmtx: The first two levels of {{NDSolve`FEM`FEMBoundaryConditionsDump`IntegratedShapeFunctionLookup[71]}.{{}[[{1,144,10297}]],{}[[{144,287,10724}]],{}[[{287,430,11009}]],{}[[{430,573,11294}]],{}[[{573,716,11579}]],{}[[{716,859,11864}]],{}[[{859,1002,12149}]],{}[[{1002,1145,12434}]],{}[[{1145,1288,12719}]],{}[[{1288,1431,13004}]],{}[[{1431,1574,13289}]],{}[[{1574,1717,13574}]],<<27>>,{}[[{5578,5721,21554}]],{}[[{5721,5864,21839}]],{}[[{5864,6007,22124}]],{}[[{6007,6150,22409}]],{}[[{6150,6293,22694}]],{}[[{6293,6436,22979}]],{}[[{6436,6579,23264}]],{}[[{6579,6722,23549}]],{}[[{6722,6865,23834}]],{}[[{6865,7008,24119}]],{}[[{7008,7151,24404}]],<<21>>}} cannot be transposed.

General::stop: Further output of Transpose::nmtx will be suppressed during this calculation.

NDSolveValue::indexss: The DAE solver failed at t = 0.`. The solver is intended for index 1 DAE systems and structural analysis indicates that the DAE is structurally singular.

As I understood, the initial errors occur because the solver tries to find the extra mesh elements above the 5000 of the basically produced. And the following errors are derived from these..

Can anybody help with a proper problem initialization?

The solver finally produces certain solution but I'm not sure that it is correct..

Answer 1

This is a bug with time dependent nonlinear NeumannValue introduced in Version 12.0. This is fixed in version 12.1 which is hopefully coming in the not to distant future at the time of this writing. There is no known top level workaround. I apologize for the inconvenience. If this is very important for you, I cloud try to implement this with the low level FEM functions - I have not tried this so I can no guarantee that this will work. Let me know what you think.

Update:

Here is a low level code.

Needs["NDSolve`FEM`"]
c = 1380*0.6 + 4200*0.4;
ρ = 1250.*0.6 + 1000*0.4;
λ = (0.111*0.6 + 0.56*0.4);
Q = 3000000;
σ = 5.675*10^-8;
ϵ = 0.33;
T0 = 300.;
R0 = 0.00007;
rr = 0.0005;
h = 0.001;
MaxT = 5;

mesh = ToElementMesh[Rectangle[{0., 0.}, {rr, h}], 
  MaxCellMeasure -> 10 10^-10, MeshQualityGoal -> 1, 
  "MeshElementType" -> "TriangleElement", "MeshOrder" -> 1]

Note that I have substantially coarsened the mesh. Try to fix the remaining issues first before switching back to a finer mesh.

Set up the PDE without the boundary conditions:

op = D[u[t, r, z], 
    t] - λ/(c ρ)*(D[r^2 D[u[t, r, z], r], r]/r^2 + 
      D[u[t, r, z], {z, 2}]);

The BCs have some issues: Split the the Piecewise into two BCs. For some reason the second BC introduces a real slow convergence. You'd need to experiment a bit with that.

Γ = {(*NeumannValue[
   Piecewise[{{Q-ϵ σ (u[t,r,z]^4-T0^4),
   0\[LessEqual]r\[LessEqual]R0},{-ϵ σ (u[t,r,z]^4-
   T0^4),R0<r\[LessEqual]rr}}],
   z\[Equal]0],*)(*NeumannValue[-ϵ σ (u[t,r,z]^4-
   T0^4),z\[Equal]h],*)
   NeumannValue[-ϵ σ (u[t, r, z]^4 - T0^4), r == rr]};

Set up FEM data:

{sdpde, sdbc, vd, sd, methodData} = 
  NDSolve`FEM`ProcessPDEEquations[{op == 0, u[0, r, z] == T0}, 
   u, {t, 0, MaxT}, {r, z} \[Element] mesh];

Now, we initialize the BCs separately:

initBCs = InitializeBoundaryConditions[vd, sd, {Γ}];
sbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd];

Set up a helper function to apply the BCs during time integration:

discretizePDEResidual[t_?NumericQ, u_?VectorQ, dudt_?VectorQ] :=

 Module[{l, s, d, tdpde, tbcs, nldpde, nlbcs, sdTemp},

  NDSolve`SetSolutionDataComponent[sd, "Time", t];
  NDSolve`SetSolutionDataComponent[sd, "DependentVariables", u];

  l = sdpde["LoadVector"];
  s = sdpde["StiffnessMatrix"];
  d = sdpde["DampingMatrix"];

  tbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd, 
    "Transient"];

  nlbcs = 
   DiscretizeBoundaryConditions[initBCs, methodData, sd, 
    "Nonlinear"];

  DeployBoundaryConditions[{l, s, d}, nlbcs];
  DeployBoundaryConditions[{l, s, d}, tbcs];
  DeployBoundaryConditions[{l, s, d}, sbcs];

  d.dudt + s.u - l
  ]

Set up the initial conditions and the sparsity pattern:

initT0 = T0 & /@ mesh["Coordinates"];
sparsity = sdpde["DampingMatrix"]["PatternArray"];

Do the time integration. Because there were may step rejections (with the other BCs) I added the "IDA" option to reduce those. This might not be necessary once the NeumannValue issues are understood.

Monitor[tufun = 
  NDSolveValue[{discretizePDEResidual[t, u[t], u'[ t]] == 0, 
    u[0] == initT0}, u, {t, 0, MaxT}
   , Method -> {
     "TimeIntegration" -> {"IDA", "MaxDifferenceOrder" -> 2}
     , "EquationSimplification" -> "Residual"}
   , Jacobian -> {Automatic, Sparse -> sparsity}
   , EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])
   (*,MaxStepFraction\[Rule]0.01*)
   ], monitor]

From the result you can then construct an interpolating function with:

ufun = ElementMeshInterpolation[{tufun["Coordinates"][[1]], 
   methodData["ElementMesh"]}, Partition[tufun["ValuesOnGrid"], 1]]

Answer 2

We can use a linear FEM solver. It is necessary to normalize u by T0, Q and $\sigma u^4$ by $\rho c T0$. Then the code is

Needs["NDSolve`FEM`"]
c = 1380*0.6 + 4200*0.4;
\[Rho] = 1250.*0.6 + 1000*0.4;
\[Lambda] = (0.111*0.6 + 0.56*0.4); T0 = 300.;
Q = 3000000/(c \[Rho] T0);
qr = 0.33 T0^3 (5.675*10^-8)/(c \[Rho] );
k = \[Lambda]/(c \[Rho] );

R0 = 0.00007;
rr = 0.0005;
h = 0.001;
tm = 5 ; tau = 1/20; nmax = Round[tm/tau];
\[CapitalOmega] = 
 ToElementMesh[Rectangle[{0., 0.}, {rr, h}], 
  MeshRefinementFunction -> 
   Function[{vertices, area}, 
    area > 0.0000005 (0.00001 + 10 Norm[Mean[vertices]])]];
\[CapitalOmega]["Wireframe"]
U[0][r_, z_] := 1; Do[
 U[i] = NDSolveValue[{(u[r, z] - U[i - 1][r, z])/tau - 
       k (D[u[r, z], r, r] + 2 D[u[r, z], r]/r + 
          D[u[r, z], {z, 2}]) == (NeumannValue[
         Q - qr  (U[i - 1][r, z]^3 u[r, z] - 1), 
         z == 0 && 0 < r <= R0] + 
        NeumannValue[-qr (U[i - 1][r, z]^3 u[r, z] - 1), 
         z == 0 && R0 < r <= rr] + 
        NeumannValue[-qr  (U[i - 1][r, z]^3 u[r, z] - 1), z == h] + 
        NeumannValue[-qr (U[i - 1][r, z]^3 u[r, z] - 1), r == rr])}, 
    u, {r, z} \[Element] \[CapitalOmega]];, {i, 1, nmax}];

Temperature visualization

{DensityPlot[T0 U[nmax][r, z], {r, z} \[Element] \[CapitalOmega], 
  ColorFunction -> "TemperatureMap", PlotLegends -> Automatic, 
  PlotRange -> All, FrameLabel -> Automatic], 
 DensityPlot[T0 U[nmax][r, z], {r, z} \[Element] \[CapitalOmega], 
  ColorFunction -> "TemperatureMap", PlotLegends -> Automatic, 
  FrameLabel -> Automatic], ListLinePlot[Table[{i tau 1., T0 U[i][.0, .0]}, {i, 0, nmax}], 
 AxesLabel -> {t, T}, PlotRange -> All]}