Skip to content

Commit

Permalink
preparing to add checks of all allocate() statements using a script: …
Browse files Browse the repository at this point in the history 
…step 3
  • Loading branch information
@komatits
komatits committed Jun 16, 2018
1 parent 936429a commit 70954a5
Show file tree
Hide file tree
Showing 25 changed files with 782 additions and 425 deletions.
6 changes: 4 additions & 2 deletions src/decompose_mesh/program_decompose_mesh_mpi.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,7 +111,9 @@ program xdecompose_mesh_mpi
endif

call bcast_all_singlei(nspec_glob)
if (myrank > 0) allocate(ipart(nspec_glob))
if (myrank > 0) then
allocate(ipart(nspec_glob))
endif
call bcast_all_i(ipart, nspec_glob)
call send_partition_mesh_to_all(myrank, ipart, npartX*npartY*npartZ)
call send_mesh_to_all(myrank)
Expand Down Expand Up @@ -205,7 +207,7 @@ subroutine send_partition_mesh_to_all(myrank, ipart, npart)
! compute elments by node table
max_elmnts_by_node = maxval(nelmnts_by_node_glob)
nelmnts_by_node_glob(:)=0
allocate(elmnts_by_node_glob(max_elmnts_by_node, nnodes_glob)) !!
allocate(elmnts_by_node_glob(max_elmnts_by_node, nnodes_glob))
elmnts_by_node_glob(:,:)=-1
do iE=1,nE
do inode=1,NGNOD
Expand Down
39 changes: 29 additions & 10 deletions src/inverse_problem_for_model/adjoint_source/adjoint_source_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -157,10 +157,17 @@ end subroutine deallocate_adjoint_source_working_arrays
!----------------------------------------------------------------------------------------------------------------------------------
subroutine allocate_adjoint_source_working_arrays()

allocate(residuals(nstep_data), raw_residuals(nstep_data), fil_residuals(nstep_data), filfil_residuals(nstep_data), &
w_tap(nstep_data), signal(nstep_data), residuals_for_cost(nstep_data), &
elastic_adjoint_source(NDIM,nstep_data), elastic_misfit(NDIM,nstep_data), &
data_trace_to_use(nstep_data), wkstmp(nstep_data))
allocate(residuals(nstep_data))
allocate(raw_residuals(nstep_data))
allocate(fil_residuals(nstep_data))
allocate(filfil_residuals(nstep_data))
allocate(w_tap(nstep_data))
allocate(signal(nstep_data))
allocate(residuals_for_cost(nstep_data))
allocate(elastic_adjoint_source(NDIM,nstep_data))
allocate(elastic_misfit(NDIM,nstep_data))
allocate(data_trace_to_use(nstep_data))
allocate(wkstmp(nstep_data))

end subroutine allocate_adjoint_source_working_arrays

Expand Down Expand Up @@ -219,10 +226,18 @@ subroutine compute_elastic_adjoint_source_displacement(irec_local, ievent, curre
case ('L2_FWI_TELESEISMIC')

! Define temporary trace vector
if (.not. allocated(trace_cal_1)) allocate(trace_cal_1(3,nstep_data))
if (.not. allocated(trace_cal_2)) allocate(trace_cal_2(3,nstep_data))
if (.not. allocated(trace_obs_1)) allocate(trace_obs_1(3,nstep_data))
if (.not. allocated(trace_obs_2)) allocate(trace_obs_2(3,nstep_data))
if (.not. allocated(trace_cal_1)) then
allocate(trace_cal_1(3,nstep_data))
endif
if (.not. allocated(trace_cal_2)) then
allocate(trace_cal_2(3,nstep_data))
endif
if (.not. allocated(trace_obs_1)) then
allocate(trace_obs_1(3,nstep_data))
endif
if (.not. allocated(trace_obs_2)) then
allocate(trace_obs_2(3,nstep_data))
endif
lat0 = acqui_simu(ievent)%Origin_chunk_lat
lon0 = acqui_simu(ievent)%Origin_chunk_lon
azi0 = acqui_simu(ievent)%Origin_chunk_azi
Expand All @@ -239,7 +254,9 @@ subroutine compute_elastic_adjoint_source_displacement(irec_local, ievent, curre

! Convolve synthetic data with wavelet
if (inversion_param%convolution_by_wavelet) then
if (.not. allocated(wavelet)) allocate(wavelet(nstep_data))
if (.not. allocated(wavelet)) then
allocate(wavelet(nstep_data))
endif
wavelet = acqui_simu(ievent)%user_source_time_function(1,:)
call myconvolution(trace_cal_2(idim,:),wavelet,nstep_data,nstep_data,tmpl,0)
trace_cal_1(idim,:) = tmpl * dt_data
Expand Down Expand Up @@ -322,7 +339,9 @@ subroutine compute_elastic_adjoint_source_displacement(irec_local, ievent, curre

! Finally cross-correlate residuals with wavelet
if (inversion_param%convolution_by_wavelet) then
if (.not. allocated(filfil_residuals_tmp)) allocate(filfil_residuals_tmp(nstep_data))
if (.not. allocated(filfil_residuals_tmp)) then
allocate(filfil_residuals_tmp(nstep_data))
endif
filfil_residuals_tmp(:) = filfil_residuals(:)
call mycorrelation(filfil_residuals_tmp,wavelet,nstep_data,nstep_data,tmpl,0)
filfil_residuals = tmpl * dt_data
Expand Down
237 changes: 176 additions & 61 deletions src/inverse_problem_for_model/adjoint_source/instantaneous_phase_envelope_misfit_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -99,42 +99,83 @@ subroutine compute_instantaneous_phase_and_envelope_data(giter)

endif

! call MPI_allreduce(MPI_IN_PLACE,amp_wlx,1,MPI_REAL,MPI_SUM,comm%mpi_comm_1,ierr_mpi)
! call MPI_allreduce(MPI_IN_PLACE,amp_wly,1,MPI_REAL,MPI_SUM,comm%mpi_comm_1,ierr_mpi)
! call MPI_allreduce(MPI_IN_PLACE,amp_wlx,1,MPI_REAL,MPI_SUM,comm%mpi_comm_1,ierr_mpi)
! call MPI_allreduce(MPI_IN_PLACE,amp_wly,1,MPI_REAL,MPI_SUM,comm%mpi_comm_1,ierr_mpi)
call MPI_allreduce(MPI_IN_PLACE,amp_wlz,1,MPI_REAL,MPI_SUM,comm%mpi_comm_1,ierr_mpi)

endif

if (nrecloc > 0) then

if (.not. allocated(envelo_vx)) allocate(envelo_vx(nrecloc,nt))
if (.not. allocated(envelo_vy)) allocate(envelo_vy(nrecloc,nt))
if (.not. allocated(envelo_vz)) allocate(envelo_vz(nrecloc,nt))

if (.not. allocated(envelc_vx)) allocate(envelc_vx(nrecloc,nt))
if (.not. allocated(envelc_vy)) allocate(envelc_vy(nrecloc,nt))
if (.not. allocated(envelc_vz)) allocate(envelc_vz(nrecloc,nt))
if (.not. allocated(envelo_vx)) then
allocate(envelo_vx(nrecloc,nt))
endif
if (.not. allocated(envelo_vy)) then
allocate(envelo_vy(nrecloc,nt))
endif
if (.not. allocated(envelo_vz)) then
allocate(envelo_vz(nrecloc,nt))
endif

if (.not. allocated(denvel_vx)) allocate(denvel_vx(nrecloc,nt))
if (.not. allocated(denvel_vy)) allocate(denvel_vy(nrecloc,nt))
if (.not. allocated(denvel_vz)) allocate(denvel_vz(nrecloc,nt))
if (.not. allocated(envelc_vx)) then
allocate(envelc_vx(nrecloc,nt))
endif
if (.not. allocated(envelc_vy)) then
allocate(envelc_vy(nrecloc,nt))
endif
if (.not. allocated(envelc_vz)) then
allocate(envelc_vz(nrecloc,nt))
endif

if (.not. allocated(dphase_vx)) allocate(dphase_vx(nrecloc,nt))
if (.not. allocated(dphase_vy)) allocate(dphase_vy(nrecloc,nt))
if (.not. allocated(dphase_vz)) allocate(dphase_vz(nrecloc,nt))
if (.not. allocated(denvel_vx)) then
allocate(denvel_vx(nrecloc,nt))
endif
if (.not. allocated(denvel_vy)) then
allocate(denvel_vy(nrecloc,nt))
endif
if (.not. allocated(denvel_vz)) then
allocate(denvel_vz(nrecloc,nt))
endif

if (.not. allocated(danalytic_vx)) allocate(danalytic_vx(nrecloc,nt))
if (.not. allocated(danalytic_vy)) allocate(danalytic_vy(nrecloc,nt))
if (.not. allocated(danalytic_vz)) allocate(danalytic_vz(nrecloc,nt))
if (.not. allocated(dphase_vx)) then
allocate(dphase_vx(nrecloc,nt))
endif
if (.not. allocated(dphase_vy)) then
allocate(dphase_vy(nrecloc,nt))
endif
if (.not. allocated(dphase_vz)) then
allocate(dphase_vz(nrecloc,nt))
endif

if (.not. allocated(an_dobs_vx)) allocate(an_dobs_vx(nrecloc,nt))
if (.not. allocated(an_dobs_vy)) allocate(an_dobs_vy(nrecloc,nt))
if (.not. allocated(an_dobs_vz)) allocate(an_dobs_vz(nrecloc,nt))
if (.not. allocated(danalytic_vx)) then
allocate(danalytic_vx(nrecloc,nt))
endif
if (.not. allocated(danalytic_vy)) then
allocate(danalytic_vy(nrecloc,nt))
endif
if (.not. allocated(danalytic_vz)) then
allocate(danalytic_vz(nrecloc,nt))
endif

if (.not. allocated(an_dcal_vx)) allocate(an_dcal_vx(nrecloc,nt))
if (.not. allocated(an_dcal_vy)) allocate(an_dcal_vy(nrecloc,nt))
if (.not. allocated(an_dcal_vz)) allocate(an_dcal_vz(nrecloc,nt))
if (.not. allocated(an_dobs_vx)) then
allocate(an_dobs_vx(nrecloc,nt))
endif
if (.not. allocated(an_dobs_vy)) then
allocate(an_dobs_vy(nrecloc,nt))
endif
if (.not. allocated(an_dobs_vz)) then
allocate(an_dobs_vz(nrecloc,nt))
endif

if (.not. allocated(an_dcal_vx)) then
allocate(an_dcal_vx(nrecloc,nt))
endif
if (.not. allocated(an_dcal_vy)) then
allocate(an_dcal_vy(nrecloc,nt))
endif
if (.not. allocated(an_dcal_vz)) then
allocate(an_dcal_vz(nrecloc,nt))
endif

!*** 2. Get their analytic signal
call get_analytic_signal
Expand Down Expand Up @@ -196,25 +237,55 @@ subroutine compute_instanteneous_phase_adjoint_source_term

if (nrecloc > 0) then

if (.not. allocated(srcterm1_vx)) allocate(srcterm1_vx(nrecloc,nt))
if (.not. allocated(srcterm1_vy)) allocate(srcterm1_vy(nrecloc,nt))
if (.not. allocated(srcterm1_vz)) allocate(srcterm1_vz(nrecloc,nt))
if (.not. allocated(srcterm1_vx)) then
allocate(srcterm1_vx(nrecloc,nt))
endif
if (.not. allocated(srcterm1_vy)) then
allocate(srcterm1_vy(nrecloc,nt))
endif
if (.not. allocated(srcterm1_vz)) then
allocate(srcterm1_vz(nrecloc,nt))
endif

if (.not. allocated(srcterm2_vx)) allocate(srcterm2_vx(nrecloc,nt))
if (.not. allocated(srcterm2_vy)) allocate(srcterm2_vy(nrecloc,nt))
if (.not. allocated(srcterm2_vz)) allocate(srcterm2_vz(nrecloc,nt))
if (.not. allocated(srcterm2_vx)) then
allocate(srcterm2_vx(nrecloc,nt))
endif
if (.not. allocated(srcterm2_vy)) then
allocate(srcterm2_vy(nrecloc,nt))
endif
if (.not. allocated(srcterm2_vz)) then
allocate(srcterm2_vz(nrecloc,nt))
endif

if (.not. allocated(srcterm2tmp_vx)) allocate(srcterm2tmp_vx(nrecloc,nt))
if (.not. allocated(srcterm2tmp_vy)) allocate(srcterm2tmp_vy(nrecloc,nt))
if (.not. allocated(srcterm2tmp_vz)) allocate(srcterm2tmp_vz(nrecloc,nt))
if (.not. allocated(srcterm2tmp_vx)) then
allocate(srcterm2tmp_vx(nrecloc,nt))
endif
if (.not. allocated(srcterm2tmp_vy)) then
allocate(srcterm2tmp_vy(nrecloc,nt))
endif
if (.not. allocated(srcterm2tmp_vz)) then
allocate(srcterm2tmp_vz(nrecloc,nt))
endif

if (.not. allocated(ft_tmp_vx)) allocate(ft_tmp_vx(nrecloc,nt))
if (.not. allocated(ft_tmp_vy)) allocate(ft_tmp_vy(nrecloc,nt))
if (.not. allocated(ft_tmp_vz)) allocate(ft_tmp_vz(nrecloc,nt))
if (.not. allocated(ft_tmp_vx)) then
allocate(ft_tmp_vx(nrecloc,nt))
endif
if (.not. allocated(ft_tmp_vy)) then
allocate(ft_tmp_vy(nrecloc,nt))
endif
if (.not. allocated(ft_tmp_vz)) then
allocate(ft_tmp_vz(nrecloc,nt))
endif

if (.not. allocated(IP_adjt_vx)) allocate(IP_adjt_vx(nrecloc,nt))
if (.not. allocated(IP_adjt_vy)) allocate(IP_adjt_vy(nrecloc,nt))
if (.not. allocated(IP_adjt_vz)) allocate(IP_adjt_vz(nrecloc,nt))
if (.not. allocated(IP_adjt_vx)) then
allocate(IP_adjt_vx(nrecloc,nt))
endif
if (.not. allocated(IP_adjt_vy)) then
allocate(IP_adjt_vy(nrecloc,nt))
endif
if (.not. allocated(IP_adjt_vz)) then
allocate(IP_adjt_vz(nrecloc,nt))
endif

!*** 1st term (easy...) remeber, imag(analytic_sig) = hilbert transform
srcterm1_vx = dphase_vx * aimag(danalytic_vx) / (envelc_vx**2 + wlx**2)
Expand Down Expand Up @@ -323,25 +394,55 @@ subroutine compute_envelope_adjoint_source_term

if (nrecloc > 0) then

if (.not. allocated(srcterm1_vx)) allocate(srcterm1_vx(nrecloc,nt))
if (.not. allocated(srcterm1_vy)) allocate(srcterm1_vy(nrecloc,nt))
if (.not. allocated(srcterm1_vz)) allocate(srcterm1_vz(nrecloc,nt))
if (.not. allocated(srcterm1_vx)) then
allocate(srcterm1_vx(nrecloc,nt))
endif
if (.not. allocated(srcterm1_vy)) then
allocate(srcterm1_vy(nrecloc,nt))
endif
if (.not. allocated(srcterm1_vz)) then
allocate(srcterm1_vz(nrecloc,nt))
endif

if (.not. allocated(srcterm2_vx)) allocate(srcterm2_vx(nrecloc,nt))
if (.not. allocated(srcterm2_vy)) allocate(srcterm2_vy(nrecloc,nt))
if (.not. allocated(srcterm2_vz)) allocate(srcterm2_vz(nrecloc,nt))
if (.not. allocated(srcterm2_vx)) then
allocate(srcterm2_vx(nrecloc,nt))
endif
if (.not. allocated(srcterm2_vy)) then
allocate(srcterm2_vy(nrecloc,nt))
endif
if (.not. allocated(srcterm2_vz)) then
allocate(srcterm2_vz(nrecloc,nt))
endif

if (.not. allocated(srcterm2tmp_vx)) allocate(srcterm2tmp_vx(nrecloc,nt))
if (.not. allocated(srcterm2tmp_vy)) allocate(srcterm2tmp_vy(nrecloc,nt))
if (.not. allocated(srcterm2tmp_vz)) allocate(srcterm2tmp_vz(nrecloc,nt))
if (.not. allocated(srcterm2tmp_vx)) then
allocate(srcterm2tmp_vx(nrecloc,nt))
endif
if (.not. allocated(srcterm2tmp_vy)) then
allocate(srcterm2tmp_vy(nrecloc,nt))
endif
if (.not. allocated(srcterm2tmp_vz)) then
allocate(srcterm2tmp_vz(nrecloc,nt))
endif

if (.not. allocated(ft_tmp_vx)) allocate(ft_tmp_vx(nrecloc,nt))
if (.not. allocated(ft_tmp_vy)) allocate(ft_tmp_vy(nrecloc,nt))
if (.not. allocated(ft_tmp_vz)) allocate(ft_tmp_vz(nrecloc,nt))
if (.not. allocated(ft_tmp_vx)) then
allocate(ft_tmp_vx(nrecloc,nt))
endif
if (.not. allocated(ft_tmp_vy)) then
allocate(ft_tmp_vy(nrecloc,nt))
endif
if (.not. allocated(ft_tmp_vz)) then
allocate(ft_tmp_vz(nrecloc,nt))
endif

if (.not. allocated(EN_adjt_vx)) allocate(EN_adjt_vx(nrecloc,nt))
if (.not. allocated(EN_adjt_vy)) allocate(EN_adjt_vy(nrecloc,nt))
if (.not. allocated(EN_adjt_vz)) allocate(EN_adjt_vz(nrecloc,nt))
if (.not. allocated(EN_adjt_vx)) then
allocate(EN_adjt_vx(nrecloc,nt))
endif
if (.not. allocated(EN_adjt_vy)) then
allocate(EN_adjt_vy(nrecloc,nt))
endif
if (.not. allocated(EN_adjt_vz)) then
allocate(EN_adjt_vz(nrecloc,nt))
endif

!*** 1st term (easy...) remeber, imag(analytic_sig) = hilbert transform
srcterm1_vx = denvel_vx * real(danalytic_vx) / (envelc_vx**2 + wlx**2)
Expand Down Expand Up @@ -445,13 +546,25 @@ subroutine get_analytic_signal

integer(kind=si) :: ff, irec, tt

if (.not. allocated(ft_dobs_vx)) allocate(ft_dobs_vx(nrecloc,nt))
if (.not. allocated(ft_dobs_vy)) allocate(ft_dobs_vy(nrecloc,nt))
if (.not. allocated(ft_dobs_vz)) allocate(ft_dobs_vz(nrecloc,nt))
if (.not. allocated(ft_dobs_vx)) then
allocate(ft_dobs_vx(nrecloc,nt))
endif
if (.not. allocated(ft_dobs_vy)) then
allocate(ft_dobs_vy(nrecloc,nt))
endif
if (.not. allocated(ft_dobs_vz)) then
allocate(ft_dobs_vz(nrecloc,nt))
endif

if (.not. allocated(ft_dcal_vx)) allocate(ft_dcal_vx(nrecloc,nt))
if (.not. allocated(ft_dcal_vy)) allocate(ft_dcal_vy(nrecloc,nt))
if (.not. allocated(ft_dcal_vz)) allocate(ft_dcal_vz(nrecloc,nt))
if (.not. allocated(ft_dcal_vx)) then
allocate(ft_dcal_vx(nrecloc,nt))
endif
if (.not. allocated(ft_dcal_vy)) then
allocate(ft_dcal_vy(nrecloc,nt))
endif
if (.not. allocated(ft_dcal_vz)) then
allocate(ft_dcal_vz(nrecloc,nt))
endif

!*** 1. Compute DFT
!* Put to zero
Expand Down Expand Up @@ -488,7 +601,9 @@ subroutine get_analytic_signal
enddo

!*** 2. Remove negative frequencies
if (.not. allocated(hh)) allocate(hh(nt))
if (.not. allocated(hh)) then
allocate(hh(nt))
endif
hh(1) = 1.
hh(floor(nt/2.)+1) = 1.

Expand Down
Loading

0 comments on commit 70954a5

Please sign in to comment.