Nonlocal KKS: a convergent nonlocal extension of the KKS phase-field corrosion model

Companion code for the paper

A. Hermann, M. Fritz, T. Köppl, A. Shojaei, S. Silling, C. J. Cyron, Unifying local and nonlocal corrosion frameworks: A convergent nonlocal extension of the KKS phase-field model.

This repository implements an explicit, matrix-free solver for a nonlocal Kim–Kim–Suzuki (KKS) phase-field model of corrosion, in which the classical gradient operators are replaced by integral operators over a finite interaction horizon δ (peridynamic-style). The nonlocal integrals are evaluated with a one-point quadrature on a uniform grid. The code reproduces all numerical results of the paper: the asymptotic-compatibility study, the diffuse-to-sharp interface transition, the three-dimensional Mg–10Gd implant-screw application, and the parallel/GPU scaling.

Simulated biodegradation of a three-dimensional Mg–10Gd implant screw.

At fixed horizon and equal volume loss, increasing the double-well coefficient (ω = 1, 2, 4) sharpens the corrosion front from diffuse to sharp.

Requirements

• A C++17 compiler with OpenMP (e.g. g++ ≥ 7) — 2D/3D study drivers.
• An MPI compiler (e.g. Intel MPI mpicxx) — distributed-memory 3D solver.
• CUDA ≥ 11 and an NVIDIA V100-class GPU (sm_70) — optional GPU solver.
• Python ≥ 3.8 with numpy and matplotlib — figure generation.
• Blender ≥ 2.93 (optional) — 3D surface renderings of the screw.

Repository layout

include/        header-only solvers (2D and 3D) and quadrature rules
studies/        small C++ drivers, one per 2D numerical experiment
mpi/            distributed-memory 3D solver (nlkks_mpi3d.cpp) and CUDA port (nlkks_gpu3d.cu)
tools/          voxelize_screw.cpp — voxelize an implant-screw point cloud onto the grid
scripts/        Python post-processing that produces the paper figures
results/        bundled *sample* numerical output (small) used by the scripts
figures/        output directory for the generated PDFs

Building

make            # 2D/3D study drivers -> build/
make mpi        # distributed-memory 3D solver -> build/nlkks_mpi3d
make gpu        # CUDA 3D solver (set I_MPI_ROOT first)

Reproducing the figures

The small numerical outputs needed for most figures are bundled under results/, so the figures regenerate directly:

pip install numpy matplotlib      # if not already available
make figures                      # writes PDFs into figures/

This reproduces the operator- and system-level convergence plots, the terminal-state field maps, the screw volume-loss / degradation-rate plots, and the parallel-scaling and throughput plots.

The diffuse-to-sharp 2D maps and the screw cross-section slices depend on large raw fields (2D φ/c arrays and 3D VTK volumes) that are not bundled here; regenerate them with the drivers below, then run scripts/make_sharp_plots.py and scripts/make_screw_slices.py.

Regenerating the raw data from scratch

The 2D experiments are driven by the programs in studies/; each writes its output under results/ and is then consumed by the matching script in scripts/. For example:

make
./build/study_operator_conv      # operator-level symbol convergence
./build/study_mconv_2d           # full-system m- and delta-convergence
./build/study_snapshots          # terminal-state phi/c fields
./build/study_kks_sharp_2d       # diffuse-to-sharp interface maps

The three-dimensional Mg–10Gd screw is run with the distributed-memory solver. First voxelize the screw geometry, then run the solver with the passivation exponent as the last argument (M(t) = M0 * 10^(-l * VLmax%)):

make mpi
./build/voxelize_screw 0.008 0.3                 # -> screw mask (8 um voxels)
NLKKS_MASKFILE=screw_mask_fine.bin OMP_NUM_THREADS=48 \
  mpirun -n 16 -ppn 1 ./build/nlkks_mpi3d \
    326 326 576  0.33333333  1.0  999  0.1215  2600  0.02  4  0.0713  0.10
# args: Nx Ny Nz  h  delta  (unused)  dt  steps  M0  omega  cL  l_passivation

A single-GPU run uses build/nlkks_gpu3d with the same argument layout.

The large raw datasets (screw point cloud and voxel mask, full 2D fields, 3D VTK volumes, Blender renders) are archived separately as a data record on Zenodo (see Citing below).

Model and method (summary)

• Nonlocal Allen–Cahn equation for the phase field φ and a nonlocal Cahn–Hilliard-type equation for the concentration c, coupled through the KKS free energy.
• Nonlocal operator L_δ u(x) = K ∑_q ω_q [u(x) − u(x+ξ_q)] with a top-hat kernel of horizon δ; three one-point quadratures are provided — full-area (FA) and the partial-area corrections PA-HHB and PA-AC — and only the partial-area schemes are asymptotically compatible.
• Explicit (forward-Euler) time stepping; the stable step is governed by the horizon, Δt ≲ δ²/(4 H_m D), not by the mesh.

Citing

If you use this code, please cite the paper above and this software release:

A. Hermann et al., Nonlocal KKS solver (version 1.0.0), Zenodo, 2026, DOI: 10.5281/zenodo.20800551.

License

MIT — see LICENSE.