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

made new spectral solver the default (no need for SOLVER=NEW anymore) and added additional warnings switche to be used with gfortran >= 4.6 

added comments
This commit is contained in:
Martin Diehl 2012-12-15 16:21:10 +00:00
parent ec9e7908a6
commit 566b680319
4 changed files with 166 additions and 150 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

154
code/DAMASK_spectral_driver.f90
View File

@ -7,7 +7,8 @@
!> @brief Driver controlling inner and outer load case looping of the various spectral solvers !> @brief Driver controlling inner and outer load case looping of the various spectral solvers
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
program DAMASK_spectral_Driver program DAMASK_spectral_Driver
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 DAMASK_interface, only: & use DAMASK_interface, only: &
DAMASK_interface_init, & DAMASK_interface_init, &
loadCaseFile, & loadCaseFile, &
@ -17,8 +18,7 @@ program DAMASK_spectral_Driver
appendToOutFile appendToOutFile
use prec, only: & use prec, only: &
pInt, & pInt, &
pReal, & pReal
DAMASK_NaN
use IO, only: & use IO, only: &
IO_isBlank, & IO_isBlank, &
IO_open_file, & IO_open_file, &
@ -78,33 +78,49 @@ program DAMASK_spectral_Driver
! variables related to information from load case and geom file ! variables related to information from load case and geom file
real(pReal), dimension(9) :: temp_valueVector = 0.0_pReal !< temporarily from loadcase file when reading in tensors (initialize to 0.0) real(pReal), dimension(9) :: temp_valueVector = 0.0_pReal !< temporarily from loadcase file when reading in tensors (initialize to 0.0)
logical, dimension(9) :: temp_maskVector = .false. !< temporarily from loadcase file when reading in tensors logical, dimension(9) :: temp_maskVector = .false. !< temporarily from loadcase file when reading in tensors
integer(pInt), parameter :: maxNchunksLoadcase = (1_pInt + 9_pInt)*3_pInt +& ! deformation, rotation, and stress integer(pInt), parameter :: maxNchunks = (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 + 1_pInt)*5_pInt + & ! time, (log)incs, temp, restartfrequency, and outputfrequency
1_pInt ! dropguessing 1_pInt, & ! dropguessing
integer(pInt), dimension(1_pInt + maxNchunksLoadcase*2_pInt) :: positions ! this is longer than needed for geometry parsing myUnit = 234_pInt !< file unit, DAMASK IO does not support newunit feature
integer(pInt), dimension(1_pInt + maxNchunks*2_pInt) :: positions ! this is longer than needed for geometry parsing
integer(pInt) :: & integer(pInt) :: &
N_l = 0_pInt, & N_l = 0_pInt, &
N_t = 0_pInt, & N_t = 0_pInt, &
N_n = 0_pInt, & N_n = 0_pInt, &
N_Fdot = 0_pInt, & ! N_Fdot = 0_pInt
myUnit = 234_pInt
character(len=1024) :: & character(len=1024) :: &
line line
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! loop variables, convergence etc. ! loop variables, convergence etc.
real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeros = 0.0_pReal real(pReal), dimension(3,3), parameter :: &
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 ones = 1.0_pReal, &
logical :: guess zeros = 0.0_pReal
integer(pInt) :: i, j, k, l, errorID, cutBackLevel = 0_pInt, stepFraction = 0_pInt real(pReal) :: &
integer(pInt) :: currentLoadcase = 0_pInt, inc, & time = 0.0_pReal, & !< elapsed time
totalIncsCounter = 0_pInt,& time0 = 0.0_pReal, & !< begin of interval
notConvergedCounter = 0_pInt, convergedCounter = 0_pInt, & timeinc = 1.0_pReal, & !< current time interval
resUnit = 0_pInt, statUnit = 0_pInt, lastRestartWritten = 0_pInt timeinc_old = 0.0_pReal ! !< previous time interval
logical :: &
guess !< guess along former trajectory
integer(pInt) :: &
i, j, k, l, &
errorID, &
cutBackLevel = 0_pInt, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$
stepFraction = 0_pInt !< fraction of current time interval
integer(pInt) :: &
currentLoadcase = 0_pInt, & !< current load case
inc, & !< current increment in current load case
totalIncsCounter = 0_pInt, & !< total No. of increments
convergedCounter = 0_pInt, & !< No. of converged increments
notConvergedCounter = 0_pInt, & !< No. of non-converged increments
resUnit = 0_pInt, & !< file unit for results writing
statUnit = 0_pInt, & !< file unit for statistics output
lastRestartWritten = 0_pInt !< total increment No. at which last restart information was written
character(len=6) :: loadcase_string character(len=6) :: loadcase_string
character(len=1024) :: incInfo character(len=1024) :: incInfo !< string parsed to solution with information about current load case
type(tLoadCase), allocatable, dimension(:) :: loadCases type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases
type(tSolutionState) solres type(tSolutionState) solres
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -122,8 +138,8 @@ program DAMASK_spectral_Driver
do do
read(myUnit,'(a1024)',END = 100) line read(myUnit,'(a1024)',END = 100) line
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
positions = IO_stringPos(line,maxNchunksLoadcase) positions = IO_stringPos(line,maxNchunks)
do i = 1_pInt, maxNchunksLoadcase, 1_pInt ! reading compulsory parameters for loadcase do i = 1_pInt, maxNchunks, 1_pInt ! reading compulsory parameters for loadcase
select case (IO_lc(IO_stringValue(line,positions,i))) select case (IO_lc(IO_stringValue(line,positions,i)))
case('l','velocitygrad','velgrad','velocitygradient') case('l','velocitygrad','velgrad','velocitygradient')
N_l = N_l + 1_pInt N_l = N_l + 1_pInt
@ -139,7 +155,7 @@ program DAMASK_spectral_Driver
100 if ((N_l + N_Fdot /= N_n) .or. (N_n /= N_t)) & ! sanity check 100 if ((N_l + N_Fdot /= N_n) .or. (N_n /= N_t)) & ! sanity check
call IO_error(error_ID=837_pInt,ext_msg = trim(loadCaseFile)) ! error message for incomplete loadcase call IO_error(error_ID=837_pInt,ext_msg = trim(loadCaseFile)) ! error message for incomplete loadcase
allocate (loadCases(N_n)) allocate (loadCases(N_n)) ! array of load cases
loadCases%P%myType='p' loadCases%P%myType='p'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -149,8 +165,8 @@ program DAMASK_spectral_Driver
read(myUnit,'(a1024)',END = 101) line read(myUnit,'(a1024)',END = 101) line
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
currentLoadCase = currentLoadCase + 1_pInt currentLoadCase = currentLoadCase + 1_pInt
positions = IO_stringPos(line,maxNchunksLoadcase) positions = IO_stringPos(line,maxNchunks)
do i = 1_pInt,maxNchunksLoadcase do i = 1_pInt,maxNchunks
select case (IO_lc(IO_stringValue(line,positions,i))) 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
temp_valueVector = 0.0_pReal temp_valueVector = 0.0_pReal
@ -160,19 +176,20 @@ program DAMASK_spectral_Driver
else else
loadCases(currentLoadCase)%deformation%myType = 'l' loadCases(currentLoadCase)%deformation%myType = 'l'
endif endif
forall (j = 1_pInt:9_pInt) temp_maskVector(j) = IO_stringValue(line,positions,i+j) /= '*' forall (j = 1_pInt:9_pInt) temp_maskVector(j) = IO_stringValue(line,positions,i+j) /= '*' ! true if not a *
do j = 1_pInt,9_pInt do j = 1_pInt,9_pInt
if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,positions,i+j) if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,positions,i+j) ! read value where applicable
enddo enddo
loadCases(currentLoadCase)%deformation%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) loadCases(currentLoadCase)%deformation%maskLogical = & ! logical mask in 3x3 notation
loadCases(currentLoadCase)%deformation%maskFloat = merge(ones,zeros,& transpose(reshape(temp_maskVector,[ 3,3]))
loadCases(currentLoadCase)%deformation%maskLogical) loadCases(currentLoadCase)%deformation%maskFloat = & ! float (1.0/0.0) mask in 3x3 notation
loadCases(currentLoadCase)%deformation%values = math_plain9to33(temp_valueVector) merge(ones,zeros,loadCases(currentLoadCase)%deformation%maskLogical)
loadCases(currentLoadCase)%deformation%values = math_plain9to33(temp_valueVector) ! values in 3x3 notation
case('p','pk1','piolakirchhoff','stress') case('p','pk1','piolakirchhoff','stress')
temp_valueVector = 0.0_pReal temp_valueVector = 0.0_pReal
forall (j = 1_pInt:9_pInt) temp_maskVector(j) = IO_stringValue(line,positions,i+j) /= '*' forall (j = 1_pInt:9_pInt) temp_maskVector(j) = IO_stringValue(line,positions,i+j) /= '*' ! true if not a *
do j = 1_pInt,9_pInt do j = 1_pInt,9_pInt
if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,positions,i+j) if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,positions,i+j) ! read value where applicable
enddo enddo
loadCases(currentLoadCase)%P%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) loadCases(currentLoadCase)%P%maskLogical = transpose(reshape(temp_maskVector,[ 3,3]))
loadCases(currentLoadCase)%P%maskFloat = merge(ones,zeros,& loadCases(currentLoadCase)%P%maskFloat = merge(ones,zeros,&
@ -190,7 +207,8 @@ program DAMASK_spectral_Driver
case('f','freq','frequency','outputfreq') ! frequency of result writings case('f','freq','frequency','outputfreq') ! frequency of result writings
loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,positions,i+1_pInt) loadCases(currentLoadCase)%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
loadCases(currentLoadCase)%restartfrequency = max(0_pInt,IO_intValue(line,positions,i+1_pInt)) loadCases(currentLoadCase)%restartfrequency = &
max(0_pInt,IO_intValue(line,positions,i+1_pInt))
case('guessreset','dropguessing') case('guessreset','dropguessing')
loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory
case('euler') ! rotation of currentLoadCase given in euler angles case('euler') ! rotation of currentLoadCase given in euler angles
@ -222,7 +240,6 @@ program DAMASK_spectral_Driver
checkLoadcases: do currentLoadCase = 1_pInt, size(loadCases) checkLoadcases: do currentLoadCase = 1_pInt, size(loadCases)
write (loadcase_string, '(i6)' ) currentLoadCase write (loadcase_string, '(i6)' ) currentLoadCase
write(6,'(1x,a,i6)') 'load case: ', currentLoadCase write(6,'(1x,a,i6)') 'load case: ', currentLoadCase
if (.not. loadCases(currentLoadCase)%followFormerTrajectory) & if (.not. loadCases(currentLoadCase)%followFormerTrajectory) &
write(6,'(2x,a)') 'drop guessing along trajectory' write(6,'(2x,a)') 'drop guessing along trajectory'
if (loadCases(currentLoadCase)%deformation%myType=='l') then if (loadCases(currentLoadCase)%deformation%myType=='l') then
@ -237,7 +254,7 @@ program DAMASK_spectral_Driver
endif endif
write(6,'(3(3(3x,f12.7,1x)/))',advance='no') & write(6,'(3(3(3x,f12.7,1x)/))',advance='no') &
merge(math_transpose33(loadCases(currentLoadCase)%deformation%values), & merge(math_transpose33(loadCases(currentLoadCase)%deformation%values), &
reshape(spread(huge(1.0_pReal),1,9),[ 3,3]), & reshape(spread(huge(1.0_pReal),1,9),[ 3,3]), & ! print *** (huge) for undefined
transpose(loadCases(currentLoadCase)%deformation%maskLogical)) transpose(loadCases(currentLoadCase)%deformation%maskLogical))
if (any(loadCases(currentLoadCase)%P%maskLogical .eqv. & if (any(loadCases(currentLoadCase)%P%maskLogical .eqv. &
loadCases(currentLoadCase)%deformation%maskLogical)) errorID = 831_pInt ! exclusive or masking only loadCases(currentLoadCase)%deformation%maskLogical)) errorID = 831_pInt ! exclusive or masking only
@ -267,9 +284,11 @@ program DAMASK_spectral_Driver
loadCases(currentLoadCase)%outputfrequency loadCases(currentLoadCase)%outputfrequency
write(6,'(2x,a,i5,/)') 'restart frequency: ', & write(6,'(2x,a,i5,/)') 'restart frequency: ', &
loadCases(currentLoadCase)%restartfrequency loadCases(currentLoadCase)%restartfrequency
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) ! exit with error message
enddo checkLoadcases enddo checkLoadcases
!--------------------------------------------------------------------------------------------------
! doing initialization depending on selected solver
select case (myspectralsolver) select case (myspectralsolver)
case (DAMASK_spectral_SolverBasic_label) case (DAMASK_spectral_SolverBasic_label)
call basic_init() call basic_init()
@ -285,15 +304,15 @@ program DAMASK_spectral_Driver
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! write header of output file ! write header of output file
if (appendToOutFile) then if (appendToOutFile) then ! after restart, append to existing results file
open(newunit=resUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& open(newunit=resUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
'.spectralOut',form='UNFORMATTED', position='APPEND', status='OLD') '.spectralOut',form='UNFORMATTED', position='APPEND', status='OLD')
open(newunit=statUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& open(newunit=statUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
'.sta',form='FORMATTED', position='APPEND', status='OLD') '.sta',form='FORMATTED', position='APPEND', status='OLD')
else else ! open new files ...
open(newunit=resUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& open(newunit=resUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
'.spectralOut',form='UNFORMATTED',status='REPLACE') '.spectralOut',form='UNFORMATTED',status='REPLACE')
write(resUnit) 'load', trim(loadCaseFile) write(resUnit) 'load', trim(loadCaseFile) ! ... and write header
write(resUnit) 'workingdir', trim(getSolverWorkingDirectoryName()) write(resUnit) 'workingdir', trim(getSolverWorkingDirectoryName())
write(resUnit) 'geometry', trim(geometryFile) write(resUnit) 'geometry', trim(geometryFile)
write(resUnit) 'resolution', res write(resUnit) 'resolution', res
@ -309,7 +328,7 @@ program DAMASK_spectral_Driver
write(resUnit) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:mesh_NcpElems) ! initial (non-deformed or read-in) results write(resUnit) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:mesh_NcpElems) ! initial (non-deformed or read-in) results
open(newunit=statUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& open(newunit=statUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
'.sta',form='FORMATTED',status='REPLACE') '.sta',form='FORMATTED',status='REPLACE')
write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file
if (debugGeneral) write(6,'(a)') 'Header of result file written out' if (debugGeneral) write(6,'(a)') 'Header of result file written out'
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -340,21 +359,24 @@ program DAMASK_spectral_Driver
timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1_pInt-loadCases(1)%incs ,pReal)) timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1_pInt-loadCases(1)%incs ,pReal))
endif endif
else ! not-1st currentLoadCase of logarithmic scale else ! not-1st currentLoadCase of logarithmic scale
timeinc = time0 *( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc,pReal)/& timeinc = time0 * &
real(loadCases(currentLoadCase)%incs ,pReal))& ( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc,pReal)/&
-(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( (inc-1_pInt),pReal)/& real(loadCases(currentLoadCase)%incs ,pReal))&
real(loadCases(currentLoadCase)%incs ,pReal)) ) -(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( (inc-1_pInt),pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal)))
endif endif
endif endif
timeinc = timeinc / 2.0_pReal**real(cutBackLevel,pReal) timeinc = timeinc / 2.0_pReal**real(cutBackLevel,pReal) ! depending on cut back level, decrease time step
if(totalIncsCounter >= restartInc) then ! do calculations (otherwise just forwarding) if(totalIncsCounter >= restartInc) then ! do calculations (otherwise just forwarding)
stepFraction = 0_pInt stepFraction = 0_pInt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new increment ! loop over sub incs
subIncLooping: do while (stepFraction/2_pInt**cutBackLevel <1_pInt) subIncLooping: do while (stepFraction/2_pInt**cutBackLevel <1_pInt)
time = time + timeinc time = time + timeinc ! forward time
stepFraction = stepFraction + 1_pInt stepFraction = stepFraction + 1_pInt
!--------------------------------------------------------------------------------------------------
! report begin of new increment
write(6,'(1/,a)') '###########################################################################' write(6,'(1/,a)') '###########################################################################'
write(6,'(a,es12.5'//& write(6,'(a,es12.5'//&
',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//& ',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//&
@ -370,6 +392,8 @@ program DAMASK_spectral_Driver
'-',stepFraction, '/', 2_pInt**cutBackLevel '-',stepFraction, '/', 2_pInt**cutBackLevel
select case(myspectralsolver) select case(myspectralsolver)
!--------------------------------------------------------------------------------------------------
! calculate solution
case (DAMASK_spectral_SolverBasic_label) case (DAMASK_spectral_SolverBasic_label)
solres = basic_solution (& solres = basic_solution (&
incInfo, guess,timeinc,timeinc_old, & incInfo, guess,timeinc,timeinc_old, &
@ -395,50 +419,52 @@ program DAMASK_spectral_Driver
rotation_BC = loadCases(currentLoadCase)%rotation) rotation_BC = loadCases(currentLoadCase)%rotation)
#endif #endif
end select end select
!--------------------------------------------------------------------------------------------------
! check solution
cutBack = .False. cutBack = .False.
if(solres%termIll .or. .not. solres%converged) then ! no solution found if(solres%termIll .or. .not. solres%converged) then ! no solution found
if (cutBackLevel < maxCutBack) then ! do cut back if (cutBackLevel < maxCutBack) then ! do cut back
write(6,'(/,a)') 'cut back detected' write(6,'(/,a)') 'cut back detected'
cutBack = .True. cutBack = .True.
stepFraction = (stepFraction - 1_pInt) * 2_pInt**cutBackLevel stepFraction = (stepFraction - 1_pInt) * 2_pInt**cutBackLevel ! adjust to new denominator
cutBackLevel = cutBackLevel + 1_pInt cutBackLevel = cutBackLevel + 1_pInt
time = time - timeinc time = time - timeinc ! rewind time
timeinc_old = timeinc timeinc_old = timeinc
timeinc = timeinc/2.0_pReal timeinc = timeinc/2.0_pReal
elseif (solres%termIll) then ! material point model cannot find a solution elseif (solres%termIll) then ! material point model cannot find a solution
if(regridMode > 0_pInt) call quit(-1*(lastRestartWritten+1)) if(regridMode > 0_pInt) call quit(-1*(lastRestartWritten+1)) ! regrid requested (mode 1 or 2)
call IO_error(850_pInt) call IO_error(850_pInt) ! no regrid (give up)
else else
if(regridMode == 2_pInt) call quit(-1*(lastRestartWritten+1)) if(regridMode == 2_pInt) call quit(-1*(lastRestartWritten+1)) ! regrid also if BVP solver do not converge
guess = .true. ! start guessing after first accepted (not converged) (sub)inc guess = .true. ! continue from non-converged solution and start guessing after accepted (sub)inc
endif endif
else else
guess = .true. ! start guessing after first converged (sub)inc guess = .true. ! start guessing after first converged (sub)inc
endif endif
if(guess) & if(guess) & ! write statistics about accepted solution
write(statUnit,*) inc, time, cutBackLevel, solres%converged, solres%iterationsNeeded write(statUnit,*) inc, time, cutBackLevel, solres%converged, solres%iterationsNeeded
enddo subIncLooping enddo subIncLooping
cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half subincs next inc cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half number of subincs next inc
if(solres%converged) then if(solres%converged) then ! report converged inc
convergedCounter = convergedCounter + 1_pInt convergedCounter = convergedCounter + 1_pInt
write(6,'(A,'//IO_intOut(totalIncsCounter)//',A)') & write(6,'(A,'//IO_intOut(totalIncsCounter)//',A)') &
'increment ', totalIncsCounter, ' converged' 'increment ', totalIncsCounter, ' converged'
else else
write(6,'(A,'//IO_intOut(totalIncsCounter)//',A)') & write(6,'(A,'//IO_intOut(totalIncsCounter)//',A)') & ! report non-converged inc
'increment ', totalIncsCounter, ' NOT converged' 'increment ', totalIncsCounter, ' NOT converged'
notConvergedCounter = notConvergedCounter + 1_pInt notConvergedCounter = notConvergedCounter + 1_pInt
endif endif
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency
write(6,'(1/,a)') '... writing results to file ......................................' write(6,'(1/,a)') '... writing results to file ......................................'
write(resUnit) materialpoint_results ! write result to file write(resUnit) materialpoint_results ! write result to file
endif endif
if( loadCases(currentLoadCase)%restartFrequency > 0_pInt .and. & if( loadCases(currentLoadCase)%restartFrequency > 0_pInt .and. &
mod(inc,loadCases(currentLoadCase)%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,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! at frequency of writing restart information set restart parameter for FEsolving (first call to CPFEM_general will write ToDo: true?)
restartWrite = .true. restartWrite = .true.
lastRestartWritten = inc lastRestartWritten = inc
endif endif
else !just time forwarding else !just time forwarding
time = time + timeinc time = time + timeinc
guess = .true. guess = .true.
endif ! end calculation/forwarding endif ! end calculation/forwarding
@ -459,6 +485,8 @@ program DAMASK_spectral_Driver
#endif #endif
end select end select
!--------------------------------------------------------------------------------------------------
! done report summary
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 ', &
@ -467,8 +495,8 @@ program DAMASK_spectral_Driver
real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, & real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, &
' %) increments converged!' ' %) increments converged!'
close(resUnit) close(resUnit)
if (notConvergedCounter > 0_pInt) call quit(3_pInt) if (notConvergedCounter > 0_pInt) call quit(3_pInt) ! error if some are not converged
call quit(0_pInt) call quit(0_pInt) ! no complains ;)
end program DAMASK_spectral_Driver end program DAMASK_spectral_Driver

4
code/DAMASK_spectral_utilities.f90
View File

@ -51,14 +51,14 @@ module DAMASK_spectral_utilities
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! derived types ! derived types
type tSolutionState !< return type of solution from spectral solver variants type tSolutionState !< return type of solution from spectral solver variants
logical :: converged = .true. logical :: converged = .true.
logical :: regrid = .false. logical :: regrid = .false.
logical :: termIll = .false. logical :: termIll = .false.
integer(pInt) :: iterationsNeeded = 0_pInt integer(pInt) :: iterationsNeeded = 0_pInt
end type tSolutionState end type tSolutionState
type tBoundaryCondition !< set of parameters defining a boundary condition type tBoundaryCondition !< set of parameters defining a boundary condition
real(pReal), dimension(3,3) :: values = 0.0_pReal real(pReal), dimension(3,3) :: values = 0.0_pReal
real(pReal), dimension(3,3) :: maskFloat = 0.0_pReal real(pReal), dimension(3,3) :: maskFloat = 0.0_pReal
logical, dimension(3,3) :: maskLogical = .false. logical, dimension(3,3) :: maskLogical = .false.

144
code/Makefile
View File

@ -9,17 +9,17 @@ SHELL = /bin/sh
# #
# Install fftw3 (v3.3 is tested): # Install fftw3 (v3.3 is tested):
# + execute # + execute
# ./configure --enable-threads --enable-sse2 --enable-shared [-enable-float] # ./configure --enable-threads --enable-sse2 --enable-shared
# make # make
# make install # make install
# + specify in the "pathinfo:FFTW" where FFTW was installed. # + specify in the "pathinfo:FFTW" where FFTW was installed.
# We essentially look for two library files "lib/libfftw3_threads" and "lib/libfftw3", # We essentially look for two library files "lib/libfftw3_threads" and "lib/libfftw3",
# so you can copy those, for instance, # so you can copy those, for instance,
# into DAMASK_ROOT/lib/fftw/lib/ and specify "./fftw/" as pathinfo:FFTW # into DAMASK_ROOT/lib/fftw/lib/ and specify "./fftw/" as pathinfo:FFTW
# Use --enable-float in above configure for single precision...
# Uses linux threads to parallelize fftw3 # Uses linux threads to parallelize fftw3
# #
# Instead of the AMD Core Math Library a standard "lib(64)/liblapack.a/dylib/etc." can be used. # Instead of the AMD Core Math Library or the Intel Kernel Math Library
# a standard "lib(64)/liblapack.a/dylib/etc." can be used.
# leave pathinfo:ACML and pathinfo:IKML blank, but specify as pathinfo:LAPACK where the library folder is located # leave pathinfo:ACML and pathinfo:IKML blank, but specify as pathinfo:LAPACK where the library folder is located
######################################################################################## ########################################################################################
# OPTIONS = standard (alternative): meaning # OPTIONS = standard (alternative): meaning
@ -111,13 +111,11 @@ endif
endif endif
endif endif
ifeq "$(SOLVER)" "NEW"
ifdef PETSC_DIR ifdef PETSC_DIR
include ${PETSC_DIR}/conf/variables include ${PETSC_DIR}/conf/variables
INCLUDE_DIRS +=${PETSC_FC_INCLUDES} -DPETSc INCLUDE_DIRS +=${PETSC_FC_INCLUDES} -DPETSc
LIBRARIES +=${PETSC_WITH_EXTERNAL_LIB} LIBRARIES +=${PETSC_WITH_EXTERNAL_LIB}
endif endif
endif
ifdef STANDARD_CHECK ifdef STANDARD_CHECK
STANDARD_CHECK_ifort =$(STANDARD_CHECK) STANDARD_CHECK_ifort =$(STANDARD_CHECK)
@ -127,9 +125,9 @@ endif
ifneq "$(FASTBUILD)" "YES" ifneq "$(FASTBUILD)" "YES"
STANDARD_CHECK_ifort ?=-stand f08 -standard-semantics -warn stderrors STANDARD_CHECK_ifort ?=-stand f08 -standard-semantics -warn stderrors
STANDARD_CHECK_gfortran ?=-std=f2008 -fall-intrinsics STANDARD_CHECK_gfortran ?=-std=f2008 -fall-intrinsics
#-std=2008ts for newer gfortran
endif endif
#-fall-intrinsics: all intrinsic procedures (including the GNU-specific extensions) are accepted. This can be useful with -std=f95 to force standard-compliance #-fall-intrinsics: all intrinsic procedures (including the GNU-specific extensions) are accepted. This can be useful with -std=f95 to force standard-compliance
#-std=f2008ts: for newer gfortran
# but get access to the full range of intrinsics available with gfortran. As a consequence, -Wintrinsics-std will be ignored and no user-defined # but get access to the full range of intrinsics available with gfortran. As a consequence, -Wintrinsics-std will be ignored and no user-defined
# procedure with the same name as any intrinsic will be called except when it is explicitly declared external # procedure with the same name as any intrinsic will be called except when it is explicitly declared external
#-Wintrinsics-std: warnings because of "flush" is not longer in the standard, but still an intrinsic fuction of the compilers #-Wintrinsics-std: warnings because of "flush" is not longer in the standard, but still an intrinsic fuction of the compilers
@ -221,7 +219,12 @@ COMPILE_OPTIONS_gfortran :=$(COMPILE_OPTIONS_gfortran)\
-Wunsafe-loop-optimizations\ -Wunsafe-loop-optimizations\
-Wunused\ -Wunused\
-Wall\ -Wall\
-Wextra -Wextra\
-Wsuggest-attribute=const\
-Wsuggest-attribute=pure\
-Wsuggest-attribute=noreturn\
-Wreal-q-constant\
-pedantic-errors
endif endif
#-xf95-cpp-input: preprocessor #-xf95-cpp-input: preprocessor
@ -246,13 +249,8 @@ endif
# -value: # -value:
# -parameter: find usused variables with "parameter" attribute # -parameter: find usused variables with "parameter" attribute
#-Wextra: #-Wextra:
#-Wreal-q-constant: warn about real-literal-constants with 'q' exponent-letter
################################################################################################### ###################################################################################################
#OPTIONS FOR GFORTRAN 4.6
#-Wsuggest-attribute=const:
#-Wsuggest-attribute=noreturn:
#-Wsuggest-attribute=pure:
#-Wreal-q-constant: Warn about real-literal-constants with 'q' exponent-letter
#
#MORE OPTIONS FOR DEBUGGING DURING COMPILING #MORE OPTIONS FOR DEBUGGING DURING COMPILING
#-Wline-truncation: too many warnings because we have comments beyond character 132 #-Wline-truncation: too many warnings because we have comments beyond character 132
#-Warray-temporarieswarnings: because we have many temporary arrays (performance issue?): #-Warray-temporarieswarnings: because we have many temporary arrays (performance issue?):
@ -286,15 +284,14 @@ COMPILED_FILES = prec.o DAMASK_spectral_interface.o IO.o numerics.o debug.o math
FEsolving.o mesh.o material.o lattice.o \ FEsolving.o mesh.o material.o lattice.o \
constitutive_dislotwin.o constitutive_j2.o constitutive_phenopowerlaw.o \ constitutive_dislotwin.o constitutive_j2.o constitutive_phenopowerlaw.o \
constitutive_titanmod.o constitutive_nonlocal.o constitutive_none.o constitutive.o crystallite.o \ constitutive_titanmod.o constitutive_nonlocal.o constitutive_none.o constitutive.o crystallite.o \
homogenization_RGC.o homogenization_isostrain.o homogenization.o CPFEM.o homogenization_RGC.o homogenization_isostrain.o homogenization.o CPFEM.o \
DAMASK_spectral_utilities.o DAMASK_spectral_solverBasic.o
ifeq "$(SOLVER)" "NEW"
ifdef PETSC_DIR ifdef PETSC_DIR
PETSC_FILES = DAMASK_spectral_solverAL.o DAMASK_spectral_solverBasicPETSc.o PETSC_FILES = DAMASK_spectral_solverAL.o DAMASK_spectral_solverBasicPETSc.o
COMPILED_FILES += $(PETSC_FILES)
endif endif
COMPILED_FILES += DAMASK_spectral_utilities.o DAMASK_spectral_solverBasic.o $(PETSC_FILES)
DAMASK_spectral.exe: DAMASK_spectral_driver.o DAMASK_spectral.exe: DAMASK_spectral_driver.o
$(PREFIX) $(COMPILERNAME) $(OPENMP_FLAG_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) $(STANDARD_CHECK_$(F90)) \ $(PREFIX) $(COMPILERNAME) $(OPENMP_FLAG_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) $(STANDARD_CHECK_$(F90)) \
@ -315,89 +312,80 @@ DAMASK_spectral_solverBasicPETSc.o: DAMASK_spectral_solverBasicPETSc.f90\
DAMASK_spectral_utilities.o: DAMASK_spectral_utilities.f90\ DAMASK_spectral_utilities.o: DAMASK_spectral_utilities.f90\
CPFEM.o CPFEM.o
else CPFEM.o: CPFEM.f90\
DAMASK_spectral.exe: DAMASK_spectral.o homogenization.o
$(PREFIX) $(COMPILERNAME) $(OPENMP_FLAG_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) $(STANDARD_CHECK_$(F90)) \
-o DAMASK_spectral.exe DAMASK_spectral.o \
$(COMPILED_FILES) $(LIB_DIRS) $(LIBRARIES) $(SUFFIX)
DAMASK_spectral.o: DAMASK_spectral.f90 CPFEM.o homogenization.o: homogenization.f90\
$(PREFIX) $(COMPILERNAME) $(COMPILE_MAXOPTI) -c DAMASK_spectral.f90 $(SUFFIX) homogenization_RGC.o \
endif homogenization_isostrain.o
CPFEM.o: CPFEM.f90\
homogenization.o
homogenization.o: homogenization.f90\ homogenization_RGC.o: homogenization_RGC.f90 \
homogenization_RGC.o \ crystallite.o
homogenization_isostrain.o
homogenization_RGC.o: homogenization_RGC.f90 \ homogenization_isostrain.o: homogenization_isostrain.f90 \
crystallite.o crystallite.o
homogenization_isostrain.o: homogenization_isostrain.f90 \ crystallite.o: crystallite.f90 \
crystallite.o constitutive.o
crystallite.o: crystallite.f90 \
constitutive.o
constitutive.o: constitutive.f90 \ constitutive.o: constitutive.f90 \
constitutive_nonlocal.o \ constitutive_nonlocal.o \
constitutive_titanmod.o \ constitutive_titanmod.o \
constitutive_dislotwin.o \ constitutive_dislotwin.o \
constitutive_phenopowerlaw.o \ constitutive_phenopowerlaw.o \
constitutive_j2.o \ constitutive_j2.o \
constitutive_none.o constitutive_none.o
constitutive_nonlocal.o: constitutive_nonlocal.f90 \ constitutive_nonlocal.o: constitutive_nonlocal.f90 \
lattice.o lattice.o
constitutive_titanmod.o: constitutive_titanmod.f90 \ constitutive_titanmod.o: constitutive_titanmod.f90 \
lattice.o lattice.o
constitutive_dislotwin.o: constitutive_dislotwin.f90 \ constitutive_dislotwin.o: constitutive_dislotwin.f90 \
lattice.o lattice.o
constitutive_phenopowerlaw.o: constitutive_phenopowerlaw.f90 \ constitutive_phenopowerlaw.o: constitutive_phenopowerlaw.f90 \
lattice.o lattice.o
constitutive_j2.o: constitutive_j2.f90 \ constitutive_j2.o: constitutive_j2.f90 \
lattice.o lattice.o
constitutive_none.o: constitutive_none.f90 \ constitutive_none.o: constitutive_none.f90 \
lattice.o lattice.o
lattice.o: lattice.f90 \ lattice.o: lattice.f90 \
material.o material.o
material.o: material.f90 \ material.o: material.f90 \
mesh.o mesh.o
mesh.o: mesh.f90 \ mesh.o: mesh.f90 \
FEsolving.o \ FEsolving.o \
math.o math.o
FEsolving.o: FEsolving.f90 \ FEsolving.o: FEsolving.f90 \
debug.o debug.o
math.o: math.f90 \ math.o: math.f90 \
debug.o debug.o
debug.o: debug.f90 \ debug.o: debug.f90 \
numerics.o numerics.o
numerics.o: numerics.f90 \ numerics.o: numerics.f90 \
IO.o IO.o
IO.o: IO.f90 \ IO.o: IO.f90 \
DAMASK_spectral_interface.o DAMASK_spectral_interface.o
DAMASK_spectral_interface.o: DAMASK_spectral_interface.f90 \ DAMASK_spectral_interface.o: DAMASK_spectral_interface.f90 \
prec.o prec.o
ifeq "$(F90)" "gfortran" # fno-range-check because NaN is defined in prec ifeq "$(F90)" "gfortran" # fno-range-check because NaN is defined in prec
prec.o: prec.f90 prec.o: prec.f90
$(PREFIX) $(COMPILERNAME) $(COMPILE) -c -fno-range-check prec.f90 $(SUFFIX) $(PREFIX) $(COMPILERNAME) $(COMPILE) -c -fno-range-check prec.f90 $(SUFFIX)
else else
prec.o: prec.f90 prec.o: prec.f90
$(PREFIX) $(COMPILERNAME) $(COMPILE) -c prec.f90 $(SUFFIX) $(PREFIX) $(COMPILERNAME) $(COMPILE) -c prec.f90 $(SUFFIX)
endif endif

14
code/numerics.f90
View File

@ -131,16 +131,16 @@ subroutine numerics_init
IO_floatValue, & IO_floatValue, &
IO_intValue, & IO_intValue, &
IO_warning IO_warning
#ifndef Marc ! Use the standard conforming module file for omp if using the spectral solver #ifndef Marc
!$ use OMP_LIB, only: omp_set_num_threads !$ use OMP_LIB, only: omp_set_num_threads ! Use the standard conforming module file for omp if not using MSC.Marc
#endif #endif
implicit none implicit none
#ifdef Marc ! use the non F90 standard include file because some versions of Marc and Abaqus crash when using the module #ifdef Marc
!$ include "omp_lib.h" !$ include "omp_lib.h" ! use the non F90 standard include file to prevent crashes with some versions of MSC.Marc
#endif #endif
#ifdef PETSc #ifdef PETSc
PetscErrorCode :: ierr PetscErrorCode :: ierr
#endif #endif
integer(pInt), parameter :: fileunit = 300_pInt ,& integer(pInt), parameter :: fileunit = 300_pInt ,&
maxNchunks = 2_pInt maxNchunks = 2_pInt
!$ integer :: gotDAMASK_NUM_THREADS = 1 !$ integer :: gotDAMASK_NUM_THREADS = 1