Code 1: M2D - Monte Carlo, Square Geometry c program m2d c-----===========------------------------------------------------------- c c monte carlo program, square geometry for gs580 c c originally from c mark christon c gs-580 c supercomputing applications c spring semester, 1987 c c modified by: pat burns c c 10/27/88 c c----------------------------------------------------------------------- c2 4 6 8(1)2 4 6 8(2)2 4 6 8(3)2 4 6 8(4)2 4 6 8(5)2 4 6 8(6)2 4 6 8(7)2 c parameter ( np = 10000, npsurf = 100000, ns = 4) c c character*8 adate,date c205 c c define common blocks c common /surf/ x(ns,2),y(ns,2),delx(ns),dely(ns),s12(ns), 2 cross(ns),enx(ns),eny(ns),x0(ns),y0(ns),e(ns,ns), 3 delxp(ns),delyp(ns),errsumi(ns),err(ns,ns), 4 etrue(ns,ns),ie(ns,ns) common /mat/ thetar(101),delthr(101),pi,twopi,pi2 common /photon/ xe(np),ye(np),ex(np),ey(np) c c initialization c t0=second() adate=date() c c input calculations c ti1 = second() call input ti2 = second() c nstop=npsurf/np if (np*nstop.ne.npsurf) then write (6,3000) write (6,3001) np,npsurf,nstop stop endif c c loop over surfaces c ts1 = second() do 100 le = 1,ns c c loop over photon emissions c do 100 n = 1,nstop c c...........emit 'np' photons c call emit (le,n) c c...........intersection loop c call intsec (le) 100 continue ts2 = second() c c.....print results, evaluate error c rnp = float(npsurf) do 150 i=1, ns do 150 j=1, ns e(i,j) = float( ie(i,j) )/rnp 150 continue c do 200 i=1, ns do 200 j=1, ns err(i,j) = e(i,j)-etrue(i,j) 200 continue c do 300 i=1, ns errsumi(i) = 0.0 do 300 j=1, ns errsumi(i) = errsumi(i)+abs(err(i,j)) 300 continue c errsum = 0. do 400 i=1,ns 400 errsum = errsum + errsumi(i) c tend = second() ti = ti2 - ti1 tmc = ts2 - ts1 ttot = tend - t0 write(6,5000) adate do 490 i=1,ns 490 write(6,5010) i,x(i,1), y(i,1), x(i,2), y(i,2) write(6,5020) np, npsurf, nstop, ns c write(6,5025) do 500 i=1,ns write(6,5026) i 500 write(6,5030) (e(i,j), j=1, ns) c write(6,5035) do 510 i=1,ns write(6,5026) i 510 write(6,5040) (err(i,j), j=1, ns) c write(6,5041) do 520 i=1,ns 520 write(6,5042) i,errsumi(i) c write(6,5045) errsum c write(6,5050) ti, tmc, ttot c if (tmc.le.0.) tmc = 1. rnflts = npsurf*ns fltper = rnflts / tmc write (6,5060) tmc,rnflts,fltper c stop c 3000 format (/ 2 /1x,'--- error termination ---'/) 3001 format (/ 2 / 1x,'bad number for photons, np, npsurf, nstop ='/5x,3i10/ 3 / 10x, '--- stop ---'/) 5000 format('1',//,5x, 2 '**************************************************************' 3 ,//,10x,'monte 2-d polygonal geometry',//,10x,a8,//, 4 10x,'geometry definition',//,7x,'no.',6x,'x1',12x,'y1',12x, 5 'x2',12x,'y2',//) 5010 format(5x,i5,'.',4(2x,e12.5)) 5020 format(//,10x,'*** parameter data ***',//,10x,'np= ',i5,/, 2 10x,'npsurf= ',i6,/,10x,'nstop= ',i5,/,10x,'ns= ',i5/) 5025 format(//,15x,'exchange factor matrix follows, row-by-row'//) 5026 format(5x,'****** row =',i5) 5030 format(10(1x,f11.7)) 5035 format(//,15x,'error matrix follows, row-by-row',//) 5040 format(10(1x,e11.4)) 5041 format(//,15x,'row sums of error follow, row-by-row'/) 5042 format(5x,i6,'.',5x,e15.8) 5045 format(//10x,'total error, summed abs = ',e15.8//) 5050 format(//,10x,'*** timing summary ***',/, 2 10x,'input time= ',e10.4,' [sec]',/, 3 10x,'calc. time= ',e10.4,' [sec]',/, 4 10x,'total time= ',e10.4,' [sec]') 5060 format(//5x,'totals follow'/ * / 5x,'solution phase time =',e12.5,' [sec]', * / 5x,'no. of flights =',e12.5 * / 5x,'flight rate =',e12.5,' [1/sec]') c end subroutine input c-----================-------------------------------------------------- c c monte carlo program - input routine for square geometry c c----------------------------------------------------------------------- c parameter ( np = 10000, npsurf = 100000, ns = 4) c common /surf/ x(ns,2),y(ns,2),delx(ns),dely(ns),s12(ns), 2 cross(ns),enx(ns),eny(ns),x0(ns),y0(ns),e(ns,ns), 3 delxp(ns),delyp(ns),errsumi(ns),err(ns,ns), 4 etrue(ns,ns),ie(ns,ns) common /mat/ thetar(101),delthr(101),pi,twopi,pi2 c enorm(x1,x2,y1,y2) = sqrt((x1-x2)**2+(y1-y2)**2) c c deifne constants c pi = 3.1415926535 twopi = 2.*pi pi2 = 0.5*pi c c define surfaces clockwise from x-axis c alpha = twopi/ns c x(1,1) = 1. x(1,2) = 1. y(1,1) = tan(alpha*0.5) y(1,2) = -y(1,1) c radius = enorm(x(1,1),0.,y(1,1),0.) c beta = -0.5*alpha do 10 i=2,ns beta = beta - alpha x(i,1) = x(i-1,2) y(i,1) = y(i-1,2) x(i,2) = radius*cos(beta) 10 y(i,2) = radius*sin(beta) c c last surface - to eliminate 'leaks' c x(ns,2) = x(1,1) y(ns,2) = y(1,1) c do 100 i = 1,ns delx(i) = x(i,2)-x(i,1) dely(i) = y(i,2)-y(i,1) s12(i) = (delx(i))**2+(dely(i))**2 cross(i)= x(i,1)*y(i,2)-x(i,2)*y(i,1) el = sqrt(s12(i)) enx(i) = dely(i)/el eny(i) = -delx(i)/el delxp(i) = delx(i)/(npsurf) delyp(i) = dely(i)/(npsurf) x0(i) = x(i,1)+0.5*delxp(i) y0(i) = y(i,1)+0.5*delyp(i) 100 continue c c set material property arrays c do 200 i= 1,101 r = float(i-1)/100. thetar(i) = asin(sqrt(r)) 200 continue c do 300 i = 1,100 delthr(i)= thetar(i+1)-thetar(i) 300 continue c delthr(101) = 0. c c initialize arrays c do 400 i = 1,ns do 400 j = 1,ns 400 ie(i,j) = 0.0 c c initialize true answers - Hottel's crossed-string method c el = enorm (x(1,1),x(1,2),y(1,1),y(1,2)) elinv = 0.5/el c do 500 i=1,ns 500 etrue(i,i) = 0.0 c do 550 i=2,ns cross1 = enorm (x(1,1),x(i,1),y(1,1),y(i,1)) cross2 = enorm (x(1,2),x(i,2),y(1,2),y(i,2)) strch1 = enorm (x(1,1),x(i,2),y(1,1),y(i,2)) strch2 = enorm (x(1,2),x(i,1),y(1,2),y(i,1)) etrue(1,i) = elinv*(cross1+cross2-strch1-strch2) 550 continue c do 580 i=2,ns do 580 j=1,ns if (j.ne.1) then etrue(i,j) = etrue(i-1,j-1) else etrue(i,j) = etrue(i-1,ns) endif 580 continue c c initialize random number generator c c call ranset( 0 ) cy205 c call setrvn( 65533 ) cy205 c return c end subroutine emit (l,n) c-----=====================--------------------------------------------- c c monte carlo emission routine c c----------------------------------------------------------------------- c parameter ( np = 10000, npsurf = 100000, ns = 4) c common /surf/ x(ns,2),y(ns,2),delx(ns),dely(ns),s12(ns), 2 cross(ns),enx(ns),eny(ns),x0(ns),y0(ns),e(ns,ns), 3 delxp(ns),delyp(ns),errsumi(ns),err(ns,ns), 4 etrue(ns,ns),ie(ns,ns) common /mat/ thetar(101),delthr(101),pi,twopi,pi2 common /photon/ xe(np),ye(np),ex(np),ey(np) c dimension ran(np),iran(np),frac(np),theta(np),phi(np),exloc(np), 2 eyloc(np),scratch(np),iran1(np) c c define starting coordinates c rnm1 = float(n-1) rnp = float(np) xstrt = x0(l) + rnm1*delxp(l)*rnp ystrt = y0(l) + rnm1*delyp(l)*rnp c c get 'np' emission coordinates - xe,ye c dx = delxp(l) dy = delyp(l) do 100 i=1,np xe(i) = ... 100 ye(i) = ... c c get 'np' local spherical angles - theta, phi c do 200 i=1,np 200 ran(i) = ranf() do 210 i=1,np 210 theta(i) = ... c do 300 i=1,np 300 ran(i) = ranf() do 310 i=1,np 310 phi(i) = ... c c get 'np' local cartesian angles - exloc, eyloc c c do 400 i=1,np exloc(i) = ... 400 eyloc(i) = ... c c get 'np' global cartesian angles - ex, ey c enyl = eny(l) enxl = enx(l) do 500 i=1,np ex(i) = ... 500 ey(i) = ... c c write(6,1000) c write(6,1010) (xe(i),ye(i),theta(i),phi(i),exloc(i),eyloc(i), c 2 ex(i),ey(i), i=1, np) c return c c1000 format(2x,'xe',5x,'ye',5x,'theta',5x,'phi',5x,'exloc',5x, c 2 'eyloc',5x,'ex',5x,'ey',//) c1010 format(8(1x,e12.6)) c1020 format(5x,'ran(i)= ',e12.6) end subroutine intsec (le) c-----======================-------------------------------------------- c c monte carlo intersection routine c c----------------------------------------------------------------------- c parameter ( np = 10000, npsurf = 100000, ns = 4) c common /surf/ x(ns,2),y(ns,2),delx(ns),dely(ns),s12(ns), 2 cross(ns),enx(ns),eny(ns),x0(ns),y0(ns),e(ns,ns), 3 delxp(ns),delyp(ns),errsumi(ns),err(ns,ns), 4 etrue(ns,ns),ie(ns,ns) common /photon/ xe(np),ye(np),ex(np),ey(np) c dimension crossp(np),crosspc(np),exc(np),eyc(np),detc(np), 2 detinvc(np),xic(np),yic(np),scratch(np),dotp(np), 3 dotpc(np),s1c(np),s2c(np),stestc(np),indexc(np) c data deteps/ 1.0e-12 / c c define photon quantities constant over surface loop c do 10 i=1,np 10 crossp(i) = ... c c loop over surfaces to determine intersection points c do 600 l= 1,ns c if( l .eq. le ) go to 600 c c........constants within loop - store in registers c delxl = delx(l) delyl = dely(l) crossl = cross(l) x1l = x(l,1) x2l = x(l,2) y1l = y(l,1) y2l = y(l,2) s12l = s12(l) enxl = enx(l) enyl = eny(l) c loop over photons c do 100 i=1,np c c........evaluate dot product between photon and surface normal c 100 dotp(i) = ... c c c2 4 6 8(1)2 4 6 8(2)2 4 6 8(3)2 4 6 8(4)2 4 6 8(5)2 4 6 8(6)2 4 6 8(7)2 c c........obtain vector (array) of indices where dotp is negative c call whenflt (np,dotp,1,0.,indexc,lenc1) cray2 c c write(6,1000) lenc1 c if( lenc1 .ge. 1 ) then c c........compress the photons incident from the front c do 200 i=1,lenc1 exc(i) = ... eyc(i) = ... 200 crosspc(i)= ... c c........calculate intersection points c do 300 i=1,lenc1 detc(i) = ... detinvc(i) = ... xic(i) = ... yic(i) = ... c write(6,1020) (xic(i), yic(i), i=1,lenc1) 300 continue c c........test for intersection between endpoints of surfaces c do 400 i=1,lenc1 s1c(i) = ... s2c(i) = ... 400 stestc(i) = ... c c........increment e(le,l) by photons absorbed c c npa = 0 c2 4 6 8(1)2 4 6 8(2)2 4 6 8(3)2 4 6 8(4)2 4 6 8(5)2 4 6 8(6)2 4 6 8(7)2 c call whenflt (lenc1,stestc,1,0.,indexc,npa) cray2 c write(6,1010) l,le,npa ie(le,l) = ie(le,l) + npa c endif c 600 continue c return c 1000 format(5x,'lenc1= ',i8) 1010 format(5x,'l= ',i5,5x,'le= ',i5,5x,'npa= ',i8) 1020 format(5x,'xic(i)= ',e12.6,5x,'yic(i)= ',e12.6) 1030 format(///,5x,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!',/, 2 5x,'!!!!! error in determinant !!!!!',/, 3 5x,'!!!!! check geometry !!!!!',/, 4 5x,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') end