# Contents of /trunk/finley/test/python/FCT_test2.py

Revision 1400 - (show annotations)
Thu Jan 24 06:04:31 2008 UTC (13 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 5364 byte(s)
```better test example for upwinding added
```
 1 # 2 # upwinding test moving a Gaussian hill around 3 # 4 # we solve U_,t - E *u_,ii + v_i u_,i =0 (E is small) 5 # 6 # the solution is given as u(x,t)=1/(4*pi*E*t)^{dim/2} * exp ( - |x-x_0(t)|^2/(4*E*t) ) 7 # 8 # where x_0(t) = [ cos(OMEGA0*T0)*0.5,-sin(OMEGA0*T0)*0.5 ] and v=[-y,x]*OMEGA0 for dim=2 and 9 # 10 # x_0(t) = [ cos(OMEGA0*T0)*0.5,-sin(OMEGA0*T0)*0.5 ] and v=[-y,x]*OMEGA0 for dim=3 11 # 12 # the solution is started from some time T0>0. 13 # 14 # We are using five quality messurements for u_h 15 # 16 # - inf(u_h) > 0 17 # - sup(u_h)/sup(u(x,t)) = sup(u_h)*(4*pi*E*t)^{dim/2} ~ 1 18 # - integrate(u_h) ~ 1 19 # - | x_0h-x_0 | ~ 0 where x_0h = integrate(x*u_h) 20 # - sigma_h/4*E*t ~ 1 where sigma_h=sqrt(integrate(length(x-x0h)**2 * u_h) * (DIM==3 ? sqrt(2./3.) :1 ) 21 # 22 # 23 from math import pi, ceil 24 NE=128/2 25 DIM=2 26 THETA=0.5 27 OMEGA0=1. 28 ALPHA=pi/4 29 T0=0.5*pi 30 E=5e-4 * 2 31 TEST_SUPG=False 32 33 def getCenter(t): 34 if DIM==2: 35 x0=[cos(OMEGA0*t)*0.5,-sin(OMEGA0*t)*0.5] 36 else: 37 x0=[cos(ALPHA)*cos(OMEGA0*t)*0.5,-sin(OMEGA0*t)*0.5,-sin(ALPHA)*cos(OMEGA0*t)*0.5] 38 return x0 39 def QUALITY(t,u_h): 40 dom=u_h.getDomain() 41 x=dom.getX() 42 a=inf(u_h) 43 b=sup(u_h)*(4*pi*E*t)**(DIM/2.)-1. 44 c=integrate(u_h,Function(dom))-1. 45 x0=getCenter(t) 46 x0h=integrate(x*u_h,Function(dom)) 47 d=length(x0-x0h) 48 sigma_h2=sqrt(integrate(length(x-x0h)**2 * u_h, Function(dom))) 49 if DIM == 3: sigma_h2*=sqrt(2./3.) 50 e=sigma_h2/sqrt(4*E*t)-1 51 return a,b,c,d,e 52 53 54 from esys.escript import * 55 from esys.escript.linearPDEs import LinearSinglePDE 56 57 class TransportPDE(object): 58 def __init__(self,domain,num_equations=1,theta=0.,dt_max=-1.,trace=True): 59 self.__domain=domain 60 self.__num_equations=num_equations 61 self.__theta=theta 62 self.__dt_max=dt_max 63 self.__transport_problem=None 64 self.__trace=trace 65 self.__matrix_type=0 66 def trace(self,text): 67 if self.__trace: print text 68 69 def getDomain(self): 70 return self.__domain 71 def getTheta(self): 72 return self.__theta 73 def getDt_max(self): 74 return self.__dt_max 75 def getNumEquations(self): 76 return self.__num_equations 77 def reduced(self): 78 return False 79 def getFunctionSpace(self): 80 if self.reduced(): 81 return ReducedSolution(self.getDomain()) 82 else: 83 return Solution(self.getDomain()) 84 85 86 def __getNewTransportProblem(self): 87 """ 88 returns an instance of a new operator 89 """ 90 self.trace("New Transport problem is allocated.") 91 return self.getDomain().newTransportProblem( \ 92 self.getTheta(), 93 self.getDt_max(), 94 self.getNumEquations(), \ 95 self.getFunctionSpace(), \ 96 self.__matrix_type) 97 98 def setValue(self,M=Data(),A=Data(),B=Data(),C=Data(),D=Data(),X=Data(),Y=Data(), 99 d=Data(),y=Data(),d_contact=Data(),y_contact=Data()): 100 self.__transport_problem=self.__getNewTransportProblem() 101 if self.getNumEquations() ==1 : 102 self.__source=Data(0.0,(),self.getFunctionSpace()) 103 else: 104 self.__source=Data(0.0,(self.getNumEquations(),),self.getFunctionSpace()) 105 self.getDomain().addPDEToTransportProblem( 106 self.__transport_problem, 107 self.__source, 108 M,A,B,C,D,X,Y,d,y,d_contact,y_contact) 109 110 def setInitialSolution(self,u): 111 self.__transport_problem.setInitialValue(interpolate(u,self.getFunctionSpace())) 112 def solve(self,dt): 113 return self.__transport_problem.solve(self.__source,dt,{"verbose" : True , "tolerance" : 1.e-8}) 114 115 116 from esys.finley import Rectangle, Brick 117 118 119 if DIM==2: 120 dom=Rectangle(NE,NE) 121 else: 122 dom=Brick(NE,NE,NE) 123 dom.setX(2*dom.getX()-1) 124 125 # set initial value 126 x=dom.getX() 127 u0=1/(4.*pi*E*T0)**(DIM/2.)*exp(-length(dom.getX()-getCenter(T0))**2/(4.*E*T0)) 128 129 print "QUALITY ",QUALITY(T0,u0) 130 131 x=Function(dom).getX() 132 if DIM == 2: 133 V=OMEGA0*(x[0]*[0,1]+x[1]*[-1,0]) 134 else: 135 V=OMEGA0*(x[0]*[0,cos(ALPHA),0]+x[1]*[-cos(ALPHA),0,sin(ALPHA)]+x[2]*[0.,-sin(ALPHA),0.]) 136 #=================== 137 fc=TransportPDE(dom,num_equations=1,theta=THETA) 138 x=Function(dom).getX() 139 fc.setValue(M=Scalar(1.,Function(dom)),C=V,A=-Scalar(E,Function(dom))*kronecker(dom)) 140 #============== 141 if TEST_SUPG: 142 supg=LinearSinglePDE(dom) 143 supg.setValue(D=1.) 144 supg.setSolverMethod(supg.LUMPING) 145 dt_supg=1./(1./inf(dom.getSize()/length(V))+1./inf(dom.getSize()**2/E))*0.3 146 u_supg=u0*1. 147 148 c=0 149 saveVTK("u.%s.xml"%c,u=u0) 150 fc.setInitialSolution(u0) 151 dt=1.266539e-02*10 152 t=T0 153 while t