Homec4science
Diffusion FENNECS: Finite Element Non-Neutral Electron Cloud Simulator (stable) Browse

fennecs/src/magnet_mod.f90
2887dfe1f6b9stable

Last Change
Raw File

magnet_mod.f90
View Options

module magnet
use constants
IMPLICIT NONE
type elongatedcoil
real(kind=db):: axialdims(2)
real(kind=db):: radialdims(2)
real(kind=db):: center(2)
integer:: nbsubcoils(2)
real(kind=db):: current
real(kind=db):: Nturns
end type elongatedcoil
type(elongatedcoil), allocatable:: the_coils(:)
contains
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Initialize the magnetic field from the standard input
!> @param[in] fileid file identifier of the input corresponding to the run parameters
!> @param[in] cstep current simulation step number
!---------------------------------------------------------------------------
subroutine magnet_init(fileid,cstep)
use basic, ONLY:B0,rnorm,bnorm,rgrid,zgrid,width,Rcurv,mpirank,Athet,Br,Bz,bscaling,magnetfile
use mpi
IMPLICIT NONE
integer:: fileid, cstep, istat, ierr
integer :: magfiletype = 0
character(256):: line=''
NAMELIST /magnetparams/ magnetfile, magfiletype, bscaling
! read the input parameters from file
Rewind(fileid)
READ(fileid, magnetparams, iostat=istat)
! save the parameters on output
IF(mpirank .eq. 0) THEN
WRITE(*, magnetparams)
END IF
if (istat.gt.0) then
backspace(fileid)
read(fileid,fmt='(A)') line
write(*,'(A)') &
'Invalid line in magnetparams: '//trim(line)
call MPI_Abort(MPI_COMM_WORLD, -1, ierr)
stop
end if
if(magfiletype .eq. 1)then
call read_magfile_txt(magnetfile,the_coils)
call compute_B_on_grid(Athet,Br,Bz,rgrid*rnorm,zgrid*rnorm)
else if(magfiletype .eq. 0 .and. len_trim(magnetfile) .lt. 1)then
call magmirror(Athet,Br,Bz,rgrid,zgrid,rnorm,B0,Rcurv,width)
else
CALL load_mag_from_h5(magnetfile,Athet,Br,Bz,rgrid,zgrid, B0, rnorm,bscaling)
end if
! We normalise the magnetic fields
Br=Br/bnorm
Bz=Bz/bnorm
end subroutine
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Loads the magnetic field defined in the .h5 file at location magfile
!> @param[in] magfile filname of .h5 file containing the definitions of A and B
!---------------------------------------------------------------------------
SUBROUTINE load_mag_from_h5(magfile, Athet,Br,Bz,rgrid,zgrid, B0, rnorm,bscaling)
USE constants, ONLY: Pi
USE futils
USE bsplines
CHARACTER(LEN=*), INTENT(IN):: magfile
REAL(kind=db), ALLOCATABLE :: magr(:), magz(:)
REAL(kind=db), ALLOCATABLE :: tempBr(:, :), tempBz(:, :), tempAthet(:, :)
real(kind=db), allocatable:: c(:,:)
real(kind=db)::B0,rnorm
Integer:: bscaling, nr, nz,i
real(kind=db):: Athet(:), Br(:), Bz(:)
Real(kind=db):: zgrid(0:), rgrid(0:)
Real(kind=db), ALLOCATABLE:: vec1(:), vec2(:)
type(spline2d):: Maginterpolation
REAL(kind=db) :: maxB
INTEGER :: magfid, dims(2)
LOGICAL:: B_is_saved
INTEGER :: magn(2), magrank
CALL openf(trim(magfile), magfid, 'r', real_prec='d')
CALL getdims(magfid, '/mag/Athet', magrank, magn)
nr=size(rgrid,1)-1
nz=size(zgrid,1)-1
! Auxiliary vectors
ALLOCATE(vec1((nz+1)*(nr+1)),vec2((nr+1)*(nz+1)))
DO i=0,nr
vec1(i*(nz+1)+1:(i+1)*(nz+1))=zgrid!(0:nz)
vec2(i*(nz+1)+1:(i+1)*(nz+1))=rgrid(i)
END DO
ALLOCATE (magr(magn(2)), magz(magn(1)))
ALLOCATE (tempAthet(magn(1), magn(2)), tempBr(magn(1), magn(2)), tempBz(magn(1), magn(2)))
! Read r and z coordinates for the definition of A_\thet, and B
CALL getarr(magfid, '/mag/r', magr)
CALL getarr(magfid, '/mag/z', magz)
CALL getarr(magfid, '/mag/Athet', tempAthet)
IF (isdataset(magfid, '/mag/Br') .and. isdataset(magfid, '/mag/Bz')) THEN
CALL getarr(magfid, '/mag/Br', tempBr)
CALL getarr(magfid, '/mag/Bz', tempBz)
IF(bscaling .gt. 0) then
maxB=sqrt(maxval(tempBr**2+tempBz**2))
tempBr=tempBr/maxB*B0
tempBz=tempBz/maxB*B0
end if
B_is_saved = .true.
ELSE
B_is_saved = .false.
END IF
magz=magz/rnorm
magr=magr/rnorm
CALL set_splcoef((/3,3/),magz,magr,Maginterpolation)
call get_dim(Maginterpolation,dims)
! Interpolation of the magnetic potential vector
allocate(c(dims(1),dims(2)))
call get_splcoef(Maginterpolation,tempAthet, c)
CALL gridval(Maginterpolation,vec1,vec2, Athet ,(/0,0/),c)
if(B_is_saved .eqv. .true.)then
! Interpolation of the Axial magnetic field
call get_splcoef(Maginterpolation,tempBz, c)
CALL gridval(Maginterpolation,vec1,vec2, Bz ,(/0,0/),c)
! Interpolation of the radial magnetic field
call get_splcoef(Maginterpolation,tempBr, c)
CALL gridval(Maginterpolation,vec1,vec2, Br ,(/0,0/),c)
else
CALL gridval(Maginterpolation,vec1,vec2, Br,(/1,0/))
Br=-Br
CALL gridval(Maginterpolation,vec1,vec2, Bz,(/0,1/))
Bz=Bz+Athet/vec2
end if
if( bscaling .lt. 0 ) then
maxB = maxval(sqrt(Bz**2 + Br**2))
Bz = Bz/maxB*B0
Br = Br/maxB*B0
end if
CALL closef(magfid)
deallocate(c)
call destroy_SP(Maginterpolation)
END SUBROUTINE load_mag_from_h5
!---------------------------------------------------------------------------
!> @author
!> Patryk kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Computes the magnetic field on the grid according to a magnetic mirror
!> @param[out] Athet magnetic potential vector
!> @param[out] Br radial magnetic field on the grid
!> @param[out] Bz Axial magnetic field on the grid
!> @param[in] rgrid radial grid points
!> @param[in] zgrid axial grid points
!> @param[in] rnorm normalisation distance constant
!> @param[in] B0 maximum mag field amplitude on axis
!> @param[in] Rcurv magnetic mirror R=max(B)/min(B) parameter
!> @param[in] width distance between the two mirror coils [m]
!>
!---------------------------------------------------------------------------
subroutine magmirror(Athet,Br,Bz,rgrid,zgrid,rnorm,B0,Rcurv,width)
USE constants, ONLY: Pi
real(kind=db):: Athet(:), Br(:), Bz(:)
REAL(kind=db) :: rg, zg, halfLz, rnorm, Rcurv, MirrorRatio, width, B0
Real(kind=db):: zgrid(:), rgrid(:)
INTEGER :: i, rindex, nr, nz
nr=size(rgrid,1)-1
nz=size(zgrid,1)-1
halfLz = (zgrid(nz) + zgrid(0))/2
MirrorRatio = (Rcurv - 1)/(Rcurv + 1)
DO i = 1, (nr + 1)*(nz + 1)
rindex = (i - 1)/(nz + 1)
rg = rgrid(rindex)
zg = zgrid(i - rindex*(nz + 1) - 1) - halfLz
Br(i) = -B0*MirrorRatio*SIN(2*pi*zg/width*rnorm)*bessi1(2*pi*rg/width*rnorm)
Bz(i) = B0*(1 - MirrorRatio*COS(2*pi*zg/width*rnorm)*bessi0(2*pi*rg/width*rnorm))
Athet(i) = 0.5*B0*(rg*rnorm - width/pi*MirrorRatio*bessi1(2*pi*rg/width*rnorm)*COS(2*pi*zg/width*rnorm))
END DO
end subroutine
!________________________________________________________________________________
!Modified Bessel functions of the first kind of the zero order
FUNCTION bessi0(x)
REAL(kind=db) :: bessi0, x
REAL(kind=db) :: ax
REAL(kind=db) p1, p2, p3, p4, p5, p6, p7, q1, q2, q3, q4, q5, q6, q7, q8, q9, y
SAVE p1, p2, p3, p4, p5, p6, p7, q1, q2, q3, q4, q5, q6, q7, q8, q9
DATA p1, p2, p3, p4, p5, p6, p7/1.0d0, 3.5156229d0, 3.0899424d0, 1.2067492d0, 0.2659732d0, 0.360768d-1, 0.45813d-2/
DATA q1, q2, q3, q4, q5, q6, q7, q8, q9/0.39894228d0, 0.1328592d-1, 0.225319d-2, -0.157565d-2, 0.916281d-2, &
& -0.2057706d-1, 0.2635537d-1, -0.1647633d-1, 0.392377d-2/
if (abs(x) .lt. 3.75) then
y = (x/3.75)**2
bessi0 = p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*(p6 + y*p7)))))
else
ax = abs(x)
y = 3.75/ax
bessi0 = (exp(ax)/sqrt(ax))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*q9))))))))
end if
return
END FUNCTION bessi0
!________________________________________________________________________________
!Modified Bessel functions of the first kind of the first order
FUNCTION bessi1(x)
REAL(kind=db) :: bessi1, x
REAL(kind=db) :: ax
REAL(kind=db) p1, p2, p3, p4, p5, p6, p7, q1, q2, q3, q4, q5, q6, q7, q8, q9, y
SAVE p1, p2, p3, p4, p5, p6, p7, q1, q2, q3, q4, q5, q6, q7, q8, q9
DATA p1, p2, p3, p4, p5, p6, p7/0.5d0, 0.87890594d0, 0.51498869d0, 0.15084934d0, 0.2658733d-1, 0.301532d-2, 0.32411d-3/
DATA q1, q2, q3, q4, q5, q6, q7, q8, q9/0.39894228d0, -0.3988024d-1, -0.362018d-2, 0.163801d-2, -0.1031555d-1, &
& 0.2282967d-1, -0.2895312d-1, 0.1787654d-1, -0.420059d-2/
if (abs(x) .lt. 3.75D0) then
y = (x/3.75D0)**2
bessi1 = x*(p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*(p6 + y*p7))))))
else
ax = abs(x)
y = 3.75D0/ax
bessi1 = (exp(ax)/sqrt(ax))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*q9))))))))
if (x .lt. 0.) bessi1 = -bessi1
end if
return
END FUNCTION bessi1
! calculation of the magnetic potential vector and magnetic fields
! using elliptic integrals
! see Smythe "static and dynamic electricity" McGraw-Hill, third edition p.290
subroutine compute_B_on_grid(A,Br,Bz,r,z)
use elliptic
use constants
implicit none
real(kind=db):: A(:), Br(:), Bz(:)
real(kind=db):: r(:), z(:)
integer:: i,j,k, l, n,linindex, nr,nz, ncoils
Real(kind=db):: kappa2, rho, distz, dr, dz
Real(kind=db):: r_wire, z_wire
Real(kind=db):: r_cen, z_cen, u, v, w, t
Real(kind=db):: EE, EK
Real(kind=db):: Isubcoil, coilsurf, Itot, currfac
nr=size(r,1)
nz=size(z,1)
ncoils=size(the_coils,1)
!$OMP PARALLEL DO PRIVATE(j,i,linindex,k,dr,dz,z_cen,r_cen,coilsurf,Itot,Isubcoil,l,z_wire,n,r_wire,kappa2,EE,EK, u,v,w,t,distz, currfac)
Do j=1,nr
Do i=1,nz
linindex=i+(j-1)*size(z,1)
A(linindex)=0.0_db
Bz(linindex)=0.0_db
Br(linindex)=0.0_db
Do k=1,ncoils
dr = (the_coils(k)%radialdims(2)-the_coils(k)%radialdims(1))
dz = (the_coils(k)%axialdims(2)-the_coils(k)%axialdims(1))
r_cen = (the_coils(k)%radialdims(2)+the_coils(k)%radialdims(1))/2
z_cen = (the_coils(k)%axialdims(2)+the_coils(k)%axialdims(1))/2
coilsurf = dr * dz;
Itot = the_coils(k)%current*the_coils(k)%Nturns
Isubcoil = Itot/(the_coils(k)%nbsubcoils(1)*the_coils(k)%nbsubcoils(2))
currfac = mu_0*Isubcoil/pi
Do l=1,the_coils(k)%nbsubcoils(1)
z_wire = (z_cen - dz/2) + (dz/the_coils(k)%nbsubcoils(1))*(l-0.5)
distz = (z(i)-z_wire)
Do n=1,the_coils(k)%nbsubcoils(2)
r_wire=(r_cen - dr/2) + (dr/the_coils(k)%nbsubcoils(2))*(n-0.5);
u=(r_wire+r(j))**2+distz**2
v=(r_wire-r(j))**2+distz**2
w=r_wire**2+r(j)**2+distz**2
t=r_wire**2-r(j)**2-distz**2
kappa2=4*r_wire*r(j)/u
EK=elliptic_fk(sqrt(kappa2)) ! complete Elliptic integral of the first kind
EE=elliptic_ek(sqrt(kappa2)) ! complete Elliptic integral of the second kind
if(r(j).gt.0.0_db)then
A(linindex)=A(linindex)+currfac*sqrt(r_wire/(r(j)*kappa2))*&
& ((1-kappa2/2)*EK-EE)
Br(linindex)=Br(linindex)+currfac/2* &
& distz/(r(j)*sqrt(u))*(-EK+w/v*EE)
end if
Bz(linindex)=Bz(linindex)+currfac/2* &
& 1/sqrt(u)*(EK+t/v*EE)
end do
end do
end do
end do
!WRITE(*,*)"id done",j,'over',nr
end do
!$OMP END PARALLEL DO
end subroutine
subroutine read_magfile_txt(magnetfile,the_coils)
type(elongatedcoil), allocatable:: the_coils(:)
character(128):: magnetfile
INTEGER:: nb_coils
INTEGER :: lu_magfile=999
INTEGER:: i, openerr, reason
CHARACTER(len=256) :: header
real(kind=db):: z1,z2, r1, r2, current
integer:: Nr, Nz
real(kind=db):: Nturns
if(allocated(the_coils)) deallocate(the_coils)
nb_coils=0
OPEN(UNIT=lu_magfile,FILE=trim(magnetfile),ACTION='READ',IOSTAT=openerr)
header=' '
IF(openerr .ne. 0) THEN
CLOSE(unit=lu_magfile)
RETURN
END IF
! The magnet table is defined as a 8 column (zmin zmax rmin rmax nturns Nz Nr current)
DO WHILE(.true.)
READ(lu_magfile,'(a)',IOSTAT=reason) header
header=adjustl(header)
if(reason .lt. 0 ) exit ! We reached end of file
if( header(1:1) .ne. '!'.and. len_trim(header) .gt. 1) then
READ(header,*,IOSTAT=reason) z1,z2, r1, r2, nturns, Nz, Nr, current
if(reason .lt. 0 ) then
WRITE(*,*) "Error in magnet definition of file: ",magnetfile,"\n with input ", header
exit
end if
nb_coils=nb_coils+1
end if
END DO
allocate(the_coils(nb_coils))
REWIND(lu_magfile)
! The magnet table is defined as a 8 column (zmin zmax rmin rmax nturns Nz Nr current)
i=1
write(*,*) "number of found coils", nb_coils
DO WHILE(i .le. nb_coils)
READ(lu_magfile,'(a)',IOSTAT=reason) header
header=adjustl(header)
!Write(*,*) i,' ',trim(header)
if(reason .lt. 0 ) exit ! We reached end of file
if( header(1:1) .ne. '!' .and. len_trim(header) .gt. 1) then
READ(header,*,IOSTAT=reason) z1,z2, r1, r2, nturns, Nz, Nr, current
the_coils(i)%axialdims=(/z1,z2/)
the_coils(i)%radialdims=(/r1,r2/)
the_coils(i)%Nturns=nturns
the_coils(i)%nbsubcoils=(/Nz,Nr/)
the_coils(i)%current=current
!Write(*,*) i,' ',the_coils(i)%axialdims, the_coils(i)%radialdims, the_coils(i)%Nturns,the_coils(i)%nbsubcoils,the_coils(i)%current
i=i+1
end if
END DO
CLOSE(unit=lu_magfile)
end subroutine
end module magnet
c4science · Help