Skip to content

Commit 328c807

Browse files
sbryngelsonsbryngelson
committed
remove surface tension
1 parent ecd0c1c commit 328c807

File tree

9 files changed

+41
-2214
lines changed

9 files changed

+41
-2214
lines changed

‎src/post_process/m_derived_variables.f90‎

Lines changed: 0 additions & 179 deletions
Original file line numberDiff line numberDiff line change
@@ -21,14 +21,12 @@ module m_derived_variables
2121

2222
private; public :: s_initialize_derived_variables_module, &
2323
s_compute_finite_difference_coefficients, &
24-
s_derive_unadvected_volume_fraction, &
2524
s_derive_specific_heat_ratio, &
2625
s_derive_liquid_stiffness, &
2726
s_derive_sound_speed, &
2827
s_derive_flux_limiter, &
2928
s_derive_vorticity_component, &
3029
s_derive_numerical_schlieren_function, &
31-
s_derive_curvature, &
3230
s_finalize_derived_variables_module
3331

3432
real(kind(0d0)), allocatable, dimension(:, :, :) :: gm_rho_sf !<
@@ -168,43 +166,6 @@ subroutine s_compute_finite_difference_coefficients(q, offset_s, &
168166

169167
end subroutine s_compute_finite_difference_coefficients ! --------------
170168

171-
!> This subroutine is used together with the volume fraction
172-
!! model and when called upon, it computes the values of the
173-
!! unadvected volume fraction from the inputted conservative
174-
!! variables, q_cons_vf. The calculated values are stored in
175-
!! the derived flow quantity storage variable, q_sf.
176-
!! @param q_cons_vf Conservative variables
177-
!! @param q_sf Unadvected volume fraction
178-
subroutine s_derive_unadvected_volume_fraction(q_cons_vf, q_sf) ! ------
179-
180-
type(scalar_field), &
181-
dimension(sys_size), &
182-
intent(IN) :: q_cons_vf
183-
184-
real(kind(0d0)), &
185-
dimension(-offset_x%beg:m + offset_x%end, &
186-
-offset_y%beg:n + offset_y%end, &
187-
-offset_z%beg:p + offset_z%end), &
188-
intent(INOUT) :: q_sf
189-
190-
integer :: i, j, k, l !< Generic loop iterators
191-
192-
! Computing unadvected volume fraction from conservative variables
193-
do l = -offset_z%beg, p + offset_z%end
194-
do k = -offset_y%beg, n + offset_y%end
195-
do j = -offset_x%beg, m + offset_x%end
196-
197-
q_sf(j, k, l) = 1d0
198-
199-
do i = adv_idx%beg, adv_idx%end
200-
q_sf(j, k, l) = q_sf(j, k, l) - q_cons_vf(i)%sf(j, k, l)
201-
end do
202-
203-
end do
204-
end do
205-
end do
206-
207-
end subroutine s_derive_unadvected_volume_fraction ! -------------------
208169

209170
!> This subroutine receives as input the specific heat ratio
210171
!! function, gamma_sf, and derives from it the specific heat
@@ -438,146 +399,6 @@ subroutine s_derive_flux_limiter(i, q_prim_vf, q_sf) ! -----------------
438399
end do
439400
end subroutine s_derive_flux_limiter ! ---------------------------------
440401

441-
!> Computes the local curvatures
442-
!! @param i Fluid indicator
443-
!! @param q_prim_vf Primitive variables
444-
!! @param q_sf Curvature
445-
subroutine s_derive_curvature(i, q_prim_vf, q_sf) ! --------------------
446-
447-
integer, intent(IN) :: i
448-
449-
type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf
450-
451-
real(kind(0d0)), dimension(-offset_x%beg:m + offset_x%end, &
452-
-offset_y%beg:n + offset_y%end, &
453-
-offset_z%beg:p + offset_z%end), &
454-
intent(INOUT) :: q_sf
455-
456-
real(kind(0d0)), allocatable, dimension(:, :, :) :: alpha_unadv_sf
457-
458-
real(kind(0d0)), allocatable, dimension(:, :) :: A
459-
real(kind(0d0)), allocatable, dimension(:) :: sol, b, AA
460-
real(kind(0d0)) :: norm
461-
real(kind(0d0)) :: xloc, yloc, zloc
462-
integer :: ndim
463-
integer :: j, k, l, jj, kk, ll, i1, i2
464-
integer :: stencil_j_min, stencil_j_max
465-
integer :: stencil_k_min, stencil_k_max
466-
integer :: stencil_l_min, stencil_l_max
467-
type(int_bounds_info) :: iz1
468-
469-
if (p > 0) then ! 3D simulation
470-
allocate (A(10, 10))
471-
allocate (sol(10)); allocate (b(10)); allocate (AA(10))
472-
else ! 2D simulation
473-
allocate (A(6, 6))
474-
allocate (sol(6)); allocate (b(6)); allocate (AA(6))
475-
end if
476-
477-
if ((i == num_fluids) .and. (adv_alphan .neqv. .true.)) then
478-
allocate (alpha_unadv_sf(-offset_x%beg:m + offset_x%end, &
479-
-offset_y%beg:n + offset_y%end, &
480-
-offset_z%beg:p + offset_z%end))
481-
call s_derive_unadvected_volume_fraction(q_prim_vf, alpha_unadv_sf)
482-
end if
483-
484-
if (p > 0) then
485-
iz1%beg = -offset_z%beg; iz1%end = p + offset_z%end
486-
else
487-
iz1%beg = -1; iz1%end = 1
488-
end if
489-
490-
! Parabolic fitting
491-
ndim = size(sol, 1)
492-
493-
do l = iz1%beg + 1, iz1%end - 1
494-
do k = -offset_y%beg + 1, n + offset_y%end - 1
495-
do j = -offset_x%beg + 1, m + offset_x%end - 1
496-
A(:, :) = 0d0
497-
b(:) = 0d0
498-
sol(:) = 0d0
499-
AA(:) = 0d0
500-
501-
stencil_j_min = j - 1; stencil_j_max = j + 1
502-
stencil_k_min = k - 1; stencil_k_max = k + 1
503-
if (p > 0) then
504-
stencil_l_min = l - 1; stencil_l_max = l + 1
505-
else
506-
stencil_l_min = 0; stencil_l_max = 0
507-
end if
508-
509-
do ll = stencil_l_min, stencil_l_max
510-
do kk = stencil_k_min, stencil_k_max
511-
do jj = stencil_j_min, stencil_j_max
512-
513-
! Ignore corner points in 3D stencil
514-
if ((p > 0) .and. (abs(jj - j) == 1) .and. (abs(kk - k) == 1) .and. (abs(ll - l) == 1)) cycle
515-
516-
! Find distance between cell centers
517-
xloc = x_cc(jj) - x_cc(j)
518-
yloc = y_cc(kk) - y_cc(k)
519-
if (p > 0) then
520-
zloc = z_cc(ll) - z_cc(l)
521-
end if
522-
523-
! Compute operator
524-
AA(1) = 1d0
525-
AA(2) = xloc
526-
AA(3) = yloc
527-
AA(4) = 5d-1*xloc**2d0
528-
AA(5) = 5d-1*yloc**2d0
529-
AA(6) = xloc*yloc
530-
if (p > 0) then
531-
AA(7) = zloc
532-
AA(8) = 5d-1*zloc**2d0
533-
AA(9) = yloc*zloc
534-
AA(10) = zloc*xloc
535-
end if
536-
537-
do i1 = 1, ndim
538-
do i2 = 1, ndim
539-
A(i1, i2) = A(i1, i2) + AA(i1)*AA(i2)
540-
end do
541-
end do
542-
543-
! Form RHS vector
544-
do i1 = 1, ndim
545-
if ((i == num_fluids) .and. (adv_alphan .neqv. .true.)) then
546-
b(i1) = b(i1) + alpha_unadv_sf(jj, kk, ll)*AA(i1)
547-
else
548-
b(i1) = b(i1) + q_prim_vf(E_idx + i)%sf(jj, kk, ll)*AA(i1)
549-
end if
550-
end do
551-
end do
552-
end do
553-
end do
554-
555-
call s_solve_linear_system(A, b, sol, ndim)
556-
557-
if (p > 0) then
558-
norm = sqrt(sol(2)**2d0 + sol(3)**2d0 + sol(7)**2d0)
559-
else
560-
norm = sqrt(sol(2)**2d0 + sol(3)**2d0)
561-
end if
562-
563-
if (p > 0) then
564-
q_sf(j, k, l) = -(sol(2)**2d0*sol(5) - 2d0*sol(2)*sol(3)*sol(6) + sol(2)**2d0*sol(8) + &
565-
sol(3)**2d0*sol(8) - 2d0*sol(3)*sol(7)*sol(9) + sol(3)**2d0*sol(4) + &
566-
sol(7)**2d0*sol(4) - 2d0*sol(7)*sol(2)*sol(10) + sol(7)**2d0*sol(5))/ &
567-
max(norm, sgm_eps)**3d0
568-
else
569-
q_sf(j, k, l) = -(sol(2)**2d0*sol(5) - 2d0*sol(2)*sol(3)*sol(6) + sol(3)**2d0*sol(4))/ &
570-
max(norm, sgm_eps)**3d0
571-
end if
572-
573-
end do
574-
end do
575-
end do
576-
577-
deallocate (A, sol, b, AA)
578-
if ((i == num_fluids) .and. (adv_alphan .neqv. .true.)) deallocate (alpha_unadv_sf)
579-
580-
end subroutine s_derive_curvature ! ------------------------------------
581402

582403
!> Computes the solution to the linear system Ax=b w/ sol = x
583404
!! @param A Input matrix

0 commit comments

Comments
 (0)