[cig-commits] [commit] master: Added sources obtained by reading sac-files (49ba0d2)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Dec 9 15:16:06 PST 2014
Repository : https://github.com/geodynamics/sw4
On branch : master
Link : https://github.com/geodynamics/sw4/compare/7a86b87d6a91ee96d2cf83382bac526772c949c4...82097c4d1dd428ac74d19be8aed92737fbef1246
>---------------------------------------------------------------
commit 49ba0d2528cb42f44d8fd008e6155612a6c82950
Author: Bjorn Sjogreen <sjogreen2 at llnl.gov>
Date: Fri Dec 5 16:07:12 2014 -0800
Added sources obtained by reading sac-files
>---------------------------------------------------------------
49ba0d2528cb42f44d8fd008e6155612a6c82950
Makefile | 2 +-
src/EW.h | 2 +-
src/GridPointSource.C | 90 +++++++++-
src/Source.C | 34 +++-
src/TimeDep.h | 2 +-
src/bcfort.f | 2 +-
src/lamb_one_point.f | 462 --------------------------------------------------
src/main.C | 2 +-
src/parallelStuff.C | 2 +-
src/parseInputFile.C | 94 +++++++++-
src/sacutils.C | 190 +++++++++++++++++++++
src/sacutils.h | 15 ++
12 files changed, 417 insertions(+), 480 deletions(-)
diff --git a/Makefile b/Makefile
index b1d922a..f959187 100644
--- a/Makefile
+++ b/Makefile
@@ -157,7 +157,7 @@ OBJ = EW.o Sarray.o version.o parseInputFile.o ForcingTwilight.o \
lamb_exact_numquad.o twilightsgfort.o EtreeFile.o MaterialIfile.o GeographicProjection.o \
rhs4curvilinear.o curvilinear4.o rhs4curvilinearsg.o curvilinear4sg.o gradients.o Image3D.o \
MaterialVolimagefile.o MaterialRfile.o randomfield3d.o innerloop-ani-sgstr-vc.o bcfortanisg.o \
- AnisotropicMaterialBlock.o checkanisomtrl.o computedtaniso.o
+ AnisotropicMaterialBlock.o checkanisomtrl.o computedtaniso.o sacutils.o
# prefix object files with build directory
FSW4 = $(addprefix $(builddir)/,$(OBJSW4))
diff --git a/src/EW.h b/src/EW.h
index e661f6d..c3a41f8 100644
--- a/src/EW.h
+++ b/src/EW.h
@@ -689,7 +689,7 @@ vector<Sarray> mC; // Anisotropic material parameters
private:
void preprocessSources( vector<Source*> & a_GlobalSources );
-
+void revvector( int npts, double* v );
// epicenter
double m_epi_lat, m_epi_lon, m_epi_depth, m_epi_t0;
diff --git a/src/GridPointSource.C b/src/GridPointSource.C
index 8729281..0f5070d 100644
--- a/src/GridPointSource.C
+++ b/src/GridPointSource.C
@@ -235,6 +235,14 @@ void GridPointSource::initializeTimeFunction()
mTimeFunc_om = Discrete_om;
mTimeFunc_omtt = Discrete_omtt;
break;
+ case iDiscrete6moments :
+ mTimeFunc = Discrete;
+ mTimeFunc_t = Discrete_t;
+ mTimeFunc_tt = Discrete_tt;
+ mTimeFunc_ttt = Discrete_ttt;
+ mTimeFunc_om = Discrete_om;
+ mTimeFunc_omtt = Discrete_omtt;
+ break;
case iC6SmoothBump :
mTimeFunc = C6SmoothBump;
mTimeFunc_t = C6SmoothBump_t;
@@ -283,6 +291,13 @@ void GridPointSource::initializeTimeFunction()
mTimeFunc_tom = Discrete_tom;
mTimeFunc_omom = Discrete_omom;
break;
+ case iDiscrete6moments :
+ mTimeFunc_tttt = Discrete_tttt;
+ mTimeFunc_tttom = Discrete_tttom;
+ mTimeFunc_ttomom = Discrete_ttomom;
+ mTimeFunc_tom = Discrete_tom;
+ mTimeFunc_omom = Discrete_omom;
+ break;
default:
// tmp
// std::cout << "High derivatives not implemented for time fuction:" << mTimeDependence <<
@@ -298,12 +313,41 @@ void GridPointSource::initializeTimeFunction()
//-----------------------------------------------------------------------
void GridPointSource::getFxyz( double t, double* fxyz ) const
{
- double afun = mTimeFunc(mFreq,t-mT0,mPar, mNpar, mIpar, mNipar );
+ double afun, afunv[6];
+ if( mTimeDependence != iDiscrete6moments )
+ afun= mTimeFunc(mFreq,t-mT0,mPar, mNpar, mIpar, mNipar );
+ else
+ {
+ int npts = mIpar[0];
+ int size = 6*(npts-1)+1;
+ size_t pos = 0;
+ afunv[0] = mTimeFunc(mFreq,t-mT0,mPar+pos, mNpar, mIpar, mNipar );
+ pos += size;
+ afunv[1] = mTimeFunc(mFreq,t-mT0,mPar+pos, mNpar, mIpar, mNipar );
+ pos += size;
+ afunv[2] = mTimeFunc(mFreq,t-mT0,mPar+pos, mNpar, mIpar, mNipar );
+ pos += size;
+ afunv[3] = mTimeFunc(mFreq,t-mT0,mPar+pos, mNpar, mIpar, mNipar );
+ pos += size;
+ afunv[4] = mTimeFunc(mFreq,t-mT0,mPar+pos, mNpar, mIpar, mNipar );
+ pos += size;
+ afunv[5] = mTimeFunc(mFreq,t-mT0,mPar+pos, mNpar, mIpar, mNipar );
+ }
+
if( m_derivative==-1)
{
- fxyz[0] = mForces[0]*afun;
- fxyz[1] = mForces[1]*afun;
- fxyz[2] = mForces[2]*afun;
+ if( mTimeDependence != iDiscrete6moments )
+ {
+ fxyz[0] = mForces[0]*afun;
+ fxyz[1] = mForces[1]*afun;
+ fxyz[2] = mForces[2]*afun;
+ }
+ else
+ {
+ fxyz[0] = mForces[0]*afunv[0]+mForces[1]*afunv[1]+mForces[2]*afunv[2];
+ fxyz[1] = mForces[0]*afunv[1]+mForces[1]*afunv[3]+mForces[2]*afunv[4];
+ fxyz[2] = mForces[0]*afunv[2]+mForces[1]*afunv[4]+mForces[2]*afunv[5];
+ }
}
else if( m_derivative >= 0 && m_derivative <= 8 )
{
@@ -360,12 +404,42 @@ void GridPointSource::getFxyz_notime( double* fxyz ) const
//-----------------------------------------------------------------------
void GridPointSource::getFxyztt( double t, double* fxyz ) const
{
- double afun = mTimeFunc_tt(mFreq,t-mT0,mPar, mNpar, mIpar, mNipar);
+ double afun, afunv[6];
+ if( mTimeDependence != iDiscrete6moments )
+ afun= mTimeFunc_tt(mFreq,t-mT0,mPar, mNpar, mIpar, mNipar );
+ else
+ {
+ int npts = mIpar[0];
+ int size = 6*(npts-1)+1;
+ size_t pos = 0;
+ afunv[0] = mTimeFunc_tt(mFreq,t-mT0,mPar+pos, mNpar, mIpar, mNipar );
+ pos += size;
+ afunv[1] = mTimeFunc_tt(mFreq,t-mT0,mPar+pos, mNpar, mIpar, mNipar );
+ pos += size;
+ afunv[2] = mTimeFunc_tt(mFreq,t-mT0,mPar+pos, mNpar, mIpar, mNipar );
+ pos += size;
+ afunv[3] = mTimeFunc_tt(mFreq,t-mT0,mPar+pos, mNpar, mIpar, mNipar );
+ pos += size;
+ afunv[4] = mTimeFunc_tt(mFreq,t-mT0,mPar+pos, mNpar, mIpar, mNipar );
+ pos += size;
+ afunv[5] = mTimeFunc_tt(mFreq,t-mT0,mPar+pos, mNpar, mIpar, mNipar );
+ }
+
+ // double afun = mTimeFunc_tt(mFreq,t-mT0,mPar, mNpar, mIpar, mNipar);
if( m_derivative==-1)
{
- fxyz[0] = mForces[0]*afun;
- fxyz[1] = mForces[1]*afun;
- fxyz[2] = mForces[2]*afun;
+ if( mTimeDependence != iDiscrete6moments )
+ {
+ fxyz[0] = mForces[0]*afun;
+ fxyz[1] = mForces[1]*afun;
+ fxyz[2] = mForces[2]*afun;
+ }
+ else
+ {
+ fxyz[0] = mForces[0]*afunv[0]+mForces[1]*afunv[1]+mForces[2]*afunv[2];
+ fxyz[1] = mForces[0]*afunv[1]+mForces[1]*afunv[3]+mForces[2]*afunv[4];
+ fxyz[2] = mForces[0]*afunv[2]+mForces[1]*afunv[4]+mForces[2]*afunv[5];
+ }
}
else if( m_derivative >= 0 && m_derivative <= 8 )
{
diff --git a/src/Source.C b/src/Source.C
index 439dcf3..6d76457 100644
--- a/src/Source.C
+++ b/src/Source.C
@@ -127,7 +127,7 @@ Source::Source(EW *a_ew,
mIpar = new int[1];
}
- if( mTimeDependence == iDiscrete )
+ if( mTimeDependence == iDiscrete || mTimeDependence == iDiscrete6moments )
spline_interpolation();
else
{
@@ -198,7 +198,7 @@ Source::Source(EW *a_ew, double frequency, double t0,
mNipar = 1;
mIpar = new int[1];
}
- if( mTimeDependence == iDiscrete )
+ if( mTimeDependence == iDiscrete || mTimeDependence == iDiscrete6moments )
spline_interpolation();
else
{
@@ -1712,6 +1712,36 @@ int Source::spline_interpolation( )
return 1;
}
+ else if( mTimeDependence == iDiscrete6moments )
+ {
+ int npts = mIpar[0];
+ // Six different time function, one for each momentum component.
+ // Store sequentially in mPar, i.e.,
+ // mPar = [tstart, first time func, tstart, second time func, ...]
+ // tstart before each function ---> Can reuse time function iDiscrete.
+
+ double* parin = new double[(npts+1)*6];
+ for( int i=0 ; i < (npts+1)*6; i++ )
+ parin[i] = mPar[i];
+ double tstart = mPar[0];
+ delete[] mPar;
+ mNpar = 6*(6*(npts-1)+1);
+ mPar = new double[mNpar];
+
+ size_t pos_in=0, pos_out=0;
+ for( int tf=0 ; tf < 6 ; tf++ )
+ {
+ Qspline quinticspline( npts, &parin[pos_in+1], tstart, 1/mFreq );
+ pos_in += npts+1;
+ mPar[pos_out] = tstart;
+ double* qsppt = quinticspline.get_polycof_ptr();
+ for( int i=0 ; i < 6*(npts-1) ; i++ )
+ mPar[pos_out+i+1] = qsppt[i];
+ pos_out += 6*(npts-1)+1;
+ }
+ delete[] parin;
+ return 1;
+ }
else
return 0;
}
diff --git a/src/TimeDep.h b/src/TimeDep.h
index 39dd3a4..1d693d3 100644
--- a/src/TimeDep.h
+++ b/src/TimeDep.h
@@ -33,6 +33,6 @@
#ifndef TIMEDEP_H
#define TIMEDEP_H
-enum timeDep { iRicker, iGaussian, iRamp, iTriangle, iSawtooth, iSmoothWave, iErf, iRickerInt, iVerySmoothBump, iBrune, iBruneSmoothed, iGaussianWindow, iLiu, iDiscrete, iDBrune, iDirac, iC6SmoothBump };
+enum timeDep { iRicker, iGaussian, iRamp, iTriangle, iSawtooth, iSmoothWave, iErf, iRickerInt, iVerySmoothBump, iBrune, iBruneSmoothed, iGaussianWindow, iLiu, iDiscrete, iDBrune, iDirac, iC6SmoothBump, iDiscrete6moments };
#endif
diff --git a/src/bcfort.f b/src/bcfort.f
index 8bde97a..321dd43 100644
--- a/src/bcfort.f
+++ b/src/bcfort.f
@@ -1100,7 +1100,7 @@ c-----------------------------------------------------------------------
integer i, j, kz
real*8 muu, lambdaa, x, y, z, t, om, c, ph
doubleprecision mu(ifirst:ilast,jfirst:jlast,kfirst:klast)
- doubleprecision lambda(ifirst:ilast,jfirst:jlast,kfirst:klast) doubleprecision x
+ doubleprecision lambda(ifirst:ilast,jfirst:jlast,kfirst:klast)
doubleprecision omstrx
doubleprecision omstry
diff --git a/src/lamb_one_point.f b/src/lamb_one_point.f
deleted file mode 100644
index 7f6d3d0..0000000
--- a/src/lamb_one_point.f
+++ /dev/null
@@ -1,462 +0,0 @@
-! SW4 LICENSE
-! # ----------------------------------------------------------------------
-! # SW4 - Seismic Waves, 4th order
-! # ----------------------------------------------------------------------
-! # Copyright (c) 2013, Lawrence Livermore National Security, LLC.
-! # Produced at the Lawrence Livermore National Laboratory.
-! #
-! # Written by:
-! # N. Anders Petersson (petersson1 at llnl.gov)
-! # Bjorn Sjogreen (sjogreen2 at llnl.gov)
-! #
-! # LLNL-CODE-643337
-! #
-! # All rights reserved.
-! #
-! # This file is part of SW4, Version: 1.0
-! #
-! # Please also read LICENCE.txt, which contains "Our Notice and GNU General Public License"
-! #
-! # This program is free software; you can redistribute it and/or modify
-! # it under the terms of the GNU General Public License (as published by
-! # the Free Software Foundation) version 2, dated June 1991.
-! #
-! # This program is distributed in the hope that it will be useful, but
-! # WITHOUT ANY WARRANTY; without even the IMPLIED WARRANTY OF
-! # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the terms and
-! # conditions of the GNU General Public License for more details.
-! #
-! # You should have received a copy of the GNU General Public License
-! # along with this program; if not, write to the Free Software
-! # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
- program lambmain
- implicit none
- real*8 t, r, uex3, mu, cs, fz, rho, tmax, dt, rocs
- integer tfun, k, nt, i
- character(360) buf
- common /funpars/ t, rocs
-
-c VerySmoothBump with tfun=1
-c C6SmoothBump with tfun=2
- tfun = 1
- cs = 1
- rho = 1.5
- mu = cs*cs*rho
- fz = 1d0
-c max time and number of time step
- tmax=30.0
-c nt = 1032
-c nt = 1677
- nt = 4695
-
-c source to receiver distance
- r=10d0
-
-c read command line arguments
- i = 1
- do while( i .lt. IARGC() )
- call GETARG( i, buf )
- if( buf .eq. '-tmax' )then
- call GETARG(i+1,buf)
- read(buf,'(f16.0)') tmax
- i = i+2
- elseif( buf.eq.'-nsteps' )then
- call GETARG(i+1,buf)
- read(buf,'(i10)') nt
- i = i+2
- elseif( buf.eq.'-dist' )then
- call GETARG(i+1,buf)
- read(buf,'(f16.0)') r
- i = i+2
- elseif( buf.eq.'-tfunc' )then
- call GETARG(i+1,buf)
- read(buf,'(i10)') tfun
- i = i+2
- else
- i = i+ 1
- endif
- enddo
-
- write(*,100)'Fz: ', fz
- 100 format(' ', a, es12.5)
- write(*,101)'Shear speed: ',cs,' Density: ',rho,
- + ' Time function: ', tfun,' (1=VerySmoothBump, 2=C6SmoothBmp)'
- write(*,101)'Source-rec dist: ', r,' Max time: ',tmax,
- + ' # time steps: ', nt
- 101 format(' ', a, es12.5, a, es12.5, a, i7, a)
-
-c basic checks
- if (.not.(tfun .eq. 1 .or. tfun .eq. 2)) then
- write(*,*)'Unknown time function, tfun=', tfun
- stop
- endif
-
- if (.not.(nt .gt. 0)) then
- write(*,*)'# time steps must be positive, not:', nt
- stop
- endif
-
- if (.not.(tmax .gt. 0)) then
- write(*,*)'End time must be positive, not:', tmax
- stop
- endif
-
- if (.not.(r .gt. 0)) then
- write(*,*)'src-rec distance must be positive, not:', r
- stop
- endif
-
-
-c name of file
- open(10,file='uzex.dat',status='unknown')
-
-c should not need to change anything beyond this point
- dt = tmax/(nt)
- rocs = r/cs
-
- do k=0,nt
- t = dt*k
- call LAMBONEPOINT(r,uex3, mu, cs, fz, tfun)
- write(10,102) t, uex3
- 102 format(' ',e23.16,e23.16)
-c testing
-c write(*,*) t, uex3
- enddo
- close(10)
- end
-
- subroutine LAMBONEPOINT(r, uex3, mu, cs, fz, tfun )
- implicit none
-
- external G1FUN, G2FUN, G2FUNNW
- external G1FUNS, G2FUNS, G2FUNNWS
-
- real*8 cpocs, gamma, gammac, pi, c21, c22, c12, c11, uex3
- parameter( cpocs=1.73205080756888d0 )
- parameter( gamma=1.08766387358054d0 )
- parameter( gammac=0.563016250305247d0 )
- parameter( pi= 3.14159265358979d0 )
-*** pi/8
- parameter( c21=0.392699081698724d0 )
-*** pi/48*sqrt(3*sqrt(3)+5)
- parameter( c22=0.208990620247103d0 )
-*** pi/96*sqrt(3)
- parameter( c12=0.0566812301323193d0)
-*** pi/96*sqrt(3*sqrt(3)-5)
- parameter( c11=0.0144935735221071d0)
-
- real*8 G1FUN, G2FUN, G2FUNNW, INTFCN3
- real*8 G1FUNS, G2FUNS, G2FUNNWS, INTFCN3S
-
- integer tfun
- real*8 t, x0, y0, cs, uzex, r, cp
- real*8 a1, b1, a2, b2, a3, b3, al1, al2, al3
- real*8 fz, mu, csor, rocs, tim
-
- integer neval, ier, limit, lenw, last, integr
- integer, allocatable, dimension(:) :: iwork
- real*8 a, b, epsabs, epsrel, result
- real*8 abserr, exact, alfa, beta
- real*8, allocatable, dimension(:) :: work
-
- common /funpars/ t, rocs
-c save /funpars/
-
- if( tfun.ne.1 .and. tfun.ne.2 )then
- write(*,*) 'ERROR LAMBEXACTNUMQUAD_ only implemented '//
- * ' for VerySmoothBump and C6SmoothBump '
- return
- endif
-c now passing t through the common block
-c tim = t
- limit = 100
- lenw = 4*limit
- epsabs = 1d-12
- epsrel = 1d-12
- allocate( iwork(limit+10), work(lenw+10))
-
-c$$$ do j=jfirst,jlast
-c$$$ do i=ifirst,ilast
-c$$$ x = (i-1)*h
-c$$$ y = (j-1)*h
-c$$$ r = SQRT( (x-x0)*(x-x0)+(y-y0)*(y-y0) )
- cp = cpocs*cs
- if( r.ge.cp*t )then
- uex3 = 0
- elseif( r.ge.1d-10 )then
- csor = cs/r
-c rocs = r/cs
- a1 = MAX(csor*(t-1),1/cpocs)
- b1 = MIN(csor*t,1d0)
- a2 = MAX(csor*(t-1),1d0)
- b2 = MIN(csor*t,gamma)
- a3 = MAX(csor*(t-1),gamma)
- b3 = csor*t
- uzex = 0
- if( b1.gt.a1+1e-12 )then
-c write(*,*) 'Doing 1'
- if( tfun.eq.1 )then
- call DQAGS( G1FUN, a1, b1, epsabs, epsrel, result,
- * abserr, neval, ier, limit, lenw, last, iwork, work )
- else
- call DQAGS( G1FUNS, a1, b1, epsabs, epsrel, result,
- * abserr, neval, ier, limit, lenw, last, iwork, work )
- endif
- if( ier.ne.0 )then
- write(*,*)'ERROR: call to DQAGS returns ier = ',ier
- endif
- if( abserr.gt.10*epsabs )then
- write(*,*) 'WARNING: DQAGS returns error ',abserr,
- * ' input tolerance ', epsabs
- endif
- if( tfun.eq.1 )then
- uzex = rocs*result +
- * 0.5d0*c21*INTFCN3(t-a1*rocs,t-b1*rocs)
- else
- uzex = rocs*result +
- * 0.5d0*c21*INTFCN3S(t-a1*rocs,t-b1*rocs)
- endif
- endif
- if( b2.gt.a2+1e-12 )then
-c write(*,*) 'Doing 2'
- alfa = 0
- beta = -0.5d0
- integr = 1
- if( abs(b2-gamma).lt.1d-12 )then
-*** Singular integral
- if( tfun.eq.1 )then
- call DQAWS( G2FUN, a2, b2, alfa, beta, integr,
- * epsabs, epsrel, result, abserr, neval, ier, limit,
- * lenw, last, iwork, work )
- else
- call DQAWS( G2FUNS, a2, b2, alfa, beta, integr,
- * epsabs, epsrel, result, abserr, neval, ier, limit,
- * lenw, last, iwork, work )
- endif
- else
-*** Away from singularity
- if( tfun.eq.1 )then
- call DQAGS( G2FUNNW, a2, b2, epsabs, epsrel,result,
- * abserr, neval, ier, limit, lenw, last, iwork,
- * work )
- else
- call DQAGS( G2FUNNWS,a2, b2, epsabs, epsrel,result,
- * abserr, neval, ier, limit, lenw, last, iwork,
- * work )
- endif
- endif
- if( ier.ne.0 )then
- write(*,*)'ERROR: call to DQAWS returns ier = ',ier
- endif
- if( abserr.gt.10*epsabs )then
- write(*,*) 'WARNING: DQAWS returns error ',abserr,
- * ' input tolerance ', epsabs
- endif
- uzex = uzex + rocs*result
- if( tfun.eq.1 )then
- uzex = uzex - c21*INTFCN3(t-b2*rocs,t-a2*rocs)
- else
- uzex = uzex - c21*INTFCN3S(t-b2*rocs,t-a2*rocs)
- endif
- endif
- if( b3.gt.a3 )then
-c write(*,*) 'Doing 3'
- if( tfun.eq.1 )then
- uzex = uzex - c21*INTFCN3(0d0,t-a3*rocs)
- else
- uzex = uzex - c21*INTFCN3S(0d0,t-a3*rocs)
- endif
- endif
- if( r.ge.1e-10 )then
- uex3 = -uzex*fz/(pi*pi*mu*r)*cpocs*cpocs
- else
- uex3=0
- endif
- endif
-c$$$ enddo
-c$$$ enddo
-
- deallocate( iwork, work )
-
- return
- end
-
-c-----------------------------------------------------------------------
- real*8 function INTFCN3( a, b )
- implicit none
- real*8 cof
- parameter( cof = 1024d0 )
- real*8 a, b
- INTFCN3 = cof*(b**5*(1-b)**5-a**5*(1-a)**5)
- end
-
-c-----------------------------------------------------------------------
- real*8 function INTFCN3S( a, b )
- implicit none
- real*8 cof
-c parameter( cof = 16384d0 )
- parameter( cof = 51480d0 )
- real*8 a, b
- INTFCN3S = cof*(b**7*(1-b)**7-a**7*(1-a)**7)
- end
-
-c-----------------------------------------------------------------------
- real*8 function G1FUN( x )
- implicit none
- real*8 x
- real*8 cof
- parameter( cof = 1024d0 )
- real*8 c11, c12, c21, c22, gamma, gammac
- parameter( gamma=1.08766387358054d0 )
- parameter( gammac=0.563016250305247d0 )
-*** pi/8
- parameter( c21=0.392699081698724d0 )
-*** pi/48*sqrt(3*sqrt(3)+5)
- parameter( c22=0.208990620247103d0 )
-*** pi/96*sqrt(3)
- parameter( c12=0.0566812301323193d0)
-*** pi/96*sqrt(3*sqrt(3)-5)
- parameter( c11=0.0144935735221071d0)
-
- real*8 t, rocs
- common /funpars/ t, rocs
-c save /funpars/
-
- G1FUN = 5*cof*(t-rocs*x)**4*(1-t+rocs*x)**4*(1-2*(t-rocs*x))*
- * ( 0.5d0*c22/SQRT(gamma**2-x*x)-c11/SQRT(x*x-gammac**2)
- * + c12/sqrt(x*x-0.25d0) )
-
- end
-
-c-----------------------------------------------------------------------
- real*8 function G1FUNS( x )
- implicit none
- real*8 x
- real*8 cof
-c parameter( cof = 16384d0 )
- parameter( cof = 51480d0 )
- real*8 c11, c12, c21, c22, gamma, gammac
- parameter( gamma=1.08766387358054d0 )
- parameter( gammac=0.563016250305247d0 )
-*** pi/8
- parameter( c21=0.392699081698724d0 )
-*** pi/48*sqrt(3*sqrt(3)+5)
- parameter( c22=0.208990620247103d0 )
-*** pi/96*sqrt(3)
- parameter( c12=0.0566812301323193d0)
-*** pi/96*sqrt(3*sqrt(3)-5)
- parameter( c11=0.0144935735221071d0)
-
- real*8 t, rocs
- common /funpars/ t, rocs
-c save /funpars/
-
- G1FUNS = 7*cof*(t-rocs*x)**6*(1-t+rocs*x)**6*(1-2*(t-rocs*x))*
- * ( 0.5d0*c22/SQRT(gamma**2-x*x)-c11/SQRT(x*x-gammac**2)
- * + c12/sqrt(x*x-0.25d0) )
-
- end
-
-c-----------------------------------------------------------------------
- real*8 function G2FUN( x )
- implicit none
- real*8 x
- real*8 cof
- parameter( cof = 1024d0 )
- real*8 c11, c12, c21, c22, gamma, gammac
- parameter( gamma=1.08766387358054d0 )
- parameter( gammac=0.563016250305247d0 )
-*** pi/8
- parameter( c21=0.392699081698724d0 )
-*** pi/48*sqrt(3*sqrt(3)+5)
- parameter( c22=0.208990620247103d0 )
-*** pi/96*sqrt(3)
- parameter( c12=0.0566812301323193d0)
-*** pi/96*sqrt(3*sqrt(3)-5)
- parameter( c11=0.0144935735221071d0)
-
- real*8 t, rocs
- common /funpars/ t, rocs
-c save /funpars/
-
- G2FUN = 5*cof*(t-rocs*x)**4*(1-t+rocs*x)**4*(1-2*(t-rocs*x))*
- * c22/SQRT(gamma+x)
- end
-
-c-----------------------------------------------------------------------
- real*8 function G2FUNS( x )
- implicit none
- real*8 x
- real*8 cof
-c parameter( cof = 16384d0 )
- parameter( cof = 51480d0 )
- real*8 c11, c12, c21, c22, gamma, gammac
- parameter( gamma=1.08766387358054d0 )
- parameter( gammac=0.563016250305247d0 )
-*** pi/8
- parameter( c21=0.392699081698724d0 )
-*** pi/48*sqrt(3*sqrt(3)+5)
- parameter( c22=0.208990620247103d0 )
-*** pi/96*sqrt(3)
- parameter( c12=0.0566812301323193d0)
-*** pi/96*sqrt(3*sqrt(3)-5)
- parameter( c11=0.0144935735221071d0)
-
- real*8 t, rocs
- common /funpars/ t, rocs
-c save /funpars/
-
- G2FUNS = 7*cof*(t-rocs*x)**6*(1-t+rocs*x)**6*(1-2*(t-rocs*x))*
- * c22/SQRT(gamma+x)
- end
-
-c-----------------------------------------------------------------------
- real*8 function G2FUNNW( x )
- implicit none
- real*8 x
- real*8 cof
- parameter( cof = 1024d0 )
- real*8 c11, c12, c21, c22, gamma, gammac
- parameter( gamma=1.08766387358054d0 )
- parameter( gammac=0.563016250305247d0 )
-*** pi/8
- parameter( c21=0.392699081698724d0 )
-*** pi/48*sqrt(3*sqrt(3)+5)
- parameter( c22=0.208990620247103d0 )
-*** pi/96*sqrt(3)
- parameter( c12=0.0566812301323193d0)
-*** pi/96*sqrt(3*sqrt(3)-5)
- parameter( c11=0.0144935735221071d0)
-
- real*8 t, rocs
- common /funpars/ t, rocs
-c save /funpars/
- G2FUNNW = 5*cof*(t-rocs*x)**4*(1-t+rocs*x)**4*(1-2*(t-rocs*x))*
- * c22/SQRT((gamma-x)*(gamma+x))
- end
-
-c-----------------------------------------------------------------------
- real*8 function G2FUNNWS( x )
- implicit none
- real*8 x
- real*8 cof
-c parameter( cof = 16384d0 )
- parameter( cof = 51480d0 )
- real*8 c11, c12, c21, c22, gamma, gammac
- parameter( gamma=1.08766387358054d0 )
- parameter( gammac=0.563016250305247d0 )
-*** pi/8
- parameter( c21=0.392699081698724d0 )
-*** pi/48*sqrt(3*sqrt(3)+5)
- parameter( c22=0.208990620247103d0 )
-*** pi/96*sqrt(3)
- parameter( c12=0.0566812301323193d0)
-*** pi/96*sqrt(3*sqrt(3)-5)
- parameter( c11=0.0144935735221071d0)
-
- real*8 t, rocs
- common /funpars/ t, rocs
-c save /funpars/
- G2FUNNWS = 7*cof*(t-rocs*x)**6*(1-t+rocs*x)**6*(1-2*(t-rocs*x))*
- * c22/SQRT((gamma-x)*(gamma+x))
- end
diff --git a/src/main.C b/src/main.C
index 1634f55..6053ecf 100644
--- a/src/main.C
+++ b/src/main.C
@@ -29,7 +29,7 @@
// # You should have received a copy of the GNU General Public License
// # along with this program; if not, write to the Free Software
// # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
-#include "mpi.h"
+//#include "mpi.h"
#include "EW.h"
diff --git a/src/parallelStuff.C b/src/parallelStuff.C
index 83aafdb..6ed379e 100644
--- a/src/parallelStuff.C
+++ b/src/parallelStuff.C
@@ -413,7 +413,7 @@ void EW::communicate_array_2d( Sarray& u, int g, int k )
//-----------------------------------------------------------------------
void EW::communicate_array_2d_ext( Sarray& u )
{
- REQUIRE2( u.m_nc == 1, "Communicate array 2d ext, only implemented for three-component arrays" );
+ REQUIRE2( u.m_nc == 1, "Communicate array 2d ext, only implemented for one-component arrays" );
int g = mNumberOfGrids-1;
int ie = m_iEnd[g]+m_ext_ghost_points, ib=m_iStart[g]-m_ext_ghost_points;
int je = m_jEnd[g]+m_ext_ghost_points, jb=m_jStart[g]-m_ext_ghost_points;
diff --git a/src/parseInputFile.C b/src/parseInputFile.C
index 07bbc96..f5e57f1 100644
--- a/src/parseInputFile.C
+++ b/src/parseInputFile.C
@@ -47,6 +47,7 @@
#include "TimeSeries.h"
#include "Filter.h"
#include "Image3D.h"
+#include "sacutils.h"
#ifdef ENABLE_OPT
#include "MaterialInvtest.h"
@@ -77,6 +78,16 @@ using namespace std;
// return a;
//}
+void EW::revvector( int npts, double* v )
+{
+ for( int i=0 ; i < npts/2; i++ )
+ {
+ double sl=v[i];
+ v[i]=v[npts-1-i];
+ v[npts-1-i]=sl;
+ }
+}
+
int computeEndGridPoint( double maxval, double dh )
{
// We round up one, so that the end point
@@ -4564,13 +4575,14 @@ void EW::processSource(char* buffer, vector<Source*> & a_GlobalUniqueSources )
bool geoCoordSet = false;
bool strikeDipRake = false;
bool dfileset=false;
+ bool sacbaseset = false;
int ncyc = 0;
bool ncyc_set = false;
timeDep tDep = iRickerInt;
char formstring[100];
- char dfile[100];
+ char dfile[1000];
strcpy(formstring, "Ricker");
@@ -4830,9 +4842,16 @@ void EW::processSource(char* buffer, vector<Source*> & a_GlobalUniqueSources )
else if (startswith("dfile=",token))
{
token += 6;
- strncpy(dfile, token,100);
+ strncpy(dfile, token,1000);
dfileset = true;
}
+ else if (startswith("sacbase=",token))
+ {
+ token += 8;
+ strncpy(dfile, token,1000);
+ sacbaseset = true;
+ isMomentType = 1;
+ }
else
{
badOption("source", token);
@@ -4879,6 +4898,77 @@ void EW::processSource(char* buffer, vector<Source*> & a_GlobalUniqueSources )
// cout << "Read disc source: t0=" << t0 << " dt=" << dt << " npts= " << npts << endl;
fclose(fd);
}
+ if( sacbaseset )
+ {
+ // Read moment tensor components from sac files.
+ // Set constant tensor to identity. m0 can give some scaling.
+ mxx = myy = mzz = 1;
+ mxy = mxz = myz = 0;
+
+ bool timereverse = false; // Set true for testing purpose only, the users want to do this themselves outside SW4.
+
+ tDep = iDiscrete6moments;
+ double dt, t0, latsac, lonsac,cmpazsac, cmpincsac;
+ int utcsac[7], npts;
+ string basename = dfile;
+ string fname;
+ fname = basename + ".xx";
+ bool byteswap;
+ readSACheader( fname.c_str(), dt, t0, latsac, lonsac, cmpazsac, cmpincsac, utcsac, npts, byteswap );
+ if( geoCoordSet )
+ {
+ double laterr = fabs((latsac-lat)/lat);
+ double lonerr = fabs((lonsac-lon)/lon);
+ if( laterr > 1e-6 || lonerr < 1e-6 )
+ {
+ if( proc_zero() )
+ cout << "WARNING in processSource: reading sac files: (lat,lon) location on sac file different from (lat,lon) on command line" << endl;
+ }
+ }
+ freq = 1/dt;
+ npar = 6*(npts+1);
+ nipar = 1;
+ par = new double [npar];
+ ipar = new int[1];
+ ipar[0] = npts;
+ size_t offset = 0;
+ par[offset] = t0;
+ fname = basename + ".xx";
+ readSACdata( fname.c_str(), npts, &par[offset+1], byteswap );
+ if( timereverse )
+ revvector( npts, &par[offset+1]);
+ offset += npts+1;
+ par[offset] = t0;
+ fname = basename + ".xy";
+ readSACdata( fname.c_str(), npts, &par[offset+1], byteswap );
+ if( timereverse )
+ revvector( npts, &par[offset+1]);
+ offset += npts+1;
+ par[offset] = t0;
+ fname = basename + ".xz";
+ readSACdata( fname.c_str(), npts, &par[offset+1], byteswap );
+ if( timereverse )
+ revvector( npts, &par[offset+1]);
+ offset += npts+1;
+ par[offset] = t0;
+ fname = basename + ".yy";
+ readSACdata( fname.c_str(), npts, &par[offset+1], byteswap );
+ if( timereverse )
+ revvector( npts, &par[offset+1]);
+ offset += npts+1;
+ par[offset] = t0;
+ fname = basename + ".yz";
+ readSACdata( fname.c_str(), npts, &par[offset+1], byteswap );
+ if( timereverse )
+ revvector( npts, &par[offset+1]);
+ offset += npts+1;
+ par[offset] = t0;
+ fname = basename + ".zz";
+ readSACdata( fname.c_str(), npts, &par[offset+1], byteswap );
+ if( timereverse )
+ revvector( npts, &par[offset+1]);
+ }
+
// ---------------------------------------------------------------------------
// Regardless of how the location for the source was specified, we are going to
// find the grid points associated with the location. (i.e., assign
diff --git a/src/sacutils.C b/src/sacutils.C
new file mode 100644
index 0000000..72a43ee
--- /dev/null
+++ b/src/sacutils.C
@@ -0,0 +1,190 @@
+#include <iostream>
+#include <cstdio>
+
+#include "Require.h"
+#include "Byteswapper.h"
+
+using namespace std;
+//-----------------------------------------------------------------------
+int lastofmonth( int year, int month )
+{
+ int days;
+ int leapyear=0;
+ leapyear = ( year % 400 == 0 ) || ( (year % 4 == 0) && !(year % 100 == 0) );
+ if( month == 2 )
+ days = 28 + leapyear;
+ else if( month==4 || month==6 || month==9 || month==11 )
+ days = 30;
+ else
+ days = 31;
+ return days;
+}
+
+//-----------------------------------------------------------------------
+void convertjday( int jday, int year, int& day, int& month )
+{
+ if( jday > 0 && jday < 367 )
+ {
+ day = 1;
+ int jd = 1;
+ month = 1;
+ while( jd < jday )
+ {
+ jd++;
+ day++;
+ if( day > lastofmonth(year,month) )
+ {
+ day = 1;
+ month++;
+ }
+ }
+ }
+ else
+ cout << "convertjday: ERROR, jday is outside range " << endl;
+}
+
+//-----------------------------------------------------------------------
+void readSACheader( const char* fname, double& dt, double& t0,
+ double& lat, double& lon, double& cmpaz,
+ double& cmpinc, int utc[7], int& npts,
+ bool& need_byte_reversal )
+{
+ float float70[70];
+ int int35[35], logical[5];
+ char kvalues[192];
+ const int maxversion = 6;
+
+ need_byte_reversal = false;
+ if( !(sizeof(float)==4) || !(sizeof(int)==4) || !(sizeof(char)==1) )
+ {
+ cout << "readSACheader: ERROR, size of datatypes do not match the SAC specification. Can not read SAC file "
+ << fname << endl;
+ return;
+ }
+
+// Open SAC file
+ FILE* fd=fopen(fname,"r");
+ if( fd == NULL )
+ {
+ cout << "readSACheader: ERROR, observed data file " << fname << " could not be opened" << endl;
+ return;
+ }
+
+// Read header data blocks
+ size_t nr = fread(float70, sizeof(float), 70, fd );
+ if( nr != 70 )
+ {
+ cout << "readSACheader: ERROR, could not read float part of header of " << fname << endl;
+ fclose(fd);
+ return;
+ }
+ nr = fread(int35, sizeof(int), 35, fd );
+ if( nr != 35 )
+ {
+ cout << "readSACheader: ERROR, could not read int part of header of " << fname << endl;
+ fclose(fd);
+ return;
+ }
+ nr = fread(logical, sizeof(int), 5, fd );
+ if( nr != 5 )
+ {
+ cout << "readSACheader: ERROR, could not read bool part of header of " << fname << endl;
+ fclose(fd);
+ return;
+ }
+ nr = fread(kvalues, sizeof(char), 192, fd );
+ if( nr != 192 )
+ {
+ cout << "readSACheader: ERROR, could not read character part of header of " << fname << endl;
+ fclose(fd);
+ return;
+ }
+ int nvhdr = int35[6];
+ if( nvhdr <0 || nvhdr > maxversion )
+ {
+ // Something is wrong, try byte reversal
+ int nvhdrold = nvhdr;
+ Byteswapper bswap;
+ bswap.byte_rev( &nvhdr, 1, "int");
+ if( nvhdr >= 0 && nvhdr <= maxversion )
+ {
+ need_byte_reversal = true;
+ bswap.byte_rev( float70, 70, "float");
+ bswap.byte_rev( int35, 35, "int" );
+ bswap.byte_rev( logical, 5, "int" );
+ }
+ else
+ CHECK_INPUT(false,"readSACheader: ERROR reading sac header of file " << fname << " header version = "
+ << nvhdrold << " byte swapped header version = " << nvhdr );
+ }
+
+// Take out wanted information
+ dt = float70[0];
+ t0 = float70[5];
+ lat = float70[31];
+ lon = float70[32];
+ cmpaz = float70[57];
+ cmpinc = float70[58];
+ utc[0] = int35[0];
+ int jday=int35[1];
+
+ convertjday( jday, utc[0], utc[2], utc[1] );
+ utc[3] = int35[2];
+ utc[4] = int35[3];
+ utc[5] = int35[4];
+ utc[6] = int35[5];
+ npts = int35[9];
+}
+
+//-----------------------------------------------------------------------
+void readSACdata( const char* fname, int npts, double* u, bool need_byte_reversal )
+{
+ if( !(sizeof(float)==4) || !(sizeof(int)==4) || !(sizeof(char)==1) )
+ {
+ cout << "readSACdata: ERROR, size of datatypes do not match the SAC specification. Can not read SAC file "
+ << fname << endl;
+ return;
+ }
+
+// Open SAC file
+ FILE* fd=fopen(fname,"r");
+ if( fd == NULL )
+ {
+ cout << "readSACdata: ERROR, observed data file " << fname << " could not be opened" << endl;
+ return;
+ }
+
+// Skip header
+ int retcode = fseek(fd,158*4,SEEK_SET);
+ if( retcode != 0 )
+ {
+ cout << "readSACdata: ERROR, could not skip header in file " << fname << endl;
+ fclose(fd);
+ return;
+ }
+
+// Read data
+ float* uf = new float[npts];
+ size_t nr = fread( uf, sizeof(float), npts, fd );
+ if( nr != npts )
+ {
+ cout << "readSACdata: ERROR, could not read float array of " << fname << endl;
+ delete[] uf;
+ fclose(fd);
+ return;
+ }
+ if( need_byte_reversal )
+ {
+ Byteswapper bswap;
+ bswap.byte_rev( uf, npts, "float" );
+ }
+
+// Return floats as doubles
+ for( int i=0 ; i < npts ; i++ )
+ u[i] = static_cast<double>(uf[i]);
+ delete[] uf;
+ fclose(fd);
+}
+
+
+
diff --git a/src/sacutils.h b/src/sacutils.h
new file mode 100644
index 0000000..d9fd11c
--- /dev/null
+++ b/src/sacutils.h
@@ -0,0 +1,15 @@
+#ifndef SW4_SACUTILS
+#define SACUTILS
+
+int lastofmonth( int year, int month );
+
+void convertjday( int jday, int year, int& day, int& month );
+
+void readSACheader( const char* fname, double& dt, double& t0,
+ double& lat, double& lon, double& cmpaz,
+ double& cmpinc, int utc[7], int& npts,
+ bool& need_byte_reversal );
+
+void readSACdata( const char* fname, int npts, double* u, bool need_byte_reversal=false );
+
+#endif
More information about the CIG-COMMITS
mailing list