Current state of the art supercomputers perform several petaflops and thus enable large-scale simulations of physical phenomena. The enormous computational power is particularly valuable when large computational domains must be resolved with a fine mesh. For minimizing the resource consumption in, e.g., compute time, processor numbers, memory, or energy, we must proceed beyond the demonstration of asymptotic optimality for the algorithms and beyond establishing the scalability of the software. Modern fast solvers must additionally be implemented with data structures that permit an efficient execution on multicore architectures exploiting their instruction level parallelism, and they must be also adapted to the complex hierarchical memory structure.This necessitates a meticulous efficiency-aware software design and an advanced performance aware implementation technology.

MLMC methods are characterized by three algorithmic levels that are potential candidates for parallel execution. As in standard Monte Carlo methods, a sequence of classical deterministic problems (samples) are solved that can be computed in parallel. Additionally one can distinguish between parallelism within an MLMC level and parallelism across MLMC levels. The third algorithmic parallelization level is the solver for the deterministic PDE problem. Indeed, the total number of samples on finer MLMC levels is typically moderate, so that the first two levels of parallelism will not suffice to exploit all available processors.We analyzed how these different levels of parallelism can be combined and how efficient parallel MLMC strategies can be designed.

In the case of a fault, the redundant interface values can be easily recovered while the lost inner unknowns are re-computed approximately with recovery algorithms using multigrid cycles for solving a local Dirichlet problem. Different strategies are compared and evaluated with respect to performance, computational cost, and speed up. Especially effective are asynchronous strategies combining global solves with accelerated local recovery.

The computations in both the faulty and healthy subdomains must be coordinated in a sensitive way, in particular, both under- and over-solving must be avoided. Both of these waste computational resources and will therefore increase the overall time-to-solution. To control the local recovery and guarantee an optimal re-coupling, we introduce a stopping criterion based on a mathematical error estimator. It involves hierarchically weighted sums of residuals within the context of uniformly refined meshes and is well-suited in the context of parallel high-performance computing. The re-coupling process is steered by local contributions of the error estimator before the fault.

The input mesh in our first experiment is the unit cube discretized by 24576 tetrahedrons, which results after uniform refinement in a linear system with 1.3·10⁸ DOFs. To load-balance the problem, we run these simulations with 48 MPI-processes and assign to each process 512(= 83) input mesh tetrahedrons. We use in the re-coupling criterion the bound σLRB and κ=1.0 for the two superman acceleration factors η_s = 2 and η _s = 4. In he fowollowing Figure 8, we present the estimated algebraic error within a faulty solution process. The plots left and right refer to ηs=2 and ηs=4, respectively. In order to satisfy the re-coupling criterion, 6 local multigrid iterations are necessary in the faulty subdomain which can be performed simultaneously while 3.7 and 2.0 iterations are executed in the healthy domain, respectively, depending on the superman acceleration ηs. In the faulty domain, the error estimator is applied in addition to the V-cycles. Asynchronously, the computations in the healthy domain are executed. The recovery process is highlighted in light green in the figures. Note, that fractional numbers imply that the healthy domain receives the stopping signal for re-coupling when the multigrid iteration in the healthy domain has not completed the most recent cycle. So for example, when the healthy domain has executed the downward branch of a symmetric V-cycle, we will have executed one half of a V-cycle and denote this in the following tables and graphics by 0.5 iterations. Summarizing our experiments, we find that the recovery with η_s=2 results in a delay of only up to three iterations and with η_s=4 in up a delay of two iterations. Thus, the adaptive re-coupling strategy reduces significantly the costs in comparison to applying no recovery strategy.

We performed strong scaling experiments of MLMC using homogeneous bulk synchronous scheduling with up to 131072 of processes. Increasing the number of processes while keeping the number of samples per level fixed results in an increased load imbalance and the total parallel efficiency decreases. Nevertheless, even with 32768 processes, we obtain a parallel efficiency of over 60%. For 131072 processes, the efficiency drops below 40%, since we have reached the limits of the scalability window. This is inevitable in any strong scaling scenario. Overall, the compute time for the MLMC estimator can be reduced from 616 to 22 seconds, resulting in an excellent parallel efficiency of the MLMC implementation.