/[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 874 - (show annotations)
Tue Oct 17 12:06:11 2006 UTC (13 years ago) by elspeth
File MIME type: text/x-python
File size: 31714 byte(s)
Changed numarray xml structure slightly, to fit better with the rest of 
the XML. 
Added lists of booleans as a data type in modelframe.
Fixed bug with bools in modelframe.

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 # import modellib temporarily removed!!!
31
32 # import the 'set' module if it's not defined (python2.3/2.4 difference)
33 try:
34 set
35 except NameError:
36 from sets import Set as set
37
38 from xml.dom import minidom
39
40 def dataNode(document, tagName, data):
41 """
42 C{dataNode}s are the building blocks of the xml documents constructed in
43 this module.
44
45 @param document: the current xml document
46 @param tagName: the associated xml tag
47 @param data: the values in the tag
48 """
49 t = document.createTextNode(str(data))
50 n = document.createElement(tagName)
51 n.appendChild(t)
52 return n
53
54 def esysDoc():
55 """
56 Global method for creating an instance of an EsysXML document.
57 """
58 doc = minidom.Document()
59 esys = doc.createElement('ESys')
60 doc.appendChild(esys)
61 return doc, esys
62
63 def all(seq):
64 for x in seq:
65 if not x:
66 return False
67 return True
68
69 def any(seq):
70 for x in seq:
71 if x:
72 return True
73 return False
74
75 LinkableObjectRegistry = {}
76
77 def registerLinkableObject(obj_id, o):
78 LinkableObjectRegistry[obj_id] = o
79
80 LinkRegistry = []
81
82 def registerLink(obj_id, l):
83 LinkRegistry.append((obj_id,l))
84
85 def parse(xml):
86 """
87 Generic parse method for EsysXML. Without this, Links don't work.
88 """
89 global LinkRegistry, LinkableObjectRegistry
90 LinkRegistry = []
91 LinkableObjectRegistry = {}
92
93 doc = minidom.parseString(xml)
94 sim = getComponent(doc.firstChild)
95 for obj_id, link in LinkRegistry:
96 link.target = LinkableObjectRegistry[obj_id]
97
98 return sim
99
100 def importName(modulename, name):
101 """ Import a named object from a module in the context of this function,
102 which means you should use fully qualified module paths.
103
104 Return None on failure.
105
106 This function from: http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/52241
107 """
108 module = __import__(modulename, globals(), locals(), [name])
109
110 try:
111 return vars(module)[name]
112 except KeyError:
113 raise ImportError("Could not import %s from %s" % (name, modulename))
114
115 def getComponent(doc):
116 """
117 Used to get components of Simualtions, Models.
118 """
119 for node in doc.childNodes:
120
121 if isinstance(node, minidom.Element):
122 if node.tagName == 'Simulation':
123 if node.getAttribute("type") == 'Simulation':
124 return Simulation.fromDom(node)
125 if node.tagName == 'Model':
126 if (node.getAttribute("module")):
127 model_module = node.getAttribute("module")
128 model_type = node.getAttribute("type")
129 return importName(model_module, model_type).fromDom(node)
130 else:
131 model_type = node.getAttribute("type")
132 model_subclasses = Model.__subclasses__()
133 for model in model_subclasses:
134 if model_type == model.__name__:
135 return Model.fromDom(node)
136 if node.tagName == 'ParameterSet':
137 parameter_type = node.getAttribute("type")
138 return ParameterSet.fromDom(node)
139 raise "Invalid simulation type, %r" % node.getAttribute("type")
140
141
142 raise ValueError("No Simulation Found")
143
144
145 class Link:
146 """
147 A Link makes an attribute of an object callable::
148
149 o.object()
150 o.a=8
151 l=Link(o,"a")
152 assert l()==8
153 """
154
155 def __init__(self,target,attribute=None):
156 """
157 Creates a link to the object target. If attribute is given, the link is
158 establised to this attribute of the target. Otherwise the attribute is
159 undefined.
160 """
161 self.target = target
162 self.attribute = None
163 self.setAttributeName(attribute)
164
165 def setAttributeName(self,attribute):
166 """
167 Set a new attribute name to be collected from the target object. The
168 target object must have the attribute with name attribute.
169 """
170 if attribute and self.target:
171 if isinstance(self.target,LinkableObject):
172 if not self.target.hasAttribute(attribute):
173 raise AttributeError("%s: target %s has no attribute %s."%(self, self.target, attribute))
174 else:
175 if not hasattr(self.target,attribute):
176 raise AttributeError("%s: target %s has no attribute %s."%(self, self.target, attribute))
177 self.attribute = attribute
178
179 def hasDefinedAttributeName(self):
180 """
181 Returns true if an attribute name is set.
182 """
183 return self.attribute != None
184
185 def __repr__(self):
186 """
187 Returns a string representation of the link.
188 """
189 if self.hasDefinedAttributeName():
190 return "<Link to attribute %s of %s>" % (self.attribute, self.target)
191 else:
192 return "<Link to target %s>" % self.target
193
194 def __call__(self,name=None):
195 """
196 Returns the value of the attribute of the target object. If the
197 atrribute is callable then the return value of the call is returned.
198 """
199 if name:
200 out=getattr(self.target, name)
201 else:
202 out=getattr(self.target, self.attribute)
203
204 if callable(out):
205 return out()
206 else:
207 return out
208
209 def toDom(self, document, node):
210 """
211 C{toDom} method of Link. Creates a Link node and appends it to the
212 current XML document.
213 """
214 link = document.createElement('Link')
215 assert (self.target != None), ("Target was none, name was %r" % self.attribute)
216 link.appendChild(dataNode(document, 'Target', self.target.id))
217 # this use of id will not work for purposes of being able to retrieve the intended
218 # target from the xml later. I need a better unique identifier.
219 assert self.attribute, "You can't xmlify a Link without a target attribute"
220 link.appendChild(dataNode(document, 'Attribute', self.attribute))
221 node.appendChild(link)
222
223 def fromDom(cls, doc):
224 targetid = doc.getElementsByTagName("Target")[0].firstChild.nodeValue.strip()
225 attribute = doc.getElementsByTagName("Attribute")[0].firstChild.nodeValue.strip()
226 l = cls(None, attribute)
227 registerLink(targetid, l)
228 return l
229
230 fromDom = classmethod(fromDom)
231
232 def writeXML(self,ostream=stdout):
233 """
234 Writes an XML representation of self to the output stream ostream.
235 If ostream is nor present the standart output stream is used. If
236 esysheader==True the esys XML header is written
237 """
238 print 'I got to the Link writeXML method'
239 document, rootnode = esysDoc()
240 self.toDom(document, rootnode)
241
242 ostream.write(document.toprettyxml())
243
244 class LinkableObject(object):
245 """
246 An object that allows to link its attributes to attributes of other objects
247 via a Link object. For instance::
248
249 p = LinkableObject()
250 p.x = Link(o,"name")
251 print p.x
252
253 links attribute C{x} of C{p} to the attribute name of object C{o}.
254
255 C{p.x} will contain the current value of attribute C{name} of object
256 C{o}.
257
258 If the value of C{getattr(o, "name")} is callable, C{p.x} will return
259 the return value of the call.
260 """
261
262 number_sequence = itertools.count(100)
263
264 def __init__(self, debug=False):
265 """
266 Initializes LinkableObject so that we can operate on Links
267 """
268 self.debug = debug
269 self.__linked_attributes={}
270 self.id = self.number_sequence.next()
271 registerLinkableObject(self.id, self)
272
273 def trace(self, msg):
274 """
275 If debugging is on, print the message, otherwise do nothing
276 """
277 if self.debug:
278 print "%s: %s"%(str(self),msg)
279
280 def __getattr__(self,name):
281 """
282 Returns the value of attribute name. If the value is a Link object the
283 object is called and the return value is returned.
284 """
285 out = self.getAttributeObject(name)
286 if isinstance(out,Link):
287 return out()
288 else:
289 return out
290
291 def getAttributeObject(self,name):
292 """
293 Return the object stored for attribute name.
294 """
295
296 if self.__dict__.has_key(name):
297 return self.__dict__[name]
298
299 if self.__linked_attributes.has_key(name):
300 return self.__linked_attributes[name]
301
302 if self.__class__.__dict__.has_key(name):
303 return self.__class.__dict__[name]
304
305 raise AttributeError,"No attribute %s."%name
306
307 def hasAttribute(self,name):
308 """
309 Returns True if self as attribute name.
310 """
311 return self.__dict__.has_key(name) or self.__linked_attributes.has_key(name) or self.__class__.__dict__.has_key(name)
312
313 def __setattr__(self,name,value):
314 """
315 Sets the value for attribute name. If value is a Link the target
316 attribute is set to name if no attribute has been specified.
317 """
318
319 if self.__dict__.has_key(name):
320 del self.__dict__[name]
321
322 if isinstance(value,Link):
323 if not value.hasDefinedAttributeName():
324 value.setAttributeName(name)
325 self.__linked_attributes[name] = value
326
327 self.trace("attribute %s is now linked by %s."%(name,value))
328 else:
329 self.__dict__[name] = value
330
331 def __delattr__(self,name):
332 """
333 Removes the attribute name.
334 """
335
336 if self.__linked_attributes.has_key[name]:
337 del self.__linked_attributes[name]
338 elif self.__dict__.has_key(name):
339 del self.__dict__[name]
340 else:
341 raise AttributeError,"No attribute %s."%name
342
343 class _ParameterIterator:
344 def __init__(self,parameterset):
345
346 self.__set=parameterset
347 self.__iter=iter(parameterset.parameters)
348
349 def next(self):
350 o=self.__iter.next()
351 return (o,self.__set.getAttributeObject(o))
352
353 def __iter__(self):
354 return self
355
356 class ParameterSet(LinkableObject):
357 """
358 A class which allows to emphazise attributes to be written and read to XML
359
360 Leaves of an ESySParameters object can be:
361
362 - a real number
363 - a integer number
364 - a string
365 - a boolean value
366 - a ParameterSet object
367 - a Simulation object
368 - a Model object
369 - a numarray object
370 - a list of booleans
371 - any other object (not considered by writeESySXML and writeXML)
372
373 Example how to create an ESySParameters object::
374
375 p11=ParameterSet(gamma1=1.,gamma2=2.,gamma3=3.)
376 p1=ParameterSet(dim=2,tol_v=0.001,output_file="/tmp/u.%3.3d.dx",runFlag=True,parm11=p11)
377 parm=ParameterSet(parm1=p1,parm2=ParameterSet(alpha=Link(p11,"gamma1")))
378
379 This can be accessed as::
380
381 parm.parm1.gamma=0.
382 parm.parm1.dim=2
383 parm.parm1.tol_v=0.001
384 parm.parm1.output_file="/tmp/u.%3.3d.dx"
385 parm.parm1.runFlag=True
386 parm.parm1.parm11.gamma1=1.
387 parm.parm1.parm11.gamma2=2.
388 parm.parm1.parm11.gamma3=3.
389 parm.parm2.alpha=1. (value of parm.parm1.parm11.gamma1)
390 """
391 def __init__(self, parameters=[], **kwargs):
392 """
393 Creates a ParameterSet with parameters parameters.
394 """
395 LinkableObject.__init__(self, **kwargs)
396 self.parameters = set()
397 self.declareParameters(parameters)
398
399 def __repr__(self):
400 return "<%s %r>" % (self.__class__.__name__,
401 [(p, getattr(self, p, None)) for p in self.parameters])
402
403 def declareParameter(self,**parameters):
404 """
405 Declares a new parameter(s) and its (their) initial value.
406 """
407 self.declareParameters(parameters)
408
409 def declareParameters(self,parameters):
410 """
411 Declares a set of parameters. parameters can be a list, a dictionary
412 or a ParameterSet.
413 """
414 if isinstance(parameters,ListType):
415 parameters = zip(parameters, itertools.repeat(None))
416 if isinstance(parameters,DictType):
417 parameters = parameters.iteritems()
418
419 for prm, value in parameters:
420 setattr(self,prm,value)
421 self.parameters.add(prm)
422
423 self.trace("parameter %s has been declared."%prm)
424
425 def releaseParameters(self,name):
426 """
427 Removes parameter name from the paramameters.
428 """
429 if self.isParameter(name):
430 self.parameters.remove(name)
431 self.trace("parameter %s has been removed."%name)
432
433 def __iter__(self):
434 """
435 Creates an iterator over the parameter and their values.
436 """
437 return _ParameterIterator(self)
438
439 def showParameters(self):
440 """
441 Returns a descrition of the parameters.
442 """
443 out="{"
444 notfirst=False
445 for i,v in self:
446 if notfirst: out=out+","
447 notfirst=True
448 if isinstance(v,ParameterSet):
449 out="%s\"%s\" : %s"%(out,i,v.showParameters())
450 else:
451 out="%s\"%s\" : %s"%(out,i,v)
452 return out+"}"
453
454 def __delattr__(self,name):
455 """
456 Removes the attribute name.
457 """
458 LinkableObject.__delattr__(self,name)
459 try:
460 self.releaseParameter(name)
461 except:
462 pass
463
464 def toDom(self, document, node):
465 """
466 C{toDom} method of ParameterSet class.
467 """
468 pset = document.createElement('ParameterSet')
469 node.appendChild(pset)
470 self._parametersToDom(document, pset)
471
472 def _parametersToDom(self, document, node):
473 node.setAttribute('id', str(self.id))
474 node.setIdAttribute("id")
475 for name,value in self:
476 param = document.createElement('Parameter')
477 param.setAttribute('type', value.__class__.__name__)
478
479 param.appendChild(dataNode(document, 'Name', name))
480
481 val = document.createElement('Value')
482
483 if isinstance(value,(ParameterSet,Link)):
484 value.toDom(document, val)
485 param.appendChild(val)
486 elif isinstance(value, numarray.NumArray):
487 shape = value.getshape()
488 if isinstance(shape, tuple):
489 size = reduce(operator.mul, shape)
490 shape = ' '.join(map(str, shape))
491 else:
492 size = shape
493 shape = str(shape)
494
495 arraytype = value.type()
496 numarrayElement = document.createElement('NumArray')
497 numarrayElement.appendChild(dataNode(document, 'ArrayType', str(arraytype)))
498 numarrayElement.appendChild(dataNode(document, 'Shape', shape))
499 numarrayElement.appendChild(dataNode(document, 'Data', ' '.join(
500 [str(x) for x in numarray.reshape(value, size)])))
501 val.appendChild(numarrayElement)
502 param.appendChild(val)
503 elif isinstance (value, list):
504 param.appendChild(dataNode(document, 'Value', ' '.join(
505 [str(x) for x in value])
506 ))
507 else:
508 param.appendChild(dataNode(document, 'Value', str(value)))
509
510 node.appendChild(param)
511
512 def fromDom(cls, doc):
513
514 # Define a host of helper functions to assist us.
515 def _children(node):
516 """
517 Remove the empty nodes from the children of this node.
518 """
519 ret = []
520 for x in node.childNodes:
521 if isinstance(x, minidom.Text):
522 if x.nodeValue.strip():
523 ret.append(x)
524 else:
525 ret.append(x)
526 return ret
527
528 def _floatfromValue(doc):
529 return float(doc.nodeValue.strip())
530
531 def _stringfromValue(doc):
532 return str(doc.nodeValue.strip())
533
534 def _intfromValue(doc):
535 return int(doc.nodeValue.strip())
536
537 def _boolfromValue(doc):
538 return _boolfromstring(doc.nodeValue.strip())
539
540 def _nonefromValue(doc):
541 return None
542
543 def _numarrayfromValue(doc):
544 for node in _children(doc):
545 if node.tagName == 'ArrayType':
546 arraytype = node.firstChild.nodeValue.strip()
547 if node.tagName == 'Shape':
548 shape = node.firstChild.nodeValue.strip()
549 shape = [int(x) for x in shape.split()]
550 if node.tagName == 'Data':
551 data = node.firstChild.nodeValue.strip()
552 data = [float(x) for x in data.split()]
553 return numarray.reshape(numarray.array(data, type=getattr(numarray, arraytype)),
554 shape)
555
556 def _listfromValue(doc):
557 return [_boolfromstring(x) for x in doc.nodeValue.split()]
558
559
560 def _boolfromstring(s):
561 if s == 'True':
562 return True
563 else:
564 return False
565 # Mapping from text types in the xml to methods used to process trees of that type
566 ptypemap = {"Simulation": Simulation.fromDom,
567 "Model":Model.fromDom,
568 "ParameterSet":ParameterSet.fromDom,
569 "Link":Link.fromDom,
570 "float":_floatfromValue,
571 "int":_intfromValue,
572 "str":_stringfromValue,
573 "bool":_boolfromValue,
574 "list":_listfromValue,
575 "NumArray":_numarrayfromValue,
576 "NoneType":_nonefromValue,
577 }
578
579 # print doc.toxml()
580
581 parameters = {}
582 for node in _children(doc):
583 ptype = node.getAttribute("type")
584
585 pname = pvalue = None
586 for childnode in _children(node):
587
588 if childnode.tagName == "Name":
589 pname = childnode.firstChild.nodeValue.strip()
590
591 if childnode.tagName == "Value":
592 nodes = _children(childnode)
593 # if ptype == 'NumArray':
594 # pvalue = _numarrayfromValue(nodes)
595 # else:
596 pvalue = ptypemap[ptype](nodes[0])
597
598 parameters[pname] = pvalue
599
600 # Create the instance of ParameterSet
601 o = cls()
602 o.declareParameters(parameters)
603 registerLinkableObject(doc.getAttribute("id"), o)
604 return o
605
606 fromDom = classmethod(fromDom)
607
608 def writeXML(self,ostream=stdout):
609 """
610 Writes the object as an XML object into an output stream.
611 """
612 # ParameterSet(d) with d[Name]=Value
613 document, node = esysDoc()
614 self.toDom(document, node)
615 ostream.write(document.toprettyxml())
616
617 class Model(ParameterSet):
618 """
619 A Model object represents a processess marching over time until a
620 finalizing condition is fullfilled. At each time step an iterative
621 process can be performed and the time step size can be controlled. A
622 Model has the following work flow::
623
624 doInitialization()
625 while not finalize():
626 dt=getSafeTimeStepSize(dt)
627 doStepPreprocessing(dt)
628 while not terminateIteration(): doStep(dt)
629 doStepPostprocessing(dt)
630 doFinalization()
631
632 where C{doInitialization}, C{finalize}, C{getSafeTimeStepSize},
633 C{doStepPreprocessing}, C{terminateIteration}, C{doStepPostprocessing},
634 C{doFinalization} are methods of the particular instance of a Model. The
635 default implementations of these methods have to be overwritten by the
636 subclass implementing a Model.
637 """
638
639 UNDEF_DT=1.e300
640
641 def __init__(self,parameters=[],**kwarg):
642 """
643 Creates a model.
644
645 Just calls the parent constructor.
646 """
647 ParameterSet.__init__(self, parameters=parameters,**kwarg)
648
649 def __str__(self):
650 return "<%s %d>"%(self.__class__,id(self))
651
652 def toDom(self, document, node):
653 """
654 C{toDom} method of Model class
655 """
656 pset = document.createElement('Model')
657 pset.setAttribute('type', self.__class__.__name__)
658 if not self.__class__.__module__.startswith('esys.escript'):
659 pset.setAttribute('module', self.__class__.__module__)
660 node.appendChild(pset)
661 self._parametersToDom(document, pset)
662
663 def doInitialization(self):
664 """
665 Initializes the time stepping scheme.
666
667 This function may be overwritten.
668 """
669 pass
670
671 def getSafeTimeStepSize(self,dt):
672 """
673 Returns a time step size which can safely be used.
674
675 C{dt} gives the previously used step size.
676
677 This function may be overwritten.
678 """
679 return self.UNDEF_DT
680
681 def finalize(self):
682 """
683 Returns False if the time stepping is finalized.
684
685 This function may be overwritten.
686 """
687 return False
688
689 def doFinalization(self):
690 """
691 Finalizes the time stepping.
692
693 This function may be overwritten.
694 """
695 pass
696
697 def doStepPreprocessing(self,dt):
698 """
699 Sets up a time step of step size dt.
700
701 This function may be overwritten.
702 """
703 pass
704
705 def doStep(self,dt):
706 """
707 Executes an iteration step at a time step.
708
709 C{dt} is the currently used time step size.
710
711 This function may be overwritten.
712 """
713 pass
714
715 def terminateIteration(self):
716 """
717 Returns True if iteration on a time step is terminated.
718 """
719 return True
720
721 def doStepPostprocessing(self,dt):
722 """
723 Finalalizes the time step.
724
725 dt is the currently used time step size.
726
727 This function may be overwritten.
728 """
729 pass
730
731 def writeXML(self, ostream=stdout):
732 document, node = esysDoc()
733 self.toDom(document, node)
734 ostream.write(document.toprettyxml())
735
736
737 class Simulation(Model):
738 """
739 A Simulation object is special Model which runs a sequence of Models.
740
741 The methods C{doInitialization}, C{finalize}, C{getSafeTimeStepSize},
742 C{doStepPreprocessing}, C{terminateIteration}, C{doStepPostprocessing},
743 C{doFinalization} are executing the corresponding methods of the models in
744 the simulation.
745 """
746
747 FAILED_TIME_STEPS_MAX=20
748 MAX_ITER_STEPS=50
749 MAX_CHANGE_OF_DT=2.
750
751 def __init__(self, models=[], **kwargs):
752 """
753 Initiates a simulation from a list of models.
754 """
755 Model.__init__(self, **kwargs)
756 self.__models=[]
757
758 for i in range(len(models)):
759 self[i] = models[i]
760
761
762 def __repr__(self):
763 """
764 Returns a string representation of the Simulation.
765 """
766 return "<Simulation %r>" % self.__models
767
768 def __str__(self):
769 """
770 Returning Simulation as a string.
771 """
772 return "<Simulation %d>"%id(self)
773
774 def iterModels(self):
775 """
776 Returns an iterator over the models.
777 """
778 return self.__models
779
780 def __getitem__(self,i):
781 """
782 Returns the i-th model.
783 """
784 return self.__models[i]
785
786 def __setitem__(self,i,value):
787 """
788 Sets the i-th model.
789 """
790 if not isinstance(value,Model):
791 raise ValueError,"assigned value is not a Model but instance of %s"%(value.__class__.__name__,)
792 for j in range(max(i-len(self.__models)+1,0)):
793 self.__models.append(None)
794 self.__models[i]=value
795
796 def __len__(self):
797 """
798 Returns the number of models.
799 """
800 return len(self.__models)
801
802 def toDom(self, document, node):
803 """
804 C{toDom} method of Simulation class.
805 """
806 simulation = document.createElement('Simulation')
807 simulation.setAttribute('type', self.__class__.__name__)
808
809 for rank, sim in enumerate(self.iterModels()):
810 component = document.createElement('Component')
811 component.setAttribute('rank', str(rank))
812
813 sim.toDom(document, component)
814
815 simulation.appendChild(component)
816
817 node.appendChild(simulation)
818
819 def writeXML(self,ostream=stdout):
820 """
821 Writes the object as an XML object into an output stream.
822 """
823 document, rootnode = esysDoc()
824 self.toDom(document, rootnode)
825 targetsList = document.getElementsByTagName('Target')
826
827 for element in targetsList:
828 targetId = int(element.firstChild.nodeValue.strip())
829 if document.getElementById(str(targetId)):
830 continue
831 targetObj = LinkableObjectRegistry[targetId]
832 targetObj.toDom(document, rootnode)
833 ostream.write(document.toprettyxml())
834
835 def getSafeTimeStepSize(self,dt):
836 """
837 Returns a time step size which can safely be used by all models.
838
839 This is the minimum over the time step sizes of all models.
840 """
841 out=min([o.getSafeTimeStepSize(dt) for o in self.iterModels()])
842 #print "%s: safe step size is %e."%(str(self),out)
843 return out
844
845 def doInitialization(self):
846 """
847 Initializes all models.
848 """
849 self.n=0
850 self.tn=0.
851 for o in self.iterModels():
852 o.doInitialization()
853
854 def finalize(self):
855 """
856 Returns True if any of the models is to be finalized.
857 """
858 return any([o.finalize() for o in self.iterModels()])
859
860 def doFinalization(self):
861 """
862 Finalalizes the time stepping for all models.
863 """
864 for i in self.iterModels(): i.doFinalization()
865 self.trace("end of time integation.")
866
867 def doStepPreprocessing(self,dt):
868 """
869 Initializes the time step for all models.
870 """
871 for o in self.iterModels():
872 o.doStepPreprocessing(dt)
873
874 def terminateIteration(self):
875 """
876 Returns True if all iterations for all models are terminated.
877 """
878 out=all([o.terminateIteration() for o in self.iterModels()])
879 return out
880
881 def doStepPostprocessing(self,dt):
882 """
883 Finalalizes the iteration process for all models.
884 """
885 for o in self.iterModels():
886 o.doStepPostprocessing(dt)
887 self.n+=1
888 self.tn+=dt
889
890 def doStep(self,dt):
891 """
892 Executes the iteration step at a time step for all model::
893
894 self.doStepPreprocessing(dt)
895 while not self.terminateIteration():
896 for all models:
897 self.doStep(dt)
898 self.doStepPostprocessing(dt)
899 """
900 self.iter=0
901 while not self.terminateIteration():
902 if self.iter==0: self.trace("iteration at %d-th time step %e starts"%(self.n+1,self.tn+dt))
903 self.iter+=1
904 self.trace("iteration step %d"%(self.iter))
905 for o in self.iterModels():
906 o.doStep(dt)
907 if self.iter>0: self.trace("iteration at %d-th time step %e finalized."%(self.n+1,self.tn+dt))
908
909 def run(self,check_point=None):
910 """
911 Run the simulation by performing essentially::
912
913 self.doInitialization()
914 while not self.finalize():
915 dt=self.getSafeTimeStepSize()
916 self.doStep(dt)
917 if n%check_point==0:
918 self.writeXML()
919 self.doFinalization()
920
921 If one of the models in throws a C{FailedTimeStepError} exception a
922 new time step size is computed through getSafeTimeStepSize() and the
923 time step is repeated.
924
925 If one of the models in throws a C{IterationDivergenceError}
926 exception the time step size is halved and the time step is repeated.
927
928 In both cases the time integration is given up after
929 C{Simulation.FAILED_TIME_STEPS_MAX} attempts.
930 """
931 dt=self.UNDEF_DT
932 self.doInitialization()
933 while not self.finalize():
934 step_fail_counter=0
935 iteration_fail_counter=0
936 if self.n==0:
937 dt_new=self.getSafeTimeStepSize(dt)
938 else:
939 dt_new=min(max(self.getSafeTimeStepSize(dt),dt/self.MAX_CHANGE_OF_DT),dt*self.MAX_CHANGE_OF_DT)
940 self.trace("%d. time step %e (step size %e.)" % (self.n+1,self.tn+dt_new,dt_new))
941 end_of_step=False
942 while not end_of_step:
943 end_of_step=True
944 if not dt_new>0:
945 raise NonPositiveStepSizeError("non-positive step size in step %d"%(self.n+1))
946 try:
947 self.doStepPreprocessing(dt_new)
948 self.doStep(dt_new)
949 self.doStepPostprocessing(dt_new)
950 except IterationDivergenceError:
951 dt_new*=0.5
952 end_of_step=False
953 iteration_fail_counter+=1
954 if iteration_fail_counter>self.FAILED_TIME_STEPS_MAX:
955 raise SimulationBreakDownError("reduction of time step to achieve convergence failed after %s steps."%self.FAILED_TIME_STEPS_MAX)
956 self.trace("Iteration failed. Time step is repeated with new step size %s."%dt_new)
957 except FailedTimeStepError:
958 dt_new=self.getSafeTimeStepSize(dt)
959 end_of_step=False
960 step_fail_counter+=1
961 self.trace("Time step is repeated with new time step size %s."%dt_new)
962 if step_fail_counter>self.FAILED_TIME_STEPS_MAX:
963 raise SimulationBreakDownError("Time integration is given up after %d attempts."%step_fail_counter)
964 dt=dt_new
965 if not check_point==None:
966 if n%check_point==0:
967 self.trace("check point is created.")
968 self.writeXML()
969 self.doFinalization()
970
971 def fromDom(cls, doc):
972 sims = []
973 for node in doc.childNodes:
974 if isinstance(node, minidom.Text):
975 continue
976
977 sims.append(getComponent(node))
978
979 return cls(sims)
980
981 fromDom = classmethod(fromDom)
982
983
984 class IterationDivergenceError(Exception):
985 """
986 Exception which is thrown if there is no convergence of the iteration
987 process at a time step.
988
989 But there is a chance that a smaller step could help to reach convergence.
990 """
991 pass
992
993 class FailedTimeStepError(Exception):
994 """
995 Exception which is thrown if the time step fails because of a step
996 size that have been choosen to be too large.
997 """
998 pass
999
1000 class SimulationBreakDownError(Exception):
1001 """
1002 Exception which is thrown if the simulation does not manage to
1003 progress in time.
1004 """
1005 pass
1006
1007 class NonPositiveStepSizeError(Exception):
1008 """
1009 Exception which is thrown if the step size is not positive.
1010 """
1011 pass
1012
1013 # 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