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

corrected reporting of basic PETSc variant and unified reporting of all spectral solvers 

improved warning and error in IO, now able to report correctly elements up to 9 digits
This commit is contained in:
Martin Diehl 2013-01-09 22:19:32 +00:00
parent b8106429f6
commit 0d5e91ac87
8 changed files with 273 additions and 1335 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1104
code/DAMASK_spectral.f90
View File

File diff suppressed because it is too large Load Diff

20
code/DAMASK_spectral_driver.f90
View File

@ -331,7 +331,8 @@ program DAMASK_spectral_Driver
open(newunit=statUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& open(newunit=statUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
'.sta',form='FORMATTED',status='REPLACE') '.sta',form='FORMATTED',status='REPLACE')
write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file
if (debugGeneral) write(6,'(a)') 'Header of result file written out' if (debugGeneral) write(6,'(/,a)') ' header of result file written out'
flush(6)
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! loopping over loadcases ! loopping over loadcases
@ -379,8 +380,8 @@ program DAMASK_spectral_Driver
stepFraction = stepFraction + 1_pInt stepFraction = stepFraction + 1_pInt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new increment ! report begin of new increment
write(6,'(1/,a)') '###########################################################################' write(6,'(/,a)') ' ###########################################################################'
write(6,'(a,es12.5'//& write(6,'(1x,a,es12.5'//&
',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//& ',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//&
',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//& ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//&
',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') & ',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') &
@ -388,9 +389,10 @@ program DAMASK_spectral_Driver
's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,& 's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,&
'-', stepFraction, '/', subStepFactor**cutBackLevel,& '-', stepFraction, '/', subStepFactor**cutBackLevel,&
' of load case ', currentLoadCase,'/',size(loadCases) ' of load case ', currentLoadCase,'/',size(loadCases)
flush(6)
write(incInfo,'(a,'//IO_intOut(totalIncsCounter)//',a,'//IO_intOut(sum(loadCases(:)%incs))//& write(incInfo,'(a,'//IO_intOut(totalIncsCounter)//',a,'//IO_intOut(sum(loadCases(:)%incs))//&
',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') & ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') &
'Inc. ',totalIncsCounter,'/',sum(loadCases(:)%incs),& 'Increment ',totalIncsCounter,'/',sum(loadCases(:)%incs),&
'-',stepFraction, '/', subStepFactor**cutBackLevel '-',stepFraction, '/', subStepFactor**cutBackLevel
select case(myspectralsolver) select case(myspectralsolver)
@ -449,13 +451,14 @@ program DAMASK_spectral_Driver
cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half number of subincs next inc cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half number of subincs next inc
if(solres%converged) then ! report converged inc if(solres%converged) then ! report converged inc
convergedCounter = convergedCounter + 1_pInt convergedCounter = convergedCounter + 1_pInt
write(6,'(A,'//IO_intOut(totalIncsCounter)//',A)') & write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',A)') &
' increment ', totalIncsCounter, ' converged' ' increment ', totalIncsCounter, ' converged'
else else
write(6,'(A,'//IO_intOut(totalIncsCounter)//',A)') & ! report non-converged inc write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',A)') & ! report non-converged inc
' increment ', totalIncsCounter, ' NOT converged' ' increment ', totalIncsCounter, ' NOT converged'
notConvergedCounter = notConvergedCounter + 1_pInt notConvergedCounter = notConvergedCounter + 1_pInt
endif endif
flush(6)
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency
write(6,'(1/,a)') ' ... writing results to file ......................................' write(6,'(1/,a)') ' ... writing results to file ......................................'
@ -489,9 +492,8 @@ program DAMASK_spectral_Driver
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! done report summary ! done report summary
write(6,'(a)') '' write(6,'(/,a)') ' ##################################################################'
write(6,'(a)') '##################################################################' write(6,'(1x,i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', &
write(6,'(i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', &
notConvergedCounter + convergedCounter, ' (', & notConvergedCounter + convergedCounter, ' (', &
real(convergedCounter, pReal)/& real(convergedCounter, pReal)/&
real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, & real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, &

75
code/DAMASK_spectral_solverAL.f90
View File

@ -68,7 +68,7 @@ module DAMASK_spectral_solverAL
err_f, & !< difference between compatible and incompatible deformation gradient err_f, & !< difference between compatible and incompatible deformation gradient
err_p !< difference of stress resulting from compatible and incompatible F err_p !< difference of stress resulting from compatible and incompatible F
logical, private :: ForwardData logical, private :: ForwardData
integer(pInt) :: reportIter = 0_pInt integer(pInt), private :: reportIter = 0_pInt
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -120,7 +120,6 @@ subroutine AL_init(temperature)
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverAL init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverAL init -+>>>'
write(6,'(a)') ' $Id: DAMASK_spectral_SolverAL.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $' write(6,'(a)') ' $Id: DAMASK_spectral_SolverAL.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $'
#include "compilation_info.f90" #include "compilation_info.f90"
write(6,'(a)') ''
allocate (F_lastInc (3,3, res(1), res(2),res(3)), source = 0.0_pReal) allocate (F_lastInc (3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (Fdot (3,3, res(1), res(2),res(3)), source = 0.0_pReal) allocate (Fdot (3,3, res(1), res(2),res(3)), source = 0.0_pReal)
@ -154,7 +153,7 @@ subroutine AL_init(temperature)
F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) F = reshape(F_lastInc,[9,res(1),res(2),res(3)])
F_lambda = F F_lambda = F
elseif (restartInc > 1_pInt) then ! using old values from file elseif (restartInc > 1_pInt) then ! using old values from file
if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',& if (debugRestart) write(6,'(/,a,i6,a)') ' reading values of increment ',&
restartInc - 1_pInt,' from file' restartInc - 1_pInt,' from file'
flush(6) flush(6)
call IO_read_jobBinaryFile(777,'F',& call IO_read_jobBinaryFile(777,'F',&
@ -269,7 +268,8 @@ type(tSolutionState) function &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! restart information for spectral solver ! restart information for spectral solver
if (restartWrite) then if (restartWrite) then
write(6,'(a)') 'writing converged results for restart' write(6,'(/,a)') ' writing converged results for restart'
flush(6)
call IO_write_jobBinaryFile(777,'F',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 write (777,rec=1) F
close (777) close (777)
@ -406,10 +406,8 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
residual_F => f_scal(:,:,1,:,:,:) residual_F => f_scal(:,:,1,:,:,:)
residual_F_lambda => f_scal(:,:,2,:,:,:) residual_F_lambda => f_scal(:,:,2,:,:,:)
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr) call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
CHKERRQ(ierr) call SNESGetIterationNumber(snes,iter,ierr); CHKERRQ(ierr)
call SNESGetIterationNumber(snes,iter,ierr)
CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new iteration ! report begin of new iteration
@ -418,13 +416,14 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
reportIter = 0_pInt reportIter = 0_pInt
endif endif
if (callNo == 0 .or. mod(callNo,2) == 1_pInt) then if (callNo == 0 .or. mod(callNo,2) == 1_pInt) then
write(6,'(/,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), &
' @ Iter. ', itmin, '≤',reportIter, '≤', itmax ' @ Iteration ', itmin, '≤',reportIter, '≤', itmax
if (debugRotation) & if (debugRotation) &
write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim (lab)=', & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', &
math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC))
write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =', & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', &
math_transpose33(F_aim) math_transpose33(F_aim)
flush(6)
reportIter = reportIter + 1_pInt reportIter = reportIter + 1_pInt
endif endif
callNo = callNo +1_pInt callNo = callNo +1_pInt
@ -473,7 +472,6 @@ end subroutine AL_formResidual
!> @brief convergence check !> @brief convergence check
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine AL_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ierr) subroutine AL_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ierr)
use numerics, only: & use numerics, only: &
itmax, & itmax, &
itmin, & itmin, &
@ -483,58 +481,59 @@ subroutine AL_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ierr)
err_stress_tolabs err_stress_tolabs
implicit none implicit none
SNES :: snes_local SNES :: snes_local
PetscInt :: it PetscInt :: it
PetscReal :: xnorm, snorm, fnorm PetscReal :: &
xnorm, &
snorm, &
fnorm
SNESConvergedReason :: reason SNESConvergedReason :: reason
PetscObject :: dummy PetscObject :: dummy
PetscErrorCode ::ierr PetscErrorCode ::ierr
logical :: Converged logical :: Converged
real(pReal) :: err_stress_tol
err_stress_tol = min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs)
Converged = (it > itmin .and. & Converged = (it > itmin .and. &
all([ err_f/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_f_tol, & all([ err_f/sqrt(sum((F_aim-math_I3)**2.0_pReal))/err_f_tol, &
err_p/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_p_tol, & err_p/sqrt(sum((F_aim-math_I3)**2.0_pReal))/err_p_tol, &
err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs)] < 1.0_pReal)) err_stress/err_stress_tol] < 1.0_pReal))
if (Converged) then if (Converged) then
reason = 1 reason = 1
elseif (it > itmax) then elseif (it >= itmax) then
reason = -1 reason = -1
else else
reason = 0 reason = 0
endif endif
write(6,'(1/,a)') ' ... reporting ....................................................'
write(6,'(a,f12.7,1x,1a,1x,es9.3)') 'error stress BC = ', & write(6,'(/,a,f8.2,a,es11.5,a,es11.4,a)') ' mismatch F = ', &
err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs),& err_f/sqrt(sum((F_aim-math_I3)**2.0_pReal))/err_f_tol, &
'@',min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs) ' (',err_f/sqrt(sum((F_aim-math_I3)**2.0_pReal)),' -, tol =',err_f_tol,')'
write(6,'(a,f12.7,1x,1a,1x,es9.3)') 'error F = ',& write(6,'(a,f8.2,a,es11.5,a,es11.4,a)') ' mismatch P = ', &
err_f/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_f_tol,& err_p/sqrt(sum((F_aim-math_I3)**2.0_pReal))/err_p_tol, &
'@',err_f_tol ' (',err_p/sqrt(sum((F_aim-math_I3)**2.0_pReal)),' -, tol =',err_p_tol,')'
write(6,'(a,f12.7,1x,1a,1x,es9.3)') 'error P = ', & write(6,'(a,f8.2,a,es11.5,a,es11.4,a)') ' error stress BC = ', &
err_p/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_p_tol,& err_stress/err_stress_tol, ' (',err_stress, ' Pa, tol =',err_stress_tol,')'
'@',err_p_tol
write(6,'(/,a)') ' ==========================================================================' write(6,'(/,a)') ' =========================================================================='
flush(6)
end subroutine AL_converged end subroutine AL_converged
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief destroy routine !> @brief destroy routine
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine AL_destroy() subroutine AL_destroy()
use DAMASK_spectral_Utilities, only: & use DAMASK_spectral_Utilities, only: &
Utilities_destroy Utilities_destroy
implicit none implicit none
PetscErrorCode :: ierr PetscErrorCode :: ierr
call VecDestroy(solution_vec,ierr) call VecDestroy(solution_vec,ierr); CHKERRQ(ierr)
CHKERRQ(ierr) call SNESDestroy(snes,ierr); CHKERRQ(ierr)
call SNESDestroy(snes,ierr) call DMDestroy(da,ierr); CHKERRQ(ierr)
CHKERRQ(ierr) call PetscFinalize(ierr); CHKERRQ(ierr)
call DMDestroy(da,ierr)
CHKERRQ(ierr)
call PetscFinalize(ierr)
CHKERRQ(ierr)
call Utilities_destroy() call Utilities_destroy()
end subroutine AL_destroy end subroutine AL_destroy

28
code/DAMASK_spectral_solverBasic.f90
View File

@ -82,7 +82,6 @@ subroutine basic_init(temperature)
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>'
write(6,'(a)') ' $Id$' write(6,'(a)') ' $Id$'
#include "compilation_info.f90" #include "compilation_info.f90"
write(6,'(a)') ''
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate global fields ! allocate global fields
@ -96,8 +95,9 @@ subroutine basic_init(temperature)
F = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity F = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity
F_lastInc = F F_lastInc = F
elseif (restartInc > 1_pInt) then ! using old values from file elseif (restartInc > 1_pInt) then ! using old values from file
if (debugRestart) write(6,'(a,'//IO_intOut(restartInc-1_pInt)//',a)') & if (debugRestart) write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'Reading values of increment', restartInc - 1_pInt, 'from file' 'reading values of increment', restartInc - 1_pInt, 'from file'
flush(6)
call IO_read_jobBinaryFile(777,'F',& call IO_read_jobBinaryFile(777,'F',&
trim(getSolverJobName()),size(F)) trim(getSolverJobName()),size(F))
read (777,rec=1) F read (777,rec=1) F
@ -216,7 +216,8 @@ type(tSolutionState) function &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! write restart information for spectral solver ! write restart information for spectral solver
if (restartWrite) then if (restartWrite) then
write(6,'(a)') 'writing converged results for restart' write(6,'(/,a)') ' writing converged results for restart'
flush(6)
call IO_write_jobBinaryFile(777,'F',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 write (777,rec=1) F
close (777) close (777)
@ -277,13 +278,14 @@ type(tSolutionState) function &
iter = iter + 1_pInt iter = iter + 1_pInt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new iteration ! report begin of new iteration
write(6,'(/,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), &
' @ Iter. ', itmin, '≤',iter, '≤', itmax ' @ Iteration ', itmin, '≤',iter, '≤', itmax
if (debugRotation) & if (debugRotation) &
write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim (lab)=', & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab)=', &
math_transpose33(math_rotate_backward33(F_aim,rotation_BC)) math_transpose33(math_rotate_backward33(F_aim,rotation_BC))
write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =', & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', &
math_transpose33(F_aim) math_transpose33(F_aim)
flush(6)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
F_aim_lastIter = F_aim F_aim_lastIter = F_aim
@ -311,6 +313,7 @@ type(tSolutionState) function &
F = F - reshape(field_real(1:res(1),1:res(2),1:res(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 F = F - reshape(field_real(1:res(1),1:res(2),1:res(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) basic_solution%converged = basic_Converged(err_div,P_av,err_stress,P_av)
write(6,'(/,a)') ' ==========================================================================' write(6,'(/,a)') ' =========================================================================='
flush(6)
if ((basic_solution%converged .and. iter >= itmin) .or. basic_solution%termIll) then if ((basic_solution%converged .and. iter >= itmin) .or. basic_solution%termIll) then
basic_solution%iterationsNeeded = iter basic_solution%iterationsNeeded = iter
exit exit
@ -353,10 +356,11 @@ logical function basic_Converged(err_div,pAvgDiv,err_stress,pAvgStress)
basic_Converged = all([ err_div/pAvgDivL2/err_div_tol,& basic_Converged = all([ err_div/pAvgDivL2/err_div_tol,&
err_stress/err_stress_tol ] < 1.0_pReal) err_stress/err_stress_tol ] < 1.0_pReal)
write(6,'(a,f6.2,a,es11.4,a)') 'error divergence = ', err_div/pAvgDivL2/err_div_tol,& write(6,'(/,a,f8.2,a,es11.5,a,es11.4,a)') ' error divergence = ', &
' (',err_div,' N/m³)' err_div/pAvgDivL2/err_div_tol, ' (',err_div/pAvgDivL2,' / m, tol =',err_div_tol,')'
write(6,'(a,f6.2,a,es11.4,a)') 'error stress = ', err_stress/err_stress_tol, & write(6,'(a,f8.2,a,es11.5,a,es11.4,a)') ' error stress BC = ', &
' (',err_stress,' Pa)' err_stress/err_stress_tol, ' (',err_stress, ' Pa , tol =',err_stress_tol,')'
flush(6)
end function basic_Converged end function basic_Converged

118
code/DAMASK_spectral_solverBasicPETSc.f90
View File

@ -58,12 +58,15 @@ module DAMASK_spectral_SolverBasicPETSc
real(pReal), private :: err_stress, err_div real(pReal), private :: err_stress, err_div
logical, private :: ForwardData logical, private :: ForwardData
integer(pInt), private :: reportIter = 0_pInt
real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal
public :: basicPETSc_init, & public :: &
basicPETSc_init, &
basicPETSc_solution ,& basicPETSc_solution ,&
basicPETSc_destroy basicPETSc_destroy
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -114,7 +117,6 @@ subroutine basicPETSc_init(temperature)
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasicPETSc init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasicPETSc init -+>>>'
write(6,'(a)') ' $Id: DAMASK_spectral_SolverBasicPETSC.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $' write(6,'(a)') ' $Id: DAMASK_spectral_SolverBasicPETSC.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $'
#include "compilation_info.f90" #include "compilation_info.f90"
write(6,'(a)') ''
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate global fields ! allocate global fields
@ -123,33 +125,27 @@ subroutine basicPETSc_init(temperature)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,snes,ierr) call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
CHKERRQ(ierr)
call DMDACreate3d(PETSC_COMM_WORLD, & call DMDACreate3d(PETSC_COMM_WORLD, &
DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, & DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, &
DMDA_STENCIL_BOX,res(1),res(2),res(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, & DMDA_STENCIL_BOX,res(1),res(2),res(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, &
9,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr) 9,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr); CHKERRQ(ierr)
CHKERRQ(ierr) call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr)
call DMCreateGlobalVector(da,solution_vec,ierr) call DMDASetLocalFunction(da,BasicPETSC_formResidual,ierr); CHKERRQ(ierr)
CHKERRQ(ierr) call SNESSetDM(snes,da,ierr); CHKERRQ(ierr)
call DMDASetLocalFunction(da,BasicPETSC_formResidual,ierr)
CHKERRQ(ierr)
call SNESSetDM(snes,da,ierr)
CHKERRQ(ierr)
call SNESSetConvergenceTest(snes,BasicPETSC_converged,dummy,PETSC_NULL_FUNCTION,ierr) call SNESSetConvergenceTest(snes,BasicPETSC_converged,dummy,PETSC_NULL_FUNCTION,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetFromOptions(snes,ierr) call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr)
CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
call DMDAVecGetArrayF90(da,solution_vec,F,ierr) ! get the data out of PETSc to work with call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with
CHKERRQ(ierr)
if (restartInc == 1_pInt) then ! no deformation (no restart) if (restartInc == 1_pInt) then ! no deformation (no restart)
F_lastInc = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity F_lastInc = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity
F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) F = reshape(F_lastInc,[9,res(1),res(2),res(3)])
elseif (restartInc > 1_pInt) then ! using old values from file elseif (restartInc > 1_pInt) then ! using old values from file
if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',& if (debugRestart) write(6,'(/,a,i6,a)') ' reading values of increment ',&
restartInc - 1_pInt,' from file' restartInc - 1_pInt,' from file'
flush(6) flush(6)
call IO_read_jobBinaryFile(777,'F',& call IO_read_jobBinaryFile(777,'F',&
@ -220,6 +216,7 @@ type(tSolutionState) function &
use FEsolving, only: & use FEsolving, only: &
restartWrite, & restartWrite, &
terminallyIll terminallyIll
implicit none implicit none
#include <finclude/petscdmda.h90> #include <finclude/petscdmda.h90>
#include <finclude/petscsnes.h90> #include <finclude/petscsnes.h90>
@ -240,15 +237,14 @@ type(tSolutionState) function &
PetscScalar, pointer :: F(:,:,:,:) PetscScalar, pointer :: F(:,:,:,:)
PetscErrorCode :: ierr PetscErrorCode :: ierr
SNESConvergedReason :: reason SNESConvergedReason :: reason
incInfo = incInfoIn incInfo = incInfoIn
call DMDAVecGetArrayF90(da,solution_vec,F,ierr) call DMDAVecGetArrayF90(da,solution_vec,F,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! write restart information for spectral solver ! write restart information for spectral solver
if (restartWrite) then if (restartWrite) then
write(6,'(a)') 'writing converged results for restart' write(6,'(/,a)') ' writing converged results for restart'
flush(6)
call IO_write_jobBinaryFile(777,'F',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 write (777,rec=1) F
close (777) close (777)
@ -296,8 +292,7 @@ type(tSolutionState) function &
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot,math_rotate_backward33(F_aim, & F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot,math_rotate_backward33(F_aim, &
rotation_BC)),[9,res(1),res(2),res(3)]) rotation_BC)),[9,res(1),res(2),res(3)])
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator) ! update stiffness (and gamma operator)
@ -320,10 +315,12 @@ type(tSolutionState) function &
BasicPETSC_solution%converged =.false. BasicPETSC_solution%converged =.false.
if (reason > 0 ) then if (reason > 0 ) then
BasicPETSC_solution%converged = .true. BasicPETSC_solution%converged = .true.
BasicPETSC_solution%iterationsNeeded = reportIter
endif endif
end function BasicPETSc_solution end function BasicPETSc_solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forms the AL residual vector !> @brief forms the AL residual vector
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -350,26 +347,36 @@ subroutine BasicPETSC_formResidual(myIn,x_scal,f_scal,dummy,ierr)
implicit none implicit none
DMDALocalInfo, dimension(*) :: myIn DMDALocalInfo, dimension(*) :: myIn
PetscScalar, dimension(3,3,res(1),res(2),res(3)) :: x_scal PetscScalar, dimension(3,3,res(1),res(2),res(3)) :: &
PetscScalar, dimension(3,3,res(1),res(2),res(3)):: f_scal x_scal, &
f_scal
PetscInt :: iter, nfuncs PetscInt :: iter, nfuncs
PetscObject :: dummy PetscObject :: dummy
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer(pInt), save :: callNo = 3_pInt
logical :: report
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr) call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
CHKERRQ(ierr) call SNESGetIterationNumber(snes,iter,ierr); CHKERRQ(ierr)
call SNESGetIterationNumber(snes,iter,ierr)
CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new iteration ! report begin of new iteration
write(6,'(/,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & if (iter == 0 .and. callNo>2) then
' @ Iter. ', itmin, '≤', iter, '≤', itmax callNo = 0_pInt
reportIter = 0_pInt
endif
if (callNo == 0 .or. mod(callNo,2) == 1_pInt) then
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), &
' @ Iteration ', itmin, '≤',reportIter, '≤', itmax
if (debugRotation) & if (debugRotation) &
write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim (lab)=', & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab)=', &
math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC))
write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =', & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', &
math_transpose33(F_aim) math_transpose33(F_aim)
flush(6)
reportIter = reportIter + 1_pInt
endif
callNo = callNo +1_pInt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
@ -396,7 +403,7 @@ subroutine BasicPETSC_formResidual(myIn,x_scal,f_scal,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! constructing residual
f_scal = reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),shape(x_scal),order=[3,4,5,1,2]) f_scal = reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),shape(x_scal),order=[3,4,5,1,2])
write(6,'(/,a)') '=========================================================================='
end subroutine BasicPETSc_formResidual end subroutine BasicPETSc_formResidual
@ -404,47 +411,53 @@ end subroutine BasicPETSc_formResidual
!> @brief convergence check !> @brief convergence check
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine BasicPETSc_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ierr) subroutine BasicPETSc_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ierr)
use numerics, only: & use numerics, only: &
itmax, & itmax, &
itmin, & itmin, &
err_div_tol, & err_div_tol, &
err_stress_tolrel, & err_stress_tolrel, &
err_stress_tolabs err_stress_tolabs
use math, only: & use math, only: &
math_mul33x33, & math_mul33x33, &
math_eigenvalues33, & math_eigenvalues33, &
math_transpose33 math_transpose33
implicit none implicit none
SNES :: snes_local SNES :: snes_local
PetscInt :: it PetscInt :: it
PetscReal :: xnorm, snorm, fnorm PetscReal :: &
xnorm, &
snorm, &
fnorm
SNESConvergedReason :: reason SNESConvergedReason :: reason
PetscObject :: dummy PetscObject :: dummy
PetscErrorCode :: ierr PetscErrorCode :: ierr
logical :: Converged logical :: Converged
real(pReal) :: pAvgDivL2 real(pReal) :: pAvgDivL2, &
err_stress_tol
err_stress_tol =min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs)
pAvgDivL2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(P_av,math_transpose33(P_av))))) pAvgDivL2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(P_av,math_transpose33(P_av)))))
Converged = (it > itmin .and. & Converged = (it >= itmin .and. &
all([ err_div/pAvgDivL2/err_div_tol, & all([ err_div/pAvgDivL2/err_div_tol, &
err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs)] < 1.0_pReal)) err_stress/err_stress_tol] < 1.0_pReal))
if (Converged) then if (Converged) then
reason = 1 reason = 1
elseif (it > itmax) then elseif (it >= itmax) then
reason = -1 reason = -1
else else
reason = 0 reason = 0
endif endif
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,'(a,f6.2,a,es11.4,a)') 'error divergence = ', err_div/pAvgDivL2/err_div_tol,&
' (',err_div/pAvgDivL2,' N/m³)'
write(6,'(a,f6.2,a,es11.4,a)') 'error stress = ', err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs), &
' (',err_stress,' Pa)'
end subroutine BasicPETSc_converged end subroutine BasicPETSc_converged
@ -452,23 +465,18 @@ end subroutine BasicPETSc_converged
!> @brief destroy routine !> @brief destroy routine
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine BasicPETSc_destroy() subroutine BasicPETSc_destroy()
use DAMASK_spectral_Utilities, only: & use DAMASK_spectral_Utilities, only: &
Utilities_destroy Utilities_destroy
implicit none implicit none
PetscErrorCode :: ierr PetscErrorCode :: ierr
call VecDestroy(solution_vec,ierr) call VecDestroy(solution_vec,ierr); CHKERRQ(ierr)
CHKERRQ(ierr) call SNESDestroy(snes,ierr); CHKERRQ(ierr)
call SNESDestroy(snes,ierr) call DMDestroy(da,ierr); CHKERRQ(ierr)
CHKERRQ(ierr) call PetscFinalize(ierr); CHKERRQ(ierr)
call DMDestroy(da,ierr)
CHKERRQ(ierr)
call PetscFinalize(ierr)
CHKERRQ(ierr)
call Utilities_destroy() call Utilities_destroy()
CHKERRQ(ierr)
end subroutine BasicPETSc_destroy end subroutine BasicPETSc_destroy
end module DAMASK_spectral_SolverBasicPETSc end module DAMASK_spectral_SolverBasicPETSc

59
code/DAMASK_spectral_utilities.f90
View File

@ -139,6 +139,7 @@ subroutine utilities_init()
write(6,'(a)') ' $Id$' write(6,'(a)') ' $Id$'
#include "compilation_info.f90" #include "compilation_info.f90"
write(6,'(a)') '' write(6,'(a)') ''
call flush(6)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set debugging parameters ! set debugging parameters
@ -149,14 +150,12 @@ subroutine utilities_init()
debugRotation = iand(debug_level(debug_spectral),debug_spectralRotation) /= 0 debugRotation = iand(debug_level(debug_spectral),debug_spectralRotation) /= 0
#ifdef PETSc #ifdef PETSc
debugPETSc = iand(debug_level(debug_spectral),debug_spectralPETSc) /= 0 debugPETSc = iand(debug_level(debug_spectral),debug_spectralPETSc) /= 0
if(debugPETSc) write(6,'(a)') ' Initializing PETSc with debug options: ', trim(PETScDebug), & if(debugPETSc) write(6,'(/,a)') ' Initializing PETSc with debug options: ', trim(PETScDebug), &
' add more using the PETSc_Options keyword in numerics.config ' ' add more using the PETSc_Options keyword in numerics.config '
call PetscOptionsClear(ierr) flush(6)
CHKERRQ(ierr) call PetscOptionsClear(ierr); CHKERRQ(ierr)
if(debugPETSc) call PetscOptionsInsertString(trim(PETScDebug),ierr) if(debugPETSc) call PetscOptionsInsertString(trim(PETScDebug),ierr); CHKERRQ(ierr)
CHKERRQ(ierr) call PetscOptionsInsertString(trim(petsc_options),ierr); CHKERRQ(ierr)
call PetscOptionsInsertString(trim(petsc_options),ierr)
CHKERRQ(ierr)
#endif #endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocation ! allocation
@ -214,7 +213,8 @@ subroutine utilities_init()
scalarField_fourier,scalarField_real,+1,fftw_planner_flag) ! input, output, backward (1), planner precision scalarField_fourier,scalarField_real,+1,fftw_planner_flag) ! input, output, backward (1), planner precision
endif endif
if (debugGeneral) write(6,'(a)') 'FFTW initialized' if (debugGeneral) write(6,'(/,a)') ' FFTW initialized'
flush(6)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)
@ -268,7 +268,8 @@ subroutine utilities_updateGamma(C,saveReference)
C_ref = C C_ref = C
if (saveReference) then if (saveReference) then
write(6,'(a)') 'writing reference stiffness to file' write(6,'(/,a)') ' writing reference stiffness to file'
flush(6)
call IO_write_jobBinaryFile(777,'C_ref',size(C_ref)) call IO_write_jobBinaryFile(777,'C_ref',size(C_ref))
write (777,rec=1) C_ref write (777,rec=1) C_ref
close(777) close(777)
@ -330,14 +331,16 @@ subroutine utilities_FFTforward(row,column)
! comparing 1 and 3x3 FT results ! comparing 1 and 3x3 FT results
if (debugFFTW) then if (debugFFTW) then
call fftw_execute_dft(plan_scalarField_forth,scalarField_real,scalarField_fourier) call fftw_execute_dft(plan_scalarField_forth,scalarField_real,scalarField_fourier)
write(6,'(a,i1,1x,i1)') 'checking FT results of compontent ', row, column write(6,'(/,a,i1,1x,i1,a)') ' .. checking FT results of compontent ', row, column, ' ..'
write(6,'(a,2(es11.4,1x))') 'max FT relative error = ',& ! print real and imaginary part seperately flush(6)
write(6,'(/,a,2(es11.4,1x))') ' max FT relative error = ',& ! print real and imaginary part seperately
maxval( real((scalarField_fourier(1:res1_red,1:res(2),1:res(3))-& maxval( real((scalarField_fourier(1:res1_red,1:res(2),1:res(3))-&
field_fourier(1:res1_red,1:res(2),1:res(3),row,column))/& field_fourier(1:res1_red,1:res(2),1:res(3),row,column))/&
scalarField_fourier(1:res1_red,1:res(2),1:res(3)))), & scalarField_fourier(1:res1_red,1:res(2),1:res(3)))), &
maxval(aimag((scalarField_fourier(1:res1_red,1:res(2),1:res(3))-& maxval(aimag((scalarField_fourier(1:res1_red,1:res(2),1:res(3))-&
field_fourier(1:res1_red,1:res(2),1:res(3),row,column))/& field_fourier(1:res1_red,1:res(2),1:res(3),row,column))/&
scalarField_fourier(1:res1_red,1:res(2),1:res(3)))) scalarField_fourier(1:res1_red,1:res(2),1:res(3))))
flush(6)
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -395,12 +398,14 @@ subroutine utilities_FFTbackward(row,column)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! comparing 1 and 3x3 inverse FT results ! comparing 1 and 3x3 inverse FT results
if (debugFFTW) then if (debugFFTW) then
write(6,'(a,i1,1x,i1)') 'checking iFT results of compontent ', row, column write(6,'(/,a,i1,1x,i1,a)') ' ... checking iFT results of compontent ', row, column, ' ..'
flush(6)
call fftw_execute_dft(plan_scalarField_back,scalarField_fourier,scalarField_real) call fftw_execute_dft(plan_scalarField_back,scalarField_fourier,scalarField_real)
write(6,'(a,es11.4)') 'max iFT relative error = ',& write(6,'(/,a,es11.4)') ' max iFT relative error = ',&
maxval((real(scalarField_real(1:res(1),1:res(2),1:res(3)))-& maxval((real(scalarField_real(1:res(1),1:res(2),1:res(3)))-&
field_real(1:res(1),1:res(2),1:res(3),row,column))/& field_real(1:res(1),1:res(2),1:res(3),row,column))/&
real(scalarField_real(1:res(1),1:res(2),1:res(3)))) real(scalarField_real(1:res(1),1:res(2),1:res(3))))
flush(6)
endif endif
field_real = field_real * wgt ! normalize the result by number of elements field_real = field_real * wgt ! normalize the result by number of elements
@ -454,6 +459,7 @@ subroutine utilities_fourierConvolution(fieldAim)
l, m, n, o l, m, n, o
write(6,'(/,a)') ' ... doing convolution .....................................................' write(6,'(/,a)') ' ... doing convolution .....................................................'
flush(6)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! to the actual spectral method calculation (mechanical equilibrium) ! to the actual spectral method calculation (mechanical equilibrium)
@ -506,7 +512,7 @@ real(pReal) function utilities_divergenceRMS()
complex(pReal), dimension(3) :: temp3_complex complex(pReal), dimension(3) :: temp3_complex
write(6,'(/,a)') ' ... calculating divergence ................................................' write(6,'(/,a)') ' ... calculating divergence ................................................'
flush(6)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculating RMS divergence criterion in Fourier space ! calculating RMS divergence criterion in Fourier space
utilities_divergenceRMS = 0.0_pReal utilities_divergenceRMS = 0.0_pReal
@ -549,11 +555,12 @@ real(pReal) function utilities_divergenceRMS()
err_real_div_max = sqrt(maxval(sum(divergence_real**2.0_pReal,dim=4))) ! max in real space err_real_div_max = sqrt(maxval(sum(divergence_real**2.0_pReal,dim=4))) ! max in real space
err_div_max = sqrt( err_div_max) ! max in Fourier space err_div_max = sqrt( err_div_max) ! max in Fourier space
write(6,'(1x,a,es11.4)') 'error divergence FT RMS = ',err_div_RMS write(6,'(/,1x,a,es11.4)') 'error divergence FT RMS = ',err_div_RMS
write(6,'(1x,a,es11.4)') 'error divergence Real RMS = ',err_real_div_RMS write(6,'(1x,a,es11.4)') 'error divergence Real RMS = ',err_real_div_RMS
write(6,'(1x,a,es11.4)') 'error divergence post RMS = ',err_post_div_RMS write(6,'(1x,a,es11.4)') 'error divergence post RMS = ',err_post_div_RMS
write(6,'(1x,a,es11.4)') 'error divergence FT max = ',err_div_max write(6,'(1x,a,es11.4)') 'error divergence FT max = ',err_div_max
write(6,'(1x,a,es11.4)') 'error divergence Real max = ',err_real_div_max write(6,'(1x,a,es11.4)') 'error divergence Real max = ',err_real_div_max
flush(6)
endif endif
end function utilities_divergenceRMS end function utilities_divergenceRMS
@ -596,9 +603,13 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
allocate (sTimesC(size_reduced,size_reduced), source =0.0_pReal) allocate (sTimesC(size_reduced,size_reduced), source =0.0_pReal)
temp99_Real = math_Plain3333to99(math_rotate_forward3333(C,rot_BC)) temp99_Real = math_Plain3333to99(math_rotate_forward3333(C,rot_BC))
if(debugGeneral) &
write(6,'(a,/,9(9(2x,f12.7,1x)/))',advance='no') 'Stiffness C rotated / GPa =',& if(debugGeneral) then
write(6,'(/,a)') ' ... updating masked compliance ............................................'
write(6,'(/,a,/,9(9(2x,f12.7,1x)/))',advance='no') ' Stiffness C rotated / GPa =',&
transpose(temp99_Real)/1.e9_pReal transpose(temp99_Real)/1.e9_pReal
flush(6)
endif
k = 0_pInt ! calculate reduced stiffness k = 0_pInt ! calculate reduced stiffness
do n = 1_pInt,9_pInt do n = 1_pInt,9_pInt
if(mask_stressVector(n)) then if(mask_stressVector(n)) then
@ -633,7 +644,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
enddo enddo
if(debugGeneral .or. errmatinv) then ! report if(debugGeneral .or. errmatinv) then ! report
write(formatString, '(I16.16)') size_reduced write(formatString, '(I16.16)') size_reduced
formatString = '(a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))' formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))'
write(6,trim(formatString),advance='no') ' C * S', transpose(matmul(c_reduced,s_reduced)) write(6,trim(formatString),advance='no') ' C * S', transpose(matmul(c_reduced,s_reduced))
write(6,trim(formatString),advance='no') ' S', transpose(s_reduced) write(6,trim(formatString),advance='no') ' S', transpose(s_reduced)
endif endif
@ -645,8 +656,9 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
temp99_real = 0.0_pReal temp99_real = 0.0_pReal
endif endif
if(debugGeneral) & ! report if(debugGeneral) & ! report
write(6,'(a,/,9(9(2x,f12.7,1x)/))',advance='no') 'Masked Compliance * GPa =', & write(6,'(/,a,/,9(9(2x,f12.7,1x)/),/)',advance='no') ' Masked Compliance * GPa =', &
transpose(temp99_Real*1.e9_pReal) transpose(temp99_Real*1.e9_pReal)
flush(6)
utilities_maskedCompliance = math_Plain99to3333(temp99_Real) utilities_maskedCompliance = math_Plain99to3333(temp99_Real)
end function utilities_maskedCompliance end function utilities_maskedCompliance
@ -693,7 +705,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
real(pReal), dimension(6) :: sigma !< cauchy stress in mandel notation real(pReal), dimension(6) :: sigma !< cauchy stress in mandel notation
real(pReal), dimension(6,6) :: dsde !< d sigma / d Epsilon real(pReal), dimension(6,6) :: dsde !< d sigma / d Epsilon
write(6,'(/,a,/)') '... evaluating constitutive response ......................................' write(6,'(/,a)') ' ... evaluating constitutive response ......................................'
if (forwardData) then ! aging results if (forwardData) then ! aging results
calcMode = 1_pInt calcMode = 1_pInt
collectMode = 4_pInt collectMode = 4_pInt
@ -719,7 +731,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
! write(6,'(a,1x,es11.4)') 'max determinant of deformation =', defgradDetMax ! write(6,'(a,1x,es11.4)') 'max determinant of deformation =', defgradDetMax
! write(6,'(a,1x,es11.4)') 'min determinant of deformation =', defgradDetMin ! write(6,'(a,1x,es11.4)') 'min determinant of deformation =', defgradDetMin
! endif ! endif
if (DebugGeneral) write(6,*) 'collect mode: ', collectMode,' calc mode: ', calcMode if (DebugGeneral) write(6,'(/,2(a,i1.1))') ' collect mode: ', collectMode,' calc mode: ', calcMode
flush(6) flush(6)
ielem = 0_pInt ielem = 0_pInt
@ -751,11 +763,12 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P
if (debugRotation) & if (debugRotation) &
write(6,'(a,/,3(3(2x,f12.7,1x)/))',advance='no') 'Piola--Kirchhoff stress (lab) / MPa =',& write(6,'(/,a,/,3(3(2x,f12.7,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',&
math_transpose33(P_av)/1.e6_pReal math_transpose33(P_av)/1.e6_pReal
P_av = math_rotate_forward33(P_av,rotation_BC) P_av = math_rotate_forward33(P_av,rotation_BC)
write(6,'(a,/,3(3(2x,f12.7,1x)/))',advance='no') 'Piola--Kirchhoff stress / MPa =',& write(6,'(/,a,/,3(3(2x,f12.7,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',&
math_transpose33(P_av)/1.e6_pReal math_transpose33(P_av)/1.e6_pReal
end subroutine utilities_constitutiveResponse end subroutine utilities_constitutiveResponse

107
code/IO.f90
View File

@ -1298,6 +1298,7 @@ subroutine IO_error(error_ID,e,i,g,ext_msg)
character(len=*), optional, intent(in) :: ext_msg character(len=*), optional, intent(in) :: ext_msg
character(len=1024) :: msg character(len=1024) :: msg
character(len=1024) :: formatString
select case (error_ID) select case (error_ID)
@ -1414,7 +1415,7 @@ subroutine IO_error(error_ID,e,i,g,ext_msg)
case (831_pInt) case (831_pInt)
msg = 'mask consistency violated in spectral loadcase' msg = 'mask consistency violated in spectral loadcase'
case (832_pInt) case (832_pInt)
msg = 'ill-defined L (each line should be either fully or not at all defined) in spectral loadcase' msg = 'ill-defined L (line party P) in spectral loadcase'
case (834_pInt) case (834_pInt)
msg = 'negative time increment in spectral loadcase' msg = 'negative time increment in spectral loadcase'
case (835_pInt) case (835_pInt)
@ -1438,7 +1439,7 @@ subroutine IO_error(error_ID,e,i,g,ext_msg)
case (846_pInt) case (846_pInt)
msg = 'not a rotation defined for loadcase rotation' msg = 'not a rotation defined for loadcase rotation'
case (847_pInt) case (847_pInt)
msg = 'updating of gamma operator not possible if it is pre calculated' msg = 'update of gamma operator not possible when pre-calculated'
case (850_pInt) case (850_pInt)
msg = 'max number of cut back exceeded' msg = 'max number of cut back exceeded'
case (880_pInt) case (880_pInt)
@ -1453,27 +1454,27 @@ subroutine IO_error(error_ID,e,i,g,ext_msg)
!* Error messages related to parsing of Abaqus input file !* Error messages related to parsing of Abaqus input file
case (900_pInt) case (900_pInt)
msg = 'PARSE ERROR: Improper definition of nodes in input file (Nnodes < 2)' msg = 'improper definition of nodes in input file (Nnodes < 2)'
case (901_pInt) case (901_pInt)
msg = 'PARSE ERROR: No Elements defined in input file (Nelems = 0)' msg = 'no elements defined in input file (Nelems = 0)'
case (902_pInt) case (902_pInt)
msg = 'PARSE ERROR: No Element sets defined in input file (Atleast one *Elset must exist)' msg = 'no element sets defined in input file (No *Elset exists)'
case (903_pInt) case (903_pInt)
msg = 'PARSE ERROR: No Materials defined in input file (Look into section assigments)' msg = 'no materials defined in input file (Look into section assigments)'
case (904_pInt) case (904_pInt)
msg = 'PARSE ERROR: No elements could be assigned for Elset: ' msg = 'no elements could be assigned for Elset: '
case (905_pInt) case (905_pInt)
msg = 'PARSE ERROR: Error in mesh_abaqus_map_materials' msg = 'error in mesh_abaqus_map_materials'
case (906_pInt) case (906_pInt)
msg = 'PARSE ERROR: Error in mesh_abaqus_count_cpElements' msg = 'error in mesh_abaqus_count_cpElements'
case (907_pInt) case (907_pInt)
msg = 'PARSE ERROR: Incorrect size of mesh_mapFEtoCPelem in mesh_abaqus_map_elements; Size cannot be zero' msg = 'size of mesh_mapFEtoCPelem in mesh_abaqus_map_elements'
case (908_pInt) case (908_pInt)
msg = 'PARSE ERROR: Incorrect size of mesh_mapFEtoCPnode in mesh_abaqus_map_nodes; Size cannot be zero' msg = 'size of mesh_mapFEtoCPnode in mesh_abaqus_map_nodes'
case (909_pInt) case (909_pInt)
msg = 'PARSE ERROR: Incorrect size of mesh_node in mesh_abaqus_build_nodes; must be equal to mesh_Nnodes' msg = 'size of mesh_node in mesh_abaqus_build_nodes not equal to mesh_Nnodes'
case (910_pInt) case (910_pInt)
msg = 'PARSE ERROR: Incorrect element type mapping in ' msg = 'incorrect element type mapping in '
!* general error messages !* general error messages
@ -1481,26 +1482,35 @@ subroutine IO_error(error_ID,e,i,g,ext_msg)
case (666_pInt) case (666_pInt)
msg = 'memory leak detected' msg = 'memory leak detected'
case default case default
msg = 'Unknown error number...' msg = 'unknown error number...'
end select end select
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) write(6,'(/,a)') ' +--------------------------------------------------------+'
write(6,'(a38)') '+------------------------------------+' write(6,'(a)') ' + error +'
write(6,'(a38)') '+ error +' write(6,'(a,i3,a)') ' + ',error_ID,' +'
write(6,'(a17,i3,a18)') '+ ',error_ID,' +' write(6,'(a)') ' + +'
write(6,'(a38)') '+ +' write(formatString,'(a,i6.6,a,i6.6,a)') '(1x,a2,a',len(trim(msg)),',',&
write(6,'(a2,a)') '+ ', trim(msg) max(1,60-len(trim(msg))-5),'x,a)'
if (present(ext_msg)) write(6,'(a2,a)') '+ ', trim(ext_msg) write(6,formatString) '+ ', trim(msg),'+'
if (present(ext_msg)) then
write(formatString,'(a,i6.6,a,i6.6,a)') '(1x,a2,a',len(trim(ext_msg)),',',&
max(1,60-len(trim(ext_msg))-5),'x,a)'
write(6,formatString) '+ ', trim(ext_msg),'+'
endif
if (present(e)) then if (present(e)) then
if (present(i) .and. present(g)) then if (present(i)) then
write(6,'(a13,i6,a4,i2,a7,i4,a2)') '+ at element ',e,' IP ',i,' grain ',g,' +' if (present(g)) then
write(6,'(a13,1x,i9,1x,a2,1x,i2,1x,a5,1x,i4,18x,a1)') ' + at element',e,'IP',i,'grain',g,'+'
else else
write(6,'(a18,i6,a14)') '+ at ',e,' +' write(6,'(a13,1x,i9,1x,a2,1x,i2,29x,a1)') ' + at element',e,'IP',i,'+'
endif
else
write(6,'(a13,1x,i9,35x,a1)') ' + at element',e,'+'
endif endif
endif endif
write(6,'(a38)') '+------------------------------------+' write(6,'(a)') ' +--------------------------------------------------------+'
flush(6) flush(6)
call quit(9000_pInt+error_ID) call quit(9000_pInt+error_ID)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
@ -1521,6 +1531,7 @@ subroutine IO_warning(warning_ID,e,i,g,ext_msg)
character(len=*), optional, intent(in) :: ext_msg character(len=*), optional, intent(in) :: ext_msg
character(len=1024) :: msg character(len=1024) :: msg
character(len=1024) :: formatString
select case (warning_ID) select case (warning_ID)
case (34_pInt) case (34_pInt)
@ -1528,47 +1539,52 @@ subroutine IO_warning(warning_ID,e,i,g,ext_msg)
case (35_pInt) case (35_pInt)
msg = 'could not get $DAMASK_NUM_THREADS' msg = 'could not get $DAMASK_NUM_THREADS'
case (40_pInt) case (40_pInt)
msg = 'Found Spectral solver parameter ' msg = 'found spectral solver parameter'
case (41_pInt) case (41_pInt)
msg = 'Found PETSc solver parameter ' msg = 'found PETSc solver parameter'
case (42_pInt) case (42_pInt)
msg = 'parameter has no effect' msg = 'parameter has no effect'
case (47_pInt) case (47_pInt)
msg = 'No valid parameter for FFTW given, using FFTW_PATIENT' msg = 'no valid parameter for FFTW, using FFTW_PATIENT'
case (101_pInt) case (101_pInt)
msg = '+ crystallite debugging off... +' msg = 'crystallite debugging off'
case (600_pInt) case (600_pInt)
msg = '+ crystallite responds elastically +' msg = 'crystallite responds elastically'
case (601_pInt) case (601_pInt)
msg = '+ stiffness close to zero +' msg = 'stiffness close to zero'
case (650_pInt) case (650_pInt)
msg = '+ polar decomposition failed +' msg = 'polar decomposition failed'
case (700_pInt) case (700_pInt)
msg = '+ unknown crystal symmetry +' msg = 'unknown crystal symmetry'
case default case default
msg = '+ unknown warning number... +' msg = 'unknown warning number'
end select end select
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) write(6,'(/,a)') ' +--------------------------------------------------------+'
write(6,'(a38)') '+------------------------------------+' write(6,'(a)') ' + warning +'
write(6,'(a38)') '+ warning +' write(6,'(a,i3,a)') ' + ',warning_ID,' +'
write(6,'(a38)') '+ +' write(6,'(a)') ' + +'
write(6,'(a17,i3,a18)') '+ ',warning_ID,' +' write(formatString,'(a,i6.6,a,i6.6,a)') '(1x,a2,a',len(trim(msg)),',',&
write(6,'(a2,a)') '+ ', trim(msg) max(1,60-len(trim(msg))-5),'x,a)'
if (present(ext_msg)) write(6,'(a2,a)') '+ ', trim(ext_msg) write(6,formatString) '+ ', trim(msg),'+'
if (present(ext_msg)) then
write(formatString,'(a,i6.6,a,i6.6,a)') '(1x,a2,a',len(trim(ext_msg)),',',&
max(1,60-len(trim(ext_msg))-5),'x,a)'
write(6,formatString) '+ ', trim(ext_msg),'+'
endif
if (present(e)) then if (present(e)) then
if (present(i)) then if (present(i)) then
if (present(g)) then if (present(g)) then
write(6,'(a12,1x,i6,1x,a2,1x,i2,1x,a5,1x,i4,a2)') '+ at element',e,'IP',i,'grain',g,' +' write(6,'(a13,1x,i9,1x,a2,1x,i2,1x,a5,1x,i4,18x,a1)') ' + at element',e,'IP',i,'grain',g,'+'
else else
write(6,'(a12,1x,i6,1x,a2,1x,i2,a13)') '+ at element',e,'IP',i,' +' write(6,'(a13,1x,i9,1x,a2,1x,i2,29x,a1)') ' + at element',e,'IP',i,'+'
endif endif
else else
write(6,'(a12,1x,i6,a19)') '+ at element',e,' +' write(6,'(a13,1x,i9,35x,a1)') ' + at element',e,'+'
endif endif
endif endif
write(6,'(a38)') '+------------------------------------+' write(6,'(a)') ' +--------------------------------------------------------+'
flush(6) flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
@ -1597,7 +1613,6 @@ recursive function abaqus_assembleInputFile(unit1,unit2) result(createSuccess)
logical :: createSuccess,fexist logical :: createSuccess,fexist
do do
read(unit2,'(A300)',END=220) line read(unit2,'(A300)',END=220) line
positions = IO_stringPos(line,maxNchunks) positions = IO_stringPos(line,maxNchunks)

9
code/numerics.f90
View File

@ -21,10 +21,12 @@
!############################################################## !##############################################################
module numerics module numerics
!############################################################## !##############################################################
use prec, only: &
use prec, only: pInt, pReal pInt, &
pReal
implicit none implicit none
private
character(len=64), parameter, private :: & character(len=64), parameter, private :: &
numerics_configFile = 'numerics.config' !< name of configuration file numerics_configFile = 'numerics.config' !< name of configuration file
@ -37,7 +39,7 @@ integer(pInt), protected, public :: &
nState = 10_pInt, & !< state loop limit nState = 10_pInt, & !< state loop limit
nStress = 40_pInt, & !< stress loop limit nStress = 40_pInt, & !< stress loop limit
pert_method = 1_pInt !< method used in perturbation technique for tangent pert_method = 1_pInt !< method used in perturbation technique for tangent
integer(pInt) :: numerics_integrationMode = 0_pInt !< integrationMode 1 = central solution ; integrationMode 2 = perturbation, Default 0: undefined, is not read from file integer(pInt), public :: numerics_integrationMode = 0_pInt !< integrationMode 1 = central solution ; integrationMode 2 = perturbation, Default 0: undefined, is not read from file
integer(pInt), dimension(2) , protected, public :: & integer(pInt), dimension(2) , protected, public :: &
numerics_integrator = 1_pInt !< method used for state integration (central & perturbed state), Default 1: fix-point iteration for both states numerics_integrator = 1_pInt !< method used for state integration (central & perturbed state), Default 1: fix-point iteration for both states
real(pReal), protected, public :: & real(pReal), protected, public :: &
@ -151,7 +153,6 @@ subroutine numerics_init
write(6,*) '$Id$' write(6,*) '$Id$'
#include "compilation_info.f90" #include "compilation_info.f90"
!$ call GET_ENVIRONMENT_VARIABLE(NAME='DAMASK_NUM_THREADS',VALUE=DAMASK_NumThreadsString,STATUS=gotDAMASK_NUM_THREADS) ! get environment variable DAMASK_NUM_THREADS... !$ call GET_ENVIRONMENT_VARIABLE(NAME='DAMASK_NUM_THREADS',VALUE=DAMASK_NumThreadsString,STATUS=gotDAMASK_NUM_THREADS) ! get environment variable DAMASK_NUM_THREADS...
!$ if(gotDAMASK_NUM_THREADS /= 0) call IO_warning(35_pInt,ext_msg=DAMASK_NumThreadsString) !$ if(gotDAMASK_NUM_THREADS /= 0) call IO_warning(35_pInt,ext_msg=DAMASK_NumThreadsString)
!$ read(DAMASK_NumThreadsString,'(i6)') DAMASK_NumThreadsInt ! ...convert it to integer... !$ read(DAMASK_NumThreadsString,'(i6)') DAMASK_NumThreadsInt ! ...convert it to integer...