/[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

revision 142 by jgs, Mon Jul 25 05:28:20 2005 UTC revision 147 by jgs, Fri Aug 12 01:45:47 2005 UTC
# Line 39  def all(seq): Line 39  def all(seq):
39              return False              return False
40      return True      return True
41    
42    def any(seq):
43        for x in seq:
44            if x:
45                return True
46        return False
47    
48  LinkableObjectRegistry = {}  LinkableObjectRegistry = {}
49    
50  def registerLinkableObject(obj_id, o):  def registerLinkableObject(obj_id, o):
# Line 50  def registerLink(obj_id, l): Line 56  def registerLink(obj_id, l):
56      LinkRegistry.append((obj_id,l))      LinkRegistry.append((obj_id,l))
57    
58  def parse(xml):  def parse(xml):
59        """
60        Generic parse method for EsysXML. Without this, Links don't work.
61        """
62      global LinkRegistry, LinkableObjectRegistry      global LinkRegistry, LinkableObjectRegistry
63      LinkRegistry = []      LinkRegistry = []
64      LinkableObjectRegistry = {}      LinkableObjectRegistry = {}
65    
66      doc = minidom.parseString(xml)      doc = minidom.parseString(xml)
   
67      sim = getComponent(doc.firstChild)      sim = getComponent(doc.firstChild)
68      for obj_id, link in LinkRegistry:      for obj_id, link in LinkRegistry:
69          link.target = LinkableObjectRegistry[obj_id]          link.target = LinkableObjectRegistry[obj_id]
# Line 63  def parse(xml): Line 71  def parse(xml):
71      return sim      return sim
72    
73  def getComponent(doc):  def getComponent(doc):
74        """
75        Used to get components of Simualtions, Models.
76        """
77      for node in doc.childNodes:      for node in doc.childNodes:
78            
79          if isinstance(node, minidom.Element):          if isinstance(node, minidom.Element):
80              if node.tagName == 'Simulation':              if node.tagName == 'Simulation':
81                  if node.getAttribute("type") == 'Simulation':                  if node.getAttribute("type") == 'Simulation':
82                      return Simulation.fromDom(node)                      return Simulation.fromDom(node)
                 elif node.getAttribute("type") == 'ExplicitSimulation':  
                     return ExplicitSimulation.fromDom(node)  
83              if node.tagName == 'Model':              if node.tagName == 'Model':
84                  model_type = node.getAttribute("type")                  model_type = node.getAttribute("type")
85                  model_subclasses = Model.__subclasses__()                  model_subclasses = Model.__subclasses__()
86                  for model in model_subclasses:                  for model in model_subclasses:
87                      if model_type == model.__name__:                      if model_type == model.__name__:
88                          return model.fromDom(node)                          return Model.fromDom(node)
89                                    if node.tagName == 'ParameterSet':
90                    parameter_type = node.getAttribute("type")
91                    return ParameterSet.fromDom(node)
92              raise "Invalid simulation type, %r" % node.getAttribute("type")              raise "Invalid simulation type, %r" % node.getAttribute("type")
93            
94    
95      raise ValueError("No Simulation Found")      raise ValueError("No Simulation Found")
96                            
# Line 106  class Link: Line 119  class Link:
119          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
120          target object must have the attribute with name attribute.          target object must have the attribute with name attribute.
121          """          """
122          if attribute and self.target and not hasattr(self.target, attribute):          if attribute and self.target:
123              raise AttributeError("%s: target %s has no attribute %s."%(self, self.target, attribute))              if isinstance(self.target,LinkableObject):
124                   if not self.target.hasAttribute(attribute):
125                      raise AttributeError("%s: target %s has no attribute %s."%(self, self.target, attribute))
126                else:
127                   if not hasattr(self.target,attribute):
128                      raise AttributeError("%s: target %s has no attribute %s."%(self, self.target, attribute))
129          self.attribute = attribute          self.attribute = attribute
130            
131      def hasDefinedAttributeName(self):      def hasDefinedAttributeName(self):
# Line 146  class Link: Line 164  class Link:
164          document          document
165          """          """
166          link = document.createElement('Link')          link = document.createElement('Link')
167            assert (self.target != None), ("Target was none, name was %r" % self.attribute)
168          link.appendChild(dataNode(document, 'Target', self.target.id))          link.appendChild(dataNode(document, 'Target', self.target.id))
169          # 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
170          # target from the xml later. I need a better unique identifier.          # target from the xml later. I need a better unique identifier.
# Line 168  class Link: Line 187  class Link:
187          If ostream is nor present the standart output stream is used.  If          If ostream is nor present the standart output stream is used.  If
188          esysheader==True the esys XML header is written          esysheader==True the esys XML header is written
189          """          """
190            print 'I got to the Link writeXML method'
191          document, rootnode = esysDoc()          document, rootnode = esysDoc()
192          self.toDom(document, rootnode)          self.toDom(document, rootnode)
193    
# Line 203  class LinkableObject(object): Line 222  class LinkableObject(object):
222          """ If debugging is on, print the message, otherwise do nothing          """ If debugging is on, print the message, otherwise do nothing
223          """          """
224          if self.debug:          if self.debug:
225              print msg              print "%s: %s"%(str(self),msg)
226            
227      def __getattr__(self,name):      def __getattr__(self,name):
228          """returns the value of attribute name. If the value is a Link object the          """returns the value of attribute name. If the value is a Link object the
# Line 223  class LinkableObject(object): Line 242  class LinkableObject(object):
242          if self.__linked_attributes.has_key(name):          if self.__linked_attributes.has_key(name):
243              return self.__linked_attributes[name]              return self.__linked_attributes[name]
244    
245            if self.__class__.__dict__.has_key(name):
246                return self.__class.__dict__[name]
247    
248          raise AttributeError,"No attribute %s."%name          raise AttributeError,"No attribute %s."%name
249            
250        def hasAttribute(self,name):
251            """returns True if self as attribute name"""
252            return self.__dict__.has_key(name) or self.__linked_attributes.has_key(name) or  self.__class__.__dict__.has_key(name)
253    
254      def __setattr__(self,name,value):      def __setattr__(self,name,value):
255          """sets the value for attribute name. If value is a Link the target          """sets the value for attribute name. If value is a Link the target
256          attribute is set to name if no attribute has been specified."""          attribute is set to name if no attribute has been specified."""
# Line 238  class LinkableObject(object): Line 264  class LinkableObject(object):
264                  value.setAttributeName(name)                  value.setAttributeName(name)
265              self.__linked_attributes[name] = value              self.__linked_attributes[name] = value
266    
267              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))
268          else:          else:
269              self.__dict__[name] = value              self.__dict__[name] = value
270            
# Line 252  class LinkableObject(object): Line 278  class LinkableObject(object):
278          else:          else:
279              raise AttributeError,"No attribute %s."%name              raise AttributeError,"No attribute %s."%name
280    
 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)  
   
281  class _ParameterIterator:  class _ParameterIterator:
282      def __init__(self,parameterset):      def __init__(self,parameterset):
283    
# Line 561  class ParameterSet(LinkableObject): Line 329  class ParameterSet(LinkableObject):
329          LinkableObject.__init__(self, **kwargs)          LinkableObject.__init__(self, **kwargs)
330          self.parameters = set()          self.parameters = set()
331          self.declareParameters(parameters)          self.declareParameters(parameters)
332    
333        def __repr__(self):
334            return "<%s %r>" % (self.__class__.__name__,
335                                [(p, getattr(self, p, None)) for p in self.parameters])
336            
337      def declareParameter(self,**parameters):      def declareParameter(self,**parameters):
338          """declares a new parameter(s) and its (their) inital value."""          """declares a new parameter(s) and its (their) inital value."""
# Line 577  class ParameterSet(LinkableObject): Line 349  class ParameterSet(LinkableObject):
349              setattr(self,prm,value)              setattr(self,prm,value)
350              self.parameters.add(prm)              self.parameters.add(prm)
351    
352              self.trace("%s: parameter %s has been declared."%(self,prm))              self.trace("parameter %s has been declared."%prm)
353    
354      def releaseParameters(self,name):      def releaseParameters(self,name):
355          """removes parameter name from the paramameters"""          """removes parameter name from the paramameters"""
356          if self.isParameter(name):          if self.isParameter(name):
357              self.parameters.remove(name)              self.parameters.remove(name)
358              self.debug("%s: parameter %s has been removed."%(self, name))              self.trace("parameter %s has been removed."%name)
359            
360      def __iter__(self):      def __iter__(self):
361          """creates an iterator over the parameter and their values"""          """creates an iterator over the parameter and their values"""
# Line 639  class ParameterSet(LinkableObject): Line 411  class ParameterSet(LinkableObject):
411    
412              node.appendChild(param)              node.appendChild(param)
413    
     
414      def fromDom(cls, doc):      def fromDom(cls, doc):
415    
416          # Define a host of helper functions to assist us.          # Define a host of helper functions to assist us.
# Line 673  class ParameterSet(LinkableObject): Line 444  class ParameterSet(LinkableObject):
444                      "bool":_boolfromValue                      "bool":_boolfromValue
445                      }                      }
446    
447    #        print doc.toxml()
448    
449          parameters = {}          parameters = {}
450          for node in _children(doc):          for node in _children(doc):
451              ptype = node.getAttribute("type")              ptype = node.getAttribute("type")
# Line 704  class ParameterSet(LinkableObject): Line 477  class ParameterSet(LinkableObject):
477          self.toDom(document, node)          self.toDom(document, node)
478          ostream.write(document.toprettyxml())          ostream.write(document.toprettyxml())
479    
480  class Model(ParameterSet,SimulationFrame):  class Model(ParameterSet):
481      """a Model is a SimulationFrame which is also a ParameterSet."""      """
482    
483        A Model object represents a processess marching over time
484        until a finalizing condition is fullfilled. At each time step an iterative
485        process can be performed and the time step size can be controlled. A Model has
486        the following work flow:
487    
488              doInitialization()
489              while not finalize():
490                   dt=getSafeTimeStepSize(dt)
491                   doStepPreprocessing(dt)
492                   while not terminateIteration(): doStep(dt)
493                   doStepPostprocessing(dt)
494              doFinalization()
495    
496              where doInitialization,finalize, getSafeTimeStepSize, doStepPreprocessing, terminateIteration, doStepPostprocessing, doFinalization
497              are methods of the particular instance of a Model. The default implementations of these methods have to be overwritten by
498              the subclass implementinf a Model.
499    
500        """
501    
502        UNDEF_DT=1.e300
503    
504      def __init__(self,parameters=[],**kwargs):      def __init__(self,parameters=[],**kwarg):
505          """creates a model"""          """creates a model
506          ParameterSet.__init__(self, parameters=parameters)  
507          SimulationFrame.__init__(self,**kwargs)              Just calls the parent constructor.
508            """
509            ParameterSet.__init__(self, parameters=parameters,**kwarg)
510    
511        def __str__(self):
512           return "<%s %d>"%(self.__class__,id(self))
513    
514      def toDom(self, document, node):      def toDom(self, document, node):
515          """ toDom method of Model class """          """ toDom method of Model class """
# Line 719  class Model(ParameterSet,SimulationFrame Line 518  class Model(ParameterSet,SimulationFrame
518          node.appendChild(pset)          node.appendChild(pset)
519          self._parametersToDom(document, pset)          self._parametersToDom(document, pset)
520    
521        def doInitialization(self):
522  class IterationDivergenceError(Exception):          """initializes the time stepping scheme. This function may be overwritten."""
523      """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          pass
524         to reach convergence."""      
525      pass      def getSafeTimeStepSize(self,dt):
526            """returns a time step size which can safely be used.
527  class IterationBreakDownError(Exception):             dt gives the previously used step size.
528      """excpetion which should be thrown if there is no conevregence and there is no chance that a time step reduction would help"""             This function may be overwritten."""
529      pass          return self.UNDEF_DT
530        
531  class FailedTimeStepError(Exception):      def finalize(self):
532      """excpetion which should be thrown if the time step fails because of a step size that have been choosen to be to large"""          """returns False if the time stepping is finalized. This function may be
533      pass          overwritten."""
534            return False
 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:  
   
         class Messenger(Model):  
             def __init__(self):  
                Model.__init__(self,parameters={"message" : "none" })  
   
         m=MyModel()  
   
        There are various ways how model parameters can be changed:  
   
        1) use object attributes:  
   
           m.message="Hello World!"  
   
        2) use setParamter method  
           
             
           m.setParameters(message="Hello World!")  
   
        3) or dictonaries  
     
            d={ message : "Hello World!" }  
            m.setParameters(**d)  
   
   
        A model executed buy staring the run method of the model:  
   
           m=Messenger()  
           m.run()  
   
        The run methods marches through time. It first calls the  
        doInitialization() method of the Model to set up the process. In each  
        time step the doStep() method is called to get from the current to the  
        next time step. The step size is defined by calling the  
        getSafeTimeStepSize() method.  the time integration process is  
        terminated when the finalize() methods return true. Final the  
        doFinalization() method is called to finalize the process. To implement  
        a particular model a subclass of the Model class is defined. The  
        subclass overwrites the default methods of Model.  
   
        The following class defines a messenger printing in the doStep method  
        what ever the current value of its parameter message is:  
   
        class Messenger(Model):  
             def __init__(self):  
                Model.__init__(self,parameters={"message" : "none" })  
   
             def doInitialization(self):  
                print "I start talking now!"  
   
             def doStep(self,t):  
                print "Message (time %e) : %s "%(t,self.message)  
   
             def doFinalization(self):  
                print "I have no more to say!"  
535                
536         If a instance of the Messenger class is run, it will print the      def doFinalization(self):
537         initialization and finalization message only.  This is because the          """finalizes the time stepping. This function may be overwritten."""
538         default method for finalize() does always returns True and therefore the          pass
539         transition is terminated startcht away.      
540                def doStepPreprocessing(self,dt):
541         Following example for solving the ODE using a forward euler scheme:          """sets up a time step of step size dt. This function may be overwritten."""
542            pass
543                  u(t=0)=u0      
544                  u_t=a*u**2       for all 0<t<=ten      def doStep(self,dt):
545            """executes an iteration step at a time step.
546        exact solution is given by u(t)=1/(1/u0-a*t)             dt is the currently used time step size.
547               This function may be overwritten."""
548            pass
549        
550        def terminateIteration(self):
551            """returns True if iteration on a time step is terminated."""
552            return True
553          
554        def doStepPostprocessing(self,dt):
555            """finalalizes the time step.
556               dt is the currently used time step size.
557               This function may be overwritten."""
558            pass
559        
560        def writeXML(self, ostream=stdout):
561            document, node = esysDoc()
562            self.toDom(document, node)
563            ostream.write(document.toprettyxml())
564        
565        
566    
567        class  Ode1(Model):  class Simulation(Model):
568           def __init__(self,**args):      """
569              Model.__init__(self,parameters={"tend" : 1., "dt" : 0.0001 ,"a" : 0.1 ,"u" : 1. },name="test",debug=True)            A Simulation object is special Model which runs a sequence of Models.
   
          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.  
570    
571         following example is a modification of the example above. Here an            The methods doInitialization,finalize, getSafeTimeStepSize, doStepPreprocessing,
572         implicit euler scheme is used. in each time step the problem            terminateIteration, doStepPostprocessing, doFinalization
573              are executing the corresponding methods of the models in the simulation.
574                        
575             0= u_{n+1}-u_{n}+a*dt*u_{n+1}**2      """
576        
577         has to be solved for u_{n+1}. The Newton scheme is used to solve this non-linear problem.      FAILED_TIME_STEPS_MAX=20
578        MAX_ITER_STEPS=50
579        
580        class  Ode2(Model):      def __init__(self, models=[], **kwargs):
581            """initiates a simulation from a list of models. """
582         def __init__(self,**args):          Model.__init__(self, **kwargs)
583             Model.__init__(self,{"tend" : 1., "dt" : 0.1 ,"a" : 10. ,"u" : 1. , "tol " : 1.e-8},"test","bla",None,True)          self.__models=[]
   
        def doInitialization(self):  
            self.__tn=0  
   
        def doIterationInitialization(self,t):  
             self.__iter=0  
             self.u_last=self.u              
             self.current_dt=t-self.tn  
             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.  
   
        A parameter of a Model can be linked to an attribute of onother object,  
        typically an parameter of another Model object.      
         
         
        which is comprised by a set of subsimulations.  
        The simulation is run through its run method which in the simplest case has the form:  
   
           s=Model()  
           s.run()  
584    
585         The run has an initializion and finalization phase. The latter is called          for i in range(len(models)):
586         if all subsimulations are to be finalized. The simulation is processing              self[i] = models[i]
        in time through calling the stepForward methods which updates the  
        observables of each subsimulation.  A time steps size which is save for  
        all subsimulation is choosen.  
587    
588         At given time step an iterative process may be performed to make sure      def __repr__(self):
589         that all observables are consistent across all subsimulations.  In this          """
590         case, similar the time dependence, an initialization and finalization of          returns a string representation of the Simulation
591         the iteration is performed.          """
592            return "<Simulation %r>" % self.__models
593    
594         A Model has input and output parameters where each input parameter can      def __str__(self):
595         be constant, time dependent or may depend on an output parameter of          """
596         another model or the model itself. To create a parameter name of a model          returning Simulation as a string
597         and to assign a value to it one can use the statement          """
598            return "<Simulation %d>"%id(self)
599        
600        def iterModels(self):
601            """returns an iterator over the models"""
602            return self.__models
603        
604        def __getitem__(self,i):
605            """returns the i-th model"""
606            return self.__models[i]
607        
608        def __setitem__(self,i,value):
609            """sets the i-th model"""
610            if not isinstance(value,Model):
611                raise ValueError("assigned value is not a Model")
612            for j in range(max(i-len(self.__models)+1,0)): self.__models.append(None)
613            self.__models[i]=value
614        
615        def __len__(self):
616            """returns the number of models"""
617            return len(self.__models)
618    
619             model.name=object      def toDom(self, document, node):
620            """ toDom method of Simulation class """
621            simulation = document.createElement('Simulation')
622            simulation.setAttribute('type', self.__class__.__name__)
623    
624            for rank, sim in enumerate(self.iterModels()):
625                component = document.createElement('Component')
626                component.setAttribute('rank', str(rank))
627    
628         At any time the current value of the parameter name can be obtained by              sim.toDom(document, component)
629    
630                 value=model.name              simulation.appendChild(component)
631    
632         If the object that has been assigned to the paramter/attribute name has          node.appendChild(simulation)
        the attribute/parameter name isself the current value of this attribute  
        of the object is returned (e.g. for model.name=object where object has  
        an attribute name, the statement value=model.name whould assign the  
        value object.name to value.). If the name of the parameters of a model  
        and an object don't match the setParameter method of model can be used.  
        So  
633    
634             model.setParameter(name,object,name_for_object)      def writeXML(self,ostream=stdout):
635            """writes the object as an XML object into an output stream"""
636            document, rootnode = esysDoc()
637            self.toDom(document, rootnode)
638            ostream.write(document.toprettyxml())
639        
640        def getSafeTimeStepSize(self,dt):
641            """returns a time step size which can safely be used by all models.
642               This is the minimum over the time step sizes of all models."""
643            out=min([o.getSafeTimeStepSize(dt) for o in self.iterModels()])
644            print "%s: safe step size is %e."%(str(self),out)
645            return out
646        
647        def doInitialization(self):
648            """initializes all models """
649            self.n=0
650            self.tn=0.
651            for o in self.iterModels():
652                o.doInitialization()
653        
654        def finalize(self):
655            """returns True if any of the models is to be finalized"""
656            return any([o.finalize() for o in self.iterModels()])
657          
658        def doFinalization(self):
659            """finalalizes the time stepping for all models."""
660            for i in self.iterModels(): i.doFinalization()
661            self.trace("end of time integation.")
662        
663        def doStepPreprocessing(self,dt):
664            """initializes the time step for all models."""
665            for o in self.iterModels():
666                o.doStepPreprocessing(dt)
667        
668        def terminateIteration(self):
669            """returns True if all iterations for all models are terminated."""
670            out=all([o.terminateIteration() for o in self.iterModels()])
671            return out
672          
673        def doStepPostprocessing(self,dt):
674            """finalalizes the iteration process for all models."""
675            for o in self.iterModels():
676                o.doStepPostprocessing(dt)
677            self.n+=1
678            self.tn+=dt
679        
680        def doStep(self,dt):
681            """
682                 executes the iteration step at a time step for all model:
683    
684                      self.doStepPreprocessing(dt)
685                      while not self.terminateIteration(): for all models: self.doStep(dt)
686                      self.doStepPostprocessing(dt)
687    
688         links the parameter name of model with the parameter name_for_object of          """
689         object.          self.iter=0
690            while not self.terminateIteration():
691                if self.iter==0: self.trace("iteration at %d-th time step %e starts"%(self.n+1,self.tn+dt))
692                self.iter+=1
693                self.trace("iteration step %d"%(self.iter))
694                for o in self.iterModels():
695                      o.doStep(dt)
696            if self.iter>0: self.trace("iteration at %d-th time step %e finalized."%(self.n+1,self.tn+dt))
697    
698         The run method initiates checkpointing (it is not clear how to do this      def run(self,check_point=None):
699         yet)          """
    =====  
               
    """  
700    
701               run the simulation by performing essentially
702        
703                   self.doInitialization()
704                   while not self.finalize():
705                      dt=self.getSafeTimeStepSize()
706                      self.doStep(dt)
707                      if n%check_point==0: self.writeXML()
708                   self.doFinalization()
709    
710               If one of the models in throws a FailedTimeStepError exception a new time step size is
711               computed through getSafeTimeStepSize() and the time step is repeated.
712      
713               If one of the models in throws a IterationDivergenceError exception the time step size
714               is halved and the time step is repeated.
715    
716  if __name__=="__main__":             In both cases the time integration is given up after Simulation.FAILED_TIME_STEPS_MAX
717      import math             attempts.
     #  
     # 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)  
718    
719          def doFinalization(self):      
720              print "I have no more to say!"          """
721              dt=self.UNDEF_DT
722      class ODETEST(Model):          self.doInitialization()
723          """ implements a solver for the ODE          while not self.finalize():
724                step_fail_counter=0
725                iteration_fail_counter=0
726                dt_new=self.getSafeTimeStepSize(dt)
727                self.trace("%d. time step %e (step size %e.)" % (self.n+1,self.tn+dt_new,dt_new))
728                end_of_step=False
729                while not end_of_step:
730                   end_of_step=True
731                   if not dt_new>0:
732                      raise NonPositiveStepSizeError("non-positive step size in step %d",self.n+1)
733                   try:
734                      self.doStepPreprocessing(dt_new)
735                      self.doStep(dt_new)
736                      self.doStepPostprocessing(dt_new)
737                   except IterationDivergenceError:
738                      dt_new*=0.5
739                      end_of_step=False
740                      iteration_fail_counter+=1
741                      if iteration_fail_counter>self.FAILED_TIME_STEPS_MAX:
742                               raise SimulationBreakDownError("reduction of time step to achieve convergence failed.")
743                      self.trace("iteration fails. time step is repeated with new step size.")
744                   except FailedTimeStepError:
745                      dt_new=self.getSafeTimeStepSize(dt)
746                      end_of_step=False
747                      step_fail_counter+=1
748                      self.trace("time step is repeated.")
749                      if step_fail_counter>self.FAILED_TIME_STEPS_MAX:
750                            raise SimulationBreakDownError("time integration is given up after %d attempts."%step_fail_counter)
751                dt=dt_new
752                if not check_point==None:
753                    if n%check_point==0:
754                        self.trace("check point is created.")
755                        self.writeXML()
756            self.doFinalization()
757    
758                du/dt=a*u+f(t)      def fromDom(cls, doc):
759            sims = []
760            for node in doc.childNodes:
761                if isinstance(node, minidom.Text):
762                    continue
763    
764             we use a implicit euler scheme :              sims.append(getComponent(node))
765    
766                u_n-u_{n-1}= dt*a u_n + st*f(t_n)          return cls(sims)
767    
768             to get u_n we run an iterative process      fromDom = classmethod(fromDom)
   
                u_{n.k}=u_{n-1}+dt*(a u_{n.i-1} + f(t_n))  
769    
770    
771             input for this model are step size dt, end time tend and a value for  class IterationDivergenceError(Exception):
772             a, f and initial value for u. we need also a tolerance tol for a      """
773             stopping criterion.         excpetion which is thrown if there is no convergence of the iteration process at a time step
774           but there is a chance taht a smaller step could help to reach convergence.
775    
776          """      """
777        pass
778    
779          def __init__(self):  class FailedTimeStepError(Exception):
780              Model.__init__(self,debug=True)      """excpetion which is thrown if the time step fails because of a step size that have been choosen to be too large"""
781              self.declareParameter(tend=1.,dt=0.1,a=0.9,u=10.,f=0.,message="",tol=1.e-8)      pass
782    
783          def doInitialization(self,t):  class SimulationBreakDownError(Exception):
784              self.__tn=t      """excpetion which is thrown if the simulation does not manage to progress in time."""
785              self.__iter=0      pass
             
         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)  
786    
787          def doIterationFinalization(self,dt):  class NonPositiveStepSizeError(Exception):
788              self.__tn+=dt      """excpetion which is thrown if the step size is not positive"""
789              self.message="current error = %e"%abs(self.u-10.*math.exp((self.a-1.)*self.__tn))      pass
   
         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.147

  ViewVC Help
Powered by ViewVC 1.1.26