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

all in one loop

This commit is contained in:
Martin Diehl 2020-03-25 11:07:47 +01:00
parent 9a188784e2
commit 32b9a5ab15
1 changed files with 32 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

55
src/crystallite.f90
View File

@ -1404,7 +1404,7 @@ end subroutine integrateStateRK4
!> adaptive step size (use 5th order solution to advance = "local extrapolation") !> adaptive step size (use 5th order solution to advance = "local extrapolation")
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateRKCK45 subroutine integrateStateRKCK45
real(pReal), dimension(5,5), parameter :: & real(pReal), dimension(5,5), parameter :: &
A = reshape([& A = reshape([&
.2_pReal, .075_pReal, .3_pReal, -11.0_pReal/54.0_pReal, 1631.0_pReal/55296.0_pReal, & .2_pReal, .075_pReal, .3_pReal, -11.0_pReal/54.0_pReal, 1631.0_pReal/55296.0_pReal, &
@ -1413,7 +1413,7 @@ subroutine integrateStateRKCK45
.0_pReal, .0_pReal, .0_pReal, 35.0_pReal/27.0_pReal, 44275.0_pReal/110592.0_pReal, & .0_pReal, .0_pReal, .0_pReal, 35.0_pReal/27.0_pReal, 44275.0_pReal/110592.0_pReal, &
.0_pReal, .0_pReal, .0_pReal, .0_pReal, 253.0_pReal/4096.0_pReal], & .0_pReal, .0_pReal, .0_pReal, .0_pReal, 253.0_pReal/4096.0_pReal], &
[5,5], order=[2,1]) !< coefficients in Butcher tableau (used for preliminary integration in stages 2 to 6) [5,5], order=[2,1]) !< coefficients in Butcher tableau (used for preliminary integration in stages 2 to 6)
real(pReal), dimension(6), parameter :: & real(pReal), dimension(6), parameter :: &
B = & B = &
[37.0_pReal/378.0_pReal, .0_pReal, 250.0_pReal/621.0_pReal, & [37.0_pReal/378.0_pReal, .0_pReal, 250.0_pReal/621.0_pReal, &
@ -1421,10 +1421,10 @@ subroutine integrateStateRKCK45
DB = B - & DB = B - &
[2825.0_pReal/27648.0_pReal, .0_pReal, 18575.0_pReal/48384.0_pReal,& [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, 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 :: &
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) 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
i, & ! integration point index in ip loop i, & ! integration point index in ip loop
@ -1437,7 +1437,7 @@ subroutine integrateStateRKCK45
sizeDotState sizeDotState
logical :: & logical :: &
nonlocalBroken 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 residuum_plastic
@ -1445,30 +1445,30 @@ subroutine integrateStateRKCK45
maxval(phase_Nsources), & maxval(phase_Nsources), &
homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
residuum_source residuum_source
call update_dotState(1.0_pReal) call update_dotState(1.0_pReal)
nonlocalBroken = .false. nonlocalBroken = .false.
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,c) !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2) do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
do stage = 1,5 do stage = 1,5
plasticState(p)%RKCK45dotState(stage,:,c) = plasticState(p)%dotState(:,c) plasticState(p)%RKCK45dotState(stage,:,c) = plasticState(p)%dotState(:,c)
plasticState(p)%dotState(:,c) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,c) plasticState(p)%dotState(:,c) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,c)
do s = 1, phase_Nsources(p) do s = 1, phase_Nsources(p)
sourceState(p)%p(s)%RKCK45dotState(stage,:,c) = sourceState(p)%p(s)%dotState(:,c) sourceState(p)%p(s)%RKCK45dotState(stage,:,c) = sourceState(p)%p(s)%dotState(:,c)
sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * sourceState(p)%p(s)%RKCK45dotState(1,:,c) 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(:,c) = plasticState(p)%dotState(:,c) & plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) &
+ A(n,stage) * plasticState(p)%RKCK45dotState(n,:,c) + A(n,stage) * plasticState(p)%RKCK45dotState(n,:,c)
@ -1477,7 +1477,7 @@ subroutine integrateStateRKCK45
+ A(n,stage) * sourceState(p)%p(s)%RKCK45dotState(n,:,c) + A(n,stage) * sourceState(p)%p(s)%RKCK45dotState(n,:,c)
enddo enddo
enddo enddo
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
+ plasticState(p)%dotState (1:sizeDotState,c) & + plasticState(p)%dotState (1:sizeDotState,c) &
@ -1488,16 +1488,16 @@ subroutine integrateStateRKCK45
+ sourceState(p)%p(s)%dotState (1:sizeDotState,c) & + sourceState(p)%p(s)%dotState (1:sizeDotState,c) &
* crystallite_subdt(g,i,e) * crystallite_subdt(g,i,e)
enddo enddo
call constitutive_dependentState(crystallite_Fe(1:3,1:3,g,i,e), & call constitutive_dependentState(crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), &
g, i, e) g, i, e)
crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage)) crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage))
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true. nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle if(.not. crystallite_todo(g,i,e)) cycle
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, & crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
@ -1510,9 +1510,9 @@ subroutine integrateStateRKCK45
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true. nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle if(.not. crystallite_todo(g,i,e)) cycle
enddo enddo
if(.not. crystallite_todo(g,i,e)) exit if(.not. crystallite_todo(g,i,e)) exit
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
@ -1537,14 +1537,23 @@ subroutine integrateStateRKCK45
matmul(B,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)) matmul(B,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c))
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
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false. if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
call update_state(1.0_pReal)
! --- relative residui and state convergence --- ! --- relative residui and state convergence ---
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,c) !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c)
@ -1575,7 +1584,7 @@ subroutine integrateStateRKCK45
call update_deltaState call update_deltaState
call update_dependentState call update_dependentState
call update_stress(1.0_pReal) call update_stress(1.0_pReal)
call setConvergenceFlag call setConvergenceFlag
if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck