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

numerics got some protected statements and is now reading in the new keywords for the solver selection in small letters. 

spectral solver got cut back facilities + improved output to screen
This commit is contained in:
Martin Diehl 2012-10-02 15:26:56 +00:00
parent 91b7883c2a
commit d5ce49c471
6 changed files with 477 additions and 339 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

170
code/DAMASK_spectral_Driver.f90
View File

@ -49,6 +49,7 @@ program DAMASK_spectral_Driver
restartInc
use numerics, only: &
maxCutBack, &
rotation_tol, &
mySpectralSolver
@ -57,9 +58,10 @@ program DAMASK_spectral_Driver
materialpoint_results
use DAMASK_spectral_Utilities, only: &
boundaryCondition, &
solutionState, &
debugGeneral
tBoundaryCondition, &
tSolutionState, &
debugGeneral, &
cutBack
use DAMASK_spectral_SolverBasic
#ifdef PETSc
@ -68,10 +70,9 @@ program DAMASK_spectral_Driver
#endif
implicit none
type loadCase
type tLoadCase
real(pReal), dimension (3,3) :: rotation = math_I3 ! rotation of BC
type(boundaryCondition) :: P, & ! stress BC
type(tBoundaryCondition) :: P, & ! stress BC
deformation ! deformation BC (Fdot or L)
real(pReal) :: time = 0.0_pReal, & ! length of increment
temperature = 300.0_pReal ! isothermal starting conditions
@ -80,7 +81,7 @@ program DAMASK_spectral_Driver
restartfrequency = 0_pInt, & ! frequency of restart writes
logscale = 0_pInt ! linear/logaritmic time inc flag
logical :: followFormerTrajectory = .true. ! follow trajectory of former loadcase
end type loadCase
end type tLoadCase
!--------------------------------------------------------------------------------------------------
! variables related to information from load case and geom file
@ -88,17 +89,15 @@ program DAMASK_spectral_Driver
logical, dimension(9) :: temp_maskVector !> temporarily from loadcase file when reading in tensors
integer(pInt), parameter :: maxNchunksLoadcase = (1_pInt + 9_pInt)*3_pInt +& ! deformation, rotation, and stress
(1_pInt + 1_pInt)*5_pInt +& ! time, (log)incs, temp, restartfrequency, and outputfrequency
1_pInt, & ! dropguessing
maxNchunksGeom = 7_pInt, & ! 4 identifiers, 3 values
myUnit = 234_pInt
1_pInt ! dropguessing
integer(pInt), dimension(1_pInt + maxNchunksLoadcase*2_pInt) :: positions ! this is longer than needed for geometry parsing
integer(pInt) :: &
N_l = 0_pInt, &
N_t = 0_pInt, &
N_n = 0_pInt, &
N_Fdot = 0_pInt ! number of Fourier points
N_Fdot = 0_pInt, & ! number of Fourier points
myUnit = 234_pInt
character(len=1024) :: &
line
@ -108,15 +107,14 @@ program DAMASK_spectral_Driver
real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal
real(pReal) :: time = 0.0_pReal, time0 = 0.0_pReal, timeinc = 1.0_pReal, timeinc_old = 0.0_pReal ! elapsed time, begin of interval, time interval, previous time interval
real(pReal) :: guessmode
real(pReal), dimension(3,3) :: temp33_Real
integer(pInt) :: i, j, k, l, errorID
integer(pInt) :: i, j, k, l, errorID, cutBackLevel = 0_pInt, stepFraction = 0_pInt
integer(pInt) :: currentLoadcase = 0_pInt, inc, &
totalIncsCounter = 0_pInt,&
notConvergedCounter = 0_pInt, convergedCounter = 0_pInt
character(len=6) :: loadcase_string
type(loadCase), allocatable, dimension(:) :: loadCases
type(solutionState) solres
character(len=1024) :: incInfo
type(tLoadCase), allocatable, dimension(:) :: loadCases
type(tSolutionState) solres
!--------------------------------------------------------------------------------------------------
! init DAMASK (all modules)
@ -126,7 +124,6 @@ program DAMASK_spectral_Driver
write(6,'(a)') ' <<<+- DAMASK_spectral_Driver init -+>>>'
write(6,'(a)') ' $Id$'
#include "compilation_info.f90"
write(6,'(a)') ''
!--------------------------------------------------------------------------------------------------
! reading basic information from load case file and allocate data structure containing load cases
call IO_open_file(myUnit,trim(loadCaseFile))
@ -234,8 +231,7 @@ program DAMASK_spectral_Driver
errorID = 0_pInt
checkLoadcases: do currentLoadCase = 1_pInt, size(loadCases)
write (loadcase_string, '(i6)' ) currentLoadCase
write(6,'(2x,a,i6)') 'load case: ', currentLoadCase
write(6,'(1x,a,i6)') 'load case: ', currentLoadCase
if (.not. loadCases(currentLoadCase)%followFormerTrajectory) write(6,'(2x,a)') 'drop guessing along trajectory'
if (loadCases(currentLoadCase)%deformation%myType=='l') then
@ -247,19 +243,19 @@ program DAMASK_spectral_Driver
else
write(6,'(2x,a)') 'deformation gradient rate:'
endif
write (6,'(3(3(3x,f12.7,1x)/))',advance='no') merge(math_transpose33(loadCases(currentLoadCase)%deformation%values),&
reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(loadCases(currentLoadCase)%deformation%maskLogical))
write (6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'stress / GPa:',&
write(6,'(3(3(3x,f12.7,1x)/))',advance='no') merge(math_transpose33(loadCases(currentLoadCase)%deformation%values),&
reshape(spread(huge(1.0_pReal),1,9),[ 3,3]),transpose(loadCases(currentLoadCase)%deformation%maskLogical))
write(6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'stress / GPa:',&
1e-9_pReal*merge(math_transpose33(loadCases(currentLoadCase)%P%values),&
reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(loadCases(currentLoadCase)%P%maskLogical))
reshape(spread(huge(1.0_pReal),1,9),[ 3,3]),transpose(loadCases(currentLoadCase)%P%maskLogical))
if (any(loadCases(currentLoadCase)%rotation /= math_I3)) &
write (6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',&
write(6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',&
math_transpose33(loadCases(currentLoadCase)%rotation)
write(6,'(2x,a,f12.6)') 'temperature:', loadCases(currentLoadCase)%temperature
write(6,'(2x,a,f12.6)') 'time: ', loadCases(currentLoadCase)%time
write(6,'(2x,a,i5)') 'increments: ', loadCases(currentLoadCase)%incs
write(6,'(2x,a,i5)') 'output frequency: ', loadCases(currentLoadCase)%outputfrequency
write(6,'(2x,a,i5)') 'restart frequency: ', loadCases(currentLoadCase)%restartfrequency
write(6,'(2x,a,i5,/)') 'restart frequency: ', loadCases(currentLoadCase)%restartfrequency
if (any(loadCases(currentLoadCase)%P%maskLogical .eqv. loadCases(currentLoadCase)%deformation%maskLogical)) errorID = 831_pInt ! exclusive or masking only
if (any(loadCases(currentLoadCase)%P%maskLogical .and. transpose(loadCases(currentLoadCase)%P%maskLogical) .and. &
@ -281,8 +277,8 @@ program DAMASK_spectral_Driver
call basic_init()
#ifdef PETSc
case (DAMASK_spectral_SolverBasicPETSC_label)
call BasicPETSC_init()
case (DAMASK_spectral_SolverBasicPETSc_label)
call basicPETSc_init()
case (DAMASK_spectral_SolverAL_label)
call AL_init()
@ -316,7 +312,6 @@ program DAMASK_spectral_Driver
if (debugGeneral) write(6,'(a)') 'Header of result file written out'
endif
!--------------------------------------------------------------------------------------------------
! loopping over loadcases
loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases)
@ -335,12 +330,12 @@ program DAMASK_spectral_Driver
!--------------------------------------------------------------------------------------------------
! forwarding time
timeinc_old = timeinc
if (loadCases(currentLoadCase)%logscale == 0_pInt) then ! linear scale
timeinc = loadCases(currentLoadCase)%time/loadCases(currentLoadCase)%incs ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used
if (loadCases(currentLoadCase)%logscale == 0_pInt) then ! linear scale
timeinc = loadCases(currentLoadCase)%time/loadCases(currentLoadCase)%incs ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used
else
if (currentLoadCase == 1_pInt) then ! 1st currentLoadCase of logarithmic scale
if (currentLoadCase == 1_pInt) then ! 1st currentLoadCase of logarithmic scale
if (inc == 1_pInt) then ! 1st inc of 1st currentLoadCase of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real( 1_pInt-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd
timeinc = loadCases(1)%time*(2.0_pReal**real( 1_pInt-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd
else ! not-1st inc of 1st currentLoadCase of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1_pInt-loadCases(1)%incs ,pReal))
endif
@ -351,63 +346,92 @@ program DAMASK_spectral_Driver
real(loadCases(currentLoadCase)%incs ,pReal)) )
endif
endif
time = time + timeinc
timeinc = timeinc / 2.0_pReal**real(cutBackLevel,pReal)
if(totalIncsCounter >= restartInc) then ! do calculations (otherwise just forwarding)
stepFraction = 0_pInt
!--------------------------------------------------------------------------------------------------
! report begin of new increment
write(6,'(a)') '##################################################################'
write(6,'(A,I5.5,A,es12.5)') 'Increment ', totalIncsCounter, ' Time ',time
subIncLooping: do while (stepFraction/2_pInt**cutBackLevel <1_pInt)
time = time + timeinc
stepFraction = stepFraction + 1_pInt
write(6,'(1/,a)') '###########################################################################'
write(6,'(a,es12.5'//&
',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//&
',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(2_pInt**cutBackLevel)//&
',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') &
'Time', time, &
's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,&
'-', stepFraction, '/', 2_pInt**cutBackLevel,&
' of load case ', currentLoadCase,'/',size(loadCases)
write(incInfo,'(a,'//IO_intOut(totalIncsCounter)//',a,'//IO_intOut(sum(loadCases(:)%incs))//&
',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(2_pInt**cutBackLevel)//')') &
'Inc. ',totalIncsCounter,'/',sum(loadCases(:)%incs),&
'-',stepFraction, '/', 2_pInt**cutBackLevel
select case(myspectralsolver)
select case(myspectralsolver)
case (DAMASK_spectral_SolverBasic_label)
solres = basic_solution (&
guessmode,timeinc,timeinc_old, &
P_BC = loadCases(currentLoadCase)%P, &
F_BC = loadCases(currentLoadCase)%deformation, &
temperature_bc = loadCases(currentLoadCase)%temperature, &
rotation_BC = loadCases(currentLoadCase)%rotation)
case (DAMASK_spectral_SolverBasic_label)
solres = basic_solution (&
incInfo, guessmode,timeinc,timeinc_old, &
P_BC = loadCases(currentLoadCase)%P, &
F_BC = loadCases(currentLoadCase)%deformation, &
temperature_bc = loadCases(currentLoadCase)%temperature, &
rotation_BC = loadCases(currentLoadCase)%rotation)
#ifdef PETSc
case (DAMASK_spectral_SolverBasicPETSC_label)
solres = BasicPETSC_solution (&
guessmode,timeinc,timeinc_old, &
P_BC = loadCases(currentLoadCase)%P, &
F_BC = loadCases(currentLoadCase)%deformation, &
temperature_bc = loadCases(currentLoadCase)%temperature, &
rotation_BC = loadCases(currentLoadCase)%rotation)
case (DAMASK_spectral_SolverBasicPETSC_label)
solres = BasicPETSC_solution (&
incInfo, guessmode,timeinc,timeinc_old, &
P_BC = loadCases(currentLoadCase)%P, &
F_BC = loadCases(currentLoadCase)%deformation, &
temperature_bc = loadCases(currentLoadCase)%temperature, &
rotation_BC = loadCases(currentLoadCase)%rotation)
case (DAMASK_spectral_SolverAL_label)
solres = AL_solution (&
guessmode,timeinc,timeinc_old, &
P_BC = loadCases(currentLoadCase)%P, &
F_BC = loadCases(currentLoadCase)%deformation, &
temperature_bc = loadCases(currentLoadCase)%temperature, &
rotation_BC = loadCases(currentLoadCase)%rotation)
case (DAMASK_spectral_SolverAL_label)
solres = AL_solution (&
incInfo, guessmode,timeinc,timeinc_old, &
P_BC = loadCases(currentLoadCase)%P, &
F_BC = loadCases(currentLoadCase)%deformation, &
temperature_bc = loadCases(currentLoadCase)%temperature, &
rotation_BC = loadCases(currentLoadCase)%rotation)
#endif
end select
write(6,'(a)') ''
write(6,'(a)') '=================================================================='
end select
cutBack = .False.
if(solres%termIll .or. .not. solres%converged) then ! no solution found
if (cutBackLevel < maxCutBack) then ! do cut back
write(6,'(/,a)') 'cut back detected'
cutBack = .True.
stepFraction = (stepFraction - 1_pInt) * 2_pInt**cutBackLevel
cutBackLevel = cutBackLevel + 1_pInt
time = time - timeinc
timeinc_old = timeinc
timeinc = timeinc/2.0_pReal
elseif (solres%termIll) then ! material point model cannot find a solution
call IO_error(850_pInt)
else
guessmode=1.0_pReal ! start guessing after first accepted (not converged) (sub)inc
endif
else
guessmode = 1.0_pReal ! start guessing after first converged (sub)inc
endif
enddo subIncLooping
cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half subincs next inc
if(solres%converged) then
convergedCounter = convergedCounter + 1_pInt
write(6,'(A,'//IO_intOut(totalIncsCounter)//',A)') &
'increment', totalIncsCounter, 'converged'
'increment ', totalIncsCounter, ' converged'
else
write(6,'(A,'//IO_intOut(totalIncsCounter)//',A)') &
'increment', totalIncsCounter, 'NOT converged'
'increment ', totalIncsCounter, ' NOT converged'
notConvergedCounter = notConvergedCounter + 1_pInt
endif
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency
write(6,'(a)') ''
write(6,'(a)') '... writing results to file ......................................'
write(538) materialpoint_results ! write result to file
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency
write(6,'(1/,a)') '... writing results to file ......................................'
write(538) materialpoint_results ! write result to file
endif
endif ! end calculation/forwarding
guessmode = 1.0_pReal ! keep guessing along former trajectory during same currentLoadCase
else !just time forwarding
time = time + timeinc
endif ! end calculation/forwarding
enddo incLooping
enddo loadCaseLooping

177
code/DAMASK_spectral_SolverAL.f90
View File

@ -15,7 +15,7 @@ module DAMASK_spectral_SolverAL
math_I3
use DAMASK_spectral_Utilities, only: &
solutionState
tSolutionState
implicit none
#include <finclude/petscsys.h>
@ -31,47 +31,48 @@ module DAMASK_spectral_SolverAL
#include <finclude/petscsnes.h90>
character (len=*), parameter, public :: &
DAMASK_spectral_SolverAL_label = 'AL'
DAMASK_spectral_SolverAL_label = 'al'
!--------------------------------------------------------------------------------------------------
! derived types
type solutionParams
type tSolutionParams
real(pReal), dimension(3,3) :: P_BC, rotation_BC
real(pReal) :: timeinc
end type solutionParams
end type tSolutionParams
type(solutionParams), private :: params
type(tSolutionParams), private :: params
!--------------------------------------------------------------------------------------------------
! PETSc data
DM, private :: da
SNES, private :: snes
Vec, private :: solution_vec
DM, private :: da
SNES, private :: snes
Vec, private :: solution_vec
!--------------------------------------------------------------------------------------------------
! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, F_lambda_lastInc
real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates
real(pReal), private, dimension(:,:,:), allocatable :: temperature
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, F_lambda_lastInc, F_lambdaDot, Fdot
real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates
real(pReal), private, dimension(:,:,:), allocatable :: temperature
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
real(pReal), private, dimension(3,3) :: &
F_aim = math_I3, &
F_aim_lastInc = math_I3, &
P_av
real(pReal), private, dimension(3,3) :: &
F_aimDot, &
F_aim = math_I3, &
F_aim_lastInc = math_I3, &
P_av
character(len=1024), private :: incInfo
real(pReal), private, dimension(3,3,3,3) :: &
C = 0.0_pReal, C_lastInc = 0.0_pReal, &
S = 0.0_pReal, &
C_scale = 0.0_pReal, &
S_scale = 0.0_pReal
real(pReal), private, dimension(3,3,3,3) :: &
C = 0.0_pReal, &
S = 0.0_pReal, &
C_scale = 0.0_pReal, &
S_scale = 0.0_pReal
real(pReal), private :: err_stress, err_f, err_p
logical, private :: ForwardData
real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal
real(pReal), private :: err_stress, err_f, err_p
logical, private :: ForwardData
real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal
contains
contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields and fills them with data, potentially from restart info
@ -93,7 +94,7 @@ subroutine AL_init()
Utilities_init, &
Utilities_constitutiveResponse, &
Utilities_updateGamma, &
debugrestart
debugRestart
use numerics, only: &
petsc_options
@ -118,14 +119,16 @@ subroutine AL_init()
call Utilities_init()
write(6,'(a)') ''
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 $'
#include "compilation_info.f90"
write(6,'(a)') ''
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,source = F_lastInc) somethin like that should be possible
allocate (F_lambda_lastInc(3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (F_lambdaDot(3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (P (3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal)
allocate (temperature( res(1), res(2),res(3)), source = 0.0_pReal)
@ -214,7 +217,8 @@ subroutine AL_init()
!--------------------------------------------------------------------------------------------------
!> @brief solution for the AL scheme with internal iterations
!--------------------------------------------------------------------------------------------------
type(solutionState) function AL_solution(guessmode,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC)
type(tSolutionState) function &
AL_solution(incInfoIn,guessmode,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC)
use numerics, only: &
update_gamma
@ -227,24 +231,27 @@ subroutine AL_init()
deformed_fft
use IO, only: &
IO_write_JobBinaryFile
use DAMASK_spectral_Utilities, only: &
boundaryCondition, &
tBoundaryCondition, &
Utilities_forwardField, &
Utilities_calculateRate, &
Utilities_maskedCompliance, &
Utilities_updateGamma
Utilities_updateGamma, &
cutBack
use FEsolving, only: &
restartWrite
restartWrite, &
terminallyIll
implicit none
!--------------------------------------------------------------------------------------------------
! input data for solution
real(pReal), intent(in) :: timeinc, timeinc_old, temperature_bc, guessmode
type(boundaryCondition), intent(in) :: P_BC,F_BC
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) :: deltaF_aim, &
F_aim_lab
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
@ -257,6 +264,7 @@ subroutine AL_init()
!--------------------------------------------------------------------------------------------------
! 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))
@ -266,31 +274,56 @@ subroutine AL_init()
write (777,rec=1) C
close(777)
endif
AL_solution%converged =.false.
AL_solution%converged =.false.
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr_psc)
F => xx_psc(0:8,:,:,:)
F_lambda => xx_psc(9:17,:,:,:)
if ( cutBack) then
F_aim = F_aim_lastInc
F_lambda= reshape(F_lambda_lastInc,[9,res(1),res(2),res(3)])
F = reshape(F_lastInc,[9,res(1),res(2),res(3)])
C = C_lastInc
else
!--------------------------------------------------------------------------------------------------
! winding forward of deformation aim in loadcase system
if (F_BC%myType=='l') then ! calculate deltaF_aim from given L and current F
deltaF_aim = timeinc * F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim)
elseif(F_BC%myType=='fdot') then ! deltaF_aim = fDot *timeinc where applicable
deltaF_aim = timeinc * F_BC%maskFloat * F_BC%values
C_lastInc = C
!--------------------------------------------------------------------------------------------------
! calculate rate for aim
if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F
f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim)
elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed
f_aimDot = F_BC%maskFloat * F_BC%values
endif
temp33_Real = F_aim
F_aim = F_aim &
+ guessmode * P_BC%maskFloat * (F_aim - F_aim_lastInc)*timeinc/timeinc_old &
+ deltaF_aim
F_aim_lastInc = temp33_Real
F_aim_lab = math_rotate_backward33(F_aim,rotation_BC) ! boundary conditions from load frame into lab (Fourier) frame
f_aimDot = f_aimDot &
+ guessmode * P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old
F_aim_lastInc = F_aim
!--------------------------------------------------------------------------------------------------
! update coordinates and rate and forward last inc
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc,timeinc_old,guessmode,F_lastInc,reshape(F,[3,3,res(1),res(2),res(3)]))
F_lambdaDot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc,timeinc_old,guessmode,F_lambda_lastInc,reshape(F_lambda,[3,3,res(1),res(2),res(3)]))
F_lastInc = reshape(F,[3,3,res(1),res(2),res(3)])
F_lambda_lastInc = reshape(F_lambda,[3,3,res(1),res(2),res(3)])
call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), &
1.0_pReal,F_lastInc,coordinates)
endif
F_aim = F_aim + f_aimDot * timeinc
!--------------------------------------------------------------------------------------------------
! update local deformation gradient and coordinates
deltaF_aim = math_rotate_backward33(deltaF_aim,rotation_BC)
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr_psc)
F => xx_psc(0:8,:,:,:)
F_lambda => xx_psc(9:17,:,:,:)
call Utilities_forwardField(deltaF_aim,timeinc,timeinc_old,guessmode,F_lastInc,F)
call Utilities_forwardField(deltaF_aim,timeinc,timeinc_old,guessmode,F_lambda_lastInc,F_lambda)
! deltaF_aim = math_rotate_backward33(deltaF_aim,rotation_BC)
F = reshape(Utilities_forwardField(timeinc,F_aim,F_lastInc,Fdot),[9,res(1),res(2),res(3)])
F_lambda = reshape(Utilities_forwardField(timeinc,F_aim,F_lambda_lastInc,F_lambdadot),[9,res(1),res(2),res(3)])
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr_psc)
call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates)
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
@ -305,6 +338,11 @@ subroutine AL_init()
call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr_psc)
call SNESGetConvergedReason(snes,reason,ierr_psc)
AL_solution%termIll = terminallyIll
terminallyIll = .false.
if (reason > 0 ) AL_solution%converged = .true.
end function AL_solution
@ -326,15 +364,18 @@ subroutine AL_init()
wgt
use DAMASK_spectral_Utilities, only: &
field_real, &
Utilities_forwardFFT, &
Utilities_FFTforward, &
Utilities_fourierConvolution, &
Utilities_backwardFFT, &
Utilities_FFTbackward, &
Utilities_constitutiveResponse
use IO, only: IO_intOut
implicit none
integer(pInt) :: i,j,k
integer(pInt), save :: callNo = 3_pInt, reportIter = 0_pInt
real(pReal), dimension(3,3) :: temp33_Real
logical :: report
DMDALocalInfo :: in(DMDA_LOCAL_INFO_SIZE)
PetscScalar, target :: x_scal(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE)
@ -355,14 +396,21 @@ subroutine AL_init()
!--------------------------------------------------------------------------------------------------
! report begin of new iteration
write(6,'(a)') ''
write(6,'(a)') '=================================================================='
write(6,'(4(a,i6.6))') ' @ Iter. ',itmin,' < ',iter,' < ',itmax, ' | # Func. calls = ',nfuncs
if (iter == 0 .and. callNo>2) then
callNo = 0_pInt
reportIter = 0_pInt
endif
if (callNo == 0 .or. mod(callNo,2) == 1_pInt) then
write(6,'(/,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), &
' @ Iter. ', itmin, '<',reportIter, '≤', itmax
write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =',&
math_transpose33(F_aim)
reportIter = reportIter + 1_pInt
endif
callNo = callNo +1_pInt
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,params%timeinc, &
residual_F,C,P_av,ForwardData,params%rotation_BC)
ForwardData = .False.
@ -384,9 +432,9 @@ subroutine AL_init()
!--------------------------------------------------------------------------------------------------
! doing Fourier transform
call Utilities_forwardFFT()
call Utilities_FFTforward()
call Utilities_fourierConvolution(math_rotate_backward33(F_aim,params%rotation_BC))
call Utilities_backwardFFT()
call Utilities_FFTbackward()
!--------------------------------------------------------------------------------------------------
! constructing residual
@ -440,8 +488,7 @@ subroutine AL_init()
write(6,'(a,es14.7)') 'error stress BC = ', err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs)
write(6,'(a,es14.7)') 'error F = ', err_f/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_f_tol
write(6,'(a,es14.7)') 'error P = ', err_p/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_p_tol
return
write(6,'(/,a)') '=========================================================================='
end subroutine AL_converged
!--------------------------------------------------------------------------------------------------

109
code/DAMASK_spectral_SolverBasic.f90
View File

@ -17,7 +17,7 @@ module DAMASK_spectral_SolverBasic
math_I3
use DAMASK_spectral_Utilities, only: &
solutionState
tSolutionState
implicit none
character (len=*), parameter, public :: &
@ -25,7 +25,7 @@ module DAMASK_spectral_SolverBasic
!--------------------------------------------------------------------------------------------------
! pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F, F_lastInc, P
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F, F_lastInc, Fdot, P
real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates
real(pReal), private, dimension(:,:,:), allocatable :: temperature
@ -35,7 +35,7 @@ module DAMASK_spectral_SolverBasic
F_aim = math_I3, &
F_aim_lastInc = math_I3
real(pReal), private,dimension(3,3,3,3) :: &
C = 0.0_pReal
C = 0.0_pReal, C_lastInc = 0.0_pReal
contains
@ -60,28 +60,26 @@ subroutine basic_init()
Utilities_init, &
Utilities_constitutiveResponse, &
Utilities_updateGamma, &
debugrestart
debugRestart
use mesh, only: &
res, &
geomdim
implicit none
integer(pInt) :: i,j,k
real(pReal), dimension(3,3) :: temp33_Real
call Utilities_Init()
write(6,'(a)') ''
write(6,'(a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>'
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>'
write(6,'(a)') ' $Id$'
#include "compilation_info.f90"
write(6,'(a)') ''
allocate (F ( 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 (P ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal)
allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal)
allocate (temperature( res(1), res(2),res(3)), source = 0.0_pReal)
@ -139,7 +137,8 @@ end subroutine basic_init
!--------------------------------------------------------------------------------------------------
!> @brief solution for the basic scheme with internal iterations
!--------------------------------------------------------------------------------------------------
type(solutionState) function basic_solution(guessmode,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC)
type(tSolutionState) function &
basic_solution(incInfo,guessmode,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC)
use numerics, only: &
itmax, &
@ -153,39 +152,40 @@ type(solutionState) function basic_solution(guessmode,timeinc,timeinc_old,P_BC,F
use mesh, only: &
res,&
geomdim, &
deformed_fft
deformed_fft, &
wgt
use IO, only: &
IO_write_JobBinaryFile, &
IO_intOut
use DAMASK_spectral_Utilities, only: &
boundaryCondition, &
tBoundaryCondition, &
field_real, &
Utilities_forwardField, &
Utilities_maskedCompliance, &
Utilities_forwardFFT, &
Utilities_FFTforward, &
Utilities_divergenceRMS, &
Utilities_fourierConvolution, &
Utilities_backwardFFT, &
Utilities_FFTbackward, &
Utilities_updateGamma, &
Utilities_constitutiveResponse
Utilities_constitutiveResponse, &
Utilities_calculateRate
use FEsolving, only: &
restartWrite, &
terminallyIll
use DAMASK_spectral_Utilities, only: cutBack
implicit none
!--------------------------------------------------------------------------------------------------
! input data for solution
real(pReal), intent(in) :: timeinc, timeinc_old, temperature_bc, guessmode
type(boundaryCondition), intent(in) :: P_BC,F_BC
type(tBoundaryCondition), intent(in) :: P_BC,F_BC
character(len=*), intent(in) :: incInfo
real(pReal), dimension(3,3), intent(in) :: rotation_BC
real(pReal), dimension(3,3,3,3) :: S
real(pReal), dimension(3,3) :: deltaF_aim, &
F_aim_lab, &
real(pReal), dimension(3,3), save :: f_aimDot = 0.0_pReal
real(pReal), dimension(3,3) :: F_aim_lab, &
F_aim_lab_lastIter, &
P_av
!--------------------------------------------------------------------------------------------------
@ -208,69 +208,80 @@ type(solutionState) function basic_solution(guessmode,timeinc,timeinc_old,P_BC,F
close(777)
endif
if ( cutBack) then
F_aim = F_aim_lastInc
F = F_lastInc
C = C_lastInc
else
C_lastInc = C
!--------------------------------------------------------------------------------------------------
! winding forward of deformation aim in loadcase system
if (F_BC%myType=='l') then ! calculate deltaF_aim from given L and current F
deltaF_aim = timeinc * F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim)
elseif(F_BC%myType=='fdot') then ! deltaF_aim = fDot *timeinc where applicable
deltaF_aim = timeinc * F_BC%maskFloat * F_BC%values
! calculate rate for aim
if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F
f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim)
elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed
f_aimDot = F_BC%maskFloat * F_BC%values
endif
f_aimDot = f_aimDot &
+ guessmode * P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old
F_aim_lastInc = F_aim
!--------------------------------------------------------------------------------------------------
! update coordinates and rate and forward last inc
call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), &
1.0_pReal,F_lastInc,coordinates)
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc,timeinc_old,guessmode,F_lastInc,F)
F_lastInc = F
endif
temp33_Real = F_aim
F_aim = F_aim &
+ guessmode * P_BC%maskFloat * (F_aim - F_aim_lastInc)*timeinc/timeinc_old &
+ deltaF_aim
F_aim_lastInc = temp33_Real
F_aim_lab = math_rotate_backward33(F_aim,rotation_BC) ! boundary conditions from load frame into lab (Fourier) frame
!--------------------------------------------------------------------------------------------------
! update local deformation gradient and coordinates
deltaF_aim = math_rotate_backward33(deltaF_aim,rotation_BC)
call Utilities_forwardField(deltaF_aim,timeinc,timeinc_old,guessmode,F_lastInc,F)
call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates)
F_aim = F_aim + f_aimDot * timeinc
F = Utilities_forwardField(timeinc,F_aim,F_lastInc,Fdot) !I thin F aim should be rotated here
print*, 'F', sum(sum(sum(F,dim=5),dim=4),dim=3)*wgt
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C)
if (update_gamma) call Utilities_updateGamma(C)
iter = 0_pInt
ForwardData = .True.
convergenceLoop: do while(iter < itmax)
iter = iter + 1_pInt
!--------------------------------------------------------------------------------------------------
! report begin of new iteration
write(6,'(a)') ''
write(6,'(a)') '=================================================================='
write(6,'(3(a,'//IO_intOut(itmax)//'))') ' Iter.', itmin, '<',iter, '<', itmax + 1_pInt
write(6,'(/,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), &
' @ Iter. ', itmin, '<',iter, '≤', itmax
write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =', &
math_transpose33(F_aim)
F_aim_lab_lastIter = math_rotate_backward33(F_aim,rotation_BC)
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
basic_solution%termIll = .false.
call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,timeinc,&
P,C,P_av,ForwardData,rotation_BC)
basic_solution%termIll = terminallyIll
terminallyIll = .False.
ForwardData = .False.
!--------------------------------------------------------------------------------------------------
! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, ((P_av - P_BC%values))) !S = 0.0 for no bc
err_stress = maxval(abs(P_BC%maskFloat * (P_av - P_BC%values))) ! mask = 0.0 for no bc
F_aim_lab = math_rotate_backward33(F_aim,rotation_BC) ! boundary conditions from load frame into lab (Fourier) frame
F_aim = F_aim - math_mul3333xx33(S, ((P_av - P_BC%values))) ! S = 0.0 for no bc
err_stress = maxval(abs(P_BC%maskFloat * (P_av - P_BC%values))) ! mask = 0.0 for no bc
F_aim_lab = math_rotate_backward33(F_aim,rotation_BC) ! boundary conditions from load frame into lab (Fourier) frame
!--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme
field_real = 0.0_pReal
field_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = reshape(P,[res(1),res(2),res(3),3,3],&
order=[4,5,1,2,3]) ! field real has a different order
call Utilities_forwardFFT()
call Utilities_FFTforward()
err_div = Utilities_divergenceRMS()
call Utilities_fourierConvolution(F_aim_lab_lastIter - F_aim_lab)
call Utilities_backwardFFT()
call Utilities_FFTbackward()
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)
write(6,'(/,a)') '=========================================================================='
if ((basic_solution%converged .and. iter > itmin) .or. basic_solution%termIll) exit
enddo convergenceLoop
@ -320,6 +331,7 @@ logical function basic_Converged(err_div,pAvgDiv,err_stress,pAvgStress)
end function basic_Converged
subroutine basic_destroy()
use DAMASK_spectral_Utilities, only: &
@ -330,6 +342,9 @@ subroutine basic_destroy()
end subroutine basic_destroy
end module DAMASK_spectral_SolverBasic

179
code/DAMASK_spectral_SolverBasicPETSC.f90
View File

@ -15,7 +15,7 @@ module DAMASK_spectral_SolverBasicPETSC
math_I3
use DAMASK_spectral_Utilities, only: &
solutionState
tSolutionState
implicit none
#include <finclude/petscsys.h>
@ -35,12 +35,12 @@ module DAMASK_spectral_SolverBasicPETSC
!--------------------------------------------------------------------------------------------------
! derived types
type solutionParams
type tSolutionParams
real(pReal), dimension(3,3) :: P_BC, rotation_BC
real(pReal) :: timeinc
end type solutionParams
end type tSolutionParams
type(solutionParams), private :: params
type(tSolutionParams), private :: params
!--------------------------------------------------------------------------------------------------
! PETSc data
@ -50,7 +50,7 @@ module DAMASK_spectral_SolverBasicPETSC
!--------------------------------------------------------------------------------------------------
! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot
real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates
real(pReal), private, dimension(:,:,:), allocatable :: temperature
@ -60,9 +60,9 @@ module DAMASK_spectral_SolverBasicPETSC
F_aim = math_I3, &
F_aim_lastInc = math_I3, &
P_av
character(len=1024), private :: incInfo
real(pReal), private, dimension(3,3,3,3) :: &
C = 0.0_pReal, &
C = 0.0_pReal, C_lastInc= 0.0_pReal, &
S = 0.0_pReal
real(pReal), private :: err_stress, err_div
@ -78,52 +78,51 @@ subroutine BasicPETSC_init()
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_read_JobBinaryFile, &
IO_write_JobBinaryFile
use IO, only: &
IO_read_JobBinaryFile, &
IO_write_JobBinaryFile
use FEsolving, only: &
restartInc
use FEsolving, only: &
restartInc
use DAMASK_interface, only: &
getSolverJobName
use DAMASK_interface, only: &
getSolverJobName
use DAMASK_spectral_Utilities, only: &
Utilities_init, &
Utilities_constitutiveResponse, &
Utilities_updateGamma, &
debugrestart
use DAMASK_spectral_Utilities, only: &
Utilities_init, &
Utilities_constitutiveResponse, &
Utilities_updateGamma, &
debugRestart
use numerics, only: &
petsc_options
use numerics, only: &
petsc_options
use mesh, only: &
res, &
geomdim, &
mesh_NcpElems
use mesh, only: &
res, &
geomdim, &
mesh_NcpElems
use math, only: &
math_invSym3333
use math, only: &
math_invSym3333
implicit none
implicit none
integer(pInt) :: i,j,k
real(pReal), dimension(:,:,:,:,:), allocatable :: P
integer(pInt) :: i,j,k
real(pReal), dimension(:,:,:,:,:), allocatable :: P
PetscErrorCode :: ierr_psc
PetscObject :: dummy
PetscMPIInt :: rank
PetscScalar, pointer :: xx_psc(:,:,:,:), F(:,:,:,:)
PetscErrorCode :: ierr_psc
PetscObject :: dummy
PetscMPIInt :: rank
PetscScalar, pointer :: xx_psc(:,:,:,:), F(:,:,:,:)
call Utilities_init()
write(6,'(a)') ''
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 $'
#include "compilation_info.f90"
write(6,'(a)') ''
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 (P (3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal)
allocate (temperature( res(1), res(2),res(3)), source = 0.0_pReal)
@ -199,8 +198,8 @@ subroutine BasicPETSC_init()
!--------------------------------------------------------------------------------------------------
!> @brief solution for the Basic PETSC scheme with internal iterations
!--------------------------------------------------------------------------------------------------
type(solutionState) function BasicPETSC_solution(guessmode,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC)
type(tSolutionState) function &
basicPETSc_solution(incInfoIn,guessmode,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC)
use numerics, only: &
update_gamma
use math, only: &
@ -212,23 +211,25 @@ subroutine BasicPETSC_init()
deformed_fft
use IO, only: &
IO_write_JobBinaryFile
use DAMASK_spectral_Utilities, only: &
boundaryCondition, &
tBoundaryCondition, &
Utilities_calculateRate, &
Utilities_forwardField, &
Utilities_maskedCompliance, &
Utilities_updateGamma
Utilities_updateGamma, &
cutBack
use FEsolving, only: &
restartWrite
restartWrite, &
terminallyIll
implicit none
!--------------------------------------------------------------------------------------------------
! input data for solution
real(pReal), intent(in) :: timeinc, timeinc_old, temperature_bc, guessmode
type(boundaryCondition), intent(in) :: P_BC,F_BC
type(tBoundaryCondition), intent(in) :: P_BC,F_BC
real(pReal), dimension(3,3), intent(in) :: rotation_BC
real(pReal), dimension(3,3) :: deltaF_aim, &
character(len=*), intent(in) :: incInfoIn
real(pReal), dimension(3,3),save :: F_aimDot=0.0_pReal
real(pReal), dimension(3,3) :: &
F_aim_lab
!--------------------------------------------------------------------------------------------------
! loop variables, convergence etc.
@ -237,11 +238,12 @@ subroutine BasicPETSC_init()
!--------------------------------------------------------------------------------------------------
!
PetscScalar, pointer :: xx_psc(:,:,:,:), F(:,:,:,:)
PetscErrorCode ierr_psc
SNESConvergedReason reason
PetscErrorCode :: ierr_psc
SNESConvergedReason :: reason
!--------------------------------------------------------------------------------------------------
! 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))
@ -251,29 +253,43 @@ subroutine BasicPETSC_init()
write (777,rec=1) C
close(777)
endif
BasicPETSC_solution%converged =.false.
!--------------------------------------------------------------------------------------------------
! winding forward of deformation aim in loadcase system
if (F_BC%myType=='l') then ! calculate deltaF_aim from given L and current F
deltaF_aim = timeinc * F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim)
elseif(F_BC%myType=='fdot') then ! deltaF_aim = fDot *timeinc where applicable
deltaF_aim = timeinc * F_BC%maskFloat * F_BC%values
endif
temp33_Real = F_aim
F_aim = F_aim &
+ guessmode * P_BC%maskFloat * (F_aim - F_aim_lastInc)*timeinc/timeinc_old &
+ deltaF_aim
F_aim_lastInc = temp33_Real
F_aim_lab = math_rotate_backward33(F_aim,rotation_BC) ! boundary conditions from load frame into lab (Fourier) frame
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr_psc)
F => xx_psc(0:8,:,:,:)
if ( cutBack) then
F_aim = F_aim_lastInc
F = reshape(F_lastInc,[9,res(1),res(2),res(3)])
C = C_lastInc
else
C_lastInc = C
!--------------------------------------------------------------------------------------------------
! update local deformation gradient and coordinates
deltaF_aim = math_rotate_backward33(deltaF_aim,rotation_BC)
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr_psc)
F => xx_psc(0:8,:,:,:)
call Utilities_forwardField(deltaF_aim,timeinc,timeinc_old,guessmode,F_lastInc,F)
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr_psc)
call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates)
! calculate rate for aim
if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F
f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim)
elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed
f_aimDot = F_BC%maskFloat * F_BC%values
endif
f_aimDot = f_aimDot &
+ guessmode * P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old
F_aim_lastInc = F_aim
!--------------------------------------------------------------------------------------------------
! update coordinates and rate and forward last inc
call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), &
1.0_pReal,F_lastInc,coordinates)
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc,timeinc_old,guessmode,F_lastInc,reshape(F,[3,3,res(1),res(2),res(3)]))
F_lastInc = reshape(F,[3,3,res(1),res(2),res(3)])
endif
F_aim = F_aim + f_aimDot * timeinc
F = reshape(Utilities_forwardField(timeinc,F_aim,F_lastInc,Fdot),[9,res(1),res(2),res(3)])
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr_psc)
call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates)
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
@ -288,6 +304,9 @@ subroutine BasicPETSC_init()
call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr_psc)
call SNESGetConvergedReason(snes,reason,ierr_psc)
basicPETSc_solution%termIll = terminallyIll
terminallyIll = .false.
BasicPETSC_solution%converged =.false.
if (reason > 0 ) BasicPETSC_solution%converged = .true.
end function BasicPETSC_solution
@ -309,12 +328,12 @@ subroutine BasicPETSC_init()
wgt
use DAMASK_spectral_Utilities, only: &
field_real, &
Utilities_forwardFFT, &
Utilities_FFTforward, &
Utilities_FFTbackward, &
Utilities_fourierConvolution, &
Utilities_backwardFFT, &
Utilities_constitutiveResponse, &
Utilities_divergenceRMS
use IO, only : IO_intOut
implicit none
real(pReal), dimension(3,3) :: F_aim_lab_lastIter, F_aim_lab
@ -327,7 +346,6 @@ subroutine BasicPETSC_init()
PetscInt :: iter, nfuncs
PetscObject :: dummy
PetscErrorCode :: ierr_psc
F => x_scal(:,:,:,:,:)
residual_F => f_scal(:,:,:,:,:)
@ -336,9 +354,8 @@ subroutine BasicPETSC_init()
!--------------------------------------------------------------------------------------------------
! report begin of new iteration
write(6,'(a)') ''
write(6,'(a)') '=================================================================='
write(6,'(4(a,i6.6))') ' @ Iter. ',itmin,' < ',iter,' < ',itmax, ' | # Func. calls = ',nfuncs
write(6,'(/,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), &
' @ Iter. ', itmin, '<',iter, '≤', itmax
write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =',&
math_transpose33(F_aim)
F_aim_lab_lastIter = math_rotate_backward33(F_aim,params%rotation_BC)
@ -347,7 +364,7 @@ subroutine BasicPETSC_init()
! evaluate constitutive response
call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,params%timeinc, &
residual_F,C,P_av,ForwardData,params%rotation_BC)
ForwardData = .False.
ForwardData = .false.
!--------------------------------------------------------------------------------------------------
! stress BC handling
@ -360,15 +377,15 @@ subroutine BasicPETSC_init()
field_real = 0.0_pReal
field_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = reshape(residual_F,[res(1),res(2),res(3),3,3],&
order=[4,5,1,2,3]) ! field real has a different order
call Utilities_forwardFFT()
call Utilities_FFTforward()
err_div = Utilities_divergenceRMS()
call Utilities_fourierConvolution(F_aim_lab_lastIter - F_aim_lab)
call Utilities_backwardFFT()
call Utilities_FFTbackward()
!--------------------------------------------------------------------------------------------------
! constructing residual
residual_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])
write(6,'(/,a)') '=========================================================================='
end subroutine BasicPETSC_formResidual
!--------------------------------------------------------------------------------------------------

112
code/DAMASK_spectral_Utilities.f90
View File

@ -33,7 +33,7 @@ module DAMASK_spectral_Utilities
real(pReal), private, dimension(3,3,3,3) :: C_ref
real(pReal), public, dimension(:,:,:,:,:), pointer :: field_real
logical, public :: cutBack =.false.
!--------------------------------------------------------------------------------------------------
! debug fftw
type(C_PTR), private :: plan_scalarField_forth, plan_scalarField_back
@ -45,7 +45,7 @@ module DAMASK_spectral_Utilities
type(C_PTR), private :: plan_divergence
real(pReal), private, dimension(:,:,:,:), pointer :: divergence_real
complex(pReal), private, dimension(:,:,:,:), pointer :: divergence_fourier
real(pReal), dimension(:,:,:,:), allocatable :: divergence_post
real(pReal), private, dimension(:,:,:,:), allocatable :: divergence_post
!--------------------------------------------------------------------------------------------------
!variables controlling debugging
@ -53,18 +53,18 @@ module DAMASK_spectral_Utilities
!--------------------------------------------------------------------------------------------------
! derived types
type solutionState
logical :: converged = .false.
type tSolutionState
logical :: converged = .true.
logical :: regrid = .false.
logical :: termIll = .false.
end type solutionState
end type tSolutionState
type boundaryCondition
type tBoundaryCondition
real(pReal), dimension(3,3) :: values = 0.0_pReal
real(pReal), dimension(3,3) :: maskFloat = 0.0_pReal
logical, dimension(3,3) :: maskLogical = .false.
character(len=64) :: myType = 'None'
end type boundaryCondition
end type tBoundaryCondition
contains
@ -105,7 +105,7 @@ subroutine Utilities_init()
write(6,'(a)') ''
write(6,'(a)') ' <<<+- DAMASK_spectralSolver Utilities init -+>>>'
write(6,'(a)') ' <<<+- DAMASK_Utilities init -+>>>'
write(6,'(a)') ' $Id$'
#include "compilation_info.f90"
write(6,'(a)') ''
@ -238,7 +238,7 @@ end subroutine Utilities_updateGamma
!> In case of debugging the FFT, also one component of the tensor (specified by row and column)
!> is independetly transformed complex to complex and compared to the whole tensor transform
!--------------------------------------------------------------------------------------------------
subroutine Utilities_forwardFFT(row,column)
subroutine Utilities_FFTforward(row,column)
use mesh, only : &
virt_dim
use math, only: &
@ -286,7 +286,7 @@ subroutine Utilities_forwardFFT(row,column)
if(res(3)>1_pInt) &
field_fourier (1:res1_red,1:res(2), res(3)/2_pInt+1_pInt,1:3,1:3)&
= cmplx(0.0_pReal,0.0_pReal,pReal)
end subroutine Utilities_forwardFFT
end subroutine Utilities_FFTforward
!--------------------------------------------------------------------------------------------------
@ -296,7 +296,7 @@ end subroutine Utilities_forwardFFT
!> is independetly transformed complex to complex and compared to the whole tensor transform
!> results is weighted by number of points stored in wgt
!--------------------------------------------------------------------------------------------------
subroutine Utilities_backwardFFT(row,column)
subroutine Utilities_FFTbackward(row,column)
implicit none
integer(pInt), intent(in), optional :: row, column
@ -335,7 +335,7 @@ subroutine Utilities_backwardFFT(row,column)
endif
field_real = field_real * wgt
end subroutine Utilities_backwardFFT
end subroutine Utilities_FFTbackward
!--------------------------------------------------------------------------------------------------
@ -353,9 +353,8 @@ subroutine Utilities_fourierConvolution(fieldAim)
integer(pInt) :: i, j, k, l, m, n, o
complex(pReal), dimension(3,3) :: temp33_complex
write(6,'(a)') ''
write(6,'(a)') '... doing convolution .................'
write(6,'(a)') ''
write(6,'(/,a)') '... doing convolution .....................................................'
!--------------------------------------------------------------------------------------------------
! to the actual spectral method calculation (mechanical equilibrium)
@ -398,8 +397,7 @@ real(pReal) function Utilities_divergenceRMS()
err_div_max, err_real_div_max
complex(pReal), dimension(3) :: temp3_complex
write(6,'(a)') ''
write(6,'(a)') '... calculating divergence .................'
write(6,'(/,a)') '... calculating divergence ................................................'
!--------------------------------------------------------------------------------------------------
! calculating RMS divergence criterion in Fourier space
@ -477,7 +475,7 @@ function Utilities_maskedCompliance(rot_BC,mask_stressVector,C)
allocate (s_reduced(size_reduced,size_reduced), source =0.0_pReal)
allocate (sTimesC(size_reduced,size_reduced), source =0.0_pReal)
C_lastInc = math_rotate_forward3333(C*0.95_pReal+0.5_pReal*C_ref,rot_BC) ! calculate stiffness from former inc
C_lastInc = math_rotate_forward3333(C*0.75_pReal+0.05_pReal*C_ref,rot_BC) ! calculate stiffness from former inc
temp99_Real = math_Plain3333to99(C_lastInc)
k = 0_pInt ! build reduced stiffness
do n = 1_pInt,9_pInt
@ -537,76 +535,98 @@ subroutine Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti
implicit none
real(pReal), dimension(res(1),res(2),res(3)) :: temperature
real(pReal), dimension(res(1),res(2),res(3),3) :: coordinates
real(pReal), dimension(res(1),res(2),res(3),3), intent(in) :: coordinates
real(pReal), dimension(3,3,res(1),res(2),res(3)) :: F,F_lastInc, P
real(pReal), dimension(3,3,res(1),res(2),res(3)), intent(in) :: F,F_lastInc
real(pReal), dimension(3,3,res(1),res(2),res(3)) :: P
real(pReal) :: timeinc
logical :: ForwardData
logical, intent(in) :: forwardData
integer(pInt) :: i, j, k, ielem
integer(pInt) :: CPFEM_mode
integer(pInt) :: calcMode, collectMode
real(pReal), dimension(3,3,3,3) :: dPdF, C
real(pReal), dimension(6) :: sigma ! cauchy stress
real(pReal), dimension(6,6) :: dsde
real(pReal), dimension(3,3) :: P_av, rotation_BC
real(pReal), dimension(3,3), intent(in) :: rotation_BC
real(pReal), dimension(3,3) :: P_av
write(6,'(a)') ''
write(6,'(a)') '... evaluating constitutive response .................'
write(6,'(a)') ''
if (ForwardData) then
CPFEM_mode = 1_pInt
write(6,'(/,a,/)') '... evaluating constitutive response ......................................'
if (forwardData) then
calcMode = 1_pInt
collectMode = 4_pInt
else
CPFEM_mode = 2_pInt
calcMode = 2_pInt
collectMode = 3_pInt
endif
if (cutBack) then
calcMode = 2_pInt
collectMode = 5_pInt
endif
print*, 'collect mode', collectMode
print*, 'calc mode', calcMode
ielem = 0_pInt
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem + 1_pInt
call CPFEM_general(3_pInt,& ! collect cycle
call CPFEM_general(collectMode,& ! collect cycle
coordinates(i,j,k,1:3), F_lastInc(1:3,1:3,i,j,k),F(1:3,1:3,i,j,k), &
temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF)
collectMode = 3_pInt
enddo; enddo; enddo
P = 0.0_pReal ! needed because of the padding for FFTW
P = 0.0_pReal ! needed because of the padding for FFTW
C = 0.0_pReal
ielem = 0_pInt
call debug_reset()
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem + 1_pInt
call CPFEM_general(CPFEM_mode,& ! first element in first iteration retains CPFEM_mode 1,
coordinates(i,j,k,1:3),F_lastInc(1:3,1:3,i,j,k), F(1:3,1:3,i,j,k), & ! others get 2 (saves winding forward effort)
call CPFEM_general(calcMode,& ! first element in first iteration retains CPFEM_mode 1,
coordinates(i,j,k,1:3),F_lastInc(1:3,1:3,i,j,k), F(1:3,1:3,i,j,k), & ! others get 2 (saves winding forward effort)
temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF)
CPFEM_mode = 2_pInt
calcMode = 2_pInt
C = C + dPdF
enddo; enddo; enddo
call debug_info()
P_av = math_rotate_forward33(sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt,rotation_BC) !average of P rotated
restartWrite = .false.
cutBack = .false.
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
C = C * wgt
end subroutine Utilities_constitutiveResponse
subroutine Utilities_forwardField(delta_aim,timeinc,timeinc_old,guessmode,field_lastInc,field)
function Utilities_calculateRate(delta_aim,timeinc,timeinc_old,guessmode,field_lastInc,field)
implicit none
real(pReal), intent(in), dimension(3,3) :: delta_aim
real(pReal), intent(in) :: timeinc, timeinc_old, guessmode
real(pReal), intent(inout), dimension(3,3,res(1),res(2),res(3)) :: field_lastInc,field
real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: field_lastInc,field
real(pReal), dimension(3,3,res(1),res(2),res(3)) :: Utilities_calculateRate
if (guessmode == 1.0_pReal) then
field = field + (field-field_lastInc) * timeinc/timeinc_old
field_lastInc = (field + field_lastInc * timeinc/timeinc_old) /(1.0_pReal + timeinc/timeinc_old)
Utilities_calculateRate = (field-field_lastInc) / timeinc_old
else
field_lastInc = field
field = field + spread(spread(spread(delta_aim,3,res(1)),4,res(2)),5,res(3))
Utilities_calculateRate = spread(spread(spread(delta_aim,3,res(1)),4,res(2)),5,res(3))/timeinc
endif
end subroutine Utilities_forwardField
end function Utilities_calculateRate
function Utilities_forwardField(timeinc,aim,field_lastInc,rate)
implicit none
real(pReal), intent(in) :: timeinc
real(pReal), intent(in), dimension(3,3) :: aim
real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: field_lastInc,rate
real(pReal), dimension(3,3,res(1),res(2),res(3)) :: Utilities_forwardField
real(pReal), dimension(3,3) :: fieldDiff
Utilities_forwardField = field_lastInc + rate*timeinc
fieldDiff = sum(sum(sum(Utilities_forwardField,dim=5),dim=4),dim=3)*wgt - aim
Utilities_forwardField = Utilities_forwardField - &
spread(spread(spread(fieldDiff,3,res(1)),4,res(2)),5,res(3))
end function Utilities_forwardField
real(pReal) function Utilities_getFilter(k)
@ -634,6 +654,7 @@ real(pReal) function Utilities_getFilter(k)
end function Utilities_getFilter
subroutine Utilities_destroy()
implicit none
@ -649,4 +670,5 @@ subroutine Utilities_destroy()
end subroutine Utilities_destroy
end module DAMASK_spectral_Utilities

51
code/numerics.f90
View File

@ -25,20 +25,23 @@ module numerics
use prec, only: pInt, pReal
implicit none
character(len=64), parameter, private ::&
character(len=64), parameter, private :: &
numerics_configFile = 'numerics.config' !< name of configuration file
integer(pInt) :: iJacoStiffness = 1_pInt, & !< frequency of stiffness update
integer(pInt), protected, public :: &
iJacoStiffness = 1_pInt, & !< frequency of stiffness update
iJacoLpresiduum = 1_pInt, & !< frequency of Jacobian update of residuum in Lp
nHomog = 20_pInt, & !< homogenization loop limit (only for debugging info, loop limit is determined by "subStepMinHomog")
nMPstate = 10_pInt, & !< materialpoint state loop limit
nCryst = 20_pInt, & !< crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst")
nState = 10_pInt, & !< state loop limit
nStress = 40_pInt, & !< stress loop limit
pert_method = 1_pInt, & !< method used in perturbation technique for tangent
numerics_integrationMode = 0_pInt !< integrationMode 1 = central solution ; integrationMode 2 = perturbation, Default 0: undefined, is not read from file
integer(pInt), dimension(2) :: numerics_integrator = 1_pInt !< method used for state integration (central & perturbed state), Default 1: fix-point iteration for both states
real(pReal) :: relevantStrain = 1.0e-7_pReal, & !< strain increment considered significant (used by crystallite to determine whether strain inc is considered significant)
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), dimension(2) , protected, public :: &
numerics_integrator = 1_pInt !< method used for state integration (central & perturbed state), Default 1: fix-point iteration for both states
real(pReal), protected, public :: &
relevantStrain = 1.0e-7_pReal, & !< strain increment considered significant (used by crystallite to determine whether strain inc is considered significant)
defgradTolerance = 1.0e-7_pReal, & !< deviation of deformation gradient that is still allowed (used by CPFEM to determine outdated ffn1)
pert_Fg = 1.0e-7_pReal, & !< strain perturbation for FEM Jacobi
subStepMinCryst = 1.0e-3_pReal, & !< minimum (relative) size of sub-step allowed during cutback in crystallite
@ -65,38 +68,48 @@ real(pReal) :: relevantStrain = 1.0e-7_pReal, &
maxVolDiscr_RGC = 1.0e-5_pReal, & !< threshold of maximum volume discrepancy allowed
volDiscrMod_RGC = 1.0e+12_pReal, & !< stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint)
volDiscrPow_RGC = 5.0_pReal !< powerlaw penalty for volume discrepancy
logical :: analyticJaco = .false. !< use analytic Jacobian or perturbation, Default .false.: calculate Jacobian using perturbations
logical, protected, public :: &
analyticJaco = .false. !< use analytic Jacobian or perturbation, Default .false.: calculate Jacobian using perturbations
!* Random seeding parameters
integer(pInt) :: fixedSeed = 0_pInt !< fixed seeding for pseudo-random number generator, Default 0: use random seed
integer(pInt), protected, public :: &
fixedSeed = 0_pInt !< fixed seeding for pseudo-random number generator, Default 0: use random seed
!* OpenMP variable
integer(pInt) :: DAMASK_NumThreadsInt = 0_pInt !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive
integer(pInt), protected, public :: &
DAMASK_NumThreadsInt = 0_pInt !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive
!* spectral parameters:
#ifdef Spectral
real(pReal) :: err_div_tol = 0.1_pReal, & !< Div(P)/avg(P)*meter
real(pReal), protected, public :: &
err_div_tol = 0.1_pReal, & !< Div(P)/avg(P)*meter
err_stress_tolrel = 0.01_pReal, & !< relative tolerance for fullfillment of stress BC, Default: 0.01 allowing deviation of 1% of maximum stress
err_stress_tolabs = huge(1.0_pReal), & !< absolute tolerance for fullfillment of stress BC, Default: 0.01 allowing deviation of 1% of maximum stress
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
character(len=64) :: 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), 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 :: &
myspectralsolver = 'basic' , & !< spectral solution method
myfilter = 'none' !< spectral filtering method
character(len=1024) :: petsc_options = '-snes_type ngmres -snes_ngmres_anderson -snes_view'
integer(pInt) :: fftw_planner_flag = 32_pInt, & !< conversion of fftw_plan_mode to integer, basically what is usually done in the include file of fftw
character(len=1024), protected, public :: &
petsc_options = '-snes_type ngmres -snes_ngmres_anderson -snes_view'
integer(pInt), protected, public :: &
fftw_planner_flag = 32_pInt, & !< conversion of fftw_plan_mode to integer, basically what is usually done in the include file of fftw
itmax = 20_pInt, & !< maximum number of iterations
itmin = 2_pInt, & !< minimum number of iterations
maxCutBack = 3_pInt !< max number of cut backs
logical :: memory_efficient = .true., & !< for fast execution (pre calculation of gamma_hat), Default .true.: do not precalculate
logical, protected , public :: &
memory_efficient = .true., & !< for fast execution (pre calculation of gamma_hat), Default .true.: do not precalculate
divergence_correction = .false., & !< correct divergence calculation in fourier space, Default .false.: no correction
update_gamma = .false. !< update gamma operator with current stiffness, Default .false.: use initial stiffness
#endif
public :: numerics_init
CONTAINS
contains
!*******************************************
! initialization subroutine
@ -256,9 +269,9 @@ subroutine numerics_init
case ('fftw_timelimit')
fftw_timelimit = IO_floatValue(line,positions,2_pInt)
case ('fftw_plan_mode')
fftw_plan_mode = IO_stringValue(line,positions,2_pInt)
fftw_plan_mode = IO_lc(IO_stringValue(line,positions,2_pInt))
case ('myfilter')
myfilter = IO_stringValue(line,positions,2_pInt)
myfilter = IO_lc(IO_stringValue(line,positions,2_pInt))
case ('rotation_tol')
rotation_tol = IO_floatValue(line,positions,2_pInt)
case ('divergence_correction')
@ -269,7 +282,7 @@ subroutine numerics_init
case ('petsc_options')
petsc_options = trim(line(positions(4):))
case ('myspectralsolver')
myspectralsolver = IO_stringValue(line,positions,2_pInt)
myspectralsolver = IO_lc(IO_stringValue(line,positions,2_pInt))
case ('err_f_tol')
err_f_tol = IO_floatValue(line,positions,2_pInt)
case ('err_p_tol')