/[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 6523 - (hide annotations)
Tue Mar 7 06:50:34 2017 UTC (2 years, 9 months ago) by jfenwick
File MIME type: text/x-python
File size: 5840 byte(s)
UPDATE THE COPYRIGHT DATES

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

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26