From e44bbb6a8cde65d2213c8379d757d9a31927f928 Mon Sep 17 00:00:00 2001 From: Mijian Xu Date: Mon, 6 Jan 2025 22:14:43 +0800 Subject: [PATCH] debug alpha_beta_rho --- CMakeLists.txt | 2 +- examples/00_checkerboard_iso/input_params.yml | 3 +++ examples/04_hawaii_topo/input_params.yml | 3 +++ examples/04_hawaii_topo/run_this_example.sh | 2 +- src/model.f90 | 16 +++++++++++++--- src/optimize.f90 | 9 ++++++--- src/shared/argparse.f90 | 12 +++++++++--- src/surfatt_cb_fwd.f90 | 5 +++-- src/tomo.f90 | 7 +++++-- 9 files changed, 44 insertions(+), 15 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 373559f..c2678c9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.10) # set the project name set(CMAKE_PROJECT_NAME "SurfATT") -project(${CMAKE_PROJECT_NAME} VERSION 1.5.1 LANGUAGES Fortran CXX) +project(${CMAKE_PROJECT_NAME} VERSION 1.6.0 LANGUAGES Fortran CXX) # set the version number and git commit hash execute_process( diff --git a/examples/00_checkerboard_iso/input_params.yml b/examples/00_checkerboard_iso/input_params.yml index 867ddb8..fe8dd21 100644 --- a/examples/00_checkerboard_iso/input_params.yml +++ b/examples/00_checkerboard_iso/input_params.yml @@ -28,6 +28,9 @@ topo: wavelen_factor: 2.5 inversion: + # ------- Inversion parameters -------- + use_alpha_beta_rho: false + rho_scaling: false # ------- Initial model -------- init_model_type: 0 # 0: increase from v0 to v1; 1: 1D inversion for average tt; 2: use a 3D vs model vel_range: [1.8, 4.2] diff --git a/examples/04_hawaii_topo/input_params.yml b/examples/04_hawaii_topo/input_params.yml index a3d998f..d33949c 100644 --- a/examples/04_hawaii_topo/input_params.yml +++ b/examples/04_hawaii_topo/input_params.yml @@ -22,6 +22,9 @@ topo: wavelen_factor: 2.5 inversion: + # ------- Inversion parameters -------- + use_alpha_beta_rho: true + rho_scaling: true # ------- Initial model -------- init_model_type: 1 # 0: increase from v0 to v1; 1: 1D inversion for average tt; 2: use a 3D vs model vel_range: [1.8, 4.2] diff --git a/examples/04_hawaii_topo/run_this_example.sh b/examples/04_hawaii_topo/run_this_example.sh index 0861a2e..3e76b4f 100644 --- a/examples/04_hawaii_topo/run_this_example.sh +++ b/examples/04_hawaii_topo/run_this_example.sh @@ -6,7 +6,7 @@ pos_str=19.5/-155.5 angle=-30 # download topography -gmt grdcut @earth_relief_01m -R-157/-152/18/21 -Ghawaii.nc +# gmt grdcut @earth_relief_01m -R-157/-152/18/21 -Ghawaii.nc # rotate source receiver file ../../bin/surfatt_rotate_src_rec -i src_rec_file_raw.csv -a $angle -c $pos_str -o src_rec_file_rotated.csv diff --git a/src/model.f90 b/src/model.f90 index 0353e62..0b2fbb4 100644 --- a/src/model.f90 +++ b/src/model.f90 @@ -236,14 +236,16 @@ subroutine inv1d(this, vsinv, niter, misfits) niter = iter end subroutine inv1d - subroutine add_pert(this,nx,ny,nz,pert_vel,hmarg,anom_size) + subroutine add_pert(this,nx,ny,nz,pert_vel,hmarg,anom_size,only_vs) class(att_model), intent(inout) :: this integer, intent(in) :: nx, ny, nz real(kind=dp), intent(in), optional :: pert_vel,hmarg,anom_size + logical, intent(in), optional :: only_vs real(kind=dp) :: ashmarg,aspert_vel,amp, asanom_size real(kind=dp), dimension(:), allocatable :: x_pert, y_pert, z_pert real(kind=dp), dimension(:,:,:), allocatable :: vs_pert real(kind=dp), dimension(3) :: para + logical :: asonly_vs integer :: ntaperx, ntapery, i, j, k if (present(pert_vel)) then @@ -264,6 +266,12 @@ subroutine add_pert(this,nx,ny,nz,pert_vel,hmarg,anom_size) asanom_size = 0. endif + if (present(only_vs)) then + asonly_vs = only_vs + else + asonly_vs = .false. + endif + if (myrank == 0) then ntaperx = ashmarg/this%d_xyz(1) ntapery = ashmarg/this%d_xyz(2) @@ -289,8 +297,10 @@ subroutine add_pert(this,nx,ny,nz,pert_vel,hmarg,anom_size) enddo enddo this%vs3d = this%vs3d*(1+vs_pert) - this%vp3d = empirical_vp(this%vs3d) - this%rho3d = empirical_rho(this%vp3d) + if (.not. asonly_vs) then + this%vp3d = empirical_vp(this%vs3d) + this%rho3d = empirical_rho(this%vp3d) + endif endif call synchronize_all() call sync_from_main_rank(this%vs3d, this%n_xyz(1), this%n_xyz(2), this%n_xyz(3)) diff --git a/src/optimize.f90 b/src/optimize.f90 index f0d9278..777dc38 100644 --- a/src/optimize.f90 +++ b/src/optimize.f90 @@ -8,6 +8,9 @@ ! (c) October 2023 ! ! Changing History: Apl 2024, Initialize Codes +! May 2024, Add L-BFGS and CG optimization +! Jun 2024, Add Hager-Zhang formula for CG +! Jan 2025, change dims for Alpha_beta_rho ! !===================================================================== module optimize @@ -35,8 +38,8 @@ subroutine get_lbfgs_direction(iter, direction) call get_gradient(iter, q_vector) dims = shape(q_vector) - allocate(gradient_diff(nstore, nker, dims(1), dims(2), dims(3))) - allocate(model_diff(nstore, nker, dims(1), dims(2), dims(3))) + allocate(gradient_diff(nstore, nker, dims(2), dims(3), dims(4))) + allocate(model_diff(nstore, nker, dims(2), dims(3), dims(4))) idx_iter = zeros(nstore) p = zeros(nstore) a = zeros(nstore) @@ -97,7 +100,7 @@ subroutine get_model(iter, model) write(key_name, '("/vs_",I3.3)') iter call h5read(model_fname, key_name, tmp_model) tmp_model = transpose_3(tmp_model) - model = zeros( nker,size(tmp_model, 1), size(tmp_model, 3), size(tmp_model, 3)) + model = zeros( nker,size(tmp_model, 1), size(tmp_model, 2), size(tmp_model, 3)) model(1,:,:,:) = tmp_model if (ap%inversion%use_alpha_beta_rho) then ! read vp model diff --git a/src/shared/argparse.f90 b/src/shared/argparse.f90 index ad2e842..c768e5c 100644 --- a/src/shared/argparse.f90 +++ b/src/shared/argparse.f90 @@ -61,18 +61,20 @@ subroutine argparse_tomo(fname, isfwd) endif end subroutine argparse_tomo - subroutine argparse_cb_fwd(fname, ncb, pert_vel, hmarg, anom_size, max_noise) + subroutine argparse_cb_fwd(fname, ncb, pert_vel, hmarg, anom_size, max_noise, only_vs) character(len=MAX_STRING_LEN),dimension(:), allocatable :: args character(len=MAX_STRING_LEN) :: arg, value character(len=MAX_STRING_LEN),intent(out) :: fname integer, dimension(3), intent(out) :: ncb real(kind=dp), intent(out) :: pert_vel, hmarg, anom_size, max_noise - integer :: i,nargs,m,ier, nopt=6 + logical, intent(out) :: only_vs + integer :: i,nargs,m,ier, nopt=7 pert_vel = 0.08 hmarg = 0. anom_size = 0. max_noise = 0. + only_vs = .false. nargs = command_argument_count() allocate(args(nargs)) do i = 1, nargs @@ -91,6 +93,7 @@ subroutine argparse_cb_fwd(fname, ncb, pert_vel, hmarg, anom_size, max_noise) write(*,*)'optional arguments:' write(*,*)' -h Print help message' write(*,*)' -e tt_noise Add random noise to travel time data, defaults to 0' + write(*,*)' -v Only add perturbations on Vs model' write(*,*)' -m margin_degree Margin area in degree between anomaly as boundary, defaults to 0' write(*,*)' -p pert_vel Magnitude of velocity perturbations, defaults to 0.08' write(*,*)' -s anom_size_km size of anomalies at the top in km, default to uniform anomaly size in Z direction' @@ -121,7 +124,10 @@ subroutine argparse_cb_fwd(fname, ncb, pert_vel, hmarg, anom_size, max_noise) elseif(arg(1:2) == '-e') then m = m+1 read(args(i+1),*,iostat=ier) max_noise - if(ier/=0) stop 'Cannot parse argument -e' + if(ier/=0) stop 'Cannot parse argument -e' + elseif(arg(1:2) == '-v') then + m = m+1 + only_vs = .true. endif enddo if (m<1) then diff --git a/src/surfatt_cb_fwd.f90 b/src/surfatt_cb_fwd.f90 index 9150943..c6d7c47 100644 --- a/src/surfatt_cb_fwd.f90 +++ b/src/surfatt_cb_fwd.f90 @@ -28,12 +28,13 @@ program surfatt_cb_fwd character(len=MAX_STRING_LEN) :: fname integer, dimension(3) :: ncb real(kind=dp) :: pert, hmarg, anom_size, max_noise + logical :: only_vs ! initialize MPI call init_mpi() ! read command line arguments - call argparse_cb_fwd(fname, ncb, pert, hmarg, anom_size, max_noise) + call argparse_cb_fwd(fname, ncb, pert, hmarg, anom_size, max_noise, only_vs) ! read parameter file call ap%read(fname) @@ -54,7 +55,7 @@ program surfatt_cb_fwd ! add perturbations call am%add_pert(ncb(1),ncb(2),ncb(3), pert_vel=pert,& - hmarg=hmarg, anom_size=anom_size) + hmarg=hmarg, anom_size=anom_size,only_vs=only_vs) ! initialize grid if (ap%data%vel_type(1)) then diff --git a/src/tomo.f90 b/src/tomo.f90 index ae52adb..2be3642 100644 --- a/src/tomo.f90 +++ b/src/tomo.f90 @@ -8,6 +8,7 @@ ! (c) October 2023 ! ! Changing History: Oct 2023, Initialize Codes +! Jan 2025, change dims for Alpha_beta_rho ! !===================================================================== module tomo @@ -84,6 +85,10 @@ subroutine initialize_inv(this) call h%add('/lat', am%ygrids) call h%add('/dep', am%zgrids) call h%add('/vs_000', transpose_3(am%vs3d)) + if (ap%inversion%use_alpha_beta_rho) then + call h%add('/vp_000', transpose_3(am%vp3d)) + call h%add('/rho_000', transpose_3(am%rho3d)) + endif call h%close(finalize=.true.) endif endif @@ -229,7 +234,6 @@ subroutine steepest_descent(this) am%vp3d = am%vp3d * (1 + gradient_s(2,:,:,:)) am%rho3d = am%rho3d * (1 + gradient_s(3,:,:,:)) endif - !!!!!!!!!!!!!!!!!!!!!!! call write_tmp_model() endif call synchronize_all() @@ -305,7 +309,6 @@ subroutine prepare_fwd_linesearch() am%vp3d_opt = am%vp3d * (1 + gradient_ls(2,:,:,:)) am%rho3d_opt = am%rho3d * (1 + gradient_ls(3,:,:,:)) else - gradient_ls = updatemax * direction / max_gk am%vs3d_opt = am%vs3d * (1 + gradient_ls(1,:,:,:)) am%vp3d_opt = empirical_vp(am%vs3d_opt) am%rho3d_opt = empirical_rho(am%vp3d_opt)