Step 1: determining \Delta \boldsymbol{V}^*

In this step the first correction to the momentum is determined. To this end the Right Hand Side (RHS) of the equation system is calculated. This part is parallelized (multithreading) since it involves a loop over all elements which can be nicely cut into pieces. The equations are solved in routine solveeq in an iterative way. This is necessary, since the LHS has been approximated by lumping. The number of iterations is set in solveeq.f and is called maxit. Right now, it has the value 1, which means that no iterations take place. If the user wishes to change this, the source code has to be recompiled. After solving the equations the variables are moved from a degree of freedom representation to a nodal based storage in resultsv1 (and stored in field v(1,*),v(2,*) and v(3,*)). Notice that \Delta \boldsymbol{V}^* is not added to \boldsymbol{V} at this point. Next, \Delta \boldsymbol{V}^* is changed such that the velocity boundary conditions are matched. These conditions are applied in the form \rho_i \boldsymbol{V}_{i+1}, where \rho_i is the density at the start of the increment and \boldsymbol{V}_{i+1} is the velocity boundary condition corresponding to the time at the end of the increment ( \rho_{i+1} is not known at this point).