3
$\begingroup$

I have a workable function with using NDSolve with predefined method for solving. But for some reasons I should use the parallel computation. The main problem is that the NDSolve can't build the structure for solving this problem with defined method, but if we remove the definition of method - it works.

The code of solving the problem is below:

 DOPRI54[] := 
  Module[{DOPRIamat, DOPRIbvec, DOPRIcvec, DOPRIevec, 
    DOPRICoefficients},
   DOPRIamat := {{1/5}, {3/40, 9/40}, {44/45, -56/15, 
      32/9}, {19372/6561, -25360/2187, 
      64448/6561, -212/729}, {9017/3168, -355/33, 46732/5247, 
      49/176, -5103/18656}, {35/384, 0, 500/1113, 125/192, -2187/6784,
       11/84}}; 
   DOPRIbvec := {35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0}; 
   DOPRIcvec := {1/5, 3/10, 4/5, 8/9, 1, 1}; 
   DOPRIevec := {71/57600, 0, -71/16695, 71/1920, -17253/339200, 
     22/525, -1/40}; 
   DOPRICoefficients[5, p_] := 
    N[{DOPRIamat, DOPRIbvec, DOPRIcvec, DOPRIevec}, p];
   Return[{"ExplicitRungeKutta", "DifferenceOrder" -> 5, 
     "Coefficients" -> DOPRICoefficients, "StiffnessTest" -> False}]
   ];

Fehlberg45[] := 
  Module[{Fehlbergamat, Fehlbergbvec, Fehlbergcvec, Fehlbergevec, 
    FehlbergCoefficients},
   Fehlbergamat := {{1/4}, {3/32, 9/32}, {1932/2197, -7200/2197, 
      7296/2197}, {439/216, -8, 3680/513, -845/4104}, {-8/27, 
      2, -3544/2565, 1859/4104, -11/40}}; 
   Fehlbergbvec := {25/216, 0, 1408/2565, 2197/4104, -1/5, 0}; 
   Fehlbergcvec := {1/4, 3/8, 12/13, 1, 1/2}; 
   Fehlbergevec := {-1/360, 0, 128/4275, 2197/75240, -1/50, -2/55}; 
   FehlbergCoefficients[4, p_] := 
    N[{Fehlbergamat, Fehlbergbvec, Fehlbergcvec, Fehlbergevec}, p];
   Return[{"ExplicitRungeKutta", 
     "Coefficients" -> FehlbergCoefficients, "DifferenceOrder" -> 4, 
     "EmbeddedDifferenceOrder" -> 5, "StiffnessTest" -> False}]
   ];

BS54[] := 
  Module[{}, 
   Return[{"ExplicitRungeKutta", 
     "Coefficients" -> "EmbeddedExplicitRungeKuttaCoefficients", 
     "DifferenceOrder" -> 5, "StiffnessTest" -> False}]];

ClassicalRK4[] := 
 Module[{crkamat, crkbvec, crkcvec, ClassicalRungeKuttaCoefficients},
  crkamat = {{1/2}, {0, 1/2}, {0, 0, 1}};
  crkbvec = {1/6, 1/3, 1/3, 1/6};
  crkcvec = {1/2, 1/2, 1};
  ClassicalRungeKuttaCoefficients[4, p_] := 
   N[{crkamat, crkbvec, crkcvec}, p];
  Return[{"ExplicitRungeKutta", "DifferenceOrder" -> 4, 
    "Coefficients" -> ClassicalRungeKuttaCoefficients}]
  ]

getMethods[methodname_] := Module[{method = ""},
  Switch[methodname, "DOPRI54", method = DOPRI54[], "Fehlberg45", 
   method = Fehlberg45[], "BS54", method = BS54[], "ClassicalRK4", 
   method = ClassicalRK4[]];
  Return[method]
  ]

getEquilibriumPoints[system_] := Block[{sol},
  sol = NSolve[{system == {0., 0., 0.}}, {x[t], y[t], z[t]}];
  Return[{x[t], y[t], z[t]} /. sol]
  ]

The problem part of the code (it has no problem for solving without parallelization):

gettingInLocality[system_, pointsCenters_, radius_, T_, u0_, 
  OptionsPattern[methodname -> "DOPRI54"]] := 
 Module[{set = {}, events, method},
  method = getMethods[OptionValue[methodname]];
  events = {(x[t] - #[[1]])^2 + (y[t] - #[[2]])^2 + (z[t] - #[[
            3]])^2 <= radius^2} & /@ pointsCenters;

  NDSolve[{{x'[t], y'[t], z'[t]} == 
       system, {x[0], y[0], z[0]} == #}, {x[t], y[t], z[t]}, {t, 0, 
      T}, Method -> {"EventLocator", 
       "Event" -> First /@ ## &[events], 
       "EventAction" :> 
        Throw[set = Append[set, {#, t}], "StopIntegration"], 
       "Method" -> method}] & /@ u0;
  Return[set]
  ]

The version of code with using parallelization:

gettingInLocality[system_, pointsCenters_, radius_, T_, u0_, 
  OptionsPattern[methodname -> "DOPRI54"]] := 
 Module[{set = {}, events, method},
  SetSharedVariable[set];
  method = getMethods[OptionValue[methodname]];
  events = {(x[t] - #[[1]])^2 + (y[t] - #[[2]])^2 + (z[t] - #[[
            3]])^2 <= radius^2} & /@ pointsCenters;

  ParallelMap[
   CriticalSection[{setClock}, 
     NDSolve[{{x'[t], y'[t], z'[t]} == 
        system, {x[0], y[0], z[0]} == #}, {x[t], y[t], z[t]}, {t, 0, 
       T}, Method -> {"EventLocator", 
        "Event" -> First /@ ## &[events], 
        "EventAction" :> 
         Throw[set = Append[set, {#, t, $KernelID}], 
          "StopIntegration"], "Method" -> method (*If we remove that part of code with choosing method - it works*)}]] &, u0];
  Return[set]
  ]

And the last part of code for solving the problem:

sigma  = 10; rho = 28 ; beta = N[8/3];
initdata = {19., -25., 10.};
system = {sigma (y[t] - x[t]), x[t] (rho - z[t]) - y[t], 
   x[t] y[t] - beta z[t]};
pointsCenters = getEquilibriumPoints[system][[{1, 2}]]
radius = 10;
T = 10.^3;
u0 = Flatten[
  Table[{x, y, z}, {x, -25., 25., 25.}, {y, -25., 25., 30.}, {z, 0., 
    50., 15.}], 2]
set = {};
set = gettingInLocality[system, pointsCenters, radius, T, u0];

The main problem is an errors about not correct structure of coefficients and initialization of the method failed. But if we remove the marked part of code (look parallel version of getting local) it works. In this case the choosing of the appropriate method is very important. How I can fix this?

The image of errors is below:

enter image description here

$\endgroup$

closed as off-topic by Michael E2, MarcoB, m_goldberg, bbgodfrey, LCarvalho Feb 26 at 16:47

  • The question does not concern the technical computing software Mathematica by Wolfram Research. Please see the help center to find out about the topics that can be asked here.
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 1
    $\begingroup$ I did not find any errors. Both versions of the code work, but when using ParallelMap list set probably contains kernel number. $\endgroup$ – Alex Trounev Feb 23 at 15:19
  • 4
    $\begingroup$ I'm voting to close this question as off-topic because there is no error when executed and both methods produce equivalent output. $\endgroup$ – Michael E2 Feb 23 at 20:44
4
$\begingroup$

DOPRI54[] returns a list of rules containing a symbol DOPRICoefficients$<<some number>> that does not exist anymore because it was scoped by Module that has the attribute Temporary. Maybe temporary symbols are not automatically shared by Parallel and its ilk.

I believe that it would be a much better strategy to define DOPRI54 (and many of the other symbols but I focus on this one) in a more encapsulated way.

In the followng, the actual function that you want to submit as option is a pure function (defined by the infix form \[Function] of Function). In order to inline the definitions DOPRIamat, DOPRIbvec, DOPRIcvec, and DOPRIevec, you may use With. Look at this:

ClearAll[DOPRI54];
DOPRI54 = With[{
   DOPRIamat = {{1/5}, {3/40, 9/40}, {44/45, -56/15, 32/9}, {19372/6561, -25360/2187, 64448/6561, -212/729}, {9017/3168, -355/33, 46732/5247, 49/176, -5103/18656}, {35/384, 0, 500/1113, 125/192, -2187/6784, 11/84}},
   DOPRIbvec = {35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0},
   DOPRIcvec = {1/5, 3/10, 4/5, 8/9, 1, 1},
   DOPRIevec = {71/57600, 0, -71/16695, 71/1920, -17253/339200, 22/525, -1/40}
   },
   {
    "ExplicitRungeKutta",
    "DifferenceOrder" -> 5,
    "Coefficients" -> ({order, p} \[Function] N[{DOPRIamat, DOPRIbvec, DOPRIcvec, DOPRIevec}, p]),
   "StiffnessTest" -> False
    }]

I have not tested this but I believe that this is much better encapsulated. This is also essential for preventing the parallel kernels from calling the main kernel (which is essential for scalability).

Moreover, make sure to learn the difference between SetDelayed (:=) and Set (=). You really want to use Set here.

$\endgroup$
  • $\begingroup$ I suspect the Module-localized variables, while Temporary, stick around for a while in the main kernel (IIRC, they persist while there are references to them) ; however, they are not distributed to the subkernels. $\endgroup$ – Michael E2 Feb 23 at 3:40
  • $\begingroup$ You are right! I totally forgot about that. $\endgroup$ – Henrik Schumacher Feb 23 at 8:20
  • $\begingroup$ Thank you @MarcoB for your edit. I really appreciate that! $\endgroup$ – Henrik Schumacher Feb 23 at 17:28
  • $\begingroup$ Does the OP's code generate errors for you? I finally tried it, and it works fine as is. $\endgroup$ – Michael E2 Feb 23 at 20:45
  • $\begingroup$ @MichaelE2 Huh, indeed, it does not! I could have promised that I had tried it yesterday... $\endgroup$ – Henrik Schumacher Feb 23 at 20:58

Not the answer you're looking for? Browse other questions tagged or ask your own question.