I am trying to calculate an algebraically complicated area integral over a 2D area. Motivated by Stokes' theorem, I break the area into tiny squares using CoordinateBoundingBoxArray[], and evaluate appropriate versions of the algebraic quantity's line integral version across each edge. I already made this calculation work efficiently using MATLAB using arrays, but recently tried the same on Mathematica because its capability to handle algebraic expressions is extremely convenient.

However, my Mathematica implementation is very slow (takes 15-30 minutes on a supercomputer's 120 GB core), whereas my MATLAB code takes under a minute on a 16 GB desktop (also single core). I did my best to look at optimization tricks and hints others have suggested online, and was able to bring computing time down from an hour. But I feel like I can do better, given how Wolfram seems like such an advanced language.

Therefore, does anyone have any hints to offer for my Mathematica implementation below? I am already substituting constant values in advance, and using ?NumericQ. Processing the algebra (derivatives, etc) takes less than a minute; i.e. the time-consuming part is evaluating function magicP[] for variable berrysingsP. I apologize for the messy algebra - but this is a minimally functioning example, where all defined quantities are required for a run. The goal is to be able to use a dense-enough mesh (meshdims > 100) to identify singular points on the grid using ListPlot3D[] at the end.

Thanks in advance.

(*** CONSTANTS ***)
kxOFFSET = 0; kyOFFSET = 0;Dirac1x = 4 Pi/(3 Sqrt[3]); Dirac2x = 8 Pi/(3 Sqrt[3]);BZx1 = Pi/Sqrt[3]; BZx2 = 3 Pi/Sqrt[3]; BZy1 = -2 Pi/3; BZy2 = 2 Pi/3;MM = 0.01; phii = Pi/2; tt = 1; t22 = 1;
meshdims = 5;(* 100 takes about 950 secs , ie ~ 15 mins *)
maxprecision = 15;

(*** Discretize mesh: I am using 'center' to try to make sure all iterations of magicP[] have valid edge points to evaluate over. ***)
bounds = Round[N[{{BZx1 - kxOFFSET, BZx2 - kxOFFSET}, {BZy1 - kyOFFSET, BZy2 - kyOFFSET}}]]
mesh = Block[{$MaxExtraPrecision = maxprecision}, N[CoordinateBoundingBoxArray[Transpose[bounds], Into[meshdims]]]];
center = Block[{$MaxExtraPrecision = maxprecision}, N[CoordinateBoundingBoxArray[Transpose[bounds], Into[meshdims], Center]]];
i1 = Dimensions[center][[1]]; (* to iterate over *)
i2 = Dimensions[center][[2]] ;

(**** required functions here ****)
b1 = {{-Sqrt[3]/2}, {3/2}}; b2 = {{-Sqrt[3]/2}, {-3/2}}; b3 = {{Sqrt[3]}, {0}};
h0[kx_, ky_, M_, phi_, t_, t2_] = {0};
hx[kx_, ky_, M_, phi_, t_, t2_] = t (1 + Cos[{kx, ky}.b1] + Cos[{kx, ky}.b2]);
hy[kx_, ky_, M_, phi_, t_, t2_] = t (Sin[{kx, ky}.b1] - Sin[{kx, ky}.b2]);
hz[kx_, ky_, M_, phi_, t_, t2_] = M - 2 t2 Sin[phi] (Sin[{kx, ky}.b1] + Sin[{kx, ky}.b2] + Sin[{kx, ky}.b3]);
H[kx_, ky_, M_, phi_, t_, t2_] = {{h0[kx, ky, M, phi, t, t2][[1]] + hz[kx, ky, M, phi, t, t2][[1]], hx[kx, ky, M, phi, t, t2][[1]] - I hy[kx, ky, M, phi, t, t2][[1]]}, {hx[kx, ky, M, phi, t, t2][[1]] + I hy[kx, ky, M, phi, t, t2][[1]], h0[kx, ky, M, phi, t, t2][[1]] - hz[kx, ky, M, phi, t, t2][[1]]}};
eP[kx_, ky_, M_, phi_, t_, t2_] = h0[kx, ky, M, phi, t, t2] + Sqrt[(hx[kx, ky, M, phi, t, t2]^2) + (hy[kx, ky, M, phi, t, t2]^2) + (hz[kx, ky, M, phi, t, t2]^2)];
eM[kx_, ky_, M_, phi_, t_, t2_] = h0[kx, ky, M, phi, t, t2] - Sqrt[(hx[kx, ky, M, phi, t, t2]^2) + (hy[kx, ky, M, phi, t, t2]^2) + (hz[kx, ky, M, phi, t, t2]^2)];
uP[kx_, ky_, M_, phi_, t_, t2_] :=Flatten[(((eP[kx, ky, M, phi, t, t2][[1]] + hz[kx, ky, M, phi, t, t2][[1]])/(2 eP[kx, ky, M, phi, t, t2][[1]]))^(1/2)) {{(hx[kx, ky, M, phi, t, t2][[1]] - I hy[kx, ky, M, phi, t, t2][[1]])/(eP[kx, ky, M, phi, t, t2][[1]] + hz[kx, ky, M, phi, t, t2][[1]])}, {-1}}, 1];
uPConj[kx_, ky_, M_, phi_, t_, t2_] := Flatten[(((eP[kx, ky, M, phi, t, t2][[1]] + hz[kx, ky, M, phi, t, t2][[1]])/(2 eP[kx, ky, M, phi, t, t2][[1]]))^(1/2)) {{(hx[kx, ky, M, phi, t, t2][[1]] + I hy[kx, ky, M, phi, t, t2][[1]])/(eP[kx, ky, M, phi, t, t2][[1]] + hz[kx, ky, M, phi, t, t2][[1]])}, {-1}}, 1];

(**** line integrals: ****)
uPAkxPre[kx_, ky_, M_, phi_, t_, t2_] = I  uPConj[kx, ky, M, phi, t, t2].D[uP[kx, ky, M, phi, t, t2], kx];
uPAkyPre[kx_, ky_, M_, phi_, t_, t2_] = I  uPConj[kx, ky, M, phi, t, t2].D[uP[kx, ky, M, phi, t, t2], ky];
(* substitute constants now to maximize numerical effectiveness *)
uPAkx[kx_?NumericQ, ky_?NumericQ] := uPAkxPre[kx, ky, MM, phii, tt, t22];
uPAky[kx_?NumericQ, ky_?NumericQ] := uPAkyPre[kx, ky, MM, phii, tt, t22];

(* evaluate above along mesh cell edges *)
magicP[k1_, k2_] := 
  Block[{x1 = k1[[1]], x2 = k2[[1]], y1 = k1[[2]], y2 = k2[[2]], 
    int12, int23, int34, int41, berryphaseP},
   int12 = NIntegrate[uPAkx[kx, y1], {kx, x1, x2}];
   int23 = NIntegrate[uPAky[x2, ky], {ky, y1, y2}];
   int34 = NIntegrate[uPAkx[kx, y2], {kx, x2, x1}];
   int41 = NIntegrate[uPAky[x1, ky], {ky, y2, y1}];
   berryphaseP = int12 + int23 + int34 + int41;
   berryphaseP];
berrysingsP = ParallelTable[magicP[mesh[[p, q]], mesh[[p + 1, q + 1]]], {p, 1, i1}, {q, 1, i2}]; // AbsoluteTiming
(* round so not identically zero *)
roundberrysingsP = Parallelize[Round[berrysingsP, 0.01]];

(* now bring back to space coordinates and visualize *)
singsP = Parallelize[MapThread[Append, {Flatten[center, 1], Flatten[roundberrysingsP, 1]}]];
ListPlot3D[singsP]