@@ -55,40 +55,52 @@ contains
5555
5656 end subroutine s_finalize_chemistry_module
5757
58- #:for NORM_DIR, XYZ in [(1 , ' x' 2 , ' y' 3 , ' z'
59- subroutine s_compute_chemistry_rhs_ ${XYZ}$ (flux_n, rhs_vf, flux_src_vf, q_prim_vf)
58+ subroutine s_compute_chemistry_advection_flux (flux_n , rhs_vf )
6059
61- type(vector_field), dimension (:), intent (IN ) :: flux_n
62- type(scalar_field), dimension (sys_size), intent (INOUT ) :: rhs_vf, flux_src_vf, q_prim_vf
63- type(int_bounds_info) :: ix, iy, iz
60+ type(vector_field), dimension (:), intent (IN ) :: flux_n
61+ type(scalar_field), dimension (sys_size), intent (INOUT ) :: rhs_vf
62+ type(int_bounds_info) :: ix, iy, iz
6463
65- integer :: x, y, z
66- integer :: eqn
67-
68- integer , parameter :: mx = ${1 if NORM_DIR == 1 else 0 }$
69- integer , parameter :: my = ${1 if NORM_DIR == 2 else 0 }$
70- integer , parameter :: mz = ${1 if NORM_DIR == 3 else 0 }$
71-
72- !$acc parallel loop collapse(4 ) present (rhs_vf, flux_n)
73- do x = 0 , m
74- do y = 0 , n
75- do z = 0 , p
76-
77- do eqn = chemxb, chemxe
78-
79- ! \nabla \cdot (F)
80- rhs_vf(eqn)%sf(x, y, z) = rhs_vf(eqn)%sf(x, y, z) + &
81- (flux_n(${NORM_DIR}$)%vf(eqn)%sf(x - mx, y - my, z - mz) - &
82- flux_n(${NORM_DIR}$)%vf(eqn)%sf(x, y, z))/ d${XYZ}$ (${XYZ}$)
64+ integer :: x, y, z
65+ integer :: eqn
8366
67+ real (kind (0d0 )) :: flux_x, flux_y, flux_z
68+
69+ #:for num_dims in range (1 , 4 )
70+ if (num_dims == ${num_dims}$) then
71+ !$acc parallel loop collapse(4 ) gang vector default(present) &
72+ !$acc private(flux_x, flux_y, flux_z)
73+ do x = 0 , m
74+ do y = 0 , n
75+ do z = 0 , p
76+ do eqn = chemxb, chemxe
77+ ! \nabla \cdot (F)
78+ flux_x = (flux_n(1 )%vf(eqn)%sf(x - 1 , y, z) - &
79+ flux_n(1 )%vf(eqn)%sf(x, y, z))/ dx(x)
80+
81+ #:if num_dims >= 2
82+ flux_x = (flux_n(2 )%vf(eqn)%sf(x, y - 1 , z) - &
83+ flux_n(2 )%vf(eqn)%sf(x, y, z))/ dy(y)
84+ #:else
85+ flux_y = 0d0
86+ #:endif
87+
88+ #:if num_dims == 3
89+ flux_z = (flux_n(3 )%vf(eqn)%sf(x, y, z - 1 ) - &
90+ flux_n(3 )%vf(eqn)%sf(x, y, z))/ dz(z)
91+ #:else
92+ flux_z = 0d0
93+ #:endif
94+
95+ rhs_vf(eqn)%sf(x, y, z) = flux_x + flux_y + flux_z
96+ end do
8497 end do
85-
8698 end do
8799 end do
88- end do
100+ end if
101+ #:endfor
89102
90- end subroutine s_compute_chemistry_rhs_ ${XYZ}$
91- #:endfor
103+ end subroutine s_compute_chemistry_advection_flux
92104
93105 subroutine s_compute_chemistry_reaction_flux (rhs_vf , q_cons_qp , q_prim_qp )
94106
@@ -104,50 +116,47 @@ contains
104116 real (kind (0d0 )) :: dyn_pres
105117 real (kind (0d0 )) :: E
106118
107- real (kind (0d0 )) :: rho
108- real (kind (1.d0 )), dimension (num_species) :: Ys
109- real (kind (1.d0 )), dimension (num_species) :: omega
119+ real (kind (0d0 )) :: rho, omega_m, omega_T
120+ real (kind (0d0 )), dimension (num_species) :: Ys
121+ real (kind (0d0 )), dimension (num_species) :: omega
110122 real (kind (0d0 )), dimension (num_species) :: enthalpies
111123 real (kind (0d0 )) :: cp_mix
112124
113125 #:if chemistry
114126
115- !$acc parallel loop collapse(4 ) private(rho)
127+ !$acc parallel loop collapse(3 ) private(rho)
116128 do x = 0 , m
117129 do y = 0 , n
118130 do z = 0 , p
119131
120- ! Maybe use q_prim_vf instead?
121132 rho = 0d0
122133 do eqn = chemxb, chemxe
123134 rho = rho + q_cons_qp(eqn)%sf(x, y, z)
124135 end do
125136
126137 do eqn = chemxb, chemxe
127- Ys(eqn - chemxb + 1 ) = q_cons_qp (eqn)%sf(x, y, z)/ rho
138+ Ys(eqn - chemxb + 1 ) = q_prim_qp (eqn)%sf(x, y, z)
128139 end do
129140
130141 dyn_pres = 0d0
131-
132142 do i = momxb, momxe
133- dyn_pres = dyn_pres + rho * q_cons_qp(i)%sf(x, y, z)* &
134- q_cons_qp (i)%sf(x, y, z)/ 2d0
143+ dyn_pres = dyn_pres + q_cons_qp(i)%sf(x, y, z)* &
144+ q_prim_qp (i)%sf(x, y, z)/ 2d0
135145 end do
136146
137- call get_temperature(.true. , q_cons_qp(E_idx)%sf(x, y, z) - dyn_pres, &
138- & q_prim_qp(tempxb)%sf(x, y, z), Ys , T)
147+ call get_temperature(q_cons_qp(E_idx)%sf(x, y, z)/ rho - dyn_pres/ rho , &
148+ & 1200d0 , Ys, .true. , T)
139149
140150 call get_net_production_rates(rho, T, Ys, omega)
141151
142- q_cons_qp(tempxb)%sf(x, y, z) = T
143- q_prim_qp(tempxb)%sf(x, y, z) = T
144-
145- !print * , x, y, z, T, rho, Ys, omega, q_cons_qp(E_idx)%sf(x, y, z), dyn_pres
152+ call get_species_enthalpies_rt(T, enthalpies)
146153
147154 do eqn = chemxb, chemxe
148155
149- rhs_vf(eqn)%sf(x, y, z) = rhs_vf(eqn)%sf(x, y, z) + &
150- mol_weights(eqn - chemxb + 1 )* omega(eqn - chemxb + 1 )
156+ omega_m = mol_weights(eqn - chemxb + 1 )* omega(eqn - chemxb + 1 )
157+ omega_T = omega_m* enthalpies(eqn - chemxb + 1 )* gas_constant* T
158+
159+ rhs_vf(eqn)%sf(x, y, z) = rhs_vf(eqn)%sf(x, y, z) + omega_m
151160
152161 end do
153162
@@ -171,9 +180,9 @@ contains
171180 integer :: eqn
172181
173182 !$acc parallel loop collapse(4 )
174- do x = ix%beg, ix% end
175- do y = iy%beg, iy%end
176- do z = iz%beg, iz%end
183+ do x = 0 , m
184+ do y = 0 , n
185+ do z = 0 , p
177186
178187 do eqn = chemxb, chemxe
179188 q_cons_qp(eqn)%sf(x, y, z) = max (0d0 , q_cons_qp(eqn)%sf(x, y, z))
0 commit comments