//
// WaterDroplets.cs
//

using System;
using System.Collections.Generic;
using System.Text;

class WaterDroplets
{
    static public void setCmdLine(int argc, string[] argv, ref int customArgc, ref string[] customArgv)
    {
        int i=0;
        for (;i<argc;i++)
        {
	        int n = argv[i].Length;
	        customArgv[i] = argv[i];
	    }

        // Set -c ClientSDK, -? for help.
        customArgv[i] = "-c";
        customArgv[i+1] = "ClientSDK";
    }

    static void Main(string[] args)
    {
    	int    customArgc = args.Length+2;
        string[] customArgv = new string[customArgc];

        setCmdLine(args.Length, args, ref customArgc, ref customArgv);

        Console.Write("*************************\n\n");
        Console.Write("Water droplets simulation\n\n");
        Console.Write("*************************\n\n");

        if (libra.Init(customArgc, customArgv) != 0)
	        return;

        //libra.SetCurrentBackend(gBackend.CPU_BACKEND);

        run();
        
        libra.Shutdown();
    }

    static public void run()
    {
        // 
        libra.SetDefaultDataType(gDataType.GFLOAT32);

        // Setup data.
        const int gridSizeX = 512;
        const int gridSizeY = 512;
        float alpha = 1.0f;

        gVar A = null, x = null, b = null;
        gVar g1 = null, g2 = null, g3 = null, g4 = null;

        setupWaveEquationSystem(ref alpha, ref A, ref x, ref b,
                                ref g1, ref g2, ref g3, ref g4,
						        gridSizeX, gridSizeY);

        // Solve system
        solveRenderWaveEquation(ref alpha, ref A, ref x, ref b,
                                ref g1, ref g2, ref g3, ref g4,
						        gridSizeX, gridSizeY);
    }

    static public void solveRenderWaveEquation(ref float alpha, ref gVar A, ref gVar x, ref gVar b,
						        ref gVar g1, ref gVar g2, ref gVar g3, ref gVar g4,
						         int gridSizeX, int gridSizeY)
    {
        gVar M = libra.diag(A);

        gVar s = new gVar(gridSizeX, gridSizeY);
        gVar u = libra.zeros(x.rowCount(), x.columnCount());
        b = libra.zeros(x.rowCount(), x.columnCount());
        int rowCount = x.rowCount();
        double totalTime = 0.0, averageIterationTime;
        const int timeSteps = 5000;
        double startTime = libra.GetTime();

        // for each step in time...
        for (int i=0;i<timeSteps;++i)
        { 
	        if (i == 10)
		        startTime = libra.GetTime();

	        // Create new drops by scattering a random value
	        // into x at the random position
	        // Scatter heavy rain
	        libra.scatter(libra.rand(15, 1) * (float)(rowCount - 1.0f), libra.rand(15, 1)*0.2f-0.1f, x);
            
            	// Setup right hand side (b) of equation
	        b.assign (alpha * (libra.gather(g1, x) + 
				                libra.gather(g2, x) + 
				                 libra.gather(g3, x) + 
				                  libra.gather(g4, x)) + 
                                  (2.0f - 4.0f * alpha) * x - u);
            
            	// u <-- old solution
	        u.assign(x);
            
	        // x <-- solve for new solution
	        for (int j=0;j<3;++j)
		        x.assign(libra.jacobiIter(A, x, b));

	        // render solution x to screen
	        libra.render(libra.log10(u + 1.0f)*15.0f, (gRenderMode.GHEIGHTFIELD2D | gRenderMode.GSTATS), gridSizeX);

            GC.Collect(GC.MaxGeneration, GCCollectionMode.Forced);
            GC.WaitForFullGCComplete();
            GC.WaitForPendingFinalizers();
        }

        totalTime = libra.GetTime() - startTime;
        averageIterationTime = totalTime / (timeSteps-10);

        Console.Write("Total time (s) : " + Convert.ToString(String.Format("{0:0.00}", totalTime)) + "\n");
        Console.Write("Average iteration time (ms) : " + Convert.ToString(String.Format("{0:0.00}", averageIterationTime * 1000.0)) + "\n");
        Console.Write("\n");
    }

    unsafe static public void setupWaveEquationSystem(ref float alpha, ref gVar A, ref gVar x, ref gVar b,
							        ref gVar g1, ref gVar g2, ref gVar g3, ref gVar g4,
							        int gridSizeX, int gridSizeY)
    {
        //
        // 2D Wave equation
        // 

        const float C = 0.5F;
        const float h = 1.0F;
        const float dT = 1.0F;
    	
        int vectorSize = gridSizeX*gridSizeY;
        int matrixSizeX = vectorSize;
        int matrixSizeY = vectorSize;

        alpha = (dT*dT*C*C) / (2*h*h);
        float mainDiag = 4*alpha + 1;

        int i, offset = 0;
        const int bandCount = 5;
     
        float[] bandValues = new float[vectorSize*bandCount];
      
        // setup data main diagonal-Y
        for (i = 0; i < vectorSize; i++)
            bandValues[offset+i] = (i >= gridSizeX) ? -alpha : 0;

        offset += i;

        // setup data main diagonal-1
        for (i = 0; i < vectorSize; i++)
            bandValues[offset + i] = (i % gridSizeY > 0) ? -alpha : 0;

        offset += i;

        // setup data main diagonal
        for (i = 0; i < vectorSize; i++)
            bandValues[offset + i] = mainDiag;

        offset += i;

        for (i = 0; i < vectorSize; i++)
            bandValues[offset + i] = (i >= gridSizeX) ? -alpha : 0;

        offset += i;

        // setup data main diagonal-1
        for (i = 0; i < vectorSize; i++)
            bandValues[offset + i] = (i % gridSizeY > 0) ? -alpha : 0;

        offset += i;

        int[] bandOffsetDatas = new int[5]{-gridSizeX, -1, 0, 1, gridSizeX};
        
        // create the sparse matrix with input on band form
        A = libra.sparseBanded(matrixSizeX, matrixSizeY, bandValues, bandOffsetDatas, bandCount);

        gVar u = libra.floor(libra.step(0.0f, vectorSize - 1.0f, 1.0f) % (float)gridSizeX);
        gVar v = libra.floor(libra.div(libra.step(0.0f, vectorSize - 1.0f, 1.0f), (float)gridSizeX)); 
        g1 = libra.max(u - 1.0f, 0.0f) + v * (float)gridSizeX;
        g2 = libra.min(u + 1.0f, gridSizeX - 1.0f) + v * (float)gridSizeX;
        g3 = u + libra.max(v - 1.0f, 0.0f) * (float)gridSizeX;
        g4 = u + libra.min(v + 1.0f, gridSizeX - 1.0f) * (float)gridSizeX;

        // setup initial solution, x = 0
        x = libra.zeros(vectorSize, 1);

        // create drops by scattering a few random values into x
        int dropCount = (int)(gridSizeX / 20.0f);
        libra.scatter(libra.rand(dropCount, 1) * (float)(x.rowCount() - 1), libra.rand(dropCount, 1)*0.5f, x);
    }
}