DAMASK_grid.f90

Go to the documentation of this file.

# 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/grid/DAMASK_grid.f90"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/grid/DAMASK_grid.f90"
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
program damask_spectral
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h" 1
!
!
!  Part of the base include file for Fortran use of 1.
!  Note: This file should contain only define statements and
!  not the declaration of variables.
 
! No spaces for #defines as some compilers (PGI) also adds
! those additional spaces during preprocessing - bad for fixed format
!
 
 
 
# 1 "/opt/petsc-3.10.3/Intel-18.4-IntelMPI-2018/include/petscconf.h" 1
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
# 13 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h" 2
 
 
 
 
# 1 "/opt/petsc-3.10.3/include/petscversion.h" 1
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
# 17 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscviewer.h" 1
!
!  Include file for Fortran use of the PetscViewer package in 1
!
 
 
 
 
 
 
 
 
 
 
 
 
 
# 31 "/opt/petsc-3.10.3/include/petsc/finclude/petscviewer.h"
 
# 18 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscerror.h" 1
 
!
!  Include file for Fortran error codes
!    These are also in include/petscerror.h
!
 
 
 
# 24 "/opt/petsc-3.10.3/include/petsc/finclude/petscerror.h"
 
# 38 "/opt/petsc-3.10.3/include/petsc/finclude/petscerror.h"
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
# 19 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petsclog.h" 1
!
!  No includes needed for logging
# 20 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscbag.h" 1
!
!
!  Include file for Fortran use of the Bag package in 1
!
 
 
 
 
 
!
!  End of Fortran include file for the IS package in 1
 
# 21 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h" 2
 
!
! The real*8,complex*16 notatiton is used so that the
! 1 double/complex variables are not affected by
! compiler options like -r4,-r8, sometimes invoked
! by the user. NAG compiler does not like integer*4,real*8
 
# 41 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h"
 
 
 
 
 
 
 
 
 
# 63 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h"
 
 
 
 
 
 
!
 
 
 
 
 
!
# 85 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h"
!
 
 
 
 
 
 
!
 
 
 
 
!
 
!
 
 
!
# 128 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h"
 
# 150 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h"
!
!     Macro for templating between real and complex
!
# 174 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h"
 
 
 
 
 
 
 
 
 
!
!    Allows the matrix Fortran Kernels to work with single precision
!    matrix data structures
!
 
!
!     PetscLogDouble variables are used to contain double precision numbers
!     that are not used in the numerical computations, but rather in logging,
!     timing etc.
!
 
 
!
!     Macros for error checking
!
 
 
 
 
 
 
 
 
# 215 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h"
 
 
 
 
 
 
 
 
 
 
 
 
# 11 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/grid/DAMASK_grid.f90" 2
 use petscsys
 use prec
 use damask_interface
 use io
 use config
 use debug
 use math
 use cpfem2
 use material
 use spectral_utilities
 use grid_mech_spectral_basic
 use grid_mech_spectral_polarisation
 use grid_mech_fem
 use grid_damage_spectral
 use grid_thermal_spectral
 use results
 
 implicit none
 
!--------------------------------------------------------------------------------------------------
! variables related to information from load case and geom file
 real(preal), dimension(9) :: temp_valuevector = 0.0_preal                                          
 logical,     dimension(9) :: temp_maskvector  = .false.                                            
 integer, allocatable, dimension(:) :: chunkpos
 integer :: &
   n_t   = 0, &                                                                                     !< # of time indicators found in load case file
   n_n   = 0, &                                                                                     
   n_def = 0                                                                                        
 character(len=pStringLen) :: &
   line
 
!--------------------------------------------------------------------------------------------------
! loop variables, convergence etc.
 real(preal), dimension(3,3), parameter :: &
   ones  = 1.0_preal, &
   zeros = 0.0_preal
 integer, parameter :: &
   substepfactor = 2
 real(preal) :: &
   time = 0.0_preal, &                                                                              
   time0 = 0.0_preal, &                                                                             
   timeinc = 1.0_preal, &                                                                           
   timeincold = 0.0_preal, &                                                                        
   remainingloadcasetime = 0.0_preal                                                                
 logical :: &
   guess, &                                                                                         !< guess along former trajectory
   stagiterate, &
   cutback = .false.
 integer :: &
   i, j, k, l, field, &
   errorid = 0, &
   cutbacklevel = 0, &                                                                              
   stepfraction = 0                                                                                 
 integer :: &
   currentloadcase = 0, &                                                                           !< current load case
   inc, &                                                                                           
   totalincscounter = 0, &                                                                          
   statunit = 0, &                                                                                  
   stagiter, &
   nactivefields = 0
 character(len=pStringLen), dimension(:), allocatable :: filecontent
 character(len=pStringLen) :: &
   incinfo, &
   loadcase_string
 type(tloadcase), allocatable, dimension(:) :: loadcases
 type(tloadcase) :: newloadcase
 type(tsolutionstate), allocatable, dimension(:) :: solres
 procedure(grid_mech_spectral_basic_init), pointer :: &
   mech_init
 procedure(grid_mech_spectral_basic_forward), pointer :: &
   mech_forward
 procedure(grid_mech_spectral_basic_solution), pointer :: &
   mech_solution
 procedure(grid_mech_spectral_basic_updatecoords), pointer :: &
   mech_updatecoords
 procedure(grid_mech_spectral_basic_restartwrite), pointer :: &
   mech_restartwrite
 
 external :: &
   quit
 
!--------------------------------------------------------------------------------------------------
! init DAMASK (all modules)
 call cpfem_initall
 write(6,'(/,a)')   ' <<<+-  DAMASK_spectral init  -+>>>'; flush(6)
 
 write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, 2019'
 write(6,'(a)')   ' https://doi.org/10.1007/978-981-10-6855-3_80'
 
!--------------------------------------------------------------------------------------------------
! initialize field solver information
 nactivefields = 1
 if (any(thermal_type  == thermal_conduction_id  )) nactivefields = nactivefields + 1
 if (any(damage_type   == damage_nonlocal_id     )) nactivefields = nactivefields + 1
 allocate(solres(nactivefields))
 allocate(newloadcase%ID(nactivefields))
 
!--------------------------------------------------------------------------------------------------
! assign mechanics solver depending on selected type
 select case (trim(config_numerics%getString('spectral_solver',defaultval='basic')))
   case ('basic')
     mech_init         => grid_mech_spectral_basic_init
     mech_forward      => grid_mech_spectral_basic_forward
     mech_solution     => grid_mech_spectral_basic_solution
     mech_updatecoords => grid_mech_spectral_basic_updatecoords
     mech_restartwrite => grid_mech_spectral_basic_restartwrite
 
   case ('polarisation')
     if(iand(debug_level(debug_spectral),debug_levelbasic)/= 0) &
       call io_warning(42, ext_msg='debug Divergence')
     mech_init         => grid_mech_spectral_polarisation_init
     mech_forward      => grid_mech_spectral_polarisation_forward
     mech_solution     => grid_mech_spectral_polarisation_solution
     mech_updatecoords => grid_mech_spectral_polarisation_updatecoords
     mech_restartwrite => grid_mech_spectral_polarisation_restartwrite
 
   case ('fem')
     if(iand(debug_level(debug_spectral),debug_levelbasic)/= 0) &
       call io_warning(42, ext_msg='debug Divergence')
     mech_init         => grid_mech_fem_init
     mech_forward      => grid_mech_fem_forward
     mech_solution     => grid_mech_fem_solution
     mech_updatecoords => grid_mech_fem_updatecoords
     mech_restartwrite => grid_mech_fem_restartwrite
 
   case default
     call io_error(error_id = 891, ext_msg = config_numerics%getString('spectral_solver'))
 
 end select
 
!--------------------------------------------------------------------------------------------------
! reading information from load case file and to sanity checks 
 filecontent = io_read_ascii(trim(loadcasefile))
 
 allocate (loadcases(0))                                                                            ! array of load cases
 do currentloadcase = 1, size(filecontent)
   line = filecontent(currentloadcase)
   if (io_isblank(line)) cycle
   chunkpos = io_stringpos(line)
 
   do i = 1, chunkpos(1)                                                                            ! reading compulsory parameters for loadcase
     select case (io_lc(io_stringvalue(line,chunkpos,i)))
       case('l','fdot','dotf','f')
         n_def = n_def + 1
       case('t','time','delta')
         n_t = n_t + 1
       case('n','incs','increments','logincs','logincrements')
         n_n = n_n + 1
     end select
   enddo
   if ((n_def /= n_n) .or. (n_n /= n_t) .or. n_n < 1) &                                             ! sanity check
     call io_error(error_id=837,el=currentloadcase,ext_msg = trim(loadcasefile))                    ! error message for incomplete loadcase
 
   newloadcase%stress%myType='stress'
   field = 1
   newloadcase%ID(field) = field_mech_id                                                            ! mechanical active by default
   thermalactive: if (any(thermal_type  == thermal_conduction_id)) then
     field = field + 1
     newloadcase%ID(field) = field_thermal_id
   endif thermalactive
   damageactive: if (any(damage_type   == damage_nonlocal_id)) then
     field = field + 1
     newloadcase%ID(field) = field_damage_id
   endif damageactive
 
   call newloadcase%rot%fromEulers(real([0.0,0.0,0.0],preal))
   readin: do i = 1, chunkpos(1)
     select case (io_lc(io_stringvalue(line,chunkpos,i)))
       case('fdot','dotf','l','f')                                                                  ! assign values for the deformation BC matrix
         temp_valuevector = 0.0_preal
         if (io_lc(io_stringvalue(line,chunkpos,i)) == 'fdot'.or. &                                 ! in case of Fdot, set type to fdot
             io_lc(io_stringvalue(line,chunkpos,i)) == 'dotf') then
           newloadcase%deformation%myType = 'fdot'
         else if (io_lc(io_stringvalue(line,chunkpos,i)) == 'f') then
           newloadcase%deformation%myType = 'f'
         else
           newloadcase%deformation%myType = 'l'
         endif
         do j = 1, 9
           temp_maskvector(j) = io_stringvalue(line,chunkpos,i+j) /= '*'                            ! true if not a *
           if (temp_maskvector(j)) temp_valuevector(j) = io_floatvalue(line,chunkpos,i+j)           ! read value where applicable
         enddo
         newloadcase%deformation%maskLogical = transpose(reshape(temp_maskvector,[ 3,3]))           ! logical mask in 3x3 notation
         newloadcase%deformation%maskFloat   = merge(ones,zeros,newloadcase%deformation%maskLogical)! float (1.0/0.0) mask in 3x3 notation
         newloadcase%deformation%values      = math_9to33(temp_valuevector)                         ! values in 3x3 notation
       case('p','stress', 's')
         temp_valuevector = 0.0_preal
         do j = 1, 9
           temp_maskvector(j) = io_stringvalue(line,chunkpos,i+j) /= '*'                            ! true if not an asterisk
           if (temp_maskvector(j)) temp_valuevector(j) = io_floatvalue(line,chunkpos,i+j)           ! read value where applicable
         enddo
         newloadcase%stress%maskLogical = transpose(reshape(temp_maskvector,[ 3,3]))
         newloadcase%stress%maskFloat   = merge(ones,zeros,newloadcase%stress%maskLogical)
         newloadcase%stress%values      = math_9to33(temp_valuevector)
       case('t','time','delta')                                                                     ! increment time
         newloadcase%time = io_floatvalue(line,chunkpos,i+1)
       case('n','incs','increments')                                                                ! number of increments
         newloadcase%incs = io_intvalue(line,chunkpos,i+1)
       case('logincs','logincrements')                                                              ! number of increments (switch to log time scaling)
         newloadcase%incs = io_intvalue(line,chunkpos,i+1)
         newloadcase%logscale = 1
       case('freq','frequency','outputfreq')                                                        ! frequency of result writings
         newloadcase%outputfrequency = io_intvalue(line,chunkpos,i+1)
       case('r','restart','restartwrite')                                                           ! frequency of writing restart information
         newloadcase%restartfrequency = io_intvalue(line,chunkpos,i+1)
       case('guessreset','dropguessing')
         newloadcase%followFormerTrajectory = .false.                                               ! do not continue to predict deformation along former trajectory
       case('euler')                                                                                ! rotation of load case given in euler angles
         temp_valuevector = 0.0_preal
         l = 1                                                                                      ! assuming values given in degrees
         k = 1                                                                                      ! assuming keyword indicating degree/radians present
         select case (io_lc(io_stringvalue(line,chunkpos,i+1)))
           case('deg','degree')
           case('rad','radian')                                                                     ! don't convert from degree to radian
             l = 0
           case default
             k = 0
         end select
         do j = 1, 3
           temp_valuevector(j) = io_floatvalue(line,chunkpos,i+k+j)
         enddo
         call newloadcase%rot%fromEulers(temp_valuevector(1:3),degrees=(l==1))
       case('rotation','rot')                                                                       ! assign values for the rotation  matrix
         temp_valuevector = 0.0_preal
         do j = 1, 9
           temp_valuevector(j) = io_floatvalue(line,chunkpos,i+j)
         enddo
         call newloadcase%rot%fromMatrix(math_9to33(temp_valuevector))
     end select
   enddo readin
 
   newloadcase%followFormerTrajectory = merge(.true.,.false.,currentloadcase > 1)                   ! by default, guess from previous load case
 
   reportandcheck: if (worldrank == 0) then
     write (loadcase_string, '(i0)' ) currentloadcase
     write(6,'(/,1x,a,i0)') 'load case: ', currentloadcase
     if (.not. newloadcase%followFormerTrajectory) write(6,'(2x,a)') 'drop guessing along trajectory'
     if (newloadcase%deformation%myType == 'l') then
       do j = 1, 3
         if (any(newloadcase%deformation%maskLogical(j,1:3) .eqv. .true.) .and. &
             any(newloadcase%deformation%maskLogical(j,1:3) .eqv. .false.)) errorid = 832           ! each row should be either fully or not at all defined
       enddo
       write(6,'(2x,a)') 'velocity gradient:'
     else if (newloadcase%deformation%myType == 'f') then
       write(6,'(2x,a)') 'deformation gradient at end of load case:'
     else
       write(6,'(2x,a)') 'deformation gradient rate:'
     endif
     do i = 1, 3; do j = 1, 3
       if(newloadcase%deformation%maskLogical(i,j)) then
         write(6,'(2x,f12.7)',advance='no') newloadcase%deformation%values(i,j)
       else
         write(6,'(2x,12a)',advance='no') '     *      '
         endif
       enddo; write(6,'(/)',advance='no')
     enddo
     if (any(newloadcase%stress%maskLogical .eqv. &
             newloadcase%deformation%maskLogical)) errorid = 831                                    ! exclusive or masking only
     if (any(newloadcase%stress%maskLogical .and. transpose(newloadcase%stress%maskLogical) &
         .and. (math_i3<1))) errorid = 838                                                          ! no rotation is allowed by stress BC
     write(6,'(2x,a)') 'stress / GPa:'
     do i = 1, 3; do j = 1, 3
       if(newloadcase%stress%maskLogical(i,j)) then
         write(6,'(2x,f12.7)',advance='no') newloadcase%stress%values(i,j)*1e-9_preal
       else
         write(6,'(2x,12a)',advance='no') '     *      '
       endif
       enddo; write(6,'(/)',advance='no')
     enddo
     if (any(abs(matmul(newloadcase%rot%asMatrix(), &
                        transpose(newloadcase%rot%asMatrix()))-math_i3) > &
                reshape(spread(tol_math_check,1,9),[ 3,3]))) errorid = 846                          ! given rotation matrix contains strain
     if (any(dneq(newloadcase%rot%asMatrix(), math_i3))) &
       write(6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',&
                transpose(newloadcase%rot%asMatrix())
     if (newloadcase%time < 0.0_preal) errorid = 834                                                ! negative time increment
     write(6,'(2x,a,f0.3)') 'time: ', newloadcase%time
     if (newloadcase%incs < 1)    errorid = 835                                                     ! non-positive incs count
     write(6,'(2x,a,i0)')  'increments: ', newloadcase%incs
     if (newloadcase%outputfrequency < 1)  errorid = 836                                            ! non-positive result frequency
     write(6,'(2x,a,i0)')  'output frequency: ', newloadcase%outputfrequency
     if (newloadcase%restartfrequency < 1)  errorid = 839                                           ! non-positive restart frequency
     if (newloadcase%restartfrequency < huge(0)) &
       write(6,'(2x,a,i0)')  'restart frequency: ', newloadcase%restartfrequency
     if (errorid > 0) call io_error(error_id = errorid, ext_msg = loadcase_string)                  ! exit with error message
   endif reportandcheck
   loadcases = [loadcases,newloadcase]                                                              ! load case is ok, append it
 enddo
 
 
!--------------------------------------------------------------------------------------------------
! doing initialization depending on active solvers
 call utilities_init
 do field = 1, nactivefields
   select case (loadcases(1)%ID(field))
     case(field_mech_id)
       call mech_init
     
     case(field_thermal_id)
       call grid_thermal_spectral_init
 
     case(field_damage_id)
       call grid_damage_spectral_init
 
   end select
 enddo
 
!--------------------------------------------------------------------------------------------------
! write header of output file
 if (worldrank == 0) then
   writeheader: if (interface_restartinc < 1) then
     open(newunit=statunit,file=trim(getsolverjobname())//'.sta',form='FORMATTED',status='REPLACE')
     write(statunit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded'                 ! statistics file
     if (iand(debug_level(debug_spectral),debug_levelbasic) /= 0) &
       write(6,'(/,a)') ' header of statistics file written out'
     flush(6)
   else writeheader
     open(newunit=statunit,file=trim(getsolverjobname())//&
                                 '.sta',form='FORMATTED', position='APPEND', status='OLD')
   endif writeheader
 endif
 
 writeundeformed: if (interface_restartinc < 1) then
   write(6,'(1/,a)') ' ... writing initial configuration to file ........................'
   call cpfem_results(0,0.0_preal)
 endif writeundeformed
 
 loadcaselooping: do currentloadcase = 1, size(loadcases)
   time0 = time                                                                                     ! load case start time
   guess = loadcases(currentloadcase)%followFormerTrajectory                                        ! change of load case? homogeneous guess for the first inc
 
   inclooping: do inc = 1, loadcases(currentloadcase)%incs
     totalincscounter = totalincscounter + 1
 
!--------------------------------------------------------------------------------------------------
! forwarding time
     timeincold = timeinc                                                                           ! last timeinc that brought former inc to an end
     if (loadcases(currentloadcase)%logscale == 0) then                                             ! linear scale
       timeinc = loadcases(currentloadcase)%time/real(loadcases(currentloadcase)%incs,preal)
     else
       if (currentloadcase == 1) then                                                               ! 1st load case of logarithmic scale
         if (inc == 1) then                                                                         ! 1st inc of 1st load case of logarithmic scale
           timeinc = loadcases(1)%time*(2.0_preal**real(    1-loadcases(1)%incs ,preal))            ! assume 1st inc is equal to 2nd
         else                                                                                       ! not-1st inc of 1st load case of logarithmic scale
           timeinc = loadcases(1)%time*(2.0_preal**real(inc-1-loadcases(1)%incs ,preal))
         endif
       else                                                                                         ! not-1st load case 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  ,preal)/&
                                                    real(loadcases(currentloadcase)%incs ,preal)))
       endif
     endif
     timeinc = timeinc * real(substepfactor,preal)**real(-cutbacklevel,preal)                       ! depending on cut back level, decrease time step
 
     skipping: if (totalincscounter <= interface_restartinc) then                                   ! not yet at restart inc?
       time = time + timeinc                                                                        ! just advance time, skip already performed calculation
       guess = .true.                                                                               ! QUESTION:why forced guessing instead of inheriting loadcase preference
     else skipping
       stepfraction = 0                                                                             ! fraction scaled by stepFactor**cutLevel
 
       substeplooping: do while (stepfraction < substepfactor**cutbacklevel)
         remainingloadcasetime = loadcases(currentloadcase)%time+time0 - time
         time = time + timeinc                                                                      ! forward target time
         stepfraction = stepfraction + 1                                                            ! count step
 
!--------------------------------------------------------------------------------------------------
! report begin of new step
         write(6,'(/,a)') ' ###########################################################################'
         write(6,'(1x,a,es12.5,6(a,i0))') &
                 'Time', time, &
                 's: Increment ', inc,'/',loadcases(currentloadcase)%incs,&
                 '-', stepfraction,'/',substepfactor**cutbacklevel,&
                 ' of load case ', currentloadcase,'/',size(loadcases)
         write(incinfo,'(4(a,i0))') &
                 'Increment ',totalincscounter,'/',sum(loadcases%incs),&
                 '-', stepfraction,'/',substepfactor**cutbacklevel
         flush(6)
 
!--------------------------------------------------------------------------------------------------
! forward fields
         do field = 1, nactivefields
           select case(loadcases(currentloadcase)%ID(field))
             case(field_mech_id)
               call mech_forward (&
                       cutback,guess,timeinc,timeincold,remainingloadcasetime, &
                       deformation_bc     = loadcases(currentloadcase)%deformation, &
                       stress_bc          = loadcases(currentloadcase)%stress, &
                       rotation_bc        = loadcases(currentloadcase)%rot)
 
             case(field_thermal_id); call grid_thermal_spectral_forward(cutback)
             case(field_damage_id);  call grid_damage_spectral_forward(cutback)
           end select
         enddo
         if(.not. cutback) call cpfem_forward
 
!--------------------------------------------------------------------------------------------------
! solve fields
         stagiter = 0
         stagiterate = .true.
         do while (stagiterate)
           do field = 1, nactivefields
             select case(loadcases(currentloadcase)%ID(field))
               case(field_mech_id)
                 solres(field) = mech_solution(&
                                        incinfo,timeinc,timeincold, &
                                        stress_bc   = loadcases(currentloadcase)%stress, &
                                        rotation_bc = loadcases(currentloadcase)%rot)
 
               case(field_thermal_id)
                 solres(field) = grid_thermal_spectral_solution(timeinc,timeincold)
 
               case(field_damage_id)
                 solres(field) = grid_damage_spectral_solution(timeinc,timeincold)
 
             end select
 
             if (.not. solres(field)%converged) exit                                                ! no solution found
 
           enddo
           stagiter = stagiter + 1
           stagiterate =            stagiter < stagitmax &
                        .and.       all(solres(:)%converged) &
                        .and. .not. all(solres(:)%stagConverged)                                    ! stationary with respect to staggered iteration
         enddo
 
!--------------------------------------------------------------------------------------------------
! check solution for either advance or retry
 
         if ( (all(solres(:)%converged .and. solres(:)%stagConverged)) &                            ! converged
              .and. .not. solres(1)%termIll) then                                                   ! and acceptable solution found
           call mech_updatecoords
           timeincold = timeinc
           cutback = .false.
           guess = .true.                                                                           ! start guessing after first converged (sub)inc
           if (worldrank == 0) then
             write(statunit,*) totalincscounter, time, cutbacklevel, &
                               solres%converged, solres%iterationsNeeded
             flush(statunit)
           endif
         elseif (cutbacklevel < maxcutback) then                                                    ! further cutbacking tolerated?
           cutback = .true.
           stepfraction = (stepfraction - 1) * substepfactor                                        ! adjust to new denominator
           cutbacklevel = cutbacklevel + 1
           time    = time - timeinc                                                                 ! rewind time
           timeinc = timeinc/real(substepfactor,preal)                                              ! cut timestep
           write(6,'(/,a)') ' cutting back '
         else                                                                                       ! no more options to continue
           call io_warning(850)
           if (worldrank == 0) close(statunit)
           call quit(0)                                                                             ! quit
         endif
 
       enddo substeplooping
 
       cutbacklevel = max(0, cutbacklevel - 1)                                                      ! try half number of subincs next inc
 
       if (all(solres(:)%converged)) then
         write(6,'(/,a,i0,a)') ' increment ', totalincscounter, ' converged'
       else
         write(6,'(/,a,i0,a)') ' increment ', totalincscounter, ' NOT converged'
       endif; flush(6)
 
       if (mod(inc,loadcases(currentloadcase)%outputFrequency) == 0) then                           ! at output frequency
         write(6,'(1/,a)') ' ... writing results to file ......................................'
         flush(6)
         call cpfem_results(totalincscounter,time)
       endif
       if (mod(inc,loadcases(currentloadcase)%restartFrequency) == 0) then
         call mech_restartwrite
         call cpfem_restartwrite
       endif
     endif skipping
 
    enddo inclooping
 
 enddo loadcaselooping
 
 
!--------------------------------------------------------------------------------------------------
! report summary of whole calculation
 write(6,'(/,a)') ' ###########################################################################'
 if (worldrank == 0) close(statunit)
 
 call quit(0)                                                                                       ! no complains ;)
 
end program damask_spectral