NDSolve cannot consistently handle a simple 2D reaction-diffusion PDE with spiky (and random) initial conditions

Question

I am still fighting with a modified version of the 2D reaction-diffusion PDE that I have already posted in a previous question ( Reaction-diffusion PDE with NDSolve: either very slow or very inaccurate ). Now I use a slightly different set of initial conditions, and I cannot convince NDSolve to do the task for me properly and within acceptable time. This is the system I use now:

    SetDirectory[NotebookDirectory[]];


    Needs["DifferentialEquations`NDSolveProblems`"]
    Needs["DifferentialEquations`NDSolveUtilities`"]
    Needs["FunctionApproximations`"]
    a = AbsoluteTime[];

    rs = 0.1;
    rn = 0.1;
    Diff = 1;
    Difn = 0.001;

    sigmas0 = 2;
    sigmans0 = 1;

    s0 = 1000.;           (* Initial biomass of Starters *)
    cs =   10;            (* Number of starter colonies *)
    ns0 = .2;             (* Initial biomass of Non-starters *)
    cns = 10;             (* Number of non-starter colonies *)

    time = 500;
    size = 50;
    rep = 1;

    k = 2 cs + 1;
    l = 2 cs + cns;


    rnd = RandomReal[2 size, 2 cs + 2 cns];
    rnd = rnd - size;

    soln = NDSolve[{
        D[S[t, x, y], t] == -rs S[t, x, y],
        D[F[t, x, y], t] == 
         Diff (D[F[t, x, y], x, x] + D[F[t, x, y], y, y]) + 
          rs S[t, x, y] - rn F[t, x, y] NS[t, x, y],
        D[NS[t, x, y], t] == 
         Difn (D[NS[t, x, y], x, x] + D[NS[t, x, y], y, y]) + 
          rn F[t, x, y] NS[t, x, y],

        S[t, -size, y] == S[t, size, y], S[t, x, -size] == S[t, x, size],
        F[t, -size, y] == F[t, size, y], F[t, x, -size] == F[t, x, size],
        NS[t, -size, y] == NS[t, size, y], 
        NS[t, x, -size] == NS[t, x, size],

        S[0, x, y] == 
         Sum[s0/cs Exp[-((x - rnd[[n]])^2/(2 sigmas0^2) + (y - 
          rnd[[cs + n]])^2/(2 sigmas0^2))], {n, 1, cs}],
        F[0, x, y] == 0.,
        NS[0, x, y] == 
         Sum[ns0/cns Exp[-((x - rnd[[nn]])^2/(2 sigmans0^2) + (y - 
                  rnd[[cns + nn]])^2/(2 sigmans0^2))], {nn, k, l}]},

       {S, F, NS}, {t, 0, time}, {x, -size, size}, {y, -size, size}, 
       Method -> {"StiffnessSwitching", "NonstiffTest" -> False}]; 

Since the culprit is most probably the spiky initial condition both in S[0,x,y] and in NS[0,x,y], I have tried the "StiffnessSwitching" method in this attempt, but even then the kernel stops with the error message

NDSolve::ndsz: At t == 28.71100997366924`, step size is effectively zero; singularity or stiff system suspected. >>

Singularity is out of the question (at least in principle it is), as the initial conditions are the sums of Gaussian functions (for S and NS) and the constant zero function (for F).

Does anyone have an idea how to make this work consistently (i.e., regardless of the actual parameters defining the initial conditions)? I have also tried setting the MaxStepSize option to 0.1 - then I got no error message, but no result either, for more than 10 hours...

Answer 1

NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent.

So there is a problem from the very beginning. What can I recommend? We set

sigmas0 = 1/2;sigmans0 = 1/2;rnd = RandomReal[{-size/3, size/3}, 2 cs + 2 cns];

Accordingly, remove the following line

rnd = rnd - size;

After that, we have nice pictures not for every random distribution, but for many of them

    Table[Plot3D[F[t, x, y] /. soln, {x, -size, size}, {y, -size, size}, 
  PlotRange -> All, ColorFunction -> Hue, Mesh -> None, 
  AxesLabel -> {"x", "y", ""}, PlotLabel -> t], {t, .2*time, 
  time, .2*time}]

    Table[Plot3D[NS[t, x, y] /. soln, {x, -size, size}, {y, -size, size}, 
  PlotRange -> All, ColorFunction -> Hue, Mesh -> None, 
  AxesLabel -> {"x", "y", ""}, PlotLabel -> t], {t, 0, time, .2*time}]