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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/esys2/escript/py_src/modelframe.py revision 142 by jgs, Mon Jul 25 05:28:20 2005 UTC trunk/escript/py_src/modelframe.py revision 939 by gross, Thu Jan 25 04:23:38 2007 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
2    
3    """
4    Environment for implementing models in escript
5    
6    @var __author__: name of author
7    @var __copyright__: copyrights
8    @var __license__: licence agreement
9    @var __url__: url entry point on documentation
10    @var __version__: version
11    @var __date__: date of the version
12    """
13    
14    __author__="Lutz Gross, l.gross@uq.edu.au"
15    __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
16                        http://www.access.edu.au
17                    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19                 http://www.opensource.org/licenses/osl-3.0.php"""
20    __url__="http://www.iservo.edu.au/esys"
21    __version__="$Revision$"
22    __date__="$Date$"
23    
24    
25  from types import StringType,IntType,FloatType,BooleanType,ListType,DictType  from types import StringType,IntType,FloatType,BooleanType,ListType,DictType
26  from sys import stdout  from sys import stdout
27    import numarray
28    import operator
29  import itertools  import itertools
 # import modellib  temporarily removed!!!  
30    
31  # import the 'set' module if it's not defined (python2.3/2.4 difference)  # import the 'set' module if it's not defined (python2.3/2.4 difference)
32  try:  try:
# Line 13  except NameError: Line 36  except NameError:
36    
37  from xml.dom import minidom  from xml.dom import minidom
38    
 def dataNode(document, tagName, data):  
     """  
     dataNodes are the building blocks of the xml documents constructed in  
     this module. document is the current xml document, tagName is the  
     associated xml tag, and data is the values in the tag.  
     """  
     t = document.createTextNode(str(data))  
     n = document.createElement(tagName)  
     n.appendChild(t)  
     return n  
   
 def esysDoc():  
     """  
     Global method for creating an instance of an EsysXML document.  
     """  
     doc = minidom.Document()  
     esys = doc.createElement('ESys')  
     doc.appendChild(esys)  
     return doc, esys  
39    
40  def all(seq):  def all(seq):
41      for x in seq:      for x in seq:
# Line 39  def all(seq): Line 43  def all(seq):
43              return False              return False
44      return True      return True
45    
46  LinkableObjectRegistry = {}  def any(seq):
47        for x in seq:
48  def registerLinkableObject(obj_id, o):          if x:
49      LinkableObjectRegistry[obj_id] = o              return True
50        return False
51    
52    def importName(modulename, name):
53        """ Import a named object from a module in the context of this function,
54            which means you should use fully qualified module paths.
55            Return None on failure.
56    
57  LinkRegistry = []          This function from: http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/52241
58        """
59        module = __import__(modulename, globals(), locals(), [name])
60            
61        try:
62            return vars(module)[name]
63        except KeyError:
64            raise ImportError("Could not import %s from %s" % (name, modulename))
65    
66  def registerLink(obj_id, l):  class ESySXMLParser(object):
67      LinkRegistry.append((obj_id,l))      """
68        parser for ESysXML file
69  def parse(xml):      """
70      global LinkRegistry, LinkableObjectRegistry      def __init__(self,xml, debug=False):
71      LinkRegistry = []        self.__dom = minidom.parseString(xml)
72      LinkableObjectRegistry = {}        self.__linkable_object_registry= {}
73          self.__link_registry=  []
74      doc = minidom.parseString(xml)        self.__esys=self.__dom.firstChild
75          self.debug=debug
76      sim = getComponent(doc.firstChild)    
77      for obj_id, link in LinkRegistry:      def getClassPath(self, node):
78          link.target = LinkableObjectRegistry[obj_id]          type = node.getAttribute("type")
79            if (node.getAttribute("module")):
80      return sim              module = node.getAttribute("module")
81                return importName(module, type)
82  def getComponent(doc):          else:
83      for node in doc.childNodes:              return importName("__main__", type)
         if isinstance(node, minidom.Element):  
             if node.tagName == 'Simulation':  
                 if node.getAttribute("type") == 'Simulation':  
                     return Simulation.fromDom(node)  
                 elif node.getAttribute("type") == 'ExplicitSimulation':  
                     return ExplicitSimulation.fromDom(node)  
             if node.tagName == 'Model':  
                 model_type = node.getAttribute("type")  
                 model_subclasses = Model.__subclasses__()  
                 for model in model_subclasses:  
                     if model_type == model.__name__:  
                         return model.fromDom(node)  
                       
             raise "Invalid simulation type, %r" % node.getAttribute("type")  
84    
85      raise ValueError("No Simulation Found")      def setLinks(self):
86                        for obj_id, link in self.__link_registry:
87                link.target = self.__linkable_object_registry[obj_id]
88    
89        def parse(self):
90           """
91           parser method for EsysXML and returns the list of generating ParameterSets
92           """
93           found=[]
94           for node in self.__esys.childNodes:
95               if isinstance(node, minidom.Element):
96                   if node.tagName == 'Simulation':
97                            found.append(Simulation.fromDom(self, node))
98                   elif node.tagName == 'Model':
99                            found.append(self.getClassPath(node).fromDom(self, node))
100                   elif node.tagName == 'ParameterSet':
101                            found.append(self.getClassPath(node).fromDom(self, node))
102                   else:
103                      raise "Invalid type, %r" % node.getAttribute("type")
104           self.setLinks()
105           return found
106    
107        def registerLink(self,obj_id, link):
108            self.__link_registry.append((int(obj_id),link))
109    
110        def registerLinkableObject(self,obj, node):
111            id_str=node.getAttribute('id').strip()
112            if len(id_str)>0:
113               id=int(id_str)
114               if self.__linkable_object_registry.has_key(id):
115                   raise ValueError("Object id %s already exists."%id)
116               else:
117                   self.__linkable_object_registry[id]=obj
118    
119        def getComponent(self, node):
120           """
121           returns a single component + rank from a simulation
122           parser method for EsysXML and returns the list of generating ParameterSets
123           """
124           rank  = int(node.getAttribute("rank"))
125           for n in node.childNodes:
126               if isinstance(n, minidom.Element):
127                   if n.tagName == 'Simulation':
128                            return (rank, Simulation.fromDom(self, n))
129                   elif n.tagName == 'Model':
130                            return (rank, self.getClassPath(n).fromDom(self, n))
131                   elif n.tagName == 'ParameterSet':
132                            return (rank, self.getClassPath(n).fromDom(self, n))
133                   else:
134                     raise ValueError("illegal component type %s"%n.tagName)
135           raise ValueError("cannot resolve Component")
136    
137    class ESySXMLCreator(object):
138        """
139        creates an XML Dom representation
140        """
141        def __init__(self):
142           self.__dom=minidom.Document()
143           self.__esys =self.__dom.createElement('ESys')
144           self.__dom.appendChild(self.__esys)
145           self.__linkable_object_registry={}
146           self.__number_sequence = itertools.count(100)
147        def getRoot(self):
148           return self.__esys
149        def createElement(self,name):
150          return self.__dom.createElement(name)
151        def createTextNode(self,name):
152          return self.__dom.createTextNode(name)
153        def getElementById(self,name):
154          return self.__dom.getElementById(name)
155        def createDataNode(self, tagName, data):
156              """
157              C{createDataNode}s are the building blocks of the xml documents constructed in
158              this module.  
159        
160              @param tagName: the associated xml tag
161              @param data: the values in the tag
162              """
163              n = self.createElement(tagName)
164              n.appendChild(self.createTextNode(str(data)))
165              return n
166        def getLinkableObjectId(self, obj):
167            for id, o in self.__linkable_object_registry.items():
168                if o == obj: return id
169            id =self.__number_sequence.next()
170            self.__linkable_object_registry[id]=obj
171            return id
172            
173        def registerLinkableObject(self, obj, node):
174            """
175            returns a unique object id for object obj
176            """
177            id=self.getLinkableObjectId(obj)
178            node.setAttribute('id',str(id))
179            node.setIdAttribute("id")
180    
181        def includeTargets(self):
182            target_written=True
183            while target_written:
184                targetsList =self.__dom.getElementsByTagName('Target')
185                target_written=False
186                for element in targetsList:
187                   targetId = int(element.firstChild.nodeValue.strip())
188                   if self.getElementById(str(targetId)): continue
189                   targetObj = self.__linkable_object_registry[targetId]
190                   targetObj.toDom(self, self.__esys)
191                   target_written=True
192    
193        def toprettyxml(self):
194            self.includeTargets()
195            return self.__dom.toprettyxml()
196    
197  class Link:  class Link:
198      """      """
199      a Link makes an attribute of an object callable:      A Link makes an attribute of an object callable::
200    
201            o.object()            o.object()
202            o.a=8            o.a=8
203            l=Link(o,"a")            l=Link(o,"a")
# Line 93  class Link: Line 206  class Link:
206            
207      def __init__(self,target,attribute=None):      def __init__(self,target,attribute=None):
208          """          """
209          creates a link to the object target. If attribute is given, the link is          Creates a link to the object target. If attribute is given, the link is
210          establised to this attribute of the target.  Otherwise the attribute is          establised to this attribute of the target.  Otherwise the attribute is
211          undefined.          undefined.
212          """          """
213          self.target = target          self.target = target
214          self.attribute = None          self.attribute = None
215          self.setAttributeName(attribute)          self.setAttributeName(attribute)
216    
217        def getTarget(self):
218            """
219            returns the target
220            """
221            return self.target
222        def getAttributeName(self):
223            """
224            returns the name of the attribute the link is pointing to
225            """
226            return self.attribute
227            
228      def setAttributeName(self,attribute):      def setAttributeName(self,attribute):
229          """          """
230          set a new attribute name to be collected from the target object. The          Set a new attribute name to be collected from the target object. The
231          target object must have the attribute with name attribute.          target object must have the attribute with name attribute.
232          """          """
233          if attribute and self.target and not hasattr(self.target, attribute):          if attribute and self.target:
234              raise AttributeError("%s: target %s has no attribute %s."%(self, self.target, attribute))              if isinstance(self.target,LinkableObject):
235                   if not self.target.hasAttribute(attribute):
236                      raise AttributeError("%s: target %s has no attribute %s."%(self, self.target, attribute))
237                else:
238                   if not hasattr(self.target,attribute):
239                      raise AttributeError("%s: target %s has no attribute %s."%(self, self.target, attribute))
240          self.attribute = attribute          self.attribute = attribute
241            
242      def hasDefinedAttributeName(self):      def hasDefinedAttributeName(self):
243          """          """
244          returns true if an attribute name is set          Returns true if an attribute name is set.
245          """          """
246          return self.attribute != None          return self.attribute != None
247            
248      def __repr__(self):      def __repr__(self):
249          """          """
250          returns a string representation of the link          Returns a string representation of the link.
251          """          """
252          if self.hasDefinedAttributeName():          if self.hasDefinedAttributeName():
253              return "<Link to attribute %s of %s>" % (self.attribute, self.target)              return "<Link to attribute %s of %s>" % (self.attribute, self.target)
# Line 127  class Link: Line 256  class Link:
256            
257      def __call__(self,name=None):      def __call__(self,name=None):
258          """          """
259          returns the value of the attribute of the target object. If the          Returns the value of the attribute of the target object. If the
260          atrribute is callable then the return value of the call is returned.          atrribute is callable then the return value of the call is returned.
261          """          """
262          if name:          if name:
# Line 140  class Link: Line 269  class Link:
269          else:          else:
270              return out              return out
271    
272      def toDom(self, document, node):      def toDom(self, esysxml, node):
273          """          """
274          toDom method of Link. Creates a Link node and appends it to the current XML          C{toDom} method of Link. Creates a Link node and appends it to the
275          document      current XML esysxml.
276          """          """
277          link = document.createElement('Link')          link = esysxml.createElement('Link')
278          link.appendChild(dataNode(document, 'Target', self.target.id))          assert (self.target != None), ("Target was none, name was %r" % self.attribute)
279            link.appendChild(esysxml.createDataNode('Target', esysxml.getLinkableObjectId(self.target)))
280          # this use of id will not work for purposes of being able to retrieve the intended          # this use of id will not work for purposes of being able to retrieve the intended
281          # target from the xml later. I need a better unique identifier.          # target from the xml later. I need a better unique identifier.
282          assert self.attribute, "You can't xmlify a Link without a target attribute"          assert self.attribute, "You can't xmlify a Link without a target attribute"
283          link.appendChild(dataNode(document, 'Attribute', self.attribute))          link.appendChild(esysxml.createDataNode('Attribute', self.attribute))
284          node.appendChild(link)          node.appendChild(link)
285    
286      def fromDom(cls, doc):      def fromDom(cls, esysxml, node):
287          targetid = doc.getElementsByTagName("Target")[0].firstChild.nodeValue.strip()          targetid = int(node.getElementsByTagName("Target")[0].firstChild.nodeValue.strip())
288          attribute = doc.getElementsByTagName("Attribute")[0].firstChild.nodeValue.strip()          attribute =str(node.getElementsByTagName("Attribute")[0].firstChild.nodeValue.strip())
289          l = cls(None, attribute)          l = cls(None, attribute)
290          registerLink(targetid, l)          esysxml.registerLink(targetid, l)
291          return l          return l
292    
293      fromDom = classmethod(fromDom)      fromDom = classmethod(fromDom)
294            
     def writeXML(self,ostream=stdout):  
         """  
         writes an XML representation of self to the output stream ostream.  
         If ostream is nor present the standart output stream is used.  If  
         esysheader==True the esys XML header is written  
         """  
   
         document, rootnode = esysDoc()  
         self.toDom(document, rootnode)  
   
         ostream.write(document.toprettyxml())  
   
295  class LinkableObject(object):  class LinkableObject(object):
296      """      """
297      An object that allows to link its attributes to attributes of other objects      An object that allows to link its attributes to attributes of other objects
298      via a Link object. For instance      via a Link object. For instance::
299                        
300             p = LinkableObject()             p = LinkableObject()
301             p.x = Link(o,"name")             p.x = Link(o,"name")
302             print p.x             print p.x
303            
304      links attribute x of p to the attribute name of object o.      links attribute C{x} of C{p} to the attribute name of object C{o}.
305    
306      p.x will contain the current value of attribute name of object o.        C{p.x} will contain the current value of attribute C{name} of object
307        C{o}.  
308    
309      If the value of getattr(o, "name") is callable, p.x will rturn the return      If the value of C{getattr(o, "name")} is callable, C{p.x} will return
310      value of the call.      the return value of the call.
311      """      """
312        
     number_sequence = itertools.count(100)  
313            
314      def __init__(self, debug=False):      def __init__(self, id = None, debug=False):
315          """ initializes LinkableObject so that we can operate on Links """          """
316        Initializes LinkableObject so that we can operate on Links
317        """
318          self.debug = debug          self.debug = debug
319          self.__linked_attributes={}          self.__linked_attributes={}
320          self.id = self.number_sequence.next()          
   
321      def trace(self, msg):      def trace(self, msg):
322          """ If debugging is on, print the message, otherwise do nothing          """
323        If debugging is on, print the message, otherwise do nothing
324          """          """
325          if self.debug:          if self.debug:
326              print msg              print "%s: %s"%(str(self),msg)
327            
328      def __getattr__(self,name):      def __getattr__(self,name):
329          """returns the value of attribute name. If the value is a Link object the          """
330          object is called and the return value is returned."""      Returns the value of attribute name. If the value is a Link object the
331            object is called and the return value is returned.
332        """
333          out = self.getAttributeObject(name)          out = self.getAttributeObject(name)
334          if isinstance(out,Link):          if isinstance(out,Link):
335              return out()              return out()
# Line 215  class LinkableObject(object): Line 337  class LinkableObject(object):
337              return out              return out
338            
339      def getAttributeObject(self,name):      def getAttributeObject(self,name):
340          """return the object stored for attribute name."""          """
341        Return the object stored for attribute name.
342        """
343    
344          if self.__dict__.has_key(name):          if self.__dict__.has_key(name):
345              return self.__dict__[name]              return self.__dict__[name]
# Line 223  class LinkableObject(object): Line 347  class LinkableObject(object):
347          if self.__linked_attributes.has_key(name):          if self.__linked_attributes.has_key(name):
348              return self.__linked_attributes[name]              return self.__linked_attributes[name]
349    
350            if self.__class__.__dict__.has_key(name):
351                return self.__class.__dict__[name]
352    
353          raise AttributeError,"No attribute %s."%name          raise AttributeError,"No attribute %s."%name
354            
355      def __setattr__(self,name,value):      def hasAttribute(self,name):
356          """sets the value for attribute name. If value is a Link the target          """
357          attribute is set to name if no attribute has been specified."""      Returns True if self as attribute name.
358        """
359            return self.__dict__.has_key(name) or self.__linked_attributes.has_key(name) or  self.__class__.__dict__.has_key(name)
360    
361        def __setattr__(self,name,value):
362            """
363        Sets the value for attribute name. If value is a Link the target
364            attribute is set to name if no attribute has been specified.
365        """
366    
367          if self.__dict__.has_key(name):          if self.__dict__.has_key(name):
368              del self.__dict__[name]              del self.__dict__[name]
# Line 238  class LinkableObject(object): Line 372  class LinkableObject(object):
372                  value.setAttributeName(name)                  value.setAttributeName(name)
373              self.__linked_attributes[name] = value              self.__linked_attributes[name] = value
374    
375              self.trace("DEBUG: %s: attribute %s is now linked by %s."%(self,name,value))              self.trace("attribute %s is now linked by %s."%(name,value))
376          else:          else:
377              self.__dict__[name] = value              self.__dict__[name] = value
378            
379      def __delattr__(self,name):      def __delattr__(self,name):
380          """removes the attribute name."""          """
381        Removes the attribute name.
382        """
383    
384          if self.__linked_attributes.has_key[name]:          if self.__linked_attributes.has_key[name]:
385              del self.__linked_attributes[name]              del self.__linked_attributes[name]
# Line 252  class LinkableObject(object): Line 388  class LinkableObject(object):
388          else:          else:
389              raise AttributeError,"No attribute %s."%name              raise AttributeError,"No attribute %s."%name
390    
 class SimulationFrame(LinkableObject):  
     """A SimulationFrame objects represents a processess marching over time  
     until a finalizing condition is fullfilled.  At each time step an iterative  
     process can be performed and the time step size can be controlled  
     """  
     UNDEF_DT=1.e300  
     MAX_TIME_STEP_REDUCTION=20  
     MAX_ITER_STEPS=50  
       
     def __init__(self,**kwargs):  
         """  
         Initialises a simulation  
           
         Just calls the parent constructor.  
         """  
         LinkableObject.__init__(self,**kwargs)  
       
     def doInitialization(self,t):  
         """initializes the time stepping scheme. This function may be  
         overwritten."""  
         pass  
       
     def getSafeTimeStepSize(self,dt):  
         """returns a time step size which can savely be used. This function may  
         be overwritten."""  
         return self.UNDEF_DT  
       
     def finalize(self):  
         """returns True if the time stepping is finalized. This function may be  
         overwritten."""  
         return True  
         
     def doFinalization(self):  
         """finalizes the time stepping.  This function may be overwritten."""  
         pass  
       
     def doIterationInitialization(self,dt):  
         """initializes the iteration at time step t. This function may be  
         overwritten. (only called if doStep is not overwritten)"""  
         pass  
       
     def doIterationStep(self,dt):  
         """executes the iteration step. This function may be overwritten. (only  
         called if doStep is not overwritten)"""  
         pass  
       
     def terminate(self):  
         """returns True if iteration on a time step is terminated. (only called  
         if doStep is not overwritten)"""  
         return True  
         
     def doIterationFinalization(self,dt):  
         """finalalizes the iteration process. (only called if doStep is not  
         overwritten)"""  
         pass  
       
     def run(self,check_point=None):  
         """run the simulation by performing essentially  
       
             self.doInitialization()  
             while not self.finalize():  
                dt=self.getSafeTimeStepSize()  
                self.doStep(dt)  
                if n%check_point==0: self.writeXML()  
             self.doFinalization()  
       
         """  
         self.__tn=0.  
         self.__n=0  
         self.__dt=None  
         self.doInitialization(self.__tn)  
         while not self.finalize():  
             self.__n+=1  
             self.__dt=self.getSafeTimeStepSize(self.__dt)  
             if self.__dt==None: self.__dt=self.UNDEF_DT  
             if not self.__dt>0:  
                  raise NonPositiveStepSizeError("non-positive step size in step %d",self.__n)  
             self.trace("%s: %d. time step %e (step size %e.)" %  
                           (self,self.__n,self.__tn+self.__dt,self.__dt))  
             endoftimestep=False  
             failcounter=0  
             while not endoftimestep:  
                 endoftimestep=True  
                 try:  
                     self.doStep(self.__dt)  
                 except FailedTimeStepError:  
                     self.__dt=self.getSafeTimeStepSize(self.__dt)  
                     if self.__dt==None: self.__dt=self.UNDEF_DT  
                     endoftimestep=False  
                     self.trace("%s: time step is repeated with new step size %e."%(self,self.__dt))  
                 except IterationDivergenceError:  
                     self.__dt*=0.5  
                     endoftimestep=False  
                     failcounter+=1  
                     if failcounter>self.MAX_TIME_STEP_REDUCTION:  
                         raise IterationBreakDownError("reduction of time step to achieve convergence failed.")  
   
                     self.trace("%s: iteration failes. time step is repeated with new step size %e."  
                         % (self,self.__dt))  
             self.__tn+=self.__dt  
             if not check_point==None:  
                 if self.__n%check_point==0:  
                     self.trace("%s: check point is created."%self)  
                     self.writeXML()  
         self.doFinalization()  
   
     def writeXML(self):  
         raise RuntimeError, "Not implemented"  
       
     def doStep(self,dt):  
         """executes a time step by iteration. This function may be overwritten.  
             
            basicly it performs :  
       
             self.doIterationInitialization(dt)  
             while not self.terminate(): self.doIterationStep(dt)  
             self.doIterationFinalization(dt)  
       
         """  
         self.doIterationInitialization(dt)  
         self.__iter=0  
         while not self.terminate():  
             self.trace("%s: iteration step %d"%(self,self.__iter))  
             self.doIterationStep(dt)  
             self.__iter+=1  
             if self.__iter>self.MAX_ITER_STEPS:  
                 raise IterationDivergenceError("%s: divergence in step %d"%(self,self.__iter))  
         self.doIterationFinalization(dt)  
   
 class Simulation(SimulationFrame):  
     """A Simulation object is comprised by SimulationFrame(s) called subsimulations."""  
       
     def __init__(self, subsimulations=[], **kwargs):  
         """initiates a simulation from a list of subsimulations. """  
         SimulationFrame.__init__(self, **kwargs)  
         self.__subsimulations=[]  
   
         for i in range(len(subsimulations)):  
             self[i] = subsimulations[i]  
       
     def iterSubsimulations(self):  
         """returns an iterator over the subsimulations"""  
         return self.__subsimulations  
       
     def __getitem__(self,i):  
         """returns the i-th subsimulation"""  
         return self.__subsimulations[i]  
       
     def __setitem__(self,i,value):  
         """sets the i-th subsimulation"""  
         if not isinstance(value,SimulationFrame):  
             raise ValueError("assigned value is not a Simulation")  
         for j in range(max(i-len(self.__subsimulations)+1,0)): self.__subsimulations.append(None)  
         self.__subsimulations[i]=value  
       
     def __len__(self):  
         """returns the number of subsimulations"""  
         return len(self.__subsimulations)  
   
     def toDom(self, document, node):  
         """ toDom method of Simulation class """  
         simulation = document.createElement('Simulation')  
         simulation.setAttribute('type', self.__class__.__name__)  
   
         for rank, sim in enumerate(self.iterSubsimulations()):  
             component = document.createElement('Component')  
             component.setAttribute('rank', str(rank))  
   
             sim.toDom(document, component)  
   
             simulation.appendChild(component)  
   
         node.appendChild(simulation)  
   
     def writeXML(self,ostream=stdout):  
         """writes the object as an XML object into an output stream"""  
         document, rootnode = esysDoc()  
         self.toDom(document, rootnode)  
         ostream.write(document.toprettyxml())  
       
     def getSafeTimeStepSize(self,dt):  
         """returns a time step size which can safely be used by all subsimulations"""  
         out=self.UNDEF_DT  
         for o in self.iterSubsimulations():  
             dt_new = o.getSafeTimeStepSize(dt)  
             if dt_new != None:  
                 out = min(out,dt_new)  
   
         return out  
       
     def doInitialization(self,dt):  
         """initializes all subsimulations """  
         for o in self.iterSubsimulations():  
             o.doInitialization(dt)  
       
     def finalize(self):  
         """returns True if all subsimulations are finalized"""  
         return all([o.finalize() for o in self.iterSubsimulations()])  
         
     def doFinalization(self):  
         """finalalizes the time stepping for all subsimulations."""  
         for i in self.iterSubsimulations(): i.doFinalization()  
       
     def doIterationInitialization(self,dt):  
         """initializes the iteration at time t for all subsimulations."""  
         self.__iter=0  
         self.trace("%s: iteration starts"%self)  
   
         for o in self.iterSubsimulations():  
             o.doIterationInitialization(dt)  
       
     def terminate(self):  
         """returns True if all iterations for all subsimulations are terminated."""  
         return all([o.terminate() for o in self.iterSubsimulations()])  
         
     def doIterationFinalization(self,dt):  
         """finalalizes the iteration process for each of the subsimulations."""  
         for o in self.iterSubsimulations():  
             o.doIterationFinalization(dt)  
         self.trace("%s: iteration finalized after %s steps"%(self,self.__iter+1))  
       
     def doIterationStep(self,dt):  
         """executes the iteration step at time step for each subsimulation"""  
         self.__iter+=1  
         self.trace("%s: iteration step %d"%(self,self.__iter))  
         for o in self.iterSubsimulations():  
             o.doIterationStep(dt)  
   
     def fromDom(cls, doc):  
         """  
         Needs to be filled out.  
         """  
         sims = []  
         for node in doc.childNodes:  
             if isinstance(node, minidom.Text):  
                 continue  
   
             sims.append(getComponent(node))  
   
   
         return cls(sims)  
   
   
               
   
   
     fromDom = classmethod(fromDom)  
   
 class ExplicitSimulation(Simulation):  
     """This is a modified form of the Simulation class. In fact it overwrites  
     the doStep method by executing the doStep method of all subsimulations  
     rather then iterating over all subsimulations."""  
   
     def doStep(self,dt):  
         """executes the time step for all subsimulation"""  
         for i in self.iterSubsimulations():  
             i.doStep(dt)  
   
391  class _ParameterIterator:  class _ParameterIterator:
392      def __init__(self,parameterset):      def __init__(self,parameterset):
393    
# Line 524  class _ParameterIterator: Line 402  class _ParameterIterator:
402          return self          return self
403    
404  class ParameterSet(LinkableObject):  class ParameterSet(LinkableObject):
405      """a class which allows to emphazise attributes to be written and read to XML      """
406        A class which allows to emphazise attributes to be written and read to XML
407                
408         Leaves of  an ESySParameters objects can be      Leaves of an ESySParameters object can be:
409            
410              a real number       - a real number
411              a integer number       - a integer number
412              a string       - a string
413              a boolean value       - a boolean value
414              a ParameterSet object       - a ParameterSet object
415              a Simulation object       - a Simulation object
416              a Model object       - a Model object
417              any other object (not considered by writeESySXML and writeXML)       - a numarray object
418                 - a list of booleans
419             Example how to create an ESySParameters object:          - any other object (not considered by writeESySXML and writeXML)
420            
421                   p11=ParameterSet(gamma1=1.,gamma2=2.,gamma3=3.)      Example how to create an ESySParameters object::
422                   p1=ParameterSet(dim=2,tol_v=0.001,output_file="/tmp/u.%3.3d.dx",runFlag=True,parm11=p11)      
423                   parm=ParameterSet(parm1=p1,parm2=ParameterSet(alpha=Link(p11,"gamma1")))          p11=ParameterSet(gamma1=1.,gamma2=2.,gamma3=3.)
424                p1=ParameterSet(dim=2,tol_v=0.001,output_file="/tmp/u.%3.3d.dx",runFlag=True,parm11=p11)
425             This can be accessed as          parm=ParameterSet(parm1=p1,parm2=ParameterSet(alpha=Link(p11,"gamma1")))
426            
427                   parm.parm1.gamma=0.      This can be accessed as::
428                   parm.parm1.dim=2      
429                   parm.parm1.tol_v=0.001      parm.parm1.gamma=0.
430                   parm.parm1.output_file="/tmp/u.%3.3d.dx"          parm.parm1.dim=2
431                   parm.parm1.runFlag=True          parm.parm1.tol_v=0.001
432                   parm.parm1.parm11.gamma1=1.          parm.parm1.output_file="/tmp/u.%3.3d.dx"
433                   parm.parm1.parm11.gamma2=2.          parm.parm1.runFlag=True
434                   parm.parm1.parm11.gamma3=3.          parm.parm1.parm11.gamma1=1.
435                   parm.parm2.alpha=1. (value of parm.parm1.parm11.gamma1)          parm.parm1.parm11.gamma2=2.
436                        parm.parm1.parm11.gamma3=3.
437            parm.parm2.alpha=1. (value of parm.parm1.parm11.gamma1)
438      """      """
439      def __init__(self, parameters=[], **kwargs):      def __init__(self, parameters=[], **kwargs):
440          """creates a ParameterSet with parameters parameters"""          """
441        Creates a ParameterSet with parameters parameters.
442        """
443          LinkableObject.__init__(self, **kwargs)          LinkableObject.__init__(self, **kwargs)
444          self.parameters = set()          self.parameters = set()
445          self.declareParameters(parameters)          self.declareParameters(parameters)
446    
447        def __repr__(self):
448            return "<%s %d>"%(self.__class__.__name__,id(self))
449            
450      def declareParameter(self,**parameters):      def declareParameter(self,**parameters):
451          """declares a new parameter(s) and its (their) inital value."""          """
452        Declares a new parameter(s) and its (their) initial value.
453        """
454          self.declareParameters(parameters)          self.declareParameters(parameters)
455            
456      def declareParameters(self,parameters):      def declareParameters(self,parameters):
457          """declares a set of parameters. parameters can be a list, a dictonary or a ParameterSet."""          """
458        Declares a set of parameters. parameters can be a list, a dictionary
459        or a ParameterSet.
460        """
461          if isinstance(parameters,ListType):          if isinstance(parameters,ListType):
462              parameters = zip(parameters, itertools.repeat(None))              parameters = zip(parameters, itertools.repeat(None))
463          if isinstance(parameters,DictType):          if isinstance(parameters,DictType):
# Line 577  class ParameterSet(LinkableObject): Line 467  class ParameterSet(LinkableObject):
467              setattr(self,prm,value)              setattr(self,prm,value)
468              self.parameters.add(prm)              self.parameters.add(prm)
469    
470              self.trace("%s: parameter %s has been declared."%(self,prm))              self.trace("parameter %s has been declared."%prm)
471    
472      def releaseParameters(self,name):      def releaseParameters(self,name):
473          """removes parameter name from the paramameters"""          """
474        Removes parameter name from the paramameters.
475        """
476          if self.isParameter(name):          if self.isParameter(name):
477              self.parameters.remove(name)              self.parameters.remove(name)
478              self.debug("%s: parameter %s has been removed."%(self, name))              self.trace("parameter %s has been removed."%name)
479            
480      def __iter__(self):      def __iter__(self):
481          """creates an iterator over the parameter and their values"""          """
482        Creates an iterator over the parameter and their values.
483        """
484          return _ParameterIterator(self)          return _ParameterIterator(self)
485            
486      def showParameters(self):      def showParameters(self):
487          """returns a descrition of the parameters"""                  """
488        Returns a descrition of the parameters.
489        """        
490          out="{"          out="{"
491          notfirst=False          notfirst=False
492          for i,v in self:          for i,v in self:
# Line 603  class ParameterSet(LinkableObject): Line 499  class ParameterSet(LinkableObject):
499          return out+"}"          return out+"}"
500            
501      def __delattr__(self,name):      def __delattr__(self,name):
502          """removes the attribute name."""          """
503        Removes the attribute name.
504        """
505          LinkableObject.__delattr__(self,name)          LinkableObject.__delattr__(self,name)
506          try:          try:
507              self.releaseParameter(name)              self.releaseParameter(name)
508          except:          except:
509              pass              pass
510    
511      def toDom(self, document, node):      def toDom(self, esysxml, node):
512          """ toDom method of ParameterSet class """          """
513          pset = document.createElement('ParameterSet')      C{toDom} method of Model class
514        """
515            pset = esysxml.createElement('ParameterSet')
516            pset.setAttribute('type', self.__class__.__name__)
517            pset.setAttribute('module', self.__class__.__module__)
518            esysxml.registerLinkableObject(self, pset)
519            self._parametersToDom(esysxml, pset)
520          node.appendChild(pset)          node.appendChild(pset)
         self._parametersToDom(document, pset)  
521    
522      def _parametersToDom(self, document, node):      def _parametersToDom(self, esysxml, node):
         node.setAttribute ('id', str(self.id))  
523          for name,value in self:          for name,value in self:
524              param = document.createElement('Parameter')              # convert list to numarray when possible:
525              param.setAttribute('type', value.__class__.__name__)              if isinstance (value, list):
526                    elem_type=-1
527                    for i in value:
528                        if isinstance(i, bool):
529                            elem_type = max(elem_type,0)
530                        if isinstance(i, int) and not isinstance(i, bool):
531                            elem_type = max(elem_type,1)
532                        if isinstance(i, float):
533                            elem_type = max(elem_type,2)
534                    if elem_type == 0: value = numarray.array(value,numarray.Bool)
535                    if elem_type == 1: value = numarray.array(value,numarray.Int)
536                    if elem_type == 2: value = numarray.array(value,numarray.Float)
537    
538              param.appendChild(dataNode(document, 'Name', name))              param = esysxml.createElement('Parameter')
539                param.setAttribute('type', value.__class__.__name__)
540    
541              val = document.createElement('Value')              param.appendChild(esysxml.createDataNode('Name', name))
542    
543              if isinstance(value,ParameterSet):              val = esysxml.createElement('Value')
544                  value.toDom(document, val)              if isinstance(value,(ParameterSet,Link,DataSource)):
545                    value.toDom(esysxml, val)
546                  param.appendChild(val)                  param.appendChild(val)
547              elif isinstance(value, Link):              elif isinstance(value, numarray.NumArray):
548                  value.toDom(document, val)                  shape = value.getshape()
549                    if isinstance(shape, tuple):
550                        size = reduce(operator.mul, shape)
551                        shape = ' '.join(map(str, shape))
552                    else:
553                        size = shape
554                        shape = str(shape)
555    
556                    arraytype = value.type()
557                    if isinstance(arraytype, numarray.BooleanType):
558                          arraytype_str="Bool"
559                    elif isinstance(arraytype, numarray.IntegralType):
560                          arraytype_str="Int"
561                    elif isinstance(arraytype, numarray.FloatingType):
562                          arraytype_str="Float"
563                    elif isinstance(arraytype, numarray.ComplexType):
564                          arraytype_str="Complex"
565                    else:
566                          arraytype_str=str(arraytype)
567                    numarrayElement = esysxml.createElement('NumArray')
568                    numarrayElement.appendChild(esysxml.createDataNode('ArrayType', arraytype_str))
569                    numarrayElement.appendChild(esysxml.createDataNode('Shape', shape))
570                    numarrayElement.appendChild(esysxml.createDataNode('Data', ' '.join(
571                        [str(x) for x in numarray.reshape(value, size)])))
572                    val.appendChild(numarrayElement)
573                  param.appendChild(val)                  param.appendChild(val)
574              elif isinstance(value,StringType):              elif isinstance(value, list):
575                  param.appendChild(dataNode(document, 'Value', value))                  param.appendChild(esysxml.createDataNode('Value', ' '.join([str(x) for x in value]) ))
576                elif isinstance(value, (str, bool, int, float, type(None))):
577                    param.appendChild(esysxml.createDataNode('Value', str(value)))
578                elif isinstance(value, dict):
579                     dic = esysxml.createElement('dictionary')
580                     if len(value.keys())>0:
581                         dic.setAttribute('key_type', value.keys()[0].__class__.__name__)
582                         dic.setAttribute('value_type', value[value.keys()[0]].__class__.__name__)
583                     for k,v in value.items():
584                        i=esysxml.createElement('item')
585                        i.appendChild(esysxml.createDataNode('key', k))
586                        i.appendChild(esysxml.createDataNode('value', v))
587                        dic.appendChild(i)
588                     param.appendChild(dic)
589              else:              else:
590                  param.appendChild(dataNode(document, 'Value', str(value)))                  print value
591                    raise ValueError("cannot serialize %s type to XML."%str(value.__class__))
592    
593              node.appendChild(param)              node.appendChild(param)
594    
595          def fromDom(cls, esysxml, node):
     def fromDom(cls, doc):  
   
596          # Define a host of helper functions to assist us.          # Define a host of helper functions to assist us.
597          def _children(node):          def _children(node):
598              """              """
599              Remove the empty nodes from the children of this node              Remove the empty nodes from the children of this node.
600              """              """
601              return [x for x in node.childNodes              ret = []
602                      if not isinstance(x, minidom.Text) or x.nodeValue.strip()]              for x in node.childNodes:
603                    if isinstance(x, minidom.Text):
604                        if x.nodeValue.strip():
605                            ret.append(x)
606                    else:
607                        ret.append(x)
608                return ret
609    
610          def _floatfromValue(doc):          def _floatfromValue(esysxml, node):
611              return float(doc.nodeValue.strip())              return float(node.nodeValue.strip())
612    
613          def _stringfromValue(doc):          def _stringfromValue(esysxml, node):
614              return str(doc.nodeValue.strip())              return str(node.nodeValue.strip())
615                
616          def _intfromValue(doc):          def _intfromValue(esysxml, node):
617              return int(doc.nodeValue.strip())              return int(node.nodeValue.strip())
618    
619          def _boolfromValue(doc):          def _boolfromValue(esysxml, node):
620              return bool(doc.nodeValue.strip())              return _boolfromstring(node.nodeValue.strip())
621          
622            def _nonefromValue(esysxml, node):
623                return None
624    
625            def _numarrayfromValue(esysxml, node):
626                for node in _children(node):
627                    if node.tagName == 'ArrayType':
628                        arraytype = node.firstChild.nodeValue.strip()
629                    if node.tagName == 'Shape':
630                        shape = node.firstChild.nodeValue.strip()
631                        shape = [int(x) for x in shape.split()]
632                    if node.tagName == 'Data':
633                        data = node.firstChild.nodeValue.strip()
634                        data = [float(x) for x in data.split()]
635                return numarray.reshape(numarray.array(data, type=getattr(numarray, arraytype)),
636                                        shape)
637          
638            def _listfromValue(esysxml, node):
639                return [x for x in node.nodeValue.split()]
640    
641            def _boolfromstring(s):
642                if s == 'True':
643                    return True
644                else:
645                    return False
646          # Mapping from text types in the xml to methods used to process trees of that type          # Mapping from text types in the xml to methods used to process trees of that type
647          ptypemap = {"Simulation": Simulation.fromDom,          ptypemap = {"Simulation": Simulation.fromDom,
648                      "Model":Model.fromDom,                      "Model":Model.fromDom,
649                      "ParameterSet":ParameterSet.fromDom,                      "ParameterSet":ParameterSet.fromDom,
650                      "Link":Link.fromDom,                      "Link":Link.fromDom,
651                        "DataSource":DataSource.fromDom,
652                      "float":_floatfromValue,                      "float":_floatfromValue,
653                      "int":_intfromValue,                      "int":_intfromValue,
654                      "str":_stringfromValue,                      "str":_stringfromValue,
655                      "bool":_boolfromValue                      "bool":_boolfromValue,
656                        "list":_listfromValue,
657                        "NumArray":_numarrayfromValue,
658                        "NoneType":_nonefromValue,
659                      }                      }
660    
661          parameters = {}          parameters = {}
662          for node in _children(doc):          for n in _children(node):
663              ptype = node.getAttribute("type")              ptype = n.getAttribute("type")
664                if not ptypemap.has_key(ptype):
665                   raise KeyError("cannot handle parameter type %s."%ptype)
666    
667              pname = pvalue = None              pname = pvalue = None
668              for childnode in _children(node):              for childnode in _children(n):
   
669                  if childnode.tagName == "Name":                  if childnode.tagName == "Name":
670                      pname = childnode.firstChild.nodeValue.strip()                      pname = childnode.firstChild.nodeValue.strip()
671    
672                  if childnode.tagName == "Value":                  if childnode.tagName == "Value":
673                      nodes = _children(childnode)                      nodes = _children(childnode)
674                      pvalue = ptypemap[ptype](nodes[0])                      pvalue = ptypemap[ptype](esysxml, nodes[0])
675    
676              parameters[pname] = pvalue              parameters[pname] = pvalue
677    
678          # Create the instance of ParameterSet          # Create the instance of ParameterSet
679          o = cls()          o = cls(debug=esysxml.debug)
680          o.declareParameters(parameters)          o.declareParameters(parameters)
681          registerLinkableObject(doc.getAttribute("id"), o)          esysxml.registerLinkableObject(o, node)
682          return o          return o
683            
684      fromDom = classmethod(fromDom)      fromDom = classmethod(fromDom)
685        
686      def writeXML(self,ostream=stdout):      def writeXML(self,ostream=stdout):
687          """writes the object as an XML object into an output stream"""          """
688          # ParameterSet(d) with d[Name]=Value      Writes the object as an XML object into an output stream.
689          document, node = esysDoc()      """
690          self.toDom(document, node)          esysxml=ESySXMLCreator()
691          ostream.write(document.toprettyxml())          self.toDom(esysxml, esysxml.getRoot())
692            ostream.write(esysxml.toprettyxml())
693        
694    class Model(ParameterSet):
695        """
696        A Model object represents a processess marching over time until a
697        finalizing condition is fullfilled. At each time step an iterative
698        process can be performed and the time step size can be controlled. A
699        Model has the following work flow::
700              
701              doInitialization()
702              while not terminateInitialIteration(): doInitializationiStep()
703              doInitialPostprocessing()
704              while not finalize():
705                   dt=getSafeTimeStepSize(dt)
706                   doStepPreprocessing(dt)
707                   while not terminateIteration(): doStep(dt)
708                   doStepPostprocessing(dt)
709              doFinalization()
710    
711        where C{doInitialization}, C{finalize}, C{getSafeTimeStepSize},
712        C{doStepPreprocessing}, C{terminateIteration}, C{doStepPostprocessing},
713        C{doFinalization} are methods of the particular instance of a Model. The
714        default implementations of these methods have to be overwritten by the
715        subclass implementing a Model.
716        """
717    
718  class Model(ParameterSet,SimulationFrame):      UNDEF_DT=1.e300
     """a Model is a SimulationFrame which is also a ParameterSet."""  
719    
720      def __init__(self,parameters=[],**kwargs):      def __init__(self,parameters=[],**kwargs):
721          """creates a model"""          """
722          ParameterSet.__init__(self, parameters=parameters)      Creates a model.
         SimulationFrame.__init__(self,**kwargs)  
   
     def toDom(self, document, node):  
         """ toDom method of Model class """  
         pset = document.createElement('Model')  
         pset.setAttribute('type', self.__class__.__name__)  
         node.appendChild(pset)  
         self._parametersToDom(document, pset)  
   
   
 class IterationDivergenceError(Exception):  
     """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  
        to reach convergence."""  
     pass  
   
 class IterationBreakDownError(Exception):  
     """excpetion which should be thrown if there is no conevregence and there is no chance that a time step reduction would help"""  
     pass  
   
 class FailedTimeStepError(Exception):  
     """excpetion which should be thrown if the time step fails because of a step size that have been choosen to be to large"""  
     pass  
   
 class NonPositiveStepSizeError(Exception):  
     """excpetion which is thrown if the step size is not positive"""  
     pass  
   
 #  
 #   ignore this text:  
 #  
 """  
 the Model class provides a framework to run a time-dependent simulation. A  
 Model has a set of parameter which may be fixed or altered by the Model itself  
 or other Models over time.    
   
        The parameters of a models are declared at instantion, e.g.  
   
            m=Model({"message" : "none" })  
   
        creates a Model with parameters p1 and p2 with inital values 1 and 2.  
        Typically a particular model is defined as a subclass of Model:  
723    
724          class Messenger(Model):          Just calls the parent constructor.
725              def __init__(self):          """
726                 Model.__init__(self,parameters={"message" : "none" })          ParameterSet.__init__(self, parameters=parameters,**kwargs)
727    
728          m=MyModel()      def __str__(self):
729           return "<%s %d>"%(self.__class__.__name__,id(self))
730    
        There are various ways how model parameters can be changed:  
731    
732         1) use object attributes:      def doInitialization(self):
733            """
734        Initializes the time stepping scheme.  
735        
736        This function may be overwritten.
737        """
738            pass
739        def doInitialStep(self):
740            """
741        performs an iteration step in the initialization phase
742    
743            m.message="Hello World!"      This function may be overwritten.
744        """
745            pass
746    
747         2) use setParamter method      def terminateInitialIteration(self):
748                    """
749                  Returns True if iteration at the inital phase is terminated.
750            m.setParameters(message="Hello World!")      """
751            return True
752    
753         3) or dictonaries      def doInitialPostprocessing(self):
754              """
755             d={ message : "Hello World!" }      finalises the initialization iteration process
            m.setParameters(**d)  
756    
757        This function may be overwritten.
758        """
759            pass
760        
761        def getSafeTimeStepSize(self,dt):
762            """
763        Returns a time step size which can safely be used.
764    
765         A model executed buy staring the run method of the model:          C{dt} gives the previously used step size.
   
           m=Messenger()  
           m.run()  
766    
767         The run methods marches through time. It first calls the          This function may be overwritten.
768         doInitialization() method of the Model to set up the process. In each      """
769         time step the doStep() method is called to get from the current to the          return self.UNDEF_DT
770         next time step. The step size is defined by calling the      
771         getSafeTimeStepSize() method.  the time integration process is      def finalize(self):
772         terminated when the finalize() methods return true. Final the          """
773         doFinalization() method is called to finalize the process. To implement      Returns False if the time stepping is finalized.
774         a particular model a subclass of the Model class is defined. The      
775         subclass overwrites the default methods of Model.      This function may be overwritten.
776        """
777         The following class defines a messenger printing in the doStep method          return False
778         what ever the current value of its parameter message is:        
779        def doFinalization(self):
780         class Messenger(Model):          """
781              def __init__(self):      Finalizes the time stepping.
782                 Model.__init__(self,parameters={"message" : "none" })      
783        This function may be overwritten.
784        """
785            pass
786        
787        def doStepPreprocessing(self,dt):
788            """
789        Sets up a time step of step size dt.
790        
791        This function may be overwritten.
792        """
793            pass
794        
795        def doStep(self,dt):
796            """
797        Executes an iteration step at a time step.
798    
799              def doInitialization(self):          C{dt} is the currently used time step size.
                print "I start talking now!"  
800    
801              def doStep(self,t):          This function may be overwritten.
802                 print "Message (time %e) : %s "%(t,self.message)      """
803            pass
804        
805        def terminateIteration(self):
806            """
807        Returns True if iteration on a time step is terminated.
808        """
809            return True
810    
             def doFinalization(self):  
                print "I have no more to say!"  
811                
812         If a instance of the Messenger class is run, it will print the      def doStepPostprocessing(self,dt):
813         initialization and finalization message only.  This is because the          """
814         default method for finalize() does always returns True and therefore the      finalises the time step.
        transition is terminated startcht away.  
           
        Following example for solving the ODE using a forward euler scheme:  
   
                 u(t=0)=u0  
                 u_t=a*u**2       for all 0<t<=ten  
815    
816        exact solution is given by u(t)=1/(1/u0-a*t)          dt is the currently used time step size.
817    
818        class  Ode1(Model):          This function may be overwritten.
819           def __init__(self,**args):      """
820              Model.__init__(self,parameters={"tend" : 1., "dt" : 0.0001 ,"a" : 0.1 ,"u" : 1. },name="test",debug=True)          pass
   
          def doInitialization(self):  
              self._tn=0  
   
          def doStep(self,t):  
              self.u=self.u+(t-self._tn)*self.a*self.u**2  
              self._tn=t  
   
          def doFinalization(self):  
              print "all done final error = ",abs(self.u-1./(1./3.-self.a*self._tn))  
   
          def getSafeTimeStepSize(self):  
              return self.dt  
   
          def finalize(self):  
              return self._tn>=self.tend  
   
        In some cases at a given time step an iteration process has to be  
        performed to get the state of the Model for the next time step. ` In  
        this case the doStep() method is replaced by a sequance of methods which  
        implements this iterative process.  The method then will control the  
        iteration process by initializing the iteration through calling the  
        doIterationInitialization() method. The iteration is preformed by  
        calling the doIterationStep() method until the terminate() method  
        returns True. The doIterationFinalization() method is called to end the  
        iteration.  
        For a particular model these methods have to overwritten by a suitable  
        subclass without touching the doStep() method.  
821    
822         following example is a modification of the example above. Here an      def toDom(self, esysxml, node):
823         implicit euler scheme is used. in each time step the problem          """
824                  C{toDom} method of Model class
825             0= u_{n+1}-u_{n}+a*dt*u_{n+1}**2      """
826            pset = esysxml.createElement('Model')
827            pset.setAttribute('type', self.__class__.__name__)
828            pset.setAttribute('module', self.__class__.__module__)
829            esysxml.registerLinkableObject(self, pset)
830            node.appendChild(pset)
831            self._parametersToDom(esysxml, pset)
832        
833    class Simulation(Model):
834        """
835        A Simulation object is special Model which runs a sequence of Models.
836    
837         has to be solved for u_{n+1}. The Newton scheme is used to solve this non-linear problem.      The methods C{doInitialization}, C{finalize}, C{getSafeTimeStepSize},
838        C{doStepPreprocessing}, C{terminateIteration}, C{doStepPostprocessing},
839        C{doFinalization} are executing the corresponding methods of the models in
840        the simulation.
841        """
842        
843        FAILED_TIME_STEPS_MAX=20
844        MAX_ITER_STEPS=50
845        MAX_CHANGE_OF_DT=2.
846        
847        def __init__(self, models=[], **kwargs):
848            """
849        Initiates a simulation from a list of models.
850        """
851            super(Simulation, self).__init__(**kwargs)
852            for m in models:
853                if not isinstance(m, Model):
854                     raise TypeError("%s is not a subclass of Model."%m)
855            self.__models=[]
856            for i in range(len(models)):
857                self[i] = models[i]
858                
859    
860        def __repr__(self):
861        class  Ode2(Model):          """
862            Returns a string representation of the Simulation.
863            """
864            return "<Simulation %r>" % self.__models
865    
866         def __init__(self,**args):      def __str__(self):
867             Model.__init__(self,{"tend" : 1., "dt" : 0.1 ,"a" : 10. ,"u" : 1. , "tol " : 1.e-8},"test","bla",None,True)          """
868            Returning Simulation as a string.
869            """
870            return "<Simulation %d>"%id(self)
871        
872        def iterModels(self):
873            """
874        Returns an iterator over the models.
875        """
876            return self.__models
877        
878        def __getitem__(self,i):
879            """
880        Returns the i-th model.
881        """
882            return self.__models[i]
883        
884        def __setitem__(self,i,value):
885            """
886        Sets the i-th model.
887        """
888            if not isinstance(value,Model):
889                raise ValueError,"assigned value is not a Model but instance of %s"%(value.__class__.__name__,)
890            for j in range(max(i-len(self.__models)+1,0)):
891                self.__models.append(None)
892            self.__models[i]=value
893        
894        def __len__(self):
895            """
896        Returns the number of models.
897        """
898            return len(self.__models)
899    
900        def getAllModels(self):
901            """
902            returns a list of all models used in the Simulation including subsimulations
903            """
904            out=[]
905            for m in self.iterModels():
906                if isinstance(m, Simulation):
907                   out+=m.getAllModels()
908                else:
909                   out.append(m)
910            return list(set(out))
911    
912         def doInitialization(self):      def checkModelLinks(self, models):
913             self.__tn=0          """
914            returns a list of (model,parameter, target model ) if the the parameter of model
915            is linking to the target_model which is not in list of models.
916            """
917            out=[]
918            for m in self.iterModels():
919                if isinstance(m, Simulation):
920                   out+=[ (m,) + f for f in  m.checkModelLinks(models) ]
921                else:
922                  for p in m:
923                     if isinstance(p[1], Link):
924                        l_m=p[1].getTarget()
925                        if isinstance(l_m, Model) and not l_m in models: out.append( (m,p[0],l_m) )
926            return out
927    
928         def doIterationInitialization(self,t):      
929              self.__iter=0      def getSafeTimeStepSize(self,dt):
930              self.u_last=self.u                      """
931              self.current_dt=t-self.tn      Returns a time step size which can safely be used by all models.
             self.__tn=t  
   
        def doIterationStep(self):  
           self.__iter+=1  
           self.u_old=self.u  
           self.u=(self.current_dt*self.a*self.u**2-self.u_last)/(2*self.current_dt*self.a*self.u-1.)  
   
        def terminate(self):  
            return abs(self.u_old-self.u)<self.tol*abs(self.u)  
   
        def doIterationFinalization(self)  
            print "all done"  
   
        def getSafeTimeStepSize(self):  
            return self.dt  
   
        def finalize(self):  
             return self.__tn>self.tend  
   
        A model can be composed from subsimulations. Subsimulations are treated  
        as model parameters. If a model parameter is set or a value of a model  
        parameter is requested, the model will search for this parameter its  
        subsimulations in the case the model does not have this parameter  
        itself. The order in which the subsimulations are searched is critical.  
        By default a Model initializes all its subsimulations, is finalized when  
        all its subsimulations are finalized and finalizes all its  
        subsimulations. In the case an iterative process is applied on a  
        particular time step the iteration is initialized for all  
        subsimulations, then the iteration step is performed for each  
        subsimulation until all subsimulations indicate termination. Then the  
        iteration is finalized for all subsimulations. Finally teh doStop()  
        method for all submethods is called.  
   
        Here we are creating a model which groups ab instantiation of the Ode2 and the Messenger Model  
   
        o=Ode2()  
        m=Messenger()  
        om=Model(subsimulations=[o,m],debug=True)  
        om.dt=0.01  
        om.u=1.  
        m.message="it's me!"  
        om.run()  
   
        Notice that dt and u are parameters of class Ode2 and message is a  
        parameter of the Messenger class. The Model formed from these models  
        automatically hand the assignment of new values down to the  
        subsimulation. om.run() starts this combined model where now the  
        soStep() method of the Messenger object printing the value of its  
        parameter message together with a time stamp is executed in each time  
        step introduced by the Ode2 model.  
932    
933         A parameter of a Model can be linked to an attribute of onother object,          This is the minimum over the time step sizes of all models.
934         typically an parameter of another Model object.          """
935                  out=min([o.getSafeTimeStepSize(dt) for o in self.iterModels()])
936            return out
937        
938        def doInitialization(self):
939            """
940        Initializes all models.
941        """
942            self.n=0
943            self.tn=0.
944            for o in self.iterModels():
945                 o.doInitialization()
946        def doInitialStep(self):
947            """
948        performs an iteration step in the initialization step for all models
949        """
950            iter=0
951            while not self.terminateInitialIteration():
952                if iter==0: self.trace("iteration for initialization starts")
953                iter+=1
954                self.trace("iteration step %d"%(iter))
955                for o in self.iterModels():
956                     o.doInitialStep()
957                if iter>self.MAX_ITER_STEPS:
958                     raise IterationDivergenceError("initial iteration did not converge after %s steps."%iter)
959            self.trace("Initialization finalized after %s iteration steps."%iter)
960    
961        def doInitialPostprocessing(self):
962            """
963        finalises the initialization iteration process for all models.
964        """
965            for o in self.iterModels():
966                o.doInitialPostprocessing()
967        def finalize(self):
968            """
969        Returns True if any of the models is to be finalized.
970        """
971            return any([o.finalize() for o in self.iterModels()])
972                
973         which is comprised by a set of subsimulations.      def doFinalization(self):
974         The simulation is run through its run method which in the simplest case has the form:          """
975        finalises the time stepping for all models.
976            s=Model()      """
977            s.run()          for i in self.iterModels(): i.doFinalization()
978            self.trace("end of time integation.")
979         The run has an initializion and finalization phase. The latter is called      
980         if all subsimulations are to be finalized. The simulation is processing      def doStepPreprocessing(self,dt):
981         in time through calling the stepForward methods which updates the          """
982         observables of each subsimulation.  A time steps size which is save for      Initializes the time step for all models.
983         all subsimulation is choosen.      """
984            for o in self.iterModels():
985         At given time step an iterative process may be performed to make sure              o.doStepPreprocessing(dt)
986         that all observables are consistent across all subsimulations.  In this      
987         case, similar the time dependence, an initialization and finalization of      def terminateIteration(self):
988         the iteration is performed.          """
989        Returns True if all iterations for all models are terminated.
990         A Model has input and output parameters where each input parameter can      """
991         be constant, time dependent or may depend on an output parameter of          out=all([o.terminateIteration() for o in self.iterModels()])
992         another model or the model itself. To create a parameter name of a model          return out
        and to assign a value to it one can use the statement  
   
            model.name=object  
   
   
        At any time the current value of the parameter name can be obtained by  
993    
994                 value=model.name      def terminateInitialIteration(self):
995            """
996        Returns True if all initial iterations for all models are terminated.
997        """
998            out=all([o.terminateInitialIteration() for o in self.iterModels()])
999            return out
1000          
1001        def doStepPostprocessing(self,dt):
1002            """
1003        finalises the iteration process for all models.
1004        """
1005            for o in self.iterModels():
1006                o.doStepPostprocessing(dt)
1007            self.n+=1
1008            self.tn+=dt
1009        
1010        def doStep(self,dt):
1011            """
1012        Executes the iteration step at a time step for all model::
1013    
1014                self.doStepPreprocessing(dt)
1015                while not self.terminateIteration():
1016                for all models:
1017                self.doStep(dt)
1018                    self.doStepPostprocessing(dt)
1019            """
1020            self.iter=0
1021            while not self.terminateIteration():
1022                if self.iter==0: self.trace("iteration at %d-th time step %e starts"%(self.n+1,self.tn+dt))
1023                self.iter+=1
1024                self.trace("iteration step %d"%(self.iter))
1025                for o in self.iterModels():
1026                      o.doStep(dt)
1027            if self.iter>0: self.trace("iteration at %d-th time step %e finalized."%(self.n+1,self.tn+dt))
1028    
1029         If the object that has been assigned to the paramter/attribute name has      def run(self,check_point=None):
1030         the attribute/parameter name isself the current value of this attribute          """
1031         of the object is returned (e.g. for model.name=object where object has      Run the simulation by performing essentially::
1032         an attribute name, the statement value=model.name whould assign the      
1033         value object.name to value.). If the name of the parameters of a model          self.doInitialization()
1034         and an object don't match the setParameter method of model can be used.              while not self.terminateInitialIteration(): self.doInitialStep()
1035         So              self.doInitialPostprocessing()
1036            while not self.finalize():
1037                dt=self.getSafeTimeStepSize()
1038                self.doStep(dt)
1039                if n%check_point==0:
1040                self.writeXML()
1041            self.doFinalization()
1042    
1043            If one of the models in throws a C{FailedTimeStepError} exception a
1044        new time step size is computed through getSafeTimeStepSize() and the
1045        time step is repeated.
1046      
1047            If one of the models in throws a C{IterationDivergenceError}
1048        exception the time step size is halved and the time step is repeated.
1049    
1050             model.setParameter(name,object,name_for_object)          In both cases the time integration is given up after
1051        C{Simulation.FAILED_TIME_STEPS_MAX} attempts.
1052            """
1053            # check the completness of the models:
1054            # first a list of all the models involved in the simulation including subsimulations:
1055            #
1056            missing=self.checkModelLinks(self.getAllModels())
1057            if len(missing)>0:
1058                msg=""
1059                for l in missing:
1060                     msg+="\n\t"+str(l[-1])+" at "+str(self)
1061                     for i in xrange(len(l)-1): msg+="."+str(l[i])
1062                raise MissingLink("link targets missing in the Simulation: %s"%msg)
1063            #==============================
1064            self.doInitialization()
1065            self.doInitialStep()
1066            self.doInitialPostprocessing()
1067            dt=self.UNDEF_DT
1068            while not self.finalize():
1069                step_fail_counter=0
1070                iteration_fail_counter=0
1071                if self.n==0:
1072                    dt_new=self.getSafeTimeStepSize(dt)
1073                else:
1074                    dt_new=min(max(self.getSafeTimeStepSize(dt),dt/self.MAX_CHANGE_OF_DT),dt*self.MAX_CHANGE_OF_DT)
1075                self.trace("%d. time step %e (step size %e.)" % (self.n+1,self.tn+dt_new,dt_new))
1076                end_of_step=False
1077                while not end_of_step:
1078                   end_of_step=True
1079                   if not dt_new>0:
1080                      raise NonPositiveStepSizeError("non-positive step size in step %d"%(self.n+1))
1081                   try:
1082                      self.doStepPreprocessing(dt_new)
1083                      self.doStep(dt_new)
1084                      self.doStepPostprocessing(dt_new)
1085                   except IterationDivergenceError:
1086                      dt_new*=0.5
1087                      end_of_step=False
1088                      iteration_fail_counter+=1
1089                      if iteration_fail_counter>self.FAILED_TIME_STEPS_MAX:
1090                               raise SimulationBreakDownError("reduction of time step to achieve convergence failed after %s steps."%self.FAILED_TIME_STEPS_MAX)
1091                      self.trace("Iteration failed. Time step is repeated with new step size %s."%dt_new)
1092                   except FailedTimeStepError:
1093                      dt_new=self.getSafeTimeStepSize(dt)
1094                      end_of_step=False
1095                      step_fail_counter+=1
1096                      self.trace("Time step is repeated with new time step size %s."%dt_new)
1097                      if step_fail_counter>self.FAILED_TIME_STEPS_MAX:
1098                            raise SimulationBreakDownError("Time integration is given up after %d attempts."%step_fail_counter)
1099                dt=dt_new
1100                if not check_point==None:
1101                    if n%check_point==0:
1102                        self.trace("check point is created.")
1103                        self.writeXML()
1104            self.doFinalization()
1105    
        links the parameter name of model with the parameter name_for_object of  
        object.  
1106    
1107         The run method initiates checkpointing (it is not clear how to do this      def toDom(self, esysxml, node):
1108         yet)          """
1109     =====      C{toDom} method of Simulation class.
1110                    """
1111     """          simulation = esysxml.createElement('Simulation')
1112            esysxml.registerLinkableObject(self, simulation)
1113            for rank, sim in enumerate(self.iterModels()):
1114                component = esysxml.createElement('Component')
1115                component.setAttribute('rank', str(rank))
1116                sim.toDom(esysxml, component)
1117                simulation.appendChild(component)
1118            node.appendChild(simulation)
1119    
1120    
1121        def fromDom(cls, esysxml, node):
1122            sims = []
1123            for n in node.childNodes:
1124                if isinstance(n, minidom.Text):
1125                    continue
1126                sims.append(esysxml.getComponent(n))
1127            sims.sort(cmp=_comp)
1128            sim=cls([s[1] for s in sims], debug=esysxml.debug)
1129            esysxml.registerLinkableObject(sim, node)
1130            return sim
1131    
1132  if __name__=="__main__":      fromDom = classmethod(fromDom)
     import math  
     #  
     # test for parameter set  
     #  
     p11=ParameterSet()  
     p11.declareParameter(gamma1=1.,gamma2=2.,gamma3=3.)  
     p1=ParameterSet()  
     p1.declareParameter(dim=2,tol_v=0.001,output_file="/tmp/u.%3.3d.dx",runFlag=True,parm11=p11)  
     parm=ParameterSet({ "parm1" : p1 , "parm2" : ParameterSet(["alpha"])})  
     parm.parm2.alpha=Link(p11,"gamma1")  
     parm.x="that should not be here!"  
     print parm.showParameters()  
     # should be something like: {"parm2" : {"alpha" : reference to attribute  
     # gamma1 of <__main__.ParameterSet instance at 0xf6db51cc>},"parm1" : {"dim"  
     # : 2,"runFlag" : True,"tol_v": 0.001,"parm11" : {"gamma3" : 3.0,"gamma2" :  
     # 2.0,"gamma1" : 1.0},"output_file" : /tmp/u.%3.3d.dx}}  
     assert parm.parm2.alpha==1.  
     parm.writeXML()  
       
     #=======================  
     class Messenger(Model):  
         def __init__(self):  
             Model.__init__(self)  
             self.declareParameter(message="none")  
   
         def doInitialization(self,t):  
             self.__t=t  
             print "I start talking now!"  
   
         def doStep(self,dt):  
             self.__t+=dt  
             print "Message (time %e) : %s "%(self.__t,self.message)  
1133    
1134          def doFinalization(self):  def _comp(a,b):
1135              print "I have no more to say!"      if a[0]<a[1]:
1136            return 1
1137      class ODETEST(Model):      elif a[0]>a[1]:
1138          """ implements a solver for the ODE        return -1
1139        else:
1140          return 0
1141    
1142                du/dt=a*u+f(t)  class IterationDivergenceError(Exception):
1143        """
1144        Exception which is thrown if there is no convergence of the iteration
1145        process at a time step.
1146    
1147             we use a implicit euler scheme :      But there is a chance that a smaller step could help to reach convergence.
1148        """
1149        pass
1150    
1151                u_n-u_{n-1}= dt*a u_n + st*f(t_n)  class FailedTimeStepError(Exception):
1152        """
1153        Exception which is thrown if the time step fails because of a step
1154        size that have been choosen to be too large.
1155        """
1156        pass
1157    
1158             to get u_n we run an iterative process  class SimulationBreakDownError(Exception):
1159        """
1160                 u_{n.k}=u_{n-1}+dt*(a u_{n.i-1} + f(t_n))      Exception which is thrown if the simulation does not manage to
1161        progress in time.
1162        """
1163        pass
1164    
1165    class NonPositiveStepSizeError(Exception):
1166        """
1167        Exception which is thrown if the step size is not positive.
1168        """
1169        pass
1170    
1171             input for this model are step size dt, end time tend and a value for  class MissingLink(Exception):
1172             a, f and initial value for u. we need also a tolerance tol for a      """
1173             stopping criterion.      Exception thrown when a link is missing
1174        """
1175        pass
1176    
1177          """  class DataSource(object):
1178        """
1179        Class for handling data sources, including local and remote files. This class is under development.
1180        """
1181    
1182          def __init__(self):      def __init__(self, uri="file.ext", fileformat="unknown"):
1183              Model.__init__(self,debug=True)          self.uri = uri
1184              self.declareParameter(tend=1.,dt=0.1,a=0.9,u=10.,f=0.,message="",tol=1.e-8)          self.fileformat = fileformat
1185    
1186        def toDom(self, esysxml, node):
1187            """
1188            C{toDom} method of DataSource. Creates a DataSource node and appends it to the
1189        current XML esysxml.
1190            """
1191            ds = esysxml.createElement('DataSource')
1192            ds.appendChild(esysxml.createDataNode('URI', self.uri))
1193            ds.appendChild(esysxml.createDataNode('FileFormat', self.fileformat))
1194            node.appendChild(ds)
1195    
1196        def fromDom(cls, esysxml, node):
1197            uri= str(node.getElementsByTagName("URI")[0].firstChild.nodeValue.strip())
1198            fileformat= str(node.getElementsByTagName("FileFormat")[0].firstChild.nodeValue.strip())
1199            ds = cls(uri, fileformat)
1200            return ds
1201    
1202          def doInitialization(self,t):      def getLocalFileName(self):
1203              self.__tn=t          return self.uri
             self.__iter=0  
             
         def doIterationInitialization(self,dt):  
             self.__iter=0  
             self.__u_last=self.u              
   
         def doIterationStep(self,dt):  
             self.__iter+=1  
             self.__u_old=self.u  
             self.u=self.__u_last+dt*(self.a*self.__u_old+self.f)  
           
         def terminate(self):  
             if self.__iter<1:  
                 return False  
             else:  
                 return abs(self.__u_old-self.u)<self.tol*abs(self.u)  
1204    
1205          def doIterationFinalization(self,dt):      fromDom = classmethod(fromDom)
1206              self.__tn+=dt      
1207              self.message="current error = %e"%abs(self.u-10.*math.exp((self.a-1.)*self.__tn))  # vim: expandtab shiftwidth=4:
   
         def getSafeTimeStepSize(self,dt):  
             return min(self.dt,1./(abs(self.a)+1.))  
   
         def finalize(self):  
             return self.__tn>=self.tend  
   
     #  
     #   s solves the coupled ODE:  
     #  
     #     du/dt=a*u+  v  
     #     dv/dt=  u+a*v  
     #    
     #   each equation is treated through the ODETEST class. The equations are  
     #   linked and iteration over each time step is performed.  the current  
     #   error of v is reported by the Messenger class.  
     #  
     o1=ODETEST()  
     o1.u=10  
     o2=ODETEST()  
     o2.u=-10.  
     o1.f=Link(o2,"u")  
     o2.f=Link(o1,"u")  
     m=Messenger()  
     o1.dt=0.01  
     m.message=Link(o1)  
     s=ExplicitSimulation([Simulation([o1,o2],debug=True),m],debug=True)  
     s.run()  
     s.writeXML()  

Legend:
Removed from v.142  
changed lines
  Added in v.939

  ViewVC Help
Powered by ViewVC 1.1.26