# -*- coding: utf-8 -*-
"""
"""
#==========================================================================
# Vector field tools.
#==========================================================================
# ----------------------------------------------------------------------------
# IMPORTS
# ----------------------------------------------------------------------------
# Standard Python Dependencies
import os
import sys
import time
import scipy.io as sio
import numpy as np
# Non-Standard Python Dependencies
# Local Module Dependencies
from fastwater.general.fileTools import generateGlobalAttributes
from fastwater.telemac.extraction import Date2Num, Date2Str, genTimeStamps, default_epoch
from fastwater.telemac.extraction import getNodesWeights, enhanceMesh, slfObj
from fastwater.general.similitude import courantNumber, froudeNumber, ekmanNumber, rossbyNumber
from fastwater.general.meshTools import getNodesAroundNode, gradientAtNode, weightsMatrix
# Other Dependencies
#--------------------------------------------------------------------------
# GLOBAL CONSTANTS
#--------------------------------------------------------------------------
# ----------------------------------------------------------------------------
# CLASS DEFINITIONS
# ----------------------------------------------------------------------------
#--------------------------------------------------------------------------
# FUNCTION DEFINITIONS
#--------------------------------------------------------------------------
[docs]def initialiseVectorField(slf):
slf.set_mpl_tri(True)
slf.set_kd_tree(True)
times = slf.tags['times']
idate = slf.datetime
nsteps = len(times)
nnodes = slf.npoin2
meshx = slf.meshx
meshy = slf.meshy
elems = slf.ikle3
edges = slf.edges #getEdgesSLF(slf.ikle3,slf.meshx,slf.meshy,False)
return times,idate,nsteps,nnodes,meshx,meshy,elems,edges
[docs]def getVariableAtNodes(slf,nodes,varIndex):
var = np.squeeze(slf.get_series(nodes+1,[varIndex],False))
return var
[docs]def getVariablesAtTime(slf,frame,varIndex):
vars = slf.getVariablesAt(frame,varIndex)
return vars
[docs]def vectorFieldAtLocation(slf,Loc,d_off=0.0,rho=1025.):
NNghb, NWghts, Locs = getNodesWeights(slf,Loc[0],Loc[1])
npts = len(NNghb)
times,idate,nsteps,nnodes,meshx,meshy,elems,edges = initialiseVectorField(slf)
elev = np.zeros([nsteps,npts])
velx = np.zeros([nsteps,npts])
vely = np.zeros([nsteps,npts])
dpth = np.zeros([nsteps,npts])
dhdx = np.zeros([nsteps,npts])
dhdy = np.zeros([nsteps,npts])
dudx = np.zeros([nsteps,npts])
dudy = np.zeros([nsteps,npts])
dvdx = np.zeros([nsteps,npts])
dvdy = np.zeros([nsteps,npts])
for inode in np.arange(npts):
node = NNghb[inode]
nodes = getNodesAroundNode(node,edges)
nodesX = meshx[nodes]
nodesY = meshy[nodes]
indx, = np.where(nodes == node)
u = getVariableAtNodes(slf,nodes,0)
v = getVariableAtNodes(slf,nodes,1)
h = getVariableAtNodes(slf,nodes,3)
b = getVariableAtNodes(slf,nodes,4)
dpth[:,inode] = d_off - b[indx,:]
elev[:,inode] = h[indx,:]
velx[:,inode] = u[indx,:]
vely[:,inode] = v[indx,:]
for istep in np.arange(nsteps):
#grad,res,rnk,sng = gradientAtNode(h,meshx,meshy,node,nodes,istep)
val = h[:,istep]
grad,res,rnk,sng = gradientAtNode(node,nodes,nodesX,nodesY,val)
dhdx[istep,inode] = grad[0]
dhdy[istep,inode] = grad[1]
#grad,res,rnk,sng = gradientAtNode(u,meshx,meshy,node,nodes,istep)
val = u[:,istep]
grad,res,rnk,sng = gradientAtNode(node,nodes,nodesX,nodesY,val)
dudx[istep,inode] = grad[0]
dudy[istep,inode] = grad[1]
#grad,res,rnk,sng = gradientAtNode(v,meshx,meshy,node,nodes,istep)
val = v[:,istep]
grad,res,rnk,sng = gradientAtNode(node,nodes,nodesX,nodesY,val)
dvdx[istep,inode] = grad[0]
dvdy[istep,inode] = grad[1]
del node, nodes, u, v, h
vort = dvdx - dudy
# Interpolate onto location of interest
wgts = weightsMatrix(NWghts,nsteps)
epochStr = Date2Str(default_epoch)
vect = {}
vect['loc'] = Loc
vect['nearestNodes'] = NNghb
vect['nodeWeights'] = NWghts
vect['nodeLocs'] = Locs
vect['timeUnits'] = 'days since '+epochStr
vect['epoch'] = epochStr
vect['times'] = genTimeStamps(idate, times, nsteps)
vect['elev'] = elev
vect['dpth'] = dpth
vect['U'] = velx
vect['V'] = vely
vect['dhdx'] = dhdx
vect['dhdy'] = dhdy
vect['dudx'] = dudx
vect['dudy'] = dudy
vect['dvdx'] = dvdx
vect['dvdy'] = dvdy
vect['vort'] = vort
vect['elev_loc'] = np.sum(elev.transpose()*wgts,0)
vect['dpth_loc'] = np.sum(dpth.transpose()*wgts,0)
vect['U_loc'] = np.sum(velx.transpose()*wgts,0)
vect['V_loc'] = np.sum(vely.transpose()*wgts,0)
vect['dhdx_loc'] = np.sum(dhdx.transpose()*wgts,0)
vect['dhdy_loc'] = np.sum(dhdy.transpose()*wgts,0)
vect['dudx_loc'] = np.sum(dudx.transpose()*wgts,0)
vect['dudy_loc'] = np.sum(dudy.transpose()*wgts,0)
vect['dvdx_loc'] = np.sum(dvdx.transpose()*wgts,0)
vect['dvdy_loc'] = np.sum(dvdy.transpose()*wgts,0)
vect['vort_loc'] = np.sum(vort.transpose()*wgts,0)
return vect
[docs]def characteriseAtNode(slf,node):
times,idate,nsteps,nnodes,meshx,meshy,elems,edges = initialiseVectorField(slf)
Cr = np.zeros([nsteps,])
Fr = np.zeros([nsteps,])
Ek = np.zeros([nsteps,])
Ro = np.zeros([nsteps,])
nodes = getNodesAroundNode(node,edges)
nodesX = meshx[nodes]
nodesY = meshy[nodes]
indx, = np.where(nodes == node)
u = getVariableAtNodes(slf,nodes,0)
v = getVariableAtNodes(slf,nodes,1)
d = getVariableAtNodes(slf,nodes,2)
dpth = np.squeeze(d[indx,:])
velx = np.squeeze(u[indx,:])
vely = np.squeeze(v[indx,:])
vel = np.sqrt(np.square(velx)+np.square(vely))
for istep in np.arange(nsteps):
val = u[:,istep]
grad,res,rnk,sng = gradientAtNode(node,nodes,nodesX,nodesY,val)
dudy = grad[1]
val = v[:,istep]
grad,res,rnk,sng = gradientAtNode(node,nodes,nodesX,nodesY,val)
dvdx = grad[0]
omega = 0.5 * (dvdx - dudy)
Cr[istep] = courantNumber(node, nodes, nodesX, nodesY, velx[istep], vely[istep], 10.0)
Fr[istep] = froudeNumber(dpth[istep], vel[istep])
Ek[istep] = ekmanNumber(dpth[istep], omega)
Ro[istep] = rossbyNumber(dpth[istep], omega)
# Store results in dictionary structure
vect = {}
vect['node'] = node
vect['nearestNodes'] = nodes
vect['nodesX'] = nodesX
vect['nodesY'] = nodesY
vect['times'] = np.reshape(Date2Num(idate) + (times / (24.*60.0*60.)),[nsteps,1])
vect['Cr'] = Cr
vect['Fr'] = Fr
vect['Ek'] = Ek
vect['Ro'] = Ro
return vect
[docs]def getVelocityField(slfHnd, dispFreq=5000):
# Construct file name parts
fparts = slfHnd.file['name'].replace('\\','/').split('/')
pidx = fparts.index('RESULTS')
fpath = "/".join(fparts[0:pidx])+'/VELOCITY'
fname = fparts[-1].split('.')[0]
# Get data dimensions
nNode = slfHnd.npoin2
nElem = slfHnd.nelem3
elems = slfHnd.ikle3.copy()
nStep = len(slfHnd.tags['times'])
# Initialise data structures
print('Initializing variables...')
velx_node = np.zeros([nNode,nStep],dtype=np.single)
vely_node = np.zeros([nNode,nStep],dtype=np.single)
vmag_node = np.zeros([nNode,nStep],dtype=np.single)
# Load surface elevation data
print('Loading velocity data...')
for tidx in np.arange(nStep, dtype=np.int64):
velx_node[:,tidx] = slfHnd.get_variables_at(tidx,[0])[0]
vely_node[:,tidx] = slfHnd.get_variables_at(tidx,[1])[0]
vmag_node[:,tidx] = np.sqrt(np.square(velx_node[:,tidx])+
np.square(vely_node[:,tidx]))
# Construct standardised time stamps
epoch = Date2Num(slfHnd.datetime)
times = slfHnd.tags['times']
t = np.reshape(epoch + (times / (24.*60.0*60.)),[nStep,1])
# Extract mesh
print('Extracting model mesh...')
bathy = slfHnd.get_variables_at(0,[4])[0]
mesh = {}
mesh['times'] = t.copy()
mesh['nNodes'] = nNode
mesh['meshx'] = slfHnd.meshx.copy()
mesh['meshy'] = slfHnd.meshy.copy()
mesh['nElems'] = nElem
mesh['elems'] = elems.copy()
mesh['bathy'] = bathy.copy()
mesh = enhanceMesh(mesh)
cWghts = mesh['cWghts'].copy()
sio.savemat(os.path.join(fpath,fname+"_mesh.mat"),mesh)
for istp in np.arange(nStep, dtype=np.int64):
oname = fname+"_velocities_nodes_step"+str(istp).zfill(4)+".mat"
vel_nodes = {}
vel_nodes['velx'] = velx_node[:,istp]
vel_nodes['vely'] = vely_node[:,istp]
vel_nodes['vmag'] = vmag_node[:,istp]
sio.savemat(os.path.join(fpath,oname),vel_nodes)
# Interploate elevation residuals onto element centres
velx_cntr = np.zeros([nElem,nStep],dtype=np.single)
vely_cntr = np.zeros([nElem,nStep],dtype=np.single)
vmag_cntr = np.zeros([nElem,nStep],dtype=np.single)
print('Interpolate onto cell centres...')
print('Number of cells to process = ',nElem)
for ielm in np.arange(nElem, dtype=np.int64):
cw = cWghts[ielm,:]
v = velx_node[elems[ielm,:],:]
velx_cntr[ielm,:] = np.sum((v.transpose()*cw),1)
v = vely_node[elems[ielm,:],:]
vely_cntr[ielm,:] = np.sum((v.transpose()*cw),1)
v = vmag_node[elems[ielm,:],:]
vmag_cntr[ielm,:] = np.sum((v.transpose()*cw),1)
if not(np.mod(ielm+1,dispFreq)):
sys.stdout.write("\r"+"cells processed = "+str(ielm+1))
sys.stdout.flush()
print("\r"+"cells processed = "+str(nElem), flush=True)
for istp in np.arange(nStep, dtype=np.int64):
oname = fname+"_velocities_cntrs_step"+str(istp).zfill(4)+".mat"
vel_cntrs = {}
vel_cntrs['velx'] = velx_cntr[:,istp]
vel_cntrs['vely'] = vely_cntr[:,istp]
vel_cntrs['vmag'] = vmag_cntr[:,istp]
sio.savemat(os.path.join(fpath,oname),vel_cntrs)
print('Process complete.')
return fname
[docs]def processNodeSet(fPath,fName,nCores,setNum):
slf = slfObj(os.path.join(fPath,fName))
times,idate,nsteps,nnodes,meshx,meshy,elems,edges = initialiseVectorField(slf)
t = np.reshape(Date2Num(idate) + (times / (24.*60.0*60.)),[nsteps,1])
npts = np.int64(np.floor(slf.npoin3/nCores))
node0 = setNum*npts
# Account for extra points to last core
if np.mod(slf.npoin3,nCores) != 0:
if (setNum+1) == nCores:
npts = npts + np.mod(slf.npoin3,nCores)
print(' Number of Nodes to Process: '+str(npts))
print(' First Node in Set: '+str(node0))
matFile = slf.file['name'].split('/')[-1].split('.')[0]+'_gradients_'+'SET'+str(setNum).zfill(4)+'.mat'
elev = np.zeros([nsteps,npts])
velx = np.zeros([nsteps,npts])
vely = np.zeros([nsteps,npts])
dpth = np.zeros([nsteps,npts])
dhdx = np.zeros([nsteps,npts])
dhdy = np.zeros([nsteps,npts])
dudx = np.zeros([nsteps,npts])
dudy = np.zeros([nsteps,npts])
dvdx = np.zeros([nsteps,npts])
dvdy = np.zeros([nsteps,npts])
tic = time.time()
for inode in np.arange(npts):
node = node0 + inode
nodes = getNodesAroundNode(node,edges)
nodesX = meshx[nodes]
nodesY = meshy[nodes]
indx = np.where(nodes == node)
u = getVariableAtNodes(slf,nodes,0)
v = getVariableAtNodes(slf,nodes,1)
h = getVariableAtNodes(slf,nodes,3)
b = getVariableAtNodes(slf,nodes,4)
dpth[:,inode] = b[indx,:]
elev[:,inode] = h[indx,:]
velx[:,inode] = u[indx,:]
vely[:,inode] = v[indx,:]
for istep in np.arange(nsteps):
#grad,res,rnk,sng = gradientAtNode(h,meshx,meshy,node,nodes,istep)
grad,res,rnk,sng = gradientAtNode(node,nodes,nodesX,nodesY,h[:,istep])
dhdx[istep,inode] = grad[0]
dhdy[istep,inode] = grad[1]
#grad,res,rnk,sng = gradientAtNode(u,meshx,meshy,node,nodes,istep)
grad,res,rnk,sng = gradientAtNode(node,nodes,nodesX,nodesY,u[:,istep])
dudx[istep,inode] = grad[0]
dudy[istep,inode] = grad[1]
#grad,res,rnk,sng = gradientAtNode(v,meshx,meshy,node,nodes,istep)
grad,res,rnk,sng = gradientAtNode(node,nodes,nodesX,nodesY,v[:,istep])
dvdx[istep,inode] = grad[0]
dvdy[istep,inode] = grad[1]
del node, nodes, u, v, h
vort = dvdx - dudy
toc = time.time() - tic
print('Time Elasped = ',toc,' s')
dict = {}
dict['times'] = t
dict['nsteps'] = nsteps
dict['nodes'] = node0+np.arange(npts)
dict['meshx'] = meshx[node0+np.arange(npts)].copy()
dict['meshy'] = meshy[node0+np.arange(npts)].copy()
dict['elev'] = elev
dict['dpth'] = dpth
dict['U'] = velx
dict['V'] = vely
dict['dhdx'] = dhdx
dict['dhdy'] = dhdy
dict['dudx'] = dudx
dict['dudy'] = dudy
dict['dvdx'] = dvdx
dict['dvdy'] = dvdy
dict['vort'] = vort
extractPath = fPath.replace('RESULTS','EXTRACT')
sio.savemat(os.path.join(extractPath,matFile),dict)
return
[docs]def processLocation(filePath,runStr,suffix,modelVer,Loc,depthOffset:float=0.0):
fName = runStr+'_2D.slf'
slf = slfObj(os.path.join(filePath,fName))
vect = vectorFieldAtLocation(slf,Loc,d_off=depthOffset)
glbattr = generateGlobalAttributes('Telemac - Extraction of dynamics at a location')
vect['globalAttributes'] = glbattr
mvStr = 'v'+str(modelVer).zfill(2)
matFile = runStr+'_'+suffix+'_DYN_'+mvStr+'.mat'
extractPath = filePath.replace('RESULTS','EXTRACT')
sio.savemat(os.path.join(extractPath,matFile),vect)
return vect