/[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 871 - (show annotations)
Sat Oct 14 08:25:54 2006 UTC (12 years, 11 months ago) by elspeth
File MIME type: text/x-python
File size: 31201 byte(s)
Numarray support added to modelframe and tests for numarray added to run_xml.
Numarrays will look like this in ESysXML:

A 2 x 3 numarray will look like this:

<Parameter type="NumArray">
                        <Name>
                                numtest
                        </Name>
                        <Value>
                                <ArrayType>
                                        Float64
                                </ArrayType>
                                <Shape>
                                        2 3
                                </Shape>
                                <Data>
                                        1.0 2.0 3.0 3.0 4.0 5.0
                                </Data>
                        </Value>
                </Parameter>

A numarray vector will look like this:

<Parameter type="NumArray">
                        <Name>
                                numtest
                        </Name>
                        <Value>
                                <ArrayType>
                                        Float64
                                </ArrayType>
                                <Shape>
                                        3
                                </Shape>
                                <Data>
                                        3.0 4.0 5.0
                                </Data>
                        </Value>
                </Parameter>



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 - any other object (not considered by writeESySXML and writeXML)
370
371 Example how to create an ESySParameters object::
372
373 p11=ParameterSet(gamma1=1.,gamma2=2.,gamma3=3.)
374 p1=ParameterSet(dim=2,tol_v=0.001,output_file="/tmp/u.%3.3d.dx",runFlag=True,parm11=p11)
375 parm=ParameterSet(parm1=p1,parm2=ParameterSet(alpha=Link(p11,"gamma1")))
376
377 This can be accessed as::
378
379 parm.parm1.gamma=0.
380 parm.parm1.dim=2
381 parm.parm1.tol_v=0.001
382 parm.parm1.output_file="/tmp/u.%3.3d.dx"
383 parm.parm1.runFlag=True
384 parm.parm1.parm11.gamma1=1.
385 parm.parm1.parm11.gamma2=2.
386 parm.parm1.parm11.gamma3=3.
387 parm.parm2.alpha=1. (value of parm.parm1.parm11.gamma1)
388 """
389 def __init__(self, parameters=[], **kwargs):
390 """
391 Creates a ParameterSet with parameters parameters.
392 """
393 LinkableObject.__init__(self, **kwargs)
394 self.parameters = set()
395 self.declareParameters(parameters)
396
397 def __repr__(self):
398 return "<%s %r>" % (self.__class__.__name__,
399 [(p, getattr(self, p, None)) for p in self.parameters])
400
401 def declareParameter(self,**parameters):
402 """
403 Declares a new parameter(s) and its (their) initial value.
404 """
405 self.declareParameters(parameters)
406
407 def declareParameters(self,parameters):
408 """
409 Declares a set of parameters. parameters can be a list, a dictionary
410 or a ParameterSet.
411 """
412 if isinstance(parameters,ListType):
413 parameters = zip(parameters, itertools.repeat(None))
414 if isinstance(parameters,DictType):
415 parameters = parameters.iteritems()
416
417 for prm, value in parameters:
418 setattr(self,prm,value)
419 self.parameters.add(prm)
420
421 self.trace("parameter %s has been declared."%prm)
422
423 def releaseParameters(self,name):
424 """
425 Removes parameter name from the paramameters.
426 """
427 if self.isParameter(name):
428 self.parameters.remove(name)
429 self.trace("parameter %s has been removed."%name)
430
431 def __iter__(self):
432 """
433 Creates an iterator over the parameter and their values.
434 """
435 return _ParameterIterator(self)
436
437 def showParameters(self):
438 """
439 Returns a descrition of the parameters.
440 """
441 out="{"
442 notfirst=False
443 for i,v in self:
444 if notfirst: out=out+","
445 notfirst=True
446 if isinstance(v,ParameterSet):
447 out="%s\"%s\" : %s"%(out,i,v.showParameters())
448 else:
449 out="%s\"%s\" : %s"%(out,i,v)
450 return out+"}"
451
452 def __delattr__(self,name):
453 """
454 Removes the attribute name.
455 """
456 LinkableObject.__delattr__(self,name)
457 try:
458 self.releaseParameter(name)
459 except:
460 pass
461
462 def toDom(self, document, node):
463 """
464 C{toDom} method of ParameterSet class.
465 """
466 pset = document.createElement('ParameterSet')
467 node.appendChild(pset)
468 self._parametersToDom(document, pset)
469
470 def _parametersToDom(self, document, node):
471 node.setAttribute('id', str(self.id))
472 node.setIdAttribute("id")
473 for name,value in self:
474 param = document.createElement('Parameter')
475 param.setAttribute('type', value.__class__.__name__)
476
477 param.appendChild(dataNode(document, 'Name', name))
478
479 val = document.createElement('Value')
480
481 if isinstance(value,ParameterSet):
482 value.toDom(document, val)
483 param.appendChild(val)
484 elif isinstance(value, Link):
485 value.toDom(document, val)
486 param.appendChild(val)
487 elif isinstance(value, numarray.NumArray):
488 shape = value.getshape()
489 if isinstance(shape, tuple):
490 size = reduce(operator.mul, shape)
491 shape = ' '.join(map(str, shape))
492 else:
493 size = shape
494 shape = str(shape)
495
496 arraytype = value.type()
497 val.appendChild(dataNode(document, 'ArrayType', str(arraytype)))
498 val.appendChild(dataNode(document, 'Shape', shape))
499 val.appendChild(dataNode(document, 'Data', ' '.join(
500 [str(x) for x in numarray.reshape(value, size)])))
501 param.appendChild(val)
502 elif isinstance(value,StringType):
503 param.appendChild(dataNode(document, 'Value', value))
504 else:
505 param.appendChild(dataNode(document, 'Value', str(value)))
506
507 node.appendChild(param)
508
509 def fromDom(cls, doc):
510
511 # Define a host of helper functions to assist us.
512 def _children(node):
513 """
514 Remove the empty nodes from the children of this node.
515 """
516 ret = []
517 for x in node.childNodes:
518 if isinstance(x, minidom.Text):
519 if x.nodeValue.strip():
520 ret.append(x)
521 else:
522 ret.append(x)
523 return ret
524
525 def _floatfromValue(doc):
526 return float(doc.nodeValue.strip())
527
528 def _stringfromValue(doc):
529 return str(doc.nodeValue.strip())
530
531 def _intfromValue(doc):
532 return int(doc.nodeValue.strip())
533
534 def _boolfromValue(doc):
535 return bool(doc.nodeValue.strip())
536
537 def _nonefromValue(doc):
538 return None
539
540 def _numarrayfromValue(doc):
541 for node in doc:
542 if node.tagName == 'ArrayType':
543 arraytype = node.firstChild.nodeValue.strip()
544 if node.tagName == 'Shape':
545 shape = node.firstChild.nodeValue.strip()
546 shape = [int(x) for x in shape.split()]
547 if node.tagName == 'Data':
548 data = node.firstChild.nodeValue.strip()
549 data = [float(x) for x in data.split()]
550 return numarray.reshape(numarray.array(data, type=getattr(numarray, arraytype)),
551 shape)
552
553 # Mapping from text types in the xml to methods used to process trees of that type
554 ptypemap = {"Simulation": Simulation.fromDom,
555 "Model":Model.fromDom,
556 "ParameterSet":ParameterSet.fromDom,
557 "Link":Link.fromDom,
558 "float":_floatfromValue,
559 "int":_intfromValue,
560 "str":_stringfromValue,
561 "bool":_boolfromValue,
562 "NoneType":_nonefromValue,
563 }
564
565 # print doc.toxml()
566
567 parameters = {}
568 for node in _children(doc):
569 ptype = node.getAttribute("type")
570
571 pname = pvalue = None
572 for childnode in _children(node):
573
574 if childnode.tagName == "Name":
575 pname = childnode.firstChild.nodeValue.strip()
576
577 if childnode.tagName == "Value":
578 nodes = _children(childnode)
579 if ptype == 'NumArray':
580 pvalue = _numarrayfromValue(nodes)
581 else:
582 pvalue = ptypemap[ptype](nodes[0])
583
584 parameters[pname] = pvalue
585
586 # Create the instance of ParameterSet
587 o = cls()
588 o.declareParameters(parameters)
589 registerLinkableObject(doc.getAttribute("id"), o)
590 return o
591
592 fromDom = classmethod(fromDom)
593
594 def writeXML(self,ostream=stdout):
595 """
596 Writes the object as an XML object into an output stream.
597 """
598 # ParameterSet(d) with d[Name]=Value
599 document, node = esysDoc()
600 self.toDom(document, node)
601 ostream.write(document.toprettyxml())
602
603 class Model(ParameterSet):
604 """
605 A Model object represents a processess marching over time until a
606 finalizing condition is fullfilled. At each time step an iterative
607 process can be performed and the time step size can be controlled. A
608 Model has the following work flow::
609
610 doInitialization()
611 while not finalize():
612 dt=getSafeTimeStepSize(dt)
613 doStepPreprocessing(dt)
614 while not terminateIteration(): doStep(dt)
615 doStepPostprocessing(dt)
616 doFinalization()
617
618 where C{doInitialization}, C{finalize}, C{getSafeTimeStepSize},
619 C{doStepPreprocessing}, C{terminateIteration}, C{doStepPostprocessing},
620 C{doFinalization} are methods of the particular instance of a Model. The
621 default implementations of these methods have to be overwritten by the
622 subclass implementing a Model.
623 """
624
625 UNDEF_DT=1.e300
626
627 def __init__(self,parameters=[],**kwarg):
628 """
629 Creates a model.
630
631 Just calls the parent constructor.
632 """
633 ParameterSet.__init__(self, parameters=parameters,**kwarg)
634
635 def __str__(self):
636 return "<%s %d>"%(self.__class__,id(self))
637
638 def toDom(self, document, node):
639 """
640 C{toDom} method of Model class
641 """
642 pset = document.createElement('Model')
643 pset.setAttribute('type', self.__class__.__name__)
644 if not self.__class__.__module__.startswith('esys.escript'):
645 pset.setAttribute('module', self.__class__.__module__)
646 node.appendChild(pset)
647 self._parametersToDom(document, pset)
648
649 def doInitialization(self):
650 """
651 Initializes the time stepping scheme.
652
653 This function may be overwritten.
654 """
655 pass
656
657 def getSafeTimeStepSize(self,dt):
658 """
659 Returns a time step size which can safely be used.
660
661 C{dt} gives the previously used step size.
662
663 This function may be overwritten.
664 """
665 return self.UNDEF_DT
666
667 def finalize(self):
668 """
669 Returns False if the time stepping is finalized.
670
671 This function may be overwritten.
672 """
673 return False
674
675 def doFinalization(self):
676 """
677 Finalizes the time stepping.
678
679 This function may be overwritten.
680 """
681 pass
682
683 def doStepPreprocessing(self,dt):
684 """
685 Sets up a time step of step size dt.
686
687 This function may be overwritten.
688 """
689 pass
690
691 def doStep(self,dt):
692 """
693 Executes an iteration step at a time step.
694
695 C{dt} is the currently used time step size.
696
697 This function may be overwritten.
698 """
699 pass
700
701 def terminateIteration(self):
702 """
703 Returns True if iteration on a time step is terminated.
704 """
705 return True
706
707 def doStepPostprocessing(self,dt):
708 """
709 Finalalizes the time step.
710
711 dt is the currently used time step size.
712
713 This function may be overwritten.
714 """
715 pass
716
717 def writeXML(self, ostream=stdout):
718 document, node = esysDoc()
719 self.toDom(document, node)
720 ostream.write(document.toprettyxml())
721
722
723 class Simulation(Model):
724 """
725 A Simulation object is special Model which runs a sequence of Models.
726
727 The methods C{doInitialization}, C{finalize}, C{getSafeTimeStepSize},
728 C{doStepPreprocessing}, C{terminateIteration}, C{doStepPostprocessing},
729 C{doFinalization} are executing the corresponding methods of the models in
730 the simulation.
731 """
732
733 FAILED_TIME_STEPS_MAX=20
734 MAX_ITER_STEPS=50
735 MAX_CHANGE_OF_DT=2.
736
737 def __init__(self, models=[], **kwargs):
738 """
739 Initiates a simulation from a list of models.
740 """
741 Model.__init__(self, **kwargs)
742 self.__models=[]
743
744 for i in range(len(models)):
745 self[i] = models[i]
746
747
748 def __repr__(self):
749 """
750 Returns a string representation of the Simulation.
751 """
752 return "<Simulation %r>" % self.__models
753
754 def __str__(self):
755 """
756 Returning Simulation as a string.
757 """
758 return "<Simulation %d>"%id(self)
759
760 def iterModels(self):
761 """
762 Returns an iterator over the models.
763 """
764 return self.__models
765
766 def __getitem__(self,i):
767 """
768 Returns the i-th model.
769 """
770 return self.__models[i]
771
772 def __setitem__(self,i,value):
773 """
774 Sets the i-th model.
775 """
776 if not isinstance(value,Model):
777 raise ValueError,"assigned value is not a Model but instance of %s"%(value.__class__.__name__,)
778 for j in range(max(i-len(self.__models)+1,0)):
779 self.__models.append(None)
780 self.__models[i]=value
781
782 def __len__(self):
783 """
784 Returns the number of models.
785 """
786 return len(self.__models)
787
788 def toDom(self, document, node):
789 """
790 C{toDom} method of Simulation class.
791 """
792 simulation = document.createElement('Simulation')
793 simulation.setAttribute('type', self.__class__.__name__)
794
795 for rank, sim in enumerate(self.iterModels()):
796 component = document.createElement('Component')
797 component.setAttribute('rank', str(rank))
798
799 sim.toDom(document, component)
800
801 simulation.appendChild(component)
802
803 node.appendChild(simulation)
804
805 def writeXML(self,ostream=stdout):
806 """
807 Writes the object as an XML object into an output stream.
808 """
809 document, rootnode = esysDoc()
810 self.toDom(document, rootnode)
811 targetsList = document.getElementsByTagName('Target')
812
813 for element in targetsList:
814 targetId = int(element.firstChild.nodeValue.strip())
815 if document.getElementById(str(targetId)):
816 continue
817 targetObj = LinkableObjectRegistry[targetId]
818 targetObj.toDom(document, rootnode)
819 ostream.write(document.toprettyxml())
820
821 def getSafeTimeStepSize(self,dt):
822 """
823 Returns a time step size which can safely be used by all models.
824
825 This is the minimum over the time step sizes of all models.
826 """
827 out=min([o.getSafeTimeStepSize(dt) for o in self.iterModels()])
828 #print "%s: safe step size is %e."%(str(self),out)
829 return out
830
831 def doInitialization(self):
832 """
833 Initializes all models.
834 """
835 self.n=0
836 self.tn=0.
837 for o in self.iterModels():
838 o.doInitialization()
839
840 def finalize(self):
841 """
842 Returns True if any of the models is to be finalized.
843 """
844 return any([o.finalize() for o in self.iterModels()])
845
846 def doFinalization(self):
847 """
848 Finalalizes the time stepping for all models.
849 """
850 for i in self.iterModels(): i.doFinalization()
851 self.trace("end of time integation.")
852
853 def doStepPreprocessing(self,dt):
854 """
855 Initializes the time step for all models.
856 """
857 for o in self.iterModels():
858 o.doStepPreprocessing(dt)
859
860 def terminateIteration(self):
861 """
862 Returns True if all iterations for all models are terminated.
863 """
864 out=all([o.terminateIteration() for o in self.iterModels()])
865 return out
866
867 def doStepPostprocessing(self,dt):
868 """
869 Finalalizes the iteration process for all models.
870 """
871 for o in self.iterModels():
872 o.doStepPostprocessing(dt)
873 self.n+=1
874 self.tn+=dt
875
876 def doStep(self,dt):
877 """
878 Executes the iteration step at a time step for all model::
879
880 self.doStepPreprocessing(dt)
881 while not self.terminateIteration():
882 for all models:
883 self.doStep(dt)
884 self.doStepPostprocessing(dt)
885 """
886 self.iter=0
887 while not self.terminateIteration():
888 if self.iter==0: self.trace("iteration at %d-th time step %e starts"%(self.n+1,self.tn+dt))
889 self.iter+=1
890 self.trace("iteration step %d"%(self.iter))
891 for o in self.iterModels():
892 o.doStep(dt)
893 if self.iter>0: self.trace("iteration at %d-th time step %e finalized."%(self.n+1,self.tn+dt))
894
895 def run(self,check_point=None):
896 """
897 Run the simulation by performing essentially::
898
899 self.doInitialization()
900 while not self.finalize():
901 dt=self.getSafeTimeStepSize()
902 self.doStep(dt)
903 if n%check_point==0:
904 self.writeXML()
905 self.doFinalization()
906
907 If one of the models in throws a C{FailedTimeStepError} exception a
908 new time step size is computed through getSafeTimeStepSize() and the
909 time step is repeated.
910
911 If one of the models in throws a C{IterationDivergenceError}
912 exception the time step size is halved and the time step is repeated.
913
914 In both cases the time integration is given up after
915 C{Simulation.FAILED_TIME_STEPS_MAX} attempts.
916 """
917 dt=self.UNDEF_DT
918 self.doInitialization()
919 while not self.finalize():
920 step_fail_counter=0
921 iteration_fail_counter=0
922 if self.n==0:
923 dt_new=self.getSafeTimeStepSize(dt)
924 else:
925 dt_new=min(max(self.getSafeTimeStepSize(dt),dt/self.MAX_CHANGE_OF_DT),dt*self.MAX_CHANGE_OF_DT)
926 self.trace("%d. time step %e (step size %e.)" % (self.n+1,self.tn+dt_new,dt_new))
927 end_of_step=False
928 while not end_of_step:
929 end_of_step=True
930 if not dt_new>0:
931 raise NonPositiveStepSizeError("non-positive step size in step %d"%(self.n+1))
932 try:
933 self.doStepPreprocessing(dt_new)
934 self.doStep(dt_new)
935 self.doStepPostprocessing(dt_new)
936 except IterationDivergenceError:
937 dt_new*=0.5
938 end_of_step=False
939 iteration_fail_counter+=1
940 if iteration_fail_counter>self.FAILED_TIME_STEPS_MAX:
941 raise SimulationBreakDownError("reduction of time step to achieve convergence failed after %s steps."%self.FAILED_TIME_STEPS_MAX)
942 self.trace("Iteration failed. Time step is repeated with new step size %s."%dt_new)
943 except FailedTimeStepError:
944 dt_new=self.getSafeTimeStepSize(dt)
945 end_of_step=False
946 step_fail_counter+=1
947 self.trace("Time step is repeated with new time step size %s."%dt_new)
948 if step_fail_counter>self.FAILED_TIME_STEPS_MAX:
949 raise SimulationBreakDownError("Time integration is given up after %d attempts."%step_fail_counter)
950 dt=dt_new
951 if not check_point==None:
952 if n%check_point==0:
953 self.trace("check point is created.")
954 self.writeXML()
955 self.doFinalization()
956
957 def fromDom(cls, doc):
958 sims = []
959 for node in doc.childNodes:
960 if isinstance(node, minidom.Text):
961 continue
962
963 sims.append(getComponent(node))
964
965 return cls(sims)
966
967 fromDom = classmethod(fromDom)
968
969
970 class IterationDivergenceError(Exception):
971 """
972 Exception which is thrown if there is no convergence of the iteration
973 process at a time step.
974
975 But there is a chance that a smaller step could help to reach convergence.
976 """
977 pass
978
979 class FailedTimeStepError(Exception):
980 """
981 Exception which is thrown if the time step fails because of a step
982 size that have been choosen to be too large.
983 """
984 pass
985
986 class SimulationBreakDownError(Exception):
987 """
988 Exception which is thrown if the simulation does not manage to
989 progress in time.
990 """
991 pass
992
993 class NonPositiveStepSizeError(Exception):
994 """
995 Exception which is thrown if the step size is not positive.
996 """
997 pass
998
999 # 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