How can I speed up NDSolve

Question

I have written a simulation that solve a pde, I used a method I found on this site to speed up the calculation but it is still slow. Is there any way to speed up the prosses?

L = 25;(*length of square*)
T = 70;(*Time integration*)
pts = 1000;
f = 0.2; (*The infiltration contrast between vegetated and bare-soil \
areas*)
A = 40.0; (*The maximal infiltration rate, obtained in densely \
vegetated areas*)
n = 30.0; (*the soil-water evaporation rate*)
R = 30.0; (*the shading effect*)
E1 = 2.0; (*root to shoot ratio of B1*)
E2 = 0.0;(*root to shoot ratio of B2*)
M1 = 0.75; (*the mortality rate of B1*)
M2 = 7.5; (*the mortality rate of B2*)
K1 = 5.0; (*the maximum standing biomass of B1*)
K2 = 0.5; (*the maximum standing biomass of B2*)
Λ1 = 0.05; (*growth rate of B1*)
Λ2 = 1.0; (*growth rate of B2*)
Γ1 = 10.0; (*the water uptake rate of B1*)
Γ2 = 15.0; (*the water uptake rate of B2*)
q = 5.0; (*reference biomass value beyond which he effect of \
vegetation on infiltration becomes significant*)
Y1 = 1.0; (*relative contribution of the woody specie (B1)*)
Y2 = 100; (*relative contribution of the herbaceous specie (B2)*)
DB1 = 0.1; (*the biomass expansion rate of B1*)
DB2 = 0.05; (*the biomass expansion rate of B2*)
DW = 3.0; (*the soil-water diffusion*)
DH = 1.0; (*the rate of surface-water spread*)
i = A (Y1 B1[t, x, y] + Y2 B2[t, x, y] + q f)/(
   Y1 B1[t, x, y] + Y2 B2[t, x, y] + q );
P = 1.7;

(*system of nonlinear PDE*)
pde = {D[B1[t, x, y], 
     t] == Λ1 W[t, x, y] B1[t, x, 
       y] (1 - B1[t, x, y]/K1) (1 + E1 B1[t, x, y])^2 - 
     M1 B1[t, x, y] + 
     DB1 (D[B1[t, x, y], x, x] + D[B1[t, x, y], y, y]), 
   D[B2[t, x, y], 
     t] == Λ2 W[t, x, y] B2[t, x, 
       y] (1 - B2[t, x, y]/K2) (1 + E2 B2[t, x, y])^2 - 
     M2 B2[t, x, y] + 
     DB2 (D[B2[t, x, y], x, x] + D[B2[t, x, y], y, y]),
   D[W[t, x, y], t] == 
    i H[t, x, y] - (n W[t, x, y])/(
     1 + R B1[t, x, y]/K1 ) - Γ1 W[t, x, y] B1[t, x, 
       y] (1 + E1 B1[t, x, y])^2 - Γ2 W[t, x, y] B2[t, 
       x, y] (1 + E2 B2[t, x, y])^2 + 
     DW (D[W[t, x, y], x, x] + D[W[t, x, y], y, y]),
   D[H[t, x, y], t] == 
    P - i H[t, x, y] + 
     DH (D[H[t, x, y]^2, x, x] + D[H[t, x, y]^2, y, y])};
(*Newman boundary condition*)
bc = {(D[B1[t, x, y], x] /. x -> -L) == 
    0, (D[B1[t, x, y], x] /. x -> L) == 
    0, (D[B1[t, x, y], y] /. y -> -L) == 
    0, (D[B1[t, x, y], y] /. y -> L) == 
    0, (D[B2[t, x, y], x] /. x -> -L) == 
    0, (D[B2[t, x, y], x] /. x -> L) == 
    0, (D[B2[t, x, y], y] /. y -> -L) == 
    0, (D[B2[t, x, y], y] /. y -> L) == 
    0, (D[W[t, x, y], x] /. x -> -L) == 
    0, (D[W[t, x, y], x] /. x -> L) == 
    0, (D[W[t, x, y], y] /. y -> -L) == 
    0, (D[W[t, x, y], y] /. y -> L) == 0,
   (D[H[t, x, y], x] /. x -> -L) == 0, (D[H[t, x, y], x] /. x -> L) ==
     0, (D[H[t, x, y], y] /. y -> -L) == 
    0, (D[H[t, x, y], y] /. y -> L) == 0};

ic = {B1[0, x, y] == If[y < 0, 1, 0],
   B2[0, x, y] == If[Abs[y] < 1, 1, 0],
   W[0, x, y] == If[y < 0, 9, 12],
   H[0, x, y] == If[y < 0, 3, 5]};

eqns = Flatten@{pde, bc, ic};

sol = NDSolve[eqns, {B1, B2, W, H}, {t, 0, T}, {x, -L, L}, {y, -L, L},
    Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> pts, "MaxPoints" -> pts}, 
     Method -> {"IDA", "ImplicitSolver" -> {"GMRES"}}}];

maybe using multiple CPU?

Thanks in advance.

(the ecology tag related because this simulation is for a model that calculate the time and space growth of plants in desert like area.

Answer 1

There are no gradients in x and y at L> 10. Nothing happens when T> 10. The number pts can be taken 81. Therefore, we can take L=10, T=10, pts=81. The initial data can be expressed using UnitStep[]. Then we get

L = 10;(*length of square*)T = 10;(*Time integration*)pts = 81;
f = 0.2;(*The infiltration contrast between vegetated and bare-soil \
areas*)A = 40.0;(*The maximal infiltration rate,obtained in densely \
vegetated areas*)n = 30.0;(*the soil-water evaporation rate*)R = \
30.0;(*the shading effect*)E1 = 2.0;(*root to shoot ratio of B1*)E2 = \
0.0;(*root to shoot ratio of B2*)M1 = 0.75;(*the mortality rate of \
B1*)M2 = 7.5;(*the mortality rate of B2*)K1 = 5.0;(*the maximum \
standing biomass of B1*)K2 = 0.5;(*the maximum standing biomass of \
B2*)Λ1 = 0.05;(*growth rate of B1*)Λ2 = \
1.0;(*growth rate of B2*)Γ1 = 10.0;(*the water uptake \
rate of B1*)Γ2 = 15.0;(*the water uptake rate of B2*)q \
= 5.0;(*reference biomass value beyond which he effect of vegetation \
on infiltration becomes significant*)Y1 = 1.0;(*relative contribution \
of the woody specie (B1)*)Y2 = 100;(*relative contribution of the \
herbaceous specie (B2)*)DB1 = 0.1;(*the biomass expansion rate of \
B1*)DB2 = 0.05;(*the biomass expansion rate of B2*)DW = 3.0;(*the \
soil-water diffusion*)DH = 1.0;(*the rate of surface-water spread*)i =
  A (Y1 B1[t, x, y] + Y2 B2[t, x, y] + q f)/(Y1 B1[t, x, y] + 
     Y2 B2[t, x, y] + q);
P = 1.7;
us[y_] := 2 (UnitStep[y] - 1/2)
(*system of nonlinear PDE*)
pde = {D[B1[t, x, y], 
     t] == Λ1 W[t, x, y] B1[t, x, 
       y] (1 - B1[t, x, y]/K1) (1 + E1 B1[t, x, y])^2 - 
     M1 B1[t, x, y] + 
     DB1 (D[B1[t, x, y], x, x] + D[B1[t, x, y], y, y]), 
   D[B2[t, x, y], 
     t] == Λ2 W[t, x, y] B2[t, x, 
       y] (1 - B2[t, x, y]/K2) (1 + E2 B2[t, x, y])^2 - 
     M2 B2[t, x, y] + 
     DB2 (D[B2[t, x, y], x, x] + D[B2[t, x, y], y, y]), 
   D[W[t, x, y], t] == 
    i H[t, x, 
       y] - (n W[t, x, y])/(1 + 
        R B1[t, x, y]/K1) - Γ1 W[t, x, y] B1[t, x, 
       y] (1 + E1 B1[t, x, y])^2 - Γ2 W[t, x, y] B2[t, 
       x, y] (1 + E2 B2[t, x, y])^2 + 
     DW (D[W[t, x, y], x, x] + D[W[t, x, y], y, y]), 
   D[H[t, x, y], t] == 
    P - i H[t, x, y] + 
     DH (D[H[t, x, y]^2, x, x] + D[H[t, x, y]^2, y, y])};
(*Newman boundary condition*)
bc = {(D[B1[t, x, y], x] /. x -> -L) == 
    0, (D[B1[t, x, y], x] /. x -> L) == 
    0, (D[B1[t, x, y], y] /. y -> -L) == 
    0, (D[B1[t, x, y], y] /. y -> L) == 
    0, (D[B2[t, x, y], x] /. x -> -L) == 
    0, (D[B2[t, x, y], x] /. x -> L) == 
    0, (D[B2[t, x, y], y] /. y -> -L) == 
    0, (D[B2[t, x, y], y] /. y -> L) == 
    0, (D[W[t, x, y], x] /. x -> -L) == 
    0, (D[W[t, x, y], x] /. x -> L) == 
    0, (D[W[t, x, y], y] /. y -> -L) == 
    0, (D[W[t, x, y], y] /. y -> L) == 
    0, (D[H[t, x, y], x] /. x -> -L) == 
    0, (D[H[t, x, y], x] /. x -> L) == 
    0, (D[H[t, x, y], y] /. y -> -L) == 
    0, (D[H[t, x, y], y] /. y -> L) == 0};

ic = {B1[0, x, y] == -(us[y] - 1)/2, 
   B2[0, x, y] == (-us[y - 1] + 1) (us[y + 1] + 1)/4, 
   W[0, x, y] == 9 + (12 - 9) (us[y] + 1)/2, 
   H[0, x, y] == 3 + (5 - 3) (us[y] + 1)/2};

eqns = Flatten@{pde, bc, ic};
sol = NDSolve[eqns, {B1, B2, W, H}, {t, 0, T}, {x, -L, L}, {y, -L, L},
     Method -> {"IndexReduction" -> Automatic, 
      "EquationSimplification" -> "Residual", 
      "PDEDiscretization" -> {"MethodOfLines", 
        "SpatialDiscretization" -> {"TensorProductGrid", 
          "MinPoints" -> pts, "MaxPoints" -> pts, 
          "DifferenceOrder" -> 2}}}] // Quiet;

{DensityPlot[B1[T, x, y] /. sol, {x, -L, L}, {y, -L, L}, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic], 
 DensityPlot[B2[T, x, y] /. sol, {x, -L, L}, {y, -L, L}, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  PlotRange -> All], 
 DensityPlot[W[T, x, y] /. sol, {x, -L, L}, {y, -L, L}, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic], 
 DensityPlot[H[T, x, y] /. sol, {x, -L, L}, {y, -L, L}, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic]}