diff --git a/.gitmodules b/.gitmodules index 5d17eb0cb..0587fff4c 100644 --- a/.gitmodules +++ b/.gitmodules @@ -2,3 +2,4 @@ path = PRIVATE url = ../PRIVATE.git branch = master + shallow = true diff --git a/PRIVATE b/PRIVATE index 313dd5de6..45ef93dbf 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 313dd5de618c996cdf9ace95a096f25e757386d9 +Subproject commit 45ef93dbfa3e0e6fa830914b3632e188c308a099 diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index a19a70432..fe8c7d1b3 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -13,7 +13,6 @@ module CPFEM use discretization_marc use material use config - use crystallite use homogenization use IO use discretization @@ -262,6 +261,7 @@ end subroutine CPFEM_general subroutine CPFEM_forward call crystallite_forward + call constitutive_forward end subroutine CPFEM_forward diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 54e381d34..636962948 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -21,7 +21,6 @@ module CPFEM2 use HDF5_utilities use homogenization use constitutive - use crystallite #if defined(Mesh) use FEM_quadrature use discretization_mesh @@ -99,6 +98,7 @@ end subroutine CPFEM_restartWrite subroutine CPFEM_forward call crystallite_forward + call constitutive_forward end subroutine CPFEM_forward diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index beabfcae1..08e7b9c1c 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -44,12 +44,9 @@ #include "source_damage_anisoDuctile.f90" #include "kinematics_cleavage_opening.f90" #include "kinematics_slipplane_opening.f90" -#include "crystallite.f90" #include "thermal_isothermal.f90" -#include "thermal_adiabatic.f90" #include "thermal_conduction.f90" #include "damage_none.f90" -#include "damage_local.f90" #include "damage_nonlocal.f90" #include "homogenization.f90" #include "homogenization_mech.f90" diff --git a/src/constitutive.f90 b/src/constitutive.f90 index a60352f81..f76b94579 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -13,11 +13,103 @@ module constitutive use results use lattice use discretization - use geometry_plastic_nonlocal, only: & - geometry_plastic_nonlocal_disable + use parallelization + use HDF5_utilities + use DAMASK_interface + use FEsolving + use results implicit none private + real(pReal), dimension(:,:,:), allocatable, public :: & + crystallite_dt !< requested time increment of each grain + real(pReal), dimension(:,:,:), allocatable :: & + crystallite_subdt, & !< substepped time increment of each grain + crystallite_subStep !< size of next integration step + type(rotation), dimension(:,:,:), allocatable :: & + crystallite_orientation !< current orientation + real(pReal), dimension(:,:,:,:,:), allocatable :: & + crystallite_F0, & !< def grad at start of FE inc + crystallite_subF, & !< def grad to be reached at end of crystallite inc + crystallite_subF0, & !< def grad at start of crystallite inc + ! + crystallite_Fe, & !< current "elastic" def grad (end of converged time step) + ! + crystallite_Fp, & !< current plastic def grad (end of converged time step) + crystallite_Fp0, & !< plastic def grad at start of FE inc + crystallite_partitionedFp0,& !< plastic def grad at start of homog inc + crystallite_subFp0,& !< plastic def grad at start of crystallite inc + ! + crystallite_subFi0,& !< intermediate def grad at start of crystallite inc + ! + crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc + crystallite_partitionedLp0, & !< plastic velocity grad at start of homog inc + ! + crystallite_S0, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc + crystallite_partitionedS0 !< 2nd Piola-Kirchhoff stress vector at start of homog inc + real(pReal), dimension(:,:,:,:,:), allocatable, public :: & + crystallite_P, & !< 1st Piola-Kirchhoff stress per grain + crystallite_Lp, & !< current plastic velocitiy grad (end of converged time step) + crystallite_S, & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step) + crystallite_partitionedF0 !< def grad at start of homog inc + real(pReal), dimension(:,:,:,:,:), allocatable, public :: & + crystallite_partitionedF !< def grad to be reached at end of homog inc + + logical, dimension(:,:,:), allocatable, public :: & + crystallite_requested !< used by upper level (homogenization) to request crystallite calculation + logical, dimension(:,:,:), allocatable :: & + crystallite_converged !< convergence flag + + type :: tOutput !< new requested output (per phase) + character(len=pStringLen), allocatable, dimension(:) :: & + label + end type tOutput + type(tOutput), allocatable, dimension(:) :: output_constituent + + type :: tTensorContainer + real(pReal), dimension(:,:,:), allocatable :: data + end type + + type(tTensorContainer), dimension(:), allocatable :: & + constitutive_mech_Fi, & + constitutive_mech_Fi0, & + constitutive_mech_partionedFi0, & + constitutive_mech_Li, & + constitutive_mech_Li0, & + constitutive_mech_partionedLi0 + + type :: tNumerics + integer :: & + iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp + nState, & !< state loop limit + nStress !< stress loop limit + real(pReal) :: & + subStepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback + subStepSizeCryst, & !< size of first substep when cutback + subStepSizeLp, & !< size of first substep when cutback in Lp calculation + subStepSizeLi, & !< size of first substep when cutback in Li calculation + stepIncreaseCryst, & !< increase of next substep size when previous substep converged + rtol_crystalliteState, & !< relative tolerance in state loop + rtol_crystalliteStress, & !< relative tolerance in stress loop + atol_crystalliteStress !< absolute tolerance in stress loop + end type tNumerics + + type(tNumerics) :: num ! numerics parameters. Better name? + + type :: tDebugOptions + logical :: & + basic, & + extensive, & + selective + integer :: & + element, & + ip, & + grain + end type tDebugOptions + + type(tDebugOptions) :: debugCrystallite + + procedure(integrateStateFPI), pointer :: integrateState integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: & phase_plasticity !< plasticity of each phase @@ -57,6 +149,45 @@ module constitutive module subroutine thermal_init end subroutine thermal_init + module function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el, & !< element + phase, & + of + real(pReal), intent(in) :: & + subdt !< timestep + real(pReal), intent(in), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems) :: & + FArray, & !< elastic deformation gradient + FpArray !< plastic deformation gradient + real(pReal), intent(in), dimension(3,3) :: & + Fi !< intermediate deformation gradient + real(pReal), intent(in), dimension(3,3) :: & + S !< 2nd Piola Kirchhoff stress (vector notation) + + logical :: broken +end function constitutive_collectDotState + + +module function constitutive_deltaState(S, Fi, ipc, ip, el, phase, of) result(broken) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el, & !< element + phase, & + of + real(pReal), intent(in), dimension(3,3) :: & + S, & !< 2nd Piola Kirchhoff stress + Fi !< intermediate deformation gradient + logical :: & + broken + + +end function constitutive_deltaState + module function plastic_active(plastic_label) result(active_plastic) character(len=*), intent(in) :: plastic_label @@ -75,67 +206,6 @@ module constitutive logical, dimension(:,:), allocatable :: active_kinematics end function kinematics_active - module subroutine plastic_isotropic_dotState(Mp,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_isotropic_dotState - - module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_phenopowerlaw_dotState - - module subroutine plastic_kinehardening_dotState(Mp,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_kinehardening_dotState - - module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T - integer, intent(in) :: & - instance, & - of - end subroutine plastic_dislotwin_dotState - - module subroutine plastic_disloTungsten_dotState(Mp,T,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T - integer, intent(in) :: & - instance, & - of - end subroutine plastic_disloTungsten_dotState - - module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & - instance,of,ip,el) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< MandelStress - real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: & - F, & !< deformation gradient - Fp !< plastic deformation gradient - real(pReal), intent(in) :: & - Temperature, & !< temperature - timestep !< substepped crystallite time increment - integer, intent(in) :: & - instance, & - of, & - ip, & !< current integration point - el !< current element number - end subroutine plastic_nonlocal_dotState - module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) integer, intent(in) :: & @@ -200,15 +270,6 @@ module constitutive el !< element end function plastic_dislotwin_homogenizedC - pure module function kinematics_thermal_expansion_initialStrain(homog,phase,offset) result(initialStrain) - integer, intent(in) :: & - phase, & - homog, & - offset - real(pReal), dimension(3,3) :: & - initialStrain - end function kinematics_thermal_expansion_initialStrain - module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) integer, intent(in) :: & instance, & @@ -268,24 +329,6 @@ module constitutive end subroutine kinematics_thermal_expansion_LiAndItsTangent - module subroutine plastic_kinehardening_deltaState(Mp,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_kinehardening_deltaState - - module subroutine plastic_nonlocal_deltaState(Mp,instance,of,ip,el) - real(pReal), dimension(3,3), intent(in) :: & - Mp - integer, intent(in) :: & - instance, & - of, & - ip, & - el - end subroutine plastic_nonlocal_deltaState - module subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -303,9 +346,7 @@ module constitutive module subroutine damage_results end subroutine damage_results - end interface - interface constitutive_LpAndItsTangents module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & S, Fi, ipc, ip, el) @@ -323,9 +364,6 @@ module constitutive dLp_dFi !< derivative of Lp with respect to Fi end subroutine constitutive_plastic_LpAndItsTangents - end interface constitutive_LpAndItsTangents - - interface constitutive_dependentState module subroutine constitutive_plastic_dependentState(F, Fp, ipc, ip, el) integer, intent(in) :: & @@ -337,10 +375,8 @@ module constitutive Fp !< plastic deformation gradient end subroutine constitutive_plastic_dependentState - end interface constitutive_dependentState - interface constitutive_SandItsTangents - + module subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -356,40 +392,63 @@ module constitutive dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient end subroutine constitutive_hooke_SandItsTangents - end interface constitutive_SandItsTangents + module subroutine integrateStateFPI(g,i,e) + integer, intent(in) :: e, i, g + end subroutine integrateStateFPI + module subroutine integrateStateEuler(g,i,e) + integer, intent(in) :: e, i, g + end subroutine integrateStateEuler + + module subroutine integrateStateAdaptiveEuler(g,i,e) + integer, intent(in) :: e, i, g + end subroutine integrateStateAdaptiveEuler + + module subroutine integrateStateRK4(g,i,e) + integer, intent(in) :: e, i, g + end subroutine integrateStateRK4 + + module subroutine integrateStateRKCK45(g,i,e) + integer, intent(in) :: e, i, g + end subroutine integrateStateRKCK45 + + end interface - type :: tDebugOptions - logical :: & - basic, & - extensive, & - selective - integer :: & - element, & - ip, & - grain - end type tDebugOptions type(tDebugOptions) :: debugConstitutive public :: & constitutive_init, & constitutive_homogenizedC, & - constitutive_LpAndItsTangents, & - constitutive_dependentState, & constitutive_LiAndItsTangents, & - constitutive_initialFi, & - constitutive_SandItsTangents, & constitutive_collectDotState, & + constitutive_collectDotState_source, & constitutive_deltaState, & + constitutive_deltaState_source, & constitutive_damage_getRateAndItsTangents, & constitutive_thermal_getRateAndItsTangents, & constitutive_results, & constitutive_allocateState, & + constitutive_forward, & + constitutive_restore, & plastic_nonlocal_updateCompatibility, & plastic_active, & source_active, & - kinematics_active + kinematics_active, & + converged, & + crystallite_init, & + crystallite_stress, & + crystallite_stressTangent, & + crystallite_orientations, & + crystallite_push33ToRef, & + crystallite_results, & + crystallite_restartWrite, & + crystallite_restartRead, & + crystallite_forward, & + crystallite_initializeRestorationPoints, & + crystallite_windForward, & + crystallite_restore + contains @@ -421,7 +480,7 @@ subroutine constitutive_init print'(/,a)', ' <<<+- constitutive init -+>>>'; flush(IO_STDOUT) - phases => config_material%get('phase') + phases => config_material%get('phase') constitutive_source_maxSizeDotState = 0 PhaseLoop2:do p = 1,phases%length !-------------------------------------------------------------------------------------------------- @@ -607,43 +666,13 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & end subroutine constitutive_LiAndItsTangents -!-------------------------------------------------------------------------------------------------- -!> @brief collects initial intermediate deformation gradient -!-------------------------------------------------------------------------------------------------- -pure function constitutive_initialFi(ipc, ip, el) - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), dimension(3,3) :: & - constitutive_initialFi !< composite initial intermediate deformation gradient - integer :: & - k !< counter in kinematics loop - integer :: & - phase, & - homog, offset - - constitutive_initialFi = math_I3 - phase = material_phaseAt(ipc,el) - - KinematicsLoop: do k = 1, phase_Nkinematics(phase) !< Warning: small initial strain assumption - kinematicsType: select case (phase_kinematics(k,phase)) - case (KINEMATICS_thermal_expansion_ID) kinematicsType - homog = material_homogenizationAt(el) - offset = material_homogenizationMemberAt(ip,el) - constitutive_initialFi = constitutive_initialFi & - + kinematics_thermal_expansion_initialStrain(homog,phase,offset) - end select kinematicsType - enddo KinematicsLoop - -end function constitutive_initialFi !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken) +function constitutive_collectDotState_source(S, ipc, ip, el,phase,of) result(broken) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -651,52 +680,14 @@ function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el el, & !< element phase, & of - real(pReal), intent(in) :: & - subdt !< timestep - real(pReal), intent(in), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems) :: & - FArray, & !< elastic deformation gradient - FpArray !< plastic deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - Fi !< intermediate deformation gradient real(pReal), intent(in), dimension(3,3) :: & S !< 2nd Piola Kirchhoff stress (vector notation) - real(pReal), dimension(3,3) :: & - Mp integer :: & - ho, & !< homogenization - tme, & !< thermal member position - i, & !< counter in source loop - instance + i !< counter in source loop logical :: broken - ho = material_homogenizationAt(el) - tme = material_homogenizationMemberAt(ip,el) - instance = phase_plasticityInstance(phase) - Mp = matmul(matmul(transpose(Fi),Fi),S) - - plasticityType: select case (phase_plasticity(phase)) - - case (PLASTICITY_ISOTROPIC_ID) plasticityType - call plastic_isotropic_dotState(Mp,instance,of) - - case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call plastic_phenopowerlaw_dotState(Mp,instance,of) - - case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_dotState(Mp,instance,of) - - case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_dotState(Mp,temperature(ho)%p(tme),instance,of) - - case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType - call plastic_disloTungsten_dotState(Mp,temperature(ho)%p(tme),instance,of) - - case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dotState(Mp,FArray,FpArray,temperature(ho)%p(tme),subdt, & - instance,of,ip,el) - end select plasticityType - broken = any(IEEE_is_NaN(plasticState(phase)%dotState(:,of))) + broken = .false. SourceLoop: do i = 1, phase_Nsources(phase) @@ -720,14 +711,14 @@ function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el enddo SourceLoop -end function constitutive_collectDotState +end function constitutive_collectDotState_source !-------------------------------------------------------------------------------------------------- !> @brief for constitutive models having an instantaneous change of state !> will return false if delta state is not needed/supported by the constitutive model !-------------------------------------------------------------------------------------------------- -function constitutive_deltaState(S, Fe, Fi, ipc, ip, el, phase, of) result(broken) +function constitutive_deltaState_source(Fe, ipc, ip, el, phase, of) result(broken) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -736,48 +727,16 @@ function constitutive_deltaState(S, Fe, Fi, ipc, ip, el, phase, of) result(broke phase, & of real(pReal), intent(in), dimension(3,3) :: & - S, & !< 2nd Piola Kirchhoff stress - Fe, & !< elastic deformation gradient - Fi !< intermediate deformation gradient - real(pReal), dimension(3,3) :: & - Mp + Fe !< elastic deformation gradient integer :: & i, & - instance, & myOffset, & mySize logical :: & broken - Mp = matmul(matmul(transpose(Fi),Fi),S) - instance = phase_plasticityInstance(phase) - - plasticityType: select case (phase_plasticity(phase)) - - case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_deltaState(Mp,instance,of) - broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of))) - - case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_deltaState(Mp,instance,of,ip,el) - broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of))) - - case default - broken = .false. - - end select plasticityType - - if(.not. broken) then - select case(phase_plasticity(phase)) - case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID) - - myOffset = plasticState(phase)%offsetDeltaState - mySize = plasticState(phase)%sizeDeltaState - plasticState(phase)%state(myOffset + 1:myOffset + mySize,of) = & - plasticState(phase)%state(myOffset + 1:myOffset + mySize,of) + plasticState(phase)%deltaState(1:mySize,of) - end select - endif + broken = .false. sourceLoop: do i = 1, phase_Nsources(phase) @@ -786,7 +745,7 @@ function constitutive_deltaState(S, Fe, Fi, ipc, ip, el, phase, of) result(broke case (SOURCE_damage_isoBrittle_ID) sourceType call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, & ipc, ip, el) - broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%deltaState(:,of))) + broken = any(IEEE_is_NaN(sourceState(phase)%p(i)%deltaState(:,of))) if(.not. broken) then myOffset = sourceState(phase)%p(i)%offsetDeltaState mySize = sourceState(phase)%p(i)%sizeDeltaState @@ -798,7 +757,7 @@ function constitutive_deltaState(S, Fe, Fi, ipc, ip, el, phase, of) result(broke enddo SourceLoop -end function constitutive_deltaState +end function constitutive_deltaState_source !-------------------------------------------------------------------------------------------------- @@ -834,6 +793,44 @@ subroutine constitutive_allocateState(state, & end subroutine constitutive_allocateState +!-------------------------------------------------------------------------------------------------- +!> @brief Restore data after homog cutback. +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_restore(i,e) + + integer, intent(in) :: & + i, & !< integration point number + e !< element number + integer :: & + c, & !< constituent number + s + + do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) + do s = 1, phase_Nsources(material_phaseAt(c,e)) + sourceState(material_phaseAt(c,e))%p(s)%state( :,material_phasememberAt(c,i,e)) = & + sourceState(material_phaseAt(c,e))%p(s)%partitionedState0(:,material_phasememberAt(c,i,e)) + enddo + enddo + +end subroutine constitutive_restore + + +!-------------------------------------------------------------------------------------------------- +!> @brief Forward data after successful increment. +! ToDo: Any guessing for the current states possible? +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_forward + + integer :: i, j + + do i = 1, size(sourceState) + do j = 1,phase_Nsources(i) + sourceState(i)%p(j)%state0 = sourceState(i)%p(j)%state + enddo; enddo + +end subroutine constitutive_forward + + !-------------------------------------------------------------------------------------------------- !> @brief writes constitutive results to HDF5 output file !-------------------------------------------------------------------------------------------------- @@ -845,4 +842,972 @@ subroutine constitutive_results end subroutine constitutive_results +!-------------------------------------------------------------------------------------------------- +!> @brief allocates and initialize per grain variables +!-------------------------------------------------------------------------------------------------- +subroutine crystallite_init + + integer :: & + Nconstituents, & + p, & + m, & + c, & !< counter in integration point component loop + i, & !< counter in integration point loop + e, & !< counter in element loop + cMax, & !< maximum number of integration point components + iMax, & !< maximum number of integration points + eMax !< maximum number of elements + + + class(tNode), pointer :: & + num_crystallite, & + debug_crystallite, & ! pointer to debug options for crystallite + phases, & + phase, & + mech + + print'(/,a)', ' <<<+- crystallite init -+>>>' + + debug_crystallite => config_debug%get('crystallite', defaultVal=emptyList) + debugCrystallite%extensive = debug_crystallite%contains('extensive') + + cMax = homogenization_maxNconstituents + iMax = discretization_nIPs + eMax = discretization_Nelems + + allocate(crystallite_partitionedF(3,3,cMax,iMax,eMax),source=0.0_pReal) + + allocate(crystallite_S0, & + crystallite_F0,crystallite_Fp0,crystallite_Lp0, & + crystallite_partitionedS0, & + crystallite_partitionedF0,crystallite_partitionedFp0,& + crystallite_partitionedLp0, & + crystallite_S,crystallite_P, & + crystallite_Fe,crystallite_Fp,crystallite_Lp, & + crystallite_subF,crystallite_subF0, & + crystallite_subFp0,crystallite_subFi0, & + source = crystallite_partitionedF) + + allocate(crystallite_dt(cMax,iMax,eMax),source=0.0_pReal) + allocate(crystallite_subdt,crystallite_subStep, & + source = crystallite_dt) + + allocate(crystallite_orientation(cMax,iMax,eMax)) + + allocate(crystallite_requested(cMax,iMax,eMax), source=.false.) + allocate(crystallite_converged(cMax,iMax,eMax), source=.true.) + + num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) + + num%subStepMinCryst = num_crystallite%get_asFloat ('subStepMin', defaultVal=1.0e-3_pReal) + num%subStepSizeCryst = num_crystallite%get_asFloat ('subStepSize', defaultVal=0.25_pReal) + num%stepIncreaseCryst = num_crystallite%get_asFloat ('stepIncrease', defaultVal=1.5_pReal) + num%subStepSizeLp = num_crystallite%get_asFloat ('subStepSizeLp', defaultVal=0.5_pReal) + num%subStepSizeLi = num_crystallite%get_asFloat ('subStepSizeLi', defaultVal=0.5_pReal) + num%rtol_crystalliteState = num_crystallite%get_asFloat ('rtol_State', defaultVal=1.0e-6_pReal) + num%rtol_crystalliteStress = num_crystallite%get_asFloat ('rtol_Stress', defaultVal=1.0e-6_pReal) + num%atol_crystalliteStress = num_crystallite%get_asFloat ('atol_Stress', defaultVal=1.0e-8_pReal) + num%iJacoLpresiduum = num_crystallite%get_asInt ('iJacoLpresiduum', defaultVal=1) + num%nState = num_crystallite%get_asInt ('nState', defaultVal=20) + num%nStress = num_crystallite%get_asInt ('nStress', defaultVal=40) + + if(num%subStepMinCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepMinCryst') + if(num%subStepSizeCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeCryst') + if(num%stepIncreaseCryst <= 0.0_pReal) call IO_error(301,ext_msg='stepIncreaseCryst') + + if(num%subStepSizeLp <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeLp') + if(num%subStepSizeLi <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeLi') + + if(num%rtol_crystalliteState <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteState') + if(num%rtol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteStress') + if(num%atol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='atol_crystalliteStress') + + if(num%iJacoLpresiduum < 1) call IO_error(301,ext_msg='iJacoLpresiduum') + + if(num%nState < 1) call IO_error(301,ext_msg='nState') + if(num%nStress< 1) call IO_error(301,ext_msg='nStress') + + select case(num_crystallite%get_asString('integrator',defaultVal='FPI')) + case('FPI') + integrateState => integrateStateFPI + case('Euler') + integrateState => integrateStateEuler + case('AdaptiveEuler') + integrateState => integrateStateAdaptiveEuler + case('RK4') + integrateState => integrateStateRK4 + case('RKCK45') + integrateState => integrateStateRKCK45 + case default + call IO_error(301,ext_msg='integrator') + end select + + phases => config_material%get('phase') + + allocate(output_constituent(phases%length)) + allocate(constitutive_mech_Fi(phases%length)) + allocate(constitutive_mech_Fi0(phases%length)) + allocate(constitutive_mech_partionedFi0(phases%length)) + allocate(constitutive_mech_Li(phases%length)) + allocate(constitutive_mech_Li0(phases%length)) + allocate(constitutive_mech_partionedLi0(phases%length)) + do p = 1, phases%length + Nconstituents = count(material_phaseAt == p) * discretization_nIPs + phase => phases%get(p) + mech => phase%get('mechanics',defaultVal = emptyDict) +#if defined(__GFORTRAN__) + output_constituent(p)%label = output_asStrings(mech) +#else + output_constituent(p)%label = mech%get_asStrings('output',defaultVal=emptyStringArray) +#endif + allocate(constitutive_mech_Fi(p)%data(3,3,Nconstituents)) + allocate(constitutive_mech_Fi0(p)%data(3,3,Nconstituents)) + allocate(constitutive_mech_partionedFi0(p)%data(3,3,Nconstituents)) + allocate(constitutive_mech_Li(p)%data(3,3,Nconstituents)) + allocate(constitutive_mech_Li0(p)%data(3,3,Nconstituents)) + allocate(constitutive_mech_partionedLi0(p)%data(3,3,Nconstituents)) + enddo + + print'(a42,1x,i10)', ' # of elements: ', eMax + print'(a42,1x,i10)', ' # of integration points/element: ', iMax + print'(a42,1x,i10)', 'max # of constituents/integration point: ', cMax + flush(IO_STDOUT) + + !$OMP PARALLEL DO PRIVATE(p,m) + do e = FEsolving_execElem(1),FEsolving_execElem(2) + do i = FEsolving_execIP(1), FEsolving_execIP(2); do c = 1, homogenization_Nconstituents(material_homogenizationAt(e)) + + p = material_phaseAt(c,e) + m = material_phaseMemberAt(c,i,e) + crystallite_Fp0(1:3,1:3,c,i,e) = material_orientation0(c,i,e)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) + crystallite_Fp0(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) & + / math_det33(crystallite_Fp0(1:3,1:3,c,i,e))**(1.0_pReal/3.0_pReal) + constitutive_mech_Fi0(p)%data(1:3,1:3,m) = math_I3 + + crystallite_F0(1:3,1:3,c,i,e) = math_I3 + + crystallite_Fe(1:3,1:3,c,i,e) = math_inv33(matmul(constitutive_mech_Fi0(p)%data(1:3,1:3,m), & + crystallite_Fp0(1:3,1:3,c,i,e))) ! assuming that euler angles are given in internal strain free configuration + crystallite_Fp(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) + constitutive_mech_Fi(p)%data(1:3,1:3,m) = constitutive_mech_Fi0(p)%data(1:3,1:3,m) + + constitutive_mech_partionedFi0(p)%data(1:3,1:3,m) = constitutive_mech_Fi0(p)%data(1:3,1:3,m) + + + crystallite_requested(c,i,e) = .true. + enddo; enddo + enddo + !$OMP END PARALLEL DO + + + crystallite_partitionedFp0 = crystallite_Fp0 + crystallite_partitionedF0 = crystallite_F0 + crystallite_partitionedF = crystallite_F0 + + call crystallite_orientations() + + !$OMP PARALLEL DO + do e = FEsolving_execElem(1),FEsolving_execElem(2) + do i = FEsolving_execIP(1),FEsolving_execIP(2) + do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) + call constitutive_plastic_dependentState(crystallite_partitionedF0(1:3,1:3,c,i,e), & + crystallite_partitionedFp0(1:3,1:3,c,i,e), & + c,i,e) ! update dependent state variables to be consistent with basic states + enddo + enddo + enddo + !$OMP END PARALLEL DO + + +end subroutine crystallite_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculate stress (P) +!-------------------------------------------------------------------------------------------------- +function crystallite_stress() + + logical, dimension(discretization_nIPs,discretization_Nelems) :: crystallite_stress + real(pReal) :: & + formerSubStep + integer :: & + NiterationCrystallite, & ! number of iterations in crystallite loop + c, & !< counter in integration point component loop + i, & !< counter in integration point loop + e, & !< counter in element loop + s, p, m + logical, dimension(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems) :: todo !ToDo: need to set some values to false for different Ngrains + real(pReal), dimension(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems) :: subFrac !ToDo: need to set some values to false for different Ngrains + real(pReal), dimension(:,:,:,:,:), allocatable :: & + subLp0,& !< plastic velocity grad at start of crystallite inc + subLi0 !< intermediate velocity grad at start of crystallite inc + + todo = .false. + + allocate(subLi0(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems)) + subLp0 = crystallite_partitionedLp0 + +!-------------------------------------------------------------------------------------------------- +! initialize to starting condition + crystallite_subStep = 0.0_pReal + !$OMP PARALLEL DO PRIVATE(p,m) + elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2) + do i = FEsolving_execIP(1),FEsolving_execIP(2); do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) + p = material_phaseAt(c,e) + m = material_phaseMemberAt(c,i,e) + subLi0(1:3,1:3,c,i,e) = constitutive_mech_partionedLi0(p)%data(1:3,1:3,m) + homogenizationRequestsCalculation: if (crystallite_requested(c,i,e)) then + plasticState (material_phaseAt(c,e))%subState0( :,material_phaseMemberAt(c,i,e)) = & + plasticState (material_phaseAt(c,e))%partitionedState0(:,material_phaseMemberAt(c,i,e)) + + do s = 1, phase_Nsources(material_phaseAt(c,e)) + sourceState(material_phaseAt(c,e))%p(s)%subState0( :,material_phaseMemberAt(c,i,e)) = & + sourceState(material_phaseAt(c,e))%p(s)%partitionedState0(:,material_phaseMemberAt(c,i,e)) + enddo + crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_partitionedFp0(1:3,1:3,c,i,e) + crystallite_subFi0(1:3,1:3,c,i,e) = constitutive_mech_partionedFi0(p)%data(1:3,1:3,m) + crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partitionedF0(1:3,1:3,c,i,e) + subFrac(c,i,e) = 0.0_pReal + crystallite_subStep(c,i,e) = 1.0_pReal/num%subStepSizeCryst + todo(c,i,e) = .true. + crystallite_converged(c,i,e) = .false. ! pretend failed step of 1/subStepSizeCryst + endif homogenizationRequestsCalculation + enddo; enddo + enddo elementLooping1 + !$OMP END PARALLEL DO + + NiterationCrystallite = 0 + cutbackLooping: do while (any(todo(:,FEsolving_execIP(1):FEsolving_execIP(2),FEsolving_execELem(1):FEsolving_execElem(2)))) + NiterationCrystallite = NiterationCrystallite + 1 + +#ifdef DEBUG + if (debugCrystallite%extensive) & + print'(a,i6)', '<< CRYST stress >> crystallite iteration ',NiterationCrystallite +#endif + !$OMP PARALLEL DO PRIVATE(formerSubStep,p,m) + elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2) + do i = FEsolving_execIP(1),FEsolving_execIP(2) + do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) + p = material_phaseAt(c,e) + m = material_phaseMemberAt(c,i,e) +!-------------------------------------------------------------------------------------------------- +! wind forward + if (crystallite_converged(c,i,e)) then + formerSubStep = crystallite_subStep(c,i,e) + subFrac(c,i,e) = subFrac(c,i,e) + crystallite_subStep(c,i,e) + crystallite_subStep(c,i,e) = min(1.0_pReal - subFrac(c,i,e), & + num%stepIncreaseCryst * crystallite_subStep(c,i,e)) + + todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on? + + if (todo(c,i,e)) then + crystallite_subF0 (1:3,1:3,c,i,e) = crystallite_subF(1:3,1:3,c,i,e) + subLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e) + subLi0(1:3,1:3,c,i,e) = constitutive_mech_Li(p)%data(1:3,1:3,m) + crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_Fp (1:3,1:3,c,i,e) + crystallite_subFi0(1:3,1:3,c,i,e) = constitutive_mech_Fi(p)%data(1:3,1:3,m) + plasticState( material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e)) & + = plasticState(material_phaseAt(c,e))%state( :,material_phaseMemberAt(c,i,e)) + do s = 1, phase_Nsources(material_phaseAt(c,e)) + sourceState( material_phaseAt(c,e))%p(s)%subState0(:,material_phaseMemberAt(c,i,e)) & + = sourceState(material_phaseAt(c,e))%p(s)%state( :,material_phaseMemberAt(c,i,e)) + enddo + endif + +!-------------------------------------------------------------------------------------------------- +! cut back (reduced time and restore) + else + crystallite_subStep(c,i,e) = num%subStepSizeCryst * crystallite_subStep(c,i,e) + crystallite_Fp (1:3,1:3,c,i,e) = crystallite_subFp0(1:3,1:3,c,i,e) + constitutive_mech_Fi(p)%data(1:3,1:3,m) = crystallite_subFi0(1:3,1:3,c,i,e) + crystallite_S (1:3,1:3,c,i,e) = crystallite_S0 (1:3,1:3,c,i,e) + if (crystallite_subStep(c,i,e) < 1.0_pReal) then ! actual (not initial) cutback + crystallite_Lp (1:3,1:3,c,i,e) = subLp0(1:3,1:3,c,i,e) + constitutive_mech_Li(p)%data(1:3,1:3,m) = subLi0(1:3,1:3,c,i,e) + endif + plasticState (material_phaseAt(c,e))%state( :,material_phaseMemberAt(c,i,e)) & + = plasticState(material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e)) + do s = 1, phase_Nsources(material_phaseAt(c,e)) + sourceState( material_phaseAt(c,e))%p(s)%state( :,material_phaseMemberAt(c,i,e)) & + = sourceState(material_phaseAt(c,e))%p(s)%subState0(:,material_phaseMemberAt(c,i,e)) + enddo + + ! cant restore dotState here, since not yet calculated in first cutback after initialization + todo(c,i,e) = crystallite_subStep(c,i,e) > num%subStepMinCryst ! still on track or already done (beyond repair) + endif + +!-------------------------------------------------------------------------------------------------- +! prepare for integration + if (todo(c,i,e)) then + crystallite_subF(1:3,1:3,c,i,e) = crystallite_subF0(1:3,1:3,c,i,e) & + + crystallite_subStep(c,i,e) *( crystallite_partitionedF (1:3,1:3,c,i,e) & + -crystallite_partitionedF0(1:3,1:3,c,i,e)) + crystallite_Fe(1:3,1:3,c,i,e) = matmul(crystallite_subF(1:3,1:3,c,i,e), & + math_inv33(matmul(constitutive_mech_Fi(p)%data(1:3,1:3,m), & + crystallite_Fp(1:3,1:3,c,i,e)))) + crystallite_subdt(c,i,e) = crystallite_subStep(c,i,e) * crystallite_dt(c,i,e) + crystallite_converged(c,i,e) = .false. + call integrateState(c,i,e) + call integrateSourceState(c,i,e) + endif + + enddo + enddo + enddo elementLooping3 + !$OMP END PARALLEL DO + +!-------------------------------------------------------------------------------------------------- +! integrate --- requires fully defined state array (basic + dependent state) + where(.not. crystallite_converged .and. crystallite_subStep > num%subStepMinCryst) & ! do not try non-converged but fully cutbacked any further + todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation + enddo cutbackLooping + +! return whether converged or not + crystallite_stress = .false. + elementLooping5: do e = FEsolving_execElem(1),FEsolving_execElem(2) + do i = FEsolving_execIP(1),FEsolving_execIP(2) + crystallite_stress(i,e) = all(crystallite_converged(:,i,e)) + enddo + enddo elementLooping5 + +end function crystallite_stress + + +!-------------------------------------------------------------------------------------------------- +!> @brief Backup data for homog cutback. +!-------------------------------------------------------------------------------------------------- +subroutine crystallite_initializeRestorationPoints(i,e) + + integer, intent(in) :: & + i, & !< integration point number + e !< element number + integer :: & + c, & !< constituent number + s,p, m + + do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) + p = material_phaseAt(c,e) + m = material_phaseMemberAt(c,i,e) + crystallite_partitionedFp0(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) + crystallite_partitionedLp0(1:3,1:3,c,i,e) = crystallite_Lp0(1:3,1:3,c,i,e) + constitutive_mech_partionedFi0(p)%data(1:3,1:3,m) = constitutive_mech_Fi0(p)%data(1:3,1:3,m) + constitutive_mech_partionedLi0(p)%data(1:3,1:3,m) = constitutive_mech_Li0(p)%data(1:3,1:3,m) + crystallite_partitionedF0(1:3,1:3,c,i,e) = crystallite_F0(1:3,1:3,c,i,e) + crystallite_partitionedS0(1:3,1:3,c,i,e) = crystallite_S0(1:3,1:3,c,i,e) + + plasticState(material_phaseAt(c,e))%partitionedState0(:,material_phasememberAt(c,i,e)) = & + plasticState(material_phaseAt(c,e))%state0( :,material_phasememberAt(c,i,e)) + do s = 1, phase_Nsources(material_phaseAt(c,e)) + sourceState(material_phaseAt(c,e))%p(s)%partitionedState0(:,material_phasememberAt(c,i,e)) = & + sourceState(material_phaseAt(c,e))%p(s)%state0( :,material_phasememberAt(c,i,e)) + enddo + enddo + +end subroutine crystallite_initializeRestorationPoints + + +!-------------------------------------------------------------------------------------------------- +!> @brief Wind homog inc forward. +!-------------------------------------------------------------------------------------------------- +subroutine crystallite_windForward(i,e) + + integer, intent(in) :: & + i, & !< integration point number + e !< element number + integer :: & + c, & !< constituent number + s, p, m + do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) + p = material_phaseAt(c,e) + m = material_phaseMemberAt(c,i,e) + crystallite_partitionedF0 (1:3,1:3,c,i,e) = crystallite_partitionedF(1:3,1:3,c,i,e) + crystallite_partitionedFp0(1:3,1:3,c,i,e) = crystallite_Fp (1:3,1:3,c,i,e) + crystallite_partitionedLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e) + constitutive_mech_partionedFi0(p)%data(1:3,1:3,m) = constitutive_mech_Fi(p)%data(1:3,1:3,m) + constitutive_mech_partionedLi0(p)%data(1:3,1:3,m) = constitutive_mech_Li(p)%data(1:3,1:3,m) + crystallite_partitionedS0 (1:3,1:3,c,i,e) = crystallite_S (1:3,1:3,c,i,e) + + plasticState (material_phaseAt(c,e))%partitionedState0(:,material_phasememberAt(c,i,e)) = & + plasticState (material_phaseAt(c,e))%state (:,material_phasememberAt(c,i,e)) + do s = 1, phase_Nsources(material_phaseAt(c,e)) + sourceState(material_phaseAt(c,e))%p(s)%partitionedState0(:,material_phasememberAt(c,i,e)) = & + sourceState(material_phaseAt(c,e))%p(s)%state (:,material_phasememberAt(c,i,e)) + enddo + enddo + +end subroutine crystallite_windForward + + +!-------------------------------------------------------------------------------------------------- +!> @brief Restore data after homog cutback. +!-------------------------------------------------------------------------------------------------- +subroutine crystallite_restore(i,e,includeL) + + integer, intent(in) :: & + i, & !< integration point number + e !< element number + logical, intent(in) :: & + includeL !< protect agains fake cutback + integer :: & + c, p, m !< constituent number + + do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) + p = material_phaseAt(c,e) + m = material_phaseMemberAt(c,i,e) + if (includeL) then + crystallite_Lp(1:3,1:3,c,i,e) = crystallite_partitionedLp0(1:3,1:3,c,i,e) + constitutive_mech_Li(p)%data(1:3,1:3,m) = constitutive_mech_partionedLi0(p)%data(1:3,1:3,m) + endif ! maybe protecting everything from overwriting makes more sense + + crystallite_Fp(1:3,1:3,c,i,e) = crystallite_partitionedFp0(1:3,1:3,c,i,e) + constitutive_mech_Fi(p)%data(1:3,1:3,m) = constitutive_mech_partionedFi0(p)%data(1:3,1:3,m) + crystallite_S (1:3,1:3,c,i,e) = crystallite_partitionedS0 (1:3,1:3,c,i,e) + + plasticState (material_phaseAt(c,e))%state( :,material_phasememberAt(c,i,e)) = & + plasticState (material_phaseAt(c,e))%partitionedState0(:,material_phasememberAt(c,i,e)) + enddo + +end subroutine crystallite_restore + + +!-------------------------------------------------------------------------------------------------- +!> @brief Calculate tangent (dPdF). +!-------------------------------------------------------------------------------------------------- +function crystallite_stressTangent(c,i,e) result(dPdF) + + real(pReal), dimension(3,3,3,3) :: dPdF + integer, intent(in) :: & + c, & !< counter in constituent loop + i, & !< counter in integration point loop + e !< counter in element loop + integer :: & + o, & + p, pp, m + + real(pReal), dimension(3,3) :: devNull, & + invSubFp0,invSubFi0,invFp,invFi, & + temp_33_1, temp_33_2, temp_33_3, temp_33_4 + real(pReal), dimension(3,3,3,3) :: dSdFe, & + dSdF, & + dSdFi, & + dLidS, & ! tangent in lattice configuration + dLidFi, & + dLpdS, & + dLpdFi, & + dFidS, & + dFpinvdF, & + rhs_3333, & + lhs_3333, & + temp_3333 + real(pReal), dimension(9,9):: temp_99 + logical :: error + + pp = material_phaseAt(c,e) + m = material_phaseMemberAt(c,i,e) + + call constitutive_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & + crystallite_Fe(1:3,1:3,c,i,e), & + constitutive_mech_Fi(pp)%data(1:3,1:3,m),c,i,e) + call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, & + crystallite_S (1:3,1:3,c,i,e), & + constitutive_mech_Fi(pp)%data(1:3,1:3,m), & + c,i,e) + + invFp = math_inv33(crystallite_Fp(1:3,1:3,c,i,e)) + invFi = math_inv33(constitutive_mech_Fi(pp)%data(1:3,1:3,m)) + invSubFp0 = math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)) + invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,c,i,e)) + + if (sum(abs(dLidS)) < tol_math_check) then + dFidS = 0.0_pReal + else + lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal + do o=1,3; do p=1,3 + lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) & + + crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) + lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) & + + invFi*invFi(p,o) + rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) & + - crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidS(1:3,1:3,o,p)) + enddo; enddo + call math_invert(temp_99,error,math_3333to99(lhs_3333)) + if (error) then + call IO_warning(warning_ID=600,el=e,ip=i,g=c, & + ext_msg='inversion error in analytic tangent calculation') + dFidS = 0.0_pReal + else + dFidS = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333) + endif + dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS + endif + + call constitutive_plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & + crystallite_S (1:3,1:3,c,i,e), & + constitutive_mech_Fi(pp)%data(1:3,1:3,m),c,i,e) + dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS + +!-------------------------------------------------------------------------------------------------- +! calculate dSdF + temp_33_1 = transpose(matmul(invFp,invFi)) + temp_33_2 = matmul(crystallite_subF(1:3,1:3,c,i,e),invSubFp0) + temp_33_3 = matmul(matmul(crystallite_subF(1:3,1:3,c,i,e),invFp), invSubFi0) + + do o=1,3; do p=1,3 + rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1) + temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) & + + matmul(temp_33_3,dLidS(1:3,1:3,p,o)) + enddo; enddo + lhs_3333 = crystallite_subdt(c,i,e)*math_mul3333xx3333(dSdFe,temp_3333) & + + math_mul3333xx3333(dSdFi,dFidS) + + call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333)) + if (error) then + call IO_warning(warning_ID=600,el=e,ip=i,g=c, & + ext_msg='inversion error in analytic tangent calculation') + dSdF = rhs_3333 + else + dSdF = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333) + endif + +!-------------------------------------------------------------------------------------------------- +! calculate dFpinvdF + temp_3333 = math_mul3333xx3333(dLpdS,dSdF) + do o=1,3; do p=1,3 + dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(c,i,e) & + * matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) + enddo; enddo + +!-------------------------------------------------------------------------------------------------- +! assemble dPdF + temp_33_1 = matmul(crystallite_S(1:3,1:3,c,i,e),transpose(invFp)) + temp_33_2 = matmul(invFp,temp_33_1) + temp_33_3 = matmul(crystallite_subF(1:3,1:3,c,i,e),invFp) + temp_33_4 = matmul(temp_33_3,crystallite_S(1:3,1:3,c,i,e)) + + dPdF = 0.0_pReal + do p=1,3 + dPdF(p,1:3,p,1:3) = transpose(temp_33_2) + enddo + do o=1,3; do p=1,3 + dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) & + + matmul(matmul(crystallite_subF(1:3,1:3,c,i,e), & + dFpinvdF(1:3,1:3,p,o)),temp_33_1) & + + matmul(matmul(temp_33_3,dSdF(1:3,1:3,p,o)), & + transpose(invFp)) & + + matmul(temp_33_4,transpose(dFpinvdF(1:3,1:3,p,o))) + enddo; enddo + +end function crystallite_stressTangent + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates orientations +!-------------------------------------------------------------------------------------------------- +subroutine crystallite_orientations + + integer & + c, & !< counter in integration point component loop + i, & !< counter in integration point loop + e !< counter in element loop + + !$OMP PARALLEL DO + do e = FEsolving_execElem(1),FEsolving_execElem(2) + do i = FEsolving_execIP(1),FEsolving_execIP(2) + do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) + call crystallite_orientation(c,i,e)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,c,i,e)))) + enddo; enddo; enddo + !$OMP END PARALLEL DO + + nonlocalPresent: if (any(plasticState%nonlocal)) then + !$OMP PARALLEL DO + do e = FEsolving_execElem(1),FEsolving_execElem(2) + if (plasticState(material_phaseAt(1,e))%nonlocal) then + do i = FEsolving_execIP(1),FEsolving_execIP(2) + call plastic_nonlocal_updateCompatibility(crystallite_orientation, & + phase_plasticityInstance(material_phaseAt(1,e)),i,e) + enddo + endif + enddo + !$OMP END PARALLEL DO + endif nonlocalPresent + +end subroutine crystallite_orientations + + +!-------------------------------------------------------------------------------------------------- +!> @brief Map 2nd order tensor to reference config +!-------------------------------------------------------------------------------------------------- +function crystallite_push33ToRef(ipc,ip,el, tensor33) + + real(pReal), dimension(3,3) :: crystallite_push33ToRef + real(pReal), dimension(3,3), intent(in) :: tensor33 + real(pReal), dimension(3,3) :: T + integer, intent(in):: & + el, & + ip, & + ipc + + T = matmul(material_orientation0(ipc,ip,el)%asMatrix(), & ! ToDo: initial orientation correct? + transpose(math_inv33(crystallite_subF(1:3,1:3,ipc,ip,el)))) + crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) + +end function crystallite_push33ToRef + + +!-------------------------------------------------------------------------------------------------- +!> @brief writes crystallite results to HDF5 output file +!-------------------------------------------------------------------------------------------------- +subroutine crystallite_results + + integer :: p,o + real(pReal), allocatable, dimension(:,:,:) :: selected_tensors + real(pReal), allocatable, dimension(:,:) :: selected_rotations + character(len=:), allocatable :: group,structureLabel + + do p=1,size(material_name_phase) + group = trim('current/phase')//'/'//trim(material_name_phase(p))//'/mechanics' + + call results_closeGroup(results_addGroup(group)) + + do o = 1, size(output_constituent(p)%label) + select case (output_constituent(p)%label(o)) + case('F') + selected_tensors = select_tensors(crystallite_partitionedF,p) + call results_writeDataset(group,selected_tensors,output_constituent(p)%label(o),& + 'deformation gradient','1') + case('F_e') + selected_tensors = select_tensors(crystallite_Fe,p) + call results_writeDataset(group,selected_tensors,output_constituent(p)%label(o),& + 'elastic deformation gradient','1') + case('F_p') + selected_tensors = select_tensors(crystallite_Fp,p) + call results_writeDataset(group,selected_tensors,output_constituent(p)%label(o),& + 'plastic deformation gradient','1') + case('F_i') + call results_writeDataset(group,constitutive_mech_Fi(p)%data,output_constituent(p)%label(o),& + 'inelastic deformation gradient','1') + case('L_p') + selected_tensors = select_tensors(crystallite_Lp,p) + call results_writeDataset(group,selected_tensors,output_constituent(p)%label(o),& + 'plastic velocity gradient','1/s') + case('L_i') + call results_writeDataset(group,constitutive_mech_Li(p)%data,output_constituent(p)%label(o),& + 'inelastic velocity gradient','1/s') + case('P') + selected_tensors = select_tensors(crystallite_P,p) + call results_writeDataset(group,selected_tensors,output_constituent(p)%label(o),& + 'First Piola-Kirchhoff stress','Pa') + case('S') + selected_tensors = select_tensors(crystallite_S,p) + call results_writeDataset(group,selected_tensors,output_constituent(p)%label(o),& + 'Second Piola-Kirchhoff stress','Pa') + case('O') + select case(lattice_structure(p)) + case(lattice_ISO_ID) + structureLabel = 'aP' + case(lattice_FCC_ID) + structureLabel = 'cF' + case(lattice_BCC_ID) + structureLabel = 'cI' + case(lattice_BCT_ID) + structureLabel = 'tI' + case(lattice_HEX_ID) + structureLabel = 'hP' + case(lattice_ORT_ID) + structureLabel = 'oP' + end select + selected_rotations = select_rotations(crystallite_orientation,p) + call results_writeDataset(group,selected_rotations,output_constituent(p)%label(o),& + 'crystal orientation as quaternion','q_0 ') + call results_addAttribute('Lattice',structureLabel,group//'/'//output_constituent(p)%label(o)) + end select + enddo + enddo + + contains + + !------------------------------------------------------------------------------------------------ + !> @brief select tensors for output + !------------------------------------------------------------------------------------------------ + function select_tensors(dataset,instance) + + integer, intent(in) :: instance + real(pReal), dimension(:,:,:,:,:), intent(in) :: dataset + real(pReal), allocatable, dimension(:,:,:) :: select_tensors + integer :: e,i,c,j + + allocate(select_tensors(3,3,count(material_phaseAt==instance)*discretization_nIPs)) + + j=0 + do e = 1, size(material_phaseAt,2) + do i = 1, discretization_nIPs + do c = 1, size(material_phaseAt,1) !ToDo: this needs to be changed for varying Ngrains + if (material_phaseAt(c,e) == instance) then + j = j + 1 + select_tensors(1:3,1:3,j) = dataset(1:3,1:3,c,i,e) + endif + enddo + enddo + enddo + + end function select_tensors + + +!-------------------------------------------------------------------------------------------------- +!> @brief select rotations for output +!-------------------------------------------------------------------------------------------------- + function select_rotations(dataset,instance) + + integer, intent(in) :: instance + type(rotation), dimension(:,:,:), intent(in) :: dataset + real(pReal), allocatable, dimension(:,:) :: select_rotations + integer :: e,i,c,j + + allocate(select_rotations(4,count(material_phaseAt==instance)*homogenization_maxNconstituents*discretization_nIPs)) + + j=0 + do e = 1, size(material_phaseAt,2) + do i = 1, discretization_nIPs + do c = 1, size(material_phaseAt,1) !ToDo: this needs to be changed for varying Ngrains + if (material_phaseAt(c,e) == instance) then + j = j + 1 + select_rotations(1:4,j) = dataset(c,i,e)%asQuaternion() + endif + enddo + enddo + enddo + + end function select_rotations + +end subroutine crystallite_results + + +!-------------------------------------------------------------------------------------------------- +!> @brief integrate stress, state with adaptive 1st order explicit Euler method +!> using Fixed Point Iteration to adapt the stepsize +!-------------------------------------------------------------------------------------------------- +subroutine integrateSourceState(g,i,e) + + integer, intent(in) :: & + e, & !< element index in element loop + i, & !< integration point index in ip loop + g !< grain index in grain loop + integer :: & + NiterationState, & !< number of iterations in state loop + p, & + c, & + s, & + size_pl + integer, dimension(maxval(phase_Nsources)) :: & + size_so + real(pReal) :: & + zeta + real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: & + r ! state residuum + real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState + logical :: & + broken + + p = material_phaseAt(g,e) + c = material_phaseMemberAt(g,i,e) + + broken = constitutive_collectDotState_source(crystallite_S(1:3,1:3,g,i,e), g,i,e,p,c) + if(broken) return + + do s = 1, phase_Nsources(p) + size_so(s) = sourceState(p)%p(s)%sizeDotState + sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%subState0(1:size_so(s),c) & + + sourceState(p)%p(s)%dotState (1:size_so(s),c) & + * crystallite_subdt(g,i,e) + source_dotState(1:size_so(s),2,s) = 0.0_pReal + enddo + + iteration: do NiterationState = 1, num%nState + + do s = 1, phase_Nsources(p) + if(nIterationState > 1) source_dotState(1:size_so(s),2,s) = source_dotState(1:size_so(s),1,s) + source_dotState(1:size_so(s),1,s) = sourceState(p)%p(s)%dotState(:,c) + enddo + + broken = constitutive_collectDotState_source(crystallite_S(1:3,1:3,g,i,e), g,i,e,p,c) + if(broken) exit iteration + + do s = 1, phase_Nsources(p) + zeta = damper(sourceState(p)%p(s)%dotState(:,c), & + source_dotState(1:size_so(s),1,s),& + source_dotState(1:size_so(s),2,s)) + sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta & + + source_dotState(1:size_so(s),1,s)* (1.0_pReal - zeta) + r(1:size_so(s)) = sourceState(p)%p(s)%state (1:size_so(s),c) & + - sourceState(p)%p(s)%subState0(1:size_so(s),c) & + - sourceState(p)%p(s)%dotState (1:size_so(s),c) * crystallite_subdt(g,i,e) + sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%state(1:size_so(s),c) & + - r(1:size_so(s)) + crystallite_converged(g,i,e) = & + crystallite_converged(g,i,e) .and. converged(r(1:size_so(s)), & + sourceState(p)%p(s)%state(1:size_so(s),c), & + sourceState(p)%p(s)%atol(1:size_so(s))) + enddo + + if(crystallite_converged(g,i,e)) then + broken = constitutive_deltaState_source(crystallite_Fe(1:3,1:3,g,i,e),g,i,e,p,c) + exit iteration + endif + + enddo iteration + + + contains + + !-------------------------------------------------------------------------------------------------- + !> @brief calculate the damping for correction of state and dot state + !-------------------------------------------------------------------------------------------------- + real(pReal) pure function damper(current,previous,previous2) + + real(pReal), dimension(:), intent(in) ::& + current, previous, previous2 + + real(pReal) :: dot_prod12, dot_prod22 + + dot_prod12 = dot_product(current - previous, previous - previous2) + dot_prod22 = dot_product(previous - previous2, previous - previous2) + if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then + damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) + else + damper = 1.0_pReal + endif + + end function damper + +end subroutine integrateSourceState + + +!-------------------------------------------------------------------------------------------------- +!> @brief determines whether a point is converged +!-------------------------------------------------------------------------------------------------- +logical pure function converged(residuum,state,atol) + + real(pReal), intent(in), dimension(:) ::& + residuum, state, atol + real(pReal) :: & + rTol + + rTol = num%rTol_crystalliteState + + converged = all(abs(residuum) <= max(atol, rtol*abs(state))) + +end function converged + + +!-------------------------------------------------------------------------------------------------- +!> @brief Write current restart information (Field and constitutive data) to file. +! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively +!-------------------------------------------------------------------------------------------------- +subroutine crystallite_restartWrite + + integer :: i + integer(HID_T) :: fileHandle, groupHandle + character(len=pStringLen) :: fileName, datasetName + + print*, ' writing field and constitutive data required for restart to file';flush(IO_STDOUT) + + write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5' + fileHandle = HDF5_openFile(fileName,'a') + + call HDF5_write(fileHandle,crystallite_partitionedF,'F') + call HDF5_write(fileHandle,crystallite_Fp, 'F_p') + call HDF5_write(fileHandle,crystallite_Lp, 'L_p') + call HDF5_write(fileHandle,crystallite_S, 'S') + + groupHandle = HDF5_addGroup(fileHandle,'phase') + do i = 1,size(material_name_phase) + write(datasetName,'(i0,a)') i,'_omega' + call HDF5_write(groupHandle,plasticState(i)%state,datasetName) + write(datasetName,'(i0,a)') i,'_F_i' + call HDF5_write(groupHandle,constitutive_mech_Fi(i)%data,datasetName) + write(datasetName,'(i0,a)') i,'_L_i' + call HDF5_write(groupHandle,constitutive_mech_Li(i)%data,datasetName) + enddo + call HDF5_closeGroup(groupHandle) + + groupHandle = HDF5_addGroup(fileHandle,'homogenization') + do i = 1, size(material_name_homogenization) + write(datasetName,'(i0,a)') i,'_omega' + call HDF5_write(groupHandle,homogState(i)%state,datasetName) + enddo + call HDF5_closeGroup(groupHandle) + + call HDF5_closeFile(fileHandle) + +end subroutine crystallite_restartWrite + + +!-------------------------------------------------------------------------------------------------- +!> @brief Read data for restart +! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively +!-------------------------------------------------------------------------------------------------- +subroutine crystallite_restartRead + + integer :: i + integer(HID_T) :: fileHandle, groupHandle + character(len=pStringLen) :: fileName, datasetName + + print'(/,a,i0,a)', ' reading restart information of increment from file' + + write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5' + fileHandle = HDF5_openFile(fileName) + + call HDF5_read(fileHandle,crystallite_F0, 'F') + call HDF5_read(fileHandle,crystallite_Fp0,'F_p') + call HDF5_read(fileHandle,crystallite_Lp0,'L_p') + call HDF5_read(fileHandle,crystallite_S0, 'S') + + groupHandle = HDF5_openGroup(fileHandle,'phase') + do i = 1,size(material_name_phase) + write(datasetName,'(i0,a)') i,'_omega' + call HDF5_read(groupHandle,plasticState(i)%state0,datasetName) + write(datasetName,'(i0,a)') i,'_F_i' + call HDF5_read(groupHandle,constitutive_mech_Fi0(i)%data,datasetName) + write(datasetName,'(i0,a)') i,'_L_i' + call HDF5_read(groupHandle,constitutive_mech_Li0(i)%data,datasetName) + enddo + call HDF5_closeGroup(groupHandle) + + groupHandle = HDF5_openGroup(fileHandle,'homogenization') + do i = 1,size(material_name_homogenization) + write(datasetName,'(i0,a)') i,'_omega' + call HDF5_read(groupHandle,homogState(i)%state0,datasetName) + enddo + call HDF5_closeGroup(groupHandle) + + call HDF5_closeFile(fileHandle) + +end subroutine crystallite_restartRead + + +!-------------------------------------------------------------------------------------------------- +!> @brief Forward data after successful increment. +! ToDo: Any guessing for the current states possible? +!-------------------------------------------------------------------------------------------------- +subroutine crystallite_forward + + integer :: i, j + + crystallite_F0 = crystallite_partitionedF + crystallite_Fp0 = crystallite_Fp + crystallite_Lp0 = crystallite_Lp + crystallite_S0 = crystallite_S + + do i = 1, size(plasticState) + plasticState(i)%state0 = plasticState(i)%state + constitutive_mech_Fi0(i) = constitutive_mech_Fi(i) + constitutive_mech_Li0(i) = constitutive_mech_Li(i) + enddo + do i = 1,size(material_name_homogenization) + homogState (i)%state0 = homogState (i)%state + damageState (i)%state0 = damageState (i)%state + enddo + +end subroutine crystallite_forward + end module constitutive diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index 6b3c6fce6..07e48b537 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -133,6 +133,67 @@ submodule(constitutive) constitutive_mech el !< current element number end subroutine plastic_nonlocal_LpAndItsTangent + module subroutine plastic_isotropic_dotState(Mp,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_isotropic_dotState + + module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_phenopowerlaw_dotState + + module subroutine plastic_kinehardening_dotState(Mp,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_kinehardening_dotState + + module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + instance, & + of + end subroutine plastic_dislotwin_dotState + + module subroutine plastic_disloTungsten_dotState(Mp,T,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + instance, & + of + end subroutine plastic_disloTungsten_dotState + + module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & + instance,of,ip,el) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< MandelStress + real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: & + F, & !< deformation gradient + Fp !< plastic deformation gradient + real(pReal), intent(in) :: & + Temperature, & !< temperature + timestep !< substepped crystallite time increment + integer, intent(in) :: & + instance, & + of, & + ip, & !< current integration point + el !< current element number + end subroutine plastic_nonlocal_dotState + module subroutine plastic_dislotwin_dependentState(T,instance,of) integer, intent(in) :: & @@ -159,6 +220,24 @@ submodule(constitutive) constitutive_mech el !< current element number end subroutine plastic_nonlocal_dependentState + module subroutine plastic_kinehardening_deltaState(Mp,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_kinehardening_deltaState + + module subroutine plastic_nonlocal_deltaState(Mp,instance,of,ip,el) + real(pReal), dimension(3,3), intent(in) :: & + Mp + integer, intent(in) :: & + instance, & + of, & + ip, & + el + end subroutine plastic_nonlocal_deltaState + module subroutine plastic_isotropic_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group @@ -454,6 +533,124 @@ module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & end subroutine constitutive_plastic_LpAndItsTangents +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +module function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el, & !< element + phase, & + of + real(pReal), intent(in) :: & + subdt !< timestep + real(pReal), intent(in), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems) :: & + FArray, & !< elastic deformation gradient + FpArray !< plastic deformation gradient + real(pReal), intent(in), dimension(3,3) :: & + Fi !< intermediate deformation gradient + real(pReal), intent(in), dimension(3,3) :: & + S !< 2nd Piola Kirchhoff stress (vector notation) + real(pReal), dimension(3,3) :: & + Mp + integer :: & + ho, & !< homogenization + tme, & !< thermal member position + i, & !< counter in source loop + instance + logical :: broken + + ho = material_homogenizationAt(el) + tme = material_homogenizationMemberAt(ip,el) + instance = phase_plasticityInstance(phase) + + Mp = matmul(matmul(transpose(Fi),Fi),S) + + plasticityType: select case (phase_plasticity(phase)) + + case (PLASTICITY_ISOTROPIC_ID) plasticityType + call plastic_isotropic_dotState(Mp,instance,of) + + case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType + call plastic_phenopowerlaw_dotState(Mp,instance,of) + + case (PLASTICITY_KINEHARDENING_ID) plasticityType + call plastic_kinehardening_dotState(Mp,instance,of) + + case (PLASTICITY_DISLOTWIN_ID) plasticityType + call plastic_dislotwin_dotState(Mp,temperature(ho)%p(tme),instance,of) + + case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType + call plastic_disloTungsten_dotState(Mp,temperature(ho)%p(tme),instance,of) + + case (PLASTICITY_NONLOCAL_ID) plasticityType + call plastic_nonlocal_dotState(Mp,FArray,FpArray,temperature(ho)%p(tme),subdt, & + instance,of,ip,el) + end select plasticityType + broken = any(IEEE_is_NaN(plasticState(phase)%dotState(:,of))) + + +end function constitutive_collectDotState + + +!-------------------------------------------------------------------------------------------------- +!> @brief for constitutive models having an instantaneous change of state +!> will return false if delta state is not needed/supported by the constitutive model +!-------------------------------------------------------------------------------------------------- +module function constitutive_deltaState(S, Fi, ipc, ip, el, phase, of) result(broken) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el, & !< element + phase, & + of + real(pReal), intent(in), dimension(3,3) :: & + S, & !< 2nd Piola Kirchhoff stress + Fi !< intermediate deformation gradient + real(pReal), dimension(3,3) :: & + Mp + integer :: & + instance, & + myOffset, & + mySize + logical :: & + broken + + Mp = matmul(matmul(transpose(Fi),Fi),S) + instance = phase_plasticityInstance(phase) + + plasticityType: select case (phase_plasticity(phase)) + + case (PLASTICITY_KINEHARDENING_ID) plasticityType + call plastic_kinehardening_deltaState(Mp,instance,of) + broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of))) + + case (PLASTICITY_NONLOCAL_ID) plasticityType + call plastic_nonlocal_deltaState(Mp,instance,of,ip,el) + broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of))) + + case default + broken = .false. + + end select plasticityType + + if(.not. broken) then + select case(phase_plasticity(phase)) + case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID) + + myOffset = plasticState(phase)%offsetDeltaState + mySize = plasticState(phase)%sizeDeltaState + plasticState(phase)%state(myOffset + 1:myOffset + mySize,of) = & + plasticState(phase)%state(myOffset + 1:myOffset + mySize,of) + plasticState(phase)%deltaState(1:mySize,of) + end select + endif + +end function constitutive_deltaState + + !-------------------------------------------------------------------------------------------- !> @brief writes plasticity constitutive results to HDF5 output file !-------------------------------------------------------------------------------------------- @@ -494,5 +691,581 @@ module subroutine plastic_results end subroutine plastic_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(ipc,ip,el,timeFraction) result(broken) + + integer, intent(in):: el, & ! element index + ip, & ! integration point index + ipc ! grain index + real(pReal), optional, intent(in) :: timeFraction ! fraction of timestep + + real(pReal), dimension(3,3):: F, & ! deformation gradient at end of timestep + Fp_new, & ! plastic deformation gradient at end of timestep + invFp_new, & ! inverse of Fp_new + invFp_current, & ! inverse of Fp_current + Lpguess, & ! current guess for plastic velocity gradient + Lpguess_old, & ! known last good guess for plastic velocity gradient + Lp_constitutive, & ! plastic velocity gradient resulting from constitutive law + residuumLp, & ! current residuum of plastic velocity gradient + residuumLp_old, & ! last residuum of plastic velocity gradient + deltaLp, & ! direction of next guess + Fi_new, & ! gradient of intermediate deformation stages + invFi_new, & + invFi_current, & ! inverse of Fi_current + Liguess, & ! current guess for intermediate velocity gradient + Liguess_old, & ! known last good guess for intermediate velocity gradient + Li_constitutive, & ! intermediate velocity gradient resulting from constitutive law + residuumLi, & ! current residuum of intermediate velocity gradient + residuumLi_old, & ! last residuum of intermediate velocity gradient + deltaLi, & ! direction of next guess + Fe, & ! elastic deformation gradient + S, & ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration + A, & + B, & + temp_33 + real(pReal), dimension(9) :: temp_9 ! needed for matrix inversion by LAPACK + integer, dimension(9) :: devNull_9 ! needed for matrix inversion by LAPACK + real(pReal), dimension(9,9) :: dRLp_dLp, & ! partial derivative of residuum (Jacobian for Newton-Raphson scheme) + dRLi_dLi ! partial derivative of residuumI (Jacobian for Newton-Raphson scheme) + real(pReal), dimension(3,3,3,3):: dS_dFe, & ! partial derivative of 2nd Piola-Kirchhoff stress + dS_dFi, & + dFe_dLp, & ! partial derivative of elastic deformation gradient + dFe_dLi, & + dFi_dLi, & + dLp_dFi, & + dLi_dFi, & + dLp_dS, & + dLi_dS + real(pReal) steplengthLp, & + steplengthLi, & + dt, & ! time increment + atol_Lp, & + atol_Li, & + devNull + integer NiterationStressLp, & ! number of stress integrations + NiterationStressLi, & ! number of inner stress integrations + ierr, & ! error indicator for LAPACK + o, & + p, & + m, & + jacoCounterLp, & + jacoCounterLi ! counters to check for Jacobian update + logical :: error,broken + + broken = .true. + + if (present(timeFraction)) then + dt = crystallite_subdt(ipc,ip,el) * timeFraction + F = crystallite_subF0(1:3,1:3,ipc,ip,el) & + + (crystallite_subF(1:3,1:3,ipc,ip,el) - crystallite_subF0(1:3,1:3,ipc,ip,el)) * timeFraction + else + dt = crystallite_subdt(ipc,ip,el) + F = crystallite_subF(1:3,1:3,ipc,ip,el) + endif + + call constitutive_plastic_dependentState(crystallite_partitionedF(1:3,1:3,ipc,ip,el), & + crystallite_Fp(1:3,1:3,ipc,ip,el),ipc,ip,el) + + p = material_phaseAt(ipc,el) + m = material_phaseMemberAt(ipc,ip,el) + + Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! take as first guess + Liguess = constitutive_mech_Li(p)%data(1:3,1:3,m) ! take as first guess + + call math_invert33(invFp_current,devNull,error,crystallite_subFp0(1:3,1:3,ipc,ip,el)) + if (error) return ! error + call math_invert33(invFi_current,devNull,error,crystallite_subFi0(1:3,1:3,ipc,ip,el)) + if (error) return ! error + + A = matmul(F,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp + + jacoCounterLi = 0 + steplengthLi = 1.0_pReal + residuumLi_old = 0.0_pReal + Liguess_old = Liguess + + NiterationStressLi = 0 + LiLoop: do + NiterationStressLi = NiterationStressLi + 1 + if (NiterationStressLi>num%nStress) return ! error + + invFi_new = matmul(invFi_current,math_I3 - dt*Liguess) + Fi_new = math_inv33(invFi_new) + + jacoCounterLp = 0 + steplengthLp = 1.0_pReal + residuumLp_old = 0.0_pReal + Lpguess_old = Lpguess + + NiterationStressLp = 0 + LpLoop: do + NiterationStressLp = NiterationStressLp + 1 + if (NiterationStressLp>num%nStress) return ! error + + B = math_I3 - dt*Lpguess + Fe = matmul(matmul(A,B), invFi_new) + call constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & + Fe, Fi_new, ipc, ip, el) + + call constitutive_plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & + S, Fi_new, ipc, ip, el) + + !* update current residuum and check for convergence of loop + atol_Lp = max(num%rtol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error + num%atol_crystalliteStress) ! minimum lower cutoff + residuumLp = Lpguess - Lp_constitutive + + if (any(IEEE_is_NaN(residuumLp))) then + return ! error + elseif (norm2(residuumLp) < atol_Lp) then ! converged if below absolute tolerance + exit LpLoop + elseif (NiterationStressLp == 1 .or. norm2(residuumLp) < norm2(residuumLp_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)... + residuumLp_old = residuumLp ! ...remember old values and... + Lpguess_old = Lpguess + steplengthLp = 1.0_pReal ! ...proceed with normal step length (calculate new search direction) + else ! not converged and residuum not improved... + steplengthLp = num%subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction + Lpguess = Lpguess_old & + + deltaLp * stepLengthLp + cycle LpLoop + endif + + calculateJacobiLi: if (mod(jacoCounterLp, num%iJacoLpresiduum) == 0) then + jacoCounterLp = jacoCounterLp + 1 + + do o=1,3; do p=1,3 + dFe_dLp(o,1:3,p,1:3) = - dt * A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) + enddo; enddo + dRLp_dLp = math_eye(9) & + - math_3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dS,dS_dFe),dFe_dLp)) + temp_9 = math_33to9(residuumLp) + call dgesv(9,1,dRLp_dLp,9,devNull_9,temp_9,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp + if (ierr /= 0) return ! error + deltaLp = - math_9to33(temp_9) + endif calculateJacobiLi + + Lpguess = Lpguess & + + deltaLp * steplengthLp + enddo LpLoop + + call constitutive_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, & + S, Fi_new, ipc, ip, el) + + !* update current residuum and check for convergence of loop + atol_Li = max(num%rtol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error + num%atol_crystalliteStress) ! minimum lower cutoff + residuumLi = Liguess - Li_constitutive + if (any(IEEE_is_NaN(residuumLi))) then + return ! error + elseif (norm2(residuumLi) < atol_Li) then ! converged if below absolute tolerance + exit LiLoop + elseif (NiterationStressLi == 1 .or. norm2(residuumLi) < norm2(residuumLi_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)... + residuumLi_old = residuumLi ! ...remember old values and... + Liguess_old = Liguess + steplengthLi = 1.0_pReal ! ...proceed with normal step length (calculate new search direction) + else ! not converged and residuum not improved... + steplengthLi = num%subStepSizeLi * steplengthLi ! ...try with smaller step length in same direction + Liguess = Liguess_old & + + deltaLi * steplengthLi + cycle LiLoop + endif + + calculateJacobiLp: if (mod(jacoCounterLi, num%iJacoLpresiduum) == 0) then + jacoCounterLi = jacoCounterLi + 1 + + temp_33 = matmul(matmul(A,B),invFi_current) + do o=1,3; do p=1,3 + dFe_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*temp_33 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) + dFi_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*invFi_current + enddo; enddo + do o=1,3; do p=1,3 + dFi_dLi(1:3,1:3,o,p) = matmul(matmul(Fi_new,dFi_dLi(1:3,1:3,o,p)),Fi_new) + enddo; enddo + dRLi_dLi = math_eye(9) & + - math_3333to99(math_mul3333xx3333(dLi_dS, math_mul3333xx3333(dS_dFe, dFe_dLi) & + + math_mul3333xx3333(dS_dFi, dFi_dLi))) & + - math_3333to99(math_mul3333xx3333(dLi_dFi, dFi_dLi)) + temp_9 = math_33to9(residuumLi) + call dgesv(9,1,dRLi_dLi,9,devNull_9,temp_9,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li + if (ierr /= 0) return ! error + deltaLi = - math_9to33(temp_9) + endif calculateJacobiLp + + Liguess = Liguess & + + deltaLi * steplengthLi + enddo LiLoop + + invFp_new = matmul(invFp_current,B) + call math_invert33(Fp_new,devNull,error,invFp_new) + if (error) return ! error + + p = material_phaseAt(ipc,el) + m = material_phaseMemberAt(ipc,ip,el) + + crystallite_P (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) + crystallite_S (1:3,1:3,ipc,ip,el) = S + crystallite_Lp (1:3,1:3,ipc,ip,el) = Lpguess + constitutive_mech_Li(p)%data(1:3,1:3,m) = Liguess + crystallite_Fp (1:3,1:3,ipc,ip,el) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize + constitutive_mech_Fi(p)%data(1:3,1:3,m) = Fi_new + crystallite_Fe (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),invFi_new) + broken = .false. + +end function integrateStress + + +!-------------------------------------------------------------------------------------------------- +!> @brief integrate stress, state with adaptive 1st order explicit Euler method +!> using Fixed Point Iteration to adapt the stepsize +!-------------------------------------------------------------------------------------------------- +module subroutine integrateStateFPI(g,i,e) + + integer, intent(in) :: & + e, & !< element index in element loop + i, & !< integration point index in ip loop + g !< grain index in grain loop + integer :: & + NiterationState, & !< number of iterations in state loop + p, & + c, & + s, & + size_pl + integer, dimension(maxval(phase_Nsources)) :: & + size_so + real(pReal) :: & + zeta + real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: & + r ! state residuum + real(pReal), dimension(constitutive_plasticity_maxSizeDotState,2) :: & + plastic_dotState + real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState + logical :: & + broken + + p = material_phaseAt(g,e) + c = material_phaseMemberAt(g,i,e) + + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partitionedF0, & + constitutive_mech_Fi(p)%data(1:3,1:3,c), & + crystallite_partitionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) return + + size_pl = plasticState(p)%sizeDotState + plasticState(p)%state(1:size_pl,c) = plasticState(p)%subState0(1:size_pl,c) & + + plasticState(p)%dotState (1:size_pl,c) & + * crystallite_subdt(g,i,e) + plastic_dotState(1:size_pl,2) = 0.0_pReal + + iteration: do NiterationState = 1, num%nState + + if(nIterationState > 1) plastic_dotState(1:size_pl,2) = plastic_dotState(1:size_pl,1) + plastic_dotState(1:size_pl,1) = plasticState(p)%dotState(:,c) + + broken = integrateStress(g,i,e) + if(broken) exit iteration + + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partitionedF0, & + constitutive_mech_Fi(p)%data(1:3,1:3,c), & + crystallite_partitionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) exit iteration + + zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState(1:size_pl,1),& + plastic_dotState(1:size_pl,2)) + plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & + + plastic_dotState(1:size_pl,1) * (1.0_pReal - zeta) + r(1:size_pl) = plasticState(p)%state (1:size_pl,c) & + - plasticState(p)%subState0(1:size_pl,c) & + - plasticState(p)%dotState (1:size_pl,c) * crystallite_subdt(g,i,e) + plasticState(p)%state(1:size_pl,c) = plasticState(p)%state(1:size_pl,c) & + - r(1:size_pl) + crystallite_converged(g,i,e) = converged(r(1:size_pl), & + plasticState(p)%state(1:size_pl,c), & + plasticState(p)%atol(1:size_pl)) + + if(crystallite_converged(g,i,e)) then + broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & + constitutive_mech_Fi(p)%data(1:3,1:3,c),g,i,e,p,c) + exit iteration + endif + + enddo iteration + + + contains + + !-------------------------------------------------------------------------------------------------- + !> @brief calculate the damping for correction of state and dot state + !-------------------------------------------------------------------------------------------------- + real(pReal) pure function damper(current,previous,previous2) + + real(pReal), dimension(:), intent(in) ::& + current, previous, previous2 + + real(pReal) :: dot_prod12, dot_prod22 + + dot_prod12 = dot_product(current - previous, previous - previous2) + dot_prod22 = dot_product(previous - previous2, previous - previous2) + if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then + damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) + else + damper = 1.0_pReal + endif + + end function damper + +end subroutine integrateStateFPI + + +!-------------------------------------------------------------------------------------------------- +!> @brief integrate state with 1st order explicit Euler method +!-------------------------------------------------------------------------------------------------- +module subroutine integrateStateEuler(g,i,e) + + integer, intent(in) :: & + e, & !< element index in element loop + i, & !< integration point index in ip loop + g !< grain index in grain loop + integer :: & + p, & + c, & + sizeDotState + logical :: & + broken + + p = material_phaseAt(g,e) + c = material_phaseMemberAt(g,i,e) + + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partitionedF0, & + constitutive_mech_Fi(p)%data(1:3,1:3,c), & + crystallite_partitionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) return + + sizeDotState = plasticState(p)%sizeDotState + plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + + plasticState(p)%dotState (1:sizeDotState,c) & + * crystallite_subdt(g,i,e) + + broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & + constitutive_mech_Fi(p)%data(1:3,1:3,c),g,i,e,p,c) + if(broken) return + + broken = integrateStress(g,i,e) + crystallite_converged(g,i,e) = .not. broken + +end subroutine integrateStateEuler + + +!-------------------------------------------------------------------------------------------------- +!> @brief integrate stress, state with 1st order Euler method with adaptive step size +!-------------------------------------------------------------------------------------------------- +module subroutine integrateStateAdaptiveEuler(g,i,e) + + integer, intent(in) :: & + e, & !< element index in element loop + i, & !< integration point index in ip loop + g !< grain index in grain loop + integer :: & + p, & + c, & + sizeDotState + logical :: & + broken + + real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic + + + p = material_phaseAt(g,e) + c = material_phaseMemberAt(g,i,e) + + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partitionedF0, & + constitutive_mech_Fi(p)%data(1:3,1:3,c), & + crystallite_partitionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) return + + sizeDotState = plasticState(p)%sizeDotState + + residuum_plastic(1:sizeDotState) = - plasticState(p)%dotstate(1:sizeDotState,c) * 0.5_pReal * crystallite_subdt(g,i,e) + plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + + plasticState(p)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) + + broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & + constitutive_mech_Fi(p)%data(1:3,1:3,c),g,i,e,p,c) + if(broken) return + + broken = integrateStress(g,i,e) + if(broken) return + + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partitionedF0, & + constitutive_mech_Fi(p)%data(1:3,1:3,c), & + crystallite_partitionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) return + + + sizeDotState = plasticState(p)%sizeDotState + crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState) & + + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e), & + plasticState(p)%state(1:sizeDotState,c), & + plasticState(p)%atol(1:sizeDotState)) + +end subroutine integrateStateAdaptiveEuler + + +!--------------------------------------------------------------------------------------------------- +!> @brief Integrate state (including stress integration) with the classic Runge Kutta method +!--------------------------------------------------------------------------------------------------- +module subroutine integrateStateRK4(g,i,e) + + integer, intent(in) :: g,i,e + + real(pReal), dimension(3,3), parameter :: & + A = reshape([& + 0.5_pReal, 0.0_pReal, 0.0_pReal, & + 0.0_pReal, 0.5_pReal, 0.0_pReal, & + 0.0_pReal, 0.0_pReal, 1.0_pReal],& + shape(A)) + real(pReal), dimension(3), parameter :: & + C = [0.5_pReal, 0.5_pReal, 1.0_pReal] + real(pReal), dimension(4), parameter :: & + 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] + + call integrateStateRK(g,i,e,A,B,C) + +end subroutine integrateStateRK4 + + +!--------------------------------------------------------------------------------------------------- +!> @brief Integrate state (including stress integration) with the Cash-Carp method +!--------------------------------------------------------------------------------------------------- +module subroutine integrateStateRKCK45(g,i,e) + + integer, intent(in) :: g,i,e + + real(pReal), dimension(5,5), parameter :: & + A = reshape([& + 1._pReal/5._pReal, .0_pReal, .0_pReal, .0_pReal, .0_pReal, & + 3._pReal/40._pReal, 9._pReal/40._pReal, .0_pReal, .0_pReal, .0_pReal, & + 3_pReal/10._pReal, -9._pReal/10._pReal, 6._pReal/5._pReal, .0_pReal, .0_pReal, & + -11._pReal/54._pReal, 5._pReal/2._pReal, -70.0_pReal/27.0_pReal, 35.0_pReal/27.0_pReal, .0_pReal, & + 1631._pReal/55296._pReal,175._pReal/512._pReal,575._pReal/13824._pReal,44275._pReal/110592._pReal,253._pReal/4096._pReal],& + shape(A)) + real(pReal), dimension(5), parameter :: & + C = [0.2_pReal, 0.3_pReal, 0.6_pReal, 1.0_pReal, 0.875_pReal] + real(pReal), dimension(6), parameter :: & + B = & + [37.0_pReal/378.0_pReal, .0_pReal, 250.0_pReal/621.0_pReal, & + 125.0_pReal/594.0_pReal, .0_pReal, 512.0_pReal/1771.0_pReal], & + DB = B - & + [2825.0_pReal/27648.0_pReal, .0_pReal, 18575.0_pReal/48384.0_pReal,& + 13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 1._pReal/4._pReal] + + call integrateStateRK(g,i,e,A,B,C,DB) + +end subroutine integrateStateRKCK45 + + +!-------------------------------------------------------------------------------------------------- +!> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an +!! embedded explicit Runge-Kutta method +!-------------------------------------------------------------------------------------------------- +subroutine integrateStateRK(g,i,e,A,B,CC,DB) + + + real(pReal), dimension(:,:), intent(in) :: A + real(pReal), dimension(:), intent(in) :: B, CC + real(pReal), dimension(:), intent(in), optional :: DB + + integer, intent(in) :: & + e, & !< element index in element loop + i, & !< integration point index in ip loop + g !< grain index in grain loop + integer :: & + stage, & ! stage index in integration stage loop + n, & + p, & + c, & + sizeDotState + logical :: & + broken + real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState + + p = material_phaseAt(g,e) + c = material_phaseMemberAt(g,i,e) + + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partitionedF0, & + constitutive_mech_Fi(p)%data(1:3,1:3,c), & + crystallite_partitionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) return + + do stage = 1,size(A,1) + sizeDotState = plasticState(p)%sizeDotState + plastic_RKdotState(1:sizeDotState,stage) = plasticState(p)%dotState(:,c) + plasticState(p)%dotState(:,c) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1) + + do n = 2, stage + sizeDotState = plasticState(p)%sizeDotState + plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) & + + A(n,stage) * plastic_RKdotState(1:sizeDotState,n) + enddo + + sizeDotState = plasticState(p)%sizeDotState + plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + + plasticState(p)%dotState (1:sizeDotState,c) & + * crystallite_subdt(g,i,e) + + broken = integrateStress(g,i,e,CC(stage)) + if(broken) exit + + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partitionedF0, & + constitutive_mech_Fi(p)%data(1:3,1:3,c), & + crystallite_partitionedFp0, & + crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c) + if(broken) exit + + enddo + if(broken) return + + sizeDotState = plasticState(p)%sizeDotState + + plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (p)%dotState(:,c) + plasticState(p)%dotState(:,c) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B) + plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + + plasticState(p)%dotState (1:sizeDotState,c) & + * crystallite_subdt(g,i,e) + if(present(DB)) & + broken = .not. converged( matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) & + * crystallite_subdt(g,i,e), & + plasticState(p)%state(1:sizeDotState,c), & + plasticState(p)%atol(1:sizeDotState)) + + if(broken) return + + broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & + constitutive_mech_Fi(p)%data(1:3,1:3,c),g,i,e,p,c) + if(broken) return + + broken = integrateStress(g,i,e) + crystallite_converged(g,i,e) = .not. broken + + +end subroutine integrateStateRK + + end submodule constitutive_mech diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index ce9e4e391..65f74b6cc 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -10,7 +10,8 @@ submodule(constitutive:constitutive_mech) plastic_nonlocal IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, & IPvolume => geometry_plastic_nonlocal_IPvolume0, & IParea => geometry_plastic_nonlocal_IParea0, & - IPareaNormal => geometry_plastic_nonlocal_IPareaNormal0 + IPareaNormal => geometry_plastic_nonlocal_IPareaNormal0, & + geometry_plastic_nonlocal_disable real(pReal), parameter :: & kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin @@ -195,7 +196,7 @@ module function plastic_nonlocal_init() result(myPlasticity) call geometry_plastic_nonlocal_disable return endif - + print*, 'Reuber et al., Acta Materialia 71:333–348, 2014' print*, 'https://doi.org/10.1016/j.actamat.2014.03.012'//IO_EOL diff --git a/src/crystallite.f90 b/src/crystallite.f90 deleted file mode 100644 index 31a6bde2d..000000000 --- a/src/crystallite.f90 +++ /dev/null @@ -1,1638 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH -!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH -!> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH -!> @author Chen Zhang, Michigan State University -!> @brief crystallite state integration functions and reporting of results -!-------------------------------------------------------------------------------------------------- - -module crystallite - use prec - use parallelization - use IO - use HDF5_utilities - use DAMASK_interface - use config - use rotations - use math - use FEsolving - use material - use constitutive - use discretization - use lattice - use results - - implicit none - private - - real(pReal), dimension(:,:,:), allocatable, public :: & - crystallite_dt !< requested time increment of each grain - real(pReal), dimension(:,:,:), allocatable :: & - crystallite_subdt, & !< substepped time increment of each grain - crystallite_subFrac, & !< already calculated fraction of increment - crystallite_subStep !< size of next integration step - type(rotation), dimension(:,:,:), allocatable :: & - crystallite_orientation !< current orientation - real(pReal), dimension(:,:,:,:,:), allocatable :: & - crystallite_F0, & !< def grad at start of FE inc - crystallite_subF, & !< def grad to be reached at end of crystallite inc - crystallite_subF0, & !< def grad at start of crystallite inc - ! - crystallite_Fe, & !< current "elastic" def grad (end of converged time step) - ! - crystallite_Fp, & !< current plastic def grad (end of converged time step) - crystallite_Fp0, & !< plastic def grad at start of FE inc - crystallite_partitionedFp0,& !< plastic def grad at start of homog inc - crystallite_subFp0,& !< plastic def grad at start of crystallite inc - ! - crystallite_Fi, & !< current intermediate def grad (end of converged time step) - crystallite_Fi0, & !< intermediate def grad at start of FE inc - crystallite_partitionedFi0,& !< intermediate def grad at start of homog inc - crystallite_subFi0,& !< intermediate def grad at start of crystallite inc - ! - crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc - crystallite_partitionedLp0, & !< plastic velocity grad at start of homog inc - ! - crystallite_Li, & !< current intermediate velocitiy grad (end of converged time step) - crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc - crystallite_partitionedLi0, & !< intermediate velocity grad at start of homog inc - ! - crystallite_S0, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc - crystallite_partitionedS0 !< 2nd Piola-Kirchhoff stress vector at start of homog inc - real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & - crystallite_P, & !< 1st Piola-Kirchhoff stress per grain - crystallite_Lp, & !< current plastic velocitiy grad (end of converged time step) - crystallite_S, & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step) - crystallite_partitionedF0 !< def grad at start of homog inc - real(pReal), dimension(:,:,:,:,:), allocatable, public :: & - crystallite_partitionedF !< def grad to be reached at end of homog inc - - logical, dimension(:,:,:), allocatable, public :: & - crystallite_requested !< used by upper level (homogenization) to request crystallite calculation - logical, dimension(:,:,:), allocatable :: & - crystallite_converged !< convergence flag - - type :: tOutput !< new requested output (per phase) - character(len=pStringLen), allocatable, dimension(:) :: & - label - end type tOutput - type(tOutput), allocatable, dimension(:) :: output_constituent - - type :: tNumerics - integer :: & - iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp - nState, & !< state loop limit - nStress !< stress loop limit - real(pReal) :: & - subStepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback - subStepSizeCryst, & !< size of first substep when cutback - subStepSizeLp, & !< size of first substep when cutback in Lp calculation - subStepSizeLi, & !< size of first substep when cutback in Li calculation - stepIncreaseCryst, & !< increase of next substep size when previous substep converged - rtol_crystalliteState, & !< relative tolerance in state loop - rtol_crystalliteStress, & !< relative tolerance in stress loop - atol_crystalliteStress !< absolute tolerance in stress loop - end type tNumerics - - type(tNumerics) :: num ! numerics parameters. Better name? - - type :: tDebugOptions - logical :: & - basic, & - extensive, & - selective - integer :: & - element, & - ip, & - grain - end type tDebugOptions - - type(tDebugOptions) :: debugCrystallite - - procedure(integrateStateFPI), pointer :: integrateState - - public :: & - crystallite_init, & - crystallite_stress, & - crystallite_stressTangent, & - crystallite_orientations, & - crystallite_push33ToRef, & - crystallite_results, & - crystallite_restartWrite, & - crystallite_restartRead, & - crystallite_forward, & - crystallite_initializeRestorationPoints, & - crystallite_windForward, & - crystallite_restore - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief allocates and initialize per grain variables -!-------------------------------------------------------------------------------------------------- -subroutine crystallite_init - - integer :: & - p, & - c, & !< counter in integration point component loop - i, & !< counter in integration point loop - e, & !< counter in element loop - cMax, & !< maximum number of integration point components - iMax, & !< maximum number of integration points - eMax !< maximum number of elements - - - class(tNode), pointer :: & - num_crystallite, & - debug_crystallite, & ! pointer to debug options for crystallite - phases, & - phase, & - mech - - print'(/,a)', ' <<<+- crystallite init -+>>>' - - debug_crystallite => config_debug%get('crystallite', defaultVal=emptyList) - debugCrystallite%basic = debug_crystallite%contains('basic') - debugCrystallite%extensive = debug_crystallite%contains('extensive') - debugCrystallite%selective = debug_crystallite%contains('selective') - debugCrystallite%element = config_debug%get_asInt('element', defaultVal=1) - debugCrystallite%ip = config_debug%get_asInt('integrationpoint', defaultVal=1) - debugCrystallite%grain = config_debug%get_asInt('grain', defaultVal=1) - - cMax = homogenization_maxNconstituents - iMax = discretization_nIPs - eMax = discretization_Nelems - - allocate(crystallite_partitionedF(3,3,cMax,iMax,eMax),source=0.0_pReal) - - allocate(crystallite_S0, & - crystallite_F0, crystallite_Fi0,crystallite_Fp0, & - crystallite_Li0,crystallite_Lp0, & - crystallite_partitionedS0, & - crystallite_partitionedF0,crystallite_partitionedFp0,crystallite_partitionedFi0, & - crystallite_partitionedLp0,crystallite_partitionedLi0, & - crystallite_S,crystallite_P, & - crystallite_Fe,crystallite_Fi,crystallite_Fp, & - crystallite_Li,crystallite_Lp, & - crystallite_subF,crystallite_subF0, & - crystallite_subFp0,crystallite_subFi0, & - source = crystallite_partitionedF) - - allocate(crystallite_dt(cMax,iMax,eMax),source=0.0_pReal) - allocate(crystallite_subdt,crystallite_subFrac,crystallite_subStep, & - source = crystallite_dt) - - allocate(crystallite_orientation(cMax,iMax,eMax)) - - allocate(crystallite_requested(cMax,iMax,eMax), source=.false.) - allocate(crystallite_converged(cMax,iMax,eMax), source=.true.) - - num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) - - num%subStepMinCryst = num_crystallite%get_asFloat ('subStepMin', defaultVal=1.0e-3_pReal) - num%subStepSizeCryst = num_crystallite%get_asFloat ('subStepSize', defaultVal=0.25_pReal) - num%stepIncreaseCryst = num_crystallite%get_asFloat ('stepIncrease', defaultVal=1.5_pReal) - num%subStepSizeLp = num_crystallite%get_asFloat ('subStepSizeLp', defaultVal=0.5_pReal) - num%subStepSizeLi = num_crystallite%get_asFloat ('subStepSizeLi', defaultVal=0.5_pReal) - num%rtol_crystalliteState = num_crystallite%get_asFloat ('rtol_State', defaultVal=1.0e-6_pReal) - num%rtol_crystalliteStress = num_crystallite%get_asFloat ('rtol_Stress', defaultVal=1.0e-6_pReal) - num%atol_crystalliteStress = num_crystallite%get_asFloat ('atol_Stress', defaultVal=1.0e-8_pReal) - num%iJacoLpresiduum = num_crystallite%get_asInt ('iJacoLpresiduum', defaultVal=1) - num%nState = num_crystallite%get_asInt ('nState', defaultVal=20) - num%nStress = num_crystallite%get_asInt ('nStress', defaultVal=40) - - if(num%subStepMinCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepMinCryst') - if(num%subStepSizeCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeCryst') - if(num%stepIncreaseCryst <= 0.0_pReal) call IO_error(301,ext_msg='stepIncreaseCryst') - - if(num%subStepSizeLp <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeLp') - if(num%subStepSizeLi <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeLi') - - if(num%rtol_crystalliteState <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteState') - if(num%rtol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteStress') - if(num%atol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='atol_crystalliteStress') - - if(num%iJacoLpresiduum < 1) call IO_error(301,ext_msg='iJacoLpresiduum') - - if(num%nState < 1) call IO_error(301,ext_msg='nState') - if(num%nStress< 1) call IO_error(301,ext_msg='nStress') - - select case(num_crystallite%get_asString('integrator',defaultVal='FPI')) - case('FPI') - integrateState => integrateStateFPI - case('Euler') - integrateState => integrateStateEuler - case('AdaptiveEuler') - integrateState => integrateStateAdaptiveEuler - case('RK4') - integrateState => integrateStateRK4 - case('RKCK45') - integrateState => integrateStateRKCK45 - case default - call IO_error(301,ext_msg='integrator') - end select - - phases => config_material%get('phase') - - allocate(output_constituent(phases%length)) - do p = 1, phases%length - phase => phases%get(p) - mech => phase%get('mechanics',defaultVal = emptyDict) -#if defined(__GFORTRAN__) - output_constituent(p)%label = output_asStrings(mech) -#else - output_constituent(p)%label = mech%get_asStrings('output',defaultVal=emptyStringArray) -#endif - enddo - -#ifdef DEBUG - if (debugCrystallite%basic) then - print'(a42,1x,i10)', ' # of elements: ', eMax - print'(a42,1x,i10)', ' # of integration points/element: ', iMax - print'(a42,1x,i10)', 'max # of constituents/integration point: ', cMax - flush(IO_STDOUT) - endif -#endif - - !$OMP PARALLEL DO PRIVATE(i,c) - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1), FEsolving_execIP(2); do c = 1, homogenization_Nconstituents(material_homogenizationAt(e)) - crystallite_Fp0(1:3,1:3,c,i,e) = material_orientation0(c,i,e)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) - crystallite_Fp0(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) & - / math_det33(crystallite_Fp0(1:3,1:3,c,i,e))**(1.0_pReal/3.0_pReal) - crystallite_Fi0(1:3,1:3,c,i,e) = constitutive_initialFi(c,i,e) - crystallite_F0(1:3,1:3,c,i,e) = math_I3 - crystallite_Fe(1:3,1:3,c,i,e) = math_inv33(matmul(crystallite_Fi0(1:3,1:3,c,i,e), & - crystallite_Fp0(1:3,1:3,c,i,e))) ! assuming that euler angles are given in internal strain free configuration - crystallite_Fp(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) - crystallite_Fi(1:3,1:3,c,i,e) = crystallite_Fi0(1:3,1:3,c,i,e) - crystallite_requested(c,i,e) = .true. - enddo; enddo - enddo - !$OMP END PARALLEL DO - - - crystallite_partitionedFp0 = crystallite_Fp0 - crystallite_partitionedFi0 = crystallite_Fi0 - crystallite_partitionedF0 = crystallite_F0 - crystallite_partitionedF = crystallite_F0 - - call crystallite_orientations() - - !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) - call constitutive_dependentState(crystallite_partitionedF0(1:3,1:3,c,i,e), & - crystallite_partitionedFp0(1:3,1:3,c,i,e), & - c,i,e) ! update dependent state variables to be consistent with basic states - enddo - enddo - enddo - !$OMP END PARALLEL DO - - -end subroutine crystallite_init - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculate stress (P) -!-------------------------------------------------------------------------------------------------- -function crystallite_stress() - - logical, dimension(discretization_nIPs,discretization_Nelems) :: crystallite_stress - real(pReal) :: & - formerSubStep - integer :: & - NiterationCrystallite, & ! number of iterations in crystallite loop - c, & !< counter in integration point component loop - i, & !< counter in integration point loop - e, & !< counter in element loop - s - logical, dimension(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems) :: todo !ToDo: need to set some values to false for different Ngrains - real(pReal), dimension(:,:,:,:,:), allocatable :: & - subLp0,& !< plastic velocity grad at start of crystallite inc - subLi0 !< intermediate velocity grad at start of crystallite inc - - todo = .false. - - subLp0 = crystallite_partitionedLp0 - subLi0 = crystallite_partitionedLi0 - -!-------------------------------------------------------------------------------------------------- -! initialize to starting condition - crystallite_subStep = 0.0_pReal - !$OMP PARALLEL DO - elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2); do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) - homogenizationRequestsCalculation: if (crystallite_requested(c,i,e)) then - plasticState (material_phaseAt(c,e))%subState0( :,material_phaseMemberAt(c,i,e)) = & - plasticState (material_phaseAt(c,e))%partitionedState0(:,material_phaseMemberAt(c,i,e)) - - do s = 1, phase_Nsources(material_phaseAt(c,e)) - sourceState(material_phaseAt(c,e))%p(s)%subState0( :,material_phaseMemberAt(c,i,e)) = & - sourceState(material_phaseAt(c,e))%p(s)%partitionedState0(:,material_phaseMemberAt(c,i,e)) - enddo - crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_partitionedFp0(1:3,1:3,c,i,e) - crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_partitionedFi0(1:3,1:3,c,i,e) - crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partitionedF0(1:3,1:3,c,i,e) - crystallite_subFrac(c,i,e) = 0.0_pReal - crystallite_subStep(c,i,e) = 1.0_pReal/num%subStepSizeCryst - todo(c,i,e) = .true. - crystallite_converged(c,i,e) = .false. ! pretend failed step of 1/subStepSizeCryst - endif homogenizationRequestsCalculation - enddo; enddo - enddo elementLooping1 - !$OMP END PARALLEL DO - - NiterationCrystallite = 0 - cutbackLooping: do while (any(todo(:,FEsolving_execIP(1):FEsolving_execIP(2),FEsolving_execELem(1):FEsolving_execElem(2)))) - NiterationCrystallite = NiterationCrystallite + 1 - -#ifdef DEBUG - if (debugCrystallite%extensive) & - print'(a,i6)', '<< CRYST stress >> crystallite iteration ',NiterationCrystallite -#endif - !$OMP PARALLEL DO PRIVATE(formerSubStep) - elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) -!-------------------------------------------------------------------------------------------------- -! wind forward - if (crystallite_converged(c,i,e)) then - formerSubStep = crystallite_subStep(c,i,e) - crystallite_subFrac(c,i,e) = crystallite_subFrac(c,i,e) + crystallite_subStep(c,i,e) - crystallite_subStep(c,i,e) = min(1.0_pReal - crystallite_subFrac(c,i,e), & - num%stepIncreaseCryst * crystallite_subStep(c,i,e)) - - todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on? - if (todo(c,i,e)) then - crystallite_subF0 (1:3,1:3,c,i,e) = crystallite_subF(1:3,1:3,c,i,e) - subLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e) - subLi0(1:3,1:3,c,i,e) = crystallite_Li (1:3,1:3,c,i,e) - crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_Fp (1:3,1:3,c,i,e) - crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_Fi (1:3,1:3,c,i,e) - plasticState( material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e)) & - = plasticState(material_phaseAt(c,e))%state( :,material_phaseMemberAt(c,i,e)) - do s = 1, phase_Nsources(material_phaseAt(c,e)) - sourceState( material_phaseAt(c,e))%p(s)%subState0(:,material_phaseMemberAt(c,i,e)) & - = sourceState(material_phaseAt(c,e))%p(s)%state( :,material_phaseMemberAt(c,i,e)) - enddo - endif - -!-------------------------------------------------------------------------------------------------- -! cut back (reduced time and restore) - else - crystallite_subStep(c,i,e) = num%subStepSizeCryst * crystallite_subStep(c,i,e) - crystallite_Fp (1:3,1:3,c,i,e) = crystallite_subFp0(1:3,1:3,c,i,e) - crystallite_Fi (1:3,1:3,c,i,e) = crystallite_subFi0(1:3,1:3,c,i,e) - crystallite_S (1:3,1:3,c,i,e) = crystallite_S0 (1:3,1:3,c,i,e) - if (crystallite_subStep(c,i,e) < 1.0_pReal) then ! actual (not initial) cutback - crystallite_Lp (1:3,1:3,c,i,e) = subLp0(1:3,1:3,c,i,e) - crystallite_Li (1:3,1:3,c,i,e) = subLi0(1:3,1:3,c,i,e) - endif - plasticState (material_phaseAt(c,e))%state( :,material_phaseMemberAt(c,i,e)) & - = plasticState(material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e)) - do s = 1, phase_Nsources(material_phaseAt(c,e)) - sourceState( material_phaseAt(c,e))%p(s)%state( :,material_phaseMemberAt(c,i,e)) & - = sourceState(material_phaseAt(c,e))%p(s)%subState0(:,material_phaseMemberAt(c,i,e)) - enddo - - ! cant restore dotState here, since not yet calculated in first cutback after initialization - todo(c,i,e) = crystallite_subStep(c,i,e) > num%subStepMinCryst ! still on track or already done (beyond repair) - endif - -!-------------------------------------------------------------------------------------------------- -! prepare for integration - if (todo(c,i,e)) then - crystallite_subF(1:3,1:3,c,i,e) = crystallite_subF0(1:3,1:3,c,i,e) & - + crystallite_subStep(c,i,e) *( crystallite_partitionedF (1:3,1:3,c,i,e) & - -crystallite_partitionedF0(1:3,1:3,c,i,e)) - crystallite_Fe(1:3,1:3,c,i,e) = matmul(crystallite_subF(1:3,1:3,c,i,e), & - math_inv33(matmul(crystallite_Fi(1:3,1:3,c,i,e), & - crystallite_Fp(1:3,1:3,c,i,e)))) - crystallite_subdt(c,i,e) = crystallite_subStep(c,i,e) * crystallite_dt(c,i,e) - crystallite_converged(c,i,e) = .false. - call integrateState(c,i,e) - endif - - enddo - enddo - enddo elementLooping3 - !$OMP END PARALLEL DO - -!-------------------------------------------------------------------------------------------------- -! integrate --- requires fully defined state array (basic + dependent state) - where(.not. crystallite_converged .and. crystallite_subStep > num%subStepMinCryst) & ! do not try non-converged but fully cutbacked any further - todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation - enddo cutbackLooping - -! return whether converged or not - crystallite_stress = .false. - elementLooping5: do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - crystallite_stress(i,e) = all(crystallite_converged(:,i,e)) - enddo - enddo elementLooping5 - -end function crystallite_stress - - -!-------------------------------------------------------------------------------------------------- -!> @brief Backup data for homog cutback. -!-------------------------------------------------------------------------------------------------- -subroutine crystallite_initializeRestorationPoints(i,e) - - integer, intent(in) :: & - i, & !< integration point number - e !< element number - integer :: & - c, & !< constituent number - s - - do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) - crystallite_partitionedFp0(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) - crystallite_partitionedLp0(1:3,1:3,c,i,e) = crystallite_Lp0(1:3,1:3,c,i,e) - crystallite_partitionedFi0(1:3,1:3,c,i,e) = crystallite_Fi0(1:3,1:3,c,i,e) - crystallite_partitionedLi0(1:3,1:3,c,i,e) = crystallite_Li0(1:3,1:3,c,i,e) - crystallite_partitionedF0(1:3,1:3,c,i,e) = crystallite_F0(1:3,1:3,c,i,e) - crystallite_partitionedS0(1:3,1:3,c,i,e) = crystallite_S0(1:3,1:3,c,i,e) - - plasticState(material_phaseAt(c,e))%partitionedState0(:,material_phasememberAt(c,i,e)) = & - plasticState(material_phaseAt(c,e))%state0( :,material_phasememberAt(c,i,e)) - do s = 1, phase_Nsources(material_phaseAt(c,e)) - sourceState(material_phaseAt(c,e))%p(s)%partitionedState0(:,material_phasememberAt(c,i,e)) = & - sourceState(material_phaseAt(c,e))%p(s)%state0( :,material_phasememberAt(c,i,e)) - enddo - enddo - -end subroutine crystallite_initializeRestorationPoints - - -!-------------------------------------------------------------------------------------------------- -!> @brief Wind homog inc forward. -!-------------------------------------------------------------------------------------------------- -subroutine crystallite_windForward(i,e) - - integer, intent(in) :: & - i, & !< integration point number - e !< element number - integer :: & - c, & !< constituent number - s - - do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) - crystallite_partitionedF0 (1:3,1:3,c,i,e) = crystallite_partitionedF(1:3,1:3,c,i,e) - crystallite_partitionedFp0(1:3,1:3,c,i,e) = crystallite_Fp (1:3,1:3,c,i,e) - crystallite_partitionedLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e) - crystallite_partitionedFi0(1:3,1:3,c,i,e) = crystallite_Fi (1:3,1:3,c,i,e) - crystallite_partitionedLi0(1:3,1:3,c,i,e) = crystallite_Li (1:3,1:3,c,i,e) - crystallite_partitionedS0 (1:3,1:3,c,i,e) = crystallite_S (1:3,1:3,c,i,e) - - plasticState (material_phaseAt(c,e))%partitionedState0(:,material_phasememberAt(c,i,e)) = & - plasticState (material_phaseAt(c,e))%state (:,material_phasememberAt(c,i,e)) - do s = 1, phase_Nsources(material_phaseAt(c,e)) - sourceState(material_phaseAt(c,e))%p(s)%partitionedState0(:,material_phasememberAt(c,i,e)) = & - sourceState(material_phaseAt(c,e))%p(s)%state (:,material_phasememberAt(c,i,e)) - enddo - enddo - -end subroutine crystallite_windForward - - -!-------------------------------------------------------------------------------------------------- -!> @brief Restore data after homog cutback. -!-------------------------------------------------------------------------------------------------- -subroutine crystallite_restore(i,e,includeL) - - integer, intent(in) :: & - i, & !< integration point number - e !< element number - logical, intent(in) :: & - includeL !< protect agains fake cutback - integer :: & - c, & !< constituent number - s - - do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) - if (includeL) then - crystallite_Lp(1:3,1:3,c,i,e) = crystallite_partitionedLp0(1:3,1:3,c,i,e) - crystallite_Li(1:3,1:3,c,i,e) = crystallite_partitionedLi0(1:3,1:3,c,i,e) - endif ! maybe protecting everything from overwriting makes more sense - crystallite_Fp(1:3,1:3,c,i,e) = crystallite_partitionedFp0(1:3,1:3,c,i,e) - crystallite_Fi(1:3,1:3,c,i,e) = crystallite_partitionedFi0(1:3,1:3,c,i,e) - crystallite_S (1:3,1:3,c,i,e) = crystallite_partitionedS0 (1:3,1:3,c,i,e) - - plasticState (material_phaseAt(c,e))%state( :,material_phasememberAt(c,i,e)) = & - plasticState (material_phaseAt(c,e))%partitionedState0(:,material_phasememberAt(c,i,e)) - do s = 1, phase_Nsources(material_phaseAt(c,e)) - sourceState(material_phaseAt(c,e))%p(s)%state( :,material_phasememberAt(c,i,e)) = & - sourceState(material_phaseAt(c,e))%p(s)%partitionedState0(:,material_phasememberAt(c,i,e)) - enddo - enddo - -end subroutine crystallite_restore - - -!-------------------------------------------------------------------------------------------------- -!> @brief Calculate tangent (dPdF). -!-------------------------------------------------------------------------------------------------- -function crystallite_stressTangent(c,i,e) result(dPdF) - - real(pReal), dimension(3,3,3,3) :: dPdF - integer, intent(in) :: & - c, & !< counter in constituent loop - i, & !< counter in integration point loop - e !< counter in element loop - integer :: & - o, & - p - - real(pReal), dimension(3,3) :: devNull, & - invSubFp0,invSubFi0,invFp,invFi, & - temp_33_1, temp_33_2, temp_33_3, temp_33_4 - real(pReal), dimension(3,3,3,3) :: dSdFe, & - dSdF, & - dSdFi, & - dLidS, & ! tangent in lattice configuration - dLidFi, & - dLpdS, & - dLpdFi, & - dFidS, & - dFpinvdF, & - rhs_3333, & - lhs_3333, & - temp_3333 - real(pReal), dimension(9,9):: temp_99 - logical :: error - - - call constitutive_SandItsTangents(devNull,dSdFe,dSdFi, & - crystallite_Fe(1:3,1:3,c,i,e), & - crystallite_Fi(1:3,1:3,c,i,e),c,i,e) - call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, & - crystallite_S (1:3,1:3,c,i,e), & - crystallite_Fi(1:3,1:3,c,i,e), & - c,i,e) - - invFp = math_inv33(crystallite_Fp(1:3,1:3,c,i,e)) - invFi = math_inv33(crystallite_Fi(1:3,1:3,c,i,e)) - invSubFp0 = math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)) - invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,c,i,e)) - - if (sum(abs(dLidS)) < tol_math_check) then - dFidS = 0.0_pReal - else - lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal - do o=1,3; do p=1,3 - lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) & - + crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) - lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) & - + invFi*invFi(p,o) - rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) & - - crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidS(1:3,1:3,o,p)) - enddo; enddo - call math_invert(temp_99,error,math_3333to99(lhs_3333)) - if (error) then - call IO_warning(warning_ID=600,el=e,ip=i,g=c, & - ext_msg='inversion error in analytic tangent calculation') - dFidS = 0.0_pReal - else - dFidS = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333) - endif - dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS - endif - - call constitutive_LpAndItsTangents(devNull,dLpdS,dLpdFi, & - crystallite_S (1:3,1:3,c,i,e), & - crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration - dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS - -!-------------------------------------------------------------------------------------------------- -! calculate dSdF - temp_33_1 = transpose(matmul(invFp,invFi)) - temp_33_2 = matmul(crystallite_subF(1:3,1:3,c,i,e),invSubFp0) - temp_33_3 = matmul(matmul(crystallite_subF(1:3,1:3,c,i,e),invFp), invSubFi0) - - do o=1,3; do p=1,3 - rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1) - temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) & - + matmul(temp_33_3,dLidS(1:3,1:3,p,o)) - enddo; enddo - lhs_3333 = crystallite_subdt(c,i,e)*math_mul3333xx3333(dSdFe,temp_3333) & - + math_mul3333xx3333(dSdFi,dFidS) - - call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333)) - if (error) then - call IO_warning(warning_ID=600,el=e,ip=i,g=c, & - ext_msg='inversion error in analytic tangent calculation') - dSdF = rhs_3333 - else - dSdF = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333) - endif - -!-------------------------------------------------------------------------------------------------- -! calculate dFpinvdF - temp_3333 = math_mul3333xx3333(dLpdS,dSdF) - do o=1,3; do p=1,3 - dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(c,i,e) & - * matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) - enddo; enddo - -!-------------------------------------------------------------------------------------------------- -! assemble dPdF - temp_33_1 = matmul(crystallite_S(1:3,1:3,c,i,e),transpose(invFp)) - temp_33_2 = matmul(invFp,temp_33_1) - temp_33_3 = matmul(crystallite_subF(1:3,1:3,c,i,e),invFp) - temp_33_4 = matmul(temp_33_3,crystallite_S(1:3,1:3,c,i,e)) - - dPdF = 0.0_pReal - do p=1,3 - dPdF(p,1:3,p,1:3) = transpose(temp_33_2) - enddo - do o=1,3; do p=1,3 - dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) & - + matmul(matmul(crystallite_subF(1:3,1:3,c,i,e), & - dFpinvdF(1:3,1:3,p,o)),temp_33_1) & - + matmul(matmul(temp_33_3,dSdF(1:3,1:3,p,o)), & - transpose(invFp)) & - + matmul(temp_33_4,transpose(dFpinvdF(1:3,1:3,p,o))) - enddo; enddo - -end function crystallite_stressTangent - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates orientations -!-------------------------------------------------------------------------------------------------- -subroutine crystallite_orientations - - integer & - c, & !< counter in integration point component loop - i, & !< counter in integration point loop - e !< counter in element loop - - !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) - call crystallite_orientation(c,i,e)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,c,i,e)))) - enddo; enddo; enddo - !$OMP END PARALLEL DO - - nonlocalPresent: if (any(plasticState%nonlocal)) then - !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) - if (plasticState(material_phaseAt(1,e))%nonlocal) then - do i = FEsolving_execIP(1),FEsolving_execIP(2) - call plastic_nonlocal_updateCompatibility(crystallite_orientation, & - phase_plasticityInstance(material_phaseAt(1,e)),i,e) - enddo - endif - enddo - !$OMP END PARALLEL DO - endif nonlocalPresent - -end subroutine crystallite_orientations - - -!-------------------------------------------------------------------------------------------------- -!> @brief Map 2nd order tensor to reference config -!-------------------------------------------------------------------------------------------------- -function crystallite_push33ToRef(ipc,ip,el, tensor33) - - real(pReal), dimension(3,3) :: crystallite_push33ToRef - real(pReal), dimension(3,3), intent(in) :: tensor33 - real(pReal), dimension(3,3) :: T - integer, intent(in):: & - el, & - ip, & - ipc - - T = matmul(material_orientation0(ipc,ip,el)%asMatrix(), & ! ToDo: initial orientation correct? - transpose(math_inv33(crystallite_subF(1:3,1:3,ipc,ip,el)))) - crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) - -end function crystallite_push33ToRef - - -!-------------------------------------------------------------------------------------------------- -!> @brief writes crystallite results to HDF5 output file -!-------------------------------------------------------------------------------------------------- -subroutine crystallite_results - - integer :: p,o - real(pReal), allocatable, dimension(:,:,:) :: selected_tensors - real(pReal), allocatable, dimension(:,:) :: selected_rotations - character(len=:), allocatable :: group,structureLabel - - do p=1,size(material_name_phase) - group = trim('current/phase')//'/'//trim(material_name_phase(p))//'/mechanics' - - call results_closeGroup(results_addGroup(group)) - - do o = 1, size(output_constituent(p)%label) - select case (output_constituent(p)%label(o)) - case('F') - selected_tensors = select_tensors(crystallite_partitionedF,p) - call results_writeDataset(group,selected_tensors,output_constituent(p)%label(o),& - 'deformation gradient','1') - case('F_e') - selected_tensors = select_tensors(crystallite_Fe,p) - call results_writeDataset(group,selected_tensors,output_constituent(p)%label(o),& - 'elastic deformation gradient','1') - case('F_p') - selected_tensors = select_tensors(crystallite_Fp,p) - call results_writeDataset(group,selected_tensors,output_constituent(p)%label(o),& - 'plastic deformation gradient','1') - case('F_i') - selected_tensors = select_tensors(crystallite_Fi,p) - call results_writeDataset(group,selected_tensors,output_constituent(p)%label(o),& - 'inelastic deformation gradient','1') - case('L_p') - selected_tensors = select_tensors(crystallite_Lp,p) - call results_writeDataset(group,selected_tensors,output_constituent(p)%label(o),& - 'plastic velocity gradient','1/s') - case('L_i') - selected_tensors = select_tensors(crystallite_Li,p) - call results_writeDataset(group,selected_tensors,output_constituent(p)%label(o),& - 'inelastic velocity gradient','1/s') - case('P') - selected_tensors = select_tensors(crystallite_P,p) - call results_writeDataset(group,selected_tensors,output_constituent(p)%label(o),& - 'First Piola-Kirchhoff stress','Pa') - case('S') - selected_tensors = select_tensors(crystallite_S,p) - call results_writeDataset(group,selected_tensors,output_constituent(p)%label(o),& - 'Second Piola-Kirchhoff stress','Pa') - case('O') - select case(lattice_structure(p)) - case(lattice_ISO_ID) - structureLabel = 'aP' - case(lattice_FCC_ID) - structureLabel = 'cF' - case(lattice_BCC_ID) - structureLabel = 'cI' - case(lattice_BCT_ID) - structureLabel = 'tI' - case(lattice_HEX_ID) - structureLabel = 'hP' - case(lattice_ORT_ID) - structureLabel = 'oP' - end select - selected_rotations = select_rotations(crystallite_orientation,p) - call results_writeDataset(group,selected_rotations,output_constituent(p)%label(o),& - 'crystal orientation as quaternion','q_0 ') - call results_addAttribute('Lattice',structureLabel,group//'/'//output_constituent(p)%label(o)) - end select - enddo - enddo - - contains - - !------------------------------------------------------------------------------------------------ - !> @brief select tensors for output - !------------------------------------------------------------------------------------------------ - function select_tensors(dataset,instance) - - integer, intent(in) :: instance - real(pReal), dimension(:,:,:,:,:), intent(in) :: dataset - real(pReal), allocatable, dimension(:,:,:) :: select_tensors - integer :: e,i,c,j - - allocate(select_tensors(3,3,count(material_phaseAt==instance)*discretization_nIPs)) - - j=0 - do e = 1, size(material_phaseAt,2) - do i = 1, discretization_nIPs - do c = 1, size(material_phaseAt,1) !ToDo: this needs to be changed for varying Ngrains - if (material_phaseAt(c,e) == instance) then - j = j + 1 - select_tensors(1:3,1:3,j) = dataset(1:3,1:3,c,i,e) - endif - enddo - enddo - enddo - - end function select_tensors - - -!-------------------------------------------------------------------------------------------------- -!> @brief select rotations for output -!-------------------------------------------------------------------------------------------------- - function select_rotations(dataset,instance) - - integer, intent(in) :: instance - type(rotation), dimension(:,:,:), intent(in) :: dataset - real(pReal), allocatable, dimension(:,:) :: select_rotations - integer :: e,i,c,j - - allocate(select_rotations(4,count(material_phaseAt==instance)*homogenization_maxNconstituents*discretization_nIPs)) - - j=0 - do e = 1, size(material_phaseAt,2) - do i = 1, discretization_nIPs - do c = 1, size(material_phaseAt,1) !ToDo: this needs to be changed for varying Ngrains - if (material_phaseAt(c,e) == instance) then - j = j + 1 - select_rotations(1:4,j) = dataset(c,i,e)%asQuaternion() - endif - enddo - enddo - enddo - - end function select_rotations - -end subroutine crystallite_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(ipc,ip,el,timeFraction) result(broken) - - integer, intent(in):: el, & ! element index - ip, & ! integration point index - ipc ! grain index - real(pReal), optional, intent(in) :: timeFraction ! fraction of timestep - - real(pReal), dimension(3,3):: F, & ! deformation gradient at end of timestep - Fp_new, & ! plastic deformation gradient at end of timestep - invFp_new, & ! inverse of Fp_new - invFp_current, & ! inverse of Fp_current - Lpguess, & ! current guess for plastic velocity gradient - Lpguess_old, & ! known last good guess for plastic velocity gradient - Lp_constitutive, & ! plastic velocity gradient resulting from constitutive law - residuumLp, & ! current residuum of plastic velocity gradient - residuumLp_old, & ! last residuum of plastic velocity gradient - deltaLp, & ! direction of next guess - Fi_new, & ! gradient of intermediate deformation stages - invFi_new, & - invFi_current, & ! inverse of Fi_current - Liguess, & ! current guess for intermediate velocity gradient - Liguess_old, & ! known last good guess for intermediate velocity gradient - Li_constitutive, & ! intermediate velocity gradient resulting from constitutive law - residuumLi, & ! current residuum of intermediate velocity gradient - residuumLi_old, & ! last residuum of intermediate velocity gradient - deltaLi, & ! direction of next guess - Fe, & ! elastic deformation gradient - S, & ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration - A, & - B, & - temp_33 - real(pReal), dimension(9) :: temp_9 ! needed for matrix inversion by LAPACK - integer, dimension(9) :: devNull_9 ! needed for matrix inversion by LAPACK - real(pReal), dimension(9,9) :: dRLp_dLp, & ! partial derivative of residuum (Jacobian for Newton-Raphson scheme) - dRLi_dLi ! partial derivative of residuumI (Jacobian for Newton-Raphson scheme) - real(pReal), dimension(3,3,3,3):: dS_dFe, & ! partial derivative of 2nd Piola-Kirchhoff stress - dS_dFi, & - dFe_dLp, & ! partial derivative of elastic deformation gradient - dFe_dLi, & - dFi_dLi, & - dLp_dFi, & - dLi_dFi, & - dLp_dS, & - dLi_dS - real(pReal) steplengthLp, & - steplengthLi, & - dt, & ! time increment - atol_Lp, & - atol_Li, & - devNull - integer NiterationStressLp, & ! number of stress integrations - NiterationStressLi, & ! number of inner stress integrations - ierr, & ! error indicator for LAPACK - o, & - p, & - jacoCounterLp, & - jacoCounterLi ! counters to check for Jacobian update - logical :: error,broken - - broken = .true. - - if (present(timeFraction)) then - dt = crystallite_subdt(ipc,ip,el) * timeFraction - F = crystallite_subF0(1:3,1:3,ipc,ip,el) & - + (crystallite_subF(1:3,1:3,ipc,ip,el) - crystallite_subF0(1:3,1:3,ipc,ip,el)) * timeFraction - else - dt = crystallite_subdt(ipc,ip,el) - F = crystallite_subF(1:3,1:3,ipc,ip,el) - endif - - call constitutive_dependentState(crystallite_partitionedF(1:3,1:3,ipc,ip,el), & - crystallite_Fp(1:3,1:3,ipc,ip,el),ipc,ip,el) - - Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! take as first guess - Liguess = crystallite_Li(1:3,1:3,ipc,ip,el) ! take as first guess - - call math_invert33(invFp_current,devNull,error,crystallite_subFp0(1:3,1:3,ipc,ip,el)) - if (error) return ! error - call math_invert33(invFi_current,devNull,error,crystallite_subFi0(1:3,1:3,ipc,ip,el)) - if (error) return ! error - - A = matmul(F,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp - - jacoCounterLi = 0 - steplengthLi = 1.0_pReal - residuumLi_old = 0.0_pReal - Liguess_old = Liguess - - NiterationStressLi = 0 - LiLoop: do - NiterationStressLi = NiterationStressLi + 1 - if (NiterationStressLi>num%nStress) return ! error - - invFi_new = matmul(invFi_current,math_I3 - dt*Liguess) - Fi_new = math_inv33(invFi_new) - - jacoCounterLp = 0 - steplengthLp = 1.0_pReal - residuumLp_old = 0.0_pReal - Lpguess_old = Lpguess - - NiterationStressLp = 0 - LpLoop: do - NiterationStressLp = NiterationStressLp + 1 - if (NiterationStressLp>num%nStress) return ! error - - B = math_I3 - dt*Lpguess - Fe = matmul(matmul(A,B), invFi_new) - call constitutive_SandItsTangents(S, dS_dFe, dS_dFi, & - Fe, Fi_new, ipc, ip, el) - - call constitutive_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & - S, Fi_new, ipc, ip, el) - - !* update current residuum and check for convergence of loop - atol_Lp = max(num%rtol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error - num%atol_crystalliteStress) ! minimum lower cutoff - residuumLp = Lpguess - Lp_constitutive - - if (any(IEEE_is_NaN(residuumLp))) then - return ! error - elseif (norm2(residuumLp) < atol_Lp) then ! converged if below absolute tolerance - exit LpLoop - elseif (NiterationStressLp == 1 .or. norm2(residuumLp) < norm2(residuumLp_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)... - residuumLp_old = residuumLp ! ...remember old values and... - Lpguess_old = Lpguess - steplengthLp = 1.0_pReal ! ...proceed with normal step length (calculate new search direction) - else ! not converged and residuum not improved... - steplengthLp = num%subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction - Lpguess = Lpguess_old & - + deltaLp * stepLengthLp - cycle LpLoop - endif - - calculateJacobiLi: if (mod(jacoCounterLp, num%iJacoLpresiduum) == 0) then - jacoCounterLp = jacoCounterLp + 1 - - do o=1,3; do p=1,3 - dFe_dLp(o,1:3,p,1:3) = - dt * A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) - enddo; enddo - dRLp_dLp = math_eye(9) & - - math_3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dS,dS_dFe),dFe_dLp)) - temp_9 = math_33to9(residuumLp) - call dgesv(9,1,dRLp_dLp,9,devNull_9,temp_9,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp - if (ierr /= 0) return ! error - deltaLp = - math_9to33(temp_9) - endif calculateJacobiLi - - Lpguess = Lpguess & - + deltaLp * steplengthLp - enddo LpLoop - - call constitutive_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, & - S, Fi_new, ipc, ip, el) - - !* update current residuum and check for convergence of loop - atol_Li = max(num%rtol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error - num%atol_crystalliteStress) ! minimum lower cutoff - residuumLi = Liguess - Li_constitutive - if (any(IEEE_is_NaN(residuumLi))) then - return ! error - elseif (norm2(residuumLi) < atol_Li) then ! converged if below absolute tolerance - exit LiLoop - elseif (NiterationStressLi == 1 .or. norm2(residuumLi) < norm2(residuumLi_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)... - residuumLi_old = residuumLi ! ...remember old values and... - Liguess_old = Liguess - steplengthLi = 1.0_pReal ! ...proceed with normal step length (calculate new search direction) - else ! not converged and residuum not improved... - steplengthLi = num%subStepSizeLi * steplengthLi ! ...try with smaller step length in same direction - Liguess = Liguess_old & - + deltaLi * steplengthLi - cycle LiLoop - endif - - calculateJacobiLp: if (mod(jacoCounterLi, num%iJacoLpresiduum) == 0) then - jacoCounterLi = jacoCounterLi + 1 - - temp_33 = matmul(matmul(A,B),invFi_current) - do o=1,3; do p=1,3 - dFe_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*temp_33 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) - dFi_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*invFi_current - enddo; enddo - do o=1,3; do p=1,3 - dFi_dLi(1:3,1:3,o,p) = matmul(matmul(Fi_new,dFi_dLi(1:3,1:3,o,p)),Fi_new) - enddo; enddo - dRLi_dLi = math_eye(9) & - - math_3333to99(math_mul3333xx3333(dLi_dS, math_mul3333xx3333(dS_dFe, dFe_dLi) & - + math_mul3333xx3333(dS_dFi, dFi_dLi))) & - - math_3333to99(math_mul3333xx3333(dLi_dFi, dFi_dLi)) - temp_9 = math_33to9(residuumLi) - call dgesv(9,1,dRLi_dLi,9,devNull_9,temp_9,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li - if (ierr /= 0) return ! error - deltaLi = - math_9to33(temp_9) - endif calculateJacobiLp - - Liguess = Liguess & - + deltaLi * steplengthLi - enddo LiLoop - - invFp_new = matmul(invFp_current,B) - call math_invert33(Fp_new,devNull,error,invFp_new) - if (error) return ! error - - crystallite_P (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) - crystallite_S (1:3,1:3,ipc,ip,el) = S - crystallite_Lp (1:3,1:3,ipc,ip,el) = Lpguess - crystallite_Li (1:3,1:3,ipc,ip,el) = Liguess - crystallite_Fp (1:3,1:3,ipc,ip,el) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize - crystallite_Fi (1:3,1:3,ipc,ip,el) = Fi_new - crystallite_Fe (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),invFi_new) - broken = .false. - -end function integrateStress - - -!-------------------------------------------------------------------------------------------------- -!> @brief integrate stress, state with adaptive 1st order explicit Euler method -!> using Fixed Point Iteration to adapt the stepsize -!-------------------------------------------------------------------------------------------------- -subroutine integrateStateFPI(g,i,e) - - integer, intent(in) :: & - e, & !< element index in element loop - i, & !< integration point index in ip loop - g !< grain index in grain loop - integer :: & - NiterationState, & !< number of iterations in state loop - p, & - c, & - s, & - size_pl - integer, dimension(maxval(phase_Nsources)) :: & - size_so - real(pReal) :: & - zeta - real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: & - r ! state residuum - real(pReal), dimension(constitutive_plasticity_maxSizeDotState,2) :: & - plastic_dotState - real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState - logical :: & - broken - - p = material_phaseAt(g,e) - c = material_phaseMemberAt(g,i,e) - - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partitionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partitionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) return - - size_pl = plasticState(p)%sizeDotState - plasticState(p)%state(1:size_pl,c) = plasticState(p)%subState0(1:size_pl,c) & - + plasticState(p)%dotState (1:size_pl,c) & - * crystallite_subdt(g,i,e) - plastic_dotState(1:size_pl,2) = 0.0_pReal - do s = 1, phase_Nsources(p) - size_so(s) = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%subState0(1:size_so(s),c) & - + sourceState(p)%p(s)%dotState (1:size_so(s),c) & - * crystallite_subdt(g,i,e) - source_dotState(1:size_so(s),2,s) = 0.0_pReal - enddo - - iteration: do NiterationState = 1, num%nState - - if(nIterationState > 1) plastic_dotState(1:size_pl,2) = plastic_dotState(1:size_pl,1) - plastic_dotState(1:size_pl,1) = plasticState(p)%dotState(:,c) - do s = 1, phase_Nsources(p) - if(nIterationState > 1) source_dotState(1:size_so(s),2,s) = source_dotState(1:size_so(s),1,s) - source_dotState(1:size_so(s),1,s) = sourceState(p)%p(s)%dotState(:,c) - enddo - - broken = integrateStress(g,i,e) - if(broken) exit iteration - - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partitionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partitionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) exit iteration - - zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState(1:size_pl,1),& - plastic_dotState(1:size_pl,2)) - plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & - + plastic_dotState(1:size_pl,1) * (1.0_pReal - zeta) - r(1:size_pl) = plasticState(p)%state (1:size_pl,c) & - - plasticState(p)%subState0(1:size_pl,c) & - - plasticState(p)%dotState (1:size_pl,c) * crystallite_subdt(g,i,e) - plasticState(p)%state(1:size_pl,c) = plasticState(p)%state(1:size_pl,c) & - - r(1:size_pl) - crystallite_converged(g,i,e) = converged(r(1:size_pl), & - plasticState(p)%state(1:size_pl,c), & - plasticState(p)%atol(1:size_pl)) - do s = 1, phase_Nsources(p) - zeta = damper(sourceState(p)%p(s)%dotState(:,c), & - source_dotState(1:size_so(s),1,s),& - source_dotState(1:size_so(s),2,s)) - sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta & - + source_dotState(1:size_so(s),1,s)* (1.0_pReal - zeta) - r(1:size_so(s)) = sourceState(p)%p(s)%state (1:size_so(s),c) & - - sourceState(p)%p(s)%subState0(1:size_so(s),c) & - - sourceState(p)%p(s)%dotState (1:size_so(s),c) * crystallite_subdt(g,i,e) - sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%state(1:size_so(s),c) & - - r(1:size_so(s)) - crystallite_converged(g,i,e) = & - crystallite_converged(g,i,e) .and. converged(r(1:size_so(s)), & - sourceState(p)%p(s)%state(1:size_so(s),c), & - sourceState(p)%p(s)%atol(1:size_so(s))) - enddo - - if(crystallite_converged(g,i,e)) then - broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - exit iteration - endif - - enddo iteration - - - contains - - !-------------------------------------------------------------------------------------------------- - !> @brief calculate the damping for correction of state and dot state - !-------------------------------------------------------------------------------------------------- - real(pReal) pure function damper(current,previous,previous2) - - real(pReal), dimension(:), intent(in) ::& - current, previous, previous2 - - real(pReal) :: dot_prod12, dot_prod22 - - dot_prod12 = dot_product(current - previous, previous - previous2) - dot_prod22 = dot_product(previous - previous2, previous - previous2) - if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then - damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) - else - damper = 1.0_pReal - endif - - end function damper - -end subroutine integrateStateFPI - - -!-------------------------------------------------------------------------------------------------- -!> @brief integrate state with 1st order explicit Euler method -!-------------------------------------------------------------------------------------------------- -subroutine integrateStateEuler(g,i,e) - - integer, intent(in) :: & - e, & !< element index in element loop - i, & !< integration point index in ip loop - g !< grain index in grain loop - integer :: & - p, & - c, & - s, & - sizeDotState - logical :: & - broken - - p = material_phaseAt(g,e) - c = material_phaseMemberAt(g,i,e) - - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partitionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partitionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) return - - sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - enddo - - broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - if(broken) return - - broken = integrateStress(g,i,e) - crystallite_converged(g,i,e) = .not. broken - -end subroutine integrateStateEuler - - -!-------------------------------------------------------------------------------------------------- -!> @brief integrate stress, state with 1st order Euler method with adaptive step size -!-------------------------------------------------------------------------------------------------- -subroutine integrateStateAdaptiveEuler(g,i,e) - - integer, intent(in) :: & - e, & !< element index in element loop - i, & !< integration point index in ip loop - g !< grain index in grain loop - integer :: & - p, & - c, & - s, & - sizeDotState - logical :: & - broken - - real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic - real(pReal), dimension(constitutive_source_maxSizeDotState,maxval(phase_Nsources)) :: residuum_source - - - p = material_phaseAt(g,e) - c = material_phaseMemberAt(g,i,e) - - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partitionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partitionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) return - - sizeDotState = plasticState(p)%sizeDotState - - residuum_plastic(1:sizeDotState) = - plasticState(p)%dotstate(1:sizeDotState,c) * 0.5_pReal * crystallite_subdt(g,i,e) - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - - residuum_source(1:sizeDotState,s) = - sourceState(p)%p(s)%dotstate(1:sizeDotState,c) & - * 0.5_pReal * crystallite_subdt(g,i,e) - sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) - enddo - - broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - if(broken) return - - broken = integrateStress(g,i,e) - if(broken) return - - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partitionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partitionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) return - - - sizeDotState = plasticState(p)%sizeDotState - crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState) & - + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e), & - plasticState(p)%state(1:sizeDotState,c), & - plasticState(p)%atol(1:sizeDotState)) - - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - crystallite_converged(g,i,e) = & - crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState,s) & - + 0.5_pReal*sourceState(p)%p(s)%dotState(:,c)*crystallite_subdt(g,i,e), & - sourceState(p)%p(s)%state(1:sizeDotState,c), & - sourceState(p)%p(s)%atol(1:sizeDotState)) - enddo - -end subroutine integrateStateAdaptiveEuler - - -!--------------------------------------------------------------------------------------------------- -!> @brief Integrate state (including stress integration) with the classic Runge Kutta method -!--------------------------------------------------------------------------------------------------- -subroutine integrateStateRK4(g,i,e) - - integer, intent(in) :: g,i,e - - real(pReal), dimension(3,3), parameter :: & - A = reshape([& - 0.5_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.5_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal],& - shape(A)) - real(pReal), dimension(3), parameter :: & - C = [0.5_pReal, 0.5_pReal, 1.0_pReal] - real(pReal), dimension(4), parameter :: & - 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] - - call integrateStateRK(g,i,e,A,B,C) - -end subroutine integrateStateRK4 - - -!--------------------------------------------------------------------------------------------------- -!> @brief Integrate state (including stress integration) with the Cash-Carp method -!--------------------------------------------------------------------------------------------------- -subroutine integrateStateRKCK45(g,i,e) - - integer, intent(in) :: g,i,e - - real(pReal), dimension(5,5), parameter :: & - A = reshape([& - 1._pReal/5._pReal, .0_pReal, .0_pReal, .0_pReal, .0_pReal, & - 3._pReal/40._pReal, 9._pReal/40._pReal, .0_pReal, .0_pReal, .0_pReal, & - 3_pReal/10._pReal, -9._pReal/10._pReal, 6._pReal/5._pReal, .0_pReal, .0_pReal, & - -11._pReal/54._pReal, 5._pReal/2._pReal, -70.0_pReal/27.0_pReal, 35.0_pReal/27.0_pReal, .0_pReal, & - 1631._pReal/55296._pReal,175._pReal/512._pReal,575._pReal/13824._pReal,44275._pReal/110592._pReal,253._pReal/4096._pReal],& - shape(A)) - real(pReal), dimension(5), parameter :: & - C = [0.2_pReal, 0.3_pReal, 0.6_pReal, 1.0_pReal, 0.875_pReal] - real(pReal), dimension(6), parameter :: & - B = & - [37.0_pReal/378.0_pReal, .0_pReal, 250.0_pReal/621.0_pReal, & - 125.0_pReal/594.0_pReal, .0_pReal, 512.0_pReal/1771.0_pReal], & - DB = B - & - [2825.0_pReal/27648.0_pReal, .0_pReal, 18575.0_pReal/48384.0_pReal,& - 13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 1._pReal/4._pReal] - - call integrateStateRK(g,i,e,A,B,C,DB) - -end subroutine integrateStateRKCK45 - - -!-------------------------------------------------------------------------------------------------- -!> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an -!! embedded explicit Runge-Kutta method -!-------------------------------------------------------------------------------------------------- -subroutine integrateStateRK(g,i,e,A,B,CC,DB) - - - real(pReal), dimension(:,:), intent(in) :: A - real(pReal), dimension(:), intent(in) :: B, CC - real(pReal), dimension(:), intent(in), optional :: DB - - integer, intent(in) :: & - e, & !< element index in element loop - i, & !< integration point index in ip loop - g !< grain index in grain loop - integer :: & - stage, & ! stage index in integration stage loop - n, & - p, & - c, & - s, & - sizeDotState - logical :: & - broken - real(pReal), dimension(constitutive_source_maxSizeDotState,size(B),maxval(phase_Nsources)) :: source_RKdotState - real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState - - p = material_phaseAt(g,e) - c = material_phaseMemberAt(g,i,e) - - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partitionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partitionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) return - - do stage = 1,size(A,1) - sizeDotState = plasticState(p)%sizeDotState - plastic_RKdotState(1:sizeDotState,stage) = plasticState(p)%dotState(:,c) - plasticState(p)%dotState(:,c) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - source_RKdotState(1:sizeDotState,stage,s) = sourceState(p)%p(s)%dotState(:,c) - sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * source_RKdotState(1:sizeDotState,1,s) - enddo - - do n = 2, stage - sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) & - + A(n,stage) * plastic_RKdotState(1:sizeDotState,n) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) & - + A(n,stage) * source_RKdotState(1:sizeDotState,n,s) - enddo - enddo - - sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - enddo - - broken = integrateStress(g,i,e,CC(stage)) - if(broken) exit - - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partitionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partitionedFp0, & - crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c) - if(broken) exit - - enddo - if(broken) return - - sizeDotState = plasticState(p)%sizeDotState - - plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (p)%dotState(:,c) - plasticState(p)%dotState(:,c) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B) - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - if(present(DB)) & - broken = .not. converged( matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) & - * crystallite_subdt(g,i,e), & - plasticState(p)%state(1:sizeDotState,c), & - plasticState(p)%atol(1:sizeDotState)) - - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - - source_RKdotState(1:sizeDotState,size(B),s) = sourceState(p)%p(s)%dotState(:,c) - sourceState(p)%p(s)%dotState(:,c) = matmul(source_RKdotState(1:sizeDotState,1:size(B),s),B) - sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - if(present(DB)) & - broken = broken .or. .not. converged(matmul(source_RKdotState(1:sizeDotState,1:size(DB),s),DB) & - * crystallite_subdt(g,i,e), & - sourceState(p)%p(s)%state(1:sizeDotState,c), & - sourceState(p)%p(s)%atol(1:sizeDotState)) - enddo - if(broken) return - - broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - if(broken) return - - broken = integrateStress(g,i,e) - crystallite_converged(g,i,e) = .not. broken - - -end subroutine integrateStateRK - - -!-------------------------------------------------------------------------------------------------- -!> @brief determines whether a point is converged -!-------------------------------------------------------------------------------------------------- -logical pure function converged(residuum,state,atol) - - real(pReal), intent(in), dimension(:) ::& - residuum, state, atol - real(pReal) :: & - rTol - - rTol = num%rTol_crystalliteState - - converged = all(abs(residuum) <= max(atol, rtol*abs(state))) - -end function converged - - -!-------------------------------------------------------------------------------------------------- -!> @brief Write current restart information (Field and constitutive data) to file. -! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively -!-------------------------------------------------------------------------------------------------- -subroutine crystallite_restartWrite - - integer :: i - integer(HID_T) :: fileHandle, groupHandle - character(len=pStringLen) :: fileName, datasetName - - print*, ' writing field and constitutive data required for restart to file';flush(IO_STDOUT) - - write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5' - fileHandle = HDF5_openFile(fileName,'a') - - call HDF5_write(fileHandle,crystallite_partitionedF,'F') - call HDF5_write(fileHandle,crystallite_Fp, 'F_p') - call HDF5_write(fileHandle,crystallite_Fi, 'F_i') - call HDF5_write(fileHandle,crystallite_Lp, 'L_p') - call HDF5_write(fileHandle,crystallite_Li, 'L_i') - call HDF5_write(fileHandle,crystallite_S, 'S') - - groupHandle = HDF5_addGroup(fileHandle,'phase') - do i = 1,size(material_name_phase) - write(datasetName,'(i0,a)') i,'_omega' - call HDF5_write(groupHandle,plasticState(i)%state,datasetName) - enddo - call HDF5_closeGroup(groupHandle) - - groupHandle = HDF5_addGroup(fileHandle,'homogenization') - do i = 1, size(material_name_homogenization) - write(datasetName,'(i0,a)') i,'_omega' - call HDF5_write(groupHandle,homogState(i)%state,datasetName) - enddo - call HDF5_closeGroup(groupHandle) - - call HDF5_closeFile(fileHandle) - -end subroutine crystallite_restartWrite - - -!-------------------------------------------------------------------------------------------------- -!> @brief Read data for restart -! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively -!-------------------------------------------------------------------------------------------------- -subroutine crystallite_restartRead - - integer :: i - integer(HID_T) :: fileHandle, groupHandle - character(len=pStringLen) :: fileName, datasetName - - print'(/,a,i0,a)', ' reading restart information of increment from file' - - write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5' - fileHandle = HDF5_openFile(fileName) - - call HDF5_read(fileHandle,crystallite_F0, 'F') - call HDF5_read(fileHandle,crystallite_Fp0,'F_p') - call HDF5_read(fileHandle,crystallite_Fi0,'F_i') - call HDF5_read(fileHandle,crystallite_Lp0,'L_p') - call HDF5_read(fileHandle,crystallite_Li0,'L_i') - call HDF5_read(fileHandle,crystallite_S0, 'S') - - groupHandle = HDF5_openGroup(fileHandle,'phase') - do i = 1,size(material_name_phase) - write(datasetName,'(i0,a)') i,'_omega' - call HDF5_read(groupHandle,plasticState(i)%state0,datasetName) - enddo - call HDF5_closeGroup(groupHandle) - - groupHandle = HDF5_openGroup(fileHandle,'homogenization') - do i = 1,size(material_name_homogenization) - write(datasetName,'(i0,a)') i,'_omega' - call HDF5_read(groupHandle,homogState(i)%state0,datasetName) - enddo - call HDF5_closeGroup(groupHandle) - - call HDF5_closeFile(fileHandle) - -end subroutine crystallite_restartRead - - -!-------------------------------------------------------------------------------------------------- -!> @brief Forward data after successful increment. -! ToDo: Any guessing for the current states possible? -!-------------------------------------------------------------------------------------------------- -subroutine crystallite_forward - - integer :: i, j - - crystallite_F0 = crystallite_partitionedF - crystallite_Fp0 = crystallite_Fp - crystallite_Lp0 = crystallite_Lp - crystallite_Fi0 = crystallite_Fi - crystallite_Li0 = crystallite_Li - crystallite_S0 = crystallite_S - - do i = 1, size(plasticState) - plasticState(i)%state0 = plasticState(i)%state - enddo - do i = 1, size(sourceState) - do j = 1,phase_Nsources(i) - sourceState(i)%p(j)%state0 = sourceState(i)%p(j)%state - enddo; enddo - do i = 1,size(material_name_homogenization) - homogState (i)%state0 = homogState (i)%state - thermalState(i)%state0 = thermalState(i)%state - damageState (i)%state0 = damageState (i)%state - enddo - -end subroutine crystallite_forward - -end module crystallite diff --git a/src/damage_local.f90 b/src/damage_local.f90 deleted file mode 100644 index 97eaf9a8c..000000000 --- a/src/damage_local.f90 +++ /dev/null @@ -1,172 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for locally evolving damage field -!-------------------------------------------------------------------------------------------------- -module damage_local - use prec - use IO - use material - use config - use YAML_types - use constitutive - use results - - implicit none - private - - type :: tParameters - character(len=pStringLen), allocatable, dimension(:) :: & - output - end type tParameters - - type, private :: tNumerics - real(pReal) :: & - residualStiffness !< non-zero residual damage - end type tNumerics - - type(tparameters), dimension(:), allocatable :: & - param - - type(tNumerics), private :: num - - public :: & - damage_local_init, & - damage_local_updateState, & - damage_local_results - -contains - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine damage_local_init - - integer :: Ninstances,Nmaterialpoints,h - class(tNode), pointer :: & - num_generic, & - material_homogenization, & - homog, & - homogDamage - - print'(/,a)', ' <<<+- damage_local init -+>>>'; flush(IO_STDOUT) - -!---------------------------------------------------------------------------------------------- -! read numerics parameter and do sanity check - num_generic => config_numerics%get('generic',defaultVal=emptyDict) - num%residualStiffness = num_generic%get_asFloat('residualStiffness', defaultVal=1.0e-6_pReal) - if (num%residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness') - - Ninstances = count(damage_type == DAMAGE_local_ID) - allocate(param(Ninstances)) - - material_homogenization => config_material%get('homogenization') - do h = 1, material_homogenization%length - if (damage_type(h) /= DAMAGE_LOCAL_ID) cycle - homog => material_homogenization%get(h) - homogDamage => homog%get('damage') - associate(prm => param(damage_typeInstance(h))) - -#if defined (__GFORTRAN__) - prm%output = output_asStrings(homogDamage) -#else - prm%output = homogDamage%get_asStrings('output',defaultVal=emptyStringArray) -#endif - - Nmaterialpoints = count(material_homogenizationAt == h) - damageState(h)%sizeState = 1 - allocate(damageState(h)%state0 (1,Nmaterialpoints), source=1.0_pReal) - allocate(damageState(h)%subState0(1,Nmaterialpoints), source=1.0_pReal) - allocate(damageState(h)%state (1,Nmaterialpoints), source=1.0_pReal) - - damage(h)%p => damageState(h)%state(1,:) - - end associate - enddo - -end subroutine damage_local_init - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates local change in damage field -!-------------------------------------------------------------------------------------------------- -function damage_local_updateState(subdt, ip, el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - subdt - logical, dimension(2) :: & - damage_local_updateState - integer :: & - homog, & - offset - real(pReal) :: & - phi, phiDot, dPhiDot_dPhi - - homog = material_homogenizationAt(el) - offset = material_homogenizationMemberAt(ip,el) - phi = damageState(homog)%subState0(1,offset) - call damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) - phi = max(num%residualStiffness,min(1.0_pReal,phi + subdt*phiDot)) - - damage_local_updateState = [ abs(phi - damageState(homog)%state(1,offset)) & - <= 1.0e-2_pReal & - .or. abs(phi - damageState(homog)%state(1,offset)) & - <= 1.0e-6_pReal*abs(damageState(homog)%state(1,offset)), & - .true.] - - damageState(homog)%state(1,offset) = phi - -end function damage_local_updateState - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates homogenized local damage driving forces -!-------------------------------------------------------------------------------------------------- -subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - phi - real(pReal) :: & - phiDot, dPhiDot_dPhi - - phiDot = 0.0_pReal - dPhiDot_dPhi = 0.0_pReal - - call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) - - phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) - dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) - -end subroutine damage_local_getSourceAndItsTangent - - -!-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file -!-------------------------------------------------------------------------------------------------- -subroutine damage_local_results(homog,group) - - integer, intent(in) :: homog - character(len=*), intent(in) :: group - - integer :: o - - associate(prm => param(damage_typeInstance(homog))) - outputsLoop: do o = 1,size(prm%output) - select case(prm%output(o)) - case ('phi') - call results_writeDataset(group,damage(homog)%p,prm%output(o),& - 'damage indicator','-') - end select - enddo outputsLoop - end associate - -end subroutine damage_local_results - - -end module damage_local diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index c4426f185..ac4d8636a 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -7,7 +7,6 @@ module damage_nonlocal use material use config use YAML_types - use crystallite use lattice use constitutive use results diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 57478e039..cbab8e468 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -11,14 +11,11 @@ module homogenization use math use material use constitutive - use crystallite use FEsolving use discretization use thermal_isothermal - use thermal_adiabatic use thermal_conduction use damage_none - use damage_local use damage_nonlocal use results @@ -162,11 +159,9 @@ subroutine homogenization_init call mech_init(num_homog) if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init - if (any(thermal_type == THERMAL_adiabatic_ID)) call thermal_adiabatic_init if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init if (any(damage_type == DAMAGE_none_ID)) call damage_none_init - if (any(damage_type == DAMAGE_local_ID)) call damage_local_init if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init @@ -212,10 +207,6 @@ subroutine materialpoint_stressAndItsTangent(dt) homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & homogState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e)) - if (thermalState(material_homogenizationAt(e))%sizeState > 0) & - thermalState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & - thermalState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e)) - if (damageState(material_homogenizationAt(e))%sizeState > 0) & damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & damageState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e)) @@ -245,9 +236,6 @@ subroutine materialpoint_stressAndItsTangent(dt) if(homogState(material_homogenizationAt(e))%sizeState > 0) & homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & homogState(material_homogenizationAt(e))%State (:,material_homogenizationMemberAt(i,e)) - if(thermalState(material_homogenizationAt(e))%sizeState > 0) & - thermalState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & - thermalState(material_homogenizationAt(e))%State (:,material_homogenizationMemberAt(i,e)) if(damageState(material_homogenizationAt(e))%sizeState > 0) & damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & damageState(material_homogenizationAt(e))%State (:,material_homogenizationMemberAt(i,e)) @@ -266,13 +254,11 @@ subroutine materialpoint_stressAndItsTangent(dt) subStep(i,e) = num%subStepSizeHomog * subStep(i,e) ! crystallite had severe trouble, so do a significant cutback call crystallite_restore(i,e,subStep(i,e) < 1.0_pReal) + call constitutive_restore(i,e) if(homogState(material_homogenizationAt(e))%sizeState > 0) & homogState(material_homogenizationAt(e))%State( :,material_homogenizationMemberAt(i,e)) = & homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) - if(thermalState(material_homogenizationAt(e))%sizeState > 0) & - thermalState(material_homogenizationAt(e))%State( :,material_homogenizationMemberAt(i,e)) = & - thermalState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) if(damageState(material_homogenizationAt(e))%sizeState > 0) & damageState(material_homogenizationAt(e))%State( :,material_homogenizationMemberAt(i,e)) = & damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) @@ -400,24 +386,6 @@ function updateState(subdt,subF,ip,el) el) end select chosenHomogenization - chosenThermal: select case (thermal_type(material_homogenizationAt(el))) - case (THERMAL_adiabatic_ID) chosenThermal - updateState = & - updateState .and. & - thermal_adiabatic_updateState(subdt, & - ip, & - el) - end select chosenThermal - - chosenDamage: select case (damage_type(material_homogenizationAt(el))) - case (DAMAGE_local_ID) chosenDamage - updateState = & - updateState .and. & - damage_local_updateState(subdt, & - ip, & - el) - end select chosenDamage - end function updateState @@ -441,8 +409,6 @@ subroutine homogenization_results group = trim(group_base)//'/damage' call results_closeGroup(results_addGroup(group)) select case(damage_type(p)) - case(DAMAGE_LOCAL_ID) - call damage_local_results(p,group) case(DAMAGE_NONLOCAL_ID) call damage_nonlocal_results(p,group) end select @@ -450,8 +416,6 @@ subroutine homogenization_results group = trim(group_base)//'/thermal' call results_closeGroup(results_addGroup(group)) select case(thermal_type(p)) - case(THERMAL_ADIABATIC_ID) - call thermal_adiabatic_results(p,group) case(THERMAL_CONDUCTION_ID) call thermal_conduction_results(p,group) end select diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index 36a882a48..5265d6172 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -26,7 +26,7 @@ contains !-------------------------------------------------------------------------------------------------- module function kinematics_thermal_expansion_init(kinematics_length) result(myKinematics) - integer, intent(in) :: kinematics_length + integer, intent(in) :: kinematics_length logical, dimension(:,:), allocatable :: myKinematics integer :: Ninstances,p,i,k @@ -35,8 +35,8 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi phases, & phase, & kinematics, & - kinematic_type - + kinematic_type + print'(/,a)', ' <<<+- kinematics_thermal_expansion init -+>>>' myKinematics = kinematics_active('thermal_expansion',kinematics_length) @@ -50,13 +50,13 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi do p = 1, phases%length if(any(myKinematics(:,p))) kinematics_thermal_expansion_instance(p) = count(myKinematics(:,1:p)) - phase => phases%get(p) + phase => phases%get(p) if(count(myKinematics(:,p)) == 0) cycle kinematics => phase%get('kinematics') do k = 1, kinematics%length if(myKinematics(k,p)) then associate(prm => param(kinematics_thermal_expansion_instance(p))) - kinematic_type => kinematics%get(k) + kinematic_type => kinematics%get(k) prm%T_ref = kinematic_type%get_asFloat('T_ref', defaultVal=0.0_pReal) @@ -81,29 +81,6 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi end function kinematics_thermal_expansion_init -!-------------------------------------------------------------------------------------------------- -!> @brief report initial thermal strain based on current temperature deviation from reference -!-------------------------------------------------------------------------------------------------- -pure module function kinematics_thermal_expansion_initialStrain(homog,phase,offset) result(initialStrain) - - integer, intent(in) :: & - phase, & - homog, & - offset - - real(pReal), dimension(3,3) :: & - initialStrain !< initial thermal strain (should be small strain, though) - - associate(prm => param(kinematics_thermal_expansion_instance(phase))) - initialStrain = & - (temperature(homog)%p(offset) - prm%T_ref)**1 / 1. * prm%A(1:3,1:3,1) + & ! constant coefficient - (temperature(homog)%p(offset) - prm%T_ref)**2 / 2. * prm%A(1:3,1:3,2) + & ! linear coefficient - (temperature(homog)%p(offset) - prm%T_ref)**3 / 3. * prm%A(1:3,1:3,3) ! quadratic coefficient - end associate - -end function kinematics_thermal_expansion_initialStrain - - !-------------------------------------------------------------------------------------------------- !> @brief constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- diff --git a/src/material.f90 b/src/material.f90 index bb5f484f6..1f2437ad3 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -41,10 +41,8 @@ module material STIFFNESS_DEGRADATION_UNDEFINED_ID, & STIFFNESS_DEGRADATION_DAMAGE_ID, & THERMAL_ISOTHERMAL_ID, & - THERMAL_ADIABATIC_ID, & THERMAL_CONDUCTION_ID, & DAMAGE_NONE_ID, & - DAMAGE_LOCAL_ID, & DAMAGE_NONLOCAL_ID, & HOMOGENIZATION_UNDEFINED_ID, & HOMOGENIZATION_NONE_ID, & @@ -86,7 +84,6 @@ module material type(tState), allocatable, dimension(:), public :: & homogState, & - thermalState, & damageState type(Rotation), dimension(:,:,:), allocatable, public, protected :: & @@ -123,10 +120,8 @@ module material STIFFNESS_DEGRADATION_UNDEFINED_ID, & STIFFNESS_DEGRADATION_DAMAGE_ID, & THERMAL_ISOTHERMAL_ID, & - THERMAL_ADIABATIC_ID, & THERMAL_CONDUCTION_ID, & DAMAGE_NONE_ID, & - DAMAGE_LOCAL_ID, & DAMAGE_NONLOCAL_ID, & HOMOGENIZATION_NONE_ID, & HOMOGENIZATION_ISOSTRAIN_ID, & @@ -152,7 +147,6 @@ subroutine material_init(restart) allocate(homogState (size(material_name_homogenization))) - allocate(thermalState (size(material_name_homogenization))) allocate(damageState (size(material_name_homogenization))) allocate(temperature (size(material_name_homogenization))) @@ -218,8 +212,6 @@ subroutine material_parseHomogenization select case (homogThermal%get_asString('type')) case('isothermal') thermal_type(h) = THERMAL_isothermal_ID - case('adiabatic') - thermal_type(h) = THERMAL_adiabatic_ID case('conduction') thermal_type(h) = THERMAL_conduction_ID case default @@ -232,8 +224,6 @@ subroutine material_parseHomogenization select case (homogDamage%get_asString('type')) case('none') damage_type(h) = DAMAGE_none_ID - case('local') - damage_type(h) = DAMAGE_local_ID case('nonlocal') damage_type(h) = DAMAGE_nonlocal_ID case default diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 deleted file mode 100644 index c67d004bf..000000000 --- a/src/thermal_adiabatic.f90 +++ /dev/null @@ -1,226 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for adiabatic temperature evolution -!-------------------------------------------------------------------------------------------------- -module thermal_adiabatic - use prec - use config - use material - use results - use constitutive - use YAML_types - use crystallite - use lattice - - implicit none - private - - type :: tParameters - character(len=pStringLen), allocatable, dimension(:) :: & - output - end type tParameters - - type(tparameters), dimension(:), allocatable :: & - param - - public :: & - thermal_adiabatic_init, & - thermal_adiabatic_updateState, & - thermal_adiabatic_getSourceAndItsTangent, & - thermal_adiabatic_getSpecificHeat, & - thermal_adiabatic_getMassDensity, & - thermal_adiabatic_results - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine thermal_adiabatic_init - - integer :: maxNinstances,h,Nmaterialpoints - class(tNode), pointer :: & - material_homogenization, & - homog, & - homogThermal - - print'(/,a)', ' <<<+- thermal_adiabatic init -+>>>'; flush(6) - - maxNinstances = count(thermal_type == THERMAL_adiabatic_ID) - if (maxNinstances == 0) return - - allocate(param(maxNinstances)) - - material_homogenization => config_material%get('homogenization') - do h = 1, size(material_name_homogenization) - if (thermal_type(h) /= THERMAL_adiabatic_ID) cycle - homog => material_homogenization%get(h) - homogThermal => homog%get('thermal') - - associate(prm => param(thermal_typeInstance(h))) - -#if defined (__GFORTRAN__) - prm%output = output_asStrings(homogThermal) -#else - prm%output = homogThermal%get_asStrings('output',defaultVal=emptyStringArray) -#endif - - Nmaterialpoints=count(material_homogenizationAt==h) - thermalState(h)%sizeState = 1 - allocate(thermalState(h)%state0 (1,Nmaterialpoints), source=thermal_initialT(h)) - allocate(thermalState(h)%subState0(1,Nmaterialpoints), source=thermal_initialT(h)) - allocate(thermalState(h)%state (1,Nmaterialpoints), source=thermal_initialT(h)) - - temperature(h)%p => thermalState(h)%state(1,:) - allocate(temperatureRate(h)%p(Nmaterialpoints),source = 0.0_pReal) - - end associate - enddo - -end subroutine thermal_adiabatic_init - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates adiabatic change in temperature based on local heat generation model -!-------------------------------------------------------------------------------------------------- -function thermal_adiabatic_updateState(subdt, ip, el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - subdt - - logical, dimension(2) :: & - thermal_adiabatic_updateState - integer :: & - homog, & - offset - real(pReal) :: & - T, Tdot, dTdot_dT - - homog = material_homogenizationAt(el) - offset = material_homogenizationMemberAt(ip,el) - - T = thermalState(homog)%subState0(1,offset) - call thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) - T = T + subdt*Tdot/(thermal_adiabatic_getSpecificHeat(ip,el)*thermal_adiabatic_getMassDensity(ip,el)) - - thermal_adiabatic_updateState = [ abs(T - thermalState(homog)%state(1,offset)) & - <= 1.0e-2_pReal & - .or. abs(T - thermalState(homog)%state(1,offset)) & - <= 1.0e-6_pReal*abs(thermalState(homog)%state(1,offset)), & - .true.] - - temperature (homog)%p(material_homogenizationMemberAt(ip,el)) = T - temperatureRate(homog)%p(material_homogenizationMemberAt(ip,el)) = & - (thermalState(homog)%state(1,offset) - thermalState(homog)%subState0(1,offset))/(subdt+tiny(0.0_pReal)) - -end function thermal_adiabatic_updateState - - -!-------------------------------------------------------------------------------------------------- -!> @brief returns heat generation rate -!-------------------------------------------------------------------------------------------------- -subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - T - real(pReal), intent(out) :: & - Tdot, dTdot_dT - integer :: & - homog - - Tdot = 0.0_pReal - dTdot_dT = 0.0_pReal - - homog = material_homogenizationAt(el) - call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S, crystallite_Lp, ip, el) - - Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal) - dTdot_dT = dTdot_dT/real(homogenization_Nconstituents(homog),pReal) - -end subroutine thermal_adiabatic_getSourceAndItsTangent - - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized specific heat capacity -!-------------------------------------------------------------------------------------------------- -function thermal_adiabatic_getSpecificHeat(ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - - real(pReal) :: & - thermal_adiabatic_getSpecificHeat - integer :: & - grain - - thermal_adiabatic_getSpecificHeat = 0.0_pReal - - do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat & - + lattice_c_p(material_phaseAt(grain,el)) - enddo - - thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat & - / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) - -end function thermal_adiabatic_getSpecificHeat - - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized mass density -!-------------------------------------------------------------------------------------------------- -function thermal_adiabatic_getMassDensity(ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal) :: & - thermal_adiabatic_getMassDensity - integer :: & - grain - - thermal_adiabatic_getMassDensity = 0.0_pReal - - do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity & - + lattice_rho(material_phaseAt(grain,el)) - enddo - - thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity & - / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) - -end function thermal_adiabatic_getMassDensity - - -!-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file -!-------------------------------------------------------------------------------------------------- -subroutine thermal_adiabatic_results(homog,group) - - integer, intent(in) :: homog - character(len=*), intent(in) :: group - - integer :: o - - associate(prm => param(damage_typeInstance(homog))) - outputsLoop: do o = 1,size(prm%output) - select case(trim(prm%output(o))) - case('T') - call results_writeDataset(group,temperature(homog)%p,'T',& - 'temperature','K') - end select - enddo outputsLoop - end associate - -end subroutine thermal_adiabatic_results - -end module thermal_adiabatic diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 602bdab35..d30e50677 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -8,7 +8,6 @@ module thermal_conduction use config use lattice use results - use crystallite use constitutive use YAML_types @@ -66,10 +65,6 @@ subroutine thermal_conduction_init #endif Nmaterialpoints=count(material_homogenizationAt==h) - thermalState(h)%sizeState = 0 - allocate(thermalState(h)%state0 (0,Nmaterialpoints)) - allocate(thermalState(h)%subState0(0,Nmaterialpoints)) - allocate(thermalState(h)%state (0,Nmaterialpoints)) allocate (temperature (h)%p(Nmaterialpoints), source=thermal_initialT(h)) allocate (temperatureRate(h)%p(Nmaterialpoints), source=0.0_pReal) diff --git a/src/thermal_isothermal.f90 b/src/thermal_isothermal.f90 index adf2257de..2a41ada49 100644 --- a/src/thermal_isothermal.f90 +++ b/src/thermal_isothermal.f90 @@ -25,10 +25,6 @@ subroutine thermal_isothermal_init if (thermal_type(h) /= THERMAL_isothermal_ID) cycle Nmaterialpoints = count(material_homogenizationAt == h) - thermalState(h)%sizeState = 0 - allocate(thermalState(h)%state0 (0,Nmaterialpoints)) - allocate(thermalState(h)%subState0(0,Nmaterialpoints)) - allocate(thermalState(h)%state (0,Nmaterialpoints)) allocate(temperature (h)%p(Nmaterialpoints),source=thermal_initialT(h)) allocate(temperatureRate(h)%p(Nmaterialpoints),source = 0.0_pReal)