Skip to content

Commit bce7426

Browse files
Dimitrios AdamDimitrios Adam
Dimitrios Adam
authored and
Dimitrios Adam
committed
Temperature stuff
1 parent 4748026 commit bce7426

6 files changed

Lines changed: 255 additions & 38 deletions

File tree

‎examples/1D_multispecies_diffusion/case.py‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -29,9 +29,9 @@
2929
"p": 0,
3030
"dt": float(dt),
3131
"t_step_start": 0,
32-
"t_step_stop": NT,
33-
"t_step_save": NS,
34-
"t_step_print": NS,
32+
"t_step_stop": 1,
33+
"t_step_save": 1,
34+
"t_step_print": 1,
3535
"parallel_io": "F",
3636
"model_eqns": 2,
3737
"num_fluids": 1,
@@ -47,8 +47,8 @@
4747
"riemann_solver": 2,
4848
"wave_speeds": 2,
4949
"avg_state": 1,
50-
"bc_x%beg": -1,
51-
"bc_x%end": -1,
50+
"bc_x%beg": -6,
51+
"bc_x%end": -6,
5252
"viscous": "F",
5353
"chemistry": "T",
5454
"chem_params%diffusion": "T",

‎examples/1D_reactive_shocktube/case.py‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -38,7 +38,7 @@
3838
u_r = -487.34
3939

4040
L = 0.12
41-
Nx = 400 * args.scale
41+
Nx = 100 * args.scale
4242
dx = L / Nx
4343
dt = dx / abs(u_r) * 0.02
4444
Tend = 230e-6
@@ -58,9 +58,9 @@
5858
"p": 0,
5959
"dt": float(dt),
6060
"t_step_start": 0,
61-
"t_step_stop": NT,
62-
"t_step_save": NS,
63-
"t_step_print": NS,
61+
"t_step_stop": 1,
62+
"t_step_save": 1,
63+
"t_step_print": 1,
6464
"parallel_io": "F" if args.mfc.get("mpi", True) else "F",
6565
# Simulation Algorithm Parameters
6666
"model_eqns": 2,
@@ -81,7 +81,7 @@
8181
"bc_x%end": -3,
8282
# Chemistry
8383
"chemistry": "F" if not args.chemistry else "T",
84-
"chem_params%diffusion": "F",
84+
"chem_params%diffusion": "T",
8585
"chem_params%reactions": "T",
8686
"cantera_file": ctfile,
8787
# Formatted Database Files Structure Parameters

‎src/common/m_boundary_common.fpp‎

Lines changed: 111 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,10 @@ module m_boundary_common
2121

2222
use m_compile_specific
2323

24+
use m_thermochem, only: &
25+
num_species, get_mixture_molecular_weight, gas_constant
26+
27+
2428
implicit none
2529

2630
type(scalar_field), dimension(:, :), allocatable :: bc_buffers
@@ -86,28 +90,29 @@ contains
8690
!> The purpose of this procedure is to populate the buffers
8791
!! of the primitive variables, depending on the selected
8892
!! boundary conditions.
89-
impure subroutine s_populate_variables_buffers(bc_type, q_prim_vf, pb_in, mv_in)
93+
impure subroutine s_populate_variables_buffers(bc_type, q_prim_vf, pb_in, mv_in, q_T_sf)
9094

9195
type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
96+
type(scalar_field), optional, intent(inout) :: q_T_sf
9297
real(stp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb_in, mv_in
9398
type(integer_field), dimension(1:num_dims, 1:2), intent(in) :: bc_type
9499

95100
integer :: k, l
96101

97102
! Population of Buffers in x-direction
98103
if (bc_x%beg >= 0) then
99-
call s_mpi_sendrecv_variables_buffers(q_prim_vf, 1, -1, sys_size, pb_in, mv_in)
104+
call s_mpi_sendrecv_variables_buffers(q_prim_vf, 1, -1, sys_size, pb_in, mv_in, q_T_sf = q_T_sf)
100105
else
101106
$:GPU_PARALLEL_LOOP(private='[l,k]', collapse=2)
102107
do l = 0, p
103108
do k = 0, n
104109
select case (int(bc_type(1, 1)%sf(0, k, l)))
105110
case (BC_CHAR_SUP_OUTFLOW:BC_GHOST_EXTRAP)
106-
call s_ghost_cell_extrapolation(q_prim_vf, 1, -1, k, l)
111+
call s_ghost_cell_extrapolation(q_prim_vf, 1, -1, k, l, q_T_sf)
107112
case (BC_REFLECTIVE)
108113
call s_symmetry(q_prim_vf, 1, -1, k, l, pb_in, mv_in)
109114
case (BC_PERIODIC)
110-
call s_periodic(q_prim_vf, 1, -1, k, l, pb_in, mv_in)
115+
call s_periodic(q_prim_vf, 1, -1, k, l, pb_in, mv_in, q_T_sf)
111116
case (BC_SLIP_WALL)
112117
call s_slip_wall(q_prim_vf, 1, -1, k, l)
113118
case (BC_NO_SLIP_WALL)
@@ -126,18 +131,18 @@ contains
126131
end if
127132

128133
if (bc_x%end >= 0) then
129-
call s_mpi_sendrecv_variables_buffers(q_prim_vf, 1, 1, sys_size, pb_in, mv_in)
134+
call s_mpi_sendrecv_variables_buffers(q_prim_vf, 1, 1, sys_size, pb_in, mv_in, q_T_sf = q_T_sf)
130135
else
131136
$:GPU_PARALLEL_LOOP(private='[l,k]', collapse=2)
132137
do l = 0, p
133138
do k = 0, n
134139
select case (int(bc_type(1, 2)%sf(0, k, l)))
135140
case (BC_CHAR_SUP_OUTFLOW:BC_GHOST_EXTRAP) ! Ghost-cell extrap. BC at end
136-
call s_ghost_cell_extrapolation(q_prim_vf, 1, 1, k, l)
141+
call s_ghost_cell_extrapolation(q_prim_vf, 1, 1, k, l, q_T_sf)
137142
case (BC_REFLECTIVE)
138143
call s_symmetry(q_prim_vf, 1, 1, k, l, pb_in, mv_in)
139144
case (BC_PERIODIC)
140-
call s_periodic(q_prim_vf, 1, 1, k, l, pb_in, mv_in)
145+
call s_periodic(q_prim_vf, 1, 1, k, l, pb_in, mv_in, q_T_sf)
141146
case (BC_SLIP_WALL)
142147
call s_slip_wall(q_prim_vf, 1, 1, k, l)
143148
case (BC_NO_SLIP_WALL)
@@ -162,20 +167,20 @@ contains
162167
#:if not MFC_CASE_OPTIMIZATION or num_dims > 1
163168

164169
if (bc_y%beg >= 0) then
165-
call s_mpi_sendrecv_variables_buffers(q_prim_vf, 2, -1, sys_size, pb_in, mv_in)
170+
call s_mpi_sendrecv_variables_buffers(q_prim_vf, 2, -1, sys_size, pb_in, mv_in, q_T_sf = q_T_sf)
166171
else
167172
$:GPU_PARALLEL_LOOP(private='[l,k]', collapse=2)
168173
do l = 0, p
169174
do k = -buff_size, m + buff_size
170175
select case (int(bc_type(2, 1)%sf(k, 0, l)))
171176
case (BC_CHAR_SUP_OUTFLOW:BC_GHOST_EXTRAP)
172-
call s_ghost_cell_extrapolation(q_prim_vf, 2, -1, k, l)
177+
call s_ghost_cell_extrapolation(q_prim_vf, 2, -1, k, l, q_T_sf)
173178
case (BC_AXIS)
174179
call s_axis(q_prim_vf, pb_in, mv_in, k, l)
175180
case (BC_REFLECTIVE)
176181
call s_symmetry(q_prim_vf, 2, -1, k, l, pb_in, mv_in)
177182
case (BC_PERIODIC)
178-
call s_periodic(q_prim_vf, 2, -1, k, l, pb_in, mv_in)
183+
call s_periodic(q_prim_vf, 2, -1, k, l, pb_in, mv_in, q_T_sf)
179184
case (BC_SLIP_WALL)
180185
call s_slip_wall(q_prim_vf, 2, -1, k, l)
181186
case (BC_NO_SLIP_WALL)
@@ -195,18 +200,18 @@ contains
195200
end if
196201

197202
if (bc_y%end >= 0) then
198-
call s_mpi_sendrecv_variables_buffers(q_prim_vf, 2, 1, sys_size, pb_in, mv_in)
203+
call s_mpi_sendrecv_variables_buffers(q_prim_vf, 2, 1, sys_size, pb_in, mv_in, q_T_sf = q_T_sf)
199204
else
200205
$:GPU_PARALLEL_LOOP(private='[l,k]', collapse=2)
201206
do l = 0, p
202207
do k = -buff_size, m + buff_size
203208
select case (int(bc_type(2, 2)%sf(k, 0, l)))
204209
case (BC_CHAR_SUP_OUTFLOW:BC_GHOST_EXTRAP)
205-
call s_ghost_cell_extrapolation(q_prim_vf, 2, 1, k, l)
210+
call s_ghost_cell_extrapolation(q_prim_vf, 2, 1, k, l, q_T_sf)
206211
case (BC_REFLECTIVE)
207212
call s_symmetry(q_prim_vf, 2, 1, k, l, pb_in, mv_in)
208213
case (BC_PERIODIC)
209-
call s_periodic(q_prim_vf, 2, 1, k, l, pb_in, mv_in)
214+
call s_periodic(q_prim_vf, 2, 1, k, l, pb_in, mv_in, q_T_sf)
210215
case (BC_SLIP_WALL)
211216
call s_slip_wall(q_prim_vf, 2, 1, k, l)
212217
case (BC_NO_SLIP_WALL)
@@ -233,18 +238,18 @@ contains
233238
#:if not MFC_CASE_OPTIMIZATION or num_dims > 2
234239

235240
if (bc_z%beg >= 0) then
236-
call s_mpi_sendrecv_variables_buffers(q_prim_vf, 3, -1, sys_size, pb_in, mv_in)
241+
call s_mpi_sendrecv_variables_buffers(q_prim_vf, 3, -1, sys_size, pb_in, mv_in, q_T_sf = q_T_sf)
237242
else
238243
$:GPU_PARALLEL_LOOP(private='[l,k]', collapse=2)
239244
do l = -buff_size, n + buff_size
240245
do k = -buff_size, m + buff_size
241246
select case (int(bc_type(3, 1)%sf(k, l, 0)))
242247
case (BC_CHAR_SUP_OUTFLOW:BC_GHOST_EXTRAP)
243-
call s_ghost_cell_extrapolation(q_prim_vf, 3, -1, k, l)
248+
call s_ghost_cell_extrapolation(q_prim_vf, 3, -1, k, l, q_T_sf)
244249
case (BC_REFLECTIVE)
245250
call s_symmetry(q_prim_vf, 3, -1, k, l, pb_in, mv_in)
246251
case (BC_PERIODIC)
247-
call s_periodic(q_prim_vf, 3, -1, k, l, pb_in, mv_in)
252+
call s_periodic(q_prim_vf, 3, -1, k, l, pb_in, mv_in, q_T_sf)
248253
case (BC_SLIP_WALL)
249254
call s_slip_wall(q_prim_vf, 3, -1, k, l)
250255
case (BC_NO_SLIP_WALL)
@@ -263,18 +268,18 @@ contains
263268
end if
264269

265270
if (bc_z%end >= 0) then
266-
call s_mpi_sendrecv_variables_buffers(q_prim_vf, 3, 1, sys_size, pb_in, mv_in)
271+
call s_mpi_sendrecv_variables_buffers(q_prim_vf, 3, 1, sys_size, pb_in, mv_in, q_T_sf = q_T_sf)
267272
else
268273
$:GPU_PARALLEL_LOOP(private='[l,k]', collapse=2)
269274
do l = -buff_size, n + buff_size
270275
do k = -buff_size, m + buff_size
271276
select case (int(bc_type(3, 2)%sf(k, l, 0)))
272277
case (BC_CHAR_SUP_OUTFLOW:BC_GHOST_EXTRAP)
273-
call s_ghost_cell_extrapolation(q_prim_vf, 3, 1, k, l)
278+
call s_ghost_cell_extrapolation(q_prim_vf, 3, 1, k, l, q_T_sf)
274279
case (BC_REFLECTIVE)
275280
call s_symmetry(q_prim_vf, 3, 1, k, l, pb_in, mv_in)
276281
case (BC_PERIODIC)
277-
call s_periodic(q_prim_vf, 3, 1, k, l, pb_in, mv_in)
282+
call s_periodic(q_prim_vf, 3, 1, k, l, pb_in, mv_in, q_T_sf)
278283
case (BC_SlIP_WALL)
279284
call s_slip_wall(q_prim_vf, 3, 1, k, l)
280285
case (BC_NO_SLIP_WALL)
@@ -296,12 +301,16 @@ contains
296301

297302
end subroutine s_populate_variables_buffers
298303

299-
subroutine s_ghost_cell_extrapolation(q_prim_vf, bc_dir, bc_loc, k, l)
304+
subroutine s_ghost_cell_extrapolation(q_prim_vf, bc_dir, bc_loc, k, l, q_T_sf)
300305
$:GPU_ROUTINE(function_name='s_ghost_cell_extrapolation', &
301306
& parallelism='[seq]', cray_inline=True)
302307
type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
308+
type(scalar_field), optional, intent(inout) :: q_T_sf
303309
integer, intent(in) :: bc_dir, bc_loc
304310
integer, intent(in) :: k, l
311+
real(wp), dimension(num_species) :: Ys_in
312+
real(wp) :: mix_mol_weight
313+
305314

306315
integer :: j, i
307316

@@ -313,13 +322,32 @@ contains
313322
q_prim_vf(i)%sf(0, k, l)
314323
end do
315324
end do
325+
if (chemistry) then
326+
do j=1, buff_size
327+
do i = chemxb,chemxe
328+
Ys_in(i - chemxb+1) = q_prim_vf(i)%sf(0, k, l)
329+
end do
330+
call get_mixture_molecular_weight(Ys_in, mix_mol_weight)
331+
q_T_sf%sf(-j,k,l) = q_prim_vf(E_idx)%sf(0, k, l)*mix_mol_weight/(gas_constant*q_prim_vf(contxb)%sf(0, k, l))
332+
end do
333+
end if
316334
else !< bc_x%end
317335
do i = 1, sys_size
318336
do j = 1, buff_size
319337
q_prim_vf(i)%sf(m + j, k, l) = &
320338
q_prim_vf(i)%sf(m, k, l)
321339
end do
322340
end do
341+
342+
if (chemistry) then
343+
do j=1, buff_size
344+
do i = chemxb,chemxe
345+
Ys_in(i - chemxb+1) = q_prim_vf(i)%sf(m, k, l)
346+
end do
347+
call get_mixture_molecular_weight(Ys_in, mix_mol_weight)
348+
q_T_sf%sf(m+j,k,l) = q_prim_vf(E_idx)%sf(m, k, l)*mix_mol_weight/(gas_constant*q_prim_vf(contxb)%sf(m, k, l))
349+
end do
350+
end if
323351
end if
324352
elseif (bc_dir == 2) then !< y-direction
325353
if (bc_loc == -1) then !< bc_y%beg
@@ -329,13 +357,31 @@ contains
329357
q_prim_vf(i)%sf(k, 0, l)
330358
end do
331359
end do
360+
if (chemistry) then
361+
do j=1, buff_size
362+
do i = chemxb,chemxe
363+
Ys_in(i - chemxb+1) = q_prim_vf(i)%sf(k, 0, l)
364+
end do
365+
call get_mixture_molecular_weight(Ys_in, mix_mol_weight)
366+
q_T_sf%sf(k,-j,l) = q_prim_vf(E_idx)%sf(k, 0, l)*mix_mol_weight/(gas_constant*q_prim_vf(contxb)%sf(k, 0, l))
367+
end do
368+
end if
332369
else !< bc_y%end
333370
do i = 1, sys_size
334371
do j = 1, buff_size
335372
q_prim_vf(i)%sf(k, n + j, l) = &
336373
q_prim_vf(i)%sf(k, n, l)
337374
end do
338375
end do
376+
if (chemistry) then
377+
do j=1, buff_size
378+
do i = chemxb,chemxe
379+
Ys_in(i - chemxb+1) = q_prim_vf(i)%sf(k, n, l)
380+
end do
381+
call get_mixture_molecular_weight(Ys_in, mix_mol_weight)
382+
q_T_sf%sf(k,n+j,l) = q_prim_vf(E_idx)%sf(k, n, l)*mix_mol_weight/(gas_constant*q_prim_vf(contxb)%sf(k, n, l))
383+
end do
384+
end if
339385
end if
340386
elseif (bc_dir == 3) then !< z-direction
341387
if (bc_loc == -1) then !< bc_z%beg
@@ -345,13 +391,31 @@ contains
345391
q_prim_vf(i)%sf(k, l, 0)
346392
end do
347393
end do
394+
if (chemistry) then
395+
do j=1, buff_size
396+
do i = chemxb,chemxe
397+
Ys_in(i - chemxb+1) = q_prim_vf(i)%sf(k, l, 0)
398+
end do
399+
call get_mixture_molecular_weight(Ys_in, mix_mol_weight)
400+
q_T_sf%sf(k, l, -j)= q_prim_vf(E_idx)%sf(k, l, 0)*mix_mol_weight/(gas_constant*q_prim_vf(contxb)%sf(k, l, 0))
401+
end do
402+
end if
348403
else !< bc_z%end
349404
do i = 1, sys_size
350405
do j = 1, buff_size
351406
q_prim_vf(i)%sf(k, l, p + j) = &
352407
q_prim_vf(i)%sf(k, l, p)
353408
end do
354409
end do
410+
if (chemistry) then
411+
do j=1, buff_size
412+
do i = chemxb,chemxe
413+
Ys_in(i - chemxb+1) = q_prim_vf(i)%sf(k, l, p)
414+
end do
415+
call get_mixture_molecular_weight(Ys_in, mix_mol_weight)
416+
q_T_sf%sf(k, l, p+j)= q_prim_vf(E_idx)%sf(k, l, p)*mix_mol_weight/(gas_constant*q_prim_vf(contxb)%sf(k, l, p))
417+
end do
418+
end if
355419
end if
356420
end if
357421

@@ -617,12 +681,15 @@ contains
617681

618682
end subroutine s_symmetry
619683

620-
subroutine s_periodic(q_prim_vf, bc_dir, bc_loc, k, l, pb_in, mv_in)
684+
subroutine s_periodic(q_prim_vf, bc_dir, bc_loc, k, l, pb_in, mv_in, q_T_sf)
621685
$:GPU_ROUTINE(parallelism='[seq]')
622686
type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
687+
type(scalar_field), intent(inout) :: q_T_sf
623688
real(stp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb_in, mv_in
624689
integer, intent(in) :: bc_dir, bc_loc
625690
integer, intent(in) :: k, l
691+
real(wp), dimension(num_species) :: Ys_in
692+
real(wp) :: mix_mol_weight
626693

627694
integer :: j, q, i
628695

@@ -635,6 +702,16 @@ contains
635702
end do
636703
end do
637704

705+
if (chemistry) then
706+
do j=1, buff_size
707+
do i = chemxb,chemxe
708+
Ys_in(i - chemxb+1) = q_prim_vf(i)%sf(m - (j - 1), k, l)
709+
end do
710+
call get_mixture_molecular_weight(Ys_in, mix_mol_weight)
711+
q_T_sf%sf(-j,k,l) = q_prim_vf(E_idx)%sf(m - (j - 1), k, l)*mix_mol_weight/(gas_constant*q_prim_vf(contxb)%sf(m - (j - 1), k, l))
712+
end do
713+
end if
714+
638715
if (qbmm .and. .not. polytropic) then
639716
do i = 1, nb
640717
do q = 1, nnode
@@ -655,6 +732,19 @@ contains
655732
end do
656733
end do
657734

735+
if (chemistry) then
736+
737+
do j=1, buff_size
738+
do i = chemxb,chemxe
739+
Ys_in(i - chemxb+1) = q_prim_vf(i)%sf(j - 1, k, l)
740+
end do
741+
call get_mixture_molecular_weight(Ys_in, mix_mol_weight)
742+
q_T_sf%sf(m + j,k,l) = q_prim_vf(E_idx)%sf(j - 1, k, l) *mix_mol_weight/(gas_constant*q_prim_vf(contxb)%sf(j - 1, k, l))
743+
end do
744+
745+
end if
746+
747+
658748
if (qbmm .and. .not. polytropic) then
659749
do i = 1, nb
660750
do q = 1, nnode

0 commit comments

Comments
 (0)