1 | include(dom.inc)
2 |
3 | SUBROUTINE SPATIAL_SCHEME(I_SCHEME, Lb, Gi, G, k_abs, k_scat, Qr, &
4 | & H, alpha, Qr_p, a_Lo, ifreq, iquad, &
5 | & nthread, npart)
6 |
7 | USE mod_slave
8 | USE mod_pmm
9 | !$ USE OMP_LIB
10 |
11 | IMPLICIT NONE
12 |
13 | include 'dom_constants.h'
14 |
15 | ! IN
16 | DOM_INT :: I_SCHEME, iquad, ifreq
17 | DOM_INT :: nthread, npart, ndir
18 | DOM_REAL :: alpha
19 | DOM_REAL,DIMENSION(is_ncells) :: Lb, Gi, k_abs, k_scat
20 | DOM_REAL,DIMENSION(is_nbfaces) :: a_Lo
21 |
22 | ! OUT
23 | DOM_REAL,DIMENSION(is_ncells) :: G
24 | DOM_REAL,DIMENSION(3,is_ncells) :: Qr
25 | DOM_REAL,DIMENSION(is_nbfaces) :: H
26 | DOM_REAL,DIMENSION(3,is_nprobes) :: Qr_p
27 |
28 | ! LOCAL
29 | DOM_INT :: ibnd, ivirt, ierr, idir_tot, ithread
30 | DOM_INT :: idir, icell, i, j, k
31 | DOM_INT :: neighfaceid, neighcellid
32 | DOM_REAL :: nmin, Lpi, dot_prod
33 | DOM_REAL,DIMENSION(is_nfacesmax) :: Dm_ij, Le
34 | DOM_REAL,DIMENSION(is_nbfaces) :: localH
35 | DOM_REAL,DIMENSION(is_ntot_vfaces, is_ntotdir) :: local_Lvirt
36 | DOM_REAL,DIMENSION(is_nfacesmax,is_ncells,nthread):: Li
37 |
38 | DOM_REAL :: L_IN, AREA_IN
39 | DOM_REAL :: SD
40 | DOM_REAL :: Y0, Y2
41 | DOM_REAL :: beta, omega, tau, KSI
42 |
43 | IF((i_inter.eq.1).and.(pmm_rank.eq.0)) THEN
44 | print*
45 | WRITE(*,'(A14XI2XA9I4)') "---> iquad:", iquad, &
46 | & " - ifreq:", ifreq
47 | ENDIF
48 |
49 | ! --------------------------------------------------!
50 | ! Update faces "in" lumminance for reflection !
51 | ! --------------------------------------------------!
52 |
53 | DO ibnd = 1,is_nbfaces
54 |
55 | j = ts_boundary(ibnd)%icell
56 | i = ts_boundary(ibnd)%iface
57 | Li(i,j,:) = a_Lo(ibnd) + (1.-ts_boundary(ibnd)%epsil)*H(ibnd)/pi
58 |
59 | ENDDO
60 |
61 | ! ----------------------------------------------!
62 | ! Initialise vectors and rounding error catcher !
63 | ! ----------------------------------------------!
64 |
65 | H = 0.
66 | G = 0.
67 | Qr = 0.
68 | Qr_p = 0.
69 |
70 | localH = 0.
71 | local_Lvirt = 0.
72 |
73 | nmin = 0.
74 | dot_prod = 0.
75 | ithread = 1
76 |
77 | ! ----------------------------------!
78 | ! Looping over all slave directions !
79 | ! ----------------------------------!
80 |
81 | !$OMP PARALLEL DO &
82 | !$OMP& PRIVATE(idir, idir_tot, icell ,i, j ,k, dot_prod, Dm_ij, &
83 | !$OMP& neighcellid, neighfaceid, Lpi, Le, ivirt, ibnd, ithread,&
84 | !$OMP& L_IN, AREA_IN, SD, Y0, Y2, beta, omega, tau, KSI) &
85 | !$OMP& SHARED(I_SCHEME, iquad, ifreq, alpha, localH,local_Lvirt, &
86 | !$OMP& Lb, Gi, k_abs, k_scat, G, Qr, Qr_p, nmin, Li)
87 |
88 | DO idir=is_dirbeg, is_dirend, 1
89 |
90 | !$ ithread=OMP_GET_THREAD_NUM()+1
91 |
92 | ! ------------------------------!
93 | ! Update virtual boundary faces !
94 | ! ------------------------------!
95 | idir_tot = is_dird + idir - 1
96 |
97 | DO ivirt = 1, is_ntot_vfaces
98 |
99 | j = ts_virtbound(ivirt)%icell_next
100 | i = ts_virtbound(ivirt)%iface_next
101 | IF(j.ne.0) THEN
102 | Li(i,j, ithread) = s_Lvirt(ivirt, idir_tot, iquad, ifreq)
103 | ENDIF
104 | ENDDO
105 |
106 | ndir = idir_tot * is_ntask
107 | IF(is_ntask.gt.24) ndir = 24
108 | IF((i_inter.eq.1).and.(pmm_rank.eq.0)) &
109 | & WRITE(*,'(A13X1I3X3F10.6)') "---- dir: ", ndir, &
110 | & s_mu(idir), s_eta(idir), s_ksi(idir)
111 | ! IF (pmm_rank.eq.0) print*, " pathw:", is_pathway(:,idir)
112 |
113 | DO icell=1, is_ncells
114 |
115 | ibnd = 1
116 | ivirt = 1
117 |
118 | ! j = icell
119 | j = is_pathway(icell,idir)
120 |
121 | ! IF (pmm_rank.eq.0) print*, " ++ Cell",icell,"/",is_ncells, &
122 | ! & ":",j
123 |
124 | ! --------------------------------------------!
125 | ! Calculate the direction.normal at each face !
126 | ! --------------------------------------------!
127 |
128 | DO i=1,is_nfcelt(j)
129 |
130 | Dm_ij(i) = s_norm(1,i,j)*s_mu(idir) + &
131 | & s_norm(2,i,j)*s_eta(idir) + &
132 | & s_norm(3,i,j)*s_ksi(idir)
133 | ENDDO
134 |
135 | ! ----------------------------------------------!
136 | ! Calculate intensity at exit faces of the cell !
137 | ! ----------------------------------------------!
138 |
139 | SELECT CASE (I_SCHEME)
140 |
141 | ! -----------------!
142 | ! MEAN FLUX SCHEME !
143 | ! -----------------!
144 | CASE(1)
145 | ! print*," (",pmm_rank,") Doing DMFS spatial discretization"
146 |
147 | L_IN = 0.
148 | AREA_IN = 0.
149 |
150 | DO i=1,is_nfcelt(j)
151 | IF (Dm_ij(i)<0.) THEN
152 | SD = s_S(i,j)*Dm_ij(i)
153 | L_IN = L_IN + SD*Li(i,j,ithread)
154 | AREA_IN = AREA_IN + SD
155 | ENDIF
156 | ENDDO
157 |
158 | Lpi = ( alpha*s_V(j)*(k_abs(j)*Lb(j) + k_scat(j)/ &
159 | & (4*pi)*Gi(j)) - L_IN ) / &
160 | & ( alpha*s_V(j)*(k_abs(j)+k_scat(j)) - AREA_IN )
161 |
162 | DO i=1,is_nfcelt(j)
163 | IF (Dm_ij(i)>0.) THEN
164 | Le(i)=( Lpi - (1.-alpha)*L_IN/AREA_IN ) / alpha
165 | IF ( Le(i).lt.0. ) Le(i) = 0.
166 | ENDIF
167 | ENDDO
168 |
169 | ! -------------------!
170 | ! EXPONENTIAL SCHEME !
171 | ! -------------------!
172 | CASE(2)
173 | ! print*," (",pmm_rank,") Doing EXP spatial discretization"
174 |
175 | Lpi = 0.
176 |
177 | beta = k_scat(j) + k_abs(j)
178 | omega = k_scat(j) / beta
179 | tau = beta*s_maxlen(icell,idir)
180 | KSI = (2./tau)*(1.-(1./tau)*(1.-exp(-tau)))
181 |
182 | Y0 = (1.-omega)*Lb(j)+omega*Gi(j)/(4*pi)
183 | Y2 = KSI*SUM(Li(:,j,ithread)*s_ss(:,icell,idir)) + &
184 | & Y0*(1-KSI)
185 |
186 | DO i=1, is_nfcelt(j)
187 | IF (Dm_ij(i)>0.) THEN
188 | Lpi = Lpi + Y2/(beta*s_V(j))*s_S(i,j)*Dm_ij(i) ! Faces de sorties
189 | ELSEIF (Dm_ij(i)<0.) THEN
190 | Lpi = Lpi + Li(i,j,ithread)/(beta*s_V(j))*s_S(i,j)* & ! Faces d'entrees
191 | & Dm_ij(i)
192 | ENDIF
193 | ENDDO
194 |
195 | Lpi = Y0 - Lpi
196 |
197 | END SELECT
198 |
199 |
200 | G(j) = G(j) + Lpi*s_w(idir)
201 |
202 | Qr(1,j) = Qr(1,j) + Lpi*s_w(idir)*s_mu (idir)
203 | Qr(2,j) = Qr(2,j) + Lpi*s_w(idir)*s_eta(idir)
204 | Qr(3,j) = Qr(3,j) + Lpi*s_w(idir)*s_ksi(idir)
205 |
206 | ! IF (pmm_rank.eq.0) print*, " Ip =",Lpi
207 |
208 | ! -----------------!
209 | ! Probe processing !
210 | ! -----------------!
211 |
212 | IF (is_nprobes.gt.0) THEN
213 | IF (is_pcells(j).gt.0 ) THEN
214 | DO k=1,is_nprobes
215 | IF(is_norm_probe(k) .eq.j ) THEN
216 | dot_prod = (s_mu(idir)*s_norm_probe(1,k)+s_eta(idir)* &
217 | & s_norm_probe(2,k)+s_ksi(idir)*s_norm_probe(3,k))
218 |
219 | IF(dot_prod.ge.nmin) THEN
220 | Qr_p(1,k) = Qr_p(1,k) + Lpi*s_w(idir)*dot_prod
221 | ELSE
222 | Qr_p(2,k) = Qr_p(2,k) + Lpi*s_w(idir)*dot_prod
223 | ENDIF
224 | Qr_p(3,k) = Qr(1,j)
225 | ENDIF
226 | ENDDO
227 | ENDIF
228 | ENDIF
229 |
230 | ! -----------------------------------------------------------!
231 | ! Update incident intensities at the wall and neighbor cells !
232 | ! -----------------------------------------------------------!
233 |
234 | DO i=1,is_nfcelt(j)
235 |
236 | IF (Dm_ij(i).ge.nmin) THEN
237 |
238 | ! ----------------------------------!
239 | ! Detect if face is in the boundary !
240 | ! ----------------------------------!
241 |
242 | neighcellid = is_neighs(2*i-1,j)
243 |
244 | ! ----------------------------------------------!
245 | ! Apply luminance to the neighbor/boundary face !
246 | ! ----------------------------------------------!
247 |
248 | IF (neighcellid.gt.0) THEN
249 |
250 | neighfaceid = is_neighs(2*i,j)
251 | Li(neighfaceid, neighcellid, ithread) = Le(i)
252 |
253 | ELSEIF (neighcellid.eq.-1) THEN
254 |
255 | CALL FIND_IVIRT(j,i,ivirt)
256 | local_Lvirt(ivirt, idir_tot) = Le(i)
257 |
258 | ELSEIF (neighcellid.eq.0) THEN
259 |
260 | ! IF (pmm_rank.eq.0) print*, " face:", i
261 | CALL FIND_IBND(j,i,ibnd)
262 | localH(ibnd)= localH(ibnd) + Le(i)*s_w(idir)*Dm_ij(i)
263 |
264 | ENDIF
265 |
266 | ENDIF
267 |
268 | ENDDO
269 | ENDDO
270 | ENDDO
271 | !$OMP END PARALLEL DO
272 |
273 | IF(is_ntask.gt.1) THEN
274 | CALL MPI_ALLREDUCE(localH, H, is_nbfaces, &
275 | & MPI_DOUBLE_PRECISION, MPI_SUM, SUB_COMM, &
276 | & ierr)
277 | ELSE
278 | H = localH
279 | ENDIF
280 |
281 | IF(npart.gt.1) &
282 | & CALL MPI_ALLREDUCE(local_Lvirt, s_Lvirt(:,:,iquad, ifreq), &
283 | & is_ntot_vfaces*is_ntotdir, &
284 | & MPI_DOUBLE_PRECISION, MPI_SUM, COMM_PARA, &
285 | & ierr)
286 |
287 | END SUBROUTINE SPATIAL_SCHEME
spatial_scheme.F could be called by: