# Optimizing code that numerically evaluates line integrals on Mesh/Grid

Posted 1 month ago
311 Views
|
4 Replies
|
6 Total Likes
|
 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] 
4 Replies
Sort By:
Posted 1 month ago
 Tharindu,With just a quick glance at the code, I see one big place to save time. You evaluate the line integrals but do not simplify them. I got more than an order of magnitude speedup by simplifying the complicated expressions. uPAkxPre[kx_, ky_, M_, phi_, t_, t2_] = I uPConj[kx, ky, M, phi, t, t2].D[uP[kx, ky, M, phi, t, t2], kx] // Simplify; uPAkyPre[kx_, ky_, M_, phi_, t_, t2_] = I uPConj[kx, ky, M, phi, t, t2].D[uP[kx, ky, M, phi, t, t2], ky] // Simplify; Doing this one change and adding an accuracy goal on the integrals (to avoid errors when the integrals are zero) int12 = NIntegrate[uPAkx[kx, y1], {kx, x1, x2}, AccuracyGoal -> 5]; I get 41 seconds of runtime (on a very fast MacBook Pro EDIT for meshdims = 100). The next step might be to try different integration options to see if you can optimize the integrals. I hope this helps.BTW, I appreciated the nicely commented code -- It makes it much easier to help.Regards,Neil
Posted 1 month ago
 Thank you for your insight, Neil - your suggestions did speed my code up significantly! I will still keep trying different integration options.Best,Tharindu
 Thirindu,I woke up with an additional thought that helps. The problem that I see is your line integral is complicated so it takes time to numerically evaluate and integrate. I tried a FullSimplify on it and I gave up after 20 min. The code below, however, yields a MUCH simpler expression that then drops the time for meshdim =100 from 41 seconds to 15 seconds! The idea is to simplify the intermediate expressions so the large expression does not "blow up" and then FullSimplify can run in a reasonable time. uPAkxPre[kx_, ky_, M_, phi_, t_, t2_] = (I uPConj[kx, ky, M, phi, t, t2] // Simplify).(D[uP[kx, ky, M, phi, t, t2], kx] // Simplify) // FullSimplify; uPAkyPre[kx_, ky_, M_, phi_, t_, t2_] = (I uPConj[kx, ky, M, phi, t, t2] // Simplify).(D[uP[kx, ky, M, phi, t, t2], ky] // Simplify) // FullSimplify; Using Simplify instead of FullSimplify takes less time on this step but runs in 18 seconds (vs the 15). However, since the FullSimplification runs only once, its not so bad (maybe 5 min). Either way you are in much better shape. How does this compare to your MATLAB computations? If you are far behind, I will think a bit more about it.Regards,Neil