DAMASK_EICMD/code/numerics.f90

! Copyright 2011-13 Max-Planck-Institut für Eisenforschung GmbH
!
! This file is part of DAMASK,
! the Düsseldorf Advanced MAterial Simulation Kit.
!
! DAMASK is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! DAMASK is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with DAMASK. If not, see <http://www.gnu.org/licenses/>.
!
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Managing of parameters related to numerics
!--------------------------------------------------------------------------------------------------
module numerics
 use prec, only: &
   pInt, &
   pReal

 implicit none
 private
 character(len=64), parameter, private :: &
   numerics_CONFIGFILE        = 'numerics.config'                                                   !< name of configuration file

 integer(pInt), protected, public :: &
   iJacoStiffness             =  1_pInt, &                                                          !< frequency of stiffness update
   iJacoLpresiduum            =  1_pInt, &                                                          !< frequency of Jacobian update of residuum in Lp
   nHomog                     = 20_pInt, &                                                          !< homogenization loop limit (only for debugging info, loop limit is determined by "subStepMinHomog")
   nMPstate                   = 10_pInt, &                                                          !< materialpoint state loop limit
   nCryst                     = 20_pInt, &                                                          !< crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst")
   nState                     = 10_pInt, &                                                          !< state loop limit
   nStress                    = 40_pInt, &                                                          !< stress loop limit
   pert_method                =  1_pInt, &                                                          !< method used in perturbation technique for tangent
   fixedSeed                  =  0_pInt, &                                                          !< fixed seeding for pseudo-random number generator, Default 0: use random seed
   DAMASK_NumThreadsInt       =  0_pInt                                                             !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive
 integer(pInt), public :: &
   numerics_integrationMode   =  0_pInt                                                             !< integrationMode 1 = central solution ; integrationMode 2 = perturbation, Default 0: undefined, is not read from file
 integer(pInt), dimension(2) , protected, public :: &
   numerics_integrator        =  1_pInt                                                             !< method used for state integration (central & perturbed state), Default 1: fix-point iteration for both states
 real(pReal), protected, public :: &
   relevantStrain             =  1.0e-7_pReal, &                                                    !< strain increment considered significant (used by crystallite to determine whether strain inc is considered significant)
   defgradTolerance           =  1.0e-7_pReal, &                                                    !< deviation of deformation gradient that is still allowed (used by CPFEM to determine outdated ffn1)
   pert_Fg                    =  1.0e-7_pReal, &                                                    !< strain perturbation for FEM Jacobi
   subStepMinCryst            =  1.0e-3_pReal, &                                                    !< minimum (relative) size of sub-step allowed during cutback in crystallite
   subStepMinHomog            =  1.0e-3_pReal, &                                                    !< minimum (relative) size of sub-step allowed during cutback in homogenization
   subStepSizeCryst           =  0.25_pReal, &                                                      !< size of first substep when cutback in crystallite
   subStepSizeHomog           =  0.25_pReal, &                                                      !< size of first substep when cutback in homogenization
   stepIncreaseCryst          =  1.5_pReal, &                                                       !< increase of next substep size when previous substep converged in crystallite
   stepIncreaseHomog          =  1.5_pReal, &                                                       !< increase of next substep size when previous substep converged in homogenization
   rTol_crystalliteState      =  1.0e-6_pReal, &                                                    !< relative tolerance in crystallite state loop 
   rTol_crystalliteTemperature=  1.0e-6_pReal, &                                                    !< relative tolerance in crystallite temperature loop 
   rTol_crystalliteStress     =  1.0e-6_pReal, &                                                    !< relative tolerance in crystallite stress loop
   aTol_crystalliteStress     =  1.0e-8_pReal, &                                                    !< absolute tolerance in crystallite stress loop, Default 1.0e-8: residuum is in Lp and hence strain is on this order
   numerics_unitlength        =  1.0_pReal, &                                                       !< determines the physical length of one computational length unit
   absTol_RGC                 =  1.0e+4_pReal, &                                                    !< absolute tolerance of RGC residuum
   relTol_RGC                 =  1.0e-3_pReal, &                                                    !< relative tolerance of RGC residuum
   absMax_RGC                 =  1.0e+10_pReal, &                                                   !< absolute maximum of RGC residuum
   relMax_RGC                 =  1.0e+2_pReal, &                                                    !< relative maximum of RGC residuum
   pPert_RGC                  =  1.0e-7_pReal, &                                                    !< perturbation for computing RGC penalty tangent
   xSmoo_RGC                  =  1.0e-5_pReal, &                                                    !< RGC penalty smoothing parameter (hyperbolic tangent)
   viscPower_RGC              =  1.0e+0_pReal, &                                                    !< power (sensitivity rate) of numerical viscosity in RGC scheme, Default 1.0e0: Newton viscosity (linear model)
   viscModus_RGC              =  0.0e+0_pReal, &                                                    !< stress modulus of RGC numerical viscosity, Default 0.0e0: No viscosity is applied
   refRelaxRate_RGC           =  1.0e-3_pReal, &                                                    !< reference relaxation rate in RGC viscosity
   maxdRelax_RGC              =  1.0e+0_pReal, &                                                    !< threshold of maximum relaxation vector increment (if exceed this then cutback)
   maxVolDiscr_RGC            =  1.0e-5_pReal, &                                                    !< threshold of maximum volume discrepancy allowed
   volDiscrMod_RGC            =  1.0e+12_pReal, &                                                   !< stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint)
   volDiscrPow_RGC            =  5.0_pReal                                                          !< powerlaw penalty for volume discrepancy
 logical, protected, public :: &                                                   
   analyticJaco               = .false., &                                                         !< use analytic Jacobian or perturbation, Default .false.: calculate Jacobian using perturbations
   usePingPong                = .true., & 
   numerics_timeSyncing       = .false.                                                            !< flag indicating if time synchronization in crystallite is used for nonlocal plasticity

!--------------------------------------------------------------------------------------------------
! spectral parameters:
#ifdef Spectral
 real(pReal), protected, public :: &
   err_div_tolAbs             =  1.0e-10_pReal, &                                                   !< absolute tolerance for equilibrium
   err_div_tolRel             =  5.0e-4_pReal, &                                                    !< relative tolerance for equilibrium
   err_curl_tolAbs            =  1.0e-10_pReal, &                                                   !< absolute tolerance for compatibility
   err_curl_tolRel            =  5.0e-4_pReal, &                                                    !< relative tolerance for compatibility
   err_stress_tolAbs          =  1.0e3_pReal,  &                                                    !< absolute tolerance for fullfillment of stress BC
   err_stress_tolRel          =  0.01_pReal, &                                                      !< relative tolerance for fullfillment of stress BC
   fftw_timelimit             = -1.0_pReal, &                                                       !< sets the timelimit of plan creation for FFTW, see manual on www.fftw.org, Default -1.0: disable timelimit
   rotation_tol               =  1.0e-12_pReal, &                                                   !< tolerance of rotation specified in loadcase, Default 1.0e-12: first guess
   polarAlpha                 =  1.0_pReal, &                                                       !< polarization scheme parameter 0.0 < alpha < 2.0. alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme 
   polarBeta                  =  1.0_pReal                                                          !< polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme 
 character(len=64), private :: &
   fftw_plan_mode             = 'FFTW_PATIENT'                                                      !< reads the planing-rigor flag, see manual on www.fftw.org, Default FFTW_PATIENT: use patient planner flag
 character(len=64), protected, public :: & 
   myspectralsolver           = 'basic'  , &                                                        !< spectral solution method 
   myfilter                   = 'none'                                                              !< spectral filtering method
 character(len=1024), protected, public :: &
   petsc_options              = '-snes_type ngmres &
                                &-snes_ngmres_anderson '
 integer(pInt), protected, public :: &
   fftw_planner_flag          =  32_pInt, &                                                         !< conversion of fftw_plan_mode to integer, basically what is usually done in the include file of fftw
   itmax                      =  250_pInt, &                                                        !< maximum number of iterations
   itmin                      =  2_pInt, &                                                          !< minimum number of iterations
   maxCutBack                 =  3_pInt, &                                                          !< max number of cut backs
   regridMode                 =  0_pInt, &                                                          !< 0: no regrid; 1: regrid if DAMASK doesn't converge; 2: regrid if DAMASK or BVP Solver doesn't converge 
   divergence_correction      =  2_pInt                                                             !< correct divergence calculation in fourier space 0: no correction, 1: size scaled to 1, 2: size scaled to Npoints
 logical, protected, public :: &
   memory_efficient           = .true., &                                                           !< for fast execution (pre calculation of gamma_hat), Default .true.: do not precalculate
   update_gamma               = .false.                                                             !< update gamma operator with current stiffness, Default .false.: use initial stiffness 
#endif

 public :: numerics_init
  
contains


!--------------------------------------------------------------------------------------------------
!> @brief reads in parameters from numerics.config and sets openMP related parameters. Also does
! a sanity check
!--------------------------------------------------------------------------------------------------
subroutine numerics_init
 use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
 use IO, only: &
   IO_read, &
   IO_error, &
   IO_open_file_stat, &
   IO_isBlank, &
   IO_stringPos, &
   IO_stringValue, &
   IO_lc, &
   IO_floatValue, &
   IO_intValue, &
   IO_warning, &
   IO_timeStamp

#ifdef Spectral
!$ use OMP_LIB, only: omp_set_num_threads                                                           ! Use the standard conforming module file for omp if not using MSC.Marc
#endif  
 implicit none
#ifndef Spectral
!$ include "omp_lib.h"                                                                              ! use the not F90 standard conforming include file to prevent crashes with some versions of MSC.Marc
#endif  
 integer(pInt), parameter ::                 FILEUNIT = 300_pInt ,&
                                             maxNchunks = 2_pInt
!$ integer ::                                gotDAMASK_NUM_THREADS = 1
 integer :: i                                                                                       ! no pInt
 integer(pInt), dimension(1+2*maxNchunks) :: positions
 character(len=65536) :: &
   tag ,&
   line
!$ character(len=6) DAMASK_NumThreadsString                                                         ! environment variable DAMASK_NUM_THREADS

 write(6,'(/,a)') ' <<<+-  numerics init  -+>>>'
 write(6,'(a)')   ' $Id$'
 write(6,'(a16,a)')   ' Current time : ',IO_timeStamp()
#include "compilation_info.f90"

!$ call GET_ENVIRONMENT_VARIABLE(NAME='DAMASK_NUM_THREADS',VALUE=DAMASK_NumThreadsString,STATUS=gotDAMASK_NUM_THREADS)   ! get environment variable DAMASK_NUM_THREADS...
!$ if(gotDAMASK_NUM_THREADS /= 0) & 
!$   call IO_warning(35_pInt,ext_msg='BEGIN:'//DAMASK_NumThreadsString//':END')
!$ read(DAMASK_NumThreadsString,'(i6)') DAMASK_NumThreadsInt                                        ! ...convert it to integer...
!$ if (DAMASK_NumThreadsInt < 1_pInt) DAMASK_NumThreadsInt = 1_pInt                                 ! ...ensure that its at least one...
!$ call omp_set_num_threads(DAMASK_NumThreadsInt)                                                   ! ...and use it as number of threads for parallel execution

!--------------------------------------------------------------------------------------------------
! try to open the config file
 fileExists: if(IO_open_file_stat(fileunit,numerics_configFile)) then 
   write(6,'(a,/)') ' using values from config file'
    
!--------------------------------------------------------------------------------------------------
! read variables from config file and overwrite default parameters if keyword is present
   line = ''
   do while (trim(line) /= '#EOF#')                                                                 ! read thru sections of phase part
     line = IO_read(fileunit)
     do i=1,len(line)
       if(line(i:i) == '=') line(i:i) = ' '                                                         ! also allow keyword = value version
     enddo
     if (IO_isBlank(line)) cycle                                                                    ! skip empty lines
     positions = IO_stringPos(line,maxNchunks)
     tag = IO_lc(IO_stringValue(line,positions,1_pInt))                                             ! extract key
     select case(tag)
       case ('relevantstrain')
         relevantStrain = IO_floatValue(line,positions,2_pInt)
       case ('defgradtolerance')
         defgradTolerance = IO_floatValue(line,positions,2_pInt)
       case ('ijacostiffness')
         iJacoStiffness = IO_intValue(line,positions,2_pInt)
       case ('ijacolpresiduum')
         iJacoLpresiduum = IO_intValue(line,positions,2_pInt)
       case ('pert_fg')
         pert_Fg = IO_floatValue(line,positions,2_pInt)
       case ('pert_method')
         pert_method = IO_intValue(line,positions,2_pInt)
       case ('nhomog')
         nHomog = IO_intValue(line,positions,2_pInt)
       case ('nmpstate')
         nMPstate = IO_intValue(line,positions,2_pInt)
       case ('ncryst')
         nCryst = IO_intValue(line,positions,2_pInt)
       case ('nstate')
         nState = IO_intValue(line,positions,2_pInt)
       case ('nstress')
         nStress = IO_intValue(line,positions,2_pInt)
       case ('substepmincryst')
         subStepMinCryst = IO_floatValue(line,positions,2_pInt)
       case ('substepsizecryst')
         subStepSizeCryst = IO_floatValue(line,positions,2_pInt)
       case ('stepincreasecryst')
         stepIncreaseCryst = IO_floatValue(line,positions,2_pInt)
       case ('substepminhomog')
         subStepMinHomog = IO_floatValue(line,positions,2_pInt)
       case ('substepsizehomog')
         subStepSizeHomog = IO_floatValue(line,positions,2_pInt)
       case ('stepincreasehomog')
         stepIncreaseHomog = IO_floatValue(line,positions,2_pInt)
       case ('rtol_crystallitestate')
         rTol_crystalliteState = IO_floatValue(line,positions,2_pInt)
       case ('rtol_crystallitetemperature')
         rTol_crystalliteTemperature = IO_floatValue(line,positions,2_pInt)
       case ('rtol_crystallitestress')
         rTol_crystalliteStress = IO_floatValue(line,positions,2_pInt)
       case ('atol_crystallitestress')
         aTol_crystalliteStress = IO_floatValue(line,positions,2_pInt)
       case ('integrator')
         numerics_integrator(1) = IO_intValue(line,positions,2_pInt)
       case ('integratorstiffness')
         numerics_integrator(2) = IO_intValue(line,positions,2_pInt)
       case ('analyticjaco')
         analyticJaco = IO_intValue(line,positions,2_pInt) > 0_pInt
       case ('usepingpong')
         usepingpong = IO_intValue(line,positions,2_pInt) > 0_pInt
       case ('timesyncing')
         numerics_timeSyncing = IO_intValue(line,positions,2_pInt) > 0_pInt
       case ('unitlength')
         numerics_unitlength = IO_floatValue(line,positions,2_pInt)

!--------------------------------------------------------------------------------------------------
! RGC parameters
       case ('atol_rgc')
         absTol_RGC = IO_floatValue(line,positions,2_pInt)
       case ('rtol_rgc')
         relTol_RGC = IO_floatValue(line,positions,2_pInt)
       case ('amax_rgc')
         absMax_RGC = IO_floatValue(line,positions,2_pInt)
       case ('rmax_rgc')
         relMax_RGC = IO_floatValue(line,positions,2_pInt)
       case ('perturbpenalty_rgc')
         pPert_RGC = IO_floatValue(line,positions,2_pInt)
       case ('relevantmismatch_rgc')
         xSmoo_RGC = IO_floatValue(line,positions,2_pInt)
       case ('viscositypower_rgc')
         viscPower_RGC = IO_floatValue(line,positions,2_pInt)
       case ('viscositymodulus_rgc')
         viscModus_RGC = IO_floatValue(line,positions,2_pInt)
       case ('refrelaxationrate_rgc')
         refRelaxRate_RGC = IO_floatValue(line,positions,2_pInt)
       case ('maxrelaxation_rgc')
         maxdRelax_RGC = IO_floatValue(line,positions,2_pInt)
       case ('maxvoldiscrepancy_rgc')
         maxVolDiscr_RGC = IO_floatValue(line,positions,2_pInt)
       case ('voldiscrepancymod_rgc')
         volDiscrMod_RGC = IO_floatValue(line,positions,2_pInt)
       case ('discrepancypower_rgc')
         volDiscrPow_RGC = IO_floatValue(line,positions,2_pInt)

!--------------------------------------------------------------------------------------------------
! random seeding parameter
       case ('fixed_seed')
         fixedSeed = IO_intValue(line,positions,2_pInt)

!--------------------------------------------------------------------------------------------------
! spectral parameters
#ifdef Spectral
       case ('err_div_tolabs')
         err_div_tolAbs = IO_floatValue(line,positions,2_pInt)
       case ('err_div_tolrel')
         err_div_tolRel = IO_floatValue(line,positions,2_pInt)
       case ('err_stress_tolrel')
         err_stress_tolrel = IO_floatValue(line,positions,2_pInt)
       case ('err_stress_tolabs')
         err_stress_tolabs = IO_floatValue(line,positions,2_pInt)
       case ('itmax')
         itmax = IO_intValue(line,positions,2_pInt)
       case ('itmin')
         itmin = IO_intValue(line,positions,2_pInt)
       case ('maxcutback')
         maxCutBack = IO_intValue(line,positions,2_pInt)
       case ('regridmode')
         regridMode = IO_intValue(line,positions,2_pInt)
       case ('memory_efficient')
         memory_efficient = IO_intValue(line,positions,2_pInt)  > 0_pInt
       case ('fftw_timelimit')
         fftw_timelimit = IO_floatValue(line,positions,2_pInt)
       case ('fftw_plan_mode')
         fftw_plan_mode = IO_lc(IO_stringValue(line,positions,2_pInt))
       case ('myfilter')
         myfilter = IO_lc(IO_stringValue(line,positions,2_pInt))
       case ('rotation_tol')
         rotation_tol = IO_floatValue(line,positions,2_pInt)
       case ('divergence_correction')
         divergence_correction = IO_intValue(line,positions,2_pInt)
       case ('update_gamma')
         update_gamma = IO_intValue(line,positions,2_pInt)  > 0_pInt
#ifdef PETSc
       case ('petsc_options')
         petsc_options = trim(line(positions(4):))
       case ('myspectralsolver')
         myspectralsolver = IO_lc(IO_stringValue(line,positions,2_pInt))
       case ('err_curl_tolabs')
         err_curl_tolAbs = IO_floatValue(line,positions,2_pInt)
       case ('err_curl_tolrel')
         err_curl_tolRel = IO_floatValue(line,positions,2_pInt)
       case ('polaralpha')
         polarAlpha = IO_floatValue(line,positions,2_pInt)
       case ('polarbeta')
         polarBeta = IO_floatValue(line,positions,2_pInt)
#endif 
#ifndef PETSc
       case ('myspectralsolver', 'petsc_options', &
             'err_curl_tolabs','err_curl_tolrel','polaralpha','polarBeta')                          ! found PETSc parameter, but compiled without PETSc
         call IO_warning(41_pInt,ext_msg=tag)
#endif
#endif
#ifndef Spectral
      case ('err_div_tolabs','err_div_tolrel','err_stress_tolrel','err_stress_tolabs',&             ! found spectral parameter for FEM build
            'itmax', 'itmin','memory_efficient','fftw_timelimit','fftw_plan_mode', &
            'rotation_tol','divergence_correction','update_gamma','myfilter', &
            'err_curl_tolabs','err_curl_tolrel', &
            'maxcutback','polaralpha','polarbeta')
         call IO_warning(40_pInt,ext_msg=tag)
#endif
       case default                                                                                ! found unknown keyword
         call IO_error(300_pInt,ext_msg=tag)
     endselect
   enddo
   close(fileunit)

 else fileExists
   write(6,'(a,/)') ' using standard values'
 endif fileExists

#ifdef Spectral
 select case(IO_lc(fftw_plan_mode))                                                                ! setting parameters for the plan creation of FFTW. Basically a translation from fftw3.f
   case('estimate','fftw_estimate')                                                                 ! ordered from slow execution (but fast plan creation) to fast execution
     fftw_planner_flag = 64_pInt
   case('measure','fftw_measure')
     fftw_planner_flag = 0_pInt
   case('patient','fftw_patient')
     fftw_planner_flag= 32_pInt
   case('exhaustive','fftw_exhaustive')
     fftw_planner_flag = 8_pInt 
   case default
     call IO_warning(warning_ID=47_pInt,ext_msg=trim(IO_lc(fftw_plan_mode)))
     fftw_planner_flag = 32_pInt
 end select
#endif

 numerics_timeSyncing = numerics_timeSyncing .and. all(numerics_integrator==2_pInt)                 ! timeSyncing only allowed for explicit Euler integrator

!--------------------------------------------------------------------------------------------------
! writing parameters to output
 write(6,'(a24,1x,es8.1)')  ' relevantStrain:         ',relevantStrain
 write(6,'(a24,1x,es8.1)')  ' defgradTolerance:       ',defgradTolerance
 write(6,'(a24,1x,i8)')     ' iJacoStiffness:         ',iJacoStiffness
 write(6,'(a24,1x,i8)')     ' iJacoLpresiduum:        ',iJacoLpresiduum
 write(6,'(a24,1x,es8.1)')  ' pert_Fg:                ',pert_Fg
 write(6,'(a24,1x,i8)')     ' pert_method:            ',pert_method
 write(6,'(a24,1x,i8)')     ' nCryst:                 ',nCryst
 write(6,'(a24,1x,es8.1)')  ' subStepMinCryst:        ',subStepMinCryst
 write(6,'(a24,1x,es8.1)')  ' subStepSizeCryst:       ',subStepSizeCryst
 write(6,'(a24,1x,es8.1)')  ' stepIncreaseCryst:      ',stepIncreaseCryst
 write(6,'(a24,1x,i8)')     ' nState:                 ',nState
 write(6,'(a24,1x,i8)')     ' nStress:                ',nStress
 write(6,'(a24,1x,es8.1)')  ' rTol_crystalliteState:  ',rTol_crystalliteState
 write(6,'(a24,1x,es8.1)')  ' rTol_crystalliteTemp:   ',rTol_crystalliteTemperature
 write(6,'(a24,1x,es8.1)')  ' rTol_crystalliteStress: ',rTol_crystalliteStress
 write(6,'(a24,1x,es8.1)')  ' aTol_crystalliteStress: ',aTol_crystalliteStress
 write(6,'(a24,2(1x,i8))')  ' integrator:             ',numerics_integrator
 write(6,'(a24,1x,L8)')     ' timeSyncing:            ',numerics_timeSyncing
 write(6,'(a24,1x,L8)')     ' analytic Jacobian:      ',analyticJaco
 write(6,'(a24,1x,L8)')     ' use ping pong scheme:   ',usepingpong
 write(6,'(a24,1x,es8.1,/)')' unitlength:             ',numerics_unitlength

 write(6,'(a24,1x,i8)')     ' nHomog:                 ',nHomog
 write(6,'(a24,1x,es8.1)')  ' subStepMinHomog:        ',subStepMinHomog
 write(6,'(a24,1x,es8.1)')  ' subStepSizeHomog:       ',subStepSizeHomog
 write(6,'(a24,1x,es8.1)')  ' stepIncreaseHomog:      ',stepIncreaseHomog
 write(6,'(a24,1x,i8,/)')   ' nMPstate:               ',nMPstate

!--------------------------------------------------------------------------------------------------
! RGC parameters
 write(6,'(a24,1x,es8.1)')   ' aTol_RGC:               ',absTol_RGC
 write(6,'(a24,1x,es8.1)')   ' rTol_RGC:               ',relTol_RGC
 write(6,'(a24,1x,es8.1)')   ' aMax_RGC:               ',absMax_RGC
 write(6,'(a24,1x,es8.1)')   ' rMax_RGC:               ',relMax_RGC
 write(6,'(a24,1x,es8.1)')   ' perturbPenalty_RGC:     ',pPert_RGC
 write(6,'(a24,1x,es8.1)')   ' relevantMismatch_RGC:   ',xSmoo_RGC
 write(6,'(a24,1x,es8.1)')   ' viscosityrate_RGC:      ',viscPower_RGC
 write(6,'(a24,1x,es8.1)')   ' viscositymodulus_RGC:   ',viscModus_RGC
 write(6,'(a24,1x,es8.1)')   ' maxrelaxation_RGC:      ',maxdRelax_RGC
 write(6,'(a24,1x,es8.1)')   ' maxVolDiscrepancy_RGC:  ',maxVolDiscr_RGC
 write(6,'(a24,1x,es8.1)')   ' volDiscrepancyMod_RGC:  ',volDiscrMod_RGC
 write(6,'(a24,1x,es8.1,/)') ' discrepancyPower_RGC:   ',volDiscrPow_RGC

!--------------------------------------------------------------------------------------------------
! Random seeding parameter
 write(6,'(a24,1x,i16,/)')    ' fixed_seed:             ',fixedSeed
 if (fixedSeed <= 0_pInt) &
   write(6,'(a,/)') ' No fixed Seed: Random is random!'
!--------------------------------------------------------------------------------------------------
! openMP parameter
  !$  write(6,'(a24,1x,i8,/)')   ' number of threads:      ',DAMASK_NumThreadsInt

!--------------------------------------------------------------------------------------------------
! spectral parameters
#ifdef Spectral
 write(6,'(a24,1x,es8.1)')   ' err_div_tolAbs:         ',err_div_tolAbs
 write(6,'(a24,1x,es8.1)')   ' err_div_tolRel:         ',err_div_tolRel
 write(6,'(a24,1x,es8.1)')   ' err_stress_tolrel:      ',err_stress_tolrel
 write(6,'(a24,1x,es8.1)')   ' err_stress_tolabs:      ',err_stress_tolabs

 write(6,'(a24,1x,i8)')      ' itmax:                  ',itmax
 write(6,'(a24,1x,i8)')      ' itmin:                  ',itmin
 write(6,'(a24,1x,i8)')      ' maxCutBack:             ',maxCutBack
 write(6,'(a24,1x,i8)')      ' regridMode:             ',regridMode
 write(6,'(a24,1x,L8)')      ' memory_efficient:       ',memory_efficient
 if(fftw_timelimit<0.0_pReal) then
   write(6,'(a24,1x,L8)')    ' fftw_timelimit:         ',.false.
 else    
   write(6,'(a24,1x,es8.1)') ' fftw_timelimit:         ',fftw_timelimit
 endif
 write(6,'(a24,1x,a)')       ' fftw_plan_mode:         ',trim(fftw_plan_mode)

 write(6,'(a24,1x,a)')       ' myfilter:               ',trim(myfilter)
 write(6,'(a24,1x,i8)')      ' fftw_planner_flag:      ',fftw_planner_flag
 write(6,'(a24,1x,es8.1)')   ' rotation_tol:           ',rotation_tol
 write(6,'(a24,1x,i8)')      ' divergence_correction:  ',divergence_correction
 write(6,'(a24,1x,L8,/)')    ' update_gamma:           ',update_gamma
#ifdef PETSc
 write(6,'(a24,1x,es8.1)')   ' err_curl_tolAbs:        ',err_curl_tolAbs
 write(6,'(a24,1x,es8.1)')   ' err_curl_tolRel:        ',err_curl_tolRel
 write(6,'(a24,1x,es8.1)')   ' polarAlpha:             ',polarAlpha
 write(6,'(a24,1x,es8.1)')   ' polarBeta:              ',polarBeta
 write(6,'(a24,1x,a)')       ' myspectralsolver:       ',trim(myspectralsolver)
 write(6,'(a24,1x,a)')       ' PETSc_options:          ',trim(petsc_options)
#endif
#endif

!--------------------------------------------------------------------------------------------------
! sanity checks
 if (relevantStrain <= 0.0_pReal)          call IO_error(301_pInt,ext_msg='relevantStrain')
 if (defgradTolerance <= 0.0_pReal)        call IO_error(301_pInt,ext_msg='defgradTolerance')
 if (iJacoStiffness < 1_pInt)              call IO_error(301_pInt,ext_msg='iJacoStiffness')
 if (iJacoLpresiduum < 1_pInt)             call IO_error(301_pInt,ext_msg='iJacoLpresiduum')
 if (pert_Fg <= 0.0_pReal)                 call IO_error(301_pInt,ext_msg='pert_Fg')
 if (pert_method <= 0_pInt .or. pert_method >= 4_pInt) &
                                           call IO_error(301_pInt,ext_msg='pert_method')
 if (nHomog < 1_pInt)                      call IO_error(301_pInt,ext_msg='nHomog')
 if (nMPstate < 1_pInt)                    call IO_error(301_pInt,ext_msg='nMPstate')
 if (nCryst < 1_pInt)                      call IO_error(301_pInt,ext_msg='nCryst')
 if (nState < 1_pInt)                      call IO_error(301_pInt,ext_msg='nState')
 if (nStress < 1_pInt)                     call IO_error(301_pInt,ext_msg='nStress')
 if (subStepMinCryst <= 0.0_pReal)         call IO_error(301_pInt,ext_msg='subStepMinCryst')
 if (subStepSizeCryst <= 0.0_pReal)        call IO_error(301_pInt,ext_msg='subStepSizeCryst')
 if (stepIncreaseCryst <= 0.0_pReal)       call IO_error(301_pInt,ext_msg='stepIncreaseCryst')
 if (subStepMinHomog <= 0.0_pReal)         call IO_error(301_pInt,ext_msg='subStepMinHomog')
 if (subStepSizeHomog <= 0.0_pReal)        call IO_error(301_pInt,ext_msg='subStepSizeHomog')
 if (stepIncreaseHomog <= 0.0_pReal)       call IO_error(301_pInt,ext_msg='stepIncreaseHomog')
 if (rTol_crystalliteState <= 0.0_pReal)   call IO_error(301_pInt,ext_msg='rTol_crystalliteState')
 if (rTol_crystalliteTemperature <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteTemperature')
 if (rTol_crystalliteStress <= 0.0_pReal)  call IO_error(301_pInt,ext_msg='rTol_crystalliteStress')
 if (aTol_crystalliteStress <= 0.0_pReal)  call IO_error(301_pInt,ext_msg='aTol_crystalliteStress')
 if (any(numerics_integrator <= 0_pInt) .or. any(numerics_integrator >= 6_pInt)) &
                                           call IO_error(301_pInt,ext_msg='integrator')
 if (numerics_unitlength <= 0.0_pReal)     call IO_error(301_pInt,ext_msg='unitlength')
 if (absTol_RGC <= 0.0_pReal)              call IO_error(301_pInt,ext_msg='absTol_RGC')
 if (relTol_RGC <= 0.0_pReal)              call IO_error(301_pInt,ext_msg='relTol_RGC')
 if (absMax_RGC <= 0.0_pReal)              call IO_error(301_pInt,ext_msg='absMax_RGC')
 if (relMax_RGC <= 0.0_pReal)              call IO_error(301_pInt,ext_msg='relMax_RGC')
 if (pPert_RGC <= 0.0_pReal)               call IO_error(301_pInt,ext_msg='pPert_RGC')
 if (xSmoo_RGC <= 0.0_pReal)               call IO_error(301_pInt,ext_msg='xSmoo_RGC')
 if (viscPower_RGC < 0.0_pReal)            call IO_error(301_pInt,ext_msg='viscPower_RGC')
 if (viscModus_RGC < 0.0_pReal)            call IO_error(301_pInt,ext_msg='viscModus_RGC')
 if (refRelaxRate_RGC <= 0.0_pReal)        call IO_error(301_pInt,ext_msg='refRelaxRate_RGC')
 if (maxdRelax_RGC <= 0.0_pReal)           call IO_error(301_pInt,ext_msg='maxdRelax_RGC')
 if (maxVolDiscr_RGC <= 0.0_pReal)         call IO_error(301_pInt,ext_msg='maxVolDiscr_RGC')
 if (volDiscrMod_RGC < 0.0_pReal)          call IO_error(301_pInt,ext_msg='volDiscrMod_RGC')
 if (volDiscrPow_RGC <= 0.0_pReal)         call IO_error(301_pInt,ext_msg='volDiscrPw_RGC')
#ifdef Spectral
 if (err_div_tolRel <= 0.0_pReal)          call IO_error(301_pInt,ext_msg='err_div_tolRel')
 if (err_div_tolAbs <= 0.0_pReal)          call IO_error(301_pInt,ext_msg='err_div_tolAbs')
 if (err_stress_tolrel <= 0.0_pReal)       call IO_error(301_pInt,ext_msg='err_stress_tolrel')
 if (err_stress_tolabs <= 0.0_pReal)       call IO_error(301_pInt,ext_msg='err_stress_tolabs')
 if (itmax <= 1_pInt)                      call IO_error(301_pInt,ext_msg='itmax')
 if (itmin > itmax .or. itmin < 1_pInt)   call IO_error(301_pInt,ext_msg='itmin')
 if (divergence_correction < 0_pInt .or. &
     divergence_correction > 2_pInt)       call IO_error(301_pInt,ext_msg='divergence_correction')
 if (maxCutBack < 0_pInt)                  call IO_error(301_pInt,ext_msg='maxCutBack')
 if (update_gamma .and. &
                   .not. memory_efficient) call IO_error(error_ID = 847_pInt)
#ifdef PETSc
 if (err_curl_tolRel <= 0.0_pReal)         call IO_error(301_pInt,ext_msg='err_curl_tolRel')
 if (err_curl_tolAbs <= 0.0_pReal)         call IO_error(301_pInt,ext_msg='err_curl_tolAbs')
 if (polarAlpha <= 0.0_pReal .or. &
     polarAlpha >  2.0_pReal)              call IO_error(301_pInt,ext_msg='polarAlpha')
 if (polarBeta < 0.0_pReal .or. &
     polarBeta >  2.0_pReal)               call IO_error(301_pInt,ext_msg='polarBeta')
#endif
#endif

end subroutine numerics_init

end module numerics