diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index 4505fb533..4c19bad4d 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -1,4 +1,4 @@ -! Copyright 2011-13 Max-Planck-Institut für Eisenforschung GmbH +!nd Copyright 2011-13 Max-Planck-Institut für Eisenforschung GmbH ! ! This file is part of DAMASK, ! the Düsseldorf Advanced Material Simulation Kit. @@ -37,8 +37,8 @@ module CPFEM real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_dcsdE_knownGood !< known good tangent logical :: CPFEM_init_done = .false., & !< remember whether init has been done already - CPFEM_init_inProgress = .false., & !< remember whether first IP is currently performing init - CPFEM_calc_done = .false. !< remember whether first IP has already calced the results + CPFEM_init_inProgress = .false., & !< remember whether first ip is currently performing init + CPFEM_calc_done = .false. !< remember whether first ip has already calced the results integer(pInt), parameter, public :: & CPFEM_CALCRESULTS = 2_pInt**0_pInt, & CPFEM_AGERESULTS = 2_pInt**1_pInt, & @@ -55,7 +55,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief call (thread safe) all module initializations !-------------------------------------------------------------------------------------------------- -subroutine CPFEM_initAll(Temperature,element,IP) +subroutine CPFEM_initAll(temperature,el,ip) use prec, only: & prec_init use numerics, only: & @@ -83,9 +83,9 @@ subroutine CPFEM_initAll(Temperature,element,IP) use DAMASK_interface implicit none - integer(pInt), intent(in) :: element, & ! FE element number - IP ! FE integration point number - real(pReal), intent(in) :: Temperature ! temperature + integer(pInt), intent(in) :: el, & ! FE el number + ip ! FE integration point number + real(pReal), intent(in) :: temperature ! temperature real(pReal) rnd integer(pInt) i,n @@ -108,12 +108,12 @@ subroutine CPFEM_initAll(Temperature,element,IP) call debug_init call math_init call FE_init - call mesh_init(IP, element) ! pass on coordinates to alter calcMode of first ip + call mesh_init(ip, el) ! pass on coordinates to alter calcMode of first ip call lattice_init call material_init call constitutive_init - call crystallite_init(Temperature) ! (have to) use temperature of first IP for whole model - call homogenization_init(Temperature) + call crystallite_init(temperature) ! (have to) use temperature of first ip for whole model + call homogenization_init call CPFEM_init #ifndef Spectral call DAMASK_interface_init() ! Spectral solver init is already done @@ -172,7 +172,7 @@ subroutine CPFEM_init homogenization_state0 implicit none - integer(pInt) i,j,k,l,m + integer(pInt) :: i,j,k,l,m write(6,'(/,a)') ' <<<+- CPFEM init -+>>>' write(6,'(a)') ' $Id$' @@ -243,10 +243,9 @@ subroutine CPFEM_init endif if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) then - write(6,'(a32,1x,6(i8,1x))') 'CPFEM_cs: ', shape(CPFEM_cs) - write(6,'(a32,1x,6(i8,1x))') 'CPFEM_dcsdE: ', shape(CPFEM_dcsdE) - write(6,'(a32,1x,6(i8,1x))') 'CPFEM_dcsdE_knownGood: ', shape(CPFEM_dcsdE_knownGood) - write(6,*) + write(6,'(a32,1x,6(i8,1x))') 'CPFEM_cs: ', shape(CPFEM_cs) + write(6,'(a32,1x,6(i8,1x))') 'CPFEM_dcsdE: ', shape(CPFEM_dcsdE) + write(6,'(a32,1x,6(i8,1x),/)') 'CPFEM_dcsdE_knownGood: ', shape(CPFEM_dcsdE_knownGood) write(6,*) 'symmetricSolver: ', symmetricSolver endif flush(6) @@ -257,81 +256,90 @@ end subroutine CPFEM_init !-------------------------------------------------------------------------------------------------- !> @brief perform initialization at first call, update variables and call the actual material model !-------------------------------------------------------------------------------------------------- -subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, element, IP, cauchyStress, jacobian) - ! note: cauchyStress = Cauchy stress cs(6) and jacobian = Consistent tangent dcs/dE - use numerics, only: defgradTolerance, & - iJacoStiffness - use debug, only: debug_level, & - debug_CPFEM, & - debug_levelBasic, & - debug_levelExtensive, & - debug_levelSelective, & - debug_e, & - debug_i, & - debug_stressMaxLocation, & - debug_stressMinLocation, & - debug_jacobianMaxLocation, & - debug_jacobianMinLocation, & - debug_stressMax, & - debug_stressMin, & - debug_jacobianMax, & - debug_jacobianMin - use FEsolving, only: outdatedFFN1, & - terminallyIll, & - cycleCounter, & - theInc, & - theTime, & - theDelta, & - FEsolving_execElem, & - FEsolving_execIP, & - restartWrite - use math, only: math_identity2nd, & - math_mul33x33, & - math_det33, & - math_transpose33, & - math_I3, & - math_Mandel3333to66, & - math_Mandel66to3333, & - math_Mandel33to6, & - math_Mandel6to33 - use mesh, only: mesh_FEasCP, & - mesh_NcpElems, & - mesh_maxNips, & - mesh_element - use material, only: homogenization_maxNgrains, & - microstructure_elemhomo, & - material_phase - use constitutive, only: constitutive_state0,constitutive_state - use crystallite, only: crystallite_partionedF,& - crystallite_F0, & - crystallite_Fp0, & - crystallite_Fp, & - crystallite_Lp0, & - crystallite_Lp, & - crystallite_dPdF0, & - crystallite_dPdF, & - crystallite_Tstar0_v, & - crystallite_Tstar_v - use homogenization, only: homogenization_sizeState, & - homogenization_state, & - homogenization_state0, & - materialpoint_F, & - materialpoint_F0, & - materialpoint_P, & - materialpoint_dPdF, & - materialpoint_results, & - materialpoint_sizeResults, & - materialpoint_Temperature, & - materialpoint_stressAndItsTangent, & - materialpoint_postResults - use IO, only: IO_write_jobRealFile, & - IO_warning +subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature, dt, elFE, ip, cauchyStress, jacobian) + use numerics, only: & + defgradTolerance, & + iJacoStiffness + use debug, only: & + debug_level, & + debug_CPFEM, & + debug_levelBasic, & + debug_levelExtensive, & + debug_levelSelective, & + debug_e, & + debug_i, & + debug_stressMaxLocation, & + debug_stressMinLocation, & + debug_jacobianMaxLocation, & + debug_jacobianMinLocation, & + debug_stressMax, & + debug_stressMin, & + debug_jacobianMax, & + debug_jacobianMin + use FEsolving, only: & + outdatedFFN1, & + terminallyIll, & + cycleCounter, & + theInc, & + theTime, & + theDelta, & + FEsolving_execElem, & + FEsolving_execIP, & + restartWrite + use math, only: & + math_identity2nd, & + math_mul33x33, & + math_det33, & + math_transpose33, & + math_I3, & + math_Mandel3333to66, & + math_Mandel66to3333, & + math_Mandel33to6, & + math_Mandel6to33 + use mesh, only: & + mesh_FEasCP, & + mesh_NcpElems, & + mesh_maxNips, & + mesh_element + use material, only: & + homogenization_maxNgrains, & + microstructure_elemhomo, & + material_phase + use constitutive, only: & + constitutive_state0,constitutive_state + use crystallite, only: & + crystallite_partionedF,& + crystallite_F0, & + crystallite_Fp0, & + crystallite_Fp, & + crystallite_Lp0, & + crystallite_Lp, & + crystallite_dPdF0, & + crystallite_dPdF, & + crystallite_Tstar0_v, & + crystallite_Tstar_v, & + crystallite_temperature + use homogenization, only: & + homogenization_sizeState, & + homogenization_state, & + homogenization_state0, & + materialpoint_F, & + materialpoint_F0, & + materialpoint_P, & + materialpoint_dPdF, & + materialpoint_results, & + materialpoint_sizeResults, & + materialpoint_stressAndItsTangent, & + materialpoint_postResults + use IO, only: & + IO_write_jobRealFile, & + IO_warning use DAMASK_interface - + implicit none - integer(pInt), intent(in) :: element, & !< FE element number - IP !< FE integration point number - real(pReal), intent(inout) :: Temperature !< temperature + integer(pInt), intent(in) :: elFE, & !< FE element number + ip !< integration point number + real(pReal), intent(in) :: temperature !< temperature real(pReal), intent(in) :: dt !< time increment real(pReal), dimension (3,3), intent(in) :: ffn, & !< deformation gradient for t=t0 ffn1 !< deformation gradient for t=t1 @@ -339,7 +347,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, el logical, intent(in) :: parallelExecution !< flag indicating parallel computation of requested IPs real(pReal), dimension(6), intent(out), optional :: cauchyStress !< stress vector in Mandel notation - real(pReal), dimension(6,6), intent(out), optional :: jacobian !< jacobian in Mandel notation + real(pReal), dimension(6,6), intent(out), optional :: jacobian !< jacobian in Mandel notation (Consistent tangent dcs/dE) real(pReal) J_inverse, & ! inverse of Jacobian rnd @@ -348,19 +356,20 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, el real(pReal), dimension (3,3,3,3) :: H_sym, & H, & jacobian3333 ! jacobian in Matrix notation - integer(pInt) cp_en, & ! crystal plasticity element number + integer(pInt) elCP, & ! crystal plasticity element number i, j, k, l, m, n logical updateJaco ! flag indicating if JAcobian has to be updated - cp_en = mesh_FEasCP('elem',element) + elCP = mesh_FEasCP('elem',elFE) +!if elCP = 0_pInt return ToDo: needed? if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt & - .and. cp_en == debug_e .and. IP == debug_i) then + .and. elCP == debug_e .and. ip == debug_i) then !$OMP CRITICAL (write2out) write(6,'(/,a)') '#############################################' - write(6,'(a1,a22,1x,i8,a13)') '#','element', cp_en, '#' - write(6,'(a1,a22,1x,i8,a13)') '#','IP', IP, '#' + write(6,'(a1,a22,1x,i8,a13)') '#','element', elCP, '#' + write(6,'(a1,a22,1x,i8,a13)') '#','ip', ip, '#' write(6,'(a1,a22,1x,f15.7,a6)') '#','theTime', theTime, '#' write(6,'(a1,a22,1x,f15.7,a6)') '#','theDelta', theDelta, '#' write(6,'(a1,a22,1x,i8,a13)') '#','theInc', theInc, '#' @@ -398,9 +407,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, el !$OMP CRITICAL (write2out) write(6,'(a)') '<< CPFEM >> aging states' if (debug_e <= mesh_NcpElems .and. debug_i <= mesh_maxNips) then - write(6,'(a,1x,i8,1x,i2,1x,i4,/,(12x,6(e20.8,1x)))') '<< CPFEM >> aged state of el ip grain',& + write(6,'(a,1x,i8,1x,i2,1x,i4,/,(12x,6(e20.8,1x)),/)') '<< CPFEM >> aged state of elFE ip grain',& debug_e, debug_i, 1, constitutive_state(1,debug_i,debug_e)%p - write(6,*) endif !$OMP END CRITICAL (write2out) endif @@ -479,18 +487,18 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, el !* In case that no parallel execution is required, there is no need to collect FEM input if (.not. parallelExecution) then - materialpoint_Temperature(IP,cp_en) = Temperature - materialpoint_F0(1:3,1:3,IP,cp_en) = ffn - materialpoint_F(1:3,1:3,IP,cp_en) = ffn1 + crystallite_temperature(ip,elCP) = temperature + materialpoint_F0(1:3,1:3,ip,elCP) = ffn + materialpoint_F(1:3,1:3,ip,elCP) = ffn1 elseif (iand(mode, CPFEM_COLLECT) /= 0_pInt) then call random_number(rnd) if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal - materialpoint_Temperature(IP,cp_en) = Temperature - materialpoint_F0(1:3,1:3,IP,cp_en) = ffn - materialpoint_F(1:3,1:3,IP,cp_en) = ffn1 - CPFEM_cs(1:6,IP,cp_en) = rnd * CPFEM_odd_stress - CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian * math_identity2nd(6) + crystallite_temperature(ip,elCP) = temperature + materialpoint_F0(1:3,1:3,ip,elCP) = ffn + materialpoint_F(1:3,1:3,ip,elCP) = ffn1 + CPFEM_cs(1:6,ip,elCP) = rnd * CPFEM_odd_stress + CPFEM_dcsde(1:6,1:6,ip,elCP) = CPFEM_odd_jacobian * math_identity2nd(6) CPFEM_calc_done = .false. endif @@ -502,13 +510,13 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, el !*** deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one - if (terminallyIll .or. outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(1:3,1:3,IP,cp_en)) > defgradTolerance)) then + if (terminallyIll .or. outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(1:3,1:3,ip,elCP)) > defgradTolerance)) then ! if (.not. terminallyIll .and. .not. outdatedFFN1) then - if (any(abs(ffn1 - materialpoint_F(1:3,1:3,IP,cp_en)) > defgradTolerance)) then + if (any(abs(ffn1 - materialpoint_F(1:3,1:3,ip,elCP)) > defgradTolerance)) then if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then !$OMP CRITICAL (write2out) - write(6,'(a,1x,i8,1x,i2)') '<< CPFEM >> OUTDATED at el ip',cp_en,IP - write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 old:',math_transpose33(materialpoint_F(1:3,1:3,IP,cp_en)) + write(6,'(a,1x,i8,1x,i2)') '<< CPFEM >> OUTDATED at elFE ip',elCP,ip + write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 old:',math_transpose33(materialpoint_F(1:3,1:3,ip,elCP)) write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 now:',math_transpose33(ffn1) !$OMP END CRITICAL (write2out) endif @@ -516,8 +524,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, el endif call random_number(rnd) if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal - CPFEM_cs(1:6,IP,cp_en) = rnd*CPFEM_odd_stress - CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian*math_identity2nd(6) + CPFEM_cs(1:6,ip,elCP) = rnd*CPFEM_odd_stress + CPFEM_dcsde(1:6,1:6,ip,elCP) = CPFEM_odd_jacobian*math_identity2nd(6) !*** deformation gradient is not outdated @@ -525,18 +533,18 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, el else updateJaco = mod(cycleCounter,iJacoStiffness) == 0 - !* no parallel computation, so we use just one single element and IP for computation + !* no parallel computation, so we use just one single elFE and ip for computation if (.not. parallelExecution) then - FEsolving_execElem(1) = cp_en - FEsolving_execElem(2) = cp_en - if (.not. microstructure_elemhomo(mesh_element(4,cp_en)) .or. & ! calculate unless homogeneous - (microstructure_elemhomo(mesh_element(4,cp_en)) .and. IP == 1_pInt)) then ! and then only first IP - FEsolving_execIP(1,cp_en) = IP - FEsolving_execIP(2,cp_en) = IP + FEsolving_execElem(1) = elCP + FEsolving_execElem(2) = elCP + if (.not. microstructure_elemhomo(mesh_element(4,elCP)) .or. & ! calculate unless homogeneous + (microstructure_elemhomo(mesh_element(4,elCP)) .and. ip == 1_pInt)) then ! and then only first ip + FEsolving_execIP(1,elCP) = ip + FEsolving_execIP(2,elCP) = ip if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then !$OMP CRITICAL (write2out) - write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for el ip ',cp_en,IP + write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for elFE ip ',elCP,ip !$OMP END CRITICAL (write2out) endif call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent @@ -561,98 +569,81 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, el if ( terminallyIll ) then call random_number(rnd) if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal - CPFEM_cs(1:6,IP,cp_en) = rnd * CPFEM_odd_stress - CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian * math_identity2nd(6) + CPFEM_cs(1:6,ip,elCP) = rnd * CPFEM_odd_stress + CPFEM_dcsde(1:6,1:6,ip,elCP) = CPFEM_odd_jacobian * math_identity2nd(6) else - if (microstructure_elemhomo(mesh_element(4,cp_en)) .and. IP > 1_pInt) then ! me homogenous? --> copy from first IP - materialpoint_P(1:3,1:3,IP,cp_en) = materialpoint_P(1:3,1:3,1,cp_en) - materialpoint_F(1:3,1:3,IP,cp_en) = materialpoint_F(1:3,1:3,1,cp_en) - materialpoint_dPdF(1:3,1:3,1:3,1:3,IP,cp_en) = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,cp_en) - materialpoint_results(1:materialpoint_sizeResults,IP,cp_en) = materialpoint_results(1:materialpoint_sizeResults,1,cp_en) + if (microstructure_elemhomo(mesh_element(4,elCP)) .and. ip > 1_pInt) then ! me homogenous? --> copy from first ip + materialpoint_P(1:3,1:3,ip,elCP) = materialpoint_P(1:3,1:3,1,elCP) + materialpoint_F(1:3,1:3,ip,elCP) = materialpoint_F(1:3,1:3,1,elCP) + materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,elCP) = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,elCP) + materialpoint_results(1:materialpoint_sizeResults,ip,elCP) = materialpoint_results(1:materialpoint_sizeResults,1,elCP) endif ! translate from P to CS - Kirchhoff = math_mul33x33(materialpoint_P(1:3,1:3,IP,cp_en), math_transpose33(materialpoint_F(1:3,1:3,IP,cp_en))) - J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,IP,cp_en)) - CPFEM_cs(1:6,IP,cp_en) = math_Mandel33to6(J_inverse * Kirchhoff) + Kirchhoff = math_mul33x33(materialpoint_P(1:3,1:3,ip,elCP), math_transpose33(materialpoint_F(1:3,1:3,ip,elCP))) + J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,ip,elCP)) + CPFEM_cs(1:6,ip,elCP) = math_Mandel33to6(J_inverse * Kirchhoff) ! translate from dP/dF to dCS/dE H = 0.0_pReal do i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3 H(i,j,k,l) = H(i,j,k,l) + & - materialpoint_F(j,m,IP,cp_en) * & - materialpoint_F(l,n,IP,cp_en) * & - materialpoint_dPdF(i,m,k,n,IP,cp_en) - & - math_I3(j,l) * materialpoint_F(i,m,IP,cp_en) * materialpoint_P(k,m,IP,cp_en) + & + materialpoint_F(j,m,ip,elCP) * & + materialpoint_F(l,n,ip,elCP) * & + materialpoint_dPdF(i,m,k,n,ip,elCP) - & + math_I3(j,l) * materialpoint_F(i,m,ip,elCP) * materialpoint_P(k,m,ip,elCP) + & 0.5_pReal * (math_I3(i,k) * Kirchhoff(j,l) + math_I3(j,l) * Kirchhoff(i,k) + & math_I3(i,l) * Kirchhoff(j,k) + math_I3(j,k) * Kirchhoff(i,l)) enddo; enddo; enddo; enddo; enddo; enddo - do i=1,3; do j=1,3; do k=1,3; do l=1,3 !< @ToDo use forall + forall(i=1:3, j=1:3,k=1:3,l=1:3) & H_sym(i,j,k,l) = 0.25_pReal * (H(i,j,k,l) + H(j,i,k,l) + H(i,j,l,k) + H(j,i,l,k)) - enddo; enddo; enddo; enddo - CPFEM_dcsde(1:6,1:6,IP,cp_en) = math_Mandel3333to66(J_inverse * H_sym) + CPFEM_dcsde(1:6,1:6,ip,elCP) = math_Mandel3333to66(J_inverse * H_sym) endif endif - !* remember extreme values of stress and jacobian - - cauchyStress33 = math_Mandel6to33(CPFEM_cs(1:6,IP,cp_en)) + cauchyStress33 = math_Mandel6to33(CPFEM_cs(1:6,ip,elCP)) if (maxval(cauchyStress33) > debug_stressMax) then - debug_stressMaxLocation = [cp_en, IP] + debug_stressMaxLocation = [elCP, ip] debug_stressMax = maxval(cauchyStress33) endif if (minval(cauchyStress33) < debug_stressMin) then - debug_stressMinLocation = [cp_en, IP] + debug_stressMinLocation = [elCP, ip] debug_stressMin = minval(cauchyStress33) endif - jacobian3333 = math_Mandel66to3333(CPFEM_dcsdE(1:6,1:6,IP,cp_en)) + jacobian3333 = math_Mandel66to3333(CPFEM_dcsdE(1:6,1:6,ip,elCP)) if (maxval(jacobian3333) > debug_jacobianMax) then - debug_jacobianMaxLocation = [cp_en, IP] + debug_jacobianMaxLocation = [elCP, ip] debug_jacobianMax = maxval(jacobian3333) endif if (minval(jacobian3333) < debug_jacobianMin) then - debug_jacobianMinLocation = [cp_en, IP] + debug_jacobianMinLocation = [elCP, ip] debug_jacobianMin = minval(jacobian3333) endif !* report stress and stiffness - if ((iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) & - .and. ((debug_e == cp_en .and. debug_i == IP) & + .and. ((debug_e == elCP .and. debug_i == ip) & .or. .not. iand(debug_level(debug_CPFEM), debug_levelSelective) /= 0_pInt)) then !$OMP CRITICAL (write2out) - write(6,'(a,i8,1x,i2,/,12x,6(f10.3,1x)/)') '<< CPFEM >> stress/MPa at el ip ', & - cp_en, IP, CPFEM_cs(1:6,IP,cp_en)/1.0e6_pReal - write(6,'(a,i8,1x,i2,/,6(12x,6(f10.3,1x)/))') '<< CPFEM >> Jacobian/GPa at el ip ', & - cp_en, IP, transpose(CPFEM_dcsdE(1:6,1:6,IP,cp_en))/1.0e9_pReal + write(6,'(a,i8,1x,i2,/,12x,6(f10.3,1x)/)') '<< CPFEM >> stress/MPa at elFE ip ', & + elCP, ip, CPFEM_cs(1:6,ip,elCP)/1.0e6_pReal + write(6,'(a,i8,1x,i2,/,6(12x,6(f10.3,1x)/))') '<< CPFEM >> Jacobian/GPa at elFE ip ', & + elCP, ip, transpose(CPFEM_dcsdE(1:6,1:6,ip,elCP))/1.0e9_pReal flush(6) !$OMP END CRITICAL (write2out) endif endif - - - - !*** homogenized result except for potentially non-isothermal starting condition - - if (theTime > 0.0_pReal) then - Temperature = materialpoint_Temperature(IP,cp_en) - endif - - !*** warn if stiffness close to zero + if (all(abs(CPFEM_dcsdE(1:6,1:6,ip,elCP)) < 1e-10_pReal)) call IO_warning(601,elCP,ip) - if (all(abs(CPFEM_dcsdE(1:6,1:6,IP,cp_en)) < 1e-10_pReal)) then - call IO_warning(601,cp_en,IP) - endif !*** copy to output if required (FEM solver) - - if(present(cauchyStress)) cauchyStress = CPFEM_cs(1:6,IP,cp_en) - if(present(jacobian)) jacobian = CPFEM_dcsdE(1:6,1:6,IP,cp_en) + if(present(cauchyStress)) cauchyStress = CPFEM_cs(1:6,ip,elCP) + if(present(jacobian)) jacobian = CPFEM_dcsdE(1:6,1:6,ip,elCP) end subroutine CPFEM_general diff --git a/code/DAMASK_spectral_driver.f90 b/code/DAMASK_spectral_driver.f90 index ab68d5843..e6bfa3196 100644 --- a/code/DAMASK_spectral_driver.f90 +++ b/code/DAMASK_spectral_driver.f90 @@ -151,7 +151,7 @@ program DAMASK_spectral_Driver external :: quit !-------------------------------------------------------------------------------------------------- ! init DAMASK (all modules) - call CPFEM_initAll(temperature = 300.0_pReal, element = 1_pInt, IP= 1_pInt) + call CPFEM_initAll(temperature = 300.0_pReal, el = 1_pInt, ip = 1_pInt) write(6,'(/,a)') ' <<<+- DAMASK_spectral_driver init -+>>>' write(6,'(a)') ' $Id$' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() diff --git a/code/DAMASK_spectral_utilities.f90 b/code/DAMASK_spectral_utilities.f90 index 653f27b6a..ab9fab231 100644 --- a/code/DAMASK_spectral_utilities.f90 +++ b/code/DAMASK_spectral_utilities.f90 @@ -819,10 +819,11 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& CPFEM_AGERESULTS, & CPFEM_BACKUPJACOBIAN, & CPFEM_RESTOREJACOBIAN + use crystallite, only: & + crystallite_temperature use homogenization, only: & materialpoint_F0, & materialpoint_F, & - materialpoint_Temperature, & materialpoint_P, & materialpoint_dPdF @@ -864,7 +865,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid)]) materialpoint_F = reshape(F, [3,3,1,product(grid)]) - materialpoint_Temperature = temperature + crystallite_temperature = temperature call debug_reset() diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 96c883932..c0d72feab 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -66,7 +66,7 @@ module constitutive constitutive_TandItsTangent, & constitutive_collectDotState, & constitutive_collectDeltaState, & - constitutive_dotTemperature, & + constitutive_heat, & constitutive_postResults private :: & @@ -463,12 +463,24 @@ pure function constitutive_homogenizedC(ipc,ip,el) use material, only: & phase_plasticity, & material_phase - use constitutive_none - use constitutive_j2 - use constitutive_phenopowerlaw - use constitutive_titanmod - use constitutive_dislotwin - use constitutive_nonlocal + use constitutive_none, only: & + CONSTITUTIVE_NONE_label, & + constitutive_none_homogenizedC + use constitutive_j2, only: & + CONSTITUTIVE_J2_label, & + constitutive_j2_homogenizedC + use constitutive_phenopowerlaw, only: & + CONSTITUTIVE_PHENOPOWERLAW_label, & + constitutive_phenopowerlaw_homogenizedC + use constitutive_titanmod, only: & + CONSTITUTIVE_TITANMOD_label, & + constitutive_titanmod_homogenizedC + use constitutive_dislotwin, only: & + CONSTITUTIVE_DISLOTWIN_label, & + constitutive_dislotwin_homogenizedC + use constitutive_nonlocal, only: & + CONSTITUTIVE_NONLOCAL_label, & + constitutive_nonlocal_homogenizedC implicit none real(pReal), dimension(6,6) :: constitutive_homogenizedC @@ -480,23 +492,23 @@ pure function constitutive_homogenizedC(ipc,ip,el) select case (phase_plasticity(material_phase(ipc,ip,el))) case (constitutive_none_label) - constitutive_homogenizedC = constitutive_none_homogenizedC(constitutive_state,ipc,ip,el) + constitutive_homogenizedC = constitutive_none_homogenizedC(ipc,ip,el) case (constitutive_j2_label) - constitutive_homogenizedC = constitutive_j2_homogenizedC(constitutive_state,ipc,ip,el) + constitutive_homogenizedC = constitutive_j2_homogenizedC(ipc,ip,el) case (constitutive_phenopowerlaw_label) - constitutive_homogenizedC = constitutive_phenopowerlaw_homogenizedC(constitutive_state,ipc,ip,el) + constitutive_homogenizedC = constitutive_phenopowerlaw_homogenizedC(ipc,ip,el) + case (constitutive_nonlocal_label) + constitutive_homogenizedC = constitutive_nonlocal_homogenizedC(ipc,ip,el) + case (constitutive_titanmod_label) constitutive_homogenizedC = constitutive_titanmod_homogenizedC(constitutive_state,ipc,ip,el) case (constitutive_dislotwin_label) constitutive_homogenizedC = constitutive_dislotwin_homogenizedC(constitutive_state,ipc,ip,el) - case (constitutive_nonlocal_label) - constitutive_homogenizedC = constitutive_nonlocal_homogenizedC(constitutive_state,ipc,ip,el) - end select end function constitutive_homogenizedC @@ -505,7 +517,7 @@ end function constitutive_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different constitutive models !-------------------------------------------------------------------------------------------------- -subroutine constitutive_microstructure(Temperature, Fe, Fp, ipc, ip, el) +subroutine constitutive_microstructure(temperature, Fe, Fp, ipc, ip, el) use material, only: & phase_plasticity, & material_phase @@ -525,7 +537,7 @@ subroutine constitutive_microstructure(Temperature, Fe, Fp, ipc, ip, el) ip, & !< integration point number el !< element number real(pReal), intent(in) :: & - Temperature + temperature real(pReal), intent(in), dimension(3,3) :: & Fe, & !< elastic deformation gradient Fp !< plastic deformation gradient @@ -533,13 +545,13 @@ subroutine constitutive_microstructure(Temperature, Fe, Fp, ipc, ip, el) select case (phase_plasticity(material_phase(ipc,ip,el))) case (constitutive_titanmod_label) - call constitutive_titanmod_microstructure(Temperature,constitutive_state,ipc,ip,el) + call constitutive_titanmod_microstructure(temperature,constitutive_state,ipc,ip,el) case (constitutive_dislotwin_label) - call constitutive_dislotwin_microstructure(Temperature,constitutive_state,ipc,ip,el) + call constitutive_dislotwin_microstructure(temperature,constitutive_state,ipc,ip,el) case (constitutive_nonlocal_label) - call constitutive_nonlocal_microstructure(constitutive_state, Temperature, Fe, Fp, ipc,ip,el) + call constitutive_nonlocal_microstructure(constitutive_state,Fe,Fp,ipc,ip,el) end select @@ -549,13 +561,14 @@ end subroutine constitutive_microstructure !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, ipc, ip, el) +subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, ipc, ip, el) + use math, only: & + math_identity2nd use material, only: & phase_plasticity, & material_phase use constitutive_none, only: & - constitutive_none_label, & - constitutive_none_LpAndItsTangent + constitutive_none_label use constitutive_j2, only: & constitutive_j2_label, & constitutive_j2_LpAndItsTangent @@ -589,22 +602,23 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, ip select case (phase_plasticity(material_phase(ipc,ip,el))) case (constitutive_none_label) - call constitutive_none_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el) + Lp = 0.0_pReal + dLp_dTstar = math_identity2nd(9) case (constitutive_j2_label) - call constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el) + call constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,constitutive_state,ipc,ip,el) case (constitutive_phenopowerlaw_label) - call constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el) + call constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,constitutive_state,ipc,ip,el) case (constitutive_titanmod_label) - call constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el) + call constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,temperature,constitutive_state,ipc,ip,el) case (constitutive_dislotwin_label) - call constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el) + call constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,temperature,constitutive_state,ipc,ip,el) case (constitutive_nonlocal_label) - call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, constitutive_state(ipc,ip,el), ipc,ip,el) + call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, constitutive_state(ipc,ip,el), ipc,ip,el) end select @@ -744,10 +758,10 @@ subroutine constitutive_collectDotState(Tstar_v, Fe, Fp, Temperature, subdt, sub constitutive_dotState(ipc,ip,el)%p = 0.0_pReal !ToDo: needed or will it remain zero anyway? case (constitutive_j2_label) - constitutive_dotState(ipc,ip,el)%p = constitutive_j2_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) + constitutive_dotState(ipc,ip,el)%p = constitutive_j2_dotState(Tstar_v,constitutive_state,ipc,ip,el) case (constitutive_phenopowerlaw_label) - constitutive_dotState(ipc,ip,el)%p = constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) + constitutive_dotState(ipc,ip,el)%p = constitutive_phenopowerlaw_dotState(Tstar_v,constitutive_state,ipc,ip,el) case (constitutive_titanmod_label) constitutive_dotState(ipc,ip,el)%p = constitutive_titanmod_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) @@ -778,7 +792,7 @@ end subroutine constitutive_collectDotState !> @brief contains the constitutive equation for calculating the incremental change of !> microstructure based on the current stress and state !-------------------------------------------------------------------------------------------------- -subroutine constitutive_collectDeltaState(Tstar_v, Temperature, ipc, ip, el) +subroutine constitutive_collectDeltaState(Tstar_v, ipc, ip, el) use prec, only: & pLongInt use debug, only: & @@ -799,8 +813,6 @@ subroutine constitutive_collectDeltaState(Tstar_v, Temperature, ipc, ip, el) ipc, & !< grain number ip, & !< integration point number el !< element number - real(pReal), intent(in) :: & - Temperature real(pReal), intent(in), dimension(6) :: & Tstar_v !< 2nd Piola-Kirchhoff stress integer(pLongInt) :: & @@ -814,7 +826,7 @@ subroutine constitutive_collectDeltaState(Tstar_v, Temperature, ipc, ip, el) select case (phase_plasticity(material_phase(ipc,ip,el))) case (constitutive_nonlocal_label) - call constitutive_nonlocal_deltaState(constitutive_deltaState(ipc,ip,el),constitutive_state, Tstar_v,Temperature,ipc,ip,el) + call constitutive_nonlocal_deltaState(constitutive_deltaState(ipc,ip,el),constitutive_state, Tstar_v,ipc,ip,el) case default constitutive_deltaState(ipc,ip,el)%p = 0.0_pReal !ToDo: needed or will it remain zero anyway? @@ -837,7 +849,7 @@ end subroutine constitutive_collectDeltaState !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -real(pReal) function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el) +real(pReal) function constitutive_heat(Tstar_v,Temperature,ipc,ip,el) implicit none integer(pInt), intent(in) :: & ipc, & !< grain number @@ -847,14 +859,16 @@ real(pReal) function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el) Temperature real(pReal), intent(in), dimension(6) :: & Tstar_v !< 2nd Piola-Kirchhoff stress - constitutive_dotTemperature = 0.0_pReal -end function constitutive_dotTemperature + + constitutive_heat = 0.0_pReal + +end function constitutive_heat !-------------------------------------------------------------------------------------------------- !> @brief returns array of constitutive results !-------------------------------------------------------------------------------------------------- -function constitutive_postResults(Tstar_v, Fe, Temperature, dt, ipc, ip, el) +function constitutive_postResults(Tstar_v, Fe, temperature, ipc, ip, el) use mesh, only: & mesh_NcpElems, & mesh_maxNips @@ -888,8 +902,7 @@ function constitutive_postResults(Tstar_v, Fe, Temperature, dt, ipc, ip, el) real(pReal), dimension(constitutive_sizePostResults(ipc,ip,el)) :: & constitutive_postResults real(pReal), intent(in) :: & - Temperature, & - dt !< timestep + temperature real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & Fe !< elastic deformation gradient real(pReal), intent(in), dimension(6) :: & @@ -902,20 +915,20 @@ function constitutive_postResults(Tstar_v, Fe, Temperature, dt, ipc, ip, el) case (constitutive_none_label) constitutive_postResults = 0.0_pReal + case (constitutive_titanmod_label) + constitutive_postResults = constitutive_titanmod_postResults(constitutive_state,ipc,ip,el) + case (constitutive_j2_label) - constitutive_postResults = constitutive_j2_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el) + constitutive_postResults = constitutive_j2_postResults(Tstar_v,constitutive_state,ipc,ip,el) case (constitutive_phenopowerlaw_label) - constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el) + constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,constitutive_state,ipc,ip,el) - case (constitutive_titanmod_label) - constitutive_postResults = constitutive_titanmod_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el) - case (constitutive_dislotwin_label) - constitutive_postResults = constitutive_dislotwin_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el) + constitutive_postResults = constitutive_dislotwin_postResults(Tstar_v,Temperature,constitutive_state,ipc,ip,el) case (constitutive_nonlocal_label) - constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, Fe, Temperature, dt, constitutive_state, & + constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, Fe, constitutive_state, & constitutive_dotstate(ipc,ip,el), ipc, ip, el) end select diff --git a/code/constitutive_dislotwin.f90 b/code/constitutive_dislotwin.f90 index 2672001b5..0d6a129bb 100644 --- a/code/constitutive_dislotwin.f90 +++ b/code/constitutive_dislotwin.f90 @@ -934,10 +934,10 @@ pure function constitutive_dislotwin_homogenizedC(state,ipc,ip,el) sumf = sum(state(ipc,ip,el)%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt))) ! safe for nt == 0 !* Homogenized elasticity matrix - constitutive_dislotwin_homogenizedC = (1.0_pReal-sumf)*constitutive_dislotwin_Cslip_66(:,:,matID) + constitutive_dislotwin_homogenizedC = (1.0_pReal-sumf)*constitutive_dislotwin_Cslip_66(1:6,1:6,matID) do i=1_pInt,nt constitutive_dislotwin_homogenizedC = & - constitutive_dislotwin_homogenizedC + state(ipc,ip,el)%p(3_pInt*ns+i)*constitutive_dislotwin_Ctwin_66(:,:,i,matID) + constitutive_dislotwin_homogenizedC + state(ipc,ip,el)%p(3_pInt*ns+i)*constitutive_dislotwin_Ctwin_66(1:6,1:6,i,matID) enddo end function constitutive_dislotwin_homogenizedC @@ -945,7 +945,7 @@ pure function constitutive_dislotwin_homogenizedC(state,ipc,ip,el) !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -subroutine constitutive_dislotwin_microstructure(Temperature,state,ipc,ip,el) +subroutine constitutive_dislotwin_microstructure(temperature,state,ipc,ip,el) use prec, only: & p_vec use math, only: & @@ -1529,7 +1529,7 @@ end function constitutive_dislotwin_dotState !-------------------------------------------------------------------------------------------------- !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- -function constitutive_dislotwin_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el) +function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el) use prec, only: & p_vec use math, only: & @@ -1558,8 +1558,7 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,dt,state,ipc,ip, real(pReal), dimension(6), intent(in) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation real(pReal), intent(in) :: & - temperature, & !< temperature at integration point - dt + temperature !< temperature at integration point integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point diff --git a/code/constitutive_j2.f90 b/code/constitutive_j2.f90 index 827211d9d..0cd0e7d43 100644 --- a/code/constitutive_j2.f90 +++ b/code/constitutive_j2.f90 @@ -340,9 +340,7 @@ end function constitutive_j2_aTolState !-------------------------------------------------------------------------------------------------- !> @brief returns the homogenized elasticity matrix !-------------------------------------------------------------------------------------------------- -pure function constitutive_j2_homogenizedC(state,ipc,ip,el) - use prec, only: & - p_vec +pure function constitutive_j2_homogenizedC(ipc,ip,el) use mesh, only: & mesh_NcpElems, & mesh_maxNips @@ -358,8 +356,6 @@ pure function constitutive_j2_homogenizedC(state,ipc,ip,el) ipc, & !< component-ID of integration point ip, & !< integration point el !< element - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & - state !< microstructure state constitutive_j2_homogenizedC = constitutive_j2_Cslip_66(1:6,1:6,& phase_plasticityInstance(material_phase(ipc,ip,el))) @@ -370,8 +366,7 @@ end function constitutive_j2_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -pure subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,& - temperature,state,ipc,ip,el) +pure subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,state,ipc,ip,el) use prec, only: & p_vec use math, only: & @@ -396,8 +391,6 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,& real(pReal), dimension(6), intent(in) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal), intent(in) :: & - temperature !< temperature at IP integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -449,7 +442,7 @@ end subroutine constitutive_j2_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -pure function constitutive_j2_dotState(Tstar_v,temperature,state,ipc,ip,el) +pure function constitutive_j2_dotState(Tstar_v,state,ipc,ip,el) use prec, only: & p_vec use math, only: & @@ -467,8 +460,6 @@ pure function constitutive_j2_dotState(Tstar_v,temperature,state,ipc,ip,el) constitutive_j2_dotState real(pReal), dimension(6), intent(in):: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal), intent(in) :: & - temperature !< temperature at integration point integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -533,7 +524,7 @@ end function constitutive_j2_dotState !-------------------------------------------------------------------------------------------------- !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- -pure function constitutive_j2_postResults(Tstar_v,temperature,dt,state,ipc,ip,el) +pure function constitutive_j2_postResults(Tstar_v,state,ipc,ip,el) use prec, only: & p_vec use math, only: & @@ -550,9 +541,6 @@ pure function constitutive_j2_postResults(Tstar_v,temperature,dt,state,ipc,ip,el implicit none real(pReal), dimension(6), intent(in) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal), intent(in) :: & - temperature, & !< temperature at integration point - dt integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point diff --git a/code/constitutive_none.f90 b/code/constitutive_none.f90 index 17416023f..afb07f948 100644 --- a/code/constitutive_none.f90 +++ b/code/constitutive_none.f90 @@ -49,8 +49,7 @@ module constitutive_none public :: & constitutive_none_init, & - constitutive_none_homogenizedC, & - constitutive_none_LpAndItsTangent + constitutive_none_homogenizedC contains @@ -185,7 +184,7 @@ end subroutine constitutive_none_init !-------------------------------------------------------------------------------------------------- !> @brief returns the homogenized elasticity matrix !-------------------------------------------------------------------------------------------------- -pure function constitutive_none_homogenizedC(state,ipc,ip,el) +pure function constitutive_none_homogenizedC(ipc,ip,el) use prec, only: & p_vec use mesh, only: & @@ -203,51 +202,10 @@ pure function constitutive_none_homogenizedC(state,ipc,ip,el) ipc, & !< component-ID of integration point ip, & !< integration point el !< element - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & - state !< microstructure state constitutive_none_homogenizedC = constitutive_none_Cslip_66(1:6,1:6,& phase_plasticityInstance(material_phase(ipc,ip,el))) end function constitutive_none_homogenizedC - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent -!> @details dummy function, returns 0.0 and Identity -!-------------------------------------------------------------------------------------------------- -pure subroutine constitutive_none_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_dev_v, & - temperature, state, ipc, ip, el) - use prec, only: & - p_vec - use math, only: & - math_identity2nd - use mesh, only: & - mesh_NcpElems, & - mesh_maxNips - use material, only: & - homogenization_maxNgrains, & - 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_dev_v !< deviatoric part of 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal), intent(in) :: & - temperature !< temperature at IP - integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & - state !< microstructure state - - Lp = 0.0_pReal ! set Lp to zero - dLp_dTstar99 = math_identity2nd(9) ! set dLp_dTstar to Identity - -end subroutine constitutive_none_LpAndItsTangent - end module constitutive_none diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index c086c399e..f031a80a6 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -45,30 +45,30 @@ character (len=*), parameter, public :: & CONSTITUTIVE_NONLOCAL_LABEL = 'nonlocal' character(len=22), dimension(11), parameter, private :: & -BASICSTATES = (/'rhoSglEdgePosMobile ', & - 'rhoSglEdgeNegMobile ', & - 'rhoSglScrewPosMobile ', & - 'rhoSglScrewNegMobile ', & - 'rhoSglEdgePosImmobile ', & - 'rhoSglEdgeNegImmobile ', & - 'rhoSglScrewPosImmobile', & - 'rhoSglScrewNegImmobile', & - 'rhoDipEdge ', & - 'rhoDipScrew ', & - 'accumulatedshear ' /) !< list of "basic" microstructural state variables that are independent from other state variables +BASICSTATES = ['rhoSglEdgePosMobile ', & + 'rhoSglEdgeNegMobile ', & + 'rhoSglScrewPosMobile ', & + 'rhoSglScrewNegMobile ', & + 'rhoSglEdgePosImmobile ', & + 'rhoSglEdgeNegImmobile ', & + 'rhoSglScrewPosImmobile', & + 'rhoSglScrewNegImmobile', & + 'rhoDipEdge ', & + 'rhoDipScrew ', & + 'accumulatedshear ' ] !< list of "basic" microstructural state variables that are independent from other state variables character(len=16), dimension(3), parameter, private :: & -DEPENDENTSTATES = (/'rhoForest ', & - 'tauThreshold ', & - 'tauBack ' /) !< list of microstructural state variables that depend on other state variables +DEPENDENTSTATES = ['rhoForest ', & + 'tauThreshold ', & + 'tauBack ' ] !< list of microstructural state variables that depend on other state variables character(len=20), dimension(6), parameter, private :: & -OTHERSTATES = (/'velocityEdgePos ', & - 'velocityEdgeNeg ', & - 'velocityScrewPos ', & - 'velocityScrewNeg ', & - 'maxDipoleHeightEdge ', & - 'maxDipoleHeightScrew' /) !< list of other dependent state variables that are not updated by microstructure +OTHERSTATES = ['velocityEdgePos ', & + 'velocityEdgeNeg ', & + 'velocityScrewPos ', & + 'velocityScrewNeg ', & + 'maxDipoleHeightEdge ', & + 'maxDipoleHeightScrew' ] !< list of other dependent state variables that are not updated by microstructure real(pReal), parameter, private :: & KB = 1.38e-23_pReal !< Physical parameter, Boltzmann constant in J/Kelvin @@ -1200,38 +1200,33 @@ end function constitutive_nonlocal_aTolState !-------------------------------------------------------------------------------------------------- !> @brief returns the homogenized elasticity matrix !-------------------------------------------------------------------------------------------------- -pure function constitutive_nonlocal_homogenizedC(state,g,ip,el) - -use mesh, only: mesh_NcpElems, & - mesh_maxNips -use material, only: homogenization_maxNgrains, & - material_phase, & - phase_plasticityInstance -implicit none - -!*** input variables -integer(pInt), intent(in) :: g, & ! current grain ID - ip, & ! current integration point - el ! current element -type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state ! microstructural state - -!*** output variables -real(pReal), dimension(6,6) :: constitutive_nonlocal_homogenizedC ! homogenized elasticity matrix - -!*** local variables -integer(pInt) :: matID ! current instance of this plasticity - -matID = phase_plasticityInstance(material_phase(g,ip,el)) - -constitutive_nonlocal_homogenizedC = Cslip66(1:6,1:6,matID) +pure function constitutive_nonlocal_homogenizedC(ipc,ip,el) + use mesh, only: & + mesh_NcpElems, & + mesh_maxNips + use material, only: & + homogenization_maxNgrains, & + material_phase, & + phase_plasticityInstance + implicit none + integer(pInt), intent(in) :: & + ipc, & ! current grain ID + ip, & ! current integration point + el ! current element + real(pReal), dimension(6,6) :: & + constitutive_nonlocal_homogenizedC + + constitutive_nonlocal_homogenizedC = & + Cslip66(1:6,1:6,phase_plasticityInstance(material_phase(ipc,ip,el))) + end function constitutive_nonlocal_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief calculates quantities characterizing the microstructure !-------------------------------------------------------------------------------------------------- -subroutine constitutive_nonlocal_microstructure(state, Temperature, Fe, Fp, gr, ip, el) +subroutine constitutive_nonlocal_microstructure(state, Fe, Fp, gr, ip, el) use IO, only: & IO_error @@ -1279,7 +1274,6 @@ implicit none integer(pInt), intent(in) :: gr, & ! current grain ID ip, & ! current integration point el ! current element -real(pReal), intent(in) :: Temperature ! temperature real(pReal), dimension(3,3), intent(in) :: & Fe, & ! elastic deformation gradient Fp ! elastic deformation gradient @@ -1547,7 +1541,7 @@ end subroutine constitutive_nonlocal_microstructure !> @brief calculates kinetics !-------------------------------------------------------------------------------------------------- subroutine constitutive_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, & - tauThreshold, c, Temperature, g, ip, el) + tauThreshold, c, Temperature, ipc, ip, el) use debug, only: debug_level, & debug_constitutive, & @@ -1563,18 +1557,18 @@ use material, only: material_phase, & implicit none !*** input variables -integer(pInt), intent(in) :: g, & !< current grain number +integer(pInt), intent(in) :: ipc, & !< current grain number ip, & !< current integration point el, & !< current element number c !< dislocation character (1:edge, 2:screw) real(pReal), intent(in) :: Temperature !< temperature -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))), & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))), & intent(in) :: tau, & !< resolved external shear stress (without non Schmid effects) tauNS, & !< resolved external shear stress (including non Schmid effects) tauThreshold !< threshold shear stress !*** output variables -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))), & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))), & intent(out) :: v, & !< velocity dv_dtau, & !< velocity derivative with respect to resolved shear stress (without non Schmid contributions) dv_dtauNS !< velocity derivative with respect to resolved shear stress (including non Schmid contributions) @@ -1606,7 +1600,7 @@ real(pReal) tauRel_P, & mobility !< dislocation mobility -instance = phase_plasticityInstance(material_phase(g,ip,el)) +instance = phase_plasticityInstance(material_phase(ipc,ip,el)) ns = totalNslip(instance) v = 0.0_pReal @@ -1689,10 +1683,10 @@ endif #ifndef _OPENMP if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & - .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)& + .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then write(6,*) - write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_kinetics at el ip g',el,ip,g + write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_kinetics at el ip ipc',el,ip,ipc write(6,*) write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tauThreshold / MPa', tauThreshold / 1e6_pReal write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tau / MPa', tau / 1e6_pReal @@ -1709,7 +1703,7 @@ end subroutine constitutive_nonlocal_kinetics !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -subroutine constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, state, g, ip, el) +subroutine constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, state, ipc, ip, el) use math, only: math_Plain3333to99, & math_mul6x6, & @@ -1733,7 +1727,7 @@ use mesh, only: mesh_ipVolume implicit none !*** input variables -integer(pInt), intent(in) :: g, & !< current grain number +integer(pInt), intent(in) :: ipc, & !< current grain number ip, & !< current integration point el !< current element number real(pReal), intent(in) :: Temperature !< temperature @@ -1758,14 +1752,14 @@ integer(pInt) matID, & !< current in s, & !< index of my current slip system sLattice !< index of my current slip system according to lattice order real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 !< derivative of Lp with respect to Tstar (3x3x3x3 matrix) -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),8) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: & rhoSgl !< single dislocation densities (including blocked) -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),4) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: & v, & !< velocity tauNS, & !< resolved shear stress including non Schmid and backstress terms dv_dtau, & !< velocity derivative with respect to the shear stress dv_dtauNS !< velocity derivative with respect to the shear stress -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & tau, & !< resolved shear stress including backstress terms gdotTotal, & !< shear rate tauBack, & !< back stress from dislocation gradients on same slip system @@ -1777,7 +1771,7 @@ real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,e Lp = 0.0_pReal dLp_dTstar3333 = 0.0_pReal -matID = phase_plasticityInstance(material_phase(g,ip,el)) +matID = phase_plasticityInstance(material_phase(ipc,ip,el)) structID = constitutive_nonlocal_structure(matID) ns = totalNslip(matID) @@ -1823,7 +1817,7 @@ tau = tau + tauBack ! edges call constitutive_nonlocal_kinetics(v(1:ns,1), dv_dtau(1:ns,1), dv_dtauNS(1:ns,1), & tau(1:ns), tauNS(1:ns,1), tauThreshold(1:ns), & - 1_pInt, Temperature, g, ip, el) + 1_pInt, Temperature, ipc, ip, el) v(1:ns,2) = v(1:ns,1) dv_dtau(1:ns,2) = dv_dtau(1:ns,1) dv_dtauNS(1:ns,2) = dv_dtauNS(1:ns,1) @@ -1839,7 +1833,7 @@ else do t = 3_pInt,4_pInt call constitutive_nonlocal_kinetics(v(1:ns,t), dv_dtau(1:ns,t), dv_dtauNS(1:ns,t), & tau(1:ns), tauNS(1:ns,t), tauThreshold(1:ns), & - 2_pInt , Temperature, g, ip, el) + 2_pInt , Temperature, ipc, ip, el) enddo endif @@ -1892,10 +1886,10 @@ dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) #ifndef _OPENMP if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & - .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)& + .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then write(6,*) - write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_LpandItsTangent at el ip g ',el,ip,g + write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_LpandItsTangent at el ip ipc ',el,ip,ipc write(6,*) write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> gdot total / 1e-3',gdotTotal*1e3_pReal write(6,'(a,/,3(12x,3(f12.7,1x),/))') '<< CONST >> Lp',transpose(Lp) @@ -1908,7 +1902,7 @@ end subroutine constitutive_nonlocal_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief (instantaneous) incremental change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine constitutive_nonlocal_deltaState(deltaState, state, Tstar_v, Temperature, g,ip,el) +subroutine constitutive_nonlocal_deltaState(deltaState, state, Tstar_v, ipc,ip,el) use debug, only: debug_level, & debug_constitutive, & @@ -1931,10 +1925,9 @@ use material, only: homogenization_maxNgrains, & implicit none !*** input variables -integer(pInt), intent(in) :: g, & ! current grain number +integer(pInt), intent(in) :: ipc, & ! current grain number ip, & ! current integration point el ! current element number -real(pReal), intent(in) :: Temperature ! temperature real(pReal), dimension(6), intent(in) :: Tstar_v ! current 2nd Piola-Kirchhoff stress in Mandel notation !*** input/output variables @@ -1952,18 +1945,18 @@ integer(pInt) matID, & ! current insta t, & ! type of dislocation s, & ! index of my current slip system sLattice ! index of my current slip system according to lattice order -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),10) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),10) :: & deltaRho, & ! density increment deltaRhoRemobilization, & ! density increment by remobilization deltaRhoDipole2SingleStress ! density increment by dipole dissociation (by stress change) -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),8) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: & rhoSgl ! current single dislocation densities (positive/negative screw and edge without dipoles) -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),4) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: & v ! dislocation glide velocity -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & tau, & ! current resolved shear stress tauBack ! current back stress from pileups on same slip system -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),2) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: & rhoDip, & ! current dipole dislocation densities (screw and edge dipoles) dLower, & ! minimum stable dipole distance for edges and screws dUpper, & ! current maximum stable dipole distance for edges and screws @@ -1973,15 +1966,15 @@ real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,e #ifndef _OPENMP if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt & - .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)& + .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then write(6,*) - write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_deltaState at el ip g ',el,ip,g + write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_deltaState at el ip ipc ',el,ip,ipc write(6,*) endif #endif -matID = phase_plasticityInstance(material_phase(g,ip,el)) +matID = phase_plasticityInstance(material_phase(ipc,ip,el)) structID = constitutive_nonlocal_structure(matID) ns = totalNslip(matID) @@ -1990,15 +1983,15 @@ ns = totalNslip(matID) forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) - rhoSgl(s,t) = max(state(g,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal) ! ensure positive single mobile densities - rhoSgl(s,t+4_pInt) = state(g,ip,el)%p(iRhoB(s,t,matID)) - v(s,t) = state(g,ip,el)%p(iV(s,t,matID)) + rhoSgl(s,t) = max(state(ipc,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal) ! ensure positive single mobile densities + rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,matID)) + v(s,t) = state(ipc,ip,el)%p(iV(s,t,matID)) endforall forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) - rhoDip(s,c) = max(state(g,ip,el)%p(iRhoD(s,c,matID)), 0.0_pReal) ! ensure positive dipole densities - dUpperOld(s,c) = state(g,ip,el)%p(iD(s,c,matID)) + rhoDip(s,c) = max(state(ipc,ip,el)%p(iRhoD(s,c,matID)), 0.0_pReal) ! ensure positive dipole densities + dUpperOld(s,c) = state(ipc,ip,el)%p(iD(s,c,matID)) endforall -tauBack = state(g,ip,el)%p(iTauB(1:ns,matID)) +tauBack = state(ipc,ip,el)%p(iTauB(1:ns,matID)) where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(matID) & .or. abs(rhoSgl) < significantRho(matID)) & @@ -2063,7 +2056,7 @@ forall (t=1_pInt:4_pInt) & !*** store new maximum dipole height in state forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & - state(g,ip,el)%p(iD(s,c,matID)) = dUpper(s,c) + state(ipc,ip,el)%p(iD(s,c,matID)) = dUpper(s,c) @@ -2084,7 +2077,7 @@ forall (s = 1:ns, c = 1_pInt:2_pInt) & #ifndef _OPENMP if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & - .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)& + .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation remobilization', deltaRhoRemobilization(1:ns,1:8) write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> dipole dissociation by stress increase', deltaRhoDipole2SingleStress @@ -2099,7 +2092,7 @@ end subroutine constitutive_nonlocal_deltaState !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -function constitutive_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, state, state0, timestep, subfrac, g,ip,el) +function constitutive_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, state, state0, timestep, subfrac, ipc,ip,el) use prec, only: DAMASK_NaN use numerics, only: numerics_integrationMode, & @@ -2144,7 +2137,7 @@ use lattice, only: lattice_Sslip_v, & implicit none !*** input variables -integer(pInt), intent(in) :: g, & !< current grain number +integer(pInt), intent(in) :: ipc, & !< current grain number ip, & !< current integration point el !< current element number real(pReal), intent(in) :: Temperature, & !< temperature @@ -2162,7 +2155,7 @@ type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), in !*** input/output variables !*** output variables -real(pReal), dimension(constitutive_nonlocal_sizeDotState(phase_plasticityInstance(material_phase(g,ip,el)))) :: & +real(pReal), dimension(constitutive_nonlocal_sizeDotState(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & constitutive_nonlocal_dotState !< evolution of state variables / microstructure !*** local variables @@ -2184,38 +2177,38 @@ integer(pInt) matID, & !< current inst s, & !< index of my current slip system sLattice, & !< index of my current slip system according to lattice order deads -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),10) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),10) :: & rhoDot, & !< density evolution rhoDotMultiplication, & !< density evolution by multiplication rhoDotFlux, & !< density evolution by flux rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide) rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation rhoDotThermalAnnihilation !< density evolution by thermal annihilation -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),8) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: & rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) rhoSglOriginal, & neighbor_rhoSgl, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) rhoSgl0, & !< single dislocation densities at start of cryst inc (positive/negative screw and edge without dipoles) my_rhoSgl !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),4) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: & v, & !< current dislocation glide velocity v0, & !< dislocation glide velocity at start of cryst inc my_v, & !< dislocation glide velocity of central ip neighbor_v, & !< dislocation glide velocity of enighboring ip gdot !< shear rates -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & rhoForest, & !< forest dislocation density tauThreshold, & !< threshold shear stress tau, & !< current resolved shear stress tauBack, & !< current back stress from pileups on same slip system vClimb, & !< climb velocity of edge dipoles nSources -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),2) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: & rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) rhoDipOriginal, & dLower, & !< minimum stable dipole distance for edges and screws dUpper !< current maximum stable dipole distance for edges and screws -real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),4) :: & +real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: & m !< direction of dislocation motion real(pReal), dimension(3,3) :: my_F, & !< my total deformation gradient neighbor_F, & !< total deformation gradient of my neighbor @@ -2237,16 +2230,16 @@ logical considerEnteringFlux, & #ifndef _OPENMP if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt & - .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)& + .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then write(6,*) - write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_dotState at el ip g ',el,ip,g + write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_dotState at el ip ipc ',el,ip,ipc write(6,*) endif #endif -matID = phase_plasticityInstance(material_phase(g,ip,el)) +matID = phase_plasticityInstance(material_phase(ipc,ip,el)) structID = constitutive_nonlocal_structure(matID) ns = totalNslip(matID) @@ -2258,16 +2251,16 @@ gdot = 0.0_pReal forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) - rhoSgl(s,t) = max(state(g,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal) ! ensure positive single mobile densities - rhoSgl(s,t+4_pInt) = state(g,ip,el)%p(iRhoB(s,t,matID)) - v(s,t) = state(g,ip,el)%p(iV(s,t,matID)) + rhoSgl(s,t) = max(state(ipc,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal) ! ensure positive single mobile densities + rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,matID)) + v(s,t) = state(ipc,ip,el)%p(iV(s,t,matID)) endforall forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) - rhoDip(s,c) = max(state(g,ip,el)%p(iRhoD(s,c,matID)), 0.0_pReal) ! ensure positive dipole densities + rhoDip(s,c) = max(state(ipc,ip,el)%p(iRhoD(s,c,matID)), 0.0_pReal) ! ensure positive dipole densities endforall -rhoForest = state(g,ip,el)%p(iRhoF(1:ns,matID)) -tauThreshold = state(g,ip,el)%p(iTauF(1:ns,matID)) -tauBack = state(g,ip,el)%p(iTauB(1:ns,matID)) +rhoForest = state(ipc,ip,el)%p(iRhoF(1:ns,matID)) +tauThreshold = state(ipc,ip,el)%p(iTauF(1:ns,matID)) +tauBack = state(ipc,ip,el)%p(iTauB(1:ns,matID)) rhoSglOriginal = rhoSgl rhoDipOriginal = rhoDip @@ -2280,9 +2273,9 @@ where (abs(rhoDip) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(matID) & if (numerics_timeSyncing) then forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) - rhoSgl0(s,t) = max(state0(g,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal) - rhoSgl0(s,t+4_pInt) = state0(g,ip,el)%p(iRhoB(s,t,matID)) - v0(s,t) = state0(g,ip,el)%p(iV(s,t,matID)) + rhoSgl0(s,t) = max(state0(ipc,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal) + rhoSgl0(s,t+4_pInt) = state0(ipc,ip,el)%p(iRhoB(s,t,matID)) + v0(s,t) = state0(ipc,ip,el)%p(iV(s,t,matID)) endforall where (abs(rhoSgl0) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(matID) & .or. abs(rhoSgl0) < significantRho(matID)) & @@ -2308,7 +2301,7 @@ forall (t = 1_pInt:4_pInt) & #ifndef _OPENMP if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt & - .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)& + .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> gdot / 1/s',gdot @@ -2364,16 +2357,16 @@ else endwhere do s = 1_pInt,ns if (nSources(s) < 1.0_pReal) then - if (sourceProbability(s,g,ip,el) > 1.0_pReal) then + if (sourceProbability(s,ipc,ip,el) > 1.0_pReal) then call random_number(rnd) - sourceProbability(s,g,ip,el) = rnd + sourceProbability(s,ipc,ip,el) = rnd !$OMP FLUSH(sourceProbability) endif - if (sourceProbability(s,g,ip,el) > 1.0_pReal - nSources(s)) then + if (sourceProbability(s,ipc,ip,el) > 1.0_pReal - nSources(s)) then rhoDotMultiplication(s,1:4) = sum(rhoSglOriginal(s,1:4) * abs(v(s,1:4))) / meshlength endif else - sourceProbability(s,g,ip,el) = 2.0_pReal + sourceProbability(s,ipc,ip,el) = 2.0_pReal rhoDotMultiplication(s,1:4) = & (sum(abs(gdot(s,1:2))) * fEdgeMultiplication(matID) + sum(abs(gdot(s,3:4)))) & / burgers(s,matID) * sqrt(rhoForest(s)) / lambda0(s,matID) @@ -2381,7 +2374,7 @@ else enddo #ifndef _OPENMP if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & - .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)& + .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then write(6,'(a,/,4(12x,12(f12.5,1x),/))') '<< CONST >> sources', nSources write(6,*) @@ -2401,7 +2394,7 @@ endif rhoDotFlux = 0.0_pReal -if (.not. phase_localPlasticity(material_phase(g,ip,el))) then ! only for nonlocal plasticity +if (.not. phase_localPlasticity(material_phase(ipc,ip,el))) then ! only for nonlocal plasticity !*** check CFL (Courant-Friedrichs-Lewy) condition for flux @@ -2433,8 +2426,8 @@ if (.not. phase_localPlasticity(material_phase(g,ip,el))) then m(1:3,1:ns,3) = -lattice_st(1:3, slipSystemLattice(1:ns,matID), structID) m(1:3,1:ns,4) = lattice_st(1:3, slipSystemLattice(1:ns,matID), structID) - my_Fe = Fe(1:3,1:3,g,ip,el) - my_F = math_mul33x33(my_Fe, Fp(1:3,1:3,g,ip,el)) + my_Fe = Fe(1:3,1:3,ipc,ip,el) + my_F = math_mul33x33(my_Fe, Fp(1:3,1:3,ipc,ip,el)) do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) ! loop through my neighbors neighbor_el = mesh_ipNeighborhood(1,n,ip,el) @@ -2447,9 +2440,9 @@ if (.not. phase_localPlasticity(material_phase(g,ip,el))) then opposite_n = mesh_ipNeighborhood(3,opposite_neighbor,ip,el) if (neighbor_n > 0_pInt) then ! if neighbor exists, average deformation gradient - neighbor_instance = phase_plasticityInstance(material_phase(g,neighbor_ip,neighbor_el)) - neighbor_Fe = Fe(1:3,1:3,g,neighbor_ip,neighbor_el) - neighbor_F = math_mul33x33(neighbor_Fe, Fp(1:3,1:3,g,neighbor_ip,neighbor_el)) + neighbor_instance = phase_plasticityInstance(material_phase(ipc,neighbor_ip,neighbor_el)) + neighbor_Fe = Fe(1:3,1:3,ipc,neighbor_ip,neighbor_el) + neighbor_F = math_mul33x33(neighbor_Fe, Fp(1:3,1:3,ipc,neighbor_ip,neighbor_el)) Favg = 0.5_pReal * (my_F + neighbor_F) else ! if no neighbor, take my value as average Favg = my_F @@ -2471,17 +2464,17 @@ if (.not. phase_localPlasticity(material_phase(g,ip,el))) then endif if (considerEnteringFlux) then - if(numerics_timeSyncing .and. (subfrac(g,neighbor_ip,neighbor_el) /= subfrac(g,ip,el))) then ! for timesyncing: in case of a timestep at the interface we have to use "state0" to make sure that fluxes n both sides are equal + if(numerics_timeSyncing .and. (subfrac(ipc,neighbor_ip,neighbor_el) /= subfrac(ipc,ip,el))) then ! for timesyncing: in case of a timestep at the interface we have to use "state0" to make sure that fluxes n both sides are equal forall (s = 1:ns, t = 1_pInt:4_pInt) - neighbor_v(s,t) = state0(g,neighbor_ip,neighbor_el)%p(iV(s,t,neighbor_instance)) - neighbor_rhoSgl(s,t) = max(state0(g,neighbor_ip,neighbor_el)%p(iRhoU(s,t,neighbor_instance)), 0.0_pReal) - neighbor_rhoSgl(s,t+4_pInt) = state0(g,neighbor_ip,neighbor_el)%p(iRhoB(s,t,neighbor_instance)) + neighbor_v(s,t) = state0(ipc,neighbor_ip,neighbor_el)%p(iV(s,t,neighbor_instance)) + neighbor_rhoSgl(s,t) = max(state0(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,t,neighbor_instance)), 0.0_pReal) + neighbor_rhoSgl(s,t+4_pInt) = state0(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,t,neighbor_instance)) endforall else forall (s = 1:ns, t = 1_pInt:4_pInt) - neighbor_v(s,t) = state(g,neighbor_ip,neighbor_el)%p(iV(s,t,neighbor_instance)) - neighbor_rhoSgl(s,t) = max(state(g,neighbor_ip,neighbor_el)%p(iRhoU(s,t,neighbor_instance)), 0.0_pReal) - neighbor_rhoSgl(s,t+4_pInt) = state(g,neighbor_ip,neighbor_el)%p(iRhoB(s,t,neighbor_instance)) + neighbor_v(s,t) = state(ipc,neighbor_ip,neighbor_el)%p(iV(s,t,neighbor_instance)) + neighbor_rhoSgl(s,t) = max(state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,t,neighbor_instance)), 0.0_pReal) + neighbor_rhoSgl(s,t+4_pInt) = state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,t,neighbor_instance)) endforall endif where (abs(neighbor_rhoSgl) * mesh_ipVolume(neighbor_ip,neighbor_el) ** 0.667_pReal & @@ -2541,11 +2534,11 @@ if (.not. phase_localPlasticity(material_phase(g,ip,el))) then my_rhoSgl = rhoSgl my_v = v if(numerics_timeSyncing) then - if (subfrac(g,ip,el) == 0.0_pReal) then + if (subfrac(ipc,ip,el) == 0.0_pReal) then my_rhoSgl = rhoSgl0 my_v = v0 elseif (neighbor_n > 0_pInt) then - if (subfrac(g,neighbor_ip,neighbor_el) == 0.0_pReal) then + if (subfrac(ipc,neighbor_ip,neighbor_el) == 0.0_pReal) then my_rhoSgl = rhoSgl0 my_v = v0 endif @@ -2658,18 +2651,18 @@ rhoDot = rhoDotFlux & + rhoDotThermalAnnihilation if (numerics_integrationMode == 1_pInt) then ! save rates for output if in central integration mode - rhoDotFluxOutput(1:ns,1:8,g,ip,el) = rhoDotFlux(1:ns,1:8) - rhoDotMultiplicationOutput(1:ns,1:2,g,ip,el) = rhoDotMultiplication(1:ns,[1,3]) - rhoDotSingle2DipoleGlideOutput(1:ns,1:2,g,ip,el) = rhoDotSingle2DipoleGlide(1:ns,9:10) - rhoDotAthermalAnnihilationOutput(1:ns,1:2,g,ip,el) = rhoDotAthermalAnnihilation(1:ns,9:10) - rhoDotThermalAnnihilationOutput(1:ns,1:2,g,ip,el) = rhoDotThermalAnnihilation(1:ns,9:10) - rhoDotEdgeJogsOutput(1:ns,g,ip,el) = 2.0_pReal * rhoDotThermalAnnihilation(1:ns,1) + rhoDotFluxOutput(1:ns,1:8,ipc,ip,el) = rhoDotFlux(1:ns,1:8) + rhoDotMultiplicationOutput(1:ns,1:2,ipc,ip,el) = rhoDotMultiplication(1:ns,[1,3]) + rhoDotSingle2DipoleGlideOutput(1:ns,1:2,ipc,ip,el) = rhoDotSingle2DipoleGlide(1:ns,9:10) + rhoDotAthermalAnnihilationOutput(1:ns,1:2,ipc,ip,el) = rhoDotAthermalAnnihilation(1:ns,9:10) + rhoDotThermalAnnihilationOutput(1:ns,1:2,ipc,ip,el) = rhoDotThermalAnnihilation(1:ns,9:10) + rhoDotEdgeJogsOutput(1:ns,ipc,ip,el) = 2.0_pReal * rhoDotThermalAnnihilation(1:ns,1) endif #ifndef _OPENMP if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & - .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)& + .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> dislocation multiplication', rhoDotMultiplication(1:ns,1:4) * timestep write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation flux', rhoDotFlux(1:ns,1:8) * timestep @@ -2892,7 +2885,7 @@ end subroutine constitutive_nonlocal_updateCompatibility !********************************************************************* !* calculates quantities characterizing the microstructure * !********************************************************************* -function constitutive_nonlocal_dislocationstress(state, Fe, g, ip, el) +pure function constitutive_nonlocal_dislocationstress(state, Fe, ipc, ip, el) use math, only: math_mul33x33, & math_mul33x3, & @@ -2917,7 +2910,7 @@ implicit none !*** input variables -integer(pInt), intent(in) :: g, & ! current grain ID +integer(pInt), intent(in) :: ipc, & ! current grain ID ip, & ! current integration point el ! current element real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & @@ -2975,11 +2968,11 @@ real(pReal), dimension(2,2,maxval(totalNslip)) :: & neighbor_rhoExcess ! excess density at neighbor material point (edge/screw,mobile/dead,slipsystem) real(pReal), dimension(2,maxval(totalNslip)) :: & rhoExcessDead -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),8) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: & rhoSgl ! single dislocation density (edge+, edge-, screw+, screw-, used edge+, used edge-, used screw+, used screw-) logical inversionError -phase = material_phase(g,ip,el) +phase = material_phase(ipc,ip,el) instance = phase_plasticityInstance(phase) latticeStruct = constitutive_nonlocal_structure(instance) ns = totalNslip(instance) @@ -2989,8 +2982,8 @@ ns = totalNslip(instance) !*** get basic states forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) - rhoSgl(s,t) = max(state(g,ip,el)%p(iRhoU(s,t,instance)), 0.0_pReal) ! ensure positive single mobile densities - rhoSgl(s,t+4_pInt) = state(g,ip,el)%p(iRhoB(s,t,instance)) + rhoSgl(s,t) = max(state(ipc,ip,el)%p(iRhoU(s,t,instance)), 0.0_pReal) ! ensure positive single mobile densities + rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,instance)) endforall @@ -3001,7 +2994,7 @@ endforall constitutive_nonlocal_dislocationstress = 0.0_pReal if (.not. phase_localPlasticity(phase)) then - call math_invert33(Fe(1:3,1:3,g,ip,el), invFe, detFe, inversionError) + call math_invert33(Fe(1:3,1:3,ipc,ip,el), invFe, detFe, inversionError) !* in case of periodic surfaces we have to find out how many periodic images in each direction we need @@ -3025,7 +3018,7 @@ if (.not. phase_localPlasticity(phase)) then do neighbor_el = 1_pInt,mesh_NcpElems ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))) - neighbor_phase = material_phase(g,neighbor_ip,neighbor_el) + neighbor_phase = material_phase(ipc,neighbor_ip,neighbor_el) if (phase_localPlasticity(neighbor_phase)) then cycle endif @@ -3035,10 +3028,10 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el)) call math_invert33(Fe(1:3,1:3,1,neighbor_ip,neighbor_el), neighbor_invFe, detFe, inversionError) neighbor_ipVolumeSideLength = mesh_ipVolume(neighbor_ip,neighbor_el) ** (1.0_pReal/3.0_pReal) ! reference volume used here forall (s = 1_pInt:neighbor_ns, c = 1_pInt:2_pInt) - neighbor_rhoExcess(c,1,s) = state(g,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_instance)) & ! positive mobiles - - state(g,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_instance)) ! negative mobiles - neighbor_rhoExcess(c,2,s) = abs(state(g,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c-1,neighbor_instance))) & ! positive deads - - abs(state(g,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c,neighbor_instance))) ! negative deads + neighbor_rhoExcess(c,1,s) = state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_instance)) & ! positive mobiles + - state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_instance)) ! negative mobiles + neighbor_rhoExcess(c,2,s) = abs(state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c-1,neighbor_instance))) & ! positive deads + - abs(state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c,neighbor_instance))) ! negative deads endforall Tdislo_neighborLattice = 0.0_pReal do deltaX = periodicImages(1,1),periodicImages(2,1) @@ -3097,8 +3090,8 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el)) cycle elseif (j > 1_pInt) then x = connection_neighborSlip(1) + sign(0.5_pReal * segmentLength, & - state(g,neighbor_ip,neighbor_el)%p(iRhoB(s,1,neighbor_instance)) & - - state(g,neighbor_ip,neighbor_el)%p(iRhoB(s,2,neighbor_instance))) + state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,1,neighbor_instance)) & + - state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2,neighbor_instance))) xsquare = x * x endif @@ -3143,8 +3136,8 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el)) cycle elseif (j > 1_pInt) then y = connection_neighborSlip(2) + sign(0.5_pReal * segmentLength, & - state(g,neighbor_ip,neighbor_el)%p(iRhoB(s,3,neighbor_instance)) & - - state(g,neighbor_ip,neighbor_el)%p(iRhoB(s,4,neighbor_instance))) + state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,3,neighbor_instance)) & + - state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,4,neighbor_instance))) ysquare = y * y endif @@ -3197,8 +3190,8 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el)) else forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & - rhoExcessDead(c,s) = state(g,ip,el)%p(iRhoB(s,2*c-1,instance)) & ! positive deads (here we use symmetry: if this has negative sign it is treated as negative density at positive position instead of positive density at negative position) - + state(g,ip,el)%p(iRhoB(s,2*c,instance)) ! negative deads (here we use symmetry: if this has negative sign it is treated as positive density at positive position instead of negative density at negative position) + rhoExcessDead(c,s) = state(ipc,ip,el)%p(iRhoB(s,2*c-1,instance)) & ! positive deads (here we use symmetry: if this has negative sign it is treated as negative density at positive position instead of positive density at negative position) + + state(ipc,ip,el)%p(iRhoB(s,2*c,instance)) ! negative deads (here we use symmetry: if this has negative sign it is treated as positive density at positive position instead of negative density at negative position) do s = 1_pInt,ns if (all(abs(rhoExcessDead(:,s)) < significantRho(instance))) then @@ -3240,80 +3233,81 @@ endif end function constitutive_nonlocal_dislocationstress -!********************************************************************* -!* return array of constitutive results * -!********************************************************************* -function constitutive_nonlocal_postResults(Tstar_v, Fe, Temperature, dt, state, dotState, g,ip,el) +!-------------------------------------------------------------------------------------------------- +!> @brief return array of constitutive results +!-------------------------------------------------------------------------------------------------- +pure function constitutive_nonlocal_postResults(Tstar_v,Fe,state,dotState,ipc,ip,el) + use math, only: & + math_mul6x6, & + math_mul33x3, & + math_mul33x33, & + pi + use mesh, only: & + mesh_NcpElems, & + mesh_maxNips + use material, only: & + homogenization_maxNgrains, & + material_phase, & + phase_plasticityInstance, & + phase_Noutput + use lattice, only: & + lattice_Sslip_v, & + lattice_sd, & + lattice_st, & + lattice_sn + + implicit none + real(pReal), dimension(6), intent(in) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + Fe !< elastic deformation gradient + type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + state !< microstructure state + type(p_vec), intent(in) :: dotState ! evolution rate of microstructural state + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + + real(pReal), dimension(constitutive_nonlocal_sizePostResults(& + phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + constitutive_nonlocal_postResults + + integer(pInt) :: & + matID, & !< current instance of this plasticity + structID, & !< current lattice structure + ns, & !< short notation for the total number of active slip systems + c, & !< character of dislocation + cs, & !< constitutive result index + o, & !< index of current output + t, & !< type of dislocation + s, & !< index of my current slip system + sLattice !< index of my current slip system according to lattice order + real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: & + rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) + rhoDotSgl !< evolution rate of single dislocation densities (positive/negative screw and edge without dipoles) + real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: & + gdot, & !< shear rates + v !< velocities + real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + rhoForest, & !< forest dislocation density + tauThreshold, & !< threshold shear stress + tau, & !< current resolved shear stress + tauBack !< back stress from pileups on same slip system + real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: & + rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) + rhoDotDip, & !< evolution rate of dipole dislocation densities (screw and edge dipoles) + dLower, & !< minimum stable dipole distance for edges and screws + dUpper !< current maximum stable dipole distance for edges and screws + real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: & + m, & !< direction of dislocation motion for edge and screw (unit vector) + m_currentconf !< direction of dislocation motion for edge and screw (unit vector) in current configuration + real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + n_currentconf !< slip system normal (unit vector) in current configuration + real(pReal), dimension(3,3) :: & + sigma -use math, only: math_mul6x6, & - math_mul33x3, & - math_mul33x33, & - pi -use mesh, only: mesh_NcpElems, & - mesh_maxNips -use material, only: homogenization_maxNgrains, & - material_phase, & - phase_plasticityInstance, & - phase_Noutput -use lattice, only: lattice_Sslip_v, & - lattice_sd, & - lattice_st, & - lattice_sn - -implicit none - -!*** input variables -integer(pInt), intent(in) :: g, & ! current grain number - ip, & ! current integration point - el ! current element number - -real(pReal), intent(in) :: Temperature, & ! temperature - dt ! time increment -real(pReal), dimension(6), intent(in) :: Tstar_v ! current 2nd Piola-Kirchhoff stress in Mandel notation -real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & - Fe ! elastic deformation gradient -type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & - state ! current microstructural state -type(p_vec), intent(in) :: dotState ! evolution rate of microstructural state - -!*** output variables -real(pReal), dimension(constitutive_nonlocal_sizePostResults(phase_plasticityInstance(material_phase(g,ip,el)))) :: & - constitutive_nonlocal_postResults - -!*** local variables -integer(pInt) matID, & ! current instance of this plasticity - structID, & ! current lattice structure - ns, & ! short notation for the total number of active slip systems - c, & ! character of dislocation - cs, & ! constitutive result index - o, & ! index of current output - t, & ! type of dislocation - s, & ! index of my current slip system - sLattice ! index of my current slip system according to lattice order -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),8) :: & - rhoSgl, & ! current single dislocation densities (positive/negative screw and edge without dipoles) - rhoDotSgl ! evolution rate of single dislocation densities (positive/negative screw and edge without dipoles) -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),4) :: & - gdot, & ! shear rates - v ! velocities -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))) :: & - rhoForest, & ! forest dislocation density - tauThreshold, & ! threshold shear stress - tau, & ! current resolved shear stress - tauBack ! back stress from pileups on same slip system -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),2) :: & - rhoDip, & ! current dipole dislocation densities (screw and edge dipoles) - rhoDotDip, & ! evolution rate of dipole dislocation densities (screw and edge dipoles) - dLower, & ! minimum stable dipole distance for edges and screws - dUpper ! current maximum stable dipole distance for edges and screws -real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),2) :: & - m, & ! direction of dislocation motion for edge and screw (unit vector) - m_currentconf ! direction of dislocation motion for edge and screw (unit vector) in current configuration -real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))) :: & - n_currentconf ! slip system normal (unit vector) in current configuration -real(pReal), dimension(3,3) :: sigma - -matID = phase_plasticityInstance(material_phase(g,ip,el)) +matID = phase_plasticityInstance(material_phase(ipc,ip,el)) structID = constitutive_nonlocal_structure(matID) ns = totalNslip(matID) @@ -3324,19 +3318,19 @@ constitutive_nonlocal_postResults = 0.0_pReal !* short hand notations for state variables forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) - rhoSgl(s,t) = state(g,ip,el)%p(iRhoU(s,t,matID)) - rhoSgl(s,t+4_pInt) = state(g,ip,el)%p(iRhoB(s,t,matID)) - v(s,t) = state(g,ip,el)%p(iV(s,t,matID)) + rhoSgl(s,t) = state(ipc,ip,el)%p(iRhoU(s,t,matID)) + rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,matID)) + v(s,t) = state(ipc,ip,el)%p(iV(s,t,matID)) rhoDotSgl(s,t) = dotState%p(iRhoU(s,t,matID)) rhoDotSgl(s,t+4_pInt) = dotState%p(iRhoB(s,t,matID)) endforall forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) - rhoDip(s,c) = state(g,ip,el)%p(iRhoD(s,c,matID)) + rhoDip(s,c) = state(ipc,ip,el)%p(iRhoD(s,c,matID)) rhoDotDip(s,c) = dotState%p(iRhoD(s,c,matID)) endforall -rhoForest = state(g,ip,el)%p(iRhoF(1:ns,matID)) -tauThreshold = state(g,ip,el)%p(iTauF(1:ns,matID)) -tauBack = state(g,ip,el)%p(iTauB(1:ns,matID)) +rhoForest = state(ipc,ip,el)%p(iRhoF(1:ns,matID)) +tauThreshold = state(ipc,ip,el)%p(iTauF(1:ns,matID)) +tauBack = state(ipc,ip,el)%p(iTauB(1:ns,matID)) @@ -3371,13 +3365,13 @@ dUpper = max(dUpper,dLower) m(1:3,1:ns,1) = lattice_sd(1:3,slipSystemLattice(1:ns,matID),structID) m(1:3,1:ns,2) = -lattice_st(1:3,slipSystemLattice(1:ns,matID),structID) forall (c = 1_pInt:2_pInt, s = 1_pInt:ns) & - m_currentconf(1:3,s,c) = math_mul33x3(Fe(1:3,1:3,g,ip,el), m(1:3,s,c)) + m_currentconf(1:3,s,c) = math_mul33x3(Fe(1:3,1:3,ipc,ip,el), m(1:3,s,c)) forall (s = 1_pInt:ns) & - n_currentconf(1:3,s) = math_mul33x3(Fe(1:3,1:3,g,ip,el), & + n_currentconf(1:3,s) = math_mul33x3(Fe(1:3,1:3,ipc,ip,el), & lattice_sn(1:3,slipSystemLattice(s,matID),structID)) -do o = 1_pInt,phase_Noutput(material_phase(g,ip,el)) +do o = 1_pInt,phase_Noutput(material_phase(ipc,ip,el)) select case(constitutive_nonlocal_output(o,matID)) case ('rho') @@ -3557,66 +3551,66 @@ do o = 1_pInt,phase_Noutput(material_phase(g,ip,el)) cs = cs + ns case ('rho_dot_gen') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,1,g,ip,el) & - + rhoDotMultiplicationOutput(1:ns,2,g,ip,el) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,1,ipc,ip,el) & + + rhoDotMultiplicationOutput(1:ns,2,ipc,ip,el) cs = cs + ns case ('rho_dot_gen_edge') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,1,g,ip,el) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,1,ipc,ip,el) cs = cs + ns case ('rho_dot_gen_screw') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,2,g,ip,el) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,2,ipc,ip,el) cs = cs + ns case ('rho_dot_sgl2dip') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,1,g,ip,el) & - + rhoDotSingle2DipoleGlideOutput(1:ns,2,g,ip,el) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,1,ipc,ip,el) & + + rhoDotSingle2DipoleGlideOutput(1:ns,2,ipc,ip,el) cs = cs + ns case ('rho_dot_sgl2dip_edge') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,1,g,ip,el) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,1,ipc,ip,el) cs = cs + ns case ('rho_dot_sgl2dip_screw') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,2,g,ip,el) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,2,ipc,ip,el) cs = cs + ns case ('rho_dot_ann_ath') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotAthermalAnnihilationOutput(1:ns,1,g,ip,el) & - + rhoDotAthermalAnnihilationOutput(1:ns,2,g,ip,el) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotAthermalAnnihilationOutput(1:ns,1,ipc,ip,el) & + + rhoDotAthermalAnnihilationOutput(1:ns,2,ipc,ip,el) cs = cs + ns case ('rho_dot_ann_the') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,1,g,ip,el) & - + rhoDotThermalAnnihilationOutput(1:ns,2,g,ip,el) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,1,ipc,ip,el) & + + rhoDotThermalAnnihilationOutput(1:ns,2,ipc,ip,el) cs = cs + ns case ('rho_dot_ann_the_edge') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,1,g,ip,el) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,1,ipc,ip,el) cs = cs + ns case ('rho_dot_ann_the_screw') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,2,g,ip,el) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,2,ipc,ip,el) cs = cs + ns case ('rho_dot_edgejogs') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotEdgeJogsOutput(1:ns,g,ip,el) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotEdgeJogsOutput(1:ns,ipc,ip,el) cs = cs + ns case ('rho_dot_flux') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,1:4,g,ip,el),2) & - + sum(abs(rhoDotFluxOutput(1:ns,5:8,g,ip,el)),2) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,1:4,ipc,ip,el),2) & + + sum(abs(rhoDotFluxOutput(1:ns,5:8,ipc,ip,el)),2) cs = cs + ns case ('rho_dot_flux_edge') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,1:2,g,ip,el),2) & - + sum(abs(rhoDotFluxOutput(1:ns,5:6,g,ip,el)),2) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,1:2,ipc,ip,el),2) & + + sum(abs(rhoDotFluxOutput(1:ns,5:6,ipc,ip,el)),2) cs = cs + ns case ('rho_dot_flux_screw') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,3:4,g,ip,el),2) & - + sum(abs(rhoDotFluxOutput(1:ns,7:8,g,ip,el)),2) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,3:4,ipc,ip,el),2) & + + sum(abs(rhoDotFluxOutput(1:ns,7:8,ipc,ip,el)),2) cs = cs + ns case ('velocity_edge_pos') @@ -3716,7 +3710,7 @@ do o = 1_pInt,phase_Noutput(material_phase(g,ip,el)) cs = cs + ns case('dislocationstress') - sigma = constitutive_nonlocal_dislocationstress(state, Fe, g, ip, el) + sigma = constitutive_nonlocal_dislocationstress(state, Fe, ipc, ip, el) constitutive_nonlocal_postResults(cs+1_pInt) = sigma(1,1) constitutive_nonlocal_postResults(cs+2_pInt) = sigma(2,2) constitutive_nonlocal_postResults(cs+3_pInt) = sigma(3,3) @@ -3726,7 +3720,7 @@ do o = 1_pInt,phase_Noutput(material_phase(g,ip,el)) cs = cs + 6_pInt case('accumulatedshear') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = state(g,ip,el)%p(iGamma(1:ns,matID)) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = state(ipc,ip,el)%p(iGamma(1:ns,matID)) cs = cs + ns case('boundarylayer') @@ -3742,6 +3736,6 @@ do o = 1_pInt,phase_Noutput(material_phase(g,ip,el)) end select enddo -endfunction +end function constitutive_nonlocal_postResults -END MODULE +end module constitutive_nonlocal diff --git a/code/constitutive_phenopowerlaw.f90 b/code/constitutive_phenopowerlaw.f90 index 006645a6e..588fba5c9 100644 --- a/code/constitutive_phenopowerlaw.f90 +++ b/code/constitutive_phenopowerlaw.f90 @@ -646,7 +646,7 @@ end function constitutive_phenopowerlaw_aTolState !-------------------------------------------------------------------------------------------------- !> @brief returns the homogenized elasticity matrix !-------------------------------------------------------------------------------------------------- -pure function constitutive_phenopowerlaw_homogenizedC(state,ipc,ip,el) +pure function constitutive_phenopowerlaw_homogenizedC(ipc,ip,el) use prec, only: & p_vec use mesh, only: & @@ -664,8 +664,6 @@ pure function constitutive_phenopowerlaw_homogenizedC(state,ipc,ip,el) ipc, & !< component-ID of integration point ip, & !< integration point el !< element - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & - state !< microstructure state constitutive_phenopowerlaw_homogenizedC = constitutive_phenopowerlaw_Cslip_66(1:6,1:6,& phase_plasticityInstance(material_phase(ipc,ip,el))) @@ -676,8 +674,7 @@ end function constitutive_phenopowerlaw_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,& - temperature,state,ipc,ip,el) +pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,state,ipc,ip,el) use prec, only: & p_vec use math, only: & @@ -709,8 +706,6 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar real(pReal), dimension(6), intent(in) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal), intent(in) :: & - temperature !< temperature at IP integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -831,7 +826,7 @@ end subroutine constitutive_phenopowerlaw_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -function constitutive_phenopowerlaw_dotState(Tstar_v,temperature,state,ipc,ip,el) +function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el) use prec, only: & p_vec use lattice, only: & @@ -854,8 +849,6 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,temperature,state,ipc,ip,el implicit none real(pReal), dimension(6), intent(in) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal), intent(in) :: & - temperature !< temperature at integration point integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -1003,7 +996,7 @@ end function constitutive_phenopowerlaw_dotState !-------------------------------------------------------------------------------------------------- !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- -pure function constitutive_phenopowerlaw_postResults(Tstar_v,temperature,dt,state,ipc,ip,el) +pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el) use prec, only: & p_vec use mesh, only: & @@ -1029,9 +1022,6 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,temperature,dt,stat implicit none real(pReal), dimension(6), intent(in) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal), intent(in) :: & - temperature, & !< temperature at integration point - dt integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point diff --git a/code/constitutive_titanmod.f90 b/code/constitutive_titanmod.f90 index fdc9ff458..86790d37d 100644 --- a/code/constitutive_titanmod.f90 +++ b/code/constitutive_titanmod.f90 @@ -1144,11 +1144,11 @@ real(pReal), dimension(constitutive_titanmod_totalNtwin(phase_plasticityInstance !-------------------------------------------------------------------------------------------------- ! homogenized elasticity matrix - constitutive_titanmod_homogenizedC = (1.0_pReal-sumf)*constitutive_titanmod_Cslip_66(:,:,matID) + constitutive_titanmod_homogenizedC = (1.0_pReal-sumf)*constitutive_titanmod_Cslip_66(1:6,1:6,matID) do i=1_pInt,nt - constitutive_titanmod_homogenizedC = & - constitutive_titanmod_homogenizedC + volumefraction_PerTwinSys(i)*constitutive_titanmod_Ctwin_66(:,:,i,matID) - + constitutive_titanmod_homogenizedC = constitutive_titanmod_homogenizedC & + + volumefraction_PerTwinSys(i)*& + constitutive_titanmod_Ctwin_66(1:6,1:6,i,matID) enddo end function constitutive_titanmod_homogenizedC @@ -1717,7 +1717,7 @@ end function constitutive_titanmod_dotState !-------------------------------------------------------------------------------------------------- !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- -pure function constitutive_titanmod_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el) +pure function constitutive_titanmod_postResults(state,ipc,ip,el) use prec, only: & p_vec use mesh, only: & @@ -1730,24 +1730,21 @@ pure function constitutive_titanmod_postResults(Tstar_v,Temperature,dt,state,ipc phase_Noutput implicit none - real(pReal), dimension(6), intent(in) :: & - Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal), intent(in) :: & - temperature, & !< temperature at integration point - dt integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & state !< microstructure state + real(pReal), dimension(constitutive_titanmod_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + constitutive_titanmod_postResults + integer(pInt) :: & matID, structID,& ns,nt,& o,i,c real(pReal) :: sumf - real(pReal), dimension(constitutive_titanmod_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & - constitutive_titanmod_postResults + real(pReal), dimension(constitutive_titanmod_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & volumefraction_PerTwinSys diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 9fc7512c1..ebf0083a2 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -42,17 +42,17 @@ module crystallite crystallite_sizePostResult !< description not available integer(pInt), dimension(:,:,:), allocatable, private :: & crystallite_symmetryID !< crystallographic symmetry 1=cubic 2=hexagonal, needed in all orientation calcs - + + real(pReal), dimension(:,:), allocatable, public :: & + crystallite_temperature !< temperature (same on all components on one IP) + real(pReal), dimension(:,:,:), allocatable, public, protected :: & + crystallite_heat !< heat source real(pReal), dimension(:,:,:), allocatable, public :: & - crystallite_dt, & !< requested time increment of each grain - crystallite_Temperature, & !< Temp of each grain - crystallite_partionedTemperature0 !< Temp of each grain at start of homog inc + crystallite_dt !< requested time increment of each grain real(pReal), dimension(:,:,:), allocatable, private :: & crystallite_subdt, & !< substepped time increment of each grain crystallite_subFrac, & !< already calculated fraction of increment - crystallite_subStep, & !< size of next integration step - crystallite_subTemperature0, & !< Temp of each grain at start of crystallite inc - crystallite_dotTemperature !< evolution of Temperature of each grain + crystallite_subStep !< size of next integration step real(pReal), dimension(:,:,:,:), allocatable, public :: & crystallite_Tstar_v, & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step) crystallite_Tstar0_v, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc @@ -122,7 +122,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates and initialize per grain variables !-------------------------------------------------------------------------------------------------- -subroutine crystallite_init(Temperature) +subroutine crystallite_init(temperature) 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_info, & @@ -178,11 +178,12 @@ subroutine crystallite_init(Temperature) constitutive_nonlocal_structureName implicit none - real(pReal), intent(in) :: Temperature - integer(pInt), parameter :: myFile = 200_pInt, & - maxNchunks = 2_pInt + real(pReal), intent(in) :: temperature + integer(pInt), parameter :: & + MYUNIT = 200_pInt, & + MAXNCHUNKS = 2_pInt - integer(pInt), dimension(1+2*maxNchunks) :: positions + integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions integer(pInt) :: & g, & !< grain number i, & !< integration point number @@ -212,10 +213,8 @@ subroutine crystallite_init(Temperature) eMax = mesh_NcpElems nMax = mesh_maxNipNeighbors - allocate(crystallite_Temperature(gMax,iMax,eMax)); crystallite_Temperature = Temperature - allocate(crystallite_partionedTemperature0(gMax,iMax,eMax)); crystallite_partionedTemperature0 = 0.0_pReal - allocate(crystallite_subTemperature0(gMax,iMax,eMax)); crystallite_subTemperature0 = 0.0_pReal - allocate(crystallite_dotTemperature(gMax,iMax,eMax)); crystallite_dotTemperature = 0.0_pReal + allocate(crystallite_temperature(iMax,eMax)); crystallite_temperature = temperature + allocate(crystallite_heat(gMax,iMax,eMax)); crystallite_heat = 0.0_pReal allocate(crystallite_Tstar0_v(6,gMax,iMax,eMax)); crystallite_Tstar0_v = 0.0_pReal allocate(crystallite_partionedTstar0_v(6,gMax,iMax,eMax)); crystallite_partionedTstar0_v = 0.0_pReal allocate(crystallite_subTstar0_v(6,gMax,iMax,eMax)); crystallite_subTstar0_v = 0.0_pReal @@ -359,7 +358,6 @@ subroutine crystallite_init(Temperature) if(any(.not. crystallite_localPlasticity) .and. .not. usePingPong) call IO_error(601_pInt) ! exit if nonlocal but no ping-pong - crystallite_partionedTemperature0 = Temperature ! isothermal assumption crystallite_partionedFp0 = crystallite_Fp0 crystallite_partionedF0 = crystallite_F0 crystallite_partionedF = crystallite_F0 @@ -401,7 +399,7 @@ subroutine crystallite_init(Temperature) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1_pInt,myNgrains - call constitutive_microstructure(crystallite_Temperature(g,i,e), & + call constitutive_microstructure(temperature, & crystallite_Fe(1:3,1:3,g,i,e), crystallite_Fp(1:3,1:3,g,i,e),g,i,e) ! update dependent state variables to be consistent with basic states enddo enddo @@ -414,8 +412,8 @@ subroutine crystallite_init(Temperature) !-------------------------------------------------------------------------------------------------- ! debug output if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Temperature: ', shape(crystallite_Temperature) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_dotTemperature: ', shape(crystallite_dotTemperature) + write(6,'(a35,1x,7(i8,1x))') 'crystallite_Temperature: ', shape(crystallite_temperature) + write(6,'(a35,1x,7(i8,1x))') 'crystallite_heat: ', shape(crystallite_heat) write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fe: ', shape(crystallite_Fe) write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fp: ', shape(crystallite_Fp) write(6,'(a35,1x,7(i8,1x))') 'crystallite_Lp: ', shape(crystallite_Lp) @@ -423,12 +421,10 @@ subroutine crystallite_init(Temperature) write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fp0: ', shape(crystallite_Fp0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_Lp0: ', shape(crystallite_Lp0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedF: ', shape(crystallite_partionedF) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedTemp0: ', shape(crystallite_partionedTemperature0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedF0: ', shape(crystallite_partionedF0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedFp0: ', shape(crystallite_partionedFp0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedLp0: ', shape(crystallite_partionedLp0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_subF: ', shape(crystallite_subF) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subTemperature0: ', shape(crystallite_subTemperature0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_symmetryID: ', shape(crystallite_symmetryID) write(6,'(a35,1x,7(i8,1x))') 'crystallite_subF0: ', shape(crystallite_subF0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFe0: ', shape(crystallite_subFe0) @@ -554,8 +550,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) P_backup real(pReal), dimension(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & Tstar_v_backup - real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - Temperature_backup integer(pInt) :: & NiterationCrystallite, & ! number of iterations in crystallite loop e, & ! element index @@ -586,13 +580,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) real(pReal) :: counter - if(iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt& - .and. debug_e > 0 .and. debug_e <= mesh_NcpElems & - .and. debug_i > 0 .and. debug_i <= mesh_maxNips & - .and. debug_g > 0 .and. debug_g <= homogenization_maxNgrains) then + if(iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then !$OMP CRITICAL (write2out) write(6,'(/,a,i8,1x,i2,1x,i3)') '<< CRYST >> crystallite start at el ip g ', debug_e, debug_i, debug_g - write(6,'(a,/,12x,f14.9)') '<< CRYST >> Temp0', crystallite_partionedTemperature0(debug_g,debug_i,debug_e) write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> F0 ', & math_transpose33(crystallite_partionedF0(1:3,1:3,debug_g,debug_i,debug_e)) write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Fp0', & @@ -610,7 +600,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) myNgrains = homogenization_Ngrains(mesh_element(3,e)) forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & g = 1_pInt:myNgrains, crystallite_requested(g,i,e)) - crystallite_subTemperature0(g,i,e) = crystallite_partionedTemperature0(g,i,e) ! ...temperature constitutive_subState0(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructure crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_partionedFp0(1:3,1:3,g,i,e) ! ...plastic def grad crystallite_subLp0(1:3,1:3,g,i,e) = crystallite_partionedLp0(1:3,1:3,g,i,e) ! ...plastic velocity grad @@ -908,7 +897,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) stepIncreaseCryst * crystallite_subStep(g,i,e)) !$OMP FLUSH(crystallite_subStep) if (crystallite_subStep(g,i,e) > 0.0_pReal) then - crystallite_subTemperature0(g,i,e) = crystallite_Temperature(g,i,e) ! wind forward... crystallite_subF0(1:3,1:3,g,i,e) = crystallite_subF(1:3,1:3,g,i,e) ! ...def grad !$OMP FLUSH(crystallite_subF0) crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_Fp(1:3,1:3,g,i,e) ! ...plastic def grad @@ -956,7 +944,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) crystallite_subStep(g,i,e) = subStepSizeCryst * crystallite_subStep(g,i,e) ! cut step in half and restore... endif !$OMP FLUSH(crystallite_subStep) - crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) ! ...temperature crystallite_Fp(1:3,1:3,g,i,e) = crystallite_subFp0(1:3,1:3,g,i,e) ! ...plastic def grad !$OMP FLUSH(crystallite_Fp) crystallite_invFp(1:3,1:3,g,i,e) = math_inv33(crystallite_Fp(1:3,1:3,g,i,e)) @@ -1115,7 +1102,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) constitutive_dotState(g,i,e)%p(1:constitutive_sizeDotState(g,i,e)) ! ... dotStates, ... endforall enddo - Temperature_backup = crystallite_Temperature ! ... Temperature, ... F_backup = crystallite_subF ! ... and kinematics Fp_backup = crystallite_Fp InvFp_backup = crystallite_invFp @@ -1155,7 +1141,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) constitutive_dotState_backup(g,i,e)%p(1:constitutive_sizeDotState(g,i,e)) endforall enddo - crystallite_Temperature = Temperature_backup crystallite_Fp = Fp_backup crystallite_invFp = InvFp_backup crystallite_Fe = Fe_backup @@ -1172,7 +1157,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) constitutive_dotState_backup(g,i,e)%p(1:constitutive_sizeDotState(g,i,e)) endforall enddo - crystallite_Temperature = crystallite_subTemperature0 crystallite_Fp = crystallite_subFp0 crystallite_Fe = crystallite_subFe0 crystallite_Lp = crystallite_subLp0 @@ -1258,7 +1242,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) constitutive_dotState_backup(g,i,e)%p(1:constitutive_sizeDotState(g,i,e)) endforall enddo - crystallite_Temperature = Temperature_backup crystallite_subF = F_backup crystallite_Fp = Fp_backup crystallite_invFp = InvFp_backup @@ -1340,7 +1323,7 @@ end subroutine crystallite_stressAndItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief integrate stress, state and Temperature with 4th order explicit Runge Kutta method +!> @brief integrate stress, state with 4th order explicit Runge Kutta method !-------------------------------------------------------------------------------------------------- subroutine crystallite_integrateStateRK4() use prec, only: pInt, & @@ -1370,7 +1353,6 @@ subroutine crystallite_integrateStateRK4() constitutive_collectDotState, & constitutive_deltaState, & constitutive_collectDeltaState, & - constitutive_dotTemperature, & constitutive_microstructure implicit none @@ -1385,8 +1367,6 @@ subroutine crystallite_integrateStateRK4() integer(pInt), dimension(2) :: eIter ! bounds for element iteration integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration gIter ! bounds for grain iteration - real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - RK4dotTemperature ! evolution of Temperature of each grain for Runge Kutta integration logical :: singleRun ! flag indicating computation for single (g,i,e) triple eIter = FEsolving_execElem(1:2) @@ -1399,16 +1379,13 @@ subroutine crystallite_integrateStateRK4() ! --- FIRST RUNGE KUTTA STEP --- - RK4dotTemperature = 0.0_pReal ! initialize Runge-Kutta dotTemperature !$OMP PARALLEL PRIVATE(mySizeDotState) !$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 constitutive_RK4dotState(g,i,e)%p = 0.0_pReal ! initialize Runge-Kutta dotState if (crystallite_todo(g,i,e)) then call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & - crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) - crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & - crystallite_Temperature(g,i,e),g,i,e) + crystallite_Temperature(i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) endif enddo; enddo; enddo !$OMP ENDDO @@ -1416,8 +1393,7 @@ subroutine crystallite_integrateStateRK4() 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 - if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState - .or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature + if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) )then ! NaN occured in 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 @@ -1445,11 +1421,9 @@ subroutine crystallite_integrateStateRK4() mySizeDotState = constitutive_sizeDotState(g,i,e) if (n < 4) then constitutive_RK4dotState(g,i,e)%p = constitutive_RK4dotState(g,i,e)%p + weight(n)*constitutive_dotState(g,i,e)%p - RK4dotTemperature(g,i,e) = RK4dotTemperature(g,i,e) + weight(n)*crystallite_dotTemperature(g,i,e) elseif (n == 4) then constitutive_dotState(g,i,e)%p = (constitutive_RK4dotState(g,i,e)%p + & weight(n)*constitutive_dotState(g,i,e)%p) / 6.0_pReal ! use weighted RKdotState for final integration - crystallite_dotTemperature(g,i,e) = (RK4dotTemperature(g,i,e) + weight(n)*crystallite_dotTemperature(g,i,e)) / 6.0_pReal endif endif enddo; enddo; enddo @@ -1460,8 +1434,6 @@ subroutine crystallite_integrateStateRK4() mySizeDotState = constitutive_sizeDotState(g,i,e) constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) & + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e) * timeStepFraction(n) - crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) & - + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) * timeStepFraction(n) if (n == 4) then ! final integration step #ifndef _OPENMP if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & @@ -1502,7 +1474,7 @@ subroutine crystallite_integrateStateRK4() !$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 if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), & + call constitutive_microstructure(crystallite_Temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states endif enddo; enddo; enddo @@ -1534,10 +1506,8 @@ subroutine crystallite_integrateStateRK4() 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 call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & - crystallite_Temperature(g,i,e), timeStepFraction(n)*crystallite_subdt(g,i,e), & ! fraction of original timestep + crystallite_Temperature(i,e), timeStepFraction(n)*crystallite_subdt(g,i,e), & ! fraction of original timestep crystallite_subFrac, g,i,e) - crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & - crystallite_Temperature(g,i,e),g,i,e) endif enddo; enddo; enddo !$OMP ENDDO @@ -1545,8 +1515,7 @@ subroutine crystallite_integrateStateRK4() 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 - if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState - .or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature + if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p)) then ! NaN occured in 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 @@ -1591,7 +1560,7 @@ end subroutine crystallite_integrateStateRK4 !-------------------------------------------------------------------------------------------------- -!> @brief integrate stress, state and Temperature with 5th order Runge-Kutta Cash-Karp method with +!> @brief integrate stress, state with 5th order Runge-Kutta Cash-Karp method with !> adaptive step size (use 5th order solution to advance = "local extrapolation") !-------------------------------------------------------------------------------------------------- subroutine crystallite_integrateStateRKCK45() @@ -1605,7 +1574,6 @@ subroutine crystallite_integrateStateRKCK45() debug_g, & debug_StateLoopDistribution use numerics, only: rTol_crystalliteState, & - rTol_crystalliteTemperature, & numerics_integrationMode use FEsolving, only: FEsolving_execElem, & FEsolving_execIP @@ -1624,7 +1592,6 @@ subroutine crystallite_integrateStateRKCK45() constitutive_collectDotState, & constitutive_deltaState, & constitutive_collectDeltaState, & - constitutive_dotTemperature, & constitutive_microstructure implicit none @@ -1638,8 +1605,6 @@ subroutine crystallite_integrateStateRKCK45() integer(pInt), dimension(2) :: eIter ! bounds for element iteration integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration gIter ! bounds for grain iteration - real(pReal), dimension(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - RKCK45dotTemperature ! evolution of Temperature of each grain for Runge Kutta Cash Karp integration real(pReal), dimension(5,5), parameter :: A = reshape([& .2_pReal, .075_pReal, .3_pReal, -11.0_pReal/54.0_pReal, 1631.0_pReal/55296.0_pReal, & .0_pReal, .225_pReal, -.9_pReal, 2.5_pReal, 175.0_pReal/512.0_pReal, & @@ -1661,9 +1626,6 @@ subroutine crystallite_integrateStateRKCK45() real(pReal), dimension(constitutive_maxSizeDotState,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & stateResiduum, & ! residuum from evolution in micrstructure relStateResiduum ! relative residuum from evolution in microstructure - real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - temperatureResiduum, & ! residuum from evolution in temperature - relTemperatureResiduum ! relative residuum from evolution in temperature logical :: singleRun ! flag indicating computation for single (g,i,e) triple ! --- LOOP ITERATOR FOR ELEMENT, GRAIN, IP --- @@ -1689,9 +1651,7 @@ subroutine crystallite_integrateStateRKCK45() 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 call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & - crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) - crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & - crystallite_Temperature(g,i,e),g,i,e) + crystallite_Temperature(i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) endif enddo; enddo; enddo !$OMP ENDDO @@ -1699,8 +1659,7 @@ subroutine crystallite_integrateStateRKCK45() 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 - if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState - .or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature + if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p)) then ! NaN occured in 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 @@ -1727,7 +1686,6 @@ subroutine crystallite_integrateStateRKCK45() if (crystallite_todo(g,i,e)) then mySizeDotState = constitutive_sizeDotState(g,i,e) constitutive_RKCK45dotState(n,g,i,e)%p = constitutive_dotState(g,i,e)%p ! store Runge-Kutta dotState - RKCK45dotTemperature(n,g,i,e) = crystallite_dotTemperature(g,i,e) ! store Runge-Kutta dotTemperature endif enddo; enddo; enddo !$OMP ENDDO @@ -1735,40 +1693,25 @@ subroutine crystallite_integrateStateRKCK45() 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 if (n == 1) then ! NEED TO DO THE ADDITION IN THIS LENGTHY WAY BECAUSE OF PARALLELIZATION (CAN'T USE A REDUCTION CLAUSE ON A POINTER OR USER DEFINED TYPE) - constitutive_dotState(g,i,e)%p = a(1,1) * constitutive_RKCK45dotState(1,g,i,e)%p - crystallite_dotTemperature(g,i,e) = a(1,1) * RKCK45dotTemperature(1,g,i,e) + constitutive_dotState(g,i,e)%p = a(1,1) * constitutive_RKCK45dotState(1,g,i,e)%p elseif (n == 2) then constitutive_dotState(g,i,e)%p = a(1,2) * constitutive_RKCK45dotState(1,g,i,e)%p & + a(2,2) * constitutive_RKCK45dotState(2,g,i,e)%p - crystallite_dotTemperature(g,i,e) = a(1,2) * RKCK45dotTemperature(1,g,i,e) & - + a(2,2) * RKCK45dotTemperature(2,g,i,e) elseif (n == 3) then constitutive_dotState(g,i,e)%p = a(1,3) * constitutive_RKCK45dotState(1,g,i,e)%p & + a(2,3) * constitutive_RKCK45dotState(2,g,i,e)%p & + a(3,3) * constitutive_RKCK45dotState(3,g,i,e)%p - crystallite_dotTemperature(g,i,e) = a(1,3) * RKCK45dotTemperature(1,g,i,e) & - + a(2,3) * RKCK45dotTemperature(2,g,i,e) & - + a(3,3) * RKCK45dotTemperature(3,g,i,e) elseif (n == 4) then constitutive_dotState(g,i,e)%p = a(1,4) * constitutive_RKCK45dotState(1,g,i,e)%p & + a(2,4) * constitutive_RKCK45dotState(2,g,i,e)%p & + a(3,4) * constitutive_RKCK45dotState(3,g,i,e)%p & + a(4,4) * constitutive_RKCK45dotState(4,g,i,e)%p - crystallite_dotTemperature(g,i,e) = a(1,4) * RKCK45dotTemperature(1,g,i,e) & - + a(2,4) * RKCK45dotTemperature(2,g,i,e) & - + a(3,4) * RKCK45dotTemperature(3,g,i,e) & - + a(4,4) * RKCK45dotTemperature(4,g,i,e) elseif (n == 5) then constitutive_dotState(g,i,e)%p = a(1,5) * constitutive_RKCK45dotState(1,g,i,e)%p & + a(2,5) * constitutive_RKCK45dotState(2,g,i,e)%p & + a(3,5) * constitutive_RKCK45dotState(3,g,i,e)%p & + a(4,5) * constitutive_RKCK45dotState(4,g,i,e)%p & + a(5,5) * constitutive_RKCK45dotState(5,g,i,e)%p - crystallite_dotTemperature(g,i,e) = a(1,5) * RKCK45dotTemperature(1,g,i,e) & - + a(2,5) * RKCK45dotTemperature(2,g,i,e) & - + a(3,5) * RKCK45dotTemperature(3,g,i,e) & - + a(4,5) * RKCK45dotTemperature(4,g,i,e) & - + a(5,5) * RKCK45dotTemperature(5,g,i,e) endif endif enddo; enddo; enddo @@ -1779,8 +1722,6 @@ subroutine crystallite_integrateStateRKCK45() mySizeDotState = constitutive_sizeDotState(g,i,e) constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) & + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e) - crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) & - + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) endif enddo; enddo; enddo !$OMP ENDDO @@ -1809,7 +1750,7 @@ subroutine crystallite_integrateStateRKCK45() !$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 if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), & + call constitutive_microstructure(crystallite_Temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states endif enddo; enddo; enddo @@ -1844,10 +1785,8 @@ subroutine crystallite_integrateStateRKCK45() 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 call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & - crystallite_Temperature(g,i,e), c(n)*crystallite_subdt(g,i,e), & ! fraction of original timestep + crystallite_Temperature(i,e), c(n)*crystallite_subdt(g,i,e), & ! fraction of original timestep crystallite_subFrac, g,i,e) - crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & - crystallite_Temperature(g,i,e),g,i,e) endif enddo; enddo; enddo !$OMP ENDDO @@ -1855,8 +1794,7 @@ subroutine crystallite_integrateStateRKCK45() 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 - if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState - .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature + if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p)) then ! NaN occured in 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 @@ -1874,17 +1812,15 @@ subroutine crystallite_integrateStateRKCK45() !-------------------------------------------------------------------------------------------------- - ! --- STATE UPDATE WITH ERROR ESTIMATE FOR STATE AND TEMPERATURE --- + ! --- STATE UPDATE WITH ERROR ESTIMATE FOR STATE --- relStateResiduum = 0.0_pReal - relTemperatureResiduum = 0.0_pReal !$OMP PARALLEL PRIVATE(mySizeDotState) !$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 if (crystallite_todo(g,i,e)) then mySizeDotState = constitutive_sizeDotState(g,i,e) constitutive_RKCK45dotState(6,g,i,e)%p = constitutive_dotState(g,i,e)%p ! store Runge-Kutta dotState - RKCK45dotTemperature(6,g,i,e) = crystallite_dotTemperature(g,i,e) ! store Runge-Kutta dotTemperature endif enddo; enddo; enddo !$OMP ENDDO @@ -1894,7 +1830,7 @@ subroutine crystallite_integrateStateRKCK45() if (crystallite_todo(g,i,e)) then mySizeDotState = constitutive_sizeDotState(g,i,e) - ! --- absolute residuum in state and temperature --- + ! --- absolute residuum in state --- ! NEED TO DO THE ADDITION IN THIS LENGTHY WAY BECAUSE OF PARALLELIZATION ! CAN'T USE A REDUCTION CLAUSE ON A POINTER OR USER DEFINED TYPE @@ -1906,15 +1842,8 @@ subroutine crystallite_integrateStateRKCK45() + db(5) * constitutive_RKCK45dotState(5,g,i,e)%p(1:mySizeDotState) & + db(6) * constitutive_RKCK45dotState(6,g,i,e)%p(1:mySizeDotState)) & * crystallite_subdt(g,i,e) - temperatureResiduum(g,i,e) = ( db(1) * RKCK45dotTemperature(1,g,i,e) & - + db(2) * RKCK45dotTemperature(2,g,i,e) & - + db(3) * RKCK45dotTemperature(3,g,i,e) & - + db(4) * RKCK45dotTemperature(4,g,i,e) & - + db(5) * RKCK45dotTemperature(5,g,i,e) & - + db(6) * RKCK45dotTemperature(6,g,i,e)) & - * crystallite_subdt(g,i,e) - ! --- dot state and dot temperature --- + ! --- dot state --- constitutive_dotState(g,i,e)%p = b(1) * constitutive_RKCK45dotState(1,g,i,e)%p & + b(2) * constitutive_RKCK45dotState(2,g,i,e)%p & @@ -1922,17 +1851,11 @@ subroutine crystallite_integrateStateRKCK45() + b(4) * constitutive_RKCK45dotState(4,g,i,e)%p & + b(5) * constitutive_RKCK45dotState(5,g,i,e)%p & + b(6) * constitutive_RKCK45dotState(6,g,i,e)%p - crystallite_dotTemperature(g,i,e) = b(1) * RKCK45dotTemperature(1,g,i,e) & - + b(2) * RKCK45dotTemperature(2,g,i,e) & - + b(3) * RKCK45dotTemperature(3,g,i,e) & - + b(4) * RKCK45dotTemperature(4,g,i,e) & - + b(5) * RKCK45dotTemperature(5,g,i,e) & - + b(6) * RKCK45dotTemperature(6,g,i,e) endif enddo; enddo; enddo !$OMP ENDDO - ! --- state and temperature update --- + ! --- state and update --- !$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 @@ -1940,8 +1863,6 @@ subroutine crystallite_integrateStateRKCK45() mySizeDotState = constitutive_sizeDotState(g,i,e) constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) & + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e) - crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) & - + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) endif enddo; enddo; enddo !$OMP ENDDO @@ -1954,15 +1875,11 @@ subroutine crystallite_integrateStateRKCK45() mySizeDotState = constitutive_sizeDotState(g,i,e) forall (s = 1_pInt:mySizeDotState, abs(constitutive_state(g,i,e)%p(s)) > 0.0_pReal) & relStateResiduum(s,g,i,e) = stateResiduum(s,g,i,e) / constitutive_state(g,i,e)%p(s) - if (crystallite_Temperature(g,i,e) > 0) then - relTemperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) / crystallite_Temperature(g,i,e) - endif - !$OMP FLUSH(relStateResiduum,relTemperatureResiduum) + !$OMP FLUSH(relStateResiduum) crystallite_todo(g,i,e) = & ( all( abs(relStateResiduum(:,g,i,e)) < rTol_crystalliteState & - .or. abs(stateResiduum(1:mySizeDotState,g,i,e)) < constitutive_aTolState(g,i,e)%p(1:mySizeDotState) ) & - .and. abs(relTemperatureResiduum(g,i,e)) < rTol_crystalliteTemperature ) + .or. abs(stateResiduum(1:mySizeDotState,g,i,e)) < constitutive_aTolState(g,i,e)%p(1:mySizeDotState) )) #ifndef _OPENMP if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt& @@ -2005,7 +1922,7 @@ subroutine crystallite_integrateStateRKCK45() !$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 if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), & + call constitutive_microstructure(crystallite_Temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states endif enddo; enddo; enddo @@ -2067,7 +1984,7 @@ end subroutine crystallite_integrateStateRKCK45 !-------------------------------------------------------------------------------------------------- -!> @brief integrate stress, state and Temperature with 1st order Euler method with adaptive step size +!> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- subroutine crystallite_integrateStateAdaptiveEuler() use debug, only: debug_level, & @@ -2080,7 +1997,6 @@ subroutine crystallite_integrateStateAdaptiveEuler() debug_g, & debug_StateLoopDistribution use numerics, only: rTol_crystalliteState, & - rTol_crystalliteTemperature, & numerics_integrationMode use FEsolving, only: FEsolving_execElem, & FEsolving_execIP @@ -2096,7 +2012,6 @@ subroutine crystallite_integrateStateAdaptiveEuler() constitutive_subState0, & constitutive_dotState, & constitutive_collectDotState, & - constitutive_dotTemperature, & constitutive_microstructure implicit none @@ -2112,9 +2027,6 @@ subroutine crystallite_integrateStateAdaptiveEuler() real(pReal), dimension(constitutive_maxSizeDotState,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & stateResiduum, & ! residuum from evolution in micrstructure relStateResiduum ! relative residuum from evolution in microstructure - real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - temperatureResiduum, & ! residuum from evolution in temperature - relTemperatureResiduum ! relative residuum from evolution in temperature logical singleRun ! flag indicating computation for single (g,i,e) triple @@ -2134,15 +2046,13 @@ subroutine crystallite_integrateStateAdaptiveEuler() if (numerics_integrationMode == 1_pInt) then - ! --- DOT STATE AND TEMPERATURE (EULER INTEGRATION) --- + ! --- DOT STATE (EULER INTEGRATION) --- !$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 if (crystallite_todo(g,i,e)) then call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & - crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) - crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & - crystallite_Temperature(g,i,e),g,i,e) + crystallite_Temperature(i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) endif enddo; enddo; enddo !$OMP ENDDO @@ -2150,8 +2060,7 @@ subroutine crystallite_integrateStateAdaptiveEuler() 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 - if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState - .or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature + if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p)) then ! NaN occured in 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 @@ -2171,12 +2080,9 @@ subroutine crystallite_integrateStateAdaptiveEuler() 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 mySizeDotState = constitutive_sizeDotState(g,i,e) - stateResiduum(1:mySizeDotState,g,i,e) = - 0.5_pReal * constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state and temperature - temperatureResiduum(g,i,e) = - 0.5_pReal * crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) + stateResiduum(1:mySizeDotState,g,i,e) = - 0.5_pReal * constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) & + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e) - crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) & - + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) endif enddo; enddo; enddo !$OMP ENDDO @@ -2205,7 +2111,7 @@ subroutine crystallite_integrateStateAdaptiveEuler() !$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 if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), & + call constitutive_microstructure(crystallite_Temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states endif enddo; enddo; enddo @@ -2233,15 +2139,13 @@ subroutine crystallite_integrateStateAdaptiveEuler() if (numerics_integrationMode == 1_pInt) then - ! --- DOT STATE AND TEMPERATURE (HEUN METHOD) --- + ! --- DOT STATE (HEUN METHOD) --- !$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 if (crystallite_todo(g,i,e)) then call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & - crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) - crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & - crystallite_Temperature(g,i,e),g,i,e) + crystallite_Temperature(i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) endif enddo; enddo; enddo !$OMP ENDDO @@ -2249,8 +2153,7 @@ subroutine crystallite_integrateStateAdaptiveEuler() 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 - if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState - .or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature + if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) ) then ! NaN occured in 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 @@ -2264,11 +2167,10 @@ subroutine crystallite_integrateStateAdaptiveEuler() !$OMP ENDDO - ! --- ERROR ESTIMATE FOR STATE AND TEMPERATURE (HEUN METHOD) --- + ! --- ERROR ESTIMATE FOR STATE (HEUN METHOD) --- !$OMP SINGLE relStateResiduum = 0.0_pReal - relTemperatureResiduum = 0.0_pReal !$OMP END SINGLE !$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 @@ -2279,19 +2181,14 @@ subroutine crystallite_integrateStateAdaptiveEuler() ! --- contribution of heun step to absolute residui --- stateResiduum(1:mySizeDotState,g,i,e) = stateResiduum(1:mySizeDotState,g,i,e) & - + 0.5_pReal * constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state and temperature - temperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) & - + 0.5_pReal * crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) - !$OMP FLUSH(stateResiduum,temperatureResiduum) + + 0.5_pReal * constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state + !$OMP FLUSH(stateResiduum) ! --- relative residui --- forall (s = 1_pInt:mySizeDotState, abs(constitutive_state(g,i,e)%p(s)) > 0.0_pReal) & relStateResiduum(s,g,i,e) = stateResiduum(s,g,i,e) / constitutive_state(g,i,e)%p(s) - if (crystallite_Temperature(g,i,e) > 0_pInt) then - relTemperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) / crystallite_Temperature(g,i,e) - endif - !$OMP FLUSH(relStateResiduum,relTemperatureResiduum) + !$OMP FLUSH(relStateResiduum) #ifndef _OPENMP if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & @@ -2316,8 +2213,7 @@ subroutine crystallite_integrateStateAdaptiveEuler() ! --- converged ? --- if ( all( abs(relStateResiduum(:,g,i,e)) < rTol_crystalliteState & - .or. abs(stateResiduum(1:mySizeDotState,g,i,e)) < constitutive_aTolState(g,i,e)%p(1:mySizeDotState)) & - .and. abs(relTemperatureResiduum(g,i,e)) < rTol_crystalliteTemperature ) then + .or. abs(stateResiduum(1:mySizeDotState,g,i,e)) < constitutive_aTolState(g,i,e)%p(1:mySizeDotState))) then crystallite_converged(g,i,e) = .true. ! ... converged per definitionem if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then !$OMP CRITICAL (distributionState) @@ -2369,7 +2265,7 @@ end subroutine crystallite_integrateStateAdaptiveEuler !-------------------------------------------------------------------------------------------------- -!> @brief integrate stress, state and Temperature with 1st order explicit Euler method +!> @brief integrate stress, and state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- subroutine crystallite_integrateStateEuler() use numerics, only: numerics_integrationMode, & @@ -2393,7 +2289,6 @@ subroutine crystallite_integrateStateEuler() constitutive_subState0, & constitutive_dotState, & constitutive_collectDotState, & - constitutive_dotTemperature, & constitutive_microstructure implicit none @@ -2420,15 +2315,13 @@ eIter = FEsolving_execElem(1:2) if (numerics_integrationMode == 1_pInt) then - ! --- DOT STATE AND TEMPERATURE --- + ! --- DOT STATE --- !$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 if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & - crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) - crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & - crystallite_Temperature(g,i,e),g,i,e) + crystallite_Temperature(i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) endif enddo; enddo; enddo !$OMP ENDDO @@ -2436,9 +2329,8 @@ eIter = FEsolving_execElem(1:2) 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) .and. .not. crystallite_converged(g,i,e)) then - if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState - .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature - if (.not. crystallite_localPlasticity(g,i,e) .and. .not. numerics_timeSyncing) then ! if broken non-local... + if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) ) then ! NaN occured in dotState + if (.not. crystallite_localPlasticity(g,i,e) .and. .not. numerics_timeSyncing) then ! if broken non-local... !$OMP CRITICAL (checkTodo) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped !$OMP END CRITICAL (checkTodo) @@ -2451,7 +2343,7 @@ eIter = FEsolving_execElem(1:2) !$OMP ENDDO - ! --- UPDATE STATE AND TEMPERATURE --- + ! --- UPDATE STATE --- !$OMP DO PRIVATE(mySizeDotState) 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 @@ -2459,8 +2351,6 @@ eIter = FEsolving_execElem(1:2) mySizeDotState = constitutive_sizeDotState(g,i,e) constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) & + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e) - crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) & - + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) #ifndef _OPENMP if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) & @@ -2502,7 +2392,7 @@ eIter = FEsolving_execElem(1:2) !$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 if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), & + call constitutive_microstructure(crystallite_Temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states endif enddo; enddo; enddo @@ -2562,7 +2452,7 @@ end subroutine crystallite_integrateStateEuler !-------------------------------------------------------------------------------------------------- -!> @brief integrate stress, state and Temperature with adaptive 1st order explicit Euler method +!> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- subroutine crystallite_integrateStateFPI() @@ -2577,8 +2467,7 @@ subroutine crystallite_integrateStateFPI() debug_StateLoopDistribution use numerics, only: nState, & numerics_integrationMode, & - rTol_crystalliteState, & - rTol_crystalliteTemperature + rTol_crystalliteState use FEsolving, only: FEsolving_execElem, & FEsolving_execIP use mesh, only: mesh_element, & @@ -2590,7 +2479,6 @@ subroutine crystallite_integrateStateFPI() constitutive_maxSizeDotState, & constitutive_dotState, & constitutive_collectDotState, & - constitutive_dotTemperature, & constitutive_microstructure, & constitutive_previousDotState, & constitutive_previousDotState2, & @@ -2609,8 +2497,7 @@ subroutine crystallite_integrateStateFPI() gIter ! bounds for grain iteration real(pReal) dot_prod12, & dot_prod22, & - stateDamper, & ! damper for integration of state - temperatureResiduum + stateDamper ! damper for integration of state real(pReal), dimension(constitutive_maxSizeDotState) :: & stateResiduum, & tempState @@ -2643,9 +2530,7 @@ subroutine crystallite_integrateStateFPI() 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 !< @ToDo: Put in loop above? call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & - crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) - crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & - crystallite_Temperature(g,i,e),g,i,e) + crystallite_Temperature(i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) endif enddo; enddo; enddo !$OMP ENDDO @@ -2653,8 +2538,7 @@ subroutine crystallite_integrateStateFPI() 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 - if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState - .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature + if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) ) then ! NaN occured in 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 !< @ToDo: no (g,i,e) index here? ...all non-locals done (and broken) @@ -2668,7 +2552,7 @@ subroutine crystallite_integrateStateFPI() !$OMP ENDDO - ! --- UPDATE STATE AND TEMPERATURE --- + ! --- UPDATE STATE --- !$OMP DO PRIVATE(mySizeDotState) 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 @@ -2676,8 +2560,6 @@ subroutine crystallite_integrateStateFPI() mySizeDotState = constitutive_sizeDotState(g,i,e) constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) & + constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e) - crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) & - + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) endif enddo; enddo; enddo !$OMP ENDDO @@ -2699,7 +2581,7 @@ subroutine crystallite_integrateStateFPI() !$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 if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), & + call constitutive_microstructure(crystallite_Temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states endif constitutive_previousDotState2(g,i,e)%p = constitutive_previousDotState(g,i,e)%p ! remember previous dotState @@ -2735,15 +2617,13 @@ subroutine crystallite_integrateStateFPI() !$OMP END SINGLE - ! --- DOT STATE AND TEMPERATURE --- + ! --- DOT STATE --- !$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 if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & - crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) - crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & - crystallite_Temperature(g,i,e),g,i,e) + crystallite_Temperature(i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) endif enddo; enddo; enddo !$OMP ENDDO @@ -2751,8 +2631,7 @@ subroutine crystallite_integrateStateFPI() 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) .and. .not. crystallite_converged(g,i,e)) then - if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState - .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature + if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) ) then ! NaN occured in dotState crystallite_todo(g,i,e) = .false. ! ... skip me next time if (.not. crystallite_localPlasticity(g,i,e)) then ! if me is non-local... !$OMP CRITICAL (checkTodo) @@ -2765,9 +2644,9 @@ subroutine crystallite_integrateStateFPI() !$OMP ENDDO - ! --- UPDATE STATE AND TEMPERATURE --- + ! --- UPDATE STATE --- - !$OMP DO PRIVATE(dot_prod12,dot_prod22,statedamper,mySizeDotState,stateResiduum,temperatureResiduum,tempState) + !$OMP DO PRIVATE(dot_prod12,dot_prod22,statedamper,mySizeDotState,stateResiduum,tempState) 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) .and. .not. crystallite_converged(g,i,e)) then @@ -2792,15 +2671,10 @@ subroutine crystallite_integrateStateFPI() - constitutive_subState0(g,i,e)%p(1:mySizeDotState) & - (constitutive_dotState(g,i,e)%p * statedamper & + constitutive_previousDotState(g,i,e)%p * (1.0_pReal - statedamper)) * crystallite_subdt(g,i,e) - temperatureResiduum = crystallite_Temperature(g,i,e) & - - crystallite_subTemperature0(g,i,e) & - - crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) - ! --- correct state and temperature with residuum --- - + ! --- correct state with residuum --- tempState(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) - stateResiduum(1:mySizeDotState) ! need to copy to local variable, since we cant flush a pointer in openmp - crystallite_Temperature(g,i,e) = crystallite_Temperature(g,i,e) - temperatureResiduum - !$OMP FLUSH(crystallite_Temperature) + #ifndef _OPENMP if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) & @@ -2822,9 +2696,7 @@ subroutine crystallite_integrateStateFPI() if ( all( abs(stateResiduum(1:mySizeDotState)) < constitutive_aTolState(g,i,e)%p(1:mySizeDotState) & .or. abs(stateResiduum(1:mySizeDotState)) < rTol_crystalliteState & - * abs(tempState(1:mySizeDotState)) ) & - .and. (abs(temperatureResiduum) < rTol_crystalliteTemperature * crystallite_Temperature(g,i,e) & - .or. crystallite_Temperature(g,i,e) == 0.0_pReal) ) then + * abs(tempState(1:mySizeDotState)) ) ) then crystallite_converged(g,i,e) = .true. ! ... converged per definitionem if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then !$OMP CRITICAL (distributionState) @@ -2925,12 +2797,11 @@ logical function crystallite_stateJump(g,i,e) crystallite_stateJump = .false. - call constitutive_collectDeltaState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Temperature(g,i,e), g,i,e) + call constitutive_collectDeltaState(crystallite_Tstar_v(1:6,g,i,e), g,i,e) mySizeDotState = constitutive_sizeDotState(g,i,e) - if (any(constitutive_deltaState(g,i,e)%p(1:mySizeDotState) /= constitutive_deltaState(g,i,e)%p(1:mySizeDotState))) then + if (any(constitutive_deltaState(g,i,e)%p(1:mySizeDotState) /= constitutive_deltaState(g,i,e)%p(1:mySizeDotState))) & return - endif constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) & + constitutive_deltaState(g,i,e)%p(1:mySizeDotState) @@ -3145,7 +3016,7 @@ logical function crystallite_integrateStress(& if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) & call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) - call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT_constitutive, Tstar_v, crystallite_Temperature(g,i,e), g, i, e) + call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT_constitutive, Tstar_v, crystallite_Temperature(i,e), g, i, e) if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then call system_clock(count=tock,count_rate=tickrate,count_max=maxticks) @@ -3432,7 +3303,7 @@ end subroutine crystallite_orientations !-------------------------------------------------------------------------------------------------- !> @brief return results of particular grain !-------------------------------------------------------------------------------------------------- -function crystallite_postResults(dt, gr, ip, el) +function crystallite_postResults(ipc, ip, el) use math, only: & math_qToEuler, & math_qToEulerAxisAngle, & @@ -3467,16 +3338,15 @@ function crystallite_postResults(dt, gr, ip, el) implicit none integer(pInt), intent(in):: el, & !< element index ip, & !< integration point index - gr !< grain index - real(pReal), intent(in):: dt !< time increment + ipc !< grain index real(pReal), dimension(1+crystallite_sizePostResults(microstructure_crystallite(mesh_element(4,el)))+ & - 1+constitutive_sizePostResults(gr,ip,el)) :: crystallite_postResults + 1+constitutive_sizePostResults(ipc,ip,el)) :: crystallite_postResults real(pReal), dimension(3,3) :: Ee real(pReal), dimension(4) :: rotation - real(pReal) detF - integer(pInt) o,c,crystID,mySize,n + real(pReal) :: detF + integer(pInt) :: o,c,crystID,mySize,n crystID = microstructure_crystallite(mesh_element(4,el)) @@ -3490,38 +3360,38 @@ function crystallite_postResults(dt, gr, ip, el) select case(crystallite_output(o,crystID)) case ('phase') mySize = 1_pInt - crystallite_postResults(c+1) = real(material_phase(gr,ip,el),pReal) ! phaseID of grain + crystallite_postResults(c+1) = real(material_phase(ipc,ip,el),pReal) ! phaseID of grain case ('texture') mySize = 1_pInt - crystallite_postResults(c+1) = real(material_texture(gr,ip,el),pReal) ! textureID of grain + crystallite_postResults(c+1) = real(material_texture(ipc,ip,el),pReal) ! textureID of grain case ('volume') mySize = 1_pInt - detF = math_det33(crystallite_partionedF(1:3,1:3,gr,ip,el)) ! V_current = det(F) * V_reference + detF = math_det33(crystallite_partionedF(1:3,1:3,ipc,ip,el)) ! V_current = det(F) * V_reference crystallite_postResults(c+1) = detF * mesh_ipVolume(ip,el) / & homogenization_Ngrains(mesh_element(3,el)) ! grain volume (not fraction but absolute) case ('orientation') mySize = 4_pInt - crystallite_postResults(c+1:c+mySize) = crystallite_orientation(1:4,gr,ip,el) ! grain orientation as quaternion + crystallite_postResults(c+1:c+mySize) = crystallite_orientation(1:4,ipc,ip,el) ! grain orientation as quaternion case ('eulerangles') mySize = 3_pInt crystallite_postResults(c+1:c+mySize) = inDeg * & - math_qToEuler(crystallite_orientation(1:4,gr,ip,el)) ! grain orientation as Euler angles in degree + math_qToEuler(crystallite_orientation(1:4,ipc,ip,el)) ! grain orientation as Euler angles in degree case ('grainrotation') mySize = 4_pInt crystallite_postResults(c+1:c+mySize) = & - math_qToEulerAxisAngle(crystallite_rotation(1:4,gr,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates + math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates crystallite_postResults(c+4) = inDeg * crystallite_postResults(c+4) ! angle in degree case ('grainrotationx') mySize = 1_pInt - rotation = math_qToEulerAxisAngle(crystallite_rotation(1:4,gr,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates + rotation = math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates crystallite_postResults(c+1) = inDeg * rotation(1) * rotation(4) ! angle in degree case ('grainrotationy') mySize = 1_pInt - rotation = math_qToEulerAxisAngle(crystallite_rotation(1:4,gr,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates + rotation = math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates crystallite_postResults(c+1) = inDeg * rotation(2) * rotation(4) ! angle in degree case ('grainrotationz') mySize = 1_pInt - rotation = math_qToEulerAxisAngle(crystallite_rotation(1:4,gr,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates + rotation = math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates crystallite_postResults(c+1) = inDeg * rotation(3) * rotation(4) ! angle in degree case ('ipcoords') mySize = 3_pInt @@ -3533,40 +3403,40 @@ function crystallite_postResults(dt, gr, ip, el) case ('defgrad','f') mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = & - reshape(math_transpose33(crystallite_partionedF(1:3,1:3,gr,ip,el)),[mySize]) + reshape(math_transpose33(crystallite_partionedF(1:3,1:3,ipc,ip,el)),[mySize]) case ('e') mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = 0.5_pReal * reshape((math_mul33x33( & - math_transpose33(crystallite_partionedF(1:3,1:3,gr,ip,el)), & - crystallite_partionedF(1:3,1:3,gr,ip,el)) - math_I3),[mySize]) + math_transpose33(crystallite_partionedF(1:3,1:3,ipc,ip,el)), & + crystallite_partionedF(1:3,1:3,ipc,ip,el)) - math_I3),[mySize]) case ('fe') mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = & - reshape(math_transpose33(crystallite_Fe(1:3,1:3,gr,ip,el)),[mySize]) + reshape(math_transpose33(crystallite_Fe(1:3,1:3,ipc,ip,el)),[mySize]) case ('ee') - Ee = 0.5_pReal *(math_mul33x33(math_transpose33(crystallite_Fe(1:3,1:3,gr,ip,el)), & - crystallite_Fe(1:3,1:3,gr,ip,el)) - math_I3) + Ee = 0.5_pReal *(math_mul33x33(math_transpose33(crystallite_Fe(1:3,1:3,ipc,ip,el)), & + crystallite_Fe(1:3,1:3,ipc,ip,el)) - math_I3) mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = reshape(Ee,[mySize]) case ('fp') mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = & - reshape(math_transpose33(crystallite_Fp(1:3,1:3,gr,ip,el)),[mySize]) + reshape(math_transpose33(crystallite_Fp(1:3,1:3,ipc,ip,el)),[mySize]) case ('lp') mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = & - reshape(math_transpose33(crystallite_Lp(1:3,1:3,gr,ip,el)),[mySize]) + reshape(math_transpose33(crystallite_Lp(1:3,1:3,ipc,ip,el)),[mySize]) case ('p','firstpiola','1stpiola') mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = & - reshape(math_transpose33(crystallite_P(1:3,1:3,gr,ip,el)),[mySize]) + reshape(math_transpose33(crystallite_P(1:3,1:3,ipc,ip,el)),[mySize]) case ('s','tstar','secondpiola','2ndpiola') mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = & - reshape(math_Mandel6to33(crystallite_Tstar_v(1:6,gr,ip,el)),[mySize]) + reshape(math_Mandel6to33(crystallite_Tstar_v(1:6,ipc,ip,el)),[mySize]) case ('elasmatrix') mySize = 36_pInt - crystallite_postResults(c+1:c+mySize) = reshape(constitutive_homogenizedC(gr,ip,el),[mySize]) + crystallite_postResults(c+1:c+mySize) = reshape(constitutive_homogenizedC(ipc,ip,el),[mySize]) case('neighboringelement') mySize = mesh_maxNipNeighbors crystallite_postResults(c+1:c+mySize) = 0.0_pReal @@ -3581,13 +3451,13 @@ function crystallite_postResults(dt, gr, ip, el) c = c + mySize enddo - crystallite_postResults(c+1) = real(constitutive_sizePostResults(gr,ip,el),pReal) ! size of constitutive results + crystallite_postResults(c+1) = real(constitutive_sizePostResults(ipc,ip,el),pReal) ! size of constitutive results c = c + 1_pInt - if (constitutive_sizePostResults(gr,ip,el) > 0_pInt) & - crystallite_postResults(c+1:c+constitutive_sizePostResults(gr,ip,el)) = & - constitutive_postResults(crystallite_Tstar_v(1:6,gr,ip,el), crystallite_Fe, & - crystallite_Temperature(gr,ip,el), dt, gr, ip, el) - c = c + constitutive_sizePostResults(gr,ip,el) + if (constitutive_sizePostResults(ipc,ip,el) > 0_pInt) & + crystallite_postResults(c+1:c+constitutive_sizePostResults(ipc,ip,el)) = & + constitutive_postResults(crystallite_Tstar_v(1:6,ipc,ip,el), crystallite_Fe, & + crystallite_Temperature(ip,el), ipc, ip, el) + c = c + constitutive_sizePostResults(ipc,ip,el) end function crystallite_postResults diff --git a/code/debug.f90 b/code/debug.f90 index d77673d42..f8be27a68 100644 --- a/code/debug.f90 +++ b/code/debug.f90 @@ -61,25 +61,23 @@ module debug debug_MARC = 12_pInt, & debug_ABAQUS = 13_pInt integer(pInt), parameter, private :: & - debug_MAXNTYPE = debug_ABAQUS ! must be set to the maximum defined debug type + debug_MAXNTYPE = debug_ABAQUS !< must be set to the maximum defined debug type integer(pInt),protected, dimension(debug_maxNtype+2_pInt), public :: & ! specific ones, and 2 for "all" and "other" debug_level = 0_pInt integer(pInt), public :: & - debug_cumLpCalls = 0_pInt, & ! total number of calls to LpAndItsTangent - debug_cumDeltaStateCalls = 0_pInt, & ! total number of calls to deltaState - debug_cumDotStateCalls = 0_pInt, & ! total number of calls to dotState - debug_cumDotTemperatureCalls = 0_pInt, & ! total number of calls to dotTemprature + debug_cumLpCalls = 0_pInt, & !< total number of calls to LpAndItsTangent + debug_cumDeltaStateCalls = 0_pInt, & !< total number of calls to deltaState + debug_cumDotStateCalls = 0_pInt, & !< total number of calls to dotState debug_e = 1_pInt, & debug_i = 1_pInt, & debug_g = 1_pInt integer(pLongInt), public :: & - debug_cumLpTicks = 0_pLongInt, & ! total cpu ticks spent in LpAndItsTangent - debug_cumDeltaStateTicks = 0_pLongInt, & ! total cpu ticks spent in deltaState - debug_cumDotStateTicks = 0_pLongInt, & ! total cpu ticks spent in dotState - debug_cumDotTemperatureTicks = 0_pLongInt ! total cpu ticks spent in dotTemperature + debug_cumLpTicks = 0_pLongInt, & !< total cpu ticks spent in LpAndItsTangent + debug_cumDeltaStateTicks = 0_pLongInt, & !< total cpu ticks spent in deltaState + debug_cumDotStateTicks = 0_pLongInt !< total cpu ticks spent in dotState integer(pInt), dimension(2), public :: & debug_stressMaxLocation = 0_pInt, & @@ -88,13 +86,13 @@ module debug debug_jacobianMinLocation = 0_pInt integer(pInt), dimension(:), allocatable, public :: & - debug_CrystalliteLoopDistribution, & ! distribution of crystallite cutbacks + debug_CrystalliteLoopDistribution, & !< distribution of crystallite cutbacks debug_MaterialpointStateLoopDistribution, & debug_MaterialpointLoopDistribution integer(pInt), dimension(:,:), allocatable, public :: & - debug_StressLoopDistribution, & ! distribution of stress iterations until convergence - debug_StateLoopDistribution ! distribution of state iterations until convergence + debug_StressLoopDistribution, & !< distribution of stress iterations until convergence + debug_StateLoopDistribution !< distribution of state iterations until convergence real(pReal), public :: & debug_stressMax = -huge(1.0_pReal), & @@ -103,7 +101,7 @@ module debug debug_jacobianMin = huge(1.0_pReal) character(len=64), parameter, private :: & - debug_CONFIGFILE = 'debug.config' ! name of configuration file + debug_CONFIGFILE = 'debug.config' !< name of configuration file #ifdef PETSc character(len=1024), parameter, public :: & @@ -335,11 +333,9 @@ subroutine debug_reset debug_cumLpTicks = 0_pLongInt debug_cumDeltaStateTicks = 0_pLongInt debug_cumDotStateTicks = 0_pLongInt - debug_cumDotTemperatureTicks = 0_pLongInt debug_cumLpCalls = 0_pInt debug_cumDeltaStateCalls = 0_pInt debug_cumDotStateCalls = 0_pInt - debug_cumDotTemperatureCalls = 0_pInt debug_stressMaxLocation = 0_pInt debug_stressMinLocation = 0_pInt debug_jacobianMaxLocation = 0_pInt @@ -394,13 +390,6 @@ subroutine debug_info write(6,'(a33,1x,f12.6)') 'avg CPU time/microsecs per call :',& real(debug_cumDeltaStateTicks,pReal)*1.0e6_pReal/real(tickrate*debug_cumDeltaStateCalls,pReal) endif - write(6,'(/,a33,1x,i12)') 'total calls to dotTemperature :',debug_cumDotTemperatureCalls - if (debug_cumdotTemperatureCalls > 0_pInt) then - write(6,'(a33,1x,f12.3)') 'total CPU time/s :',& - real(debug_cumDotTemperatureTicks,pReal)/real(tickrate,pReal) - write(6,'(a33,1x,f12.6)') 'avg CPU time/microsecs per call :',& - real(debug_cumDotTemperatureTicks,pReal)*1.0e6_pReal/real(tickrate*debug_cumDotTemperatureCalls,pReal) - endif integral = 0_pInt write(6,'(3/,a)') 'distribution_StressLoop : stress stiffness' diff --git a/code/homogenization.f90 b/code/homogenization.f90 index 95be41da4..caa5cfe4b 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -36,8 +36,6 @@ module homogenization private type(p_vec), dimension(:,:), allocatable, public :: & homogenization_state0 !< pointer array to homogenization state at start of FE increment - real(pReal), dimension(:,:), allocatable, public :: & - materialpoint_Temperature !< temperature at IP real(pReal), dimension(:,:,:,:), allocatable, public :: & materialpoint_F0, & !< def grad of IP at start of FE increment materialpoint_F, & !< def grad of IP to be reached at end of FE increment @@ -63,7 +61,8 @@ module homogenization real(pReal), dimension(:,:), allocatable, private :: & materialpoint_subFrac, & materialpoint_subStep, & - materialpoint_subdt + materialpoint_subdt, & + materialpoint_heat integer(pInt), dimension(:,:), allocatable, private :: & homogenization_sizePostResults !< size of postResults array per material point integer(pInt), private :: & @@ -82,7 +81,7 @@ module homogenization homogenization_partitionDeformation, & homogenization_updateState, & homogenization_averageStressAndItsTangent, & - homogenization_averageTemperature, & + homogenization_averageHeat, & homogenization_postResults contains @@ -91,14 +90,16 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief module initialization !-------------------------------------------------------------------------------------------------- -subroutine homogenization_init(Temperature) +subroutine homogenization_init() use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use math, only: & math_I3 use debug, only: & debug_level, & debug_homogenization, & - debug_levelBasic + debug_levelBasic, & + debug_e, & + debug_g use IO, only: & IO_error, & IO_open_file, & @@ -121,18 +122,16 @@ subroutine homogenization_init(Temperature) use homogenization_RGC implicit none - real(pReal) Temperature integer(pInt), parameter :: fileunit = 200_pInt - integer(pInt) e,i,p,myInstance + integer(pInt) :: e,i,p,myInstance integer(pInt), dimension(:,:), pointer :: thisSize character(len=64), dimension(:,:), pointer :: thisOutput logical :: knownHomogenization !-------------------------------------------------------------------------------------------------- ! parse homogenization from config file - if (.not. IO_open_jobFile_stat(fileunit,material_localFileExt)) then ! no local material configuration present... + if (.not. IO_open_jobFile_stat(fileunit,material_localFileExt)) & ! no local material configuration present... call IO_open_file(fileunit,material_configFile) ! ... open material.config file - endif call homogenization_isostrain_init(fileunit) call homogenization_RGC_init(fileunit) close(fileunit) @@ -173,7 +172,8 @@ subroutine homogenization_init(Temperature) homogenization_sizeState = 0_pInt allocate(homogenization_sizePostResults(mesh_maxNips,mesh_NcpElems)) homogenization_sizePostResults = 0_pInt - + allocate(materialpoint_heat(mesh_maxNips,mesh_NcpElems)) + materialpoint_heat = 0.0_pReal allocate(materialpoint_dPdF(3,3,3,3,mesh_maxNips,mesh_NcpElems)) materialpoint_dPdF = 0.0_pReal allocate(materialpoint_F0(3,3,mesh_maxNips,mesh_NcpElems)) @@ -185,8 +185,6 @@ subroutine homogenization_init(Temperature) materialpoint_subF = 0.0_pReal allocate(materialpoint_P(3,3,mesh_maxNips,mesh_NcpElems)) materialpoint_P = 0.0_pReal - allocate(materialpoint_Temperature(mesh_maxNips,mesh_NcpElems)) - materialpoint_Temperature = Temperature allocate(materialpoint_subFrac(mesh_maxNips,mesh_NcpElems)) materialpoint_subFrac = 0.0_pReal allocate(materialpoint_subStep(mesh_maxNips,mesh_NcpElems)) @@ -204,7 +202,7 @@ subroutine homogenization_init(Temperature) materialpoint_F = materialpoint_F0 !-------------------------------------------------------------------------------------------------- -! allocate and initialize global state and postrestuls variables +! allocate and initialize global state and postresutls variables elementLooping: do e = 1,mesh_NcpElems myInstance = homogenization_typeInstance(mesh_element(3,e)) IpLooping: do i = 1,FE_Nips(FE_geomtype(mesh_element(2,e))) @@ -247,38 +245,38 @@ subroutine homogenization_init(Temperature) + 1 + constitutive_maxSizePostResults) ! constitutive size & constitutive results allocate(materialpoint_results(materialpoint_sizeResults,mesh_maxNips,mesh_NcpElems)) - write(6,'(/,a)') ' <<<+- homogenization init -+>>>' - write(6,'(a)') ' $Id$' - write(6,'(a16,a)') ' Current time : ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- homogenization init -+>>>' + write(6,'(a)') ' $Id$' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then - write(6,'(a32,1x,7(i8,1x))') 'homogenization_state0: ', shape(homogenization_state0) - write(6,'(a32,1x,7(i8,1x))') 'homogenization_subState0: ', shape(homogenization_subState0) - write(6,'(a32,1x,7(i8,1x))') 'homogenization_state: ', shape(homogenization_state) - write(6,'(a32,1x,7(i8,1x))') 'homogenization_sizeState: ', shape(homogenization_sizeState) - write(6,'(a32,1x,7(i8,1x))') 'homogenization_sizePostResults: ', shape(homogenization_sizePostResults) - write(6,*) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F0: ', shape(materialpoint_F0) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F: ', shape(materialpoint_F) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subF0: ', shape(materialpoint_subF0) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subF: ', shape(materialpoint_subF) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_P: ', shape(materialpoint_P) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_Temperature: ', shape(materialpoint_Temperature) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subFrac: ', shape(materialpoint_subFrac) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subStep: ', shape(materialpoint_subStep) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subdt: ', shape(materialpoint_subdt) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_requested: ', shape(materialpoint_requested) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_converged: ', shape(materialpoint_converged) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_doneAndHappy: ', shape(materialpoint_doneAndHappy) - write(6,*) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_results: ', shape(materialpoint_results) - write(6,*) - write(6,'(a32,1x,7(i8,1x))') 'maxSizeState: ', homogenization_maxSizeState - write(6,'(a32,1x,7(i8,1x))') 'maxSizePostResults: ', homogenization_maxSizePostResults + write(6,'(a32,1x,7(i8,1x))') 'homogenization_state0: ', shape(homogenization_state0) + write(6,'(a32,1x,7(i8,1x))') 'homogenization_subState0: ', shape(homogenization_subState0) + write(6,'(a32,1x,7(i8,1x))') 'homogenization_state: ', shape(homogenization_state) + write(6,'(a32,1x,7(i8,1x))') 'homogenization_sizeState: ', shape(homogenization_sizeState) + write(6,'(a32,1x,7(i8,1x),/)') 'homogenization_sizePostResults: ', shape(homogenization_sizePostResults) + write(6,'(a32,1x,7(i8,1x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF) + write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F0: ', shape(materialpoint_F0) + write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F: ', shape(materialpoint_F) + write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subF0: ', shape(materialpoint_subF0) + write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subF: ', shape(materialpoint_subF) + write(6,'(a32,1x,7(i8,1x))') 'materialpoint_P: ', shape(materialpoint_P) + write(6,'(a32,1x,7(i8,1x))') 'materialpoint_heat: ', shape(materialpoint_heat) + write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subFrac: ', shape(materialpoint_subFrac) + write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subStep: ', shape(materialpoint_subStep) + write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subdt: ', shape(materialpoint_subdt) + write(6,'(a32,1x,7(i8,1x))') 'materialpoint_requested: ', shape(materialpoint_requested) + write(6,'(a32,1x,7(i8,1x))') 'materialpoint_converged: ', shape(materialpoint_converged) + write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_doneAndHappy: ', shape(materialpoint_doneAndHappy) + write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_results: ', shape(materialpoint_results) + write(6,'(a32,1x,7(i8,1x))') 'maxSizeState: ', homogenization_maxSizeState + write(6,'(a32,1x,7(i8,1x))') 'maxSizePostResults: ', homogenization_maxSizePostResults endif flush(6) - + + if (debug_g < 1 .or. debug_g > homogenization_Ngrains(mesh_element(3,debug_e))) & + call IO_error(602_pInt,ext_msg='component (grain)') + end subroutine homogenization_init @@ -309,7 +307,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) constitutive_partionedState0, & constitutive_state use crystallite, only: & - crystallite_Temperature, & + crystallite_heat, & crystallite_F0, & crystallite_Fp0, & crystallite_Fp, & @@ -319,7 +317,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) crystallite_dPdF0, & crystallite_Tstar0_v, & crystallite_Tstar_v, & - crystallite_partionedTemperature0, & crystallite_partionedF0, & crystallite_partionedF, & crystallite_partionedFp0, & @@ -357,13 +354,10 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) !-------------------------------------------------------------------------------------------------- ! initialize to starting condition - if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt .and. & - debug_e > 0 .and. debug_e <= mesh_NcpElems .and. debug_i > 0 .and. debug_i <= mesh_maxNips) then + if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then !$OMP CRITICAL (write2out) - write(6,*) - write(6,'(a,i5,1x,i2)') '<< HOMOG >> Material Point start at el ip ', debug_e, debug_i - write(6,'(a,/,12x,f14.9)') '<< HOMOG >> Temp0', & - materialpoint_Temperature(debug_i,debug_e) + write(6,'(/a,i5,1x,i2)') '<< HOMOG >> Material Point start at el ip ', debug_e, debug_i + write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F0', & math_transpose33(materialpoint_F0(1:3,1:3,debug_i,debug_e)) write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F', & @@ -377,7 +371,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) myNgrains = homogenization_Ngrains(mesh_element(3,e)) forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains) constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p ! ...microstructures - crystallite_partionedTemperature0(g,i,e) = materialpoint_Temperature(i,e) ! ...temperatures crystallite_partionedFp0(1:3,1:3,g,i,e) = crystallite_Fp0(1:3,1:3,g,i,e) ! ...plastic def grads crystallite_partionedLp0(1:3,1:3,g,i,e) = crystallite_Lp0(1:3,1:3,g,i,e) ! ...plastic velocity grads crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,g,i,e) = crystallite_dPdF0(1:3,1:3,1:3,1:3,g,i,e) ! ...stiffness @@ -427,7 +420,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) steppingNeeded: if (materialpoint_subStep(i,e) > subStepMinHomog) then ! wind forward grain starting point of... - crystallite_partionedTemperature0(1:myNgrains,i,e) = crystallite_Temperature(1:myNgrains,i,e) ! ...temperatures crystallite_partionedF0(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedF(1:3,1:3,1:myNgrains,i,e) ! ...def grads crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Fp(1:3,1:3,1:myNgrains,i,e) ! ...plastic def grads crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Lp(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads @@ -476,8 +468,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) !-------------------------------------------------------------------------------------------------- ! restore... - crystallite_Temperature(1:myNgrains,i,e) = crystallite_partionedTemperature0(1:myNgrains,i,e) ! ...temperatures - ! ...initial def grad unchanged crystallite_Fp(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic def grads crystallite_Lp(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads crystallite_dPdF(1:3,1:3,1:3,1:3,1:myNgrains,i,e) = crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,1:myNgrains,i,e) ! ...stiffness @@ -575,7 +565,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) elementLooping4: do e = FEsolving_execElem(1),FEsolving_execElem(2) IpLooping4: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) call homogenization_averageStressAndItsTangent(i,e) - materialpoint_Temperature(i,e) = homogenization_averageTemperature(i,e) + materialpoint_heat(i,e) = homogenization_averageHeat(i,e) enddo IpLooping4 enddo elementLooping4 !$OMP END PARALLEL DO @@ -639,7 +629,7 @@ subroutine materialpoint_postResults(dt) grainLooping :do g = 1,myNgrains theSize = (1 + crystallite_sizePostResults(myCrystallite)) + (1 + constitutive_sizePostResults(g,i,e)) - materialpoint_results(thePos+1:thePos+theSize,i,e) = crystallite_postResults(dt,g,i,e) ! tell crystallite results + materialpoint_results(thePos+1:thePos+theSize,i,e) = crystallite_postResults(g,i,e) ! tell crystallite results thePos = thePos + theSize enddo grainLooping enddo IpLooping @@ -677,18 +667,15 @@ subroutine homogenization_partitionDeformation(ip,el) case (homogenization_isostrain_label) chosenHomogenization call homogenization_isostrain_partitionDeformation(& crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), & - crystallite_partionedF0(1:3,1:3,1:homogenization_maxNgrains,ip,el),& materialpoint_subF(1:3,1:3,ip,el),& - homogenization_state(ip,el), & - ip, & el) case (homogenization_RGC_label) chosenHomogenization - call homogenization_RGC_partitionDeformation(crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), & - crystallite_partionedF0(1:3,1:3,1:homogenization_maxNgrains,ip,el),& - materialpoint_subF(1:3,1:3,ip,el),& - homogenization_state(ip,el), & - ip, & - el) + call homogenization_RGC_partitionDeformation(& + crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), & + materialpoint_subF(1:3,1:3,ip,el),& + homogenization_state(ip,el), & + ip, & + el) end select chosenHomogenization end subroutine homogenization_partitionDeformation @@ -765,35 +752,35 @@ subroutine homogenization_averageStressAndItsTangent(ip,el) chosenHomogenization: select case(homogenization_type(mesh_element(3,el))) case (homogenization_isostrain_label) chosenHomogenization - call homogenization_isostrain_averageStressAndItsTangent(materialpoint_P(1:3,1:3,ip,el), & - materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& - crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), & - crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), & - ip, & - el) + call homogenization_isostrain_averageStressAndItsTangent(& + materialpoint_P(1:3,1:3,ip,el), & + materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& + crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), & + crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), & + el) case (homogenization_RGC_label) chosenHomogenization - call homogenization_RGC_averageStressAndItsTangent( materialpoint_P(1:3,1:3,ip,el), & - materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& - crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), & - crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), & - ip, & - el) + call homogenization_RGC_averageStressAndItsTangent(& + materialpoint_P(1:3,1:3,ip,el), & + materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& + crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), & + crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), & + el) end select chosenHomogenization end subroutine homogenization_averageStressAndItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief derive average temperature from constituent quantities (does not depend on choosen +!> @brief derive average heat from constituent quantities (does not depend on choosen !! homogenization scheme) !-------------------------------------------------------------------------------------------------- -real(pReal) function homogenization_averageTemperature(ip,el) +real(pReal) function homogenization_averageHeat(ip,el) use mesh, only: & mesh_element use material, only: & homogenization_Ngrains use crystallite, only: & - crystallite_Temperature + crystallite_heat implicit none integer(pInt), intent(in) :: & @@ -803,11 +790,11 @@ real(pReal) function homogenization_averageTemperature(ip,el) Ngrains !-------------------------------------------------------------------------------------------------- -! computing the average temperature +! computing the average heat Ngrains = homogenization_Ngrains(mesh_element(3,el)) - homogenization_averageTemperature= sum(crystallite_Temperature(1:Ngrains,ip,el))/real(Ngrains,pReal) + homogenization_averageHeat= sum(crystallite_heat(1:Ngrains,ip,el))/real(Ngrains,pReal) -end function homogenization_averageTemperature +end function homogenization_averageHeat !-------------------------------------------------------------------------------------------------- @@ -835,9 +822,9 @@ function homogenization_postResults(ip,el) homogenization_postResults = 0.0_pReal chosenHomogenization: select case (homogenization_type(mesh_element(3,el))) case (homogenization_isostrain_label) chosenHomogenization - homogenization_postResults = homogenization_isostrain_postResults(homogenization_state(ip,el),ip,el) + homogenization_postResults = homogenization_isostrain_postResults(el) case (homogenization_RGC_label) chosenHomogenization - homogenization_postResults = homogenization_RGC_postResults(homogenization_state(ip,el),ip,el) + homogenization_postResults = homogenization_RGC_postResults(homogenization_state(ip,el),el) end select chosenHomogenization end function homogenization_postResults diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index 75c057f29..bdf68a63f 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -77,7 +77,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -subroutine homogenization_RGC_init(myFile) +subroutine homogenization_RGC_init(myUnit) 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, & @@ -101,7 +101,7 @@ subroutine homogenization_RGC_init(myFile) use material implicit none - integer(pInt), intent(in) :: myFile !< file pointer to material configuration + integer(pInt), intent(in) :: myUnit !< file pointer to material configuration integer(pInt), parameter :: MAXNCHUNKS = 4_pInt integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions integer(pInt) ::section=0_pInt, maxNinstance, i,j,e, output=-1_pInt, mySize, myInstance @@ -130,14 +130,14 @@ subroutine homogenization_RGC_init(myFile) allocate(homogenization_RGC_orientation(3,3,mesh_maxNips,mesh_NcpElems)) homogenization_RGC_orientation = spread(spread(math_I3,3,mesh_maxNips),4,mesh_NcpElems) ! initialize to identity - rewind(myFile) + rewind(myUnit) do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization) ! wind forward to - line = IO_read(myFile) + line = IO_read(myUnit) enddo do while (trim(line) /= '#EOF#') - line = IO_read(myFile) + line = IO_read(myUnit) if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'[',']') /= '') then ! next section @@ -199,8 +199,7 @@ subroutine homogenization_RGC_init(myFile) if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then do i = 1_pInt,maxNinstance - write(6,'(a15,1x,i4)') 'instance: ', i - write(6,*) + write(6,'(a15,1x,i4,/)') 'instance: ', i write(6,'(a25,3(1x,i8))') 'cluster size: ',(homogenization_RGC_Ngrains(j,i),j=1_pInt,3_pInt) write(6,'(a25,1x,e10.3)') 'scaling parameter: ', homogenization_RGC_xiAlpha(i) write(6,'(a25,1x,e10.3)') 'over-proportionality: ', homogenization_RGC_ciAlpha(i) @@ -228,11 +227,11 @@ subroutine homogenization_RGC_init(myFile) mySize = 0_pInt end select - if (mySize > 0_pInt) then ! any meaningful output found + outputFound: if (mySize > 0_pInt) then homogenization_RGC_sizePostResult(j,i) = mySize homogenization_RGC_sizePostResults(i) = & homogenization_RGC_sizePostResults(i) + mySize - endif + endif outputFound enddo homogenization_RGC_sizeState(i) & @@ -249,7 +248,7 @@ end subroutine homogenization_RGC_init !-------------------------------------------------------------------------------------------------- !> @brief partitions the deformation gradient onto the constituents !-------------------------------------------------------------------------------------------------- -subroutine homogenization_RGC_partitionDeformation(F,F0,avgF,state,ip,el) +subroutine homogenization_RGC_partitionDeformation(F,avgF,state,ip,el) use prec, only: & p_vec use debug, only: & @@ -268,7 +267,6 @@ subroutine homogenization_RGC_partitionDeformation(F,F0,avgF,state,ip,el) implicit none real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned F per grain - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: F0 !< initial partioned F per grain real(pReal), dimension (3,3), intent(in) :: avgF !< averaged F type(p_vec), intent(in) :: state integer(pInt), intent(in) :: & @@ -427,7 +425,7 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el !-------------------------------------------------------------------------------------------------- ! calculating volume discrepancy and stress penalty related to overall volume discrepancy - call homogenization_RGC_volumePenalty(D,volDiscrep,F,avgF,ip,el,homID) + call homogenization_RGC_volumePenalty(D,volDiscrep,F,avgF,ip,el) !-------------------------------------------------------------------------------------------------- ! debugging the mismatch, stress and penalties of grains @@ -513,6 +511,7 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el endif homogenization_RGC_updateState = .false. + !-------------------------------------------------------------------------------------------------- ! If convergence reached => done and happy if (residMax < relTol_RGC*stresMax .or. residMax < absTol_RGC) then @@ -521,8 +520,7 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & .and. debug_e == el .and. debug_i == ip) then !$OMP CRITICAL (write2out) - write(6,'(1x,a55)')'... done and happy' - write(6,*)' ' + write(6,'(1x,a55,/)')'... done and happy' flush(6) !$OMP END CRITICAL (write2out) endif @@ -552,16 +550,14 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & .and. debug_e == el .and. debug_i == ip) then !$OMP CRITICAL (write2out) - write(6,'(1x,a30,1x,e15.8)')'Constitutive work: ',constitutiveWork + write(6,'(1x,a30,1x,e15.8)') 'Constitutive work: ',constitutiveWork write(6,'(1x,a30,3(1x,e15.8))')'Magnitude mismatch: ',sum(NN(1,:))/real(nGrain,pReal), & sum(NN(2,:))/real(nGrain,pReal), & sum(NN(3,:))/real(nGrain,pReal) - write(6,'(1x,a30,1x,e15.8)')'Penalty energy: ',penaltyEnergy - write(6,'(1x,a30,1x,e15.8)')'Volume discrepancy: ',volDiscrep - write(6,*)'' - write(6,'(1x,a30,1x,e15.8)')'Maximum relaxation rate: ',maxval(abs(drelax))/dt - write(6,'(1x,a30,1x,e15.8)')'Average relaxation rate: ',sum(abs(drelax))/dt/real(3_pInt*nIntFaceTot,pReal) - write(6,*)'' + write(6,'(1x,a30,1x,e15.8)') 'Penalty energy: ',penaltyEnergy + write(6,'(1x,a30,1x,e15.8,/)') 'Volume discrepancy: ',volDiscrep + write(6,'(1x,a30,1x,e15.8)') 'Maximum relaxation rate: ',maxval(abs(drelax))/dt + write(6,'(1x,a30,1x,e15.8,/)') 'Average relaxation rate: ',sum(abs(drelax))/dt/real(3_pInt*nIntFaceTot,pReal) flush(6) !$OMP END CRITICAL (write2out) endif @@ -577,8 +573,7 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & .and. debug_e == el .and. debug_i == ip) then !$OMP CRITICAL (write2out) - write(6,'(1x,a55)')'... broken' - write(6,*)' ' + write(6,'(1x,a55,/)')'... broken' flush(6) !$OMP END CRITICAL (write2out) endif @@ -589,8 +584,7 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & .and. debug_e == el .and. debug_i == ip) then !$OMP CRITICAL (write2out) - write(6,'(1x,a55)')'... not yet done' - write(6,*)' ' + write(6,'(1x,a55,/)')'... not yet done' flush(6) !$OMP END CRITICAL (write2out) endif @@ -668,9 +662,9 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el p_relax = relax p_relax(ipert) = relax(ipert) + pPert_RGC ! perturb the relaxation vector state%p(1:3*nIntFaceTot) = p_relax - call homogenization_RGC_grainDeformation(pF,F0,avgF,state,ip,el) ! compute the grains deformation from perturbed state + call homogenization_RGC_grainDeformation(pF,avgF,state,ip,el) ! compute the grains deformation from perturbed state call homogenization_RGC_stressPenalty(pR,pNN,avgF,pF,ip,el,homID) ! compute stress penalty due to interface mismatch from perturbed state - call homogenization_RGC_volumePenalty(pD,volDiscrep,pF,avgF,ip,el,homID) ! compute stress penalty due to volume discrepancy from perturbed state + call homogenization_RGC_volumePenalty(pD,volDiscrep,pF,avgF,ip,el) ! compute stress penalty due to volume discrepancy from perturbed state !-------------------------------------------------------------------------------------------------- ! computing the global stress residual array from the perturbed state @@ -814,9 +808,7 @@ end function homogenization_RGC_updateState !-------------------------------------------------------------------------------------------------- !> @brief derive average stress and stiffness from constituent quantities !-------------------------------------------------------------------------------------------------- -subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ip,el ) - use prec, only: & - p_vec +subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,el) use debug, only: & debug_level, & debug_homogenization,& @@ -826,15 +818,14 @@ subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF, use math, only: math_Plain3333to99 implicit none - real(pReal), dimension (3,3), intent(out) :: avgP ! average stress at material point - real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF ! average stiffness at material point - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P ! array of current grain stresses - real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF ! array of current grain stiffnesses + real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point + real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point + real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P !< array of current grain stresses + real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF !< array of current grain stiffnesses + integer(pInt), intent(in) :: el !< element number real(pReal), dimension (9,9) :: dPdF99 - integer(pInt), intent(in) :: & - ip, & ! integration point number - el ! element number - integer(pInt) homID, i, j, Ngrains, iGrain + + integer(pInt) :: homID, i, j, Ngrains, iGrain homID = homogenization_typeInstance(mesh_element(3,el)) Ngrains = homogenization_Ngrains(mesh_element(3,el)) @@ -866,7 +857,7 @@ end subroutine homogenization_RGC_averageStressAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief return array of homogenization results for post file inclusion !-------------------------------------------------------------------------------------------------- -pure function homogenization_RGC_postResults(state,ip,el) +pure function homogenization_RGC_postResults(state,el) use prec, only: & p_vec use mesh, only: & @@ -877,9 +868,7 @@ pure function homogenization_RGC_postResults(state,ip,el) implicit none type(p_vec), intent(in) :: state ! my State - integer(pInt), intent(in) :: & - ip, & ! integration point number - el ! element number + integer(pInt), intent(in) :: el ! element number integer(pInt) homID,o,c,nIntFaceTot real(pReal), dimension(homogenization_RGC_sizePostResults(homogenization_typeInstance(mesh_element(3,el)))) :: & homogenization_RGC_postResults @@ -943,10 +932,10 @@ subroutine homogenization_RGC_stressPenalty(rPen,nMis,avgF,fDef,ip,el,homID) xSmoo_RGC implicit none - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: rPen ! stress-like penalty - real(pReal), dimension (3,homogenization_maxNgrains), intent(out) :: nMis ! total amount of mismatch - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: fDef ! deformation gradients - real(pReal), dimension (3,3), intent(in) :: avgF ! initial effective stretch tensor + real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: rPen !< stress-like penalty + real(pReal), dimension (3,homogenization_maxNgrains), intent(out) :: nMis !< total amount of mismatch + real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: fDef !< deformation gradients + real(pReal), dimension (3,3), intent(in) :: avgF !< initial effective stretch tensor integer(pInt), intent(in) :: ip,el integer(pInt), dimension (4) :: intFace integer(pInt), dimension (3) :: iGrain3,iGNghb3,nGDim @@ -960,7 +949,6 @@ subroutine homogenization_RGC_stressPenalty(rPen,nMis,avgF,fDef,ip,el,homID) real(pReal), parameter :: nDefToler = 1.0e-10_pReal nGDim = homogenization_RGC_Ngrains(1:3,homID) - rPen = 0.0_pReal nMis = 0.0_pReal @@ -1064,7 +1052,7 @@ end subroutine homogenization_RGC_stressPenalty !-------------------------------------------------------------------------------------------------- !> @brief calculate stress-like penalty due to volume discrepancy !-------------------------------------------------------------------------------------------------- -subroutine homogenization_RGC_volumePenalty(vPen,vDiscrep,fDef,fAvg,ip,el, homID) +subroutine homogenization_RGC_volumePenalty(vPen,vDiscrep,fDef,fAvg,ip,el) use debug, only: & debug_level, & debug_homogenization,& @@ -1092,7 +1080,7 @@ subroutine homogenization_RGC_volumePenalty(vPen,vDiscrep,fDef,fAvg,ip,el, homID integer(pInt), intent(in) :: ip,& ! integration point el real(pReal), dimension (homogenization_maxNgrains) :: gVol - integer(pInt) :: homID,iGrain,nGrain,i,j + integer(pInt) :: iGrain,nGrain,i,j nGrain = homogenization_Ngrains(mesh_element(3,el)) @@ -1434,7 +1422,7 @@ end function homogenization_RGC_interface1to4 !> @brief calculating the grain deformation gradient (the same with ! homogenization_RGC_partionDeformation, but used only for perturbation scheme) !-------------------------------------------------------------------------------------------------- -subroutine homogenization_RGC_grainDeformation(F, F0, avgF, state, ip, el) +subroutine homogenization_RGC_grainDeformation(F, avgF, state, ip, el) use prec, only: & p_vec use mesh, only: & @@ -1446,7 +1434,6 @@ subroutine homogenization_RGC_grainDeformation(F, F0, avgF, state, ip, el) implicit none real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned F per grain - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: F0 !< initiatial partioned F per grain real(pReal), dimension (3,3), intent(in) :: avgF !< type(p_vec), intent(in) :: state integer(pInt), intent(in) :: & diff --git a/code/homogenization_isostrain.f90 b/code/homogenization_isostrain.f90 index f863a9f20..a548e9de6 100644 --- a/code/homogenization_isostrain.f90 +++ b/code/homogenization_isostrain.f90 @@ -55,7 +55,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -subroutine homogenization_isostrain_init(myFile) +subroutine homogenization_isostrain_init(myUnit) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use math, only: & math_Mandel3333to66, & @@ -64,7 +64,7 @@ subroutine homogenization_isostrain_init(myFile) use material implicit none - integer(pInt), intent(in) :: myFile + integer(pInt), intent(in) :: myUnit integer(pInt), parameter :: MAXNCHUNKS = 2_pInt integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions integer(pInt) :: & @@ -96,15 +96,15 @@ subroutine homogenization_isostrain_init(myFile) allocate(homogenization_isostrain_output(maxval(homogenization_Noutput),maxNinstance)) homogenization_isostrain_output = '' - rewind(myFile) + rewind(myUnit) section = 0_pInt do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization) ! wind forward to - line = IO_read(myFile) + line = IO_read(myUnit) enddo do while (trim(line) /= '#EOF#') - line = IO_read(myFile) + line = IO_read(myUnit) if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'[',']') /= '') then ! next section @@ -120,7 +120,7 @@ subroutine homogenization_isostrain_init(myFile) case ('(output)') output = output + 1_pInt homogenization_isostrain_output(output,i) = IO_lc(IO_stringValue(line,positions,2_pInt)) - case ('ngrains') + case ('ngrains','ncomponents') homogenization_isostrain_Ngrains(i) = IO_intValue(line,positions,2_pInt) case ('mapping') homogenization_isostrain_mapping(i) = IO_lc(IO_stringValue(line,positions,2_pInt)) @@ -134,17 +134,17 @@ subroutine homogenization_isostrain_init(myFile) do j = 1_pInt,maxval(homogenization_Noutput) select case(homogenization_isostrain_output(j,i)) - case('ngrains') + case('ngrains','ncomponents') mySize = 1_pInt case default mySize = 0_pInt end select - if (mySize > 0_pInt) then ! any meaningful output found + outputFound: if (mySize > 0_pInt) then homogenization_isostrain_sizePostResult(j,i) = mySize homogenization_isostrain_sizePostResults(i) = & - homogenization_isostrain_sizePostResults(i) + mySize - endif + homogenization_isostrain_sizePostResults(i) + mySize + endif outputFound enddo enddo @@ -154,10 +154,9 @@ end subroutine homogenization_isostrain_init !-------------------------------------------------------------------------------------------------- !> @brief partitions the deformation gradient onto the constituents !-------------------------------------------------------------------------------------------------- -subroutine homogenization_isostrain_partitionDeformation(F,F0,avgF,state,ip,el) +subroutine homogenization_isostrain_partitionDeformation(F,avgF,el) use prec, only: & - pReal, & - p_vec + pReal use mesh, only: & mesh_element use material, only: & @@ -165,12 +164,9 @@ subroutine homogenization_isostrain_partitionDeformation(F,F0,avgF,state,ip,el) homogenization_Ngrains implicit none - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F ! partioned def grad per grain - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: F0 ! initial partioned def grad per grain - real(pReal), dimension (3,3), intent(in) :: avgF ! my average def grad - type(p_vec), intent(in) :: state ! my state + real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned def grad per grain + real(pReal), dimension (3,3), intent(in) :: avgF !< my average def grad integer(pInt), intent(in) :: & - ip, & !< integration point number el !< element number F = spread(avgF,3,homogenization_Ngrains(mesh_element(3,el))) @@ -181,7 +177,7 @@ end subroutine homogenization_isostrain_partitionDeformation !-------------------------------------------------------------------------------------------------- !> @brief derive average stress and stiffness from constituent quantities !-------------------------------------------------------------------------------------------------- -subroutine homogenization_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ip,el) +subroutine homogenization_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,el) use prec, only: & pReal use mesh, only: & @@ -196,9 +192,7 @@ subroutine homogenization_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P !< array of current grain stresses real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF !< array of current grain stiffnesses - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number + integer(pInt), intent(in) :: el !< element number integer(pInt) :: & homID, & Ngrains @@ -224,10 +218,9 @@ end subroutine homogenization_isostrain_averageStressAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief return array of homogenization results for post file inclusion !-------------------------------------------------------------------------------------------------- -pure function homogenization_isostrain_postResults(state,ip,el) +pure function homogenization_isostrain_postResults(el) use prec, only: & - pReal,& - p_vec + pReal use mesh, only: & mesh_element use material, only: & @@ -235,10 +228,7 @@ pure function homogenization_isostrain_postResults(state,ip,el) homogenization_Noutput implicit none - type(p_vec), intent(in) :: state - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number + integer(pInt), intent(in) :: el !< element number real(pReal), dimension(homogenization_isostrain_sizePostResults & (homogenization_typeInstance(mesh_element(3,el)))) :: & homogenization_isostrain_postResults @@ -253,7 +243,7 @@ pure function homogenization_isostrain_postResults(state,ip,el) do o = 1_pInt,homogenization_Noutput(mesh_element(3,el)) select case(homogenization_isostrain_output(o,homID)) - case ('ngrains') + case ('ngrains','ncomponents') homogenization_isostrain_postResults(c+1_pInt) = real(homogenization_isostrain_Ngrains(homID),pReal) c = c + 1_pInt end select diff --git a/code/material.f90 b/code/material.f90 index 4bfb0e411..05106352c 100644 --- a/code/material.f90 +++ b/code/material.f90 @@ -122,14 +122,17 @@ module material homogenization_active - public :: material_init + public :: & + material_init - private :: material_parseHomogenization, & - material_parseMicrostructure, & - material_parseCrystallite, & - material_parsePhase, & - material_parseTexture, & - material_populateGrains + private :: & + material_parseHomogenization, & + material_parseMicrostructure, & + material_parseCrystallite, & + material_parsePhase, & + material_parseTexture, & + material_populateGrains + contains diff --git a/code/mesh.f90 b/code/mesh.f90 index 955ae6fa6..73a253644 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -669,7 +669,7 @@ integer(pInt) function mesh_FEasCP(what,myID) return endif - do while (upper-lower > 1_pInt) ! binary search in between bounds + binarySearch: do while (upper-lower > 1_pInt) center = (lower+upper)/2_pInt if (lookupMap(1_pInt,center) < myID) then lower = center @@ -679,7 +679,7 @@ integer(pInt) function mesh_FEasCP(what,myID) mesh_FEasCP = lookupMap(2_pInt,center) exit endif - enddo + enddo binarySearch end function mesh_FEasCP diff --git a/code/numerics.f90 b/code/numerics.f90 index 588e647ab..0428f8d52 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -60,7 +60,6 @@ module numerics 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 - rTol_crystalliteTemperature= 1.0e-6_pReal, & !< relative tolerance in crystallite temperature loop rTol_crystalliteStress = 1.0e-6_pReal, & !< relative tolerance in crystallite stress loop aTol_crystalliteStress = 1.0e-8_pReal, & !< absolute tolerance in crystallite stress loop, Default 1.0e-8: residuum is in Lp and hence strain is on this order numerics_unitlength = 1.0_pReal, & !< determines the physical length of one computational length unit @@ -222,8 +221,6 @@ subroutine numerics_init stepIncreaseHomog = IO_floatValue(line,positions,2_pInt) case ('rtol_crystallitestate') rTol_crystalliteState = IO_floatValue(line,positions,2_pInt) - case ('rtol_crystallitetemperature') - rTol_crystalliteTemperature = IO_floatValue(line,positions,2_pInt) case ('rtol_crystallitestress') rTol_crystalliteStress = IO_floatValue(line,positions,2_pInt) case ('atol_crystallitestress') @@ -377,7 +374,6 @@ subroutine numerics_init write(6,'(a24,1x,i8)') ' nState: ',nState write(6,'(a24,1x,i8)') ' nStress: ',nStress write(6,'(a24,1x,es8.1)') ' rTol_crystalliteState: ',rTol_crystalliteState - write(6,'(a24,1x,es8.1)') ' rTol_crystalliteTemp: ',rTol_crystalliteTemperature write(6,'(a24,1x,es8.1)') ' rTol_crystalliteStress: ',rTol_crystalliteStress write(6,'(a24,1x,es8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress write(6,'(a24,2(1x,i8))') ' integrator: ',numerics_integrator @@ -469,7 +465,6 @@ subroutine numerics_init 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') if (rTol_crystalliteState <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteState') - if (rTol_crystalliteTemperature <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteTemperature') if (rTol_crystalliteStress <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteStress') if (aTol_crystalliteStress <= 0.0_pReal) call IO_error(301_pInt,ext_msg='aTol_crystalliteStress') if (any(numerics_integrator <= 0_pInt) .or. any(numerics_integrator >= 6_pInt)) &