diff --git a/code/DAMASK_spectral_solverAL.f90 b/code/DAMASK_spectral_solverAL.f90 index d393bf515..885022a38 100644 --- a/code/DAMASK_spectral_solverAL.f90 +++ b/code/DAMASK_spectral_solverAL.f90 @@ -106,7 +106,6 @@ subroutine AL_init(temperature) temperature #include #include - integer(pInt) :: i,j,k real(pReal), dimension(3,3, res(1), res(2),res(3)) :: P real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal @@ -376,7 +375,9 @@ end function AL_solution subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) use numerics, only: & itmax, & - itmin + itmin, & + polarAlpha, & + polarBeta use math, only: & math_rotate_backward33, & math_transpose33, & @@ -395,7 +396,6 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) implicit none integer(pInt), save :: callNo = 3_pInt - real(pReal), dimension(3,3) :: temp33_Real !-------------------------------------------------------------------------------------------------- ! strange syntax in the next line because otherwise macros expand beyond 132 character limit @@ -451,10 +451,29 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) reportIter = reportIter + 1_pInt endif callNo = callNo +1_pInt + +!-------------------------------------------------------------------------------------------------- +! + field_real = 0.0_pReal + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(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_lambda(1:3,1:3,i,j,k)) + enddo; enddo; enddo + +!-------------------------------------------------------------------------------------------------- +! doing convolution in Fourier space + call Utilities_FFTforward() + call Utilities_fourierConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) + call Utilities_FFTbackward() + +!-------------------------------------------------------------------------------------------------- +! constructing residual + residual_F_lambda = polarBeta*F - reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),& + [3,3,res(1),res(2),res(3)],order=[3,4,5,1,2]) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call Utilities_constitutiveResponse(F_lastInc,F,params%temperature,params%timeinc, & + call Utilities_constitutiveResponse(F_lastInc,F - residual_F_lambda/polarBeta,params%temperature,params%timeinc, & residual_F,C,P_av,ForwardData,params%rotation_BC) ForwardData = .False. @@ -462,31 +481,19 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) ! 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 - + !-------------------------------------------------------------------------------------------------- -! - field_real = 0.0_pReal +! constructing residual do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - temp33_Real = math_mul3333xx33(S_scale,residual_F(1:3,1:3,i,j,k)) + math_I3 - residual_F(1:3,1:3,i,j,k) = temp33_Real - field_real(i,j,k,1:3,1:3) = -math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k)-F(1:3,1:3,i,j,k)) + residual_F(1:3,1:3,i,j,k) = math_mul3333xx33(S_scale,residual_F(1:3,1:3,i,j,k)) - & + F_lambda(1:3,1:3,i,j,k) + & + F(1:3,1:3,i,j,k) + & + residual_F_lambda(1:3,1:3,i,j,k) enddo; enddo; enddo -!-------------------------------------------------------------------------------------------------- -! doing convolution in Fourier space - call Utilities_FFTforward() - call Utilities_fourierConvolution(math_rotate_backward33(F_aim,params%rotation_BC)) - call Utilities_FFTbackward() - -!-------------------------------------------------------------------------------------------------- -! constructing residual - residual_F_lambda = F - reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),& - [3,3,res(1),res(2),res(3)],order=[3,4,5,1,2]) - residual_F = residual_F - F_lambda + residual_F_lambda - !-------------------------------------------------------------------------------------------------- ! calculating errors - err_f = wgt*sqrt(sum(residual_F_lambda**2.0_pReal)) + err_f = wgt*sqrt(sum(residual_F_lambda**2.0_pReal))/polarBeta err_p = wgt*sqrt(sum((residual_F - residual_F_lambda)**2.0_pReal)) end subroutine AL_formResidual diff --git a/code/numerics.f90 b/code/numerics.f90 index 3a3aba694..177efd996 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -92,7 +92,9 @@ module numerics err_f_tol = 1e-6_pReal, & err_p_tol = 1e-5_pReal, & 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 + 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 character(len=64), private :: & fftw_plan_mode = 'FFTW_PATIENT' !< reads the planing-rigor flag, see manual on www.fftw.org, Default FFTW_PATIENT: use patient planner flag character(len=64), protected, public :: & @@ -305,9 +307,14 @@ subroutine numerics_init err_f_tol = IO_floatValue(line,positions,2_pInt) case ('err_p_tol') err_p_tol = 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_tol', 'err_p_tol') ! found PETSc parameter, but compiled without PETSc + case ('myspectralsolver', 'petsc_options','err_f_tol', 'err_p_tol', & + 'polaralpha','polarBeta') ! found PETSc parameter, but compiled without PETSc call IO_warning(41_pInt,ext_msg=tag) #endif #endif @@ -315,7 +322,7 @@ subroutine numerics_init 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_tol', 'err_p_tol', 'maxcutback') + 'err_f_tol', 'err_p_tol', 'maxcutback','polaralpha','polarbeta') call IO_warning(40_pInt,ext_msg=tag) #endif case default ! found unknown keyword @@ -425,6 +432,8 @@ subroutine numerics_init #ifdef PETSc write(6,'(a24,1x,es8.1)') ' err_f_tol: ',err_f_tol write(6,'(a24,1x,es8.1)') ' err_p_tol: ',err_p_tol + write(6,'(a24,1x,es8.1)') ' polarAlpha: ',polarAlpha + write(6,'(a24,1x,es8.1)') ' polarBeta: ',polarBeta write(6,'(a24,1x,a)') ' myspectralsolver: ',trim(myspectralsolver) write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_options) #endif @@ -484,6 +493,11 @@ subroutine numerics_init #ifdef PETSc if (err_f_tol <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_f_tol') if (err_p_tol <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_p_tol') + if (polarAlpha <= 0.0_pReal .or. & + polarAlpha > 2.0_pReal) call IO_error(301_pInt,ext_msg='polarAlpha') + if (polarBeta <= 0.0_pReal .or. & + polarBeta > 2.0_pReal) call IO_error(301_pInt,ext_msg='polarBeta') + if (err_p_tol <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_p_tol') #endif #endif if (fixedSeed <= 0_pInt) &