/[escript]/trunk/escript/py_src/modelframe.py
ViewVC logotype

Contents of /trunk/escript/py_src/modelframe.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 950 - (show annotations)
Tue Feb 6 07:01:11 2007 UTC (12 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 40764 byte(s)
escript data objects can now be saved to netCDF files, see http://www.unidata.ucar.edu/software/netcdf/.
Currently only constant data are implemented with expanded and tagged data to follow.
There are two new functions to dump a data object

   s=Data(...)
   s.dump(<filename>)

and to recover it

   s=load(<filename>, domain)

Notice that the function space of s is recovered but domain is still need. 

dump and load will replace archive and extract.

The installation needs now the netCDF installed. 


1 # $Id$
2
3 """
4 Environment for implementing models in escript
5
6 @var __author__: name of author
7 @var __copyright__: copyrights
8 @var __license__: licence agreement
9 @var __url__: url entry point on documentation
10 @var __version__: version
11 @var __date__: date of the version
12 """
13
14 __author__="Lutz Gross, l.gross@uq.edu.au"
15 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
16 http://www.access.edu.au
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="http://www.iservo.edu.au/esys"
21 __version__="$Revision$"
22 __date__="$Date$"
23
24
25 from types import StringType,IntType,FloatType,BooleanType,ListType,DictType
26 from sys import stdout
27 import numarray
28 import operator
29 import itertools
30
31 # import the 'set' module if it's not defined (python2.3/2.4 difference)
32 try:
33 set
34 except NameError:
35 from sets import Set as set
36
37 from xml.dom import minidom
38
39
40 def all(seq):
41 for x in seq:
42 if not x:
43 return False
44 return True
45
46 def any(seq):
47 for x in seq:
48 if x:
49 return True
50 return False
51
52 def importName(modulename, name):
53 """ Import a named object from a module in the context of this function,
54 which means you should use fully qualified module paths.
55 Return None on failure.
56
57 This function from: http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/52241
58 """
59 module = __import__(modulename, globals(), locals(), [name])
60
61 try:
62 return vars(module)[name]
63 except KeyError:
64 raise ImportError("Could not import %s from %s" % (name, modulename))
65
66 class ESySXMLParser(object):
67 """
68 parser for ESysXML file
69 """
70 def __init__(self,xml, debug=False):
71 self.__dom = minidom.parseString(xml)
72 self.__linkable_object_registry= {}
73 self.__link_registry= []
74 self.__esys=self.__dom.getElementsByTagName('ESys')[0]
75 self.debug=debug
76
77 def getClassPath(self, node):
78 type = node.getAttribute("type")
79 if (node.getAttribute("module")):
80 module = node.getAttribute("module")
81 return importName(module, type)
82 else:
83 return importName("__main__", type)
84
85 def setLinks(self):
86 for obj_id, link in self.__link_registry:
87 link.target = self.__linkable_object_registry[obj_id]
88
89 def parse(self):
90 """
91 parser method for EsysXML and returns the list of generating ParameterSets
92 """
93 found=[]
94 for node in self.__esys.childNodes:
95 if isinstance(node, minidom.Element):
96 if node.tagName == 'Simulation':
97 found.append(Simulation.fromDom(self, node))
98 elif node.tagName == 'Model':
99 found.append(self.getClassPath(node).fromDom(self, node))
100 elif node.tagName == 'ParameterSet':
101 found.append(self.getClassPath(node).fromDom(self, node))
102 else:
103 raise "Invalid type, %r" % node.getAttribute("type")
104 self.setLinks()
105 return found
106
107 def registerLink(self,obj_id, link):
108 self.__link_registry.append((int(obj_id),link))
109
110 def registerLinkableObject(self,obj, node):
111 id_str=node.getAttribute('id').strip()
112 if len(id_str)>0:
113 id=int(id_str)
114 if self.__linkable_object_registry.has_key(id):
115 raise ValueError("Object id %s already exists."%id)
116 else:
117 self.__linkable_object_registry[id]=obj
118
119 def getComponent(self, node):
120 """
121 returns a single component + rank from a simulation
122 parser method for EsysXML and returns the list of generating ParameterSets
123 """
124 rank = int(node.getAttribute("rank"))
125 for n in node.childNodes:
126 if isinstance(n, minidom.Element):
127 if n.tagName == 'Simulation':
128 return (rank, Simulation.fromDom(self, n))
129 elif n.tagName == 'Model':
130 return (rank, self.getClassPath(n).fromDom(self, n))
131 elif n.tagName == 'ParameterSet':
132 return (rank, self.getClassPath(n).fromDom(self, n))
133 else:
134 raise ValueError("illegal component type %s"%n.tagName)
135 raise ValueError("cannot resolve Component")
136
137 class ESySXMLCreator(object):
138 """
139 creates an XML Dom representation
140 """
141 def __init__(self):
142 self.__dom=minidom.Document()
143 self.__esys =self.__dom.createElement('ESys')
144 self.__dom.appendChild(self.__esys)
145 self.__linkable_object_registry={}
146 self.__number_sequence = itertools.count(100)
147 def getRoot(self):
148 return self.__esys
149 def createElement(self,name):
150 return self.__dom.createElement(name)
151 def createTextNode(self,name):
152 return self.__dom.createTextNode(name)
153 def getElementById(self,name):
154 return self.__dom.getElementById(name)
155 def createDataNode(self, tagName, data):
156 """
157 C{createDataNode}s are the building blocks of the xml documents constructed in
158 this module.
159
160 @param tagName: the associated xml tag
161 @param data: the values in the tag
162 """
163 n = self.createElement(tagName)
164 n.appendChild(self.createTextNode(str(data)))
165 return n
166 def getLinkableObjectId(self, obj):
167 for id, o in self.__linkable_object_registry.items():
168 if o == obj: return id
169 id =self.__number_sequence.next()
170 self.__linkable_object_registry[id]=obj
171 return id
172
173 def registerLinkableObject(self, obj, node):
174 """
175 returns a unique object id for object obj
176 """
177 id=self.getLinkableObjectId(obj)
178 node.setAttribute('id',str(id))
179 node.setIdAttribute("id")
180
181 def includeTargets(self):
182 target_written=True
183 while target_written:
184 targetsList =self.__dom.getElementsByTagName('Target')
185 target_written=False
186 for element in targetsList:
187 targetId = int(element.firstChild.nodeValue.strip())
188 if self.getElementById(str(targetId)): continue
189 targetObj = self.__linkable_object_registry[targetId]
190 targetObj.toDom(self, self.__esys)
191 target_written=True
192
193 def toprettyxml(self):
194 self.includeTargets()
195 return self.__dom.toprettyxml()
196
197 class Link:
198 """
199 A Link makes an attribute of an object callable::
200
201 o.object()
202 o.a=8
203 l=Link(o,"a")
204 assert l()==8
205 """
206
207 def __init__(self,target,attribute=None):
208 """
209 Creates a link to the object target. If attribute is given, the link is
210 establised to this attribute of the target. Otherwise the attribute is
211 undefined.
212 """
213 self.target = target
214 self.attribute = None
215 self.setAttributeName(attribute)
216
217 def getTarget(self):
218 """
219 returns the target
220 """
221 return self.target
222 def getAttributeName(self):
223 """
224 returns the name of the attribute the link is pointing to
225 """
226 return self.attribute
227
228 def setAttributeName(self,attribute):
229 """
230 Set a new attribute name to be collected from the target object. The
231 target object must have the attribute with name attribute.
232 """
233 if attribute and self.target:
234 if isinstance(self.target,LinkableObject):
235 if not self.target.hasAttribute(attribute):
236 raise AttributeError("%s: target %s has no attribute %s."%(self, self.target, attribute))
237 else:
238 if not hasattr(self.target,attribute):
239 raise AttributeError("%s: target %s has no attribute %s."%(self, self.target, attribute))
240 self.attribute = attribute
241
242 def hasDefinedAttributeName(self):
243 """
244 Returns true if an attribute name is set.
245 """
246 return self.attribute != None
247
248 def __repr__(self):
249 """
250 Returns a string representation of the link.
251 """
252 if self.hasDefinedAttributeName():
253 return "<Link to attribute %s of %s>" % (self.attribute, self.target)
254 else:
255 return "<Link to target %s>" % self.target
256
257 def __call__(self,name=None):
258 """
259 Returns the value of the attribute of the target object. If the
260 atrribute is callable then the return value of the call is returned.
261 """
262 if name:
263 out=getattr(self.target, name)
264 else:
265 out=getattr(self.target, self.attribute)
266
267 if callable(out):
268 return out()
269 else:
270 return out
271
272 def toDom(self, esysxml, node):
273 """
274 C{toDom} method of Link. Creates a Link node and appends it to the
275 current XML esysxml.
276 """
277 link = esysxml.createElement('Link')
278 assert (self.target != None), ("Target was none, name was %r" % self.attribute)
279 link.appendChild(esysxml.createDataNode('Target', esysxml.getLinkableObjectId(self.target)))
280 # this use of id will not work for purposes of being able to retrieve the intended
281 # target from the xml later. I need a better unique identifier.
282 assert self.attribute, "You can't xmlify a Link without a target attribute"
283 link.appendChild(esysxml.createDataNode('Attribute', self.attribute))
284 node.appendChild(link)
285
286 def fromDom(cls, esysxml, node):
287 targetid = int(node.getElementsByTagName("Target")[0].firstChild.nodeValue.strip())
288 attribute =str(node.getElementsByTagName("Attribute")[0].firstChild.nodeValue.strip())
289 l = cls(None, attribute)
290 esysxml.registerLink(targetid, l)
291 return l
292
293 fromDom = classmethod(fromDom)
294
295 class LinkableObject(object):
296 """
297 An object that allows to link its attributes to attributes of other objects
298 via a Link object. For instance::
299
300 p = LinkableObject()
301 p.x = Link(o,"name")
302 print p.x
303
304 links attribute C{x} of C{p} to the attribute name of object C{o}.
305
306 C{p.x} will contain the current value of attribute C{name} of object
307 C{o}.
308
309 If the value of C{getattr(o, "name")} is callable, C{p.x} will return
310 the return value of the call.
311 """
312
313
314 def __init__(self, id = None, debug=False):
315 """
316 Initializes LinkableObject so that we can operate on Links
317 """
318 self.debug = debug
319 self.__linked_attributes={}
320
321 def trace(self, msg):
322 """
323 If debugging is on, print the message, otherwise do nothing
324 """
325 if self.debug:
326 print "%s: %s"%(str(self),msg)
327
328 def __getattr__(self,name):
329 """
330 Returns the value of attribute name. If the value is a Link object the
331 object is called and the return value is returned.
332 """
333 out = self.getAttributeObject(name)
334 if isinstance(out,Link):
335 return out()
336 else:
337 return out
338
339 def getAttributeObject(self,name):
340 """
341 Return the object stored for attribute name.
342 """
343
344 if self.__dict__.has_key(name):
345 return self.__dict__[name]
346
347 if self.__linked_attributes.has_key(name):
348 return self.__linked_attributes[name]
349
350 if self.__class__.__dict__.has_key(name):
351 return self.__class.__dict__[name]
352
353 raise AttributeError,"No attribute %s."%name
354
355 def hasAttribute(self,name):
356 """
357 Returns True if self as attribute name.
358 """
359 return self.__dict__.has_key(name) or self.__linked_attributes.has_key(name) or self.__class__.__dict__.has_key(name)
360
361 def __setattr__(self,name,value):
362 """
363 Sets the value for attribute name. If value is a Link the target
364 attribute is set to name if no attribute has been specified.
365 """
366
367 if self.__dict__.has_key(name):
368 del self.__dict__[name]
369
370 if isinstance(value,Link):
371 if not value.hasDefinedAttributeName():
372 value.setAttributeName(name)
373 self.__linked_attributes[name] = value
374
375 self.trace("attribute %s is now linked by %s."%(name,value))
376 else:
377 self.__dict__[name] = value
378
379 def __delattr__(self,name):
380 """
381 Removes the attribute name.
382 """
383
384 if self.__linked_attributes.has_key[name]:
385 del self.__linked_attributes[name]
386 elif self.__dict__.has_key(name):
387 del self.__dict__[name]
388 else:
389 raise AttributeError,"No attribute %s."%name
390
391 class _ParameterIterator:
392 def __init__(self,parameterset):
393
394 self.__set=parameterset
395 self.__iter=iter(parameterset.parameters)
396
397 def next(self):
398 o=self.__iter.next()
399 return (o,self.__set.getAttributeObject(o))
400
401 def __iter__(self):
402 return self
403
404 class ParameterSet(LinkableObject):
405 """
406 A class which allows to emphazise attributes to be written and read to XML
407
408 Leaves of an ESySParameters object can be:
409
410 - a real number
411 - a integer number
412 - a string
413 - a boolean value
414 - a ParameterSet object
415 - a Simulation object
416 - a Model object
417 - a numarray object
418 - a list of booleans
419 - any other object (not considered by writeESySXML and writeXML)
420
421 Example how to create an ESySParameters object::
422
423 p11=ParameterSet(gamma1=1.,gamma2=2.,gamma3=3.)
424 p1=ParameterSet(dim=2,tol_v=0.001,output_file="/tmp/u.%3.3d.dx",runFlag=True,parm11=p11)
425 parm=ParameterSet(parm1=p1,parm2=ParameterSet(alpha=Link(p11,"gamma1")))
426
427 This can be accessed as::
428
429 parm.parm1.gamma=0.
430 parm.parm1.dim=2
431 parm.parm1.tol_v=0.001
432 parm.parm1.output_file="/tmp/u.%3.3d.dx"
433 parm.parm1.runFlag=True
434 parm.parm1.parm11.gamma1=1.
435 parm.parm1.parm11.gamma2=2.
436 parm.parm1.parm11.gamma3=3.
437 parm.parm2.alpha=1. (value of parm.parm1.parm11.gamma1)
438 """
439 def __init__(self, parameters=[], **kwargs):
440 """
441 Creates a ParameterSet with parameters parameters.
442 """
443 LinkableObject.__init__(self, **kwargs)
444 self.parameters = set()
445 self.declareParameters(parameters)
446
447 def __repr__(self):
448 return "<%s %d>"%(self.__class__.__name__,id(self))
449
450 def declareParameter(self,**parameters):
451 """
452 Declares a new parameter(s) and its (their) initial value.
453 """
454 self.declareParameters(parameters)
455
456 def declareParameters(self,parameters):
457 """
458 Declares a set of parameters. parameters can be a list, a dictionary
459 or a ParameterSet.
460 """
461 if isinstance(parameters,ListType):
462 parameters = zip(parameters, itertools.repeat(None))
463 if isinstance(parameters,DictType):
464 parameters = parameters.iteritems()
465
466 for prm, value in parameters:
467 setattr(self,prm,value)
468 self.parameters.add(prm)
469
470 self.trace("parameter %s has been declared."%prm)
471
472 def releaseParameters(self,name):
473 """
474 Removes parameter name from the paramameters.
475 """
476 if self.isParameter(name):
477 self.parameters.remove(name)
478 self.trace("parameter %s has been removed."%name)
479
480 def checkLinkTargets(self, models, hash):
481 """
482 returns a list of (paramter set,parameter, target model ) if the the parameter of model
483 is linking to the target_model which is not in list of models.
484 """
485 out=[]
486 for name, value in self:
487 if isinstance(value, Link):
488 print str(self), name, value
489 m=value.getTarget()
490 if (m,name) in hash:
491 raise RuntimeError("recursive link: %s"%(hash+[ (m,name) ]))
492 if isinstance(m, Model) and not m in models: out.append( (self,name,m) )
493 if isinstance(m, ParameterSet):
494 try:
495 out+=m.checkLinkTargets(models, hash+[ (self,name) ] )
496 except RuntimeError, e:
497 print str(e), str(self), name
498 raise e
499 return list(set(out))
500
501 def __iter__(self):
502 """
503 Creates an iterator over the parameter and their values.
504 """
505 return _ParameterIterator(self)
506
507 def showParameters(self):
508 """
509 Returns a descrition of the parameters.
510 """
511 out="{"
512 notfirst=False
513 for i,v in self:
514 if notfirst: out=out+","
515 notfirst=True
516 if isinstance(v,ParameterSet):
517 out="%s\"%s\" : %s"%(out,i,v.showParameters())
518 else:
519 out="%s\"%s\" : %s"%(out,i,v)
520 return out+"}"
521
522 def __delattr__(self,name):
523 """
524 Removes the attribute name.
525 """
526 LinkableObject.__delattr__(self,name)
527 try:
528 self.releaseParameter(name)
529 except:
530 pass
531
532 def toDom(self, esysxml, node):
533 """
534 C{toDom} method of Model class
535 """
536 pset = esysxml.createElement('ParameterSet')
537 pset.setAttribute('type', self.__class__.__name__)
538 pset.setAttribute('module', self.__class__.__module__)
539 esysxml.registerLinkableObject(self, pset)
540 self._parametersToDom(esysxml, pset)
541 node.appendChild(pset)
542
543 def _parametersToDom(self, esysxml, node):
544 for name,value in self:
545 # convert list to numarray when possible:
546 if isinstance (value, list):
547 elem_type=-1
548 for i in value:
549 if isinstance(i, bool):
550 elem_type = max(elem_type,0)
551 if isinstance(i, int) and not isinstance(i, bool):
552 elem_type = max(elem_type,1)
553 if isinstance(i, float):
554 elem_type = max(elem_type,2)
555 if elem_type == 0: value = numarray.array(value,numarray.Bool)
556 if elem_type == 1: value = numarray.array(value,numarray.Int)
557 if elem_type == 2: value = numarray.array(value,numarray.Float)
558
559 param = esysxml.createElement('Parameter')
560 param.setAttribute('type', value.__class__.__name__)
561
562 param.appendChild(esysxml.createDataNode('Name', name))
563
564 val = esysxml.createElement('Value')
565 if isinstance(value,(ParameterSet,Link,DataSource)):
566 value.toDom(esysxml, val)
567 param.appendChild(val)
568 elif isinstance(value, numarray.NumArray):
569 shape = value.getshape()
570 if isinstance(shape, tuple):
571 size = reduce(operator.mul, shape)
572 shape = ' '.join(map(str, shape))
573 else:
574 size = shape
575 shape = str(shape)
576
577 arraytype = value.type()
578 if isinstance(arraytype, numarray.BooleanType):
579 arraytype_str="Bool"
580 elif isinstance(arraytype, numarray.IntegralType):
581 arraytype_str="Int"
582 elif isinstance(arraytype, numarray.FloatingType):
583 arraytype_str="Float"
584 elif isinstance(arraytype, numarray.ComplexType):
585 arraytype_str="Complex"
586 else:
587 arraytype_str=str(arraytype)
588 numarrayElement = esysxml.createElement('NumArray')
589 numarrayElement.appendChild(esysxml.createDataNode('ArrayType', arraytype_str))
590 numarrayElement.appendChild(esysxml.createDataNode('Shape', shape))
591 numarrayElement.appendChild(esysxml.createDataNode('Data', ' '.join(
592 [str(x) for x in numarray.reshape(value, size)])))
593 val.appendChild(numarrayElement)
594 param.appendChild(val)
595 elif isinstance(value, list):
596 param.appendChild(esysxml.createDataNode('Value', ' '.join([str(x) for x in value]) ))
597 elif isinstance(value, (str, bool, int, float, type(None))):
598 param.appendChild(esysxml.createDataNode('Value', str(value)))
599 elif isinstance(value, dict):
600 dic = esysxml.createElement('dictionary')
601 if len(value.keys())>0:
602 dic.setAttribute('key_type', value.keys()[0].__class__.__name__)
603 dic.setAttribute('value_type', value[value.keys()[0]].__class__.__name__)
604 for k,v in value.items():
605 i=esysxml.createElement('item')
606 i.appendChild(esysxml.createDataNode('key', k))
607 i.appendChild(esysxml.createDataNode('value', v))
608 dic.appendChild(i)
609 param.appendChild(dic)
610 else:
611 print value
612 raise ValueError("cannot serialize %s type to XML."%str(value.__class__))
613
614 node.appendChild(param)
615
616 def fromDom(cls, esysxml, node):
617 # Define a host of helper functions to assist us.
618 def _children(node):
619 """
620 Remove the empty nodes from the children of this node.
621 """
622 ret = []
623 for x in node.childNodes:
624 if isinstance(x, minidom.Text):
625 if x.nodeValue.strip():
626 ret.append(x)
627 else:
628 ret.append(x)
629 return ret
630
631 def _floatfromValue(esysxml, node):
632 return float(node.nodeValue.strip())
633
634 def _stringfromValue(esysxml, node):
635 return str(node.nodeValue.strip())
636
637 def _intfromValue(esysxml, node):
638 return int(node.nodeValue.strip())
639
640 def _boolfromValue(esysxml, node):
641 return _boolfromstring(node.nodeValue.strip())
642
643 def _nonefromValue(esysxml, node):
644 return None
645
646 def _numarrayfromValue(esysxml, node):
647 for node in _children(node):
648 if node.tagName == 'ArrayType':
649 arraytype = node.firstChild.nodeValue.strip()
650 if node.tagName == 'Shape':
651 shape = node.firstChild.nodeValue.strip()
652 shape = [int(x) for x in shape.split()]
653 if node.tagName == 'Data':
654 data = node.firstChild.nodeValue.strip()
655 data = [float(x) for x in data.split()]
656 return numarray.reshape(numarray.array(data, type=getattr(numarray, arraytype)),
657 shape)
658
659 def _listfromValue(esysxml, node):
660 return [x for x in node.nodeValue.split()]
661
662 def _boolfromstring(s):
663 if s == 'True':
664 return True
665 else:
666 return False
667 # Mapping from text types in the xml to methods used to process trees of that type
668 ptypemap = {"Simulation": Simulation.fromDom,
669 "Model":Model.fromDom,
670 "ParameterSet":ParameterSet.fromDom,
671 "Link":Link.fromDom,
672 "DataSource":DataSource.fromDom,
673 "float":_floatfromValue,
674 "int":_intfromValue,
675 "str":_stringfromValue,
676 "bool":_boolfromValue,
677 "list":_listfromValue,
678 "NumArray":_numarrayfromValue,
679 "NoneType":_nonefromValue,
680 }
681
682 parameters = {}
683 for n in _children(node):
684 ptype = n.getAttribute("type")
685 if not ptypemap.has_key(ptype):
686 raise KeyError("cannot handle parameter type %s."%ptype)
687
688 pname = pvalue = None
689 for childnode in _children(n):
690 if childnode.tagName == "Name":
691 pname = childnode.firstChild.nodeValue.strip()
692
693 if childnode.tagName == "Value":
694 nodes = _children(childnode)
695 pvalue = ptypemap[ptype](esysxml, nodes[0])
696
697 parameters[pname] = pvalue
698
699 # Create the instance of ParameterSet
700 o = cls(debug=esysxml.debug)
701 o.declareParameters(parameters)
702 esysxml.registerLinkableObject(o, node)
703 return o
704
705 fromDom = classmethod(fromDom)
706
707 def writeXML(self,ostream=stdout):
708 """
709 Writes the object as an XML object into an output stream.
710 """
711 esysxml=ESySXMLCreator()
712 self.toDom(esysxml, esysxml.getRoot())
713 ostream.write(esysxml.toprettyxml())
714
715 class Model(ParameterSet):
716 """
717 A Model object represents a processess marching over time until a
718 finalizing condition is fullfilled. At each time step an iterative
719 process can be performed and the time step size can be controlled. A
720 Model has the following work flow::
721
722 doInitialization()
723 while not terminateInitialIteration(): doInitializationiStep()
724 doInitialPostprocessing()
725 while not finalize():
726 dt=getSafeTimeStepSize(dt)
727 doStepPreprocessing(dt)
728 while not terminateIteration(): doStep(dt)
729 doStepPostprocessing(dt)
730 doFinalization()
731
732 where C{doInitialization}, C{finalize}, C{getSafeTimeStepSize},
733 C{doStepPreprocessing}, C{terminateIteration}, C{doStepPostprocessing},
734 C{doFinalization} are methods of the particular instance of a Model. The
735 default implementations of these methods have to be overwritten by the
736 subclass implementing a Model.
737 """
738
739 UNDEF_DT=1.e300
740
741 def __init__(self,parameters=[],**kwargs):
742 """
743 Creates a model.
744
745 Just calls the parent constructor.
746 """
747 ParameterSet.__init__(self, parameters=parameters,**kwargs)
748
749 def __str__(self):
750 return "<%s %d>"%(self.__class__.__name__,id(self))
751
752
753 def doInitialization(self):
754 """
755 Initializes the time stepping scheme.
756
757 This function may be overwritten.
758 """
759 pass
760 def doInitialStep(self):
761 """
762 performs an iteration step in the initialization phase
763
764 This function may be overwritten.
765 """
766 pass
767
768 def terminateInitialIteration(self):
769 """
770 Returns True if iteration at the inital phase is terminated.
771 """
772 return True
773
774 def doInitialPostprocessing(self):
775 """
776 finalises the initialization iteration process
777
778 This function may be overwritten.
779 """
780 pass
781
782 def getSafeTimeStepSize(self,dt):
783 """
784 Returns a time step size which can safely be used.
785
786 C{dt} gives the previously used step size.
787
788 This function may be overwritten.
789 """
790 return self.UNDEF_DT
791
792 def finalize(self):
793 """
794 Returns False if the time stepping is finalized.
795
796 This function may be overwritten.
797 """
798 return False
799
800 def doFinalization(self):
801 """
802 Finalizes the time stepping.
803
804 This function may be overwritten.
805 """
806 pass
807
808 def doStepPreprocessing(self,dt):
809 """
810 Sets up a time step of step size dt.
811
812 This function may be overwritten.
813 """
814 pass
815
816 def doStep(self,dt):
817 """
818 Executes an iteration step at a time step.
819
820 C{dt} is the currently used time step size.
821
822 This function may be overwritten.
823 """
824 pass
825
826 def terminateIteration(self):
827 """
828 Returns True if iteration on a time step is terminated.
829 """
830 return True
831
832
833 def doStepPostprocessing(self,dt):
834 """
835 finalises the time step.
836
837 dt is the currently used time step size.
838
839 This function may be overwritten.
840 """
841 pass
842
843 def toDom(self, esysxml, node):
844 """
845 C{toDom} method of Model class
846 """
847 pset = esysxml.createElement('Model')
848 pset.setAttribute('type', self.__class__.__name__)
849 pset.setAttribute('module', self.__class__.__module__)
850 esysxml.registerLinkableObject(self, pset)
851 node.appendChild(pset)
852 self._parametersToDom(esysxml, pset)
853
854 class Simulation(Model):
855 """
856 A Simulation object is special Model which runs a sequence of Models.
857
858 The methods C{doInitialization}, C{finalize}, C{getSafeTimeStepSize},
859 C{doStepPreprocessing}, C{terminateIteration}, C{doStepPostprocessing},
860 C{doFinalization} are executing the corresponding methods of the models in
861 the simulation.
862 """
863
864 FAILED_TIME_STEPS_MAX=20
865 MAX_ITER_STEPS=50
866 MAX_CHANGE_OF_DT=2.
867
868 def __init__(self, models=[], **kwargs):
869 """
870 Initiates a simulation from a list of models.
871 """
872 super(Simulation, self).__init__(**kwargs)
873 for m in models:
874 if not isinstance(m, Model):
875 raise TypeError("%s is not a subclass of Model."%m)
876 self.__models=[]
877 for i in range(len(models)):
878 self[i] = models[i]
879
880
881 def __repr__(self):
882 """
883 Returns a string representation of the Simulation.
884 """
885 return "<Simulation %r>" % self.__models
886
887 def __str__(self):
888 """
889 Returning Simulation as a string.
890 """
891 return "<Simulation %d>"%id(self)
892
893 def iterModels(self):
894 """
895 Returns an iterator over the models.
896 """
897 return self.__models
898
899 def __getitem__(self,i):
900 """
901 Returns the i-th model.
902 """
903 return self.__models[i]
904
905 def __setitem__(self,i,value):
906 """
907 Sets the i-th model.
908 """
909 if not isinstance(value,Model):
910 raise ValueError,"assigned value is not a Model but instance of %s"%(value.__class__.__name__,)
911 for j in range(max(i-len(self.__models)+1,0)):
912 self.__models.append(None)
913 self.__models[i]=value
914
915 def __len__(self):
916 """
917 Returns the number of models.
918 """
919 return len(self.__models)
920
921 def getAllModels(self):
922 """
923 returns a list of all models used in the Simulation including subsimulations
924 """
925 out=[]
926 for m in self.iterModels():
927 if isinstance(m, Simulation):
928 out+=m.getAllModels()
929 else:
930 out.append(m)
931 return list(set(out))
932
933 def checkModels(self, models, hash):
934 """
935 returns a list of (model,parameter, target model ) if the the parameter of model
936 is linking to the target_model which is not in list of models.
937 """
938 out=self.checkLinkTargets(models, hash)
939 for m in self.iterModels():
940 if isinstance(m, Simulation):
941 out+=[ (m,) + f for f in m.checkModels(models, hash) ]
942 else:
943 out+=[ (m,) + f for f in m.checkLinkTargets(models, hash) ]
944 return out
945
946
947 def getSafeTimeStepSize(self,dt):
948 """
949 Returns a time step size which can safely be used by all models.
950
951 This is the minimum over the time step sizes of all models.
952 """
953 out=min([o.getSafeTimeStepSize(dt) for o in self.iterModels()])
954 return out
955
956 def doInitialization(self):
957 """
958 Initializes all models.
959 """
960 self.n=0
961 self.tn=0.
962 for o in self.iterModels():
963 o.doInitialization()
964 def doInitialStep(self):
965 """
966 performs an iteration step in the initialization step for all models
967 """
968 iter=0
969 while not self.terminateInitialIteration():
970 if iter==0: self.trace("iteration for initialization starts")
971 iter+=1
972 self.trace("iteration step %d"%(iter))
973 for o in self.iterModels():
974 o.doInitialStep()
975 if iter>self.MAX_ITER_STEPS:
976 raise IterationDivergenceError("initial iteration did not converge after %s steps."%iter)
977 self.trace("Initialization finalized after %s iteration steps."%iter)
978
979 def doInitialPostprocessing(self):
980 """
981 finalises the initialization iteration process for all models.
982 """
983 for o in self.iterModels():
984 o.doInitialPostprocessing()
985 def finalize(self):
986 """
987 Returns True if any of the models is to be finalized.
988 """
989 return any([o.finalize() for o in self.iterModels()])
990
991 def doFinalization(self):
992 """
993 finalises the time stepping for all models.
994 """
995 for i in self.iterModels(): i.doFinalization()
996 self.trace("end of time integation.")
997
998 def doStepPreprocessing(self,dt):
999 """
1000 Initializes the time step for all models.
1001 """
1002 for o in self.iterModels():
1003 o.doStepPreprocessing(dt)
1004
1005 def terminateIteration(self):
1006 """
1007 Returns True if all iterations for all models are terminated.
1008 """
1009 out=all([o.terminateIteration() for o in self.iterModels()])
1010 return out
1011
1012 def terminateInitialIteration(self):
1013 """
1014 Returns True if all initial iterations for all models are terminated.
1015 """
1016 out=all([o.terminateInitialIteration() for o in self.iterModels()])
1017 return out
1018
1019 def doStepPostprocessing(self,dt):
1020 """
1021 finalises the iteration process for all models.
1022 """
1023 for o in self.iterModels():
1024 o.doStepPostprocessing(dt)
1025 self.n+=1
1026 self.tn+=dt
1027
1028 def doStep(self,dt):
1029 """
1030 Executes the iteration step at a time step for all model::
1031
1032 self.doStepPreprocessing(dt)
1033 while not self.terminateIteration():
1034 for all models:
1035 self.doStep(dt)
1036 self.doStepPostprocessing(dt)
1037 """
1038 self.iter=0
1039 while not self.terminateIteration():
1040 if self.iter==0: self.trace("iteration at %d-th time step %e starts"%(self.n+1,self.tn+dt))
1041 self.iter+=1
1042 self.trace("iteration step %d"%(self.iter))
1043 for o in self.iterModels():
1044 o.doStep(dt)
1045 if self.iter>0: self.trace("iteration at %d-th time step %e finalized."%(self.n+1,self.tn+dt))
1046
1047 def run(self,check_point=None):
1048 """
1049 Run the simulation by performing essentially::
1050
1051 self.doInitialization()
1052 while not self.terminateInitialIteration(): self.doInitialStep()
1053 self.doInitialPostprocessing()
1054 while not self.finalize():
1055 dt=self.getSafeTimeStepSize()
1056 self.doStep(dt)
1057 if n%check_point==0:
1058 self.writeXML()
1059 self.doFinalization()
1060
1061 If one of the models in throws a C{FailedTimeStepError} exception a
1062 new time step size is computed through getSafeTimeStepSize() and the
1063 time step is repeated.
1064
1065 If one of the models in throws a C{IterationDivergenceError}
1066 exception the time step size is halved and the time step is repeated.
1067
1068 In both cases the time integration is given up after
1069 C{Simulation.FAILED_TIME_STEPS_MAX} attempts.
1070 """
1071 # check the completness of the models:
1072 # first a list of all the models involved in the simulation including subsimulations:
1073 #
1074 missing=self.checkModels(self.getAllModels(), [])
1075 if len(missing)>0:
1076 msg=""
1077 for l in missing:
1078 msg+="\n\t"+str(l[-1])+" at "+str(self)
1079 for i in xrange(len(l)-1): msg+="."+str(l[i])
1080 raise MissingLink("link targets missing in the Simulation: %s"%msg)
1081 #==============================
1082 self.doInitialization()
1083 self.doInitialStep()
1084 self.doInitialPostprocessing()
1085 dt=self.UNDEF_DT
1086 while not self.finalize():
1087 step_fail_counter=0
1088 iteration_fail_counter=0
1089 if self.n==0:
1090 dt_new=self.getSafeTimeStepSize(dt)
1091 else:
1092 dt_new=min(max(self.getSafeTimeStepSize(dt),dt/self.MAX_CHANGE_OF_DT),dt*self.MAX_CHANGE_OF_DT)
1093 self.trace("%d. time step %e (step size %e.)" % (self.n+1,self.tn+dt_new,dt_new))
1094 end_of_step=False
1095 while not end_of_step:
1096 end_of_step=True
1097 if not dt_new>0:
1098 raise NonPositiveStepSizeError("non-positive step size in step %d"%(self.n+1))
1099 try:
1100 self.doStepPreprocessing(dt_new)
1101 self.doStep(dt_new)
1102 self.doStepPostprocessing(dt_new)
1103 except IterationDivergenceError:
1104 dt_new*=0.5
1105 end_of_step=False
1106 iteration_fail_counter+=1
1107 if iteration_fail_counter>self.FAILED_TIME_STEPS_MAX:
1108 raise SimulationBreakDownError("reduction of time step to achieve convergence failed after %s steps."%self.FAILED_TIME_STEPS_MAX)
1109 self.trace("Iteration failed. Time step is repeated with new step size %s."%dt_new)
1110 except FailedTimeStepError:
1111 dt_new=self.getSafeTimeStepSize(dt)
1112 end_of_step=False
1113 step_fail_counter+=1
1114 self.trace("Time step is repeated with new time step size %s."%dt_new)
1115 if step_fail_counter>self.FAILED_TIME_STEPS_MAX:
1116 raise SimulationBreakDownError("Time integration is given up after %d attempts."%step_fail_counter)
1117 dt=dt_new
1118 if not check_point==None:
1119 if n%check_point==0:
1120 self.trace("check point is created.")
1121 self.writeXML()
1122 self.doFinalization()
1123
1124
1125 def toDom(self, esysxml, node):
1126 """
1127 C{toDom} method of Simulation class.
1128 """
1129 simulation = esysxml.createElement('Simulation')
1130 esysxml.registerLinkableObject(self, simulation)
1131 for rank, sim in enumerate(self.iterModels()):
1132 component = esysxml.createElement('Component')
1133 component.setAttribute('rank', str(rank))
1134 sim.toDom(esysxml, component)
1135 simulation.appendChild(component)
1136 node.appendChild(simulation)
1137
1138
1139 def fromDom(cls, esysxml, node):
1140 sims = []
1141 for n in node.childNodes:
1142 if isinstance(n, minidom.Text):
1143 continue
1144 sims.append(esysxml.getComponent(n))
1145 sims.sort(cmp=_comp)
1146 sim=cls([s[1] for s in sims], debug=esysxml.debug)
1147 esysxml.registerLinkableObject(sim, node)
1148 return sim
1149
1150 fromDom = classmethod(fromDom)
1151
1152 def _comp(a,b):
1153 if a[0]<a[1]:
1154 return 1
1155 elif a[0]>a[1]:
1156 return -1
1157 else:
1158 return 0
1159
1160 class IterationDivergenceError(Exception):
1161 """
1162 Exception which is thrown if there is no convergence of the iteration
1163 process at a time step.
1164
1165 But there is a chance that a smaller step could help to reach convergence.
1166 """
1167 pass
1168
1169 class FailedTimeStepError(Exception):
1170 """
1171 Exception which is thrown if the time step fails because of a step
1172 size that have been choosen to be too large.
1173 """
1174 pass
1175
1176 class SimulationBreakDownError(Exception):
1177 """
1178 Exception which is thrown if the simulation does not manage to
1179 progress in time.
1180 """
1181 pass
1182
1183 class NonPositiveStepSizeError(Exception):
1184 """
1185 Exception which is thrown if the step size is not positive.
1186 """
1187 pass
1188
1189 class MissingLink(Exception):
1190 """
1191 Exception thrown when a link is missing
1192 """
1193 pass
1194
1195 class DataSource(object):
1196 """
1197 Class for handling data sources, including local and remote files. This class is under development.
1198 """
1199
1200 def __init__(self, uri="file.ext", fileformat="unknown"):
1201 self.uri = uri
1202 self.fileformat = fileformat
1203
1204 def toDom(self, esysxml, node):
1205 """
1206 C{toDom} method of DataSource. Creates a DataSource node and appends it to the
1207 current XML esysxml.
1208 """
1209 ds = esysxml.createElement('DataSource')
1210 ds.appendChild(esysxml.createDataNode('URI', self.uri))
1211 ds.appendChild(esysxml.createDataNode('FileFormat', self.fileformat))
1212 node.appendChild(ds)
1213
1214 def fromDom(cls, esysxml, node):
1215 uri= str(node.getElementsByTagName("URI")[0].firstChild.nodeValue.strip())
1216 fileformat= str(node.getElementsByTagName("FileFormat")[0].firstChild.nodeValue.strip())
1217 ds = cls(uri, fileformat)
1218 return ds
1219
1220 def getLocalFileName(self):
1221 return self.uri
1222
1223 fromDom = classmethod(fromDom)
1224
1225 # vim: expandtab shiftwidth=4:

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26