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

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26