Discretize nonlinear PDE-system fails

Question

Inspired by this post Solving the system of nonlinear PDE

I would like to discretize these two stationary nonlinear pde's for further use with FEMsolver of NDSolve.

parameters used:

{beta, chiu, chiv, kappa, Du, Dv} = {9, 0.3, 2.4, 0.001, 0.733, 0.733}
{a, b, c, d} = {-0.1, 0.4, 0, 0.2}
region = Rectangle[{-1, -1}, {1, 1}];
u0 = RandomVariate[NormalDistribution[22.222, 0.0001]];
v0 = RandomVariate[NormalDistribution[44.444, 0.0001]];

discretization (my unsuccessfull apporach...)

vd = NDSolve`VariableData[{"DependentVariables","Space" } -> {{u, v}, {x, y} }] ;
sd = NDSolve`SolutionData[{"Space"} -> {Discretize[region]} ] 

c11 = Du - beta chiu (a - b) (u[x, y] v[x, y])/(u[x, y] + v[x, y])^2;
c12 = beta chiu (a - b) (u[x, y] u[x, y])/(u[x, y] + v[x, y])^2;
c21 = -beta chiv (c - d) (v[x, y] v[x, y])/(u[x, y] + v[x, y])^2;
c22 = Dv + beta chiv (c - d) (u[x, y] v[x, y])/(u[x, y] + v[x, y])^2;
cdata = InitializePDECoefficients[vd, sd, 
"DiffusionCoefficients" -> { { -{{c11, 0}, {0, c11}}, -{{c12,0}, {0, c12}}}, {-{{c21 , 0}, {0, c21 }}, -{{c22, 0}, {0, c22}}} }, 
"ReactionCoefficients" -> {{ (a u[ x, y] + b v[ x, y])/(u[ x, y] + v[ x, y]) - 
  kappa*(u[ x, y] + v[ x, y]), 0}, {0, (c u[ x, y] + d v[ x, y])/(u[ x, y] + v[ x, y]) - 
  kappa*(u[ x, y] + v[ x, y]), 0} } ]

bcdata = InitializeBoundaryConditions[vd,sd, { {DirichletCondition[u[ x, y] == u0, True],DirichletCondition[v[ x, y] == v0, True]}  }]

mdata = InitializePDEMethodData[vd, sd] 

dpde = DiscretizePDE[cdata, mdata, sd, "SaveFiniteElements" -> True,"AssembleSystemMatrices" -> True] 

Unfortunately the last step

{load, stiffness, damping, mass} = dpde["All"]

fails and gives the error "Set::shape: Lists {load,stiffness,damping,mass} and <<1>> are not the same shape."

Any ideas what went wrong here and how to correct it to make the nonlinear FEMmodel run? Thanks!

Answer 1

Here is the system:

{beta, chiu, chiv, kappa, Du, Dv} = {9, 0.3, 2.4, 0.001, 0.733, 0.733};
{a, b, c, d} = {-0.1, 0.4, 0, 0.2};
region = Rectangle[{-1, -1}, {1, 1}];
u0 = RandomVariate[NormalDistribution[22.222, 0.0001]];
v0 = RandomVariate[NormalDistribution[44.444, 0.0001]];

c11 = Du - 
   beta  chiu  (a - b)  (u[x, y]  v[x, y])/(u[x, y] + v[x, y])^2;
c12 = beta  chiu  (a - b)  (u[x, y]  u[x, y])/(u[x, y] + v[x, y])^2;
c21 = -beta  chiv  (c - d)  (v[x, y]  v[x, y])/(u[x, y] + v[x, y])^2;
c22 = Dv + 
   beta  chiv  (c - d)  (u[x, y]  v[x, y])/(u[x, y] + v[x, y])^2;

vars = {{u[x, y], v[x, y]}, {x, y}};
op = DiffusionPDETerm[vars, {{c11, c12}, {c21, c22}}] + 
   ReactionPDETerm[
    vars, {{(a  u[x, y] + b  v[x, y])/(u[x, y] + 
          v[x, y]) - kappa*(u[x, y] + v[x, y]), 
      0}, {0, (c  u[x, y] + d  v[x, y])/(u[x, y] + 
          v[x, y]) - kappa*(u[x, y] + v[x, y])}}];

NDSolveValue[{op == {0, 0}, DirichletCondition[u[x, y] == u0, True], 
  DirichletCondition[v[x, y] == v0, True]}, {u, v}, {x, y} \[Element] 
  region, "InitialSeeding" -> {u[x, y] == 1, v[x, y] == 1}]

but you need a decent initial seed, otherwise it will fail.

FindRoot::dfmin: The minimal damping factor of 1/10000 has been reached.

And here is the low-level approach:

Needs["NDSolve`FEM`"]

vd = NDSolve`VariableData[{"DependentVariables", 
     "Space"} -> {{u, v}, {x, y}}];
sd = NDSolve`SolutionData[{"Space"} -> {ToElementMesh[region]}]

cdata = InitializePDECoefficients[vd, sd, 
   "DiffusionCoefficients" -> {{-{{c11, 0}, {0, c11}}, -{{c12, 0}, {0,
          c12}}}, {-{{c21, 0}, {0, c21}}, -{{c22, 0}, {0, c22}}}}, 
   "ReactionCoefficients" -> {{(a  u[x, y] + b  v[x, y])/(u[x, y] + 
          v[x, y]) - kappa*(u[x, y] + v[x, y]), 
      0}, {0, (c  u[x, y] + d  v[x, y])/(u[x, y] + v[x, y]) - 
       kappa*(u[x, y] + v[x, y])}}, 
   "VerificationData" -> {"DependentVariables" -> {1, 1}}];

bcdata = 
 InitializeBoundaryConditions[vd, 
  sd, {{DirichletCondition[u[x, y] == u0, True]}, {DirichletCondition[
     v[x, y] == v0, True]}}]
mdata = InitializePDEMethodData[vd, sd]

dpde = DiscretizePDE[cdata, mdata, sd]