/[escript]/trunk/doc/examples/geotutorial/steadystate_variablek.py
ViewVC logotype

Contents of /trunk/doc/examples/geotutorial/steadystate_variablek.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (2 years, 3 months ago) by jfenwick
File MIME type: text/x-python
File size: 2054 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1 from __future__ import division
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2018 by The University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Apache License, version 2.0
9 # http://www.apache.org/licenses/LICENSE-2.0
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-2018 by The University of Queensland
18 http://www.uq.edu.au
19 Primary Business: Queensland, Australia"""
20 __license__="""Licensed under the Apache License, version 2.0
21 http://www.apache.org/licenses/LICENSE-2.0"""
22 __url__="https://launchpad.net/escript-finley"
23
24 # import tools
25 from esys.escript import *
26 from esys.escript.linearPDEs import LinearPDE
27 try:
28 from esys.dudley import Rectangle
29 HAVE_DUDLEY = True
30 except ImportError:
31 HAVE_DUDLEY = False
32 from esys.weipa import saveVTK
33
34 if not HAVE_DUDLEY:
35 print("Dudley module not available")
36 else:
37 # dimensions:
38 L0=1.;L1=1.;
39 # height of k interface:
40 H=L1*0.75
41 # bottom temperature:
42 T_bot=100
43 # location, size and value of heat source
44 xc=[0.3,0.4]; r=0.1; Qc=3000
45 # two values for k
46 k0=1; k1=10
47 # create domain:
48 mydomain=Rectangle(l0=L0,l1=L1,n0=20,n1=20)
49 x=mydomain.getX()
50 # set variable k
51 k=k0+(k1-k0)*wherePositive(x[1]-H)
52 # boundary temperature
53 T_D=T_bot/L1*(L1-x[1])
54 # heat source
55 Q=Qc*whereNegative(length(x-xc)-r)
56 # create PDE and set coefficients:
57 mypde=LinearPDE(mydomain)
58 mypde.setSymmetryOn()
59 # set PDE coefficients:
60 mypde.setValue(A=k*kronecker(mydomain),Y=Q, r=T_D, \
61 q=whereZero(x[1])+whereZero(x[1]-L1))
62 # get temperature:
63 T=mypde.getSolution()
64 # save as VTK for visualisation:
65 saveVTK("u.vtu",T=T)
66

  ViewVC Help
Powered by ViewVC 1.1.26