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

* Introduced preprocessor directives in order to suppress compilation of most write statements when using openmp. This tremendously improves efficiency of parallelization. 

* Also added some more openmp directives to increase percentage of parallelized code.

* "implicit none" was missing in two subroutines of homogenization and constitutive.
This commit is contained in:
Christoph Kords 2011-03-29 07:27:19 +00:00
parent 1b3dba4113
commit 3d51dd36fa
8 changed files with 842 additions and 821 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
code/CPFEM.f90
View File

@ -373,12 +373,14 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
endif
!$OMP END CRITICAL (write2out)
endif
do k = 1,mesh_NcpElems
do j = 1,mesh_maxNips
if (homogenization_sizeState(j,k) > 0_pInt) &
homogenization_state0(j,k)%p = homogenization_state(j,k)%p ! internal state of homogenization scheme
!$OMP PARALLEL DO
do k = 1,mesh_NcpElems
do j = 1,mesh_maxNips
if (homogenization_sizeState(j,k) > 0_pInt) &
homogenization_state0(j,k)%p = homogenization_state(j,k)%p ! internal state of homogenization scheme
enddo
enddo
enddo
!$OMP END PARALLEL DO
! *** dump the last converged values of each essential variable to a binary file
@ -492,16 +494,18 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
endif
call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent (parallel execution inside)
call materialpoint_postResults(dt) ! post results
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! loop over all parallely processed elements
if (microstructure_elemhomo(mesh_element(4,e))) then ! dealing with homogeneous element?
forall (i = 2:FE_Nips(mesh_element(2,e))) ! copy results of first IP to all others
materialpoint_P(1:3,1:3,i,e) = materialpoint_P(1:3,1:3,1,e)
materialpoint_F(1:3,1:3,i,e) = materialpoint_F(1:3,1:3,1,e)
materialpoint_dPdF(1:3,1:3,1:3,1:3,i,e) = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e)
materialpoint_results(1:materialpoint_sizeResults,i,e) = materialpoint_results(1:materialpoint_sizeResults,1,e)
end forall
endif
enddo
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! loop over all parallely processed elements
if (microstructure_elemhomo(mesh_element(4,e))) then ! dealing with homogeneous element?
forall (i = 2:FE_Nips(mesh_element(2,e))) ! copy results of first IP to all others
materialpoint_P(1:3,1:3,i,e) = materialpoint_P(1:3,1:3,1,e)
materialpoint_F(1:3,1:3,i,e) = materialpoint_F(1:3,1:3,1,e)
materialpoint_dPdF(1:3,1:3,1:3,1:3,i,e) = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e)
materialpoint_results(1:materialpoint_sizeResults,i,e) = materialpoint_results(1:materialpoint_sizeResults,1,e)
end forall
endif
enddo
!$OMP END PARALLEL DO
CPFEM_calc_done = .true.
endif
@ -511,7 +515,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
CPFEM_cs(1:6,IP,cp_en) = rnd * CPFEM_odd_stress
CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian * math_identity2nd(6)
else
! translate from P to CS
! translate from P to CS
Kirchhoff = math_mul33x33(materialpoint_P(1:3,1:3,IP, cp_en), math_transpose3x3(materialpoint_F(1:3,1:3,IP,cp_en)))
J_inverse = 1.0_pReal / math_det3x3(materialpoint_F(1:3,1:3,IP,cp_en))
CPFEM_cs(1:6,IP,cp_en) = math_Mandel33to6(J_inverse * Kirchhoff)

552
code/constitutive.f90
View File

@ -10,27 +10,27 @@
MODULE constitutive
!*** Include other modules ***
use prec
implicit none
use prec
implicit none
type(p_vec), dimension(:,:,:), allocatable :: constitutive_state0, & ! pointer array to microstructure at start of FE inc
constitutive_partionedState0, & ! pointer array to microstructure at start of homogenization inc
constitutive_subState0, & ! pointer array to microstructure at start of crystallite inc
constitutive_state, & ! pointer array to current microstructure (end of converged time step)
constitutive_state_backup, & ! pointer array to backed up microstructure (end of converged time step)
constitutive_dotState, & ! pointer array to evolution of current microstructure
constitutive_previousDotState,& ! pointer array to previous evolution of current microstructure
constitutive_previousDotState2,&! pointer array to 2nd previous evolution of current microstructure
constitutive_dotState_backup, & ! pointer array to backed up evolution of current microstructure
constitutive_RK4dotState, & ! pointer array to evolution of microstructure defined by classical Runge-Kutta method
constitutive_aTolState ! pointer array to absolute state tolerance
type(p_vec), dimension(:,:,:,:), allocatable :: constitutive_RKCK45dotState ! pointer array to evolution of microstructure used by Cash-Karp Runge-Kutta method
integer(pInt), dimension(:,:,:), allocatable :: constitutive_sizeDotState, & ! size of dotState array
constitutive_sizeState, & ! size of state array per grain
constitutive_sizePostResults ! size of postResults array per grain
integer(pInt) constitutive_maxSizeDotState, &
constitutive_maxSizeState, &
constitutive_maxSizePostResults
type(p_vec), dimension(:,:,:), allocatable :: constitutive_state0, & ! pointer array to microstructure at start of FE inc
constitutive_partionedState0, & ! pointer array to microstructure at start of homogenization inc
constitutive_subState0, & ! pointer array to microstructure at start of crystallite inc
constitutive_state, & ! pointer array to current microstructure (end of converged time step)
constitutive_state_backup, & ! pointer array to backed up microstructure (end of converged time step)
constitutive_dotState, & ! pointer array to evolution of current microstructure
constitutive_previousDotState,& ! pointer array to previous evolution of current microstructure
constitutive_previousDotState2,&! pointer array to 2nd previous evolution of current microstructure
constitutive_dotState_backup, & ! pointer array to backed up evolution of current microstructure
constitutive_RK4dotState, & ! pointer array to evolution of microstructure defined by classical Runge-Kutta method
constitutive_aTolState ! pointer array to absolute state tolerance
type(p_vec), dimension(:,:,:,:), allocatable :: constitutive_RKCK45dotState ! pointer array to evolution of microstructure used by Cash-Karp Runge-Kutta method
integer(pInt), dimension(:,:,:), allocatable :: constitutive_sizeDotState, & ! size of dotState array
constitutive_sizeState, & ! size of state array per grain
constitutive_sizePostResults ! size of postResults array per grain
integer(pInt) constitutive_maxSizeDotState, &
constitutive_maxSizeState, &
constitutive_maxSizePostResults
CONTAINS
!****************************************
@ -44,273 +44,299 @@ CONTAINS
!* - constitutive_postResults
!****************************************
subroutine constitutive_init()
!**************************************
!* Module initialization *
!**************************************
use prec, only: pReal,pInt
use debug, only: debug_verbosity, debug_selectiveDebugger, debug_e, debug_i, debug_g
use numerics, only: numerics_integrator
use IO, only: IO_error, IO_open_file, IO_open_jobFile
use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips
use material
use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_titanmod
use constitutive_dislotwin
use constitutive_nonlocal
subroutine constitutive_init()
use prec, only: pReal,pInt
use debug, only: debug_verbosity, debug_selectiveDebugger, debug_e, debug_i, debug_g
use numerics, only: numerics_integrator
use IO, only: IO_error, IO_open_file, IO_open_jobFile
use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips
use material
use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_titanmod
use constitutive_dislotwin
use constitutive_nonlocal
integer(pInt), parameter :: fileunit = 200
integer(pInt) e,i,g,p,s,myInstance,myNgrains
integer(pInt), dimension(:,:), pointer :: thisSize
character(len=64), dimension(:,:), pointer :: thisOutput
logical knownConstitution
implicit none
if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file
integer(pInt), parameter :: fileunit = 200
integer(pInt) g, & ! grain number
i, & ! integration point number
e, & ! element number
gMax, & ! maximum number of grains
iMax, & ! maximum number of integration points
eMax, & ! maximum number of elements
p, &
s, &
myInstance,&
myNgrains
integer(pInt), dimension(:,:), pointer :: thisSize
character(len=64), dimension(:,:), pointer :: thisOutput
logical knownConstitution
call constitutive_j2_init(fileunit) ! parse all phases of this constitution
call constitutive_phenopowerlaw_init(fileunit)
call constitutive_titanmod_init(fileunit)
call constitutive_dislotwin_init(fileunit)
call constitutive_nonlocal_init(fileunit)
close(fileunit)
! --- PARSE CONSTITUTIONS FROM CONFIG FILE ---
! write description file for constitutive phase output
if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file
call constitutive_j2_init(fileunit)
call constitutive_phenopowerlaw_init(fileunit)
call constitutive_titanmod_init(fileunit)
call constitutive_dislotwin_init(fileunit)
call constitutive_nonlocal_init(fileunit)
close(fileunit)
if(.not. IO_open_jobFile(fileunit,'outputConstitutive')) call IO_error (50) ! problems in writing file
do p = 1,material_Nphase
i = phase_constitutionInstance(p) ! which instance of a constitution is present phase
knownConstitution = .true. ! assume valid
select case(phase_constitution(p)) ! split per constitiution
case (constitutive_j2_label)
thisOutput => constitutive_j2_output
thisSize => constitutive_j2_sizePostResult
case (constitutive_phenopowerlaw_label)
thisOutput => constitutive_phenopowerlaw_output
thisSize => constitutive_phenopowerlaw_sizePostResult
case (constitutive_titanmod_label)
thisOutput => constitutive_titanmod_output
thisSize => constitutive_titanmod_sizePostResult
case (constitutive_dislotwin_label)
thisOutput => constitutive_dislotwin_output
thisSize => constitutive_dislotwin_sizePostResult
case (constitutive_nonlocal_label)
thisOutput => constitutive_nonlocal_output
thisSize => constitutive_nonlocal_sizePostResult
case default
knownConstitution = .false.
end select
! --- WRITE DESCRIPTION FILE FOR CONSTITUTIVE PHASE OUTPUT ---
write(fileunit,*)
write(fileunit,'(a)') '['//trim(phase_name(p))//']'
write(fileunit,*)
if (knownConstitution) then
write(fileunit,'(a)') '(constitution)'//char(9)//trim(phase_constitution(p))
do e = 1,phase_Noutput(p)
write(fileunit,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
enddo
endif
enddo
if(.not. IO_open_jobFile(fileunit,'outputConstitutive')) then ! problems in writing file
call IO_error (50)
endif
do p = 1,material_Nphase
i = phase_constitutionInstance(p) ! which instance of a constitution is present phase
knownConstitution = .true. ! assume valid
select case(phase_constitution(p)) ! split per constitiution
case (constitutive_j2_label)
thisOutput => constitutive_j2_output
thisSize => constitutive_j2_sizePostResult
case (constitutive_phenopowerlaw_label)
thisOutput => constitutive_phenopowerlaw_output
thisSize => constitutive_phenopowerlaw_sizePostResult
case (constitutive_titanmod_label)
thisOutput => constitutive_titanmod_output
thisSize => constitutive_titanmod_sizePostResult
case (constitutive_dislotwin_label)
thisOutput => constitutive_dislotwin_output
thisSize => constitutive_dislotwin_sizePostResult
case (constitutive_nonlocal_label)
thisOutput => constitutive_nonlocal_output
thisSize => constitutive_nonlocal_sizePostResult
case default
knownConstitution = .false.
end select
write(fileunit,*)
write(fileunit,'(a)') '['//trim(phase_name(p))//']'
write(fileunit,*)
if (knownConstitution) then
write(fileunit,'(a)') '(constitution)'//char(9)//trim(phase_constitution(p))
do e = 1,phase_Noutput(p)
write(fileunit,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
enddo
endif
enddo
close(fileunit)
close(fileunit)
! allocate memory for state management
! --- ALLOCATION OF STATES ---
allocate(constitutive_state0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_partionedState0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_subState0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_state(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_state_backup(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_dotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_dotState_backup(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_aTolState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_sizeDotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeDotState = 0_pInt
allocate(constitutive_sizeState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeState = 0_pInt
allocate(constitutive_sizePostResults(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)); constitutive_sizePostResults = 0_pInt
if (any(numerics_integrator == 1)) then
allocate(constitutive_previousDotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_previousDotState2(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
endif
if (any(numerics_integrator == 4)) &
allocate(constitutive_RK4dotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
if (any(numerics_integrator == 5)) &
allocate(constitutive_RKCK45dotState(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
gMax = homogenization_maxNgrains
iMax = mesh_maxNips
eMax = mesh_NcpElems
do e = 1,mesh_NcpElems ! loop over elements
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs
do g = 1,myNgrains ! loop over grains
myInstance = phase_constitutionInstance(material_phase(g,i,e))
select case(phase_constitution(material_phase(g,i,e)))
allocate(constitutive_state0(gMax,iMax,eMax))
allocate(constitutive_partionedState0(gMax,iMax,eMax))
allocate(constitutive_subState0(gMax,iMax,eMax))
allocate(constitutive_state(gMax,iMax,eMax))
allocate(constitutive_state_backup(gMax,iMax,eMax))
allocate(constitutive_dotState(gMax,iMax,eMax))
allocate(constitutive_dotState_backup(gMax,iMax,eMax))
allocate(constitutive_aTolState(gMax,iMax,eMax))
allocate(constitutive_sizeDotState(gMax,iMax,eMax)) ; constitutive_sizeDotState = 0_pInt
allocate(constitutive_sizeState(gMax,iMax,eMax)) ; constitutive_sizeState = 0_pInt
allocate(constitutive_sizePostResults(gMax,iMax,eMax)); constitutive_sizePostResults = 0_pInt
if (any(numerics_integrator == 1)) then
allocate(constitutive_previousDotState(gMax,iMax,eMax))
allocate(constitutive_previousDotState2(gMax,iMax,eMax))
endif
if (any(numerics_integrator == 4)) then
allocate(constitutive_RK4dotState(gMax,iMax,eMax))
endif
if (any(numerics_integrator == 5)) then
allocate(constitutive_RKCK45dotState(6,gMax,iMax,eMax))
endif
case (constitutive_j2_label)
allocate(constitutive_state0(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_state(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_dotState(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)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 4)) &
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
if (any(numerics_integrator == 5)) then
do s = 1,6
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_j2_stateInit(myInstance)
constitutive_aTolState(g,i,e)%p = constitutive_j2_aTolState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_j2_sizeState(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_j2_sizeDotState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_j2_sizePostResults(myInstance)
!$OMP PARALLEL DO PRIVATE(myNgrains,myInstance)
do e = 1,mesh_NcpElems ! loop over elements
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs
do g = 1,myNgrains ! loop over grains
myInstance = phase_constitutionInstance(material_phase(g,i,e))
select case(phase_constitution(material_phase(g,i,e)))
case (constitutive_phenopowerlaw_label)
allocate(constitutive_state0(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_state(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_dotState(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)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 4)) &
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
if (any(numerics_integrator == 5)) then
do s = 1,6
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_phenopowerlaw_stateInit(myInstance)
constitutive_aTolState(g,i,e)%p = constitutive_phenopowerlaw_aTolState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_phenopowerlaw_sizeState(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_phenopowerlaw_sizeDotState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_phenopowerlaw_sizePostResults(myInstance)
case (constitutive_j2_label)
allocate(constitutive_state0(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_state(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_dotState(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)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 4)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 5)) then
do s = 1,6
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_j2_stateInit(myInstance)
constitutive_aTolState(g,i,e)%p = constitutive_j2_aTolState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_j2_sizeState(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_j2_sizeDotState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_j2_sizePostResults(myInstance)
case (constitutive_titanmod_label)
allocate(constitutive_state0(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_state(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_dotState(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)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 4)) &
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
if (any(numerics_integrator == 5)) then
do s = 1,6
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_titanmod_stateInit(myInstance)
constitutive_aTolState(g,i,e)%p = constitutive_titanmod_aTolState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_titanmod_sizeState(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_titanmod_sizeDotState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_titanmod_sizePostResults(myInstance)
case (constitutive_phenopowerlaw_label)
allocate(constitutive_state0(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_state(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_dotState(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)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 4)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 5)) then
do s = 1,6
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_phenopowerlaw_stateInit(myInstance)
constitutive_aTolState(g,i,e)%p = constitutive_phenopowerlaw_aTolState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_phenopowerlaw_sizeState(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_phenopowerlaw_sizeDotState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_phenopowerlaw_sizePostResults(myInstance)
case (constitutive_dislotwin_label)
allocate(constitutive_state0(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_state(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_dotState(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)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 4)) &
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
if (any(numerics_integrator == 5)) then
do s = 1,6
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_dislotwin_stateInit(myInstance)
constitutive_aTolState(g,i,e)%p = constitutive_dislotwin_aTolState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_dislotwin_sizeState(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_dislotwin_sizeDotState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_dislotwin_sizePostResults(myInstance)
case (constitutive_titanmod_label)
allocate(constitutive_state0(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_state(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_dotState(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)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 4)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 5)) then
do s = 1,6
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_titanmod_stateInit(myInstance)
constitutive_aTolState(g,i,e)%p = constitutive_titanmod_aTolState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_titanmod_sizeState(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_titanmod_sizeDotState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_titanmod_sizePostResults(myInstance)
case (constitutive_nonlocal_label)
allocate(constitutive_state0(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_state(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_dotState(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)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 4)) &
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
if (any(numerics_integrator == 5)) then
do s = 1,6
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_nonlocal_stateInit(myInstance)
constitutive_aTolState(g,i,e)%p = constitutive_nonlocal_aTolState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_nonlocal_sizeState(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_nonlocal_sizeDotState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_nonlocal_sizePostResults(myInstance)
case (constitutive_dislotwin_label)
allocate(constitutive_state0(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_state(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_dotState(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)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 4)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 5)) then
do s = 1,6
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_dislotwin_stateInit(myInstance)
constitutive_aTolState(g,i,e)%p = constitutive_dislotwin_aTolState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_dislotwin_sizeState(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_dislotwin_sizeDotState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_dislotwin_sizePostResults(myInstance)
case default
call IO_error(200,material_phase(g,i,e)) ! unknown constitution
case (constitutive_nonlocal_label)
allocate(constitutive_state0(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_state(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_dotState(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)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 4)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 5)) then
do s = 1,6
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_nonlocal_stateInit(myInstance)
constitutive_aTolState(g,i,e)%p = constitutive_nonlocal_aTolState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_nonlocal_sizeState(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_nonlocal_sizeDotState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_nonlocal_sizePostResults(myInstance)
end select
constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p
enddo
enddo
enddo
case default
call IO_error(200,material_phase(g,i,e)) ! unknown constitution
constitutive_maxSizeState = maxval(constitutive_sizeState)
constitutive_maxSizeDotState = maxval(constitutive_sizeDotState)
constitutive_maxSizePostResults = maxval(constitutive_sizePostResults)
end select
constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p
enddo
enddo
enddo
!$OMP END PARALLEL DO
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) '<<<+- constitutive init -+>>>'
write(6,*) '$Id$'
write(6,*)
if (debug_verbosity > 0) then
write(6,'(a32,x,7(i5,x))') 'constitutive_state0: ', shape(constitutive_state0)
write(6,'(a32,x,7(i5,x))') 'constitutive_partionedState0: ', shape(constitutive_partionedState0)
write(6,'(a32,x,7(i5,x))') 'constitutive_subState0: ', shape(constitutive_subState0)
write(6,'(a32,x,7(i5,x))') 'constitutive_state: ', shape(constitutive_state)
write(6,'(a32,x,7(i5,x))') 'constitutive_aTolState: ', shape(constitutive_aTolState)
write(6,'(a32,x,7(i5,x))') 'constitutive_dotState: ', shape(constitutive_dotState)
write(6,'(a32,x,7(i5,x))') 'constitutive_sizeState: ', shape(constitutive_sizeState)
write(6,'(a32,x,7(i5,x))') 'constitutive_sizeDotState: ', shape(constitutive_sizeDotState)
write(6,'(a32,x,7(i5,x))') 'constitutive_sizePostResults: ', shape(constitutive_sizePostResults)
write(6,*)
write(6,'(a32,x,7(i5,x))') 'maxSizeState: ', constitutive_maxSizeState
write(6,'(a32,x,7(i5,x))') 'maxSizeDotState: ', constitutive_maxSizeDotState
write(6,'(a32,x,7(i5,x))') 'maxSizePostResults: ', constitutive_maxSizePostResults
endif
call flush(6)
!$OMP END CRITICAL (write2out)
return
constitutive_maxSizeState = maxval(constitutive_sizeState)
constitutive_maxSizeDotState = maxval(constitutive_sizeDotState)
constitutive_maxSizePostResults = maxval(constitutive_sizePostResults)
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) '<<<+- constitutive init -+>>>'
write(6,*) '$Id$'
write(6,*)
if (debug_verbosity > 0) then
write(6,'(a32,x,7(i5,x))') 'constitutive_state0: ', shape(constitutive_state0)
write(6,'(a32,x,7(i5,x))') 'constitutive_partionedState0: ', shape(constitutive_partionedState0)
write(6,'(a32,x,7(i5,x))') 'constitutive_subState0: ', shape(constitutive_subState0)
write(6,'(a32,x,7(i5,x))') 'constitutive_state: ', shape(constitutive_state)
write(6,'(a32,x,7(i5,x))') 'constitutive_aTolState: ', shape(constitutive_aTolState)
write(6,'(a32,x,7(i5,x))') 'constitutive_dotState: ', shape(constitutive_dotState)
write(6,'(a32,x,7(i5,x))') 'constitutive_sizeState: ', shape(constitutive_sizeState)
write(6,'(a32,x,7(i5,x))') 'constitutive_sizeDotState: ', shape(constitutive_sizeDotState)
write(6,'(a32,x,7(i5,x))') 'constitutive_sizePostResults: ', shape(constitutive_sizePostResults)
write(6,*)
write(6,'(a32,x,7(i5,x))') 'maxSizeState: ', constitutive_maxSizeState
write(6,'(a32,x,7(i5,x))') 'maxSizeDotState: ', constitutive_maxSizeDotState
write(6,'(a32,x,7(i5,x))') 'maxSizePostResults: ', constitutive_maxSizePostResults
endif
call flush(6)
!$OMP END CRITICAL (write2out)
endsubroutine

127
code/constitutive_nonlocal.f90
View File

@ -280,6 +280,7 @@ constitutive_nonlocal_interactionSlipSlip = 0.0_pReal
constitutive_nonlocal_dLowerEdgePerSlipFamily = 0.0_pReal
constitutive_nonlocal_dLowerScrewPerSlipFamily = 0.0_pReal
!*** readout data from material.config file
rewind(file)
@ -384,7 +385,8 @@ enddo
lattice_initializeStructure(constitutive_nonlocal_structureName(i), constitutive_nonlocal_CoverA(i)) ! our lattice structure is defined in the material.config file by the structureName (and the c/a ratio)
myStructure = constitutive_nonlocal_structure(i)
!*** sanity checks
!*** sanity checks
if (myStructure < 1 .or. myStructure > 3) call IO_error(205)
if (sum(constitutive_nonlocal_Nslip(:,i)) <= 0_pInt) call IO_error(235,ext_msg='Nslip')
@ -421,7 +423,7 @@ enddo
.or. constitutive_nonlocal_surfaceTransmissivity(i) > 1.0_pReal) call IO_error(235,ext_msg='surfaceTransmissivity')
!*** determine total number of active slip systems
!*** determine total number of active slip systems
constitutive_nonlocal_Nslip(1:lattice_maxNslipFamily,i) = min( lattice_NslipSystem(1:lattice_maxNslipFamily, myStructure), &
constitutive_nonlocal_Nslip(1:lattice_maxNslipFamily,i) ) ! we can't use more slip systems per family than specified in lattice
@ -471,7 +473,8 @@ do i = 1,maxNinstance
myStructure = constitutive_nonlocal_structure(i) ! lattice structure of this instance
!*** Inverse lookup of my slip system family and the slip system in lattice
!*** Inverse lookup of my slip system family and the slip system in lattice
l = 0_pInt
do f = 1,lattice_maxNslipFamily
@ -482,7 +485,7 @@ do i = 1,maxNinstance
enddo; enddo
!*** determine size of state array
!*** determine size of state array
ns = constitutive_nonlocal_totalNslip(i)
constitutive_nonlocal_sizeState(i) = size(constitutive_nonlocal_listBasicStates) * ns &
@ -490,7 +493,7 @@ do i = 1,maxNinstance
constitutive_nonlocal_sizeDotState(i) = size(constitutive_nonlocal_listBasicStates) * ns
!*** determine size of postResults array
!*** determine size of postResults array
do o = 1,maxval(phase_Noutput)
select case(constitutive_nonlocal_output(o,i))
@ -572,7 +575,7 @@ do i = 1,maxNinstance
enddo
!*** elasticity matrix and shear modulus according to material.config
!*** elasticity matrix and shear modulus according to material.config
select case (myStructure)
case(1:2) ! cubic(s)
@ -605,12 +608,10 @@ do i = 1,maxNinstance
/ ( 4.0_pReal*constitutive_nonlocal_C11(i) + 6.0_pReal*constitutive_nonlocal_C12(i) &
+ 2.0_pReal*constitutive_nonlocal_C44(i) ) ! C12iso/(C11iso+C12iso) with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5
do s1 = 1,ns
f = constitutive_nonlocal_slipFamily(s1,i)
!*** burgers vector, mean free path prefactor and minimum dipole distance for each slip system
!*** burgers vector, mean free path prefactor and minimum dipole distance for each slip system
constitutive_nonlocal_burgersPerSlipSystem(s1,i) = constitutive_nonlocal_burgersPerSlipFamily(f,i)
constitutive_nonlocal_lambda0PerSlipSystem(s1,i) = constitutive_nonlocal_lambda0PerSlipFamily(f,i)
@ -619,31 +620,33 @@ do i = 1,maxNinstance
do s2 = 1,ns
!*** calculation of forest projections for edge and screw dislocations. s2 acts as forest to s1
!*** calculation of forest projections for edge and screw dislocations. s2 acts as forest for s1
constitutive_nonlocal_forestProjectionEdge(s1, s2, i) &
= abs(math_mul3x3(lattice_sn(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), &
lattice_st(1:3, constitutive_nonlocal_slipSystemLattice(s2,i), myStructure))) ! forest projection of edge dislocations is the projection of (t = b x n) onto the slip normal of the respective slip plane
constitutive_nonlocal_forestProjectionEdge(s1,s2,i) &
= abs(math_mul3x3(lattice_sn(1:3,constitutive_nonlocal_slipSystemLattice(s1,i),myStructure), &
lattice_st(1:3,constitutive_nonlocal_slipSystemLattice(s2,i),myStructure))) ! forest projection of edge dislocations is the projection of (t = b x n) onto the slip normal of the respective slip plane
constitutive_nonlocal_forestProjectionScrew(s1, s2, i) &
= abs(math_mul3x3(lattice_sn(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), &
lattice_sd(1:3, constitutive_nonlocal_slipSystemLattice(s2,i), myStructure))) ! forest projection of screw dislocations is the projection of b onto the slip normal of the respective splip plane
constitutive_nonlocal_forestProjectionScrew(s1,s2,i) &
= abs(math_mul3x3(lattice_sn(1:3,constitutive_nonlocal_slipSystemLattice(s1,i),myStructure), &
lattice_sd(1:3,constitutive_nonlocal_slipSystemLattice(s2,i),myStructure))) ! forest projection of screw dislocations is the projection of b onto the slip normal of the respective splip plane
!*** calculation of interaction matrices
!*** calculation of interaction matrices
constitutive_nonlocal_interactionMatrixSlipSlip(s1, s2, i) &
= constitutive_nonlocal_interactionSlipSlip( lattice_interactionSlipSlip(constitutive_nonlocal_slipSystemLattice(s1,i), &
constitutive_nonlocal_slipSystemLattice(s2,i), &
myStructure), &
i )
constitutive_nonlocal_interactionMatrixSlipSlip(s1,s2,i) &
= constitutive_nonlocal_interactionSlipSlip(lattice_interactionSlipSlip(constitutive_nonlocal_slipSystemLattice(s1,i), &
constitutive_nonlocal_slipSystemLattice(s2,i), &
myStructure), &
i)
enddo; enddo
enddo
enddo
!*** calculation of prefactor for activation enthalpy for dislocation glide
!*** calculation of prefactor for activation enthalpy for dislocation glide
constitutive_nonlocal_Qeff0(1:ns,i) = constitutive_nonlocal_burgersPerSlipSystem(1:ns,i) ** 3.0_pReal &
* sqrt(0.5_pReal * constitutive_nonlocal_d0(i) ** 3.0_pReal &
* constitutive_nonlocal_Gmod(i) * constitutive_nonlocal_tauObs(i))
* sqrt(0.5_pReal * constitutive_nonlocal_d0(i) ** 3.0_pReal &
* constitutive_nonlocal_Gmod(i) * constitutive_nonlocal_tauObs(i))
enddo
@ -1029,16 +1032,6 @@ if (.not. phase_localConstitution(myPhase)) then
endif
if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
!$OMP CRITICAL (write2out)
write(6,*)
write(6,'(a,i5,x,i2,x,i1)') '<< CONST >> nonlocal_microstructure at el ip g',el,ip,g
write(6,*)
write(6,'(a,/,12(x),12(e10.3,x))') '<< CONST >> rhoForest', rhoForest
write(6,'(a,/,12(x),12(f10.5,x))') '<< CONST >> tauThreshold / MPa', tauThreshold/1e6
write(6,'(a,/,3(12(x),3(f10.5,x),/))') '<< CONST >> Tdislocation / MPa', math_Mandel6to33(Tdislocation_v)/1e6
!$OMP END CRITICAL (write2out)
endif
!**********************************************************************
!*** set states
@ -1049,6 +1042,18 @@ state(g,ip,el)%p(10*ns+1:11*ns) = rhoForest
state(g,ip,el)%p(11*ns+1:12*ns) = tauThreshold
state(g,ip,el)%p(12*ns+1:12*ns+6) = Tdislocation_v
#ifndef _OPENMP
if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,i5,x,i2,x,i1)') '<< CONST >> nonlocal_microstructure at el ip g',el,ip,g
write(6,*)
write(6,'(a,/,12(x),12(e10.3,x))') '<< CONST >> rhoForest', rhoForest
write(6,'(a,/,12(x),12(f10.5,x))') '<< CONST >> tauThreshold / MPa', tauThreshold/1e6
write(6,'(a,/,3(12(x),3(f10.5,x),/))') '<< CONST >> Tdislocation / MPa', math_Mandel6to33(Tdislocation_v)/1e6
endif
#endif
endsubroutine
@ -1166,15 +1171,15 @@ endif
constitutive_nonlocal_v(1:ns,1:4,g,ip,el) = v
!$OMP FLUSH(constitutive_nonlocal_v)
if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
!$OMP CRITICAL (write2out)
#ifndef _OPENMP
if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,i5,x,i2,x,i1)') '<< CONST >> nonlocal_kinetics at el ip g',el,ip,g
write(6,*)
write(6,'(a,/,12(x),12(f12.5,x))') '<< CONST >> tau / MPa', tau/1e6_pReal
write(6,'(a,/,4(12(x),12(f12.5,x),/))') '<< CONST >> v / 1e-3m/s', constitutive_nonlocal_v(:,:,g,ip,el)*1e3
!$OMP END CRITICAL (write2out)
endif
endif
#endif
endsubroutine
@ -1289,16 +1294,16 @@ do s = 1,ns
enddo
dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333)
if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
!$OMP CRITICAL (write2out)
#ifndef _OPENMP
if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,i5,x,i2,x,i1)') '<< CONST >> nonlocal_LpandItsTangent at el ip g ',el,ip,g
write(6,*)
write(6,'(a,/,4(12(x),12(f12.5,x)),/)') '<< CONST >> gdot / 1e-3',gdot*1e3_pReal
write(6,'(a,/,12(x),12(f12.5,x))') '<< CONST >> gdot total / 1e-3',gdotTotal*1e3_pReal
write(6,'(a,/,3(12(x),3(f12.7,x),/))') '<< CONST >> Lp',Lp
!$OMP END CRITICAL (write2out)
endif
endif
#endif
endsubroutine
@ -1441,21 +1446,19 @@ real(pReal) area, & ! area
logical considerEnteringFlux, &
considerLeavingFlux
if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
!$OMP CRITICAL (write2out)
#ifndef _OPENMP
if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,i5,x,i2,x,i1)') '<< CONST >> nonlocal_dotState at el ip g ',el,ip,g
write(6,*)
!$OMP END CRITICAL (write2out)
endif
endif
#endif
select case(mesh_element(2,el))
case (1,6,7,8,9)
! all fine
case default
call IO_error(-1,el,ip,g,'element type not supported for nonlocal constitution')
end select
myInstance = phase_constitutionInstance(material_phase(g,ip,el))
@ -1489,7 +1492,7 @@ endif
!****************************************************************************
!*** Calculate shear rate
call constitutive_nonlocal_kinetics(Tstar_v, Temperature, state(g,ip,el), g, ip, el) ! get velocities
call constitutive_nonlocal_kinetics(Tstar_v, Temperature, state(g,ip,el), g, ip, el) ! velocities
forall (t = 1:4) &
gdot(1:ns,t) = rhoSgl(1:ns,t) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
@ -1498,12 +1501,12 @@ forall (s = 1:ns, t = 1:4, rhoSgl(s,t+4) * constitutive_nonlocal_v(s,t,g,ip,el)
gdot(s,t) = gdot(s,t) + abs(rhoSgl(s,t+4)) * constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) &
* constitutive_nonlocal_v(s,t,g,ip,el)
if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
!$OMP CRITICAL (write2out)
#ifndef _OPENMP
if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
write(6,'(a,/,10(12(x),12(e12.5,x),/))') '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip
write(6,'(a,/,4(12(x),12(e12.5,x),/))') '<< CONST >> gdot / 1/s',gdot
!$OMP END CRITICAL (write2out)
endif
endif
#endif
@ -1751,8 +1754,10 @@ forall (t = 1:10) &
+ rhoDotAthermalAnnihilation(1:ns,t) &
+ rhoDotThermalAnnihilation(1:ns,t)
if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
!$OMP CRITICAL (write2out)
dotState%p(1:10*ns) = dotState%p(1:10*ns) + reshape(rhoDot,(/10*ns/))
#ifndef _OPENMP
if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
write(6,'(a,/,8(12(x),12(e12.5,x),/))') '<< CONST >> dislocation remobilization', rhoDotRemobilization(1:ns,1:8) * timestep
write(6,'(a,/,4(12(x),12(e12.5,x),/))') '<< CONST >> dislocation multiplication', rhoDotMultiplication(1:ns,1:4) * timestep
write(6,'(a,/,8(12(x),12(e12.5,x),/))') '<< CONST >> dislocation flux', rhoDotFlux(1:ns,1:8) * timestep
@ -1765,10 +1770,8 @@ if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g
write(6,'(a,/,10(12(x),12(f12.7,x),/))') '<< CONST >> relative density change', &
rhoDot(1:ns,1:8) * timestep / (abs(rhoSgl)+1.0e-10), &
rhoDot(1:ns,9:10) * timestep / (rhoDip+1.0e-10)
!$OMP END CRITICAL (write2out)
endif
dotState%p(1:10*ns) = dotState%p(1:10*ns) + reshape(rhoDot,(/10*ns/))
endif
#endif
endsubroutine
@ -2089,7 +2092,7 @@ forall (t = 1:4) &
do s = 1,ns
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
tau(s) = math_mul6x6( Tstar_v + Tdislocation_v, lattice_Sslip_v(1:6,sLattice,myStructure) )
tau(s) = math_mul6x6(Tstar_v + Tdislocation_v, lattice_Sslip_v(1:6,sLattice,myStructure))
enddo
dLower(1:ns,1) = constitutive_nonlocal_dLowerEdgePerSlipSystem(1:ns,myInstance)

434
code/crystallite.f90
View File

@ -15,7 +15,7 @@ MODULE crystallite
use prec, only: pReal, pInt
implicit none
!
! ****************************************************************
! *** General variables for the crystallite calculation ***
! ****************************************************************
@ -76,6 +76,7 @@ logical, dimension (:,:,:), allocatable :: &
CONTAINS
!********************************************************************
! allocate and initialize per grain variables
!********************************************************************
@ -213,7 +214,7 @@ allocate(crystallite_output(maxval(crystallite_Noutput), &
material_Ncrystallite)) ; crystallite_output = ''
allocate(crystallite_sizePostResults(material_Ncrystallite)) ; crystallite_sizePostResults = 0_pInt
allocate(crystallite_sizePostResult(maxval(crystallite_Noutput), &
material_Ncrystallite)) ; crystallite_sizePostResult = 0_pInt
material_Ncrystallite)) ; crystallite_sizePostResult = 0_pInt
if(.not. IO_open_file(file,material_configFile)) call IO_error (100) ! corrupt config file
@ -450,17 +451,13 @@ use math, only: math_inv3x3, &
math_transpose3x3, &
math_I3
use FEsolving, only: FEsolving_execElem, &
FEsolving_execIP, &
theInc, &
cycleCounter
FEsolving_execIP
use mesh, only: mesh_element, &
mesh_NcpElems, &
mesh_maxNips
use material, only: homogenization_Ngrains, &
homogenization_maxNgrains
use constitutive, only: constitutive_maxSizeState, &
constitutive_maxSizeDotState, &
constitutive_sizeState, &
use constitutive, only: constitutive_sizeState, &
constitutive_sizeDotState, &
constitutive_state, &
constitutive_state_backup, &
@ -468,10 +465,7 @@ use constitutive, only: constitutive_maxSizeState,
constitutive_partionedState0, &
constitutive_homogenizedC, &
constitutive_dotState, &
constitutive_dotState_backup, &
constitutive_collectDotState, &
constitutive_dotTemperature, &
constitutive_microstructure
constitutive_dotState_backup
implicit none
@ -517,21 +511,21 @@ logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
! --+>> INITIALIZE TO STARTING CONDITION <<+--
if (debug_verbosity > 4 .and. debug_e > 0 .and. debug_e <= mesh_NcpElems &
.and. debug_i > 0 .and. debug_i <= mesh_maxNips &
.and. debug_g > 0 .and. debug_g <= homogenization_maxNgrains) then
!$OMP CRITICAL (write2out)
write (6,*)
write (6,'(a,i5,x,i2,x,i3)') '<< CRYST >> crystallite start at el ip g ', debug_e, debug_i, debug_g
write (6,'(a,/,12(x),f14.9)') '<< CRYST >> Temp0', crystallite_partionedTemperature0(debug_g,debug_i,debug_e)
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> F0 ', &
math_transpose3x3(crystallite_partionedF0(1:3,1:3,debug_g,debug_i,debug_e))
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Fp0', &
math_transpose3x3(crystallite_partionedFp0(1:3,1:3,debug_g,debug_i,debug_e))
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Lp0', &
math_transpose3x3(crystallite_partionedLp0(1:3,1:3,debug_g,debug_i,debug_e))
!$OMP END CRITICAL (write2out)
endif
if (debug_verbosity > 4 .and. debug_e > 0 .and. debug_e <= mesh_NcpElems &
.and. debug_i > 0 .and. debug_i <= mesh_maxNips &
.and. debug_g > 0 .and. debug_g <= homogenization_maxNgrains) then
!$OMP CRITICAL (write2out)
write (6,*)
write (6,'(a,i5,x,i2,x,i3)') '<< CRYST >> crystallite start at el ip g ', debug_e, debug_i, debug_g
write (6,'(a,/,12(x),f14.9)') '<< CRYST >> Temp0', crystallite_partionedTemperature0(debug_g,debug_i,debug_e)
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> F0 ', &
math_transpose3x3(crystallite_partionedF0(1:3,1:3,debug_g,debug_i,debug_e))
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Fp0', &
math_transpose3x3(crystallite_partionedFp0(1:3,1:3,debug_g,debug_i,debug_e))
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Lp0', &
math_transpose3x3(crystallite_partionedLp0(1:3,1:3,debug_g,debug_i,debug_e))
!$OMP END CRITICAL (write2out)
endif
crystallite_subStep = 0.0_pReal
@ -575,15 +569,15 @@ do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2
! --- wind forward ---
if (crystallite_converged(g,i,e)) then
#ifndef _OPENMP
if (debug_verbosity > 4 &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
!$OMP CRITICAL (write2out)
write(6,'(a,f10.8,a,f10.8,a)') '<< CRYST >> winding forward from ', &
crystallite_subFrac(g,i,e),' to current crystallite_subfrac ', &
crystallite_subFrac(g,i,e)+crystallite_subStep(g,i,e),' in crystallite_stressAndItsTangent'
write(6,*)
!$OMP END CRITICAL (write2out)
write(6,'(a,f10.8,a,f10.8,a)') '<< CRYST >> winding forward from ', &
crystallite_subFrac(g,i,e),' to current crystallite_subfrac ', &
crystallite_subFrac(g,i,e)+crystallite_subStep(g,i,e),' in crystallite_stressAndItsTangent'
write(6,*)
endif
#endif
crystallite_subFrac(g,i,e) = crystallite_subFrac(g,i,e) + crystallite_subStep(g,i,e)
formerSubStep = crystallite_subStep(g,i,e)
!$OMP FLUSH(crystallite_subFrac)
@ -618,15 +612,15 @@ do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2
constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure
crystallite_Tstar_v(1:6,g,i,e) = crystallite_subTstar0_v(1:6,g,i,e) ! ...2nd PK stress
! cant restore dotState here, since not yet calculated in first cutback after initialization
!$OMP FLUSH(crystallite_subStep,crystallite_invFp)
!$OMP FLUSH(crystallite_invFp)
#ifndef _OPENMP
if (debug_verbosity > 4 &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
!$OMP CRITICAL (write2out)
write(6,'(a,f10.8)') '<< CRYST >> cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',&
crystallite_subStep(g,i,e)
write(6,*)
!$OMP END CRITICAL (write2out)
write(6,'(a,f10.8)') '<< CRYST >> cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',&
crystallite_subStep(g,i,e)
write(6,*)
endif
#endif
endif
! --- prepare for integration ---
@ -684,18 +678,17 @@ enddo
math_Mandel33to6( math_mul33x33(transpose(Fe_guess),Fe_guess) - math_I3 ) ) )
crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(Fe_guess,math_mul33x33(Tstar,transpose(invFp)))
endif
#ifndef _OPENMP
if (debug_verbosity > 4 &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
!$OMP FLUSH(crystallite_P)
!$OMP CRITICAL (write2out)
write (6,'(a,i5,x,i2,x,i3)') '<< CRYST >> central solution of cryst_StressAndTangent at el ip g ',e,i,g
write (6,*)
write (6,'(a,/,3(12(x),3(f12.4,x)/))') '<< CRYST >> P / MPa', math_transpose3x3(crystallite_P(1:3,1:3,g,i,e)) / 1e6
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Fp', math_transpose3x3(crystallite_Fp(1:3,1:3,g,i,e))
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Lp', math_transpose3x3(crystallite_Lp(1:3,1:3,g,i,e))
write (6,*)
!$OMP END CRITICAL (write2out)
write (6,'(a,i5,x,i2,x,i3)') '<< CRYST >> central solution of cryst_StressAndTangent at el ip g ',e,i,g
write (6,*)
write (6,'(a,/,3(12(x),3(f12.4,x)/))') '<< CRYST >> P / MPa', math_transpose3x3(crystallite_P(1:3,1:3,g,i,e)) / 1e6
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Fp', math_transpose3x3(crystallite_Fp(1:3,1:3,g,i,e))
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Lp', math_transpose3x3(crystallite_Lp(1:3,1:3,g,i,e))
write (6,*)
endif
#endif
enddo
enddo
enddo
@ -1007,18 +1000,18 @@ do n = 1,4
if (crystallite_todo(g,i,e)) then
if (crystallite_integrateStress(g,i,e,timeStepFraction(n))) then ! fraction of original times step
if (n == 4) then ! final integration step
#ifndef _OPENMP
if (debug_verbosity > 5 &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
mySizeDotState = constitutive_sizeDotState(g,i,e)
!$OMP CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState)
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
write(6,*)
!$OMP END CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState)
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
write(6,*)
endif
#endif
crystallite_converged(g,i,e) = .true. ! ... converged per definition
crystallite_todo(g,i,e) = .false. ! ... integration done
if (debug_verbosity > 0) then
@ -1238,13 +1231,11 @@ endif
! --- FIRST RUNGE KUTTA STEP ---
#ifndef _OPENMP
if (debug_verbosity > 5) then
!$OMP SINGLE
!$OMP CRITICAL (write2out)
write(6,'(a,x,i1)') '<< CRYST >> RUNGE KUTTA STEP',1
!$OMP END CRITICAL (write2out)
!$OMP END SINGLE
write(6,'(a,x,i1)') '<< CRYST >> RUNGE KUTTA STEP',1
endif
#endif
!$OMP DO
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
@ -1253,7 +1244,7 @@ endif
crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
crystallite_Temperature(g,i,e),g,i,e)
endif
enddo; enddo; enddo
enddo; enddo; enddo
!$OMP ENDDO
!$OMP DO
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
@ -1376,13 +1367,11 @@ do n = 1,5
! --- dot state and RK dot state---
#ifndef _OPENMP
if (debug_verbosity > 5) then
!$OMP SINGLE
!$OMP CRITICAL (write2out)
write(6,'(a,x,i1)') '<< CRYST >> RUNGE KUTTA STEP',n+1
!$OMP END CRITICAL (write2out)
!$OMP END SINGLE
write(6,'(a,x,i1)') '<< CRYST >> RUNGE KUTTA STEP',n+1
endif
#endif
!$OMP DO
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
@ -1502,24 +1491,23 @@ relTemperatureResiduum = 0.0_pReal
.or. abs(stateResiduum(1:mySizeDotState,g,i,e)) < constitutive_aTolState(g,i,e)%p(1:mySizeDotState) ) &
.and. abs(relTemperatureResiduum(g,i,e)) < rTol_crystalliteTemperature )
#ifndef _OPENMP
if (debug_verbosity > 5 &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
!$OMP CRITICAL (write2out)
write(6,'(a,i5,x,i3,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,*)
write(6,'(a,/,(12(x),12(f12.1,x)))') '<< CRYST >> absolute residuum tolerance', &
stateResiduum(1:mySizeDotState,g,i,e) / constitutive_aTolState(g,i,e)%p(1:mySizeDotState)
write(6,*)
write(6,'(a,/,(12(x),12(f12.1,x)))') '<< CRYST >> relative residuum tolerance', &
relStateResiduum(1:mySizeDotState,g,i,e) / rTol_crystalliteState
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState)
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
write(6,*)
!$OMP END CRITICAL (write2out)
write(6,'(a,i5,x,i3,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,*)
write(6,'(a,/,(12(x),12(f12.1,x)))') '<< CRYST >> absolute residuum tolerance', &
stateResiduum(1:mySizeDotState,g,i,e) / constitutive_aTolState(g,i,e)%p(1:mySizeDotState)
write(6,*)
write(6,'(a,/,(12(x),12(f12.1,x)))') '<< CRYST >> relative residuum tolerance', &
relStateResiduum(1:mySizeDotState,g,i,e) / rTol_crystalliteState
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState)
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
write(6,*)
endif
#endif
endif
enddo; enddo; enddo
!$OMP ENDDO
@ -1565,11 +1553,11 @@ relTemperatureResiduum = 0.0_pReal
! --- nonlocal convergence check ---
#ifndef _OPENMP
if (debug_verbosity > 5) then
!$OMP CRITICAL (write2out)
write(6,'(a,L)') '<< CRYST >> crystallite_converged',crystallite_converged
!$OMP END CRITICAL (write2out)
write(6,'(a,L)') '<< CRYST >> crystallite_converged',crystallite_converged
endif
#endif
if (.not. singleRun) then ! if not requesting Integration of just a single IP
if ( any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) then ! any non-local not yet converged (or broken)...
crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! ...restart all non-local as not converged
@ -1812,25 +1800,24 @@ relTemperatureResiduum = 0.0_pReal
relTemperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) / crystallite_Temperature(g,i,e)
!$OMP FLUSH(relStateResiduum,relTemperatureResiduum)
#ifndef _OPENMP
if (debug_verbosity > 5 &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
!$OMP CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,*)
write(6,'(a,/,(12(x),12(f12.1,x)))') '<< CRYST >> absolute residuum tolerance', &
stateResiduum(1:mySizeDotState,g,i,e) / constitutive_aTolState(g,i,e)%p(1:mySizeDotState)
write(6,*)
write(6,'(a,/,(12(x),12(f12.1,x)))') '<< CRYST >> relative residuum tolerance', &
relStateResiduum(1:mySizeDotState,g,i,e) / rTol_crystalliteState
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState) &
- 2.0_pReal * stateResiduum(1:mySizeDotState,g,i,e) / crystallite_subdt(g,i,e) ! calculate former dotstate from higher order solution and state residuum
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
write(6,*)
!$OMP END CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,*)
write(6,'(a,/,(12(x),12(f12.1,x)))') '<< CRYST >> absolute residuum tolerance', &
stateResiduum(1:mySizeDotState,g,i,e) / constitutive_aTolState(g,i,e)%p(1:mySizeDotState)
write(6,*)
write(6,'(a,/,(12(x),12(f12.1,x)))') '<< CRYST >> relative residuum tolerance', &
relStateResiduum(1:mySizeDotState,g,i,e) / rTol_crystalliteState
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState) &
- 2.0_pReal * stateResiduum(1:mySizeDotState,g,i,e) / crystallite_subdt(g,i,e) ! calculate former dotstate from higher order solution and state residuum
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
write(6,*)
endif
#endif
! --- converged ? ---
@ -1854,11 +1841,11 @@ relTemperatureResiduum = 0.0_pReal
! --- NONLOCAL CONVERGENCE CHECK ---
#ifndef _OPENMP
if (debug_verbosity > 5) then
!$OMP CRITICAL (write2out)
write(6,'(a,L)') '<< CRYST >> crystallite_converged',crystallite_converged
!$OMP END CRITICAL (write2out)
write(6,'(a,L)') '<< CRYST >> crystallite_converged',crystallite_converged
endif
#endif
if (.not. singleRun) then ! if not requesting Integration of just a single IP
if ( any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) then ! any non-local not yet converged (or broken)...
crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! ...restart all non-local as not converged
@ -1988,16 +1975,16 @@ endif
endif
enddo; enddo; enddo
!$OMP ENDDO
#ifndef _OPENMP
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
!$OMP CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(debug_g,debug_i,debug_e)%p(1:mySizeDotState)
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(debug_g,debug_i,debug_e)%p(1:mySizeDotState)
write(6,*)
!$OMP END CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(debug_g,debug_i,debug_e)%p(1:mySizeDotState)
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(debug_g,debug_i,debug_e)%p(1:mySizeDotState)
write(6,*)
endif
#endif
! --- UPDATE DEPENDENT STATES ---
@ -2104,7 +2091,9 @@ real(pReal) dot_prod12, &
dot_prod22
logical singleRun, & ! flag indicating computation for single (g,i,e) triple
stateConverged, & ! flag indicating convergence of state integration
temperatureConverged ! flag indicating convergence of temperature integration
temperatureConverged, & ! flag indicating convergence of temperature integration
stateUpdateDone, & ! flag indicating successfull state update
temperatureUpdateDone ! flag indicating successfull temperature update
if (present(ee) .and. present(ii) .and. present(gg)) then
@ -2126,7 +2115,7 @@ endif
! --- RESET DOTSTATE ---
!$OMP PARALLEL PRIVATE(stateConverged,temperatureConverged)
!$OMP PARALLEL
!$OMP DO
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
@ -2154,14 +2143,14 @@ endif
!$OMP SINGLE
crystallite_statedamper = 1.0_pReal
!$OMP END SINGLE
!$OMP DO
!$OMP DO PRIVATE(stateUpdateDone,temperatureUpdateDone,stateConverged,temperatureConverged)
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
call crystallite_updateState(crystallite_todo(g,i,e), stateConverged, g,i,e) ! update state
!$OMP FLUSH(crystallite_todo)
call crystallite_updateTemperature(crystallite_todo(g,i,e), temperatureConverged, g,i,e) ! update temperature
!$OMP FLUSH(crystallite_todo)
if ( .not. crystallite_localConstitution(g,i,e) .and. .not. crystallite_todo(g,i,e) ) then ! if updateState or updateTemperature signals broken non-local...
call crystallite_updateState(stateUpdateDone, stateConverged, g,i,e) ! update state
call crystallite_updateTemperature(temperatureUpdateDone, temperatureConverged, g,i,e) ! update temperature
crystallite_todo(g,i,e) = stateUpdateDone .and. temperatureUpdateDone
if ( (.not. stateUpdateDone .or. .not. temperatureUpdateDone) &
.and. .not. crystallite_localConstitution(g,i,e) ) then ! if updateState or updateTemperature signals broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
@ -2184,16 +2173,17 @@ endif
constitutive_dotState(g,i,e)%p = 0.0_pReal ! reset dotState to zero
enddo; enddo; enddo
!$OMP ENDDO
!$OMP END PARALLEL
! --+>> STATE LOOP <<+--
NiterationState = 0_pInt
do while (any(crystallite_todo) .and. NiterationState < nState ) ! convergence loop for crystallite
NiterationState = NiterationState + 1_pInt
!$OMP PARALLEL PRIVATE(stateConverged,temperatureConverged)
!$OMP PARALLEL
! --- STRESS INTEGRATION ---
@ -2213,13 +2203,11 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
enddo; enddo; enddo
!$OMP ENDDO
#ifndef _OPENMP
if (debug_verbosity > 5) then
!$OMP SINGLE
!$OMP CRITICAL (write2out)
write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after stress integration'
!$OMP END CRITICAL (write2out)
!$OMP END SINGLE
write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after stress integration'
endif
#endif
! --- DOT STATES ---
@ -2239,7 +2227,7 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
!$OMP SINGLE
crystallite_statedamper = 1.0_pReal
!$OMP END SINGLE
!$OMP DO PRIVATE(dot_prod12,dot_prod22)
!$OMP DO PRIVATE(dot_prod12,dot_prod22,stateUpdateDone,temperatureUpdateDone,stateConverged,temperatureConverged)
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
@ -2258,12 +2246,12 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
! --- updates ---
call crystallite_updateState(crystallite_todo(g,i,e), stateConverged, g,i,e) ! update state
!$OMP FLUSH(crystallite_todo)
call crystallite_updateTemperature(crystallite_todo(g,i,e), temperatureConverged, g,i,e) ! update temperature
!$OMP FLUSH(crystallite_todo)
call crystallite_updateState(stateUpdateDone, stateConverged, g,i,e) ! update state
call crystallite_updateTemperature(temperatureUpdateDone, temperatureConverged, g,i,e) ! update temperature
crystallite_todo(g,i,e) = stateUpdateDone .and. temperatureUpdateDone
crystallite_converged(g,i,e) = stateConverged .and. temperatureConverged
if ( .not. crystallite_localConstitution(g,i,e) .and. .not. crystallite_todo(g,i,e)) then ! if updateState or updateTemperature signals broken non-local...
if ( (.not. stateUpdateDone .or. .not. temperatureUpdateDone) &
.and. .not. crystallite_localConstitution(g,i,e) ) then ! if updateState or updateTemperature signals broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
@ -2296,13 +2284,13 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
!$OMP END PARALLEL
#ifndef _OPENMP
if (debug_verbosity > 5) then
!$OMP CRITICAL (write2out)
write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), &
' grains converged after state integration no. ', NiterationState
write(6,*)
!$OMP END CRITICAL (write2out)
write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), &
' grains converged after state integration no. ', NiterationState
write(6,*)
endif
#endif
! --- CONVERGENCE CHECK ---
@ -2314,14 +2302,14 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
endif
crystallite_todo = crystallite_todo .and. .not. crystallite_converged ! skip all converged
#ifndef _OPENMP
if (debug_verbosity > 5) then
!$OMP CRITICAL (write2out)
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. ',&
NiterationState
write(6,*)
!$OMP END CRITICAL (write2out)
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. ',&
NiterationState
write(6,*)
endif
#endif
enddo ! crystallite convergence loop
@ -2388,11 +2376,11 @@ dotState(1:mySize) = constitutive_dotState(g,i,e)%p(1:mySize) * crystallite_stat
residuum = constitutive_state(g,i,e)%p(1:mySize) - constitutive_subState0(g,i,e)%p(1:mySize) &
- dotState(1:mySize) * crystallite_subdt(g,i,e)
if (any(residuum /= residuum)) then ! if NaN occured then return without changing the state
#ifndef _OPENMP
if (debug_verbosity > 4) then
!$OMP CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState encountered NaN at el ip g ',e,i,g
!$OMP END CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState encountered NaN at el ip g ',e,i,g
endif
#endif
return
endif
@ -2404,22 +2392,22 @@ done = .true.
converged = all( abs(residuum) < constitutive_aTolState(g,i,e)%p(1:mySize) &
.or. abs(residuum) < rTol_crystalliteState * abs(state(1:mySize)) )
#ifndef _OPENMP
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
!$OMP CRITICAL (write2out)
if (converged) then
write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState converged at el ip g ',e,i,g
else
write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState did not converge at el ip g ',e,i,g
endif
write(6,*)
write(6,'(a,f6.1)') '<< CRYST >> crystallite_statedamper ',crystallite_statedamper(g,i,e)
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState',dotState(1:mySize)
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state',state(1:mySize)
write(6,*)
!$OMP END CRITICAL (write2out)
if (converged) then
write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState converged at el ip g ',e,i,g
else
write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState did not converge at el ip g ',e,i,g
endif
write(6,*)
write(6,'(a,f6.1)') '<< CRYST >> crystallite_statedamper ',crystallite_statedamper(g,i,e)
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState',dotState(1:mySize)
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state',state(1:mySize)
write(6,*)
endif
#endif
!* sync global state and dotState with local copies
@ -2471,11 +2459,11 @@ residuum = crystallite_Temperature(g,i,e) - crystallite_subTemperature0(g,i,e) &
- constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e),crystallite_Temperature(g,i,e),g,i,e) &
* crystallite_subdt(g,i,e)
if (residuum /= residuum) then
#ifndef _OPENMP
if (debug_verbosity > 4) then
!$OMP CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateTemperature encountered NaN at el ip g ',e,i,g
!$OMP END CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateTemperature encountered NaN at el ip g ',e,i,g
endif
#endif
return
endif
@ -2598,11 +2586,11 @@ integer(pLongInt) tick, &
!* be pessimistic
crystallite_integrateStress = .false.
#ifndef _OPENMP
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
!$OMP CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> integrateStress at el ip g ',e,i,g
!$OMP END CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> integrateStress at el ip g ',e,i,g
endif
#endif
!* only integrate over fraction of timestep?
@ -2628,15 +2616,15 @@ Lpguess = crystallite_Lp(1:3,1:3,g,i,e) ! ... and tak
invFp_current = math_inv3x3(Fp_current)
if (all(invFp_current == 0.0_pReal)) then ! ... failed?
#ifndef _OPENMP
if (debug_verbosity > 4) then
!$OMP CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> integrateStress failed on invFp_current inversion at el ip g ',e,i,g
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> invFp_new',math_transpose3x3(invFp_new(1:3,1:3))
endif
!$OMP END CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> integrateStress failed on invFp_current inversion at el ip g ',e,i,g
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> invFp_new',math_transpose3x3(invFp_new(1:3,1:3))
endif
endif
#endif
return
endif
A = math_mul33x33(transpose(invFp_current), math_mul33x33(transpose(Fg_new),math_mul33x33(Fg_new,invFp_current)))
@ -2662,12 +2650,12 @@ LpLoop: do
!* too many loops required ?
if (NiterationStress > nStress) then
#ifndef _OPENMP
if (debug_verbosity > 4) then
!$OMP CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> integrateStress reached loop limit at el ip g ',e,i,g
write(6,*)
!$OMP END CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> integrateStress reached loop limit at el ip g ',e,i,g
write(6,*)
endif
#endif
return
endif
@ -2700,15 +2688,15 @@ LpLoop: do
!$OMP END CRITICAL (debugTimingLpTangent)
endif
#ifndef _OPENMP
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger) &
.and. numerics_integrationMode == 1_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(a,i3)') '<< CRYST >> iteration ', NiterationStress
write(6,*)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lp_constitutive', math_transpose3x3(Lp_constitutive)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lpguess', math_transpose3x3(Lpguess)
!$OMP END CRITICAL (write2out)
write(6,'(a,i3)') '<< CRYST >> iteration ', NiterationStress
write(6,*)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lp_constitutive', math_transpose3x3(Lp_constitutive)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lpguess', math_transpose3x3(Lpguess)
endif
#endif
!* update current residuum and check for convergence of loop
@ -2728,13 +2716,13 @@ LpLoop: do
!* NaN occured at regular speed? -> return
if (any(residuum/=residuum) .and. leapfrog == 1.0) then
#ifndef _OPENMP
if (debug_verbosity > 4) then
!$OMP CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3,a,i3,a)') '<< CRYST >> integrateStress encountered NaN at el ip g ',e,i,g,&
' ; iteration ', NiterationStress,&
' >> returning..!'
!$OMP END CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3,a,i3,a)') '<< CRYST >> integrateStress encountered NaN at el ip g ',e,i,g,&
' ; iteration ', NiterationStress,&
' >> returning..!'
endif
#endif
return
@ -2746,12 +2734,12 @@ LpLoop: do
any(residuum/=residuum) & ! NaN occured
) &
) then
#ifndef _OPENMP
if (debug_verbosity > 5) then
!$OMP CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3,x,a,i3)') '<< CRYST >> integrateStress encountered high-speed crash at el ip g ',e,i,g,&
'; iteration ', NiterationStress
!$OMP END CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3,x,a,i3)') '<< CRYST >> integrateStress encountered high-speed crash at el ip g ',e,i,g,&
'; iteration ', NiterationStress
endif
#endif
maxleap = 0.5_pReal * leapfrog ! limit next acceleration
leapfrog = 1.0_pReal ! grinding halt
jacoCounter = 0_pInt ! reset counter for Jacobian update (we want to do an update next time!)
@ -2778,21 +2766,21 @@ LpLoop: do
invdRdLp = 0.0_pReal
call math_invert(9,dRdLp,invdRdLp,dummy,error) ! invert dR/dLp --> dLp/dR
if (error) then
#ifndef _OPENMP
if (debug_verbosity > 4) then
!$OMP CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLp inversion at el ip g ',e,i,g,&
' ; iteration ', NiterationStress
if (debug_verbosity > 5 &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,/,9(12(x),9(e15.3,x)/))') '<< CRYST >> dRdLp',transpose(dRdLp)
write(6,'(a,/,9(12(x),9(e15.3,x)/))') '<< CRYST >> dLpdT_constitutive',transpose(dLpdT_constitutive)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lp_constitutive',math_transpose3x3(Lp_constitutive)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lpguess',math_transpose3x3(Lpguess)
endif
!$OMP END CRITICAL (write2out)
endif
return
write(6,'(a,i5,x,i2,x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLp inversion at el ip g ',e,i,g,&
' ; iteration ', NiterationStress
if (debug_verbosity > 5 &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,/,9(12(x),9(e15.3,x)/))') '<< CRYST >> dRdLp',transpose(dRdLp)
write(6,'(a,/,9(12(x),9(e15.3,x)/))') '<< CRYST >> dLpdT_constitutive',transpose(dLpdT_constitutive)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lp_constitutive',math_transpose3x3(Lp_constitutive)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lpguess',math_transpose3x3(Lpguess)
endif
endif
#endif
return
endif
endif
jacoCounter = jacoCounter + 1_pInt ! increase counter for jaco update
@ -2822,16 +2810,16 @@ invFp_new = math_mul33x33(invFp_current,B)
invFp_new = invFp_new/math_det3x3(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize by det
call math_invert3x3(invFp_new,Fp_new,det,error)
if (error) then
#ifndef _OPENMP
if (debug_verbosity > 4) then
!$OMP CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip g ',e,i,g, &
' ; iteration ', NiterationStress
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> invFp_new',math_transpose3x3(invFp_new)
endif
!$OMP END CRITICAL (write2out)
write(6,'(a,i5,x,i2,x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip g ',e,i,g, &
' ; iteration ', NiterationStress
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> invFp_new',math_transpose3x3(invFp_new)
endif
endif
#endif
return
endif
Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resulting Fe
@ -2855,17 +2843,17 @@ crystallite_invFp(1:3,1:3,g,i,e) = invFp_new
!* set return flag to true
crystallite_integrateStress = .true.
#ifndef _OPENMP
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger) &
.and. numerics_integrationMode == 1_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> P / MPa',math_transpose3x3(crystallite_P(1:3,1:3,g,i,e))/1e6
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> Cauchy / MPa', &
math_mul33x33(crystallite_P(1:3,1:3,g,i,e), math_transpose3x3(Fg_new)) / 1e6 / math_det3x3(Fg_new)
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> Fe Lp Fe^-1', &
math_transpose3x3(math_mul33x33(Fe_new, math_mul33x33(crystallite_Lp(1:3,1:3,g,i,e), math_inv3x3(Fe_new)))) ! transpose to get correct print out order
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> Fp',math_transpose3x3(crystallite_Fp(1:3,1:3,g,i,e))
!$OMP END CRITICAL (write2out)
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> P / MPa',math_transpose3x3(crystallite_P(1:3,1:3,g,i,e))/1e6
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> Cauchy / MPa', &
math_mul33x33(crystallite_P(1:3,1:3,g,i,e), math_transpose3x3(Fg_new)) / 1e6 / math_det3x3(Fg_new)
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> Fe Lp Fe^-1', &
math_transpose3x3(math_mul33x33(Fe_new, math_mul33x33(crystallite_Lp(1:3,1:3,g,i,e), math_inv3x3(Fe_new)))) ! transpose to get correct print out order
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> Fp',math_transpose3x3(crystallite_Fp(1:3,1:3,g,i,e))
endif
#endif
if (debug_verbosity > 0) then
!$OMP CRITICAL (distributionStress)

366
code/homogenization.f90
View File

@ -53,165 +53,169 @@ CONTAINS
!* Module initialization *
!**************************************
subroutine homogenization_init(Temperature)
use prec, only: pReal,pInt
use math, only: math_I3
use debug, only: debug_verbosity
use IO, only: IO_error, IO_open_file, IO_open_jobFile
use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips
use material
use constitutive, only: constitutive_maxSizePostResults
use crystallite, only: crystallite_maxSizePostResults
use homogenization_isostrain
use homogenization_RGC
use prec, only: pReal,pInt
use math, only: math_I3
use debug, only: debug_verbosity
use IO, only: IO_error, IO_open_file, IO_open_jobFile
use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips
use material
use constitutive, only: constitutive_maxSizePostResults
use crystallite, only: crystallite_maxSizePostResults
use homogenization_isostrain
use homogenization_RGC
real(pReal) Temperature
integer(pInt), parameter :: fileunit = 200
integer(pInt) e,i,g,p,myInstance,j
integer(pInt), dimension(:,:), pointer :: thisSize
character(len=64), dimension(:,:), pointer :: thisOutput
logical knownHomogenization
implicit none
if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file
call homogenization_isostrain_init(fileunit) ! parse all homogenizations of this type
call homogenization_RGC_init(fileunit)
close(fileunit)
! write description file for homogenization output
if(.not. IO_open_jobFile(fileunit,'outputHomogenization')) call IO_error (50) ! problems in writing file
do p = 1,material_Nhomogenization
i = homogenization_typeInstance(p) ! which instance of this homogenization type
knownHomogenization = .true. ! assume valid
select case(homogenization_type(p)) ! split per homogenization type
case (homogenization_isostrain_label)
thisOutput => homogenization_isostrain_output
thisSize => homogenization_isostrain_sizePostResult
case (homogenization_RGC_label)
thisOutput => homogenization_RGC_output
thisSize => homogenization_RGC_sizePostResult
case default
knownHomogenization = .false.
end select
write(fileunit,*)
write(fileunit,'(a)') '['//trim(homogenization_name(p))//']'
write(fileunit,*)
if (knownHomogenization) then
write(fileunit,'(a)') '(type)'//char(9)//trim(homogenization_type(p))
write(fileunit,'(a,i)') '(ngrains)'//char(9),homogenization_Ngrains(p)
do e = 1,homogenization_Noutput(p)
write(fileunit,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
enddo
endif
enddo
close(fileunit)
real(pReal) Temperature
integer(pInt), parameter :: fileunit = 200
integer(pInt) e,i,g,p,myInstance,j
integer(pInt), dimension(:,:), pointer :: thisSize
character(len=64), dimension(:,:), pointer :: thisOutput
logical knownHomogenization
allocate(homogenization_state0(mesh_maxNips,mesh_NcpElems))
allocate(homogenization_subState0(mesh_maxNips,mesh_NcpElems))
allocate(homogenization_state(mesh_maxNips,mesh_NcpElems))
allocate(homogenization_sizeState(mesh_maxNips,mesh_NcpElems)); homogenization_sizeState = 0_pInt
allocate(homogenization_sizePostResults(mesh_maxNips,mesh_NcpElems)); homogenization_sizePostResults = 0_pInt
! --- PARSE HOMOGENIZATIONS FROM CONFIG FILE ---
allocate(materialpoint_dPdF(3,3,3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_dPdF = 0.0_pReal
allocate(materialpoint_F0(3,3,mesh_maxNips,mesh_NcpElems));
allocate(materialpoint_F(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_F = 0.0_pReal
allocate(materialpoint_subF0(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_subF0 = 0.0_pReal
allocate(materialpoint_subF(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_subF = 0.0_pReal
allocate(materialpoint_P(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_P = 0.0_pReal
allocate(materialpoint_Temperature(mesh_maxNips,mesh_NcpElems)); materialpoint_Temperature = Temperature
allocate(materialpoint_subFrac(mesh_maxNips,mesh_NcpElems)); materialpoint_subFrac = 0.0_pReal
allocate(materialpoint_subStep(mesh_maxNips,mesh_NcpElems)); materialpoint_subStep = 0.0_pReal
allocate(materialpoint_subdt(mesh_maxNips,mesh_NcpElems)); materialpoint_subdt = 0.0_pReal
allocate(materialpoint_requested(mesh_maxNips,mesh_NcpElems)); materialpoint_requested = .false.
allocate(materialpoint_converged(mesh_maxNips,mesh_NcpElems)); materialpoint_converged = .true.
allocate(materialpoint_doneAndHappy(2,mesh_maxNips,mesh_NcpElems)); materialpoint_doneAndHappy = .true.
forall (i = 1:mesh_maxNips,e = 1:mesh_NcpElems)
materialpoint_F0(1:3,1:3,i,e) = math_I3
materialpoint_F(1:3,1:3,i,e) = math_I3
end forall
do e = 1,mesh_NcpElems ! loop over elements
myInstance = homogenization_typeInstance(mesh_element(3,e))
do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs
select case(homogenization_type(mesh_element(3,e)))
!* isostrain
case (homogenization_isostrain_label)
if (homogenization_isostrain_sizeState(myInstance) > 0_pInt) then
allocate(homogenization_state0(i,e)%p(homogenization_isostrain_sizeState(myInstance)))
allocate(homogenization_subState0(i,e)%p(homogenization_isostrain_sizeState(myInstance)))
allocate(homogenization_state(i,e)%p(homogenization_isostrain_sizeState(myInstance)))
homogenization_state0(i,e)%p = homogenization_isostrain_stateInit(myInstance)
homogenization_sizeState(i,e) = homogenization_isostrain_sizeState(myInstance)
endif
homogenization_sizePostResults(i,e) = homogenization_isostrain_sizePostResults(myInstance)
!* RGC homogenization
case (homogenization_RGC_label)
if (homogenization_RGC_sizeState(myInstance) > 0_pInt) then
allocate(homogenization_state0(i,e)%p(homogenization_RGC_sizeState(myInstance)))
allocate(homogenization_subState0(i,e)%p(homogenization_RGC_sizeState(myInstance)))
allocate(homogenization_state(i,e)%p(homogenization_RGC_sizeState(myInstance)))
homogenization_state0(i,e)%p = homogenization_RGC_stateInit(myInstance)
homogenization_sizeState(i,e) = homogenization_RGC_sizeState(myInstance)
endif
homogenization_sizePostResults(i,e) = homogenization_RGC_sizePostResults(myInstance)
case default
call IO_error(201,ext_msg=homogenization_type(mesh_element(3,e))) ! unknown type 201 is homogenization!
end select
enddo
enddo
homogenization_maxSizeState = maxval(homogenization_sizeState)
homogenization_maxSizePostResults = maxval(homogenization_sizePostResults)
materialpoint_sizeResults = 1+ 1+homogenization_maxSizePostResults + & ! grain count, homogSize, homogResult
homogenization_maxNgrains*(1+crystallite_maxSizePostResults+ & ! results count, cryst results
1+constitutive_maxSizePostResults) ! results count, constitutive results
allocate(materialpoint_results(materialpoint_sizeResults, mesh_maxNips,mesh_NcpElems))
if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file
call homogenization_isostrain_init(fileunit)
call homogenization_RGC_init(fileunit)
close(fileunit)
! --- WRITE DESCRIPTION FILE FOR HOMOGENIZATION OUTPUT ---
if(.not. IO_open_jobFile(fileunit,'outputHomogenization')) then ! problems in writing file
call IO_error (50)
endif
do p = 1,material_Nhomogenization
i = homogenization_typeInstance(p) ! which instance of this homogenization type
knownHomogenization = .true. ! assume valid
select case(homogenization_type(p)) ! split per homogenization type
case (homogenization_isostrain_label)
thisOutput => homogenization_isostrain_output
thisSize => homogenization_isostrain_sizePostResult
case (homogenization_RGC_label)
thisOutput => homogenization_RGC_output
thisSize => homogenization_RGC_sizePostResult
case default
knownHomogenization = .false.
end select
write(fileunit,*)
write(fileunit,'(a)') '['//trim(homogenization_name(p))//']'
write(fileunit,*)
if (knownHomogenization) then
write(fileunit,'(a)') '(type)'//char(9)//trim(homogenization_type(p))
write(fileunit,'(a,i)') '(ngrains)'//char(9),homogenization_Ngrains(p)
do e = 1,homogenization_Noutput(p)
write(fileunit,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
enddo
endif
enddo
close(fileunit)
! --- ALLOCATE AND INITIALIZE GLOBAL VARIABLES ---
allocate(homogenization_state0(mesh_maxNips,mesh_NcpElems))
allocate(homogenization_subState0(mesh_maxNips,mesh_NcpElems))
allocate(homogenization_state(mesh_maxNips,mesh_NcpElems))
allocate(homogenization_sizeState(mesh_maxNips,mesh_NcpElems)); homogenization_sizeState = 0_pInt
allocate(homogenization_sizePostResults(mesh_maxNips,mesh_NcpElems)); homogenization_sizePostResults = 0_pInt
allocate(materialpoint_dPdF(3,3,3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_dPdF = 0.0_pReal
allocate(materialpoint_F0(3,3,mesh_maxNips,mesh_NcpElems));
allocate(materialpoint_F(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_F = 0.0_pReal
allocate(materialpoint_subF0(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_subF0 = 0.0_pReal
allocate(materialpoint_subF(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_subF = 0.0_pReal
allocate(materialpoint_P(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_P = 0.0_pReal
allocate(materialpoint_Temperature(mesh_maxNips,mesh_NcpElems)); materialpoint_Temperature = Temperature
allocate(materialpoint_subFrac(mesh_maxNips,mesh_NcpElems)); materialpoint_subFrac = 0.0_pReal
allocate(materialpoint_subStep(mesh_maxNips,mesh_NcpElems)); materialpoint_subStep = 0.0_pReal
allocate(materialpoint_subdt(mesh_maxNips,mesh_NcpElems)); materialpoint_subdt = 0.0_pReal
allocate(materialpoint_requested(mesh_maxNips,mesh_NcpElems)); materialpoint_requested = .false.
allocate(materialpoint_converged(mesh_maxNips,mesh_NcpElems)); materialpoint_converged = .true.
allocate(materialpoint_doneAndHappy(2,mesh_maxNips,mesh_NcpElems)); materialpoint_doneAndHappy = .true.
forall (i = 1:mesh_maxNips,e = 1:mesh_NcpElems)
materialpoint_F0(1:3,1:3,i,e) = math_I3
materialpoint_F(1:3,1:3,i,e) = math_I3
end forall
! --- ALLOCATE AND INITIALIZE GLOBAL STATE AND POSTRESULTS VARIABLES ---
!$OMP PARALLEL DO PRIVATE(myInstance)
do e = 1,mesh_NcpElems ! loop over elements
myInstance = homogenization_typeInstance(mesh_element(3,e))
do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs
select case(homogenization_type(mesh_element(3,e)))
case (homogenization_isostrain_label)
if (homogenization_isostrain_sizeState(myInstance) > 0_pInt) then
allocate(homogenization_state0(i,e)%p(homogenization_isostrain_sizeState(myInstance)))
allocate(homogenization_subState0(i,e)%p(homogenization_isostrain_sizeState(myInstance)))
allocate(homogenization_state(i,e)%p(homogenization_isostrain_sizeState(myInstance)))
homogenization_state0(i,e)%p = homogenization_isostrain_stateInit(myInstance)
homogenization_sizeState(i,e) = homogenization_isostrain_sizeState(myInstance)
endif
homogenization_sizePostResults(i,e) = homogenization_isostrain_sizePostResults(myInstance)
case (homogenization_RGC_label)
if (homogenization_RGC_sizeState(myInstance) > 0_pInt) then
allocate(homogenization_state0(i,e)%p(homogenization_RGC_sizeState(myInstance)))
allocate(homogenization_subState0(i,e)%p(homogenization_RGC_sizeState(myInstance)))
allocate(homogenization_state(i,e)%p(homogenization_RGC_sizeState(myInstance)))
homogenization_state0(i,e)%p = homogenization_RGC_stateInit(myInstance)
homogenization_sizeState(i,e) = homogenization_RGC_sizeState(myInstance)
endif
homogenization_sizePostResults(i,e) = homogenization_RGC_sizePostResults(myInstance)
case default
call IO_error(201,ext_msg=homogenization_type(mesh_element(3,e))) ! unknown type 201 is homogenization!
end select
enddo
enddo
!$OMP END PARALLEL DO
homogenization_maxSizeState = maxval(homogenization_sizeState)
homogenization_maxSizePostResults = maxval(homogenization_sizePostResults)
materialpoint_sizeResults = 1 & ! grain count
+ 1 + homogenization_maxSizePostResults & ! homogSize & homogResult
+ homogenization_maxNgrains * (1 + crystallite_maxSizePostResults & ! crystallite size & crystallite results
+ 1 + constitutive_maxSizePostResults) ! constitutive size & constitutive results
allocate(materialpoint_results(materialpoint_sizeResults,mesh_maxNips,mesh_NcpElems))
! *** Output to MARC output file ***
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) '<<<+- homogenization init -+>>>'
write(6,*) '$Id$'
write(6,*)
if (debug_verbosity > 0) then
write(6,'(a32,x,7(i5,x))') 'homogenization_state0: ', shape(homogenization_state0)
write(6,'(a32,x,7(i5,x))') 'homogenization_subState0: ', shape(homogenization_subState0)
write(6,'(a32,x,7(i5,x))') 'homogenization_state: ', shape(homogenization_state)
write(6,'(a32,x,7(i5,x))') 'homogenization_sizeState: ', shape(homogenization_sizeState)
write(6,'(a32,x,7(i5,x))') 'homogenization_sizePostResults: ', shape(homogenization_sizePostResults)
write(6,*)
write(6,'(a32,x,7(i5,x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF)
write(6,'(a32,x,7(i5,x))') 'materialpoint_F0: ', shape(materialpoint_F0)
write(6,'(a32,x,7(i5,x))') 'materialpoint_F: ', shape(materialpoint_F)
write(6,'(a32,x,7(i5,x))') 'materialpoint_subF0: ', shape(materialpoint_subF0)
write(6,'(a32,x,7(i5,x))') 'materialpoint_subF: ', shape(materialpoint_subF)
write(6,'(a32,x,7(i5,x))') 'materialpoint_P: ', shape(materialpoint_P)
write(6,'(a32,x,7(i5,x))') 'materialpoint_Temperature: ', shape(materialpoint_Temperature)
write(6,'(a32,x,7(i5,x))') 'materialpoint_subFrac: ', shape(materialpoint_subFrac)
write(6,'(a32,x,7(i5,x))') 'materialpoint_subStep: ', shape(materialpoint_subStep)
write(6,'(a32,x,7(i5,x))') 'materialpoint_subdt: ', shape(materialpoint_subdt)
write(6,'(a32,x,7(i5,x))') 'materialpoint_requested: ', shape(materialpoint_requested)
write(6,'(a32,x,7(i5,x))') 'materialpoint_converged: ', shape(materialpoint_converged)
write(6,'(a32,x,7(i5,x))') 'materialpoint_doneAndHappy: ', shape(materialpoint_doneAndHappy)
write(6,*)
write(6,'(a32,x,7(i5,x))') 'materialpoint_results: ', shape(materialpoint_results)
write(6,*)
write(6,'(a32,x,7(i5,x))') 'maxSizeState: ', homogenization_maxSizeState
write(6,'(a32,x,7(i5,x))') 'maxSizePostResults: ', homogenization_maxSizePostResults
endif
call flush(6)
write(6,*)
write(6,*) '<<<+- homogenization init -+>>>'
write(6,*) '$Id$'
write(6,*)
if (debug_verbosity > 0) then
write(6,'(a32,x,7(i5,x))') 'homogenization_state0: ', shape(homogenization_state0)
write(6,'(a32,x,7(i5,x))') 'homogenization_subState0: ', shape(homogenization_subState0)
write(6,'(a32,x,7(i5,x))') 'homogenization_state: ', shape(homogenization_state)
write(6,'(a32,x,7(i5,x))') 'homogenization_sizeState: ', shape(homogenization_sizeState)
write(6,'(a32,x,7(i5,x))') 'homogenization_sizePostResults: ', shape(homogenization_sizePostResults)
write(6,*)
write(6,'(a32,x,7(i5,x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF)
write(6,'(a32,x,7(i5,x))') 'materialpoint_F0: ', shape(materialpoint_F0)
write(6,'(a32,x,7(i5,x))') 'materialpoint_F: ', shape(materialpoint_F)
write(6,'(a32,x,7(i5,x))') 'materialpoint_subF0: ', shape(materialpoint_subF0)
write(6,'(a32,x,7(i5,x))') 'materialpoint_subF: ', shape(materialpoint_subF)
write(6,'(a32,x,7(i5,x))') 'materialpoint_P: ', shape(materialpoint_P)
write(6,'(a32,x,7(i5,x))') 'materialpoint_Temperature: ', shape(materialpoint_Temperature)
write(6,'(a32,x,7(i5,x))') 'materialpoint_subFrac: ', shape(materialpoint_subFrac)
write(6,'(a32,x,7(i5,x))') 'materialpoint_subStep: ', shape(materialpoint_subStep)
write(6,'(a32,x,7(i5,x))') 'materialpoint_subdt: ', shape(materialpoint_subdt)
write(6,'(a32,x,7(i5,x))') 'materialpoint_requested: ', shape(materialpoint_requested)
write(6,'(a32,x,7(i5,x))') 'materialpoint_converged: ', shape(materialpoint_converged)
write(6,'(a32,x,7(i5,x))') 'materialpoint_doneAndHappy: ', shape(materialpoint_doneAndHappy)
write(6,*)
write(6,'(a32,x,7(i5,x))') 'materialpoint_results: ', shape(materialpoint_results)
write(6,*)
write(6,'(a32,x,7(i5,x))') 'maxSizeState: ', homogenization_maxSizeState
write(6,'(a32,x,7(i5,x))') 'maxSizePostResults: ', homogenization_maxSizePostResults
endif
call flush(6)
!$OMP END CRITICAL (write2out)
return
endsubroutine
@ -335,13 +339,13 @@ subroutine materialpoint_stressAndItsTangent(&
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
if ( materialpoint_converged(i,e) ) then
#ifndef _OPENMP
if (debug_verbosity > 2 .and. ((e == debug_e .and. i == debug_i) .or. .not. debug_selectiveDebugger)) then
!$OMP CRITICAL (write2out)
write(6,'(a,x,f10.8,x,a,x,f10.8,x,a,/)') '<< HOMOG >> winding forward from', &
materialpoint_subFrac(i,e), 'to current materialpoint_subFrac', &
materialpoint_subFrac(i,e)+materialpoint_subStep(i,e),'in materialpoint_stressAndItsTangent'
!$OMP END CRITICAL (write2out)
write(6,'(a,x,f10.8,x,a,x,f10.8,x,a,/)') '<< HOMOG >> winding forward from', &
materialpoint_subFrac(i,e), 'to current materialpoint_subFrac', &
materialpoint_subFrac(i,e)+materialpoint_subStep(i,e),'in materialpoint_stressAndItsTangent'
endif
#endif
! calculate new subStep and new subFrac
materialpoint_subFrac(i,e) = materialpoint_subFrac(i,e) + materialpoint_subStep(i,e)
@ -386,12 +390,12 @@ subroutine materialpoint_stressAndItsTangent(&
materialpoint_subStep(i,e) = subStepSizeHomog * materialpoint_subStep(i,e) ! crystallite had severe trouble, so do a significant cutback
!$OMP FLUSH(materialpoint_subStep)
#ifndef _OPENMP
if (debug_verbosity > 2 .and. ((e == debug_e .and. i == debug_i) .or. .not. debug_selectiveDebugger)) then
!$OMP CRITICAL (write2out)
write(6,'(a,x,f10.8,/)') '<< HOMOG >> cutback step in materialpoint_stressAndItsTangent with new materialpoint_subStep:',&
materialpoint_subStep(i,e)
!$OMP END CRITICAL (write2out)
write(6,'(a,x,f10.8,/)') '<< HOMOG >> cutback step in materialpoint_stressAndItsTangent with new materialpoint_subStep:',&
materialpoint_subStep(i,e)
endif
#endif
! restore...
crystallite_Temperature(1:myNgrains,i,e) = crystallite_partionedTemperature0(1:myNgrains,i,e) ! ...temperatures
@ -496,7 +500,6 @@ subroutine materialpoint_stressAndItsTangent(&
if (.not. terminallyIll ) then
call crystallite_orientations() ! calculate crystal orientations
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
@ -505,21 +508,12 @@ subroutine materialpoint_stressAndItsTangent(&
call homogenization_averageTemperature(i,e)
enddo; enddo
!$OMP END PARALLEL DO
if (debug_verbosity > 2) then
!$OMP CRITICAL (write2out)
write (6,*)
write (6,'(a)') '<< HOMOG >> Material Point end'
write (6,*)
!$OMP END CRITICAL (write2out)
endif
else
!$OMP CRITICAL (write2out)
write (6,*)
write (6,'(a)') '<< HOMOG >> Material Point terminally ill'
write (6,*)
write (6,*)
write (6,'(a)') '<< HOMOG >> Material Point terminally ill'
write (6,*)
!$OMP END CRITICAL (write2out)
endif
return
@ -540,25 +534,31 @@ subroutine materialpoint_postResults(dt)
implicit none
real(pReal), intent(in) :: dt
integer(pInt) g,i,e,c,d,myNgrains,myCrystallite
integer(pInt) g,i,e,pos,size,myNgrains,myCrystallite
!$OMP PARALLEL DO PRIVATE(myNgrains,myCrystallite,c,d)
!$OMP PARALLEL DO PRIVATE(myNgrains,myCrystallite,pos,size)
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e))
myCrystallite = microstructure_crystallite(mesh_element(4,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
c = 0_pInt
materialpoint_results(c+1,i,e) = myNgrains; c = c+1_pInt ! tell number of grains at materialpoint
d = homogenization_sizePostResults(i,e)
materialpoint_results(c+1,i,e) = d; c = c+1_pInt ! tell size of homogenization results
if (d > 0_pInt) then ! any homogenization results to mention?
materialpoint_results(c+1:c+d,i,e) = & ! tell homogenization results
homogenization_postResults(i,e); c = c+d
pos = 0_pInt
materialpoint_results(pos+1,i,e) = myNgrains ! tell number of grains at materialpoint
pos = pos + 1_pInt
size = homogenization_sizePostResults(i,e)
materialpoint_results(pos+1,i,e) = size ! tell size of homogenization results
pos = pos + 1_pInt
if (size > 0_pInt) then ! any homogenization results to mention?
materialpoint_results(pos+1:pos+size,i,e) = homogenization_postResults(i,e) ! tell homogenization results
pos = pos + size
endif
do g = 1,myNgrains ! loop over all grains
d = 1+crystallite_sizePostResults(myCrystallite) + 1+constitutive_sizePostResults(g,i,e)
materialpoint_results(c+1:c+d,i,e) = & ! tell crystallite results
crystallite_postResults(dt,g,i,e); c = c+d
size = (1 + crystallite_sizePostResults(myCrystallite)) + (1 + constitutive_sizePostResults(g,i,e))
materialpoint_results(pos+1:pos+size,i,e) = crystallite_postResults(dt,g,i,e) ! tell crystallite results
pos = pos + size
enddo
enddo
enddo

40
code/mpie_cpfem_abaqus_exp.f
View File

@ -14,7 +14,7 @@
!
!********************************************************************
include "prec.f90" ! uses nothing else
#include "prec.f90"
MODULE mpie_interface
@ -73,25 +73,25 @@ end function
END MODULE
include "IO.f90" ! uses prec
include "numerics.f90" ! uses prec, IO
include "debug.f90" ! uses prec, numerics
include "math.f90" ! uses prec, numerics, debug
include "FEsolving.f90" ! uses prec, IO, debug
include "mesh.f90" ! uses prec, math, IO, FEsolving, debug
include "material.f90" ! uses prec, math, IO, mesh, debug
include "lattice.f90" ! uses prec, math, IO, material, debug
include "constitutive_j2.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_phenopowerlaw.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_titanmod.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_dislotwin.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_nonlocal.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug
include "crystallite.f90" ! uses prec, math, IO, numerics, Fesolving, material, mesh, constitutive, debug
include "homogenization_isostrain.f90" ! uses prec, math, IO, debug
include "homogenization_RGC.f90" ! uses prec, math, IO, numerics, mesh, debug
include "homogenization.f90" ! uses prec, math, IO, numerics, debug
include "CPFEM.f90" ! uses prec, math, IO, numerics, debug, FEsolving, mesh, lattice, constitutive, crystallite, debug
#include "IO.f90"
#include "numerics.f90"
#include "debug.f90"
#include "math.f90"
#include "FEsolving.f90"
#include "mesh.f90"
#include "material.f90"
#include "lattice.f90"
#include "constitutive_j2.f90"
#include "constitutive_phenopowerlaw.f90"
#include "constitutive_titanmod.f90"
#include "constitutive_dislotwin.f90"
#include "constitutive_nonlocal.f90"
#include "constitutive.f90"
#include "crystallite.f90"
#include "homogenization_isostrain.f90"
#include "homogenization_RGC.f90"
#include "homogenization.f90"
#include "CPFEM.f90"
subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, &
stepTime, totalTime, dt, cmname, coordMp, charLength, &

40
code/mpie_cpfem_abaqus_std.f
View File

@ -14,7 +14,7 @@
!
!********************************************************************
include "prec.f90" ! uses nothing else
#include "prec.f90"
MODULE mpie_interface
@ -75,25 +75,25 @@ end function
END MODULE
include "IO.f90" ! uses prec
include "numerics.f90" ! uses prec, IO
include "debug.f90" ! uses prec, numerics
include "math.f90" ! uses prec, numerics, debug
include "FEsolving.f90" ! uses prec, IO, debug
include "mesh.f90" ! uses prec, math, IO, FEsolving, debug
include "material.f90" ! uses prec, math, IO, mesh, debug
include "lattice.f90" ! uses prec, math, IO, material, debug
include "constitutive_j2.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_phenopowerlaw.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_titanmod.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_dislotwin.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_nonlocal.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug
include "crystallite.f90" ! uses prec, math, IO, numerics, Fesolving, material, mesh, constitutive, debug
include "homogenization_isostrain.f90" ! uses prec, math, IO, debug
include "homogenization_RGC.f90" ! uses prec, math, IO, numerics, mesh, debug
include "homogenization.f90" ! uses prec, math, IO, numerics, debug
include "CPFEM.f90" ! uses prec, math, IO, numerics, debug, FEsolving, mesh, lattice, constitutive, crystallite, debug
#include "IO.f90"
#include "numerics.f90"
#include "debug.f90"
#include "math.f90"
#include "FEsolving.f90"
#include "mesh.f90"
#include "material.f90"
#include "lattice.f90"
#include "constitutive_j2.f90"
#include "constitutive_phenopowerlaw.f90"
#include "constitutive_titanmod.f90"
#include "constitutive_dislotwin.f90"
#include "constitutive_nonlocal.f90"
#include "constitutive.f90"
#include "crystallite.f90"
#include "homogenization_isostrain.f90"
#include "homogenization_RGC.f90"
#include "homogenization.f90"
#include "CPFEM.f90"
subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,&

40
code/mpie_cpfem_marc.f90
View File

@ -36,7 +36,7 @@
! - creeps: timinc
!********************************************************************
!
include "prec.f90" ! uses nothing else
#include "prec.f90"
MODULE mpie_interface
@ -100,25 +100,25 @@ end function
END MODULE
include "IO.f90" ! uses prec
include "numerics.f90" ! uses prec, IO
include "debug.f90" ! uses prec, numerics
include "math.f90" ! uses prec, numerics, debug
include "FEsolving.f90" ! uses prec, IO, debug
include "mesh.f90" ! uses prec, math, IO, FEsolving, debug
include "material.f90" ! uses prec, math, IO, mesh, debug
include "lattice.f90" ! uses prec, math, IO, material, debug
include "constitutive_j2.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_phenopowerlaw.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_titanmod.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_dislotwin.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_nonlocal.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug
include "crystallite.f90" ! uses prec, math, IO, numerics, Fesolving, material, mesh, constitutive, debug
include "homogenization_isostrain.f90" ! uses prec, math, IO, debug
include "homogenization_RGC.f90" ! uses prec, math, IO, numerics, mesh, debug
include "homogenization.f90" ! uses prec, math, IO, numerics, debug
include "CPFEM.f90" ! uses prec, math, IO, numerics, debug, FEsolving, mesh, lattice, constitutive, crystallite, debug
#include "IO.f90"
#include "numerics.f90"
#include "debug.f90"
#include "math.f90"
#include "FEsolving.f90"
#include "mesh.f90"
#include "material.f90"
#include "lattice.f90"
#include "constitutive_j2.f90"
#include "constitutive_phenopowerlaw.f90"
#include "constitutive_titanmod.f90"
#include "constitutive_dislotwin.f90"
#include "constitutive_nonlocal.f90"
#include "constitutive.f90"
#include "crystallite.f90"
#include "homogenization_isostrain.f90"
#include "homogenization_RGC.f90"
#include "homogenization.f90"
#include "CPFEM.f90"
!********************************************************************