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

intermediate, not working state of new solver structure

This commit is contained in:
Martin Diehl 2012-07-20 15:33:13 +00:00
parent 9887f42383
commit bd9e81fbec
3 changed files with 430 additions and 548 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

330
code/DAMASK_spectralDriver.f90
View File

@ -34,37 +34,113 @@
! !
! MPI fuer Eisenforschung, Duesseldorf ! MPI fuer Eisenforschung, Duesseldorf
#include "spectral_quit.f90"
program DAMASK_spectral program DAMASK_spectralDriver
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
use DAMASK_interface, only: &
DAMASK_interface_init, &
loadCaseFile, &
geometryFile, &
getSolverWorkingDirectoryName, &
getSolverJobName, &
appendToOutFile
use prec, only: & use prec, only: &
pInt, & pInt, &
pReal pReal, &
DAMASK_NaN
use IO, only: & use IO, only: &
IO_error,& IO_isBlank, &
IO_open_file, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_error, &
IO_lc, &
IO_read_jobBinaryFile, &
IO_write_jobBinaryFile IO_write_jobBinaryFile
use math use math
use mesh, only : &
mesh_spectral_getResolution, &
mesh_spectral_getDimension, &
mesh_spectral_getHomogenization
use CPFEM, only: &
CPFEM_initAll
use FEsolving, only: & use FEsolving, only: &
restartWrite, & restartWrite, &
restartInc restartInc
use numerics, only: &
rotation_tol
use homogenization, only: & use homogenization, only: &
materialpoint_sizeResults, & materialpoint_sizeResults, &
materialpoint_results materialpoint_results
use DAMASK_spectralSovler use DAMASK_spectralSolver !, only: &
!solution, &
!solution_t
implicit none implicit none
!--------------------------------------------------------------------------------------------------
! variables related to information from load case and geom file
real(pReal), dimension(9) :: &
temp_valueVector !> temporarily from loadcase file when reading in tensors
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
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, &
Npoints ! number of Fourier points
character(len=1024) :: &
line
type(bc_type), allocatable, dimension(:) :: bc
type(solution_t) solres
type(init) initres
!--------------------------------------------------------------------------------------------------
! BC related information
real(pReal), dimension(3,3) :: &
F_aim = math_I3, &
F_aim_lastInc = math_I3, &
mask_stress, &
mask_defgrad, &
deltaF_aim
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! loop variables, convergence etc. ! loop variables, convergence etc.
integer(pInt) :: i, j, k, l, m, n, p, errorID 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
real(pReal) :: guessmode
real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal
real(pReal), dimension(3,3) :: temp33_Real
integer(pInt) :: i, j, k, p, errorID
integer(pInt) :: N_Loadcases, loadcase = 0_pInt, inc, &
totalIncsCounter = 0_pInt,&
notConvergedCounter = 0_pInt, convergedCounter = 0_pInt
character(len=6) :: loadcase_string
call DAMASK_interface_init
call DAMASK_interface_init
write(6,'(a)') '' write(6,'(a)') ''
write(6,'(a)') ' <<<+- DAMASK_spectral init -+>>>' write(6,'(a)') ' <<<+- DAMASK_spectral init -+>>>'
write(6,'(a)') ' $Id$' write(6,'(a)') ' $Id$'
@ -72,6 +148,7 @@ program DAMASK_spectral
write(6,'(a)') ' Working Directory: ',trim(getSolverWorkingDirectoryName()) write(6,'(a)') ' Working Directory: ',trim(getSolverWorkingDirectoryName())
write(6,'(a)') ' Solver Job Name: ',trim(getSolverJobName()) write(6,'(a)') ' Solver Job Name: ',trim(getSolverJobName())
write(6,'(a)') '' write(6,'(a)') ''
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reading the load case file and allocate data structure containing load cases ! reading the load case file and allocate data structure containing load cases
call IO_open_file(myUnit,trim(loadCaseFile)) call IO_open_file(myUnit,trim(loadCaseFile))
@ -108,63 +185,63 @@ program DAMASK_spectral
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
loadcase = loadcase + 1_pInt loadcase = loadcase + 1_pInt
positions = IO_stringPos(line,maxNchunksLoadcase) positions = IO_stringPos(line,maxNchunksLoadcase)
do j = 1_pInt,maxNchunksLoadcase do i = 1_pInt,maxNchunksLoadcase
select case (IO_lc(IO_stringValue(line,positions,j))) select case (IO_lc(IO_stringValue(line,positions,i)))
case('fdot','dotf','l','velocitygrad','velgrad','velocitygradient') ! assign values for the deformation BC matrix case('fdot','dotf','l','velocitygrad','velgrad','velocitygradient') ! assign values for the deformation BC matrix
bc(loadcase)%velGradApplied = & bc(loadcase)%velGradApplied = &
(IO_lc(IO_stringValue(line,positions,j)) == 'l'.or. & ! in case of given L, set flag to true (IO_lc(IO_stringValue(line,positions,i)) == 'l'.or. & ! in case of given L, set flag to true
IO_lc(IO_stringValue(line,positions,j)) == 'velocitygrad'.or.& IO_lc(IO_stringValue(line,positions,i)) == 'velocitygrad'.or.&
IO_lc(IO_stringValue(line,positions,j)) == 'velgrad'.or.& IO_lc(IO_stringValue(line,positions,i)) == 'velgrad'.or.&
IO_lc(IO_stringValue(line,positions,j)) == 'velocitygradient') IO_lc(IO_stringValue(line,positions,i)) == 'velocitygradient')
temp_valueVector = 0.0_pReal temp_valueVector = 0.0_pReal
temp_maskVector = .false. temp_maskVector = .false.
forall (k = 1_pInt:9_pInt) temp_maskVector(k) = IO_stringValue(line,positions,j+k) /= '*' forall (j = 1_pInt:9_pInt) temp_maskVector(j) = IO_stringValue(line,positions,i+j) /= '*'
do k = 1_pInt,9_pInt do j = 1_pInt,9_pInt
if (temp_maskVector(k)) temp_valueVector(k) = IO_floatValue(line,positions,j+k) if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,positions,i+j)
enddo enddo
bc(loadcase)%maskDeformation = transpose(reshape(temp_maskVector,[ 3,3])) bc(loadcase)%maskDeformation = transpose(reshape(temp_maskVector,[ 3,3]))
bc(loadcase)%deformation = math_plain9to33(temp_valueVector) bc(loadcase)%deformation = math_plain9to33(temp_valueVector)
case('p','pk1','piolakirchhoff','stress') case('p','pk1','piolakirchhoff','stress')
temp_valueVector = 0.0_pReal temp_valueVector = 0.0_pReal
forall (k = 1_pInt:9_pInt) bc(loadcase)%maskStressVector(k) =& forall (j = 1_pInt:9_pInt) bc(loadcase)%maskStressVector(j) =&
IO_stringValue(line,positions,j+k) /= '*' IO_stringValue(line,positions,i+j) /= '*'
do k = 1_pInt,9_pInt do j = 1_pInt,9_pInt
if (bc(loadcase)%maskStressVector(k)) temp_valueVector(k) =& if (bc(loadcase)%maskStressVector(j)) temp_valueVector(j) =&
IO_floatValue(line,positions,j+k) ! assign values for the bc(loadcase)%stress matrix IO_floatValue(line,positions,i+j) ! assign values for the bc(loadcase)%stress matrix
enddo enddo
bc(loadcase)%maskStress = transpose(reshape(bc(loadcase)%maskStressVector,[ 3,3])) bc(loadcase)%maskStress = transpose(reshape(bc(loadcase)%maskStressVector,[ 3,3]))
bc(loadcase)%stress = math_plain9to33(temp_valueVector) bc(loadcase)%stress = math_plain9to33(temp_valueVector)
case('t','time','delta') ! increment time case('t','time','delta') ! increment time
bc(loadcase)%time = IO_floatValue(line,positions,j+1_pInt) bc(loadcase)%time = IO_floatValue(line,positions,i+1_pInt)
case('temp','temperature') ! starting temperature case('temp','temperature') ! starting temperature
bc(loadcase)%temperature = IO_floatValue(line,positions,j+1_pInt) bc(loadcase)%temperature = IO_floatValue(line,positions,i+1_pInt)
case('n','incs','increments','steps') ! number of increments case('n','incs','increments','steps') ! number of increments
bc(loadcase)%incs = IO_intValue(line,positions,j+1_pInt) bc(loadcase)%incs = IO_intValue(line,positions,i+1_pInt)
case('logincs','logincrements','logsteps') ! number of increments (switch to log time scaling) case('logincs','logincrements','logsteps') ! number of increments (switch to log time scaling)
bc(loadcase)%incs = IO_intValue(line,positions,j+1_pInt) bc(loadcase)%incs = IO_intValue(line,positions,i+1_pInt)
bc(loadcase)%logscale = 1_pInt bc(loadcase)%logscale = 1_pInt
case('f','freq','frequency','outputfreq') ! frequency of result writings case('f','freq','frequency','outputfreq') ! frequency of result writings
bc(loadcase)%outputfrequency = IO_intValue(line,positions,j+1_pInt) bc(loadcase)%outputfrequency = IO_intValue(line,positions,i+1_pInt)
case('r','restart','restartwrite') ! frequency of writing restart information case('r','restart','restartwrite') ! frequency of writing restart information
bc(loadcase)%restartfrequency = max(0_pInt,IO_intValue(line,positions,j+1_pInt)) bc(loadcase)%restartfrequency = max(0_pInt,IO_intValue(line,positions,i+1_pInt))
case('guessreset','dropguessing') case('guessreset','dropguessing')
bc(loadcase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory bc(loadcase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory
case('euler') ! rotation of loadcase given in euler angles case('euler') ! rotation of loadcase given in euler angles
p = 0_pInt ! assuming values given in radians p = 0_pInt ! assuming values given in radians
l = 1_pInt ! assuming keyword indicating degree/radians k = 1_pInt ! assuming keyword indicating degree/radians
select case (IO_lc(IO_stringValue(line,positions,j+1_pInt))) select case (IO_lc(IO_stringValue(line,positions,i+1_pInt)))
case('deg','degree') case('deg','degree')
p = 1_pInt ! for conversion from degree to radian p = 1_pInt ! for conversion from degree to radian
case('rad','radian') case('rad','radian')
case default case default
l = 0_pInt ! immediately reading in angles, assuming radians k = 0_pInt ! immediately reading in angles, assuming radians
end select end select
forall(k = 1_pInt:3_pInt) temp33_Real(k,1) = & forall(j = 1_pInt:3_pInt) temp33_Real(j,1) = &
IO_floatValue(line,positions,j+l+k) * real(p,pReal) * inRad IO_floatValue(line,positions,i+k+j) * real(p,pReal) * inRad
bc(loadcase)%rotation = math_EulerToR(temp33_Real(:,1)) bc(loadcase)%rotation = math_EulerToR(temp33_Real(:,1))
case('rotation','rot') ! assign values for the rotation of loadcase matrix case('rotation','rot') ! assign values for the rotation of loadcase matrix
temp_valueVector = 0.0_pReal temp_valueVector = 0.0_pReal
forall (k = 1_pInt:9_pInt) temp_valueVector(k) = IO_floatValue(line,positions,j+k) forall (j = 1_pInt:9_pInt) temp_valueVector(j) = IO_floatValue(line,positions,i+j)
bc(loadcase)%rotation = math_plain9to33(temp_valueVector) bc(loadcase)%rotation = math_plain9to33(temp_valueVector)
end select end select
enddo; enddo enddo; enddo
@ -175,16 +252,7 @@ program DAMASK_spectral
call CPFEM_initAll(bc(1)%temperature,1_pInt,1_pInt) call CPFEM_initAll(bc(1)%temperature,1_pInt,1_pInt)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! get resolution, dimension, homogenization and variables derived from resolution ! output of geometry information
res = mesh_spectral_getResolution()
geomdim = mesh_spectral_getDimension()
homog = mesh_spectral_getHomogenization()
res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
Npoints = res(1)*res(2)*res(3)
wgt = 1.0_pReal/real(Npoints, pReal)
!--------------------------------------------------------------------------------------------------
! output of geometry
write(6,'(a)') '' write(6,'(a)') ''
write(6,'(a)') '#############################################################' write(6,'(a)') '#############################################################'
write(6,'(a)') 'DAMASK spectral:' write(6,'(a)') 'DAMASK spectral:'
@ -193,9 +261,9 @@ program DAMASK_spectral
write(6,'(a)') '#############################################################' write(6,'(a)') '#############################################################'
write(6,'(a)') 'geometry file: ',trim(geometryFile) write(6,'(a)') 'geometry file: ',trim(geometryFile)
write(6,'(a)') '=============================================================' write(6,'(a)') '============================================================='
write(6,'(a,3(i12 ))') 'resolution a b c:', res write(6,'(a,3(i12 ))') 'resolution a b c:', mesh_spectral_getResolution()
write(6,'(a,3(f12.5))') 'dimension x y z:', geomdim write(6,'(a,3(f12.5))') 'dimension x y z:', mesh_spectral_getDimension()
write(6,'(a,i5)') 'homogenization: ',homog write(6,'(a,i5)') 'homogenization: ', mesh_spectral_getHomogenization()
write(6,'(a)') '#############################################################' write(6,'(a)') '#############################################################'
write(6,'(a)') 'loadcase file: ',trim(loadCaseFile) write(6,'(a)') 'loadcase file: ',trim(loadCaseFile)
@ -246,27 +314,50 @@ program DAMASK_spectral
if (bc(loadcase)%outputfrequency < 1_pInt) errorID = 836_pInt ! non-positive result frequency if (bc(loadcase)%outputfrequency < 1_pInt) errorID = 836_pInt ! non-positive result frequency
if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string)
enddo enddo
!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!
initres = solverInit('AL')
F_aim = initres%F_init
!!!!!!!!!!!!!!
!!!!!!!!!!!!!!
!--------------------------------------------------------------------------------------------------
! write header of output file
if (appendToOutFile) then
open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.spectralOut',&
form='UNFORMATTED', position='APPEND', status='OLD')
else
open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.spectralOut',&
form='UNFORMATTED',status='REPLACE')
write(538) 'load', trim(loadCaseFile)
write(538) 'workingdir', trim(getSolverWorkingDirectoryName())
write(538) 'geometry', trim(geometryFile)
write(538) 'resolution', mesh_spectral_getResolution()
write(538) 'dimension', mesh_spectral_getDimension()
write(538) 'materialpoint_sizeResults', materialpoint_sizeResults
write(538) 'loadcases', N_Loadcases
write(538) 'frequencies', bc(1:N_Loadcases)%outputfrequency ! one entry per loadcase
write(538) 'times', bc(1:N_Loadcases)%time ! one entry per loadcase
write(538) 'logscales', bc(1:N_Loadcases)%logscale
write(538) 'increments', bc(1:N_Loadcases)%incs ! one entry per loadcase
write(538) 'startingIncrement', restartInc - 1_pInt ! start with writing out the previous inc
write(538) 'eoh' ! end of header
write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! initial (non-deformed or read-in) results
if (debugGeneral) write(6,'(a)') 'Header of result file written out'
endif
!################################################################################################## !##################################################################################################
! Loop over loadcases defined in the loadcase file ! Loop over loadcases defined in the loadcase file
!################################################################################################## !##################################################################################################
do loadcase = 1_pInt, N_Loadcases do loadcase = 1_pInt, N_Loadcases
time0 = time ! loadcase start time time0 = time ! loadcase start time
if (bc(loadcase)%followFormerTrajectory .and. & if (bc(loadcase)%followFormerTrajectory) then
(restartInc < totalIncsCounter .or. &
restartInc > totalIncsCounter+bc(loadcase)%incs) ) then ! continue to guess along former trajectory where applicable
guessmode = 1.0_pReal guessmode = 1.0_pReal
else else
guessmode = 0.0_pReal ! change of load case, homogeneous guess for the first inc guessmode = 0.0_pReal ! change of load case, homogeneous guess for the first inc
endif endif
!--------------------------------------------------------------------------------------------------
! arrays for mixed boundary conditions
mask_defgrad = merge(ones,zeroes,bc(loadcase)%maskDeformation)
mask_stress = merge(ones,zeroes,bc(loadcase)%maskStress)
size_reduced = int(count(bc(loadcase)%maskStressVector), pInt)
allocate (c_reduced(size_reduced,size_reduced), source =0.0_pReal)
allocate (s_reduced(size_reduced,size_reduced), source =0.0_pReal)
!################################################################################################## !##################################################################################################
! loop oper incs defined in input file for current loadcase ! loop oper incs defined in input file for current loadcase
!################################################################################################## !##################################################################################################
@ -309,101 +400,70 @@ program DAMASK_spectral
+ deltaF_aim + deltaF_aim
F_aim_lastInc = temp33_Real F_aim_lastInc = temp33_Real
!--------------------------------------------------------------------------------------------------
! update local deformation gradient and coordinates
deltaF_aim = math_rotate_backward33(deltaF_aim,bc(loadcase)%rotation)
call
call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,bc(loadcase)%rotation),& ! calculate current coordinates
1.0_pReal,F_lastInc,coordinates)
!--------------------------------------------------------------------------------------------------
! calculate reduced compliance
if(size_reduced > 0_pInt) then ! calculate compliance in case stress BC is applied
C_lastInc = math_rotate_forward3333(C,bc(loadcase)%rotation) ! calculate stiffness from former inc
temp99_Real = math_Plain3333to99(C_lastInc)
k = 0_pInt ! build reduced stiffness
do n = 1_pInt,9_pInt
if(bc(loadcase)%maskStressVector(n)) then
k = k + 1_pInt
j = 0_pInt
do m = 1_pInt,9_pInt
if(bc(loadcase)%maskStressVector(m)) then
j = j + 1_pInt
c_reduced(k,j) = temp99_Real(n,m)
endif; enddo; endif; enddo
call math_invert(size_reduced, c_reduced, s_reduced, i, errmatinv) ! invert reduced stiffness
if(errmatinv) call IO_error(error_ID=400_pInt)
temp99_Real = 0.0_pReal ! build full compliance
k = 0_pInt
do n = 1_pInt,9_pInt
if(bc(loadcase)%maskStressVector(n)) then
k = k + 1_pInt
j = 0_pInt
do m = 1_pInt,9_pInt
if(bc(loadcase)%maskStressVector(m)) then
j = j + 1_pInt
temp99_Real(n,m) = s_reduced(k,j)
endif; enddo; endif; enddo
S_lastInc = (math_Plain99to3333(temp99_Real))
endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new increment ! report begin of new increment
write(6,'(a)') '##################################################################' write(6,'(a)') '##################################################################'
write(6,'(A,I5.5,A,es12.5)') 'Increment ', totalIncsCounter, ' Time ',time write(6,'(A,I5.5,A,es12.5)') 'Increment ', totalIncsCounter, ' Time ',time
guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase
iter = 0_pInt solres = solution('AL')
err_div = huge(err_div_tol) ! go into loop !converged = solution(mySolver,ForwardFields(solver,deltaF_aim,timeinc/timeinc_old,guessmode, restartWrite))
converged = solution(mySolver,ForwardFields(solver,deltaF_aim,timeinc/timeinc_old,guessmode))
CPFEM_mode = 1_pInt ! winding forward
C = C * wgt
write(6,'(a)') '' write(6,'(a)') ''
write(6,'(a)') '==================================================================' write(6,'(a)') '=================================================================='
if(err_div > err_div_tol .or. err_stress > err_stress_tol) then if(solres%converged) then
write(6,'(A,I5.5,A)') 'increment ', totalIncsCounter, ' NOT converged'
notConvergedCounter = notConvergedCounter + 1_pInt
else
convergedCounter = convergedCounter + 1_pInt convergedCounter = convergedCounter + 1_pInt
write(6,'(A,I5.5,A)') 'increment ', totalIncsCounter, ' converged' write(6,'(A,I5.5,A)') 'increment ', totalIncsCounter, ' converged'
else
write(6,'(A,I5.5,A)') 'increment ', totalIncsCounter, ' NOT converged'
notConvergedCounter = notConvergedCounter + 1_pInt
endif endif
if (mod(inc,bc(loadcase)%outputFrequency) == 0_pInt) then ! at output frequency if (mod(inc,bc(loadcase)%outputFrequency) == 0_pInt) then ! at output frequency
write(6,'(a)') '' write(6,'(a)') ''
write(6,'(a)') '... writing results to file ......................................' write(6,'(a)') '... writing results to file ......................................'
write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! write result to file write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! write result to file
flush(538)
endif endif
if( bc(loadcase)%restartFrequency > 0_pInt .and. &
mod(inc,bc(loadcase)%restartFrequency) == 0_pInt) then ! at frequency of writing restart information set restart parameter for FEsolving (first call to CPFEM_general will write ToDo: true?)
restartInc=totalIncsCounter
restartWrite = .true.
write(6,'(a)') 'writing converged results for restart'
call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F)) ! writing deformation gradient field to file
write (777,rec=1) F
close (777)
call IO_write_jobBinaryFile(777,'C',size(C))
write (777,rec=1) C
close(777)
endif
endif ! end calculation/forwarding endif ! end calculation/forwarding
enddo ! end looping over incs in current loadcase guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase
deallocate(c_reduced)
deallocate(s_reduced) enddo ! end looping over incs in current loadcase
enddo ! end looping over loadcases enddo ! end looping over loadcases
write(6,'(a)') '' write(6,'(a)') ''
write(6,'(a)') '##################################################################' write(6,'(a)') '##################################################################'
write(6,'(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, &
' %) increments converged!' ' %) increments converged!'
close(538) close(538)
if (notConvergedCounter > 0_pInt) call quit(3_pInt) if (notConvergedCounter > 0_pInt) call quit(3_pInt)
call quit(0_pInt) call quit(0_pInt)
end program DAMASK_spectral
end program DAMASK_spectralDriver
subroutine quit(stop_id)
use prec, only: &
pInt
implicit none
integer(pInt), intent(in) :: stop_id
integer, dimension(8) :: dateAndTime ! type default integer
call date_and_time(values = dateAndTime)
write(6,'(/,a)') 'DAMASK terminated on:'
write(6,'(a,2(i2.2,a),i4.4)') 'Date: ',dateAndTime(3),'/',&
dateAndTime(2),'/',&
dateAndTime(1)
write(6,'(a,2(i2.2,a),i2.2)') 'Time: ',dateAndTime(5),':',&
dateAndTime(6),':',&
dateAndTime(7)
if (stop_id == 0_pInt) stop 0 ! normal termination
if (stop_id < 0_pInt) then ! trigger regridding
write(0,'(a,i6)') 'restart at ', stop_id*(-1_pInt)
stop 2
endif
if (stop_id == 3_pInt) stop 3 ! not all steps converged
stop 1 ! error (message from IO_error)
end subroutine

548
code/DAMASK_spectralSolver.f90
View File

@ -34,113 +34,36 @@
! !
! MPI fuer Eisenforschung, Duesseldorf ! MPI fuer Eisenforschung, Duesseldorf
#include "spectral_quit.f90" module DAMASK_spectralSolver
program DAMASK_spectral
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
use prec, only: pReal, pInt
use DAMASK_interface, only: & use math
use DAMASK_interface, only: &
DAMASK_interface_init, & DAMASK_interface_init, &
loadCaseFile, & loadCaseFile, &
geometryFile, & geometryFile, &
getSolverWorkingDirectoryName, & getSolverWorkingDirectoryName, &
getSolverJobName, & getSolverJobName, &
appendToOutFile appendToOutFile
use debug, only: &
use prec, only: &
pInt, &
pReal, &
DAMASK_NaN
use IO, only: &
IO_isBlank, &
IO_open_file, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_error, &
IO_lc, &
IO_read_jobBinaryFile, &
IO_write_jobBinaryFile
use debug, only: &
debug_level, & debug_level, &
debug_spectral, & debug_spectral, &
debug_levelBasic, & debug_levelBasic, &
debug_spectralDivergence, &
debug_spectralRestart, & debug_spectralRestart, &
debug_spectralFFTW, & debug_spectralFFTW
debug_reset, & use IO
debug_info
use math
use mesh, only : &
mesh_spectral_getResolution, &
mesh_spectral_getDimension, &
mesh_spectral_getHomogenization
use CPFEM, only: &
CPFEM_general, &
CPFEM_initAll
use FEsolving, only: &
restartWrite, &
restartInc
use numerics, only: &
err_div_tol, &
err_stress_tolrel, &
err_stress_tolabs, &
rotation_tol, &
itmax,&
itmin, &
memory_efficient, &
divergence_correction, &
DAMASK_NumThreadsInt, &
fftw_planner_flag, &
fftw_timelimit
use homogenization, only: &
materialpoint_sizeResults, &
materialpoint_results
implicit none implicit none
type solution_t ! mask of stress boundary conditions
#ifdef PETSC logical :: converged = .false.
#include <finclude/petscsys.h> logical :: regrid = .false.
#include <finclude/petscvec.h> end type solution_t
#include <finclude/petscsnes.h>
#include <finclude/petscvec.h90>
#include <finclude/petscsnes.h90>
#endif
!--------------------------------------------------------------------------------------------------
! variables related to information from load case and geom file
real(pReal), dimension(9) :: &
temp_valueVector !> temporarily from loadcase file when reading in tensors
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
integer(pInt), dimension(1_pInt + maxNchunksLoadcase*2_pInt) :: positions ! this is longer than needed for geometry parsing
integer(pInt) :: & type init
N_l = 0_pInt, & real(pReal), dimension(3,3) :: F_init
N_t = 0_pInt, & end type init
N_n = 0_pInt, &
N_Fdot = 0_pInt, &
Npoints,& ! number of Fourier points
homog, & ! homogenization scheme used
res1_red ! to store res(1)/2 +1
character(len=1024) :: & type bc_type
line
type bc_type
real(pReal), dimension (3,3) :: deformation = 0.0_pReal, & ! applied velocity gradient or time derivative of deformation gradient real(pReal), dimension (3,3) :: deformation = 0.0_pReal, & ! applied velocity gradient or time derivative of deformation gradient
stress = 0.0_pReal, & ! stress BC (if applicable) stress = 0.0_pReal, & ! stress BC (if applicable)
rotation = math_I3 ! rotation of BC (if applicable) rotation = math_I3 ! rotation of BC (if applicable)
@ -156,15 +79,74 @@ program DAMASK_spectral
maskStress = .false. ! mask of stress boundary conditions maskStress = .false. ! mask of stress boundary conditions
logical, dimension(9) :: maskStressVector = .false. ! linear mask of boundary conditions logical, dimension(9) :: maskStressVector = .false. ! linear mask of boundary conditions
end type end type
real(pReal), dimension(:,:,:,:,:), allocatable :: F, F_lastInc
real(pReal), dimension(:,:,:,:), allocatable :: coordinates
real(pReal), dimension(:,:,:), allocatable :: temperature
type(bc_type), allocatable, dimension(:) :: bc real(pReal), dimension(:,:,:,:,:), pointer :: P_real, deltaF_real ! field in real space (pointer)
complex(pReal), dimension(:,:,:,:,:), pointer :: P_fourier,deltaF_fourier ! field in fourier space (pointer)
complex(pReal), dimension(:,:,:), pointer :: scalarField_real
complex(pReal), dimension(:,:,:), pointer :: scalarField_fourier
!--------------------------------------------------------------------------------------------------
! variables for additional output of divergence calculations
type(C_PTR) :: divergence, plan_divergence, plan_correction
real(pReal), dimension(:,:,:,:), pointer :: divergence_real
complex(pReal), dimension(:,:,:,:), pointer :: divergence_fourier
real(pReal), dimension(:,:,:,:), allocatable :: divergence_post
contains
type(init) function solverInit(solver,restartInc,loadcase)
real(pReal) :: wgt use mesh, only : &
real(pReal), dimension(3) :: geomdim = 0.0_pReal, virt_dim = 0.0_pReal ! physical dimension of volume element per direction mesh_spectral_getResolution, &
integer(pInt), dimension(3) :: res = 1_pInt ! resolution (number of Fourier points) in each direction mesh_spectral_getDimension
use CPFEM, only: &
CPFEM_general
!-------------------------------------------------------------------------------------------------- use numerics, only: &
memory_efficient, &
divergence_correction, &
DAMASK_NumThreadsInt, &
fftw_planner_flag, &
fftw_timelimit
use debug, only: &
debug_level, &
debug_spectral, &
debug_levelBasic, &
debug_spectralDivergence, &
debug_spectralRestart, &
debug_spectralFFTW, &
debug_reset, &
debug_info
implicit none
real(pReal) :: restartInc
character(len=*) :: solver
integer(pInt) :: &
Npoints,& ! number of Fourier points
homog, & ! homogenization scheme used
res1_red
!--------------------------------------------------------------------------------------------------
! loop variables, convergence etc.
real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal
complex(pReal), dimension(3) :: temp3_Complex
complex(pReal), dimension(3,3) :: temp33_Complex
real(pReal), dimension(3,3) :: temp33_Real
integer(pInt) :: i, j, k, l, m, n, p, errorID
integer(pInt) :: inc, iter, ielem, CPFEM_mode=1_pInt, &
ierr
logical :: errmatinv
real(pReal) :: defgradDet
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc. ! stress, stiffness and compliance average etc.
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
P_av, & P_av, &
@ -176,71 +158,46 @@ program DAMASK_spectral
F_aim_lab, & F_aim_lab, &
F_aim_lab_lastIter, & F_aim_lab_lastIter, &
P_av_lab P_av_lab
real(pReal), dimension(3,3,3,3) :: & real(pReal), dimension(3,3,3,3) :: &
dPdF, & dPdF, &
C_ref = 0.0_pReal, & C_ref = 0.0_pReal, &
C = 0.0_pReal, & C = 0.0_pReal, &
S_lastInc, & S_lastInc, &
C_lastInc ! stiffness and compliance C_lastInc
real(pReal), dimension(6) :: sigma ! cauchy stress
real(pReal), dimension(6,6) :: dsde
real(pReal), dimension(9,9) :: temp99_Real ! compliance and stiffness in matrix notation
real(pReal), dimension(:,:), allocatable :: s_reduced, c_reduced ! reduced compliance and stiffness (only for stress BC)
integer(pInt) :: size_reduced = 0_pInt ! number of stress BCs
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! pointwise data ! pointwise data
type(C_PTR) :: tensorField ! field in real an fourier space type(C_PTR) :: tensorField ! field in real an fourier space
real(pReal), dimension(:,:,:,:,:), pointer :: P_real, deltaF_real ! field in real space (pointer) type(bc_type) :: loadcase ! field in real an fourier space
complex(pReal), dimension(:,:,:,:,:), pointer :: P_fourier,deltaF_fourier ! field in fourier space (pointer)
real(pReal), dimension(:,:,:,:,:), allocatable :: F, F_lastInc
real(pReal), dimension(:,:,:,:), allocatable :: coordinates
real(pReal), dimension(:,:,:), allocatable :: temperature
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! variables storing information for spectral method and FFTW ! variables storing information for spectral method and FFTW
type(C_PTR) :: plan_stress, plan_correction ! plans for fftw type(C_PTR) :: plan_stress, plan_backward ! plans for fftw
real(pReal), dimension(3,3) :: xiDyad ! product of wave vectors real(pReal), dimension(3,3) :: xiDyad ! product of wave vectors
real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat ! gamma operator (field) for spectral method real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat ! gamma operator (field) for spectral method
real(pReal), dimension(:,:,:,:), allocatable :: xi ! wave vector field for divergence and for gamma operator real(pReal), dimension(:,:,:,:), allocatable :: xi ! wave vector field for divergence and for gamma operator
integer(pInt), dimension(3) :: k_s integer(pInt), dimension(3) :: k_s
real(pReal), dimension(6) :: sigma ! cauchy stress
real(pReal), dimension(6,6) :: dsde
real(pReal) :: wgt
real(pReal), dimension(3) :: geomdim = 0.0_pReal, virt_dim = 0.0_pReal ! physical dimension of volume element per direction
integer(pInt), dimension(3) :: res = 1_pInt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! loop variables, convergence etc. ! pointwise data ! field in real an fourier space
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 real(pReal), dimension(:,:,:,:,:), pointer :: P_real, deltaF_real ! field in real space (pointer)
real(pReal) :: guessmode, err_div, err_stress, err_stress_tol
real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal complex(pReal), dimension(:,:,:,:,:), pointer :: P_fourier,deltaF_fourier ! field in fourier space (pointer)
complex(pReal), dimension(3) :: temp3_Complex
complex(pReal), dimension(3,3) :: temp33_Complex
real(pReal), dimension(3,3) :: temp33_Real
integer(pInt) :: i, j, k, l, m, n, p, errorID
integer(pInt) :: N_Loadcases, loadcase = 0_pInt, inc, iter, ielem, CPFEM_mode=1_pInt, &
ierr, totalIncsCounter = 0_pInt,&
notConvergedCounter = 0_pInt, convergedCounter = 0_pInt
logical :: errmatinv
real(pReal) :: defgradDet
character(len=6) :: loadcase_string
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!variables controlling debugging !variables controlling debugging
logical :: debugGeneral, debugDivergence, debugRestart, debugFFTW logical :: debugGeneral, debugDivergence, debugRestart, debugFFTW
!--------------------------------------------------------------------------------------------------
!variables for additional output due to general debugging
real(pReal) :: defgradDetMax, defgradDetMin, maxCorrectionSym, maxCorrectionSkew
!--------------------------------------------------------------------------------------------------
! variables for additional output of divergence calculations
type(C_PTR) :: divergence, plan_divergence
real(pReal), dimension(:,:,:,:), pointer :: divergence_real
complex(pReal), dimension(:,:,:,:), pointer :: divergence_fourier
real(pReal), dimension(:,:,:,:), allocatable :: divergence_post
real(pReal) :: pstress_av_L2, err_div_RMS, err_real_div_RMS, err_post_div_RMS,&
err_div_max, err_real_div_max
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! variables for debugging fft using a scalar field ! variables for debugging fft using a scalar field
@ -250,16 +207,10 @@ program DAMASK_spectral
complex(pReal), dimension(:,:,:), pointer :: scalarField_fourier complex(pReal), dimension(:,:,:), pointer :: scalarField_fourier
integer(pInt) :: row, column integer(pInt) :: row, column
!################################################################################################## if (solver == 'AL') solverInit%F_init=1.0_pReal
! reading of information from load case file and geometry file
!##################################################################################################
#ifdef PETSC
integer :: ierr_psc
call PetscInitialize(PETSC_NULL_CHARACTER, ierr_psc)
#endif
call DAMASK_interface_init
write(6,'(a)') '' write(6,'(a)') ''
write(6,'(a)') ' <<<+- DAMASK_spectral init -+>>>' write(6,'(a)') ' <<<+- DAMASK_spectralSolver init -+>>>'
write(6,'(a)') ' $Id$' write(6,'(a)') ' $Id$'
#include "compilation_info.f90" #include "compilation_info.f90"
write(6,'(a)') '' write(6,'(a)') ''
@ -279,12 +230,8 @@ program DAMASK_spectral
allocate (F_lastInc ( res(1), res(2),res(3),3,3), source = 0.0_pReal) allocate (F_lastInc ( res(1), res(2),res(3),3,3), source = 0.0_pReal)
allocate (xi (3,res1_red,res(2),res(3)), source = 0.0_pReal) allocate (xi (3,res1_red,res(2),res(3)), source = 0.0_pReal)
allocate (coordinates( res(1), res(2),res(3),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 = bc(1)%temperature) ! start out isothermally allocate (temperature( res(1), res(2),res(3)), source = loadcase%temperature) ! start out isothermally
tensorField = fftw_alloc_complex(int(res1_red*res(2)*res(3)*9_pInt,C_SIZE_T)) ! allocate continous data using a C function, C_SIZE_T is of type integer(8) tensorField = fftw_alloc_complex(int(res1_red*res(2)*res(3)*9_pInt,C_SIZE_T)) ! allocate continous data using a C function, C_SIZE_T is of type integer(8)
call c_f_pointer(tensorField, P_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for a real representation on tensorField
call c_f_pointer(tensorField, deltaF_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for a real representation on tensorField
call c_f_pointer(tensorField, P_fourier, [ res1_red, res(2),res(3),3,3]) ! place a pointer for a complex representation on tensorField
call c_f_pointer(tensorField, deltaF_fourier, [ res1_red, res(2),res(3),3,3]) ! place a pointer for a complex representation on tensorField
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! general initialization of fftw (see manual on fftw.org for more details) ! general initialization of fftw (see manual on fftw.org for more details)
@ -423,121 +370,129 @@ program DAMASK_spectral
gamma_hat(1,1,1, 1:3,1:3,1:3,1:3) = 0.0_pReal ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 gamma_hat(1,1,1, 1:3,1:3,1:3,1:3) = 0.0_pReal ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
endif endif
!-------------------------------------------------------------------------------------------------- end function solverInit
! write header of output file
if (appendToOutFile) then
open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.spectralOut',&
form='UNFORMATTED', position='APPEND', status='OLD') type(solution_t) function solution(solver,load,restartWrite)
else
open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.spectralOut',& use numerics, only: &
form='UNFORMATTED',status='REPLACE') err_div_tol, &
write(538) 'load', trim(loadCaseFile) err_stress_tolrel, &
write(538) 'workingdir', trim(getSolverWorkingDirectoryName()) err_stress_tolabs, &
write(538) 'geometry', trim(geometryFile) rotation_tol, &
write(538) 'resolution', res itmax,&
write(538) 'dimension', geomdim itmin, &
write(538) 'materialpoint_sizeResults', materialpoint_sizeResults memory_efficient, &
write(538) 'loadcases', N_Loadcases divergence_correction, &
write(538) 'frequencies', bc(1:N_Loadcases)%outputfrequency ! one entry per loadcase DAMASK_NumThreadsInt, &
write(538) 'times', bc(1:N_Loadcases)%time ! one entry per loadcase fftw_planner_flag, &
write(538) 'logscales', bc(1:N_Loadcases)%logscale fftw_timelimit
write(538) 'increments', bc(1:N_Loadcases)%incs ! one entry per loadcase
write(538) 'startingIncrement', restartInc - 1_pInt ! start with writing out the previous inc !--------------------------------------------------------------------------------------------------
write(538) 'eoh' ! end of header
write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! initial (non-deformed or read-in) results
if (debugGeneral) write(6,'(a)') 'Header of result file written out'
endif
!##################################################################################################
! Loop over loadcases defined in the loadcase file
!##################################################################################################
do loadcase = 1_pInt, N_Loadcases
time0 = time ! loadcase start time
if (bc(loadcase)%followFormerTrajectory .and. &
(restartInc < totalIncsCounter .or. &
restartInc > totalIncsCounter+bc(loadcase)%incs) ) then ! continue to guess along former trajectory where applicable
guessmode = 1.0_pReal
else
guessmode = 0.0_pReal ! change of load case, homogeneous guess for the first inc
endif
!--------------------------------------------------------------------------------------------------
! arrays for mixed boundary conditions ! arrays for mixed boundary conditions
mask_defgrad = merge(ones,zeroes,bc(loadcase)%maskDeformation) logical restartWrite
mask_stress = merge(ones,zeroes,bc(loadcase)%maskStress) character(len=*) :: solver
size_reduced = int(count(bc(loadcase)%maskStressVector), pInt) type(bc_type) :: load
real(pReal) :: pstress_av_L2, err_div_RMS, err_real_div_RMS, err_post_div_RMS,&
err_div_max, err_real_div_max
real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal
complex(pReal), dimension(3) :: temp3_complex
complex(pReal), dimension(3,3) :: temp33_complex
integer(pInt) :: size_reduced = 0_pInt
real(pReal), dimension(:,:), allocatable :: s_reduced, c_reduced ! reduced compliance and stiffness (only for stress BC)
real(pReal), dimension(6) :: sigma ! cauchy stress
real(pReal), dimension(6,6) :: dsde
real(pReal), dimension(9,9) :: temp99_Real ! compliance and stiffness in matrix notation
integer(pInt) :: Npoints
!--------------------------------------------------------------------------------------------------
! pointwise data
type(C_PTR) :: tensorField ! field in real an fourier space
real(pReal), dimension(:,:,:,:,:), pointer :: P_real, deltaF_real ! field in real space (pointer)
complex(pReal), dimension(:,:,:,:,:), pointer :: P_fourier,deltaF_fourier ! field in fourier space (pointer)
!--------------------------------------------------------------------------------------------------
! variables storing information for spectral method and FFTW
type(C_PTR) :: plan_stress, plan_correction ! plans for fftw
real(pReal), dimension(3,3) :: xiDyad ! product of wave vectors
real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat ! gamma operator (field) for spectral method
real(pReal), dimension(:,:,:,:), allocatable :: xi ! wave vector field for divergence and for gamma operator
integer(pInt), dimension(3) :: k_s, res
!--------------------------------------------------------------------------------------------------
! loop variables, convergence etc.
real(pReal) :: time0 = 0.0_pReal, timeinc = 1.0_pReal, timeinc_old = 0.0_pReal ! elapsed time, begin of interval, time interval
real(pReal) :: guessmode, err_div, err_stress, err_stress_tol
real(pReal), dimension(3,3) :: temp33_Real
integer(pInt) :: i, j, k, l, m, n, p, errorID
integer(pInt) :: N_Loadcases, loadcase = 0_pInt, inc, iter, ielem, CPFEM_mode=1_pInt, &
ierr, totalIncsCounter = 0_pInt,&
notConvergedCounter = 0_pInt, convergedCounter = 0_pInt
logical :: errmatinv
real(pReal) :: defgradDet
character(len=6) :: loadcase_string
!--------------------------------------------------------------------------------------------------
!variables controlling debugging
logical :: debugGeneral, debugDivergence, debugRestart, debugFFTW
!--------------------------------------------------------------------------------------------------
!variables for additional output due to general debugging
real(pReal) :: defgradDetMax, defgradDetMin, maxCorrectionSym, maxCorrectionSkew
!--------------------------------------------------------------------------------------------------
! variables for additional output of divergence calculations
type(C_PTR) :: divergence, plan_divergence
type(C_PTR) :: scalarField_realC, scalarField_fourierC,&
plan_scalarField_forth, plan_scalarField_back
real(pReal), dimension(:,:,:,:), pointer :: divergence_real
complex(pReal), dimension(:,:,:,:), pointer :: divergence_fourier
real(pReal), dimension(:,:,:,:), allocatable :: divergence_post
integer(pInt) :: row, column
real(pReal), dimension(3,3) :: &
P_av, &
F_aim = math_I3, &
F_aim_lastInc = math_I3, &
mask_stress, &
mask_defgrad, &
deltaF_aim, &
F_aim_lab, &
F_aim_lab_lastIter, &
P_av_lab
real(pReal), dimension(3,3,3,3) :: &
dPdF, &
C_ref = 0.0_pReal, &
C = 0.0_pReal, &
S_lastInc, &
C_lastInc
real(pReal), dimension(3) :: geomdim = 0.0_pReal, virt_dim = 0.0_pReal
integer(pInt) :: &
res1_red
real(pReal) :: wgt
if (solver == 'AL') solution%converged=.true.
mask_defgrad = merge(ones,zeroes,load%maskDeformation)
mask_stress = merge(ones,zeroes,load%maskStress)
size_reduced = int(count(load%maskStressVector), pInt)
allocate (c_reduced(size_reduced,size_reduced), source =0.0_pReal) allocate (c_reduced(size_reduced,size_reduced), source =0.0_pReal)
allocate (s_reduced(size_reduced,size_reduced), source =0.0_pReal) allocate (s_reduced(size_reduced,size_reduced), source =0.0_pReal)
!################################################################################################## !--------------------------------------------------------------------------------------------------
! loop oper incs defined in input file for current loadcase
!##################################################################################################
do inc = 1_pInt, bc(loadcase)%incs
totalIncsCounter = totalIncsCounter + 1_pInt
!--------------------------------------------------------------------------------------------------
! forwarding time
timeinc_old = timeinc
if (bc(loadcase)%logscale == 0_pInt) then ! linear scale
timeinc = bc(loadcase)%time/bc(loadcase)%incs ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used
else
if (loadcase == 1_pInt) then ! 1st loadcase of logarithmic scale
if (inc == 1_pInt) then ! 1st inc of 1st loadcase of logarithmic scale
timeinc = bc(1)%time*(2.0_pReal**real( 1_pInt-bc(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd
else ! not-1st inc of 1st loadcase of logarithmic scale
timeinc = bc(1)%time*(2.0_pReal**real(inc-1_pInt-bc(1)%incs ,pReal))
endif
else ! not-1st loadcase of logarithmic scale
timeinc = time0 *( (1.0_pReal + bc(loadcase)%time/time0 )**(real( inc,pReal)/&
real(bc(loadcase)%incs ,pReal))&
-(1.0_pReal + bc(loadcase)%time/time0 )**(real( (inc-1_pInt),pReal)/&
real(bc(loadcase)%incs ,pReal)) )
endif
endif
time = time + timeinc
if(totalIncsCounter >= restartInc) then ! do calculations (otherwise just forwarding)
if (bc(loadcase)%velGradApplied) then ! calculate deltaF_aim from given L and current F
deltaF_aim = timeinc * mask_defgrad * math_mul33x33(bc(loadcase)%deformation, F_aim)
else ! deltaF_aim = fDot *timeinc where applicable
deltaF_aim = timeinc * mask_defgrad * bc(loadcase)%deformation
endif
!--------------------------------------------------------------------------------------------------
! winding forward of deformation aim in loadcase system
temp33_Real = F_aim
F_aim = F_aim &
+ guessmode * mask_stress * (F_aim - F_aim_lastInc)*timeinc/timeinc_old &
+ deltaF_aim
F_aim_lastInc = temp33_Real
!--------------------------------------------------------------------------------------------------
! update local deformation gradient and coordinates
deltaF_aim = math_rotate_backward33(deltaF_aim,bc(loadcase)%rotation)
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
temp33_Real = F(i,j,k,1:3,1:3)
F(i,j,k,1:3,1:3) = F(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon
+ guessmode * (F(i,j,k,1:3,1:3) - F_lastInc(i,j,k,1:3,1:3))& ! guessing...
*timeinc/timeinc_old &
+ (1.0_pReal-guessmode) * deltaF_aim ! if not guessing, use prescribed average deformation where applicable
F_lastInc(i,j,k,1:3,1:3) = temp33_Real
enddo; enddo; enddo
call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,bc(loadcase)%rotation),& ! calculate current coordinates
1.0_pReal,F_lastInc,coordinates)
!--------------------------------------------------------------------------------------------------
! calculate reduced compliance ! calculate reduced compliance
if(size_reduced > 0_pInt) then ! calculate compliance in case stress BC is applied if(size_reduced > 0_pInt) then ! calculate compliance in case stress BC is applied
C_lastInc = math_rotate_forward3333(C,bc(loadcase)%rotation) ! calculate stiffness from former inc C_lastInc = math_rotate_forward3333(C,load%rotation) ! calculate stiffness from former inc
temp99_Real = math_Plain3333to99(C_lastInc) temp99_Real = math_Plain3333to99(C_lastInc)
k = 0_pInt ! build reduced stiffness k = 0_pInt ! build reduced stiffness
do n = 1_pInt,9_pInt do n = 1_pInt,9_pInt
if(bc(loadcase)%maskStressVector(n)) then if(load%maskStressVector(n)) then
k = k + 1_pInt k = k + 1_pInt
j = 0_pInt j = 0_pInt
do m = 1_pInt,9_pInt do m = 1_pInt,9_pInt
if(bc(loadcase)%maskStressVector(m)) then if(load%maskStressVector(m)) then
j = j + 1_pInt j = j + 1_pInt
c_reduced(k,j) = temp99_Real(n,m) c_reduced(k,j) = temp99_Real(n,m)
endif; enddo; endif; enddo endif; enddo; endif; enddo
@ -546,29 +501,23 @@ program DAMASK_spectral
temp99_Real = 0.0_pReal ! build full compliance temp99_Real = 0.0_pReal ! build full compliance
k = 0_pInt k = 0_pInt
do n = 1_pInt,9_pInt do n = 1_pInt,9_pInt
if(bc(loadcase)%maskStressVector(n)) then if(load%maskStressVector(n)) then
k = k + 1_pInt k = k + 1_pInt
j = 0_pInt j = 0_pInt
do m = 1_pInt,9_pInt do m = 1_pInt,9_pInt
if(bc(loadcase)%maskStressVector(m)) then if(load%maskStressVector(m)) then
j = j + 1_pInt j = j + 1_pInt
temp99_Real(n,m) = s_reduced(k,j) temp99_Real(n,m) = s_reduced(k,j)
endif; enddo; endif; enddo endif; enddo; endif; enddo
S_lastInc = (math_Plain99to3333(temp99_Real)) S_lastInc = (math_Plain99to3333(temp99_Real))
endif endif
guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase
!--------------------------------------------------------------------------------------------------
! report begin of new increment
write(6,'(a)') '##################################################################'
write(6,'(A,I5.5,A,es12.5)') 'Increment ', totalIncsCounter, ' Time ',time
guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase
iter = 0_pInt iter = 0_pInt
err_div = huge(err_div_tol) ! go into loop err_div = huge(err_div_tol) ! go into loop
!##################################################################################################
!##################################################################################################
! convergence loop (looping over iterations) ! convergence loop (looping over iterations)
!################################################################################################## !##################################################################################################
do while((iter < itmax .and. (err_div > err_div_tol .or. err_stress > err_stress_tol))& do while((iter < itmax .and. (err_div > err_div_tol .or. err_stress > err_stress_tol))&
.or. iter < itmin) .or. iter < itmin)
iter = iter + 1_pInt iter = iter + 1_pInt
@ -577,14 +526,14 @@ program DAMASK_spectral
! report begin of new iteration ! report begin of new iteration
write(6,'(a)') '' write(6,'(a)') ''
write(6,'(a)') '==================================================================' write(6,'(a)') '=================================================================='
write(6,'(6(a,i6.6))') 'Loadcase ',loadcase,' Inc. ',inc,'/',bc(loadcase)%incs,& write(6,'(6(a,i6.6))') 'Loadcase ',loadcase,' Inc. ',inc,'/',load%incs,&
' @ Iter. ',itmin,' < ',iter,' < ',itmax ' @ Iter. ',itmin,' < ',iter,' < ',itmax
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)
write(6,'(a)') '' write(6,'(a)') ''
write(6,'(a)') '... update stress field P(F) .....................................' write(6,'(a)') '... update stress field P(F) .....................................'
if (restartWrite) write(6,'(a)') 'writing restart info for last increment' if (restartWrite) write(6,'(a)') 'writing restart info for last increment'
F_aim_lab_lastIter = math_rotate_backward33(F_aim,bc(loadcase)%rotation) F_aim_lab_lastIter = math_rotate_backward33(F_aim,load%rotation)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
ielem = 0_pInt ielem = 0_pInt
@ -635,7 +584,7 @@ program DAMASK_spectral
call fftw_execute_dft_r2c(plan_stress,P_real,P_fourier) call fftw_execute_dft_r2c(plan_stress,P_real,P_fourier)
P_av_lab = real(P_fourier(1,1,1,1:3,1:3),pReal)*wgt P_av_lab = real(P_fourier(1,1,1,1:3,1:3),pReal)*wgt
P_av = math_rotate_forward33(P_av_lab,bc(loadcase)%rotation) P_av = math_rotate_forward33(P_av_lab,load%rotation)
write (6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'Piola-Kirchhoff stress / MPa =',& write (6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'Piola-Kirchhoff stress / MPa =',&
math_transpose33(P_av)/1.e6_pReal math_transpose33(P_av)/1.e6_pReal
@ -666,19 +615,19 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
if(size_reduced > 0_pInt) then ! calculate stress BC if applied if(size_reduced > 0_pInt) then ! calculate stress BC if applied
err_stress = maxval(abs(mask_stress * (P_av - bc(loadcase)%stress))) ! maximum deviaton (tensor norm not applicable) err_stress = maxval(abs(mask_stress * (P_av - load%stress))) ! maximum deviaton (tensor norm not applicable)
err_stress_tol = min(maxval(abs(P_av)) * err_stress_tolrel,err_stress_tolabs) ! don't use any tensor norm for the relative criterion because the comparison should be coherent err_stress_tol = min(maxval(abs(P_av)) * err_stress_tolrel,err_stress_tolabs) ! don't use any tensor norm for the relative criterion because the comparison should be coherent
write(6,'(a)') '' write(6,'(a)') ''
write(6,'(a)') '... correcting deformation gradient to fulfill BCs ...............' write(6,'(a)') '... correcting deformation gradient to fulfill BCs ...............'
write(6,'(a,f6.2,a,es11.4,a)') 'error stress = ', err_stress/err_stress_tol, & write(6,'(a,f6.2,a,es11.4,a)') 'error stress = ', err_stress/err_stress_tol, &
' (',err_stress,' Pa)' ' (',err_stress,' Pa)'
F_aim = F_aim - math_mul3333xx33(S_lastInc, ((P_av - bc(loadcase)%stress))) ! residual on given stress components F_aim = F_aim - math_mul3333xx33(S_lastInc, ((P_av - load%stress))) ! residual on given stress components
write(6,'(a,1x,es11.4)')'determinant of new deformation = ',math_det33(F_aim) write(6,'(a,1x,es11.4)')'determinant of new deformation = ',math_det33(F_aim)
else else
err_stress_tol = +huge(1.0_pReal) err_stress_tol = +huge(1.0_pReal)
endif endif
F_aim_lab = math_rotate_backward33(F_aim,bc(loadcase)%rotation) ! boundary conditions from load frame into lab (Fourier) frame F_aim_lab = math_rotate_backward33(F_aim,load%rotation) ! boundary conditions from load frame into lab (Fourier) frame
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! actual spectral method ! actual spectral method
@ -870,26 +819,11 @@ program DAMASK_spectral
CPFEM_mode = 1_pInt ! winding forward CPFEM_mode = 1_pInt ! winding forward
C = C * wgt C = C * wgt
write(6,'(a)') ''
write(6,'(a)') '=================================================================='
if(err_div > err_div_tol .or. err_stress > err_stress_tol) then
write(6,'(A,I5.5,A)') 'increment ', totalIncsCounter, ' NOT converged'
notConvergedCounter = notConvergedCounter + 1_pInt
else
convergedCounter = convergedCounter + 1_pInt
write(6,'(A,I5.5,A)') 'increment ', totalIncsCounter, ' converged'
endif
if (mod(inc,bc(loadcase)%outputFrequency) == 0_pInt) then ! at output frequency
write(6,'(a)') ''
write(6,'(a)') '... writing results to file ......................................'
write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! write result to file
flush(538)
endif
if( bc(loadcase)%restartFrequency > 0_pInt .and. & if( load%restartFrequency > 0_pInt .and. &
mod(inc,bc(loadcase)%restartFrequency) == 0_pInt) then ! at frequency of writing restart information set restart parameter for FEsolving (first call to CPFEM_general will write ToDo: true?) mod(inc,load%restartFrequency) == 0_pInt) then ! at frequency of writing restart information set restart parameter for FEsolving (first call to CPFEM_general will write ToDo: true?)
restartInc=totalIncsCounter !restartInc=totalIncsCounter
restartWrite = .true. restartWrite = .true.
write(6,'(a)') 'writing converged results for restart' write(6,'(a)') 'writing converged results for restart'
call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F)) ! writing deformation gradient field to file call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F)) ! writing deformation gradient field to file
@ -900,25 +834,11 @@ program DAMASK_spectral
close(777) close(777)
endif endif
endif ! end calculation/forwarding
enddo ! end looping over incs in current loadcase
deallocate(c_reduced) deallocate(c_reduced)
deallocate(s_reduced) deallocate(s_reduced)
enddo ! end looping over loadcases end function solution
write(6,'(a)') ''
write(6,'(a)') '##################################################################'
write(6,'(i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', &
notConvergedCounter + convergedCounter, ' (', & end module DAMASK_spectralSolver
real(convergedCounter, pReal)/&
real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, &
' %) increments converged!'
close(538)
call fftw_destroy_plan(plan_stress); call fftw_destroy_plan(plan_correction)
if (debugDivergence) call fftw_destroy_plan(plan_divergence)
if (debugFFTW) then
call fftw_destroy_plan(plan_scalarField_forth)
call fftw_destroy_plan(plan_scalarField_back)
endif
if (notConvergedCounter > 0_pInt) call quit(3_pInt)
call quit(0_pInt)
end program DAMASK_spectral

100
code/DAMASK_spectralSolver_Basic.f90
View File

@ -107,108 +107,10 @@ program DAMASK_spectral
materialpoint_results materialpoint_results
implicit none implicit none
! field in real an fourier space
#ifdef PETSC
#include <finclude/petscsys.h>
#include <finclude/petscvec.h>
#include <finclude/petscsnes.h>
#include <finclude/petscvec.h90>
#include <finclude/petscsnes.h90>
#endif
!--------------------------------------------------------------------------------------------------
! variables related to information from load case and geom file
real(pReal), dimension(9) :: &
temp_valueVector !> temporarily from loadcase file when reading in tensors
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
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, &
Npoints,& ! number of Fourier points
homog, & ! homogenization scheme used
res1_red ! to store res(1)/2 +1
character(len=1024) :: &
line
type bc_type
real(pReal), dimension (3,3) :: deformation = 0.0_pReal, & ! applied velocity gradient or time derivative of deformation gradient
stress = 0.0_pReal, & ! stress BC (if applicable)
rotation = math_I3 ! rotation of BC (if applicable)
real(pReal) :: time = 0.0_pReal, & ! length of increment
temperature = 300.0_pReal ! isothermal starting conditions
integer(pInt) :: incs = 0_pInt, & ! number of increments
outputfrequency = 1_pInt, & ! frequency of result writes
restartfrequency = 0_pInt, & ! frequency of restart writes
logscale = 0_pInt ! linear/logaritmic time inc flag
logical :: followFormerTrajectory = .true., & ! follow trajectory of former loadcase
velGradApplied = .false. ! decide wether velocity gradient or fdot is given
logical, dimension(3,3) :: maskDeformation = .false., & ! mask of deformation boundary conditions
maskStress = .false. ! mask of stress boundary conditions
logical, dimension(9) :: maskStressVector = .false. ! linear mask of boundary conditions
end type
type(bc_type), allocatable, dimension(:) :: bc
real(pReal) :: wgt
real(pReal), dimension(3) :: geomdim = 0.0_pReal, virt_dim = 0.0_pReal ! physical dimension of volume element per direction
integer(pInt), dimension(3) :: res = 1_pInt ! resolution (number of Fourier points) in each direction
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
real(pReal), dimension(3,3) :: &
P_av, &
F_aim = math_I3, &
F_aim_lastInc = math_I3, &
mask_stress, &
mask_defgrad, &
deltaF_aim, &
F_aim_lab, &
F_aim_lab_lastIter, &
P_av_lab
real(pReal), dimension(3,3,3,3) :: &
dPdF, &
C_ref = 0.0_pReal, &
C = 0.0_pReal, &
S_lastInc, &
C_lastInc ! stiffness and compliance
real(pReal), dimension(6) :: sigma ! cauchy stress
real(pReal), dimension(6,6) :: dsde
real(pReal), dimension(9,9) :: temp99_Real ! compliance and stiffness in matrix notation
real(pReal), dimension(:,:), allocatable :: s_reduced, c_reduced ! reduced compliance and stiffness (only for stress BC)
integer(pInt) :: size_reduced = 0_pInt ! number of stress BCs
!--------------------------------------------------------------------------------------------------
! pointwise data
type(C_PTR) :: tensorField ! field in real an fourier space
real(pReal), dimension(:,:,:,:,:), pointer :: P_real, deltaF_real ! field in real space (pointer) real(pReal), dimension(:,:,:,:,:), pointer :: P_real, deltaF_real ! field in real space (pointer)
complex(pReal), dimension(:,:,:,:,:), pointer :: P_fourier,deltaF_fourier ! field in fourier space (pointer) complex(pReal), dimension(:,:,:,:,:), pointer :: P_fourier,deltaF_fourier ! field in fourier space (pointer)
real(pReal), dimension(:,:,:,:,:), allocatable :: F, F_lastInc
real(pReal), dimension(:,:,:,:), allocatable :: coordinates
real(pReal), dimension(:,:,:), allocatable :: temperature
!--------------------------------------------------------------------------------------------------
! variables storing information for spectral method and FFTW
type(C_PTR) :: plan_stress, plan_correction ! plans for fftw
real(pReal), dimension(3,3) :: xiDyad ! product of wave vectors
real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat ! gamma operator (field) for spectral method
real(pReal), dimension(:,:,:,:), allocatable :: xi ! wave vector field for divergence and for gamma operator
integer(pInt), dimension(3) :: k_s
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! loop variables, convergence etc. ! loop variables, convergence etc.
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 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