From b35ff67f99d56aeda4605e225f408c06aabe55c1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 4 May 2017 00:32:44 +0200 Subject: [PATCH 01/62] using IEEE_is_NAN and IEEE_quiet_NaN instead of hand-written solution, will not work for gfortran < 5 --- src/crystallite.f90 | 120 +++++++++++++++++++------------------ src/lattice.f90 | 6 +- src/math.f90 | 36 ++++------- src/plastic_isotropic.f90 | 17 +++--- src/plastic_nonlocal.f90 | 10 ++-- src/prec.f90 | 36 +---------- src/spectral_utilities.f90 | 6 +- 7 files changed, 92 insertions(+), 139 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index efeab9dcd..8d2ae5f23 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1232,8 +1232,8 @@ end subroutine crystallite_stressAndItsTangent !> @brief integrate stress, state with 4th order explicit Runge Kutta method !-------------------------------------------------------------------------------------------------- subroutine crystallite_integrateStateRK4() - use prec, only: & - prec_isNaN + use, intrinsic :: & + IEEE_arithmetic use numerics, only: & numerics_integrationMode use debug, only: & @@ -1316,7 +1316,7 @@ subroutine crystallite_integrateStateRK4() ! first Runge-Kutta step !$OMP PARALLEL !$OMP DO - do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) & call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, & @@ -1326,22 +1326,22 @@ subroutine crystallite_integrateStateRK4() !$OMP ENDDO !$OMP DO PRIVATE(p,c,NaN) - do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then c = phasememberAt(g,i,e) p = phaseAt(g,i,e) - NaN = any(prec_isNaN(plasticState(p)%dotState(:,c))) + NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c))) do mySource = 1_pInt, phase_Nsources(p) - NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c))) + NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c))) enddo - if (NaN) then ! NaN occured in any dotState - if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... + if (NaN) then ! NaN occured in any dotState + if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) - crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped + crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped !$OMP END CRITICAL (checkTodo) - else ! if broken local... - crystallite_todo(g,i,e) = .false. ! ... skip this one next time + else ! if broken local... + crystallite_todo(g,i,e) = .false. ! ... skip this one next time endif endif endif @@ -1356,7 +1356,7 @@ subroutine crystallite_integrateStateRK4() !$OMP PARALLEL !$OMP DO PRIVATE(p,c) - do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then p = phaseAt(g,i,e) c = phasememberAt(g,i,e) @@ -1371,7 +1371,7 @@ subroutine crystallite_integrateStateRK4() !$OMP ENDDO !$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,c) - do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then p = phaseAt(g,i,e) @@ -1475,9 +1475,9 @@ subroutine crystallite_integrateStateRK4() p = phaseAt(g,i,e) c = phasememberAt(g,i,e) - NaN = any(prec_isNaN(plasticState(p)%dotState(:,c))) + NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c))) do mySource = 1_pInt, phase_Nsources(p) - NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c))) + NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c))) enddo if (NaN) then ! NaN occured in any dotState if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... @@ -1528,8 +1528,8 @@ end subroutine crystallite_integrateStateRK4 !> adaptive step size (use 5th order solution to advance = "local extrapolation") !-------------------------------------------------------------------------------------------------- subroutine crystallite_integrateStateRKCK45() - use prec, only: & - prec_isNaN + use, intrinsic :: & + IEEE_arithmetic use debug, only: & debug_level, & debug_crystallite, & @@ -1647,9 +1647,9 @@ subroutine crystallite_integrateStateRKCK45() if (crystallite_todo(g,i,e)) then cc = phasememberAt(g,i,e) p = phaseAt(g,i,e) - NaN = any(prec_isNaN(plasticState(p)%dotState(:,cc))) + NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,cc))) do mySource = 1_pInt, phase_Nsources(p) - NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,cc))) + NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,cc))) enddo if (NaN) then ! NaN occured in any dotState if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... @@ -1801,9 +1801,9 @@ subroutine crystallite_integrateStateRKCK45() p = phaseAt(g,i,e) cc = phasememberAt(g,i,e) - NaN = any(prec_isNaN(plasticState(p)%dotState(:,cc))) + NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,cc))) do mySource = 1_pInt, phase_Nsources(p) - NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,cc))) + NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,cc))) enddo if (NaN) then ! NaN occured in any dotState if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... @@ -2031,8 +2031,8 @@ end subroutine crystallite_integrateStateRKCK45 !> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- subroutine crystallite_integrateStateAdaptiveEuler() - use prec, only: & - prec_isNaN + use, intrinsic :: & + IEEE_arithmetic use debug, only: & debug_level, & debug_crystallite, & @@ -2133,9 +2133,9 @@ subroutine crystallite_integrateStateAdaptiveEuler() if (crystallite_todo(g,i,e)) then p = phaseAt(g,i,e) c = phasememberAt(g,i,e) - NaN = any(prec_isNaN(plasticState(p)%dotState(:,c))) + NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c))) do mySource = 1_pInt, phase_Nsources(p) - NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c))) + NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c))) enddo if (NaN) then ! NaN occured in any dotState if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... @@ -2254,9 +2254,9 @@ subroutine crystallite_integrateStateAdaptiveEuler() if (crystallite_todo(g,i,e)) then p = phaseAt(g,i,e) c = phasememberAt(g,i,e) - NaN = any(prec_isNaN(plasticState(p)%dotState(:,c))) + NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c))) do mySource = 1_pInt, phase_Nsources(p) - NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c))) + NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c))) enddo if (NaN) then ! NaN occured in any dotState if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... @@ -2391,8 +2391,8 @@ end subroutine crystallite_integrateStateAdaptiveEuler !> @brief integrate stress, and state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- subroutine crystallite_integrateStateEuler() - use prec, only: & - prec_isNaN + use, intrinsic :: & + IEEE_arithmetic use debug, only: & debug_level, & debug_crystallite, & @@ -2471,9 +2471,9 @@ eIter = FEsolving_execElem(1:2) if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then c = phasememberAt(g,i,e) p = phaseAt(g,i,e) - NaN = any(prec_isNaN(plasticState(p)%dotState(:,c))) + NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c))) do mySource = 1_pInt, phase_Nsources(p) - NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c))) + NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c))) enddo if (NaN) then ! NaN occured in any dotState if (.not. crystallite_localPlasticity(g,i,e) .and. .not. numerics_timeSyncing) then ! if broken non-local... @@ -2614,8 +2614,8 @@ end subroutine crystallite_integrateStateEuler !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- subroutine crystallite_integrateStateFPI() - use prec, only: & - prec_isNaN + use, intrinsic :: & + IEEE_arithmetic use debug, only: & debug_e, & debug_i, & @@ -2732,22 +2732,22 @@ subroutine crystallite_integrateStateFPI() !$OMP ENDDO !$OMP DO PRIVATE(p,c,NaN) - do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then p = phaseAt(g,i,e) c = phasememberAt(g,i,e) - NaN = any(prec_isNaN(plasticState(p)%dotState(:,c))) + NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c))) do mySource = 1_pInt, phase_Nsources(p) - NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c))) + NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c))) enddo - if (NaN) then ! NaN occured in any dotState - if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a non-local... + if (NaN) then ! NaN occured in any dotState + if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a non-local... !$OMP CRITICAL (checkTodo) - crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals done (and broken) + crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals done (and broken) !$OMP END CRITICAL (checkTodo) - else ! broken one was local... - crystallite_todo(g,i,e) = .false. ! ... done (and broken) + else ! broken one was local... + crystallite_todo(g,i,e) = .false. ! ... done (and broken) endif endif endif @@ -2851,9 +2851,9 @@ subroutine crystallite_integrateStateFPI() if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then p = phaseAt(g,i,e) c = phasememberAt(g,i,e) - NaN = any(prec_isNaN(plasticState(p)%dotState(:,c))) + NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c))) do mySource = 1_pInt, phase_Nsources(p) - NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c))) + NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c))) enddo if (NaN) then ! NaN occured in any dotState crystallite_todo(g,i,e) = .false. ! ... skip me next time @@ -3062,8 +3062,9 @@ end subroutine crystallite_integrateStateFPI !> returns true, if state jump was successfull or not needed. false indicates NaN in delta state !-------------------------------------------------------------------------------------------------- logical function crystallite_stateJump(ipc,ip,el) + use, intrinsic :: & + IEEE_arithmetic use prec, only: & - prec_isNaN, & dNeq0 use debug, only: & debug_level, & @@ -3098,7 +3099,7 @@ logical function crystallite_stateJump(ipc,ip,el) p = phaseAt(ipc,ip,el) call constitutive_collectDeltaState(crystallite_Tstar_v(1:6,ipc,ip,el), crystallite_Fe(1:3,1:3,ipc,ip,el), ipc,ip,el) mySizePlasticDeltaState = plasticState(p)%sizeDeltaState - if( any(prec_isNaN(plasticState(p)%deltaState(:,c)))) then ! NaN occured in deltaState + if( any(IEEE_is_NaN(plasticState(p)%deltaState(:,c)))) then ! NaN occured in deltaState crystallite_stateJump = .false. return endif @@ -3106,7 +3107,7 @@ logical function crystallite_stateJump(ipc,ip,el) plasticState(p)%deltaState(1:mySizePlasticDeltaState,c) do mySource = 1_pInt, phase_Nsources(p) mySizeSourceDeltaState = sourceState(p)%p(mySource)%sizeDeltaState - if( any(prec_isNaN(sourceState(p)%p(mySource)%deltaState(:,c)))) then ! NaN occured in deltaState + if( any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(:,c)))) then ! NaN occured in deltaState crystallite_stateJump = .false. return endif @@ -3169,9 +3170,10 @@ logical function crystallite_integrateStress(& el,& ! element number timeFraction & ) + use, intrinsic :: & + IEEE_arithmetic use prec, only: pLongInt, & tol_math_check, & - prec_isNaN, & dEq0 use numerics, only: nStress, & aTol_crystalliteStress, & @@ -3426,11 +3428,11 @@ logical function crystallite_integrateStress(& !* update current residuum and check for convergence of loop - aTolLp = max(rTol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error - aTol_crystalliteStress) ! minimum lower cutoff + aTolLp = max(rTol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error + aTol_crystalliteStress) ! minimum lower cutoff residuumLp = Lpguess - Lp_constitutive - if (any(prec_isNaN(residuumLp))) then ! NaN in residuum... + if (any(IEEE_is_NaN(residuumLp))) then ! NaN in residuum... #ifndef _OPENMP if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) & write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3,a)') '<< CRYST >> integrateStress encountered NaN at el (elFE) ip ipc ', & @@ -3438,16 +3440,16 @@ logical function crystallite_integrateStress(& ' ; iteration ', NiterationStressLp,& ' >> returning..!' #endif - return ! ...me = .false. to inform integrator about problem - elseif (norm2(residuumLp) < aTolLp) then ! converged if below absolute tolerance - exit LpLoop ! ...leave iteration loop + return ! ...me = .false. to inform integrator about problem + elseif (norm2(residuumLp) < aTolLp) then ! converged if below absolute tolerance + exit LpLoop ! ...leave iteration loop elseif ( NiterationStressLp == 1_pInt & - .or. norm2(residuumLp) < norm2(residuumLp_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)... - residuumLp_old = residuumLp ! ...remember old values and... + .or. norm2(residuumLp) < norm2(residuumLp_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)... + residuumLp_old = residuumLp ! ...remember old values and... Lpguess_old = Lpguess - steplengthLp = steplengthLp0 ! ...proceed with normal step length (calculate new search direction) - else ! not converged and residuum not improved... - steplengthLp = 0.5_pReal * steplengthLp ! ...try with smaller step length in same direction + steplengthLp = steplengthLp0 ! ...proceed with normal step length (calculate new search direction) + else ! not converged and residuum not improved... + steplengthLp = 0.5_pReal * steplengthLp ! ...try with smaller step length in same direction Lpguess = Lpguess_old + steplengthLp * deltaLp cycle LpLoop endif @@ -3520,7 +3522,7 @@ logical function crystallite_integrateStress(& aTolLi = max(rTol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error aTol_crystalliteStress) ! minimum lower cutoff residuumLi = Liguess - Li_constitutive - if (any(prec_isNaN(residuumLi))) then ! NaN in residuum... + if (any(IEEE_is_NaN(residuumLi))) then ! NaN in residuum... return ! ...me = .false. to inform integrator about problem elseif (norm2(residuumLi) < aTolLi) then ! converged if below absolute tolerance exit LiLoop ! ...leave iteration loop diff --git a/src/lattice.f90 b/src/lattice.f90 index 7726c772d..a970ed85a 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2170,8 +2170,8 @@ pure function lattice_symmetrize33(struct,T33) !> @brief figures whether unit quat falls into stereographic standard triangle !-------------------------------------------------------------------------------------------------- logical pure function lattice_qInSST(Q, struct) - use prec, only: & - prec_isNaN + use, intrinsic :: & + IEEE_arithmetic use math, only: & math_qToRodrig @@ -2181,7 +2181,7 @@ logical pure function lattice_qInSST(Q, struct) real(pReal), dimension(3) :: Rodrig ! Rodrigues vector of Q Rodrig = math_qToRodrig(Q) - if (any(prec_isNaN(Rodrig))) then + if (any(IEEE_is_NaN(Rodrig))) then lattice_qInSST = .false. else select case (struct) diff --git a/src/math.f90 b/src/math.f90 index 9629401ca..3a0533ddf 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1666,15 +1666,16 @@ end function math_qToEulerAxisAngle !> @brief Rodrigues vector (x, y, z) from unit quaternion (w+ix+jy+kz) !-------------------------------------------------------------------------------------------------- pure function math_qToRodrig(Q) + use, intrinsic :: & + IEEE_arithmetic use prec, only: & - DAMASK_NaN, & tol_math_check implicit none real(pReal), dimension(4), intent(in) :: Q real(pReal), dimension(3) :: math_qToRodrig - math_qToRodrig = merge(Q(2:4)/Q(1),DAMASK_NaN,abs(Q(1)) > tol_math_check) ! NaN for 180 deg since Rodrig is unbound + math_qToRodrig = merge(Q(2:4)/Q(1),IEEE_value(1.0_pReal,IEEE_quiet_NaN),abs(Q(1)) > tol_math_check)! NaN for 180 deg since Rodrig is unbound end function math_qToRodrig @@ -2092,11 +2093,11 @@ end function math_rotationalPart33 !-------------------------------------------------------------------------------------------------- !> @brief Eigenvalues of symmetric matrix m -! will return NaN on error +! will return NaN on error !-------------------------------------------------------------------------------------------------- function math_eigenvaluesSym(m) - use prec, only: & - DAMASK_NaN + use, intrinsic :: & + IEEE_arithmetic implicit none real(pReal), dimension(:,:), intent(in) :: m @@ -2109,7 +2110,7 @@ function math_eigenvaluesSym(m) vectors = m ! copy matrix to input (doubles as output) array call dsyev('N','U',size(m,1),vectors,size(m,1),math_eigenvaluesSym,work,(64+2)*size(m,1),info) - if (info /= 0_pInt) math_eigenvaluesSym = DAMASK_NaN + if (info /= 0_pInt) math_eigenvaluesSym = IEEE_value(1.0_pReal,IEEE_quiet_NaN) end function math_eigenvaluesSym @@ -2701,29 +2702,13 @@ pure function math_rotate_forward3333(tensor,rot_tensor) end function math_rotate_forward3333 -!-------------------------------------------------------------------------------------------------- -!> @brief calculate average of tensor field -!-------------------------------------------------------------------------------------------------- -function math_tensorAvg(field) - - implicit none - real(pReal), dimension(3,3) :: math_tensorAvg - real(pReal), intent(in), dimension(:,:,:,:,:) :: field - real(pReal) :: wgt - - wgt = 1.0_pReal/real(size(field,3)*size(field,4)*size(field,5), pReal) - math_tensorAvg = sum(sum(sum(field,dim=5),dim=4),dim=3)*wgt - -end function math_tensorAvg - - !-------------------------------------------------------------------------------------------------- !> @brief limits a scalar value to a certain range (either one or two sided) ! Will return NaN if left > right !-------------------------------------------------------------------------------------------------- real(pReal) pure function math_limit(a, left, right) - use prec, only: & - DAMASK_NaN + use, intrinsic :: & + IEEE_arithmetic implicit none real(pReal), intent(in) :: a @@ -2735,7 +2720,8 @@ real(pReal) pure function math_limit(a, left, right) merge(right, huge(a), present(right)) & ) - if (present(left) .and. present(right)) math_limit = merge (DAMASK_NaN,math_limit, left>right) + if (present(left) .and. present(right)) & + math_limit = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_limit, left>right) end function math_limit diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90 index 5e3722e6a..def4c5a35 100644 --- a/src/plastic_isotropic.f90 +++ b/src/plastic_isotropic.f90 @@ -9,8 +9,7 @@ module plastic_isotropic use prec, only: & pReal,& - pInt, & - DAMASK_NaN + pInt implicit none private @@ -36,14 +35,14 @@ module plastic_isotropic integer(kind(undefined_ID)), allocatable, dimension(:) :: & outputID real(pReal) :: & - fTaylor = DAMASK_NaN, & - tau0 = DAMASK_NaN, & - gdot0 = DAMASK_NaN, & - n = DAMASK_NaN, & - h0 = DAMASK_NaN, & + fTaylor, & + tau0, & + gdot0, & + n, & + h0, & h0_slopeLnRate = 0.0_pReal, & - tausat = DAMASK_NaN, & - a = DAMASK_NaN, & + tausat, & + a, & aTolFlowstress = 1.0_pReal, & aTolShear = 1.0e-6_pReal, & tausat_SinhFitA= 0.0_pReal, & diff --git a/src/plastic_nonlocal.f90 b/src/plastic_nonlocal.f90 index 0c0a48523..0a8c4c3f9 100644 --- a/src/plastic_nonlocal.f90 +++ b/src/plastic_nonlocal.f90 @@ -2377,9 +2377,9 @@ end subroutine plastic_nonlocal_deltaState !--------------------------------------------------------------------------------------------------- subroutine plastic_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, & timestep,subfrac, ip,el) - -use prec, only: DAMASK_NaN, & - dNeq0, & +use, intrinsic :: & + IEEE_arithmetic +use prec, only: dNeq0, & dNeq, & dEq0 use numerics, only: numerics_integrationMode, & @@ -2701,7 +2701,7 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then write(6,'(a)') '<< CONST >> enforcing cutback !!!' endif #endif - plasticState(p)%dotState = DAMASK_NaN ! -> return NaN and, hence, enforce cutback + plasticState(p)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) ! -> return NaN and, hence, enforce cutback return endif @@ -2984,7 +2984,7 @@ if ( any(rhoSglOriginal(1:ns,1:4) + rhoDot(1:ns,1:4) * timestep < -aTolRho(in write(6,'(a)') '<< CONST >> enforcing cutback !!!' endif #endif - plasticState(p)%dotState = DAMASK_NaN + plasticState(p)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) return else forall (s = 1:ns, t = 1_pInt:4_pInt) diff --git a/src/prec.f90 b/src/prec.f90 index c6ae656b8..671e15990 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -5,22 +5,12 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH !> @brief setting precision for real and int type -!> @details setting precision for real and int type and for DAMASK_NaN. Definition is made -!! depending on makro "INT" defined during compilation -!! for details on NaN see https://software.intel.com/en-us/forums/topic/294680 !-------------------------------------------------------------------------------------------------- module prec - -#if !(defined(__GFORTRAN__) && __GNUC__ < 5) - use, intrinsic :: & ! unfortunately not avialable in gfortran <= 5 - IEEE_arithmetic -#endif - implicit none private #if (FLOAT==8) integer, parameter, public :: pReal = 8 !< floating point double precision (was selected_real_kind(15,300), number with 15 significant digits, up to 1e+-300) - real(pReal), parameter, public :: DAMASK_NaN = real(Z'7FF8000000000000',pReal) !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html) #else NO SUITABLE PRECISION FOR REAL SELECTED, STOPPING COMPILATION #endif @@ -106,7 +96,6 @@ module prec public :: & prec_init, & - prec_isNaN, & dEq, & dEq0, & cEq, & @@ -118,7 +107,7 @@ contains !-------------------------------------------------------------------------------------------------- -!> @brief reporting precision and checking if DAMASK_NaN is set correctly +!> @brief reporting precision !-------------------------------------------------------------------------------------------------- subroutine prec_init use, intrinsic :: & @@ -133,36 +122,13 @@ subroutine prec_init write(6,'(a,i3)') ' Bytes for pReal: ',pReal write(6,'(a,i3)') ' Bytes for pInt: ',pInt write(6,'(a,i3)') ' Bytes for pLongInt: ',pLongInt - write(6,'(a,e10.3)') ' NaN: ', DAMASK_NaN - write(6,'(a,l3)') ' NaN != NaN: ',DAMASK_NaN /= DAMASK_NaN - write(6,'(a,l3,/)') ' NaN check passed ',prec_isNAN(DAMASK_NaN) - if ((.not. prec_isNaN(DAMASK_NaN)) .or. (DAMASK_NaN == DAMASK_NaN)) call quit(9000) realloc_lhs_test = [1_pInt,2_pInt] if (realloc_lhs_test(2)/=2_pInt) call quit(9000) - end subroutine prec_init -!-------------------------------------------------------------------------------------------------- -!> @brief figures out if a floating point number is NaN -! basically just a small wrapper, because gfortran < 5.0 does not have the IEEE module -!-------------------------------------------------------------------------------------------------- -logical elemental pure function prec_isNaN(a) - - implicit none - real(pReal), intent(in) :: a - -#if (defined(__GFORTRAN__) && __GNUC__ < 5) - intrinsic :: isNaN - prec_isNaN = isNaN(a) -#else - prec_isNaN = IEEE_is_NaN(a) -#endif -end function prec_isNaN - - !-------------------------------------------------------------------------------------------------- !> @brief equality comparison for float with double precision ! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90 index a33c7a9f5..0773a9065 100644 --- a/src/spectral_utilities.f90 +++ b/src/spectral_utilities.f90 @@ -741,8 +741,8 @@ end function utilities_curlRMS !> @brief calculates mask compliance tensor used to adjust F to fullfill stress BC !-------------------------------------------------------------------------------------------------- function utilities_maskedCompliance(rot_BC,mask_stress,C) - use prec, only: & - prec_isNaN + use, intrinsic :: & + IEEE_arithmetic use IO, only: & IO_error use math, only: & @@ -794,7 +794,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) endif; enddo; endif; enddo call math_invert(size_reduced, c_reduced, s_reduced, errmatinv) ! invert reduced stiffness - if (any(prec_isNaN(s_reduced))) errmatinv = .true. + if (any(IEEE_is_NaN(s_reduced))) errmatinv = .true. if(errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance') temp99_Real = 0.0_pReal ! fill up compliance with zeros k = 0_pInt From f52e62e6b9d5c67d62d6e6a19e9f8bb76b397c96 Mon Sep 17 00:00:00 2001 From: Test User Date: Tue, 16 May 2017 20:39:45 +0200 Subject: [PATCH 02/62] [skip ci] updated version information after successful test of v2.0.1-734-g64fde94 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 5a76f3ec4..8484dacfc 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-728-g66efd90 +v2.0.1-734-g64fde94 From 60b7870abea949456f8706f74c0a06586ff7647e Mon Sep 17 00:00:00 2001 From: Test User Date: Tue, 16 May 2017 23:00:49 +0200 Subject: [PATCH 03/62] [skip ci] updated version information after successful test of v2.0.1-735-gf52e62e --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 8484dacfc..5dabd576c 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-734-g64fde94 +v2.0.1-735-gf52e62e From eff2b4ca1cb92221757d3c732e165352f25dd019 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 17 May 2017 08:53:37 +0200 Subject: [PATCH 04/62] scripts now work if sourced from new location --- DAMASK_env.csh | 2 +- DAMASK_env.sh | 2 +- DAMASK_env.zsh | 2 +- env/{DAMASK_env.csh => DAMASK.csh} | 21 +++++++++++++++------ env/{DAMASK_env.sh => DAMASK.sh} | 27 +++++++++++++-------------- env/{DAMASK_env.zsh => DAMASK.zsh} | 20 ++++++++++---------- 6 files changed, 41 insertions(+), 33 deletions(-) rename env/{DAMASK_env.csh => DAMASK.csh} (84%) rename env/{DAMASK_env.sh => DAMASK.sh} (88%) rename env/{DAMASK_env.zsh => DAMASK.zsh} (89%) diff --git a/DAMASK_env.csh b/DAMASK_env.csh index 4535015d8..73a008b34 120000 --- a/DAMASK_env.csh +++ b/DAMASK_env.csh @@ -1 +1 @@ -env/DAMASK_env.csh \ No newline at end of file +./env/DAMASK.csh \ No newline at end of file diff --git a/DAMASK_env.sh b/DAMASK_env.sh index e49b7df7c..2ec9efd4a 120000 --- a/DAMASK_env.sh +++ b/DAMASK_env.sh @@ -1 +1 @@ -env/DAMASK_env.sh \ No newline at end of file +./env/DAMASK.sh \ No newline at end of file diff --git a/DAMASK_env.zsh b/DAMASK_env.zsh index 22da3bd4d..7af0ffdbc 120000 --- a/DAMASK_env.zsh +++ b/DAMASK_env.zsh @@ -1 +1 @@ -env/DAMASK_env.zsh \ No newline at end of file +./env/DAMASK.zsh \ No newline at end of file diff --git a/env/DAMASK_env.csh b/env/DAMASK.csh similarity index 84% rename from env/DAMASK_env.csh rename to env/DAMASK.csh index 81d4de421..819b263ef 100644 --- a/env/DAMASK_env.csh +++ b/env/DAMASK.csh @@ -3,16 +3,21 @@ set CALLED=($_) set DIRNAME=`dirname $CALLED[2]` + +# transition compatibility (renamed $DAMASK_ROOT/DAMASK_env.csh to $DAMASK_ROOT/env/DAMASK.csh) +set FILENAME=`basename $CALLED[2]` +if ($FILENAME == "DAMASK.csh") then + set DIRNAME=$DIRNAME"/../" +endif + set DAMASK_ROOT=`python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" $DIRNAME` + source $DAMASK_ROOT/CONFIG -# if DAMASK_BIN is present and not in $PATH, add it +# if DAMASK_BIN is present if ( $?DAMASK_BIN) then - set MATCH=`echo :${PATH}: | grep ${DAMASK_BIN}:` - if ( "x$MATCH" == "x" ) then - set PATH=${DAMASK_BIN}:${PATH} - endif + set path = ($DAMASK_BIN $path) endif set SOLVER=`which DAMASK_spectral` @@ -62,4 +67,8 @@ if ( $?prompt ) then endif setenv DAMASK_NUM_THREADS $DAMASK_NUM_THREADS -setenv PYTHONPATH $DAMASK_ROOT/lib:$PYTHONPATH +if ( ! $?PYTHONPATH ) then + setenv PYTHONPATH $DAMASK_ROOT/lib +else + setenv PYTHONPATH $DAMASK_ROOT/lib:$PYTHONPATH +endif diff --git a/env/DAMASK_env.sh b/env/DAMASK.sh similarity index 88% rename from env/DAMASK_env.sh rename to env/DAMASK.sh index 264ae9c94..2674505f7 100644 --- a/env/DAMASK_env.sh +++ b/env/DAMASK.sh @@ -1,5 +1,5 @@ # sets up an environment for DAMASK on bash -# usage: source DAMASK_env.sh +# usage: source DAMASK.sh if [ "$OSTYPE" == "linux-gnu" ] || [ "$OSTYPE" == 'linux' ]; then @@ -10,6 +10,11 @@ else DAMASK_ROOT=${STAT##* } fi +# transition compatibility (renamed $DAMASK_ROOT/DAMASK_env.sh to $DAMASK_ROOT/env/DAMASK.sh) +if [ ${BASH_SOURCE##*/} == "DAMASK.sh" ]; then + DAMASK_ROOT=$DAMASK_ROOT'/..' +fi + # shorthand command to change to DAMASK_ROOT directory eval "function damask() { cd $DAMASK_ROOT; }" @@ -20,22 +25,16 @@ set() { source $DAMASK_ROOT/CONFIG unset -f set -# add DAMASK_BIN if present but not yet in $PATH -if [[ "x$DAMASK_BIN" != "x" && ! $(echo ":$PATH:" | grep $DAMASK_BIN:) ]]; then - export PATH=$DAMASK_BIN:$PATH -fi +# add DAMASK_BIN if present +[ "x$DAMASK_BIN" != "x" ] && PATH=$DAMASK_BIN:$PATH SOLVER=$(which DAMASK_spectral || true 2>/dev/null) -if [ "x$SOLVER" == "x" ]; then - SOLVER='Not found!' -fi +[ "x$SOLVER" == "x" ] && SOLVER='Not found!' + PROCESSING=$(which postResults || true 2>/dev/null) -if [ "x$PROCESSING" == "x" ]; then - PROCESSING='Not found!' -fi -if [ "x$DAMASK_NUM_THREADS" == "x" ]; then - DAMASK_NUM_THREADS=1 -fi +[ "x$PROCESSING" == "x" ] && PROCESSING='Not found!' + +[ "x$DAMASK_NUM_THREADS" == "x" ] && DAMASK_NUM_THREADS=1 # according to http://software.intel.com/en-us/forums/topic/501500 # this seems to make sense for the stack size diff --git a/env/DAMASK_env.zsh b/env/DAMASK.zsh similarity index 89% rename from env/DAMASK_env.zsh rename to env/DAMASK.zsh index 9098352bf..23b057d68 100644 --- a/env/DAMASK_env.zsh +++ b/env/DAMASK.zsh @@ -1,7 +1,12 @@ # sets up an environment for DAMASK on zsh -# usage: source DAMASK_env.zsh +# usage: source DAMASK.zsh -DAMASK_ROOT=${0:a:h} +# transition compatibility (renamed $DAMASK_ROOT/DAMASK_env.zsh to $DAMASK_ROOT/env/DAMASK.zsh) +if [ ${0:t:r} = 'DAMASK' ]; then + DAMASK_ROOT=${0:a:h}'/..' +else + DAMASK_ROOT=${0:a:h} +fi # shorthand command to change to DAMASK_ROOT directory eval "function damask() { cd $DAMASK_ROOT; }" @@ -13,17 +18,12 @@ set() { source $DAMASK_ROOT/CONFIG unset -f set -# add DAMASK_BIN if present but not yet in $PATH -MATCH=`echo ":$PATH:" | grep $DAMASK_BIN:` -if [[ ( "x$DAMASK_BIN" != "x" ) && ( "x$MATCH" = "x" ) ]]; then - export PATH=$DAMASK_BIN:$PATH -fi +# add DAMASK_BIN if present +[ ( "x$DAMASK_BIN" != "x" ) ] && PATH=$DAMASK_BIN:$PATH SOLVER=`which DAMASK_spectral || True 2>/dev/null` PROCESSING=`which postResults || True 2>/dev/null` -if [ "x$DAMASK_NUM_THREADS" = "x" ]; then - DAMASK_NUM_THREADS=1 -fi +[ "x$DAMASK_NUM_THREADS" = "x" ] && DAMASK_NUM_THREADS=1 # according to http://software.intel.com/en-us/forums/topic/501500 # this seems to make sense for the stack size From ee8afc33cf7600ff63d720698a8c6a7154fed805 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 17 May 2017 10:15:30 +0200 Subject: [PATCH 05/62] following yamllint suggestions --- .gitlab-ci.yml | 53 ++++++++++++++++++++++++++++---------------------- 1 file changed, 30 insertions(+), 23 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 709655863..cc541307f 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -1,3 +1,4 @@ +--- stages: - prepareAll - preprocessing @@ -23,26 +24,30 @@ stages: ################################################################################################### before_script: - - if [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue)x == 'x' ]; then echo $CI_PIPELINE_ID >> $TESTROOT/GitLabCI.queue; fi - - while [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue) -ne 1 ];do sleep 5m; done + - if [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue)x == 'x' ] + then echo $CI_PIPELINE_ID >> $TESTROOT/GitLabCI.queue + fi + - while [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue) != 1 ] + do sleep 5m + done - source $DAMASKROOT/DAMASK_env.sh - cd $DAMASKROOT/PRIVATE/testing ################################################################################################### variables: - #================================================================================================ + # =============================================================================================== # GitLab Settings - #================================================================================================ + # =============================================================================================== GIT_SUBMODULE_STRATEGY: none - #================================================================================================ + # =============================================================================================== # Shortcut names - #================================================================================================ - DAMASKROOT: "$TESTROOT/GitLabCI_Pipeline_$CI_PIPELINE_ID/DAMASK" + # =============================================================================================== + DAMASKROOT: "$TESTROOT/GitLabCI_Pipeline_$CI_PIPELINE_ID/DAMASK" - #================================================================================================ + # =============================================================================================== # Names of module files to load - #================================================================================================ + # =============================================================================================== # ++++++++++++ Compiler ++++++++++++++++++++++++++++++++++++++++++++++ IntelCompiler16_0: "Compiler/Intel/16.0 Libraries/IMKL/2016" IntelCompiler17_0: "Compiler/Intel/17.0 Libraries/IMKL/2017" @@ -81,11 +86,13 @@ variables: ################################################################################################### -checkout: +checkout: stage: prepareAll - before_script: + before_script: - echo $CI_PIPELINE_ID >> $TESTROOT/GitLabCI.queue - - while [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue) -ne 1 ];do sleep 5m; done + - while [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue) != 1 ] + do sleep 5m + done script: - mkdir -p $DAMASKROOT - cd $DAMASKROOT @@ -186,7 +193,7 @@ Post_ParaviewRelated: ################################################################################################### Compile_Intel: - stage: compileSpectralIntel + stage: compileSpectralIntel script: - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel - SpectralAll_compile/test.py @@ -214,7 +221,7 @@ Compile_Intel_Prepare: except: - master - release - + ################################################################################################### Spectral_PackedGeometry: stage: spectral @@ -442,9 +449,9 @@ SpectralRuntime: except: - master - release - + ################################################################################################### -AbaqusExp: +AbaqusExp: stage: createDocumentation script: - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel $Doxygen @@ -453,7 +460,7 @@ AbaqusExp: - master - release -AbaqusStd: +AbaqusStd: stage: createDocumentation script: - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel $Doxygen @@ -462,7 +469,7 @@ AbaqusStd: - master - release -Marc: +Marc: stage: createDocumentation script: - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel $Doxygen @@ -471,7 +478,7 @@ Marc: - master - release -Spectral: +Spectral: stage: createDocumentation script: - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel $Doxygen @@ -481,7 +488,7 @@ Spectral: - release ################################################################################################## -backupData: +backupData: stage: saveDocumentation script: - cd $TESTROOT/performance # location of new runtime results @@ -496,7 +503,7 @@ backupData: - development ################################################################################################## -mergeIntoMaster: +mergeIntoMaster: stage: updateMaster script: - cd $DAMASKROOT @@ -527,9 +534,9 @@ removeData: - release ################################################################################################### -removeLock: +removeLock: stage: releaseLock - before_script: + before_script: - echo 'Do nothing' when: always script: sed -i "/$CI_PIPELINE_ID/d" $TESTROOT/GitLabCI.queue From 965bee5889889f28471e4166a63f2897d7185024 Mon Sep 17 00:00:00 2001 From: Franz Roters Date: Wed, 17 May 2017 15:00:16 +0200 Subject: [PATCH 06/62] corrected thermal expansion coefficient --- examples/ConfigFiles/Kinematics_Thermal_Expansion.config | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/ConfigFiles/Kinematics_Thermal_Expansion.config b/examples/ConfigFiles/Kinematics_Thermal_Expansion.config index 6c778136f..0cedff049 100644 --- a/examples/ConfigFiles/Kinematics_Thermal_Expansion.config +++ b/examples/ConfigFiles/Kinematics_Thermal_Expansion.config @@ -1,2 +1,2 @@ (kinematics) thermal_expansion -thermal_expansion11 0.00231 +thermal_expansion11 0.0000231 From a774fffc2476fa1738890853ef021a5b8cfc0ad7 Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 17 May 2017 16:31:04 +0200 Subject: [PATCH 07/62] [skip ci] updated version information after successful test of v2.0.1-737-geff2b4c --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 5dabd576c..f37d06efd 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-735-gf52e62e +v2.0.1-737-geff2b4c From ea87e8e53cbe661f3258f3204dba7f1edf29848f Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 17 May 2017 21:33:40 +0200 Subject: [PATCH 08/62] [skip ci] updated version information after successful test of v2.0.1-740-gae28ca9 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 5dabd576c..229dd031b 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-735-gf52e62e +v2.0.1-740-gae28ca9 From 13def4ade7350f149d6a377dd65a84a91cc2a8e8 Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 18 May 2017 02:37:58 +0200 Subject: [PATCH 09/62] [skip ci] updated version information after successful test of v2.0.1-741-g965bee5 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 5dabd576c..406fd3675 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-735-gf52e62e +v2.0.1-741-g965bee5 From d5af1a7018f4fd64f9407c54475aea95c016b03d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 May 2017 05:32:08 +0200 Subject: [PATCH 10/62] superfluous brackets --- env/DAMASK.zsh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/env/DAMASK.zsh b/env/DAMASK.zsh index 23b057d68..e4a371db4 100644 --- a/env/DAMASK.zsh +++ b/env/DAMASK.zsh @@ -19,7 +19,7 @@ source $DAMASK_ROOT/CONFIG unset -f set # add DAMASK_BIN if present -[ ( "x$DAMASK_BIN" != "x" ) ] && PATH=$DAMASK_BIN:$PATH +[ "x$DAMASK_BIN" != "x" ] && PATH=$DAMASK_BIN:$PATH SOLVER=`which DAMASK_spectral || True 2>/dev/null` PROCESSING=`which postResults || True 2>/dev/null` From 92fefea3cc7dda2d5cacd904f533d6a2a484884b Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 18 May 2017 10:34:24 +0200 Subject: [PATCH 11/62] [skip ci] updated version information after successful test of v2.0.1-750-gd5af1a7 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index f37d06efd..f39a8c9e8 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-737-geff2b4c +v2.0.1-750-gd5af1a7 From 295bcd20f0ecf1fe6da94ec27149dec552312576 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 26 Apr 2017 19:18:47 +0200 Subject: [PATCH 12/62] variable cut back factor for Lp as suggested by Duancheng --- src/crystallite.f90 | 55 ++++++++++++++++++++++----------------------- src/numerics.f90 | 5 +++++ 2 files changed, 32 insertions(+), 28 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 8d2ae5f23..33b717508 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1316,7 +1316,7 @@ subroutine crystallite_integrateStateRK4() ! first Runge-Kutta step !$OMP PARALLEL !$OMP DO - do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) & call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, & @@ -1326,7 +1326,7 @@ subroutine crystallite_integrateStateRK4() !$OMP ENDDO !$OMP DO PRIVATE(p,c,NaN) - do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then c = phasememberAt(g,i,e) @@ -1335,13 +1335,13 @@ subroutine crystallite_integrateStateRK4() do mySource = 1_pInt, phase_Nsources(p) NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c))) enddo - if (NaN) then ! NaN occured in any dotState - if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... + if (NaN) then ! NaN occured in any dotState + if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) - crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped + crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped !$OMP END CRITICAL (checkTodo) - else ! if broken local... - crystallite_todo(g,i,e) = .false. ! ... skip this one next time + else ! if broken local... + crystallite_todo(g,i,e) = .false. ! ... skip this one next time endif endif endif @@ -1356,7 +1356,7 @@ subroutine crystallite_integrateStateRK4() !$OMP PARALLEL !$OMP DO PRIVATE(p,c) - do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then p = phaseAt(g,i,e) c = phasememberAt(g,i,e) @@ -1371,7 +1371,7 @@ subroutine crystallite_integrateStateRK4() !$OMP ENDDO !$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,c) - do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then p = phaseAt(g,i,e) @@ -2732,7 +2732,7 @@ subroutine crystallite_integrateStateFPI() !$OMP ENDDO !$OMP DO PRIVATE(p,c,NaN) - do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then p = phaseAt(g,i,e) @@ -2741,13 +2741,13 @@ subroutine crystallite_integrateStateFPI() do mySource = 1_pInt, phase_Nsources(p) NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c))) enddo - if (NaN) then ! NaN occured in any dotState - if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a non-local... + if (NaN) then ! NaN occured in any dotState + if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a non-local... !$OMP CRITICAL (checkTodo) - crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals done (and broken) + crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals done (and broken) !$OMP END CRITICAL (checkTodo) - else ! broken one was local... - crystallite_todo(g,i,e) = .false. ! ... done (and broken) + else ! broken one was local... + crystallite_todo(g,i,e) = .false. ! ... done (and broken) endif endif endif @@ -3107,7 +3107,7 @@ logical function crystallite_stateJump(ipc,ip,el) plasticState(p)%deltaState(1:mySizePlasticDeltaState,c) do mySource = 1_pInt, phase_Nsources(p) mySizeSourceDeltaState = sourceState(p)%p(mySource)%sizeDeltaState - if( any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(:,c)))) then ! NaN occured in deltaState + if( any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(:,c)))) then ! NaN occured in deltaState crystallite_stateJump = .false. return endif @@ -3374,8 +3374,7 @@ logical function crystallite_integrateStress(& NiterationStressLp = 0_pInt jacoCounterLp = 0_pInt - steplengthLp0 = 1.0_pReal - steplengthLp = steplengthLp0 + steplengthLp = 1.0_pReal residuumLp_old = 0.0_pReal Lpguess_old = Lpguess @@ -3428,8 +3427,8 @@ logical function crystallite_integrateStress(& !* update current residuum and check for convergence of loop - aTolLp = max(rTol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error - aTol_crystalliteStress) ! minimum lower cutoff + aTolLp = max(rTol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error + aTol_crystalliteStress) ! minimum lower cutoff residuumLp = Lpguess - Lp_constitutive if (any(IEEE_is_NaN(residuumLp))) then ! NaN in residuum... @@ -3440,16 +3439,16 @@ logical function crystallite_integrateStress(& ' ; iteration ', NiterationStressLp,& ' >> returning..!' #endif - return ! ...me = .false. to inform integrator about problem - elseif (norm2(residuumLp) < aTolLp) then ! converged if below absolute tolerance - exit LpLoop ! ...leave iteration loop + return ! ...me = .false. to inform integrator about problem + elseif (norm2(residuumLp) < aTolLp) then ! converged if below absolute tolerance + exit LpLoop ! ...leave iteration loop elseif ( NiterationStressLp == 1_pInt & - .or. norm2(residuumLp) < norm2(residuumLp_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)... - residuumLp_old = residuumLp ! ...remember old values and... + .or. norm2(residuumLp) < norm2(residuumLp_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)... + residuumLp_old = residuumLp ! ...remember old values and... Lpguess_old = Lpguess - steplengthLp = steplengthLp0 ! ...proceed with normal step length (calculate new search direction) - else ! not converged and residuum not improved... - steplengthLp = 0.5_pReal * steplengthLp ! ...try with smaller step length in same direction + steplengthLp = 1.0_pReal ! ...proceed with normal step length (calculate new search direction) + else ! not converged and residuum not improved... + steplengthLp = subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction Lpguess = Lpguess_old + steplengthLp * deltaLp cycle LpLoop endif diff --git a/src/numerics.f90 b/src/numerics.f90 index eb974b3c4..108ec22d1 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -42,6 +42,7 @@ module numerics subStepMinHomog = 1.0e-3_pReal, & !< minimum (relative) size of sub-step allowed during cutback in homogenization subStepSizeCryst = 0.25_pReal, & !< size of first substep when cutback in crystallite subStepSizeHomog = 0.25_pReal, & !< size of first substep when cutback in homogenization + subStepSizeLp = 0.5_pReal, & !< size of first substep when cutback in Lp calculation stepIncreaseCryst = 1.5_pReal, & !< increase of next substep size when previous substep converged in crystallite stepIncreaseHomog = 1.5_pReal, & !< increase of next substep size when previous substep converged in homogenization rTol_crystalliteState = 1.0e-6_pReal, & !< relative tolerance in crystallite state loop @@ -295,6 +296,8 @@ subroutine numerics_init subStepSizeCryst = IO_floatValue(line,chunkPos,2_pInt) case ('stepincreasecryst') stepIncreaseCryst = IO_floatValue(line,chunkPos,2_pInt) + case ('substepsizelp') + subStepSizeLp = IO_floatValue(line,chunkPos,2_pInt) case ('substepminhomog') subStepMinHomog = IO_floatValue(line,chunkPos,2_pInt) case ('substepsizehomog') @@ -515,6 +518,7 @@ subroutine numerics_init write(6,'(a24,1x,es8.1)') ' subStepMinCryst: ',subStepMinCryst write(6,'(a24,1x,es8.1)') ' subStepSizeCryst: ',subStepSizeCryst write(6,'(a24,1x,es8.1)') ' stepIncreaseCryst: ',stepIncreaseCryst + write(6,'(a24,1x,es8.1)') ' subStepSizeLp: ',subStepSizeLp write(6,'(a24,1x,i8)') ' nState: ',nState write(6,'(a24,1x,i8)') ' nStress: ',nStress write(6,'(a24,1x,es8.1)') ' rTol_crystalliteState: ',rTol_crystalliteState @@ -643,6 +647,7 @@ subroutine numerics_init if (subStepMinCryst <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepMinCryst') if (subStepSizeCryst <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeCryst') if (stepIncreaseCryst <= 0.0_pReal) call IO_error(301_pInt,ext_msg='stepIncreaseCryst') + if (subStepSizeLp <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeLp') if (subStepMinHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepMinHomog') if (subStepSizeHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeHomog') if (stepIncreaseHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='stepIncreaseHomog') From 590a5c8b91192474fefc4bd8251fdf7b132d3985 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 28 Apr 2017 14:31:03 +0200 Subject: [PATCH 13/62] using also variable cutback factor for Li --- src/crystallite.f90 | 10 +++++----- src/numerics.f90 | 5 +++++ 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 33b717508..8601697ba 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -3179,7 +3179,9 @@ logical function crystallite_integrateStress(& aTol_crystalliteStress, & rTol_crystalliteStress, & iJacoLpresiduum, & - numerics_integrationMode + numerics_integrationMode, & + subStepSizeLp, & + subStepSizeLi use debug, only: debug_level, & debug_crystallite, & debug_levelBasic, & @@ -3265,9 +3267,7 @@ logical function crystallite_integrateStress(& dLp_dT3333, & dLi_dT3333 real(pReal) detInvFi, & ! determinant of InvFi - steplengthLp0, & steplengthLp, & - steplengthLi0, & steplengthLi, & dt, & ! time increment aTolLp, & @@ -3529,9 +3529,9 @@ logical function crystallite_integrateStress(& .or. norm2(residuumLi) < norm2(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) + steplengthLi = 1.0_pReal ! ...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 = subStepSizeLi * steplengthLi ! ...try with smaller step length in same direction Liguess = Liguess_old + steplengthLi * deltaLi cycle LiLoop endif diff --git a/src/numerics.f90 b/src/numerics.f90 index 108ec22d1..db7bf0fe4 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -43,6 +43,7 @@ module numerics subStepSizeCryst = 0.25_pReal, & !< size of first substep when cutback in crystallite subStepSizeHomog = 0.25_pReal, & !< size of first substep when cutback in homogenization subStepSizeLp = 0.5_pReal, & !< size of first substep when cutback in Lp calculation + subStepSizeLi = 0.5_pReal, & !< size of first substep when cutback in Li calculation stepIncreaseCryst = 1.5_pReal, & !< increase of next substep size when previous substep converged in crystallite stepIncreaseHomog = 1.5_pReal, & !< increase of next substep size when previous substep converged in homogenization rTol_crystalliteState = 1.0e-6_pReal, & !< relative tolerance in crystallite state loop @@ -298,6 +299,8 @@ subroutine numerics_init stepIncreaseCryst = IO_floatValue(line,chunkPos,2_pInt) case ('substepsizelp') subStepSizeLp = IO_floatValue(line,chunkPos,2_pInt) + case ('substepsizeli') + subStepSizeLi = IO_floatValue(line,chunkPos,2_pInt) case ('substepminhomog') subStepMinHomog = IO_floatValue(line,chunkPos,2_pInt) case ('substepsizehomog') @@ -519,6 +522,7 @@ subroutine numerics_init write(6,'(a24,1x,es8.1)') ' subStepSizeCryst: ',subStepSizeCryst write(6,'(a24,1x,es8.1)') ' stepIncreaseCryst: ',stepIncreaseCryst write(6,'(a24,1x,es8.1)') ' subStepSizeLp: ',subStepSizeLp + write(6,'(a24,1x,es8.1)') ' subStepSizeLi: ',subStepSizeLi write(6,'(a24,1x,i8)') ' nState: ',nState write(6,'(a24,1x,i8)') ' nStress: ',nStress write(6,'(a24,1x,es8.1)') ' rTol_crystalliteState: ',rTol_crystalliteState @@ -648,6 +652,7 @@ subroutine numerics_init if (subStepSizeCryst <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeCryst') if (stepIncreaseCryst <= 0.0_pReal) call IO_error(301_pInt,ext_msg='stepIncreaseCryst') if (subStepSizeLp <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeLp') + if (subStepSizeLi <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeLi') if (subStepMinHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepMinHomog') if (subStepSizeHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeHomog') if (stepIncreaseHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='stepIncreaseHomog') From 9be3cac9472662face5e841e5b0a722b59487610 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 15 May 2017 13:06:15 +0200 Subject: [PATCH 14/62] unused variable --- src/homogenization.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 34bc93af0..b70eacc56 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -928,7 +928,6 @@ subroutine materialpoint_postResults constitutive_plasticity_maxSizePostResults, & constitutive_source_maxSizePostResults, & #endif - constitutive_postResults use crystallite, only: & #ifdef FEM crystallite_maxSizePostResults, & From 12f66fd806fc6a97617b9a73c1d1e2feec4320b2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 May 2017 11:33:11 +0200 Subject: [PATCH 15/62] bug introduced during merge --- src/crystallite.f90 | 3 +-- src/homogenization.f90 | 4 ++-- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 8601697ba..c5bd4d979 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -3353,8 +3353,7 @@ logical function crystallite_integrateStress(& NiterationStressLi = 0_pInt jacoCounterLi = 0_pInt - steplengthLi0 = 1.0_pReal - steplengthLi = steplengthLi0 + steplengthLi = 1.0_pReal residuumLi_old = 0.0_pReal LiLoop: do diff --git a/src/homogenization.f90 b/src/homogenization.f90 index b70eacc56..c8c5fad01 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -923,10 +923,10 @@ subroutine materialpoint_postResults material_phase, & homogenization_Ngrains, & microstructure_crystallite - use constitutive, only: & #ifdef FEM + use constitutive, only: & constitutive_plasticity_maxSizePostResults, & - constitutive_source_maxSizePostResults, & + constitutive_source_maxSizePostResults #endif use crystallite, only: & #ifdef FEM From 1bfdb4d6f35b920d68dbaaf33d83e6d4dd50a2f3 Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 18 May 2017 16:38:34 +0200 Subject: [PATCH 16/62] [skip ci] updated version information after successful test of v2.0.1-755-g12f66fd --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index f39a8c9e8..a5f8d7ad5 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-750-gd5af1a7 +v2.0.1-755-g12f66fd From 5ec4c8c8b6af1b682287bbdc2acc7eec22b78400 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 May 2017 17:07:19 +0200 Subject: [PATCH 17/62] line delimiters missing, still valid yaml but does not run --- .gitlab-ci.yml | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index cc541307f..0542a5252 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -24,11 +24,11 @@ stages: ################################################################################################### before_script: - - if [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue)x == 'x' ] - then echo $CI_PIPELINE_ID >> $TESTROOT/GitLabCI.queue + - if [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue)x == 'x' ]; + then echo $CI_PIPELINE_ID >> $TESTROOT/GitLabCI.queue; fi - - while [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue) != 1 ] - do sleep 5m + - while [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue) != 1 ]; + do sleep 5m; done - source $DAMASKROOT/DAMASK_env.sh - cd $DAMASKROOT/PRIVATE/testing @@ -90,8 +90,8 @@ checkout: stage: prepareAll before_script: - echo $CI_PIPELINE_ID >> $TESTROOT/GitLabCI.queue - - while [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue) != 1 ] - do sleep 5m + - while [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue) != 1 ]; + do sleep 5m; done script: - mkdir -p $DAMASKROOT From 867a7f2ee2773ad510e99ab4b608e8d0275eff51 Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 18 May 2017 22:09:38 +0200 Subject: [PATCH 18/62] [skip ci] updated version information after successful test of v2.0.1-760-g5ec4c8c --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index a5f8d7ad5..909a4dad5 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-755-g12f66fd +v2.0.1-760-g5ec4c8c From 691f338f5defa6a8c6d22af2e1be87e4f0ce07f4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 19 May 2017 09:33:58 +0200 Subject: [PATCH 19/62] executable bit got lost --- processing/pre/mentat_pbcOnBoxMesh.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100644 => 100755 processing/pre/mentat_pbcOnBoxMesh.py diff --git a/processing/pre/mentat_pbcOnBoxMesh.py b/processing/pre/mentat_pbcOnBoxMesh.py old mode 100644 new mode 100755 From b1ff3b533158579c2c35235fa266457fea621e9b Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 19 May 2017 14:30:49 +0200 Subject: [PATCH 20/62] [skip ci] updated version information after successful test of v2.0.1-762-g691f338 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 909a4dad5..0a6808d61 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-760-g5ec4c8c +v2.0.1-762-g691f338 From 8e86a35492f762e0b93ce7a5edf3e256cde6c4c7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 20 May 2017 15:36:56 +0200 Subject: [PATCH 21/62] location of env files has changed --- .gitlab-ci.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 0542a5252..4455c9cf5 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -30,7 +30,7 @@ before_script: - while [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue) != 1 ]; do sleep 5m; done - - source $DAMASKROOT/DAMASK_env.sh + - source $DAMASKROOT/env/DAMASK.sh - cd $DAMASKROOT/PRIVATE/testing ################################################################################################### @@ -99,7 +99,7 @@ checkout: - git clone -q git@magit1.mpie.de:damask/DAMASK.git . - git checkout $CI_COMMIT_SHA - git submodule update --init - - source DAMASK_env.sh + - source env/DAMASK.sh - make processing except: - master From 00a3c9751a018a1b8680bc2b1f9edb108b78d88a Mon Sep 17 00:00:00 2001 From: Test User Date: Sat, 20 May 2017 20:32:53 +0200 Subject: [PATCH 22/62] [skip ci] updated version information after successful test of v2.0.1-764-g8e86a35 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 0a6808d61..d4ba8a1cb 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-762-g691f338 +v2.0.1-764-g8e86a35 From d3467705cad7bab00edb84363532df5389c915b5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 22 May 2017 10:08:16 +0200 Subject: [PATCH 23/62] compilation exception not needed any more for intrinsic NaN function --- src/CMakeLists.txt | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 83647d8d0..a3625ef1a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,16 +1,5 @@ # special flags for some files if (${CMAKE_Fortran_COMPILER_ID} STREQUAL "GNU") - - SET_SOURCE_FILES_PROPERTIES( "prec.f90" PROPERTIES - COMPILE_FLAGS "-fno-range-check -fall-intrinsics -fno-fast-math") - # fno-range-check: Disable range checking on results of simplification of constant expressions during compilation - # --> allows the definition of DAMASK_NaN - #-fall-intrinsics: all intrinsic procedures (including the GNU-specific extensions) are accepted. -Wintrinsics-std will be ignored - # and no user-defined procedure with the same name as any intrinsic will be called except when it is explicitly declared external - # --> allows the use of 'isnan' - #-fno-fast-math: - # --> otherwise, when setting -ffast-math, isnan always evaluates to false (I would call it a bug) - SET_SOURCE_FILES_PROPERTIES( "lattice.f90" PROPERTIES COMPILE_FLAGS "-ffree-line-length-240") # long lines for interaction matrix From 4c0b2e85980e74251978d3e4e276f4647b28ccff Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 22 May 2017 15:01:29 +0200 Subject: [PATCH 24/62] [skip ci] updated version information after successful test of v2.0.1-766-gd346770 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index d4ba8a1cb..7cc702b88 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-764-g8e86a35 +v2.0.1-766-gd346770 From 2f3c0ec0c6860caff6acc3ed220d2b03a9a15782 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 24 May 2017 06:31:58 +0200 Subject: [PATCH 25/62] consistent naming --- env/DAMASK.sh | 2 +- env/DAMASK.zsh | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/env/DAMASK.sh b/env/DAMASK.sh index 2674505f7..3e9889fec 100644 --- a/env/DAMASK.sh +++ b/env/DAMASK.sh @@ -16,7 +16,7 @@ if [ ${BASH_SOURCE##*/} == "DAMASK.sh" ]; then fi # shorthand command to change to DAMASK_ROOT directory -eval "function damask() { cd $DAMASK_ROOT; }" +eval "function DAMASK_root() { cd $DAMASK_ROOT; }" # defining set() allows to source the same file for tcsh and bash, with and without space around = set() { diff --git a/env/DAMASK.zsh b/env/DAMASK.zsh index e4a371db4..c662e6aee 100644 --- a/env/DAMASK.zsh +++ b/env/DAMASK.zsh @@ -9,7 +9,7 @@ else fi # shorthand command to change to DAMASK_ROOT directory -eval "function damask() { cd $DAMASK_ROOT; }" +eval "function DAMASK_root() { cd $DAMASK_ROOT; }" # defining set() allows to source the same file for tcsh and zsh, with and without space around = set() { From b9e986783c3f26db98f4f3cdcc8dde5b290c561d Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 24 May 2017 11:30:47 +0200 Subject: [PATCH 26/62] [skip ci] updated version information after successful test of v2.0.1-768-g2f3c0ec --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 7cc702b88..2ad87936f 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-766-gd346770 +v2.0.1-768-g2f3c0ec From 8b529d8b0431918696fa120ce016734e2ee00f33 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 24 May 2017 18:12:36 +0200 Subject: [PATCH 27/62] cleaner finalizing in case of interrupted simulation --- src/DAMASK_spectral.f90 | 133 ++++++++++++++++++++-------------------- 1 file changed, 68 insertions(+), 65 deletions(-) mode change 100644 => 100755 src/DAMASK_spectral.f90 diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90 old mode 100644 new mode 100755 index 6f1794e35..6a932fd82 --- a/src/DAMASK_spectral.f90 +++ b/src/DAMASK_spectral.f90 @@ -364,12 +364,12 @@ program DAMASK_spectral case (DAMASK_spectral_SolverBasicPETSc_label) call basicPETSc_init case (DAMASK_spectral_SolverAL_label) - if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0 .and. worldrank == 0_pInt) & + if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) & call IO_warning(42_pInt, ext_msg='debug Divergence') call AL_init case (DAMASK_spectral_SolverPolarisation_label) - if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0 .and. worldrank == 0_pInt) & + if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) & call IO_warning(42_pInt, ext_msg='debug Divergence') call Polarisation_init @@ -450,8 +450,7 @@ program DAMASK_spectral if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_write') enddo fileOffset = fileOffset + sum(outputSize) ! forward to current file position - if (worldrank == 0) & - write(6,'(1/,a)') ' ... writing initial configuration to file ........................' + write(6,'(1/,a)') ' ... writing initial configuration to file ........................' endif !-------------------------------------------------------------------------------------------------- ! loopping over loadcases @@ -498,21 +497,20 @@ program DAMASK_spectral !-------------------------------------------------------------------------------------------------- ! report begin of new increment - if (worldrank == 0) then - write(6,'(/,a)') ' ###########################################################################' - write(6,'(1x,a,es12.5'//& - ',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//& - ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//& - ',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') & - 'Time', time, & - 's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,& - '-', stepFraction, '/', subStepFactor**cutBackLevel,& - ' of load case ', currentLoadCase,'/',size(loadCases) - flush(6) - write(incInfo,'(a,'//IO_intOut(totalIncsCounter)//',a,'//IO_intOut(sum(loadCases%incs))//& - ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') & - 'Increment ',totalIncsCounter,'/',sum(loadCases%incs),& - '-',stepFraction, '/', subStepFactor**cutBackLevel + write(6,'(/,a)') ' ###########################################################################' + write(6,'(1x,a,es12.5'//& + ',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//& + ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//& + ',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') & + 'Time', time, & + 's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,& + '-', stepFraction, '/', subStepFactor**cutBackLevel,& + ' of load case ', currentLoadCase,'/',size(loadCases) + flush(6) + write(incInfo,'(a,'//IO_intOut(totalIncsCounter)//',a,'//IO_intOut(sum(loadCases%incs))//& + ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') & + 'Increment ',totalIncsCounter,'/',sum(loadCases%incs),& + '-',stepFraction, '/', subStepFactor**cutBackLevel endif !-------------------------------------------------------------------------------------------------- @@ -595,7 +593,7 @@ program DAMASK_spectral cutBack = .False. if(solres(1)%termIll .or. .not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found if (cutBackLevel < maxCutBack) then ! do cut back - if (worldrank == 0) write(6,'(/,a)') ' cut back detected' + write(6,'(/,a)') ' cut back detected' cutBack = .True. stepFraction = (stepFraction - 1_pInt) * subStepFactor ! adjust to new denominator cutBackLevel = cutBackLevel + 1_pInt @@ -603,11 +601,15 @@ program DAMASK_spectral timeinc = timeinc/2.0_pReal elseif (solres(1)%termIll) then ! material point model cannot find a solution, exit in any casy call IO_warning(850_pInt) + call MPI_file_close(resUnit,ierr) + close(statUnit) call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written elseif (continueCalculation == 1_pInt) then guess = .true. ! accept non converged BVP solution else ! default behavior, exit if spectral solver does not converge call IO_warning(850_pInt) + call MPI_file_close(resUnit,ierr) + close(statUnit) call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written endif else @@ -624,13 +626,11 @@ program DAMASK_spectral cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half number of subincs next inc if(all(solres(:)%converged)) then ! report converged inc convergedCounter = convergedCounter + 1_pInt - if (worldrank == 0) & - write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & - ' increment ', totalIncsCounter, ' converged' + write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & + ' increment ', totalIncsCounter, ' converged' else - if (worldrank == 0) & - write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc - ' increment ', totalIncsCounter, ' NOT converged' + write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc + ' increment ', totalIncsCounter, ' NOT converged' notConvergedCounter = notConvergedCounter + 1_pInt endif; flush(6) if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency @@ -665,45 +665,15 @@ program DAMASK_spectral !-------------------------------------------------------------------------------------------------- ! report summary of whole calculation - if (worldrank == 0) then - write(6,'(/,a)') ' ###########################################################################' - write(6,'(1x,i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', & - notConvergedCounter + convergedCounter, ' (', & - real(convergedCounter, pReal)/& - real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, & - ' %) increments converged!' - endif + write(6,'(/,a)') ' ###########################################################################' + write(6,'(1x,i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', & + notConvergedCounter + convergedCounter, ' (', & + real(convergedCounter, pReal)/& + real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, & + ' %) increments converged!' call MPI_file_close(resUnit,ierr) close(statUnit) - do field = 1, nActiveFields - select case(loadCases(1)%ID(field)) - case(FIELD_MECH_ID) - select case (spectral_solver) - case (DAMASK_spectral_SolverBasicPETSc_label) - call BasicPETSC_destroy() - case (DAMASK_spectral_SolverAL_label) - call AL_destroy() - case (DAMASK_spectral_SolverPolarisation_label) - call Polarisation_destroy() - end select - case(FIELD_THERMAL_ID) - call spectral_thermal_destroy() - case(FIELD_DAMAGE_ID) - call spectral_damage_destroy() - end select - enddo - call utilities_destroy() - - call PETScFinalize(ierr); CHKERRQ(ierr) - -#ifdef _OPENMP - call MPI_finalize(i) - if (i /= 0_pInt) then - call IO_error(error_ID=894, el=i, ext_msg="Finalize()") - endif -#endif - if (notConvergedCounter > 0_pInt) call quit(3_pInt) ! error if some are not converged call quit(0_pInt) ! no complains ;) @@ -721,11 +691,43 @@ end program DAMASK_spectral subroutine quit(stop_id) use prec, only: & pInt + use spectral_mech_Basic, only: & + BasicPETSC_destroy + use spectral_mech_AL, only: & + AL_destroy + use spectral_mech_Polarisation, only: & + Polarisation_destroy + use spectral_damage, only: & + spectral_damage_destroy + use spectral_thermal, only: & + spectral_thermal_destroy + use spectral_utilities, only: & + utilities_destroy implicit none + + #include integer(pInt), intent(in) :: stop_id integer, dimension(8) :: dateAndTime ! type default integer + integer(pInt) :: error = 0_pInt + PetscErrorCode :: ierr = 0 + logical :: ErrorInQuit + + call BasicPETSC_destroy() + call AL_destroy() + call Polarisation_destroy() + call spectral_damage_destroy() + call spectral_thermal_destroy() + call utilities_destroy() + call PETScFinalize(ierr) + if(ierr /= 0) write(6,'(a)') ' Error in PETScFinalize' +#ifdef _OPENMP + call MPI_finalize(error) + if(error /= 0) write(6,'(a)') ' Error in MPI_finalize' +#endif + ErrorInQuit = (ierr /= 0 .or. error /= 0_pInt) + call date_and_time(values = dateAndTime) write(6,'(/,a)') 'DAMASK terminated on:' write(6,'(a,2(i2.2,a),i4.4)') 'Date: ',dateAndTime(3),'/',& @@ -735,12 +737,13 @@ subroutine quit(stop_id) dateAndTime(6),':',& dateAndTime(7) - if (stop_id == 0_pInt) stop 0 ! normal termination - if (stop_id < 0_pInt) then ! terminally ill, restart might help + if (stop_id == 0_pInt .and. .not. ErrorInQuit) stop 0 ! normal termination + if (stop_id < 0_pInt .and. .not. ErrorInQuit) then ! terminally ill, restart might help write(0,'(a,i6)') 'restart information available at ', stop_id*(-1_pInt) stop 2 endif - if (stop_id == 3_pInt) stop 3 ! not all incs converged + if (stop_id == 3_pInt .and. .not. ErrorInQuit) stop 3 ! not all incs converged + stop 1 ! error (message from IO_error) end subroutine quit From 618bf95a439c561798937d2154d2e0e5882411dc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 24 May 2017 22:16:35 +0200 Subject: [PATCH 28/62] did not compile, adjusted dummy compilation routine to detect such errors --- src/CMakeLists.txt | 4 +++- src/DAMASK_spectral.f90 | 7 +++++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a3625ef1a..eb1448028 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -174,8 +174,10 @@ if ("${PROJECT_NAME}" STREQUAL "DAMASK_spectral") list(APPEND OBJECTFILES $) if(NOT "${CMAKE_BUILD_TYPE}" STREQUAL "SYNTAXONLY") add_executable(DAMASK_spectral "DAMASK_spectral.f90" ${OBJECTFILES}) - add_dependencies(DAMASK_spectral SPECTRAL_SOLVER) + else() + add_library(DAMASK_spectral OBJECT "DAMASK_spectral.f90") endif() + add_dependencies(DAMASK_spectral SPECTRAL_SOLVER) elseif ("${PROJECT_NAME}" STREQUAL "DAMASK_FEM") add_library(FEM_UTILITIES OBJECT "FEM_utilities.f90") add_dependencies(FEM_UTILITIES DAMASK_CPFE) diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90 index 6a932fd82..dfa1746b2 100755 --- a/src/DAMASK_spectral.f90 +++ b/src/DAMASK_spectral.f90 @@ -511,7 +511,6 @@ program DAMASK_spectral ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') & 'Increment ',totalIncsCounter,'/',sum(loadCases%incs),& '-',stepFraction, '/', subStepFactor**cutBackLevel - endif !-------------------------------------------------------------------------------------------------- ! forward fields @@ -706,13 +705,17 @@ subroutine quit(stop_id) implicit none - #include +#include integer(pInt), intent(in) :: stop_id integer, dimension(8) :: dateAndTime ! type default integer integer(pInt) :: error = 0_pInt PetscErrorCode :: ierr = 0 logical :: ErrorInQuit + external :: & + PETScFinalize, & + MPI_finalize + call BasicPETSC_destroy() call AL_destroy() call Polarisation_destroy() From 83f8795638877ee437d86cbdb886c18cf0bf7f26 Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 25 May 2017 03:10:53 +0200 Subject: [PATCH 29/62] [skip ci] updated version information after successful test of v2.0.1-771-g618bf95 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 2ad87936f..e4cefa5d6 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-768-g2f3c0ec +v2.0.1-771-g618bf95 From e5af0630fe46f81c9f20adb6471cb39aec270eb7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 29 May 2017 10:32:22 +0200 Subject: [PATCH 30/62] gaussian filter (wrapper no ndimage) --- PRIVATE | 2 +- processing/post/addGaussian.py | 137 +++++++++++++++++++++++++++++++ processing/pre/geom_fromTable.py | 3 +- 3 files changed, 140 insertions(+), 2 deletions(-) create mode 100755 processing/post/addGaussian.py diff --git a/PRIVATE b/PRIVATE index 19a53f622..596ec41a5 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 19a53f6229603aeafb2466b58679a1cd04fc0142 +Subproject commit 596ec41a5117c5f2a09356ffeee8cd8ce9a149d7 diff --git a/processing/post/addGaussian.py b/processing/post/addGaussian.py new file mode 100755 index 000000000..91d3e4667 --- /dev/null +++ b/processing/post/addGaussian.py @@ -0,0 +1,137 @@ +#!/usr/bin/env python2.7 +# -*- coding: UTF-8 no BOM -*- + +import os,sys +import numpy as np +from optparse import OptionParser +from scipy import ndimage +import damask + +scriptName = os.path.splitext(os.path.basename(__file__))[0] +scriptID = ' '.join([scriptName,damask.version]) + + +# -------------------------------------------------------------------- +# MAIN +# -------------------------------------------------------------------- + +parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [ASCIItable(s)]', description = """ +Add column(s) containing gradient of requested column(s). +Operates on periodic ordered three-dimensional data sets. +Deals with both vector- and scalar fields. + +""", version = scriptID) + +parser.add_option('-p','--pos','--periodiccellcenter', + dest = 'pos', + type = 'string', metavar = 'string', + help = 'label of coordinates [%default]') +parser.add_option('-s','--scalar', + dest = 'scalar', + action = 'extend', metavar = '', + help = 'label(s) of scalar field values') +parser.add_option('-o','--order', + dest = 'order', + type = int, + metavar = 'int', + help = 'order of the filter') +parser.add_option('--sigma', + dest = 'sigma', + type = float, + metavar = 'float', + help = 'standard deviation') +parser.add_option('--periodic', + dest = 'periodic', + action = 'store_true', + help = 'assume periodic grain structure' + ) + + + +parser.set_defaults(pos = 'pos', + order = 0, + sigma = 1, + periodic = False + ) + +(options,filenames) = parser.parse_args() + +if options.scalar is None: + parser.error('no data column specified.') + +# --- loop over input files ------------------------------------------------------------------------ + +if filenames == []: filenames = [None] + +for name in filenames: + try: table = damask.ASCIItable(name = name,buffered = False) + except: continue + damask.util.report(scriptName,name) + +# ------------------------------------------ read header ------------------------------------------ + + table.head_read() + +# ------------------------------------------ sanity checks ---------------------------------------- + + items = { + 'scalar': {'dim': 1, 'shape': [1], 'labels':options.scalar, 'active':[], 'column': []}, + } + errors = [] + remarks = [] + column = {} + + if table.label_dimension(options.pos) != 3: errors.append('coordinates {} are not a vector.'.format(options.pos)) + else: colCoord = table.label_index(options.pos) + + for type, data in items.iteritems(): + for what in (data['labels'] if data['labels'] is not None else []): + dim = table.label_dimension(what) + if dim != data['dim']: remarks.append('column {} is not a {}.'.format(what,type)) + else: + items[type]['active'].append(what) + items[type]['column'].append(table.label_index(what)) + + if remarks != []: damask.util.croak(remarks) + if errors != []: + damask.util.croak(errors) + table.close(dismiss = True) + continue + +# ------------------------------------------ assemble header -------------------------------------- + + table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) + for type, data in items.iteritems(): + for label in data['active']: + table.labels_append(['Gauss{}({})'.format(options.sigma,label)]) # extend ASCII header with new labels + table.head_write() + +# --------------- figure out size and grid --------------------------------------------------------- + + table.data_readArray() + + coords = [np.unique(table.data[:,colCoord+i]) for i in range(3)] + mincorner = np.array(map(min,coords)) + maxcorner = np.array(map(max,coords)) + grid = np.array(map(len,coords),'i') + size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1) + size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 equal to smallest among other ones + +# ------------------------------------------ process value field ----------------------------------- + + stack = [table.data] + for type, data in items.iteritems(): + for i,label in enumerate(data['active']): + stack.append(ndimage.filters.gaussian_filter(table.data[:,data['column'][i]], + options.sigma,options.order, + mode = 'wrap' if options.periodic else 'nearest' + ).reshape([table.data.shape[0],1]) + ) + +# ------------------------------------------ output result ----------------------------------------- + if len(stack) > 1: table.data = np.hstack(tuple(stack)) + table.data_writeArray('%.12g') + +# ------------------------------------------ output finalization ----------------------------------- + + table.close() # close input ASCII table (works for stdin) diff --git a/processing/pre/geom_fromTable.py b/processing/pre/geom_fromTable.py index 9af92888d..b10bc9f88 100755 --- a/processing/pre/geom_fromTable.py +++ b/processing/pre/geom_fromTable.py @@ -23,7 +23,7 @@ Generate geometry description and material configuration from position, phase, a parser.add_option('--coordinates', dest = 'pos', type = 'string', metavar = 'string', - help = 'coordinates label') + help = 'coordinates label (%default)') parser.add_option('--phase', dest = 'phase', type = 'string', metavar = 'string', @@ -90,6 +90,7 @@ parser.set_defaults(symmetry = [damask.Symmetry.lattices[-1]], homogenization = 1, crystallite = 1, verbose = False, + pos = 'pos', ) (options,filenames) = parser.parse_args() From 5bd93e2117f5508c9f0390ebf9fb5577e39c6e27 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 29 May 2017 10:46:26 +0200 Subject: [PATCH 31/62] trouble when restarting pipeline --- .gitlab-ci.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 4455c9cf5..df8991900 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -96,6 +96,7 @@ checkout: script: - mkdir -p $DAMASKROOT - cd $DAMASKROOT + - if [ -d DAMASK ]; then rm -rf DAMASK; fi # there might be some leftovers from a failed clone - git clone -q git@magit1.mpie.de:damask/DAMASK.git . - git checkout $CI_COMMIT_SHA - git submodule update --init From d7b6c2005ee075fa979cdb3f7e443fed6fe9cf52 Mon Sep 17 00:00:00 2001 From: Test User Date: Tue, 30 May 2017 05:40:47 +0200 Subject: [PATCH 32/62] [skip ci] updated version information after successful test of v2.0.1-774-g5bd93e2 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index e4cefa5d6..64ca177c7 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-771-g618bf95 +v2.0.1-774-g5bd93e2 From 6599f7299da0172307919f6ed14bff0b0ed40ce5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 1 Jun 2017 10:00:21 +0200 Subject: [PATCH 33/62] correct description --- PRIVATE | 2 +- processing/post/addGaussian.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/PRIVATE b/PRIVATE index 596ec41a5..19a53f622 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 596ec41a5117c5f2a09356ffeee8cd8ce9a149d7 +Subproject commit 19a53f6229603aeafb2466b58679a1cd04fc0142 diff --git a/processing/post/addGaussian.py b/processing/post/addGaussian.py index 91d3e4667..bc5599d49 100755 --- a/processing/post/addGaussian.py +++ b/processing/post/addGaussian.py @@ -16,9 +16,9 @@ scriptID = ' '.join([scriptName,damask.version]) # -------------------------------------------------------------------- parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [ASCIItable(s)]', description = """ -Add column(s) containing gradient of requested column(s). -Operates on periodic ordered three-dimensional data sets. -Deals with both vector- and scalar fields. +Add column(s) containing Gaussian filtered values of requested column(s). +Operates on periodic and non-periodic ordered three-dimensional data sets. +For Details see scipy.ndimage documentation. """, version = scriptID) From 31859c08363c460542549c553a3227f2514c8dbc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 1 Jun 2017 16:08:29 +0200 Subject: [PATCH 34/62] improve runtime test --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 19a53f622..057371b82 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 19a53f6229603aeafb2466b58679a1cd04fc0142 +Subproject commit 057371b82e3f5e880271b9631ace46c54280a753 From 6ce0888864e2165248ba6bf25eabed3cd105e64b Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 1 Jun 2017 21:27:05 +0200 Subject: [PATCH 35/62] [skip ci] updated version information after successful test of v2.0.1-777-g31859c0 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 64ca177c7..1ad92c82c 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-774-g5bd93e2 +v2.0.1-777-g31859c0 From f085f61c4dac5ba7149386082334f51203b71d8c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jun 2017 09:25:58 +0200 Subject: [PATCH 36/62] testing hook that checks for executable bit --- processing/post/addGrainID.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processing/post/addGrainID.py b/processing/post/addGrainID.py index c39d67983..45034034b 100755 --- a/processing/post/addGrainID.py +++ b/processing/post/addGrainID.py @@ -148,7 +148,7 @@ for name in filenames: bg.set_message('reading positions...') - table.data_readArray(options.pos) # read position vectors + table.data_readArray(options.pos) # read position vectors grainID = -np.ones(len(table.data),dtype=int) start = tick = time.clock() From d9c8072a77dd7845ac1d0504b0718f4e1e1d0da1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jun 2017 09:54:54 +0200 Subject: [PATCH 37/62] checking file mode --- installation/symlink_Processing.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/installation/symlink_Processing.py b/installation/symlink_Processing.py index d10b5af55..50ea86e68 100755 --- a/installation/symlink_Processing.py +++ b/installation/symlink_Processing.py @@ -30,7 +30,7 @@ for subDir in processing_subDirs: for theFile in os.listdir(theDir): theName,theExt = os.path.splitext(theFile) - if theExt in processing_extensions: # only consider files with proper extensions + if theExt in processing_extensions: # only consider files with proper extensions src = os.path.abspath(os.path.join(theDir,theFile)) sym_link = os.path.abspath(os.path.join(binDir,theName)) From 931eee496df64ffd433e49dde07931456c249701 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jun 2017 09:58:44 +0200 Subject: [PATCH 38/62] trying executable check --- installation/mods_MarcMentat/apply_DAMASK_modifications.sh | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/installation/mods_MarcMentat/apply_DAMASK_modifications.sh b/installation/mods_MarcMentat/apply_DAMASK_modifications.sh index d1c1e37dd..8cd987ba2 100755 --- a/installation/mods_MarcMentat/apply_DAMASK_modifications.sh +++ b/installation/mods_MarcMentat/apply_DAMASK_modifications.sh @@ -163,5 +163,4 @@ echo '' echo 'precompiling $VERSION HYPELA2 user subroutine...' echo 'not yet implemented..!' -echo '' -echo 'done.' +echo -e '\ndone.' From 4d43b05a155dddc85ec55be070954808d1f9dbf6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jun 2017 09:59:54 +0200 Subject: [PATCH 39/62] still working on executable detection --- installation/mods_MarcMentat/apply_DAMASK_modifications.sh | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/installation/mods_MarcMentat/apply_DAMASK_modifications.sh b/installation/mods_MarcMentat/apply_DAMASK_modifications.sh index 8cd987ba2..7eef9c729 100755 --- a/installation/mods_MarcMentat/apply_DAMASK_modifications.sh +++ b/installation/mods_MarcMentat/apply_DAMASK_modifications.sh @@ -102,8 +102,7 @@ for filename in 'edit_window' \ done # Mentat scripts -echo '' -echo 'adapting Mentat menus...' +echo -e '\nadapting Mentat menus...' theDIR=$INSTALLDIR/mentat$VERSION/menus for filename in 'job_run.ms'; do cp $SCRIPTLOCATION/$VERSION/Mentat_menus/$filename $theDIR From 8ddb43b6db9f0748330b529a29adaa7a21474fef Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jun 2017 10:00:53 +0200 Subject: [PATCH 40/62] testing executable test --- installation/mods_MarcMentat/test.py | 1 + 1 file changed, 1 insertion(+) create mode 100644 installation/mods_MarcMentat/test.py diff --git a/installation/mods_MarcMentat/test.py b/installation/mods_MarcMentat/test.py new file mode 100644 index 000000000..e904651a4 --- /dev/null +++ b/installation/mods_MarcMentat/test.py @@ -0,0 +1 @@ +print('test') From 331cda5da7d04a32b7a806381efd5b96d624d0df Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jun 2017 10:03:07 +0200 Subject: [PATCH 41/62] testing ... --- installation/mods_MarcMentat/test.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100644 => 100755 installation/mods_MarcMentat/test.py diff --git a/installation/mods_MarcMentat/test.py b/installation/mods_MarcMentat/test.py old mode 100644 new mode 100755 From b0d4fe45515ffbf9d9c591d5503f397e4deb405a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jun 2017 10:10:21 +0200 Subject: [PATCH 42/62] still trying --- installation/mods_MarcMentat/test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/installation/mods_MarcMentat/test.py b/installation/mods_MarcMentat/test.py index e904651a4..c419ddb69 100755 --- a/installation/mods_MarcMentat/test.py +++ b/installation/mods_MarcMentat/test.py @@ -1 +1 @@ -print('test') +print('test2') From c9a31274b29bd32b2ff7fd0f27adb78da70570b6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jun 2017 10:10:44 +0200 Subject: [PATCH 43/62] not executable for testing --- installation/mods_MarcMentat/test.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100755 => 100644 installation/mods_MarcMentat/test.py diff --git a/installation/mods_MarcMentat/test.py b/installation/mods_MarcMentat/test.py old mode 100755 new mode 100644 From 10b5864c660f65b5dce5f2e24d918790a0c98d16 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jun 2017 10:12:01 +0200 Subject: [PATCH 44/62] running again for test --- installation/mods_MarcMentat/test.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100644 => 100755 installation/mods_MarcMentat/test.py diff --git a/installation/mods_MarcMentat/test.py b/installation/mods_MarcMentat/test.py old mode 100644 new mode 100755 From eed5adf6cc0417717acf48e9f40abe40d8959a9b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jun 2017 10:14:38 +0200 Subject: [PATCH 45/62] off again --- installation/mods_MarcMentat/test.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100755 => 100644 installation/mods_MarcMentat/test.py diff --git a/installation/mods_MarcMentat/test.py b/installation/mods_MarcMentat/test.py old mode 100755 new mode 100644 From b84e8be7a163ceba4cf1d5f0875cb08dd22e4ed0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jun 2017 10:17:26 +0200 Subject: [PATCH 46/62] on again --- installation/mods_MarcMentat/test.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100644 => 100755 installation/mods_MarcMentat/test.py diff --git a/installation/mods_MarcMentat/test.py b/installation/mods_MarcMentat/test.py old mode 100644 new mode 100755 From 1faea4a0b48918c1c9e30d76493b380d646e41ab Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jun 2017 10:20:26 +0200 Subject: [PATCH 47/62] off again:wq: --- installation/mods_MarcMentat/test.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100755 => 100644 installation/mods_MarcMentat/test.py diff --git a/installation/mods_MarcMentat/test.py b/installation/mods_MarcMentat/test.py old mode 100755 new mode 100644 From 3c56cb9100160bcdcd2ac09d73d2862001e01bb5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jun 2017 10:29:47 +0200 Subject: [PATCH 48/62] on again --- installation/mods_MarcMentat/test.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100644 => 100755 installation/mods_MarcMentat/test.py diff --git a/installation/mods_MarcMentat/test.py b/installation/mods_MarcMentat/test.py old mode 100644 new mode 100755 From ec254f745d140d94adbab098af7ab80a07a66683 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jun 2017 10:31:27 +0200 Subject: [PATCH 49/62] off --- installation/mods_MarcMentat/test.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100755 => 100644 installation/mods_MarcMentat/test.py diff --git a/installation/mods_MarcMentat/test.py b/installation/mods_MarcMentat/test.py old mode 100755 new mode 100644 From 3ffeebfc3a68cdd009d73705324fe07be5f52310 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jun 2017 10:38:42 +0200 Subject: [PATCH 50/62] on --- installation/mods_MarcMentat/test.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100644 => 100755 installation/mods_MarcMentat/test.py diff --git a/installation/mods_MarcMentat/test.py b/installation/mods_MarcMentat/test.py old mode 100644 new mode 100755 From c66bda1e40e2ec67f10ab018c2f36a4fb4eaa756 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jun 2017 10:42:34 +0200 Subject: [PATCH 51/62] off --- installation/mods_MarcMentat/test.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100755 => 100644 installation/mods_MarcMentat/test.py diff --git a/installation/mods_MarcMentat/test.py b/installation/mods_MarcMentat/test.py old mode 100755 new mode 100644 From 28e78ef7761602883174daf1442ea34d113deea4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jun 2017 10:43:31 +0200 Subject: [PATCH 52/62] on --- installation/mods_MarcMentat/test.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100644 => 100755 installation/mods_MarcMentat/test.py diff --git a/installation/mods_MarcMentat/test.py b/installation/mods_MarcMentat/test.py old mode 100644 new mode 100755 From 031df004fa23767dcf3e89604e7bfc5a97244ac1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jun 2017 10:44:11 +0200 Subject: [PATCH 53/62] not needed any more --- installation/mods_MarcMentat/test.py | 1 - 1 file changed, 1 deletion(-) delete mode 100755 installation/mods_MarcMentat/test.py diff --git a/installation/mods_MarcMentat/test.py b/installation/mods_MarcMentat/test.py deleted file mode 100755 index c419ddb69..000000000 --- a/installation/mods_MarcMentat/test.py +++ /dev/null @@ -1 +0,0 @@ -print('test2') From 7a737c149e7f1effb06e9ca10d753c7047bd840f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jun 2017 10:44:56 +0200 Subject: [PATCH 54/62] non executable fiel --- test.sh | 1 + 1 file changed, 1 insertion(+) create mode 100644 test.sh diff --git a/test.sh b/test.sh new file mode 100644 index 000000000..539e9e8e9 --- /dev/null +++ b/test.sh @@ -0,0 +1 @@ +should not be executable From 9792f66a863a9ea4a9a19026615030b2000ad7a2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jun 2017 10:45:20 +0200 Subject: [PATCH 55/62] now executable, should war --- test.sh | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100644 => 100755 test.sh diff --git a/test.sh b/test.sh old mode 100644 new mode 100755 From 47ca37cc28987601a45f3caaf756fb7f67972e65 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jun 2017 10:45:42 +0200 Subject: [PATCH 56/62] everythin seems to work --- test.sh | 1 - 1 file changed, 1 deletion(-) delete mode 100755 test.sh diff --git a/test.sh b/test.sh deleted file mode 100755 index 539e9e8e9..000000000 --- a/test.sh +++ /dev/null @@ -1 +0,0 @@ -should not be executable From 930dd7c70e94e112d16403ed0efd81648e4f26a1 Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 2 Jun 2017 14:21:54 +0200 Subject: [PATCH 57/62] [skip ci] updated version information after successful test of v2.0.1-779-gf085f61 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 1ad92c82c..7272b2d06 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-777-g31859c0 +v2.0.1-779-gf085f61 From 9894af74ca19200956881c095ecb9366d4038852 Mon Sep 17 00:00:00 2001 From: Tias Maiti Date: Tue, 6 Jun 2017 20:37:23 -0400 Subject: [PATCH 58/62] =?UTF-8?q?new=20material=20subroutine=20implementin?= =?UTF-8?q?g=20the=20diagonal=20hardening=20concept=20outlined=20by=20"Z.?= =?UTF-8?q?=20Zhao=20et=20al.=20/=20International=20Journal=20of=20Plastic?= =?UTF-8?q?ity=2024=20(2008)=202278=E2=80=932297"?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/CMakeLists.txt | 1 + src/constitutive.f90 | 27 +- src/material.f90 | 26 +- src/plastic_MITdiag.f90 | 716 ++++++++++++++++++++++++++++++++++++++++ 4 files changed, 755 insertions(+), 15 deletions(-) create mode 100644 src/plastic_MITdiag.f90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index eb1448028..6b32e39ec 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -74,6 +74,7 @@ add_library (PLASTIC OBJECT "plastic_disloUCLA.f90" "plastic_isotropic.f90" "plastic_phenopowerlaw.f90" + "plastic_MITdiag.f90" "plastic_titanmod.f90" "plastic_nonlocal.f90" "plastic_none.f90" diff --git a/src/constitutive.f90 b/src/constitutive.f90 index de8f61c2a..99594a846 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -70,6 +70,7 @@ subroutine constitutive_init() PLASTICITY_none_ID, & PLASTICITY_isotropic_ID, & PLASTICITY_phenopowerlaw_ID, & + PLASTICITY_mitdiag_ID ,& PLASTICITY_phenoplus_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_disloucla_ID, & @@ -93,6 +94,7 @@ subroutine constitutive_init() PLASTICITY_NONE_label, & PLASTICITY_ISOTROPIC_label, & PLASTICITY_PHENOPOWERLAW_label, & + PLASTICITY_MITDIAG_label, & PLASTICITY_PHENOPLUS_label, & PLASTICITY_DISLOTWIN_label, & PLASTICITY_DISLOUCLA_label, & @@ -113,6 +115,7 @@ subroutine constitutive_init() use plastic_none use plastic_isotropic use plastic_phenopowerlaw + use plastic_mitdiag use plastic_phenoplus use plastic_dislotwin use plastic_disloucla @@ -158,6 +161,7 @@ subroutine constitutive_init() if (any(phase_plasticity == PLASTICITY_NONE_ID)) call plastic_none_init if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init(FILEUNIT) + if (any(phase_plasticity == PLASTICITY_MITDIAG_ID)) call plastic_mitdiag_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_PHENOPLUS_ID)) call plastic_phenoplus_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID)) call plastic_disloucla_init(FILEUNIT) @@ -207,7 +211,7 @@ subroutine constitutive_init() outputName = PLASTICITY_NONE_label thisNoutput => null() thisOutput => null() - thisSize => null() + thisSize => null() case (PLASTICITY_ISOTROPIC_ID) plasticityType outputName = PLASTICITY_ISOTROPIC_label thisNoutput => plastic_isotropic_Noutput @@ -218,6 +222,11 @@ subroutine constitutive_init() thisNoutput => plastic_phenopowerlaw_Noutput thisOutput => plastic_phenopowerlaw_output thisSize => plastic_phenopowerlaw_sizePostResult + case (PLASTICITY_MITDIAG_ID) plasticityType + outputName = PLASTICITY_MITDIAG_label + thisNoutput => plastic_mitdiag_Noutput + thisOutput => plastic_mitdiag_output + thisSize => plastic_mitdiag_sizePostResult case (PLASTICITY_PHENOPLUS_ID) plasticityType outputName = PLASTICITY_PHENOPLUS_label thisNoutput => plastic_phenoplus_Noutput @@ -501,6 +510,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v PLASTICITY_NONE_ID, & PLASTICITY_ISOTROPIC_ID, & PLASTICITY_PHENOPOWERLAW_ID, & + PLASTICITY_MITDIAG_ID, & PLASTICITY_PHENOPLUS_ID, & PLASTICITY_DISLOTWIN_ID, & PLASTICITY_DISLOUCLA_ID, & @@ -510,6 +520,8 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v plastic_isotropic_LpAndItsTangent use plastic_phenopowerlaw, only: & plastic_phenopowerlaw_LpAndItsTangent + use plastic_mitdiag, only: & + plastic_mitdiag_LpAndItsTangent use plastic_phenoplus, only: & plastic_phenoplus_LpAndItsTangent use plastic_dislotwin, only: & @@ -560,6 +572,8 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v call plastic_isotropic_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el) + case (PLASTICITY_MITDIAG_ID) plasticityType + call plastic_mitdiag_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el) case (PLASTICITY_PHENOPLUS_ID) plasticityType call plastic_phenoplus_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) plasticityType @@ -884,6 +898,7 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra PLASTICITY_none_ID, & PLASTICITY_isotropic_ID, & PLASTICITY_phenopowerlaw_ID, & + PLASTICITY_mitdiag_ID, & PLASTICITY_phenoplus_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_disloucla_ID, & @@ -897,6 +912,8 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra plastic_isotropic_dotState use plastic_phenopowerlaw, only: & plastic_phenopowerlaw_dotState + use plastic_mitdiag, only: & + plastic_mitdiag_dotState use plastic_phenoplus, only: & plastic_phenoplus_dotState use plastic_dislotwin, only: & @@ -950,6 +967,8 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra call plastic_isotropic_dotState (Tstar_v,ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType call plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) + case (PLASTICITY_MITDIAG_ID) plasticityType + call plastic_mitdiag_dotState(Tstar_v,ipc,ip,el) case (PLASTICITY_PHENOPLUS_ID) plasticityType call plastic_phenoplus_dotState (Tstar_v,ipc,ip,el) case (PLASTICITY_DISLOTWIN_ID) plasticityType @@ -1093,6 +1112,7 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) PLASTICITY_NONE_ID, & PLASTICITY_ISOTROPIC_ID, & PLASTICITY_PHENOPOWERLAW_ID, & + PLASTICITY_MITDIAG_ID, & PLASTICITY_PHENOPLUS_ID, & PLASTICITY_DISLOTWIN_ID, & PLASTICITY_DISLOUCLA_ID, & @@ -1108,6 +1128,8 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) plastic_phenopowerlaw_postResults use plastic_phenoplus, only: & plastic_phenoplus_postResults + use plastic_mitdiag, only: & + plastic_mitdiag_postResults use plastic_dislotwin, only: & plastic_dislotwin_postResults use plastic_disloucla, only: & @@ -1160,6 +1182,9 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType constitutive_postResults(startPos:endPos) = & plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el) + case (PLASTICITY_MITDIAG_ID) plasticityType + constitutive_postResults(startPos:endPos) = & + plastic_mitdiag_postResults(Tstar_v,ipc,ip,el) case (PLASTICITY_PHENOPLUS_ID) plasticityType constitutive_postResults(startPos:endPos) = & plastic_phenoplus_postResults(Tstar_v,ipc,ip,el) diff --git a/src/material.f90 b/src/material.f90 index a77c4871a..1acff4f93 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -25,6 +25,7 @@ module material PLASTICITY_none_label = 'none', & PLASTICITY_isotropic_label = 'isotropic', & PLASTICITY_phenopowerlaw_label = 'phenopowerlaw', & + PLASTICITY_mitdiag_label = 'mitdiag', & PLASTICITY_phenoplus_label = 'phenoplus', & PLASTICITY_dislotwin_label = 'dislotwin', & PLASTICITY_disloucla_label = 'disloucla', & @@ -73,6 +74,7 @@ module material enumerator :: PLASTICITY_undefined_ID, & PLASTICITY_none_ID, & PLASTICITY_isotropic_ID, & + PLASTICITY_mitdiag_ID, & PLASTICITY_phenopowerlaw_ID, & PLASTICITY_phenoplus_ID, & PLASTICITY_dislotwin_ID, & @@ -312,6 +314,7 @@ module material PLASTICITY_none_ID, & PLASTICITY_isotropic_ID, & PLASTICITY_phenopowerlaw_ID, & + PLASTICITY_mitdiag_ID, & PLASTICITY_phenoplus_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_disloucla_ID, & @@ -671,7 +674,7 @@ subroutine material_parseHomogenization(fileUnit,myPart) case ('porosity') select case (IO_lc(IO_stringValue(line,chunkPos,2_pInt))) - case(POROSITY_NONE_label) + case(POROSITY_none_label) porosity_type(section) = POROSITY_none_ID case(POROSITY_phasefield_label) porosity_type(section) = POROSITY_phasefield_ID @@ -729,8 +732,6 @@ end subroutine material_parseHomogenization !> @brief parses the microstructure part in the material configuration file !-------------------------------------------------------------------------------------------------- subroutine material_parseMicrostructure(fileUnit,myPart) - use prec, only: & - dNeq use IO use mesh, only: & mesh_element, & @@ -740,6 +741,7 @@ subroutine material_parseMicrostructure(fileUnit,myPart) character(len=*), intent(in) :: myPart integer(pInt), intent(in) :: fileUnit + integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt) :: Nsections, section, constituent, e, i character(len=65536) :: & @@ -759,7 +761,7 @@ subroutine material_parseMicrostructure(fileUnit,myPart) allocate(microstructure_elemhomo(Nsections), source=.false.) if(any(mesh_element(4,1:mesh_NcpElems) > Nsections)) & - call IO_error(155_pInt,ext_msg='More microstructures in geometry than sections in material.config') + call IO_error(155_pInt,ext_msg='Microstructure in geometry > Sections in material.config') forall (e = 1_pInt:mesh_NcpElems) microstructure_active(mesh_element(4,e)) = .true. ! current microstructure used in model? Elementwise view, maximum N operations for N elements @@ -802,7 +804,7 @@ subroutine material_parseMicrostructure(fileUnit,myPart) microstructure_crystallite(section) = IO_intValue(line,chunkPos,2_pInt) case ('(constituent)') constituent = constituent + 1_pInt - do i = 2_pInt,6_pInt,2_pInt + do i=2_pInt,6_pInt,2_pInt tag = IO_lc(IO_stringValue(line,chunkPos,i)) select case (tag) case('phase') @@ -817,12 +819,6 @@ subroutine material_parseMicrostructure(fileUnit,myPart) endif enddo - !sanity check -do section = 1_pInt, Nsections - if (dNeq(sum(microstructure_fraction(:,section)),1.0_pReal)) & - call IO_error(153_pInt,ext_msg=microstructure_name(section)) -enddo - end subroutine material_parseMicrostructure @@ -979,17 +975,19 @@ subroutine material_parsePhase(fileUnit,myPart) end select case ('plasticity') select case (IO_lc(IO_stringValue(line,chunkPos,2_pInt))) - case (PLASTICITY_NONE_label) - phase_plasticity(section) = PLASTICITY_NONE_ID + case (PLASTICITY_none_label) + phase_plasticity(section) = PLASTICITY_none_ID case (PLASTICITY_ISOTROPIC_label) phase_plasticity(section) = PLASTICITY_ISOTROPIC_ID + case (PLASTICITY_mitdiag_label) + phase_plasticity(section) = PLASTICITY_MITDIAG_ID case (PLASTICITY_PHENOPOWERLAW_label) phase_plasticity(section) = PLASTICITY_PHENOPOWERLAW_ID case (PLASTICITY_PHENOPLUS_label) phase_plasticity(section) = PLASTICITY_PHENOPLUS_ID case (PLASTICITY_DISLOTWIN_label) phase_plasticity(section) = PLASTICITY_DISLOTWIN_ID - case (PLASTICITY_DISLOUCLA_label) + case (PLASTICITY_disloUCLA_label) phase_plasticity(section) = PLASTICITY_DISLOUCLA_ID case (PLASTICITY_TITANMOD_label) phase_plasticity(section) = PLASTICITY_TITANMOD_ID diff --git a/src/plastic_MITdiag.f90 b/src/plastic_MITdiag.f90 new file mode 100644 index 000000000..28d9bc3d7 --- /dev/null +++ b/src/plastic_MITdiag.f90 @@ -0,0 +1,716 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Tias Maiti, Michigan State University +!> @author Philip Eisenlohr, Michigan State University +!> @brief material subroutine implementing the diagonal hardening concept outlined by +!> Raul Radovitzky in the paper +!-------------------------------------------------------------------------------------------------- + +module plastic_MITdiag + use prec, only: & + pReal,& + pInt + use lattice, only: & + lattice_maxNinteraction, & + lattice_maxNslipFamily + + implicit none + private + integer(pInt), dimension(:), allocatable, public, protected :: & + plastic_MITdiag_sizePostResults !< cumulative size of post results + + integer(pInt), dimension(:,:), allocatable, target, public :: & + plastic_MITdiag_sizePostResult !< size of each post result output + + character(len=64), dimension(:,:), allocatable, target, public :: & + plastic_MITdiag_output !< name of each post result output + + integer(pInt), dimension(:), allocatable, target, public :: & + plastic_MITdiag_Noutput !< number of outputs per instance + + integer(pInt), dimension(:), allocatable, public, protected :: & + plastic_MITdiag_totalNslip !< no. of slip system used in simulation + + real(pReal), dimension(:,:,:), allocatable, private :: & + plastic_MITdiag_hardeningMatrix + + enum, bind(c) + enumerator :: undefined_ID, & + resistance_ID, & + accumulatedshear_ID, & + shearrate_ID, & + resolvedstress_ID + end enum + + type, private :: tParameters !< container type for internal constitutive parameters + integer(kind(undefined_ID)), allocatable, dimension(:) :: & + outputID + real(pReal), dimension(lattice_maxNinteraction) :: & + interaction + real(pReal), dimension(:,:), allocatable :: & + hardeningMatrix + integer(pInt), dimension(lattice_maxNslipFamily) :: & + n_slip = 0_pInt + real(pReal), dimension(lattice_maxNslipFamily) :: & + b, & + g0, & + rho0, & + rhosat, & + gammasat + real(pReal) :: & + aTolResistance = 1.0_pReal, & + aTolShear = 1.0e-6_pReal, & + gdot0, & + n + end type + + type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + + type, private :: tMITdiagState !< internal state aliases + real(pReal), pointer, dimension(:,:) :: & ! scalars along NipcMyInstance + resistance, & + accumulatedShear + end type + type, private :: tMITdiagAbsTol !< internal alias for abs tolerance in state + real(pReal), pointer :: & ! scalars along NipcMyInstance + resistance, & + accumulatedShear + end type + type(tMITdiagState), allocatable, dimension(:), private :: & !< state aliases per instance + state, & + state0, & + dotState + type(tMITdiagAbsTol), allocatable, dimension(:), private :: & !< state aliases per instance + stateAbsTol + + public :: & + plastic_MITdiag_init, & + plastic_MITdiag_LpAndItsTangent, & + plastic_MITdiag_dotState, & + plastic_MITdiag_postResults + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +subroutine plastic_MITdiag_init(fileUnit) + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use debug, only: & + debug_level, & + debug_constitutive, & + debug_levelBasic + use numerics, only: & + numerics_integrator + use math, only: & + math_Mandel3333to66, & + math_Voigt66to3333 + use IO, only: & + IO_read, & + IO_lc, & + IO_getTag, & + IO_isBlank, & + IO_stringPos, & + IO_stringValue, & + IO_intValue, & + IO_floatValue, & + IO_error, & + IO_timeStamp, & + IO_EOF, & + IO_warning + use lattice + use material, only: & + phase_plasticity, & + phase_plasticityInstance, & + phase_Noutput, & + PLASTICITY_MITdiag_label, & + PLASTICITY_MITdiag_ID, & + material_phase, & + plasticState, & + MATERIAL_partPhase + + use lattice + + implicit none + integer(pInt), intent(in) :: fileUnit + + + integer(pInt), allocatable, dimension(:) :: chunkPos + integer(pInt) :: & + f,i,j,o,k, & + phase, & + instance, & + maxNinstance, & + mySize = 0_pInt, & + sizeDotState, & + sizeState, & + sizeDeltaState, & + totalNslip, & + index_myFamily, & + index_otherFamily, & + Nchunks_SlipSlip = 0_pInt, & + Nchunks_SlipFamilies = 0_pInt + character(len=65536) :: & + tag = '', & + line = '', & + extmsg = '' + real(pReal), dimension(:), allocatable :: tempPerSlip + character(len=64) :: & + outputtag = '' + integer(pInt) :: NipcMyPhase + + write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_MITdiag_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + + maxNinstance = int(count(phase_plasticity == PLASTICITY_MITdiag_ID),pInt) + if (maxNinstance == 0_pInt) return + + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & + write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + + allocate(plastic_MITdiag_sizePostResults(maxNinstance), source=0_pInt) + allocate(plastic_MITdiag_sizePostResult(maxval(phase_Noutput), maxNinstance),source=0_pInt) + allocate(plastic_MITdiag_output(maxval(phase_Noutput), maxNinstance)) + plastic_MITdiag_output = '' + allocate(plastic_MITdiag_Noutput(maxNinstance), source=0_pInt) + allocate(plastic_MITdiag_totalNslip(maxNinstance), source=0_pInt) + allocate(param(maxNinstance)) ! one container of parameters per instance + + rewind(fileUnit) + phase = 0_pInt + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partPhase) ! wind forward to + line = IO_read(fileUnit) + enddo + + parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part + line = IO_read(fileUnit) + if (IO_isBlank(line)) cycle ! skip empty lines + if (IO_getTag(line,'<','>') /= '') then ! stop at next part + line = IO_read(fileUnit, .true.) ! reset IO_read + exit + endif + if (IO_getTag(line,'[',']') /= '') then ! next section + phase = phase + 1_pInt ! advance section counter + if (phase_plasticity(phase) == PLASTICITY_MITdiag_ID) then + instance = phase_plasticityInstance(phase) ! count instances of my constitutive law + allocate(param(instance)%outputID(phase_Noutput(phase))) ! allocate space for IDs of every requested output + Nchunks_SlipFamilies = count(lattice_NslipSystem(:,phase) > 0_pInt) ! maximum number of slip families according to lattice type of current phase + Nchunks_SlipSlip = maxval(lattice_interactionSlipSlip(:,:,phase)) + if(allocated(tempPerSlip)) deallocate(tempPerSlip) + allocate(tempPerSlip(Nchunks_SlipFamilies)) + endif + cycle ! skip to next line + endif + + if (phase > 0_pInt) then; if (phase_plasticity(phase) == PLASTICITY_MITdiag_ID) then ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran + instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase + chunkPos = IO_stringPos(line) + tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key + extmsg = trim(tag)//' ('//PLASTICITY_MITdiag_label//')' ! prepare error message identifier + + select case(tag) + case ('(output)') + outputtag = IO_lc(IO_stringValue(line,chunkPos,2_pInt)) + select case(outputtag) + case ('resistance') + plastic_MITdiag_Noutput(instance) = plastic_MITdiag_Noutput(instance) + 1_pInt + param(instance)%outputID(plastic_MITdiag_Noutput(instance)) = resistance_ID + plastic_MITdiag_output(plastic_MITdiag_Noutput(instance),instance) = outputtag + case ('accumulatedshear') + plastic_MITdiag_Noutput(instance) = plastic_MITdiag_Noutput(instance) + 1_pInt + param(instance)%outputID(plastic_MITdiag_Noutput(instance)) = accumulatedshear_ID + plastic_MITdiag_output(plastic_MITdiag_Noutput(instance),instance) = outputtag + case ('shearrate') + plastic_MITdiag_Noutput(instance) = plastic_MITdiag_Noutput(instance) + 1_pInt + param(instance)%outputID(plastic_MITdiag_Noutput(instance)) = shearrate_ID + plastic_MITdiag_output(plastic_MITdiag_Noutput(instance),instance) = outputtag + case ('resolvedstress') + plastic_MITdiag_Noutput(instance) = plastic_MITdiag_Noutput(instance) + 1_pInt + param(instance)%outputID(plastic_MITdiag_Noutput(instance)) = resolvedstress_ID + plastic_MITdiag_output(plastic_MITdiag_Noutput(instance),instance) = outputtag + + end select + +!-------------------------------------------------------------------------------------------------- +! parameters depending on number of slip families + case ('nslip') + if (chunkPos(1) < Nchunks_SlipFamilies + 1_pInt) & + call IO_warning(50_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_MITdiag_label//')') + if (chunkPos(1) > Nchunks_SlipFamilies + 1_pInt) & + call IO_error(150_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_MITdiag_label//')') + Nchunks_SlipFamilies = chunkPos(1) - 1_pInt ! user specified number of (possibly) active slip families (e.g. 6 0 6 --> 3) + do j = 1_pInt, Nchunks_SlipFamilies + param(instance)%n_slip(j) = IO_intValue(line,chunkPos,1_pInt+j) + enddo + case ('b','g0','rho0','rhosat','gammasat') + tempPerSlip = 0.0_pReal + do j = 1_pInt, Nchunks_SlipFamilies + if (param(instance)%n_slip(j) > 0_pInt) & + tempPerSlip(j) = IO_floatValue(line,chunkPos,1_pInt+j) + enddo + select case(tag) + case ('b') + param(instance)%b(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('g0') + param(instance)%g0(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('rho0') + param(instance)%rho0(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('rhosat') + param(instance)%rhosat(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('gammasat') + param(instance)%gammasat(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies) + end select + + case ('interaction_slipslip') + if (chunkPos(1) < 1_pInt + Nchunks_SlipSlip) & + call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_MITDIAG_label//')') + do j = 1_pInt, Nchunks_SlipSlip + param(instance)%interaction(j) = IO_floatValue(line,chunkPos,1_pInt+j) + enddo + + case ('gdot0') + param(instance)%gdot0 = IO_floatValue(line,chunkPos,2_pInt) + if (param(instance)%gdot0 <= 0.0_pReal) call IO_error(211_pInt,ext_msg=extmsg) + + case ('n') + param(instance)%n = IO_floatValue(line,chunkPos,2_pInt) + if (param(instance)%n <= 0.0_pReal) call IO_error(211_pInt,ext_msg=extmsg) + + case ('atol_resistance') + param(instance)%aTolResistance = IO_floatValue(line,chunkPos,2_pInt) + if (param(instance)%aTolResistance <= 0.0_pReal) call IO_error(211_pInt,ext_msg=extmsg) + + case ('atol_shear') + param(instance)%aTolShear = IO_floatValue(line,chunkPos,2_pInt) + + case default + + end select + endif; endif + enddo parsingFile + + allocate(state(maxNinstance)) ! internal state aliases + allocate(state0(maxNinstance)) + allocate(dotState(maxNinstance)) + allocate(stateAbsTol(maxNinstance)) + + initializeInstances: do phase = 1_pInt, size(phase_plasticity) ! loop over every plasticity + myPhase: if (phase_plasticity(phase) == PLASTICITY_MITdiag_ID) then ! isolate instances of own constitutive description + NipcMyPhase = count(material_phase == phase) ! number of own material points (including point components ipc) + instance = phase_plasticityInstance(phase) +!-------------------------------------------------------------------------------------------------- +! sanity checks + if (param(instance)%aTolShear <= 0.0_pReal) & + param(instance)%aTolShear = 1.0e-6_pReal ! default absolute tolerance 1e-6 + param(instance)%n_slip(1:lattice_maxNslipFamily) = & + min(lattice_NslipSystem(1:lattice_maxNslipFamily,phase),& ! limit active slip systems per family to min of available and requested + param(instance)%n_slip(1:lattice_maxNslipFamily)) + plastic_MITdiag_totalNslip(instance) = sum(param(instance)%n_slip(1:lattice_maxNslipFamily)) ! how many slip systems altogether + totalNslip = plastic_MITdiag_totalNslip(instance) + + allocate(plastic_MITdiag_hardeningMatrix(maxval(plastic_MITdiag_totalNslip),& ! slip resistance from slip activity + maxval(plastic_MITdiag_totalNslip),& + maxNinstance), source=0.0_pReal) + +!-------------------------------------------------------------------------------------------------- +! Determine size of postResults array + outputsLoop: do o = 1_pInt,plastic_MITdiag_Noutput(instance) + select case(param(instance)%outputID(o)) + case(resistance_ID, & + accumulatedshear_ID, & + shearrate_ID, & + resolvedstress_ID & + ) + mySize = plastic_MITdiag_totalNslip(instance) + case default + end select + + outputFound: if (mySize > 0_pInt) then + plastic_MITdiag_sizePostResult(o,instance) = mySize + plastic_MITdiag_sizePostResults(instance) = & + plastic_MITdiag_sizePostResults(instance) + mySize + endif outputFound + enddo outputsLoop + +!-------------------------------------------------------------------------------------------------- +! allocate state arrays + sizeState = totalNslip & ! g (resistance) + + totalNslip ! accshear + sizeDotState = sizeState ! both evolve + sizeDeltaState = 0_pInt ! no sudden jumps in state + plasticState(phase)%sizeState = sizeState + plasticState(phase)%sizeDotState = sizeDotState + plasticState(phase)%sizeDeltaState = sizeDeltaState + plasticState(phase)%sizePostResults = plastic_MITdiag_sizePostResults(instance) + plasticState(phase)%nSlip = totalNslip + plasticState(phase)%nTwin = 0 + plasticState(phase)%nTrans= 0 + allocate(plasticState(phase)%aTolState ( sizeState)) + + allocate(plasticState(phase)%state0 ( sizeState,NipcMyPhase),source=0.0_pReal) + + allocate(plasticState(phase)%partionedState0 ( sizeState,NipcMyPhase),source=0.0_pReal) + allocate(plasticState(phase)%subState0 ( sizeState,NipcMyPhase),source=0.0_pReal) + allocate(plasticState(phase)%state ( sizeState,NipcMyPhase),source=0.0_pReal) + allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) + allocate(plasticState(phase)%deltaState (sizeDeltaState,NipcMyPhase),source=0.0_pReal) + if (any(numerics_integrator == 1_pInt)) then + allocate(plasticState(phase)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal) + allocate(plasticState(phase)%previousDotState2(sizeDotState,NipcMyPhase),source=0.0_pReal) + endif + if (any(numerics_integrator == 4_pInt)) & + allocate(plasticState(phase)%RK4dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) + if (any(numerics_integrator == 5_pInt)) & + allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NipcMyPhase),source=0.0_pReal) + + do f = 1_pInt,lattice_maxNslipFamily ! >>> interaction slip -- X + index_myFamily = sum(param(instance)%n_slip(1:f-1_pInt)) + do j = 1_pInt,param(instance)%n_slip(f) ! loop over (active) systems in my family (slip) + do o = 1_pInt,lattice_maxNslipFamily + index_otherFamily = sum(param(instance)%n_slip(1:o-1_pInt)) + do k = 1_pInt,param(instance)%n_slip(o) ! loop over (active) systems in other family (slip) + plastic_MITdiag_hardeningMatrix(index_myFamily+j, index_otherFamily+k, instance) = & + param(instance)%interaction(lattice_interactionSlipSlip( & + sum(lattice_NslipSystem(1:f-1,phase))+j, & + sum(lattice_NslipSystem(1:o-1,phase))+k, & + phase)) + enddo + enddo + enddo + enddo + +!-------------------------------------------------------------------------------------------------- +! globally required state aliases + plasticState(phase)%slipRate => plasticState(phase)%dotState(totalNslip+1:2*totalNslip,1:NipcMyPhase) + plasticState(phase)%accumulatedSlip => plasticState(phase)%state (totalNslip+1:2*totalNslip,1:NipcMyPhase) + +!-------------------------------------------------------------------------------------------------- +! locally defined state aliases + state(instance)%resistance => plasticState(phase)%state (1:totalNslip,1:NipcMyPhase) + state0(instance)%resistance => plasticState(phase)%state0 (1:totalNslip,1:NipcMyPhase) + dotState(instance)%resistance => plasticState(phase)%dotState (1:totalNslip,1:NipcMyPhase) + stateAbsTol(instance)%resistance => plasticState(phase)%aTolState(1) + + state(instance)%accumulatedShear => plasticState(phase)%state (totalNslip+1:2*totalNslip,1:NipcMyPhase) + state0(instance)%accumulatedShear => plasticState(phase)%state0 (totalNslip+1:2*totalNslip,1:NipcMyPhase) + dotState(instance)%accumulatedShear => plasticState(phase)%dotState (totalNslip+1:2*totalNslip,1:NipcMyPhase) + stateAbsTol(instance)%accumulatedShear => plasticState(phase)%aTolState(2) + +!-------------------------------------------------------------------------------------------------- +! init state + state0(instance)%resistance = param(instance)%g0(1) + state0(instance)%accumulatedShear = 1.0e-150_pReal + +!-------------------------------------------------------------------------------------------------- +! init absolute state tolerances + stateAbsTol(instance)%resistance = param(instance)%aTolResistance + stateAbsTol(instance)%accumulatedShear = param(instance)%aTolShear + + endif myPhase + enddo initializeInstances + +end subroutine plastic_MITdiag_init + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates plastic velocity gradient and its tangent +!-------------------------------------------------------------------------------------------------- +subroutine plastic_MITdiag_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) + use prec, only: & + dNeq0 + use debug, only: & + debug_level, & + debug_constitutive, & + debug_levelBasic, & + debug_levelExtensive, & + debug_levelSelective, & + debug_e, & + debug_i, & + debug_g + use math, only: & + math_mul6x6, & + math_Mandel6to33, & + math_Plain3333to99, & + math_deviatoric33, & + math_mul33xx33, & + math_transpose33 + use lattice, only: & + lattice_Sslip, & + lattice_Sslip_v, & + lattice_maxNslipFamily, & + lattice_NslipSystem + use material, only: & + phaseAt, phasememberAt, & + material_phase, & + phase_plasticityInstance + + implicit none + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(9,9), intent(out) :: & + dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress + + real(pReal), dimension(6), intent(in) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + + integer(pInt) :: & + instance, & + index_myFamily, & + f,i,j,k,l,m,n, & + of, & + ph + + real(pReal) :: & + tau_slip, & + gdot_slip, & + dgdot_dtauslip + + real(pReal), dimension(3,3,3,3) :: & + dLp_dTstar3333 !< derivative of Lp with respect to Tstar as 4th order tensor + + of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember + ph = phaseAt(ipc,ip,el) + instance = phase_plasticityInstance(material_phase(ipc,ip,el)) + + gdot_slip = 0.0_pReal + Lp = 0.0_pReal + dLp_dTstar3333 = 0.0_pReal + dLp_dTstar99 = 0.0_pReal + +!-------------------------------------------------------------------------------------------------- +! Slip part + j = 0_pInt + slipFamilies: do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family + slipSystems: do i = 1_pInt,param(instance)%n_slip(f) + j = j+1_pInt + + ! Calculation of Lp + tau_slip = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) + + if (dNeq0(abs(tau_slip))) then + gdot_slip = param(instance)%gdot0 * & + ((abs(tau_slip)/(state(instance)%resistance(j,of))) & + **param(instance)%n)*sign(1.0_pReal,tau_slip) + + Lp = Lp + gdot_slip*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph) + end if + + ! Calculation of the tangent of Lp + if (dNeq0(abs(tau_slip))) then + dgdot_dtauslip = gdot_slip * param(instance)%n / tau_slip + 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) + & + dgdot_dtauslip * lattice_Sslip(k,l,1,index_myFamily+i,ph)* & + lattice_Sslip(m,n,1,index_myFamily+i,ph) + endif + + enddo slipSystems + enddo slipFamilies + + dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) + + if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt & + .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & + .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then + write(6,'(/,a,/,3(12x,3(f12.4,1x)/))') '<< CONST isotropic >> Lp', & + math_transpose33(Lp(1:3,1:3)) + write(6,'(/,a,/,f12.5)') '<< CONST phenopowerlaw >> gdot', gdot_slip + end if + +end subroutine plastic_MITdiag_LpAndItsTangent + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +subroutine plastic_MITdiag_dotState(Tstar_v,ipc,ip,el) + use prec, only: & + dNeq0 + use math, only: & + pi + use material, only: & + phaseAt, phasememberAt, & + material_phase, & + phase_plasticityInstance + use lattice, only: & + lattice_mu, & + lattice_Sslip_v, & + lattice_maxNslipFamily, & + lattice_NslipSystem + + implicit none + real(pReal), dimension(6), intent(in):: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + + integer(pInt) :: & + instance,ph, & + nSlip, & + f,i,j, & + index_myFamily, & + of + real(pReal) :: & + tau_slip, & + forresthardening + + real(pReal), dimension(plastic_MITdiag_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + gdot_slip, & + tau_c, & + gamma_c, & + h, rho + + of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember + ph = phaseAt(ipc,ip,el) + instance = phase_plasticityInstance(material_phase(ipc,ip,el)) + nSlip = plastic_MITdiag_totalNslip(instance) + +!-------------------------------------------------------------------------------------------------- +! calculate left and right vectors and calculate dot gammas + gdot_slip = 0.0_pReal + tau_slip = 0.0_pReal + rho = 0.0_pReal + + j = 0_pInt + slipFamilies1: do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family + slipSystems1: do i = 1_pInt,param(instance)%n_slip(f) + j = j+1_pInt + + tau_slip = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) + + if (dNeq0(abs(tau_slip))) then + gdot_slip(j) = param(instance)%gdot0 * & + ((abs(tau_slip)/(state(instance)%resistance(j,of))) & + **param(instance)%n * sign(1.0_pReal,tau_slip)) + + + endif + rho(j) = param(instance)%rhosat(f)* & + (1.0_pReal - (1.0_pReal - (param(instance)%rho0(f) / param(instance)%rhosat(f))) * & + exp(-state(instance)%accumulatedShear(j,of) / param(instance)%gammasat(f))) + + dotState(ph)%accumulatedShear(j,of) = abs(gdot_slip(j)) + + enddo slipSystems1 + enddo slipFamilies1 + + tau_c = 0.0_pReal + gamma_c = 0.0_pReal + h = 0.0_pReal + + + j = 0_pInt + slipFamilies2: do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family + slipSystems2: do i = 1_pInt,param(instance)%n_slip(f) + j = j+1_pInt + + forresthardening = & + dot_product(plastic_MITdiag_hardeningMatrix(j,1:nSlip,instance), rho) + tau_c(j) = 0.3_pReal * lattice_mu(ph) * param(instance)%b(f) * & + sqrt(pi * forresthardening) + gamma_c(j) = param(instance)%b(f) * rho(j) / & + (2.0_pReal * sqrt(forresthardening)) + h(j) = (tau_c(j) / gamma_c(j)) * (state(instance)%resistance(j,of) / tau_c(j))**3.0_pReal * & + (cosh((tau_c(j) / state(instance)%resistance(j,of))**2.0_pReal) - 1.0_pReal) + dotState(ph)%resistance(j,of) = h(j) * abs(gdot_slip(j)) + + enddo slipSystems2 + enddo slipFamilies2 + +end subroutine plastic_MITdiag_dotState + +!-------------------------------------------------------------------------------------------------- +!> @brief return array of constitutive results +!-------------------------------------------------------------------------------------------------- +function plastic_MITdiag_postResults(Tstar_v,ipc,ip,el) + use math, only: & + math_mul6x6 + use material, only: & + material_phase, & + phaseAt, phasememberAt, & + phase_plasticityInstance + use lattice, only: & + lattice_NslipSystem, & + lattice_maxNslipFamily, & + lattice_Sslip_v + + implicit none + real(pReal), dimension(6), intent(in) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), dimension(plastic_MITdiag_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + plastic_MITdiag_postResults + + integer(pInt) :: & + instance,ph, of, & + nSlip, & + o,f,i,c,j, & + index_myFamily + real(pReal) :: & + tau_slip + + of = phasememberAt(ipc,ip,el) + ph = phaseAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember + instance = phase_plasticityInstance(material_phase(ipc,ip,el)) + + nSlip = plastic_MITdiag_totalNslip(instance) + + c = 0_pInt + plastic_MITdiag_postResults = 0.0_pReal + + outputsLoop: do o = 1_pInt,plastic_MITdiag_Noutput(instance) + select case(param(instance)%outputID(o)) + case (resistance_ID) + plastic_MITdiag_postResults(c+1_pInt:c+nSlip) = state(instance)%resistance(:,of) + c = c + nSlip + + case (accumulatedshear_ID) + plastic_MITdiag_postResults(c+1_pInt:c+nSlip) = state(instance)%accumulatedShear(:,of) + c = c + nSlip + + case (shearrate_ID) + j = 0_pInt + slipFamilies1: do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family + slipSystems1: do i = 1_pInt,param(instance)%n_slip(f) + j = j + 1_pInt + tau_slip = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) + plastic_MITdiag_postResults(c+j) = param(instance)%gdot0 * & + ((abs(tau_slip) / state(instance)%resistance(j,of)) ** param(instance)%n & + *sign(1.0_pReal,tau_slip)) + enddo slipSystems1 + enddo slipFamilies1 + c = c + nSlip + + case (resolvedstress_ID) + j = 0_pInt + slipFamilies2: do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family + slipSystems2: do i = 1_pInt,param(instance)%n_slip(f) + j = j + 1_pInt + plastic_MITdiag_postResults(c+j) = & + dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) + enddo slipSystems2 + enddo slipFamilies2 + c = c + nSlip + + end select + enddo outputsLoop + +end function plastic_MITdiag_postResults + + +end module plastic_MITdiag From dd68374afd452c2fe3da1a8c2970b87adf93de85 Mon Sep 17 00:00:00 2001 From: Tias Maiti Date: Sun, 11 Jun 2017 18:47:21 -0400 Subject: [PATCH 59/62] moved new constitutive law to new branch for further testing --- src/CMakeLists.txt | 1 - src/constitutive.f90 | 27 +- src/material.f90 | 26 +- src/plastic_MITdiag.f90 | 716 ---------------------------------------- 4 files changed, 15 insertions(+), 755 deletions(-) delete mode 100644 src/plastic_MITdiag.f90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 6b32e39ec..eb1448028 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -74,7 +74,6 @@ add_library (PLASTIC OBJECT "plastic_disloUCLA.f90" "plastic_isotropic.f90" "plastic_phenopowerlaw.f90" - "plastic_MITdiag.f90" "plastic_titanmod.f90" "plastic_nonlocal.f90" "plastic_none.f90" diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 99594a846..de8f61c2a 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -70,7 +70,6 @@ subroutine constitutive_init() PLASTICITY_none_ID, & PLASTICITY_isotropic_ID, & PLASTICITY_phenopowerlaw_ID, & - PLASTICITY_mitdiag_ID ,& PLASTICITY_phenoplus_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_disloucla_ID, & @@ -94,7 +93,6 @@ subroutine constitutive_init() PLASTICITY_NONE_label, & PLASTICITY_ISOTROPIC_label, & PLASTICITY_PHENOPOWERLAW_label, & - PLASTICITY_MITDIAG_label, & PLASTICITY_PHENOPLUS_label, & PLASTICITY_DISLOTWIN_label, & PLASTICITY_DISLOUCLA_label, & @@ -115,7 +113,6 @@ subroutine constitutive_init() use plastic_none use plastic_isotropic use plastic_phenopowerlaw - use plastic_mitdiag use plastic_phenoplus use plastic_dislotwin use plastic_disloucla @@ -161,7 +158,6 @@ subroutine constitutive_init() if (any(phase_plasticity == PLASTICITY_NONE_ID)) call plastic_none_init if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init(FILEUNIT) - if (any(phase_plasticity == PLASTICITY_MITDIAG_ID)) call plastic_mitdiag_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_PHENOPLUS_ID)) call plastic_phenoplus_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID)) call plastic_disloucla_init(FILEUNIT) @@ -211,7 +207,7 @@ subroutine constitutive_init() outputName = PLASTICITY_NONE_label thisNoutput => null() thisOutput => null() - thisSize => null() + thisSize => null() case (PLASTICITY_ISOTROPIC_ID) plasticityType outputName = PLASTICITY_ISOTROPIC_label thisNoutput => plastic_isotropic_Noutput @@ -222,11 +218,6 @@ subroutine constitutive_init() thisNoutput => plastic_phenopowerlaw_Noutput thisOutput => plastic_phenopowerlaw_output thisSize => plastic_phenopowerlaw_sizePostResult - case (PLASTICITY_MITDIAG_ID) plasticityType - outputName = PLASTICITY_MITDIAG_label - thisNoutput => plastic_mitdiag_Noutput - thisOutput => plastic_mitdiag_output - thisSize => plastic_mitdiag_sizePostResult case (PLASTICITY_PHENOPLUS_ID) plasticityType outputName = PLASTICITY_PHENOPLUS_label thisNoutput => plastic_phenoplus_Noutput @@ -510,7 +501,6 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v PLASTICITY_NONE_ID, & PLASTICITY_ISOTROPIC_ID, & PLASTICITY_PHENOPOWERLAW_ID, & - PLASTICITY_MITDIAG_ID, & PLASTICITY_PHENOPLUS_ID, & PLASTICITY_DISLOTWIN_ID, & PLASTICITY_DISLOUCLA_ID, & @@ -520,8 +510,6 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v plastic_isotropic_LpAndItsTangent use plastic_phenopowerlaw, only: & plastic_phenopowerlaw_LpAndItsTangent - use plastic_mitdiag, only: & - plastic_mitdiag_LpAndItsTangent use plastic_phenoplus, only: & plastic_phenoplus_LpAndItsTangent use plastic_dislotwin, only: & @@ -572,8 +560,6 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v call plastic_isotropic_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el) - case (PLASTICITY_MITDIAG_ID) plasticityType - call plastic_mitdiag_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el) case (PLASTICITY_PHENOPLUS_ID) plasticityType call plastic_phenoplus_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) plasticityType @@ -898,7 +884,6 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra PLASTICITY_none_ID, & PLASTICITY_isotropic_ID, & PLASTICITY_phenopowerlaw_ID, & - PLASTICITY_mitdiag_ID, & PLASTICITY_phenoplus_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_disloucla_ID, & @@ -912,8 +897,6 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra plastic_isotropic_dotState use plastic_phenopowerlaw, only: & plastic_phenopowerlaw_dotState - use plastic_mitdiag, only: & - plastic_mitdiag_dotState use plastic_phenoplus, only: & plastic_phenoplus_dotState use plastic_dislotwin, only: & @@ -967,8 +950,6 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra call plastic_isotropic_dotState (Tstar_v,ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType call plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) - case (PLASTICITY_MITDIAG_ID) plasticityType - call plastic_mitdiag_dotState(Tstar_v,ipc,ip,el) case (PLASTICITY_PHENOPLUS_ID) plasticityType call plastic_phenoplus_dotState (Tstar_v,ipc,ip,el) case (PLASTICITY_DISLOTWIN_ID) plasticityType @@ -1112,7 +1093,6 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) PLASTICITY_NONE_ID, & PLASTICITY_ISOTROPIC_ID, & PLASTICITY_PHENOPOWERLAW_ID, & - PLASTICITY_MITDIAG_ID, & PLASTICITY_PHENOPLUS_ID, & PLASTICITY_DISLOTWIN_ID, & PLASTICITY_DISLOUCLA_ID, & @@ -1128,8 +1108,6 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) plastic_phenopowerlaw_postResults use plastic_phenoplus, only: & plastic_phenoplus_postResults - use plastic_mitdiag, only: & - plastic_mitdiag_postResults use plastic_dislotwin, only: & plastic_dislotwin_postResults use plastic_disloucla, only: & @@ -1182,9 +1160,6 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType constitutive_postResults(startPos:endPos) = & plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el) - case (PLASTICITY_MITDIAG_ID) plasticityType - constitutive_postResults(startPos:endPos) = & - plastic_mitdiag_postResults(Tstar_v,ipc,ip,el) case (PLASTICITY_PHENOPLUS_ID) plasticityType constitutive_postResults(startPos:endPos) = & plastic_phenoplus_postResults(Tstar_v,ipc,ip,el) diff --git a/src/material.f90 b/src/material.f90 index 1acff4f93..a77c4871a 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -25,7 +25,6 @@ module material PLASTICITY_none_label = 'none', & PLASTICITY_isotropic_label = 'isotropic', & PLASTICITY_phenopowerlaw_label = 'phenopowerlaw', & - PLASTICITY_mitdiag_label = 'mitdiag', & PLASTICITY_phenoplus_label = 'phenoplus', & PLASTICITY_dislotwin_label = 'dislotwin', & PLASTICITY_disloucla_label = 'disloucla', & @@ -74,7 +73,6 @@ module material enumerator :: PLASTICITY_undefined_ID, & PLASTICITY_none_ID, & PLASTICITY_isotropic_ID, & - PLASTICITY_mitdiag_ID, & PLASTICITY_phenopowerlaw_ID, & PLASTICITY_phenoplus_ID, & PLASTICITY_dislotwin_ID, & @@ -314,7 +312,6 @@ module material PLASTICITY_none_ID, & PLASTICITY_isotropic_ID, & PLASTICITY_phenopowerlaw_ID, & - PLASTICITY_mitdiag_ID, & PLASTICITY_phenoplus_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_disloucla_ID, & @@ -674,7 +671,7 @@ subroutine material_parseHomogenization(fileUnit,myPart) case ('porosity') select case (IO_lc(IO_stringValue(line,chunkPos,2_pInt))) - case(POROSITY_none_label) + case(POROSITY_NONE_label) porosity_type(section) = POROSITY_none_ID case(POROSITY_phasefield_label) porosity_type(section) = POROSITY_phasefield_ID @@ -732,6 +729,8 @@ end subroutine material_parseHomogenization !> @brief parses the microstructure part in the material configuration file !-------------------------------------------------------------------------------------------------- subroutine material_parseMicrostructure(fileUnit,myPart) + use prec, only: & + dNeq use IO use mesh, only: & mesh_element, & @@ -741,7 +740,6 @@ subroutine material_parseMicrostructure(fileUnit,myPart) character(len=*), intent(in) :: myPart integer(pInt), intent(in) :: fileUnit - integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt) :: Nsections, section, constituent, e, i character(len=65536) :: & @@ -761,7 +759,7 @@ subroutine material_parseMicrostructure(fileUnit,myPart) allocate(microstructure_elemhomo(Nsections), source=.false.) if(any(mesh_element(4,1:mesh_NcpElems) > Nsections)) & - call IO_error(155_pInt,ext_msg='Microstructure in geometry > Sections in material.config') + call IO_error(155_pInt,ext_msg='More microstructures in geometry than sections in material.config') forall (e = 1_pInt:mesh_NcpElems) microstructure_active(mesh_element(4,e)) = .true. ! current microstructure used in model? Elementwise view, maximum N operations for N elements @@ -804,7 +802,7 @@ subroutine material_parseMicrostructure(fileUnit,myPart) microstructure_crystallite(section) = IO_intValue(line,chunkPos,2_pInt) case ('(constituent)') constituent = constituent + 1_pInt - do i=2_pInt,6_pInt,2_pInt + do i = 2_pInt,6_pInt,2_pInt tag = IO_lc(IO_stringValue(line,chunkPos,i)) select case (tag) case('phase') @@ -819,6 +817,12 @@ subroutine material_parseMicrostructure(fileUnit,myPart) endif enddo + !sanity check +do section = 1_pInt, Nsections + if (dNeq(sum(microstructure_fraction(:,section)),1.0_pReal)) & + call IO_error(153_pInt,ext_msg=microstructure_name(section)) +enddo + end subroutine material_parseMicrostructure @@ -975,19 +979,17 @@ subroutine material_parsePhase(fileUnit,myPart) end select case ('plasticity') select case (IO_lc(IO_stringValue(line,chunkPos,2_pInt))) - case (PLASTICITY_none_label) - phase_plasticity(section) = PLASTICITY_none_ID + case (PLASTICITY_NONE_label) + phase_plasticity(section) = PLASTICITY_NONE_ID case (PLASTICITY_ISOTROPIC_label) phase_plasticity(section) = PLASTICITY_ISOTROPIC_ID - case (PLASTICITY_mitdiag_label) - phase_plasticity(section) = PLASTICITY_MITDIAG_ID case (PLASTICITY_PHENOPOWERLAW_label) phase_plasticity(section) = PLASTICITY_PHENOPOWERLAW_ID case (PLASTICITY_PHENOPLUS_label) phase_plasticity(section) = PLASTICITY_PHENOPLUS_ID case (PLASTICITY_DISLOTWIN_label) phase_plasticity(section) = PLASTICITY_DISLOTWIN_ID - case (PLASTICITY_disloUCLA_label) + case (PLASTICITY_DISLOUCLA_label) phase_plasticity(section) = PLASTICITY_DISLOUCLA_ID case (PLASTICITY_TITANMOD_label) phase_plasticity(section) = PLASTICITY_TITANMOD_ID diff --git a/src/plastic_MITdiag.f90 b/src/plastic_MITdiag.f90 deleted file mode 100644 index 28d9bc3d7..000000000 --- a/src/plastic_MITdiag.f90 +++ /dev/null @@ -1,716 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Tias Maiti, Michigan State University -!> @author Philip Eisenlohr, Michigan State University -!> @brief material subroutine implementing the diagonal hardening concept outlined by -!> Raul Radovitzky in the paper -!-------------------------------------------------------------------------------------------------- - -module plastic_MITdiag - use prec, only: & - pReal,& - pInt - use lattice, only: & - lattice_maxNinteraction, & - lattice_maxNslipFamily - - implicit none - private - integer(pInt), dimension(:), allocatable, public, protected :: & - plastic_MITdiag_sizePostResults !< cumulative size of post results - - integer(pInt), dimension(:,:), allocatable, target, public :: & - plastic_MITdiag_sizePostResult !< size of each post result output - - character(len=64), dimension(:,:), allocatable, target, public :: & - plastic_MITdiag_output !< name of each post result output - - integer(pInt), dimension(:), allocatable, target, public :: & - plastic_MITdiag_Noutput !< number of outputs per instance - - integer(pInt), dimension(:), allocatable, public, protected :: & - plastic_MITdiag_totalNslip !< no. of slip system used in simulation - - real(pReal), dimension(:,:,:), allocatable, private :: & - plastic_MITdiag_hardeningMatrix - - enum, bind(c) - enumerator :: undefined_ID, & - resistance_ID, & - accumulatedshear_ID, & - shearrate_ID, & - resolvedstress_ID - end enum - - type, private :: tParameters !< container type for internal constitutive parameters - integer(kind(undefined_ID)), allocatable, dimension(:) :: & - outputID - real(pReal), dimension(lattice_maxNinteraction) :: & - interaction - real(pReal), dimension(:,:), allocatable :: & - hardeningMatrix - integer(pInt), dimension(lattice_maxNslipFamily) :: & - n_slip = 0_pInt - real(pReal), dimension(lattice_maxNslipFamily) :: & - b, & - g0, & - rho0, & - rhosat, & - gammasat - real(pReal) :: & - aTolResistance = 1.0_pReal, & - aTolShear = 1.0e-6_pReal, & - gdot0, & - n - end type - - type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) - - type, private :: tMITdiagState !< internal state aliases - real(pReal), pointer, dimension(:,:) :: & ! scalars along NipcMyInstance - resistance, & - accumulatedShear - end type - type, private :: tMITdiagAbsTol !< internal alias for abs tolerance in state - real(pReal), pointer :: & ! scalars along NipcMyInstance - resistance, & - accumulatedShear - end type - type(tMITdiagState), allocatable, dimension(:), private :: & !< state aliases per instance - state, & - state0, & - dotState - type(tMITdiagAbsTol), allocatable, dimension(:), private :: & !< state aliases per instance - stateAbsTol - - public :: & - plastic_MITdiag_init, & - plastic_MITdiag_LpAndItsTangent, & - plastic_MITdiag_dotState, & - plastic_MITdiag_postResults - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine plastic_MITdiag_init(fileUnit) - use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) - use debug, only: & - debug_level, & - debug_constitutive, & - debug_levelBasic - use numerics, only: & - numerics_integrator - use math, only: & - math_Mandel3333to66, & - math_Voigt66to3333 - use IO, only: & - IO_read, & - IO_lc, & - IO_getTag, & - IO_isBlank, & - IO_stringPos, & - IO_stringValue, & - IO_intValue, & - IO_floatValue, & - IO_error, & - IO_timeStamp, & - IO_EOF, & - IO_warning - use lattice - use material, only: & - phase_plasticity, & - phase_plasticityInstance, & - phase_Noutput, & - PLASTICITY_MITdiag_label, & - PLASTICITY_MITdiag_ID, & - material_phase, & - plasticState, & - MATERIAL_partPhase - - use lattice - - implicit none - integer(pInt), intent(in) :: fileUnit - - - integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: & - f,i,j,o,k, & - phase, & - instance, & - maxNinstance, & - mySize = 0_pInt, & - sizeDotState, & - sizeState, & - sizeDeltaState, & - totalNslip, & - index_myFamily, & - index_otherFamily, & - Nchunks_SlipSlip = 0_pInt, & - Nchunks_SlipFamilies = 0_pInt - character(len=65536) :: & - tag = '', & - line = '', & - extmsg = '' - real(pReal), dimension(:), allocatable :: tempPerSlip - character(len=64) :: & - outputtag = '' - integer(pInt) :: NipcMyPhase - - write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_MITdiag_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" - - maxNinstance = int(count(phase_plasticity == PLASTICITY_MITdiag_ID),pInt) - if (maxNinstance == 0_pInt) return - - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance - - allocate(plastic_MITdiag_sizePostResults(maxNinstance), source=0_pInt) - allocate(plastic_MITdiag_sizePostResult(maxval(phase_Noutput), maxNinstance),source=0_pInt) - allocate(plastic_MITdiag_output(maxval(phase_Noutput), maxNinstance)) - plastic_MITdiag_output = '' - allocate(plastic_MITdiag_Noutput(maxNinstance), source=0_pInt) - allocate(plastic_MITdiag_totalNslip(maxNinstance), source=0_pInt) - allocate(param(maxNinstance)) ! one container of parameters per instance - - rewind(fileUnit) - phase = 0_pInt - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partPhase) ! wind forward to - line = IO_read(fileUnit) - enddo - - parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part - line = IO_read(fileUnit) - if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') then ! stop at next part - line = IO_read(fileUnit, .true.) ! reset IO_read - exit - endif - if (IO_getTag(line,'[',']') /= '') then ! next section - phase = phase + 1_pInt ! advance section counter - if (phase_plasticity(phase) == PLASTICITY_MITdiag_ID) then - instance = phase_plasticityInstance(phase) ! count instances of my constitutive law - allocate(param(instance)%outputID(phase_Noutput(phase))) ! allocate space for IDs of every requested output - Nchunks_SlipFamilies = count(lattice_NslipSystem(:,phase) > 0_pInt) ! maximum number of slip families according to lattice type of current phase - Nchunks_SlipSlip = maxval(lattice_interactionSlipSlip(:,:,phase)) - if(allocated(tempPerSlip)) deallocate(tempPerSlip) - allocate(tempPerSlip(Nchunks_SlipFamilies)) - endif - cycle ! skip to next line - endif - - if (phase > 0_pInt) then; if (phase_plasticity(phase) == PLASTICITY_MITdiag_ID) then ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran - instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase - chunkPos = IO_stringPos(line) - tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key - extmsg = trim(tag)//' ('//PLASTICITY_MITdiag_label//')' ! prepare error message identifier - - select case(tag) - case ('(output)') - outputtag = IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - select case(outputtag) - case ('resistance') - plastic_MITdiag_Noutput(instance) = plastic_MITdiag_Noutput(instance) + 1_pInt - param(instance)%outputID(plastic_MITdiag_Noutput(instance)) = resistance_ID - plastic_MITdiag_output(plastic_MITdiag_Noutput(instance),instance) = outputtag - case ('accumulatedshear') - plastic_MITdiag_Noutput(instance) = plastic_MITdiag_Noutput(instance) + 1_pInt - param(instance)%outputID(plastic_MITdiag_Noutput(instance)) = accumulatedshear_ID - plastic_MITdiag_output(plastic_MITdiag_Noutput(instance),instance) = outputtag - case ('shearrate') - plastic_MITdiag_Noutput(instance) = plastic_MITdiag_Noutput(instance) + 1_pInt - param(instance)%outputID(plastic_MITdiag_Noutput(instance)) = shearrate_ID - plastic_MITdiag_output(plastic_MITdiag_Noutput(instance),instance) = outputtag - case ('resolvedstress') - plastic_MITdiag_Noutput(instance) = plastic_MITdiag_Noutput(instance) + 1_pInt - param(instance)%outputID(plastic_MITdiag_Noutput(instance)) = resolvedstress_ID - plastic_MITdiag_output(plastic_MITdiag_Noutput(instance),instance) = outputtag - - end select - -!-------------------------------------------------------------------------------------------------- -! parameters depending on number of slip families - case ('nslip') - if (chunkPos(1) < Nchunks_SlipFamilies + 1_pInt) & - call IO_warning(50_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_MITdiag_label//')') - if (chunkPos(1) > Nchunks_SlipFamilies + 1_pInt) & - call IO_error(150_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_MITdiag_label//')') - Nchunks_SlipFamilies = chunkPos(1) - 1_pInt ! user specified number of (possibly) active slip families (e.g. 6 0 6 --> 3) - do j = 1_pInt, Nchunks_SlipFamilies - param(instance)%n_slip(j) = IO_intValue(line,chunkPos,1_pInt+j) - enddo - case ('b','g0','rho0','rhosat','gammasat') - tempPerSlip = 0.0_pReal - do j = 1_pInt, Nchunks_SlipFamilies - if (param(instance)%n_slip(j) > 0_pInt) & - tempPerSlip(j) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo - select case(tag) - case ('b') - param(instance)%b(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies) - case ('g0') - param(instance)%g0(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies) - case ('rho0') - param(instance)%rho0(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies) - case ('rhosat') - param(instance)%rhosat(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies) - case ('gammasat') - param(instance)%gammasat(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies) - end select - - case ('interaction_slipslip') - if (chunkPos(1) < 1_pInt + Nchunks_SlipSlip) & - call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_MITDIAG_label//')') - do j = 1_pInt, Nchunks_SlipSlip - param(instance)%interaction(j) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo - - case ('gdot0') - param(instance)%gdot0 = IO_floatValue(line,chunkPos,2_pInt) - if (param(instance)%gdot0 <= 0.0_pReal) call IO_error(211_pInt,ext_msg=extmsg) - - case ('n') - param(instance)%n = IO_floatValue(line,chunkPos,2_pInt) - if (param(instance)%n <= 0.0_pReal) call IO_error(211_pInt,ext_msg=extmsg) - - case ('atol_resistance') - param(instance)%aTolResistance = IO_floatValue(line,chunkPos,2_pInt) - if (param(instance)%aTolResistance <= 0.0_pReal) call IO_error(211_pInt,ext_msg=extmsg) - - case ('atol_shear') - param(instance)%aTolShear = IO_floatValue(line,chunkPos,2_pInt) - - case default - - end select - endif; endif - enddo parsingFile - - allocate(state(maxNinstance)) ! internal state aliases - allocate(state0(maxNinstance)) - allocate(dotState(maxNinstance)) - allocate(stateAbsTol(maxNinstance)) - - initializeInstances: do phase = 1_pInt, size(phase_plasticity) ! loop over every plasticity - myPhase: if (phase_plasticity(phase) == PLASTICITY_MITdiag_ID) then ! isolate instances of own constitutive description - NipcMyPhase = count(material_phase == phase) ! number of own material points (including point components ipc) - instance = phase_plasticityInstance(phase) -!-------------------------------------------------------------------------------------------------- -! sanity checks - if (param(instance)%aTolShear <= 0.0_pReal) & - param(instance)%aTolShear = 1.0e-6_pReal ! default absolute tolerance 1e-6 - param(instance)%n_slip(1:lattice_maxNslipFamily) = & - min(lattice_NslipSystem(1:lattice_maxNslipFamily,phase),& ! limit active slip systems per family to min of available and requested - param(instance)%n_slip(1:lattice_maxNslipFamily)) - plastic_MITdiag_totalNslip(instance) = sum(param(instance)%n_slip(1:lattice_maxNslipFamily)) ! how many slip systems altogether - totalNslip = plastic_MITdiag_totalNslip(instance) - - allocate(plastic_MITdiag_hardeningMatrix(maxval(plastic_MITdiag_totalNslip),& ! slip resistance from slip activity - maxval(plastic_MITdiag_totalNslip),& - maxNinstance), source=0.0_pReal) - -!-------------------------------------------------------------------------------------------------- -! Determine size of postResults array - outputsLoop: do o = 1_pInt,plastic_MITdiag_Noutput(instance) - select case(param(instance)%outputID(o)) - case(resistance_ID, & - accumulatedshear_ID, & - shearrate_ID, & - resolvedstress_ID & - ) - mySize = plastic_MITdiag_totalNslip(instance) - case default - end select - - outputFound: if (mySize > 0_pInt) then - plastic_MITdiag_sizePostResult(o,instance) = mySize - plastic_MITdiag_sizePostResults(instance) = & - plastic_MITdiag_sizePostResults(instance) + mySize - endif outputFound - enddo outputsLoop - -!-------------------------------------------------------------------------------------------------- -! allocate state arrays - sizeState = totalNslip & ! g (resistance) - + totalNslip ! accshear - sizeDotState = sizeState ! both evolve - sizeDeltaState = 0_pInt ! no sudden jumps in state - plasticState(phase)%sizeState = sizeState - plasticState(phase)%sizeDotState = sizeDotState - plasticState(phase)%sizeDeltaState = sizeDeltaState - plasticState(phase)%sizePostResults = plastic_MITdiag_sizePostResults(instance) - plasticState(phase)%nSlip = totalNslip - plasticState(phase)%nTwin = 0 - plasticState(phase)%nTrans= 0 - allocate(plasticState(phase)%aTolState ( sizeState)) - - allocate(plasticState(phase)%state0 ( sizeState,NipcMyPhase),source=0.0_pReal) - - allocate(plasticState(phase)%partionedState0 ( sizeState,NipcMyPhase),source=0.0_pReal) - allocate(plasticState(phase)%subState0 ( sizeState,NipcMyPhase),source=0.0_pReal) - allocate(plasticState(phase)%state ( sizeState,NipcMyPhase),source=0.0_pReal) - allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - allocate(plasticState(phase)%deltaState (sizeDeltaState,NipcMyPhase),source=0.0_pReal) - if (any(numerics_integrator == 1_pInt)) then - allocate(plasticState(phase)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - allocate(plasticState(phase)%previousDotState2(sizeDotState,NipcMyPhase),source=0.0_pReal) - endif - if (any(numerics_integrator == 4_pInt)) & - allocate(plasticState(phase)%RK4dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - if (any(numerics_integrator == 5_pInt)) & - allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NipcMyPhase),source=0.0_pReal) - - do f = 1_pInt,lattice_maxNslipFamily ! >>> interaction slip -- X - index_myFamily = sum(param(instance)%n_slip(1:f-1_pInt)) - do j = 1_pInt,param(instance)%n_slip(f) ! loop over (active) systems in my family (slip) - do o = 1_pInt,lattice_maxNslipFamily - index_otherFamily = sum(param(instance)%n_slip(1:o-1_pInt)) - do k = 1_pInt,param(instance)%n_slip(o) ! loop over (active) systems in other family (slip) - plastic_MITdiag_hardeningMatrix(index_myFamily+j, index_otherFamily+k, instance) = & - param(instance)%interaction(lattice_interactionSlipSlip( & - sum(lattice_NslipSystem(1:f-1,phase))+j, & - sum(lattice_NslipSystem(1:o-1,phase))+k, & - phase)) - enddo - enddo - enddo - enddo - -!-------------------------------------------------------------------------------------------------- -! globally required state aliases - plasticState(phase)%slipRate => plasticState(phase)%dotState(totalNslip+1:2*totalNslip,1:NipcMyPhase) - plasticState(phase)%accumulatedSlip => plasticState(phase)%state (totalNslip+1:2*totalNslip,1:NipcMyPhase) - -!-------------------------------------------------------------------------------------------------- -! locally defined state aliases - state(instance)%resistance => plasticState(phase)%state (1:totalNslip,1:NipcMyPhase) - state0(instance)%resistance => plasticState(phase)%state0 (1:totalNslip,1:NipcMyPhase) - dotState(instance)%resistance => plasticState(phase)%dotState (1:totalNslip,1:NipcMyPhase) - stateAbsTol(instance)%resistance => plasticState(phase)%aTolState(1) - - state(instance)%accumulatedShear => plasticState(phase)%state (totalNslip+1:2*totalNslip,1:NipcMyPhase) - state0(instance)%accumulatedShear => plasticState(phase)%state0 (totalNslip+1:2*totalNslip,1:NipcMyPhase) - dotState(instance)%accumulatedShear => plasticState(phase)%dotState (totalNslip+1:2*totalNslip,1:NipcMyPhase) - stateAbsTol(instance)%accumulatedShear => plasticState(phase)%aTolState(2) - -!-------------------------------------------------------------------------------------------------- -! init state - state0(instance)%resistance = param(instance)%g0(1) - state0(instance)%accumulatedShear = 1.0e-150_pReal - -!-------------------------------------------------------------------------------------------------- -! init absolute state tolerances - stateAbsTol(instance)%resistance = param(instance)%aTolResistance - stateAbsTol(instance)%accumulatedShear = param(instance)%aTolShear - - endif myPhase - enddo initializeInstances - -end subroutine plastic_MITdiag_init - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent -!-------------------------------------------------------------------------------------------------- -subroutine plastic_MITdiag_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) - use prec, only: & - dNeq0 - use debug, only: & - debug_level, & - debug_constitutive, & - debug_levelBasic, & - debug_levelExtensive, & - debug_levelSelective, & - debug_e, & - debug_i, & - debug_g - use math, only: & - math_mul6x6, & - math_Mandel6to33, & - math_Plain3333to99, & - math_deviatoric33, & - math_mul33xx33, & - math_transpose33 - use lattice, only: & - lattice_Sslip, & - lattice_Sslip_v, & - lattice_maxNslipFamily, & - lattice_NslipSystem - use material, only: & - phaseAt, phasememberAt, & - material_phase, & - phase_plasticityInstance - - implicit none - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(9,9), intent(out) :: & - dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress - - real(pReal), dimension(6), intent(in) :: & - Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - - integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - - integer(pInt) :: & - instance, & - index_myFamily, & - f,i,j,k,l,m,n, & - of, & - ph - - real(pReal) :: & - tau_slip, & - gdot_slip, & - dgdot_dtauslip - - real(pReal), dimension(3,3,3,3) :: & - dLp_dTstar3333 !< derivative of Lp with respect to Tstar as 4th order tensor - - of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember - ph = phaseAt(ipc,ip,el) - instance = phase_plasticityInstance(material_phase(ipc,ip,el)) - - gdot_slip = 0.0_pReal - Lp = 0.0_pReal - dLp_dTstar3333 = 0.0_pReal - dLp_dTstar99 = 0.0_pReal - -!-------------------------------------------------------------------------------------------------- -! Slip part - j = 0_pInt - slipFamilies: do f = 1_pInt,lattice_maxNslipFamily - index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family - slipSystems: do i = 1_pInt,param(instance)%n_slip(f) - j = j+1_pInt - - ! Calculation of Lp - tau_slip = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) - - if (dNeq0(abs(tau_slip))) then - gdot_slip = param(instance)%gdot0 * & - ((abs(tau_slip)/(state(instance)%resistance(j,of))) & - **param(instance)%n)*sign(1.0_pReal,tau_slip) - - Lp = Lp + gdot_slip*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph) - end if - - ! Calculation of the tangent of Lp - if (dNeq0(abs(tau_slip))) then - dgdot_dtauslip = gdot_slip * param(instance)%n / tau_slip - 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) + & - dgdot_dtauslip * lattice_Sslip(k,l,1,index_myFamily+i,ph)* & - lattice_Sslip(m,n,1,index_myFamily+i,ph) - endif - - enddo slipSystems - enddo slipFamilies - - dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) - - if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt & - .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & - .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then - write(6,'(/,a,/,3(12x,3(f12.4,1x)/))') '<< CONST isotropic >> Lp', & - math_transpose33(Lp(1:3,1:3)) - write(6,'(/,a,/,f12.5)') '<< CONST phenopowerlaw >> gdot', gdot_slip - end if - -end subroutine plastic_MITdiag_LpAndItsTangent - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates the rate of change of microstructure -!-------------------------------------------------------------------------------------------------- -subroutine plastic_MITdiag_dotState(Tstar_v,ipc,ip,el) - use prec, only: & - dNeq0 - use math, only: & - pi - use material, only: & - phaseAt, phasememberAt, & - material_phase, & - phase_plasticityInstance - use lattice, only: & - lattice_mu, & - lattice_Sslip_v, & - lattice_maxNslipFamily, & - lattice_NslipSystem - - implicit none - real(pReal), dimension(6), intent(in):: & - Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - - integer(pInt) :: & - instance,ph, & - nSlip, & - f,i,j, & - index_myFamily, & - of - real(pReal) :: & - tau_slip, & - forresthardening - - real(pReal), dimension(plastic_MITdiag_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & - gdot_slip, & - tau_c, & - gamma_c, & - h, rho - - of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember - ph = phaseAt(ipc,ip,el) - instance = phase_plasticityInstance(material_phase(ipc,ip,el)) - nSlip = plastic_MITdiag_totalNslip(instance) - -!-------------------------------------------------------------------------------------------------- -! calculate left and right vectors and calculate dot gammas - gdot_slip = 0.0_pReal - tau_slip = 0.0_pReal - rho = 0.0_pReal - - j = 0_pInt - slipFamilies1: do f = 1_pInt,lattice_maxNslipFamily - index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family - slipSystems1: do i = 1_pInt,param(instance)%n_slip(f) - j = j+1_pInt - - tau_slip = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) - - if (dNeq0(abs(tau_slip))) then - gdot_slip(j) = param(instance)%gdot0 * & - ((abs(tau_slip)/(state(instance)%resistance(j,of))) & - **param(instance)%n * sign(1.0_pReal,tau_slip)) - - - endif - rho(j) = param(instance)%rhosat(f)* & - (1.0_pReal - (1.0_pReal - (param(instance)%rho0(f) / param(instance)%rhosat(f))) * & - exp(-state(instance)%accumulatedShear(j,of) / param(instance)%gammasat(f))) - - dotState(ph)%accumulatedShear(j,of) = abs(gdot_slip(j)) - - enddo slipSystems1 - enddo slipFamilies1 - - tau_c = 0.0_pReal - gamma_c = 0.0_pReal - h = 0.0_pReal - - - j = 0_pInt - slipFamilies2: do f = 1_pInt,lattice_maxNslipFamily - index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family - slipSystems2: do i = 1_pInt,param(instance)%n_slip(f) - j = j+1_pInt - - forresthardening = & - dot_product(plastic_MITdiag_hardeningMatrix(j,1:nSlip,instance), rho) - tau_c(j) = 0.3_pReal * lattice_mu(ph) * param(instance)%b(f) * & - sqrt(pi * forresthardening) - gamma_c(j) = param(instance)%b(f) * rho(j) / & - (2.0_pReal * sqrt(forresthardening)) - h(j) = (tau_c(j) / gamma_c(j)) * (state(instance)%resistance(j,of) / tau_c(j))**3.0_pReal * & - (cosh((tau_c(j) / state(instance)%resistance(j,of))**2.0_pReal) - 1.0_pReal) - dotState(ph)%resistance(j,of) = h(j) * abs(gdot_slip(j)) - - enddo slipSystems2 - enddo slipFamilies2 - -end subroutine plastic_MITdiag_dotState - -!-------------------------------------------------------------------------------------------------- -!> @brief return array of constitutive results -!-------------------------------------------------------------------------------------------------- -function plastic_MITdiag_postResults(Tstar_v,ipc,ip,el) - use math, only: & - math_mul6x6 - use material, only: & - material_phase, & - phaseAt, phasememberAt, & - phase_plasticityInstance - use lattice, only: & - lattice_NslipSystem, & - lattice_maxNslipFamily, & - lattice_Sslip_v - - implicit none - real(pReal), dimension(6), intent(in) :: & - Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), dimension(plastic_MITdiag_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & - plastic_MITdiag_postResults - - integer(pInt) :: & - instance,ph, of, & - nSlip, & - o,f,i,c,j, & - index_myFamily - real(pReal) :: & - tau_slip - - of = phasememberAt(ipc,ip,el) - ph = phaseAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember - instance = phase_plasticityInstance(material_phase(ipc,ip,el)) - - nSlip = plastic_MITdiag_totalNslip(instance) - - c = 0_pInt - plastic_MITdiag_postResults = 0.0_pReal - - outputsLoop: do o = 1_pInt,plastic_MITdiag_Noutput(instance) - select case(param(instance)%outputID(o)) - case (resistance_ID) - plastic_MITdiag_postResults(c+1_pInt:c+nSlip) = state(instance)%resistance(:,of) - c = c + nSlip - - case (accumulatedshear_ID) - plastic_MITdiag_postResults(c+1_pInt:c+nSlip) = state(instance)%accumulatedShear(:,of) - c = c + nSlip - - case (shearrate_ID) - j = 0_pInt - slipFamilies1: do f = 1_pInt,lattice_maxNslipFamily - index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family - slipSystems1: do i = 1_pInt,param(instance)%n_slip(f) - j = j + 1_pInt - tau_slip = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) - plastic_MITdiag_postResults(c+j) = param(instance)%gdot0 * & - ((abs(tau_slip) / state(instance)%resistance(j,of)) ** param(instance)%n & - *sign(1.0_pReal,tau_slip)) - enddo slipSystems1 - enddo slipFamilies1 - c = c + nSlip - - case (resolvedstress_ID) - j = 0_pInt - slipFamilies2: do f = 1_pInt,lattice_maxNslipFamily - index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family - slipSystems2: do i = 1_pInt,param(instance)%n_slip(f) - j = j + 1_pInt - plastic_MITdiag_postResults(c+j) = & - dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) - enddo slipSystems2 - enddo slipFamilies2 - c = c + nSlip - - end select - enddo outputsLoop - -end function plastic_MITdiag_postResults - - -end module plastic_MITdiag From 68d0b5c6b4aff6c63644cb5a58cb04eda0b866ea Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 12 Jun 2017 05:40:03 +0200 Subject: [PATCH 60/62] [skip ci] updated version information after successful test of v2.0.1-803-gdd68374 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 7272b2d06..85d08b940 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-779-gf085f61 +v2.0.1-803-gdd68374 From dc3eda336d0e871575e4b7e52ba7a2ad81b208d0 Mon Sep 17 00:00:00 2001 From: Franz Roters Date: Fri, 14 Jul 2017 15:28:09 +0200 Subject: [PATCH 61/62] corrected unit of atomic volume to b^3 --- examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.config | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.config b/examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.config index e7c9d4e19..0330e6f4a 100644 --- a/examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.config +++ b/examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.config @@ -43,7 +43,7 @@ q_slip 1.0 # q-exponent in glide velocity CLambdaSlip 10.0 # Adj. parameter controlling dislocation mean free path D0 4.0e-5 # Vacancy diffusion prefactor [m**2/s] Qsd 4.5e-19 # Activation energy for climb [J] -Catomicvolume 1.0 # Adj. parameter controlling the atomic volume [in b] +Catomicvolume 1.0 # Adj. parameter controlling the atomic volume [in b^3] Cedgedipmindistance 1.0 # Adj. parameter controlling the minimum dipole distance [in b] atol_rho 1.0 interactionSlipSlip 0.122 0.122 0.625 0.07 0.137 0.122 # Interaction coefficients (Kubin et al. 2008) From a4ceebc00b582d342c27cf38886e57d07e9e6821 Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 14 Jul 2017 20:26:55 +0200 Subject: [PATCH 62/62] [skip ci] updated version information after successful test of v2.0.1-805-gdc3eda3 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 85d08b940..12e0d6df0 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-803-gdd68374 +v2.0.1-805-gdc3eda3