diff options
Diffstat (limited to 'Dragon/src/g2s_generatingMC.f90')
| -rw-r--r-- | Dragon/src/g2s_generatingMC.f90 | 825 |
1 files changed, 825 insertions, 0 deletions
diff --git a/Dragon/src/g2s_generatingMC.f90 b/Dragon/src/g2s_generatingMC.f90 new file mode 100644 index 0000000..286fe2e --- /dev/null +++ b/Dragon/src/g2s_generatingMC.f90 @@ -0,0 +1,825 @@ +! +!----------------------------------------------------------------------- +! +!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 |