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

Merge commit 'v2.0.1-741-gea87e8e'

This commit is contained in:
Test User 2017-05-17 21:33:40 +02:00
parent a774fffc24 ea87e8e53c
commit 85917b05ba
7 changed files with 92 additions and 139 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

120
src/crystallite.f90
View File

@ -1232,8 +1232,8 @@ end subroutine crystallite_stressAndItsTangent
!> @brief integrate stress, state with 4th order explicit Runge Kutta method !> @brief integrate stress, state with 4th order explicit Runge Kutta method
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateRK4() subroutine crystallite_integrateStateRK4()
use prec, only: & use, intrinsic :: &
prec_isNaN IEEE_arithmetic
use numerics, only: & use numerics, only: &
numerics_integrationMode numerics_integrationMode
use debug, only: & use debug, only: &
@ -1316,7 +1316,7 @@ subroutine crystallite_integrateStateRK4()
! first Runge-Kutta step ! first Runge-Kutta step
!$OMP PARALLEL !$OMP PARALLEL
!$OMP DO !$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) & if (crystallite_todo(g,i,e)) &
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), & call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), &
crystallite_Fe, & crystallite_Fe, &
@ -1326,22 +1326,22 @@ subroutine crystallite_integrateStateRK4()
!$OMP ENDDO !$OMP ENDDO
!$OMP DO PRIVATE(p,c,NaN) !$OMP DO PRIVATE(p,c,NaN)
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo) !$OMP FLUSH(crystallite_todo)
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
c = phasememberAt(g,i,e) c = phasememberAt(g,i,e)
p = phaseAt(g,i,e) p = phaseAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,c))) NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p) do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c))) NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo enddo
if (NaN) then ! NaN occured in any dotState if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo) !$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo) !$OMP END CRITICAL (checkTodo)
else ! if broken local... else ! if broken local...
crystallite_todo(g,i,e) = .false. ! ... skip this one next time crystallite_todo(g,i,e) = .false. ! ... skip this one next time
endif endif
endif endif
endif endif
@ -1356,7 +1356,7 @@ subroutine crystallite_integrateStateRK4()
!$OMP PARALLEL !$OMP PARALLEL
!$OMP DO PRIVATE(p,c) !$OMP DO PRIVATE(p,c)
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
p = phaseAt(g,i,e) p = phaseAt(g,i,e)
c = phasememberAt(g,i,e) c = phasememberAt(g,i,e)
@ -1371,7 +1371,7 @@ subroutine crystallite_integrateStateRK4()
!$OMP ENDDO !$OMP ENDDO
!$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,c) !$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,c)
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
p = phaseAt(g,i,e) p = phaseAt(g,i,e)
@ -1475,9 +1475,9 @@ subroutine crystallite_integrateStateRK4()
p = phaseAt(g,i,e) p = phaseAt(g,i,e)
c = phasememberAt(g,i,e) c = phasememberAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,c))) NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p) do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c))) NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo enddo
if (NaN) then ! NaN occured in any dotState if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
@ -1528,8 +1528,8 @@ end subroutine crystallite_integrateStateRK4
!> adaptive step size (use 5th order solution to advance = "local extrapolation") !> adaptive step size (use 5th order solution to advance = "local extrapolation")
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateRKCK45() subroutine crystallite_integrateStateRKCK45()
use prec, only: & use, intrinsic :: &
prec_isNaN IEEE_arithmetic
use debug, only: & use debug, only: &
debug_level, & debug_level, &
debug_crystallite, & debug_crystallite, &
@ -1647,9 +1647,9 @@ subroutine crystallite_integrateStateRKCK45()
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
cc = phasememberAt(g,i,e) cc = phasememberAt(g,i,e)
p = phaseAt(g,i,e) p = phaseAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,cc))) NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,cc)))
do mySource = 1_pInt, phase_Nsources(p) do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,cc))) NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,cc)))
enddo enddo
if (NaN) then ! NaN occured in any dotState if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
@ -1801,9 +1801,9 @@ subroutine crystallite_integrateStateRKCK45()
p = phaseAt(g,i,e) p = phaseAt(g,i,e)
cc = phasememberAt(g,i,e) cc = phasememberAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,cc))) NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,cc)))
do mySource = 1_pInt, phase_Nsources(p) do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,cc))) NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,cc)))
enddo enddo
if (NaN) then ! NaN occured in any dotState if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
@ -2031,8 +2031,8 @@ end subroutine crystallite_integrateStateRKCK45
!> @brief integrate stress, state with 1st order Euler method with adaptive step size !> @brief integrate stress, state with 1st order Euler method with adaptive step size
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateAdaptiveEuler() subroutine crystallite_integrateStateAdaptiveEuler()
use prec, only: & use, intrinsic :: &
prec_isNaN IEEE_arithmetic
use debug, only: & use debug, only: &
debug_level, & debug_level, &
debug_crystallite, & debug_crystallite, &
@ -2133,9 +2133,9 @@ subroutine crystallite_integrateStateAdaptiveEuler()
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
p = phaseAt(g,i,e) p = phaseAt(g,i,e)
c = phasememberAt(g,i,e) c = phasememberAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,c))) NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p) do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c))) NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo enddo
if (NaN) then ! NaN occured in any dotState if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
@ -2254,9 +2254,9 @@ subroutine crystallite_integrateStateAdaptiveEuler()
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
p = phaseAt(g,i,e) p = phaseAt(g,i,e)
c = phasememberAt(g,i,e) c = phasememberAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,c))) NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p) do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c))) NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo enddo
if (NaN) then ! NaN occured in any dotState if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
@ -2391,8 +2391,8 @@ end subroutine crystallite_integrateStateAdaptiveEuler
!> @brief integrate stress, and state with 1st order explicit Euler method !> @brief integrate stress, and state with 1st order explicit Euler method
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateEuler() subroutine crystallite_integrateStateEuler()
use prec, only: & use, intrinsic :: &
prec_isNaN IEEE_arithmetic
use debug, only: & use debug, only: &
debug_level, & debug_level, &
debug_crystallite, & debug_crystallite, &
@ -2471,9 +2471,9 @@ eIter = FEsolving_execElem(1:2)
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
c = phasememberAt(g,i,e) c = phasememberAt(g,i,e)
p = phaseAt(g,i,e) p = phaseAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,c))) NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p) do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c))) NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo enddo
if (NaN) then ! NaN occured in any dotState if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e) .and. .not. numerics_timeSyncing) then ! if broken non-local... if (.not. crystallite_localPlasticity(g,i,e) .and. .not. numerics_timeSyncing) then ! if broken non-local...
@ -2614,8 +2614,8 @@ end subroutine crystallite_integrateStateEuler
!> using Fixed Point Iteration to adapt the stepsize !> using Fixed Point Iteration to adapt the stepsize
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateFPI() subroutine crystallite_integrateStateFPI()
use prec, only: & use, intrinsic :: &
prec_isNaN IEEE_arithmetic
use debug, only: & use debug, only: &
debug_e, & debug_e, &
debug_i, & debug_i, &
@ -2732,22 +2732,22 @@ subroutine crystallite_integrateStateFPI()
!$OMP ENDDO !$OMP ENDDO
!$OMP DO PRIVATE(p,c,NaN) !$OMP DO PRIVATE(p,c,NaN)
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo) !$OMP FLUSH(crystallite_todo)
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
p = phaseAt(g,i,e) p = phaseAt(g,i,e)
c = phasememberAt(g,i,e) c = phasememberAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,c))) NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p) do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c))) NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo enddo
if (NaN) then ! NaN occured in any dotState if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a non-local... if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a non-local...
!$OMP CRITICAL (checkTodo) !$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals done (and broken) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals done (and broken)
!$OMP END CRITICAL (checkTodo) !$OMP END CRITICAL (checkTodo)
else ! broken one was local... else ! broken one was local...
crystallite_todo(g,i,e) = .false. ! ... done (and broken) crystallite_todo(g,i,e) = .false. ! ... done (and broken)
endif endif
endif endif
endif endif
@ -2851,9 +2851,9 @@ subroutine crystallite_integrateStateFPI()
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
p = phaseAt(g,i,e) p = phaseAt(g,i,e)
c = phasememberAt(g,i,e) c = phasememberAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,c))) NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p) do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c))) NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo enddo
if (NaN) then ! NaN occured in any dotState if (NaN) then ! NaN occured in any dotState
crystallite_todo(g,i,e) = .false. ! ... skip me next time crystallite_todo(g,i,e) = .false. ! ... skip me next time
@ -3062,8 +3062,9 @@ end subroutine crystallite_integrateStateFPI
!> returns true, if state jump was successfull or not needed. false indicates NaN in delta state !> returns true, if state jump was successfull or not needed. false indicates NaN in delta state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
logical function crystallite_stateJump(ipc,ip,el) logical function crystallite_stateJump(ipc,ip,el)
use, intrinsic :: &
IEEE_arithmetic
use prec, only: & use prec, only: &
prec_isNaN, &
dNeq0 dNeq0
use debug, only: & use debug, only: &
debug_level, & debug_level, &
@ -3098,7 +3099,7 @@ logical function crystallite_stateJump(ipc,ip,el)
p = phaseAt(ipc,ip,el) p = phaseAt(ipc,ip,el)
call constitutive_collectDeltaState(crystallite_Tstar_v(1:6,ipc,ip,el), crystallite_Fe(1:3,1:3,ipc,ip,el), ipc,ip,el) call constitutive_collectDeltaState(crystallite_Tstar_v(1:6,ipc,ip,el), crystallite_Fe(1:3,1:3,ipc,ip,el), ipc,ip,el)
mySizePlasticDeltaState = plasticState(p)%sizeDeltaState mySizePlasticDeltaState = plasticState(p)%sizeDeltaState
if( any(prec_isNaN(plasticState(p)%deltaState(:,c)))) then ! NaN occured in deltaState if( any(IEEE_is_NaN(plasticState(p)%deltaState(:,c)))) then ! NaN occured in deltaState
crystallite_stateJump = .false. crystallite_stateJump = .false.
return return
endif endif
@ -3106,7 +3107,7 @@ logical function crystallite_stateJump(ipc,ip,el)
plasticState(p)%deltaState(1:mySizePlasticDeltaState,c) plasticState(p)%deltaState(1:mySizePlasticDeltaState,c)
do mySource = 1_pInt, phase_Nsources(p) do mySource = 1_pInt, phase_Nsources(p)
mySizeSourceDeltaState = sourceState(p)%p(mySource)%sizeDeltaState mySizeSourceDeltaState = sourceState(p)%p(mySource)%sizeDeltaState
if( any(prec_isNaN(sourceState(p)%p(mySource)%deltaState(:,c)))) then ! NaN occured in deltaState if( any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(:,c)))) then ! NaN occured in deltaState
crystallite_stateJump = .false. crystallite_stateJump = .false.
return return
endif endif
@ -3169,9 +3170,10 @@ logical function crystallite_integrateStress(&
el,& ! element number el,& ! element number
timeFraction & timeFraction &
) )
use, intrinsic :: &
IEEE_arithmetic
use prec, only: pLongInt, & use prec, only: pLongInt, &
tol_math_check, & tol_math_check, &
prec_isNaN, &
dEq0 dEq0
use numerics, only: nStress, & use numerics, only: nStress, &
aTol_crystalliteStress, & aTol_crystalliteStress, &
@ -3426,11 +3428,11 @@ logical function crystallite_integrateStress(&
!* update current residuum and check for convergence of loop !* update current residuum and check for convergence of loop
aTolLp = max(rTol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error aTolLp = max(rTol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error
aTol_crystalliteStress) ! minimum lower cutoff aTol_crystalliteStress) ! minimum lower cutoff
residuumLp = Lpguess - Lp_constitutive residuumLp = Lpguess - Lp_constitutive
if (any(prec_isNaN(residuumLp))) then ! NaN in residuum... if (any(IEEE_is_NaN(residuumLp))) then ! NaN in residuum...
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3,a)') '<< CRYST >> integrateStress encountered NaN at el (elFE) ip ipc ', & write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3,a)') '<< CRYST >> integrateStress encountered NaN at el (elFE) ip ipc ', &
@ -3438,16 +3440,16 @@ logical function crystallite_integrateStress(&
' ; iteration ', NiterationStressLp,& ' ; iteration ', NiterationStressLp,&
' >> returning..!' ' >> returning..!'
#endif #endif
return ! ...me = .false. to inform integrator about problem return ! ...me = .false. to inform integrator about problem
elseif (norm2(residuumLp) < aTolLp) then ! converged if below absolute tolerance elseif (norm2(residuumLp) < aTolLp) then ! converged if below absolute tolerance
exit LpLoop ! ...leave iteration loop exit LpLoop ! ...leave iteration loop
elseif ( NiterationStressLp == 1_pInt & elseif ( NiterationStressLp == 1_pInt &
.or. norm2(residuumLp) < norm2(residuumLp_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)... .or. norm2(residuumLp) < norm2(residuumLp_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
residuumLp_old = residuumLp ! ...remember old values and... residuumLp_old = residuumLp ! ...remember old values and...
Lpguess_old = Lpguess Lpguess_old = Lpguess
steplengthLp = steplengthLp0 ! ...proceed with normal step length (calculate new search direction) steplengthLp = steplengthLp0 ! ...proceed with normal step length (calculate new search direction)
else ! not converged and residuum not improved... else ! not converged and residuum not improved...
steplengthLp = 0.5_pReal * steplengthLp ! ...try with smaller step length in same direction steplengthLp = 0.5_pReal * steplengthLp ! ...try with smaller step length in same direction
Lpguess = Lpguess_old + steplengthLp * deltaLp Lpguess = Lpguess_old + steplengthLp * deltaLp
cycle LpLoop cycle LpLoop
endif endif
@ -3520,7 +3522,7 @@ logical function crystallite_integrateStress(&
aTolLi = max(rTol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error aTolLi = max(rTol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error
aTol_crystalliteStress) ! minimum lower cutoff aTol_crystalliteStress) ! minimum lower cutoff
residuumLi = Liguess - Li_constitutive residuumLi = Liguess - Li_constitutive
if (any(prec_isNaN(residuumLi))) then ! NaN in residuum... if (any(IEEE_is_NaN(residuumLi))) then ! NaN in residuum...
return ! ...me = .false. to inform integrator about problem return ! ...me = .false. to inform integrator about problem
elseif (norm2(residuumLi) < aTolLi) then ! converged if below absolute tolerance elseif (norm2(residuumLi) < aTolLi) then ! converged if below absolute tolerance
exit LiLoop ! ...leave iteration loop exit LiLoop ! ...leave iteration loop

6
src/lattice.f90
View File

@ -2170,8 +2170,8 @@ pure function lattice_symmetrize33(struct,T33)
!> @brief figures whether unit quat falls into stereographic standard triangle !> @brief figures whether unit quat falls into stereographic standard triangle
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
logical pure function lattice_qInSST(Q, struct) logical pure function lattice_qInSST(Q, struct)
use prec, only: & use, intrinsic :: &
prec_isNaN IEEE_arithmetic
use math, only: & use math, only: &
math_qToRodrig math_qToRodrig
@ -2181,7 +2181,7 @@ logical pure function lattice_qInSST(Q, struct)
real(pReal), dimension(3) :: Rodrig ! Rodrigues vector of Q real(pReal), dimension(3) :: Rodrig ! Rodrigues vector of Q
Rodrig = math_qToRodrig(Q) Rodrig = math_qToRodrig(Q)
if (any(prec_isNaN(Rodrig))) then if (any(IEEE_is_NaN(Rodrig))) then
lattice_qInSST = .false. lattice_qInSST = .false.
else else
select case (struct) select case (struct)

36
src/math.f90
View File

@ -1666,15 +1666,16 @@ end function math_qToEulerAxisAngle
!> @brief Rodrigues vector (x, y, z) from unit quaternion (w+ix+jy+kz) !> @brief Rodrigues vector (x, y, z) from unit quaternion (w+ix+jy+kz)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_qToRodrig(Q) pure function math_qToRodrig(Q)
use, intrinsic :: &
IEEE_arithmetic
use prec, only: & use prec, only: &
DAMASK_NaN, &
tol_math_check tol_math_check
implicit none implicit none
real(pReal), dimension(4), intent(in) :: Q real(pReal), dimension(4), intent(in) :: Q
real(pReal), dimension(3) :: math_qToRodrig real(pReal), dimension(3) :: math_qToRodrig
math_qToRodrig = merge(Q(2:4)/Q(1),DAMASK_NaN,abs(Q(1)) > tol_math_check) ! NaN for 180 deg since Rodrig is unbound math_qToRodrig = merge(Q(2:4)/Q(1),IEEE_value(1.0_pReal,IEEE_quiet_NaN),abs(Q(1)) > tol_math_check)! NaN for 180 deg since Rodrig is unbound
end function math_qToRodrig end function math_qToRodrig
@ -2092,11 +2093,11 @@ end function math_rotationalPart33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Eigenvalues of symmetric matrix m !> @brief Eigenvalues of symmetric matrix m
! will return NaN on error ! will return NaN on error
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function math_eigenvaluesSym(m) function math_eigenvaluesSym(m)
use prec, only: & use, intrinsic :: &
DAMASK_NaN IEEE_arithmetic
implicit none implicit none
real(pReal), dimension(:,:), intent(in) :: m real(pReal), dimension(:,:), intent(in) :: m
@ -2109,7 +2110,7 @@ function math_eigenvaluesSym(m)
vectors = m ! copy matrix to input (doubles as output) array vectors = m ! copy matrix to input (doubles as output) array
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)
if (info /= 0_pInt) math_eigenvaluesSym = DAMASK_NaN if (info /= 0_pInt) math_eigenvaluesSym = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
end function math_eigenvaluesSym end function math_eigenvaluesSym
@ -2701,29 +2702,13 @@ pure function math_rotate_forward3333(tensor,rot_tensor)
end function math_rotate_forward3333 end function math_rotate_forward3333
!--------------------------------------------------------------------------------------------------
!> @brief calculate average of tensor field
!--------------------------------------------------------------------------------------------------
function math_tensorAvg(field)
implicit none
real(pReal), dimension(3,3) :: math_tensorAvg
real(pReal), intent(in), dimension(:,:,:,:,:) :: field
real(pReal) :: wgt
wgt = 1.0_pReal/real(size(field,3)*size(field,4)*size(field,5), pReal)
math_tensorAvg = sum(sum(sum(field,dim=5),dim=4),dim=3)*wgt
end function math_tensorAvg
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief limits a scalar value to a certain range (either one or two sided) !> @brief limits a scalar value to a certain range (either one or two sided)
! Will return NaN if left > right ! Will return NaN if left > right
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
real(pReal) pure function math_limit(a, left, right) real(pReal) pure function math_limit(a, left, right)
use prec, only: & use, intrinsic :: &
DAMASK_NaN IEEE_arithmetic
implicit none implicit none
real(pReal), intent(in) :: a real(pReal), intent(in) :: a
@ -2735,7 +2720,8 @@ real(pReal) pure function math_limit(a, left, right)
merge(right, huge(a), present(right)) & merge(right, huge(a), present(right)) &
) )
if (present(left) .and. present(right)) math_limit = merge (DAMASK_NaN,math_limit, left>right) if (present(left) .and. present(right)) &
math_limit = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_limit, left>right)
end function math_limit end function math_limit

17
src/plastic_isotropic.f90
View File

@ -9,8 +9,7 @@
module plastic_isotropic module plastic_isotropic
use prec, only: & use prec, only: &
pReal,& pReal,&
pInt, & pInt
DAMASK_NaN
implicit none implicit none
private private
@ -36,14 +35,14 @@ module plastic_isotropic
integer(kind(undefined_ID)), allocatable, dimension(:) :: & integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID outputID
real(pReal) :: & real(pReal) :: &
fTaylor = DAMASK_NaN, & fTaylor, &
tau0 = DAMASK_NaN, & tau0, &
gdot0 = DAMASK_NaN, & gdot0, &
n = DAMASK_NaN, & n, &
h0 = DAMASK_NaN, & h0, &
h0_slopeLnRate = 0.0_pReal, & h0_slopeLnRate = 0.0_pReal, &
tausat = DAMASK_NaN, & tausat, &
a = DAMASK_NaN, & a, &
aTolFlowstress = 1.0_pReal, & aTolFlowstress = 1.0_pReal, &
aTolShear = 1.0e-6_pReal, & aTolShear = 1.0e-6_pReal, &
tausat_SinhFitA= 0.0_pReal, & tausat_SinhFitA= 0.0_pReal, &

10
src/plastic_nonlocal.f90
View File

@ -2377,9 +2377,9 @@ end subroutine plastic_nonlocal_deltaState
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
subroutine plastic_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, & subroutine plastic_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, &
timestep,subfrac, ip,el) timestep,subfrac, ip,el)
use, intrinsic :: &
use prec, only: DAMASK_NaN, & IEEE_arithmetic
dNeq0, & use prec, only: dNeq0, &
dNeq, & dNeq, &
dEq0 dEq0
use numerics, only: numerics_integrationMode, & use numerics, only: numerics_integrationMode, &
@ -2701,7 +2701,7 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then
write(6,'(a)') '<< CONST >> enforcing cutback !!!' write(6,'(a)') '<< CONST >> enforcing cutback !!!'
endif endif
#endif #endif
plasticState(p)%dotState = DAMASK_NaN ! -> return NaN and, hence, enforce cutback plasticState(p)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) ! -> return NaN and, hence, enforce cutback
return return
endif endif
@ -2984,7 +2984,7 @@ if ( any(rhoSglOriginal(1:ns,1:4) + rhoDot(1:ns,1:4) * timestep < -aTolRho(in
write(6,'(a)') '<< CONST >> enforcing cutback !!!' write(6,'(a)') '<< CONST >> enforcing cutback !!!'
endif endif
#endif #endif
plasticState(p)%dotState = DAMASK_NaN plasticState(p)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
return return
else else
forall (s = 1:ns, t = 1_pInt:4_pInt) forall (s = 1:ns, t = 1_pInt:4_pInt)

36
src/prec.f90
View File

@ -5,22 +5,12 @@
!> @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 !> @brief setting precision for real and int type
!> @details setting precision for real and int type and for DAMASK_NaN. Definition is made
!! depending on makro "INT" defined during compilation
!! for details on NaN see https://software.intel.com/en-us/forums/topic/294680
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module prec module prec
#if !(defined(__GFORTRAN__) && __GNUC__ < 5)
use, intrinsic :: & ! unfortunately not avialable in gfortran <= 5
IEEE_arithmetic
#endif
implicit none implicit none
private private
#if (FLOAT==8) #if (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)
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)
#else #else
NO SUITABLE PRECISION FOR REAL SELECTED, STOPPING COMPILATION NO SUITABLE PRECISION FOR REAL SELECTED, STOPPING COMPILATION
#endif #endif
@ -106,7 +96,6 @@ module prec
public :: & public :: &
prec_init, & prec_init, &
prec_isNaN, &
dEq, & dEq, &
dEq0, & dEq0, &
cEq, & cEq, &
@ -118,7 +107,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief reporting precision and checking if DAMASK_NaN is set correctly !> @brief reporting precision
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine prec_init subroutine prec_init
use, intrinsic :: & use, intrinsic :: &
@ -133,36 +122,13 @@ subroutine prec_init
write(6,'(a,i3)') ' Bytes for pReal: ',pReal write(6,'(a,i3)') ' Bytes for pReal: ',pReal
write(6,'(a,i3)') ' Bytes for pInt: ',pInt write(6,'(a,i3)') ' Bytes for pInt: ',pInt
write(6,'(a,i3)') ' Bytes for pLongInt: ',pLongInt write(6,'(a,i3)') ' Bytes for pLongInt: ',pLongInt
write(6,'(a,e10.3)') ' NaN: ', DAMASK_NaN
write(6,'(a,l3)') ' NaN != NaN: ',DAMASK_NaN /= DAMASK_NaN
write(6,'(a,l3,/)') ' NaN check passed ',prec_isNAN(DAMASK_NaN)
if ((.not. prec_isNaN(DAMASK_NaN)) .or. (DAMASK_NaN == DAMASK_NaN)) call quit(9000)
realloc_lhs_test = [1_pInt,2_pInt] realloc_lhs_test = [1_pInt,2_pInt]
if (realloc_lhs_test(2)/=2_pInt) call quit(9000) if (realloc_lhs_test(2)/=2_pInt) call quit(9000)
end subroutine prec_init end subroutine prec_init
!--------------------------------------------------------------------------------------------------
!> @brief figures out if a floating point number is NaN
! basically just a small wrapper, because gfortran < 5.0 does not have the IEEE module
!--------------------------------------------------------------------------------------------------
logical elemental pure function prec_isNaN(a)
implicit none
real(pReal), intent(in) :: a
#if (defined(__GFORTRAN__) && __GNUC__ < 5)
intrinsic :: isNaN
prec_isNaN = isNaN(a)
#else
prec_isNaN = IEEE_is_NaN(a)
#endif
end function prec_isNaN
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief equality comparison for float with double precision !> @brief equality comparison for float with double precision
! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq ! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq

6
src/spectral_utilities.f90
View File

@ -741,8 +741,8 @@ end function utilities_curlRMS
!> @brief calculates mask compliance tensor used to adjust F to fullfill stress BC !> @brief calculates mask compliance tensor used to adjust F to fullfill stress BC
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function utilities_maskedCompliance(rot_BC,mask_stress,C) function utilities_maskedCompliance(rot_BC,mask_stress,C)
use prec, only: & use, intrinsic :: &
prec_isNaN IEEE_arithmetic
use IO, only: & use IO, only: &
IO_error IO_error
use math, only: & use math, only: &
@ -794,7 +794,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
endif; enddo; endif; enddo endif; enddo; endif; enddo
call math_invert(size_reduced, c_reduced, s_reduced, errmatinv) ! invert reduced stiffness call math_invert(size_reduced, c_reduced, s_reduced, errmatinv) ! invert reduced stiffness
if (any(prec_isNaN(s_reduced))) errmatinv = .true. if (any(IEEE_is_NaN(s_reduced))) errmatinv = .true.
if(errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance') if(errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance')
temp99_Real = 0.0_pReal ! fill up compliance with zeros temp99_Real = 0.0_pReal ! fill up compliance with zeros
k = 0_pInt k = 0_pInt