diff --git a/code/DAMASK_spectral_solverAL.f90 b/code/DAMASK_spectral_solverAL.f90 index b00ee815d..d061d6045 100644 --- a/code/DAMASK_spectral_solverAL.f90 +++ b/code/DAMASK_spectral_solverAL.f90 @@ -90,8 +90,8 @@ module DAMASK_spectral_solverAL real(pReal), private :: & err_stress, & !< deviation from stress BC - err_f, & !< difference between compatible and incompatible deformation gradient - err_p !< difference of stress resulting from compatible and incompatible F + err_curl, & !< RMS of curl of F + err_div !< RMS of div of P logical, private :: ForwardData integer(pInt), private :: & totalIter = 0_pInt !< total iteration in current increment @@ -470,6 +470,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) math_transpose33, & math_mul3333xx33, & math_invSym3333, & + math_mul33x33, & PI use DAMASK_spectral_Utilities, only: & grid, & @@ -481,7 +482,9 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) Utilities_fourierConvolution, & Utilities_inverseLaplace, & Utilities_FFTbackward, & - Utilities_constitutiveResponse + Utilities_constitutiveResponse, & + Utilities_divergenceRMS, & + Utilities_curlRMS use debug, only: & debug_level, & debug_spectral, & @@ -513,7 +516,6 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) PetscErrorCode :: ierr integer(pInt) :: & i, j, k, e - real(pReal) :: correctionFactor F => x_scal(1:3,1:3,1,& XG_RANGE,YG_RANGE,ZG_RANGE) @@ -550,7 +552,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) residual_F = params%density*product(geomSize/grid)*residual_F field_real = 0.0_pReal field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3) = reshape(residual_F,[grid(1),grid(2),grid(3),3,3],& - order=[4,5,1,2,3]) ! field real has a different order + order=[4,5,1,2,3]) ! field real has a different order call Utilities_FFTforward() call Utilities_inverseLaplace() inertiaField_fourier = field_fourier @@ -562,8 +564,11 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) ! field_real = 0.0_pReal do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) - field_real(i,j,k,1:3,1:3) = math_mul3333xx33(C_scale,polarBeta*F(1:3,1:3,i,j,k) - & - polarAlpha*F_lambda(1:3,1:3,i,j,k)) + field_real(i,j,k,1:3,1:3) = & + polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& + polarAlpha*math_mul33x33(F(1:3,1:3,i,j,k), & + math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k) - math_I3)) + enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- @@ -584,34 +589,35 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) call Utilities_constitutiveResponse(F_lastInc,F - residual_F_lambda/polarBeta,params%temperature,params%timeinc, & residual_F,C_volAvg,C_minMaxAvg,P_av,ForwardData,params%rotation_BC) ForwardData = .False. - + !-------------------------------------------------------------------------------------------------- -! stress BC handling - write(6,'(/,a)') ' ... correcting F to fullfill stress BC ....................................' - correctionFactor = (cos((1.0-10000.0_pReal**(-sum((P_av-P_avLastEval)**2.0_pReal)/& ! only correct when averages stress of last two calls doesn't strongly deviate - sum(P_av**2.0_pReal)))*PI)+1.0)/2.0_pReal - write(6,'(/,a,f10.4)') ' stress BC correction factor = ', correctionFactor - F_aim = F_aim - correctionFactor *math_mul3333xx33(S, ((P_av - params%P_BC))) ! S = 0.0 for no bc - err_stress = maxval(abs(mask_stress * (P_av - params%P_BC))) ! mask = 0.0 for no bc +! calculate divergence + field_real = 0.0_pReal + field_real = reshape(residual_F,[grid(1),grid(2),grid(3),3,3],order=[4,5,1,2,3]) + call Utilities_FFTforward() + err_div = Utilities_divergenceRMS() + call Utilities_FFTbackward() !-------------------------------------------------------------------------------------------------- ! constructing residual e = 0_pInt - err_p = 0.0_pReal do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) e = e + 1_pInt residual_F(1:3,1:3,i,j,k) = math_mul3333xx33(math_invSym3333(materialpoint_dPdF(:,:,:,:,1,e) + C_scale), & residual_F(1:3,1:3,i,j,k) - & - math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k) - math_I3)) & - + residual_F_lambda(1:3,1:3,i,j,k) - err_p = err_p + sum((math_mul3333xx33(C_scale,residual_F(1:3,1:3,i,j,k) - residual_F_lambda(1:3,1:3,i,j,k)))**2.0_pReal) + math_mul33x33(F(1:3,1:3,i,j,k), & + math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k) - math_I3))) & + + residual_F_lambda(1:3,1:3,i,j,k) enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- -! calculating errors - err_f = wgt*sqrt(sum(residual_F_lambda**2.0_pReal))/polarBeta - err_p = wgt*sqrt(err_p) - +! calculating curl + field_real = 0.0_pReal + field_real = reshape(F,[grid(1),grid(2),grid(3),3,3],order=[4,5,1,2,3]) + call Utilities_FFTforward() + err_curl = Utilities_curlRMS() + call Utilities_FFTbackward() + end subroutine AL_formResidual @@ -619,14 +625,17 @@ end subroutine AL_formResidual !> @brief convergence check !-------------------------------------------------------------------------------------------------- subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) -use numerics, only: & - itmax, & - itmin, & - err_f_tolabs, & - err_p_tolabs, & - err_f_p_tolrel, & - err_stress_tolabs, & - err_stress_tolrel + use math, only: & + math_mul3333xx33 + use numerics, only: & + itmax, & + itmin, & + err_div_tolRel, & + err_div_tolAbs, & + err_curl_tolRel, & + err_curl_tolAbs, & + err_stress_tolabs, & + err_stress_tolrel use FEsolving, only: & terminallyIll @@ -641,31 +650,23 @@ use numerics, only: & PetscObject :: dummy PetscErrorCode ::ierr real(pReal) :: & - mismatch_f_tol, & - mismatch_p_tol, & - stressBC_tol - - mismatch_f_tol = max(maxval(abs(F_aim-math_I3))*err_f_p_tolrel,err_f_tolabs) - mismatch_p_tol = max(maxval(abs(P_av)) *err_f_p_tolrel,err_p_tolabs) - stressBC_tol = max(maxval(abs(P_av)) *err_stress_tolrel,err_stress_tolabs) - - write(6,'(1/,a)') ' ... reporting .............................................................' - write(6,'(/,a,f8.2,a,es11.5,a,es11.4,a)') ' mismatch F = ', & - err_f/mismatch_f_tol, & - ' (',err_f,' -, tol =',mismatch_f_tol,')' - write(6,'(a,f8.2,a,es11.5,a,es11.4,a)') ' mismatch P = ', & - err_p/mismatch_p_tol, & - ' (',err_p,' -, tol =',mismatch_p_tol,')' - write(6,'(a,f8.2,a,es11.5,a,es11.4,a)') ' error stress BC = ', & - err_stress/stressBC_tol, ' (',err_stress, ' Pa, tol =',stressBC_tol,')' - write(6,'(/,a)') ' ===========================================================================' - flush(6) - + curlTol, & + divTol, & + stressTol + +!-------------------------------------------------------------------------------------------------- +! stress BC handling + F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%P_BC))) ! S = 0.0 for no bc + err_stress = maxval(abs(mask_stress * (P_av - params%P_BC))) ! mask = 0.0 for no bc + curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel,err_curl_tolAbs) + divTol = max(maxval(abs(P_av)) *err_div_tolRel,err_div_tolAbs) + stressTol = max(maxval(abs(P_av)) *err_stress_tolrel,err_stress_tolabs) + converged: if ((totalIter >= itmin .and. & - all([ err_f/mismatch_f_tol, & - err_p/mismatch_p_tol, & - err_stress/stressBC_tol] < 1.0_pReal)) & - .or. terminallyIll) then + all([ err_div/divTol, & + err_curl/curlTol, & + err_stress/err_stress_tol ] < 1.0_pReal)) & + .or. terminallyIll) reason = 1 elseif (totalIter >= itmax) then converged reason = -1 @@ -673,6 +674,18 @@ use numerics, only: & reason = 0 endif converged +!-------------------------------------------------------------------------------------------------- +! report + write(6,'(1/,a)') ' ... reporting .............................................................' + write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & + err_curl/curlTol,' (',err_curl,' -, tol =',curlTol,')' + write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + err_div/divTol, ' (',err_div,' / m, tol =',divTol,')' + write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & + err_stress/stressTol, ' (',err_stress, ' Pa, tol =',stressTol,')' + write(6,'(/,a)') ' ===========================================================================' + flush(6) + end subroutine AL_converged !-------------------------------------------------------------------------------------------------- diff --git a/code/DAMASK_spectral_solverBasic.f90 b/code/DAMASK_spectral_solverBasic.f90 index 328779890..724e7fd89 100644 --- a/code/DAMASK_spectral_solverBasic.f90 +++ b/code/DAMASK_spectral_solverBasic.f90 @@ -337,8 +337,8 @@ type(tSolutionState) function basic_solution(& call Utilities_fourierConvolution(math_rotate_backward33(F_aim_lastIter-F_aim,rotation_BC)) call Utilities_FFTbackward() F = F - reshape(field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3),shape(F),order=[3,4,5,1,2]) ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization - basic_solution%converged = basic_Converged(err_div,P_av,err_stress,P_av) - write(6,'(/,a)') ' ==========================================================================' + basic_solution%converged = basic_Converged(err_div,P_av,err_stress) + write(6,'(/,a)') ' ===========================================================================' flush(6) if ((basic_solution%converged .and. iter >= itmin) .or. basic_solution%termIll) then basic_solution%iterationsNeeded = iter @@ -352,10 +352,11 @@ end function basic_solution !-------------------------------------------------------------------------------------------------- !> @brief convergence check for basic scheme based on div of P and deviation from stress aim !-------------------------------------------------------------------------------------------------- -logical function basic_Converged(err_div,pAvgDiv,err_stress,pAvgStress) +logical function basic_Converged(err_div,P_av,err_stress) use numerics, only: & itmin, & - err_div_tol, & + err_div_tolRel, & + err_div_tolAbs, & err_stress_tolrel, & err_stress_tolabs use math, only: & @@ -365,28 +366,30 @@ logical function basic_Converged(err_div,pAvgDiv,err_stress,pAvgStress) implicit none real(pReal), dimension(3,3), intent(in) :: & - pAvgDiv,& - pAvgStress + P_av real(pReal), intent(in) :: & err_div, & err_stress real(pReal) :: & - err_stress_tol, & - pAvgDivL2 + stressTol, & + divTol - pAvgDivL2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(pAvgDiv,math_transpose33(pAvgDiv))))) ! L_2 norm of average stress (http://mathworld.wolfram.com/SpectralNorm.html) - err_stress_tol = max(maxval(abs(pAvgStress))*err_stress_tolrel,err_stress_tolabs) + divTol = max(maxval(abs(P_av))*err_div_tolRel,err_div_tolAbs) + stressTol = max(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs) - basic_Converged = all([ err_div/pAvgDivL2/err_div_tol,& - err_stress/err_stress_tol ] < 1.0_pReal) - - write(6,'(/,a,f10.2,a,es11.5,a,es11.4,a)') ' error divergence = ', & - err_div/pAvgDivL2/err_div_tol, ' (',err_div/pAvgDivL2,' / m, tol =',err_div_tol,')' - write(6, '(a,f10.2,a,es11.5,a,es11.4,a)') ' error stress BC = ', & - err_stress/err_stress_tol, ' (',err_stress, ' Pa, tol =',err_stress_tol,')' - flush(6) + basic_Converged = all([ err_div/divTol,& + err_stress/stressTol ] < 1.0_pReal) + +!-------------------------------------------------------------------------------------------------- +! report + write(6,'(1/,a)') ' ... reporting .............................................................' + write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + err_div/divTol, ' (',err_div,' / m, tol =',divTol,')' + write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & + err_stress/stressTol, ' (',err_stress, ' Pa, tol =',stressTol,')' + flush(6) end function basic_Converged diff --git a/code/DAMASK_spectral_solverBasicPETSc.f90 b/code/DAMASK_spectral_solverBasicPETSc.f90 index cdaa7e0e2..5fb3863b7 100644 --- a/code/DAMASK_spectral_solverBasicPETSc.f90 +++ b/code/DAMASK_spectral_solverBasicPETSc.f90 @@ -495,9 +495,10 @@ subroutine BasicPETSc_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ier use numerics, only: & itmax, & itmin, & - err_div_tol, & - err_stress_tolrel, & - err_stress_tolabs + err_div_TolAbs, & + err_div_TolRel, & + err_stress_tolRel, & + err_stress_tolAbs use math, only: & math_mul33x33, & math_eigenvalues33, & @@ -515,32 +516,33 @@ subroutine BasicPETSc_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ier SNESConvergedReason :: reason PetscObject :: dummy PetscErrorCode :: ierr - logical :: Converged real(pReal) :: & - pAvgDivL2, & - err_stress_tol + divTol, & + stressTol - err_stress_tol =max(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs) - pAvgDivL2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(P_av,math_transpose33(P_av))))) - Converged = (it >= itmin .and. & - all([ err_div/pAvgDivL2/err_div_tol, & + divTol = max(maxval(abs(P_av))*err_div_tolRel,err_div_tolAbs) + stressTol = max(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs) + + converged: if ((it >= itmin .and. & + all([ err_div/divTol, & err_stress/err_stress_tol ] < 1.0_pReal)) & - .or. terminallyIll - - if (Converged) then + .or. terminallyIll) reason = 1 - elseif (it >= itmax) then + elseif (totalIter >= itmax) then converged reason = -1 - else + else converged reason = 0 - endif + endif converged + +!-------------------------------------------------------------------------------------------------- +! report write(6,'(1/,a)') ' ... reporting .............................................................' - write(6,'(/,a,f8.2,a,es11.5,a,es11.4,a)') ' error divergence = ', & - err_div/pAvgDivL2/err_div_tol, ' (',err_div/pAvgDivL2,' / m, tol =',err_div_tol,')' - write(6,'(a,f8.2,a,es11.5,a,es11.4,a)') ' error stress BC = ', & - err_stress/err_stress_tol, ' (',err_stress, ' Pa , tol =',err_stress_tol,')' - write(6,'(/,a)') ' ===========================================================================' - flush(6) + write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + err_div/divTol, ' (',err_div,' / m, tol =',divTol,')' + write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & + err_stress/stressTol, ' (',err_stress, ' Pa, tol =',stressTol,')' + write(6,'(/,a)') ' ===========================================================================' + flush(6) end subroutine BasicPETSc_converged diff --git a/code/DAMASK_spectral_solverPolarisation.f90 b/code/DAMASK_spectral_solverPolarisation.f90 index 23b265f64..5e71f3c3c 100644 --- a/code/DAMASK_spectral_solverPolarisation.f90 +++ b/code/DAMASK_spectral_solverPolarisation.f90 @@ -79,7 +79,7 @@ module DAMASK_spectral_solverPolarisation F_aim_lastInc = math_I3, & !< previous average deformation gradient P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress P_avLastEval = 0.0_pReal !< average 1st Piola--Kirchhoff stress last call of CPFEM_general - character(len=1024), private :: incInfo !< time and increment information + character(len=1024), private :: incInfo !< time and increment information real(pReal), private, dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness @@ -90,8 +90,8 @@ module DAMASK_spectral_solverPolarisation real(pReal), private :: & err_stress, & !< deviation from stress BC - err_f, & !< difference between compatible and incompatible deformation gradient - err_p !< difference of stress resulting from compatible and incompatible F + err_curl, & !< RMS of curl of F + err_div !< RMS of div of P logical, private :: ForwardData integer(pInt), private :: & totalIter = 0_pInt !< total iteration in current increment @@ -470,6 +470,7 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) math_transpose33, & math_mul3333xx33, & math_invSym3333, & + math_mul33x33, & PI use DAMASK_spectral_Utilities, only: & grid, & @@ -481,7 +482,9 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) Utilities_fourierConvolution, & Utilities_inverseLaplace, & Utilities_FFTbackward, & - Utilities_constitutiveResponse + Utilities_constitutiveResponse, & + Utilities_divergenceRMS, & + Utilities_curlRMS use debug, only: & debug_level, & debug_spectral, & @@ -513,7 +516,6 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) PetscErrorCode :: ierr integer(pInt) :: & i, j, k, e - real(pReal) :: correctionFactor F => x_scal(1:3,1:3,1,& XG_RANGE,YG_RANGE,ZG_RANGE) @@ -562,8 +564,11 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) ! field_real = 0.0_pReal do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) - field_real(i,j,k,1:3,1:3) = math_mul3333xx33(C_scale,(polarAlpha + polarBeta)*F(1:3,1:3,i,j,k) - & - (polarAlpha)*F_tau(1:3,1:3,i,j,k)) + field_real(i,j,k,1:3,1:3) = & + polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& + polarAlpha*math_mul33x33(F(1:3,1:3,i,j,k), & + math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3)) + enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- @@ -572,7 +577,7 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) field_fourier = field_fourier + polarAlpha*inertiaField_fourier call Utilities_fourierConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) call Utilities_FFTbackward() - + !-------------------------------------------------------------------------------------------------- ! constructing residual residual_F_tau = polarBeta*F - reshape(field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3),& @@ -584,15 +589,14 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) call Utilities_constitutiveResponse(F_lastInc,F - residual_F_tau/polarBeta,params%temperature,params%timeinc, & residual_F,C_volAvg,C_minMaxAvg,P_av,ForwardData,params%rotation_BC) ForwardData = .False. - + !-------------------------------------------------------------------------------------------------- -! stress BC handling - write(6,'(/,a)') ' ... correcting F to fullfill stress BC ....................................' - correctionFactor = (cos((1.0-10000.0_pReal**(-sum((P_av-P_avLastEval)**2.0_pReal)/& ! only correct when averages stress of last two calls doesn't strongly deviate - sum(P_av**2.0_pReal)))*PI)+1.0)/2.0_pReal - write(6,'(/,a,f10.4)') ' stress BC correction factor = ', correctionFactor - F_aim = F_aim - correctionFactor *math_mul3333xx33(S, ((P_av - params%P_BC))) ! S = 0.0 for no bc - err_stress = maxval(abs(mask_stress * (P_av - params%P_BC))) ! mask = 0.0 for no bc +! calculate divergence + field_real = 0.0_pReal + field_real = reshape(residual_F,[grid(1),grid(2),grid(3),3,3],order=[4,5,1,2,3]) + call Utilities_FFTforward() + err_div = Utilities_divergenceRMS() + call Utilities_FFTbackward() !-------------------------------------------------------------------------------------------------- ! constructing residual @@ -602,15 +606,18 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) e = e + 1_pInt residual_F(1:3,1:3,i,j,k) = math_mul3333xx33(math_invSym3333(materialpoint_dPdF(:,:,:,:,1,e) + C_scale), & residual_F(1:3,1:3,i,j,k) - & - math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3)) & - + residual_F_tau(1:3,1:3,i,j,k) - err_p = err_p + sum((math_mul3333xx33(C_scale,residual_F(1:3,1:3,i,j,k) - residual_F_tau(1:3,1:3,i,j,k)))**2.0_pReal) + math_mul33x33(F(1:3,1:3,i,j,k), & + math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) & + + residual_F_tau(1:3,1:3,i,j,k) enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- -! calculating errors - err_f = wgt*sqrt(sum(residual_F_tau**2.0_pReal))/polarBeta - err_p = wgt*sqrt(err_p) +! calculating curl + field_real = 0.0_pReal + field_real = reshape(F,[grid(1),grid(2),grid(3),3,3],order=[4,5,1,2,3]) + call Utilities_FFTforward() + err_curl = Utilities_curlRMS() + call Utilities_FFTbackward() end subroutine Polarisation_formResidual @@ -622,9 +629,10 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason, use numerics, only: & itmax, & itmin, & - err_f_tolabs, & - err_p_tolabs, & - err_f_p_tolrel, & + err_div_tolRel, & + err_div_tolAbs, & + err_curl_tolRel, & + err_curl_tolAbs, & err_stress_tolabs, & err_stress_tolrel use FEsolving, only: & @@ -641,31 +649,23 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason, PetscObject :: dummy PetscErrorCode ::ierr real(pReal) :: & - mismatch_f_tol, & - mismatch_p_tol, & + curlTol, & + divTol, & stressBC_tol - - mismatch_f_tol = max(maxval(abs(F_aim-math_I3))*err_f_p_tolrel,err_f_tolabs) - mismatch_p_tol = max(maxval(abs(P_av)) *err_f_p_tolrel,err_p_tolabs) - stressBC_tol = max(maxval(abs(P_av)) *err_stress_tolrel,err_stress_tolabs) - - write(6,'(1/,a)') ' ... reporting .............................................................' - write(6,'(/,a,f8.2,a,es11.5,a,es11.4,a)') ' mismatch F = ', & - err_f/mismatch_f_tol, & - ' (',err_f,' -, tol =',mismatch_f_tol,')' - write(6,'(a,f8.2,a,es11.5,a,es11.4,a)') ' mismatch P = ', & - err_p/mismatch_p_tol, & - ' (',err_p,' -, tol =',mismatch_p_tol,')' - write(6,'(a,f8.2,a,es11.5,a,es11.4,a)') ' error stress BC = ', & - err_stress/stressBC_tol, ' (',err_stress, ' Pa, tol =',stressBC_tol,')' - write(6,'(/,a)') ' ===========================================================================' - flush(6) - + +!-------------------------------------------------------------------------------------------------- +! stress BC handling + F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%P_BC))) ! S = 0.0 for no bc + err_stress = maxval(abs(mask_stress * (P_av - params%P_BC))) ! mask = 0.0 for no bc + curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel,err_curl_tolAbs) + divTol = max(maxval(abs(P_av)) *err_div_tolRel,err_div_tolAbs) + stressBC_tol = max(maxval(abs(P_av)) *err_stress_tolrel,err_stress_tolabs) + converged: if ((totalIter >= itmin .and. & - all([ err_f/mismatch_f_tol, & - err_p/mismatch_p_tol, & - err_stress/stressBC_tol] < 1.0_pReal)) & - .or. terminallyIll) then + all([ err_div/divTol, & + err_curl/curlTol, & + err_stress/err_stress_tol ] < 1.0_pReal)) & + .or. terminallyIll) reason = 1 elseif (totalIter >= itmax) then converged reason = -1 @@ -673,6 +673,17 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason, reason = 0 endif converged +!-------------------------------------------------------------------------------------------------- +! report + write(6,'(1/,a)') ' ... reporting .............................................................' + write(6,'(/,a,f12.2,a,es12.2,a,es12.2,a)') ' error curl = ', & + err_curl/curlTol,' (',err_curl,' -, tol =',curlTol,')' + write(6,'(a,f12.2,a,es12.8,a,es12.8,a)') ' error divergence = ', & + err_div/divTol, ' (',err_div,' / m, tol =',divTol,')' + write(6,'(a,f8.2,a,es11.5,a,es11.4,a)') ' error stress BC = ', & + err_stress/stressBC_tol, ' (',err_stress, ' Pa, tol =',stressBC_tol,')' + write(6,'(/,a)') ' ===========================================================================' + flush(6) end subroutine Polarisation_converged diff --git a/code/DAMASK_spectral_utilities.f90 b/code/DAMASK_spectral_utilities.f90 index aec69b4c4..1298cd8b7 100644 --- a/code/DAMASK_spectral_utilities.f90 +++ b/code/DAMASK_spectral_utilities.f90 @@ -690,7 +690,7 @@ real(pReal) function utilities_curlRMS() 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 + utilities_curlRMS = sqrt(utilities_curlRMS) * wgt end function utilities_curlRMS @@ -709,10 +709,10 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) math_invert implicit none - real(pReal), dimension(3,3,3,3) :: utilities_maskedCompliance !< masked compliance - real(pReal), intent(in) , dimension(3,3,3,3) :: C !< current average stiffness - real(pReal), intent(in) , dimension(3,3) :: rot_BC !< rotation of load frame - logical, intent(in), dimension(3,3) :: mask_stress !< mask of stress BC + real(pReal), dimension(3,3,3,3) :: utilities_maskedCompliance !< masked compliance + real(pReal), intent(in) , dimension(3,3,3,3) :: C !< current average stiffness + real(pReal), intent(in) , dimension(3,3) :: rot_BC !< rotation of load frame + logical, intent(in), dimension(3,3) :: mask_stress !< mask of stress BC integer(pInt) :: j, k, m, n logical, dimension(9) :: mask_stressVector real(pReal), dimension(9,9) :: temp99_Real @@ -912,10 +912,10 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P if (debugRotation) & write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',& - math_transpose33(P_av)/1.e6_pReal + math_transpose33(P_av)*1.e-6_pReal P_av = math_rotate_forward33(P_av,rotation_BC) write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',& - math_transpose33(P_av)/1.e6_pReal + math_transpose33(P_av)*1.e-6_pReal end subroutine utilities_constitutiveResponse diff --git a/code/math.f90 b/code/math.f90 index 92d43bb98..e445631b6 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -582,7 +582,7 @@ end function math_mul6x6 !-------------------------------------------------------------------------------------------------- -!> @brief matrix multiplication 33x33 = 1 (double contraction --> ij * ij) +!> @brief matrix multiplication 33xx33 = 1 (double contraction --> ij * ij) !-------------------------------------------------------------------------------------------------- real(pReal) pure function math_mul33xx33(A,B) @@ -2320,7 +2320,6 @@ pure subroutine math_hi(M,HI1M,HI2M,HI3M) HI2M=HI1M**2.0_pReal/2.0_pReal- (M(1,1)**2.0_pReal+M(2,2)**2.0_pReal+M(3,3)**2.0_pReal)& /2.0_pReal-M(1,2)*M(2,1)-M(1,3)*M(3,1)-M(2,3)*M(3,2) HI3M=math_det33(M) -! QUESTION: is 3rd equiv det(M) ?? if yes, use function math_det !agreed on YES end subroutine math_hi diff --git a/code/numerics.f90 b/code/numerics.f90 index a70bb2436..cc5755d29 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -77,20 +77,20 @@ module numerics 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 logical, protected, public :: & - analyticJaco = .false., & !< use analytic Jacobian or perturbation, Default .false.: calculate Jacobian using perturbations + analyticJaco = .false., & !< use analytic Jacobian or perturbation, Default .false.: calculate Jacobian using perturbations usePingPong = .true., & - numerics_timeSyncing = .false. !< flag indicating if time synchronization in crystallite is used for nonlocal plasticity + numerics_timeSyncing = .false. !< flag indicating if time synchronization in crystallite is used for nonlocal plasticity !-------------------------------------------------------------------------------------------------- ! spectral parameters: #ifdef Spectral real(pReal), protected, public :: & - err_div_tol = 5.0e-4_pReal, & !< Div(P)/avg(P)*meter - err_stress_tolrel = 0.01_pReal, & !< relative tolerance for fullfillment of stress BC - err_stress_tolabs = 1.0e3_pReal, & !< absolute tolerance for fullfillment of stress BC - err_f_tolabs = 1.0e-8_pReal, & !< absolute tolerance mismatch F - err_p_tolabs = 1.0e2_pReal, & !< absolute tolerance mismatch P - err_f_p_tolrel = 1.0e-5_pReal, & !< relative tolerance for mismatch F and P + err_div_tolAbs = 1.0e-10_pReal, & !< absolute tolerance for equilibrium + err_div_tolRel = 5.0e-4_pReal, & !< relative tolerance for equilibrium + err_curl_tolAbs = 1.0e-10_pReal, & !< absolute tolerance for compatibility + err_curl_tolRel = 5.0e-4_pReal, & !< relative tolerance for compatibility + 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 fftw_timelimit = -1.0_pReal, & !< sets the timelimit of plan creation for FFTW, see manual on www.fftw.org, Default -1.0: disable timelimit 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 @@ -277,8 +277,10 @@ subroutine numerics_init !-------------------------------------------------------------------------------------------------- ! spectral parameters #ifdef Spectral - case ('err_div_tol') - err_div_tol = IO_floatValue(line,positions,2_pInt) + case ('err_div_tolabs') + err_div_tolAbs = IO_floatValue(line,positions,2_pInt) + case ('err_div_tolrel') + err_div_tolRel = IO_floatValue(line,positions,2_pInt) case ('err_stress_tolrel') err_stress_tolrel = IO_floatValue(line,positions,2_pInt) case ('err_stress_tolabs') @@ -310,28 +312,27 @@ subroutine numerics_init petsc_options = trim(line(positions(4):)) case ('myspectralsolver') myspectralsolver = IO_lc(IO_stringValue(line,positions,2_pInt)) - case ('err_f_tolabs') - err_f_tolabs = IO_floatValue(line,positions,2_pInt) - case ('err_p_tolabs') - err_p_tolabs = IO_floatValue(line,positions,2_pInt) - case ('err_f_p_tolrel') - err_f_p_tolrel = IO_floatValue(line,positions,2_pInt) + case ('err_curl_tolabs') + err_curl_tolAbs = IO_floatValue(line,positions,2_pInt) + case ('err_curl_tolrel') + err_curl_tolRel = IO_floatValue(line,positions,2_pInt) case ('polaralpha') polarAlpha = IO_floatValue(line,positions,2_pInt) case ('polarbeta') polarBeta = IO_floatValue(line,positions,2_pInt) #endif #ifndef PETSc - case ('myspectralsolver', 'petsc_options','err_f_tolabs', 'err_p_tolabs', & - 'err_f_p_tolrel','polaralpha','polarBeta') ! found PETSc parameter, but compiled without PETSc + case ('myspectralsolver', 'petsc_options', & + 'err_curl_tolabs','err_curl_tolrel','polaralpha','polarBeta') ! found PETSc parameter, but compiled without PETSc call IO_warning(41_pInt,ext_msg=tag) #endif #endif #ifndef Spectral - case ('err_div_tol','err_stress_tolrel','err_stress_tolabs',& ! found spectral parameter for FEM build - 'itmax', 'itmin','memory_efficient','fftw_timelimit','fftw_plan_mode','myspectralsolver', & - 'rotation_tol','divergence_correction','update_gamma','petsc_options','myfilter', & - 'err_f_tolabs', 'err_p_tolabs', 'err_f_p_tolrel', 'maxcutback','polaralpha','polarbeta') + case ('err_div_tolabs','err_div_tolrel','err_stress_tolrel','err_stress_tolabs',& ! found spectral parameter for FEM build + 'itmax', 'itmin','memory_efficient','fftw_timelimit','fftw_plan_mode', & + 'rotation_tol','divergence_correction','update_gamma','myfilter', & + 'err_div_tolabs','err_div_tolrel', & + 'maxcutback','polaralpha','polarbeta') call IO_warning(40_pInt,ext_msg=tag) #endif case default ! found unknown keyword @@ -419,7 +420,8 @@ subroutine numerics_init !-------------------------------------------------------------------------------------------------- ! spectral parameters #ifdef Spectral - write(6,'(a24,1x,es8.1)') ' err_div_tol: ',err_div_tol + write(6,'(a24,1x,es8.1)') ' err_div_tolAbs: ',err_div_tolAbs + write(6,'(a24,1x,es8.1)') ' err_div_tolRel: ',err_div_tolRel write(6,'(a24,1x,es8.1)') ' err_stress_tolrel: ',err_stress_tolrel write(6,'(a24,1x,es8.1)') ' err_stress_tolabs: ',err_stress_tolabs @@ -441,9 +443,8 @@ subroutine numerics_init write(6,'(a24,1x,i8)') ' divergence_correction: ',divergence_correction write(6,'(a24,1x,L8,/)') ' update_gamma: ',update_gamma #ifdef PETSc - write(6,'(a24,1x,es8.1)') ' err_f_tolabs: ',err_f_tolabs - write(6,'(a24,1x,es8.1)') ' err_p_tolabs: ',err_p_tolabs - write(6,'(a24,1x,es8.1)') ' err_f_p_tolrel: ',err_f_p_tolrel + write(6,'(a24,1x,es8.1)') ' err_curl_tolAbs: ',err_curl_tolAbs + write(6,'(a24,1x,es8.1)') ' err_curl_tolRel: ',err_curl_tolRel write(6,'(a24,1x,es8.1)') ' polarAlpha: ',polarAlpha write(6,'(a24,1x,es8.1)') ' polarBeta: ',polarBeta write(6,'(a24,1x,a)') ' myspectralsolver: ',trim(myspectralsolver) @@ -492,23 +493,23 @@ subroutine numerics_init if (volDiscrMod_RGC < 0.0_pReal) call IO_error(301_pInt,ext_msg='volDiscrMod_RGC') if (volDiscrPow_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='volDiscrPw_RGC') #ifdef Spectral - if (err_div_tol <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_div_tol') + if (err_div_tolRel <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_div_tolRel') + if (err_div_tolAbs <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_div_tolAbs') if (err_stress_tolrel <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_stress_tolrel') if (err_stress_tolabs <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_stress_tolabs') if (itmax <= 1_pInt) call IO_error(301_pInt,ext_msg='itmax') - if (itmin > itmax .or. itmin < 1_pInt) call IO_error(301_pInt,ext_msg='itmin') + if (itmin > itmax .or. itmin < 1_pInt) call IO_error(301_pInt,ext_msg='itmin') if (divergence_correction < 0_pInt .or. & divergence_correction > 2_pInt) call IO_error(301_pInt,ext_msg='divergence_correction') if (maxCutBack < 0_pInt) call IO_error(301_pInt,ext_msg='maxCutBack') if (update_gamma .and. & .not. memory_efficient) call IO_error(error_ID = 847_pInt) #ifdef PETSc - if (err_f_tolabs <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_f_tolabs') - if (err_p_tolabs <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_p_tolabs') - if (err_f_p_tolrel <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_f_p_tolrel') + if (err_curl_tolRel <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_curl_tolRel') + if (err_curl_tolAbs <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_curl_tolAbs') if (polarAlpha <= 0.0_pReal .or. & polarAlpha > 2.0_pReal) call IO_error(301_pInt,ext_msg='polarAlpha') - if (polarBeta <= 0.0_pReal .or. & + if (polarBeta < 0.0_pReal .or. & polarBeta > 2.0_pReal) call IO_error(301_pInt,ext_msg='polarBeta') #endif #endif