! program tau-dis.f90 ! IRI = Irradiation Instability in a Keplerian, optically thick disk ! ! includes dislin graphics, ! compiles with ifort (all options not necessary..) ! ifort -qopenmp -O3 -no-prec-div -fp-model fast=2 -xHost \ ! -L$DISLIN10 -ldislin10 -I$DISLIN10/ifc \ ! -L/opt/intel/composer_xe_2015.3.187/compiler/lib/intel64 -lirc !*.f90 -o !*.x ! ! Experiment with changing: ! pgf90 compiler ! schedule(static,Np/Nthreads) -> schedule(dynamic,Np/Nthreads) or none ! Nthreads = 1, 6, 24 ! $omp atomic, $omp critical updates to Map -> watch the speed drop ! ! dtau = density map of particles. Size: Nr by Nphi. ! tau = optical thickness of the disk (map) converted to beta(r,phi) ! partic = table of particle positions, L, and radial speeds: ! 1=r, 2=phi, 3=beta, 4=v_phi, 5= L, 6 = v_r ! ! Disk outer radius ~ Rout ! Np = number of particles module dataset use omp_lib implicit none ! constants integer, parameter :: Tau_inf = 5, Beta = 0.2 integer, parameter :: Nthreads = 12, Nj = 8, Mil = 1000000 integer, parameter :: Nr = 400/12*13, Nphi = Nr integer, parameter :: Np = 1*Mil /Nthreads*(1+Nthreads) real, parameter :: pi = atan(1.)*4., Rout = 2.5 real, parameter :: T_max = pi real, parameter :: dr = Rout/Nr ! arrays and variables real :: partic(Nj,Np) ! particle data: r,phi,vr,L real :: dtau(0:Nr-1,0:Nphi-1), tau(0:Nr-1,0:Nphi-1) ! dtau,tau (r,phi) real :: rand_r(2,Np*3) real :: rad, c = 1.0 real :: time = 0., dt = 0.004 integer :: ix, iy, mode, i, j contains !___________________________________________________________________ real function Sigma(x) real :: arg, x arg = max( min(x-0.666,0.666),0.) ! 0..1 for x = 2/3...4/3 Sigma = sin(arg*pi/2.) ! sine ramp if (x > Rout*0.95) Sigma = 0. return end !___________________________________________________________________ subroutine rad_pressure (mode) integer :: ir, iphi, mode ! compute the optical depth from particles positions ! contribute to the map of dtau, then integrate dtau = 0. if(mode==0) c = 1. !$omp parallel do num_threads(Nthreads) private(ir,iphi) do i = 1, Np ir = partic(1,i)/Rout*Nr ir = max(0,min(ir, Nr-1)) iphi = max(0,min(int(partic(2,i)/2./pi*Nphi), Nphi-1)) ! Histogram update can cause race condition. We do diagnose the ! magnitude of the problem by comparing the total of Map entries ! with the number of particles thrown in it. If P(race condition) ! is censiderable, then we can activate either the critical region ! (guaranteed to be done by single thread) or the atomic update ! (guaranteed non-overlapping writes to memory by threads). ! In our program, race occurs with P < 0.0003, and we do not activate ! those OpenMP safeguards, since they slow down the execution of the ! program 3x (atomic summation) and 30 times (critical region). !///$omp critical !///$omp atomic dtau(ir,iphi) = dtau(ir,iphi) + c ! partic(1,i)**(-2) !///$omp end critical end do ! summarize accurate the mapping is if(mode==0) print*,'mapping accuracy',sum(dtau)/Np ! integrate along radius tau(0,:) = 0. !$omp parallel do num_threads(Nthreads) schedule(static,Nphi/12) & !$omp private(i) do j = 0, Nphi-1 do i = 1, Nr-1 tau(i,j) = tau(i-1,j) + dr*dtau(i,j)/(i*i) ! dtau ~ 1/r^2 end do end do ! get proper Tau normalization in mode=0 from tau(outer edge) if (mode==0) then c = Tau_inf * Nphi/sum(tau(Nr-1,:)) print*,' mode 0, c=',c dtau = dtau * c tau = tau * c end if ! convert tau to beta tau = exp(-tau) * Beta ! tau now holds beta(r,phi) ! now pull the resulting betas out of the array !$omp parallel do num_threads(Nthreads) private(ir,iphi) do i = 1, Np ir = partic(1,i)/Rout*Nr -1 ir = max(0,min(ir,Nr-1)) iphi = max(0,min(int(partic(2,i)/2./pi*Nphi), Nphi-1)) partic(3,i) = tau(ir,iphi) ! beta(t) of particle i end do end subroutine rad_pressure !_____________________________________________________________________ subroutine initialize() integer :: n real :: force(Np), ran ! randomize arrays call random_number(rand_r) call random_number(partic) ! initialize r, phi ! try up to 3*Np different r's and reject some ! in order to produce radial distribution with density Sigma(r) n = 0 !$omp critical do i = 1, 3*Np ran = rand_r(1,i) ! random 0..1 rad = rand_r(2,i) * Rout ! trial radius, 0..Rout if (ran < Sigma(rad)) then ! accept particle's r n = n + 1; if (n > Np) exit ! when all r,phi chosen partic(2,n) = pi*2.*partic(2,n) ! phi = 0...pi*2 partic(1,n) = rad ! radius if(n<100) then print*,n,' r,p',rad,partic(2,n) end if end if end do !$omp end critical ! compute tau(r) and radiation pressure force print*,'r,phi established, calling rad_press.' call rad_pressure(0) ! finish initialization of velocity, based on beta partic(7,:) = (1.-partic(3,:))/partic(1,:)**2 ! GM(1-beta)/r^2 partic(4,:) = sqrt(partic(7,:)*partic(1,:)) ! v_phi from smooth force partic(4,:) = partic(4,:) *(1.+0.04*(partic(5,:)-0.5)) ! randomize partic(5,:) = partic(4,:)*partic(1,:) ! store vphi * r = L = const. partic(6,:) = 0.02 * partic(1,:)**(-0.5) ! vr = +-0.01 * vK end subroutine initialize !________________________________________________________________________ subroutine timestep(step) integer :: step time = time + dt ! evaluate optical thickness, beta = F_rad/F_grav, ! and radial forces (accelerations) call rad_pressure (1) ! evolve particles using angular momentum conservation partic(4,:) = partic(5,:)/partic(1,:) ! L/r = v_phi partic(2,:) = mod(partic(2,:)+dt*partic(4,:), 2*pi) ! phi+=v_phi*dt -1 partic(6,:) = partic(6,:) + dt*(-partic(7,:) + & ! r += (dr/dt)*dt +partic(4,:)**2/partic(1,:)) partic(1,:) = partic(1,:) + dt*partic(6,:) ! r += (dr/dt)*dt end subroutine timestep end module dataset !_______________________________________________________________________ ! MAIN DRIVER ! program IRI use dataset real (kind=8) :: t integer :: step ! init call initialize() print*,'disk initialized' ! main time loop do step = 0, int(T_max/dt) t = omp_get_wtime() call timestep(step) t = omp_get_wtime() -t if(mod(step,40)==0) print*,'time',time,' t[s]',real(t) end do ! sometimes do i = 1,Np,Np/40 print('(7e11.3)'),partic(1:7,i) end do print*,'max tau()',maxval(tau) print*, tau(398:403, 399) print*, tau(398:403, 400) print*, tau(398:403, 401) print*, ' Nr-1',tau(nr_1, 0:Nphi-1:10) ! display particle map call display() end program IRI !______________________________________________________________________ !______________________________________________________________________ ! GRAPHICS subroutine display() use dataset use dislin real :: arr(Nr,Nphi) real :: a0, a1, v0, v1 ! prepare normalized data v1 = maxval(dtau) print*,' min, max dtau',minval(dtau),v1 arr = alog10(0.01+dtau) ! arr = arr/(1.+(arr/2)**4)**0.25 a0 = minval(arr) a1 = maxval(arr) print*,' min/max val arr =',a0,a1 print*,' array slice',arr(396:800:4,400) ! plot !call scrmod('REVERS') call metafl('CONS') call disini() call pagera() call hwfont() call titlin('Map of particles in polar coordinates',4) call titlin('Disk with total optical depth = 5',2) call name('R ','x') call name('Phi ','y') call name(' ~ log (surface density)','z') call intax() call autres(Nr,Nphi) call axspos(300,1850) ! (300,1850) call ax3len(1400,1400,1400) ! (2200,1400,1400) a0 = 1. call graf3(0.,2.5,0.,0.2,0.,2*pi,0.,pi/4,a0,a1,a0,(a1-a0)/10) call crvmat(arr,Nr,Nphi,1,1) call height(50) call title() !call mpaepl(3) call disfin() stop end subroutine display !_________________________________________________________________________