2012-07-30 19:36:22 +05:30
!--------------------------------------------------------------------------------------------------
2012-07-19 22:54:56 +05:30
!* $Id$
2012-07-30 19:36:22 +05:30
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Driver controlling inner and outer load case looping of the various spectral solvers
!--------------------------------------------------------------------------------------------------
2012-07-24 22:37:10 +05:30
program DAMASK_spectral_Driver
2012-07-20 21:03:13 +05:30
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
2012-07-19 22:54:56 +05:30
use prec , only : &
pInt , &
2012-07-20 21:03:13 +05:30
pReal , &
DAMASK_NaN
2012-07-19 22:54:56 +05:30
use IO , only : &
2012-07-20 21:03:13 +05:30
IO_isBlank , &
IO_open_file , &
IO_stringPos , &
IO_stringValue , &
IO_floatValue , &
IO_intValue , &
IO_error , &
IO_lc , &
IO_read_jobBinaryFile , &
2012-07-19 22:54:56 +05:30
IO_write_jobBinaryFile
use math
2012-07-20 21:03:13 +05:30
use mesh , only : &
2012-08-03 14:55:48 +05:30
res , &
geomdim , &
mesh_NcpElems
2012-07-20 21:03:13 +05:30
use CPFEM , only : &
CPFEM_initAll
2012-07-19 22:54:56 +05:30
use FEsolving , only : &
restartWrite , &
restartInc
2012-07-20 21:03:13 +05:30
use numerics , only : &
2012-07-25 19:31:39 +05:30
rotation_tol , &
2012-08-03 14:55:48 +05:30
mySpectralSolver
2012-07-20 21:03:13 +05:30
2012-07-19 22:54:56 +05:30
use homogenization , only : &
materialpoint_sizeResults , &
materialpoint_results
2012-08-03 14:55:48 +05:30
use DAMASK_spectral_Utilities , only : &
boundaryCondition , &
solutionState , &
debugGeneral
2012-07-24 22:37:10 +05:30
use DAMASK_spectral_SolverBasic
2012-08-14 22:28:23 +05:30
use DAMASK_spectral_SolverBasicPETSC
2012-08-06 14:23:12 +05:30
use DAMASK_spectral_SolverAL
2012-07-19 22:54:56 +05:30
implicit none
2012-07-20 21:03:13 +05:30
2012-08-03 14:55:48 +05:30
type loadCase
real ( pReal ) , dimension ( 3 , 3 ) :: rotation = math_I3 ! rotation of BC
type ( boundaryCondition ) :: P , & ! stress BC
deformation ! deformation BC (Fdot or L)
2012-07-23 15:42:31 +05:30
real ( pReal ) :: time = 0.0_pReal , & ! length of increment
temperature = 30 0.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
2012-08-03 14:55:48 +05:30
logical :: followFormerTrajectory = . true . ! follow trajectory of former loadcase
end type loadCase
2012-07-23 15:42:31 +05:30
2012-07-20 21:03:13 +05:30
!--------------------------------------------------------------------------------------------------
! variables related to information from load case and geom file
2012-07-30 19:36:22 +05:30
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
2012-07-20 21:03:13 +05:30
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 , &
2012-07-23 15:42:31 +05:30
N_Fdot = 0_pInt ! number of Fourier points
2012-07-20 21:03:13 +05:30
character ( len = 1024 ) :: &
line
2012-08-03 14:55:48 +05:30
2012-07-19 22:54:56 +05:30
!--------------------------------------------------------------------------------------------------
! loop variables, convergence etc.
2012-08-03 14:55:48 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , parameter :: ones = 1.0_pReal , zeroes = 0.0_pReal
real ( pReal ) :: time = 0.0_pReal , time0 = 0.0_pReal , timeinc = 1.0_pReal , timeinc_old = 0.0_pReal ! elapsed time, begin of interval, time interval, previous time interval
2012-07-20 21:03:13 +05:30
real ( pReal ) :: guessmode
2012-07-30 19:36:22 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: temp33_Real
2012-07-24 22:37:10 +05:30
integer ( pInt ) :: i , j , k , l , errorID
integer ( pInt ) :: currentLoadcase = 0_pInt , inc , &
2012-07-20 21:03:13 +05:30
totalIncsCounter = 0_pInt , &
notConvergedCounter = 0_pInt , convergedCounter = 0_pInt
2012-07-23 15:42:31 +05:30
character ( len = 6 ) :: loadcase_string
2012-08-03 14:55:48 +05:30
type ( loadCase ) , allocatable , dimension ( : ) :: loadCases
type ( solutionState ) solres
2012-07-19 22:54:56 +05:30
2012-08-03 14:55:48 +05:30
!--------------------------------------------------------------------------------------------------
! init DAMASK (all modules)
2012-07-30 19:36:22 +05:30
call CPFEM_initAll ( temperature = 30 0.0_pReal , element = 1_pInt , IP = 1_pInt )
2012-07-19 22:54:56 +05:30
write ( 6 , '(a)' ) ''
2012-07-25 19:31:39 +05:30
write ( 6 , '(a)' ) ' <<<+- DAMASK_spectral_Driver init -+>>>'
2012-07-19 22:54:56 +05:30
write ( 6 , '(a)' ) ' $Id$'
#include "compilation_info.f90"
write ( 6 , '(a)' ) ''
!--------------------------------------------------------------------------------------------------
2012-08-03 14:55:48 +05:30
! reading basic information from load case file and allocate data structure containing load cases
2012-07-19 22:54:56 +05:30
call IO_open_file ( myUnit , trim ( loadCaseFile ) )
rewind ( myUnit )
do
read ( myUnit , '(a1024)' , END = 100 ) line
if ( IO_isBlank ( line ) ) cycle ! skip empty lines
positions = IO_stringPos ( line , maxNchunksLoadcase )
do i = 1_pInt , maxNchunksLoadcase , 1_pInt ! reading compulsory parameters for loadcase
select case ( IO_lc ( IO_stringValue ( line , positions , i ) ) )
case ( 'l' , 'velocitygrad' , 'velgrad' , 'velocitygradient' )
N_l = N_l + 1_pInt
case ( 'fdot' , 'dotf' )
N_Fdot = N_Fdot + 1_pInt
case ( 't' , 'time' , 'delta' )
N_t = N_t + 1_pInt
case ( 'n' , 'incs' , 'increments' , 'steps' , 'logincs' , 'logincrements' , 'logsteps' )
N_n = N_n + 1_pInt
end select
enddo ! count all identifiers to allocate memory and do sanity check
enddo
2012-07-24 22:37:10 +05:30
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
2012-08-03 14:55:48 +05:30
allocate ( loadCases ( N_n ) )
loadCases % P % myType = 'p'
2012-07-19 22:54:56 +05:30
2012-08-03 14:55:48 +05:30
!--------------------------------------------------------------------------------------------------
2012-07-19 22:54:56 +05:30
! reading the load case and assign values to the allocated data structure
rewind ( myUnit )
do
read ( myUnit , '(a1024)' , END = 101 ) line
if ( IO_isBlank ( line ) ) cycle ! skip empty lines
2012-08-03 14:55:48 +05:30
currentLoadCase = currentLoadCase + 1_pInt
2012-07-19 22:54:56 +05:30
positions = IO_stringPos ( line , maxNchunksLoadcase )
2012-07-20 21:03:13 +05:30
do i = 1_pInt , maxNchunksLoadcase
select case ( IO_lc ( IO_stringValue ( line , positions , i ) ) )
2012-07-19 22:54:56 +05:30
case ( 'fdot' , 'dotf' , 'l' , 'velocitygrad' , 'velgrad' , 'velocitygradient' ) ! assign values for the deformation BC matrix
2012-08-03 14:55:48 +05:30
if ( IO_lc ( IO_stringValue ( line , positions , i ) ) == 'l' . or . & ! in case of given L, set flag to true
IO_lc ( IO_stringValue ( line , positions , i ) ) == 'velocitygrad' . or . &
IO_lc ( IO_stringValue ( line , positions , i ) ) == 'velgrad' . or . &
IO_lc ( IO_stringValue ( line , positions , i ) ) == 'velocitygradient' ) then
loadCases ( currentLoadCase ) % deformation % myType = 'l'
else
loadCases ( currentLoadCase ) % deformation % myType = 'fdot'
endif
2012-07-20 21:03:13 +05:30
forall ( j = 1_pInt : 9_pInt ) temp_maskVector ( j ) = IO_stringValue ( line , positions , i + j ) / = '*'
do j = 1_pInt , 9_pInt
if ( temp_maskVector ( j ) ) temp_valueVector ( j ) = IO_floatValue ( line , positions , i + j )
2012-07-19 22:54:56 +05:30
enddo
2012-08-03 14:55:48 +05:30
loadCases ( currentLoadCase ) % deformation % maskLogical = transpose ( reshape ( temp_maskVector , [ 3 , 3 ] ) )
loadCases ( currentLoadCase ) % deformation % maskFloat = merge ( ones , zeroes , &
loadCases ( currentLoadCase ) % deformation % maskLogical )
loadCases ( currentLoadCase ) % deformation % values = math_plain9to33 ( temp_valueVector )
2012-07-19 22:54:56 +05:30
case ( 'p' , 'pk1' , 'piolakirchhoff' , 'stress' )
temp_valueVector = 0.0_pReal
2012-08-03 14:55:48 +05:30
forall ( j = 1_pInt : 9_pInt ) temp_maskVector ( j ) = IO_stringValue ( line , positions , i + j ) / = '*'
2012-07-20 21:03:13 +05:30
do j = 1_pInt , 9_pInt
2012-08-03 14:55:48 +05:30
if ( temp_maskVector ( j ) ) temp_valueVector ( j ) = IO_floatValue ( line , positions , i + j )
2012-07-19 22:54:56 +05:30
enddo
2012-08-03 14:55:48 +05:30
loadCases ( currentLoadCase ) % P % maskLogical = transpose ( reshape ( temp_maskVector , [ 3 , 3 ] ) )
loadCases ( currentLoadCase ) % P % maskFloat = merge ( ones , zeroes , &
loadCases ( currentLoadCase ) % P % maskLogical )
loadCases ( currentLoadCase ) % P % values = math_plain9to33 ( temp_valueVector )
2012-07-19 22:54:56 +05:30
case ( 't' , 'time' , 'delta' ) ! increment time
2012-08-03 14:55:48 +05:30
loadCases ( currentLoadCase ) % time = IO_floatValue ( line , positions , i + 1_pInt )
2012-07-19 22:54:56 +05:30
case ( 'temp' , 'temperature' ) ! starting temperature
2012-08-03 14:55:48 +05:30
loadCases ( currentLoadCase ) % temperature = IO_floatValue ( line , positions , i + 1_pInt )
2012-07-19 22:54:56 +05:30
case ( 'n' , 'incs' , 'increments' , 'steps' ) ! number of increments
2012-08-03 14:55:48 +05:30
loadCases ( currentLoadCase ) % incs = IO_intValue ( line , positions , i + 1_pInt )
2012-07-19 22:54:56 +05:30
case ( 'logincs' , 'logincrements' , 'logsteps' ) ! number of increments (switch to log time scaling)
2012-08-03 14:55:48 +05:30
loadCases ( currentLoadCase ) % incs = IO_intValue ( line , positions , i + 1_pInt )
loadCases ( currentLoadCase ) % logscale = 1_pInt
2012-07-19 22:54:56 +05:30
case ( 'f' , 'freq' , 'frequency' , 'outputfreq' ) ! frequency of result writings
2012-08-03 14:55:48 +05:30
loadCases ( currentLoadCase ) % outputfrequency = IO_intValue ( line , positions , i + 1_pInt )
2012-07-19 22:54:56 +05:30
case ( 'r' , 'restart' , 'restartwrite' ) ! frequency of writing restart information
2012-08-03 14:55:48 +05:30
loadCases ( currentLoadCase ) % restartfrequency = max ( 0_pInt , IO_intValue ( line , positions , i + 1_pInt ) )
2012-07-19 22:54:56 +05:30
case ( 'guessreset' , 'dropguessing' )
2012-08-03 14:55:48 +05:30
loadCases ( currentLoadCase ) % followFormerTrajectory = . false . ! do not continue to predict deformation along former trajectory
case ( 'euler' ) ! rotation of currentLoadCase given in euler angles
2012-07-24 22:37:10 +05:30
l = 0_pInt ! assuming values given in radians
2012-07-20 21:03:13 +05:30
k = 1_pInt ! assuming keyword indicating degree/radians
select case ( IO_lc ( IO_stringValue ( line , positions , i + 1_pInt ) ) )
2012-07-19 22:54:56 +05:30
case ( 'deg' , 'degree' )
2012-07-24 22:37:10 +05:30
l = 1_pInt ! for conversion from degree to radian
2012-07-19 22:54:56 +05:30
case ( 'rad' , 'radian' )
case default
2012-08-03 14:55:48 +05:30
k = 0_pInt ! immediately readingk in angles, assuming radians
2012-07-19 22:54:56 +05:30
end select
2012-07-20 21:03:13 +05:30
forall ( j = 1_pInt : 3_pInt ) temp33_Real ( j , 1 ) = &
2012-07-24 22:37:10 +05:30
IO_floatValue ( line , positions , i + k + j ) * real ( l , pReal ) * inRad
2012-08-03 14:55:48 +05:30
loadCases ( currentLoadCase ) % rotation = math_EulerToR ( temp33_Real ( : , 1 ) )
case ( 'rotation' , 'rot' ) ! assign values for the rotation of currentLoadCase matrix
2012-07-19 22:54:56 +05:30
temp_valueVector = 0.0_pReal
2012-07-20 21:03:13 +05:30
forall ( j = 1_pInt : 9_pInt ) temp_valueVector ( j ) = IO_floatValue ( line , positions , i + j )
2012-08-03 14:55:48 +05:30
loadCases ( currentLoadCase ) % rotation = math_plain9to33 ( temp_valueVector )
2012-07-19 22:54:56 +05:30
end select
enddo ; enddo
101 close ( myUnit )
2012-08-03 14:55:48 +05:30
2012-07-19 22:54:56 +05:30
!--------------------------------------------------------------------------------------------------
! consistency checks and output of load case
2012-08-03 14:55:48 +05:30
loadCases ( 1 ) % followFormerTrajectory = . false . ! cannot guess along trajectory for first inc of first currentLoadCase
2012-07-19 22:54:56 +05:30
errorID = 0_pInt
2012-08-03 14:55:48 +05:30
checkLoadcases : do currentLoadCase = 1_pInt , size ( loadCases )
write ( loadcase_string , '(i6)' ) currentLoadCase
2012-07-19 22:54:56 +05:30
2012-08-03 14:55:48 +05:30
write ( 6 , '(2x,a,i6)' ) 'load case: ' , currentLoadCase
2012-07-19 22:54:56 +05:30
2012-08-03 14:55:48 +05:30
if ( . not . loadCases ( currentLoadCase ) % followFormerTrajectory ) write ( 6 , '(2x,a)' ) 'drop guessing along trajectory'
if ( loadCases ( currentLoadCase ) % deformation % myType == 'l' ) then
2012-07-19 22:54:56 +05:30
do j = 1_pInt , 3_pInt
2012-08-03 14:55:48 +05:30
if ( any ( loadCases ( currentLoadCase ) % deformation % maskLogical ( j , 1 : 3 ) . eqv . . true . ) . and . &
any ( loadCases ( currentLoadCase ) % deformation % maskLogical ( j , 1 : 3 ) . eqv . . false . ) ) errorID = 832_pInt ! each row should be either fully or not at all defined
2012-07-19 22:54:56 +05:30
enddo
2012-08-03 14:55:48 +05:30
write ( 6 , '(2x,a)' ) 'velocity gradient:'
2012-07-19 22:54:56 +05:30
else
2012-08-03 14:55:48 +05:30
write ( 6 , '(2x,a)' ) 'deformation gradient rate:'
2012-07-19 22:54:56 +05:30
endif
2012-08-03 14:55:48 +05:30
write ( 6 , '(3(3(3x,f12.7,1x)/))' , advance = 'no' ) merge ( math_transpose33 ( loadCases ( currentLoadCase ) % deformation % values ) , &
reshape ( spread ( DAMASK_NaN , 1 , 9 ) , [ 3 , 3 ] ) , transpose ( loadCases ( currentLoadCase ) % deformation % maskLogical ) )
write ( 6 , '(2x,a,/,3(3(3x,f12.7,1x)/))' , advance = 'no' ) 'stress / GPa:' , &
1e-9_pReal * merge ( math_transpose33 ( loadCases ( currentLoadCase ) % P % values ) , &
reshape ( spread ( DAMASK_NaN , 1 , 9 ) , [ 3 , 3 ] ) , transpose ( loadCases ( currentLoadCase ) % P % maskLogical ) )
if ( any ( loadCases ( currentLoadCase ) % rotation / = math_I3 ) ) &
write ( 6 , '(2x,a,/,3(3(3x,f12.7,1x)/))' , advance = 'no' ) 'rotation of loadframe:' , &
math_transpose33 ( loadCases ( currentLoadCase ) % rotation )
write ( 6 , '(2x,a,f12.6)' ) 'temperature:' , loadCases ( currentLoadCase ) % temperature
write ( 6 , '(2x,a,f12.6)' ) 'time: ' , loadCases ( currentLoadCase ) % time
write ( 6 , '(2x,a,i5)' ) 'increments: ' , loadCases ( currentLoadCase ) % incs
write ( 6 , '(2x,a,i5)' ) 'output frequency: ' , loadCases ( currentLoadCase ) % outputfrequency
write ( 6 , '(2x,a,i5)' ) 'restart frequency: ' , loadCases ( currentLoadCase ) % restartfrequency
2012-07-19 22:54:56 +05:30
2012-08-03 14:55:48 +05:30
if ( any ( loadCases ( currentLoadCase ) % P % maskLogical . eqv . loadCases ( currentLoadCase ) % deformation % maskLogical ) ) errorID = 831_pInt ! exclusive or masking only
if ( any ( loadCases ( currentLoadCase ) % P % maskLogical . and . transpose ( loadCases ( currentLoadCase ) % P % maskLogical ) . and . &
2012-07-19 22:54:56 +05:30
reshape ( [ . false . , . true . , . true . , . true . , . false . , . true . , . true . , . true . , . false . ] , [ 3 , 3 ] ) ) ) &
errorID = 838_pInt ! no rotation is allowed by stress BC
2012-08-03 14:55:48 +05:30
if ( any ( abs ( math_mul33x33 ( loadCases ( currentLoadCase ) % rotation , math_transpose33 ( loadCases ( currentLoadCase ) % rotation ) ) &
2012-07-19 22:54:56 +05:30
- math_I3 ) > reshape ( spread ( rotation_tol , 1 , 9 ) , [ 3 , 3 ] ) ) &
2012-08-03 14:55:48 +05:30
. or . abs ( math_det33 ( loadCases ( currentLoadCase ) % rotation ) ) > 1.0_pReal + rotation_tol ) &
2012-07-19 22:54:56 +05:30
errorID = 846_pInt ! given rotation matrix contains strain
2012-08-03 14:55:48 +05:30
if ( loadCases ( currentLoadCase ) % time < 0.0_pReal ) errorID = 834_pInt ! negative time increment
if ( loadCases ( currentLoadCase ) % incs < 1_pInt ) errorID = 835_pInt ! non-positive incs count
if ( loadCases ( currentLoadCase ) % outputfrequency < 1_pInt ) errorID = 836_pInt ! non-positive result frequency
2012-07-19 22:54:56 +05:30
if ( errorID > 0_pInt ) call IO_error ( error_ID = errorID , ext_msg = loadcase_string )
2012-07-24 22:37:10 +05:30
enddo checkLoadcases
2012-07-20 21:03:13 +05:30
2012-08-03 14:55:48 +05:30
select case ( myspectralsolver )
case ( DAMASK_spectral_SolverBasic_label )
call basic_init ( )
2012-08-14 22:28:23 +05:30
case ( DAMASK_spectral_SolverBasicPETSC_label )
call BasicPETSC_init ( )
2012-08-06 14:23:12 +05:30
case ( DAMASK_spectral_SolverAL_label )
call AL_init ( )
2012-08-03 14:55:48 +05:30
end select
2012-07-20 21:03:13 +05:30
!--------------------------------------------------------------------------------------------------
! 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 )
2012-08-03 14:55:48 +05:30
write ( 538 ) 'resolution' , res
write ( 538 ) 'dimension' , geomdim
2012-07-20 21:03:13 +05:30
write ( 538 ) 'materialpoint_sizeResults' , materialpoint_sizeResults
2012-08-03 14:55:48 +05:30
write ( 538 ) 'loadcases' , size ( loadCases )
write ( 538 ) 'frequencies' , loadCases % outputfrequency ! one entry per currentLoadCase
write ( 538 ) 'times' , loadCases % time ! one entry per currentLoadCase
write ( 538 ) 'logscales' , loadCases % logscale
write ( 538 ) 'increments' , loadCases % incs ! one entry per currentLoadCase
2012-07-20 21:03:13 +05:30
write ( 538 ) 'startingIncrement' , restartInc - 1_pInt ! start with writing out the previous inc
write ( 538 ) 'eoh' ! end of header
2012-08-03 14:55:48 +05:30
write ( 538 ) materialpoint_results ( 1_pInt : materialpoint_sizeResults , 1 , 1_pInt : mesh_NcpElems ) ! initial (non-deformed or read-in) results
2012-07-20 21:03:13 +05:30
if ( debugGeneral ) write ( 6 , '(a)' ) 'Header of result file written out'
endif
2012-07-23 15:42:31 +05:30
2012-08-03 14:55:48 +05:30
!--------------------------------------------------------------------------------------------------
! loopping over loadcases
loadCaseLooping : do currentLoadCase = 1_pInt , size ( loadCases )
time0 = time ! currentLoadCase start time
if ( loadCases ( currentLoadCase ) % followFormerTrajectory ) then
2012-07-19 22:54:56 +05:30
guessmode = 1.0_pReal
else
guessmode = 0.0_pReal ! change of load case, homogeneous guess for the first inc
endif
2012-08-03 14:55:48 +05:30
!--------------------------------------------------------------------------------------------------
! loop oper incs defined in input file for current currentLoadCase
incLooping : do inc = 1_pInt , loadCases ( currentLoadCase ) % incs
2012-07-19 22:54:56 +05:30
totalIncsCounter = totalIncsCounter + 1_pInt
!--------------------------------------------------------------------------------------------------
! forwarding time
timeinc_old = timeinc
2012-08-03 14:55:48 +05:30
if ( loadCases ( currentLoadCase ) % logscale == 0_pInt ) then ! linear scale
timeinc = loadCases ( currentLoadCase ) % time / loadCases ( currentLoadCase ) % incs ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used
2012-07-19 22:54:56 +05:30
else
2012-08-03 14:55:48 +05:30
if ( currentLoadCase == 1_pInt ) then ! 1st currentLoadCase of logarithmic scale
if ( inc == 1_pInt ) then ! 1st inc of 1st currentLoadCase of logarithmic scale
timeinc = loadCases ( 1 ) % time * ( 2.0_pReal ** real ( 1_pInt - loadCases ( 1 ) % incs , pReal ) ) ! assume 1st inc is equal to 2nd
else ! not-1st inc of 1st currentLoadCase of logarithmic scale
timeinc = loadCases ( 1 ) % time * ( 2.0_pReal ** real ( inc - 1_pInt - loadCases ( 1 ) % incs , pReal ) )
2012-07-19 22:54:56 +05:30
endif
2012-08-03 14:55:48 +05:30
else ! not-1st currentLoadCase of logarithmic scale
timeinc = time0 * ( ( 1.0_pReal + loadCases ( currentLoadCase ) % time / time0 ) ** ( real ( inc , pReal ) / &
real ( loadCases ( currentLoadCase ) % incs , pReal ) ) &
- ( 1.0_pReal + loadCases ( currentLoadCase ) % time / time0 ) ** ( real ( ( inc - 1_pInt ) , pReal ) / &
real ( loadCases ( currentLoadCase ) % incs , pReal ) ) )
2012-07-19 22:54:56 +05:30
endif
endif
time = time + timeinc
if ( totalIncsCounter > = restartInc ) then ! do calculations (otherwise just forwarding)
!--------------------------------------------------------------------------------------------------
! report begin of new increment
write ( 6 , '(a)' ) '##################################################################'
write ( 6 , '(A,I5.5,A,es12.5)' ) 'Increment ' , totalIncsCounter , ' Time ' , time
2012-07-25 19:31:39 +05:30
select case ( myspectralsolver )
case ( DAMASK_spectral_SolverBasic_label )
solres = basic_solution ( &
guessmode , timeinc , timeinc_old , &
2012-08-03 14:55:48 +05:30
P_BC = loadCases ( currentLoadCase ) % P , &
F_BC = loadCases ( currentLoadCase ) % deformation , &
temperature_bc = loadCases ( currentLoadCase ) % temperature , &
rotation_BC = loadCases ( currentLoadCase ) % rotation )
2012-07-25 19:31:39 +05:30
2012-08-14 22:28:23 +05:30
case ( DAMASK_spectral_SolverBasicPETSC_label )
solres = BasicPETSC_solution ( &
guessmode , timeinc , timeinc_old , &
P_BC = loadCases ( currentLoadCase ) % P , &
F_BC = loadCases ( currentLoadCase ) % deformation , &
temperature_bc = loadCases ( currentLoadCase ) % temperature , &
rotation_BC = loadCases ( currentLoadCase ) % rotation )
2012-08-06 14:23:12 +05:30
case ( DAMASK_spectral_SolverAL_label )
solres = AL_solution ( &
guessmode , timeinc , timeinc_old , &
P_BC = loadCases ( currentLoadCase ) % P , &
F_BC = loadCases ( currentLoadCase ) % deformation , &
temperature_bc = loadCases ( currentLoadCase ) % temperature , &
rotation_BC = loadCases ( currentLoadCase ) % rotation )
2012-07-25 19:31:39 +05:30
end select
2012-07-20 21:03:13 +05:30
2012-07-19 22:54:56 +05:30
write ( 6 , '(a)' ) ''
write ( 6 , '(a)' ) '=================================================================='
2012-07-20 21:03:13 +05:30
if ( solres % converged ) then
2012-07-19 22:54:56 +05:30
convergedCounter = convergedCounter + 1_pInt
write ( 6 , '(A,I5.5,A)' ) 'increment ' , totalIncsCounter , ' converged'
2012-07-20 21:03:13 +05:30
else
write ( 6 , '(A,I5.5,A)' ) 'increment ' , totalIncsCounter , ' NOT converged'
notConvergedCounter = notConvergedCounter + 1_pInt
2012-07-19 22:54:56 +05:30
endif
2012-08-03 14:55:48 +05:30
if ( mod ( inc , loadCases ( currentLoadCase ) % outputFrequency ) == 0_pInt ) then ! at output frequency
2012-07-19 22:54:56 +05:30
write ( 6 , '(a)' ) ''
write ( 6 , '(a)' ) '... writing results to file ......................................'
2012-08-03 14:55:48 +05:30
write ( 538 ) materialpoint_results ! write result to file
2012-07-19 22:54:56 +05:30
endif
endif ! end calculation/forwarding
2012-08-03 14:55:48 +05:30
guessmode = 1.0_pReal ! keep guessing along former trajectory during same currentLoadCase
2012-07-20 21:03:13 +05:30
2012-07-24 22:37:10 +05:30
enddo incLooping
enddo loadCaseLooping
2012-08-03 14:55:48 +05:30
2012-07-25 19:31:39 +05:30
select case ( myspectralsolver )
case ( DAMASK_spectral_SolverBasic_label )
call basic_destroy ( )
2012-08-14 22:28:23 +05:30
case ( DAMASK_spectral_SolverBasicPETSC_label )
call BasicPETSC_destroy ( )
2012-08-06 14:23:12 +05:30
case ( DAMASK_spectral_SolverAL_label )
call AL_destroy ( )
2012-08-03 14:55:48 +05:30
end select
2012-07-20 21:03:13 +05:30
write ( 6 , '(a)' ) ''
write ( 6 , '(a)' ) '##################################################################'
write ( 6 , '(i6.6,a,i6.6,a,f5.1,a)' ) convergedCounter , ' out of ' , &
notConvergedCounter + convergedCounter , ' (' , &
real ( convergedCounter , pReal ) / &
real ( notConvergedCounter + convergedCounter , pReal ) * 10 0.0_pReal , &
' %) increments converged!'
2012-07-19 22:54:56 +05:30
close ( 538 )
if ( notConvergedCounter > 0_pInt ) call quit ( 3_pInt )
call quit ( 0_pInt )
2012-07-20 21:03:13 +05:30
2012-07-24 22:37:10 +05:30
end program DAMASK_spectral_Driver
2012-07-20 21:03:13 +05:30
2012-07-23 15:42:31 +05:30
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
2012-07-20 21:03:13 +05:30
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)
2012-07-23 15:42:31 +05:30
end subroutine
