From 60f6511dc5c6f8884320f5c3e07d3c1707b38853 Mon Sep 17 00:00:00 2001 From: daniel peter Date: Thu, 18 Jun 2015 00:18:19 +0200 Subject: [PATCH] updates tomo model to be usable with additional Qp and Qs values and so for gll model to read in also Q values; updates attenuation routines and inserts a criteria from Anderson & Hart (1978) to limit Qs/Qp ratios for bulk attenuation and Olson scaling --- .../create_regions_mesh.f90 | 4 +- src/generate_databases/get_model.f90 | 14 +- src/generate_databases/model_default.f90 | 7 +- src/generate_databases/model_gll.f90 | 68 ++- src/generate_databases/model_gll_adios.F90 | 31 +- src/generate_databases/model_tomography.f90 | 290 ++++++++-- .../read_partition_files.f90 | 2 + src/generate_databases/save_arrays_solver.f90 | 35 +- .../save_arrays_solver_adios.F90 | 15 + src/meshfem3D/check_mesh_quality.f90 | 230 ++++---- src/meshfem3D/create_regions_mesh.f90 | 2 + src/meshfem3D/meshfem3D.f90 | 16 +- src/meshfem3D/save_databases.f90 | 58 +- src/shared/get_attenuation_model.f90 | 502 +++++++++--------- src/shared/rules.mk | 1 + 15 files changed, 826 insertions(+), 449 deletions(-) diff --git a/src/generate_databases/create_regions_mesh.f90 b/src/generate_databases/create_regions_mesh.f90 index 452922a3e..5dc0ac718 100644 --- a/src/generate_databases/create_regions_mesh.f90 +++ b/src/generate_databases/create_regions_mesh.f90 @@ -343,7 +343,8 @@ subroutine create_regions_mesh() deallocate(xixstore,xiystore,xizstore,& etaxstore,etaystore,etazstore,& gammaxstore,gammaystore,gammazstore) - deallocate(jacobianstore,qkappa_attenuation_store,qmu_attenuation_store) + deallocate(jacobianstore) + deallocate(qkappa_attenuation_store,qmu_attenuation_store) deallocate(kappastore,mustore,rho_vp,rho_vs) deallocate(rho_vpI,rho_vpII,rho_vsI) deallocate(rhoarraystore,kappaarraystore,etastore,phistore,tortstore,permstore) @@ -403,6 +404,7 @@ subroutine crm_ext_allocate_arrays(nspec,LOCAL_PATH,myrank, & allocate(xelm(NGNOD),yelm(NGNOD),zelm(NGNOD),stat=ier) if (ier /= 0) stop 'error allocating array xelm etc.' +! attenuation allocate(qkappa_attenuation_store(NGLLX,NGLLY,NGLLZ,nspec),stat=ier) if (ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays') diff --git a/src/generate_databases/get_model.f90 b/src/generate_databases/get_model.f90 index 51b356909..748c13ce9 100644 --- a/src/generate_databases/get_model.f90 +++ b/src/generate_databases/get_model.f90 @@ -433,6 +433,12 @@ subroutine get_model_values(materials_ext_mesh,nmat_ext_mesh, & iflag_aniso = 0 idomain_id = IDOMAIN_ELASTIC + ! attenuation + ! shear attenuation: arbitrary value, see maximum in constants.h + qmu_atten = ATTENUATION_COMP_MAXIMUM + ! bulk attenuation: arbitrary (high) default value + qkappa_atten = 9999.0_CUSTOM_REAL + ! selects chosen velocity model select case (IMODEL) @@ -450,33 +456,27 @@ subroutine get_model_values(materials_ext_mesh,nmat_ext_mesh, & case (IMODEL_1D_PREM ) ! 1D model profile from PREM call model_1D_prem_iso(xmesh,ymesh,zmesh,rho,vp,vs,qmu_atten) - qkappa_atten = 9999. ! undefined in this model case (IMODEL_1D_PREM_PB ) ! 1D model profile from PREM modified by Piero imaterial_PB = abs(imaterial_id) call model_1D_PREM_routine_PB(xmesh,ymesh,zmesh,rho,vp,vs,imaterial_PB) - ! attenuation: arbitrary value, see maximum in constants.h - qmu_atten = ATTENUATION_COMP_MAXIMUM + ! sets acoustic/elastic domain as given in materials properties iundef = - imaterial_id ! iundef must be positive read(undef_mat_prop(6,iundef),*) idomain_id - qkappa_atten = 9999. ! undefined in this model case (IMODEL_1D_CASCADIA ) ! 1D model profile for Cascadia region call model_1D_cascadia(xmesh,ymesh,zmesh,rho,vp,vs,qmu_atten) - qkappa_atten = 9999. ! undefined in this model case (IMODEL_1D_SOCAL ) ! 1D model profile for Southern California call model_1D_socal(xmesh,ymesh,zmesh,rho,vp,vs,qmu_atten) - qkappa_atten = 9999. ! undefined in this model case (IMODEL_SALTON_TROUGH ) ! gets model values from tomography file call model_salton_trough(xmesh,ymesh,zmesh,rho,vp,vs,qmu_atten) - qkappa_atten = 9999. ! undefined in this model case (IMODEL_TOMO ) ! gets model values from tomography file diff --git a/src/generate_databases/model_default.f90 b/src/generate_databases/model_default.f90 index a4d42ec67..bb38591e4 100644 --- a/src/generate_databases/model_default.f90 +++ b/src/generate_databases/model_default.f90 @@ -35,14 +35,15 @@ subroutine model_default(materials_ext_mesh,nmat_ext_mesh, & undef_mat_prop,nundefMat_ext_mesh, & imaterial_id,imaterial_def, & xmesh,ymesh,zmesh, & - rho,vp,vs,iflag_aniso,qkappa_atten,qmu_atten,idomain_id, & + rho,vp,vs, & + iflag_aniso,qkappa_atten,qmu_atten,idomain_id, & rho_s,kappa_s,rho_f,kappa_f,eta_f,kappa_fr,mu_fr, & phi,tort,kxx,kxy,kxz,kyy,kyz,kzz) ! takes model values specified by mesh properties use generate_databases_par, only: IDOMAIN_ACOUSTIC,IDOMAIN_ELASTIC,IDOMAIN_POROELASTIC - use create_regions_mesh_ext_par + use create_regions_mesh_ext_par, only: CUSTOM_REAL,MAX_STRING_LEN implicit none @@ -71,6 +72,8 @@ subroutine model_default(materials_ext_mesh,nmat_ext_mesh, & character(len=MAX_STRING_LEN) :: str_domain integer :: ier + !print *,'model defaults: ',imaterial_id,imaterial_def + ! check if the material is known or unknown if (imaterial_id > 0) then ! gets velocity model as specified by (cubit) mesh files for elastic & acoustic diff --git a/src/generate_databases/model_gll.f90 b/src/generate_databases/model_gll.f90 index f3521cd44..2a49b1227 100644 --- a/src/generate_databases/model_gll.f90 +++ b/src/generate_databases/model_gll.f90 @@ -38,8 +38,12 @@ subroutine model_gll(myrank,nspec,LOCAL_PATH) - use generate_databases_par,only: NGLLX,NGLLY,NGLLZ,FOUR_THIRDS,IMAIN - use create_regions_mesh_ext_par + use generate_databases_par,only: NGLLX,NGLLY,NGLLZ,FOUR_THIRDS,IMAIN,MAX_STRING_LEN, & + ATTENUATION,FULL_ATTENUATION_SOLID + + use create_regions_mesh_ext_par,only: rhostore,kappastore,mustore,rho_vp,rho_vs, & + qkappa_attenuation_store,qmu_attenuation_store + implicit none integer, intent(in) :: myrank,nspec @@ -63,9 +67,12 @@ subroutine model_gll(myrank,nspec,LOCAL_PATH) !!! use lines for vp only ! density - allocate( rho_read(NGLLX,NGLLY,NGLLZ,nspec),stat=ier) + allocate(rho_read(NGLLX,NGLLY,NGLLZ,nspec),stat=ier) if (ier /= 0) stop 'error allocating array rho_read' + ! user output + if (myrank==0) write(IMAIN,*) ' reading in: rho.bin' + filename = prname_lp(1:len_trim(prname_lp))//'rho.bin' open(unit=28,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier) if (ier /= 0) then @@ -77,9 +84,12 @@ subroutine model_gll(myrank,nspec,LOCAL_PATH) close(28) ! vp - allocate( vp_read(NGLLX,NGLLY,NGLLZ,nspec),stat=ier) + allocate(vp_read(NGLLX,NGLLY,NGLLZ,nspec),stat=ier) if (ier /= 0) stop 'error allocating array vp_read' + ! user output + if (myrank==0) write(IMAIN,*) ' reading in: vp.bin' + filename = prname_lp(1:len_trim(prname_lp))//'vp.bin' open(unit=28,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier) if (ier /= 0) then @@ -91,9 +101,12 @@ subroutine model_gll(myrank,nspec,LOCAL_PATH) close(28) ! vs - allocate( vs_read(NGLLX,NGLLY,NGLLZ,nspec),stat=ier) + allocate(vs_read(NGLLX,NGLLY,NGLLZ,nspec),stat=ier) if (ier /= 0) stop 'error allocating array vs_read' + ! user output + if (myrank==0) write(IMAIN,*) ' reading in: vs.bin' + filename = prname_lp(1:len_trim(prname_lp))//'vs.bin' open(unit=28,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier) if (ier /= 0) then @@ -127,11 +140,46 @@ subroutine model_gll(myrank,nspec,LOCAL_PATH) !!! update arrays that will be saved and used in the solver xspecfem3D !!! the following part is neccessary if you uncommented something above - rhostore = rho_read - kappastore = rhostore * ( vp_read * vp_read - FOUR_THIRDS * vs_read * vs_read ) - mustore = rhostore * vs_read * vs_read - rho_vp = rhostore * vp_read - rho_vs = rhostore * vs_read + rhostore(:,:,:,:) = rho_read(:,:,:,:) + kappastore(:,:,:,:) = rhostore(:,:,:,:) * ( vp_read(:,:,:,:) * vp_read(:,:,:,:) & + - FOUR_THIRDS * vs_read(:,:,:,:) * vs_read(:,:,:,:) ) + mustore(:,:,:,:) = rhostore(:,:,:,:) * vs_read(:,:,:,:) * vs_read(:,:,:,:) + + rho_vp(:,:,:,:) = rhostore(:,:,:,:) * vp_read(:,:,:,:) + rho_vs(:,:,:,:) = rhostore(:,:,:,:) * vs_read(:,:,:,:) + + ! gets attenuation arrays from files + if (ATTENUATION) then + ! shear attenuation + ! user output + if (myrank==0) write(IMAIN,*) ' reading in: qmu.bin' + + filename = prname_lp(1:len_trim(prname_lp))//'qmu.bin' + open(unit=28,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier) + if (ier /= 0) then + print *,'Error opening file: ',trim(filename) + stop 'Error reading qmu.bin file' + endif + + read(28) qmu_attenuation_store + close(28) + + ! bulk attenuation + if (FULL_ATTENUATION_SOLID) then + ! user output + if (myrank==0) write(IMAIN,*) ' reading in: qkappa.bin' + + filename = prname_lp(1:len_trim(prname_lp))//'qkappa.bin' + open(unit=28,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier) + if (ier /= 0) then + print *,'error opening file: ',trim(filename) + stop 'error reading qkappa.bin file' + endif + + read(28) qkappa_attenuation_store + close(28) + endif + endif ! free memory deallocate( rho_read,vp_read,vs_read) diff --git a/src/generate_databases/model_gll_adios.F90 b/src/generate_databases/model_gll_adios.F90 index 8f7682c91..69f970be6 100644 --- a/src/generate_databases/model_gll_adios.F90 +++ b/src/generate_databases/model_gll_adios.F90 @@ -38,7 +38,10 @@ subroutine model_gll_adios(myrank,nspec,LOCAL_PATH) use adios_read_mod - use generate_databases_par,only: NGLLX,NGLLY,NGLLZ,FOUR_THIRDS,IMAIN + + use generate_databases_par,only: NGLLX,NGLLY,NGLLZ,FOUR_THIRDS,IMAIN, & + ATTENUATION,FULL_ATTENUATION_SOLID + use create_regions_mesh_ext_par implicit none @@ -94,6 +97,8 @@ subroutine model_gll_adios(myrank,nspec,LOCAL_PATH) start(1) = local_dim_rho * myrank count_ad(1) = NGLLX * NGLLY * NGLLZ * nspec call adios_selection_boundingbox(sel, 1, start, count_ad) + + ! wave speeds call adios_schedule_read(handle, sel, "rho/array", 0, 1, & rho_read, ier) call adios_schedule_read(handle, sel, "vp/array", 0, 1, & @@ -101,11 +106,24 @@ subroutine model_gll_adios(myrank,nspec,LOCAL_PATH) call adios_schedule_read(handle, sel, "vp/array", 0, 1, & vp_read, ier) + ! attenuation + if (ATTENUATION) then + ! shear attenuation + call adios_schedule_read(handle, sel, "qmu/array", 0, 1, & + qmu_attenuation_store, ier) + ! bulk attenuation + if (FULL_ATTENUATION_SOLID) then + call adios_schedule_read(handle, sel, "qkappa/array", 0, 1, & + qkappa_attenuation_store, ier) + endif + endif + !---------------------------------------. ! Perform read and close the adios file | !---------------------------------------' call adios_perform_reads(handle, ier) if (ier /= 0) call stop_all() + call adios_read_close(handle,ier) call adios_read_finalize_method(ADIOS_READ_METHOD_BP, ier) @@ -132,11 +150,12 @@ subroutine model_gll_adios(myrank,nspec,LOCAL_PATH) !!! update arrays that will be saved and used in the solver xspecfem3D !!! the following part is neccessary if you uncommented something above - rhostore = rho_read - kappastore = rhostore * ( vp_read * vp_read - FOUR_THIRDS * vs_read * vs_read ) - mustore = rhostore * vs_read * vs_read - rho_vp = rhostore * vp_read - rho_vs = rhostore * vs_read + rhostore(:,:,:,:) = rho_read(:,:,:,:) + kappastore(:,:,:,:) = rhostore(:,:,:,:) * ( vp_read(:,:,:,:) * vp_read(:,:,:,:) & + - FOUR_THIRDS * vs_read(:,:,:,:) * vs_read(:,:,:,:) ) + mustore(:,:,:,:) = rhostore(:,:,:,:) * vs_read(:,:,:,:) * vs_read(:,:,:,:) + rho_vp(:,:,:,:) = rhostore(:,:,:,:) * vp_read(:,:,:,:) + rho_vs(:,:,:,:) = rhostore(:,:,:,:) * vs_read(:,:,:,:) ! free memory deallocate( rho_read,vp_read,vs_read) diff --git a/src/generate_databases/model_tomography.f90 b/src/generate_databases/model_tomography.f90 index 7ed1a05ca..cff363e01 100644 --- a/src/generate_databases/model_tomography.f90 +++ b/src/generate_databases/model_tomography.f90 @@ -55,6 +55,7 @@ module model_tomography_par ! models parameter records real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: vp_tomography,vs_tomography,rho_tomography,z_tomography + real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: qp_tomography,qs_tomography ! models entries integer, dimension(:), allocatable :: NX,NY,NZ @@ -66,6 +67,9 @@ module model_tomography_par ! process rank integer :: myrank_tomo + ! data format + logical,dimension(:),allocatable :: tomo_has_q_values + end module model_tomography_par ! @@ -138,6 +142,10 @@ subroutine init_tomography_files() character(len=MAX_STRING_LEN) :: filename character(len=MAX_STRING_LEN) :: string_read integer :: nmaterials + ! data format + logical :: has_q_values + integer :: ntokens + logical,dimension(:),allocatable :: materials_with_q ! sets number of materials to loop over nmaterials = nundefMat_ext_mesh @@ -151,6 +159,11 @@ subroutine init_tomography_files() nmaterials = 1 endif + ! data format flag + allocate(materials_with_q(nmaterials),stat=ier) + if (ier /= 0) stop 'Error allocating array materials_with_q' + materials_with_q(:) = .false. + ! loops over number of undefined materials do iundef = 1, nmaterials @@ -192,13 +205,17 @@ subroutine init_tomography_files() call exit_MPI(myrank_tomo,'error reading tomography file') endif + ! header infos + ! format: #origin_x #origin_y #origin_z #end_x #end_y #end_z call tomo_read_next_line(IIN,string_read) read(string_read,*) dummy,dummy,dummy,dummy,dummy,dummy + ! format: #dx #dy #dz call tomo_read_next_line(IIN,string_read) read(string_read,*) dummy,dummy,dummy ! reads in models entries + ! format: #nx #ny #nz call tomo_read_next_line(IIN,string_read) read(string_read,*) temp_x,temp_y,temp_z @@ -206,11 +223,33 @@ subroutine init_tomography_files() nrec = int(temp_x*temp_y*temp_z) nrecord_max = max(nrecord_max,nrec) + ! format: #vp_min #vp_max #vs_min #vs_max #density_min #density_max call tomo_read_next_line(IIN,string_read) read(string_read,*) dummy,dummy,dummy,dummy,dummy,dummy + ! data records ! checks the number of records for points definition - nlines = 0 + ! first record + nlines = 1 + call tomo_read_next_line(IIN,string_read) ! allows to skip comment lines before data section + + ! checks number of entries of first data line + call tomo_get_number_of_tokens(string_read,ntokens) + !print *,'tomography file: number of tokens on first data line: ',ntokens,' line: ',trim(string_read) + if (ntokens /= 6 .and. ntokens /= 8) then + print *,'Error reading tomography file, data line has wrong number of entries: ',trim(string_read) + stop 'Error reading tomography file' + endif + + ! determines data format + if (ntokens == 8) then + has_q_values = .true. + else + has_q_values = .false. + endif + materials_with_q(ifiles_tomo) = has_q_values + + ! counts remaining records do while (ier == 0) read(IIN,*,iostat=ier) if (ier == 0) nlines = nlines + 1 @@ -269,6 +308,23 @@ subroutine init_tomography_files() VP_MAX(NFILES_TOMO),VS_MAX(NFILES_TOMO),RHO_MAX(NFILES_TOMO),stat=ier) if (ier /= 0) call exit_MPI(myrank_tomo,'not enough memory to allocate tomo arrays') + ! q values + allocate(tomo_has_q_values(NFILES_TOMO),stat=ier) + if (ier /= 0) call exit_MPI(myrank_tomo,'not enough memory to allocate tomo q-flag array') + tomo_has_q_values(:) = .false. + ! stores data format flag + do ifiles_tomo = 1,NFILES_TOMO + tomo_has_q_values(ifiles_tomo) = materials_with_q(ifiles_tomo) + enddo + deallocate(materials_with_q) + + ! only allocate q arrays if needed + if (any(tomo_has_q_values)) then + allocate(qp_tomography(NFILES_TOMO,nrecord_max), & + qs_tomography(NFILES_TOMO,nrecord_max),stat=ier) + if (ier /= 0) call exit_MPI(myrank_tomo,'not enough memory to allocate tomo q-value arrays') + endif + end subroutine init_tomography_files ! @@ -289,11 +345,13 @@ subroutine read_model_tomography() ! local parameters real(kind=CUSTOM_REAL) :: x_tomo,y_tomo,z_tomo,vp_tomo,vs_tomo,rho_tomo + real(kind=CUSTOM_REAL) :: qp_tomo,qs_tomo integer :: irecord,ier,iundef,imat character(len=MAX_STRING_LEN*2) :: tomo_filename character(len=MAX_STRING_LEN) :: filename character(len=MAX_STRING_LEN) :: string_read integer :: nmaterials + logical :: has_q_values ! sets number of materials to loop over nmaterials = nundefMat_ext_mesh @@ -375,7 +433,28 @@ subroutine read_model_tomography() ! first record ! format: #x #y #z #vp #vs #density call tomo_read_next_line(IIN,string_read) - read(string_read,*) x_tomo,y_tomo,z_tomo,vp_tomo,vs_tomo,rho_tomo + + ! determines data format + has_q_values = tomo_has_q_values(imat) + + ! user output + if (myrank_tomo == 0) then + if (has_q_values) then + write(IMAIN,*) ' data format: #x #y #z #vp #vs #density #Q_p #Q_s' + else + write(IMAIN,*) ' data format: #x #y #z #vp #vs #density' + endif + call flush_IMAIN() + endif + + ! reads in first data values + if (has_q_values) then + ! format: #x #y #z #vp #vs #density #Q_p #Q_s + read(string_read,*) x_tomo,y_tomo,z_tomo,vp_tomo,vs_tomo,rho_tomo,qp_tomo,qs_tomo + else + ! format: #x #y #z #vp #vs #density + read(string_read,*) x_tomo,y_tomo,z_tomo,vp_tomo,vs_tomo,rho_tomo + endif ! stores record values vp_tomography(imat,1) = vp_tomo @@ -384,22 +463,40 @@ subroutine read_model_tomography() z_tomography(imat,1) = z_tomo ! reads in record sections - do irecord = 2,nrecord(imat) - read(IIN,*) x_tomo,y_tomo,z_tomo,vp_tomo,vs_tomo,rho_tomo - - ! stores record values - vp_tomography(imat,irecord) = vp_tomo - vs_tomography(imat,irecord) = vs_tomo - rho_tomography(imat,irecord) = rho_tomo - z_tomography(imat,irecord) = z_tomo - enddo - + if (has_q_values) then + do irecord = 2,nrecord(imat) + ! format: #x #y #z #vp #vs #density #Q_p #Q_s + read(IIN,*,iostat=ier) x_tomo,y_tomo,z_tomo,vp_tomo,vs_tomo,rho_tomo,qp_tomo,qs_tomo + if (ier /= 0) stop 'Error reading tomo file line format with q values' + + ! stores record values + vp_tomography(imat,irecord) = vp_tomo + vs_tomography(imat,irecord) = vs_tomo + rho_tomography(imat,irecord) = rho_tomo + z_tomography(imat,irecord) = z_tomo + qp_tomography(imat,irecord) = qp_tomo + qs_tomography(imat,irecord) = qs_tomo + enddo + else + do irecord = 2,nrecord(imat) + ! format: #x #y #z #vp #vs #density + read(IIN,*,iostat=ier) x_tomo,y_tomo,z_tomo,vp_tomo,vs_tomo,rho_tomo + if (ier /= 0) stop 'Error reading tomo file line format' + + ! stores record values + vp_tomography(imat,irecord) = vp_tomo + vs_tomography(imat,irecord) = vs_tomo + rho_tomography(imat,irecord) = rho_tomo + z_tomography(imat,irecord) = z_tomo + enddo + endif close(IIN) ! user output if (myrank_tomo == 0) then write(IMAIN,*) ' number of grid points = NX*NY*NZ:',nrecord(imat) write(IMAIN,*) + call flush_IMAIN() endif enddo @@ -429,8 +526,10 @@ subroutine tomo_read_next_line(unit_in,string_read) ! suppress trailing carriage return (ASCII code 13) if any (e.g. if input text file coming from Windows/DOS) if (index(string_read,achar(13)) > 0) string_read = string_read(1:index(string_read,achar(13))-1) - ! exit loop when we find the first line that is not a comment or a white line + ! reads next line if empty if (len_trim(string_read) == 0) cycle + + ! exit loop when we find the first line that is not a comment or a white line if (string_read(1:1) /= '#') exit enddo @@ -450,7 +549,50 @@ end subroutine tomo_read_next_line !------------------------------------------------------------------------------------------------- ! - subroutine model_tomography(xmesh,ymesh,zmesh,rho_model,vp_model,vs_model,qkappa_atten,qmu_atten,imaterial_id, & + subroutine tomo_get_number_of_tokens(string_read,ntokens) + + use constants, only: MAX_STRING_LEN + implicit none + + character(len=MAX_STRING_LEN),intent(in) :: string_read + integer,intent(out) :: ntokens + + ! local parameters + integer :: i + logical :: previous_is_delim + character(len=1), parameter :: delim_space = ' ' + character(len=1), parameter :: delim_tab = achar(9) ! tab delimiter + + ! initializes + ntokens = 0 + + ! checks if anything to do + if (len_trim(string_read) == 0) return + + ! counts tokens + ntokens = 1 + previous_is_delim = .true. + do i = 1, len_trim(string_read) + ! finds next delimiter (space or tabular) + if (string_read(i:i) == delim_space .or. string_read(i:i) == delim_tab) then + if (.not. previous_is_delim) then + ntokens = ntokens + 1 + previous_is_delim = .true. + endif + else + if (previous_is_delim) previous_is_delim = .false. + endif + enddo + + end subroutine tomo_get_number_of_tokens + + +! +!------------------------------------------------------------------------------------------------- +! + + subroutine model_tomography(xmesh,ymesh,zmesh,rho_model,vp_model,vs_model, & + qkappa_atten,qmu_atten,imaterial_id, & has_tomo_value) use generate_databases_par, only: undef_mat_prop,nundefMat_ext_mesh,IMODEL,ATTENUATION_COMP_MAXIMUM @@ -477,10 +619,17 @@ subroutine model_tomography(xmesh,ymesh,zmesh,rho_model,vp_model,vs_model,qkappa double precision :: gamma_interp_z1,gamma_interp_z2,gamma_interp_z3, & gamma_interp_z4,gamma_interp_z5,gamma_interp_z6,gamma_interp_z7,gamma_interp_z8 - real(kind=CUSTOM_REAL) :: vp1,vp2,vp3,vp4,vp5,vp6,vp7,vp8, & - vs1,vs2,vs3,vs4,vs5,vs6,vs7,vs8,rho1,rho2,rho3,rho4,rho5,rho6,rho7,rho8 + real(kind=CUSTOM_REAL) :: vp1,vp2,vp3,vp4,vp5,vp6,vp7,vp8 + real(kind=CUSTOM_REAL) :: vs1,vs2,vs3,vs4,vs5,vs6,vs7,vs8 + real(kind=CUSTOM_REAL) :: rho1,rho2,rho3,rho4,rho5,rho6,rho7,rho8 - real(kind=CUSTOM_REAL), dimension(NFILES_TOMO) :: vp_final,vs_final,rho_final + real(kind=CUSTOM_REAL) :: vp_final,vs_final,rho_final + + ! attenuation + real(kind=CUSTOM_REAL) :: qp1,qp2,qp3,qp4,qp5,qp6,qp7,qp8 + real(kind=CUSTOM_REAL) :: qs1,qs2,qs3,qs4,qs5,qs6,qs7,qs8 + real(kind=CUSTOM_REAL) :: qp_final,qs_final + real(kind=CUSTOM_REAL) :: L_val ! initializes flag has_tomo_value = .false. @@ -650,7 +799,7 @@ subroutine model_tomography(xmesh,ymesh,zmesh,rho_model,vp_model,vs_model,qkappa rho8 = rho_tomography(imat,p7+1) ! use trilinear interpolation in cell to define Vp Vs and rho - vp_final(imat) = & + vp_final = & vp1*(1.-gamma_interp_x)*(1.-gamma_interp_y)*(1.-gamma_interp_z1) + & vp2*gamma_interp_x*(1.-gamma_interp_y)*(1.-gamma_interp_z2) + & vp3*gamma_interp_x*gamma_interp_y*(1.-gamma_interp_z3) + & @@ -660,7 +809,7 @@ subroutine model_tomography(xmesh,ymesh,zmesh,rho_model,vp_model,vs_model,qkappa vp7*gamma_interp_x*gamma_interp_y*gamma_interp_z3 + & vp8*(1.-gamma_interp_x)*gamma_interp_y*gamma_interp_z4 - vs_final(imat) = & + vs_final = & vs1*(1.-gamma_interp_x)*(1.-gamma_interp_y)*(1.-gamma_interp_z1) + & vs2*gamma_interp_x*(1.-gamma_interp_y)*(1.-gamma_interp_z2) + & vs3*gamma_interp_x*gamma_interp_y*(1.-gamma_interp_z3) + & @@ -670,7 +819,7 @@ subroutine model_tomography(xmesh,ymesh,zmesh,rho_model,vp_model,vs_model,qkappa vs7*gamma_interp_x*gamma_interp_y*gamma_interp_z3 + & vs8*(1.-gamma_interp_x)*gamma_interp_y*gamma_interp_z4 - rho_final(imat) = & + rho_final = & rho1*(1.-gamma_interp_x)*(1.-gamma_interp_y)*(1.-gamma_interp_z1) + & rho2*gamma_interp_x*(1.-gamma_interp_y)*(1.-gamma_interp_z2) + & rho3*gamma_interp_x*gamma_interp_y*(1.-gamma_interp_z3) + & @@ -681,25 +830,91 @@ subroutine model_tomography(xmesh,ymesh,zmesh,rho_model,vp_model,vs_model,qkappa rho8*(1.-gamma_interp_x)*gamma_interp_y*gamma_interp_z4 ! impose minimum and maximum velocity and density if needed - if (vp_final(imat) < VP_MIN(imat)) vp_final(imat) = VP_MIN(imat) - if (vp_final(imat) > VP_MAX(imat)) vp_final(imat) = VP_MAX(imat) + if (vp_final < VP_MIN(imat)) vp_final = VP_MIN(imat) + if (vp_final > VP_MAX(imat)) vp_final = VP_MAX(imat) - if (vs_final(imat) < VS_MIN(imat)) vs_final(imat) = VS_MIN(imat) - if (vs_final(imat) > VS_MAX(imat)) vs_final(imat) = VS_MAX(imat) + if (vs_final < VS_MIN(imat)) vs_final = VS_MIN(imat) + if (vs_final > VS_MAX(imat)) vs_final = VS_MAX(imat) - if (rho_final(imat) > RHO_MAX(imat)) rho_final(imat) = RHO_MAX(imat) - if (rho_final(imat) < RHO_MIN(imat)) rho_final(imat) = RHO_MIN(imat) + if (rho_final > RHO_MAX(imat)) rho_final = RHO_MAX(imat) + if (rho_final < RHO_MIN(imat)) rho_final = RHO_MIN(imat) ! model parameters for the associated negative imaterial_id index in materials file - rho_model = rho_final(imat) - vp_model = vp_final(imat) - vs_model = vs_final(imat) - - ! attenuation: arbitrary value, see maximum in constants.h - qmu_atten = ATTENUATION_COMP_MAXIMUM + rho_model = rho_final + vp_model = vp_final + vs_model = vs_final + + ! attenuation + if (tomo_has_q_values(imat)) then + qp1 = qp_tomography(imat,p0+1) + qp2 = qp_tomography(imat,p1+1) + qp3 = qp_tomography(imat,p2+1) + qp4 = qp_tomography(imat,p3+1) + qp5 = qp_tomography(imat,p4+1) + qp6 = qp_tomography(imat,p5+1) + qp7 = qp_tomography(imat,p6+1) + qp8 = qp_tomography(imat,p7+1) + + qs1 = qs_tomography(imat,p0+1) + qs2 = qs_tomography(imat,p1+1) + qs3 = qs_tomography(imat,p2+1) + qs4 = qs_tomography(imat,p3+1) + qs5 = qs_tomography(imat,p4+1) + qs6 = qs_tomography(imat,p5+1) + qs7 = qs_tomography(imat,p6+1) + qs8 = qs_tomography(imat,p7+1) + + ! use trilinear interpolation in cell + qp_final = & + qp1*(1.-gamma_interp_x)*(1.-gamma_interp_y)*(1.-gamma_interp_z1) + & + qp2*gamma_interp_x*(1.-gamma_interp_y)*(1.-gamma_interp_z2) + & + qp3*gamma_interp_x*gamma_interp_y*(1.-gamma_interp_z3) + & + qp4*(1.-gamma_interp_x)*gamma_interp_y*(1.-gamma_interp_z4) + & + qp5*(1.-gamma_interp_x)*(1.-gamma_interp_y)*gamma_interp_z1 + & + qp6*gamma_interp_x*(1.-gamma_interp_y)*gamma_interp_z2 + & + qp7*gamma_interp_x*gamma_interp_y*gamma_interp_z3 + & + qp8*(1.-gamma_interp_x)*gamma_interp_y*gamma_interp_z4 + + qs_final = & + qs1*(1.-gamma_interp_x)*(1.-gamma_interp_y)*(1.-gamma_interp_z1) + & + qs2*gamma_interp_x*(1.-gamma_interp_y)*(1.-gamma_interp_z2) + & + qs3*gamma_interp_x*gamma_interp_y*(1.-gamma_interp_z3) + & + qs4*(1.-gamma_interp_x)*gamma_interp_y*(1.-gamma_interp_z4) + & + qs5*(1.-gamma_interp_x)*(1.-gamma_interp_y)*gamma_interp_z1 + & + qs6*gamma_interp_x*(1.-gamma_interp_y)*gamma_interp_z2 + & + qs7*gamma_interp_x*gamma_interp_y*gamma_interp_z3 + & + qs8*(1.-gamma_interp_x)*gamma_interp_y*gamma_interp_z4 + + ! attenuation zero (means negligible attenuation) + if (qs_final <= 1.e-5) qs_final = ATTENUATION_COMP_MAXIMUM + if (qp_final <= 1.e-5) qp_final = ATTENUATION_COMP_MAXIMUM + + ! Anderson & Hart (1978) conversion between (Qp,Qs) and (Qkappa,Qmu) + ! factor L + L_val = 4.0/3.0 * (vs_model/vp_model)**2 + + ! shear attenuation + qmu_atten = qs_final + + ! converts to bulk attenuation + qkappa_atten = (1.0 - L_val) / (1.0/qp_final - L_val/qs_final) + + ! attenuation zero (means negligible attenuation) + if (qmu_atten <= 1.e-5) qmu_atten = ATTENUATION_COMP_MAXIMUM + if (qkappa_atten <= 1.e-5) qkappa_atten = ATTENUATION_COMP_MAXIMUM + + ! limits Q values + if (qmu_atten < 1.0d0) qmu_atten = 1.0d0 + if (qmu_atten > ATTENUATION_COMP_MAXIMUM) qmu_atten = ATTENUATION_COMP_MAXIMUM + if (qkappa_atten < 1.0d0) qkappa_atten = 1.0d0 + if (qkappa_atten > ATTENUATION_COMP_MAXIMUM) qkappa_atten = ATTENUATION_COMP_MAXIMUM -!! DK DK Q_kappa is not implemented in this model_tomography routine yet, thus set it to dummy value - qkappa_atten = 9999. + else + ! attenuation: arbitrary value, see maximum in constants.h + qmu_atten = ATTENUATION_COMP_MAXIMUM + ! Q_kappa is not implemented in this model_tomography routine yet, thus set it to dummy value + qkappa_atten = ATTENUATION_COMP_MAXIMUM + endif ! value found has_tomo_value = .true. @@ -734,4 +949,11 @@ subroutine deallocate_tomography_files() deallocate(VP_MIN,VS_MIN,RHO_MIN) deallocate(VP_MAX,VS_MAX,RHO_MAX) + ! q values + if (any(tomo_has_q_values)) then + deallocate(qp_tomography) + deallocate(qs_tomography) + endif + deallocate(tomo_has_q_values) + end subroutine deallocate_tomography_files diff --git a/src/generate_databases/read_partition_files.f90 b/src/generate_databases/read_partition_files.f90 index 67cc2ad39..d6bf2547a 100644 --- a/src/generate_databases/read_partition_files.f90 +++ b/src/generate_databases/read_partition_files.f90 @@ -72,6 +72,8 @@ subroutine read_partition_files allocate(materials_ext_mesh(16,nmat_ext_mesh),stat=ier) if (ier /= 0) stop 'Error allocating array materials_ext_mesh' + materials_ext_mesh(:,:) = 0.d0 + do imat = 1, nmat_ext_mesh ! (visco)elastic or acoustic format: ! #(1) rho #(2) vp #(3) vs #(4) Q_mu #(5) anisotropy_flag #(6) material_domain_id #(7) Q_kappa diff --git a/src/generate_databases/save_arrays_solver.f90 b/src/generate_databases/save_arrays_solver.f90 index 41918d100..f0f24879c 100644 --- a/src/generate_databases/save_arrays_solver.f90 +++ b/src/generate_databases/save_arrays_solver.f90 @@ -429,8 +429,7 @@ subroutine save_arrays_solver_files(nspec,nglob,ibool) write(IOUT) v_tmp close(IOUT) - ! VTK file output - ! vp values + ! vp values - VTK file output filename = prname(1:len_trim(prname))//'vp' call write_VTK_data_gll_cr(nspec,nglob, & xstore_dummy,ystore_dummy,zstore_dummy,ibool, & @@ -451,8 +450,7 @@ subroutine save_arrays_solver_files(nspec,nglob,ibool) write(IOUT) v_tmp close(IOUT) - ! VTK file output - ! vs values + ! vs values - VTK file output filename = prname(1:len_trim(prname))//'vs' call write_VTK_data_gll_cr(nspec,nglob, & xstore_dummy,ystore_dummy,zstore_dummy,ibool, & @@ -466,18 +464,37 @@ subroutine save_arrays_solver_files(nspec,nglob,ibool) write(IOUT) v_tmp close(IOUT) - ! VTK file output - ! saves attenuation flag assigned on each gll point into a vtk file - filename = prname(1:len_trim(prname))//'attenuation' + ! attenuation + ! shear attenuation Qmu + open(unit=IOUT,file=prname(1:len_trim(prname))//'qmu.bin',status='unknown',form='unformatted',iostat=ier) + if (ier /= 0) stop 'error opening file qmu.bin' + write(IOUT) qmu_attenuation_store + close(IOUT) + + ! shear attenuation - VTK file output + filename = prname(1:len_trim(prname))//'qmu' call write_VTK_data_gll_cr(nspec,nglob, & xstore_dummy,ystore_dummy,zstore_dummy,ibool, & qmu_attenuation_store,filename) + ! bulk attenuation Qkappa + open(unit=IOUT,file=prname(1:len_trim(prname))//'qkappa.bin',status='unknown',form='unformatted',iostat=ier) + if (ier /= 0) stop 'error opening file qkappa.bin' + write(IOUT) qkappa_attenuation_store + close(IOUT) + + ! bulk attenuation - VTK file output + filename = prname(1:len_trim(prname))//'qkappa' + call write_VTK_data_gll_cr(nspec,nglob, & + xstore_dummy,ystore_dummy,zstore_dummy,ibool, & + qkappa_attenuation_store,filename) + + ! frees temporary array deallocate(v_tmp) - ! VTK file output + ! additional VTK file output if (DEBUG) then - + ! user output call synchronize_all() if (myrank == 0) then write(IMAIN,*) ' saving debugging mesh files' diff --git a/src/generate_databases/save_arrays_solver_adios.F90 b/src/generate_databases/save_arrays_solver_adios.F90 index f94e265b4..fcbde9a85 100644 --- a/src/generate_databases/save_arrays_solver_adios.F90 +++ b/src/generate_databases/save_arrays_solver_adios.F90 @@ -1447,6 +1447,8 @@ subroutine save_arrays_solver_files_adios(nspec,nglob,ibool, nspec_wmax, & call define_adios_scalar(group, groupsize, "", STRINGIFY_VAR(nglob)) local_dim = NGLLX * NGLLY * NGLLZ * nspec_wmax + + ! wave speeds & density call define_adios_global_array1D(group, groupsize, local_dim, & "", "vp", vp_tmp) call define_adios_global_array1D(group, groupsize, local_dim, & @@ -1454,6 +1456,12 @@ subroutine save_arrays_solver_files_adios(nspec,nglob,ibool, nspec_wmax, & call define_adios_global_array1D(group, groupsize, local_dim, & "", "rho", rho_tmp) + ! attenuation + call define_adios_global_array1D(group, groupsize, local_dim, & + "", "qmu", qmu_attenuation_store) + call define_adios_global_array1D(group, groupsize, local_dim, & + "", "qkappa", qkappa_attenuation_store) + !------------------------------------------------------------. ! Open an handler to the ADIOS file and setup the group size | !------------------------------------------------------------' @@ -1472,6 +1480,7 @@ subroutine save_arrays_solver_files_adios(nspec,nglob,ibool, nspec_wmax, & call adios_write(handle, STRINGIFY_VAR(nglob), ier) local_dim = NGLLX * NGLLY * NGLLZ * nspec_wmax + ! wave speeds & density call write_adios_global_1d_array(handle, myrank, sizeprocs, local_dim, & "vp", vp_tmp) call write_adios_global_1d_array(handle, myrank, sizeprocs, local_dim, & @@ -1479,6 +1488,12 @@ subroutine save_arrays_solver_files_adios(nspec,nglob,ibool, nspec_wmax, & call write_adios_global_1d_array(handle, myrank, sizeprocs, local_dim, & "rho", rho_tmp) + ! attenuation + call write_adios_global_1d_array(handle, myrank, sizeprocs, local_dim, & + "qmu", qmu_attenuation_store) + call write_adios_global_1d_array(handle, myrank, sizeprocs, local_dim, & + "qkappa", qkappa_attenuation_store) + !----------------------------------. ! Perform the actual write to disk | !----------------------------------' diff --git a/src/meshfem3D/check_mesh_quality.f90 b/src/meshfem3D/check_mesh_quality.f90 index 7086ff870..58173a760 100644 --- a/src/meshfem3D/check_mesh_quality.f90 +++ b/src/meshfem3D/check_mesh_quality.f90 @@ -96,6 +96,7 @@ subroutine check_mesh_quality(myrank,VP_MAX,NGLOB,NSPEC,x,y,z,ibool, & write(IMAIN,*) '**************************' write(IMAIN,*) write(IMAIN,*) 'start computing the minimum and maximum edge size' + call flush_IMAIN() endif if (NGNOD /= 8) stop 'error: check_mesh_quality only supports NGNOD == 8 for now' @@ -189,71 +190,73 @@ subroutine check_mesh_quality(myrank,VP_MAX,NGLOB,NSPEC,x,y,z,ibool, & if (myrank == 0) then - if (skewness_max_rank /= myrank) then - call recv_i_t(tmp_ispec_max_skewness_MPI,1,skewness_max_rank) - ispec_max_skewness_MPI = tmp_ispec_max_skewness_MPI(1) - else - ispec_max_skewness_MPI = ispec_max_skewness - endif - - write(IMAIN,*) 'done processing ' - - write(IMAIN,*) - write(IMAIN,*) '------------' - write(IMAIN,*) 'mesh quality parameter definitions:' - write(IMAIN,*) - write(IMAIN,*) 'equiangle skewness: 0. perfect, 1. bad' - write(IMAIN,*) 'skewness max deviation angle: 0. perfect, 90. bad' - write(IMAIN,*) 'edge aspect ratio: 1. perfect, above 1. gives stretching factor' - write(IMAIN,*) 'diagonal aspect ratio: 1. perfect, above 1. gives stretching factor' - write(IMAIN,*) '------------' - - write(IMAIN,*) - write(IMAIN,*) 'minimum length of an edge in the whole mesh (m) = ',distance_min_MPI!,' in element ',ispec_min_edge_length - write(IMAIN,*) - write(IMAIN,*) 'maximum length of an edge in the whole mesh (m) = ',distance_max_MPI!,' in element ',ispec_max_edge_length - write(IMAIN,*) - write(IMAIN,*) '***' - write(IMAIN,*) '*** max equiangle skewness = ',equiangle_skewness_max_MPI,' in element ',ispec_max_skewness_MPI, & - ' of slice ',skewness_max_rank - write(IMAIN,*) '***' - write(IMAIN,*) - write(IMAIN,*) 'max deviation angle from a right angle (90 degrees) is therefore = ',90.*equiangle_skewness_max_MPI - write(IMAIN,*) - write(IMAIN,*) 'worst angle in the mesh is therefore ',90.*(1. - equiangle_skewness_max_MPI) - write(IMAIN,*) 'or ',180. - 90.*(1. - equiangle_skewness_max_MPI),' degrees' - write(IMAIN,*) - write(IMAIN,*) 'max edge aspect ratio = ',edge_aspect_ratio_max_MPI - write(IMAIN,*) - write(IMAIN,*) 'max diagonal aspect ratio = ',diagonal_aspect_ratio_max_MPI - write(IMAIN,*) - write(IMAIN,*) '***' - write(IMAIN,'(a50,f13.8)') ' *** Maximum suggested time step for simulation = ',dt_suggested_max_MPI - write(IMAIN,*) '***' - write(IMAIN,*) '*** Max CFL stability condition of the time scheme (must be below about 0.55 or so) = ',stability_max_MPI - write(IMAIN,*) '*** computed using the maximum P wave velocity = ',VP_MAX - write(IMAIN,*) '***' - - ! max stability CFL value -!! DK DK we could probably increase this, now that the Stacey conditions have been fixed - max_CFL_stability_limit = 0.55d0 !! DK DK increased this 0.48d0 - - if (stability_max_MPI >= max_CFL_stability_limit) then - write(IMAIN,*) '*********************************************' - write(IMAIN,*) '*********************************************' - write(IMAIN,*) ' WARNING, that value is above the upper CFL limit of ',max_CFL_stability_limit - write(IMAIN,*) 'therefore the run should be unstable' - write(IMAIN,*) 'You can try to reduce the time step' - write(IMAIN,*) '*********************************************' - write(IMAIN,*) '*********************************************' - else - write(IMAIN,*) 'that value is below the upper CFL limit of ',max_CFL_stability_limit - write(IMAIN,*) 'therefore the run should be stable' - endif - write(IMAIN,*) - - ! create statistics about mesh quality - write(IMAIN,*) 'creating histogram and statistics of mesh quality' + if (skewness_max_rank /= myrank) then + call recv_i_t(tmp_ispec_max_skewness_MPI,1,skewness_max_rank) + ispec_max_skewness_MPI = tmp_ispec_max_skewness_MPI(1) + else + ispec_max_skewness_MPI = ispec_max_skewness + endif + + write(IMAIN,*) 'done processing ' + + write(IMAIN,*) + write(IMAIN,*) '------------' + write(IMAIN,*) 'mesh quality parameter definitions:' + write(IMAIN,*) + write(IMAIN,*) 'equiangle skewness: 0. perfect, 1. bad' + write(IMAIN,*) 'skewness max deviation angle: 0. perfect, 90. bad' + write(IMAIN,*) 'edge aspect ratio: 1. perfect, above 1. gives stretching factor' + write(IMAIN,*) 'diagonal aspect ratio: 1. perfect, above 1. gives stretching factor' + write(IMAIN,*) '------------' + + write(IMAIN,*) + write(IMAIN,*) 'minimum length of an edge in the whole mesh (m) = ',distance_min_MPI!,' in element ',ispec_min_edge_length + write(IMAIN,*) + write(IMAIN,*) 'maximum length of an edge in the whole mesh (m) = ',distance_max_MPI!,' in element ',ispec_max_edge_length + write(IMAIN,*) + write(IMAIN,*) '***' + write(IMAIN,*) '*** max equiangle skewness = ',equiangle_skewness_max_MPI,' in element ',ispec_max_skewness_MPI, & + ' of slice ',skewness_max_rank + write(IMAIN,*) '***' + write(IMAIN,*) + write(IMAIN,*) 'max deviation angle from a right angle (90 degrees) is therefore = ',90.*equiangle_skewness_max_MPI + write(IMAIN,*) + write(IMAIN,*) 'worst angle in the mesh is therefore ',90.*(1. - equiangle_skewness_max_MPI) + write(IMAIN,*) 'or ',180. - 90.*(1. - equiangle_skewness_max_MPI),' degrees' + write(IMAIN,*) + write(IMAIN,*) 'max edge aspect ratio = ',edge_aspect_ratio_max_MPI + write(IMAIN,*) + write(IMAIN,*) 'max diagonal aspect ratio = ',diagonal_aspect_ratio_max_MPI + write(IMAIN,*) + write(IMAIN,*) '***' + write(IMAIN,'(a50,f13.8)') ' *** Maximum suggested time step for simulation = ',dt_suggested_max_MPI + write(IMAIN,*) '***' + write(IMAIN,*) '*** Max CFL stability condition of the time scheme (must be below about 0.55 or so) = ',stability_max_MPI + write(IMAIN,*) '*** computed using the maximum P wave velocity = ',VP_MAX + write(IMAIN,*) '***' + call flush_IMAIN() + + ! max stability CFL value + !! DK DK we could probably increase this, now that the Stacey conditions have been fixed + max_CFL_stability_limit = 0.55d0 !! DK DK increased this 0.48d0 + + if (stability_max_MPI >= max_CFL_stability_limit) then + write(IMAIN,*) '*********************************************' + write(IMAIN,*) '*********************************************' + write(IMAIN,*) ' WARNING, that value is above the upper CFL limit of ',max_CFL_stability_limit + write(IMAIN,*) 'therefore the run should be unstable' + write(IMAIN,*) 'You can try to reduce the time step' + write(IMAIN,*) '*********************************************' + write(IMAIN,*) '*********************************************' + else + write(IMAIN,*) 'that value is below the upper CFL limit of ',max_CFL_stability_limit + write(IMAIN,*) 'therefore the run should be stable' + endif + write(IMAIN,*) + + ! create statistics about mesh quality + write(IMAIN,*) 'creating histogram and statistics of mesh quality' + call flush_IMAIN() endif ! erase histogram of skewness @@ -282,53 +285,56 @@ subroutine check_mesh_quality(myrank,VP_MAX,NGLOB,NSPEC,x,y,z,ibool, & call sum_all_i(NSPEC,NSPEC_ALL_SLICES) if (myrank == 0) then - ! create histogram of skewness and save in Gnuplot file - write(IMAIN,*) - write(IMAIN,*) 'histogram of skewness (0. good - 1. bad):' - write(IMAIN,*) - total_percent = 0. - open(unit=14,file=OUTPUT_FILES(1:len_trim(OUTPUT_FILES))//'mesh_quality_histogram.txt',status='unknown') - do iclass = 0,NCLASS-1 - current_percent = 100.*dble(classes_skewnessMPI(iclass))/dble(NSPEC_ALL_SLICES) - total_percent = total_percent + current_percent - write(IMAIN,*) real(iclass/dble(NCLASS)),' - ',real((iclass+1)/dble(NCLASS)),classes_skewnessMPI(iclass),& - ' ',sngl(current_percent),' %' - write(14,*) 0.5*(real(iclass/dble(NCLASS)) + real((iclass+1)/dble(NCLASS))),' ',sngl(current_percent) - enddo - close(14) - - ! create script for Gnuplot histogram file - open(unit=14,file=OUTPUT_FILES(1:len_trim(OUTPUT_FILES))// & - 'plot_mesh_quality_histogram.gnu',status='unknown') - write(14,*) 'set term wxt' - write(14,*) '#set term gif' - write(14,*) '#set output "mesh_quality_histogram.gif"' - write(14,*) - write(14,*) 'set xrange [0:1]' - write(14,*) 'set xtics 0,0.1,1' - write(14,*) 'set boxwidth ',1./real(NCLASS) - write(14,*) 'set xlabel "Skewness range"' - write(14,*) 'set ylabel "Percentage of elements (%)"' - write(14,*) 'plot "mesh_quality_histogram.txt" with boxes' - write(14,*) 'pause -1 "hit any key..."' - close(14) - - ! display warning if maximum skewness is too high - if (equiangle_skewness_max >= 0.75d0) then - write(IMAIN,*) - write(IMAIN,*) '*********************************************' - write(IMAIN,*) '*********************************************' - write(IMAIN,*) ' WARNING, mesh is bad (max skewness >= 0.75)' - write(IMAIN,*) '*********************************************' - write(IMAIN,*) '*********************************************' - write(IMAIN,*) - endif - - if (total_percent < 99.9d0 .or. total_percent > 100.1d0) then - write(IMAIN,*) 'total percentage = ',total_percent,' %' - stop 'total percentage should be 100%' - endif - + ! create histogram of skewness and save in Gnuplot file + write(IMAIN,*) + write(IMAIN,*) 'histogram of skewness (0. good - 1. bad):' + write(IMAIN,*) + call flush_IMAIN() + + total_percent = 0. + open(unit=14,file=OUTPUT_FILES(1:len_trim(OUTPUT_FILES))//'mesh_quality_histogram.txt',status='unknown') + do iclass = 0,NCLASS-1 + current_percent = 100.*dble(classes_skewnessMPI(iclass))/dble(NSPEC_ALL_SLICES) + total_percent = total_percent + current_percent + write(IMAIN,*) real(iclass/dble(NCLASS)),' - ',real((iclass+1)/dble(NCLASS)),classes_skewnessMPI(iclass),& + ' ',sngl(current_percent),' %' + write(14,*) 0.5*(real(iclass/dble(NCLASS)) + real((iclass+1)/dble(NCLASS))),' ',sngl(current_percent) + enddo + close(14) + + ! create script for Gnuplot histogram file + open(unit=14,file=OUTPUT_FILES(1:len_trim(OUTPUT_FILES))// & + 'plot_mesh_quality_histogram.gnu',status='unknown') + write(14,*) 'set term wxt' + write(14,*) '#set term gif' + write(14,*) '#set output "mesh_quality_histogram.gif"' + write(14,*) + write(14,*) 'set xrange [0:1]' + write(14,*) 'set xtics 0,0.1,1' + write(14,*) 'set boxwidth ',1./real(NCLASS) + write(14,*) 'set xlabel "Skewness range"' + write(14,*) 'set ylabel "Percentage of elements (%)"' + write(14,*) 'plot "mesh_quality_histogram.txt" with boxes' + write(14,*) 'pause -1 "hit any key..."' + close(14) + + ! display warning if maximum skewness is too high + if (equiangle_skewness_max >= 0.75d0) then + write(IMAIN,*) + write(IMAIN,*) '*********************************************' + write(IMAIN,*) '*********************************************' + write(IMAIN,*) ' WARNING, mesh is bad (max skewness >= 0.75)' + write(IMAIN,*) '*********************************************' + write(IMAIN,*) '*********************************************' + write(IMAIN,*) + call flush_IMAIN() + endif + + if (total_percent < 99.9d0 .or. total_percent > 100.1d0) then + write(IMAIN,*) 'total percentage = ',total_percent,' %' + stop 'total percentage should be 100%' + endif + call flush_IMAIN() endif ! debug: for vtk output diff --git a/src/meshfem3D/create_regions_mesh.f90 b/src/meshfem3D/create_regions_mesh.f90 index 6dadef1f1..3728349a3 100644 --- a/src/meshfem3D/create_regions_mesh.f90 +++ b/src/meshfem3D/create_regions_mesh.f90 @@ -391,7 +391,9 @@ subroutine create_regions_mesh(xgrid,ygrid,zgrid,ibool, & write(IMAIN,*) write(IMAIN,*) 'writing out problematic mesh: mesh_debug.vtk' write(IMAIN,*) 'please check your mesh setup...' + call flush_IMAIN() endif + ! debug file output ! vtk file format output open(66,file=prname(1:len_trim(prname))//'mesh_debug.vtk',status='unknown') diff --git a/src/meshfem3D/meshfem3D.f90 b/src/meshfem3D/meshfem3D.f90 index 8db94f955..47dba3528 100644 --- a/src/meshfem3D/meshfem3D.f90 +++ b/src/meshfem3D/meshfem3D.f90 @@ -542,7 +542,9 @@ subroutine meshfem3D if (myrank == 0) then write(IMAIN,*) 'creating global slice addressing' write(IMAIN,*) + call flush_IMAIN() endif + do iproc_eta=0,NPROC_ETA-1 do iproc_xi=0,NPROC_XI-1 iprocnum = iproc_eta * NPROC_XI + iproc_xi @@ -560,6 +562,7 @@ subroutine meshfem3D enddo write(IMAIN,'(a1)',advance='yes') ' ' enddo + call flush_IMAIN() endif if (myrank == 0) then @@ -583,6 +586,7 @@ subroutine meshfem3D write(IMAIN,*) 'Surface shape functions defined by NGNOD2D = ',NGNOD2D_FOUR_CORNERS,' control nodes' write(IMAIN,*) 'Beware! Curvature (i.e. HEX27 elements) is not handled by our internal mesher' write(IMAIN,*) + call flush_IMAIN() endif ! check that the constants.h file is correct @@ -647,6 +651,7 @@ subroutine meshfem3D write(IMAIN,*) 'using UTM projection in region ',UTM_PROJECTION_ZONE endif write(IMAIN,*) + call flush_IMAIN() endif ! get addressing for this process @@ -664,6 +669,7 @@ subroutine meshfem3D write(IMAIN,*) write(IMAIN,*) 'Reading interface data from file ',trim(MF_IN_DATA_FILES)//trim(INTERFACES_FILE) write(IMAIN,*) + call flush_IMAIN() endif open(unit=IIN,file=trim(MF_IN_DATA_FILES)//trim(INTERFACES_FILE),status='old',iostat=ier) @@ -826,6 +832,7 @@ subroutine meshfem3D write(IMAIN,*) 'creating mesh in the model' write(IMAIN,*) '**************************' write(IMAIN,*) + call flush_IMAIN() endif ! assign theoretical number of elements @@ -864,6 +871,7 @@ subroutine meshfem3D if (myrank == 0) then ! compare to exact theoretical value (bottom is always flat) write(IMAIN,*) ' exact area: ',(UTM_Y_MAX-UTM_Y_MIN)*(UTM_X_MAX-UTM_X_MIN) + call flush_IMAIN() endif ! make sure everybody is synchronized @@ -898,6 +906,7 @@ subroutine meshfem3D write(IMAIN,*) write(IMAIN,*) 'smallest and largest possible floating-point numbers are: ',tiny(1._CUSTOM_REAL),huge(1._CUSTOM_REAL) write(IMAIN,*) + call flush_IMAIN() endif ! end of section executed by main process only ! elapsed time since beginning of mesh generation @@ -907,12 +916,11 @@ subroutine meshfem3D write(IMAIN,*) 'Elapsed time for mesh generation and buffer creation in seconds = ',tCPU write(IMAIN,*) 'End of mesh generation' write(IMAIN,*) - endif - -! close main output file - if (myrank == 0) then write(IMAIN,*) 'done' write(IMAIN,*) + call flush_IMAIN() + + ! close main output file close(IMAIN) endif diff --git a/src/meshfem3D/save_databases.f90 b/src/meshfem3D/save_databases.f90 index 9d0cdf0e7..68171f4cb 100644 --- a/src/meshfem3D/save_databases.f90 +++ b/src/meshfem3D/save_databases.f90 @@ -93,6 +93,15 @@ subroutine save_databases(prname,nspec,nglob,iproc_xi,iproc_eta, & integer,dimension(2,nspec) :: material_index character(len=MAX_STRING_LEN), dimension(6,1) :: undef_mat_prop + !------------------------------------------------------------------ + ! user parameter + + ! stores mesh files as cubit for single process run + ! todo: we could put this parameter into the Mesh_Par_file + logical, parameter :: SAVE_MESH_AS_CUBIT = .false. + + !------------------------------------------------------------------ + ! assignes material index ! format: (1,ispec) = #material_id , (2,ispec) = #material_definition material_index (:,:) = 0 @@ -116,6 +125,7 @@ subroutine save_databases(prname,nspec,nglob,iproc_xi,iproc_eta, & ndef = 0 nundef = 0 do i = 1,NMATERIALS + ! material properties format: #rho #vp #vs #Q_flag #anisotropy_flag #domain_id #material_id mat_id = material_properties(i,7) if (mat_id > 0) ndef = ndef + 1 if (mat_id < 0) nundef = nundef + 1 @@ -140,16 +150,22 @@ subroutine save_databases(prname,nspec,nglob,iproc_xi,iproc_eta, & ! materials ! format: #number of defined materials #number of undefined materials write(IIN_database) ndef, nundef + ! writes out defined materials do i = 1,NMATERIALS + ! material properties format: #rho #vp #vs #Q_flag #anisotropy_flag #domain_id #material_id mat_id = material_properties(i,7) if (mat_id > 0) then ! pad dummy zeros to fill up 16 entries (poroelastic medium not allowed) matpropl(:) = 0.d0 + ! material properties format: #rho #vp #vs #Q_flag #anisotropy_flag #domain_id matpropl(1:6) = material_properties(i,1:6) - write(IIN_database) matpropl + ! fills adds arbitrary value for Q_kappa + matpropl(7) = 9999.0 + write(IIN_database) matpropl(:) endif enddo + ! writes out undefined materials do i = 1,NMATERIALS domain_id = material_properties(i,6) @@ -360,11 +376,12 @@ subroutine save_databases(prname,nspec,nglob,iproc_xi,iproc_eta, & write(IIN_database) 0,0 !! VM VM add outputs as cubit - call save_output_mesh_files_as_cubit(nspec,nglob, ibool,nodes_coords, ispec_material_id, & - nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,NSPEC2D_BOTTOM,NSPEC2D_TOP, & - NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NMATERIALS,material_properties, & - ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top) - + if (SAVE_MESH_AS_CUBIT) then + call save_output_mesh_files_as_cubit(nspec,nglob, ibool,nodes_coords, ispec_material_id, & + nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,NSPEC2D_BOTTOM,NSPEC2D_TOP, & + NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NMATERIALS,material_properties, & + ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top) + endif endif close(IIN_database) @@ -387,28 +404,28 @@ subroutine save_output_mesh_files_as_cubit(nspec,nglob, ibool,nodes_coords, ispe integer, parameter :: IIN_database = 15 ! number of spectral elements in each block - integer nspec + integer :: nspec ! number of vertices in each block - integer nglob + integer :: nglob ! MPI cartesian topology ! E for East (= XI_MIN), W for West (= XI_MAX), S for South (= ETA_MIN), N for North (= ETA_MAX) integer, parameter :: W=1,E=2,S=3,N=4,NW=5,NE=6,SE=7,SW=8 ! arrays with the mesh - integer ibool(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) + integer :: ibool(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) double precision :: nodes_coords(nglob,3) - integer ispec_material_id(nspec) + integer :: ispec_material_id(nspec) ! boundary parameters locator - integer NSPEC2D_BOTTOM,NSPEC2D_TOP,NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX - integer nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax - integer ibelm_xmin(NSPEC2DMAX_XMIN_XMAX),ibelm_xmax(NSPEC2DMAX_XMIN_XMAX) - integer ibelm_ymin(NSPEC2DMAX_YMIN_YMAX),ibelm_ymax(NSPEC2DMAX_YMIN_YMAX) - integer ibelm_bottom(NSPEC2D_BOTTOM) - integer ibelm_top(NSPEC2D_TOP) + integer :: NSPEC2D_BOTTOM,NSPEC2D_TOP,NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX + integer :: nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax + integer :: ibelm_xmin(NSPEC2DMAX_XMIN_XMAX),ibelm_xmax(NSPEC2DMAX_XMIN_XMAX) + integer :: ibelm_ymin(NSPEC2DMAX_YMIN_YMAX),ibelm_ymax(NSPEC2DMAX_YMIN_YMAX) + integer :: ibelm_bottom(NSPEC2D_BOTTOM) + integer :: ibelm_top(NSPEC2D_TOP) ! material properties integer :: NMATERIALS @@ -416,13 +433,18 @@ subroutine save_output_mesh_files_as_cubit(nspec,nglob, ibool,nodes_coords, ispe ! second dimension : #rho #vp #vs #Q_flag #anisotropy_flag #domain_id #material_id double precision , dimension(NMATERIALS,7) :: material_properties - integer :: i,ispec,iglob + integer :: i,ispec,iglob,ier ! name of the database files integer :: mat_id,domain_id - open(IIN_database, file = 'MESH/nummaterial_velocity_file') + open(IIN_database, file = 'MESH/nummaterial_velocity_file',status='unknown',action='write',iostat=ier) + if (ier /= 0) then + print *,'Error opening file: ','MESH/nummaterial_velocity_file' + print *,'Please check if directory MESH/ exists for saving mesh files as cubit...' + stop 'Error opening mesh file' + endif do i = 1,NMATERIALS domain_id = material_properties(i,6) diff --git a/src/shared/get_attenuation_model.f90 b/src/shared/get_attenuation_model.f90 index efc73651a..e1aa9eda5 100644 --- a/src/shared/get_attenuation_model.f90 +++ b/src/shared/get_attenuation_model.f90 @@ -25,6 +25,40 @@ ! !===================================================================== + module attenuation_model + + implicit none + + ! model_attenuation_storage_var + type model_attenuation_storage_var + sequence + double precision, dimension(:,:), pointer :: tau_eps_storage + double precision, dimension(:), pointer :: Qmu_storage + integer Q_resolution + integer Q_max + end type model_attenuation_storage_var + type (model_attenuation_storage_var) AM_S + + ! attenuation_simplex_variables + type attenuation_simplex_variables + sequence + double precision Q ! Q = Desired Value of Attenuation or Q + double precision iQ ! iQ = 1/Q + double precision, dimension(:), pointer :: f + ! f = Frequencies at which to evaluate the solution + double precision, dimension(:), pointer :: tau_s + ! tau_s = Tau_sigma defined by the frequency range and + ! number of standard linear solids + integer nf ! nf = Number of Frequencies + integer nsls ! nsls = Number of Standard Linear Solids + end type attenuation_simplex_variables + type(attenuation_simplex_variables) AS_V + + end module attenuation_model + +! +!------------------------------------------------------------------------ +! subroutine get_attenuation_model_olsen(vs_val,Q_mu,OLSEN_ATTENUATION_RATIO) @@ -56,50 +90,27 @@ subroutine get_attenuation_model_olsen(vs_val,Q_mu,OLSEN_ATTENUATION_RATIO) ! v_s in m/s Q_mu = OLSEN_ATTENUATION_RATIO * vs_val - ! uses a simple, 2-constant model mentioned in Olsen et al. (2003) + ! scaling rule if (USE_SIMPLE_OLSEN) then + ! uses a simple, 2-constant model mentioned in Olsen et al. (2003) ! vs (in m/s) if (vs_val < 2000.0_CUSTOM_REAL) then Q_mu = 0.02 * vs_val else Q_mu = 0.1 * vs_val endif - endif - - ! uses discrete values in sediment range - if (USE_DISCRETE_OLSEN) then + else if (USE_DISCRETE_OLSEN) then + ! uses discrete values in sediment range int_Q_mu = 10 * nint(Q_mu / 10.) + ! limits Q to sediment range values if (int_Q_mu < 40) int_Q_mu = 40 if (int_Q_mu > 150) int_Q_mu = 150 - if (int_Q_mu == 40) then - Q_mu = 40.0d0 - else if (int_Q_mu == 50) then - Q_mu = 50.0d0 - else if (int_Q_mu == 60) then - Q_mu = 60.0d0 - else if (int_Q_mu == 70) then - Q_mu = 70.0d0 - else if (int_Q_mu == 80) then - Q_mu = 80.0d0 - else if (int_Q_mu == 90) then - Q_mu = 90.0d0 - else if (int_Q_mu == 100) then - Q_mu = 100.0d0 - else if (int_Q_mu == 110) then - Q_mu = 110.0d0 - else if (int_Q_mu == 120) then - Q_mu = 120.0d0 - else if (int_Q_mu == 130) then - Q_mu = 130.0d0 - else if (int_Q_mu == 140) then - Q_mu = 140.0d0 - else if (int_Q_mu == 150) then - Q_mu = 150.0d0 - else - stop 'incorrect attenuation coefficient' - endif + ! converts to double precision value + Q_mu = dble(int_Q_mu) + else + stop 'Error no Olsen model specified, please set rule in get_attenuation_model.f90' endif ! limits Q_mu value range @@ -123,21 +134,23 @@ subroutine get_attenuation_model(myrank,nspec,USE_OLSEN_ATTENUATION,OLSEN_ATTENU implicit none - double precision :: OLSEN_ATTENUATION_RATIO - integer :: myrank,nspec - logical :: USE_OLSEN_ATTENUATION + double precision,intent(in) :: OLSEN_ATTENUATION_RATIO + integer,intent(in) :: myrank,nspec + logical,intent(in) :: USE_OLSEN_ATTENUATION - real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: mustore - real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: rho_vs - real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: qkappa_attenuation_store - real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: qmu_attenuation_store + real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec),intent(in) :: mustore + real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec),intent(in) :: rho_vs + real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec),intent(in) :: qkappa_attenuation_store + real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec),intent(in) :: qmu_attenuation_store - real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: kappastore - real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: rho_vp + real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec),intent(in) :: kappastore + real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec),intent(in) :: rho_vp - logical, dimension(nspec) :: ispec_is_elastic - real(kind=CUSTOM_REAL) :: min_resolved_period - character(len=MAX_STRING_LEN) :: prname + logical, dimension(nspec),intent(in) :: ispec_is_elastic + real(kind=CUSTOM_REAL),intent(in) :: min_resolved_period + character(len=MAX_STRING_LEN),intent(in) :: prname + + logical,intent(in) :: FULL_ATTENUATION_SOLID ! local parameters real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: one_minus_sum_beta @@ -149,6 +162,7 @@ subroutine get_attenuation_model(myrank,nspec,USE_OLSEN_ATTENUATION,OLSEN_ATTENU double precision factor_scale_dble,one_minus_sum_beta_dble,& factor_scale_dble_kappa,one_minus_sum_beta_dble_kappa double precision :: Q_mu,Q_kappa,Q_p,Q_s + double precision :: L_val double precision :: f_c_source double precision :: MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD real(kind=CUSTOM_REAL), dimension(N_SLS) :: tau_sigma @@ -157,7 +171,15 @@ subroutine get_attenuation_model(myrank,nspec,USE_OLSEN_ATTENUATION,OLSEN_ATTENU real(kind=CUSTOM_REAL):: vs_val,vp_val integer :: i,j,k,ispec,ier double precision :: qmin,qmax,qmin_all,qmax_all - logical :: FULL_ATTENUATION_SOLID + double precision :: qmin_kappa,qmax_kappa,qmin_kappa_all,qmax_kappa_all + + !----------------------------------------------------- + ! user parameter + + ! enforces ratio Qs/Qp >= L factor from Anderson & Hart (1978) + logical, parameter :: USE_ANDERSON_CRITERIA = .true. + + !----------------------------------------------------- ! initializes arrays allocate(one_minus_sum_beta(NGLLX,NGLLY,NGLLZ,nspec), & @@ -189,17 +211,27 @@ subroutine get_attenuation_model(myrank,nspec,USE_OLSEN_ATTENUATION,OLSEN_ATTENU ! precalculates tau_sigma depending on period band (constant for all Q_mu), and ! determines central frequency f_c_source of attenuation period band call get_attenuation_constants(min_resolved_period,tau_sigma_dble,& - f_c_source,MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD) + f_c_source,MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD) ! user output if (myrank == 0) then write(IMAIN,*) write(IMAIN,*) "attenuation: " write(IMAIN,*) " reference period (s) : ",sngl(1.0/ATTENUATION_f0_REFERENCE), & - " frequency: ",sngl(ATTENUATION_f0_REFERENCE) + " frequency: ",sngl(ATTENUATION_f0_REFERENCE) write(IMAIN,*) " period band min/max (s): ",sngl(MIN_ATTENUATION_PERIOD),sngl(MAX_ATTENUATION_PERIOD) write(IMAIN,*) " central period (s) : ",sngl(1.0/f_c_source), & - " frequency: ",sngl(f_c_source) + " frequency: ",sngl(f_c_source) + if (FULL_ATTENUATION_SOLID) then + write(IMAIN,*) " using full attenuation (Q_kappa & Q_mu)" + endif + if (USE_OLSEN_ATTENUATION) then + write(IMAIN,*) " using Olsen scaling with attenuation ratio Qmu/vs = ",sngl(OLSEN_ATTENUATION_RATIO) + if (FULL_ATTENUATION_SOLID .and. USE_ANDERSON_CRITERIA) then + write(IMAIN,*) " using Anderson & Hart criteria for ratio Qs/Qp" + endif + endif + call flush_IMAIN() endif ! determines inverse of tau_sigma @@ -214,6 +246,9 @@ subroutine get_attenuation_model(myrank,nspec,USE_OLSEN_ATTENUATION,OLSEN_ATTENU ! precalculates factors for shear modulus scaling according to attenuation model qmin = HUGEVAL qmax = 0.0 + qmin_kappa = HUGEVAL + qmax_kappa = 0.0 + do ispec = 1,nspec ! skips non elastic elements @@ -224,49 +259,107 @@ subroutine get_attenuation_model(myrank,nspec,USE_OLSEN_ATTENUATION,OLSEN_ATTENU do j=1,NGLLY do i=1,NGLLX + ! initializes Q + Q_mu = 0.d0 + Q_kappa = 0.d0 + ! shear moduli attenuation ! gets Q_mu value if (USE_OLSEN_ATTENUATION) then + ! shear attenuation ! use scaling rule similar to Olsen et al. (2003) vs_val = mustore(i,j,k,ispec) / rho_vs(i,j,k,ispec) call get_attenuation_model_olsen(vs_val,Q_mu,OLSEN_ATTENUATION_RATIO) - if (FULL_ATTENUATION_SOLID) then + else + ! takes Q set in (CUBIT) mesh + Q_mu = qmu_attenuation_store(i,j,k,ispec) + endif + + ! attenuation zero (skips point if zero value) + if (Q_mu <= 1.e-5) cycle + + ! limits Q + if (Q_mu < 1.0d0) Q_mu = 1.0d0 + if (Q_mu > ATTENUATION_COMP_MAXIMUM) Q_mu = ATTENUATION_COMP_MAXIMUM + + ! statistics on Q_mu + if (Q_mu < qmin) qmin = Q_mu + if (Q_mu > qmax) qmax = Q_mu + + ! bulk moduli attenuation + ! gets Q_kappa value + if (FULL_ATTENUATION_SOLID) then + ! gets Q_kappa value + if (USE_OLSEN_ATTENUATION) then + ! bulk attenuation + ! compressional wave speed vp vp_val = (kappastore(i,j,k,ispec) + 2.0d0 * mustore(i,j,k,ispec) / 3.0d0) / rho_vp(i,j,k,ispec) + + ! Anderson & Hart (1978), Q of the Earth, JGR, 83, No. B12 + ! conversion between (Qp,Qs) and (Qkappa,Qmu) + ! factor L + L_val = 4.0d0/3.d0 * (vs_val/vp_val)**2 + + ! attenuation Qs (eq.1) Q_s = Q_mu + + ! scales Qp from Qs (scaling factor introduced by Zhinan?) + ! todo: should we scale Qkappa directly? e.g. Q_kappa = 10.d0 * Q_mu Q_p = 1.5d0 * Q_s - Q_kappa = 1.0d0 / ((1.0/Q_p - 4.0d0/3.0d0*(vp_val/vs_val)**2*(1.d0/Q_mu)) /(1.0d0 - 4.0d0/3.0d0*(vp_val/vs_val)**2)) - if (Q_kappa < 1.0d0) Q_kappa = 1.0d0 - if (Q_kappa > ATTENUATION_COMP_MAXIMUM) Q_kappa = ATTENUATION_COMP_MAXIMUM + + ! Anderson & Hart criteria: Qs/Qp >= L (eq. 4) since Qmu and Qkappa must be positive + if (USE_ANDERSON_CRITERIA) then + ! enforces (eq. 4) from Anderson & Hart + ! note: this might lead to Q_p < Q_s + if ((Q_s - L_val * Q_p) <= 0.d0 ) then + ! negligible bulk attenuation (1/Q_kappa -> zero) + Q_kappa = ATTENUATION_COMP_MAXIMUM + else + ! converts to bulk attenuation (eq. 3) + Q_kappa = (1.0d0 - L_val) * Q_p * Q_s / (Q_s - L_val * Q_p) + endif + else + ! note: this case might lead to: Q_kappa < Q_mu + ! thus having a solid with stronger bulk attenuation than shear attenuation? + + ! avoids division by zero + if (abs(Q_s - L_val * Q_p) <= 1.d-5 ) then + Q_kappa = ATTENUATION_COMP_MAXIMUM + else + ! converts to bulk attenuation (eq. 3) + Q_kappa = (1.0d0 - L_val) * Q_p * Q_s / (Q_s - L_val * Q_p) + endif + endif + + ! debug + !print *,'Qs,Qp:',Q_s,Q_p,'vs,vp:',vs_val,vp_val,'L:',L_val,'crit=',Q_s - L_val * Q_p + !print *,'Qmu,Qkappa:',Q_mu,Q_kappa,'Q_s,Q_p:',Q_s,1./(L_val/Q_mu + (1.0 - L_val)/Q_kappa),'Qp scaled:',Q_p + else + ! takes Q set in (CUBIT) mesh + Q_kappa = qkappa_attenuation_store(i,j,k,ispec) endif - else - ! takes Q set in (CUBIT) mesh - Q_mu = qmu_attenuation_store(i,j,k,ispec) - Q_kappa = qkappa_attenuation_store(i,j,k,ispec) - ! attenuation zero - if (Q_mu <= 1.e-5) cycle - if (Q_kappa <= 1.e-5) cycle + ! attenuation zero (means negligible attenuation) + if (Q_kappa <= 1.e-5) Q_kappa = ATTENUATION_COMP_MAXIMUM ! limits Q - if (Q_mu < 1.0d0) Q_mu = 1.0d0 - if (Q_mu > ATTENUATION_COMP_MAXIMUM) Q_mu = ATTENUATION_COMP_MAXIMUM - if (Q_kappa < 1.0d0) Q_kappa = 1.0d0 if (Q_kappa > ATTENUATION_COMP_MAXIMUM) Q_kappa = ATTENUATION_COMP_MAXIMUM + ! statistics on Q_kappa + if (Q_kappa < qmin_kappa) qmin_kappa = Q_kappa + if (Q_kappa > qmax_kappa) qmax_kappa = Q_kappa endif - ! statistics on Q_mu - if (Q_mu < qmin) qmin = Q_mu - if (Q_mu > qmax) qmax = Q_mu - ! gets beta, on_minus_sum_beta and factor_scale ! based on calculation of strain relaxation times tau_eps call get_attenuation_factors(myrank,Q_mu,MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD, & - f_c_source,tau_sigma_dble, & - beta_dble,one_minus_sum_beta_dble,factor_scale_dble, & - Q_kappa,beta_dble_kappa,one_minus_sum_beta_dble_kappa,factor_scale_dble_kappa,FULL_ATTENUATION_SOLID) + f_c_source,tau_sigma_dble, & + beta_dble,one_minus_sum_beta_dble,factor_scale_dble, & + Q_kappa,beta_dble_kappa,one_minus_sum_beta_dble_kappa,factor_scale_dble_kappa, & + FULL_ATTENUATION_SOLID) + ! shear attenuation ! stores factor for unrelaxed parameter one_minus_sum_beta(i,j,k,ispec) = one_minus_sum_beta_dble @@ -280,6 +373,7 @@ subroutine get_attenuation_model(myrank,nspec,USE_OLSEN_ATTENUATION,OLSEN_ATTENU ! stores scale factor for mu moduli scale_factor(i,j,k,ispec) = factor_scale_dble + ! bulk attenuation if (FULL_ATTENUATION_SOLID) then one_minus_sum_beta_kappa(i,j,k,ispec) = one_minus_sum_beta_dble_kappa beta_kappa(:) = beta_dble_kappa(:) @@ -294,6 +388,23 @@ subroutine get_attenuation_model(myrank,nspec,USE_OLSEN_ATTENUATION,OLSEN_ATTENU enddo enddo + ! statistics + call min_all_dp(qmin,qmin_all) + call max_all_dp(qmax,qmax_all) + if (FULL_ATTENUATION_SOLID) then + call min_all_dp(qmin_kappa,qmin_kappa_all) + call max_all_dp(qmax_kappa,qmax_kappa_all) + endif + + ! user output + if (myrank == 0) then + write(IMAIN,*) " Q_mu min/max : ",sngl(qmin_all),sngl(qmax_all) + if (FULL_ATTENUATION_SOLID) then + write(IMAIN,*) " Q_kappa min/max : ",sngl(qmin_kappa_all),sngl(qmax_kappa_all) + endif + write(IMAIN,*) + endif + ! stores attenuation arrays into files open(unit=27, file=prname(1:len_trim(prname))//'attenuation.bin', & status='unknown',action='write',form='unformatted',iostat=ier) @@ -302,10 +413,11 @@ subroutine get_attenuation_model(myrank,nspec,USE_OLSEN_ATTENUATION,OLSEN_ATTENU call exit_mpi(myrank,'error opening attenuation.bin file') endif write(27) nspec + ! shear attenuation write(27) one_minus_sum_beta write(27) factor_common write(27) scale_factor - + ! bulk attenuation if (FULL_ATTENUATION_SOLID) then write(27) one_minus_sum_beta_kappa write(27) factor_common_kappa @@ -314,21 +426,12 @@ subroutine get_attenuation_model(myrank,nspec,USE_OLSEN_ATTENUATION,OLSEN_ATTENU close(27) + ! frees memory deallocate(one_minus_sum_beta,factor_common,scale_factor) - if (FULL_ATTENUATION_SOLID) then deallocate(one_minus_sum_beta_kappa,factor_common_kappa,scale_factor_kappa) endif - ! statistics - call min_all_dp(qmin,qmin_all) - call max_all_dp(qmax,qmax_all) - ! user output - if (myrank == 0) then - write(IMAIN,*) " Q_mu min/max : ",sngl(qmin_all),sngl(qmax_all) - write(IMAIN,*) - endif - end subroutine get_attenuation_model ! @@ -372,7 +475,7 @@ end subroutine get_attenuation_memory_values ! subroutine get_attenuation_constants(min_resolved_period,tau_sigma, & - f_c_source,MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD) + f_c_source,MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD) ! returns: period band constants tau_sigma and center frequency f_c_source @@ -416,9 +519,10 @@ end subroutine get_attenuation_constants ! subroutine get_attenuation_factors(myrank,Q_mu,MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD, & - f_c_source,tau_sigma, & - beta,one_minus_sum_beta,factor_scale, & - Q_kappa,beta_kappa,one_minus_sum_beta_kappa,factor_scale_kappa,FULL_ATTENUATION_SOLID) + f_c_source,tau_sigma, & + beta,one_minus_sum_beta,factor_scale, & + Q_kappa,beta_kappa,one_minus_sum_beta_kappa,factor_scale_kappa, & + FULL_ATTENUATION_SOLID) ! returns: attenuation mechanisms beta,one_minus_sum_beta,factor_scale @@ -447,7 +551,7 @@ subroutine get_attenuation_factors(myrank,Q_mu,MIN_ATTENUATION_PERIOD,MAX_ATTENU ! determines tau_eps for Q_mu call get_attenuation_tau_eps(Q_mu,tau_sigma,tau_eps, & - MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD) + MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD) ! determines one_minus_sum_beta call get_attenuation_property_values(tau_sigma,tau_eps,beta,one_minus_sum_beta) @@ -458,7 +562,7 @@ subroutine get_attenuation_factors(myrank,Q_mu,MIN_ATTENUATION_PERIOD,MAX_ATTENU if (FULL_ATTENUATION_SOLID) then ! determines tau_eps for Q_kappa call get_attenuation_tau_eps(Q_kappa,tau_sigma,tau_eps_kappa, & - MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD) + MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD) ! determines one_minus_sum_beta call get_attenuation_property_values(tau_sigma,tau_eps_kappa,beta_kappa,one_minus_sum_beta_kappa) @@ -512,7 +616,7 @@ end subroutine get_attenuation_property_values !------------------------------------------------------------------------------------------------- ! - subroutine get_attenuation_scale_factor(myrank, f_c_source, tau_eps, tau_sigma, Q_mu, scale_factor) + subroutine get_attenuation_scale_factor(myrank, f_c_source, tau_eps, tau_sigma, Q_val, scale_factor) ! returns: physical dispersion scaling factor scale_factor @@ -521,7 +625,7 @@ subroutine get_attenuation_scale_factor(myrank, f_c_source, tau_eps, tau_sigma, implicit none integer :: myrank - double precision :: scale_factor, Q_mu, f_c_source + double precision :: scale_factor, Q_val, f_c_source ! strain and stress relaxation times double precision, dimension(N_SLS) :: tau_eps, tau_sigma @@ -544,7 +648,7 @@ subroutine get_attenuation_scale_factor(myrank, f_c_source, tau_eps, tau_sigma, ! Geophys. J. R. Astron. Soc., vol. 47, pp. 41-58 (1976) ! and in Aki, K. and Richards, P. G., Quantitative seismology, theory and methods, ! W. H. Freeman, (1980), second edition, sections 5.5 and 5.5.2, eq. (5.81) p. 170 - factor_scale_mu0 = ONE + TWO * log(f_c_source / ATTENUATION_f0_REFERENCE ) / (PI * Q_mu) + factor_scale_mu0 = ONE + TWO * log(f_c_source / ATTENUATION_f0_REFERENCE ) / (PI * Q_val) !--- compute a, b and Omega parameters ! see e.g.: @@ -577,8 +681,9 @@ subroutine get_attenuation_scale_factor(myrank, f_c_source, tau_eps, tau_sigma, if (scale_factor < 0.7 .or. scale_factor > 1.3) then write(*,*) "error : in get_attenuation_scale_factor() " write(*,*) " scale factor: ", scale_factor, " should be between 0.7 and 1.3" + write(*,*) " Q value = ", Q_val, " central frequency = ",f_c_source write(*,*) " please check your reference frequency ATTENUATION_f0_REFERENCE in constants.h" - call exit_MPI(myrank,'incorrect correction factor in attenuation model') + call exit_MPI(myrank,'unreliable correction factor in attenuation model') endif end subroutine get_attenuation_scale_factor @@ -753,40 +858,17 @@ subroutine get_attenuation_tau_eps(Qmu_in,tau_s,tau_eps,min_period,max_period) ! local parameters integer :: rw - ! model_attenuation_storage_var - type model_attenuation_storage_var - sequence - double precision, dimension(:,:), pointer :: tau_eps_storage - double precision, dimension(:), pointer :: Qmu_storage - integer Q_resolution - integer Q_max - end type model_attenuation_storage_var - type (model_attenuation_storage_var) AM_S - ! attenuation_simplex_variables - type attenuation_simplex_variables - sequence - double precision Q ! Q = Desired Value of Attenuation or Q - double precision iQ ! iQ = 1/Q - double precision, dimension(:), pointer :: f - ! f = Frequencies at which to evaluate the solution - double precision, dimension(:), pointer :: tau_s - ! tau_s = Tau_sigma defined by the frequency range and - ! number of standard linear solids - integer nf ! nf = Number of Frequencies - integer nsls ! nsls = Number of Standard Linear Solids - end type attenuation_simplex_variables - type(attenuation_simplex_variables) AS_V ! READ rw = 1 - call model_attenuation_storage(Qmu_in, tau_eps, rw, AM_S) + call model_attenuation_storage(Qmu_in, tau_eps, rw) if (rw > 0) return - call attenuation_invert_by_simplex(min_period, max_period, N_SLS, Qmu_in, tau_s, tau_eps, AS_V) + call attenuation_invert_by_simplex(min_period, max_period, N_SLS, Qmu_in, tau_s, tau_eps) ! WRITE rw = -1 - call model_attenuation_storage(Qmu_in, tau_eps, rw, AM_S) + call model_attenuation_storage(Qmu_in, tau_eps, rw) end subroutine get_attenuation_tau_eps @@ -795,50 +877,46 @@ end subroutine get_attenuation_tau_eps ! - subroutine model_attenuation_storage(Qmu, tau_eps, rw, AM_S) + subroutine model_attenuation_storage(Qmu, tau_eps, rw) use constants + use attenuation_model,only: AM_S + implicit none -! model_attenuation_storage_var - type model_attenuation_storage_var - sequence - double precision, dimension(:,:), pointer :: tau_eps_storage - double precision, dimension(:), pointer :: Qmu_storage - integer Q_resolution - integer Q_max - end type model_attenuation_storage_var + double precision :: Qmu, Qmu_new + double precision, dimension(N_SLS) :: tau_eps + integer :: rw - type (model_attenuation_storage_var) AM_S -! model_attenuation_storage_var + integer :: Qtmp + integer :: ier - double precision Qmu, Qmu_new - double precision, dimension(N_SLS) :: tau_eps - integer rw + double precision, parameter :: ZERO_TOL = 1.e-5 - integer Qtmp integer, save :: first_time_called = 1 - double precision, parameter :: ZERO_TOL = 1.e-5 - integer ier + ! allocates arrays when first called if (first_time_called == 1) then - first_time_called = 0 - AM_S%Q_resolution = 10**ATTENUATION_COMP_RESOLUTION - AM_S%Q_max = ATTENUATION_COMP_MAXIMUM - Qtmp = AM_S%Q_resolution * AM_S%Q_max + first_time_called = 0 + AM_S%Q_resolution = 10**ATTENUATION_COMP_RESOLUTION + AM_S%Q_max = ATTENUATION_COMP_MAXIMUM + Qtmp = AM_S%Q_resolution * AM_S%Q_max - allocate(AM_S%tau_eps_storage(N_SLS, Qtmp), & + allocate(AM_S%tau_eps_storage(N_SLS, Qtmp), & AM_S%Qmu_storage(Qtmp),stat=ier) - if (ier /= 0) stop 'error allocating arrays for attenuation storage' - AM_S%Qmu_storage(:) = -1 + if (ier /= 0) stop 'error allocating arrays for attenuation storage' + AM_S%Qmu_storage(:) = -1 + + ! debug + !print*,'allocated attenuation arrays: q storage size = ',Qtmp,AM_S%Q_resolution,AM_S%Q_max endif if (Qmu < 0.0d0 .OR. Qmu > AM_S%Q_max) then - write(IMAIN,*) 'Error attenuation_storage()' - write(IMAIN,*) 'Attenuation Value out of Range: ', Qmu - write(IMAIN,*) 'Attenuation Value out of Range: Min, Max ', 0, AM_S%Q_max - call exit_MPI(0, 'Attenuation Value out of Range') + print *,'Error attenuation_storage()' + print *,'Attenuation Value out of Range: ', Qmu + print *,'Attenuation Value out of Range: Min, Max ', 0, AM_S%Q_max + stop 'Attenuation Value out of Range' endif if (rw > 0 .AND. Qmu <= ZERO_TOL) then @@ -855,7 +933,7 @@ subroutine model_attenuation_storage(Qmu, tau_eps, rw, AM_S) ! by default: resolution is Q_resolution = 10 ! converts Qmu to an array integer index: ! e.g. Qmu = 150.31 -> Qtmp = 150.31 * 10 = int( 1503.10 ) = 1503 - Qtmp = int(Qmu * dble(AM_S%Q_resolution)) + Qtmp = int(Qmu * dble(AM_S%Q_resolution)) ! rounds to corresponding double value: ! e.g. Qmu_new = dble( 1503 ) / dble(10) = 150.30 @@ -863,11 +941,19 @@ subroutine model_attenuation_storage(Qmu, tau_eps, rw, AM_S) Qmu_new = dble(Qtmp) / dble(AM_S%Q_resolution) if (rw > 0) then + ! checks + if (first_time_called == 0) then + if (.not. associated(AM_S%Qmu_storage)) & + stop 'error calling model_attenuation_storage() routine without AM_S array' + else + stop 'error calling model_attenuation_storage() routine with first_time_called value invalid' + endif + ! READ if (AM_S%Qmu_storage(Qtmp) > 0) then ! READ SUCCESSFUL - tau_eps(:) = AM_S%tau_eps_storage(:, Qtmp) - Qmu = AM_S%Qmu_storage(Qtmp) + tau_eps(:) = AM_S%tau_eps_storage(:,Qtmp) + Qmu = AM_S%Qmu_storage(Qtmp) rw = 1 else ! READ NOT SUCCESSFUL @@ -875,8 +961,8 @@ subroutine model_attenuation_storage(Qmu, tau_eps, rw, AM_S) endif else ! WRITE SUCCESSFUL - AM_S%tau_eps_storage(:,Qtmp) = tau_eps(:) - AM_S%Qmu_storage(Qtmp) = Qmu + AM_S%tau_eps_storage(:,Qtmp) = tau_eps(:) + AM_S%Qmu_storage(Qtmp) = Qmu rw = 1 endif @@ -886,26 +972,10 @@ end subroutine model_attenuation_storage !------------------------------------------------------------------------------------------------- ! - subroutine attenuation_invert_by_simplex(t2, t1, n, Q_real, tau_s, tau_eps, AS_V) + subroutine attenuation_invert_by_simplex(t2, t1, n, Q_real, tau_s, tau_eps) implicit none - ! attenuation_simplex_variables - type attenuation_simplex_variables - sequence - double precision Q ! Q = Desired Value of Attenuation or Q - double precision iQ ! iQ = 1/Q - double precision, dimension(:), pointer :: f - ! f = Frequencies at which to evaluate the solution - double precision, dimension(:), pointer :: tau_s - ! tau_s = Tau_sigma defined by the frequency range and - ! number of standard linear solids - integer nf ! nf = Number of Frequencies - integer nsls ! nsls = Number of Standard Linear Solids - end type attenuation_simplex_variables - type(attenuation_simplex_variables) AS_V - ! attenuation_simplex_variables - ! Input / Output double precision t1, t2 double precision Q_real @@ -959,7 +1029,7 @@ subroutine attenuation_invert_by_simplex(t2, t1, n, Q_real, tau_s, tau_eps, AS_V ! Shove the paramters into the module - call attenuation_simplex_setup(nf,n,f,Q_real,tau_s,AS_V) + call attenuation_simplex_setup(nf,n,f,Q_real,tau_s) ! Set the Tau_epsilon (tau_eps) to an initial value at omega*tau = 1 ! tan_delta = 1/Q = (tau_eps - tau_s)/(2 * sqrt(tau e*tau_s)) @@ -970,7 +1040,7 @@ subroutine attenuation_invert_by_simplex(t2, t1, n, Q_real, tau_s, tau_eps, AS_V enddo ! Run a simplex search to determine the optimum values of tau_eps - call fminsearch(attenuation_eval, tau_eps, n, iterations, min_value, prnt, err,AS_V) + call fminsearch(attenuation_eval, tau_eps, n, iterations, min_value, prnt, err) if (err > 0) then write(*,*)'Search did not converge for an attenuation of ', Q_real write(*,*)' Iterations: ', iterations @@ -981,7 +1051,7 @@ subroutine attenuation_invert_by_simplex(t2, t1, n, Q_real, tau_s, tau_eps, AS_V !deallocate(f) - call attenuation_simplex_finish(AS_V) + call attenuation_simplex_finish() end subroutine attenuation_invert_by_simplex @@ -990,28 +1060,14 @@ end subroutine attenuation_invert_by_simplex ! - subroutine attenuation_simplex_setup(nf_in,nsls_in,f_in,Q_in,tau_s_in,AS_V) + subroutine attenuation_simplex_setup(nf_in,nsls_in,f_in,Q_in,tau_s_in) ! - Inserts necessary parameters into the module attenuation_simplex_variables ! - See module for explaination - implicit none + use attenuation_model,only: AS_V - ! attenuation_simplex_variables - type attenuation_simplex_variables - sequence - double precision Q ! Q = Desired Value of Attenuation or Q - double precision iQ ! iQ = 1/Q - double precision, dimension(:), pointer :: f - ! f = Frequencies at which to evaluate the solution - double precision, dimension(:), pointer :: tau_s - ! tau_s = Tau_sigma defined by the frequency range and - ! number of standard linear solids - integer nf ! nf = Number of Frequencies - integer nsls ! nsls = Number of Standard Linear Solids - end type attenuation_simplex_variables - type(attenuation_simplex_variables) AS_V - ! attenuation_simplex_variables + implicit none integer nf_in, nsls_in double precision Q_in @@ -1020,7 +1076,7 @@ subroutine attenuation_simplex_setup(nf_in,nsls_in,f_in,Q_in,tau_s_in,AS_V) integer ier allocate(AS_V%f(nf_in), & - AS_V%tau_s(nsls_in),stat=ier) + AS_V%tau_s(nsls_in),stat=ier) if (ier /= 0) stop 'error allocating arrays for attenuation simplex' AS_V%nf = nf_in @@ -1037,7 +1093,7 @@ end subroutine attenuation_simplex_setup ! - double precision function attenuation_eval(Xin,AS_V) + double precision function attenuation_eval(Xin) ! - Computes the misfit from a set of relaxation paramters ! given a set of frequencies and target attenuation @@ -1059,25 +1115,10 @@ double precision function attenuation_eval(Xin,AS_V) ! ! See atteunation_simplex_setup ! + use attenuation_model,only: AS_V implicit none - ! attenuation_simplex_variables - type attenuation_simplex_variables - sequence - double precision Q ! Q = Desired Value of Attenuation or Q - double precision iQ ! iQ = 1/Q - double precision, dimension(:), pointer :: f - ! f = Frequencies at which to evaluate the solution - double precision, dimension(:), pointer :: tau_s - ! tau_s = Tau_sigma defined by the frequency range and - ! number of standard linear solids - integer nf ! nf = Number of Frequencies - integer nsls ! nsls = Number of Standard Linear Solids - end type attenuation_simplex_variables - type(attenuation_simplex_variables) AS_V - ! attenuation_simplex_variables - ! Input double precision, dimension(AS_V%nsls) :: Xin double precision, dimension(AS_V%nsls) :: tau_eps @@ -1171,7 +1212,7 @@ end subroutine attenuation_maxwell ! - subroutine fminsearch(funk, x, n, itercount, tolf, prnt, err, AS_V) + subroutine fminsearch(funk, x, n, itercount, tolf, prnt, err) ! subroutine fminsearch ! - Computes the minimization of funk(x(n)) using the simplex method @@ -1205,23 +1246,6 @@ subroutine fminsearch(funk, x, n, itercount, tolf, prnt, err, AS_V) implicit none -! attenuation_simplex_variables - type attenuation_simplex_variables - sequence - double precision Q ! Q = Desired Value of Attenuation or Q - double precision iQ ! iQ = 1/Q - double precision, dimension(:), pointer :: f - ! f = Frequencies at which to evaluate the solution - double precision, dimension(:), pointer :: tau_s - ! tau_s = Tau_sigma defined by the frequency range and - ! number of standard linear solids - integer nf ! nf = Number of Frequencies - integer nsls ! nsls = Number of Standard Linear Solids - end type attenuation_simplex_variables - - type(attenuation_simplex_variables) AS_V -! attenuation_simplex_variables - ! Input double precision, external :: funk @@ -1283,7 +1307,7 @@ subroutine fminsearch(funk, x, n, itercount, tolf, prnt, err, AS_V) v(:,1) = xin x = xin - fv(1) = funk(xin,AS_V) + fv(1) = funk(xin) usual_delta = 0.05 zero_term_delta = 0.00025 @@ -1297,7 +1321,7 @@ subroutine fminsearch(funk, x, n, itercount, tolf, prnt, err, AS_V) endif v(:,j+1) = y x(:) = y - fv(j+1) = funk(x,AS_V) + fv(j+1) = funk(x) enddo call qsort_local(fv,n+1,place) @@ -1340,13 +1364,13 @@ subroutine fminsearch(funk, x, n, itercount, tolf, prnt, err, AS_V) enddo xr = (1 + rho)*xbar - rho*v(:,n+1) x(:) = xr - fxr = funk(x,AS_V) + fxr = funk(x) func_evals = func_evals + 1 if (fxr < fv(1)) then ! Calculate the expansion point xe = (1 + rho*chi)*xbar - rho*chi*v(:,n+1) x = xe - fxe = funk(x,AS_V) + fxe = funk(x) func_evals = func_evals+1 if (fxe < fxr) then v(:,n+1) = xe @@ -1368,7 +1392,7 @@ subroutine fminsearch(funk, x, n, itercount, tolf, prnt, err, AS_V) ! Perform an outside contraction xc = (1 + psi*rho)*xbar - psi*rho*v(:,n+1) x(:) = xc - fxc = funk(x,AS_V) + fxc = funk(x) func_evals = func_evals+1 if (fxc <= fxr) then @@ -1383,7 +1407,7 @@ subroutine fminsearch(funk, x, n, itercount, tolf, prnt, err, AS_V) ! Perform an inside contraction xcc = (1-psi)*xbar + psi*v(:,n+1) x(:) = xcc - fxcc = funk(x,AS_V) + fxcc = funk(x) func_evals = func_evals+1 if (fxcc < fv(n+1)) then @@ -1399,7 +1423,7 @@ subroutine fminsearch(funk, x, n, itercount, tolf, prnt, err, AS_V) do j=2,n+1 v(:,j)=v(:,1)+sigma*(v(:,j) - v(:,1)) x(:) = v(:,j) - fv(j) = funk(x,AS_V) + fv(j) = funk(x) enddo func_evals = func_evals + n endif @@ -1574,25 +1598,11 @@ end subroutine qsort_local !------------------------------------------------------------------------------------------------- ! - subroutine attenuation_simplex_finish(AS_V) + subroutine attenuation_simplex_finish() - implicit none + use attenuation_model,only: AS_V - ! attenuation_simplex_variables - type attenuation_simplex_variables - sequence - double precision Q ! Q = Desired Value of Attenuation or Q - double precision iQ ! iQ = 1/Q - double precision, dimension(:), pointer :: f - ! f = Frequencies at which to evaluate the solution - double precision, dimension(:), pointer :: tau_s - ! tau_s = Tau_sigma defined by the frequency range and - ! number of standard linear solids - integer nf ! nf = Number of Frequencies - integer nsls ! nsls = Number of Standard Linear Solids - end type attenuation_simplex_variables - type(attenuation_simplex_variables) AS_V - ! attenuation_simplex_variables + implicit none deallocate(AS_V%f) deallocate(AS_V%tau_s) diff --git a/src/shared/rules.mk b/src/shared/rules.mk index 20a5282d4..584c99e3a 100644 --- a/src/shared/rules.mk +++ b/src/shared/rules.mk @@ -78,6 +78,7 @@ shared_MODULES = \ $(FC_MODDIR)/shared_input_parameters.$(FC_MODEXT) \ $(FC_MODDIR)/shared_compute_parameters.$(FC_MODEXT) \ $(FC_MODDIR)/shared_parameters.$(FC_MODEXT) \ + $(FC_MODDIR)/attenuation_model.$(FC_MODEXT) \ $(EMPTY_MACRO)