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

Go to the source code of this file.

Functions/Subroutines

subroutine ns_fvm_diffusion (ale_connect, multi_fvm, time_step, ebcs_tab, diffusion, ipm, pm, iparg, elbuf_tab, nercvois, nesdvois, lercvois, lesdvois, ixs, func_value)

Function/Subroutine Documentation

◆ ns_fvm_diffusion()

subroutine ns_fvm_diffusion ( type(t_ale_connectivity), intent(in) ale_connect,
type(multi_fvm_struct), intent(inout) multi_fvm,
intent(in) time_step,
type(t_ebcs_tab), intent(inout), target ebcs_tab,
type(t_diffusion), intent(inout) diffusion,
integer, dimension(npropmi, *), intent(in) ipm,
dimension(npropm, *), intent(in) pm,
integer, dimension(nparg, *), intent(in) iparg,
type(elbuf_struct_), dimension(ngroup), intent(inout), target elbuf_tab,
integer, dimension(*), intent(in) nercvois,
integer, dimension(*), intent(in) nesdvois,
integer, dimension(*), intent(in) lercvois,
integer, dimension(*), intent(in) lesdvois,
integer, dimension(nixs, numels), intent(in) ixs,
dimension(*), intent(in) func_value )

Definition at line 36 of file ns_fvm_diffusion.F.

39 USE ale_connectivity_mod
40 USE multi_fvm_mod
41 USE diffusion_mod
42 USE ebcs_mod
43 USE elbufdef_mod
44c-----------------------------------------------
45c i m p l i c i t t y p e s
46c-----------------------------------------------
47#include "implicit_f.inc"
48! nspmd, ngroup
49#include "com01_c.inc"
50! numels
51#include "com04_c.inc"
52! nsvois
53! npropmi, npropm, nparg
54#include "param_c.inc"
55c-----------------------------------------------
56c m e s s a g e p a s s i n g
57c-----------------------------------------------
58#ifdef mpi
59#endif
60! dummy
61 type(t_ale_connectivity), intent(in) :: ale_connect
62 type(multi_fvm_struct), intent(inout) :: multi_fvm
63 my_real, intent(in) :: time_step
64 type(t_diffusion), intent(inout) :: diffusion
65 type(t_ebcs_tab), intent(inout), target :: ebcs_tab
66 integer, intent(in) :: ipm(npropmi, *)
67 my_real, intent(in) :: pm(npropm, *)
68 integer, intent(in) :: iparg(nparg, *)
69 integer, intent(in) :: nercvois(*), nesdvois(*), lercvois(*), lesdvois(*)
70 integer, intent(in) :: ixs(nixs, numels)
71 my_real, intent(in) :: func_value(*)
72 type(elbuf_struct_), target, dimension(ngroup), intent(inout) :: elbuf_tab
73c----------------------------------------------
74c l o c a l v a r i a b l e s
75c-----------------------------------------------
76 integer :: i, ii, jj, kk, ng, iebcs, typ, nelem, ielem
77 my_real :: mass, surf, ux, uy, uz, nu, dist, xk(3), xl(3), xf(3)
78 double precision, dimension(:), pointer :: sol
79 my_real :: vel(3), nx, ny, nz
80 integer :: mtn, nel, nft, imat, submatlaw, submatid, lencom, ifunc
81 type(g_bufel_), pointer :: gbuf
82 integer :: mat_nz, glob_dim, ierr, max_id, icount1, icount2
83 class (t_ebcs), pointer :: ebcs
84 integer :: iad, lgth
85 integer :: kface
86
87! diffusion coefficient
88 diffusion%nu(:) = zero
89 if (multi_fvm%nbmat == 1) then
90 do ng = 1, ngroup
91 mtn = iparg(1, ng)
92 if (mtn == 151) then
93 nel = iparg(2, ng)
94 nft = iparg(3, ng)
95 submatid = ipm(20 + 1, ixs(1, 1 + nft))
96 submatlaw = ipm(2, submatid)
97 if (submatlaw == 6) then
98 do ii = 1, nel
99 i = ii + nft
100 diffusion%nu(i) = pm(24, submatid) * multi_fvm%rho(i)
101 enddo
102 endif
103 endif
104 enddo
105 else
106 do ng = 1, ngroup
107 mtn = iparg(1, ng)
108 if (mtn == 151) then
109 nel = iparg(2, ng)
110 nft = iparg(3, ng)
111 do imat = 1, multi_fvm%nbmat
112 submatid = ipm(20 + imat, ixs(1, 1 + nft))
113 submatlaw = ipm(2, submatid)
114 if (submatlaw == 6) then
115 do ii = 1, nel
116 i = ii + nft
117 diffusion%nu(i) = diffusion%nu(i) + multi_fvm%phase_alpha(imat, i) *
118 . pm(24, submatid) * multi_fvm%phase_rho(imat, i)
119 enddo
120 endif
121 enddo
122 endif
123 enddo
124 endif
125! mpi comm
126 if (nspmd > 1) then
127 lencom = nercvois(nspmd + 1) + nesdvois(nspmd + 1)
128 call spmd_e1vois(diffusion%nu,
129 . nercvois, nesdvois, lercvois, lesdvois, lencom)
130 endif
131
132
133! ebcs / fluxout and ebcs / inlet : newman boundary conditions
134 if (ebcs_tab%nebcs_fvm > 0) then
135! flag_outlet = 0 for standard dirichlet condition: u_boundary = 0 un Euler, or w_grid in ALE. Rhs needs be modified
136! flag_outlet = -1 for newman boundary condition
137! flag_outlet = iebcs for dirichlet boundary condition, value of velocity is then taken as the one given in the inlet card.
138 if (.not. diffusion%outlet_flagged) then
139 do iebcs = 1, ebcs_tab%nebcs_fvm
140 typ = ebcs_tab%tab(iebcs)%poly%type
141 nelem = ebcs_tab%tab(iebcs)%poly%nb_elem
142 ebcs => ebcs_tab%tab(iebcs)%poly
143 select type (ebcs)
144 type is (t_ebcs_inlet)
145! inlet
146 do ielem = 1, nelem
147 ii = ebcs_tab%tab(iebcs)%poly%ielem(ielem)
148 jj = ebcs_tab%tab(iebcs)%poly%iface(ielem)
149 diffusion%flag_outlet(6 * (ii - 1) + jj) = iebcs
150 enddo
151 type is (t_ebcs_fluxout)
152! fluxout
153 do ielem = 1, nelem
154 ii = ebcs_tab%tab(iebcs)%poly%ielem(ielem)
155 jj = ebcs_tab%tab(iebcs)%poly%iface(ielem)
156 diffusion%flag_outlet(6 * (ii - 1) + jj) = -1
157 enddo
158 end select
159 enddo
160 diffusion%outlet_flagged = .true.
161 endif
162 endif
163
164 diffusion%rhs%val(1:3 * numels) = zero
165! fill in matrix + rhs
166 icount1 = 0
167 icount2 = numels
168 do ii = 1, numels
169 iad = ale_connect%ee_connect%iad_connect(ii)
170 lgth = ale_connect%ee_connect%iad_connect(ii+1) - iad
171 icount1 = icount1 + 1
172 diffusion%mat%irow(icount1) = ale_connect%idglob%id(ii)
173 diffusion%mat%jcol(icount1) = ale_connect%idglob%id(ii)
174 xk(1:3) = multi_fvm%elem_data%centroid(1:3, ii)
175 mass = multi_fvm%vol(ii) * multi_fvm%rho(ii)
176 diffusion%mat%val(icount1) = mass
177! diagonal part
178 do jj = 1, lgth
179 kk = ale_connect%ee_connect%connected(iad + jj - 1)
180 if (kk > 0) then
181 surf = multi_fvm%face_data%surf(jj, ii)
182 xl(1:3) = multi_fvm%elem_data%centroid(1:3, kk)
183 dist = sqrt((xk(1) - xl(1))**2 + (xk(2) - xl(2))**2 + (xk(3) - xl(3))**2)
184 nu = 0.5 * (diffusion%nu(ii) + diffusion%nu(kk))
185 diffusion%mat%val(icount1) = diffusion%mat%val(icount1) + time_step * surf * nu / dist
186 else if (diffusion%flag_outlet(6 * (ii - 1) + jj) == -1) then
187! Neumann boundary condition
188 else if (diffusion%flag_outlet(6 * (ii - 1) + jj) >= 0) then
189! Dirichlet boundary condition
190 if (diffusion%flag_outlet(6 * (ii - 1) + jj) == 0) then
191 vel(1:3) = zero
192 else
193 iebcs = diffusion%flag_outlet(6 * (ii - 1) + jj)
194 ebcs => ebcs_tab%tab(iebcs)%poly
195 nx = multi_fvm%face_data%normal(1, jj, ii)
196 ny = multi_fvm%face_data%normal(2, jj, ii)
197 nz = multi_fvm%face_data%normal(3, jj, ii)
198 select type (ebcs)
199 type is (t_ebcs_inlet)
200 if (ebcs%fvm_inlet_data%vector_velocity == 0) then
201! normal velocity imposed
202 ifunc = ebcs%fvm_inlet_data%func_vel(1)
203 if (ifunc > 0) then
204 vel(1) = -ebcs%fvm_inlet_data%val_vel(1) * func_value(ifunc) * nx
205 vel(2) = -ebcs%fvm_inlet_data%val_vel(2) * func_value(ifunc) * ny
206 vel(3) = -ebcs%fvm_inlet_data%val_vel(3) * func_value(ifunc) * nz
207 else
208 vel(1) = -ebcs%fvm_inlet_data%val_vel(1) * nx
209 vel(2) = -ebcs%fvm_inlet_data%val_vel(2) * ny
210 vel(3) = -ebcs%fvm_inlet_data%val_vel(3) * nz
211 endif
212 else
213! all components of velocity are imposed
214 ifunc = ebcs%fvm_inlet_data%func_vel(1)
215 if (ifunc > 0) then
216 vel(1) = ebcs%fvm_inlet_data%val_vel(1) * func_value(ifunc)
217 else
218 vel(1) = ebcs%fvm_inlet_data%val_vel(1)
219 endif
220 ifunc = ebcs%fvm_inlet_data%func_vel(2)
221 if (ifunc > 0) then
222 vel(2) = ebcs%fvm_inlet_data%val_vel(2) * func_value(ifunc)
223 else
224 vel(2) = ebcs%fvm_inlet_data%val_vel(2)
225 endif
226 ifunc = ebcs%fvm_inlet_data%func_vel(3)
227 if (ifunc > 0) then
228 vel(3) = ebcs%fvm_inlet_data%val_vel(3) * func_value(ifunc)
229 else
230 vel(3) = ebcs%fvm_inlet_data%val_vel(3)
231 endif
232 endif
233 end select
234 endif
235 surf = multi_fvm%face_data%surf(jj, ii)
236 xf(1:3) = multi_fvm%face_data%centroid(1:3, jj, ii)
237 xl(1:3) = two * xf(1:3) - xk(1:3)
238 dist = sqrt((xk(1) - xl(1))**2 + (xk(2) - xl(2))**2 + (xk(3) - xl(3))**2)
239 nu = diffusion%nu(ii)
240 diffusion%mat%val(icount1) = diffusion%mat%val(icount1) + time_step * surf * nu / dist
241 diffusion%rhs%val(ii + (1 - 1) * numels) = diffusion%rhs%val(ii + (1 - 1) * numels) +
242 . time_step * surf * nu * vel(1) / dist
243 diffusion%rhs%val(ii + (2 - 1) * numels) = diffusion%rhs%val(ii + (2 - 1) * numels) +
244 . time_step * surf * nu * vel(2) / dist
245 diffusion%rhs%val(ii + (3 - 1) * numels) = diffusion%rhs%val(ii + (3 - 1) * numels) +
246 . time_step * surf * nu * vel(3) / dist
247 endif
248 enddo
249! non diagonal part
250 do jj = 1, lgth
251 kk = ale_connect%ee_connect%connected(iad + jj - 1)
252 if (kk > 0) then
253 icount2 = icount2 + 1
254 diffusion%mat%irow(icount2) = ale_connect%idglob%id(ii)
255 diffusion%mat%jcol(icount2) = ale_connect%idglob%id(kk)
256 surf = multi_fvm%face_data%surf(jj, ii)
257 xl(1:3) = multi_fvm%elem_data%centroid(1:3, kk)
258 dist = sqrt((xk(1) - xl(1))**2 + (xk(2) - xl(2))**2 + (xk(3) - xl(3))**2)
259 nu = 0.5 * (diffusion%nu(ii) + diffusion%nu(kk))
260 diffusion%mat%val(icount2) = - time_step * surf * nu / dist
261 endif
262 enddo
263 enddo
264
265! right hand side
266 icount1 = 0
267 do ii = 1, numels
268 ux = multi_fvm%vel(1, ii)
269 uy = multi_fvm%vel(2, ii)
270 uz = multi_fvm%vel(3, ii)
271 mass = multi_fvm%vol(ii) * multi_fvm%rho(ii)
272 icount1 = icount1 + 1
273 diffusion%rhs%irow(icount1) = ale_connect%idglob%id(ii)
274 diffusion%rhs%val(ii + (1 - 1) * numels) = diffusion%rhs%val(ii + (1 - 1) * numels) + mass * ux
275 diffusion%rhs%val(ii + (2 - 1) * numels) = diffusion%rhs%val(ii + (2 - 1) * numels) + mass * uy
276 diffusion%rhs%val(ii + (3 - 1) * numels) = diffusion%rhs%val(ii + (3 - 1) * numels) + mass * uz
277 enddo
278
279 nullify(sol)
280 call diffusion%solve_diffusion()
281 call diffusion%get_solution(sol, glob_dim)
282
283 if (associated(sol)) then
284 do ii = 1, numels
285 ux = multi_fvm%vel(1, ii)
286 uy = multi_fvm%vel(2, ii)
287 uz = multi_fvm%vel(3, ii)
288 multi_fvm%eint(ii) = multi_fvm%eint(ii) + 0.5 * multi_fvm%rho(ii) * (ux**2 + uy**2 + uz**2)
289 ux = diffusion%sol(ale_connect%idglob%id(ii) + 0 * glob_dim)
290 uy = diffusion%sol(ale_connect%idglob%id(ii) + 1 * glob_dim)
291 uz = diffusion%sol(ale_connect%idglob%id(ii) + 2 * glob_dim)
292 multi_fvm%vel(1, ii) = ux
293 multi_fvm%vel(2, ii) = uy
294 multi_fvm%vel(3, ii) = uz
295 multi_fvm%eint(ii) = multi_fvm%eint(ii) - 0.5 * multi_fvm%rho(ii) * (ux**2 + uy**2 + uz**2)
296 enddo
297 endif
my_real
#define my_real
Definition cppsort.cpp:32
ale_connectivity_mod
Definition ale_connectivity_mod.F:223
diffusion_mod
Definition diffusion_mod.F:36
spmd_e1vois
subroutine spmd_e1vois(phi, nercvois, nesdvois, lercvois, lesdvois, lencom)
Definition spmd_cfd.F:375