diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 1abb04377..d64b44f6d 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -636,7 +636,7 @@ select case (phase_constitution(material_phase(ipc,ip,el))) case (constitutive_nonlocal_label) call constitutive_nonlocal_dotState(constitutive_dotState(ipc,ip,el), Tstar_v, Fe, Fp, Temperature, constitutive_state, & - constitutive_aTolState, subdt, orientation, ipc, ip, el) + subdt, orientation, ipc, ip, el) end select diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index 6f8a00c88..3f6ba1c89 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -1516,7 +1516,7 @@ endsubroutine !********************************************************************* !* rate of change of microstructure * !********************************************************************* -subroutine constitutive_nonlocal_dotState(dotState, Tstar_v, Fe, Fp, Temperature, state, aTolState, timestep, orientation, g,ip,el) +subroutine constitutive_nonlocal_dotState(dotState, Tstar_v, Fe, Fp, Temperature, state, timestep, orientation, g,ip,el) use prec, only: pReal, & pInt, & @@ -1579,8 +1579,7 @@ real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & orientation ! crystal lattice orientation type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & - state, & ! current microstructural state - aTolState ! absolute state tolerance + state ! current microstructural state !*** input/output variables type(p_vec), intent(inout) :: dotState ! evolution of state variables / microstructure @@ -1972,6 +1971,7 @@ rhoDotThermalAnnihilation(1:ns,10) = 0.0_pReal !**************************************************************************** !*** assign the rates of dislocation densities to my dotState +!*** if evolution rates lead to negative densities, a cutback is enforced rhoDot = 0.0_pReal rhoDot = rhoDotFlux & @@ -1981,7 +1981,21 @@ rhoDot = rhoDotFlux & + rhoDotAthermalAnnihilation & + rhoDotThermalAnnihilation -dotState%p(1:10*ns) = dotState%p(1:10*ns) + reshape(rhoDot,(/10*ns/)) +if ( any(rhoSgl(1:ns,1:4) + rhoDot(1:ns,1:4) * timestep < - constitutive_nonlocal_aTolRho(myInstance)) & + .or. any(rhoDip(1:ns,1:2) + rhoDot(1:ns,9:10) * timestep < - constitutive_nonlocal_aTolRho(myInstance))) then +#ifndef _OPENMP + if (debug_verbosity > 6) then + write(6,'(a,i5,a,i2)') '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip + write(6,'(a)') '<< CONST >> enforcing cutback !!!' + endif +#endif + dotState%p = DAMASK_NaN + return +else + dotState%p(1:10*ns) = dotState%p(1:10*ns) + reshape(rhoDot,(/10*ns/)) +endif + + #ifndef _OPENMP if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then