Partial Differential Equation in Parallel

Question

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?

Answer 1

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.

Answer 2

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

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.