From 6378cb1ff143c92da1433c6669b25512ac824fa9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 29 Oct 2016 09:39:08 +0200 Subject: [PATCH] using improved dEq0/dNeq instead of dEq/dNeq with 0.0 argument --- code/crystallite.f90 | 18 ++++++++--------- code/homogenization_RGC.f90 | 4 ++-- code/plastic_disloUCLA.f90 | 10 +++++----- code/plastic_dislotwin.f90 | 30 ++++++++++++++--------------- code/plastic_isotropic.f90 | 4 ++-- code/plastic_nonlocal.f90 | 35 +++++++++++++++++----------------- code/plastic_phenoplus.f90 | 13 ++++++------- code/plastic_phenopowerlaw.f90 | 13 ++++++------- 8 files changed, 63 insertions(+), 64 deletions(-) diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 2212ba9a1..a570ad714 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -508,7 +508,7 @@ end subroutine crystallite_init subroutine crystallite_stressAndItsTangent(updateJaco) use prec, only: & tol_math_check, & - dNeq + dNeq0 use numerics, only: & subStepMinCryst, & subStepSizeCryst, & @@ -776,7 +776,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) endif else subFracIntermediate = maxval(crystallite_subFrac, mask=.not.crystallite_localPlasticity) - if (dNeq(subFracIntermediate,0.0_pReal)) then + if (dNeq0(subFracIntermediate)) then crystallite_neighborEnforcedCutback = .false. ! look for ips that require a cutback because of a nonconverged neighbor !$OMP PARALLEL !$OMP DO PRIVATE(neighboring_e,neighboring_i) @@ -817,7 +817,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. dNeq(crystallite_subFrac(1,i,e),0.0_pReal)) then + if (.not. crystallite_localPlasticity(1,i,e) .and. dNeq0(crystallite_subFrac(1,i,e))) 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) @@ -3070,7 +3070,7 @@ end subroutine crystallite_integrateStateFPI logical function crystallite_stateJump(ipc,ip,el) use prec, only: & prec_isNaN, & - dNeq + dNeq0 use debug, only: & debug_level, & debug_crystallite, & @@ -3122,7 +3122,7 @@ logical function crystallite_stateJump(ipc,ip,el) enddo #ifndef _OPENMP - if (any(dNeq(plasticState(p)%deltaState(1:mySizePlasticDeltaState,c),0.0_pReal)) & + if (any(dNeq0(plasticState(p)%deltaState(1:mySizePlasticDeltaState,c))) & .and. iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then @@ -3178,7 +3178,7 @@ logical function crystallite_integrateStress(& use prec, only: pLongInt, & tol_math_check, & prec_isNaN, & - dEq + dEq0 use numerics, only: nStress, & aTol_crystalliteStress, & rTol_crystalliteStress, & @@ -3325,7 +3325,7 @@ logical function crystallite_integrateStress(& !* inversion of Fp_current... invFp_current = math_inv33(Fp_current) - failedInversionFp: if (all(dEq(invFp_current,0.0_pReal))) then + failedInversionFp: if (all(dEq0(invFp_current))) 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)') '<< CRYST >> integrateStress failed on inversion of Fp_current at el (elFE) ip g ',& @@ -3341,7 +3341,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 + failedInversionFi: if (all(dEq0(invFi_current))) 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)') '<< CRYST >> integrateStress failed on inversion of Fi_current at el (elFE) ip ipc ',& @@ -3600,7 +3600,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) - failedInversionInvFp: if (all(dEq(Fp_new,0.0_pReal))) then + failedInversionInvFp: if (all(dEq0(Fp_new))) 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 ',& diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index 3b75ec016..ef293fc22 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -387,7 +387,7 @@ end subroutine homogenization_RGC_partitionDeformation !-------------------------------------------------------------------------------------------------- function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) use prec, only: & - dEq + dEq0 use debug, only: & debug_level, & debug_homogenization,& @@ -443,7 +443,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 - zeroTimeStep: if(dEq(dt,0.0_pReal)) then + zeroTimeStep: if(dEq0(dt)) then homogenization_RGC_updateState = .true. ! pretend everything is fine and return return endif zeroTimeStep diff --git a/code/plastic_disloUCLA.f90 b/code/plastic_disloUCLA.f90 index e91883173..147be2f54 100644 --- a/code/plastic_disloUCLA.f90 +++ b/code/plastic_disloUCLA.f90 @@ -974,7 +974,7 @@ end subroutine plastic_disloUCLA_LpAndItsTangent subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,ipc,ip,el) use prec, only: & tol_math_check, & - dEq + dEq0 use math, only: & pi use material, only: & @@ -1112,7 +1112,7 @@ subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,ipc,ip,el) !* Dipole formation EdgeDipMinDistance = & plastic_disloUCLA_CEdgeDipMinDistance(instance)*plastic_disloUCLA_burgersPerSlipSystem(j,instance) - if (dEq(tau_slip_pos,0.0_pReal)) then + if (dEq0(tau_slip_pos)) then DotRhoDipFormation = 0.0_pReal else EdgeDipDistance = & @@ -1140,7 +1140,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 (dEq(tau_slip_pos,0.0_pReal)) then + if (dEq0(tau_slip_pos)) then DotRhoEdgeDipClimb = 0.0_pReal else ClimbVelocity = & @@ -1174,7 +1174,7 @@ end subroutine plastic_disloUCLA_dotState function plastic_disloUCLA_postResults(Tstar_v,Temperature,ipc,ip,el) use prec, only: & tol_math_check, & - dEq + dEq0 use math, only: & pi use material, only: & @@ -1402,7 +1402,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 (dEq(gdot_slip_pos(j)+gdot_slip_neg(j),0.0_pReal)) then + if (dEq0(gdot_slip_pos(j)+gdot_slip_neg(j))) 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 e205f5abb..43fe0c6d8 100644 --- a/code/plastic_dislotwin.f90 +++ b/code/plastic_dislotwin.f90 @@ -200,7 +200,8 @@ 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, & + dEq0, & + dNeq0, & dNeq use debug, only: & debug_level,& @@ -748,8 +749,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 (dEq(plastic_dislotwin_SFE_0K(instance), 0.0_pReal) .and. & - dEq(plastic_dislotwin_dSFE_dT(instance),0.0_pReal) .and. & + if (dEq0(plastic_dislotwin_SFE_0K(instance)) .and. & + dEq0(plastic_dislotwin_dSFE_dT(instance)) .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) & @@ -758,8 +759,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 (dEq(plastic_dislotwin_SFE_0K(instance), 0.0_pReal) .and. & - dEq(plastic_dislotwin_dSFE_dT(instance),0.0_pReal) .and. & + if (dEq0(plastic_dislotwin_SFE_0K(instance)) .and. & + dEq0(plastic_dislotwin_dSFE_dT(instance)) .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) & @@ -772,7 +773,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 (dNeq(plastic_dislotwin_dipoleFormationFactor(instance), 0.0_pReal) .and. & + if (dNeq0(plastic_dislotwin_dipoleFormationFactor(instance)) .and. & dNeq(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. & @@ -1624,7 +1625,7 @@ end subroutine plastic_dislotwin_microstructure subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,ipc,ip,el) use prec, only: & tol_math_check, & - dNeq + dNeq0 use math, only: & math_Plain3333to99, & math_Mandel6to33, & @@ -1771,8 +1772,7 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature !-------------------------------------------------------------------------------------------------- ! Shear banding (shearband) part - if(dNeq(plastic_dislotwin_sbVelocity(instance), 0.0_pReal) .and. & - dNeq(plastic_dislotwin_sbResistance(instance),0.0_pReal)) then + if(dNeq0(plastic_dislotwin_sbVelocity(instance)) .and. dNeq0(plastic_dislotwin_sbResistance(instance))) then gdot_sb = 0.0_pReal dgdot_dtausb = 0.0_pReal call math_eigenValuesVectorsSym(math_Mandel6to33(Tstar_v),eigValues,eigVectors,error) @@ -1939,7 +1939,7 @@ end subroutine plastic_dislotwin_LpAndItsTangent subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el) use prec, only: & tol_math_check, & - dEq + dEq0 use math, only: & pi use material, only: & @@ -2040,7 +2040,7 @@ subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el) !* Dipole formation EdgeDipMinDistance = & plastic_dislotwin_CEdgeDipMinDistance(instance)*plastic_dislotwin_burgersPerSlipSystem(j,instance) - if (dEq(tau_slip(j),0.0_pReal)) then + if (dEq0(tau_slip(j))) then DotRhoDipFormation = 0.0_pReal else EdgeDipDistance = & @@ -2068,10 +2068,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 (dEq(tau_slip(j),0.0_pReal)) then + if (dEq0(tau_slip(j))) then DotRhoEdgeDipClimb = 0.0_pReal else - if (dEq(EdgeDipDistance-EdgeDipMinDistance,0.0_pReal)) then + if (dEq0(EdgeDipDistance-EdgeDipMinDistance)) then DotRhoEdgeDipClimb = 0.0_pReal else ClimbVelocity = 3.0_pReal*lattice_mu(ph)*VacancyDiffusion*AtomicVolume/ & @@ -2187,7 +2187,7 @@ end subroutine plastic_dislotwin_dotState function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el) use prec, only: & tol_math_check, & - dEq + dEq0 use math, only: & pi, & math_Mandel6to33, & @@ -2503,7 +2503,7 @@ function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el) !* Stress exponent plastic_dislotwin_postResults(c+j) = & - merge(0.0_pReal,(tau/gdot_slip(j))*dgdot_dtauslip,dEq(gdot_slip(j),0.0_pReal)) + merge(0.0_pReal,(tau/gdot_slip(j))*dgdot_dtauslip,dEq0(gdot_slip(j))) enddo ; enddo c = c + ns case (sb_eigenvalues_ID) diff --git a/code/plastic_isotropic.f90 b/code/plastic_isotropic.f90 index 9f8129fac..c26b63cfe 100644 --- a/code/plastic_isotropic.f90 +++ b/code/plastic_isotropic.f90 @@ -517,7 +517,7 @@ end subroutine plastic_isotropic_LiAndItsTangent !-------------------------------------------------------------------------------------------------- subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el) use prec, only: & - dEq + dEq0 use math, only: & math_mul6x6 use material, only: & @@ -564,7 +564,7 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el) !-------------------------------------------------------------------------------------------------- ! hardening coefficient if (abs(gamma_dot) > 1e-12_pReal) then - if (dEq(param(instance)%tausat_SinhFitA,0.0_pReal)) then + if (dEq0(param(instance)%tausat_SinhFitA)) then saturation = param(instance)%tausat else saturation = ( param(instance)%tausat & diff --git a/code/plastic_nonlocal.f90 b/code/plastic_nonlocal.f90 index f7b36897a..0c0a48523 100644 --- a/code/plastic_nonlocal.f90 +++ b/code/plastic_nonlocal.f90 @@ -1542,7 +1542,7 @@ end subroutine plastic_nonlocal_aTolState !-------------------------------------------------------------------------------------------------- subroutine plastic_nonlocal_microstructure(Fe, Fp, ip, el) use prec, only: & - dEq + dEq0 use IO, only: & IO_error use math, only: & @@ -1786,7 +1786,7 @@ if (.not. phase_localPlasticity(ph) .and. shortRangeStressCorrection(instance)) - neighbor_rhoExcess(c,s,neighbors(2)) enddo invConnections = math_inv33(connections) - if (all(dEq(invConnections,0.0_pReal))) & + if (all(dEq0(invConnections))) & 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)) @@ -2195,7 +2195,7 @@ end subroutine plastic_nonlocal_LpAndItsTangent !-------------------------------------------------------------------------------------------------- subroutine plastic_nonlocal_deltaState(Tstar_v,ip,el) use prec, only: & - dNeq + dNeq0 use debug, only: debug_level, & debug_constitutive, & debug_levelBasic, & @@ -2318,8 +2318,8 @@ dUpper(1:ns,2) = lattice_mu(ph) * burgers(1:ns,instance) / (4.0_pReal * pi * abs forall (c = 1_pInt:2_pInt) - 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)) & + where(dNeq0(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) = 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)) @@ -2331,7 +2331,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. & - dNeq(dUpperOld(s,c) - dLower(s,c),0.0_pReal)) & + dNeq0(dUpperOld(s,c) - dLower(s,c))) & deltaRhoDipole2SingleStress(s,8_pInt+c) = rhoDip(s,c) * deltaDUpper(s,c) & / (dUpperOld(s,c) - dLower(s,c)) @@ -2379,8 +2379,9 @@ subroutine plastic_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, & timestep,subfrac, ip,el) use prec, only: DAMASK_NaN, & + dNeq0, & dNeq, & - dEq + dEq0 use numerics, only: numerics_integrationMode, & numerics_timeSyncing use IO, only: IO_error @@ -2614,8 +2615,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(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)) & + where(dNeq0(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) = 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)) @@ -2827,11 +2828,11 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then my_rhoSgl = rhoSgl my_v = v if(numerics_timeSyncing) then - if (dEq(subfrac(1_pInt,ip,el),0.0_pReal)) then + if (dEq0(subfrac(1_pInt,ip,el))) then my_rhoSgl = rhoSgl0 my_v = v0 elseif (neighbor_n > 0_pInt) then - if (dEq(subfrac(1_pInt,neighbor_ip,neighbor_el),0.0_pReal)) then + if (dEq0(subfrac(1_pInt,neighbor_ip,neighbor_el))) then my_rhoSgl = rhoSgl0 my_v = v0 endif @@ -3170,7 +3171,7 @@ end subroutine plastic_nonlocal_updateCompatibility !********************************************************************* function plastic_nonlocal_dislocationstress(Fe, ip, el) use prec, only: & - dEq + dEq0 use math, only: math_mul33x33, & math_mul33x3, & math_inv33, & @@ -3383,7 +3384,7 @@ if (.not. phase_localPlasticity(ph)) then Rsquare = R * R Rcube = Rsquare * R denominator = R * (R + flipSign * lambda) - if (dEq(denominator,0.0_pReal)) exit ipLoop + if (dEq0(denominator)) exit ipLoop sigma(1,1) = sigma(1,1) - real(side,pReal) & * flipSign * z / denominator & @@ -3428,7 +3429,7 @@ if (.not. phase_localPlasticity(ph)) then Rsquare = R * R Rcube = Rsquare * R denominator = R * (R + flipSign * lambda) - if (dEq(denominator,0.0_pReal)) exit ipLoop + if (dEq0(denominator)) exit ipLoop sigma(1,2) = sigma(1,2) - real(side,pReal) * flipSign * z & * (1.0_pReal - lattice_nu(ph)) / denominator & @@ -3518,7 +3519,7 @@ end function plastic_nonlocal_dislocationstress !-------------------------------------------------------------------------------------------------- function plastic_nonlocal_postResults(Tstar_v,Fe,ip,el) use prec, only: & - dNeq + dNeq0 use math, only: & math_mul6x6, & math_mul33x3, & @@ -3635,8 +3636,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(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)) & + where(dNeq0(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) = 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 b9d09fb3d..105a64b34 100644 --- a/code/plastic_phenoplus.f90 +++ b/code/plastic_phenoplus.f90 @@ -113,7 +113,7 @@ 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 + dEq0 use debug, only: & debug_level, & debug_constitutive,& @@ -479,8 +479,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(dEq(plastic_phenoplus_a_slip(instance),0.0_pReal) .and. & - plastic_phenoplus_Nslip(:,instance) > 0)) & + if (any(dEq0(plastic_phenoplus_a_slip(instance)) .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. & plastic_phenoplus_Ntwin(:,instance) > 0)) & @@ -918,7 +917,7 @@ end subroutine plastic_phenoplus_microstructure !-------------------------------------------------------------------------------------------------- subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) use prec, only: & - dNeq + dNeq0 use math, only: & math_Plain3333to99, & math_Mandel6to33 @@ -1034,7 +1033,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 (dNeq(gdot_slip_pos,0.0_pReal)) then + if (dNeq0(gdot_slip_pos)) 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) + & @@ -1042,7 +1041,7 @@ subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) nonSchmid_tensor(m,n,1) endif - if (dNeq(gdot_slip_neg,0.0_pReal)) then + if (dNeq0(gdot_slip_neg)) 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) + & @@ -1069,7 +1068,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 (dNeq(gdot_twin,0.0_pReal)) then + if (dNeq0(gdot_twin)) 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 0db6a6072..a7c7b10e6 100644 --- a/code/plastic_phenopowerlaw.f90 +++ b/code/plastic_phenopowerlaw.f90 @@ -125,7 +125,7 @@ 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 + dEq0 use debug, only: & debug_level, & debug_constitutive,& @@ -485,8 +485,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(dEq(plastic_phenopowerlaw_a_slip(instance),0.0_pReal) .and. & - plastic_phenopowerlaw_Nslip(:,instance) > 0)) & + if (any(dEq0(plastic_phenopowerlaw_a_slip(instance)) .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. & plastic_phenopowerlaw_Ntwin(:,instance) > 0)) & @@ -769,7 +768,7 @@ end subroutine plastic_phenopowerlaw_aTolState !-------------------------------------------------------------------------------------------------- subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) use prec, only: & - dNeq + dNeq0 use math, only: & math_Plain3333to99, & math_Mandel6to33 @@ -859,7 +858,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 (dNeq(gdot_slip_pos,0.0_pReal)) then + if (dNeq0(gdot_slip_pos)) 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) + & @@ -867,7 +866,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip, nonSchmid_tensor(m,n,1) endif - if (dNeq(gdot_slip_neg,0.0_pReal)) then + if (dNeq0(gdot_slip_neg)) 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) + & @@ -894,7 +893,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 (dNeq(gdot_twin,0.0_pReal)) then + if (dNeq0(gdot_twin)) 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) + &