We solve the wind flow across entire offshore wind farms. These huge simulations require highly parallel-efficient tools that can give valuable insights into physical phenomena. The filtered incompressible Navier-Stokes equations are numerically solved using SOD2D (Spectral high-Order coDe 2 solve partial Differential equations), a low- dissipation spectral element method (SEM) open-source code. This code has been developed in-house to efficiently leverage the computational power deployed worldwide. It runs on GPU or CPU architectures. The code is written in Fortran and uses MPI and OpenACC directives to provide parallelism at both coarse and fine-grained levels. For the temporal discretisation, a BDF-EXT3 high-order operator splitting approach is used to solve the velocity-pressure coupling. A projection stabilisation precludes numerical oscillations due to dominant convection while introducing low numerical dissipation. The aliasing effects of the reduced order integration caused by employing SEM integration for convective terms are countered with a skew-symmetric splitting. The inflow conditions are generated using a concurrent-precursor method, which simulates the ABL using lateral periodic boundaries. A proportional integral derivative control algorithm is implemented, to obtain the desired average wind speed and direction at hub height. The goal is to achieve the target wind velocity as soon as possible; minimising inertial oscillations. The rotation of the velocity field in the precursor simulation is enforced by a rotation of the frame of reference, without affecting the wind profile. A novel initial condition is used to reduce the cpu time to statistically steady. A conventionally neutral boundary layer (CNBL) is modelled, with a strong capping inversion above the ABL. The capping inversion can be displaced upwards by the presence of the wind farm, generating gravity waves that affect the distribution of pressure across the wind farm. The subgrid turbulence models ILSA [3] and Vreman are used, and compared, to obtain the wind field. Gravity waves are reflected on the boundaries. These reflections are minimised adding damping layers close to the boundaries. We follow the work of Khan et al. with some modifications. The code and methodology are validated by the simulation of a wind farm immersed in a shallow CNBL of 300m, and by the simulation of an offshore wind farm with hundreds of wind turbines operated by OceanWinds.