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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5288 - (show annotations)
Tue Dec 2 23:18:40 2014 UTC (4 years, 9 months ago) by sshaw
File MIME type: text/x-python
File size: 2540 byte(s)
fixing tests for cases where required domains not built
1 from __future__ import division
2 from __future__ import print_function
3 ##############################################################################
4 #
5 # Copyright (c) 2003-2014 by University of Queensland
6 # http://www.uq.edu.au
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
13 # Development 2012-2013 by School of Earth Sciences
14 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
15 #
16 ##############################################################################
17
18 __copyright__="""Copyright (c) 2003-2014 by University of Queensland
19 http://www.uq.edu.au
20 Primary Business: Queensland, Australia"""
21 __license__="""Licensed under the Open Software License version 3.0
22 http://www.opensource.org/licenses/osl-3.0.php"""
23 __url__="https://launchpad.net/escript-finley"
24
25 # import tools
26 from esys.escript import *
27 from esys.escript.linearPDEs import LinearPDE
28 try:
29 from esys.dudley import Rectangle
30 HAVE_DUDLEY = True
31 except ImportError:
32 HAVE_DUDLEY = False
33 from esys.weipa import saveVTK
34
35 if not HAVE_DUDLEY:
36 print("Dudley module not available")
37 else:
38 # end of simulation time
39 t_end=0.1
40 # dimensions:
41 L0=1.;L1=1.
42 # location, size and value of heat source
43 xc=[0.3,0.4]; r=0.1; Qc=3000
44 # material parameter
45 k=1.; rhocp=100;
46 # bottom temperature:
47 T_bot=100
48 # generate domain:
49 mydomain=Rectangle(l0=L0,l1=L1,n0=20,n1=20)
50 x=mydomain.getX()
51 # set boundray temperature:
52 T_D=T_bot/L1*(L1-x[1])
53 # set heat source:
54 Q=Qc*whereNegative(length(x-xc)-r)
55 # time step size:
56 dt=0.01
57 # or use adaptive choice dt=0.05*inf(rhocp*mydomain.getSize()**2/k)
58 print("time step size = %s"%dt)
59 # generate domain:
60 mypde=LinearPDE(mydomain)
61 mypde.setSymmetryOn()
62 # set PDE coefficients:
63 mypde.setValue(D=rhocp,
64 r=T_D, q=whereZero(x[1])+whereZero(x[1]-L1))
65 # initial temperature
66 T=T_D
67 # step counter and time marker:
68 N=0; t=0
69 # stop when t_end is reached:
70 while t<t_end:
71 print("time step %d, t=%s, T_max=%s"%(N, t, Lsup(T)))
72 # update PDE coefficient:
73 mypde.setValue(Y=rhocp*T+dt*Q, X=-k*dt*grad(T))
74 # new temperature:
75 T=mypde.getSolution()
76 # save as VTK for visualisation:
77 # saveVTK("u.%s.vtu"%N,T=T)
78 # increase counter and marker:
79 N+=1; t+=dt
80

  ViewVC Help
Powered by ViewVC 1.1.26