LinearSolve problem occurs as solving PDEs with NDSolve

Question

I am new to Mathematica, please forgive me asking naive questions. I tried to solve PDEs numerically using NDSolve, but failed to go through due to errors. Two of three PDEs are time-dependent and another one is time-independent (or implicitly time-dependent). Here is the code:

Needs["NDSolve`FEM`"]
dr = 0.5; rx = 2.5/2; xs = -(dr + rx); ys = 0.0;
ec = RegionUnion[Disk[{-rx, 0}, dr], Rectangle[{-rx, -dr}, {rx, dr}], 
   Disk[{rx, 0}, dr]];
Region[ec]
bmesh = ToBoundaryMesh[ec, "MaxBoundaryCellMeasure" -> .05];
bmesh["Wireframe"]
f = Function[{vertices, area}, 
   area > 0.004 (0.9 - 0.5 Norm[Mean[vertices]])];
(mesh = ToElementMesh[bmesh, MeshRefinementFunction -> f])["Wireframe"]
u0 = 11.2; v0 = 1; d0 = 0.97; V2 = -0.2;
Eqs = {D[u[t, x, y], t] == 
    d0 (Laplacian[u[t, x, y], {x, y}] + w[t, x, y]),
   D[v[t, x, y], t] == Laplacian[v[t, x, y], {x, y}] - w[t, x, y],
   Laplacian[
     w[t, x, y], {x, y}] == -4 \[Pi] (v[t, x, y] - u[t, x, y])};
ics = {u[0, x, y] == u0, v[0, x, y] == v0, w[0, x, y] == -0.0};
bcs = DirichletCondition[w[t, x, y] == V2, {x, y} \[Element] bmesh];
{usol, vsol, wsol} = 
  NDSolveValue[{Eqs, ics, bcs}, {u, v, w}, {t, 0, 
    10^6}, {x, y} \[Element] mesh];

The result shows errors

LinearSolve::parpiv: Zero pivot was detected during the numerical factorization or there was a problem in the iterative refinement process. It is possible that the matrix is ill-conditioned or singular.
Set::shape: Lists {usol,vsol,wsol} and NDSolveValue[<<1>>] are not the same shape.
LinearSolve::parpiv: Zero pivot was detected during the numerical factorization or there was a problem in the iterative refinement process. It is possible that the matrix is ill-conditioned or singular.

Totally have no ideas where errors come from. Thanks for any instructions and suggestions.

Answer 1

There are several problems discussed here and here, and we can resolve these problems by adding small regularization term to the last equation as follows

Needs["NDSolve`FEM`"]
dr = 0.5; rx = 2.5/2; xs = -(dr + rx); ys = 0.0;
ec = RegionUnion[Disk[{-rx, 0}, dr], Rectangle[{-rx, -dr}, {rx, dr}], 
   Disk[{rx, 0}, dr]];
Region[ec]
bmesh = ToBoundaryMesh[ec, "MaxBoundaryCellMeasure" -> .05];
bmesh["Wireframe"]
f = Function[{vertices, area}, 
   area > 0.004 (0.9 - 0.5 Norm[Mean[vertices]])];
(mesh = ToElementMesh[bmesh, MeshRefinementFunction -> f])["Wireframe"]
u0 = 11.2; v0 = 1; d0 = 0.97; V2 = -0.2;
Eqs = Inactivate[{D[u[t, x, y], t] - 
     d0 (Laplacian[u[t, x, y], {x, y}] + w[t, x, y]), 
    D[v[t, x, y], t] - Laplacian[v[t, x, y], {x, y}] + w[t, x, y], 
    0.001 D[w[t, x, y], t] - Laplacian[w[t, x, y], {x, y}] - 
     4 \[Pi] (v[t, x, y] - u[t, x, y])}, Laplacian | D];
ics = {u[0, x, y] == u0, v[0, x, y] == v0, w[0, x, y] == 0.0};
bcs = DirichletCondition[w[t, x, y] == V2, True];
{usol, vsol, wsol} = 
 NDSolveValue[{Activate[Eqs] == {NeumannValue[0, True], 
     NeumannValue[0, True], 0}, ics, bcs}, {u, v, 
   w}, {x, y} \[Element] mesh, {t, 0, 10}];

Note, in this case we don't need to integrate up to t=10^6, even t=5 could be fine. Visualization

Plot[{usol[t, 0, 0], vsol[t, 0, 0], wsol[t, 0, 0]}, {t, 0, 10}, 
 PlotLegends -> {u, v, w}, Frame -> True, FrameLabel -> Automatic, 
 PlotRange -> All]

{DensityPlot[usol[10, x, y], {x, y} \[Element] mesh, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  AspectRatio -> Automatic], 
 DensityPlot[vsol[10, x, y], {x, y} \[Element] mesh, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  AspectRatio -> Automatic], 
 DensityPlot[wsol[10, x, y], {x, y} \[Element] mesh, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  AspectRatio -> Automatic]}

Answer 2

Thanks to Alex's prompt answer and instruction about using regularization term. The PDEs in the last question is a simplified set of equations for pilot study because I encountered the same problem in the PDEs to be solved. Then I tried the method of regularization to the following:

Needs["NDSolve`FEM`"]
dr = 0.5; rx = 2.5/2; xs = -(dr + rx); ys = 0.0;
ec = RegionUnion[Disk[{-rx, 0}, dr], Rectangle[{-rx, -dr}, {rx, dr}], 
   Disk[{rx, 0}, dr]];
Region[ec]
bmesh = ToBoundaryMesh[ec, "MaxBoundaryCellMeasure" -> .05];
bmesh["Wireframe"]
f = Function[{vertices, area}, 
   area > 0.004 (0.9 - 0.5 Norm[Mean[vertices]])];
(mesh = ToElementMesh[bmesh, MeshRefinementFunction -> f])["Wireframe"]
u0 = 11.2; v0 = 1; d0 = 0.97; V2 = -0.2;
Eqs = Inactivate[{D[u[t, x, y], t] - 
     d0 (Laplacian[u[t, x, y], {x, y}] + 
        Div[(u[t, x, y]) Grad[w[t, x, y], {x, y}], {x, y}]), 
    D[v[t, x, y], t] - Laplacian[v[t, x, y], {x, y}] - 
     Div[(v[t, x, y]) Grad[w[t, x, y], {x, y}], {x, y}], 
    0.001 D[w[t, x, y], t] - Laplacian[w[t, x, y], {x, y}] - 
     4 \[Pi] (v[t, x, y] - u[t, x, y])}, Laplacian | D];
ics = {u[0, x, y] == u0, v[0, x, y] == v0, w[0, x, y] == 0.0};
bcs = DirichletCondition[w[t, x, y] == V2, True];
{usol, vsol, wsol} = 
  NDSolveValue[{Activate[Eqs] == {NeumannValue[0, True], 
      NeumannValue[0, True], 0}, ics, bcs}, {u, v, 
    w}, {x, y} \[Element] mesh, {t, 0, 10}];

You may find out the differences are the divergence of u or v times the gradient of w following the Laplacian term in the 1st and 2nd PDEs. Small regularization term is applied to the 3rd PDE as the prior setting. Soon I encounter the error as:

NDSolveValue::ndsz: At t == 0.009774811738543655`, step size is effectively zero; singularity or stiff system suspected.

I serached stackexchange to see if it could be solved and modified NDSolve as

  {usol, vsol, wsol} = 
  NDSolveValue[{Activate[Eqs] == {NeumannValue[0, True], 
      NeumannValue[0, True], 0}, ics, bcs}, {u, v, 
    w}, {x, y} \[Element] mesh, {t, 0, 1}, 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, w -> 1}, 
     "MeshOptions" -> {"MaxCellMeasure" -> 0.0005}}];

However, it ran over a day and not finished yet. I guess it's not the best way to remove sigularity or stiffness here. I'd like to learn form you or anone if any approach to solve this type of problem? Thanks for any help and instruction.