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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 123 - (hide 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 jgs 121 # $Id$
2    
3 jgs 122 from types import StringType,IntType,FloatType,BooleanType,ListType,DictType
4     from sys import stdout
5 jgs 123 import itertools
6 jgs 121
7 jgs 123 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 jgs 122 from xml.dom import minidom
16 jgs 121
17 jgs 122 def dataNode(document, tagName, data):
18 jgs 123 t = document.createTextNode(str(data))
19     n = document.createElement(tagName)
20     n.appendChild(t)
21     return n
22 jgs 121
23 jgs 122 def esysDoc():
24 jgs 123 doc = minidom.Document()
25     esys = doc.createElement('ESys')
26     doc.appendChild(esys)
27     return doc, esys
28 jgs 121
29 jgs 123 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 jgs 122 class Link:
79 jgs 123 """
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 jgs 121
121 jgs 123 if callable(out):
122     return out()
123     else:
124     return out
125 jgs 121
126 jgs 123 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 jgs 121
137 jgs 123 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 jgs 121
144 jgs 123 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 jgs 122 class LinkableObject(object):
157 jgs 123 """
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 jgs 121
167 jgs 123 p.x will contain the current value of attribute name of object o.
168 jgs 121
169 jgs 123 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 jgs 122 class SimulationFrame(LinkableObject):
235 jgs 123 """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 jgs 121
331 jgs 123 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 jgs 121
340 jgs 123 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 jgs 121
363 jgs 122 class Simulation(SimulationFrame):
364 jgs 123 """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 jgs 121
371 jgs 123 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 jgs 121
393 jgs 123 def toDom(self, document, node):
394     """ toDom method of Simulation class """
395     simulation = document.createElement('Simulation')
396     simulation.setAttribute('type', self.__class__.__name__)
397 jgs 121
398 jgs 123 for rank, sim in enumerate(self.iterSubsimulations()):
399     component = document.createElement('Component')
400     component.setAttribute('rank', str(rank))
401 jgs 121
402 jgs 123 sim.toDom(document, component)
403 jgs 121
404 jgs 123 simulation.appendChild(component)
405 jgs 121
406 jgs 123 node.appendChild(simulation)
407 jgs 121
408 jgs 123 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 jgs 122 class ExplicitSimulation(Simulation):
483 jgs 123 """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 jgs 121
487 jgs 123 def doStep(self,dt):
488     """executes the time step for all subsimulation"""
489     for i in self.iterSubsimulations():
490     i.doStep(dt)
491 jgs 121
492 jgs 122 class _ParameterIterator:
493 jgs 123 def __init__(self,parameterset):
494 jgs 121
495 jgs 123 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 jgs 122 class ParameterSet(LinkableObject):
506 jgs 123 """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 jgs 121
555 jgs 123 for prm, value in parameters:
556     setattr(self,prm,value)
557     self.parameters.add(prm)
558 jgs 121
559 jgs 123 self.trace("%s: parameter %s has been declared."%(self,prm))
560 jgs 121
561 jgs 123 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 jgs 121
592 jgs 123 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 jgs 121
598 jgs 123 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 jgs 121
604 jgs 123 param.appendChild(dataNode(document, 'Name', name))
605 jgs 121
606 jgs 123 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 jgs 122 class Model(ParameterSet,SimulationFrame):
681 jgs 123 """a Model is a SimulationFrame which is also a ParameterSet."""
682 jgs 121
683 jgs 123 def __init__(self,parameters=[],*args,**kwargs):
684     """creates a model"""
685     ParameterSet.__init__(self, parameters=parameters)
686     SimulationFrame.__init__(self,*args,**kwargs)
687 jgs 121
688 jgs 123 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 jgs 121
695 jgs 123
696 jgs 122 class IterationDivergenceError(Exception):
697 jgs 123 """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 jgs 121
701 jgs 122 class IterationBreakDownError(Exception):
702 jgs 123 """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 jgs 121
705 jgs 122 class FailedTimeStepError(Exception):
706 jgs 123 """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 jgs 121
709 jgs 122 class NonPositiveStepSizeError(Exception):
710 jgs 123 """excpetion which is thrown if the step size is not positive"""
711     pass
712 jgs 121
713 jgs 122 #
714 jgs 123 # ignore this text:
715 jgs 122 #
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 jgs 123 or other Models over time.
719 jgs 121
720 jgs 123 The parameters of a models are declared at instantion, e.g.
721 jgs 121
722 jgs 123 m=Model({"message" : "none" })
723 jgs 121
724 jgs 123 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 jgs 121
727 jgs 123 class Messenger(Model):
728     def __init__(self):
729     Model.__init__(self,parameters={"message" : "none" })
730 jgs 121
731 jgs 123 m=MyModel()
732 jgs 121
733 jgs 123 There are various ways how model parameters can be changed:
734 jgs 121
735 jgs 123 1) use object attributes:
736 jgs 121
737 jgs 123 m.message="Hello World!"
738 jgs 121
739 jgs 123 2) use setParamter method
740    
741    
742     m.setParameters(message="Hello World!")
743 jgs 121
744 jgs 123 3) or dictonaries
745 jgs 122
746 jgs 123 d={ message : "Hello World!" }
747     m.setParameters(**d)
748 jgs 121
749    
750 jgs 123 A model executed buy staring the run method of the model:
751 jgs 122
752 jgs 123 m=Messenger()
753     m.run()
754 jgs 121
755 jgs 123 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 jgs 121
765 jgs 123 The following class defines a messenger printing in the doStep method
766     what ever the current value of its parameter message is:
767 jgs 121
768 jgs 123 class Messenger(Model):
769     def __init__(self):
770     Model.__init__(self,parameters={"message" : "none" })
771 jgs 121
772 jgs 123 def doInitialization(self):
773     print "I start talking now!"
774 jgs 121
775 jgs 123 def doStep(self,t):
776     print "Message (time %e) : %s "%(t,self.message)
777 jgs 121
778 jgs 123 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 jgs 121
788 jgs 123 u(t=0)=u0
789     u_t=a*u**2 for all 0<t<=ten
790 jgs 121
791 jgs 123 exact solution is given by u(t)=1/(1/u0-a*t)
792 jgs 121
793 jgs 123 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 jgs 121
797 jgs 123 def doInitialization(self):
798     self._tn=0
799 jgs 121
800 jgs 123 def doStep(self,t):
801     self.u=self.u+(t-self._tn)*self.a*self.u**2
802     self._tn=t
803 jgs 121
804 jgs 123 def doFinalization(self):
805     print "all done final error = ",abs(self.u-1./(1./3.-self.a*self._tn))
806 jgs 121
807 jgs 123 def getSafeTimeStepSize(self):
808     return self.dt
809 jgs 121
810 jgs 123 def finalize(self):
811     return self._tn>=self.tend
812 jgs 121
813 jgs 123 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 jgs 121
825 jgs 123 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 jgs 121
830 jgs 123 has to be solved for u_{n+1}. The Newton scheme is used to solve this non-linear problem.
831 jgs 121
832 jgs 122
833 jgs 123 class Ode2(Model):
834 jgs 121
835 jgs 123 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 jgs 121
838 jgs 123 def doInitialization(self):
839     self.__tn=0
840 jgs 121
841 jgs 123 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 jgs 121
847 jgs 123 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 jgs 121
852 jgs 123 def terminate(self):
853     return abs(self.u_old-self.u)<self.tol*abs(self.u)
854 jgs 121
855 jgs 123 def doIterationFinalization(self)
856     print "all done"
857 jgs 121
858 jgs 123 def getSafeTimeStepSize(self):
859     return self.dt
860 jgs 121
861 jgs 123 def finalize(self):
862     return self.__tn>self.tend
863 jgs 121
864 jgs 123 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 jgs 121
878 jgs 123 Here we are creating a model which groups ab instantiation of the Ode2 and the Messenger Model
879 jgs 121
880 jgs 123 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 jgs 121
888 jgs 123 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 jgs 121
896 jgs 123 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 jgs 121
903 jgs 123 s=Model()
904     s.run()
905 jgs 121
906 jgs 123 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 jgs 121
912 jgs 123 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 jgs 121
917 jgs 123 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 jgs 121
922 jgs 123 model.name=object
923 jgs 121
924    
925 jgs 123 At any time the current value of the parameter name can be obtained by
926 jgs 121
927 jgs 123 value=model.name
928 jgs 121
929 jgs 123 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 jgs 121
937 jgs 123 model.setParameter(name,object,name_for_object)
938 jgs 121
939 jgs 123 links the parameter name of model with the parameter name_for_object of
940     object.
941 jgs 121
942 jgs 123 The run method initiates checkpointing (it is not clear how to do this
943     yet)
944 jgs 122 =====
945 jgs 123
946 jgs 122 """
947 jgs 121
948    
949 jgs 122
950 jgs 121 if __name__=="__main__":
951 jgs 123 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 jgs 121
976 jgs 123 def doInitialization(self,t):
977     self.__t=t
978     print "I start talking now!"
979 jgs 121
980 jgs 123 def doStep(self,dt):
981     self.__t+=dt
982     print "Message (time %e) : %s "%(self.__t,self.message)
983 jgs 121
984 jgs 123 def doFinalization(self):
985     print "I have no more to say!"
986 jgs 121
987 jgs 123 class ODETEST(Model):
988     """ implements a solver for the ODE
989 jgs 121
990 jgs 123 du/dt=a*u+f(t)
991 jgs 121
992 jgs 123 we use a implicit euler scheme :
993 jgs 121
994 jgs 123 u_n-u_{n-1}= dt*a u_n + st*f(t_n)
995 jgs 121
996 jgs 123 to get u_n we run an iterative process
997 jgs 122
998 jgs 123 u_{n.k}=u_{n-1}+dt*(a u_{n.i-1} + f(t_n))
999 jgs 121
1000    
1001 jgs 123 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 jgs 121
1005 jgs 123 """
1006 jgs 121
1007 jgs 123 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 jgs 121
1011 jgs 123 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 jgs 121
1019 jgs 123 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 jgs 121
1030 jgs 123 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 jgs 121
1034 jgs 123 def getSafeTimeStepSize(self,dt):
1035     return min(self.dt,1./(abs(self.a)+1.))
1036 jgs 121
1037 jgs 123 def finalize(self):
1038     return self.__tn>=self.tend
1039 jgs 121
1040 jgs 123 #
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