2012-07-30 19:36:22 +05:30
!--------------------------------------------------------------------------------------------------
2012-10-24 17:01:40 +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-08-31 01:56:28 +05:30
IO_write_jobBinaryFile , &
IO_intOut
2012-10-24 17:01:40 +05:30
use math ! need to include the whole module for FFTW
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-10-02 20:56:56 +05:30
maxCutBack , &
2012-07-25 19:31:39 +05:30
rotation_tol , &
2012-08-03 14:55:48 +05:30
mySpectralSolver
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 : &
2012-10-02 20:56:56 +05:30
tBoundaryCondition , &
tSolutionState , &
debugGeneral , &
cutBack
2012-07-24 22:37:10 +05:30
use DAMASK_spectral_SolverBasic
2012-08-29 00:49:47 +05:30
#ifdef PETSc
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-08-29 00:49:47 +05:30
#endif
2012-07-19 22:54:56 +05:30
implicit none
2012-10-02 20:56:56 +05:30
type tLoadCase
2012-10-19 14:14:21 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: rotation = math_I3 !< rotation of BC
type ( tBoundaryCondition ) :: P , & !< stress BC
deformation !< deformation BC (Fdot or L)
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
logical :: followFormerTrajectory = . true . !< follow trajectory of former loadcase
2012-10-02 20:56:56 +05:30
end type tLoadCase
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-11-09 03:03:58 +05:30
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
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
2012-10-02 20:56:56 +05:30
1_pInt ! dropguessing
2012-07-20 21:03:13 +05:30
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-10-24 17:01:40 +05:30
N_Fdot = 0_pInt , & !
2012-10-02 20:56:56 +05:30
myUnit = 234_pInt
2012-07-20 21:03:13 +05:30
character ( len = 1024 ) :: &
line
2012-07-19 22:54:56 +05:30
!--------------------------------------------------------------------------------------------------
! loop variables, convergence etc.
2012-10-19 14:14:21 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , parameter :: ones = 1.0_pReal , zeros = 0.0_pReal
2012-08-03 14:55:48 +05:30
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-11-09 01:02:00 +05:30
logical :: guess
2012-10-02 20:56:56 +05:30
integer ( pInt ) :: i , j , k , l , errorID , cutBackLevel = 0_pInt , stepFraction = 0_pInt
2012-07-24 22:37:10 +05:30
integer ( pInt ) :: currentLoadcase = 0_pInt , inc , &
2012-07-20 21:03:13 +05:30
totalIncsCounter = 0_pInt , &
2012-10-19 14:14:21 +05:30
notConvergedCounter = 0_pInt , convergedCounter = 0_pInt , &
resUnit = 0_pInt , statUnit = 0_pInt
2012-07-23 15:42:31 +05:30
character ( len = 6 ) :: loadcase_string
2012-10-02 20:56:56 +05:30
character ( len = 1024 ) :: incInfo
type ( tLoadCase ) , allocatable , dimension ( : ) :: loadCases
type ( tSolutionState ) 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-10-19 14:14:21 +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"
2012-10-24 17:01:40 +05:30
2012-07-19 22:54:56 +05:30
!--------------------------------------------------------------------------------------------------
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-10-24 17:01:40 +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-10-24 17:01:40 +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-11-09 02:05:31 +05:30
temp_valueVector = 0.0_pReal
2012-10-24 17:01:40 +05:30
if ( IO_lc ( IO_stringValue ( line , positions , i ) ) == 'l' . or . & ! in case of given L, set flag to true
2012-08-03 14:55:48 +05:30
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 ] ) )
2012-10-19 14:14:21 +05:30
loadCases ( currentLoadCase ) % deformation % maskFloat = merge ( ones , zeros , &
2012-08-03 14:55:48 +05:30
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 ] ) )
2012-10-19 14:14:21 +05:30
loadCases ( currentLoadCase ) % P % maskFloat = merge ( ones , zeros , &
2012-08-03 14:55:48 +05:30
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-10-24 17:01:40 +05:30
loadCases ( currentLoadCase ) % followFormerTrajectory = . false . ! do not continue to predict deformation along former trajectory
2012-08-03 14:55:48 +05:30
case ( 'euler' ) ! rotation of currentLoadCase given in euler angles
2012-11-09 03:03:58 +05:30
temp_valueVector = 0.0_pReal
2012-09-06 19:35:28 +05:30
l = 0_pInt ! assuming values given in degrees
k = 0_pInt ! assuming keyword indicating degree/radians
2012-07-20 21:03:13 +05:30
select case ( IO_lc ( IO_stringValue ( line , positions , i + 1_pInt ) ) )
2012-09-06 19:35:28 +05:30
case ( 'rad' , 'radian' )
l = 1_pInt
k = 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 default
end select
2012-09-06 19:35:28 +05:30
forall ( j = 1_pInt : 3_pInt ) temp_valueVector ( j ) = IO_floatValue ( line , positions , i + k + j )
if ( k == 1_pInt ) temp_valueVector ( 1 : 3 ) = temp_valueVector ( 1 : 3 ) * inRad ! convert to rad
loadCases ( currentLoadCase ) % rotation = math_EulerToR ( temp_valueVector ( 1 : 3 ) ) ! convert rad Eulers to rotation matrix
2012-08-03 14:55:48 +05:30
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-10-24 17:01:40 +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-10-02 20:56:56 +05:30
write ( 6 , '(1x,a,i6)' ) 'load case: ' , currentLoadCase
2012-07-19 22:54:56 +05:30
2012-10-24 17:01:40 +05:30
if ( . not . loadCases ( currentLoadCase ) % followFormerTrajectory ) &
write ( 6 , '(2x,a)' ) 'drop guessing along trajectory'
2012-08-03 14:55:48 +05:30
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 . &
2012-10-24 17:01:40 +05:30
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-10-24 17:01:40 +05:30
write ( 6 , '(3(3(3x,f12.7,1x)/))' , advance = 'no' ) &
merge ( math_transpose33 ( loadCases ( currentLoadCase ) % deformation % values ) , &
reshape ( spread ( huge ( 1.0_pReal ) , 1 , 9 ) , [ 3 , 3 ] ) , &
transpose ( loadCases ( currentLoadCase ) % deformation % maskLogical ) )
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 . &
reshape ( [ . false . , . true . , . true . , . true . , . false . , . true . , . true . , . true . , . false . ] , [ 3 , 3 ] ) ) ) &
errorID = 838_pInt ! no rotation is allowed by stress BC
2012-10-02 20:56:56 +05:30
write ( 6 , '(2x,a,/,3(3(3x,f12.7,1x)/))' , advance = 'no' ) 'stress / GPa:' , &
2012-10-24 17:01:40 +05:30
1e-9_pReal * merge ( math_transpose33 ( loadCases ( currentLoadCase ) % P % values ) , &
reshape ( spread ( huge ( 1.0_pReal ) , 1 , 9 ) , [ 3 , 3 ] ) , &
transpose ( loadCases ( currentLoadCase ) % P % maskLogical ) )
if ( any ( abs ( math_mul33x33 ( loadCases ( currentLoadCase ) % rotation , &
math_transpose33 ( loadCases ( currentLoadCase ) % rotation ) ) - math_I3 ) > &
reshape ( spread ( rotation_tol , 1 , 9 ) , [ 3 , 3 ] ) ) &
. or . abs ( math_det33 ( loadCases ( currentLoadCase ) % rotation ) ) > &
1.0_pReal + rotation_tol ) errorID = 846_pInt ! given rotation matrix contains strain
2012-08-03 14:55:48 +05:30
if ( any ( loadCases ( currentLoadCase ) % rotation / = math_I3 ) ) &
2012-10-02 20:56:56 +05:30
write ( 6 , '(2x,a,/,3(3(3x,f12.7,1x)/))' , advance = 'no' ) 'rotation of loadframe:' , &
2012-10-24 17:01:40 +05:30
math_transpose33 ( loadCases ( currentLoadCase ) % rotation )
2012-08-03 14:55:48 +05:30
write ( 6 , '(2x,a,f12.6)' ) 'temperature:' , loadCases ( currentLoadCase ) % temperature
2012-10-24 17:01:40 +05:30
if ( loadCases ( currentLoadCase ) % time < 0.0_pReal ) errorID = 834_pInt ! negative time increment
2012-08-03 14:55:48 +05:30
write ( 6 , '(2x,a,f12.6)' ) 'time: ' , loadCases ( currentLoadCase ) % time
2012-10-24 17:01:40 +05:30
if ( loadCases ( currentLoadCase ) % incs < 1_pInt ) errorID = 835_pInt ! non-positive incs count
2012-08-03 14:55:48 +05:30
write ( 6 , '(2x,a,i5)' ) 'increments: ' , loadCases ( currentLoadCase ) % incs
2012-10-24 17:01:40 +05:30
if ( loadCases ( currentLoadCase ) % outputfrequency < 1_pInt ) errorID = 836_pInt ! non-positive result frequency
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
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-29 00:49:47 +05:30
#ifdef PETSc
2012-10-02 20:56:56 +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-29 00:49:47 +05:30
#endif
case default
call IO_error ( error_ID = 891 , ext_msg = trim ( myspectralsolver ) )
2012-08-03 14:55:48 +05:30
end select
2012-08-29 00:49:47 +05:30
2012-07-20 21:03:13 +05:30
!--------------------------------------------------------------------------------------------------
! write header of output file
if ( appendToOutFile ) then
2012-10-19 14:14:21 +05:30
open ( newunit = resUnit , file = trim ( getSolverWorkingDirectoryName ( ) ) / / trim ( getSolverJobName ( ) ) / / &
2012-10-24 17:01:40 +05:30
'.spectralOut' , form = 'UNFORMATTED' , position = 'APPEND' , status = 'OLD' )
2012-10-19 14:14:21 +05:30
open ( newunit = statUnit , file = trim ( getSolverWorkingDirectoryName ( ) ) / / trim ( getSolverJobName ( ) ) / / &
2012-10-24 17:01:40 +05:30
'.sta' , form = 'FORMATTED' , position = 'APPEND' , status = 'OLD' )
2012-07-20 21:03:13 +05:30
else
2012-10-19 14:14:21 +05:30
open ( newunit = resUnit , file = trim ( getSolverWorkingDirectoryName ( ) ) / / trim ( getSolverJobName ( ) ) / / &
2012-10-24 17:01:40 +05:30
'.spectralOut' , form = 'UNFORMATTED' , status = 'REPLACE' )
2012-10-19 14:14:21 +05:30
write ( resUnit ) 'load' , trim ( loadCaseFile )
write ( resUnit ) 'workingdir' , trim ( getSolverWorkingDirectoryName ( ) )
write ( resUnit ) 'geometry' , trim ( geometryFile )
write ( resUnit ) 'resolution' , res
write ( resUnit ) 'dimension' , geomdim
write ( resUnit ) 'materialpoint_sizeResults' , materialpoint_sizeResults
write ( resUnit ) 'loadcases' , size ( loadCases )
2012-10-24 17:01:40 +05:30
write ( resUnit ) 'frequencies' , loadCases % outputfrequency ! one entry per currentLoadCase
write ( resUnit ) 'times' , loadCases % time ! one entry per currentLoadCase
2012-10-19 14:14:21 +05:30
write ( resUnit ) 'logscales' , loadCases % logscale
2012-10-24 17:01:40 +05:30
write ( resUnit ) 'increments' , loadCases % incs ! one entry per currentLoadCase
write ( resUnit ) 'startingIncrement' , restartInc - 1_pInt ! start with writing out the previous inc
write ( resUnit ) 'eoh' ! end of header
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 ( ) ) / / &
'.sta' , form = 'FORMATTED' , status = 'REPLACE' )
write ( statUnit , '(a)' ) 'Increment Time CutbackLevel Converged IterationsNeeded'
2012-07-20 21:03:13 +05:30
if ( debugGeneral ) write ( 6 , '(a)' ) 'Header of result file written out'
endif
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-11-09 01:02:00 +05:30
guess = . true .
2012-07-19 22:54:56 +05:30
else
2012-11-09 01:02:00 +05:30
guess = . false . ! change of load case, homogeneous guess for the first inc
2012-07-19 22:54:56 +05:30
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-10-02 20:56:56 +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-10-02 20:56:56 +05:30
if ( currentLoadCase == 1_pInt ) then ! 1st currentLoadCase of logarithmic scale
2012-08-03 14:55:48 +05:30
if ( inc == 1_pInt ) then ! 1st inc of 1st currentLoadCase of logarithmic scale
2012-10-02 20:56:56 +05:30
timeinc = loadCases ( 1 ) % time * ( 2.0_pReal ** real ( 1_pInt - loadCases ( 1 ) % incs , pReal ) ) ! assume 1st inc is equal to 2nd
2012-08-03 14:55:48 +05:30
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
2012-10-24 17:01:40 +05:30
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
2012-10-02 20:56:56 +05:30
timeinc = timeinc / 2.0_pReal ** real ( cutBackLevel , pReal )
2012-07-19 22:54:56 +05:30
if ( totalIncsCounter > = restartInc ) then ! do calculations (otherwise just forwarding)
2012-10-02 20:56:56 +05:30
stepFraction = 0_pInt
2012-07-19 22:54:56 +05:30
!--------------------------------------------------------------------------------------------------
! report begin of new increment
2012-10-02 20:56:56 +05:30
subIncLooping : do while ( stepFraction / 2_pInt ** cutBackLevel < 1_pInt )
time = time + timeinc
stepFraction = stepFraction + 1_pInt
write ( 6 , '(1/,a)' ) '###########################################################################'
write ( 6 , '(a,es12.5' / / &
',a,' / / IO_intOut ( inc ) / / ',a,' / / IO_intOut ( loadCases ( currentLoadCase ) % incs ) / / &
',a,' / / IO_intOut ( stepFraction ) / / ',a,' / / IO_intOut ( 2_pInt ** cutBackLevel ) / / &
',a,' / / IO_intOut ( currentLoadCase ) / / ',a,' / / IO_intOut ( size ( loadCases ) ) / / ')' ) &
2012-10-24 17:01:40 +05:30
'Time' , time , &
's: Increment ' , inc , '/' , loadCases ( currentLoadCase ) % incs , &
'-' , stepFraction , '/' , 2_pInt ** cutBackLevel , &
' of load case ' , currentLoadCase , '/' , size ( loadCases )
2012-10-02 20:56:56 +05:30
write ( incInfo , '(a,' / / IO_intOut ( totalIncsCounter ) / / ',a,' / / IO_intOut ( sum ( loadCases ( : ) % incs ) ) / / &
2012-10-24 17:01:40 +05:30
',a,' / / IO_intOut ( stepFraction ) / / ',a,' / / IO_intOut ( 2_pInt ** cutBackLevel ) / / ')' ) &
2012-10-02 20:56:56 +05:30
'Inc. ' , totalIncsCounter , '/' , sum ( loadCases ( : ) % incs ) , &
'-' , stepFraction , '/' , 2_pInt ** cutBackLevel
select case ( myspectralsolver )
case ( DAMASK_spectral_SolverBasic_label )
solres = basic_solution ( &
2012-11-09 01:02:00 +05:30
incInfo , guess , timeinc , timeinc_old , &
2012-10-02 20:56:56 +05:30
P_BC = loadCases ( currentLoadCase ) % P , &
F_BC = loadCases ( currentLoadCase ) % deformation , &
temperature_bc = loadCases ( currentLoadCase ) % temperature , &
rotation_BC = loadCases ( currentLoadCase ) % rotation )
2012-08-29 00:49:47 +05:30
#ifdef PETSc
2012-10-02 20:56:56 +05:30
case ( DAMASK_spectral_SolverBasicPETSC_label )
solres = BasicPETSC_solution ( &
2012-11-09 01:02:00 +05:30
incInfo , guess , timeinc , timeinc_old , &
2012-10-02 20:56:56 +05:30
P_BC = loadCases ( currentLoadCase ) % P , &
F_BC = loadCases ( currentLoadCase ) % deformation , &
temperature_bc = loadCases ( currentLoadCase ) % temperature , &
rotation_BC = loadCases ( currentLoadCase ) % rotation )
case ( DAMASK_spectral_SolverAL_label )
solres = AL_solution ( &
2012-11-09 01:02:00 +05:30
incInfo , guess , timeinc , timeinc_old , &
2012-10-02 20:56:56 +05:30
P_BC = loadCases ( currentLoadCase ) % P , &
F_BC = loadCases ( currentLoadCase ) % deformation , &
temperature_bc = loadCases ( currentLoadCase ) % temperature , &
rotation_BC = loadCases ( currentLoadCase ) % rotation )
2012-08-29 00:49:47 +05:30
#endif
2012-10-02 20:56:56 +05:30
end select
cutBack = . False .
if ( solres % termIll . or . . not . solres % converged ) then ! no solution found
if ( cutBackLevel < maxCutBack ) then ! do cut back
write ( 6 , '(/,a)' ) 'cut back detected'
cutBack = . True .
stepFraction = ( stepFraction - 1_pInt ) * 2_pInt ** cutBackLevel
cutBackLevel = cutBackLevel + 1_pInt
time = time - timeinc
timeinc_old = timeinc
timeinc = timeinc / 2.0_pReal
elseif ( solres % termIll ) then ! material point model cannot find a solution
call IO_error ( 850_pInt )
else
2012-11-09 01:02:00 +05:30
guess = . true . ! start guessing after first accepted (not converged) (sub)inc
2012-10-02 20:56:56 +05:30
endif
else
2012-11-09 01:02:00 +05:30
guess = . true . ! start guessing after first converged (sub)inc
2012-10-02 20:56:56 +05:30
endif
2012-11-09 01:02:00 +05:30
if ( guess ) &
2012-10-19 14:14:21 +05:30
write ( statUnit , * ) inc , time , cutBackLevel , solres % converged , solres % iterationsNeeded
2012-10-02 20:56:56 +05:30
enddo subIncLooping
cutBackLevel = max ( 0_pInt , cutBackLevel - 1_pInt ) ! try half subincs next inc
2012-07-20 21:03:13 +05:30
if ( solres % converged ) then
2012-07-19 22:54:56 +05:30
convergedCounter = convergedCounter + 1_pInt
2012-08-31 01:56:28 +05:30
write ( 6 , '(A,' / / IO_intOut ( totalIncsCounter ) / / ',A)' ) &
2012-10-02 20:56:56 +05:30
'increment ' , totalIncsCounter , ' converged'
2012-07-20 21:03:13 +05:30
else
2012-08-31 01:56:28 +05:30
write ( 6 , '(A,' / / IO_intOut ( totalIncsCounter ) / / ',A)' ) &
2012-10-02 20:56:56 +05:30
'increment ' , totalIncsCounter , ' NOT converged'
2012-07-20 21:03:13 +05:30
notConvergedCounter = notConvergedCounter + 1_pInt
2012-07-19 22:54:56 +05:30
endif
2012-10-02 20:56:56 +05:30
if ( mod ( inc , loadCases ( currentLoadCase ) % outputFrequency ) == 0_pInt ) then ! at output frequency
write ( 6 , '(1/,a)' ) '... writing results to file ......................................'
2012-10-19 14:14:21 +05:30
write ( resUnit ) materialpoint_results ! write result to file
2012-07-19 22:54:56 +05:30
endif
2012-10-24 17:01:40 +05:30
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?)
restartWrite = . true .
endif
2012-10-02 20:56:56 +05:30
else !just time forwarding
2012-10-24 17:01:40 +05:30
time = time + timeinc
2012-11-09 01:02:00 +05:30
guess = . true .
2012-10-02 20:56:56 +05:30
endif ! end calculation/forwarding
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-29 00:49:47 +05:30
#ifdef PETSc
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-29 00:49:47 +05:30
#endif
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-10-19 14:14:21 +05:30
close ( resUnit )
2012-07-19 22:54:56 +05:30
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-08-20 18:52:34 +05:30
#include "spectral_quit.f90"
