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)
5 |
6 | USE mod_slave
7 | USE mod_pmm
8 |
9 | IMPLICIT NONE
10 |
11 | include 'dom_constants.h'
12 |
13 | ! IN
14 | DOM_INT :: I_SCHEME, t,u
15 | DOM_REAL :: alpha
16 | DOM_REAL,DIMENSION(is_ncells) :: Lb
17 | DOM_REAL,DIMENSION(is_ncells) :: k_abs, k_scat
18 | DOM_REAL,DIMENSION(is_ncells) :: Gi
19 | DOM_REAL,DIMENSION(is_nbfaces) :: a_Lo
20 |
21 |
22 | ! OUT
23 | DOM_REAL,DIMENSION(is_ncells) :: G
24 | DOM_REAL,DIMENSION(3,is_ncells) :: Qr
25 | DOM_REAL,DIMENSION(is_nfacesmax,is_ncells) :: Li
26 | DOM_REAL,DIMENSION(is_nbfaces) :: H
27 | DOM_REAL,DIMENSION(3,is_nprobes) :: Qr_p
28 |
29 | ! LOCAL
30 | DOM_INT :: inode, ibnd, ierr
31 | DOM_INT :: idir,icell,i,j,k
32 | DOM_INT :: neighfaceid, neighcellid
33 | DOM_REAL :: nmin
34 | DOM_REAL :: Lpi
35 | DOM_REAL :: dot_prod
36 | DOM_REAL,DIMENSION(is_nfacesmax) :: Dm_ij
37 | DOM_REAL,DIMENSION(is_nbfaces) :: localH
38 | DOM_REAL,DIMENSION(is_nfacesmax) :: Le
39 |
40 | ! --------------------------------------------------!
41 | ! Update faces "in" lumminance for reflection !
42 | ! --------------------------------------------------!
43 |
44 | DO ibnd = 1,is_nbfaces
45 |
46 | j = is_bcell(ibnd)
47 | i = is_bface(ibnd)
48 | Li(i,j) = a_Lo(ibnd) + (1.-s_epsil(ibnd))*H(ibnd)/pi
49 |
50 | ENDDO
51 |
52 | ! ----------------------------------------------!
53 | ! Initialise vectors and rounding error catcher !
54 | ! ----------------------------------------------!
55 |
56 | H = 0.
57 | G = 0.
58 | Qr = 0.
59 | nmin = 0.
60 |
61 | Qr_p = 0.
62 | dot_prod = 0.
63 | localH = 0.
64 |
65 | ! ----------------------------------!
66 | ! Looping over all slave directions !
67 | ! ----------------------------------!
68 |
69 | DO idir=is_dirbeg, is_dirend, 1
70 |
71 | ! IF (pmm_rank.eq.0) print*, "--- dir:", &
72 | ! & s_mu(idir),s_eta(idir),s_ksi(idir)
73 | ! IF (pmm_rank.eq.0) print*, " pathw:", is_pathway(:,idir)
74 |
75 | DO icell=1, is_ncells
76 |
77 | ibnd = 1
78 | j = is_pathway(icell,idir)
79 |
80 | ! IF (pmm_rank.eq.0) print*, " ++ Cell",icell,"/",is_ncells, &
81 | ! & ":",j
82 |
83 | ! --------------------------------------------!
84 | ! Calculate the direction.normal at each face !
85 | ! --------------------------------------------!
86 |
87 | DO i=1,is_nfcelt(j)
88 |
89 | Dm_ij(i) = s_norm(1,i,j)*s_mu(idir) + &
90 | & s_norm(2,i,j)*s_eta(idir) + &
91 | & s_norm(3,i,j)*s_ksi(idir)
92 | ENDDO
93 |
94 | ! ----------------------------------------------!
95 | ! Calculate intensity at exit faces of the cell !
96 | ! ----------------------------------------------!
97 |
98 | SELECT CASE (I_SCHEME)
99 |
100 | ! -----------------!
101 | ! MEAN FLUX SCHEME !
102 | ! -----------------!
103 | CASE(1)
104 | ! print*," (",pmm_rank,") Doing DMFS spatial discretization"
105 |
106 | CALL MFSCHEME(k_abs(j), k_scat(j), s_V(j), is_nfcelt(j), &
107 | & s_S(:,j), Dm_ij, Lb(j), Gi(j), Li(:,j), &
108 | & Le(:), Lpi, alpha, is_nfacesmax)
109 |
110 | ! -------------------!
111 | ! EXPONENTIAL SCHEME !
112 | ! -------------------!
113 | CASE(2)
114 | ! print*," (",pmm_rank,") Doing EXP spatial discretization"
115 |
116 | CALL EXPOSCHEME(k_abs(j), k_scat(j), s_V(j), &
117 | & is_nfcelt(j), s_S(:,j), Dm_ij, Lb(j), &
118 | & s_maxlen(icell,idir), s_ss(:,icell,idir), &
119 | & Gi(j), Li(:,j), Le(:), Lpi,is_nfacesmax)
120 |
121 | END SELECT
122 |
123 |
124 | G(j) = G(j) + Lpi*s_w(idir)
125 |
126 | Qr(1,j) = Qr(1,j) + Lpi*s_w(idir)*s_mu (idir)
127 | Qr(2,j) = Qr(2,j) + Lpi*s_w(idir)*s_eta(idir)
128 | Qr(3,j) = Qr(3,j) + Lpi*s_w(idir)*s_ksi(idir)
129 |
130 | ! IF (pmm_rank.eq.0) print*, " Ip =",Lpi
131 |
132 | ! -----------------!
133 | ! Probe processing !
134 | ! -----------------!
135 |
136 | IF (is_nprobes.gt.0) THEN
137 | IF (is_pcells(j).gt.0 ) THEN
138 | DO k=1,is_nprobes
139 | IF(is_norm_probe(k) .eq.j ) THEN
140 | dot_prod = (s_mu(idir)*s_norm_probe(1,k)+s_eta(idir)* &
141 | & s_norm_probe(2,k)+s_ksi(idir)*s_norm_probe(3,k))
142 |
143 | IF(dot_prod.ge.nmin) THEN
144 | Qr_p(1,k) = Qr_p(1,k) + Lpi*s_w(idir)*dot_prod
145 | ELSE
146 | Qr_p(2,k) = Qr_p(2,k) + Lpi*s_w(idir)*dot_prod
147 | ENDIF
148 | Qr_p(3,k) = Qr(1,j)
149 | ENDIF
150 | ENDDO
151 | ENDIF
152 | ENDIF
153 |
154 | ! -----------------------------------------------------------!
155 | ! Update incident intensities at the wall and neighbor cells !
156 | ! -----------------------------------------------------------!
157 |
158 | DO i=1,is_nfcelt(j)
159 |
160 | IF (Dm_ij(i).ge.nmin) THEN
161 |
162 | ! ----------------------------------!
163 | ! Detect if face is in the boundary !
164 | ! ----------------------------------!
165 |
166 | neighcellid = is_neighs(2*i-1,j)
167 |
168 | ! ----------------------------------------------!
169 | ! Apply luminance to the neighbor/boundary face !
170 | ! ----------------------------------------------!
171 |
172 | IF (neighcellid.ne.0) THEN
173 |
174 | neighfaceid = is_neighs(2*i,j)
175 | Li(neighfaceid, neighcellid) = Le(i)
176 |
177 | ELSE
178 |
179 | ! IF (pmm_rank.eq.0) print*, " face:", i
180 | CALL FIND_IBND(j,i,ibnd)
181 | localH(ibnd)= localH(ibnd) + Le(i)*s_w(idir)*Dm_ij(i)
182 |
183 | ENDIF
184 |
185 | ENDIF
186 |
187 | ENDDO
188 | ENDDO
189 | ENDDO
190 |
191 | CALL MPI_ALLREDUCE(localH, H, is_nbfaces, &
192 | & MPI_DOUBLE_PRECISION, MPI_SUM, COMM_PARA, &
193 | & ierr)
194 |
195 |
196 | END SUBROUTINE SPATIAL_SCHEME
spatial_scheme.F could be called by: