/[escript]/trunk/modellib/py_src/input.py
ViewVC logotype

Contents of /trunk/modellib/py_src/input.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4657 - (show annotations)
Thu Feb 6 06:12:20 2014 UTC (5 years, 8 months ago) by jfenwick
File MIME type: text/x-python
File size: 20698 byte(s)
I changed some files.
Updated copyright notices, added GeoComp.



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