Deep learning for synthetic microstructure generation in a materials-by-design framework for heterogeneous energetic materials

The sensitivity of heterogeneous energetic (HE) materials (propellants, explosives, and pyrotechnics) is critically dependent on their microstructure. Initiation of chemical reactions occurs at hot spots due to energy localization at sites of porosities and other defects. Emerging multi-scale predictive models of HE response to loads account for the physics at the meso-scale, i.e. at the scale of statistically representative clusters of particles and other features in the microstructure. Meso-scale physics is infused in machine-learned closure models informed by resolved meso-scale simulations. Since microstructures are stochastic, ensembles of meso-scale simulations are required to quantify hot spot ignition and growth and to develop models for microstructure-dependent energy deposition rates. We propose utilizing generative adversarial networks (GAN) to spawn ensembles of synthetic heterogeneous energetic material microstructures. The method generates qualitatively and quantitatively realistic microstructures by learning from images of HE microstructures. We show that the proposed GAN method also permits the generation of new morphologies, where the porosity distribution can be controlled and spatially manipulated. Such control paves the way for the design of novel microstructures to engineer HE materials for targeted performance in a materials-by-design framework.

Scientific RepoRtS | (2020) 10:13307 | https://doi.org/10.1038/s41598-020-70149-0 www.nature.com/scientificreports/ in this paper is directed towards establishing a computational approach to relate stochastic microstructures of energetic materials to their observed performance. To achieve this overarching goal, in silico experiments on an ensemble of stochastic microstructures have been performed in the previous work 23,24 , to extract quantitative data from resolved meso-scale simulations in the form of surrogate models 25,26 . These models are used to close the macro-scale governing equations in a multi-scale predictive framework. The work in this paper presents a deep learning approach to generate ensembles of synthetic microstructures that can be used for simulations; the methodology also allows for designing new microstructures, paving the way for materials-by-design of energetic materials 13,27 . Simulations on an ensemble of microstructures also facilitate uncertainty quantification (UQ) due to microstructural variability 24 i.e. aleatory uncertainties associated with the inherently stochastic microstructure. The propagation of uncertainty across scales influences the overall prediction uncertainty of HE sensitivity at the macro-scale. Therefore, ensemble simulations performed on a sufficiently large set of synthetic microstructures will be a key enabling tool for the reliability-based design [28][29][30] of next-generation HEs.
need for ensemble simulations on synthetic microstructures. In recent years, multi-scale, multiphysics models that establish structure-property linkages have begun to be developed for the response of HEs to shock and impact loading 31 . A key task in the multi-scale modeling workflow is the performance of ensembles of high resolution, high-fidelity simulations of the meso-scale reactive mechanics 23 . The goal of such simulations is to capture the essential physical ingredients and to quantify energy localization at micro-scale morphological features such as pores and crystal-crystal interfaces. A predominant mechanism for energy localization in the microstructure is the creation of hot spots which are formed due to the collapse of pores 32 . Pore collapse is a well-studied problem, both theoretically [33][34][35][36] and experimentally 37 . Reactive simulations of the dynamics of pore collapse in microstructures extracted from images taken using scanning electron microscopy (SEM), X-ray computed tomography (XCT), etc have also been performed 23 . However, since the microstructure is stochastic, a large ensemble of meso-scale calculations, using a large enough number of statistically representative microstructures is needed to extract statistically meaningful information from such simulations. However, computational modelers typically do not have access to large sets of imaged data for a variety of types of energetic materials, or even for various micro-morphologies of a single type of energetic material. Image acquisition is expensive, and imaged data may be distribution-sensitive and specific to limited formulations and material types. From a computational mechanics standpoint, it is extremely useful to have access to an ensemble of microstructural geometries so that simulations can be performed and the statistics of micro-morphologies can be correlated with measures of sensitivity. In the absence of a large database of microstructural images, computational scientists rely on generated synthetic microstructures 19,[38][39][40][41][42][43] (hereinafter abbreviated as "SynµS") as proxies to the real microstructures (hereinafter abbreviated as "RealµS"); a large array of stochastic Synµ S that closely mimic the real sample must be created and used in in silico experiments. Generation of realistic Synµ S, with adequate statistical representation and ability to control global and local features in a versatile and flexible computational framework will prove to be a key component in the materials-by-design process for precise and controlled performance of energetic materials. previous approaches. Generating an ensemble of stochastic Synµ S that stand in for the Realµ S is a challenging task. In general, Synµ S can be generated using several different approaches 42,44 . Of these, approaches based on shape descriptors have been used in the past; objects are inserted into a computational domain using shape packing algorithms constrained by global shape descriptors, such as volume fractions and particle size distributions. While packing algorithms 44 can be used to generate microstructures with specified morphometric characteristics (e.g., porosity, particle size distributions), these methods are limited to regular/analytical shapes such as spheres, ellipsoids, and polygons 45 . It is also difficult to pack regular shapes for theoretical maximum densities (TMDs) significantly higher than the close-packing limit; typical pressed energetic materials of the type simulated in the current work can have TMD values greater than 90% 46 . Another popular approach, particularly to achieve high TMDs, is to start with space-filling polygons via tessellation, followed by superposition of voids or other phases in the mixture 47 . In both packing and tessellation approaches, however, the microstructure fails to mimic the morphology of real samples for a broad class of heterogeneous materials; i.e. flexibility as well as realism are lacking. In the particular context of energetic materials, both shape-packing and tessellation-based approaches 19,47 have been used to perform meso-scale simulations of the shock response of stochastic microstructures. However, the Synµ S generated by these previous approaches are still too "ideal", in that the range of possible void/defect/interface shapes are not well represented in these two approaches. This shortcoming is significant for two reasons: 1. Realµ S possess features that contain non-ideal distributions of shape features, including outlier features such as large cracks and elongated, tortuous void structures. It has been shown experimentally that such structures possessing large surface-to-volume ratios and other shape characteristics play a significant role in energy localization and sensitivity 48 . 2. Computational studies 34,35,49 have shown that local features in microstructures, such as inter-void distances, void shapes and orientations are in fact key aspects of sensitivity that distinguish different classes of the same material. Therefore, Synµ S must mimic not only global features but also local structural characteristics to be useful in building predictive models of material performance and materials-by-design frameworks.
Machine learning approaches are promising alternatives to overcome the limitations of shape-descriptor based approaches. In recent years, convolutional neural networks 50  www.nature.com/scientificreports/ proposed a method for Synµ S generation based on a general-purpose texture synthesis method in computer vision that uses transfer learning 52 . For a given Realµ S input, their method generates a Synµ S having the same "style" as the input Realµ S by minimizing the style difference, where the style of a microstructure is defined by the Gram matrix of feature maps produced by the CNN. However, the transfer learning method requires a Realµ S as a reference to generate Synµ S. Hence, from the materials-by-design standpoint, material morphology can be explored only in the "neighborhood" of the existing Realµ S samples. This poses a critical limitation in generating a large ensemble of microstructures that can span the space of candidate material morphologies. In another deep learning approach, Cang et al. 53 and Guo et al. 54 employed an encoder-decoder architecture to generate Synµ S. The encoder-decoder architecture develops a codified representation of micro-morphology by learning to compress the image pixels (encoder) and reconstruct it back to the original one (decoder). The code values learned by encoder-decoder networks parameterize micro-morphology, allowing the generation and manipulation of Synµ S by "turning knobs, " where the code values act as the knob control parameters. However, these networks tend to generate blurry images. On the other hand, another deep learning-based method, the patch-based generative adversarial networks (GAN), can generate much sharper and crisper images as we will demonstrate later. Unlike the original GAN method 55 , patch-based GAN approaches evaluates the quality of generated images at different local regions, enforcing the details to be clearer and more realistic. Yang et al. 56 and Mosser et al. 57 demonstrated that GANs are not only capable of generating realistic Synµ S but can also be used to continuously parameterize the micromorphology; this paves the way for smoothly varying the morphology to produce new microstructures. Fokina et al. 58 also confirmed the good performance of GANs in generating Synµ S of ALPORAS aluminum foam. However, in these previous works, the stochastic variations in micro-morphology were not investigated or quantified. In addition, the output microstructure was fixed at a certain size and was not scalable to arbitrary sizes. Furthermore, the above-mentioned methods lacked the capacity to control micro-morphology at different local regions to produce spatially varying morphologies in a single Synµ S sample. In this paper, we develop a flexible and versatile algorithm for generating realistic microstructures using GAN. The new algorithm allows control of micromorphology in different regions and can be scaled to arbitrary image dimensions seamlessly.
A method for microstructure generation using GAn. Here, we employ a patch-based, fully convolutional GAN architecture as illustrated in Fig. 1. The generator takes two input vectors ρ ∈ R r and ∈ R l defined at each of h × w grid locations, forming an h × w × (r + l) input tensor. The input tensor is then up-convolved five times, each time scaling the dimension by a factor of 2, to produce an H × W microstructure image. The generator is trained together with a detector network symmetric to the generator, in which the upconvolution layers are replaced by regular convolution layers of the same kernel size and stride. During training, the detector network is presented with an arbitrarily chosen image, either Realµ S or Synµ S, and is tasked to determine if the presented image is real or synthetic. The determination of whether a presented image is real or synthetic occurs at each of the h × w grid locations of the generator. The patch-wise feedback on the image quality promotes details to be more closely captured 59 . Furthermore, since the receptive fields of the detector network overlap by 32 pixels between adjacent patches, smooth and seamless connection between the patches is naturally enforced. Finally, it is worthwhile to note that the proposed architecture is fully convolutional so that arbitrary-sized images can be produced without stitching 58 by varying the size of the input tensor. www.nature.com/scientificreports/ Two input vectors shown in Fig. 1-the local stochasticity parameters ρ and the global morphology parameters -play a critical role in the proposed GAN model. The role of the local stochasticity parameters ρ is the same as the "noise tensors" in the standard GAN implementation; they serve as seeds for adding stochastic variations. The global morphology parameters , on the other hand, control the overall morphological characteristics of the generated image, such as grain (or void) sizes, orientations, and aspect ratios. We achieve such a control of the global morphology by setting to be constant across different grid locations during training, while ρ varies randomly across grid locations. Both ρ and are uniformly distributed in the range [−1, 1] and hence parameterize the morphological variations of a material in the domain [−1, 1] (l+r) . In the following section, we demonstrate the fidelity and versatility of the proposed GAN approach for Synµ S generation in comparison with the current state-of-the-art transfer learning (TL) approach presented Li et al. 51 .

Results
Qualitative and quantitative comparison of real and synthetic microstructures. When mechanical pressing techniques are used to produce HE materials, defects such as voids, cracks, and inclusions are created 60 . Figure 2 (Top row labeled 'ground truth') shows examples of one such microstructure 61 -subsampled from an image obtained using SEM-illustrating the distribution of features. In this type of pressed HE, both inter-crystal and intra-crystal voids are nearly uniformly presented within the microstructure, i.e. crystals and voids do not appear to vary significantly in size and concentration within each sampled image nor across the different samples. As shown in previous work 23 , when the pressed HE material is subject to shock loading in the range of 10-20 GPa, voids in the microstructure collapse, leading to localization of energy and formation of high temperature hot spots. The ignition and subsequent growth of these high temperature hot spots depends on the physicochemical properties of the crystalline material and on the morphology of the void and crystal distributions in the microstructure 33,48,62 . In the current paper, to produce Synµ S that mimic the behavior of the real (i.e. imaged) ones, we focus on comparing the distributions of key morphological metrics of void and crystal phases and the shock response of the Realµ S and SynµS.
The Realµ S (top row), Synµ S generated by the TL method of Li et al. 51 (middle row), and Synµ S generated by the proposed GAN method (bottom row) are presented in Fig. 2 for qualitative comparison. Overall, the GAN generated Synµ S are visually more similar to the Realµ S, whereas artifacts and blurry crystal boundaries are observed in TL generated Synµ S as highlighted in closed-up images. In addition, complex geometry of natural voids and crystals in Realµ S are well emulated in the GAN generated Synµ S making them almost indistinguishable from the Realµ S. Quantitatively, we compare the statistical distributions of void shape descriptors such as void diameter D void , void aspect ratio AR, and void orientation θ , which are shown in Fig. 3. From each of the categories (i.e., Realµ S, transfer learned (TL) Synµ S, and GAN-SynµS), 25 random samples were drawn. Morphometric analysis was performed using methods described in Roy et al. 31 . The distribution of the shape descriptors of the Realµ S and Synµ S were generally in good agreement, both between the GAN-Synµ S and Realµ S as well as between the TL-Synµ S and Realµ S. The whiskers in the figures indicate the standard deviation among the 25 microstructures while the curves correspond to the means. In general, the generated Synµ S have www.nature.com/scientificreports/ voids of sizes in the same range as the Realµ S; the peak of the distribution is shifted slightly in the generated microstructures. As demonstrated below, the observed differences in the size distribution have negligible effects on the computed quantity of interest, viz. the hot spot ignition and growth rate. The void aspect ratio distribution of the real and imaged microstructures is likewise in good agreement, albeit with a shift in the peak of the distribution. On the other hand, the void orientation θ distribution plots show good agreement; a small peak is observed at θ = 45 • but the overall distribution of void orientations is fairly uniform, i.e. there is no strong orientational preference of the voids. As seen from the figure, transfer learning also shows overall good agreement with the Realµ S void distribution. For porous energetic materials, the sensitivity is strongly dependent on the characteristics of the void field in the microstructure. The void fields are interstices between packed crystals. Both the GAN and TL approaches produce crystals in the microstructure that are in good agreement with those in the real microstructures. To demonstrate this, statistical distribution of crystal shape descriptors have been obtained and are shown in the online supplementary materials ( Supplementary Fig. 4).
In addition to the above shape descriptors, two-point correlation functions 63 were obtained to quantify the void phase morphology. The computed two-dimensional two-point correlation functions for the Realµ S and Synµ S are displayed in Fig. 4 where ϕ(r) indicates the volume fraction of void phase at distance r. The result shows that all three microstructures have similar volume fraction of the void phase. However, unlike the transfer learning approach, the GAN generated microstructure shows a weak correlation that persists across the plot. Nevertheless, the magnitude of the correlation is small, and does not exert any discernable influence on the hot spot dynamics simulated below. Furthermore, the horizontal cross-sections of ϕ 2 (r) extracted from the center of the two-dimensional plots show that both the Synµ S are within the standard deviation of the Realµ S. In addition, both plots stabilize to ϕ 2 at r = 1µ m, showing strong agreement of the correlation length with the Realµ S. Therefore, the two-point correlation function plots reveal that the generated microstructures closely resemble the statistics of the void phase in the Realµ S. Figure 4 indicates that the TL approach produces a 2-point correlation in better agreement with the real image that does GAN. However, we remark that the performance of TL comes at the expense of solving an optimization problem for each Synµ S generation, as TL generates a Synµ S by minimizing the style difference of the generated Synµ S from a reference Realµ S input. Furthermore, it should also be noted that a TL-Synµ S tends to hew close to the immediate neighborhood of the reference Realµ S in the morphology space due to the style-difference minimization scheme. Whereas, the proposed GAN method can explore the entire morphology space by parameterizing the space by varying the input parameter , as will be demonstrated below. On the other hand, it can be argued that TL-Synµ S display slightly better agreement with the Realµ S in morphometry due to the style difference minimization by which the TL method generates images: the TL approach generates Synµ S only within a narrow margin of style difference and, therefore, TL-Synµ S are biased towards the Realµ S samples; whereas the proposed GAN method generates a broader range of novel  www.nature.com/scientificreports/ microstructures. Hence, for the challenge of generating a large ensemble of Synµ S from a small set of Realµ S images, the proposed GAN method offers greater practical benefits than the TL approach. Aside from the lack of artifacts and blurring in the GAN-Synµ S compared to the TL-Synµ S samples as shown in Fig. 2, the capacity to generate a broader range of novel microstructures using only a small number of Realµ S samples in the case of GAN is a meaningful advantage given the typical paucity of available microstructure images for a general HE. The GAN approach therefore provides greater flexibility and versatility in the generated microstructures and is better aligned with our goal of creating a pathway to materials-by-design.
Simulations of reactive dynamics in real and synthetic microstructures. The void microstructure of a pressed energetic material strongly affects the meso-scale reactive dynamics of the material under shock loading. Both experiments 6,48,64 and simulations 33,34,65 have shown that crystal and void size distributions, void volume fractions, and void morphology affect the meso-scale sensitivity of the material. Here, reactive computations are used to compare the dynamics of the void collapse process and subsequent chemical decomposition of the solid HMX material for real and GAN generated microstructures. The methods for computation of void collapse process and reaction initiation from imaged microstructures are described in several previous publications 23,36,65 and only briefly outlined in the methods section of this paper. The simulations are performed by applying shock loading from the left end of the domain at a shock pressure P s = 9.5 GPa. The size of the Realµ S and GAN-generated Synµ S are both 25 µm × 25 µm . The sample is placed in the middle of the computational domain to avoid edge effects from the domain boundaries. Extra padding regions are provided surrounding the imaged sample to account for the translation of the material during the shock passage through the material. In Fig. 5, snapshots of the evolving temperature and pressure fields are shown at several instants of time. The temperature field measures the intensity of hot spots that resulted from the process of void collapse and the pressure plot shows the blast waves that emanate from the void collapse events. As seen in the figure the temperature and pressure fields generated due to shock passage are nearly identical for the Realµ S and GAN-Synµ S. Thus, the dynamics of the void collapse and hot spot evolution is well represented by the GAN-generated microstructure. A more quantitative assessment of the physically realistic response of the Synµ S is obtained from Fig. 6, which shows the time evolution of a quantity of interest in meso-scale simulations, which is the total area occupied by hot spots in the domain, calculated using the approach described in the methods section. As seen in the figure the total hot spot area A HS calculated for the Realµ S and Synµ S cases is in good agreement over the course of the evolution of the hot spots. In summary, from Fig. 5 the temperature and pressure fields as well as the quantity of interest (QoI) obtained from simulations using the synthetic microstructure captures the dynamics of real, imaged microstructures with good fidelity. Therefore the GAN-Synµ S can serve as suitable proxies in ensemble simulations to extract the physics of this type of pressed HMX material. towards materials-by-design: controlling the micro-morphology. Controlling morphological features by changing and ρ. As discussed above, the GAN model encodes global micro-morphology characteristics into a parameter vector . The numerical values of the elements of the vector are the "knobs, " which    Fig. 8 illustrates subtle variations in micro-morphology for a specific value of the morphology parameter , obtained by changing the values of ρ . It is noticeable that the overall morphology looks similar across the columns, while the layout of crystals/voids and the subtle details are different. Note that significant variations in the void and crystal characteristics are achieved by varying the (comparing across rows), while subtle variation, particularly stochasticity, is achieved by varying ρ (comparing across columns).
Using to create "in between" and spatially varying microstructures. Another useful property of the GAN morphology parameter is that changes in lead to smooth and linearly proportional changes in micro-morphology, as demonstrated in Fig. 9. In the case of GAN, to span the space of morphological parameters (void/crystal size, shape, and orientation), the vector was linearly interpolated between two book-ending images to produce a smoothly varying range of microstructures. In the case of TL, there was no trivial method for linearly interpolating the microstructure, but we attempted to modify the TL algorithm by making it possible to interpolate the Gram matrices (i.e. the representation of "styles"). Figure 9 shows that an interpolation between two different microstructures can be achieved by both TL and GAN approaches. However, the TL method was not able to generate a smooth, continuously varying interpolation, while the proposed GAN method demonstrated the opposite: e.g. grain sizes GAN-Synµ S in Fig. 9 grow almost linearly from α = 0 to α = 1 . This is a clear advantage compared to the current state-of-the-art machine learning based Synµ S generation approaches, where such smooth and continuous morphology changes cannot be easily accomplished. The parameter linearly maps the space of micro-morphology, so that continuously varying from one value (say 0 ) to another (say 1 ) produces smoothly and linearly varying micro-morphologies.
The ability to generate "in-between" morphologies using the GAN approach is of key importance for the realization of a materials-by-design framework. For example, it is well-known that the sensitivity of pressed HE samples is affected by the crystal (and concomitant void) sizes 48 . However, in previous numerical simulations of the meso-scale shock response 23 , it has not been possible to thoroughly quantify the behavior of a range of real, i.e. imaged microstructures of a specific material nor the change in sensitivity due to variabilities in the microstructure; this is due to the lack of sufficiently large numbers and types of imaged samples and the lack of ability to control micro-morphology in image-based shock simulations. Using the present capability, it becomes possible to conduct shock simulations with ensembles of stochastic, morphology-controlled Synµ S. Figure 10 shows results from shock simulations conducted on two samples of Synµ S with significantly different crystal/ void sizes; in Fig. 10a the temperature field for a sample with small crystal/void sizes is shown on the left and for a sample with larger crystal/void sizes is shown on the right. In Fig. 10a, while the difference in the mean crystal size is rather modest (around 1 µ m for the left column and 2.5 µ m for the right column), the simulation results show significant differences in the hot spot development and therefore in the overall material sensitivity. The sample with smaller crystals/voids (left column) shows a higher density of hot spots in the control volume when compared to the sample with larger crystals/voids (right column). Figure 10b shows that there is a noticeable difference also in the pressure field between the two microstructures; similar to the hot spot temperature field in Fig. 10a, it is observed that the pressure field in the case of the small crystal size sample (left column) achieves higher overall magnitude in a larger part of the domain when compared to the pressure field for the sample with larger crystal (right column). The relatively higher sensitivity of the sample with small crystals/voids  www.nature.com/scientificreports/ is quantitatively depicted in Fig. 11 which shows the evolution of the hot spot area with time. Clearly, the small crystal microstructure has a larger hot spot growth rate compared to the large crystal microstructure. In addition to its ability to create "in between" microstructures by smoothly interpolating between two microstructures that define a range of possible morphologies, the proposed GAN method can exercise spatial control of the micro-morphology, providing the ability to spatially grade morphologies within a sample. This is achieved by spatially varying the morphology parameter . Although is constrained to be constant across the spatial grid during training, at the inference time, different locations can have different values of . Through this, for example, one can generate a smoothly graded Synµ S by allowing to vary continuously from one location to the other, as in Fig. 12a. Here, we show that the generated Synµ S has crystal and void sizes that vary spatially, with small voids/ crystals on the left end of the sample transitioning smoothly to larger voids/crystals on the right end. This type of spatial variation in the microstructure can allow for more precise control on the rate of burning of the energetic material. Furthermore, one can "paint" micro-morphologies by placing different "shades" of , as demonstrated in Fig. 12b, where a "Hawkeye" pattern shown on the left is imprinted into the crystal size distribution, resulting in a microstructure on the right with a morphology that reflects this pattern. In both examples, it is noteworthy that there are no observable stitching boundaries or awkward image transitions. Using this facility, combined Figure 11. Hot spot area evolution with respect to time. It is noticeable that the crystal size affects the growth rate of the hot spot area. Figure 12. (a) A graded microstructure generated by linearly interpolating two global morphology parameters. A global morphology parameter 0 corresponding to small crystal sizes and a parameter 1 corresponding to large crystals were obtained manually. The image was then generated by spatially grading 0 and 1 linearly: www.nature.com/scientificreports/ with the on-going micro-scale additive manufacturing (or micro 3D printing) techniques 13,21,27,66 , the present controlled micro-morphology generation technique can lead to the realization of microstructure-engineered materials, which is being pursued by the authors in ongoing work.

Discussion
A novel GAN model for Synµ S generation has been developed in this work. The method is shown to produce microstructures that qualitatively and quantitatively replicate real microstructures obtained from images. The approach also provides control over the generated morphology, including the ability to spatially vary the morphological properties seamlessly. Compared to the recent deep learning based material synthesis methods 51,53,54,[56][57][58] , the advantage of the present method is the ability to scale to arbitrary size without stitching or quilting, to produce a linear and continuous control of morphology, and an ability to generate a microstructure with spatially controlled morphology. There are several avenues for further development of the current approach. For example, the current morphology parameter is at a single scale, and cannot yet capture a large variation of crystal sizes, which occurs frequently in natural images. As described in depth in the methods section, an input grid location of the generator controls 125 × 125 pixel area of the generated image. We found that the GAN was not able to generate crystals bigger than 125 pixels in size, which currently limits our capacity to cover the full range of morphometry. To this end, multiscale GAN architectures recently proposed in GAN literature (e.g., Karras et al. 67 ), in which the input parameters are injected at different scales of convolution are being explored. Another interesting direction of study being currently pursued is to unravel the intuitive, semantic meaning of the global morphology parameters . This can potentially be achieved via correlation analyses between each dimension of the morphology parameters and the conventional morphometric measures, such as void/crystal diameters, aspect ratios, etc. Our preliminary study, in which we conducted a multivariate linear regression to test the correlation between GAN parameter and morphometric measures (Supplementary Figs. 5, 6 and Supplementary Tables 1, 2), revealed statistically significant correlations between morphometry parameters and GAN parameters. Hence, a further investigation may allow the GAN be reparameterized with morphometry parameters, permitting more intuitive handles to the users. It is also remains to establish the optimal dimension of . Unfortunately, theoretical works on this issue are lacking in the literature, other than naïve search methods for finding the minimal loss by varying the dimensions. Therefore, more fundamental research should be directed on this avenue.
With regard to design and analysis of energetic materials, an immediate extension of this work would be to use the GAN method to extensively study the sensitivity of heterogeneous energetic materials to applied loads; for this purpose we will characterize the effects of micromorphology on different physical QoIs, such as hot spot ignition and growth rates 26 . This can be accomplished by performing computational simulations on a large ensemble of morphology-controlled microstructures, as briefly demonstrated in Figs. 10 and 11. This will lead to a comprehensive quantitative characterization of the detonation behavior of HE materials as a function of micro-morphology, and will provide surrogate models for bridging meso-scale behaviors to macro-scale simulations 16 . Further development will lead to a materials-by-design framework where a desired performance or property of a HE material can be engineered through optimization of material morphology, analogous to topology optimization in mechanical component design. The ability of the proposed GAN method to spatially control micro-morphology will play a key role in the realization of such a framework.

Methods
Data. We used a raw image of a class V cyclotetramethylene-tetranitramine (HMX) pressed energetic material 60,61 obtained using SEM for the training of the GAN model. The original image was 3, 000 × 3, 000 in size, from which we sampled 161 × 161 image patches by cropping the image at a random position. The physical resolution of the image is 52 nm per pixel, or 19 pixels per 1 µ m. The 161 × 161 sample patches, therefore, correspond to 8.5 µm × 8.5 µm in physical dimension. For training, a total of 12,500 of such samples were generated and the pixel values were re-scaled from [0, 255] to [− 1, 1].
GAn. GAN 55 is a generative model capitalizes an adversarial training of two or more competing neural networks. Typically, an image generator network G and a real/synthetic discriminator D are trained together, which poses a minimax game where G tries to maximize the misdetection rate while D tries to minimize it: where P data is the distribution of real images and z is the input parameter to G . The convergence is reached at the equilibrium, where the discriminator D is no longer able to distinguish synthetic images generated by G from real images drawn from the data set.
In our work, the GAN loss in Eq. (1) is evaluated multiple times at grid patches and averaged over grid locations, resulting the following loss function: where D j,i is the prediction of D at jth row and i-th column on the h × w grid in Fig. 1.
The input tensor z ∈ R h×w×(r+l) is composed of the local stochasticity parameter ρ j,i ∈ R r and the global morphology parameter vector j,i ∈ R l defined at each grid location (j, i). During the training, we set the global www.nature.com/scientificreports/ morphology parameter j,i = for all i, j for some constant randomly drawn from the uniform distribution in [− 1, 1], while the local stochasticity parameter ρ j,i is independently drawn from the uniform random distribution in [− 1, 1]. At the inference time, j,i can vary across the grid, as demonstrated in Fig. 12 and controlled by the user. The architectures of the generator and the discriminator are symmetric. The discriminator contains a stack of five convolution layers of kernel size 5 × 5 and stride 2. The generator has the mirrored architecture, where the convolutions are replaced by upconvolutions of the same kernel size and the stride 1/2. Such architecture results in a 125 × 125 receptive field, which determines the window size that one grid parameter j,i controls.
For the training, the adaptive moment estimation (ADAM) optimizer with the learning rate of 0.0002 and the batch size of 10 was used. We set r = 30 and l = 15 as the dimensions of ρ and , respectively. During the training, h = w = 5 was used.
Morphometry. Morphometry refers to the quantification of the microstructural features using shape descriptors or correlation functions. The shape descriptors calculated in this paper are used to quantify the size, shape and orientation of voids and crystals, the two primary features in the microstructure of pressed HMX. In the present work, the shape descriptors are obtained using level-set based morphometry techniques described in detail in Roy et al. 31 Briefly, the gray-scale microstructural images are segmented in the image (pixel) space using active contouring to obtain level-set fields 68 . With the zero level set value defining the crystal shapes, various approaches have been developed in Roy et al. 31 to calculate the sizes of voids and crystals, their aspect ratios and their orientations. The orientation of voids ( θ ) are defined as the mean of the orientations of the voids relative to the direction of the shock propagation (from left to right across this paper). It is calculated by averaging over the local orientations within each void, where the local orientation is calculated at each Cartesian grid point in the interior of the voids using the volume orientation method 69 . These three quantities are then computed for each void/crystal and histograms (pdfs) are developed to obtain the distribution of void sizes, aspect ratios and orientations of all voids/crystals in the domain. The second approach to morphometry employed here is the use of correlation functions; n-point correlation functions provide rich information on the structures embedded in the domain. Typically, n < 3 in the interest of computational cost. Two-point correlation functions 63 determine the probability of finding 2 random points with position vectors p and p + r in a given phase. Geometrically, the two-point correlation function S 2 (r) can be interpreted as the probability of having both ends of the chord vector r in the same phase. The two-dimensional two-point correlation function S 2 (r) is calculated as the following probability function: S 2 (r) = PI(p) = 1, I(p + r) = 1 where, I(p) and I(p + r) are the indicator function values at locations p and p + r respectively, r is the magnitude of the length of the chord vector in the normal direction r/ r . The indicator function I(p) is used to identify phase at a particular point in the image space. In this work since the void phase is of interest, I(p) is equal to 1 if the position vector p indicates void phase and 0 otherwise. In the present work, the two-point correlations are obtained using the open source software developed by de Geus et al. 70 .
Meso-scale simulations of reactive mechanics. Meso-scale simulations are performed to simulate the collapse of voids in the microstructure due to the passage of a shock wave. The shock is parameterized by the pressure P s . Reactive calculations are performed using methods discussed extensively in previous publications 65 . Several validation exercises have been demonstrated in previous work, providing high confidence in the physical correctness of the shock computations 36 . By performing the meso-scale simulations, temperature, pressure and species field data are utilized to quantify the response of the pressed material to the imposed shock, as shown in previous work 23 . In the present context, the QoIs used to quantify the effect of microstructure on the sensitivity of the material, are calculated by following the evolution of the temperature field and the reaction product mass fraction in the sample. The temperature field T(x, t) in the domain measures the intensity of a hot spot resulted from the process of void collapse. Higher temperature hot spots formed due the collapse of voids in the material lead to higher chemical decomposition rates. The reaction zone defines the hot spot in the domain, which is defined as the region where the temperature of the material exceeds the value of the temperature ( T bulk )reached after the passage of a planar shock wave. The hot spot area A hs is a significant quantity of interest for determining sensitivity and is calculated as the area of the domain where the temperature T(x, t) > T bulk . The hot spot area is recorded throughout the simulation to track the evolution of the hot spots.
Scientific RepoRtS | (2020) 10:13307 | https://doi.org/10.1038/s41598-020-70149-0 www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creat iveco mmons .org/licen ses/by/4.0/.