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
Attachments:
