The TELEMAC-MASCARET system  trunk
section_khione.f
Go to the documentation of this file.
1 ! *************************
2  SUBROUTINE section_khione
3 ! *************************
4  &(mesh,ikle,nelmax,ifabor,f,u,v,qvc,ca,cv,h,s,lt)
5 !
6 !***********************************************************************
7 ! TELEMAC2D V8P2 23/06/2020
8 !***********************************************************************
9 !
10 !brief Computes caracteristics of the clogged section
11 !+ Inspired by the algorithm of flusec_telemac2d
12 !
13 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
14 !| CA |-->| WET AREA OF THE SECTION
15 !| CV |-->| AVERAGE CONCENTRATION ON THE SECTION
16 !| F |-->| TRACER VALUES
17 !| H |-->| WATER DEPTH
18 !| IFABOR |-->| TABLE OF ELEMENT NEIGHBORS
19 !| IKLE |-->| CONNECTIVITY TABLE
20 !| LT |-->| NUMBER OF TIME STEPS
21 !| MESH |-->| MESH STRUCTURE
22 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENT
23 !| QVC |-->| COMPUTED FRAZIL FLUX ACROSS THE SECTION
24 !| S |-->| VOID STRUCTURE
25 !| U |-->| X-COMPONENT OF THE VELOCITY
26 !| V |-->| Y-COMPONENT OF THE VELOCITY
27 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
28 !
29  USE bief
30 !
31  USE declarations_khione, ONLY: nseclog,seclog,msksec,t1,t2,t3,
32  & bm1,bm2,cv1,cv2,cv3,cv4,ind_fra,
33  & clog_tlgth,nfrclog,t4,t5,un1,un2,
34  & nseg_fs,liste_fs
35  USE declarations_special
36  USE interface_parallel, ONLY : p_dsum,p_dmin,p_dmax,p_imin
37  IMPLICIT NONE
38  INTEGER, INTENT(IN) :: NELMAX,LT
39  INTEGER, INTENT(IN) :: IKLE(nelmax,*),IFABOR(nelmax,*)
40  DOUBLE PRECISION, INTENT(INOUT) :: QVC(nseclog)
41  DOUBLE PRECISION, INTENT(INOUT) :: CA(nseclog),CV(nseclog)
42  TYPE(bief_obj), INTENT(IN) :: F,U,V,H,S
43  TYPE(bief_mesh), INTENT(INOUT) :: MESH
44 !
45 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
46 !
47  INTEGER, PARAMETER :: NSEMAX = 200
48 !
49  INTEGER NCP,ISEC,DEP,ARR,PT,IELEM,I1,I2,I3,ELBEST,IGBEST,ILBEST
50  INTEGER ILPREC,ISEG,J1,J2,J3,ERR,IEL,IELMH,IELMU,II
51  DOUBLE PRECISION SUR6,DIST,DIST1,DIST2,DIST3,X1,X2
52  DOUBLE PRECISION Y1,Y2,NORM,L2(nseclog),DTMP1,DTMP2
53 !
54  LOGICAL OK
55 !
56 !-----------------------------------------------------------------------
57 !
58  CALL os('X=C ', x=t3, c=1.d0)
59  CALL os('X=Y ', x=t2, y=h )
60  CALL os('X=YZ ', x=t1, y=h, z=f%ADR(ind_fra)%P )
61  sur6=1.d0/6.d0
62  ncp=nseclog*2
63 !
64 ! SEARCHES PATHS JOINING POINT PAIRS:
65 !
66 !
67  IF(lt.EQ.1) THEN
68 
69  ALLOCATE(nseg_fs(ncp) ,stat=err)
70  ALLOCATE(liste_fs(ncp,nsemax,2) ,stat=err)
71  CALL check_allocate(err, "SEC_KHIONE:LISTE_FS")
72 !
73  DO isec =1,nseclog
74 ! IN THE SERIAL CASE, OR "CLASSICAL" IN PARALLEL,
75 ! FOLLOW THE ALGORITHM OF FINDING SEGMENT CHAINS
76 !
77 ! NOTE: IF YOU CHANGE THE ALGORITHM, CHANGE IT IN PARTEL AS WELL
78 !
79  dep = seclog(1+2*(isec-1))
80  arr = seclog(2+2*(isec-1))
81  IF(ncsize.GT.1) THEN
82  dep=global_to_local_point(dep,mesh)
83  arr=global_to_local_point(arr,mesh)
84  ENDIF
85  IF(dep.GT.0) THEN
86  x1 = mesh%X%R(dep)
87  y1 = mesh%Y%R(dep)
88  ELSE
89  x1 = 0.d0
90  y1 = 0.d0
91  ENDIF
92  IF(arr.GT.0) THEN
93  x2 = mesh%X%R(arr)
94  y2 = mesh%Y%R(arr)
95  ELSE
96  x2 = 0.d0
97  y2 = 0.d0
98  ENDIF
99  x1 = p_dsum(x1)
100  x2 = p_dsum(x2)
101  y1 = p_dsum(y1)
102  y2 = p_dsum(y2)
103  norm = sqrt((x1-x2)**2+(y1-y2)**2)
104  IF(norm.LE.0.d0) THEN
105  WRITE(lu,*) 'SECTION_KHIONE: SECTION TOO SHORT '
106  WRITE(lu,*) ' PROBABLY BECAUSE N_START=N_END'
107  CALL plante(1)
108  stop
109  ENDIF
110  clog_tlgth(isec+nfrclog) = norm
111  un1(isec) = (y1 - y2) / norm
112  un2(isec) = (x2 - x1) / norm
113  IF(ncsize.GT.1) THEN
114  IF(dep.EQ.0.AND.arr.EQ.0) THEN
115  nseg_fs(isec)=0
116  cycle
117  ENDIF
118  IF((dep.EQ.0.AND.arr.NE.0).OR.(dep.NE.0.AND.arr.EQ.0)) THEN
119  nseg_fs(isec)=-1
120  cycle
121  ENDIF
122  ENDIF
123  pt = dep
124  iseg = 0
125  dist=(mesh%X%R(dep)-mesh%X%R(arr))**2+
126  & (mesh%Y%R(dep)-mesh%Y%R(arr))**2
127 !
128 10 CONTINUE
129 !
130 !
131  DO ielem =1,mesh%NELEM
132 !
133  i1 = ikle(ielem,1)
134  i2 = ikle(ielem,2)
135  i3 = ikle(ielem,3)
136 ! IF THE ELEMENT CONTAINS THE SELECTED NODE:
137  IF(pt.EQ.i1.OR.pt.EQ.i2.OR.pt.EQ.i3) THEN
138  dist1 = (mesh%X%R(i1)-mesh%X%R(arr))**2 +
139  & (mesh%Y%R(i1)-mesh%Y%R(arr))**2
140  dist2 = (mesh%X%R(i2)-mesh%X%R(arr))**2 +
141  & (mesh%Y%R(i2)-mesh%Y%R(arr))**2
142  dist3 = (mesh%X%R(i3)-mesh%X%R(arr))**2 +
143  & (mesh%Y%R(i3)-mesh%Y%R(arr))**2
144  IF(dist1.LT.dist) THEN
145  dist = dist1
146  elbest = ielem
147  igbest = i1
148  ilbest = 1
149  IF(i1.EQ.pt) ilprec = 1
150  IF(i2.EQ.pt) ilprec = 2
151  IF(i3.EQ.pt) ilprec = 3
152  ENDIF
153  IF(dist2.LT.dist) THEN
154  dist = dist2
155  elbest = ielem
156  igbest = i2
157  ilbest = 2
158  IF(i1.EQ.pt) ilprec = 1
159  IF(i2.EQ.pt) ilprec = 2
160  IF(i3.EQ.pt) ilprec = 3
161  ENDIF
162  IF(dist3.LT.dist) THEN
163  dist = dist3
164  elbest = ielem
165  igbest = i3
166  ilbest = 3
167  IF(i1.EQ.pt) ilprec = 1
168  IF(i2.EQ.pt) ilprec = 2
169  IF(i3.EQ.pt) ilprec = 3
170  ENDIF
171  ENDIF
172 !
173  ENDDO ! IELEM
174 !
175  IF(igbest.EQ.pt) THEN
176  WRITE(lu,33)
177 33 FORMAT(1x,'SECTION_KHIONE : ALGORITHM FAILED')
178  CALL plante(1)
179  stop
180  ELSE
181  pt = igbest
182  iseg = iseg + 1
183  IF(iseg.GT.nsemax) THEN
184  WRITE(lu,*) 'SECTION_KHIONE: TOO MANY SEGMENTS IN A '
185  WRITE(lu,*) ' SECTION. INCREASE NSEMAX'
186  CALL plante(1)
187  stop
188  ENDIF
189  liste_fs(isec,iseg,1) = ikle(elbest,ilprec)
190  liste_fs(isec,iseg,2) = ikle(elbest,ilbest)
191  IF(igbest.NE.arr) GO TO 10
192  ENDIF
193 !
194  nseg_fs(isec) = iseg
195 !
196 ! IF COMPATIBLE COMPUTATION OF FLUXES=YES
197  IF(nseg_fs(isec).GE.1) THEN
198 !
199 ! LOOKING AT ALL ELEMENTS TOUCHING THE SECTION WITH 2 POINTS
200 ! MARKING THEM WITH 1 ON ONE SIDE AND -1 ON THE OTHER SIDE
201 !
202  DO iel=1,mesh%NELEM
203  msksec%ADR(isec)%P%R(iel)=0.d0
204  j1=ikle(iel,1)
205  j2=ikle(iel,2)
206  j3=ikle(iel,3)
207  DO iseg=1,nseg_fs(isec)
208  i1 = liste_fs(isec,iseg,1)
209  i2 = liste_fs(isec,iseg,2)
210 ! LEFT SIDE
211  IF ( (j1.EQ.i1.AND.j2.EQ.i2) .OR.
212  & (j2.EQ.i1.AND.j3.EQ.i2) .OR.
213  & (j3.EQ.i1.AND.j1.EQ.i2) ) THEN
214  msksec%ADR(isec)%P%R(iel)=1.d0
215 ! RIGHT SIDE
216  ELSEIF( (j1.EQ.i2.AND.j2.EQ.i1) .OR.
217  & (j2.EQ.i2.AND.j3.EQ.i1) .OR.
218  & (j3.EQ.i2.AND.j1.EQ.i1) ) THEN
219  msksec%ADR(isec)%P%R(iel)=-1.d0
220  ENDIF
221  ENDDO
222  ENDDO
223 !
224 ! OTHER TRIANGLES WITH ONLY 1 POINT TOUCHING THE SECTION
225 ! LOOKING AT NEIGHBOURS TO FIND THEIR SIDE
226 !
227  DO
228  ok=.true.
229  DO iel=1,mesh%NELEM
230  j1=ikle(iel,1)
231  j2=ikle(iel,2)
232  j3=ikle(iel,3)
233  DO iseg=1,nseg_fs(isec)
234  i1 = liste_fs(isec,iseg,1)
235  i2 = liste_fs(isec,iseg,2)
236  IF((j1.EQ.i1.OR.j2.EQ.i1.OR.j3.EQ.i1.OR.
237  & j1.EQ.i2.OR.j2.EQ.i2.OR.j3.EQ.i2).AND.
238  & abs(msksec%ADR(isec)%P%R(iel)).LT.0.5d0) THEN
239 ! LOOKING AT NEIGHBOURS
240  IF(ifabor(iel,1).GT.0) THEN
241  ielem=ifabor(iel,1)
242  IF(abs(msksec%ADR(isec)%P%R(ielem)).GT.0.5d0) THEN
243  msksec%ADR(isec)%P%R(iel)=
244  & msksec%ADR(isec)%P%R(ielem)
245  ENDIF
246  ENDIF
247  IF(ifabor(iel,2).GT.0) THEN
248  ielem=ifabor(iel,2)
249  IF(abs(msksec%ADR(isec)%P%R(ielem)).GT.0.5d0) THEN
250  msksec%ADR(isec)%P%R(iel)=
251  & msksec%ADR(isec)%P%R(ielem)
252  ENDIF
253  ENDIF
254  IF(ifabor(iel,3).GT.0) THEN
255  ielem=ifabor(iel,3)
256  IF(abs(msksec%ADR(isec)%P%R(ielem)).GT.0.5d0) THEN
257  msksec%ADR(isec)%P%R(iel)=
258  & msksec%ADR(isec)%P%R(ielem)
259  ENDIF
260  ENDIF
261  IF(abs(msksec%ADR(isec)%P%R(iel)).LT.0.5d0) ok=.false.
262  ENDIF
263  ENDDO
264  ENDDO
265  IF(ok) EXIT
266  ENDDO
267 !
268 ! IF(COMFLU.AND.NSEG_FS(ISEC).GE.1) THEN
269  ENDIF
270 !
271  ENDDO ! ISEC
272  ENDIF
273 !
274 !-----------------------------------------------------------------------
275 !
276  ielmh=t2%ELM
277  ielmu=u%ELM
278 !
279  DO isec = 1 , nseclog
280 !
281  qvc(isec) = 0.d0
282  ca(isec) = 0.d0
283  l2(isec)= 0.d0
284  cv(isec) = 0.d0
285 !
286  IF(nseg_fs(isec).GE.1) THEN
287 !
288  CALL os('X=C ', x=t4, c=un1(isec))
289  CALL os('X=C ', x=t5, c=un2(isec))
290 ! COMPUTING THE FLUX AS IN THE CONTINUITY EQUATION
291 ! (HOWEVER IMPLICITATION SHOULD PERHAPS BE ALSO USED)
292 !
293  CALL matrix(bm1,'M=N ','MATFGR X',ielmh,ielmu,
294  & 1.d0,t1,s,s,s,s,s,mesh,.true.,msksec%ADR(isec)%P)
295  CALL matrix(bm2,'M=N ','MATFGR Y',ielmh,ielmu,
296  & 1.d0,t1,s,s,s,s,s,mesh,.true.,msksec%ADR(isec)%P)
297 !
298  CALL matvec( 'X=AY ',cv1,bm1,u,0.d0,mesh)
299  CALL matvec( 'X=X+AY ',cv1,bm2,v,0.d0,mesh)
300  CALL matvec( 'X=AY ',cv4,bm1,t4,0.d0,mesh)
301  CALL matvec( 'X=X+AY ',cv4,bm2,t5,0.d0,mesh)
302 !
303  CALL matrix(bm1,'M=N ','MATFGR X',ielmh,ielmu,
304  & 1.d0,t2,s,s,s,s,s,mesh,.true.,msksec%ADR(isec)%P)
305  CALL matrix(bm2,'M=N ','MATFGR Y',ielmh,ielmu,
306  & 1.d0,t2,s,s,s,s,s,mesh,.true.,msksec%ADR(isec)%P)
307 !
308  CALL matvec( 'X=AY ',cv2,bm1,t4,0.d0,mesh)
309  CALL matvec( 'X=X+AY ',cv2,bm2,t5,0.d0,mesh)
310 !
311  CALL matrix(bm1,'M=N ','MATFGR X',ielmh,ielmu,
312  & 1.d0,t3,s,s,s,s,s,mesh,.true.,msksec%ADR(isec)%P)
313  CALL matrix(bm2,'M=N ','MATFGR Y',ielmh,ielmu,
314  & 1.d0,t3,s,s,s,s,s,mesh,.true.,msksec%ADR(isec)%P)
315 !
316  CALL matvec( 'X=AY ',cv3,bm1,t4,0.d0,mesh)
317  CALL matvec( 'X=X+AY ',cv3,bm2,t5,0.d0,mesh)
318 !
319 ! SUMMING CV1 FOR ALL POINTS OF THE SECTION, THIS IS THE FLUX !
320 ! (BTAINED BY CONTINUITY EQUATION AND AN INTEGRATION BY PARTS)
321 ! CV1 (flux frazil), CV2 (wet area), CV3 (length), CV4
322 ! (concentration)
323 !
324  DO iseg = 1 , nseg_fs(isec)
325  i1 = liste_fs(isec,iseg,1)
326  qvc(isec) = qvc(isec) + cv1%R(i1)
327  ca(isec) = ca(isec) + cv2%R(i1)
328  l2(isec) = l2(isec) + cv3%R(i1)
329  cv(isec) = cv(isec) + cv4%R(i1)
330  ENDDO
331 ! LAST SEGMENT, ADDING THE LAST POINT
332  i2 = liste_fs(isec,nseg_fs(isec),2)
333  qvc(isec) = qvc(isec) + cv1%R(i2)
334  ca(isec) = ca(isec) + cv2%R(i2)
335  l2(isec) = l2(isec) + cv3%R(i2)
336  cv(isec) = cv(isec) + cv4%R(i2)
337 !
338 ! WHEN BOTH UPWIND AND DOWNSTREAM ELEMENTS ARE TAKEN INTO ACCOUNT
339 ! WITH DIFFERENT SIGNS, WE GET TWICE THE FLUX
340 !
341  qvc(isec)=abs(qvc(isec)*0.5d0)/clog_tlgth(isec+nfrclog)
342  l2(isec) = abs(l2(isec)*0.5d0)
343  ca(isec)=abs(ca(isec)*0.5d0)/l2(isec)*clog_tlgth(isec+nfrclog)
344  cv(isec)=abs(cv(isec)*0.5d0)/l2(isec)*clog_tlgth(isec+nfrclog)
345 !
346  ENDIF
347  IF (ncsize.GT.1) THEN
348  ii = p_imin(nseg_fs(isec))
349  IF(ii.GE.0) THEN
350  dtmp1 = p_dmin(qvc(isec))
351  dtmp2 = p_dmax(qvc(isec))
352  qvc(isec) = dtmp1+dtmp2
353  dtmp1 = p_dmin(ca(isec))
354  dtmp2 = p_dmax(ca(isec))
355  ca(isec) = dtmp1+dtmp2
356  dtmp1 = p_dmin(cv(isec))
357  dtmp2 = p_dmax(cv(isec))
358  cv(isec) = dtmp1+dtmp2
359  ENDIF
360  ENDIF
361 !
362  ENDDO ! ISEC
363 !
364 !-----------------------------------------------------------------------
365 !
366  RETURN
367  END
interface_parallel
Definition: interface_parallel.f:3
declarations_khione::t2
type(bief_obj), pointer t2
Definition: declarations_khione.f:1090
declarations_special
Definition: declarations_special.F:3
p_dmax
double precision function p_dmax(MYPART)
Definition: p_dmax.F:7
declarations_khione::seclog
integer, dimension(:), allocatable seclog
Definition: declarations_khione.f:537
section_khione
subroutine section_khione(MESH, IKLE, NELMAX, IFABOR, F, U, V, QVC, CA, CV, H, S, LT)
Definition: section_khione.f:6
declarations_special::lu
integer lu
Definition: declarations_special.F:45
declarations_khione::nseclog
integer nseclog
Definition: declarations_khione.f:536
global_to_local_point
integer function global_to_local_point(IPOIN, MESH)
Definition: global_to_local_point.f:7
declarations_khione::msksec
type(bief_obj), target msksec
Definition: declarations_khione.f:551
declarations_khione
Definition: declarations_khione.f:3
declarations_khione::t1
type(bief_obj), pointer t1
Definition: declarations_khione.f:1090
matrix
subroutine matrix(M, OP, FORMUL, IELM1, IELM2, XMUL, F, G, H, U, V, W, MESH, MSK, MASKEL)
Definition: matrix.f:7
p_dsum
double precision function p_dsum(MYPART)
Definition: p_dsum.F:7
p_dmin
double precision function p_dmin(MYPART)
Definition: p_dmin.F:7
os
subroutine os(OP, X, Y, Z, C, IOPT, INFINI, ZERO)
Definition: os.f:7
matvec
subroutine matvec(OP, X, A, Y, C, MESH, LEGO)
Definition: matvec.f:7
declarations_khione::t3
type(bief_obj), pointer t3
Definition: declarations_khione.f:1090
p_imin
integer function p_imin(MYPART)
Definition: p_imin.F:7
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:34:45 for The TELEMAC-MASCARET system by doxygen 1.8.13