Detection of magnetohydrodynamic waves by using convolutional neural networks

Shock waves are common features in high-speed compressible flows, where the fluid state (density, pressure, etc.) changes drastically over a thin region whose thickness is proportional to the mean free path of the gas. Because of the importance of shock waves in several applications, such as aircraft aerodynamics, astrophysical phenomena¹ and inertial confinement fusion (ICF),² etc., the detection of shock waves is an important and challenging problem in computational fluid dynamics (CFD).³ Several shock detection methods have been developed for hydrodynamic shocks. Pagendarm and Seitz⁴ proposed a shock detection method based on finding the density gradient maxima. Although the method is easy to understand, proper filters must be carefully designed to remove false results. Samtaney⁵ proposed the zero-crossing of the Laplacian of the density field to detect shocks. Since the normal direction of shock wave is parallel to a local pressure gradient and, thus, normal Mach number can be obtained from a pressure distribution, Lovely and Haimes⁶ proposed a method using the iso-surface of unity normal Mach number to represent the detected shock wave surface. The above methods are relatively easily implemented, but their effectiveness and accuracy require further improvement. Kanamori and Suzuki^7,8 proposed a shock detection method by calculating critical lines for the vector field of the characteristics. This approach can be applied in steady and unsteady two- (2D) and three-dimensional (3D) flows. A shock-wave-detection technique for continuum and rarefied-gas flows has been proposed using Schlieren imaging.⁹ The scheme is applicable for any existing 2D flow fields obtained experimentally or numerically. Samtaney et al.¹⁰ proposed a method to extract shocklets in turbulence simulations based on finding the minimum for a cost function defined in terms of local Rankine–Hugoniot jump conditions. Fujimoto et al.^11,12 recently offered an alternative shock detection method by integrating canny edge detection, which is one of the image processing methods to detect edges, and Rankine–Hugoniot relations. These two methods provide good accuracy, but implementation is somewhat complicated.

In magnetohydrodynamics (MHD), a wave is considered physical only if it satisfies both viscosity admissibility condition and evolutionary condition.¹³ In a strongly planar system, the fast and slow waves (shocks or expansion fans), contact discontinuities, $1→3$ and $2→4$ intermediate shocks, slow and fast compound waves are considered admissible, while fast and slow waves (shocks or expansion fans), contact discontinuities, and 180° rotational discontinuity (RD) are admissible in the planar system. Note that a flow is planar if there is no derivatives in the out-of-plane (z) direction and strongly planar if there is also a reference frame in which there is no vector component in the z-direction. Therefore, shock wave detection in MHD becomes more challenging due to the variety of wave types.^13,14 Wheatley et al.¹⁵ discussed in detail the variety of waves resulting from shock refraction in MHD. Most shock detection methods developed for hydrodynamics have not yet been adopted or extended to MHD due to myriad MHD wave types, which cause these methods to fail classification and detection. MHD wave classification is traditionally based on the Rankine–Hugoniot relations, producing a set of nonlinear algebraic relations and compatibility relationships for a complex wave structure. Snow et al.¹⁶ proposed a detection method for MHD shocks based on upstream and downstream velocities relative to characteristic speeds of the system. This method was specifically applied to shock identification and classification in 2D MHD compressible turbulence Orszag–Tang vortex evolution. Another example is that of chromospheric detections of intermediate shocks for MHD nonlinearities in sunspot atmospheres,¹⁷ and a technique to identify the high-frequency MHD waves of an oscillating coronal loop¹⁸ has been developed. It is claimed that current detection methods for MHD wave detection remain somewhat limited, and a thorough investigation of MHD wave detection remains largely unexplored.

This study developed two simple methods to detect MHD waves based on convolutional neural network (CNN). The proposed approach was applied to MHD shock refraction at an inclined density interface (e.g., Refs. 15 and 19). Figure 1(a) shows a canonical physical setup to investigate shock refraction. Flow is characterized by incident hydrodynamic shock sonic Mach number M (or fast magnetosonic Mach number for fast-mode MHD shocks); density ratio of the interface $η=ρb/ρ0$⁠; ratio of specific heats γ; angle between incident shock normal and interface α; and nondimensional strength of an initially applied magnetic field $β−1=B2/2p0$⁠, where B and p₀ denote dimensionless magnitudes of the applied magnetic field and gas pressure. Figure 1(b) shows that shock refraction produces a pair of reflected and a pair of transmitted waves. RF and TF are the fast reflected and transmitted magnetosonic shocks, respectively, whereas the wave type of TS and RS is strongly dependent on the chosen parameters (e.g., $β,α$⁠, initial magnetic field orientation whether parallel or perpendicular to the motion of incident shock, etc.).^15,20 Figure 2 shows an example of RS and TS transitions with increasing β (B_x is initially present) by fixing a set of parameters $(M=2,η=3,γ=1.4)$⁠. With increasing β in the strongly planar system (denoted as Ic), RS and TS transit from a slow shock to an $2→4$ intermediate shock and then into a slow-mode compound wave C₁. In the planar system (noted as Ir), with increasing β, RS and TS transit from a slow shock to RD followed by slow shock and then into a RD followed by slow-mode expansion. Further details are provided elsewhere.^15,20 Although Fig. 2 is somewhat difficult to understand, it illustrates that there are many different MHD shock refraction structures with different parameters that are too complicated to identify by simply extending the existing methods used for hydrodynamics.

From the point of view of machine learning (ML) application to fluid mechanics, deep learning has achieved considerable recent success in CFD in recent years.^21–24 For the application of CNN, Lee and You²⁵ have shown that developed CNN could predict vortex dynamics at distinctive flow regimes with flow structures at different scales. Dai et al.²⁶ proposed a CNN-based machine learning model to predict energy conversion efficiency for multi-droplet jumping with different initial distribution angles and radius ratios. The applicability of the machine learning-based reduced order model (ML-ROM) to three-dimensional complex flows has been examined.²⁷ Specially, some shock detection methods based on CNN have been recently proposed for hydrodynamic cases,^28,29 confirming that CNN methods can produce better detection results than some traditional methods. However, deep learning applications for MHD remain limited to date. Kates-Harbeck et al. developed a fusion recurrent neural network (FRNN) to predict disruptive instabilities in controlled fusion plasmas,³⁰ achieving good efficiency to predict disruptions. The current study aims to develop a simple CNN to detect MHD waves, constraining these waves to admissible types for the planar and strongly planar systems.

The outline of this paper is as follows: In Sec. II, we first introduce the fixed output model in which the wave number of each training datum is fixed, followed by the detection model without fixing the number of waves for each sample. Section III presents results and discussions. The two proposed models are first verified via hydrodynamic cases to demonstrate their robustness and efficiency in Sec. III A. The application of MHD wave detection in 1D and the reconstruction of 2D shock refraction are discussed in Secs. III B and III C, respectively. Finally, some conclusions are presented in Sec. IV.

In this section, we present two CNN algorithms to detect MHD waves resulting from the shock refraction process. The two algorithms are developed based on supervised learning, using a training dataset to teach models to yield the desired output. Hence, we need to provide a labeled dataset to train algorithms to accurately classify data and predict outcomes. The labeled dataset consists of input training datasets and labels. Primitive variables $U=(ρ,p,V,B)$ (each such state vector comprises six variables) are employed as the input training dataset, where ρ and p are the density and the gas pressure, respectively, while V and B are two-component velocity vector and the magnetic induction, respectively. Each primitive variable is specified in a training sample, comprising N_t points. Presently, we choose N_t = 100 points in each training sample. Furthermore, each training sample is obtained from a 1D numerical simulation or a straight line containing all five waves $RF,RS,SC,TS$⁠, and TF resulting from the shock refraction process. Here, we only consider these five waves in MHD shock refraction since IS type and unshocked density interface don't vary as the parameters change, e.g., $β,α$⁠, γ, etc. Thus, each training datum has dimension $(6×100)$⁠. The labels of each sample consist of labeled class $c∈C$ and bounds (s₁, s₂) for each wave. It leads to 3w dimension of labels for a sample that consists of w waves. Presently, we consider C classes for seven MHD waves, viz., contact discontinuity (CD), fast shock (FS), slow-mode expansion (Sexp), slow-mode compound wave (Scw), slow shock (SS), rotational discontinuity + slow shock $(RD+SS)$⁠, and $2→4$ intermediate shock $(I24)$⁠. The C classes comprise three waves or hydrodynamics cases: contact discontinuity (CD), shock wave (S), and expansion fans (Exp). We first test and verify the proposed algorithms in hydrodynamics (see Sec. III A), for which an example of training sample is shown in Fig. 3(a).

Figure 3(b) shows 1D CNN architecture used in MHD wave detection. The developed two CNN algorithms have the same main architecture except for the output layer, which is different between the two proposed algorithms. The details of the two CNN algorithms are presented in the Subsections II A and II B. The neural network architecture is composed of the input layer, convolutional layer followed by max-pooling layer (three levels), two fully connected layers, and output layer. The set of some main hyper-parameters of two algorithms are common, e.g., the optimizer “Adam” is chosen since gradient descent algorithms are basic optimizations to attain the optimized value, and the Adam algorithm is one of the most traditional optimizers for its robustness and effectiveness. Another hyper-parameter is the activation function “ReLU” (Rectified Linear Unit) that is set to achieve the nonlinear connection between the fully connected layers. The set of the other hyper-parameters of the two algorithms and their details are presented in the following subsections. The neural networks are implemented using the popular TensorFlow machine learning library.

We now present the first CNN model, where the number of waves w of each training sample is fixed (w = 5 for MHD, and w = 3 for hydrodynamic shock refraction in 1D). The fixed output system models both classification and regression problems. For each wave, the system leads to a $(2+C)$ tensor to predict the wave location and classify its type. The output layer is then encoded as a $w×(2+C)$ tensor. Here, we briefly introduce the hyper-parameter setting for this model:

• Prediction of w wave locations leads to 2w neurons at the output layer, i.e., 10 (6) neurons for MHD (hydrodynamic) wave detection. This is a regression problem, and the activation function connected to the previous layer is set as the ReLU type. The output from this part is evaluated by minimizing the mean square error (MSE) cost function. The metric to monitor the training process is root mean square error (RMSE).

• Classification of w wave types. This is a multi-classification problem (C classes in total) resulting in $w*C$ neurons due to the application of the “softmax” activation function at the final output layer. The number of neurons from the C-classification is 35 for MHD cases, while it is nine for hydrodynamic cases. The traditional cost function for multi-classification problem “categorical-crossentropy” is adopted, and the accuracy is used as a metric.

 Appendix A provides the remaining hyper-parameter settings for this method. The model can be used for some special cases where the number of waves is known, e.g., Sod shock tube and shock refraction problems, but it is not a true detection algorithm since the number of waves is fixed a priori rather than remaining unspecified. However, this model can serve as a verification process due to similarity of its architecture as in the real detection model, which will be presented in Subsection II B.

Recently, several mature algorithms have been developed for real-time object detection in images, including Fast R-CNN,³¹ Mask-R-CNN[R24],³² YOLO,³³ etc. The proposed detection model was inspired by the YOLO algorithm. There are two main differences between the present detection model and the fixed output model. First, we employ divided grids over the input data, leaving the number of waves in the samples unspecified rather than fixing their number. Hence, this algorithm models detection as a regression problem, whereas the fixed output model regards this detection as two problems, namely classification and regression. The main concept is as follows:

• The algorithm divides the input data into N_grid grids. If the center of a wave falls into a grid cell, then this grid cell is responsible for detecting this wave. Each grid cell predicts Bo bounding boxes that predict wave locations, with correlated confidence scores κ. Presently, we set Bo = 1 as this is deemed sufficient for our 1D CNN algorithm.

• The confidence score κ reflects how confident the box contains a wave and also how accurate it considers the box that it predicts. We define $κ=Pr(wave)×IOUpredtrue$⁠, where Pr(wave) is the probability for a wave being present, and $IOUpredtrue$ is the intersection over union (IOU) between the predicted box and the ground truth. κ = 0 if no wave exists in that grid cell and $κ=IOU$ otherwise.

• Each bounding box comprises three predictions: two predictions are the bounds for the predicted bounding box, denoted as (s₁, s₂), and the third prediction is the confidence score κ, representing the IOU between the predicted bounding box and any ground truth box.

• Each grid cell also predicts C conditional class probabilities, $Pr(Classi|wave)$⁠. These probabilities are conditioned on the grid cell containing a wave. We predict one set of class probabilities per grid cell.

• We multiply the conditional class probabilities and the individual box confidence predictions,

  $Pr(Classi|wave)×Pr(wave)×IOUpredtrue=Pr(Classi)×IOUpredtrue,$

  (1)

  which encodes both the probability of that class appearing in the predicted box and how well the box fits the wave. Thus, we consider a prediction to be correct if the classification $Pr(Classi)$ is correct and IOU > 0.5, for example.

This algorithm leads to the final output layer encoded as an $Ngrid×(3 Bo+C)$ tensor. For the chosen default set of parameters, Bo = 1 and N_grid = 10, the number of neurons in the output layer is 100 (60) for the MHD (hydrodynamic) case. Since the problem is modeled as a regression algorithm, the activation function for both the fully connected and output layers is set as ReLU type; loss function and metric are “MSE” and “accuracy,” respectively.  Appendix B details the remaining hyper-parameters for the detection algorithm.

This section presents verifications using the two proposed models for known hydrodynamic cases. In hydrodynamics, training datasets are selected from numerical simulations of 1D shock tube,³⁴ shock contact interaction problems,³⁵ etc. To generate the training dataset, we set output time as 1 and randomly generated the ratio of specific heats, initial ratio of left/right density or pressure, Mach number, etc. Only data with three waves are acceptable for the fixed model, but this is not a limitation for the detection model. Hence, we select the training data with only single shock and isolated contact discontinuity for the second model. Moreover, we use at least the second-order precision numerical scheme or self-similar methods, with more than 200 meshes per unit length. The 5000 sample dataset is divided into training and evaluation sets at 80:20 ratio. In each sample, the number of waves is fixed at three waves, and the results from left to right in the 1D solution profiles comprise a reflected shock or expansion fan R, a contact discontinuity CD, and a transmitted shock T. Detected waves include a combination of S, CD, and Exp. Moreover, the main metrics to analyze the predicted results are accuracy, precision, recall, and F₁ score,

respectively. TP and TN are true positive and negative, respectively, i.e., the number of predictions where the model correctly predicts the positive (negative) class as positive (negative); FP and FN are false positive and negative, respectively, i.e., the number of predictions where the model incorrectly predicts the negative (positive) class as positive (negative).

Figure 4 shows the loss function history resulting from the fixed output model. The loss function history of wave position prediction is almost identical to the total loss function, implying that the main loss function comes from the regression, i.e., prediction of the position for each wave [see Fig. 4(a)]. Figure 4(b) shows that the loss function history resulting from the classification is relatively small compared with the location prediction one, confirming that the fixed output model predicts wave type with very high accuracy [see Fig. 5(a)]. The accuracy of classification approaches unity for all three waves. We also examine loss function and accuracy histories for the detection model in Fig. 5(b). In contrast with the fixed output model, loss and accuracy are shown as their totals since the detection model considers the whole system as a regression problem rather than combined regression and classification problems. Compared with the fixed output model, accuracy of the detection model is not always exactly unity, although it approaches 0.9995 for epoch >30. Relative lower performance for the detection model is attributed to the algorithm that decides if a prediction is correct or not. In the detection model, conditions for deciding whether a prediction is correct or not is monitored by the confidence score within each box and the class probability at each grid cell. Thus, some waves can fail to be detected if the confidence score threshold is set too high, thereby reducing recall. Conversely, some waves may be incorrectly detected if the class probability is set too low, and it leads to reduction of precision. Therefore, these two thresholds are chosen as trade-offs between precision and recall. Presently, we consider each prediction to be correct only if $κ>0.5$ and the probability of class i $Pr(Classi|wave)>0.5$⁠.

After training the models, we employ the two trained models to predict an untrained dataset. Figure 6 shows the confusion matrices for the trained models predicting on an untrained dataset. The fixed model achieves all metrics to unity, whereas the detection model achieves mean accuracy = 0.99, precision = 1, recall = 0.99, and $F1=0.995$⁠. The detection model metrics are only slightly reduced, and hence, the performance of the detection model remains excellent. Thus, we conclude that the detection model is a promising method for MHD wave detection.

This section investigates MHD wave detection for the models verified above using hydrodynamic cases. Data from MHD shock refraction in 1D are employed as the training dataset. We consider seven MHD wave types with five fixed waves in each sample for the fixed output model (RF, RS, CD, TS, and TF), whereas the number of waves in each sample remained unfixed for the detection model. We first examine the performance of the detection model with different values for N_grid, as shown in Fig. 7(a). The evolution of average metrics is not linear as N_grid increases. In the studied range of N_grid, all mean metrics exhibit a local minimum at N_grid = 20, with mean recall, precision, and F₁ being symmetry about N_grid = 20. Hence, the detection model achieves good performance by considering the trade-off between accuracy (0.9435) and F₁ (0.9696) at N_grid = 15. If a high precision >0.99 is expected, N_grid = 10 is also a good option with a relative lower accuracy of 0.9112 and $F1=0.9469$⁠. We next investigate the influence of thresholds of class probability on the model by fixing the confidence as 0.5, see in Fig. 7(b). The model has good performance at several thresholds of class probability, and to retain the stability of the model, we finally choose the threshold of class probability of 0.5. Figure 7(c) shows that there is no significant performance difference for $κ∈(0.2,0.4)$⁠, and F₁ and accuracy reduce as confidence increases above $κ>0.4$⁠. Hence, $κ=0.4$ is a reasonable choice for the confidence threshold. The detection model achieves good performance for MHD wave detection by setting suitable parameters, i.e., N_grid, and thresholds of confidence and class probability. Therefore, we set N_grid = 10, $κ=0.5$⁠, and class probability as 0.5 thresholds for subsequent detection model investigations.

Figure 8(a) confirms that wave location prediction is the main contribution to the total loss function for the fixed output model on MHD cases, consistent with hydrodynamics outcomes. The evolution of the total loss function with increasing epoch is very similar to the regression component (predicting of wave locations). Figure 8(b) shows that classification loss function approximated to $10−7$ for the middle wave (CD) when the model is close to being well trained (epoch >20). This is relatively small loss compared with the other losses since CD is shocked contact for all the training datasets, and hence, it is reasonable that CD loss is the smallest component. In contrast, Fig. 8(c) shows classification loss for RF and TF waves $≈10−2$ when the model is close to being well trained. CD is the fixed type, but each RF and TF is not and could be fast shock (FS) or slow-mode expansion fan $(S exp )$⁠, causing slightly higher loss function for each training epoch compared with CD loss function. Figure 8(d) shows that the same reasoning applies to classification loss evolution for TS and RS waves. TS and RS could be any of SS, I24, S_cw, or RD depending on the initial conditions. TS and RS accuracy was not exactly unity after training due to the prevalence of multiple wave types, whereas accuracy is precisely unity for the other three waves, RF, TF, and CD, as shown in Fig. 9. Figure 10 shows that detection model overall loss tends to be stable and accuracy approximated to 0.999 after training for MHD cases. Thus, the detection model exhibits good performance for MHD cases.

Figures 11(a) and 11(b) show confusion matrices for predicting an untrained dataset using the trained fixed output and detection models, respectively. The fixed model achieves accuracy around 0.999 with similar high recall = 0.998, precision = 0.999, and $F1=0.998$⁠. This high F₁ indicates that the fixed model has excellent robustness for MHD wave detection due to the use of shock refraction in 1D. However, the detection model fails to find several waves, particularly CD type, whereas the fixed model detects CD-type waves with accuracy = 1 and smallest loss. The parameter set is the same as for hydrodynamics: N_grid = 10, confidence $=0.5$⁠, and class probability threshold = 0.5. For the detection model, the accuracy is 0.9113, the recall and precision are 0.9133 and 0.9975, respectively, and the F₁ score is thus 0.9469. The relative low recall may be attributed to N_grid and confidence and class probability threshold settings. Thus, the detection model shows promise for MHD wave detection since its simple algorithm uses only 1D input data.

Finally, we applied the best trained detection model based on 1D input dataset to reconstruct wave structures for 2D MHD shock refraction. The 2D dataset is obtained from numerical simulations using the developed self-similar method^35,36 and subsequently split into numerous 1D datasets by taking vertical or horizontal slices from the 2D data. Each 1D dataset contains R, CD, and T waves for hydrodynamics and RF, RS, CD, TS, and TF waves for MHD. Figures 12 and 13 show reconstructed 2D wave patterns for hydrodynamics and MHD, respectively. The hydrodynamic case fails to detect some waves, including S and CD waves in the triple-point vicinity, and CD waves are erroneously detected as S. Non-detection and incorrect classification are caused by the waves being in close physical proximity, with consequential numerical smearing due to the underlying numerical method. The same reasoning applies for non-detection and incorrect classification near reflected waves from the bottom wall. However, the simple 1D detection model still successfully reproduces the general hydrodynamic shock refraction structure. Non-detected waves and incorrect classification occur in the quintuple-point vicinity for MHD cases, i.e., where all the waves intersect. In this case, RS and TS are SS and I24 waves, respectively, and incorrect wave-type classification occurs for these two types. Similar to the hydrodynamics case, non-detection and incorrect classification in the reflection region from the bottom wall are more severe since the wave patterns are complicated and in physical proximity with each other. All the waves are accurately detected just a small distance away from the triple or quintuple points for hydrodynamics and MHD, respectively. One way to further develop a method is to connect the dots to reproduce the entire wave front for each wave. We successfully reproduce the main 2D shock refraction structure, particularly in vicinity where all the waves interact. Although the proposed model classification is not as accurate as traditional MHD wave classification methods, i.e., Rankine–Hugoniot relations, it offers a reasonable alternative that is somewhat easier to implement and more straightforward to use.

In summary, the present work represents the first instance detecting MHD waves based on a machine learning approach. MHD shock refraction includes many different wave types that interact with each other in a complex manner. Detecting such waves by traditional methods, such as those to detect hydrodynamic shocks, is cumbersome and error prone. The proposed concept comprises neural networks that employ simple 1D training datasets to detect multi-class MHD waves, such as intermediate shocks, compound waves, rotational discontinuities, etc. We developed separate fixed output and detection models. The fixed output model achieved accuracy, recall, and precision all up 0.99. The main shortcoming for this model is that the wave number needs to be fixed a priori for all training datasets. Hence, the model has limited applicability for general wave detection cases. However, the underlying CNN can be adapted to the detection model, which considered wave classification and position prediction as a single regression problem. The essential concept is that we divide the input dataset into N_grid grid cells, which enables the model to be applied to more general wave detection problems with the accompanying advantage that the number of waves need not be fixed. Each grid predicts the confidence score and class probability which are then used as conditions to decide whether the detection is correct or not. The detection model is sensitively affected by the parameter settings, such as N_grid, and confidence score and class probability thresholds, thereby achieving good accuracy approximated to 0.9 and $F1≈0.95$⁠. However, the proposed detection model does not achieve good performance in the vicinity where many waves are in close physical proximity. This was attributed to numerical smearing of these waves, and that the model was trained with the one-dimensional datasets. The proposed CNN methods provide a viable alternative to the traditional Rankine–Hugoniot approach if acceptable accuracy is not extremely high.

This research was supported by funding from King Abdullah University of Science and Technology (KAUST) under Grant No. BAS/1/1349–01-01.

The authors have no conflicts to disclose.

Fang Chen: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Software (equal); Writing – original draft (equal). Ravi Samtaney: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Resources (equal); Software (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.