/[escript]/branches/3.4.1/modellib/py_src/input.py
ViewVC logotype

Annotation of /branches/3.4.1/modellib/py_src/input.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4587 - (hide annotations)
Wed Dec 11 06:17:09 2013 UTC (5 years, 4 months ago) by jfenwick
File MIME type: text/x-python
File size: 20630 byte(s)
Preparation begins

1 ksteube 1809
2 jfenwick 3981 ##############################################################################
3 ksteube 1312 #
4 jfenwick 4154 # Copyright (c) 2003-2013 by University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 ksteube 1312 #
7 ksteube 1809 # Primary Business: Queensland, Australia
8     # Licensed under the Open Software License version 3.0
9     # http://www.opensource.org/licenses/osl-3.0.php
10 ksteube 1312 #
11 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     # Development since 2012 by School of Earth Sciences
13     #
14     ##############################################################################
15 elspeth 628
16 jfenwick 4154 __copyright__="""Copyright (c) 2003-2013 by University of Queensland
17 jfenwick 3981 http://www.uq.edu.au
18 ksteube 1809 Primary Business: Queensland, Australia"""
19 elspeth 628 __license__="""Licensed under the Open Software License version 3.0
20 ksteube 1809 http://www.opensource.org/licenses/osl-3.0.php"""
21 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
22 elspeth 628
23 jfenwick 3432 from esys.escript import length, wherePositive, whereNegative, exp, inf, sup
24 jgs 149 from esys.escript.modelframe import Model,ParameterSet
25 gross 927 from esys.escript.linearPDEs import LinearPDE
26 jgs 147 from math import log
27 jfenwick 3432 import numpy
28 jgs 127
29 jgs 147 class Sequencer(Model):
30 jgs 148 """
31 gross 814 Runs through time until t_end is reached.
32 gross 720
33 jfenwick 2625 :ivar t_end: model is terminated when t_end is passed, default 1 (in).
34     :type t_end: ``float``
35     :ivar dt_max: maximum time step size, default `Model.UNDEF_DT` (in)
36     :type dt_max: ``float``
37     :ivar t: current time stamp (in/out). By default it is initialized with zero.
38     :type t: ``float``
39 gross 814
40 jgs 148 """
41 gross 918 def __init__(self,**kwargs):
42 jgs 148 """
43 jgs 149 """
44 gross 918 super(Sequencer,self).__init__(**kwargs)
45 gross 423 self.declareParameter(t=0.,
46     t_end=1.,
47     dt_max=Model.UNDEF_DT)
48 jgs 127
49 jgs 147 def doInitialization(self):
50 jgs 148 """
51 jgs 149 initialize time integration
52 jgs 148 """
53     self.__t_old = self.t
54 jgs 147
55 jgs 148 def doStepPreprocessing(self, dt):
56     self.t = self.__t_old+dt
57 jgs 147
58 jgs 148 def doStepPostprocessing(self, dt):
59     self.__t_old = self.t
60 jgs 147
61     def finalize(self):
62 jgs 148 """
63 jfenwick 2625 returns true when `t` has reached `t_end`
64 jgs 148 """
65     return self.t >= self.t_end
66 jgs 147
67 jgs 148 def getSafeTimeStepSize(self, dt):
68     """
69 jfenwick 2625 returns `dt_max`
70 jgs 148 """
71     return self.dt_max
72 jgs 147
73 jgs 148 class GaussianProfile(ParameterSet):
74     """
75 jgs 149 Generates a Gaussian profile at center x_c, width width and height A
76 jgs 148 over a domain
77 jgs 127
78 jfenwick 2625 :ivar domain: domain
79     :ivar x_c: center of the Gaussian profile (default [0.,0.,0.])
80     :ivar A: (in) height of the profile. A maybe a vector. (default 1.)
81     :ivar width: (in) width of the profile (default 0.1)
82     :ivar r: (in) radius of the circle (default = 0)
83 jgs 127
84 jgs 148 In the case that the spatial dimension is two, The third component of
85 gross 814 x_c is dropped.
86 jgs 127 """
87 gross 918 def __init__(self,**kwargs):
88     super(GaussianProfile, self).__init__(**kwargs)
89 jgs 147 self.declareParameter(domain=None,
90 jfenwick 2455 x_c=numpy.zeros([3]),
91 jgs 147 A=1.,
92     width=0.1,
93     r=0)
94 jgs 127
95     def out(self):
96 jgs 149 """
97     Generate the Gaussian profile
98 gross 814
99     Link against this method to get the output of this model.
100 jgs 149 """
101 jgs 148 x = self.domain.getX()
102     dim = self.domain.getDim()
103     l = length(x-self.x_c[:dim])
104 gross 323 m = whereNegative(l-self.r)
105 jgs 148
106 jgs 127 return (m+(1.-m)*exp(-log(2.)*(l/self.width)**2))*self.A
107    
108 jgs 148 class InterpolateOverBox(ParameterSet):
109 jgs 147 """
110 jgs 148 Returns values at each time. The values are defined through given values
111 gross 814 at time node. For two dimensional domains back values are ignored.
112 jgs 147
113 jfenwick 2625 :ivar domain: domain
114     :ivar value_left_bottom_front: (in) value at left,bottom,front corner
115     :ivar value_right_bottom_front: (in) value at right, bottom, front corner
116     :ivar value_left_top_front: (in) value at left,top,front corner
117     :ivar value_right_top_front: (in) value at right,top,front corner
118     :ivar value_left_bottom_back: (in) value at left,bottom,back corner
119     :ivar value_right_bottom_back: (in) value at right,bottom,back corner
120     :ivar value_left_top_back: (in) value at left,top,back corner
121     :ivar value_right_top_back: (in) value at right,top,back corner
122 jgs 147 """
123    
124 gross 918 def __init__(self, **kwargs):
125 gross 911 super(InterpolateOverBox, self).__init__(self)
126 jgs 147 self.declareParameter(domain=None,
127     value_left_bottom_front=0.,
128     value_right_bottom_front=0.,
129     value_left_top_front=0.,
130     value_right_top_front=0.,
131     value_left_bottom_back=0.,
132     value_right_bottom_back=0.,
133     value_left_top_back=0.,
134     value_right_top_back=0.)
135    
136    
137     def out(self):
138 gross 814 """
139     values at domain locations by bilinear interpolation of the given values.
140    
141     Link against this method to get the output of this model.
142     """
143 jgs 148 x = self.domain.getX()
144     if self.domain.getDim() == 2:
145 gross 819 x0,x1=x[0],x[1]
146     left_bottom_front0,right_top_back0=inf(x0),sup(x0)
147 gross 820 left_bottom_front1,right_top_back1=inf(x1),sup(x1)
148     f_right = (x0 - left_bottom_front0)/(right_top_back0 -left_bottom_front0)
149 jgs 148 f_left = 1. - f_right
150 gross 820 f_top = (x1 - left_bottom_front1)/(right_top_back1 - left_bottom_front1)
151 jgs 148 f_bottom = 1. - f_top
152 gross 816 out = f_left * f_bottom * self.value_left_bottom_front \
153     + f_right * f_bottom * self.value_right_bottom_front \
154     + f_left * f_top * self.value_left_top_front \
155     + f_right * f_top * self.value_right_top_front
156 jgs 147 else:
157 gross 819 x0,x1,x2=x[0],x[1],x[2]
158     left_bottom_front0,right_top_back0=inf(x0),sup(x0)
159 gross 820 left_bottom_front1,right_top_back1=inf(x1),sup(x1)
160     left_bottom_front2,right_top_back2=inf(x2),sup(x2)
161     f_right = (x0 - left_bottom_front0)/(right_top_back0 - left_bottom_front0)
162 jgs 148 f_left = 1. - f_right
163 gross 820 f_top = (x1 - left_bottom_front1)/(right_top_back1 - left_bottom_front1)
164 jgs 148 f_bottom = 1. - f_top
165 gross 820 f_back = (x2 - left_bottom_front1)/(right_top_back2 - left_bottom_front2)
166 jgs 148 f_front = 1. - f_back
167 gross 816 out = f_left * f_bottom * f_front * self.value_left_bottom_front\
168     + f_right * f_bottom * f_front * self.value_right_bottom_front\
169     + f_left * f_top * f_front * self.value_left_top_front\
170     + f_right * f_top * f_front * self.value_right_top_front\
171     + f_left * f_bottom * f_back * self.value_left_bottom_back\
172     + f_right * f_bottom * f_back * self.value_right_bottom_back\
173     + f_left * f_top * f_back * self.value_left_top_back\
174     + f_right * f_top * f_back * self.value_right_top_back
175 jgs 147 return out
176    
177    
178 jgs 148 class InterpolatedTimeProfile(ParameterSet):
179     """
180    
181     Returns values at each time. The values are defined through given
182     values at time node.
183 jgs 147
184 jgs 148 value[i] defines the value at time nodes[i]. Between nodes linear
185     interpolation is used.
186    
187     For time t<nodes[0], value[0] is used and for t>nodes[l], values[l]
188     is used where l=len(nodes)-1.
189 jgs 147
190 jfenwick 2625 :ivar t: (in) current time
191     :ivar node: (in) list of time nodes
192     :ivar values: (in) list of values at time nodes
193 jgs 148 """
194 jgs 127
195 gross 918 def __init__(self,**kwargs):
196     super( InterpolatedTimeProfile, self).__init__(**kwargs)
197 jgs 147 self.declareParameter(t=0., \
198     nodes=[0.,1.],\
199     values=[1.,1.])
200     def out(self):
201 gross 814 """
202     current value
203    
204     Link against this method to get the output of this model.
205     """
206 jgs 148 l = len(self.nodes) - 1
207     t = self.t
208     if t <= self.nodes[0]:
209     return self.values[0]
210 jgs 147 else:
211     for i in range(1,l):
212 jgs 148 if t < self.nodes[i]:
213     m = (self.values[i-1] - self.values[i])/\
214     (self.nodes[i-1] - self.nodes[i])
215     return m*(t-self.nodes[i-1]) + self.values[i-1]
216 jgs 147 return self.values[l]
217 jgs 127
218 gross 903 class ScalarDistributionFromTags(ParameterSet):
219 jgs 148 """
220 gross 953 creates a scalar distribution on a domain from tags, If tag_map is given
221     the tags can be given a names and tag_map is used to map it into domain tags.
222 gross 903
223 jfenwick 2625 :ivar domain: domain
224     :type domain: `esys.escript.Domain`
225     :ivar default: default value
226     :ivar tag0: tag 0
227     :type tag0: ``int``
228     :ivar value0: value for tag 0
229     :type value0: ``float``
230     :ivar tag1: tag 1
231     :type tag1: ``int``
232     :ivar value1: value for tag 1
233     :type value1: ``float``
234     :ivar tag2: tag 2
235     :type tag2: ``int``
236     :ivar value2: value for tag 2
237     :type value2: ``float``
238     :ivar tag3: tag 3
239     :type tag3: ``int``
240     :ivar value3: value for tag 3
241     :type value3: ``float``
242     :ivar tag4: tag 4
243     :type tag4: ``int``
244     :ivar value4: value for tag 4
245     :type value4: ``float``
246     :ivar tag5: tag 5
247     :type tag5: ``int``
248     :ivar value5: value for tag 5
249     :type value5: ``float``
250     :ivar tag6: tag 6
251     :type tag6: ``int``
252     :ivar value6: value for tag 6
253     :type value6: ``float``
254     :ivar tag7: tag 7
255     :type tag7: ``int``
256     :ivar value7: value for tag 7
257     :type value7: ``float``
258     :ivar tag8: tag 8
259     :type tag8: ``int``
260     :ivar value8: value for tag 8
261     :type value8: ``float``
262     :ivar tag9: tag 9
263     :type tag9: ``int``
264     :ivar value9: value for tag 9
265     :type value9: ``float``
266 gross 903 """
267 gross 918 def __init__(self,**kwargs):
268     super(ScalarDistributionFromTags, self).__init__(**kwargs)
269 gross 903 self.declareParameter(domain=None,
270     default=0.,
271     tag0=None,
272     value0=0.,
273     tag1=None,
274     value1=0.,
275     tag2=None,
276     value2=0.,
277     tag3=None,
278     value3=0.,
279     tag4=None,
280     value4=0.,
281     tag5=None,
282     value5=0.,
283     tag6=None,
284     value6=0.,
285     tag7=None,
286     value7=0.,
287     tag8=None,
288     value8=0.,
289     tag9=None,
290     value9=0.)
291    
292    
293     def out(self):
294     """
295 jfenwick 2625 returns a `esys.escript.Data` object
296 gross 903 Link against this method to get the output of this model.
297     """
298     d=Scalar(self.default,Function(self.domain))
299 gross 1098 if not self.tag0 == None: d.setTaggedValue(self.tag0,self.value0)
300     if not self.tag1 == None: d.setTaggedValue(self.tag1,self.value1)
301     if not self.tag2 == None: d.setTaggedValue(self.tag2,self.value2)
302     if not self.tag3 == None: d.setTaggedValue(self.tag3,self.value3)
303     if not self.tag4 == None: d.setTaggedValue(self.tag4,self.value4)
304     if not self.tag5 == None: d.setTaggedValue(self.tag5,self.value5)
305     if not self.tag6 == None: d.setTaggedValue(self.tag6,self.value6)
306     if not self.tag7 == None: d.setTaggedValue(self.tag7,self.value7)
307     if not self.tag8 == None: d.setTaggedValue(self.tag8,self.value8)
308     if not self.tag9 == None: d.setTaggedValue(self.tag9,self.value9)
309 gross 903 return d
310    
311 gross 927 class SmoothScalarDistributionFromTags(ParameterSet):
312     """
313     creates a smooth scalar distribution on a domain from region tags
314    
315 jfenwick 2625 :ivar domain: domain
316     :type domain: `esys.escript.Domain`
317     :ivar default: default value
318     :ivar tag0: tag 0
319     :type tag0: ``int``
320     :ivar value0: value for tag 0
321     :type value0: ``float``
322     :ivar tag1: tag 1
323     :type tag1: ``int``
324     :ivar value1: value for tag 1
325     :type value1: ``float``
326     :ivar tag2: tag 2
327     :type tag2: ``int``
328     :ivar value2: value for tag 2
329     :type value2: ``float``
330     :ivar tag3: tag 3
331     :type tag3: ``int``
332     :ivar value3: value for tag 3
333     :type value3: ``float``
334     :ivar tag4: tag 4
335     :type tag4: ``int``
336     :ivar value4: value for tag 4
337     :type value4: ``float``
338     :ivar tag5: tag 5
339     :type tag5: ``int``
340     :ivar value5: value for tag 5
341     :type value5: ``float``
342     :ivar tag6: tag 6
343     :type tag6: ``int``
344     :ivar value6: value for tag 6
345     :type value6: ``float``
346     :ivar tag7: tag 7
347     :type tag7: ``int``
348     :ivar value7: value for tag 7
349     :type value7: ``float``
350     :ivar tag8: tag 8
351     :type tag8: ``int``
352     :ivar value8: value for tag 8
353     :type value8: ``float``
354     :ivar tag9: tag 9
355     :type tag9: ``int``
356     :ivar value9: value for tag 9
357     :type value9: ``float``
358 gross 927 """
359     def __init__(self,**kwargs):
360     super(SmoothScalarDistributionFromTags, self).__init__(**kwargs)
361     self.declareParameter(domain=None,
362     default=0.,
363     tag0=None,
364     value0=0.,
365     tag1=None,
366     value1=0.,
367     tag2=None,
368     value2=0.,
369     tag3=None,
370     value3=0.,
371     tag4=None,
372     value4=0.,
373     tag5=None,
374     value5=0.,
375     tag6=None,
376     value6=0.,
377     tag7=None,
378     value7=0.,
379     tag8=None,
380     value8=0.,
381     tag9=None,
382     value9=0.)
383    
384    
385     def __update(self,tag,tag_value,value):
386     if self.__pde==None:
387     self.__pde=LinearPDE(self.domain,numSolutions=1)
388     mask=Scalar(0.,Function(self.domain))
389 gross 1098 mask.setTaggedValue(tag,1.)
390 gross 927 self.__pde.setValue(Y=mask)
391     mask=wherePositive(abs(self.__pde.getRightHandSide()))
392     value*=(1.-mask)
393     value+=tag_value*mask
394     return value
395    
396     def out(self):
397     """
398 jfenwick 2625 returns a `esys.escript.Data` object
399 gross 927 Link against this method to get the output of this model.
400     """
401     d=Scalar(self.default,Solution(self.domain))
402     self.__pde=None
403     if not self.tag0 == None: d=self.__update(self.tag0,self.value0,d)
404     if not self.tag1 == None: d=self.__update(self.tag1,self.value1,d)
405     if not self.tag2 == None: d=self.__update(self.tag2,self.value2,d)
406     if not self.tag3 == None: d=self.__update(self.tag3,self.value3,d)
407     if not self.tag4 == None: d=self.__update(self.tag4,self.value4,d)
408     if not self.tag5 == None: d=self.__update(self.tag5,self.value5,d)
409     if not self.tag6 == None: d=self.__update(self.tag6,self.value6,d)
410     if not self.tag7 == None: d=self.__update(self.tag7,self.value7,d)
411     if not self.tag8 == None: d=self.__update(self.tag8,self.value8,d)
412     if not self.tag9 == None: d=self.__update(self.tag9,self.value9,d)
413     return d
414    
415 gross 903 class LinearCombination(ParameterSet):
416     """
417 jgs 148 Returns a linear combination of the f0*v0+f1*v1+f2*v2+f3*v3+f4*v4
418 jgs 147
419 jfenwick 2625 :ivar f0: numerical object or None, default=None (in)
420     :ivar v0: numerical object or None, default=None (in)
421     :ivar f1: numerical object or None, default=None (in)
422     :ivar v1: numerical object or None, default=None (in)
423     :ivar f2: numerical object or None, default=None (in)
424     :ivar v2: numerical object or None, default=None (in)
425     :ivar f3: numerical object or None, default=None (in)
426     :ivar v3: numerical object or None, default=None (in)
427     :ivar f4: numerical object or None, default=None (in)
428     :ivar v4: numerical object or None, default=None (in)
429 jgs 148 """
430 gross 918 def __init__(self,**kwargs):
431     super(LinearCombination, self).__init__(**kwargs)
432 jgs 148 self.declareParameter(f0=None, \
433     v0=None, \
434     f1=None, \
435     v1=None, \
436     f2=None, \
437     v2=None, \
438     f3=None, \
439     v3=None, \
440     f4=None, \
441     v4=None)
442 jgs 147
443 jgs 148 def out(self):
444 gross 814 """
445     returns f0*v0+f1*v1+f2*v2+f3*v3+f4*v4.
446     Link against this method to get the output of this model.
447     """
448 jgs 148 if not self.f0 == None and not self.v0 == None:
449     fv0 = self.f0*self.v0
450     else:
451     fv0 = None
452 jgs 147
453 jgs 148 if not self.f1 == None and not self.v1 == None:
454     fv1 = self.f1*self.v1
455     else:
456     fv1 = None
457 jgs 147
458 jgs 148 if not self.f2 == None and not self.v2 == None:
459     fv2 = f2*v2
460     else:
461     fv2 = None
462 jgs 147
463 jgs 148 if not self.f3 == None and not self.v3 == None:
464     fv3 = self.f3*self.v3
465     else:
466     fv3 = None
467 jgs 147
468 jgs 148 if not self.f4 == None and not self.v4 == None:
469     fv4 = self.f4*self.v4
470     else:
471     fv4 = None
472 jgs 147
473 jgs 148 if fv0 == None:
474     out = 0.
475     else:
476     out = fv0
477     if not fv1 == None:
478     out += fv1
479     if not fv2 == None:
480     out += fv2
481     if not fv3 == None:
482     out += fv3
483     return out
484 jgs 147
485 gross 904 class MergeConstraints(ParameterSet):
486     """
487     Returns a linear combination of the f0*v0+f1*v1+f2*v2+f3*v3+f4*v4
488     """
489 gross 918 def __init__(self,**kwargs):
490     super(MergeConstraints, self).__init__(**kwargs)
491 gross 904 self.declareParameter(location_of_constraint0=None, \
492     value_of_constraint0=None, \
493     location_of_constraint1=None, \
494     value_of_constraint1=None, \
495     location_of_constraint2=None, \
496     value_of_constraint2=None, \
497     location_of_constraint3=None, \
498     value_of_constraint3=None, \
499     location_of_constraint4=None, \
500     value_of_constraint4=None)
501     def location_of_constraint(self):
502     """
503     return the values used to constrain a solution
504    
505 jfenwick 2625 :return: the mask marking the locations of the constraints
506     :rtype: `escript.Scalar`
507 gross 904 """
508     out_loc=0
509 gross 906 if not self.location_of_constraint0 == None:
510 gross 904 out_loc=wherePositive(out_loc+wherePositive(self.location_of_constraint0))
511 gross 906 if not self.location_of_constraint1 == None:
512 gross 904 out_loc=wherePositive(out_loc+wherePositive(self.location_of_constraint1))
513 gross 906 if not self.location_of_constraint2 == None:
514 gross 904 out_loc=wherePositive(out_loc+wherePositive(self.location_of_constraint2))
515 gross 906 if not self.location_of_constraint3 == None:
516 gross 904 out_loc=wherePositive(out_loc+wherePositive(self.location_of_constraint3))
517     return out_loc
518    
519     def value_of_constraint(self):
520     """
521     return the values used to constrain a solution
522    
523 jfenwick 2625 :return: values to be used at the locations of the constraints. If
524     ``value`` is not given ``None`` is rerturned.
525     :rtype: `escript.Scalar`
526 gross 904 """
527     out_loc=0
528     out=0
529 gross 906 if not self.location_of_constraint0 == None:
530 gross 904 tmp=wherePositive(self.location_of_constraint0)
531 gross 906 out=out*(1.-tmp)+self.value_of_constraint0*tmp
532 gross 904 out_loc=wherePositive(out_loc+tmp)
533 gross 906 if not self.location_of_constraint1 == None:
534 gross 904 tmp=wherePositive(self.location_of_constraint1)
535 gross 906 out=out*(1.-tmp)+self.value_of_constraint1*tmp
536 gross 904 out_loc=wherePositive(out_loc+tmp)
537 gross 906 if not self.location_of_constraint2 == None:
538 gross 904 tmp=wherePositive(self.location_of_constraint2)
539 gross 906 out=out*(1.-tmp)+self.value_of_constraint2*tmp
540 gross 904 out_loc=wherePositive(out_loc+tmp)
541 gross 906 if not self.location_of_constraint3 == None:
542 gross 904 tmp=wherePositive(self.location_of_constraint3)
543 gross 906 out=out*(1.-tmp)+self.value_of_constraint3*tmp
544 gross 904 out_loc=wherePositive(out_loc+tmp)
545     return out
546 jgs 148 # vim: expandtab shiftwidth=4:

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26