# Annotation of /trunk/finley/test/python/FCT_test1.py

Revision 3981 - (hide annotations)
Fri Sep 21 02:47:54 2012 UTC (7 years, 4 months ago) by jfenwick
File MIME type: text/x-python
File size: 3562 byte(s)
```First pass of updating copyright notices
```
 1 ksteube 1811 2 jfenwick 3981 ############################################################################## 3 gross 1411 # 4 jfenwick 3911 # Copyright (c) 2003-2012 by University of Queensland 5 jfenwick 3981 6 ksteube 1811 # 7 # Primary Business: Queensland, Australia 8 # Licensed under the Open Software License version 3.0 9 10 # 11 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC) 12 # Development since 2012 by School of Earth Sciences 13 # 14 ############################################################################## 15 ksteube 1811 16 jfenwick 3911 __copyright__="""Copyright (c) 2003-2012 by University of Queensland 17 jfenwick 3981 http://www.uq.edu.au 18 ksteube 1811 Primary Business: Queensland, Australia""" 19 __license__="""Licensed under the Open Software License version 3.0 20 21 jfenwick 2344 __url__= 22 ksteube 1811 23 # 24 gross 1411 # 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 gross 3793 from esys.escript.linearPDEs import TransportPDE 47 gross 3808 from esys.finley import Rectangle, Brick 48 #from esys.ripley import Rectangle, Brick 49 caltinay 3346 from esys.weipa import saveVTK 50 gross 1411 from math import pi, ceil 51 NE=128 52 gross 3793 #NE=4 53 gross 1411 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 gross 3808 dom=Rectangle(NE,NE) 64 gross 1562 u0=dom.getX()[0] 65 caltinay 2534 # saveVTK("u.%s.vtu"%0,u=u0) 66 gross 1562 # print "XX"*80 67 gross 3808 68 # set initial value 69 gross 3793 #dom.setX(2*dom.getX()-1) 70 gross 3808 #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 gross 1411 74 gross 3808 #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 gross 1411 x=dom.getX() 81 82 gross 3808 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 gross 1411 x=Function(dom).getX() 100 gross 3808 V=OMEGA0*((0.5-x[0])*[0,1]+(0.5-x[1])*[-1,0]) 101 gross 1411 #=================== 102 gross 3822 103 gross 3793 fc=TransportPDE(dom,numEquations=1) 104 fc.getSolverOptions().setVerbosityOn() 105 gross 3822 #fc.getSolverOptions().setODESolver(fc.getSolverOptions().BACKWARD_EULER) 106 gross 3793 fc.getSolverOptions().setODESolver(fc.getSolverOptions().LINEAR_CRANK_NICOLSON) 107 fc.getSolverOptions().setODESolver(fc.getSolverOptions().CRANK_NICOLSON) 108 gross 1411 x=Function(dom).getX() 109 gross 3793 fc.setValue(M=1,C=V) 110 gross 1411 111 c=0 112 gross 3793 saveVTK("u.%s.vtu"%c,u=u0) 113 gross 1411 fc.setInitialSolution(u0) 114 gross 3808 dt=fc.getSafeTimeStepSize() 115 #dt=1.e-3 116 gross 3793 print "dt = ",dt 117 gross 1411 t=T0 118 jfenwick 3772 print("QUALITY FCT: time = %s pi"%(t/pi),inf(u0),sup(u0),integrate(u0)) 119 gross 3808 #T_END=200*dt 120 gross 1411 while t