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

Merge branch 'clean-constitutive' into 'development' 

Clean constitutive

See merge request damask/DAMASK!309
This commit is contained in:
Sharan Roongta 2021-01-06 14:22:56 +01:00
parent 2791f3d192 5f5d9ed908
commit f8dd5df0cc
20 changed files with 1260 additions and 1528 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
src/CPFEM.f90
View File

@ -260,7 +260,7 @@ end subroutine CPFEM_general
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine CPFEM_forward subroutine CPFEM_forward
call crystallite_forward call homogenization_forward
call constitutive_forward call constitutive_forward
end subroutine CPFEM_forward end subroutine CPFEM_forward
@ -277,7 +277,6 @@ subroutine CPFEM_results(inc,time)
call results_openJobFile call results_openJobFile
call results_addIncrement(inc,time) call results_addIncrement(inc,time)
call constitutive_results call constitutive_results
call crystallite_results
call homogenization_results call homogenization_results
call discretization_results call discretization_results
call results_finalizeIncrement call results_finalizeIncrement

3
src/CPFEM2.f90
View File

@ -97,7 +97,7 @@ end subroutine CPFEM_restartWrite
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine CPFEM_forward subroutine CPFEM_forward
call crystallite_forward call homogenization_forward
call constitutive_forward call constitutive_forward
end subroutine CPFEM_forward end subroutine CPFEM_forward
@ -114,7 +114,6 @@ subroutine CPFEM_results(inc,time)
call results_openJobFile call results_openJobFile
call results_addIncrement(inc,time) call results_addIncrement(inc,time)
call constitutive_results call constitutive_results
call crystallite_results
call homogenization_results call homogenization_results
call discretization_results call discretization_results
call results_finalizeIncrement call results_finalizeIncrement

1244
src/constitutive.f90
View File

File diff suppressed because it is too large Load Diff

69
src/constitutive_damage.f90
View File

@ -1,5 +1,5 @@
!---------------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------------
!> @brief internal microstructure state for all damage sources and kinematics constitutive models !> @brief internal microstructure state for all damage sources and kinematics constitutive models
!---------------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------------
submodule(constitutive) constitutive_damage submodule(constitutive) constitutive_damage
@ -8,7 +8,7 @@ submodule(constitutive) constitutive_damage
module function source_damage_anisoBrittle_init(source_length) result(mySources) module function source_damage_anisoBrittle_init(source_length) result(mySources)
integer, intent(in) :: source_length integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources logical, dimension(:,:), allocatable :: mySources
end function source_damage_anisoBrittle_init end function source_damage_anisoBrittle_init
module function source_damage_anisoDuctile_init(source_length) result(mySources) module function source_damage_anisoDuctile_init(source_length) result(mySources)
integer, intent(in) :: source_length integer, intent(in) :: source_length
@ -23,7 +23,7 @@ submodule(constitutive) constitutive_damage
module function source_damage_isoDuctile_init(source_length) result(mySources) module function source_damage_isoDuctile_init(source_length) result(mySources)
integer, intent(in) :: source_length integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources logical, dimension(:,:), allocatable :: mySources
end function source_damage_isoDuctile_init end function source_damage_isoDuctile_init
module function kinematics_cleavage_opening_init(kinematics_length) result(myKinematics) module function kinematics_cleavage_opening_init(kinematics_length) result(myKinematics)
integer, intent(in) :: kinematics_length integer, intent(in) :: kinematics_length
@ -39,14 +39,14 @@ submodule(constitutive) constitutive_damage
module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
integer, intent(in) :: & integer, intent(in) :: &
phase, & !< phase ID of element phase, & !< phase ID of element
constituent !< position of element within its phase instance constituent !< position of element within its phase instance
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi !< damage parameter phi !< damage parameter
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
localphiDot, & localphiDot, &
dLocalphiDot_dPhi dLocalphiDot_dPhi
end subroutine source_damage_anisoBrittle_getRateAndItsTangent end subroutine source_damage_anisoBrittle_getRateAndItsTangent
module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
integer, intent(in) :: & integer, intent(in) :: &
phase, & !< phase ID of element phase, & !< phase ID of element
@ -129,7 +129,7 @@ module subroutine damage_init
allocate(sourceState(ph)%p(phase_Nsources(ph))) allocate(sourceState(ph)%p(phase_Nsources(ph)))
enddo enddo
allocate(phase_source(maxval(phase_Nsources),phases%length), source = SOURCE_undefined_ID) allocate(phase_source(maxval(phase_Nsources),phases%length), source = SOURCE_undefined_ID)
! initialize source mechanisms ! initialize source mechanisms
if(maxval(phase_Nsources) /= 0) then if(maxval(phase_Nsources) /= 0) then
@ -141,19 +141,19 @@ module subroutine damage_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize kinematic mechanisms ! initialize kinematic mechanisms
allocate(phase_Nkinematics(phases%length),source = 0) allocate(phase_Nkinematics(phases%length),source = 0)
do ph = 1,phases%length do ph = 1,phases%length
phase => phases%get(ph) phase => phases%get(ph)
kinematics => phase%get('kinematics',defaultVal=emptyList) kinematics => phase%get('kinematics',defaultVal=emptyList)
phase_Nkinematics(ph) = kinematics%length phase_Nkinematics(ph) = kinematics%length
enddo enddo
allocate(phase_kinematics(maxval(phase_Nkinematics),phases%length), source = KINEMATICS_undefined_ID) allocate(phase_kinematics(maxval(phase_Nkinematics),phases%length), source = KINEMATICS_undefined_ID)
if(maxval(phase_Nkinematics) /= 0) then if(maxval(phase_Nkinematics) /= 0) then
where(kinematics_cleavage_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_cleavage_opening_ID where(kinematics_cleavage_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_cleavage_opening_ID
where(kinematics_slipplane_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_slipplane_opening_ID where(kinematics_slipplane_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_slipplane_opening_ID
endif endif
end subroutine damage_init end subroutine damage_init
@ -167,7 +167,7 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi !< damage parameter phi !< damage parameter
real(pReal), intent(inout) :: & real(pReal), intent(inout) :: &
phiDot, & phiDot, &
dPhiDot_dPhi dPhiDot_dPhi
@ -183,7 +183,7 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi
phiDot = 0.0_pReal phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal dPhiDot_dPhi = 0.0_pReal
do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el)) do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
phase = material_phaseAt(grain,el) phase = material_phaseAt(grain,el)
constituent = material_phasememberAt(grain,ip,el) constituent = material_phasememberAt(grain,ip,el)
@ -217,32 +217,35 @@ end subroutine constitutive_damage_getRateAndItsTangents
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief writes damage sources results to HDF5 output file !< @brief writes damage sources results to HDF5 output file
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module subroutine damage_results module subroutine damage_results(group,ph)
integer :: p,i character(len=*), intent(in) :: group
character(len=pStringLen) :: group integer, intent(in) :: ph
do p = 1, size(material_name_phase) integer :: so
sourceLoop: do i = 1, phase_Nsources(p) sourceLoop: do so = 1, phase_Nsources(ph)
group = trim('current/phase')//'/'//trim(material_name_phase(p))
group = trim(group)//'/sources'
call results_closeGroup(results_addGroup(group))
sourceType: select case (phase_source(i,p)) if (phase_source(so,ph) /= SOURCE_UNDEFINED_ID) &
call results_closeGroup(results_addGroup(group//'sources/')) ! should be 'damage'
case (SOURCE_damage_anisoBrittle_ID) sourceType sourceType: select case (phase_source(so,ph))
call source_damage_anisoBrittle_results(p,group)
case (SOURCE_damage_anisoDuctile_ID) sourceType
call source_damage_anisoDuctile_results(p,group)
case (SOURCE_damage_isoBrittle_ID) sourceType
call source_damage_isoBrittle_results(p,group)
case (SOURCE_damage_isoDuctile_ID) sourceType
call source_damage_isoDuctile_results(p,group)
end select sourceType
enddo SourceLoop case (SOURCE_damage_anisoBrittle_ID) sourceType
enddo call source_damage_anisoBrittle_results(ph,group//'sources/')
case (SOURCE_damage_anisoDuctile_ID) sourceType
call source_damage_anisoDuctile_results(ph,group//'sources/')
case (SOURCE_damage_isoBrittle_ID) sourceType
call source_damage_isoBrittle_results(ph,group//'sources/')
case (SOURCE_damage_isoDuctile_ID) sourceType
call source_damage_isoDuctile_results(ph,group//'sources/')
end select sourceType
enddo SourceLoop
end subroutine damage_results end subroutine damage_results

897
src/constitutive_mech.f90
View File

File diff suppressed because it is too large Load Diff

10
src/constitutive_plastic_dislotwin.f90
View File

@ -485,12 +485,12 @@ end function plastic_dislotwin_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Return the homogenized elasticity matrix. !> @brief Return the homogenized elasticity matrix.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) module function plastic_dislotwin_homogenizedC(co,ip,el) result(homogenizedC)
real(pReal), dimension(6,6) :: & real(pReal), dimension(6,6) :: &
homogenizedC homogenizedC
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< component-ID of integration point co, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
@ -498,9 +498,9 @@ module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC)
of of
real(pReal) :: f_unrotated real(pReal) :: f_unrotated
of = material_phasememberAt(ipc,ip,el) of = material_phasememberAt(co,ip,el)
associate(prm => param(phase_plasticityInstance(material_phaseAt(ipc,el))),& associate(prm => param(phase_plasticityInstance(material_phaseAt(co,el))),&
stt => state(phase_plasticityInstance(material_phaseAT(ipc,el)))) stt => state(phase_plasticityInstance(material_phaseAT(co,el))))
f_unrotated = 1.0_pReal & f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,of)) & - sum(stt%f_tw(1:prm%sum_N_tw,of)) &

29
src/constitutive_plastic_nonlocal.f90
View File

@ -552,11 +552,10 @@ end function plastic_nonlocal_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates quantities characterizing the microstructure !> @brief calculates quantities characterizing the microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el) module subroutine plastic_nonlocal_dependentState(F, instance, of, ip, el)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
F, & F
Fp
integer, intent(in) :: & integer, intent(in) :: &
instance, & instance, &
of, & of, &
@ -564,6 +563,8 @@ module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el)
el el
integer :: & integer :: &
ph, &
me, &
no, & !< neighbor offset no, & !< neighbor offset
neighbor_el, & ! element number of neighboring material point neighbor_el, & ! element number of neighboring material point
neighbor_ip, & ! integration point of neighboring material point neighbor_ip, & ! integration point of neighboring material point
@ -643,8 +644,10 @@ module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el)
rho0 = getRho0(instance,of,ip,el) rho0 = getRho0(instance,of,ip,el)
if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then
invFp = math_inv33(Fp) ph = material_phaseAt(1,el)
invFe = matmul(Fp,math_inv33(F)) me = material_phaseMemberAt(1,ip,el)
invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me))
invFe = matmul(constitutive_mech_Fp(ph)%data(1:3,1:3,me),math_inv33(F))
rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg) rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg)
rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg) rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg)
@ -973,14 +976,13 @@ end subroutine plastic_nonlocal_deltaState
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure !> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & module subroutine plastic_nonlocal_dotState(Mp, F, Temperature,timestep, &
instance,of,ip,el) instance,of,ip,el)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress Mp !< MandelStress
real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: & real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: &
F, & !< elastic deformation gradient F !< Deformation gradient
Fp !< plastic deformation gradient
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
Temperature, & !< temperature Temperature, & !< temperature
timestep !< substepped crystallite time increment timestep !< substepped crystallite time increment
@ -1147,7 +1149,7 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
- rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
- rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have
rhoDot = rhoDotFlux(F,Fp,timestep, instance,of,ip,el) & rhoDot = rhoDotFlux(F,timestep, instance,of,ip,el) &
+ rhoDotMultiplication & + rhoDotMultiplication &
+ rhoDotSingle2DipoleGlide & + rhoDotSingle2DipoleGlide &
+ rhoDotAthermalAnnihilation & + rhoDotAthermalAnnihilation &
@ -1176,11 +1178,10 @@ end subroutine plastic_nonlocal_dotState
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure !> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
function rhoDotFlux(F,Fp,timestep, instance,of,ip,el) function rhoDotFlux(F,timestep, instance,of,ip,el)
real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: & real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: &
F, & !< elastic deformation gradient F !< Deformation gradient
Fp !< plastic deformation gradient
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timestep !< substepped crystallite time increment timestep !< substepped crystallite time increment
integer, intent(in) :: & integer, intent(in) :: &
@ -1293,7 +1294,7 @@ function rhoDotFlux(F,Fp,timestep, instance,of,ip,el)
m(1:3,:,4) = prm%slip_transverse m(1:3,:,4) = prm%slip_transverse
my_F = F(1:3,1:3,1,ip,el) my_F = F(1:3,1:3,1,ip,el)
my_Fe = matmul(my_F, math_inv33(Fp(1:3,1:3,1,ip,el))) my_Fe = matmul(my_F, math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,of)))
neighbors: do n = 1,nIPneighbors neighbors: do n = 1,nIPneighbors
@ -1311,7 +1312,7 @@ function rhoDotFlux(F,Fp,timestep, instance,of,ip,el)
if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient
neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el)) neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el))
neighbor_F = F(1:3,1:3,1,neighbor_ip,neighbor_el) neighbor_F = F(1:3,1:3,1,neighbor_ip,neighbor_el)
neighbor_Fe = matmul(neighbor_F, math_inv33(Fp(1:3,1:3,1,neighbor_ip,neighbor_el))) neighbor_Fe = matmul(neighbor_F, math_inv33(constitutive_mech_Fp(np)%data(1:3,1:3,no)))
Favg = 0.5_pReal * (my_F + neighbor_F) Favg = 0.5_pReal * (my_F + neighbor_F)
else ! if no neighbor, take my value as average else ! if no neighbor, take my value as average
Favg = my_F Favg = my_F

6
src/damage_nonlocal.f90
View File

@ -148,12 +148,12 @@ real(pReal) function damage_nonlocal_getMobility(ip,el)
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
integer :: & integer :: &
ipc co
damage_nonlocal_getMobility = 0.0_pReal damage_nonlocal_getMobility = 0.0_pReal
do ipc = 1, homogenization_Nconstituents(material_homogenizationAt(el)) do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_M(material_phaseAt(ipc,el)) damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_M(material_phaseAt(co,el))
enddo enddo
damage_nonlocal_getMobility = damage_nonlocal_getMobility/& damage_nonlocal_getMobility = damage_nonlocal_getMobility/&

250
src/homogenization.f90
View File

@ -48,20 +48,6 @@ module homogenization
type(tNumerics) :: num type(tNumerics) :: num
type :: tDebugOptions
logical :: &
basic, &
extensive, &
selective
integer :: &
element, &
ip, &
grain
end type tDebugOptions
type(tDebugOptions) :: debugHomog
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
interface interface
@ -112,6 +98,7 @@ module homogenization
public :: & public :: &
homogenization_init, & homogenization_init, &
materialpoint_stressAndItsTangent, & materialpoint_stressAndItsTangent, &
homogenization_forward, &
homogenization_results homogenization_results
contains contains
@ -124,24 +111,10 @@ subroutine homogenization_init
class (tNode) , pointer :: & class (tNode) , pointer :: &
num_homog, & num_homog, &
num_homogGeneric, & num_homogGeneric
debug_homogenization
print'(/,a)', ' <<<+- homogenization init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- homogenization init -+>>>'; flush(IO_STDOUT)
debug_homogenization => config_debug%get('homogenization', defaultVal=emptyList)
debugHomog%basic = debug_homogenization%contains('basic')
debugHomog%extensive = debug_homogenization%contains('extensive')
debugHomog%selective = debug_homogenization%contains('selective')
debugHomog%element = config_debug%get_asInt('element',defaultVal = 1)
debugHomog%ip = config_debug%get_asInt('integrationpoint',defaultVal = 1)
debugHomog%grain = config_debug%get_asInt('grain',defaultVal = 1)
if (debugHomog%grain < 1 &
.or. debugHomog%grain > homogenization_Nconstituents(material_homogenizationAt(debugHomog%element))) &
call IO_error(602,ext_msg='constituent', el=debugHomog%element, g=debugHomog%grain)
num_homog => config_numerics%get('homogenization',defaultVal=emptyDict) num_homog => config_numerics%get('homogenization',defaultVal=emptyDict)
num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict) num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict)
@ -177,172 +150,144 @@ subroutine materialpoint_stressAndItsTangent(dt)
integer :: & integer :: &
NiterationHomog, & NiterationHomog, &
NiterationMPstate, & NiterationMPstate, &
i, & !< integration point number ip, & !< integration point number
e, & !< element number el, & !< element number
myNgrains myNgrains, co, ce
real(pReal), dimension(discretization_nIPs,discretization_Nelems) :: & real(pReal) :: &
subFrac, & subFrac, &
subStep subStep
logical, dimension(discretization_nIPs,discretization_Nelems) :: & logical :: &
requested, & requested, &
converged converged
logical, dimension(2,discretization_nIPs,discretization_Nelems) :: & logical, dimension(2) :: &
doneAndHappy doneAndHappy
integer :: m
!$OMP PARALLEL DO PRIVATE(ce,myNgrains,NiterationMPstate,NiterationHomog,subFrac,converged,subStep,requested,doneAndHappy)
do el = FEsolving_execElem(1),FEsolving_execElem(2)
do ip = FEsolving_execIP(1),FEsolving_execIP(2)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize restoration points ! initialize restoration points
do e = FEsolving_execElem(1),FEsolving_execElem(2) call constitutive_initializeRestorationPoints(ip,el)
do i = FEsolving_execIP(1),FEsolving_execIP(2);
call crystallite_initializeRestorationPoints(i,e) subFrac = 0.0_pReal
converged = .false. ! pretend failed step ...
subStep = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation
requested = .true. ! everybody requires calculation
subFrac(i,e) = 0.0_pReal if (homogState(material_homogenizationAt(el))%sizeState > 0) &
converged(i,e) = .false. ! pretend failed step ... homogState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = &
subStep(i,e) = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation homogState(material_homogenizationAt(el))%State0( :,material_homogenizationMemberAt(ip,el))
requested(i,e) = .true. ! everybody requires calculation
if (homogState(material_homogenizationAt(e))%sizeState > 0) & if (damageState(material_homogenizationAt(el))%sizeState > 0) &
homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & damageState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = &
homogState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e)) damageState(material_homogenizationAt(el))%State0( :,material_homogenizationMemberAt(ip,el))
if (damageState(material_homogenizationAt(e))%sizeState > 0) & NiterationHomog = 0
damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & cutBackLooping: do while (.not. terminallyIll .and. subStep > num%subStepMinHomog)
damageState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e))
enddo
enddo
NiterationHomog = 0 myNgrains = homogenization_Nconstituents(material_homogenizationAt(el))
cutBackLooping: do while (.not. terminallyIll .and. & if (converged) then
any(subStep(FEsolving_execIP(1):FEsolving_execIP(2),& subFrac = subFrac + subStep
FEsolving_execElem(1):FEsolving_execElem(2)) > num%subStepMinHomog)) subStep = min(1.0_pReal-subFrac,num%stepIncreaseHomog*subStep) ! introduce flexibility for step increase/acceleration
!$OMP PARALLEL DO PRIVATE(m) steppingNeeded: if (subStep > num%subStepMinHomog) then
elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Nconstituents(material_homogenizationAt(e))
IpLooping1: do i = FEsolving_execIP(1),FEsolving_execIP(2)
if (converged(i,e)) then
subFrac(i,e) = subFrac(i,e) + subStep(i,e)
subStep(i,e) = min(1.0_pReal-subFrac(i,e),num%stepIncreaseHomog*subStep(i,e)) ! introduce flexibility for step increase/acceleration
steppingNeeded: if (subStep(i,e) > num%subStepMinHomog) then
! wind forward grain starting point ! wind forward grain starting point
call crystallite_windForward(i,e) call constitutive_windForward(ip,el)
if(homogState(material_homogenizationAt(e))%sizeState > 0) & if(homogState(material_homogenizationAt(el))%sizeState > 0) &
homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & homogState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = &
homogState(material_homogenizationAt(e))%State (:,material_homogenizationMemberAt(i,e)) homogState(material_homogenizationAt(el))%State (:,material_homogenizationMemberAt(ip,el))
if(damageState(material_homogenizationAt(e))%sizeState > 0) & if(damageState(material_homogenizationAt(el))%sizeState > 0) &
damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & damageState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = &
damageState(material_homogenizationAt(e))%State (:,material_homogenizationMemberAt(i,e)) damageState(material_homogenizationAt(el))%State (:,material_homogenizationMemberAt(ip,el))
endif steppingNeeded endif steppingNeeded
else else
if ( (myNgrains == 1 .and. subStep(i,e) <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite if ( (myNgrains == 1 .and. subStep <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite
num%subStepSizeHomog * subStep(i,e) <= num%subStepMinHomog ) then ! would require too small subStep num%subStepSizeHomog * subStep <= num%subStepMinHomog ) then ! would require too small subStep
! cutback makes no sense ! cutback makes no sense
if (.not. terminallyIll) then ! so first signals terminally ill... if (.not. terminallyIll) then ! so first signals terminally ill...
print*, ' Integration point ', i,' at element ', e, ' terminally ill' print*, ' Integration point ', ip,' at element ', el, ' terminally ill'
endif endif
terminallyIll = .true. ! ...and kills all others terminallyIll = .true. ! ...and kills all others
else ! cutback makes sense else ! cutback makes sense
subStep(i,e) = num%subStepSizeHomog * subStep(i,e) ! crystallite had severe trouble, so do a significant cutback subStep = num%subStepSizeHomog * subStep ! crystallite had severe trouble, so do a significant cutback
call crystallite_restore(i,e,subStep(i,e) < 1.0_pReal) call crystallite_restore(ip,el,subStep < 1.0_pReal)
call constitutive_restore(i,e) call constitutive_restore(ip,el)
if(homogState(material_homogenizationAt(e))%sizeState > 0) & if(homogState(material_homogenizationAt(el))%sizeState > 0) &
homogState(material_homogenizationAt(e))%State( :,material_homogenizationMemberAt(i,e)) = & homogState(material_homogenizationAt(el))%State( :,material_homogenizationMemberAt(ip,el)) = &
homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) homogState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el))
if(damageState(material_homogenizationAt(e))%sizeState > 0) & if(damageState(material_homogenizationAt(el))%sizeState > 0) &
damageState(material_homogenizationAt(e))%State( :,material_homogenizationMemberAt(i,e)) = & damageState(material_homogenizationAt(el))%State( :,material_homogenizationMemberAt(ip,el)) = &
damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) damageState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el))
endif endif
endif endif
if (subStep(i,e) > num%subStepMinHomog) then if (subStep > num%subStepMinHomog) then
requested(i,e) = .true. requested = .true.
doneAndHappy(1:2,i,e) = [.false.,.true.] doneAndHappy = [.false.,.true.]
endif endif
enddo IpLooping1
enddo elementLooping1
!$OMP END PARALLEL DO
NiterationMPstate = 0
convergenceLooping: do while (.not. terminallyIll .and. & NiterationMPstate = 0
any( requested(:,FEsolving_execELem(1):FEsolving_execElem(2)) & convergenceLooping: do while (.not. terminallyIll .and. requested &
.and. .not. doneAndHappy(1,:,FEsolving_execELem(1):FEsolving_execElem(2)) & .and. .not. doneAndHappy(1) &
) .and. & .and. NiterationMPstate < num%nMPstate)
NiterationMPstate < num%nMPstate) NiterationMPstate = NiterationMPstate + 1
NiterationMPstate = NiterationMPstate + 1
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! deformation partitioning ! deformation partitioning
!$OMP PARALLEL DO PRIVATE(myNgrains,m)
elementLooping2: do e = FEsolving_execElem(1),FEsolving_execElem(2) if(requested .and. .not. doneAndHappy(1)) then ! requested but not yet done
myNgrains = homogenization_Nconstituents(material_homogenizationAt(e)) ce = (el-1)*discretization_nIPs + ip
IpLooping2: do i = FEsolving_execIP(1),FEsolving_execIP(2) call mech_partition(homogenization_F0(1:3,1:3,ce) &
if(requested(i,e) .and. .not. doneAndHappy(1,i,e)) then ! requested but not yet done + (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce))&
m = (e-1)*discretization_nIPs + i *(subStep+subFrac), &
call mech_partition(homogenization_F0(1:3,1:3,m) & ip,el)
+ (homogenization_F(1:3,1:3,m)-homogenization_F0(1:3,1:3,m))& converged = .true.
*(subStep(i,e)+subFrac(i,e)), & do co = 1, myNgrains
i,e) converged = converged .and. crystallite_stress(dt*subStep,co,ip,el)
crystallite_dt(1:myNgrains,i,e) = dt*subStep(i,e) ! propagate materialpoint dt to grains enddo
crystallite_requested(1:myNgrains,i,e) = .true. ! request calculation for constituents
else
crystallite_requested(1:myNgrains,i,e) = .false. ! calculation for constituents not required anymore
endif endif
enddo IpLooping2
enddo elementLooping2
!$OMP END PARALLEL DO
!--------------------------------------------------------------------------------------------------
! crystallite integration
converged = crystallite_stress() !ToDo: MD not sure if that is the best logic
!-------------------------------------------------------------------------------------------------- if (requested .and. .not. doneAndHappy(1)) then
! state update if (.not. converged) then
!$OMP PARALLEL DO PRIVATE(m) doneAndHappy = [.true.,.false.]
elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2)
IpLooping3: do i = FEsolving_execIP(1),FEsolving_execIP(2)
if (requested(i,e) .and. .not. doneAndHappy(1,i,e)) then
if (.not. converged(i,e)) then
doneAndHappy(1:2,i,e) = [.true.,.false.]
else else
m = (e-1)*discretization_nIPs + i ce = (el-1)*discretization_nIPs + ip
doneAndHappy(1:2,i,e) = updateState(dt*subStep(i,e), & doneAndHappy = updateState(dt*subStep, &
homogenization_F0(1:3,1:3,m) & homogenization_F0(1:3,1:3,ce) &
+ (homogenization_F(1:3,1:3,m)-homogenization_F0(1:3,1:3,m)) & + (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce)) &
*(subStep(i,e)+subFrac(i,e)), & *(subStep+subFrac), &
i,e) ip,el)
converged(i,e) = all(doneAndHappy(1:2,i,e)) ! converged if done and happy converged = all(doneAndHappy)
endif endif
endif endif
enddo IpLooping3
enddo elementLooping3
!$OMP END PARALLEL DO
enddo convergenceLooping enddo convergenceLooping
NiterationHomog = NiterationHomog + 1
NiterationHomog = NiterationHomog + 1 enddo cutBackLooping
enddo
enddo cutBackLooping enddo
!$OMP END PARALLEL DO
if (.not. terminallyIll ) then if (.not. terminallyIll ) then
call crystallite_orientations() ! calculate crystal orientations call crystallite_orientations() ! calculate crystal orientations
!$OMP PARALLEL DO !$OMP PARALLEL DO
elementLooping4: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
IpLooping4: do i = FEsolving_execIP(1),FEsolving_execIP(2) IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2)
call mech_homogenize(i,e) call mech_homogenize(ip,el)
enddo IpLooping4 enddo IpLooping3
enddo elementLooping4 enddo elementLooping3
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
else else
print'(/,a,/)', ' << HOMOG >> Material Point terminally ill' print'(/,a,/)', ' << HOMOG >> Material Point terminally ill'
@ -399,6 +344,7 @@ subroutine homogenization_results
integer :: p integer :: p
character(len=:), allocatable :: group_base,group character(len=:), allocatable :: group_base,group
call results_closeGroup(results_addGroup('current/homogenization/'))
do p=1,size(material_name_homogenization) do p=1,size(material_name_homogenization)
group_base = 'current/homogenization/'//trim(material_name_homogenization(p)) group_base = 'current/homogenization/'//trim(material_name_homogenization(p))
@ -424,4 +370,20 @@ subroutine homogenization_results
end subroutine homogenization_results end subroutine homogenization_results
!--------------------------------------------------------------------------------------------------
!> @brief Forward data after successful increment.
! ToDo: Any guessing for the current states possible?
!--------------------------------------------------------------------------------------------------
subroutine homogenization_forward
integer :: ho
do ho = 1, size(material_name_homogenization)
homogState (ho)%state0 = homogState (ho)%state
damageState(ho)%state0 = damageState(ho)%state
enddo
end subroutine homogenization_forward
end module homogenization end module homogenization

1
src/homogenization_mech.f90
View File

@ -4,6 +4,7 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(homogenization) homogenization_mech submodule(homogenization) homogenization_mech
interface interface
module subroutine mech_none_init module subroutine mech_none_init

178
src/homogenization_mech_RGC.f90
View File

@ -18,8 +18,6 @@ submodule(homogenization:homogenization_mech) homogenization_mech_RGC
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
D_alpha, & D_alpha, &
a_g a_g
integer :: &
of_debug = 0
character(len=pStringLen), allocatable, dimension(:) :: & character(len=pStringLen), allocatable, dimension(:) :: &
output output
end type tParameters end type tParameters
@ -151,12 +149,6 @@ module subroutine mech_RGC_init(num_homogMech)
st0 => state0(homogenization_typeInstance(h)), & st0 => state0(homogenization_typeInstance(h)), &
dst => dependentState(homogenization_typeInstance(h))) dst => dependentState(homogenization_typeInstance(h)))
#ifdef DEBUG
if (h==material_homogenizationAt(debugHomog%element)) then
prm%of_debug = material_homogenizationMemberAt(debugHomog%ip,debugHomog%element)
endif
#endif
#if defined (__GFORTRAN__) #if defined (__GFORTRAN__)
prm%output = output_asStrings(homogMech) prm%output = output_asStrings(homogMech)
#else #else
@ -239,17 +231,6 @@ module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
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
enddo enddo
F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! resulting relaxed deformation gradient F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! resulting relaxed deformation gradient
#ifdef DEBUG
if (debugHomog%extensive) then
print'(a,i3)',' Deformation gradient of grain: ',iGrain
do i = 1,3
print'(1x,3(e15.8,1x))',(F(i,j,iGrain), j = 1,3)
enddo
print*,' '
flush(IO_STDOUT)
endif
#endif
enddo enddo
end associate end associate
@ -273,10 +254,6 @@ module procedure mech_RGC_updateState
logical :: error logical :: error
real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix
real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax
#ifdef DEBUG
integer, dimension(3) :: stresLoc
integer, dimension(2) :: residLoc
#endif
zeroTimeStep: if(dEq0(dt)) then zeroTimeStep: if(dEq0(dt)) then
mech_RGC_updateState = .true. ! pretend everything is fine and return mech_RGC_updateState = .true. ! pretend everything is fine and return
@ -303,16 +280,6 @@ module procedure mech_RGC_updateState
relax = stt%relaxationVector(:,of) relax = stt%relaxationVector(:,of)
drelax = stt%relaxationVector(:,of) - st0%relaxationVector(:,of) drelax = stt%relaxationVector(:,of) - st0%relaxationVector(:,of)
#ifdef DEBUG
if (debugHomog%extensive) then
print*, 'Obtained state: '
do i = 1,size(stt%relaxationVector(:,of))
print'(1x,2(e15.8,1x))', stt%relaxationVector(i,of)
enddo
print*,' '
endif
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! computing interface mismatch and stress penalty tensor for all interfaces of all grains ! computing interface mismatch and stress penalty tensor for all interfaces of all grains
call stressPenalty(R,NN,avgF,F,ip,el,instance,of) call stressPenalty(R,NN,avgF,F,ip,el,instance,of)
@ -353,13 +320,6 @@ module procedure mech_RGC_updateState
enddo enddo
enddo enddo
#ifdef DEBUG
if (debugHomog%extensive) then
print'(a,i3)',' Traction at interface: ',iNum
print'(1x,3(e15.8,1x))',(tract(iNum,j), j = 1,3)
print*,' '
endif
#endif
enddo enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -367,29 +327,12 @@ module procedure mech_RGC_updateState
stresMax = maxval(abs(P)) ! get the maximum of first Piola-Kirchhoff (material) stress stresMax = maxval(abs(P)) ! get the maximum of first Piola-Kirchhoff (material) stress
residMax = maxval(abs(tract)) ! get the maximum of the residual residMax = maxval(abs(tract)) ! get the maximum of the residual
#ifdef DEBUG
if (debugHomog%extensive .and. prm%of_debug == of) then
stresLoc = maxloc(abs(P))
residLoc = maxloc(abs(tract))
print'(a,i2,1x,i4)',' RGC residual check ... ',ip,el
print'(a,e15.8,a,i3,a,i2,i2)', ' Max stress: ',stresMax, &
'@ grain ',stresLoc(3),' in component ',stresLoc(1),stresLoc(2)
print'(a,e15.8,a,i3,a,i2)',' Max residual: ',residMax, &
' @ iface ',residLoc(1),' in direction ',residLoc(2)
flush(IO_STDOUT)
endif
#endif
mech_RGC_updateState = .false. mech_RGC_updateState = .false.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! If convergence reached => done and happy ! If convergence reached => done and happy
if (residMax < num%rtol*stresMax .or. residMax < num%atol) then if (residMax < num%rtol*stresMax .or. residMax < num%atol) then
mech_RGC_updateState = .true. mech_RGC_updateState = .true.
#ifdef DEBUG
if (debugHomog%extensive .and. prm%of_debug == of) &
print*, '... done and happy'; flush(IO_STDOUT)
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 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
@ -406,41 +349,14 @@ module procedure mech_RGC_updateState
dst%relaxationRate_avg(of) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal) dst%relaxationRate_avg(of) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal)
dst%relaxationRate_max(of) = maxval(abs(drelax))/dt dst%relaxationRate_max(of) = maxval(abs(drelax))/dt
#ifdef DEBUG
if (debugHomog%extensive .and. prm%of_debug == of) then
print'(a,e15.8)', ' Constitutive work: ',stt%work(of)
print'(a,3(1x,e15.8))', ' Magnitude mismatch: ',dst%mismatch(1,of), &
dst%mismatch(2,of), &
dst%mismatch(3,of)
print'(a,e15.8)', ' Penalty energy: ', stt%penaltyEnergy(of)
print'(a,e15.8,/)', ' Volume discrepancy: ', dst%volumeDiscrepancy(of)
print'(a,e15.8)', ' Maximum relaxation rate: ', dst%relaxationRate_max(of)
print'(a,e15.8,/)', ' Average relaxation rate: ', dst%relaxationRate_avg(of)
flush(IO_STDOUT)
endif
#endif
return return
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! if residual blows-up => done but unhappy ! if residual blows-up => done but unhappy
elseif (residMax > num%relMax*stresMax .or. residMax > num%absMax) then ! try to restart when residual blows up exceeding maximum bound elseif (residMax > num%relMax*stresMax .or. residMax > num%absMax) then ! try to restart when residual blows up exceeding maximum bound
mech_RGC_updateState = [.true.,.false.] ! with direct cut-back mech_RGC_updateState = [.true.,.false.] ! with direct cut-back
#ifdef DEBUG
if (debugHomog%extensive .and. prm%of_debug == of) &
print'(a,/)', ' ... broken'; flush(IO_STDOUT)
#endif
return return
endif
else ! proceed with computing the Jacobian and state update
#ifdef DEBUG
if (debugHomog%extensive .and. prm%of_debug == of) &
print'(a,/)', ' ... not yet done'; flush(IO_STDOUT)
#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
@ -492,17 +408,6 @@ module procedure mech_RGC_updateState
enddo enddo
enddo enddo
#ifdef DEBUG
if (debugHomog%extensive) then
print*, 'Jacobian matrix of stress'
do i = 1,3*nIntFaceTot
print'(1x,100(e11.4,1x))',(smatrix(i,j), j = 1,3*nIntFaceTot)
enddo
print*,' '
flush(IO_STDOUT)
endif
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! ... of the stress penalty tangent (mismatch penalty and volume penalty, computed using numerical ! ... of the stress penalty tangent (mismatch penalty and volume penalty, computed using numerical
! perturbation method) "pmatrix" ! perturbation method) "pmatrix"
@ -552,16 +457,6 @@ module procedure mech_RGC_updateState
pmatrix(:,ipert) = p_resid/num%pPert pmatrix(:,ipert) = p_resid/num%pPert
enddo enddo
#ifdef DEBUG
if (debugHomog%extensive) then
print*, 'Jacobian matrix of penalty'
do i = 1,3*nIntFaceTot
print'(1x,100(e11.4,1x))',(pmatrix(i,j), j = 1,3*nIntFaceTot)
enddo
print*,' '
flush(IO_STDOUT)
endif
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! ... of the numerical viscosity traction "rmatrix" ! ... of the numerical viscosity traction "rmatrix"
@ -571,48 +466,16 @@ module procedure mech_RGC_updateState
(abs(drelax(i))/(num%refRelaxRate*dt))**(num%viscPower - 1.0_pReal) ! only in the main diagonal term (abs(drelax(i))/(num%refRelaxRate*dt))**(num%viscPower - 1.0_pReal) ! only in the main diagonal term
enddo enddo
#ifdef DEBUG
if (debugHomog%extensive) then
print*, 'Jacobian matrix of penalty'
do i = 1,3*nIntFaceTot
print'(1x,100(e11.4,1x))',(rmatrix(i,j), j = 1,3*nIntFaceTot)
enddo
print*,' '
flush(IO_STDOUT)
endif
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! The overall Jacobian matrix summarizing contributions of smatrix, pmatrix, rmatrix ! The overall Jacobian matrix summarizing contributions of smatrix, pmatrix, rmatrix
allocate(jmatrix(3*nIntFaceTot,3*nIntFaceTot)); jmatrix = smatrix + pmatrix + rmatrix allocate(jmatrix(3*nIntFaceTot,3*nIntFaceTot)); jmatrix = smatrix + pmatrix + rmatrix
#ifdef DEBUG
if (debugHomog%extensive) then
print*, 'Jacobian matrix (total)'
do i = 1,3*nIntFaceTot
print'(1x,100(e11.4,1x))',(jmatrix(i,j), j = 1,3*nIntFaceTot)
enddo
print*,' '
flush(IO_STDOUT)
endif
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! computing the update of the state variable (relaxation vectors) using the Jacobian matrix ! computing the update of the state variable (relaxation vectors) using the Jacobian matrix
allocate(jnverse(3*nIntFaceTot,3*nIntFaceTot),source=0.0_pReal) allocate(jnverse(3*nIntFaceTot,3*nIntFaceTot),source=0.0_pReal)
call math_invert(jnverse,error,jmatrix) call math_invert(jnverse,error,jmatrix)
#ifdef DEBUG
if (debugHomog%extensive) then
print*, 'Jacobian inverse'
do i = 1,3*nIntFaceTot
print'(1x,100(e11.4,1x))',(jnverse(i,j), j = 1,3*nIntFaceTot)
enddo
print*,' '
flush(IO_STDOUT)
endif
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate the state update (global relaxation vectors) for the next Newton-Raphson iteration ! calculate the state update (global relaxation vectors) for the next Newton-Raphson iteration
drelax = 0.0_pReal drelax = 0.0_pReal
@ -629,17 +492,6 @@ module procedure mech_RGC_updateState
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
#ifdef DEBUG
if (debugHomog%extensive) then
print*, 'Returned state: '
do i = 1,size(stt%relaxationVector(:,of))
print'(1x,2(e15.8,1x))', stt%relaxationVector(i,of)
enddo
print*,' '
flush(IO_STDOUT)
endif
#endif
end associate end associate
contains contains
@ -676,12 +528,6 @@ module procedure mech_RGC_updateState
associate(prm => param(instance)) associate(prm => param(instance))
#ifdef DEBUG
if (debugHomog%extensive .and. prm%of_debug == of) then
print'(a,2(1x,i3))', ' Correction factor: ',ip,el
print*, surfCorr
endif
#endif
!----------------------------------------------------------------------------------------------- !-----------------------------------------------------------------------------------------------
! computing the mismatch and penalty stress tensor of all grains ! computing the mismatch and penalty stress tensor of all grains
@ -717,13 +563,7 @@ module procedure mech_RGC_updateState
enddo; enddo enddo; enddo
nDefNorm = max(nDefToler,sqrt(nDefNorm)) ! approximation to zero mismatch if mismatch is zero (singularity) nDefNorm = max(nDefToler,sqrt(nDefNorm)) ! approximation to zero mismatch if mismatch is zero (singularity)
nMis(abs(intFace(1)),iGrain) = nMis(abs(intFace(1)),iGrain) + nDefNorm ! total amount of mismatch experienced by the grain (at all six interfaces) nMis(abs(intFace(1)),iGrain) = nMis(abs(intFace(1)),iGrain) + nDefNorm ! total amount of mismatch experienced by the grain (at all six interfaces)
#ifdef DEBUG
if (debugHomog%extensive .and. prm%of_debug == of) then
print'(a,i2,a,i3)',' Mismatch to face: ',intFace(1),' neighbor grain: ',iGNghb
print*, transpose(nDef)
print'(a,e11.4)', ' with magnitude: ',nDefNorm
endif
#endif
!------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------
! compute the stress penalty of all interfaces ! compute the stress penalty of all interfaces
@ -735,12 +575,7 @@ module procedure mech_RGC_updateState
*tanh(nDefNorm/num%xSmoo) *tanh(nDefNorm/num%xSmoo)
enddo; enddo;enddo; enddo enddo; enddo;enddo; enddo
enddo interfaceLoop enddo interfaceLoop
#ifdef DEBUG
if (debugHomog%extensive .and. prm%of_debug == of) then
print'(a,i2)', ' Penalty of grain: ',iGrain
print*, transpose(rPen(1:3,1:3,iGrain))
endif
#endif
enddo grainLoop enddo grainLoop
@ -783,13 +618,6 @@ module procedure mech_RGC_updateState
vPen(:,:,i) = -1.0_pReal/real(nGrain,pReal)*num%volDiscrMod*num%volDiscrPow/num%maxVolDiscr* & vPen(:,:,i) = -1.0_pReal/real(nGrain,pReal)*num%volDiscrMod*num%volDiscrPow/num%maxVolDiscr* &
sign((abs(vDiscrep)/num%maxVolDiscr)**(num%volDiscrPow - 1.0),vDiscrep)* & sign((abs(vDiscrep)/num%maxVolDiscr)**(num%volDiscrPow - 1.0),vDiscrep)* &
gVol(i)*transpose(math_inv33(fDef(:,:,i))) gVol(i)*transpose(math_inv33(fDef(:,:,i)))
#ifdef DEBUG
if (debugHomog%extensive .and. param(instance)%of_debug == of) then
print'(a,i2)',' Volume penalty of grain: ',i
print*, transpose(vPen(:,:,i))
endif
#endif
enddo enddo
end subroutine volumePenalty end subroutine volumePenalty

6
src/kinematics_cleavage_opening.f90
View File

@ -99,10 +99,10 @@ end function kinematics_cleavage_opening_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< grain number co, & !< grain number
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
@ -124,7 +124,7 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S,
Ld = 0.0_pReal Ld = 0.0_pReal
dLd_dTstar = 0.0_pReal dLd_dTstar = 0.0_pReal
associate(prm => param(kinematics_cleavage_opening_instance(material_phaseAt(ipc,el)))) associate(prm => param(kinematics_cleavage_opening_instance(material_phaseAt(co,el))))
do i = 1,prm%sum_N_cl do i = 1,prm%sum_N_cl
traction_crit = prm%g_crit(i)* damage(homog)%p(damageOffset)**2.0_pReal traction_crit = prm%g_crit(i)* damage(homog)%p(damageOffset)**2.0_pReal

6
src/kinematics_slipplane_opening.f90
View File

@ -117,10 +117,10 @@ end function kinematics_slipplane_opening_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< grain number co, & !< grain number
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
@ -138,7 +138,7 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S
traction_d, traction_t, traction_n, traction_crit, & traction_d, traction_t, traction_n, traction_crit, &
udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
phase = material_phaseAt(ipc,el) phase = material_phaseAt(co,el)
instance = kinematics_slipplane_opening_instance(phase) instance = kinematics_slipplane_opening_instance(phase)
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
damageOffset = material_homogenizationMemberAt(ip,el) damageOffset = material_homogenizationMemberAt(ip,el)

6
src/kinematics_thermal_expansion.f90
View File

@ -84,10 +84,10 @@ end function kinematics_thermal_expansion_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief constitutive equation for calculating the velocity gradient !> @brief constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el) module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, co, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< grain number co, & !< grain number
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), intent(out), dimension(3,3) :: & real(pReal), intent(out), dimension(3,3) :: &
@ -101,7 +101,7 @@ module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, i
real(pReal) :: & real(pReal) :: &
T, TDot T, TDot
phase = material_phaseAt(ipc,el) phase = material_phaseAt(co,el)
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
T = temperature(homog)%p(material_homogenizationMemberAt(ip,el)) T = temperature(homog)%p(material_homogenizationMemberAt(ip,el))
TDot = temperatureRate(homog)%p(material_homogenizationMemberAt(ip,el)) TDot = temperatureRate(homog)%p(material_homogenizationMemberAt(ip,el))

46
src/material.f90
View File

@ -17,29 +17,6 @@ module material
private private
enum, bind(c); enumerator :: & enum, bind(c); enumerator :: &
ELASTICITY_UNDEFINED_ID, &
ELASTICITY_HOOKE_ID, &
PLASTICITY_UNDEFINED_ID, &
PLASTICITY_NONE_ID, &
PLASTICITY_ISOTROPIC_ID, &
PLASTICITY_PHENOPOWERLAW_ID, &
PLASTICITY_KINEHARDENING_ID, &
PLASTICITY_DISLOTWIN_ID, &
PLASTICITY_DISLOTUNGSTEN_ID, &
PLASTICITY_NONLOCAL_ID, &
SOURCE_UNDEFINED_ID ,&
SOURCE_THERMAL_DISSIPATION_ID, &
SOURCE_THERMAL_EXTERNALHEAT_ID, &
SOURCE_DAMAGE_ISOBRITTLE_ID, &
SOURCE_DAMAGE_ISODUCTILE_ID, &
SOURCE_DAMAGE_ANISOBRITTLE_ID, &
SOURCE_DAMAGE_ANISODUCTILE_ID, &
KINEMATICS_UNDEFINED_ID ,&
KINEMATICS_CLEAVAGE_OPENING_ID, &
KINEMATICS_SLIPPLANE_OPENING_ID, &
KINEMATICS_THERMAL_EXPANSION_ID, &
STIFFNESS_DEGRADATION_UNDEFINED_ID, &
STIFFNESS_DEGRADATION_DAMAGE_ID, &
THERMAL_ISOTHERMAL_ID, & THERMAL_ISOTHERMAL_ID, &
THERMAL_CONDUCTION_ID, & THERMAL_CONDUCTION_ID, &
DAMAGE_NONE_ID, & DAMAGE_NONE_ID, &
@ -96,29 +73,6 @@ module material
public :: & public :: &
material_init, & material_init, &
ELASTICITY_UNDEFINED_ID, &
ELASTICITY_HOOKE_ID, &
PLASTICITY_UNDEFINED_ID, &
PLASTICITY_NONE_ID, &
PLASTICITY_ISOTROPIC_ID, &
PLASTICITY_PHENOPOWERLAW_ID, &
PLASTICITY_KINEHARDENING_ID, &
PLASTICITY_DISLOTWIN_ID, &
PLASTICITY_DISLOTUNGSTEN_ID, &
PLASTICITY_NONLOCAL_ID, &
SOURCE_UNDEFINED_ID ,&
SOURCE_THERMAL_DISSIPATION_ID, &
SOURCE_THERMAL_EXTERNALHEAT_ID, &
SOURCE_DAMAGE_ISOBRITTLE_ID, &
SOURCE_DAMAGE_ISODUCTILE_ID, &
SOURCE_DAMAGE_ANISOBRITTLE_ID, &
SOURCE_DAMAGE_ANISODUCTILE_ID, &
KINEMATICS_UNDEFINED_ID ,&
KINEMATICS_CLEAVAGE_OPENING_ID, &
KINEMATICS_SLIPPLANE_OPENING_ID, &
KINEMATICS_THERMAL_EXPANSION_ID, &
STIFFNESS_DEGRADATION_UNDEFINED_ID, &
STIFFNESS_DEGRADATION_DAMAGE_ID, &
THERMAL_ISOTHERMAL_ID, & THERMAL_ISOTHERMAL_ID, &
THERMAL_CONDUCTION_ID, & THERMAL_CONDUCTION_ID, &
DAMAGE_NONE_ID, & DAMAGE_NONE_ID, &

2
src/results.f90
View File

@ -111,8 +111,6 @@ subroutine results_addIncrement(inc,time)
call results_closeGroup(results_addGroup(trim('inc'//trim(adjustl(incChar))))) call results_closeGroup(results_addGroup(trim('inc'//trim(adjustl(incChar)))))
call results_setLink(trim('inc'//trim(adjustl(incChar))),'current') call results_setLink(trim('inc'//trim(adjustl(incChar))),'current')
call results_addAttribute('time/s',time,trim('inc'//trim(adjustl(incChar)))) call results_addAttribute('time/s',time,trim('inc'//trim(adjustl(incChar))))
call results_closeGroup(results_addGroup('current/phase'))
call results_closeGroup(results_addGroup('current/homogenization'))
end subroutine results_addIncrement end subroutine results_addIncrement

8
src/source_damage_anisoBrittle.f90
View File

@ -120,10 +120,10 @@ end function source_damage_anisoBrittle_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state !> @brief calculates derived quantities from state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< component-ID of integration point co, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
@ -139,8 +139,8 @@ module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
real(pReal) :: & real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit traction_d, traction_t, traction_n, traction_crit
phase = material_phaseAt(ipc,el) phase = material_phaseAt(co,el)
constituent = material_phasememberAt(ipc,ip,el) constituent = material_phasememberAt(co,ip,el)
sourceOffset = source_damage_anisoBrittle_offset(phase) sourceOffset = source_damage_anisoBrittle_offset(phase)
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
damageOffset = material_homogenizationMemberAt(ip,el) damageOffset = material_homogenizationMemberAt(ip,el)

8
src/source_damage_anisoDuctile.f90
View File

@ -107,10 +107,10 @@ end function source_damage_anisoDuctile_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state !> @brief calculates derived quantities from state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) module subroutine source_damage_anisoDuctile_dotState(co, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< component-ID of integration point co, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
@ -121,8 +121,8 @@ module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el)
damageOffset, & damageOffset, &
homog homog
phase = material_phaseAt(ipc,el) phase = material_phaseAt(co,el)
constituent = material_phasememberAt(ipc,ip,el) constituent = material_phasememberAt(co,ip,el)
sourceOffset = source_damage_anisoDuctile_offset(phase) sourceOffset = source_damage_anisoDuctile_offset(phase)
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
damageOffset = material_homogenizationMemberAt(ip,el) damageOffset = material_homogenizationMemberAt(ip,el)

8
src/source_damage_isoBrittle.f90
View File

@ -94,10 +94,10 @@ end function source_damage_isoBrittle_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state !> @brief calculates derived quantities from state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) module subroutine source_damage_isoBrittle_deltaState(C, Fe, co, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< component-ID of integration point co, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
@ -114,8 +114,8 @@ module subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
real(pReal) :: & real(pReal) :: &
strainenergy strainenergy
phase = material_phaseAt(ipc,el) !< phase ID at ipc,ip,el phase = material_phaseAt(co,el) !< phase ID at co,ip,el
constituent = material_phasememberAt(ipc,ip,el) !< state array offset for phase ID at ipc,ip,el constituent = material_phasememberAt(co,ip,el) !< state array offset for phase ID at co,ip,el
sourceOffset = source_damage_isoBrittle_offset(phase) sourceOffset = source_damage_isoBrittle_offset(phase)
strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3) strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3)

8
src/source_damage_isoDuctile.f90
View File

@ -98,10 +98,10 @@ end function source_damage_isoDuctile_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state !> @brief calculates derived quantities from state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine source_damage_isoDuctile_dotState(ipc, ip, el) module subroutine source_damage_isoDuctile_dotState(co, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< component-ID of integration point co, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
@ -112,8 +112,8 @@ module subroutine source_damage_isoDuctile_dotState(ipc, ip, el)
damageOffset, & damageOffset, &
homog homog
phase = material_phaseAt(ipc,el) phase = material_phaseAt(co,el)
constituent = material_phasememberAt(ipc,ip,el) constituent = material_phasememberAt(co,ip,el)
sourceOffset = source_damage_isoDuctile_offset(phase) sourceOffset = source_damage_isoDuctile_offset(phase)
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
damageOffset = material_homogenizationMemberAt(ip,el) damageOffset = material_homogenizationMemberAt(ip,el)