Commit b761309a authored by Christophe Geuzaine's avatar Christophe Geuzaine

Merge branch 'PeweClean' into 'master'

Pewe clean

See merge request !39
parents e430729c 30fda5bc
Pipeline #1466 passed with stage
in 12 minutes and 12 seconds
......@@ -1170,7 +1170,11 @@ struct StringXFunction2Nbr F_Function[] = { /* #Par #Arg */
{"ElastodynamicsCylinderWallOut", (CAST)F_ElastodynamicsCylinderWallOut, 5, 1},
{"ElastodynamicsCylinderWallsOut", (CAST)F_ElastodynamicsCylinderWallsOut, 5, 1},
{"ElastoCylinderWallOutAbc",(CAST)F_ElastoCylinderWallOutAbc, 6, 1},
{"ElastoCylinderWallsOutAbc",(CAST)F_ElastoCylinderWallsOutAbc, 6, 1},
{"ElastoCylinderWallOutAbc2",(CAST)F_ElastoCylinderWallOutAbc2, 6, 1},
{"ElastoCylinderWallOutAbc2Pade",(CAST)F_ElastoCylinderWallOutAbc2Pade, 10, 1},
{"ElastoCylinderWallsOutAbc2Pade",(CAST)F_ElastoCylinderWallsOutAbc2Pade, 10, 1},
// F_Raytracing: ray tracing functions
{"CylinderPhase", (CAST)F_CylinderPhase, 0, 1 },
......
......@@ -78,7 +78,10 @@ void F_ElastodynamicsCylinderWallS(F_ARG);
void F_ElastodynamicsCylinderWallOut(F_ARG);
void F_ElastodynamicsCylinderWallsOut(F_ARG);
void F_ElastoCylinderWallOutAbc(F_ARG);
void F_ElastoCylinderWallsOutAbc(F_ARG);
void F_ElastoCylinderWallOutAbc2(F_ARG);
void F_ElastoCylinderWallOutAbc2Pade(F_ARG);
void F_ElastoCylinderWallsOutAbc2Pade(F_ARG);
/* F_Geometry */
......
......@@ -1686,7 +1686,7 @@ void F_AcousticFieldSoftCylinder(F_ARG)
std::complex<double> I(0,1);
double vr=0, vi=0;
#if defined(_OPENMP)
#pragma omp parallel for reduction(+: vr,vi)
//#pragma omp parallel for reduction(+: vr,vi)
#endif
for(int n = nStart ; n < nEnd ; n++){
std::complex<double> HnkR( jn(n,kR), yn(n,kR) );
......
......@@ -33,10 +33,26 @@ extern "C" {
double *lambda, double *mu, double *rho, double *a,
double *b);
extern void cylindrical_wallsoutabc_(double *du, double *dv, double *dut, double *dvt,
double *X, double *Y, double *t, double *omega,
double *lambda, double *mu, double *rho, double *a,
double *b);
extern void cylindrical_walloutabc2_(double *du, double *dv, double *dut, double *dvt,
double *X, double *Y, double *t, double *omega,
double *lambda, double *mu, double *rho, double *a,
double *b);
extern void cylindrical_walloutabc2pade_(double *du, double *dv, double *dut, double *dvt,
double *X, double *Y, double *t, double *omega,
double *lambda, double *mu, double *rho, double *a,
double *b, double *L, double *alpha, double *eps_p, double *eps_s);
extern void cylindrical_wallsoutabc2pade_(double *du, double *dv, double *dut, double *dvt,
double *X, double *Y, double *t, double *omega,
double *lambda, double *mu, double *rho, double *a,
double *b, double *L, double *alpha, double *eps_p, double *eps_s);
}
void F_ElastodynamicsCylinderCavity(F_ARG)
......@@ -120,7 +136,8 @@ void F_ElastodynamicsCylinderWallOut(F_ARG)
double mu = Fct->Para[2];
double rho = Fct->Para[3];
double a = Fct->Para[4];
cylindrical_wallout_(&du_re,&dv_re,&du_im,&dv_im,&X,&Y,&t,&omega,&lambda,&mu,&rho,&a);
cylindrical_wallout_(&du_re,&dv_re,&du_im,&dv_im,&X,&Y,&t,
&omega,&lambda,&mu,&rho,&a);
#else
Message::Error("ElastodynamicsCylinderWallOut requires PeWe");
#endif
......@@ -179,6 +196,31 @@ void F_ElastoCylinderWallOutAbc(F_ARG)
V->Type = VECTOR ;
}
void F_ElastoCylinderWallsOutAbc(F_ARG)
{
double du_re = 0., dv_re = 0., du_im = 0., dv_im = 0.;
#if defined(HAVE_PEWE)
double X = A->Val[0];
double Y = A->Val[1];
double t = 0.;
double omega = Fct->Para[0];
double lambda = Fct->Para[1];
double mu = Fct->Para[2];
double rho = Fct->Para[3];
double a = Fct->Para[4];
double b = Fct->Para[5];
cylindrical_wallsoutabc_(&du_re,&dv_re,&du_im,&dv_im,&X,&Y,&t,&omega,
&lambda,&mu,&rho,&a,&b);
#else
Message::Error("ElastodynamicsCylinderWallsOutABC requires PeWe");
#endif
V->Val[0] = du_re;
V->Val[1] = dv_re;
V->Val[MAX_DIM] = du_im;
V->Val[MAX_DIM+1] = dv_im;
V->Type = VECTOR ;
}
void F_ElastoCylinderWallOutAbc2(F_ARG)
{
double du_re = 0., dv_re = 0., du_im = 0., dv_im = 0.;
......@@ -204,5 +246,62 @@ void F_ElastoCylinderWallOutAbc2(F_ARG)
V->Type = VECTOR ;
}
void F_ElastoCylinderWallOutAbc2Pade(F_ARG)
{
double du_re = 0., dv_re = 0., du_im = 0., dv_im = 0.;
#if defined(HAVE_PEWE)
double X = A->Val[0];
double Y = A->Val[1];
double t = 0.;
double omega = Fct->Para[0];
double lambda = Fct->Para[1];
double mu = Fct->Para[2];
double rho = Fct->Para[3];
double a = Fct->Para[4];
double b = Fct->Para[5];
double L = Fct->Para[6];
double alpha = Fct->Para[7];
double eps_p = Fct->Para[8];
double eps_s = Fct->Para[9];
cylindrical_walloutabc2pade_(&du_re,&dv_re,&du_im,&dv_im,&X,&Y,&t,&omega,
&lambda,&mu,&rho,&a,&b,&L,&alpha,&eps_p,&eps_s);
#else
Message::Error("ElastodynamicsCylinderWallOutABC2_Pade requires PeWe");
#endif
V->Val[0] = du_re;
V->Val[1] = dv_re;
V->Val[MAX_DIM] = du_im;
V->Val[MAX_DIM+1] = dv_im;
V->Type = VECTOR ;
}
void F_ElastoCylinderWallsOutAbc2Pade(F_ARG)
{
double du_re = 0., dv_re = 0., du_im = 0., dv_im = 0.;
#if defined(HAVE_PEWE)
double X = A->Val[0];
double Y = A->Val[1];
double t = 0.;
double omega = Fct->Para[0];
double lambda = Fct->Para[1];
double mu = Fct->Para[2];
double rho = Fct->Para[3];
double a = Fct->Para[4];
double b = Fct->Para[5];
double L = Fct->Para[6];
double alpha = Fct->Para[7];
double eps_p = Fct->Para[8];
double eps_s = Fct->Para[9];
cylindrical_wallsoutabc2pade_(&du_re,&dv_re,&du_im,&dv_im,&X,&Y,&t,&omega,
&lambda,&mu,&rho,&a,&b,&L,&alpha,&eps_p,&eps_s);
#else
Message::Error("ElastodynamicsCylinderWallsOutABC2_Pade requires PeWe");
#endif
V->Val[0] = du_re;
V->Val[1] = dv_re;
V->Val[MAX_DIM] = du_im;
V->Val[MAX_DIM+1] = dv_im;
V->Type = VECTOR ;
}
#undef F_ARG
......@@ -11,8 +11,10 @@ set(SRC
fortran/cylindrical_wallout.f90
fortran/cylindrical_wallsout.f90
fortran/cylindrical_walloutabc.f90
fortran/cylindrical_wallsoutabc.f90
fortran/cylindrical_walloutabc2.f90
fortran/cylindrical_walloutabc2pade.f90
fortran/cylindrical_wallsoutabc2pade.f90
)
file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
......
......@@ -4,7 +4,6 @@
!
! Copyright (C) 2015 Kristoffer Virta & Daniel Appelo
!
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
......@@ -13,17 +12,15 @@
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
! GNU General Public License for more dksils.
!
! You should have received a copy of the GNU General Public License
! along with this program. If not, see <http://www.gnu.org/licenses/>.
!
! Modified by Steven Roman for inclusion in GetDP:
! Modified by Vanessa Mattesi for inclusion in GetDP:
!
! Exact solution for a cylindrical wall (zero displacement on the
! boundary).
! boundary) for a pressure incident wave.
subroutine cylindrical_wallout(du,dv,dut,dvt,X,Y,t,omega,lambda,mu,rho,a)
......@@ -31,149 +28,65 @@ subroutine cylindrical_wallout(du,dv,dut,dvt,X,Y,t,omega,lambda,mu,rho,a)
double precision :: x,y,du,dv,dvt,dut,t
double precision :: r,theta
double precision :: lambda,mu,rho
double precision :: cp,cs,omega,gamma,eta,a
double precision :: phi_0,f_00,f_02,epsilon_0,epsilon_1,c01,c02
double precision :: epsilon_n
double precision :: f_10,f_11,f_12,f_13,f_14,f_15
double complex :: a_00,a_02,b_01,b_04,m11,m12,m21,m22,ab0(2),q
double complex :: a_10,a_11,a_12,a_13,a_14,a_15
double complex :: b_10,b_11,b_12,b_13,b_14,b_15
double complex :: c11,c12,ab1(2),v,cn1,cn2,abn(2)
double complex :: f_n0,f_n1,f_n2,f_n3,f_n4,f_n5
double complex :: a_n0,a_n1,a_n2,a_n3,a_n4,a_n5
double complex :: b_n0,b_n1,b_n2,b_n3,b_n4,b_n5
integer :: n
double precision , parameter :: pi = acos(-1.d0)
! for GetDP
integer :: ns
double precision :: cp,cs,omega,kp,ks,a
double precision :: phi_0,epsilon_n
double complex :: m11,m12,m21,m22,An,Bn,f_1,f_2
double complex :: p,q
integer :: n,ns
double precision , parameter :: pi = 4d0*datan(1d0)
double complex, external :: besselh
! Computes displacement field (du(x,y,t),dy(x,y,t))
! of incoming P waves, scattered P waves and scattered S waves at
! time t = 0. Solution at other times are given by
! du(x,y,t) = du(x,y,0)*exp(1i*omega*t),
! dv(x,y,t) = dv(x,y,0)*exp(1i*omega*t).
double complex, external :: besselh !besselh(order n,kind 1 or 2,k*r) Hankel function of order n of the first or second kind of argument (kr)
double complex, external :: dr_h !dr_h(kind 1 or 2,k,r,n) derivative of the Hankel function of order n of the first or second kind of argument (kr)
double complex, external :: dr_j !dr_j(k,r,n) derivative of the bessel function of first kind of order n of argument (kr)
! Compute radius r and angle theta
r = sqrt(X**2+Y**2)
theta = atan2(Y,X)
r = sqrt(X**2+Y**2);
theta = atan2(Y,X);
! P and S wave speeds
cp = sqrt((lambda+2.d0*mu)/rho)
cs = sqrt(mu/rho)
cp = sqrt((lambda+2.d0*mu)/rho);
cs = sqrt(mu/rho);
! To satisfy elastic wave equation
gamma = omega/cp
eta = omega/cs
! Radius of cylinder
!a = 1.d0
kp = omega/cp;
ks = omega/cs;
! Amplitude of incomming wave
phi_0 = 1.d0
! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Series expansion of solution%%%%%%%%%%%
! disp('COMPUTING COEFICIENTS')
n = 0
f_00 = -gamma*besjn(1,gamma*a)
a_00 = -gamma*besselh(1,1,gamma*a)
b_01 = eta*besselh(1,1,eta*a)
epsilon_0 = 1.d0
m11 = a_00
m12 = 0.d0
m21 = 0.d0
m22 = b_01
c01 = -phi_0*epsilon_0*f_00
c02 = 0.d0
AB0(1) = 1.d0/(m11*m22 - m12*m21)*(m22*c01-m12*c02)
AB0(2) = 1.d0/(m11*m22 - m12*m21)*(-m21*c01+m11*c02)
q = phi_0*epsilon_0*(-gamma)*besjn(1,gamma*r) &
+ AB0(1)*(-gamma)*(besselh(1,1,gamma*r))
n = 1
f_10 = gamma/2.d0*(besjn(0,gamma*a)-besjn(2,gamma*a))
f_11 = -1.d0/a*besjn(1,gamma*a)
a_10 = gamma/2.d0*(besselh(0,1,gamma*a)-besselh(2,1,gamma*a))
a_11 = -1.d0/a*besselh(1,1,gamma*a)
b_10 = 1.d0/a*besselh(1,1,eta*a)
b_11 = -eta/2.d0*(besselh(0,1,eta*a)-besselh(2,1,eta*a))
epsilon_1 = 2.d0
m11 = a_10
m12 = b_10
m21 = a_11
m22 = b_11
c11 = -phi_0*epsilon_1*(0.d0,-1.d0)*f_10
c12 = -phi_0*epsilon_1*(0.d0,-1.d0)*f_11
AB1(1) = 1.d0/(m11*m22 - m12*m21)*(m22*c11-m12*c12)
AB1(2) = 1.d0/(m11*m22 - m12*m21)*(-m21*c11+m11*c12)
q = q + (phi_0*epsilon_1*(0.d0,-1.d0)*gamma/2.d0*(besjn(0,gamma*r)&
-besjn(2,gamma*r)) + AB1(1)*gamma/2.d0*(besselh(0,1,gamma*r)&
-besselh(2,1,gamma*r)) + AB1(2)*besselh(1,1,eta*r)/r)*cos(theta)
v = (-phi_0*epsilon_1*(0.d0,-1.d0)*besjn(1,gamma*r)/r&
- AB1(1)*besselh(1,1,gamma*r)/r &
- AB1(2)*eta/2.d0*(besselh(0,1,eta*r)&
-besselh(2,1,eta*r)))*sin(theta)
! for GetDP
ns = max(eta*a+30, 2*eta);
! for GetDP: instead of 2:24
!$OMP PARALLEL DO PRIVATE(f_n0,f_n1,a_n0,a_n1,b_n0,b_n1,epsilon_n,m11,m12,m21,m22,cn1,cn2,ABn) REDUCTION(+:q,v)
do n = 2,ns
f_n0 = gamma/2.d0*(besjn(n-1,gamma*a)-besjn(n+1,gamma*a))
f_n1 = -dble(n)/a*besjn(n,gamma*a)
a_n0 = gamma/2.d0*(besselh(n-1,1,gamma*a)-besselh(n+1,1,gamma*a))
a_n1 = -dble(n)/a*besselh(n,1,gamma*a)
b_n0 = dble(n)/a*besselh(n,1,eta*a)
b_n1 = -eta/2.d0*(besselh(n-1,1,eta*a)-besselh(n+1,1,eta*a))
epsilon_n = 2.d0
m11 = a_n0
m12 = b_n0
m21 = a_n1
m22 = b_n1
cn1 = -phi_0*epsilon_n*(0.d0,-1.d0)**n * f_n0
cn2 = -phi_0*epsilon_n*(0.d0,-1.d0)**n * f_n1
ABn(1) = 1.d0/(m11*m22 - m12*m21)*(m22*cn1-m12*cn2)
ABn(2) = 1.d0/(m11*m22 - m12*m21)*(-m21*cn1+m11*cn2)
! Make sure magnitude of coefficients are sufficiently small
! max(abs(ABn))
q = q+(phi_0*epsilon_n*(0.d0,-1.d0)**n*gamma/2.d0*(besjn(n-1,gamma*r)&
-besjn(n+1,gamma*r)) + ABn(1)*gamma/2.d0*(besselh(n-1,1,gamma*r)&
-besselh(n+1,1,gamma*r)) &
+ ABn(2)*dble(n)*besselh(n,1,eta*r)/r)*cos(dble(n)*theta);
v = v + (-dble(n)*phi_0*epsilon_n*(0.d0,-1.d0)**n*besjn(n,gamma*r)/r &
- ABn(1)*dble(n)*besselh(n,1,gamma*r)/r &
- ABn(2)*eta/2.d0*(besselh(n-1,1,eta*r)&
-besselh(n+1,1,eta*r)))*sin(dble(n)*theta);
phi_0 = 1.d0;
! size of the series
ns = 2*FLOOR(omega);
!-------------------------------
! Series expansion of solution -
!-------------------------------
! initialisation of the radial and azimutal part of the solution
p=0.d0;
q=0.d0;
!$OMP PARALLEL DO PRIVATE(f_1,f_2,epsilon_n,n,m11,m12,m21,m22,An,Bn) REDUCTION(+:p,q)
do n = 0,ns
If(n.eq.0) Then
epsilon_n = 1.d0;
elseIf (n>0) Then
epsilon_n = 2.d0;
Endif
m11 = dr_h(1,kp,a,n);
m12 = (dble(n)/a)*besselh(int(n),1,ks*a);
m21 =-(dble(n)/a)*besselh(int(n),1,kp*a);
m22 =-dr_h(1,ks,a,n);
f_1 =-phi_0*epsilon_n*((0.d0,-1.d0)**n) *dr_j(kp,a,n);
f_2 = phi_0*epsilon_n*((0.d0,-1.d0)**n) *(dble(n)/a)*besjn(n,kp*a);
An = 1.d0/(m11*m22 - m12*m21)*( m22*f_1-m12*f_2);
Bn = 1.d0/(m11*m22 - m12*m21)*(-m21*f_1+m11*f_2);
p = p+ ( An*dr_h(1,kp,r,n) + Bn*(dble(n)/r)*besselh(int(n),1,ks*r)) *cos(n*theta);
q = q+ (-An*(dble(n)/r)*besselh(int(n),1,kp*r) - Bn*dr_h(1,ks,r,n)) *sin(n*theta);
! if(maxval(abs(abn)) .lt. 1.0d-16) exit
!if(maxval(abs(abn)) .lt. 1.0d-16) exit
end do
!$OMP END PARALLEL DO
! disp('DONE COMPUTING COEFICIENTS')
du = dreal(exp(omega*(0.d0,1.d0)*t)*(cos(theta)*q-sin(theta)*v))
dv = dreal(exp(omega*(0.d0,1.d0)*t)*(sin(theta)*q+cos(theta)*v))
!dut = dreal(omega*(0.d0,1.d0)*exp(omega*(0.d0,1.d0)*t)*(cos(theta)*q-sin(theta)*v))
!dvt = dreal(omega*(0.d0,1.d0)*exp(omega*(0.d0,1.d0)*t)*(sin(theta)*q+cos(theta)*v))
! for GetDP: return imaginary part in dut, dvt
dut = dimag(exp(omega*(0.d0,1.d0)*t)*(cos(theta)*q-sin(theta)*v))
dvt = dimag(exp(omega*(0.d0,1.d0)*t)*(sin(theta)*q+cos(theta)*v))
! return real part in du, dv
du = dreal(exp(omega*(0.d0,1.d0)*t)*(cos(theta)*p-sin(theta)*q))
dv = dreal(exp(omega*(0.d0,1.d0)*t)*(sin(theta)*p+cos(theta)*q))
! return imaginary part in dut, dvt
dut = dimag(exp(omega*(0.d0,1.d0)*t)*(cos(theta)*p-sin(theta)*q))
dvt = dimag(exp(omega*(0.d0,1.d0)*t)*(sin(theta)*p+cos(theta)*q))
end subroutine cylindrical_wallout
This diff is collapsed.
......@@ -88,7 +88,7 @@ subroutine cylindrical_walloutabc2(du,dv,dut,dvt,X,Y,t,omega,lambda,mu,rho,a,b)
q=0
epsilon_1 = 2.0d0
! for GetDP
ns = max(ks*a+30, 3*ks);
ns = 55;
do n = 1,1
......
This diff is collapsed.
......@@ -4,7 +4,6 @@
!
! Copyright (C) 2015 Kristoffer Virta & Daniel Appelo
!
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
......@@ -18,8 +17,6 @@
! You should have received a copy of the GNU General Public License
! along with this program. If not, see <http://www.gnu.org/licenses/>.
!
! Modified by Vanessa Mattesi for inclusion in GetDP:
!
! Exact solution for a cylindrical wall (zero displacement on the
......@@ -31,130 +28,66 @@ subroutine cylindrical_wallsout(du,dv,dut,dvt,X,Y,t,omega,lambda,mu,rho,a)
double precision :: x,y,du,dv,dvt,dut,t
double precision :: r,theta
double precision :: lambda,mu,rho
double precision :: cp,cs,omega,gamma,eta,a
double precision :: psi_0,f_01,f_02,epsilon_0,epsilon_1,c01,c02
double precision :: epsilon_n
double precision :: f_10,f_11,f_12,f_13,f_14,f_15
double complex :: a_00,a_02,b_01,b_04,m11,m12,m21,m22,ab0(2),q
double complex :: a_10,a_11,a_12,a_13,a_14,a_15
double complex :: b_10,b_11,b_12,b_13,b_14,b_15
double complex :: c11,c12,ab1(2),v,cn1,cn2,abn(2)
double complex :: f_n0,f_n1,f_n2,f_n3,f_n4,f_n5
double complex :: a_n0,a_n1,a_n2,a_n3,a_n4,a_n5
double complex :: b_n0,b_n1,b_n2,b_n3,b_n4,b_n5
integer :: n
double precision , parameter :: pi = acos(-1.d0)
! for GetDP
integer :: ns
double precision :: cp,cs,omega,kp,ks,a
double precision :: psi_0,epsilon_n
double complex :: dr_j_p,dr_h_p,dr_h_s,dr_hr_p,dr_hr_s
double complex :: m11,m12,m21,m22,An,Bn,f_1,f_2
double complex :: p,q
integer :: n,ns
double precision , parameter :: pi = 4d0*datan(1d0)
double complex, external :: besselh
! Computes displacement field (du(x,y,t),dy(x,y,t))
! of incoming S waves, scattered P waves and scattered S waves at
! time t = 0. Solution at other times are given by
! du(x,y,t) = du(x,y,0)*exp(1i*omega*t),
! dv(x,y,t) = dv(x,y,0)*exp(1i*omega*t).
double complex, external :: besselh !besselh(order n,kind 1 or 2,k*r) Hankel function of order n of the first or second kind of argument (kr)
double complex, external :: dr_h !dr_h(kind 1 or 2,k,r,n) derivative of the Hankel function of order n of the first or second kind of argument (kr)
double complex, external :: dr_j !dr_j(k,r,n) derivative of the bessel function of first kind of order n of argument (kr)
! Compute radius r and angle theta
r = sqrt(X**2+Y**2)
theta = atan2(Y,X)
r = sqrt(X**2+Y**2);
theta = atan2(Y,X);
! P and S wave speeds
cp = sqrt((lambda+2.d0*mu)/rho)
cs = sqrt(mu/rho)
cp = sqrt((lambda+2.d0*mu)/rho);
cs = sqrt(mu/rho);
! To satisfy elastic wave equation
gamma = omega/cp
eta = omega/cs
! Radius of cylinder
!a = 1.d0
kp = omega/cp;
ks = omega/cs;
! Amplitude of incomming wave
psi_0 = 1.d0
! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Series expansion of solution%%%%%%%%%%%
! disp('COMPUTING COEFICIENTS')
n = 0
f_01= eta*besjn(1,eta*a)
a_00 = -gamma*besselh(1,1,gamma*a)
b_01 = eta*besselh(1,1,eta*a)
epsilon_0 = 1.d0
m11 = a_00
m12 = 0.d0
m21 = 0.d0
m22 = b_01
c01 = 0.d0
c02 = -psi_0*epsilon_0*f_01
psi_0 = 1.d0;
! size of the series
ns = 2*FLOOR(omega);
!-------------------------------
! Series expansion of solution -
!-------------------------------
! initialisation of the radial and azimutal part of the solution
p=0.d0;
q=0.d0;
!$OMP PARALLEL DO PRIVATE(n,f_1,f_2,epsilon_n,m11,m12,m21,m22,An,Bn) REDUCTION(+:p,q)
do n = 0,ns
If (n==0) then
epsilon_n=1.d0;
else if (n>0) then
epsilon_n=2.d0;
endIf
m11 = dr_h(1,kp,a,n);
m12 =-(dble(n)/a)*besselh(n,1,ks*a);
m21 = (dble(n)/a)*besselh(n,1,kp*a);
m22 =-dr_h(1,ks,a,n);
f_1 = psi_0*epsilon_n*((0.d0,-1.d0)**n) *(dble(n)/a)*besjn(n,ks*a);
f_2 = psi_0*epsilon_n*((0.d0,-1.d0)**n) *dr_j(ks,a,n);
An = 1.d0/(m11*m22 - m12*m21)*( m22*f_1-m12*f_2);
Bn = 1.d0/(m11*m22 - m12*m21)*(-m21*f_1+m11*f_2);
p = p+ ( An*dr_h(1,kp,r,n) - Bn*(dble(n)/r)*besselh(n,1,ks*r) ) *sin(n*theta);
q = q+ ( An*(dble(n)/r)*besselh(n,1,kp*r) - Bn*dr_h(1,ks,r,n) ) *cos(n*theta);
AB0(1) = 1.d0/(m11*m22 - m12*m21)*(m22*c01-m12*c02)
AB0(2) = 1.d0/(m11*m22 - m12*m21)*(-m21*c01+m11*c02)
v = AB0(2) * eta*(besselh(1,1,eta*r)) + psi_0*epsilon_0* eta *besjn(1,eta*r)
! for GetDP
ns = max(eta*a+30, 2*eta);
! for GetDP: instead of 2:24
epsilon_1 = 2.d0
do n = 1,1
f_n0 = psi_0*epsilon_1*((0.d0,-1.d0)**n)*(dble(n)/a)*besjn(n,eta*a)
f_n1 = psi_0*epsilon_1*((0.d0,-1.d0)**n)*(eta/2.d0)*(besjn(n-1,eta*a)-besjn(n+1,eta*a))
m11 = (gamma/2.d0)*(besselh(n-1,1,gamma*a)-besselh(n+1,1,gamma*a))
m12 = -(dble(n)/a)*(besselh(int(n),1,eta*a))
m21 = (dble(n)/a)*besselh(int(n),1,gamma*a)
m22 = -(eta/2.d0)*(besselh(n-1,1,eta*a)-besselh(n+1,1,eta*a))
ABn(1) = (1.d0/(m11*m22 - m12*m21))*(m22*f_n0-m12*f_n1)
ABn(2) = (1.d0/(m11*m22 - m12*m21))*(-m21*f_n0+m11*f_n1)
q = q+(-psi_0*epsilon_1*((0.d0,-1.d0)**n)*(dble(n)/r)*besjn(n,eta*r)&
+ ABn(1)*(gamma/2.d0)*(besselh(n-1,1,gamma*r)&
-besselh(n+1,1,gamma*r)) &
- ABn(2)*(dble(n)/r)*besselh(int(n),1,eta*r))*sin(dble(n)*theta)
v = v + (-psi_0*epsilon_1*((0.d0,-1.d0)**n)*(eta/2.d0)*( besjn(n-1,eta*r) &
-besjn(n+1,eta*r)) &
+ ABn(1)*(dble(n)/r)*besselh(int(n),1,gamma*r) &
- ABn(2)*(eta/2.d0)*(besselh(n-1,1,eta*r)&
-besselh(n+1,1,eta*r)))*cos(dble(n)*theta)
end do
!$OMP PARALLEL DO PRIVATE(f_n0,f_n1,a_n0,a_n1,b_n0,b_n1,epsilon_n,m11,m12,m21,m22,cn1,cn2,ABn) REDUCTION(+:q,v)
do n = 2,ns
f_n0 = psi_0*epsilon_1*((0.d0,-1.d0)**n)*(dble(n)/a)*besjn(n,eta*a)
f_n1 = psi_0*epsilon_1*((0.d0,-1.d0)**n)*(eta/2.d0)*(besjn(n-1,eta*a)-besjn(n+1,eta*a))
m11 = (gamma/2.d0)*(besselh(n-1,1,gamma*a)-besselh(n+1,1,gamma*a))
m12 = -(dble(n)/a)*(besselh(int(n),1,eta*a))
m21 = (dble(n)/a)*besselh(int(n),1,gamma*a)
m22 = -(eta/2.d0)*(besselh(n-1,1,eta*a)-besselh(n+1,1,eta*a))
ABn(1) = (1.d0/(m11*m22 - m12*m21))*(m22*f_n0-m12*f_n1)
ABn(2) = (1.d0/(m11*m22 - m12*m21))*(-m21*f_n0+m11*f_n1)
q = q+(-psi_0*epsilon_1*((0.d0,-1.d0)**n)*(dble(n)/r)*besjn(n,eta*r)&
+ ABn(1)*(gamma/2.d0)*(besselh(n-1,1,gamma*r)&
-besselh(n+1,1,gamma*r)) &
- ABn(2)*(dble(n)/r)*besselh(int(n),1,eta*r))*sin(dble(n)*theta)
v = v + (-psi_0*epsilon_1*((0.d0,-1.d0)**n)*(eta/2.d0)*( besjn(n-1,eta*r) &
-besjn(n+1,eta*r)) &
+ ABn(1)*(dble(n)/r)*besselh(int(n),1,gamma*r) &
- ABn(2)*(eta/2.d0)*(besselh(n-1,1,eta*r)&
-besselh(n+1,1,eta*r)))*cos(dble(n)*theta)
end do
!$OMP END PARALLEL DO
! disp('DONE COMPUTING COEFICIENTS')
du = dreal(exp(omega*(0.d0,1.d0)*t)*(cos(theta)*q-sin(theta)*v))
dv = dreal(exp(omega*(0.d0,1.d0)*t)*(sin(theta)*q+cos(theta)*v))
! for GetDP: return imaginary part in dut, dvt
dut = dimag(exp(omega*(0.d0,1.d0)*t)*(cos(theta)*q-sin(theta)*v))
dvt = dimag(exp(omega*(0.d0,1.d0)*t)*(sin(theta)*q+cos(theta)*v))
! return real part in du,dv
du = dreal(exp(omega*(0.d0,1.d0)*t)*(cos(theta)*p-sin(theta)*q))
dv = dreal(exp(omega*(0.d0,1.d0)*t)*(sin(theta)*p+cos(theta)*q))
! return imaginary part in dut, dvt
dut = dimag(exp(omega*(0.d0,1.d0)*t)*(cos(theta)*p-sin(theta)*q))
dvt = dimag(exp(omega*(0.d0,1.d0)*t)*(sin(theta)*p+cos(theta)*q))
end subroutine cylindrical_wallsout
! cylindrical_cavity This program computes particular solutions to
! the elastic wave equation in cylindrical geometries,
! see: https://bitbucket.org/appelo/pewe
!
! Copyright (C) 2015 Kristoffer Virta & Daniel Appelo
!
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program. If not, see <http://www.gnu.org/licenses/>.