Neuman boundary conditions for system of PDE in 2D

Question

The system of equations for three variables: n(x,y,t), b(x,y,t) and a(x,y,t)

The rest of letters are coefficients. System is defined over a square region with the lower left corner at the origin.

Boundary conditions:

Bold n is the normal vector to the boundary (Neuman condition).

And initial conditions:

My code for solution is the following

(*region*)
\[CapitalOmega] = Rectangle[{0, 0}, {1, 1}];

(*constants*)
\[Mu]n = 0.001;
\[Mu]b = 0.001;
\[Chi] = 0.1;
bHat = 1.5;
\[Lambda]0 = 100;
\[Lambda]1 = 100;
\[Lambda]2 = 10;
\[Lambda]4 = 100;
\[Lambda]5 = 10;
nHat = 1;
\[Lambda]7 = 10;
\[Delta] = 0.01;
\[Alpha] = 2.5;

aNeuman = 
  NeumannValue[\[Lambda]7*bHat*a[t, x, y], {x, y} \[Element]RegionBoundary[\[CapitalOmega]]];

(*equations*)
eq = {-D[n[t, x, y], t] + \[Mu]n*Laplacian[n[t, x, y], {x, y}] - \[Chi] *
      Div[n[t, x, y]*Grad[a[t, x, y], {x, y}], {x, y}] + \[Lambda]1*
      a[t, x, y]*b[t, x, y] - \[Lambda]2*n[t, x, y] - \[Lambda]0*(n[t, x, y])^2 == 0,

     -D[a[t, x, y], t] + Laplacian[a[t, x, y], {x, y}] + 
      \[Lambda]4 /2*(1 + Tanh[(1 - b[t, x, y])/\[Delta]]) - (\[Lambda]4 + b[t, x, y])*a[t, x, y] == aNeuman, 

     -D[b[t, x, y], t] + \[Mu]b*Div[n[t, x, y]*Grad[b[t, x, y], {x, y}], {x, y}] - 
     Norm[\[Mu]n*Grad[n[t, x, y], {x, y}] - \[Chi] *n[t, x, y]*Grad[a[t, x, y], {x, y}]] == 0};

(*initial conditions*)
incs = {n[0, x, y] == If[x == 0 || y == 0 || x == 1 || y == 1, nHat, 0],
        a[0, x, y] == 0,
        b[0, x, y] == If[x == 0 || y == 0 || x == 1 || y == 1, bHat, 0]};

(*boundary conditions*)
bcs = {DirichletCondition[n[t, x, y] == nHat*E^(-\[Alpha]*t), True],
       DirichletCondition[b[t, x, y] == bHat, True]};

(*solution*)
{nfun, afun, bfun} = NDSolveValue[Join[eq, bcs, incs], {n, a, b}, {t, 0, 2}, {x, y} \[Element] \[CapitalOmega]];

This code gives the error (apparently because of the wrong way of adding Neumann condition)

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.

I actually have two questions.

1. How to properly implement the boundary condition for a(x,y,t)?

2. When I replace this Neuman condition with the Dirichlet one (which is wrong, but for testing puprposes), the solver starts to work, but then I get another error:

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

This error goes away when I set all coeffitienst of the equations equal to 1. But could you, please, help me with solving the system with original coeffitients?

P.S. My Mathematica version is 12.1

P.P.S. Sorry for that Greek letters. I tried make them consistent with the original equations. If they really distract you, please, tell, I will rename.

Answer 1

We can solve this system with using FEM. But we need to make some modification concerning Norm. This code has some message for the last point t=2

 (*region*)\[CapitalOmega] = Rectangle[{0, 0}, {1, 1}];

(*constants*)
\[Mu]n = 0.001;
\[Mu]b = 0.001;
\[Chi] = 0.1;
bHat = 1.5;
\[Lambda]0 = 100;
\[Lambda]1 = 100;
\[Lambda]2 = 10;
\[Lambda]4 = 100;
\[Lambda]5 = 10;
nHat = 1;
\[Lambda]7 = 10;
\[Delta] = 0.01;
\[Alpha] = 2.5;
eps = 10^-5;
aNeuman = 
  NeumannValue[-\[Lambda]7*bHat*a[t, x, y], {x, y} \[Element] 
    RegionBoundary[\[CapitalOmega]]];

(*equations*)
eq = {-D[n[t, x, y], t] + \[Mu]n*
      Laplacian[n[t, x, y], {x, y}] - \[Chi]*
      Div[n[t, x, y]*Grad[a[t, x, y], {x, y}], {x, y}] + \[Lambda]1*
      a[t, x, y]*b[t, x, y] - \[Lambda]2*
      n[t, x, y] - \[Lambda]0*(n[t, x, y])^2 == 
    0, -D[a[t, x, y], t] + 
     Laplacian[
      a[t, x, y], {x, y}] + \[Lambda]4/
       2*(1 + Tanh[(1 - b[t, x, y])/\[Delta]]) - (\[Lambda]4 + 
        b[t, x, y])*a[t, x, y] == 
    aNeuman, -D[b[t, x, y], t] + \[Mu]b*
      Div[n[t, x, y]*Grad[b[t, x, y], {x, y}], {x, y}] - 
     Sqrt[(\[Mu]n*Grad[n[t, x, y], {x, y}] - \[Chi]*n[t, x, y]*
           Grad[a[t, x, y], {x, y}]) . (\[Mu]n*
           Grad[n[t, x, y], {x, y}] - \[Chi]*n[t, x, y]*
           Grad[a[t, x, y], {x, y}]) + eps^2] == 0};

(*initial conditions*)
incs = {n[0, x, y] == 
    If[x == 0 || y == 0 || x == 1 || y == 1, nHat, 0], 
   a[0, x, y] == 0, 
   b[0, x, y] == If[x == 0 || y == 0 || x == 1 || y == 1, bHat, 0]};

(*boundary conditions*)
bcs = {DirichletCondition[n[t, x, y] == nHat*E^(-\[Alpha]*t), True], 
   DirichletCondition[b[t, x, y] == bHat, True]};


(*solution*)
{nfun, afun, bfun} = 
  NDSolveValue[
    Join[eq, bcs, incs], {n, a, b}, {t, 0, 
     2}, {x, y} \[Element] \[CapitalOmega], 
    Method -> {"FiniteElement", 
      InterpolationOrder -> {n -> 2, a -> 2, b -> 2}, 
      "MeshOptions" -> {"MaxCellMeasure" -> 0.001}}, 
    PrecisionGoal -> 5, AccuracyGoal -> 5] // Quiet; 

Visualization

Table[{Row[{"t =", t}], 
  Plot3D[nfun[t, x, y], {x, y} \[Element] \[CapitalOmega], 
   PlotRange -> All, AxesLabel -> Automatic, PlotLabel -> "n", 
   Mesh -> None, PlotPoints -> 50], 
  Plot3D[afun[t, x, y], {x, y} \[Element] \[CapitalOmega], 
   PlotRange -> All, AxesLabel -> Automatic, PlotLabel -> "a", 
   Mesh -> None, PlotPoints -> 50], 
  Plot3D[bfun[t, x, y], {x, y} \[Element] \[CapitalOmega], 
   PlotRange -> All, AxesLabel -> Automatic, PlotLabel -> "b", 
   Mesh -> None, PlotPoints -> 50]}, {t, .1, 1.9, .3}]