Merge branch 'development' of magit1.mpie.de:damask/DAMASK into development

This commit is contained in:
Philip Eisenlohr 2017-09-19 08:13:08 -04:00
commit 67ac5c7b84
3 changed files with 80 additions and 25 deletions

View File

@ -1 +1 @@
v2.0.1-909-g9e502c7
v2.0.1-918-gae868d3

View File

@ -1557,13 +1557,7 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
case (400_pInt)
msg = 'matrix inversion error'
case (401_pInt)
msg = 'math_check: quat -> axisAngle -> quat failed'
case (402_pInt)
msg = 'math_check: quat -> R -> quat failed'
case (403_pInt)
msg = 'math_check: quat -> euler -> quat failed'
case (404_pInt)
msg = 'math_check: R -> euler -> R failed'
msg = 'math_check failed'
case (405_pInt)
msg = 'I_TO_HALTON-error: an input base BASE is <= 1'
case (406_pInt)

View File

@ -180,15 +180,11 @@ subroutine math_init
implicit none
integer(pInt) :: i
real(pReal), dimension(3,3) :: R,R2
real(pReal), dimension(3) :: Eulers,v
real(pReal), dimension(4) :: q,q2,axisangle,randTest
real(pReal), dimension(4) :: randTest
! the following variables are system dependend and shound NOT be pInt
integer :: randSize ! gfortran requires a variable length to compile
integer, dimension(:), allocatable :: randInit ! if recalculations of former randomness (with given seed) is necessary
! comment the first random_seed call out, set randSize to 1, and use ifort
character(len=64) :: error_msg
write(6,'(/,a)') ' <<<+- math init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
@ -211,7 +207,7 @@ subroutine math_init
enddo
write(6,'(a,I2)') ' size of random seed: ', randSize
do i =1, randSize
do i = 1_pInt,randSize
write(6,'(a,I2,I14)') ' value of random seed: ', i, randInit(i)
enddo
write(6,'(a,4(/,26x,f17.14),/)') ' start of random sequence: ', randTest
@ -221,6 +217,27 @@ subroutine math_init
call halton_seed_set(int(randInit(1), pInt))
call halton_ndim_set(3_pInt)
call math_check()
end subroutine math_init
!--------------------------------------------------------------------------------------------------
!> @brief check correctness of (some) math functions
!--------------------------------------------------------------------------------------------------
subroutine math_check
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: tol_math_check
use numerics, only: &
fixedSeed
use IO, only: IO_error, IO_timeStamp
character(len=64) :: error_msg
real(pReal), dimension(3,3) :: R,R2
real(pReal), dimension(3) :: Eulers,v
real(pReal), dimension(4) :: q,q2,axisangle
! --- check rotation dictionary ---
q = math_qRand() ! random quaternion
@ -230,7 +247,8 @@ subroutine math_init
q2 = math_axisAngleToQ(axisangle(1:3),axisangle(4))
if ( any(abs( q-q2) > tol_math_check) .and. &
any(abs(-q-q2) > tol_math_check) ) then
write (error_msg, '(a,e14.6)' ) 'maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
write (error_msg, '(a,e14.6)' ) &
'quat -> axisAngle -> quat maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
call IO_error(401_pInt,ext_msg=error_msg)
endif
@ -239,8 +257,9 @@ subroutine math_init
q2 = math_RtoQ(R)
if ( any(abs( q-q2) > tol_math_check) .and. &
any(abs(-q-q2) > tol_math_check) ) then
write (error_msg, '(a,e14.6)' ) 'maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
call IO_error(402_pInt,ext_msg=error_msg)
write (error_msg, '(a,e14.6)' ) &
'quat -> R -> quat maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
call IO_error(401_pInt,ext_msg=error_msg)
endif
! +++ q -> euler -> q +++
@ -248,28 +267,46 @@ subroutine math_init
q2 = math_EulerToQ(Eulers)
if ( any(abs( q-q2) > tol_math_check) .and. &
any(abs(-q-q2) > tol_math_check) ) then
write (error_msg, '(a,e14.6)' ) 'maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
call IO_error(403_pInt,ext_msg=error_msg)
write (error_msg, '(a,e14.6)' ) &
'quat -> euler -> quatmaximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
call IO_error(401_pInt,ext_msg=error_msg)
endif
! +++ R -> euler -> R +++
Eulers = math_RtoEuler(R)
R2 = math_EulerToR(Eulers)
if ( any(abs( R-R2) > tol_math_check) ) then
write (error_msg, '(a,e14.6)' ) 'maximum deviation ',maxval(abs( R-R2))
call IO_error(404_pInt,ext_msg=error_msg)
write (error_msg, '(a,e14.6)' ) &
'R -> euler -> R maximum deviation ',maxval(abs( R-R2))
call IO_error(401_pInt,ext_msg=error_msg)
endif
! +++ check rotation sense of q and R +++
q = math_qRand() ! random quaternion
call halton(3_pInt,v) ! random vector
R = math_qToR(q)
if (any(abs(math_mul33x3(R,v) - math_qRot(q,v)) > tol_math_check)) then
write(6,'(a,4(f8.3,1x))') 'q',q
call IO_error(409_pInt)
write (error_msg, '(a)' ) 'R(q)*v has different sense than q*v'
call IO_error(401_pInt,ext_msg=error_msg)
endif
end subroutine math_init
! +++ check vector expansion +++
if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,3.0_pReal,3.0_pReal,3.0_pReal] - &
math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1_pInt,2_pInt,3_pInt,0_pInt])) > tol_math_check)) then
write (error_msg, '(a)' ) 'math_expand [1,2,3] by [1,2,3,0] => [1,2,2,3,3,3]'
call IO_error(401_pInt,ext_msg=error_msg)
endif
if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal] - &
math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1_pInt,2_pInt])) > tol_math_check)) then
write (error_msg, '(a)' ) 'math_expand [1,2,3] by [1,2] => [1,2,2]'
call IO_error(401_pInt,ext_msg=error_msg)
endif
if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,1.0_pReal,1.0_pReal,1.0_pReal] - &
math_expand([1.0_pReal,2.0_pReal],[1_pInt,2_pInt,3_pInt])) > tol_math_check)) then
write (error_msg, '(a)' ) 'math_expand [1,2] by [1,2,3] => [1,2,2,1,1,1]'
call IO_error(401_pInt,ext_msg=error_msg)
endif
end subroutine math_check
!--------------------------------------------------------------------------------------------------
@ -334,6 +371,30 @@ recursive subroutine math_qsort(a, istart, iend)
end subroutine math_qsort
!--------------------------------------------------------------------------------------------------
!> @brief vector expansion
!> @details takes a set of numbers (a,b,c,...) and corresponding multipliers (x,y,z,...)
!> to return a vector of x times a, y times b, z times c, ...
!--------------------------------------------------------------------------------------------------
pure function math_expand(what,how)
implicit none
real(pReal), dimension(:), intent(in) :: what
integer(pInt), dimension(:), intent(in) :: how
real(pReal), dimension(sum(how)) :: math_expand
integer(pInt) :: i,j,o
o = 0_pInt
do i = 1, size(how)
do j = 1, how(i)
o = o + 1_pInt
math_expand(o) = what(1+mod(i-1,size(what)))
enddo
enddo
end function math_expand
!--------------------------------------------------------------------------------------------------
!> @brief range of integers starting at one
!--------------------------------------------------------------------------------------------------