/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | faSavageHutterFOAM
    \\  /    A nd           | Copyright (C) 2017 Matthias Rauter
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Application
    faSavageHutterFoam

Description
    A depth-integrated solver for shallow granular flows.
    The solver is based on the Finite Area Method.
	Experimetnal version of an explicit solver

Author
    Matthias Rauter matthias.rauter@uibk.ac.at

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "faCFD.H"
#include "frictionModel.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
#   include "setRootCase.H"
#   include "createTime.H"
#   include "createMesh.H"
#   include "createFaMesh.H"
#   include "readGravitationalAcceleration.H"
#   include "createFaFields.H"
#   include "createFvFields.H"
#   include "readTransportProperties.H"
#   include "createTimeControls.H"

    Info << "\nStarting time loop\n" << endl;


    Switch initDeltaT =
        runTime.controlDict().lookup("initDeltaT");



    if (initDeltaT)
    {
#       include "readTimeControls.H"
#       include "surfaceCourantNo.H"
        runTime.setDeltaT(
            min(maxCo/(CoNum + SMALL)*runTime.deltaT().value(), maxDeltaT)
        );
    }

#       include "readSolutionControls.H"
#       include "readTimeControls.H"
#       include "surfaceCourantNo.H"
#       include "setDeltaT.H"

    scalar oldTime = runTime.value();

    bool adams = (ddtScheme == "adams");
    bool rungekutta = (ddtScheme == "rungekutta");
    if (adams)
        Info << "Using explicit Adams-Bashorth method" << endl;
    else if (rungekutta)
        Info << "Using explicit Runge-Kutta method" << endl;
    else
        Info << "Using semi-implicit Runge-Kutta method" << endl;

    scalar totVolIni = sum(h*aMesh.S()).value();

    while (runTime.run())
    {

        runTime++;
        scalar totVol = sum(h*aMesh.S()).value();
        Info << "Time = " << runTime.timeName() << nl <<
                "Total volume = " << totVol << "(" << (totVol/totVolIni*100) << "%)" << nl << endl;


        if (adams)
        {
            // factors for non-uniform time-stepping with Adams-Method
            scalar c = 0.;
            if (h.timeIndex() != h.oldTime().timeIndex())
            {
                scalar dtOld = runTime.value()-oldTime;
                scalar dtNew = runTime.deltaT().value();

                c = dtNew/dtOld;
            }
            else    // c = 0 leads to explicit euler
            {
                c = 0;
            }


            oldTime = runTime.value();

            ddthUs.storeOldTime();
            ddth.storeOldTime();
            tau.storeOldTime();
#           include "calcBasalstress.H"
#           include "calcFluxes.h"

            h = h+runTime.deltaT()/2.*(ddth*(2.+c) - ddth.oldTime()*c);
            hUs = hUs+runTime.deltaT()/2.*((ddthUs+tau)*(2.+c) - (ddthUs.oldTime()+tau.oldTime())*c);
        }
        else if (rungekutta)
        {
#           include "calcBasalstress.H"
#           include "calcFluxes.h"

            h.storeOldTime();
            hUs.storeOldTime();

            h = h+runTime.deltaT()/2.*(ddth);
            hUs = hUs+runTime.deltaT()/2.*ddthUs
                     ;//+runTime.deltaT()/2.*tau;

            vectorField pred = hUs+runTime.deltaT()/2.*tau;
            forAll(hUs, i)
            {
                for(int j = 0; j<3; j++)
                {
                    if ((pred[i][j]*tau[i][j]) > 0)//(mag(hUs[i][j]) < mag(tau[i][j])*runTime.deltaT().value()/2.)
                        hUs[i][j] = 0;
                    else
                        hUs[i][j] = pred[i][j];
                }
            }

            hUs -= n*(n & hUs);

            forAll(hUs, i)
            {
                if(h[i] < hminflux.value())
                    hUs[i] = vector(0,0,0);
            }

            h = max(h, dimensionedScalar("0", h.dimensions(), 0));

            Us = hUs/max(h, hminflux);

#           include "calcBasalstress.H"
#           include "calcFluxes.h"

            h = h.oldTime()+runTime.deltaT()*(ddth);
            hUs = hUs.oldTime()+runTime.deltaT()*ddthUs
                               ;//+runTime.deltaT()*tau;


            pred = hUs+runTime.deltaT()*tau;
            forAll(hUs, i)
            {
                for(int j = 0; j<3; j++)
                {
                    if ((pred[i][j]*tau[i][j]) > 0)//(mag(hUs[i][j]) < mag(tau[i][j])*runTime.deltaT().value()/2.)
                        hUs[i][j] = 0;
                    else
                        hUs[i][j] = pred[i][j];
                }
            }

            hUs -= n*(n & hUs);

            forAll(hUs, i)
            {
                if(h[i] < hminflux.value())
                    hUs[i] = vector(0,0,0);
            }

        }

        h = max(h, dimensionedScalar("0", h.dimensions(), 0));

        Us = hUs/max(h, hminflux);

        hUs.correctBoundaryConditions();
        h.correctBoundaryConditions();


        if (runTime.outputTime())
        {
            vsm.mapToVolume(hUs, HU.boundaryField());
            vsm.mapToVolume(Us, U.boundaryField());
            vsm.mapToVolume(h, H.boundaryField());
            vsm.mapToVolume(pb, Pb.boundaryField());

            runTime.write();
        }


        Info << "ExecutionTime = "
            << scalar(runTime.elapsedCpuTime())
            << " s\n" << endl << endl;


#       include "readSolutionControls.H"
#       include "readTimeControls.H"
#       include "surfaceCourantNo.H"
#       include "setDeltaT.H"

    }

    return(0);
}

// ************************************************************************* //
