Context Navigation

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

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

Line
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
18
19	class Quantity:
20
21
22	def __init__(self, domain, vertex_values=None):
23	#Initialise Quantity using optional vertex values.
24
25	from domain import Domain
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)
54	self.centroid_backup_values = zeros(N, Float)
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
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
81
82	def __len__(self):
83	"""
84	Returns number of intervals.
85	"""
86	return self.centroid_values.shape[0]
87
88	def interpolate(self):
89	"""
90	Compute interpolated values at centroid
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
111	of elements in the mesh. Otherwise it must be of dimension Nx2
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
124	if location not in ['vertices', 'centroids']:
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
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
155
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':
166
167	P = self.domain.centroids
168	self.set_values(f(P), location)
169	else:
170	#Vertices
171
172	P = self.domain.get_vertices()
173
174	for i in range(2):
175
176	self.vertex_values[:,i] = f(P[:,i])
177
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.
228	Permissible options are: vertices, centroid
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:
263	cells = self.domain.vertexlist[unique_vert_id]
264
265	#In case there are unused points
266	if cells is None:
267	msg = 'Unique vertex not associated with cells'
268	raise msg
269
270	# Go through all cells, vertex pairs
271	# Average the values
272	sum = 0
273	for cell_id, vertex_id in cells:
274	sum += self.vertex_values[cell_id, vertex_id]
275	vert_values.append(sum/len(cells))
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
282
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.
291	If requested the coordinates will be returned in 1D arrays X.
292
293	The connectivity is represented as an integer array, V, of dimension
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.
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
304	vertices 0, vertices 1
305
306
307	Calling convention
308	if x is True:
309	X,A,V = get_vertex_values
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:
381	X = self.domain.get_vertex_coordinates()
382
383	#X = C[:,0:6:2].copy()
384	#Y = C[:,1:6:2].copy()
385
386	return X.flat, A, V
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
430	#print 'compute_gradient'
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
450	k2 = k+2
451
452	q0 = Q[k0]
453	q1 = Q[k1]
454	q2 = Q[k2]
455
456	x0 = X[k0] #V0 centroid
457	x1 = X[k1] #V1 centroid
458	x2 = X[k2]
459
460	#Gradient
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)
465
466	elif k == N-1:
467
468	#Get data
469	k0 = k
470	k1 = k-1
471	k2 = k-2
472
473	q0 = Q[k0]
474	q1 = Q[k1]
475	q2 = Q[k2]
476
477	x0 = X[k0] #V0 centroid
478	x1 = X[k1] #V1 centroid
479	x2 = X[k2]
480
481	#Gradient
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)
486
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
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
521	#print 'compute_minmod_gradients'
522
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
549	k2 = k+2
550
551	q0 = Q[k0]
552	q1 = Q[k1]
553	q2 = Q[k2]
554
555	x0 = X[k0] #V0 centroid
556	x1 = X[k1] #V1 centroid
557	x2 = X[k2]
558
559	#Gradient
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)
564
565	elif k == N-1:
566
567	#Get data
568	k0 = k
569	k1 = k-1
570	k2 = k-2
571
572	q0 = Q[k0]
573	q1 = Q[k1]
574	q2 = Q[k2]
575
576	x0 = X[k0] #V0 centroid
577	x1 = X[k1] #V1 centroid
578	x2 = X[k2]
579
580	#Gradient
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)
585
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
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()
673	elif self.domain.limiter == "minmod_steve":
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
773	#print 'limit_range'
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])
806	#elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
807	# phi = 0.0
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]
897
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
903	self.centroid_backup_values[:] = (self.centroid_values).astype('f')
904
905	def saxpy_centroid_values(self,a,b):
906	# Call correct module function
907	# (either from this module or C-extension)
908	self.centroid_values[:] = (aself.centroid_values + bself.centroid_backup_values).astype('f')
909
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
923
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))
937
938	def newLinePlot(title='Simple Plot'):
939	import pylab as g
940	g.ion()
941	g.hold(False)
942	g.title(title)
943	g.xlabel('x')
944	g.ylabel('y')
945
946
947	def linePlot(x,y):
948	import pylab as g
949	g.plot(x.flat,y.flat)
950
951
952	def closePlots():
953	import pylab as g
954	g.close('all')
955
956	if __name__ == "__main__":
957	#from domain import Domain
958	from shallow_water_domain import Domain
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
966	Q1 = Quantity(D1, vertex_values)
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()
1006	#Q1.limit_minmod()
1007
1008	newLinePlot('plots')
1009	linePlot(Xc,Qc)
1010	raw_input('press return')
1011
1012	points2 = arange(10)
1013	D2 = Domain(points2)
1014
1015	Q2 = Quantity(D2)
1016	Q2.set_values(fun,'vertices')
1017	Xc = Q2.domain.vertices
1018	Qc = Q2.vertex_values
1019	linePlot(Xc,Qc)
1020	raw_input('press return')
1021
1022
1023	Q2.extrapolate_second_order()
1024	#Q2.limit_minmod()
1025	Xc = Q2.domain.vertices
1026	Qc = Q2.vertex_values
1027	print Q2.centroid_values
1028	print Qc
1029	linePlot(Xc,Qc)
1030	raw_input('press return')
1031
1032
1033	for i in range(10):
1034	import pylab as g
1035	g.hold(True)
1036	fun = FunClass(i/10.0)
1037	Q2.set_values(fun,'centroids')
1038	Q2.extrapolate_second_order()
1039	#Q2.limit_minmod()
1040	Qc = Q2.vertex_values
1041	linePlot(Xc,Qc)
1042	raw_input('press return')
1043
1044	raw_input('press return to quit')
1045	closePlots()

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

Download in other formats: