We present a simulation tool, TETHYS (i.e., Two-dimensional Emitter of THz, Hydrodynamic Simulation), developed with the aim of studying hydrodynamic models for electronic transport in graphene, being particularly tailored for the simulation of the configuration of top-gated graphene field-effect transistors (GFET) near room temperature. We present the model governing the system and briefly discuss its validity and physical parameters; then we proceed to discuss the algorithms and numerical methods used in its implementation, namely, two-step Richtmyer and weighted forward-time centred-space schemes, which address the hyperbolic and parabolic parts of the model, respectively. An abridged usage guide and description of the output are also provided, as well as two illustrative simulation cases: the Dyakonov–Shur instability and the boundary layer of pulsatile electric flow. The simulation program here brought forward is uncomplicated and lightweight, yet robust and accurate, having a wide variety of applicable scenarios and configurations to explore. Program summaryProgram Title: TETHYSCPC Library link to program files:https://doi.org/10.17632/s28zwgpnpw.1Developer's repository link:github.com/pcosme/TETHYS-Graphene-Hydrodynamic-SimulationLicensing provisions: MITProgramming language: C++Nature of problem: The main objective of this code is to simulate and analyse the electronic conduction on gated (mono-layer) graphene resorting to a hydrodynamic model [1]. Note that, given the Dirac-like dispersion relation of such material, the mass of the carriers is ill-defined; and, in fact, the mass of the fluid element varies with the number density, which sets this system apart from other 2D electron fluids [2]. The model assumes enough doping, so that only one kind of carriers, electrons or holes, participates in the dynamics of the system. Moreover, the simulation is limited to the fully degenerate limit.Thus, these assumptions correspond to the situation of a regular GFET and this application can be used to explore atypical phenomena arising from the effects of uniform and static magnetic field and Hall viscosity, and even spatially non-uniform gate capacitance.Furthermore, this software was designed with the intent of examining unstable nonlinear modes that may lead to the emission of radiation in the THz range [3].Solution method: The fluid equations (continuity, momentum and energy transport) are solved with a split-operator technique, with a second-order finite-volume Lax–Wendroff method [4] for the hyperbolic part and a forward-time centred-space of second-order stencil to tackle the dissipation and heat diffusion. That is, the fields of number density, velocity/momentum density and temperature are calculated over a rectangular grid and stored at a HDF5 file for further processing.Then, the obtained field profiles are used to compute the relevant macroscopic electronic quantities, such as drain-to-source current or dissipated power, and, particularly, the electric dipole moment of the charge distribution, which can subsequently be used to characterize the emitted far-field radiation. This is mostly done by numerically integrating the data in space, with the necessary time derivatives obtained by Gaussian convolution.