223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

fixed al and pol solver, now checking for div(p) = curl(f) = 0

This commit is contained in:
Martin Diehl 2013-08-07 17:20:05 +00:00
parent cc6c524740
commit ac92b90e0b
7 changed files with 213 additions and 184 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

123
code/DAMASK_spectral_solverAL.f90
View File

@ -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
!--------------------------------------------------------------------------------------------------

39
code/DAMASK_spectral_solverBasic.f90
View File

@ -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

46
code/DAMASK_spectral_solverBasicPETSc.f90
View File

@ -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

107
code/DAMASK_spectral_solverPolarisation.f90
View File

@ -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

14
code/DAMASK_spectral_utilities.f90
View File

@ -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

3
code/math.f90
View File

@ -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

65
code/numerics.f90
View File

@ -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