Message Boards Message Boards

How do I tune my NDSolve code to start fast, so it can be run hundreds of times?

Hi!

I'm simulating a monolayer epithelia biological system, in which the cells behave fluid-like and frequently exchange neighbors. This neighbor exchange causes a problem, because it changes the differential equations which govern the tissue. So, I have a WhenEvent set up for each edge within the tissue, to track if it becomes short enough for the cells to exchange neighbors. Then, when such an event occurs, the code stops the NDSolve block, runs the neighbor exchange, then restarts.

But, this introduces huge overhead that I can't figure out the source of; I have a tracker keeping track of the time step that the solver is on, which stalls for significant amounts of time on each startup of NDSolve. I'm not sure what the program is doing at this time. Here's the code (only necessary modifications are changing the directories for the attatchments which I've included on this post)

SetDirectory[NotebookDirectory[]];
a = 225;
b = 225;
x0 = 0;
y0 = 0;
d = Sqrt[a^2 - b^2];
center = {x0, y0};
focus1 = {x0 - d, y0};
focus2 = {x0 + d, y0};
nCells = 50;
scale = Sqrt[
   41.3083/(Pi*a*b/nCells)]; (*one length unit equals scale*micron*)
originalNCells = nCells;
ellipsePointList = 
  Table[{a*Cos[\[Theta]], b*Sin[\[Theta]]}, {\[Theta], 0, 2 \[Pi], 
    2 \[Pi]/200}];
cellularVoronoiDiagram = 
  Import["C:\\Users\\kaden\\Downloads\\savedvoronoidiagram.m"];
tAnneal = 750;
tHealing = 750;
\[Mu] = 30;
timeConstant = 
  0.8969476207347036;(*experimentally determined!!!!!!!!!! time \
should now be measured in seconds*)
timeDifference = 
  timeConstant/
   8; (*this is a meaningless variable but I'm affraid if I delete it \
something will break*)
intercolationProportion = 0.00001;
nonDimensionalizedContractileConstant = 0.01771;
outerTension = 0.14;
lineTensionVariability = 0.01;
decayRate = 150;
Import["C:\\Users\\kaden\\Downloads\\voronoiDiagramFunctions.m"]
contractileTension = 
  Table[nonDimensionalizedContractileConstant, {i, 
    Length[cellularVoronoiDiagram[[2]]]}];
meanTensionCell = 
  Table[(contractileTension[[i]]/
      originalScale)*(originalPerimeterScale - P0s[[i]]), {i, 
    Length[P0s]}];
meanTensionEdge = 
  Table[Total[
    meanTensionCell[[getCellsAdjacentToEdge[edgeData[[All, 2]][[i]], 
       cellularVoronoiDiagram]]]], {i, Length[edgeData[[All, 2]]]}];

tMax = Max[{tHealing, tAnneal}] + 200;

lineTensions = 
  Table[proc = 
    ItoProcess[\[DifferentialD]lineTension[
        t] == -1/
         decayRate*(lineTension[t] - 
          meanTensionEdge[[i]]) \[DifferentialD]t + \[DifferentialD]w[
         t], lineTension[t], {lineTension, meanTensionEdge[[i]]}, t, 
     w \[Distributed] WienerProcess[0, lineTensionVariability]];
   Interpolation[
     MovingAverage[RandomFunction[proc, {0, tMax, 0.5}], 5]][
    t + 10], {i, Length[edgeData]}];

healths = Table[1, {i, Length[cellularVoronoiDiagram[[2]]]}];
Print["linetensionscreated"];
nVertices = Length[vertexList];
tCurrent = 0;
Subscript[vertexLists, 0] = vertexList;
Subscript[timeStop, 0] = 0;
Subscript[cellularVoronoiDiagrams, 0] = cellularVoronoiDiagram;
Subscript[edgeDatas, 0] = edgeData;
stepSize = 0.0000001;
SetSystemOptions[
  "NDSolveOptions" -> "DefaultSolveTimeConstraint" -> 10.`];
solutions = {};
dummies = Table[Subscript[L, i], {i, Length[lineTensions]}];
substitutionsLineTensions = 
  Table[Subscript[L, i] -> lineTensions[[i]], {i, 
    Length[lineTensions]}];
Timing[Do[
  (*odeList=D[Flatten[vertexList],t] == Flatten[
  ParallelTable[(1/\[Mu])*forceOnVertex[i,cellularVoronoiDiagram,
  vertexList,A0s,contractileTension,Table[0,{i,Length[edgeData]}],
  edgeData,healths],{i,Length[vertexList]}]];*)
  diffEQs = 
   D[Flatten[vertexList], t] == 
    Flatten[Table[(1/\[Mu])*
        forceOnVertexSolved[i, cellularVoronoiDiagram, vertexList, 
         A0s, contractileTension,(*lineTensions*)(*meanTensionEdge*)
         dummies, edgeData, healths] /. substitutionsLineTensions, {i,
        Length[vertexList]}]];
  bcList = 
   Flatten[vertexList /. t -> tCurrent] == 
    Flatten[cellularVoronoiDiagram[[1]]];
  (*Check working precision*)
  Timing[times = {};
   Monitor[
    s = NDSolve[{diffEQs, bcList, 
       With[{edgeData = edgeData, vertexList = vertexList}, 
        WhenEvent[# - intercolationLength < 0, 
           (cellularVoronoiDiagram[[1]] = 
             Table[{vertexList[[i, 1]], vertexList[[i, 2]]}, {i, 
               Length[vertexList]}];

            intercolatedEdge = 
             Position[edgeData[[All, 1]], #, Heads -> False][[1]][[1]];
            Print["intercolating ", intercolatedEdge];
            tCurrent = t;
            "StopIntegration"
            )] & /@ edgeData[[All, 1]]]}, 
      Flatten[vertexList], {t, tCurrent, tAnneal}, 
      StartingStepSize -> stepSize, 
      EvaluationMonitor :> AppendTo[times, t], AccuracyGoal -> 2, 
      PrecisionGoal -> 2], times[[-1]] ];
   ];
  stepSize = (Subscript[x, 1][t] /. s[[1]] /. {t -> "Coordinates"} // 
       First // Differences)[[-1]];
  AppendTo[solutions, s];
  numIntercolations = j - 1;
  Subscript[timeStop, j] = times[[-1]];
  If[tAnneal - 0.01 < times[[-1]] < tAnneal + 0.01, Break[]];
  tempCenter = 
   Mean[{cellularVoronoiDiagram[[1]][[edgeData[[intercolatedEdge]][[2]\
][[1]]]], 
     cellularVoronoiDiagram[[1]][[edgeData[[intercolatedEdge]][[2]][[\
2]]]]}]; 
  newEdge = 
   r[((cellularVoronoiDiagram[[1]][[edgeData[[intercolatedEdge]][[2]][\
[1]]]] - 
         cellularVoronoiDiagram[[1]][[edgeData[[intercolatedEdge]][[2]\
][[2]]]])/
       Norm[cellularVoronoiDiagram[[1]][[edgeData[[intercolatedEdge]][\
[2]][[1]]]] - 
         cellularVoronoiDiagram[[1]][[edgeData[[intercolatedEdge]][[2]\
][[2]]]]])*(intercolationLength + 1)];
  intercolate[edgeData[[intercolatedEdge]]]; 
  cellularVoronoiDiagram[[1]][[edgeData[[intercolatedEdge]][[2]][[1]]]\
] = (tempCenter + 0.5*newEdge ); 
  cellularVoronoiDiagram[[1]][[edgeData[[intercolatedEdge]][[2]][[2]]]\
] = (tempCenter - 0.5*newEdge);
  Subscript[cellularVoronoiDiagrams, j] = cellularVoronoiDiagram;
  Subscript[edgeDatas, j] = edgeData;
  Subscript[vertexLists, j] = vertexList,

  {j, 40000}]]

nCells = Length[cellularVoronoiDiagram[[2]]];
solutionTable = 
  Table[vertexList /. solutions[[i]][[1]], {i, Length[solutions]}];
cellularVoronoiDiagram[[1]] = 
  solutionTable[[numIntercolations + 1]] /. t -> tAnneal;

Thanks for any and all help!!

Best, Kaden

POSTED BY: Kaden Tro
3 Replies

I haven't run your code, but NDSolve is a very high level function that does a large amount of pre-evaluation of the problem before actually integrating, in which it is picking the best method, etc. I notice that you don't specify a method, so you may be triggering that whole process every time you restart. I suggest you experiment with specifying as many of the NDSolve options that you can, but especially the method.

POSTED BY: Gareth Russell

This worked great, thank you so much!

Another thing I tried was compiling the interpolated random walks which are coefficients in the differential equations. It seems to make things faster, but in trials with the exact same seed I end up with less triggers on my whenevent code, and I have no idea why that would be--I don't think the compilation should affect the solution at all, but it seems to be doing just that. Any ideas why that might be?

POSTED BY: Kaden Tro

Sorry, I am not sure about that. Generally I would look at precision issues, as the compiled version will be limited to machine precision.

POSTED BY: Gareth Russell
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract