Context Navigation

source: anuga_work/development/anuga_1d/quantity.py @ 6806

Last change on this file since 6806 was 6453, checked in by steve, 16 years ago
Added Padarn's (2008/2009 summer student) files
File size: 30.7 KB

Rev	Line
[5535]	1	"""Class Quantity - Implements values at each 1d element
	2
	3	To create:
	4
	5	Quantity(domain, vertex_values)
	6
	7	domain: Associated domain structure. Required.
	8
	9	vertex_values: N x 2 array of values at each vertex for each element.
	10	Default None
	11
	12	If vertex_values are None Create array of zeros compatible with domain.
	13	Otherwise check that it is compatible with dimenions of domain.
	14	Otherwise raise an exception
	15	"""
	16
	17
[5858]	18
[5535]	19	class Quantity:
	20
[5858]	21
[5535]	22	def __init__(self, domain, vertex_values=None):
[5858]	23	#Initialise Quantity using optional vertex values.
	24
[5536]	25	from domain import Domain
[5535]	26	from Numeric import array, zeros, Float
	27
	28	msg = 'First argument in Quantity.__init__ '
	29	msg += 'must be of class Domain (or a subclass thereof)'
	30	assert isinstance(domain, Domain), msg
	31
	32	if vertex_values is None:
	33	N = domain.number_of_elements
	34	self.vertex_values = zeros((N, 2), Float)
	35	else:
	36	self.vertex_values = array(vertex_values, Float)
	37
	38	N, V = self.vertex_values.shape
	39	assert V == 2,\
	40	'Two vertex values per element must be specified'
	41
	42
	43	msg = 'Number of vertex values (%d) must be consistent with'\
	44	%N
	45	msg += 'number of elements in specified domain (%d).'\
	46	%domain.number_of_elements
	47
	48	assert N == domain.number_of_elements, msg
	49
	50	self.domain = domain
	51
	52	#Allocate space for other quantities
	53	self.centroid_values = zeros(N, Float)
[5738]	54	self.centroid_backup_values = zeros(N, Float)
[5535]	55	#self.edge_values = zeros((N, 2), Float)
	56	#edge values are values of the ends of each interval
	57
	58	#Intialise centroid values
	59	self.interpolate()
	60
[5536]	61
	62	from Numeric import zeros, Float
	63
	64	#Allocate space for boundary values
	65	#L = len(domain.boundary)
	66	self.boundary_values = zeros(2, Float) #assumes no parrellism
	67
	68	#Allocate space for updates of conserved quantities by
	69	#flux calculations and forcing functions
	70
	71	N = domain.number_of_elements
	72	self.explicit_update = zeros(N, Float )
	73	self.semi_implicit_update = zeros(N, Float )
	74
	75	self.gradients = zeros(N, Float)
	76	self.qmax = zeros(self.centroid_values.shape, Float)
	77	self.qmin = zeros(self.centroid_values.shape, Float)
	78
	79	self.beta = domain.beta
	80
[5858]	81
[5535]	82	def __len__(self):
[5858]	83	"""
	84	Returns number of intervals.
	85	"""
[5535]	86	return self.centroid_values.shape[0]
	87
	88	def interpolate(self):
[5858]	89	"""
	90	Compute interpolated values at centroid
[5535]	91	Pre-condition: vertex_values have been set
	92	"""
	93
	94	N = self.vertex_values.shape[0]
	95	for i in range(N):
	96	v0 = self.vertex_values[i, 0]
	97	v1 = self.vertex_values[i, 1]
	98
	99	self.centroid_values[i] = (v0 + v1)/2.0
	100
	101	def set_values(self, X, location='vertices'):
	102	"""Set values for quantity
	103
	104	X: Compatible list, Numeric array (see below), constant or function
	105	location: Where values are to be stored.
	106	Permissible options are: vertices, centroid
	107	Default is "vertices"
	108
	109	In case of location == 'centroid' the dimension values must
	110	be a list of a Numerical array of length N, N being the number
[5858]	111	of elements in the mesh. Otherwise it must be of dimension Nx2
[5535]	112
	113	The values will be stored in elements following their
	114	internal ordering.
	115
	116	If values are described a function, it will be evaluated at specified points
	117
	118	If selected location is vertices, values for centroid and edges
	119	will be assigned interpolated values.
	120	In any other case, only values for the specified locations
	121	will be assigned and the others will be left undefined.
	122	"""
	123
[5858]	124	if location not in ['vertices', 'centroids']:
[5535]	125	msg = 'Invalid location: %s, (possible choices vertices, centroids)' %location
	126	raise msg
	127
	128	if X is None:
	129	msg = 'Given values are None'
	130	raise msg
	131
	132	import types
	133
	134	if callable(X):
	135	#Use function specific method
	136	self.set_function_values(X, location)
[6453]	137
[5535]	138	elif type(X) in [types.FloatType, types.IntType, types.LongType]:
	139	if location == 'centroids':
	140	self.centroid_values[:] = X
	141	else:
	142	self.vertex_values[:] = X
	143
	144	else:
	145	#Use array specific method
	146	self.set_array_values(X, location)
	147
	148	if location == 'vertices':
	149	#Intialise centroid and edge values
	150	self.interpolate()
	151
	152
	153
	154
[5858]	155
[5535]	156	def set_function_values(self, f, location='vertices'):
	157	"""Set values for quantity using specified function
	158
	159	f: x -> z Function where x and z are arrays
	160	location: Where values are to be stored.
	161	Permissible options are: vertices, centroid
	162	Default is "vertices"
	163	"""
	164
	165	if location == 'centroids':
[6453]	166
[5535]	167	P = self.domain.centroids
	168	self.set_values(f(P), location)
	169	else:
	170	#Vertices
[6453]	171
[5535]	172	P = self.domain.get_vertices()
	173
	174	for i in range(2):
[6453]	175
	176	self.vertex_values[:,i] = f(P[:,i])
	177
[5535]	178	def set_array_values(self, values, location='vertices'):
	179	"""Set values for quantity
	180
	181	values: Numeric array
	182	location: Where values are to be stored.
	183	Permissible options are: vertices, centroid, edges
	184	Default is "vertices"
	185
	186	In case of location == 'centroid' the dimension values must
	187	be a list of a Numerical array of length N, N being the number
	188	of elements in the mesh. Otherwise it must be of dimension Nx2
	189
	190	The values will be stored in elements following their
	191	internal ordering.
	192
	193	If selected location is vertices, values for centroid
	194	will be assigned interpolated values.
	195	In any other case, only values for the specified locations
	196	will be assigned and the others will be left undefined.
	197	"""
	198
	199	from Numeric import array, Float
	200
	201	values = array(values).astype(Float)
	202
	203	N = self.centroid_values.shape[0]
	204
	205	msg = 'Number of values must match number of elements'
	206	assert values.shape[0] == N, msg
	207
	208	if location == 'centroids':
	209	assert len(values.shape) == 1, 'Values array must be 1d'
	210	self.centroid_values = values
	211	#elif location == 'edges':
	212	# assert len(values.shape) == 2, 'Values array must be 2d'
	213	# msg = 'Array must be N x 2'
	214	# self.edge_values = values
	215	else:
	216	assert len(values.shape) == 2, 'Values array must be 2d'
	217	msg = 'Array must be N x 2'
	218	assert values.shape[1] == 2, msg
	219
	220	self.vertex_values = values
	221
	222
	223	def get_values(self, location='vertices', indices = None):
	224	"""get values for quantity
	225
	226	return X, Compatible list, Numeric array (see below)
	227	location: Where values are to be stored.
[5858]	228	Permissible options are: vertices, centroid
[5535]	229	and unique vertices. Default is 'vertices'
	230
	231	In case of location == 'centroids' the dimension values must
	232	be a list of a Numerical array of length N, N being the number
	233	of elements. Otherwise it must be of dimension Nx3
	234
	235	The returned values with be a list the length of indices
	236	(N if indices = None). Each value will be a list of the three
	237	vertex values for this quantity.
	238
	239	Indices is the set of element ids that the operation applies to.
	240
	241	"""
	242	from Numeric import take
	243
	244	if location not in ['vertices', 'centroids', 'unique vertices']:
	245	msg = 'Invalid location: %s' %location
	246	raise msg
	247
	248	import types, Numeric
	249	assert type(indices) in [types.ListType, types.NoneType,
	250	Numeric.ArrayType],\
	251	'Indices must be a list or None'
	252
	253	if location == 'centroids':
	254	if (indices == None):
	255	indices = range(len(self))
	256	return take(self.centroid_values,indices)
	257	elif location == 'unique vertices':
	258	if (indices == None):
	259	indices=range(self.domain.coordinates.shape[0])
	260	vert_values = []
	261	#Go through list of unique vertices
	262	for unique_vert_id in indices:
[5858]	263	cells = self.domain.vertexlist[unique_vert_id]
[5535]	264
	265	#In case there are unused points
[5858]	266	if cells is None:
	267	msg = 'Unique vertex not associated with cells'
[5535]	268	raise msg
	269
[5858]	270	# Go through all cells, vertex pairs
[5535]	271	# Average the values
	272	sum = 0
[5858]	273	for cell_id, vertex_id in cells:
	274	sum += self.vertex_values[cell_id, vertex_id]
	275	vert_values.append(sum/len(cells))
[5535]	276	return Numeric.array(vert_values)
	277	else:
	278	if (indices == None):
	279	indices = range(len(self))
	280	return take(self.vertex_values,indices)
	281
[5858]	282
[5535]	283	def get_vertex_values(self,
	284	x=True,
	285	smooth = None,
	286	precision = None,
	287	reduction = None):
	288	"""Return vertex values like an OBJ format
	289
	290	The vertex values are returned as one sequence in the 1D float array A.
[5858]	291	If requested the coordinates will be returned in 1D arrays X.
[5535]	292
	293	The connectivity is represented as an integer array, V, of dimension
[5858]	294	M x 2, where M is the number of volumes. Each row has two indices
	295	into the X, A arrays defining the element.
[5535]	296
	297	if smooth is True, vertex values corresponding to one common
	298	coordinate set will be smoothed according to the given
	299	reduction operator. In this case vertex coordinates will be
	300	de-duplicated.
	301
	302	If no smoothings is required, vertex coordinates and values will
	303	be aggregated as a concatenation of values at
[5858]	304	vertices 0, vertices 1
[5535]	305
	306
	307	Calling convention
[5858]	308	if x is True:
	309	X,A,V = get_vertex_values
[5535]	310	else:
	311	A,V = get_vertex_values
	312
	313	"""
	314
	315	from Numeric import concatenate, zeros, Float, Int, array, reshape
	316
	317
	318	if smooth is None:
	319	smooth = self.domain.smooth
	320
	321	if precision is None:
	322	precision = Float
	323
	324	if reduction is None:
	325	reduction = self.domain.reduction
	326
	327	#Create connectivity
	328
	329	if smooth == True:
	330
	331	V = self.domain.get_vertices()
	332	N = len(self.domain.vertexlist)
	333	#N = len(self.domain.vertices)
	334	A = zeros(N, precision)
	335
	336	#Smoothing loop
	337	for k in range(N):
	338	L = self.domain.vertexlist[k]
	339	#L = self.domain.vertices[k]
	340
	341	#Go through all triangle, vertex pairs
	342	#contributing to vertex k and register vertex value
	343
	344	if L is None: continue #In case there are unused points
	345
	346	contributions = []
	347	for volume_id, vertex_id in L:
	348	v = self.vertex_values[volume_id, vertex_id]
	349	contributions.append(v)
	350
	351	A[k] = reduction(contributions)
	352
	353	if x is True:
	354	#X = self.domain.coordinates[:,0].astype(precision)
	355	X = self.domain.coordinates[:].astype(precision)
	356	#Y = self.domain.coordinates[:,1].astype(precision)
	357
	358	#return X, Y, A, V
	359	return X, A, V
	360
	361	#else:
	362	return A, V
	363	else:
	364	#Don't smooth
	365	#obj machinery moved to general_mesh
	366
	367	# Create a V like [[0 1 2], [3 4 5]....[3m-2 3m-1 3*m]]
	368	# These vert_id's will relate to the verts created below
	369	#m = len(self.domain) #Number of volumes
	370	#M = 3*m #Total number of unique vertices
	371	#V = reshape(array(range(M)).astype(Int), (m,3))
	372
	373	#V = self.domain.get_triangles(obj=True)
	374	V = self.domain.get_vertices
	375	#FIXME use get_vertices, when ready
	376
	377	A = self.vertex_values.flat
	378
	379	#Do vertex coordinates
	380	if x is True:
[5858]	381	X = self.domain.get_vertex_coordinates()
[5535]	382
[5858]	383	#X = C[:,0:6:2].copy()
	384	#Y = C[:,1:6:2].copy()
[5535]	385
[5858]	386	return X.flat, A, V
[5535]	387	else:
	388	return A, V
	389
	390	def get_integral(self):
	391	"""Compute the integral of quantity across entire domain
	392	"""
	393	integral = 0
	394	for k in range(self.domain.number_of_elements):
	395	area = self.domain.areas[k]
	396	qc = self.centroid_values[k]
	397	integral += qc*area
	398
	399	return integral
	400
	401
	402	def update(self, timestep):
	403	"""Update centroid values based on values stored in
	404	explicit_update and semi_implicit_update as well as given timestep
	405	"""
	406
	407	from Numeric import sum, equal, ones, Float
	408
	409	N = self.centroid_values.shape[0]
	410
	411	#Explicit updates
	412	self.centroid_values += timestep*self.explicit_update
	413
	414	#Semi implicit updates
	415	denominator = ones(N, Float)-timestep*self.semi_implicit_update
	416
	417	if sum(equal(denominator, 0.0)) > 0.0:
	418	msg = 'Zero division in semi implicit update. Call Stephen :-)'
	419	raise msg
	420	else:
	421	#Update conserved_quantities from semi implicit updates
	422	self.centroid_values /= denominator
	423
	424
	425	def compute_gradients(self):
	426	"""Compute gradients of piecewise linear function defined by centroids of
	427	neighbouring volumes.
	428	"""
	429
[5858]	430	#print 'compute_gradient'
[5535]	431
	432	from Numeric import array, zeros, Float
	433
	434	N = self.centroid_values.shape[0]
	435
	436
	437	G = self.gradients
	438	Q = self.centroid_values
	439	X = self.domain.centroids
	440
	441	for k in range(N):
	442
	443	# first and last elements have boundaries
	444
	445	if k == 0:
	446
	447	#Get data
	448	k0 = k
	449	k1 = k+1
[5858]	450	k2 = k+2
[5535]	451
	452	q0 = Q[k0]
	453	q1 = Q[k1]
[5858]	454	q2 = Q[k2]
[5535]	455
	456	x0 = X[k0] #V0 centroid
	457	x1 = X[k1] #V1 centroid
[5858]	458	x2 = X[k2]
[5535]	459
	460	#Gradient
[5858]	461	#G[k] = (q1 - q0)/(x1 - x0)
	462
	463	G[k] = (q1 - q0)(x2 - x0)(x2 - x0) - (q2 - q0)(x1 - x0)(x1 - x0)
	464	G[k] /= (x1 - x0)(x2 - x0)(x2 - x1)
[5535]	465
	466	elif k == N-1:
	467
	468	#Get data
	469	k0 = k
	470	k1 = k-1
[5858]	471	k2 = k-2
[5535]	472
	473	q0 = Q[k0]
	474	q1 = Q[k1]
[5858]	475	q2 = Q[k2]
[5535]	476
	477	x0 = X[k0] #V0 centroid
	478	x1 = X[k1] #V1 centroid
[5858]	479	x2 = X[k2]
[5535]	480
	481	#Gradient
[5858]	482	#G[k] = (q1 - q0)/(x1 - x0)
	483
	484	G[k] = (q1 - q0)(x2 - x0)(x2 - x0) - (q2 - q0)(x1 - x0)(x1 - x0)
	485	G[k] /= (x1 - x0)(x2 - x0)(x2 - x1)
[5535]	486
[5858]	487	## q0 = Q[k0]
	488	## q1 = Q[k1]
	489	##
	490	## x0 = X[k0] #V0 centroid
	491	## x1 = X[k1] #V1 centroid
	492	##
	493	## #Gradient
	494	## G[k] = (q1 - q0)/(x1 - x0)
	495
[5535]	496	else:
	497	#Interior Volume (2 neighbours)
	498
	499	#Get data
	500	k0 = k-1
	501	k2 = k+1
	502
	503	q0 = Q[k0]
	504	q1 = Q[k]
	505	q2 = Q[k2]
	506
	507	x0 = X[k0] #V0 centroid
	508	x1 = X[k] #V1 centroid (Self)
	509	x2 = X[k2] #V2 centroid
	510
	511	#Gradient
	512	#G[k] = (q2-q0)/(x2-x0)
	513	G[k] = ((q0-q1)/(x0-x1)(x2-x1) - (q2-q1)/(x2-x1)(x0-x1))/(x2-x0)
	514
	515
	516	def compute_minmod_gradients(self):
	517	"""Compute gradients of piecewise linear function defined by centroids of
	518	neighbouring volumes.
	519	"""
	520
[5858]	521	#print 'compute_minmod_gradients'
	522
[5535]	523	from Numeric import array, zeros, Float,sign
	524
	525	def xmin(a,b):
	526	return 0.5(sign(a)+sign(b))min(abs(a),abs(b))
	527
	528	def xmic(t,a,b):
	529	return xmin(txmin(a,b), 0.50(a+b) )
	530
	531
	532
	533	N = self.centroid_values.shape[0]
	534
	535
	536	G = self.gradients
	537	Q = self.centroid_values
	538	X = self.domain.centroids
	539
	540	for k in range(N):
	541
	542	# first and last elements have boundaries
	543
	544	if k == 0:
	545
	546	#Get data
	547	k0 = k
	548	k1 = k+1
[5858]	549	k2 = k+2
[5535]	550
	551	q0 = Q[k0]
	552	q1 = Q[k1]
[5858]	553	q2 = Q[k2]
[5535]	554
	555	x0 = X[k0] #V0 centroid
	556	x1 = X[k1] #V1 centroid
[5858]	557	x2 = X[k2]
[5535]	558
	559	#Gradient
[5858]	560	#G[k] = (q1 - q0)/(x1 - x0)
	561
	562	G[k] = (q1 - q0)(x2 - x0)(x2 - x0) - (q2 - q0)(x1 - x0)(x1 - x0)
	563	G[k] /= (x1 - x0)(x2 - x0)(x2 - x1)
[5535]	564
	565	elif k == N-1:
	566
	567	#Get data
	568	k0 = k
	569	k1 = k-1
[5858]	570	k2 = k-2
[5535]	571
	572	q0 = Q[k0]
	573	q1 = Q[k1]
[5858]	574	q2 = Q[k2]
[5535]	575
	576	x0 = X[k0] #V0 centroid
	577	x1 = X[k1] #V1 centroid
[5858]	578	x2 = X[k2]
[5535]	579
	580	#Gradient
[5858]	581	#G[k] = (q1 - q0)/(x1 - x0)
	582
	583	G[k] = (q1 - q0)(x2 - x0)(x2 - x0) - (q2 - q0)(x1 - x0)(x1 - x0)
	584	G[k] /= (x1 - x0)(x2 - x0)(x2 - x1)
[5535]	585
[5858]	586	## #Get data
	587	## k0 = k
	588	## k1 = k-1
	589	##
	590	## q0 = Q[k0]
	591	## q1 = Q[k1]
	592	##
	593	## x0 = X[k0] #V0 centroid
	594	## x1 = X[k1] #V1 centroid
	595	##
	596	## #Gradient
	597	## G[k] = (q1 - q0)/(x1 - x0)
	598
[5535]	599	elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
	600	G[k] = 0.0
	601
	602	else:
	603	#Interior Volume (2 neighbours)
	604
	605	#Get data
	606	k0 = k-1
	607	k2 = k+1
	608
	609	q0 = Q[k0]
	610	q1 = Q[k]
	611	q2 = Q[k2]
	612
	613	x0 = X[k0] #V0 centroid
	614	x1 = X[k] #V1 centroid (Self)
	615	x2 = X[k2] #V2 centroid
	616
	617	# assuming uniform grid
	618	d1 = (q1 - q0)/(x1-x0)
	619	d2 = (q2 - q1)/(x2-x1)
	620
	621	#Gradient
	622	#G[k] = (d1+d2)*0.5
	623	#G[k] = (d1(x2-x1) - d2(x0-x1))/(x2-x0)
	624	G[k] = xmic( self.domain.beta, d1, d2 )
	625
	626
	627	def extrapolate_first_order(self):
	628	"""Extrapolate conserved quantities from centroid to
	629	vertices for each volume using
	630	first order scheme.
	631	"""
	632
	633	qc = self.centroid_values
	634	qv = self.vertex_values
	635
	636	for i in range(2):
	637	qv[:,i] = qc
	638
	639
	640	def extrapolate_second_order(self):
	641	"""Extrapolate conserved quantities from centroid to
	642	vertices for each volume using
	643	second order scheme.
	644	"""
	645	if self.domain.limiter == "pyvolution":
	646	#Z = self.gradients
	647	#print "gradients 1",Z
	648	self.compute_gradients()
	649	#print "gradients 2",Z
	650
	651	#Z = self.gradients
	652	#print "gradients 1",Z
	653	#self.compute_minmod_gradients()
	654	#print "gradients 2", Z
	655
	656	G = self.gradients
	657	V = self.domain.vertices
	658	qc = self.centroid_values
	659	qv = self.vertex_values
	660
	661	#Check each triangle
	662	for k in range(self.domain.number_of_elements):
	663	#Centroid coordinates
	664	x = self.domain.centroids[k]
	665
	666	#vertex coordinates
	667	x0, x1 = V[k,:]
	668
	669	#Extrapolate
	670	qv[k,0] = qc[k] + G[k]*(x0-x)
	671	qv[k,1] = qc[k] + G[k]*(x1-x)
	672	self.limit_pyvolution()
[5844]	673	elif self.domain.limiter == "minmod_steve":
[5535]	674	self.limit_minmod()
	675	else:
	676	self.limit_range()
	677
	678
	679
	680	def limit_minmod(self):
	681	#Z = self.gradients
	682	#print "gradients 1",Z
	683	self.compute_minmod_gradients()
	684	#print "gradients 2", Z
	685
	686	G = self.gradients
	687	V = self.domain.vertices
	688	qc = self.centroid_values
	689	qv = self.vertex_values
	690
	691	#Check each triangle
	692	for k in range(self.domain.number_of_elements):
	693	#Centroid coordinates
	694	x = self.domain.centroids[k]
	695
	696	#vertex coordinates
	697	x0, x1 = V[k,:]
	698
	699	#Extrapolate
	700	qv[k,0] = qc[k] + G[k]*(x0-x)
	701	qv[k,1] = qc[k] + G[k]*(x1-x)
	702
	703
	704	def limit_pyvolution(self):
	705	"""
	706	Limit slopes for each volume to eliminate artificial variance
	707	introduced by e.g. second order extrapolator
	708
	709	This is an unsophisticated limiter as it does not take into
	710	account dependencies among quantities.
	711
	712	precondition:
	713	vertex values are estimated from gradient
	714	postcondition:
	715	vertex values are updated
	716	"""
	717	from Numeric import zeros, Float
	718
	719	N = self.domain.number_of_elements
	720	beta = self.domain.beta
	721	#beta = 0.8
	722
	723	qc = self.centroid_values
	724	qv = self.vertex_values
	725
	726	#Find min and max of this and neighbour's centroid values
	727	qmax = self.qmax
	728	qmin = self.qmin
	729
	730	for k in range(N):
	731	qmax[k] = qmin[k] = qc[k]
	732	for i in range(2):
	733	n = self.domain.neighbours[k,i]
	734	if n >= 0:
	735	qn = qc[n] #Neighbour's centroid value
	736
	737	qmin[k] = min(qmin[k], qn)
	738	qmax[k] = max(qmax[k], qn)
	739
	740
	741	#Diffences between centroids and maxima/minima
	742	dqmax = qmax - qc
	743	dqmin = qmin - qc
	744
	745	#Deltas between vertex and centroid values
	746	dq = zeros(qv.shape, Float)
	747	for i in range(2):
	748	dq[:,i] = qv[:,i] - qc
	749
	750	#Phi limiter
	751	for k in range(N):
	752
	753	#Find the gradient limiter (phi) across vertices
	754	phi = 1.0
	755	for i in range(2):
	756	r = 1.0
	757	if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
	758	if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
	759
	760	phi = min( min(r*beta, 1), phi )
	761
	762	#Then update using phi limiter
	763	for i in range(2):
	764	qv[k,i] = qc[k] + phi*dq[k,i]
	765
	766	def limit_range(self):
	767	import sys
	768	from Numeric import zeros, Float
	769	from util import minmod, minmod_kurganov, maxmod, vanleer, vanalbada
	770	limiter = self.domain.limiter
	771	#print limiter
	772
[5858]	773	#print 'limit_range'
[5535]	774	N = self.domain.number_of_elements
	775	qc = self.centroid_values
	776	qv = self.vertex_values
	777	C = self.domain.centroids
	778	X = self.domain.vertices
	779	beta_p = zeros(N,Float)
	780	beta_m = zeros(N,Float)
	781	beta_x = zeros(N,Float)
	782
	783	for k in range(N):
	784
	785	n0 = self.domain.neighbours[k,0]
	786	n1 = self.domain.neighbours[k,1]
	787
	788	if ( n0 >= 0) & (n1 >= 0):
	789	#SLOPE DERIVATIVE LIMIT
	790	beta_p[k] = (qc[k]-qc[k-1])/(C[k]-C[k-1])
	791	beta_m[k] = (qc[k+1]-qc[k])/(C[k+1]-C[k])
	792	beta_x[k] = (qc[k+1]-qc[k-1])/(C[k+1]-C[k-1])
	793
	794	dq = zeros(qv.shape, Float)
	795	for i in range(2):
	796	dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
	797
	798	#Phi limiter
	799	for k in range(N):
	800	n0 = self.domain.neighbours[k,0]
	801	n1 = self.domain.neighbours[k,1]
	802	if n0 < 0:
	803	phi = (qc[k+1] - qc[k])/(C[k+1] - C[k])
	804	elif n1 < 0:
	805	phi = (qc[k] - qc[k-1])/(C[k] - C[k-1])
[5827]	806	#elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
	807	# phi = 0.0
[5535]	808	else:
	809	if limiter == "minmod":
	810	phi = minmod(beta_p[k],beta_m[k])
	811
	812	elif limiter == "minmod_kurganov":#Change this
	813	# Also known as monotonized central difference limiter
	814	# if theta = 2.0
	815	theta = 2.0
	816	phi = minmod_kurganov(thetabeta_p[k],thetabeta_m[k],beta_x[k])
	817	elif limiter == "superbee":
	818	slope1 = minmod(beta_m[k],2.0*beta_p[k])
	819	slope2 = minmod(2.0*beta_m[k],beta_p[k])
	820	phi = maxmod(slope1,slope2)
	821
	822	elif limiter == "vanleer":
	823	phi = vanleer(beta_p[k],beta_m[k])
	824
	825	elif limiter == "vanalbada":
	826	phi = vanalbada(beta_m[k],beta_p[k])
	827
	828	for i in range(2):
	829	qv[k,i] = qc[k] + phi*dq[k,i]
	830
	831	def limit_steve_slope(self):
	832
	833	import sys
	834	from Numeric import zeros, Float
	835	from util import minmod, minmod_kurganov, maxmod, vanleer
	836
	837	N = self.domain.number_of_elements
	838	limiter = self.domain.limiter
	839	limiter_type = self.domain.limiter_type
	840
	841	qc = self.centroid_values
	842	qv = self.vertex_values
	843
	844	#Find min and max of this and neighbour's centroid values
	845	beta_p = zeros(N,Float)
	846	beta_m = zeros(N,Float)
	847	beta_x = zeros(N,Float)
	848	C = self.domain.centroids
	849	X = self.domain.vertices
	850
	851	for k in range(N):
	852
	853	n0 = self.domain.neighbours[k,0]
	854	n1 = self.domain.neighbours[k,1]
	855
	856	if (n0 >= 0) & (n1 >= 0):
	857	# Check denominator not zero
	858	if (qc[k+1]-qc[k]) == 0.0:
	859	beta_p[k] = float(sys.maxint)
	860	beta_m[k] = float(sys.maxint)
	861	else:
	862	#STEVE LIMIT
	863	beta_p[k] = (qc[k]-qc[k-1])/(qc[k+1]-qc[k])
	864	beta_m[k] = (qc[k+2]-qc[k+1])/(qc[k+1]-qc[k])
	865
	866	#Deltas between vertex and centroid values
	867	dq = zeros(qv.shape, Float)
	868	for i in range(2):
	869	dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
	870
	871	#Phi limiter
	872	for k in range(N):
	873
	874	phi = 0.0
	875	if limiter == "flux_minmod":
	876	#FLUX MINMOD
	877	phi = minmod_kurganov(1.0,beta_m[k],beta_p[k])
	878	elif limiter == "flux_superbee":
	879	#FLUX SUPERBEE
	880	phi = max(0.0,min(1.0,2.0beta_m[k]),min(2.0,beta_m[k]))+max(0.0,min(1.0,2.0beta_p[k]),min(2.0,beta_p[k]))-1.0
	881	elif limiter == "flux_muscl":
	882	#FLUX MUSCL
	883	phi = max(0.0,min(2.0,2.0beta_m[k],2.0beta_p[k],0.5*(beta_m[k]+beta_p[k])))
	884	elif limiter == "flux_vanleer":
	885	#FLUX VAN LEER
	886	phi = (beta_m[k]+abs(beta_m[k]))/(1.0+abs(beta_m[k]))+(beta_p[k]+abs(beta_p[k]))/(1.0+abs(beta_p[k]))-1.0
	887
	888	#Then update using phi limiter
	889	n = self.domain.neighbours[k,1]
	890	if n>=0:
	891	#qv[k,0] = qc[k] - 0.5phi(qc[k+1]-qc[k])
	892	#qv[k,1] = qc[k] + 0.5phi(qc[k+1]-qc[k])
	893	qv[k,0] = qc[k] + 0.5phi(qv[k,0]-qc[k])
	894	qv[k,1] = qc[k] + 0.5phi(qv[k,1]-qc[k])
	895	else:
	896	qv[k,i] = qc[k]
[5536]	897
[5738]	898	def backup_centroid_values(self):
	899	# Call correct module function
	900	# (either from this module or C-extension)
	901	#backup_centroid_values(self)
	902
[5742]	903	self.centroid_backup_values[:] = (self.centroid_values).astype('f')
	904
[5738]	905	def saxpy_centroid_values(self,a,b):
	906	# Call correct module function
	907	# (either from this module or C-extension)
[5743]	908	self.centroid_values[:] = (aself.centroid_values + bself.centroid_backup_values).astype('f')
[5742]	909
[5536]	910	class Conserved_quantity(Quantity):
	911	"""Class conserved quantity adds to Quantity:
	912
	913	storage and method for updating, and
	914	methods for extrapolation from centropid to vertices inluding
	915	gradients and limiters
	916	"""
	917
	918	def __init__(self, domain, vertex_values=None):
	919	Quantity.__init__(self, domain, vertex_values)
	920
	921	print "Use Quantity instead of Conserved_quantity"
	922
[5535]	923
[5858]	924	##
	925	##def newLinePlot(title='Simple Plot'):
	926	## import Gnuplot
	927	## g = Gnuplot.Gnuplot()
	928	## g.title(title)
	929	## g('set data style linespoints')
	930	## g.xlabel('x')
	931	## g.ylabel('y')
	932	## return g
	933	##
	934	##def linePlot(g,x,y):
	935	## import Gnuplot
	936	## g.plot(Gnuplot.PlotItems.Data(x.flat,y.flat))
[5738]	937
[5535]	938	def newLinePlot(title='Simple Plot'):
[5858]	939	import pylab as g
	940	g.ion()
	941	g.hold(False)
[5535]	942	g.title(title)
	943	g.xlabel('x')
	944	g.ylabel('y')
[5858]	945
[5535]	946
[5858]	947	def linePlot(x,y):
	948	import pylab as g
	949	g.plot(x.flat,y.flat)
[5535]	950
	951
[5858]	952	def closePlots():
	953	import pylab as g
	954	g.close('all')
	955
[5535]	956	if __name__ == "__main__":
	957	#from domain import Domain
[5536]	958	from shallow_water_domain import Domain
[5535]	959	from Numeric import arange
	960
	961	points1 = [0.0, 1.0, 2.0, 3.0]
	962	vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
	963
	964	D1 = Domain(points1)
	965
[5536]	966	Q1 = Quantity(D1, vertex_values)
[5535]	967
	968	print Q1.vertex_values
	969	print Q1.centroid_values
	970
	971	new_vertex_values = [[2.0,1.0],[3.0,4.0],[-2.0,4.0]]
	972
	973	Q1.set_values(new_vertex_values)
	974
	975	print Q1.vertex_values
	976	print Q1.centroid_values
	977
	978	new_centroid_values = [20,30,40]
	979	Q1.set_values(new_centroid_values,'centroids')
	980
	981	print Q1.vertex_values
	982	print Q1.centroid_values
	983
	984	class FunClass:
	985	def __init__(self,value):
	986	self.value = value
	987
	988	def __call__(self,x):
	989	return self.value(x*2)
	990
	991
	992	fun = FunClass(1.0)
	993	Q1.set_values(fun,'vertices')
	994
	995	print Q1.vertex_values
	996	print Q1.centroid_values
	997
	998	Xc = Q1.domain.vertices
	999	Qc = Q1.vertex_values
	1000	print Xc
	1001	print Qc
	1002
	1003	Qc[1,0] = 3
	1004
	1005	Q1.extrapolate_second_order()
[5858]	1006	#Q1.limit_minmod()
[5535]	1007
[5858]	1008	newLinePlot('plots')
	1009	linePlot(Xc,Qc)
	1010	raw_input('press return')
[5535]	1011
	1012	points2 = arange(10)
	1013	D2 = Domain(points2)
	1014
[5536]	1015	Q2 = Quantity(D2)
[5535]	1016	Q2.set_values(fun,'vertices')
	1017	Xc = Q2.domain.vertices
	1018	Qc = Q2.vertex_values
[5858]	1019	linePlot(Xc,Qc)
	1020	raw_input('press return')
[5535]	1021
	1022
	1023	Q2.extrapolate_second_order()
[5858]	1024	#Q2.limit_minmod()
[5535]	1025	Xc = Q2.domain.vertices
	1026	Qc = Q2.vertex_values
	1027	print Q2.centroid_values
	1028	print Qc
[5858]	1029	linePlot(Xc,Qc)
[5535]	1030	raw_input('press return')
	1031
	1032
	1033	for i in range(10):
[5858]	1034	import pylab as g
	1035	g.hold(True)
[5535]	1036	fun = FunClass(i/10.0)
[5858]	1037	Q2.set_values(fun,'centroids')
	1038	Q2.extrapolate_second_order()
	1039	#Q2.limit_minmod()
[5535]	1040	Qc = Q2.vertex_values
[5858]	1041	linePlot(Xc,Qc)
	1042	raw_input('press return')
[5535]	1043
[5858]	1044	raw_input('press return to quit')
[6453]	1045	closePlots()

Note: See TracBrowser for help on using the repository browser.

Download in other formats: