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

fixed floating point exception in nonlocal (calculation of dUpper at several places)

This commit is contained in:
Martin Diehl 2014-06-26 20:01:38 +00:00
parent 4c7d6c4464
commit 9be7c0d5f3
3 changed files with 34 additions and 38 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
code/constitutive.f90
View File

@ -1020,7 +1020,7 @@ logical function constitutive_collectDeltaState(Tstar_v, ipc, ip, el)
#ifdef NEWSTATE #ifdef NEWSTATE
call constitutive_nonlocal_deltaState(Tstar_v,ip,el) call constitutive_nonlocal_deltaState(Tstar_v,ip,el)
#else #else
call constitutive_nonlocal_deltaState(constitutive_deltaState(ipc,ip,el),& call constitutive_nonlocal_deltaState(constitutive_deltaState(ipc,ip,el)%p,&
constitutive_state(ipc,ip,el), Tstar_v,ipc,ip,el) constitutive_state(ipc,ip,el), Tstar_v,ipc,ip,el)
#endif #endif
case default case default

54
code/constitutive_nonlocal.f90
View File

@ -2335,7 +2335,7 @@ real(pReal), dimension(6), intent(in) :: Tstar_v ! current
#ifndef NEWSTATE #ifndef NEWSTATE
type(p_vec), intent(inout) :: & type(p_vec), intent(inout) :: &
state ! current microstructural state state ! current microstructural state
type(p_vec), intent(out) :: deltaState ! change of state variables / microstructure real(pReal), dimension(:), intent(out) :: deltaState ! change of state variables / microstructure
#else #else
integer(pInt) :: & integer(pInt) :: &
p, & !< phase p, & !< phase
@ -2370,15 +2370,11 @@ real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,e
dUpperOld, & ! old maximum stable dipole distance for edges and screws dUpperOld, & ! old maximum stable dipole distance for edges and screws
deltaDUpper ! change in maximum stable dipole distance for edges and screws deltaDUpper ! change in maximum stable dipole distance for edges and screws
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt &
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == 1)& .and. ((debug_e == el .and. debug_i == ip .and. debug_g == 1)&
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) &
write(6,*) write(6,'(/,a,i8,1x,i2,1x,i1,/)') '<< CONST >> nonlocal_deltaState at el ip ipc ',el,ip,1
write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_deltaState at el ip ipc ',el,ip,1
write(6,*)
endif
#endif #endif
phase = material_phase(1,ip,el) phase = material_phase(1,ip,el)
@ -2457,19 +2453,21 @@ dUpper(1:ns,1) = lattice_mu(phase) * burgers(1:ns,instance) &
dUpper(1:ns,2) = lattice_mu(phase) * burgers(1:ns,instance) / (4.0_pReal * pi * abs(tau)) dUpper(1:ns,2) = lattice_mu(phase) * burgers(1:ns,instance) / (4.0_pReal * pi * abs(tau))
!in the test there is an FPE exception. Division by zero? forall (c = 1_pInt:2_pInt)
forall (c = 1_pInt:2_pInt) & where(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)) >= tiny(0.0_pReal)) &
dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) &
+ abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), & + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
dUpper(1:ns,c)) dUpper(1:ns,c))
end forall
dUpper = max(dUpper,dLower) dUpper = max(dUpper,dLower)
deltaDUpper = dUpper - dUpperOld deltaDUpper = dUpper - dUpperOld
!*** dissociation by stress increase !*** dissociation by stress increase
deltaRhoDipole2SingleStress = 0.0_pReal deltaRhoDipole2SingleStress = 0.0_pReal
forall (c=1_pInt:2_pInt, s=1_pInt:ns, deltaDUpper(s,c) < 0.0_pReal) & forall (c=1_pInt:2_pInt, s=1_pInt:ns, deltaDUpper(s,c) < 0.0_pReal .and. &
(dUpperOld(s,c) - dLower(s,c))/= 0.0_pReal) &
deltaRhoDipole2SingleStress(s,8_pInt+c) = rhoDip(s,c) * deltaDUpper(s,c) & deltaRhoDipole2SingleStress(s,8_pInt+c) = rhoDip(s,c) * deltaDUpper(s,c) &
/ (dUpperOld(s,c) - dLower(s,c)) / (dUpperOld(s,c) - dLower(s,c))
@ -2478,7 +2476,6 @@ forall (t=1_pInt:4_pInt) &
* deltaRhoDipole2SingleStress(1_pInt:ns,(t-1_pInt)/2_pInt+9_pInt) * deltaRhoDipole2SingleStress(1_pInt:ns,(t-1_pInt)/2_pInt+9_pInt)
!*** store new maximum dipole height in state !*** store new maximum dipole height in state
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) &
@ -2504,16 +2501,15 @@ endforall
forall (s = 1:ns, c = 1_pInt:2_pInt) & forall (s = 1:ns, c = 1_pInt:2_pInt) &
plasticState(p)%deltaState(iRhoD(s,c,instance),o) = deltaRho(s,c+8_pInt) plasticState(p)%deltaState(iRhoD(s,c,instance),o) = deltaRho(s,c+8_pInt)
#else #else
deltaState%p = 0.0_pReal deltaState = 0.0_pReal
forall (s = 1:ns, t = 1_pInt:4_pInt) forall (s = 1:ns, t = 1_pInt:4_pInt)
deltaState%p(iRhoU(s,t,instance)) = deltaRho(s,t) deltaState(iRhoU(s,t,instance)) = deltaRho(s,t)
deltaState%p(iRhoB(s,t,instance)) = deltaRho(s,t+4_pInt) deltaState(iRhoB(s,t,instance)) = deltaRho(s,t+4_pInt)
endforall endforall
forall (s = 1:ns, c = 1_pInt:2_pInt) & forall (s = 1:ns, c = 1_pInt:2_pInt) &
deltaState%p(iRhoD(s,c,instance)) = deltaRho(s,c+8_pInt) deltaState(iRhoD(s,c,instance)) = deltaRho(s,c+8_pInt)
#endif #endif
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)&
@ -2698,11 +2694,8 @@ allocate(constitutive_nonlocal_dotState(plasticState(p)%sizeDotState), source =
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt &
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)&
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) &
write(6,*) write(6,'(/,a,i8,1x,i2,1x,i1,/)') '<< CONST >> nonlocal_dotState at el ip ipc ',el,ip,ipc
write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_dotState at el ip ipc ',el,ip,ipc
write(6,*)
endif
#endif #endif
phase = material_phase(ipc,ip,el) phase = material_phase(ipc,ip,el)
@ -2811,14 +2804,15 @@ dUpper(1:ns,1) = lattice_mu(phase) * burgers(1:ns,instance) &
/ (8.0_pReal * pi * (1.0_pReal - lattice_nu(phase)) * abs(tau)) / (8.0_pReal * pi * (1.0_pReal - lattice_nu(phase)) * abs(tau))
dUpper(1:ns,2) = lattice_mu(phase) * burgers(1:ns,instance) & dUpper(1:ns,2) = lattice_mu(phase) * burgers(1:ns,instance) &
/ (4.0_pReal * pi * abs(tau)) / (4.0_pReal * pi * abs(tau))
forall (c = 1_pInt:2_pInt) & forall (c = 1_pInt:2_pInt)
where(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)) >= tiny(0.0_pReal)) &
dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) &
+ abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), & + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
dUpper(1:ns,c)) dUpper(1:ns,c))
end forall
dUpper = max(dUpper,dLower) dUpper = max(dUpper,dLower)
!**************************************************************************** !****************************************************************************
!*** calculate dislocation multiplication !*** calculate dislocation multiplication
@ -3920,10 +3914,13 @@ dUpper(1:ns,1) = lattice_mu(phase) * burgers(1:ns,instance) &
/ (8.0_pReal * pi * (1.0_pReal - lattice_nu(phase)) * abs(tau)) / (8.0_pReal * pi * (1.0_pReal - lattice_nu(phase)) * abs(tau))
dUpper(1:ns,2) = lattice_mu(phase) * burgers(1:ns,instance) & dUpper(1:ns,2) = lattice_mu(phase) * burgers(1:ns,instance) &
/ (4.0_pReal * pi * abs(tau)) / (4.0_pReal * pi * abs(tau))
forall (c = 1_pInt:2_pInt) & forall (c = 1_pInt:2_pInt)
where(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)) >= tiny(0.0_pReal)) &
dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) &
+ abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), & + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
dUpper(1:ns,c)) dUpper(1:ns,c))
end forall
dUpper = max(dUpper,dLower) dUpper = max(dUpper,dLower)
!*** dislocation motion !*** dislocation motion
@ -4433,10 +4430,13 @@ dUpper(1:ns,1) = lattice_mu(phase) * burgers(1:ns,instance) &
/ (8.0_pReal * pi * (1.0_pReal - lattice_nu(phase)) * abs(tau)) / (8.0_pReal * pi * (1.0_pReal - lattice_nu(phase)) * abs(tau))
dUpper(1:ns,2) = lattice_mu(phase) * burgers(1:ns,instance) & dUpper(1:ns,2) = lattice_mu(phase) * burgers(1:ns,instance) &
/ (4.0_pReal * pi * abs(tau)) / (4.0_pReal * pi * abs(tau))
forall (c = 1_pInt:2_pInt) & forall (c = 1_pInt:2_pInt)
where(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)) >= tiny(0.0_pReal)) &
dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) &
+ abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), & + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
dUpper(1:ns,c)) dUpper(1:ns,c))
end forall
dUpper = max(dUpper,dLower) dUpper = max(dUpper,dLower)

8
code/prec.f90
View File

@ -16,11 +16,8 @@ module prec
implicit none implicit none
private private
#if (FLOAT==4) #if (FLOAT==4)
#ifdef Spectral #if defined(Spectral) || defined(FEM)
SPECTRAL SOLVER DOES NOT SUPPORT SINGLE PRECISION, STOPPING COMPILATION SPECTRAL SOLVER AND OWN FEM DO NOT SUPPORT SINGLE PRECISION, STOPPING COMPILATION
#endif
#ifdef FEM
SPECTRAL SOLVER DOES NOT SUPPORT SINGLE PRECISION, STOPPING COMPILATION
#endif #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 __INTEL_COMPILER #ifdef __INTEL_COMPILER
@ -81,7 +78,6 @@ module prec
end type end type
#endif #endif
public :: & public :: &
prec_init prec_init