/[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 4657 - (hide annotations)
Thu Feb 6 06:12:20 2014 UTC (5 years, 10 months ago) by jfenwick
File MIME type: text/x-python
File size: 5847 byte(s)
I changed some files.
Updated copyright notices, added GeoComp.



1 gross 3379 """
2    
3     a simple comparison for row-sum and HRZ lumping in case of the wave equation
4    
5     """
6    
7 jfenwick 3981 ##############################################################################
8 gross 3379 #
9 jfenwick 4657 # Copyright (c) 2003-2014 by University of Queensland
10 jfenwick 3981 # http://www.uq.edu.au
11 gross 3379 #
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 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
17 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
18     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
19 jfenwick 3981 #
20     ##############################################################################
21 gross 3379
22 jfenwick 4657 __copyright__="""Copyright (c) 2003-2014 by University of Queensland
23 jfenwick 3981 http://www.uq.edu.au
24 gross 3379 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
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(mypde_HRZ.getSolverOptions().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(mypde_RS.getSolverOptions().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 jfenwick 3772 print(t_list[-1], u_list[-1], f_list[-1], HRZ_list[-1] , RS_list[-1])
80 gross 3379
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 jfenwick 3772 print(t_list[-1], u_list[-1], f_list[-1], HRZ_list[-1] , RS_list[-1])
102 gross 3379
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(mypde_HRZ.getSolverOptions().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(mypde_RS.getSolverOptions().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 jfenwick 3772 print(t_list[-1], u_list[-1], f_list[-1], HRZ_list[-1] , RS_list[-1])
154 gross 3379
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 jfenwick 3772 print(t_list[-1], u_list[-1], f_list[-1], HRZ_list[-1] , RS_list[-1])
172 gross 3379
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