avaframe/AvaFrame

View on GitHub
avaframe/ana1Tests/damBreak.py

Summary

Maintainability
B
6 hrs
Test Coverage
A
95%
"""
Simple python script to reproduce analytic solution for a Riemann problem,
following the derivations in Faccanoni and Mangeney (2012), Test 2, Case 1.2.
but scaled up in size.

Here the instantanous release of fluid from rest is described using incompressible,
thickness-avaeraged mass and momentum conservation equations and a Coulomb-tpye friction law.
"""

import numpy as np
import logging

# local imports
import avaframe.com1DFA.com1DFA as com1DFA
from avaframe.com1DFA import DFAtools
import avaframe.ana1Tests.analysisTools as anaTools
import avaframe.out3Plot.outAna1Plots as outAna1Plots

# create local logger
# change log level in calling module to DEBUG to see log messages
log = logging.getLogger(__name__)


def damBreakSol(avaDir, cfg, cfgC, outDirTest):
    """ Compute analytical Flow thickness and velocity for dam break test case

        for a granular flow over a dry rough sloping bed with the Savage Hutter model

        Parameters
        -----------
        avaDir: str
            path to avalanche directory
        cfg: configParser object
            main avaframe configuration - here showPlot flag used
        cfgC: configParser object
            configuration setting for avalanche simulation including DAMBREAK section
        outDirTest: pathlib path
            path to output directory

        Returns
        -------
        solDam: dict
            analytic solution dictionary:
                tAna: 1D numpy array
                    time array
                h0: float
                    release thickness
                hAna: 2D numpy array
                    Flow thickness (rows for x and columns for time)
                uAna: 2D numpy array
                    flow velocity (rows for x and columns for time)
                xAna: 2D numpy array
                    extent of domain in the horizontal plane coordinate system (rows for x and columns for time)
                xMidAna: 1D numpy array
                    middle of the material in x dir in the horizontal plane coordinate system
                    (used to compute the error)
    """

    # Set Parameters
    # Coordinate system chosen in the direction of the inclined plane
    g = cfgC['GENERAL'].getfloat('gravAcc')       # acceleration due to gravity [ms-2]
    phiDeg = cfgC['DAMBREAK'].getfloat('phi')       # slope angle [°]
    deltaDeg = cfgC['DAMBREAK'].getfloat('delta')       # friction angle [°]
    phi = np.radians(phiDeg)                          # slope angle [rad]
    delta = np.radians(deltaDeg)                        # bed friction angle [rad]
    gz = g * np.cos(phi)                          # projection of g perpendicular to the inclined plane
    m0 = gz * (np.tan(phi) - np.tan(delta))       # constant x-acceleration resulting from gravity and friction force
    hL = cfgC['GENERAL'].getfloat('relTh')        # initial height [m] in Riemann problem in state 1 (x<0), hR (x>0)=0
    cL = np.sqrt(gz * hL)
    x0R = cfgC['DAMBREAK'].getfloat('xBack') / np.cos(phi)
    tSave = cfgC['DAMBREAK'].getfloat('tSave')
    # Define time [0-1] seconds and space [-2,2] meters domains multiplied times 100
    tAna = np.arange(0, cfgC['DAMBREAK'].getfloat('tEnd'), cfgC['DAMBREAK'].getfloat('dt'))
    x = np.arange(cfgC['DAMBREAK'].getfloat('xStart') / np.cos(phi) , cfgC['DAMBREAK'].getfloat('xEnd') / np.cos(phi),
                  cfgC['DAMBREAK'].getfloat('dx'))
    nt = len(tAna)
    nx = len(x)
    # Initialise Flow thickness solution and velocity
    h = np.zeros((nx, nt))
    u = np.zeros((nx, nt))
    xMiddle = np.zeros(nt)

    # Compute exact solution for case: 'dry bed' - including three different states
    for m in range(nt):
        cond1R = x0R + ((m0*tAna[m]) / 2.0 - 2*cL) * tAna[m]
        cond2R = x0R + ((m0*tAna[m]) / 2.0 + cL) * tAna[m]
        cond1 = ((m0*tAna[m]) / 2.0 - cL) * tAna[m]
        cond2 = (2.0 *cL + ((m0*tAna[m]) / 2.0)) * tAna[m]
        # elif x[k] > cond2:
        u[:, m] = 0
        h[:, m] = 0
        # elif cond1 < x[k] <= cond2:
        if tAna[m] > 0:
            u[:, m] = np.where(cond2 >= x, (2./3.) * (cL + (x / tAna[m]) + m0 * tAna[m]), u[:, m])
            h[:, m] = np.where(cond2 >= x, ((2.* cL - (x / tAna[m]) + ((m0 * tAna[m]) / 2.))**2) / (9. * gz), h[:, m])
        # if x[k] <= cond1:
        u[:, m] = np.where(cond1 >= x, m0 * tAna[m], u[:, m])
        h[:, m] = np.where(cond1 >= x, hL, h[:, m])
        # uncomment this if you also want the rear part of the dam break
        # to get the analytical solution on the back part (not working very well, theory is probably wrong)
        # # elif cond1R < x[k] <= cond2R:
        # if tAna[m] > 0:
        #     u[:, m] = np.where(cond2R >= x, (2./3.) * (-cL + ((x-x0R) / tAna[m]) + m0 * tAna[m]), u[:, m])
        #     h[:, m] = np.where(cond2R >= x, ((2.* cL + ((x-x0R) / tAna[m]) -
        #                        ((m0 * tAna[m]) / 2.))**2) / (9. * gz), h[:, m])
        # # if x[k] <= cond1R:
        # u[:, m] = np.where(cond1R >= x, 0, u[:, m])
        # h[:, m] = np.where(cond1R >= x, 0, h[:, m])


        xMiddle[m] = (2*cond2R + cond1)/3

    #-----------------------------Plot results --------------
    # Reproduce figure 6, case 1.2.1 - Test 2
    outAna1Plots.plotDamAnaResults(tAna, x, xMiddle, h, u, tSave, cfg, outDirTest)

    x = x * np.cos(phi)  # projected on the horizontal plane
    xMiddle = xMiddle * np.cos(phi)  # projected on the horizontal plane
    solDam = {'tAna': tAna, 'h0': hL, 'xAna': x, 'hAna': h, 'uAna': u, 'xMidAna': xMiddle}
    return solDam


def postProcessDamBreak(avalancheDir, cfgMain, cfgDam, simDF, solDam, outDirTest):
    """ loop on all DFA simulations and compare then to the anlytic solution

    Parameters
    -----------
    avalancheDir: str or pathlib path
        avalanche directory
    cfgMain: confiparser
        avaframeCfg configuration
    cfgDam: configParser object
        configuration setting for avalanche simulation including DAMBREAK section
    simDF: pandas dataFrame
        configuration DF
    solDam: dict
        analytic solution dictionary
    outDirTest: pathlib path
        path to output directory

    Returns
    --------
    simDF: pandas dataFrame
        configuration DF appended with the analysis results
    """
    # loop on all the simulations and make the comparison to reference
    for simHash, simDFrow in simDF.iterrows():
        simName = simDFrow['simName']
        # fetch the simulation results
        fieldsList, fieldHeader, timeList = com1DFA.readFields(avalancheDir, ['FT', 'FV', 'Vx', 'Vy', 'Vz'],
                                                               simName=simName, flagAvaDir=True, comModule='com1DFA')
        # analyze and compare results
        if cfgDam['DAMBREAK']['tSave'] == '':
            indTime = -1
            tSave = timeList[-1]
        else:
            tSave = cfgDam['DAMBREAK'].getfloat('tSave')
            indTime = min(np.searchsorted(timeList, tSave), min(len(timeList)-1, len(fieldsList)-1))
            if cfgDam['DAMBREAK'].getboolean('onlyLast'):
                fieldsList = [fieldsList[0], fieldsList[indTime]]
                timeList = [0, timeList[indTime]]
                indTime = 1

        # computeAndPlotGradient(avalancheDir, particlesList, timeList, solDam, cfgDam, outDirTest, simHash, simDFrow)
        hEL2Array, hELMaxArray, vhEL2Array, vhELMaxArray, t = analyzeResults(avalancheDir, fieldsList, timeList, solDam,
                                                                             fieldHeader, cfgDam, outDirTest, simHash,
                                                                             simDFrow)
        # add result of error analysis
        # save results in the simDF
        simDF.loc[simHash, 'timeError'] = t
        simDF.loc[simHash, 'hErrorL2'] = hEL2Array[indTime]
        simDF.loc[simHash, 'vhErrorL2'] = vhEL2Array[indTime]
        simDF.loc[simHash, 'hErrorLMax'] = hELMaxArray[indTime]
        simDF.loc[simHash, 'vhErrorLMax'] = vhELMaxArray[indTime]
        # +++++++++POSTPROCESS++++++++++++++++++++++++
        # -------------------------------
    name = 'results' + str(round(t)) + '.p'
    simDF.to_pickle(outDirTest / name)

    return simDF


def analyzeResults(avalancheDir, fieldsList, timeList, solDam, fieldHeader, cfg, outDirTest, simHash, simDFrow):
    """Compare analytical and com1DFA results, compute error

        Parameters
        -----------
        avalancheDir: pathlib path
            path to avalanche directory
        fieldsList: list
            list of fields dictionaries
        timeList: list
            list of time steps
        solDam: dict
            analytic solution dictionary
        fieldHeader: dict
            header dictionary with info about the extend and cell size
        cfg: configParser object
            configuration setting for avalanche simulation including DAMBREAK section
        outDirTest: pathlib path
            path output directory (where to save the figures)
        simHash: str
            com1DFA simulation id
        simDFrow: pandas object
            com1DFA simulation rox coresponding to simHash

        Returns
        --------
        hErrorL2Array: numpy array
            L2 error on Flow thickness for saved time steps
        hErrorLMaxArray: numpy array
            LMax error on Flow thickness for saved time steps
        vErrorL2Array: numpy array
            L2 error on flow velocity for saved time steps
        vErrorLMaxArray: numpy array
            LMax error on flow velocity for saved time steps
        tSave: float
            time corresponding the errors

    """
    cfgDam = cfg['DAMBREAK']
    phi = cfgDam.getfloat('phi')
    phiRad = np.radians(phi)
    cellSize = fieldHeader['cellsize']
    xMiddle = solDam['xMidAna']
    tAna = solDam['tAna']
    xAna = solDam['xAna']
    hAna = solDam['hAna']
    uAna = solDam['uAna']
    # get plots limits
    limits = outAna1Plots.getPlotLimits(cfgDam, fieldsList, fieldHeader)
    hErrorL2Array = np.zeros((len(fieldsList)))
    vhErrorL2Array = np.zeros((len(fieldsList)))
    hErrorLMaxArray = np.zeros((len(fieldsList)))
    vhErrorLMaxArray = np.zeros((len(fieldsList)))
    count = 0
    # run the comparison routine for each saved time step
    for t, field in zip(timeList, fieldsList):
        # get similartiy solution h, u at required time step
        indTime = np.searchsorted(tAna, t)
        xM = xMiddle[indTime]
        # get extent where the sol should be compared
        xDamPlus, nColMid, nColMax, nRowMin, nRowMax = getDamExtend(fieldHeader, xM, cfgDam)
        # get analytical solution (in same format as the simulation results)
        hDamPlus = np.interp(xDamPlus, xAna, hAna[:, indTime])
        uDamPlus = np.interp(xDamPlus, xAna, uAna[:, indTime])
        # get numirical sol on good interval
        hNumPlus = field['FT'][nRowMin:nRowMax, nColMid:nColMax]
        dataNumV = DFAtools.scalProd(field['Vx'], field['Vy'], field['Vz'], np.cos(phiRad), 0, -np.sin(phiRad))
        uNumPlus = dataNumV[nRowMin:nRowMax, nColMid:nColMax]

        # make analytical results 2D
        nrows, ncols = np.shape(hNumPlus)
        O = np.ones((nrows, ncols))
        hDamPlus = O*hDamPlus
        uDamPlus = O*uDamPlus

        hEL2Plus, hL2RelPlus, hLmaxPlus, hLmaxRelPlus = anaTools.normL2Scal(hDamPlus, hNumPlus, cellSize,
                                                                            np.cos(phiRad))
        vhL2Plus, vhL2RelPlus, vhLmaxPlus, vhLmaxRelPlus = anaTools.normL2Scal(hDamPlus*uDamPlus, hNumPlus*uNumPlus,
                                                                               cellSize, np.cos(phiRad))
        if cfgDam.getboolean('relativError'):
            hErrorL2Array[count] = hL2RelPlus
            hErrorLMaxArray[count] = hLmaxRelPlus
            vhErrorL2Array[count] = vhL2RelPlus
            vhErrorLMaxArray[count] = vhLmaxRelPlus
        else:
            hErrorL2Array[count] = hEL2Plus
            hErrorLMaxArray[count] = hLmaxPlus
            vhErrorL2Array[count] = vhL2Plus
            vhErrorLMaxArray[count] = vhLmaxPlus
        title = outAna1Plots.getTitleError(cfgDam.getboolean('relativError'))
        log.debug("L2 %s error on the Flow thickness at t=%.2f s is : %.4f" % (title, t, hEL2Plus))
        log.debug("L2 %s error on the momentum at t=%.2f s is : %.4f" % (title, t, vhL2Plus))
        # Make all individual time step comparison plot
        if cfgDam.getboolean('plotSequence'):
            outAna1Plots.plotComparisonDam(cfgDam, simHash, fieldsList[0], field, fieldHeader, solDam, t, limits,
                                           outDirTest)
        count = count + 1

    tSave = cfgDam.getfloat('tSave')
    indT = min(np.searchsorted(timeList, tSave), min(len(timeList)-1, len(fieldsList)-1))
    tSave = timeList[indT]
    if cfgDam.getboolean('plotErrorTime') and len(timeList)>1:
        # Create result plots

        outAna1Plots.plotErrorTime(timeList, hErrorL2Array, hErrorLMaxArray, vhErrorL2Array, vhErrorLMaxArray,
                                   cfgDam.getboolean('relativError'), tSave, simHash, outDirTest)

    outAna1Plots.plotDamBreakSummary(avalancheDir, timeList, fieldsList, fieldHeader, solDam, hErrorL2Array,
                                     hErrorLMaxArray, vhErrorL2Array, vhErrorLMaxArray, outDirTest, simDFrow, simHash,
                                     cfg)

    return hErrorL2Array, hErrorLMaxArray, vhErrorL2Array, vhErrorLMaxArray, tSave


def getDamExtend(fieldHeader, xM, cfgDam):
    """ Get the extend where the analytic and simulation should be compared

    Parameters
    -----------
    fieldHeader: dict
        header dictionary with info about the extend and cell size
    xM: float
        x coordinate ot the start of the comparizon domain
    cfgDam: configParser object
        configuration setting for the DAMBREAK section

    Returns
    --------
    xDamPlus: 1D numpy array
        x array corresponding to the comparizon domain
    nColMid: int
        index of the column corresponding to xM
    nColMax: int
        index of the column corresponding to xEnd
    nRowMin: int
        index of the row corresponding to yStart
    nRowMax: int
        index of the row corresponding to yEnd
    """
    cellSize = fieldHeader['cellsize']
    xllc = fieldHeader['xllcenter']
    yllc = fieldHeader['yllcenter']
    nCols = fieldHeader['ncols']
    xEnd = cfgDam.getfloat('xEnd')
    yStart = cfgDam.getfloat('yStart')
    yEnd = cfgDam.getfloat('yEnd')
    xArray = np.arange(nCols)*cellSize + xllc
    nColMid = round((xM-xllc)/cellSize)
    nColMax = round((xEnd-xllc)/cellSize)
    nRowMin = round((yStart-yllc)/cellSize)
    nRowMax = round((yEnd-yllc)/cellSize)

    xDamPlus = xArray[nColMid:nColMax]

    return xDamPlus, nColMid, nColMax, nRowMin, nRowMax