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

  ViewVC Help
Powered by ViewVC 1.1.26