How to solve a reaction-diffusion?

Question

I would like to solve a PDE system reaction-diffusion type (2D spatial + 1 temporal) coupled as described below. Another question of this same system was solved here:

System of nonlinear PDE 2D (Reaction-Diffusion type) with periodic (in laterals) and null flux (top and bottom) boundary condition

Also, initial conditions is assumed with a noise.

The system is shown below:

The code:

L = 5;(*length of square*)
pts = 100;
T = 300;(*Time integration*)
α = 1.5;
k = 0.1;
ρ = 18.5;
a = 92;
b = 64;
d = 10;
γ = 9;
(*system of nonlinear PDE*)
pde = {D[u[t, x, y], 
     t] == (D[u[t, x, y], x, x] + 
       D[u[t, x, y], y, y]) + γ (α - 
        u[t, x, y] - (ρ u[t, x, y] v[t, x, y])/(
        1 + u[t, x, y] + k u[t, x, y] u[t, x, y])), 
   D[v[t, x, y], t] == 
    d (D[v[t, x, y], x, x] + 
        D[v[t, x, y], y, 
         y]) + γ (α (b - v[t, x, y]) - (ρ u [t, x, 
          y] v[t, x, y])/(1 + u[t, x, y] + k u[t, x, y] u[t, x, y]))};
(*Newman boundary condition*)
bc = {(D[u[t, x, y], x] /. x -> 0) == 
    0, (D[u[t, x, y], x] /. x -> 2 L) == 0, 
   u[t, x, 0] == u[t, x, 2 L] , (D[v[t, x, y], x] /. x -> 0) == 
    0, (D[v[t, x, y], x] /. x -> 2 L) == 0, 
   v[t, x, 0] == v[t, x, 2 L]};
(*initial condition*)
ic = {u[0, x, y] == 
    Interpolation[
      Flatten[Table[{x, y, RandomReal[{0.01, 1}]}, {x, 0, 2 L, 
         2/pts}, {y, 0, 2 L, 2/pts}], 1]][x, y], 
   v[0, x, y] == 
    Interpolation[
      Flatten[Table[{x, y, RandomReal[{0.01, 1}]}, {x, 0, 2 L, 
         2/pts}, {y, 0, 2 L, 2/pts}], 1]][x, y]};
eqns = Flatten@{pde, bc, ic};
sol = NDSolve[eqns, {u, v}, {t, 0, T}, {x, 0, 2 L}, {y, 0, 2 L}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> pts, "MaxPoints" -> pts}}];

But, there are errors:

NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent. NDSolve::eerri: Warning: estimated initial error on the specified spatial grid in the direction of independent variable x exceeds prescribed error tolerance. NDSolve::eerri: Warning: estimated initial error on the specified spatial grid in the direction of independent variable y exceeds prescribed error tolerance.

Can anybody help me?

Answer 1

I must admit, I haven't tried to debug what's wrong with your implementation (a first guess could be ensuring periodicity in your initial conditions, see below) - but here's some working FEM code to get you started:

Equations

h[\[Rho]_,k_][u_,v_]=(\[Rho] u v)/(1+u+k u^2);
f[a_,\[Rho]_,k_][u_,v_]=a-u-h[\[Rho],k][u,v];
g[\[Alpha]_,b_,\[Rho]_,k_][u_,v_]=\[Alpha](b-v)-h[\[Rho],k][u,v];
eqns[d_,[Gamma],[Alpha],[Rho],k,a_,b_]={
Derivative[1,0,0][u][t,x,y]-[Gamma] f[a,[Rho],k][u[t,x,y],v[t,x,y]]-Inactive[Div][( Inactive[Grad][u[t,x,y],{x,y}]),{x,y}]==0,
Derivative[1,0,0][v][t,x,y]-[Gamma] g[[Alpha],b,[Rho],k][u[t,x,y],v[t,x,y]]-Inactive[Div][(d Inactive[Grad][v[t,x,y],{x,y}]),{x,y}]==0
};

Initial Conditions

horizontalPeriodicRandomField[{min_,max_},pts_,L_]:=Block[
{rand=RandomReal[{min,max},{pts,pts}]},
rand[[-1]]=rand[[1]];
ListInterpolation[rand,{{-L,L},{-L,L}},PeriodicInterpolation->{True,False}]]
ics={
u[0,x,y]==horizontalPeriodicRandomField[{0.01,1},100,5][x,y],
v[0,x,y]==horizontalPeriodicRandomField[{0.01,1},100,5][x,y]
};

Boundary Conditions

bcs={
PeriodicBoundaryCondition[u[t,x,y],x==-5,TranslationTransform[{10,0}]],
PeriodicBoundaryCondition[v[t,x,y],x==-5,TranslationTransform[{10,0}]]
};

Solution

Monitor[{ufun,vfun}=NDSolveValue[
{eqns[10,9,3/2,37/2,1/10,92,64],bcs,ics},
{u,v},{x,y}\[Element]Rectangle[{-5,-5},{5,5}],{t,0,1},
Method->{"PDEDiscretization"->{"MethodOfLines","SpatialDiscretization"->{"FiniteElement","MeshOptions"->{"MaxCellMeasure"->0.01}}}},
EvaluationMonitor:>(currentTime=Row[{"t = ",CForm[t]}])];,currentTime];

Visualization

This seems to work, but it's quite slow to converge with the parameters and initial conditions you want. FWIW, here's the output of $v$ for 10 evenly-spaced time frames between $t=(0,1]$ (Note the periodic boundary conditions are correctly respected along x):

vFrames=With[{sol=vfun,reg=Rectangle[{-5,-5},{5,5}]},
Table[ContourPlot[sol[t,x,y],{x,y}\[Element]reg,PlotRange->All,Frame->None,Axes->None,PlotPoints->25],{t,Rest@Subdivide[10]}]];
Rasterize[Multicolumn[vFrames,5,Appearance->"Horizontal"],ImageSize->800]

Spectral Methods (?)

Seems like you only have periodic boundary conditions along the horizontal direction, but if this isn't crucial to your problem - you can consider using FFT-based spectral methods. Here's an example walk-through (click on the Wolfram buttons on the top to see the Wolfram code) with code using the Gray-Scott and Cahn-Hilliard reaction-diffusion equations.