debug alpha_beta_rho

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(

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]

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]

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

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))

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

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

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

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)