OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
alefvm_aflux3.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| alefvm_aflux3 ../engine/source/ale/alefvm/alefvm_aflux3.F
25!||--- called by ------------------------------------------------------
26!|| aflux0 ../engine/source/ale/aflux0.F
27!||--- uses -----------------------------------------------------
28!|| ale_connectivity_mod ../common_source/modules/ale/ale_connectivity_mod.F
29!|| alefvm_mod ../common_source/modules/ale/alefvm_mod.F
30!|| i22tri_mod ../common_source/modules/interfaces/cut-cell-search_mod.F
31!||====================================================================
32 SUBROUTINE alefvm_aflux3(PM ,IXS ,W ,FLUX ,
33 2 FLU1 ,ALE_CONNECT ,
34 3 IPM ,NV46 ,X ,
35 4 NEL )
36C-----------------------------------------------
37C D e s c r i p t i o n
38C-----------------------------------------------
39C 'alefvm' is related to a collocated scheme (built from FVM and based on Godunov scheme)
40C which was temporarily introduced for experimental option /INTER/TYPE22 (FSI coupling with cut cell method)
41C This cut cell method is not completed, abandoned, and is not an official option.
42C There is no other use for this scheme which is automatically enabled when /INTER/TYPE22 is defined (INT22>0 => IALEFVM=1).
43C
44C This subroutine is treating an uncut cell.
45C-----------------------------------------------
46C M o d u l e s
47C-----------------------------------------------
48 USE i22tri_mod
49 USE alefvm_mod
50 USE ale_connectivity_mod
51C-----------------------------------------------
52C I m p l i c i t T y p e s
53C-----------------------------------------------
54#include "implicit_f.inc"
55C-----------------------------------------------
56C G l o b a l P a r a m e t e r s
57C-----------------------------------------------
58#include "mvsiz_p.inc"
59C-----------------------------------------------
60C C o m m o n B l o c k s
61C-----------------------------------------------
62#include "vect01_c.inc"
63#include "param_c.inc"
64#include "com01_c.inc"
65#include "com08_c.inc"
66C-----------------------------------------------
67C D u m m y A r g u m e n t s
68C-----------------------------------------------
69 INTEGER :: IXS(NIXS,*), NV46, IPM(NPROPMI,*),NEL
70 my_real :: PM(NPROPM,*), FLUX(6,*), FLU1(*), X(3,*),W(3,*)
71 TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECT
72C-----------------------------------------------
73C L o c a l V a r i a b l e s
74C-----------------------------------------------
75 INTEGER MAT(MVSIZ), NC1(MVSIZ), NC2(MVSIZ), NC3(MVSIZ), NC4(MVSIZ), NC5(MVSIZ), NC6(MVSIZ),
76 . NC7(MVSIZ), NC8(MVSIZ), I, II,J,IMAT,IALEFVM_FLG,IDV, ICF(4,6),
77 . IX1,IX2,IX3,IX4
78 my_real
79 . x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz), x5(mvsiz), x6(mvsiz), x7(mvsiz), x8(mvsiz),
80 . y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz), y5(mvsiz), y6(mvsiz), y7(mvsiz), y8(mvsiz), z1(mvsiz),
81 . z2(mvsiz), z3(mvsiz), z4(mvsiz), z5(mvsiz), z6(mvsiz), z7(mvsiz), z8(mvsiz), n1x(mvsiz),
82 . n2x(mvsiz), n3x(mvsiz), n4x(mvsiz), n5x(mvsiz), n6x(mvsiz), n1y(mvsiz), n2y(mvsiz), n3y(mvsiz),
83 . n4y(mvsiz), n5y(mvsiz), n6y(mvsiz), n1z(mvsiz), n2z(mvsiz), n3z(mvsiz), n4z(mvsiz), n5z(mvsiz),
84 . n6z(mvsiz), flux1(mvsiz), flux2(mvsiz), flux3(mvsiz), flux4(mvsiz), flux5(mvsiz),
85 . flux6(mvsiz), vx1(mvsiz), vx2(mvsiz), vx3(mvsiz), vx4(mvsiz), vx5(mvsiz), vx6(mvsiz),
86 . vy1(mvsiz), vy2(mvsiz), vy3(mvsiz), vy4(mvsiz), vy5(mvsiz), vy6(mvsiz), vz1(mvsiz), vz2(mvsiz),
87 . vz3(mvsiz), vz4(mvsiz), vz5(mvsiz), vz6(mvsiz),
88 . reduc,upwl(6,mvsiz),valvois(6,nv46,mvsiz), valel(6,mvsiz),vec(3,6),cf(3,mvsiz),
89 . z(3),zadj(3),zcf(3),zzadj(3),zcf_,zzadj_,cos, ps, denom, denom1,denom2,lambda,xsi,wface(3,6,mvsiz),sr1,sr2,
90 . surf(6), term2, n(3,6,mvsiz)
91 LOGICAL debug_outp
92 INTEGER :: idbf, idbl, NC(8,MVSIZ), IAD2, LGTH
93C-----------------------------------------------
94 DATA icf/1,4,3,2,3,4,8,7,5,6,7,8,1,2,6,5,2,3,7,6,1,5,8,4/
95C-----------------------------------------------
96C D e s c r i p t i o n
97C-----------------------------------------------
98C This subroutine enables to compute face velocities
99C from centroid values [rho*U] stored in Fcell vector
100C (later used for centroid forces)
101C Face velocities are stored in vectors F_Face (ALEFVM_MOD)
102C This same vector is used later to store forces on faces
103C (used to compute centroid force)
104 !IPM(251,MAT(1)) = I_ALE_SOLVER
105 ! 0: Default = 1 expect if /ALE/MAT or /EULER/MAT has IFROM flag defined.
106 ! 1 : FEM
107 ! 2 : FVM U average
108 ! 3 : FVM rho.U average
109 ! 4 : FVM rho.c.U average
110 ! 5 : Godunov Acoustic
111 ! 6 : experimental
112C-----------------------------------------------
113C P r e - C o n d i t i o n s
114C-----------------------------------------------
115 IF(alefvm_param%IEnabled == 0) RETURN
116 imat = ixs(1,1+nft)
117 ialefvm_flg = ipm(251,imat)
118 IF(ialefvm_flg <= 1)RETURN
119C-----------------------------------------------
120C S o u r c e L i n e s
121C-----------------------------------------------
122
123 DO i=1,nel
124 ii = i + nft
125 mat(i) = ixs(1,ii)
126 !---8 local node numbers NC1 TO NC8 for solid element I ---!
127 nc1(i)=ixs(2,ii)
128 nc2(i)=ixs(3,ii)
129 nc3(i)=ixs(4,ii)
130 nc4(i)=ixs(5,ii)
131 nc5(i)=ixs(6,ii)
132 nc6(i)=ixs(7,ii)
133 nc7(i)=ixs(8,ii)
134 nc8(i)=ixs(9,ii)
135 !
136 !---Coordinates of the 8 nodes
137 x1(i)=x(1,nc1(i))
138 y1(i)=x(2,nc1(i))
139 z1(i)=x(3,nc1(i))
140 !
141 x2(i)=x(1,nc2(i))
142 y2(i)=x(2,nc2(i))
143 z2(i)=x(3,nc2(i))
144 !
145 x3(i)=x(1,nc3(i))
146 y3(i)=x(2,nc3(i))
147 z3(i)=x(3,nc3(i))
148 !
149 x4(i)=x(1,nc4(i))
150 y4(i)=x(2,nc4(i))
151 z4(i)=x(3,nc4(i))
152 !
153 x5(i)=x(1,nc5(i))
154 y5(i)=x(2,nc5(i))
155 z5(i)=x(3,nc5(i))
156 !
157 x6(i)=x(1,nc6(i))
158 y6(i)=x(2,nc6(i))
159 z6(i)=x(3,nc6(i))
160 !
161 x7(i)=x(1,nc7(i))
162 y7(i)=x(2,nc7(i))
163 z7(i)=x(3,nc7(i))
164 !
165 x8(i)=x(1,nc8(i))
166 y8(i)=x(2,nc8(i))
167 z8(i)=x(3,nc8(i))
168 ENDDO
169
170 !======================================================!
171 ! FACE VELOCITIES !
172 !======================================================!
173 DO i=1,nel
174 wface(1,1,i) = fourth*(w(1,nc1(i))+w(1,nc2(i))+w(1,nc3(i))+w(1,nc4(i)))
175 wface(2,1,i) = fourth*(w(2,nc1(i))+w(2,nc2(i))+w(2,nc3(i))+w(2,nc4(i)))
176 wface(3,1,i) = fourth*(w(3,nc1(i))+w(3,nc2(i))+w(3,nc3(i))+w(3,nc4(i)))
177
178 wface(1,2,i) = fourth*(w(1,nc3(i))+w(1,nc4(i))+w(1,nc7(i))+w(1,nc8(i)))
179 wface(2,2,i) = fourth*(w(2,nc3(i))+w(2,nc4(i))+w(2,nc7(i))+w(2,nc8(i)))
180 wface(3,2,i) = fourth*(w(3,nc3(i))+w(3,nc4(i))+w(3,nc7(i))+w(3,nc8(i)))
181
182
183 wface(1,3,i) = fourth*(w(1,nc5(i))+w(1,nc6(i))+w(1,nc7(i))+w(1,nc8(i)))
184 wface(2,3,i) = fourth*(w(2,nc5(i))+w(2,nc6(i))+w(2,nc7(i))+w(2,nc8(i)))
185 wface(3,3,i) = fourth*(w(3,nc5(i))+w(3,nc6(i))+w(3,nc7(i))+w(3,nc8(i)))
186
187 wface(1,4,i) = fourth*(w(1,nc1(i))+w(1,nc2(i))+w(1,nc5(i))+w(1,nc6(i)))
188 wface(2,4,i) = fourth*(w(2,nc1(i))+w(2,nc2(i))+w(2,nc5(i))+w(2,nc6(i)))
189 wface(3,4,i) = fourth*(w(3,nc1(i))+w(3,nc2(i))+w(3,nc5(i))+w(3,nc6(i)))
190
191 wface(1,5,i) = fourth*(w(1,nc2(i))+w(1,nc3(i))+w(1,nc6(i))+w(1,nc7(i)))
192 wface(2,5,i) = fourth*(w(2,nc2(i))+w(2,nc3(i))+w(2,nc6(i))+w(2,nc7(i)))
193 wface(3,5,i) = fourth*(w(3,nc2(i))+w(3,nc3(i))+w(3,nc6(i))+w(3,nc7(i)))
194
195 wface(1,6,i) = fourth*(w(1,nc1(i))+w(1,nc4(i))+w(1,nc5(i))+w(1,nc8(i)))
196 wface(2,6,i) = fourth*(w(2,nc1(i))+w(2,nc4(i))+w(2,nc5(i))+w(2,nc8(i)))
197 wface(3,6,i) = fourth*(w(3,nc1(i))+w(3,nc4(i))+w(3,nc5(i))+w(3,nc8(i)))
198 ENDDO
199
200 !======================================================!
201 ! INITIALIZATION : ADJAENT VALUES !
202 !======================================================!
203 !FCELL = [rhoU]. this quantity must be know for all groups
204 DO i=1,nel
205 ii = i + nft
206 valel(1,i) = alefvm_buffer%FCELL(1,ii) !rho.ux
207 valel(2,i) = alefvm_buffer%FCELL(2,ii) !rho.uy
208 valel(3,i) = alefvm_buffer%FCELL(3,ii) !rho.uz
209 valel(4,i) = alefvm_buffer%FCELL(4,ii) !rho
210 valel(5,i) = alefvm_buffer%FCELL(5,ii) !ssp
211 valel(6,i) = alefvm_buffer%FCELL(6,ii) !P
212 iad2 = ale_connect%ee_connect%iad_connect(ii)
213 lgth = ale_connect%ee_connect%iad_connect(ii+1) - iad2
214 DO j=1,lgth
215 idv = ale_connect%ee_connect%connected(iad2 + j - 1)
216 IF(idv <= 0) THEN
217 valvois(1,j,i) = -alefvm_buffer%FCELL(1,ii)
218 valvois(2,j,i) = -alefvm_buffer%FCELL(2,ii)
219 valvois(3,j,i) = -alefvm_buffer%FCELL(3,ii)
220 valvois(4,j,i) = alefvm_buffer%FCELL(4,ii)
221 valvois(5,j,i) = alefvm_buffer%FCELL(5,ii)
222 valvois(6,j,i) = alefvm_buffer%FCELL(6,ii)
223 ELSE
224 valvois(1,j,i) = alefvm_buffer%FCELL(1,idv)
225 valvois(2,j,i) = alefvm_buffer%FCELL(2,idv)
226 valvois(3,j,i) = alefvm_buffer%FCELL(3,idv)
227 valvois(4,j,i) = alefvm_buffer%FCELL(4,idv)
228 valvois(5,j,i) = alefvm_buffer%FCELL(5,idv)
229 valvois(6,j,i) = alefvm_buffer%FCELL(6,idv)
230 ENDIF
231 enddo!next J
232 enddo!next I
233
234 !======================================================!
235 ! NORMAL VECTORS ON EACH DACE !
236 ! 2S[n] = [diag1] x [diag2] !
237 ! where !
238 ! [n] : unitary normal vector on face !
239 !======================================================!
240 DO i=1,nel
241 ! Face-1
242 n1x(i)=(y3(i)-y1(i))*(z2(i)-z4(i)) - (z3(i)-z1(i))*(y2(i)-y4(i))
243 n1y(i)=(z3(i)-z1(i))*(x2(i)-x4(i)) - (x3(i)-x1(i))*(z2(i)-z4(i))
244 n1z(i)=(x3(i)-x1(i))*(y2(i)-y4(i)) - (y3(i)-y1(i))*(x2(i)-x4(i))
245 ! Face-2
246 n2x(i)=(y7(i)-y4(i))*(z3(i)-z8(i)) - (z7(i)-z4(i))*(y3(i)-y8(i))
247 n2y(i)=(z7(i)-z4(i))*(x3(i)-x8(i)) - (x7(i)-x4(i))*(z3(i)-z8(i))
248 n2z(i)=(x7(i)-x4(i))*(y3(i)-y8(i)) - (y7(i)-y4(i))*(x3(i)-x8(i))
249 ! Face-3
250 n3x(i)=(y6(i)-y8(i))*(z7(i)-z5(i)) - (z6(i)-z8(i))*(y7(i)-y5(i))
251 n3y(i)=(z6(i)-z8(i))*(x7(i)-x5(i)) - (x6(i)-x8(i))*(z7(i)-z5(i))
252 n3z(i)=(x6(i)-x8(i))*(y7(i)-y5(i)) - (y6(i)-y8(i))*(x7(i)-x5(i))
253 ! Face-4
254 n4x(i)=(y2(i)-y5(i))*(z6(i)-z1(i)) - (z2(i)-z5(i))*(y6(i)-y1(i))
255 n4y(i)=(z2(i)-z5(i))*(x6(i)-x1(i)) - (x2(i)-x5(i))*(z6(i)-z1(i))
256 n4z(i)=(x2(i)-x5(i))*(y6(i)-y1(i)) - (y2(i)-y5(i))*(x6(i)-x1(i))
257 ! Face-5
258 n5x(i)=(y7(i)-y2(i))*(z6(i)-z3(i)) - (z7(i)-z2(i))*(y6(i)-y3(i))
259 n5y(i)=(z7(i)-z2(i))*(x6(i)-x3(i)) - (x7(i)-x2(i))*(z6(i)-z3(i))
260 n5z(i)=(x7(i)-x2(i))*(y6(i)-y3(i)) - (y7(i)-y2(i))*(x6(i)-x3(i))
261 ! Face-6
262 n6x(i)=(y8(i)-y1(i))*(z4(i)-z5(i)) - (z8(i)-z1(i))*(y4(i)-y5(i))
263 n6y(i)=(z8(i)-z1(i))*(x4(i)-x5(i)) - (x8(i)-x1(i))*(z4(i)-z5(i))
264 n6z(i)=(x8(i)-x1(i))*(y4(i)-y5(i)) - (y8(i)-y1(i))*(x4(i)-x5(i))
265 ! Stack
266 n(1,1,i) = n1x(i)
267 n(2,1,i) = n1y(i)
268 n(3,1,i) = n1z(i)
269 !
270 n(1,2,i) = n2x(i)
271 n(2,2,i) = n2y(i)
272 n(3,2,i) = n2z(i)
273 !
274 n(1,3,i) = n3x(i)
275 n(2,3,i) = n3y(i)
276 n(3,3,i) = n3z(i)
277 !
278 n(1,4,i) = n4x(i)
279 n(2,4,i) = n4y(i)
280 n(3,4,i) = n4z(i)
281 !
282 n(1,5,i) = n5x(i)
283 n(2,5,i) = n5y(i)
284 n(3,5,i) = n5z(i)
285 !
286 n(1,6,i) = n6x(i)
287 n(2,6,i) = n6y(i)
288 n(3,6,i) = n6z(i)
289 END DO
290
291 !======================================================!
292 ! RELATIVE VELOCITIES ON EACH FACE !
293 !======================================================!
294 !UPDATE LATER WITH LINEAR INTERPOLATION INSTEAD OF MEAN VALUE
295 !MEAN VALUE OK FOR CARTESIAN REGULAR MESH
296
297 IF(alefvm_param%IFORM/=0)ialefvm_flg = alefvm_param%IFORM !for debug purpose if set in alefvm_init.F
298
299 SELECT CASE(alefvm_param%ISOLVER)
300
301
302 CASE(1)
303 !centered fvm - u average
304 DO i=1,nel
305 ii = i + nft
306 DO j=1,nv46
307 vec(1,j) = half * (valel(1,i)/valel(4,i) + valvois(1,j,i)/valvois(4,j,i))
308 vec(2,j) = half * (valel(2,i)/valel(4,i) + valvois(2,j,i)/valvois(4,j,i))
309 vec(3,j) = half * (valel(3,i)/valel(4,i) + valvois(3,j,i)/valvois(4,j,i))
310 vec(1,j) = vec(1,j) - wface(1,j,i)
311 vec(2,j) = vec(2,j) - wface(2,j,i)
312 vec(3,j) = vec(3,j) - wface(3,j,i)
313 enddo!next J
314 DO j=1,nv46
315 alefvm_buffer%F_FACE(1,j,ii) = vec(1,j)
316 alefvm_buffer%F_FACE(2,j,ii) = vec(2,j)
317 alefvm_buffer%F_FACE(3,j,ii) = vec(3,j)
318 ENDDO
319 !---face-1
320 vx1(i) = vec(1,1)
321 vy1(i) = vec(2,1)
322 vz1(i) = vec(3,1)
323 !---face-2
324 vx2(i) = vec(1,2)
325 vy2(i) = vec(2,2)
326 vz2(i) = vec(3,2)
327 !---face-3
328 vx3(i) = vec(1,3)
329 vy3(i) = vec(2,3)
330 vz3(i) = vec(3,3)
331 !---face-4
332 vx4(i) = vec(1,4)
333 vy4(i) = vec(2,4)
334 vz4(i) = vec(3,4)
335 !---face-5
336 vx5(i) = vec(1,5)
337 vy5(i) = vec(2,5)
338 vz5(i) = vec(3,5)
339 !---face-6
340 vx6(i) = vec(1,6)
341 vy6(i) = vec(2,6)
342 vz6(i) = vec(3,6)
343 enddo!next I
344
345 CASE(2)
346 !CENTERED FVM - rho.U average
347 DO i=1,nel
348 ii = i + nft
349 DO j=1,nv46
350 vec(1,j) = (valel(1,i) + valvois(1,j,i)) / (valel(4,i)+valvois(4,j,i))
351 vec(2,j) = (valel(2,i) + valvois(2,j,i)) / (valel(4,i)+valvois(4,j,i))
352 vec(3,j) = (valel(3,i) + valvois(3,j,i)) / (valel(4,i)+valvois(4,j,i))
353 vec(1,j) = vec(1,j) - wface(1,j,i)
354 vec(2,j) = vec(2,j) - wface(2,j,i)
355 vec(3,j) = vec(3,j) - wface(3,j,i)
356 enddo!next J
357 DO j=1,nv46
358 alefvm_buffer%F_FACE(1,j,ii) = vec(1,j)
359 alefvm_buffer%F_FACE(2,j,ii) = vec(2,j)
360 alefvm_buffer%F_FACE(3,j,ii) = vec(3,j)
361 ENDDO
362 !---face-1
363 vx1(i) = vec(1,1)
364 vy1(i) = vec(2,1)
365 vz1(i) = vec(3,1)
366 !---face-2
367 vx2(i) = vec(1,2)
368 vy2(i) = vec(2,2)
369 vz2(i) = vec(3,2)
370 !---face-3
371 vx3(i) = vec(1,3)
372 vy3(i) = vec(2,3)
373 vz3(i) = vec(3,3)
374 !---face-4
375 vx4(i) = vec(1,4)
376 vy4(i) = vec(2,4)
377 vz4(i) = vec(3,4)
378 !---face-5
379 vx5(i) = vec(1,5)
380 vy5(i) = vec(2,5)
381 vz5(i) = vec(3,5)
382 !---face-6
383 vx6(i) = vec(1,6)
384 vy6(i) = vec(2,6)
385 vz6(i) = vec(3,6)
386 enddo!next I
387
388 CASE(3)
389 !CENTERED FVM - Roe-average
390 DO i=1,nel
391 ii = i + nft
392 sr1 = sqrt(valel(4,i))
393 DO j=1,nv46
394 sr2 = sqrt(valvois(4,j,i))
395 vec(1,j) = (valel(1,i)/sr1 + valvois(1,j,i)/sr2) / (sr1+sr2)
396 vec(2,j) = (valel(2,i)/sr1 + valvois(2,j,i)/sr2) / (sr1+sr2)
397 vec(3,j) = (valel(3,i)/sr1 + valvois(3,j,i)/sr2) / (sr1+sr2)
398 vec(1,j) = vec(1,j) - wface(1,j,i)
399 vec(2,j) = vec(2,j) - wface(2,j,i)
400 vec(3,j) = vec(3,j) - wface(3,j,i)
401 enddo!next J
402 DO j=1,nv46
403 alefvm_buffer%F_FACE(1,j,ii) = vec(1,j)
404 alefvm_buffer%F_FACE(2,j,ii) = vec(2,j)
405 alefvm_buffer%F_FACE(3,j,ii) = vec(3,j)
406 ENDDO
407 !---face-1
408 vx1(i) = vec(1,1)
409 vy1(i) = vec(2,1)
410 vz1(i) = vec(3,1)
411 !---face-2
412 vx2(i) = vec(1,2)
413 vy2(i) = vec(2,2)
414 vz2(i) = vec(3,2)
415 !---face-3
416 vx3(i) = vec(1,3)
417 vy3(i) = vec(2,3)
418 vz3(i) = vec(3,3)
419 !---face-4
420 vx4(i) = vec(1,4)
421 vy4(i) = vec(2,4)
422 vz4(i) = vec(3,4)
423 !---face-5
424 vx5(i) = vec(1,5)
425 vy5(i) = vec(2,5)
426 vz5(i) = vec(3,5)
427 !---face-6
428 vx6(i) = vec(1,6)
429 vy6(i) = vec(2,6)
430 vz6(i) = vec(3,6)
431 enddo!next I
432
433 CASE(4)
434 !CENTERED FVM - rho.c.U average
435 DO i=1,nel
436 ii = i + nft
437 DO j=1,nv46
438 IF(dt1==zero)THEN
439 vec(1,j) = (valel(1,i) + valvois(1,j,i) )
440 . / (valel(4,i) + valvois(4,j,i) )
441 vec(2,j) = (valel(2,i) + valvois(2,j,i) )
442 . / (valel(4,i) + valvois(4,j,i) )
443 vec(3,j) = (valel(3,i) + valvois(3,j,i) )
444 . / (valel(4,i) + valvois(4,j,i) )
445 ELSE
446 vec(1,j) = (valel(1,i)*valel(5,i) + valvois(1,j,i)*valvois(5,j,i))
447 . / (valel(4,i)*valel(5,i) + valvois(4,j,i)*valvois(5,j,i))
448 vec(2,j) = (valel(2,i)*valel(5,i) + valvois(2,j,i)*valvois(5,j,i))
449 . / (valel(4,i)*valel(5,i) + valvois(4,j,i)*valvois(5,j,i))
450 vec(3,j) = (valel(3,i)*valel(5,i) + valvois(3,j,i)*valvois(5,j,i))
451 . / (valel(4,i)*valel(5,i) + valvois(4,j,i)*valvois(5,j,i))
452 ENDIF
453 vec(1,j) = vec(1,j) - wface(1,j,i)
454 vec(2,j) = vec(2,j) - wface(2,j,i)
455 vec(3,j) = vec(3,j) - wface(3,j,i)
456 enddo!next J
457 DO j=1,nv46
458 alefvm_buffer%F_FACE(1,j,ii) = vec(1,j)
459 alefvm_buffer%F_FACE(2,j,ii) = vec(2,j)
460 alefvm_buffer%F_FACE(3,j,ii) = vec(3,j)
461 ENDDO
462 !---face-1
463 vx1(i) = vec(1,1)
464 vy1(i) = vec(2,1)
465 vz1(i) = vec(3,1)
466 !---face-2
467 vx2(i) = vec(1,2)
468 vy2(i) = vec(2,2)
469 vz2(i) = vec(3,2)
470 !---face-3
471 vx3(i) = vec(1,3)
472 vy3(i) = vec(2,3)
473 vz3(i) = vec(3,3)
474 !---face-4
475 vx4(i) = vec(1,4)
476 vy4(i) = vec(2,4)
477 vz4(i) = vec(3,4)
478 !---face-5
479 vx5(i) = vec(1,5)
480 vy5(i) = vec(2,5)
481 vz5(i) = vec(3,5)
482 !---face-6
483 vx6(i) = vec(1,6)
484 vy6(i) = vec(2,6)
485 vz6(i) = vec(3,6)
486 enddo!next I
487
488 CASE(5)
489 !Godunov Acoustic for Riemann problem
490 DO i=1,nel
491 ii = i + nft
492 IF(dt1==zero)THEN
493 DO j=1,nv46
494 vec(1,j) = (valel(1,i) + valvois(1,j,i))
495 . / (valel(4,i) + valvois(4,j,i))
496 vec(2,j) = (valel(2,i) + valvois(2,j,i))
497 . / (valel(4,i) + valvois(4,j,i))
498 vec(3,j) = (valel(3,i) + valvois(3,j,i))
499 . / (valel(4,i) + valvois(4,j,i))
500 vec(1,j) = vec(1,j) - wface(1,j,i)
501 vec(2,j) = vec(2,j) - wface(2,j,i)
502 vec(3,j) = vec(3,j) - wface(3,j,i)
503 enddo!next J
504 ELSE
505 surf(1) = one/sqrt(n1x(i)*n1x(i)+n1y(i)*n1y(i)+n1z(i)*n1z(i))
506 surf(2) = one/sqrt(n2x(i)*n2x(i)+n2y(i)*n2y(i)+n2z(i)*n2z(i))
507 surf(3) = one/sqrt(n3x(i)*n3x(i)+n3y(i)*n3y(i)+n3z(i)*n3z(i))
508 surf(4) = one/sqrt(n4x(i)*n4x(i)+n4y(i)*n4y(i)+n4z(i)*n4z(i))
509 surf(5) = one/sqrt(n5x(i)*n5x(i)+n5y(i)*n5y(i)+n5z(i)*n5z(i))
510 surf(6) = one/sqrt(n6x(i)*n6x(i)+n6y(i)*n6y(i)+n6z(i)*n6z(i))
511 DO j=1,nv46
512 term2 = surf(j) * ( valel(6,i)-valvois(6,j,i) ) / (valel(4,i)*valel(5,i)+valvois(4,j,i)*valvois(5,j,i))
513 !TERM2 = TERM2 * ZERO
514 vec(1,j) = (valel(1,i)*valel(5,i) + valvois(1,j,i)*valvois(5,j,i))
515 . / (valel(4,i)*valel(5,i)+valvois(4,j,i)*valvois(5,j,i)) + term2*n(1,j,i)
516 vec(2,j) = (valel(2,i)*valel(5,i) + valvois(2,j,i)*valvois(5,j,i))
517 . / (valel(4,i)*valel(5,i)+valvois(4,j,i)*valvois(5,j,i)) + term2*n(2,j,i)
518 vec(3,j) = (valel(3,i)*valel(5,i) + valvois(3,j,i)*valvois(5,j,i))
519 . / (valel(4,i)*valel(5,i)+valvois(4,j,i)*valvois(5,j,i)) + term2*n(3,j,i)
520 vec(1,j) = vec(1,j) - wface(1,j,i)
521 vec(2,j) = vec(2,j) - wface(2,j,i)
522 vec(3,j) = vec(3,j) - wface(3,j,i)
523 enddo!next J
524 ENDIF
525 DO j=1,nv46
526 alefvm_buffer%F_FACE(1,j,ii) = vec(1,j)
527 alefvm_buffer%F_FACE(2,j,ii) = vec(2,j)
528 alefvm_buffer%F_FACE(3,j,ii) = vec(3,j)
529 ENDDO
530 !---face-1
531 vx1(i) = vec(1,1)
532 vy1(i) = vec(2,1)
533 vz1(i) = vec(3,1)
534 !---face-2
535 vx2(i) = vec(1,2)
536 vy2(i) = vec(2,2)
537 vz2(i) = vec(3,2)
538 !---face-3
539 vx3(i) = vec(1,3)
540 vy3(i) = vec(2,3)
541 vz3(i) = vec(3,3)
542 !---face-4
543 vx4(i) = vec(1,4)
544 vy4(i) = vec(2,4)
545 vz4(i) = vec(3,4)
546 !---face-5
547 vx5(i) = vec(1,5)
548 vy5(i) = vec(2,5)
549 vz5(i) = vec(3,5)
550 !---face-6
551 vx6(i) = vec(1,6)
552 vy6(i) = vec(2,6)
553 vz6(i) = vec(3,6)
554 enddo!next I
555
556 CASE(6)
557 !INTERPOLATED
558 !WARNING : ADD SYMMETRY LATER IF FORMULATION IS KEPT (mean value for identical value on both sides)
559 DO i=1,nel
560 ii = i + nft
561 nc(1,i) = ixs(2,ii)
562 nc(2,i) = ixs(3,ii)
563 nc(3,i) = ixs(4,ii)
564 nc(4,i) = ixs(5,ii)
565 nc(5,i) = ixs(6,ii)
566 nc(6,i) = ixs(7,ii)
567 nc(7,i) = ixs(8,ii)
568 nc(8,i) = ixs(9,ii)
569 z(1) = one_over_8*sum(x(1,nc(1:8,i)))
570 z(2) = one_over_8*sum(x(2,nc(1:8,i)))
571 z(3) = one_over_8*sum(x(3,nc(1:8,i)))
572 iad2 = ale_connect%ee_connect%iad_connect(ii)
573 lgth = ale_connect%ee_connect%iad_connect(ii + 1) - iad2
574 DO j=1,lgth
575 idv = ale_connect%ee_connect%connected(iad2 + j - 1)
576 IF(idv<=0)THEN
577 vec(1:3,j) = zero
578 cycle
579 ENDIF
580 ix1=ixs(icf(1,j)+1,ii)
581 ix2=ixs(icf(2,j)+1,ii)
582 ix3=ixs(icf(3,j)+1,ii)
583 ix4=ixs(icf(4,j)+1,ii)
584 cf(1,i) = fourth*(x(1,ix1)+x(1,ix2)+x(1,ix3)+x(1,ix4))
585 cf(2,i) = fourth*(x(2,ix1)+x(2,ix2)+x(2,ix3)+x(2,ix4))
586 cf(3,i) = fourth*(x(3,ix1)+x(3,ix2)+x(3,ix3)+x(3,ix4))
587 nc(1,i)=ixs(2,idv)
588 nc(2,i)=ixs(3,idv)
589 nc(3,i)=ixs(4,idv)
590 nc(4,i)=ixs(5,idv)
591 nc(5,i)=ixs(6,idv)
592 nc(6,i)=ixs(7,idv)
593 nc(7,i)=ixs(8,idv)
594 nc(8,i)=ixs(9,idv)
595 zadj(1) = one_over_8*sum(x(1,nc(1:8,i)))
596 zadj(2) = one_over_8*sum(x(2,nc(1:8,i)))
597 zadj(3) = one_over_8*sum(x(3,nc(1:8,i)))
598 zzadj(1) = zadj(1)-z(1)
599 zzadj(2) = zadj(2)-z(2)
600 zzadj(3) = zadj(3)-z(3)
601 zcf(1) = cf(1,i) - z(1)
602 zcf(2) = cf(2,i) - z(2)
603 zcf(3) = cf(3,i) - z(3)
604 ps = zcf(1)*zzadj(1) + zcf(2)*zzadj(2) + zcf(3)*zzadj(3)
605 zcf_ = zcf(1)**2 + zcf(2)**2 + zcf(3)**2
606 zzadj_ = zzadj(1)**2 + zzadj(2)**2 + zzadj(3)**2
607 denom1 = sqrt(zcf_)
608 denom2 = sqrt(zzadj_)
609 denom = denom1*denom2
610 cos = ps/denom
611 xsi = sqrt(zcf_) * cos
612 lambda = xsi / denom2
613 !interpoler VALEL -> VALVOIS utilisant lambda
614 vec(1,j) = valel(1,i) + lambda*(valvois(1,j,i)/valvois(4,j,i)-valel(1,i)/valel(4,i)) - wface(1,j,i)
615 vec(2,j) = valel(2,i) + lambda*(valvois(2,j,i)/valvois(4,j,i)-valel(2,i)/valel(4,i)) - wface(2,j,i)
616 vec(3,j) = valel(3,i) + lambda*(valvois(3,j,i)/valvois(4,j,i)-valel(3,i)/valel(4,i)) - wface(3,j,i)
617 enddo!next J
618 DO j=1,nv46
619 alefvm_buffer%F_FACE(1,j,ii) = vec(1,j)
620 alefvm_buffer%F_FACE(2,j,ii) = vec(2,j)
621 alefvm_buffer%F_FACE(3,j,ii) = vec(3,j)
622 ENDDO
623 !---face-1
624 vx1(i) = vec(1,1)
625 vy1(i) = vec(2,1)
626 vz1(i) = vec(3,1)
627 !---face-2
628 vx2(i) = vec(1,2)
629 vy2(i) = vec(2,2)
630 vz2(i) = vec(3,2)
631 !---face-3
632 vx3(i) = vec(1,3)
633 vy3(i) = vec(2,3)
634 vz3(i) = vec(3,3)
635 !---face-4
636 vx4(i) = vec(1,4)
637 vy4(i) = vec(2,4)
638 vz4(i) = vec(3,4)
639 !---face-5
640 vx5(i) = vec(1,5)
641 vy5(i) = vec(2,5)
642 vz5(i) = vec(3,5)
643 !---face-6
644 vx6(i) = vec(1,6)
645 vy6(i) = vec(2,6)
646 vz6(i) = vec(3,6)
647 enddo!next I
648
649 END SELECT !I_ALE_SOLVER
650
651 !======================================================!
652 ! FLUXES CALCULATION ON EACH FACE !
653 ! FLUX_face = [0.5*V_face] . [2S*n] !
654 !======================================================!
655 DO i=1,nel
656 flux1(i)=half*(vx1(i)*n1x(i)+vy1(i)*n1y(i)+vz1(i)*n1z(i))
657 flux2(i)=half*(vx2(i)*n2x(i)+vy2(i)*n2y(i)+vz2(i)*n2z(i))
658 flux3(i)=half*(vx3(i)*n3x(i)+vy3(i)*n3y(i)+vz3(i)*n3z(i))
659 flux4(i)=half*(vx4(i)*n4x(i)+vy4(i)*n4y(i)+vz4(i)*n4z(i))
660 flux5(i)=half*(vx5(i)*n5x(i)+vy5(i)*n5y(i)+vz5(i)*n5z(i))
661 flux6(i)=half*(vx6(i)*n6x(i)+vy6(i)*n6y(i)+vz6(i)*n6z(i))
662 ENDDO
663
664 !DEBUG-OUTPUT---------------!
665 if(alefvm_param%IOUTP_FLUX /= 0)then
666 debug_outp = .false.
667 if(alefvm_param%IOUTP_FLUX>0)then
668 do i=lft,llt
669 ii = nft + i
670 if(ixs(11,ii)==alefvm_param%IOUTP_FLUX)THEN
671 debug_outp = .true.
672 idbf = i
673 idbl = i
674 EXIT
675 endif
676 enddo
677 elseif(alefvm_param%IOUTP_FLUX==-1)then
678 debug_outp=.true.
679 idbf = lft
680 idbl = llt
681 endif
682!#!include "lockon.inc"
683 if(debug_outp)then
684!#!include "lockon.inc"
685 print *, " |----alefvm_aflux3.F-----|"
686 print *, " | THREAD INFORMATION |"
687 print *, " |------------------------|"
688 print *, " NCYCLE =", ncycle
689 do i=idbf,idbl
690 ii = nft + i
691 print *, " brique=", ixs(11,nft+i)
692 write (*,fmt='(A,6E26.14)') " Flux(1:6) =", flux1(i),flux2(i),flux3(i),flux4(i),flux5(i),flux6(i)
693 write(*,fmt='(A24,1A20)') " ",
694 . "#--------- cell-----------"
695 write (*,fmt='(A,1E26.14)') " V_cell-X =", alefvm_buffer%FCELL(1,ii)
696 write (*,fmt='(A,1E26.14)') " V_cell-Y =", alefvm_buffer%FCELL(2,ii)
697 write (*,fmt='(A,1E26.14)') " V_cell-Z =", alefvm_buffer%FCELL(3,ii)
698 write(*,fmt='(A24,6A26)') " ",
699 . "#--------- adj_1 ---------","#--------- adj_2 ---------",
700 . "#--------- adj_3 ---------","#--------- adj_4 ---------",
701 . "#--------- adj_5 ---------","#--------- adj_6 --------#"
702 write (*,fmt='(A,6E26.14)') " V_faces-X =", alefvm_buffer%F_FACE(1,1:6,ii)
703 write (*,fmt='(A,6E26.14)') " V_faces-Y =", alefvm_buffer%F_FACE(2,1:6,ii)
704 write (*,fmt='(A,6E26.14)') " V_faces-Z =", alefvm_buffer%F_FACE(3,1:6,ii)
705 print *, " "
706 enddo
707!#!include "lockoff.inc"
708 endif
709 endif
710 !-----------------------------------------!
711
712
713 !======================================================!
714 ! TRIMATERIAL CASE INITIALIZATION (LAW51) !
715 ! -->RETURN !
716 !======================================================!
717 IF(nint(pm(19,mat(1)))==51)THEN
718 DO i=1,nel
719 flux(1,i)=flux1(i)
720 flux(2,i)=flux2(i)
721 flux(3,i)=flux3(i)
722 flux(4,i)=flux4(i)
723 flux(5,i)=flux5(i)
724 flux(6,i)=flux6(i)
725 ENDDO !next I
726 RETURN
727 ENDIF
728
729 !======================================================!
730 ! BOUNDARY FACE : no volume flux by default !
731 ! slip wall bc !
732 !======================================================!
733 DO j=1,6
734 DO i=1,nel
735 upwl(j,i)=pm(16,mat(i))
736 ENDDO !next I
737 ENDDO !next J
738
739 DO i=1,nel
740 reduc=pm(92,mat(i))
741 iad2 = ale_connect%ee_connect%iad_connect(i + nft)
742 !---face1---!
743 ii=ale_connect%ee_connect%connected(iad2 + 1 - 1)
744 IF(ii==0)THEN
745 flux1(i)=flux1(i)*reduc
746 ENDIF
747 !---face2---!
748 ii=ale_connect%ee_connect%connected(iad2 + 2 - 1)
749 IF(ii==0)THEN
750 flux2(i)=flux2(i)*reduc
751 ENDIF
752 !---face3---!
753 ii=ale_connect%ee_connect%connected(iad2 + 3 - 1)
754 IF(ii==0)THEN
755 flux3(i)=flux3(i)*reduc
756 ENDIF
757 !---face4---!
758 ii=ale_connect%ee_connect%connected(iad2 + 4 - 1)
759 IF(ii==0)THEN
760 flux4(i)=flux4(i)*reduc
761 ENDIF
762 !---face5---!
763 ii=ale_connect%ee_connect%connected(iad2 + 5 - 1)
764 IF(ii==0)THEN
765 flux5(i)=flux5(i)*reduc
766 ENDIF
767 !---face6---!
768 ii=ale_connect%ee_connect%connected(iad2 + 6 - 1)
769 IF(ii==0)THEN
770 flux6(i)=flux6(i)*reduc
771 ENDIF
772 ENDDO !next I
773
774 DO i=1,nel
775 flux(1,i)=flux1(i)-upwl(1,i)*abs(flux1(i))
776 flux(2,i)=flux2(i)-upwl(2,i)*abs(flux2(i))
777 flux(3,i)=flux3(i)-upwl(3,i)*abs(flux3(i))
778 flux(4,i)=flux4(i)-upwl(4,i)*abs(flux4(i))
779 flux(5,i)=flux5(i)-upwl(5,i)*abs(flux5(i))
780 flux(6,i)=flux6(i)-upwl(6,i)*abs(flux6(i))
781
782 flu1(i) =flux1(i)+upwl(1,i)*abs(flux1(i))
783 . +flux2(i)+upwl(2,i)*abs(flux2(i))
784 . +flux3(i)+upwl(3,i)*abs(flux3(i))
785 . +flux4(i)+upwl(4,i)*abs(flux4(i))
786 . +flux5(i)+upwl(5,i)*abs(flux5(i))
787 . +flux6(i)+upwl(6,i)*abs(flux6(i))
788 ENDDO !next I
789
790
791
792
793
794 RETURN
795 END
alefvm_aflux3
subroutine alefvm_aflux3(pm, ixs, w, flux, flu1, ale_connect, ipm, nv46, x, nel)
Definition alefvm_aflux3.F:36
ale_connectivity_mod
Definition ale_connectivity_mod.F:223
alefvm_mod
Definition alefvm_mod.F:88
alefvm_mod::alefvm_buffer
type(alefvm_buffer_), target alefvm_buffer
Definition alefvm_mod.F:120
alefvm_mod::alefvm_param
type(alefvm_param_), target alefvm_param
Definition alefvm_mod.F:121
i22tri_mod
Definition cut-cell-search_mod.F:212