I solved the following ideal magnetohydrodynamic (MHD) equations numerically in two dimensions to prepare input images to GANs: ∂ ρ ∂ t + ∇ ⋅ ρ v = 0 (1) ∂ ∂ t ρ v + ∇ ⋅ ρ v v + p T I − B B = 0 (2) ∂ B ∂ t + ∇ ⋅ v B − B v = 0 (3) ∂ e ∂ t + ∇ ⋅ e + p T v − B v ⋅ B = 0 (4) p T = p + | B | 2 2 (5) e = p γ − 1 + ρ | v | 2 2 + | B | 2 2 (6) where ρ, p, and v are the density, pressure, and velocity of the gas; B is the magnetic field; γ represents the heat capacity ratio and is equal to 5/3 in this paper; p _T and e represent the total pressure and the internal energy density; I is the unit matrix.

One of typical test problems for MHD, the so-called Orszag-Tang vortex problem (Orszag and Tang, 1979), was solved by the Roe scheme (Roe 1981) with MUSCL [monotonic upstream-centered scheme for conservation laws; (van Leer 1979)]. The initial conditions are summarized in Table 1. B ₀ is a parameter for controlling the magnetic field strength. The compuational domain is 0 ≤ x ≤ 1 and 0 ≤ y ≤ 1. The periodic boundary condition is applied in both x- and y- directions. Simulations for each condition were performed twice on computational grids with different resolutions. The number of grid points is N x × N y = 51 × 51 or 251 × 251 . In the case of N x × N y = 251 × 251 , the calculation time is more than 70 times longer than in the other case though the obtained solution is expected to be close to the true solution.

After the original concept of GANs was proposed by Goodfellow et al. (2014), various GANs have been researched. Among such networks, I focused on pix2pix, which is a type of conditional GAN and a network for learning the relationship between the input and output images. The feasibility of translating from the results of low-resolution grid simulations to those of high-resolution grid simulations has been investigated in this research. Furthermore, in order to enable the translation across two time-series, the architecture combined pix2pix and LSTM has been constructed.

Figure 1 shows the schematic picture of the architecture of the generator in this research. The role of the LSTM layer is to adjust the image translation dependent on the physical time of the simulation; for the initial state of the simulation T = 0 , no translation is needed at all, but as physical time passes, progressively larger translations are needed. Note that the weights of the encoder (decoder) before (after) the LSTM layer are the same in the time direction. Plots of the time evolution of the density in the low-resolution simulations are input to the generator (plots are read as single-channel images). The input images are converted to vectors by the first-half of a U-shaped network (U-Net). In Figure 2, I denote the architecture of the first-half of U-Net in detail. It consists of eight convolutional blocks with a kernel size of 4 × 4 or 2 × 2 . The instance normalization (Ulyanov et al. (2017)) is applied except for the first and last blocks. The activation function is a leaky rectified linear unit [leaky ReLU; Maas et al. (2013)] with a slope of 0.2 for all blocks. A 512-dimensional vector is generated at the end of this architecture.

A series of 512-dimensional vectors converted from the time-series plots is input to the LSTM layer. An input vector x _t originated from the plot at time = t is calculated with the hidden state h _t−1 and memory cell c _t−1. A forget gate f , an input gate i , an output gate o , and part of the term to be added to the memory cell z in Figure 3 are calculated as follows: f = σ W f x t + R f h t − 1 + b f (7) i = σ W i x t + R i h t − 1 + b i (8) o = σ W o x t + R o h t − 1 + b o (9) z = tanh W z x t + R z h t − 1 + b z (10)

where σ is the sigmoid function and tanh is the hyperbolic tangent function; W _⋅ and R _⋅ are the input-to-hidden weight matrices and the recurrent weight matrices; b _⋅ are bias vectors. The hidden state and memory cell are updated by: c t = f ☉ c t − 1 + i ☉ z (11) h t = o ☉ tanh c t (12) The hidden state h _t is reshaped as 1,1,512 . The reshaped hidden state h _t ’ is passed to the latter-half of U-Net and is decoded to the image data (Figure 4). This part consists of eight deconvolutional blocks with an upsampling of the feature map, convolution with a kernel size of 2 × 2 or 4 × 4 (the size of the feature map does not change because the stride of convolution is 1), the instance normalization and activation by ReLU function except the last block. As seen in Figure 1, the generator outputs synthetic time-series plots of density distribution.

The authenticity of the images is judged by the discriminator. Figure 5 shows the details of the architecture of the discriminator in this research. A real image (plot of the density distribution in a high-resolution simulation) or a synthesis image is input to the discriminator. It consists of five convolutional blocks with a kernel size of 4 × 4 . The instance normalization is applied except for the first and last blocks. Except for the last block, the leaky ReLU function with a slope of 0.2 is applied as the activation function. The 16 × 16 patch is eventually output. The discriminator classifies each patch into real or synthetic. We call its architecture the patchGAN (Isola et al., 2017).

The objective of the network is the same as the regular pix2pix as follows: G * = arg min G max D L c G A N G , D + λ L L 1 G (13) L c G A N G , D = E log D x y + E log 1 − D x , G x (14) L L 1 G = E ‖ y − G x ‖ (15) where G and D denote the generator and discriminator, λ is the weighted sum parameter and equal to 100 in this research, and x and y mean the source and target images. G x returns a synthesis image and D x , y or D x , G x returns the probability that y or G x is a real target image. L L 1 G is the mean absolute error (L1 loss) calculated from the pixel-wise comparison between the real image and the synthetic image. The optimizer is Adam with a learning rate of 0.0002.

The architecture is implemented using Keras 2.5.0 and TensorFlow as a backend. The model was trained on Google Colaboratory with Tesla P100-PCIe GPU. For applying the convolution and deconvolution to the sequential data, sets of operations as shown in Figures 2,4,5 are passed to the TimeDistributed layer. The skip connections are implemented by feeding the outputs of the previous upsampling block in the latter-half of U-Net and the same-level (it means that the size of the feature map is the same) convolutional block in the first-half of U-Net to the Concatenate layer. The “return_sequences” and “stateful” parameters in the LSTM layer are set to True and False, respectively.

Figure 6 shows two examples of the time-evolution of density distribution for the training datasets. The top and bottom images of Figures 6A,B show the simulation results, and the middle images are synthesis ones generated from the top ones (the results of low-resolution grid simulations) through the generator. Compared to the high-resolution grid cases, the density distributions in the low-resolution grid cases show less fine structure and become closer to the uniform. Figure 7 displays the comparison of the inhomogeneity of the density between the high-resolution grid cases and the low-resolution grid cases. The inhomogeneity is defined by α = σ ρ / ρ ̄ , where σ _ρ and ρ ̄ are the standard deviation and the average of the density. In the low-resolution grid, the numerical diffusion is larger than in the high-resolution grid, and therefore the inhomogeneity of the density tends to be smaller especially from the middle stage of the vortex development and in the relatively strong magnetic field (Figure 7B). The synthesis images reproduce the fine structures of the density distributions and appear to be well consistent with the high-resolution grid results.

To quantitatively evaluate the quality of the synthesis images, I estimated the density inhomogeneity from the distribution map. When calculating the density inhomogeneity from the simulation result, we can use the value of the density on each grid; however, the density distribution maps (including synthesis images in this research) have only the information of the RGB values. Therefore, to estimate the density inhomogeneity from the distribution map, I trained a three-layer fully connected neural network with 196,608 (256pixel × 256pixel × 3) inputs, two hidden layers of 1,024 and 128 neurons and one output layer. Figure 8 shows the result of the inhomogeneity prediction from the density distribution maps. The horizontal axis is the inhomogeneity calculated from the density values on the grids, and the vertical axis is the inhomogeneity predicted from the distribution maps by the trained neural network. The coefficient of determination R 2 is equal to 0.999. Thus we conclude that the trained neural network provides an accurate estimation of the density inhomogeneity from the distribution maps and the synthesis images.

We can quantitatively evaluate the quality of the synthesis images by inputting those into the neural network and comparing the outputted inhomogeneity with the inhomogeneity calculated from the high-resolution grid simulation results. Figure 9 shows that the inhomogeneity predicted from the synthesis images matches that calculated from the high-resolution grid simulation results with good accuracy; therefore the quality of the synthesis images is definitely good for the training datasets.

In the previous subsection, I have shown that the results for the training datasets are pretty good. However, the generalization ability needs to be investigated for practical use. The testing datasets (the magnetic field strength is different from the training datasets as shown in Table 2) that were not used for training are input to the trained model, and the synthesis images are output from the generator. Figures 10A,B show the comparison of the simulation results and the synthesis images for two example cases. From the 19 cases in the testing datasets, the results for the cases with B ₀ = 0.75 and 1.7 were selected for presentation. The B ₀ = 1.7 case is especially suitable for verifying the generalization ability because there is no training data between B ₀ = 1.5 and 2.0. The top images show the time evolution of the density distribution of low-resolution grid simulations, which are input for the generator; the bottom images show that of high-resolution grid simulations, which are compared with the synthesis images; the middle images are synthesis ones generated through the generator. As with the cases for the training datasets, the synthesis images qualitatively reproduce the density distributions of the high-resolution grid simulations. Even in the B ₀ = 1.7 case, the synthesis images show the fine structure of the density distribution very similar to that in the ground truth images, as shown in the zoomed-in image in Figure 10B.

Figure 11 is almost the same as Figure 9 but for the testing datasets. The density inhomogeneity predicted from the synthesis images through the fully connected neural network (explained in the previous subsection) is in good agreement with the inhomogeneity calculated from the results of high-resolution grid simulations. This result indicates that the method in this research is capable of obtaining high generalization ability.

To demonstrate the effectiveness of the proposed method and the quality of the generated images, I compare the results with those obtained by conventional super-resolution algorithms. The algorithms investigated here are a bicubic interpolation, a Lanczos interpolation, and Laplacian Pyramid Super-Resolution Network [LapSRN; Lai et al. (2017)]. The pixel size of the image to be used as the basis of the super-resolution is 64 × 64, and each algorithm quadruples the pixel size. These results were compared qualitatively and quantitatively with the result of high-resolution grid simulation and the image generated by the proposed method. Plots of the density distribution in high-resolution simulations in the training datasets were used to train LapSRN.

I performed the super-resolution algorithms to the testing datasets (380 images). As an example, the results for the B ₀ = 1.7 and T = 0.38 case are compared in Figure 12. In this case, none of the three conventional super-resolution algorithms can work with a quality comparable to the method proposed in this research. To compare the proposed method with the others quantitatively, the pixel-wise mean squared error (MSE) and the structural similarity index measure [SSIM; Wang et al. (2004)] are calculated between the ground truth image and the synthesis image or the result of super-resolution. Figure 13 shows that the quality of the synthesis images by the proposed method is significantly high compared to that of the results by the conventional super-resolution algorithms.

In this subsection, I discuss an application of this research. As mentioned above, results of high computational cost simulations can be estimated from those of low-cost simulations by the method in this paper. However, it is important to note that simulation results of quite a few cases are needed to train the network ¹ . Therefore, it is not beneficial for a small number of simulations. The more simulations are required, the greater the benefits arise. One such case is optimization based on CFD simulations. As the number of objective variables to be optimized increases, the number of calculations required to obtain the desired performance is expected to increase; in some cases, it takes several thousand cases to evaluate. In such multi-objective optimization simulations, for example, the first dozens to several hundred cases are simulated on both high- and low-resolution grids, and the results are used to train the GANs. After the GANs are trained, low-resolution grid simulations are run, the results are input to the GANs to reproduce the results of high-resolution grid simulations, and objective variables are estimated from synthesis images by, for example, a neural network.

I demonstrate the estimation of computational cost reduction. If the number of simulations required originally and that to train the GANs are N (several thousands in some cases) and N _t (N > N _t ), the calculation times of the high- and low-resolution grid simulations are T _h and T _l (T _h > T _l ), and the computational cost to train the GANs is T _t , the computational cost reduction is roughly equal to N × T h − N t × T h + T t + N × T l (16) where the first term corresponds to the computational cost in the case that all simulations are run on the high-resolution grid, and the second term corresponds to that in the case that the method in this research is applied (the cost to reproduce the results of high-resolution grid simulations by the GANs is negligible compare to performing the simulations). In this way, by substituting low-resolution grid simulations and the result conversion by the GANs for quite a part of high-resolution grid simulations, a great reduction of the computational cost should be achieved.