avaframe/AvaFrame

View on GitHub
avaframe/com1DFA/testSPH.py

Summary

Maintainability
F
3 days
Test Coverage
F
0%

import numpy as np
import time
import math
import matplotlib.pyplot as plt

# Local imports
import avaframe.out3Plot.plotUtils as pU
import avaframe.com1DFA.DFAtools as DFAtls
import avaframe.com1DFA.com1DFA as com1DFA
# import avaframe.com1DFA.SPHfunctions as SPH
from avaframe.in3Utils import cfgUtils
import avaframe.com1DFA.DFAfunctionsCython as DFAfunC

cfg = cfgUtils.getModuleConfig(com1DFA)['GENERAL']
cfgFull = cfgUtils.getModuleConfig(com1DFA)

########################################################################
# CHOOSE YOUR SETUP
##########################################################################
# Choose the snow thickness you want to use (h function)


def Hfunction(x, y, z):
    h = np.ones(np.shape(x))
    GHx = np.zeros(np.shape(x))
    GHy = np.zeros(np.shape(x))
    # h = x*x*y/10000 + 1
    GHx = -4/(Lx*Lx)*(2*x-Lx)*np.sin(math.pi*y/Lx)
    GHy = -4/(Lx*Lx)*(x-Lx)*x*np.cos(math.pi*y/Lx)*math.pi/Lx
    h = -4/(Lx*Lx)*(x-Lx)*x*np.sin(math.pi*y/Lx)
    # r = np.sqrt((x-12.5)*(x-12.5)+(y-12.5)*(y-12.5))
    # H0 = 1
    # h = H0 * (1 - (r/12.5) * (r/12.5))
    GHz = np.zeros(np.shape(x))
    return h, GHx, GHy, GHz


# Choose the surface shape you want to use
def Sfunction(x, y, Lx, Ly):
    # plane
    Z = x*np.tan(slopeAnglex) + y*np.tan(slopeAngley)
    sx = np.tan(slopeAnglex)*np.ones(np.shape(x))
    sy = np.tan(slopeAngley)*np.ones(np.shape(x))
    area = np.sqrt(1 + (sx*sx) + (sy*sy))

    # plane
    # Z = 200-np.sqrt(200*200-x*x)
    # sx = -200*x/np.sqrt(200*200-x*x)
    # sy = np.zeros(np.shape(x))
    # area = np.sqrt(1 + (sx*sx) + (sy*sy))
    return Z, sx, sy, area


slopeAnglex = 20*math.pi/180
slopeAngley = 30*math.pi/180
# set the size of the mesh grid [m]
NDX = 5
# choose the number of particles per DX and DY
# if you choose 3, you will have 3*3 = 9 particles per grid cell
nPartPerDList = [16]

# choose if the particles should be randomly distributed.
# 0 no random part, up to 1, random fluctuation of dx/2 and dy/2
coef = 0.
rho = 200
##############################################################################
# END CHOOSE SETUP
###############################################################################


def definePart(dx, dy, Lx, Ly):
    # define particles
    nx = np.int(Lx/dx)-1
    ny = np.int(Ly/dy)-1
    nPart = nx*ny
    x = np.linspace(dx, Lx-dx, nx)
    y = np.linspace(dy, Ly-dy, ny)
    xx, yy = np.meshgrid(x, y)
    xx = xx.flatten()
    yy = yy.flatten()
    # Xpart = np.array([50., 54.])
    # Ypart = np.array([51., 53.])
    # nPart = 2
    Xpart = xx + (np.random.rand(nPart) - 0.5) * dx * coef
    Ypart = yy + (np.random.rand(nPart) - 0.5) * dy * coef
    # adding z component
    Zpart, sx, sy, area = Sfunction(Xpart, Ypart, Lx, Ly)
    Hpart, _, _, _ = Hfunction(Xpart, Ypart, Zpart)
    Mpart = Hpart * dx * dy * area * rho
    # Mpart = np.array([2000., 3000.])
    # create dictionary to store particles properties
    particles = {}
    particles['nPart'] = nPart
    particles['mTot'] = np.sum(Mpart)
    particles['x'] = Xpart
    particles['y'] = Ypart
    particles['z'] = Zpart
    particles['trajectoryLengthXY'] = np.zeros(np.shape(Ypart))
    particles['ux'] = 2*np.ones(np.shape(Ypart))
    particles['uy'] = np.ones(np.shape(Ypart))
    particles['uz'] = np.zeros(np.shape(Ypart))
    particles['m'] = Mpart
    particles['h'] = Hpart
    return particles


def defineGrid(Lx, Ly, csz):
    # define grid
    NX = np.int(Lx/csz + 1)
    NY = np.int(Ly/csz + 1)
    header = {}
    header['ncols'] = NX
    header['nrows'] = NY
    header['cellsize'] = csz
    dem = {}
    dem['header'] = header
    X = np.linspace(0, Lx, NX)
    Y = np.linspace(0, Ly, NY)
    XX, YY = np.meshgrid(X, Y)
    ZZ, _, _, _ = Sfunction(XX, YY, Lx, Ly)
    dem['rasterData'] = ZZ
    # Initialize mesh
    dem = com1DFA.initializeMesh(dem, num=4)

    Dummy = np.zeros((NY, NX))
    fields = {}
    fields['pfv'] = Dummy
    fields['ppr'] = Dummy
    fields['pft'] = Dummy
    fields['V'] = Dummy
    fields['P'] = Dummy
    fields['FT'] = Dummy
    return dem, fields


def plotFT(ax, x, xx, h, particles, ind, mark, count):
    if count == 1:
        ax.plot(xx, h, color='r', linestyle='-', label='real flow thickness')
    ax.plot(particles[x][ind], particles['h2'][ind], color='b',
             marker=mark, linestyle='None', label='corrected flow thickness N = ' + str(nPartPerDList*nPartPerDList))
    ax.plot(particles[x][ind], particles['h1'][ind], color='g',
             marker=mark, linestyle='None', label='flow thickness')
    # ax.plot(particles[x][ind], particles['h2'][ind], color='c',
    #          marker=mark, linestyle='None', label='half corrected flow thickness')
    ax.plot(particles[x][ind], particles['h'][ind], color='k',
             marker=mark, linestyle='None', label='flow thickness bilinear')
    return ax


def plotGrad(ax, x, xx, particles, ind, mark, count):
    # if count == 1:
    #     ax.plot(xx, gx, color='r', linestyle='-', label='real gradHX')
    #     ax.plot(xx, gy, color='g', linestyle='-', label='real gradHY')
    #     ax.plot(xx, gz, color='k', linestyle='-', label='real gradHZ')

    ax.plot(particles[x][ind], GHX[ind], color='m', marker=mark, markersize=5,
             linestyle='None', label='SPH N = ' + str(nPartPerDList*nPartPerDList))
    ax.plot(particles[x][ind], GHY[ind], color='m', marker=mark, markersize=5,
             linestyle='None')
    ax.plot(particles[x][ind], GHZ[ind], color='m', marker=mark, markersize=5,
             linestyle='None')
    ax.plot(particles[x][ind], GHX2[ind], color='c', marker=mark, markersize=5,
             linestyle='None', label='Corrected 1 SPH N = ' + str(nPartPerDList*nPartPerDList))
    ax.plot(particles[x][ind], GHY2[ind], color='c', marker=mark, markersize=5,
             linestyle='None')
    ax.plot(particles[x][ind], GHZ2[ind], color='c', marker=mark, markersize=5,
             linestyle='None')
    # ax.plot(particles[x][ind], GHX4[ind], color='y', marker=mark, markersize=5,
    #          linestyle='None', label='Corrected 2 SPH N = ' + str(nPartPerDList*nPartPerDList))
    # ax.plot(particles[x][ind], GHY4[ind], color='y', marker=mark, markersize=5,
    #          linestyle='None')
    # ax.plot(particles[x][ind], GHZ4[ind], color='y', marker=mark, markersize=5,
    #          linestyle='None')
    # ax.plot(particles[x][ind], GHX5[ind], color='b', marker=mark, markersize=5,
    #          linestyle='None', label='Corrected 3 SPH N = ' + str(nPartPerDList*nPartPerDList))
    # ax.plot(particles[x][ind], GHY5[ind], color='b', marker=mark, markersize=5,
    #          linestyle='None')
    # ax.plot(particles[x][ind], GHZ5[ind], color='b', marker=mark, markersize=5,
    #          linestyle='None')
    return ax


cfgPlot = cfgUtils.getModuleConfig(pU)
figW = 2.0 * float(cfg.getfloat("figW"))
figH = 2.0 * float(cfg.getfloat("figH"))

fig1, ax1 = plt.subplots(figsize=(figW, figH))
ax1.set_title('h(x)')
ax1.set_xlabel('x [m]')
ax1.set_ylabel('h [m]')
fig2, ax2 = plt.subplots(figsize=(figW, figH))
ax2.set_title('h(y)')
ax2.set_xlabel('y [m]')
ax2.set_ylabel('h [m]')
fig3, ax3 = plt.subplots(figsize=(figW, figH))
ax3.set_title('Gradh(x)')
ax3.set_xlabel('x [m]')
ax3.set_ylabel('Gradh []')
ax3.set_ylim([-0.02, 0.02])
# ax3.set_ylim([-0.016, 0.016])
fig4, ax4 = plt.subplots(figsize=(figW, figH))
ax4.set_title('Gradh(y)')
ax4.set_xlabel('y [m]')
ax4.set_ylabel('Gradh []')
# ax4.set_ylim([-0.02, 0.02])
ax4.set_ylim([-0.016, 0.016])

markers = ['o', 's', 'd', '*', 'p', 'P', '^', '>', '<', 'X', 'h']
count = 0

csz = NDX
nCells = 30

# set the extend of your mesh
Lx = nCells*csz
Ly = nCells*csz

for nPartPerD in nPartPerDList:
    dx = csz/nPartPerD
    dy = csz/nPartPerD

    # ------------------------------------------
    # define particles
    particles = definePart(dx, dy, Lx, Ly)
    Htrue = particles['h']
    # ------------------------------------------
    # define grid
    dem, fields = defineGrid(Lx, Ly, csz)

    # ------------------------------------------
    # find neighbours
    particles = DFAfunC.getNeighboursC(particles, dem)
    particles, fields = DFAfunC.updateFieldsC(cfg, particles, dem, fields)

    # ------------------------------------------
    # Compute SPH gradient
    header = dem['header']
    Nx = dem['Nx']
    Ny = dem['Ny']
    Nz = dem['Nz']
    m = particles['m']
    x = particles['x']
    y = particles['y']
    ux = particles['x']
    uy = particles['y']
    uz = particles['z']
    nx, ny, nz = DFAtls.getNormalArray(x, y, Nx, Ny, Nz, csz)
    # normal component of the velocity
    uN = ux*nx + uy*ny + uz*nz
    # print(nx, ny, nz)
    # print(norm(ux, uy, uz), uN)
    # remove normal component of the velocity
    particles['ux'] = ux - uN * nx
    particles['uy'] = uy - uN * ny
    particles['uz'] = uz - uN * nz
    indPartInCell = (particles['indPartInCell']).astype('int')
    partInCell = (particles['partInCell']).astype('int')
    indX = particles['indX'].astype('int')
    indY = particles['indY'].astype('int')

    startTime = time.time()
    # Compute sph FT
    H, C, W = DFAfunC.computeFTC(cfg, particles, header, Nx, Ny, Nz, indX, indY)
    H = np.asarray(H)
    W = np.asarray(W)
    C = np.asarray(C)
    particles['hSPH'] = (H-C)/W
    particles['W'] = W
    particles['h1'] = H
    particles['h2'] = H/W
    particles['h3'] = (H-C)/W
    tottime = time.time() - startTime
    print('Time FT: ', tottime)

    startTime = time.time()
    GHX, GHY, GHZ = DFAfunC.computeGradC(cfg, particles, header, Nx, Ny, Nz, indX, indY, SPHOption=1, gradient=0)
    GHX = np.asarray(GHX)
    GHY = np.asarray(GHY)
    GHZ = np.asarray(GHZ)
    tottime = time.time() - startTime
    print(GHX, GHY, GHZ)
    print('Time SPHOption 1: ', tottime)

    startTime = time.time()
    GHX2, GHY2, GHZ2 = DFAfunC.computeGradC(cfg, particles, header, Nx, Ny, Nz, indX, indY, SPHOption=2, gradient=0)
    GHX2 = np.asarray(GHX2)
    GHY2 = np.asarray(GHY2)
    GHZ2 = np.asarray(GHZ2)
    particles['gradx'] = GHX2
    particles['grady'] = GHY2
    particles['gradz'] = GHZ2
    tottime = time.time() - startTime
    print(GHX2, GHY2, GHZ2)
    print('Time SPHOption 2: ', tottime)

    # startTime = time.time()
    # GHX4, GHY4, GHZ4 = DFAfunC.computeGradC(cfg, particles, header, Nx, Ny, Nz, indX, indY, SPHOption=4, gradient=1)
    # GHX4 = np.asarray(GHX4)
    # GHY4 = np.asarray(GHY4)
    # GHZ4 = np.asarray(GHZ4)
    # tottime = time.time() - startTime
    # print('Time SPHOption 4: ', tottime)
    #
    # startTime = time.time()
    # GHX5, GHY5, GHZ5 = DFAfunC.computeGradC(cfg, particles, header, Nx, Ny, Nz, indX, indY, SPHOption=5, gradient=1)
    # GHX5 = np.asarray(GHX5)
    # GHY5 = np.asarray(GHY5)
    # GHZ5 = np.asarray(GHZ5)
    # tottime = time.time() - startTime
    # print('Time SPHOption 5: ', tottime)


    # ------------------------------------------------------
    # Post processing

    m = particles['m']
    x = particles['x']
    y = particles['y']
    z = particles['z']
    h, Ghx, Ghy, Ghz = Hfunction(particles['x'], particles['y'], particles['z'])

    count = count + 1
    mark = markers[count-1]
    ind = np.where(((y > Ly/4-0.5*dy) & (y < Ly/4+0.5*dy)))
    xx = np.linspace(0, Lx, 100)
    yy = Ly/4*np.ones(100)
    zz = np.zeros(100)
    h, gx, gy, gz = Hfunction(xx, yy, zz)
    ax1 = plotFT(ax1, 'x', xx, h, particles, ind, mark, count)
    ax3 = plotGrad(ax3, 'x', xx, particles, ind, mark, count)

    ind = np.where(((x > Lx/4-0.5*dx) & (x < Lx/4+0.5*dx)))
    yy = np.linspace(0, Ly, 100)
    xx = Lx/4*np.ones(100)
    zz = np.zeros(100)
    h, gx, gy, gz = Hfunction(xx, yy, zz)
    ax2 = plotFT(ax2, 'y', yy, h, particles, ind, mark, count)
    ax4 = plotGrad(ax4, 'y', yy, particles, ind, mark, count)

fig1.legend()
fig2.legend()
fig3.legend()
fig4.legend()
# Save figure to file
saveName = '_exact_semirand'
fig1.savefig(('FT_x' + saveName))
fig2.savefig(('FT_y' + saveName))
fig3.savefig(('Grad_x' + saveName))
fig4.savefig(('Grad_y' + saveName))
plt.show()