MODULE korc_simple_equilibrium_pdf USE korc_types USE korc_constants USE korc_HDF5 USE korc_hpc USE special_functions IMPLICIT NONE TYPE, PRIVATE :: PARAMS REAL(rp) :: E ! Parallel electric field normalized using the critical electric field REAL(rp) :: Zeff ! Effective atomic number of impurities REAL(rp) :: max_pitch_angle ! Maximum pitch angle of sampled PDF in degrees REAL(rp) :: min_pitch_angle ! Minimum pitch angle of sampled PDF in degrees REAL(rp) :: po ! Momentum of sampled PDF in units of mc REAL(rp) :: Bo REAL(rp) :: lambda END TYPE PARAMS TYPE(PARAMS), PRIVATE :: pdf_params REAL(rp), PRIVATE, PARAMETER :: xo = (C_ME*C_C**2/C_E)/1.0E6 REAL(rp), PRIVATE, PARAMETER :: Tol = 1.0E-5_rp REAL(rp), PRIVATE, PARAMETER :: minmax_buffer_size = 10.0_rp PUBLIC :: get_equilibrium_distribution PRIVATE :: initialize_params,& save_params,& sample_distribution,& deg2rad,& rad2deg,& fRE,& random_norm,& PR,& P_integral,& IntK,& IntBesselK CONTAINS SUBROUTINE get_equilibrium_distribution(params,eta,go,etao) TYPE(KORC_PARAMS), INTENT(IN) :: params REAL(rp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: eta REAL(rp), INTENT(IN) :: go REAL(rp), INTENT(OUT) :: etao call initialize_params(params,go) call save_params(params) call sample_distribution(params,eta,etao) END SUBROUTINE get_equilibrium_distribution SUBROUTINE initialize_params(params,go) TYPE(KORC_PARAMS), INTENT(IN) :: params REAL(rp), INTENT(IN) :: go REAL(rp) :: max_pitch_angle REAL(rp) :: min_pitch_angle REAL(rp) :: Zeff REAL(rp) :: E REAL(rp) :: Bo REAL(rp) :: lambda NAMELIST /SimpleEquilibriumPDF/ max_pitch_angle,min_pitch_angle,Zeff,E,Bo,lambda open(unit=default_unit_open,file=TRIM(params%path_to_inputs),status='OLD',form='formatted') read(default_unit_open,nml=SimpleEquilibriumPDF) close(default_unit_open) pdf_params%max_pitch_angle = max_pitch_angle pdf_params%min_pitch_angle = min_pitch_angle pdf_params%Zeff = Zeff pdf_params%E = E pdf_params%Bo = Bo pdf_params%lambda = lambda pdf_params%po = sqrt(go**2 - 1.0_rp) END SUBROUTINE initialize_params FUNCTION deg2rad(x) REAL(rp), INTENT(IN) :: x REAL(rp) :: deg2rad deg2rad = C_PI*x/180.0_rp END FUNCTION FUNCTION rad2deg(x) REAL(rp), INTENT(IN) :: x REAL(rp) :: rad2deg rad2deg = 180.0_rp*x/C_PI END FUNCTION FUNCTION fRE(eta,p) REAL(rp), INTENT(IN) :: eta ! pitch angle in degrees REAL(rp), INTENT(IN) :: p ! momentum in units of mc REAL(rp) :: fRE REAL(rp) :: A A = (2.0_rp*pdf_params%E/(pdf_params%Zeff + 1.0_rp))*(p**2/SQRT(p**2.0_rp + 1.0_rp)) fRE = 0.5_rp*A*EXP(A*COS(deg2rad(eta)))/SINH(A) ! fRE = fRE*PR(eta,p,pdf_params%Bo,pdf_params%lambda) END FUNCTION fRE FUNCTION random_norm(mean,sigma) REAL(rp), INTENT(IN) :: mean REAL(rp), INTENT(IN) :: sigma REAL(rp) :: random_norm REAL(rp) :: rand1, rand2 call RANDOM_NUMBER(rand1) call RANDOM_NUMBER(rand2) random_norm = SQRT(-2.0_rp*LOG(1.0_rp-rand1))*COS(2.0_rp*C_PI*rand2); END FUNCTION random_norm FUNCTION IntK(v,x) IMPLICIT NONE REAL(rp) :: IntK REAL(rp), INTENT(IN) :: v REAL(rp), INTENT(IN) :: x IntK = (C_PI/SQRT(2.0_rp))*(1.0_rp - 0.25_rp*(4.0_rp*v**2 - 1.0_rp))*ERFC(SQRT(x))& + 0.25_rp*(4.0_rp*v**2 - 1.0_rp)*SQRT(0.5_rp*C_PI/x)*EXP(-x) END FUNCTION IntK FUNCTION besselk(v,x) IMPLICIT NONE REAL(rp) :: besselk REAL(rp), INTENT(IN) :: x REAL(rp), INTENT(IN) :: v REAL(4) :: ri,rk,rip,rkp call bessik(REAL(x,4),REAL(v,4),ri,rk,rip,rkp) besselk = REAL(rk,rp) END FUNCTION besselk FUNCTION IntBesselK(a,b) REAL(rp), INTENT(IN) :: a REAL(rp), INTENT(IN) :: b REAL(rp) :: IntBesselK REAL(rp) :: Iold REAL(rp) :: Inew REAL(rp) :: rerr REAL(rp) :: sum_f REAL(rp) :: v,h,z INTEGER :: ii,jj,npoints LOGICAL :: flag v = 5.0_rp/3.0_rp h = b - a sum_f = 0.5*(besselk(v,a) + besselk(v,b)) Iold = 0.0_rp Inew = sum_f*h ii = 1_idef flag = .TRUE. do while (flag) Iold = Inew ii = ii + 1_idef npoints = 2_idef**(ii-2_idef) h = 0.5_rp*(b-a)/REAL(npoints,rp) sum_f = 0.0_rp do jj=1_idef,npoints z = a + h + 2.0_rp*(REAL(jj,rp) - 1.0_rp)*h sum_f = sum_f + besselk(v,z) end do Inew = 0.5_rp*Iold + sum_f*h rerr = ABS((Inew - Iold)/Iold) flag = .NOT.(rerr.LT.Tol) end do IntBesselK = Inew END FUNCTION IntBesselK SUBROUTINE P_integral(z,P) REAL(rp), INTENT(OUT) :: P REAL(rp), INTENT(IN) :: z REAL(rp) :: a P = 0.0_rp IF (z .LT. 0.5_rp) THEN a = (2.16_rp/2.0_rp**(2.0_rp/3.0_rp))*z**(1.0_rp/3.0_rp) P = IntBesselK(z,a) + IntK(5.0_rp/3.0_rp,a) ELSE IF ((z .GE. 0.5_rp).AND.(z .LT. 2.5_rp)) THEN a = 0.72_rp*(z + 1.0_rp) P = IntBesselK(z,a) + IntK(5.0_rp/3.0_rp,a) ELSE P = IntK(5.0_rp/3.0_rp,z) END IF END SUBROUTINE P_integral FUNCTION PR(eta,p,Bo,l) REAL(rp), INTENT(IN) :: eta ! in radians REAL(rp), INTENT(IN) :: p ! dimensionless (in units of mc) REAL(rp), INTENT(IN) :: Bo REAL(rp), INTENT(IN) :: l REAL(rp) :: PR REAL(rp) :: g REAL(rp) :: v REAL(rp) :: k REAL(rp) :: lc REAL(rp) :: z REAL(rp) :: Pi g = SQRT(p**2 + 1.0_rp) v = C_C*SQRT(1.0_rp - 1.0_rp/g**2) k = C_E*Bo*SIN(deg2rad(eta))/(g*C_ME*v) lc = (4.0_rp*C_PI/3.0_rp)/(k*g**3) ! Critical wavelength z = lc/l call P_integral(z,Pi) PR = (C_C*C_E**2)*Pi/(SQRT(3.0_rp)*C_E0*g**2*l**3) END FUNCTION PR SUBROUTINE sample_distribution(params,eta,etao) TYPE(KORC_PARAMS), INTENT(IN) :: params REAL(rp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: eta REAL(rp), INTENT(OUT) :: etao REAL(rp) :: go_root REAL(rp) :: etao_root REAL(rp) :: eta_buffer, eta_test REAL(rp) :: ratio, rand_unif REAL(rp), DIMENSION(:), ALLOCATABLE :: eta_samples REAL(rp), DIMENSION(:), ALLOCATABLE :: eta_tmp REAL(rp) :: minmax, min_pitch_angle, max_pitch_angle REAL(rp) :: deta LOGICAL :: leta INTEGER :: num_accepted INTEGER :: ii,jj,ppp,nsamples INTEGER :: mpierr ppp = SIZE(eta) nsamples = ppp*params%mpi_params%nmpi deta = (pdf_params%max_pitch_angle - pdf_params%min_pitch_angle)/100.0_rp if (pdf_params%min_pitch_angle.GE.korc_zero) then do jj=1_idef,INT(minmax_buffer_size,idef) minmax = pdf_params%min_pitch_angle - REAL(jj,rp)*deta if (minmax.GT.0.0_rp) then min_pitch_angle = minmax end if end do else min_pitch_angle = pdf_params%min_pitch_angle end if do jj=1_idef,INT(minmax_buffer_size,idef) minmax = pdf_params%max_pitch_angle + REAL(jj,rp)*deta if (minmax.LE.90.0_rp) then max_pitch_angle = minmax end if end do ! write(6,*) min_pitch_angle,max_pitch_angle if (params%mpi_params%rank.EQ.0_idef) then ALLOCATE(eta_samples(nsamples))! Number of samples to distribute among all MPI processes ALLOCATE(eta_tmp(nsamples))! Number of samples to distribute among all MPI processes !* * * Transient * * *! call RANDOM_SEED() call RANDOM_NUMBER(rand_unif) eta_buffer = pdf_params%min_pitch_angle + (pdf_params%max_pitch_angle - pdf_params%min_pitch_angle)*rand_unif ii=2_idef do while (ii .LE. 1000_idef) eta_test = eta_buffer + random_norm(0.0_rp,deta) do while ((ABS(eta_test) .GT. pdf_params%max_pitch_angle).OR.(ABS(eta_test) .LT. pdf_params%min_pitch_angle)) eta_test = eta_buffer + random_norm(0.0_rp,deta) end do ratio = fRE(eta_test,pdf_params%po)/fRE(eta_buffer,pdf_params%po) if (ratio .GE. 1.0_rp) then eta_buffer = eta_test ii = ii + 1_idef else call RANDOM_NUMBER(rand_unif) if (rand_unif .LT. ratio) then eta_buffer = eta_test ii = ii + 1_idef end if end if end do !* * * Transient * * *! eta_tmp(1) = eta_buffer call RANDOM_SEED() call RANDOM_NUMBER(rand_unif) num_accepted = 0_idef do while(num_accepted.LT.nsamples) ii=2_idef do while (ii .LE. nsamples) eta_test = eta_tmp(ii-1) + random_norm(0.0_rp,deta) do while ((ABS(eta_test) .GT. max_pitch_angle).OR.(ABS(eta_test) .LT. min_pitch_angle)) eta_test = eta_tmp(ii-1) + random_norm(0.0_rp,deta) end do ratio = fRE(eta_test,pdf_params%po)/fRE(eta_tmp(ii-1),pdf_params%po) if (ratio .GE. 1.0_rp) then eta_tmp(ii) = eta_test ii = ii + 1_idef else call RANDOM_NUMBER(rand_unif) if (rand_unif .LT. ratio) then eta_tmp(ii) = eta_test ii = ii + 1_idef end if end if end do eta_tmp = ABS(eta_tmp) ii = 1_idef do while ( (ii.LT.nsamples).AND.(num_accepted.LT.nsamples) ) leta = (eta_tmp(ii).LE.pdf_params%max_pitch_angle).AND.(eta_tmp(ii).GE.pdf_params%min_pitch_angle) if (leta) then num_accepted = num_accepted + 1_idef eta_samples(num_accepted) = eta_tmp(ii) end if ii = ii + 1_idef end do end do etao = SUM(eta_samples)/nsamples end if CALL MPI_SCATTER(eta_samples,ppp,MPI_REAL8,eta,ppp,MPI_REAL8,0,MPI_COMM_WORLD,mpierr) CALL MPI_BCAST(etao,1,MPI_REAL8,0,MPI_COMM_WORLD,mpierr) call MPI_BARRIER(MPI_COMM_WORLD,mpierr) if (params%mpi_params%rank.EQ.0_idef) then DEALLOCATE(eta_samples) DEALLOCATE(eta_tmp) end if END SUBROUTINE sample_distribution SUBROUTINE save_params(params) TYPE(KORC_PARAMS), INTENT(IN) :: params CHARACTER(MAX_STRING_LENGTH) :: filename CHARACTER(MAX_STRING_LENGTH) :: gname CHARACTER(MAX_STRING_LENGTH), DIMENSION(:), ALLOCATABLE :: attr_array CHARACTER(MAX_STRING_LENGTH) :: dset CHARACTER(MAX_STRING_LENGTH) :: attr INTEGER(HID_T) :: h5file_id INTEGER(HID_T) :: group_id INTEGER :: h5error REAL(rp) :: units if (params%mpi_params%rank .EQ. 0) then filename = TRIM(params%path_to_outputs) // "simple_equilibrium_pdf.h5" call h5fcreate_f(TRIM(filename), H5F_ACC_TRUNC_F, h5file_id, h5error) gname = "pdf_params" call h5gcreate_f(h5file_id, TRIM(gname), group_id, h5error) dset = TRIM(gname) // "/max_pitch_angle" attr = "Maximum pitch angle in avalanche PDF (degrees)" call save_to_hdf5(h5file_id,dset,pdf_params%max_pitch_angle,attr) dset = TRIM(gname) // "/min_pitch_angle" attr = "Minimum pitch angle in avalanche PDF (degrees)" call save_to_hdf5(h5file_id,dset,pdf_params%min_pitch_angle,attr) dset = TRIM(gname) // "/Zeff" attr = "Effective atomic number of ions." call save_to_hdf5(h5file_id,dset,pdf_params%Zeff,attr) dset = TRIM(gname) // "/E" attr = "Parallel electric field in (Ec)" call save_to_hdf5(h5file_id,dset,pdf_params%E,attr) dset = TRIM(gname) // "/Bo" attr = "Characteristic magnetic field in T (in case of using PR)" call save_to_hdf5(h5file_id,dset,pdf_params%Bo,attr) dset = TRIM(gname) // "/lambda" attr = "Characteristic wavelength in m (in case of using PR)" call save_to_hdf5(h5file_id,dset,pdf_params%lambda,attr) call h5gclose_f(group_id, h5error) call h5fclose_f(h5file_id, h5error) end if END SUBROUTINE save_params END MODULE korc_simple_equilibrium_pdf