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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3379 - (show annotations)
Wed Nov 24 04:48:49 2010 UTC (8 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 5713 byte(s)
some clarification on lumping
1 """
2
3 a simple comparison for row-sum and HRZ lumping in case of the wave equation
4
5 """
6
7 ########################################################
8 #
9 # Copyright (c) 2003-2010 by University of Queensland
10 # Earth Systems Science Computational Center (ESSCC)
11 # http://www.uq.edu.au/esscc
12 #
13 # Primary Business: Queensland, Australia
14 # Licensed under the Open Software License version 3.0
15 # http://www.opensource.org/licenses/osl-3.0.php
16 #
17 ########################################################
18
19 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
20 Earth Systems Science Computational Center (ESSCC)
21 http://www.uq.edu.au/esscc
22 Primary Business: Queensland, Australia"""
23 __license__="""Licensed under the Open Software License version 3.0
24 http://www.opensource.org/licenses/osl-3.0.php"""
25 __url__="https://launchpad.net/escript-finley"
26
27 from esys.escript import *
28 from esys.escript.linearPDEs import LinearSinglePDE
29 from esys.escript.pdetools import Locator
30 from esys import finley
31 from math import pi
32
33
34
35 c=1.
36 n=5.
37
38 def ref_u(x,t):
39 return sin(pi*n*(x[0]-c*t))
40
41 def runValetAcceleration(order):
42 domain=finley.Rectangle(100,10,order)
43 x=domain.getX()
44
45
46 # test Velet scheme
47 dt=inf(domain.getSize()/c)*(1./6.)
48 q=whereZero(x[0])+whereZero(x[0]-1.)
49
50 mypde_f=LinearSinglePDE(domain)
51 mypde_f.setSymmetryOn()
52 mypde_f.setValue(D=1,q=q)
53 u_f_old=ref_u(x,-dt)
54 u_f=ref_u(x,0)
55
56 mypde_HRZ=LinearSinglePDE(domain)
57 mypde_HRZ.getSolverOptions().setSolverMethod(mypde_HRZ.getSolverOptions().HRZ_LUMPING)
58 mypde_HRZ.setValue(D=1,q=q)
59 u_HRZ_old=ref_u(x,-dt)
60 u_HRZ=ref_u(x,0)
61
62 mypde_RS=LinearSinglePDE(domain)
63 mypde_RS.getSolverOptions().setSolverMethod(mypde_RS.getSolverOptions().ROWSUM_LUMPING)
64 mypde_RS.setValue(D=1,q=q)
65 u_RS_old=ref_u(x,-dt)
66 u_RS=ref_u(x,0)
67
68 l=Locator(domain,[0.5,0.5])
69 t=0
70
71 u=ref_u(x,t)
72 t_list=[t]
73 u_list=[l(u)]
74 f_list=[l(u_f)]
75 HRZ_list=[l(u_HRZ)]
76 RS_list=[l(u_RS)]
77 print t_list[-1], u_list[-1], f_list[-1], HRZ_list[-1] , RS_list[-1]
78
79 while t< 4/n/c:
80 t+=dt
81 u=ref_u(x,t)
82 mypde_f.setValue(X=-c**2*grad(u_f), r=-c**2*u)
83 mypde_HRZ.setValue(X=-c**2*grad(u_HRZ), r=-c**2*u)
84 mypde_RS.setValue(X=-c**2*grad(u_RS), r=-c**2*u)
85
86 a_f=mypde_f.getSolution()
87 a_HRZ=mypde_HRZ.getSolution()
88 a_RS=mypde_RS.getSolution()
89
90 u_f, u_f_old = 2*u_f-u_f_old + dt**2*a_f , u_f
91 u_HRZ, u_HRZ_old = 2*u_HRZ-u_HRZ_old + dt**2*a_HRZ , u_HRZ
92 u_RS, u_RS_old = 2*u_RS-u_RS_old + dt**2*a_RS , u_RS
93
94 t_list.append(t)
95 u_list.append(l(u))
96 f_list.append(l(u_f))
97 HRZ_list.append(l(u_HRZ))
98 RS_list.append(l(u_RS))
99 print t_list[-1], u_list[-1], f_list[-1], HRZ_list[-1] , RS_list[-1]
100
101
102 import matplotlib.pyplot as plt
103 if getMPIRankWorld() == 0:
104 plt.clf()
105 plt.plot(t_list, u_list, '-', label="exact", linewidth=1)
106 plt.plot(t_list, f_list, '-', label="full", linewidth=1)
107 plt.plot(t_list, HRZ_list, '-', label="HRZ lumping", linewidth=1)
108 plt.plot(t_list, RS_list, '-', label="row sum lumping", linewidth=1)
109 plt.axis([0.,max(t_list),-1.3,1.3])
110 plt.xlabel('time')
111 plt.ylabel('displacement')
112 plt.legend()
113 plt.savefig('lumping_valet_a_%d.png'%order, format='png')
114
115 def runValetDisplacement(order):
116 domain=finley.Rectangle(100,10,order)
117 x=domain.getX()
118
119
120 # test Velet scheme
121 dt=inf(domain.getSize()/c)*(1./6.)
122 q=whereZero(x[0])+whereZero(x[0]-1.)
123
124 mypde_f=LinearSinglePDE(domain)
125 mypde_f.setSymmetryOn()
126 mypde_f.setValue(D=1,q=q)
127 u_f_old=ref_u(x,-dt)
128 u_f=ref_u(x,0)
129
130 mypde_HRZ=LinearSinglePDE(domain)
131 mypde_HRZ.getSolverOptions().setSolverMethod(mypde_HRZ.getSolverOptions().HRZ_LUMPING)
132 mypde_HRZ.setValue(D=1,q=q)
133 u_HRZ_old=ref_u(x,-dt)
134 u_HRZ=ref_u(x,0)
135
136 mypde_RS=LinearSinglePDE(domain)
137 mypde_RS.getSolverOptions().setSolverMethod(mypde_RS.getSolverOptions().ROWSUM_LUMPING)
138 mypde_RS.setValue(D=1,q=q)
139 u_RS_old=ref_u(x,-dt)
140 u_RS=ref_u(x,0)
141
142 l=Locator(domain,[0.5,0.5])
143 t=0
144
145 u=ref_u(x,t)
146 t_list=[t]
147 u_list=[l(u)]
148 f_list=[l(u_f)]
149 HRZ_list=[l(u_HRZ)]
150 RS_list=[l(u_RS)]
151 print t_list[-1], u_list[-1], f_list[-1], HRZ_list[-1] , RS_list[-1]
152
153 while t< 4/n/c:
154 t+=dt
155 u=ref_u(x,t)
156 mypde_f.setValue(X=-(c*dt)**2*grad(u_f), r=u, Y=2*u_f-u_f_old )
157 mypde_HRZ.setValue(X=-(c*dt)**2*grad(u_HRZ), r=u, Y= 2*u_HRZ-u_HRZ_old )
158 mypde_RS.setValue(X=-(c*dt)**2*grad(u_RS), r=u, Y= 2*u_RS-u_RS_old )
159
160 u_f, u_f_old =mypde_f.getSolution(), u_f
161 u_HRZ, u_HRZ_old =mypde_HRZ.getSolution(), u_HRZ
162 u_RS, u_RS_old =mypde_RS.getSolution(), u_RS
163
164 t_list.append(t)
165 u_list.append(l(u))
166 f_list.append(l(u_f))
167 HRZ_list.append(l(u_HRZ))
168 RS_list.append(l(u_RS))
169 print t_list[-1], u_list[-1], f_list[-1], HRZ_list[-1] , RS_list[-1]
170
171
172 import matplotlib.pyplot as plt
173 if getMPIRankWorld() == 0:
174 plt.clf()
175 plt.plot(t_list, u_list, '-', label="exact", linewidth=1)
176 plt.plot(t_list, f_list, '-', label="full", linewidth=1)
177 plt.plot(t_list, HRZ_list, '-', label="HRZ lumping", linewidth=1)
178 plt.plot(t_list, RS_list, '-', label="row sum lumping", linewidth=1)
179 plt.axis([0.,max(t_list),-1.3,1.3])
180 plt.xlabel('time')
181 plt.ylabel('displacement')
182 plt.legend()
183 plt.savefig('lumping_valet_u_%d.png'%order, format='png')
184
185
186 #runValetAcceleration(1)
187 #runValetAcceleration(2)
188 runValetDisplacement(1)
189 runValetDisplacement(2)
190

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26