How to parallelize this code?

Question

Initialize a matrix a with value of $\alpha$ in central cell as "ice", and $\beta$ in every other cell.

L = 599; alpha = 0.501; beta = 0.3; gamma = 0.0001; 
a = ConstantArray[beta, {L, L}];
a[[Ceiling[L/2], Ceiling[L/2]]] = alpha;

Store all ice cell in a1, and store all water cell in a2.

a1 = ParallelMap[If[# >= alpha, # , 0] &, a, {2}];
a2 = ParallelMap[If[# < alpha, #, 0] &, a, {2}];

Now scan all the cells in a, copy all the ice cells and their neighbors into a1, put zeros in a2 in these positions. Then add $\gamma$(growth rate) to non-zero cell in a1.

Do[If[Total[Boole[{a[[i, j]] >= alpha, a[[i, j + 1]] >= alpha, 
  a[[i, j - 1]] >= alpha, a[[i + 1, j]] >= alpha, 
  a[[i + 1, j + 1]] >= alpha, a[[i + 1, j - 1]] >= alpha, 
  a[[i - 1, j]] >= alpha, a[[i - 1, j + 1]] >= alpha, 
  a[[i - 1, j - 1]] >= alpha}]] >= 1, 
a1[[i, j]] = a[[i, j]] + gamma; a2[[i, j]] = 0;, 
a2[[i, j]] = a[[i, j]]; a1[[i, j]] = 0;], {i, 2, L - 1}, {j, 2, L - 1}];

And I got(For L=9)

a1 // MatrixForm
0   0   0   0   0   0   0   0   0
0   0   0   0   0   0   0   0   0
0   0   0   0   0   0   0   0   0
0   0   0   0.3001  0.3001  0.3001  0   0   0
0   0   0   0.3001  0.5011  0.3001  0   0   0
0   0   0   0.3001  0.3001  0.3001  0   0   0
0   0   0   0   0   0   0   0   0
0   0   0   0   0   0   0   0   0
0   0   0   0   0   0   0   0   0
a2 // MatrixForm
0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
0.3 0.3 0.3 0   0   0   0.3 0.3 0.3
0.3 0.3 0.3 0   0   0   0.3 0.3 0.3
0.3 0.3 0.3 0   0   0   0.3 0.3 0.3
0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3

The result is correct, but it's very slow when L is large.I tried to use ParallelDo instead of Do, but there is side effect.
I was wondering if it's possible to parallelize the above code?

Answer 1

Seems that you're just trying to make the code faster, then try this:

L = 599; alpha = 0.501; beta = 0.3; gamma = 0.0001;

a = ConstantArray[beta, {L, L}];
a[[Ceiling[L/2], Ceiling[L/2]]] = alpha;
{a1, a2} = {a #, a (1 - #)} &@ UnitStep[a - alpha];

{Function[{x, y}, a1[[x - 1 ;; x + 1, y - 1 ;; y + 1]] = 
                  a[[x - 1 ;; x + 1, y - 1 ;; y + 1]]] @@@#, 
 Function[{x, y}, a2[[x - 1 ;; x + 1, y - 1 ;; y + 1]] = 0.] @@@#}&@Position[a, alpha];
a1 = (1 - UnitStep[-#]) gamma + # &@a1;

Answer 2

Like xzczd I suggest you create the initial arrays using UnitStepto define a binary mask:

L = 599; alpha = 0.501; beta = 0.3; gamma = 0.0001;
a = ConstantArray[beta, {L, L}];
a[[Ceiling[L/2], Ceiling[L/2]]] = alpha;

With[{mask = UnitStep[a - alpha]},
  a1 = mask a;
  a2 = (1 - mask) a];

A similar approach can be used to replace the Do loop. The loop goes through each element of a and works out if that element or any of its neighbours is greater than alpha. This is equivalent to asking if the maximum value of the element and its neighbours is greater than alpha. You can use MaxFilterfor that, and then UnitStep again to create a binary mask:

With[{mask = UnitStep[MaxFilter[a, 1] - alpha]},
  a1 = mask (a + gamma);
  a2 = (1 - mask) a];

This is much faster than the procedural Do loop.