In 1995, two groups <Kageyama et al., 1995; Glatzmaier and Roberts, 1995a, 1995b> reported results of numerical integrations of fully three-dimensional, fully nonlinear dynamos. Their papers were precursors of a stream of such models that have focused particularly on the geodynamo. They provide us, in unprecedented detail, with spectacular realizations of interesting geomagnetic field behaviors, such as secular variation and even polarity reversals. The proliferation of models has, however, created some confusion and apparently conflicting results. This can be partly attributed to the different ways in which different groups have modeled the core, normalized their equations, defined their dimensionless parameters, chosen their boundary conditions, and selected their energy sources. This has made it difficult to compare the results of different simulations directly. In this paper, we first try, as far as possible, to overcome this difficulty, so that all reported results can be compared on common ground. We then review the results, emphasizing three major topics: (1) onset and evolution of convection, (2) character of the magnetic field generated, and (3) comparison with the observed geomagnetic field. Although there are large differences in the way that the simulations are defined, the magnetic fields that they generate have some surprising similarities. The fields are dominated by the axial dipole. In some models they are most strongly generated in shear layers near the upper and lower boundaries and near the tangent cylinder, an imaginary surface touching the inner core on its equator. Convection rolls occur within which a type of the α effect distorts the toroidal field lines to create poloidal magnetic field. Some features of the models are found to strongly affect the fields that they produce. In particular, the boundary conditions defining the energy flow (e.g., an inhomogeneous heat flux or distribution of buoyancy sources) are very influential and have been extensively studied. They change the frequency and the mode of magnetic polarity reversals as well as the ratio in strengths of the dipole and nondipole moments. As the ultimate goal of geodynamo simulations is to explain the features of the real geomagnetic field, it is essential that proper comparisons be made between simulation results and observations. It is remarkable that polarity reversals reminiscent of the paleomagnetically observed field reversals have already been simulated by some of the models. Other features such as drift of the field, its secular variation, and statistical properties of Gauss coefficients are discussed in this paper and are compared with observations. These comparisons are ratherprimitive, not only because self-consistent dynamo models are still too new and too few but also because many of the observations (and especially the paleomagnetic data) are themselves not yet reliable or decisive enough. The aim of the third part of this paper is therefore more to demonstrate the potential use of simulations than to elucidate the nature of geomagnetic field generation.