Merge remote-tracking branch 'remotes/origin/improve-Lp-guessing' into development

This commit is contained in:
Franz Roters 2019-05-15 08:10:46 +02:00
commit d29967d8b2
4 changed files with 109 additions and 171 deletions

@ -1 +1 @@
Subproject commit aadf2d82a7e04646e3f20c3d526f27a6f90acef0
Subproject commit 183d7a3a3bafa0a308c3eac858ca03c08fc03d50

View File

@ -21,8 +21,7 @@ module crystallite
homogenization_Ngrains
use future
implicit none
implicit none
private
character(len=64), dimension(:,:), allocatable, private :: &
crystallite_output !< name of each post result output
@ -111,6 +110,24 @@ module crystallite
end type tOutput
type(tOutput), allocatable, dimension(:), private :: output_constituent
type, private :: 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?
procedure(), pointer :: integrateState
public :: &
@ -165,14 +182,13 @@ subroutine crystallite_init
use config, only: &
config_deallocate, &
config_crystallite, &
config_numerics, &
config_phase, &
crystallite_name
use constitutive, only: &
constitutive_initialFi, &
constitutive_microstructure ! derived (shortcut) quantities of given state
implicit none
integer, parameter :: FILEUNIT=434
logical, dimension(:,:), allocatable :: devNull
integer :: &
@ -242,6 +258,38 @@ subroutine crystallite_init
allocate(crystallite_sizePostResults(size(config_crystallite)),source=0)
allocate(crystallite_sizePostResult(maxval(crystallite_Noutput), &
size(config_crystallite)), source=0)
num%subStepMinCryst = config_numerics%getFloat('substepmincryst', defaultVal=1.0e-3_pReal)
num%subStepSizeCryst = config_numerics%getFloat('substepsizecryst', defaultVal=0.25_pReal)
num%stepIncreaseCryst = config_numerics%getFloat('stepincreasecryst', defaultVal=1.5_pReal)
num%subStepSizeLp = config_numerics%getFloat('substepsizelp', defaultVal=0.5_pReal)
num%subStepSizeLi = config_numerics%getFloat('substepsizeli', defaultVal=0.5_pReal)
num%rTol_crystalliteState = config_numerics%getFloat('rtol_crystallitestate', defaultVal=1.0e-6_pReal)
num%rTol_crystalliteStress = config_numerics%getFloat('rtol_crystallitestress',defaultVal=1.0e-6_pReal)
num%aTol_crystalliteStress = config_numerics%getFloat('atol_crystallitestress',defaultVal=1.0e-8_pReal)
num%iJacoLpresiduum = config_numerics%getInt ('ijacolpresiduum', defaultVal=1)
num%nState = config_numerics%getInt ('nstate', defaultVal=20)
num%nStress = config_numerics%getInt ('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(numerics_integrator)
case(1)
@ -433,10 +481,6 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
use prec, only: &
tol_math_check, &
dNeq0
use numerics, only: &
subStepMinCryst, &
subStepSizeCryst, &
stepIncreaseCryst
#ifdef DEBUG
use debug, only: &
debug_level, &
@ -463,7 +507,6 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
phase_Nsources, &
phaseAt, phasememberAt
implicit none
logical, dimension(theMesh%elem%nIPs,theMesh%Nelems) :: crystallite_stress
real(pReal), intent(in), optional :: &
dummyArgumentToPreventInternalCompilerErrorWithGCC
@ -519,7 +562,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partionedF0(1:3,1:3,c,i,e)
crystallite_subS0(1:3,1:3,c,i,e) = crystallite_partionedS0(1:3,1:3,c,i,e)
crystallite_subFrac(c,i,e) = 0.0_pReal
crystallite_subStep(c,i,e) = 1.0_pReal/subStepSizeCryst
crystallite_subStep(c,i,e) = 1.0_pReal/num%subStepSizeCryst
crystallite_todo(c,i,e) = .true.
crystallite_converged(c,i,e) = .false. ! pretend failed step of 1/subStepSizeCryst
endif homogenizationRequestsCalculation
@ -554,7 +597,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
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), &
stepIncreaseCryst * crystallite_subStep(c,i,e))
num%stepIncreaseCryst * crystallite_subStep(c,i,e))
crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on?
if (crystallite_todo(c,i,e)) then
@ -584,7 +627,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
!--------------------------------------------------------------------------------------------------
! cut back (reduced time and restore)
else
crystallite_subStep(c,i,e) = subStepSizeCryst * crystallite_subStep(c,i,e)
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_invFp(1:3,1:3,c,i,e) = math_inv33(crystallite_Fp (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)
@ -602,7 +645,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
enddo
! cant restore dotState here, since not yet calculated in first cutback after initialization
crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > subStepMinCryst ! still on track or already done (beyond repair)
crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > num%subStepMinCryst ! still on track or already done (beyond repair)
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 &
.and. ((e == debug_e .and. i == debug_i .and. c == debug_g) &
@ -652,7 +695,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
!--------------------------------------------------------------------------------------------------
! integrate --- requires fully defined state array (basic + dependent state)
if (any(crystallite_todo)) call integrateState ! TODO: unroll into proper elementloop to avoid N^2 for single point evaluation
where(.not. crystallite_converged .and. crystallite_subStep > subStepMinCryst) & ! do not try non-converged but fully cutbacked any further
where(.not. crystallite_converged .and. crystallite_subStep > num%subStepMinCryst) & ! do not try non-converged but fully cutbacked any further
crystallite_todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation
@ -702,7 +745,7 @@ end function crystallite_stress
!--------------------------------------------------------------------------------------------------
!> @brief calculate tangent (dPdF)
!--------------------------------------------------------------------------------------------------
subroutine crystallite_stressTangent()
subroutine crystallite_stressTangent
use prec, only: &
tol_math_check, &
dNeq0
@ -728,7 +771,6 @@ subroutine crystallite_stressTangent()
constitutive_LpAndItsTangents, &
constitutive_LiAndItsTangents
implicit none
integer :: &
c, & !< counter in integration point component loop
i, & !< counter in integration point loop
@ -882,7 +924,6 @@ subroutine crystallite_orientations
use plastic_nonlocal, only: &
plastic_nonlocal_updateCompatibility
implicit none
integer &
c, & !< counter in integration point component loop
i, & !< counter in integration point loop
@ -919,7 +960,6 @@ function crystallite_push33ToRef(ipc,ip,el, tensor33)
use material, only: &
material_EulerAngles ! ToDo: Why stored? We also have crystallite_orientation0
implicit none
real(pReal), dimension(3,3) :: crystallite_push33ToRef
real(pReal), dimension(3,3), intent(in) :: tensor33
real(pReal), dimension(3,3) :: T
@ -962,7 +1002,6 @@ function crystallite_postResults(ipc, ip, el)
use rotations, only: &
rotation
implicit none
integer, intent(in):: &
el, & !< element index
ip, & !< integration point index
@ -1089,7 +1128,6 @@ subroutine crystallite_results
use material, only: &
material_phase_plasticity_type => phase_plasticity
implicit none
integer :: p,o
real(pReal), allocatable, dimension(:,:,:) :: selected_tensors
type(rotation), allocatable, dimension(:) :: selected_rotations
@ -1233,12 +1271,6 @@ logical function integrateStress(ipc,ip,el,timeFraction)
IEEE_arithmetic
use prec, only: tol_math_check, &
dEq0
use numerics, only: nStress, &
aTol_crystalliteStress, &
rTol_crystalliteStress, &
iJacoLpresiduum, &
subStepSizeLp, &
subStepSizeLi
#ifdef DEBUG
use debug, only: debug_level, &
debug_e, &
@ -1263,7 +1295,6 @@ logical function integrateStress(ipc,ip,el,timeFraction)
math_33to9, &
math_9to33
implicit none
integer, intent(in):: el, & ! element index
ip, & ! integration point index
ipc ! grain index
@ -1382,10 +1413,10 @@ logical function integrateStress(ipc,ip,el,timeFraction)
LiLoop: do
NiterationStressLi = NiterationStressLi + 1
LiLoopLimit: if (NiterationStressLi > nStress) then
LiLoopLimit: if (NiterationStressLi > num%nStress) then
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) &
write(6,'(a,i3,a,i8,1x,i2,1x,i3,/)') '<< CRYST integrateStress >> reached Li loop limit',nStress, &
write(6,'(a,i3,a,i8,1x,i2,1x,i3,/)') '<< CRYST integrateStress >> reached Li loop limit',num%nStress, &
' at el ip ipc ', el,ip,ipc
#endif
return
@ -1404,10 +1435,10 @@ logical function integrateStress(ipc,ip,el,timeFraction)
LpLoop: do
NiterationStressLp = NiterationStressLp + 1
LpLoopLimit: if (NiterationStressLp > nStress) then
LpLoopLimit: if (NiterationStressLp > num%nStress) then
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) &
write(6,'(a,i3,a,i8,1x,i2,1x,i3,/)') '<< CRYST integrateStress >> reached Lp loop limit',nStress, &
write(6,'(a,i3,a,i8,1x,i2,1x,i3,/)') '<< CRYST integrateStress >> reached Lp loop limit',num%nStress, &
' at el ip ipc ', el,ip,ipc
#endif
return
@ -1438,8 +1469,8 @@ logical function integrateStress(ipc,ip,el,timeFraction)
#endif
!* update current residuum and check for convergence of loop
aTolLp = max(rTol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error
aTol_crystalliteStress) ! minimum lower cutoff
aTolLp = 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
@ -1459,7 +1490,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
Lpguess_old = Lpguess
steplengthLp = 1.0_pReal ! ...proceed with normal step length (calculate new search direction)
else ! not converged and residuum not improved...
steplengthLp = subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction
steplengthLp = num%subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction
Lpguess = Lpguess_old + steplengthLp * deltaLp
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 &
@ -1473,7 +1504,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
!* calculate Jacobian for correction term
if (mod(jacoCounterLp, iJacoLpresiduum) == 0) then
if (mod(jacoCounterLp, num%iJacoLpresiduum) == 0) then
do o=1,3; do p=1,3
dFe_dLp(o,1:3,p,1:3) = A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
enddo; enddo
@ -1537,8 +1568,8 @@ logical function integrateStress(ipc,ip,el,timeFraction)
#endif
!* update current residuum and check for convergence of loop
aTolLi = max(rTol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error
aTol_crystalliteStress) ! minimum lower cutoff
aTolLi = 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 ! NaN in residuum...
#ifdef DEBUG
@ -1557,13 +1588,13 @@ logical function integrateStress(ipc,ip,el,timeFraction)
Liguess_old = Liguess
steplengthLi = 1.0_pReal ! ...proceed with normal step length (calculate new search direction)
else ! not converged and residuum not improved...
steplengthLi = subStepSizeLi * steplengthLi ! ...try with smaller step length in same direction
steplengthLi = num%subStepSizeLi * steplengthLi ! ...try with smaller step length in same direction
Liguess = Liguess_old + steplengthLi * deltaLi
cycle LiLoop
endif
!* calculate Jacobian for correction term
if (mod(jacoCounterLi, iJacoLpresiduum) == 0) then
if (mod(jacoCounterLi, num%iJacoLpresiduum) == 0) then
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)
@ -1661,7 +1692,7 @@ end function integrateStress
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
subroutine integrateStateFPI()
subroutine integrateStateFPI
#ifdef DEBUG
use debug, only: debug_level, &
debug_e, &
@ -1672,8 +1703,6 @@ subroutine integrateStateFPI()
debug_levelExtensive, &
debug_levelSelective
#endif
use numerics, only: &
nState
use mesh, only: &
mesh_element
use material, only: &
@ -1686,8 +1715,6 @@ subroutine integrateStateFPI()
constitutive_plasticity_maxSizeDotState, &
constitutive_source_maxSizeDotState
implicit none
integer :: &
NiterationState, & !< number of iterations in state loop
e, & !< element index in element loop
@ -1712,7 +1739,7 @@ subroutine integrateStateFPI()
NiterationState = 0
doneWithIntegration = .false.
crystalliteLooping: do while (.not. doneWithIntegration .and. NiterationState < nState)
crystalliteLooping: do while (.not. doneWithIntegration .and. NiterationState < num%nState)
NiterationState = NiterationState + 1
#ifdef DEBUG
@ -1852,7 +1879,6 @@ subroutine integrateStateFPI()
!--------------------------------------------------------------------------------------------------
real(pReal) pure function damper(current,previous,previous2)
implicit none
real(pReal), dimension(:), intent(in) ::&
current, previous, previous2
@ -1874,12 +1900,10 @@ end subroutine integrateStateFPI
!--------------------------------------------------------------------------------------------------
!> @brief integrate state with 1st order explicit Euler method
!--------------------------------------------------------------------------------------------------
subroutine integrateStateEuler()
subroutine integrateStateEuler
use material, only: &
plasticState
implicit none
call update_dotState(1.0_pReal)
call update_state(1.0_pReal)
call update_deltaState
@ -1894,7 +1918,7 @@ end subroutine integrateStateEuler
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with 1st order Euler method with adaptive step size
!--------------------------------------------------------------------------------------------------
subroutine integrateStateAdaptiveEuler()
subroutine integrateStateAdaptiveEuler
use mesh, only: &
theMesh, &
mesh_element
@ -1909,7 +1933,6 @@ subroutine integrateStateAdaptiveEuler()
constitutive_plasticity_maxSizeDotState, &
constitutive_source_maxSizeDotState
implicit none
integer :: &
e, & ! element index in element loop
i, & ! integration point index in ip loop
@ -2001,7 +2024,7 @@ end subroutine integrateStateAdaptiveEuler
!> @brief integrate stress, state with 4th order explicit Runge Kutta method
! ToDo: This is totally BROKEN: RK4dotState is never used!!!
!--------------------------------------------------------------------------------------------------
subroutine integrateStateRK4()
subroutine integrateStateRK4
use mesh, only: &
mesh_element
use material, only: &
@ -2011,7 +2034,6 @@ subroutine integrateStateRK4()
phase_Nsources, &
phaseAt, phasememberAt
implicit none
real(pReal), dimension(4), parameter :: &
TIMESTEPFRACTION = [0.5_pReal, 0.5_pReal, 1.0_pReal, 1.0_pReal] ! factor giving the fraction of the original timestep used for Runge Kutta Integration
real(pReal), dimension(4), parameter :: &
@ -2069,7 +2091,7 @@ end subroutine integrateStateRK4
!> @brief integrate stress, state with 5th order Runge-Kutta Cash-Karp method with
!> adaptive step size (use 5th order solution to advance = "local extrapolation")
!--------------------------------------------------------------------------------------------------
subroutine integrateStateRKCK45()
subroutine integrateStateRKCK45
use mesh, only: &
mesh_element, &
theMesh
@ -2084,7 +2106,6 @@ subroutine integrateStateRKCK45()
constitutive_plasticity_maxSizeDotState, &
constitutive_source_maxSizeDotState
implicit none
real(pReal), dimension(5,5), parameter :: &
A = reshape([&
.2_pReal, .075_pReal, .3_pReal, -11.0_pReal/54.0_pReal, 1631.0_pReal/55296.0_pReal, &
@ -2252,9 +2273,7 @@ end subroutine integrateStateRKCK45
!> @brief sets convergence flag for nonlocal calculations
!> @detail one non-converged nonlocal sets all other nonlocals to non-converged to trigger cut back
!--------------------------------------------------------------------------------------------------
subroutine nonlocalConvergenceCheck()
implicit none
subroutine nonlocalConvergenceCheck
if (any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) & ! any non-local not yet converged (or broken)...
where( .not. crystallite_localPlasticity) crystallite_converged = .false.
@ -2267,10 +2286,10 @@ end subroutine nonlocalConvergenceCheck
! still .true. is considered as converged
!> @details: For explicitEuler, RK4 and RKCK45, adaptive Euler and FPI have their on criteria
!--------------------------------------------------------------------------------------------------
subroutine setConvergenceFlag()
subroutine setConvergenceFlag
use mesh, only: &
mesh_element
implicit none
integer :: &
e, & !< element index in element loop
i, & !< integration point index in ip loop
@ -2291,14 +2310,13 @@ end subroutine setConvergenceFlag
!> @brief determines whether a point is converged
!--------------------------------------------------------------------------------------------------
logical pure function converged(residuum,state,aTol)
use prec, only: &
dEq0
use numerics, only: &
rTol => rTol_crystalliteState
implicit none
real(pReal), intent(in), dimension(:) ::&
residuum, state, aTol
real(pReal) :: &
rTol
rTol = num%rTol_crystalliteState
converged = all(abs(residuum) <= max(aTol, rTol*abs(state)))
@ -2311,7 +2329,7 @@ end subroutine setConvergenceFlag
subroutine update_stress(timeFraction)
use mesh, only: &
mesh_element
implicit none
real(pReal), intent(in) :: &
timeFraction
integer :: &
@ -2341,13 +2359,12 @@ end subroutine update_stress
!--------------------------------------------------------------------------------------------------
!> @brief tbd
!--------------------------------------------------------------------------------------------------
subroutine update_dependentState()
subroutine update_dependentState
use mesh, only: &
mesh_element
use constitutive, only: &
constitutive_dependentState => constitutive_microstructure
implicit none
integer :: e, & ! element index in element loop
i, & ! integration point index in ip loop
g ! grain index in grain loop
@ -2378,7 +2395,6 @@ subroutine update_state(timeFraction)
use mesh, only: &
mesh_element
implicit none
real(pReal), intent(in) :: &
timeFraction
integer :: &
@ -2431,7 +2447,6 @@ subroutine update_dotState(timeFraction)
use constitutive, only: &
constitutive_collectDotState
implicit none
real(pReal), intent(in) :: &
timeFraction
integer :: &
@ -2490,7 +2505,6 @@ subroutine update_deltaState
phaseAt, phasememberAt
use constitutive, only: &
constitutive_collectDeltaState
implicit none
integer :: &
e, & !< element index in element loop
i, & !< integration point index in ip loop
@ -2579,7 +2593,6 @@ logical function stateJump(ipc,ip,el)
use constitutive, only: &
constitutive_collectDeltaState
implicit none
integer, intent(in):: &
el, & ! element index
ip, & ! integration point index

View File

@ -445,45 +445,45 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
steppingNeeded: if (materialpoint_subStep(i,e) > subStepMinHomog) then
! wind forward grain starting point of...
crystallite_partionedF0(1:3,1:3,1:myNgrains,i,e) = &
crystallite_partionedF0 (1:3,1:3,1:myNgrains,i,e) = &
crystallite_partionedF(1:3,1:3,1:myNgrains,i,e) ! ...def grads
crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) = &
crystallite_Fp(1:3,1:3,1:myNgrains,i,e) ! ...plastic def grads
crystallite_partionedFp0 (1:3,1:3,1:myNgrains,i,e) = &
crystallite_Fp (1:3,1:3,1:myNgrains,i,e) ! ...plastic def grads
crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) = &
crystallite_Lp(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_partionedLp0 (1:3,1:3,1:myNgrains,i,e) = &
crystallite_Lp (1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_partionedFi0(1:3,1:3,1:myNgrains,i,e) = &
crystallite_Fi(1:3,1:3,1:myNgrains,i,e) ! ...intermediate def grads
crystallite_partionedFi0 (1:3,1:3,1:myNgrains,i,e) = &
crystallite_Fi (1:3,1:3,1:myNgrains,i,e) ! ...intermediate def grads
crystallite_partionedLi0(1:3,1:3,1:myNgrains,i,e) = &
crystallite_Li(1:3,1:3,1:myNgrains,i,e) ! ...intermediate velocity grads
crystallite_partionedLi0 (1:3,1:3,1:myNgrains,i,e) = &
crystallite_Li (1:3,1:3,1:myNgrains,i,e) ! ...intermediate velocity grads
crystallite_partionedS0(1:3,1:3,1:myNgrains,i,e) = &
crystallite_S(1:3,1:3,1:myNgrains,i,e) ! ...2nd PK stress
crystallite_partionedS0 (1:3,1:3,1:myNgrains,i,e) = &
crystallite_S (1:3,1:3,1:myNgrains,i,e) ! ...2nd PK stress
do g = 1,myNgrains
plasticState (phaseAt(g,i,e))%partionedState0(:,phasememberAt(g,i,e)) = &
plasticState (phaseAt(g,i,e))%state( :,phasememberAt(g,i,e))
plasticState (phaseAt(g,i,e))%state (:,phasememberAt(g,i,e))
do mySource = 1, phase_Nsources(phaseAt(g,i,e))
sourceState(phaseAt(g,i,e))%p(mySource)%partionedState0(:,phasememberAt(g,i,e)) = &
sourceState(phaseAt(g,i,e))%p(mySource)%state( :,phasememberAt(g,i,e))
sourceState(phaseAt(g,i,e))%p(mySource)%state (:,phasememberAt(g,i,e))
enddo
enddo
forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), &
homogState(material_homogenizationAt(e))%sizeState > 0) &
homogState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = &
homogState(material_homogenizationAt(e))%State( :,mappingHomogenization(1,i,e)) ! ...internal homogenization state
homogState(material_homogenizationAt(e))%State (:,mappingHomogenization(1,i,e)) ! ...internal homogenization state
forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), &
thermalState(material_homogenizationAt(e))%sizeState > 0) &
thermalState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = &
thermalState(material_homogenizationAt(e))%State( :,mappingHomogenization(1,i,e)) ! ...internal thermal state
thermalState(material_homogenizationAt(e))%State (:,mappingHomogenization(1,i,e)) ! ...internal thermal state
forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), &
damageState(material_homogenizationAt(e))%sizeState > 0) &
damageState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = &
damageState(material_homogenizationAt(e))%State( :,mappingHomogenization(1,i,e)) ! ...internal damage state
damageState(material_homogenizationAt(e))%State (:,mappingHomogenization(1,i,e)) ! ...internal damage state
materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e) ! ...def grad
endif steppingNeeded
@ -515,14 +515,16 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
!--------------------------------------------------------------------------------------------------
! restore...
if (materialpoint_subStep(i,e) < 1.0_pReal) then ! protect against fake cutback from \Delta t = 2 to 1. Maybe that "trick" is not necessary anymore at all? I.e. start with \Delta t = 1
crystallite_Lp(1:3,1:3,1:myNgrains,i,e) = &
crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_Li(1:3,1:3,1:myNgrains,i,e) = &
crystallite_partionedLi0(1:3,1:3,1:myNgrains,i,e) ! ...intermediate velocity grads
endif ! maybe protecting everything from overwriting (not only L) makes even more sense
crystallite_Fp(1:3,1:3,1:myNgrains,i,e) = &
crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic def grads
crystallite_Lp(1:3,1:3,1:myNgrains,i,e) = &
crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_Fi(1:3,1:3,1:myNgrains,i,e) = &
crystallite_partionedFi0(1:3,1:3,1:myNgrains,i,e) ! ...intermediate def grads
crystallite_Li(1:3,1:3,1:myNgrains,i,e) = &
crystallite_partionedLi0(1:3,1:3,1:myNgrains,i,e) ! ...intermediate velocity grads
crystallite_S(1:3,1:3,1:myNgrains,i,e) = &
crystallite_partionedS0(1:3,1:3,1:myNgrains,i,e) ! ...2nd PK stress
do g = 1, myNgrains

View File

@ -13,12 +13,7 @@ module numerics
integer(pInt), protected, public :: &
iJacoStiffness = 1_pInt, & !< frequency of stiffness update
iJacoLpresiduum = 1_pInt, & !< frequency of Jacobian update of residuum in Lp
nMPstate = 10_pInt, & !< materialpoint state loop limit
nCryst = 20_pInt, & !< crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst")
nState = 10_pInt, & !< state loop limit
nStress = 40_pInt, & !< stress loop limit
pert_method = 1_pInt, & !< method used in perturbation technique for tangent
randomSeed = 0_pInt, & !< fixed seeding for pseudo-random number generator, Default 0: use random seed
worldrank = 0_pInt, & !< MPI worldrank (/=0 for MPI simulations only)
worldsize = 1_pInt, & !< MPI worldsize (/=1 for MPI simulations only)
@ -26,20 +21,10 @@ module numerics
integer(4), protected, public :: &
DAMASK_NumThreadsInt = 0 !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive
real(pReal), protected, public :: &
relevantStrain = 1.0e-7_pReal, & !< strain increment considered significant (used by crystallite to determine whether strain inc is considered significant)
defgradTolerance = 1.0e-7_pReal, & !< deviation of deformation gradient that is still allowed (used by CPFEM to determine outdated ffn1)
pert_Fg = 1.0e-7_pReal, & !< strain perturbation for FEM Jacobi
subStepMinCryst = 1.0e-3_pReal, & !< minimum (relative) size of sub-step allowed during cutback in crystallite
subStepMinHomog = 1.0e-3_pReal, & !< minimum (relative) size of sub-step allowed during cutback in homogenization
subStepSizeCryst = 0.25_pReal, & !< size of first substep when cutback in crystallite
subStepSizeHomog = 0.25_pReal, & !< size of first substep when cutback in homogenization
subStepSizeLp = 0.5_pReal, & !< size of first substep when cutback in Lp calculation
subStepSizeLi = 0.5_pReal, & !< size of first substep when cutback in Li calculation
stepIncreaseCryst = 1.5_pReal, & !< increase of next substep size when previous substep converged in crystallite
stepIncreaseHomog = 1.5_pReal, & !< increase of next substep size when previous substep converged in homogenization
rTol_crystalliteState = 1.0e-6_pReal, & !< relative tolerance in crystallite state loop
rTol_crystalliteStress = 1.0e-6_pReal, & !< relative tolerance in crystallite stress loop
aTol_crystalliteStress = 1.0e-8_pReal, & !< absolute tolerance in crystallite stress loop, Default 1.0e-8: residuum is in Lp and hence strain is on this order
numerics_unitlength = 1.0_pReal, & !< determines the physical length of one computational length unit
absTol_RGC = 1.0e+4_pReal, & !< absolute tolerance of RGC residuum
relTol_RGC = 1.0e-3_pReal, & !< relative tolerance of RGC residuum
@ -57,8 +42,7 @@ module numerics
charLength = 1.0_pReal, & !< characteristic length scale for gradient problems
residualStiffness = 1.0e-6_pReal !< non-zero residual damage
logical, protected, public :: &
usePingPong = .true., &
numerics_timeSyncing = .false. !< flag indicating if time synchronization in crystallite is used for nonlocal plasticity
usePingPong = .true.
!--------------------------------------------------------------------------------------------------
! field parameters:
@ -194,48 +178,18 @@ subroutine numerics_init
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('relevantstrain')
relevantStrain = IO_floatValue(line,chunkPos,2_pInt)
case ('defgradtolerance')
defgradTolerance = IO_floatValue(line,chunkPos,2_pInt)
case ('ijacostiffness')
iJacoStiffness = IO_intValue(line,chunkPos,2_pInt)
case ('ijacolpresiduum')
iJacoLpresiduum = IO_intValue(line,chunkPos,2_pInt)
case ('pert_fg')
pert_Fg = IO_floatValue(line,chunkPos,2_pInt)
case ('pert_method')
pert_method = IO_intValue(line,chunkPos,2_pInt)
case ('nmpstate')
nMPstate = IO_intValue(line,chunkPos,2_pInt)
case ('ncryst')
nCryst = IO_intValue(line,chunkPos,2_pInt)
case ('nstate')
nState = IO_intValue(line,chunkPos,2_pInt)
case ('nstress')
nStress = IO_intValue(line,chunkPos,2_pInt)
case ('substepmincryst')
subStepMinCryst = IO_floatValue(line,chunkPos,2_pInt)
case ('substepsizecryst')
subStepSizeCryst = IO_floatValue(line,chunkPos,2_pInt)
case ('stepincreasecryst')
stepIncreaseCryst = IO_floatValue(line,chunkPos,2_pInt)
case ('substepsizelp')
subStepSizeLp = IO_floatValue(line,chunkPos,2_pInt)
case ('substepsizeli')
subStepSizeLi = IO_floatValue(line,chunkPos,2_pInt)
case ('substepminhomog')
subStepMinHomog = IO_floatValue(line,chunkPos,2_pInt)
case ('substepsizehomog')
subStepSizeHomog = IO_floatValue(line,chunkPos,2_pInt)
case ('stepincreasehomog')
stepIncreaseHomog = IO_floatValue(line,chunkPos,2_pInt)
case ('rtol_crystallitestate')
rTol_crystalliteState = IO_floatValue(line,chunkPos,2_pInt)
case ('rtol_crystallitestress')
rTol_crystalliteStress = IO_floatValue(line,chunkPos,2_pInt)
case ('atol_crystallitestress')
aTol_crystalliteStress = IO_floatValue(line,chunkPos,2_pInt)
case ('integrator')
numerics_integrator = IO_intValue(line,chunkPos,2_pInt)
case ('usepingpong')
@ -356,23 +310,8 @@ subroutine numerics_init
!--------------------------------------------------------------------------------------------------
! writing parameters to output
write(6,'(a24,1x,es8.1)') ' relevantStrain: ',relevantStrain
write(6,'(a24,1x,es8.1)') ' defgradTolerance: ',defgradTolerance
write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness
write(6,'(a24,1x,i8)') ' iJacoLpresiduum: ',iJacoLpresiduum
write(6,'(a24,1x,es8.1)') ' pert_Fg: ',pert_Fg
write(6,'(a24,1x,i8)') ' pert_method: ',pert_method
write(6,'(a24,1x,i8)') ' nCryst: ',nCryst
write(6,'(a24,1x,es8.1)') ' subStepMinCryst: ',subStepMinCryst
write(6,'(a24,1x,es8.1)') ' subStepSizeCryst: ',subStepSizeCryst
write(6,'(a24,1x,es8.1)') ' stepIncreaseCryst: ',stepIncreaseCryst
write(6,'(a24,1x,es8.1)') ' subStepSizeLp: ',subStepSizeLp
write(6,'(a24,1x,es8.1)') ' subStepSizeLi: ',subStepSizeLi
write(6,'(a24,1x,i8)') ' nState: ',nState
write(6,'(a24,1x,i8)') ' nStress: ',nStress
write(6,'(a24,1x,es8.1)') ' rTol_crystalliteState: ',rTol_crystalliteState
write(6,'(a24,1x,es8.1)') ' rTol_crystalliteStress: ',rTol_crystalliteStress
write(6,'(a24,1x,es8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress
write(6,'(a24,1x,i8)') ' integrator: ',numerics_integrator
write(6,'(a24,1x,L8)') ' use ping pong scheme: ',usepingpong
write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength
@ -452,28 +391,12 @@ subroutine numerics_init
!--------------------------------------------------------------------------------------------------
! sanity checks
if (relevantStrain <= 0.0_pReal) call IO_error(301_pInt,ext_msg='relevantStrain')
if (defgradTolerance <= 0.0_pReal) call IO_error(301_pInt,ext_msg='defgradTolerance')
if (iJacoStiffness < 1_pInt) call IO_error(301_pInt,ext_msg='iJacoStiffness')
if (iJacoLpresiduum < 1_pInt) call IO_error(301_pInt,ext_msg='iJacoLpresiduum')
if (pert_Fg <= 0.0_pReal) call IO_error(301_pInt,ext_msg='pert_Fg')
if (pert_method <= 0_pInt .or. pert_method >= 4_pInt) &
call IO_error(301_pInt,ext_msg='pert_method')
if (nMPstate < 1_pInt) call IO_error(301_pInt,ext_msg='nMPstate')
if (nCryst < 1_pInt) call IO_error(301_pInt,ext_msg='nCryst')
if (nState < 1_pInt) call IO_error(301_pInt,ext_msg='nState')
if (nStress < 1_pInt) call IO_error(301_pInt,ext_msg='nStress')
if (subStepMinCryst <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepMinCryst')
if (subStepSizeCryst <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeCryst')
if (stepIncreaseCryst <= 0.0_pReal) call IO_error(301_pInt,ext_msg='stepIncreaseCryst')
if (subStepSizeLp <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeLp')
if (subStepSizeLi <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeLi')
if (subStepMinHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepMinHomog')
if (subStepSizeHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeHomog')
if (stepIncreaseHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='stepIncreaseHomog')
if (rTol_crystalliteState <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteState')
if (rTol_crystalliteStress <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteStress')
if (aTol_crystalliteStress <= 0.0_pReal) call IO_error(301_pInt,ext_msg='aTol_crystalliteStress')
if (numerics_integrator <= 0_pInt .or. numerics_integrator >= 6_pInt) &
call IO_error(301_pInt,ext_msg='integrator')
if (numerics_unitlength <= 0.0_pReal) call IO_error(301_pInt,ext_msg='unitlength')