Is it possible to speed up ContourPlot on multi-core machines?

I second @Verbeia's suggestion: compute the function on a mesh of points and use ListContourPlot. The disadvantage is that ListContourPlot has no adaptive sampling, so it'd be preferable if we could do our own adaptive sampling somehow. Adaptive sampling can give you a much better result while needing to compute the function in far less points---and the problem here is indeed computation time. So ContourPlot with its adaptive sampling might give a better result in less time on a single CPU than ListContourPlot will with a high resolution mesh computed on many CPUs.

Adaptive sampling is what I asked about (and solved) here: Adaptive sampling for slow to compute functions in 2D

The method I implemented there is usable (I am using it for something very similar to what you describe) but it is not nearly as good as ContourPlot's own. So one might still try to somehow make use of it. I'm quoting one suggestion I received from Leonid Shifrin there (in a comment):

You probably can control the DensityPlot, although not directly. Since it calls your function, you can simply Sow the values until some criteria (which you define) is violated (or satisfied). Then, you stop via throwing an exception, and catching it in the outer function, but still inside Reap. Alternatively, you could just start fooling DensityPlot by supplying faked values (perhaps, interpolated, or whatever), and it will stop by itself, I guess. Not sure this will work for you, but it may be worth trying.

I have not tried to implement this before, but I think it could work if your function is sufficiently smooth (which mine is definitely not, but yours may be).

Here's a quick sample implementation of how it could work:

First, let's define a sample function to plot:

fun[{x_, y_}] := 1/(1 + Exp[10 (Norm[{x, y}] - 3)])

Let's divide both the $x$ and $y$ axes into 5 parts on the interval $[0,5]$ and generate a mesh of points:

initialDivision = Range[0, 5];

points = N@Tuples[initialDivision, {2}];

Calculate function values on the intial mesh. This can be parallelized (just use ParallelMap)

values = fun /@ points;

This counter i will be used to control the maximal subdivisions in ContourPlot:

i = 0;

Now put the following code into a single cell, and evaluate it several times. Each time a finer and finer approximation will be computed. The points where function values have been computed will also be visualized. Note that I fixed the plot points in ContourPlot to force it to use the same initial mesh that I used, and I also fixed the number of contours.

if = Interpolation@ArrayFlatten[{{points, List /@ values}}]

{plot, {newpoints}} = Reap[
   ContourPlot[if[x, y], {x, 0, 5}, {y, 0, 5}, 
    Contours -> Range[0, 1, .1], MaxRecursion -> (++i), 
    PlotPoints -> Length[initialDivision], 
    EvaluationMonitor :> Sow[{x, y}]]
   ];
plot

newpoints = Complement[newpoints, points];
newvalues = fun /@ newpoints;  (* <-- this can be parallelized *)
points = Join[points, newpoints];
values = Join[values, newvalues];

Graphics[Point[points]]

After a few iterations the contour plot and the point mesh will look like this (note that the code above only plots the contours for the previous step, not the current results):

After 3 iterations, this method has computed the function value in 3809 points for this particular function.

Let's compare this with a plain ContourPlot using the same parameters:

ContourPlot[fun[{x, y}], {x, 0, 5}, {y, 0, 5}, 
    PlotPoints -> 6, MaxRecursion -> 3]

The quality of the plot is about the same with a plain ContourPlot as well.

How many points did the plain CountoutPlot use?

Reap[ContourPlot[fun[{x, y}], {x, 0, 5}, {y, 0, 5}, PlotPoints -> 6, 
    MaxRecursion -> 3, EvaluationMonitor :> Sow[{x, y}]]][[2, 1]] // Length

(* ==> 3790 *)

It uses almost the same number of points, so if the bottleneck is computing f, the method I described is going to be almost as fast as ContourPlot on a single core, with the advantage that it is parallelizable for multiple cores.

The next step would be packaging this up into a self-contained function, but seeing how the quality improves step by step is also valuable as you can make decisions about when to stop calculating (and avoid excessive computation times).

---

I find it quite disappointing that all those nice and fast algorithms that plotting functions use (fast Voronoi cells, Delaunay trinagulation, adaptive sampling) are not directly accessible by users. We either have to use hacks to access these algorithms or reimplement them.

One idea, which was originally appealing to me, was to simply split the plot domain in four quarter size pieces and have these pieces calculated in parallel using ParallelTable and then combined using Show. It appears though that the overhead of getting the graphics data back to Show is pretty large, so this only yields some extra speed if the function is computationally slow.

a = 36;
(parPlot = Show[
     ParallelTable[
      RegionPlot[
       Mod[Sqrt[x^2 + y^2] - 7/2 ArcTan[x, y] +Sin[x] +Cos[y], \[Pi]] < \[Pi]/2, 
       {x, -a + i a, a + i a}, {y, -a + i a, a + i a}, 
       PlotRange -> {{-a, a}, {-a, a}}, PlotPoints -> 50], 
       {i, 0, 1}, {j, 0, 1}]
     ]); // AbsoluteTiming

{25.5344605, Null}

parPlot 

The original serial calculation:

a = 36;
g = RegionPlot[
    Mod[Sqrt[x^2 + y^2] - 7/2 ArcTan[x, y] + Sin[x] + 
       Cos[y], \[Pi]] < \[Pi]/2, {x, -a, a}, {y, -a, a}, 
    PlotPoints -> 100]; // AbsoluteTiming

(*
==> {51.2839332, Null}
*)

g

It's twice as slow, but it looks slightly better in places. RegionPlot is worse than ContourPlot in this sense.

Try generating a mesh of points using ParallelMap and then plot them using ListContourPlot.

You might find the techniques in this question on adaptive sampling for hard-to-compute functions useful.