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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3911 - (hide annotations)
Thu Jun 14 01:01:03 2012 UTC (7 years, 5 months ago) by jfenwick
File MIME type: text/x-python
File size: 5717 byte(s)
Copyright changes
1 gross 3379 """
2    
3     a simple comparison for row-sum and HRZ lumping in case of the wave equation
4    
5     """
6    
7     ########################################################
8     #
9 jfenwick 3911 # Copyright (c) 2003-2012 by University of Queensland
10 gross 3379 # 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 jfenwick 3911 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
20 gross 3379 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 jfenwick 3772 print(t_list[-1], u_list[-1], f_list[-1], HRZ_list[-1] , RS_list[-1])
78 gross 3379
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 jfenwick 3772 print(t_list[-1], u_list[-1], f_list[-1], HRZ_list[-1] , RS_list[-1])
100 gross 3379
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 jfenwick 3772 print(t_list[-1], u_list[-1], f_list[-1], HRZ_list[-1] , RS_list[-1])
152 gross 3379
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 jfenwick 3772 print(t_list[-1], u_list[-1], f_list[-1], HRZ_list[-1] , RS_list[-1])
170 gross 3379
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