Entry 10 (Second Place)

../../_images/entry10.png

Authors

  • Kristen Thyng

Description

This figure is from a paper for the European Wave and Tidal Energy Conference (EWTEC) that I am working on with a co-author. The paper describes a study we did together combining a simulation of an idealized tidal channel with a symmetric headland in the middle and an array of ten turbines near the headland. This combination of efforts allows us to start to understand the effects of turbines on increasingly realistic flows, both near the turbine array and farther field, kilometers away down the channel.

The plot shows the difference in maximum turbulent kinetic energy (TKE) between a simulation with no turbines and a simulation with an array of turbines. Positive (red) values indicate areas in which the initial case has larger maximum values and negative (grey) values indicate areas in which the turbine array case has larger values. Line plots to the left (bottom) of the main plot area show the along- (across-) channel average of the plot values, with coloring indicating position greater than (red) or less than (grey) zero.

There is a distinction in relative behavior in the cases in the near- and far-field. Very near to the turbines and in their wake, the TKE is stronger in the turbine case, which is expected given that several terms have been added to the turbulence equations in order to model the effect of the turbines. Further downstream, the initial, no turbine case TKE is larger than the turbine case, probably due to the dissipation due to the turbines in the regular array case. These trends are reflected in the along- and across-channel averages shown alongside the plot: the initial case tends to be larger except near the turbines themselves.

Products

Source

import plot_gen as pg
import numpy as np

# # Some subtraction plots
# # vorticity (da) 
# # levels = np.linspace(-.09,.09,10)
# smax: np.linspace(-.9,.9,10), smean: np.linspace(-.45,.45,10)
# zetamax: np.linspace(-.01,.01,9)
var = {'vort':{'levels':np.linspace(-.0045,.0045,10),'colormap':'RdGy_r','title':r'$\Delta$Max vertical vorticity $(1/s)$'},
	   'tke':{'levels':np.linspace(-.09,.09,10),'colormap':'RdGy_r','title':r'$\Delta$Max TKE ($m^2/s^2$)'},
	   'w':{'levels':np.linspace(-.9,.9,10),'colormap':'RdGy_r','title':r'$\Delta$Max vertical velocity $(m/s)$'},
	   's':{'levels':np.linspace(-.9,.9,10),'colormap':'RdGy_r','title':r'$\Delta$Max horizontal speed $(m/s)$'},
	   'zeta':{'levels':np.linspace(-.0001,.0001,9),'colormap':'RdGy_r','title':r'$\Delta$Mean zeta $(m)$'},
	   'mkpd':{'levels':np.linspace(-.9,.9,10),'colormap':'RdGy_r','title':r'$\Delta$ mean kinetic power density $(kW/m^2)$'},
	   'abi':{'levels':np.linspace(-27,27,10),'colormap':'RdGy_r','title':r'$\Delta$ asymmetry (degrees)'},
	   'stdbi':{'levels':np.linspace(-27,27,10),'colormap':'RdGy_r','title':r'$\Delta$ directional deviation (degrees)'}}

name = 'tke' # when i want to do more do a loop
oper = 'max' # operation (max, mean, metric, min,gradmean)
zoom = 'zmid'
# vlevel = 'hubheight' #'surface'#'hubheight' (depth-averaged)
data = np.load('diff.npz')
diff = data['diff']
x = data['x']
y = data['y']
mask = data['mask']
data.close()

### TEMP
# ind = np.arange(0,50)#15,15+25) # half the cycle

# diff = d.max(axis=0)-dr.max(axis=0)

pg.plt(x/1000.,y/1000.,diff,mask,zoom,
	var[name]['colormap'],var[name]['levels'],
	name + oper + '.pdf',
	var[name]['title'])
'''
troc_plot_gen.py
Kristen M. Thyng
November 2012

Plot properties that were saved and calculated in troc_calc.py, at 
multiple depths. More general than troc_plot.py.

Movies:
	Density
	U Velocity
	V Velocity
	TKE (5 m down)
	Vertical Velocity (5 m down)
	Free Surface
Plots:
	Asymmetry properties:
		bidirectionality
		directional deviation
	Vorticity (need to calculate in sigma coordinates)
	Mean speed
	Mean Density
	Mean TKE
	Mean w
	Mean vorticity
	Mean kinetic power density
	Mean free surface
	Free surface signal at edge of domain

'''

# import matplotlib
# matplotlib.use("Agg") # set matplotlib to use the backend that does not require a windowing system
import os
import numpy as np
from matplotlib.pyplot import *
import matplotlib.colors as colors
from matplotlib import rc, ticker, cm
import troc_calc as tr
import pdb
import op
import netCDF4 as netCDF
from pylab import *
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes
from mpl_toolkits.axes_grid1.inset_locator import mark_inset

#rc('text', usetex=True) # for latex rendering
##rc('font', family='serif')

# # User stuff
# var = 'mssurface'
# vtype = 'calculate'
# zoom = 'zout'
# xlims = (0,40)
# ylims = (0,7)

hgrey = '#737373'
hred = '#c94741'

def plt(x,y,z,mask,zoom,colormap,levels,fname,Title):

	if zoom == 'zout':
		xlims = (0,40)
		ylims = (0,7)
		dy = 24
		dx = 18
		ms = 1.5
	elif zoom == 'zmid':
		xlims = (10,30)
		ylims = (0,7)
		dy = 24
		dx = 14
		ms = 1.5 # min shaft length
	elif zoom == 'zin':
		xlims = (17,23)
		ylims = (3,5)
		dy = 12#6
		dx = 8#4
		ms = .25

	# Make plots
	# rc('text',usetex=True)
	matplotlib.rcParams.update({'font.size': 16})
	# if vtype == 'calculate':
		# if var == 'abisurface' or var == 'abihubheight' or var == 'stdbisurface' or var == 'stdbihubheight':
		# 	figure(figsize=(18,6))
		# 	contourf(x,y,d,[5,10,15,20,25,30,35,40],cmap=colormap)#,vmin=dmin,vmax=dmax)
		# 	xlim(xlims)
		# 	ylim(ylims)
		# 	colorbar()
		# 	xlabel('Along-channel distance (km)')
		# 	ylabel('Across-channel distance (km)')
		# 	title(Title)
		# 	savefig('figures/calculate/' + var + zoom + '.png',bbox_inches='tight')
		# 	close()
		# else:
	# figure(figsize=(18,6))
	
	# definitions for axes, a la http://matplotlib.org/examples/pylab_examples/scatter_hist.html
	left, bottom = 0.05, 0.1 # for very left and very bottom
	thickx, thicky = 0.03, 0.125 # thickness of the left and bottom plots
	width, height = 0.75, 0.6 # width and height of main plot
	dpx = .025 # space between plots
	dpy = .08

	# axes for each of the three plots: [left,bottom,width,height]
	rect_main = [left+thickx+dpx,bottom+thicky+dpy,width+.153,height] # main plot
	rect_left = [left,bottom+thicky+dpy,thickx,height] # left side small plot
	rect_bottom = [left+thickx+dpx,bottom,width,thicky] # bottom small plot
	nullfmt   = NullFormatter()         # no labels
	nullloc   = NullLocator()         

	# start with a rectangular figure
	figure(1,figsize=(18,10))#18,6))

	axMain = axes(rect_main)

	# main plot

	# pcolormesh(x,y,z,cmap=colormap)#,vmin=dmin,vmax=dmax)
	# cmap = cm.PRGn
	# levels=[-0.08      , -0.06222222, -0.04444444, -0.02666667, -0.00888889,
 #        0.00888889,  0.02666667,  0.04444444,  0.06222222,  0.08]
 	# dmax = np.max(z.max(),abs(z.min()))
	# this is to prevent the whiting out of larger values in the contour maps
	dmax = np.max(levels)
	ind = (z>=dmax)
	z[ind] = dmax
	ind = (z<=-dmax)
	z[ind] = -dmax
	# pdb.set_trace()
	contourf(x,y,z,levels,cmap=get_cmap(colormap))#levels=np.linspace(-dmax,dmax,10))#,levels=np.linspace(-.09,.09,10))#cmap=cm.get_cmap(colormap, len(levels)-1))#cmap=colormap)#,vmin=dmin,vmax=dmax#levels=np.linspace(dmin,dmax,9))
	# pcolormesh(x,y,z,cmap=colormap,vmin=-dmax,vmax=dmax)
	colorbar(pad=.02)
	#set_cmap(colormap)

	axMain.set_xlim(xlims)
	axMain.set_ylim(ylims)
	axMain.set_title(Title)
	axMain.set_ylabel('Across-channel distance (km)')
	axMain.set_xlabel('Along-channel distance (km)')
	# axMain.xaxis.set_major_formatter(nullfmt) # no labels
	# axMain.yaxis.set_major_formatter(nullfmt) # no labels

	# Draw in headland to distinguish from the flow in whiter plots
	hold('on')
	contour(x,y,mask,levels=[0],colors='grey',alpha=.5)
	# ind = (mask == 0)
	# pcolor(x,y,ind*.5)

	# Try including various averages along the edges of the plots: left and bottom
	# Left plot: averaged along-channel
	axLeft = axes(rect_left,frameon=False)#,sharey=axMain)
	ind = ~np.isnan(z)
	# pdb.set_trace()
	# zacross = z.mean(axis=1)
	# poor man's nanmean
	zacross = np.nansum(z,axis=1)/ind.sum(axis=1)
	indp = (zacross>0)
	indn = (zacross<=0)
	zacrossp = zacross.copy()
	zacrossp[indn] = np.nan
	zacrossn = zacross.copy()
	zacrossn[indp] = np.nan
	axLeft.plot(zacrossp,y[:,0],color=hred,alpha=1,linewidth=2)
	hold('on')
	# ind = (zacross<=0)
	axLeft.plot(zacrossn,y[:,0],color=hgrey,alpha=.7,linewidth=2)
	# axLeft.plot(np.zeros(y.shape[0]),y[:,0],':k')
	xmax = max(zacross.max(),abs(zacross.min()))
	axLeft.set_xlim((xmax,-xmax))	
	axLeft.set_ylim(ylims)
	axLeft.set_xlabel('\n\n Channel\n Directional\n averages',color=hgrey)
	# locs,labels = xticks()
	# setp(labels,visible=False)#rotation=90,color='grey')
	# setp(locs,visible=False)#rotation=90,color='grey')
	axLeft.yaxis.set_major_locator(nullloc) # no labels
	axLeft.xaxis.set_major_locator(nullloc) # no labels

	# Bottom plot: averaged across-channel
	axBottom = axes(rect_bottom,frameon=False)#,sharex=axMain)
	zalong = z.mean(axis=0)
	indp = (zalong>0)
	indn = (zalong<=0)
	zalongp = zalong.copy()
	zalongp[indn] = np.nan
	zalongn = zalong.copy()
	zalongn[indp] = np.nan
	axBottom.plot(x[0,:],zalongn,color=hgrey,alpha=.7,linewidth=2)
	hold('on')
	axBottom.plot(x[0,:],zalongp,color=hred,alpha=1,linewidth=2)
	# axBottom.plot(x[0,:],np.zeros(x.shape[1]),':k')
	axBottom.set_xlim(xlims)
	ymax = max(zalong.max(),abs(zalong.min()))
	axBottom.set_ylim((-ymax,ymax))	
	axBottom.yaxis.set_major_locator(nullloc) # no labels
	axBottom.xaxis.set_major_locator(nullloc) # no labels
	# axBottom.yaxis.tick_right()

	# Inset magnified plot
	# Inset image
	axins = zoomed_inset_axes(axMain,1.5,loc=1) # zoom=6
	# pcolormesh(x,y,d[i,:,:],cmap=colormap,vmin=dmin,vmax=dmax)
	# pcolormesh(x,y,z,cmap=colormap,vmin=-dmax,vmax=dmax)
	contourf(x,y,z,levels,cmap=get_cmap(colormap))
	hold('on')
	contour(x,y,mask,levels=[0],colors='grey',alpha=.5)
	# clabel(C,fontsize=14,fmt='%2.1f')#,manual=True)
	# subregion of the original image
	x1,x2,y1,y2 = 17,23,2.75,4.5
	axins.set_xlim(x1,x2)
	axins.set_ylim(y1,y2)
	xticks(visible=False)
	yticks(visible=False)
	setp(axins,xticks=[],yticks=[])
	# draw a bbox of the region of the inset axes in the parent axes and
	# connecting lines between the bbox and the inset axes area
	mark_inset(axMain, axins, loc1=2, loc2=4, fc="none", ec="0.5")
	draw()
	show()



	# title(Title)
	savefig(fname,bbox_inches='tight')
	# savefig('figures/' + var + zoom + '.png',bbox_inches='tight')
	show()
	# close()