/[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 123 - (show annotations)
Fri Jul 8 04:08:13 2005 UTC (13 years, 9 months ago) by jgs
Original Path: trunk/esys2/escript/py_src/modelframe.py
File MIME type: text/x-python
File size: 37492 byte(s)
Merge of development branch back to main trunk on 2005-07-08

1 # $Id$
2
3 from types import StringType,IntType,FloatType,BooleanType,ListType,DictType
4 from sys import stdout
5 import itertools
6
7 import modellib
8
9 # import the 'set' module if it's not defined (python2.3/2.4 difference)
10 try:
11 set
12 except NameError:
13 from sets import Set as set
14
15 from xml.dom import minidom
16
17 def dataNode(document, tagName, data):
18 t = document.createTextNode(str(data))
19 n = document.createElement(tagName)
20 n.appendChild(t)
21 return n
22
23 def esysDoc():
24 doc = minidom.Document()
25 esys = doc.createElement('ESys')
26 doc.appendChild(esys)
27 return doc, esys
28
29 def all(seq):
30 for x in seq:
31 if not x:
32 return False
33 return True
34
35 LinkableObjectRegistry = {}
36
37 def registerLinkableObject(obj_id, o):
38 LinkableObjectRegistry[obj_id] = o
39
40 LinkRegistry = []
41
42 def registerLink(obj_id, l):
43 LinkRegistry.append((obj_id,l))
44
45 def parse(xml):
46 global LinkRegistry, LinkableObjectRegistry
47 LinkRegistry = []
48 LinkableObjectRegistry = {}
49
50 doc = minidom.parseString(xml)
51
52 sim = getComponent(doc.firstChild)
53 for obj_id, link in LinkRegistry:
54 link.target = LinkableObjectRegistry[obj_id]
55
56 return sim
57
58 def getComponent(doc):
59 for node in doc.childNodes:
60 if isinstance(node, minidom.Element):
61 if node.tagName == 'Simulation':
62 if node.getAttribute("type") == 'Simulation':
63 return Simulation.fromDom(node)
64 elif node.getAttribute("type") == 'ExplicitSimulation':
65 return ExplicitSimulation.fromDom(node)
66 if node.tagName == 'Model':
67 model_type = node.getAttribute("type")
68 model_subclasses = Model.__subclasses__()
69 for model in model_subclasses:
70 if model_type == model.__name__:
71 return model.fromDom(node)
72
73 raise "Invalid simulation type, %r" % node.getAttribute("type")
74
75 raise ValueError("No Simulation Found")
76
77
78 class Link:
79 """
80 a Link makes an attribute of an object callable:
81 o.object()
82 o.a=8
83 l=Link(o,"a")
84 assert l()==8
85 """
86
87 def __init__(self,target,attribute=None):
88 """creates a link to the object target. If attribute is given, the link is
89 establised to this attribute of the target. Otherwise the attribute is
90 undefined."""
91 self.target = target
92 self.attribute = None
93 self.setAttributeName(attribute)
94
95 def setAttributeName(self,attribute):
96 """set a new attribute name to be collected from the target object. The
97 target object must have the attribute with name attribute."""
98 if attribute and self.target and not hasattr(self.target, attribute):
99 raise AttributeError("%s: target %s has no attribute %s."%(self, self.target, attribute))
100 self.attribute = attribute
101
102 def hasDefinedAttributeName(self):
103 """returns true if an attribute name is set"""
104 return self.attribute != None
105
106 def __repr__(self):
107 """returns a string representation of the link"""
108 if self.hasDefinedAttributeName():
109 return "<Link to attribute %s of %s>" % (self.attribute, self.target)
110 else:
111 return "<Link to target %s>" % self.target
112
113 def __call__(self,name=None):
114 """returns the value of the attribute of the target object. If the
115 atrribute is callable then the return value of the call is returned."""
116 if name:
117 out=getattr(self.target, name)
118 else:
119 out=getattr(self.target, self.attribute)
120
121 if callable(out):
122 return out()
123 else:
124 return out
125
126 def toDom(self, document, node):
127 """ toDom method of Link. Creates a Link node and appends it to the current XML
128 document """
129 link = document.createElement('Link')
130 link.appendChild(dataNode(document, 'Target', self.target.id))
131 # this use of id will not work for purposes of being able to retrieve the intended
132 # target from the xml later. I need a better unique identifier.
133 assert self.attribute, "You can't xmlify a Link without a target attribute"
134 link.appendChild(dataNode(document, 'Attribute', self.attribute))
135 node.appendChild(link)
136
137 def fromDom(cls, doc):
138 targetid = doc.getElementsByTagName("Target")[0].firstChild.nodeValue.strip()
139 attribute = doc.getElementsByTagName("Attribute")[0].firstChild.nodeValue.strip()
140 l = cls(None, attribute)
141 registerLink(targetid, l)
142 return l
143
144 fromDom = classmethod(fromDom)
145
146 def writeXML(self,ostream=stdout):
147 """writes an XML representation of self to the output stream ostream.
148 If ostream is nor present the standart output stream is used. If
149 esysheader==True the esys XML header is written"""
150
151 document, rootnode = esysDoc()
152 self.toDom(document, rootnode)
153
154 ostream.write(document.toprettyxml())
155
156 class LinkableObject(object):
157 """
158 An object that allows to link its attributes to attributes of other objects
159 via a Link object. For instance
160
161 p = LinkableObject()
162 p.x = Link(o,"name")
163 print p.x
164
165 links attribute x of p to the attribute name of object o.
166
167 p.x will contain the current value of attribute name of object o.
168
169 If the value of getattr(o, "name") is callable, p.x will rturn the return
170 value of the call.
171 """
172
173 number_sequence = itertools.count(100)
174
175 def __init__(self, debug=False):
176 """ initializes LinkableObject so that we can operate on Links """
177 self.debug = debug
178 self.__linked_attributes={}
179 self.id = self.number_sequence.next()
180
181 def trace(self, msg):
182 """ If debugging is on, print the message, otherwise do nothing
183 """
184 if self.debug:
185 print msg
186
187 def __getattr__(self,name):
188 """returns the value of attribute name. If the value is a Link object the
189 object is called and the return value is returned."""
190 out = self.getAttributeObject(name)
191 if isinstance(out,Link):
192 return out()
193 else:
194 return out
195
196 def getAttributeObject(self,name):
197 """return the object stored for attribute name."""
198
199 if self.__dict__.has_key(name):
200 return self.__dict__[name]
201
202 if self.__linked_attributes.has_key(name):
203 return self.__linked_attributes[name]
204
205 raise AttributeError,"No attribute %s."%name
206
207 def __setattr__(self,name,value):
208 """sets the value for attribute name. If value is a Link the target
209 attribute is set to name if no attribute has been specified."""
210
211
212 if self.__dict__.has_key(name):
213 del self.__dict__[name]
214
215 if isinstance(value,Link):
216 if not value.hasDefinedAttributeName():
217 value.setAttributeName(name)
218 self.__linked_attributes[name] = value
219
220 self.trace("DEBUG: %s: attribute %s is now linked by %s."%(self,name,value))
221 else:
222 self.__dict__[name] = value
223
224 def __delattr__(self,name):
225 """removes the attribute name."""
226
227 if self.__linked_attributes.has_key[name]:
228 del self.__linked_attributes[name]
229 elif self.__dict__.has_key(name):
230 del self.__dict__[name]
231 else:
232 raise AttributeError,"No attribute %s."%name
233
234 class SimulationFrame(LinkableObject):
235 """A SimulationFrame objects represents a processess marching over time
236 until a finalizing condition is fullfilled. At each time step an iterative
237 process can be performed and the time step size can be controlled
238 """
239 UNDEF_DT=1.e300
240 MAX_TIME_STEP_REDUCTION=20
241 MAX_ITER_STEPS=50
242
243 def __init__(self,*args, **kwargs):
244 """
245 Initialises a simulation
246
247 Just calls the parent constructor.
248 """
249 LinkableObject.__init__(self,*args, **kwargs)
250
251 def doInitialization(self,t):
252 """initializes the time stepping scheme. This function may be
253 overwritten."""
254 pass
255
256 def getSafeTimeStepSize(self,dt):
257 """returns a time step size which can savely be used. This function may
258 be overwritten."""
259 return self.UNDEF_DT
260
261 def finalize(self):
262 """returns True if the time stepping is finalized. This function may be
263 overwritten."""
264 return True
265
266 def doFinalization(self):
267 """finalizes the time stepping. This function may be overwritten."""
268 pass
269
270 def doIterationInitialization(self,dt):
271 """initializes the iteration at time step t. This function may be
272 overwritten. (only called if doStep is not overwritten)"""
273 pass
274
275 def doIterationStep(self,dt):
276 """executes the iteration step. This function may be overwritten. (only
277 called if doStep is not overwritten)"""
278 pass
279
280 def terminate(self):
281 """returns True if iteration on a time step is terminated. (only called
282 if doStep is not overwritten)"""
283 return True
284
285 def doIterationFinalization(self,dt):
286 """finalalizes the iteration process. (only called if doStep is not
287 overwritten)"""
288 pass
289
290 def run(self,check_point=None):
291 """run the simulation by performing essentially
292
293 self.doInitialization()
294 while not self.finalize():
295 dt=self.getSafeTimeStepSize()
296 self.doStep(dt)
297 if n%check_point==0: self.writeXML()
298 self.doFinalization()
299
300 """
301 self.__tn=0.
302 self.__n=0
303 self.__dt=None
304 self.doInitialization(self.__tn)
305 while not self.finalize():
306 self.__n+=1
307 self.__dt=self.getSafeTimeStepSize(self.__dt)
308 if self.__dt==None: self.__dt=self.UNDEF_DT
309 if not self.__dt>0:
310 raise NonPositiveStepSizeError("non-positive step size in step %d",self.__n)
311 self.trace("%s: %d. time step %e (step size %e.)" %
312 (self,self.__n,self.__tn+self.__dt,self.__dt))
313 endoftimestep=False
314 failcounter=0
315 while not endoftimestep:
316 endoftimestep=True
317 try:
318 self.doStep(self.__dt)
319 except FailedTimeStepError:
320 self.__dt=self.getSafeTimeStepSize(self.__dt)
321 if self.__dt==None: self.__dt=self.UNDEF_DT
322 endoftimestep=False
323 self.trace("%s: time step is repeated with new step size %e."%(self,self.__dt))
324 except IterationDivergenceError:
325 self.__dt*=0.5
326 endoftimestep=False
327 failcounter+=1
328 if failcounter>self.MAX_TIME_STEP_REDUCTION:
329 raise IterationBreakDownError("reduction of time step to achieve convergence failed.")
330
331 self.trace("%s: iteration failes. time step is repeated with new step size %e."
332 % (self,self.__dt))
333 self.__tn+=self.__dt
334 if not check_point==None:
335 if self.__n%check_point==0:
336 self.trace("%s: check point is created."%self)
337 self.writeXML()
338 self.doFinalization()
339
340 def writeXML(self):
341 raise RuntimeError, "Not implemented"
342
343 def doStep(self,dt):
344 """executes a time step by iteration. This function may be overwritten.
345
346 basicly it performs :
347
348 self.doIterationInitialization(dt)
349 while not self.terminate(): self.doIterationStep(dt)
350 self.doIterationFinalization(dt)
351
352 """
353 self.doIterationInitialization(dt)
354 self.__iter=0
355 while not self.terminate():
356 self.trace("%s: iteration step %d"%(self,self.__iter))
357 self.doIterationStep(dt)
358 self.__iter+=1
359 if self.__iter>self.MAX_ITER_STEPS:
360 raise IterationDivergenceError("%s: divergence in step %d"%(self,self.__iter))
361 self.doIterationFinalization(dt)
362
363 class Simulation(SimulationFrame):
364 """A Simulation object is comprised by SimulationFrame(s) called subsimulations."""
365
366 def __init__(self, subsimulations=[], *args, **kwargs):
367 """initiates a simulation from a list of subsimulations. """
368 SimulationFrame.__init__(self, *args, **kwargs)
369 self.__subsimulations=[]
370
371 for i in range(len(subsimulations)):
372 self[i] = subsimulations[i]
373
374 def iterSubsimulations(self):
375 """returns an iterator over the subsimulations"""
376 return self.__subsimulations
377
378 def __getitem__(self,i):
379 """returns the i-th subsimulation"""
380 return self.__subsimulations[i]
381
382 def __setitem__(self,i,value):
383 """sets the i-th subsimulation"""
384 if not isinstance(value,SimulationFrame):
385 raise ValueError("assigned value is not a Simulation")
386 for j in range(max(i-len(self.__subsimulations)+1,0)): self.__subsimulations.append(None)
387 self.__subsimulations[i]=value
388
389 def __len__(self):
390 """returns the number of subsimulations"""
391 return len(self.__subsimulations)
392
393 def toDom(self, document, node):
394 """ toDom method of Simulation class """
395 simulation = document.createElement('Simulation')
396 simulation.setAttribute('type', self.__class__.__name__)
397
398 for rank, sim in enumerate(self.iterSubsimulations()):
399 component = document.createElement('Component')
400 component.setAttribute('rank', str(rank))
401
402 sim.toDom(document, component)
403
404 simulation.appendChild(component)
405
406 node.appendChild(simulation)
407
408 def writeXML(self,ostream=stdout):
409 """writes the object as an XML object into an output stream"""
410 document, rootnode = esysDoc()
411 self.toDom(document, rootnode)
412 ostream.write(document.toprettyxml())
413
414 def getSafeTimeStepSize(self,dt):
415 """returns a time step size which can safely be used by all subsimulations"""
416 out=self.UNDEF_DT
417 for o in self.iterSubsimulations():
418 dt_new = o.getSafeTimeStepSize(dt)
419 if dt_new != None:
420 out = min(out,dt_new)
421
422 return out
423
424 def doInitialization(self,dt):
425 """initializes all subsimulations """
426 for o in self.iterSubsimulations():
427 o.doInitialization(dt)
428
429 def finalize(self):
430 """returns True if all subsimulations are finalized"""
431 return all([o.finalize() for o in self.iterSubsimulations()])
432
433 def doFinalization(self):
434 """finalalizes the time stepping for all subsimulations."""
435 for i in self.iterSubsimulations(): i.doFinalization()
436
437 def doIterationInitialization(self,dt):
438 """initializes the iteration at time t for all subsimulations."""
439 self.__iter=0
440 self.trace("%s: iteration starts"%self)
441
442 for o in self.iterSubsimulations():
443 o.doIterationInitialization(dt)
444
445 def terminate(self):
446 """returns True if all iterations for all subsimulations are terminated."""
447 return all([o.terminate() for o in self.iterSubsimulations()])
448
449 def doIterationFinalization(self,dt):
450 """finalalizes the iteration process for each of the subsimulations."""
451 for o in self.iterSubsimulations():
452 o.doIterationFinalization(dt)
453 self.trace("%s: iteration finalized after %s steps"%(self,self.__iter+1))
454
455 def doIterationStep(self,dt):
456 """executes the iteration step at time step for each subsimulation"""
457 self.__iter+=1
458 self.trace("%s: iteration step %d"%(self,self.__iter))
459 for o in self.iterSubsimulations():
460 o.doIterationStep(dt)
461
462 def fromDom(cls, doc):
463 """
464 Needs to be filled out.
465 """
466 sims = []
467 for node in doc.childNodes:
468 if isinstance(node, minidom.Text):
469 continue
470
471 sims.append(getComponent(node))
472
473
474 return cls(sims)
475
476
477
478
479
480 fromDom = classmethod(fromDom)
481
482 class ExplicitSimulation(Simulation):
483 """This is a modified form of the Simulation class. In fact it overwrites
484 the doStep method by executing the doStep method of all subsimulations
485 rather then iterating over all subsimulations."""
486
487 def doStep(self,dt):
488 """executes the time step for all subsimulation"""
489 for i in self.iterSubsimulations():
490 i.doStep(dt)
491
492 class _ParameterIterator:
493 def __init__(self,parameterset):
494
495 self.__set=parameterset
496 self.__iter=iter(parameterset.parameters)
497
498 def next(self):
499 o=self.__iter.next()
500 return (o,self.__set.getAttributeObject(o))
501
502 def __iter__(self):
503 return self
504
505 class ParameterSet(LinkableObject):
506 """a class which allows to emphazise attributes to be written and read to XML
507
508 Leaves of an ESySParameters objects can be
509
510 a real number
511 a integer number
512 a string
513 a boolean value
514 a ParameterSet object
515 a Simulation object
516 a Model object
517 any other object (not considered by writeESySXML and writeXML)
518
519 Example how to create an ESySParameters object:
520
521 p11=ParameterSet(gamma1=1.,gamma2=2.,gamma3=3.)
522 p1=ParameterSet(dim=2,tol_v=0.001,output_file="/tmp/u.%3.3d.dx",runFlag=True,parm11=p11)
523 parm=ParameterSet(parm1=p1,parm2=ParameterSet(alpha=Link(p11,"gamma1")))
524
525 This can be accessed as
526
527 parm.parm1.gamma=0.
528 parm.parm1.dim=2
529 parm.parm1.tol_v=0.001
530 parm.parm1.output_file="/tmp/u.%3.3d.dx"
531 parm.parm1.runFlag=True
532 parm.parm1.parm11.gamma1=1.
533 parm.parm1.parm11.gamma2=2.
534 parm.parm1.parm11.gamma3=3.
535 parm.parm2.alpha=1. (value of parm.parm1.parm11.gamma1)
536
537 """
538 def __init__(self, parameters=[], **kwargs):
539 """creates a ParameterSet with parameters parameters"""
540 LinkableObject.__init__(self, **kwargs)
541 self.parameters = set()
542 self.declareParameters(parameters)
543
544 def declareParameter(self,**parameters):
545 """declares a new parameter(s) and its (their) inital value."""
546 self.declareParameters(parameters)
547
548 def declareParameters(self,parameters):
549 """declares a set of parameters. parameters can be a list, a dictonary or a ParameterSet."""
550 if isinstance(parameters,ListType):
551 parameters = zip(parameters, itertools.repeat(None))
552 if isinstance(parameters,DictType):
553 parameters = parameters.iteritems()
554
555 for prm, value in parameters:
556 setattr(self,prm,value)
557 self.parameters.add(prm)
558
559 self.trace("%s: parameter %s has been declared."%(self,prm))
560
561 def releaseParameters(self,name):
562 """removes parameter name from the paramameters"""
563 if self.isParameter(name):
564 self.parameters.remove(name)
565 self.debug("%s: parameter %s has been removed."%(self, name))
566
567 def __iter__(self):
568 """creates an iterator over the parameter and their values"""
569 return _ParameterIterator(self)
570
571 def showParameters(self):
572 """returns a descrition of the parameters"""
573 out="{"
574 notfirst=False
575 for i,v in self:
576 if notfirst: out=out+","
577 notfirst=True
578 if isinstance(v,ParameterSet):
579 out="%s\"%s\" : %s"%(out,i,v.showParameters())
580 else:
581 out="%s\"%s\" : %s"%(out,i,v)
582 return out+"}"
583
584 def __delattr__(self,name):
585 """removes the attribute name."""
586 LinkableObject.__delattr__(self,name)
587 try:
588 self.releaseParameter(name)
589 except:
590 pass
591
592 def toDom(self, document, node):
593 """ toDom method of ParameterSet class """
594 pset = document.createElement('ParameterSet')
595 node.appendChild(pset)
596 self._parametersToDom(document, pset)
597
598 def _parametersToDom(self, document, node):
599 node.setAttribute ('id', str(self.id))
600 for name,value in self:
601 param = document.createElement('Parameter')
602 param.setAttribute('type', value.__class__.__name__)
603
604 param.appendChild(dataNode(document, 'Name', name))
605
606 val = document.createElement('Value')
607
608 if isinstance(value,ParameterSet):
609 value.toDom(document, val)
610 param.appendChild(val)
611 elif isinstance(value, Link):
612 value.toDom(document, val)
613 param.appendChild(val)
614 elif isinstance(value,StringType):
615 param.appendChild(dataNode(document, 'Value', value))
616 else:
617 param.appendChild(dataNode(document, 'Value', str(value)))
618
619 node.appendChild(param)
620
621
622 def fromDom(cls, doc):
623
624 # Define a host of helper functions to assist us.
625 def _children(node):
626 """
627 Remove the empty nodes from the children of this node
628 """
629 return [x for x in node.childNodes
630 if not isinstance(x, minidom.Text) or x.nodeValue.strip()]
631
632 def _floatfromValue(doc):
633 return float(doc.nodeValue.strip())
634
635 def _stringfromValue(doc):
636 return str(doc.nodeValue.strip())
637
638 # Mapping from text types in the xml to methods used to process trees of that type
639 ptypemap = {"Simulation": Simulation.fromDom,
640 "Model":Model.fromDom,
641 "ParameterSet":ParameterSet.fromDom,
642 "Link":Link.fromDom,
643 "float":_floatfromValue,
644 #"int":_intfromValue,
645 "str":_stringfromValue,
646 #"bool":_boolfromValue
647 }
648
649 parameters = {}
650 for node in _children(doc):
651 ptype = node.getAttribute("type")
652
653 pname = pvalue = None
654 for childnode in _children(node):
655
656 if childnode.tagName == "Name":
657 pname = childnode.firstChild.nodeValue.strip()
658
659 if childnode.tagName == "Value":
660 nodes = _children(childnode)
661 pvalue = ptypemap[ptype](nodes[0])
662
663 parameters[pname] = pvalue
664
665 # Create the instance of ParameterSet
666 o = cls()
667 o.declareParameters(parameters)
668 registerLinkableObject(doc.getAttribute("id"), o)
669 return o
670
671 fromDom = classmethod(fromDom)
672
673 def writeXML(self,ostream=stdout):
674 """writes the object as an XML object into an output stream"""
675 # ParameterSet(d) with d[Name]=Value
676 document, node = esysDoc()
677 self.toDom(document, node)
678 ostream.write(document.toprettyxml())
679
680 class Model(ParameterSet,SimulationFrame):
681 """a Model is a SimulationFrame which is also a ParameterSet."""
682
683 def __init__(self,parameters=[],*args,**kwargs):
684 """creates a model"""
685 ParameterSet.__init__(self, parameters=parameters)
686 SimulationFrame.__init__(self,*args,**kwargs)
687
688 def toDom(self, document, node):
689 """ toDom method of Model class """
690 pset = document.createElement('Model')
691 pset.setAttribute('type', self.__class__.__name__)
692 node.appendChild(pset)
693 self._parametersToDom(document, pset)
694
695
696 class IterationDivergenceError(Exception):
697 """excpetion which should be thrown if there is no convergence of the iteration process at a time step but there is a chance taht a smaller step could help
698 to reach convergence."""
699 pass
700
701 class IterationBreakDownError(Exception):
702 """excpetion which should be thrown if there is no conevregence and there is no chance that a time step reduction would help"""
703 pass
704
705 class FailedTimeStepError(Exception):
706 """excpetion which should be thrown if the time step fails because of a step size that have been choosen to be to large"""
707 pass
708
709 class NonPositiveStepSizeError(Exception):
710 """excpetion which is thrown if the step size is not positive"""
711 pass
712
713 #
714 # ignore this text:
715 #
716 """ the Model class provides a framework to run a time-dependent simulation. A
717 Model has a set of parameter which may be fixed or altered by the Model itself
718 or other Models over time.
719
720 The parameters of a models are declared at instantion, e.g.
721
722 m=Model({"message" : "none" })
723
724 creates a Model with parameters p1 and p2 with inital values 1 and 2.
725 Typically a particular model is defined as a subclass of Model:
726
727 class Messenger(Model):
728 def __init__(self):
729 Model.__init__(self,parameters={"message" : "none" })
730
731 m=MyModel()
732
733 There are various ways how model parameters can be changed:
734
735 1) use object attributes:
736
737 m.message="Hello World!"
738
739 2) use setParamter method
740
741
742 m.setParameters(message="Hello World!")
743
744 3) or dictonaries
745
746 d={ message : "Hello World!" }
747 m.setParameters(**d)
748
749
750 A model executed buy staring the run method of the model:
751
752 m=Messenger()
753 m.run()
754
755 The run methods marches through time. It first calls the
756 doInitialization() method of the Model to set up the process. In each
757 time step the doStep() method is called to get from the current to the
758 next time step. The step size is defined by calling the
759 getSafeTimeStepSize() method. the time integration process is
760 terminated when the finalize() methods return true. Final the
761 doFinalization() method is called to finalize the process. To implement
762 a particular model a subclass of the Model class is defined. The
763 subclass overwrites the default methods of Model.
764
765 The following class defines a messenger printing in the doStep method
766 what ever the current value of its parameter message is:
767
768 class Messenger(Model):
769 def __init__(self):
770 Model.__init__(self,parameters={"message" : "none" })
771
772 def doInitialization(self):
773 print "I start talking now!"
774
775 def doStep(self,t):
776 print "Message (time %e) : %s "%(t,self.message)
777
778 def doFinalization(self):
779 print "I have no more to say!"
780
781 If a instance of the Messenger class is run, it will print the
782 initialization and finalization message only. This is because the
783 default method for finalize() does always returns True and therefore the
784 transition is terminated startcht away.
785
786 Following example for solving the ODE using a forward euler scheme:
787
788 u(t=0)=u0
789 u_t=a*u**2 for all 0<t<=ten
790
791 exact solution is given by u(t)=1/(1/u0-a*t)
792
793 class Ode1(Model):
794 def __init__(self,**args):
795 Model.__init__(self,parameters={"tend" : 1., "dt" : 0.0001 ,"a" : 0.1 ,"u" : 1. },name="test",debug=True)
796
797 def doInitialization(self):
798 self._tn=0
799
800 def doStep(self,t):
801 self.u=self.u+(t-self._tn)*self.a*self.u**2
802 self._tn=t
803
804 def doFinalization(self):
805 print "all done final error = ",abs(self.u-1./(1./3.-self.a*self._tn))
806
807 def getSafeTimeStepSize(self):
808 return self.dt
809
810 def finalize(self):
811 return self._tn>=self.tend
812
813 In some cases at a given time step an iteration process has to be
814 performed to get the state of the Model for the next time step. ` In
815 this case the doStep() method is replaced by a sequance of methods which
816 implements this iterative process. The method then will control the
817 iteration process by initializing the iteration through calling the
818 doIterationInitialization() method. The iteration is preformed by
819 calling the doIterationStep() method until the terminate() method
820 returns True. The doIterationFinalization() method is called to end the
821 iteration.
822 For a particular model these methods have to overwritten by a suitable
823 subclass without touching the doStep() method.
824
825 following example is a modification of the example above. Here an
826 implicit euler scheme is used. in each time step the problem
827
828 0= u_{n+1}-u_{n}+a*dt*u_{n+1}**2
829
830 has to be solved for u_{n+1}. The Newton scheme is used to solve this non-linear problem.
831
832
833 class Ode2(Model):
834
835 def __init__(self,**args):
836 Model.__init__(self,{"tend" : 1., "dt" : 0.1 ,"a" : 10. ,"u" : 1. , "tol " : 1.e-8},"test","bla",None,True)
837
838 def doInitialization(self):
839 self.__tn=0
840
841 def doIterationInitialization(self,t):
842 self.__iter=0
843 self.u_last=self.u
844 self.current_dt=t-self.tn
845 self.__tn=t
846
847 def doIterationStep(self):
848 self.__iter+=1
849 self.u_old=self.u
850 self.u=(self.current_dt*self.a*self.u**2-self.u_last)/(2*self.current_dt*self.a*self.u-1.)
851
852 def terminate(self):
853 return abs(self.u_old-self.u)<self.tol*abs(self.u)
854
855 def doIterationFinalization(self)
856 print "all done"
857
858 def getSafeTimeStepSize(self):
859 return self.dt
860
861 def finalize(self):
862 return self.__tn>self.tend
863
864 A model can be composed from subsimulations. Subsimulations are treated
865 as model parameters. If a model parameter is set or a value of a model
866 parameter is requested, the model will search for this parameter its
867 subsimulations in the case the model does not have this parameter
868 itself. The order in which the subsimulations are searched is critical.
869 By default a Model initializes all its subsimulations, is finalized when
870 all its subsimulations are finalized and finalizes all its
871 subsimulations. In the case an iterative process is applied on a
872 particular time step the iteration is initialized for all
873 subsimulations, then the iteration step is performed for each
874 subsimulation until all subsimulations indicate termination. Then the
875 iteration is finalized for all subsimulations. Finally teh doStop()
876 method for all submethods is called.
877
878 Here we are creating a model which groups ab instantiation of the Ode2 and the Messenger Model
879
880 o=Ode2()
881 m=Messenger()
882 om=Model(subsimulations=[o,m],debug=True)
883 om.dt=0.01
884 om.u=1.
885 m.message="it's me!"
886 om.run()
887
888 Notice that dt and u are parameters of class Ode2 and message is a
889 parameter of the Messenger class. The Model formed from these models
890 automatically hand the assignment of new values down to the
891 subsimulation. om.run() starts this combined model where now the
892 soStep() method of the Messenger object printing the value of its
893 parameter message together with a time stamp is executed in each time
894 step introduced by the Ode2 model.
895
896 A parameter of a Model can be linked to an attribute of onother object,
897 typically an parameter of another Model object.
898
899
900 which is comprised by a set of subsimulations.
901 The simulation is run through its run method which in the simplest case has the form:
902
903 s=Model()
904 s.run()
905
906 The run has an initializion and finalization phase. The latter is called
907 if all subsimulations are to be finalized. The simulation is processing
908 in time through calling the stepForward methods which updates the
909 observables of each subsimulation. A time steps size which is save for
910 all subsimulation is choosen.
911
912 At given time step an iterative process may be performed to make sure
913 that all observables are consistent across all subsimulations. In this
914 case, similar the time dependence, an initialization and finalization of
915 the iteration is performed.
916
917 A Model has input and output parameters where each input parameter can
918 be constant, time dependent or may depend on an output parameter of
919 another model or the model itself. To create a parameter name of a model
920 and to assign a value to it one can use the statement
921
922 model.name=object
923
924
925 At any time the current value of the parameter name can be obtained by
926
927 value=model.name
928
929 If the object that has been assigned to the paramter/attribute name has
930 the attribute/parameter name isself the current value of this attribute
931 of the object is returned (e.g. for model.name=object where object has
932 an attribute name, the statement value=model.name whould assign the
933 value object.name to value.). If the name of the parameters of a model
934 and an object don't match the setParameter method of model can be used.
935 So
936
937 model.setParameter(name,object,name_for_object)
938
939 links the parameter name of model with the parameter name_for_object of
940 object.
941
942 The run method initiates checkpointing (it is not clear how to do this
943 yet)
944 =====
945
946 """
947
948
949
950 if __name__=="__main__":
951 import math
952 #
953 # test for parameter set
954 #
955 p11=ParameterSet()
956 p11.declareParameter(gamma1=1.,gamma2=2.,gamma3=3.)
957 p1=ParameterSet()
958 p1.declareParameter(dim=2,tol_v=0.001,output_file="/tmp/u.%3.3d.dx",runFlag=True,parm11=p11)
959 parm=ParameterSet({ "parm1" : p1 , "parm2" : ParameterSet(["alpha"])})
960 parm.parm2.alpha=Link(p11,"gamma1")
961 parm.x="that should not be here!"
962 print parm.showParameters()
963 # should be something like: {"parm2" : {"alpha" : reference to attribute
964 # gamma1 of <__main__.ParameterSet instance at 0xf6db51cc>},"parm1" : {"dim"
965 # : 2,"runFlag" : True,"tol_v": 0.001,"parm11" : {"gamma3" : 3.0,"gamma2" :
966 # 2.0,"gamma1" : 1.0},"output_file" : /tmp/u.%3.3d.dx}}
967 assert parm.parm2.alpha==1.
968 parm.writeXML()
969
970 #=======================
971 class Messenger(Model):
972 def __init__(self):
973 Model.__init__(self)
974 self.declareParameter(message="none")
975
976 def doInitialization(self,t):
977 self.__t=t
978 print "I start talking now!"
979
980 def doStep(self,dt):
981 self.__t+=dt
982 print "Message (time %e) : %s "%(self.__t,self.message)
983
984 def doFinalization(self):
985 print "I have no more to say!"
986
987 class ODETEST(Model):
988 """ implements a solver for the ODE
989
990 du/dt=a*u+f(t)
991
992 we use a implicit euler scheme :
993
994 u_n-u_{n-1}= dt*a u_n + st*f(t_n)
995
996 to get u_n we run an iterative process
997
998 u_{n.k}=u_{n-1}+dt*(a u_{n.i-1} + f(t_n))
999
1000
1001 input for this model are step size dt, end time tend and a value for
1002 a, f and initial value for u. we need also a tolerance tol for a
1003 stopping criterion.
1004
1005 """
1006
1007 def __init__(self):
1008 Model.__init__(self,debug=True)
1009 self.declareParameter(tend=1.,dt=0.1,a=0.9,u=10.,f=0.,message="",tol=1.e-8)
1010
1011 def doInitialization(self,t):
1012 self.__tn=t
1013 self.__iter=0
1014
1015 def doIterationInitialization(self,dt):
1016 self.__iter=0
1017 self.__u_last=self.u
1018
1019 def doIterationStep(self,dt):
1020 self.__iter+=1
1021 self.__u_old=self.u
1022 self.u=self.__u_last+dt*(self.a*self.__u_old+self.f)
1023
1024 def terminate(self):
1025 if self.__iter<1:
1026 return False
1027 else:
1028 return abs(self.__u_old-self.u)<self.tol*abs(self.u)
1029
1030 def doIterationFinalization(self,dt):
1031 self.__tn+=dt
1032 self.message="current error = %e"%abs(self.u-10.*math.exp((self.a-1.)*self.__tn))
1033
1034 def getSafeTimeStepSize(self,dt):
1035 return min(self.dt,1./(abs(self.a)+1.))
1036
1037 def finalize(self):
1038 return self.__tn>=self.tend
1039
1040 #
1041 # s solves the coupled ODE:
1042 #
1043 # du/dt=a*u+ v
1044 # dv/dt= u+a*v
1045 #
1046 # each equation is treated through the ODETEST class. The equations are
1047 # linked and iteration over each time step is performed. the current
1048 # error of v is reported by the Messenger class.
1049 #
1050 o1=ODETEST()
1051 o1.u=10
1052 o2=ODETEST()
1053 o2.u=-10.
1054 o1.f=Link(o2,"u")
1055 o2.f=Link(o1,"u")
1056 m=Messenger()
1057 o1.dt=0.01
1058 m.message=Link(o1)
1059 s=ExplicitSimulation([Simulation([o1,o2],debug=True),m],debug=True)
1060 s.run()
1061 s.writeXML()

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26