In metals, degradation occurs through the corrosion process, which usually consists of electrochemical reactions, including anodic and cathodic reactions as well as the formation of side products [299].

For Mg, the corrosion reactions comprise the following steps [299]: first, the material is released as metallic ions and free electrons, which causes the volume of the bulk material to be reduced:

The free electron reduces water to hydrogen gas and hydroxide ions:

Then, with the combination of the metallic and hydroxide ions, a porous film is formed on the surface, slowing down the degradation rate by protecting the material underneath:

With the presence of some specific ions in the surrounding medium, such as chloride ions in a saline solution, the protective film might be broken partially, which contributes to an increase of the rate of degradation:

The degradation process of metals is a continuous repetition of the above reactions.

In the example used in this study, the state variable in Eq. 7.5 is the concentration of effective chemical components involved in the degradation process, namely magnesium ions and the protective layer, denoted by $C_{\mathrm{Mg}}$ and $C_{\mathrm{Film}}$ respectively.

The defined saturation term acts as a function of space that varies between 0 and 1 in each point. By adding this term to the concentration of Mg ions, we can write

Since the film is a porous layer and allows the ions to diffuse through it, the diffusion coefficient in Eq. 7.10 is a function of space and not a constant value (which is the reason for being denoted as $D_{\mathrm{Mg}}^{e}$). We can calculate this effective diffusion function by interpolating two values at any point: 1) $D_{\mathrm{Mg}}^{e}=D_{\mathrm{Mg}}$ when $C_{\mathrm{Film}}=0$, and 2) $D_{\mathrm{Mg}}^{e}=\frac{\epsilon}{\tau}D_{\mathrm{Mg}}$ when $C_{\mathrm{Film}}=[\mathrm{Film}]_{\max}$, in which $\epsilon$ and $\tau$ are the porosity and tortuosity of the protective film, respectively. The interpolation leads to the effective diffusion function:

The level set method is a methodology that allows moving interfaces to be described by an implicit function. In other words, the boundaries of domains can be tracked as a function instead of being explicitly defined. In the level set method, a signed distance function, $\phi=\phi(x,y,z,t)$, describes the distance of each point in space to the interface, and the zero iso-contour of this function implies the interface [219]. In the current study, this function was defined in a way that divides the domain into two subdomains: 1) the bulk material, in which the implicit function is positive ($\phi>0$), and 2) the medium, in which the function is negative ($\phi<0$). The interface is defined as the points in space where $\phi=0$. Fig. 7.2 shows a schematic representation of the solid-medium interface in the current study, in which the interface moves as the material degrades over time.

Inserting the obtained velocity of Eq. 7.16 into Eq. 7.14 and considering the direction of the shrinkage velocity, which is in the opposite direction of the surface normal vector, yields

Eq. 7.17 is the final formulation of the level set equation in the current study, which alongside Eqs. 7.10 and 7.11 forms the mathematical model of degradation of Mg with a moving interface. Eq. 7.17 contributes indirectly to the evolution of Eqs. 7.10 and 7.11 as it defines the boundary, the zero iso-contour of the $\phi$ function, on which the boundary conditions of the equations are applied.

In order to solve Eq. 7.10 numerically, we used a finite difference scheme for the temporal term and a finite element formulation for the spatial terms. For simplicity of writing, notations of variables are changed, so $C_{\mathrm{Mg}}$ is represented as $u$ (the main unknown state variable to find), $C_{\mathrm{Film}}$ is denoted by $p$, $[\mathrm{Cl}]$ is denoted by $q$, and the saturation term $(1-\frac{F}{F_{\max}})$ is denoted by $s$. By doing this, Eq. 7.10 can be written as

To obtain the finite element formulation, the weak form of derived PDE is required. In order to get this, we define a space of test functions and then, multiply each term of the PDE by any arbitrary function as a member of this space. The test function space is

in which the $\Omega$ is the domain of interest, $\Gamma$ is the boundary of $\Omega$, and $\mathcal{H}^{1}$ denotes the Sobolev space of the domain $\Omega$, which is a space of functions whose derivatives are square-integrable functions in $\Omega$. The solution of the PDE belongs to a trial function space, which is similarly defined as

Then, we multiply Eq. 7.18 to an arbitrary function $v\in\mathcal{V}$:

One can approximate the unknown function $u$ in Eq. 7.29 by $u(x)\approx\sum_{i=0}^{N}c_{i}\psi_{i}(x)$, where the $\psi_{i}$ are the basis functions used to discretize the function space, and $c_{0},\ldots,c_{N}$ are the unknown coefficients. The finite element method uses Lagrange polynomials as the basis function and discretizes the computational domain using a new function space $\mathcal{V}_{h}$ spanned by the basis functions $\left\{\psi_{i}\right\}_{i\in\mathcal{I}_{s}}$, in which $\mathcal{I}_{s}$ is defined as $\mathcal{I}_{s}=\{0,\ldots,N\}$, where $N$ denotes the degrees of freedom in the computational mesh. The computational mesh discretizes the space into a finite number of elements, in each of which the $\psi_{i}$ is non-zero inside the $i$th element and zero everywhere else. In this study, 1st order Lagrange polynomials were used as the basis functions to define the finite element space.

For 1D elements, a 1st order Lagrange polynomial for the $i$th element with the width of $h$ can be written as:

A similar approach can be applied to define the basis function space in 2D and 3D spaces.

Eq. 7.32 is a linear system:

with

$M$ (which traditionally is called the mass matrix), $K$ (which traditionally is called the stiffness matrix), $f$, $c$, and $c_{1}$ are defined as:

By solving Eq. 7.33 and substituting the obtained $c$ in Eq. 7.31, $u$ ($C_{\mathrm{Mg}}$ in the example in this study) can be calculated in the current time step. As stated before, the same approach can be applied to Eq. 7.11 and Eq. 7.17 to get $C_{\mathrm{Film}}$ and $\phi$. This procedure is repeated in each time step to compute the values of $C_{\mathrm{Mg}}$, $C_{\mathrm{Film}}$, and $\phi$ over time.

A common practice to save time for solving Eq. 7.33 for a constant time step size is to compute the left-hand side matrix ($A$ in Eq. 7.34) once and compute only the right-hand side vector of the equation at each time iteration. But in this case, although the time step size is fixed, due to the presence of the $\alpha$ coefficient, the matrix changes along the time. The $\alpha$ coefficient is not constant and should be updated in each time step because it depends on the saturation term $s$ (which is a function of the concentration of the film as can be seen by comparing Eq. 7.10 and Eq. 7.18). In addition to this, the diffusion coefficient is not constant (Eq. 7.12), making the second term in Eq. 7.34 non-constant even in the absence of $\alpha$ coefficient. Consequently, the left-hand side matrix of the Eq. 7.33 cannot be computed before the start of the main time loop, and computing it in each time step is an extra but inevitable computational task in comparison to similar efficient and high-performance finite element implementations. This contributes to a slower algorithm for solving the aforementioned PDEs.

The model was implemented in FreeFEM [112], which is an open-source PDE solver to facilitate converting the weak formulation (Eq. 7.27) to a linear system $Ax=b$ (with $A$ from Eq. 7.34 and $b$ from Eq. 7.35). The computational mesh was generated using Netgen [233] in the SALOME platform [217] by a set of linear tetrahedral elements, and all the other preprocessing steps were performed in FreeFEM. The mesh was adaptively refined on the material-medium interface in order to increase the accuracy of the level set model. Postprocessing of the results was carried out using Paraview [8].

Computing the diffusion solely in the medium domain causes oscillations close to the interface, and to prevent this, the mass lumping feature of FreeFEM was employed. In this technique, the desired mass matrix is handled node-wise and not element-wise. Technically speaking, this means that the state variable is stored in the mesh nodes, and although this is the natural formulation in the finite difference method, it requires artificial modification in the standard finite element formulation [274]. The mass lumping feature of FreeFEM applies a quadratic formula at the vertices of elements to make the mass matrix diagonal, which contributes positively to the convergence of the solution.

The main parallelization approach for the current study was domain decomposition, in which the mesh is split into smaller domains (can be overlapping or non-overlapping), and the global solution of the linear system is achieved by solving the problem on each smaller local mesh. What really matters in this approach is providing virtual boundary conditions to the smaller sub-domains by ghost elements, transferring neighboring sub-domain solutions [21]. As a result, a high-performance parallelism is feasible by assigning each sub-domain to one processing unit.

In computational science, preconditioning is widely used to enhance the convergence, which means instead of directly working with a linear system $Ax=b$, one can consider the preconditioned system [62]:

in which the $M^{-1}$ is the preconditioner. In the current study, we considered this approach for both the domain composition and the solution of the linear system. We opted to use an overlapping Schwarz method for domain decomposition, in which the mesh is first divided into a graph of $N$ non-overlapping meshes using METIS (or ParMETIS) [144]. Then, by defining a positive number $\delta$, the overlapping decomposition $\left\{\mathcal{T}_{i}^{\delta}\right\}_{1\leqslant i\leqslant N}$ can be created recursively for each sub-mesh $\left\{\mathcal{T}_{i}\right\}_{1\leqslant i\leqslant N}$ by adding all adjacent elements of $\mathcal{T}_{i}^{\delta-1}$ to it. Then, the finite element space $\mathcal{V}_{h}$ (Eq. 7.19) can be mapped to the local space $\left\{\mathcal{V}_{i}^{\delta}\right\}_{1\leqslant i\leqslant N}$ by considering the restrictions $\left\{R_{i}\right\}_{1\leqslant i\leqslant N}$ and a local partition of unity $\left\{D_{i}\right\}_{1\leqslant i\leqslant N}$ such that:

where $I$ and $n$ denote identity matrix and the global number of unknowns, respectively [68].

In this study, we decomposed the mesh by using the one-level preconditioner Restricted Additive Schwarz (RAS):

in which $\left\{A_{i}\right\}_{1\leqslant i\leqslant N}$ is the local operator of the sub-matrices [68]. For this purpose, we took advantage of the HPDDM (high-performance domain decomposition methods) package interface in FreeFEM [140]. The partitioned mesh is shown in Fig. 7.3. The effect of the construction of these local sub-domains on the sparsity pattern of the global matrix is also depicted in Fig. 7.4. The global matrix is a sparse matrix according to Eq. 7.34 and the definition of the basis function $\psi$.

Generally, two categories of methods have been used to solve a large linear system of equations on parallel machines: direct solvers (e.g. Multifrontal Massively Parallel Sparse, MUMPS [14]) and iterative solvers (e.g. Generalized Minimal Residual Method, GMRES [223]). While direct solvers are quite robust, they suffer from the memory requirement problem on large systems. Inversely, iterative solvers are quite efficient on memory consumption, but similar to other iterative approaches, they are not very reliable in some cases [224]. Direct solvers modify the matrix by factorization (e.g. Cholesky decomposition), but an iterative solver does not manipulate the matrix and works solely using basic algebraic operations. However, for an efficient usage of iterative solvers, a proper preconditioner is crucial [224]. By evaluating and comparing the performance of the aforementioned methods for the current model, we decided to use an iterative approach using the Krylov subspaces (KSP) method, in which we preconditioned the equation using a proper preconditioner (Eq. 7.38) and then solved it with an iterative solver.

Krylov methods have been frequently used by researchers as robust iterative approaches to parallelism [131]. What matters in this regard is ensuring proper scaling of the parallelized algorithm for both the assembling of the matrices and the solution of the linear system of equations. One good solution to this challenge is taking advantage of HPC-ready mathematical libraries to achieve efficient distributed-memory parallelism through the Message Passing Interface (MPI). In the current study, we used the PETSc (Portable, Extensible Toolkit for Scientific Computation) library [25], which provides a collection of high-performance preconditioners and solvers for this purpose.

In order to yield the highest performance, a variety of different combinations of KSP types and preconditioners were evaluated, such as Conjugate Gradients (CG) [115], Successive Over-Relaxation (SOR) [104], block Jacobi, and Algebraic Multigrid (AMG) [187], to name a few. The performance tests results are presented in the appendix section of this chapter. The best performance for the reaction-diffusion system model was achieved using the HYPRE preconditioner [78] and the GMRES solver [223]. This was the combination used for all the performance analysis tests.

As mentioned before, in order to track the interface of the bulk material and the surrounding fluid, an implicit signed distance function is defined as the solution of Eq. 7.17. This equation can be solved using the aforementioned finite element discretization, but in a practical implementation, there are usually a couple of problems associated with this PDE.

The first issue is defining $D_{\mathrm{Mg}}^{e}$ and $\nabla_{n}C_{\mathrm{Mg}}$ on the moving interface. To ensure correct boundary conditions for Eq. 7.16, the value of $C_{\mathrm{Mg}}$ is set constant on the whole bulk material by using the penalty method. As a result, the implicit interface is not necessarily aligned on the computational mesh. Although this is a beneficial fact for the interface tracking, it inserts the problem of overestimation of $C_{\mathrm{Mg}}$ on the nodes close to the interface, which makes it difficult to calculate $\nabla_{n}C_{\mathrm{Mg}}$ on these nodes correctly. The same problem exists for calculating $D_{\mathrm{Mg}}^{e}$. To overcome this issue, the values of $C_{\mathrm{Mg}}$ and $D_{\mathrm{Mg}}^{e}$ are calculated at the distance $h$ from the interface in the normal direction (towards the medium), where $h$ is the edge size of the smallest element of the computational mesh.

The next issue is a well-known problem of the level set method: if the velocity of the interface is not constant (as in Eq. 7.13), the level set function $\phi$ may become distorted by having too flat or too steep gradients close to the moving front. This could cause unwanted movements of the interface. The problem becomes even worse when the distance function is advected. A solution to this issue is re-initializing the distance function in each time step (re-distancing), but this operation requires solving a new PDE. From numerical investigations, it has been observed that this operation inserts new errors in the numerical computation of the level set equation [222]. This can be resolved by improving the method of reconstruction of the distance function [222].

However, re-initialization results in another issue on a massively parallel implementation: as the mesh is partitioned into smaller sub-meshes, it is not feasible anymore to evaluate the distance to the interface globally on each sub-domain. As a result, the inverse process of domain decomposition should be taken to assemble the mesh again. This can be done by the restriction matrix and the partition of unity (defined in Eqs. 7.39 and 7.40), but it is rather a very inefficient procedure regarding the parallelization of the simulation and results in a long execution time in each time step.

In the current study, the distance function $\phi$ was initialized only once at the beginning of the simulation. The re-initialization process was unnecessary in this case because according to Eq. 7.17, the distance function is advected only in the regions where there is a gradient of the concentration of Mg ions, which means that advection is applied only on the regions close to the interface in the medium. This prevented the whole distance function of being distorted, and as a result, it was not required to re-initialize it in each time step. This also removed the need for inverting the decomposition process.

In order to verify the performance of the developed model, a degradation experiment was reconstructed in-silico, in which the degradation of a block of Mg (with the size of $13mm\times 13mm\times 4mm$) was investigated in a simulated body fluid solution. All the experimental parameter data (used to setup the simulation), as well as the degradation rates (used to calibrate and validate the numerical model) were extracted from [192].

As can be seen in Eqs. 7.1 and 7.2, each mole removed from the Mg block corresponds to one mole of the produced hydrogen. As a result, instead of a direct measurement of mass loss, one can collect and measure the amount of produced hydrogen to monitor the degradation rate. This is a common way of reporting degradation in this type of studies [4]. In order to get this quantity out of the developed model, we used the level set model output. The total mass loss of Mg at each desired time can be calculated based on the movement of the corrosion front:

where $\Omega_{+}(t)=\{\mathbf{x}:\phi(\mathbf{x},t)\geq 0\}$. It is worth noting that this integration should be performed by ignoring the ghost elements generated in the mesh partitioning process, otherwise the calculated material loss will be higher than the real value. Then, the amount of formed hydrogen gas can be calculated based on the ideal gas law:

in which $R$ is the universal gas constant, $P$ is the pressure, $T$ is the solution temperature, $A$ is the exposed corrosion surface area (which can be computed using the level set function), and $\mathrm{Mg}_{\mathrm{mol}}$ is the molar mass of Mg. Plotting a comparison of the predicted and experimentally obtained values of hydrogen can show the overall validity of the mathematical model because both the diffusion-reaction equations and the level set equation contribute to the prediction made by the computational model.

The geometry of the simulation experiment is depicted in Fig. 7.5. Based on this geometry, an Eulerian computational mesh was constructed by generating tetrahedral elements on the whole domain, including the Mg block and the medium. This resulted in ${830808}$ elements with a total of ${143719}$ DOFs for each PDE (Eqs. 7.10, 7.11, and 7.17), which indicates the size of matrix $A$ in Eq. 7.34. Model parameters and material properties were obtained from [23]. The diffusion coefficient of Mg was calculated using an inverse problem setup in which a Bayesian optimization process [196] was used to run the simulation code multiple times and minimize the difference of the model output and the experimental data reported by [192]. A time step convergence study was performed to measure the sensitivity of the model to the time stepping parameter, and based on the results, the time step value was set to 0.025 hours.

To investigate the performance and scaling behavior of the implemented parallel code, we conducted a set of weak-scaling and strong-scaling tests on the computational model. To do this, the time required to solve each PDE in each time step was measured in a simulation. This acted as a rough estimation of the time required in each time step because it ignores all the other factors contributing to speedup results such as communication costs, load imbalance, limited memory bandwidth, and parallelization-caused overhead.

Weak-scaling was evaluated by dividing the computational domain into smaller sub-domains (each of which was $\frac{1}{16}$ of the whole domain, Fig. 7.6) and conducting simulation experiments with 1, 2, 4, and 8 computational cores in a way that the number of processors corresponded to the number of employed sub-domains. In Fig. 7.6 the upper row shows different domains as an accumulation of the smaller divisions, and the lower row shows the corresponding domain decomposition for parallel computing by depicting each processing unit in a different color. In fact, it demonstrates the concept of increasing the number of MPI processing units as we increase the size of the problem.

After calculating the speedup of each test (by comparing the differences in execution time), we can use Gustafson’s law [100] to calculate the sequential and parallelizable portion of computation in the current implementation in weak-scaling evaluation:

where $N$ is the total number of computational cores, $f$ is the fraction of operations in the computation that are sequential, and as a result, $1-f$ is the fraction of the execution time spent on the parallelizable part.

The strong-scaling evaluation was performed using the entire domain. The evaluation was done using 1, 8, 16, 40, 60, 90, 200, and 300 MPI cores. In strong-scaling, Amdahl’s law [13] is used to calculate the portion of the algorithm that runs in parallel:

in which the parameters are the same as Eq. 7.43.

Simulations were conducted on the VSC (Flemish Supercomputer Center) supercomputer with the availability of Intel CPUs in three different micro-architectures: Ivy Bridge, Haswell, and Skylake. Due to a better performance, the strong and weak-scaling measurements were solely performed on the Skylake nodes. On this supercomputer, we made use of 3 nodes, 36 cores each, with 576 GB of the total memory, each node holding 2 Intel Xeon Gold 6132 CPUs with a base clock speed of 2.6 GHz. The nodes in the supercomputer are connected using an InfiniBand EDR network (bandwidth 25 Gb/s). For interprocess communication, Intel’s MPI implementation 2018 was used.

The performed numerical simulation produces the output of three main quantities: the concentration of the Mg ions in the medium (as the solution of Eq. 7.10), the concentration of the protective film (as the solution of Eq. 7.11), and the level set function values at each element (as the solution of Eq. 7.17). In addition to this, a quantitative prediction of the mass loss is also generated according to Eqs. 7.41 and 7.42.

In order to have quantitative predictions, the coefficients of Eqs. 7.10 and 7.11 (diffusion rates and reaction rates) should be calibrated using an inverse problem. Fig. 7.7 shows the results produced by the computational model after this parameter estimation stage. A narrow layer of the protective film is formed on the surface of the Mg block, and the volume of produced hydrogen gas is compared with values obtained from experiments. Additionally, by plotting the zero iso-contour of the level set function, we can obtain the shape of the material block as it degrades during the degradation process (i.e. tracking the moving corrosion front). This is depicted by the grey surface in Fig. 7.7.

Weak-scaling results are plotted in Fig. 7.8, in which the execution time of each time step is broken down into the time spent on each PDE. The results show good scalability of the parallel implementation.

Speedup and parallel efficiency of the weak-scaling experiment is plotted in Fig. 7.9. By fitting a curve based on the Gustafson equation (Eq. 7.43) on the obtained results (Fig. 7.9), the sequential proportion of the current implementation was calculated to be $18\%$, which means that $82$ percent of the code can be parallelized, which is a proper but not an ideal scalability.

The strong-scaling results are plotted in Fig. 7.10, which shows a better scalability in comparison to the weak-scaling test. For a better representation, exact measured values are presented in Table 7.1.

Similar to weak-scaling results, Fig. 7.11 demonstrates the speedup and parallel efficiency of the developed code for strong-scaling evaluation. From the results, it is obvious that increasing the number of cores leads to a better performance but a lower efficiency. By fitting Amdahl’s equation (Eq. 7.44) on the obtained speedup results (Fig. 7.11), $f$ was obtained as $0.01$, which means in strong-scaling terms that $99\%$ of the code can run in parallel.