Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix hypoelastic instability #773

Merged
merged 3 commits into from
Dec 26, 2024
Merged

Fix hypoelastic instability #773

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
a88a8d1
fix hypo
ChrisZYJ Dec 24, 2024
98787cc
add tests
ChrisZYJ Dec 24, 2024
3779fe3
regenerate tests
ChrisZYJ Dec 24, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
173 changes: 59 additions & 114 deletions src/simulation/m_hypoelastic.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,14 +13,13 @@

use m_global_parameters !< Definitions of the global parameters

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

! ==========================================================================

implicit none

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

real(wp), allocatable, dimension(:) :: Gs
Expand All @@ -34,16 +33,11 @@
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, k, r
integer :: i

@: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 @@ -61,29 +55,6 @@
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 @@ -99,7 +70,7 @@

real(wp) :: rho_K, G_K

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

ndirs = 1; if (n > 0) ndirs = 2; if (p > 0) ndirs = 3
Expand All @@ -112,91 +83,82 @@
do q = 0, p
do l = 0, n
do k = 0, m
du_dx(k, l, q) = 0._wp
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))

Check warning on line 91 in src/simulation/m_hypoelastic.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_hypoelastic.fpp#L86-L91 

Added lines #L86 - L91 were not covered by tests
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) = 0._wp; dv_dx(k, l, q) = 0._wp; dv_dy(k, l, q) = 0._wp
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)) &

Check warning on line 117 in src/simulation/m_hypoelastic.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_hypoelastic.fpp#L101-L117 

Added lines #L101 - L117 were not covered by tests
/(12._wp*dy(l))
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) = 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
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)) &

Check warning on line 139 in src/simulation/m_hypoelastic.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_hypoelastic.fpp#L129-L139 

Added lines #L129 - L139 were not covered by tests
/(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)) &

Check warning on line 157 in src/simulation/m_hypoelastic.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_hypoelastic.fpp#L141-L157 

Added lines #L141 - L157 were not covered by tests
/(12._wp*dz(q))
end do
end do
end do
!$acc end parallel loop
end if
end if

Expand All @@ -213,7 +175,7 @@
G_K_field(k, l, q) = G_K

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

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: 0 additions & 1 deletion src/simulation/m_start_up.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1636,7 +1636,6 @@ 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
Loading