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