12 |
#import urllib |
#import urllib |
13 |
import urllib2 |
import urllib2 |
14 |
|
|
15 |
EPS=1.e-5 |
EPS=1.e-8 |
16 |
|
|
17 |
#==================================================================================================================================================== |
#==================================================================================================================================================== |
18 |
class GridData: |
class GridData: |
354 |
if pl.getDiameter()>self.region.resolution: |
if pl.getDiameter()>self.region.resolution: |
355 |
short_pl=[pl[0]] |
short_pl=[pl[0]] |
356 |
for i in range(1,len(pl)): |
for i in range(1,len(pl)): |
357 |
if pl[i]-short_pl[-1]>self.region.resolution/5: |
if pl[i]-short_pl[-1]>0*EPS+self.region.resolution/10: |
358 |
short_pl.append(pl[i]) |
short_pl.append(pl[i]) |
359 |
elif i==len(pl)-1: |
elif i==len(pl)-1: |
360 |
short_pl[-1]=pl[i] |
short_pl[-1]=pl[i] |
422 |
vertices_on_face.append(q) |
vertices_on_face.append(q) |
423 |
vertices_on_face.sort(self.region.comparePointsOnFace) |
vertices_on_face.sort(self.region.comparePointsOnFace) |
424 |
index=0 |
index=0 |
425 |
walking_on_water=west_south_is_water |
walking_on_water=True |
426 |
l=len(vertices_on_face) |
l=len(vertices_on_face) |
427 |
while index<l: |
while index<l: |
428 |
p1=vertices_on_face[(index+1)%l] |
p1=vertices_on_face[(index+1)%l] |
610 |
finley_file.close() |
finley_file.close() |
611 |
# get the mesh |
# get the mesh |
612 |
out=ReadMesh("%s.msh"%self.fn) |
out=ReadMesh("%s.msh"%self.fn) |
613 |
# os.remove("%s.msh"%self.fn) |
os.remove("%s.msh"%self.fn) |
614 |
return out |
return out |
615 |
|
|
616 |
|
|
692 |
dx=x_grd[1]-ox |
dx=x_grd[1]-ox |
693 |
nx=1 |
nx=1 |
694 |
nx=1 |
nx=1 |
695 |
while abs(x_grd[nx]-ox)>1.e-10*diam: |
while abs(x_grd[nx]-ox)>EPS*diam: |
696 |
nx+=1 |
nx+=1 |
697 |
dy=y_grd[nx]-oy |
dy=y_grd[nx]-oy |
698 |
ny=len(x_grd)/nx |
ny=len(x_grd)/nx |
789 |
x_lat_hat=-b*(x_long-mid_long)+a*(x_lat-mid_lat) |
x_lat_hat=-b*(x_long-mid_long)+a*(x_lat-mid_lat) |
790 |
# x_lat-direction |
# x_lat-direction |
791 |
s=abs(x_lat_hat)-self.width |
s=abs(x_lat_hat)-self.width |
792 |
m=s.whereNegative() |
m=whereNegative(s) |
793 |
f1=(1.-m)*exp(-(s*beta)**2)+m |
f1=(1.-m)*exp(-(s*beta)**2)+m |
794 |
# x_long-direction |
# x_long-direction |
795 |
s=abs(x_long_hat)-dist |
s=abs(x_long_hat)-dist |
796 |
m=s.whereNegative() |
m=whereNegative(s) |
797 |
f0=(1.-m)*exp(-(s*beta)**2)+m |
f0=(1.-m)*exp(-(s*beta)**2)+m |
798 |
self.wave_height=f1*f0*self.amplitude |
self.wave_height=f1*f0*self.amplitude |
799 |
|
|
827 |
self.__pde=LinearPDE(self.domain) |
self.__pde=LinearPDE(self.domain) |
828 |
self.__pde.setSolverMethod(self.__pde.LUMPING) |
self.__pde.setSolverMethod(self.__pde.LUMPING) |
829 |
self.__pde.setValue(D=1.) |
self.__pde.setValue(D=1.) |
830 |
self.__c2=RegionOnEarthSurface.GRAVITY*self.bathymetry/RegionOnEarthSurface.RADIUS**2 |
self.__c2=RegionOnEarthSurface.GRAVITY*self.bathymetry/(RegionOnEarthSurface.RADIUS*2*numarray.pi/360.)**2 |
831 |
c_max=math.sqrt(Lsup(self.__c2)) |
c_max=math.sqrt(Lsup(self.__c2)) |
832 |
self.__dt=self.safety_factor*inf(self.domain.getSize()/(sqrt(self.__c2)+EPS*c_max)) |
self.__dt=self.safety_factor*inf(self.domain.getSize()/(sqrt(self.__c2)+EPS*c_max)) |
833 |
|
print c_max,inf(self.domain.getSize()) |
834 |
if self.initial_time_step==None: self.initial_time_step=self.__dt |
if self.initial_time_step==None: self.initial_time_step=self.__dt |
835 |
self.trace("maximum wave velocity %s m/sec"%c_max) |
self.trace("maximum wave velocity %s m/sec"%c_max) |
836 |
self.trace("Time step size is %s sec"%self.__dt) |
self.trace("Time step size is %s sec"%self.__dt) |
1120 |
@param dt: |
@param dt: |
1121 |
""" |
""" |
1122 |
if self.t>=self.__next_t: |
if self.t>=self.__next_t: |
1123 |
print self.t,"write ",Lsup(self.wave_height) |
print self.t,"write ",Lsup(self.wave_height)," to ",self.__fn |
1124 |
saveVTK(self.__fn,h=self.wave_height) |
saveVTK(self.__fn,h=self.wave_height) |
1125 |
# vtkobj=... |
# vtkobj=... |
1126 |
# save(self.__frame_name) |
# save(self.__frame_name) |
1141 |
waveMapper = vtk.vtkDataSetMapper() |
waveMapper = vtk.vtkDataSetMapper() |
1142 |
waveMapper.SetInput(waveGrid) |
waveMapper.SetInput(waveGrid) |
1143 |
waveMapper.SetLookupTable(self.lutTrans) |
waveMapper.SetLookupTable(self.lutTrans) |
1144 |
waveMapper.SetScalarRange(-self.max_height, self.max_height) |
waveMapper.SetScalarRange(-self.max_height*0.3, self.max_height*0.3) |
1145 |
|
|
1146 |
self.waveActor.SetMapper(waveMapper) |
self.waveActor.SetMapper(waveMapper) |
1147 |
|
|
1219 |
from esys.modellib.input import Sequencer |
from esys.modellib.input import Sequencer |
1220 |
|
|
1221 |
oc=OceanRegionCollector() |
oc=OceanRegionCollector() |
1222 |
oc.resolution=2 |
oc.north= 26.7 |
1223 |
oc.south=-45.5 |
oc.west= 102.9 |
1224 |
oc.north=-5.4 |
oc.range360= True |
1225 |
oc.east=161.1 |
oc.east= 232.6 |
1226 |
oc.west=108.2 |
oc.resolution= 1. |
1227 |
oc.range360=True |
oc.south= -71.3 |
1228 |
|
|
1229 |
|
|
1230 |
b=Bathymetry() |
b=Bathymetry() |
1241 |
|
|
1242 |
src=TsunamiSource() |
src=TsunamiSource() |
1243 |
src.domain=Link(oreg,"domain") |
src.domain=Link(oreg,"domain") |
1244 |
src.start_lat=-10. |
src.decay_zone= 0.01 |
1245 |
src.end_lat=-12. |
src.end_long= 185. |
1246 |
src.start_long=110. |
src.end_lat= -37. |
1247 |
src.end_long=120. |
src.width=0.5 |
1248 |
src.width=1. |
src.start_long= 174. |
1249 |
src.decay_zone=0.01 |
src.start_lat= -15. |
1250 |
src.amplitude=1. |
src.amplitude= 3 |
1251 |
|
|
1252 |
ts=TsunamiInDeepWater() |
ts=TsunamiInDeepWater() |
1253 |
ts.domain=Link(oreg,"domain") |
ts.domain=Link(oreg,"domain") |
1256 |
ts.bathymetry=Link(oreg,"bathymetry") |
ts.bathymetry=Link(oreg,"bathymetry") |
1257 |
|
|
1258 |
sq=Sequencer() |
sq=Sequencer() |
1259 |
sq.t_end=100000. |
sq.t_end=15000. |
1260 |
|
|
1261 |
sm=SurfMovie() |
sm=SurfMovie() |
1262 |
sm.bathymetry=Link(b,"bathymetry") |
sm.bathymetry=Link(b,"bathymetry") |
1263 |
sm.wave_height=Link(ts,"wave_height") |
sm.wave_height=Link(ts,"wave_height") |
1264 |
sm.coastline=Link(oreg,"coastline") |
sm.coastline=Link(oreg,"coastline") |
1265 |
sm.t=Link(sq,"t") |
sm.t=Link(sq,"t") |
|
sm.dt=5000. |
|
1266 |
sm.filename="movie.mpg" |
sm.filename="movie.mpg" |
1267 |
|
sm.north= 8.7 |
1268 |
|
sm.west= 138.9 |
1269 |
|
sm.dt= 50. |
1270 |
|
sm.east= 196.6 |
1271 |
|
sm.south= -53.3 |
1272 |
sm.max_height=Link(src,"amplitude") |
sm.max_height=Link(src,"amplitude") |
1273 |
|
|
1274 |
s=Simulation([sq,oc,b,oreg,src,ts,sm]) |
s=Simulation([sq,oc,b,oreg,src,ts,sm]) |