223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch 'RKCK45-one-loop' into 'development' 

Rkck45 one loop

See merge request damask/DAMASK!147
This commit is contained in:
Franz Roters 2020-03-27 13:27:05 +01:00
parent 7030f8a923 adbf6739a9
commit fc0f2db8eb
1 changed files with 148 additions and 132 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

186
src/crystallite.f90
View File

@ -1423,7 +1423,7 @@ subroutine integrateStateRKCK45
13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 0.25_pReal] !< coefficients in Butcher tableau (used for final integration and error estimate) 13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 0.25_pReal] !< coefficients in Butcher tableau (used for final integration and error estimate)
real(pReal), dimension(5), parameter :: & real(pReal), dimension(5), parameter :: &
C = [0.2_pReal, 0.3_pReal, 0.6_pReal, 1.0_pReal, 0.875_pReal] !< coefficients in Butcher tableau (fractions of original time step in stages 2 to 6) CC = [0.2_pReal, 0.3_pReal, 0.6_pReal, 1.0_pReal, 0.875_pReal] !< coefficients in Butcher tableau (fractions of original time step in stages 2 to 6)
integer :: & integer :: &
e, & ! element index in element loop e, & ! element index in element loop
@ -1432,137 +1432,153 @@ subroutine integrateStateRKCK45
stage, & ! stage index in integration stage loop stage, & ! stage index in integration stage loop
n, & n, &
p, & p, &
cc, & c, &
s, & s, &
sizeDotState sizeDotState
logical :: &
! ToDo: MD: once all constitutives use allocate state, attach residuum arrays to the state in case of RKCK45 nonlocalBroken
real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & real(pReal), dimension(constitutive_plasticity_maxSizeDotState, &
homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
residuum_plastic ! relative residuum from evolution in microstructure residuum_plastic
real(pReal), dimension(constitutive_source_maxSizeDotState, & real(pReal), dimension(constitutive_source_maxSizeDotState, &
maxval(phase_Nsources), & maxval(phase_Nsources), &
homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
residuum_source ! relative residuum from evolution in microstructure residuum_source
call update_dotState(1.0_pReal) nonlocalBroken = .false.
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,c)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then
! --- SECOND TO SIXTH RUNGE KUTTA STEP --- p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do s = 1, phase_Nsources(p)
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
do stage = 1,5 do stage = 1,5
! --- state update --- plasticState(p)%RKCK45dotState(stage,:,c) = plasticState(p)%dotState(:,c)
plasticState(p)%dotState(:,c) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,c)
!$OMP PARALLEL DO PRIVATE(p,cc)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e)) then
p = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e)
plasticState(p)%RKCK45dotState(stage,:,cc) = plasticState(p)%dotState(:,cc)
plasticState(p)%dotState(:,cc) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,cc)
do s = 1, phase_Nsources(p) do s = 1, phase_Nsources(p)
sourceState(p)%p(s)%RKCK45dotState(stage,:,cc) = sourceState(p)%p(s)%dotState(:,cc) sourceState(p)%p(s)%RKCK45dotState(stage,:,c) = sourceState(p)%p(s)%dotState(:,c)
sourceState(p)%p(s)%dotState(:,cc) = A(1,stage) * sourceState(p)%p(s)%RKCK45dotState(1,:,cc) sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * sourceState(p)%p(s)%RKCK45dotState(1,:,c)
enddo enddo
do n = 2, stage do n = 2, stage
plasticState(p)%dotState(:,cc) = plasticState(p)%dotState(:,cc) & plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) &
+ A(n,stage) * plasticState(p)%RKCK45dotState(n,:,cc) + A(n,stage) * plasticState(p)%RKCK45dotState(n,:,c)
do s = 1, phase_Nsources(p) do s = 1, phase_Nsources(p)
sourceState(p)%p(s)%dotState(:,cc) = sourceState(p)%p(s)%dotState(:,cc) & sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) &
+ A(n,stage) * sourceState(p)%p(s)%RKCK45dotState(n,:,cc) + A(n,stage) * sourceState(p)%p(s)%RKCK45dotState(n,:,c)
enddo enddo
enddo enddo
endif
enddo; enddo; enddo
!$OMP END PARALLEL DO
call update_state(1.0_pReal) !MD: 1.0 correct?
call update_deltaState
call update_dependentState
call update_stress(C(stage))
call update_dotState(C(stage))
enddo
!--------------------------------------------------------------------------------------------------
! --- STATE UPDATE WITH ERROR ESTIMATE FOR STATE ---
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,cc)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e)) then
p = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e)
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
plasticState(p)%RKCK45dotState(6,:,cc) = plasticState (p)%dotState(:,cc) + plasticState(p)%dotState (1:sizeDotState,c) &
residuum_plastic(1:sizeDotState,g,i,e) = &
matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,cc)),DB) & ! why transpose? Better to transpose constant DB
* crystallite_subdt(g,i,e) * crystallite_subdt(g,i,e)
plasticState(p)%dotState(:,cc) = &
matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,cc)), B) ! why transpose? Better to transpose constant B
do s = 1, phase_Nsources(p) do s = 1, phase_Nsources(p)
sizeDotState = sourceState(p)%p(s)%sizeDotState 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)%RKCK45dotState(6,:,cc) = sourceState(p)%p(s)%dotState(:,cc) + sourceState(p)%p(s)%dotState (1:sizeDotState,c) &
residuum_source(1:sizeDotState,s,g,i,e) = &
matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,cc)),DB) &
* crystallite_subdt(g,i,e) * crystallite_subdt(g,i,e)
sourceState(p)%p(s)%dotState(:,cc) = &
matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,cc)),B)
enddo enddo
endif call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), &
enddo; enddo; enddo crystallite_Fp(1:3,1:3,g,i,e), &
!$OMP END PARALLEL DO g, i, e)
call update_state(1.0_pReal) crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage))
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) exit
! --- relative residui and state convergence --- call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, &
crystallite_subdt(g,i,e)*CC(stage), g,i,e)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do s = 1, phase_Nsources(p)
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) exit
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,cc) enddo
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2) if(.not. crystallite_todo(g,i,e)) cycle
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e)) then
p = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e)
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
plasticState(p)%RKCK45dotState(6,:,c) = plasticState (p)%dotState(:,c)
residuum_plastic(1:sizeDotState,g,i,e) = matmul(DB,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c)) &
* crystallite_subdt(g,i,e)
plasticState(p)%dotState(:,c) = matmul(B,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c))
plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
+ plasticState(p)%dotState (1:sizeDotState,c) &
* crystallite_subdt(g,i,e)
crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), & crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), &
plasticState(p)%state(1:sizeDotState,cc), & plasticState(p)%state(1:sizeDotState,c), &
plasticState(p)%atol(1:sizeDotState)) plasticState(p)%atol(1:sizeDotState))
do s = 1, phase_Nsources(p) do s = 1, phase_Nsources(p)
sizeDotState = sourceState(p)%p(s)%sizeDotState sizeDotState = sourceState(p)%p(s)%sizeDotState
crystallite_todo(g,i,e) = & sourceState(p)%p(s)%RKCK45dotState(6,:,c) = sourceState(p)%p(s)%dotState(:,c)
crystallite_todo(g,i,e) .and. converged(residuum_source(1:sizeDotState,s,g,i,e), & residuum_source(1:sizeDotState,s,g,i,e) = matmul(DB,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)) &
sourceState(p)%p(s)%state(1:sizeDotState,cc), & * crystallite_subdt(g,i,e)
sourceState(p)%p(s)%dotState(:,c) = matmul(B,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c))
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)
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. &
converged(residuum_source(1:sizeDotState,s,g,i,e), &
sourceState(p)%p(s)%state(1:sizeDotState,c), &
sourceState(p)%p(s)%atol(1:sizeDotState)) sourceState(p)%p(s)%atol(1:sizeDotState))
enddo enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
crystallite_todo(g,i,e) = stateJump(g,i,e)
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), &
g, i, e)
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
crystallite_todo(g,i,e) = integrateStress(g,i,e)
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
crystallite_converged(g,i,e) = crystallite_todo(g,i,e) ! consider converged if not broken
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
call update_deltaState if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
call update_dependentState
call update_stress(1.0_pReal)
call setConvergenceFlag
if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck
end subroutine integrateStateRKCK45 end subroutine integrateStateRKCK45