! !----------------------------------------------------------------------- ! !Purpose: ! Generate a dataset for a Monte Carlo code. ! !Copyright: ! Copyright (C) 2001 Ecole Polytechnique de Montreal. ! This library is free software; you can redistribute it and/or ! modify it under the terms of the GNU Lesser General Public ! License as published by the Free Software Foundation; either ! version 2.1 of the License, or (at your option) any later version. ! !Author(s): ! G. Civario (CS-SI) ! !Comments: ! fonctions: ! - initMC : initialisation des tableaux ! - destroyMC : liberation des pointeurs ! - createElemGene : constructeur d'un element generique ! - isSameEG : teste l'equivalence de deux elements generiques ! - isSameWayEG : teste si deux elements generiques ont meme orientation ! - prepareMCData : preparation des donnes necessaires et remplissage ! des tableaux globaux ! - generateTripoliFile : ecriture du fichier de donnees Tripoli4 ! - generateMCNPFile : ecriture du fichier de donnees MCNP ! - generateSerpentFile : ecriture du fichier de donnees Serpent ! - putOn80col : formattage sur 80 colonnes style MCNP, d'un buffer a afficher ! - findParalleleWithTrans : trouve l'element geometrique a associer avec un ! autre, dans le cas d'une translation ! !----------------------------------------------------------------------- ! module monteCarlo use boundCond use cast use cellulePlaced use constType use constUtiles use segArc use generSAL implicit none type t_elemGene logical :: isPlan !plan -> true ; cylindre -> false double precision :: x,y ! -> origine ; -> centre double precision :: dx,dy ! -> extremite ; -> 0. double precision :: r ! -> 0. ; -> rayon integer :: limite ! ==Tri_XXXX end type t_elemGene type(t_elemGene),dimension(:),allocatable,save :: tabEG type t_volume integer :: lg !longueur des tableaux integer :: mix !numero du milieu integer,dimension(:), allocatable :: indElem !indice des elements generiques !du contour logical,dimension(:), allocatable :: side !true -> cote + ; false -> cote - integer,dimension(:), allocatable :: typCL !condition limite de la face end type t_volume type(t_volume),dimension(:),allocatable,save :: tabVolume integer,parameter :: Tri_Not=-1,Tri_Void=0,Tri_Refl=1,Tri_Trans=2,Tri_Cos=3 contains subroutine initMC(nbNode) integer,intent(in) :: nbNode allocate(tabVolume(nbNode),stat=alloc_ok) if (alloc_ok /= 0) call XABORT("G2S: initMC => allocation pb") tabVolume(1:nbNode)%lg = 0 end subroutine initMC subroutine destroyMC() integer :: i do i = 1,size(tabVolume) deallocate(tabVolume(i)%indElem) deallocate(tabVolume(i)%side) deallocate(tabVolume(i)%typCL) end do deallocate(tabEG,tabVolume) end subroutine destroyMC function createElemGene(sa,strDCL) type(t_segArc),intent(in) :: sa character(len=4),intent(in) :: strDCL type(t_elemGene) :: createElemGene integer :: sunsetCL createElemGene%isPlan = (sa%typ==tseg) createElemGene%x = sa%x createElemGene%y = sa%y if (createElemGene%isPlan) then createElemGene%dx = sa%dx createElemGene%dy = sa%dy createElemGene%r = 0.d0 else createElemGene%dx = 0.d0 createElemGene%dy = 0.d0 createElemGene%r = sa%r end if createElemGene%limite = Tri_Not sunsetCL = minval((/sa%nodeg,sa%noded/)) if (sunsetCL<=0) then sunsetCL = mod(-sunsetCL,100) select case(sunsetCL) case(B_Void,B_Zero) createElemGene%limite = Tri_Void case(B_Refl,B_Ssym,B_Syme,B_Diag) createElemGene%limite = Tri_Refl case(B_Tran) createElemGene%limite = Tri_Trans case(B_Albe) createElemGene%limite = Tri_Cos case(-fooMix,0) !il faut utiliser la CL par defaut if (strDCL=='ALBE') then createElemGene%limite = Tri_Cos else if (strDCL=='REFL') then createElemGene%limite = Tri_Refl else ! => VOID createElemGene%limite = Tri_Void end if case default call XABORT("G2MC: boundary condition not allowed") end select end if end function createElemGene function isSameEG(eg1,eg2) type(t_elemGene),intent(in) :: eg1,eg2 logical :: isSameEG if (eg1%isPlan.neqv.eg2%isPlan) then isSameEG = .false. else if (eg1%isPlan) then isSameEG = pointsAlignes(eg1%x,eg1%y,eg1%dx,eg1%dy,eg2%x,eg2%y) & .and. pointsAlignes(eg1%x,eg1%y,eg1%dx,eg1%dy,eg2%dx,eg2%dy) else isSameEG = isEqualConst(eg1%x,eg2%x) & .and. isEqualConst(eg1%y,eg2%y) & .and. isEqualConst(eg1%r,eg2%r) end if isSameEG = isSameEG .and. (eg1%limite==eg2%limite) end function isSameEG function isSameWayEG(eg1,eg2) type(t_elemGene),intent(in) :: eg1,eg2 logical :: isSameWayEG if (eg1%isPlan .and. eg2%isPlan) then isSameWayEG = & (((eg1%dx-eg1%x)*(eg2%dx-eg2%x)+(eg1%dy-eg1%y)*(eg2%dy-eg2%y))>0) else isSameWayEG = .true. end if end function isSameWayEG subroutine prepareMCData(szSA,nbNode,szEG) integer,intent(in) :: szSA,nbNode integer,intent(inout) :: szEG integer :: i,j,volNbr,defautCl,lgMax,indFictive logical :: found,toPut real :: albedo character*4 :: strDCL double precision :: ptx,pty,rxx,ryy type(t_elemGene) :: eg integer,dimension(:),allocatable :: indEG,nbPoints logical,dimension(:),allocatable :: isGoodWay,withFictiveEG,inACircle double precision,dimension(:,:),allocatable :: xx,yy type(t_elemGene),dimension(:),allocatable :: tmpTabEG !recuperation des donnees de conditions aux limites call calculDefaultCl(defautCl,albedo,strDCL) !preparation des donnees de remplissage des tableaux globaux allocate(indEG(szSA),isGoodWay(szSA),tmpTabEG(szSA+nbNode),stat=alloc_ok) if (alloc_ok /= 0) call XABORT("G2S: prepareMCData(1) => allocation pb") do i = 1,szSA eg = createElemGene(tabSegArc(i),strDCL) found = .false. do j = 1,szEG found = isSameEG(eg,tmpTabEG(j)) if (found) exit end do if (found) then indEG(i) = j isGoodWay(i) = isSameWayEG(eg,tmpTabEG(j)) else szEG = szEG + 1 tmpTabEG(szEG) = eg indEG(i) = szEG isGoodWay(i) = .true. end if end do !preparation des donnes permettant de savoir si un element fictif !englobant est necessaire allocate(inACircle(nbNode),withFictiveEG(nbNode),stat=alloc_ok) if (alloc_ok /= 0) call XABORT("G2S: prepareMCData(2) => allocation pb") withFictiveEG(:nbNode) = .false. ; inACircle(:nbNode) = .false. !calcul des dimensions des pointeurs du tableau des volumes do i = 1,szSA !cote gauche volNbr = tabSegArc(i)%nodeg if (volNbr>0) then tabVolume(volNbr)%lg = tabVolume(volNbr)%lg + 1 if (tabSegArc(i)%typ==tcer .and. .not. inACircle(volNbr)) then !le node est a l'interieur d'un cercle => on n'a pas besoin de !l'englober inACircle(volNbr) = .true. if (withFictiveEG(volNbr)) then !on avait commencer a englober fictivement => on laisse tomber tabVolume(volNbr)%lg = tabVolume(volNbr)%lg - 1 withFictiveEG(volNbr) = .false. end if end if end if !cote droit volNbr = tabSegArc(i)%noded if (volNbr>0) then tabVolume(volNbr)%lg = tabVolume(volNbr)%lg + 1 if (tabSegArc(i)%typ==tarc .and. .not. withFictiveEG(volNbr) & .and. .not. inACircle(volNbr)) then !le node est a l'exterieur d'un arc => on englobe si pas encore !fait tabVolume(volNbr)%lg = tabVolume(volNbr)%lg + 1 withFictiveEG(volNbr) = .true. end if end if end do deallocate(inACircle) lgMax = maxval(tabVolume(:nbNode)%lg) !allocation des pointeurs du tableau des volumes do i = 1,nbNode allocate(tabVolume(i)%indElem(tabVolume(i)%lg),stat=alloc_ok) if (alloc_ok /= 0) call XABORT("G2S: prepareMCData(3) => allocation pb") tabVolume(i)%indElem = 0 allocate(tabVolume(i)%side(tabVolume(i)%lg),stat=alloc_ok) if (alloc_ok /= 0) call XABORT("G2S: prepareMCData(4) => allocation pb") allocate(tabVolume(i)%typCL(tabVolume(i)%lg),stat=alloc_ok) if (alloc_ok /= 0) call XABORT("G2S: prepareMCData(5) => allocation pb") tabVolume(i)%typCL = 0 !remise a zero du compteur, pour gerer les doublons tabVolume(i)%lg = 0 end do !remplissage des pointeurs, avec elimination des doublons allocate(xx(nbNode,lgMax),yy(nbNode,lgMax)) allocate(nbPoints(nbNode),stat=alloc_ok) if (alloc_ok /= 0) call XABORT("G2S: prepareMCData(6) => allocation pb") nbPoints(:nbNode) = 0 ; xx(:nbNode,:lgMax) = 0.d0 ; yy(:nbNode,:lgMax) = 0.d0 do i = 1,szSA !cote gauche volNbr = tabSegArc(i)%nodeg if (volNbr>0) then toPut = ( count(tabVolume(volNbr)%indElem(:)==indEG(i)) == 0 ) if (toPut) then tabVolume(volNbr)%mix = tabSegArc(i)%neutronicMixg tabVolume(volNbr)%lg = tabVolume(volNbr)%lg + 1 tabVolume(volNbr)%indElem(tabVolume(volNbr)%lg) = indEG(i) tabVolume(volNbr)%typCL(tabVolume(volNbr)%lg) = & tmpTabEG(indEG(i))%limite tabVolume(volNbr)%side(tabVolume(volNbr)%lg) = .not.isGoodWay(i) end if !preparation du calcul des coordonnees du cercle englobant if (withFictiveEG(volNbr).and.tabSegArc(i)%typ/=tcer) then call giveOrigine(tabSegArc(i),ptx,pty) nbPoints(volNbr) = nbPoints(volNbr) + 1 xx(volNbr,nbPoints(volNbr)) = ptx yy(volNbr,nbPoints(volNbr)) = pty end if end if !cote droit volNbr = tabSegArc(i)%noded if (volNbr>0) then toPut = ( count(tabVolume(volNbr)%indElem(:)==indEG(i)) == 0 ) if (toPut) then tabVolume(volNbr)%mix = tabSegArc(i)%neutronicMixd tabVolume(volNbr)%lg = tabVolume(volNbr)%lg + 1 tabVolume(volNbr)%indElem(tabVolume(volNbr)%lg) = indEG(i) tabVolume(volNbr)%typCL(tabVolume(volNbr)%lg) = & tmpTabEG(indEG(i))%limite tabVolume(volNbr)%side(tabVolume(volNbr)%lg) = isGoodWay(i) end if !preparation du calcul des coordonnees du cercle englobant if (withFictiveEG(volNbr).and.tabSegArc(i)%typ/=tcer) then call giveExtremite(tabSegArc(i),ptx,pty) nbPoints(volNbr) = nbPoints(volNbr) + 1 xx(volNbr,nbPoints(volNbr)) = ptx yy(volNbr,nbPoints(volNbr)) = pty end if end if end do !calcul des cercles englobants do i = 1,nbNode if (.not.withFictiveEG(i)) cycle !creation du cercle eg%isPlan = .false. eg%x = sum(xx(i,:nbPoints(i)))/nbPoints(i) eg%y = sum(yy(i,:nbPoints(i)))/nbPoints(i) eg%dx = 0.d0 ; eg%dy = 0.d0 ; eg%r = 0.d0 ; eg%limite = Tri_Not do j = 1,nbPoints(i) rxx = xx(i,j)-eg%x ; ryy = yy(i,j)-eg%y eg%r = max(eg%r,longVect(rxx,ryy)) end do !integration eventuelle dans le tableau des elements geometriques found = .false. do j = 1,szEG found = isSameEG(eg,tmpTabEG(j)) if (found) exit end do if (found) then indFictive = j else szEG = szEG + 1 tmpTabEG(szEG) = eg indFictive = szEG end if !ajout eventuel dans la liste des elements geometriques constitutifs !de la zone toPut = ( count(tabVolume(i)%indElem(:)==indFictive) == 0 ) if (toPut) then tabVolume(i)%lg = tabVolume(i)%lg + 1 tabVolume(i)%indElem(tabVolume(i)%lg) = indFictive tabVolume(i)%typCL(tabVolume(i)%lg) = Tri_Not tabVolume(i)%side(tabVolume(i)%lg) = .false. end if end do deallocate(nbPoints,xx,yy,indEG,isGoodWay,withFictiveEG) !recopie du tableau global des elements geometriques, et liberation !du temporaire allocate(tabEG(szEG),stat=alloc_ok) if (alloc_ok /= 0) call XABORT("G2S: prepareMCData(7) => allocation pb") tabEG(1:szEG) = tmpTabEG(1:szEG) deallocate(tmpTabEG) end subroutine prepareMCData subroutine generateTripoliFile(fileTripoli,szSA,nbNode) integer,intent(in) :: fileTripoli,szSA,nbNode real :: a,b,c,d,xx,yy,rr integer :: i,j,szEG,nbPlus,nbMoins,nbMix,nbCLTot,nbCL integer,dimension(:),allocatable :: sp,sm,listMix,volmix,clList szEG = 0 call initMC(nbNode) call prepareMCData(szSA,nbNode,szEG) write(fileTripoli,'(/"GEOMETRY"/)') write(fileTripoli,'("//jeu de donnees geometriques TRIPOLI")') write(fileTripoli,'("//genere pour DRAGON par le module G2MC"/)') write(fileTripoli,'("TITRE geometrie DRAGON pour TRIPOLI"/)') !definition des surfaces write(fileTripoli,'("//definition des surfaces"/)') c = 0. do i = 1,szEG if (tabEG(i)%isPlan) then a = real(tabEG(i)%dy - tabEG(i)%y) b = real(tabEG(i)%x - tabEG(i)%dx) d = real(tabEG(i)%y*tabEG(i)%dx - tabEG(i)%x*tabEG(i)%dy) write(fileTripoli,*) " SURF ",i," PLAN ",a,b,c,d else xx = real(tabEG(i)%x) yy = real(tabEG(i)%y) rr = real(tabEG(i)%r) write(fileTripoli,*) " SURF ",i," CYLZ ",xx,yy,rr end if end do write(fileTripoli,*) " SURF ",szEG+1," PLANZ 0. //plan inferieur" write(fileTripoli,*) " SURF ",szEG+2," PLANZ 30. //plan superieur" !definition des volumes write(fileTripoli,'(/"//definition des volumes"/)') do i = 1,nbNode nbPlus = count(tabVolume(i)%side(1:tabVolume(i)%lg)) nbMoins = tabVolume(i)%lg - nbPlus allocate(sp(nbPlus),sm(nbMoins),stat=alloc_ok) if (alloc_ok /= 0) call XABORT("G2S: generateTripoliFile(1) => allocation pb") sp = pack(tabVolume(i)%indElem(1:tabVolume(i)%lg), & tabVolume(i)%side(1:tabVolume(i)%lg)) sm = pack(tabVolume(i)%indElem(1:tabVolume(i)%lg), & .not.tabVolume(i)%side(1:tabVolume(i)%lg)) write(fileTripoli,*) " VOLU" write(fileTripoli,*) " ",i write(fileTripoli,*) " EQUA" write(fileTripoli,*) " PLUS",nbPlus+1 write(fileTripoli,*) " ",szEG+1,sp(:) !plan inf en premier write(fileTripoli,*) " MOINS",nbMoins+1 write(fileTripoli,*) " ",sm(:),szEG+2 !plan sup en dernier write(fileTripoli,*) " FINV" write(fileTripoli,*) deallocate(sp,sm) end do write(fileTripoli,'("FINGEOM")') !conditions aux limites write(fileTripoli,'(/"//conditions aux limites")') write(fileTripoli,'("LIMIT")') !calcul du nombre total de lignes de conditions limites particulieres nbCLTot = 2*nbNode do i = 1,nbNode nbCLTot = nbCLTot & + count(tabVolume(i)%typCL(1:tabVolume(i)%lg)==Tri_Refl) & + count(tabVolume(i)%typCL(1:tabVolume(i)%lg)==Tri_Trans) & + count(tabVolume(i)%typCL(1:tabVolume(i)%lg)==Tri_Cos) end do write(fileTripoli,*) nbCLTot !affichage des conditions limites do i = 1,nbNode !Reflexions write(fileTripoli,*) " ",i," REFLECTION ",szEG+1 write(fileTripoli,*) " ",i," REFLECTION ",szEG+2 nbCL = count(tabVolume(i)%typCL(1:tabVolume(i)%lg)==Tri_Refl) if (nbCL/=0) then allocate(clList(nbCL),stat=alloc_ok) if (alloc_ok /= 0) call XABORT("G2S: generateTripoliFile(2) => allocation pb") clList = pack(tabVolume(i)%indElem(1:tabVolume(i)%lg), & tabVolume(i)%typCL(1:tabVolume(i)%lg)==Tri_Refl) do j = 1,nbCL write(fileTripoli,*) " ",i," REFLECTION ",clList(j) end do deallocate(clList) end if !Translations nbCL = count(tabVolume(i)%typCL(1:tabVolume(i)%lg)==Tri_Trans) if (nbCL/=0) then allocate(clList(nbCL),stat=alloc_ok) if (alloc_ok /= 0) call XABORT("G2S: generateTripoliFile(3) => allocation pb") clList = pack(tabVolume(i)%indElem(1:tabVolume(i)%lg), & tabVolume(i)%typCL(1:tabVolume(i)%lg)==Tri_Trans) do j = 1,nbCL write(fileTripoli,*) " ",i," TRANSLATION ",clList(j) end do deallocate(clList) end if !Cosinus nbCL = count(tabVolume(i)%typCL(1:tabVolume(i)%lg)==Tri_Cos) if (nbCL/=0) then allocate(clList(nbCL),stat=alloc_ok) if (alloc_ok /= 0) call XABORT("G2S: generateTripoliFile(4) => allocation pb") clList = pack(tabVolume(i)%indElem(1:tabVolume(i)%lg), & tabVolume(i)%typCL(1:tabVolume(i)%lg)==Tri_Cos) do j = 1,nbCL write(fileTripoli,*) " ",i," COSINUS ",clList(j) end do deallocate(clList) end if end do write(fileTripoli,'("FIN_LIMIT"/)') ! insertion des donnees neutroniques des materiaux ! le milieu ne correspond au materiaux nomme "MIX_n" write(fileTripoli,'("//insertion des donnees neutroniques")') write(fileTripoli,'("//des matreriaux")') write(fileTripoli,'("FILE")') write(fileTripoli,'(" mix_data"/)') !le nom est fixe, mais pourra etre !change si besoin ###### !definition du milieu par volume write(fileTripoli,'("//definition du milieu par volume")') write(fileTripoli,'("GEOMCOMP")') allocate(listMix(nbNode),stat=alloc_ok) if (alloc_ok /= 0) call XABORT("G2S: generateTripoliFile(5) => allocation pb") listMix(:nbNode) = 0 nbMix = 0 do i = 1,nbNode if (count(listMix(:nbNode)==tabVolume(i)%mix)==0) then nbMix = nbMix + 1 listMix(nbMix) = tabVolume(i)%mix end if end do do i = 1,nbMix allocate(volMix(count(tabVolume(:)%mix==listMix(i))),stat=alloc_ok) if (alloc_ok /= 0) call XABORT("G2S: generateTripoliFile(6) => allocation pb") volMix = pack((/(j,j=1,nbNode)/),tabVolume(:)%mix==listMix(i)) write(fileTripoli,*) " MIX_"//trim(i2s(listMix(i))),size(volMix) write(fileTripoli,*) " ",volMix deallocate(volMix) end do deallocate(listMix) write(fileTripoli,'("FIN_GEOMCOMP"/)') call destroyMC() end subroutine generateTripoliFile subroutine generateMCNPFile(fileMCNP,szSA,nbNode) integer,intent(in) :: fileMCNP,szSA,nbNode integer :: i,j,k,szEG,ind,transInd,nbEGCL real :: a,b,c,d,xx,yy,rr character(len=1) :: sgn character(len=5) :: nbr,mix character(len=6) :: other,signedNbr character(len=12) :: text12 character(len=1000) :: buffer integer,dimension(:),allocatable :: tabEGCL szEG = 0 call initMC(nbNode) call prepareMCData(szSA,nbNode,szEG) !title card write(fileMCNP,'("Jeu de donnees geometriques MCNP genere pour DRAGON & &par le module G2MC")') !message en commentaire write(fileMCNP,'("C Attention ! Les densites des milieux n''etant pas & &disponibles dans les")') write(fileMCNP,'("C jeux de donnees geometriques DRAGON, il faut passer & &une moulinette pour")') write(fileMCNP,'("C remplacer les ""DENSITE_I"" presents dans ce fichier, & &par leur valeur."/)') !cell card write(fileMCNP,'("C cell card")') do i = 1,nbNode text12 = i2s(i) nbr = text12(:5) text12 = i2s(tabVolume(i)%mix) mix = text12(:5) buffer =nbr//" "//trim(mix)//" DENSITE_"//trim(mix) ind = len_trim(buffer)+1 write(buffer(ind:),*) merge(tabVolume(i)%indElem(1:tabVolume(i)%lg), & -tabVolume(i)%indElem(1:tabVolume(i)%lg), & tabVolume(i)%side(1:tabVolume(i)%lg)) , & szEG+1 , -(szEG+2) call putOn80col(fileMCNP,buffer) end do text12=i2s(nbNode+1) nbr = text12(:5) buffer = nbr//" 0" ind = 9 if (geomTyp==HexTyp.or.geomTyp==TriaTyp) then !on a des risques de geometrie non convexe ! => on prend les complementaires do i = 1,nbNode text12 = i2s(i) nbr = text12(:5) buffer(ind:) = "#"//nbr ind = ind + len_trim(nbr) + 2 end do else !la geometrie est convexe ! on prend l'union des complementaires des cotes a condition limite nbEGCL = count(tabEG(:)%limite/=Tri_Not) allocate(tabEGCL(nbEGCL),stat=alloc_ok) if (alloc_ok /= 0) call XABORT("G2S: generateTripoliFile(7) => allocation pb") tabEGCL = pack((/(i,i=1,szEG)/),tabEG(:)%limite/=Tri_Not) do i = 1,nbEGCL do j = 1,nbNode do k = 1,tabVolume(j)%lg if (tabVolume(j)%indElem(k)==tabEGCL(i)) then if (tabVolume(j)%side(k)) tabEGCL(i) = -tabEGCL(i) exit end if end do end do end do do i = 1,nbEGCL text12 = i2s(tabEGCL(i)) signedNbr = text12(:6) buffer(ind:) = trim(signedNbr)//" :" ind = ind + len_trim(signedNbr) + 3 end do buffer(ind:) = trim(i2s(-(szEG+1)))//" : "//i2s(szEG+2) deallocate(tabEGCL) end if call putOn80col(fileMCNP,buffer) !surface card write(fileMCNP,'(/"C surface card")') c = 0. do i = 1,szEG text12 = i2s(i) nbr = text12(:5) sgn = ' ' other = ' ' select case(tabEG(i)%limite) case(Tri_Refl) sgn = '*' case(Tri_Trans) transInd = findParalleleWithTrans(i,szEG) text12=i2s(transInd) if (isSameWayEG(tabEG(i),tabEG(transInd))) then other = '-'//text12(:6) else other = text12(:6) end if case(Tri_Cos) sgn = '+' end select if (tabEG(i)%isPlan) then a = real(tabEG(i)%dy - tabEG(i)%y) b = real(tabEG(i)%x - tabEG(i)%dx) d = real(tabEG(i)%y*tabEG(i)%dx - tabEG(i)%x*tabEG(i)%dy) buffer = adjustl(sgn//nbr//trim(" "//other)//" P") ind = len_trim(buffer) + 1 write(buffer(ind:),*) a,b,c,d else xx = real(tabEG(i)%x) yy = real(tabEG(i)%y) rr = real(tabEG(i)%r) buffer = adjustl(sgn//nbr//" C/Z") ind = len_trim(buffer) + 1 write(buffer(ind:),*) xx,yy,rr end if call putOn80col(fileMCNP,buffer) end do text12 = i2s(szEG+1) nbr = text12(:5) buffer = "*"//nbr//" PZ 0. $plan inferieur" call putOn80col(fileMCNP,buffer) text12 = i2s(szEG+2) nbr = text12(:5) buffer = "*"//nbr//" PZ 30. $plan superieur" call putOn80col(fileMCNP,buffer) !data card call destroyMC() end subroutine generateMCNPFile subroutine generateSerpentFile(fileSerpent,szSA,nbNode) integer,intent(in) :: fileSerpent,szSA,nbNode integer :: i,j,k,szEG,ind,transInd,nbEGCL,nb2,isurf,jjj real :: a,b,c,d,xx,yy,rr,cuboid(6) character(len=1) :: sgn character(len=5) :: nbr,mix character(len=6) :: other character(len=12) :: text12 character(len=1000) :: buffer integer,dimension(:),allocatable :: tabEGCL,tmpnod,tmpnod2 logical :: lcuboid szEG = 0 call initMC(nbNode) call prepareMCData(szSA,nbNode,szEG) !title card write(fileSerpent,'("% Serpent combinatorial geometry generated by module G2MC:")') !cell card nbEGCL = count(tabEG(:)%limite/=Tri_Not) allocate(tabEGCL(nbEGCL),stat=alloc_ok) if (alloc_ok /= 0) call XABORT("G2S: generateSerpentFile(1) => allocation pb") tabEGCL = pack((/(i,i=1,szEG)/),tabEG(:)%limite/=Tri_Not) do i = 1,nbEGCL do j = 1,nbNode do k = 1,tabVolume(j)%lg if (tabVolume(j)%indElem(k)==tabEGCL(i)) then if (tabVolume(j)%side(k)) tabEGCL(i) = -tabEGCL(i) exit end if end do end do end do write(fileSerpent,'("% cell cards")') nb2=0 cuboid=(/ 0.0, 0.0, 0.0, 0.0, -1.0E10, 1.0E10 /) do i = 1,nbNode allocate(tmpnod(tabVolume(i)%lg),tmpnod2(tabVolume(i)%lg),stat=alloc_ok) if (alloc_ok /= 0) call XABORT("G2S: generateSerpentFile(2) => allocation pb") text12 = i2s(i) nbr = text12(:5) text12 = i2s(tabVolume(i)%mix) mix = text12(:5) buffer ="cell "//nbr//" 0 MIX_"//trim(mix) ind = len_trim(buffer)+1 tmpnod(:tabVolume(i)%lg)=merge(-tabVolume(i)%indElem(1:tabVolume(i)%lg), & tabVolume(i)%indElem(1:tabVolume(i)%lg), & tabVolume(i)%side(1:tabVolume(i)%lg)) jjj=0 lcuboid=.true. do j = 1,tabVolume(i)%lg do k =1,nbEGCL if(tmpnod(j) == tabEGCL(k)) then isurf=abs(tabEGCL(k)) if (.not.tabEG(isurf)%isPlan) call XABORT('g2s_generatingMC: unsupported BC') a = real(tabEG(isurf)%dy - tabEG(isurf)%y) b = real(tabEG(isurf)%x - tabEG(isurf)%dx) d = real(tabEG(isurf)%y*tabEG(isurf)%dx - tabEG(isurf)%x*tabEG(isurf)%dy) if ((a == 0.0).and.(tabEGCL(k) < 0)) then cuboid(4)=-d/b else if ((a == 0.0).and.(tabEGCL(k) > 0)) then cuboid(3)=-d/b else if ((b == 0.0).and.(tabEGCL(k) < 0)) then cuboid(1)=-d/a else if ((b == 0.0).and.(tabEGCL(k) > 0)) then cuboid(2)=-d/a endif if(lcuboid) then jjj=jjj+1 tmpnod2(jjj)=-szEG-1 lcuboid=.false. endif go to 10 endif end do jjj=jjj+1 tmpnod2(jjj)=-tmpnod(j) 10 continue end do if(jjj.gt.0) then write(buffer(ind:),*) tmpnod2(:jjj) write(fileSerpent,'(a)') trim(buffer) endif deallocate(tmpnod2,tmpnod) end do if (nbEGCL > 0) then text12 = i2s(nbNode+1) nbr = text12(:5) text12 = i2s(szEG+1) buffer ="cell "//nbr//" 0 outside "//trim(text12) write(fileSerpent,'(a)') trim(buffer) endif !surface card write(fileSerpent,'(/"% surface cards")') c = 0. do i = 1,szEG do k = 1,nbEGCL if (i == abs(tabEGCL(k))) go to 20 end do text12 = i2s(i) nbr = text12(:5) sgn = ' ' other = ' ' select case(tabEG(i)%limite) case(Tri_Refl) sgn = '*' case(Tri_Trans) transInd = findParalleleWithTrans(i,szEG) text12=i2s(transInd) if (isSameWayEG(tabEG(i),tabEG(transInd))) then other = '-'//text12(:6) else other = text12(:6) end if case(Tri_Cos) sgn = '+' end select if (tabEG(i)%isPlan) then a = real(tabEG(i)%dy - tabEG(i)%y) b = real(tabEG(i)%x - tabEG(i)%dx) d = -real(tabEG(i)%y*tabEG(i)%dx - tabEG(i)%x*tabEG(i)%dy) buffer = adjustl("surf "//sgn//nbr//trim(" "//other)//" plane") ind = len_trim(buffer) + 1 write(buffer(ind:),*) a,b,c,d else xx = real(tabEG(i)%x) yy = real(tabEG(i)%y) rr = real(tabEG(i)%r) buffer = adjustl("surf "//sgn//nbr//" cyl") ind = len_trim(buffer) + 1 write(buffer(ind:),*) xx,yy,rr end if write(fileSerpent,'(a)') trim(buffer) 20 continue end do if (nbEGCL > 0) then text12 = i2s(szEG+1) nbr = text12(:5) buffer = adjustl("surf "//sgn//nbr//" cuboid ") ind = len_trim(buffer) + 1 write(buffer(ind:),*) (cuboid(i),i=1,6) write(fileSerpent,'(a)') trim(buffer) endif deallocate(tabEGCL) !data card call destroyMC() end subroutine generateSerpentFile subroutine putOn80col(outFile,strIn) integer,intent(in) :: outFile character(len=*),intent(in) :: strIn integer :: lg,lastBlank character(len=15) :: frm character(len=80) :: truncStr character(len=len(strIn)) :: str str = strIn do lg = len_trim(str) if (lg <= 80) then frm = "(a"//i2s(lg)//")" write(outFile,frm) str return end if truncStr = str(1:78) do lastBlank = 78,1,-1 if (truncStr(lastBlank:lastBlank)==' ') exit end do truncStr = trim(str(1:lastBlank-1))//' &' frm = "(a"//i2s(len_trim(truncStr))//")" write(outFile,frm) truncStr str = ' '//adjustl(str(lastBlank:)) end do end subroutine putOn80col function findParalleleWithTrans(ind,szEG) integer,intent(in) :: ind,szEG integer :: findParalleleWithTrans integer :: i findParalleleWithTrans = 0 do i = 1,szEG if (i==ind) cycle if (.not.tabEG(i)%isPlan) cycle if (tabEG(i)%limite/=Tri_Trans) cycle if (.not.estColi(tabEG(i)%dx-tabEG(i)%x,tabEG(i)%dy-tabEG(i)%y, & tabEG(ind)%dx-tabEG(ind)%x,tabEG(ind)%dy-tabEG(ind)%y)) cycle findParalleleWithTrans = i return end do call XABORT("G2MC: error, no parallel side found for translation & &boundary condition") end function findParalleleWithTrans end module monteCarlo