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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 126 - (hide annotations)
Fri Jul 22 03:53:08 2005 UTC (14 years, 4 months ago) by jgs
File MIME type: text/x-python
File size: 38044 byte(s)
Merge of development branch back to main trunk on 2005-07-22

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26