# ParallelDo errors using Physical Constant and a list

I am new to using Mathematica in parallel and need help with the following code. When it is not run in parallel as follow, it works.

Remove["Global*"]
Needs["PhysicalConstants"]

num = 5;
lista = {};

sis := (
x0 = RandomReal[];
v0 = First[SpeedOfSound];

sol = First[x /. NDSolve[{x''[t] + x[t] == 0, x[0] == x0, x'[0] == v0}, x, {t, 1}]];
inte = Integrate[sol[t], {t, 0, 1}];
AppendTo[lista, {inte}]
)

Do[sis, {num}];
lista


But, when it is run in parallel using ParallelDo and SetSharedVariable[lista] it generates errors.

Remove["Global*"]
Needs["PhysicalConstants"]

num = 5;
lista = {};

sis := (
x0 = RandomReal[];
v0 = First[SpeedOfSound];

sol = First[
x /. NDSolve[{x''[t] + x[t] == 0, x[0] == x0, x'[0] == v0},
x, {t, 1}]];
inte = Integrate[sol[t], {t, 0, 1}];

AppendTo[lista, {inte}])

SetSharedVariable[lista];
ParallelDo[sis, {num}];
lista


The error has several lines. Could somebody tell me what I need to change in the code above to make it work?

• PhysicalConstants is obsolete. You can use the new built-in Quantity["SpeedOfSound"] instead of PhysicalConstantsSpeedOfSound. See the transition guide here. – MarcoB Aug 12 '15 at 3:43
• Your parallel code runs fine the first time you execute it on a fresh kernel, but fails when I try to re-run it (Set::wrsym: Symbol lista is Protected.). Is that what you observe as well? It would be useful if you could confirm or add the exact error you are experiencing to your question. – MarcoB Aug 12 '15 at 3:44
• See here for why it is a bad idea in general to use AppendTo and SetSharedVariable. – Szabolcs Aug 12 '15 at 7:55

To start, if you need a package for your parallel computation, you need to load that package on all kernels using ParallelNeeds instead of Needs. As a start, in your case using the following at least takes care of the undefined symbol errors you were receiving.

ParallelNeeds["PhysicalConstants"]


Furthermore, the speed of sound and other physical constants can be obtained from the built-in Quantity system since version 9, (see Upgrading From Physical Constants). Admittedly, the built-in system is still clunkier than I would like, but it does give access to the numerical values of physical constants in an almost limitless variety of units of measurement:

speedofsound = QuantityMagnitude@UnitConvert[Quantity["SpeedOfSound"], "m/s"];


I assume that you are aware of the fact that your differential equation can be solved analytically, so I suppose that this might be just a simplified example similar, but not identical, to your application.

If that is not the case, then you should definitely consider carrying out most of the math symbolically here:

Clear[f]
f[x0_] = Integrate[
x[t] /. First@DSolve[{x''[t] + x[t] == 0, x[0] == x0, x'[0] == v0}, x[t], t],
{t, 0, 1.}
] /. v0 -> speedofsound

Map[f, RandomReal[1, 5]]; // AbsoluteTiming

(* Out: {0.0000349718, Null} *)


In this case, the calculation is so trivial that I saw a loss in performance when attempting to parallelize execution, even for very large numbers (e.g. 5 million random values), at least on my two-kernel laptop.

In those cases in which you have to solve the differential equation numerically, then it is still simpler than your approach to use ParametricNDSolveValue with $x_0$ as a parameter.

Clear[integral]

integral[param_] :=
NIntegrate[
ParametricNDSolveValue[
{x''[t] + x[t] == 0, x[0] == x0, x'[0] == speedofsound},
x, {t, 0, 1}, x0
][param][t],
{t, 0, 1}
]

Map[integral, RandomReal[1, 5000]]; // AbsoluteTiming
ParallelMap[integral, RandomReal[1, 5000]]; // AbsoluteTiming

(* Out:
{19.5816, Null}
{9.95256, Null}
*)


The speedup from parallelization is appreciable, even on my laptop with only two kernels. I assume that it would get better if your ran more examples. In order to fully appreciate its effects, you will want to time a parallel run after your first, so you don't include the time it takes the system to start the parallel kernels.

In order to answer the question as you asked it, though, it is worth mentioning that AppendTo is generally a bad idea performance-wise (e.g. it is slow, and it creates a copy of the list to be appended to). If you want to collect an undefined number of results from your computation, it is more effective to use a combination of Sow and Reap. See for example:

There is a small complication when using Reap and Sow in parallel, as explained in this tutorial from WRI: Sow, Reap, and Parallel Programming, and in questions on this site (e.g. Reap, Sow with Parallelize: bad performance, why?, and How to collect result continuously (interruptible calculation) when running parallel calculations?). In short, one always needs to Sow on the master kernel.

This is accomplished below by following the advice of that tutorial and creating ParallelSow, a wrapper for Sow that is specified to always run on the master kernel using SetSharedFunction.

Here is a version of your code that uses Sow and Reap and will run in parallel:

Clear[ParallelSow]
ParallelSow[expr_] := Sow[expr]
SetSharedFunction[ParallelSow]

num = 15;

lista = ParallelDo[
Module[
{x0, sol, inte},
x0 = RandomReal[];
sol = First[x /. NDSolve[{x''[t] + x[t] == 0, x[0] == x0, x'[0] == speedofsound}, x, {t, 1}]];
inte = Integrate[sol[t], {t, 0, 1}];
ParallelSow[inte]
],
{num}
] // Reap;

lista[[2, 1]]

(* Out:
{156.49, 156.987, 156.926, 156.996, 156.684, 156.564, 156.436, 157.219, 156.972, 156.847,
157.153, 156.851, 156.522, 156.784, 157.088}
*)
`