Context Navigation

source: inundation/pyvolution/domain.py @ 2866

Last change on this file since 2866 was 2866, checked in by ole, 19 years ago
Implemented Domain(meshfile) variant widely, updated documentation and added more statistics
File size: 29.6 KB

Line
1	"""Class Domain - 2D triangular domains for finite-volume computations of
2	conservation laws.
3
4
5	Copyright 2004
6	Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
7	Geoscience Australia
8	"""
9
10	from pyvolution.mesh import Mesh
11	from pyvolution.generic_boundary_conditions import Boundary
12	from pyvolution.generic_boundary_conditions import File_boundary
13	from pyvolution.generic_boundary_conditions import Dirichlet_boundary
14	from pyvolution.generic_boundary_conditions import Time_boundary
15	from pyvolution.generic_boundary_conditions import Transmissive_boundary
16	from pyvolution.pmesh2domain import pmesh_to_domain
17
18	import types
19
20	class Domain(Mesh):
21
22
23	def __init__(self,
24	source=None,
25	triangles=None,
26	boundary=None,
27	conserved_quantities=None,
28	other_quantities=None,
29	tagged_elements=None,
30	geo_reference=None,
31	use_inscribed_circle=False,
32	mesh_filename=None,
33	use_cache=False,
34	verbose=False,
35	full_send_dict=None,
36	ghost_recv_dict=None,
37	processor=0,
38	numproc=1):
39
40
41	"""Instantiate generic computational Domain.
42
43	Input:
44	source: Either a mesh filename or coordinates of mesh vertices.
45	If it is a filename values specified for triangles will
46	be overridden.
47	triangles: Mesh connectivity (see mesh.py for more information)
48	boundary: See mesh.py for more information
49
50	conserved_quantities: List of quantity names entering the
51	conservation equations
52	other_quantities: List of other quantity names
53
54	tagged_elements:
55	...
56
57
58	"""
59
60	# Determine whether source is a mesh filename or coordinates
61	if type(source) == types.StringType:
62	mesh_filename = source
63	else:
64	coordinates = source
65
66
67	# In case a filename has been specified, extract content
68	if mesh_filename is not None:
69	coordinates, triangles, boundary, vertex_quantity_dict, \
70	tagged_elements, geo_reference = \
71	pmesh_to_domain(file_name=mesh_filename,
72	use_cache=use_cache,
73	verbose=verbose)
74
75
76	# Initialise underlying mesh structure
77	Mesh.__init__(self, coordinates, triangles, boundary,
78	tagged_elements, geo_reference, use_inscribed_circle,
79	verbose=verbose)
80
81	if verbose: print 'Initialising Domain'
82	from Numeric import zeros, Float, Int
83	from quantity import Quantity, Conserved_quantity
84
85	# List of quantity names entering
86	# the conservation equations
87	if conserved_quantities is None:
88	self.conserved_quantities = []
89	else:
90	self.conserved_quantities = conserved_quantities
91
92	# List of other quantity names
93	if other_quantities is None:
94	self.other_quantities = []
95	else:
96	self.other_quantities = other_quantities
97
98
99	#Build dictionary of Quantity instances keyed by quantity names
100	self.quantities = {}
101
102	#FIXME: remove later - maybe OK, though....
103	for name in self.conserved_quantities:
104	self.quantities[name] = Conserved_quantity(self)
105	for name in self.other_quantities:
106	self.quantities[name] = Quantity(self)
107
108	#Create an empty list for explicit forcing terms
109	self.forcing_terms = []
110
111	#Setup the ghost cell communication
112	if full_send_dict is None:
113	self.full_send_dict = {}
114	else:
115	self.full_send_dict = full_send_dict
116
117	# List of other quantity names
118	if ghost_recv_dict is None:
119	self.ghost_recv_dict = []
120	else:
121	self.ghost_recv_dict = ghost_recv_dict
122
123	self.processor = processor
124	self.numproc = numproc
125
126	# Setup Communication Buffers
127
128	if verbose: print 'Domain: Set up communication buffers (parallel)'
129	self.nsys = len(self.conserved_quantities)
130	for key in self.full_send_dict:
131	buffer_shape = self.full_send_dict[key][0].shape[0]
132	self.full_send_dict[key].append(zeros( (buffer_shape,self.nsys) ,Float))
133
134
135	for key in self.ghost_recv_dict:
136	buffer_shape = self.ghost_recv_dict[key][0].shape[0]
137	self.ghost_recv_dict[key].append(zeros( (buffer_shape,self.nsys) ,Float))
138
139	#Defaults
140	from config import max_smallsteps, beta_w, beta_h, epsilon, CFL
141	self.beta_w = beta_w
142	self.beta_h = beta_h
143	self.epsilon = epsilon
144
145	#FIXME: Maybe have separate orders for h-limiter and w-limiter?
146	#Or maybe get rid of order altogether and use beta_w and beta_h
147	self.set_default_order(1)
148	#self.default_order = 1
149	#self.order = self.default_order
150
151	self.smallsteps = 0
152	self.max_smallsteps = max_smallsteps
153	self.number_of_steps = 0
154	self.number_of_first_order_steps = 0
155	self.CFL = CFL
156
157	#Model time
158	self.time = 0.0
159	self.finaltime = None
160	self.min_timestep = self.max_timestep = 0.0
161	self.starttime = 0 #Physical starttime if any (0 is 1 Jan 1970 00:00:00)
162
163	######OBSOLETE
164	#Origin in UTM coordinates
165	#FIXME: This should be set if read by a msh file
166	#self.zone = zone
167	#self.xllcorner = xllcorner
168	#self.yllcorner = yllcorner
169
170
171	#Checkpointing and storage
172	from config import default_datadir
173	self.datadir = default_datadir
174	self.filename = 'domain'
175	self.checkpoint = False
176
177	#MH310505 To avoid calculating the flux across each edge twice, keep an integer (boolean) array,
178	#to be used during the flux calculation
179	N=self.number_of_elements
180	self.already_computed_flux = zeros((N, 3), Int)
181
182	if mesh_filename is not None:
183	# If the mesh file passed any quantity values
184	# , initialise with these values.
185	if verbose: print 'Domain: Initialising quantity values'
186	self.set_quantity_vertices_dict(vertex_quantity_dict)
187
188
189	if verbose: print 'Domain: Done'
190
191
192
193
194	def set_default_order(self, n):
195	"""Set default (spatial) order to either 1 or 2
196	"""
197
198	msg = 'Default order must be either 1 or 2. I got %s' %n
199	assert n in [1,2], msg
200
201	self.default_order = n
202	self.order = self.default_order
203
204
205	#Public interface to Domain
206	def get_conserved_quantities(self, vol_id, vertex=None, edge=None):
207	"""Get conserved quantities at volume vol_id
208
209	If vertex is specified use it as index for vertex values
210	If edge is specified use it as index for edge values
211	If neither are specified use centroid values
212	If both are specified an exeception is raised
213
214	Return value: Vector of length == number_of_conserved quantities
215
216	"""
217
218	from Numeric import zeros, Float
219
220	if not (vertex is None or edge is None):
221	msg = 'Values for both vertex and edge was specified.'
222	msg += 'Only one (or none) is allowed.'
223	raise msg
224
225	q = zeros( len(self.conserved_quantities), Float)
226
227	for i, name in enumerate(self.conserved_quantities):
228	Q = self.quantities[name]
229	if vertex is not None:
230	q[i] = Q.vertex_values[vol_id, vertex]
231	elif edge is not None:
232	q[i] = Q.edge_values[vol_id, edge]
233	else:
234	q[i] = Q.centroid_values[vol_id]
235
236	return q
237
238	def set_time(self, time=0.0):
239	"""Set the model time (seconds)"""
240
241	self.time = time
242
243	def set_quantity_vertices_dict(self, quantity_dict):
244	"""Set values for named quantities.
245	The index is the quantity
246
247	name: Name of quantity
248	X: Compatible list, Numeric array, const or function (see below)
249
250	The values will be stored in elements following their
251	internal ordering.
252
253	"""
254	for key in quantity_dict.keys():
255	self.set_quantity(key, quantity_dict[key], location='vertices')
256
257
258	def set_quantity(self, name, args, *kwargs):
259	"""Set values for named quantity
260
261
262	One keyword argument is documented here:
263	expression = None, # Arbitrary expression
264
265	expression:
266	Arbitrary expression involving quantity names
267
268	See Quantity.set_values for further documentation.
269	"""
270
271	#FIXME (Ole): Allow new quantities here
272	#from quantity import Quantity, Conserved_quantity
273	#Create appropriate quantity object
274	##if name in self.conserved_quantities:
275	## self.quantities[name] = Conserved_quantity(self)
276	##else:
277	## self.quantities[name] = Quantity(self)
278
279
280	#Do the expression stuff
281	if kwargs.has_key('expression'):
282	expression = kwargs['expression']
283	del kwargs['expression']
284
285	Q = self.create_quantity_from_expression(expression)
286	kwargs['quantity'] = Q
287
288	#Assign values
289	self.quantities[name].set_values(args, *kwargs)
290
291
292	def get_quantity(self, name, location='vertices', indices = None):
293	"""Get quantity object.
294
295	name: Name of quantity
296
297	See methods inside the quantity object for more options
298	"""
299
300	return self.quantities[name] #.get_values( location, indices = indices)
301
302
303	def get_quantity_object(self, name):
304	"""Get object for named quantity
305
306	name: Name of quantity
307
308	FIXME: Obsolete
309	"""
310
311	print 'get_quantity_object has been deprecated. Please use get_quantity'
312	return self.quantities[name]
313
314
315	def create_quantity_from_expression(self, expression):
316	"""Create new quantity from other quantities using arbitrary expression
317
318	Combine existing quantities in domain using expression and return
319	result as a new quantity.
320
321	Note, the new quantity could e.g. be used in set_quantity
322
323	Valid expressions are limited to operators defined in class Quantity
324
325	Example:
326
327
328	"""
329
330	from util import apply_expression_to_dictionary
331	return apply_expression_to_dictionary(expression, self.quantities)
332
333
334
335
336	def set_boundary(self, boundary_map):
337	"""Associate boundary objects with tagged boundary segments.
338
339	Input boundary_map is a dictionary of boundary objects keyed
340	by symbolic tags to matched against tags in the internal dictionary
341	self.boundary.
342
343	As result one pointer to a boundary object is stored for each vertex
344	in the list self.boundary_objects.
345	More entries may point to the same boundary object
346
347	Schematically the mapping is from two dictionaries to one list
348	where the index is used as pointer to the boundary_values arrays
349	within each quantity.
350
351	self.boundary: (vol_id, edge_id): tag
352	boundary_map (input): tag: boundary_object
353	----------------------------------------------
354	self.boundary_objects: ((vol_id, edge_id), boundary_object)
355
356
357	Pre-condition:
358	self.boundary has been built.
359
360	Post-condition:
361	self.boundary_objects is built
362
363	If a tag from the domain doesn't appear in the input dictionary an
364	exception is raised.
365	However, if a tag is not used to the domain, no error is thrown.
366	FIXME: This would lead to implementation of a
367	default boundary condition
368
369	Note: If a segment is listed in the boundary dictionary and if it is
370	not None, it will become a boundary -
371	even if there is a neighbouring triangle.
372	This would be the case for internal boundaries
373
374	Boundary objects that are None will be skipped.
375
376	FIXME: If set_boundary is called multiple times and if Boundary
377	object is changed into None, the neighbour structure will not be
378	restored!!!
379
380
381	"""
382
383	self.boundary_objects = []
384	self.boundary_map = boundary_map #Store for use with eg. boundary_stats.
385
386	#FIXME: Try to remove the sorting and fix test_mesh.py
387	x = self.boundary.keys()
388	x.sort()
389
390	#Loop through edges that lie on the boundary and associate them with
391	#callable boundary objects depending on their tags
392	for k, (vol_id, edge_id) in enumerate(x):
393	tag = self.boundary[ (vol_id, edge_id) ]
394
395	if boundary_map.has_key(tag):
396	B = boundary_map[tag] #Get callable boundary object
397
398	if B is not None:
399	self.boundary_objects.append( ((vol_id, edge_id), B) )
400	self.neighbours[vol_id, edge_id] = -len(self.boundary_objects)
401	else:
402	pass
403	#FIXME: Check and perhaps fix neighbour structure
404
405
406	else:
407	msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag
408	msg += 'bound to a boundary object.\n'
409	msg += 'All boundary tags defined in domain must appear '
410	msg += 'in the supplied dictionary.\n'
411	msg += 'The tags are: %s' %self.get_boundary_tags()
412	raise msg
413
414
415	def set_region(self, functions):
416	# The order of functions in the list is used.
417	if type(functions) not in [types.ListType,types.TupleType]:
418	functions = [functions]
419	for function in functions:
420	for tag in self.tagged_elements.keys():
421	function(tag, self.tagged_elements[tag], self)
422
423	#Do we need to do this sort of thing?
424	#self = function(tag, self.tagged_elements[tag], self)
425
426	#MISC
427	def check_integrity(self):
428	Mesh.check_integrity(self)
429
430	for quantity in self.conserved_quantities:
431	msg = 'Conserved quantities must be a subset of all quantities'
432	assert quantity in self.quantities, msg
433
434	##assert hasattr(self, 'boundary_objects')
435
436	def write_time(self):
437	print self.timestepping_statistics()
438
439	#Old version
440	#if self.min_timestep == self.max_timestep:
441	# print 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
442	# %(self.time, self.min_timestep, self.number_of_steps,
443	# self.number_of_first_order_steps)
444	#elif self.min_timestep > self.max_timestep:
445	# print 'Time = %.4f, steps=%d (%d)'\
446	# %(self.time, self.number_of_steps,
447	# self.number_of_first_order_steps)
448	#else:
449	# print 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
450	# %(self.time, self.min_timestep,
451	# self.max_timestep, self.number_of_steps,
452	# self.number_of_first_order_steps)
453
454	def timestepping_statistics(self):
455	"""Return string with time stepping statistics for printing or logging
456	"""
457
458	msg = ''
459	if self.min_timestep == self.max_timestep:
460	msg += 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
461	%(self.time, self.min_timestep, self.number_of_steps,
462	self.number_of_first_order_steps)
463	elif self.min_timestep > self.max_timestep:
464	msg += 'Time = %.4f, steps=%d (%d)'\
465	%(self.time, self.number_of_steps,
466	self.number_of_first_order_steps)
467	else:
468	msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
469	%(self.time, self.min_timestep,
470	self.max_timestep, self.number_of_steps,
471	self.number_of_first_order_steps)
472
473	return msg
474
475
476	def write_boundary_statistics(self, quantities = None, tags = None):
477	print self.boundary_statistics(quantities, tags)
478
479	def boundary_statistics(self, quantities = None, tags = None):
480	"""Output statistics about boundary forcing at each timestep
481
482
483	Input:
484	quantities: either None, a string or a list of strings naming the quantities to be reported
485	tags: either None, a string or a list of strings naming the tags to be reported
486
487
488	Example output:
489	Tag 'wall':
490	stage in [2, 5.5]
491	xmomentum in []
492	ymomentum in []
493	Tag 'ocean'
494
495
496	If quantities are specified only report on those. Otherwise take all conserved quantities.
497	If tags are specified only report on those, otherwise take all tags.
498
499	"""
500
501	#Input checks
502	import types, string
503
504	if quantities is None:
505	quantities = self.conserved_quantities
506	elif type(quantities) == types.StringType:
507	quantities = [quantities] #Turn it into a list
508
509	msg = 'Keyword argument quantities must be either None, '
510	msg += 'string or list. I got %s' %str(quantities)
511	assert type(quantities) == types.ListType, msg
512
513
514	if tags is None:
515	tags = self.get_boundary_tags()
516	elif type(tags) == types.StringType:
517	tags = [tags] #Turn it into a list
518
519	msg = 'Keyword argument tags must be either None, '
520	msg += 'string or list. I got %s' %str(tags)
521	assert type(tags) == types.ListType, msg
522
523	#Determine width of longest quantity name (for cosmetic purposes)
524	maxwidth = 0
525	for name in quantities:
526	w = len(name)
527	if w > maxwidth:
528	maxwidth = w
529
530	#Output stats
531	msg = 'Boundary values at time %.4f:\n' %self.time
532	for tag in tags:
533	msg += ' %s:\n' %tag
534
535	for name in quantities:
536	q = self.quantities[name]
537
538	#Find range of boundary values for tag and q
539	maxval = minval = None
540	for i, ((vol_id, edge_id), B) in\
541	enumerate(self.boundary_objects):
542	if self.boundary[(vol_id, edge_id)] == tag:
543	v = q.boundary_values[i]
544	if minval is None or v < minval: minval = v
545	if maxval is None or v > maxval: maxval = v
546
547	if minval is None or maxval is None:
548	msg += ' Sorry no information available about' +\
549	' tag %s and quantity %s\n' %(tag, name)
550	else:
551	msg += ' %s in [%12.8f, %12.8f]\n'\
552	%(string.ljust(name, maxwidth), minval, maxval)
553
554
555	return msg
556
557
558	def get_name(self):
559	return self.filename
560
561	def set_name(self, name):
562	self.filename = name
563
564	def get_datadir(self):
565	return self.datadir
566
567	def set_datadir(self, name):
568	self.datadir = name
569
570
571
572	#def set_defaults(self):
573	# """Set default values for uninitialised quantities.
574	# Should be overridden or specialised by specific modules
575	# """#
576	#
577	# for name in self.conserved_quantities + self.other_quantities:
578	# self.set_quantity(name, 0.0)
579
580
581	###########################
582	#Main components of evolve
583
584	def evolve(self,
585	yieldstep = None,
586	finaltime = None,
587	duration = None,
588	skip_initial_step = False):
589	"""Evolve model through time starting from self.starttime.
590
591
592	yieldstep: Interval between yields where results are stored,
593	statistics written and domain inspected or
594	possibly modified. If omitted the internal predefined
595	max timestep is used.
596	Internally, smaller timesteps may be taken.
597
598	duration: Duration of simulation
599
600	finaltime: Time where simulation should end
601
602	If both duration and finaltime are given an exception is thrown.
603
604
605	skip_initial_step: Boolean flag that decides whether the first
606	yield step is skipped or not. This is useful for example to avoid
607	duplicate steps when multiple evolve processes are dove tailed.
608
609
610	Evolve is implemented as a generator and is to be called as such, e.g.
611
612	for t in domain.evolve(yieldstep, finaltime):
613	<Do something with domain and t>
614
615
616	All times are given in seconds
617
618	"""
619
620	from config import min_timestep, max_timestep, epsilon
621
622	#FIXME: Maybe lump into a larger check prior to evolving
623	msg = 'Boundary tags must be bound to boundary objects before evolving system, '
624	msg += 'e.g. using the method set_boundary.\n'
625	msg += 'This system has the boundary tags %s '\
626	%self.get_boundary_tags()
627	assert hasattr(self, 'boundary_objects'), msg
628
629	##self.set_defaults()
630
631	if yieldstep is None:
632	yieldstep = max_timestep
633	else:
634	yieldstep = float(yieldstep)
635
636	self.order = self.default_order
637
638
639	if finaltime is not None and duration is not None:
640	print 'F', finaltime, duration
641	msg = 'Only one of finaltime and duration may be specified'
642	raise msg
643	else:
644	if finaltime is not None:
645	self.finaltime = float(finaltime)
646	if duration is not None:
647	self.finaltime = self.starttime + float(duration)
648
649
650
651
652	self.yieldtime = 0.0 #Time between 'yields'
653
654	#Initialise interval of timestep sizes (for reporting only)
655	self.min_timestep = max_timestep
656	self.max_timestep = min_timestep
657	self.number_of_steps = 0
658	self.number_of_first_order_steps = 0
659
660	#update ghosts
661	self.update_ghosts()
662
663	#Initial update of vertex and edge values
664	self.distribute_to_vertices_and_edges()
665
666	#Initial update boundary values
667	self.update_boundary()
668
669	#Or maybe restore from latest checkpoint
670	if self.checkpoint is True:
671	self.goto_latest_checkpoint()
672
673	if skip_initial_step is False:
674	yield(self.time) #Yield initial values
675
676	while True:
677
678	#Compute fluxes across each element edge
679	self.compute_fluxes()
680
681	#Update timestep to fit yieldstep and finaltime
682	self.update_timestep(yieldstep, finaltime)
683
684	#Update conserved quantities
685	self.update_conserved_quantities()
686
687	#update ghosts
688	self.update_ghosts()
689
690	#Update vertex and edge values
691	self.distribute_to_vertices_and_edges()
692
693	#Update boundary values
694	self.update_boundary()
695
696	#Update time
697	self.time += self.timestep
698	self.yieldtime += self.timestep
699	self.number_of_steps += 1
700	if self.order == 1:
701	self.number_of_first_order_steps += 1
702
703	#Yield results
704	if finaltime is not None and self.time >= finaltime:
705
706	if self.time > finaltime:
707	#FIXME (Ole, 30 April 2006): Do we need this check?
708	print 'WARNING (domain.py): time overshot finaltime. Contact Ole.Nielsen@ga.gov.au'
709	self.time = finaltime
710
711	# Yield final time and stop
712	yield(self.time)
713	break
714
715
716	if self.yieldtime >= yieldstep:
717	# Yield (intermediate) time and allow inspection of domain
718
719	if self.checkpoint is True:
720	self.store_checkpoint()
721	self.delete_old_checkpoints()
722
723	#Pass control on to outer loop for more specific actions
724	yield(self.time)
725
726	# Reinitialise
727	self.yieldtime = 0.0
728	self.min_timestep = max_timestep
729	self.max_timestep = min_timestep
730	self.number_of_steps = 0
731	self.number_of_first_order_steps = 0
732
733
734	def evolve_to_end(self, finaltime = 1.0):
735	"""Iterate evolve all the way to the end
736	"""
737
738	for _ in self.evolve(yieldstep=None, finaltime=finaltime):
739	pass
740
741
742
743	def update_boundary(self):
744	"""Go through list of boundary objects and update boundary values
745	for all conserved quantities on boundary.
746	"""
747
748	#FIXME: Update only those that change (if that can be worked out)
749	#FIXME: Boundary objects should not include ghost nodes.
750	for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
751	if B is None:
752	print 'WARNING: Ignored boundary segment %d (None)'
753	else:
754	q = B.evaluate(vol_id, edge_id)
755
756	for j, name in enumerate(self.conserved_quantities):
757	Q = self.quantities[name]
758	Q.boundary_values[i] = q[j]
759
760
761	def compute_fluxes(self):
762	msg = 'Method compute_fluxes must be overridden by Domain subclass'
763	raise msg
764
765
766	def update_timestep(self, yieldstep, finaltime):
767
768	from config import min_timestep, max_timestep
769
770	# self.timestep is calculated from speed of characteristics
771	# Apply CFL condition here
772	timestep = min(self.CFL*self.timestep,max_timestep)
773
774	#Record maximal and minimal values of timestep for reporting
775	self.max_timestep = max(timestep, self.max_timestep)
776	self.min_timestep = min(timestep, self.min_timestep)
777
778	#Protect against degenerate time steps
779	if timestep < min_timestep:
780
781	#Number of consecutive small steps taken b4 taking action
782	self.smallsteps += 1
783
784	if self.smallsteps > self.max_smallsteps:
785	self.smallsteps = 0 #Reset
786
787	if self.order == 1:
788	msg = 'WARNING: Too small timestep %.16f reached '\
789	%timestep
790	msg += 'even after %d steps of 1 order scheme'\
791	%self.max_smallsteps
792	print msg
793	timestep = min_timestep #Try enforcing min_step
794
795	#raise msg
796	else:
797	#Try to overcome situation by switching to 1 order
798	self.order = 1
799
800	else:
801	self.smallsteps = 0
802	if self.order == 1 and self.default_order == 2:
803	self.order = 2
804
805
806	#Ensure that final time is not exceeded
807	if finaltime is not None and self.time + timestep > finaltime:
808	timestep = finaltime-self.time
809
810	#Ensure that model time is aligned with yieldsteps
811	if self.yieldtime + timestep > yieldstep:
812	timestep = yieldstep-self.yieldtime
813
814	self.timestep = timestep
815
816
817
818	def compute_forcing_terms(self):
819	"""If there are any forcing functions driving the system
820	they should be defined in Domain subclass and appended to
821	the list self.forcing_terms
822	"""
823
824	for f in self.forcing_terms:
825	f(self)
826
827
828
829	def update_conserved_quantities(self):
830	"""Update vectors of conserved quantities using previously
831	computed fluxes specified forcing functions.
832	"""
833
834	from Numeric import ones, sum, equal, Float
835
836	N = self.number_of_elements
837	d = len(self.conserved_quantities)
838
839	timestep = self.timestep
840
841	#Compute forcing terms
842	self.compute_forcing_terms()
843
844	#Update conserved_quantities
845	for name in self.conserved_quantities:
846	Q = self.quantities[name]
847	Q.update(timestep)
848
849	#Clean up
850	#Note that Q.explicit_update is reset by compute_fluxes
851
852	#MH090605 commented out the following since semi_implicit_update is now re-initialized
853	#at the end of the _update function in quantity_ext.c (This is called by the
854	#preceeding Q.update(timestep) statement above).
855	#For run_profile.py with N=128, the time of update_conserved_quantities is cut from 14.00 secs
856	#to 8.35 secs
857
858	#Q.semi_implicit_update[:] = 0.0
859
860	def update_ghosts(self):
861	pass
862
863	def distribute_to_vertices_and_edges(self):
864	"""Extrapolate conserved quantities from centroid to
865	vertices and edge-midpoints for each volume
866
867	Default implementation is straight first order,
868	i.e. constant values throughout each element and
869	no reference to non-conserved quantities.
870	"""
871
872	for name in self.conserved_quantities:
873	Q = self.quantities[name]
874	if self.order == 1:
875	Q.extrapolate_first_order()
876	elif self.order == 2:
877	Q.extrapolate_second_order()
878	Q.limit()
879	else:
880	raise 'Unknown order'
881	Q.interpolate_from_vertices_to_edges()
882
883
884	def centroid_norm(self, quantity, normfunc):
885	"""Calculate the norm of the centroid values
886	of a specific quantity, using normfunc.
887
888	normfunc should take a list to a float.
889
890	common normfuncs are provided in the module utilities.norms
891	"""
892	return normfunc(self.quantities[quantity].centroid_values)
893
894
895
896	##############################################
897	#Initialise module
898
899	#Optimisation with psyco
900	from config import use_psyco
901	if use_psyco:
902	try:
903	import psyco
904	except:
905	import os
906	if os.name == 'posix' and os.uname()[4] == 'x86_64':
907	pass
908	#Psyco isn't supported on 64 bit systems, but it doesn't matter
909	else:
910	msg = 'WARNING: psyco (speedup) could not import'+\
911	', you may want to consider installing it'
912	print msg
913	else:
914	psyco.bind(Domain.update_boundary)
915	#psyco.bind(Domain.update_timestep) #Not worth it
916	psyco.bind(Domain.update_conserved_quantities)
917	psyco.bind(Domain.distribute_to_vertices_and_edges)
918
919
920	if __name__ == "__main__":
921	pass

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

Download in other formats: