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 35 of file alefvm_sfint3_int22.F.

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