diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index c8b7e3cca..37e4ba737 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -192,11 +192,10 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS else validCalculation if (debugCPFEM%extensive) print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip - call homogenization_mechanical_response(dt,[ip,ip],[elCP,elCP]) + call homogenization_mechanical_response(dt,(elCP-1)*discretization_nIPs + ip,(elCP-1)*discretization_nIPs + ip) if (.not. terminallyIll) & call homogenization_mechanical_response2(dt,[ip,ip],[elCP,elCP]) - terminalIllness: if (terminallyIll) then call random_number(rnd) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index d95580145..34a976644 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -812,9 +812,9 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& homogenization_F = reshape(F,[3,3,product(cells(1:2))*cells3]) ! set materialpoint target F to estimated field - call homogenization_mechanical_response(Delta_t,[1,1],[1,product(cells(1:2))*cells3]) ! calculate P field + call homogenization_mechanical_response(Delta_t,1,product(cells(1:2))*cells3) ! calculate P field if (.not. terminallyIll) & - call homogenization_thermal_response(Delta_t,[1,1],[1,product(cells(1:2))*cells3]) + call homogenization_thermal_response(Delta_t,1,product(cells(1:2))*cells3) if (.not. terminallyIll) & call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(cells(1:2))*cells3]) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 3d96b007f..8539df994 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -222,14 +222,13 @@ end subroutine homogenization_init !-------------------------------------------------------------------------------------------------- !> @brief !-------------------------------------------------------------------------------------------------- -subroutine homogenization_mechanical_response(Delta_t,FEsolving_execIP,FEsolving_execElem) +subroutine homogenization_mechanical_response(Delta_t,cell_start,cell_end) real(pReal), intent(in) :: Delta_t !< time increment - integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP + integer, intent(in) :: & + cell_start, cell_end integer :: & NiterationMPstate, & - ip, & !< integration point number - el, & !< element number co, ce, ho, en logical :: & converged @@ -237,45 +236,42 @@ subroutine homogenization_mechanical_response(Delta_t,FEsolving_execIP,FEsolving doneAndHappy - !$OMP PARALLEL DO PRIVATE(ce,en,ho,NiterationMPstate,converged,doneAndHappy) - do el = FEsolving_execElem(1),FEsolving_execElem(2) + !$OMP PARALLEL DO PRIVATE(en,ho,co,NiterationMPstate,converged,doneAndHappy) + do ce = cell_start, cell_end - do ip = FEsolving_execIP(1),FEsolving_execIP(2) - ce = (el-1)*discretization_nIPs + ip - en = material_homogenizationEntry(ce) - ho = material_homogenizationID(ce) + en = material_homogenizationEntry(ce) + ho = material_homogenizationID(ce) - call phase_restore(ce,.false.) ! wrong name (is more a forward function) + call phase_restore(ce,.false.) ! wrong name (is more a forward function) - if(homogState(ho)%sizeState > 0) homogState(ho)%state(:,en) = homogState(ho)%state0(:,en) - if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%state(:,en) = damageState_h(ho)%state0(:,en) - call damage_partition(ce) + if(homogState(ho)%sizeState > 0) homogState(ho)%state(:,en) = homogState(ho)%state0(:,en) + if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%state(:,en) = damageState_h(ho)%state0(:,en) + call damage_partition(ce) - doneAndHappy = [.false.,.true.] + doneAndHappy = [.false.,.true.] - NiterationMPstate = 0 - convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)) & - .and. NiterationMPstate < num%nMPstate) - NiterationMPstate = NiterationMPstate + 1 + NiterationMPstate = 0 + convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)) & + .and. NiterationMPstate < num%nMPstate) + NiterationMPstate = NiterationMPstate + 1 - call mechanical_partition(homogenization_F(1:3,1:3,ce),ce) - converged = all([(phase_mechanical_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))]) - if (converged) then - doneAndHappy = mechanical_updateState(Delta_t,homogenization_F(1:3,1:3,ce),ce) - converged = all(doneAndHappy) - else - doneAndHappy = [.true.,.false.] - endif - enddo convergenceLooping + call mechanical_partition(homogenization_F(1:3,1:3,ce),ce) + converged = all([(phase_mechanical_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))]) + if (converged) then + doneAndHappy = mechanical_updateState(Delta_t,homogenization_F(1:3,1:3,ce),ce) + converged = all(doneAndHappy) + else + doneAndHappy = [.true.,.false.] + end if + end do convergenceLooping - converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))]) + converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))]) - if (.not. converged) then - if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill' - terminallyIll = .true. - endif - enddo - enddo + if (.not. converged) then + if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill' + terminallyIll = .true. + end if + end do !$OMP END PARALLEL DO end subroutine homogenization_mechanical_response @@ -284,31 +280,27 @@ end subroutine homogenization_mechanical_response !-------------------------------------------------------------------------------------------------- !> @brief !-------------------------------------------------------------------------------------------------- -subroutine homogenization_thermal_response(Delta_t,FEsolving_execIP,FEsolving_execElem) +subroutine homogenization_thermal_response(Delta_t,cell_start,cell_end) real(pReal), intent(in) :: Delta_t !< time increment - integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP + integer, intent(in) :: & + cell_start, cell_end integer :: & - ip, & !< integration point number - el, & !< element number co, ce, ho - !$OMP PARALLEL DO PRIVATE(ho,ce) - do el = FEsolving_execElem(1),FEsolving_execElem(2) + !$OMP PARALLEL DO PRIVATE(ho) + do ce = cell_start, cell_end if (terminallyIll) continue - do ip = FEsolving_execIP(1),FEsolving_execIP(2) - ce = (el-1)*discretization_nIPs + ip - ho = material_homogenizationID(ce) - call thermal_partition(ce) - do co = 1, homogenization_Nconstituents(ho) - if (.not. phase_thermal_constitutive(Delta_t,material_phaseID(co,ce),material_phaseEntry(co,ce))) then - if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill' - terminallyIll = .true. - endif - enddo - enddo - enddo + ho = material_homogenizationID(ce) + call thermal_partition(ce) + do co = 1, homogenization_Nconstituents(ho) + if (.not. phase_thermal_constitutive(Delta_t,material_phaseID(co,ce),material_phaseEntry(co,ce))) then + if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill' + terminallyIll = .true. + end if + end do + end do !$OMP END PARALLEL DO end subroutine homogenization_thermal_response @@ -331,11 +323,11 @@ subroutine homogenization_mechanical_response2(Delta_t,FEsolving_execIP,FEsolvin elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2) IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2) ce = (el-1)*discretization_nIPs + ip - ho = material_homogenizationID(ce) - do co = 1, homogenization_Nconstituents(ho) - call crystallite_orientations(co,ip,el) + ho = material_homogenizationID(ce) + do co = 1, homogenization_Nconstituents(ho) + call crystallite_orientations(co,ip,el) enddo - call mechanical_homogenize(Delta_t,ce) + call mechanical_homogenize(Delta_t,ce) enddo IpLooping3 enddo elementLooping3 !$OMP END PARALLEL DO diff --git a/src/material.f90 b/src/material.f90 index fd1531964..bf78a77a6 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -17,16 +17,17 @@ module material implicit none private - type :: tRotationContainer - type(tRotation), dimension(:), allocatable :: data - end type - type :: tTensorContainer + type, public :: tRotationContainer + type(tRotation), dimension(:), allocatable :: data + end type tRotationContainer + + type, public :: tTensorContainer real(pReal), dimension(:,:,:), allocatable :: data - end type + end type tTensorContainer - type(tRotationContainer), dimension(:), allocatable :: material_O_0 - type(tTensorContainer), dimension(:), allocatable :: material_F_i_0 + type(tRotationContainer), dimension(:), allocatable, public, protected :: material_O_0 + type(tTensorContainer), dimension(:), allocatable, public, protected :: material_F_i_0 integer, dimension(:), allocatable, public, protected :: & homogenization_Nconstituents !< number of grains in each homogenization @@ -37,20 +38,14 @@ module material material_name_phase, & !< name of each phase material_name_homogenization !< name of each homogenization - integer, dimension(:), allocatable, public, protected :: & ! (elem) - material_homogenizationID, & !< per cell TODO: material_ID_homogenization - material_homogenizationEntry !< per cell TODO: material_entry_homogenization - integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem) - material_phaseAt, & !< phase ID of each element TODO: remove - material_phaseID, & !< per (constituent,cell) TODO: material_ID_phase - material_phaseEntry !< per (constituent,cell) TODO: material_entry_phase - integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem) - material_phaseMemberAt !TODO: remove + integer, dimension(:), allocatable, public, protected :: & ! (cell) + material_homogenizationID, & ! TODO: rename to material_ID_homogenization + material_homogenizationEntry ! TODO: rename to material_entry_homogenization + integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,cell) + material_phaseID, & ! TODO: rename to material_ID_phase + material_phaseEntry ! TODO: rename to material_entry_phase + public :: & - tTensorContainer, & - tRotationContainer, & - material_F_i_0, & - material_O_0, & material_init contains @@ -97,11 +92,12 @@ subroutine parse() counterPhase, & counterHomogenization - real(pReal) :: & - frac + real(pReal) :: v integer :: & - el, ip, co, ma, & - h, ce + el, ip, & + ho, ph, & + co, ce, & + ma materials => config_material%get('material') phases => config_material%get('phase') @@ -118,51 +114,48 @@ subroutine parse() #endif allocate(homogenization_Nconstituents(homogenizations%length)) - do h=1, homogenizations%length - homogenization => homogenizations%get(h) - homogenization_Nconstituents(h) = homogenization%get_asInt('N_constituents') + do ho=1, homogenizations%length + homogenization => homogenizations%get(ho) + homogenization_Nconstituents(ho) = homogenization%get_asInt('N_constituents') end do homogenization_maxNconstituents = maxval(homogenization_Nconstituents) allocate(counterPhase(phases%length),source=0) allocate(counterHomogenization(homogenizations%length),source=0) - allocate(material_phaseAt(homogenization_maxNconstituents,discretization_Nelems),source=0) - allocate(material_phaseMemberAt(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems),source=0) - - allocate(material_homogenizationID(discretization_nIPs*discretization_Nelems),source=0) allocate(material_homogenizationEntry(discretization_nIPs*discretization_Nelems),source=0) + allocate(material_phaseID(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) allocate(material_phaseEntry(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) do el = 1, discretization_Nelems - material => materials%get(discretization_materialAt(el)) - constituents => material%get('constituents') + material => materials%get(discretization_materialAt(el)) + ho = homogenizations%getIndex(material%get_asString('homogenization')) do ip = 1, discretization_nIPs ce = (el-1)*discretization_nIPs + ip - material_homogenizationID(ce) = homogenizations%getIndex(material%get_asString('homogenization')) - counterHomogenization(material_homogenizationID(ce)) = counterHomogenization(material_homogenizationID(ce)) + 1 - material_homogenizationEntry(ce) = counterHomogenization(material_homogenizationID(ce)) + material_homogenizationID(ce) = ho + counterHomogenization(ho) = counterHomogenization(ho) + 1 + material_homogenizationEntry(ce) = counterHomogenization(ho) end do - frac = 0.0_pReal + v = 0.0_pReal + constituents => material%get('constituents') do co = 1, constituents%length constituent => constituents%get(co) - frac = frac + constituent%get_asFloat('v') + v = v + constituent%get_asFloat('v') - material_phaseAt(co,el) = phases%getIndex(constituent%get_asString('phase')) + ph = phases%getIndex(constituent%get_asString('phase')) do ip = 1, discretization_nIPs ce = (el-1)*discretization_nIPs + ip - counterPhase(material_phaseAt(co,el)) = counterPhase(material_phaseAt(co,el)) + 1 - material_phaseMemberAt(co,ip,el) = counterPhase(material_phaseAt(co,el)) - material_phaseEntry(co,ce) = counterPhase(material_phaseAt(co,el)) - material_phaseID(co,ce) = material_phaseAt(co,el) + material_phaseID(co,ce) = ph + counterPhase(ph) = counterPhase(ph) + 1 + material_phaseEntry(co,ce) = counterPhase(ph) end do end do - if (dNeq(frac,1.0_pReal,1.e-12_pReal)) call IO_error(153,ext_msg='constituent') + if (dNeq(v,1.0_pReal,1.e-12_pReal)) call IO_error(153,ext_msg='constituent') end do diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 703c6c818..4e31a0475 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -150,7 +150,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) print'(/,1x,a)', '... evaluating constitutive response ......................................' - call homogenization_mechanical_response(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field + call homogenization_mechanical_response(timeinc,1,mesh_maxNips*mesh_NcpElems) ! calculate P field if (.not. terminallyIll) & call homogenization_mechanical_response2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) cutBack = .false. diff --git a/src/phase.f90 b/src/phase.f90 index 3de551131..c161fb48c 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -228,15 +228,15 @@ module phase end function phase_thermal_constitutive - module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_) + module function phase_damage_constitutive(Delta_t,co,ce) result(converged_) real(pReal), intent(in) :: Delta_t - integer, intent(in) :: co, ip, el + integer, intent(in) :: co, ce logical :: converged_ end function phase_damage_constitutive - module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged_) + module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) real(pReal), intent(in) :: Delta_t - integer, intent(in) :: co, ip, el + integer, intent(in) :: co, ce logical :: converged_ end function phase_mechanical_constitutive @@ -264,19 +264,18 @@ module phase real(pReal) :: f end function phase_f_T - module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) + module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el) integer, intent(in) :: & ph, & - i, & - e + ip, & + el type(tRotationContainer), dimension(:), intent(in) :: orientation end subroutine plastic_nonlocal_updateCompatibility - module subroutine plastic_dependentState(co,ip,el) + module subroutine plastic_dependentState(ph,en) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + ph, & + en end subroutine plastic_dependentState @@ -503,8 +502,8 @@ subroutine crystallite_init() ce, & co, & !< counter in integration point component loop ip, & !< counter in integration point loop - el !< counter in element loop - + el, & !< counter in element loop + en, ph class(tNode), pointer :: & num_crystallite, & phases @@ -543,13 +542,15 @@ subroutine crystallite_init() phases => config_material%get('phase') - !$OMP PARALLEL DO PRIVATE(ce) + !$OMP PARALLEL DO PRIVATE(ce,ph,en) do el = 1, discretization_Nelems do ip = 1, discretization_nIPs ce = (el-1)*discretization_nIPs + ip do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) + en = material_phaseEntry(co,ce) + ph = material_phaseID(co,ce) call crystallite_orientations(co,ip,el) - call plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states + call plastic_dependentState(ph,en) ! update dependent state variables to be consistent with basic states end do end do end do @@ -577,8 +578,8 @@ subroutine crystallite_orientations(co,ip,el) call phase_O(ph)%data(en)%fromMatrix(transpose(math_rotationalPart(mechanical_F_e(ph,en)))) - if (plasticState(material_phaseAt(1,el))%nonlocal) & - call plastic_nonlocal_updateCompatibility(phase_O,material_phaseAt(1,el),ip,el) + if (plasticState(material_phaseID(1,(el-1)*discretization_nIPs + ip))%nonlocal) & + call plastic_nonlocal_updateCompatibility(phase_O,material_phaseID(1,(el-1)*discretization_nIPs + ip),ip,el) end subroutine crystallite_orientations diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 2fa5ab9ea..e4700938f 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -127,21 +127,20 @@ end subroutine damage_init !-------------------------------------------------------------------------------------------------- !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- -module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_) +module function phase_damage_constitutive(Delta_t,co,ce) result(converged_) real(pReal), intent(in) :: Delta_t integer, intent(in) :: & co, & - ip, & - el + ce logical :: converged_ integer :: & ph, en - ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) - en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) + ph = material_phaseID(co,ce) + en = material_phaseEntry(co,ce) converged_ = .not. integrateDamageState(Delta_t,ph,en) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index eb758d8ed..05866a3e1 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -79,15 +79,12 @@ submodule(phase) mechanical en end subroutine plastic_isotropic_LiAndItsTangent - module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState) + module function plastic_dotState(subdt,ph,en) result(dotState) integer, intent(in) :: & - co, & !< constituent - ip, & !< integration point - el, & !< element ph, & en real(pReal), intent(in) :: & - subdt !< timestep + subdt !< timestep real(pReal), dimension(plasticState(ph)%sizeDotState) :: & dotState end function plastic_dotState @@ -365,13 +362,11 @@ end subroutine mechanical_results !> @brief calculation of stress (P) with time integration based on a residuum in Lp and !> intermediate acceleration of the Newton-Raphson correction !-------------------------------------------------------------------------------------------------- -function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) +function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken) real(pReal), dimension(3,3), intent(in) :: F,subFp0,subFi0 real(pReal), intent(in) :: Delta_t - integer, intent(in):: el, & ! element index - ip, & ! integration point index - co ! grain index + integer, intent(in) :: ph, en real(pReal), dimension(3,3):: Fp_new, & ! plastic deformation gradient at end of timestep invFp_new, & ! inverse of Fp_new @@ -419,19 +414,13 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) ierr, & ! error indicator for LAPACK o, & p, & - ph, & - en, & jacoCounterLp, & jacoCounterLi ! counters to check for Jacobian update logical :: error,broken broken = .true. - - ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) - en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) - - call plastic_dependentState(co,ip,el) + call plastic_dependentState(ph,en) Lpguess = phase_mechanical_Lp(ph)%data(1:3,1:3,en) ! take as first guess Liguess = phase_mechanical_Li(ph)%data(1:3,1:3,en) ! take as first guess @@ -579,37 +568,32 @@ end function integrateStress !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken) +function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t integer, intent(in) :: & - el, & !< element index in element loop - ip, & !< integration point index in ip loop - co !< grain index in grain loop + ph, & + en logical :: & broken integer :: & NiterationState, & !< number of iterations in state loop - ph, & - en, & sizeDotState real(pReal) :: & zeta - real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & r, & ! state residuum dotState - real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState,2) :: & + real(pReal), dimension(plasticState(ph)%sizeDotState,2) :: & dotState_last - ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) - en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) broken = .true. - dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) return sizeDotState = plasticState(ph)%sizeDotState @@ -620,10 +604,10 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul dotState_last(1:sizeDotState,2) = merge(dotState_last(1:sizeDotState,1),0.0, nIterationState > 1) dotState_last(1:sizeDotState,1) = dotState - broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) + broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en) if(broken) exit iteration - dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) exit iteration zeta = damper(dotState,dotState_last(1:sizeDotState,1),dotState_last(1:sizeDotState,2)) @@ -672,31 +656,26 @@ end function integrateStateFPI !-------------------------------------------------------------------------------------------------- !> @brief integrate state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- -function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken) +function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t integer, intent(in) :: & - el, & !< element index in element loop - ip, & !< integration point index in ip loop - co !< grain index in grain loop + ph, & + en !< grain index in grain loop logical :: & broken - real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & dotState integer :: & - ph, & - en, & sizeDotState - ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) - en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) broken = .true. - dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) return sizeDotState = plasticState(ph)%sizeDotState @@ -706,7 +685,7 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res broken = plastic_deltaState(ph,en) if(broken) return - broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) + broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en) end function integrateStateEuler @@ -714,32 +693,27 @@ end function integrateStateEuler !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- -function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken) +function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t integer, intent(in) :: & - el, & !< element index in element loop - ip, & !< integration point index in ip loop - co !< grain index in grain loop + ph, & + en logical :: & broken integer :: & - ph, & - en, & sizeDotState - real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & r, & dotState - ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) - en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) broken = .true. - dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) return sizeDotState = plasticState(ph)%sizeDotState @@ -751,10 +725,10 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip broken = plastic_deltaState(ph,en) if(broken) return - broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) + broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en) if(broken) return - dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) return broken = .not. converged(r + 0.5_pReal * dotState * Delta_t, & @@ -767,12 +741,12 @@ end function integrateStateAdaptiveEuler !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the classic Runge Kutta method !--------------------------------------------------------------------------------------------------- -function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken) +function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t - integer, intent(in) :: co,ip,el + integer, intent(in) :: ph, en logical :: broken real(pReal), dimension(3,3), parameter :: & @@ -787,7 +761,7 @@ function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal] - broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C) + broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C) end function integrateStateRK4 @@ -795,12 +769,12 @@ end function integrateStateRK4 !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the Cash-Carp method !--------------------------------------------------------------------------------------------------- -function integrateStateRKCK45(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken) +function integrateStateRKCK45(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t - integer, intent(in) :: co,ip,el + integer, intent(in) :: ph, en logical :: broken real(pReal), dimension(5,5), parameter :: & @@ -822,7 +796,7 @@ function integrateStateRKCK45(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) re 13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 1._pReal/4._pReal] - broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,DB) + broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB) end function integrateStateRKCK45 @@ -831,7 +805,7 @@ end function integrateStateRKCK45 !> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an !! embedded explicit Runge-Kutta method !-------------------------------------------------------------------------------------------------- -function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,DB) result(broken) +function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 real(pReal), intent(in),dimension(:) :: subState0 @@ -840,28 +814,23 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D real(pReal), dimension(:), intent(in) :: B, C real(pReal), dimension(:), intent(in), optional :: DB integer, intent(in) :: & - el, & !< element index in element loop - ip, & !< integration point index in ip loop - co !< grain index in grain loop + ph, & + en logical :: broken integer :: & stage, & ! stage index in integration stage loop n, & - ph, & - en, & sizeDotState - real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & dotState - real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState,size(B)) :: & + real(pReal), dimension(plasticState(ph)%sizeDotState,size(B)) :: & plastic_RKdotState - ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) - en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) broken = .true. - dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) return sizeDotState = plasticState(ph)%sizeDotState @@ -879,10 +848,10 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D plasticState(ph)%state(1:sizeDotState,en) = subState0 & + dotState * Delta_t - broken = integrateStress(F_0 + (F-F_0) * Delta_t*C(stage),subFp0,subFi0,Delta_t*C(stage),co,ip,el) + broken = integrateStress(F_0+(F-F_0)*Delta_t*C(stage),subFp0,subFi0,Delta_t*C(stage), ph,en) if(broken) exit - dotState = plastic_dotState(Delta_t*C(stage), co,ip,el,ph,en) + dotState = plastic_dotState(Delta_t*C(stage), ph,en) if (any(IEEE_is_NaN(dotState))) exit enddo @@ -904,7 +873,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D broken = plastic_deltaState(ph,en) if(broken) return - broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) + broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en) end function integrateStateRK @@ -1006,13 +975,12 @@ end subroutine mechanical_forward !-------------------------------------------------------------------------------------------------- !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- -module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged_) +module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) real(pReal), intent(in) :: Delta_t integer, intent(in) :: & co, & - ip, & - el + ce logical :: converged_ real(pReal) :: & @@ -1028,16 +996,15 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged subLi0, & subF0, & subF - real(pReal), dimension(:), allocatable :: subState0 + real(pReal), dimension(plasticState(material_phaseID(co,ce))%sizeState) :: subState0 - ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) - en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) - sizeDotState = plasticState(ph)%sizeDotState + ph = material_phaseID(co,ce) + en = material_phaseEntry(co,ce) + subState0 = plasticState(ph)%state0(:,en) subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en) subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en) - allocate(subState0,source=plasticState(ph)%State0(:,en)) subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en) subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,en) @@ -1070,7 +1037,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged subStep = num%subStepSizeCryst * subStep phase_mechanical_Fp(ph)%data(1:3,1:3,en) = subFp0 phase_mechanical_Fi(ph)%data(1:3,1:3,en) = subFi0 - phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en) ! why no subS0 ? is S0 of any use? + phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en) if (subStep < 1.0_pReal) then ! actual (not initial) cutback phase_mechanical_Lp(ph)%data(1:3,1:3,en) = subLp0 phase_mechanical_Li(ph)%data(1:3,1:3,en) = subLi0 @@ -1082,9 +1049,10 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged !-------------------------------------------------------------------------------------------------- ! prepare for integration if (todo) then + sizeDotState = plasticState(ph)%sizeDotState subF = subF0 & + subStep * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en)) - converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,co,ip,el) + converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,ph,en) endif enddo cutbackLooping diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 70181d724..72b67ef64 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -160,16 +160,14 @@ submodule(phase:mechanical) plastic dotState end function dislotungsten_dotState - module subroutine nonlocal_dotState(Mp,timestep,ph,en,ip,el) + module subroutine nonlocal_dotState(Mp,timestep,ph,en) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress real(pReal), intent(in) :: & timestep !< substepped crystallite time increment integer, intent(in) :: & ph, & - en, & - ip, & !< current integration point - el !< current element number + en end subroutine nonlocal_dotState module subroutine dislotwin_dependentState(T,ph,en) @@ -186,12 +184,10 @@ submodule(phase:mechanical) plastic en end subroutine dislotungsten_dependentState - module subroutine nonlocal_dependentState(ph, en, ip, el) + module subroutine nonlocal_dependentState(ph,en) integer, intent(in) :: & ph, & - en, & - ip, & !< current integration point - el !< current element number + en end subroutine nonlocal_dependentState module subroutine plastic_kinehardening_deltaState(Mp,ph,en) @@ -301,12 +297,9 @@ end subroutine plastic_LpAndItsTangents !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState) +module function plastic_dotState(subdt,ph,en) result(dotState) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el, & !< element ph, & en real(pReal), intent(in) :: & @@ -340,7 +333,7 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState) dotState = dislotungsten_dotState(Mp,ph,en) case (PLASTIC_NONLOCAL_ID) plasticType - call nonlocal_dotState(Mp,subdt,ph,en,ip,el) + call nonlocal_dotState(Mp,subdt,ph,en) dotState = plasticState(ph)%dotState(:,en) end select plasticType @@ -352,21 +345,13 @@ end function plastic_dotState !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different plasticity constitutive models !-------------------------------------------------------------------------------------------------- -module subroutine plastic_dependentState(co, ip, el) +module subroutine plastic_dependentState(ph,en) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element - - integer :: & ph, & en - ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) - en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) - plasticType: select case (phase_plasticity(ph)) case (PLASTIC_DISLOTWIN_ID) plasticType @@ -376,7 +361,7 @@ module subroutine plastic_dependentState(co, ip, el) call dislotungsten_dependentState(ph,en) case (PLASTIC_NONLOCAL_ID) plasticType - call nonlocal_dependentState(ph,en,ip,el) + call nonlocal_dependentState(ph,en) end select plasticType @@ -422,7 +407,7 @@ module function plastic_deltaState(ph, en) result(broken) broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en))) if (.not. broken) then - mySize = plasticState(ph)%sizeDeltaState + mySize = plasticState(ph)%sizeDeltaState plasticState(ph)%deltaState2(1:mySize,en) = plasticState(ph)%deltaState2(1:mySize,en) & + plasticState(ph)%deltaState(1:mySize,en) end if diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index d668e40fe..dbc02b7a4 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -15,6 +15,9 @@ submodule(phase:plastic) nonlocal type :: tGeometry real(pReal), dimension(:), allocatable :: V_0 + integer, dimension(:,:,:), allocatable :: IPneighborhood + real(pReal), dimension(:,:), allocatable :: IParea, IPcoordinates + real(pReal), dimension(:,:,:), allocatable :: IPareaNormal end type tGeometry type(tGeometry), dimension(:), allocatable :: geom @@ -123,8 +126,10 @@ submodule(phase:plastic) nonlocal type :: tNonlocalDependentState real(pReal), allocatable, dimension(:,:) :: & - tau_pass, & - tau_Back + tau_pass, & + tau_Back + real(pReal), allocatable, dimension(:,:,:,:,:) :: & + compatibility end type tNonlocalDependentState type :: tNonlocalState @@ -160,7 +165,6 @@ submodule(phase:plastic) nonlocal state0 type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters - type(tNonlocalDependentState), dimension(:), allocatable :: dependentState contains @@ -406,6 +410,10 @@ module function plastic_nonlocal_init() result(myPlasticity) call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState,0) ! ToDo: state structure does not follow convention allocate(geom(ph)%V_0(Nmembers)) + allocate(geom(ph)%IPneighborhood(3,nIPneighbors,Nmembers)) + allocate(geom(ph)%IPareaNormal(3,nIPneighbors,Nmembers)) + allocate(geom(ph)%IParea(nIPneighbors,Nmembers)) + allocate(geom(ph)%IPcoordinates(3,Nmembers)) call storeGeometry(ph) if(plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) & @@ -489,6 +497,7 @@ module function plastic_nonlocal_init() result(myPlasticity) allocate(dst%tau_pass(prm%sum_N_sl,Nmembers),source=0.0_pReal) allocate(dst%tau_back(prm%sum_N_sl,Nmembers),source=0.0_pReal) + allocate(dst%compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,Nmembers),source=0.0_pReal) end associate if (Nmembers > 0) call stateInit(ini,ph,Nmembers) @@ -543,13 +552,11 @@ end function plastic_nonlocal_init !-------------------------------------------------------------------------------------------------- !> @brief calculates quantities characterizing the microstructure !-------------------------------------------------------------------------------------------------- -module subroutine nonlocal_dependentState(ph, en, ip, el) +module subroutine nonlocal_dependentState(ph, en) integer, intent(in) :: & ph, & - en, & - ip, & - el + en integer :: & no, & !< neighbor offset @@ -650,12 +657,12 @@ module subroutine nonlocal_dependentState(ph, en, ip, el) nRealNeighbors = 0.0_pReal neighbor_rhoTotal = 0.0_pReal do n = 1,nIPneighbors - neighbor_el = IPneighborhood(1,n,ip,el) - neighbor_ip = IPneighborhood(2,n,ip,el) - no = material_phasememberAt(1,neighbor_ip,neighbor_el) - if (neighbor_el > 0 .and. neighbor_ip > 0) then - if (material_phaseAt(1,neighbor_el) == ph) then + neighbor_el = geom(ph)%IPneighborhood(1,n,en) + neighbor_ip = geom(ph)%IPneighborhood(2,n,en) + if (neighbor_el > 0 .and. neighbor_ip > 0) then + if (material_phaseID(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip) == ph) then + no = material_phaseEntry(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip) nRealNeighbors = nRealNeighbors + 1.0_pReal rho_neighbor0 = getRho0(ph,no) @@ -665,11 +672,11 @@ module subroutine nonlocal_dependentState(ph, en, ip, el) neighbor_rhoTotal(1,:,n) = sum(abs(rho_neighbor0(:,edg)),2) neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor0(:,scr)),2) - connection_latticeConf(1:3,n) = matmul(invFe, discretization_IPcoords(1:3,neighbor_el+neighbor_ip-1) & - - discretization_IPcoords(1:3,el+neighbor_ip-1)) - normal_latticeConf = matmul(transpose(invFp), IPareaNormal(1:3,n,ip,el)) + connection_latticeConf(1:3,n) = matmul(invFe, geom(ph)%IPcoordinates(1:3,no) & + - geom(ph)%IPcoordinates(1:3,en)) + normal_latticeConf = matmul(transpose(invFp), geom(ph)%IPareaNormal(1:3,n,en)) if (math_inner(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) & ! neighboring connection points in opposite direction to face normal: must be periodic image - connection_latticeConf(1:3,n) = normal_latticeConf * IPvolume(ip,el)/IParea(n,ip,el) ! instead take the surface normal scaled with the diameter of the cell + connection_latticeConf(1:3,n) = normal_latticeConf * geom(ph)%V_0(en)/geom(ph)%IParea(n,en) ! instead take the surface normal scaled with the diameter of the cell else ! local neighbor or different lattice structure or different constitution instance -> use central values instead connection_latticeConf(1:3,n) = 0.0_pReal @@ -939,7 +946,7 @@ end subroutine plastic_nonlocal_deltaState !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- module subroutine nonlocal_dotState(Mp,timestep, & - ph,en,ip,el) + ph,en) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress @@ -947,9 +954,7 @@ module subroutine nonlocal_dotState(Mp,timestep, & timestep !< substepped crystallite time increment integer, intent(in) :: & ph, & - en, & - ip, & !< current integration point - el !< current element number + en integer :: & c, & !< character of dislocation @@ -1099,7 +1104,7 @@ module subroutine nonlocal_dotState(Mp,timestep, & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - rhoDot = rhoDotFlux(timestep, ph,en,ip,el) & + rhoDot = rhoDotFlux(timestep, ph,en) & + rhoDotMultiplication & + rhoDotSingle2DipoleGlide & + rhoDotAthermalAnnihilation & @@ -1110,7 +1115,7 @@ module subroutine nonlocal_dotState(Mp,timestep, & .or. any(rho(:,dip) + rhoDot(:,9:10) * timestep < -prm%atol_rho)) then #ifdef DEBUG if (debugConstitutive%extensive) then - print'(a,i5,a,i2)', '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip + print'(a,i5,a,i2)', '<< CONST >> evolution rate leads to negative density at ph ',ph,' en ',en print'(a)', '<< CONST >> enforcing cutback !!!' end if #endif @@ -1129,20 +1134,17 @@ end subroutine nonlocal_dotState !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- #if __INTEL_COMPILER >= 2020 -non_recursive function rhoDotFlux(timestep,ph,en,ip,el) +non_recursive function rhoDotFlux(timestep,ph,en) #else -function rhoDotFlux(timestep,ph,en,ip,el) +function rhoDotFlux(timestep,ph,en) #endif real(pReal), intent(in) :: & timestep !< substepped crystallite time increment integer, intent(in) :: & ph, & - en, & - ip, & !< current integration point - el !< current element number + en integer :: & - neighbor_ph, & !< phase of my neighbor's plasticity ns, & !< short notation for the total number of active slip systems c, & !< character of dislocation n, & !< index of my current neighbor @@ -1215,14 +1217,14 @@ function rhoDotFlux(timestep,ph,en,ip,el) !*** check CFL (Courant-Friedrichs-Lewy) condition for flux if (any( abs(dot_gamma) > 0.0_pReal & ! any active slip system ... .and. prm%C_CFL * abs(v0) * timestep & - > IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) + > geom(ph)%V_0(en)/ maxval(geom(ph)%IParea(:,en)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) #ifdef DEBUG if (debugConstitutive%extensive) then - print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip + print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at ph ',ph,' en ',en print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at ', & maxval(abs(v0), abs(dot_gamma) > 0.0_pReal & .and. prm%C_CFL * abs(v0) * timestep & - > IPvolume(ip,el) / maxval(IParea(:,ip,el))), & + > geom(ph)%V_0(en) / maxval(geom(ph)%IParea(:,en))), & ' at a timestep of ',timestep print*, '<< CONST >> enforcing cutback !!!' end if @@ -1245,19 +1247,18 @@ function rhoDotFlux(timestep,ph,en,ip,el) neighbors: do n = 1,nIPneighbors - neighbor_el = IPneighborhood(1,n,ip,el) - neighbor_ip = IPneighborhood(2,n,ip,el) - neighbor_n = IPneighborhood(3,n,ip,el) - np = material_phaseAt(1,neighbor_el) - no = material_phasememberAt(1,neighbor_ip,neighbor_el) + neighbor_el = geom(ph)%IPneighborhood(1,n,en) + neighbor_ip = geom(ph)%IPneighborhood(2,n,en) + neighbor_n = geom(ph)%IPneighborhood(3,n,en) + np = material_phaseID(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip) + no = material_phaseEntry(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip) opposite_neighbor = n + mod(n,2) - mod(n+1,2) - opposite_el = IPneighborhood(1,opposite_neighbor,ip,el) - opposite_ip = IPneighborhood(2,opposite_neighbor,ip,el) - opposite_n = IPneighborhood(3,opposite_neighbor,ip,el) + opposite_el = geom(ph)%IPneighborhood(1,opposite_neighbor,en) + opposite_ip = geom(ph)%IPneighborhood(2,opposite_neighbor,en) + opposite_n = geom(ph)%IPneighborhood(3,opposite_neighbor,en) if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient - neighbor_ph = material_phaseAt(1,neighbor_el) neighbor_F = phase_mechanical_F(np)%data(1:3,1:3,no) neighbor_Fe = matmul(neighbor_F, math_inv33(phase_mechanical_Fp(np)%data(1:3,1:3,no))) Favg = 0.5_pReal * (my_F + neighbor_F) @@ -1276,12 +1277,12 @@ function rhoDotFlux(timestep,ph,en,ip,el) !* The entering flux from my neighbor will be distributed on my slip systems according to the !* compatibility if (neighbor_n > 0) then - if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTIC_NONLOCAL_ID .and. & - any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) then + if (phase_plasticity(np) == PLASTIC_NONLOCAL_ID .and. & + any(dependentState(ph)%compatibility(:,:,:,n,en) > 0.0_pReal)) then forall (s = 1:ns, t = 1:4) - neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,neighbor_ph),no) - neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_ph),no),0.0_pReal) + neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,np),no) + neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,np),no),0.0_pReal) endforall where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%rho_min & @@ -1301,12 +1302,12 @@ function rhoDotFlux(timestep,ph,en,ip,el) .and. v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density lineLength = neighbor_rhoSgl0(s,t) * neighbor_v0(s,t) & * math_inner(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface - where (compatibility(c,:,s,n,ip,el) > 0.0_pReal) & + where (dependentState(ph)%compatibility(c,:,s,n,en) > 0.0_pReal) & rhoDotFlux(:,t) = rhoDotFlux(1:ns,t) & - + lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2 ! transferring to equally signed mobile dislocation type - where (compatibility(c,:,s,n,ip,el) < 0.0_pReal) & + + lineLength/geom(ph)%V_0(en)*dependentState(ph)%compatibility(c,:,s,n,en)**2 ! transferring to equally signed mobile dislocation type + where (dependentState(ph)%compatibility(c,:,s,n,en) < 0.0_pReal) & rhoDotFlux(:,topp) = rhoDotFlux(:,topp) & - + lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2 ! transferring to opposite signed mobile dislocation type + + lineLength/geom(ph)%V_0(en)*dependentState(ph)%compatibility(c,:,s,n,en)**2 ! transferring to opposite signed mobile dislocation type end if end do @@ -1322,28 +1323,28 @@ function rhoDotFlux(timestep,ph,en,ip,el) !* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density. !* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations. if (opposite_n > 0) then - if (phase_plasticity(material_phaseAt(1,opposite_el)) == PLASTIC_NONLOCAL_ID) then + if (phase_plasticity(np) == PLASTIC_NONLOCAL_ID) then normal_me2neighbor_defConf = math_det33(Favg) & - * matmul(math_inv33(transpose(Favg)),IPareaNormal(1:3,n,ip,el)) ! normal of the interface in (average) deformed configuration (pointing en => neighbor) + * matmul(math_inv33(transpose(Favg)),geom(ph)%IPareaNormal(1:3,n,en)) ! normal of the interface in (average) deformed configuration (pointing en => neighbor) normal_me2neighbor = matmul(transpose(my_Fe), normal_me2neighbor_defConf) & / math_det33(my_Fe) ! interface normal in my lattice configuration - area = IParea(n,ip,el) * norm2(normal_me2neighbor) + area = geom(ph)%IParea(n,en) * norm2(normal_me2neighbor) normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor) ! normalize the surface normal to unit length do s = 1,ns do t = 1,4 c = (t + 1) / 2 if (v0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from en to my neighbor == leaving flux for en (might also be a pure flux from my mobile density to dead density if interface not at all transmissive) if (v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal) then ! no sign change in flux density - transmissivity = sum(compatibility(c,:,s,n,ip,el)**2) ! overall transmissivity from this slip system to my neighbor + transmissivity = sum(dependentState(ph)%compatibility(c,:,s,n,en)**2) ! overall transmissivity from this slip system to my neighbor else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor transmissivity = 0.0_pReal end if lineLength = my_rhoSgl0(s,t) * v0(s,t) & * math_inner(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface - rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / IPvolume(ip,el) ! subtract dislocation flux from current type + rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / geom(ph)%V_0(en) ! subtract dislocation flux from current type rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) & - + lineLength / IPvolume(ip,el) * (1.0_pReal - transmissivity) & + + lineLength / geom(ph)%V_0(en) * (1.0_pReal - transmissivity) & * sign(1.0_pReal, v0(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point end if end do @@ -1364,14 +1365,14 @@ end function rhoDotFlux ! plane normals and signed cosine of the angle between the slip directions. Only the largest values ! that sum up to a total of 1 are considered, all others are set to zero. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) +module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el) type(tRotationContainer), dimension(:), intent(in) :: & orientation ! crystal orientation integer, intent(in) :: & ph, & - i, & - e + ip, & + el integer :: & n, & ! neighbor index @@ -1397,16 +1398,16 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) associate(prm => param(ph)) ns = prm%sum_N_sl - en = material_phaseMemberAt(1,i,e) + en = material_phaseEntry(1,(el-1)*discretization_nIPs + ip) !*** start out fully compatible my_compatibility = 0.0_pReal forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pReal neighbors: do n = 1,nIPneighbors - neighbor_e = IPneighborhood(1,n,i,e) - neighbor_i = IPneighborhood(2,n,i,e) - neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e) - neighbor_phase = material_phaseAt(1,neighbor_e) + neighbor_e = IPneighborhood(1,n,ip,el) + neighbor_i = IPneighborhood(2,n,ip,el) + neighbor_me = material_phaseEntry(1,(neighbor_e-1)*discretization_nIPs + neighbor_i) + neighbor_phase = material_phaseID(1,(neighbor_e-1)*discretization_nIPs + neighbor_i) if (neighbor_e <= 0 .or. neighbor_i <= 0) then !* FREE SURFACE @@ -1465,7 +1466,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) end do neighbors - compatibility(:,:,:,:,i,e) = my_compatibility + dependentState(ph)%compatibility(:,:,:,:,material_phaseEntry(1,(el-1)*discretization_nIPs + ip)) = my_compatibility end associate @@ -1756,14 +1757,29 @@ subroutine storeGeometry(ph) integer, intent(in) :: ph - integer :: ce, co + integer :: ce, co, nCell real(pReal), dimension(:), allocatable :: V + integer, dimension(:,:,:), allocatable :: neighborhood + real(pReal), dimension(:,:), allocatable :: area, coords + real(pReal), dimension(:,:,:), allocatable :: areaNormal + nCell = product(shape(IPvolume)) + + V = reshape(IPvolume,[nCell]) + neighborhood = reshape(IPneighborhood,[3,nIPneighbors,nCell]) + area = reshape(IParea,[nIPneighbors,nCell]) + areaNormal = reshape(IPareaNormal,[3,nIPneighbors,nCell]) + coords = reshape(discretization_IPcoords,[3,nCell]) - V = reshape(IPvolume,[product(shape(IPvolume))]) do ce = 1, size(material_homogenizationEntry,1) do co = 1, homogenization_maxNconstituents - if (material_phaseID(co,ce) == ph) geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce) + if (material_phaseID(co,ce) == ph) then + geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce) + geom(ph)%IPneighborhood(:,:,material_phaseEntry(co,ce)) = neighborhood(:,:,ce) + geom(ph)%IParea(:,material_phaseEntry(co,ce)) = area(:,ce) + geom(ph)%IPareaNormal(:,:,material_phaseEntry(co,ce)) = areaNormal(:,:,ce) + geom(ph)%IPcoordinates(:,material_phaseEntry(co,ce)) = coords(:,ce) + end if end do end do