223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

switched dNeq and dEq to relative tolerance, removed single precision (makes things complicated 

and was never used anyway)
This commit is contained in:
Martin Diehl 2016-03-20 23:20:58 +01:00
parent 97b52f60e7
commit 28259b2c46
3 changed files with 35 additions and 75 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
code/crystallite.f90
View File

@ -3569,11 +3569,7 @@ logical function crystallite_integrateStress(&
maxticks maxticks
external :: & external :: &
#if(FLOAT==8)
dgesv dgesv
#elif(FLOAT==4)
sgesv
#endif
!* be pessimistic !* be pessimistic
crystallite_integrateStress = .false. crystallite_integrateStress = .false.
@ -3756,11 +3752,7 @@ logical function crystallite_integrateStress(&
- math_Plain3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dT3333,dT_dFe3333),dFe_dLp3333)) - math_Plain3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dT3333,dT_dFe3333),dFe_dLp3333))
dRLp_dLp2 = dRLp_dLp ! will be overwritten in first call to LAPACK routine dRLp_dLp2 = dRLp_dLp ! will be overwritten in first call to LAPACK routine
work = math_plain33to9(residuumLp) work = math_plain33to9(residuumLp)
#if(FLOAT==8)
call dgesv(9,1,dRLp_dLp2,9,ipiv,work,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp call dgesv(9,1,dRLp_dLp2,9,ipiv,work,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp
#elif(FLOAT==4)
call sgesv(9,1,dRLp_dLp2,9,ipiv,work,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp
#endif
if (ierr /= 0_pInt) then if (ierr /= 0_pInt) then
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
@ -3849,31 +3841,27 @@ logical function crystallite_integrateStress(&
math_mul3333xx3333(dT_dFi3333, dFi_dLi3333))) & math_mul3333xx3333(dT_dFi3333, dFi_dLi3333))) &
- math_Plain3333to99(math_mul3333xx3333(dLi_dFi3333, dFi_dLi3333)) - math_Plain3333to99(math_mul3333xx3333(dLi_dFi3333, dFi_dLi3333))
work = math_plain33to9(residuumLi) work = math_plain33to9(residuumLi)
#if(FLOAT==8)
call dgesv(9,1,dRLi_dLi,9,ipiv,work,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li call dgesv(9,1,dRLi_dLi,9,ipiv,work,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li
#elif(FLOAT==4) if (ierr /= 0_pInt) then
call sgesv(9,1,dRLi_dLi,9,ipiv,work,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li
#endif
if (ierr /= 0_pInt) then
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLi inversion at el ip ipc ', & write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLi inversion at el ip ipc ', &
el,mesh_element(1,el),ip,ipc el,mesh_element(1,el),ip,ipc
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g)& .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g)&
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,*) write(6,*)
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dR_dLi',transpose(dRLi_dLi) write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dR_dLi',transpose(dRLi_dLi)
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLi',transpose(math_Plain3333to99(dFe_dLi3333)) write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLi',transpose(math_Plain3333to99(dFe_dLi3333))
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dFi_constitutive',transpose(math_Plain3333to99(dT_dFi3333)) write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dFi_constitutive',transpose(math_Plain3333to99(dT_dFi3333))
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLi_dT_constitutive',transpose(math_Plain3333to99(dLi_dT3333)) write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLi_dT_constitutive',transpose(math_Plain3333to99(dLi_dT3333))
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Li_constitutive',math_transpose33(Li_constitutive) write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Li_constitutive',math_transpose33(Li_constitutive)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Liguess',math_transpose33(Liguess) write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Liguess',math_transpose33(Liguess)
endif
endif endif
#endif
return
endif endif
#endif
return
endif
deltaLi = - math_plain9to33(work) deltaLi = - math_plain9to33(work)
endif endif

31
code/math.f90
View File

@ -186,10 +186,6 @@ module math
halton_seed_set, & halton_seed_set, &
i_to_halton, & i_to_halton, &
prime prime
external :: &
dsyev, &
dgetrf, &
dgetri
contains contains
@ -811,15 +807,13 @@ function math_invSym3333(A)
integer(pInt), dimension(6) :: ipiv6 integer(pInt), dimension(6) :: ipiv6
real(pReal), dimension(6,6) :: temp66_Real real(pReal), dimension(6,6) :: temp66_Real
real(pReal), dimension(6) :: work6 real(pReal), dimension(6) :: work6
external :: &
dgetrf, &
dgetri
temp66_real = math_Mandel3333to66(A) temp66_real = math_Mandel3333to66(A)
#if(FLOAT==8)
call dgetrf(6,6,temp66_real,6,ipiv6,ierr) call dgetrf(6,6,temp66_real,6,ipiv6,ierr)
call dgetri(6,temp66_real,6,ipiv6,work6,6,ierr) call dgetri(6,temp66_real,6,ipiv6,work6,6,ierr)
#elif(FLOAT==4)
call sgetrf(6,6,temp66_real,6,ipiv6,ierr)
call sgetri(6,temp66_real,6,ipiv6,work6,6,ierr)
#endif
if (ierr == 0_pInt) then if (ierr == 0_pInt) then
math_invSym3333 = math_Mandel66to3333(temp66_real) math_invSym3333 = math_Mandel66to3333(temp66_real)
else else
@ -847,13 +841,8 @@ subroutine math_invert(myDim,A, InvA, error)
logical, intent(out) :: error logical, intent(out) :: error
invA = A invA = A
#if(FLOAT==8)
call dgetrf(myDim,myDim,invA,myDim,ipiv,ierr) call dgetrf(myDim,myDim,invA,myDim,ipiv,ierr)
call dgetri(myDim,InvA,myDim,ipiv,work,myDim,ierr) call dgetri(myDim,InvA,myDim,ipiv,work,myDim,ierr)
#elif(FLOAT==4)
call sgetrf(myDim,myDim,invA,myDim,ipiv,ierr)
call sgetri(myDim,InvA,myDim,ipiv,work,myDim,ierr)
#endif
error = merge(.true.,.false., ierr /= 0_pInt) ! http://fortraninacworld.blogspot.de/2012/12/ternary-operator.html error = merge(.true.,.false., ierr /= 0_pInt) ! http://fortraninacworld.blogspot.de/2012/12/ternary-operator.html
end subroutine math_invert end subroutine math_invert
@ -1937,16 +1926,13 @@ subroutine math_eigenValuesVectorsSym(m,values,vectors,error)
real(pReal), dimension(size(m,1)), intent(out) :: values real(pReal), dimension(size(m,1)), intent(out) :: values
real(pReal), dimension(size(m,1),size(m,1)), intent(out) :: vectors real(pReal), dimension(size(m,1),size(m,1)), intent(out) :: vectors
logical, intent(out) :: error logical, intent(out) :: error
integer(pInt) :: info integer(pInt) :: info
real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f
external :: &
dsyev
vectors = m ! copy matrix to input (doubles as output) array vectors = m ! copy matrix to input (doubles as output) array
#if(FLOAT==8)
call dsyev('V','U',size(m,1),vectors,size(m,1),values,work,(64+2)*size(m,1),info) call dsyev('V','U',size(m,1),vectors,size(m,1),values,work,(64+2)*size(m,1),info)
#elif(FLOAT==4)
call ssyev('V','U',size(m,1),vectors,size(m,1),values,work,(64+2)*size(m,1),info)
#endif
error = (info == 0_pInt) error = (info == 0_pInt)
end subroutine math_eigenValuesVectorsSym end subroutine math_eigenValuesVectorsSym
@ -2135,16 +2121,13 @@ function math_eigenvaluesSym(m)
real(pReal), dimension(:,:), intent(in) :: m real(pReal), dimension(:,:), intent(in) :: m
real(pReal), dimension(size(m,1)) :: math_eigenvaluesSym real(pReal), dimension(size(m,1)) :: math_eigenvaluesSym
real(pReal), dimension(size(m,1),size(m,1)) :: vectors real(pReal), dimension(size(m,1),size(m,1)) :: vectors
integer(pInt) :: info integer(pInt) :: info
real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f
external :: &
dsyev
vectors = m ! copy matrix to input (doubles as output) array vectors = m ! copy matrix to input (doubles as output) array
#if(FLOAT==8)
call dsyev('N','U',size(m,1),vectors,size(m,1),math_eigenvaluesSym,work,(64+2)*size(m,1),info) call dsyev('N','U',size(m,1),vectors,size(m,1),math_eigenvaluesSym,work,(64+2)*size(m,1),info)
#elif(FLOAT==4)
call ssyev('N','U',size(m,1),vectors,size(m,1),math_eigenvaluesSym,work,(64+2)*size(m,1),info)
#endif
if (info /= 0_pInt) math_eigenvaluesSym = DAMASK_NaN if (info /= 0_pInt) math_eigenvaluesSym = DAMASK_NaN
end function math_eigenvaluesSym end function math_eigenvaluesSym

33
code/prec.f90
View File

@ -4,9 +4,9 @@
!> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH !> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH !> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH
!> @brief setting precision for real and int type depending on makros "FLOAT" and "INT" !> @brief setting precision for real and int type
!> @details setting precision for real and int type and for DAMASK_NaN. Definition is made !> @details setting precision for real and int type and for DAMASK_NaN. Definition is made
!! depending on makros "FLOAT" and "INT" defined during compilation !! depending on makro "INT" defined during compilation
!! for details on NaN see https://software.intel.com/en-us/forums/topic/294680 !! for details on NaN see https://software.intel.com/en-us/forums/topic/294680
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module prec module prec
@ -18,18 +18,7 @@ module prec
implicit none implicit none
private private
#if (FLOAT==4) #if (FLOAT==8)
#if defined(Spectral) || defined(FEM)
SPECTRAL SOLVER AND OWN FEM DO NOT SUPPORT SINGLE PRECISION, STOPPING COMPILATION
#endif
integer, parameter, public :: pReal = 4 !< floating point single precition (was selected_real_kind(6,37), number with 6 significant digits, up to 1e+-37)
#ifdef __INTEL_COMPILER
real(pReal), parameter, public :: DAMASK_NaN = Z'7F800001' !< quiet NaN for single precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran)
#endif
#ifdef __GFORTRAN__
real(pReal), parameter, public :: DAMASK_NaN = real(Z'7F800001', pReal) !< quiet NaN for single precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran)
#endif
#elif (FLOAT==8)
integer, parameter, public :: pReal = 8 !< floating point double precision (was selected_real_kind(15,300), number with 15 significant digits, up to 1e+-300) integer, parameter, public :: pReal = 8 !< floating point double precision (was selected_real_kind(15,300), number with 15 significant digits, up to 1e+-300)
#ifdef __INTEL_COMPILER #ifdef __INTEL_COMPILER
real(pReal), parameter, public :: DAMASK_NaN = Z'7FF8000000000000' !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran) real(pReal), parameter, public :: DAMASK_NaN = Z'7FF8000000000000' !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran)
@ -188,29 +177,29 @@ logical elemental pure function prec_isNaN(a)
end function prec_isNaN end function prec_isNaN
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief equality comparison for double precision !> @brief equality comparison for double precision
! replaces "==" but for certain (absolute) tolerance. Counterpart to dNeq ! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq
! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
logical elemental pure function dEq(a,b,tol) logical elemental pure function dEq(a,b,tol)
real(pReal), intent(in) :: a,b real(pReal), intent(in) :: a,b
real(pReal), intent(in), optional :: tol real(pReal), intent(in), optional :: tol
real(pReal), parameter :: eps = 1.0e-15_pReal real(pReal), parameter :: eps = 2.2204460492503131E-16 ! DBL_EPSILON in C
dEq = merge(.True., .False.,abs(a-b) < merge(tol,eps,present(tol))) dEq = merge(.True., .False.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b])))
end function dEq end function dEq
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief inequality comparison for double precision !> @brief inequality comparison for double precision
! replaces "!=" but for certain (absolute) tolerance. Counterpart to dEq ! replaces "!=" but for certain (relative) tolerance. Counterpart to dEq
! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
logical elemental pure function dNeq(a,b,tol) logical elemental pure function dNeq(a,b,tol)
real(pReal), intent(in) :: a,b real(pReal), intent(in) :: a,b
real(pReal), intent(in), optional :: tol real(pReal), intent(in), optional :: tol
real(pReal), parameter :: eps = 1.0e-15_pReal real(pReal), parameter :: eps = 2.2204460492503131E-16 ! DBL_EPSILON in C
dNeq = merge(.True., .False.,abs(a-b) < merge(tol,eps,present(tol))) dNeq = merge(.False., .True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b])))
end function dNeq end function dNeq
end module prec end module prec