The implementation of three-dimensional magnetic fields, such as the ones of stellarators, in the GBS code (Ricci et al 2012 Plasma Phys. Control. Fusion54 124047; Giacomin et al 2022 J. Comput. Phys.464 111294) is presented, and simulation results are discussed. The geometrical operators appearing in the drif-reduced Braginskii equations evolved by GBS are expanded considering the typical parameter ordering of stellarator configurations. It turns out that most of the operators have a similar structure as the one implemented in the tokamak axisymmetric version of the code. In particular, the perpendicular laplacian only acts on the poloidal plane, which avoids the need of a three-dimensional solver for the electrostatic potential. The simulation of an island divertor stellarator is then presented, showing the derivation of the magnetic equilibrium in detail and extending the results in (Coelho et al 2022 Nucl. Fusion62 074004). Although the island magnetic field-lines divert the plasma towards the strike points of the walls, the islands do not seem to have a significant impact on the turbulence properties. The dominant mode, identified as interchange-driven, is field-aligned and breaks the stellarator toroidal symmetry. The radial and poloidal extensions of the mode are of the same order, in contrast to typical tokamak simulations. This has consequences on the poloidal dependence of turbulent transport.