| 1 |
######################################################## |
| 2 |
# |
| 3 |
# Copyright (c) 2003-2010 by University of Queensland |
| 4 |
# Earth Systems Science Computational Center (ESSCC) |
| 5 |
# http://www.uq.edu.au/esscc |
| 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 |
######################################################## |
| 12 |
|
| 13 |
__copyright__="""Copyright (c) 2003-2010 by University of Queensland |
| 14 |
Earth Systems Science Computational Center (ESSCC) |
| 15 |
http://www.uq.edu.au/esscc |
| 16 |
Primary Business: Queensland, Australia""" |
| 17 |
__license__="""Licensed under the Open Software License version 3.0 |
| 18 |
http://www.opensource.org/licenses/osl-3.0.php""" |
| 19 |
__url__="https://launchpad.net/escript-finley" |
| 20 |
|
| 21 |
# import tools |
| 22 |
from esys.escript import * |
| 23 |
from esys.escript.linearPDEs import LinearPDE |
| 24 |
from esys.dudley import Rectangle |
| 25 |
from esys.weipa import saveVTK |
| 26 |
# end of simulation time |
| 27 |
t_end=0.1 |
| 28 |
# dimensions: |
| 29 |
L0=1.;L1=1. |
| 30 |
# location, size and value of heat source |
| 31 |
xc=[0.3,0.4]; r=0.1; Qc=3000 |
| 32 |
# material parameter |
| 33 |
k=1.; rhocp=100; |
| 34 |
# bottom temperature: |
| 35 |
T_bot=100 |
| 36 |
# generate domain: |
| 37 |
mydomain=Rectangle(l0=L0,l1=L1,n0=20,n1=20) |
| 38 |
x=mydomain.getX() |
| 39 |
# set boundray temperature: |
| 40 |
T_D=T_bot/L1*(L1-x[1]) |
| 41 |
# set heat source: |
| 42 |
Q=Qc*whereNegative(length(x-xc)-r) |
| 43 |
# time step size: |
| 44 |
dt=0.01 |
| 45 |
# or use adaptive choice dt=0.05*inf(rhocp*mydomain.getSize()**2/k) |
| 46 |
print "time step size = ",dt |
| 47 |
# generate domain: |
| 48 |
mypde=LinearPDE(mydomain) |
| 49 |
mypde.setSymmetryOn() |
| 50 |
# set PDE coefficients: |
| 51 |
mypde.setValue(D=rhocp, |
| 52 |
r=T_D, q=whereZero(x[1])+whereZero(x[1]-L1)) |
| 53 |
# initial temperature |
| 54 |
T=T_D |
| 55 |
# step counter and time marker: |
| 56 |
N=0; t=0 |
| 57 |
# stop when t_end is reached: |
| 58 |
while t<t_end: |
| 59 |
print N,"-th time step t=",t," T_max=", Lsup(T) |
| 60 |
# update PDE coefficient: |
| 61 |
mypde.setValue(Y=rhocp*T+dt*Q, X=-k*dt*grad(T)) |
| 62 |
# new temperature: |
| 63 |
T=mypde.getSolution() |
| 64 |
# save as VTK for visualisation: |
| 65 |
# saveVTK("u.%s.vtu"%N,T=T) |
| 66 |
# increase counter and marker: |
| 67 |
N+=1; t+=dt |
Parent Directory
| 
Revision Log