From e98495b639ab128db61c6e0a6c32311738711828 Mon Sep 17 00:00:00 2001 From: Daniel Peter Date: Wed, 25 Sep 2024 20:15:18 +0200 Subject: [PATCH] updates attenuation info output --- src/specfem2D/attenuation_model.f90 | 56 +++ src/specfem2D/prepare_attenuation.f90 | 574 +++++++++++++++----------- 2 files changed, 385 insertions(+), 245 deletions(-) diff --git a/src/specfem2D/attenuation_model.f90 b/src/specfem2D/attenuation_model.f90 index f0432b437..744559e78 100644 --- a/src/specfem2D/attenuation_model.f90 +++ b/src/specfem2D/attenuation_model.f90 @@ -567,6 +567,10 @@ subroutine decomposition_LU(a,i_min,n,indx,d) end subroutine decomposition_LU +! +!-------------------------------------------------------------------------------- +! + subroutine LUbksb(a,i_min,n,indx,b,m) implicit none @@ -608,6 +612,10 @@ subroutine LUbksb(a,i_min,n,indx,b,m) end subroutine LUbksb +! +!-------------------------------------------------------------------------------- +! + subroutine syst_LU(a,i_min,n,b,m) implicit none @@ -633,6 +641,10 @@ subroutine syst_LU(a,i_min,n,b,m) end subroutine syst_LU +! +!-------------------------------------------------------------------------------- +! + subroutine lfit_zener(x,y,sig,ndat,poids,ia,covar,chisq,ma,Qref,point) ! ma = nombre de variable diffusive ! ndat = m = K nombre d'abcisse freq_k @@ -715,6 +727,10 @@ subroutine lfit_zener(x,y,sig,ndat,poids,ia,covar,chisq,ma,Qref,point) end subroutine lfit_zener +! +!-------------------------------------------------------------------------------- +! + subroutine func_zener(x,afunc,N,Qref,point) implicit none @@ -735,6 +751,10 @@ subroutine func_zener(x,afunc,N,Qref,point) end subroutine func_zener +! +!-------------------------------------------------------------------------------- +! + subroutine remplit_point(fmin,fmax,N,point) use constants, only: TWO_PI @@ -759,6 +779,10 @@ subroutine remplit_point(fmin,fmax,N,point) end subroutine remplit_point +! +!-------------------------------------------------------------------------------- +! + subroutine classical_linear_least_squares(Qref,poids,point,N,fmin,fmax) use constants, only: TWO_PI @@ -798,6 +822,10 @@ subroutine classical_linear_least_squares(Qref,poids,point,N,fmin,fmax) end subroutine classical_linear_least_squares +! +!-------------------------------------------------------------------------------- +! + ! Calcul des coefficients par optimization non-lineaire avec contraintes subroutine solvopt(n,x,f,fun,flg,grad,options,flfc,func,flgc,gradc,Qref,Kopt,theta_min,theta_max,f_min,f_max) @@ -2052,6 +2080,10 @@ subroutine solvopt(n,x,f,fun,flg,grad,options,flfc,func,flgc,gradc,Qref,Kopt,the end subroutine solvopt +! +!-------------------------------------------------------------------------------- +! + subroutine soptions(default) ! SOPTIONS returns the default values for the optional parameters ! used by SolvOpt. @@ -2076,6 +2108,10 @@ subroutine soptions(default) end subroutine soptions +! +!-------------------------------------------------------------------------------- +! + subroutine func_objective(x,res,freq,Qref,N,Nopt) implicit none @@ -2097,6 +2133,10 @@ subroutine func_objective(x,res,freq,Qref,N,Nopt) end subroutine func_objective +! +!-------------------------------------------------------------------------------- +! + subroutine func_mini(x,res,Qref,N,Nopt,K,f_min,f_max) ! Nopt = 2*N : nombre de coefficients a optimiser @@ -2124,6 +2164,10 @@ subroutine func_mini(x,res,Qref,N,Nopt,K,f_min,f_max) end subroutine func_mini +! +!-------------------------------------------------------------------------------- +! + subroutine grad_func_mini(x,grad,Qref,N,Nopt,K,f_min,f_max) use constants, only: TWO_PI @@ -2180,6 +2224,10 @@ subroutine grad_func_mini(x,grad,Qref,N,Nopt,K,f_min,f_max) end subroutine grad_func_mini +! +!-------------------------------------------------------------------------------- +! + subroutine max_residu(x,res,N,Nopt,theta_min,theta_max) implicit none @@ -2203,6 +2251,10 @@ subroutine max_residu(x,res,N,Nopt,theta_min,theta_max) end subroutine max_residu +! +!-------------------------------------------------------------------------------- +! + subroutine grad_max_residu(x,grad,N,Nopt,theta_min,theta_max) implicit none @@ -2249,6 +2301,10 @@ subroutine grad_max_residu(x,grad,N,Nopt,theta_min,theta_max) end subroutine grad_max_residu +! +!-------------------------------------------------------------------------------- +! + subroutine nonlinear_optimization(N,Qref,f0,point,poids,f_min,f_max) use shared_input_parameters, only: USE_SOLVOPT diff --git a/src/specfem2D/prepare_attenuation.f90 b/src/specfem2D/prepare_attenuation.f90 index fafca8414..6837931a8 100644 --- a/src/specfem2D/prepare_attenuation.f90 +++ b/src/specfem2D/prepare_attenuation.f90 @@ -40,13 +40,13 @@ subroutine prepare_attenuation() implicit none ! local parameters - integer :: i,j,ispec,n,ier + integer :: i,j,ispec,ier ! for shifting of velocities if needed in the case of viscoelasticity double precision :: vp,vs,rhol,mul,kappal double precision :: qkappal,qmul - ! attenuation factors + ! temporary attenuation array factors real(kind=CUSTOM_REAL) :: Mu_nu1_sent,Mu_nu2_sent real(kind=CUSTOM_REAL), dimension(:), allocatable :: tau_epsilon_nu1_sent,tau_epsilon_nu2_sent real(kind=CUSTOM_REAL), dimension(:), allocatable :: inv_tau_sigma_nu1_sent,inv_tau_sigma_nu2_sent, & @@ -60,15 +60,321 @@ subroutine prepare_attenuation() ! user output call synchronize_all() if (myrank == 0) then - write(IMAIN,*) write(IMAIN,*) 'Attenuation:' - write(IMAIN,*) ' viscoelastic attenuation:',ATTENUATION_VISCOELASTIC,'(shear & bulk attenuation in elastic domains)' - write(IMAIN,*) ' viscoacoustic attenuation:',ATTENUATION_VISCOACOUSTIC,'(bulk attenuation in acoustic domains)' + write(IMAIN,*) ' viscoelastic attenuation: ',ATTENUATION_VISCOELASTIC, & + '(shear & bulk attenuation in elastic domains)' + write(IMAIN,*) ' viscoacoustic attenuation: ',ATTENUATION_VISCOACOUSTIC, & + '(bulk attenuation in acoustic domains)' + write(IMAIN,*) ' poro-fluid attenuation: ',ATTENUATION_PORO_FLUID_PART, & + '(viscous attenuation in poroelastic domains)' + write(IMAIN,*) ' permittivity attenuation: ',ATTENUATION_PERMITTIVITY, & + '(for electromagnetic materials)' write(IMAIN,*) call flush_IMAIN() endif ! attenuation array allocations + call prepare_attenuation_global_arrays() + + ! checks if anything further to do + if (.not. ATTENUATION_VISCOELASTIC .and. & + .not. ATTENUATION_VISCOACOUSTIC .and. & + .not. ATTENUATION_PORO_FLUID_PART .and. & + .not. ATTENUATION_PERMITTIVITY) return + + ! physical dispersion: scales moduli from reference frequency to simulation (source) center frequency + ! + ! if attenuation is on, shift the velocity model to right frequency; + ! rescale mu to average frequency for attenuation + ! + ! the formulas to implement the scaling can be found for instance in + ! Liu, H. P., Anderson, D. L. and Kanamori, H., Velocity dispersion due to + ! anelasticity: implications for seismology and mantle composition, + ! 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. + ! + ! Beware that in the book of Aki and Richards eq. (5.81) is given for velocities + ! while we need an equation for "mu" and thus we have an additional factor of 2 + ! in the scaling factor below and in equation (49) of Komatitsch and Tromp, Geophys. J. Int. (2002) 149, 390-412, + ! because "mu" is related to the square of velocity. + ! + ! mu(omega_c) = mu(omega_0)[ 1 + 2/(pi Q_mu) ln(omega_c / omega_0) ] + ! + ! if source is not a Dirac or Heavyside then ATTENUATION_f0_REFERENCE is f0 of the first source + if (.not. (time_function_type(1) == 4 .or. time_function_type(1) == 5)) then + ATTENUATION_f0_REFERENCE = f0_source(1) + + if (ATTENUATION_PERMITTIVITY) f0_electromagnetic = f0_source(1) + endif + + ! user output + if (myrank == 0) then + write(IMAIN,*) "Preparing attenuation" + write(IMAIN,*) " The code uses a constant Q quality factor, but approximated" + write(IMAIN,*) " based on a series of Zener standard linear solids (SLS)." + write(IMAIN,*) " Approximation is performed in the following frequency band:" + write(IMAIN,*) + write(IMAIN,*) " number of SLS bodies: ",N_SLS + write(IMAIN,*) + call flush_IMAIN() + endif + + ! setup attenuation + if (ATTENUATION_VISCOELASTIC .or. ATTENUATION_VISCOACOUSTIC) then + ! user output + if (myrank == 0) then + write(IMAIN,*) ' Attenuation in viscoelastic or viscoacoustic parts of the model:' + write(IMAIN,*) " reference frequency (Hz) : ",sngl(ATTENUATION_f0_REFERENCE) + write(IMAIN,*) " period (s) : ",sngl(1.d0/ATTENUATION_f0_REFERENCE) + write(IMAIN,*) + if (READ_VELOCITIES_AT_f0) then + write(IMAIN,*) ' reading velocity at f0 : ',READ_VELOCITIES_AT_f0 + write(IMAIN,*) ' assuming velocity model given at reference frequency' + else + write(IMAIN,*) ' assuming velocity model given at unrelaxed state (infinite frequency)' + endif + write(IMAIN,*) + call flush_IMAIN() + endif + + ! temporary arrays for function argument + allocate(tau_epsilon_nu1_sent(N_SLS), & + tau_epsilon_nu2_sent(N_SLS), & + inv_tau_sigma_nu1_sent(N_SLS), & + inv_tau_sigma_nu2_sent(N_SLS), & + phi_nu1_sent(N_SLS), & + phi_nu2_sent(N_SLS),stat=ier) + if (ier /= 0) call stop_the_code('Error allocating attenuation coefficient arrays') + tau_epsilon_nu1_sent(:) = 0.0_CUSTOM_REAL + tau_epsilon_nu2_sent(:) = 0.0_CUSTOM_REAL + inv_tau_sigma_nu1_sent(:) = 0.0_CUSTOM_REAL + inv_tau_sigma_nu2_sent(:) = 0.0_CUSTOM_REAL + phi_nu1_sent(:) = 0.0_CUSTOM_REAL + phi_nu2_sent(:) = 0.0_CUSTOM_REAL + + ! define the attenuation quality factors. + do ispec = 1,nspec + + ! checks element type for velocity shifts + if (READ_VELOCITIES_AT_f0) then + ! safety check + if (ispec_is_anisotropic(ispec) .or. ispec_is_poroelastic(ispec)) & + call stop_the_code('READ_VELOCITIES_AT_f0 only implemented for non anisotropic, non-poroelastic materials for now') + + if (ispec_is_acoustic(ispec)) then + print * + print *,'******************************' + print *,'WARNING: READ_VELOCITIES_AT_f0 in viscoacoustic elements may imply having to rebuild the mass matrix & + &with the shifted velocities, since the fluid mass matrix contains Kappa; not implemented yet, BEWARE!!' + print * + print *,'viscoacoustic element: ',ispec + print *,'******************************' + print * + endif + endif + + do j = 1,NGLLZ + do i = 1,NGLLX + ! determines relaxation factors + ! bulk attenuation + qkappal = qkappa_attenuation_store(i,j,ispec) + ! shear attenuation + qmul = qmu_attenuation_store(i,j,ispec) + + ! if no attenuation in that elastic element + if (qkappal > 9998.999d0 .and. qmul > 9998.999d0) cycle + + ! determines attenuation factors + call attenuation_model(qkappal,qmul,ATTENUATION_f0_REFERENCE,N_SLS, & + tau_epsilon_nu1_sent,inv_tau_sigma_nu1_sent,phi_nu1_sent,Mu_nu1_sent, & + tau_epsilon_nu2_sent,inv_tau_sigma_nu2_sent,phi_nu2_sent,Mu_nu2_sent) + + ! stores attenuation values + ! bulk attenuation (Qkappa) + inv_tau_sigma_nu1(i,j,ispec,:) = inv_tau_sigma_nu1_sent(:) + tau_epsilon_nu1(i,j,ispec,:) = tau_epsilon_nu1_sent(:) + phi_nu1(i,j,ispec,:) = phi_nu1_sent(:) + + ! shear attenuation (Qmu) + inv_tau_sigma_nu2(i,j,ispec,:) = inv_tau_sigma_nu2_sent(:) + tau_epsilon_nu2(i,j,ispec,:) = tau_epsilon_nu2_sent(:) + phi_nu2(i,j,ispec,:) = phi_nu2_sent(:) + + Mu_nu1(i,j,ispec) = Mu_nu1_sent + Mu_nu2(i,j,ispec) = Mu_nu2_sent + + ! acoustic domains + if (ATTENUATION_VISCOACOUSTIC .and. USE_A_STRONG_FORMULATION_FOR_E1 .and. time_stepping_scheme == 1 ) then + ! bulk attenuation (Qkappa) + phinu(:) = phi_nu1(i,j,ispec,:) + tauinvnu(:) = inv_tau_sigma_nu1(i,j,ispec,:) + temp(:) = exp(- 0.5d0 * tauinvnu(:) * DT) + coef(:) = (1.d0 - temp(:)) / tauinvnu(:) + A_newmark_e1_sf(:,i,j,ispec) = temp(:) + B_newmark_e1_sf(:,i,j,ispec) = phinu(:) * coef(:) + endif + + ! elastic domains + if (ATTENUATION_VISCOELASTIC .and. time_stepping_scheme == 1 ) then + ! bulk attenuation (Qkappa) + phinu(:) = phi_nu1(i,j,ispec,:) + tauinvnu(:) = inv_tau_sigma_nu1(i,j,ispec,:) + temp(:) = exp(- 0.5d0 * tauinvnu(:) * DT) + coef(:) = (1.d0 - temp(:)) / tauinvnu(:) + A_newmark_nu1(:,i,j,ispec) = temp(:) + B_newmark_nu1(:,i,j,ispec) = phinu(:) * coef(:) + + ! shear attenuation (Qmu) + phinu(:) = phi_nu2(i,j,ispec,:) + tauinvnu(:) = inv_tau_sigma_nu2(i,j,ispec,:) + temp(:) = exp(- 0.5d0 * tauinvnu(:) * DT) + coef(:) = (1.d0 - temp(:)) / tauinvnu(:) + A_newmark_nu2(:,i,j,ispec) = temp(:) + B_newmark_nu2(:,i,j,ispec) = phinu(:) * coef(:) + endif + + ! shifts velocities + if (READ_VELOCITIES_AT_f0) then + ! shifts velocity model + rhol = dble(rhostore(i,j,ispec)) + vp = dble(rho_vpstore(i,j,ispec)/rhol) + vs = dble(rho_vsstore(i,j,ispec)/rhol) + + ! shifts vp and vs (according to f0 and attenuation band) + call shift_velocities_from_f0(vp,vs,rhol, & + ATTENUATION_f0_REFERENCE,N_SLS, & + tau_epsilon_nu1_sent,tau_epsilon_nu2_sent, & + inv_tau_sigma_nu1_sent,inv_tau_sigma_nu2_sent) + + ! stores shifted values + ! determines mu and kappa + mul = rhol * vs * vs + if (AXISYM) then ! CHECK kappa + kappal = rhol * vp * vp - FOUR_THIRDS * mul + else + kappal = rhol * vp * vp - mul + endif + ! to compare: + !lambdal = rhol * vp*vp - TWO * mul + !if (AXISYM) then ! CHECK kappa + ! kappal = lambdal + TWO_THIRDS * mul + ! vp = sqrt((kappal + FOUR_THIRDS * mul)/rhol) + !else + ! kappal = lambdal + mul + ! vp = sqrt((kappal + mul)/rhol) + !endif + + ! stores unrelaxed moduli + mustore(i,j,ispec) = mul + kappastore(i,j,ispec) = kappal + + ! stores density times vp and vs + rho_vpstore(i,j,ispec) = rhol * vp + rho_vsstore(i,j,ispec) = rhol * vs + endif + enddo + enddo + enddo + + ! for PMLs + if (PML_BOUNDARY_CONDITIONS) call prepare_attenuation_with_PML() + + endif ! of if (ATTENUATION_VISCOELASTIC .or. ATTENUATION_VISCOACOUSTIC) + + ! allocate memory variables for viscous attenuation (poroelastic media) + if (ATTENUATION_PORO_FLUID_PART) then + ! user output + if (myrank == 0) then + write(IMAIN,*) ' Attenuation in poroelastic parts of the model:' + write(IMAIN,*) ' f0 poroelastic : ',freq0_poroelastic + write(IMAIN,*) ' Q0 poroelastic : ',Q0_poroelastic + write(IMAIN,*) + call flush_IMAIN() + endif + + ! precompute Runge Kutta coefficients if viscous attenuation + ! viscous attenuation is implemented following the memory variable formulation of + ! J. M. Carcione Wave fields in real media: wave propagation in anisotropic, + ! anelastic and porous media, Elsevier, p. 304-305, 2007 + theta_e = (sqrt(Q0_poroelastic**2+1.d0) +1.d0)/(2.d0*pi*freq0_poroelastic*Q0_poroelastic) + theta_s = (sqrt(Q0_poroelastic**2+1.d0) -1.d0)/(2.d0*pi*freq0_poroelastic*Q0_poroelastic) + + thetainv = - 1.d0 / theta_s + alphaval = 1.d0 + DT * thetainv + DT**2 * thetainv**2 / 2.d0 & + + DT**3 * thetainv**3 / 6.d0 + DT**4 * thetainv**4 / 24.d0 + betaval = DT / 2.d0 + DT**2 * thetainv / 3.d0 + DT**3 * thetainv**2 / 8.d0 + DT**4 * thetainv**3 / 24.d0 + gammaval = DT / 2.d0 + DT**2 * thetainv / 6.d0 + DT**3 * thetainv**2 / 24.d0 + endif + + ! allocate memory variables for permittivity attenuation (electromagnetic media) + if (ATTENUATION_PERMITTIVITY) then + ! if source is not a Dirac or Heavyside then f0_attenuation is f0 of the first + ! source + if (.not. (time_function_type(1) == 4 .or. time_function_type(1) == 5)) then + f0_electromagnetic = f0_source(1) + endif + + ! user output + if (myrank == 0) then + write(IMAIN,*) ' Attenuation permittivity in electromagnetic parts of the model:' + write(IMAIN,*) ' f0 electromagnetic : ',f0_electromagnetic + write(IMAIN,*) + call flush_IMAIN() + endif + + ! precompute Runge Kutta coefficients if viscous attenuation + ! viscous attenuation is implemented following the memory variable + ! formulation of + ! J. M. Carcione Wave fields in real media: wave propagation in anisotropic, + ! anelastic and porous media, Elsevier, p. 304-305, 2007 + do ispec = 1,nspec + do j = 1,NGLLZ + do i = 1,NGLLX + + tau_e(i,j,ispec,1) = (sqrt(Qe11_electromagnetic(kmato(ispec))**2+1.d0) +1.d0)/ & + (2.d0*pi*f0_electromagnetic*Qe11_electromagnetic(kmato(ispec))) + tau_e(i,j,ispec,2) = (sqrt(Qe33_electromagnetic(kmato(ispec))**2+1.d0) +1.d0)/ & + (2.d0*pi*f0_electromagnetic*Qe33_electromagnetic(kmato(ispec))) + tau_d(i,j,ispec,1) = (sqrt(Qe11_electromagnetic(kmato(ispec))**2+1.d0) -1.d0)/ & + (2.d0*pi*f0_electromagnetic*Qe11_electromagnetic(kmato(ispec))) + tau_d(i,j,ispec,2) = (sqrt(Qe33_electromagnetic(kmato(ispec))**2+1.d0) -1.d0)/ & + (2.d0*pi*f0_electromagnetic*Qe33_electromagnetic(kmato(ispec))) + + tauinv(:) = - 1.d0 / tau_d(i,j,ispec,:) + alphaval_em(i,j,ispec,:) = 1.d0 + deltat*tauinv(:) + deltat**2*tauinv(:)**2 / 2.d0 & + + deltat**3*tauinv(:)**3 / 6.d0 + deltat**4*tauinv(:)**4 / 24.d0 + betaval_em(i,j,ispec,:) = deltat / 2.d0 + deltat**2*tauinv(:) / 3.d0 + deltat**3*tauinv(:)**2 / 8.d0 & + + deltat**4*tauinv(:)**3 / 24.d0 + gammaval_em(i,j,ispec,:) = deltat / 2.d0 + deltat**2*tauinv(:) / 6.d0 + deltat**3*tauinv(:)**2 / 24.d0 + + enddo + enddo + enddo + endif + + ! synchronizes all processes + call synchronize_all() + + end subroutine prepare_attenuation + +! +!------------------------------------------------------------------------------------- +! + + subroutine prepare_attenuation_global_arrays() + +! allocates all attenuation arrays (required with and without attenuation) + + use constants, only: IMAIN,TWO,PI,FOUR_THIRDS,TWO_THIRDS,USE_A_STRONG_FORMULATION_FOR_E1 + use specfem_par + + implicit none + + ! local parameters + integer :: ier + ! elastic domains if (ATTENUATION_VISCOELASTIC) then nspec_ATT_el = nspec @@ -88,6 +394,13 @@ subroutine prepare_attenuation() A_newmark_nu2(N_SLS,NGLLX,NGLLZ,nspec_ATT_el), & B_newmark_nu2(N_SLS,NGLLX,NGLLZ,nspec_ATT_el), stat=ier) + e1(:,:,:,:) = 0._CUSTOM_REAL + e11(:,:,:,:) = 0._CUSTOM_REAL + e13(:,:,:,:) = 0._CUSTOM_REAL + dux_dxl_old(:,:,:) = 0._CUSTOM_REAL + duz_dzl_old(:,:,:) = 0._CUSTOM_REAL + dux_dzl_plus_duz_dxl_old(:,:,:) = 0._CUSTOM_REAL + ! acoustic domains if (ATTENUATION_VISCOACOUSTIC) then nglob_ATT = nglob @@ -101,7 +414,6 @@ subroutine prepare_attenuation() sum_forces_old(NGLLX,NGLLZ,nspec_ATT_ac), stat=ier) if (ATTENUATION_VISCOACOUSTIC .and. .not. USE_A_STRONG_FORMULATION_FOR_E1) then - allocate(e1_acous(nglob_acoustic,N_SLS), & dot_e1(nglob_acoustic,N_SLS), & A_newmark_e1_sf(1,1,1,1), & @@ -127,7 +439,6 @@ subroutine prepare_attenuation() endif else if (ATTENUATION_VISCOACOUSTIC .and. USE_A_STRONG_FORMULATION_FOR_E1) then - allocate(e1_acous(1,N_SLS), & e1_acous_temp(1,N_SLS), & dot_e1(1,N_SLS), & @@ -151,13 +462,6 @@ subroutine prepare_attenuation() endif if (ier /= 0) call stop_the_code('Error allocating attenuation arrays') - e1(:,:,:,:) = 0._CUSTOM_REAL - e11(:,:,:,:) = 0._CUSTOM_REAL - e13(:,:,:,:) = 0._CUSTOM_REAL - dux_dxl_old(:,:,:) = 0._CUSTOM_REAL - duz_dzl_old(:,:,:) = 0._CUSTOM_REAL - dux_dzl_plus_duz_dxl_old(:,:,:) = 0._CUSTOM_REAL - e1_acous(:,:) = 0._CUSTOM_REAL e1_acous_sf(:,:,:,:) = 0._CUSTOM_REAL @@ -165,6 +469,7 @@ subroutine prepare_attenuation() dot_e1 = 0._CUSTOM_REAL sum_forces_old = 0._CUSTOM_REAL + ! backward/kernel simulations if (SIMULATION_TYPE == 3) then ! acoustic domains if (any_acoustic) then @@ -192,6 +497,7 @@ subroutine prepare_attenuation() endif endif + ! LDDRK time scheme if (time_stepping_scheme == 2) then ! LDDRK ! elastic domains @@ -223,6 +529,7 @@ subroutine prepare_attenuation() e1_LDDRK_acous(:,:) = 0._CUSTOM_REAL + ! RK time scheme if (time_stepping_scheme == 3) then ! RK scheme ! elastic domains @@ -274,15 +581,6 @@ subroutine prepare_attenuation() tau_epsilon_nu2(NGLLX,NGLLZ,max(nspec_ATT_el,nspec_ATT_ac),N_SLS),stat=ier) if (ier /= 0) call stop_the_code('Error allocating attenuation arrays') - ! temporary arrays for function argument - allocate(tau_epsilon_nu1_sent(N_SLS), & - tau_epsilon_nu2_sent(N_SLS), & - inv_tau_sigma_nu1_sent(N_SLS), & - inv_tau_sigma_nu2_sent(N_SLS), & - phi_nu1_sent(N_SLS), & - phi_nu2_sent(N_SLS),stat=ier) - if (ier /= 0) call stop_the_code('Error allocating attenuation coefficient arrays') - ! initialize to dummy values ! convention to indicate that Q = 9999 in that element i.e. that there is no viscoelasticity in that element inv_tau_sigma_nu1(:,:,:,:) = -1._CUSTOM_REAL @@ -299,166 +597,6 @@ subroutine prepare_attenuation() Mu_nu1(:,:,:) = +1._CUSTOM_REAL Mu_nu2(:,:,:) = +1._CUSTOM_REAL - ! physical dispersion: scales moduli from reference frequency to simulation (source) center frequency - ! - ! if attenuation is on, shift the velocity model to right frequency; - ! rescale mu to average frequency for attenuation - ! - ! the formulas to implement the scaling can be found for instance in - ! Liu, H. P., Anderson, D. L. and Kanamori, H., Velocity dispersion due to - ! anelasticity: implications for seismology and mantle composition, - ! 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. - ! - ! Beware that in the book of Aki and Richards eq. (5.81) is given for velocities - ! while we need an equation for "mu" and thus we have an additional factor of 2 - ! in the scaling factor below and in equation (49) of Komatitsch and Tromp, Geophys. J. Int. (2002) 149, 390-412, - ! because "mu" is related to the square of velocity. - ! - ! mu(omega_c) = mu(omega_0)[ 1 + 2/(pi Q_mu) ln(omega_c / omega_0) ] - ! - ! if source is not a Dirac or Heavyside then ATTENUATION_f0_REFERENCE is f0 of the first source - if (.not. (time_function_type(1) == 4 .or. time_function_type(1) == 5)) then - ATTENUATION_f0_REFERENCE = f0_source(1) - endif - - ! setup attenuation - if (ATTENUATION_VISCOELASTIC .or. ATTENUATION_VISCOACOUSTIC) then - ! user output - if (myrank == 0) then - write(IMAIN,*) 'Preparing attenuation in viscoelastic or viscoacoustic parts of the model:' - write(IMAIN,*) ' reading velocity at f0 : ',READ_VELOCITIES_AT_f0 - write(IMAIN,*) - write(IMAIN,*) ' using an attenuation reference frequency of ',ATTENUATION_f0_REFERENCE,'Hz' - write(IMAIN,*) - call flush_IMAIN() - endif - - ! define the attenuation quality factors. - do ispec = 1,nspec - - do j = 1,NGLLZ - do i = 1,NGLLX - ! determines relaxation factors - ! bulk attenuation - qkappal = qkappa_attenuation_store(i,j,ispec) - ! shear attenuation - qmul = qmu_attenuation_store(i,j,ispec) - - ! if no attenuation in that elastic element - if (qkappal > 9998.999d0 .and. qmul > 9998.999d0) cycle - - ! determines attenuation factors - call attenuation_model(qkappal,qmul,ATTENUATION_f0_REFERENCE,N_SLS, & - tau_epsilon_nu1_sent,inv_tau_sigma_nu1_sent,phi_nu1_sent,Mu_nu1_sent, & - tau_epsilon_nu2_sent,inv_tau_sigma_nu2_sent,phi_nu2_sent,Mu_nu2_sent) - - ! stores attenuation values - ! bulk attenuation (Qkappa) - inv_tau_sigma_nu1(i,j,ispec,:) = inv_tau_sigma_nu1_sent(:) - tau_epsilon_nu1(i,j,ispec,:) = tau_epsilon_nu1_sent(:) - phi_nu1(i,j,ispec,:) = phi_nu1_sent(:) - - ! shear attenuation (Qmu) - inv_tau_sigma_nu2(i,j,ispec,:) = inv_tau_sigma_nu2_sent(:) - tau_epsilon_nu2(i,j,ispec,:) = tau_epsilon_nu2_sent(:) - phi_nu2(i,j,ispec,:) = phi_nu2_sent(:) - - Mu_nu1(i,j,ispec) = Mu_nu1_sent - Mu_nu2(i,j,ispec) = Mu_nu2_sent - - ! acoustic domains - if (ATTENUATION_VISCOACOUSTIC .and. USE_A_STRONG_FORMULATION_FOR_E1 .and. time_stepping_scheme == 1 ) then - ! bulk attenuation (Qkappa) - phinu(:) = phi_nu1(i,j,ispec,:) - tauinvnu(:) = inv_tau_sigma_nu1(i,j,ispec,:) - temp(:) = exp(- 0.5d0 * tauinvnu(:) * DT) - coef(:) = (1.d0 - temp(:)) / tauinvnu(:) - A_newmark_e1_sf(:,i,j,ispec) = temp(:) - B_newmark_e1_sf(:,i,j,ispec) = phinu(:) * coef(:) - endif - - ! elastic domains - if (ATTENUATION_VISCOELASTIC .and. time_stepping_scheme == 1 ) then - ! bulk attenuation (Qkappa) - phinu(:) = phi_nu1(i,j,ispec,:) - tauinvnu(:) = inv_tau_sigma_nu1(i,j,ispec,:) - temp(:) = exp(- 0.5d0 * tauinvnu(:) * DT) - coef(:) = (1.d0 - temp(:)) / tauinvnu(:) - A_newmark_nu1(:,i,j,ispec) = temp(:) - B_newmark_nu1(:,i,j,ispec) = phinu(:) * coef(:) - - ! shear attenuation (Qmu) - phinu(:) = phi_nu2(i,j,ispec,:) - tauinvnu(:) = inv_tau_sigma_nu2(i,j,ispec,:) - temp(:) = exp(- 0.5d0 * tauinvnu(:) * DT) - coef(:) = (1.d0 - temp(:)) / tauinvnu(:) - A_newmark_nu2(:,i,j,ispec) = temp(:) - B_newmark_nu2(:,i,j,ispec) = phinu(:) * coef(:) - endif - - ! shifts velocities - if (READ_VELOCITIES_AT_f0) then - - ! safety check - if (ispec_is_anisotropic(ispec) .or. ispec_is_poroelastic(ispec)) & - call stop_the_code('READ_VELOCITIES_AT_f0 only implemented for non anisotropic, non poroelastic materials for now') - - if (ispec_is_acoustic(ispec)) then - do n = 1,100 - print *,'WARNING: READ_VELOCITIES_AT_f0 in viscoacoustic elements may imply having to rebuild the mass matrix & - &with the shifted velocities, since the fluid mass matrix contains Kappa; not implemented yet, BEWARE!!' - enddo - endif - - ! shifts velocity model - rhol = dble(rhostore(i,j,ispec)) - vp = dble(rho_vpstore(i,j,ispec)/rhol) - vs = dble(rho_vsstore(i,j,ispec)/rhol) - - ! shifts vp and vs (according to f0 and attenuation band) - call shift_velocities_from_f0(vp,vs,rhol, & - ATTENUATION_f0_REFERENCE,N_SLS, & - tau_epsilon_nu1_sent,tau_epsilon_nu2_sent, & - inv_tau_sigma_nu1_sent,inv_tau_sigma_nu2_sent) - - ! stores shifted values - ! determines mu and kappa - mul = rhol * vs * vs - if (AXISYM) then ! CHECK kappa - kappal = rhol * vp * vp - FOUR_THIRDS * mul - else - kappal = rhol * vp * vp - mul - endif - ! to compare: - !lambdal = rhol * vp*vp - TWO * mul - !if (AXISYM) then ! CHECK kappa - ! kappal = lambdal + TWO_THIRDS * mul - ! vp = sqrt((kappal + FOUR_THIRDS * mul)/rhol) - !else - ! kappal = lambdal + mul - ! vp = sqrt((kappal + mul)/rhol) - !endif - - ! stores unrelaxed moduli - mustore(i,j,ispec) = mul - kappastore(i,j,ispec) = kappal - - ! stores density times vp and vs - rho_vpstore(i,j,ispec) = rhol * vp - rho_vsstore(i,j,ispec) = rhol * vs - endif - enddo - enddo - enddo - - ! for PMLs - if (PML_BOUNDARY_CONDITIONS) call prepare_attenuation_with_PML() - - endif ! of if (ATTENUATION_VISCOELASTIC .or. ATTENUATION_VISCOACOUSTIC) - ! allocate memory variables for viscous attenuation (poroelastic media) if (ATTENUATION_PORO_FLUID_PART) then allocate(rx_viscous(NGLLX,NGLLZ,nspec)) @@ -488,19 +626,6 @@ subroutine prepare_attenuation() rx_viscous_force_RK(:,:,:,:) = 0.d0 rz_viscous_force_RK(:,:,:,:) = 0.d0 endif - - ! precompute Runge Kutta coefficients if viscous attenuation - ! viscous attenuation is implemented following the memory variable formulation of - ! J. M. Carcione Wave fields in real media: wave propagation in anisotropic, - ! anelastic and porous media, Elsevier, p. 304-305, 2007 - theta_e = (sqrt(Q0_poroelastic**2+1.d0) +1.d0)/(2.d0*pi*freq0_poroelastic*Q0_poroelastic) - theta_s = (sqrt(Q0_poroelastic**2+1.d0) -1.d0)/(2.d0*pi*freq0_poroelastic*Q0_poroelastic) - - thetainv = - 1.d0 / theta_s - alphaval = 1.d0 + DT * thetainv + DT**2 * thetainv**2 / 2.d0 & - + DT**3 * thetainv**3 / 6.d0 + DT**4 * thetainv**4 / 24.d0 - betaval = DT / 2.d0 + DT**2 * thetainv / 3.d0 + DT**3 * thetainv**2 / 8.d0 + DT**4 * thetainv**3 / 24.d0 - gammaval = DT / 2.d0 + DT**2 * thetainv / 6.d0 + DT**3 * thetainv**2 / 24.d0 endif ! allocate memory variables for permittivity attenuation (electromagnetic media) @@ -514,48 +639,6 @@ subroutine prepare_attenuation() alphaval_em(NGLLX,NGLLZ,nspec,2), & betaval_em(NGLLX,NGLLZ,nspec,2), & gammaval_em(NGLLX,NGLLZ,nspec,2),stat=ier) - if (ier /= 0) stop 'Error allocating rx_permattenuation,.. arrays' - - ! initialize memory variables for attenuation - rx_permattenuation(:,:,:) = 0.d0 - rz_permattenuation(:,:,:) = 0.d0 - permx(:,:,:) = 0.d0 - permz(:,:,:) = 0.d0 - - ! if source is not a Dirac or Heavyside then f0_attenuation is f0 of the first - ! source - if (.not. (time_function_type(1) == 4 .or. time_function_type(1) == 5)) then - f0_electromagnetic = f0_source(1) - endif - - ! precompute Runge Kutta coefficients if viscous attenuation - ! viscous attenuation is implemented following the memory variable - ! formulation of - ! J. M. Carcione Wave fields in real media: wave propagation in anisotropic, - ! anelastic and porous media, Elsevier, p. 304-305, 2007 - do ispec = 1,nspec - do j = 1,NGLLZ - do i = 1,NGLLX - - tau_e(i,j,ispec,1) = (sqrt(Qe11_electromagnetic(kmato(ispec))**2+1.d0) +1.d0)/ & - (2.d0*pi*f0_electromagnetic*Qe11_electromagnetic(kmato(ispec))) - tau_e(i,j,ispec,2) = (sqrt(Qe33_electromagnetic(kmato(ispec))**2+1.d0) +1.d0)/ & - (2.d0*pi*f0_electromagnetic*Qe33_electromagnetic(kmato(ispec))) - tau_d(i,j,ispec,1) = (sqrt(Qe11_electromagnetic(kmato(ispec))**2+1.d0) -1.d0)/ & - (2.d0*pi*f0_electromagnetic*Qe11_electromagnetic(kmato(ispec))) - tau_d(i,j,ispec,2) = (sqrt(Qe33_electromagnetic(kmato(ispec))**2+1.d0) -1.d0)/ & - (2.d0*pi*f0_electromagnetic*Qe33_electromagnetic(kmato(ispec))) - - tauinv(:) = - 1.d0 / tau_d(i,j,ispec,:) - alphaval_em(i,j,ispec,:) = 1.d0 + deltat*tauinv(:) + deltat**2*tauinv(:)**2 / 2.d0 & - + deltat**3*tauinv(:)**3 / 6.d0 + deltat**4*tauinv(:)**4 / 24.d0 - betaval_em(i,j,ispec,:) = deltat / 2.d0 + deltat**2*tauinv(:) / 3.d0 + deltat**3*tauinv(:)**2 / 8.d0 & - + deltat**4*tauinv(:)**3 / 24.d0 - gammaval_em(i,j,ispec,:) = deltat / 2.d0 + deltat**2*tauinv(:) / 6.d0 + deltat**3*tauinv(:)**2 / 24.d0 - - enddo - enddo - enddo else ! dummy allocation ! (rx_permattenuation,rz_permattenuation allocation needed for subroutine call) @@ -567,14 +650,16 @@ subroutine prepare_attenuation() tau_d(1,1,1,1), & alphaval_em(1,1,1,1), & betaval_em(1,1,1,1), & - gammaval_em(1,1,1,1)) + gammaval_em(1,1,1,1),stat=ier) endif + if (ier /= 0) stop 'Error allocating rx_permattenuation,.. arrays' + ! initialize memory variables for attenuation + rx_permattenuation(:,:,:) = 0.d0 + rz_permattenuation(:,:,:) = 0.d0 + permx(:,:,:) = 0.d0 + permz(:,:,:) = 0.d0 - ! synchronizes all processes - call synchronize_all() - - end subroutine prepare_attenuation - + end subroutine prepare_attenuation_global_arrays ! !------------------------------------------------------------------------------------- @@ -610,7 +695,6 @@ subroutine prepare_attenuation_with_PML() ! user output if (myrank == 0) then write(IMAIN,*) ' preparing attenuation within PML region' - write(IMAIN,*) call flush_IMAIN() endif