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

No crystallite _converged

This commit is contained in:
Martin Diehl 2021-01-06 14:24:52 +01:00 committed by Sharan Roongta
parent f8dd5df0cc
commit 959c18c85e
18 changed files with 113 additions and 190 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
src/CPFEM.f90
View File

@ -5,7 +5,6 @@
!--------------------------------------------------------------------------------------------------
module CPFEM
use prec
use FEsolving
use math
use rotations
use YAML_types
@ -197,11 +196,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6)
else validCalculation
FEsolving_execElem = elCP
FEsolving_execIP = ip
if (debugCPFEM%extensive) &
print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip
call materialpoint_stressAndItsTangent(dt)
call materialpoint_stressAndItsTangent(dt,[ip,ip],[elCP,elCP])
terminalIllness: if (terminallyIll) then

1
src/CPFEM2.f90
View File

@ -6,7 +6,6 @@
module CPFEM2
use prec
use config
use FEsolving
use math
use rotations
use YAML_types

1
src/DAMASK_marc.f90
View File

@ -176,7 +176,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
use DAMASK_interface
use config
use YAML_types
use FEsolving
use discretization_marc
use homogenization
use CPFEM

15
src/FEsolving.f90
View File

@ -1,15 +0,0 @@
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief global variables for flow control
!--------------------------------------------------------------------------------------------------
module FEsolving
implicit none
public
integer, dimension(2) :: &
FEsolving_execElem, & !< for ping-pong scheme always whole range, otherwise one specific element
FEsolving_execIP !< for ping-pong scheme always range to max IP, otherwise one specific IP
end module FEsolving

1
src/commercialFEM_fileList.f90
View File

@ -13,7 +13,6 @@
#include "math.f90"
#include "quaternions.f90"
#include "rotations.f90"
#include "FEsolving.f90"
#include "element.f90"
#include "HDF5_utilities.f90"
#include "results.f90"

83
src/constitutive.f90
View File

@ -16,7 +16,6 @@ module constitutive
use parallelization
use HDF5_utilities
use DAMASK_interface
use FEsolving
use results
implicit none
@ -65,10 +64,6 @@ module constitutive
real(pReal), dimension(:,:,:,:,:), allocatable, public :: &
crystallite_partitionedF !< def grad to be reached at end of homog inc
logical, dimension(:,:,:), allocatable :: &
crystallite_converged !< convergence flag
type :: tTensorContainer
real(pReal), dimension(:,:,:), allocatable :: data
end type
@ -186,10 +181,10 @@ module constitutive
! == cleaned:end ===================================================================================
module function crystallite_stress(dt,co,ip,el)
module function crystallite_stress(dt,co,ip,el) result(converged_)
real(pReal), intent(in) :: dt
integer, intent(in) :: co, ip, el
logical :: crystallite_stress
logical :: converged_
end function crystallite_stress
module function constitutive_homogenizedC(co,ip,el) result(C)
@ -873,10 +868,8 @@ subroutine crystallite_init
source = crystallite_partitionedF)
allocate(crystallite_subdt(cMax,iMax,eMax),source=0.0_pReal)
allocate(crystallite_orientation(cMax,iMax,eMax))
allocate(crystallite_converged(cMax,iMax,eMax), source=.true.)
num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
@ -940,8 +933,8 @@ subroutine crystallite_init
flush(IO_STDOUT)
!$OMP PARALLEL DO PRIVATE(ph,me)
do el = FEsolving_execElem(1),FEsolving_execElem(2)
do ip = FEsolving_execIP(1), FEsolving_execIP(2); do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
do el = 1, size(material_phaseMemberAt,3)
do ip = 1, size(material_phaseMemberAt,2); do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
@ -967,14 +960,14 @@ subroutine crystallite_init
crystallite_partitionedF0 = crystallite_F0
crystallite_partitionedF = crystallite_F0
call crystallite_orientations()
!$OMP PARALLEL DO PRIVATE(ph,me)
do el = FEsolving_execElem(1),FEsolving_execElem(2)
do ip = FEsolving_execIP(1),FEsolving_execIP(2)
do el = 1, size(material_phaseMemberAt,3)
do ip = 1, size(material_phaseMemberAt,2)
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
call crystallite_orientations(co,ip,el)
call constitutive_plastic_dependentState(crystallite_partitionedF0(1:3,1:3,co,ip,el),co,ip,el) ! update dependent state variables to be consistent with basic states
enddo
enddo
@ -1089,7 +1082,7 @@ function crystallite_stressTangent(co,ip,el) result(dPdF)
el !< counter in element loop
integer :: &
o, &
p, pp, m
p, ph, me
real(pReal), dimension(3,3) :: devNull, &
invSubFp0,invSubFi0,invFp,invFi, &
@ -1109,19 +1102,19 @@ function crystallite_stressTangent(co,ip,el) result(dPdF)
real(pReal), dimension(9,9):: temp_99
logical :: error
pp = material_phaseAt(co,el)
m = material_phaseMemberAt(co,ip,el)
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
call constitutive_hooke_SandItsTangents(devNull,dSdFe,dSdFi, &
crystallite_Fe(1:3,1:3,co,ip,el), &
constitutive_mech_Fi(pp)%data(1:3,1:3,m),co,ip,el)
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el)
call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, &
crystallite_S (1:3,1:3,co,ip,el), &
constitutive_mech_Fi(pp)%data(1:3,1:3,m), &
constitutive_mech_Fi(ph)%data(1:3,1:3,me), &
co,ip,el)
invFp = math_inv33(constitutive_mech_Fp(pp)%data(1:3,1:3,m))
invFi = math_inv33(constitutive_mech_Fi(pp)%data(1:3,1:3,m))
invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me))
invFi = math_inv33(constitutive_mech_Fi(ph)%data(1:3,1:3,me))
invSubFp0 = math_inv33(crystallite_subFp0(1:3,1:3,co,ip,el))
invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,co,ip,el))
@ -1150,7 +1143,7 @@ function crystallite_stressTangent(co,ip,el) result(dPdF)
call constitutive_plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, &
crystallite_S (1:3,1:3,co,ip,el), &
constitutive_mech_Fi(pp)%data(1:3,1:3,m),co,ip,el)
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el)
dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS
!--------------------------------------------------------------------------------------------------
@ -1210,34 +1203,20 @@ end function crystallite_stressTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates orientations
!--------------------------------------------------------------------------------------------------
subroutine crystallite_orientations
subroutine crystallite_orientations(co,ip,el)
integer &
integer, intent(in) :: &
co, & !< counter in integration point component loop
ip, & !< counter in integration point loop
el !< counter in element loop
!$OMP PARALLEL DO
do el = FEsolving_execElem(1),FEsolving_execElem(2)
do ip = FEsolving_execIP(1),FEsolving_execIP(2)
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,co,ip,el))))
enddo; enddo; enddo
!$OMP END PARALLEL DO
call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,co,ip,el))))
if (plasticState(material_phaseAt(1,el))%nonlocal) &
call plastic_nonlocal_updateCompatibility(crystallite_orientation, &
phase_plasticityInstance(material_phaseAt(1,el)),ip,el)
nonlocalPresent: if (any(plasticState%nonlocal)) then
!$OMP PARALLEL DO
do el = FEsolving_execElem(1),FEsolving_execElem(2)
if (plasticState(material_phaseAt(1,el))%nonlocal) then
do ip = FEsolving_execIP(1),FEsolving_execIP(2)
call plastic_nonlocal_updateCompatibility(crystallite_orientation, &
phase_plasticityInstance(material_phaseAt(1,el)),ip,el)
enddo
endif
enddo
!$OMP END PARALLEL DO
endif nonlocalPresent
end subroutine crystallite_orientations
@ -1268,7 +1247,7 @@ end function crystallite_push33ToRef
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
subroutine integrateSourceState(co,ip,el)
function integrateSourceState(co,ip,el) result(broken)
integer, intent(in) :: &
el, & !< element index in element loop
@ -1288,12 +1267,13 @@ subroutine integrateSourceState(co,ip,el)
r ! state residuum
real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState
logical :: &
broken
broken, converged_
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
converged_ = .true.
broken = constitutive_thermal_collectDotState(ph,me)
broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,co,ip,el), co,ip,el,ph,me)
if(broken) return
@ -1328,19 +1308,20 @@ subroutine integrateSourceState(co,ip,el)
- sourceState(ph)%p(so)%dotState (1:size_so(so),me) * crystallite_subdt(co,ip,el)
sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%state(1:size_so(so),me) &
- r(1:size_so(so))
crystallite_converged(co,ip,el) = &
crystallite_converged(co,ip,el) .and. converged(r(1:size_so(so)), &
sourceState(ph)%p(so)%state(1:size_so(so),me), &
sourceState(ph)%p(so)%atol(1:size_so(so)))
converged_ = converged_ .and. converged(r(1:size_so(so)), &
sourceState(ph)%p(so)%state(1:size_so(so),me), &
sourceState(ph)%p(so)%atol(1:size_so(so)))
enddo
if(crystallite_converged(co,ip,el)) then
if(converged_) then
broken = constitutive_damage_deltaState(crystallite_Fe(1:3,1:3,co,ip,el),co,ip,el,ph,me)
exit iteration
endif
enddo iteration
broken = broken .or. .not. converged_
contains
@ -1364,7 +1345,7 @@ subroutine integrateSourceState(co,ip,el)
end function damper
end subroutine integrateSourceState
end function integrateSourceState
!--------------------------------------------------------------------------------------------------

73
src/constitutive_mech.f90
View File

@ -951,7 +951,7 @@ end function integrateStress
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
subroutine integrateStateFPI(F_0,F,Delta_t,co,ip,el)
function integrateStateFPI(F_0,F,Delta_t,co,ip,el) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F
real(pReal), intent(in) :: Delta_t
@ -1004,11 +1004,7 @@ subroutine integrateStateFPI(F_0,F,Delta_t,co,ip,el)
- plasticState(ph)%dotState (1:size_pl,me) * Delta_t
plasticState(ph)%state(1:size_pl,me) = plasticState(ph)%state(1:size_pl,me) &
- r(1:size_pl)
crystallite_converged(co,ip,el) = converged(r(1:size_pl), &
plasticState(ph)%state(1:size_pl,me), &
plasticState(ph)%atol(1:size_pl))
if(crystallite_converged(co,ip,el)) then
if (converged(r(1:size_pl),plasticState(ph)%state(1:size_pl,me),plasticState(ph)%atol(1:size_pl))) then
broken = constitutive_deltaState(crystallite_S(1:3,1:3,co,ip,el), &
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el,ph,me)
exit iteration
@ -1016,7 +1012,6 @@ subroutine integrateStateFPI(F_0,F,Delta_t,co,ip,el)
enddo iteration
contains
!--------------------------------------------------------------------------------------------------
@ -1039,13 +1034,13 @@ subroutine integrateStateFPI(F_0,F,Delta_t,co,ip,el)
end function damper
end subroutine integrateStateFPI
end function integrateStateFPI
!--------------------------------------------------------------------------------------------------
!> @brief integrate state with 1st order explicit Euler method
!--------------------------------------------------------------------------------------------------
subroutine integrateStateEuler(F_0,F,Delta_t,co,ip,el)
function integrateStateEuler(F_0,F,Delta_t,co,ip,el) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F
real(pReal), intent(in) :: Delta_t
@ -1075,15 +1070,14 @@ subroutine integrateStateEuler(F_0,F,Delta_t,co,ip,el)
if(broken) return
broken = integrateStress(F,Delta_t,co,ip,el)
crystallite_converged(co,ip,el) = .not. broken
end subroutine integrateStateEuler
end function integrateStateEuler
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with 1st order Euler method with adaptive step size
!--------------------------------------------------------------------------------------------------
subroutine integrateStateAdaptiveEuler(F_0,F,Delta_t,co,ip,el)
function integrateStateAdaptiveEuler(F_0,F,Delta_t,co,ip,el) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F
real(pReal), intent(in) :: Delta_t
@ -1123,24 +1117,22 @@ subroutine integrateStateAdaptiveEuler(F_0,F,Delta_t,co,ip,el)
broken = mech_collectDotState(Delta_t, co,ip,el,ph,me)
if(broken) return
broken = .not. converged(residuum_plastic(1:sizeDotState) + 0.5_pReal * plasticState(ph)%dotState(:,me) * Delta_t, &
plasticState(ph)%state(1:sizeDotState,me), &
plasticState(ph)%atol(1:sizeDotState))
sizeDotState = plasticState(ph)%sizeDotState
crystallite_converged(co,ip,el) = converged(residuum_plastic(1:sizeDotState) &
+ 0.5_pReal * plasticState(ph)%dotState(:,me) * Delta_t, &
plasticState(ph)%state(1:sizeDotState,me), &
plasticState(ph)%atol(1:sizeDotState))
end subroutine integrateStateAdaptiveEuler
end function integrateStateAdaptiveEuler
!---------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with the classic Runge Kutta method
!---------------------------------------------------------------------------------------------------
subroutine integrateStateRK4(F_0,F,Delta_t,co,ip,el)
function integrateStateRK4(F_0,F,Delta_t,co,ip,el) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: co,ip,el
logical :: broken
real(pReal), dimension(3,3), parameter :: &
A = reshape([&
@ -1153,19 +1145,20 @@ subroutine integrateStateRK4(F_0,F,Delta_t,co,ip,el)
real(pReal), dimension(4), parameter :: &
B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal]
call integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C)
broken = integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C)
end subroutine integrateStateRK4
end function integrateStateRK4
!---------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with the Cash-Carp method
!---------------------------------------------------------------------------------------------------
subroutine integrateStateRKCK45(F_0,F,Delta_t,co,ip,el)
function integrateStateRKCK45(F_0,F,Delta_t,co,ip,el) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: co,ip,el
logical :: broken
real(pReal), dimension(5,5), parameter :: &
A = reshape([&
@ -1185,16 +1178,16 @@ subroutine integrateStateRKCK45(F_0,F,Delta_t,co,ip,el)
[2825.0_pReal/27648.0_pReal, .0_pReal, 18575.0_pReal/48384.0_pReal,&
13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 1._pReal/4._pReal]
call integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB)
broken = integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB)
end subroutine integrateStateRKCK45
end function integrateStateRKCK45
!--------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an
!! embedded explicit Runge-Kutta method
!--------------------------------------------------------------------------------------------------
subroutine integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB)
function integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F
real(pReal), intent(in) :: Delta_t
@ -1205,15 +1198,14 @@ subroutine integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB)
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
logical :: broken
integer :: &
integer :: &
stage, & ! stage index in integration stage loop
n, &
ph, &
me, &
sizeDotState
logical :: &
broken
real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState
@ -1266,10 +1258,8 @@ subroutine integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB)
if(broken) return
broken = integrateStress(F,Delta_t,co,ip,el)
crystallite_converged(co,ip,el) = .not. broken
end subroutine integrateStateRK
end function integrateStateRK
!--------------------------------------------------------------------------------------------------
@ -1479,15 +1469,14 @@ end function constitutive_homogenizedC
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress (P)
!--------------------------------------------------------------------------------------------------
module function crystallite_stress(dt,co,ip,el)
module function crystallite_stress(dt,co,ip,el) result(converged_)
real(pReal), intent(in) :: dt
integer, intent(in) :: &
co, &
ip, &
el
logical :: crystallite_stress
logical :: converged_
real(pReal) :: &
formerSubStep
@ -1519,7 +1508,7 @@ module function crystallite_stress(dt,co,ip,el)
subFrac = 0.0_pReal
subStep = 1.0_pReal/num%subStepSizeCryst
todo = .true.
crystallite_converged(co,ip,el) = .false. ! pretend failed step of 1/subStepSizeCryst
converged_ = .false. ! pretend failed step of 1/subStepSizeCryst
todo = .true.
NiterationCrystallite = 0
@ -1528,7 +1517,7 @@ module function crystallite_stress(dt,co,ip,el)
!--------------------------------------------------------------------------------------------------
! wind forward
if (crystallite_converged(co,ip,el)) then
if (converged_) then
formerSubStep = subStep
subFrac = subFrac + subStep
subStep = min(1.0_pReal - subFrac, num%stepIncreaseCryst * subStep)
@ -1579,17 +1568,13 @@ module function crystallite_stress(dt,co,ip,el)
math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), &
constitutive_mech_Fp(ph)%data(1:3,1:3,me))))
crystallite_subdt(co,ip,el) = subStep * dt
crystallite_converged(co,ip,el) = .false.
call integrateState(subF0,crystallite_subF(1:3,1:3,co,ip,el),&
crystallite_subdt(co,ip,el),co,ip,el)
call integrateSourceState(co,ip,el)
converged_ = .not. integrateState(subF0,crystallite_subF(1:3,1:3,co,ip,el),&
crystallite_subdt(co,ip,el),co,ip,el)
converged_ = converged_ .and. .not. integrateSourceState(co,ip,el)
endif
enddo cutbackLooping
! return whether converged or not
crystallite_stress = crystallite_converged(co,ip,el)
end function crystallite_stress
end submodule constitutive_mech

4
src/grid/discretization_grid.f90
View File

@ -19,7 +19,6 @@ module discretization_grid
use results
use discretization
use geometry_plastic_nonlocal
use FEsolving
implicit none
private
@ -117,9 +116,6 @@ subroutine discretization_grid_init(restart)
(grid(1)+1) * (grid(2)+1) * grid3,& ! ...unless not last process
worldrank+1==worldsize))
FEsolving_execElem = [1,product(myGrid)] ! parallel loop bounds set to comprise all elements
FEsolving_execIP = [1,1] ! parallel loop bounds set to comprise the only IP
!--------------------------------------------------------------------------------------------------
! store geometry information for post processing
if(.not. restart) then

1
src/grid/grid_mech_FEM.f90
View File

@ -18,7 +18,6 @@ module grid_mech_FEM
use math
use rotations
use spectral_utilities
use FEsolving
use config
use homogenization
use discretization

1
src/grid/grid_mech_spectral_basic.f90
View File

@ -18,7 +18,6 @@ module grid_mech_spectral_basic
use math
use rotations
use spectral_utilities
use FEsolving
use config
use homogenization
use discretization_grid

1
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -18,7 +18,6 @@ module grid_mech_spectral_polarisation
use math
use rotations
use spectral_utilities
use FEsolving
use config
use homogenization
use discretization_grid

4
src/grid/spectral_utilities.f90
View File

@ -810,9 +810,9 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
print'(/,a)', ' ... evaluating constitutive response ......................................'
flush(IO_STDOUT)
homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field
homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field
call materialpoint_stressAndItsTangent(timeinc) ! calculate P field
call materialpoint_stressAndItsTangent(timeinc,[1,1],[1,product(grid(1:2))*grid3]) ! calculate P field
P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3])
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P

76
src/homogenization.f90
View File

@ -11,7 +11,6 @@ module homogenization
use math
use material
use constitutive
use FEsolving
use discretization
use thermal_isothermal
use thermal_conduction
@ -144,27 +143,29 @@ end subroutine homogenization_init
!--------------------------------------------------------------------------------------------------
!> @brief parallelized calculation of stress and corresponding tangent at material points
!--------------------------------------------------------------------------------------------------
subroutine materialpoint_stressAndItsTangent(dt)
subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execElem)
real(pReal), intent(in) :: dt !< time increment
integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP
integer :: &
NiterationHomog, &
NiterationMPstate, &
ip, & !< integration point number
el, & !< element number
myNgrains, co, ce
myNgrains, co, ce, ho
real(pReal) :: &
subFrac, &
subStep
logical :: &
requested, &
converged
logical, dimension(2) :: &
doneAndHappy
!$OMP PARALLEL DO PRIVATE(ce,myNgrains,NiterationMPstate,NiterationHomog,subFrac,converged,subStep,requested,doneAndHappy)
!$OMP PARALLEL DO PRIVATE(ce,ho,myNgrains,NiterationMPstate,NiterationHomog,subFrac,converged,subStep,doneAndHappy)
do el = FEsolving_execElem(1),FEsolving_execElem(2)
ho = material_homogenizationAt(el)
myNgrains = homogenization_Nconstituents(ho)
do ip = FEsolving_execIP(1),FEsolving_execIP(2)
!--------------------------------------------------------------------------------------------------
@ -174,21 +175,19 @@ subroutine materialpoint_stressAndItsTangent(dt)
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
if (homogState(material_homogenizationAt(el))%sizeState > 0) &
homogState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = &
homogState(material_homogenizationAt(el))%State0( :,material_homogenizationMemberAt(ip,el))
if (homogState(ho)%sizeState > 0) &
homogState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = &
homogState(ho)%State0( :,material_homogenizationMemberAt(ip,el))
if (damageState(ho)%sizeState > 0) &
damageState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = &
damageState(ho)%State0( :,material_homogenizationMemberAt(ip,el))
if (damageState(material_homogenizationAt(el))%sizeState > 0) &
damageState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = &
damageState(material_homogenizationAt(el))%State0( :,material_homogenizationMemberAt(ip,el))
NiterationHomog = 0
cutBackLooping: do while (.not. terminallyIll .and. subStep > num%subStepMinHomog)
myNgrains = homogenization_Nconstituents(material_homogenizationAt(el))
if (converged) then
subFrac = subFrac + subStep
subStep = min(1.0_pReal-subFrac,num%stepIncreaseHomog*subStep) ! introduce flexibility for step increase/acceleration
@ -198,22 +197,20 @@ subroutine materialpoint_stressAndItsTangent(dt)
! wind forward grain starting point
call constitutive_windForward(ip,el)
if(homogState(material_homogenizationAt(el))%sizeState > 0) &
homogState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = &
homogState(material_homogenizationAt(el))%State (:,material_homogenizationMemberAt(ip,el))
if(damageState(material_homogenizationAt(el))%sizeState > 0) &
damageState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = &
damageState(material_homogenizationAt(el))%State (:,material_homogenizationMemberAt(ip,el))
if(homogState(ho)%sizeState > 0) &
homogState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = &
homogState(ho)%State (:,material_homogenizationMemberAt(ip,el))
if(damageState(ho)%sizeState > 0) &
damageState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = &
damageState(ho)%State (:,material_homogenizationMemberAt(ip,el))
endif steppingNeeded
else
if ( (myNgrains == 1 .and. subStep <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite
num%subStepSizeHomog * subStep <= num%subStepMinHomog ) then ! would require too small subStep
! cutback makes no sense
if (.not. terminallyIll) then ! so first signals terminally ill...
if (.not. terminallyIll) & ! so first signals terminally ill...
print*, ' Integration point ', ip,' at element ', el, ' terminally ill'
endif
terminallyIll = .true. ! ...and kills all others
else ! cutback makes sense
subStep = num%subStepSizeHomog * subStep ! crystallite had severe trouble, so do a significant cutback
@ -221,23 +218,19 @@ subroutine materialpoint_stressAndItsTangent(dt)
call crystallite_restore(ip,el,subStep < 1.0_pReal)
call constitutive_restore(ip,el)
if(homogState(material_homogenizationAt(el))%sizeState > 0) &
homogState(material_homogenizationAt(el))%State( :,material_homogenizationMemberAt(ip,el)) = &
homogState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el))
if(damageState(material_homogenizationAt(el))%sizeState > 0) &
damageState(material_homogenizationAt(el))%State( :,material_homogenizationMemberAt(ip,el)) = &
damageState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el))
if(homogState(ho)%sizeState > 0) &
homogState(ho)%State( :,material_homogenizationMemberAt(ip,el)) = &
homogState(ho)%subState0(:,material_homogenizationMemberAt(ip,el))
if(damageState(ho)%sizeState > 0) &
damageState(ho)%State( :,material_homogenizationMemberAt(ip,el)) = &
damageState(ho)%subState0(:,material_homogenizationMemberAt(ip,el))
endif
endif
if (subStep > num%subStepMinHomog) then
requested = .true.
doneAndHappy = [.false.,.true.]
endif
if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.]
NiterationMPstate = 0
convergenceLooping: do while (.not. terminallyIll .and. requested &
convergenceLooping: do while (.not. terminallyIll &
.and. .not. doneAndHappy(1) &
.and. NiterationMPstate < num%nMPstate)
NiterationMPstate = NiterationMPstate + 1
@ -245,7 +238,7 @@ subroutine materialpoint_stressAndItsTangent(dt)
!--------------------------------------------------------------------------------------------------
! deformation partitioning
if(requested .and. .not. doneAndHappy(1)) then ! requested but not yet done
if (.not. doneAndHappy(1)) then
ce = (el-1)*discretization_nIPs + ip
call mech_partition(homogenization_F0(1:3,1:3,ce) &
+ (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce))&
@ -255,10 +248,7 @@ subroutine materialpoint_stressAndItsTangent(dt)
do co = 1, myNgrains
converged = converged .and. crystallite_stress(dt*subStep,co,ip,el)
enddo
endif
if (requested .and. .not. doneAndHappy(1)) then
if (.not. converged) then
doneAndHappy = [.true.,.false.]
else
@ -281,10 +271,14 @@ subroutine materialpoint_stressAndItsTangent(dt)
!$OMP END PARALLEL DO
if (.not. terminallyIll ) then
call crystallite_orientations() ! calculate crystal orientations
!$OMP PARALLEL DO
!$OMP PARALLEL DO PRIVATE(ho,myNgrains)
elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
ho = material_homogenizationAt(el)
myNgrains = homogenization_Nconstituents(ho)
IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2)
do co = 1, myNgrains
call crystallite_orientations(co,ip,el)
enddo
call mech_homogenize(ip,el)
enddo IpLooping3
enddo elementLooping3

24
src/homogenization_mech.f90
View File

@ -128,35 +128,35 @@ module subroutine mech_homogenize(ip,el)
integer, intent(in) :: &
ip, & !< integration point
el !< element number
integer :: c,m
integer :: co,ce
real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
m = (el-1)* discretization_nIPs + ip
ce = (el-1)* discretization_nIPs + ip
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization
homogenization_P(1:3,1:3,m) = crystallite_P(1:3,1:3,1,ip,el)
homogenization_dPdF(1:3,1:3,1:3,1:3,m) = crystallite_stressTangent(1,ip,el)
homogenization_P(1:3,1:3,ce) = crystallite_P(1:3,1:3,1,ip,el)
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = crystallite_stressTangent(1,ip,el)
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
do c = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el)
enddo
call mech_isostrain_averageStressAndItsTangent(&
homogenization_P(1:3,1:3,m), &
homogenization_dPdF(1:3,1:3,1:3,1:3,m),&
homogenization_P(1:3,1:3,ce), &
homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
dPdFs, &
homogenization_typeInstance(material_homogenizationAt(el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
do c = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el)
enddo
call mech_RGC_averageStressAndItsTangent(&
homogenization_P(1:3,1:3,m), &
homogenization_dPdF(1:3,1:3,1:3,1:3,m),&
homogenization_P(1:3,1:3,ce), &
homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
dPdFs, &
homogenization_typeInstance(material_homogenizationAt(el)))

4
src/marc/discretization_marc.f90
View File

@ -12,7 +12,6 @@ module discretization_marc
use DAMASK_interface
use IO
use config
use FEsolving
use element
use discretization
use geometry_plastic_nonlocal
@ -89,9 +88,6 @@ subroutine discretization_marc_init
if (debug_e < 1 .or. debug_e > nElems) call IO_error(602,ext_msg='element')
if (debug_i < 1 .or. debug_i > elem%nIPs) call IO_error(602,ext_msg='IP')
FEsolving_execElem = [1,nElems]
FEsolving_execIP = [1,elem%nIPs]
allocate(cellNodeDefinition(elem%nNodes-1))
allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,nElems))
call buildCells(connectivity_cell,cellNodeDefinition,&

1
src/mesh/DAMASK_mesh.f90
View File

@ -15,7 +15,6 @@ program DAMASK_mesh
use IO
use math
use CPFEM2
use FEsolving
use config
use discretization_mesh
use FEM_Utilities

2
src/mesh/FEM_utilities.f90
View File

@ -160,7 +160,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
print'(/,a)', ' ... evaluating constitutive response ......................................'
call materialpoint_stressAndItsTangent(timeinc) ! calculate P field
call materialpoint_stressAndItsTangent(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field
cutBack = .false. ! reset cutBack status

6
src/mesh/discretization_mesh.f90
View File

@ -18,7 +18,6 @@ module discretization_mesh
use config
use discretization
use results
use FEsolving
use FEM_quadrature
use YAML_types
use prec
@ -30,7 +29,7 @@ module discretization_mesh
mesh_Nboundaries, &
mesh_NcpElemsGlobal
integer :: &
integer, public, protected :: &
mesh_NcpElems !< total number of CP elements in mesh
!!!! BEGIN DEPRECATED !!!!!
@ -174,9 +173,6 @@ subroutine discretization_mesh_init(restart)
if (debug_element < 1 .or. debug_element > mesh_NcpElems) call IO_error(602,ext_msg='element')
if (debug_ip < 1 .or. debug_ip > mesh_maxNips) call IO_error(602,ext_msg='IP')
FEsolving_execElem = [1,mesh_NcpElems] ! parallel loop bounds set to comprise all DAMASK elements
FEsolving_execIP = [1,mesh_maxNips]
allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal)
call discretization_init(materialAt,&