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 37 of file ns_fvm_diffusion.F.

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