From 1f4d7c2ca46facee494a58e391387c8369c2d0bb Mon Sep 17 00:00:00 2001 From: Pratheek Shanthraj Date: Wed, 6 Mar 2013 14:31:13 +0000 Subject: [PATCH] changed reference stiffness to min max avg from volume avg. using terminally ill information in convergence check --- code/DAMASK_spectral_solverAL.f90 | 32 +++++---- code/DAMASK_spectral_solverBasic.f90 | 10 +-- code/DAMASK_spectral_solverBasicPETSc.f90 | 15 ++-- code/DAMASK_spectral_utilities.f90 | 86 ++++++++++++++++++++++- 4 files changed, 117 insertions(+), 26 deletions(-) diff --git a/code/DAMASK_spectral_solverAL.f90 b/code/DAMASK_spectral_solverAL.f90 index 377fb2884..f8c0feb35 100644 --- a/code/DAMASK_spectral_solverAL.f90 +++ b/code/DAMASK_spectral_solverAL.f90 @@ -57,7 +57,7 @@ module DAMASK_spectral_solverAL P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress character(len=1024), private :: incInfo !< time and increment information real(pReal), private, dimension(3,3,3,3) :: & - C = 0.0_pReal, & !< current average stiffness + C = 0.0_pReal, C_minmaxAvg = 0.0_pReal, & !< current average stiffness C_lastInc = 0.0_pReal, & !< previous average stiffness S = 0.0_pReal, & !< current compliance (filled up with zeros) C_scale = 0.0_pReal, & @@ -186,25 +186,25 @@ subroutine AL_init(temperature) read (777,rec=1) C_lastInc close (777) call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(temp3333_Real)) - read (777,rec=1) temp3333_Real + read (777,rec=1) C_minmaxAvg close (777) endif mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomdim,reshape(& F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems]) - call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,temp3333_Real2,& + call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,temp3333_Real,temp3333_Real2,& temp33_Real,.false.,math_I3) call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! reference stiffness if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness - temp3333_Real = temp3333_Real2 - C = temp3333_Real2 + C_minmaxAvg = temp3333_Real2 + C = temp3333_Real endif - call Utilities_updateGamma(temp3333_Real,.True.) - C_scale = temp3333_Real - S_scale = math_invSym3333(temp3333_Real) + call Utilities_updateGamma(temp3333_Real2,.True.) + C_scale = temp3333_Real2 + S_scale = math_invSym3333(temp3333_Real2) end subroutine AL_init @@ -219,7 +219,8 @@ type(tSolutionState) function & itmax use math, only: & math_mul33x33 ,& - math_rotate_backward33 + math_rotate_backward33, & + math_invSym3333 use mesh, only: & res, & geomdim, & @@ -342,7 +343,11 @@ use mesh, only: & !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C) - if (update_gamma) call Utilities_updateGamma(C,restartWrite) + if (update_gamma) then + call Utilities_updateGamma(C_minmaxAvg,restartWrite) + C_scale = C_minmaxAvg + S_scale = math_invSym3333(C_minmaxAvg) + endif ForwardData = .True. mask_stress = P_BC%maskFloat @@ -478,7 +483,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response call Utilities_constitutiveResponse(F_lastInc,F - residual_F_tau/polarBeta,params%temperature,params%timeinc, & - residual_F,C,P_av,ForwardData,params%rotation_BC) + residual_F,C,C_minmaxAvg,P_av,ForwardData,params%rotation_BC) ForwardData = .False. !-------------------------------------------------------------------------------------------------- @@ -520,6 +525,8 @@ subroutine AL_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ierr) err_p_tol, & err_stress_tolrel, & err_stress_tolabs + use FEsolving, only: & + terminallyIll implicit none SNES :: snes_local @@ -538,7 +545,8 @@ subroutine AL_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ierr) Converged = (it > itmin .and. & all([ err_f/sqrt(sum((F_aim-math_I3)**2.0_pReal))/err_f_tol, & err_p/sqrt(sum((F_aim-math_I3)**2.0_pReal))/err_p_tol, & - err_stress/err_stress_tol] < 1.0_pReal)) + err_stress/err_stress_tol] < 1.0_pReal)) & + .or. terminallyIll if (Converged) then reason = 1 diff --git a/code/DAMASK_spectral_solverBasic.f90 b/code/DAMASK_spectral_solverBasic.f90 index dbfb26118..6078e383d 100644 --- a/code/DAMASK_spectral_solverBasic.f90 +++ b/code/DAMASK_spectral_solverBasic.f90 @@ -35,7 +35,7 @@ module DAMASK_spectral_SolverBasic F_aim_lastInc = math_I3, & !< deformation gradient aim last increment F_aimDot = 0.0_pReal !< assumed rate real(pReal), private,dimension(3,3,3,3) :: & - C = 0.0_pReal, & !< average stiffness + C = 0.0_pReal, C_minmaxAvg = 0.0_pReal, & !< average stiffness C_lastInc = 0.0_pReal !< average stiffness last increment contains @@ -128,9 +128,9 @@ subroutine basic_init(temperature) close (777) endif mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomdim,F),[3,1,mesh_NcpElems]) - call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,C,temp33_Real,.false.,math_I3) ! constitutive response with no deformation in no time to get reference stiffness + call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,C,C_minmaxAvg,temp33_Real,.false.,math_I3) ! constitutive response with no deformation in no time to get reference stiffness if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness - temp3333_Real = C + temp3333_Real = C_minmaxAvg endif call Utilities_updateGamma(temp3333_Real,.True.) @@ -270,7 +270,7 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C) - if (update_gamma) call Utilities_updateGamma(C,restartWrite) + if (update_gamma) call Utilities_updateGamma(C_minmaxAvg,restartWrite) !-------------------------------------------------------------------------------------------------- ! iteration till converged @@ -294,7 +294,7 @@ type(tSolutionState) function & F_aim_lastIter = F_aim basic_solution%termIll = .false. call Utilities_constitutiveResponse(F_lastInc,F,temperature_bc,timeinc,& - P,C,P_av,ForwardData,rotation_BC) + P,C,C_minmaxAvg,P_av,ForwardData,rotation_BC) basic_solution%termIll = terminallyIll terminallyIll = .false. ForwardData = .false. diff --git a/code/DAMASK_spectral_solverBasicPETSc.f90 b/code/DAMASK_spectral_solverBasicPETSc.f90 index e7f8c19a6..1e13a1910 100644 --- a/code/DAMASK_spectral_solverBasicPETSc.f90 +++ b/code/DAMASK_spectral_solverBasicPETSc.f90 @@ -53,7 +53,7 @@ module DAMASK_spectral_SolverBasicPETSc F_aimDot=0.0_pReal character(len=1024), private :: incInfo real(pReal), private, dimension(3,3,3,3) :: & - C = 0.0_pReal, C_lastInc= 0.0_pReal, & + C = 0.0_pReal, C_minmaxAvg = 0.0_pReal, C_lastInc= 0.0_pReal, & S = 0.0_pReal real(pReal), private :: err_stress, err_div @@ -182,10 +182,10 @@ subroutine basicPETSc_init(temperature) call Utilities_constitutiveResponse(& reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),& reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),& - temperature,0.0_pReal,P,C,temp33_Real,.false.,math_I3) + temperature,0.0_pReal,P,C,C_minmaxAvg,temp33_Real,.false.,math_I3) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back into PETSc if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness - temp3333_Real = C + temp3333_Real = C_minmaxAvg endif call Utilities_updateGamma(temp3333_Real,.True.) @@ -302,7 +302,7 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C) - if (update_gamma) call Utilities_updateGamma(C,restartWrite) + if (update_gamma) call Utilities_updateGamma(C_minmaxAvg,restartWrite) ForwardData = .True. @@ -393,7 +393,7 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response call Utilities_constitutiveResponse(F_lastInc,x_scal,params%temperature,params%timeinc, & - f_scal,C,P_av,ForwardData,params%rotation_BC) + f_scal,C,C_minmaxAvg,P_av,ForwardData,params%rotation_BC) ForwardData = .false. !-------------------------------------------------------------------------------------------------- @@ -433,6 +433,8 @@ subroutine BasicPETSc_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ier math_mul33x33, & math_eigenvalues33, & math_transpose33 + use FEsolving, only: & + terminallyIll implicit none SNES :: snes_local @@ -453,7 +455,8 @@ subroutine BasicPETSc_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ier pAvgDivL2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(P_av,math_transpose33(P_av))))) Converged = (it >= itmin .and. & all([ err_div/pAvgDivL2/err_div_tol, & - err_stress/err_stress_tol ] < 1.0_pReal)) + err_stress/err_stress_tol ] < 1.0_pReal)) & + .or. terminallyIll if (Converged) then reason = 1 diff --git a/code/DAMASK_spectral_utilities.f90 b/code/DAMASK_spectral_utilities.f90 index da666df0a..fff140dcf 100644 --- a/code/DAMASK_spectral_utilities.f90 +++ b/code/DAMASK_spectral_utilities.f90 @@ -573,6 +573,66 @@ real(pReal) function utilities_divergenceRMS() endif end function utilities_divergenceRMS + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculate max of curl of field_fourier +!-------------------------------------------------------------------------------------------------- +real(pReal) function utilities_curlRMS() + use math !< must use the whole module for use of FFTW + use mesh, only: & + res, & + res1_red, & + wgt + + implicit none + integer(pInt) :: i, j, k, l + complex(pReal), dimension(3,3) :: curl_fourier + real(pReal) :: curl_abs + + write(6,'(/,a)') ' ... calculating curl ................................................' + flush(6) +!-------------------------------------------------------------------------------------------------- +! calculating max curl criterion in Fourier space + utilities_curlRMS = 0.0_pReal + + do k = 1_pInt, res(3); do j = 1_pInt, res(2); + do i = 2_pInt, res1_red - 1_pInt + do l = 1_pInt, 3_pInt + curl_fourier(l,1) = (field_fourier(i,j,k,l,3)*xi(2,i,j,k)& + - field_fourier(i,j,k,l,2)*xi(3,i,j,k))*TWOPIIMG + curl_fourier(l,2) = (-field_fourier(i,j,k,l,3)*xi(1,i,j,k)& + +field_fourier(i,j,k,l,1)*xi(3,i,j,k) )*TWOPIIMG + curl_fourier(l,3) = ( field_fourier(i,j,k,l,2)*xi(1,i,j,k)& + -field_fourier(i,j,k,l,1)*xi(2,i,j,k) )*TWOPIIMG + enddo + utilities_curlRMS = utilities_curlRMS + & + 2.0_pReal*sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal) + enddo + do l = 1_pInt, 3_pInt + curl_fourier = (field_fourier(1,j,k,l,3)*xi(2,1,j,k)& + - field_fourier(1,j,k,l,2)*xi(3,1,j,k))*TWOPIIMG + curl_fourier = (-field_fourier(1,j,k,l,3)*xi(1,1,j,k)& + +field_fourier(1,j,k,l,1)*xi(3,1,j,k) )*TWOPIIMG + curl_fourier = ( field_fourier(1,j,k,l,2)*xi(1,1,j,k)& + -field_fourier(1,j,k,l,1)*xi(2,1,j,k) )*TWOPIIMG + enddo + utilities_curlRMS = utilities_curlRMS + & + 2.0_pReal*sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal) + do l = 1_pInt, 3_pInt + curl_fourier = (field_fourier(res1_red,j,k,l,3)*xi(2,res1_red,j,k)& + - field_fourier(res1_red,j,k,l,2)*xi(3,res1_red,j,k))*TWOPIIMG + curl_fourier = (-field_fourier(res1_red,j,k,l,3)*xi(1,res1_red,j,k)& + +field_fourier(res1_red,j,k,l,1)*xi(3,res1_red,j,k) )*TWOPIIMG + curl_fourier = ( field_fourier(res1_red,j,k,l,2)*xi(1,res1_red,j,k)& + -field_fourier(res1_red,j,k,l,1)*xi(2,res1_red,j,k) )*TWOPIIMG + enddo + utilities_curlRMS = utilities_curlRMS + & + 2.0_pReal*sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal) + enddo; enddo + utilities_curlRMS = sqrt(utilities_curlRMS) *wgt + +end function utilities_curlRMS !-------------------------------------------------------------------------------------------------- @@ -677,7 +737,7 @@ end function utilities_maskedCompliance !> @brief calculates constitutive response !-------------------------------------------------------------------------------------------------- subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& - P,C,P_av,forwardData,rotation_BC) + P,C_volAvg,C_minmaxAvg,P_av,forwardData,rotation_BC) use debug, only: & debug_reset, & debug_info @@ -713,7 +773,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& logical, intent(in) :: forwardData !< age results real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame - real(pReal),intent(out), dimension(3,3,3,3) :: C !< average stiffness + real(pReal),intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress real(pReal),intent(out), dimension(3,3,res(1),res(2),res(3)) :: P !< PK stress @@ -723,6 +783,9 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& real(pReal), dimension(3,3,3,3) :: dPdF !< d P / d F real(pReal), dimension(6) :: sigma !< cauchy stress in mandel notation real(pReal), dimension(6,6) :: dsde !< d sigma / d Epsilon + real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF + real(pReal) :: max_dPdF_norm, min_dPdF_norm + integer(pInt) :: k write(6,'(/,a)') ' ... evaluating constitutive response ......................................' calcMode = CPFEM_CALCRESULTS @@ -764,9 +827,26 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& call CPFEM_general(calcMode,F_lastInc(1:3,1:3,1,1,1), F(1:3,1:3,1,1,1), & ! first call calculates everything temperature,timeinc,1_pInt,1_pInt,sigma,dsde,P(1:3,1:3,1,1,1),dPdF) + + max_dPdF = 0.0_pReal + max_dPdF_norm = 0.0_pReal + min_dPdF = huge(1.0_pReal) + min_dPdF_norm = huge(1.0_pReal) + do k = 1_pInt, res(3)*res(2)*res(3) + if (max_dPdF_norm < sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal)) then + max_dPdF = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k) + max_dPdF_norm = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal) + endif + if (min_dPdF_norm > sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal)) then + min_dPdF = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k) + min_dPdF_norm = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal) + endif + enddo P = reshape(materialpoint_P, [3,3,res(1),res(2),res(3)]) - C = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt + C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt + C_minmaxAvg = 0.5_pReal*(max_dPdF + min_dPdF) + call debug_info() restartWrite = .false. ! reset restartWrite status