diff --git a/code/crystallite.f90 b/code/crystallite.f90 index b50c71272..450da5685 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -843,7 +843,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !$OMP DO PRIVATE(neighboring_e,neighboring_i) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - if (.not. crystallite_localPlasticity(1,i,e) .and. abs(crystallite_subFrac(1,i,e)) > tiny(0.0_pReal)) then + if (.not. crystallite_localPlasticity(1,i,e) .and. dNeq(crystallite_subFrac(1,i,e),0.0_pReal)) then do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e)))) neighboring_e = mesh_ipNeighborhood(1,n,i,e) neighboring_i = mesh_ipNeighborhood(2,n,i,e) @@ -3623,6 +3623,7 @@ logical function crystallite_integrateStress(& !* inversion of Fi_current... + invFi_current = math_inv33(Fi_current) failedInversionFi: if (all(dEq(invFi_current,0.0_pReal))) then #ifndef _OPENMP if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then @@ -3883,7 +3884,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) - failedInversionFp2: if (all(dEq(invFp_new,0.0_pReal))) then + failedInversionInvFp: if (all(dEq(Fp_new,0.0_pReal))) then #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 ipc ',& @@ -3895,7 +3896,7 @@ logical function crystallite_integrateStress(& endif #endif return - endif failedInversionFp2 + endif failedInversionInvFp Fe_new = math_mul33x33(math_mul33x33(Fg_new,invFp_new),invFi_new) ! calc resulting Fe !* calculate 1st Piola-Kirchhoff stress diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index 0919b1e5e..3b75ec016 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -386,6 +386,8 @@ end subroutine homogenization_RGC_partitionDeformation ! "happy" with result !-------------------------------------------------------------------------------------------------- function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) + use prec, only: & + dEq use debug, only: & debug_level, & debug_homogenization,& @@ -441,10 +443,10 @@ 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(abs(dt) < tiny(0.0_pReal)) then ! zero time step - homogenization_RGC_updateState = .true. ! pretend everything is fine and return + zeroTimeStep: if(dEq(dt,0.0_pReal)) then + homogenization_RGC_updateState = .true. ! pretend everything is fine and return return - endif + endif zeroTimeStep !-------------------------------------------------------------------------------------------------- ! get the dimension of the cluster (grains and interfaces) diff --git a/code/math.f90 b/code/math.f90 index a5a525461..12803ba3d 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -761,6 +761,8 @@ end function math_inv33 ! returns error if not possible, i.e. if det close to zero !-------------------------------------------------------------------------------------------------- pure subroutine math_invert33(A, InvA, DetA, error) + use prec, only: & + dEq implicit none logical, intent(out) :: error @@ -774,7 +776,7 @@ pure subroutine math_invert33(A, InvA, DetA, error) DetA = A(1,1) * InvA(1,1) + A(1,2) * InvA(2,1) + A(1,3) * InvA(3,1) - if (abs(DetA) <= tiny(DetA)) then + if (dEq(DetA,0.0_pReal)) then InvA = 0.0_pReal error = .true. else @@ -1279,6 +1281,8 @@ end function math_qNorm !> @brief quaternion inversion !-------------------------------------------------------------------------------------------------- pure function math_qInv(Q) + use prec, only: & + dNeq implicit none real(pReal), dimension(4), intent(in) :: Q @@ -1288,8 +1292,7 @@ pure function math_qInv(Q) math_qInv = 0.0_pReal squareNorm = math_qDot(Q,Q) - if (abs(squareNorm) > tiny(squareNorm)) & - math_qInv = math_qConj(Q) / squareNorm + if (dNeq(squareNorm,0.0_pReal)) math_qInv = math_qConj(Q) / squareNorm end function math_qInv @@ -2093,6 +2096,8 @@ end function math_eigenvectorBasisSym33 !> @brief rotational part from polar decomposition of 33 tensor m !-------------------------------------------------------------------------------------------------- function math_rotationalPart33(m) + use prec, only: & + dEq use IO, only: & IO_warning @@ -2104,12 +2109,12 @@ function math_rotationalPart33(m) U = math_eigenvectorBasisSym33(math_mul33x33(transpose(m),m)) Uinv = math_inv33(U) - if (all(abs(Uinv) <= tiny(Uinv))) then ! math_inv33 returns zero when failed, avoid floating point equality comparison + inversionFailed: if (all(dEq(Uinv,0.0_pReal))) then math_rotationalPart33 = math_I3 call IO_warning(650_pInt) - else + else inversionFailed math_rotationalPart33 = math_mul33x33(m,Uinv) - endif + endif inversionFailed end function math_rotationalPart33 diff --git a/code/plastic_disloUCLA.f90 b/code/plastic_disloUCLA.f90 index 11863c5c4..0e7105c70 100644 --- a/code/plastic_disloUCLA.f90 +++ b/code/plastic_disloUCLA.f90 @@ -981,7 +981,8 @@ end subroutine plastic_disloUCLA_LpAndItsTangent !-------------------------------------------------------------------------------------------------- subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,ipc,ip,el) use prec, only: & - tol_math_check + tol_math_check, & + dEq use math, only: & pi use material, only: & @@ -1119,7 +1120,7 @@ subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,ipc,ip,el) !* Dipole formation EdgeDipMinDistance = & plastic_disloUCLA_CEdgeDipMinDistance(instance)*plastic_disloUCLA_burgersPerSlipSystem(j,instance) - if (abs(tau_slip_pos) <= tiny(0.0_pReal)) then + if (dEq(tau_slip_pos,0.0_pReal)) then DotRhoDipFormation = 0.0_pReal else EdgeDipDistance = & @@ -1147,7 +1148,7 @@ subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,ipc,ip,el) plastic_disloUCLA_CAtomicVolume(instance)*plastic_disloUCLA_burgersPerSlipSystem(j,instance)**(3.0_pReal) VacancyDiffusion = & plastic_disloUCLA_D0(instance)*exp(-plastic_disloUCLA_Qsd(instance)/(kB*Temperature)) - if (abs(tau_slip_pos) <= tiny(0.0_pReal)) then + if (dEq(tau_slip_pos,0.0_pReal)) then DotRhoEdgeDipClimb = 0.0_pReal else ClimbVelocity = & @@ -1180,7 +1181,8 @@ end subroutine plastic_disloUCLA_dotState !-------------------------------------------------------------------------------------------------- function plastic_disloUCLA_postResults(Tstar_v,Temperature,ipc,ip,el) use prec, only: & - tol_math_check + tol_math_check, & + dEq use math, only: & pi use material, only: & @@ -1408,7 +1410,7 @@ function plastic_disloUCLA_postResults(Tstar_v,Temperature,ipc,ip,el) c = c + ns elseif(plastic_disloUCLA_outputID(o,instance) == stress_exponent_ID) then do j = 1_pInt, ns - if (abs(gdot_slip_pos(j)+gdot_slip_neg(j))<=tiny(0.0_pReal)) then + if (dEq(gdot_slip_pos(j)+gdot_slip_neg(j),0.0_pReal)) then plastic_disloUCLA_postResults(c+j) = 0.0_pReal else plastic_disloUCLA_postResults(c+j) = (tau_slip_pos(j)+tau_slip_neg(j))/& diff --git a/code/plastic_dislotwin.f90 b/code/plastic_dislotwin.f90 index ebee2738a..0902ecad6 100644 --- a/code/plastic_dislotwin.f90 +++ b/code/plastic_dislotwin.f90 @@ -200,6 +200,7 @@ contains subroutine plastic_dislotwin_init(fileUnit) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use prec, only: & + dEq, & dNeq use debug, only: & debug_level,& @@ -751,8 +752,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 (abs(plastic_dislotwin_SFE_0K(instance)) <= tiny(0.0_pReal) .and. & - abs(plastic_dislotwin_dSFE_dT(instance)) <= tiny(0.0_pReal) .and. & + if (dEq(plastic_dislotwin_SFE_0K(instance), 0.0_pReal) .and. & + dEq(plastic_dislotwin_dSFE_dT(instance),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) & @@ -761,8 +762,8 @@ subroutine plastic_dislotwin_init(fileUnit) call IO_error(211_pInt,el=instance,ext_msg='aTolTwinFrac ('//PLASTICITY_DISLOTWIN_label//')') endif if (sum(plastic_dislotwin_Ntrans(:,instance)) > 0_pInt) then - if (abs(plastic_dislotwin_SFE_0K(instance)) <= tiny(0.0_pReal) .and. & - abs(plastic_dislotwin_dSFE_dT(instance)) <= tiny(0.0_pReal) .and. & + if (dEq(plastic_dislotwin_SFE_0K(instance), 0.0_pReal) .and. & + dEq(plastic_dislotwin_dSFE_dT(instance),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_aTolTransFrac(instance) <= 0.0_pReal) & @@ -1630,7 +1631,8 @@ end subroutine plastic_dislotwin_microstructure !-------------------------------------------------------------------------------------------------- subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,ipc,ip,el) use prec, only: & - tol_math_check + tol_math_check, & + dNeq use math, only: & math_Plain3333to99, & math_Mandel6to33, & @@ -1777,8 +1779,8 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature !-------------------------------------------------------------------------------------------------- ! Shear banding (shearband) part - if(abs(plastic_dislotwin_sbVelocity(instance)) > tiny(0.0_pReal) .and. & - abs(plastic_dislotwin_sbResistance(instance)) > tiny(0.0_pReal)) then + if(dNeq(plastic_dislotwin_sbVelocity(instance), 0.0_pReal) .and. & + dNeq(plastic_dislotwin_sbResistance(instance),0.0_pReal)) then gdot_sb = 0.0_pReal dgdot_dtausb = 0.0_pReal call math_eigenValuesVectorsSym(math_Mandel6to33(Tstar_v),eigValues,eigVectors,error) @@ -1944,7 +1946,8 @@ end subroutine plastic_dislotwin_LpAndItsTangent !-------------------------------------------------------------------------------------------------- subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el) use prec, only: & - tol_math_check + tol_math_check, & + dEq use math, only: & pi use material, only: & @@ -2045,7 +2048,7 @@ subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el) !* Dipole formation EdgeDipMinDistance = & plastic_dislotwin_CEdgeDipMinDistance(instance)*plastic_dislotwin_burgersPerSlipSystem(j,instance) - if (abs(tau_slip(j)) <= tiny(0.0_pReal)) then + if (dEq(tau_slip(j),0.0_pReal)) then DotRhoDipFormation = 0.0_pReal else EdgeDipDistance = & @@ -2073,10 +2076,10 @@ 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 (abs(tau_slip(j)) <= tiny(0.0_pReal)) then + if (dEq(tau_slip(j),0.0_pReal)) then DotRhoEdgeDipClimb = 0.0_pReal else - if (EdgeDipDistance-EdgeDipMinDistance <= tiny(0.0_pReal)) then + if (dEq(EdgeDipDistance-EdgeDipMinDistance,0.0_pReal)) then DotRhoEdgeDipClimb = 0.0_pReal else ClimbVelocity = 3.0_pReal*lattice_mu(ph)*VacancyDiffusion*AtomicVolume/ & @@ -2191,7 +2194,8 @@ end subroutine plastic_dislotwin_dotState !-------------------------------------------------------------------------------------------------- function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el) use prec, only: & - tol_math_check + tol_math_check, & + dEq use math, only: & pi, & math_Mandel6to33, & @@ -2506,11 +2510,8 @@ function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el) endif !* Stress exponent - 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 - endif + plastic_dislotwin_postResults(c+j) = & + merge(0.0_pReal,(tau/gdot_slip(j))*dgdot_dtauslip,dEq(gdot_slip(j),0.0_pReal)) enddo ; enddo c = c + ns case (sb_eigenvalues_ID) diff --git a/code/plastic_isotropic.f90 b/code/plastic_isotropic.f90 index a074c3862..7d3e1aa7c 100644 --- a/code/plastic_isotropic.f90 +++ b/code/plastic_isotropic.f90 @@ -524,6 +524,8 @@ end subroutine plastic_isotropic_LiAndItsTangent !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el) + use prec, only: & + dEq use math, only: & math_mul6x6 use material, only: & @@ -570,7 +572,7 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el) !-------------------------------------------------------------------------------------------------- ! hardening coefficient if (abs(gamma_dot) > 1e-12_pReal) then - if (abs(param(instance)%tausat_SinhFitA) <= tiny(0.0_pReal)) then + if (dEq(param(instance)%tausat_SinhFitA,0.0_pReal)) then saturation = param(instance)%tausat else saturation = ( param(instance)%tausat & diff --git a/code/plastic_nonlocal.f90 b/code/plastic_nonlocal.f90 index 4557ca237..a6c7acad1 100644 --- a/code/plastic_nonlocal.f90 +++ b/code/plastic_nonlocal.f90 @@ -1549,6 +1549,8 @@ end subroutine plastic_nonlocal_aTolState !> @brief calculates quantities characterizing the microstructure !-------------------------------------------------------------------------------------------------- subroutine plastic_nonlocal_microstructure(Fe, Fp, ip, el) +use prec, only: & + dEq use IO, only: & IO_error use math, only: & @@ -1792,7 +1794,7 @@ if (.not. phase_localPlasticity(ph) .and. shortRangeStressCorrection(instance)) - neighbor_rhoExcess(c,s,neighbors(2)) enddo invConnections = math_inv33(connections) - if (all(abs(invConnections) <= tiny(0.0_pReal))) & ! check for failed in version (math_inv33 returns 0) and avoid floating point equality comparison + if (all(dEq(invConnections,0.0_pReal))) & 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)) @@ -2200,6 +2202,8 @@ end subroutine plastic_nonlocal_LpAndItsTangent !> @brief (instantaneous) incremental change of microstructure !-------------------------------------------------------------------------------------------------- subroutine plastic_nonlocal_deltaState(Tstar_v,ip,el) +use prec, only: & + dNeq use debug, only: debug_level, & debug_constitutive, & debug_levelBasic, & @@ -2322,8 +2326,8 @@ dUpper(1:ns,2) = lattice_mu(ph) * burgers(1:ns,instance) / (4.0_pReal * pi * abs 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)) & + where(dNeq(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)),0.0_pReal)) & 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)), & dUpper(1:ns,c)) @@ -2335,7 +2339,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. & - abs(dUpperOld(s,c) - dLower(s,c)) > tiny(0.0_pReal)) & + dNeq(dUpperOld(s,c) - dLower(s,c),0.0_pReal)) & deltaRhoDipole2SingleStress(s,8_pInt+c) = rhoDip(s,c) * deltaDUpper(s,c) & / (dUpperOld(s,c) - dLower(s,c)) @@ -2383,7 +2387,8 @@ subroutine plastic_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, & timestep,subfrac, ip,el) use prec, only: DAMASK_NaN, & - dNeq + dNeq, & + dEq use numerics, only: numerics_integrationMode, & numerics_timeSyncing use IO, only: IO_error @@ -2617,8 +2622,8 @@ dUpper(1:ns,1) = lattice_mu(ph) * burgers(1:ns,instance) & dUpper(1:ns,2) = lattice_mu(ph) * burgers(1:ns,instance) & / (4.0_pReal * pi * abs(tau)) 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)) & + where(dNeq(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)),0.0_pReal)) & 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)), & dUpper(1:ns,c)) @@ -2830,11 +2835,11 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then my_rhoSgl = rhoSgl my_v = v if(numerics_timeSyncing) then - if (abs(subfrac(1_pInt,ip,el))<= tiny(0.0_pReal)) then + if (dEq(subfrac(1_pInt,ip,el),0.0_pReal)) then my_rhoSgl = rhoSgl0 my_v = v0 elseif (neighbor_n > 0_pInt) then - if (abs(subfrac(1_pInt,neighbor_ip,neighbor_el))<= tiny(0.0_pReal)) then + if (dEq(subfrac(1_pInt,neighbor_ip,neighbor_el),0.0_pReal)) then my_rhoSgl = rhoSgl0 my_v = v0 endif @@ -3172,6 +3177,8 @@ end subroutine plastic_nonlocal_updateCompatibility !* calculates quantities characterizing the microstructure * !********************************************************************* function plastic_nonlocal_dislocationstress(Fe, ip, el) +use prec, only: & + dEq use math, only: math_mul33x33, & math_mul33x3, & math_inv33, & @@ -3384,7 +3391,7 @@ if (.not. phase_localPlasticity(ph)) then Rsquare = R * R Rcube = Rsquare * R denominator = R * (R + flipSign * lambda) - if (abs(denominator)<= tiny(0.0_pReal)) exit ipLoop + if (dEq(denominator,0.0_pReal)) exit ipLoop sigma(1,1) = sigma(1,1) - real(side,pReal) & * flipSign * z / denominator & @@ -3429,7 +3436,7 @@ if (.not. phase_localPlasticity(ph)) then Rsquare = R * R Rcube = Rsquare * R denominator = R * (R + flipSign * lambda) - if (abs(denominator)<= tiny(0.0_pReal)) exit ipLoop + if (dEq(denominator,0.0_pReal)) exit ipLoop sigma(1,2) = sigma(1,2) - real(side,pReal) * flipSign * z & * (1.0_pReal - lattice_nu(ph)) / denominator & @@ -3518,6 +3525,8 @@ end function plastic_nonlocal_dislocationstress !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- function plastic_nonlocal_postResults(Tstar_v,Fe,ip,el) + use prec, only: & + dNeq use math, only: & math_mul6x6, & math_mul33x3, & @@ -3634,8 +3643,8 @@ dUpper(1:ns,1) = lattice_mu(ph) * burgers(1:ns,instance) & dUpper(1:ns,2) = lattice_mu(ph) * burgers(1:ns,instance) & / (4.0_pReal * pi * abs(tau)) 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)) & + where(dNeq(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)),0.0_pReal)) & 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)), & dUpper(1:ns,c)) diff --git a/code/plastic_phenoplus.f90 b/code/plastic_phenoplus.f90 index 96ad0d647..5df3d490b 100644 --- a/code/plastic_phenoplus.f90 +++ b/code/plastic_phenoplus.f90 @@ -112,6 +112,8 @@ contains !-------------------------------------------------------------------------------------------------- subroutine plastic_phenoplus_init(fileUnit) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use prec, only: & + dEq use debug, only: & debug_level, & debug_constitutive,& @@ -481,7 +483,7 @@ subroutine plastic_phenoplus_init(fileUnit) if (any(plastic_phenoplus_tausat_slip(:,instance) <= 0.0_pReal .and. & plastic_phenoplus_Nslip(:,instance) > 0)) & call IO_error(211_pInt,el=instance,ext_msg='tausat_slip ('//PLASTICITY_PHENOPLUS_label//')') - if (any(abs(plastic_phenoplus_a_slip(instance)) <= tiny(0.0_pReal) .and. & + if (any(dEq(plastic_phenoplus_a_slip(instance),0.0_pReal) .and. & plastic_phenoplus_Nslip(:,instance) > 0)) & call IO_error(211_pInt,el=instance,ext_msg='a_slip ('//PLASTICITY_PHENOPLUS_label//')') if (any(plastic_phenoplus_tau0_twin(:,instance) < 0.0_pReal .and. & @@ -923,6 +925,8 @@ end subroutine plastic_phenoplus_microstructure !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) + use prec, only: & + dNeq use math, only: & math_Plain3333to99, & math_Mandel6to33 @@ -1038,7 +1042,7 @@ subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) (gdot_slip_pos+gdot_slip_neg)*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph) ! Calculation of the tangent of Lp - if (abs(gdot_slip_pos) > tiny(0.0_pReal)) then + if (dNeq(gdot_slip_pos,0.0_pReal)) then dgdot_dtauslip_pos = gdot_slip_pos*plastic_phenoplus_n_slip(instance)/tau_slip_pos forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & @@ -1046,7 +1050,7 @@ subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) nonSchmid_tensor(m,n,1) endif - if (abs(gdot_slip_neg) > tiny(0.0_pReal)) then + if (dNeq(gdot_slip_neg,0.0_pReal)) then dgdot_dtauslip_neg = gdot_slip_neg*plastic_phenoplus_n_slip(instance)/tau_slip_neg forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & @@ -1073,7 +1077,7 @@ subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) Lp = Lp + gdot_twin*lattice_Stwin(1:3,1:3,index_myFamily+i,ph) ! Calculation of the tangent of Lp - if (abs(gdot_twin) > tiny(0.0_pReal)) then + if (dNeq(gdot_twin,0.0_pReal)) then dgdot_dtautwin = gdot_twin*plastic_phenoplus_n_twin(instance)/tau_twin forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & diff --git a/code/plastic_phenopowerlaw.f90 b/code/plastic_phenopowerlaw.f90 index 8bc86c97d..15b7aa78b 100644 --- a/code/plastic_phenopowerlaw.f90 +++ b/code/plastic_phenopowerlaw.f90 @@ -124,6 +124,8 @@ contains !-------------------------------------------------------------------------------------------------- subroutine plastic_phenopowerlaw_init(fileUnit) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use prec, only: & + dEq use debug, only: & debug_level, & debug_constitutive,& @@ -487,7 +489,7 @@ subroutine plastic_phenopowerlaw_init(fileUnit) if (any(plastic_phenopowerlaw_tausat_slip(:,instance) <= 0.0_pReal .and. & plastic_phenopowerlaw_Nslip(:,instance) > 0)) & call IO_error(211_pInt,el=instance,ext_msg='tausat_slip ('//PLASTICITY_PHENOPOWERLAW_label//')') - if (any(abs(plastic_phenopowerlaw_a_slip(instance)) <= tiny(0.0_pReal) .and. & + if (any(dEq(plastic_phenopowerlaw_a_slip(instance),0.0_pReal) .and. & plastic_phenopowerlaw_Nslip(:,instance) > 0)) & call IO_error(211_pInt,el=instance,ext_msg='a_slip ('//PLASTICITY_PHENOPOWERLAW_label//')') if (any(plastic_phenopowerlaw_tau0_twin(:,instance) < 0.0_pReal .and. & @@ -774,6 +776,8 @@ end subroutine plastic_phenopowerlaw_aTolState !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) + use prec, only: & + dNeq use math, only: & math_Plain3333to99, & math_Mandel6to33 @@ -863,7 +867,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip, (gdot_slip_pos+gdot_slip_neg)*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph) ! Calculation of the tangent of Lp - if (abs(gdot_slip_pos) > tiny(0.0_pReal)) then + if (dNeq(gdot_slip_pos,0.0_pReal)) then dgdot_dtauslip_pos = gdot_slip_pos*plastic_phenopowerlaw_n_slip(instance)/tau_slip_pos forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & @@ -871,7 +875,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip, nonSchmid_tensor(m,n,1) endif - if (abs(gdot_slip_neg) > tiny(0.0_pReal)) then + if (dNeq(gdot_slip_neg,0.0_pReal)) then dgdot_dtauslip_neg = gdot_slip_neg*plastic_phenopowerlaw_n_slip(instance)/tau_slip_neg forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & @@ -898,7 +902,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip, Lp = Lp + gdot_twin*lattice_Stwin(1:3,1:3,index_myFamily+i,ph) ! Calculation of the tangent of Lp - if (abs(gdot_twin) > tiny(0.0_pReal)) then + if (dNeq(gdot_twin,0.0_pReal)) then dgdot_dtautwin = gdot_twin*plastic_phenopowerlaw_n_twin(instance)/tau_twin forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & diff --git a/code/prec.f90 b/code/prec.f90 index d9375bafa..594201bee 100644 --- a/code/prec.f90 +++ b/code/prec.f90 @@ -115,7 +115,9 @@ module prec prec_init, & prec_isNaN, & dEq, & - dNeq + cEq, & + dNeq, & + cNeq contains @@ -180,23 +182,23 @@ end function prec_isNaN !-------------------------------------------------------------------------------------------------- -!> @brief equality comparison for double precision +!> @brief equality comparison for float with double precision ! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq ! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm !-------------------------------------------------------------------------------------------------- logical elemental pure function dEq(a,b,tol) implicit none - real(pReal), intent(in) :: a,b + real(pReal), intent(in) :: a,b real(pReal), intent(in), optional :: tol - real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C + real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C dEq = merge(.True., .False.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b]))) end function dEq !-------------------------------------------------------------------------------------------------- -!> @brief inequality comparison for double precision +!> @brief inequality comparison for float with double precision ! replaces "!=" but for certain (relative) tolerance. Counterpart to dEq ! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm !-------------------------------------------------------------------------------------------------- @@ -205,9 +207,44 @@ logical elemental pure function dNeq(a,b,tol) implicit none real(pReal), intent(in) :: a,b real(pReal), intent(in), optional :: tol - real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C + real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C dNeq = merge(.False., .True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b]))) end function dNeq +!-------------------------------------------------------------------------------------------------- +!> @brief equality comparison for complex with double precision +! replaces "==" but for certain (relative) tolerance. Counterpart to cNeq +! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm +! probably a component wise comparison would be more accurate than the comparsion of the absolute +! value +!-------------------------------------------------------------------------------------------------- +logical elemental pure function cEq(a,b,tol) + + implicit none + complex(pReal), intent(in) :: a,b + real(pReal), intent(in), optional :: tol + real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C + + cEq = merge(.True., .False.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b]))) +end function cEq + + +!-------------------------------------------------------------------------------------------------- +!> @brief inequality comparison for complex with double precision +! replaces "!=" but for certain (relative) tolerance. Counterpart to cEq +! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm +! probably a component wise comparison would be more accurate than the comparsion of the absolute +! value +!-------------------------------------------------------------------------------------------------- +logical elemental pure function cNeq(a,b,tol) + + implicit none + complex(pReal), intent(in) :: a,b + real(pReal), intent(in), optional :: tol + real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C + + cNeq = merge(.False., .True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b]))) +end function cNeq + end module prec diff --git a/code/spectral_utilities.f90 b/code/spectral_utilities.f90 index 2979fb52e..1bf322843 100644 --- a/code/spectral_utilities.f90 +++ b/code/spectral_utilities.f90 @@ -402,7 +402,7 @@ subroutine utilities_updateGamma(C,saveReference) integer(pInt) :: & i, j, k, & l, m, n, o - logical :: ierr + logical :: err C_ref = C if (saveReference) then @@ -426,7 +426,7 @@ subroutine utilities_updateGamma(C,saveReference) matA(1:3,1:3) = real(temp33_complex); matA(4:6,4:6) = real(temp33_complex) matA(1:3,4:6) = aimag(temp33_complex); matA(4:6,1:3) = -aimag(temp33_complex) if (abs(math_det33(matA(1:3,1:3))) > 1e-16) then - call math_invert(6_pInt, matA, matInvA, ierr) + call math_invert(6_pInt, matA, matInvA, err) temp33_complex = cmplx(matInvA(1:3,1:3),matInvA(1:3,4:6),pReal) forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt) & gamma_hat(l,m,n,o,i,j,k-grid3Offset) = temp33_complex(l,n)* & @@ -938,6 +938,10 @@ end subroutine utilities_fourierTensorDivergence !-------------------------------------------------------------------------------------------------- subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, & P,C_volAvg,C_minmaxAvg,P_av,forwardData,rotation_BC) + use prec, only: & + dNeq + use IO, only: & + IO_error use debug, only: & debug_reset, & debug_info @@ -976,10 +980,9 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, & age integer(pInt) :: & - j,k + j,k,ierr real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF real(pReal) :: max_dPdF_norm, min_dPdF_norm, defgradDetMin, defgradDetMax, defgradDet - PetscErrorCode :: ierr external :: & MPI_Reduce, & @@ -1013,7 +1016,9 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, & defgradDetMin = min(defgradDetMin,defgradDet) end do call MPI_reduce(MPI_IN_PLACE,defgradDetMax,1,MPI_DOUBLE,MPI_MAX,0,PETSC_COMM_WORLD,ierr) + if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce max') call MPI_reduce(MPI_IN_PLACE,defgradDetMin,1,MPI_DOUBLE,MPI_MIN,0,PETSC_COMM_WORLD,ierr) + if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce min') if (worldrank == 0_pInt) then write(6,'(a,1x,es11.4)') ' max determinant of deformation =', defgradDetMax write(6,'(a,1x,es11.4)') ' min determinant of deformation =', defgradDetMin @@ -1039,7 +1044,9 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, & end do call MPI_Allreduce(MPI_IN_PLACE,max_dPdF,81,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) + if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce max') call MPI_Allreduce(MPI_IN_PLACE,min_dPdF,81,MPI_DOUBLE,MPI_MIN,PETSC_COMM_WORLD,ierr) + if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce min') C_minmaxAvg = 0.5_pReal*(max_dPdF + min_dPdF) C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt @@ -1194,6 +1201,10 @@ end function utilities_getFreqDerivative ! convolution !-------------------------------------------------------------------------------------------------- subroutine utilities_updateIPcoords(F) + use prec, only: & + cNeq + use IO, only: & + IO_error use math, only: & math_mul33x3 use mesh, only: & @@ -1205,10 +1216,9 @@ subroutine utilities_updateIPcoords(F) implicit none real(pReal), dimension(3,3,grid(1),grid(2),grid3), intent(in) :: F - integer(pInt) :: i, j, k, m + integer(pInt) :: i, j, k, m, ierr real(pReal), dimension(3) :: step, offset_coords real(pReal), dimension(3,3) :: Favg - PetscErrorCode :: ierr external & MPI_Bcast @@ -1219,8 +1229,8 @@ subroutine utilities_updateIPcoords(F) call utilities_FFTtensorForward() call utilities_fourierTensorDivergence() - do k = 1_pInt, grid3; do j = 1_pInt, grid(2) ;do i = 1_pInt, grid1Red - if (any(abs(xi1st(1:3,i,j,k)) > tiny(0.0_pReal))) & + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid1Red + if (any(cNeq(xi1st(1:3,i,j,k),cmplx(0.0_pReal,0.0_pReal)))) & vectorField_fourier(1:3,i,j,k) = vectorField_fourier(1:3,i,j,k)/ & sum(conjg(-xi1st(1:3,i,j,k))*xi1st(1:3,i,j,k)) enddo; enddo; enddo @@ -1230,12 +1240,14 @@ subroutine utilities_updateIPcoords(F) ! average F if (grid3Offset == 0_pInt) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt call MPI_Bcast(Favg,9,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr) + if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='update_IPcoords') !-------------------------------------------------------------------------------------------------- ! add average to fluctuation and put (0,0,0) on (0,0,0) step = geomSize/real(grid, pReal) if (grid3Offset == 0_pInt) offset_coords = vectorField_real(1:3,1,1,1) call MPI_Bcast(offset_coords,3,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr) + if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='update_IPcoords') offset_coords = math_mul33x3(Favg,step/2.0_pReal) - offset_coords m = 1_pInt do k = 1_pInt,grid3; do j = 1_pInt,grid(2); do i = 1_pInt,grid(1)