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

fixed bug in crystallite_FPI: stateDamper always has to be defined for each grain 

added some OMP FLUSH statements were necessary 
replaced openmp do by forall construct where possible; this is much safer and perhaps even as fast for small loops
This commit is contained in:
Christoph Kords 2012-11-07 15:43:29 +00:00
parent bbcffa668b
commit dad9922f54
6 changed files with 421 additions and 451 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

439
code/constitutive.f90
View File

@ -219,215 +219,209 @@ if (any(numerics_integrator == 5_pInt)) then
allocate(constitutive_RKCK45dotState(6,gMax,iMax,eMax)) allocate(constitutive_RKCK45dotState(6,gMax,iMax,eMax))
endif endif
!$OMP PARALLEL DO PRIVATE(myNgrains,myInstance) do e = 1_pInt,mesh_NcpElems ! loop over elements
do e = 1_pInt,mesh_NcpElems ! loop over elements myNgrains = homogenization_Ngrains(mesh_element(3,e))
myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = 1_pInt,FE_Nips(mesh_element(2,e)) ! loop over IPs
do i = 1_pInt,FE_Nips(mesh_element(2,e)) ! loop over IPs do g = 1_pInt,myNgrains ! loop over grains
do g = 1_pInt,myNgrains ! loop over grains select case(phase_elasticity(material_phase(g,i,e)))
select case(phase_elasticity(material_phase(g,i,e)))
case (constitutive_hooke_label) case (constitutive_hooke_label)
! valid elasticity but nothing to do ! valid elasticity but nothing to do
case default case default
call IO_error(200_pInt,ext_msg=trim(phase_elasticity(material_phase(g,i,e)))) ! unknown elasticity call IO_error(200_pInt,ext_msg=trim(phase_elasticity(material_phase(g,i,e)))) ! unknown elasticity
end select end select
myInstance = phase_plasticityInstance(material_phase(g,i,e)) myInstance = phase_plasticityInstance(material_phase(g,i,e))
select case(phase_plasticity(material_phase(g,i,e))) select case(phase_plasticity(material_phase(g,i,e)))
case (constitutive_none_label) case (constitutive_none_label)
allocate(constitutive_state0(g,i,e)%p(constitutive_none_sizeState(myInstance))) allocate(constitutive_state0(g,i,e)%p(constitutive_none_sizeState(myInstance)))
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_none_sizeState(myInstance))) allocate(constitutive_partionedState0(g,i,e)%p(constitutive_none_sizeState(myInstance)))
allocate(constitutive_subState0(g,i,e)%p(constitutive_none_sizeState(myInstance))) allocate(constitutive_subState0(g,i,e)%p(constitutive_none_sizeState(myInstance)))
allocate(constitutive_state(g,i,e)%p(constitutive_none_sizeState(myInstance))) allocate(constitutive_state(g,i,e)%p(constitutive_none_sizeState(myInstance)))
allocate(constitutive_state_backup(g,i,e)%p(constitutive_none_sizeState(myInstance))) allocate(constitutive_state_backup(g,i,e)%p(constitutive_none_sizeState(myInstance)))
allocate(constitutive_aTolState(g,i,e)%p(constitutive_none_sizeState(myInstance))) allocate(constitutive_aTolState(g,i,e)%p(constitutive_none_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_none_sizeDotState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_none_sizeDotState(myInstance)))
allocate(constitutive_deltaState(g,i,e)%p(constitutive_none_sizeDotState(myInstance))) allocate(constitutive_deltaState(g,i,e)%p(constitutive_none_sizeDotState(myInstance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_none_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_none_sizeDotState(myInstance)))
if (any(numerics_integrator == 1_pInt)) then if (any(numerics_integrator == 1_pInt)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_none_sizeDotState(myInstance))) allocate(constitutive_previousDotState(g,i,e)%p(constitutive_none_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_none_sizeDotState(myInstance))) allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_none_sizeDotState(myInstance)))
endif endif
if (any(numerics_integrator == 4_pInt)) then if (any(numerics_integrator == 4_pInt)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_none_sizeDotState(myInstance))) allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_none_sizeDotState(myInstance)))
endif endif
if (any(numerics_integrator == 5_pInt)) then if (any(numerics_integrator == 5_pInt)) then
do s = 1_pInt,6_pInt do s = 1_pInt,6_pInt
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_none_sizeDotState(myInstance))) allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_none_sizeDotState(myInstance)))
enddo enddo
endif endif
constitutive_state0(g,i,e)%p = constitutive_none_stateInit(myInstance) constitutive_state0(g,i,e)%p = constitutive_none_stateInit(myInstance)
constitutive_aTolState(g,i,e)%p = constitutive_none_aTolState(myInstance) constitutive_aTolState(g,i,e)%p = constitutive_none_aTolState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_none_sizeState(myInstance) constitutive_sizeState(g,i,e) = constitutive_none_sizeState(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_none_sizeDotState(myInstance) constitutive_sizeDotState(g,i,e) = constitutive_none_sizeDotState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_none_sizePostResults(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_none_sizePostResults(myInstance)
case (constitutive_j2_label) case (constitutive_j2_label)
allocate(constitutive_state0(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_state0(g,i,e)%p(constitutive_j2_sizeState(myInstance)))
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_partionedState0(g,i,e)%p(constitutive_j2_sizeState(myInstance)))
allocate(constitutive_subState0(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_subState0(g,i,e)%p(constitutive_j2_sizeState(myInstance)))
allocate(constitutive_state(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_state(g,i,e)%p(constitutive_j2_sizeState(myInstance)))
allocate(constitutive_state_backup(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_state_backup(g,i,e)%p(constitutive_j2_sizeState(myInstance)))
allocate(constitutive_aTolState(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_aTolState(g,i,e)%p(constitutive_j2_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
allocate(constitutive_deltaState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) allocate(constitutive_deltaState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
if (any(numerics_integrator == 1_pInt)) then if (any(numerics_integrator == 1_pInt)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) allocate(constitutive_previousDotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
endif endif
if (any(numerics_integrator == 4_pInt)) then if (any(numerics_integrator == 4_pInt)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
endif endif
if (any(numerics_integrator == 5_pInt)) then if (any(numerics_integrator == 5_pInt)) then
do s = 1_pInt,6_pInt do s = 1_pInt,6_pInt
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
enddo enddo
endif endif
constitutive_state0(g,i,e)%p = constitutive_j2_stateInit(myInstance) constitutive_state0(g,i,e)%p = constitutive_j2_stateInit(myInstance)
constitutive_aTolState(g,i,e)%p = constitutive_j2_aTolState(myInstance) constitutive_aTolState(g,i,e)%p = constitutive_j2_aTolState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_j2_sizeState(myInstance) constitutive_sizeState(g,i,e) = constitutive_j2_sizeState(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_j2_sizeDotState(myInstance) constitutive_sizeDotState(g,i,e) = constitutive_j2_sizeDotState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_j2_sizePostResults(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_j2_sizePostResults(myInstance)
case (constitutive_phenopowerlaw_label) case (constitutive_phenopowerlaw_label)
allocate(constitutive_state0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) allocate(constitutive_state0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance)))
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) allocate(constitutive_partionedState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance)))
allocate(constitutive_subState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) allocate(constitutive_subState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance)))
allocate(constitutive_state(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) allocate(constitutive_state(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance)))
allocate(constitutive_state_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) allocate(constitutive_state_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance)))
allocate(constitutive_aTolState(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) allocate(constitutive_aTolState(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
allocate(constitutive_deltaState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) allocate(constitutive_deltaState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
if (any(numerics_integrator == 1_pInt)) then if (any(numerics_integrator == 1_pInt)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) allocate(constitutive_previousDotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
endif endif
if (any(numerics_integrator == 4_pInt)) then if (any(numerics_integrator == 4_pInt)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
endif endif
if (any(numerics_integrator == 5_pInt)) then if (any(numerics_integrator == 5_pInt)) then
do s = 1_pInt,6_pInt do s = 1_pInt,6_pInt
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
enddo enddo
endif endif
constitutive_state0(g,i,e)%p = constitutive_phenopowerlaw_stateInit(myInstance) constitutive_state0(g,i,e)%p = constitutive_phenopowerlaw_stateInit(myInstance)
constitutive_aTolState(g,i,e)%p = constitutive_phenopowerlaw_aTolState(myInstance) constitutive_aTolState(g,i,e)%p = constitutive_phenopowerlaw_aTolState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_phenopowerlaw_sizeState(myInstance) constitutive_sizeState(g,i,e) = constitutive_phenopowerlaw_sizeState(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_phenopowerlaw_sizeDotState(myInstance) constitutive_sizeDotState(g,i,e) = constitutive_phenopowerlaw_sizeDotState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_phenopowerlaw_sizePostResults(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_phenopowerlaw_sizePostResults(myInstance)
case (constitutive_titanmod_label) case (constitutive_titanmod_label)
allocate(constitutive_state0(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) allocate(constitutive_state0(g,i,e)%p(constitutive_titanmod_sizeState(myInstance)))
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) allocate(constitutive_partionedState0(g,i,e)%p(constitutive_titanmod_sizeState(myInstance)))
allocate(constitutive_subState0(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) allocate(constitutive_subState0(g,i,e)%p(constitutive_titanmod_sizeState(myInstance)))
allocate(constitutive_state(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) allocate(constitutive_state(g,i,e)%p(constitutive_titanmod_sizeState(myInstance)))
allocate(constitutive_state_backup(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) allocate(constitutive_state_backup(g,i,e)%p(constitutive_titanmod_sizeState(myInstance)))
allocate(constitutive_aTolState(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) allocate(constitutive_aTolState(g,i,e)%p(constitutive_titanmod_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
allocate(constitutive_deltaState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) allocate(constitutive_deltaState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
if (any(numerics_integrator == 1_pInt)) then if (any(numerics_integrator == 1_pInt)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) allocate(constitutive_previousDotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
endif endif
if (any(numerics_integrator == 4_pInt)) then if (any(numerics_integrator == 4_pInt)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
endif endif
if (any(numerics_integrator == 5_pInt)) then if (any(numerics_integrator == 5_pInt)) then
do s = 1_pInt,6_pInt do s = 1_pInt,6_pInt
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
enddo enddo
endif endif
constitutive_state0(g,i,e)%p = constitutive_titanmod_stateInit(myInstance) constitutive_state0(g,i,e)%p = constitutive_titanmod_stateInit(myInstance)
constitutive_aTolState(g,i,e)%p = constitutive_titanmod_aTolState(myInstance) constitutive_aTolState(g,i,e)%p = constitutive_titanmod_aTolState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_titanmod_sizeState(myInstance) constitutive_sizeState(g,i,e) = constitutive_titanmod_sizeState(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_titanmod_sizeDotState(myInstance) constitutive_sizeDotState(g,i,e) = constitutive_titanmod_sizeDotState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_titanmod_sizePostResults(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_titanmod_sizePostResults(myInstance)
case (constitutive_dislotwin_label) case (constitutive_dislotwin_label)
allocate(constitutive_state0(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) allocate(constitutive_state0(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance)))
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) allocate(constitutive_partionedState0(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance)))
allocate(constitutive_subState0(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) allocate(constitutive_subState0(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance)))
allocate(constitutive_state(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) allocate(constitutive_state(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance)))
allocate(constitutive_state_backup(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) allocate(constitutive_state_backup(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance)))
allocate(constitutive_aTolState(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) allocate(constitutive_aTolState(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
allocate(constitutive_deltaState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) allocate(constitutive_deltaState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
if (any(numerics_integrator == 1_pInt)) then if (any(numerics_integrator == 1_pInt)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) allocate(constitutive_previousDotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
endif endif
if (any(numerics_integrator == 4_pInt)) then if (any(numerics_integrator == 4_pInt)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
endif endif
if (any(numerics_integrator == 5_pInt)) then if (any(numerics_integrator == 5_pInt)) then
do s = 1_pInt,6_pInt do s = 1_pInt,6_pInt
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
enddo enddo
endif endif
constitutive_state0(g,i,e)%p = constitutive_dislotwin_stateInit(myInstance) constitutive_state0(g,i,e)%p = constitutive_dislotwin_stateInit(myInstance)
constitutive_aTolState(g,i,e)%p = constitutive_dislotwin_aTolState(myInstance) constitutive_aTolState(g,i,e)%p = constitutive_dislotwin_aTolState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_dislotwin_sizeState(myInstance) constitutive_sizeState(g,i,e) = constitutive_dislotwin_sizeState(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_dislotwin_sizeDotState(myInstance) constitutive_sizeDotState(g,i,e) = constitutive_dislotwin_sizeDotState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_dislotwin_sizePostResults(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_dislotwin_sizePostResults(myInstance)
case (constitutive_nonlocal_label) case (constitutive_nonlocal_label)
select case(mesh_element(2,e)) select case(mesh_element(2,e))
case (1_pInt,6_pInt,7_pInt,8_pInt,9_pInt) case (1_pInt,6_pInt,7_pInt,8_pInt,9_pInt)
! all fine ! all fine
case default case default
call IO_error(253_pInt,e,i,g) call IO_error(253_pInt,e,i,g)
end select end select
allocate(constitutive_state0(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) allocate(constitutive_state0(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance)))
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) allocate(constitutive_partionedState0(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance)))
allocate(constitutive_subState0(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) allocate(constitutive_subState0(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance)))
allocate(constitutive_state(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) allocate(constitutive_state(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance)))
allocate(constitutive_state_backup(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) allocate(constitutive_state_backup(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance)))
allocate(constitutive_aTolState(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) allocate(constitutive_aTolState(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
allocate(constitutive_deltaState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) allocate(constitutive_deltaState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
if (any(numerics_integrator == 1_pInt)) then if (any(numerics_integrator == 1_pInt)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) allocate(constitutive_previousDotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
endif endif
if (any(numerics_integrator == 4_pInt)) then if (any(numerics_integrator == 4_pInt)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
endif endif
if (any(numerics_integrator == 5_pInt)) then if (any(numerics_integrator == 5_pInt)) then
do s = 1_pInt,6_pInt do s = 1_pInt,6_pInt
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
enddo enddo
endif endif
constitutive_aTolState(g,i,e)%p = constitutive_nonlocal_aTolState(myInstance) constitutive_aTolState(g,i,e)%p = constitutive_nonlocal_aTolState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_nonlocal_sizeState(myInstance) constitutive_sizeState(g,i,e) = constitutive_nonlocal_sizeState(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_nonlocal_sizeDotState(myInstance) constitutive_sizeDotState(g,i,e) = constitutive_nonlocal_sizeDotState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_nonlocal_sizePostResults(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_nonlocal_sizePostResults(myInstance)
case default case default
call IO_error(201_pInt,ext_msg=trim(phase_plasticity(material_phase(g,i,e)))) ! unknown plasticity call IO_error(201_pInt,ext_msg=trim(phase_plasticity(material_phase(g,i,e)))) ! unknown plasticity
end select end select
enddo
enddo enddo
enddo enddo
!$OMP END PARALLEL DO enddo
call constitutive_nonlocal_stateInit(constitutive_state0(1,1:iMax,1:eMax)) call constitutive_nonlocal_stateInit(constitutive_state0(1,1:iMax,1:eMax))
!$OMP PARALLEL DO PRIVATE(myNgrains) do e = 1_pInt,mesh_NcpElems ! loop over elements
do e = 1_pInt,mesh_NcpElems ! loop over elements myNgrains = homogenization_Ngrains(mesh_element(3,e))
myNgrains = homogenization_Ngrains(mesh_element(3,e)) forall(i = 1_pInt:FE_Nips(mesh_element(2,e)), g = 1_pInt:myNgrains)
do i = 1_pInt,FE_Nips(mesh_element(2,e)) ! loop over IPs constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p
do g = 1_pInt,myNgrains ! loop over grains constitutive_state(g,i,e)%p = constitutive_state0(g,i,e)%p ! need to be defined for first call of constitutive_microstructure in crystallite_init
constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p endforall
constitutive_state(g,i,e)%p = constitutive_state0(g,i,e)%p ! need to be defined for first call of constitutive_microstructure in crystallite_init enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
!----- write out state size file---------------- !----- write out state size file----------------
call IO_write_jobBinaryIntFile(777,'sizeStateConst', size(constitutive_sizeState)) call IO_write_jobBinaryIntFile(777,'sizeStateConst', size(constitutive_sizeState))
@ -467,7 +461,7 @@ end subroutine constitutive_init
!* - ip : current integration point * !* - ip : current integration point *
!* - el : current element * !* - el : current element *
!********************************************************************* !*********************************************************************
function constitutive_homogenizedC(ipc,ip,el) pure function constitutive_homogenizedC(ipc,ip,el)
use prec, only: pReal use prec, only: pReal
use material, only: phase_plasticity,material_phase use material, only: phase_plasticity,material_phase
@ -479,7 +473,7 @@ function constitutive_homogenizedC(ipc,ip,el)
use constitutive_nonlocal use constitutive_nonlocal
implicit none implicit none
integer(pInt) :: ipc,ip,el integer(pInt), intent(in) :: ipc,ip,el
real(pReal), dimension(6,6) :: constitutive_homogenizedC real(pReal), dimension(6,6) :: constitutive_homogenizedC
select case (phase_plasticity(material_phase(ipc,ip,el))) select case (phase_plasticity(material_phase(ipc,ip,el)))
@ -693,22 +687,26 @@ end subroutine constitutive_LpAndItsTangent
!* - ip : current integration point * !* - ip : current integration point *
!* - el : current element * !* - el : current element *
!************************************************************************ !************************************************************************
subroutine constitutive_TandItsTangent(T, dT_dFe, Fe, ipc, ip, el) pure subroutine constitutive_TandItsTangent(T, dT_dFe, Fe, ipc, ip, el)
use prec, only: pReal use prec, only: pReal
use material, only: phase_elasticity,material_phase use material, only: phase_elasticity,material_phase
implicit none implicit none
integer(pInt) :: ipc,ip,el
real(pReal), dimension(3,3) :: T, Fe
real(pReal), dimension(3,3,3,3) :: dT_dFe
select case (phase_elasticity(material_phase(ipc,ip,el))) integer(pInt), intent(in) :: ipc,ip,el
real(pReal), dimension(3,3), intent(in) :: Fe
case (constitutive_hooke_label) real(pReal), dimension(3,3), intent(out) :: T
call constitutive_hooke_TAndItsTangent(T, dT_dFe, Fe, ipc, ip, el) real(pReal), dimension(3,3,3,3), intent(out) :: dT_dFe
end select
select case (phase_elasticity(material_phase(ipc,ip,el)))
case (constitutive_hooke_label)
call constitutive_hooke_TandItsTangent(T, dT_dFe, Fe, ipc, ip, el)
end select
end subroutine constitutive_TandItsTangent end subroutine constitutive_TandItsTangent
@ -727,17 +725,24 @@ end subroutine constitutive_TandItsTangent
!* - ip : current integration point * !* - ip : current integration point *
!* - el : current element * !* - el : current element *
!************************************************************************ !************************************************************************
subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, Fe, g, i, e) pure subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, Fe, g, i, e)
use prec, only: p_vec use prec, only: p_vec
use math use math
implicit none
implicit none
!* Definition of variables !* Definition of variables
integer(pInt) g, i, e, p, o
real(pReal), dimension(3,3) :: T, Fe integer(pInt), intent(in) :: g, i, e
real(pReal), dimension(6,6) :: C_66 real(pReal), dimension(3,3), intent(in) :: Fe
real(pReal), dimension(3,3,3,3) :: dT_dFe, C
integer(pInt) p, o
real(pReal), dimension(3,3), intent(out) :: T
real(pReal), dimension(3,3,3,3), intent(out) :: dT_dFe
real(pReal), dimension(6,6) :: C_66
real(pReal), dimension(3,3,3,3) :: C
!* get elasticity tensor !* get elasticity tensor

6
code/constitutive_j2.f90
View File

@ -338,7 +338,7 @@ pure function constitutive_j2_aTolState(myInstance)
end function constitutive_j2_aTolState end function constitutive_j2_aTolState
function constitutive_j2_homogenizedC(state,ipc,ip,el) pure function constitutive_j2_homogenizedC(state,ipc,ip,el)
!********************************************************************* !*********************************************************************
!* homogenized elacticity matrix * !* homogenized elacticity matrix *
!* INPUT: * !* INPUT: *
@ -353,9 +353,9 @@ function constitutive_j2_homogenizedC(state,ipc,ip,el)
implicit none implicit none
integer(pInt), intent(in) :: ipc,ip,el integer(pInt), intent(in) :: ipc,ip,el
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state
integer(pInt) :: matID integer(pInt) :: matID
real(pReal), dimension(6,6) :: constitutive_j2_homogenizedC real(pReal), dimension(6,6) :: constitutive_j2_homogenizedC
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
matID = phase_plasticityInstance(material_phase(ipc,ip,el)) matID = phase_plasticityInstance(material_phase(ipc,ip,el))
constitutive_j2_homogenizedC = constitutive_j2_Cslip_66(1:6,1:6,matID) constitutive_j2_homogenizedC = constitutive_j2_Cslip_66(1:6,1:6,matID)
@ -537,7 +537,7 @@ pure function constitutive_j2_dotState(Tstar_v, Temperature, state, g, ip, el)
endif endif
! dotState ! dotState
constitutive_j2_dotState = hardening * gamma_dot constitutive_j2_dotState = hardening * gamma_dot
end function constitutive_j2_dotState end function constitutive_j2_dotState

4
code/constitutive_none.f90
View File

@ -270,7 +270,7 @@ pure function constitutive_none_aTolState(myInstance)
end function constitutive_none_aTolState end function constitutive_none_aTolState
function constitutive_none_homogenizedC(state,ipc,ip,el) pure function constitutive_none_homogenizedC(state,ipc,ip,el)
!********************************************************************* !*********************************************************************
!* homogenized elacticity matrix * !* homogenized elacticity matrix *
!* INPUT: * !* INPUT: *
@ -285,9 +285,9 @@ function constitutive_none_homogenizedC(state,ipc,ip,el)
implicit none implicit none
integer(pInt), intent(in) :: ipc,ip,el integer(pInt), intent(in) :: ipc,ip,el
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state
integer(pInt) :: matID integer(pInt) :: matID
real(pReal), dimension(6,6) :: constitutive_none_homogenizedC real(pReal), dimension(6,6) :: constitutive_none_homogenizedC
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
matID = phase_plasticityInstance(material_phase(ipc,ip,el)) matID = phase_plasticityInstance(material_phase(ipc,ip,el))
constitutive_none_homogenizedC = constitutive_none_Cslip_66(1:6,1:6,matID) constitutive_none_homogenizedC = constitutive_none_Cslip_66(1:6,1:6,matID)

2
code/constitutive_phenopowerlaw.f90
View File

@ -593,7 +593,7 @@ end function constitutive_phenopowerlaw_aTolState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief homogenized elacticity matrix !> @brief homogenized elacticity matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_phenopowerlaw_homogenizedC(state,ipc,ip,el) pure function constitutive_phenopowerlaw_homogenizedC(state,ipc,ip,el)
use prec, only: p_vec use prec, only: p_vec
use mesh, only: mesh_NcpElems,mesh_maxNips use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_plasticityInstance use material, only: homogenization_maxNgrains,material_phase, phase_plasticityInstance

362
code/crystallite.f90
View File

@ -312,23 +312,17 @@ subroutine crystallite_init(Temperature)
close(myFile) close(myFile)
!$OMP PARALLEL PRIVATE(myNgrains,myPhase,myMat,myStructure)
!$OMP DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over all cp elements
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over all cp elements myNgrains = homogenization_Ngrains(mesh_element(3,e)) ! look up homogenization-->grainCount
myNgrains = homogenization_Ngrains(mesh_element(3,e)) ! look up homogenization-->grainCount forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1_pInt:myNgrains)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element crystallite_Fp0(1:3,1:3,g,i,e) = math_EulerToR(material_EulerAngles(1:3,g,i,e)) ! plastic def gradient reflects init orientation
do g = 1_pInt,myNgrains crystallite_F0(1:3,1:3,g,i,e) = math_I3
crystallite_Fp0(1:3,1:3,g,i,e) = math_EulerToR(material_EulerAngles(1:3,g,i,e)) ! plastic def gradient reflects init orientation crystallite_localPlasticity(g,i,e) = phase_localPlasticity(material_phase(g,i,e))
crystallite_F0(1:3,1:3,g,i,e) = math_I3 crystallite_Fe(1:3,1:3,g,i,e) = math_transpose33(crystallite_Fp0(1:3,1:3,g,i,e))
crystallite_localPlasticity(g,i,e) = phase_localPlasticity(material_phase(g,i,e)) crystallite_Fp(1:3,1:3,g,i,e) = crystallite_Fp0(1:3,1:3,g,i,e)
!$OMP FLUSH(crystallite_Fp0) endforall
crystallite_Fe(1:3,1:3,g,i,e) = math_transpose33(crystallite_Fp0(1:3,1:3,g,i,e)) enddo
crystallite_Fp(1:3,1:3,g,i,e) = crystallite_Fp0(1:3,1:3,g,i,e)
enddo
enddo
enddo
!$OMP ENDDO
crystallite_partionedTemperature0 = Temperature ! isothermal assumption crystallite_partionedTemperature0 = Temperature ! isothermal assumption
crystallite_partionedFp0 = crystallite_Fp0 crystallite_partionedFp0 = crystallite_Fp0
crystallite_partionedF0 = crystallite_F0 crystallite_partionedF0 = crystallite_F0
@ -338,34 +332,31 @@ crystallite_requested = .true.
! Initialize crystallite_symmetryID(g,i,e) ! Initialize crystallite_symmetryID(g,i,e)
!$OMP DO do e = FEsolving_execElem(1),FEsolving_execElem(2)
do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Ngrains(mesh_element(3,e))
myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1_pInt,myNgrains
do g = 1_pInt,myNgrains myPhase = material_phase(g,i,e)
myPhase = material_phase(g,i,e) myMat = phase_plasticityInstance(myPhase)
myMat = phase_plasticityInstance(myPhase) select case (phase_plasticity(myPhase))
select case (phase_plasticity(myPhase)) case (constitutive_phenopowerlaw_label)
case (constitutive_phenopowerlaw_label) myStructure = constitutive_phenopowerlaw_structure(myMat)
myStructure = constitutive_phenopowerlaw_structure(myMat) case (constitutive_titanmod_label)
case (constitutive_titanmod_label) myStructure = constitutive_titanmod_structure(myMat)
myStructure = constitutive_titanmod_structure(myMat) case (constitutive_dislotwin_label)
case (constitutive_dislotwin_label) myStructure = constitutive_dislotwin_structure(myMat)
myStructure = constitutive_dislotwin_structure(myMat) case (constitutive_nonlocal_label)
case (constitutive_nonlocal_label) myStructure = constitutive_nonlocal_structure(myMat)
myStructure = constitutive_nonlocal_structure(myMat) case default
case default myStructure = -1_pInt ! does this happen for j2 material?
myStructure = -1_pInt ! does this happen for j2 material? end select
end select if (myStructure > 0_pInt) then
if (myStructure > 0_pInt) then crystallite_symmetryID(g,i,e) = lattice_symmetryType(myStructure) ! structure = 1(fcc) or 2(bcc) => 1; 3(hex)=>2
crystallite_symmetryID(g,i,e) = lattice_symmetryType(myStructure) ! structure = 1(fcc) or 2(bcc) => 1; 3(hex)=>2 endif
endif
enddo
enddo enddo
enddo enddo
!$OMP ENDDO enddo
!$OMP END PARALLEL
call crystallite_orientations() call crystallite_orientations()
crystallite_orientation0 = crystallite_orientation ! Store initial orientations for calculation of grain rotations crystallite_orientation0 = crystallite_orientation ! Store initial orientations for calculation of grain rotations
@ -529,9 +520,7 @@ integer(pInt) NiterationCrystallite, &
o, & o, &
p, & p, &
perturbation , & ! loop counter for forward,backward perturbation mode perturbation , & ! loop counter for forward,backward perturbation mode
myNgrains, & myNgrains
mySizeState, &
mySizeDotState
logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
convergenceFlag_backup convergenceFlag_backup
! local variables used for calculating analytic Jacobian ! local variables used for calculating analytic Jacobian
@ -568,31 +557,24 @@ endif
crystallite_subStep = 0.0_pReal crystallite_subStep = 0.0_pReal
!$OMP PARALLEL DO PRIVATE(myNgrains) do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e))
myNgrains = homogenization_Ngrains(mesh_element(3,e)) forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1_pInt:myNgrains, crystallite_requested(g,i,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed crystallite_subTemperature0(g,i,e) = crystallite_partionedTemperature0(g,i,e) ! ...temperature
do g = 1_pInt,myNgrains constitutive_subState0(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructure
if (crystallite_requested(g,i,e)) then ! initialize restoration point of ... crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_partionedFp0(1:3,1:3,g,i,e) ! ...plastic def grad
crystallite_subTemperature0(g,i,e) = crystallite_partionedTemperature0(g,i,e) ! ...temperature crystallite_subLp0(1:3,1:3,g,i,e) = crystallite_partionedLp0(1:3,1:3,g,i,e) ! ...plastic velocity grad
constitutive_subState0(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructure crystallite_dPdF0(1:3,1:3,1:3,1:3,g,i,e) = crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,g,i,e) ! ...stiffness
crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_partionedFp0(1:3,1:3,g,i,e) ! ...plastic def grad crystallite_subF0(1:3,1:3,g,i,e) = crystallite_partionedF0(1:3,1:3,g,i,e) ! ...def grad
crystallite_subLp0(1:3,1:3,g,i,e) = crystallite_partionedLp0(1:3,1:3,g,i,e) ! ...plastic velocity grad crystallite_subTstar0_v(1:6,g,i,e) = crystallite_partionedTstar0_v(1:6,g,i,e) !...2nd PK stress
crystallite_dPdF0(1:3,1:3,1:3,1:3,g,i,e) = crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,g,i,e) ! ...stiffness crystallite_subFe0(1:3,1:3,g,i,e) = math_mul33x33(crystallite_subF0(1:3,1:3,g,i,e), &
crystallite_subF0(1:3,1:3,g,i,e) = crystallite_partionedF0(1:3,1:3,g,i,e) ! ...def grad math_inv33(crystallite_subFp0(1:3,1:3,g,i,e))) ! only needed later on for stiffness calculation
crystallite_subTstar0_v(1:6,g,i,e) = crystallite_partionedTstar0_v(1:6,g,i,e) !...2nd PK stress crystallite_subFrac(g,i,e) = 0.0_pReal
crystallite_subFe0(1:3,1:3,g,i,e) = math_mul33x33(crystallite_subF0(1:3,1:3,g,i,e), & crystallite_subStep(g,i,e) = 1.0_pReal/subStepSizeCryst
math_inv33(crystallite_subFp0(1:3,1:3,g,i,e))) ! only needed later on for stiffness calculation crystallite_todo(g,i,e) = .true.
crystallite_converged(g,i,e) = .false. ! pretend failed step of twice the required size
crystallite_subFrac(g,i,e) = 0.0_pReal endforall
crystallite_subStep(g,i,e) = 1.0_pReal/subStepSizeCryst enddo
crystallite_todo(g,i,e) = .true.
crystallite_converged(g,i,e) = .false. ! pretend failed step of twice the required size
endif
enddo
enddo
enddo
!$OMP END PARALLEL DO
! --+>> CRYSTALLITE CUTBACK LOOP <<+-- ! --+>> CRYSTALLITE CUTBACK LOOP <<+--
@ -653,6 +635,7 @@ do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2
!$OMP FLUSH(crystallite_subStep) !$OMP FLUSH(crystallite_subStep)
crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) ! ...temperature crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) ! ...temperature
crystallite_Fp(1:3,1:3,g,i,e) = crystallite_subFp0(1:3,1:3,g,i,e) ! ...plastic def grad crystallite_Fp(1:3,1:3,g,i,e) = crystallite_subFp0(1:3,1:3,g,i,e) ! ...plastic def grad
!$OMP FLUSH(crystallite_Fp)
crystallite_invFp(1:3,1:3,g,i,e) = math_inv33(crystallite_Fp(1:3,1:3,g,i,e)) crystallite_invFp(1:3,1:3,g,i,e) = math_inv33(crystallite_Fp(1:3,1:3,g,i,e))
!$OMP FLUSH(crystallite_invFp) !$OMP FLUSH(crystallite_invFp)
crystallite_Lp(1:3,1:3,g,i,e) = crystallite_subLp0(1:3,1:3,g,i,e) ! ...plastic velocity grad crystallite_Lp(1:3,1:3,g,i,e) = crystallite_subLp0(1:3,1:3,g,i,e) ! ...plastic velocity grad
@ -714,39 +697,37 @@ enddo
! --+>> CHECK FOR NON-CONVERGED CRYSTALLITES <<+-- ! --+>> CHECK FOR NON-CONVERGED CRYSTALLITES <<+--
!$OMP PARALLEL DO PRIVATE(myNgrains,invFp,Fe_guess,Tstar) do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e))
myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains
do g = 1,myNgrains if (.not. crystallite_converged(g,i,e)) then ! respond fully elastically (might be not required due to becoming terminally ill anyway)
if (.not. crystallite_converged(g,i,e)) then ! respond fully elastically (might be not required due to becoming terminally ill anyway) if(iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
#ifndef _OPENMP !$OMP CRITICAL (write2out)
if(iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> no convergence : respond fully elastic at el ip g ',e,i,g
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> no convergence : respond fully elastic at el ip g ',e,i,g
write(6,*)
endif
#endif
invFp = math_inv33(crystallite_partionedFp0(1:3,1:3,g,i,e))
Fe_guess = math_mul33x33(crystallite_partionedF(1:3,1:3,g,i,e), invFp)
call constitutive_TandItsTangent(Tstar, junk2, Fe_guess,g,i,e)
crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(Fe_guess,math_mul33x33(Tstar,transpose(invFp)))
endif
#ifndef _OPENMP
if(iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_level(debug_crystallite),debug_levelSelective) /= 0_pInt)) then
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> central solution of cryst_StressAndTangent at el ip g ',e,i,g
write(6,*)
write(6,'(a,/,3(12x,3(f12.4,1x)/))') '<< CRYST >> P / MPa', math_transpose33(crystallite_P(1:3,1:3,g,i,e))/1.0e6_pReal
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Fp', math_transpose33(crystallite_Fp(1:3,1:3,g,i,e))
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Lp', math_transpose33(crystallite_Lp(1:3,1:3,g,i,e))
write(6,*) write(6,*)
!$OMP END CRITICAL (write2out)
endif endif
#endif invFp = math_inv33(crystallite_partionedFp0(1:3,1:3,g,i,e))
enddo Fe_guess = math_mul33x33(crystallite_partionedF(1:3,1:3,g,i,e), invFp)
call constitutive_TandItsTangent(Tstar, junk2, Fe_guess,g,i,e)
crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(Fe_guess,math_mul33x33(Tstar,transpose(invFp)))
endif
if(iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_level(debug_crystallite),debug_levelSelective) /= 0_pInt)) then
!$OMP CRITICAL (write2out)
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> central solution of cryst_StressAndTangent at el ip g ',e,i,g
write(6,*)
write(6,'(a,/,3(12x,3(f12.4,1x)/))') '<< CRYST >> P / MPa', math_transpose33(crystallite_P(1:3,1:3,g,i,e))/1.0e6_pReal
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Fp', math_transpose33(crystallite_Fp(1:3,1:3,g,i,e))
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Lp', math_transpose33(crystallite_Lp(1:3,1:3,g,i,e))
write(6,*)
!$OMP END CRITICAL (write2out)
endif
enddo enddo
enddo enddo
!$OMP END PARALLEL DO enddo
! --+>> STIFFNESS CALCULATION <<+-- ! --+>> STIFFNESS CALCULATION <<+--
@ -760,19 +741,15 @@ if(updateJaco) then
! --- BACKUP --- ! --- BACKUP ---
!$OMP PARALLEL DO PRIVATE(myNgrains,mySizeState,mySizeDotState) do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e))
myNgrains = homogenization_Ngrains(mesh_element(3,e)) forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed constitutive_state_backup(g,i,e)%p(1:constitutive_sizeState(g,i,e)) = &
do g = 1,myNgrains constitutive_state(g,i,e)%p(1:constitutive_sizeState(g,i,e)) ! remember unperturbed, converged state, ...
mySizeState = constitutive_sizeState(g,i,e) ! number of state variables for this grain constitutive_dotState_backup(g,i,e)%p(1:constitutive_sizeDotState(g,i,e)) = &
mySizeDotState = constitutive_sizeDotState(g,i,e) ! number of dotStates for this grain constitutive_dotState(g,i,e)%p(1:constitutive_sizeDotState(g,i,e)) ! ... dotStates, ...
constitutive_state_backup(g,i,e)%p(1:mySizeState) = & endforall
constitutive_state(g,i,e)%p(1:mySizeState) ! remember unperturbed, converged state, ... enddo
constitutive_dotState_backup(g,i,e)%p(1:mySizeDotState) = &
constitutive_dotState(g,i,e)%p(1:mySizeDotState) ! ... dotStates, ...
enddo; enddo; enddo
!$OMP END PARALLEL DO
Temperature_backup = crystallite_Temperature ! ... Temperature, ... Temperature_backup = crystallite_Temperature ! ... Temperature, ...
F_backup = crystallite_subF ! ... and kinematics F_backup = crystallite_subF ! ... and kinematics
Fp_backup = crystallite_Fp Fp_backup = crystallite_Fp
@ -792,31 +769,27 @@ if(updateJaco) then
if (iand(pert_method,perturbation) > 0) then ! mask for desired direction if (iand(pert_method,perturbation) > 0) then ! mask for desired direction
myPert = -pert_Fg * (-1.0_pReal)**perturbation ! set perturbation step myPert = -pert_Fg * (-1.0_pReal)**perturbation ! set perturbation step
do k = 1,3; do l = 1,3 ! ...alter individual components do k = 1,3; do l = 1,3 ! ...alter individual components
#ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a,2(1x,i1),1x,a)') '<< CRYST >> [[[[[[ Stiffness perturbation',k,l,']]]]]]' write(6,'(a,2(1x,i1),1x,a)') '<< CRYST >> [[[[[[ Stiffness perturbation',k,l,']]]]]]'
write(6,*) write(6,*)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
#endif
! --- INITIALIZE UNPERTURBED STATE --- ! --- INITIALIZE UNPERTURBED STATE ---
select case(numerics_integrator(numerics_integrationMode)) select case(numerics_integrator(numerics_integrationMode))
case(1_pInt) ! Fix-point method: restore to last converged state at end of subinc, since this is probably closest to perturbed state case(1_pInt) ! Fix-point method: restore to last converged state at end of subinc, since this is probably closest to perturbed state
!$OMP PARALLEL DO PRIVATE(myNgrains,mySizeState,mySizeDotState) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Ngrains(mesh_element(3,e))
myNgrains = homogenization_Ngrains(mesh_element(3,e)) forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) constitutive_state(g,i,e)%p(1:constitutive_sizeState(g,i,e)) = &
do g = 1_pInt,myNgrains constitutive_state_backup(g,i,e)%p(1:constitutive_sizeState(g,i,e))
mySizeState = constitutive_sizeState(g,i,e) constitutive_dotState(g,i,e)%p(1:constitutive_sizeDotState(g,i,e)) = &
mySizeDotState = constitutive_sizeDotState(g,i,e) constitutive_dotState_backup(g,i,e)%p(1:constitutive_sizeDotState(g,i,e))
constitutive_state(g,i,e)%p(1:mySizeState) = constitutive_state_backup(g,i,e)%p(1:mySizeState) endforall
constitutive_dotState(g,i,e)%p(1:mySizeDotState) = constitutive_dotState_backup(g,i,e)%p(1:mySizeDotState) enddo
enddo; enddo; enddo
!OMP END PARALLEL DO
crystallite_Temperature = Temperature_backup crystallite_Temperature = Temperature_backup
crystallite_Fp = Fp_backup crystallite_Fp = Fp_backup
crystallite_invFp = InvFp_backup crystallite_invFp = InvFp_backup
@ -825,17 +798,15 @@ if(updateJaco) then
crystallite_Tstar_v = Tstar_v_backup crystallite_Tstar_v = Tstar_v_backup
case(2_pInt,3_pInt) ! explicit Euler methods: nothing to restore (except for F), since we are only doing a stress integration step case(2_pInt,3_pInt) ! explicit Euler methods: nothing to restore (except for F), since we are only doing a stress integration step
case(4_pInt,5_pInt) ! explicit Runge-Kutta methods: restore to start of subinc, since we are doing a full integration of state and stress case(4_pInt,5_pInt) ! explicit Runge-Kutta methods: restore to start of subinc, since we are doing a full integration of state and stress
!$OMP PARALLEL DO PRIVATE(myNgrains,mySizeState,mySizeDotState) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Ngrains(mesh_element(3,e))
myNgrains = homogenization_Ngrains(mesh_element(3,e)) forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) constitutive_state(g,i,e)%p(1:constitutive_sizeState(g,i,e)) = &
do g = 1_pInt,myNgrains constitutive_subState0(g,i,e)%p(1:constitutive_sizeState(g,i,e))
mySizeState = constitutive_sizeState(g,i,e) constitutive_dotState(g,i,e)%p(1:constitutive_sizeDotState(g,i,e)) = &
mySizeDotState = constitutive_sizeDotState(g,i,e) constitutive_dotState_backup(g,i,e)%p(1:constitutive_sizeDotState(g,i,e))
constitutive_state(g,i,e)%p(1:mySizeState) = constitutive_subState0(g,i,e)%p(1:mySizeState) endforall
constitutive_dotState(g,i,e)%p(1:mySizeDotState) = constitutive_dotState_backup(g,i,e)%p(1:mySizeDotState) enddo
enddo; enddo; enddo
!OMP END PARALLEL DO
crystallite_Temperature = crystallite_subTemperature0 crystallite_Temperature = crystallite_subTemperature0
crystallite_Fp = crystallite_subFp0 crystallite_Fp = crystallite_subFp0
crystallite_Fe = crystallite_subFe0 crystallite_Fe = crystallite_subFe0
@ -865,21 +836,19 @@ if(updateJaco) then
call crystallite_integrateStateRKCK45() call crystallite_integrateStateRKCK45()
end select end select
!OMP PARALLEL DO PRIVATE(myNgrains) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Ngrains(mesh_element(3,e))
myNgrains = homogenization_Ngrains(mesh_element(3,e)) select case(perturbation)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) case(1_pInt)
do g = 1_pInt,myNgrains forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains, &
if (crystallite_requested(g,i,e) .and. crystallite_converged(g,i,e)) then ! converged state warrants stiffness update crystallite_requested(g,i,e) .and. crystallite_converged(g,i,e)) & ! converged state warrants stiffness update
select case(perturbation) dPdF_perturbation1(1:3,1:3,k,l,g,i,e) = (crystallite_P(1:3,1:3,g,i,e) - P_backup(1:3,1:3,g,i,e)) / myPert ! tangent dP_ij/dFg_kl
case(1_pInt) case(2_pInt)
dPdF_perturbation1(1:3,1:3,k,l,g,i,e) = (crystallite_P(1:3,1:3,g,i,e) - P_backup(1:3,1:3,g,i,e)) / myPert ! tangent dP_ij/dFg_kl forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains, &
case(2_pInt) crystallite_requested(g,i,e) .and. crystallite_converged(g,i,e)) & ! converged state warrants stiffness update
dPdF_perturbation2(1:3,1:3,k,l,g,i,e) = (crystallite_P(1:3,1:3,g,i,e) - P_backup(1:3,1:3,g,i,e)) / myPert ! tangent dP_ij/dFg_kl dPdF_perturbation2(1:3,1:3,k,l,g,i,e) = (crystallite_P(1:3,1:3,g,i,e) - P_backup(1:3,1:3,g,i,e)) / myPert ! tangent dP_ij/dFg_kl
end select end select
endif enddo
enddo; enddo; enddo
!OMP END PARALLEL DO
enddo; enddo ! k,l loop enddo; enddo ! k,l loop
endif endif
@ -888,43 +857,41 @@ if(updateJaco) then
! --- STIFFNESS ACCORDING TO PERTURBATION METHOD AND CONVERGENCE --- ! --- STIFFNESS ACCORDING TO PERTURBATION METHOD AND CONVERGENCE ---
!$OMP PARALLEL DO PRIVATE(myNgrains) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Ngrains(mesh_element(3,e))
myNgrains = homogenization_Ngrains(mesh_element(3,e)) select case(pert_method)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) case(1_pInt)
do g = 1_pInt,myNgrains forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains, &
if (crystallite_requested(g,i,e) .and. convergenceFlag_backup(g,i,e)) then ! central solution converged crystallite_requested(g,i,e) .and. convergenceFlag_backup(g,i,e)) & ! perturbation mode 1: central solution converged
select case(pert_method) crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = dPdF_perturbation1(1:3,1:3,1:3,1:3,g,i,e)
case(1_pInt) case(2_pInt)
crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = dPdF_perturbation1(1:3,1:3,1:3,1:3,g,i,e) forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains, &
case(2_pInt) crystallite_requested(g,i,e) .and. convergenceFlag_backup(g,i,e)) & ! perturbation mode 2: central solution converged
crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = dPdF_perturbation2(1:3,1:3,1:3,1:3,g,i,e) crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = dPdF_perturbation2(1:3,1:3,1:3,1:3,g,i,e)
case(3_pInt) case(3_pInt)
crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = 0.5_pReal* ( dPdF_perturbation1(1:3,1:3,1:3,1:3,g,i,e) & forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains, &
+ dPdF_perturbation2(1:3,1:3,1:3,1:3,g,i,e)) crystallite_requested(g,i,e) .and. convergenceFlag_backup(g,i,e)) & ! perturbation mode 3: central solution converged
end select crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = 0.5_pReal* ( dPdF_perturbation1(1:3,1:3,1:3,1:3,g,i,e) &
elseif (crystallite_requested(g,i,e) .and. .not. convergenceFlag_backup(g,i,e)) then ! central solution did not converge + dPdF_perturbation2(1:3,1:3,1:3,1:3,g,i,e))
crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = crystallite_fallbackdPdF(1:3,1:3,1:3,1:3,g,i,e) ! use (elastic) fallback end select
endif forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains, &
enddo crystallite_requested(g,i,e) .and. .not. convergenceFlag_backup(g,i,e)) ! for any pertubation mode: if central solution did not converge...
enddo crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = crystallite_fallbackdPdF(1:3,1:3,1:3,1:3,g,i,e) ! ...use (elastic) fallback
enddo endforall
!$OMP END PARALLEL DO enddo
! --- RESTORE --- ! --- RESTORE ---
!$OMP PARALLEL DO PRIVATE(myNgrains,mySizeState,mySizeDotState) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Ngrains(mesh_element(3,e))
myNgrains = homogenization_Ngrains(mesh_element(3,e)) forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) constitutive_state(g,i,e)%p(1:constitutive_sizeState(g,i,e)) = &
do g = 1_pInt,myNgrains constitutive_state_backup(g,i,e)%p(1:constitutive_sizeState(g,i,e))
mySizeState = constitutive_sizeState(g,i,e) constitutive_dotState(g,i,e)%p(1:constitutive_sizeDotState(g,i,e)) = &
mySizeDotState = constitutive_sizeDotState(g,i,e) constitutive_dotState_backup(g,i,e)%p(1:constitutive_sizeDotState(g,i,e))
constitutive_state(g,i,e)%p(1:mySizeState) = constitutive_state_backup(g,i,e)%p(1:mySizeState) endforall
constitutive_dotState(g,i,e)%p(1:mySizeDotState) = constitutive_dotState_backup(g,i,e)%p(1:mySizeDotState) enddo
enddo; enddo; enddo
!OMP END PARALLEL DO
crystallite_Temperature = Temperature_backup crystallite_Temperature = Temperature_backup
crystallite_subF = F_backup crystallite_subF = F_backup
crystallite_Fp = Fp_backup crystallite_Fp = Fp_backup
@ -2345,7 +2312,8 @@ real(pReal) dot_prod12, &
stateDamper, & ! damper for integration of state stateDamper, & ! damper for integration of state
temperatureResiduum temperatureResiduum
real(pReal), dimension(constitutive_maxSizeDotState) :: & real(pReal), dimension(constitutive_maxSizeDotState) :: &
stateResiduum stateResiduum, &
tempState
logical singleRun ! flag indicating computation for single (g,i,e) triple logical singleRun ! flag indicating computation for single (g,i,e) triple
singleRun = present(ee) .and. present(ii) .and. present(gg) singleRun = present(ee) .and. present(ii) .and. present(gg)
@ -2423,7 +2391,6 @@ endif
! --+>> STATE LOOP <<+-- ! --+>> STATE LOOP <<+--
statedamper = 1.0_pReal
NiterationState = 0_pInt NiterationState = 0_pInt
do while (any(crystallite_todo) .and. NiterationState < nState ) ! convergence loop for crystallite do while (any(crystallite_todo) .and. NiterationState < nState ) ! convergence loop for crystallite
NiterationState = NiterationState + 1_pInt NiterationState = NiterationState + 1_pInt
@ -2462,11 +2429,11 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
!$OMP ENDDO !$OMP ENDDO
#ifndef _OPENMP !$OMP CRITICAL (write2out)
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after stress integration' write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after stress integration'
endif endif
#endif !$OMP END CRITICAL (write2out)
! --- DOT STATE AND TEMPERATURE --- ! --- DOT STATE AND TEMPERATURE ---
@ -2500,7 +2467,7 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
! --- UPDATE STATE AND TEMPERATURE --- ! --- UPDATE STATE AND TEMPERATURE ---
!$OMP DO PRIVATE(dot_prod12,dot_prod22,statedamper,mySizeDotState,stateResiduum,temperatureResiduum) !$OMP DO PRIVATE(dot_prod12,dot_prod22,statedamper,mySizeDotState,stateResiduum,temperatureResiduum,tempState)
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
@ -2514,6 +2481,8 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
.and. ( dot_prod12 < 0.0_pReal & .and. ( dot_prod12 < 0.0_pReal &
.or. dot_product(constitutive_dotState(g,i,e)%p, constitutive_previousDotState(g,i,e)%p) < 0.0_pReal) ) then .or. dot_product(constitutive_dotState(g,i,e)%p, constitutive_previousDotState(g,i,e)%p) < 0.0_pReal) ) then
statedamper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) statedamper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22)
else
statedamper = 1.0_pReal
endif endif
! --- get residui --- ! --- get residui ---
@ -2529,9 +2498,9 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
! --- correct state and temperature with residuum --- ! --- correct state and temperature with residuum ---
constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) & tempState(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) - stateResiduum(1:mySizeDotState) ! need to copy to local variable, since we cant flush a pointer in openmp
- stateResiduum(1:mySizeDotState)
crystallite_Temperature(g,i,e) = crystallite_Temperature(g,i,e) - temperatureResiduum crystallite_Temperature(g,i,e) = crystallite_Temperature(g,i,e) - temperatureResiduum
!$OMP FLUSH(crystallite_Temperature)
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt & if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
@ -2542,7 +2511,7 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
write(6,*) write(6,*)
write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> state residuum',stateResiduum(1:mySizeDotState) write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> state residuum',stateResiduum(1:mySizeDotState)
write(6,*) write(6,*)
write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> new state',constitutive_state(g,i,e)%p(1:mySizeDotState) write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> new state',tempState(1:mySizeDotState)
write(6,*) write(6,*)
endif endif
#endif #endif
@ -2557,7 +2526,7 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
if ( all( abs(stateResiduum(1:mySizeDotState)) < constitutive_aTolState(g,i,e)%p(1:mySizeDotState) & if ( all( abs(stateResiduum(1:mySizeDotState)) < constitutive_aTolState(g,i,e)%p(1:mySizeDotState) &
.or. abs(stateResiduum(1:mySizeDotState)) < rTol_crystalliteState & .or. abs(stateResiduum(1:mySizeDotState)) < rTol_crystalliteState &
* abs(constitutive_state(g,i,e)%p(1:mySizeDotState)) ) & * abs(tempState(1:mySizeDotState)) ) &
.and. (abs(temperatureResiduum) < rTol_crystalliteTemperature * crystallite_Temperature(g,i,e) & .and. (abs(temperatureResiduum) < rTol_crystalliteTemperature * crystallite_Temperature(g,i,e) &
.or. crystallite_Temperature(g,i,e) == 0.0_pReal) ) then .or. crystallite_Temperature(g,i,e) == 0.0_pReal) ) then
crystallite_converged(g,i,e) = .true. ! ... converged per definitionem crystallite_converged(g,i,e) = .true. ! ... converged per definitionem
@ -2568,6 +2537,7 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
!$OMP END CRITICAL (distributionState) !$OMP END CRITICAL (distributionState)
endif endif
endif endif
constitutive_state(g,i,e)%p(1:mySizeDotState) = tempState(1:mySizeDotState) ! copy local backup to global pointer
endif endif
enddo; enddo; enddo enddo; enddo; enddo
@ -2594,13 +2564,13 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
!$OMP END PARALLEL !$OMP END PARALLEL
#ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL(write2out)
write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), & write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), &
' grains converged after state integration no. ', NiterationState ' grains converged after state integration no. ', NiterationState
write(6,*) write(6,*)
!$OMP END CRITICAL(write2out)
endif endif
#endif
! --- NON-LOCAL CONVERGENCE CHECK --- ! --- NON-LOCAL CONVERGENCE CHECK ---
@ -2612,14 +2582,14 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
endif endif
crystallite_todo = crystallite_todo .and. .not. crystallite_converged ! skip all converged crystallite_todo = crystallite_todo .and. .not. crystallite_converged ! skip all converged
#ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL(write2out)
write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_converged(:,:,:)),' grains converged after non-local check' write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_converged(:,:,:)),' grains converged after non-local check'
write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after state integration no. ',& write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after state integration no. ',&
NiterationState NiterationState
write(6,*) write(6,*)
!$OMP END CRITICAL(write2out)
endif endif
#endif
enddo ! crystallite convergence loop enddo ! crystallite convergence loop

35
code/homogenization.f90
View File

@ -334,33 +334,28 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
endif endif
!$OMP PARALLEL DO PRIVATE(myNgrains) ! initialize restoration points of ...
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains)
constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p ! ...microstructures
! initialize restoration points of grain... crystallite_partionedTemperature0(g,i,e) = materialpoint_Temperature(i,e) ! ...temperatures
forall (g = 1:myNgrains) constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p ! ...microstructures crystallite_partionedFp0(1:3,1:3,g,i,e) = crystallite_Fp0(1:3,1:3,g,i,e) ! ...plastic def grads
crystallite_partionedTemperature0(1:myNgrains,i,e) = materialpoint_Temperature(i,e) ! ...temperatures crystallite_partionedLp0(1:3,1:3,g,i,e) = crystallite_Lp0(1:3,1:3,g,i,e) ! ...plastic velocity grads
crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Fp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic def grads crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,g,i,e) = crystallite_dPdF0(1:3,1:3,1:3,1:3,g,i,e) ! ...stiffness
crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Lp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads crystallite_partionedF0(1:3,1:3,g,i,e) = crystallite_F0(1:3,1:3,g,i,e) ! ...def grads
crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,1:myNgrains,i,e) = & crystallite_partionedTstar0_v(1:6,g,i,e) = crystallite_Tstar0_v(1:6,g,i,e) ! ...2nd PK stress
crystallite_dPdF0(1:3,1:3,1:3,1:3,1:myNgrains,i,e) ! ...stiffness endforall
crystallite_partionedF0(1:3,1:3,1:myNgrains,i,e) = crystallite_F0(1:3,1:3,1:myNgrains,i,e) ! ...def grads forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e))
crystallite_partionedTstar0_v(1:6,1:myNgrains,i,e) = crystallite_Tstar0_v(1:6,1:myNgrains,i,e) ! ...2nd PK stress
! initialize restoration points of ...
if (homogenization_sizeState(i,e) > 0_pInt) &
homogenization_subState0(i,e)%p = homogenization_state0(i,e)%p ! ...internal homogenization state
materialpoint_subF0(1:3,1:3,i,e) = materialpoint_F0(1:3,1:3,i,e) ! ...def grad materialpoint_subF0(1:3,1:3,i,e) = materialpoint_F0(1:3,1:3,i,e) ! ...def grad
materialpoint_subFrac(i,e) = 0.0_pReal materialpoint_subFrac(i,e) = 0.0_pReal
materialpoint_subStep(i,e) = 1.0_pReal/subStepSizeHomog ! <<added to adopt flexibility in cutback size>> materialpoint_subStep(i,e) = 1.0_pReal/subStepSizeHomog ! <<added to adopt flexibility in cutback size>>
materialpoint_converged(i,e) = .false. ! pretend failed step of twice the required size materialpoint_converged(i,e) = .false. ! pretend failed step of twice the required size
materialpoint_requested(i,e) = .true. ! everybody requires calculation materialpoint_requested(i,e) = .true. ! everybody requires calculation
enddo endforall
forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), homogenization_sizeState(i,e) > 0_pInt) &
homogenization_subState0(i,e)%p = homogenization_state0(i,e)%p ! ...internal homogenization state
enddo enddo
!$OMP END PARALLEL DO
NiterationHomog = 0_pInt NiterationHomog = 0_pInt