examples/inversion/synthetic_VTI.py

##############################################################################
#
# Copyright (c) 2003-2018 by The University of Queensland
# http://www.uq.edu.au
#
# Primary Business: Queensland, Australia
# Licensed under the Apache License, version 2.0
# http://www.apache.org/licenses/LICENSE-2.0
#
# Development 2012-2013 by School of Earth Sciences
# Development from 2014 by Centre for Geoscience Computing (GeoComp)
#
##############################################################################
from __future__ import print_function, division

__copyright__="""Copyright (c) 2003-2018 by The University of Queensland
http://www.uq.edu.au
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Apache License, version 2.0
http://www.apache.org/licenses/LICENSE-2.0"""
__url__="https://launchpad.net/escript-finley"

from esys.escript import *
from esys.escript import unitsSI as U
from esys.escript.pdetools import Locator
from esys.weipa import saveSilo
from esys.downunder import Ricker, VTIWave, SimpleSEGYWriter
from math import ceil

try:
    from esys.speckley import Brick, Rectangle
    HAVE_SPECKLEY=True
except ImportError:
    HAVE_SPECKLEY=False

if HAVE_SPECKLEY:
    DIM=2          # spatial dimension

    depth=1*U.km    # depth 
    v_p_top=1.5*U.km/U.sec
    v_p_bottom=3*U.km/U.sec
    absorption_zone=100*U.m
    ne_z=50.

    reflector_at=0.5*depth


    t_end=0.008*U.sec #only this low for testing purposes
    frq=8.*U.Hz
    sampling_interval=4*U.msec
    numRcvPerLine=101
    rangeRcv=800*U.m

    # location of source in crossing array lines with in 0..numRcvInLine one needs to be None
    srcEW=numRcvPerLine//2
    srcNS=None

    # dommain dimension
    width_x=rangeRcv + 4*absorption_zone
    width_y=width_x
    #
    # create array 
    #
    receiver_line=[2*absorption_zone + i * (rangeRcv//(numRcvPerLine-1)) for i in range(numRcvPerLine) ]
    #
    #   set source location with tag "source""
    #
    src_tags=["source"]
    if DIM == 2:
       src_locations = [ (receiver_line[srcEW], depth)]
       src_loc_2D=(receiver_line[srcEW], 0.)
    else:
       if srcEW:
          srcNS=numRcvPerLine//2
       elif srcNS:
          srcEW=numRcvPerLine//2
       else:
           raise ValueError("on of the variables srcEW or srcNS must be None!")
       src_locations  = [ (receiver_line[srcEW], receiver_line[srcNS], depth)]
       src_loc_2D=(receiver_line[srcEW], receiver_line[srcNS])
    #
    #   create sensor arrays:
    #
    # East-west line of receiver
    rcvEW_locations=[]
    rgEW=[]
    mid_point=receiver_line[len(receiver_line)//2]

    for ix in range(len(receiver_line)):
            if DIM == 2:
                rcvEW_locations.append((receiver_line[ix], depth))
                rgEW.append( ( receiver_line[ix], 0.) ) 
            else:
               rcvEW_locations.append((receiver_line[ix], mid_point, depth))
               rgEW.append( ( receiver_line[ix], mid_point) ) 
    # North-south line of receiver
    if DIM == 3:
       rcvNS_locations=[]
       rgNS=[]

       for iy in range(len(receiver_line)):
           rcvNS_locations.append((mid_point, receiver_line[iy],  depth))
           rgNS.append( (  mid_point, receiver_line[iy]) ) 
    #
    # create domain:
    #
    order = 5
    if DIM == 2:
        domain=Rectangle(order, ceil(ne_z*width_x/depth),ne_z,l0=width_x,l1=depth, 
                    diracPoints=src_locations, diracTags=src_tags)
    else:
        domain=Brick(order, ceil(ne_z*width_x/depth), ceil(ne_z*width_y/depth),
                    ne_z, l0=width_x, l1=width_y, l2=depth,
                    diracPoints=src_locations, diracTags=src_tags)
    wl=Ricker(frq)

    #======================================================================
    # m=whereNegative(Function(domain).getX()[DIM-1]-reflector_at)
    # v_p=v_p_bottom*m+v_p_top*(1-m)
    v_p=2*U.km/U.sec
    v_s=0.9*U.km/U.sec
    vareps=0.1*0
    gamma=0.15*0
    delta=0.05*0
    rho=2000*U.kg/U.m**3
    src_dir=[0,0,1]

    sw=VTIWave(domain, v_p, v_s, wl, src_tags[0], source_vector = src_dir,
                    eps=vareps, gamma=gamma, delta=delta, rho=rho,
                    absorption_zone=300*U.m, absorption_cut=1e-2, lumping=True)

    locEW=Locator(domain,rcvEW_locations)
    tracerEW_x=SimpleSEGYWriter(receiver_group=rgEW, source=src_loc_2D, sampling_interval=sampling_interval, text='x-displacement - east-west line')
    tracerEW_z=SimpleSEGYWriter(receiver_group=rgEW, source=src_loc_2D, sampling_interval=sampling_interval, text='z-displacement - east-west line')
    if DIM==3:
        locNS=Locator(domain,rcvNS_locations)
        tracerEW_y=SimpleSEGYWriter(receiver_group=rgEW, source=src_loc_2D,
            sampling_interval=sampling_interval,
            text='x-displacement - east-west line')
        tracerNS_x=SimpleSEGYWriter(receiver_group=rgNS, source=src_loc_2D,
            sampling_interval=sampling_interval,
            text='x-displacement - north-south line')
        tracerNS_y=SimpleSEGYWriter(receiver_group=rgNS, source=src_loc_2D,
            sampling_interval=sampling_interval,
            text='y-displacement - north-south line')
        tracerNS_z=SimpleSEGYWriter(receiver_group=rgNS, source=src_loc_2D,
            sampling_interval=sampling_interval,
            text='z-displacement - north-south line')
    if not tracerEW_x.obspy_available():
        print("\nWARNING: obspy not available, SEGY files will not be written\n")
    elif getMPISizeWorld() > 1:
        print("\nWARNING: SEGY files cannot be written with multiple processes\n")

    t=0.
    mkDir('tmp')
    n=0
    while t < t_end:
            t,u = sw.update(t+sampling_interval)
            tracerEW_x.addRecord(locEW(u[0]))
            tracerEW_z.addRecord(locEW(u[DIM-1]))
            if DIM==3:
                   tracerEW_y.addRecord(locEW(u[1]))
                   tracerNS_x.addRecord(locNS(u[0]))
                   tracerNS_y.addRecord(locNS(u[1]))
                   tracerNS_z.addRecord(locNS(u[2]))
            print(t, locEW(u[DIM-1])[len(rgEW)//2-4:len(rgEW)//2+1], wl.getValue(t))
            #if n%5 == 0 : saveSilo("tmp/u_%d.silo"%(n/5,), u=u)
            saveSilo("tmp/u_%d.silo"%(n,), u=u, cycle=n, time=t)
            n+=1
    if tracerEW_x.obspy_available() and getMPISizeWorld() == 1:
        tracerEW_x.write('lineEW_x.sgy')
        tracerEW_z.write('lineEW_z.sgy')
        if DIM == 3: 
            tracerEW_y.write('lineEW_y.sgy')
            tracerNS_x.write('lineNS_x.sgy')
            tracerNS_y.write('lineNS_y.sgy')
            tracerNS_z.write('lineNS_z.sgy')

else: # no speckley
    print("The Speckley module is not available")

1	gross	4548	##############################################################################
2			#
3	jfenwick	6651	# Copyright (c) 2003-2018 by The University of Queensland
4	gross	4548	# http://www.uq.edu.au
5			#
6			# Primary Business: Queensland, Australia
7	jfenwick	6112	# Licensed under the Apache License, version 2.0
8			# http://www.apache.org/licenses/LICENSE-2.0
9	gross	4548	#
10	jfenwick	4657	# Development 2012-2013 by School of Earth Sciences
11			# Development from 2014 by Centre for Geoscience Computing (GeoComp)
12	gross	4548	#
13			##############################################################################
14	sshaw	5666	from __future__ import print_function, division
15	gross	4548
16	jfenwick	6651	__copyright__="""Copyright (c) 2003-2018 by The University of Queensland
17	gross	4548	http://www.uq.edu.au
18			Primary Business: Queensland, Australia"""
19	jfenwick	6112	__license__="""Licensed under the Apache License, version 2.0
20			http://www.apache.org/licenses/LICENSE-2.0"""
21	gross	4548	__url__="https://launchpad.net/escript-finley"
22
23			from esys.escript import *
24			from esys.escript import unitsSI as U
25			from esys.escript.pdetools import Locator
26			from esys.weipa import saveSilo
27	gross	4562	from esys.downunder import Ricker, VTIWave, SimpleSEGYWriter
28	gross	4548	from math import ceil
29
30	caltinay	6119	try:
31			from esys.speckley import Brick, Rectangle
32			HAVE_SPECKLEY=True
33			except ImportError:
34			HAVE_SPECKLEY=False
35	gross	4548
36	caltinay	6119	if HAVE_SPECKLEY:
37			DIM=2 # spatial dimension
38	gross	4548
39	caltinay	6119	depth=1*U.km # depth
40			v_p_top=1.5*U.km/U.sec
41			v_p_bottom=3*U.km/U.sec
42			absorption_zone=100*U.m
43			ne_z=50.
44	gross	4553
45	caltinay	6119	reflector_at=0.5*depth
46	gross	4548
47
48	caltinay	6119	t_end=0.008*U.sec #only this low for testing purposes
49			frq=8.*U.Hz
50			sampling_interval=4*U.msec
51			numRcvPerLine=101
52			rangeRcv=800*U.m
53	gross	4548
54	caltinay	6119	# location of source in crossing array lines with in 0..numRcvInLine one needs to be None
55			srcEW=numRcvPerLine//2
56			srcNS=None
57	gross	4548
58	caltinay	6119	# dommain dimension
59			width_x=rangeRcv + 4*absorption_zone
60			width_y=width_x
61			#
62			# create array
63			#
64			receiver_line=[2absorption_zone + i (rangeRcv//(numRcvPerLine-1)) for i in range(numRcvPerLine) ]
65			#
66			# set source location with tag "source""
67			#
68			src_tags=["source"]
69			if DIM == 2:
70			src_locations = [ (receiver_line[srcEW], depth)]
71			src_loc_2D=(receiver_line[srcEW], 0.)
72			else:
73			if srcEW:
74			srcNS=numRcvPerLine//2
75			elif srcNS:
76			srcEW=numRcvPerLine//2
77			else:
78			raise ValueError("on of the variables srcEW or srcNS must be None!")
79			src_locations = [ (receiver_line[srcEW], receiver_line[srcNS], depth)]
80			src_loc_2D=(receiver_line[srcEW], receiver_line[srcNS])
81			#
82			# create sensor arrays:
83			#
84			# East-west line of receiver
85			rcvEW_locations=[]
86			rgEW=[]
87			mid_point=receiver_line[len(receiver_line)//2]
88	gross	4553
89	caltinay	6119	for ix in range(len(receiver_line)):
90			if DIM == 2:
91			rcvEW_locations.append((receiver_line[ix], depth))
92			rgEW.append( ( receiver_line[ix], 0.) )
93			else:
94			rcvEW_locations.append((receiver_line[ix], mid_point, depth))
95			rgEW.append( ( receiver_line[ix], mid_point) )
96			# North-south line of receiver
97			if DIM == 3:
98			rcvNS_locations=[]
99			rgNS=[]
100	gross	4553
101	caltinay	6119	for iy in range(len(receiver_line)):
102			rcvNS_locations.append((mid_point, receiver_line[iy], depth))
103			rgNS.append( ( mid_point, receiver_line[iy]) )
104			#
105			# create domain:
106			#
107			order = 5
108			if DIM == 2:
109			domain=Rectangle(order, ceil(ne_z*width_x/depth),ne_z,l0=width_x,l1=depth,
110			diracPoints=src_locations, diracTags=src_tags)
111			else:
112			domain=Brick(order, ceil(ne_zwidth_x/depth), ceil(ne_zwidth_y/depth),
113			ne_z, l0=width_x, l1=width_y, l2=depth,
114			diracPoints=src_locations, diracTags=src_tags)
115			wl=Ricker(frq)
116	gross	4548
117	caltinay	6119	#======================================================================
118			# m=whereNegative(Function(domain).getX()[DIM-1]-reflector_at)
119			# v_p=v_p_bottomm+v_p_top(1-m)
120			v_p=2*U.km/U.sec
121			v_s=0.9*U.km/U.sec
122			vareps=0.1*0
123			gamma=0.15*0
124			delta=0.05*0
125			rho=2000U.kg/U.m*3
126			src_dir=[0,0,1]
127	gross	4548
128	caltinay	6119	sw=VTIWave(domain, v_p, v_s, wl, src_tags[0], source_vector = src_dir,
129			eps=vareps, gamma=gamma, delta=delta, rho=rho,
130			absorption_zone=300*U.m, absorption_cut=1e-2, lumping=True)
131	gross	4562
132	caltinay	6119	locEW=Locator(domain,rcvEW_locations)
133			tracerEW_x=SimpleSEGYWriter(receiver_group=rgEW, source=src_loc_2D, sampling_interval=sampling_interval, text='x-displacement - east-west line')
134			tracerEW_z=SimpleSEGYWriter(receiver_group=rgEW, source=src_loc_2D, sampling_interval=sampling_interval, text='z-displacement - east-west line')
135			if DIM==3:
136			locNS=Locator(domain,rcvNS_locations)
137			tracerEW_y=SimpleSEGYWriter(receiver_group=rgEW, source=src_loc_2D,
138			sampling_interval=sampling_interval,
139			text='x-displacement - east-west line')
140			tracerNS_x=SimpleSEGYWriter(receiver_group=rgNS, source=src_loc_2D,
141			sampling_interval=sampling_interval,
142			text='x-displacement - north-south line')
143			tracerNS_y=SimpleSEGYWriter(receiver_group=rgNS, source=src_loc_2D,
144			sampling_interval=sampling_interval,
145			text='y-displacement - north-south line')
146			tracerNS_z=SimpleSEGYWriter(receiver_group=rgNS, source=src_loc_2D,
147			sampling_interval=sampling_interval,
148			text='z-displacement - north-south line')
149			if not tracerEW_x.obspy_available():
150			print("\nWARNING: obspy not available, SEGY files will not be written\n")
151			elif getMPISizeWorld() > 1:
152			print("\nWARNING: SEGY files cannot be written with multiple processes\n")
153	gross	4548
154	caltinay	6119	t=0.
155			mkDir('tmp')
156			n=0
157			while t < t_end:
158			t,u = sw.update(t+sampling_interval)
159			tracerEW_x.addRecord(locEW(u[0]))
160			tracerEW_z.addRecord(locEW(u[DIM-1]))
161			if DIM==3:
162			tracerEW_y.addRecord(locEW(u[1]))
163			tracerNS_x.addRecord(locNS(u[0]))
164			tracerNS_y.addRecord(locNS(u[1]))
165			tracerNS_z.addRecord(locNS(u[2]))
166			print(t, locEW(u[DIM-1])[len(rgEW)//2-4:len(rgEW)//2+1], wl.getValue(t))
167			#if n%5 == 0 : saveSilo("tmp/u_%d.silo"%(n/5,), u=u)
168	aellery	6812	saveSilo("tmp/u_%d.silo"%(n,), u=u, cycle=n, time=t)
169	caltinay	6119	n+=1
170			if tracerEW_x.obspy_available() and getMPISizeWorld() == 1:
171			tracerEW_x.write('lineEW_x.sgy')
172			tracerEW_z.write('lineEW_z.sgy')
173			if DIM == 3:
174			tracerEW_y.write('lineEW_y.sgy')
175			tracerNS_x.write('lineNS_x.sgy')
176			tracerNS_y.write('lineNS_y.sgy')
177			tracerNS_z.write('lineNS_z.sgy')
178
179			else: # no speckley
180			print("The Speckley module is not available")
181