is there any native way to implement multi-core parallel solving of PDE in Wolfram Mathematica?

WM 10 now supports Finite Elements Method, but it is actually useless without parallelization. Usually PDE-systems are heavy and no one solves it using just one core. If there is any sort of domain decomposition or some other way parallelize the calculation?

  • 4
    $\begingroup$ Calling it "useless" is unfair. You get great convenience from Mathematica at the price of performance. If you solve the same in C or Fortran, you'll gain much more performance than parallelization would give you in Mathematica. Use the right tool for the problem. $\endgroup$
    – Szabolcs
    Commented Jul 26, 2014 at 22:55
  • 1
    $\begingroup$ Yes, yes, this is usual practice to use C/Fortran codes or sort of COMSOL, ANSYS, OpenFOAM et ct. Still, they're not posess Mathematica flexibility and I think I could find some applications for Mathematica too. $\endgroup$
    – user18790
    Commented Jul 27, 2014 at 13:10
  • $\begingroup$ Has anybody tried using Mathematica to work out the equations symbolically, and then exporting the problem description to some other solver? $\endgroup$
    – Åsmund Hj
    Commented Jan 23, 2015 at 7:48
  • 1
    $\begingroup$ Side note: At least since v11.3, TensorProductGrid method of NDSolve automatically solves PDE in parallel, which doesn't seem to be the case in v9.0.1 $\endgroup$
    – xzczd
    Commented Oct 7, 2019 at 6:45

2 Answers 2


The Finite Element solver in Mathematica does run in parallel, both element computation and the linear solve process are spread over the CPU cores available.

Additionally, the option "MeshElementBlocks" for ToElementMesh splits the mesh elements in blocks which could be used for a domain decomposition.

To get a more detailed answer you'd need to clarify you question.

  • $\begingroup$ one thing that can be addded is the fact that FE spreads over all cores available on a single machine. I can agree that, because of this single machine restriction, some applications are beyond reanosability (I myself frequently require 300 cores to remain within acceptable times). But even if NDSolve could spread over 300 cores, I don't think that the current practical use we can make of it, really requires more than, lets say, 24 cores. If it does, most likely, it is not the best tool for the job, because of other aspects besides the cores limitation... (but I may be wrong...) $\endgroup$
    – P. Fonseca
    Commented Feb 2, 2016 at 21:46
  • $\begingroup$ An example of a larger system that runs FEM code in parallel $\endgroup$
    – user21
    Commented Feb 3, 2016 at 19:02
  • $\begingroup$ @user21 Do I need to invoke Parallelize with ToElementMesh and NDSolve for it to use all available cores? $\endgroup$
    – Young
    Commented Aug 14, 2016 at 22:18
  • $\begingroup$ @user21 I can't seem to get Mathematica to use more than one core while evaluating NDSolve using the "FiniteElement" method. Hints? $\endgroup$
    – Young
    Commented Aug 15, 2016 at 0:28
  • $\begingroup$ @user21 Can you point me to an example of using "MeshElementBlocks" with NDSolve to break a problem into regions and solve each region in parallel? $\endgroup$
    – Young
    Commented Aug 16, 2016 at 5:09

The open source Wolfram Research FEMAddOns package has a domain decomposition solver called DecompositionNDSolveValue. You can install the paclet with evaluating:


That domain decomposition solver can solve stationary PDEs on a cluster of heterogeneous computers/CPUs. The package has an extensive tutorial about domain decomposition, it's pros and cons and information on how to tweak options to get a good performance. I'll give a few examples below but you should look at the tutorial. You'd probably needed a Version 11.3 for these functions to work. For installation instructions please see link above.

To use the package it needs to be loaded.


We set up a model problem:

op = -Laplacian[u[x, y], {x, y}] - 1;
bcs = DirichletCondition[u[x, y] == 0, 
   x == 0 || x == 1 || y == 0 || y == 1];
\[CapitalOmega] = MengerMesh[4];

To do some comparison we do a dry run of solving the PDE with NDSolve. This compiles and loads some code that we do not want in the comparison.

NDSolveValue[{op == 0, bcs}, u, {x, y} \[Element] \[CapitalOmega]];

Solve the PDE with NDSolve and measure memory consumption:

mem1 = MaxMemoryUsed[
  sol1 = NDSolveValue[{op == 0, bcs}, 
    u, {x, y} \[Element] \[CapitalOmega]]]

Solve the PDE with domain decomposition and measure memory consumption. This problem takes a few moments to solve:

mem2 = MaxMemoryUsed[
  sol2 = DecompositionNDSolveValue[{op == 0, bcs}, 
    u, {x, y} \[Element] \[CapitalOmega]]]

Compare the memory consumption and solutions:

{N[mem1/mem2], Norm[sol1["ValuesOnGrid"] - sol2["ValuesOnGrid"]]}
{3.12322, 1.11832*10^-7}

The memory consumption of DecompositionNDSolveValue is smaller compared to NDSolve. In large scale computations this can be the difference that makes it possible to solve a PDE at all. The amount of memory saved can be influenced through the setting of options which are discussed in the tutorial.

Domain decomposition methods work by dividing the computational region into subdomains. The PDE under consideration is then solved on these subdomains and the solution on the original region is stitched together from the solutions of the subdomains. The memory consumption when computing the solution will never be larger than what it is when solving the problem over the largest subdomain.

The domain decomposition method is also advantageous because the PDE can be solved simultaneously over the subdomains; not necessarily on the same workstation. For example, a computer cluster can be used by assigning one subdomain to each server.

kernels = LaunchKernels[4];
sol2 = DecompositionNDSolveValue[{op == 0, bcs}, 
   u, {x, y} \[Element] \[CapitalOmega], "Kernels" -> kernels];
Norm[sol1["ValuesOnGrid"] - sol2["ValuesOnGrid"]]

And because we are missing a pretty picture in this answer, here is how a mesh is partitioned for NDSolve to solve the sub problems:

mesh = ToElementMesh[MengerMesh[2]];
subdomains = 
  ElementMeshDecomposition[mesh, "Subdomains" -> 4, "Overlap" -> 6];
subdomains["Wireframe"[ImageSize -> 300]]

enter image description here

The FEMAddOns package is an open source project and contributions are very much welcome. If the finite element method or numerics in general is your field of interest then this is a change to contribute FEM and related code. It may also help Wolfram Research to see if and how customers would like to interact with the company on a very technical level.

  • 4
    $\begingroup$ That's good new stuff. I have to give it a try very soon. $\endgroup$ Commented Jul 9, 2018 at 13:49
  • $\begingroup$ @user21 , may I ask a question, how can we install FEMAddOns in windows? I give a try, however, i have no idea how to install FEMAddOns.paclet, I never find such FEMAddOns.paclet folder... $\endgroup$
    Commented Mar 16, 2019 at 21:37
  • $\begingroup$ @ABCDEMMM, you'll find the paclet on the FEMAddOns release page $\endgroup$
    – user21
    Commented Mar 18, 2019 at 6:21

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.