4
$\begingroup$

How can "MeshElementBlocks" be used with NDSolve to run a parallel kernel evaluation of a Finite Element problem?

Ω = RegionDifference[Rectangle[{0, 0}, {100, 100}], Rectangle[{40, 40}, {60, 60}]];

fea[grid_, goal_, blocks_: 1] := {"FiniteElement",
  "MeshOptions" -> {MaxCellMeasure -> {"Area" -> grid}, 
    AccuracyGoal -> goal, PrecisionGoal -> goal, 
    "MeshElementBlocks" -> blocks}}

sol = Parallelize[NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 0, 
   DirichletCondition[u[x, y] == 100., 
    x == 40 && 40 <= y <= 60 || x == 60 && 40 <= y <= 60 || 
     40 <= x <= 60 && y == 40 || 40 <= x <= 60 && y == 60], 
   u[x, 0] == u[x, 100] == u[0, y] == u[100, y] == 0},
  u, {x, y} ∈ Ω, Method -> fea[1, 8, 4]]]

Reference:

Partial Differential Equation in Parallel

$\endgroup$
  • $\begingroup$ Have a look at this answer. Hope that's helpful. $\endgroup$ – user21 Jul 9 '18 at 13:23
3
$\begingroup$

Here is part of a different approach. This leads in essence to a domain decomposition. One can use the Options "PartialSystemMatricesAssembly" of DiscretizePDE. This example is from the documentation.

Needs["NDSolve`FEM`"]
nr = ToNumericalRegion[Rectangle[{0, 0}, {1, 1/2}]];
vd = NDSolve`VariableData[{"DependentVariables", 
     "Space"} -> {{u}, {x, y}}];
sd = NDSolve`SolutionData[{"Space"} -> {nr}];
cdata = InitializePDECoefficients[vd, sd, 
   "DiffusionCoefficients" -> {{-IdentityMatrix[2]}}];
mdata2 = InitializePDEMethodData[vd, sd, 
   Method -> {"FiniteElement", 
     "MeshOptions" -> {"MeshElementBlocks" -> 5}}];
mdata2["ElementMesh"]

Partially discretize a PDE with blocks numbers 1, 2, and 5 :

dpde1 = 
 DiscretizePDE[cdata, mdata2, sd, 
  "PartialSystemMatricesAssembly" -> {1, 2, 5}]

Extract and visualize the assembled stiffness matrix :

MatrixPlot[dpde1["StiffnessMatrix"]]

enter image description here

Partially discretize a PDE with blocks numbers 3 and 4 :

dpde2 = DiscretizePDE[cdata, mdata2, sd, 
  "PartialSystemMatricesAssembly" -> {3, 4}]

Extract and visualize the assembled stiffness matrix :

MatrixPlot[dpde2["StiffnessMatrix"]]

enter image description here

Discretize a PDE over all mesh elements:

dpde = DiscretizePDE[cdata, mdata2, sd]

Verify that the sum of the partially assembled system matrices is \ equal to the system matrices assembled as a whole :

dpde["StiffnessMatrix"] == 
   dpde1["StiffnessMatrix"] + dpde2["StiffnessMatrix"]
True

Now, in a next step one would solve over the partially assembled system matrices and then construct the solution from that. I don't have code for that. Here is an example from the MATLAB PDE toolbox that does something like this, perhaps that's useful.

$\endgroup$
  • $\begingroup$ @Young, if you do figure this out. I'd be very interested in hearing about that. Good luck. $\endgroup$ – user21 Aug 16 '16 at 15:22
  • $\begingroup$ I've made some progress on this but need more free time to compete ... thank you for the bounty on the other question! $\endgroup$ – Young Aug 23 '16 at 1:09
  • $\begingroup$ @Young, if you figure this one out there will be another one ;-) $\endgroup$ – user21 Aug 23 '16 at 2:06
  • $\begingroup$ @Young did you get a chance to look at this in more detail? $\endgroup$ – user21 Aug 28 '16 at 14:32
2
$\begingroup$

For using "MeshElementBlocks" there is an example for that in the Finite Element Programming tutorial but I am not sure this is what you are looking for. The "MeshElementBlocks" are for lowering the memory consumption and/or for doing something called domain decomposition. While domain decomposition is essentially for solving FEM models in parallel I do not have an example for that.

When you run this

Ω = RegionDifference[Rectangle[{0, 0}, {100, 100}], Rectangle[{40, 40}, {60, 60}]];

sol = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 0, 
   DirichletCondition[u[x, y] == 100., 
    x == 40 && 40 <= y <= 60 || x == 60 && 40 <= y <= 60 || 
     40 <= x <= 60 && y == 40 || 40 <= x <= 60 && y == 60], 
   u[x, 0] == u[x, 100] == u[0, y] == u[100, y] == 0}, 
  u, {x, y} ∈ Ω, 
  Method -> {"FiniteElement", 
    "MeshOptions" -> {MaxCellMeasure -> 0.05}}]

you will see that some parts (the FEM element computation and the LinearSolve) will run in parallel if you look at, for example, top on Linux.

$\endgroup$
  • $\begingroup$ So when using blocks, the MaxMemoryUsed[] should be lower? $\endgroup$ – Young Aug 16 '16 at 13:22
  • $\begingroup$ It's a bit more complicated than that. Have a look at the section about solving memory intensive PDEs in the Finite Element best practice tutorial $\endgroup$ – user21 Aug 16 '16 at 13:30
  • $\begingroup$ Any chance of seeing an example of how to use blocks to parallelize a NDSolve FE problem? $\endgroup$ – Young Aug 16 '16 at 13:45
  • $\begingroup$ @Young not in the immediate future. I'd like to focus on something else before that. Sorry about that. $\endgroup$ – user21 Aug 16 '16 at 14:30
  • $\begingroup$ I understand. Can you just tell me how to break this ElementMesh[{{-0.0365, 0.0365}, {-0.0365, 0.0365}}, {TriangleElement[ "<" 36134 ">"], TriangleElement["<" 36134 ">"], TriangleElement["<" 36134 ">"], TriangleElement["<" 36134 ">"]}] into the 4 separate block regions so I can solve them independently? $\endgroup$ – Young Aug 16 '16 at 14:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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