OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
alefvm_sfint3_int22.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "com01_c.inc"
#include "com04_c.inc"
#include "inter22.inc"
#include "param_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine alefvm_sfint3_int22 (ixs, nv46, itask, nbf, nbl, nin)

Function/Subroutine Documentation

◆ alefvm_sfint3_int22()

subroutine alefvm_sfint3_int22 ( integer, dimension(nixs,numels) ixs,
integer nv46,
integer itask,
integer nbf,
integer nbl,
integer nin )

Definition at line 34 of file alefvm_sfint3_int22.F.

35C-----------------------------------------------
36C D e s c r i p t i o n
37C-----------------------------------------------
38C 'alefvm' is related to a collocated scheme (built from FVM and based on Godunov scheme)
39C which was temporarily introduced for experimental option /INTER/TYPE22 (FSI coupling with cut cell method)
40C This cut cell method is not completed, abandoned, and is not an official option.
41C There is no other use for this scheme which is automatically enabled when /INTER/TYPE22 is defined (INT22>0 => IALEFVM=1).
42C
43C This subroutine is treating an uncut cell.
44C-----------------------------------------------
45C M o d u l e s
46C-----------------------------------------------
47 USE alefvm_mod
48 USE i22bufbric_mod
49 USE i22tri_mod
50C-----------------------------------------------
51C I m p l i c i t T y p e s
52C-----------------------------------------------
53#include "implicit_f.inc"
54C-----------------------------------------------
55C G l o b a l P a r a m e t e r s
56C-----------------------------------------------
57#include "mvsiz_p.inc"
58C-----------------------------------------------
59C C o m m o n B l o c k s
60C-----------------------------------------------
61#include "com01_c.inc"
62#include "com04_c.inc"
63#include "inter22.inc"
64#include "param_c.inc"
65C-----------------------------------------------
66C D e s c r i p t i o n
67C-----------------------------------------------
68C This subroutines computes internal forces for
69C finite volume scheme (IALEFVM==1)
70C
71C If option is not detected in input file then
72C subroutine is unplugged
73C-----------------------------------------------
74C D u m m y A r g u m e n t s
75C-----------------------------------------------
76 INTEGER :: IXS(NIXS,NUMELS),NV46, ITASK
77 INTEGER :: NBF, NBL, NIN
78C-----------------------------------------------
79C L o c a l V a r i a b l e s
80C-----------------------------------------------
81 INTEGER :: IB, IBm, IBs, IADJ, MCELL, NumSECND, NADJ, ISECND
82 INTEGER :: I, II, J, IBv, NBCUT, NBCUTv, K, IV
83 INTEGER :: ICELLv, JV, ICELLs, ClosedSurface, IPRES_MOM, ISGN_NORM
84 my_real :: f0(3), surf, pf
85 my_real :: norm, z1,z2, p1,p2, un1,un2, denom
86 my_real :: n(3), area, areav, js, coef1, coef2
87 my_real :: theta, m1, m2, mf
88C-----------------------------------------------
89C P r e - C o n d i t i o n s
90C-----------------------------------------------
91 IF(alefvm_param%IEnabled==0) RETURN
92 IF(int22 == 0) RETURN
93C-----------------------------------------------
94C S o u r c e L i n e s
95C-----------------------------------------------
96
97c 0. Domaine d'integration : supercell = main u (U secnd)
98c 1. Calcul de Pf
99c 2. Calcul de Ff(1:3) = -Pf.S.n
100c 3. Assemblge de FINT(1:3) = Sum Ff(1:3)
101
102 !IPRES_MOM = 0,1,2,3,4 : (P1+P2)/2
103 !IPRES_MOM = 5 : (rho1c1P1+rho2c2P2)/(rho1c1+rho2c2) + (rho1c1*rho2c2/(rho1c1+rho2c2)) <uel-uadj,nel> acoustic solver
104
105 ipres_mom = alefvm_param%ISOLVER
106
107! print *, "IPRES_MOM, SFINT3, int22", IPRES_MOM
108
109 SELECT CASE(ipres_mom)
110 CASE (5)
111 !Godunov
112 coef1 = zero
113 coef2 = one
114 CASE DEFAULT
115 coef1 = one
116 coef2 = zero
117 END SELECT
118
119 !UN1 is <U1,N1>
120 !UN2 is <U2,N2>
121 !then (<Ucell-Uadj,ncell> = UN1+UN2
122
123 !Inner Forces
124 !------------
125 DO ib=nbf,nbl
126 !---STARTING WITH main CELL
127 !-----------------------------------
128 f0(1:3) = zero
129 ii = brick_list(nin,ib)%ID
130 mcell = brick_list(nin,ib)%MainID
131 nbcut = brick_list(nin,ib)%NBCUT
132 IF (mcell==0)cycle !no internal force to compute at this centroid
133 !init. with main cut faces
134 p1 = brick_list(nin,ib)%SIG(0)
135 z1 = brick_list(nin,ib)%RHOC
136 m1 = brick_list(nin,ib)%MACH
137
138 IF(nbcut>0)THEN
139 !---Face-0
140 !-----------------------------------
141 isgn_norm = one
142 IF(mcell==9)THEN
143 isgn_norm=-one
144 DO k=1,nbcut
145 n(1:3) = isgn_norm* brick_list(nin,ib)%PCUT(k)%N(1:3)
146 norm = sqrt(n(1)**2+n(2)**2+n(3)**2)
147 n(1) = n(1)/norm
148 n(2) = n(2)/norm
149 n(3) = n(3)/norm
150 surf = brick_list(nin,ib)%PCUT(k)%SCUT(1)
151 un1 = brick_list(nin,ib)%POLY(mcell)%FACE0%U_N(k)
152 mf = m1
153 theta = min(one,mf) !see dellacherie,omnes,raviart : fixing low mach issue
154 pf = p1 + coef2*theta*half*z1*un1
155 f0(1) = f0(1) -pf*surf*n(1)
156 f0(2) = f0(2) -pf*surf*n(2)
157 f0(3) = f0(3) -pf*surf*n(3)
158 enddo!next K
159 ELSE
160 isgn_norm = one
161 k = mcell
162 n(1:3) = isgn_norm* brick_list(nin,ib)%PCUT(k)%N(1:3)
163 norm = sqrt(n(1)**2+n(2)**2+n(3)**2)
164 n(1) = n(1)/norm
165 n(2) = n(2)/norm
166 n(3) = n(3)/norm
167 surf = brick_list(nin,ib)%PCUT(k)%SCUT(1)
168 un1 = brick_list(nin,ib)%POLY(mcell)%FACE0%U_N(k)
169 !Pf = P1
170 mf = m1
171 theta = min(one,mf) !see dellacherie,omnes,raviart : fixing low mach issue
172 pf = p1 + coef2*theta*half*z1*un1
173 f0(1) = f0(1) -pf*surf*n(1)
174 f0(2) = f0(2) -pf*surf*n(2)
175 f0(3) = f0(3) -pf*surf*n(3)
176 ENDIF
177 ENDIF
178
179 !--Face-1:6
180 !-----------------------------------
181 closedsurface = 0
182 DO j=1,nv46
183 nadj = brick_list(nin,ib)%POLY(mcell)%FACE(j)%NAdjCell
184 closedsurface = brick_list(nin,ib)%ClosedSurf(j)
185 !
186 ! +---------+---------+
187 ! + \ + +
188 ! + \+ +
189 ! + + +
190 ! + + <----------- closed surface
191 ! + + +
192 ! + /+ +
193 ! + / + +
194 ! +---------+---------+
195 !
196 !if(closedsurface == 1)then
197 ! print *, "brick_id=",ixs(11,ii)," face=",J, " has a closed surface"
198 ! print *, "internal forces treated without adjacent cell."
199 !endif
200 IF(nadj==0 .OR. closedsurface == 1)THEN
201 !sliding rigid wall boundary condition
202 un1 = brick_list(nin,ib)%POLY(mcell)%FACE(j)%U_N
203 un2 = zero
204 mf = m1
205 theta = min(one,mf) !see dellacherie,omnes,raviart : fixing low mach issue
206 !Pf = P1
207 pf = p1 + coef2*theta*half*z1*un1
208 n(1:3) = brick_list(nin,ib)%N(j,1:3)
209 surf = brick_list(nin,ib)%POLY(mcell)%FACE(j)%Surf
210 f0(1) = f0(1) -pf*surf*n(1)
211 f0(2) = f0(2) -pf*surf*n(2)
212 f0(3) = f0(3) -pf*surf*n(3)
213 ELSE
214 n(1:3) = brick_list(nin,ib)%N(j,1:3)
215 area = brick_list(nin,ib)%POLY(mcell)%FACE(j)%Surf
216 iv = brick_list(nin,ib)%Adjacent_Brick(j,1)
217 ibv = brick_list(nin,ib)%Adjacent_Brick(j,4)
218 jv = brick_list(nin,ib)%Adjacent_Brick(j,5)
219 un1 = brick_list(nin,ib)%POLY(mcell)%FACE(j)%U_N
220 IF(area<=zero)cycle
221 DO iadj=1,nadj
222 icellv = brick_list(nin,ib)%POLY(mcell)%FACE(j)%Adjacent_Cell(iadj)
223 ibm = 0
224 IF(ibv/=0)THEN
225 ibm = brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(4)
226 IF(ibm == ib)cycle ! inner face skipped
227 areav = brick_list(nin,ibv)%POLY(icellv)%FACE(jv)%Surf
228 surf = min(area,areav)
229 p2 = brick_list(nin,ibm)%SIG(0)
230 z2 = brick_list(nin,ibm)%RHOC
231 m2 = brick_list(nin,ibm)%MACH
232 un2 = brick_list(nin,ibv)%POLY(icellv)%FACE(jv)%U_N
233 denom = z1+z2
234 !Pf = COEF1*(HALF*(P1+P2)) + COEF2*((Z1*P2 + Z2*P1)/DENOM
235 mf = min(m1,m2)
236 theta = min(one,mf) !see dellacherie,omnes,raviart : fixing low mach issue
237 pf = coef1*(half*(p1+p2)) + coef2*((z1*p2 + z2*p1)/denom + theta*z1*z2*(un1+un2)/denom)
238 f0(1) = f0(1) - pf*surf*n(1)
239 f0(2) = f0(2) - pf*surf*n(2)
240 f0(3) = f0(3) - pf*surf*n(3)
241 else! THEN IBV==0
242 p2 = alefvm_buffer%F_FACE(1 ,4 ,iv)
243 z2 = alefvm_buffer%F_FACE(1 ,3 ,iv)
244 un2 = alefvm_buffer%F_FACE(3 ,jv ,iv)
245 m2 = alefvm_buffer%F_FACE(1 ,5 ,iv)
246 surf = area
247 denom = z1+z2
248 !Pf = COEF1*(HALF*(P1+P2)) + COEF2*((Z1*P2 + Z2*P1)/DENOM
249 mf = min(m1,m2)
250 theta = min(one,mf) !see dellacherie,omnes,raviart : fixing low mach issue
251 pf = coef1*(half*(p1+p2)) + coef2*((z1*p2 + z2*p1)/denom + theta*z1*z2*(un1+un2)/denom)
252 f0(1) = f0(1) - pf*surf*n(1)
253 f0(2) = f0(2) - pf*surf*n(2)
254 f0(3) = f0(3) - pf*surf*n(3)
255 ENDIF
256 enddo!next IADJ
257 ENDIF !(NADJ==0)
258
259 enddo!next J
260
261 !---GO ON WITH ATTACHED SECND CELLS
262 !-----------------------------------
263 numsecnd = brick_list(nin,ib)%SecndList%Num
264 DO isecnd=1,numsecnd
265 ibs = brick_list(nin,ib)%SecndList%IBV(isecnd)
266 icells = brick_list(nin,ib)%SecndList%ICELLv(isecnd)
267 js = brick_list(nin,ibs)%POLY(icells)%WhereIsMain(1)
268 nbcutv = brick_list(nin,ibs)%NBCUT
269 !---Face-0
270 !-----------------------------------
271
272 !non ne pas sommer sur les plans de coupes il ne sont pas forcemment tous voisins :
273 !example : inter22/1D_EXPANSION/9ELEMS/LAW51/PerfectGas/EULER/SMALL_MAIN/0.FVM/MAT51
274 !print *, "pause" ; pause
275
276 !! boucler si icells=9
277
278
279 isgn_norm = one
280 IF(icells==9)THEN
281 isgn_norm=-one
282 DO k=1,nbcutv
283 n(1:3) = isgn_norm* brick_list(nin,ibs)%PCUT(k)%N(1:3)
284 norm = sqrt(n(1)**2+n(2)**2+n(3)**2)
285 n(1) = n(1)/norm
286 n(2) = n(2)/norm
287 n(3) = n(3)/norm
288 surf = brick_list(nin,ibs)%PCUT(k)%SCUT(1)
289 un1 = brick_list(nin,ibs)%POLY(icells)%FACE0%U_N(k)
290 !Pf = P1
291 mf = m1
292 theta = min(one,mf) !see dellacherie,omnes,raviart : fixing low mach issue
293 pf = p1 + coef2*theta*half*z1*un1
294 f0(1) = f0(1) -pf*surf*n(1)
295 f0(2) = f0(2) -pf*surf*n(2)
296 f0(3) = f0(3) -pf*surf*n(3)
297 ENDDO
298 ELSE
299 k=icells
300 n(1:3) = isgn_norm* brick_list(nin,ibs)%PCUT(k)%N(1:3)
301 norm = sqrt(n(1)**2+n(2)**2+n(3)**2)
302 n(1) = n(1)/norm
303 n(2) = n(2)/norm
304 n(3) = n(3)/norm
305 surf = brick_list(nin,ibs)%PCUT(k)%SCUT(1)
306 un1 = brick_list(nin,ibs)%POLY(icells)%FACE0%U_N(k)
307 mf = m1
308 theta = min(one,mf) !see dellacherie,omnes,raviart : fixing low mach issue
309 !Pf = P1
310 pf = p1 + coef2*theta*half*z1*un1
311 f0(1) = f0(1) -pf*surf*n(1)
312 f0(2) = f0(2) -pf*surf*n(2)
313 f0(3) = f0(3) -pf*surf*n(3)
314 ENDIF
315
316 !---Face-1:6
317 !-----------------------------------
318 DO j=1,nv46
319 IF(j==js) cycle ! do no take into account internal face (inside supercell), i.e. between main cell and secnd cell
320 nadj = brick_list(nin,ibs)%POLY(icells)%FACE(j)%NAdjCell
321 z2 = zero
322 IF(nadj==0 .OR. closedsurface == 1)THEN
323 !sliding rigid wall boundary condition
324 un1 = brick_list(nin,ibs)%POLY(icells)%FACE(j)%U_N
325 !Pf = P1
326 mf = m1
327 theta = min(one,mf) !see dellacherie,omnes,raviart : fixing low mach issue
328 pf = p1 + coef2*theta*half*z1*un1
329 n(1:3) = brick_list(nin,ibs)%N(j,1:3)
330 surf = brick_list(nin,ibs)%POLY(icells)%FACE(j)%Surf
331 f0(1) = f0(1) -pf*surf*n(1)
332 f0(2) = f0(2) -pf*surf*n(2)
333 f0(3) = f0(3) -pf*surf*n(3)
334 ELSE
335 n(1:3) = brick_list(nin,ibs)%N(j,1:3)
336 area = brick_list(nin,ibs)%POLY(icells)%FACE(j)%Surf
337 iv = brick_list(nin,ibs)%Adjacent_Brick(j,1)
338 ibv = brick_list(nin,ibs)%Adjacent_Brick(j,4)
339 jv = brick_list(nin,ibs)%Adjacent_Brick(j,5)
340 un1 = brick_list(nin,ibs)%POLY(icells)%FACE(j)%U_N
341 IF(area<=zero)cycle
342 DO iadj=1,nadj
343 icellv = brick_list(nin,ibs)%POLY(icells)%FACE(j)%Adjacent_Cell(iadj)
344 ibm = 0
345 IF(ibv/=0)THEN
346 ibm = brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(4)
347 areav = brick_list(nin,ibv)%POLY(icellv)%FACE(jv)%Surf
348 surf = min(area,areav)
349 p2 = brick_list(nin,ibm)%SIG(0)
350 z2 = brick_list(nin,ibm)%RHOC
351 un2 = brick_list(nin,ibv)%POLY(icellv)%FACE(jv)%U_N
352 denom = z1+z2
353 m2 = brick_list(nin,ibm)%MACH
354 mf = min(m1,m2)
355 theta = min(one,mf) !see dellacherie,omnes,raviart : fixing low mach issue
356 ! Pf = COEF1*(HALF*(P1+P2)) + COEF2*((Z1*P2 + Z2*P1)/DENOM
357 pf = coef1*(half*(p1+p2)) + coef2*((z1*p2 + z2*p1)/denom + theta*z1*z2*(un1+un2)/denom)
358 IF(ibm == ib)cycle ! inner face skipped
359 f0(1) = f0(1) - pf*surf*n(1)
360 f0(2) = f0(2) - pf*surf*n(2)
361 f0(3) = f0(3) - pf*surf*n(3)
362 else! THEN IBV==0
363 p2 = alefvm_buffer%F_FACE(1 ,4 ,iv)
364 z2 = alefvm_buffer%F_FACE(1 ,3 ,iv)
365 un2 = alefvm_buffer%F_FACE(3 ,jv ,iv)
366 m2 = alefvm_buffer%F_FACE(1 ,5 ,iv)
367 surf = area
368 denom = z1+z2
369 !Pf = COEF1*(HALF*(P1+P2)) + COEF2*((Z1*P2 + Z2*P1)/DENOM
370 mf = min(m1,m2)
371 theta = min(one,mf) !see dellacherie,omnes,raviart : fixing low mach issue
372 pf = coef1*(half*(p1+p2)) + coef2*((z1*p2 + z2*p1)/denom + theta*z1*z2*(un1+un2)/denom)
373 f0(1) = f0(1) - pf*surf*n(1)
374 f0(2) = f0(2) - pf*surf*n(2)
375 f0(3) = f0(3) - pf*surf*n(3)
376 ENDIF
377 enddo!next IADJ
378 ENDIF !(NADJ==0)
379 enddo!next J
380 enddo! next ISECND
381
382
383 alefvm_buffer%FINT_CELL(1,ii) = f0(1)
384 alefvm_buffer%FINT_CELL(2,ii) = f0(2)
385 alefvm_buffer%FINT_CELL(3,ii) = f0(3)
386
387 ! write (*,FMT='(A,I6,A,3F30.16)') "22.brick ID=", ixs(11,II)," Fint=", F0(1:3)
388
389
390 !Fcell already contains gravity force, so do not erase but add
391 alefvm_buffer%FCELL(1,ii) = alefvm_buffer%FCELL(1,ii) + f0(1) !+ FEXT_CELL(1,II)
392 alefvm_buffer%FCELL(2,ii) = alefvm_buffer%FCELL(2,ii) + f0(2) !+ FEXT_CELL(2,II)
393 alefvm_buffer%FCELL(3,ii) = alefvm_buffer%FCELL(3,ii) + f0(3) !+ FEXT_CELL(3,II)
394
395 enddo!next IB
396
397 !-------------------------------------------------------------!
398 ! Debug Output !
399 !-------------------------------------------------------------!
400 if(alefvm_param%IOUTP_FINT /= 0)then
401
402 call my_barrier
403
404 if(itask==0)then
405 print *, " |---alefvm_sfint3_int22.F---|"
406 print *, " | THREAD INFORMATION |"
407 print *, " |---------------------------|"
408 print *, " NCYCLE =", ncycle
409
410 do i=1,nb
411 ii = brick_list(nin,i)%id
412 mcell = brick_list(nin,i)%MainID
413 IF (mcell==0)cycle
414 print *, " brique_=", ixs(11,ii)
415 write(*,fmt='(A34,6A26)') " ",
416 . "#--- internal force -----#"
417
418
419 write (*,fmt='(A,8E26.14)' ) " force-X =", alefvm_buffer%FCELL(1,ii)
420 write (*,fmt='(A,8E26.14)' ) " force-Y =", alefvm_buffer%FCELL(2,ii)
421 write (*,fmt='(A,8E26.14)' ) " force-Z =", alefvm_buffer%FCELL(3,ii)
422 write (*,fmt='(A,8E26.14)' ) " P =", brick_list(nin,i)%SIG(0)
423
424
425 print *, " "
426 enddo
427
428 endif!if(itask)
429
430 call my_barrier
431
432 endif!if(IALEFVM_OUTP_FINT /= 0)
433 !-----------------------------------------!
434
435 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
min
#define min(a, b)
Definition macros.h:20
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
i22bufbric_mod
Definition cut-cell-search_mod.F:95
i22bufbric_mod::brick_list
type(brick_entity), dimension(:,:), allocatable, target brick_list
Definition cut-cell-search_mod.F:105
i22tri_mod
Definition cut-cell-search_mod.F:212
i22tri_mod::nb
integer nb
Definition cut-cell-search_mod.F:487
my_barrier
subroutine my_barrier
Definition machine.F:31