Context Navigation

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

Last change on this file since 5858 was 5858, checked in by steve, 15 years ago

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

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

Download in other formats: