Skip to content

Commit

Permalink
Revert "Fix hypoelastic instability" (#776)
Browse files Browse the repository at this point in the history
  • Loading branch information
@sbryngelson
sbryngelson authored Dec 27, 2024
1 parent 5d8c07c commit d667ea7
Show file tree
Hide file tree
Showing 8 changed files with 167 additions and 565 deletions.
174 changes: 114 additions & 60 deletions src/simulation/m_hypoelastic.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,14 +8,14 @@
module m_hypoelastic

use m_derived_types !< Definitions of the derived types

use m_global_parameters !< Definitions of the global parameters

use m_mpi_proxy !< Message passing interface (MPI) module proxy
use m_finite_differences
use m_helper

implicit none

private; public :: s_initialize_hypoelastic_module, &
s_finalize_hypoelastic_module, &
s_compute_hypoelastic_rhs

real(wp), allocatable, dimension(:) :: Gs
Expand All @@ -29,11 +29,16 @@ module m_hypoelastic
real(wp), allocatable, dimension(:, :, :) :: rho_K_field, G_K_field
!$acc declare create(rho_K_field, G_K_field)

real(wp), allocatable, dimension(:, :) :: fd_coeff_x_h
real(wp), allocatable, dimension(:, :) :: fd_coeff_y_h
real(wp), allocatable, dimension(:, :) :: fd_coeff_z_h
!$acc declare create(fd_coeff_x_h,fd_coeff_y_h,fd_coeff_z_h)

contains

subroutine s_initialize_hypoelastic_module

integer :: i
integer :: i, k, r

@:ALLOCATE(Gs(1:num_fluids))
@:ALLOCATE(rho_K_field(0:m,0:n,0:p), G_K_field(0:m,0:n,0:p))
Expand All @@ -51,6 +56,29 @@ contains
end do
!$acc update device(Gs)

@:ALLOCATE(fd_coeff_x_h(-fd_number:fd_number, 0:m))
if (n > 0) then
@:ALLOCATE(fd_coeff_y_h(-fd_number:fd_number, 0:n))
end if
if (p > 0) then
@:ALLOCATE(fd_coeff_z_h(-fd_number:fd_number, 0:p))
end if

! Computing centered finite difference coefficients
call s_compute_finite_difference_coefficients(m, x_cc, fd_coeff_x_h, buff_size, &
fd_number, fd_order)
!$acc update device(fd_coeff_x_h)
if (n > 0) then
call s_compute_finite_difference_coefficients(n, y_cc, fd_coeff_y_h, buff_size, &
fd_number, fd_order)
!$acc update device(fd_coeff_y_h)
end if
if (p > 0) then
call s_compute_finite_difference_coefficients(p, z_cc, fd_coeff_z_h, buff_size, &
fd_number, fd_order)
!$acc update device(fd_coeff_z_h)
end if

end subroutine s_initialize_hypoelastic_module

!> The purpose of this procedure is to compute the source terms
Expand All @@ -66,7 +94,7 @@ contains

real(wp) :: rho_K, G_K

integer :: i, k, l, q !< Loop variables
integer :: i, k, l, q, r !< Loop variables
integer :: ndirs !< Number of coordinate directions

ndirs = 1; if (n > 0) ndirs = 2; if (p > 0) ndirs = 3
Expand All @@ -79,82 +107,91 @@ contains
do q = 0, p
do l = 0, n
do k = 0, m
du_dx(k, l, q) = &
(q_prim_vf(momxb)%sf(k - 2, l, q) &
- 8._wp*q_prim_vf(momxb)%sf(k - 1, l, q) &
+ 8._wp*q_prim_vf(momxb)%sf(k + 1, l, q) &
- q_prim_vf(momxb)%sf(k + 2, l, q)) &
/(12._wp*dx(k))
du_dx(k, l, q) = 0._wp
end do
end do
end do
!$acc end parallel loop

!$acc parallel loop collapse(3) gang vector default(present)
do q = 0, p
do l = 0, n
do k = 0, m
!$acc loop seq
do r = -fd_number, fd_number
du_dx(k, l, q) = du_dx(k, l, q) &
+ q_prim_vf(momxb)%sf(k + r, l, q)*fd_coeff_x_h(r, k)
end do

end do
end do
end do
!$acc end parallel loop

if (ndirs > 1) then
!$acc parallel loop collapse(3) gang vector default(present)
do q = 0, p
do l = 0, n
do k = 0, m
du_dy(k, l, q) = &
(q_prim_vf(momxb)%sf(k, l - 2, q) &
- 8._wp*q_prim_vf(momxb)%sf(k, l - 1, q) &
+ 8._wp*q_prim_vf(momxb)%sf(k, l + 1, q) &
- q_prim_vf(momxb)%sf(k, l + 2, q)) &
/(12._wp*dy(l))
dv_dx(k, l, q) = &
(q_prim_vf(momxb + 1)%sf(k - 2, l, q) &
- 8._wp*q_prim_vf(momxb + 1)%sf(k - 1, l, q) &
+ 8._wp*q_prim_vf(momxb + 1)%sf(k + 1, l, q) &
- q_prim_vf(momxb + 1)%sf(k + 2, l, q)) &
/(12._wp*dx(k))
dv_dy(k, l, q) = &
(q_prim_vf(momxb + 1)%sf(k, l - 2, q) &
- 8._wp*q_prim_vf(momxb + 1)%sf(k, l - 1, q) &
+ 8._wp*q_prim_vf(momxb + 1)%sf(k, l + 1, q) &
- q_prim_vf(momxb + 1)%sf(k, l + 2, q)) &
/(12._wp*dy(l))
du_dy(k, l, q) = 0._wp; dv_dx(k, l, q) = 0._wp; dv_dy(k, l, q) = 0._wp
end do
end do
end do
!$acc end parallel loop

!$acc parallel loop collapse(3) gang vector default(present)
do q = 0, p
do l = 0, n
do k = 0, m
!$acc loop seq
do r = -fd_number, fd_number
du_dy(k, l, q) = du_dy(k, l, q) &
+ q_prim_vf(momxb)%sf(k, l + r, q)*fd_coeff_y_h(r, l)
dv_dx(k, l, q) = dv_dx(k, l, q) &
+ q_prim_vf(momxb + 1)%sf(k + r, l, q)*fd_coeff_x_h(r, k)
dv_dy(k, l, q) = dv_dy(k, l, q) &
+ q_prim_vf(momxb + 1)%sf(k, l + r, q)*fd_coeff_y_h(r, l)
end do
end do
end do
end do
!$acc end parallel loop

! 3D
if (ndirs == 3) then

!$acc parallel loop collapse(3) gang vector default(present)
do q = 0, p
do l = 0, n
do k = 0, m
du_dz(k, l, q) = &
(q_prim_vf(momxb)%sf(k, l, q - 2) &
- 8._wp*q_prim_vf(momxb)%sf(k, l, q - 1) &
+ 8._wp*q_prim_vf(momxb)%sf(k, l, q + 1) &
- q_prim_vf(momxb)%sf(k, l, q + 2)) &
/(12._wp*dz(q))
dv_dz(k, l, q) = &
(q_prim_vf(momxb + 1)%sf(k, l, q - 2) &
- 8._wp*q_prim_vf(momxb + 1)%sf(k, l, q - 1) &
+ 8._wp*q_prim_vf(momxb + 1)%sf(k, l, q + 1) &
- q_prim_vf(momxb + 1)%sf(k, l, q + 2)) &
/(12._wp*dz(q))
dw_dx(k, l, q) = &
(q_prim_vf(momxe)%sf(k - 2, l, q) &
- 8._wp*q_prim_vf(momxe)%sf(k - 1, l, q) &
+ 8._wp*q_prim_vf(momxe)%sf(k + 1, l, q) &
- q_prim_vf(momxe)%sf(k + 2, l, q)) &
/(12._wp*dx(k))
dw_dy(k, l, q) = &
(q_prim_vf(momxe)%sf(k, l - 2, q) &
- 8._wp*q_prim_vf(momxe)%sf(k, l - 1, q) &
+ 8._wp*q_prim_vf(momxe)%sf(k, l + 1, q) &
- q_prim_vf(momxe)%sf(k, l + 2, q)) &
/(12._wp*dy(l))
dw_dz(k, l, q) = &
(q_prim_vf(momxe)%sf(k, l, q - 2) &
- 8._wp*q_prim_vf(momxe)%sf(k, l, q - 1) &
+ 8._wp*q_prim_vf(momxe)%sf(k, l, q + 1) &
- q_prim_vf(momxe)%sf(k, l, q + 2)) &
/(12._wp*dz(q))
du_dz(k, l, q) = 0_wp; dv_dz(k, l, q) = 0_wp; dw_dx(k, l, q) = 0_wp;
dw_dy(k, l, q) = 0_wp; dw_dz(k, l, q) = 0_wp;
end do
end do
end do
!$acc end parallel loop

!$acc parallel loop collapse(3) gang vector default(present)
do q = 0, p
do l = 0, n
do k = 0, m
!$acc loop seq
do r = -fd_number, fd_number
du_dz(k, l, q) = du_dz(k, l, q) &
+ q_prim_vf(momxb)%sf(k, l, q + r)*fd_coeff_z_h(r, q)
dv_dz(k, l, q) = dv_dz(k, l, q) &
+ q_prim_vf(momxb + 1)%sf(k, l, q + r)*fd_coeff_z_h(r, q)
dw_dx(k, l, q) = dw_dx(k, l, q) &
+ q_prim_vf(momxe)%sf(k + r, l, q)*fd_coeff_x_h(r, k)
dw_dy(k, l, q) = dw_dy(k, l, q) &
+ q_prim_vf(momxe)%sf(k, l + r, q)*fd_coeff_y_h(r, l)
dw_dz(k, l, q) = dw_dz(k, l, q) &
+ q_prim_vf(momxe)%sf(k, l, q + r)*fd_coeff_z_h(r, q)
end do
end do
end do
end do
!$acc end parallel loop
end if
end if

Expand All @@ -171,7 +208,7 @@ contains
G_K_field(k, l, q) = G_K

!TODO: take this out if not needed
if (G_K < 1000) then
if (G_K < verysmall) then
G_K_field(k, l, q) = 0
end if
end do
Expand Down Expand Up @@ -296,4 +333,21 @@ contains

end subroutine s_compute_hypoelastic_rhs

subroutine s_finalize_hypoelastic_module()

@:DEALLOCATE(Gs)
@:DEALLOCATE(rho_K_field, G_K_field)
@:DEALLOCATE(du_dx)
@:DEALLOCATE(fd_coeff_x_h)
if (n > 0) then
@:DEALLOCATE(du_dy,dv_dx,dv_dy)
@:DEALLOCATE(fd_coeff_y_h)
if (p > 0) then
@:DEALLOCATE(du_dz, dv_dz, dw_dx, dw_dy, dw_dz)
@:DEALLOCATE(fd_coeff_z_h)
end if
end if

end subroutine s_finalize_hypoelastic_module

end module m_hypoelastic
1 change: 1 addition & 0 deletions src/simulation/m_start_up.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1627,6 +1627,7 @@ contains
subroutine s_finalize_modules

call s_finalize_time_steppers_module()
if (hypoelasticity) call s_finalize_hypoelastic_module()
if (hyperelasticity) call s_finalize_hyperelastic_module()
call s_finalize_derived_variables_module()
call s_finalize_data_output_module()
Expand Down
Loading

0 comments on commit d667ea7

Please sign in to comment.