I need to solve several stoke's equations for very similar boundary conditions, which represent particles in a fluid with different velocities U[i]. I'm interested in computing the force acting on the particles as a function of their velocities, and in the real case will be several particles and for several configurations. This is why I need to speed up the code wherever is possible.

This is the code that I'm using:

L = 10;
cords = {0, 0};

<< NDSolve`FEM`

r = RegionDifference[Rectangle[{-L, -L}, {L, L}], Disk[cords, 1]];

mesh2 = ToElementMesh[ r, "MaxBoundaryCellMeasure" -> .1, 
   "MaxCellMeasure" -> .75];

I need the options for the mesh since the stress tensor depends on the derivatives of the fluid velocity, thus I want a smooth function.

The stokes equations with zero boundary conditions:

eqn = {
   Laplacian[u[x, y], {x, y}] - D[p[x, y], x] == 0,
   Laplacian[v[x, y], {x, y}] - D[p[x, y], y] == 0,
   D[u[x, y], x] + D[v[x, y], y] == 0};

(*with null boundary conditions in "infinite"*)

bcs = {u[x, -L] == 0, u[x, L] == 0,
   u[-L, y] == 0, u[L, y] == 0,
   v[x, -L] == 0, v[x, L] == 0,
   v[-L, y] == 0, v[L, y] == 0,
   p[L, y] == 0, p[-L, y] == 0};

And now a loop on different velocities for the particles.

plots = Table[{}, {i, 4}];
forceX =  Table[{}, {i, 4}];
 U = {1 + i/10., 0};

 pde = {eqn, bcs, 
   DirichletCondition[u[x, y] == U[[1]], Norm[{x, y}]^2 <= 1.], 
   DirichletCondition[v[x, y] == U[[2]], Norm[{x, y}]^2 <= 1.]};

 {xVel, yVel, pressure} = 
  NDSolveValue[pde, {u, v, p}, {x, y} \[Element] mesh2, 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];

 plots[[i]] = ContourPlot[yVel[x, y], {x, y} \[Element] mesh2];

 Subscript[\[Sigma], yy][x_, y_] = Evaluate[2 D[yVel[x, y], {y}]];
 Subscript[\[Sigma], xx][x_, y_] = Evaluate[2 D[xVel[x, y], {x}]];
 Subscript[\[Sigma], xy][x_, y_] = 
  Evaluate[D[xVel[x, y], {y}] + D[yVel[x, y], {x}]];

 forceX[[i]] = 
  Quiet@NIntegrate[(-pressure[Cos[x], Sin[x]] + 
        Subscript[\[Sigma], xx][Cos[x], Sin[x]]) Cos[x] + 
     Subscript[\[Sigma], xy][Cos[x], Sin[x]] Sin[x], {x, 0, 
     2 \[Pi]}], {i, 4}]

(*{-9.8274, -10.7208, -11.6142, -12.5076}*)

I obtain like this the force on the particle as a function of U. But if you look at plots, one can see that all the solutions are essentially the same, save a scale factor:

enter image description here

So the question is, can I help NDSolve to find the solution for U[i+1] giving it the solution previously found for U[i]? For what I understand of FEM this is a minimization problem, so giving it a better starting point should reduce the number of iterations significantly.

  • $\begingroup$ How do you deduce that the solutions are the same? They look different to me and not only due to scale (the scale is same in all of them, right?). Seems like the effect of boundary conditions is reflected differently at areas that you actually have the highest gradients. $\endgroup$
    – MathX
    Apr 2, 2016 at 3:57
  • $\begingroup$ same functional form, with four lobes and decaying on the borders. The force is different in each one too, last comment of the code. My point is, it's easier to start searching from one of those functions than start from a constant or a linear function. Plus in the "real" case I'm interested in small variations of U so the solutions are even more similar than in this example case. $\endgroup$ Apr 2, 2016 at 4:29
  • $\begingroup$ Yes your idea is in line with modifying the initialization. I honestly don't know if you can do it in a general case. I mean how would you know without solving the problem once? Except maybe from some kind of simplified theoretical approach that could be problem dependent. However, I'll just trow in an idea, although I guess you already know about it. $\endgroup$
    – MathX
    Apr 2, 2016 at 17:33
  • $\begingroup$ Maybe you can solve it for a range of initial values and then form a function of unknown type using FindFormula to get an estimated solution as a function of initial values and then use that? $\endgroup$
    – MathX
    Apr 2, 2016 at 17:37
  • $\begingroup$ to answer your first comment, you have the theoretical solution for U in the case you have only one particle. The assumption I will make is that the solution for several particles is a superposition of the single particle solution. Your second comment is interesting, I was thinking to use Newton's method or steepest descent but maybe the function is well behaved and has just one minimum easier to find that way. I still like to be able to give NDSolve a starter function though. $\endgroup$ Apr 2, 2016 at 23:31

1 Answer 1


Here is how I'd do it. The idea is to use the low level FEM functions and to discretize the PDE just once. All stationary boundary conditions are deployed. Additional there are two parametric Dirichlet boundary conditions that can change during the loop. Only those are re-discretized and re-deployed in every loop step. Have a look at the FEM programming tutorial and the FEM refpages.

L = 10;
\[Mu] = 1;
center = {0, 0};
r = RegionDifference[Rectangle[{-L, -L}, {L, L}], Disk[center, 1]];

(* setup of numerical region, variable and solution data; note two
parameter p1 and p2 in the variable data and initial values for those
in the solution data *)
nr = ToNumericalRegion[r];
vd = NDSolve`VariableData[{"DependentVariables", "Space", 
     "Parameter"} -> {{u, v, p}, {x, y}, {p1, p2}}];
sd = NDSolve`SolutionData[{"Space", "Parameter"} -> {nr, {1.1, 0.}}];

(* init stokes coefficients *)
mu = \[Mu]*IdentityMatrix[2];
zero = 0*IdentityMatrix[2];
iPDE = InitializePDECoefficients[vd, sd, "DiffusionCoefficients" -> {
     {-mu, zero, zero},
     {zero, -mu, zero},
     {zero, zero, zero}
   "ConvectionCoefficients" -> {
     {{{0, 0}}, {{0, 0}}, {{1, 0}}},
     {{{0, 0}}, {{0, 0}}, {{0, 1}}},
     {{{1, 0}}, {{0, 1}}, {{0, 0}}}

(* init BCs coefficients; note two DBCs with parameter p1 and p2 *)
iBCs = InitializeBoundaryConditions[vd, sd, {
    {DirichletCondition[u[x, y] == 0, 
      x == L || x == -L || y == L || y == -L], 
     DirichletCondition[u[x, y] == p1, x^2 + y^2 <= 1.]},
    {DirichletCondition[v[x, y] == 0, 
      x == L || x == -L || y == L || y == -L], 
     DirichletCondition[v[x, y] == p2, x^2 + y^2 <= 1.]},
    {DirichletCondition[p[x, y] == 0., x == L || x == -L]}

(* setup method data and possibly mesh specifications *)
md = InitializePDEMethodData[vd, sd, 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1},
     "MeshOptions" -> {"MaxBoundaryCellMeasure" -> .1, 
       "MaxCellMeasure" -> .75}}];

(* discretize the PDE and the "Stationary" boundary conditions *)
dPDE = DiscretizePDE[iPDE, md, sd];
dBCs = DiscretizeBoundaryConditions[iBCs, md, sd, "Stationary"];
load = dPDE["LoadVector"];
stiffness = dPDE["StiffnessMatrix"];
(* deploy stationary BCs *)
DeployBoundaryConditions[{load, stiffness}, dBCs];

(* pre-compute how the solution vector is split *)
spans = MapThread[
   Span, {Most[md["IncidentOffsets"] + 1], 

(* In  every loop step update the parametric solution data to contain
the current parametric values; copy the system matrices and deploy
the discretized parametric BCs; then LinearSolve *)
 NDSolve`SetSolutionDataComponent[sd, "Parametric", {1 + i/10., 0}];
 {lcopy, scopy} = {load, stiffness};
 dBCs2 = DiscretizeBoundaryConditions[iBCs, md, sd, "Parametric"];
 DeployBoundaryConditions[{lcopy, scopy}, dBCs2];
 res = LinearSolve[scopy, lcopy];
 , {i, 1, 4}];

(* visualize the last result *)
mesh = nr["ElementMesh"]; mesh1 = MeshOrderAlteration[mesh, 1];
xVel = ElementMeshInterpolation[{mesh}, res[[spans[[1]]]]];
yVel = ElementMeshInterpolation[{mesh}, res[[spans[[2]]]]];
pressure = ElementMeshInterpolation[{mesh1}, res[[spans[[3]]]]];
{ContourPlot[xVel[x, y], {x, y} \[Element] mesh], 
 ContourPlot[yVel[x, y], {x, y} \[Element] mesh], 
 ContourPlot[pressure[x, y], {x, y} \[Element] mesh1]}

enter image description here

Now, in principal one could go even further and eliminate the repeated calls to LinearSolve, but that would requite a little more thought. See how far you get with this and in case it's not sufficient that could be tried as well.

One other thing, however, that might make more sense is to talk a bit about this integral boundary (?) condition you want to have. If you could put forward an explanation for that then it might be possible to implement that in stead.

  • $\begingroup$ This is great, thanks! I saw some of those low-level programing examples but they are quite obscure, I could have never come up with something like this. Will work tomorrow on implementing the minimization algorithm with this and post it as an answer to the next question I'm making about the integral boundary condition. Thanks again. $\endgroup$ Apr 4, 2016 at 2:58
  • $\begingroup$ The reason this is obscure is that it's a complicated example to learn the low level FEM stuff. Spend a day or two and go through the low level FEM stuff which has much simpler examples that illustrate how it works and you'll find a good tool to work on all kind of PDE problems. $\endgroup$
    – user21
    Apr 4, 2016 at 3:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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