{'fraction_non_alphanumeric': 0.027969589174862246, 'fraction_numerical': 0.008125828276487411, 'mean_word_length': 4.698529411764706, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical simulation of skin transport using Parareal

In-silico investigations of skin permeation are an important but also computationally demanding problem. To resolve all involved scales in full detail will probably not only require exascale computing capacities but also suitable parallel algorithms. This article investigates the applicability of the time-parallel Parareal algorithm to a brick and mortar setup, a precursory problem to skin permeation. A C++ library implementing Parareal is combined with the UG4 simulation framework, which provides the spatial discretization and parallelization. The combination's performance is studied with respect to convergence and speedup. While limited speedup from the time parallelization is shown to be possible, load balancing is identified as an important issue in the application of Parareal with implicit integrators to complex PDEs.


Introduction
Permeation of chemical substances through human skin is an interesting and important process e.g. for the development of cosmetics or drugs. In-vitro studies with humans constitute the "gold standard" but they are expensive and limited by ethical and practical concerns. Here, in-silico studies are a viable alternative. They have been successfully used in the past (cf. the reviews in [1,20]) and can be expected to expand further in the future. They allow for hypothesis testing and may lead to experiments through which effects not known today could be discovered.
Yet, numerical simulations in this field are demanding in terms of computational resources. The problem covers vastly different physical scales and, in case a complex full-fledged model is used, massive computational parallelism needs to be exploited to reach reasonable times-to-solutions. Therefore, many interesting aspects such as substructures of lipid bilayers or networks of keratin filaments [32,33] have not yet been investigated in three spatial dimensions (3D) using numerical simulations. Finally, modern imaging techniques make resolving a spectral range of a few hundred nanometers possible, which results in "big data" for analyses. Understanding the functional mechanism of the skin is thus a potential candidate from the life sciences for applying exascale computing.
Considering the technology trend towards more and more parallelism, the application of new parallel methods to the problem investigated here becomes relevant. Promising candidates for such methods are parallelin-time integration methods that can add a direction of concurrency in addition to spatial parallelization, e.g., as done here, through parallel multi-grid. In recent years, time-parallel methods have matured from a mainly mathematical concept to an approach with demonstrated efficiency in massively parallel computations [29,26], and they have been listed as a direction for mathematical research that has the potential to help reaching exascale computing [9].
One of the most widely investigated parallel-in-time methods is Parareal, introduced in 2001 by Lions, Maday and Turinici [16]. It has been used for benchmark problems motivated by applications from fields as diverse as plasma physics [28], computational fluid dynamics [8,21] or quantum chemistry [7]. Improvements with respect to implementation are considered e.g. in [5,10]. Parareal's most appreciated aspect is probably that it is non-intrusive and rather easy to integrate into existing codes. Its drawback on the other hand is a quite severe bound on achievable parallel efficiency. However, because several other "across-the-step" timeparallel methods share similar features with Parareal (e.g. PITA [13], MGRIT [12] or PFASST [11]), studying Parareal's performance often already gives important insights.
Theoretical estimates for stability and convergence of Parareal for linear diffusive problems with constant coefficients can be found in [14]. Theory for diffusive problems with constant coefficients can also be found in [6]. For 2D diffusion with space-time dependent coefficients, numerical experiments showed only a marginal reduction in convergence speed [27]. In [4], the small impact of a time dependent viscosity on Parareal's convergence for an advection-diffusion problem is demonstrated. However, performance for a 3D diffusive problem on a complex geometry with anisotropies has not yet been studied.
In preparation for the eventual application of Parareal to skin permeation, this article provides an investigation of Parareal's performance for a 3D brick and mortar problem. From our point of view this model serves as an excellent benchmark because it features challenges resulting from a complex anisotropic geometry and from jumping coefficients, which need to be resolved adaptively over long time intervals. Although locally the mathematical formulation of the brick and mortar problem is clear, the overall picture is highly complex and linked to a real world application requiring a sound simulation infrastructure in terms of numerical methods and software. Here, we employ the simulation framework ug4 [31] for the spatial discretization and linear solvers, for which excellent parallel scaling has been demonstrated [22]. We parallelize ug4's serial temporal solvers through the C++ Parareal library Lib4PrM, which is integrated as a plugin.
The present article establishes the principle applicability of Parareal to the 3D brick and mortar prob-lem. In doing so, it identifies a set of relevant issues to be tackled in order to develop an improved parallelin-time integrator that can deliver reasonable efficiency for the skin transport problem. In particular, load balancing in time is identified as a critical issue when combining implicit methods for a complex PDE with Parareal. Because the number of iterations of the spatial solver varies in time, balancing temporal subintervals in Parareal simply by the number of time steps induces load imbalances. Among other things, we show that this has a noticeable effect on parallel efficiency.

Problem and methods
As a test case we study a simplified version of the 3D brick and mortar problem introduced earlier in [23,19]. The benchmark is defined on a biphasic domain Ω ⊂ R 3 that consists of two disjoint subsets Ω cor , Ω lip ⊂ R 3 representing the so called corneocyte and lipid phase respectively. To be more specific, Ω is the interior of the union Ω cor ∪ Ω lip of closures. On this geometry we solve a diffusion equation with space dependent diffusion coefficients. The evolution of the drug concentration c p (x, t), p ∈ {cor,lib}, is modeled by the equation with x ∈ Ω p , t ∈ [0, T ], and a phase-dependent diffusion coefficient In terms of boundary conditions we obtain an interior phase boundary Γ = Ω cor ∩Ω lip and an exterior boundary ∂Ω. For the exterior boundary we consider a mix of Dirichlet and homogeneous Neumann conditions, i.e.
where n denotes the outward pointing normal vector on ∂Ω. At the phase boundaries Γ , the flux must be continuous, i.e.
Along Γ , concentrations can be discontinuous. However, they are often assumed to be linked by a Nernst's equilibrium K cor/lip c lip = c cor . When K cor/lip is constant, the model may be reformulated, e.g., in terms of a continuous concentration [23,19]. Hence this work employs the simplifying assumption c lip = c cor . Note that more complex situations with locally fluctuating (a) with cornoecytes (yellow) (b) corneocytes removed Fig. 1 A sketch of the 3D brick and mortar problem is shown. The geometry has ten layers of corneocytes (Ω cor , yellow) that are embedded in a matrix of lipid bilayers (Ω lip , blue).
The geometry used in this article is a 3D brick and mortar configuration as depicted in figure 1. It features ten layers of corneocytes (yellow) that are embedded in a matrix of lipid bilayers (blue). It is a simplified version of more elaborate tetrakaidekahedral models with hybrid grids presented previously [19]; the brick and mortar model only consist of hexahedra with a reduced level of anisotropy. However, since jumping coefficients are present, it features some of the issues one encounters also in more complex situations. Figure 2 shows the solution at three different times. To allow for a view into the interior of the domain, a cuboid representing a quarter of the total domain has been removed in the representation. Initially, the solution is zero on the whole domain with a Dirichlet boundary condition of c p (x, t) = 1 on the top side (not shown).
Then, the tracer starts to diffuse downwards in the lipid channels and much slower in comparison through the corneocytes. In the first subfigure, diffusion has just started and filled the upper half of the domain but the concentration is still essentially zero in the lower half. In the last subfigure, approximately at a quarter of the total simulation time, the tracer has diffused down through the whole domain. Concentrations continue to change with smaller changes from step to step, eventually approaching a steady state. A prospective 3D model with a fully resolved lipid-bilayer substructure, as suggested for two dimensions in [32,33], will feature a similar effect on an even smaller scale in time and space.

Parareal
Let C and F be a "coarse" and "fine" time integration method 1 with time step size ∆t and δt ∆t, respectively. For the sake of simplicity assume that all subintervals have the same length and that a constant number of both δt and ∆t steps covers a subinterval exactly. Then, instead of serially integrating across [0, T ] with F, Parareal uses the iteration with k as iteration index. For the first subinterval, i.e. for [0, t 1 ], set for all k, where c 0 is the given initial value. Note that as the iteration converges and c k+1 n+1 − c k n+1 approaches zero for all n = 0, . . . , N t −1, the Parareal solution c k+1 n+1 converges to the serial fine solution F(c n ).
Formulation (6) introduces concurrency because as soon as the iterates c k n are known, the computationally expensive evaluation of the fine method can be done in parallel over all subintervals. That is, the time spent using the fine method in parallel equals the runtime of the fine method across a single subinterval instead of the full interval [0, T ]. However, multiple iterations are typically required and the propagation of corrections by the coarse method remains serial in time. Speedup   therefore depends on finding a cheap enough coarse integrator that still leads to convergence in a small number of iterations.

Speedup from Parareal
Denote by N c the number of coarse time steps per subinterval, by N f the number of fine time steps per subinterval and by N t the number of subintervals, which is assumed to be equal to the number of processors in the temporal direction. Further, denote by N i the number of Parareal iterations and by τ c and τ f the computational runtime for a coarse or fine time step, re-spectively. If one assumes that every time step takes the same amount of time, speedup from Parareal can be modeled by See e.g. [17] or [4] for a more detailed discussion of this bound. For the results presented below, however, this bound turns out to be very inaccurate, because, as the solution approaches a steady-state, the spatial solver requires fewer and fewer iterations per time step, making later time steps cheaper. We thus devise a somewhat modified model, dropping the assumption that each time step takes the same amount of time to compute. Here, to simplify the notation, we omit the index range for sums, maxima and minima; it is always implied to be n = 0, . . . , N t −1. Now, denote by γ c n and γ f n the cost of running the coarse and fine method across the subinterval [t n , t n+1 ]. Then, a serial run of the fine method amounts to the duration while a Parareal run with N i iterations costs as the runtime for the parallel fine solve will be dominated by the subinterval with the longest simulation time. Also, using proper pipelining, the parallel runtime of the serial coarse correction step will be governed by the most expensive subinterval for C. The resulting estimate for the speedup of Parareal is then . (12) This is a slight generalization of (8) in the sense that if γ c n = N c τ c and γ f n = N f τ f are constant for all subintervals, we get and for which (12) simplifies to (8). According to (12), in the case of imbalances in the distribution of computational load across subintervals, possible speedup is limited by the subinterval with the longest runtime for both the fine and coarse method. The optimal configuration therefore corresponds to equal computing times for all subintervals. For explicit schemes, where the cost per time step is more or less constant, this is relatively easy to achieve by making sure every subinterval handles the same number of coarse and fine steps, resulting in the simple speedup model (8). For implicit methods, however, cost per time step is mainly determined by the cost of the spatial solver (i.e. the number of iterations), which in turn depends on the a priori unknown dynamics of the solution. As demonstrated below, load balancing Parareal with implicit methods based on the number of time steps alone can lead to a significant loss in efficiency. Unfortunately, it seems that devising a proper load balancing in time for the implicit case is not straightforward and, to the best of our knowledge, has not yet been addressed in the literature. The easiest approach may be to use the information from the initial coarse run in Parareal to determine the size of the subintervals but this requires a non-negligible amount of implementation, might inhibit proper pipelining when requiring synchronization at some point and is thus left for future work.

Implementation
For Parareal we use the C++ library Lib4PrM that uses MPI for the necessary communication of volume data in time. A straightforward approach to implementing (6) is sketched as pseudo code e.g. in [4]. Here, however, we use a somewhat more elaborate implementation that is based on the following observation. After the first iteration on the first subinterval [0, t 1 ], the coarse terms in the Parareal iteration (6) cancel out, resulting in for k ≥ 1. That is, after one iteration, the first subinterval is guaranteed to have converged and the processors responsible for the first subinterval can "retire". After the second iteration, by the same argument, this will then be true for the processors handling the second subinterval and so on. After k iterations, the timeparallel fine method is guaranteed to have converged on the first k subintervals and all processors with an MPI rank (in time) smaller or equal to k could in principle be used otherwise. Put differently, Parareal converges always at least at a rate of one subinterval per iteration and when k = n the Parareal method is guaranteed to have converged at t k ≤ t n . While leaving processors idle according to this implementation does not affect runtime negatively, it has the potential to reduce the energy cost of a simulation, particularly in combination with "dynamic voltage and frequency scaling". This will be studied in a future work [15]. Also, if not enough processors are available to cover the whole interval [0, T ] by subintervals of a given size, converged processors could pick up subintervals at the end in a caterpillar-like way. Such even more involved implementations are left for future studies, though. Finally, in production runs one could also use some tolerance e.g. for the updates between iterations or a proper residual to decide when the solution at the end of a subinterval is converged [25].

Spatio-temporal discretization and solvers
Discretization in space and time is provided by the software package ug4 [31]. We employ a plain vanilla vertex centered finite volume scheme in space that is combined with an implicit Euler scheme in time. For each time step this gives rise to a large linear system of equations, where the number of degrees of freedom corresponds to the number of vertices of the grid. The solver is a multigrid method with three steps of damped (ω = 0.6) Jacobi relaxation used for pre-and post-smoothing. The coarse level solver is an ILUT-iteration with a threshold of 10 −6 . On the coarsest level at most two steps are performed. More details are provided in [30].

Numerical results
We report results from solving the brick and mortar problem described in §2 with the simulation framework also described above. For both C and F we use an implicit Euler method with the time step size ∆t for C being significantly larger than the time step size δt for F. All runs are performed on the Cray XE6 Monte Rosa supercomputer at the Swiss National Supercomputing Centre CSCS in Lugano, Switzerland. This supercomputer is equipped with 1'496 compute nodes, each of which consists of two 16-core AMD Opteron CPUs, making for a total of 47'872 compute cores. 2 Its peak performance is 402 TFlops, placing it at position 171 in the Top500 November, 2014 list. 3 As compiler we used version 4.7.0 of the the GNU compiler collection 4 and, in the following, report runtimes of simulations without I/O operations.

Convergence of Parareal
Generally speaking, convergence of Parareal is affected by a number of parameters: The time step sizes and methods used for C and F, the number of concurrently computed subintervals, the discretization used for the spatial derivatives, and the dynamics of the problem to be solved. To measure convergence, below the relative defect d k n between Parareal after k iterations and the fine solution run in serial is reported, i.e.
In order to avoid distortions through I/O times, only the final solution values are written out and the defect d k Nt−1 is reported except for in section §3.3. There, the defect is analyzed not only as a function of the number of iterations but also time. Figure 3 shows the defect d k Nt−1 versus the number of iterations k. In addition, the estimated discretization error of both the coarse and fine method resulting from a comparison against a run of F with a time step four times smaller than δt is shown. In all configurations Parareal converges reasonably well although the defect stalls for the first few iterations. As can be expected, computing more subintervals in parallel slows down convergence somewhat. For N t = 8, five iterations suffice to reduce the defect below the discretization error of the fine method, while for N t = 64 eight iterations are required. This confirms the usability of Parareal also for complex diffusion problems with anisotropic geometries and large jumps in the coefficients. Particularly the relatively mild increase in required iterations as N t is increased illustrates the potential for using larger numbers of cores (of the order of 100) to parallelize in time for this kind of problem, should a sufficiently large machine be available.

Effect of spatially varying coefficients
In the brick and mortar problem, the diffusion coefficients jump between D lip = 1 in the lipid channels and D cor = 10 −3 in the corneocytes. To assess the impact these jumps have on the convergence of Parareal, figure 4 gives a comparison of the defect for the brick and mortar problem (red) and a reference configuration with D lip = D cor = 10 −3 throughout the whole domain. Because of the smaller overall diffusivity, the  case of constant coefficients starts out with a larger defect but then converges quickly in three iterations. At k = 3 the defect is essentially zero. In contrast, the brick and mortar case -a 3D anisotropic problem with jumping coefficients - figure 4 shows a significant impact of the jumps on the rate of convergence. It has to be noted that in the simpler 2D problem studied in [27] this effect is much less pronounced.

Error over time
So far only the defect at the end of the simulation has been reported. In basic configurations, the defect typically increases monotonically with time so that controlling the defect at the end is already sufficient. For the complex problem studied here, however, the picture is more involved. Figure 5 shows both the defect d k n of Parareal for k = 1 (red) and k = 8 (green) as well as the estimated discretization error of the coarse (blue) and fine (yellow) integrator. The upper figure shows the defect for Parareal using N t = 64 concurrently computed subintervals, the lower figure for N t = 16.
Clearly, with just a single iteration, Parareal is more accurate than the coarse method only for the first couple of subintervals. Later subintervals are probably dominated by the error from the inaccurate starting value and are only as accurate as the coarse method. For eight iterations, Parareal is in both cases as accurate as the fine propagator, that is the defect is below the discretization error of F throughout the whole time interval [0, T ] (cf. the discussion in [4]). Note that the first eight subintervals are fully converged after eight iterations with a defect of the order of machine precision, and that because of the logarithmic scaling of the vertical axis some points are not showing up for the green line. It can also be observed that the defect and discretization error are high at the beginning and decay over time, as the solution approaches a steady state. However, the ordering of the lines in the beginning is roughly the same as in the end so that monitoring only the final defect seems to guarantee that Parareal reproduces the accuracy of the underlying fine method, although its particular value does not give an indication of the value at earlier times.
It is interesting to note that the computing time per subinterval decays in a similar way as the error. Where the error is high, the solution is far away from the steady state. Therefore, the solution from the previous time step is only a rough initial guess for the linear solver, which consequently requires many iterations in the beginning. As the error decays and the solution becomes more and more steady, the solution from the previous time step becomes a better and better initial guess, leading to a successive reduction in the number of required iterations and thus runtimes. This induces a load imbalance across time in Parareal (see the discussion in §3.4). F. Measured speedup follows the theoretical model reasonably well with some drop-off towards larger numbers of subintervals. It has to be noted that because of load imbalances Parareal performs significantly worse than for the ideal case described by (8). A short calculation shows that for the selected time step sizes and iteration numbers, a speedup of about S(64) ≈ 6.4 would be possible with optimal load balancing, in contrast to the actually measured value of S(64) ≈ 1.8. This clearly illustrates the need for a proper load balancing in time to achieve at least the theoretically possible speedup given by (8). Derivation and testing of such a time load balancing is left for future work.

Space-time parallel weak scaling
When doubling the spatial resolution, the resulting increase in the computational cost per time step can be compensated for by a corresponding increase in the number of cores used for the spatial parallelization - at least if both the employed method and implementation show proper weak scaling. For the spatial solver and parallelization of ug4 applied to the benchmark used here, this has been demonstrated successfully in [30]. However, when the number of fine time steps N f is also doubled, twice as many time steps have to be computed in order to keep the ratio δt/δx constant, which leads to a doubling of time-to-solution (see also the discussion in [4]). Time parallelization can provide some mitigation by also doubling the number of subintervals and thus of cores used to parallelize along time. Figure 7 compares the convergence behavior for two Parareal runs. The first one (blue) uses a time step size of ∆t = 1 /8 and δt = 1 /128 in units of T . The second one (red) on the other hand uses half the coarse and fine time step and half the spatial resolution so that δt/δx and ∆t/δx are the same in both runs. It also uses twice as many cores in time and eight times more cores in space, so that both the number of degreesof-freedom per core and the number of time steps per core remain the same, too. While the more accurate space-and time-discretization leads to a smaller defect at the beginning, convergence of Parareal in the configuration with more subintervals is significantly slower. This likely stems from the fact that Parareal is not optimal in the multi-grid sense [12], which means that the number of iterations typically goes up as the number of subintervals is increased, even if the total time interval [0, T ] remains the same. It should be noted that recent extensions of Parareal towards full time multi-grid have shown optimal behavior at least for simple diffusion problems [12]. Optimal behavior has also been demonstrated for the PFASST method [18] but again for a simple problem. As a skin permeation problem with, e.g., fully resolved lipid bi-layers will require a significant further increase in resolution, an optimal or near-optimal time-parallel method will be highly valuable.

Conclusions
Computational modeling of skin permeation is of interest for different applications. However, the full problem requires resolution of a vast range of scales, leading to enormous computational requirements. Massively parallel computers are needed but these require suitable parallel numerical methods to be used efficiently. This article investigates the applicability and performance of the time-parallel Parareal integrator to a relevant precursory problem of skin transport, i.e. a 3D brick and mortar configuration. For this the C++ Parareal library Lib4PrM is integrated as a plug-in into the simulation environment ug4 using implicit integrators and a geometric multi-grid as spatial solver. While the brick and mortar problem does not yet feature the same geometric level of detail as the skin transport problem, it already has jumps in the diffusion coefficients of several orders of magnitude on a highly anisotropic domain. The article is an extension of a previous study on a 2D problem on a domain with a simpler structure [27].
Performance of the space-time parallel solver is studied in several numerical experiments. It is confirmed that Parareal still converges reasonably fast for the brick and mortar case, although in contrast to the simple 2D case the jumps in the diffusion coefficients are found to have a noticeable impact. Furthermore, the defect of Parareal and the discretization error of the coarse and fine method are investigated. As the solution of the brick and mortar problem approaches a steady state in time, initial guesses for the geometric multi-grid become more accurate. This leads to faster convergence, which in turn induces an imbalance between the different processors in time. A basic model for the speedup in the imbalanced case is given and shown to agree reasonably well with actual speedup measurements. However, the imbalance is found to have a significant effect on parallel efficiency and further exacerbates Parareal's issues concerning maximum parallel efficiency. For Parareal with implicit integrators applied to complex PDEs, this load imbalance is an important issue that has to be addressed. In this context, weak scaling and optimality of Parareal, and their implications for the "trap of weak scaling" are briefly illustrated and discussed.
A number of possible directions for future research emerge from the experiments presented here besides load imbalances from the spatial solver and optimality. So far, coarsening in Parareal was done only in time by using a larger time step. Much better results can be expected if the spatial discretization is coarsened simultaneously. This requires a closer integration of Parareal with the spatial multi-grid solver, in order to provide interpolation and restriction routines. Also, this approach can be taken further by interweaving iterations of the time-parallel method with multi-grid cycles, as discussed e.g. for PFASST in [18]. Another important issue that is also connected to load balancing is spatial and temporal adaptivity. While both can in principle be used in Parareal, they greatly complicate the load balancing problem which, as shown here, is important to achieve reasonable efficiency. Finally, as energy consumption is becoming a more and more important issue in high-performance computing, a thorough benchmarking in terms of energy-to-solution is also an important direction for future work.