New in:
11.3​
| Modified in:
| Obsolete in:
| Excised in:
Categorization

Keywords

Details


Domain decomposition

Introduction

The amount of memory needed by
NDSolve
to solve stationary
Partial Differential Equations (PDEs)
with the Finite Element Method (FEM) is proportional to the number of equations and the number of nodes in the underlying mesh. In some cases solving a particular PDE may require more computer memory than is available. Some easy-to-use options for reducing memory consumption are presented in the section
Solving Memory-Intensive PDEs
, but these aren't always enough. This tutorial explains how to additionally use Domain Decomposition Methods (DDMs) in Wolfram Language for solving memory intensive PDEs.
Domain decomposition solvers are useful primarily for two reasons:

To parallelize PDE solving, possibly making use of remote kernels in a computer cluster.

To decrease memory consumption during PDE solving.
Being able to solve larger scale PDEs comes at a cost of a longer solution time as will be explained later in this tutorial.

Reducing memory consumption of PDE solving

Load the package:
In[2]:=
Needs["FEMAddOns`"]
The following is an example of how DecompositionNDSolveValue, which implements a domain decomposition method, uses less memory than NDSolve when solving the Poisson equation over a Menger mesh.
Set up the equations:
In[3]:=
op=-Laplacian[u[x,y],{x,y}]-1;​​bcs=DirichletCondition[u[x,y]0,x0||x1||y0||y1];​​Ω=MengerMesh[4];
Because
NDSolve
does some initialization when run the first time that affects memory usage, an unmetered call to
NDSolve
is performed to give a fair comparison.
Solve the PDE with
NDSolve
to initialize:
In[6]:=
NDSolveValue[{op0,bcs},u,{x,y}∈Ω];
Solve the PDE with
NDSolve
and measure memory consumption:
In[7]:=
mem1=MaxMemoryUsed[sol1=NDSolveValue[{op0,bcs},u,{x,y}∈Ω]]
Out[7]=
47120592
Solve the PDE with domain decomposition and measure memory consumption (this problem takes a long time to solve):
In[8]:=
mem2=MaxMemoryUsed[sol2=DecompositionNDSolveValue[{op0,bcs},u,{x,y}∈Ω]]
Out[8]=
15087696
Compare the memory consumption and solutions:
In[9]:=
{N[mem1/mem2],Norm[sol1["ValuesOnGrid"]-sol2["ValuesOnGrid"]]}
Out[9]=
{3.12311,1.11832×
-7
10
}
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 will be discussed further down.
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 savings in memory come at the price of computational time.
Compare
NDSolve
and DecompositionNDSolveValue with respect to timings:
In[10]:=
t1=AbsoluteTiming[NDSolveValue[{op0,bcs},u,{x,y}∈Ω];][[1]];​​t2=AbsoluteTiming[DecompositionNDSolveValue[{op0,bcs},u,{x,y}∈Ω];][[1]];​​t2/t1
Out[12]=
26.7834

Computing PDEs on a cluster

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.
Solve the PDE with remote kernel objects:
In[13]:=
kernels=LaunchKernels[4];​​sol2=DecompositionNDSolveValue[{op0,bcs},u,{x,y}∈Ω,"Kernels"kernels];​​Norm[sol1["ValuesOnGrid"]-sol2["ValuesOnGrid"]]
Out[15]=
5.70608×
-7
10

How domain decomposition methods work

To solve a PDE, three pieces of information need to be given. A PDE, its boundary conditions, and a region (domain) over which the PDE is to be solved. Consider the following decomposition of the one-dimensional domain
a<x<b
, and how domain decomposition methods can be used to solve a PDE with it.
Two overlapping domains
Ω
1
and
Ω
2
constitute the overall domain
a<x<b
.
Generally speaking there are two types of domain decomposition methods. There are those that make use of overlapping domains, called Schwarz methods, and those that make use of non-overlapping domains. DecompositionNDSolveValue implements Schwarz methods. Schwarz methods work by solving the PDE not over the entire domain Ω but by solving the same PDE over overlapping subsets
Ω
1
and
Ω
2
. In other words, instead of finding the solution u(x) in the entire domain the Schwarz methods try to find two solutions
u
1
(x) for x ∈ (a,
γ
1
)
and
u
2
(x) for x ∈
(
γ
2
, b). The solutions
u
1
and
u
2
are then stitched together to make the solution
u
. To be able to solve the PDE on the subdomains artificial Dirichlet boundary conditions are introduced at
γ
1
and
γ
2
.
u
1
and
u
2
are coupled by the artificial Dirichlet conditions
u
1
(
γ
1
)=
u
2
(
γ
1
)
and
u
2
(
γ
2
)=
u
1
(
γ
2
)
. The boundary conditions are satisfied through an iterative scheme which, supposing that the solution is zero and that
(0)
u
2
is given, might look like this:
Dirichlet boundary conditions
k
u
g
at
Γ
g
for iteration number
k
.
The figure gives an intuition for why, at least in this case, the convergence rate of Schwarz methods is dependent on the size of the overlap. More generally, one might have subdomains in a higher dimension. For example:
Here is a more formal description of the solution process. Consider the PDE, which can also have more general boundary conditions,

Lu=f
onΩ
u=g
on∂Ω
​where L is an elliptic differential operator, u is the dependent variable, Ω is the union of two subdomains
Ω
1
and
Ω
2
, and ∂Ω is its boundary. A Schwarz method proceeds by first solving the PDE over
Ω
1
using boundary conditions from the physical boundary
∂
Ω
1
\
Γ
1
and an artificial Dirichlet boundary on
Γ
1
, with
(0)
u
2
given.
(k+1)
Lu
1
=f
in
Ω
1
(k+1)
u
1
=
(k)
u
2
on
Γ
1
(k+1)
u
1
=0
on∂
Ω
1
\
Γ
1
We specify a model PDE that will be used in this section.
Set up a model PDE:

The multiplicative and additive Schwarz methods

DecompositionNDSolveValue implements the multiplicative Schwarz method and the additive Schwarz method for sequential evaluation. They both have robust convergence and the same decrease in the amount of memory, but multiplicative Schwarz has a better convergence rate and is therefore faster when used sequentially (i.e. no kernels have been provided). The additive Schwarz method is inherently parallel, and it is currently the only algorithm available when using the Kernels option to solve problems in parallel.
Solve the equations with additive Schwarz and multiplicative Schwarz:
Make a table comparing the solutions:
The example shows that multiplicative Schwarz and additive Schwarz consume the same amount of memory, but that multiplicative Schwarz is faster.

Adjusting the tolerance of the convergence criterion

DecompositionNDSolveValue stops iteration when the norm falls below the given tolerance:

Providing an approximate solution to speed up convergence

The computations can be sped up by making the Schwarz algorithm converge faster. There are two way to achieve a faster convergence:
The better of an approximation of the actual solution that can be provided, the fewer iterations will have to be done to satisfy the tolerance criterion. The following example shows how to provide an initial guess.
Create a mesh and a coarse version of that mesh:
Find the solution using the coarse mesh:
Count the number of iteration until convergence:
Count the number iterations it takes for the solution to converge when starting from the solution obtained using the coarse mesh:
The example shows that providing an estimate of the solution can decrease the number of iterations required for the solution to converge and therefore decrease the amount of time it takes to solve the problem. The code makes use of the function SetMaxCellMeasure, which takes a mesh and its corresponding symbolic region and creates a new, coarser mesh. It does this by pruning nodes in the body of the mesh while leaving the boundary nodes as they are, thus ensuring that the new mesh covers exactly the same space as the original mesh.

Solving problems in parallel using local or remote kernels

Solve a problem in parallel using two local kernels:
The additive Schwarz method is always used when kernels are provided, because it is very well suited as a parallel domain decomposition method. A parallel version of multiplicative Schwarz has not been implemented in this version. DecompositionNDSolveValue will divide the domain into as many subdomains as there are kernels. If subdomains are specified then the number of kernels should match the number of subdomains. DecompositionNDSolveValue cannot be used with fewer kernels than subdomains, if there are more kernels than subdomains then only the required number of kernels will be used. The subdomains are setup on one kernel each. Each kernel is responsible for solving the PDE over its subdomain once every iteration. The first kernel in the list is matched with the first subdomain in the list and so on. This could be used to place a particularly large subdomain on a computer with more memory. Using a different number of kernels than subdomains will give an error message.

Monitoring the convergence of the solution

The option EvaluationMonitor can be used to evaluate a function after each iteration. The function is given as its first argument the current estimate of the solution and as its second argument the previous estimate of the solution. This can be used to, for example, visualize the progress of the solution, or to compute the norm of the difference between the two arguments to see how far away the norm is from dropping below the specified tolerance.
Monitor convergence using EvaluationMonitor:
The example shows how the norm decreases iteration by iteration until it drops below the tolerance threshold.

Passing subdomain options to DecompositionNDSolveValue

More information can be found in the documentation for ElementMeshDecomposition below.

DecompositionNDSolveValue options

The following table summarizes the options available for DecompositionNDSolveValue.
The currently available domain decomposition methods work for stationary, elliptic partial differential equations in one, two, and three dimensions. As a rule of thumb, any steady partial differential equation that the FEM package can solve, can also be solved with the Schwarz algorithms. Currently, only additive Schwarz is supported for solving PDEs in parallel.
The following example demonstrates that domain decomposition methods can handle equations with convection.
Set up the equation:
Solve convection-diffusion equation using NDSolve:
Solve convection-diffusion equation using DecompositionNDSolveValue:
Compare solutions:
The class of supported PDEs is not limited to systems of equations with one dependent variable as in the previous examples. This makes it possible to use the Schwarz algorithms for solving problems in, for example, structural mechanics. In the following problem, a force is pulling at the edge of a metal plate:
Set up the equations:
Solve the coupled PDEs using DecompositionNDSolveValue:
The convergence rate of one-level Schwarz algorithms, which is the kind of Schwarz algorithms currently implemented, is known to deteriorate with an increasing number of subdomains. The reason is that subdomains only exchange information with their nearest neighbors.
Problem setup:
Solve the problem using NDSolve:
Solve the problem using 2 subdomains, overlap 4, and 50 iterations:
Solve the problem using 5 subdomains, overlap 4, and 50 iterations:
Solve the problem using 10 subdomains, overlap 4, and 50 iterations:
The previous examples show that 50 iterations yield a much better solution when there are fewer subdomains. This is because the convergence rate decreases with an increasing number of subdomains. The convergence rate can be improved in two ways; by increasing the amount of overlap, or by providing a better initial guess.
Solve the problem using 10 subdomains, overlap 20, and 50 iterations:
Solve the problem using 10 subdomains, overlap 4, and 50 iterations, with an initial guess:

Monitoring convergence

Specifying a method for LinearSolve

Set up equations:
Solve the PDE using default solver:
Solve the PDE using "Pardiso":
The smaller memory consumption means that it almost always makes sense to use Pardiso when trying to decrease memory consumption using domain decomposition methods.

Subdomains

The preceding sections have talked about how to solve PDEs using domain decomposition methods. The details of how domains are decomposed into parts, subdomains, has not been discussed as it has been handled by DecompositionNDSolveValue. Sometimes it can be useful to decompose the element mesh separately without using DecompositionNDSolveValue. For example, to check what the subdomains look like, to inspect their properties, or even to implement custom domain decomposition methods. ElementMeshDecomposition is a function that takes an element mesh and returns a decomposition of that element mesh into subdomains. The subdomains can be queried to find information about the subdomain mesh properties and the relationships between the subdomains in terms of for example shared incidents, artificial boundary incidents, and physical boundary incidents. Subdomains returned by ElementMeshDecomposition can be given to DecompositionNDSolveValue in place of a region.

Creating subdomains using ElementMeshDecomposition

Any element mesh can be divided into subdomains using the ElementMeshDecomposition function.
Decompose a rectangle into two subdomains and visualize them:
ElementMeshDecomposition has two options: "Subdomains" and "Overlap". Overlap is the number of layers one subdomain reaches into the neighboring subdomain. In the following example, four layers of elements overlap, which is two layers in each direction from the boundary in between the subdomains in the preceding example. That is, two layers in each direction from the boundary in between the corresponding nonoverlapping subdomains.
Specify two layers of elements of overlap:
"Subdomains" specifies how many subdomains the region should be decomposed into. If an integer is given then ElementMeshDecomposition will attempt to split the element mesh into that many subdomains. If a list of fractions is given then it will try to decompose the element mesh into as many subdomains as there are elements in the list, each with approximately the specified fraction of the total area.
Decompose a Menger mesh into four subdomains with six layers of overlap:
The figure shows four subdomains. The overlapping regions are given distinct colors as well.
Use a list of fractions to divide an element mesh into subdomains of different sizes:
In the preceding example, the blue subdomain takes up roughly 60 % of the area, the yellow subdomain takes up roughly 30 % and the green subdomain takes up roughly 10 %. Subdomains of unequal sizes can be used to distribute the computational load unevenly across kernels if, for example, one kernel has access to more memory than another.

Set overlap as a fraction of the subdomain element graph diameter

Consider an undirected graph representing the connections between elements in a subdomain element mesh. It turns out that the convergence rate is related to the overlap as a percentage of the graph diameter. It is often easier to specify the amount of overlap as a fraction of the graph diameter than as a number of layers because of this. This can be done.
Set the overlap to 10 % of the graph diameter
Use the same overlap specification for a mesh with more subdomains
Use the same overlap specification for a mesh with more mesh elements
The examples show that specifying the overlap as a fraction of the subdomain graph diameter is a good way of ensuring that the overlap will adapt if the number of elements in the subdomains changes.

Specifying overlap for subdomains individually

The subdomains used so far have all extended themselves into their neighbors by approximately same number of layers. It is also possible specify the number of layers for each subdomain individually.
Set overlap to zero for the subdomains in the middle, 20 for the ones at the edge:

One-dimensional and three-dimensional examples

It is not just two-dimensional domains which can be decomposed by ElementMeshDecomposition, but also one-dimensional and three-dimensional domains.
Decomposing a one-dimensional domain into subdomains:
Decomposing a three-dimensional domain into subdomains:

Subdomain wireframe visualization

"Wireframe" accepts the following options, which are the same for subdomains as it is for ElementMesh:
Any options accepted by Graphics will also be applied to the wireframe.
Using the ImageSize option of Graphics to set the size of the resulting graphics:
Specify an edge color for all elements regardless of subdomain:
Lists can be used to specify different options for different subdomains. The list will repeat itself cyclically if the number of stylings in the list is shorter than the number of subdomains.
Style specifications are reused cyclically:

Subdomain data structures

The structure of the object returned by ElementMeshDecomposition is SubdomainGroup[mesh, {Subdomain[...], Subdomain[...], ...}] where mesh corresponds to the global domain. The Subdomain and SubdomainGroup objects have properties. These properties can be used to implement custom domain decomposition methods. They are listed below.
The corresponding incidents in the global mesh are given by the following properties. This is the same as applying the mapping given by LocalToGlobal to the local incidents.

Global versus local incidents

Because each subdomain has fewer incidents than the full mesh the subdomain incidents do not match up with the list of coordinates in the full mesh. This means that two sets of incidents are needed to, say, map the solution vector of a subdomain onto the solution vector of the full mesh. Each type of incident is therefore provided both as local incidents and global incidents. A map is also provided which maps local incidents into global incidents.
Create two subdomains:
Find the coordinates of the artificial boundary using a global mesh:
Find the coordinates of the artificial boundary using a subdomain mesh:
The difference between the first way of getting the boundary incidents and the second is that the incidents are retrieved using the GlobalArtificialBoundaryIncidents property in the first example, and ArtificialBoundaryIncidents in the second example. As always, incidents correspond to coordinates, but in this case they correspond to coordinates in different element mesh objects. The list GlobalArtificialBoundaryIncidents corresponds to coordinates in the full mesh and ArtificialBoundaryIncidents corresponds to elements in the local mesh. In the next example it is shown how incidents retrieved using ArtificialBoundaryIncidents can be converted into incidents that can be used with the coordinates of the full mesh.
How to convert local incidents into global incidents:
Problem setup:
Extract subdomain mesh properties:
Initialize local and global solution vectors:
Extract incident data from subdomains:
Implementation of the Alternating Schwarz method:
Compare solution with NDSolve:
Visualize the solution using ElementMeshInterpolation and Plot3D:

Possible issues

DecompositionNDSolveValue may not be able to solve some PDEs with discontinuous coefficients.
Set up equations:
Solve using DecompositionNDSolveValue:
Compare solutions:

References

Related Tutorials