store data where it is needed
avoid globals, even if they are read only
This commit is contained in:
Martin Diehl
parent
0bb7fea782
commit
e4792e56fb
|
|
@ -23,7 +23,7 @@ submodule(homogenization) homogenization_mech_RGC
|
|||
character(len=pStringLen), allocatable, dimension(:) :: &
|
||||
output
|
||||
end type tParameters
|
||||
|
||||
|
||||
type :: tRGCstate
|
||||
real(pReal), pointer, dimension(:) :: &
|
||||
work, &
|
||||
|
|
@ -31,7 +31,7 @@ submodule(homogenization) homogenization_mech_RGC
|
|||
real(pReal), pointer, dimension(:,:) :: &
|
||||
relaxationVector
|
||||
end type tRGCstate
|
||||
|
||||
|
||||
type :: tRGCdependentState
|
||||
real(pReal), allocatable, dimension(:) :: &
|
||||
volumeDiscrepancy, &
|
||||
|
|
@ -42,7 +42,24 @@ submodule(homogenization) homogenization_mech_RGC
|
|||
real(pReal), allocatable, dimension(:,:,:) :: &
|
||||
orientation
|
||||
end type tRGCdependentState
|
||||
|
||||
|
||||
type :: tNumerics
|
||||
real(pReal) :: &
|
||||
atol, & !< absolute tolerance of RGC residuum
|
||||
rtol, & !< relative tolerance of RGC residuum
|
||||
absMax, & !< absolute maximum of RGC residuum
|
||||
relMax, & !< relative maximum of RGC residuum
|
||||
pPert, & !< perturbation for computing RGC penalty tangent
|
||||
xSmoo, & !< RGC penalty smoothing parameter (hyperbolic tangent)
|
||||
viscPower, & !< power (sensitivity rate) of numerical viscosity in RGC scheme, Default 1.0e0: Newton viscosity (linear model)
|
||||
viscModus, & !< stress modulus of RGC numerical viscosity, Default 0.0e0: No viscosity is applied
|
||||
refRelaxRate, & !< reference relaxation rate in RGC viscosity
|
||||
maxdRelax, & !< threshold of maximum relaxation vector increment (if exceed this then cutback)
|
||||
maxVolDiscr, & !< threshold of maximum volume discrepancy allowed
|
||||
volDiscrMod, & !< stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint)
|
||||
volDiscrPow !< powerlaw penalty for volume discrepancy
|
||||
end type tNumerics
|
||||
|
||||
type(tparameters), dimension(:), allocatable :: &
|
||||
param
|
||||
type(tRGCstate), dimension(:), allocatable :: &
|
||||
|
|
@ -50,6 +67,8 @@ submodule(homogenization) homogenization_mech_RGC
|
|||
state0
|
||||
type(tRGCdependentState), dimension(:), allocatable :: &
|
||||
dependentState
|
||||
type(tNumerics) :: &
|
||||
num ! numerics parameters. Better name?
|
||||
|
||||
contains
|
||||
|
||||
|
|
@ -63,26 +82,52 @@ module subroutine mech_RGC_init
|
|||
h, &
|
||||
NofMyHomog, &
|
||||
sizeState, nIntFaceTot
|
||||
|
||||
|
||||
|
||||
write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_RGC_label//' init -+>>>'; flush(6)
|
||||
|
||||
|
||||
write(6,'(/,a)') ' Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009'
|
||||
write(6,'(a)') ' https://doi.org/10.1007/s12289-009-0619-1'
|
||||
|
||||
|
||||
write(6,'(/,a)') ' Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010'
|
||||
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) &
|
||||
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
|
||||
|
||||
|
||||
allocate(param(Ninstance))
|
||||
allocate(state(Ninstance))
|
||||
allocate(state0(Ninstance))
|
||||
allocate(dependentState(Ninstance))
|
||||
|
||||
|
||||
|
||||
num%atol = config_numerics%getFloat('atol_rgc', defaultVal=1.0e+4_pReal)
|
||||
num%rtol = config_numerics%getFloat('rtol_rgc', defaultVal=1.0e-3_pReal)
|
||||
num%absMax = config_numerics%getFloat('amax_rgc', defaultVal=1.0e+10_pReal)
|
||||
num%relMax = config_numerics%getFloat('rmax_rgc', defaultVal=1.0e+2_pReal)
|
||||
num%pPert = config_numerics%getFloat('perturbpenalty_rgc', defaultVal=1.0e-7_pReal)
|
||||
num%xSmoo = config_numerics%getFloat('relvantmismatch_rgc', defaultVal=1.0e-5_pReal)
|
||||
num%viscPower = config_numerics%getFloat('viscositypower_rgc', defaultVal=1.0e+0_pReal)
|
||||
num%viscModus = config_numerics%getFloat('viscositymodulus_rgc', defaultVal=0.0e+0_pReal)
|
||||
num%refRelaxRate = config_numerics%getFloat('refrelaxationrate_rgc',defaultVal=1.0e-3_pReal)
|
||||
num%maxdRelax = config_numerics%getFloat('maxrelaxationrate_rgc',defaultVal=1.0e+0_pReal)
|
||||
num%maxVolDiscr = config_numerics%getFloat('maxvoldiscrepancy_rgc',defaultVal=1.0e-5_pReal)
|
||||
num%volDiscrMod = config_numerics%getFloat('voldiscrepancymod_rgc',defaultVal=1.0e+12_pReal)
|
||||
num%volDiscrPow = config_numerics%getFloat('dicrepancypower_rgc', defaultVal=5.0_pReal)
|
||||
|
||||
if (num%atol <= 0.0_pReal) call IO_error(301,ext_msg='absTol_RGC')
|
||||
if (num%rtol <= 0.0_pReal) call IO_error(301,ext_msg='relTol_RGC')
|
||||
if (num%absMax <= 0.0_pReal) call IO_error(301,ext_msg='absMax_RGC')
|
||||
if (num%relMax <= 0.0_pReal) call IO_error(301,ext_msg='relMax_RGC')
|
||||
if (num%pPert <= 0.0_pReal) call IO_error(301,ext_msg='pPert_RGC')
|
||||
if (num%xSmoo <= 0.0_pReal) call IO_error(301,ext_msg='xSmoo_RGC')
|
||||
if (num%viscPower < 0.0_pReal) call IO_error(301,ext_msg='viscPower_RGC')
|
||||
if (num%viscModus < 0.0_pReal) call IO_error(301,ext_msg='viscModus_RGC')
|
||||
if (num%refRelaxRate <= 0.0_pReal) call IO_error(301,ext_msg='refRelaxRate_RGC')
|
||||
if (num%maxdRelax <= 0.0_pReal) call IO_error(301,ext_msg='maxdRelax_RGC')
|
||||
if (num%maxVolDiscr <= 0.0_pReal) call IO_error(301,ext_msg='maxVolDiscr_RGC')
|
||||
if (num%volDiscrMod < 0.0_pReal) call IO_error(301,ext_msg='volDiscrMod_RGC')
|
||||
if (num%volDiscrPow <= 0.0_pReal) call IO_error(301,ext_msg='volDiscrPw_RGC')
|
||||
|
||||
do h = 1, size(homogenization_type)
|
||||
if (homogenization_type(h) /= HOMOGENIZATION_RGC_ID) cycle
|
||||
associate(prm => param(homogenization_typeInstance(h)), &
|
||||
|
|
@ -90,24 +135,24 @@ module subroutine mech_RGC_init
|
|||
st0 => state0(homogenization_typeInstance(h)), &
|
||||
dst => dependentState(homogenization_typeInstance(h)), &
|
||||
config => config_homogenization(h))
|
||||
|
||||
|
||||
#ifdef DEBUG
|
||||
if (h==material_homogenizationAt(debug_e)) then
|
||||
prm%of_debug = material_homogenizationMemberAt(debug_i,debug_e)
|
||||
endif
|
||||
#endif
|
||||
|
||||
prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
|
||||
|
||||
prm%Nconstituents = config%getInts('clustersize',requiredSize=3)
|
||||
if (homogenization_Ngrains(h) /= product(prm%Nconstituents)) &
|
||||
call IO_error(211,ext_msg='clustersize ('//HOMOGENIZATION_RGC_label//')')
|
||||
|
||||
|
||||
prm%xiAlpha = config%getFloat('scalingparameter')
|
||||
prm%ciAlpha = config%getFloat('overproportionality')
|
||||
|
||||
|
||||
prm%dAlpha = config%getFloats('grainsize', requiredSize=3)
|
||||
prm%angles = config%getFloats('clusterorientation',requiredSize=3)
|
||||
|
||||
prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
|
||||
|
||||
NofMyHomog = count(material_homogenizationAt == h)
|
||||
nIntFaceTot = 3*( (prm%Nconstituents(1)-1)*prm%Nconstituents(2)*prm%Nconstituents(3) &
|
||||
|
|
@ -115,17 +160,17 @@ module subroutine mech_RGC_init
|
|||
+ prm%Nconstituents(1)*prm%Nconstituents(2)*(prm%Nconstituents(3)-1))
|
||||
sizeState = nIntFaceTot &
|
||||
+ size(['avg constitutive work ','average penalty energy'])
|
||||
|
||||
|
||||
homogState(h)%sizeState = sizeState
|
||||
allocate(homogState(h)%state0 (sizeState,NofMyHomog), source=0.0_pReal)
|
||||
allocate(homogState(h)%subState0(sizeState,NofMyHomog), source=0.0_pReal)
|
||||
allocate(homogState(h)%state (sizeState,NofMyHomog), source=0.0_pReal)
|
||||
|
||||
|
||||
stt%relaxationVector => homogState(h)%state(1:nIntFaceTot,:)
|
||||
st0%relaxationVector => homogState(h)%state0(1:nIntFaceTot,:)
|
||||
stt%work => homogState(h)%state(nIntFaceTot+1,:)
|
||||
stt%penaltyEnergy => homogState(h)%state(nIntFaceTot+2,:)
|
||||
|
||||
|
||||
allocate(dst%volumeDiscrepancy( NofMyHomog))
|
||||
allocate(dst%relaxationRate_avg( NofMyHomog))
|
||||
allocate(dst%relaxationRate_max( NofMyHomog))
|
||||
|
|
@ -135,11 +180,11 @@ module subroutine mech_RGC_init
|
|||
! assigning cluster orientations
|
||||
dependentState(homogenization_typeInstance(h))%orientation = spread(eu2om(prm%angles*inRad),3,NofMyHomog)
|
||||
!dst%orientation = spread(eu2om(prm%angles*inRad),3,NofMyHomog) ifort version 18.0.1 crashes (for whatever reason)
|
||||
|
||||
|
||||
end associate
|
||||
|
||||
enddo
|
||||
|
||||
|
||||
enddo
|
||||
|
||||
end subroutine mech_RGC_init
|
||||
|
||||
|
||||
|
|
@ -149,19 +194,19 @@ end subroutine mech_RGC_init
|
|||
module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
|
||||
|
||||
real(pReal), dimension (:,:,:), intent(out) :: F !< partioned F per grain
|
||||
|
||||
|
||||
real(pReal), dimension (3,3), intent(in) :: avgF !< averaged F
|
||||
integer, intent(in) :: &
|
||||
instance, &
|
||||
of
|
||||
|
||||
|
||||
real(pReal), dimension(3) :: aVect,nVect
|
||||
integer, dimension(4) :: intFace
|
||||
integer, dimension(3) :: iGrain3
|
||||
integer :: iGrain,iFace,i,j
|
||||
|
||||
|
||||
associate(prm => param(instance))
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! compute the deformation gradient of individual grains due to relaxations
|
||||
F = 0.0_pReal
|
||||
|
|
@ -187,14 +232,14 @@ module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
|
|||
endif
|
||||
#endif
|
||||
enddo
|
||||
|
||||
|
||||
end associate
|
||||
|
||||
end subroutine mech_RGC_partitionDeformation
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief update the internal state of the homogenization scheme and tell whether "done" and
|
||||
!> @brief update the internal state of the homogenization scheme and tell whether "done" and
|
||||
! "happy" with result
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module procedure mech_RGC_updateState
|
||||
|
|
@ -213,17 +258,17 @@ module procedure mech_RGC_updateState
|
|||
integer, dimension(3) :: stresLoc
|
||||
integer, dimension(2) :: residLoc
|
||||
#endif
|
||||
|
||||
|
||||
zeroTimeStep: if(dEq0(dt)) then
|
||||
mech_RGC_updateState = .true. ! pretend everything is fine and return
|
||||
return
|
||||
return
|
||||
endif zeroTimeStep
|
||||
|
||||
instance = homogenization_typeInstance(material_homogenizationAt(el))
|
||||
of = material_homogenizationMemberAt(ip,el)
|
||||
|
||||
|
||||
associate(stt => state(instance), st0 => state0(instance), dst => dependentState(instance), prm => param(instance))
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! get the dimension of the cluster (grains and interfaces)
|
||||
nGDim = prm%Nconstituents
|
||||
|
|
@ -238,7 +283,7 @@ module procedure mech_RGC_updateState
|
|||
allocate(tract(nIntFaceTot,3), source=0.0_pReal)
|
||||
relax = stt%relaxationVector(:,of)
|
||||
drelax = stt%relaxationVector(:,of) - st0%relaxationVector(:,of)
|
||||
|
||||
|
||||
#ifdef DEBUG
|
||||
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0) then
|
||||
write(6,'(1x,a30)')'Obtained state: '
|
||||
|
|
@ -254,7 +299,7 @@ module procedure mech_RGC_updateState
|
|||
call stressPenalty(R,NN,avgF,F,ip,el,instance,of)
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! calculating volume discrepancy and stress penalty related to overall volume discrepancy
|
||||
! calculating volume discrepancy and stress penalty related to overall volume discrepancy
|
||||
call volumePenalty(D,dst%volumeDiscrepancy(of),avgF,F,nGrain,instance,of)
|
||||
|
||||
#ifdef DEBUG
|
||||
|
|
@ -284,7 +329,7 @@ module procedure mech_RGC_updateState
|
|||
iGrN = grain3to1(iGr3N,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
|
||||
intFaceN = getInterface(2*faceID(1),iGr3N)
|
||||
normN = interfaceNormal(intFaceN,instance,of)
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! identify the right/up/front grain (+|P)
|
||||
iGr3P = iGr3N
|
||||
|
|
@ -296,7 +341,7 @@ module procedure mech_RGC_updateState
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
! compute the residual of traction at the interface (in local system, 4-dimensional index)
|
||||
do i = 1,3
|
||||
tract(iNum,i) = sign(viscModus_RGC*(abs(drelax(i+3*(iNum-1)))/(refRelaxRate_RGC*dt))**viscPower_RGC, &
|
||||
tract(iNum,i) = sign(num%viscModus*(abs(drelax(i+3*(iNum-1)))/(num%refRelaxRate*dt))**num%viscPower, &
|
||||
drelax(i+3*(iNum-1))) ! contribution from the relaxation viscosity
|
||||
do j = 1,3
|
||||
tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP) + D(i,j,iGrP))*normP(j) & ! contribution from material stress P, mismatch penalty R, and volume penalty D projected into the interface
|
||||
|
|
@ -304,7 +349,7 @@ module procedure mech_RGC_updateState
|
|||
resid(i+3*(iNum-1)) = tract(iNum,i) ! translate the local residual into global 1-dimensional residual array
|
||||
enddo
|
||||
enddo
|
||||
|
||||
|
||||
#ifdef DEBUG
|
||||
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0) then
|
||||
write(6,'(1x,a30,1x,i3)')'Traction at interface: ',iNum
|
||||
|
|
@ -332,12 +377,12 @@ module procedure mech_RGC_updateState
|
|||
flush(6)
|
||||
endif
|
||||
#endif
|
||||
|
||||
|
||||
mech_RGC_updateState = .false.
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! If convergence reached => done and happy
|
||||
if (residMax < relTol_RGC*stresMax .or. residMax < absTol_RGC) then
|
||||
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) &
|
||||
|
|
@ -358,7 +403,7 @@ module procedure mech_RGC_updateState
|
|||
dst%mismatch(1:3,of) = sum(NN,2)/real(nGrain,pReal)
|
||||
dst%relaxationRate_avg(of) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal)
|
||||
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
|
||||
write(6,'(1x,a30,1x,e15.8)') 'Constitutive work: ',stt%work(of)
|
||||
|
|
@ -372,31 +417,31 @@ module procedure mech_RGC_updateState
|
|||
flush(6)
|
||||
endif
|
||||
#endif
|
||||
|
||||
|
||||
return
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! if residual blows-up => done but unhappy
|
||||
elseif (residMax > relMax_RGC*stresMax .or. residMax > absMax_RGC) 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
|
||||
|
||||
|
||||
#ifdef DEBUG
|
||||
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 .and. prm%of_debug == of) &
|
||||
write(6,'(1x,a,/)') '... broken'; flush(6)
|
||||
#endif
|
||||
|
||||
return
|
||||
|
||||
|
||||
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) &
|
||||
write(6,'(1x,a,/)') '... not yet done'; flush(6)
|
||||
#endif
|
||||
|
||||
|
||||
endif
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! construct the global Jacobian matrix for updating the global relaxation vector array when
|
||||
! construct the global Jacobian matrix for updating the global relaxation vector array when
|
||||
! convergence is not yet reached ...
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -404,7 +449,7 @@ module procedure mech_RGC_updateState
|
|||
allocate(smatrix(3*nIntFaceTot,3*nIntFaceTot), source=0.0_pReal)
|
||||
do iNum = 1,nIntFaceTot
|
||||
faceID = interface1to4(iNum,param(instance)%Nconstituents) ! assembling of local dPdF into global Jacobian matrix
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! identify the left/bottom/back grain (-|N)
|
||||
iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate sytem
|
||||
|
|
@ -424,7 +469,7 @@ module procedure mech_RGC_updateState
|
|||
! to obtain the Jacobian matrix contribution of dPdF
|
||||
endif
|
||||
enddo
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! identify the right/up/front grain (+|P)
|
||||
iGr3P = iGr3N
|
||||
|
|
@ -444,7 +489,7 @@ module procedure mech_RGC_updateState
|
|||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
||||
|
||||
#ifdef DEBUG
|
||||
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0) then
|
||||
write(6,'(1x,a30)')'Jacobian matrix of stress'
|
||||
|
|
@ -455,9 +500,9 @@ module procedure mech_RGC_updateState
|
|||
flush(6)
|
||||
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"
|
||||
allocate(pmatrix(3*nIntFaceTot,3*nIntFaceTot), source=0.0_pReal)
|
||||
allocate(p_relax(3*nIntFaceTot), source=0.0_pReal)
|
||||
|
|
@ -465,7 +510,7 @@ module procedure mech_RGC_updateState
|
|||
|
||||
do ipert = 1,3*nIntFaceTot
|
||||
p_relax = relax
|
||||
p_relax(ipert) = relax(ipert) + pPert_RGC ! perturb the relaxation vector
|
||||
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
|
||||
|
|
@ -483,7 +528,7 @@ module procedure mech_RGC_updateState
|
|||
iGrN = grain3to1(iGr3N,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
|
||||
intFaceN = getInterface(2*faceID(1),iGr3N) ! identify the interface ID of the grain
|
||||
normN = interfaceNormal(intFaceN,instance,of)
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! identify the right/up/front grain (+|P)
|
||||
iGr3P = iGr3N
|
||||
|
|
@ -491,9 +536,9 @@ module procedure mech_RGC_updateState
|
|||
iGrP = grain3to1(iGr3P,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
|
||||
intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identify the interface ID of the grain
|
||||
normP = interfaceNormal(intFaceP,instance,of)
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! compute the residual stress (contribution of mismatch and volume penalties) from perturbed state
|
||||
! compute the residual stress (contribution of mismatch and volume penalties) from perturbed state
|
||||
! at all interfaces
|
||||
do i = 1,3; do j = 1,3
|
||||
p_resid(i+3*(iNum-1)) = p_resid(i+3*(iNum-1)) + (pR(i,j,iGrP) - R(i,j,iGrP))*normP(j) &
|
||||
|
|
@ -502,9 +547,9 @@ module procedure mech_RGC_updateState
|
|||
+ (pD(i,j,iGrN) - D(i,j,iGrN))*normN(j)
|
||||
enddo; enddo
|
||||
enddo
|
||||
pmatrix(:,ipert) = p_resid/pPert_RGC
|
||||
pmatrix(:,ipert) = p_resid/num%pPert
|
||||
enddo
|
||||
|
||||
|
||||
#ifdef DEBUG
|
||||
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0) then
|
||||
write(6,'(1x,a30)')'Jacobian matrix of penalty'
|
||||
|
|
@ -515,15 +560,15 @@ module procedure mech_RGC_updateState
|
|||
flush(6)
|
||||
endif
|
||||
#endif
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! ... of the numerical viscosity traction "rmatrix"
|
||||
allocate(rmatrix(3*nIntFaceTot,3*nIntFaceTot),source=0.0_pReal)
|
||||
do i=1,3*nIntFaceTot
|
||||
rmatrix(i,i) = viscModus_RGC*viscPower_RGC/(refRelaxRate_RGC*dt)* & ! tangent due to numerical viscosity traction appears
|
||||
(abs(drelax(i))/(refRelaxRate_RGC*dt))**(viscPower_RGC - 1.0_pReal) ! only in the main diagonal term
|
||||
rmatrix(i,i) = num%viscModus*num%viscPower/(num%refRelaxRate*dt)* & ! tangent due to numerical viscosity traction appears
|
||||
(abs(drelax(i))/(num%refRelaxRate*dt))**(num%viscPower - 1.0_pReal) ! only in the main diagonal term
|
||||
enddo
|
||||
|
||||
|
||||
#ifdef DEBUG
|
||||
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0) then
|
||||
write(6,'(1x,a30)')'Jacobian matrix of penalty'
|
||||
|
|
@ -554,7 +599,7 @@ module procedure mech_RGC_updateState
|
|||
! computing the update of the state variable (relaxation vectors) using the Jacobian matrix
|
||||
allocate(jnverse(3*nIntFaceTot,3*nIntFaceTot),source=0.0_pReal)
|
||||
call math_invert(jnverse,error,jmatrix)
|
||||
|
||||
|
||||
#ifdef DEBUG
|
||||
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0) then
|
||||
write(6,'(1x,a30)')'Jacobian inverse'
|
||||
|
|
@ -573,7 +618,7 @@ module procedure mech_RGC_updateState
|
|||
drelax(i) = drelax(i) - jnverse(i,j)*resid(j) ! Calculate the correction for the state variable
|
||||
enddo; enddo
|
||||
stt%relaxationVector(:,of) = relax + drelax ! Updateing the state variable for the next iteration
|
||||
if (any(abs(drelax) > maxdRelax_RGC)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large
|
||||
if (any(abs(drelax) > num%maxdRelax)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large
|
||||
mech_RGC_updateState = [.true.,.false.]
|
||||
!$OMP CRITICAL (write2out)
|
||||
write(6,'(1x,a,1x,i3,1x,a,1x,i3,1x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback'
|
||||
|
|
@ -603,11 +648,11 @@ module procedure mech_RGC_updateState
|
|||
|
||||
real(pReal), dimension (:,:,:), intent(out) :: rPen !< stress-like penalty
|
||||
real(pReal), dimension (:,:), intent(out) :: nMis !< total amount of mismatch
|
||||
|
||||
|
||||
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
|
||||
|
||||
|
||||
integer, dimension (4) :: intFace
|
||||
integer, dimension (3) :: iGrain3,iGNghb3,nGDim
|
||||
real(pReal), dimension (3,3) :: gDef,nDef
|
||||
|
|
@ -623,13 +668,13 @@ module procedure mech_RGC_updateState
|
|||
nGDim = param(instance)%Nconstituents
|
||||
rPen = 0.0_pReal
|
||||
nMis = 0.0_pReal
|
||||
|
||||
|
||||
!----------------------------------------------------------------------------------------------
|
||||
! get the correction factor the modulus of penalty stress representing the evolution of area of
|
||||
! get the correction factor the modulus of penalty stress representing the evolution of area of
|
||||
! the interfaces due to deformations
|
||||
|
||||
surfCorr = surfaceCorrection(avgF,instance,of)
|
||||
|
||||
|
||||
associate(prm => param(instance))
|
||||
|
||||
#ifdef DEBUG
|
||||
|
|
@ -640,15 +685,15 @@ module procedure mech_RGC_updateState
|
|||
write(6,*) surfCorr
|
||||
endif
|
||||
#endif
|
||||
|
||||
|
||||
!-----------------------------------------------------------------------------------------------
|
||||
! computing the mismatch and penalty stress tensor of all grains
|
||||
! computing the mismatch and penalty stress tensor of all grains
|
||||
grainLoop: do iGrain = 1,product(prm%Nconstituents)
|
||||
Gmoduli = equivalentModuli(iGrain,ip,el)
|
||||
muGrain = Gmoduli(1) ! collecting the equivalent shear modulus of grain
|
||||
bgGrain = Gmoduli(2) ! and the lengthh of Burgers vector
|
||||
iGrain3 = grain1to3(iGrain,prm%Nconstituents) ! get the grain ID in local 3-dimensional index (x,y,z)-position
|
||||
|
||||
|
||||
interfaceLoop: do iFace = 1,6
|
||||
intFace = getInterface(iFace,iGrain3) ! get the 4-dimensional index of the interface in local numbering system of the grain
|
||||
nVect = interfaceNormal(intFace,instance,of)
|
||||
|
|
@ -662,7 +707,7 @@ module procedure mech_RGC_updateState
|
|||
muGNghb = Gmoduli(1)
|
||||
bgGNghb = Gmoduli(2)
|
||||
gDef = 0.5_pReal*(fDef(1:3,1:3,iGNghb) - fDef(1:3,1:3,iGrain)) ! difference/jump in deformation gradeint across the neighbor
|
||||
|
||||
|
||||
!-------------------------------------------------------------------------------------------
|
||||
! compute the mismatch tensor of all interfaces
|
||||
nDefNorm = 0.0_pReal
|
||||
|
|
@ -690,7 +735,7 @@ module procedure mech_RGC_updateState
|
|||
*surfCorr(abs(intFace(1)))/prm%dAlpha(abs(intFace(1))) &
|
||||
*cosh(prm%ciAlpha*nDefNorm) &
|
||||
*0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) &
|
||||
*tanh(nDefNorm/xSmoo_RGC)
|
||||
*tanh(nDefNorm/num%xSmoo)
|
||||
enddo; enddo;enddo; enddo
|
||||
enddo interfaceLoop
|
||||
#ifdef DEBUG
|
||||
|
|
@ -699,32 +744,32 @@ module procedure mech_RGC_updateState
|
|||
write(6,*) transpose(rPen(1:3,1:3,iGrain))
|
||||
endif
|
||||
#endif
|
||||
|
||||
|
||||
enddo grainLoop
|
||||
|
||||
|
||||
end associate
|
||||
|
||||
|
||||
end subroutine stressPenalty
|
||||
|
||||
|
||||
|
||||
|
||||
!------------------------------------------------------------------------------------------------
|
||||
!> @brief calculate stress-like penalty due to volume discrepancy
|
||||
!> @brief calculate stress-like penalty due to volume discrepancy
|
||||
!------------------------------------------------------------------------------------------------
|
||||
subroutine volumePenalty(vPen,vDiscrep,fAvg,fDef,nGrain,instance,of)
|
||||
|
||||
|
||||
real(pReal), dimension (:,:,:), intent(out) :: vPen ! stress-like penalty due to volume
|
||||
real(pReal), intent(out) :: vDiscrep ! total volume discrepancy
|
||||
|
||||
|
||||
real(pReal), dimension (:,:,:), intent(in) :: fDef ! deformation gradients
|
||||
real(pReal), dimension (3,3), intent(in) :: fAvg ! overall deformation gradient
|
||||
integer, intent(in) :: &
|
||||
Ngrain, &
|
||||
instance, &
|
||||
of
|
||||
|
||||
|
||||
real(pReal), dimension(size(vPen,3)) :: gVol
|
||||
integer :: i
|
||||
|
||||
|
||||
!----------------------------------------------------------------------------------------------
|
||||
! compute the volumes of grains and of cluster
|
||||
vDiscrep = math_det33(fAvg) ! compute the volume of the cluster
|
||||
|
|
@ -733,15 +778,15 @@ module procedure mech_RGC_updateState
|
|||
vDiscrep = vDiscrep - gVol(i)/real(nGrain,pReal) ! calculate the difference/dicrepancy between
|
||||
! the volume of the cluster and the the total volume of grains
|
||||
enddo
|
||||
|
||||
|
||||
!----------------------------------------------------------------------------------------------
|
||||
! calculate the stress and penalty due to volume discrepancy
|
||||
vPen = 0.0_pReal
|
||||
do i = 1,nGrain
|
||||
vPen(:,:,i) = -1.0_pReal/real(nGrain,pReal)*volDiscrMod_RGC*volDiscrPow_RGC/maxVolDiscr_RGC* &
|
||||
sign((abs(vDiscrep)/maxVolDiscr_RGC)**(volDiscrPow_RGC - 1.0),vDiscrep)* &
|
||||
vPen(:,:,i) = -1.0_pReal/real(nGrain,pReal)*num%volDiscrMod*num%volDiscrPow/num%maxVolDiscr* &
|
||||
sign((abs(vDiscrep)/num%maxVolDiscr)**(num%volDiscrPow - 1.0),vDiscrep)* &
|
||||
gVol(i)*transpose(math_inv33(fDef(:,:,i)))
|
||||
|
||||
|
||||
#ifdef DEBUG
|
||||
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 &
|
||||
.and. param(instance)%of_debug == of) then
|
||||
|
|
@ -755,11 +800,11 @@ module procedure mech_RGC_updateState
|
|||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief compute the correction factor accouted for surface evolution (area change) due to
|
||||
!> @brief compute the correction factor accouted for surface evolution (area change) due to
|
||||
! deformation
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function surfaceCorrection(avgF,instance,of)
|
||||
|
||||
|
||||
real(pReal), dimension(3) :: surfaceCorrection
|
||||
|
||||
real(pReal), dimension(3,3), intent(in) :: avgF !< average F
|
||||
|
|
@ -771,9 +816,9 @@ module procedure mech_RGC_updateState
|
|||
real(pReal) :: detF
|
||||
integer :: i,j,iBase
|
||||
logical :: error
|
||||
|
||||
|
||||
call math_invert33(invC,detF,error,matmul(transpose(avgF),avgF))
|
||||
|
||||
|
||||
surfaceCorrection = 0.0_pReal
|
||||
do iBase = 1,3
|
||||
nVect = interfaceNormal([iBase,1,1,1],instance,of)
|
||||
|
|
@ -782,7 +827,7 @@ module procedure mech_RGC_updateState
|
|||
enddo; enddo
|
||||
surfaceCorrection(iBase) = sqrt(surfaceCorrection(iBase))*detF ! get the surface correction factor (area contraction/enlargement)
|
||||
enddo
|
||||
|
||||
|
||||
end function surfaceCorrection
|
||||
|
||||
|
||||
|
|
@ -790,9 +835,9 @@ module procedure mech_RGC_updateState
|
|||
!> @brief compute the equivalent shear and bulk moduli from the elasticity tensor
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function equivalentModuli(grainID,ip,el)
|
||||
|
||||
|
||||
real(pReal), dimension(2) :: equivalentModuli
|
||||
|
||||
|
||||
integer, intent(in) :: &
|
||||
grainID,&
|
||||
ip, & !< integration point number
|
||||
|
|
@ -802,9 +847,9 @@ module procedure mech_RGC_updateState
|
|||
cEquiv_11, &
|
||||
cEquiv_12, &
|
||||
cEquiv_44
|
||||
|
||||
|
||||
elasTens = constitutive_homogenizedC(grainID,ip,el)
|
||||
|
||||
|
||||
!----------------------------------------------------------------------------------------------
|
||||
! compute the equivalent shear modulus after Turterltaub and Suiker, JMPS (2005)
|
||||
cEquiv_11 = (elasTens(1,1) + elasTens(2,2) + elasTens(3,3))/3.0_pReal
|
||||
|
|
@ -812,37 +857,37 @@ module procedure mech_RGC_updateState
|
|||
elasTens(1,3) + elasTens(2,1) + elasTens(3,2))/6.0_pReal
|
||||
cEquiv_44 = (elasTens(4,4) + elasTens(5,5) + elasTens(6,6))/3.0_pReal
|
||||
equivalentModuli(1) = 0.2_pReal*(cEquiv_11 - cEquiv_12) + 0.6_pReal*cEquiv_44
|
||||
|
||||
|
||||
!----------------------------------------------------------------------------------------------
|
||||
! obtain the length of Burgers vector (could be model dependend)
|
||||
equivalentModuli(2) = 2.5e-10_pReal
|
||||
|
||||
|
||||
end function equivalentModuli
|
||||
|
||||
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief calculating the grain deformation gradient (the same with
|
||||
!> @brief calculating the grain deformation gradient (the same with
|
||||
! homogenization_RGC_partitionDeformation, but used only for perturbation scheme)
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine grainDeformation(F, avgF, instance, of)
|
||||
|
||||
|
||||
real(pReal), dimension(:,:,:), intent(out) :: F !< partioned F per grain
|
||||
|
||||
|
||||
real(pReal), dimension(:,:), intent(in) :: avgF !< averaged F
|
||||
integer, intent(in) :: &
|
||||
instance, &
|
||||
of
|
||||
|
||||
|
||||
real(pReal), dimension(3) :: aVect,nVect
|
||||
integer, dimension(4) :: intFace
|
||||
integer, dimension(3) :: iGrain3
|
||||
integer :: iGrain,iFace,i,j
|
||||
|
||||
|
||||
!-------------------------------------------------------------------------------------------------
|
||||
! compute the deformation gradient of individual grains due to relaxations
|
||||
|
||||
|
||||
associate(prm => param(instance))
|
||||
|
||||
|
||||
F = 0.0_pReal
|
||||
do iGrain = 1,product(prm%Nconstituents)
|
||||
iGrain3 = grain1to3(iGrain,prm%Nconstituents)
|
||||
|
|
@ -855,16 +900,16 @@ module procedure mech_RGC_updateState
|
|||
enddo
|
||||
F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! relaxed deformation gradient
|
||||
enddo
|
||||
|
||||
|
||||
end associate
|
||||
|
||||
|
||||
end subroutine grainDeformation
|
||||
|
||||
|
||||
end procedure mech_RGC_updateState
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief derive average stress and stiffness from constituent quantities
|
||||
!> @brief derive average stress and stiffness from constituent quantities
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
|
||||
|
||||
|
|
@ -888,9 +933,9 @@ module subroutine mech_RGC_results(instance,group)
|
|||
|
||||
integer, intent(in) :: instance
|
||||
character(len=*), intent(in) :: group
|
||||
|
||||
|
||||
integer :: o
|
||||
|
||||
|
||||
associate(stt => state(instance), dst => dependentState(instance), prm => param(instance))
|
||||
outputsLoop: do o = 1,size(prm%output)
|
||||
select case(trim(prm%output(o)))
|
||||
|
|
@ -915,7 +960,7 @@ module subroutine mech_RGC_results(instance,group)
|
|||
end select
|
||||
enddo outputsLoop
|
||||
end associate
|
||||
|
||||
|
||||
end subroutine mech_RGC_results
|
||||
|
||||
|
||||
|
|
@ -928,7 +973,7 @@ pure function relaxationVector(intFace,instance,of)
|
|||
|
||||
integer, intent(in) :: instance,of
|
||||
integer, dimension(4), intent(in) :: intFace !< set of interface ID in 4D array (normal and position)
|
||||
|
||||
|
||||
integer :: iNum
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -945,10 +990,10 @@ end function relaxationVector
|
|||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief identify the normal of an interface
|
||||
!> @brief identify the normal of an interface
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure function interfaceNormal(intFace,instance,of)
|
||||
|
||||
|
||||
real(pReal), dimension(3) :: interfaceNormal
|
||||
|
||||
integer, dimension(4), intent(in) :: intFace !< interface ID in 4D array (normal and position)
|
||||
|
|
@ -980,10 +1025,10 @@ pure function getInterface(iFace,iGrain3)
|
|||
integer, intent(in) :: iFace !< face index (1..6) mapped like (-e1,-e2,-e3,+e1,+e2,+e3) or iDir = (-1,-2,-3,1,2,3)
|
||||
|
||||
integer :: iDir !< direction of interface normal
|
||||
|
||||
|
||||
iDir = (int(real(iFace-1,pReal)/2.0_pReal)+1)*(-1)**iFace
|
||||
getInterface(1) = iDir
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! identify the interface position by the direction of its normal
|
||||
getInterface(2:4) = iGrain3
|
||||
|
|
@ -996,7 +1041,7 @@ end function getInterface
|
|||
!> @brief map grain ID from in 1D (global array) to in 3D (local position)
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure function grain1to3(grain1,nGDim)
|
||||
|
||||
|
||||
integer, dimension(3) :: grain1to3
|
||||
|
||||
integer, intent(in) :: grain1 !< grain ID in 1D array
|
||||
|
|
@ -1016,7 +1061,7 @@ integer pure function grain3to1(grain3,nGDim)
|
|||
|
||||
integer, dimension(3), intent(in) :: grain3 !< grain ID in 3D array (pos.x,pos.y,pos.z)
|
||||
integer, dimension(3), intent(in) :: nGDim
|
||||
|
||||
|
||||
grain3to1 = grain3(1) &
|
||||
+ nGDim(1)*(grain3(2)-1) &
|
||||
+ nGDim(1)*nGDim(2)*(grain3(3)-1)
|
||||
|
|
@ -1028,11 +1073,11 @@ end function grain3to1
|
|||
!> @brief maps interface ID from 4D (normal and local position) into 1D (global array)
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
integer pure function interface4to1(iFace4D, nGDim)
|
||||
|
||||
|
||||
integer, dimension(4), intent(in) :: iFace4D !< interface ID in 4D array (n.dir,pos.x,pos.y,pos.z)
|
||||
integer, dimension(3), intent(in) :: nGDim
|
||||
|
||||
|
||||
|
||||
select case(abs(iFace4D(1)))
|
||||
|
||||
case(1)
|
||||
|
|
@ -1064,7 +1109,7 @@ integer pure function interface4to1(iFace4D, nGDim)
|
|||
|
||||
case default
|
||||
interface4to1 = -1
|
||||
|
||||
|
||||
end select
|
||||
|
||||
end function interface4to1
|
||||
|
|
@ -1074,7 +1119,7 @@ end function interface4to1
|
|||
!> @brief maps interface ID from 1D (global array) into 4D (normal and local position)
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure function interface1to4(iFace1D, nGDim)
|
||||
|
||||
|
||||
integer, dimension(4) :: interface1to4
|
||||
|
||||
integer, intent(in) :: iFace1D !< interface ID in 1D array
|
||||
|
|
|
|||
109
src/numerics.f90
109
src/numerics.f90
|
|
@ -17,12 +17,12 @@ module numerics
|
|||
private
|
||||
|
||||
integer, protected, public :: &
|
||||
iJacoStiffness = 1, & !< frequency of stiffness update
|
||||
nMPstate = 10, & !< materialpoint state loop limit
|
||||
randomSeed = 0, & !< fixed seeding for pseudo-random number generator, Default 0: use random seed
|
||||
worldrank = 0, & !< MPI worldrank (/=0 for MPI simulations only)
|
||||
worldsize = 1, & !< MPI worldsize (/=1 for MPI simulations only)
|
||||
numerics_integrator = 1 !< method used for state integration Default 1: fix-point iteration
|
||||
iJacoStiffness = 1, & !< frequency of stiffness update
|
||||
nMPstate = 10, & !< materialpoint state loop limit
|
||||
randomSeed = 0, & !< fixed seeding for pseudo-random number generator, Default 0: use random seed
|
||||
worldrank = 0, & !< MPI worldrank (/=0 for MPI simulations only)
|
||||
worldsize = 1, & !< MPI worldsize (/=1 for MPI simulations only)
|
||||
numerics_integrator = 1 !< method used for state integration Default 1: fix-point iteration
|
||||
integer(4), protected, public :: &
|
||||
DAMASK_NumThreadsInt = 0 !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive
|
||||
real(pReal), protected, public :: &
|
||||
|
|
@ -31,22 +31,9 @@ module numerics
|
|||
subStepSizeHomog = 0.25_pReal, & !< size of first substep when cutback in homogenization
|
||||
stepIncreaseHomog = 1.5_pReal, & !< increase of next substep size when previous substep converged in homogenization
|
||||
numerics_unitlength = 1.0_pReal, & !< determines the physical length of one computational length unit
|
||||
absTol_RGC = 1.0e+4_pReal, & !< absolute tolerance of RGC residuum
|
||||
relTol_RGC = 1.0e-3_pReal, & !< relative tolerance of RGC residuum
|
||||
absMax_RGC = 1.0e+10_pReal, & !< absolute maximum of RGC residuum
|
||||
relMax_RGC = 1.0e+2_pReal, & !< relative maximum of RGC residuum
|
||||
pPert_RGC = 1.0e-7_pReal, & !< perturbation for computing RGC penalty tangent
|
||||
xSmoo_RGC = 1.0e-5_pReal, & !< RGC penalty smoothing parameter (hyperbolic tangent)
|
||||
viscPower_RGC = 1.0e+0_pReal, & !< power (sensitivity rate) of numerical viscosity in RGC scheme, Default 1.0e0: Newton viscosity (linear model)
|
||||
viscModus_RGC = 0.0e+0_pReal, & !< stress modulus of RGC numerical viscosity, Default 0.0e0: No viscosity is applied
|
||||
refRelaxRate_RGC = 1.0e-3_pReal, & !< reference relaxation rate in RGC viscosity
|
||||
maxdRelax_RGC = 1.0e+0_pReal, & !< threshold of maximum relaxation vector increment (if exceed this then cutback)
|
||||
maxVolDiscr_RGC = 1.0e-5_pReal, & !< threshold of maximum volume discrepancy allowed
|
||||
volDiscrMod_RGC = 1.0e+12_pReal, & !< stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint)
|
||||
volDiscrPow_RGC = 5.0_pReal, & !< powerlaw penalty for volume discrepancy
|
||||
charLength = 1.0_pReal, & !< characteristic length scale for gradient problems
|
||||
residualStiffness = 1.0e-6_pReal !< non-zero residual damage
|
||||
logical, protected, public :: &
|
||||
residualStiffness = 1.0e-6_pReal !< non-zero residual damage
|
||||
logical, protected, public :: &
|
||||
usePingPong = .true.
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -75,8 +62,8 @@ module numerics
|
|||
err_stress_tolAbs = 1.0e3_pReal, & !< absolute tolerance for fullfillment of stress BC
|
||||
err_stress_tolRel = 0.01_pReal, & !< relative tolerance for fullfillment of stress BC
|
||||
rotation_tol = 1.0e-12_pReal, & !< tolerance of rotation specified in loadcase, Default 1.0e-12: first guess
|
||||
polarAlpha = 1.0_pReal, & !< polarization scheme parameter 0.0 < alpha < 2.0. alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme
|
||||
polarBeta = 1.0_pReal !< polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme
|
||||
polarAlpha = 1.0_pReal, & !< polarization scheme parameter 0.0 < alpha < 2.0. alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme
|
||||
polarBeta = 1.0_pReal !< polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme
|
||||
character(len=pStringLen), protected, public :: &
|
||||
petsc_options = ''
|
||||
#endif
|
||||
|
|
@ -87,8 +74,8 @@ module numerics
|
|||
integer, protected, public :: &
|
||||
integrationOrder = 2, & !< order of quadrature rule required
|
||||
structOrder = 2 !< order of displacement shape functions
|
||||
logical, protected, public :: &
|
||||
BBarStabilisation = .false.
|
||||
logical, protected, public :: &
|
||||
BBarStabilisation = .false.
|
||||
character(len=*), parameter, public :: &
|
||||
petsc_defaultOptions = '-mech_snes_type newtonls &
|
||||
&-mech_snes_linesearch_type cp &
|
||||
|
|
@ -106,7 +93,7 @@ module numerics
|
|||
#endif
|
||||
|
||||
public :: numerics_init
|
||||
|
||||
|
||||
contains
|
||||
|
||||
|
||||
|
|
@ -130,7 +117,7 @@ subroutine numerics_init
|
|||
call MPI_Comm_size(PETSC_COMM_WORLD,worldsize,ierr);CHKERRQ(ierr)
|
||||
#endif
|
||||
write(6,'(/,a)') ' <<<+- numerics init -+>>>'
|
||||
|
||||
|
||||
!$ call GET_ENVIRONMENT_VARIABLE(NAME='DAMASK_NUM_THREADS',VALUE=DAMASK_NumThreadsString,STATUS=gotDAMASK_NUM_THREADS) ! get environment variable DAMASK_NUM_THREADS...
|
||||
!$ if(gotDAMASK_NUM_THREADS /= 0) then ! could not get number of threads, set it to 1
|
||||
!$ call IO_warning(35,ext_msg='BEGIN:'//DAMASK_NumThreadsString//':END')
|
||||
|
|
@ -140,15 +127,15 @@ subroutine numerics_init
|
|||
!$ if (DAMASK_NumThreadsInt < 1_4) DAMASK_NumThreadsInt = 1_4 ! in case of string conversion fails, set it to one
|
||||
!$ endif
|
||||
!$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution
|
||||
|
||||
|
||||
inquire(file='numerics.config', exist=fexist)
|
||||
|
||||
|
||||
fileExists: if (fexist) then
|
||||
write(6,'(a,/)') ' using values from config file'
|
||||
flush(6)
|
||||
fileContent = IO_read_ASCII('numerics.config')
|
||||
fileContent = IO_read_ASCII('numerics.config')
|
||||
do j=1, size(fileContent)
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! read variables from config file and overwrite default parameters if keyword is present
|
||||
line = fileContent(j)
|
||||
|
|
@ -179,35 +166,6 @@ subroutine numerics_init
|
|||
case ('unitlength')
|
||||
numerics_unitlength = IO_floatValue(line,chunkPos,2)
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! RGC parameters
|
||||
case ('atol_rgc')
|
||||
absTol_RGC = IO_floatValue(line,chunkPos,2)
|
||||
case ('rtol_rgc')
|
||||
relTol_RGC = IO_floatValue(line,chunkPos,2)
|
||||
case ('amax_rgc')
|
||||
absMax_RGC = IO_floatValue(line,chunkPos,2)
|
||||
case ('rmax_rgc')
|
||||
relMax_RGC = IO_floatValue(line,chunkPos,2)
|
||||
case ('perturbpenalty_rgc')
|
||||
pPert_RGC = IO_floatValue(line,chunkPos,2)
|
||||
case ('relevantmismatch_rgc')
|
||||
xSmoo_RGC = IO_floatValue(line,chunkPos,2)
|
||||
case ('viscositypower_rgc')
|
||||
viscPower_RGC = IO_floatValue(line,chunkPos,2)
|
||||
case ('viscositymodulus_rgc')
|
||||
viscModus_RGC = IO_floatValue(line,chunkPos,2)
|
||||
case ('refrelaxationrate_rgc')
|
||||
refRelaxRate_RGC = IO_floatValue(line,chunkPos,2)
|
||||
case ('maxrelaxation_rgc')
|
||||
maxdRelax_RGC = IO_floatValue(line,chunkPos,2)
|
||||
case ('maxvoldiscrepancy_rgc')
|
||||
maxVolDiscr_RGC = IO_floatValue(line,chunkPos,2)
|
||||
case ('voldiscrepancymod_rgc')
|
||||
volDiscrMod_RGC = IO_floatValue(line,chunkPos,2)
|
||||
case ('discrepancypower_rgc')
|
||||
volDiscrPow_RGC = IO_floatValue(line,chunkPos,2)
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! random seeding parameter
|
||||
case ('random_seed','fixed_seed')
|
||||
|
|
@ -280,8 +238,6 @@ subroutine numerics_init
|
|||
#endif
|
||||
end select
|
||||
enddo
|
||||
|
||||
|
||||
else fileExists
|
||||
write(6,'(a,/)') ' using standard values'
|
||||
flush(6)
|
||||
|
|
@ -301,21 +257,6 @@ subroutine numerics_init
|
|||
write(6,'(a24,1x,es8.1)') ' stepIncreaseHomog: ',stepIncreaseHomog
|
||||
write(6,'(a24,1x,i8,/)') ' nMPstate: ',nMPstate
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! RGC parameters
|
||||
write(6,'(a24,1x,es8.1)') ' aTol_RGC: ',absTol_RGC
|
||||
write(6,'(a24,1x,es8.1)') ' rTol_RGC: ',relTol_RGC
|
||||
write(6,'(a24,1x,es8.1)') ' aMax_RGC: ',absMax_RGC
|
||||
write(6,'(a24,1x,es8.1)') ' rMax_RGC: ',relMax_RGC
|
||||
write(6,'(a24,1x,es8.1)') ' perturbPenalty_RGC: ',pPert_RGC
|
||||
write(6,'(a24,1x,es8.1)') ' relevantMismatch_RGC: ',xSmoo_RGC
|
||||
write(6,'(a24,1x,es8.1)') ' viscosityrate_RGC: ',viscPower_RGC
|
||||
write(6,'(a24,1x,es8.1)') ' viscositymodulus_RGC: ',viscModus_RGC
|
||||
write(6,'(a24,1x,es8.1)') ' maxrelaxation_RGC: ',maxdRelax_RGC
|
||||
write(6,'(a24,1x,es8.1)') ' maxVolDiscrepancy_RGC: ',maxVolDiscr_RGC
|
||||
write(6,'(a24,1x,es8.1)') ' volDiscrepancyMod_RGC: ',volDiscrMod_RGC
|
||||
write(6,'(a24,1x,es8.1,/)') ' discrepancyPower_RGC: ',volDiscrPow_RGC
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! Random seeding parameter
|
||||
write(6,'(a16,1x,i16,/)') ' random_seed: ',randomSeed
|
||||
|
|
@ -366,7 +307,6 @@ subroutine numerics_init
|
|||
write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_defaultOptions)//' '//trim(petsc_options)
|
||||
write(6,'(a24,1x,L8)') ' B-Bar stabilisation: ',BBarStabilisation
|
||||
#endif
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! sanity checks
|
||||
|
|
@ -379,19 +319,6 @@ subroutine numerics_init
|
|||
if (numerics_integrator <= 0 .or. numerics_integrator >= 6) &
|
||||
call IO_error(301,ext_msg='integrator')
|
||||
if (numerics_unitlength <= 0.0_pReal) call IO_error(301,ext_msg='unitlength')
|
||||
if (absTol_RGC <= 0.0_pReal) call IO_error(301,ext_msg='absTol_RGC')
|
||||
if (relTol_RGC <= 0.0_pReal) call IO_error(301,ext_msg='relTol_RGC')
|
||||
if (absMax_RGC <= 0.0_pReal) call IO_error(301,ext_msg='absMax_RGC')
|
||||
if (relMax_RGC <= 0.0_pReal) call IO_error(301,ext_msg='relMax_RGC')
|
||||
if (pPert_RGC <= 0.0_pReal) call IO_error(301,ext_msg='pPert_RGC')
|
||||
if (xSmoo_RGC <= 0.0_pReal) call IO_error(301,ext_msg='xSmoo_RGC')
|
||||
if (viscPower_RGC < 0.0_pReal) call IO_error(301,ext_msg='viscPower_RGC')
|
||||
if (viscModus_RGC < 0.0_pReal) call IO_error(301,ext_msg='viscModus_RGC')
|
||||
if (refRelaxRate_RGC <= 0.0_pReal) call IO_error(301,ext_msg='refRelaxRate_RGC')
|
||||
if (maxdRelax_RGC <= 0.0_pReal) call IO_error(301,ext_msg='maxdRelax_RGC')
|
||||
if (maxVolDiscr_RGC <= 0.0_pReal) call IO_error(301,ext_msg='maxVolDiscr_RGC')
|
||||
if (volDiscrMod_RGC < 0.0_pReal) call IO_error(301,ext_msg='volDiscrMod_RGC')
|
||||
if (volDiscrPow_RGC <= 0.0_pReal) call IO_error(301,ext_msg='volDiscrPw_RGC')
|
||||
if (residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness')
|
||||
if (itmax <= 1) call IO_error(301,ext_msg='itmax')
|
||||
if (itmin > itmax .or. itmin < 1) call IO_error(301,ext_msg='itmin')
|
||||
|
|
|
|||
Loading…
Reference in New Issue