How to force evaluation of an expression with Part left unnevaluated?

Question

My goal is to minimize a real valued multi-variable function averageFreeEnergyDensity. The function accepts an integer argument n, a set of real parameters j,d,a,h,dx and two real matrices th,ph. I want to find matrices th,ph which minimize the function. FindMinimum minimizes the function well enough but it is too slow at the moment. I am trying to compile the function to accelerate FindMinimum. If I compile the function for arbitrary n it is at least an order of magnitude more slow than if I make a set of functions each of which is defined for a particular n. My current fastest solution is

averageFreeEnergyDensityCompiled["square",n_]:=
 averageFreeEnergyDensityCompiled["square", n]= 
  Compile[{{j, _Real}, {d, _Real}, {a, _Real}, {h, _Real}, {dx,
_Real}, {th, _Real, 2}, {ph, _Real, 2}}, Evaluate[
    (*apply boundary conditions*)
    s[i1_, i2_] := Which[
      i1 == n + 1 && i2 == n + 1, s[1, 1],
      i1 == n + 1, s[1, i2],
      i2 == n + 1, s[i1, 1],
      True, {Cos[ph[[i1, i2]]] Sin[th[[i1, i2]]], 
       Sin[ph[[i1, i2]]] Sin[th[[i1, i2]]], Cos[th[[i1, i2]]]}
      ];
    (*define finite difference functions*)
    dsdx[i1_, i2_] := (s[i1 + 1, i2] - s[i1, i2])/dx;
    dsdy[i1_, i2_] := (s[i1, i2 + 1] - s[i1, i2])/dx;
    (*return average free energy density function*)
    Sum[j (dsdx[i1, i2].dsdx[i1, i2] + dsdy[i1, i2].dsdy[i1, i2]) + 
       a s[i1, i2][[3]]^2 - h s[i1, i2][[3]], {i1, 1, n}, {i2, 1, n}]/
     n^2
    ]]

I need to delay evaluation because the expression being compiled can only be evaluated for particular n. I store the function because I will evaluate it for many values of j,d,a,h,dx. This solution is fast but it throws an error the first time the function is evaluated

Block[{n = 100}, 
 AbsoluteTiming[
  averageFreeEnergyDensityCompiled["square", n][1, 1, 1, 0, 0.1, 
   Array[0 &, {n, n}], RandomReal[2 Pi, {n, n}]]]]
Part::partd: Part specification ph$[[2,1]] is longer than depth of object. >>
Part::partd: Part specification th$[[2,1]] is longer than depth of object. >>
Part::partd: Part specification ph$[[2,1]] is longer than depth of object. >>
General::stop: Further output of Part::partd will be suppressed during this calculation. >>

After several attempts to evaluate the Part function Mathematica gives up and evaluates the expression ignoring Part. The final expression can then be compiled. I want to manually tell Mathematica to ignore Part. I have tried using Inactive,Activate but this slows down the final compiled function by orders of magnitude. I would ideally give an option to Evaluate like Evaluate["insert expression here",Holding->Part] but Evaluate has no options that I know of.

Note: The boundary condition I have used here is simple and can be taken out of compile, but I will be using more complicated boundary conditions once I solve this problem. I do not need to use finite difference steps of order greater than 1. I will also be using this code on a triangular grid which will further change the finite difference.

Answer 1

I believe Block will work for you here:

mem : averageFreeEnergyDensityCompiled["square", n_] := mem =
  Block[{Part},
   Compile[{{j, _Real}, {d, _Real}, {a, _Real}, {h, _Real}, {dx, _Real}, {th, _Real, 
      2}, {ph, _Real, 2}}, 
    Evaluate[(*apply boundary conditions*)
     s[i1_, i2_] := 
      Which[i1 == n + 1 && i2 == n + 1, s[1, 1], i1 == n + 1, s[1, i2], i2 == n + 1, 
       s[i1, 1], True, {Cos[ph[[i1, i2]]] Sin[th[[i1, i2]]], 
        Sin[ph[[i1, i2]]] Sin[th[[i1, i2]]], Cos[th[[i1, i2]]]}];
     (*define finite difference functions*)
     dsdx[i1_, i2_] := (s[i1 + 1, i2] - s[i1, i2])/dx;
     dsdy[i1_, i2_] := (s[i1, i2 + 1] - s[i1, i2])/dx;
     (*return average free energy density function*)
     Sum[j (dsdx[i1, i2].dsdx[i1, i2] + dsdy[i1, i2].dsdy[i1, i2]) + a s[i1, i2][[3]]^2 - 
        h s[i1, i2][[3]], {i1, 1, n}, {i2, 1, n}]/n^2]]
  ]

However I believe Compile`GetElement will also work and perform better too. Search for "GetElement" (on this site) for more information.

mem : averageFreeEnergyDensityCompiled["square", n_] := mem =
  With[{Part = Compile`GetElement},
   Compile[{{j, _Real}, {d, _Real}, {a, _Real}, {h, _Real}, {dx, _Real}, {th, _Real, 
      2}, {ph, _Real, 2}}, 
    Evaluate[(*apply boundary conditions*)
     s[i1_, i2_] := 
      Which[i1 == n + 1 && i2 == n + 1, s[1, 1], i1 == n + 1, s[1, i2], i2 == n + 1, 
       s[i1, 1], True, {Cos[ph[[i1, i2]]] Sin[th[[i1, i2]]], 
        Sin[ph[[i1, i2]]] Sin[th[[i1, i2]]], Cos[th[[i1, i2]]]}];
     (*define finite difference functions*)
     dsdx[i1_, i2_] := (s[i1 + 1, i2] - s[i1, i2])/dx;
     dsdy[i1_, i2_] := (s[i1, i2 + 1] - s[i1, i2])/dx;
     (*return average free energy density function*)
     Sum[j (dsdx[i1, i2].dsdx[i1, i2] + dsdy[i1, i2].dsdy[i1, i2]) + a s[i1, i2][[3]]^2 - 
        h s[i1, i2][[3]], {i1, 1, n}, {i2, 1, n}]/n^2]]
  ]

(See the note at the bottom of this post for explanation of the mem thing.)