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

the latest RGC model + corrections for "element homogeneous" feature

This commit is contained in:
Denny Tjahjanto 2010-03-24 13:20:12 +00:00
parent 3aa2dd5fef
commit 40b1478dac
13 changed files with 364 additions and 162 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
code/CPFEM.f90
View File

@ -86,7 +86,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
cycleCounter, & cycleCounter, &
theInc, & theInc, &
theTime, & theTime, &
theDelta, &
FEsolving_execElem, & FEsolving_execElem, &
FEsolving_execIP FEsolving_execIP
use math, only: math_init, & use math, only: math_init, &
@ -102,11 +101,14 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
mesh_maxNips mesh_maxNips
use lattice, only: lattice_init use lattice, only: lattice_init
use material, only: material_init, & use material, only: material_init, &
homogenization_maxNgrains homogenization_maxNgrains, &
microstructure_elemhomo
use constitutive, only: constitutive_init,& use constitutive, only: constitutive_init,&
constitutive_state0,constitutive_state constitutive_state0,constitutive_state
use crystallite, only: crystallite_init, & use crystallite, only: crystallite_init, &
crystallite_F0, & crystallite_F0, &
crystallite_rexParm, &
crystallite_critVal, &
crystallite_partionedF, & crystallite_partionedF, &
crystallite_Fp0, & crystallite_Fp0, &
crystallite_Fp, & crystallite_Fp, &
@ -122,6 +124,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
materialpoint_F0, & materialpoint_F0, &
materialpoint_P, & materialpoint_P, &
materialpoint_dPdF, & materialpoint_dPdF, &
materialpoint_results, &
materialpoint_Temperature, & materialpoint_Temperature, &
materialpoint_stressAndItsTangent, & materialpoint_stressAndItsTangent, &
materialpoint_postResults materialpoint_postResults
@ -286,6 +289,16 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
elseif (.not. CPFEM_calc_done) then elseif (.not. CPFEM_calc_done) then
call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent (parallel execution inside) call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent (parallel execution inside)
call materialpoint_postResults(dt) ! post results call materialpoint_postResults(dt) ! post results
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! loop over all parallely processed elements
if (microstructure_elemhomo(mesh_element(4,e))) then ! dealing with homogeneous element?
forall (i = 2:FE_Nips(mesh_element(2,e))) ! copy results of first IP to all others
materialpoint_P(:,:,i,e) = materialpoint_P(:,:,1,e)
materialpoint_F(:,:,i,e) = materialpoint_F(:,:,1,e)
materialpoint_dPdF(:,:,:,:,i,e) = materialpoint_dPdF(:,:,:,:,1,e)
materialpoint_results(:,i,e) = materialpoint_results(:,1,e)
end forall
endif
enddo
CPFEM_calc_done = .true. CPFEM_calc_done = .true.
endif endif
@ -321,8 +334,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
else if (mode == 5) then else if (mode == 5) then
CPFEM_dcsde = CPFEM_dcsde_knownGood ! --+>> RESTORE CONSISTENT JACOBIAN FROM FORMER CONVERGED INC CPFEM_dcsde = CPFEM_dcsde_knownGood ! --+>> RESTORE CONSISTENT JACOBIAN FROM FORMER CONVERGED INC
end if end if
call random_number(rnd) call random_number(rnd)
if (rnd < 0.5_pReal) rnd = 1.0_pReal - rnd if (rnd < 0.5_pReal) rnd = 1.0_pReal - rnd
materialpoint_Temperature(IP,cp_en) = Temperature materialpoint_Temperature(IP,cp_en) = Temperature
materialpoint_F0(:,:,IP,cp_en) = ffn materialpoint_F0(:,:,IP,cp_en) = ffn
materialpoint_F(:,:,IP,cp_en) = ffn1 materialpoint_F(:,:,IP,cp_en) = ffn1

2
code/FEsolving.f90
View File

@ -1,4 +1,4 @@
!* $Id$ !* 'Id'
!############################################################## !##############################################################
MODULE FEsolving MODULE FEsolving
!############################################################## !##############################################################

45
code/constitutive.f90
View File

@ -32,6 +32,7 @@ CONTAINS
!**************************************** !****************************************
!* - constitutive_init !* - constitutive_init
!* - constitutive_homogenizedC !* - constitutive_homogenizedC
!* - constitutive_averageBurgers
!* - constitutive_microstructure !* - constitutive_microstructure
!* - constitutive_LpAndItsTangent !* - constitutive_LpAndItsTangent
!* - constitutive_collectDotState !* - constitutive_collectDotState
@ -268,8 +269,48 @@ function constitutive_homogenizedC(ipc,ip,el)
return return
endfunction endfunction
function constitutive_averageBurgers(ipc,ip,el)
!*********************************************************************
!* This function returns the average length of Burgers vector *
!* INPUT: *
!* - state : state variables *
!* - ipc : component-ID of current integration point *
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt
use material, only: phase_constitution,material_phase
use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_dislotwin
use constitutive_nonlocal
implicit none
subroutine constitutive_microstructure(Temperature,Tstar_v,Fe,Fp,ipc,ip,el) !* Definition of variables
integer(pInt) ipc,ip,el
real(pReal) :: constitutive_averageBurgers
select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_j2_label)
constitutive_averageBurgers = 2.5e-10_pReal !constitutive_j2_averageBurgers(constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label)
constitutive_averageBurgers = 2.5e-10_pReal !constitutive_phenopowerlaw_averageBurgers(constitutive_state,ipc,ip,el)
case (constitutive_dislotwin_label)
constitutive_averageBurgers = 2.5e-10_pReal !constitutive_dislotwin_averageBurgers(constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label)
constitutive_averageBurgers = 2.5e-10_pReal !constitutive_nonlocal_averageBurgers(constitutive_state,ipc,ip,el)
end select
return
endfunction
subroutine constitutive_microstructure(Temperature,Fe,Fp,ipc,ip,el)
!********************************************************************* !*********************************************************************
!* This function calculates from state needed variables * !* This function calculates from state needed variables *
!* INPUT: * !* INPUT: *
@ -538,4 +579,4 @@ return
endfunction endfunction
END MODULE END MODULE

17
code/constitutive_nonlocal.f90
View File

@ -73,10 +73,10 @@ real(pReal), dimension(:,:), allocatable :: constitutive_nonlocal_
constitutive_nonlocal_lambda0PerSlipSystem, & ! mean free path prefactor for each slip system and instance constitutive_nonlocal_lambda0PerSlipSystem, & ! mean free path prefactor for each slip system and instance
constitutive_nonlocal_burgersPerSlipFamily, & ! absolute length of burgers vector [m] for each family and instance constitutive_nonlocal_burgersPerSlipFamily, & ! absolute length of burgers vector [m] for each family and instance
constitutive_nonlocal_burgersPerSlipSystem, & ! absolute length of burgers vector [m] for each slip system and instance constitutive_nonlocal_burgersPerSlipSystem, & ! absolute length of burgers vector [m] for each slip system and instance
constitutive_nonlocal_dLowerEdgePerSlipFamily, & ! minimum stable edge dipole height for each family and instance constitutive_nonlocal_dLowerEdgePerSlipFamily, & ! minimum stable edge dipole height for each family and instance
constitutive_nonlocal_dLowerEdgePerSlipSystem, & ! minimum stable edge dipole height for each slip system and instance constitutive_nonlocal_dLowerEdgePerSlipSystem, & ! minimum stable edge dipole height for each slip system and instance
constitutive_nonlocal_dLowerScrewPerSlipFamily, & ! minimum stable screw dipole height for each family and instance constitutive_nonlocal_dLowerScrewPerSlipFamily, & ! minimum stable screw dipole height for each family and instance
constitutive_nonlocal_dLowerScrewPerSlipSystem, & ! minimum stable screw dipole height for each slip system and instance constitutive_nonlocal_dLowerScrewPerSlipSystem, & ! minimum stable screw dipole height for each slip system and instance
constitutive_nonlocal_interactionSlipSlip ! coefficients for slip-slip interaction for each interaction type and instance constitutive_nonlocal_interactionSlipSlip ! coefficients for slip-slip interaction for each interaction type and instance
real(pReal), dimension(:,:,:,:,:), allocatable :: constitutive_nonlocal_v, & ! dislocation velocity real(pReal), dimension(:,:,:,:,:), allocatable :: constitutive_nonlocal_v, & ! dislocation velocity
constitutive_nonlocal_rhoDotFlux ! dislocation convection term constitutive_nonlocal_rhoDotFlux ! dislocation convection term
@ -1230,17 +1230,9 @@ integer(pInt) myInstance, & ! current
neighboring_ip, & ! integration point of my neighbor neighboring_ip, & ! integration point of my neighbor
c, & ! character of dislocation c, & ! character of dislocation
n, & ! index of my current neighbor n, & ! index of my current neighbor
opposite_n, & ! index of my opposite neighbor
opposite_ip, & ! ip of my opposite neighbor
opposite_el, & ! element index of my opposite neighbor
t, & ! type of dislocation t, & ! type of dislocation
s, & ! index of my current slip system s, & ! index of my current slip system
sLattice ! index of my current slip system according to lattice order sLattice ! index of my current slip system according to lattice order
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),8) :: &
rhoSgl, & ! current single dislocation densities (positive/negative screw and edge without dipoles)
previousRhoSgl, & ! previous single dislocation densities (positive/negative screw and edge without dipoles)
totalRhoDotSgl, & ! total rate of change of single dislocation densities
thisRhoDotSgl ! rate of change of single dislocation densities for this mechanism
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: & real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: &
fluxdensity, & ! flux density at central material point fluxdensity, & ! flux density at central material point
neighboring_fluxdensity, &! flux density at neighbroing material point neighboring_fluxdensity, &! flux density at neighbroing material point
@ -1428,7 +1420,6 @@ fluxdensity = rhoSgl(:,1:4) * constitutive_nonlocal_v(:,:,g,ip,el)
do n = 1,FE_NipNeighbors(mesh_element(2,el)) ! loop through my neighbors do n = 1,FE_NipNeighbors(mesh_element(2,el)) ! loop through my neighbors
opposite_n = n - 1_pInt + 2_pInt*mod(n,2_pInt) opposite_n = n - 1_pInt + 2_pInt*mod(n,2_pInt)
neighboring_el = mesh_ipNeighborhood(1,n,ip,el) neighboring_el = mesh_ipNeighborhood(1,n,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,n,ip,el) neighboring_ip = mesh_ipNeighborhood(2,n,ip,el)

86
code/crystallite.f90
View File

@ -291,7 +291,8 @@ subroutine crystallite_init(Temperature)
write(6,'(a35,x,7(i5,x))') 'crystallite_subTstar0_v: ', shape(crystallite_subTstar0_v) write(6,'(a35,x,7(i5,x))') 'crystallite_subTstar0_v: ', shape(crystallite_subTstar0_v)
write(6,'(a35,x,7(i5,x))') 'crystallite_dPdF: ', shape(crystallite_dPdF) write(6,'(a35,x,7(i5,x))') 'crystallite_dPdF: ', shape(crystallite_dPdF)
write(6,'(a35,x,7(i5,x))') 'crystallite_fallbackdPdF: ', shape(crystallite_fallbackdPdF) write(6,'(a35,x,7(i5,x))') 'crystallite_fallbackdPdF: ', shape(crystallite_fallbackdPdF)
write(6,'(a35,x,7(i5,x))') 'crystallite_orientation: ', shape(crystallite_orientation) write(6,'(a35,x,7(i5,x))') 'crystallite_R: ', shape(crystallite_R)
write(6,'(a35,x,7(i5,x))') 'crystallite_eulerangles: ', shape(crystallite_eulerangles)
write(6,'(a35,x,7(i5,x))') 'crystallite_misorientation: ', shape(crystallite_misorientation) write(6,'(a35,x,7(i5,x))') 'crystallite_misorientation: ', shape(crystallite_misorientation)
write(6,'(a35,x,7(i5,x))') 'crystallite_dt: ', shape(crystallite_dt) write(6,'(a35,x,7(i5,x))') 'crystallite_dt: ', shape(crystallite_dt)
write(6,'(a35,x,7(i5,x))') 'crystallite_subdt: ', shape(crystallite_subdt) write(6,'(a35,x,7(i5,x))') 'crystallite_subdt: ', shape(crystallite_subdt)
@ -299,12 +300,11 @@ subroutine crystallite_init(Temperature)
write(6,'(a35,x,7(i5,x))') 'crystallite_subStep: ', shape(crystallite_subStep) write(6,'(a35,x,7(i5,x))') 'crystallite_subStep: ', shape(crystallite_subStep)
write(6,'(a35,x,7(i5,x))') 'crystallite_localConstitution: ', shape(crystallite_localConstitution) write(6,'(a35,x,7(i5,x))') 'crystallite_localConstitution: ', shape(crystallite_localConstitution)
write(6,'(a35,x,7(i5,x))') 'crystallite_requested: ', shape(crystallite_requested) write(6,'(a35,x,7(i5,x))') 'crystallite_requested: ', shape(crystallite_requested)
write(6,'(a35,x,7(i5,x))') 'crystallite_todo: ', shape(crystallite_todo) write(6,'(a35,x,7(i5,x))') 'crystallite_onTrack: ', shape(crystallite_onTrack)
write(6,'(a35,x,7(i5,x))') 'crystallite_converged: ', shape(crystallite_converged) write(6,'(a35,x,7(i5,x))') 'crystallite_converged: ', shape(crystallite_converged)
write(6,'(a35,x,7(i5,x))') 'crystallite_stateConverged: ', shape(crystallite_stateConverged) write(6,'(a35,x,7(i5,x))') 'crystallite_stateConverged: ', shape(crystallite_stateConverged)
write(6,'(a35,x,7(i5,x))') 'crystallite_temperatureConverged: ', shape(crystallite_temperatureConverged) write(6,'(a35,x,7(i5,x))') 'crystallite_temperatureConverged: ', shape(crystallite_temperatureConverged)
write(6,'(a35,x,7(i5,x))') 'crystallite_sizePostResults: ', shape(crystallite_sizePostResults) write(6,'(a35,x,7(i5,x))') 'crystallite_todo: ', shape(crystallite_todo)
write(6,'(a35,x,7(i5,x))') 'crystallite_sizePostResult: ', shape(crystallite_sizePostResult)
write(6,*) write(6,*)
write(6,*) 'Number of nonlocal grains: ',count(.not. crystallite_localConstitution) write(6,*) 'Number of nonlocal grains: ',count(.not. crystallite_localConstitution)
call flush(6) call flush(6)
@ -333,14 +333,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
pert_Fg, & pert_Fg, &
pert_method, & pert_method, &
nState, & nState, &
nCryst, & nCryst
iJacoStiffness
use debug, only: debugger, & use debug, only: debugger, &
selectiveDebugger, &
verboseDebugger, &
debug_e, &
debug_i, &
debug_g, &
debug_CrystalliteLoopDistribution, & debug_CrystalliteLoopDistribution, &
debug_CrystalliteStateLoopDistribution, & debug_CrystalliteStateLoopDistribution, &
debug_StiffnessStateLoopDistribution debug_StiffnessStateLoopDistribution
@ -351,7 +345,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
math_Mandel6to33, & math_Mandel6to33, &
math_Mandel33to6, & math_Mandel33to6, &
math_I3, & math_I3, &
math_Plain3333to99 math_Plain3333to99, &
math_EulerToR
use FEsolving, only: FEsolving_execElem, & use FEsolving, only: FEsolving_execElem, &
FEsolving_execIP, & FEsolving_execIP, &
theInc, & theInc, &
@ -385,7 +380,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
!*** local variables ***! !*** local variables ***!
real(pReal) myTemperature, & ! local copy of the temperature real(pReal) myTemperature, & ! local copy of the temperature
myPert ! perturbation with correct sign myPert, & ! perturbation with correct sign
rexCritStrain ! perturbation with correct sign
real(pReal), dimension(3,3) :: invFp, & ! inverse of the plastic deformation gradient real(pReal), dimension(3,3) :: invFp, & ! inverse of the plastic deformation gradient
Fe_guess, & ! guess for elastic deformation gradient Fe_guess, & ! guess for elastic deformation gradient
Tstar ! 2nd Piola-Kirchhoff stress tensor Tstar ! 2nd Piola-Kirchhoff stress tensor
@ -401,8 +397,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains, & myNgrains, &
mySizeState, & mySizeState, &
mySizeDotState mySizeDotState
integer(pInt), dimension(2,9,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & integer(pInt), dimension(2,9) :: kl
kl
logical onTrack, & ! flag indicating whether we are still on track logical onTrack, & ! flag indicating whether we are still on track
temperatureConverged, & ! flag indicating if temperature converged temperatureConverged, & ! flag indicating if temperature converged
stateConverged, & ! flag indicating if state converged stateConverged, & ! flag indicating if state converged
@ -412,8 +407,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
real(pReal), dimension(constitutive_maxSizeDotState) :: delta_dotState1, & ! difference between current and previous dotstate real(pReal), dimension(constitutive_maxSizeDotState) :: delta_dotState1, & ! difference between current and previous dotstate
delta_dotState2 ! difference between previousDotState and previousDotState2 delta_dotState2 ! difference between previousDotState and previousDotState2
real(pReal) dot_prod12, & real(pReal) dot_prod12, &
dot_prod22, & dot_prod22
formerSubStep
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
storedF, & storedF, &
storedFp, & storedFp, &
@ -432,9 +426,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
storedConvergenceFlag storedConvergenceFlag
logical, dimension(3,3) :: mask logical, dimension(3,3) :: mask
logical forceLocalStiffnessCalculation ! flag indicating that stiffness calculation is always done locally
forceLocalStiffnessCalculation = .false.
! ------ initialize to starting condition ------ ! ------ initialize to starting condition ------
@ -451,17 +442,17 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
crystallite_subStep = 0.0_pReal crystallite_subStep = 0.0_pReal
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
if (crystallite_requested(g,i,e)) then ! initialize restoration point of ... if (crystallite_requested(g,i,e)) then ! initialize restoration point of ...
crystallite_subTemperature0(g,i,e) = crystallite_partionedTemperature0(g,i,e) ! ...temperature crystallite_subTemperature0(g,i,e) = crystallite_partionedTemperature0(g,i,e) ! ...temperature
constitutive_subState0(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructure constitutive_subState0(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructure
crystallite_subFp0(:,:,g,i,e) = crystallite_partionedFp0(:,:,g,i,e) ! ...plastic def grad crystallite_subFp0(:,:,g,i,e) = crystallite_partionedFp0(:,:,g,i,e) ! ...plastic def grad
crystallite_subLp0(:,:,g,i,e) = crystallite_partionedLp0(:,:,g,i,e) ! ...plastic velocity grad crystallite_subLp0(:,:,g,i,e) = crystallite_partionedLp0(:,:,g,i,e) ! ...plastic velocity grad
crystallite_subF0(:,:,g,i,e) = crystallite_partionedF0(:,:,g,i,e) ! ...def grad crystallite_subF0(:,:,g,i,e) = crystallite_partionedF0(:,:,g,i,e) ! ...def grad
crystallite_subTstar0_v(:,g,i,e) = crystallite_partionedTstar0_v(:,g,i,e) !...2nd PK stress crystallite_subTstar0_v(:,g,i,e) = crystallite_partionedTstar0_v(:,g,i,e) ! ...2nd PK stress
crystallite_subFrac(g,i,e) = 0.0_pReal crystallite_subFrac(g,i,e) = 0.0_pReal
crystallite_subStep(g,i,e) = 1.0_pReal/subStepSizeCryst crystallite_subStep(g,i,e) = 1.0_pReal/subStepSizeCryst
@ -472,7 +463,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo enddo
enddo enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
! --+>> crystallite loop <<+-- ! --+>> crystallite loop <<+--
@ -553,7 +543,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
if (crystallite_todo(g,i,e)) then ! all undone crystallites if (crystallite_todo(g,i,e)) then ! all undone crystallites
call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, &
@ -583,7 +573,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
if (crystallite_todo(g,i,e)) then ! all undone crystallites if (crystallite_todo(g,i,e)) then ! all undone crystallites
@ -612,9 +602,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
! to account for substepping within _integrateStress ! to account for substepping within _integrateStress
! results in crystallite_Fp,.._Lp ! results in crystallite_Fp,.._Lp
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
if (crystallite_todo(g,i,e)) then ! all undone crystallites if (crystallite_todo(g,i,e)) then ! all undone crystallites
@ -645,10 +635,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
if (crystallite_todo(g,i,e)) then ! all undone crystallites if (crystallite_todo(g,i,e)) then ! all undone crystallites
constitutive_previousDotState2(g,i,e)%p = constitutive_previousDotState(g,i,e)%p constitutive_previousDotState2(g,i,e)%p = constitutive_previousDotState(g,i,e)%p
constitutive_previousDotState(g,i,e)%p = constitutive_dotState(g,i,e)%p constitutive_previousDotState(g,i,e)%p = constitutive_dotState(g,i,e)%p
constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotState constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotState
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
@ -658,7 +648,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
if (crystallite_todo(g,i,e)) then ! all undone crystallites if (crystallite_todo(g,i,e)) then ! all undone crystallites
@ -679,9 +669,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
if (crystallite_todo(g,i,e)) then ! all undone crystallites if (crystallite_todo(g,i,e)) then ! all undone crystallites
@ -728,18 +718,18 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
!$OMPEND CRITICAL (write2out) !$OMPEND CRITICAL (write2out)
endif endif
enddo ! crystallite convergence loop enddo ! crystallite convergence loop
NiterationCrystallite = NiterationCrystallite + 1 NiterationCrystallite = NiterationCrystallite + 1
enddo ! cutback loop enddo ! cutback loop
! ------ check for non-converged crystallites ------ ! ------ check for non-converged crystallites ------
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
if (.not. crystallite_converged(g,i,e)) then ! respond fully elastically (might be not required due to becoming terminally ill anyway) if (.not. crystallite_converged(g,i,e)) then ! respond fully elastically (might be not required due to becoming terminally ill anyway)
! call IO_warning(600,e,i,g) ! call IO_warning(600,e,i,g)
invFp = math_inv3x3(crystallite_partionedFp0(:,:,g,i,e)) invFp = math_inv3x3(crystallite_partionedFp0(:,:,g,i,e))
Fe_guess = math_mul33x33(crystallite_partionedF(:,:,g,i,e),invFp) Fe_guess = math_mul33x33(crystallite_partionedF(:,:,g,i,e),invFp)
@ -861,11 +851,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo ! perturbation direction enddo ! perturbation direction
select case(pert_method) select case(pert_method)
case (1) case (1)
crystallite_dPdF(:,:,:,:,g,i,e) = dPdF_perturbation(:,:,:,:,1) crystallite_dPdF(:,:,:,:,g,i,e) = dPdF_perturbation(:,:,:,:,1)
case (2) case (2)
crystallite_dPdF(:,:,:,:,g,i,e) = dPdF_perturbation(:,:,:,:,2) crystallite_dPdF(:,:,:,:,g,i,e) = dPdF_perturbation(:,:,:,:,2)
case (3) case (3)
crystallite_dPdF(:,:,:,:,g,i,e) = 0.5_pReal*(dPdF_perturbation(:,:,:,:,1)+dPdF_perturbation(:,:,:,:,2)) crystallite_dPdF(:,:,:,:,g,i,e) = 0.5_pReal*(dPdF_perturbation(:,:,:,:,1)+dPdF_perturbation(:,:,:,:,2))
end select end select
else ! grain did not converge else ! grain did not converge
crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! use (elastic) fallback crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! use (elastic) fallback
@ -876,7 +866,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo ! element loop enddo ! element loop
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
elseif (any(.not. crystallite_localConstitution)) then ! if any nonlocal grain present, we have to do a full loop over all grains after each perturbance elseif (any(.not. crystallite_localConstitution)) then ! if any nonlocal grain present, we have to do a full loop over all grains after each perturbance
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))

2
code/debug.f90
View File

@ -41,8 +41,6 @@ subroutine debug_init()
nHomog nHomog
implicit none implicit none
if (.not. debugger) verboseDebugger = .false.
write(6,*) write(6,*)
write(6,*) '<<<+- debug init -+>>>' write(6,*) '<<<+- debug init -+>>>'
write(6,*) '$Id$' write(6,*) '$Id$'

3
code/homogenization.f90
View File

@ -263,6 +263,7 @@ subroutine materialpoint_stressAndItsTangent(&
debug_i, & debug_i, &
debug_MaterialpointLoopDistribution, & debug_MaterialpointLoopDistribution, &
debug_MaterialpointStateLoopDistribution debug_MaterialpointStateLoopDistribution
use math, only: math_pDecomposition
implicit none implicit none
@ -270,6 +271,7 @@ subroutine materialpoint_stressAndItsTangent(&
logical, intent(in) :: updateJaco logical, intent(in) :: updateJaco
integer(pInt) NiterationHomog,NiterationMPstate integer(pInt) NiterationHomog,NiterationMPstate
integer(pInt) g,i,e,myNgrains integer(pInt) g,i,e,myNgrains
logical error
! ------ initialize to starting condition ------ ! ------ initialize to starting condition ------
@ -283,6 +285,7 @@ subroutine materialpoint_stressAndItsTangent(&
write (6,'(a,/,3(3(f12.7,x)/))') 'Lp0 of 1 1 1',crystallite_Lp0(1:3,:,1,1,1) write (6,'(a,/,3(3(f12.7,x)/))') 'Lp0 of 1 1 1',crystallite_Lp0(1:3,:,1,1,1)
endif endif
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))

304
code/homogenization_RGC.f90
View File

@ -18,14 +18,15 @@ MODULE homogenization_RGC
character (len=*), parameter :: homogenization_RGC_label = 'rgc' character (len=*), parameter :: homogenization_RGC_label = 'rgc'
integer(pInt), dimension(:), allocatable :: homogenization_RGC_sizeState, & integer(pInt), dimension(:), allocatable :: homogenization_RGC_sizeState, &
homogenization_RGC_sizePostResults homogenization_RGC_sizePostResults
integer(pInt), dimension(:,:), allocatable,target :: homogenization_RGC_sizePostResult integer(pInt), dimension(:,:), allocatable,target :: homogenization_RGC_sizePostResult
integer(pInt), dimension(:,:), allocatable :: homogenization_RGC_Ngrains integer(pInt), dimension(:,:), allocatable :: homogenization_RGC_Ngrains
real(pReal), dimension(:,:), allocatable :: homogenization_RGC_xiAlpha, & real(pReal), dimension(:,:), allocatable :: homogenization_RGC_dAlpha, &
homogenization_RGC_ciAlpha homogenization_RGC_angles
real(pReal), dimension(:), allocatable :: homogenization_RGC_maxVol0, & real(pReal), dimension(:,:,:,:), allocatable :: homogenization_RGC_orientation
homogenization_RGC_vPower0 real(pReal), dimension(:), allocatable :: homogenization_RGC_xiAlpha, &
homogenization_RGC_ciAlpha
character(len=64), dimension(:,:), allocatable,target :: homogenization_RGC_output ! name of each post result output character(len=64), dimension(:,:), allocatable,target :: homogenization_RGC_output ! name of each post result output
CONTAINS CONTAINS
@ -71,11 +72,15 @@ subroutine homogenization_RGC_init(&
allocate(homogenization_RGC_sizePostResult(maxval(homogenization_Noutput), & allocate(homogenization_RGC_sizePostResult(maxval(homogenization_Noutput), &
maxNinstance)); homogenization_RGC_sizePostResult = 0_pInt maxNinstance)); homogenization_RGC_sizePostResult = 0_pInt
allocate(homogenization_RGC_Ngrains(3,maxNinstance)); homogenization_RGC_Ngrains = 0_pInt allocate(homogenization_RGC_Ngrains(3,maxNinstance)); homogenization_RGC_Ngrains = 0_pInt
allocate(homogenization_RGC_ciAlpha(3,maxNinstance)); homogenization_RGC_ciAlpha = 0.0_pReal allocate(homogenization_RGC_ciAlpha(maxNinstance)); homogenization_RGC_ciAlpha = 0.0_pReal
allocate(homogenization_RGC_xiAlpha(3,maxNinstance)); homogenization_RGC_xiAlpha = 0.0_pReal allocate(homogenization_RGC_xiAlpha(maxNinstance)); homogenization_RGC_xiAlpha = 0.0_pReal
allocate(homogenization_RGC_maxVol0(maxNinstance)); homogenization_RGC_maxVol0 = 0.0_pReal allocate(homogenization_RGC_dAlpha(3,maxNinstance)); homogenization_RGC_dAlpha = 0.0_pReal
allocate(homogenization_RGC_vPower0(maxNinstance)); homogenization_RGC_vPower0 = 0.0_pReal allocate(homogenization_RGC_angles(3,maxNinstance)); homogenization_RGC_angles = 400.0_pReal
allocate(homogenization_RGC_output(maxval(homogenization_Noutput),maxNinstance)); homogenization_RGC_output = '' allocate(homogenization_RGC_output(maxval(homogenization_Noutput),maxNinstance)); homogenization_RGC_output = ''
allocate(homogenization_RGC_orientation(3,3,mesh_maxNips,mesh_NcpElems))
forall (i = 1:mesh_maxNips,e = 1:mesh_NcpElems)
homogenization_RGC_orientation(:,:,i,e) = math_I3
end forall
rewind(file) rewind(file)
line = '' line = ''
@ -105,23 +110,55 @@ subroutine homogenization_RGC_init(&
homogenization_RGC_Ngrains(1,i) = IO_intValue(line,positions,2) homogenization_RGC_Ngrains(1,i) = IO_intValue(line,positions,2)
homogenization_RGC_Ngrains(2,i) = IO_intValue(line,positions,3) homogenization_RGC_Ngrains(2,i) = IO_intValue(line,positions,3)
homogenization_RGC_Ngrains(3,i) = IO_intValue(line,positions,4) homogenization_RGC_Ngrains(3,i) = IO_intValue(line,positions,4)
case ('grainsizeparameter') case ('scalingparameter')
homogenization_RGC_xiAlpha(1,i) = IO_floatValue(line,positions,2) homogenization_RGC_xiAlpha(i) = IO_floatValue(line,positions,2)
homogenization_RGC_xiAlpha(2,i) = IO_floatValue(line,positions,3)
homogenization_RGC_xiAlpha(3,i) = IO_floatValue(line,positions,4)
case ('overproportionality') case ('overproportionality')
homogenization_RGC_ciAlpha(1,i) = IO_floatValue(line,positions,2) homogenization_RGC_ciAlpha(i) = IO_floatValue(line,positions,2)
homogenization_RGC_ciAlpha(2,i) = IO_floatValue(line,positions,3) case ('grainsize')
homogenization_RGC_ciAlpha(3,i) = IO_floatValue(line,positions,4) homogenization_RGC_dAlpha(1,i) = IO_floatValue(line,positions,2)
case ('maxvoldiscrepancy') homogenization_RGC_dAlpha(2,i) = IO_floatValue(line,positions,3)
homogenization_RGC_maxVol0(i) = IO_floatValue(line,positions,2) homogenization_RGC_dAlpha(3,i) = IO_floatValue(line,positions,4)
case ('discrepancypower') case ('clusterorientation')
homogenization_RGC_vPower0(i) = IO_floatValue(line,positions,2) homogenization_RGC_angles(1,i) = IO_floatValue(line,positions,2)
homogenization_RGC_angles(2,i) = IO_floatValue(line,positions,3)
homogenization_RGC_angles(3,i) = IO_floatValue(line,positions,4)
end select end select
endif endif
enddo enddo
100 do i = 1,maxNinstance ! sanity checks !*** assigning cluster orientations
do e = 1,mesh_NcpElems
if (homogenization_type(mesh_element(3,e)) == homogenization_RGC_label) then
myInstance = homogenization_typeInstance(mesh_element(3,e))
if (all (homogenization_RGC_angles(:,myInstance) >= 399.9_pReal)) then
homogenization_RGC_orientation(:,:,1,e) = math_EulerToR(math_sampleRandomOri())
do i = 1,FE_Nips(mesh_element(2,e))
if (microstructure_elemhomo(mesh_element(4,e))) then
homogenization_RGC_orientation(:,:,i,e) = homogenization_RGC_orientation(:,:,1,e)
else
homogenization_RGC_orientation(:,:,i,e) = math_EulerToR(math_sampleRandomOri())
endif
enddo
else
do i = 1,FE_Nips(mesh_element(2,e))
homogenization_RGC_orientation(:,:,i,e) = math_EulerToR(homogenization_RGC_angles(:,myInstance)*inRad)
enddo
endif
endif
enddo
100 do i = 1,maxNinstance ! sanity checks
write(6,*)
write(6,*) '<<<+- homogenization_RGC init -+>>>'
write(6,*) '$Id$'
write(6,*)
write(6,'(a15,x,i4)') 'instance: ', i
write(6,*)
write(6,'(a25,3(x,i8))') 'cluster size: ',(homogenization_RGC_Ngrains(j,i),j=1,3)
write(6,'(a25,x,e10.3)') 'scaling parameter: ', homogenization_RGC_xiAlpha(i)
write(6,'(a25,x,e10.3)') 'over-proportionality: ', homogenization_RGC_ciAlpha(i)
write(6,'(a25,3(x,e10.3))') 'grain size: ',(homogenization_RGC_dAlpha(j,i),j=1,3)
write(6,'(a25,3(x,e10.3))') 'cluster orientation: ',(homogenization_RGC_angles(j,i),j=1,3)
enddo enddo
do i = 1,maxNinstance do i = 1,maxNinstance
@ -130,13 +167,17 @@ subroutine homogenization_RGC_init(&
case('constitutivework') case('constitutivework')
mySize = 1 mySize = 1
case('magnitudemismatch') case('magnitudemismatch')
mySize = 1 mySize = 3
case('penaltyenergy') case('penaltyenergy')
mySize = 1 mySize = 1
case('volumediscrepancy') case('volumediscrepancy')
mySize = 1 mySize = 1
case default case default
mySize = 0 mySize = 0
case('averagerelaxrate')
mySize = 1
case('maximumrelaxrate')
mySize = 1
end select end select
if (mySize > 0_pInt) then ! any meaningful output found if (mySize > 0_pInt) then ! any meaningful output found
@ -146,13 +187,12 @@ subroutine homogenization_RGC_init(&
endif endif
enddo enddo
homogenization_RGC_sizeState(i) & homogenization_RGC_sizeState(i) &
= 3*(homogenization_RGC_Ngrains(1,i)-1)*homogenization_RGC_Ngrains(2,i)*homogenization_RGC_Ngrains(3,i) & = 3*(homogenization_RGC_Ngrains(1,i)-1)*homogenization_RGC_Ngrains(2,i)*homogenization_RGC_Ngrains(3,i) &
+ 3*homogenization_RGC_Ngrains(1,i)*(homogenization_RGC_Ngrains(2,i)-1)*homogenization_RGC_Ngrains(3,i) & + 3*homogenization_RGC_Ngrains(1,i)*(homogenization_RGC_Ngrains(2,i)-1)*homogenization_RGC_Ngrains(3,i) &
+ 3*homogenization_RGC_Ngrains(1,i)*homogenization_RGC_Ngrains(2,i)*(homogenization_RGC_Ngrains(3,i)-1) & + 3*homogenization_RGC_Ngrains(1,i)*homogenization_RGC_Ngrains(2,i)*(homogenization_RGC_Ngrains(3,i)-1) &
+ 4_pInt ! (1) Average constitutive work, (2) Overall mismatch, (3) Average penalty energy, + 8_pInt ! (1) Average constitutive work, (2-4) Overall mismatch, (5) Average penalty energy,
! (4) Volume discrepancy ! (6) Volume discrepancy, (7) Avg relaxation rate component, (8) Max relaxation rate component
enddo enddo
return return
@ -235,7 +275,7 @@ subroutine homogenization_RGC_partitionDeformation(&
do iFace = 1,nFace do iFace = 1,nFace
call homogenization_RGC_getInterface(intFace,iFace,iGrain3) call homogenization_RGC_getInterface(intFace,iFace,iGrain3)
call homogenization_RGC_relaxationVector(aVect,intFace,state,homID) call homogenization_RGC_relaxationVector(aVect,intFace,state,homID)
call homogenization_RGC_interfaceNormal(nVect,intFace) call homogenization_RGC_interfaceNormal(nVect,intFace,ip,el)
forall (i=1:3,j=1:3) & forall (i=1:3,j=1:3) &
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
enddo enddo
@ -279,7 +319,7 @@ function homogenization_RGC_updateState(&
use material, only: homogenization_maxNgrains,homogenization_typeInstance, & use material, only: homogenization_maxNgrains,homogenization_typeInstance, &
homogenization_Ngrains homogenization_Ngrains
use numerics, only: absTol_RGC,relTol_RGC,absMax_RGC,relMax_RGC,pPert_RGC, & use numerics, only: absTol_RGC,relTol_RGC,absMax_RGC,relMax_RGC,pPert_RGC, &
maxdRelax_RGC,ratePower_RGC,viscModus_RGC maxdRelax_RGC,viscPower_RGC,viscModus_RGC,refRelaxRate_RGC
use FEsolving, only: theInc,cycleCounter,theTime use FEsolving, only: theInc,cycleCounter,theTime
implicit none implicit none
@ -299,7 +339,7 @@ function homogenization_RGC_updateState(&
integer(pInt), dimension (2) :: residLoc integer(pInt), dimension (2) :: residLoc
integer(pInt) homID,i1,i2,i3,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ival,ipert,iGrain,nGrain integer(pInt) homID,i1,i2,i3,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ival,ipert,iGrain,nGrain
real(pReal), dimension (3,3,homogenization_maxNgrains) :: R,pF,pR,D,pD real(pReal), dimension (3,3,homogenization_maxNgrains) :: R,pF,pR,D,pD
real(pReal), dimension (homogenization_maxNgrains) :: NN,pNN real(pReal), dimension (3,homogenization_maxNgrains) :: NN,pNN
real(pReal), dimension (3) :: normP,normN,mornP,mornN real(pReal), dimension (3) :: normP,normN,mornP,mornN
real(pReal) residMax,stresMax,constitutiveWork,penaltyEnergy,volDiscrep,penDiscrep real(pReal) residMax,stresMax,constitutiveWork,penaltyEnergy,volDiscrep,penDiscrep
logical error,RGCdebug,RGCdebugJacobi,RGCcheck logical error,RGCdebug,RGCdebugJacobi,RGCcheck
@ -336,7 +376,7 @@ function homogenization_RGC_updateState(&
endif endif
!* Stress-like penalty related to mismatch or incompatibility at interfaces !* Stress-like penalty related to mismatch or incompatibility at interfaces
call homogenization_RGC_stressPenalty(R,NN,F,ip,el,homID) call homogenization_RGC_stressPenalty(R,NN,avgF,F,ip,el,homID)
!* Stress-like penalty related to overall volume discrepancy !* Stress-like penalty related to overall volume discrepancy
call homogenization_RGC_volumePenalty(D,volDiscrep,F,avgF,ip,el,homID) call homogenization_RGC_volumePenalty(D,volDiscrep,F,avgF,ip,el,homID)
@ -344,7 +384,7 @@ function homogenization_RGC_updateState(&
!* Debugging the mismatch, stress and penalties of grains !* Debugging the mismatch, stress and penalties of grains
if (RGCdebug) then if (RGCdebug) then
do iGrain = 1,nGrain do iGrain = 1,nGrain
write(6,'(x,a30,x,i3,x,a4,x,e14.8)')'Mismatch magnitude of grain(',iGrain,') :',NN(iGrain) write(6,'(x,a30,x,i3,x,a4,3(x,e14.8))')'Mismatch magnitude of grain(',iGrain,') :',NN(1,iGrain),NN(2,iGrain),NN(3,iGrain)
write(6,*)' ' write(6,*)' '
write(6,'(x,a30,x,i3)')'Stress and penalties of grain: ',iGrain write(6,'(x,a30,x,i3)')'Stress and penalties of grain: ',iGrain
do i = 1,3 do i = 1,3
@ -364,15 +404,15 @@ function homogenization_RGC_updateState(&
iGr3N = faceID(2:4) iGr3N = faceID(2:4)
call homogenization_RGC_grain3to1(iGrN,iGr3N,homID) call homogenization_RGC_grain3to1(iGrN,iGr3N,homID)
call homogenization_RGC_getInterface(intFaceN,2*faceID(1),iGr3N) call homogenization_RGC_getInterface(intFaceN,2*faceID(1),iGr3N)
call homogenization_RGC_interfaceNormal(normN,intFaceN) ! get the interface normal call homogenization_RGC_interfaceNormal(normN,intFaceN,ip,el) ! get the interface normal
!* Identify the right/up/front grain (+|P) !* Identify the right/up/front grain (+|P)
iGr3P = iGr3N iGr3P = iGr3N
iGr3P(faceID(1)) = iGr3N(faceID(1))+1 iGr3P(faceID(1)) = iGr3N(faceID(1))+1
call homogenization_RGC_grain3to1(iGrP,iGr3P,homID) call homogenization_RGC_grain3to1(iGrP,iGr3P,homID)
call homogenization_RGC_getInterface(intFaceP,2*faceID(1)-1,iGr3P) call homogenization_RGC_getInterface(intFaceP,2*faceID(1)-1,iGr3P)
call homogenization_RGC_interfaceNormal(normP,intFaceP) ! get the interface normal call homogenization_RGC_interfaceNormal(normP,intFaceP,ip,el) ! get the interface normal
do i = 1,3 ! compute the traction balance at the interface do i = 1,3 ! compute the traction balance at the interface
tract(iNum,i) = sign(viscModus_RGC*(abs(drelax(i+3*(iNum-1))/dt))**ratePower_RGC, & tract(iNum,i) = sign(viscModus_RGC*(abs(drelax(i+3*(iNum-1)))/(refRelaxRate_RGC*dt))**viscPower_RGC, &
drelax(i+3*(iNum-1))) drelax(i+3*(iNum-1)))
do j = 1,3 do j = 1,3
tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP) + D(i,j,iGrP))*normP(j) & tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP) + D(i,j,iGrP))*normP(j) &
@ -416,7 +456,7 @@ function homogenization_RGC_updateState(&
!* Then compute/update the state for postResult, i.e., ... !* Then compute/update the state for postResult, i.e., ...
!* ... all energy densities computed by time-integration !* ... all energy densities computed by time-integration
constitutiveWork = state%p(3*nIntFaceTot+1) constitutiveWork = state%p(3*nIntFaceTot+1)
penaltyEnergy = state%p(3*nIntFaceTot+3) penaltyEnergy = state%p(3*nIntFaceTot+5)
do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) do iGrain = 1,homogenization_Ngrains(mesh_element(3,el))
do i = 1,3 do i = 1,3
do j = 1,3 do j = 1,3
@ -428,16 +468,23 @@ function homogenization_RGC_updateState(&
!* ... the bulk mechanical/constitutive work !* ... the bulk mechanical/constitutive work
state%p(3*nIntFaceTot+1) = constitutiveWork state%p(3*nIntFaceTot+1) = constitutiveWork
!* ... the overall mismatch !* ... the overall mismatch
state%p(3*nIntFaceTot+2) = sum(NN)/dble(nGrain) state%p(3*nIntFaceTot+2) = sum(NN(1,:))/dble(nGrain)
state%p(3*nIntFaceTot+3) = penaltyEnergy state%p(3*nIntFaceTot+3) = sum(NN(2,:))/dble(nGrain)
state%p(3*nIntFaceTot+4) = sum(NN(3,:))/dble(nGrain)
state%p(3*nIntFaceTot+5) = penaltyEnergy
!* ... the volume discrepancy !* ... the volume discrepancy
state%p(3*nIntFaceTot+4) = volDiscrep state%p(3*nIntFaceTot+6) = volDiscrep
state%p(3*nIntFaceTot+7) = sum(abs(drelax))/dt/dble(3*nIntFaceTot)
state%p(3*nIntFaceTot+8) = maxval(abs(drelax))/dt
if (el == 1 .and. ip == 1) then if (el == 1 .and. ip == 1) then
write(6,'(x,a30,x,e14.8)')'Constitutive work: ',constitutiveWork write(6,'(x,a30,x,e14.8)')'Constitutive work: ',constitutiveWork
write(6,'(x,a30,x,e14.8)')'Magnitude mismatch: ',sum(NN) write(6,'(x,a30,3(x,e14.8))')'Magnitude mismatch: ',sum(NN(1,:))/dble(nGrain),sum(NN(2,:))/dble(nGrain),sum(NN(3,:))/dble(nGrain)
write(6,'(x,a30,x,e14.8)')'Penalty energy: ',penaltyEnergy write(6,'(x,a30,x,e14.8)')'Penalty energy: ',penaltyEnergy
write(6,'(x,a30,x,e14.8)')'Volume discrepancy: ',volDiscrep write(6,'(x,a30,x,e14.8)')'Volume discrepancy: ',volDiscrep
write(6,*)'' write(6,*)''
write(6,'(x,a30,x,e14.8)')'Maximum relaxation rate: ',maxval(abs(drelax))/dt
write(6,'(x,a30,x,e14.8)')'Average relaxation rate: ',sum(abs(drelax))/dt/dble(3*nIntFaceTot)
write(6,*)''
call flush(6) call flush(6)
endif endif
deallocate(tract,resid,relax,drelax) deallocate(tract,resid,relax,drelax)
@ -446,8 +493,8 @@ function homogenization_RGC_updateState(&
elseif (residMax > relMax_RGC*stresMax .or. residMax > absMax_RGC) then elseif (residMax > relMax_RGC*stresMax .or. residMax > absMax_RGC) then
!* Try to restart when residual blows up exceeding maximum bound !* Try to restart when residual blows up exceeding maximum bound
homogenization_RGC_updateState = (/.true.,.false./) ! ... with direct cut-back homogenization_RGC_updateState = (/.true.,.false./) ! ... with direct cut-back
write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback' ! write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback'
write(6,'(x,a,x,e14.8,x,a,x,e14.8)')'due to high residual =',residMax,'for stress =',stresMax ! write(6,'(x,a,x,e14.8,x,a,x,e14.8)')'due to high residual =',residMax,'for stress =',stresMax
! state%p(1:3*nIntFaceTot) = 0.0_pReal ! ... with full Taylor assumption ! state%p(1:3*nIntFaceTot) = 0.0_pReal ! ... with full Taylor assumption
! write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'relaxation vectors reset to zero' ! write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'relaxation vectors reset to zero'
if (RGCcheck) then if (RGCcheck) then
@ -475,10 +522,10 @@ function homogenization_RGC_updateState(&
iGr3N = faceID(2:4) iGr3N = faceID(2:4)
call homogenization_RGC_grain3to1(iGrN,iGr3N,homID) call homogenization_RGC_grain3to1(iGrN,iGr3N,homID)
call homogenization_RGC_getInterface(intFaceN,2*faceID(1),iGr3N) call homogenization_RGC_getInterface(intFaceN,2*faceID(1),iGr3N)
call homogenization_RGC_interfaceNormal(normN,intFaceN) ! get the interface normal call homogenization_RGC_interfaceNormal(normN,intFaceN,ip,el) ! get the interface normal
do iFace = 1,nFace do iFace = 1,nFace
call homogenization_RGC_getInterface(intFaceN,iFace,iGr3N) call homogenization_RGC_getInterface(intFaceN,iFace,iGr3N)
call homogenization_RGC_interfaceNormal(mornN,intFaceN) ! get influencing interfaces normal call homogenization_RGC_interfaceNormal(mornN,intFaceN,ip,el) ! get influencing interfaces normal
call homogenization_RGC_interface4to1(iMun,intFaceN,homID) call homogenization_RGC_interface4to1(iMun,intFaceN,homID)
if (iMun .gt. 0) then ! get the tangent if (iMun .gt. 0) then ! get the tangent
forall(i=1:3,j=1:3,k=1:3,l=1:3) & forall(i=1:3,j=1:3,k=1:3,l=1:3) &
@ -490,10 +537,10 @@ function homogenization_RGC_updateState(&
iGr3P(faceID(1)) = iGr3N(faceID(1))+1 iGr3P(faceID(1)) = iGr3N(faceID(1))+1
call homogenization_RGC_grain3to1(iGrP,iGr3P,homID) call homogenization_RGC_grain3to1(iGrP,iGr3P,homID)
call homogenization_RGC_getInterface(intFaceP,2*faceID(1)-1,iGr3P) call homogenization_RGC_getInterface(intFaceP,2*faceID(1)-1,iGr3P)
call homogenization_RGC_interfaceNormal(normP,intFaceP) ! get the interface normal call homogenization_RGC_interfaceNormal(normP,intFaceP,ip,el) ! get the interface normal
do iFace = 1,nFace do iFace = 1,nFace
call homogenization_RGC_getInterface(intFaceP,iFace,iGr3P) call homogenization_RGC_getInterface(intFaceP,iFace,iGr3P)
call homogenization_RGC_interfaceNormal(mornP,intFaceP) ! get influencing interfaces normal call homogenization_RGC_interfaceNormal(mornP,intFaceP,ip,el) ! get influencing interfaces normal
call homogenization_RGC_interface4to1(iMun,intFaceP,homID) call homogenization_RGC_interface4to1(iMun,intFaceP,homID)
if (iMun .gt. 0) then ! get the tangent if (iMun .gt. 0) then ! get the tangent
forall(i=1:3,j=1:3,k=1:3,l=1:3) & forall(i=1:3,j=1:3,k=1:3,l=1:3) &
@ -519,8 +566,8 @@ function homogenization_RGC_updateState(&
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
state%p(1:3*nIntFaceTot) = p_relax state%p(1:3*nIntFaceTot) = p_relax
call homogenization_RGC_grainDeformation(pF,F0,avgF,state,el) call homogenization_RGC_grainDeformation(pF,F0,avgF,state,ip,el)
call homogenization_RGC_stressPenalty(pR,pNN,pF,ip,el,homID) call homogenization_RGC_stressPenalty(pR,pNN,avgF,pF,ip,el,homID)
call homogenization_RGC_volumePenalty(pD,volDiscrep,pF,avgF,ip,el,homID) call homogenization_RGC_volumePenalty(pD,volDiscrep,pF,avgF,ip,el,homID)
p_resid = 0.0_pReal p_resid = 0.0_pReal
do iNum = 1,nIntFaceTot do iNum = 1,nIntFaceTot
@ -529,13 +576,13 @@ function homogenization_RGC_updateState(&
iGr3N = faceID(2:4) iGr3N = faceID(2:4)
call homogenization_RGC_grain3to1(iGrN,iGr3N,homID) call homogenization_RGC_grain3to1(iGrN,iGr3N,homID)
call homogenization_RGC_getInterface(intFaceN,2*faceID(1),iGr3N) call homogenization_RGC_getInterface(intFaceN,2*faceID(1),iGr3N)
call homogenization_RGC_interfaceNormal(normN,intFaceN) ! get the corresponding normal call homogenization_RGC_interfaceNormal(normN,intFaceN,ip,el) ! get the corresponding normal
!* Identify the right/up/front grain (+|P) !* Identify the right/up/front grain (+|P)
iGr3P = iGr3N iGr3P = iGr3N
iGr3P(faceID(1)) = iGr3N(faceID(1))+1 iGr3P(faceID(1)) = iGr3N(faceID(1))+1
call homogenization_RGC_grain3to1(iGrP,iGr3P,homID) call homogenization_RGC_grain3to1(iGrP,iGr3P,homID)
call homogenization_RGC_getInterface(intFaceP,2*faceID(1)-1,iGr3P) call homogenization_RGC_getInterface(intFaceP,2*faceID(1)-1,iGr3P)
call homogenization_RGC_interfaceNormal(normP,intFaceP) ! get the corresponding normal call homogenization_RGC_interfaceNormal(normP,intFaceP,ip,el) ! get the corresponding normal
!* Compute the perturbed traction at interface !* Compute the perturbed traction at interface
do i = 1,3 do i = 1,3
do j = 1,3 do j = 1,3
@ -561,7 +608,8 @@ function homogenization_RGC_updateState(&
!* Construct the Jacobian matrix of the numerical viscosity tangent !* Construct the Jacobian matrix of the numerical viscosity tangent
allocate(rmatrix(3*nIntFaceTot,3*nIntFaceTot)); rmatrix = 0.0_pReal allocate(rmatrix(3*nIntFaceTot,3*nIntFaceTot)); rmatrix = 0.0_pReal
forall (i=1:3*nIntFaceTot) & forall (i=1:3*nIntFaceTot) &
rmatrix(i,i) = viscModus_RGC*ratePower_RGC/dt*(abs(drelax(i)/dt))**(ratePower_RGC - 1.0_pReal) rmatrix(i,i) = viscModus_RGC*viscPower_RGC/(refRelaxRate_RGC*dt)* &
(abs(drelax(i))/(refRelaxRate_RGC*dt))**(viscPower_RGC - 1.0_pReal)
!* Debugging the global Jacobian matrix of numerical viscosity tangent !* Debugging the global Jacobian matrix of numerical viscosity tangent
if (RGCdebugJacobi) then if (RGCdebugJacobi) then
write(6,'(x,a30)')'Jacobian matrix of penalty' write(6,'(x,a30)')'Jacobian matrix of penalty'
@ -609,9 +657,9 @@ function homogenization_RGC_updateState(&
! state%p(1:3*nIntFaceTot) = 0.0_pReal ! state%p(1:3*nIntFaceTot) = 0.0_pReal
! write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'relaxation vectors reset to zero' ! write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'relaxation vectors reset to zero'
homogenization_RGC_updateState = (/.true.,.false./) homogenization_RGC_updateState = (/.true.,.false./)
write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback' ! write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback'
write(6,'(x,a,x,e14.8)')'due to large relaxation change =',maxval(abs(drelax)) ! write(6,'(x,a,x,e14.8)')'due to large relaxation change =',maxval(abs(drelax))
call flush(6) ! call flush(6)
endif endif
!* Debugging the return state !* Debugging the return state
if (RGCdebugJacobi) then if (RGCdebugJacobi) then
@ -744,13 +792,21 @@ pure function homogenization_RGC_postResults(&
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+1) homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+1)
c = c + 1 c = c + 1
case('magnitudemismatch') case('magnitudemismatch')
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+3)
c = c + 1
case('penaltyenergy')
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+2) 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)
c = c + 3
case('penaltyenergy')
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+5)
c = c + 1 c = c + 1
case('volumediscrepancy') case('volumediscrepancy')
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+4) homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+6)
c = c + 1
case('averagerelaxrate')
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+7)
c = c + 1
case('maximumrelaxrate')
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+8)
c = c + 1 c = c + 1
end select end select
enddo enddo
@ -766,6 +822,7 @@ subroutine homogenization_RGC_stressPenalty(&
rPen, & ! stress-like penalty rPen, & ! stress-like penalty
nMis, & ! total amount of mismatch nMis, & ! total amount of mismatch
! !
avgF, & ! initial effective stretch tensor
fDef, & ! deformation gradients fDef, & ! deformation gradients
ip, & ! integration point ip, & ! integration point
el, & ! element el, & ! element
@ -774,22 +831,23 @@ subroutine homogenization_RGC_stressPenalty(&
use prec, only: pReal,pInt,p_vec use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element use mesh, only: mesh_element
use constitutive, only: constitutive_homogenizedC use constitutive, only: constitutive_homogenizedC
use math, only: math_civita use math, only: math_civita,math_invert3x3
use material, only: homogenization_maxNgrains,homogenization_Ngrains use material, only: homogenization_maxNgrains,homogenization_Ngrains
use numerics, only: xSmoo_RGC use numerics, only: xSmoo_RGC
implicit none implicit none
!* Definition of variables !* Definition of variables
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: rPen real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: rPen
real(pReal), dimension (homogenization_maxNgrains), intent(out) :: nMis real(pReal), dimension (3,homogenization_maxNgrains), intent(out) :: nMis
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: fDef real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: fDef
real(pReal), dimension (3,3), intent(in) :: avgF
integer(pInt), intent(in) :: ip,el integer(pInt), intent(in) :: ip,el
integer(pInt), dimension (4) :: intFace integer(pInt), dimension (4) :: intFace
integer(pInt), dimension (3) :: iGrain3,iGNghb3,nGDim integer(pInt), dimension (3) :: iGrain3,iGNghb3,nGDim
real(pReal), dimension (3,3) :: gDef,nDef real(pReal), dimension (3,3) :: gDef,nDef,avgC
real(pReal), dimension (3) :: nVect real(pReal), dimension (3) :: nVect,surfCorr
integer(pInt) homID,iGrain,iGNghb,iFace,i,j,k,l integer(pInt) homID,iGrain,iGNghb,iFace,i,j,k,l,m
real(pReal) muGrain,muGNghb,nDefNorm real(pReal) muGrain,muGNghb,nDefNorm,xiAlpha,ciAlpha,bgGrain,bgGNghb,detF
! !
integer(pInt), parameter :: nFace = 6 integer(pInt), parameter :: nFace = 6
real(pReal), parameter :: nDefToler = 1.0e-10 real(pReal), parameter :: nDefToler = 1.0e-10
@ -798,13 +856,23 @@ subroutine homogenization_RGC_stressPenalty(&
rPen = 0.0_pReal rPen = 0.0_pReal
nMis = 0.0_pReal nMis = 0.0_pReal
!* Get the penalty correction factor representing surface evolution
call homogenization_RGC_surfaceCorrection(surfCorr,avgF,ip,el)
!* Debugging the surface correction factor
! if (ip == 1 .and. el == 1) then
! write(6,'(x,a20,2(x,i3))')'Correction factor: ',ip,el
! write(6,'(x,3(e10.4,x))')(surfCorr(i), i = 1,3)
! endif
do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) do iGrain = 1,homogenization_Ngrains(mesh_element(3,el))
call homogenization_RGC_equivalentModuli(muGrain,constitutive_homogenizedC(iGrain,ip,el)) call homogenization_RGC_equivalentModuli(muGrain,bgGrain,iGrain,ip,el)
call homogenization_RGC_grain1to3(iGrain3,iGrain,homID) call homogenization_RGC_grain1to3(iGrain3,iGrain,homID)
!* Compute the mismatch tensor at all six interfaces !* Compute the mismatch tensor at all six interfaces
do iFace = 1,nFace do iFace = 1,nFace
call homogenization_RGC_getInterface(intFace,iFace,iGrain3) call homogenization_RGC_getInterface(intFace,iFace,iGrain3)
call homogenization_RGC_interfaceNormal(nVect,intFace) ! get the interface normal call homogenization_RGC_interfaceNormal(nVect,intFace,ip,el) ! get the interface normal
iGNghb3 = iGrain3 ! identify the grain neighbor iGNghb3 = iGrain3 ! identify the grain neighbor
iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) + int(dble(intFace(1))/dble(abs(intFace(1)))) iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) + int(dble(intFace(1))/dble(abs(intFace(1))))
!* The grain periodicity along e1 !* The grain periodicity along e1
@ -817,7 +885,7 @@ subroutine homogenization_RGC_stressPenalty(&
if (iGNghb3(3) < 1) iGNghb3(3) = nGDim(3) if (iGNghb3(3) < 1) iGNghb3(3) = nGDim(3)
if (iGNghb3(3) > nGDim(3)) iGNghb3(3) = 1 if (iGNghb3(3) > nGDim(3)) iGNghb3(3) = 1
call homogenization_RGC_grain3to1(iGNghb,iGNghb3,homID) ! get the grain neighbor call homogenization_RGC_grain3to1(iGNghb,iGNghb3,homID) ! get the grain neighbor
call homogenization_RGC_equivalentModuli(muGNghb,constitutive_homogenizedC(iGNghb,ip,el)) call homogenization_RGC_equivalentModuli(muGNghb,bgGNghb,iGNghb,ip,el)
gDef = 0.5_pReal*(fDef(:,:,iGNghb) - fDef(:,:,iGrain)) ! difference in F with the neighbor gDef = 0.5_pReal*(fDef(:,:,iGNghb) - fDef(:,:,iGrain)) ! difference in F with the neighbor
nDefNorm = 0.0_pReal nDefNorm = 0.0_pReal
nDef = 0.0_pReal nDef = 0.0_pReal
@ -845,8 +913,9 @@ subroutine homogenization_RGC_stressPenalty(&
do j = 1,3 do j = 1,3
do k = 1,3 do k = 1,3
do l = 1,3 do l = 1,3
rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain + muGNghb)/homogenization_RGC_xiAlpha(abs(intFace(1)),homID) & rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*homogenization_RGC_xiAlpha(homID) &
*cosh(homogenization_RGC_ciAlpha(abs(intFace(1)),homID)*nDefNorm) & *surfCorr(abs(intFace(1)))/homogenization_RGC_dAlpha(abs(intFace(1)),homID) &
*cosh(homogenization_RGC_ciAlpha(homID)*nDefNorm) &
*0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_civita(k,l,j) & *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_civita(k,l,j) &
*tanh(nDefNorm/xSmoo_RGC) *tanh(nDefNorm/xSmoo_RGC)
enddo enddo
@ -854,7 +923,7 @@ subroutine homogenization_RGC_stressPenalty(&
enddo enddo
enddo enddo
!* Total amount of mismatch experienced by the grain (at all six interfaces) !* Total amount of mismatch experienced by the grain (at all six interfaces)
nMis(iGrain) = nMis(iGrain) + nDefNorm nMis(abs(intFace(1)),iGrain) = nMis(abs(intFace(1)),iGrain) + nDefNorm
enddo enddo
!* Debugging the stress-like penalty !* Debugging the stress-like penalty
! if (ip == 1 .and. el == 1) then ! if (ip == 1 .and. el == 1) then
@ -884,7 +953,6 @@ subroutine homogenization_RGC_volumePenalty(&
) )
use prec, only: pReal,pInt,p_vec use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element use mesh, only: mesh_element
use constitutive, only: constitutive_homogenizedC
use math, only: math_det3x3,math_inv3x3 use math, only: math_det3x3,math_inv3x3
use material, only: homogenization_maxNgrains,homogenization_Ngrains use material, only: homogenization_maxNgrains,homogenization_Ngrains
use numerics, only: maxVolDiscr_RGC,volDiscrMod_RGC,volDiscrPow_RGC use numerics, only: maxVolDiscr_RGC,volDiscrMod_RGC,volDiscrPow_RGC
@ -903,18 +971,18 @@ subroutine homogenization_RGC_volumePenalty(&
nGrain = homogenization_Ngrains(mesh_element(3,el)) nGrain = homogenization_Ngrains(mesh_element(3,el))
!* Compute volumes of grain and the effective volume and the total volume discrepancy !* Compute volumes of grain and the effective volume and the total volume discrepancy
vDiscrep = math_det3x3(fAvg) vDiscrep = math_det3x3(fAvg)
do iGrain = 1,nGrain do iGrain = 1,nGrain
gVol(iGrain) = math_det3x3(fDef(:,:,iGrain)) gVol(iGrain) = math_det3x3(fDef(:,:,iGrain))
vDiscrep = vDiscrep - gVol(iGrain)/dble(nGrain) vDiscrep = vDiscrep - gVol(iGrain)/dble(nGrain)
enddo enddo
!* Calculate the stress and penalty due to volume discrepancy !* Calculate the stress and penalty due to volume discrepancy
vPen = 0.0_pReal vPen = 0.0_pReal
do iGrain = 1,nGrain do iGrain = 1,nGrain
vPen(:,:,iGrain) = -1.0_pReal*sign(volDiscrMod_RGC*volDiscrPow_RGC/maxVolDiscr_RGC* & vPen(:,:,iGrain) = -1.0_pReal/dble(nGrain)*volDiscrMod_RGC*volDiscrPow_RGC/maxVolDiscr_RGC* &
(abs(vDiscrep)/maxVolDiscr_RGC)**(volDiscrPow_RGC - 1.0)*gVol(iGrain)/dble(nGrain),vDiscrep)* & sign((abs(vDiscrep)/maxVolDiscr_RGC)**(volDiscrPow_RGC - 1.0),vDiscrep)* &
transpose(math_inv3x3(fDef(:,:,iGrain))) gVol(iGrain)*transpose(math_inv3x3(fDef(:,:,iGrain)))
!* Debugging the stress-like penalty of volume discrepancy !* Debugging the stress-like penalty of volume discrepancy
! if (ip == 1 .and. el == 1) then ! if (ip == 1 .and. el == 1) then
@ -929,25 +997,79 @@ subroutine homogenization_RGC_volumePenalty(&
endsubroutine endsubroutine
!********************************************************************
! subroutine to compute the correction factor due to surface area evolution
!********************************************************************
subroutine homogenization_RGC_surfaceCorrection(&
vSurf, & ! surface correction factor
!
avgF, & ! average deformation gradient
ip, & ! my IP
el & ! my element
)
use prec, only: pReal,pInt,p_vec
use math, only: math_invert3x3,math_mul33x33
implicit none
!* Definition of variables
real(pReal), dimension(3,3), intent(in) :: avgF
real(pReal), dimension(3), intent(out) :: vSurf
integer(pInt), intent(in) :: ip,el
real(pReal), dimension(3,3) :: invC,avgC
real(pReal), dimension(3) :: nVect
real(pReal) detF
integer(pInt), dimension(4) :: intFace
integer(pInt) i,j,iBase
logical error
!* Compute the correction factor accouted for surface evolution (area change)
avgC = 0.0_pReal
avgC = math_mul33x33(transpose(avgF),avgF)
invC = 0.0_pReal
call math_invert3x3(avgC,invC,detF,error)
vSurf = 0.0_pReal
do iBase = 1,3
intFace = (/iBase,1_pInt,1_pInt,1_pInt/)
call homogenization_RGC_interfaceNormal(nVect,intFace,ip,el)
do i = 1,3
do j = 1,3
vSurf(iBase) = vSurf(iBase) + invC(i,j)*nVect(i)*nVect(j)
enddo
enddo
vSurf(iBase) = sqrt(vSurf(iBase))*detF
enddo
return
endsubroutine
!******************************************************************** !********************************************************************
! subroutine to compute the equivalent shear and bulk moduli from the elasticity tensor ! subroutine to compute the equivalent shear and bulk moduli from the elasticity tensor
!******************************************************************** !********************************************************************
subroutine homogenization_RGC_equivalentModuli(& subroutine homogenization_RGC_equivalentModuli(&
shearMod, & ! equivalent (isotropic) shear modulus shearMod, & ! equivalent (isotropic) shear modulus
vBurgers, & ! length of burgers vector
! !
elasTens & ! elasticity tensor in Mandel notation grainID, & ! grain ID
ip, & ! IP number
el & ! element number
) )
use prec, only: pReal,pInt,p_vec use prec, only: pReal,pInt,p_vec
use constitutive, only: constitutive_homogenizedC,constitutive_averageBurgers
use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips
use material, only: homogenization_typeInstance use material, only: homogenization_typeInstance
implicit none implicit none
!* Definition of variables !* Definition of variables
real(pReal), dimension (6,6), intent(in) :: elasTens integer(pInt), intent(in) :: grainID,ip,el
real(pReal), intent(out) :: shearMod real(pReal), dimension (6,6) :: elasTens
real(pReal), intent(out) :: shearMod,vBurgers
real(pReal) cEquiv_11,cEquiv_12,cEquiv_44 real(pReal) cEquiv_11,cEquiv_12,cEquiv_44
elasTens = constitutive_homogenizedC(grainID,ip,el)
!* Compute the equivalent shear modulus using Turterltaub and Suiker, JMPS (2005) !* Compute the equivalent shear modulus using Turterltaub and Suiker, JMPS (2005)
cEquiv_11 = (elasTens(1,1) + elasTens(2,2) + elasTens(3,3))/3.0_pReal cEquiv_11 = (elasTens(1,1) + elasTens(2,2) + elasTens(3,3))/3.0_pReal
cEquiv_12 = (elasTens(1,2) + elasTens(2,3) + elasTens(3,1) + & cEquiv_12 = (elasTens(1,2) + elasTens(2,3) + elasTens(3,1) + &
@ -955,6 +1077,8 @@ subroutine homogenization_RGC_equivalentModuli(&
cEquiv_44 = (elasTens(4,4) + elasTens(5,5) + elasTens(6,6))/3.0_pReal cEquiv_44 = (elasTens(4,4) + elasTens(5,5) + elasTens(6,6))/3.0_pReal
shearMod = 0.2_pReal*(cEquiv_11 - cEquiv_12) + 0.6_pReal*cEquiv_44 shearMod = 0.2_pReal*(cEquiv_11 - cEquiv_12) + 0.6_pReal*cEquiv_44
vBurgers = constitutive_averageBurgers(grainID,ip,el)
return return
endsubroutine endsubroutine
@ -998,23 +1122,36 @@ endsubroutine
subroutine homogenization_RGC_interfaceNormal(& subroutine homogenization_RGC_interfaceNormal(&
nVect, & ! interface normal nVect, & ! interface normal
! !
intFace & ! interface ID in 4D array (normal and position) intFace, & ! interface ID in 4D array (normal and position)
ip, & ! my IP
el & ! my element
) )
use prec, only: pReal,pInt,p_vec use prec, only: pReal,pInt,p_vec
use math, only: math_mul33x3
use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips
implicit none implicit none
!* Definition of variables !* Definition of variables
real(pReal), dimension (3), intent(out) :: nVect real(pReal), dimension (3), intent(out) :: nVect
integer(pInt), dimension (4), intent(in) :: intFace integer(pInt), dimension (4), intent(in) :: intFace
integer(pInt) nPos integer(pInt), intent(in) :: ip,el
integer(pInt) nPos,i
!* Get the normal of the interface, identified from the value of intFace(1) !* Get the normal of the interface, identified from the value of intFace(1)
nVect = 0.0_pReal nVect = 0.0_pReal
nPos = abs(intFace(1)) nPos = abs(intFace(1))
nVect(nPos) = intFace(1)/abs(intFace(1)) nVect(nPos) = intFace(1)/abs(intFace(1))
nVect = math_mul33x3(homogenization_RGC_orientation(:,:,ip,el),nVect)
! if (ip == 1 .and. el == 1) then
! write(6,'(x,a32,3(x,i3))')'Interface normal: ',intFace(1)
! write(6,'(x,3(e14.8,x))')(nVect(i), i = 1,3)
! write(6,*)' '
! call flush(6)
! endif
return return
endsubroutine endsubroutine
@ -1211,6 +1348,7 @@ subroutine homogenization_RGC_grainDeformation(&
F0, & ! initial partioned def grad per grain F0, & ! initial partioned def grad per grain
avgF, & ! my average def grad avgF, & ! my average def grad
state, & ! my state state, & ! my state
ip, & ! my IP
el & ! my element el & ! my element
) )
use prec, only: pReal,pInt,p_vec use prec, only: pReal,pInt,p_vec
@ -1223,7 +1361,7 @@ subroutine homogenization_RGC_grainDeformation(&
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: F0 real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: F0
real(pReal), dimension (3,3), intent(in) :: avgF real(pReal), dimension (3,3), intent(in) :: avgF
type(p_vec), intent(in) :: state type(p_vec), intent(in) :: state
integer(pInt), intent(in) :: el integer(pInt), intent(in) :: el,ip
! !
real(pReal), dimension (3) :: aVect,nVect real(pReal), dimension (3) :: aVect,nVect
integer(pInt), dimension (4) :: intFace integer(pInt), dimension (4) :: intFace
@ -1240,7 +1378,7 @@ subroutine homogenization_RGC_grainDeformation(&
do iFace = 1,nFace do iFace = 1,nFace
call homogenization_RGC_getInterface(intFace,iFace,iGrain3) call homogenization_RGC_getInterface(intFace,iFace,iGrain3)
call homogenization_RGC_relaxationVector(aVect,intFace,state,homID) call homogenization_RGC_relaxationVector(aVect,intFace,state,homID)
call homogenization_RGC_interfaceNormal(nVect,intFace) call homogenization_RGC_interfaceNormal(nVect,intFace,ip,el)
forall (i=1:3,j=1:3) & forall (i=1:3,j=1:3) &
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
enddo enddo

2
code/material.f90
View File

@ -535,6 +535,7 @@ subroutine material_populateGrains()
use math, only: math_sampleRandomOri, math_sampleGaussOri, math_sampleFiberOri, math_symmetricEulers, inDeg use math, only: math_sampleRandomOri, math_sampleGaussOri, math_sampleFiberOri, math_symmetricEulers, inDeg
use mesh, only: mesh_element, mesh_maxNips, mesh_NcpElems, mesh_ipVolume, FE_Nips use mesh, only: mesh_element, mesh_maxNips, mesh_NcpElems, mesh_ipVolume, FE_Nips
use IO, only: IO_error, IO_hybridIA use IO, only: IO_error, IO_hybridIA
use FEsolving, only: FEsolving_execIP
implicit none implicit none
integer(pInt), dimension (:,:), allocatable :: Ngrains integer(pInt), dimension (:,:), allocatable :: Ngrains
@ -711,6 +712,7 @@ subroutine material_populateGrains()
material_phase(g,i,e) = phaseOfGrain(grain+g) material_phase(g,i,e) = phaseOfGrain(grain+g)
material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+g) material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+g)
end forall end forall
FEsolving_execIP(2,e) = 1_pInt ! restrict calculation to first IP only, since all other results are to be copied from this
grain = grain + dGrains ! wind forward by NgrainsPerIP grain = grain + dGrains ! wind forward by NgrainsPerIP
else else
forall (i = 1:FE_Nips(mesh_element(2,e)), g = 1:dGrains) ! loop over IPs and grains forall (i = 1:FE_Nips(mesh_element(2,e)), g = 1:dGrains) ! loop over IPs and grains

26
code/math.f90
View File

@ -176,8 +176,6 @@ subroutine math_misorientation(dQ, Q1, Q2, symmetryType)
return return
endif endif
allocate(mySymOperations(4,NsymOperations(symmetryType)))
mySymOperations = symOperations(:,sum(NsymOperations(1:symmetryType-1))+1:sum(NsymOperations(1:symmetryType))) ! choose symmetry operations according to crystal symmetry
dQ(1) = -1.0_pReal ! start with maximum misorientation angle dQ(1) = -1.0_pReal ! start with maximum misorientation angle
do s = 1,NsymOperations(symmetryType) ! loop ver symmetry operations do s = 1,NsymOperations(symmetryType) ! loop ver symmetry operations
@ -999,6 +997,30 @@ pure function math_transpose3x3(A)
ENDFUNCTION ENDFUNCTION
!********************************************************************
! equivalent scalar quantity of a full strain tensor
!********************************************************************
PURE FUNCTION math_equivStrain33(m)
use prec, only: pReal,pInt
implicit none
real(pReal), dimension(3,3), intent(in) :: m
real(pReal) math_equivStrain33,e11,e22,e33,s12,s23,s31
e11 = (2.0_pReal*m(1,1)-m(2,2)-m(3,3))/3.0_pReal
e22 = (2.0_pReal*m(2,2)-m(3,3)-m(1,1))/3.0_pReal
e33 = (2.0_pReal*m(3,3)-m(1,1)-m(2,2))/3.0_pReal
s12 = 2.0_pReal*m(1,2)
s23 = 2.0_pReal*m(2,3)
s31 = 2.0_pReal*m(3,1)
math_equivStrain33 = 2.0_pReal*(1.50_pReal*(e11**2.0_pReal+e22**2.0_pReal+e33**2.0_pReal) + &
0.75_pReal*(s12**2.0_pReal+s23**2.0_pReal+s31**2.0_pReal))**(0.5_pReal)/3.0_pReal
return
ENDFUNCTION
!******************************************************************** !********************************************************************
! determinant of a 3x3 matrix ! determinant of a 3x3 matrix
!******************************************************************** !********************************************************************

2
code/mesh.f90
View File

@ -240,7 +240,7 @@
call mesh_build_ipAreas() call mesh_build_ipAreas()
call mesh_tell_statistics() call mesh_tell_statistics()
parallelExecution = (mesh_Nelems == mesh_NcpElems) ! plus potential STOP if non-local constitutive present parallelExecution = (parallelExecution .and. (mesh_Nelems == mesh_NcpElems)) ! plus potential killer from non-local constitutive
else else
call IO_error(101) ! cannot open input file call IO_error(101) ! cannot open input file
endif endif

1
code/mpie_cpfem_marc.f90
View File

@ -212,7 +212,6 @@ subroutine hypela2(&
outdatedByNewInc = .true. outdatedByNewInc = .true.
terminallyIll = .false. terminallyIll = .false.
cycleCounter = 0 cycleCounter = 0
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write (6,'(i6,x,i2,x,a)') n(1),nn,'<< hypela2 >> lastIncConverged + outdated'; call flush(6) write (6,'(i6,x,i2,x,a)') n(1),nn,'<< hypela2 >> lastIncConverged + outdated'; call flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)

15
code/numerics.f90
View File

@ -35,8 +35,9 @@ real(pReal) relevantStrain, & ! strain
relMax_RGC, & ! relative maximum of RGC residuum relMax_RGC, & ! relative maximum of RGC residuum
pPert_RGC, & ! perturbation for computing RGC penalty tangent pPert_RGC, & ! perturbation for computing RGC penalty tangent
xSmoo_RGC, & ! RGC penalty smoothing parameter (hyperbolic tangent) xSmoo_RGC, & ! RGC penalty smoothing parameter (hyperbolic tangent)
ratePower_RGC, & ! power (sensitivity rate) of numerical viscosity in RGC scheme viscPower_RGC, & ! power (sensitivity rate) of numerical viscosity in RGC scheme
viscModus_RGC, & ! stress modulus of RGC numerical viscosity viscModus_RGC, & ! stress modulus of RGC numerical viscosity
refRelaxRate_RGC, & ! reference relaxation rate in RGC viscosity
maxdRelax_RGC, & ! threshold of maximum relaxation vector increment (if exceed this then cutback) maxdRelax_RGC, & ! threshold of maximum relaxation vector increment (if exceed this then cutback)
maxVolDiscr_RGC, & ! threshold of maximum volume discrepancy allowed maxVolDiscr_RGC, & ! threshold of maximum volume discrepancy allowed
volDiscrMod_RGC, & ! stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint) volDiscrMod_RGC, & ! stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint)
@ -111,8 +112,9 @@ subroutine numerics_init()
relMax_RGC = 1.0e+2 relMax_RGC = 1.0e+2
pPert_RGC = 1.0e-7 pPert_RGC = 1.0e-7
xSmoo_RGC = 1.0e-5 xSmoo_RGC = 1.0e-5
ratePower_RGC = 1.0e+0 ! Newton viscosity (linear model) viscPower_RGC = 1.0e+0 ! Newton viscosity (linear model)
viscModus_RGC = 0.0e+0 ! No viscosity is applied viscModus_RGC = 0.0e+0 ! No viscosity is applied
refRelaxRate_RGC = 1.0e-3
maxdRelax_RGC = 1.0e+0 maxdRelax_RGC = 1.0e+0
maxVolDiscr_RGC = 1.0e-5 ! tolerance for volume discrepancy allowed maxVolDiscr_RGC = 1.0e-5 ! tolerance for volume discrepancy allowed
volDiscrMod_RGC = 1.0e+12 volDiscrMod_RGC = 1.0e+12
@ -189,10 +191,12 @@ subroutine numerics_init()
pPert_RGC = IO_floatValue(line,positions,2) pPert_RGC = IO_floatValue(line,positions,2)
case ('relevantmismatch_rgc') case ('relevantmismatch_rgc')
xSmoo_RGC = IO_floatValue(line,positions,2) xSmoo_RGC = IO_floatValue(line,positions,2)
case ('viscosityrate_rgc') case ('viscositypower_rgc')
ratePower_RGC = IO_floatValue(line,positions,2) viscPower_RGC = IO_floatValue(line,positions,2)
case ('viscositymodulus_rgc') case ('viscositymodulus_rgc')
viscModus_RGC = IO_floatValue(line,positions,2) viscModus_RGC = IO_floatValue(line,positions,2)
case ('refrelaxationrate_rgc')
refRelaxRate_RGC = IO_floatValue(line,positions,2)
case ('maxrelaxation_rgc') case ('maxrelaxation_rgc')
maxdRelax_RGC = IO_floatValue(line,positions,2) maxdRelax_RGC = IO_floatValue(line,positions,2)
case ('maxvoldiscrepancy_rgc') case ('maxvoldiscrepancy_rgc')
@ -292,8 +296,9 @@ subroutine numerics_init()
if (relMax_RGC <= 0.0_pReal) call IO_error(275) if (relMax_RGC <= 0.0_pReal) call IO_error(275)
if (pPert_RGC <= 0.0_pReal) call IO_error(276) !! oops !! if (pPert_RGC <= 0.0_pReal) call IO_error(276) !! oops !!
if (xSmoo_RGC <= 0.0_pReal) call IO_error(277) if (xSmoo_RGC <= 0.0_pReal) call IO_error(277)
if (ratePower_RGC < 0.0_pReal) call IO_error(278) if (viscPower_RGC < 0.0_pReal) call IO_error(278)
if (viscModus_RGC < 0.0_pReal) call IO_error(278) if (viscModus_RGC < 0.0_pReal) call IO_error(278)
if (refRelaxRate_RGC <= 0.0_pReal) call IO_error(278)
if (maxdRelax_RGC <= 0.0_pReal) call IO_error(288) if (maxdRelax_RGC <= 0.0_pReal) call IO_error(288)
if (maxVolDiscr_RGC <= 0.0_pReal) call IO_error(289) if (maxVolDiscr_RGC <= 0.0_pReal) call IO_error(289)
if (volDiscrMod_RGC < 0.0_pReal) call IO_error(289) if (volDiscrMod_RGC < 0.0_pReal) call IO_error(289)