DAMASK structure updated; define debug parameters when used by a module

This commit is contained in:
Sharan Roongta 2020-06-18 16:06:11 +02:00
parent 4e60d8e133
commit c987f55f69
5 changed files with 107 additions and 65 deletions

View File

@ -54,15 +54,20 @@ module homogenization
interface
module subroutine mech_none_init
module subroutine mech_none_init(debug_homogenization)
class(tNode), pointer, intent(in) :: &
debug_homogenization
end subroutine mech_none_init
module subroutine mech_isostrain_init
module subroutine mech_isostrain_init(debug_homogenization)
class(tNode), pointer, intent(in) :: &
debug_homogenization
end subroutine mech_isostrain_init
module subroutine mech_RGC_init(num_homogMech)
module subroutine mech_RGC_init(num_homogMech, debug_homogenization)
class(tNode), pointer, intent(in) :: &
num_homogMech
num_homogMech, &
debug_homogenization
end subroutine mech_RGC_init
@ -71,12 +76,15 @@ module homogenization
real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point
end subroutine mech_isostrain_partitionDeformation
module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of, &
debug_homogenization)
real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient
real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point
integer, intent(in) :: &
instance, &
of
class(tNode), pointer, intent(in) :: &
debug_homogenization
end subroutine mech_RGC_partitionDeformation
@ -98,8 +106,7 @@ module homogenization
integer, intent(in) :: instance
end subroutine mech_RGC_averageStressAndItsTangent
module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el,debug_homogenization)
logical, dimension(2) :: mech_RGC_updateState
real(pReal), dimension(:,:,:), intent(in) :: &
P,& !< partitioned stresses
@ -111,6 +118,8 @@ module homogenization
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
class(tNode), pointer, intent(in) :: &
debug_homogenization
end function mech_RGC_updateState
@ -137,15 +146,21 @@ subroutine homogenization_init
class (tNode) , pointer :: &
num_homog, &
num_homogMech, &
num_homogGeneric
num_homogGeneric, &
debug_homogenization
integer :: &
debug_g, &
debug_e
num_homog => numerics_root%get('homogenization',defaultVal=emptyDict)
num_homogMech => num_homog%get('mech',defaultVal=emptyDict)
num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict)
if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mech_none_init
if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mech_isostrain_init
if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call mech_RGC_init(num_homogMech)
debug_homogenization => debug_root%get('homogenization',defaultVal=emptyList)
if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mech_none_init(debug_homogenization)
if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mech_isostrain_init(debug_homogenization)
if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call mech_RGC_init(num_homogMech,debug_homogenization)
if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init
if (any(thermal_type == THERMAL_adiabatic_ID)) call thermal_adiabatic_init
@ -166,6 +181,8 @@ subroutine homogenization_init
write(6,'(/,a)') ' <<<+- homogenization init -+>>>'; flush(6)
debug_g = debug_root%get_asInt('grain', defaultVal=1)
debug_e = debug_root%get_asInt('element', defaultVal=1)
if (debug_g < 1 .or. debug_g > homogenization_Ngrains(material_homogenizationAt(debug_e))) &
call IO_error(602,ext_msg='constituent', el=debug_e, g=debug_g)
@ -197,7 +214,9 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
i, & !< integration point number
e, & !< element number
mySource, &
myNgrains
myNgrains, &
debug_e, &
debug_i
real(pReal), dimension(discretization_nIP,discretization_nElem) :: &
subFrac, &
subStep
@ -206,9 +225,15 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
converged
logical, dimension(2,discretization_nIP,discretization_nElem) :: &
doneAndHappy
class(tNode), pointer :: &
debug_homogenization
#ifdef DEBUG
if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0) then
debug_e = debug_root%get_asInt('element', defaultVal=1)
debug_i = debug_root%get_asInt('integrationpoint',defaultVal=1)
debug_homogenization => debug_root%get('homogenization',defaultVal=emptyList)
if (debug_homogenization%contains('basic')) then
write(6,'(/a,i5,1x,i2)') '<< HOMOG >> Material Point start at el ip ', debug_e, debug_i
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F0', &
@ -273,9 +298,9 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
if (converged(i,e)) then
#ifdef DEBUG
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0 &
if (debug_homogenization%contains('extensive') &
.and. ((e == debug_e .and. i == debug_i) &
.or. .not. iand(debug_level(debug_homogenization),debug_levelSelective) /= 0)) then
.or. .not. debug_homogenization%contains('selective'))) then
write(6,'(a,1x,f12.8,1x,a,1x,f12.8,1x,a,i8,1x,i2/)') '<< HOMOG >> winding forward from', &
subFrac(i,e), 'to current subFrac', &
subFrac(i,e)+subStep(i,e),'in materialpoint_stressAndItsTangent at el ip',e,i
@ -332,9 +357,9 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
subStep(i,e) = num%subStepSizeHomog * subStep(i,e) ! crystallite had severe trouble, so do a significant cutback
#ifdef DEBUG
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0 &
if (debug_homogenization%contains('extensive') &
.and. ((e == debug_e .and. i == debug_i) &
.or. .not. iand(debug_level(debug_homogenization), debug_levelSelective) /= 0)) then
.or. .not. debug_homogenization%contains('selective'))) then
write(6,'(a,1x,f12.8,a,i8,1x,i2/)') &
'<< HOMOG >> cutback step in materialpoint_stressAndItsTangent with new subStep:',&
subStep(i,e),' at el ip',e,i

View File

@ -75,16 +75,19 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields, reads information from material configuration file
!--------------------------------------------------------------------------------------------------
module subroutine mech_RGC_init(num_homogMech)
module subroutine mech_RGC_init(num_homogMech,debug_homogenization)
class(tNode), pointer, intent(in) :: &
num_homogMech
num_homogMech, &
debug_homogenization
integer :: &
Ninstance, &
h, &
NofMyHomog, &
sizeState, nIntFaceTot
sizeState, nIntFaceTot, &
debug_e, &
debug_i
class (tNode), pointer :: &
num_RGC
@ -98,7 +101,7 @@ module subroutine mech_RGC_init(num_homogMech)
write(6,'(a)') ' https://doi.org/10.1088/0965-0393/18/1/015006'
Ninstance = count(homogenization_type == HOMOGENIZATION_RGC_ID)
if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0) &
if (debug_homogenization%contains('basic')) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(param(Ninstance))
@ -146,6 +149,8 @@ module subroutine mech_RGC_init(num_homogMech)
config => config_homogenization(h))
#ifdef DEBUG
debug_e = debug_root%get_asInt('element',defaultVal=1)
debug_i = debug_root%get_asInt('integrationpoint',defaultVal=1)
if (h==material_homogenizationAt(debug_e)) then
prm%of_debug = material_homogenizationMemberAt(debug_i,debug_e)
endif
@ -200,7 +205,7 @@ end subroutine mech_RGC_init
!--------------------------------------------------------------------------------------------------
!> @brief partitions the deformation gradient onto the constituents
!--------------------------------------------------------------------------------------------------
module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of,debug_homogenization)
real(pReal), dimension (:,:,:), intent(out) :: F !< partioned F per grain
@ -208,6 +213,8 @@ module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
integer, intent(in) :: &
instance, &
of
class(tNode), pointer, intent(in) :: &
debug_homogenization
real(pReal), dimension(3) :: aVect,nVect
integer, dimension(4) :: intFace
@ -231,7 +238,7 @@ module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! resulting relaxed deformation gradient
#ifdef DEBUG
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0) then
if (debug_homogenization%contains('extensive')) then
write(6,'(1x,a32,1x,i3)')'Deformation gradient of grain: ',iGrain
do i = 1,3
write(6,'(1x,3(e15.8,1x))')(F(i,j,iGrain), j = 1,3)
@ -294,7 +301,7 @@ module procedure mech_RGC_updateState
drelax = stt%relaxationVector(:,of) - st0%relaxationVector(:,of)
#ifdef DEBUG
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0) then
if (debug_homogenization%contains('extensive')) then
write(6,'(1x,a30)')'Obtained state: '
do i = 1,size(stt%relaxationVector(:,of))
write(6,'(1x,2(e15.8,1x))') stt%relaxationVector(i,of)
@ -305,14 +312,14 @@ module procedure mech_RGC_updateState
!--------------------------------------------------------------------------------------------------
! computing interface mismatch and stress penalty tensor for all interfaces of all grains
call stressPenalty(R,NN,avgF,F,ip,el,instance,of)
call stressPenalty(R,NN,avgF,F,ip,el,instance,of,debug_homogenization)
!--------------------------------------------------------------------------------------------------
! calculating volume discrepancy and stress penalty related to overall volume discrepancy
call volumePenalty(D,dst%volumeDiscrepancy(of),avgF,F,nGrain,instance,of)
call volumePenalty(D,dst%volumeDiscrepancy(of),avgF,F,nGrain,instance,of,debug_homogenization)
#ifdef DEBUG
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0) then
if (debug_homogenization%contains('extensive')) then
do iGrain = 1,nGrain
write(6,'(1x,a30,1x,i3,1x,a4,3(1x,e15.8))')'Mismatch magnitude of grain(',iGrain,') :',&
NN(1,iGrain),NN(2,iGrain),NN(3,iGrain)
@ -360,7 +367,7 @@ module procedure mech_RGC_updateState
enddo
#ifdef DEBUG
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0) then
if (debug_homogenization%contains('extensive')) then
write(6,'(1x,a30,1x,i3)')'Traction at interface: ',iNum
write(6,'(1x,3(e15.8,1x))')(tract(iNum,j), j = 1,3)
write(6,*)' '
@ -374,7 +381,7 @@ module procedure mech_RGC_updateState
residMax = maxval(abs(tract)) ! get the maximum of the residual
#ifdef DEBUG
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 .and. prm%of_debug == of) then
if (debug_homogenization%contains('extensive') .and. prm%of_debug == of) then
stresLoc = maxloc(abs(P))
residLoc = maxloc(abs(tract))
write(6,'(1x,a)')' '
@ -394,7 +401,7 @@ module procedure mech_RGC_updateState
if (residMax < num%rtol*stresMax .or. residMax < num%atol) then
mech_RGC_updateState = .true.
#ifdef DEBUG
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 .and. prm%of_debug == of) &
if (debug_homogenization%contains('extensive') .and. prm%of_debug == of) &
write(6,'(1x,a55,/)')'... done and happy'; flush(6)
#endif
@ -414,7 +421,7 @@ module procedure mech_RGC_updateState
dst%relaxationRate_max(of) = maxval(abs(drelax))/dt
#ifdef DEBUG
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 .and. prm%of_debug == of) then
if (debug_homogenization%contains('extensive') .and. prm%of_debug == of) then
write(6,'(1x,a30,1x,e15.8)') 'Constitutive work: ',stt%work(of)
write(6,'(1x,a30,3(1x,e15.8))')'Magnitude mismatch: ',dst%mismatch(1,of), &
dst%mismatch(2,of), &
@ -435,7 +442,7 @@ module procedure mech_RGC_updateState
mech_RGC_updateState = [.true.,.false.] ! with direct cut-back
#ifdef DEBUG
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 .and. prm%of_debug == of) &
if (debug_homogenization%contains('extensive') .and. prm%of_debug == of) &
write(6,'(1x,a,/)') '... broken'; flush(6)
#endif
@ -443,7 +450,7 @@ module procedure mech_RGC_updateState
else ! proceed with computing the Jacobian and state update
#ifdef DEBUG
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 .and. prm%of_debug == of) &
if (debug_homogenization%contains('extensive') .and. prm%of_debug == of) &
write(6,'(1x,a,/)') '... not yet done'; flush(6)
#endif
@ -500,7 +507,7 @@ module procedure mech_RGC_updateState
enddo
#ifdef DEBUG
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0) then
if (debug_homogenization%contains('extensive')) then
write(6,'(1x,a30)')'Jacobian matrix of stress'
do i = 1,3*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(smatrix(i,j), j = 1,3*nIntFaceTot)
@ -522,8 +529,8 @@ module procedure mech_RGC_updateState
p_relax(ipert) = relax(ipert) + num%pPert ! perturb the relaxation vector
stt%relaxationVector(:,of) = p_relax
call grainDeformation(pF,avgF,instance,of) ! rain deformation from perturbed state
call stressPenalty(pR,DevNull, avgF,pF,ip,el,instance,of) ! stress penalty due to interface mismatch from perturbed state
call volumePenalty(pD,devNull(1,1), avgF,pF,nGrain,instance,of) ! stress penalty due to volume discrepancy from perturbed state
call stressPenalty(pR,DevNull, avgF,pF,ip,el,instance,of,debug_homogenization) ! stress penalty due to interface mismatch from perturbed state
call volumePenalty(pD,devNull(1,1), avgF,pF,nGrain,instance,of,debug_homogenization) ! stress penalty due to volume discrepancy from perturbed state
!--------------------------------------------------------------------------------------------------
! computing the global stress residual array from the perturbed state
@ -560,7 +567,7 @@ module procedure mech_RGC_updateState
enddo
#ifdef DEBUG
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0) then
if (debug_homogenization%contains('extensive')) then
write(6,'(1x,a30)')'Jacobian matrix of penalty'
do i = 1,3*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(pmatrix(i,j), j = 1,3*nIntFaceTot)
@ -579,7 +586,7 @@ module procedure mech_RGC_updateState
enddo
#ifdef DEBUG
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0) then
if (debug_homogenization%contains('extensive')) then
write(6,'(1x,a30)')'Jacobian matrix of penalty'
do i = 1,3*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(rmatrix(i,j), j = 1,3*nIntFaceTot)
@ -594,7 +601,7 @@ module procedure mech_RGC_updateState
allocate(jmatrix(3*nIntFaceTot,3*nIntFaceTot)); jmatrix = smatrix + pmatrix + rmatrix
#ifdef DEBUG
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0) then
if (debug_homogenization%contains('extensive')) then
write(6,'(1x,a30)')'Jacobian matrix (total)'
do i = 1,3*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(jmatrix(i,j), j = 1,3*nIntFaceTot)
@ -610,7 +617,7 @@ module procedure mech_RGC_updateState
call math_invert(jnverse,error,jmatrix)
#ifdef DEBUG
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0) then
if (debug_homogenization%contains('extensive')) then
write(6,'(1x,a30)')'Jacobian inverse'
do i = 1,3*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(jnverse(i,j), j = 1,3*nIntFaceTot)
@ -637,7 +644,7 @@ module procedure mech_RGC_updateState
endif
#ifdef DEBUG
if (iand(debug_homogenization, debug_levelExtensive) > 0) then
if (debug_homogenization%contains('extensive')) then
write(6,'(1x,a30)')'Returned state: '
do i = 1,size(stt%relaxationVector(:,of))
write(6,'(1x,2(e15.8,1x))') stt%relaxationVector(i,of)
@ -653,7 +660,7 @@ module procedure mech_RGC_updateState
!------------------------------------------------------------------------------------------------
!> @brief calculate stress-like penalty due to deformation mismatch
!------------------------------------------------------------------------------------------------
subroutine stressPenalty(rPen,nMis,avgF,fDef,ip,el,instance,of)
subroutine stressPenalty(rPen,nMis,avgF,fDef,ip,el,instance,of,debug_homogenization)
real(pReal), dimension (:,:,:), intent(out) :: rPen !< stress-like penalty
real(pReal), dimension (:,:), intent(out) :: nMis !< total amount of mismatch
@ -661,6 +668,7 @@ module procedure mech_RGC_updateState
real(pReal), dimension (:,:,:), intent(in) :: fDef !< deformation gradients
real(pReal), dimension (3,3), intent(in) :: avgF !< initial effective stretch tensor
integer, intent(in) :: ip,el,instance,of
class(tNode), pointer, intent(in) :: debug_homogenization
integer, dimension (4) :: intFace
integer, dimension (3) :: iGrain3,iGNghb3,nGDim
@ -687,7 +695,7 @@ module procedure mech_RGC_updateState
associate(prm => param(instance))
#ifdef DEBUG
debugActive = iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 .and. prm%of_debug == of
debugActive = debug_homogenization%contains('extensive') .and. prm%of_debug == of
if (debugActive) then
write(6,'(1x,a20,2(1x,i3))')'Correction factor: ',ip,el
@ -764,7 +772,7 @@ module procedure mech_RGC_updateState
!------------------------------------------------------------------------------------------------
!> @brief calculate stress-like penalty due to volume discrepancy
!------------------------------------------------------------------------------------------------
subroutine volumePenalty(vPen,vDiscrep,fAvg,fDef,nGrain,instance,of)
subroutine volumePenalty(vPen,vDiscrep,fAvg,fDef,nGrain,instance,of,debug_homogenization)
real(pReal), dimension (:,:,:), intent(out) :: vPen ! stress-like penalty due to volume
real(pReal), intent(out) :: vDiscrep ! total volume discrepancy
@ -775,6 +783,7 @@ module procedure mech_RGC_updateState
Ngrain, &
instance, &
of
class(tNode), pointer, intent(in) :: debug_homogenization
real(pReal), dimension(size(vPen,3)) :: gVol
integer :: i
@ -797,7 +806,7 @@ module procedure mech_RGC_updateState
gVol(i)*transpose(math_inv33(fDef(:,:,i)))
#ifdef DEBUG
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 &
if (debug_homogenization%contains('extensive') &
.and. param(instance)%of_debug == of) then
write(6,'(1x,a30,i2)')'Volume penalty of grain: ',i
write(6,*) transpose(vPen(:,:,i))

View File

@ -26,7 +26,10 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields, reads information from material configuration file
!--------------------------------------------------------------------------------------------------
module subroutine mech_isostrain_init
module subroutine mech_isostrain_init(debug_homogenization)
class(tNode), pointer, intent(in) :: &
debug_homogenization
integer :: &
Ninstance, &
@ -38,7 +41,7 @@ module subroutine mech_isostrain_init
write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_ISOSTRAIN_LABEL//' init -+>>>'
Ninstance = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)
if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0) &
if (debug_homogenization%contains('basic')) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(param(Ninstance)) ! one container of parameters per instance

View File

@ -11,7 +11,10 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields, reads information from material configuration file
!--------------------------------------------------------------------------------------------------
module subroutine mech_none_init
module subroutine mech_none_init(debug_homogenization)
class(tNode), pointer, intent(in) :: &
debug_homogenization
integer :: &
Ninstance, &
@ -21,7 +24,7 @@ module subroutine mech_none_init
write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_NONE_label//' init -+>>>'; flush(6)
Ninstance = count(homogenization_type == HOMOGENIZATION_NONE_ID)
if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0) &
if (debug_homogenization%contains('basic')) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
do h = 1, size(homogenization_type)

View File

@ -8,6 +8,7 @@ module material
use prec
use math
use config
use YAML_types
use results
use IO
use debug
@ -215,22 +216,23 @@ subroutine material_init(restart)
integer, dimension(:), allocatable :: &
CounterPhase, &
CounterHomogenization
myDebug = debug_level(debug_material)
class(tNode), pointer :: &
debug_material
write(6,'(/,a)') ' <<<+- material init -+>>>'; flush(6)
debug_material => debug_root%get('material',defaultVal=emptyList)
call material_parsePhase()
if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Phase parsed'; flush(6)
if (debug_material%contains('basic')) write(6,'(a)') ' Phase parsed'; flush(6)
call material_parseMicrostructure()
if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Microstructure parsed'; flush(6)
if (debug_material%contains('basic')) write(6,'(a)') ' Microstructure parsed'; flush(6)
call material_parseHomogenization()
if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Homogenization parsed'; flush(6)
if (debug_material%contains('basic')) write(6,'(a)') ' Homogenization parsed'; flush(6)
call material_parseTexture()
if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Texture parsed'; flush(6)
if (debug_material%contains('basic')) write(6,'(a)') ' Texture parsed'; flush(6)
material_Nphase = size(config_phase)
material_Nhomogenization = size(config_homogenization)
@ -266,7 +268,7 @@ subroutine material_init(restart)
enddo
if(homogenization_maxNgrains > size(microstructure_phase,1)) call IO_error(148)
debugOut: if (iand(myDebug,debug_levelExtensive) /= 0) then
debugOut: if (debug_material%contains('extensive')) then
write(6,'(/,a,/)') ' MATERIAL configuration'
write(6,'(a32,1x,a16,1x,a6)') 'homogenization ','type ','grains'
do h = 1,size(config_homogenization)