[cig-commits] r16286 - in seismo/3D/SPECFEM3D_GLOBE/trunk: . UTILS/Cluster/lsf UTILS/cut_velocity UTILS/lib
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Sat Feb 20 17:38:57 PST 2010
Author: danielpeter
Date: 2010-02-20 17:38:56 -0800 (Sat, 20 Feb 2010)
New Revision: 16286
Added:
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/Cluster/lsf/remap_lsf_machines.pl
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/Cluster/lsf/run_lsf.bash
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/Cluster/lsf/run_lsf.kernel
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/cut_velocity/Makefile
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/INSTALL
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/cut_velocity/compile_cut
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/cut_velocity/cut_velocity.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/lib/rw_ascfile_c.c
Log:
adds in directory UTILS/cut_velocity a Makefile to use `make` for compilation rather than the provided script compile_cut
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/INSTALL
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/INSTALL 2010-02-20 20:40:52 UTC (rev 16285)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/INSTALL 2010-02-21 01:38:56 UTC (rev 16286)
@@ -59,7 +59,7 @@
please make sure, your installation is working and that
the created files 'constant.h' and 'Makefile' satisfy
- you needs.
+ your needs.
more information is given in the manual provided in USERS_MANUAL.
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/Cluster/lsf/remap_lsf_machines.pl
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/Cluster/lsf/remap_lsf_machines.pl (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/Cluster/lsf/remap_lsf_machines.pl 2010-02-21 01:38:56 UTC (rev 16286)
@@ -0,0 +1,23 @@
+#!/usr/bin/perl -w
+
+# this program remaps the LSF output machine file to
+# the standard one-node-per-column machine file
+
+if (@ARGV != 1) {die("remap_lsf_machines.pl machinefile\n");}
+
+$machine = $ARGV[0];
+
+open(FILE,"$machine") or die("Error opening file $machine\n");
+(@junk) = <FILE>;
+close(FILE);
+
+for($i=0;$i<@junk;$i++) {
+ @node_array = split(" ",$junk[$i]);
+ foreach $node (@node_array) {
+ next if ( $node =~ /^[0-9]/ );
+ push(@nodes, $node);
+ }
+}
+foreach $node (@nodes) {
+ print "$node\n";
+}
Property changes on: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/Cluster/lsf/remap_lsf_machines.pl
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/Cluster/lsf/run_lsf.bash
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/Cluster/lsf/run_lsf.bash (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/Cluster/lsf/run_lsf.bash 2010-02-21 01:38:56 UTC (rev 16286)
@@ -0,0 +1,38 @@
+#!/bin/bash
+
+# this is the launch script to run a regular forward simulation
+# on CITerra at Caltech with the LSF scheduler
+# Qinya Liu, Caltech, May 2007
+
+# use the normal queue unless otherwise directed
+queue="-q normal"
+if [ $# -eq 1 ]; then
+ echo "Setting the queue to $1"
+ queue="-q $1"
+fi
+
+d=`date`
+echo "Starting compilation $d"
+make clean
+make xmeshfem3D
+make xcreate_header_file
+xcreate_header_file
+make xspecfem3D
+d=`date`
+echo "Finished compilation $d"
+
+# compute total number of nodes needed
+NPROC_XI=`grep NPROC_XI DATA/Par_file | cut -d = -f 2`
+NPROC_ETA=`grep NPROC_ETA DATA/Par_file | cut -d = -f 2`
+NCHUNKS=`grep NCHUNKS DATA/Par_file | cut -d = -f 2`
+
+# total number of nodes is the product of the values read
+numnodes=$(( $NCHUNKS * $NPROC_XI * $NPROC_ETA ))
+
+#rm -r -f OUTPUT_FILES/*
+
+echo "Submitting job"
+
+# time below is given in hh:mm
+bsub $queue -n $numnodes -W 48:00 -C 0 < go_mesher_solver_lsf_globe.bash
+
Property changes on: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/Cluster/lsf/run_lsf.bash
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/Cluster/lsf/run_lsf.kernel
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/Cluster/lsf/run_lsf.kernel (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/Cluster/lsf/run_lsf.kernel 2010-02-21 01:38:56 UTC (rev 16286)
@@ -0,0 +1,42 @@
+#!/bin/bash
+
+# this is the launch script to run a kernel simulation
+# on CITerra with the LSF scheduler
+# Qinya Liu, Caltech, May 2007
+
+# use the normal queue unless otherwise directed
+queue="-q normal"
+if [ $# -eq 1 ]; then
+ echo "Setting the queue to $1"
+ queue="-q $1"
+fi
+
+rm -f OUTPUT_FILES/*
+
+d=`date`
+echo "Starting compilation $d"
+# you have to make sure you set this first before compiling the mesher & solver
+change_simulation_type.pl -b
+make clean
+make meshfem3D
+make create_header_file
+xcreate_header_file
+make specfem3D
+d=`date`
+echo "Finished compilation $d"
+
+# compute total number of nodes needed
+NPROC_XI=`grep NPROC_XI DATA/Par_file | cut -d = -f 2`
+NPROC_ETA=`grep NPROC_ETA DATA/Par_file | cut -d = -f 2`
+NCHUNKS=`grep NCHUNKS DATA/Par_file | cut -d = -f 2`
+
+# total number of nodes is the product of the values read
+numnodes=$(( $NCHUNKS * $NPROC_XI * $NPROC_ETA ))
+
+#rm -r -f OUTPUT_FILES/*
+
+echo "Submitting job"
+
+# time below is given in hh:mm
+bsub $queue -n $numnodes -W 48:00 -C 0 < go_mesher_solver_lsf_globe.kernel
+
Property changes on: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/Cluster/lsf/run_lsf.kernel
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/cut_velocity/Makefile
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/cut_velocity/Makefile (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/cut_velocity/Makefile 2010-02-21 01:38:56 UTC (rev 16286)
@@ -0,0 +1,29 @@
+# Makefile
+
+#############################################################
+## modify to match your compiler defaults
+
+# compilers
+F90 = gfortran
+CC = gcc
+
+#############################################################
+
+
+LIB_OBJS = cut_velocity.o rw_ascfile_c.o
+
+# targets
+all: xcut_velocity
+
+xcut_velocity: $(LIB_OBJS)
+ ${F90} -Wall -o xcut_velocity $(LIB_OBJS)
+
+cut_velocity.o: cut_velocity.f90
+ ${F90} -Wall -c cut_velocity.f90
+
+rw_ascfile_c.o: ../lib/rw_ascfile_c.c
+ ${CC} -c -o rw_ascfile_c.o ../lib/rw_ascfile_c.c
+
+
+clean:
+ rm -f *.o xcut_velocity
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/cut_velocity/compile_cut
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/cut_velocity/compile_cut 2010-02-20 20:40:52 UTC (rev 16285)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/cut_velocity/compile_cut 2010-02-21 01:38:56 UTC (rev 16286)
@@ -1,3 +1,16 @@
-gcc -c -o rw_ascfile_c.o rw_ascfile_c.c
-gcc -c -o rw_fortran_wrapper.o rw_fortran_wrapper.c
-ifort -O3 -o xcut_velocity cut_velocity.f90 rw_ascfile_f.f90 rw_ascfile_c.o rw_fortran_wrapper.o
+# use make
+echo "xcut_velocity:"
+echo
+echo "please,in this directory UTILS/cut_velocity"
+echo "add your compiler specifics to: Makefile"
+echo
+echo "and run in this directory UTILS/cut_velocity:"
+echo
+echo "> make"
+echo
+echo
+
+# obsolete...
+#gcc -c -o rw_ascfile_c.o rw_ascfile_c.c
+#gcc -c -o rw_fortran_wrapper.o rw_fortran_wrapper.c
+#ifort -O3 -o xcut_velocity cut_velocity.f90 rw_ascfile_f.f90 rw_ascfile_c.o rw_fortran_wrapper.o
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/cut_velocity/cut_velocity.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/cut_velocity/cut_velocity.f90 2010-02-20 20:40:52 UTC (rev 16285)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/cut_velocity/cut_velocity.f90 2010-02-21 01:38:56 UTC (rev 16286)
@@ -1,18 +1,16 @@
program cut_velocity
-! this program cuts certain portion of the seismograms and convert them into
-! the adjoints sources for generating banana-dougnut kernels.
-! rotation conforms with the SAC 'rotate to gcp normal' command, which says
-! R = cos(th)*N+sin(th)*E
-! T = -sin(th)*N+cos(th)*E
-! where th = baz - pi
+! this program cuts certain portion of the seismograms and converts them into
+! the adjoint sources for generating banana-dougnut kernels.
+! Qinya Liu, Caltech, May 2007
!
-! Qinya Liu, Caltech, May 2007
-
+! call by: ./xcut_velocity t1 t2 ifile[0-5] E/N/Z-ascii-files [baz]
+!
implicit none
- integer :: i, nfile, is, ie, nstep, j, itime ,ifile,ios, i1, i2, nstep_old
- character(len=100) :: arg(100), file(100)
+ integer :: i, is, ie, nstep, j, itime ,ifile,ios, i1, i2, nstep_old
+ character(len=256) :: arg(100), file(100)
+ character(len=256) :: filename
integer,parameter :: NMAX = 30000
real*8, parameter :: EPS = 1.0d-17
real*8, parameter :: PI = 3.1415926d0
@@ -23,9 +21,25 @@
i = 1
lrot = .false.
+ ! reads in file arguments
do while (1 == 1)
call getarg(i,arg(i))
- if (i < 6 .and. trim(arg(i)) == '') stop 'cut_velocity t1 t2 ifile[0-5] E/N/Z-ascii-files [baz]'
+ if (i < 6 .and. trim(arg(i)) == '') then
+ print*,'Usage: '
+ print*,' xcut_velocity t1 t2 ifile[0-5] E/N/Z-ascii-files [baz]'
+ print*,'with'
+ print*,' t1: window start time'
+ print*,' t2: window end time'
+ print*,' ifile: 0 = adjoint source calculated for each seismogram component'
+ print*,' ifile: 1 = adjoint source given by East component only'
+ print*,' ifile: 2 = adjoint source given by North component'
+ print*,' ifile: 3 = adjoint source given by Z component'
+ print*,' ifile: 4 = adjoint source given by rotated transversal component (requires baz)'
+ print*,' ifile: 5 = adjoint source given by rotated radial component (requires baz)'
+ print*,' E/N/Z-ascii-files : displacement traces stored as ascii files'
+ print*,' [baz]: (optional) back-azimuth, requires ifile = 4 or ifile = 5'
+ stop 'cut_velocity t1 t2 ifile[0-5] E/N/Z-ascii-files [baz]'
+ endif
if (trim(arg(i)) == '') exit
if (i == 1) then
read(arg(i),*,iostat=ios) ts
@@ -48,6 +62,7 @@
i = i + 1
enddo
+ ! checks rotation baz and ifile parameter
i = i - 1
if (lrot) then
if (i /= 7) stop 'cut_velocity t1 t2 ifile[0-5] E/N/Z-ascii-files [baz]'
@@ -60,52 +75,73 @@
if (i /= 6) stop 'cut_velocity t1 t2 ifile[0-5] E/N/Z-ascii-files [baz]'
endif
+ ! user output
print *, 'ifile = ', ifile, ' lrot = ', lrot
print *, ' '
-
+ ! reads seismograms (ascii format)
do i = 1, 3
- print *, 'reading asc file '//trim(file(i))//' ...'
- call dread_ascfile_f(file(i),t0,dt,nstep,data(i,:))
+ filename = trim(file(i))
+ print *, 'reading asc file '//trim(filename)//' ...'
+ call dread_ascfile_c(trim(filename)//char(0),t0,dt,nstep,data(i,:))
if (nstep > NMAX) stop 'Change the data array range limit'
if (i == 1) then
t0_old = t0; dt_old = dt; nstep_old = nstep
else
- if (i > 1 .and. abs(t0_old - t0) > EPS .and. abs(dt_old - dt) > EPS .and. nstep_old /= nstep) &
+ if (i > 1 .and. abs(t0_old - t0) > EPS &
+ .and. abs(dt_old - dt) > EPS &
+ .and. nstep_old /= nstep) &
stop 'Error different t0, dt, nstep'
endif
enddo
+ print *, ' '
+ print *, 'start time:',t0
+ print *, 'time step:',dt
+ print *, 'number of steps:',nstep
+ print *, ' '
+ ! component rotation
if (lrot) then
data(4,:) = costh * data(1,:) - sinth * data(2,:)
data(5,:) = sinth * data(1,:) + costh * data(2,:)
- call dwrite_ascfile_f('t.txt',t0,dt,nstep,data(4,:))
- call dwrite_ascfile_f('r.txt',t0,dt,nstep,data(5,:))
+ call dwrite_ascfile_c('t.txt',t0,dt,nstep,data(4,:))
+ call dwrite_ascfile_c('r.txt',t0,dt,nstep,data(5,:))
i1 = 3; i2 = 5
else
i1 = 1; i2 = 3
endif
-
+ ! loops over seismogram components
do i = i1, i2
+ ! start and end index
is = (ts - t0) / dt + 1
ie = (te - t0) / dt + 1
if (is < 1 .or. ie <= is .or. ie > nstep) then
print *, 'Error in ts, te'; stop
endif
+
+ ! time window (parabola shaped)
tw(1:nstep) = 0.
+ if( i == i1 ) open(44,file='plot_time_window.txt',status='unknown')
do j = is, ie
tw(j) = 1 - (2 * (dble(j) - is)/(ie - is) - 1) ** 2
+ if( i == i1 ) write(44,*) j,tw(j)
enddo
+ if( i == i1 ) close(44)
+
+ ! calculates velocity (by finite-differences)
do itime = 2, nstep-1
out(itime) = (data(i,itime+1) - data(i,itime-1)) / (2 * dt)
enddo
out(1) = (data(i,2) - data(i,1)) / dt
out(nstep) = (data(i,nstep) - data(i,nstep-1)) /dt
+
+ ! normalization factor
norm = dt * sum( tw(1:nstep) * out(1:nstep) * out(1:nstep))
print *, 'i = ', i, 'norm = ', norm
if (ifile /= 0 .and. ifile /= i) norm = 0.0
-
+
+ ! adjoint source
if (abs(norm) > EPS) then
adj(1:nstep) = - out(1:nstep) * tw(1:nstep) / norm
else
@@ -113,18 +149,25 @@
adj(:) = 0.
endif
data(i,:) = adj(:)
+
enddo
+ print *, ' '
+ ! component rotation back to cartesian x-y-z
if (lrot) then
- call dwrite_ascfile_f('t-cut.txt',t0,dt,nstep,data(4,:))
- call dwrite_ascfile_f('r-cut.txt',t0,dt,nstep,data(5,:))
+ call dwrite_ascfile_c('t-cut.txt',t0,dt,nstep,data(4,:))
+ call dwrite_ascfile_c('r-cut.txt',t0,dt,nstep,data(5,:))
data(1,:) = costh * data(4,:) + sinth * data(5,:)
data(2,:) = -sinth * data(4,:) + costh * data(5,:)
endif
+ ! file output for component BHE/BHN/BHZ
do i = 1, 3
- print *, 'write to asc file '//trim(file(i))//'.ad'
- call dwrite_ascfile_f(trim(file(i))//'.ad',t0,dt,nstep,data(i,:))
+ filename = trim(file(i))//'.adj'
+ print *, 'write to asc file '//trim(filename)
+ call dwrite_ascfile_c(trim(filename)//char(0),t0,dt,nstep,data(i,:))
enddo
end program cut_velocity
+
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/lib/rw_ascfile_c.c
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/lib/rw_ascfile_c.c 2010-02-20 20:40:52 UTC (rev 16285)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/lib/rw_ascfile_c.c 2010-02-21 01:38:56 UTC (rev 16286)
@@ -1,7 +1,10 @@
#include <stdio.h>
-#include "drw_ascfile.h"
+#include <stdlib.h>
+//#include "drw_ascfile.h"
+// Qinya Liu, Caltech, May 2007
-void dread_ascfile(const char *ascfile,
+
+void dread_ascfile_c_(const char *ascfile,
double *t0, double *dt, int *n,
double *data)
@@ -18,33 +21,41 @@
while ( fscanf(fd,"%lf %lf\n",&junk, data+i) != EOF ) {
if (i == 0) junk1 = junk;
if (i == 1) *dt = junk - junk1;
- i++;}
+ i++;
+ }
*t0 = junk1;
*n = i;
if (fclose(fd) != 0) {
printf(" file %s cannot be closed\n",ascfile);
- exit(1);}
+ exit(1);
+ }
}
-void dwrite_ascfile(const char *ascfile,
- double t0, double dt, int n,
+void dwrite_ascfile_c_(const char *ascfile,
+ double *t0, double *dt, int *n,
const double *data)
{
FILE *fd;
int i;
+ double tmp,tmp0;
+ //printf("t0: %f ,dt: %f ,n: %i \n",*t0,*dt,*n);
+
if ((fd = fopen(ascfile,"w")) == NULL) {
printf(" file %s cannot be opened to write\n",ascfile);
exit(1);
}
i = 0;
- for (i=0; i<n; i++) {
- fprintf(fd,"%14.7g %18.7g\n", t0+i*dt, data[i]);
+ tmp = *dt;
+ tmp0 = *t0;
+ for (i=0; i< *n; i++) {
+ fprintf(fd,"%14.7g %18.7g\n", tmp0+i*tmp, data[i]);
}
if (fclose(fd) != 0) {
printf("file %s cannot be closed\n",ascfile);
- exit(1);}
+ exit(1);
+ }
}
More information about the CIG-COMMITS
mailing list