Context Navigation

source: anuga_work/development/sudi/sw_1d/periodic_waves/johns/quantity.py @ 7837

Last change on this file since 7837 was 7837, checked in by mungkasi, 14 years ago
Again, adding some codes for 1d problems on debris avalanche and periodic waves.
File size: 30.6 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	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
154
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.
225	Permissible options are: vertices, centroid
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:
260	cells = self.domain.vertexlist[unique_vert_id]
261
262	#In case there are unused points
263	if cells is None:
264	msg = 'Unique vertex not associated with cells'
265	raise msg
266
267	# Go through all cells, vertex pairs
268	# Average the values
269	sum = 0
270	for cell_id, vertex_id in cells:
271	sum += self.vertex_values[cell_id, vertex_id]
272	vert_values.append(sum/len(cells))
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
279
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.
288	If requested the coordinates will be returned in 1D arrays X.
289
290	The connectivity is represented as an integer array, V, of dimension
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.
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
301	vertices 0, vertices 1
302
303
304	Calling convention
305	if x is True:
306	X,A,V = get_vertex_values
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:
378	X = self.domain.get_vertex_coordinates()
379
380	#X = C[:,0:6:2].copy()
381	#Y = C[:,1:6:2].copy()
382
383	return X.flat, A, V
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
427	#print 'compute_gradient'
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
447	k2 = k+2
448
449	q0 = Q[k0]
450	q1 = Q[k1]
451	q2 = Q[k2]
452
453	x0 = X[k0] #V0 centroid
454	x1 = X[k1] #V1 centroid
455	x2 = X[k2]
456
457	#Gradient
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)
462
463	elif k == N-1:
464
465	#Get data
466	k0 = k
467	k1 = k-1
468	k2 = k-2
469
470	q0 = Q[k0]
471	q1 = Q[k1]
472	q2 = Q[k2]
473
474	x0 = X[k0] #V0 centroid
475	x1 = X[k1] #V1 centroid
476	x2 = X[k2]
477
478	#Gradient
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)
483
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
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
518	#print 'compute_minmod_gradients'
519
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
546	k2 = k+2
547
548	q0 = Q[k0]
549	q1 = Q[k1]
550	q2 = Q[k2]
551
552	x0 = X[k0] #V0 centroid
553	x1 = X[k1] #V1 centroid
554	x2 = X[k2]
555
556	#Gradient
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)
561
562	elif k == N-1:
563
564	#Get data
565	k0 = k
566	k1 = k-1
567	k2 = k-2
568
569	q0 = Q[k0]
570	q1 = Q[k1]
571	q2 = Q[k2]
572
573	x0 = X[k0] #V0 centroid
574	x1 = X[k1] #V1 centroid
575	x2 = X[k2]
576
577	#Gradient
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)
582
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
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()
670	elif self.domain.limiter == "minmod_steve":
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
770	#print 'limit_range'
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	#Note: neighbours=[-1, {0, 1, 2, ..., N-1}, -2]
782	n0 = self.domain.neighbours[k,0]
783	n1 = self.domain.neighbours[k,1]
784	if ( n0 >= 0) & (n1 >= 0):
785	#SLOPE DERIVATIVE LIMIT
786	beta_p[k] = (qc[k]-qc[k-1])/(C[k]-C[k-1])
787	beta_m[k] = (qc[k+1]-qc[k])/(C[k+1]-C[k])
788	beta_x[k] = (qc[k+1]-qc[k-1])/(C[k+1]-C[k-1])
789
790	dq = zeros(qv.shape, Float)
791	for i in range(2):
792	dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
793
794	#Phi limiter
795	for k in range(N):
796	#Note: neighbours=[-1, {0, 1, 2, ..., N-1}, -2]
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])
803	else:
804	if limiter == "minmod":
805	phi = minmod(beta_p[k],beta_m[k])
806	elif limiter == "minmod_kurganov":
807	# Also known as monotonized central difference limiter if theta = 2.0
808	theta = 2.0
809	phi = minmod_kurganov(thetabeta_p[k],thetabeta_m[k],beta_x[k])
810	elif limiter == "superbee":
811	slope1 = minmod(beta_m[k],2.0*beta_p[k])
812	slope2 = minmod(2.0*beta_m[k],beta_p[k])
813	phi = maxmod(slope1,slope2)
814	elif limiter == "vanleer":
815	phi = vanleer(beta_p[k],beta_m[k])
816	elif limiter == "vanalbada":
817	phi = vanalbada(beta_m[k],beta_p[k])
818	for i in range(2):
819	qv[k,i] = qc[k] + phi*dq[k,i]
820
821	def limit_steve_slope(self):
822
823	import sys
824	from Numeric import zeros, Float
825	from util import minmod, minmod_kurganov, maxmod, vanleer
826
827	N = self.domain.number_of_elements
828	limiter = self.domain.limiter
829	limiter_type = self.domain.limiter_type
830
831	qc = self.centroid_values
832	qv = self.vertex_values
833
834	#Find min and max of this and neighbour's centroid values
835	beta_p = zeros(N,Float)
836	beta_m = zeros(N,Float)
837	beta_x = zeros(N,Float)
838	C = self.domain.centroids
839	X = self.domain.vertices
840
841	for k in range(N):
842
843	n0 = self.domain.neighbours[k,0]
844	n1 = self.domain.neighbours[k,1]
845
846	if (n0 >= 0) & (n1 >= 0):
847	# Check denominator not zero
848	if (qc[k+1]-qc[k]) == 0.0:
849	beta_p[k] = float(sys.maxint)
850	beta_m[k] = float(sys.maxint)
851	else:
852	#STEVE LIMIT
853	beta_p[k] = (qc[k]-qc[k-1])/(qc[k+1]-qc[k])
854	beta_m[k] = (qc[k+2]-qc[k+1])/(qc[k+1]-qc[k])
855
856	#Deltas between vertex and centroid values
857	dq = zeros(qv.shape, Float)
858	for i in range(2):
859	dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
860
861	#Phi limiter
862	for k in range(N):
863
864	phi = 0.0
865	if limiter == "flux_minmod":
866	#FLUX MINMOD
867	phi = minmod_kurganov(1.0,beta_m[k],beta_p[k])
868	elif limiter == "flux_superbee":
869	#FLUX SUPERBEE
870	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
871	elif limiter == "flux_muscl":
872	#FLUX MUSCL
873	phi = max(0.0,min(2.0,2.0beta_m[k],2.0beta_p[k],0.5*(beta_m[k]+beta_p[k])))
874	elif limiter == "flux_vanleer":
875	#FLUX VAN LEER
876	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
877
878	#Then update using phi limiter
879	n = self.domain.neighbours[k,1]
880	if n>=0:
881	#qv[k,0] = qc[k] - 0.5phi(qc[k+1]-qc[k])
882	#qv[k,1] = qc[k] + 0.5phi(qc[k+1]-qc[k])
883	qv[k,0] = qc[k] + 0.5phi(qv[k,0]-qc[k])
884	qv[k,1] = qc[k] + 0.5phi(qv[k,1]-qc[k])
885	else:
886	qv[k,i] = qc[k]
887
888	def backup_centroid_values(self):
889	# Call correct module function
890	# (either from this module or C-extension)
891	#backup_centroid_values(self)
892
893	self.centroid_backup_values[:] = (self.centroid_values).astype('f')
894
895	def saxpy_centroid_values(self,a,b):
896	# Call correct module function
897	# (either from this module or C-extension)
898	self.centroid_values[:] = (aself.centroid_values + bself.centroid_backup_values).astype('f')
899
900	class Conserved_quantity(Quantity):
901	"""Class conserved quantity adds to Quantity:
902
903	storage and method for updating, and
904	methods for extrapolation from centropid to vertices inluding
905	gradients and limiters
906	"""
907
908	def __init__(self, domain, vertex_values=None):
909	Quantity.__init__(self, domain, vertex_values)
910
911	print "Use Quantity instead of Conserved_quantity"
912
913
914	##
915	##def newLinePlot(title='Simple Plot'):
916	## import Gnuplot
917	## g = Gnuplot.Gnuplot()
918	## g.title(title)
919	## g('set data style linespoints')
920	## g.xlabel('x')
921	## g.ylabel('y')
922	## return g
923	##
924	##def linePlot(g,x,y):
925	## import Gnuplot
926	## g.plot(Gnuplot.PlotItems.Data(x.flat,y.flat))
927
928	def newLinePlot(title='Simple Plot'):
929	import pylab as g
930	g.ion()
931	g.hold(False)
932	g.title(title)
933	g.xlabel('x')
934	g.ylabel('y')
935
936
937	def linePlot(x,y):
938	import pylab as g
939	g.plot(x.flat,y.flat)
940
941
942	def closePlots():
943	import pylab as g
944	g.close('all')
945
946	if __name__ == "__main__":
947	#from domain import Domain
948	from domain_johns import Domain
949	from Numeric import arange
950
951	points1 = [0.0, 1.0, 2.0, 3.0]
952	vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
953
954	D1 = Domain(points1)
955
956	Q1 = Quantity(D1, vertex_values)
957
958	print Q1.vertex_values
959	print Q1.centroid_values
960
961	new_vertex_values = [[2.0,1.0],[3.0,4.0],[-2.0,4.0]]
962
963	Q1.set_values(new_vertex_values)
964
965	print Q1.vertex_values
966	print Q1.centroid_values
967
968	new_centroid_values = [20,30,40]
969	Q1.set_values(new_centroid_values,'centroids')
970
971	print Q1.vertex_values
972	print Q1.centroid_values
973
974	class FunClass:
975	def __init__(self,value):
976	self.value = value
977
978	def __call__(self,x):
979	return self.value(x*2)
980
981
982	fun = FunClass(1.0)
983	Q1.set_values(fun,'vertices')
984
985	print Q1.vertex_values
986	print Q1.centroid_values
987
988	Xc = Q1.domain.vertices
989	Qc = Q1.vertex_values
990	print Xc
991	print Qc
992
993	Qc[1,0] = 3
994
995	Q1.extrapolate_second_order()
996	#Q1.limit_minmod()
997
998	newLinePlot('plots')
999	linePlot(Xc,Qc)
1000	raw_input('press return')
1001
1002	points2 = arange(10)
1003	D2 = Domain(points2)
1004
1005	Q2 = Quantity(D2)
1006	Q2.set_values(fun,'vertices')
1007	Xc = Q2.domain.vertices
1008	Qc = Q2.vertex_values
1009	linePlot(Xc,Qc)
1010	raw_input('press return')
1011
1012
1013	Q2.extrapolate_second_order()
1014	#Q2.limit_minmod()
1015	Xc = Q2.domain.vertices
1016	Qc = Q2.vertex_values
1017	print Q2.centroid_values
1018	print Qc
1019	linePlot(Xc,Qc)
1020	raw_input('press return')
1021
1022
1023	for i in range(10):
1024	import pylab as g
1025	g.hold(True)
1026	fun = FunClass(i/10.0)
1027	Q2.set_values(fun,'centroids')
1028	Q2.extrapolate_second_order()
1029	#Q2.limit_minmod()
1030	Qc = Q2.vertex_values
1031	linePlot(Xc,Qc)
1032	raw_input('press return')
1033
1034	raw_input('press return to quit')
1035	closePlots()

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

Download in other formats: