From e4792e56fb214f1c9866d4ca339515603f5bddcf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 16 Mar 2020 21:39:53 +0100 Subject: [PATCH] store data where it is needed avoid globals, even if they are read only --- src/homogenization_mech_RGC.f90 | 299 ++++++++++++++++++-------------- src/numerics.f90 | 109 ++---------- 2 files changed, 190 insertions(+), 218 deletions(-) diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index 2152211b7..ca84789b4 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -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 diff --git a/src/numerics.f90 b/src/numerics.f90 index 61297ff6f..0f4bd68e8 100644 --- a/src/numerics.f90 +++ b/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')