Source code for ana1Tests.FPtest

""" Flat plane test

"""

# imports
import numpy as np
import os
import logging
import matplotlib.pyplot as plt
import pathlib

# local imports
import avaframe.com1DFA.DFAtools as DFAtls
import avaframe.com1DFA.DFAfunctionsCython as DFAfunC
import avaframe.in2Trans.ascUtils as IOf
import avaframe.out3Plot.plotUtils as pU
import avaframe.in3Utils.fileHandlerUtils as fU


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


[docs]def getReleaseThickness(avaDir, cfg, demFile): """ define release thickness for Flat Plane solution test """ # Read dem demOri = IOf.readRaster(demFile) nrows = demOri['header']['nrows'] ncols = demOri['header']['ncols'] xllc = demOri['header']['xllcenter'] yllc = demOri['header']['yllcenter'] csz = demOri['header']['cellsize'] # define release thickness distribution cfgFP = cfg['FPSOL'] H0 = float(cfgFP['H0']) deltaX = float(cfgFP['deltaX']) slope = float(cfgFP['slope']) x = np.linspace(0, ncols-1, ncols)*csz+xllc y = np.linspace(0, nrows-1, nrows)*csz+yllc X, Y = np.meshgrid(x, y) r = np.sqrt((X*X)+(Y*Y)) relTh = H0 - (r-deltaX)*slope relTh = np.where(relTh < 0, 0, relTh) relTh = np.where(relTh > H0, H0, relTh) relDict = {'relTh': relTh, 'demOri': demOri, 'X': X, 'Y': Y} # save release thickness field to file relThPath = pathlib.Path(avaDir, 'Inputs', 'RELTH') fU.makeADir(relThPath) relThFile = relThPath / 'releaseThickness.asc' IOf.writeResultToAsc(demOri['header'], relTh, relThFile, flip=True) return relDict
[docs]def postProcessFPcom1DFA(cfgGen, particles, fields, ind_t, relDict): """ get fields and particles dictionaries for given time step """ demOri = relDict['demOri'] nrows = demOri['header']['nrows'] ncols = demOri['header']['ncols'] xllc = demOri['header']['xllcenter'] yllc = demOri['header']['yllcenter'] csz = demOri['header']['cellsize'] dem = relDict['dem'] x = particles['x'] y = particles['y'] ux = particles['ux'] uy = particles['uy'] m = particles['m'] h = particles['h'] force2 = {} particles, force2 = DFAfunC.computeForceSPHC(cfgGen, particles, force2, dem, cfgGen.getint('sphOption'), gradient=1) gradNorm = DFAtls.norm(force2['forceSPHX'], force2['forceSPHY'], force2['forceSPHZ']) x1, y1, z1, = DFAtls.normalize(x+xllc, y+yllc, 0*x) uMag = DFAtls.norm(ux, uy, 0) v = DFAtls.scalProd(ux, uy, 0, x1, y1, z1) grad = DFAtls.scalProd(force2['forceSPHX'], force2['forceSPHY'], force2['forceSPHZ'], x1, y1, z1) Grad = np.zeros((nrows, ncols)) MassBilinear = np.zeros((nrows, ncols)) MassBilinear = DFAfunC.pointsToRasterC(x, y, m, MassBilinear, csz=5) Grad = DFAfunC.pointsToRasterC(x, y, m*gradNorm, Grad, csz=5) indMass = np.where(MassBilinear > 0) Grad[indMass] = Grad[indMass]/MassBilinear[indMass] x = particles['x']+xllc y = particles['y']+yllc r = np.sqrt(x*x + y*y) com1DFASol = {'x': x, 'y': y, 'r': r, 'h': h, 'v': v, 'gradNorm': gradNorm, 'grad': grad, 'Grad': Grad, 'uMag': uMag, 'fields': fields} return com1DFASol
[docs]def plotProfilesFPtest(cfg, ind_time, relDict, comSol): """ Plot flow thickness and gradient for FlatPlane simulation results Parameters ----------- cfg: configparser ind_time: int time index for simiSol relDict: dict dictionary of release area info comSol: dict dictionary of simulation results and info (particles, fields, indices, time step) """ cfgGen = cfg['GENERAL'] mu = cfgGen.getfloat('mu' + cfgGen['frictModel'].lower()) cfgFP = cfg['FPSOL'] H0 = float(cfgFP['H0']) deltaX = float(cfgFP['deltaX']) slope = float(cfgFP['slope']) # get info from dem demOri = relDict['demOri'] relTh = relDict['relTh'] ncols = demOri['header']['ncols'] nrows = demOri['header']['nrows'] xllc = demOri['header']['xllcenter'] yllc = demOri['header']['yllcenter'] csz = demOri['header']['cellsize'] # com1DFA results fields = comSol['fields'] x = comSol['x'] y = comSol['y'] r = comSol['r'] h = comSol['h'] v = comSol['v'] gradNorm = comSol['gradNorm'] v = comSol['v'] outDirTest = comSol['outDirTest'] showPlot = comSol['showPlot'] Tsave = comSol['Tsave'] fig = plt.figure(figsize=(pU.figW, pU.figH)) ax1 = fig.add_subplot(121) ax2 = fig.add_subplot(122) ax1.plot(np.linspace(xllc, xllc+(ncols-1)*csz, ncols), fields['FT'][100,:], '--b', label='field flow thickness') ax1.plot(r, h, color='b', marker='.', linestyle='None', label='particle flow thickness') # ax2.plot(r, grad, color='b', marker='.', linestyle='None') ax2.plot(r, gradNorm, color='k', marker='o', linestyle='None', label='SPH gradient used') # ax2.plot(r, v, color='b', marker='.', linestyle='None') ax1.plot(np.linspace(xllc, xllc+(ncols-1)*csz, ncols), relTh[100, :], '--k') ax1.plot(r, H0-mu*(r-deltaX), '-k', label='initial expected flow thickness') ax1.set_xlabel('r in [m]') ax1.set_title('flow thickness, t=%.2f s' % (Tsave)) ax2.plot(r, mu*np.ones(np.shape(r)), '-k', label='friction threashold') ax2.set_xlabel('r in [m]') ax2.set_title('Gradient of the flow thickness') ax1.legend() ax2.legend() fig.savefig(os.path.join(outDirTest, 'radialCutSol.%s' % (pU.outputFormat))) if showPlot: plt.show() else: plt.close(fig)