/[escript]/trunk/finley/test/python/FCT_test1.py
ViewVC logotype

Contents of /trunk/finley/test/python/FCT_test1.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3981 - (show annotations)
Fri Sep 21 02:47:54 2012 UTC (7 years, 2 months ago) by jfenwick
File MIME type: text/x-python
File size: 3562 byte(s)
First pass of updating copyright notices
1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2012 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 since 2012 by School of Earth Sciences
13 #
14 ##############################################################################
15
16 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
17 http://www.uq.edu.au
18 Primary Business: Queensland, Australia"""
19 __license__="""Licensed under the Open Software License version 3.0
20 http://www.opensource.org/licenses/osl-3.0.php"""
21 __url__="https://launchpad.net/escript-finley"
22
23 #
24 # upwinding test moving a Gaussian hill around
25 #
26 # we solve U_,t + v_i u_,i =0
27 #
28 # the solution is given as u(x,t)=1/(4*pi*E*t)^{dim/2} * exp ( - |x-x_0(t)|^2/(4*E*t) )
29 #
30 # where x_0(t) = [ cos(OMEGA0*T0)*0.5,-sin(OMEGA0*T0)*0.5 ] and v=[-y,x]*OMEGA0 for dim=2 and
31 #
32 # x_0(t) = [ cos(OMEGA0*T0)*0.5,-sin(OMEGA0*T0)*0.5 ] and v=[-y,x]*OMEGA0 for dim=3
33 #
34 # the solution is started from some time T0>0.
35 #
36 # We are using five quality messurements for u_h
37 #
38 # - inf(u_h) > 0
39 # - sup(u_h)/sup(u(x,t)) = sup(u_h)*(4*pi*E*t)^{dim/2} ~ 1
40 # - integrate(u_h) ~ 1
41 # - | x_0h-x_0 | ~ 0 where x_0h = integrate(x*u_h)
42 # - sigma_h/4*E*t ~ 1 where sigma_h=sqrt(integrate(length(x-x0h)**2 * u_h) * (DIM==3 ? sqrt(2./3.) :1 )
43 #
44 #
45 from esys.escript import *
46 from esys.escript.linearPDEs import TransportPDE
47 from esys.finley import Rectangle, Brick
48 #from esys.ripley import Rectangle, Brick
49 from esys.weipa import saveVTK
50 from math import pi, ceil
51 NE=128
52 #NE=4
53 DIM=2
54 THETA=0.5
55 OMEGA0=1.
56 ALPHA=pi/4
57 T0=0
58 T_END=2.*pi
59 dt=1e-3*10*10
60 E=1.e-3
61
62
63 dom=Rectangle(NE,NE)
64 u0=dom.getX()[0]
65 # saveVTK("u.%s.vtu"%0,u=u0)
66 # print "XX"*80
67
68 # set initial value
69 #dom.setX(2*dom.getX()-1)
70 #x=dom.getX()
71 #r=sqrt(x[0]**2+(x[1]-1./3.)**2)
72 #u0=whereNegative(r-1./3.)*wherePositive(wherePositive(abs(x[0])-0.05)+wherePositive(x[1]-0.5))
73
74 #x=Function(dom).getX()
75 #if DIM == 2:
76 # V=OMEGA0*(x[0]*[0,-1]+x[1]*[1,0])
77 #else:
78 # V=OMEGA0*(x[0]*[0,cos(ALPHA),0]+x[1]*[-cos(ALPHA),0,sin(ALPHA)]+x[2]*[0.,-sin(ALPHA),0.])
79
80 x=dom.getX()
81
82 R0=0.15
83 #cylinder:
84 X0=0.5
85 Y0=0.75
86 r=sqrt((x[0]-X0)**2+(x[1]-Y0)**2)/R0
87 u0=whereNegative(r-1)*wherePositive(wherePositive(abs(x[0]-X0)-0.025)+wherePositive(x[1]-0.85))
88 # cone:
89 X0=0.5
90 Y0=0.25
91 r=sqrt((x[0]-X0)**2+(x[1]-Y0)**2)/R0
92 u0=u0+wherePositive(1-r)*(1-r)
93 #hump
94 X0=0.25
95 Y0=0.5
96 r=sqrt((x[0]-X0)**2+(x[1]-Y0)**2)/R0
97 u0=u0+1./4.*(1+cos(pi*clip(r,maxval=1)))
98
99 x=Function(dom).getX()
100 V=OMEGA0*((0.5-x[0])*[0,1]+(0.5-x[1])*[-1,0])
101 #===================
102
103 fc=TransportPDE(dom,numEquations=1)
104 fc.getSolverOptions().setVerbosityOn()
105 #fc.getSolverOptions().setODESolver(fc.getSolverOptions().BACKWARD_EULER)
106 fc.getSolverOptions().setODESolver(fc.getSolverOptions().LINEAR_CRANK_NICOLSON)
107 fc.getSolverOptions().setODESolver(fc.getSolverOptions().CRANK_NICOLSON)
108 x=Function(dom).getX()
109 fc.setValue(M=1,C=V)
110
111 c=0
112 saveVTK("u.%s.vtu"%c,u=u0)
113 fc.setInitialSolution(u0)
114 dt=fc.getSafeTimeStepSize()
115 #dt=1.e-3
116 print "dt = ",dt
117 t=T0
118 print("QUALITY FCT: time = %s pi"%(t/pi),inf(u0),sup(u0),integrate(u0))
119 #T_END=200*dt
120 while t<T_END:
121
122 print("time step t=",t+dt)
123 u=fc.getSolution(dt)
124 print("QUALITY FCT: time = %s pi"%(t+dt/pi),inf(u),sup(u),integrate(u))
125 saveVTK("u.%s.vtu"%(c+1,),u=u)
126 c+=1
127 t+=dt

  ViewVC Help
Powered by ViewVC 1.1.26