From 836172028d4f45201139a9b5d4eb9d5d70d14c5c Mon Sep 17 00:00:00 2001 From: Eureka Pai Kulyadi Date: Fri, 6 Sep 2019 15:44:38 -0400 Subject: [PATCH 1/3] modified plastic_isotropic_LiAndItsTangent to behave correctly --- src/plastic_isotropic.f90 | 47 +++++++++++++++++++++------------------ 1 file changed, 25 insertions(+), 22 deletions(-) diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90 index 3b5ec6e5a..5c5d134f5 100644 --- a/src/plastic_isotropic.f90 +++ b/src/plastic_isotropic.f90 @@ -283,49 +283,52 @@ end subroutine plastic_isotropic_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent -! ToDo: Rename Tstar to Mi? !-------------------------------------------------------------------------------------------------- -subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar,Tstar,instance,of) +subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) real(pReal), dimension(3,3), intent(out) :: & Li !< inleastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & - dLi_dTstar !< derivative of Li with respect to the Mandel stress + dLi_dMi !< derivative of Li with respect to the Mandel stress real(pReal), dimension(3,3), intent(in) :: & - Tstar !< Mandel stress ToDo: Mi? + Mi !< Mandel stress integer, intent(in) :: & instance, & of real(pReal), dimension(3,3) :: & - Tstar_sph !< sphiatoric part of the Mandel stress + Mi_sph !< spherical part of the Mandel stress real(pReal) :: & - dot_gamma, & !< strainrate - norm_Tstar_sph, & !< euclidean norm of Tstar_sph - squarenorm_Tstar_sph !< square of the euclidean norm of Tstar_sph + dot_gamma, & !< shear rate + tr !< pressure integer :: & k, l, m, n associate(prm => param(instance), stt => state(instance)) - Tstar_sph = math_spherical33(Tstar) - squarenorm_Tstar_sph = math_mul33xx33(Tstar_sph,Tstar_sph) - norm_Tstar_sph = sqrt(squarenorm_Tstar_sph) + tr=math_trace33(math_spherical33(Mi)) + + if (prm%dilatation .and. abs(tr) > 0.0_pReal) then ! no stress or J2 plasticity --> Li and its derivative are zero + Li = math_I3 & + * prm%dot_gamma_0/prm%M * (3.0_pReal*prm%M*stt%xi(of))**(-prm%n) & + * tr * abs(tr)**(prm%n-1.0_pReal) - if (prm%dilatation .and. norm_Tstar_sph > 0.0_pReal) then ! no stress or J2 plastitiy --> Li and its derivative are zero - dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Tstar_sph /(prm%M*stt%xi(of))) **prm%n - - Li = math_I3/sqrt(3.0_pReal) * dot_gamma/prm%M +#ifdef DEBUG + if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 & + .and. (of == prm%of_debug .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0)) then + write(6,'(/,a,/,f12.5)') '<< CONST isotropic >> pressure / MPa', tr/3.0_pReal*1.0e-6_pReal + write(6,'(/,a,/,f12.5)') '<< CONST isotropic >> gdot', prm%dot_gamma_0 * (3.0_pReal*prm%M*stt%xi(of))**(-prm%n) & + * tr * abs(tr)**(prm%n-1.0_pReal) + end if +#endif + forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLi_dTstar(k,l,m,n) = (prm%n-1.0_pReal) * Tstar_sph(k,l)*Tstar_sph(m,n) / squarenorm_Tstar_sph - forall (k=1:3,l=1:3) & - dLi_dTstar(k,l,k,l) = dLi_dTstar(k,l,k,l) + 1.0_pReal - - dLi_dTstar = dot_gamma / prm%M * dLi_dTstar / norm_Tstar_sph + dLi_dMi(k,l,m,n) = n / tr * Li(k,l) * math_I3(m,n) + else - Li = 0.0_pReal - dLi_dTstar = 0.0_pReal + Li = 0.0_pReal + dLi_dMi = 0.0_pReal endif end associate From cf37f8d405ad933231f01acd13224f4608577515 Mon Sep 17 00:00:00 2001 From: Eureka Pai Kulyadi Date: Fri, 6 Sep 2019 15:45:49 -0400 Subject: [PATCH 2/3] added some lines to help while debugging --- src/crystallite.f90 | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 246cb7092..0157db0e0 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1274,7 +1274,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then - write(6,'(a,i3,/)') '<< CRYST integrateStress >> iteration ', NiterationStressLp + write(6,'(a,i3,/)') '<< CRYST integrateStress >> Lp iteration ', NiterationStressLp write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST integrateStress >> Lpguess', transpose(Lpguess) write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST integrateStress >> Lp_constitutive', transpose(Lp_constitutive) write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST integrateStress >> Fi', transpose(Fi_new) @@ -1377,6 +1377,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then + write(6,'(a,i3,/)') '<< CRYST integrateStress >> Li iteration ', NiterationStressLi write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> Li_constitutive', transpose(Li_constitutive) write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> Liguess', transpose(Liguess) endif @@ -1405,6 +1406,13 @@ logical function integrateStress(ipc,ip,el,timeFraction) else ! not converged and residuum not improved... steplengthLi = num%subStepSizeLi * steplengthLi ! ...try with smaller step length in same direction Liguess = Liguess_old + steplengthLi * deltaLi +#ifdef DEBUG + if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & + .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & + .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then + write(6,'(a,1x,f7.4)') '<< CRYST integrateStress >> linear search for Liguess with step', steplengthLi + endif +#endif cycle LiLoop endif @@ -1450,6 +1458,13 @@ logical function integrateStress(ipc,ip,el,timeFraction) jacoCounterLi = jacoCounterLi + 1 Liguess = Liguess + steplengthLi * deltaLi +#ifdef DEBUG + if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & + .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & + .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then + write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> corrected Liguess by', transpose(deltaLi) + endif +#endif enddo LiLoop !* calculate new plastic and elastic deformation gradient From ca760f9ca18828731b6a1e4bf7cef821a9cbd8e2 Mon Sep 17 00:00:00 2001 From: Eureka Pai Kulyadi Date: Fri, 6 Sep 2019 15:47:02 -0400 Subject: [PATCH 3/3] added text lines --- src/grid/DAMASK_grid.f90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index edf26af19..0bec7606c 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -118,6 +118,7 @@ program DAMASK_spectral ! assign mechanics solver depending on selected type select case (trim(config_numerics%getString('spectral_solver',defaultVal='basic'))) case ('basic') + write(6,'(a)') ' Basic solution algorithm' mech_init => grid_mech_spectral_basic_init mech_forward => grid_mech_spectral_basic_forward mech_solution => grid_mech_spectral_basic_solution @@ -125,6 +126,7 @@ program DAMASK_spectral case ('polarisation') if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) & call IO_warning(42, ext_msg='debug Divergence') + write(6,'(a)') ' Polarisation solution algorithm' mech_init => grid_mech_spectral_polarisation_init mech_forward => grid_mech_spectral_polarisation_forward mech_solution => grid_mech_spectral_polarisation_solution @@ -132,6 +134,7 @@ program DAMASK_spectral case ('fem') if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) & call IO_warning(42, ext_msg='debug Divergence') + write(6,'(a)') ' FEM solution algorithm' mech_init => grid_mech_FEM_init mech_forward => grid_mech_FEM_forward mech_solution => grid_mech_FEM_solution