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

Finished applying newstate to Homogenization.

This commit is contained in:
Luv Sharma 2014-09-19 17:59:06 +00:00
parent 14530b0a76
commit 649ec1fca8
8 changed files with 87 additions and 502 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

89
code/CPFEM.f90
View File

@ -147,12 +147,11 @@ subroutine CPFEM_init
use material, only: & use material, only: &
homogenization_maxNgrains, & homogenization_maxNgrains, &
material_phase, & material_phase, &
#ifdef NEWSTATE
homogState, & homogState, &
mappingHomogenization, & mappingHomogenization, &
#endif
phase_plasticity, & phase_plasticity, &
plasticState plasticState, &
material_Nhomogenization
use crystallite, only: & use crystallite, only: &
crystallite_F0, & crystallite_F0, &
crystallite_Fp0, & crystallite_Fp0, &
@ -160,14 +159,9 @@ subroutine CPFEM_init
crystallite_dPdF0, & crystallite_dPdF0, &
crystallite_Tstar0_v, & crystallite_Tstar0_v, &
crystallite_localPlasticity crystallite_localPlasticity
use homogenization, only: &
#ifndef NEWSTATE
homogenization_state0, &
#endif
homogenization_sizeState
implicit none implicit none
integer(pInt) :: i,j,k,l,m,ph integer(pInt) :: i,j,k,l,m,ph,homog
write(6,'(/,a)') ' <<<+- CPFEM init -+>>>' write(6,'(/,a)') ' <<<+- CPFEM init -+>>>'
write(6,'(a)') ' $Id$' write(6,'(a)') ' $Id$'
@ -214,29 +208,26 @@ subroutine CPFEM_init
call IO_read_realFile(777,'convergedStateConst',modelName) call IO_read_realFile(777,'convergedStateConst',modelName)
m = 0_pInt m = 0_pInt
readInstances: do ph = 1_pInt, size(phase_plasticity) readPlasticityInstances: do ph = 1_pInt, size(phase_plasticity)
do k = 1_pInt, plasticState(ph)%sizeState do k = 1_pInt, plasticState(ph)%sizeState
do l = 1, size(plasticState(ph)%state0(1,:)) do l = 1, size(plasticState(ph)%state0(1,:))
m = m+1_pInt m = m+1_pInt
read(777,rec=m) plasticState(ph)%state0(k,l) read(777,rec=m) plasticState(ph)%state0(k,l)
enddo; enddo enddo; enddo
enddo readInstances enddo readPlasticityInstances
close (777) close (777)
call IO_read_realFile(777,'convergedStateHomog',modelName) call IO_read_realFile(777,'convergedStateHomog',modelName)
m = 0_pInt m = 0_pInt
do k = 1,mesh_NcpElems; do j = 1,mesh_maxNips readHomogInstances: do homog = 1_pInt, material_Nhomogenization
do l = 1,homogenization_sizeState(j,k) do k = 1_pInt, homogState(homog)%sizeState
m = m+1_pInt do l = 1, size(homogState(homog)%state0(1,:))
#ifdef NEWSTATE m = m+1_pInt
read(777,rec=m) homogState(mappingHomogenization(2,j,k))%state0(l,mappingHomogenization(1,j,k)) read(777,rec=m) homogState(homog)%state0(k,l)
#else enddo; enddo
read(777,rec=m) homogenization_state0(j,k)%p(l) enddo readHomogInstances
#endif
enddo
enddo; enddo
close (777) close (777)
#if defined(Marc4DAMASK) || defined(Abaqus) #if defined(Marc4DAMASK) || defined(Abaqus)
call IO_read_realFile(777,'convergeddcsdE',modelName,size(CPFEM_dcsdE)) call IO_read_realFile(777,'convergeddcsdE',modelName,size(CPFEM_dcsdE))
read (777,rec=1) CPFEM_dcsdE read (777,rec=1) CPFEM_dcsdE
@ -314,14 +305,12 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip)
microstructure_elemhomo, & microstructure_elemhomo, &
plasticState, & plasticState, &
damageState, & damageState, &
#ifdef NEWSTATE
homogState, & homogState, &
mappingHomogenization, &
#endif
thermalState, & thermalState, &
mappingConstitutive, & mappingConstitutive, &
material_phase, & material_phase, &
phase_plasticity phase_plasticity, &
material_Nhomogenization
use crystallite, only: & use crystallite, only: &
crystallite_partionedF,& crystallite_partionedF,&
crystallite_F0, & crystallite_F0, &
@ -334,13 +323,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip)
crystallite_Tstar0_v, & crystallite_Tstar0_v, &
crystallite_Tstar_v, & crystallite_Tstar_v, &
crystallite_temperature crystallite_temperature
use homogenization, only: & use homogenization, only: & materialpoint_F, &
homogenization_sizeState, &
#ifndef NEWSTATE
homogenization_state, &
homogenization_state0, &
#endif
materialpoint_F, &
materialpoint_F0, & materialpoint_F0, &
materialpoint_P, & materialpoint_P, &
materialpoint_dPdF, & materialpoint_dPdF, &
@ -378,7 +361,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip)
#endif #endif
integer(pInt) elCP, & ! crystal plasticity element number integer(pInt) elCP, & ! crystal plasticity element number
i, j, k, l, m, n, ph i, j, k, l, m, n, ph, homog
logical updateJaco ! flag indicating if JAcobian has to be updated logical updateJaco ! flag indicating if JAcobian has to be updated
#if defined(Marc4DAMASK) || defined(Abaqus) #if defined(Marc4DAMASK) || defined(Abaqus)
@ -433,19 +416,10 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip)
endif endif
endif endif
!$OMP PARALLEL DO do homog = 1_pInt, material_Nhomogenization
do k = 1,mesh_NcpElems homogState(homog)%state0 = homogState(homog)%state
do j = 1,mesh_maxNips enddo
if (homogenization_sizeState(j,k) > 0_pInt) &
#ifdef NEWSTATE
homogState(mappingHomogenization(2,j,k))%state0(:,mappingHomogenization(1,j,k)) = &
homogState(mappingHomogenization(2,j,k))%state(:,mappingHomogenization(1,j,k)) ! internal state of homogenization scheme
#else
homogenization_state0(j,k)%p = homogenization_state(j,k)%p ! internal state of homogenization scheme
#endif
enddo
enddo
!$OMP END PARALLEL DO
! * dump the last converged values of each essential variable to a binary file ! * dump the last converged values of each essential variable to a binary file
@ -479,27 +453,24 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip)
call IO_write_jobRealFile(777,'convergedStateConst') call IO_write_jobRealFile(777,'convergedStateConst')
m = 0_pInt m = 0_pInt
writeInstances: do ph = 1_pInt, size(phase_plasticity) writePlasticityInstances: do ph = 1_pInt, size(phase_plasticity)
do k = 1_pInt, plasticState(ph)%sizeState do k = 1_pInt, plasticState(ph)%sizeState
do l = 1, size(plasticState(ph)%state0(1,:)) do l = 1, size(plasticState(ph)%state0(1,:))
m = m+1_pInt m = m+1_pInt
write(777,rec=m) plasticState(ph)%state0(k,l) write(777,rec=m) plasticState(ph)%state0(k,l)
enddo; enddo enddo; enddo
enddo writeInstances enddo writePlasticityInstances
close (777) close (777)
call IO_write_jobRealFile(777,'convergedStateHomog') call IO_write_jobRealFile(777,'convergedStateHomog')
m = 0_pInt m = 0_pInt
do k = 1,mesh_NcpElems; do j = 1,mesh_maxNips writeHomogInstances: do homog = 1_pInt, material_Nhomogenization
do l = 1,homogenization_sizeState(j,k) do k = 1_pInt, homogState(homog)%sizeState
m = m+1_pInt do l = 1, size(homogState(homog)%state0(1,:))
#ifdef NEWSTATE m = m+1_pInt
write(777,rec=m) homogState(mappingHomogenization(2,j,k))%state0(l,mappingHomogenization(1,j,k)) write(777,rec=m) homogState(homog)%state0(k,l)
#else enddo; enddo
write(777,rec=m) homogenization_state0(j,k)%p(l) enddo writeHomogInstances
#endif
enddo
enddo; enddo
close (777) close (777)
#if defined(Marc4DAMASK) || defined(Abaqus) #if defined(Marc4DAMASK) || defined(Abaqus)

1
code/Makefile
View File

@ -458,6 +458,7 @@ FEM_utilities.o: FEM_utilities.f90 \
FEZoo.o: $(wildcard FEZoo.f90) \ FEZoo.o: $(wildcard FEZoo.f90) \
IO.o IO.o
$(IGNORE) $(PREFIX) $(COMPILERNAME) $(COMPILE) -c ../private/FEM/code/FEZoo.f90 $(SUFFIX) $(IGNORE) $(PREFIX) $(COMPILERNAME) $(COMPILE) -c ../private/FEM/code/FEZoo.f90 $(SUFFIX)
touch FEZoo.o
CPFEM.o: CPFEM.f90\ CPFEM.o: CPFEM.f90\
homogenization.o homogenization.o

7
code/constitutive.f90
View File

@ -215,17 +215,10 @@ subroutine constitutive_init
#endif #endif
#ifdef TODO #ifdef TODO
!--------------------------------------------------------------------------------------------------
! write out state size file
call IO_write_jobIntFile(777,'sizeStateConst', size(constitutive_sizeState))
write (777,rec=1) constitutive_sizeState
close(777)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report ! report
constitutive_maxSizeState = maxval(constitutive_sizeState) constitutive_maxSizeState = maxval(constitutive_sizeState)
constitutive_maxSizeDotState = maxval(constitutive_sizeDotState) constitutive_maxSizeDotState = maxval(constitutive_sizeDotState)
constitutive_maxSizePostResults = maxval(constitutive_sizePostResults)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) then
write(6,'(a32,1x,7(i8,1x))') 'constitutive_state0: ', shape(constitutive_state0) write(6,'(a32,1x,7(i8,1x))') 'constitutive_state0: ', shape(constitutive_state0)

130
code/homogenization.f90
View File

@ -9,17 +9,12 @@
module homogenization module homogenization
use prec, only: & use prec, only: &
pInt, & pInt, &
pReal, & pReal
p_vec
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! General variables for the homogenization at a material point ! General variables for the homogenization at a material point
implicit none implicit none
private private
#ifndef NEWSTATE
type(p_vec), dimension(:,:), allocatable, public :: &
homogenization_state0 !< pointer array to homogenization state at start of FE increment
#endif
real(pReal), dimension(:,:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:,:), allocatable, public :: &
materialpoint_F0, & !< def grad of IP at start of FE increment materialpoint_F0, & !< def grad of IP at start of FE increment
materialpoint_F, & !< def grad of IP to be reached at end of FE increment materialpoint_F, & !< def grad of IP to be reached at end of FE increment
@ -28,20 +23,12 @@ module homogenization
materialpoint_dPdF !< tangent of first P--K stress at IP materialpoint_dPdF !< tangent of first P--K stress at IP
real(pReal), dimension(:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:), allocatable, public :: &
materialpoint_results !< results array of material point materialpoint_results !< results array of material point
#ifndef NEWSTATE
type(p_vec), dimension(:,:), allocatable, public, protected :: &
homogenization_state !< pointer array to current homogenization state (end of converged time step)
#endif
integer(pInt), dimension(:,:), allocatable, public, protected :: &
homogenization_sizeState !< size of state array per grain
integer(pInt), public, protected :: & integer(pInt), public, protected :: &
materialpoint_sizeResults, & materialpoint_sizeResults, &
homogenization_maxSizePostResults homogenization_maxSizePostResults
real(pReal), dimension(:,:), allocatable, public, protected :: & real(pReal), dimension(:,:), allocatable, public, protected :: &
materialpoint_heat materialpoint_heat
type(p_vec), dimension(:,:), allocatable, private :: &
homogenization_subState0 !< pointer array to homogenization state at start of homogenization increment
real(pReal), dimension(:,:,:,:), allocatable, private :: & real(pReal), dimension(:,:,:,:), allocatable, private :: &
materialpoint_subF0, & !< def grad of IP at beginning of homogenization increment materialpoint_subF0, & !< def grad of IP at beginning of homogenization increment
materialpoint_subF !< def grad of IP to be reached at end of homog inc materialpoint_subF !< def grad of IP to be reached at end of homog inc
@ -233,12 +220,6 @@ subroutine homogenization_init()
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate and initialize global variables ! allocate and initialize global variables
#ifndef NEWSTATE
allocate(homogenization_state0(mesh_maxNips,mesh_NcpElems))
allocate(homogenization_subState0(mesh_maxNips,mesh_NcpElems))
allocate(homogenization_state(mesh_maxNips,mesh_NcpElems))
#endif
allocate(homogenization_sizeState(mesh_maxNips,mesh_NcpElems), source=0_pInt)
allocate(homogenization_sizePostResults(mesh_maxNips,mesh_NcpElems), source=0_pInt) allocate(homogenization_sizePostResults(mesh_maxNips,mesh_NcpElems), source=0_pInt)
allocate(materialpoint_heat(mesh_maxNips,mesh_NcpElems), source=0.0_pReal) allocate(materialpoint_heat(mesh_maxNips,mesh_NcpElems), source=0.0_pReal)
allocate(materialpoint_dPdF(3,3,3,3,mesh_maxNips,mesh_NcpElems), source=0.0_pReal) allocate(materialpoint_dPdF(3,3,3,3,mesh_maxNips,mesh_NcpElems), source=0.0_pReal)
@ -265,50 +246,13 @@ subroutine homogenization_init()
InstancePosition(myInstance) = InstancePosition(myInstance)+1_pInt InstancePosition(myInstance) = InstancePosition(myInstance)+1_pInt
mapping(e,1:4) = [instancePosition(myinstance),myinstance,e,i] mapping(e,1:4) = [instancePosition(myinstance),myinstance,e,i]
#endif #endif
select case(homogenization_type(mesh_element(3,e))) homogenization_sizePostResults(i,e) = homogState(mappingHomogenization(2,i,e))%sizePostResults
case (HOMOGENIZATION_none_ID)
#ifdef NEWSTATE
homogenization_sizePostResults(i,e) = homogState(mappingHomogenization(2,i,e))%sizePostResults
#else
homogenization_sizePostResults(i,e) = 0_pInt
#endif
case (HOMOGENIZATION_ISOSTRAIN_ID)
#ifdef NEWSTATE
homogenization_sizePostResults(i,e) = homogState(mappingHomogenization(2,i,e))%sizePostResults
#else
homogenization_sizePostResults(i,e) = homogenization_isostrain_sizePostResults(myInstance)
#endif
case (HOMOGENIZATION_RGC_ID)
if (homogenization_RGC_sizeState(myInstance) > 0_pInt) then
#ifdef NEWSTATE
homogenization_sizeState(i,e) = homogState(mappingHomogenization(2,i,e))%sizeState
#else
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 = 0.0_pReal
homogenization_sizeState(i,e) = homogenization_RGC_sizeState(myInstance)
#endif
endif
#ifdef NEWSTATE
homogenization_sizePostResults(i,e) = homogState(mappingHomogenization(2,i,e))%sizePostResults
#else
homogenization_sizePostResults(i,e) = homogenization_RGC_sizePostResults(myInstance)
#endif
end select
enddo IpLooping enddo IpLooping
enddo elementLooping enddo elementLooping
#ifdef HDF #ifdef HDF
call HDF5_mappingHomogenization(mapping) call HDF5_mappingHomogenization(mapping)
#endif #endif
!--------------------------------------------------------------------------------------------------
! write state size file out
call IO_write_jobIntFile(777,'sizeStateHomog',size(homogenization_sizeState))
write (777,rec=1) homogenization_sizeState
close(777)
homogenization_maxSizeState = maxval(homogenization_sizeState)
homogenization_maxSizePostResults = maxval(homogenization_sizePostResults) homogenization_maxSizePostResults = maxval(homogenization_sizePostResults)
materialpoint_sizeResults = 1 & ! grain count materialpoint_sizeResults = 1 & ! grain count
+ 1 + homogenization_maxSizePostResults & ! homogSize & homogResult + 1 + homogenization_maxSizePostResults & ! homogSize & homogResult
@ -323,12 +267,11 @@ subroutine homogenization_init()
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then
#ifndef NEWSTATE #ifdef TODO
write(6,'(a32,1x,7(i8,1x))') 'homogenization_state0: ', shape(homogenization_state0) write(6,'(a32,1x,7(i8,1x))') 'homogenization_state0: ', shape(homogenization_state0)
write(6,'(a32,1x,7(i8,1x))') 'homogenization_subState0: ', shape(homogenization_subState0) write(6,'(a32,1x,7(i8,1x))') 'homogenization_subState0: ', shape(homogenization_subState0)
write(6,'(a32,1x,7(i8,1x))') 'homogenization_state: ', shape(homogenization_state) write(6,'(a32,1x,7(i8,1x))') 'homogenization_state: ', shape(homogenization_state)
#endif #endif
write(6,'(a32,1x,7(i8,1x))') 'homogenization_sizeState: ', shape(homogenization_sizeState)
write(6,'(a32,1x,7(i8,1x),/)') 'homogenization_sizePostResults: ', shape(homogenization_sizePostResults) write(6,'(a32,1x,7(i8,1x),/)') 'homogenization_sizePostResults: ', shape(homogenization_sizePostResults)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F0: ', shape(materialpoint_F0) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F0: ', shape(materialpoint_F0)
@ -344,7 +287,6 @@ subroutine homogenization_init()
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_converged: ', shape(materialpoint_converged) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_converged: ', shape(materialpoint_converged)
write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_doneAndHappy: ', shape(materialpoint_doneAndHappy) write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_doneAndHappy: ', shape(materialpoint_doneAndHappy)
write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_results: ', shape(materialpoint_results) write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_results: ', shape(materialpoint_results)
write(6,'(a32,1x,7(i8,1x))') 'maxSizeState: ', homogenization_maxSizeState
write(6,'(a32,1x,7(i8,1x))') 'maxSizePostResults: ', homogenization_maxSizePostResults write(6,'(a32,1x,7(i8,1x))') 'maxSizePostResults: ', homogenization_maxSizePostResults
endif endif
flush(6) flush(6)
@ -379,10 +321,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
plasticState, & plasticState, &
damageState, & damageState, &
thermalState, & thermalState, &
#ifdef NEWSTATE
homogState, & homogState, &
mappingHomogenization, & mappingHomogenization, &
#endif
mappingConstitutive, & mappingConstitutive, &
homogenization_Ngrains homogenization_Ngrains
@ -472,14 +412,10 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
materialpoint_converged(i,e) = .false. ! pretend failed step of twice the required size materialpoint_converged(i,e) = .false. ! pretend failed step of twice the required size
materialpoint_requested(i,e) = .true. ! everybody requires calculation materialpoint_requested(i,e) = .true. ! everybody requires calculation
endforall endforall
#ifdef NEWSTATE forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), &
forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), homogenization_sizeState(i,e) > 0_pInt) & homogState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) &
homogState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = & homogState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = &
homogState(mappingHomogenization(2,i,e))%State0(:,mappingHomogenization(1,i,e)) ! ...internal homogenization state homogState(mappingHomogenization(2,i,e))%State0( :,mappingHomogenization(1,i,e)) ! ...internal homogenization state
#else
forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), homogenization_sizeState(i,e) > 0_pInt) &
homogenization_subState0(i,e)%p = homogenization_state0(i,e)%p ! ...internal homogenization state
#endif
enddo enddo
NiterationHomog = 0_pInt NiterationHomog = 0_pInt
@ -526,13 +462,9 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
thermalState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) = & thermalState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) = &
thermalState(mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e)) thermalState(mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e))
end forall end forall
if (homogenization_sizeState(i,e) > 0_pInt) & if (homogState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) &
#ifdef NEWSTATE
homogState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = & homogState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = &
homogState(mappingHomogenization(2,i,e))%state( :,mappingHomogenization(1,i,e)) homogState(mappingHomogenization(2,i,e))%state( :,mappingHomogenization(1,i,e))
#else
homogenization_subState0(i,e)%p = homogenization_state(i,e)%p ! ...internal state of homog scheme
#endif
materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e) ! ...def grad materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e) ! ...def grad
!$OMP FLUSH(materialpoint_subF0) !$OMP FLUSH(materialpoint_subF0)
elseif (materialpoint_requested(i,e)) then steppingNeeded ! already at final time (??) elseif (materialpoint_requested(i,e)) then steppingNeeded ! already at final time (??)
@ -585,13 +517,9 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
thermalState(mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e)) = & thermalState(mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e)) = &
thermalState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) thermalState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e))
end forall end forall
if (homogenization_sizeState(i,e) > 0_pInt) & if (homogState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) &
#ifdef NEWSTATE
homogState(mappingHomogenization(2,i,e))%state( :,mappingHomogenization(1,i,e)) = & homogState(mappingHomogenization(2,i,e))%state( :,mappingHomogenization(1,i,e)) = &
homogState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) homogState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e))
#else
homogenization_state(i,e)%p = homogenization_subState0(i,e)%p ! ...internal state of homog scheme
#endif
endif endif
endif converged endif converged
@ -803,21 +731,11 @@ subroutine homogenization_partitionDeformation(ip,el)
materialpoint_subF(1:3,1:3,ip,el),& materialpoint_subF(1:3,1:3,ip,el),&
el) el)
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
#ifdef NEWSTATE
call homogenization_RGC_partitionDeformation(& call homogenization_RGC_partitionDeformation(&
crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), & crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), &
materialpoint_subF(1:3,1:3,ip,el),& materialpoint_subF(1:3,1:3,ip,el),&
ip, & ip, &
el) el)
#else
call homogenization_RGC_partitionDeformation(&
crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), &
materialpoint_subF(1:3,1:3,ip,el),&
homogenization_state(ip,el), &
ip, &
el)
#endif
end select chosenHomogenization end select chosenHomogenization
end subroutine homogenization_partitionDeformation end subroutine homogenization_partitionDeformation
@ -852,7 +770,7 @@ function homogenization_updateState(ip,el)
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
homogenization_updateState = & homogenization_updateState = &
#ifdef NEWSTATE
homogenization_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), & homogenization_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), &
crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), & crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), &
crystallite_partionedF0(1:3,1:3,1:homogenization_maxNgrains,ip,el),& crystallite_partionedF0(1:3,1:3,1:homogenization_maxNgrains,ip,el),&
@ -861,18 +779,6 @@ function homogenization_updateState(ip,el)
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), & crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), &
ip, & ip, &
el) el)
#else
homogenization_RGC_updateState(homogenization_state(ip,el), &
homogenization_subState0(ip,el), &
crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), &
crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), &
crystallite_partionedF0(1:3,1:3,1:homogenization_maxNgrains,ip,el),&
materialpoint_subF(1:3,1:3,ip,el),&
materialpoint_subdt(ip,el), &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), &
ip, &
el)
#endif
case default chosenHomogenization case default chosenHomogenization
homogenization_updateState = .true. homogenization_updateState = .true.
end select chosenHomogenization end select chosenHomogenization
@ -1027,11 +933,11 @@ function field_getMassDensity(ip,el)
select case(field_thermal_type(material_homog(ip,el))) select case(field_thermal_type(material_homog(ip,el)))
case (FIELD_THERMAL_ADIABATIC_ID) case (FIELD_THERMAL_ADIABATIC_ID)
field_getMassDensity = 0.0_pReal field_getMassDensity = 0.0_pReal
case (FIELD_THERMAL_CONDUCTION_ID) case (FIELD_THERMAL_CONDUCTION_ID)
do ipc = 1, homogenization_Ngrains(mesh_element(3,el)) do ipc = 1, homogenization_Ngrains(mesh_element(3,el))
field_getMassDensity = field_getMassDensity + lattice_massDensity(material_phase(ipc,ip,el)) field_getMassDensity = field_getMassDensity + lattice_massDensity(material_phase(ipc,ip,el))
enddo enddo
end select end select
@ -1215,6 +1121,7 @@ end function field_getDAMAGE
subroutine field_putDAMAGE(ip,el,fieldDamageValue) ! naming scheme subroutine field_putDAMAGE(ip,el,fieldDamageValue) ! naming scheme
use material, only: & use material, only: &
fieldDamage, & fieldDamage, &
material_homog, &
mappingHomogenization, & mappingHomogenization, &
field_damage_type, & field_damage_type, &
FIELD_DAMAGE_NONLOCAL_ID FIELD_DAMAGE_NONLOCAL_ID
@ -1338,22 +1245,11 @@ function homogenization_postResults(ip,el)
materialpoint_P(1:3,1:3,ip,el), & materialpoint_P(1:3,1:3,ip,el), &
materialpoint_F(1:3,1:3,ip,el)) materialpoint_F(1:3,1:3,ip,el))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
#ifdef NEWSTATE
homogenization_postResults = homogenization_RGC_postResults(& homogenization_postResults = homogenization_RGC_postResults(&
ip, & ip, &
el, & el, &
materialpoint_P(1:3,1:3,ip,el), & materialpoint_P(1:3,1:3,ip,el), &
materialpoint_F(1:3,1:3,ip,el)) materialpoint_F(1:3,1:3,ip,el))
#else
homogenization_postResults = homogenization_RGC_postResults(&
homogenization_state(ip,el),&
ip, &
el, &
materialpoint_P(1:3,1:3,ip,el), &
materialpoint_F(1:3,1:3,ip,el))
#endif
end select chosenHomogenization end select chosenHomogenization
end function homogenization_postResults end function homogenization_postResults

235
code/homogenization_RGC.f90
View File

@ -21,10 +21,8 @@ module homogenization_RGC
homogenization_RGC_sizePostResult homogenization_RGC_sizePostResult
character(len=64), dimension(:,:), allocatable,target, public :: & character(len=64), dimension(:,:), allocatable,target, public :: &
homogenization_RGC_output ! name of each post result output homogenization_RGC_output ! name of each post result output
#ifdef NEWSTATE
integer(pInt), dimension(:), allocatable, private :: & integer(pInt), dimension(:), allocatable, private :: &
homogenization_RGC_Noutput !< number of outputs per homog instance homogenization_RGC_Noutput !< number of outputs per homog instance
#endif
integer(pInt), dimension(:,:), allocatable, private :: & integer(pInt), dimension(:,:), allocatable, private :: &
homogenization_RGC_Ngrains homogenization_RGC_Ngrains
real(pReal), dimension(:,:), allocatable, private :: & real(pReal), dimension(:,:), allocatable, private :: &
@ -49,7 +47,7 @@ module homogenization_RGC
avgfirstpiola_ID avgfirstpiola_ID
end enum end enum
integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: &
homogenization_RGC_outputID !< ID of each post result output homogenization_RGC_outputID !< ID of each post result output
public :: & public :: &
homogenization_RGC_init, & homogenization_RGC_init, &
@ -103,17 +101,15 @@ subroutine homogenization_RGC_init(fileUnit)
use material use material
implicit none implicit none
integer(pInt), intent(in) :: fileUnit !< file pointer to material configuration integer(pInt), intent(in) :: fileUnit !< file pointer to material configuration
integer(pInt), parameter :: MAXNCHUNKS = 4_pInt integer(pInt), parameter :: MAXNCHUNKS = 4_pInt
integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions
#ifdef NEWSTATE
integer :: & integer :: &
homog, & homog, &
NofMyHomog, & NofMyHomog, &
o, & o, &
instance, & instance, &
sizeHState sizeHState
#endif
integer(pInt) :: section=0_pInt, maxNinstance, i,j,e, output=-1_pInt, mySize, myInstance integer(pInt) :: section=0_pInt, maxNinstance, i,j,e, output=-1_pInt, mySize, myInstance
character(len=65536) :: & character(len=65536) :: &
tag = '', & tag = '', &
@ -126,15 +122,11 @@ subroutine homogenization_RGC_init(fileUnit)
maxNinstance = int(count(homogenization_type == HOMOGENIZATION_RGC_ID),pInt) maxNinstance = int(count(homogenization_type == HOMOGENIZATION_RGC_ID),pInt)
if (maxNinstance == 0_pInt) return if (maxNinstance == 0_pInt) return
#ifdef NEWSTATE
if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
#endif
allocate(homogenization_RGC_sizeState(maxNinstance), source=0_pInt) allocate(homogenization_RGC_sizeState(maxNinstance), source=0_pInt)
allocate(homogenization_RGC_sizePostResults(maxNinstance), source=0_pInt) allocate(homogenization_RGC_sizePostResults(maxNinstance), source=0_pInt)
#ifdef NEWSTATE
allocate(homogenization_RGC_Noutput(maxNinstance), source=0_pInt) allocate(homogenization_RGC_Noutput(maxNinstance), source=0_pInt)
#endif
allocate(homogenization_RGC_Ngrains(3,maxNinstance), source=0_pInt) allocate(homogenization_RGC_Ngrains(3,maxNinstance), source=0_pInt)
allocate(homogenization_RGC_ciAlpha(maxNinstance), source=0.0_pReal) allocate(homogenization_RGC_ciAlpha(maxNinstance), source=0.0_pReal)
allocate(homogenization_RGC_xiAlpha(maxNinstance), source=0.0_pReal) allocate(homogenization_RGC_xiAlpha(maxNinstance), source=0.0_pReal)
@ -177,54 +169,34 @@ subroutine homogenization_RGC_init(fileUnit)
homogenization_RGC_output(output,i) = IO_lc(IO_stringValue(line,positions,2_pInt)) homogenization_RGC_output(output,i) = IO_lc(IO_stringValue(line,positions,2_pInt))
select case(homogenization_RGC_output(output,i)) select case(homogenization_RGC_output(output,i))
case('constitutivework') case('constitutivework')
#ifdef NEWSTATE
homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) + 1_pInt homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) + 1_pInt
#endif
homogenization_RGC_outputID(output,i) = constitutivework_ID homogenization_RGC_outputID(output,i) = constitutivework_ID
case('penaltyenergy') case('penaltyenergy')
#ifdef NEWSTATE
homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) + 1_pInt homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) + 1_pInt
#endif
homogenization_RGC_outputID(output,i) = penaltyenergy_ID homogenization_RGC_outputID(output,i) = penaltyenergy_ID
case('volumediscrepancy') case('volumediscrepancy')
#ifdef NEWSTATE
homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) + 1_pInt homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) + 1_pInt
#endif
homogenization_RGC_outputID(output,i) = volumediscrepancy_ID homogenization_RGC_outputID(output,i) = volumediscrepancy_ID
case('averagerelaxrate') case('averagerelaxrate')
#ifdef NEWSTATE
homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) + 1_pInt homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) + 1_pInt
#endif
homogenization_RGC_outputID(output,i) = averagerelaxrate_ID homogenization_RGC_outputID(output,i) = averagerelaxrate_ID
case('maximumrelaxrate') case('maximumrelaxrate')
#ifdef NEWSTATE
homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) + 1_pInt homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) + 1_pInt
#endif
homogenization_RGC_outputID(output,i) = maximumrelaxrate_ID homogenization_RGC_outputID(output,i) = maximumrelaxrate_ID
case('magnitudemismatch') case('magnitudemismatch')
#ifdef NEWSTATE
homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) + 1_pInt homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) + 1_pInt
#endif
homogenization_RGC_outputID(output,i) = magnitudemismatch_ID homogenization_RGC_outputID(output,i) = magnitudemismatch_ID
case('temperature') case('temperature')
#ifdef NEWSTATE
homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) + 1_pInt homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) + 1_pInt
#endif
homogenization_RGC_outputID(output,i) = temperature_ID homogenization_RGC_outputID(output,i) = temperature_ID
case('ipcoords') case('ipcoords')
#ifdef NEWSTATE
homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) + 1_pInt homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) + 1_pInt
#endif
homogenization_RGC_outputID(output,i) = ipcoords_ID homogenization_RGC_outputID(output,i) = ipcoords_ID
case('avgdefgrad','avgf') case('avgdefgrad','avgf')
#ifdef NEWSTATE
homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) + 1_pInt homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) + 1_pInt
#endif
homogenization_RGC_outputID(output,i) = avgdefgrad_ID homogenization_RGC_outputID(output,i) = avgdefgrad_ID
case('avgp','avgfirstpiola','avg1stpiola') case('avgp','avgfirstpiola','avg1stpiola')
#ifdef NEWSTATE
homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) + 1_pInt homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) + 1_pInt
#endif
homogenization_RGC_outputID(output,i) = avgfirstpiola_ID homogenization_RGC_outputID(output,i) = avgfirstpiola_ID
case default case default
call IO_error(105_pInt,ext_msg=IO_stringValue(line,positions,2_pInt)//& call IO_error(105_pInt,ext_msg=IO_stringValue(line,positions,2_pInt)//&
@ -289,7 +261,6 @@ subroutine homogenization_RGC_init(fileUnit)
enddo enddo
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
#ifdef NEWSTATE
initializeInstances: do homog = 1_pInt, material_Nhomogenization initializeInstances: do homog = 1_pInt, material_Nhomogenization
myHomog: if (homogenization_type(homog) == HOMOGENIZATION_RGC_ID) then myHomog: if (homogenization_type(homog) == HOMOGENIZATION_RGC_ID) then
NofMyHomog = count(material_homog == homog) NofMyHomog = count(material_homog == homog)
@ -329,44 +300,14 @@ subroutine homogenization_RGC_init(fileUnit)
! allocate state arrays ! allocate state arrays
homogState(homog)%sizeState = sizeHState homogState(homog)%sizeState = sizeHState
homogState(homog)%sizePostResults = homogenization_RGC_sizePostResults(instance) homogState(homog)%sizePostResults = homogenization_RGC_sizePostResults(instance)
allocate(homogState(homog)%state0 ( sizeHState,NofMyHomog), source=0.0_pReal) allocate(homogState(homog)%state0 (sizeHState,NofMyHomog), source=0.0_pReal)
allocate(homogState(homog)%subState0 ( sizeHState,NofMyHomog), source=0.0_pReal) allocate(homogState(homog)%subState0(sizeHState,NofMyHomog), source=0.0_pReal)
allocate(homogState(homog)%state ( sizeHState,NofMyHomog), source=0.0_pReal) allocate(homogState(homog)%state (sizeHState,NofMyHomog), source=0.0_pReal)
endif myHomog endif myHomog
enddo initializeInstances enddo initializeInstances
#else
do i = 1_pInt,maxNinstance
do j = 1_pInt,maxval(homogenization_Noutput)
select case(homogenization_RGC_outputID(j,i))
case(temperature_ID,constitutivework_ID,penaltyenergy_ID,volumediscrepancy_ID, &
averagerelaxrate_ID,maximumrelaxrate_ID)
mySize = 1_pInt
case(ipcoords_ID,magnitudemismatch_ID)
mySize = 3_pInt
case(avgdefgrad_ID,avgfirstpiola_ID)
mySize = 9_pInt
case default
mySize = 0_pInt
end select
outputFound: if (mySize > 0_pInt) then
homogenization_RGC_sizePostResult(j,i) = mySize
homogenization_RGC_sizePostResults(i) = &
homogenization_RGC_sizePostResults(i) + mySize
endif outputFound
enddo
homogenization_RGC_sizeState(i) &
= 3_pInt*(homogenization_RGC_Ngrains(1,i)-1_pInt)*homogenization_RGC_Ngrains(2,i)*homogenization_RGC_Ngrains(3,i) &
+ 3_pInt*homogenization_RGC_Ngrains(1,i)*(homogenization_RGC_Ngrains(2,i)-1_pInt)*homogenization_RGC_Ngrains(3,i) &
+ 3_pInt*homogenization_RGC_Ngrains(1,i)*homogenization_RGC_Ngrains(2,i)*(homogenization_RGC_Ngrains(3,i)-1_pInt) &
+ 8_pInt ! (1) Average constitutive work, (2-4) Overall mismatch, (5) Average penalty energy,
! (6) Volume discrepancy, (7) Avg relaxation rate component, (8) Max relaxation rate component
enddo
#endif
end subroutine homogenization_RGC_init end subroutine homogenization_RGC_init
@ -374,14 +315,7 @@ end subroutine homogenization_RGC_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief partitions the deformation gradient onto the constituents !> @brief partitions the deformation gradient onto the constituents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
#ifdef NEWSTATE
subroutine homogenization_RGC_partitionDeformation(F,avgF,ip,el) subroutine homogenization_RGC_partitionDeformation(F,avgF,ip,el)
#else
subroutine homogenization_RGC_partitionDeformation(F,avgF,state,ip,el)
#endif
use prec, only: &
p_vec
use debug, only: & use debug, only: &
debug_level, & debug_level, &
debug_homogenization, & debug_homogenization, &
@ -391,10 +325,8 @@ subroutine homogenization_RGC_partitionDeformation(F,avgF,state,ip,el)
use material, only: & use material, only: &
homogenization_maxNgrains, & homogenization_maxNgrains, &
homogenization_Ngrains,& homogenization_Ngrains,&
#ifdef NEWSTATE
homogState, & homogState, &
mappingHomogenization, & mappingHomogenization, &
#endif
homogenization_typeInstance homogenization_typeInstance
use FEsolving, only: & use FEsolving, only: &
theInc,& theInc,&
@ -403,9 +335,6 @@ subroutine homogenization_RGC_partitionDeformation(F,avgF,state,ip,el)
implicit none implicit none
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned F per grain real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned F per grain
real(pReal), dimension (3,3), intent(in) :: avgF !< averaged F real(pReal), dimension (3,3), intent(in) :: avgF !< averaged F
#ifndef NEWSTATE
type(p_vec), intent(in) :: state
#endif
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
@ -437,11 +366,9 @@ subroutine homogenization_RGC_partitionDeformation(F,avgF,state,ip,el)
iGrain3 = homogenization_RGC_grain1to3(iGrain,homID) iGrain3 = homogenization_RGC_grain1to3(iGrain,homID)
do iFace = 1_pInt,nFace do iFace = 1_pInt,nFace
intFace = homogenization_RGC_getInterface(iFace,iGrain3) ! identifying 6 interfaces of each grain intFace = homogenization_RGC_getInterface(iFace,iGrain3) ! identifying 6 interfaces of each grain
#ifdef NEWSTATE
aVect = homogenization_RGC_relaxationVector(intFace,homID, ip, el) ! get the relaxation vectors for each interface from global relaxation vector array aVect = homogenization_RGC_relaxationVector(intFace,homID, ip, el) ! get the relaxation vectors for each interface from global relaxation vector array
#else
aVect = homogenization_RGC_relaxationVector(intFace,state,homID) ! get the relaxation vectors for each interface from global relaxation vector array
#endif
nVect = homogenization_RGC_interfaceNormal(intFace,ip,el) ! get the normal of each interface nVect = homogenization_RGC_interfaceNormal(intFace,ip,el) ! get the normal of each interface
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) & forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) &
F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! calculating deformation relaxations due to interface relaxation F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! calculating deformation relaxations due to interface relaxation
@ -470,13 +397,7 @@ end subroutine homogenization_RGC_partitionDeformation
!> @brief update the internal state of the homogenization scheme and tell whether "done" and !> @brief update the internal state of the homogenization scheme and tell whether "done" and
! "happy" with result ! "happy" with result
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
#ifdef NEWSTATE function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
function homogenization_RGC_updateState( P,F,F0,avgF,dt,dPdF,ip,el)
#else
function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el)
#endif
use prec, only: &
p_vec
use debug, only: & use debug, only: &
debug_level, & debug_level, &
debug_homogenization,& debug_homogenization,&
@ -490,10 +411,8 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el
use material, only: & use material, only: &
homogenization_maxNgrains, & homogenization_maxNgrains, &
homogenization_typeInstance, & homogenization_typeInstance, &
#ifdef NEWSTATE
homogState, & homogState, &
mappingHomogenization, & mappingHomogenization, &
#endif
homogenization_Ngrains homogenization_Ngrains
use numerics, only: & use numerics, only: &
absTol_RGC, & absTol_RGC, &
@ -507,10 +426,7 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el
refRelaxRate_RGC refRelaxRate_RGC
implicit none implicit none
#ifndef NEWSTATE
type(p_vec), intent(inout) :: state !< current state
type(p_vec), intent(in) :: state0 !< initial state
#endif
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: & real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: &
P,& !< array of P P,& !< array of P
F,& !< array of F F,& !< array of F
@ -554,28 +470,19 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el
! allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster ! allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster
allocate(resid(3_pInt*nIntFaceTot), source=0.0_pReal) allocate(resid(3_pInt*nIntFaceTot), source=0.0_pReal)
allocate(tract(nIntFaceTot,3), source=0.0_pReal) allocate(tract(nIntFaceTot,3), source=0.0_pReal)
#ifdef NEWSTATE
allocate(relax(3_pInt*nIntFaceTot)); relax= homogState(mappingHomogenization(2,ip,el))% & allocate(relax(3_pInt*nIntFaceTot)); relax= homogState(mappingHomogenization(2,ip,el))% &
state(1:3_pInt*nIntFaceTot,mappingHomogenization(1,ip,el)) state(1:3_pInt*nIntFaceTot,mappingHomogenization(1,ip,el))
allocate(drelax(3_pInt*nIntFaceTot)); drelax= homogState(mappingHomogenization(2,ip,el))% & allocate(drelax(3_pInt*nIntFaceTot)); drelax= homogState(mappingHomogenization(2,ip,el))% &
state(1:3_pInt*nIntFaceTot,mappingHomogenization(1,ip,el)) - & state(1:3_pInt*nIntFaceTot,mappingHomogenization(1,ip,el)) - &
homogState(mappingHomogenization(2,ip,el))% & homogState(mappingHomogenization(2,ip,el))% &
state0(1:3_pInt*nIntFaceTot,mappingHomogenization(1,ip,el)) state0(1:3_pInt*nIntFaceTot,mappingHomogenization(1,ip,el))
#else
allocate(relax(3_pInt*nIntFaceTot)); relax=state%p(1:3_pInt*nIntFaceTot)
allocate(drelax(3_pInt*nIntFaceTot)); drelax=state%p(1:3_pInt*nIntFaceTot)-state0%p(1:3_pInt*nIntFaceTot)
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! debugging the obtained state ! debugging the obtained state
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Obtained state: ' write(6,'(1x,a30)')'Obtained state: '
do i = 1_pInt,3_pInt*nIntFaceTot do i = 1_pInt,3_pInt*nIntFaceTot
#ifdef NEWSTATE
write(6,'(1x,2(e15.8,1x))')homogState(mappingHomogenization(2,ip,el))%state(i,mappingHomogenization(1,ip,el)) write(6,'(1x,2(e15.8,1x))')homogState(mappingHomogenization(2,ip,el))%state(i,mappingHomogenization(1,ip,el))
#else
write(6,'(1x,2(e15.8,1x))')state%p(i)
#endif
enddo enddo
write(6,*)' ' write(6,*)' '
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
@ -607,7 +514,7 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
!!!------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------
! computing the residual stress from the balance of traction at all (interior) interfaces ! computing the residual stress from the balance of traction at all (interior) interfaces
do iNum = 1_pInt,nIntFaceTot do iNum = 1_pInt,nIntFaceTot
faceID = homogenization_RGC_interface1to4(iNum,homID) ! identifying the interface ID in local coordinate system (4-dimensional index) faceID = homogenization_RGC_interface1to4(iNum,homID) ! identifying the interface ID in local coordinate system (4-dimensional index)
@ -689,13 +596,8 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! compute/update the state for postResult, i.e., all energy densities computed by time-integration ! compute/update the state for postResult, i.e., all energy densities computed by time-integration
#ifdef NEWSTATE
constitutiveWork = homogState(mappingHomogenization(2,ip,el))%state(3*nIntFaceTot+1,mappingHomogenization(1,ip,el)) constitutiveWork = homogState(mappingHomogenization(2,ip,el))%state(3*nIntFaceTot+1,mappingHomogenization(1,ip,el))
penaltyEnergy = homogState(mappingHomogenization(2,ip,el))%state(3*nIntFaceTot+5,mappingHomogenization(1,ip,el)) penaltyEnergy = homogState(mappingHomogenization(2,ip,el))%state(3*nIntFaceTot+5,mappingHomogenization(1,ip,el))
#else
constitutiveWork = state%p(3*nIntFaceTot+1)
penaltyEnergy = state%p(3*nIntFaceTot+5)
#endif
do iGrain = 1_pInt,homogenization_Ngrains(mesh_element(3,el)) ! time-integration loop for the calculating the work and energy do iGrain = 1_pInt,homogenization_Ngrains(mesh_element(3,el)) ! time-integration loop for the calculating the work and energy
do i = 1_pInt,3_pInt do i = 1_pInt,3_pInt
do j = 1_pInt,3_pInt do j = 1_pInt,3_pInt
@ -704,7 +606,6 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el
enddo enddo
enddo enddo
enddo enddo
#ifdef NEWSTATE
homogState(mappingHomogenization(2,ip,el))% & homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+1,mappingHomogenization(1,ip,el)) = constitutiveWork ! the bulk mechanical/constitutive work state(3*nIntFaceTot+1,mappingHomogenization(1,ip,el)) = constitutiveWork ! the bulk mechanical/constitutive work
homogState(mappingHomogenization(2,ip,el))% & homogState(mappingHomogenization(2,ip,el))% &
@ -723,17 +624,6 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el
homogState(mappingHomogenization(2,ip,el))% & homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+8,mappingHomogenization(1,ip,el)) = maxval(abs(drelax))/dt ! the maximum rate of relaxation vectors state(3*nIntFaceTot+8,mappingHomogenization(1,ip,el)) = maxval(abs(drelax))/dt ! the maximum rate of relaxation vectors
#else
state%p(3*nIntFaceTot+1) = constitutiveWork ! the bulk mechanical/constitutive work
state%p(3*nIntFaceTot+2) = sum(NN(1,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e1-direction
state%p(3*nIntFaceTot+3) = sum(NN(2,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e2-direction
state%p(3*nIntFaceTot+4) = sum(NN(3,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e3-direction
state%p(3*nIntFaceTot+5) = penaltyEnergy ! the overall penalty energy
state%p(3*nIntFaceTot+6) = volDiscrep ! the overall volume discrepancy
state%p(3*nIntFaceTot+7) = sum(abs(drelax))/dt/real(3_pInt*nIntFaceTot,pReal) ! the average rate of relaxation vectors
state%p(3*nIntFaceTot+8) = maxval(abs(drelax))/dt ! the maximum rate of relaxation vectors
#endif
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then .and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
@ -778,7 +668,7 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el
endif endif
!!!------------------------------------------------------------------------------------------------ !---------------------------------------------------------------------------------------------------
! construct the global Jacobian matrix for updating the global relaxation vector array when ! construct the global Jacobian matrix for updating the global relaxation vector array when
! convergence is not yet reached ... ! convergence is not yet reached ...
@ -848,13 +738,8 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el
do ipert = 1_pInt,3_pInt*nIntFaceTot do ipert = 1_pInt,3_pInt*nIntFaceTot
p_relax = relax p_relax = relax
p_relax(ipert) = relax(ipert) + pPert_RGC ! perturb the relaxation vector p_relax(ipert) = relax(ipert) + pPert_RGC ! perturb the relaxation vector
#ifdef NEWSTATE
homogState(mappingHomogenization(2,ip,el))%state(1:3*nIntFaceTot,mappingHomogenization(1,ip,el)) = p_relax homogState(mappingHomogenization(2,ip,el))%state(1:3*nIntFaceTot,mappingHomogenization(1,ip,el)) = p_relax
call homogenization_RGC_grainDeformation(pF,avgF,ip,el) ! compute the grains deformation from perturbed state call homogenization_RGC_grainDeformation(pF,avgF,ip,el) ! compute the grains deformation from perturbed state
#else
state%p(1:3*nIntFaceTot) = p_relax
call homogenization_RGC_grainDeformation(pF,avgF,state,ip,el) ! compute the grains deformation from perturbed state
#endif
call homogenization_RGC_stressPenalty(pR,pNN,avgF,pF,ip,el,homID) ! compute stress penalty due to interface mismatch from perturbed state call homogenization_RGC_stressPenalty(pR,pNN,avgF,pF,ip,el,homID) ! compute stress penalty due to interface mismatch from perturbed state
call homogenization_RGC_volumePenalty(pD,volDiscrep,pF,avgF,ip,el) ! compute stress penalty due to volume discrepancy from perturbed state call homogenization_RGC_volumePenalty(pD,volDiscrep,pF,avgF,ip,el) ! compute stress penalty due to volume discrepancy from perturbed state
@ -969,11 +854,7 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el
enddo enddo
enddo enddo
relax = relax + drelax ! Updateing the state variable for the next iteration relax = relax + drelax ! Updateing the state variable for the next iteration
#ifdef NEWSTATE
homogState(mappingHomogenization(2,ip,el))%state(1:3*nIntFaceTot,mappingHomogenization(1,ip,el)) = relax homogState(mappingHomogenization(2,ip,el))%state(1:3*nIntFaceTot,mappingHomogenization(1,ip,el)) = relax
#else
state%p(1:3*nIntFaceTot) = relax
#endif
if (any(abs(drelax) > maxdRelax_RGC)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large if (any(abs(drelax) > maxdRelax_RGC)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large
homogenization_RGC_updateState = [.true.,.false.] homogenization_RGC_updateState = [.true.,.false.]
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
@ -989,11 +870,7 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Returned state: ' write(6,'(1x,a30)')'Returned state: '
do i = 1_pInt,3_pInt*nIntFaceTot do i = 1_pInt,3_pInt*nIntFaceTot
#ifdef NEWSTATE
write(6,'(1x,2(e15.8,1x))')homogState(mappingHomogenization(2,ip,el))%state(i,mappingHomogenization(1,ip,el)) write(6,'(1x,2(e15.8,1x))')homogState(mappingHomogenization(2,ip,el))%state(i,mappingHomogenization(1,ip,el))
#else
write(6,'(1x,2(e15.8,1x))')state%p(i)
#endif
enddo enddo
write(6,*)' ' write(6,*)' '
flush(6) flush(6)
@ -1016,10 +893,8 @@ subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,
use mesh, only: mesh_element use mesh, only: mesh_element
use material, only: & use material, only: &
homogenization_maxNgrains, & homogenization_maxNgrains, &
#ifdef NEWSTATE
homogState, & homogState, &
mappingHomogenization, & mappingHomogenization, &
#endif
homogenization_Ngrains, & homogenization_Ngrains, &
homogenization_typeInstance homogenization_typeInstance
use math, only: math_Plain3333to99 use math, only: math_Plain3333to99
@ -1064,22 +939,14 @@ end subroutine homogenization_RGC_averageStressAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of homogenization results for post file inclusion !> @brief return array of homogenization results for post file inclusion
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
#ifdef NEWSTATE
pure function homogenization_RGC_postResults(ip,el,avgP,avgF) pure function homogenization_RGC_postResults(ip,el,avgP,avgF)
#else
pure function homogenization_RGC_postResults(state,ip,el,avgP,avgF)
#endif
use prec, only: &
p_vec
use mesh, only: & use mesh, only: &
mesh_element, & mesh_element, &
mesh_ipCoordinates mesh_ipCoordinates
use material, only: & use material, only: &
homogenization_typeInstance,& homogenization_typeInstance,&
#ifdef NEWSTATE
homogState, & homogState, &
mappingHomogenization, & mappingHomogenization, &
#endif
homogenization_Noutput homogenization_Noutput
use crystallite, only: & use crystallite, only: &
crystallite_temperature crystallite_temperature
@ -1091,10 +958,7 @@ pure function homogenization_RGC_postResults(state,ip,el,avgP,avgF)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
avgP, & !< average stress at material point avgP, & !< average stress at material point
avgF !< average deformation gradient at material point avgF !< average deformation gradient at material point
#ifndef NEWSTATE
type(p_vec), intent(in) :: &
state ! my State
#endif
integer(pInt) homID,o,c,nIntFaceTot integer(pInt) homID,o,c,nIntFaceTot
real(pReal), dimension(homogenization_RGC_sizePostResults(homogenization_typeInstance(mesh_element(3,el)))) :: & real(pReal), dimension(homogenization_RGC_sizePostResults(homogenization_typeInstance(mesh_element(3,el)))) :: &
homogenization_RGC_postResults homogenization_RGC_postResults
@ -1121,63 +985,33 @@ pure function homogenization_RGC_postResults(state,ip,el,avgP,avgF)
homogenization_RGC_postResults(c+1_pInt:c+3_pInt) = mesh_ipCoordinates(1:3,ip,el) ! current ip coordinates homogenization_RGC_postResults(c+1_pInt:c+3_pInt) = mesh_ipCoordinates(1:3,ip,el) ! current ip coordinates
c = c + 3_pInt c = c + 3_pInt
case (constitutivework_ID) case (constitutivework_ID)
#ifdef NEWSTATE
homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% & homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+1,mappingHomogenization(1,ip,el)) state(3*nIntFaceTot+1,mappingHomogenization(1,ip,el))
#else
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+1)
#endif
c = c + 1_pInt c = c + 1_pInt
case (magnitudemismatch_ID) case (magnitudemismatch_ID)
#ifdef NEWSTATE
homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% & homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+2,mappingHomogenization(1,ip,el)) state(3*nIntFaceTot+2,mappingHomogenization(1,ip,el))
homogenization_RGC_postResults(c+2) = homogState(mappingHomogenization(2,ip,el))% & homogenization_RGC_postResults(c+2) = homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+3,mappingHomogenization(1,ip,el)) state(3*nIntFaceTot+3,mappingHomogenization(1,ip,el))
homogenization_RGC_postResults(c+3) = homogState(mappingHomogenization(2,ip,el))% & homogenization_RGC_postResults(c+3) = homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+4,mappingHomogenization(1,ip,el)) state(3*nIntFaceTot+4,mappingHomogenization(1,ip,el))
#else
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+2)
homogenization_RGC_postResults(c+2) = state%p(3*nIntFaceTot+3)
homogenization_RGC_postResults(c+3) = state%p(3*nIntFaceTot+4)
#endif
c = c + 3_pInt c = c + 3_pInt
case (penaltyenergy_ID) case (penaltyenergy_ID)
#ifdef NEWSTATE
homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% & homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+5,mappingHomogenization(1,ip,el)) state(3*nIntFaceTot+5,mappingHomogenization(1,ip,el))
c = c + 1_pInt c = c + 1_pInt
#else
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+5)
c = c + 1_pInt
#endif
case (volumediscrepancy_ID) case (volumediscrepancy_ID)
#ifdef NEWSTATE
homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% & homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+6,mappingHomogenization(1,ip,el)) state(3*nIntFaceTot+6,mappingHomogenization(1,ip,el))
c = c + 1_pInt c = c + 1_pInt
#else
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+6)
c = c + 1_pInt
#endif
case (averagerelaxrate_ID) case (averagerelaxrate_ID)
#ifdef NEWSTATE
homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% & homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+7,mappingHomogenization(1,ip,el)) state(3*nIntFaceTot+7,mappingHomogenization(1,ip,el))
c = c + 1_pInt c = c + 1_pInt
#else
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+7)
c = c + 1_pInt
#endif
case (maximumrelaxrate_ID) case (maximumrelaxrate_ID)
#ifdef NEWSTATE
homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% & homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+8,mappingHomogenization(1,ip,el)) state(3*nIntFaceTot+8,mappingHomogenization(1,ip,el))
c = c + 1_pInt c = c + 1_pInt
#else
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+8)
c = c + 1_pInt
#endif
end select end select
enddo enddo
@ -1203,10 +1037,8 @@ subroutine homogenization_RGC_stressPenalty(rPen,nMis,avgF,fDef,ip,el,homID)
math_invert33 math_invert33
use material, only: & use material, only: &
homogenization_maxNgrains,& homogenization_maxNgrains,&
#ifdef NEWSTATE
homogState, & homogState, &
mappingHomogenization, & mappingHomogenization, &
#endif
homogenization_Ngrains homogenization_Ngrains
use numerics, only: & use numerics, only: &
xSmoo_RGC xSmoo_RGC
@ -1346,10 +1178,6 @@ subroutine homogenization_RGC_volumePenalty(vPen,vDiscrep,fDef,fAvg,ip,el)
math_inv33 math_inv33
use material, only: & use material, only: &
homogenization_maxNgrains,& homogenization_maxNgrains,&
#ifdef NEWSTATE
homogState, &
mappingHomogenization, &
#endif
homogenization_Ngrains homogenization_Ngrains
use numerics, only: & use numerics, only: &
maxVolDiscr_RGC,& maxVolDiscr_RGC,&
@ -1409,11 +1237,6 @@ function homogenization_RGC_surfaceCorrection(avgF,ip,el)
use math, only: & use math, only: &
math_invert33, & math_invert33, &
math_mul33x33 math_mul33x33
#ifdef NEWSTATE
use material, only: &
homogState, &
mappingHomogenization
#endif
implicit none implicit none
real(pReal), dimension(3) :: homogenization_RGC_surfaceCorrection real(pReal), dimension(3) :: homogenization_RGC_surfaceCorrection
@ -1483,28 +1306,15 @@ end function homogenization_RGC_equivalentModuli
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief collect relaxation vectors of an interface !> @brief collect relaxation vectors of an interface
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
#ifdef NEWSTATE
function homogenization_RGC_relaxationVector(intFace,homID, ip, el) function homogenization_RGC_relaxationVector(intFace,homID, ip, el)
#else
function homogenization_RGC_relaxationVector(intFace,state,homID)
#endif
use prec, only: &
p_vec
#ifdef NEWSTATE
use material, only: & use material, only: &
homogState, & homogState, &
mappingHomogenization mappingHomogenization
#endif
implicit none implicit none
#ifdef NEWSTATE
integer(pInt), intent(in) :: ip, el integer(pInt), intent(in) :: ip, el
#endif
real(pReal), dimension (3) :: homogenization_RGC_relaxationVector real(pReal), dimension (3) :: homogenization_RGC_relaxationVector
integer(pInt), dimension (4), intent(in) :: intFace !< set of interface ID in 4D array (normal and position) integer(pInt), dimension (4), intent(in) :: intFace !< set of interface ID in 4D array (normal and position)
#ifndef NEWSTATE
type(p_vec), intent(in) :: state !< set of global relaxation vectors
#endif
integer(pInt), dimension (3) :: nGDim integer(pInt), dimension (3) :: nGDim
integer(pInt) :: & integer(pInt) :: &
iNum, & iNum, &
@ -1515,12 +1325,8 @@ function homogenization_RGC_relaxationVector(intFace,state,homID)
homogenization_RGC_relaxationVector = 0.0_pReal homogenization_RGC_relaxationVector = 0.0_pReal
nGDim = homogenization_RGC_Ngrains(1:3,homID) nGDim = homogenization_RGC_Ngrains(1:3,homID)
iNum = homogenization_RGC_interface4to1(intFace,homID) ! identify the position of the interface in global state array iNum = homogenization_RGC_interface4to1(intFace,homID) ! identify the position of the interface in global state array
#ifdef NEWSTATE
if (iNum > 0_pInt) homogenization_RGC_relaxationVector = homogState(mappingHomogenization(2,ip,el))% & if (iNum > 0_pInt) homogenization_RGC_relaxationVector = homogState(mappingHomogenization(2,ip,el))% &
state((3*iNum-2):(3*iNum),mappingHomogenization(1,ip,el)) ! get the corresponding entries state((3*iNum-2):(3*iNum),mappingHomogenization(1,ip,el)) ! get the corresponding entries
#else
if (iNum > 0_pInt) homogenization_RGC_relaxationVector = state%p((3*iNum-2):(3*iNum)) ! get the corresponding entries
#endif
end function homogenization_RGC_relaxationVector end function homogenization_RGC_relaxationVector
@ -1734,30 +1540,17 @@ end function homogenization_RGC_interface1to4
!> @brief calculating the grain deformation gradient (the same with !> @brief calculating the grain deformation gradient (the same with
! homogenization_RGC_partionDeformation, but used only for perturbation scheme) ! homogenization_RGC_partionDeformation, but used only for perturbation scheme)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
#ifdef NEWSTATE
subroutine homogenization_RGC_grainDeformation(F, avgF, ip, el) subroutine homogenization_RGC_grainDeformation(F, avgF, ip, el)
#else
subroutine homogenization_RGC_grainDeformation(F, avgF, state, ip, el)
#endif
use prec, only: &
p_vec
use mesh, only: & use mesh, only: &
mesh_element mesh_element
use material, only: & use material, only: &
homogenization_maxNgrains,& homogenization_maxNgrains,&
homogenization_Ngrains, & homogenization_Ngrains, &
#ifdef NEWSTATE
homogState, &
mappingHomogenization, &
#endif
homogenization_typeInstance homogenization_typeInstance
implicit none implicit none
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned F per grain real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned F per grain
real(pReal), dimension (3,3), intent(in) :: avgF !< real(pReal), dimension (3,3), intent(in) :: avgF !<
#ifndef NEWSTATE
type(p_vec), intent(in) :: state
#endif
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
el, & !< element number el, & !< element number
ip !< integration point number ip !< integration point number
@ -1775,11 +1568,7 @@ subroutine homogenization_RGC_grainDeformation(F, avgF, state, ip, el)
iGrain3 = homogenization_RGC_grain1to3(iGrain,homID) iGrain3 = homogenization_RGC_grain1to3(iGrain,homID)
do iFace = 1_pInt,nFace do iFace = 1_pInt,nFace
intFace = homogenization_RGC_getInterface(iFace,iGrain3) intFace = homogenization_RGC_getInterface(iFace,iGrain3)
#ifdef NEWSTATE
aVect = homogenization_RGC_relaxationVector(intFace,homID, ip, el) aVect = homogenization_RGC_relaxationVector(intFace,homID, ip, el)
#else
aVect = homogenization_RGC_relaxationVector(intFace,state,homID)
#endif
nVect = homogenization_RGC_interfaceNormal(intFace,ip,el) nVect = homogenization_RGC_interfaceNormal(intFace,ip,el)
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) & forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) &
F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations

77
code/homogenization_isostrain.f90
View File

@ -18,10 +18,8 @@ module homogenization_isostrain
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
homogenization_isostrain_output !< name of each post result output homogenization_isostrain_output !< name of each post result output
#ifdef NEWSTATE integer(pInt), dimension(:), allocatable, private :: &
integer(pInt), dimension(:), allocatable, private :: &
homogenization_isostrain_Noutput !< number of outputs per homog instance homogenization_isostrain_Noutput !< number of outputs per homog instance
#endif
integer(pInt), dimension(:), allocatable, private :: & integer(pInt), dimension(:), allocatable, private :: &
homogenization_isostrain_Ngrains homogenization_isostrain_Ngrains
enum, bind(c) enum, bind(c)
@ -56,14 +54,11 @@ contains
subroutine homogenization_isostrain_init(fileUnit) subroutine homogenization_isostrain_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: & use prec, only: &
pReal, & pReal
pInt
#ifdef NEWSTATE
use debug, only: & use debug, only: &
debug_HOMOGENIZATION, & debug_HOMOGENIZATION, &
debug_level, & debug_level, &
debug_levelBasic debug_levelBasic
#endif
use IO use IO
use material use material
@ -75,14 +70,10 @@ subroutine homogenization_isostrain_init(fileUnit)
section = 0_pInt, i, j, output, mySize, o section = 0_pInt, i, j, output, mySize, o
integer :: & integer :: &
maxNinstance, & maxNinstance, &
#ifdef NEWSTATE
homog, & homog, &
NofMyHomog, & instance
instance, & integer :: &
sizeHState, & NofMyHomog ! no pInt (stores a system dependen value from 'count'
#endif
k
! no pInt (stores a system dependen value from 'count'
character(len=65536) :: & character(len=65536) :: &
tag = '', & tag = '', &
line = '' line = ''
@ -94,16 +85,13 @@ subroutine homogenization_isostrain_init(fileUnit)
maxNinstance = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID) maxNinstance = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)
if (maxNinstance == 0) return if (maxNinstance == 0) return
#ifdef NEWSTATE
if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
#endif
allocate(homogenization_isostrain_sizePostResults(maxNinstance), source=0_pInt) allocate(homogenization_isostrain_sizePostResults(maxNinstance), source=0_pInt)
allocate(homogenization_isostrain_sizePostResult(maxval(homogenization_Noutput),maxNinstance), & allocate(homogenization_isostrain_sizePostResult(maxval(homogenization_Noutput),maxNinstance), &
source=0_pInt) source=0_pInt)
#ifdef NEWSTATE allocate(homogenization_isostrain_Noutput(maxNinstance), source=0_pInt)
allocate(homogenization_isostrain_Noutput(maxNinstance), source=0_pInt)
#endif
allocate(homogenization_isostrain_Ngrains(maxNinstance), source=0_pInt) allocate(homogenization_isostrain_Ngrains(maxNinstance), source=0_pInt)
allocate(homogenization_isostrain_mapping(maxNinstance), source=average_ID) allocate(homogenization_isostrain_mapping(maxNinstance), source=average_ID)
allocate(homogenization_isostrain_output(maxval(homogenization_Noutput),maxNinstance)) allocate(homogenization_isostrain_output(maxval(homogenization_Noutput),maxNinstance))
@ -116,7 +104,7 @@ subroutine homogenization_isostrain_init(fileUnit)
line = IO_read(fileUnit) line = IO_read(fileUnit)
enddo enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of homogenization part parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of homogenization part
line = IO_read(fileUnit) line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part if (IO_getTag(line,'<','>') /= '') then ! stop at next part
@ -143,29 +131,19 @@ subroutine homogenization_isostrain_init(fileUnit)
homogenization_isostrain_output(output,i) = IO_lc(IO_stringValue(line,positions,2_pInt)) homogenization_isostrain_output(output,i) = IO_lc(IO_stringValue(line,positions,2_pInt))
select case(homogenization_isostrain_output(output,i)) select case(homogenization_isostrain_output(output,i))
case('nconstituents','ngrains') case('nconstituents','ngrains')
#ifdef NEWSTATE
homogenization_isostrain_Noutput(i) = homogenization_isostrain_Noutput(i) + 1_pInt homogenization_isostrain_Noutput(i) = homogenization_isostrain_Noutput(i) + 1_pInt
#endif
homogenization_isostrain_outputID(output,i) = nconstituents_ID homogenization_isostrain_outputID(output,i) = nconstituents_ID
case('temperature') case('temperature')
#ifdef NEWSTATE
homogenization_isostrain_Noutput(i) = homogenization_isostrain_Noutput(i) + 1_pInt homogenization_isostrain_Noutput(i) = homogenization_isostrain_Noutput(i) + 1_pInt
#endif
homogenization_isostrain_outputID(output,i) = temperature_ID homogenization_isostrain_outputID(output,i) = temperature_ID
case('ipcoords') case('ipcoords')
#ifdef NEWSTATE
homogenization_isostrain_Noutput(i) = homogenization_isostrain_Noutput(i) + 1_pInt homogenization_isostrain_Noutput(i) = homogenization_isostrain_Noutput(i) + 1_pInt
#endif
homogenization_isostrain_outputID(output,i) = ipcoords_ID homogenization_isostrain_outputID(output,i) = ipcoords_ID
case('avgdefgrad','avgf') case('avgdefgrad','avgf')
#ifdef NEWSTATE
homogenization_isostrain_Noutput(i) = homogenization_isostrain_Noutput(i) + 1_pInt homogenization_isostrain_Noutput(i) = homogenization_isostrain_Noutput(i) + 1_pInt
#endif
homogenization_isostrain_outputID(output,i) = avgdefgrad_ID homogenization_isostrain_outputID(output,i) = avgdefgrad_ID
case('avgp','avgfirstpiola','avg1stpiola') case('avgp','avgfirstpiola','avg1stpiola')
#ifdef NEWSTATE
homogenization_isostrain_Noutput(i) = homogenization_isostrain_Noutput(i) + 1_pInt homogenization_isostrain_Noutput(i) = homogenization_isostrain_Noutput(i) + 1_pInt
#endif
homogenization_isostrain_outputID(output,i) = avgfirstpiola_ID homogenization_isostrain_outputID(output,i) = avgfirstpiola_ID
case default case default
call IO_error(105_pInt,ext_msg=IO_stringValue(line,positions,2_pInt)//& call IO_error(105_pInt,ext_msg=IO_stringValue(line,positions,2_pInt)//&
@ -189,7 +167,6 @@ subroutine homogenization_isostrain_init(fileUnit)
endif endif
enddo parsingFile enddo parsingFile
#ifdef NEWSTATE
initializeInstances: do homog = 1_pInt, material_Nhomogenization initializeInstances: do homog = 1_pInt, material_Nhomogenization
myHomog: if (homogenization_type(homog) == HOMOGENIZATION_ISOSTRAIN_ID) then myHomog: if (homogenization_type(homog) == HOMOGENIZATION_ISOSTRAIN_ID) then
NofMyHomog = count(material_homog == homog) NofMyHomog = count(material_homog == homog)
@ -216,43 +193,15 @@ subroutine homogenization_isostrain_init(fileUnit)
enddo outputsLoop enddo outputsLoop
! allocate state arrays ! allocate state arrays
sizeHState = 0_pInt homogState(homog)%sizeState = 0_pInt
homogState(homog)%sizeState = sizeHState
homogState(homog)%sizePostResults = homogenization_isostrain_sizePostResults(instance) homogState(homog)%sizePostResults = homogenization_isostrain_sizePostResults(instance)
allocate(homogState(homog)%state0 ( sizeHState,NofMyHomog), source=0.0_pReal) allocate(homogState(homog)%state0 (0_pInt,NofMyHomog), source=0.0_pReal)
allocate(homogState(homog)%subState0 ( sizeHState,NofMyHomog), source=0.0_pReal) allocate(homogState(homog)%subState0(0_pInt,NofMyHomog), source=0.0_pReal)
allocate(homogState(homog)%state ( sizeHState,NofMyHomog), source=0.0_pReal) allocate(homogState(homog)%state (0_pInt,NofMyHomog), source=0.0_pReal)
endif myHomog endif myHomog
enddo initializeInstances enddo initializeInstances
#else
do i = 1,maxNinstance
do j = 1_pInt,maxval(homogenization_Noutput)
select case(homogenization_isostrain_outputID(j,i))
case(nconstituents_ID, temperature_ID)
mySize = 1_pInt
case(ipcoords_ID)
mySize = 3_pInt
case(avgdefgrad_ID, avgfirstpiola_ID)
mySize = 9_pInt
case default
mySize = 0_pInt
end select
outputFound: if (mySize > 0_pInt) then
homogenization_isostrain_sizePostResult(j,i) = mySize
homogenization_isostrain_sizePostResults(i) = &
homogenization_isostrain_sizePostResults(i) + mySize
endif outputFound
enddo
enddo
#endif
end subroutine homogenization_isostrain_init end subroutine homogenization_isostrain_init

26
code/homogenization_none.f90
View File

@ -29,36 +29,28 @@ subroutine homogenization_none_init()
use material use material
implicit none implicit none
#ifdef NEWSTATE integer(pInt) :: &
integer :: &
homog, & homog, &
NofMyHomog, & NofMyHomog
instance, &
sizeHState
#endif
write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_NONE_label//' init -+>>>' write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_NONE_label//' init -+>>>'
write(6,'(a)') ' $Id$' write(6,'(a)') ' $Id$'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
#ifdef NEWSTATE
initializeInstances: do homog = 1_pInt, material_Nhomogenization initializeInstances: do homog = 1_pInt, material_Nhomogenization
myhomog: if (homogenization_type(homog) == HOMOGENIZATION_none_ID) then myhomog: if (homogenization_type(homog) == HOMOGENIZATION_none_ID) then
NofMyHomog = count(material_homog == homog) NofMyHomog = count(material_homog == homog)
! instance = phase_plasticityInstance(phase) homogState(homog)%sizeState = 0_pInt
! allocate homogenization state arrays
sizeHState = 0_pInt
homogState(homog)%sizeState = sizeHState
homogState(homog)%sizePostResults = 0_pInt homogState(homog)%sizePostResults = 0_pInt
allocate(homogState(homog)%state0 ( sizeHState,NofMyHomog), source=0.0_pReal) allocate(homogState(homog)%state0 (0_pInt,NofMyHomog), source=0.0_pReal)
allocate(homogState(homog)%subState0 ( sizeHState,NofMyHomog), source=0.0_pReal) allocate(homogState(homog)%subState0(0_pInt,NofMyHomog), source=0.0_pReal)
allocate(homogState(homog)%state ( sizeHState,NofMyHomog), source=0.0_pReal) allocate(homogState(homog)%state (0_pInt,NofMyHomog), source=0.0_pReal)
endif myhomog endif myhomog
enddo initializeInstances enddo initializeInstances
#endif
end subroutine homogenization_none_init end subroutine homogenization_none_init

14
code/prec.f90
View File

@ -77,18 +77,12 @@ module prec
RK4dotState RK4dotState
real(pReal), allocatable, dimension(:,:,:) :: RKCK45dotState real(pReal), allocatable, dimension(:,:,:) :: RKCK45dotState
end type end type
#ifdef NEWSTATE #ifdef NEWSTATE
type, public :: hState type, public :: tFieldData
integer(pInt) :: sizeState = 0_pInt , & integer(pInt) :: sizeField = 0_pInt , &
sizePostResults = 0_pInt sizePostResults = 0_pInt
real(pReal), allocatable, dimension(:,:) :: state, & ! material points, state size real(pReal), allocatable, dimension(:,:) :: field ! material points, state size
state0, &
subState0
end type
type, public :: fState
integer(pInt) :: sizeState = 0_pInt , &
sizePostResults = 0_pInt
real(pReal), allocatable, dimension(:,:) :: state ! material points, state size
end type end type
#endif #endif