19
$\begingroup$

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?

$\endgroup$
  • 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 Jul 26 '14 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 Jul 27 '14 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 Jan 23 '15 at 7:48
  • $\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 Oct 7 at 6:45
13
$\begingroup$

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.

$\endgroup$
  • $\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 Feb 2 '16 at 21:46
  • $\begingroup$ An example of a larger system that runs FEM code in parallel $\endgroup$ – user21 Feb 3 '16 at 19:02
  • $\begingroup$ @user21 Do I need to invoke Parallelize with ToElementMesh and NDSolve for it to use all available cores? $\endgroup$ – Young Aug 14 '16 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 Aug 15 '16 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 Aug 16 '16 at 5:09
11
$\begingroup$

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

ResourceFunction["FEMAddOnsInstall"][]

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.

Needs["FEMAddOns`"]

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]]]
47120448

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]]]
15087120

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"]]
5.70608*10^-7

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.

$\endgroup$
  • 4
    $\begingroup$ That's good new stuff. I have to give it a try very soon. $\endgroup$ – Henrik Schumacher Jul 9 '18 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$ – ABCDEMMM Mar 16 at 21:37
  • $\begingroup$ @ABCDEMMM, you'll find the paclet on the FEMAddOns release page $\endgroup$ – user21 Mar 18 at 6:21

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.