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

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26