some minor improvement on precision detection: checking only once (in prec and no longer in math and crystallite), added one more 4/8 switch for LAPACK, as there is no single precision FFTW, stopping compilation for spectral method if FLOAT=4

new function in IO to print integers without leading zeros, implemented it at some places in the new spectral solver (reporting still needs some serious polishing)
updated preprocessing for documentation to handle precision correctly
This commit is contained in:
Martin Diehl 2012-08-30 20:26:28 +00:00
parent b7dc9f9944
commit 22812c9a91
6 changed files with 331 additions and 315 deletions

View File

@ -31,7 +31,8 @@ program DAMASK_spectral_Driver
IO_error, & IO_error, &
IO_lc, & IO_lc, &
IO_read_jobBinaryFile, & IO_read_jobBinaryFile, &
IO_write_jobBinaryFile IO_write_jobBinaryFile, &
IO_intOut
use math use math
@ -390,9 +391,11 @@ program DAMASK_spectral_Driver
write(6,'(a)') '==================================================================' write(6,'(a)') '=================================================================='
if(solres%converged) then if(solres%converged) then
convergedCounter = convergedCounter + 1_pInt convergedCounter = convergedCounter + 1_pInt
write(6,'(A,I5.5,A)') 'increment ', totalIncsCounter, ' converged' write(6,'(A,'//IO_intOut(totalIncsCounter)//',A)') &
'increment', totalIncsCounter, 'converged'
else else
write(6,'(A,I5.5,A)') 'increment ', totalIncsCounter, ' NOT converged' write(6,'(A,'//IO_intOut(totalIncsCounter)//',A)') &
'increment', totalIncsCounter, 'NOT converged'
notConvergedCounter = notConvergedCounter + 1_pInt notConvergedCounter = notConvergedCounter + 1_pInt
endif endif

View File

@ -47,7 +47,8 @@ subroutine basic_init()
use IO, only: & use IO, only: &
IO_read_JobBinaryFile, & IO_read_JobBinaryFile, &
IO_write_JobBinaryFile IO_write_JobBinaryFile, &
IO_intOut
use FEsolving, only: & use FEsolving, only: &
restartInc restartInc
@ -95,8 +96,8 @@ subroutine basic_init()
- geomdim/real(2_pInt*res,pReal) - geomdim/real(2_pInt*res,pReal)
enddo; enddo; enddo enddo; enddo; enddo
elseif (restartInc > 1_pInt) then ! using old values from file elseif (restartInc > 1_pInt) then ! using old values from file
if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',& if (debugRestart) write(6,'(a,'//IO_intOut(restartInc-1_pInt)//',a)') &
restartInc - 1_pInt,' from file' 'Reading values of increment', restartInc - 1_pInt, 'from file'
call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',& call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',&
trim(getSolverJobName()),size(F)) trim(getSolverJobName()),size(F))
read (777,rec=1) F read (777,rec=1) F
@ -154,7 +155,8 @@ type(solutionState) function basic_solution(guessmode,timeinc,timeinc_old,P_BC,F
geomdim, & geomdim, &
deformed_fft deformed_fft
use IO, only: & use IO, only: &
IO_write_JobBinaryFile IO_write_JobBinaryFile, &
IO_intOut
use DAMASK_spectral_Utilities, only: & use DAMASK_spectral_Utilities, only: &
boundaryCondition, & boundaryCondition, &
@ -240,7 +242,7 @@ type(solutionState) function basic_solution(guessmode,timeinc,timeinc_old,P_BC,F
! report begin of new iteration ! report begin of new iteration
write(6,'(a)') '' write(6,'(a)') ''
write(6,'(a)') '==================================================================' write(6,'(a)') '=================================================================='
write(6,'(3(a,i6.6))') ' Iter. ',itmin,' < ',iter,' < ',itmax + 1_pInt write(6,'(3(a,'//IO_intOut(itmax)//'))') ' Iter.', itmin, '<',iter, '<', itmax + 1_pInt
write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =', & write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =', &
math_transpose33(F_aim) math_transpose33(F_aim)
F_aim_lab_lastIter = math_rotate_backward33(F_aim,rotation_BC) F_aim_lab_lastIter = math_rotate_backward33(F_aim,rotation_BC)
@ -279,7 +281,6 @@ end function basic_solution
!> @brief convergence check for basic scheme based on div of P and deviation from stress aim !> @brief convergence check for basic scheme based on div of P and deviation from stress aim
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
logical function basic_Converged(err_div,pAvgDiv,err_stress,pAvgStress) logical function basic_Converged(err_div,pAvgDiv,err_stress,pAvgStress)
use numerics, only: & use numerics, only: &
itmin, & itmin, &
err_div_tol, & err_div_tol, &
@ -292,7 +293,6 @@ logical function basic_Converged(err_div,pAvgDiv,err_stress,pAvgStress)
math_transpose33 math_transpose33
implicit none implicit none
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
pAvgDiv,& pAvgDiv,&
pAvgStress pAvgStress

View File

@ -62,7 +62,8 @@ module IO
IO_countContinuousIntValues, & IO_countContinuousIntValues, &
IO_continuousIntValues, & IO_continuousIntValues, &
IO_error, & IO_error, &
IO_warning IO_warning, &
IO_intOut
#ifndef Spectral #ifndef Spectral
public :: IO_open_inputFile, & public :: IO_open_inputFile, &
IO_open_logFile IO_open_logFile
@ -1267,6 +1268,17 @@ function IO_continuousIntValues(myUnit,maxN,lookupName,lookupMap,lookupMaxN)
100 end function IO_continuousIntValues 100 end function IO_continuousIntValues
pure function IO_intOut(intToPrint)
implicit none
character(len=64) :: N_Digits
character(len=64) :: IO_intOut
integer(pInt), intent(in) :: intToPrint
write(N_Digits, '(I16.16)') 1_pInt + int(log10(real(intToPrint)),pInt)
IO_intOut = '1x,I'//trim(N_Digits)//'.'//trim(N_Digits)//',1x'
end function IO_intOut
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief write error statements to standard out and terminate the Marc/spectral run with exit #9xxx !> @brief write error statements to standard out and terminate the Marc/spectral run with exit #9xxx
!> in ABAQUS either time step is reduced or execution terminated !> in ABAQUS either time step is reduced or execution terminated
@ -1384,8 +1396,6 @@ subroutine IO_error(error_ID,e,i,g,ext_msg)
!* errors related to spectral solver !* errors related to spectral solver
case (808_pInt)
msg = 'precision not suitable for FFTW'
case (809_pInt) case (809_pInt)
msg = 'initializing FFTW' msg = 'initializing FFTW'
case (831_pInt) case (831_pInt)

View File

@ -30,7 +30,7 @@
!* - _postResults * !* - _postResults *
!*************************************** !***************************************
MODULE crystallite module crystallite
use prec, only: pReal, pInt use prec, only: pReal, pInt
@ -46,62 +46,61 @@ private :: crystallite_integrateStateFPI, &
! **************************************************************** ! ****************************************************************
! *** General variables for the crystallite calculation *** ! *** General variables for the crystallite calculation ***
! **************************************************************** ! ****************************************************************
integer(pInt) crystallite_maxSizePostResults integer(pInt) crystallite_maxSizePostResults
integer(pInt), dimension(:), allocatable :: crystallite_sizePostResults integer(pInt), dimension(:), allocatable :: crystallite_sizePostResults
integer(pInt), dimension(:,:), allocatable :: crystallite_sizePostResult integer(pInt), dimension(:,:), allocatable :: crystallite_sizePostResult
character(len=64), dimension(:,:), allocatable :: crystallite_output ! name of each post result output character(len=64), dimension(:,:), allocatable :: crystallite_output !< name of each post result output
integer(pInt), dimension (:,:,:), allocatable :: & integer(pInt), dimension (:,:,:), allocatable :: &
crystallite_symmetryID ! crystallographic symmetry 1=cubic 2=hexagonal, needed in all orientation calcs crystallite_symmetryID !< crystallographic symmetry 1=cubic 2=hexagonal, needed in all orientation calcs
real(pReal), dimension (:,:,:), allocatable :: & real(pReal), dimension (:,:,:), allocatable :: &
crystallite_dt, & ! requested time increment of each grain crystallite_dt, & !< requested time increment of each grain
crystallite_subdt, & ! substepped time increment of each grain crystallite_subdt, & !< substepped time increment of each grain
crystallite_subFrac, & ! already calculated fraction of increment crystallite_subFrac, & !< already calculated fraction of increment
crystallite_subStep, & ! size of next integration step crystallite_subStep, & !< size of next integration step
crystallite_Temperature, & ! Temp of each grain crystallite_Temperature, & !< Temp of each grain
crystallite_partionedTemperature0, & ! Temp of each grain at start of homog inc crystallite_partionedTemperature0, & !< Temp of each grain at start of homog inc
crystallite_subTemperature0, & ! Temp of each grain at start of crystallite inc crystallite_subTemperature0, & !< Temp of each grain at start of crystallite inc
crystallite_dotTemperature ! evolution of Temperature of each grain crystallite_dotTemperature !< evolution of Temperature of each grain
real(pReal), dimension (:,:,:,:), allocatable :: & real(pReal), dimension (:,:,:,:), allocatable :: &
crystallite_Tstar_v, & ! current 2nd Piola-Kirchhoff stress vector (end of converged time step) crystallite_Tstar_v, & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step)
crystallite_Tstar0_v, & ! 2nd Piola-Kirchhoff stress vector at start of FE inc crystallite_Tstar0_v, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc
crystallite_partionedTstar0_v, & ! 2nd Piola-Kirchhoff stress vector at start of homog inc crystallite_partionedTstar0_v, & !< 2nd Piola-Kirchhoff stress vector at start of homog inc
crystallite_subTstar0_v, & ! 2nd Piola-Kirchhoff stress vector at start of crystallite inc crystallite_subTstar0_v, & !< 2nd Piola-Kirchhoff stress vector at start of crystallite inc
crystallite_orientation, & ! orientation as quaternion crystallite_orientation, & !< orientation as quaternion
crystallite_orientation0, & ! initial orientation as quaternion crystallite_orientation0, & !< initial orientation as quaternion
crystallite_rotation ! grain rotation away from initial orientation as axis-angle (in degrees) crystallite_rotation !< grain rotation away from initial orientation as axis-angle (in degrees)
real(pReal), dimension (:,:,:,:,:), allocatable :: & real(pReal), dimension (:,:,:,:,:), allocatable :: &
crystallite_Fe, & ! current "elastic" def grad (end of converged time step) crystallite_Fe, & !< current "elastic" def grad (end of converged time step)
crystallite_subFe0,& ! "elastic" def grad at start of crystallite inc crystallite_subFe0,& !< "elastic" def grad at start of crystallite inc
crystallite_Fp, & ! current plastic def grad (end of converged time step) crystallite_Fp, & !< current plastic def grad (end of converged time step)
crystallite_invFp, & ! inverse of current plastic def grad (end of converged time step) crystallite_invFp, & !< inverse of current plastic def grad (end of converged time step)
crystallite_Fp0, & ! plastic def grad at start of FE inc crystallite_Fp0, & !< plastic def grad at start of FE inc
crystallite_partionedFp0,& ! plastic def grad at start of homog inc crystallite_partionedFp0,& !< plastic def grad at start of homog inc
crystallite_subFp0,& ! plastic def grad at start of crystallite inc crystallite_subFp0,& !< plastic def grad at start of crystallite inc
crystallite_F0, & ! def grad at start of FE inc crystallite_F0, & !< def grad at start of FE inc
crystallite_partionedF, & ! def grad to be reached at end of homog inc crystallite_partionedF, & !< def grad to be reached at end of homog inc
crystallite_partionedF0, & ! def grad at start of homog inc crystallite_partionedF0, & !< def grad at start of homog inc
crystallite_subF, & ! def grad to be reached at end of crystallite inc crystallite_subF, & !< def grad to be reached at end of crystallite inc
crystallite_subF0, & ! def grad at start of crystallite inc crystallite_subF0, & !< def grad at start of crystallite inc
crystallite_Lp, & ! current plastic velocitiy grad (end of converged time step) crystallite_Lp, & !< current plastic velocitiy grad (end of converged time step)
crystallite_Lp0, & ! plastic velocitiy grad at start of FE inc crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc
crystallite_partionedLp0,& ! plastic velocity grad at start of homog inc crystallite_partionedLp0,& !< plastic velocity grad at start of homog inc
crystallite_subLp0,& ! plastic velocity grad at start of crystallite inc crystallite_subLp0,& !< plastic velocity grad at start of crystallite inc
crystallite_P, & ! 1st Piola-Kirchhoff stress per grain crystallite_P, & !< 1st Piola-Kirchhoff stress per grain
crystallite_disorientation ! disorientation between two neighboring ips (only calculated for single grain IPs) crystallite_disorientation !< disorientation between two neighboring ips (only calculated for single grain IPs)
real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: & real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: &
crystallite_dPdF, & ! current individual dPdF per grain (end of converged time step) crystallite_dPdF, & !< current individual dPdF per grain (end of converged time step)
crystallite_dPdF0, & ! individual dPdF per grain at start of FE inc crystallite_dPdF0, & !< individual dPdF per grain at start of FE inc
crystallite_partioneddPdF0, & ! individual dPdF per grain at start of homog inc crystallite_partioneddPdF0, & !< individual dPdF per grain at start of homog inc
crystallite_fallbackdPdF ! dPdF fallback for non-converged grains (elastic prediction) crystallite_fallbackdPdF !< dPdF fallback for non-converged grains (elastic prediction)
logical, dimension (:,:,:), allocatable :: & logical, dimension (:,:,:), allocatable :: &
crystallite_localPlasticity, & ! indicates this grain to have purely local constitutive law crystallite_localPlasticity, & !< indicates this grain to have purely local constitutive law
crystallite_requested, & ! flag to request crystallite calculation crystallite_requested, & !< flag to request crystallite calculation
crystallite_todo, & ! flag to indicate need for further computation crystallite_todo, & !< flag to indicate need for further computation
crystallite_converged ! convergence flag crystallite_converged !< convergence flag
CONTAINS contains
!******************************************************************** !********************************************************************
@ -149,8 +148,6 @@ integer(pInt), parameter :: myFile = 200_pInt, &
!*** input variables ***! !*** input variables ***!
real(pReal) Temperature real(pReal) Temperature
!*** output variables ***!
!*** local variables ***! !*** local variables ***!
integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) g, & ! grain number integer(pInt) g, & ! grain number
@ -1048,8 +1045,8 @@ use constitutive, only: constitutive_sizeDotState, &
implicit none implicit none
real(pReal), dimension(4), parameter :: timeStepFraction = (/0.5_pReal, 0.5_pReal, 1.0_pReal, 1.0_pReal/) ! factor giving the fraction of the original timestep used for Runge Kutta Integration real(pReal), dimension(4), parameter :: timeStepFraction = [0.5_pReal, 0.5_pReal, 1.0_pReal, 1.0_pReal] ! factor giving the fraction of the original timestep used for Runge Kutta Integration
real(pReal), dimension(4), parameter :: weight = (/1.0_pReal, 2.0_pReal, 2.0_pReal, 1.0_pReal/) ! weight of slope used for Runge Kutta integration real(pReal), dimension(4), parameter :: weight = [1.0_pReal, 2.0_pReal, 2.0_pReal, 1.0_pReal] ! weight of slope used for Runge Kutta integration
!*** input variables ***! !*** input variables ***!
integer(pInt), optional, intent(in):: ee, & ! element index integer(pInt), optional, intent(in):: ee, & ! element index
@ -3040,8 +3037,6 @@ LpLoop: do
#elif(FLOAT==4) #elif(FLOAT==4)
call sgetrf(9,9,inv_dR_dLp,9,ipiv,error) ! invert dR/dLp --> dLp/dR call sgetrf(9,9,inv_dR_dLp,9,ipiv,error) ! invert dR/dLp --> dLp/dR
call sgetri(9,inv_dR_dLp,9,ipiv,work,9,error) call sgetri(9,inv_dR_dLp,9,ipiv,work,9,error)
#else
NO SUITABLE PRECISION SELECTED, COMPILATION ABORTED
#endif #endif
if (error) then if (error) then
#ifndef _OPENMP #ifndef _OPENMP

View File

@ -903,15 +903,20 @@ function math_invSym3333(A)
real(pReal),dimension(3,3,3,3),intent(in) :: A real(pReal),dimension(3,3,3,3),intent(in) :: A
integer(pInt) :: ierr1, ierr2 integer(pInt) :: ierr
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
temp66_real = math_Mandel3333to66(A) temp66_real = math_Mandel3333to66(A)
call dgetrf(6,6,temp66_real,6,ipiv6,ierr1) #if(FLOAT==8)
call dgetri(6,temp66_real,6,ipiv6,work6,6,ierr2) call dgetrf(6,6,temp66_real,6,ipiv6,ierr)
if (ierr1*ierr2 == 0_pInt) then 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
math_invSym3333 = math_Mandel66to3333(temp66_real) math_invSym3333 = math_Mandel66to3333(temp66_real)
else else
call IO_error(400_pInt, ext_msg = 'math_invSym3333') call IO_error(400_pInt, ext_msg = 'math_invSym3333')
@ -945,8 +950,6 @@ subroutine math_invert(myDim,A, InvA, error)
#elif(FLOAT==4) #elif(FLOAT==4)
call sgetrf(myDim,myDim,invA,myDim,ipiv,ierr) call sgetrf(myDim,myDim,invA,myDim,ipiv,ierr)
call sgetri(myDim,InvA,myDim,ipiv,work,myDim,ierr) call sgetri(myDim,InvA,myDim,ipiv,work,myDim,ierr)
#else
NO SUITABLE PRECISION SELECTED, COMPILATION ABORTED
#endif #endif
if (ierr == 0_pInt) then if (ierr == 0_pInt) then
error = .false. error = .false.
@ -2770,7 +2773,6 @@ function math_curlFFT(geomdim,field)
wgt = 1.0_pReal/real(res(1)*res(2)*res(3),pReal) wgt = 1.0_pReal/real(res(1)*res(2)*res(3),pReal)
res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt)
call fftw_set_timelimit(fftw_timelimit) call fftw_set_timelimit(fftw_timelimit)
field_fftw = fftw_alloc_complex(int(res1_red *res(2)*res(3)*vec_tens*3_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8) field_fftw = fftw_alloc_complex(int(res1_red *res(2)*res(3)*vec_tens*3_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8)
call c_f_pointer(field_fftw, field_real, [res(1)+2_pInt,res(2),res(3),vec_tens,3_pInt]) call c_f_pointer(field_fftw, field_real, [res(1)+2_pInt,res(2),res(3),vec_tens,3_pInt])
@ -2890,7 +2892,6 @@ function math_divergenceFFT(geomdim,field)
res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
wgt = 1.0_pReal/real(res(1)*res(2)*res(3),pReal) wgt = 1.0_pReal/real(res(1)*res(2)*res(3),pReal)
if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt)
call fftw_set_timelimit(fftw_timelimit) call fftw_set_timelimit(fftw_timelimit)
field_fftw = fftw_alloc_complex(int(res1_red*res(2)*res(3)*vec_tens*3_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8) field_fftw = fftw_alloc_complex(int(res1_red*res(2)*res(3)*vec_tens*3_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8)
call c_f_pointer(field_fftw, field_real, [res(1)+2_pInt,res(2),res(3),vec_tens,3_pInt]) call c_f_pointer(field_fftw, field_real, [res(1)+2_pInt,res(2),res(3),vec_tens,3_pInt])

View File

@ -38,6 +38,9 @@ module prec
private private
#if (FLOAT==4) #if (FLOAT==4)
#ifdef Spectral
SPECTRAL SOLVER DOES NOT SUPPORT SINGLE PRECISION, STOPING 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) 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 LEGACY_COMPILER #ifdef LEGACY_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) 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)
@ -51,12 +54,16 @@ module prec
#else #else
real(pReal), parameter, public :: DAMASK_NaN = real(Z'7FF8000000000000', pReal) !< 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 = real(Z'7FF8000000000000', pReal) !< 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)
#endif #endif
#else
NO SUITABLE PRECISION SELECTED, STOPING COMPILATION
#endif #endif
#if (INT==4) #if (INT==4)
integer, parameter, public :: pInt = 4 !< integer representation 32 bit (was selected_int_kind(9), number with at least up to +- 1e9) integer, parameter, public :: pInt = 4 !< integer representation 32 bit (was selected_int_kind(9), number with at least up to +- 1e9)
#elif (INT==8) #elif (INT==8)
integer, parameter, public :: pInt = 8 !< integer representation 64 bit (was selected_int_kind(12), number with at least up to +- 1e12) integer, parameter, public :: pInt = 8 !< integer representation 64 bit (was selected_int_kind(12), number with at least up to +- 1e12)
#else
NO SUITABLE PRECISION SELECTED, STOPING COMPILATION
#endif #endif
integer, parameter, public :: pLongInt = 8 !< integer representation 64 bit (was selected_int_kind(12), number with at least up to +- 1e12) integer, parameter, public :: pLongInt = 8 !< integer representation 64 bit (was selected_int_kind(12), number with at least up to +- 1e12)