package integrators;

import main.*;

/** BackwardEulerIntegrator is a custom-modified backward Euler method
with simple-minded adaptive stepsizing that exploits the same partially
explicit corrector formula as the SemiExplicitAPCIntegrator.  Basically
it makes an estimate of y(t+h) by first making such an estimate with
the forward Euler step y(t+h) = y(t) + hf(y(t)), then using that estimate
to calculate the backward Euler step y(t+h) = y(t) + hf(y(t+h)), then
iterating the backward step to convergence.  At high accuracy (<10E-4)
this is usually much slower than a higher-order method; however, since
it is "guaranteed" to be stable at low accuracy (10E-2 or even 10E-1) it
can often take much bigger steps than other methods, especially if the
parameters are such as to render the equations relatively stiff... of
course they are still low-accuracy steps. */
public class BackwardEulerIntegrator extends Integrator	{
	
	final String		type = "BackwardEuler";
	
	NodeValues		[]	yT;
	NodeValues		[]	fT;
	NodeValues		[]	yP;
	NodeValues		[]	fhat;
	NodeValues		[]	dfdyi;
	NodeValues		[]	yC;
	NodeValues		[]	yErr;

	final float			maxcut=10f, maxboost=4f;
	final float			D = 0.5f;
	final float			GAMMA = 1f/2f;
	final int			MAXIT=4;

	boolean				gotCurrentInfo=false, finalsSet=false;	// used in recursion
	boolean				warned = false;

	public BackwardEulerIntegrator() {}
	
	public void init(Model model)	{
		super.init(model);
		
		errorTolerance = 0.01f;
		
		yT = new NodeValues[num_cells * num_nodes];
		fT = new NodeValues[num_cells * num_nodes];
		yP = new NodeValues[num_cells * num_nodes];
		fhat = new NodeValues[num_cells * num_nodes];
		dfdyi = new NodeValues[num_cells * num_nodes];
		yC = new NodeValues[num_cells * num_nodes];
		yErr = new NodeValues[num_cells * num_nodes];
		
		for(int i = 0; i < num_cells; i++) {
			for(int j = 0; j < num_nodes; j++) {
				int k = (i * num_nodes) + j;
				yT[k] = new NodeValues();
				yT[k].numSides = model.cells[0].nodes[j].numSides;
				fT[k] = new NodeValues();
				fT[k].numSides = model.cells[0].nodes[j].numSides;
				yP[k] = new NodeValues();
				yP[k].numSides = model.cells[0].nodes[j].numSides;
				fhat[k] = new NodeValues();
				fhat[k].numSides = model.cells[0].nodes[j].numSides;
				dfdyi[k] = new NodeValues();
				dfdyi[k].numSides = model.cells[0].nodes[j].numSides;
				yC[k] = new NodeValues();
				yC[k].numSides = model.cells[0].nodes[j].numSides;
				yErr[k] = new NodeValues();
				yErr[k].numSides = model.cells[0].nodes[j].numSides;
			}
		}
	}
	
	public void reset()	{
		timestep = defaultStepsize;
		gotCurrentInfo = finalsSet = warned = false;
	}
	
	public float IntegrateOneStep()	{
		int i, j, k, l, iterCount=0;
		float yCold, delta = 0.0f, epsilon=0.0f, stepTaken = 0.0f, idealstep;
		finalsSet = false;
		boolean done = false;
		
		if(!gotCurrentInfo)	{
			getCurrentValues(yT);
			getDerivatives(fT);
			gotCurrentInfo = true;
		}
		for(i = 0; i < num_cells; i++) {
			for(j = 0; j < num_nodes; j++) {
				k = (i*num_nodes) + j;
				for(l = 0; l < yP[k].numSides; l++)	{
					yP[k].values[l] = yT[k].values[l] + (timestep * fT[k].values[l]);
				}
			}
		}
		setIntegrationValues(yP);
		copyNVarray(yP, yC);
		do	{
			getFhatAndDfdyi(fhat, dfdyi);
			epsilon = 0f;
			for(i = 0; i < num_cells; i++) {
				for(j = 0; j < num_nodes; j++) {
					k = (i*num_nodes) + j;
					for(l = 0; l < yP[k].numSides; l++)	{
						yCold = yC[k].values[l];
						yC[k].values[l] = (yT[k].values[l] + timestep * fhat[k].values[l]) /
														(1f - timestep * dfdyi[k].values[l]);
						yErr[k].values[l] = Math.abs(yC[k].values[l] - yCold);
						epsilon += yErr[k].values[l] * yErr[k].values[l];
					}
				}
			}
			setIntegrationValues(yC);
			iterCount++;
			epsilon = (float)Math.pow((double)epsilon,0.5);
			if(epsilon <= errorTolerance) done = true;
		}
		while(iterCount < MAXIT && !done);
		eEst = epsilon;	
		if(!done)	{			// failed to converge
			setIntegrationValues(yT);
			timestep /= 2.0f;
			stepTaken += IntegrateOneStep();
		}
		else	{
			delta = 0.0f;
			for(i = 0; i < num_cells; i++)	{
				for(j = 0; j < num_nodes; j++)	{
					k = (i*num_nodes)+j;
					for(l = 0; l < yP[k].numSides; l++)	{
						yErr[k].values[l] = Math.abs(yC[k].values[l] - yP[k].values[l]);
						delta += yErr[k].values[l] * yErr[k].values[l];
					}
				}
			}
			delta = D * (float)Math.pow((double)delta,0.5);
			idealstep = timestep * ((GAMMA * (errorTolerance))/delta);
			if(delta <= errorTolerance/4f && idealstep > timestep && iterCount <= MAXIT/2)	{	// accept step but try larger step next time 
				stepTaken += timestep;
				timestep = Math.max(timestep*maxboost,idealstep);
				setFinalValues(yC);
				setIntegrationValues(yC);
			}
			else if(delta > errorTolerance)	{	// reject this step and make correction
				setIntegrationValues(yT);
				timestep = Math.max(timestep/maxcut,idealstep);
				stepTaken += IntegrateOneStep();
			}
			else	{	// accept this step, leave timestep alone
				stepTaken += timestep;
				setFinalValues(yC);
				setIntegrationValues(yC);
			}
		}
		
		gotCurrentInfo = false;
		return stepTaken;

	}
	
	protected void getFhatAndDfdyi(NodeValues [] fhat_vals, NodeValues [] dfdyi_vals)	{
		for(int i = 0; i < num_cells; i++) {
			for(int j = 0; j < num_nodes; j++) {
				model.cells[i].nodes[j].getSeparatedDeriv(fhat_vals[(i*num_nodes)+j].values, dfdyi_vals[(i*num_nodes)+j].values);
			}
		}
	}
	
	
}