How to Write Efficient C++ Simulators Using Arrays of Structures Not Structures of Arrays in OpenCSMP
Stephan Konrad Matthai · 2024
Summary Multi-physics simulation of complex geologic systems needs to express dynamic behaviour manifesting as plumes, convection cells and-or fracturing and faulting. The “playing field” typically is a feature-rich geomodel. Instabilities emerge across multiple length scales, and fracturing creates complex patterns with power-law – frequency distributions and cross-cutting relationships revealing sequential formation over geologic time. The Complex Systems Modelling Platform (CSP then CSMP++, now OpenCSMP) employs a space-time adaptive FE/FV/FD discretisation to model this behaviour, with a focus on geo-energy systems. Emergent features are captured in different ways, for instance via the introduction of discontinuities / interfaces. Unstructured mesh patches are inserted into a structured mesh to transit from coarse to fine regions. Unstructured tetrahedral /triangular domains connect with hexahedral structured domains using prisms and pyramids. This discretisation is implemented with dynamic polymorphic data structures. Combined with a distributed property storage the global impact of local dynamic mesh adaptation is minimised. The current code uses policy-based class design, and elements, nodes and other data structures are implemented as class templates. Notwithstanding, memory footprint per computational-degrees-of-freedom matches commercial codes with an array-based rigid mesh representation. Here we review CSMP’s journey through alternative implementations of its data model and FEM/FVM computations, discussing pros and cons of design choices and code readability and complexity. CSMP++ now underpins multiple simulators, integrating functionality that has enabled the numeric simulation research in 250 peer-reviewed articles, including letters in Science and Nature Geoscience. We observe that computational speed is limited data access and CPU cache throughput. Yet there are very significant differences between optimised code of competing class template designs including memory footprint. During refactoring we managed to shift some dynamic dispatching into the compile time domain, benefitting from move semantics and compiler R-value optimisations for classes allocated on the stack as opposed to the heap. Regarding polymorphism we settled on a compromise using both static and dynamic polymorphism in different parts of the code, instantiating objects on the stack wherever possible, benefitting from compiler elimination of virtual function calls. In the analysis conducted here we confirm the benefit of flattening nested arrays and identify areas for improvement.