  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$.
