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

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

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

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

Legend:
Removed from v.121  
changed lines
  Added in v.126

  ViewVC Help
Powered by ViewVC 1.1.26