/[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 122 - (hide annotations)
Thu Jun 9 05:38:05 2005 UTC (14 years, 4 months ago) by jgs
Original Path: trunk/esys2/escript/py_src/modelframe.py
File MIME type: text/x-python
File size: 29908 byte(s)
Merge of development branch back to main trunk on 2005-06-09

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26