Graphics Processing Units Research Articles

GPU parallel implementation of a finite volume lattice Boltzmann method for incompressible flows

This work presents a graphics processing units (GPU) parallel algorithm of a cell-centered finite volume lattice Boltzmann method (FVLBM) on unstructured meshes. In the present GPU parallel algorithm, the parallelization is performed in the physical space. To reduce the frequency of GPU memory accesses, this algorithm develops coalesced access to GPU memory. In addition, to avoid the race for resources leading to data anomalies, such as dirty read or phantom read etc., and the double counting for flux calculation, the efficient face-based data structure often used for flux calculation in cells in the central processing unit (CPU) version of FVLBM is modified into a face-based data structure used for the fluxes on all faces, followed by a cell-based loop for the final residuals in all cells. Therefore, the proposed GPU parallel algorithm does not need to use the resource lock and retains the high efficiency of the face-based data structure in the fluxes computation to enhance its’ parallel efficiency. Additionally, to demonstrate the computational efficiency of the proposed GPU parallel algorithm, various benchmark studies are performed in this work by the proposed parallel scheme on a double precision NVIDIA GeForce RTX 3090Ti GPU card, including (a) the lid-driven flow in a two-dimensional (2D) square cavity, (b) a 2D flow past a cylinder, and (c) the lid-driven flow in a three-dimensional (3D) cubic cavity. The numerical results show that the proposed GPU parallel algorithm can be as accurate as the original CPU serial scheme with 1 to 2 orders of speedup.

High-speed turbulent flows towards the exascale: STREAmS-2 porting and performance

Exascale High Performance Computing (HPC) represents a tremendous opportunity to push the boundaries of Computational Fluid Dynamics (CFD), but despite the consolidated trend towards the use of Graphics Processing Units (GPUs), programmability is still an issue. STREAmS-2 (Bernardini et al. Comput. Phys. Commun. 285 (2023) 108644) is a compressible solver for canonical wall-bounded turbulent flows capable of harvesting the potential of NVIDIA GPUs. Here we extend the already available CUDA Fortran backend with a novel HIP backend targeting AMD GPU architectures. The main implementation strategies are discussed along with a novel Python tool that can generate the HIP and CPU code versions allowing developers to focus their attention only on the CUDA Fortran backend. Single GPU performance is analysed focusing on NVIDIA A100 and AMD MI250x cards which are currently at the core of several HPC clusters. The gap between peak GPU performance and STREAmS-2 performance is found to be generally smaller for NVIDIA cards. Roofline analysis allows tracing this behavior to unexpectedly different computational intensities of the same kernel using the two cards. Additional single-GPU comparisons are performed to assess the impact of grid size, number of parallelized loops, thread masking and thread divergence. Parallel performance is measured on the two largest EuroHPC pre-exascale systems, LUMI (AMD GPUs) and Leonardo (NVIDIA GPUs). Strong scalability reveals more than 80% efficiency up to 16 nodes for Leonardo and up to 32 for LUMI. Weak scalability shows an impressive efficiency of over 95% up to the maximum number of nodes tested (256 for LUMI and 512 for Leonardo). This analysis shows that STREAmS-2 is the perfect candidate to fully exploit the power of current pre-exascale HPC systems in Europe, allowing users to simulate flows with over a trillion mesh points, thus reducing the gap between the Reynolds numbers achievable in high-fidelity simulations and those of real engineering applications.

Graphics card-based real-time processing for dual comb interferometry.

The technique of performing interferometry with two optical frequency combs is used by an increasing number of research groups and even for field deployed commercial applications. Real-time interferogram acquisition, correction, and averaging are, however, still not broadly accessible. This limits the deployment and wider adoption of this high resolution, high sensitivity technique. We herein introduce and describe a freely available correction software performing real-time processing on a graphics processing unit.

Mesh Convolution With Continuous Filters for 3-D Surface Parsing.

Geometric feature learning for 3-D surfaces is critical for many applications in computer graphics and 3-D vision. However, deep learning currently lags in hierarchical modeling of 3-D surfaces due to the lack of required operations and/or their efficient implementations. In this article, we propose a series of modular operations for effective geometric feature learning from 3-D triangle meshes. These operations include novel mesh convolutions, efficient mesh decimation, and associated mesh (un)poolings. Our mesh convolutions exploit spherical harmonics as orthonormal bases to create continuous convolutional filters. The mesh decimation module is graphics processing unit (GPU)-accelerated and able to process batched meshes on-the-fly, while the (un)pooling operations compute features for upsampled/downsampled meshes. We provide an open-source implementation of these operations, collectively termed Picasso. Picasso supports heterogeneous mesh batching and processing. Leveraging its modular operations, we further contribute a novel hierarchical neural network for perceptual parsing of 3-D surfaces, named PicassoNet++. It achieves highly competitive performance for shape analysis and scene segmentation on prominent 3-D benchmarks. The code, data, and trained models are available at https://github.com/EnyaHermite/Picasso.

Fast interactive simulations of cardiac electrical activity in anatomically accurate heart structures by compressing sparse uniform cartesian grids

Background and Objective:Numerical simulations are valuable tools for studying cardiac arrhythmias. Not only do they complement experimental studies, but there is also an increasing expectation for their use in clinical applications to guide patient-specific procedures. However, numerical studies that solve the reaction–diffusion equations describing cardiac electrical activity remain challenging to set up, are time-consuming, and in many cases, are prohibitively computationally expensive for long studies. The computational cost of cardiac simulations of complex models on anatomically accurate structures necessitates parallel computing. Graphics processing units (GPUs), which have thousands of cores, have been introduced as a viable technology for carrying out fast cardiac simulations, sometimes including real-time interactivity. Our main objective is to increase the performance and accuracy of such GPU implementations while conserving computational resources. Methods:In this work, we present a compression algorithm that can be used to conserve GPU memory and improve efficiency by managing the sparsity that is inherent in using Cartesian grids to represent cardiac structures directly obtained from high-resolution MRI and mCT scans. Furthermore, we present a discretization scheme that includes the cross-diagonal terms in the computational cell to increase numerical accuracy, which is especially important for simulating thin tissue sections without the need for costly mesh refinement. Results:Interactive WebGL simulations of atrial/ventricular structures (on PCs, laptops, tablets, and phones) demonstrate the algorithm’s ability to reduce memory demand by an order of magnitude and achieve calculations up to 20x faster. We further showcase its superiority in slender tissues and validate results against experiments performed in live explanted human hearts. Conclusions:In this work, we present a compression algorithm that accelerates electrical activity simulations on realistic anatomies by an order of magnitude (up to 20x), thereby allowing the use of finer grid resolutions while conserving GPU memory. Additionally, improved accuracy is achieved through cross-diagonal terms, which are essential for thin tissues, often found in heart structures such as pectinate muscles and trabeculae, as well as Purkinje fibers. Our method enables interactive simulations with even interactive domain boundary manipulation (unlike finite element/volume methods). Finally, agreement with experiments and ease of mesh import into WebGL paves the way for virtual cohorts and digital twins, aiding arrhythmia analysis and personalized therapies.

Parallel and scalable AI in HPC systems for CFD applications and beyond

This manuscript presents the library AI4HPC with its architecture and components. The library enables large-scale trainings of AI models on High-Performance Computing systems. It addresses challenges in handling non-uniform datasets through data manipulation routines, model complexity through specialized ML architectures, scalability through extensive code optimizations that augment performance, HyperParameter Optimization (HPO), and performance monitoring. The scalability of the library is demonstrated by strong scaling experiments on up to 3,664 Graphical Processing Units (GPUs) resulting in a scaling efficiency of 96%, using the performance on 1 node as baseline. Furthermore, code optimizations and communication/computation bottlenecks are discussed for training a neural network on an actuated Turbulent Boundary Layer (TBL) simulation dataset (8.3 TB) on the HPC system JURECA at the Jülich Supercomputing Centre. The distributed training approach significantly influences the accuracy, which can be drastically compromised by varying mini-batch sizes. Therefore, AI4HPC implements learning rate scaling and adaptive summation algorithms, which are tested and evaluated in this work. For the TBL use case, results scaled up to 64 workers are shown. A further increase in the number of workers causes an additional overhead due to too small dataset samples per worker. Finally, the library is applied for the reconstruction of TBL flows with a convolutional autoencoder-based architecture and a diffusion model. In case of the autoencoder, a modal decomposition shows that the network provides accurate reconstructions of the underlying field and achieves a mean drag prediction error of ≈5%. With the diffusion model, a reconstruction error of ≈4% is achieved when super-resolution is applied to 5-fold coarsened velocity fields. The AI4HPC library is agnostic to the underlying network and can be adapted across various scientific and technical disciplines.

A piecewise-hierarchical particle count method suitable for the implementation of the unified gas-kinetic wave–particle method on graphics processing unit devices

A piecewise-hierarchical particle count method suitable for the implementation of the unified gas-kinetic wave–particle method on graphics processing unit devices

Deflagration-to-detonation transition and detonation propagation in supersonic flows with hydrogen injection and downstream ignition

This study investigates the mechanisms of flame acceleration and deflagration-to-detonation transition (DDT) in supersonic flows using transverse hydrogen injection and downstream ignition. Utilizing the graphics processing unit accelerated adaptive mesh refinement approach, we examine the influence of downstream ignition jet pressure on DDT through high-resolution computational simulations. Our results indicate that the transverse injection of hydrogen into the supersonic mainstream generates strong turbulence and numerous vortices due to Kelvin–Helmholtz instability, enhancing fuel mixing efficiency along the flow but deviating from the ideal premixed state. Following the injection of the downstream ignition jet into the supersonic main flow, initial flame acceleration is less effective than in the premixed state due to the non-uniformity of the incoming flow. However, within the boundary layer, the flame remains stable, and the intense turbulence fosters shock–flame interactions. The convergence of multiple compression waves into a shock wave facilitates energy deposition, coupling with the flame to trigger local detonation via the reactive gradient mechanism. The detonation wave exhibits complex wavefront structures, including vertical and oblique fronts induced by boundary layer interactions. Ignition jet pressure significantly impacts the DDT process and detonation wave characteristics, reducing ignition time and affecting the detonation temperature, pressure, and propagation speed. This study provides valuable insights into the dynamics of flame acceleration and DDT in supersonic flows with non-uniform fuel distribution and downstream jet ignition. The findings highlight the critical role of ignition jet pressure in optimizing ignition and detonation processes, offering new perspectives for achieving low-energy, rapid detonation initiation within the tube.

Hybrid finite element and laplace transform method for efficient numerical solutions of fractional PDEs on graphics processing units

Fractional Partial Differential equations (FPDEs) are essential for modeling complex systems across various scientific and engineering areas. However, efficiently solving FPDEs presents significant computational challenges due to their inherent memory effects, often leading to increased execution times for numerical solutions. This study proposes a highly parallelizable hybrid computational approach that combines the Finite Element Method (FEM) for spatial discretization with Numerical Inversion of the Laplace Transform (NILT) for time-domain solutions, optimized for execution on Graphics Processing Units (GPUs). The NILT method’s high parallelizability, stemming from the independence of its series terms, combined with the robust spatial discretization provided by FEM, enables the efficient and accurate solution of FPDEs on GPUs, demonstrating substantial performance improvements over traditional CPU-based implementations. We observe a generalized pattern in execution time behavior that accounts for both the number of nodes and the number of NILT terms. Specifically, execution time scales quadratically with the number of nodes, while showing only a logarithmic marginal increase with the number of NILT terms These behaviors not only enables consistent performance assessment but also highlights potential areas for algorithm optimization. Validation against exact solutions of fractional diffusion and wave equations, employing Caputo, modified Caputo-Fabrizio, and modified Atangana-Baleanu derivatives, demonstrates the accuracy and convergence of the hybrid FEM-NILT method. Notably, the exact solutions of wave equation are novel in literature. The results highlight the method’s potential for enabling high-precision, large-scale simulations in fractional calculus applications, thereby advancing computational capabilities and efficiency in the field.

Blood Vessel Segmentation and Classification of Diabetic Retinopathy with Machine Learning-Based Ensemble Model

The incidence of diabetes has increased in recent times due to factors such as obesity and genetic predisposition. Diabetes wears out the eye vessels over time. Diabetic retinopathy (DR) is a serious disease that leads to vision problems. DR can be diagnosed by specialists who examine the fundus images of the eye at regular intervals. With 537 million diabetics in 2021, this method can be time-consuming, costly and inadequate. Artificial intelligence algorithms can provide fast and cost-effective solutions for DR diagnosis. In this study, the noise of blood vessels in fundus images was eliminated using the LinkNet-RCB7 model, and diabetic retinopathy was categorized into five classes using a machine learning-based ensemble model. Artificial intelligence-based classification training using images as input takes a long time and requires high resource requirements such as Random Access Memory (RAM) and Graphics Processing Unit (GPU). By using Gray Level Cooccurrence Matrix (GLCM) attributes in the classification phase, a lower resource requirement was aimed for. A Dice coefficient of 85.95% was achieved for the segmentation of blood vessels in the Stare dataset, in addition to 97.46% accuracy for binary classification and 96.10% accuracy for classifying DR into five classes in the dataset APTOS 2019.

МЕТОД ОБУЧЕНИЯ ИНТЕЛЛЕКТУАЛЬНОГО АГЕНТА С ПОМОЩЬЮ СЕТЕЙ DOUBLE DQN, ПУТЕВЫХ ТОЧЕК И ФУНКЦИИ ВОЗНАГРАЖДЕНИЯ

The problems of increasing the control efficiency of autonomous vehicles are considered. The problem of reducing the required computational resources for the intelligent vehicle control module is highlighted. A neural network training algorithm for Double DQN architecture with modified reward functions is proposed. The basis of the proposed solution is the use of lane segmentation, reward function and the use of additional waypoints in training. A software model has been developed and simulation of the learning process has been performed. The results obtained from a comparative analysis with known solutions show a stable increase in episode duration, and effective training in a realistic urban simulation. The study points to the possibility of reducing the need for high computing power, which will enable the use of central processing units (CPUs) for basic functions of unmanned vehicles instead of graphics processing units (GPUs).

Efficient Training Data Acquisition Technique for Deep Learning Networks in Radar Applications

In the field of radar, deep learning techniques have shown considerably superior performance over traditional classifiers in detecting and classifying targets. However, acquiring sufficient training data for deep learning applications is often challenging and time consuming. In this study, we propose a technique for acquiring training data efficiently using a combination of synthesized data and measured background data. We utilized graphics processing unit (GPU)-based physical optics methods to obtain the backscattered field of moving targets. We then generated a virtual dataset by mixing the synthesized target signal with the background signal real data. Subsequently, we trained a convolutional neural network using the virtual dataset to identify three different classes—Bird, Drone, and Background—from a range-Doppler map. When tested using the measurement data, the trained model achieved an accuracy of over 90%, demonstrating the effectiveness of the proposed method in acquiring training data for radar-based deep learning applications.

A graphics processing unit accelerated sparse direct solver and preconditioner with block low rank compression

We present the GPU implementation efforts and challenges of the sparse solver package STRUMPACK. The code is made publicly available on github with a permissive BSD license. STRUMPACK implements an approximate multifrontal solver, a sparse LU factorization which makes use of compression methods to accelerate time to solution and reduce memory usage. Multiple compression schemes based on rank-structured and hierarchical matrix approximations are supported, including hierarchically semi-separable, hierarchically off-diagonal butterfly, and block low rank. In this paper, we present the GPU implementation of the block low rank (BLR) compression method within a multifrontal solver. Our GPU implementation relies on highly optimized vendor libraries such as cuBLAS and cuSOLVER for NVIDIA GPUs, rocBLAS and rocSOLVER for AMD GPUs and the Intel oneAPI Math Kernel Library (oneMKL) for Intel GPUs. Additionally, we rely on external open source libraries such as SLATE (Software for Linear Algebra Targeting Exascale), MAGMA (Matrix Algebra on GPU and Multi-core Architectures), and KBLAS (KAUST BLAS). SLATE is used as a GPU-capable ScaLAPACK replacement. From MAGMA we use variable sized batched dense linear algebra operations such as GEMM, TRSM and LU with partial pivoting. KBLAS provides efficient (batched) low rank matrix compression for NVIDIA GPUs using an adaptive randomized sampling scheme. The resulting sparse solver and preconditioner runs on NVIDIA, AMD and Intel GPUs. Interfaces are available from PETSc, Trilinos and MFEM, or the solver can be used directly in user code. We report results for a range of benchmark applications, using the Perlmutter system from NERSC, Frontier from ORNL, and Aurora from ALCF. For a high frequency wave equation on a regular mesh, using 32 Perlmutter compute nodes, the factorization phase of the exact GPU solver is about 6.5× faster compared to the CPU-only solver. The BLR-enabled GPU solver is about 13.8× faster than the CPU exact solver. For a collection of SuiteSparse matrices, the STRUMPACK exact factorization on a single GPU is on average 1.9× faster than NVIDIA’s cuDSS solver.

A GPU based 3D raytracing algorithm for DUED laser fusion code

Abstract These days, graphical processing units (GPUs) deliver performance
comparable to that of hundreds of CPU cores. This level of performance allows
certain classes of simulations to be run in-house on a standard consumer workstation,
eliminating the need for a cluster. In this paper, it is shown that medium-resolution,
2D radiation hydrodynamics simulations for laser-driven inertial confinement fusion
with realistic 3D laser raytracing can now be conducted on a single consumer device.
A novel raytracing module has indeed been developed for the 2D Lagrangian radiation-
hydro-nuclear code DUED to leverage the computational power of GPUs. By
employing 3D raytracing, more realistic investigations of laser-driven plasmas become
feasible, with a particular focus on perturbations resulting from non-uniform laser
irradiation.

CGH calculation algorithm for expressing reflection on a curved mirror surface

Rendering techniques are important in computer-generated holograms (CGHs) for expressing various types of holo-realistic 3D images. Rendering techniques such as hidden surface removal and reflection by a planar mirror have been proposed thus far, but reflection on a curved mirror surface has yet to be achieved. In this study, we propose a calculation algorithm that can express the reflection of the surrounding environment on a mirror surface defined as a Bézier surface. The results of optical experiments demonstrate that the proposed algorithm enables reflection on such a mirror surface and that the calculation can be accelerated by using a graphics processing unit (GPU).

Cascaded lattice Boltzmann simulation of Newtonian and non-Newtonian mixture nanofluids with variable thermophysical properties in a cavity with vertical heat radiator

The central moments-based cascaded lattice Boltzmann method (CLBM) for then Newtonian and non-Newtonian Buongiorno’s model mixture nanofluids (CuO, ZnO, Al2O3-water) has been implemented and applied in a square chamber with a vertical heat radiator using accelerated graphics processing unit (GPU) computing through compute unified device architecture (CUDA) C/C++ platform. Due to the higher numerical stability, the CLBM is a superior numerical tool to the raw moments-based MRT-LBM (multiple-relaxation-time lattice Boltzmann method). Three different models for the viscosity and thermal conductivity of the nanofluids: (i) the Binkmann model for the constant viscosity and the Maxwell model for the constant thermal conductivity (ii) Binkmann and Maxwell model with temperature dependent Brownian motion and (iii) Corcione model with non-Newtonian fluid where the temperature and strain rate determine the nanofluid effective thermal conductivity and viscosity, have been used. The enclosure’s upper and bottom walls are thermally adiabatic, but the left and right walls are uniformly cold. A vertical heater is immersed in the middle position of the cavity. The benchmark results for non-Newtonian, Newtonian, and nanofluids for the various computational domains are used to validate the current code adequately. The Bingham number (Bn), the Rayleigh number (Ra), and The volume fraction of the nanoparticles (ϕ) are the three key parameters that are varied in this investigation to demonstrate the effects of natural convection on the isotherms, streamlines, isolines of nanoparticle volume fractionation, yielded and unyeilded zone, and average Nusselt number (Nu¯). The Brownian motion effects of the nanoparticles augmented the average rate of heat transfer and the use of the Bingham nanofluids reduced the heat transfer enhancement. For the CuO-water nanofluid, the augmentation of the rate of heat transfer is 15.42% from ϕ=0 to 4% while Ra=106 and the corresponding heat transfer enhancement for the ZnO-water nanofluid is 11.11%. For the Bingham fluid, the rate of heat transfer increases 7% from Ra=105 to Ra=106 while ϕ=2% and Bn=0.3. The findings can be applied to optimize automotive radiator systems, which are crucial for maintaining engine temperature.

Enhanced seagull optimization for enhanced accuracy in CUDA-accelerated Levenberg–Marquardt backpropagation neural networks for earthquake forecasting

Hyperparameter tuning is crucial for enhancing the accuracy and reliability of artificial neural networks (ANNs). This study presents an optimization of the Levenberg–Marquardt backpropagation neural network (LM-BPNN) by integrating an improved seagull optimization algorithm (ISOA). The proposed ISOA-LM-BPNN model is designed to forecast earthquakes in the Caribbean region. The study further explores the impact of data and model parallelism, revealing that hybrid parallelism effectively mitigates the limitations of both. This leads to substantial gains in throughput and overall performance. To address computational demands, this model leverages the compute unified device architecture (CUDA) framework, enabling hybrid parallelism on graphics processing units (GPUs). This approach significantly enhances the model’s computational speed. The experimental results demonstrate that the ISOA-LM-BPNN model achieves a 20% improvement in accuracy compared to four baseline algorithms across three diverse datasets. The integration of ISOA with LM-BPNN refines the neural network’s hyperparameters, leading to more precise earthquake predictions. Additionally, the model’s computational efficiency is evidenced by a 56% speed increase when utilizing a single GPU, and an even greater acceleration with dual GPUs connected via NVLink compared to traditional CPU-based computations. The findings underscore the potential of ISOA-LM-BPNN as a robust tool for earthquake forecasting, combining high accuracy with enhanced computational speed, making it suitable for real-time applications in seismic monitoring and early warning systems.

Multi-Stream Scheduling of Inference Pipelines on Edge Devices - a DRL Approach

Low-power edge devices equipped with Graphics Processing Units (GPUs) are a popular target platform for real-time scheduling of inference pipelines. Such application-architecture combinations are popular in Advanced Driver-assistance Systems for aiding in the real-time decision-making of automotive controllers. However, the real-time throughput sustainable by such inference pipelines is limited by resource constraints of the target edge devices. Modern GPUs, both in edge devices and workstation variants, support the facility of concurrent execution of computation kernels and data transfers using the primitive of streams , also allowing for the assignment of priority to these streams. This opens up the possibility of executing computation layers of inference pipelines within a multi-priority, multi-stream environment on the GPU. However, manually co-scheduling such applications while satisfying their throughput requirement and platform memory budget may require an unmanageable number of profiling runs. In this work, we propose a Deep Reinforcement Learning (DRL)-based method for deciding the start time of various operations in each pipeline layer while optimizing the latency of execution of inference pipelines as well as memory consumption. Experimental results demonstrate the promising efficacy of the proposed DRL approach in comparison with the baseline methods, particularly in terms of real-time performance enhancements, schedulability ratio, and memory savings. We have additionally assessed the effectiveness of the proposed DRL approach using a real-time traffic simulation tool IPG CarMaker.

Correlation mining of multimodal features based on higher-order partial least squares for emotion recognition in conversations

In fields requiring an understanding of emotions, such as digital human interaction and public opinion analysis, achieving a dependable and interpretable model for mining correlations among multimodal features remains a primary objective. However, current deep learning methods often lack transparency and suffer from low interpretability. To address these challenges, we propose a novel Correlation Mining method based on Higher-Order Partial Least Squares (HOPLS) for multimodal Emotion Recognition in conversations (CMHER) in this paper. CMHER innovatively combines HOPLS with transformers and Gated Recurrent Units (GRUs) to compute correlation matrices within unimodal data streams and between cross-modal sources. HOPLS projects source data into a latent space to predict target data via correlation matrix computations, eliminating the need for Graphical Processing Unit (GPU) acceleration and making it suitable for experimental and edge systems. The integration of HOPLS with deep neural networks involves preprocessing multimodal features into suitable dimensions and latent representations, followed by HOPLS computing correlation matrices for cross-modal latent vectors and final labels through optimal joint subspace approximation, which aims at the improvements of both interpretability and reliability. Additionally, a generalization error fitting module further refines the predicted correlation matrices to improve predictive capability and overall model performance. Experiments on two public datasets validate the superiority of our proposed method.

Thermodynamic quantum Fokker–Planck equations and their application to thermostatic Stirling engine

We developed a computer code for the thermodynamic quantum Fokker–Planck equations (T-QFPE), derived from a thermodynamic system–bath model. This model consists of an anharmonic subsystem coupled to multiple Ohmic baths at different temperatures, which are connected to or disconnected from the subsystem as a function of time. The code numerically integrates the T-QFPE and their classical expression to simulate isothermal, isentropic, thermostatic, and entropic processes in both quantum and classical cases. The accuracy of the results was verified by comparing the analytical solutions of the Brownian oscillator. In addition, we illustrated a breakdown of the Markovian Lindblad-master equation in the pure quantum regime. As a demonstration, we simulated a thermostatic Stirling engine employed to develop non-equilibrium thermodynamics [S. Koyanagi and Y. Tanimura, J. Chem. Phys. 161, 114113 (2024)] under quasi-static conditions. The quasi-static thermodynamic potentials, described as intensive and extensive variables, were depicted as work diagrams. In the classical case, the work done by the external field is independent of the system–bath coupling strength. In contrast, in the quantum case, the work decreases as the coupling strength increases due to quantum entanglement between the subsystem and bath. The codes were developed for multicore processors using Open Multi-Processing (OpenMP) and for graphics processing units using the Compute Unified Device Architecture. These codes are provided in the supplementary material.