diff --git a/src_3d/a_ye_computing.f90 b/src_3d/a_ye_computing.f90 index 8aebf0bd5bad0f349492e06004ce477ca6bc4460..f38cbad5bf565e536a96c60a0d8c05f4dcdcef98 100644 --- a/src_3d/a_ye_computing.f90 +++ b/src_3d/a_ye_computing.f90 @@ -38,7 +38,7 @@ if(INFO_B .NE. 0) then endif -CALL DESCINIT(DESCC, n_w,nd_bez, NB, NB, 0, 0, CONTEXT, LDA_ye, INFO_C ) +CALL DESCINIT(DESC_ye, n_w,nd_bez, NB, NB, 0, 0, CONTEXT, LDA_ye, INFO_C ) if(INFO_C .NE. 0) then write(6,*) "Something is wrong in a_ye_computing CALL DESCINIT DESCC, INFO_C=",INFO_C stop @@ -46,7 +46,7 @@ endif CALL PDGEMM('T','N', n_w ,nd_bez, n_w , 1.0D0 , s_ww_loc ,1,1,DESCA, a_we_loc,1,1,DESCB,& - 1.0D0, a_ye_loc,1 , 1 , DESCC) + 1.0D0, a_ye_loc,1 , 1 , DESC_ye) do i=1,n_w diff --git a/src_3d/input.f90 b/src_3d/input.f90 index c4a6ae25836cc38449f058087c8f217a893d7883..88ae033c68da1b6f3723dc1d23f1226930d6c715 100644 --- a/src_3d/input.f90 +++ b/src_3d/input.f90 @@ -80,7 +80,7 @@ subroutine input call MPI_BCAST(NB, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ier) call MPI_BCAST(ORFAC, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier) call MPI_BCAST(lwork_cooficient, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ier) - !=================================================================================== + !================================================================================== if(rank==0) then write(outp,*) ' ScaLAPACK Input parameters:' write(outp,'(A,I4)') ' NB =', NB @@ -88,14 +88,16 @@ subroutine input write(outp,'(A,I4)') endif - !Only one task read file type because only one tast will write afterwards data to the file. if(rank==0) read(inp, params_output) ! Read namelist from STDIN + !================================================================================== + call MPI_BCAST(format_type, len(format_type), MPI_CHARACTER,0, MPI_COMM_WORLD, ier) + !================================================================================== if(rank==0) then write(outp,'(A,I4)') write(outp,*) 'Output starwall-response.dat:' write(outp,'(A,A)') ' file format = ',format_type endif - + ! --- Check input parameters if ( ( i_response < 0 ) .or. ( i_response > 2 ) ) then diff --git a/src_3d/mod_sca.f90 b/src_3d/mod_sca.f90 index b74707453de18c9e3bc58b425553beb8f949fa56..a010223ec494ae80943267066ce841528f783eb5 100644 --- a/src_3d/mod_sca.f90 +++ b/src_3d/mod_sca.f90 @@ -8,6 +8,7 @@ module sca LDA_we,LDA_ee,LDA_ww,LDA_rw,LDA_sww, LDA_s_ww_inv integer :: LDA_ey,LDA_ye,LDA_dee integer :: DESCA(9),DESCB(9),DESCZ(9),DESCC(9) + integer :: DESC_ye(9) integer :: INFO,INFO_A,INFO_B,INFO_Z,INFO_C real :: ORFAC integer :: lwork_cooficient diff --git a/src_3d/output.f90 b/src_3d/output.f90 index 5531e61e160dd0ef2f0fc414887594d5bc2a0892..a13ad957bf30215fd38f3b8221bf59bdfe6b89b6 100644 --- a/src_3d/output.f90 +++ b/src_3d/output.f90 @@ -17,6 +17,7 @@ if(rank==0) write(*,*) 'output begings' !====================================== !call print_starwall_m_ey !====================================== +!call print_starwall_m_ey2 !====================================== !call print_starwall_m_ee @@ -27,7 +28,8 @@ if(rank==0) write(*,*) 'output begings' !====================================== !====================================== -call print_starwall_response +!call print_starwall_response +call print_starwall_response2 !====================================== !====================================== diff --git a/src_3d/print_starwall_m_ey.f90 b/src_3d/print_starwall_m_ey.f90 index bf536c2f3368a00aa49a96e4b3dc00453f1c827b..7b5a89f27635d0c2eff23a73bc29f967a92121ba 100644 --- a/src_3d/print_starwall_m_ey.f90 +++ b/src_3d/print_starwall_m_ey.f90 @@ -34,3 +34,54 @@ real :: num if(rank==0) close(60) end subroutine print_starwall_m_ey + + + +subroutine print_starwall_m_ey2 + +use mpi_v +use resistive +use solv +use sca +!----------------------------------------------------------------------- +implicit none +include "mpif.h" + +integer :: i_loc,j_loc, IC,JC +real :: num +integer INDXL2G +EXTERNAL INDXL2G + +!----------------------------------------------------------------------- + call SCA_GRID(nd_bez,n_w) + + CALL DESCINIT(DESCA, nd_bez, n_w, NB, NB, 0, 0, CONTEXT, LDA_ey, INFO_A ) + if(INFO_A .NE. 0) then + write(6,*) "Something is wrong in print_starwall_m_ye CALL DESCINIT DESCA, INFO_A=",INFO_A + stop + endif + +! if(rank==0) then + open(60,file='starwall_m_ey2',form="formatted",iostat=ier) + ! write(60,'(2i8)') nd_bez,n_w + + ! endif + + + DO i_loc = 1,MP_A + + IC= INDXL2G( i_loc, NB, MYROW, 0, NPROW) + DO j_loc = 1,NQ_A + JC= INDXL2G( j_loc, NB, MYCOL, 0, NPCOL) + write(60,'(2i8,1pe16.8,1i8)') IC,JC,a_ey_loc(i_loc,j_loc), rank + END DO + END DO + +!if(rank==0) close(60) +close(60) + + + +end subroutine print_starwall_m_ey2 + + diff --git a/src_3d/print_starwall_response.f90 b/src_3d/print_starwall_response.f90 index f2da79833e502b1dc993f85e9f6ba945cdce5859..dffd93e0b75af942620269128c1f5463fee638c4 100644 --- a/src_3d/print_starwall_response.f90 +++ b/src_3d/print_starwall_response.f90 @@ -364,3 +364,434 @@ end if !if (format_type == 'formatted' ) then close(60) end subroutine print_starwall_response + + + + +subroutine print_starwall_response2 + +use icontr +use contr_su +use solv +use tri_w +use tri_p +use coil2d + +use sca +use mpi_v +use time +use resistive +!----------------------------------------------------------------------- +implicit none +include "mpif.h" + + +integer :: i,j,k,i_loc,j_loc +real :: num=0. +real,dimension(:,:),allocatable :: a_ye_print(:,:),a_ey_print(:,:),a_ee_print(:,:), & + s_ww_print(:,:),s_ww_inv_print(:,:) +integer :: IC, JC +integer INDXL2G +EXTERNAL INDXL2G + +!----------------------------------------------------------------------- + +! (Determine positions of wall triangle nodes) +if(rank==0) then + 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 +endif + +if (format_type == 'formatted' ) then + if(rank==0) then + open(60, file='starwall-response2.dat', form=format_type, status='replace', action='write') + 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 + + endif + + allocate(a_ye_print(n_w,nd_bez), stat=ier) + IF (IER /= 0) THEN + WRITE (*,*) "output.f90, can not allocate a_ye_print MYPROC_NUM=",MYPNUM + STOP + END IF + a_ye_print=0. + + call SCA_GRID(n_w,nd_bez) + DO i_loc = 1,MP_A + IC= INDXL2G( i_loc, NB, MYROW, 0, NPROW) + DO j_loc = 1,NQ_A + JC= INDXL2G( j_loc, NB, MYCOL, 0, NPCOL) + a_ye_print(IC,JC) = a_ye_loc(i_loc,j_loc) + END DO + END DO + + + if(rank==0) then + call MPI_REDUCE(MPI_IN_PLACE, a_ye_print, n_w*nd_bez, & + MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, IER) + else + call MPI_REDUCE(a_ye_print, a_ye_print, n_w*nd_bez, & + MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, IER) + endif + + + if(rank==0) write(60,133) a_ye_print(:,:) + deallocate (a_ye_print, a_ye_loc) + +!========================================================================================= + + allocate(a_ey_print(nd_bez,n_w),stat=ier) + IF (IER /= 0) THEN + WRITE (*,*) "output.f90, can not allocate a_ey_print MYPROC_NUM=",MYPNUM + STOP + END IF + a_ey_print=0. + + call SCA_GRID(nd_bez,n_w) + DO i_loc = 1,MP_A + IC= INDXL2G( i_loc, NB, MYROW, 0, NPROW) + DO j_loc = 1,NQ_A + JC= INDXL2G( j_loc, NB, MYCOL, 0, NPCOL) + a_ey_print(IC,JC) = a_ey_loc(i_loc,j_loc) + END DO + END DO + + if(rank==0) then + call MPI_REDUCE(MPI_IN_PLACE, a_ey_print, nd_bez*n_w, & + MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, IER) + else + call MPI_REDUCE(a_ey_print, a_ey_print, nd_bez*n_w, & + MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, IER) + + endif + + if(rank==0) then + write(60,132) 'ey', 2, 'float', nd_bez, n_w + write(60,133) a_ey_print(:,:) + endif + deallocate (a_ey_print, a_ey_loc) + +!========================================================================================= + + allocate(a_ee_print(nd_bez,nd_bez),stat=ier) + IF (IER /= 0) THEN + WRITE (*,*) "output.f90, can not allocate a_ee_print MYPROC_NUM=",MYPNUM + STOP + END IF + a_ee_print=0. + call SCA_GRID(nd_bez,nd_bez) + DO i_loc = 1,MP_A + IC= INDXL2G( i_loc, NB, MYROW, 0, NPROW) + DO j_loc = 1,NQ_A + JC= INDXL2G( j_loc, NB, MYCOL, 0, NPCOL) + a_ee_print(IC,JC) = a_ee_loc(i_loc,j_loc) + END DO + END DO + + if(rank==0) then + call MPI_REDUCE(MPI_IN_PLACE, a_ee_print, nd_bez*nd_bez, & + MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, IER) + else + call MPI_REDUCE(a_ee_print, a_ee_print, nd_bez*nd_bez, & + MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, IER) + endif + + if(rank==0) then + write(60,132) 'ee', 2, 'float', nd_bez, nd_bez + write(60,133) a_ee_print(:,:) + endif + deallocate (a_ee_print, a_ee_loc) + +!========================================================================================= + + allocate(s_ww_print(n_w,n_w),stat=ier) + IF (IER /= 0) THEN + WRITE (*,*) "output.f90, can not allocate s_ww_print MYPROC_NUM=",MYPNUM + STOP + END IF + s_ww_print=0. + + call SCA_GRID(n_w,n_w) + DO i_loc = 1,MP_A + IC= INDXL2G( i_loc, NB, MYROW, 0, NPROW) + DO j_loc = 1,NQ_A + JC= INDXL2G( j_loc, NB, MYCOL, 0, NPCOL) + s_ww_print(IC,JC) = s_ww_loc(i_loc,j_loc) + END DO + END DO + + if(rank==0) then + call MPI_REDUCE(MPI_IN_PLACE, s_ww_print, n_w*n_w, & + MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, IER) + else + call MPI_REDUCE(s_ww_print, s_ww_print, n_w*n_w, & + MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, IER) + endif + + if(rank==0) then + write(60,132) 's_ww', 2, 'float', n_w, n_w + write(60,133) s_ww_print(:,:) + endif + deallocate (s_ww_print, s_ww_loc) +!========================================================================================= + allocate(s_ww_inv_print(n_w,n_w),stat=ier) + IF (IER /= 0) THEN + WRITE (*,*) "output.f90, can not allocate s_ww_inv_print MYPROC_NUM=",MYPNUM + STOP + END IF + s_ww_inv_print=0. + call SCA_GRID(n_w,n_w) + DO i_loc = 1,MP_A + IC= INDXL2G( i_loc, NB, MYROW, 0, NPROW) + DO j_loc = 1,NQ_A + JC= INDXL2G( j_loc, NB, MYCOL, 0, NPCOL) + s_ww_inv_print(IC,JC) = s_ww_inv_loc(i_loc,j_loc) + END DO + END DO + + if(rank==0) then + call MPI_REDUCE(MPI_IN_PLACE, s_ww_inv_print, n_w*n_w, & + MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, IER) + else + call MPI_REDUCE(s_ww_inv_print, s_ww_inv_print, n_w*n_w, & + MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, IER) + endif + + if(rank==0) then + write(60,132) 's_ww_inv', 2, 'float', n_w, n_w + write(60,133) s_ww_inv_print(:,:) + endif + deallocate (s_ww_inv_print, s_ww_inv_loc) + + +!========================================================================================= + if(rank==0) then + 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(:,:) + + deallocate(xyzpot_w) + + endif + deallocate(jpot_w) +!========================================================================================= +else !if (format_type == 'unformatted' ) then + + if(rank==0) then + + open(60, file='starwall-response2.dat', form=format_type, status='replace', action='write') + + 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 + + endif + + allocate(a_ye_print(n_w,nd_bez), stat=ier) + IF (IER /= 0) THEN + WRITE (*,*) "output.f90, can not allocate a_ye_print MYPROC_NUM=",MYPNUM + STOP + END IF + a_ye_print=0. + + call SCA_GRID(n_w,nd_bez) + DO i_loc = 1,MP_A + IC= INDXL2G( i_loc, NB, MYROW, 0, NPROW) + DO j_loc = 1,NQ_A + JC= INDXL2G( j_loc, NB, MYCOL, 0, NPCOL) + a_ye_print(IC,JC) = a_ye_loc(i_loc,j_loc) + END DO + END DO + + if(rank==0) then + call MPI_REDUCE(MPI_IN_PLACE, a_ye_print, n_w*nd_bez, & + MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, IER) + else + call MPI_REDUCE(a_ye_print, a_ye_print, n_w*nd_bez, & + MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, IER) + endif + + + if(rank==0) write(60) a_ye_print(:,:) + deallocate (a_ye_print, a_ye_loc) + +!========================================================================================= + + allocate(a_ey_print(nd_bez,n_w),stat=ier) + IF (IER /= 0) THEN + WRITE (*,*) "output.f90, can not allocate a_ey_print MYPROC_NUM=",MYPNUM + STOP + END IF + a_ey_print=0. + + call SCA_GRID(nd_bez,n_w) + DO i_loc = 1,MP_A + IC= INDXL2G( i_loc, NB, MYROW, 0, NPROW) + DO j_loc = 1,NQ_A + JC= INDXL2G( j_loc, NB, MYCOL, 0, NPCOL) + a_ey_print(IC,JC) = a_ey_loc(i_loc,j_loc) + END DO + END DO + + if(rank==0) then + call MPI_REDUCE(MPI_IN_PLACE, a_ey_print, nd_bez*n_w, & + MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, IER) + else + call MPI_REDUCE(a_ey_print, a_ey_print, nd_bez*n_w, & + MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, IER) + + endif + + if(rank==0) then + write(60) '#@array ', 'ey ', 2, 'float ', nd_bez, n_w + write(60) a_ey_print(:,:) + endif + deallocate (a_ey_print, a_ey_loc) + +!========================================================================================= + + allocate(a_ee_print(nd_bez,nd_bez),stat=ier) + IF (IER /= 0) THEN + WRITE (*,*) "output.f90, can not allocate a_ee_print MYPROC_NUM=",MYPNUM + STOP + END IF + a_ee_print=0. + call SCA_GRID(nd_bez,nd_bez) + DO i_loc = 1,MP_A + IC= INDXL2G( i_loc, NB, MYROW, 0, NPROW) + DO j_loc = 1,NQ_A + JC= INDXL2G( j_loc, NB, MYCOL, 0, NPCOL) + a_ee_print(IC,JC) = a_ee_loc(i_loc,j_loc) + END DO + END DO + + if(rank==0) then + call MPI_REDUCE(MPI_IN_PLACE, a_ee_print, nd_bez*nd_bez, & + MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, IER) + else + call MPI_REDUCE(a_ee_print, a_ee_print, nd_bez*nd_bez, & + MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, IER) + endif + + if(rank==0) then + write(60) '#@array ', 'ee ', 2, 'float ', nd_bez, nd_bez + write(60) a_ee_print(:,:) + endif + deallocate (a_ee_print, a_ee_loc) +!========================================================================================= + + allocate(s_ww_print(n_w,n_w),stat=ier) + IF (IER /= 0) THEN + WRITE (*,*) "output.f90, can not allocate s_ww_print MYPROC_NUM=",MYPNUM + STOP + END IF + s_ww_print=0. + + call SCA_GRID(n_w,n_w) + DO i_loc = 1,MP_A + IC= INDXL2G( i_loc, NB, MYROW, 0, NPROW) + DO j_loc = 1,NQ_A + JC= INDXL2G( j_loc, NB, MYCOL, 0, NPCOL) + s_ww_print(IC,JC) = s_ww_loc(i_loc,j_loc) + END DO + END DO + + if(rank==0) then + call MPI_REDUCE(MPI_IN_PLACE, s_ww_print, n_w*n_w, & + MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, IER) + else + call MPI_REDUCE(s_ww_print, s_ww_print, n_w*n_w, & + MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, IER) + endif + + if(rank==0) then + write(60) '#@array ', 's_ww ', 2, 'float ', n_w, n_w + write(60) s_ww_print(:,:) + endif + deallocate (s_ww_print, s_ww_loc) +!========================================================================================= + allocate(s_ww_inv_print(n_w,n_w),stat=ier) + IF (IER /= 0) THEN + WRITE (*,*) "output.f90, can not allocate s_ww_inv_print MYPROC_NUM=",MYPNUM + STOP + END IF + s_ww_inv_print=0. + call SCA_GRID(n_w,n_w) + DO i_loc = 1,MP_A + IC= INDXL2G( i_loc, NB, MYROW, 0, NPROW) + DO j_loc = 1,NQ_A + JC= INDXL2G( j_loc, NB, MYCOL, 0, NPCOL) + s_ww_inv_print(IC,JC) = s_ww_inv_loc(i_loc,j_loc) + END DO + END DO + + if(rank==0) then + call MPI_REDUCE(MPI_IN_PLACE, s_ww_inv_print, n_w*n_w, & + MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, IER) + else + call MPI_REDUCE(s_ww_inv_print, s_ww_inv_print, n_w*n_w, & + MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, IER) + endif + + if(rank==0) then + write(60) '#@array ', 's_ww_inv ', 2, 'float ', n_w, n_w + write(60) s_ww_inv_print(:,:) + endif + deallocate (s_ww_inv_print, s_ww_inv_loc) + + +!========================================================================================= + if(rank==0) then + 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(:,:) + + deallocate (xyzpot_w) + endif + deallocate (jpot_w) +!========================================================================================= + + +end if !if (format_type == 'formatted' ) then + +close(60) + +end subroutine print_starwall_response2 diff --git a/src_3d/surface_wall.f90 b/src_3d/surface_wall.f90 index ff37aebe8fe0445a84322bec8f549bc4f6f41986..768e834cee24411eb0b27045d740ed0cc03347a2 100644 --- a/src_3d/surface_wall.f90 +++ b/src_3d/surface_wall.f90 @@ -80,6 +80,13 @@ if(rank==0) write(outp,*) ' n_w m_w rc_w rs_w zc_w end do if(rank==0) write(outp,*) + ! --- Check input parameters + if ( ( nwu*nwv > 19600 ) ) then + write(outp,*) 'ERROR: npot_w (nwu*nwv) must have a value smaller than 19600 for current output version' + stop 1 + end if + + ! Remarks: ! large Phi: independent variables ! small phi: current potential at each triangle vertex