Merge branch 'development' into adaptive-euler-substate0
This commit is contained in:
Martin Diehl
commit
cff62219b6
|
|
@ -10,8 +10,8 @@
|
|||
*.pdf binary
|
||||
|
||||
# ignore files from MSC.Marc in language statistics
|
||||
installation/mods_MarcMentat/* linguist-vendored
|
||||
src/MarcInclude/* linguist-vendored
|
||||
installation/mods_MarcMentat/20*/* linguist-vendored
|
||||
src/marc/include/* linguist-vendored
|
||||
|
||||
# ignore reference files for tests in language statistics
|
||||
python/tests/reference/* linguist-vendored
|
||||
|
|
|
|||
|
|
@ -1048,7 +1048,7 @@ subroutine integrateStateFPI
|
|||
sourceState(p)%p(s)%previousDotState (:,c) = sourceState(p)%p(s)%dotState(:,c)
|
||||
enddo
|
||||
|
||||
call constitutive_dependentState(crystallite_Fe(1:3,1:3,g,i,e), &
|
||||
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), &
|
||||
crystallite_Fp(1:3,1:3,g,i,e), &
|
||||
g, i, e)
|
||||
|
||||
|
|
@ -1199,7 +1199,7 @@ subroutine integrateStateEuler
|
|||
nonlocalBroken = .true.
|
||||
if(.not. crystallite_todo(g,i,e)) cycle
|
||||
|
||||
call constitutive_dependentState(crystallite_Fe(1:3,1:3,g,i,e), &
|
||||
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), &
|
||||
crystallite_Fp(1:3,1:3,g,i,e), &
|
||||
g, i, e)
|
||||
|
||||
|
|
@ -1275,7 +1275,7 @@ subroutine integrateStateAdaptiveEuler
|
|||
nonlocalBroken = .true.
|
||||
if(.not. crystallite_todo(g,i,e)) cycle
|
||||
|
||||
call constitutive_dependentState(crystallite_Fe(1:3,1:3,g,i,e), &
|
||||
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), &
|
||||
crystallite_Fp(1:3,1:3,g,i,e), &
|
||||
g, i, e)
|
||||
|
||||
|
|
@ -1405,165 +1405,181 @@ end subroutine integrateStateRK4
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine integrateStateRKCK45
|
||||
|
||||
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, &
|
||||
.0_pReal, .225_pReal, -.9_pReal, 2.5_pReal, 175.0_pReal/512.0_pReal, &
|
||||
.0_pReal, .0_pReal, 1.2_pReal, -70.0_pReal/27.0_pReal, 575.0_pReal/13824.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], &
|
||||
[5,5], order=[2,1]) !< coefficients in Butcher tableau (used for preliminary integration in stages 2 to 6)
|
||||
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, &
|
||||
.0_pReal, .225_pReal, -.9_pReal, 2.5_pReal, 175.0_pReal/512.0_pReal, &
|
||||
.0_pReal, .0_pReal, 1.2_pReal, -70.0_pReal/27.0_pReal, 575.0_pReal/13824.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], &
|
||||
[5,5], order=[2,1]) !< coefficients in Butcher tableau (used for preliminary integration in stages 2 to 6)
|
||||
|
||||
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], & !< coefficients in Butcher tableau (used for final integration and error estimate)
|
||||
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, 0.25_pReal] !< coefficients in Butcher tableau (used for final integration and error estimate)
|
||||
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], & !< coefficients in Butcher tableau (used for final integration and error estimate)
|
||||
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, 0.25_pReal] !< coefficients in Butcher tableau (used for final integration and error estimate)
|
||||
|
||||
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)
|
||||
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)
|
||||
|
||||
integer :: &
|
||||
e, & ! element index in element loop
|
||||
i, & ! integration point index in ip loop
|
||||
g, & ! grain index in grain loop
|
||||
stage, & ! stage index in integration stage loop
|
||||
n, &
|
||||
p, &
|
||||
cc, &
|
||||
s, &
|
||||
sizeDotState
|
||||
integer :: &
|
||||
e, & ! element index in element loop
|
||||
i, & ! integration point index in ip loop
|
||||
g, & ! grain index in grain loop
|
||||
stage, & ! stage index in integration stage loop
|
||||
n, &
|
||||
p, &
|
||||
c, &
|
||||
s, &
|
||||
sizeDotState
|
||||
logical :: &
|
||||
nonlocalBroken
|
||||
|
||||
! ToDo: MD: once all constitutives use allocate state, attach residuum arrays to the state in case of RKCK45
|
||||
|
||||
real(pReal), dimension(constitutive_plasticity_maxSizeDotState, &
|
||||
homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
|
||||
residuum_plastic ! relative residuum from evolution in microstructure
|
||||
real(pReal), dimension(constitutive_source_maxSizeDotState, &
|
||||
maxval(phase_Nsources), &
|
||||
homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
|
||||
residuum_source ! relative residuum from evolution in microstructure
|
||||
real(pReal), dimension(constitutive_plasticity_maxSizeDotState, &
|
||||
homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
|
||||
residuum_plastic
|
||||
real(pReal), dimension(constitutive_source_maxSizeDotState, &
|
||||
maxval(phase_Nsources), &
|
||||
homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
|
||||
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)
|
||||
|
||||
do stage = 1,5
|
||||
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
|
||||
|
||||
! --- state update ---
|
||||
do stage = 1,5
|
||||
|
||||
!$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,:,c) = plasticState(p)%dotState(:,c)
|
||||
plasticState(p)%dotState(:,c) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,c)
|
||||
do s = 1, phase_Nsources(p)
|
||||
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)
|
||||
enddo
|
||||
|
||||
plasticState(p)%RKCK45dotState(stage,:,cc) = plasticState(p)%dotState(:,cc)
|
||||
plasticState(p)%dotState(:,cc) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,cc)
|
||||
do n = 2, stage
|
||||
plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) &
|
||||
+ A(n,stage) * plasticState(p)%RKCK45dotState(n,:,c)
|
||||
do s = 1, phase_Nsources(p)
|
||||
sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) &
|
||||
+ A(n,stage) * sourceState(p)%p(s)%RKCK45dotState(n,:,c)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do s = 1, phase_Nsources(p)
|
||||
sourceState(p)%p(s)%RKCK45dotState(stage,:,cc) = sourceState(p)%p(s)%dotState(:,cc)
|
||||
sourceState(p)%p(s)%dotState(:,cc) = A(1,stage) * sourceState(p)%p(s)%RKCK45dotState(1,:,cc)
|
||||
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
|
||||
|
||||
do n = 2, stage
|
||||
plasticState(p)%dotState(:,cc) = plasticState(p)%dotState(:,cc) &
|
||||
+ A(n,stage) * plasticState(p)%RKCK45dotState(n,:,cc)
|
||||
do s = 1, phase_Nsources(p)
|
||||
sourceState(p)%p(s)%dotState(:,cc) = sourceState(p)%p(s)%dotState(:,cc) &
|
||||
+ A(n,stage) * sourceState(p)%p(s)%RKCK45dotState(n,:,cc)
|
||||
enddo
|
||||
enddo
|
||||
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), &
|
||||
crystallite_Fp(1:3,1:3,g,i,e), &
|
||||
g, i, e)
|
||||
|
||||
endif
|
||||
enddo; enddo; enddo
|
||||
!$OMP END PARALLEL DO
|
||||
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
|
||||
|
||||
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))
|
||||
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
|
||||
|
||||
enddo
|
||||
enddo
|
||||
|
||||
if(.not. crystallite_todo(g,i,e)) cycle
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! --- STATE UPDATE WITH ERROR ESTIMATE FOR STATE ---
|
||||
sizeDotState = plasticState(p)%sizeDotState
|
||||
|
||||
!$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)
|
||||
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), &
|
||||
plasticState(p)%state(1:sizeDotState,c), &
|
||||
plasticState(p)%atol(1:sizeDotState))
|
||||
|
||||
sizeDotState = plasticState(p)%sizeDotState
|
||||
do s = 1, phase_Nsources(p)
|
||||
sizeDotState = sourceState(p)%p(s)%sizeDotState
|
||||
|
||||
plasticState(p)%RKCK45dotState(6,:,cc) = plasticState (p)%dotState(:,cc)
|
||||
sourceState(p)%p(s)%RKCK45dotState(6,:,c) = sourceState(p)%p(s)%dotState(:,c)
|
||||
residuum_source(1:sizeDotState,s,g,i,e) = matmul(DB,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)) &
|
||||
* 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))
|
||||
enddo
|
||||
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
|
||||
nonlocalBroken = .true.
|
||||
if(.not. crystallite_todo(g,i,e)) cycle
|
||||
|
||||
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_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
|
||||
|
||||
plasticState(p)%dotState(:,cc) = &
|
||||
matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,cc)), B) ! why transpose? Better to transpose constant B
|
||||
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), &
|
||||
crystallite_Fp(1:3,1:3,g,i,e), &
|
||||
g, i, e)
|
||||
|
||||
do s = 1, phase_Nsources(p)
|
||||
sizeDotState = sourceState(p)%p(s)%sizeDotState
|
||||
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
|
||||
nonlocalBroken = .true.
|
||||
if(.not. crystallite_todo(g,i,e)) cycle
|
||||
|
||||
sourceState(p)%p(s)%RKCK45dotState(6,:,cc) = sourceState(p)%p(s)%dotState(:,cc)
|
||||
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
|
||||
|
||||
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)
|
||||
endif
|
||||
enddo; enddo; enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
sourceState(p)%p(s)%dotState(:,cc) = &
|
||||
matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,cc)),B)
|
||||
enddo
|
||||
|
||||
endif
|
||||
enddo; enddo; enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
call update_state(1.0_pReal)
|
||||
|
||||
! --- relative residui and state convergence ---
|
||||
|
||||
!$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
|
||||
|
||||
crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), &
|
||||
plasticState(p)%state(1:sizeDotState,cc), &
|
||||
plasticState(p)%atol(1:sizeDotState))
|
||||
|
||||
do s = 1, phase_Nsources(p)
|
||||
sizeDotState = sourceState(p)%p(s)%sizeDotState
|
||||
|
||||
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,cc), &
|
||||
sourceState(p)%p(s)%atol(1:sizeDotState))
|
||||
enddo
|
||||
endif
|
||||
enddo; enddo; enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
call update_deltaState
|
||||
call update_dependentState
|
||||
call update_stress(1.0_pReal)
|
||||
call setConvergenceFlag
|
||||
if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck
|
||||
if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
|
||||
if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck
|
||||
|
||||
end subroutine integrateStateRKCK45
|
||||
|
||||
|
|
@ -1666,7 +1682,7 @@ subroutine update_dependentState
|
|||
do i = FEsolving_execIP(1),FEsolving_execIP(2)
|
||||
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
|
||||
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) &
|
||||
call constitutive_dependentState(crystallite_Fe(1:3,1:3,g,i,e), &
|
||||
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), &
|
||||
crystallite_Fp(1:3,1:3,g,i,e), &
|
||||
g, i, e)
|
||||
enddo; enddo; enddo
|
||||
|
|
|
|||
Loading…
Reference in New Issue