Commit 4e75df8b authored by Serhiy Mochalskyy's avatar Serhiy Mochalskyy
Browse files

cleaning

parent e8df7704
Pipeline #2433 skipped
......@@ -240,6 +240,7 @@ do i=1,ntri_w
enddo
do iw =1,ntri_w
do ip=1,ntri_p
do kw =1,3
......@@ -267,7 +268,6 @@ do iw =1,ntri_w
enddo
enddo
enddo
!===================================================================
call ScaLAPACK_mapping_i(1+ncoil,i_loc,inside_i)
......
......@@ -53,6 +53,40 @@ subroutine ScaLAPACK_mapping_i(i,i_loc,inside_row)
end subroutine ScaLAPACK_mapping_i
!DIR$ ATTRIBUTES FORCEINLINE :: ScaLAPACK_mapping_j2
subroutine ScaLAPACK_mapping_j2(j,j_loc,inside_col)
! ----------------------------------------------------------------------
! purpose: 20/08/2015
! Check the mapping "j" index for ScaLAPCK local matrix
!! ----------------------------------------------------------------------
use mpi_v ! new module with mpi variables
use sca ! new module for ScaLAPACK variables
! ----------------------------------------------------------------------
implicit none
integer INDXG2L,INDXG2P
EXTERNAL INDXG2L,INDXG2P
integer j,j_loc
logical inside_col
inside_col = .true.
ipc = INDXG2P(j,NB,0,0,NPCOL)
if (ipc.eq.MYCOL) then
inside_col = .false.
j_loc = INDXG2L(j,NB,0,0,NPCOL)
endif
end subroutine ScaLAPACK_mapping_j2
!DIR$ ATTRIBUTES FORCEINLINE :: ScaLAPACK_mapping_j
subroutine ScaLAPACK_mapping_j(j,j_loc,inside_col)
! ----------------------------------------------------------------------
......
......@@ -194,9 +194,6 @@ if(rank==0)write(210+rank,*) time12-time11
call a_we_computing
!==================================================================
deallocate(a_pwe_loc_s)
call MPI_BARRIER(MPI_COMM_WORLD,ier)
time13=MPI_WTIME()
if(rank==0)write(220+rank,*) time13-time12
......@@ -209,7 +206,6 @@ call matrix_multiplication
deallocate(a_wp_loc,a_pwe_loc,a_ep_loc_sca)
call MPI_BARRIER(MPI_COMM_WORLD,ier)
time14=MPI_WTIME()
if(rank==0)write(230+rank,*) time14-time13
......
......@@ -330,3 +330,195 @@ subroutine tri_induct_3_2_b(ntri_p,ntri_w,k,i,xw,yw,zw,dima_sca)
end subroutine tri_induct_3_2_b
subroutine tri_induct(xp,yp,zp,ntri_p,xw,yw,zw,ntri_w,dima)
! ----------------------------------------------------------------------
! 27.06.01
! purpose:
!! ----------------------------------------------------------------------
implicit none
real,dimension (ntri_p,3) :: xp,yp,zp
real,dimension (ntri_w,3) :: xw,yw,zw
real,dimension (ntri_p,ntri_w) :: dima
integer :: ntri_p, ntri_w
! ----------------------------------------------------------------------
! local variables
real,dimension (7*ntri_p) :: xm,ym,zm
real,dimension (7) :: w
real :: pi41,sq15, f1, f2, f3, f4, f5
real :: x31, y31,z31, x21, y21, z21
real :: xp21,yp21,zp21 &
,x1,y1,z1,x2,y2,z2,x3,y3,z3,s1,s2,s3 &
,d221,d232,d213,al1,al2,al3 &
,d21,d32,d13,nx,ny,nz,area &
,tx1,ty1,tz1,tx2,ty2,tz2,tx3,ty3,tz3 &
,h,s21,s22,s23,dp1,dp2,dp3 &
,ap1,am1,ap2,am2,ap3,am3,ata1,ata2,ata3 &
,at,a1,a2,a3,ar1,ar2,ar3,ve &
,x32,y32,z32,x13,y13,z13 &
,dep1,dep2,dep3,dem1,dem2,dem3
integer :: i,i1,k,ik,l
!-----------------------------------------------------------------
pi41 = .125/asin(1.)
sq15 = sqrt(15.)
f1 = 1./3.
f2 = (6.+ sq15)/21.
f3 = (9.- 2.*sq15)/21.
f4 = (9.+ 2.*sq15)/21.
f5 = (6.- sq15)/21.
w(1) = .1125
w(2) = (155.+sq15)/2400.
w(3) = (155.+sq15)/2400.
w(4) = (155.+sq15)/2400.
w(5) = (155.-sq15)/2400.
w(6) = (155.-sq15)/2400.
w(7) = (155.-sq15)/2400.
! ----------------------------------------------------------------------
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE( &
!$OMP x31, y31,z31, x21, y21, z21 &
!$OMP ,xp31,yp31,zp31,xp21,yp21,zp21 &
!$OMP ,x1,y1,z1,x2,y2,z2,x3,y3,z3,s1,s2,s3 &
!$OMP ,d221,d232,d213,al1,al2,al3 &
!$OMP ,d21,d32,d13,nx,ny,nz,area &
!$OMP ,tx1,ty1,tz1,tx2,ty2,tz2,tx3,ty3,tz3 &
!$OMP ,h,s21,s22,s23,dp1,dp2,dp3 &
!$OMP ,ap1,am1,ap2,am2,ap3,am3,ata1,ata2,ata3 &
!$OMP ,at,a1,a2,a3,ar1,ar2,ar3,ve &
!$OMP ,x32,y32,z32,x13,y13,z13 &
!$OMP ,dep1,dep2,dep3,dem1,dem2,dem3 &
!$OMP ,i,i1,k,ik,l )
!$OMP DO
do i=1,ntri_p
i1 = 7*(i-1)
x21 = xp(i,2) - xp(i,1)
y21 = yp(i,2) - yp(i,1)
z21 = zp(i,2) - zp(i,1)
x31 = xp(i,3) - xp(i,1)
y31 = yp(i,3) - yp(i,1)
z31 = zp(i,3) - zp(i,1)
xm(1+i1) = xp(i,1)+f1*(x21+x31)
ym(1+i1) = yp(i,1)+f1*(y21+y31)
zm(1+i1) = zp(i,1)+f1*(z21+z31)
xm(2+i1) = xp(i,1)+f2*(x21+x31)
ym(2+i1) = yp(i,1)+f2*(y21+y31)
zm(2+i1) = zp(i,1)+f2*(z21+z31)
xm(5+i1) = xp(i,1)+f5*(x21+x31)
ym(5+i1) = yp(i,1)+f5*(y21+y31)
zm(5+i1) = zp(i,1)+f5*(z21+z31)
xm(3+i1) = xp(i,1)+f3*x21+f2*x31
ym(3+i1) = yp(i,1)+f3*y21+f2*y31
zm(3+i1) = zp(i,1)+f3*z21+f2*z31
xm(4+i1) = xp(i,1)+f2*x21+f3*x31
ym(4+i1) = yp(i,1)+f2*y21+f3*y31
zm(4+i1) = zp(i,1)+f2*z21+f3*z31
xm(6+i1) = xp(i,1)+f4*x21+f5*x31
ym(6+i1) = yp(i,1)+f4*y21+f5*y31
zm(6+i1) = zp(i,1)+f4*z21+f5*z31
xm(7+i1) = xp(i,1)+f5*x21+f4*x31
ym(7+i1) = yp(i,1)+f5*y21+f4*y31
zm(7+i1) = zp(i,1)+f5*z21+f4*z31
enddo
!$OMP END DO
!$OMP DO
do k=1,ntri_w
x21 = xw(k,2) - xw(k,1)
y21 = yw(k,2) - yw(k,1)
z21 = zw(k,2) - zw(k,1)
x32 = xw(k,3) - xw(k,2)
y32 = yw(k,3) - yw(k,2)
z32 = zw(k,3) - zw(k,2)
x13 = xw(k,1) - xw(k,3)
y13 = yw(k,1) - yw(k,3)
z13 = zw(k,1) - zw(k,3)
d221 = x21**2+y21**2+z21**2
d232 = x32**2+y32**2+z32**2
d213 = x13**2+y13**2+z13**2
d21 = sqrt(d221)
d32 = sqrt(d232)
d13 = sqrt(d213)
nx = -y21*z13 + z21*y13
ny = -z21*x13 + x21*z13
nz = -x21*y13 + y21*x13
area = 1./sqrt(nx*nx+ny*ny+nz*nz)
nx = nx*area
ny = ny*area
nz = nz*area
tx1 = y32*nz-z32*ny
ty1 = z32*nx-x32*nz
tz1 = x32*ny-y32*nx
tx2 = y13*nz-z13*ny
ty2 = z13*nx-x13*nz
tz2 = x13*ny-y13*nx
tx3 = y21*nz-z21*ny
ty3 = z21*nx-x21*nz
tz3 = x21*ny-y21*nx
! ----------------------------------------------------------------------
do l=1,7
do i=1,ntri_p
ik = l+7*(i-1)
x1 = xw(k,1) - xm(ik)
y1 = yw(k,1) - ym(ik)
z1 = zw(k,1) - zm(ik)
x2 = xw(k,2) - xm(ik)
y2 = yw(k,2) - ym(ik)
z2 = zw(k,2) - zm(ik)
x3 = xw(k,3) - xm(ik)
y3 = yw(k,3) - ym(ik)
z3 = zw(k,3) - zm(ik)
! ----------------------------------------------------------------------
h = abs(nx*x1+ny*y1+nz*z1)
s21 = x1**2+y1**2+z1**2
s22 = x2**2+y2**2+z2**2
s23 = x3**2+y3**2+z3**2
s1 = sqrt(s21)
s2 = sqrt(s22)
s3 = sqrt(s23)
al1 = alog((s2+s1+d21)/(s2+s1-d21))
al2 = alog((s3+s2+d32)/(s3+s2-d32))
al3 = alog((s1+s3+d13)/(s1+s3-d13))
ar1 = x1*tx3+y1*ty3+z1*tz3
ar2 = x2*tx1+y2*ty1+z2*tz1
ar3 = x3*tx2+y3*ty2+z3*tz2
dp1 = .5*(s22-s21+d221)
dp2 = .5*(s23-s22+d232)
dp3 = .5*(s21-s23+d213)
ap1 = ar1*dp1
dep1 = ar1**2+h*d221*(h+s2)
ap2 = ar2*dp2
dep2 = ar2**2+h*d232*(h+s3)
ap3 = ar3*dp3
dep3 = ar3**2+h*d213*(h+s1)
am1 = ar1*(dp1-d221)
dem1 = ar1**2+h*d221*(h+s1)
am2 = ar2*(dp2-d232)
dem2 = ar2**2+h*d232*(h+s2)
am3 = ar3*(dp3-d213)
dem3 = ar3**2+h*d213*(h+s3)
ata1 = atan2(ap1*dem1-am1*dep1,dep1*dem1+ap1*am1)
ata2 = atan2(ap2*dem2-am2*dep2,dep2*dem2+ap2*am2)
ata3 = atan2(ap3*dem3-am3*dep3,dep3*dem3+ap3*am3)
ve = pi41*area*(-h*(ata1+ata2+ata3) &
+ ar1*al1/d21+ar2*al2/d32+ar3*al3/d13)
dima(i,k)=dima(i,k) +w(l)*ve
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
end subroutine tri_induct
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment