We present a simulation-based cosmological analysis using a combination of Gaussian and non-Gaussian statistics of the weak lensing mass (convergence) maps from the first three years (Y3) of the Dark Energy Survey (DES). We implement: 1) second and third moments; 2) wavelet phase harmonics; 3) the scattering transform. Our analysis is fully based on simulations, spans a space of seven \(\nu w\)CDM cosmological parameters, and forward models the most relevant sources of systematics inherent in the data: masks, noise variations, clustering of the sources, intrinsic alignments, and shear and redshift calibration. We implement a neural network compression of the summary statistics, and we estimate the parameter posteriors using a simulation-based inference approach. Including and combining different non-Gaussian statistics is a powerful tool that strongly improves constraints over Gaussian statistics (in our case, the second moments); in particular, the Figure of Merit \(\textrm{FoM}(S_8, \Omega_{\textrm{m}})\) is improved by 70 percent (\(\Lambda\)CDM) and 90 percent (\(w\)CDM). When all the summary statistics are combined, we achieve a 2 percent constraint on the amplitude of fluctuations parameter \(S_8 \equiv \sigma_8 (\Omega_{\textrm{m}}/0.3)^{0.5}\), obtaining \(S_8 = 0.794 \pm 0.017\) (\(\Lambda\)CDM) and \(S_8 = 0.817 \pm 0.021\) (\(w\)CDM). The constraints from different statistics are shown to be internally consistent (with a \(p\)-value>0.1 for all combinations of statistics examined). We compare our results to other weak lensing results from the DES Y3 data, finding good consistency; we also compare with results from external datasets, such as \planck{} constraints from the Cosmic Microwave Background, finding statistical agreement, with discrepancies no greater than \(<2.2\sigma\).