Commit ba6990a6 authored by Alex Breuer's avatar Alex Breuer

Added comparison of forward and SGT-derived solution of the SGT test case.

parent 18b414dc
[![CC0](http://i.creativecommons.org/p/zero/1.0/88x31.png)](http://creativecommons.org/publicdomain/zero/1.0/)
To the extent possible under law, [Alexander Nikolas Breuer](http://dial3343.org) has waived all copyright and related or neighboring rights to the following work:
* `forward.tar.bz2`
* `point_sources.tar.bz2`
* `reciprocal.tar.bz2`
* `stream.pdf`
This work is published from: United States.
#!/bin/bash
##
# @file This file is part of EDGE.
#
# @author Alexander Breuer (anbreuer AT ucsd.edu)
#
# @section LICENSE
# Copyright (c) 2018, Regents of the University of California
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
#
# 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
#
# 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
# @section DESCRIPTION
# Generates the reciprocal solutions.
##
STATIONS=(VC2E MFU DFU FFU UPSAR01 FZ8 GH3E SC3E )
IDS=(1 3 5 7 9 11 13 15)
for run in `seq 0 7`
do
id=${IDS[$run]}
station=${STATIONS[$run]}
echo "processing station with id ${id}: " ${STATIONS[$run]}
python3 ./sgt.py --lambda 17136000000 \
--mu 4032000000 \
--source_deconv \
../../sources/gaussian.ascii \
--source_conv ../../sources/ricker.ascii \
--source_resample 100 \
--receiver point_sources/hypocenter.csv \
--sim_id ${id} \
--out reciprocal/${station}.csv
done
#!/usr/bin/python3
##
# @file This file is part of EDGE.
#
# @author Alexander Breuer (anbreuer AT ucsd.edu)
#
# @section LICENSE
# Copyright (c) 2018, Regents of the University of California
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
#
# 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
#
# 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
# @section DESCRIPTION
# Computes the Strain Green's Tensors
##
import logging
import argparse
import numpy
import pandas
def strainMat2D( i_lam,
i_mu ):
"""Computes the the numpy matrix, which is required for converting stresses to strain.
Arguments:
i_lam {float} -- Lame parameter lambda.
i_mu {float} -- Lame parameter mu.
"""
l_mat = numpy.zeros( (3,3) )
l_sca = -i_lam / ( 2*i_mu * ( 3*i_lam + 2*i_mu ) )
l_mat[0,0] = l_sca
l_mat[0,1] = l_sca
l_mat[1,0] = l_sca
l_mat[1,1] = l_sca
l_sca = 1 / ( 2*i_mu )
l_mat[0,0] = l_mat[0,0] + l_sca
l_mat[1,1] = l_mat[1,1] + l_sca
l_mat[2,2] = l_mat[2,2] + l_sca
return l_mat
def strainMat3D( i_lam,
i_mu ):
"""Computes the the numpy matrix, which is required for converting stresses to strain.
Arguments:
i_lam {float} -- Lame parameter lambda.
i_mu {float} -- Lame parameter mu.
"""
l_mat = numpy.zeros( (6,6) )
l_yMod = i_mu * ( 3*i_lam + 2*i_mu )
l_yMod = l_yMod / ( i_lam + i_mu )
l_pRat = i_lam / (2*i_lam + 2*i_mu)
l_mat[0,0] = 1
l_mat[0,1] = -l_pRat
l_mat[0,2] = -l_pRat
l_mat[1,0] = -l_pRat
l_mat[1,1] = 1
l_mat[1,2] = -l_pRat
l_mat[2,0] = -l_pRat
l_mat[2,1] = -l_pRat
l_mat[2,2] = 1
l_mat[3,3] = (1 + l_pRat)
l_mat[4,4] = (1 + l_pRat)
l_mat[5,5] = (1 + l_pRat)
l_mat = l_mat * 1 / l_yMod
return l_mat
# set up logger
logging.basicConfig( level=logging.INFO,
format='%(asctime)s - %(name)s - %(levelname)s - %(message)s' )
# command line arguments
l_parser = argparse.ArgumentParser( description='Computes Strain Green\'s Tensors, deconvolves the synthetics and convolves them with a new source.' )
l_parser.add_argument( '--lambda',
dest = 'lambda',
required = True,
type = float,
help = 'Lame parameter lambda' )
l_parser.add_argument( '--mu',
dest = 'mu',
required = True,
type = float,
help = 'Lame parameter mu' )
l_parser.add_argument( '--source_deconv',
dest = 'source_deconv',
required = True,
type = str,
help = 'source, which will be deconvolved from the synthetics' )
l_parser.add_argument( '--sim_id',
dest = 'sim_id',
required = True,
type = int,
help = 'id of the simulation for which the receiver data is extracted.' )
l_parser.add_argument( '--source_conv',
dest = 'source_conv',
required = True,
type = str,
help = 'source, which will be convolved to the synthetics' )
l_parser.add_argument( '--source_resample',
dest = 'source_resample',
required = True,
type = int,
help = 'Resampling of the source time function. 1 takes the given source, 2 every other value, etc' )
l_parser.add_argument( '--receiver',
dest = 'receiver',
required = True,
type = str,
help = 'receiver for which the SGT is computed' )
l_parser.add_argument( '--out',
dest = 'out',
required = True,
type = str,
help = 'output file for the SGT-derived seismograms' )
l_args = vars(l_parser.parse_args())
logging.info( 'Computing SGTs.' )
# read the CSV
l_sim = str(l_args['sim_id'])
l_stresses = ['Q0', 'Q1', 'Q2', 'Q3', 'Q4', 'Q5']
for l_qt in range(6):
l_stresses[l_qt] = l_stresses[l_qt] + '_C' + l_sim
l_recv = pandas.read_csv( l_args['receiver'],
delim_whitespace=False )
# derive the strain tensor
l_mat = strainMat3D( l_args['lambda'],
l_args['mu'] )
l_strain = l_mat * numpy.matrix( l_recv[[*l_stresses]][:].transpose() )
# read the source of the forward simulation
l_srcDec = pandas.read_csv( l_args['source_deconv'],
skiprows=0,
delim_whitespace=False )
l_srcDec = numpy.array( l_srcDec ).flatten()
l_srcDec = numpy.squeeze( l_srcDec[0:-1:l_args['source_resample']] )
# read the new source, which gets convolved
l_srcCon = pandas.read_csv( l_args['source_conv'],
skiprows=0,
delim_whitespace=False )
l_srcCon = numpy.array( l_srcCon ).flatten()
l_srcCon = l_srcCon
l_srcCon = numpy.squeeze( l_srcCon[0:-1:l_args['source_resample']] )
# transfer sources to frequency space
l_fftSrcDec = numpy.zeros( l_strain.shape[1] )
l_fftSrcCon = numpy.zeros( l_strain.shape[1] )
l_fftSrcDec = numpy.fft.fft( l_srcDec, n=l_strain.shape[1] )
l_fftSrcCon = numpy.fft.fft( l_srcCon, n=l_strain.shape[1] )
# parse the source signal
l_signalSrc = numpy.asarray( l_strain )
# transfer signal to frequency space
l_fftSignalSrc = numpy.zeros( l_signalSrc.shape )
l_fftSignalSrc = numpy.fft.fft( l_signalSrc )
# do the deconvolution, convolution in frequency space and transfer back to time domain
l_signalDes = numpy.fft.ifft( l_fftSrcCon * (l_fftSignalSrc / l_fftSrcDec) )
# write data
print( 'printing: ' + l_args['out'] )
l_signalDf = pandas.DataFrame( data=numpy.real( -2*(1E12/2800)*l_signalDes[3] ) )
l_signalDf.to_csv( l_args['out'], mode='w', index=False )
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
#!/usr/bin/python3
##
# @file This file is part of EDGE.
#
# @author Alexander Breuer (anbreuer AT ucsd.edu)
#
# @section LICENSE
# Copyright (c) 2018, Regents of the University of California
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
#
# 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
#
# 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
# @section DESCRIPTION
# Visualizes the SGT solution in comparison to double couple.
##
###
import pandas
import numpy
import math
import obspy
l_stations = pandas.read_csv( 'stations.csv',
header=None,
delim_whitespace=False )
l_stations = numpy.array( l_stations )
# get traces
l_traces = []
for l_st in l_stations:
l_dist = l_st[1]*l_st[1] + l_st[2]*l_st[2] + (l_st[3]+7622.4) * (l_st[3]+7622.4)
l_dist = math.sqrt( l_dist )
print( 'processing station: ' + l_st[0] + ' at distance ' + str(l_dist/1000) + ' km' )
l_forward = pandas.read_csv( 'forward/'+l_st[0]+'.csv',
delim_whitespace=False )
l_forward = numpy.squeeze( numpy.array( l_forward[['Q7_C0']] ) )
l_reciprocal = pandas.read_csv( 'reciprocal/'+l_st[0]+'.csv',
delim_whitespace=False )
l_reciprocal = numpy.squeeze( numpy.array( l_reciprocal ) )
l_traces.append( obspy.Trace( l_forward ) )
l_traces.append( obspy.Trace( l_reciprocal ) )
l_traces[-2].stats.distance = l_dist
l_traces[-2].stats.delta = 0.01
l_traces[-2].stats.network = 'Forward simulation with source at (0,0,-7622.4)'
l_traces[-1].stats.distance = l_dist
l_traces[-1].stats.delta = 0.01
l_traces[-1].stats.station = l_st[0]
l_traces[-1].stats.network = 'Reciprocal Strain Green\'s Tensor solution'
# plot solution
l_stream = obspy.Stream( traces=l_traces )
l_stream.plot( type='section', color='network', linewidth=1, scale=0.8, outfile='stream.pdf' )
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment