1515!! Reference: B. Xie and F. Xiao, " Toward efficient and accurate interface
1616!! capturing on arbitrary hybrid unstructured grids: The THINC method with
1717!! quadratic surface representation and Gaussian quadrature,"
18- !! J. Comput. Phys. , vol. 349 , pp. 415-440 , 2017 .
18+ !! Journal Computational Physics , vol. 349 , pp. 415-440 , 2017 .
1919module m_thinc
2020
2121 use m_derived_types
@@ -33,13 +33,13 @@ module m_thinc
3333
3434 !> 3 - point Gauss- Legendre quadrature on [- 1 / 2 , 1 / 2 ]
3535 real (wp), parameter :: gq3_pts(3 ) = [ &
36- - 5e-1_wp * 0.7745966692414834_wp , &
37- 0._wp , &
38- 5e-1_wp * 0.7745966692414834_wp ]
36+ - 5e-1_wp * 0.7745966692414834_wp , &
37+ 0._wp , &
38+ 5e-1_wp * 0.7745966692414834_wp ]
3939 real (wp), parameter :: gq3_wts(3 ) = [ &
40- 5._wp / 18._wp , &
41- 8._wp / 18._wp , &
42- 5._wp / 18._wp ]
40+ 5._wp / 18._wp , &
41+ 8._wp / 18._wp , &
42+ 5._wp / 18._wp ]
4343 real (wp), parameter :: ln2 = 0.6931471805599453_wp
4444
4545 !> MTHINC precomputed data : unit normal components and interface
@@ -51,10 +51,6 @@ module m_thinc
5151
5252contains
5353
54- ! ================================================================
55- ! Device- callable helper functions for MTHINC
56- ! ================================================================
57-
5854 !> @brief Stable computation of ln(cosh (x)).
5955 !! For large |x| the naive formula overflows; uses
6056 !! ln(cosh (x)) = |x| + ln(1 + exp (- 2 |x|)) - ln(2 ).
@@ -82,7 +78,7 @@ contains
8278 res = 5e-1_wp * (1._wp + tanh (a))
8379 else
8480 res = 5e-1_wp + (f_log_cosh(a + 5e-1_wp * b) &
85- - f_log_cosh(a - 5e-1_wp * b))/ (2._wp * b)
81+ - f_log_cosh(a - 5e-1_wp * b))/ (2._wp * b)
8682 end if
8783 end function f_thinc_integral_1d
8884
@@ -110,7 +106,7 @@ contains
110106 do q2 = 1 , 3
111107 a = beta* (n2* gq3_pts(q1) + n3* gq3_pts(q2) + d)
112108 res = res + gq3_wts(q1)* gq3_wts(q2) &
113- * f_thinc_integral_1d(a, beta* n1)
109+ * f_thinc_integral_1d(a, beta* n1)
114110 end do
115111 end do
116112 end if
@@ -136,20 +132,20 @@ contains
136132 do q1 = 1 , 3
137133 do q2 = 1 , 3
138134 th = tanh (beta* (n1* gq3_pts(q1) &
139- + n2* gq3_pts(q2) + d))
135+ + n2* gq3_pts(q2) + d))
140136 res = res + gq3_wts(q1)* gq3_wts(q2) &
141- * (1._wp - th* th)
137+ * (1._wp - th* th)
142138 end do
143139 end do
144140 else
145141 do q1 = 1 , 3
146142 do q2 = 1 , 3
147143 do q3 = 1 , 3
148144 th = tanh (beta* (n1* gq3_pts(q1) &
149- + n2* gq3_pts(q2) &
150- + n3* gq3_pts(q3) + d))
145+ + n2* gq3_pts(q2) &
146+ + n3* gq3_pts(q3) + d))
151147 res = res + gq3_wts(q1)* gq3_wts(q2) &
152- * gq3_wts(q3)* (1._wp - th* th)
148+ * gq3_wts(q3)* (1._wp - th* th)
153149 end do
154150 end do
155151 end do
@@ -183,7 +179,7 @@ contains
183179 !! The face coordinate in face_dir is fixed; remaining directions are
184180 !! integrated over [- 1 / 2 , 1 / 2 ] (analytically along one, Gauss for others).
185181 function f_mthinc_face_average (n1 , n2 , n3 , d , beta , &
186- face_dir , face_pos , ndim ) result(res)
182+ face_dir , face_pos , ndim ) result(res)
187183 $:GPU_ROUTINE(parallelism= ' [seq]' )
188184 real (wp), intent (in ) :: n1, n2, n3, d, beta, face_pos
189185 integer , intent (in ) :: face_dir, ndim
@@ -213,17 +209,13 @@ contains
213209 res = 0._wp
214210 do q = 1 , 3
215211 a = beta* (n_face* face_pos &
216- + n_t0* gq3_pts(q) + d)
212+ + n_t0* gq3_pts(q) + d)
217213 res = res + gq3_wts(q) &
218- * f_thinc_integral_1d(a, beta* n_t1)
214+ * f_thinc_integral_1d(a, beta* n_t1)
219215 end do
220216 end if
221217 end function f_mthinc_face_average
222218
223- ! ================================================================
224- ! Module initialization / finalization
225- ! ================================================================
226-
227219 subroutine s_initialize_thinc_module ()
228220
229221 if (int_comp == 2 ) then
@@ -239,10 +231,6 @@ contains
239231
240232 end subroutine s_initialize_thinc_module
241233
242- ! ================================================================
243- ! MTHINC normal + interface position precomputation
244- ! ================================================================
245-
246234 !> @brief Compute the unit normal and interface- position parameter d
247235 !! at each interface cell. The normal is computed in cell- scaled
248236 !! coordinates (central differences without physical distance
@@ -285,18 +273,18 @@ contains
285273 ! (no division by physical distance — aspect ratio
286274 ! is embedded in the normal automatically)
287275 nr_x = (v_vf(advxb)%sf(j + 1 , k, l) &
288- - v_vf(advxb)%sf(j - 1 , k, l))* 5e-1_wp
276+ - v_vf(advxb)%sf(j - 1 , k, l))* 5e-1_wp
289277
290278 nr_y = 0._wp
291279 if (n > 0 ) then
292280 nr_y = (v_vf(advxb)%sf(j, k + 1 , l) &
293- - v_vf(advxb)%sf(j, k - 1 , l))* 5e-1_wp
281+ - v_vf(advxb)%sf(j, k - 1 , l))* 5e-1_wp
294282 end if
295283
296284 nr_z = 0._wp
297285 if (p > 0 ) then
298286 nr_z = (v_vf(advxb)%sf(j, k, l + 1 ) &
299- - v_vf(advxb)%sf(j, k, l - 1 ))* 5e-1_wp
287+ - v_vf(advxb)%sf(j, k, l - 1 ))* 5e-1_wp
300288 end if
301289
302290 nmag = sqrt (nr_x* nr_x + nr_y* nr_y + nr_z* nr_z)
@@ -311,8 +299,8 @@ contains
311299 mthinc_nhat(3 , j, k, l) = nr_z
312300
313301 mthinc_d(j, k, l) = f_mthinc_solve_d( &
314- nr_x, nr_y, nr_z, &
315- ic_beta, real (ac, wp), num_dims)
302+ nr_x, nr_y, nr_z, &
303+ ic_beta, real (ac, wp), num_dims)
316304 end if
317305
318306 end if
@@ -324,10 +312,6 @@ contains
324312
325313 end subroutine s_compute_mthinc_normals
326314
327- ! ================================================================
328- ! THINC / MTHINC interface compression
329- ! ================================================================
330-
331315 !> @brief Applies THINC (int_comp= 1 ) or MTHINC (int_comp= 2 ) interface
332316 !! compression to sharpen volume- fraction reconstructions at material
333317 !! interfaces. Called after WENO/ MUSCL reconstruction per direction.
@@ -385,7 +369,7 @@ contains
385369 if (aC >= ic_eps .and. aC <= 1._wp - ic_eps) then
386370
387371 if (int_comp == 2 .and. n > 0 ) then
388- ! ---- MTHINC: multi- dimensional face averages ----
372+ ! MTHINC: multi- dimensional face averages
389373
390374 ! Map reshaped (j,k,l) to physical (ix,iy,iz)
391375 #:if REC_DIR == 1
@@ -411,8 +395,8 @@ contains
411395
412396 ! Left face (face_pos = - 0.5 )
413397 aTHINC = f_mthinc_face_average( &
414- nh1, nh2, nh3, d_local, ic_beta, &
415- ${REC_DIR}$, - 5e-1_wp , num_dims)
398+ nh1, nh2, nh3, d_local, ic_beta, &
399+ ${REC_DIR}$, - 5e-1_wp , num_dims)
416400 if (aTHINC < ic_eps) aTHINC = ic_eps
417401 if (aTHINC > 1._wp - ic_eps) aTHINC = 1._wp - ic_eps
418402 vL_rs_vf_${XYZ}$ (j, k, l, contxb) = rho1* aTHINC
@@ -422,8 +406,8 @@ contains
422406
423407 ! Right face (face_pos = + 0.5 )
424408 aTHINC = f_mthinc_face_average( &
425- nh1, nh2, nh3, d_local, ic_beta, &
426- ${REC_DIR}$, 5e-1_wp , num_dims)
409+ nh1, nh2, nh3, d_local, ic_beta, &
410+ ${REC_DIR}$, 5e-1_wp , num_dims)
427411 if (aTHINC < ic_eps) aTHINC = ic_eps
428412 if (aTHINC > 1._wp - ic_eps) aTHINC = 1._wp - ic_eps
429413 vR_rs_vf_${XYZ}$ (j, k, l, contxb) = rho1* aTHINC
@@ -434,7 +418,7 @@ contains
434418 end if
435419
436420 else
437- ! ---- 1D THINC: constant sharpness ----
421+ ! 1D THINC: constant sharpness
438422
439423 moncon = (aCR - aC)* (aC - aCL)
440424
0 commit comments