diff --git a/src/basic_mod.f90 b/src/basic_mod.f90 index 7afade0..817e1a1 100644 --- a/src/basic_mod.f90 +++ b/src/basic_mod.f90 @@ -1,311 +1,312 @@ MODULE basic ! USE hashtable USE constants USE bsplines USE mumps_bsplines USE futils USE mpihelper IMPLICIT NONE ! ! Basic module for time dependent problems ! CHARACTER(len=80) :: label1, label2, label3, label4 ! ! BASIC Namelist ! LOGICAL :: nlres = .FALSE. !< Restart flag LOGICAL :: nlsave = .TRUE. !< Checkpoint (save) flag LOGICAL :: newres=.FALSE. !< New result HDF5 file LOGICAL :: nlxg=.FALSE. !< Show graphical interface Xgrafix INTEGER :: nrun=1 !< Number of time steps to run DOUBLE PRECISION :: job_time=3600.0 !< Time allocated to this job in seconds DOUBLE PRECISION :: tmax=100000.0 !< Maximum simulation time DOUBLE PRECISION :: extra_time=60.0 !< Extra time allocated DOUBLE PRECISION :: dt=1 !< Time step DOUBLE PRECISION :: time=0 !< Current simulation time (Init from restart file) ! ! Other basic global vars and arrays ! INTEGER :: jobnum !< Job number INTEGER :: step !< Calculation step of this run INTEGER :: cstep=0 !< Current step number (Init from restart file) LOGICAL :: nlend !< Signal end of run INTEGER :: ierr !< Integer used for MPI INTEGER :: it0d=1 !< Number of iterations between 0d values writes to hdf5 INTEGER :: it2d=1 !< Number of iterations between 2d values writes to hdf5 INTEGER :: itparts=1 !< Number of iterations between particles values writes to hdf5 INTEGER :: ittext=1 !< Number of iterations between text outputs in the console INTEGER :: itrestart=10000 !< Number of iterations between save of restart.h5 file INTEGER :: itgraph !< Number of iterations between graphical interface updates INTEGER :: mpirank !< MPIrank of the current processus INTEGER :: mpisize !< Size of the MPI_COMM_WORLD communicator INTEGER :: rightproc !< Rank of next processor in the z decomposition INTEGER :: leftproc !< Rank of previous processor in the z decomposition ! ! List of logical file units INTEGER :: lu_in = 90 !< File duplicated from STDIN INTEGER :: lu_stop = 91 !< stop file, see subroutine TESEND ! ! HDF5 file CHARACTER(len=64) :: resfile = "results.h5" !< Main result file CHARACTER(len=64) :: rstfile = "restart.h5" !< Restart file INTEGER :: fidres !< File ID for resfile INTEGER :: fidrst !< File ID for restart file TYPE(BUFFER_TYPE) :: hbuf0 !< Hashtable for 0d var ! ! Plasma parameters LOGICAL :: nlPhis= .TRUE. !< Calculate self consistent electric field flag LOGICAL :: nlclassical= .FALSE. !< If true, solves the equation of motion according to classical !! dynamics LOGICAL :: nlperiod(2)=(/.false.,.false./)!< Set periodic splines on or off LOGICAL :: partperiodic= .TRUE. !< Sets if the particles boundary conditions are periodic or open INTEGER :: nplasma !< Number of macro-particles INTEGER :: nblock !< Number of slices in Z for stable distribution initialisation DOUBLE PRECISION :: potinn !< Electric potential at the inner metallic wall DOUBLE PRECISION :: potout !< Electric potential at the outer metallic wall DOUBLE PRECISION :: B0 !< Max magnitude of magnetic field DOUBLE PRECISION, allocatable :: Bz(:), Br(:) !< Magnetic field components DOUBLE PRECISION, allocatable :: Athet(:) !< Theta component of the magnetic vector potential TYPE(spline2d) :: splrz !< Spline at r and z DOUBLE PRECISION, allocatable :: Ez(:), Er(:) !< Electric field components DOUBLE PRECISION, allocatable :: pot(:) !< Electro static potential DOUBLE PRECISION :: radii(4) !< Inner and outer radius of cylinder and radii of fine mesh region DOUBLE PRECISION :: plasmadim(4) !< Zmin Zmax Rmin Rmax values for plasma particle loading INTEGER :: distribtype=1 !< Type of distribution function used to load the particles !!1: gaussian, 2: Stable as defined in 4.95 of Davidson DOUBLE PRECISION :: H0=0 !< Initial value of Hamiltonian for distribution 2 DOUBLE PRECISION :: P0=0 !< Initial canonical angular momentum for distribution 2 DOUBLE PRECISION :: temprescale = -1.0 !< Factor used for temperature rescaling in case of a restart (<0 -> no rescaling) DOUBLE PRECISION :: lz(2) !< Lower and upper cylinder limits in z direction DOUBLE PRECISION :: n0 !< Physical plasma density parameter DOUBLE PRECISION, ASYNCHRONOUS, DIMENSION(:,:), ALLOCATABLE, SAVE:: moments !< Moments of the distribution function evaluated every it2d DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, SAVE:: rhs !< right hand side of the poisson equation solver + DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, SAVE:: volume !< Volume covered by each spline for density calculation INTEGER :: nz !< Number of grid intervals in z INTEGER :: nr !< Total number of grid intervals in r INTEGER :: nnr(3) !< Number of grid intervals in r in each subdomain DOUBLE PRECISION :: dz !< Cell size in z DOUBLE PRECISION :: dr(3) !< Cell size in r for each region DOUBLE PRECISION, ALLOCATABLE :: zgrid(:) !< Nodes positions in longitudinal direction DOUBLE PRECISION, ALLOCATABLE :: rgrid(:) !< Nodes positions in radial direction DOUBLE PRECISION :: bnorm,enorm,vnorm,tnorm,rnorm,phinorm !< Normalization constants DOUBLE PRECISION :: qsim !< Charge of superparticles DOUBLE PRECISION :: msim !< Mass of superparticles INTEGER :: femorder(2) !< FEM order INTEGER :: ngauss(2) !< Number of gauss points LOGICAL :: nlppform =.TRUE. !< Argument of set_spline INTEGER, SAVE :: nrank(2) !< Number of splines in both directions DOUBLE PRECISION :: omegac !< Cyclotronic frequency DOUBLE PRECISION :: omegap !< Plasma frequency DOUBLE PRECISION :: V !< Normalised volume DOUBLE PRECISION :: temp !< Initial temperature of plasma DOUBLE PRECISION :: Rcurv !< Magnetic field curvature coefficient DOUBLE PRECISION :: Width !< Distance between two magnetic mirrors INTEGER, DIMENSION(:), ALLOCATABLE :: Zbounds !< Index of bounds for local processus in Z direction TYPE(BASICDATA) :: bdata CONTAINS ! !================================================================================ SUBROUTINE basic_data ! ! Define basic data ! use mpihelper IMPLICIT NONE ! ! Local vars and arrays CHARACTER(len=256) :: line, inputfilename ! NAMELIST /BASIC/ job_time, extra_time, nrun, tmax, dt, nlres, nlsave, newres, nlxg, & & nplasma, potinn, potout, B0, lz, n0, nz, nnr, femorder, ngauss, & & nlppform, plasmadim, radii, temp, Rcurv, width, it0d, it2d, itparts, ittext, & & resfile, rstfile, itgraph, nlPhis, distribtype, nblock, nlclassical, H0, P0, partperiodic, & & temprescale !________________________________________________________________________________ ! 1. Process Standard Input File ! IF(mpirank .eq. 0) THEN IF(COMMAND_ARGUMENT_COUNT().NE.1)THEN WRITE(*,*)'ERROR, ONE COMMAND-LINE ARGUMENT REQUIRED, STOPPING' STOP ENDIF CALL GET_COMMAND_ARGUMENT(1,inputfilename) OPEN(UNIT=lu_in,FILE=trim(inputfilename)) !________________________________________________________________________________ ! 1. Label the run ! READ(lu_in,'(a)') label1 READ(lu_in,'(a)') label2 READ(lu_in,'(a)') label3 READ(lu_in,'(a)') label4 ! WRITE(*,'(12x,a/)') label1(1:len_trim(label1)) WRITE(*,'(12x,a/)') label2(1:len_trim(label2)) WRITE(*,'(12x,a/)') label3(1:len_trim(label3)) WRITE(*,'(12x,a/)') label4(1:len_trim(label4)) !________________________________________________________________________________ ! 2. Read in basic data specific to run ! READ(lu_in,basic) WRITE(*,basic) bdata=BASICDATA( nlres, nlsave, newres, nlxg, nlppform, nlPhis, nlclassical, & & nplasma, nz, it0d, it2d, itparts, itgraph, distribtype, nblock, & & nrun, job_time, extra_time, tmax, dt, potinn, potout, B0, n0, temp, & & Rcurv, width, H0, P0, femorder, ngauss, nnr, lz, radii, plasmadim, resfile, partperiodic) END IF ! Distribute run parameters to all MPI workers CALL init_mpitypes ! initialize all mpi types that will be needed in the simulation CALL MPI_Bcast(bdata, 1, basicdata_type, 0, MPI_COMM_WORLD, ierr) CALL readbdata ! END SUBROUTINE basic_data !================================================================================ SUBROUTINE daytim(str) ! ! Print date and time ! IMPLICIT NONE ! CHARACTER(len=*), INTENT(in) :: str ! ! Local vars and arrays CHARACTER(len=16) :: d, t, dat, functime !________________________________________________________________________________ ! CALL DATE_AND_TIME(d,t) dat=d(7:8) // '/' // d(5:6) // '/' // d(1:4) functime=t(1:2) // ':' // t(3:4) // ':' // t(5:10) WRITE(*,'(a,1x,a,1x,a)') str, dat(1:10), functime(1:12) ! END SUBROUTINE daytim !================================================================================ SUBROUTINE timera(cntrl, str, eltime) ! ! Timers (cntrl=0/1 to Init/Update) ! IMPLICIT NONE INTEGER, INTENT(in) :: cntrl CHARACTER(len=*), INTENT(in) :: str DOUBLE PRECISION, OPTIONAL, INTENT(out) :: eltime ! INTEGER, PARAMETER :: ncmax=128 INTEGER, SAVE :: icall=0, nc=0 DOUBLE PRECISION, SAVE :: startt0=0.0 CHARACTER(len=16), SAVE :: which(ncmax) INTEGER :: lstr, found, i DOUBLE PRECISION :: seconds DOUBLE PRECISION, DIMENSION(ncmax), SAVE :: startt = 0.0, endt = 0.0 !________________________________________________________________________________ IF( icall .EQ. 0 ) THEN icall = icall+1 startt0 = seconds() END IF lstr = LEN_TRIM(str) IF( lstr .GT. 0 ) found = loc(str) !________________________________________________________________________________ ! SELECT CASE (cntrl) ! CASE(-1) ! Current wall time IF( PRESENT(eltime) ) THEN eltime = seconds() - startt0 ELSE WRITE(*,'(/a,a,1pe10.3/)') "++ ", ' Wall time used so far = ', seconds() - startt0 END IF ! CASE(0) ! Init Timer IF( found .EQ. 0 ) THEN ! Called for the 1st time for 'str' nc = nc+1 which(nc) = str(1:lstr) found = nc END IF startt(found) = seconds() ! CASE(1) ! Update timer endt(found) = seconds() - startt(found) IF( PRESENT(eltime) ) THEN eltime = endt(found) ELSE WRITE(*,'(/a,a,1pe10.3/)') "++ "//str, ' wall clock time = ', endt(found) END IF ! CASE(2) ! Update and reset timer endt(found) = endt(found) + seconds() - startt(found) startt(found) = seconds() IF( PRESENT(eltime) ) THEN eltime = endt(found) END IF ! CASE(9) ! Display all timers IF( nc .GT. 0 ) THEN WRITE(*,'(a)') "Timer Summary" WRITE(*,'(a)') "=============" DO i=1,nc WRITE(*,'(a20,2x,2(1pe12.3))') TRIM(which(i))//":", endt(i) END DO END IF ! END SELECT ! CONTAINS INTEGER FUNCTION loc(str) CHARACTER(len=*), INTENT(in) :: str INTEGER :: i, ind loc = 0 DO i=1,nc ind = INDEX(which(i), str(1:lstr)) IF( ind .GT. 0 .AND. LEN_TRIM(which(i)) .EQ. lstr ) THEN loc = i EXIT END IF END DO END FUNCTION loc END SUBROUTINE timera !================================================================================ SUBROUTINE readbdata nlres = bdata%nlres nlsave = bdata%nlsave newres = bdata%newres nlxg = bdata%nlxg nlppform = bdata%nlppform nlPhis = bdata%nlPhis nlclassical = bdata%nlclassical nplasma = bdata%nplasma nz = bdata%nz it0d = bdata%it0d it2d = bdata%it2d itparts = bdata%itparts itgraph = bdata%itgraph distribtype = bdata%distribtype nblock = bdata%nblock nrun = bdata%nrun job_time = bdata%job_time extra_time = bdata%extra_time tmax = bdata%tmax dt = bdata%dt potinn = bdata%potinn potout = bdata%potout B0 = bdata%B0 n0 = bdata%n0 temp = bdata%temp Rcurv = bdata%Rcurv width = bdata%width H0 = bdata%H0 P0 = bdata%P0 femorder = bdata%femorder ngauss = bdata%ngauss nnr = bdata%nnr lz = bdata%lz radii = bdata%radii plasmadim = bdata%plasmadim resfile = bdata%resfile partperiodic = bdata%partperiodic END SUBROUTINE readbdata !================================================================================ END MODULE basic diff --git a/src/diagnose.f90 b/src/diagnose.f90 index 8241134..fc37d5d 100644 --- a/src/diagnose.f90 +++ b/src/diagnose.f90 @@ -1,239 +1,240 @@ SUBROUTINE diagnose(kstep) ! ! Diagnostics ! USE basic USE futils USE hashtable USE beam, ONLY : parts, epot, ekin, etot, etot0, ireducerequest, Nplocs_all USE xg, ONLY : initw, updt_xg_var USE fields, ONLY: phi_spline USE mpi IMPLICIT NONE ! INTEGER, INTENT(in) :: kstep ! ! Local vars and arrays INTEGER, PARAMETER :: BUFSIZE = 20 CHARACTER(len=128) :: str, fname ! Only process 0 should save on file IF(mpirank .ne. 0) RETURN !________________________________________________________________________________ ! 1. Initial diagnostics IF( kstep .EQ. 0) THEN ! WRITE(*,'(a)') ' Initial diagnostics' ! ! 1.1 Initial run or when NEWRES set to .TRUE. ! IF( .NOT. nlres .OR. newres) THEN CALL creatf(resfile, fidres, 'Espic2d result file', 'd') WRITE(*,'(3x,a,a)') TRIM(resfile), ' created' ! ! Label the run IF( LEN_TRIM(label1).GT.0 ) CALL attach(fidres, "/", "label1", TRIM(label1)) IF( LEN_TRIM(label2).GT.0 ) CALL attach(fidres, "/", "label2", TRIM(label2)) IF( LEN_TRIM(label3).GT.0 ) CALL attach(fidres, "/", "label3", TRIM(label3)) IF( LEN_TRIM(label4).GT.0 ) CALL attach(fidres, "/", "label4", TRIM(label4)) ! ! Job number jobnum = 0 ! ! Data group CALL creatg(fidres, "/data", "data") CALL creatg(fidres, "/data/var0d", "0d history arrays") CALL creatg(fidres, "/data/var1d","1d history arrays") CALL creatg(fidres, "/data/part", "part phase space") CALL creatg(fidres, "/data/fields", "Electro static potential and Er Ez fields") ! ! File group CALL creatg(fidres, "/files", "files") CALL attach(fidres, "/files", "jobnum", jobnum) CALL putarr(fidres, "/data/var1d/rgrid", rgrid) CALL putarr(fidres, "/data/var1d/zgrid", zgrid) ! Part and fields vectors ! Initialize time-dependent particle variables CALL creatd(fidres, 1, SHAPE(parts%R), "/data/part/R", "R") CALL creatd(fidres, 1, SHAPE(parts%Z), "/data/part/Z", "Z") CALL creatd(fidres, 1, SHAPE(parts%Z), "/data/part/THET", "THET") CALL creatd(fidres, 1, SHAPE(parts%UZ), "/data/part/UR", "UZ") CALL creatd(fidres, 1, SHAPE(parts%UR), "/data/part/UZ", "UR") CALL creatd(fidres, 1, SHAPE(parts%UTHET), "/data/part/UTHET", "UTHET") CALL creatd(fidres, 1, SHAPE(parts%Rindex), "/data/part/Rindex", "Rindex") CALL creatd(fidres, 1, SHAPE(parts%Zindex), "/data/part/Zindex", "Zindex") CALL creatd(fidres, 1, SHAPE(parts%partindex), "/data/part/partindex", "partindex") CALL creatd(fidres, 1, SHAPE(parts%pot), "/data/part/pot", "pot") CALL creatd(fidres, 0, SHAPE(time), "/data/part/time", "time" ) CALL creatd(fidres, 1, SHAPE(pot), "/data/fields/pot", "pot") CALL creatd(fidres, 1, SHAPE(Er), "/data/fields/Er", "Er") CALL creatd(fidres, 1, SHAPE(Ez), "/data/fields/Ez", "Ez") CALL creatd(fidres, 1, SHAPE(phi_spline), "/data/fields/phi", "spline form of Phi") CALL creatd(fidres, 2, SHAPE(moments), "/data/fields/moments", "moments") !CALL creatd(fidres, 1, SHAPE(partdensity), "/data/fields/partdensity", "partdensity") !CALL creatd(fidres, 1, SHAPE(fluidur), "/data/fields/fluidur", "fluidur") !CALL creatd(fidres, 1, SHAPE(fluiduz), "/data/fields/fluiduz", "fluiduz") !CALL creatd(fidres, 1, SHAPE(fluiduthet), "/data/fields/fluiduthet", "fluiduthet") !CALL creatd(fidres, 2, SHAPE(Presstens), "/data/fields/presstens", "presstens") CALL creatd(fidres, 0, SHAPE(time), "/data/fields/time", "time" ) IF(mpisize .gt. 1) CALL creatd(fidres, 1, SHAPE(Nplocs_all), "/data/part/Nplocs_all", "Nplocs_all") CALL putarr(fidres, "/data/fields/Br", Br) CALL putarr(fidres, "/data/fields/Bz", Bz) CALL putarr(fidres, "/data/fields/Athet", Athet) + CALL putarr(fidres, "/data/fields/volume", Volume) ! ! 1.2 Restart run ! ELSE CALL cp2bk(resfile) ! backup previous result file CALL openf(resfile, fidres,'rw', 'd') WRITE(*,'(3x,a,a)') TRIM(resfile), ' open' CALL getatt(fidres, "/files", "jobnum", jobnum) jobnum = jobnum+1 WRITE(*,'(3x,a,i3)') "Current Job Number =", jobnum CALL attach(fidres, "/files", "jobnum", jobnum) END IF ! ! Add input namelist variables as attributes of /data/input WRITE(str,'(a,i2.2)') "/data/input.",jobnum CALL creatg(fidres, TRIM(str)) CALL attach(fidres, TRIM(str), "job_time", job_time) CALL attach(fidres, TRIM(str), "extra_time", extra_time) CALL attach(fidres, TRIM(str), "dt", dt) CALL attach(fidres, TRIM(str), "tmax", tmax) CALL attach(fidres, TRIM(str), "nrun", nrun) CALL attach(fidres, TRIM(str), "nlres", nlres) CALL attach(fidres, TRIM(str), "nlsave", nlsave) CALL attach(fidres, TRIM(str), "newres", newres) CALL attach(fidres, TRIM(str), "nz", nz) CALL putarr(fidres, TRIM(str)//"/lz", lz) CALL attach(fidres, TRIM(str), "nplasma", nplasma) CALL attach(fidres, TRIM(str), "potinn", potinn) CALL attach(fidres, TRIM(str), "potout", potout) CALL attach(fidres, TRIM(str), "B0", B0) CALL attach(fidres, TRIM(str), "Rcurv", Rcurv) CALL attach(fidres, TRIM(str), "width", width) CALL attach(fidres, TRIM(str), "n0", n0) CALL attach(fidres, TRIM(str), "temp", temp) CALL attach(fidres, TRIM(str), "it0d", it0d) CALL attach(fidres, TRIM(str), "it2d", it2d) CALL attach(fidres, TRIM(str), "itparts", itparts) CALL attach(fidres, TRIM(str), "nlclassical", nlclassical) CALL attach(fidres, TRIM(str), "nlPhis", nlPhis) CALL attach(fidres, TRIM(str), "qsim", qsim) CALL attach(fidres, TRIM(str), "msim", msim) CALL attach(fidres, TRIM(str), "V", V) CALL attach(fidres, TRIM(str), "startstep", cstep) CALL putarr(fidres, TRIM(str)//"/femorder", femorder) CALL putarr(fidres, TRIM(str)//"/ngauss", ngauss) CALL putarr(fidres, TRIM(str)//"/nnr", nnr) CALL putarr(fidres, TRIM(str)//"/radii", radii) CALL putarr(fidres, TRIM(str)//"/plasmadim", plasmadim) ! Save STDIN of this run WRITE(str,'(a,i2.2)') "/files/STDIN.",jobnum INQUIRE(unit=lu_in, name=fname) CALL putfile(fidres, TRIM(str), TRIM(fname)) CLOSE(lu_in) ! ! Initialize buffers for 0d history arrays CALL htable_init(hbuf0, BUFSIZE) CALL set_htable_fileid(hbuf0, fidres, "/data/var0d") ! ! ! Initialize Xgrafix ! IF(nlxg) THEN CALL initw END IF END IF !________________________________________________________________________________ ! 2. Periodic diagnostics IF( kstep .NE. -1) THEN ! IF(modulo(step,ittext).eq. 0 .or. nlend) THEN WRITE(*,'(a,1x,i8.8,a1,i8.8,20x,a,1pe10.3)') & & '*** Timestep (this run/total) =', step, '/', cstep, 'Time =', time WRITE(*,'(a,4(1pe12.4))') 'Epot, Ekin, Etot, Eerr', epot, ekin, etot, etot-etot0 END IF !________________________________________________________________________________ ! ! 2.1 0d history arrays ! IF(modulo(step,it0d).eq. 0 .or. nlend) THEN CALL add_record(hbuf0, "time", "simulation time", time) CALL add_record(hbuf0, "epot", "potential energy", epot) CALL add_record(hbuf0, "ekin", "kinetic energy", ekin) CALL add_record(hbuf0, "etot", "total energy", etot) CALL add_record(hbuf0, "etot0", "theoretical total energy", etot0) IF (mpisize .gt. 1) THEN CALL add_record(hbuf0, "nbparts", "number of remaining parts in the simulation space", REAL(sum(Nplocs_all),kind=db)) ELSE CALL add_record(hbuf0, "nbparts", "number of remaining parts in the simulation space", REAL(parts%Nptot,kind=db)) END IF CALL htable_endstep(hbuf0) END IF ! ! 2.2 1d profiles IF(modulo(step,it2d).eq. 0 .or. nlend) THEN CALL append(fidres, "/data/fields/time", time) CALL append(fidres, "/data/fields/pot", pot) CALL append(fidres, "/data/fields/Er", Er) CALL append(fidres, "/data/fields/Ez", Ez) CALL append(fidres, "/data/fields/phi", phi_spline) CALL append(fidres, "/data/fields/moments", moments) !CALL append(fidres, "/data/fields/partdensity", partdensity) !CALL append(fidres, "/data/fields/fluidur", fluidur) !CALL append(fidres, "/data/fields/fluiduz", fluiduz) !CALL append(fidres, "/data/fields/fluiduthet", fluiduthet) !CALL append(fidres, "/data/fields/presstens", Presstens) END IF ! ! 2.3 2d profiles IF(modulo(step,itparts).eq. 0 .or. nlend) THEN !PRINT*, 'write particles to file_____________________' CALL append(fidres, "/data/part/time", time) CALL append(fidres, "/data/part/R", parts%R) CALL append(fidres, "/data/part/Z", parts%Z) CALL append(fidres, "/data/part/THET", parts%THET) CALL append(fidres, "/data/part/UZ", parts%UZ) CALL append(fidres, "/data/part/UR", parts%UR) CALL append(fidres, "/data/part/UTHET", parts%UTHET) CALL append(fidres, "/data/part/pot", parts%pot) CALL append(fidres, "/data/part/Rindex", REAL(parts%Rindex,kind=db)) CALL append(fidres, "/data/part/Zindex", REAL(parts%Zindex,kind=db)) CALL append(fidres, "/data/part/partindex", REAL(parts%partindex,kind=db)) IF(mpisize .gt. 1) CALL append(fidres, "/data/part/Nplocs_all", REAL(Nplocs_all,kind=db)) ! END IF ! ! 2.4 3d profiles ! ! ! 2.5 Xgrafix ! IF(nlxg .AND. modulo(kstep,itgraph) .eq. 0) THEN call xgevent CALL updt_xg_var CALL xgupdate END IF !________________________________________________________________________________ ! 3. Final diagnostics ELSE ! ! Flush 0d history array buffers CALL htable_hdf5_flush(hbuf0) ! ! Close all diagnostic files CALL closef(fidres) !________________________________________________________________________________ END IF ! END SUBROUTINE diagnose diff --git a/src/fields_mod.f90 b/src/fields_mod.f90 index e0081dd..fef560f 100644 --- a/src/fields_mod.f90 +++ b/src/fields_mod.f90 @@ -1,533 +1,590 @@ !------------------------------------------------------------------------------ ! EPFL/Swiss Plasma Center !------------------------------------------------------------------------------ ! ! MODULE: beam ! !> @author !> Patryk Kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> Module responsible for initializing the magnetic field and solve the Poisson equation !------------------------------------------------------------------------------ MODULE fields USE constants USE basic, ONLY: nr, nz, zgrid, rgrid, Br, Bz, Er, Ez, femorder, ngauss, nlppform, pot, Athet, splrz, nlperiod, phinorm, nlPhis, nrank, mpirank, mpisize, step, it2d USE beam, ONLY: parts USE bsplines USE mumps_bsplines USE mpi IMPLICIT NONE DOUBLE PRECISION, allocatable, SAVE :: matcoef(:,:), phi_spline(:), vec1(:), vec2(:) TYPE(mumps_mat), SAVE :: femat ! Finite Element Method matrix CONTAINS !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Set-up the necessary variables for solving Poisson and computes the magnetic field on the grid ! !--------------------------------------------------------------------------- SUBROUTINE fields_init USE basic, ONLY: pot, nlperiod, nrank, rhs, moments, leftproc, rightproc, Zbounds USE bsplines USE mumps_bsplines USE mpihelper INTEGER :: nrz(2) ! Set up 2d spline splrz used in the FEM CALL set_spline(femorder, ngauss, zgrid, rgrid, splrz, nlppform=nlppform, period=nlperiod) ! Calculate dimension of splines nrz(1)=nz nrz(2)=nr CALL get_dim(splrz, nrank, nrz, femorder) ALLOCATE(matcoef(nrank(1),nrank(2))) ALLOCATE(pot((nr+1)*(nz+1))) ALLOCATE(rhs(nrank(1)*nrank(2))) ALLOCATE(moments(10,nrank(1)*nrank(2))) ALLOCATE(phi_spline(nrank(1)*nrank(2))) + ALLOCATE(volume(nrank(1)*nrank(2))) ALLOCATE(Br((nr+1)*(nz+1)),Bz((nr+1)*(nz+1))) ALLOCATE(Athet((nr+1)*(nz+1))) ALLOCATE(Er((nr+1)*(nz+1)),Ez((nr+1)*(nz+1))) ! Calculate magnetic field mirror components in grid points (Davidson analytical formula employed) CALL magnet ! Calculate and factorise FEM matrix (depends only on mesh) CALL fematrix CALL factor(femat) + CALL comp_volume IF(mpisize .gt. 1) THEN CALL init_overlaps(nrank, femorder, Zbounds(mpirank+1), moments, mpirank, leftproc, rightproc) END IF END SUBROUTINE fields_init !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Construct the right hand side vector used in the FEM Poisson solver ! !--------------------------------------------------------------------------- SUBROUTINE rhscon USE bsplines USE mumps_bsplines USE mpi USE basic, ONLY: omegap, omegac, V, potinn, potout, ierr, nrank, rhs, nplasma, moments, nlend, step, it2d USE beam, ONLY: parts USE mpihelper INTEGER ::irow, jcol, it, jw, mu, i, k, iend, nbunch INTEGER,DIMENSION(:),ALLOCATABLE::zleft, rleft DOUBLE PRECISION :: vr, vthet, vz, coeff, chargecoeff DOUBLE PRECISION, ALLOCATABLE :: fun(:,:,:), fun2(:,:,:) nbunch=64 ! Particle bunch size used when calling basfun ALLOCATE(zleft(nbunch),rleft(nbunch)) ALLOCATE(fun(1:femorder(1)+1,0:0,nbunch),fun2(1:femorder(2)+1,0:0,nbunch)) ! Arrays keeping values of b-splines at gauss node moments=0 ! Reset the moments matrix rhs=0 ! reset rhs chargecoeff=0.5/pi*omegap/omegac*V/nplasma ! Normalized charge density simulated by each macro particle ! Assemble rhs IF (parts%Nptot .ne. 0) THEN !$OMP PARALLEL DO DEFAULT(SHARED), PRIVATE(i, zleft, rleft, jw, it, iend, irow, jcol, mu, k, vr, vz, vthet, coeff) PRIVATE(fun, fun2) DO i=1,parts%Nptot,nbunch ! Avoid segmentation fault by accessing non relevant data iend=min(i+nbunch-1,parts%Nptot) ! Localize the particle CALL locintv(splrz%sp1, parts%Z(i:iend), zleft) CALL locintv(splrz%sp2, parts%R(i:iend), rleft) ! Compute the value of the splines at the particles positions CALL basfun(parts%Z(i:iend), splrz%sp1, fun, zleft+1) CALL basfun(parts%R(i:iend), splrz%sp2, fun2, rleft+1) IF(modulo(step,it2d).eq. 0 .or. nlend) THEN ! We want to compute the higher order moments !$OMP SIMD DO k=1,(iend-i+1) DO jw=1,(femorder(2)+1) DO it=1,(femorder(1)+1) irow=zleft(k)+it jcol=rleft(k)+jw mu=irow+(jcol-1)*(nrank(1)) coeff=fun(it,0,k)*fun2(jw,0,k) ! Add contribution of particle nbunch to rhs grid point mu vr=parts%UR(i+k-1)/parts%Gamma(i+k-1) vz=parts%UZ(i+k-1)/parts%Gamma(i+k-1) vthet=parts%UTHET(i+k-1)/parts%Gamma(i+k-1) !$OMP ATOMIC UPDATE moments(1, mu) = moments(1, mu) + coeff !$OMP END ATOMIC !$OMP ATOMIC UPDATE moments(2, mu) = moments(2, mu) + coeff*vr !$OMP END ATOMIC !$OMP ATOMIC UPDATE moments(3, mu) = moments(3, mu) + coeff*vthet !$OMP END ATOMIC !$OMP ATOMIC UPDATE moments(4, mu) = moments(4, mu) + coeff*vz !$OMP END ATOMIC !$OMP ATOMIC UPDATE moments(5, mu) = moments(5, mu) + coeff*vr*vr !$OMP END ATOMIC !$OMP ATOMIC UPDATE moments(6, mu) = moments(6, mu) + coeff*vr*vthet !$OMP END ATOMIC !$OMP ATOMIC UPDATE moments(7, mu) = moments(7, mu) + coeff*vr*vz !$OMP END ATOMIC !$OMP ATOMIC UPDATE moments(8, mu) = moments(8, mu) + coeff*vthet*vthet !$OMP END ATOMIC !$OMP ATOMIC UPDATE moments(9, mu) = moments(9, mu) + coeff*vthet*vz !$OMP END ATOMIC !$OMP ATOMIC UPDATE moments(10,mu) = moments(10,mu) + coeff*vz*vz !$OMP END ATOMIC END DO END DO END DO !$OMP END SIMD ELSE !$OMP SIMD DO k=1,(iend-i+1) DO jw=1,(femorder(2)+1) DO it=1,(femorder(1)+1) irow=zleft(k)+it jcol=rleft(k)+jw mu=irow+(jcol-1)*(nrank(1)) ! Add contribution of particle nbunch to rhs grid point mu !$OMP ATOMIC update moments(1, mu) = moments(1, mu) + fun(it,0,k)*fun2(jw,0,k) !$OMP END ATOMIC END DO END DO END DO !$OMP END SIMD END IF END DO !$OMP END PARALLEL DO END IF ! If we are using MPI parallelism, reduce the rhs on the root process IF(mpisize .gt. 1) THEN CALL fields_comm END IF IF(nlphis) THEN ! We consider self-consistent interactions rhs=moments(1,1:nrank(1)*nrank(2))*chargecoeff ! Normalized charge density simulated by each macro particle ELSE rhs=0 END IF ! Apply the Dirichlet boundary conditions on the coaxial insert (if present) IF(rgrid(0).ne.0.0_db) rhs(1:nrank(1))=potinn ! Apply the Dirichlet boundary conditions on the cylinder walls rhs((nrank(2)-1)*nrank(1)+1:nrank(2)*nrank(1))=potout DEALLOCATE(fun,fun2, zleft, rleft) END SUBROUTINE rhscon SUBROUTINE fields_comm USE mpihelper USE Basic, ONLY: moments, Zbounds, mpirank, nlend, it2d, step, leftproc, rightproc INTEGER:: ierr, i, j INTEGER:: displs(mpisize), counts(mpisize) INTEGER:: overlap_type displs=Zbounds(0:mpisize-1) counts=Zbounds(1:mpisize)-Zbounds(0:mpisize-1) counts(mpisize)=counts(mpisize)+femorder(1) IF(modulo(step,it2d).eq. 0 .or. nlend) THEN CALL start_momentscomm(mpirank, leftproc, rightproc, moments, nrank, femorder, Zbounds(mpirank+1)) IF(leftproc .lt. mpirank) THEN !$OMP PARALLEL DO SIMD DEFAULT(SHARED), PRIVATE(i) DO j=1,femorder(1) DO i=1,nrank(2) moments(1:10,Zbounds(mpirank)+(i-1)*nrank(1)+j) = moments(1:10,Zbounds(mpirank)+(i-1)*nrank(1)+j)& & + momentsoverlap_buffer(10*(nrank(2)*(j-1)+i-1)+1:10*(nrank(2)*(j-1)+i)) END DO END DO !$OMP END PARALLEL DO SIMD END IF ! Set communication vector type overlap_type=momentsoverlap_type ELSE CALL start_rhscomm(mpirank, leftproc, rightproc, moments, nrank, femorder, Zbounds(mpirank+1)) IF(leftproc .lt. mpirank) THEN !$OMP PARALLEL DO SIMD DEFAULT(SHARED), PRIVATE(i) DO j=1,femorder(1) DO i=1,nrank(2) moments(1, Zbounds(mpirank)+(i-1)*nrank(1)+j) = moments(1,Zbounds(mpirank)+(i-1)*nrank(1)+j)& & + rhsoverlap_buffer(nrank(2)*(j-1)+i) END DO END DO !$OMP END PARALLEL DO SIMD END IF ! Set communication vector type overlap_type=rhsoverlap_type END IF IF(mpirank .eq. 0) THEN CALL MPI_GATHERV(MPI_IN_PLACE, counts(mpirank+1), overlap_type, & & moments, counts, displs, overlap_type, 0, MPI_COMM_WORLD, ierr) ELSE CALL MPI_GATHERV(moments(1, displs(mpirank+1)+1), counts(mpirank+1), overlap_type, & & moments, counts, displs, overlap_type, 0, MPI_COMM_WORLD, ierr) END IF END SUBROUTINE fields_comm !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Solves Poisson equation using FEM. Distributes the result on all MPI workers and interpolate the electric forces !> for each particle. ! !--------------------------------------------------------------------------- SUBROUTINE poisson USE basic, ONLY: rhs, nrank, pot USE bsplines USE mumps_bsplines USE futils INTEGER:: partend, i, ierr, jder(2), iend INTEGER:: nbunch=64 jder(1)=0 jder(2)=0 partend=parts%Nptot ! Only the root process solves Poisson IF(mpirank.eq.0) CALL bsolve(femat,rhs,phi_spline) ! Broadcast the solution to all MPI workers CALL MPI_Bcast(phi_spline(1), nrank(1)*nrank(2), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) matcoef= reshape(phi_spline, (/nrank(1),nrank(2)/)) ! Evaluate the electric potential on the grid points CALL gridval(splrz, vec1, vec2, pot, jder, matcoef) ! Evaluate the electric potential and field at the particles position !$OMP PARALLEL DO SIMD DEFAULT(SHARED), PRIVATE(iend) DO i=1,parts%Nptot,nbunch ! Avoid segmentation fault by accessing non relevant data iend=min(i+nbunch-1,parts%Nptot) CALL getgrad(splrz, parts%Z(i:iend), parts%R(i:iend), parts%pot(i:iend), parts%EZ(i:iend), parts%ER(i:iend)) parts%EZ(i:iend)=-parts%Ez(i:iend) parts%ER(i:iend)=-parts%ER(i:iend) END DO !$OMP END PARALLEL DO SIMD IF( mpirank .eq. 0 .and. modulo(step,it2d) .eq. 0) THEN ! On the root process, compute the electric field for diagnostic purposes CALL gridval(splrz, vec1, vec2, Ez, (/1,0/)) CALL gridval(splrz, vec1, vec2, Er, (/0,1/)) Ez=-Ez Er=-Er END IF END SUBROUTINE poisson !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Constucts the FEM matrix using bsplines initialized in fields_init !--------------------------------------------------------------------------- SUBROUTINE fematrix USE bsplines USE mumps_bsplines DOUBLE PRECISION, ALLOCATABLE :: xgauss(:,:), wgauss(:), zg(:), rg(:), wzg(:), wrg(:) DOUBLE PRECISION, ALLOCATABLE :: f(:,:), aux(:) DOUBLE PRECISION, ALLOCATABLE :: coefs(:) DOUBLE PRECISION, ALLOCATABLE :: fun(:,:), fun2(:,:) DOUBLE PRECISION :: contrib INTEGER, ALLOCATABLE :: idert(:), iderw(:) INTEGER :: i, j, k, jt, iw, irow, jcol, mu, igauss, iterm, irow2, jcol2, mu2, kterms kterms=2 !Number of terms in weak form ALLOCATE(fun(1:femorder(1)+1,0:1),fun2(1:femorder(2)+1,0:1))!Arrays keeping values of b-splines at gauss node ALLOCATE(xgauss(ngauss(1)*ngauss(2),2), wgauss(ngauss(1)*ngauss(2)),zg(ngauss(1)),rg(ngauss(2)), wzg(ngauss(1)), wrg(ngauss(2))) !Gaussian nodes and weights arrays ALLOCATE(f((femorder(1)+1)*(femorder(2)+1),2),aux(femorder(1)+1)) !Auxiliary arrays ordering bsplines ALLOCATE(idert(kterms), iderw(kterms), coefs(kterms)) !Pointers on the order of derivatives ! Constuction of auxiliary array ordering bsplines in given interval DO i=1,(femorder(1)+1) aux(i)=i END DO DO i=1,(femorder(2)+1) f((i-1)*(femorder(1)+1)+1:i*(femorder(1)+1),1)=aux f((i-1)*(femorder(1)+1)+1:i*(femorder(1)+1),2)=i END DO ! Initialisation of FEM matrix CALL init(nrank(1)*nrank(2),2,femat) ! Assemble FEM matrix DO j=1,nr ! Loop on r position DO i=1,nz ! Loop on z position ! Computation of gauss weight and position in r and z direction for gaussian integration CALL get_gauss(splrz%sp1, ngauss(1), i, zg, wzg) CALL get_gauss(splrz%sp2, ngauss(2), j, rg, wrg) ! Construction of matrix xgauss and wgauss storing the weight and position for 2d gaussian integration DO k=1,ngauss(2) xgauss((k-1)*ngauss(1)+1:k*ngauss(1),1)=zg xgauss((k-1)*ngauss(1)+1:k*ngauss(1),2)=rg(k) wgauss((k-1)*ngauss(1)+1:k*ngauss(1))=wrg(k)*wzg END DO DO igauss=1,ngauss(1)*ngauss(2) ! Loop on gaussian weights and positions CALL basfun(xgauss(igauss,1), splrz%sp1, fun, i) CALL basfun(xgauss(igauss,2), splrz%sp2, fun2, j) CALL coefeq(xgauss(igauss,:), idert, iderw, coefs) DO iterm=1,kterms ! Loop on the two integration dimensions DO jt=1,(1+femorder(1))*(femorder(2)+1) DO iw=1,(1+femorder(1))*(femorder(2)+1) irow=i+f(jt,1)-1; jcol=j+f(jt,2)-1 irow2=i+f(iw,1)-1; jcol2=j+f(iw,2)-1 mu=irow+(jcol-1)*nrank(1) mu2=irow2+(jcol2-1)*nrank(1) contrib=fun(f(jt,1),idert(iterm))* fun(f(iw,1),idert(iterm))* & & fun2(f(jt,2),iderw(iterm))*fun2(f(iw,2),iderw(iterm))* & & wgauss(igauss)*xgauss(igauss,2) CALL updtmat(femat, mu, mu2, contrib) END DO END DO END DO END DO END DO END DO ! Impose Dirichlet boundary conditions on FEM matrix CALL fe_dirichlet DEALLOCATE(xgauss, wgauss,zg,rg, wzg, wrg) - DEALLOCATE(f, aux) + DEALLOCATE(f, aux) DEALLOCATE(idert, iderw, coefs, fun,fun2) END SUBROUTINE fematrix + SUBROUTINE comp_volume + + USE bsplines + USE mumps_bsplines + + DOUBLE PRECISION, ALLOCATABLE :: xgauss(:,:), wgauss(:), zg(:), rg(:), wzg(:), wrg(:) + DOUBLE PRECISION, ALLOCATABLE :: f(:,:), aux(:) + DOUBLE PRECISION, ALLOCATABLE :: coefs(:) + DOUBLE PRECISION, ALLOCATABLE :: fun(:,:), fun2(:,:) + INTEGER :: i, j, k, jt, irow, jcol, mu, igauss + + ALLOCATE(fun(1:femorder(1)+1,0:0),fun2(1:femorder(2)+1,0:0))!Arrays keeping values of b-splines at gauss node + ALLOCATE(xgauss(ngauss(1)*ngauss(2),2), wgauss(ngauss(1)*ngauss(2)),zg(ngauss(1)),rg(ngauss(2)), wzg(ngauss(1)), wrg(ngauss(2))) !Gaussian nodes and weights arrays + ALLOCATE(f((femorder(1)+1)*(femorder(2)+1),2),aux(femorder(1)+1)) !Auxiliary arrays ordering bsplines + +! Constuction of auxiliary array ordering bsplines in given interval + DO i=1,(femorder(1)+1) + aux(i)=i + END DO + DO i=1,(femorder(2)+1) + f((i-1)*(femorder(1)+1)+1:i*(femorder(1)+1),1)=aux + f((i-1)*(femorder(1)+1)+1:i*(femorder(1)+1),2)=i + END DO + +! Assemble Volume matrix + DO j=1,nr ! Loop on r position + DO i=1,nz ! Loop on z position + ! Computation of gauss weight and position in r and z direction for gaussian integration + CALL get_gauss(splrz%sp1, ngauss(1), i, zg, wzg) + CALL get_gauss(splrz%sp2, ngauss(2), j, rg, wrg) + ! Construction of matrix xgauss and wgauss storing the weight and position for 2d gaussian integration + DO k=1,ngauss(2) + xgauss((k-1)*ngauss(1)+1:k*ngauss(1),1)=zg + xgauss((k-1)*ngauss(1)+1:k*ngauss(1),2)=rg(k) + wgauss((k-1)*ngauss(1)+1:k*ngauss(1))=wrg(k)*wzg + END DO + DO igauss=1,ngauss(1)*ngauss(2) ! Loop on gaussian weights and positions + CALL basfun(xgauss(igauss,1), splrz%sp1, fun, i) + CALL basfun(xgauss(igauss,2), splrz%sp2, fun2, j) + DO jt=1,(1+femorder(1))*(femorder(2)+1) + irow=i+f(jt,1)-1; jcol=j+f(jt,2)-1 + mu=irow+(jcol-1)*nrank(1) + volume(mu)=volume(mu)+2*pi*fun(f(jt,1),0) * fun2(f(jt,2),0)* wgauss(igauss)*xgauss(igauss,2) + END DO + END DO + END DO + END DO + + + DEALLOCATE(xgauss, wgauss,zg,rg, wzg, wrg) + DEALLOCATE(f, aux) + DEALLOCATE(fun,fun2) + + END SUBROUTINE comp_volume + !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Imposes the dirichlet boundary conditions on the FEM matrix. !--------------------------------------------------------------------------- SUBROUTINE fe_dirichlet DOUBLE PRECISION, ALLOCATABLE :: arr(:) INTEGER :: i ALLOCATE(arr(nrank(1)*nrank(2))) DO i=1,nrank(1)*nrank(2) CALL getrow(femat,i,arr) END DO DO i=1,nrank(1) IF(rgrid(0).ne.0.0_db) THEN arr=0; arr(i)=1; CALL putrow(femat,i,arr) END IF arr=0; arr(nrank(1)*nrank(2)+1-i)=1 ; CALL putrow(femat,nrank(1)*nrank(2)+1-i,arr) END DO DEALLOCATE(arr) END SUBROUTINE fe_dirichlet !________________________________________________________________________________ SUBROUTINE coefeq(x, idt, idw, c) DOUBLE PRECISION, INTENT(in) :: x(2) INTEGER, INTENT(out) :: idt(:), idw(:) DOUBLE PRECISION, INTENT(out) :: c(SIZE(idt)) c(1) = x(2) idt(1) = 0 idw(1) = 1 c(2) = x(2) idt(2) = 1 idw(2) = 0 END SUBROUTINE coefeq !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Computes the magnetic field on the grid according to a magnetic mirror. !--------------------------------------------------------------------------- SUBROUTINE magnet USE basic, ONLY: B0, Rcurv, rgrid, zgrid, width, rnorm, nr, nz, bnorm USE constants, ONLY: Pi DOUBLE PRECISION :: rg, zg,halfLz, MirrorRatio INTEGER :: i, rindex 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)/bnorm Bz(i)=B0*(1-MirrorRatio*COS(2*pi*zg/width*rnorm)*bessi0(2*pi*rg/width*rnorm))/bnorm 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 magnet !________________________________________________________________________________ !Modified Bessel functions of the first kind of the zero order FUNCTION bessi0(x) DOUBLE PRECISION :: bessi0,x DOUBLE PRECISION :: ax DOUBLE PRECISION 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)))))))) endif return END FUNCTION bessi0 !________________________________________________________________________________ !Modified Bessel functions of the first kind of the first order FUNCTION bessi1(x) DOUBLE PRECISION :: bessi1,x DOUBLE PRECISION :: ax DOUBLE PRECISION 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.75) then y=(x/3.75)**2 bessi1=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))) else ax=abs(x) y=3.75/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 endif return END FUNCTION bessi1 !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Free the memory used by the fields module !--------------------------------------------------------------------------- SUBROUTINE clean_fields Use bsplines USE basic, ONLY: rhs, moments DEALLOCATE(matcoef) DEALLOCATE(pot) DEALLOCATE(rhs) DEALLOCATE(phi_spline) DEALLOCATE(moments) DEALLOCATE(Br,Bz) DEALLOCATE(Er,Ez) DEALLOCATE(vec1,vec2) Call DESTROY_SP(splrz) END SUBROUTINE clean_fields END MODULE fields