[cig-commits] [commit] devel: updates CPML routines for convolution and pml_compute_memory_variables_acoustic for slightly faster calculations (c9c0282)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Fri Nov 28 14:14:10 PST 2014
Repository : https://github.com/geodynamics/specfem3d
On branch : devel
Link : https://github.com/geodynamics/specfem3d/compare/9a84d06c76e869f5276019e4f84affce23830a4d...8dce71b713c1fd0b510e7538fea0f2307c7b29e8
>---------------------------------------------------------------
commit c9c0282f8874b26dee7c152e5b242b452764e885
Author: daniel peter <peterda at ethz.ch>
Date: Thu Nov 27 16:59:00 2014 +0100
updates CPML routines for convolution and pml_compute_memory_variables_acoustic for slightly faster calculations
>---------------------------------------------------------------
c9c0282f8874b26dee7c152e5b242b452764e885
src/specfem3D/pml_compute_accel_contribution.f90 | 134 ++++++++++++++---------
src/specfem3D/pml_compute_memory_variables.f90 | 96 +++++++---------
2 files changed, 126 insertions(+), 104 deletions(-)
diff --git a/src/specfem3D/pml_compute_accel_contribution.f90 b/src/specfem3D/pml_compute_accel_contribution.f90
index c4911b8..20db342 100644
--- a/src/specfem3D/pml_compute_accel_contribution.f90
+++ b/src/specfem3D/pml_compute_accel_contribution.f90
@@ -1021,81 +1021,117 @@ Subroutine compute_convolution_coef(bb, deltat, coef0, coef1, coef2, singularity
use constants, only: CUSTOM_REAL
- logical,parameter :: FIRST_ORDER_CONVOLUTION = .false.
- real(kind=CUSTOM_REAL) :: bb, deltat, coef0, coef1, coef2, time_nplus1, time_n
- integer :: singularity_type
+ implicit none
+
+ real(kind=CUSTOM_REAL),intent(in) :: bb, deltat
+ real(kind=CUSTOM_REAL),intent(out) :: coef0, coef1, coef2
+ integer,intent(in) :: singularity_type
+ real(kind=CUSTOM_REAL),intent(in) :: time_nplus1, time_n
- coef0 = exp(-bb * deltat)
+ ! local parameters
+ logical,parameter :: FIRST_ORDER_CONVOLUTION = .false.
+ real(kind=CUSTOM_REAL) :: bbpow2,bbpow3
+ real(kind=CUSTOM_REAL) :: deltatpow2,deltatpow3,deltatpow4,deltatpow5,deltatpow6,deltat_half
+ real(kind=CUSTOM_REAL) :: prod1,prod1_half
+
+ ! permanent factors (avoids divisions which are computationally expensive)
+ real(kind=CUSTOM_REAL),parameter :: ONE_OVER_8 = 0.125_CUSTOM_REAL
+ real(kind=CUSTOM_REAL),parameter :: ONE_OVER_12 = 1._CUSTOM_REAL / 12._CUSTOM_REAL
+ real(kind=CUSTOM_REAL),parameter :: ONE_OVER_24 = 1._CUSTOM_REAL / 24._CUSTOM_REAL
+ real(kind=CUSTOM_REAL),parameter :: ONE_OVER_48 = 1._CUSTOM_REAL / 48._CUSTOM_REAL
+ real(kind=CUSTOM_REAL),parameter :: ONE_OVER_128 = 0.0078125_CUSTOM_REAL
+ real(kind=CUSTOM_REAL),parameter :: ONE_OVER_384 = 1._CUSTOM_REAL / 384._CUSTOM_REAL
+ real(kind=CUSTOM_REAL),parameter :: ONE_OVER_960 = 1._CUSTOM_REAL / 960._CUSTOM_REAL
+ real(kind=CUSTOM_REAL),parameter :: ONE_OVER_1920 = 1._CUSTOM_REAL / 1920._CUSTOM_REAL
+ real(kind=CUSTOM_REAL),parameter :: SEVEN_OVER_3840 = 7._CUSTOM_REAL / 3840._CUSTOM_REAL
+ real(kind=CUSTOM_REAL),parameter :: FIVE_OVER_11520 = 5._CUSTOM_REAL/11520._CUSTOM_REAL
+
+ ! helper variables
+ bbpow2 = bb**2
+ bbpow3 = bb**3
+
+ deltatpow2 = deltat**2
+ deltatpow3 = deltat**3
+ deltatpow4 = deltat**4
+ deltatpow5 = deltat**5
+ deltatpow6 = deltat**6
+ deltat_half = deltat * 0.5_CUSTOM_REAL
+
+ prod1 = bb * deltat
+ prod1_half = prod1 * 0.5_CUSTOM_REAL
+
+ ! calculates coefficients
+ coef0 = exp(-prod1)
if (singularity_type == 0)then
if (abs(bb) >= 1.e-5_CUSTOM_REAL) then
if (FIRST_ORDER_CONVOLUTION) then
- coef1 = (1._CUSTOM_REAL - exp(-bb * deltat) ) / bb
+ coef1 = (1._CUSTOM_REAL - exp(-prod1) ) / bb
coef2 = 0._CUSTOM_REAL
else
- coef1 = (1._CUSTOM_REAL - exp(-bb * deltat/2._CUSTOM_REAL) ) / bb
- coef2 = (1._CUSTOM_REAL - exp(-bb * deltat/2._CUSTOM_REAL) ) * exp(-bb * deltat/2._CUSTOM_REAL) / bb
+ coef1 = (1._CUSTOM_REAL - exp(-prod1_half) ) / bb
+ coef2 = (1._CUSTOM_REAL - exp(-prod1_half) ) * exp(-prod1_half) / bb
endif
else
if (FIRST_ORDER_CONVOLUTION) then
coef1 = deltat
coef2 = 0._CUSTOM_REAL
else
- coef1 = deltat/2._CUSTOM_REAL + &
- (- deltat**2*bb/8._CUSTOM_REAL + &
- (deltat**3*bb**2/48._CUSTOM_REAL - &
- deltat**4*bb**3/384._CUSTOM_REAL))
- coef2 = deltat/2._CUSTOM_REAL + &
- (- 3._CUSTOM_REAL*deltat**2*bb/8._CUSTOM_REAL + &
- (7._CUSTOM_REAL*deltat**3*bb**2/48._CUSTOM_REAL - &
- 5._CUSTOM_REAL*deltat**4*bb**3/128._CUSTOM_REAL))
+ coef1 = deltat_half + &
+ (- deltatpow2*bb*ONE_OVER_8 + &
+ (deltatpow3*bbpow2*ONE_OVER_48 - &
+ deltatpow4*bbpow3*ONE_OVER_384))
+ coef2 = deltat_half + &
+ (- 3._CUSTOM_REAL*deltatpow2*bb*ONE_OVER_8 + &
+ (7._CUSTOM_REAL*deltatpow3*bbpow2*ONE_OVER_48 - &
+ 5._CUSTOM_REAL*deltatpow4*bbpow3*ONE_OVER_128))
endif
endif
else if (singularity_type == 1)then
if (abs(bb) >= 1.e-5_CUSTOM_REAL) then
- coef1 = (1._CUSTOM_REAL - exp(-bb * deltat/2._CUSTOM_REAL) ) / bb
- coef1 = time_nplus1 * coef1 + (deltat/2._CUSTOM_REAL*exp(-bb * deltat/2._CUSTOM_REAL) - coef1) / bb
+ coef1 = (1._CUSTOM_REAL - exp(-prod1_half) ) / bb
+ coef1 = time_nplus1 * coef1 + (deltat_half*exp(-prod1_half) - coef1) / bb
- coef2 = (1._CUSTOM_REAL - exp(-bb * deltat/2._CUSTOM_REAL) ) * exp(-bb * deltat/2._CUSTOM_REAL) / bb
- coef2 = time_n * coef2 + ( deltat/2._CUSTOM_REAL*exp(-bb * deltat/2._CUSTOM_REAL) - coef2) / bb
+ coef2 = (1._CUSTOM_REAL - exp(-prod1_half) ) * exp(-prod1_half) / bb
+ coef2 = time_n * coef2 + ( deltat_half*exp(-prod1_half) - coef2) / bb
else
- coef1 = -deltat**2/8._CUSTOM_REAL + (deltat**3*bb/24._CUSTOM_REAL + &
- (- deltat**4*bb**2/128._CUSTOM_REAL + deltat**5*bb**3/960._CUSTOM_REAL)) + &
- time_nplus1* (deltat/2._CUSTOM_REAL + (- deltat**2*bb/8._CUSTOM_REAL + &
- (deltat**3*bb**2/48._CUSTOM_REAL - deltat**4*bb**3/384._CUSTOM_REAL)))
- coef2 = deltat**2/8._CUSTOM_REAL + (-deltat**3*bb/12._CUSTOM_REAL + &
- (11._CUSTOM_REAL*deltat**4*bb**2/384._CUSTOM_REAL - 13._CUSTOM_REAL*deltat**5*bb**3/1920._CUSTOM_REAL)) + &
- time_n * (deltat/2._CUSTOM_REAL + (- 3._CUSTOM_REAL*deltat**2*bb/8._CUSTOM_REAL + &
- (7._CUSTOM_REAL*deltat**3*bb**2/48._CUSTOM_REAL - &
- 5._CUSTOM_REAL*deltat**4*bb**3/128._CUSTOM_REAL)))
+ coef1 = -deltatpow2*ONE_OVER_8 + (deltatpow3*bb*ONE_OVER_24 + &
+ (- deltatpow4*bbpow2*ONE_OVER_128 + deltatpow5*bbpow3*ONE_OVER_960)) + &
+ time_nplus1* (deltat_half + (- deltatpow2*bb*ONE_OVER_8 + &
+ (deltatpow3*bbpow2*ONE_OVER_48 - deltatpow4*bbpow3*ONE_OVER_384)))
+ coef2 = deltatpow2*ONE_OVER_8 + (-deltatpow3*bb*ONE_OVER_12 + &
+ (11._CUSTOM_REAL*deltatpow4*bbpow2*ONE_OVER_384 - 13._CUSTOM_REAL*deltatpow5*bbpow3*ONE_OVER_1920)) + &
+ time_n * (deltat_half + (- 3._CUSTOM_REAL*deltatpow2*bb*ONE_OVER_8 + &
+ (7._CUSTOM_REAL*deltatpow3*bbpow2*ONE_OVER_48 - &
+ 5._CUSTOM_REAL*deltatpow4*bbpow3*ONE_OVER_128)))
endif
else if (singularity_type == 2)then
if (abs(bb) >= 1.e-5_CUSTOM_REAL) then
- coef1 = (1._CUSTOM_REAL - exp(-bb * deltat/2._CUSTOM_REAL) ) / bb
- coef1 = time_nplus1**2 * coef1 + (time_nplus1*(-2._CUSTOM_REAL/bb*coef1 + deltat/bb*exp(-bb * deltat/2._CUSTOM_REAL)) + &
- ((2._CUSTOM_REAL/bb**2* coef1 - exp(-bb * deltat/2._CUSTOM_REAL)*deltat/bb**2) - &
- (deltat/2._CUSTOM_REAL)**2/bb*exp(-bb * deltat/2._CUSTOM_REAL)))
+ coef1 = (1._CUSTOM_REAL - exp(-prod1_half) ) / bb
+ coef1 = time_nplus1**2 * coef1 + (time_nplus1*(-2._CUSTOM_REAL/bb*coef1 + deltat/bb*exp(-prod1_half)) + &
+ ((2._CUSTOM_REAL/bbpow2* coef1 - exp(-prod1_half)*deltat/bbpow2) - &
+ (deltat_half)**2/bb*exp(-prod1_half)))
- coef2 = (1._CUSTOM_REAL - exp(-bb * deltat/2._CUSTOM_REAL) ) * exp(-bb * deltat/2._CUSTOM_REAL) / bb
+ coef2 = (1._CUSTOM_REAL - exp(-prod1_half) ) * exp(-prod1_half) / bb
coef2 = time_n**2 * coef2 + &
- (time_n*(-2._CUSTOM_REAL/bb*coef2 + deltat/bb*exp(-bb * deltat/2._CUSTOM_REAL)) + &
- ((2._CUSTOM_REAL/bb**2* coef2 - exp(-bb * deltat/2._CUSTOM_REAL)*deltat/bb**2) - &
- (deltat/2._CUSTOM_REAL)**2/bb*exp(-bb * deltat/2._CUSTOM_REAL)))
+ (time_n*(-2._CUSTOM_REAL/bb*coef2 + deltat/bb*exp(-prod1_half)) + &
+ ((2._CUSTOM_REAL/bbpow2* coef2 - exp(-prod1_half)*deltat/bbpow2) - &
+ (deltat_half)**2/bb*exp(-prod1_half)))
else
- coef1 = deltat**3/24._CUSTOM_REAL + (- deltat**4*bb*3._CUSTOM_REAL/192._CUSTOM_REAL &
- + (deltat**5*bb**2*3._CUSTOM_REAL/960._CUSTOM_REAL - deltat**6*bb**3*5._CUSTOM_REAL/11520._CUSTOM_REAL)) + &
- time_nplus1*2._CUSTOM_REAL*(-deltat**2/8._CUSTOM_REAL + (deltat**3*bb/24._CUSTOM_REAL + &
- (- deltat**4*bb**2/128._CUSTOM_REAL + deltat**5*bb**3/960._CUSTOM_REAL))) + &
- time_nplus1**2*(deltat/2._CUSTOM_REAL + (- deltat**2*bb/8._CUSTOM_REAL + &
- (deltat**3*bb**2/48._CUSTOM_REAL - deltat**4*bb**3/384._CUSTOM_REAL)))
- coef2 = deltat**3/24._CUSTOM_REAL + (- deltat**4*bb*5._CUSTOM_REAL/192._CUSTOM_REAL &
- + (deltat**5*bb**2*8._CUSTOM_REAL/960._CUSTOM_REAL - deltat**6*bb**3*7._CUSTOM_REAL/3840._CUSTOM_REAL)) + &
- time_n*2._CUSTOM_REAL*(deltat**2/8._CUSTOM_REAL + (-deltat**3*bb/12._CUSTOM_REAL + &
- (11._CUSTOM_REAL*deltat**4*bb**2/384._CUSTOM_REAL - 13._CUSTOM_REAL*deltat**5*bb**3/1920._CUSTOM_REAL))) + &
- time_n**2 * (deltat/2._CUSTOM_REAL + (- 3._CUSTOM_REAL*deltat**2*bb/8._CUSTOM_REAL + &
- (7._CUSTOM_REAL*deltat**3*bb**2/48._CUSTOM_REAL - &
- 5._CUSTOM_REAL*deltat**4*bb**3/128._CUSTOM_REAL)))
+ coef1 = deltatpow3*ONE_OVER_24 + (- deltatpow4*bb*3._CUSTOM_REAL/192._CUSTOM_REAL &
+ + (deltatpow5*bbpow2*3._CUSTOM_REAL*ONE_OVER_960 - deltatpow6*bbpow3*FIVE_OVER_11520)) + &
+ time_nplus1*2._CUSTOM_REAL*(-deltatpow2*ONE_OVER_8 + (deltatpow3*bb*ONE_OVER_24 + &
+ (- deltatpow4*bbpow2*ONE_OVER_128 + deltatpow5*bbpow3*ONE_OVER_960))) + &
+ time_nplus1**2*(deltat_half + (- deltatpow2*bb*ONE_OVER_8 + &
+ (deltatpow3*bbpow2*ONE_OVER_48 - deltatpow4*bbpow3*ONE_OVER_384)))
+ coef2 = deltatpow3*ONE_OVER_24 + (- deltatpow4*bb*5._CUSTOM_REAL/192._CUSTOM_REAL &
+ + (deltatpow5*bbpow2*8._CUSTOM_REAL*ONE_OVER_960 - deltatpow6*bbpow3*SEVEN_OVER_3840)) + &
+ time_n*2._CUSTOM_REAL*(deltatpow2*ONE_OVER_8 + (-deltatpow3*bb*ONE_OVER_12 + &
+ (11._CUSTOM_REAL*deltatpow4*bbpow2*ONE_OVER_384 - 13._CUSTOM_REAL*deltatpow5*bbpow3*ONE_OVER_1920))) + &
+ time_n**2 * (deltat_half + (- 3._CUSTOM_REAL*deltatpow2*bb*ONE_OVER_8 + &
+ (7._CUSTOM_REAL*deltatpow3*bbpow2*ONE_OVER_48 - &
+ 5._CUSTOM_REAL*deltatpow4*bbpow3*ONE_OVER_128)))
endif
else
stop "error in singularity_type in compute_convolution_coefficient"
diff --git a/src/specfem3D/pml_compute_memory_variables.f90 b/src/specfem3D/pml_compute_memory_variables.f90
index 6c117aa..46846d3 100644
--- a/src/specfem3D/pml_compute_memory_variables.f90
+++ b/src/specfem3D/pml_compute_memory_variables.f90
@@ -393,7 +393,9 @@ subroutine pml_compute_memory_variables_acoustic(ispec,ispec_CPML,temp1,temp2,te
real(kind=CUSTOM_REAL) :: coef0_1,coef1_1,coef2_1,coef0_2,coef1_2,coef2_2,coef0_3,coef1_3,coef2_3
integer :: CPML_region_local
integer :: singularity_type_2,singularity_type_3
- real(kind=CUSTOM_REAL) :: kappa_x,kappa_y,kappa_z,d_x,d_y,d_z,alpha_x,alpha_y,alpha_z
+ real(kind=CUSTOM_REAL) :: kappa_x,kappa_y,kappa_z
+ real(kind=CUSTOM_REAL) :: d_x,d_y,d_z,alpha_x,alpha_y,alpha_z
+
CPML_region_local = CPML_regions(ispec_CPML)
time_nplus1 = (it-1) * deltat
@@ -402,7 +404,6 @@ subroutine pml_compute_memory_variables_acoustic(ispec,ispec_CPML,temp1,temp2,te
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
- !---------------------- A6, A7, A8, A9 --------------------------
kappa_x = k_store_x(i,j,k,ispec_CPML)
kappa_y = k_store_y(i,j,k,ispec_CPML)
kappa_z = k_store_z(i,j,k,ispec_CPML)
@@ -413,6 +414,7 @@ subroutine pml_compute_memory_variables_acoustic(ispec,ispec_CPML,temp1,temp2,te
alpha_y = alpha_store_y(i,j,k,ispec_CPML)
alpha_z = alpha_store_z(i,j,k,ispec_CPML)
+ !---------------------- A6, A7, A8, A9 --------------------------
call lijk_parameter_computation(time_nplus1,deltat,&
kappa_z,d_z,alpha_z,kappa_y,d_y,alpha_y,kappa_x,d_x,alpha_x,&
CPML_region_local,231,A6,A7,A8,A9,&
@@ -429,16 +431,6 @@ subroutine pml_compute_memory_variables_acoustic(ispec,ispec_CPML,temp1,temp2,te
coef1_3 * PML_dpotential_dxl_new(i,j,k) + coef2_3 * PML_dpotential_dxl_old(i,j,k)
!---------------------- A10,A11,A12,A13 --------------------------
- kappa_x = k_store_x(i,j,k,ispec_CPML)
- kappa_y = k_store_y(i,j,k,ispec_CPML)
- kappa_z = k_store_z(i,j,k,ispec_CPML)
- d_x = d_store_x(i,j,k,ispec_CPML)
- d_y = d_store_y(i,j,k,ispec_CPML)
- d_z = d_store_z(i,j,k,ispec_CPML)
- alpha_x = alpha_store_x(i,j,k,ispec_CPML)
- alpha_y = alpha_store_y(i,j,k,ispec_CPML)
- alpha_z = alpha_store_z(i,j,k,ispec_CPML)
-
call lijk_parameter_computation(time_nplus1,deltat,&
kappa_x,d_x,alpha_x,kappa_z,d_z,alpha_z,kappa_y,d_y,alpha_y,&
CPML_region_local,132,A10,A11,A12,A13,&
@@ -455,16 +447,6 @@ subroutine pml_compute_memory_variables_acoustic(ispec,ispec_CPML,temp1,temp2,te
coef1_3 * PML_dpotential_dyl_new(i,j,k) + coef2_3 * PML_dpotential_dyl_old(i,j,k)
!---------------------- A14,A15,A16,A17 --------------------------
- kappa_x = k_store_x(i,j,k,ispec_CPML)
- kappa_y = k_store_y(i,j,k,ispec_CPML)
- kappa_z = k_store_z(i,j,k,ispec_CPML)
- d_x = d_store_x(i,j,k,ispec_CPML)
- d_y = d_store_y(i,j,k,ispec_CPML)
- d_z = d_store_z(i,j,k,ispec_CPML)
- alpha_x = alpha_store_x(i,j,k,ispec_CPML)
- alpha_y = alpha_store_y(i,j,k,ispec_CPML)
- alpha_z = alpha_store_z(i,j,k,ispec_CPML)
-
call lijk_parameter_computation(time_nplus1,deltat,&
kappa_x,d_x,alpha_x,kappa_y,d_y,alpha_y,kappa_z,d_z,alpha_z,&
CPML_region_local,123,A14,A15,A16,A17,&
@@ -480,6 +462,7 @@ subroutine pml_compute_memory_variables_acoustic(ispec,ispec_CPML,temp1,temp2,te
rmemory_dpotential_dzl(i,j,k,ispec_CPML,3) = coef0_3 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,3) + &
coef1_3 * PML_dpotential_dzl_new(i,j,k) + coef2_3 * PML_dpotential_dzl_old(i,j,k)
+ ! derivatives
dpotentialdxl = A6 * PML_dpotential_dxl(i,j,k) + &
A7 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,1) + &
A8 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,2) + &
@@ -689,9 +672,10 @@ subroutine lijk_parameter_computation(timeval,deltat,kappa_x,d_x,alpha_x,kappa_y
implicit none
real(kind=CUSTOM_REAL), intent(in) :: timeval,deltat
- real(kind=CUSTOM_REAL) :: kappa_x,d_x,alpha_x,beta_x, &
- kappa_y,d_y,alpha_y,beta_y, &
- kappa_z,d_z,alpha_z,beta_z
+ real(kind=CUSTOM_REAL), intent(in) :: kappa_x,d_x,alpha_x, &
+ kappa_y,d_y,alpha_y, &
+ kappa_z,d_z,alpha_z
+
integer, intent(in) :: CPML_region_local,index_ijk
real(kind=CUSTOM_REAL), intent(out) :: A_0,A_1,A_2,A_3
@@ -700,6 +684,8 @@ subroutine lijk_parameter_computation(timeval,deltat,kappa_x,d_x,alpha_x,kappa_y
integer, intent(out) :: singularity_type_2,singularity_type_3
!local variable
+ real(kind=CUSTOM_REAL) :: a_x,a_y
+ real(kind=CUSTOM_REAL) :: beta_x,beta_y,beta_z
real(kind=CUSTOM_REAL) :: time_nplus1,time_n
real(kind=CUSTOM_REAL) :: bar_A_0,bar_A_1,bar_A_2,bar_A_3,alpha_0
@@ -775,11 +761,11 @@ subroutine lijk_parameter_computation(timeval,deltat,kappa_x,d_x,alpha_x,kappa_y
!----------------A1,2,3-------------------------
alpha_0 = max(alpha_x,alpha_y)
- alpha_x = alpha_0
- alpha_y = alpha_0
+ a_x = alpha_0
+ a_y = alpha_0
- beta_x = alpha_x + d_x / kappa_x
- beta_y = alpha_y + d_y / kappa_y
+ beta_x = a_x + d_x / kappa_x
+ beta_y = a_y + d_y / kappa_y
beta_z = alpha_z + d_z / kappa_z
bar_A_1 = bar_A_0 * (-2._CUSTOM_REAL * alpha_0**3 + beta_x*beta_y*beta_z &
@@ -799,8 +785,8 @@ subroutine lijk_parameter_computation(timeval,deltat,kappa_x,d_x,alpha_x,kappa_y
singularity_type_2 = 1 ! 0 means no singularity, 1 means first order singularity, 2 means second order singularity
singularity_type_3 = 0 ! 0 means no singularity, 1 means first order singularity, 2 means second order singularity
- call compute_convolution_coef(alpha_x, deltat, coef0_1, coef1_1, coef2_1, 0, time_nplus1,time_n)
- call compute_convolution_coef(alpha_y, deltat, coef0_2, coef1_2, coef2_2, singularity_type_2, time_nplus1,time_n)
+ call compute_convolution_coef(a_x, deltat, coef0_1, coef1_1, coef2_1, 0, time_nplus1,time_n)
+ call compute_convolution_coef(a_y, deltat, coef0_2, coef1_2, coef2_2, singularity_type_2, time_nplus1,time_n)
call compute_convolution_coef(beta_z, deltat, coef0_3, coef1_3, coef2_3, singularity_type_3, time_nplus1,time_n)
else if (abs(alpha_x-alpha_y) >= 1.e-5_CUSTOM_REAL .and. abs(alpha_x-beta_z) < 1.e-5_CUSTOM_REAL &
@@ -808,9 +794,9 @@ subroutine lijk_parameter_computation(timeval,deltat,kappa_x,d_x,alpha_x,kappa_y
!----------------A1,2,3-------------------------
alpha_0 = max(alpha_x,beta_z)
- alpha_x = alpha_0
+ a_x = alpha_0
- beta_x = alpha_x + d_x / kappa_x
+ beta_x = a_x + d_x / kappa_x
beta_y = alpha_y + d_y / kappa_y
beta_z = alpha_0
@@ -831,7 +817,7 @@ subroutine lijk_parameter_computation(timeval,deltat,kappa_x,d_x,alpha_x,kappa_y
singularity_type_2 = 0 ! 0 means no singularity, 1 means first order singularity, 2 means second order singularity
singularity_type_3 = 1 ! 0 means no singularity, 1 means first order singularity, 2 means second order singularity
- call compute_convolution_coef(alpha_x, deltat, coef0_1, coef1_1, coef2_1, 0, time_nplus1,time_n)
+ call compute_convolution_coef(a_x, deltat, coef0_1, coef1_1, coef2_1, 0, time_nplus1,time_n)
call compute_convolution_coef(alpha_y, deltat, coef0_2, coef1_2, coef2_2, singularity_type_2, time_nplus1,time_n)
call compute_convolution_coef(beta_z, deltat, coef0_3, coef1_3, coef2_3, singularity_type_3, time_nplus1,time_n)
@@ -840,10 +826,10 @@ subroutine lijk_parameter_computation(timeval,deltat,kappa_x,d_x,alpha_x,kappa_y
!----------------A1,2,3-------------------------
alpha_0 = max(alpha_y,beta_z)
- alpha_y = alpha_0
+ a_y = alpha_0
beta_x = alpha_x + d_x / kappa_x
- beta_y = alpha_y + d_y / kappa_y
+ beta_y = a_y + d_y / kappa_y
beta_z = alpha_0
bar_A_1 = - bar_A_0 * (alpha_x-alpha_z) * (alpha_x - beta_x) * (alpha_x - beta_y) / &
@@ -864,7 +850,7 @@ subroutine lijk_parameter_computation(timeval,deltat,kappa_x,d_x,alpha_x,kappa_y
singularity_type_3 = 1 ! 0 means no singularity, 1 means first order singularity, 2 means second order singularity
call compute_convolution_coef(alpha_x, deltat, coef0_1, coef1_1, coef2_1, 0, time_nplus1,time_n)
- call compute_convolution_coef(alpha_y, deltat, coef0_2, coef1_2, coef2_2, singularity_type_2, time_nplus1,time_n)
+ call compute_convolution_coef(a_y, deltat, coef0_2, coef1_2, coef2_2, singularity_type_2, time_nplus1,time_n)
call compute_convolution_coef(beta_z, deltat, coef0_3, coef1_3, coef2_3, singularity_type_3, time_nplus1,time_n)
else if ((abs(alpha_x-alpha_y) < 1.e-5_CUSTOM_REAL .and. abs(alpha_x-beta_z) < 1.e-5_CUSTOM_REAL &
@@ -881,11 +867,11 @@ subroutine lijk_parameter_computation(timeval,deltat,kappa_x,d_x,alpha_x,kappa_y
!----------------A1,2,3-------------------------
alpha_0 = max(alpha_x,alpha_y,beta_z)
- alpha_x = alpha_0
- alpha_y = alpha_0
+ a_x = alpha_0
+ a_y = alpha_0
- beta_x = alpha_x + d_x / kappa_x
- beta_y = alpha_y + d_y / kappa_y
+ beta_x = a_x + d_x / kappa_x
+ beta_y = a_y + d_y / kappa_y
beta_z = alpha_0
bar_A_1 = bar_A_0 * (-3._CUSTOM_REAL * alpha_0 + alpha_z + beta_x + beta_y)
@@ -900,8 +886,8 @@ subroutine lijk_parameter_computation(timeval,deltat,kappa_x,d_x,alpha_x,kappa_y
singularity_type_2 = 1 ! 0 means no singularity, 1 means first order singularity, 2 means second order singularity
singularity_type_3 = 2 ! 0 means no singularity, 1 means first order singularity, 2 means second order singularity
- call compute_convolution_coef(alpha_x, deltat, coef0_1, coef1_1, coef2_1, 0, time_nplus1,time_n)
- call compute_convolution_coef(alpha_y, deltat, coef0_2, coef1_2, coef2_2, singularity_type_2, time_nplus1,time_n)
+ call compute_convolution_coef(a_x, deltat, coef0_1, coef1_1, coef2_1, 0, time_nplus1,time_n)
+ call compute_convolution_coef(a_y, deltat, coef0_2, coef1_2, coef2_2, singularity_type_2, time_nplus1,time_n)
call compute_convolution_coef(beta_z, deltat, coef0_3, coef1_3, coef2_3, singularity_type_3, time_nplus1,time_n)
else
@@ -940,10 +926,10 @@ subroutine lijk_parameter_computation(timeval,deltat,kappa_x,d_x,alpha_x,kappa_y
!----------------A1,2,3-------------------------
alpha_0 = max(alpha_y,beta_z)
- alpha_y = alpha_0
+ a_y = alpha_0
beta_x = alpha_x + d_x / kappa_x
- beta_y = alpha_y + d_y / kappa_y
+ beta_y = a_y + d_y / kappa_y
beta_z = alpha_0
bar_A_1 = 0._CUSTOM_REAL
@@ -958,7 +944,7 @@ subroutine lijk_parameter_computation(timeval,deltat,kappa_x,d_x,alpha_x,kappa_y
singularity_type_3 = 1 ! 0 means no singularity, 1 means first order singularity, 2 means second order singularity
call compute_convolution_coef(alpha_x, deltat, coef0_1, coef1_1, coef2_1, 0, time_nplus1,time_n)
- call compute_convolution_coef(alpha_y, deltat, coef0_2, coef1_2, coef2_2, singularity_type_2, time_nplus1,time_n)
+ call compute_convolution_coef(a_y, deltat, coef0_2, coef1_2, coef2_2, singularity_type_2, time_nplus1,time_n)
call compute_convolution_coef(beta_z, deltat, coef0_3, coef1_3, coef2_3, singularity_type_3, time_nplus1,time_n)
else
@@ -997,9 +983,9 @@ subroutine lijk_parameter_computation(timeval,deltat,kappa_x,d_x,alpha_x,kappa_y
!----------------A1,2,3-------------------------
alpha_0 = max(alpha_x,beta_z)
- alpha_x = alpha_0
+ a_x = alpha_0
- beta_x = alpha_x + d_x / kappa_x
+ beta_x = a_x + d_x / kappa_x
beta_y = alpha_y + d_y / kappa_y
beta_z = alpha_0
@@ -1014,7 +1000,7 @@ subroutine lijk_parameter_computation(timeval,deltat,kappa_x,d_x,alpha_x,kappa_y
singularity_type_2 = 0 ! 0 means no singularity, 1 means first order singularity, 2 means second order singularity
singularity_type_3 = 1 ! 0 means no singularity, 1 means first order singularity, 2 means second order singularity
- call compute_convolution_coef(alpha_x, deltat, coef0_1, coef1_1, coef2_1, 0, time_nplus1,time_n)
+ call compute_convolution_coef(a_x, deltat, coef0_1, coef1_1, coef2_1, 0, time_nplus1,time_n)
call compute_convolution_coef(alpha_y, deltat, coef0_2, coef1_2, coef2_2, singularity_type_2, time_nplus1,time_n)
call compute_convolution_coef(beta_z, deltat, coef0_3, coef1_3, coef2_3, singularity_type_3, time_nplus1,time_n)
@@ -1055,11 +1041,11 @@ subroutine lijk_parameter_computation(timeval,deltat,kappa_x,d_x,alpha_x,kappa_y
!----------------A1,2,3-------------------------
alpha_0 = max(alpha_x,alpha_y)
- alpha_x = alpha_0
- alpha_y = alpha_0
+ a_x = alpha_0
+ a_y = alpha_0
- beta_x = alpha_x + d_x / kappa_x
- beta_y = alpha_y + d_y / kappa_y
+ beta_x = a_x + d_x / kappa_x
+ beta_y = a_y + d_y / kappa_y
beta_z = alpha_z + d_z / kappa_z
bar_A_1 = bar_A_0 * (-2._CUSTOM_REAL * alpha_0 + (beta_x + beta_y))
@@ -1073,8 +1059,8 @@ subroutine lijk_parameter_computation(timeval,deltat,kappa_x,d_x,alpha_x,kappa_y
singularity_type_2 = 1 ! 0 means no singularity, 1 means first order singularity, 2 means second order singularity
singularity_type_3 = 0 ! 0 means no singularity, 1 means first order singularity, 2 means second order singularity
- call compute_convolution_coef(alpha_x, deltat, coef0_1, coef1_1, coef2_1, 0, time_nplus1,time_n)
- call compute_convolution_coef(alpha_y, deltat, coef0_2, coef1_2, coef2_2, singularity_type_2, time_nplus1,time_n)
+ call compute_convolution_coef(a_x, deltat, coef0_1, coef1_1, coef2_1, 0, time_nplus1,time_n)
+ call compute_convolution_coef(a_y, deltat, coef0_2, coef1_2, coef2_2, singularity_type_2, time_nplus1,time_n)
call compute_convolution_coef(beta_z, deltat, coef0_3, coef1_3, coef2_3, singularity_type_3, time_nplus1,time_n)
else
More information about the CIG-COMMITS
mailing list