""" -----------------------------------------------sanddune_erosion_operator

This script provides an Anuga operator that simulates the removal of sand
associated with over-topping of a sand dune by a tsunami. It is provided to
permit the exploration of a scenario where a dune line provides protection for
land(or water) behind the dune line that could be lost if part of the dune was
to be over topped by the first or successive waves of the tsunami.

The erosion mechanism assumes a clear water scour (no significant sediment
entrained in the approaching wave that could impact detachment rates and an
environment where the detached sediment is well within the transport capacity
of the water column in the eroding area(viz no settlement only erosion
within the target area).

The erosional relationships used are based on the work by Dr David Froelich.
Froelich (2002) IMPACT Project Field Tests 1 and 2 "Blind" Simulation by DaveF
Mo-i-Rana Norway, 2002

Erosional parameters of relevance extracted from the above are;

Bed shear  Tau_bed = Wd*G*(n**2)* (m**2) / (d**2.333) Pa (m is absolute momentum)
Where
Wd = water mass density,               (1000 kg/m3)
Sd = sediment mass density             ( 1800 kg/m3)
G = Accel due to gravity,              (9.8 m/sec/sec)
n = mannings n,                        (sand n = 0.025)
m = abs momentum                       ((mx**2+my**2)**0.5)
d = water depth                        (Anuga stage-elevation m)

Dune erosion will occur when bed shear stress > critical bed shear stress
from Froelich Table 1
Critical (detachment) bed shear stress for sand Tau_crit is 2.1 Pa
Detachment rate constant Kd is 0.0250 Kg/sec/m2/Pa and
detachment rate Sk = Kd*(tau_bed-Tau_Crit)/Sd m3/sec/m2 (Froelich Eq4)

Note That while Froelich provides erosional properties for other soils, this
operator is specific to the erosion of a sand dune.

The second process simulated by this operator is the collapse, fluidisation and
removal of sand from the dune system as a consequence of the above (vertically)
aligned erosion process creating face slopes that would be greater than the
angle of repose.  This is applied after the erosion computations for a particular
timestep have been completed and a new eroded surface computed. Each triangle within the
specified erosion area is then checked to see if the CG (elevation) of any neighbouring
triangles lie above or below the angle of repose line from the lowest neighbours
CG . If so then the neighbouring triangles CG (elevation) that is above the repose angle
is adjusted to lie at the angle of repose relative to the lowest neighbours CG.
This adjustment is also made to the current triangles CG if it also is above the repose angle
relative to the lowest neighbours CG. A check is performed to prevent any reduction in
elevation that would lower a triangles CG below the specified base level.
This is an approximate simplified lower bound approach to the full solution of surface
level adjustments based on true dip and strike.  While this algorithm allows
each triangles height  to be reviewed and adjusted up to three times in each timestep,
it does not guarantee  all triangles will be fully adjusted back to a stable face.
Given the multiple  iterations at each time step and small size of each time step
it is considered that it can however reasonably simulate the collapse of steep sand faces.

It is noted that a third removal process that involves the regressive fluidisation and
collapse of the downstream face from a rising phreatic surface within the dune from the
rising sea level. This has not been included in this model as it is considered a
significantly longer duration process than the other two processes.

NOTE - NOTE - NOTE - NOTE - NOTE - NOTE - NOTE - NOTE
Anuga is a 2D model and as such can not simulate vertical components of velocity where
they exist. Flow approaching and in particular down the downstream face of a dune will
have a significant vertical velocity component that will alter the rate of detachment
of sand. This code is therefore only an approximation of the realworld process and its
results should not be used where full simulation of the real world process is the intent.



Script Author E Rigby, ted.rigby@rienco.com.au
Version       1.00 October 2015
"""

from anuga.operators.base_operator import Operator
from anuga import Region
import pdb



class sanddune_erosion_operator(Operator, Region)  :

    def __init__(self,
                 domain,
                 threshold= 0.0,
                 base=0.0,
                 indices=None,
                 polygon=None,
                 center=None,
                 radius=None,
                 description = None,
                 label = None,
                 logging = False,
                 verbose = False):




        Operator.__init__(self, domain, description, label, logging, verbose)



        Region.__init__(self, domain,
                        indices=indices,
                        polygon=polygon,
                        center=center,
                        radius=radius,
                        verbose=verbose)

        #------------------------------------------
        # Local variables
        #------------------------------------------
        self.base = base


    def __call__(self):
        """
        Applies  removal of sand by erosion and  collapse due to slope instability to
        those triangles defined in indices by the specified polygon

        indices == [], then don't apply anywhere
        otherwise apply for the specific indices within the erosion area polygon
        """

        import math
        import numpy as num

        updated = True             # flag called OK

        if self.indices is []:     # empty list no polygon - return
            return (updated)

        #-------------------------------------------
        # set some sand and water properties
        #-------------------------------------------
        Wd       = 1000       # water mass density kg/m3
        Sd       = 1800       # sediment mass density kg/m3
        G        = 9.8        # acceleration due to gravity m/sec/sec
        n        = 0.030      # sand mannings n - mostly bare with undulations
        Tau_crit = 2.1        # critical (detachment) bed shear stress Pa
        Kd       = 0.025      # detachment factor Froelich Table 2 Kg/sec/m2/Pa
        Ra       = 34.0       # Repose angle in degrees for dry sand (note wet varies 30-45)
        Rs       = math.tan(Ra*math.pi/180)    # Repose slope

        #-------------------------------------------
        # get some useful model parameters
        #-------------------------------------------

        dt = self.get_timestep()

        #-----------------------------------------------------------------------------------------
        # Compute erosion depths during the timestep and update centroid elevations accordingly
        # Don't allow erosion below the specified base level
        #-----------------------------------------------------------------------------------------


        ind = self.indices                                                   # indices of triangles in polygon

        # store previous timestep centroid heights(depths) for later use
        h = self.stage_c[ind] - self.elev_c[ind]

        m = num.sqrt(self.xmom_c[ind]**2 + self.ymom_c[ind]**2)              # abs Momentum m2/sec vector
        d =(self.stage_c[ind]-self.elev_c[ind])                              # Depth meters vector

        # compute bed shear stress Pa
        Tau_bed  = Wd*G*(n**2)*(m**2)/((d**2.333)+ 0.000001)                 # refer Froelich 2002 vector

        # compute de (m) elevation change increment due to scour during timestep
        de = (Kd*(Tau_bed-Tau_crit)/Sd) *dt                                  # refer Froelich 2002 vector

        # no scour though if Tau_bed < Tau_crit ie if de is <= 0
        de = num.where(de > 0.0, de, 0.0)                                    # de=de whenever de>0 vector

        # Also ensure we don't erode below the base surface level
        self.elev_c[ind] = num.maximum(self.elev_c[ind] - de, num.minimum(self.elev_c[ind],self.base[ind]) )

        #-------------------------------------------------------------------------------------------
        # Reduce triangles elevations in any area where erosion has created unstable face slopes, so
        # that the final face slope lies at or about the angle of repose. This is programmed as a loop
        # over all triangles in the erosion poly - there may be speedier ways of doing this and this
        # only approximates the collapse as slopes not oriented to the dip angle.
        #--------------------------------------------------------------------------------------------
        neighbours = self.domain.surrogate_neighbours

        use_old_code = False
        if use_old_code:
            for i in range(len(ind))                                  :          # loop over all tris in poly

                neighbourindices = neighbours[ind[i]]                            # get the neighbour Indices for the current triangle
                if len(neighbourindices) == 3                         :          # only do this if all neighbours present

                    # locate lowest neighbour triangle
                    n0ind = neighbourindices[0]                                  # lowest neighbour index
                    n1ind = neighbourindices[1]                                  # other neighbours indices
                    n2ind = neighbourindices[2]
                    lowestelev = self.elev_c[n0ind]
                    if self.elev_c[neighbourindices[1]] <  lowestelev  :
                        n0ind = neighbourindices[1]
                        n1ind = neighbourindices[0]                              # other neighbours indices
                        n2ind = neighbourindices[2]
                        lowestelev = self.elev_c[n0ind]
                    if self.elev_c[neighbourindices[2]] <  lowestelev  :
                        n0ind = neighbourindices[2]
                        n1ind = neighbourindices[1]                              # other neighbours indices
                        n2ind = neighbourindices[0]
                        lowestelev = self.elev_c[n0ind]

                    # compute the plan distance from lowest neighbours CG to other neighbours CGs and associated slope
                    lxy = num.sqrt((self.coord_c[n1ind][0]-self.coord_c[n0ind][0])**2 + (self.coord_c[n1ind][1]-self.coord_c[n0ind][1])**2 )                               # dist to next neighbour in xy (plan) plane
                    s   = (self.elev_c[n1ind]-self.elev_c[n0ind])/lxy
                    if s > Rs                                        :          # unstable slope - lower other neighbours elev
                        self.elev_c[n1ind] = num.maximum(self.elev_c[n0ind]+(Rs*lxy), self.base[n1ind])
                    lxy = num.sqrt((self.coord_c[n2ind][0]-self.coord_c[n0ind][0])**2 + (self.coord_c[n2ind][1]-self.coord_c[n0ind][1])**2 )                                                           # dist to other neighbour
                    s   = (self.elev_c[n2ind]-self.elev_c[n0ind])/lxy
                    if s > Rs                                        :          # unstable slope - lower other neighbours elev
                        self.elev_c[n2ind] = num.maximum(self.elev_c[n0ind]+(Rs*lxy), self.base[n2ind])

        else:

            neighbourindices = neighbours[ind]           # get the neighbour Indices for the current triangle

            n0 = neighbourindices[:,0]
            n1 = neighbourindices[:,1]
            n2 = neighbourindices[:,2]

            k = n0.shape[0]

            import numpy
            e = numpy.zeros((k,3))

            e[:,0] = self.elev_c[n0]
            e[:,1] = self.elev_c[n1]
            e[:,2] = self.elev_c[n2]

            minid0 = numpy.argmin(e, axis=1)


            minid1 = numpy.mod(minid0+1, 3)
            minid2 = numpy.mod(minid0+2, 3)
            ident  = numpy.arange(k)

            n0ind = neighbourindices[ident,minid0]
            n1ind = neighbourindices[ident,minid1]
            n2ind = neighbourindices[ident,minid2]



            lxy = num.sqrt((self.coord_c[n1ind][:,0]-self.coord_c[n0ind][:,0])**2 + (self.coord_c[n1ind][:,1]-self.coord_c[n0ind][:,1])**2 )


            s   = (self.elev_c[n1ind]-self.elev_c[n0ind])/lxy
            self.elev_c[n1ind] = numpy.where(s>Rs, num.maximum(self.elev_c[n0ind]+(Rs*lxy), self.base[n1ind]),
                                             self.elev_c[n1ind])

            lxy = num.sqrt((self.coord_c[n2ind][:,0]-self.coord_c[n0ind][:,0])**2 + (self.coord_c[n2ind][:,1]-self.coord_c[n0ind][:,1])**2 )

            s   = (self.elev_c[n2ind]-self.elev_c[n0ind])/lxy
            self.elev_c[n2ind] = numpy.where(s>Rs, num.maximum(self.elev_c[n0ind]+(Rs*lxy), self.base[n2ind]),
                                             self.elev_c[n2ind])





        # once all triangles processed for erosion and collapse add back height to maintain volumes correctly
        self.stage_c[ind] = self.elev_c[ind] + h

        return (updated)


    def parallel_safe(self):
        """If Operator is applied independently on each cell and
        so is parallel safe.
        """

        if self.domain.get_using_discontinuous_elevation():
            return True
        else:
            return False


