 
 __device__ double RK_ErrorNorm(double* scale, double* DY) {
 	double sum = 0;
-	#pragma unroll 1
+	#pragma unroll 8
 	for (int i = 0; i < NSP; ++i){
 		sum += (scale[i] * scale[i] * DY[i] * DY[i]);
 	}

            }
            NewtonIncrementOld = NewtonIncrement;
            //! Update solution: Z4 <-- Z4 + DZ4
            #pragma unroll 1
            for (int i = 0; i < NSP; i++) {
            	Z4[i] += DZ4[i];
            }

#ifdef SDIRK_ERROR
		//!~~~>   Prepare the loop-independent part of the right-hand side
		//!       G = H*rkBgam(0)*F0 + rkTheta(1)*Z1 + rkTheta(2)*Z2 + rkTheta(3)*Z3
		#pragma unroll 8
		for (int i = 0; i < NSP; i++) {
			Z4[i] = Z3[i];
			G[i] = rkBgam[0]*F0[i]*H + rkTheta[0] * Z1[i] + rkTheta[1] * Z2[i] + rkTheta[2] * Z3[i];
