/**
 * AugRie_CUDA.h
 *
 ****
 **** Approximate Augmented Riemann Solver for the Shallow Water Equations
 **** vectorized using SIMD-Extensions
 ****
 *
 *  Created on: May 28, 2013
 *  Last Update: Jul 14, 2013
 *
 ****
 *
 *  Author: Alexander Breuer
 *    Homepage: http://www5.in.tum.de/wiki/index.php/Dipl.-Math._Alexander_Breuer
 *    E-Mail: breuera AT in.tum.de
 *
 *  Some optimizations: Martin Schreiber
 *    Homepage: http://www5.in.tum.de/wiki/index.php/Martin_Schreiber
 *    E-Mail: schreibm AT in.tum.de
 * 
 *  CUDA: Wolfgang Hölzl
 *    E-Mail: hoelzlw AT in.tum.de
 *
 ****
 *
 * (Main) Literature:
 *
 *   @phdthesis{george2006finite,
 *              Author = {George, D.L.},
 *              Title = {Finite volume methods and adaptive refinement for tsunami propagation and inundation},
 *              Year = {2006}}
 *
 *   @article{george2008augmented,
 *            Author = {George, D.L.},
 *            Journal = {Journal of Computational Physics},
 *            Number = {6},
 *            Pages = {3089--3113},
 *            Publisher = {Elsevier},
 *            Title = {Augmented Riemann solvers for the shallow water equations over variable topography with steady states and inundation},
 *            Volume = {227},
 *            Year = {2008}}
 *
 *   @book{leveque2002finite,
 *         Author = {LeVeque, R. J.},
 *         Date-Added = {2011-09-13 14:09:31 +0000},
 *         Date-Modified = {2011-10-31 09:46:40 +0000},
 *         Publisher = {Cambridge University Press},
 *         Title = {Finite Volume Methods for Hyperbolic Problems},
 *         Volume = {31},
 *         Year = {2002}}
 *
 *   @webpage{levequeclawpack,
 *            Author = {LeVeque, R. J.},
 *            Lastchecked = {January, 05, 2011},
 *            Title = {Clawpack Sofware},
 *            Url = {https://github.com/clawpack/clawpack-4.x/blob/master/geoclaw/2d/lib}}
 *
 ****
 *
 * Acknowledgments:
 *   Special thanks go to R.J. LeVeque and D.L. George for publishing their code
 *   and the corresponding documentation (-> Literature).
 */


/*
 * TODO: store nLow/nHigh variables in array[2]
 */


#ifndef AUGRIE_CUDA_HPP_
#define AUGRIE_CUDA_HPP_

#define FLOAT32

#include "SIMD_TYPES.hpp"

#include <cassert>
#include <cmath>
#include <iostream>
#include <string>
#include <vector>


const integer DryDry = 0;
const integer WetWet = 1;
const integer WetDryInundation = 2;
const integer WetDryWall = 3;
const integer WetDryWallInundation = 4;
const integer DryWetInundation = 5;
const integer DryWetWall = 6;
const integer DryWetWallInundation = 7;

const integer DrySingleRarefaction = 0;
const integer SingleRarefactionDry = 1;
const integer ShockShock = 2;
const integer ShockRarefaction = 3;
const integer RarefactionShock = 4;
const integer RarefactionRarefaction = 5;

__device__
inline void
computeMiddleState (
	const real &i_hLeft, const real &i_hRight,
	const real &i_uLeft, const real &i_uRight,
	const real dryTol, const real newtonTol, const real g, const real sqrt_g,
	const unsigned int &i_maxNumberOfNewtonIterations,
	real& o_hMiddle,
	real o_middleStateSpeeds[2]
);

/**
 * Compute net updates for the left/right cell of the edge.
 *
 * maxWaveSpeed will be set to the maximum (linearized) wave speed => CFL
 */
__device__
void
augRieComputeNetUpdates (
	const real i_hLeft, const real i_hRight,
	const real i_huLeft, const real i_huRight,
	const real i_bLeft, const real i_bRight,

	const real g, const real dryTol, const real newtonTol, const real zeroTol, const unsigned int maxNumberOfNewtonIterations,

	real o_netUpdates[5]
)
{
	real hLeft = i_hLeft;			real hRight = i_hRight;
	real uLeft = static_cast<real>(0);	real uRight = static_cast<real>(0);
	real huLeft = i_huLeft;			real huRight = i_huRight;
	real bLeft = i_bLeft;			real bRight = i_bRight;

	//declare variables which are used over and over again
	const real sqrt_g = sqrtf(g);
	real sqrt_g_hLeft;
	real sqrt_g_hRight;

	real sqrt_hLeft;
	real sqrt_hRight;

	//set speeds to zero (will be determined later)
	uLeft = uRight = 0.;

	//reset net updates and the maximum wave speed
	o_netUpdates[0] = o_netUpdates[1] = o_netUpdates[2] = o_netUpdates[3] = static_cast<real>(0);
	o_netUpdates[4] = static_cast<real>(0);

	real hMiddle = static_cast<real>(0);
	real middleStateSpeeds[2] = { static_cast<real>(0) };

	//determine the wet/dry state and compute local variables correspondingly

/***************************************************************************************
 * Determine Wet Dry State Begin
 **************************************************************************************/
	integer wetDryState;
	//compute speeds or set them to zero (dry cells)
	if (hLeft > dryTol) {
		uLeft = huLeft / hLeft;
	} else {
		bLeft += hLeft;
		hLeft = huLeft = uLeft = 0;
	}

	if (hRight > dryTol) {
		uRight = huRight / hRight;
	} else {
		bRight += hRight;
		hRight = huRight = uRight = 0;
	}

	if (hLeft >= dryTol and hRight >= dryTol) {
		//test for simple wet/wet case since this is most probably the
		//most frequently executed branch
		wetDryState = WetWet;
	} else if (hLeft < dryTol and hRight < dryTol) {
		//check for the dry/dry-case
		wetDryState = DryDry;
	} else if (hLeft < dryTol and hRight + bRight > bLeft) {
		//we have a shoreline: one cell dry, one cell wet

		//check for simple inundation problems
		// (=> dry cell lies lower than the wet cell)
		wetDryState = DryWetInundation;
	} else if (hRight < dryTol and hLeft + bLeft > bRight) {
		wetDryState = WetDryInundation;
	} else if (hLeft < dryTol) {
		//dry cell lies higher than the wet cell
		//lets check if the momentum is able to overcome the difference in height
		//  => solve homogeneous Riemann-problem to determine the middle state height
		//     which would arise if there is a wall (wall-boundary-condition)
		//       \cite[ch. 6.8.2]{george2006finite})
		//       \cite[ch. 5.2]{george2008augmented}
		computeMiddleState(
			hRight, hRight,
			-uRight, uRight,
			dryTol, newtonTol, g, sqrt_g,
			maxNumberOfNewtonIterations,
			hMiddle, middleStateSpeeds
		);

		if (hMiddle + bRight > bLeft) {
			//momentum is large enough, continue with the original values
			//          bLeft = o_hMiddle + bRight;
			wetDryState = DryWetWallInundation;
		} else {
			//momentum is not large enough, use wall-boundary-values
			hLeft = hRight;
			uLeft = -uRight;
			huLeft = -huRight;
			bLeft = bRight = static_cast<real>(0);
			wetDryState = DryWetWall;
		}
	} else if (hRight < dryTol) {
		//lets check if the momentum is able to overcome the difference in height
		//  => solve homogeneous Riemann-problem to determine the middle state height
		//     which would arise if there is a wall (wall-boundary-condition)
		//       \cite[ch. 6.8.2]{george2006finite})
		//       \cite[ch. 5.2]{george2008augmented}
		computeMiddleState(
			hLeft, hLeft,
			uLeft, -uLeft,
			dryTol, newtonTol, g, sqrt_g,
			maxNumberOfNewtonIterations,
			hMiddle, middleStateSpeeds
		);

		if (hMiddle + bLeft > bRight) {
			//momentum is large enough, continue with the original values
			//          bRight = o_hMiddle + bLeft;
			wetDryState = WetDryWallInundation;
		} else {
			hRight = hLeft;
			uRight = -uLeft;
			huRight = -huLeft;
			bRight = bLeft = static_cast<real>(0);
			wetDryState = WetDryWall;
		}
	} else {
		//done with all cases
		assert(false);
	}

	//limit the effect of the source term if there is a "wall"
	//\cite[end of ch. 5.2?]{george2008augmented}
	//\cite[rpn2ez_fast_geo.f][levequeclawpack]
	if (wetDryState == DryWetWallInundation) {
		bLeft = hRight + bRight;
	} else if (wetDryState == WetDryWallInundation) {
		bRight = hLeft + bLeft;
	}
/***************************************************************************************
 * Determine Wet Dry State End
 **************************************************************************************/

	if (wetDryState != DryDry) {
		//precompute some terms which are fixed during
		//the computation after some specific point
		sqrt_hLeft = std::sqrt(hLeft);
		sqrt_hRight = std::sqrt(hRight);

		sqrt_g_hLeft = sqrt_g*sqrt_hLeft;
		sqrt_g_hRight = sqrt_g*sqrt_hRight;


		//where to store the three waves
		real fWaves[3][2];
		//and their speeds
		real waveSpeeds[3];

		//compute the augmented decomposition
		//  (thats the place where the computational work is done..)
/***************************************************************************************
 * Compute Wave Decomposition Begin
 **************************************************************************************/
		//compute eigenvalues of the jacobian matrices in states Q_{i-1} and Q_{i} (char. speeds)
		const real characteristicSpeeds[2] = {
			uLeft - sqrt_g_hLeft,
			uRight + sqrt_g_hRight
		};

		//compute "Roe speeds"
		const real hRoe = static_cast<real>(0.5) * (hRight + hLeft);
		const real uRoe = (uLeft * sqrt_hLeft + uRight * sqrt_hRight) / (sqrt_hLeft + sqrt_hRight);

		//optimization for dumb compilers
		const real sqrt_g_hRoe = std::sqrt(g * hRoe);
		const real roeSpeeds[2] = {
			uRoe - sqrt_g_hRoe,
			uRoe + sqrt_g_hRoe
		};

		//compute the middle state of the homogeneous Riemann-Problem
		if (wetDryState != WetDryWall and wetDryState != DryWetWall) {
			//case WDW and DWW was computed in determineWetDryState already
			computeMiddleState(
				hLeft, hRight,
				uLeft, uRight,
				dryTol, newtonTol, g, sqrt_g,
				1,
				hMiddle, middleStateSpeeds
			);
		}

		//compute extended eindfeldt speeds (einfeldt speeds + middle state speeds)
		//  \cite[ch. 5.2]{george2008augmented}, \cite[ch. 6.8]{george2006finite}
		real extEinfeldtSpeeds[2] = {static_cast<real>(0), static_cast<real>(0)};
		if (wetDryState == WetWet or wetDryState == WetDryWall or wetDryState == DryWetWall) {
			extEinfeldtSpeeds[0] = fmin(characteristicSpeeds[0], roeSpeeds[0]);
			extEinfeldtSpeeds[0] = fmin(extEinfeldtSpeeds[0], middleStateSpeeds[1]);

			extEinfeldtSpeeds[1] = fmax(characteristicSpeeds[1], roeSpeeds[1]);
			extEinfeldtSpeeds[1] = fmax(extEinfeldtSpeeds[1], middleStateSpeeds[0]);
		} else if (hLeft < dryTol) {
			//ignore undefined speeds
			extEinfeldtSpeeds[0] = fmin(roeSpeeds[0], middleStateSpeeds[1]);
			extEinfeldtSpeeds[1] = fmax(characteristicSpeeds[1], roeSpeeds[1]);

			assert(middleStateSpeeds[0] < extEinfeldtSpeeds[1]);
		} else if (hRight < dryTol) {
			//ignore undefined speeds
			extEinfeldtSpeeds[0] = fmin(characteristicSpeeds[0], roeSpeeds[0]);
			extEinfeldtSpeeds[1] = fmax(roeSpeeds[1], middleStateSpeeds[0]);

			assert(middleStateSpeeds[1] > extEinfeldtSpeeds[0]);
		} else {
			assert(false);
		}

		//HLL middle state
		//  \cite[theorem 3.1]{george2006finite}, \cite[ch. 4.1]{george2008augmented}
		const real hLLMiddleHeight = fmax((huLeft - huRight + extEinfeldtSpeeds[1] * hRight - extEinfeldtSpeeds[0] * hLeft) / (extEinfeldtSpeeds[1] - extEinfeldtSpeeds[0]), static_cast<real>(0));

		//define eigenvalues
		const real eigenValues[3] = {
			extEinfeldtSpeeds[0],
			static_cast<real>(0.5) * (extEinfeldtSpeeds[0] + extEinfeldtSpeeds[1]),
			extEinfeldtSpeeds[1]
		};
		//define eigenvectors
		const real eigenVectors[3][3] = {
			{
				static_cast<real>(1),
				static_cast<real>(0),
				static_cast<real>(1)
			},
			{
				eigenValues[0],
				static_cast<real>(0),
				eigenValues[2]
			},
			{
				eigenValues[0] * eigenValues[0],
				static_cast<real>(1),
				eigenValues[2] * eigenValues[2]
			}
		};


		//compute rarefaction corrector wave
		//  \cite[ch. 6.7.2]{george2006finite}, \cite[ch. 5.1]{george2008augmented}
#pragma message "strongRarefaction set to false. Further investigations needed about senseless initialization of eigenValues[1]"

		//compute the jump in state
		real rightHandSide[3] = {
			hRight - hLeft,
			huRight - huLeft,
			huRight * uRight + static_cast<real>(0.5) * g * hRight * hRight - (huLeft * uLeft + static_cast<real>(0.5) * g * hLeft * hLeft)
		};

		//compute steady state wave
		//  \cite[ch. 4.2.1 \& app. A]{george2008augmented}, \cite[ch. 6.2 \& ch. 4.4]{george2006finite}
		const real hBar = (hLeft + hRight) * static_cast<real>(0.5);

		real steadyStateWave[2] = {
			-(bRight - bLeft),
			-g * hBar * (bRight - bLeft)
		};

		//preserve depth-positivity
		//  \cite[ch. 6.5.2]{george2006finite}, \cite[ch. 4.2.3]{george2008augmented}
		if (eigenValues[0] < -zeroTol and eigenValues[2] > zeroTol) {
			//subsonic
			steadyStateWave[0] = fmax(steadyStateWave[0], hLLMiddleHeight * (eigenValues[2] - eigenValues[0]) / eigenValues[0]);
			steadyStateWave[0] = fmin(steadyStateWave[0], hLLMiddleHeight * (eigenValues[2] - eigenValues[0]) / eigenValues[2]);
		} else if (eigenValues[0] > zeroTol) {
			//supersonic right TODO: motivation?
			steadyStateWave[0] = fmax(steadyStateWave[0], -hLeft);
			steadyStateWave[0] = fmin(steadyStateWave[0], hLLMiddleHeight * (eigenValues[2] - eigenValues[0]) / eigenValues[0]);
		} else if (eigenValues[2] < -zeroTol) {
			//supersonic left TODO: motivation?
			steadyStateWave[0] = fmax(steadyStateWave[0], hLLMiddleHeight * (eigenValues[2] - eigenValues[0]) / eigenValues[2]);
			steadyStateWave[0] = fmin(steadyStateWave[0], hRight);
		}

		//Limit the effect of the source term
		//  \cite[ch. 6.4.2]{george2006finite}
		steadyStateWave[1] = fmin(steadyStateWave[1], g * fmax(-hLeft * (bRight - bLeft), -hRight * (bRight - bLeft)));
		steadyStateWave[1] = fmax(steadyStateWave[1], g * fmin(-hLeft * (bRight - bLeft), -hRight * (bRight - bLeft)));

		rightHandSide[0] -= steadyStateWave[0];
		//rightHandSide[1]: no source term
		rightHandSide[2] -= steadyStateWave[1];

		//everything is ready, solve the equations!
/***************************************************************************************
 * Solve linear equation begin
 **************************************************************************************/
		// compute inverse of 3x3 matrix
		const real m[3][3] = {
			{
				(eigenVectors[1][1] * eigenVectors[2][2] - eigenVectors[1][2] * eigenVectors[2][1]),
				-(eigenVectors[0][1] * eigenVectors[2][2] - eigenVectors[0][2] * eigenVectors[2][1]),
				(eigenVectors[0][1] * eigenVectors[1][2] - eigenVectors[0][2] * eigenVectors[1][1])
			},
			{
				-(eigenVectors[1][0] * eigenVectors[2][2] - eigenVectors[1][2] * eigenVectors[2][0]),
				(eigenVectors[0][0] * eigenVectors[2][2] - eigenVectors[0][2] * eigenVectors[2][0]),
				-(eigenVectors[0][0] * eigenVectors[1][2] - eigenVectors[0][2] * eigenVectors[1][0])
			},
			{
				(eigenVectors[1][0] * eigenVectors[2][1] - eigenVectors[1][1] * eigenVectors[2][0]),
				-(eigenVectors[0][0] * eigenVectors[2][1] - eigenVectors[0][1] * eigenVectors[2][0]),
				(eigenVectors[0][0] * eigenVectors[1][1] - eigenVectors[0][1] * eigenVectors[1][0])
			}
		};
		const real d = (eigenVectors[0][0] * m[0][0] + eigenVectors[0][1] * m[1][0] + eigenVectors[0][2] * m[2][0]);

		// m stores not really the inverse matrix, but the inverse multiplied by d
		const real s = 1 / d;

		// compute m*rightHandSide
		const real beta[3] = {
			(m[0][0] * rightHandSide[0] + m[0][1] * rightHandSide[1] + m[0][2] * rightHandSide[2]) * s,
			(m[1][0] * rightHandSide[0] + m[1][1] * rightHandSide[1] + m[1][2] * rightHandSide[2]) * s,
			(m[2][0] * rightHandSide[0] + m[2][1] * rightHandSide[1] + m[2][2] * rightHandSide[2]) * s
		};
/***************************************************************************************
 * Solve linear equation end
 **************************************************************************************/

		//compute f-waves and wave-speeds
		if (wetDryState == WetDryWall) {
			//zero ghost updates (wall boundary)
			//care about the left going wave (0) only
			fWaves[0][0] = beta[0] * eigenVectors[1][0];
			fWaves[0][1] = beta[0] * eigenVectors[2][0];

			//set the rest to zero
			fWaves[1][0] = fWaves[1][1] = static_cast<real>(0);
			fWaves[2][0] = fWaves[2][1] = static_cast<real>(0);

			waveSpeeds[0] = eigenValues[0];
			waveSpeeds[1] = waveSpeeds[2] = static_cast<real>(0);

			assert(eigenValues[0] < zeroTol);
		} else if (wetDryState == DryWetWall) {
			//zero ghost updates (wall boundary)
			//care about the right going wave (2) only
			fWaves[2][0] = beta[2] * eigenVectors[1][2];
			fWaves[2][1] = beta[2] * eigenVectors[2][2];

			//set the rest to zero
			fWaves[0][0] = fWaves[0][1] = static_cast<real>(0);
			fWaves[1][0] = fWaves[1][1] = static_cast<real>(0);

			waveSpeeds[2] = eigenValues[2];
			waveSpeeds[0] = waveSpeeds[1] = static_cast<real>(0);

			assert(eigenValues[2] > -zeroTol);
		} else {
			//compute f-waves (default)
			for (int waveNumber = 0; waveNumber < 3; waveNumber++) {
				fWaves[waveNumber][0] = beta[waveNumber] * eigenVectors[1][waveNumber]; //select 2nd and
				fWaves[waveNumber][1] = beta[waveNumber] * eigenVectors[2][waveNumber]; //3rd component of the augmented decomposition
			}

			waveSpeeds[0] = eigenValues[0];
			waveSpeeds[1] = eigenValues[1];
			waveSpeeds[2] = eigenValues[2];
		}
/***************************************************************************************
 * Compute Wave Decomposition End
 **************************************************************************************/


		//compute the updates from the three propagating waves
		//A^-\delta Q = \sum{s[i]<0} \beta[i] * r[i] = A^-\delta Q = \sum{s[i]<0} Z^i
		//A^+\delta Q = \sum{s[i]>0} \beta[i] * r[i] = A^-\delta Q = \sum{s[i]<0} Z^i
		for (int waveNumber = 0; waveNumber < 3; waveNumber++) {
			if (waveSpeeds[waveNumber] < -zeroTol) {
				//left going
				o_netUpdates[0] += fWaves[waveNumber][0];
				o_netUpdates[2] += fWaves[waveNumber][1];
			} else if (waveSpeeds[waveNumber] > zeroTol) {
				//right going
				o_netUpdates[1] += fWaves[waveNumber][0];
				o_netUpdates[3] += fWaves[waveNumber][1];
			} else {
				//TODO: this case should not happen mathematically, but it does. Where is the bug? Machine accuracy only?
				o_netUpdates[0] += static_cast<real>(0.5) * fWaves[waveNumber][0];
				o_netUpdates[2] += static_cast<real>(0.5) * fWaves[waveNumber][1];

				o_netUpdates[1] += static_cast<real>(0.5) * fWaves[waveNumber][0];
				o_netUpdates[3] += static_cast<real>(0.5) * fWaves[waveNumber][1];
			}
		}

		//compute maximum wave speed (-> CFL-condition)
		waveSpeeds[0] = fabs(waveSpeeds[0]);
		waveSpeeds[1] = fabs(waveSpeeds[1]);
		waveSpeeds[2] = fabs(waveSpeeds[2]);

		o_netUpdates[4] = fmax(waveSpeeds[0], waveSpeeds[1]);
		o_netUpdates[4] = fmax(o_netUpdates[4], waveSpeeds[2]);
	}
}

/**
 * Computes the middle state of the homogeneous Riemann-problem.
 *   -> (\cite[ch. 13]{leveque2002finite})
 *
 * @param i_hLeft height on the left side of the edge.
 * @param i_hRight height on the right side of the edge.
 * @param i_uLeft velocity on the left side of the edge.
 * @param i_uRight velocity on the right side of the edge.
 * @param i_huLeft momentum on the left side of the edge.
 * @param i_huRight momentum on the right side of the edge.
 * @param i_maxNumberOfNewtonIterations maximum number of Newton iterations.
 */
__device__
inline void
computeMiddleState (
	const real &i_hLeft, const real &i_hRight,
	const real &i_uLeft, const real &i_uRight,
	const real dryTol, const real newtonTol, const real g, const real sqrt_g,
	const unsigned int &i_maxNumberOfNewtonIterations,
	real& o_hMiddle,
	real o_middleStateSpeeds[2]
)
{
	//set everything to zero
	o_hMiddle = static_cast<real>(0);
	o_middleStateSpeeds[0] = static_cast<real>(0);
	o_middleStateSpeeds[1] = static_cast<real>(0);

	//compute local square roots
	//(not necessarily the same ones as in computeNetUpdates!)
	real l_sqrt_g_hRight = std::sqrt(g * i_hRight);
	real l_sqrt_g_hLeft = std::sqrt(g * i_hLeft);

	//single rarefaction in the case of a wet/dry interface
	integer riemannStructure;
	if (i_hLeft < dryTol) {
		o_middleStateSpeeds[1] = o_middleStateSpeeds[0] = i_uRight - static_cast<real>(2) * l_sqrt_g_hRight;
		riemannStructure = DrySingleRarefaction;
		return;
	} else if (i_hRight < dryTol) {
		o_middleStateSpeeds[0] = o_middleStateSpeeds[1] = i_uLeft + static_cast<real>(2) * l_sqrt_g_hLeft;
		riemannStructure = SingleRarefactionDry;
		return;
	}

	//determine the wave structure of the Riemann-problem
/***************************************************************************************
 * Determine riemann structure begin
 **************************************************************************************/
	const real hMin = fmin(i_hLeft, i_hRight);
	const real hMax = fmax(i_hLeft, i_hRight);

	const real uDif = i_uRight - i_uLeft;

	if (0 <= static_cast<real>(2) * (std::sqrt(g * hMin) - std::sqrt(g * hMax)) + uDif) {
		riemannStructure = RarefactionRarefaction;
	} else if ((hMax - hMin) * std::sqrt(g * static_cast<real>(0.5) * (1 / hMax + 1 / hMin)) + uDif <= 0) {
		riemannStructure = ShockShock;
	} else if (i_hLeft < i_hRight) {
		riemannStructure = ShockRarefaction;
	} else {
		riemannStructure = RarefactionShock;
	}
/***************************************************************************************
 * Determine riemann structure end
 **************************************************************************************/

	//will be computed later
	real sqrt_g_hMiddle = static_cast<real>(0);

	if (riemannStructure == ShockShock) {
		o_hMiddle = fmin(i_hLeft, i_hRight);

		real l_sqrtTermH[2] = {0, 0};

		for (unsigned int i = 0; i < i_maxNumberOfNewtonIterations; i++) {
			l_sqrtTermH[0] = std::sqrt(static_cast<real>(0.5) * g * ((o_hMiddle + i_hLeft) / (o_hMiddle * i_hLeft)));
			l_sqrtTermH[1] = std::sqrt(static_cast<real>(0.5) * g * ((o_hMiddle + i_hRight) / (o_hMiddle * i_hRight)));

			real phi = i_uRight - i_uLeft + (o_hMiddle - i_hLeft) * l_sqrtTermH[0] + (o_hMiddle - i_hRight) * l_sqrtTermH[1];

			if (std::fabs(phi) < newtonTol) {
				break;
			}

			real derivativePhi = l_sqrtTermH[0] + l_sqrtTermH[1]
					- static_cast<real>(0.25) * g * (o_hMiddle - i_hLeft) / (l_sqrtTermH[0] * o_hMiddle * o_hMiddle)
					- static_cast<real>(0.25) * g * (o_hMiddle - i_hRight) / (l_sqrtTermH[1] * o_hMiddle * o_hMiddle);

			o_hMiddle = o_hMiddle - phi / derivativePhi; //Newton step
			assert(o_hMiddle >= dryTol);
		}

		sqrt_g_hMiddle = std::sqrt(g * o_hMiddle);
	}

	if (riemannStructure == RarefactionRarefaction) {
		o_hMiddle = fmax(static_cast<real>(0), i_uLeft - i_uRight + static_cast<real>(2) * (l_sqrt_g_hLeft + l_sqrt_g_hRight));
		o_hMiddle = static_cast<real>(1) / (static_cast<real>(16) * g) * o_hMiddle * o_hMiddle;

		sqrt_g_hMiddle = std::sqrt(g * o_hMiddle);
	}

	if (riemannStructure == ShockRarefaction or riemannStructure == RarefactionShock) {
		real hMin, hMax;
		if (riemannStructure == ShockRarefaction) {
			hMin = i_hLeft;
			hMax = i_hRight;
		} else {
			hMin = i_hRight;
			hMax = i_hLeft;
		}

		o_hMiddle = hMin;

		sqrt_g_hMiddle = std::sqrt(g * o_hMiddle);
		real sqrt_g_hMax = std::sqrt(g * hMax);
		for (unsigned int i = 0; i < i_maxNumberOfNewtonIterations; i++) {
			real sqrtTermHMin = std::sqrt(static_cast<real>(0.5) * g * ((o_hMiddle + hMin) / (o_hMiddle * hMin)));

			real phi = i_uRight - i_uLeft + (o_hMiddle - hMin) * sqrtTermHMin + static_cast<real>(2) * (sqrt_g_hMiddle - sqrt_g_hMax);

			if (std::fabs(phi) < newtonTol) {
				break;
			}

			real derivativePhi = sqrtTermHMin - static_cast<real>(0.25) * g * (o_hMiddle - hMin) / (o_hMiddle * o_hMiddle * sqrtTermHMin) + sqrt_g / sqrt_g_hMiddle;

			o_hMiddle = o_hMiddle - phi / derivativePhi; //Newton step

			sqrt_g_hMiddle = std::sqrt(g * o_hMiddle);
		}
	}

	o_middleStateSpeeds[0] = i_uLeft + static_cast<real>(2) * l_sqrt_g_hLeft - static_cast<real>(3) * sqrt_g_hMiddle;
	o_middleStateSpeeds[1] = i_uRight - static_cast<real>(2) * l_sqrt_g_hRight + static_cast<real>(3) * sqrt_g_hMiddle;

	assert(o_hMiddle >= 0);
}

#endif /* AUGRIE_CUDA_HPP_ */
