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

restructured AL solver to resemble polarization scheme and introduced alpha and beta parameters in numerics to determine what kind of solution you want. alpha = beta = 1 results in the AL scheme. alpha = beta = 2 results in the accelerated scheme. safe values are anything in between. changed the solver to use stress calculations on the campatible field which seems to result in faster convergence.

This commit is contained in:
Pratheek Shanthraj 2013-02-28 17:35:02 +00:00
parent 1ad85a350b
commit 0209a1d5a7
2 changed files with 47 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

53
code/DAMASK_spectral_solverAL.f90
View File

@ -106,7 +106,6 @@ subroutine AL_init(temperature)
temperature temperature
#include <finclude/petscdmda.h90> #include <finclude/petscdmda.h90>
#include <finclude/petscsnes.h90> #include <finclude/petscsnes.h90>
integer(pInt) :: i,j,k
real(pReal), dimension(3,3, res(1), res(2),res(3)) :: P real(pReal), dimension(3,3, res(1), res(2),res(3)) :: P
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
temp33_Real = 0.0_pReal temp33_Real = 0.0_pReal
@ -376,7 +375,9 @@ end function AL_solution
subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
use numerics, only: & use numerics, only: &
itmax, & itmax, &
itmin itmin, &
polarAlpha, &
polarBeta
use math, only: & use math, only: &
math_rotate_backward33, & math_rotate_backward33, &
math_transpose33, & math_transpose33, &
@ -395,7 +396,6 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
implicit none implicit none
integer(pInt), save :: callNo = 3_pInt 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 ! 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 reportIter = reportIter + 1_pInt
endif endif
callNo = callNo +1_pInt 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 ! 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) residual_F,C,P_av,ForwardData,params%rotation_BC)
ForwardData = .False. ForwardData = .False.
@ -462,31 +481,19 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
! stress BC handling ! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%P_BC))) ! S = 0.0 for no bc 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 err_stress = maxval(abs(mask_stress * (P_av - params%P_BC))) ! mask = 0.0 for no bc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! ! constructing residual
field_real = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) 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) = math_mul3333xx33(S_scale,residual_F(1:3,1:3,i,j,k)) - &
residual_F(1:3,1:3,i,j,k) = temp33_Real F_lambda(1:3,1:3,i,j,k) + &
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)) F(1:3,1:3,i,j,k) + &
residual_F_lambda(1:3,1:3,i,j,k)
enddo; enddo; enddo 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 ! 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)) err_p = wgt*sqrt(sum((residual_F - residual_F_lambda)**2.0_pReal))
end subroutine AL_formResidual end subroutine AL_formResidual

20
code/numerics.f90
View File

@ -92,7 +92,9 @@ module numerics
err_f_tol = 1e-6_pReal, & err_f_tol = 1e-6_pReal, &
err_p_tol = 1e-5_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 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 :: & 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 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 :: & character(len=64), protected, public :: &
@ -305,9 +307,14 @@ subroutine numerics_init
err_f_tol = IO_floatValue(line,positions,2_pInt) err_f_tol = IO_floatValue(line,positions,2_pInt)
case ('err_p_tol') case ('err_p_tol')
err_p_tol = IO_floatValue(line,positions,2_pInt) 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 #endif
#ifndef PETSc #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) call IO_warning(41_pInt,ext_msg=tag)
#endif #endif
#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 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', & 'itmax', 'itmin','memory_efficient','fftw_timelimit','fftw_plan_mode','myspectralsolver', &
'rotation_tol','divergence_correction','update_gamma','petsc_options','myfilter', & '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) call IO_warning(40_pInt,ext_msg=tag)
#endif #endif
case default ! found unknown keyword case default ! found unknown keyword
@ -425,6 +432,8 @@ subroutine numerics_init
#ifdef PETSc #ifdef PETSc
write(6,'(a24,1x,es8.1)') ' err_f_tol: ',err_f_tol 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)') ' 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)') ' myspectralsolver: ',trim(myspectralsolver)
write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_options) write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_options)
#endif #endif
@ -484,6 +493,11 @@ subroutine numerics_init
#ifdef PETSc #ifdef PETSc
if (err_f_tol <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_f_tol') 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 (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
#endif #endif
if (fixedSeed <= 0_pInt) & if (fixedSeed <= 0_pInt) &