Skip to content

updates scattering perturbations (w/ white noise simplification) #1851

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
danielpeter merged 2 commits into from
Dec 8, 2025
Merged

updates scattering perturbations (w/ white noise simplification) #1851

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6f36936
updates scattering perturbations (using white noise for correlation f…
danielpeter Dec 8, 2025
cfbe0a2
avoid warnings
danielpeter Dec 8, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
257 changes: 231 additions & 26 deletions src/generate_databases/model_scattering.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,10 @@ module model_scattering_par
integer :: num_target_materials = 0
logical :: USE_ALL_MATERIALS = .false.

! white noise (otherwise von Karman noise distribution)
logical :: USE_WHITE_NOISE = .false.

! von Karman distribution
! perturbation array (regular x/y/z grid)
real(kind=CUSTOM_REAL), dimension(:,:,:),allocatable :: perturbation_grid
integer :: pert_Nx,pert_Ny,pert_Nz
Expand Down Expand Up @@ -260,6 +264,7 @@ subroutine setup_grid()
implicit none
real(kind=CUSTOM_REAL) :: x_min,x_min_all,y_min,y_min_all,z_min,z_min_all, &
x_max,x_max_all,y_max,y_max_all,z_max,z_max_all
real(kind=CUSTOM_REAL) :: x_min_elem,y_min_elem,z_min_elem,x_max_elem,y_max_elem,z_max_elem
real(kind=CUSTOM_REAL) :: dim_x,dim_y,dim_z,dim_max
real(kind=CUSTOM_REAL) :: elemsize_min,elemsize_max,elemsize_min_slice,elemsize_min_norm
integer :: ispec,imaterial_id
Expand All @@ -268,14 +273,47 @@ subroutine setup_grid()
if (num_target_materials == 0) return

! get mesh properties
! Calculation of whole computational domain
x_min = minval(xstore(:))
x_max = maxval(xstore(:))
y_min = minval(ystore(:))
y_max = maxval(ystore(:))
z_min = minval(zstore(:))
z_max = maxval(zstore(:))
if (USE_ALL_MATERIALS) then
! Calculation of whole computational domain
x_min = minval(xstore(:))
x_max = maxval(xstore(:))
y_min = minval(ystore(:))
y_max = maxval(ystore(:))
z_min = minval(zstore(:))
z_max = maxval(zstore(:))
else
! only domains of target materials
x_min = HUGEVAL
y_min = HUGEVAL
z_min = HUGEVAL

x_max = - HUGEVAL
y_max = - HUGEVAL
z_max = - HUGEVAL

do ispec = 1,NSPEC_AB
! material id: index 1 has associated material number
imaterial_id = mat_ext_mesh(1,ispec)

! check only elements for target materials (target id 0 means to apply to all materials)
if (any(target_material_ids(:) == imaterial_id)) then
! determine min/max extent of element
call get_elem_minmax_xyz(x_min_elem,x_max_elem,y_min_elem,y_max_elem,z_min_elem,z_max_elem, &
ispec,NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore)

! domain min/max
x_min = min(x_min,x_min_elem)
y_min = min(y_min,y_min_elem)
z_min = min(z_min,z_min_elem)

x_max = max(x_max,x_max_elem)
y_max = max(y_max,y_max_elem)
z_max = max(z_max,z_max_elem)
endif
enddo
endif

! determine min/max on all processes
x_min_all = HUGEVAL
y_min_all = HUGEVAL
z_min_all = HUGEVAL
Expand All @@ -297,6 +335,10 @@ subroutine setup_grid()
dim_y = y_max_all - y_min_all
dim_z = z_max_all - z_min_all

if (dim_x < 0.0_CUSTOM_REAL) dim_x = 0.0_CUSTOM_REAL
if (dim_y < 0.0_CUSTOM_REAL) dim_y = 0.0_CUSTOM_REAL
if (dim_z < 0.0_CUSTOM_REAL) dim_z = 0.0_CUSTOM_REAL

! maximum mesh length
dim_max = max(dim_x,dim_y,dim_z)

Expand All @@ -321,6 +363,23 @@ subroutine setup_grid()
call flush_IMAIN()
endif

! check if domain volume found, otherwise target material ids won't match any in the mesh
if (dim_x == 0.0_CUSTOM_REAL .or. dim_y == 0.0_CUSTOM_REAL .or. dim_z == 0.0_CUSTOM_REAL) then
! user info
if (myrank == 0) then
write(IMAIN,*) ' no valid domain found (target material ids not found in mesh)'
write(IMAIN,*) ' no perturbations will be applied...'
write(IMAIN,*)
call flush_IMAIN()
endif
! set zero targets
num_target_materials = 0
! all done
return
endif

! determine cell size for FFT grid
!
! get minimum element size
elemsize_min_slice = HUGEVAL
do ispec = 1, NSPEC_AB
Expand Down Expand Up @@ -466,9 +525,6 @@ subroutine generate_perturbations()
integer :: num_seed
integer,dimension(:),allocatable :: myseed

! debug output
character(len=MAX_STRING_LEN) :: name

! external function
real,external :: psd_vonKarman_3D

Expand All @@ -479,6 +535,91 @@ subroutine generate_perturbations()
! checks if anything to do
if (num_target_materials == 0) return

! white noise
! in case correlation length is smaller than GLL point distance, we use white noise
!
! correlation length
! such that k * a ~ 1 (for correlation factor == 1)
! a_corr = real(grid_lambda_min / (2.d0 * PI) * SCATTERING_CORRELATION,kind=CUSTOM_REAL)
!
! therefore
! grid_lambda_min / (2.d0 * PI) * SCATTERING_CORRELATION < elemsize_min_norm / (NGLLX-1)
! using the estimates for grid_lambda_min = ( elemsize_min_norm / (NGLLX-1) * 5), this is equal to
! 5 / (2.d0 * PI) * SCATTERING_CORRELATION < 1
! and
! SCATTERING_CORRELATION < 1 / 5 * 2 * PI ~ 1.25
! let's use white noise for factors < 1
if (SCATTERING_CORRELATION < 1.0d0) then
! set flag
USE_WHITE_NOISE = .true.

! user output
if (myrank == 0) then
write(IMAIN,*) ' perturbation distribution : ','white noise (correlation factor < 1.0)'
write(IMAIN,*) ' perturbation correlation factor : ',sngl(SCATTERING_CORRELATION)
write(IMAIN,*) ' perturbation maximum amplitude : ',sngl(SCATTERING_STRENGTH)
write(IMAIN,*)
call flush_IMAIN()
endif

! no need for FFTs, just random number generation
!
! initializes random number generator
! seed with fixed value to make successive runs repeatable
call random_seed(size=num_seed)
allocate(myseed(num_seed))
myseed(1:num_seed) = 12345
call random_seed(put=myseed)

! visualization output (for debugging)
if (SAVE_MESH_FILES .and. myrank == 0) then
! file output of perturbation grid
!
! note: this perturbation grid is not actually used and allocated for white noise perturbations.
! the file output here is meant for debugging and visualization.

! determines grid dimensions
! next power of 2 for FFT
npower_of_2 = int(log2(grid_length / grid_cell_size)) + 1
N = 2**npower_of_2
! limits minimum/maximum number of points for FFT
!if (N < 1024) N = 1024 ! number of minimum points along x-direction (power of 2 for fft) 2**(10)
!if (N < 512) N = 512 ! number of minimum points along x-direction (power of 2 for fft) 2**(9)
if (N < 256) N = 256 ! number of minimum points along x-direction (power of 2 for fft) 2**(8)
!if (N > 8192) N = 8192 ! number of maximum points along x-direction (power of 2 for fft) 2**(13)
!if (N > 4096) N = 4096 ! number of maximum points along x-direction (power of 2 for fft) 2**(12)
if (N > 2048) N = 2048 ! number of minimum points along x-direction (power of 2 for fft) 2**(11)
! power of 2
npower_of_2 = ceiling(log2(real(N)))
! for mesh interpolations (uses double precision values)
grid_dx = grid_length / (N-1) ! grid space increment, dx == dy == dz
! for perturbation grid dimension
pert_Nx = int(grid_dim_x / grid_normalization_factor / grid_dx) + 1
pert_Ny = int(grid_dim_y / grid_normalization_factor / grid_dx) + 1
pert_Nz = int(grid_dim_z / grid_normalization_factor / grid_dx) + 1
! make sure perturbation grid is not bigger than wavenumber distribution array
if (pert_Nx > N) pert_Nx = N
if (pert_Ny > N) pert_Ny = N
if (pert_Nz > N) pert_Nz = N

! user output
if (myrank == 0) then
write(IMAIN,*) ' grid spacing : dx = ',sngl(grid_dx * grid_normalization_factor),'(m)'
write(IMAIN,*) ' perturbations: grid Nx/Ny/Nz = ',pert_Nx,'/',pert_Ny,'/',pert_Nz
write(IMAIN,*)
call flush_IMAIN()
endif

! output VTU file for visual inspection
call plot_grid_data()
endif

! all done
return
endif

! from here on, we generate noise with von Karman distribution

! user output
if (myrank == 0) then
write(IMAIN,*) ' perturbation distribution : ','von Karman'
Expand Down Expand Up @@ -528,10 +669,10 @@ subroutine generate_perturbations()
if (myrank == 0) then
write(IMAIN,*) ' FFT using number of grid points : N = ',N
write(IMAIN,*) ' FFT using power of 2 : npower_of_2 = ',npower_of_2
write(IMAIN,*) ' grid spacing : dx = ',sngl(grid_dx * grid_normalization_factor),'(m)'
! memory required
mb_size = dble(N) * dble(N) * dble(N) * dble(CUSTOM_REAL) / 1024.d0 / 1024.d0
write(IMAIN,*) ' grid memory size : ',sngl(mb_size),'MB'
write(IMAIN,*) ' FFT using memory size : ',sngl(mb_size),'MB'
write(IMAIN,*) ' grid spacing : dx = ',sngl(grid_dx * grid_normalization_factor),'(m)'
write(IMAIN,*)
call flush_IMAIN()
endif
Expand Down Expand Up @@ -793,8 +934,7 @@ subroutine generate_perturbations()
! visualization output
if (SAVE_MESH_FILES .and. myrank == 0) then
! output VTU file for visual inspection
name = 'perturbation_grid'
call plot_grid_data(perturbation_grid,pert_Nx,pert_Ny,pert_Nz,name)
call plot_grid_data()
endif

end subroutine generate_perturbations
Expand Down Expand Up @@ -853,17 +993,44 @@ end subroutine get_grid_average_max
!-------------------------------------------------------------------------------------------------
!

subroutine plot_grid_data(array,Nx,Ny,Nz,name)
function get_random_perturbation_value() result(pert_val)

use constants, only: myrank,IMAIN,CUSTOM_REAL,MAX_STRING_LEN
use model_scattering_par, only: grid_origin_x,grid_origin_y,grid_origin_z,grid_dx,grid_normalization_factor
use constants, only: CUSTOM_REAL
use shared_parameters, only: SCATTERING_STRENGTH

implicit none
real(kind=CUSTOM_REAL) :: pert_val
real(kind=CUSTOM_REAL) :: rand_val

integer, intent(in) :: Nx,Ny,Nz
real(kind=CUSTOM_REAL),dimension(Nx,Ny,Nz),intent(in) :: array
character(len=MAX_STRING_LEN),intent(in) :: name
! random value between [0,1]
call random_number(rand_val)

! white noise between [-1,1]
pert_val = 2.0_CUSTOM_REAL * rand_val - 1.0_CUSTOM_REAL

! scale to target strength
pert_val = pert_val * real(SCATTERING_STRENGTH,kind=CUSTOM_REAL)

end function get_random_perturbation_value

!
!-------------------------------------------------------------------------------------------------
!

subroutine plot_grid_data()

use constants, only: myrank,IMAIN,CUSTOM_REAL,MAX_STRING_LEN

use model_scattering_par, only: grid_origin_x,grid_origin_y,grid_origin_z, &
grid_dx,grid_normalization_factor,USE_WHITE_NOISE

use model_scattering_par, only: &
array => perturbation_grid, &
Nx => pert_Nx, &
Ny => pert_Ny, &
Nz => pert_Nz

implicit none
! local parameters
integer :: ne,np,ixyz,n1,n2,n3,n4,n5,n6,n7,n8
integer :: i,j,k,ier
Expand All @@ -875,6 +1042,12 @@ subroutine plot_grid_data(array,Nx,Ny,Nz,name)
integer,dimension(:,:),allocatable :: total_dat_con
character(len=MAX_STRING_LEN) :: mesh_file,var_name

! for white noise
real(kind=CUSTOM_REAL),external :: get_random_perturbation_value

! file output
character(len=MAX_STRING_LEN) :: name

! regular grid
np = Nx * Ny * Nz ! total number of points
ne = (Nx-1) * (Ny-1) * (Nz-1) ! total number of elements
Expand Down Expand Up @@ -906,12 +1079,17 @@ subroutine plot_grid_data(array,Nx,Ny,Nz,name)
if (np /= ixyz) stop 'Invalid grid point'

! array value
total_dat(np) = array(i,j,k)
if (USE_WHITE_NOISE) then
! random perturbation value
total_dat(np) = get_random_perturbation_value()
else
total_dat(np) = array(i,j,k)
endif

! grid point position [-L/2,L/2]
total_dat_xyz(1,np) = grid_origin_x + (i-1) * cell_size
total_dat_xyz(2,np) = grid_origin_y + (j-1) * cell_size
total_dat_xyz(3,np) = grid_origin_z + (k-1) * cell_size
total_dat_xyz(1,np) = real(grid_origin_x + (i-1) * cell_size,kind=CUSTOM_REAL)
total_dat_xyz(2,np) = real(grid_origin_y + (j-1) * cell_size,kind=CUSTOM_REAL)
total_dat_xyz(3,np) = real(grid_origin_z + (k-1) * cell_size,kind=CUSTOM_REAL)
enddo
enddo
enddo
Expand Down Expand Up @@ -950,6 +1128,7 @@ subroutine plot_grid_data(array,Nx,Ny,Nz,name)
enddo

! VTU binary format
name = 'perturbation_grid'
mesh_file = 'OUTPUT_FILES/' // trim(name) // '.vtu'
var_name = 'val'
call write_VTU_movie_data_binary(ne,np,total_dat_xyz,total_dat_con,total_dat,mesh_file,var_name)
Expand Down Expand Up @@ -992,20 +1171,33 @@ subroutine model_scattering_add_perturbations(imaterial_id,xmesh,ymesh,zmesh, &
logical, intent(in) :: ANISOTROPY

! local parameters
double precision :: pert_val,scaling
double precision :: pert_val
integer :: ix,iy,iz,Nx,Ny,Nz
! positioning
double precision :: offset_x,offset_y,offset_z
double precision :: spac_x,spac_y,spac_z
double precision :: gamma_interp_x,gamma_interp_y,gamma_interp_z
double precision :: val1,val2,val3,val4,val5,val6,val7,val8

! for white noise
real(kind=CUSTOM_REAL),external :: get_random_perturbation_value

! checks if anything to do
if (num_target_materials == 0) return

! check if we apply perturbation (to target materials only; target id 0 means to apply to all materials)
if (.not. (any(target_material_ids(:) == imaterial_id) .or. USE_ALL_MATERIALS)) return

! white noise
if (USE_WHITE_NOISE) then
! get random perturbation
pert_val = get_random_perturbation_value()
! apply perturbation to model parameters
call apply_perturbation(pert_val)
! all done
return
endif

! determine spacing and cell for linear interpolation
offset_x = (xmesh - grid_origin_x) / grid_normalization_factor ! position offset in mesh (normalized)
offset_y = (ymesh - grid_origin_y) / grid_normalization_factor
Expand Down Expand Up @@ -1098,8 +1290,19 @@ subroutine model_scattering_add_perturbations(imaterial_id,xmesh,ymesh,zmesh, &
! print *,'debug: loc ',sngl(xmesh),sngl(ymesh),sngl(zmesh),'ixyz',ix,iy,iz, &
! 'perturbation ',pert_val,'corners',val1,val2,val3,val4,val5,val6,val7,val8

! apply perturbation to model parameters
call apply_perturbation(pert_val)

contains

subroutine apply_perturbation(pert_val)
implicit none
double precision, intent(in) :: pert_val
! local parameters
real(kind=CUSTOM_REAL) :: scaling

! adds perturbation
scaling = 1.d0 + pert_val
scaling = real(1.d0 + pert_val,kind=CUSTOM_REAL)

! apply perturbation to model
! define new elastic parameters in the model
Expand Down Expand Up @@ -1131,4 +1334,6 @@ subroutine model_scattering_add_perturbations(imaterial_id,xmesh,ymesh,zmesh, &
c66 = scaling * c66
endif

end subroutine apply_perturbation

end subroutine model_scattering_add_perturbations
Loading
Loading