diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 6688e89ad..660599efa 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -3705,17 +3705,17 @@ logical function crystallite_integrateStress(& !* feed local variables - Fp_current = crystallite_subFp0(1:3,1:3,g,i,e) ! "Fp_current" is only used as temp var here... - Lpguess = crystallite_Lp (1:3,1:3,g,i,e) ! ... and take it as first guess - Fi_current = crystallite_subFi0(1:3,1:3,g,i,e) ! intermediate configuration, assume decomposition as F = Fe Fi Fp - Liguess = crystallite_Li (1:3,1:3,g,i,e) ! ... and take it as first guess + Fp_current = crystallite_subFp0(1:3,1:3,g,i,e) ! "Fp_current" is only used as temp var here... + Lpguess = crystallite_Lp (1:3,1:3,g,i,e) ! ... and take it as first guess + Fi_current = crystallite_subFi0(1:3,1:3,g,i,e) ! intermediate configuration, assume decomposition as F = Fe Fi Fp + Liguess = crystallite_Li (1:3,1:3,g,i,e) ! ... and take it as first guess Liguess_old = Liguess !* 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 ',& @@ -3726,12 +3726,12 @@ logical function crystallite_integrateStress(& #endif return endif - A = math_mul33x33(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp + A = math_mul33x33(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp !* 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 ',& @@ -3890,7 +3890,7 @@ logical function crystallite_integrateStress(& endif deltaLp = - math_plain9to33(work) endif - jacoCounterLp = jacoCounterLp + 1_pInt ! increase counter for jaco update + jacoCounterLp = jacoCounterLp + 1_pInt ! increase counter for jaco update Lpguess = Lpguess + steplengthLp * deltaLp @@ -3910,20 +3910,20 @@ logical function crystallite_integrateStress(& !* update current residuum and check for convergence of loop - aTolLi = max(rTol_crystalliteStress * max(math_norm33(Liguess),math_norm33(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error - aTol_crystalliteStress) ! minimum lower cutoff + aTolLi = max(rTol_crystalliteStress * max(math_norm33(Liguess),math_norm33(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error + aTol_crystalliteStress) ! minimum lower cutoff residuumLi = Liguess - Li_constitutive - if (any(residuumLi /= residuumLi)) then ! NaN in residuum... - return ! ...me = .false. to inform integrator about problem - elseif (math_norm33(residuumLi) < aTolLi) then ! converged if below absolute tolerance - exit LiLoop ! ...leave iteration loop + if (any(residuumLi /= residuumLi)) then ! NaN in residuum... + return ! ...me = .false. to inform integrator about problem + elseif (math_norm33(residuumLi) < aTolLi) then ! converged if below absolute tolerance + exit LiLoop ! ...leave iteration loop elseif ( NiterationStressLi == 1_pInt & - .or. math_norm33(residuumLi) < math_norm33(residuumLi_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)... - residuumLi_old = residuumLi ! ...remember old values and... + .or. math_norm33(residuumLi) < math_norm33(residuumLi_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)... + residuumLi_old = residuumLi ! ...remember old values and... Liguess_old = Liguess - steplengthLi = steplengthLi0 ! ...proceed with normal step length (calculate new search direction) - else ! not converged and residuum not improved... - steplengthLi = 0.5_pReal * steplengthLi ! ...try with smaller step length in same direction + steplengthLi = steplengthLi0 ! ...proceed with normal step length (calculate new search direction) + else ! not converged and residuum not improved... + steplengthLi = 0.5_pReal * steplengthLi ! ...try with smaller step length in same direction Liguess = Liguess_old + steplengthLi * deltaLi cycle LiLoop endif @@ -3935,7 +3935,7 @@ logical function crystallite_integrateStress(& dFe_dLi3333 = 0.0_pReal dFi_dLi3333 = 0.0_pReal forall(o=1_pInt:3_pInt,p=1_pInt:3_pInt) - dFe_dLi3333(1:3,o,1:3,p) = -dt*math_I3(o,p)*temp_33 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) + dFe_dLi3333(1:3,o,1:3,p) = -dt*math_I3(o,p)*temp_33 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) dFi_dLi3333(1:3,o,1:3,p) = -dt*math_I3(o,p)*invFi_current end forall forall(o=1_pInt:3_pInt,p=1_pInt:3_pInt) & @@ -3947,16 +3947,16 @@ logical function crystallite_integrateStress(& - math_Plain3333to99(math_mul3333xx3333(dLi_dFi3333, dFi_dLi3333)) 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) - call sgesv(9,1,dRLi_dLi,9,ipiv,work,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li + 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 return endif deltaLi = - math_plain9to33(work) endif - jacoCounterLi = jacoCounterLi + 1_pInt ! increase counter for jaco update + jacoCounterLi = jacoCounterLi + 1_pInt ! increase counter for jaco update Liguess = Liguess + steplengthLi * deltaLi enddo LiLoop @@ -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 ',& diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index 9d5439c1c..7fe024039 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -464,8 +464,8 @@ 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 - homogenization_RGC_updateState = .true. ! pretend everything is fine and return + if(abs(dt) < tiny(0.0_pReal)) then ! zero time step + homogenization_RGC_updateState = .true. ! pretend everything is fine and return return endif diff --git a/code/lattice.f90 b/code/lattice.f90 index acf9df0f5..77e15c6de 100644 --- a/code/lattice.f90 +++ b/code/lattice.f90 @@ -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)) diff --git a/code/math.f90 b/code/math.f90 index 8fd49d516..dff1e103e 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -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 diff --git a/code/mesh.f90 b/code/mesh.f90 index 8b7ebce8d..f887d5ac1 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -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 diff --git a/code/plastic_dislotwin.f90 b/code/plastic_dislotwin.f90 index 5189f4f5a..5d00f560f 100644 --- a/code/plastic_dislotwin.f90 +++ b/code/plastic_dislotwin.f90 @@ -682,9 +682,9 @@ 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. & - lattice_structure(phase) == LATTICE_fcc_ID) & + 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) & call IO_error(211_pInt,el=instance,ext_msg='aTolRho ('//PLASTICITY_DISLOTWIN_label//')') @@ -708,8 +708,8 @@ 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. & - plastic_dislotwin_dipoleFormationFactor(instance) /= 1.0_pReal) & + 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. & plastic_dislotwin_qShearBand(instance) <= 0.0_pReal) & @@ -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 diff --git a/code/plastic_nonlocal.f90 b/code/plastic_nonlocal.f90 index 74c7a0d58..c1b94b70b 100644 --- a/code/plastic_nonlocal.f90 +++ b/code/plastic_nonlocal.f90 @@ -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 & diff --git a/code/plastic_titanmod.f90 b/code/plastic_titanmod.f90 index 80a9adaf5..b590c4fd4 100644 --- a/code/plastic_titanmod.f90 +++ b/code/plastic_titanmod.f90 @@ -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, &