4
$\begingroup$

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.

$\endgroup$
4
  • 2
    $\begingroup$ If no option is specified the integration is faster. In other words why do you specify pts=1000? Also, this does run on multiple CPUs. $\endgroup$
    – user21
    Feb 20, 2020 at 12:40
  • $\begingroup$ no particular reason, just wanted high definition. and it does not run on multiple CPU, how can i make it so it does? $\endgroup$
    – Ben9172
    Feb 20, 2020 at 12:46
  • $\begingroup$ Works in parallel on my machine. Have you changed the setup on yours? $\endgroup$
    – user21
    Feb 20, 2020 at 15:17
  • 1
    $\begingroup$ As @user21 suggests, requiring 1000 points is a big problem. Since you have a 2D system, I think that means that NDSolve discretizes the PDE into 1000*1000 ODEs. Not surprising that's going to be slow. $\endgroup$
    – Chris K
    Feb 20, 2020 at 16:26

1 Answer 1

7
$\begingroup$

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]}

Figure 1

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.