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

comparison for zero should be abs: 

a == 0.0_pReal   becomes abs(a) <= tiny(a)
a /= 0.0_pReal   becomes abs(a) > tiny(a)
remove unused variables
This commit is contained in:
Martin Diehl 2015-04-11 08:25:23 +00:00
parent 6170209198
commit 470fc2dce3
8 changed files with 56 additions and 67 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
code/crystallite.f90
View File

@ -3715,7 +3715,7 @@ logical function crystallite_integrateStress(&
!* inversion of Fp_current...
invFp_current = math_inv33(Fp_current)
if (all(invFp_current <= tiny(0.0_pReal))) then ! math_inv33 returns zero when failed, avoid floating point comparison
if (all(abs(invFp_current) <= tiny(0.0_pReal))) then ! math_inv33 returns zero when failed, avoid floating point comparison
#ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fp_current at el (elFE) ip g ',&
@ -3731,7 +3731,7 @@ logical function crystallite_integrateStress(&
!* inversion of Fi_current...
invFi_current = math_inv33(Fi_current)
if (all(invFi_current <= tiny(0.0_pReal))) then ! math_inv33 returns zero when failed, avoid floating point comparison
if (all(abs(invFi_current) <= tiny(0.0_pReal))) then ! math_inv33 returns zero when failed, avoid floating point comparison
#ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fi_current at el (elFE) ip g ',&
@ -3973,7 +3973,7 @@ logical function crystallite_integrateStress(&
invFp_new = math_mul33x33(invFp_current,B)
invFp_new = invFp_new / math_det33(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize by det
Fp_new = math_inv33(invFp_new)
if (all(Fp_new <= tiny(0.0_pReal))) then ! math_inv33 returns zero when failed, avoid floating point comparison
if (all(abs(Fp_new)<= tiny(0.0_pReal))) then ! math_inv33 returns zero when failed, avoid floating point comparison
#ifndef _OPENMP
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 invFp_new inversion at el ip g ',&

2
code/homogenization_RGC.f90
View File

@ -464,7 +464,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix
real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax
if(dt < tiny(1.0_pReal)) then ! zero time step
if(abs(dt) < tiny(0.0_pReal)) then ! zero time step
homogenization_RGC_updateState = .true. ! pretend everything is fine and return
return
endif

12
code/lattice.f90
View File

@ -260,7 +260,7 @@ module lattice
],[ 4_pInt,LATTICE_fcc_Ntrans])
real(pReal), dimension(LATTICE_fcc_Ntrans,LATTICE_fcc_Ntrans), parameter, private :: & ! Matrix for projection of shear from slip system to fault-band (twin) systems
LATTICE_fcc_projectionTrans = reshape([& ! For ns = nt = nr
LATTICE_fcc_projectionTrans = reshape(real([& ! For ns = nt = nr
0, 1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
-1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
@ -273,7 +273,7 @@ module lattice
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,-1, &
0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 1, &
0, 0, 0, 0, 0, 0, 0, 0, 0, 1,-1, 0 &
],[LATTICE_fcc_Ntrans,LATTICE_fcc_Ntrans],order=[2,1])
],pReal),[LATTICE_fcc_Ntrans,LATTICE_fcc_Ntrans],order=[2,1])
real(pReal), parameter, private :: &
LATTICE_fcc_projectionTransFactor = sqrt(3.0_pReal/4.0_pReal)
@ -1374,12 +1374,12 @@ subroutine lattice_initializeStructure(myPhase,CoverA,aA,aM,cM)
rb(1:3,i) = lattice_fcc_bainRot(1:3,i)
ab(i) = lattice_fcc_bainRot(4,i)
xb(1:3,i) = LATTICE_fcc_bainVariant(1:3,i)
yb(1:3,i) = LATTICE_fcc_bainVariant(4:6,i)
zb(1:3,i) = LATTICE_fcc_bainVariant(7:9,i)
xb(1:3,i) = real(LATTICE_fcc_bainVariant(1:3,i),pReal)
yb(1:3,i) = real(LATTICE_fcc_bainVariant(4:6,i),pReal)
zb(1:3,i) = real(LATTICE_fcc_bainVariant(7:9,i),pReal)
ub(1:3,1:3,i) = 0.0_pReal
if ((aA > 0.0_pReal) .and. (aM > 0.0_pReal) .and. (cM == 0.0_pReal)) then
if ((aA > 0.0_pReal) .and. (aM > 0.0_pReal) .and. (abs(cM) <= tiny(0.0_pReal))) then
ub(1:3,1:3,i) = (aM/aA)*math_tensorproduct(xb(1:3,i), xb(1:3,i)) + &
sqrt(2.0_pReal)*(aM/aA)*math_tensorproduct(yb(1:3,i), yb(1:3,i)) + &
sqrt(2.0_pReal)*(aM/aA)*math_tensorproduct(zb(1:3,i), zb(1:3,i))

6
code/math.f90
View File

@ -745,7 +745,7 @@ pure function math_inv33(A)
- A(1,2) * (A(2,1) * A(3,3) - A(2,3) * A(3,1))&
+ A(1,3) * (A(2,1) * A(3,2) - A(2,2) * A(3,1))
if (abs(DetA) > tiny(abs(DetA))) then
if (abs(DetA) > tiny(DetA)) then
math_inv33(1,1) = ( A(2,2) * A(3,3) - A(2,3) * A(3,2)) / DetA
math_inv33(2,1) = (-A(2,1) * A(3,3) + A(2,3) * A(3,1)) / DetA
math_inv33(3,1) = ( A(2,1) * A(3,2) - A(2,2) * A(3,1)) / DetA
@ -783,7 +783,7 @@ pure subroutine math_invert33(A, InvA, DetA, error)
- A(1,2) * (A(2,1) * A(3,3) - A(2,3) * A(3,1))&
+ A(1,3) * (A(2,1) * A(3,2) - A(2,2) * A(3,1))
if (abs(DetA) <= tiny(abs(DetA))) then
if (abs(DetA) <= tiny(DetA)) then
error = .true.
else
InvA(1,1) = ( A(2,2) * A(3,3) - A(2,3) * A(3,2)) / DetA
@ -1318,7 +1318,7 @@ pure function math_qInv(Q)
math_qInv = 0.0_pReal
squareNorm = math_qDot(Q,Q)
if (squareNorm > tiny(squareNorm)) &
if (abs(squareNorm) > tiny(squareNorm)) &
math_qInv = math_qConj(Q) / squareNorm
end function math_qInv

14
code/mesh.f90
View File

@ -567,16 +567,13 @@ subroutine mesh_init(ip,el)
gridMPI = gridGlobal
alloc_local = fftw_mpi_local_size_3d(gridMPI(3), gridMPI(2), gridMPI(1)/2 +1, &
MPI_COMM_WORLD, local_K, local_K_offset)
gridLocal(1) = gridGlobal(1)
gridLocal(2) = gridGlobal(2)
gridLocal(3) = local_K
gridOffset = local_K_offset
gridLocal = [gridGlobal(1:2) int(local_K,pInt)]
gridOffset = int(local_K_offset,pInt)
geomSizeGlobal = mesh_spectral_getSize(fileUnit)
geomSizeLocal(1) = geomSizeGlobal(1)
geomSizeLocal(2) = geomSizeGlobal(2)
geomSizeLocal(3) = geomSizeGlobal(3)*real(gridLocal(3))/real(gridGlobal(3))
geomSizeOffset = geomSizeGlobal(3)*real(gridOffset) /real(gridGlobal(3))
geomSizeLocal = [geomSizeGlobal(1:2), &
geomSizeGlobal(3)*real(gridLocal(3),pReal)/real(gridGlobal(3),pReal)]
geomSizeOffset = geomSizeGlobal(3)*real(gridOffset,pReal) /real(gridGlobal(3),pReal)
#else
call IO_open_file(FILEUNIT,geometryFile) ! parse info from geometry file...
if (myDebug) write(6,'(a)') ' Opened geometry file'; flush(6)
@ -2206,7 +2203,6 @@ function mesh_volumeMismatch(gDim,F,nodes) result(vMismatch)
debug_level, &
debug_levelBasic
use math, only: &
PI, &
math_det33, &
math_volTetrahedron

14
code/plastic_dislotwin.f90
View File

@ -682,8 +682,8 @@ subroutine plastic_dislotwin_init(fileUnit)
if (plastic_dislotwin_Qsd(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='Qsd ('//PLASTICITY_DISLOTWIN_label//')')
if (sum(plastic_dislotwin_Ntwin(:,instance)) > 0_pInt) then
if (plastic_dislotwin_SFE_0K(instance) <= tiny(0.0_pReal) .and. &
plastic_dislotwin_dSFE_dT(instance) <= tiny(0.0_pReal) .and. &
if (abs(plastic_dislotwin_SFE_0K(instance)) <= tiny(0.0_pReal) .and. &
abs(plastic_dislotwin_dSFE_dT(instance)) <= tiny(0.0_pReal) .and. &
lattice_structure(phase) == LATTICE_fcc_ID) &
call IO_error(211_pInt,el=instance,ext_msg='SFE0K ('//PLASTICITY_DISLOTWIN_label//')')
if (plastic_dislotwin_aTolRho(instance) <= 0.0_pReal) &
@ -708,7 +708,7 @@ subroutine plastic_dislotwin_init(fileUnit)
if (plastic_dislotwin_sbVelocity(instance) > 0.0_pReal .and. &
plastic_dislotwin_pShearBand(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='pShearBand ('//PLASTICITY_DISLOTWIN_label//')')
if (plastic_dislotwin_dipoleFormationFactor(instance) > tiny(0.0_pReal) .and. &
if (abs(plastic_dislotwin_dipoleFormationFactor(instance)) > tiny(0.0_pReal) .and. &
plastic_dislotwin_dipoleFormationFactor(instance) /= 1.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='dipoleFormationFactor ('//PLASTICITY_DISLOTWIN_label//')')
if (plastic_dislotwin_sbVelocity(instance) > 0.0_pReal .and. &
@ -1507,8 +1507,8 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature
!--------------------------------------------------------------------------------------------------
! Shear banding (shearband) part
if(plastic_dislotwin_sbVelocity(instance) > tiny(0.0_pReal) .and. &
plastic_dislotwin_sbResistance(instance) > tiny(0.0_pReal)) then
if(abs(plastic_dislotwin_sbVelocity(instance)) > tiny(0.0_pReal) .and. &
abs(plastic_dislotwin_sbResistance(instance)) > tiny(0.0_pReal)) then
gdot_sb = 0.0_pReal
dgdot_dtausb = 0.0_pReal
call math_spectralDecompositionSym33(math_Mandel6to33(Tstar_v),eigValues,eigVectors, error)
@ -1792,7 +1792,7 @@ subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el)
plastic_dislotwin_CAtomicVolume(instance)*plastic_dislotwin_burgersPerSlipSystem(j,instance)**(3.0_pReal)
VacancyDiffusion = &
plastic_dislotwin_D0(instance)*exp(-plastic_dislotwin_Qsd(instance)/(kB*Temperature))
if (tau_slip(j) <= tiny(0.0_pReal)) then
if (abs(tau_slip(j)) <= tiny(0.0_pReal)) then
DotRhoEdgeDipClimb(j) = 0.0_pReal
else
ClimbVelocity(j) = &
@ -2255,7 +2255,7 @@ function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
endif
!* Stress exponent
if (gdot_slip(j)<=tiny(0.0_pReal)) then
if (abs(gdot_slip(j))<=tiny(0.0_pReal)) then
plastic_dislotwin_postResults(c+j) = 0.0_pReal
else
plastic_dislotwin_postResults(c+j) = (tau/gdot_slip(j))*dgdot_dtauslip

17
code/plastic_nonlocal.f90
View File

@ -286,14 +286,12 @@ use debug, only: debug_level, &
use mesh, only: mesh_NcpElems, &
mesh_maxNips, &
mesh_maxNipNeighbors
use material, only: homogenization_maxNgrains, &
phase_plasticity, &
use material, only: phase_plasticity, &
phase_plasticityInstance, &
phase_Noutput, &
PLASTICITY_NONLOCAL_label, &
PLASTICITY_NONLOCAL_ID, &
plasticState, &
! material_phase, &
material_Nphase, &
MATERIAL_partPhase ,&
material_phase
@ -1420,7 +1418,6 @@ use material, only: material_phase, &
phase_plasticityInstance, &
plasticState, &
mappingConstitutive, &
material_Nphase, &
phase_plasticity ,&
PLASTICITY_NONLOCAL_ID
implicit none
@ -1794,7 +1791,7 @@ if (.not. phase_localPlasticity(ph) .and. shortRangeStressCorrection(instance))
- neighbor_rhoExcess(c,s,neighbors(2))
enddo
invConnections = math_inv33(connections)
if (all(invConnections <= tiny(0.0_pReal))) & ! check for failed in version (math_inv33 returns 0) and avoid floating point equality comparison
if (all(abs(invConnections) <= tiny(0.0_pReal))) & ! check for failed in version (math_inv33 returns 0) and avoid floating point equality comparison
call IO_error(-1_pInt,ext_msg='back stress calculation: inversion error')
rhoExcessGradient(c) = math_mul3x3(m(1:3,s,c), &
math_mul33x3(invConnections,rhoExcessDifferences))
@ -2338,7 +2335,7 @@ deltaDUpper = dUpper - dUpperOld
!*** dissociation by stress increase
deltaRhoDipole2SingleStress = 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)) > tiny(0.0_pReal)) &
abs(dUpperOld(s,c) - dLower(s,c)) > tiny(0.0_pReal)) &
deltaRhoDipole2SingleStress(s,8_pInt+c) = rhoDip(s,c) * deltaDUpper(s,c) &
/ (dUpperOld(s,c) - dLower(s,c))
@ -2834,11 +2831,11 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then
my_rhoSgl = rhoSgl
my_v = v
if(numerics_timeSyncing) then
if (subfrac(1_pInt,ip,el) <= tiny(0.0_pReal)) then
if (abs(subfrac(1_pInt,ip,el))<= tiny(0.0_pReal)) then
my_rhoSgl = rhoSgl0
my_v = v0
elseif (neighbor_n > 0_pInt) then
if (subfrac(1_pInt,neighbor_ip,neighbor_el) <= tiny(0.0_pReal)) then
if (abs(subfrac(1_pInt,neighbor_ip,neighbor_el))<= tiny(0.0_pReal)) then
my_rhoSgl = rhoSgl0
my_v = v0
endif
@ -3394,7 +3391,7 @@ if (.not. phase_localPlasticity(ph)) then
Rsquare = R * R
Rcube = Rsquare * R
denominator = R * (R + flipSign * lambda)
if (denominator <= tiny(0.0_pReal)) exit ipLoop
if (abs(denominator)<= tiny(0.0_pReal)) exit ipLoop
sigma(1,1) = sigma(1,1) - real(side,pReal) &
* flipSign * z / denominator &
@ -3439,7 +3436,7 @@ if (.not. phase_localPlasticity(ph)) then
Rsquare = R * R
Rcube = Rsquare * R
denominator = R * (R + flipSign * lambda)
if (denominator <= tiny(0.0_pReal)) exit ipLoop
if (abs(denominator)<= tiny(0.0_pReal)) exit ipLoop
sigma(1,2) = sigma(1,2) - real(side,pReal) * flipSign * z &
* (1.0_pReal - lattice_nu(ph)) / denominator &

4
code/plastic_titanmod.f90
View File

@ -1225,11 +1225,7 @@ end function plastic_titanmod_homogenizedC
!--------------------------------------------------------------------------------------------------
subroutine plastic_titanmod_microstructure(temperature,ipc,ip,el)
use mesh, only: &
mesh_NcpElems, &
mesh_maxNips
use material, only: &
homogenization_maxNgrains, &
material_phase,&
phase_plasticityInstance, &
plasticState, &