made restart working for all solvers and added corresponding test

corrected error handling for FFT based function in case of odd numbers etc.
This commit is contained in:
Martin Diehl 2012-12-14 17:30:22 +00:00
parent 7b87987751
commit b1f28ba5e7
9 changed files with 109 additions and 66 deletions

View File

@ -440,7 +440,7 @@ program DAMASK_spectral
!--------------------------------------------------------------------------------------------------
! general initialization of fftw (see manual on fftw.org for more details)
if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt) ! check for correct precision in C
if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=0_pInt) ! check for correct precision in C
!$ if(DAMASK_NumThreadsInt > 0_pInt) then
!$ ierr = fftw_init_threads()
!$ if (ierr == 0_pInt) call IO_error(error_ID = 809_pInt)

View File

@ -93,6 +93,7 @@ subroutine AL_init()
use mesh, only: &
res, &
geomdim, &
wgt, &
mesh_NcpElems, &
mesh_ipCoordinates
use math, only: &
@ -103,7 +104,12 @@ subroutine AL_init()
#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) :: &
temp33_Real = 0.0_pReal
real(pReal), dimension(3,3,3,3) :: &
temp3333_Real = 0.0_pReal, &
temp3333_Real2 = 0.0_pReal
PetscErrorCode :: ierr
PetscObject :: dummy
PetscScalar, pointer, dimension(:,:,:,:) :: xx_psc, F, F_lambda
@ -153,50 +159,55 @@ subroutine AL_init()
elseif (restartInc > 1_pInt) then ! using old values from file
if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',&
restartInc - 1_pInt,' from file'
call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',&
trim(getSolverJobName()),size(F_lastInc))
flush(6)
call IO_read_jobBinaryFile(777,'F',&
trim(getSolverJobName()),size(F))
read (777,rec=1) F
close (777)
call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',&
call IO_read_jobBinaryFile(777,'F_lastInc',&
trim(getSolverJobName()),size(F_lastInc))
read (777,rec=1) F_lastInc
close (777)
call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda',&
trim(getSolverJobName()),size(F_lambda_lastInc))
F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F
F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc
call IO_read_jobBinaryFile(777,'F_lambda',&
trim(getSolverJobName()),size(F_lambda))
read (777,rec=1) F_lambda
close (777)
call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda_lastInc',&
trim(getSolverJobName()),size(F_lambda_lastInc))
call IO_read_jobBinaryFile(777,'F_lambda_lastInc',&
trim(getSolverJobName()),size(F_lambda_lastInc))
read (777,rec=1) F_lambda_lastInc
close (777)
call IO_read_jobBinaryFile(777,'F_aim',trim(getSolverJobName()),size(F_aim))
read (777,rec=1) F_aim
call IO_read_jobBinaryFile(777,'C_lastInc',trim(getSolverJobName()),size(C_lastInc))
read (777,rec=1) C_lastInc
close (777)
call IO_read_jobBinaryFile(777,'F_aim_lastInc',trim(getSolverJobName()),size(F_aim_lastInc))
read (777,rec=1) F_aim_lastInc
call IO_read_jobBinaryFile(777,'C',trim(getSolverJobName()),size(C))
read (777,rec=1) C
close (777)
call IO_read_jobBinaryFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot))
read (777,rec=1) f_aimDot
close (777)
call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(temp3333_Real))
read (777,rec=1) temp3333_Real
close (777)
endif
mesh_ipCoordinates = 0.0_pReal !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,C,P_av,.false.,math_I3)
call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,temp3333_Real2,&
temp33_Real,.false.,math_I3)
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr)
CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! reference stiffness
if (restartInc == 1_pInt) then
call IO_write_jobBinaryFile(777,'C_ref',size(C))
write (777,rec=1) C
close(777)
elseif (restartInc > 1_pInt) then
call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(C))
read (777,rec=1) C
close (777)
endif
if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness
temp3333_Real = temp3333_Real2
C = temp3333_Real2
endif
call Utilities_updateGamma(C,.True.)
C_scale = C
S_scale = math_invSym3333(C)
call Utilities_updateGamma(temp3333_Real,.True.)
C_scale = temp3333_Real
S_scale = math_invSym3333(temp3333_Real)
end subroutine AL_init
@ -239,35 +250,46 @@ type(tSolutionState) function &
type(tBoundaryCondition), intent(in) :: P_BC,F_BC
character(len=*), intent(in) :: incInfoIn
real(pReal), dimension(3,3), intent(in) :: rotation_BC
real(pReal), dimension(3,3) ,save :: F_aimDot
real(pReal), dimension(3,3) :: F_aim_lab
!--------------------------------------------------------------------------------------------------
! loop variables, convergence etc.
real(pReal), dimension(3,3) :: temp33_Real
!--------------------------------------------------------------------------------------------------
! PETSc Data
PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F, F_lambda
PetscErrorCode :: ierr
SNESConvergedReason ::reason
incInfo = incInfoIn
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr)
F => xx_psc(0:8,:,:,:)
F_lambda => xx_psc(9:17,:,:,:)
!--------------------------------------------------------------------------------------------------
! restart information for spectral solver
incInfo = incInfoIn
if (restartWrite) then
write(6,'(a)') 'writing converged results for restart'
call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F_lastInc))
write (777,rec=1) F_LastInc
call IO_write_jobBinaryFile(777,'F',size(F)) ! writing deformation gradient field to file
write (777,rec=1) F
close (777)
call IO_write_jobBinaryFile(777,'F_lastInc',size(F_lastInc)) ! writing F_lastInc field to file
write (777,rec=1) F_lastInc
close (777)
call IO_write_jobBinaryFile(777,'F_lambda',size(F_lambda)) ! writing deformation gradient field to file
write (777,rec=1) F_lambda
close (777)
call IO_write_jobBinaryFile(777,'F_lambda_lastInc',size(F_lambda_lastInc)) ! writing F_lastInc field to file
write (777,rec=1) F_lambda_lastInc
close (777)
call IO_write_jobBinaryFile(777,'C',size(C))
write (777,rec=1) C
close(777)
call IO_write_jobBinaryFile(777,'C_lastInc',size(C_lastInc))
write (777,rec=1) C_lastInc
close(777)
call IO_write_jobBinaryFile(777,'F_aimDot',size(F_aimDot))
write (777,rec=1) F_aimDot
close(777)
endif
AL_solution%converged =.false.
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr)
F => xx_psc(0:8,:,:,:)
F_lambda => xx_psc(9:17,:,:,:)
if ( cutBack) then
F_aim = F_aim_lastInc

View File

@ -97,11 +97,11 @@ subroutine basic_init()
elseif (restartInc > 1_pInt) then ! using old values from file
if (debugRestart) write(6,'(a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'Reading values of increment', restartInc - 1_pInt, 'from file'
call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',&
call IO_read_jobBinaryFile(777,'F',&
trim(getSolverJobName()),size(F))
read (777,rec=1) F
close (777)
call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',&
call IO_read_jobBinaryFile(777,'F_lastInc',&
trim(getSolverJobName()),size(F_lastInc))
read (777,rec=1) F_lastInc
close (777)
@ -216,10 +216,10 @@ type(tSolutionState) function &
! write restart information for spectral solver
if (restartWrite) then
write(6,'(a)') 'writing converged results for restart'
call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F)) ! writing deformation gradient field to file
call IO_write_jobBinaryFile(777,'F',size(F)) ! writing deformation gradient field to file
write (777,rec=1) F
close (777)
call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',size(F_lastInc)) ! writing F_lastInc field to file
call IO_write_jobBinaryFile(777,'F_lastInc',size(F_lastInc)) ! writing F_lastInc field to file
write (777,rec=1) F_lastInc
close (777)
call IO_write_jobBinaryFile(777,'C',size(C))

View File

@ -150,11 +150,11 @@ subroutine basicPETSc_init()
if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',&
restartInc - 1_pInt,' from file'
flush(6)
call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',&
call IO_read_jobBinaryFile(777,'F',&
trim(getSolverJobName()),size(F))
read (777,rec=1) F
close (777)
call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',&
call IO_read_jobBinaryFile(777,'F_lastInc',&
trim(getSolverJobName()),size(F_lastInc))
read (777,rec=1) F_lastInc
close (777)
@ -167,8 +167,8 @@ subroutine basicPETSc_init()
call IO_read_jobBinaryFile(777,'C',trim(getSolverJobName()),size(C))
read (777,rec=1) C
close (777)
call IO_read_jobBinaryFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot))
read (777,rec=1) f_aimDot
call IO_read_jobBinaryFile(777,'F_aimDot',trim(getSolverJobName()),size(F_aimDot))
read (777,rec=1) F_aimDot
close (777)
call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(temp3333_Real))
read (777,rec=1) temp3333_Real
@ -243,14 +243,15 @@ type(tSolutionState) function &
incInfo = incInfoIn
call DMDAVecGetArrayF90(da,solution_vec,F,ierr)
!--------------------------------------------------------------------------------------------------
! write restart information for spectral solver
if (restartWrite) then
write(6,'(a)') 'writing converged results for restart'
call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F)) ! writing deformation gradient field to file
call IO_write_jobBinaryFile(777,'F',size(F)) ! writing deformation gradient field to file
write (777,rec=1) F
close (777)
call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',size(F_lastInc)) ! writing F_lastInc field to file
call IO_write_jobBinaryFile(777,'F_lastInc',size(F_lastInc)) ! writing F_lastInc field to file
write (777,rec=1) F_lastInc
close (777)
call IO_write_jobBinaryFile(777,'C',size(C))
@ -259,8 +260,8 @@ type(tSolutionState) function &
call IO_write_jobBinaryFile(777,'C_lastInc',size(C_lastInc))
write (777,rec=1) C_lastInc
close(777)
call IO_write_jobBinaryFile(777,'F_aimDot',size(f_aimDot))
write (777,rec=1) f_aimDot
call IO_write_jobBinaryFile(777,'F_aimDot',size(F_aimDot))
write (777,rec=1) F_aimDot
close(777)
endif
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomdim,reshape(F,[3,3,res(1),res(2),res(3)])),&

View File

@ -145,7 +145,7 @@ subroutine utilities_init()
!--------------------------------------------------------------------------------------------------
! general initialization of FFTW (see manual on fftw.org for more details)
if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt) ! check for correct precision in C
if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(0_pInt,ext_msg='Fortran to C') ! check for correct precision in C
!$ if(DAMASK_NumThreadsInt > 0_pInt) then
!$ i = fftw_init_threads() ! returns 0 in case of problem
!$ if (i == 0_pInt) call IO_error(error_ID = 809_pInt)
@ -588,7 +588,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
c_reduced(k,j) = temp99_Real(n,m)
endif; enddo; endif; enddo
call math_invert(size_reduced, c_reduced, s_reduced, errmatinv) ! invert reduced stiffness
if(errmatinv) call IO_error(error_ID=400_pInt)
if(errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance')
temp99_Real = 0.0_pReal ! fill up compliance with zeros
k = 0_pInt
do n = 1_pInt,9_pInt
@ -615,7 +615,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
write(6,trim(formatString),advance='no') 'C * S', transpose(matmul(c_reduced,s_reduced))
write(6,trim(formatString),advance='no') 'S', transpose(s_reduced)
endif
if(errmatinv) call IO_error(error_ID=400_pInt)
if(errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance')
deallocate(c_reduced)
deallocate(s_reduced)
deallocate(sTimesC)

View File

@ -70,3 +70,4 @@ divergence_correction 0 # Use dimension-independent divergence cr
myspectralsolver basic # Type of spectral solver (basic: basic, basicPETSc: basic with PETSc, AL: augmented Lagrange)
myfilter none # Type of filtering method to mitigate Gibb's phenomenon (none, cosine, ...)
petsc_options -snes_type ngmres -snes_ngmres_anderson -snes_view # PetSc solver options
regridMode 0 # 0: no regrid; 1: regrid if DAMASK doesn't converge; 2: regrid if DAMASK or BVP Solver doesn't converge

View File

@ -2869,16 +2869,22 @@ function math_curlFFT(geomdim,field)
res = [size(field,1),size(field,2),size(field,3)]
vec_tens = size(field,4)
if (vec_tens /= 1_pInt .and. vec_tens /= 3_pInt) call IO_error(0_pInt, &
ext_msg = 'Curl: invalid data dimension')
if (iand(debug_level(debug_math),debug_levelBasic) /= 0_pInt) then
if (vec_tens == 1_pInt) print*, 'Calculating curl of vector field'
if (vec_tens == 3_pInt) print*, 'Calculating curl of tensor field'
print '(a,3(e12.5))', ' Dimension: ', geomdim
print '(a,3(i5))', ' Resolution:', res
endif
if (vec_tens /= 1_pInt .and. vec_tens /= 3_pInt) &
call IO_error(0_pInt, ext_msg = 'Invalid data type in math_curlFFT')
if ((mod(res(3),2_pInt)/=0_pInt .and. res(3) /= 1_pInt) .or. &
mod(res(2),2_pInt)/=0_pInt .or. &
mod(res(1),2_pInt)/=0_pInt) &
call IO_error(0_pInt,ext_msg='Resolution in math_curlFFT')
if (pReal /= C_DOUBLE .or. pInt /= C_INT) &
call IO_error(0_pInt,ext_msg='Fortran to C in math_curlFFT')
wgt = 1.0_pReal/real(res(1)*res(2)*res(3),pReal)
res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
@ -2995,8 +3001,6 @@ function math_divergenceFFT(geomdim,field)
res = [size(field,1),size(field,2),size(field,3)]
vec_tens = size(field,4)
if (vec_tens /= 1_pInt .and. vec_tens /= 3_pInt) call IO_error(0_pInt, &
ext_msg = 'Divergence: invalid data dimension')
if (iand(debug_level(debug_math),debug_levelBasic) /= 0_pInt) then
if (vec_tens == 1_pInt) print*, 'Calculating divergence of vector field'
@ -3004,7 +3008,15 @@ function math_divergenceFFT(geomdim,field)
print '(a,3(e12.5))', ' Dimension: ', geomdim
print '(a,3(i5))', ' Resolution:', res
endif
if (vec_tens /= 1_pInt .and. vec_tens /= 3_pInt) &
call IO_error(0_pInt, ext_msg = 'Invalid data type in math_divergenceFFT')
if ((mod(res(3),2_pInt)/=0_pInt .and. res(3) /= 1_pInt) .or. &
mod(res(2),2_pInt)/=0_pInt .or. &
mod(res(1),2_pInt)/=0_pInt) &
call IO_error(0_pInt,ext_msg='Resolution in math_divergenceFFT')
if (pReal /= C_DOUBLE .or. pInt /= C_INT) &
call IO_error(0_pInt,ext_msg='Fortran to C in math_divergenceFFT')
res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
wgt = 1.0_pReal/real(res(1)*res(2)*res(3),pReal)
@ -3081,6 +3093,7 @@ end function math_divergenceFFT
!--------------------------------------------------------------------------------------------------
function math_divergenceFDM(geomdim,order,field)
use IO, only: IO_error
use debug, only: debug_math, &
debug_level, &
debug_levelBasic
@ -3111,6 +3124,9 @@ function math_divergenceFDM(geomdim,order,field)
print '(a,3(i5))', ' Resolution:', res
endif
if (vec_tens /= 1_pInt .and. vec_tens /= 3_pInt) &
call IO_error(0_pInt, ext_msg = 'Invalid data type in math_divergenceFDM')
math_divergenceFDM = 0.0_pReal
do k = 0_pInt, res(3)-1_pInt; do j = 0_pInt, res(2)-1_pInt; do i = 0_pInt, res(1)-1_pInt
do m = 1_pInt, order + 1_pInt

View File

@ -382,7 +382,6 @@ subroutine mesh_init(ip,element)
use DAMASK_interface
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use IO, only: &
IO_error, &
#ifdef Abaqus
IO_abaqus_hasNoPart, &
#endif
@ -1045,8 +1044,6 @@ end subroutine mesh_spectral_count_cpSizes
!--------------------------------------------------------------------------------------------------
subroutine mesh_spectral_build_nodes()
use IO, only: &
IO_error
use numerics, only: numerics_unitlength
implicit none
@ -2018,8 +2015,12 @@ function mesh_deformedCoordsFFT(geomdim,F,scalingIn,FavgIn)
res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
step = geomdim/real(res, pReal)
if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt)
if ((mod(res(3),2_pInt)/=0_pInt .and. res(3) /= 1_pInt) .or. &
mod(res(2),2_pInt)/=0_pInt .or. &
mod(res(1),2_pInt)/=0_pInt) &
call IO_error(0_pInt,ext_msg='Resolution in mesh_deformedCoordsFFT')
if (pReal /= C_DOUBLE .or. pInt /= C_INT) &
call IO_error(0_pInt,ext_msg='Fortran to C in mesh_deformedCoordsFFT')
call fftw_set_timelimit(fftw_timelimit)
defgrad_fftw = fftw_alloc_complex(int(res1_red *res(2)*res(3)*9_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8)
call c_f_pointer(defgrad_fftw, F_real, [res(1)+2_pInt,res(2),res(3),3_pInt,3_pInt])

View File

@ -339,8 +339,10 @@ subroutine numerics_init
end select
#endif
#ifdef PETSc
call PetscOptionsInsertString(petsc_options,ierr)
write(6,'(a)') ' Initializing PETSc'
call PetscOptionsClear(ierr)
CHKERRQ(ierr)
call PetscOptionsInsertString(petsc_options,ierr)
CHKERRQ(ierr)
#endif