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

  ViewVC Help
Powered by ViewVC 1.1.26