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.

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

Revert "Fix hypoelastic instability" #776

Merged
merged 4 commits into from
Dec 27, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
174 changes: 114 additions & 60 deletions src/simulation/m_hypoelastic.fpp
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 @@
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 @@
end do
!$acc update device(Gs)

@:ALLOCATE(fd_coeff_x_h(-fd_number:fd_number, 0:m))

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

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_hypoelastic.fpp#L59

Added line #L59 was not covered by tests
if (n > 0) then
@:ALLOCATE(fd_coeff_y_h(-fd_number:fd_number, 0:n))

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

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_hypoelastic.fpp#L61

Added line #L61 was not covered by tests
end if
if (p > 0) then
@:ALLOCATE(fd_coeff_z_h(-fd_number:fd_number, 0:p))

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

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_hypoelastic.fpp#L64

Added line #L64 was not covered by tests
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 @@

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 @@
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

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

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_hypoelastic.fpp#L110

Added line #L110 was 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)

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

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_hypoelastic.fpp#L122-L123

Added lines #L122 - L123 were not covered by tests
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

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

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_hypoelastic.fpp#L136

Added line #L136 was 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_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)

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

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_hypoelastic.fpp#L148-L153

Added lines #L148 - L153 were not covered by tests
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;

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

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_hypoelastic.fpp#L167-L168

Added lines #L167 - L168 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_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)

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

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_hypoelastic.fpp#L180-L189

Added lines #L180 - L189 were not covered by tests
end do
end do
end do
end do
!$acc end parallel loop
end if
end if

Expand All @@ -171,7 +208,7 @@
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 @@

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)

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

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_hypoelastic.fpp#L338-L341

Added lines #L338 - L341 were not covered by tests
if (n > 0) then
@:DEALLOCATE(du_dy,dv_dx,dv_dy)
@:DEALLOCATE(fd_coeff_y_h)

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

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_hypoelastic.fpp#L343-L344

Added lines #L343 - L344 were not covered by tests
if (p > 0) then
@:DEALLOCATE(du_dz, dv_dz, dw_dx, dw_dy, dw_dz)
@:DEALLOCATE(fd_coeff_z_h)

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

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_hypoelastic.fpp#L346-L347

Added lines #L346 - L347 were not covered by tests
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
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
Loading