Commit 22502cc5 authored by Serhiy Mochalskyy's avatar Serhiy Mochalskyy
Browse files

magnetic coil version with clean code after forcheck

parent cc31c56e
......@@ -8,7 +8,7 @@ ifeq ($(DEBUG),yes)
else
FFLAGS += -O3
# FFLAGS += -O0 -warn all,nounused -check all,noarg_temp_created -debug all -debug-parameters -fstack-security-check -ftrapuv -traceback
FFLAGS += -O0 -warn all,nounused -check all,noarg_temp_created -debug all -debug-parameters -fstack-security-check -ftrapuv -traceback
endif
FPPFLAGS =
......
subroutine a_ee_computing
use solv
use tri_w
use tri_p
use mpi_v
use sca
use coil2d
!-----------------------------------------------------------------------
implicit none
include "mpif.h"
if(rank==0) write(*,*) 'a_ee_computing starts, nd_bez=', nd_bez
......
subroutine a_ew_computing
use solv
use tri_w
use tri_p
use mpi_v
use sca
use coil2d
!-----------------------------------------------------------------------
implicit none
include "mpif.h"
if(rank==0) write(*,*) "a_ew_computing start"
......
subroutine a_ey_computing
use solv
use tri_w
use mpi_v
use sca
use resistive
!-----------------------------------------------------------------------
implicit none
include "mpif.h"
if(rank==0) write(*,*) "a_ey_computing start"
......
......@@ -8,7 +8,7 @@ use coil2d
use mpi_v
implicit none
include "mpif.h"
if(rank==0) write(*,*) 'a_pwe transpose begins'
......
subroutine a_pwe_s_computing
use solv
use tri_w
use tri_p
use mpi_v
use sca
......@@ -9,7 +8,7 @@ use sca
implicit none
include "mpif.h"
integer :: i,j,k
integer :: i,j
integer BLACS_PNUM
EXTERNAL BLACS_PNUM
......
subroutine a_we_computing
use solv
use tri_w
use tri_p
use mpi_v
use sca
use coil2d
!-----------------------------------------------------------------------
implicit none
include "mpif.h"
integer :: i,j,k
if(rank==0) write(*,*) 'a_we_computing starts'
......
subroutine a_ye_computing
use solv
use tri_w
use mpi_v
use sca
use resistive
!-----------------------------------------------------------------------
implicit none
include "mpif.h"
integer :: i_loc,i
logical :: inside_i
......@@ -50,8 +49,8 @@ CALL PDGEMM('T','N', n_w ,nd_bez, n_w , 1.0D0 , s_ww_loc ,1,1,DESCA, a_we_loc,1
do i=1,n_w
call ScaLAPACK_mapping_i2(i,i_loc,inside_i)
if (inside_i /= .true.) then
call ScaLAPACK_mapping_i(i,i_loc,inside_i)
if (inside_i) then
a_ye_loc(i_loc,:) = a_ye_loc(i_loc,:)/gamma(i)
endif
enddo
......
!----------------------------------------------------------------------
subroutine bezier2real_space(bez,f,N_bnd,n_points,ndim)
! ----------------------------------------------------------------------
! purpose: 12/08/09
!
! ----------------------------------------------------------------------
implicit none
integer :: N_bnd,n_points,ndim
real,dimension(2*N_bnd, ndim) :: bez
real,dimension(N_bnd*n_points,ndim) :: f
real,dimension(n_points ) :: h11,h12,h21,h22
real,dimension(2 , 2) :: H,H_s,H_ss
real :: s
integer :: i,idim,im,im2,ip,ip2,j,k
f = 0.
do j=1,n_points
s = float(j-1)/float(n_points)
call basisfunctions1(s,H,H_s,H_ss)
h11(j) = H(1,1)
h12(j) = H(1,2)
h21(j) = H(2,1)
h22(j) = -H(2,2)
enddo
do i = 1, N_bnd
im = 2*i -1
ip = 2*i
im2 = im + 2
ip2 = ip + 2
if(im2.gt.2*N_bnd) im2 = 1
if(ip2.gt.2*N_bnd) ip2 = 2
do k=1,n_points
j = k+n_points*(i-1)
f(j,:) = f(j,:) + h11(k)*bez(im ,:) &
+ h21(k)*bez(im2,:) &
+ h12(k)*bez(ip ,:) &
+ h22(k)*bez(ip2,:)
enddo
enddo
end subroutine bezier2real_space
!----------------------------------------------------------------------
subroutine bezier_s_2real_space(bez,f,N_bnd,n_points,ndim)
! ----------------------------------------------------------------------
! purpose: 12/08/09
!
! ----------------------------------------------------------------------
implicit none
integer :: N_bnd,n_points,ndim
real,dimension(2*N_bnd, ndim) :: bez
real,dimension(N_bnd*n_points,ndim) :: f
real,dimension(n_points ) :: h11,h12,h21,h22
real,dimension(2 , 2) :: H,H_s,H_ss
real :: s
integer :: i,idim,im,im2,ip,ip2,j,k
f = 0.
do j=1,n_points
s = float(j-1)/float(n_points)
call basisfunctions1(s,H,H_s,H_ss)
h11(j) = H_s(1,1)
h12(j) = H_s(1,2)
h21(j) = H_s(2,1)
h22(j) = -H_s(2,2)
enddo
do i = 1, N_bnd
im = 2*i -1
ip = 2*i
im2 = im + 2
ip2 = ip + 2
if(im2.gt.2*N_bnd) im2 = 1
if(ip2.gt.2*N_bnd) ip2 = 2
do k=1,n_points
j = k+n_points*(i-1)
f(j,:) = f(j,:) + h11(k)*bez(im ,:) &
+ h21(k)*bez(im2,:) &
+ h12(k)*bez(ip ,:) &
+ h22(k)*bez(ip2,:)
enddo
enddo
end subroutine bezier_s_2real_space
......@@ -15,11 +15,10 @@
real,dimension (n):: 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 &
,ata1,ata2,ata3,at &
,s21,s22,s23,dp1,dm1,dp2,dm2,dp3,dm3 &
,ap1,am1,ap2,am2,ap3,am3 &
,h,ln1,ln2,ln3,ar1,ar2,ar3 &
,h,ar1,ar2,ar3 &
,x21,y21,z21,x32,y32,z32,x13,y13,z13,vx,vy,vz&
,tx1,ty1,tz1,tx2,ty2,tz2,tx3,ty3,tz3 &
,nx,ny,nz,pi41,area,d21,d32,d13,jx,jy,jz &
......
......@@ -4,9 +4,7 @@ subroutine bfield_par_par(b_par, xp, yp, zp, xu, yu, zu, n, n_b, n_e, x, y, z, n
! 23.09.09
! purpose:
!
!
! ----------------------------------------------------------------------
use mpi_v
implicit none
integer :: n,n_b,n_e,ntri
......@@ -19,11 +17,10 @@ integer,dimension (ntri,3) :: j_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 &
,ata1,ata2,ata3,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 &
,h,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 &
......
program check_bezier
real*8, allocatable :: f(:,:), bez(:,:), index(:), freco(:), Lij(:,:,:), Rij(:,:,:), Zij(:,:,:)
real*8 :: H(2,2)
integer,allocatable :: bnd_node(:,:), bnd_node_index(:,:), bnd_node_ndof(:)
open(921,file='fstar.dat',form='unformatted')
read(921) n_f,ndim
allocate(f(n_f,ndim))
read(921) f
close(921)
open(922,file='bstar.dat',form='unformatted')
read(922) n_dof_bnd,ndim
allocate(bez(n_dof_bnd,ndim))
read(922) bez
close(922)
write(*,*) ' f : ',n_f,ndim
write(*,*) ' bez : ',n_dof_bnd,ndim
!------------------------------------------------------------------------------
open(22,file='boundary.txt')
read(22,*) n_bnd, n_bnd_nodes ! the number of boundary elements, nodes
allocate(Lij(2,2,n_bnd),Rij(2,2,n_bnd), Zij(2,2,n_bnd), bnd_node(n_bnd,2))
!----- Rij : first index counts the two nodes of a boundary element
! second index counts the basis function H
! third index counts the number of the boundary element
do i=1,n_bnd
read(22,*) ibnd, bnd_node(i,1), bnd_node(i,2), &
Rij(1,1,i), Zij(1,1,i), Rij(1,2,i), Zij(1,2,i), Lij(1,1,i), Lij(1,2,i), &
Rij(2,1,i), Zij(2,1,i), Rij(2,2,i), Zij(2,2,i), Lij(2,1,i), Lij(2,2,i)
enddo
allocate(bnd_node_index(n_bnd_nodes,2),bnd_node_ndof(n_bnd_nodes))
do i=1,n_bnd_nodes
read(22,*) j,bnd_node_index(i,:),bnd_node_ndof(i) ! the index in the starwall response (as requested by JOREK)
write(*,'(A,3i5)') 'bnd_node_index : ',i,bnd_node_index(i,:)
enddo
close(22)
n_bnd_elm = n_bnd
n_points = 10
allocate(freco(n_bnd_elm*n_points))
freco = 0.
n = 32
do j=1,n_points
s = float(j-1)/10.
call basisfunctions1(s,H,H_s,H_ss)
do k=1, n_bnd_elm
im = bnd_node_index(bnd_node(k,1),1) ! node 1, first dof
ip = bnd_node_index(bnd_node(k,1),2) ! node 1, second dof
im2 = bnd_node_index(bnd_node(k,2),1) ! node 2, first dof
ip2 = bnd_node_index(bnd_node(k,2),2) ! node 2, second dof
freco((k-1)*n_points+j) = freco((k-1)*n_points+j) &
+ bez(im,n) * H(1,1) * Lij(1,1,k) &
+ bez(ip,n) * H(1,2) * Lij(1,2,k) &
+ bez(im2,n) * H(2,1) * Lij(2,1,k) &
+ bez(ip2,n) * H(2,2) * Lij(2,2,k)
enddo
enddo
do i=1,n_points*n_bnd_elm
write(*,*) i,f(i,n),freco(i)
enddo
allocate(index(ndim))
np1 = 20 !ndim/16
np2 = n_f
do i=1,np2
index(i) = i
enddo
call begplt('check.ps')
call lplot6(1,1,index(np1:np2),f(np1:np2,n),np2-np1+1,' ')
call lplot6(1,1,index(np1:np2),freco(np1:np2),-(np2-np1+1),' ')
call finplt
end
subroutine basisfunctions1(s,H,H_s,H_ss)
!--------------------------------------------------------------
! subroutine which defines the basis functions in one dimension
!
! derived from a mixed Bezier/Cubic finite element representation
! index 1 : counts the vertex
! index 2 : counts the variables (p,u,v,w)
!
! the functions are defined on the interval [0,1]
!----------------------------------------------------------------
implicit none
real*8 :: s, H(2,2), H_s(2,2), H_ss(2,2)
!---------------------------------------------------------- vertex (1)
H(1,1) =(-1.d0 + s)**2*(1.d0 + 2.d0*s)
H_s(1,1) =6.d0*(-1.d0 + s)*s
H_ss(1,1)=6.d0*(-1.d0 + 2.d0*s)
H(1,2) =3.d0*(-1.d0 + s)**2*s
H_s(1,2) =3.d0*(-1.d0 + s)*(-1.d0 + 3.d0*s)
H_ss(1,2)=6.d0*(-2.d0 + 3.d0*s)
!---------------------------------------------------------- vertex (2)
H(2,1) =-s**2*(-3.d0 + 2.d0*s)
H_s(2,1) =-6.d0*(-1.d0 + s)*s
H_ss(2,1)=-6.d0*(-1.d0 + 2.d0*s)
H(2,2) =-3.d0*(-1.d0 + s)*s**2
H_s(2,2) =-3.d0*s*(-2.d0 + 3.d0*s)
H_ss(2,2)=-6.d0*(-1.d0 + 3.d0*s)
end subroutine basisfunctions1
subroutine cholesky_solver
use solv
use tri_w
use tri_p
use mpi_v
use sca
!-----------------------------------------------------------------------
implicit none
include "mpif.h"
if(rank==0) write(*,*) 'cholesky_solver begins'
!=================================================================================
......
......@@ -4,7 +4,7 @@ use sca
use resistive
!-----------------------------------------------------------------------
implicit none
include "mpif.h"
integer, allocatable :: ipiv(:)
real, dimension(:),allocatable :: WORK
......
subroutine control_boundary
!------------------------------------------------------------------------------
use icontr
use gauss
use contr_su
use mpi_v
......@@ -10,7 +9,7 @@ implicit none
include "mpif.h"
real,dimension(:,:,:),allocatable :: Rij, Zij
real :: pi2,alv,fnv,fnp,sq,sq1,co,si,s
integer :: i, j, k, l, iv1, iv2,index,ibnd,ku,kv,nuv
integer :: i, j, k, l, index,ibnd,ku,kv,nuv
integer :: j_exists
logical :: exists
......
subroutine d_ee_computing
use solv
use tri_w
use mpi_v
use sca
use resistive
!-----------------------------------------------------------------------
implicit none
include "mpif.h"
if(rank==0) write(*,*) "d_ee_computing start"
......
F_MAIN_SRC = mod_icontr.f90\
mod_gauss.f90\
mod_contr_su.f90\
mod_tri_w.f90\
mod_tri_p.f90\
......@@ -15,7 +14,6 @@ F_MAIN_SRC = mod_icontr.f90\
input.f90\
solver.f90\
coil_2d.f90\
linequ.f90\
simil_trafo.f90\
resistive_wall_response.f90\
bfield_c.f90\
......@@ -35,7 +33,6 @@ F_MAIN_SRC = mod_icontr.f90\
tri_contr_surf.f90\
real_space2bezier.f90\
real_space2bezier_par.f90\
real_space2bezier_s.f90\
read_coil_data.f90\
matrix_ep.f90\
matrix_ew.f90\
......@@ -59,8 +56,7 @@ F_MAIN_SRC = mod_icontr.f90\
computing_s_ww_inverse.f90\
get_index_dima.f90\
output.f90\
print_starwall_response.f90\
ideal_wall_response.f90
print_starwall_response.f90
F_MAIN_OBJ = $(addprefix $(OBJ_DIR)/,$(F_MAIN_SRC:.f90=.o))
MAIN_OBJ = $(F_MAIN_OBJ)
......@@ -80,7 +76,7 @@ $(OBJ_DIR)/input.o: \
$(OBJ_DIR)/mod_mpi_v.o $(OBJ_DIR)/mod_resistive.o
$(OBJ_DIR)/control_boundary.o: \
$(OBJ_DIR)/mod_icontr.o $(OBJ_DIR)/mod_contr_su.o\
$(OBJ_DIR)/mod_gauss.o $(OBJ_DIR)/mod_mpi_v.o
$(OBJ_DIR)/mod_mpi_v.o
$(OBJ_DIR)/matrix_pp.o: \
$(OBJ_DIR)/mod_icontr.o $(OBJ_DIR)/mod_solv.o \
$(OBJ_DIR)/mod_contr_su.o $(OBJ_DIR)/mod_tri_p.o\
......
!------------------------------------------------------------------------
subroutine generate_index_fu(tx,ty,tz,ntri,j_pot,nupot,nhole)
!------------------------------------------------------------------------
implicit none
real, dimension(ntri,3) :: tx,ty,tz
integer,dimension(ntri,3) :: j_pot,nbou
integer :: nhole,nupot,j
real :: delta,dis,disp,tx1,ty1,tz1,tx2,ty2,tz2
integer :: ntri,k,i,kp,ip,k1,kp1
!------------------------------------------------------------------------
!find the boundaries of a multiply connected wall composed of triangles
! nbou(i,k) = 0 : edge (i,k)-(i,k+1) interior edge
! nbou(i,k) = -1 : edge (i,k)-(i,k+1) boundary edge
delta = 1.e-10
nbou = -1
do k=1,3
do i=1,ntri
if(nbou(i,k).lt.0) then
do kp=1,3
do ip=i+1,ntri
dis = abs(tx(ip,kp )-tx(i,k )) &
+ abs(ty(ip,kp )-ty(i,k )) &
+ abs(tz(ip,kp )-tz(i,k ))
if(dis.le.delta) then
kp1 = mod(kp+1,3) + 1
k1 = mod(k ,3) + 1
disp = abs(tx(ip,kp1)-tx(i,k1)) &
+ abs(ty(ip,kp1)-ty(i,k1)) &
+ abs(tz(ip,kp1)-tz(i,k1))
if(disp.le.delta) then
nbou(i ,k ) = 0
nbou(ip,kp1) = 0
endif
endif
enddo
enddo
endif
enddo
enddo
! find the boundary of the holes
nhole = 0
do i=1,ntri
do k=1,3
if(nbou(i,k).eq.-1) then
nhole = nhole+1
nbou(i,k) = nhole
tx1 = tx(i,k)
ty1 = ty(i,k)
tz1 = tz(i,k)
k1 = mod(k,3)+1
tx2 = tx(i,k1)
ty2 = ty(i,k1)
tz2 = tz(i,k1)
98 continue
do ip=i+1,ntri
do kp=1,3
if(nbou(ip,kp).eq.-1) then
kp1 = mod(kp,3) +1
dis = abs(tx2-tx(ip,kp)) &
+abs(ty2-ty(ip,kp)) &
+abs(tz2-tz(ip,kp))
if(dis.le.delta) then
nbou(ip,kp) = nhole
tx2 = tx(ip,kp1)
ty2 = ty(ip,kp1)
tz2 = tz(ip,kp1)
disp = abs(tx1-tx(ip,kp1)) &
+ abs(ty1-ty(ip,kp1)) &
+ abs(tz1-tz(ip,kp1))
if(disp.le.delta) goto 99
go to 98
endif
endif
enddo
enddo
write(6,*) 'boundary nr.:',nhole+1,' is not closed: stop'
stop
99 continue
write(6,*) ' boundary nr.:',nhole,'is closed'
endif
enddo
enddo
write(6,*) 'number of boundaries: ',nhole
! compute the index function
j_pot = 0
if(nhole.gt.0) then
do k=1,3
do i=1,ntri
if(nbou(i,k).gt.0.and.j_pot(i,k).eq.0) then
! if(nbou(i,k).gt.0) then
j_pot(i,k) = nbou(i,k)
do kp=1,3
do ip=1,ntri
dis = abs(tx(i,k)-tx(ip,kp)) &
+abs(ty(i,k)-ty(ip,kp)) &
+abs(tz(i,k)-tz(ip,kp))
if(dis.le.delta) &
j_pot(ip,kp) = j_pot(i,k)
enddo
enddo
endif
enddo
enddo
endif
j = nhole
do i=1,ntri
do k=1,3
if(j_pot(i,k).eq.0) then
j = j+1
j_pot(i,k) = j
do ip=i+1,ntri
do kp=1,3
dis = abs(tx(i</