don't clutter with statements that are never used

This commit is contained in:
Martin Diehl 2020-12-23 14:10:26 +01:00
parent 972e041f59
commit 8ac880c0ad
2 changed files with 4 additions and 204 deletions

View File

@ -48,20 +48,6 @@ module homogenization
type(tNumerics) :: num type(tNumerics) :: num
type :: tDebugOptions
logical :: &
basic, &
extensive, &
selective
integer :: &
element, &
ip, &
grain
end type tDebugOptions
type(tDebugOptions) :: debugHomog
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
interface interface
@ -125,24 +111,10 @@ subroutine homogenization_init
class (tNode) , pointer :: & class (tNode) , pointer :: &
num_homog, & num_homog, &
num_homogGeneric, & num_homogGeneric
debug_homogenization
print'(/,a)', ' <<<+- homogenization init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- homogenization init -+>>>'; flush(IO_STDOUT)
debug_homogenization => config_debug%get('homogenization', defaultVal=emptyList)
debugHomog%basic = debug_homogenization%contains('basic')
debugHomog%extensive = debug_homogenization%contains('extensive')
debugHomog%selective = debug_homogenization%contains('selective')
debugHomog%element = config_debug%get_asInt('element',defaultVal = 1)
debugHomog%ip = config_debug%get_asInt('integrationpoint',defaultVal = 1)
debugHomog%grain = config_debug%get_asInt('grain',defaultVal = 1)
if (debugHomog%grain < 1 &
.or. debugHomog%grain > homogenization_Nconstituents(material_homogenizationAt(debugHomog%element))) &
call IO_error(602,ext_msg='constituent', el=debugHomog%element, g=debugHomog%grain)
num_homog => config_numerics%get('homogenization',defaultVal=emptyDict) num_homog => config_numerics%get('homogenization',defaultVal=emptyDict)
num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict) num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict)

View File

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