diff --git a/src/phase.f90 b/src/phase.f90 index bea1a5543..b7bace81a 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -408,7 +408,7 @@ end subroutine phase_init !> @brief Allocate the components of the state structure for a given phase !-------------------------------------------------------------------------------------------------- subroutine phase_allocateState(state, & - NEntries,sizeState,sizeDotState,sizeDeltaState) + NEntries,sizeState,sizeDotState,sizeDeltaState,offsetDeltaState) class(tState), intent(inout) :: & state @@ -417,12 +417,17 @@ subroutine phase_allocateState(state, & sizeState, & sizeDotState, & sizeDeltaState - + integer, intent(in), optional :: & + offsetDeltaState state%sizeState = sizeState state%sizeDotState = sizeDotState state%sizeDeltaState = sizeDeltaState - state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition + if (present(offsetDeltaState)) then + state%offsetDeltaState = offsetDeltaState ! ToDo: this is a fix for broken nonlocal + else + state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition + end if allocate(state%atol (sizeState), source=0.0_pReal) allocate(state%state0 (sizeState,NEntries), source=0.0_pReal) @@ -432,7 +437,6 @@ subroutine phase_allocateState(state, & allocate(state%deltaState (sizeDeltaState,NEntries), source=0.0_pReal) - end subroutine phase_allocateState diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 234bee8e9..04cc4b946 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -88,7 +88,7 @@ submodule(phase) mechanical en real(pReal), intent(in) :: & subdt !< timestep - real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & dotState end function plastic_dotState diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 6d2aac827..912aadc03 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -154,11 +154,10 @@ submodule(phase:mechanical) plastic en end subroutine dislotungsten_dotState - module subroutine nonlocal_dotState(Mp,Temperature,timestep,ph,en,ip,el) + module subroutine nonlocal_dotState(Mp,timestep,ph,en,ip,el) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress real(pReal), intent(in) :: & - Temperature, & !< temperature timestep !< substepped crystallite time increment integer, intent(in) :: & ph, & @@ -308,7 +307,7 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState) subdt !< timestep real(pReal), dimension(3,3) :: & Mp - real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & dotState @@ -334,7 +333,7 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState) call dislotungsten_dotState(Mp,thermal_T(ph,en),ph,en) case (PLASTIC_NONLOCAL_ID) plasticType - call nonlocal_dotState(Mp,thermal_T(ph,en),subdt,ph,en,ip,el) + call nonlocal_dotState(Mp,subdt,ph,en,ip,el) end select plasticType end if @@ -394,6 +393,7 @@ module function plastic_deltaState(ph, en) result(broken) myOffset, & mySize + broken = .false. select case (phase_plasticity(ph)) diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 0dc3dc4bd..a8f56b5c5 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -403,7 +403,7 @@ module function plastic_nonlocal_init() result(myPlasticity) 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure sizeDeltaState = sizeDotState - call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState) + call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState,0) ! ToDo: state structure does not follow convention allocate(geom(ph)%V_0(Nmembers)) call storeGeometry(ph) @@ -411,9 +411,6 @@ module function plastic_nonlocal_init() result(myPlasticity) if(plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) & call IO_error(212,ext_msg='IPneighborhood does not exist') - - plasticState(ph)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention - st0%rho => plasticState(ph)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) stt%rho => plasticState(ph)%state (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) dot%rho => plasticState(ph)%dotState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) @@ -941,13 +938,12 @@ end subroutine plastic_nonlocal_deltaState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- -module subroutine nonlocal_dotState(Mp, Temperature,timestep, & +module subroutine nonlocal_dotState(Mp,timestep, & ph,en,ip,el) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress real(pReal), intent(in) :: & - Temperature, & !< temperature timestep !< substepped crystallite time increment integer, intent(in) :: & ph, & @@ -984,7 +980,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & real(pReal) :: & D_SD, & mu, & - nu + nu, Temperature if (timestep <= 0.0_pReal) then plasticState(ph)%dotState = 0.0_pReal @@ -995,6 +991,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & mu = elastic_mu(ph,en) nu = elastic_nu(ph,en) + Temperature = thermal_T(ph,en) tau = 0.0_pReal dot_gamma = 0.0_pReal @@ -1195,7 +1192,6 @@ function rhoDotFlux(timestep,ph,en,ip,el) associate(prm => param(ph), & dst => dependentState(ph), & - dot => dotState(ph), & stt => state(ph)) ns = prm%sum_N_sl @@ -1206,7 +1202,7 @@ function rhoDotFlux(timestep,ph,en,ip,el) rho0 = getRho0(ph,en) my_rhoSgl0 = rho0(:,sgl) - forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) !ToDo: MD: I think we should use state0 here + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) !ToDo: MD: I think we should use state0 here dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en) @@ -1217,7 +1213,7 @@ function rhoDotFlux(timestep,ph,en,ip,el) if (plasticState(ph)%nonlocal) then !*** check CFL (Courant-Friedrichs-Lewy) condition for flux - if (any( abs(dot_gamma) > 0.0_pReal & ! any active slip system ... + 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) #ifdef DEBUG