[cig-commits] r16135 - in seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate: PLOTTING src
carltape at geodynamics.org
carltape at geodynamics.org
Fri Jan 15 14:07:22 PST 2010
Author: carltape
Date: 2010-01-15 14:07:21 -0800 (Fri, 15 Jan 2010)
New Revision: 16135
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model_all.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex08.f90
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_solver.f90
Log:
Included a special example for reading in an initial heterogeneous model and a heterogenous target model.
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model.pl 2010-01-15 21:10:53 UTC (rev 16134)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model.pl 2010-01-15 22:07:21 UTC (rev 16135)
@@ -4,7 +4,7 @@
#
# plot_surf_model.pl
# Carl Tape
-# 18-Nov-2009
+# 15-Jan-2010
#
# Plot the source time function, the velocity model for the synthetics,
# and the velocity model for the data (optional).
@@ -23,6 +23,7 @@
# ../plot_surf_model.pl 4 500 0 1/0/5 0.1/0.05 stf_00001_1 0
# ../plot_surf_model.pl 4 600 0 1/0/5 0.1/0.05 stf_00001_1 0
# ../plot_surf_model.pl 4 700 0 1/0/5 0.1/0.05 stf_00001_1 0
+# ../plot_surf_model.pl 4 800 0 0/1/5 0.1/0.05 stf_00001_1 0
#
#==========================================================
@@ -131,7 +132,7 @@
$Dscale1 = "-D$Dx/$Dy/$Dlen/0.10h";
#$Bscale1 = sprintf("-B%2.2e:\" Phase Velocity ( km s\@+-1\@+ )\": -E10p",$ptick);
-$Bscale1 = sprintf("-B%2.2e:\" Percent Perturbation from %3.3f km s\@+-1\@+ \": -E10p",$ptick,$beta0/1000);
+$Bscale1 = sprintf("-B%2.2e:\" Perturbation from %3.3f km s\@+-1\@+ \": -E10p",$ptick,$beta0/1000);
#-------------------------
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model_all.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model_all.pl (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model_all.pl 2010-01-15 22:07:21 UTC (rev 16135)
@@ -0,0 +1,266 @@
+#!/usr/bin/perl -w
+
+#==========================================================
+#
+# plot_surf_model_all.pl
+# Carl Tape
+# 15-Jan-2010
+#
+# Plot a 6-panel figure for each model iteration.
+# Input allows for a range of models as well.
+#
+# Example:
+# ../plot_surf_model_all.pl 800 0/32 0/1/5 0.1/0.05 0
+#
+#==========================================================
+
+if (@ARGV < 5) {die("Usage: plot_surf_model_all.pl xxx\n");}
+($irun0,$imodels,$ievents,$pbar,$ifinite) = @ARGV;
+$comp = 1;
+
+$icolor = 1;
+$ijpg = 0;
+$ipdf = 1;
+
+# base directory
+$pwd = $ENV{PWD};
+$dirplot = `dirname $pwd`; chomp($dirplot);
+$basedir = `dirname $dirplot`; chomp($basedir);
+#$basedir = "/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_work";
+
+($imodelmin,$imodelmax) = split("/",$imodels);
+($pmax,$ptick) = split("/",$pbar);
+($ievent_all,$ievent_one,$ievent) = split("/",$ievents);
+$rec_file1 = "sr.txt";
+$efile_syn = "events_syn_lonlat.dat";
+$mfile_syn = "structure_syn.dat";
+$mfile_dat = "structure_dat.dat";
+
+# directories
+$plotdir = "${basedir}/PLOTTING";
+$odir = "${basedir}/OUTPUT";
+$idir = "${basedir}/INPUT";
+$figdir = "${plotdir}/FIGURES";
+if (not -e $figdir) {die("Check if figdir $figdir exist or not\n");}
+
+$stirun0 = sprintf("%4.4i",$irun0);
+
+$dir0 = "$odir/run_${stirun0}";
+$dir2 = sprintf("$dir0/event_%3.3i",$ievent);
+$rec_file = "$dir2/${rec_file1}";
+$evefile = "$dir0/${efile_syn}"; # event file from the base directory
+if (not -f $evefile) { die("Check if $evefile exist or not\n") }
+
+# read in reference values: alpha0, beta0, per
+if (1==1) {
+ $fileref1 = "$dir0/reference_values.dat";
+ $fileref2 = "$dir0/reference_period.dat";
+ open(IN,$fileref1); @lines = <IN>; ($alpha0,$beta0) = split(" ",$lines[0]);
+ open(IN,$fileref2); $per = <IN>; chomp($per);
+} else {
+ $per = 20;
+ $beta0 = 3500;
+}
+$per_lab = sprintf("T = %3.3f s",$per);
+$beta0_lab = sprintf("beta0 = %3.3f km/s",$beta0/1000);
+#print "\n-- $per -- $per_lab -- $beta0_lab --\n"; die("testing");
+
+$plabel = "${plotdir}/plot_surf_model_all.pl";
+
+#die("\ntesting\n");
+
+# plotting specifications
+$fsize0 = "14";
+$fsize1 = "11";
+$fsize2 = "9";
+$fsize3 = "5";
+$fontno = "1"; # 1 or 4
+$tick = "0.1c";
+
+$tinc = 50; # tick increment (in seconds) for source time function
+
+# plot symbols for sources, receivers, and shelf
+if($ifinite==0) {$src0 = "-W0.75p -Sa0.20 -G255/0/0"} else {$src0 = "-Sc0.05"}
+$src = "-W0.75p -Sa0.20";
+$rec = "-W0.5p/0/0/0 -St0.10";
+$rec0 = "-Sc10p -W0.5p";
+$Wshelf = "-W1.0/0/0/0tap";
+
+# resolution of color plots
+$interp = "-I0.5m/0.5m -S4m"; # key information
+$interp = "-I2m/2m -S4m"; # key information
+$grdfile = "temp.grd";
+
+#-------------------------
+# color
+
+# KEY: scaling for color
+$scale_color = 21.0;
+$colorbar = "seis";
+
+#-------------------------
+
+# write plotting scripts
+$wid = 2.5;
+$J1 = "-JM${wid}i"; # in lat-lon
+$origin = "-X1.5 -Y7.25";
+
+$B1dat = "-B1:.\" \":WeSN";
+$B1syn = "-B1:.\" \":wESN";
+
+$Dlen = 2.0;
+$Dx = $wid/2;
+$Dy = -0.35;
+
+$Dscale1 = "-D$Dx/$Dy/$Dlen/0.15h";
+
+$Bscale1 = sprintf("-B%2.2e:\" Perturbation from %3.3f km s\@+-1\@+ \": -E10p",$ptick,$beta0/1000);
+$Bscale2 = sprintf("-B%2.2e:\" Perturbation from m00 \": -E10p",$ptick);
+
+#===============================================
+$cshfile = "plot_surf_model_all.csh";
+print "\nWriting to CSH file ${cshfile}...\n";
+open(CSH,">$cshfile");
+print CSH "gmtset BASEMAP_TYPE plain PAPER_MEDIA letter TICK_LENGTH $tick LABEL_FONT_SIZE $fsize2 ANOT_FONT_SIZE $fsize2 PLOT_DEGREE_FORMAT D HEADER_FONT $fontno ANOT_FONT $fontno LABEL_FONT $fontno HEADER_FONT_SIZE $fsize1\n";
+#===============================================
+
+# make colorpoint file
+# $T1 = "-T3/4/0.1";
+# $pmax = 10;
+$dc = $pmax/10;
+$T1 = "-T-$pmax/$pmax/$dc";
+$cptfile1 = "color0.cpt";
+print CSH "makecpt -C$colorbar -D $T1 > $cptfile1\n";
+
+#$rec_file = "$edir/sr.txt";
+if (not -f $rec_file) {die("Check if rec_file $rec_file exist or not\n");}
+
+# sides to show the lat-lon tick marks
+$B0 = "-B1:.\" \":";
+ at Bopts = ("WesN","wEsN","Wesn","wEsn","Wesn","wEsn");
+
+$dX = -1.1*$wid; $dY = -1.1*$wid; $shift1 = "-X$dX -Y$dY";
+$dX = 1.1*$wid; $dY = 0; $shift2 = "-X$dX -Y$dY";
+ at shifts = ($origin,$shift2,$shift1,$shift2,$shift1,$shift2);
+
+ at finds = (7,7,7,7,7,7,7);
+
+$xlab = $wid/2;
+$ylab = -0.02*$wid;
+$olab = "-Xa$xlab -Ya$ylab";
+$J_title = "-JM1";
+$R_title = "-R0/1/0/1";
+
+#================
+
+# loop over models
+for ($m = $imodelmin; $m <= $imodelmax; $m = $m+1) {
+
+ $imodel = $m;
+ $stirun = sprintf("%4.4i",$irun0+$imodel);
+ $stmod = sprintf("m%2.2i",$imodel/2);
+
+ $dir1 = "$odir/run_${stirun}";
+
+ $name = "model_all_${stirun}";
+ $psfile = "$name.ps";
+ $jpgfile = "$name.jpg";
+
+ @files = ("$dir0/${mfile_syn}","$dir0/${mfile_syn}","$dir1/${mfile_syn}","$dir1/${mfile_syn}","$dir0/${mfile_dat}","$dir0/${mfile_dat}");
+ @stlabs = ("m00","ln(m00/m00)",$stmod,"ln($stmod/m00)","mtar","ln(mtar/m00)");
+
+ # only plot even models (m04 = iteration 8)
+ $kmax = 6*(1 - $m % 2);
+
+ # loop over subplots -- only for even m
+ for ($k = 0; $k < $kmax; $k = $k+1) {
+
+ $iB = $Bopts[$k];
+ $shift = $shifts[$k];
+ $B = "${B0}${iB}";
+ $find = $finds[$k];
+ $find2 = 7 + $find;
+
+ $file0 = $files[0];
+ $file1 = $files[$k];
+ if (not -f $file1) {die("Check if $file1 exist or not\n");}
+
+ if ($k % 2 == 1) {
+ $file2 = "ftemp";
+ print CSH "paste $file1 $file0 > $file2\n";
+ }
+
+ # phase velocity map
+ if ($k==0) {
+ # set bounds for the plotting
+ #$xmin = -121; $xmax = -115.0; $zmin = 31.75; $zmax = 36.5;
+ ($xmin,$xmax,$zmin,$zmax,$smin,$smax,$tmin,$tmax) = split(" ",`minmax -C $file1`);
+ $dinc = 0.25; # buffer around min/max
+ $xmin = $xmin-$dinc; $xmax = $xmax+$dinc;
+ $zmin = $zmin-$dinc; $zmax = $zmax+$dinc;
+ $R = "-R$xmin/$xmax/$zmin/$zmax";
+ print "\n$R\n";
+
+ print CSH "psbasemap $B $R $J1 -K -V -P $origin > $psfile\n"; # START
+ } else {
+ print CSH "psbasemap $B $R $J1 -K -O -V $shift >> $psfile\n";
+ }
+
+ if ($icolor == 1) {
+ #print CSH "awk '{print \$1,\$2,\$7}' $file1 | pscontour $R $J1 -A- -C$cptfile1 -I -K -O -V >> $psfile\n";
+ if ($k % 2 == 0) {
+ print CSH "awk '{print \$1,\$2,\$${find}}' $file1 | nearneighbor -G$grdfile $R $interp\n";
+ } else {
+ print CSH "awk '{print \$1,\$2,\$${find}-\$${find2}}' $file2 | nearneighbor -G$grdfile $R $interp\n";
+ }
+ print CSH "grdimage $grdfile -C$cptfile1 $J1 -K -O -V -Q >> $psfile\n";
+ }
+ print CSH "pscoast $J1 $R $B -W1p -Na/1p -Dh -K -O -V >> $psfile\n";
+
+ # plot receivers
+ if (0==1) { # numbered large circles
+ print CSH "awk '\$1 == \"R\" {print \$2,\$3}' $rec_file | psxy -N $J1 $R -K -O -V $rec0 >> $psfile\n";
+ $rec_file2 = text_rec; $angle = 0; $just = "CM";
+ print CSH "awk '\$1 == \"R\" {print \$2,\$3,$fsize3,$angle,$fontno,\"$just\",\$4}' $rec_file > $rec_file2\n";
+ print CSH "pstext $rec_file2 -N $J1 $R -K -O -V >> $psfile\n";
+
+ } else { # small circles
+ print CSH "awk '\$1 == \"R\" {print \$2,\$3}' $rec_file | psxy -N -Sc5p -W0.5p $J1 $R -K -O -V >> $psfile\n";
+
+ }
+
+ # plot sources
+ if ($ievent_one) {print CSH "awk '\$1 == \"S\" {print \$2,\$3}' $rec_file | psxy $src0 -N $J1 $R -K -O -V >> $psfile\n";}
+ if ($ievent_all) {print CSH "awk '{print \$1,\$2}' $evefile | psxy $src -N $J1 $R -K -O -V >> $psfile\n";}
+
+ # plot label
+ $stlab = $stlabs[$k];
+ $textinfo = "-N -C4p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
+ print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab >>$psfile<<EOF\n0 0 12 0 $fontno CM $stlab\nEOF\n";
+
+ # colorscale bars
+ if ($k==4) {print CSH "psscale -C$cptfile1 $Dscale1 $Bscale1 -K -O -V >> $psfile \n";}
+ if ($k==5) {print CSH "psscale -C$cptfile1 $Dscale1 $Bscale2 -K -O -V >> $psfile \n";}
+
+ # plot title and GMT header
+ if ($k==5) {
+ $plot_title = "Model m00, model $stmod, and target model (run_${stirun})";
+ $Utag = "-U/0/0.25/$plabel";
+ $shift = "-Xa-$dX -Ya8.5";
+ print CSH "pstext -N $J1 $R $Utag -O -V $shift >>$psfile<<EOF\n $xmin $zmin $fsize0 0 $fontno LM $plot_title\nEOF\n"; # FINISH
+
+ if($ijpg == 1) {print CSH "convert $psfile $jpgfile\n";}
+ if($ipdf == 1) {print CSH "ps2pdf $psfile\n";}
+
+ }
+ } # for loop
+
+ #-----------------------------
+
+} # loop over models
+
+close (CSH);
+system("csh -f $cshfile");
+system("gv $psfile &")
+
+#=================================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model_all.pl
___________________________________________________________________
Name: svn:executable
+ *
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d.f90 2010-01-15 21:10:53 UTC (rev 16134)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d.f90 2010-01-15 22:07:21 UTC (rev 16135)
@@ -42,7 +42,7 @@
double precision :: d,dmin
double precision, dimension(:), allocatable :: shelf_lat,shelf_lon,shelf_z,shelf_x
double precision, dimension(:), allocatable :: coast_lat,coast_lon,coast_z,coast_x
- character(len=200) :: shelf_file,coast_file
+ character(len=200) :: shelf_file,coast_file,idir
! corners of the grid (for membrane surface wave simulations)
double precision, dimension(4) :: corners_utm_x, corners_utm_z, corners_lon, corners_lat
@@ -141,7 +141,8 @@
double precision, dimension(NGLLX,NGLLZ,NSPEC) :: kbulk_smooth, kbeta_smooth ! scale_array
! double precision, dimension(NGLLX,NGLLZ,NSPEC) :: bulk_pert_syn, beta_pert_syn
double precision :: scale_bulk, scale_beta, scale_source, scale_struct, scale_struct_gji
- double precision :: sigma_bulk, sigma_beta, sigma_xs, sigma_zs, sigma_ts, sigma_checker_scale
+ double precision :: sigma_bulk0, sigma_beta0, sigma_bulk, sigma_beta, sigma_checker_scale
+ double precision :: sigma_xs, sigma_zs, sigma_ts
double precision :: cov_fac, cov_src_fac, cov_str_fac
double precision :: afac_target, coverage_str, coverage_src
double precision :: joint_str_src_grad_ratio, jsw_str, jsw_src, joint_total_weight, jfac
@@ -333,7 +334,7 @@
endif
!----------
- ! KEY: define 1D model for body waves
+ ! define 1D model for body waves
! SoCal 1D model (Dreger & Helmberger 1990)
!z_breaks(1) = 5500.0 ; z_breaks(2) = 16000.0 ; z_breaks(3) = 35000.0 ; nlayer = 4
@@ -705,9 +706,65 @@
w_scale = 2*pi/(Nfac*2*beta0*t_target)
endif
+ !----------
+ ! load in heterogeneous reference and target models
+ ! NOTE: THESE ARE ONLY FOR A SPECIFIC MESH
+ if ((IMODEL_SYN == 3).and.(IMODEL_DAT == 3)) then
+
+ idir = '/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_INPUT/random_fields/model_001/'
+
+ ! read in structure model for synthetics
+ write(filename2,'(a)') trim(idir)//'structure_syn_in.dat'
+ open(unit=19,file=filename2,status='old')
+ do ispec = 1,NSPEC
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ read(19,'(7e20.10)') temp1, temp2, temp3, temp4, temp5, temp6, temp7
+ kappa_syn(i,j,ispec) = temp3
+ mu_syn(i,j,ispec) = temp4
+ rho_syn(i,j,ispec) = temp5
+ enddo
+ enddo
+ enddo
+ close(19)
+
+ ! read in structure model for data
+ write(filename2,'(a)') trim(idir)//'structure_dat_in.dat'
+ open(unit=20,file=filename2,status='old')
+ do ispec = 1,NSPEC
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ read(20,'(7e20.10)') temp1, temp2, temp3, temp4, temp5, temp6, temp7
+ kappa_dat(i,j,ispec) = temp3
+ mu_dat(i,j,ispec) = temp4
+ rho_dat(i,j,ispec) = temp5
+ enddo
+ enddo
+ enddo
+ close(20)
+
+ beta_syn = sqrt( mu_syn / rho_syn )
+ beta_dat = sqrt( mu_dat / rho_dat )
+
+ ! read reference values to file
+ idir = '/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_INPUT/meshes/mesh_001/'
+ open(unit=16,file=trim(idir)//'reference_values.dat', status='unknown')
+ read(16,'(6e20.10)') alpha0, beta0, rho0, bulk0, kappa0, mu0
+ close(16)
+
+ ! sigma values describing the distribution
+ sigma_beta0 = 0.05
+ sigma_bulk0 = 0.05
+
+ endif
+
!=========================================
! (INITIAL) REFERENCE MODEL (synthetics)
+ ! initialize
+ !kappa_syn = 0.0 ; mu_syn = 0.0 ; rho_syn = 0.0
+ !alpha_syn = 0.0 ; beta_syn = 0.0 ; bulk_syn = 0.0
+
! set velocity field for the synthetics (e.g., homogeneous reference model)
! assigns local arrays mu_syn, kappa_syn, rho_syn
iref = 1
@@ -725,15 +782,17 @@
endif
! compute alpha0 and beta0 by finding the average of the phase velocity map
- print *, 'computing the average of the phase velocity map'
- kappa0 = sum(kappa_syn * da_local) / (LENGTH * HEIGHT)
- mu0 = sum( mu_syn * da_local) / (LENGTH * HEIGHT)
- rho0 = sum( rho_syn * da_local) / (LENGTH * HEIGHT)
- alpha0 = sqrt( (kappa0 + FOUR_THIRDS*mu0) / rho0 )
- beta0 = sqrt( mu0 / rho0 )
- bulk0 = sqrt( kappa0 / rho0 )
- !alpha0 = sum(alpha_initial_syn * da_local) / (LENGTH * HEIGHT)
- !beta0 = sum(beta_initial_syn * da_local) / (LENGTH * HEIGHT)
+ if ( (IMODEL_SYN /= 3) .or. (IMODEL_DAT /= 3) ) then
+ print *, 'computing the average of the phase velocity map'
+ kappa0 = sum(kappa_syn * da_local) / (LENGTH * HEIGHT)
+ mu0 = sum( mu_syn * da_local) / (LENGTH * HEIGHT)
+ rho0 = sum( rho_syn * da_local) / (LENGTH * HEIGHT)
+ alpha0 = sqrt( (kappa0 + FOUR_THIRDS*mu0) / rho0 )
+ beta0 = sqrt( mu0 / rho0 )
+ bulk0 = sqrt( kappa0 / rho0 )
+ !alpha0 = sum(alpha_initial_syn * da_local) / (LENGTH * HEIGHT)
+ !beta0 = sum(beta_initial_syn * da_local) / (LENGTH * HEIGHT)
+ endif
!!$ ! initial and reference velocity structure
!!$ ! From this model, we compute fractional perturbations.
@@ -743,6 +802,10 @@
!=========================================
! TARGET MODEL (data)
+ ! initialize
+ !kappa_dat = 0.0 ; mu_dat = 0.0 ; rho_dat = 0.0
+ !alpha_dat = 0.0 ; beta_dat = 0.0 ; bulk_dat = 0.0
+
! set velocity field for the data (perturbed reference model)
! assigns local arrays mu_dat, kappa_dat, rho_dat
iref = 0
@@ -753,14 +816,10 @@
! to what was used in the base directory for the CG algorithm.
if (HESSIAN .and. INV_STRUCT == 1) then
- kappa_syn = 0.0
- mu_syn = 0.0
- rho_syn = 0.0
+ kappa_syn = 0.0 ; mu_syn = 0.0 ; rho_syn = 0.0
+ alpha_syn = 0.0 ; beta_syn = 0.0 ; bulk_syn = 0.0
- ! THERE IS SOME UNRESOLVED ISSUE WITH EITHER THE MATLAB FILES WRITING OVER OLD FILES,
- ! OR THE SCRATCH DISK ON DENALI. WHATEVER THE CASE, WE MUST CHECK THAT THE OUTPUT
- ! FILE IS THE SAME AS THE ONE THAT IS READ IN.
-
+ ! read in structure model for synthetics
write(filename2,'(a,i2.2,a)') trim(out_dir2)//'structure_syn_h',hmod,'.dat'
open(unit=19,file=filename2,status='unknown')
do ispec = 1,NSPEC
@@ -798,7 +857,7 @@
bulk_syn = sqrt( kappa_syn / rho_syn )
endif
- ! unperturbed structure (set target map to the map used to compute synthetics)
+ ! unperturbed structure (set target map to be the map used to compute synthetics)
if(PERT_STRUCT == 0) then
kappa_dat = kappa_syn
mu_dat = mu_syn
@@ -1645,8 +1704,13 @@
! scaling vector for STRUCTURE parameters
sigma_checker_scale = 2.0 ! to mimic a Gaussian distibution
- m_scale_str(1) = (afac/sigma_checker_scale)/100.0 ! dimensionless : bulk* = ln(bulk/bulk0)
- m_scale_str(2) = (afac/sigma_checker_scale)/100.0 ! dimensionless : beta* = ln(beta/beta0)
+ if ((IMODEL_SYN == 3).and.(IMODEL_DAT == 3)) then
+ m_scale_str(1) = sigma_bulk0
+ m_scale_str(2) = sigma_beta0
+ else
+ m_scale_str(1) = (afac/sigma_checker_scale)/100.0 ! dimensionless : bulk* = ln(bulk/bulk0)
+ m_scale_str(2) = (afac/sigma_checker_scale)/100.0 ! dimensionless : beta* = ln(beta/beta0)
+ endif
! scaling vector for SOURCE parameters
if( GJI_PAPER == 1 ) then
@@ -1753,6 +1817,9 @@
! Thus, the norm of the target model will then be somewhat GREATER than 1.0.
coverage_str = 0.666 / 0.962
coverage_src = 0.946 / 1.018
+
+ ! If the initial and target models are from a Gaussian distribution,
+ ! then this factor is not needed.
!coverage_str = 1.0
!coverage_src = 1.0
@@ -2547,7 +2614,7 @@
! write the current models to file
! only write the files for the first event, since they are the same for each event
- if(ifirst_event == 1) then
+ if (ifirst_event == 1) then
! write velocity structure for data and synthetics to file -- LOCAL LEVEL
file_structure_dat = 'structure_dat.dat'
@@ -3315,7 +3382,7 @@
mt_vec(1 : nmod_str) = m0_vec_initial(1 : nmod_str) ! initial structure
else
mt_vec(1 : NLOCAL) = bulk0 * exp( mt(1 : NLOCAL) )
- mt_vec(NLOCAL+1 : nmod_str) = beta0 * exp( mt(NLOCAL+1 : nmod_str) )
+ mt_vec(NLOCAL+1 : nmod_str) = beta0 * exp( mt(NLOCAL+1 : nmod_str) )
endif
! SOURCE PARAMETERS
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex08.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex08.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex08.f90 2010-01-15 22:07:21 UTC (rev 16135)
@@ -0,0 +1,336 @@
+module wave2d_constants
+
+! This file is copied into the output directory when wave2d.f90 is run
+! There are two basic options: membrane (surface) waves and body waves.
+! Several lines need to be commented/uncommented to go from one to the other.
+
+! index of the reference directory for the simulation output
+ integer, parameter :: IRUNZ = 800
+
+!========================
+! GRID, TIME-STEP, AND SOURCE PARAMETERS
+
+! NFRAME : number of frames to save: membrane, 10; body, 20
+! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
+! NSTEP : number of timesteps
+ integer, parameter :: NFRAME = 10 ! membrane
+ !integer, parameter :: NFRAME = 11 ! body
+ integer, parameter :: NSAVE = 400 ! 200,400
+ integer, parameter :: NSTEP = NFRAME*NSAVE
+
+! time step in seconds
+ !double precision, parameter :: DT = 2.0e-2 ! body waves
+ double precision, parameter :: DT = 6.0e-2 ! membrane surface waves
+
+! temporal properties of source (source time function)
+ integer, parameter :: ISRC_TIME = 1 ! type (1)
+ double precision, parameter :: hdur = 10.0 ! HALF-duration (s), membrane
+ !double precision, parameter :: hdur = 2.0 ! HALF-duration (s), body
+ double precision, parameter :: tshift = 2.0*DT*dble(NSAVE) ! time shift (s)
+ !double precision, parameter :: tshift = 8.0*hdur
+ logical, parameter :: SRC_TAPER = .true. ! taper the endpoints of the time series
+
+! normalization factor of point source force
+ double precision, parameter :: FNORM = 1.0d10
+
+! forward wavefield source-time function
+! (For membrane waves, FOR_Y = 1 is all that matters.)
+ integer, parameter :: FOR_X = 1 ! boolean
+ integer, parameter :: FOR_Y = 1 ! boolean
+ integer, parameter :: FOR_Z = 1 ! boolean
+
+! adjoint wavefield source-time function
+! (For membrane waves, REV_Y = 1 is all that matters.)
+! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
+ integer, parameter :: REV_X = 1 ! boolean
+ integer, parameter :: REV_Y = 1 ! boolean
+ integer, parameter :: REV_Z = 1 ! boolean
+
+! spatial properties of sources
+! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
+! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
+! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
+ integer, parameter :: ISRC_SPACE = 6 ! see wave2d.f90
+
+! spatial properties of receivers
+! IREC_SPACE
+! (1) individual station(s)
+! (2) SoCal (used for GJI paper)
+! (3) regular mesh on land
+! (4) regular mesh
+! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
+! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
+! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
+ integer, parameter :: IREC_SPACE = 2 ! see wave2d.f90
+ integer, parameter :: NMESH_REC = 17
+ double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m ! 4km for membrane surface waves
+ double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m ! 15km for membrane surface waves
+ double precision, parameter :: STATION_COAST_BUFFER = 0.0e3 ! m
+
+! lower right corner for membrane surface waves plotting grid
+ double precision, parameter :: LAT_MIN = 32.0
+ double precision, parameter :: LON_MIN = -120.0
+ integer, parameter :: UTM_PROJECTION_ZONE = 11 ! southern California
+
+! mesh specifications: membrane surface waves
+ double precision, parameter :: LENGTH = 480.0e3 ! m
+ double precision, parameter :: HEIGHT = 480.0e3 ! m
+ double precision, parameter :: AREA = LENGTH*HEIGHT
+ integer, parameter :: NEX = 40
+ integer, parameter :: NEZ = 40
+!!$
+!!$! mesh specifications: body waves
+!!$ double precision, parameter :: LENGTH = 200.0e3 ! m ! 400 for 1D body waves
+!!$ double precision, parameter :: HEIGHT = 80.0e3 ! m
+!!$ integer, parameter :: NEX = 80 ! 160
+!!$ integer, parameter :: NEZ = 32 ! 32
+
+!========================
+! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
+
+! model perturbations for HOMOGENEOUS model (or perturbation)
+! scaling from beta to alpha
+! value is from Master et al. (2000), "The relative behavior of shear velocity..."
+ double precision, parameter :: R_BETA_OVER_ALPHA = 1.3
+ double precision, parameter :: PBETA = 10.0
+ double precision, parameter :: PALPHA = PBETA / R_BETA_OVER_ALPHA
+ double precision, parameter :: PRHO = 0.0
+
+ ! reference model and target model choice
+ integer, parameter :: IMODEL_SYN = 3
+ integer, parameter :: IMODEL_DAT = 3
+ !-----------------------------------------------------------------------------------------
+ ! 0 1 2 3
+ !-----------------------------------------------------------------------------------------
+ ! ISURFACE=1, IMODEL_SYN : homo checker het
+ ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het pert
+ !-----------------------------------------------------------------------------------------
+ ! ISURFACE=0, IMODEL_SYN : homo 1D model checker het
+ ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert het pert
+ !-----------------------------------------------------------------------------------------
+
+ ! whether you want to smooth the structure model for computing the synthetics
+ integer, parameter :: ISMOOTH_INITIAL_MODEL = 0
+ integer, parameter :: ISMOOTH_MODEL_UPDATE = 0
+ integer, parameter :: ISMOOTH_EVENT_KERNEL = 1
+ integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
+
+ ! smoothing scalelength
+ double precision, parameter :: GAMMA_SMOOTH_MODEL = 30.0e3
+
+!========================
+
+! boolean parameters
+! IKER: (0) waveform
+! (1) traveltime, cross-correlation, misfit
+! (2) amplitude, cross-correlation, misfit
+! (3) traveltime, multitaper
+! (4) amplitude, multitaper
+! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
+! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
+ integer, parameter :: IKER = 1
+ integer, parameter :: IAMP_VEL = 0 ! measure amplitudes between velocity traces
+
+! KEY: USE ONE LINE OR THE OTHER
+ integer, parameter :: ISURFACE = 1, NCOMP = 1, NABSORB = 4 ! surface waves
+! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
+
+! parameters controlling what to write to file
+! NOTE: for the tomography simulations, ALL of these can be .false.
+
+ logical, parameter :: WRITE_STF_F = .false.
+ logical, parameter :: WRITE_SEISMO_F = .false.
+! logical, parameter :: WRITE_SPECTRA_F = .false.
+! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
+
+ logical, parameter :: WRITE_SEISMO_RECONSTRUCT = .false. ! multitaper
+! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
+
+ logical, parameter :: WRITE_STF_A = .false.
+ logical, parameter :: WRITE_SEISMO_A = .false. ! source inversions
+! logical, parameter :: WRITE_SPECTRA_A = .false.
+! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
+
+ logical, parameter :: WRITE_KERNELS = .true. ! write all nine kernels
+ logical, parameter :: WRITE_KERNEL_SNAPSHOTS = .false. ! kernel snapshots
+ logical, parameter :: WRITE_WAVFIELD_SNAPSHOTS = .false. ! wavefield snapshots
+
+!--------------------------------------
+! INVERSION PARAMETERS
+
+ ! whether you want to compute kernels, or simply the misfit function
+ logical, parameter :: COMPUTE_KERNELS = .true.
+
+ ! whether to use the data subspace method, which has a Hessian
+ logical, parameter :: HESSIAN = .false.
+
+ ! stopping criteria
+ ! NITERATION : number of iterations
+ ! VAR_RED_MIN : minimum variance reduction (in percent)
+ ! SIGMA_FAC : stop if a model value exceeds SIGMA_FAC * sigma_m
+ ! CONV_STOP : stop when the misfit value is this fraction of the INITIAL misfit value
+ integer, parameter :: NITERATION = 16
+ double precision, parameter :: VAR_RED_MIN = 8.0
+
+ ! Gaussian errors containted in input file
+ ! see wave2d_sigmas.m and INPUT/sigma_0p1_pert.dat
+ double precision, parameter :: SIGMA_DT = 0.10
+ double precision, parameter :: SIGMA_DLNA = 1.0
+ double precision, parameter :: SIGMA_WAVEFORM = 1.0
+ logical, parameter :: ADD_DATA_ERRORS = .true.
+
+ !double precision, parameter :: SIGMA_FAC = 2.0
+ !double precision, parameter :: CONV_STOP = 1.0e-4
+
+ ! order of interpolating polynomial in conjugate gradient algorithm
+ ! using POLY_ORDER = 3 required computing the gradient of the test model
+ integer, parameter :: POLY_ORDER = 2 ! 2 (preferred) or 3
+
+ ! use inversion and scaling used for the GJI paper
+ !integer, parameter :: GJI_PAPER = 1
+
+ ! what to perturb, what to invert
+ integer, parameter :: PERT_STRUCT = 1
+ integer, parameter :: PERT_SOURCE_T = 0
+ integer, parameter :: PERT_SOURCE_X = 0
+
+ integer, parameter :: INV_STRUCT = 1
+ integer, parameter :: INV_SOURCE_T = 0
+ integer, parameter :: INV_SOURCE_X = 0
+
+ ! whether to include the model norm term in the misfit function, which acts like damping
+ logical, parameter :: INCLUDE_MODEL_NORM = .true.
+
+ ! log file showing source loactions
+ logical, parameter :: ISOURCE_LOG = .true.
+
+ ! DO NOT CHANGE THESE
+ integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ integer, parameter :: NVAR_SOURCE = 3 ! x position, z position, origin time
+ integer, parameter :: NVAR = NVAR_STRUCT + NVAR_SOURCE
+
+! parameterization for structure in the inversion
+! 1 : kappa-mu-rho
+! 2 : alpha-beta-rho
+! 3 : c-beta-rho
+ integer, parameter :: STRUCTURE_PARAMETER_TYPE = 2
+
+! scalelength of smoothing Gaussian
+! GJI paper : 30,60,90
+! body waves : 10
+ double precision, parameter :: GAMMA_SMOOTH = 60.0e3
+
+ ! parameters for smoothing
+ logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
+ logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
+
+! homogeneous background model (S.I. units)
+ double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
+ double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
+ double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
+
+!---------------------------------------------------------------
+! measurement windows
+
+ double precision, parameter :: HWIN1 = 2.5*hdur ! half-window width for pulse (2.5)
+ double precision, parameter :: HWIN2 = 5.0*hdur ! half-window with for MTM measurement (zero-padding)
+
+!---------------------------------------------------------------
+! CHT: do not change these
+
+! UTM zone for Southern California region
+! integer, parameter :: UTM_PROJECTION_ZONE = 11
+
+! to suppress UTM projection for SCEC benchmarks
+ logical, parameter :: SUPPRESS_UTM_PROJECTION = .false.
+
+! flag for projection from latitude/longitude to UTM, and back
+ integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
+
+! flag for projection from latitude/longitude to mesh-UTM, and back
+ integer, parameter :: ILONLAT2MESH = 0, IMESH2LONLAT = 1
+
+! max number of fake receivers
+ integer, parameter :: MAX_SR_FAKE = 1000
+
+! max number of events, receivers, and phass
+ integer, parameter :: MAX_EVENT = 50
+ integer, parameter :: MAX_SR = 1400
+ integer, parameter :: MAX_PHASE = 1
+ integer, parameter :: MAX_COMP = NCOMP
+
+!========================
+! GRID AND GLL POINTS
+
+ integer, parameter :: NELE = MAX(NEX,NEZ)
+ integer, parameter :: NSPEC = NEX*NEZ
+
+! number of GLL points (polynomial degree plus one)
+ integer, parameter :: NGLLX = 5
+ integer, parameter :: NGLLZ = 5
+ integer, parameter :: NGLL = MAX(NGLLX,NGLLZ)
+
+! number of global points
+ integer, parameter :: NGLLSQUARE = NGLLX * NGLLZ ! element GLL points
+ integer, parameter :: NGLOB = ((NGLLX-1)*NEX + 1)*((NGLLZ-1)*NEZ +1) ! global GLL points
+ integer, parameter :: NLOCAL = NGLLX * NGLLZ * NSPEC ! local GLL points
+ integer, parameter :: NSPEC_CORNER = (NEX+1) * (NEZ+1) ! element corner points
+
+! number of nodes for 2D and 3D shape functions for hexahedra
+! we use 8-node mesh bricks, which are more stable than 27-node elements
+ integer, parameter :: NGNOD = 8, NGNOD2D = 4
+
+! number of iterations to solve the system for xi and eta
+ integer, parameter :: NUM_ITER = 1
+
+! very large and very small values
+ double precision, parameter :: HUGEVAL = 1.0e30, TINYVAL = 1.0e-9
+
+! for the Gauss-Lobatto-Legendre points and weights
+ double precision, parameter :: GAUSSALPHA = 0.0,GAUSSBETA = 0.0
+
+!
+! CONSTANTS
+!
+ double precision, parameter :: PI = 3.141592653589793
+ double precision, parameter :: FOUR_THIRDS = 4.0 / 3.0
+ double precision, parameter :: ONE_THIRD = 1.0 / 3.0
+ double precision, parameter :: ONEOVERTWO = 0.5
+ double precision, parameter :: DEG = 180.0/PI
+! double precision, parameter :: EPS = 1.0e-35
+
+!!$! parameter for FFTW
+!!$ integer, parameter :: NOUT = NSTEP/2 + 1
+!!$
+!!$! bounds for bandpass filter
+!!$ double precision, parameter :: hwid = 3.0 ! HALF-width of window (s)
+!!$ double precision, parameter :: tmin = 2.0*hdur-hwid
+!!$ double precision, parameter :: tmax = 2.0*hdur+hwid
+!!$ double precision, parameter :: fmin = 1.0/tmax, fmax = 1.0/tmin
+!!$ double precision, parameter :: trbdndw = 0.3, a = 30.0
+!!$ integer, parameter :: passes = 2, iord = 4
+!!$
+!!$!
+!!$! MULTI-TAPER PARAMETERS
+!!$!
+!!$! Ying Zhou: The fit between the recovered data and the data can be improved
+!!$! by either increasing the window width (HWIN above) or by decreasing NPI.
+!!$! In her experience, NPI = 2.5 is good for noisy data.
+!!$! For synthetic data, we can use a lower NPI.
+!!$! number of tapers should be fixed as twice NPI -- see Latex notes
+!!$ double precision, parameter :: WTR = 0.02
+!!$ double precision, parameter :: NPI = 2.5
+!!$ integer, parameter :: NTAPER = int(2.0*NPI)
+!!$
+!!$! FFT parameters
+!!$! FOR SOME REASON, I CANNOT SET NDIM = npt when lnpt <= 11 -- WHY NOT?
+!!$ integer, parameter :: lnpt = 13, npt = 2**lnpt, NDIM = 20000
+!!$ !integer, parameter :: lnpt = 14, npt = 2**lnpt, NDIM = npt
+!!$
+!!$ !double precision, parameter :: ZZIGN = -1.0 ! sign convention -- WATCH OUT
+!!$ double precision, parameter :: FORWARD_FFT = 1.0
+!!$ double precision, parameter :: REVERSE_FFT = -1.0
+
+ integer, parameter :: NDIM = 20000
+
+end module wave2d_constants
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_solver.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_solver.f90 2010-01-15 21:10:53 UTC (rev 16134)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_solver.f90 2010-01-15 22:07:21 UTC (rev 16135)
@@ -375,10 +375,9 @@
stop 'check model in set_model_property.f90'
else
- stop 'check model in set_model_property.f90'
- !kappa_syn(i,j,ispec) = INCOMPRESSIBILITY ! not used
- !mu_syn(i,j,ispec) = DENSITY * beta_syn(iglob)*beta_syn(iglob)
- !rho_syn(i,j,ispec) = DENSITY
+ kappa_syn(i,j,ispec) = INCOMPRESSIBILITY ! not used
+ mu_syn(i,j,ispec) = DENSITY * beta_syn(i,j,ispec)*beta_syn(i,j,ispec)
+ rho_syn(i,j,ispec) = DENSITY
endif
More information about the CIG-COMMITS
mailing list