1 |
# $Id$ |
2 |
|
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: |
86 |
""" |
87 |
a Link makes an attribute of an object callable: |
88 |
o.object() |
89 |
o.a=8 |
90 |
l=Link(o,"a") |
91 |
assert l()==8 |
92 |
""" |
93 |
|
94 |
def __init__(self,target,attribute=None): |
95 |
""" |
96 |
creates a link to the object target. If attribute is given, the link is |
97 |
establised to this attribute of the target. Otherwise the attribute is |
98 |
undefined. |
99 |
""" |
100 |
self.target = target |
101 |
self.attribute = None |
102 |
self.setAttributeName(attribute) |
103 |
|
104 |
def setAttributeName(self,attribute): |
105 |
""" |
106 |
set a new attribute name to be collected from the target object. The |
107 |
target object must have the attribute with name attribute. |
108 |
""" |
109 |
if attribute and self.target and not hasattr(self.target, attribute): |
110 |
raise AttributeError("%s: target %s has no attribute %s."%(self, self.target, attribute)) |
111 |
self.attribute = attribute |
112 |
|
113 |
def hasDefinedAttributeName(self): |
114 |
""" |
115 |
returns true if an attribute name is set |
116 |
""" |
117 |
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 |
number_sequence = itertools.count(100) |
195 |
|
196 |
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. |
749 |
|
750 |
m=Model({"message" : "none" }) |
751 |
|
752 |
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): |
756 |
def __init__(self): |
757 |
Model.__init__(self,parameters={"message" : "none" }) |
758 |
|
759 |
m=MyModel() |
760 |
|
761 |
There are various ways how model parameters can be changed: |
762 |
|
763 |
1) use object attributes: |
764 |
|
765 |
m.message="Hello World!" |
766 |
|
767 |
2) use setParamter method |
768 |
|
769 |
|
770 |
m.setParameters(message="Hello World!") |
771 |
|
772 |
3) or dictonaries |
773 |
|
774 |
d={ message : "Hello World!" } |
775 |
m.setParameters(**d) |
776 |
|
777 |
|
778 |
A model executed buy staring the run method of the model: |
779 |
|
780 |
m=Messenger() |
781 |
m.run() |
782 |
|
783 |
The run methods marches through time. It first calls the |
784 |
doInitialization() method of the Model to set up the process. In each |
785 |
time step the doStep() method is called to get from the current to the |
786 |
next time step. The step size is defined by calling the |
787 |
getSafeTimeStepSize() method. the time integration process is |
788 |
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 |
794 |
what ever the current value of its parameter message is: |
795 |
|
796 |
class Messenger(Model): |
797 |
def __init__(self): |
798 |
Model.__init__(self,parameters={"message" : "none" }) |
799 |
|
800 |
def doInitialization(self): |
801 |
print "I start talking now!" |
802 |
|
803 |
def doStep(self,t): |
804 |
print "Message (time %e) : %s "%(t,self.message) |
805 |
|
806 |
def doFinalization(self): |
807 |
print "I have no more to say!" |
808 |
|
809 |
If a instance of the Messenger class is run, it will print the |
810 |
initialization and finalization message only. This is because the |
811 |
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: |
815 |
|
816 |
u(t=0)=u0 |
817 |
u_t=a*u**2 for all 0<t<=ten |
818 |
|
819 |
exact solution is given by u(t)=1/(1/u0-a*t) |
820 |
|
821 |
class Ode1(Model): |
822 |
def __init__(self,**args): |
823 |
Model.__init__(self,parameters={"tend" : 1., "dt" : 0.0001 ,"a" : 0.1 ,"u" : 1. },name="test",debug=True) |
824 |
|
825 |
def doInitialization(self): |
826 |
self._tn=0 |
827 |
|
828 |
def doStep(self,t): |
829 |
self.u=self.u+(t-self._tn)*self.a*self.u**2 |
830 |
self._tn=t |
831 |
|
832 |
def doFinalization(self): |
833 |
print "all done final error = ",abs(self.u-1./(1./3.-self.a*self._tn)) |
834 |
|
835 |
def getSafeTimeStepSize(self): |
836 |
return self.dt |
837 |
|
838 |
def finalize(self): |
839 |
return self._tn>=self.tend |
840 |
|
841 |
In some cases at a given time step an iteration process has to be |
842 |
performed to get the state of the Model for the next time step. ` In |
843 |
this case the doStep() method is replaced by a sequance of methods which |
844 |
implements this iterative process. The method then will control the |
845 |
iteration process by initializing the iteration through calling the |
846 |
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 |
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 |
857 |
|
858 |
has to be solved for u_{n+1}. The Newton scheme is used to solve this non-linear problem. |
859 |
|
860 |
|
861 |
class Ode2(Model): |
862 |
|
863 |
def __init__(self,**args): |
864 |
Model.__init__(self,{"tend" : 1., "dt" : 0.1 ,"a" : 10. ,"u" : 1. , "tol " : 1.e-8},"test","bla",None,True) |
865 |
|
866 |
def doInitialization(self): |
867 |
self.__tn=0 |
868 |
|
869 |
def doIterationInitialization(self,t): |
870 |
self.__iter=0 |
871 |
self.u_last=self.u |
872 |
self.current_dt=t-self.tn |
873 |
self.__tn=t |
874 |
|
875 |
def doIterationStep(self): |
876 |
self.__iter+=1 |
877 |
self.u_old=self.u |
878 |
self.u=(self.current_dt*self.a*self.u**2-self.u_last)/(2*self.current_dt*self.a*self.u-1.) |
879 |
|
880 |
def terminate(self): |
881 |
return abs(self.u_old-self.u)<self.tol*abs(self.u) |
882 |
|
883 |
def doIterationFinalization(self) |
884 |
print "all done" |
885 |
|
886 |
def getSafeTimeStepSize(self): |
887 |
return self.dt |
888 |
|
889 |
def finalize(self): |
890 |
return self.__tn>self.tend |
891 |
|
892 |
A model can be composed from subsimulations. Subsimulations are treated |
893 |
as model parameters. If a model parameter is set or a value of a model |
894 |
parameter is requested, the model will search for this parameter its |
895 |
subsimulations in the case the model does not have this parameter |
896 |
itself. The order in which the subsimulations are searched is critical. |
897 |
By default a Model initializes all its subsimulations, is finalized when |
898 |
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 |
907 |
|
908 |
o=Ode2() |
909 |
m=Messenger() |
910 |
om=Model(subsimulations=[o,m],debug=True) |
911 |
om.dt=0.01 |
912 |
om.u=1. |
913 |
m.message="it's me!" |
914 |
om.run() |
915 |
|
916 |
Notice that dt and u are parameters of class Ode2 and message is a |
917 |
parameter of the Messenger class. The Model formed from these models |
918 |
automatically hand the assignment of new values down to the |
919 |
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, |
925 |
typically an parameter of another Model object. |
926 |
|
927 |
|
928 |
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: |
930 |
|
931 |
s=Model() |
932 |
s.run() |
933 |
|
934 |
The run has an initializion and finalization phase. The latter is called |
935 |
if all subsimulations are to be finalized. The simulation is processing |
936 |
in time through calling the stepForward methods which updates the |
937 |
observables of each subsimulation. A time steps size which is save for |
938 |
all subsimulation is choosen. |
939 |
|
940 |
At given time step an iterative process may be performed to make sure |
941 |
that all observables are consistent across all subsimulations. In this |
942 |
case, similar the time dependence, an initialization and finalization of |
943 |
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 |
951 |
|
952 |
|
953 |
At any time the current value of the parameter name can be obtained by |
954 |
|
955 |
value=model.name |
956 |
|
957 |
If the object that has been assigned to the paramter/attribute name has |
958 |
the attribute/parameter name isself the current value of this attribute |
959 |
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) |
966 |
|
967 |
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 |
971 |
yet) |
972 |
===== |
973 |
|
974 |
""" |
975 |
|
976 |
|
977 |
|
978 |
if __name__=="__main__": |
979 |
import math |
980 |
# |
981 |
# test for parameter set |
982 |
# |
983 |
p11=ParameterSet() |
984 |
p11.declareParameter(gamma1=1.,gamma2=2.,gamma3=3.) |
985 |
p1=ParameterSet() |
986 |
p1.declareParameter(dim=2,tol_v=0.001,output_file="/tmp/u.%3.3d.dx",runFlag=True,parm11=p11) |
987 |
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): |
1013 |
print "I have no more to say!" |
1014 |
|
1015 |
class ODETEST(Model): |
1016 |
""" implements a solver for the ODE |
1017 |
|
1018 |
du/dt=a*u+f(t) |
1019 |
|
1020 |
we use a implicit euler scheme : |
1021 |
|
1022 |
u_n-u_{n-1}= dt*a u_n + st*f(t_n) |
1023 |
|
1024 |
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 |
|
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 |
""" |
1034 |
|
1035 |
def __init__(self): |
1036 |
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 doInitialization(self,t): |
1040 |
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 doIterationStep(self,dt): |
1048 |
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 |
def doIterationFinalization(self,dt): |
1059 |
self.__tn+=dt |
1060 |
self.message="current error = %e"%abs(self.u-10.*math.exp((self.a-1.)*self.__tn)) |
1061 |
|
1062 |
def getSafeTimeStepSize(self,dt): |
1063 |
return min(self.dt,1./(abs(self.a)+1.)) |
1064 |
|
1065 |
def finalize(self): |
1066 |
return self.__tn>=self.tend |
1067 |
|
1068 |
# |
1069 |
# s solves the coupled ODE: |
1070 |
# |
1071 |
# du/dt=a*u+ v |
1072 |
# dv/dt= u+a*v |
1073 |
# |
1074 |
# each equation is treated through the ODETEST class. The equations are |
1075 |
# linked and iteration over each time step is performed. the current |
1076 |
# error of v is reported by the Messenger class. |
1077 |
# |
1078 |
o1=ODETEST() |
1079 |
o1.u=10 |
1080 |
o2=ODETEST() |
1081 |
o2.u=-10. |
1082 |
o1.f=Link(o2,"u") |
1083 |
o2.f=Link(o1,"u") |
1084 |
m=Messenger() |
1085 |
o1.dt=0.01 |
1086 |
m.message=Link(o1) |
1087 |
s=ExplicitSimulation([Simulation([o1,o2],debug=True),m],debug=True) |
1088 |
s.run() |
1089 |
s.writeXML() |