diff --git a/PRIVATE b/PRIVATE index aadf2d82a..183d7a3a3 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit aadf2d82a7e04646e3f20c3d526f27a6f90acef0 +Subproject commit 183d7a3a3bafa0a308c3eac858ca03c08fc03d50 diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 1df0f0d9f..3116345b6 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -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 diff --git a/src/homogenization.f90 b/src/homogenization.f90 index fe15364b8..ba2384632 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -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 diff --git a/src/numerics.f90 b/src/numerics.f90 index f7c603c60..77d1c7714 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -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')