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