Commit 61332ccc authored by Serhiy Mochalskyy's avatar Serhiy Mochalskyy
Browse files

Delete svn folder

parent 55cee5b5
Pipeline #1723 skipped
# config.in file for the LINUX CLUSTER in Garching
include $(FILES_MK)
ifeq ($(DEBUG),yes)
FFLAGS = -autodouble -O0 -I$(OBJ_DIR) -module $(OBJ_DIR) -warn all -check all -traceback
else
FFLAGS = -autodouble -O3 -openmp -I$(OBJ_DIR) -module $(OBJ_DIR)
endif
FPPFLAGS =
LD = $(FC)
LDFLAGS = $(FFLAGS)
DPREF = -D
FPPFLAGS += $(DPREF)$(COMPILER) $(DPREF)$(OSTYPE)
########################################################################
MKLPATH = $(MKLROOT)/lib/intel64
LAPACK_LIB = -L$(MKLPATH) -lmkl_lapack95_lp64 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -liomp5
LIB_DIR =
LIB = $(LAPACK_LIB)
########################################################################
all: $(EXECOS_MAIN)
$(EXECOS_MAIN): $(MAIN_OBJ) $(OS_MK) $(FILES_MK)
$(LD) $(LDFLAGS) -o $@ $(MAIN_OBJ) $(addprefix -L,$(LIB_DIR)) $(LIB)
mv $@ $(OBJ_DIR)
$(OBJ_DIR)/%.o: %.f90
$(FC) -fpp $(FPPFLAGS) $(FFLAGS) -c -o $@ $<
! ----------------------------------------------------------------------
subroutine bfield_par(b_par,xp,yp,zp,xu,yu,zu,n,x,y,z,ntri,j_pot,n_pot)
! ----------------------------------------------------------------------
! 23.09.09
! purpose:
!
!
! ----------------------------------------------------------------------
implicit none
real, dimension (n) :: xp,yp,zp,xu,yu,zu
real, dimension (n,n_pot) :: b_par
real, dimension (ntri,3) :: x,y,z
integer,dimension (ntri,3) :: j_pot
integer :: n,ntri,n_pot
! ----------------------------------------------------------------------
! local variables
real :: x1,y1,z1,x2,y2,z2,x3,y3,z3,sn
real :: s1,s2,s3,d221,d232,d213,al1,al2,al3 &
,vx32,vy32,vz32,vx13,vy13,vz13 &
,ata1,ata2,ata3,v21,v32,v13,at &
,s21,s22,s23,dp1,dm1,dp2,dm2,dp3,dm3 &
,ap1,am1,ap2,am2,ap3,am3 &
,h,ln1,ln2,ln3,ar1,ar2,ar3,vx,vy,vz &
,x21,y21,z21,x32,y32,z32,x13,y13,z13 &
,tx1,ty1,tz1,tx2,ty2,tz2,tx3,ty3,tz3 &
,nx,ny,nz,pi41,area,d21,d32,d13 &
,jx1,jy1,jz1,jx2,jy2,jz2,jx3,jy3,jz3,temp3 &
,dep1,dep2,dep3,dem1,dem2,dem3,temp1,temp2
integer :: i,k,j1,j2,j3
! ----------------------------------------------------------------------
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE( &
!$OMP x1,y1,z1,x2,y2,z2,x3,y3,z3,sn &
!$OMP ,s1,s2,s3,d221,d232,d213,al1,al2,al3 &
!$OMP ,vx32,vy32,vz32,vx13,vy13,vz13 &
!$OMP ,ata1,ata2,ata3,v21,v32,v13,at &
!$OMP ,s21,s22,s23,dp1,dm1,dp2,dm2,dp3,dm3 &
!$OMP ,ap1,am1,ap2,am2,ap3,am3 &
!$OMP ,h,ln1,ln2,ln3,ar1,ar2,ar3,vx,vy,vz &
!$OMP ,x21,y21,z21,x32,y32,z32,x13,y13,z13 &
!$OMP ,tx1,ty1,tz1,tx2,ty2,tz2,tx3,ty3,tz3 &
!$OMP ,nx,ny,nz,pi41,area,d21,d32,d13 &
!$OMP ,jx1,jy1,jz1,jx2,jy2,jz2,jx3,jy3,jz3,temp3 &
!$OMP ,dep1,dep2,dep3,dem1,dem2,dem3,temp1,temp2 &
!$OMP ,i,k,j1,j2,j3 )
!$OMP DO
do k=1,ntri
x21 = x(k,2) - x(k,1)
y21 = y(k,2) - y(k,1)
z21 = z(k,2) - z(k,1)
x32 = x(k,3) - x(k,2)
y32 = y(k,3) - y(k,2)
z32 = z(k,3) - z(k,2)
x13 = x(k,1) - x(k,3)
y13 = y(k,1) - y(k,3)
z13 = z(k,1) - z(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)
pi41 = .125/asin(1.)
jx1 = x32*area*pi41
jy1 = y32*area*pi41
jz1 = z32*area*pi41
jx2 = x13*area*pi41
jy2 = y13*area*pi41
jz2 = z13*area*pi41
jx3 = x21*area*pi41
jy3 = y21*area*pi41
jz3 = z21*area*pi41
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)
j1 = j_pot(k,1)
j2 = j_pot(k,2)
j3 = j_pot(k,3)
do i=1,n
x1 = x(k,1) - xp(i)
y1 = y(k,1) - yp(i)
z1 = z(k,1) - zp(i)
x2 = x(k,2) - xp(i)
y2 = y(k,2) - yp(i)
z2 = z(k,2) - zp(i)
x3 = x(k,3) - xp(i)
y3 = y(k,3) - yp(i)
z3 = z(k,3) - zp(i)
sn = nx*x1+ny*y1+nz*z1
h = abs(sn)
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)/(s1+s2-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)
dm1 = dp1-d221
dm2 = dp2-d232
dm3 = dp3-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*dm1
dem1 = ar1**2+h*d221*(h+s1)
am2 = ar2*dm2
dem2 = ar2**2+h*d232*(h+s2)
am3 = ar3*dm3
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)
at = sign(1.,sn)*(ata1+ata2+ata3)
vx = -nx*at + al1*tx3/d21+al2*tx1/d32+al3*tx2/d13
vy = -ny*at + al1*ty3/d21+al2*ty1/d32+al3*ty2/d13
vz = -nz*at + al1*tz3/d21+al2*tz1/d32+al3*tz2/d13
temp1 = (vy*jz1-vz*jy1)*xu(i) &
+ (vz*jx1-vx*jz1)*yu(i) &
+ (vx*jy1-vy*jx1)*zu(i)
!$OMP ATOMIC
b_par(i,j1) = b_par(i,j1)+ temp1
temp2 = (vy*jz2-vz*jy2)*xu(i) &
+ (vz*jx2-vx*jz2)*yu(i) &
+ (vx*jy2-vy*jx2)*zu(i)
!$OMP ATOMIC
b_par(i,j2) = b_par(i,j2)+ temp2
temp3 = (vy*jz3-vz*jy3)*xu(i) &
+ (vz*jx3-vx*jz3)*yu(i) &
+ (vx*jy3-vy*jx3)*zu(i)
!$OMP ATOMIC
b_par(i,j3) = b_par(i,j3)+ temp3
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
end subroutine bfield_par
subroutine resistive_wall_response
use icontr
use contr_su
use solv
use tri_w
use tri_p
use coil2d
!-----------------------------------------------------------------------
implicit none
real, dimension(:,:),allocatable :: a_ye,a_ey,S_ww,d_ee
real, dimension( :),allocatable :: gamma
integer :: i,ier,j,k,nd_bez,n_w,lwork
integer :: n_tor_jorek, i_tor_jorek(999)
integer, allocatable :: ipiv(:)
real, allocatable :: s_ww_inv(:,:), xyzpot_w(:,:), work(:)
character(len=64) :: format_type
character(len=512) :: char512
!-----------------------------------------------------------------------
nd_bez = (nd_harm + nd_harm0)/2*n_dof_bnd
n_w = npot_w+ncoil
write(*,*) ' nd_harm :',nd_harm
write(*,*) ' nd_harm0 :',nd_harm0
write(*,*) ' n_bnd :',n_bnd
write(*,*) ' nd_bez :',nd_bez
allocate(gamma(n_w),S_ww(n_w,n_w),d_ee(nd_bez,nd_bez), &
a_ye(n_w,nd_bez),a_ey(nd_bez,n_w), stat=ier )
gamma =0. ; S_ww= 0.; a_ye=0.; a_ey=0.
call simil_trafo(a_ww,a_rw,n_w,gamma,S_ww)
a_ye =0.
do i=1,n_w
do k=1,nd_bez
do j=1,n_w
a_ye(i,k) = a_ye(i,k) + S_ww(j,i)*a_we(j,k)
enddo
enddo
enddo
do i=1,n_w
do k=1,nd_bez
a_ye(i,k) = a_ye(i,k)/gamma(i)
enddo
enddo
a_ey = 0.
do i=1,nd_bez
do k=1,n_w
do j=1,n_w
a_ey(i,k) = a_ey(i,k) + a_ew(i,j)*S_ww(j,k)
enddo
enddo
enddo
!------------------------------------------------------------------
d_ee = 0.
do i=1,nd_bez
do k=1,nd_bez
do j=1,n_w
d_ee(i,k)= d_ee(i,k) + a_ey(i,j)*a_ye(j,k)
enddo
enddo
enddo
do i=1,nd_bez
do k=1,nd_bez
write(123,'(2i6,1pe12.4)') i,k,d_ee(i,k)
enddo
enddo
!------------------------------------------------------------------
open(60,file='starwall_d_yy',form="formatted",iostat=ier)
write(60,'(i8)') n_w
do i=1,n_w
write(60,'(i8,1pe16.8)') i,1./gamma(i)
enddo
close(60)
open(60,file='starwall_m_ye',form="formatted",iostat=ier)
write(60,'(2i8)') n_w,nd_bez
do i=1,n_w
do k=1,nd_bez
write(60,'(2i8,1pe16.8)') i,k,a_ye(i,k)
enddo
enddo
close(60)
open(60,file='starwall_m_ey',form="formatted",iostat=ier)
write(60,'(2i8)') nd_bez,n_w
do i=1,nd_bez
do k=1,n_w
write(60,'(2i8,1pe16.8)') i,k,a_ey(i,k)
enddo
enddo
close(60)
open(60,file='starwall_m_ee',form="formatted",iostat=ier)
write(60,'(2i8)') nd_bez,nd_bez
do i=1,nd_bez
do k=1,nd_bez
write(60,'(2i8,1pe16.8)') i,k,a_ee(i,k)
enddo
enddo
close(60)
n_tor_jorek = 0
do i = 1, n_harm
if ( n_tor(i) == 0 ) then
n_tor_jorek = n_tor_jorek + 1
i_tor_jorek(n_tor_jorek) = 1
else
n_tor_jorek = n_tor_jorek + 2
i_tor_jorek(n_tor_jorek-1:n_tor_jorek) = n_tor(i)*2 + (/ 0, 1 /)
end if
end do
! --- Determine inverse of S_ww matrix
allocate( ipiv(n_w), s_ww_inv(n_w,n_w) )
s_ww_inv(:,:) = s_ww(:,:)
call dgetrf(n_w, n_w, s_ww_inv, n_w, ipiv, ier)
if ( ier /= 0 ) then
write(*,*) 'ERROR calling dgetrf when computing inverse of S_ww. info=', ier
stop
end if
lwork = 10*n_w; allocate( work(lwork) )
call dgetri(n_w, s_ww_inv, n_w, ipiv, work, lwork, ier)
if ( ier/= 0 ) then
write(*,*) 'ERROR calling dgetri when computing inverse of S_ww. info=', ier
stop
end if
deallocate( ipiv, work )
! --- Write out the STARWALL response as a single file
format_type = 'unformatted'
! (Determine positions of wall triangle nodes)
allocate( xyzpot_w(npot_w,3) )
do i = 1, ntri_w
do j = 1, 3
xyzpot_w(jpot_w(i,j),:) = (/ xw(i,j), yw(i,j), zw(i,j) /)
end do
end do
open(60, file='starwall-response.dat', form=format_type, status='replace', action='write')
if ( format_type == 'formatted' ) then
130 format(a)
131 format('#@intparam ',a24,99i12)
132 format('#@array ',a24,i12,a24,99i12)
133 format(4es24.16)
134 format(8i12)
write(60,130) '#@content STARWALL VACUUM RESPONSE DATA FOR THE JOREK CODE'
write(60,131) 'file_version', 1
write(60,131) 'n_bnd', N_bnd
write(60,131) 'nd_bez', nd_bez
write(60,131) 'ncoil', ncoil
write(60,131) 'npot_w', npot_w
write(60,131) 'n_w', n_w
write(60,131) 'ntri_w', ntri_w
write(60,131) 'n_tor', n_tor_jorek
write(60,132) 'i_tor', 1, 'int', n_tor_jorek, 0
write(60,134) i_tor_jorek(1:n_tor_jorek)
write(60,132) 'yy', 1, 'float', n_w, 0
write(60,133) 1./gamma(:)
write(60,132) 'ye', 2, 'float', n_w, nd_bez
write(60,133) a_ye(:,:)
write(60,132) 'ey', 2, 'float', nd_bez, n_w
write(60,133) a_ey(:,:)
write(60,132) 'ee', 2, 'float', nd_bez, nd_bez
write(60,133) a_ee(:,:)
write(60,132) 's_ww', 2, 'float', n_w, n_w
write(60,133) S_ww(:,:)
write(60,132) 's_ww_inv', 2, 'float', n_w, n_w
write(60,133) s_ww_inv(:,:)
write(60,132) 'xyzpot_w', 2, 'float', npot_w, 3
write(60,133) xyzpot_w(:,:)
write(60,132) 'jpot_w', 2, 'int', ntri_w, 3
write(60,134) jpot_w(:,:)
else
char512='#@content STARWALL VACUUM RESPONSE DATA FOR THE JOREK CODE'
write(60) 42, 42.d0 !### for an elementary check in JOREK
write(60) char512
write(60) '#@intparam ', 'file_version ', 1
write(60) '#@intparam ', 'n_bnd ', N_bnd
write(60) '#@intparam ', 'nd_bez ', nd_bez
write(60) '#@intparam ', 'ncoil ', ncoil
write(60) '#@intparam ', 'npot_w ', npot_w
write(60) '#@intparam ', 'n_w ', n_w
write(60) '#@intparam ', 'ntri_w ', ntri_w
write(60) '#@intparam ', 'n_tor ', n_tor_jorek
write(60) '#@array ', 'i_tor ', 1, 'int ', n_tor_jorek, 0
write(60) i_tor_jorek(1:n_tor_jorek)
write(60) '#@array ', 'yy ', 1, 'float ', n_w, 0
write(60) 1.d0/gamma(:)
write(60) '#@array ', 'ye ', 2, 'float ', n_w, nd_bez
write(60) a_ye(:,:)
write(60) '#@array ', 'ey ', 2, 'float ', nd_bez, n_w
write(60) a_ey(:,:)
write(60) '#@array ', 'ee ', 2, 'float ', nd_bez, nd_bez
write(60) a_ee(:,:)
write(60) '#@array ', 's_ww ', 2, 'float ', n_w, n_w
write(60) S_ww(:,:)
write(60) '#@array ', 's_ww_inv ', 2, 'float ', n_w, n_w
write(60) s_ww_inv(:,:)
write(60) '#@array ', 'xyzpot_w ', 2, 'float ', npot_w, 3
write(60) xyzpot_w(:,:)
write(60) '#@array ', 'jpot_w ', 2, 'int ', ntri_w, 3
write(60) jpot_w(:,:)
end if
close(60)
! --- VTK: Control surface
140 format(a)
141 format(a,i8,a)
142 format(3es24.16)
143 format(a,2i8)
144 format(4i8)
open(60, file='control.vtk')
write(60,140) '# vtk DataFile Version 2.0'
write(60,140) 'testdata'
write(60,140) 'ASCII'
write(60,140) 'DATASET POLYDATA'
write(60,141) 'POINTS', 3*ntri_p, ' float'
do i = 1, ntri_p
do j = 1, 3
write(60,142) xp(i,j), yp(i,j), zp(i,j)
end do
end do
write(60,143) 'POLYGONS', ntri_p, ntri_p*4
do i = 1, ntri_p
write(60,144) 3, 3*(i-1), 3*(i-1)+1, 3*(i-1)+2
end do
write(60,141) 'POINT_DATA', 3*ntri_p
write(60,140) 'SCALARS potentials float'
write(60,140) 'LOOKUP_TABLE default'
do i = 1, 3*ntri_p
write(60,144) mod(i/3,5) !###
end do
close(60)
! --- VTK: Wall
open(60, file='wall.vtk')
write(60,140) '# vtk DataFile Version 2.0'
write(60,140) 'testdata'
write(60,140) 'ASCII'
write(60,140) 'DATASET POLYDATA'
write(60,141) 'POINTS', npot_w, ' float'
do i = 1, npot_w
write(60,142) xyzpot_w(i,:)
end do
write(60,143) 'POLYGONS', ntri_w, ntri_w*4
do i = 1, ntri_w
write(60,144) 3, jpot_w(i,:)-1
end do
write(60,141) 'POINT_DATA', npot_w
write(60,140) 'SCALARS potentials float'
write(60,140) 'LOOKUP_TABLE default'
do i = 1, npot_w
write(60,144) mod(i,3) !###
end do
close(60)
! --- VTK: Coils
open(60, file='coils.vtk')
write(60,140) '# vtk DataFile Version 2.0'
write(60,140) 'testdata'
write(60,140) 'ASCII'
write(60,140) 'DATASET POLYDATA'
write(60,141) 'POINTS', 3*ntri_c, ' float'
do i = 1, ntri_c
write(60,142) x_coil(i,1), z_coil(i,1), y_coil(i,1)
write(60,142) x_coil(i,2), z_coil(i,2), y_coil(i,2)
write(60,142) x_coil(i,3), z_coil(i,3), y_coil(i,3)
end do
write(60,143) 'POLYGONS', ntri_c, ntri_c*4
do i = 1, ntri_c
write(60,144) 3, 3 * i - (/ 3, 2, 1 /)
end do
write(60,141) 'POINT_DATA', 3*ntri_c
write(60,140) 'SCALARS potentials float'
write(60,140) 'LOOKUP_TABLE default'
do i = 1, ntri_c
write(60,142) 0.d0
end do
close(60)
deallocate( xyzpot_w )
end subroutine resistive_wall_response
subroutine solver
use icontr
use solv
use contr_su
use tri_w
use tri_p
use coil2d
!-----------------------------------------------------------------------
implicit none
integer :: i,j,k,nd_w,nd_we,nd_bez,info
!-----------------------------------------------------------------------
write(*,*) ' SOLVER: n_dof_bnd : ',n_dof_bnd
nd_bez = (nd_harm+nd_harm0)/2 * n_dof_bnd
nd_w = ncoil + npot_w
nd_we = nd_w + nd_bez
write(*,*) ' SOLVER: nd_bez : ',nd_bez
!-----------------------------------------------------------------------
!!!!!call perfon('ma_pp')
call matrix_pp
!!!!!call perfoff
!!!!!call perfon('ma_wp')
call matrix_wp
!!!!!call perfoff
!!!!!call perfon('ma_ww')
call matrix_ww
!!!!!call perfoff
!!!!!call perfon('ma_rw')
call matrix_rw
!!!!!call perfoff
if (nu_coil.ne.0) then
!!!!!call perfon('coil')
call matrix_cc
call matrix_cp
call matrix_wc
call matrix_rc
!!!!!call perfoff
endif
!!!!!call perfon('ma_pe')
call matrix_pe
!!!!!call perfoff
!!!!!call perfon('ma_ep')
call matrix_ep
!!!!!call perfoff
!!!!!call perfon('ma_ew')
call matrix_ew
!!!!!call perfoff
if (nu_coil.ne.0) then
!!!!!call perfon('ma_ec')
call matrix_ec
!!!!!call perfoff
endif
write(*,'(A,i5)') 'SOLVER : npot_p : ',npot_p
write(*,'(A,i5)') 'SOLVER : nd_w : ',nd_w
write(*,'(A,2i5,A,2i5)') 'SOLVER : copying a_wp ', nd_w, npot_p,' to a_pwe ',npot_p,nd_w
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,k)
!$OMP DO
do i=1,npot_p
do k=1,nd_w
a_pwe(i,k) = a_wp(k,i)
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
!!!!!call perfon('dpof')
call dpotrf('U',npot_p,a_pp,npot_p,info)
call dpotrs('U',npot_p,nd_we,a_pp,npot_p,a_pwe,npot_p,info)
!!!!!call perfoff
deallocate(a_pp)
!-----------------------------------------------------------------------
! compute parallel component of magnetic field
allocate(a_we(nd_w,nd_bez),a_ee(nd_bez,nd_bez))
a_we = 0.; a_ee = 0.
! compute a_ee = a_ep*a_pp^(-1)*a_pe a_ee = vacuum_response
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k)
!$OMP DO
do i=1,nd_bez
do k=1,nd_bez
do j=1,npot_p
a_ee(i,k) = a_ee(i,k)+a_ep(i,j)*a_pwe(j,k+nd_w)
enddo
enddo
enddo
!$OMP END DO
! compute a_ew = a_ew - a_ep*a_pp^(-1)*a_pw
!$OMP DO
do i=1,nd_bez
do k=1,nd_w
do j=1,npot_p