Skip to content

Commit

Permalink
Merge branch 'Galerking-index-out-of-bounds' into 'development'
Browse files Browse the repository at this point in the history
fix invalid array access for worldrank !=0

See merge request damask/DAMASK!984
  • Loading branch information
eisenlohr committed Oct 21, 2024
2 parents 96bc8c3 + 667db84 commit 907e0ba
Show file tree
Hide file tree
Showing 4 changed files with 39 additions and 24 deletions.
3 changes: 0 additions & 3 deletions src/grid/grid_mech_spectral_Galerkin.f90
Original file line number Diff line number Diff line change
Expand Up @@ -179,9 +179,6 @@ subroutine grid_mechanical_spectral_Galerkin_init(num_grid)
extmsg, &
petsc_options

KSP :: ksp
PC :: pc


print'(/,1x,a)', '<<<+- grid_mechanical_spectral_Galerkin init -+>>>'; flush(IO_STDOUT)

Expand Down
18 changes: 5 additions & 13 deletions src/grid/spectral_utilities.f90
Original file line number Diff line number Diff line change
Expand Up @@ -548,21 +548,11 @@ function utilities_G_Convolution(field,stress_mask,fieldAim) result(G_Field)
i, j, k, &
l, m

complex(pREAL), dimension(3,3) :: field_zero_freq !< zero freq of field

! print'(/,1x,a)', '... doing G convolution ...............................................'
! flush(IO_STDOUT)

tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pREAL
tensorField_real(1:3,1:3,1:cells(1), 1:cells(2),1:cells3) = field
call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier)

if (present(fieldAim)) then
field_zero_freq = cmplx(fieldAim/wgt,0.0_pREAL,pREAL)
else
field_zero_freq = tensorField_fourier(1:3,1:3,1,1,1) ! f_hat(k=0) = f_ave * vol = f_ave / wgt
endif

!$OMP PARALLEL DO PRIVATE(l,m,temp33_cmplx)
do j = 1, cells2; do k = 1, cells(3); do i = 1, cells1Red
if (any([i,j+cells2Offset,k] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
Expand All @@ -581,9 +571,11 @@ function utilities_G_Convolution(field,stress_mask,fieldAim) result(G_Field)
!$OMP END PARALLEL DO

! Apply stress bounday conditions
if (cells3Offset == 0) &
tensorField_fourier(1:3,1:3,1,1,1) = merge(cmplx(0.0_pREAL,0.0_pREAL,pREAL),field_zero_freq,stress_mask)

if (cells3Offset == 0) then
if (present(fieldAim)) &
tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pREAL,pREAL)
where (stress_mask) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(0.0_pREAL,0.0_pREAL,pREAL)
end if

call fftw_mpi_execute_dft_c2r(planTensorBack,tensorField_fourier,tensorField_real)
G_Field = tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3)
Expand Down
8 changes: 4 additions & 4 deletions src/math.f90
Original file line number Diff line number Diff line change
Expand Up @@ -353,7 +353,7 @@ end function math_cross


!--------------------------------------------------------------------------------------------------
!> @brief outer product of arbitrary sized vectors (A ⊗ B / i,j)
!> @brief Outer product of arbitrary sized vectors (A ⊗ B / i,j).
!--------------------------------------------------------------------------------------------------
pure function math_outer(A,B)

Expand All @@ -375,7 +375,7 @@ end function math_outer


!--------------------------------------------------------------------------------------------------
!> @brief inner product of arbitrary sized vectors (A · B / i,i)
!> @brief Inner product of arbitrary sized vectors (A · B / i,i).
!--------------------------------------------------------------------------------------------------
real(pREAL) pure function math_inner(A,B)

Expand All @@ -389,7 +389,7 @@ end function math_inner


!--------------------------------------------------------------------------------------------------
!> @brief double contraction of 3x3 matrices (A : B / ij,ij)
!> @brief Double contraction of 3x3 matrices (A : B / ij,ij).
!--------------------------------------------------------------------------------------------------
real(pREAL) pure function math_tensordot(A,B)

Expand All @@ -402,7 +402,7 @@ end function math_tensordot


!--------------------------------------------------------------------------------------------------
!> @brief matrix double contraction 3333x33 = 33 (ijkl,kl)
!> @brief Matrix double contraction 3333x33 = 33 (ijkl,kl).
!--------------------------------------------------------------------------------------------------
pure function math_mul3333xx33(A,B)

Expand Down
34 changes: 30 additions & 4 deletions src/misc.f90
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ module misc
module procedure misc_optional_pI32
module procedure misc_optional_pI64
module procedure misc_optional_real
module procedure misc_optional_complex
module procedure misc_optional_str
end interface misc_optional

Expand Down Expand Up @@ -44,7 +45,7 @@ end subroutine misc_init
!--------------------------------------------------------------------------------------------------
!> @brief Return bool value if given, otherwise default.
!--------------------------------------------------------------------------------------------------
pure function misc_optional_bool(given,default) result(var)
pure elemental function misc_optional_bool(given,default) result(var)

logical, intent(in), optional :: given
logical, intent(in) :: default
Expand All @@ -63,7 +64,7 @@ end function misc_optional_bool
!--------------------------------------------------------------------------------------------------
!> @brief Return integer(pI32) value if given, otherwise default.
!--------------------------------------------------------------------------------------------------
pure function misc_optional_pI32(given,default) result(var)
pure elemental function misc_optional_pI32(given,default) result(var)

integer(pI32), intent(in), optional :: given
integer(pI32), intent(in) :: default
Expand All @@ -82,7 +83,7 @@ end function misc_optional_pI32
!--------------------------------------------------------------------------------------------------
!> @brief Return integer(pI64) value if given, otherwise default.
!--------------------------------------------------------------------------------------------------
pure function misc_optional_pI64(given,default) result(var)
pure elemental function misc_optional_pI64(given,default) result(var)

integer(pI64), intent(in), optional :: given
integer(pI64), intent(in) :: default
Expand All @@ -101,7 +102,7 @@ end function misc_optional_pI64
!--------------------------------------------------------------------------------------------------
!> @brief Return real value if given, otherwise default.
!--------------------------------------------------------------------------------------------------
pure function misc_optional_real(given,default) result(var)
pure elemental function misc_optional_real(given,default) result(var)

real(pREAL), intent(in), optional :: given
real(pREAL), intent(in) :: default
Expand All @@ -117,6 +118,24 @@ pure function misc_optional_real(given,default) result(var)
end function misc_optional_real


!--------------------------------------------------------------------------------------------------
!> @brief Return complex value if given, otherwise default.
!--------------------------------------------------------------------------------------------------
pure elemental function misc_optional_complex(given,default) result(var)

complex(pREAL), intent(in), optional :: given
complex(pREAL), intent(in) :: default
complex(pREAL) :: var


if (present(given)) then
var = given
else
var = default
end if

end function misc_optional_complex

!--------------------------------------------------------------------------------------------------
!> @brief Return string value if given, otherwise default.
!--------------------------------------------------------------------------------------------------
Expand Down Expand Up @@ -192,10 +211,16 @@ end function misc_ones
subroutine misc_selfTest()

real(pREAL) :: r
real(pREAL), dimension(:), allocatable :: rN
character(len=:), allocatable :: str,out
integer :: N


call random_number(r)
N = int(r*99._pREAL)
allocate(rN(N),source=0.0_pREAL)
call random_number(rN)

if (test_str('DAMASK') /= 'DAMASK') error stop 'optional_str, present'
if (test_str() /= 'default') error stop 'optional_str, not present'
if (misc_optional(default='default') /= 'default') error stop 'optional_str, default only'
Expand All @@ -205,6 +230,7 @@ subroutine misc_selfTest()
if (dNeq(test_real(r),r)) error stop 'optional_real, present'
if (dNeq(test_real(),0.0_pREAL)) error stop 'optional_real, not present'
if (dNeq(misc_optional(default=r),r)) error stop 'optional_real, default only'
if (any(dNeq(misc_optional(default=rN),rN))) error stop 'optional_real array, default only'
if (test_bool(r<0.5_pREAL) .neqv. r<0.5_pREAL) error stop 'optional_bool, present'
if (.not. test_bool()) error stop 'optional_bool, not present'
if (misc_optional(default=r>0.5_pREAL) .neqv. r>0.5_pREAL) error stop 'optional_bool, default only'
Expand Down

0 comments on commit 907e0ba

Please sign in to comment.