DAMASK_EICMD/src/numerics.f90

!--------------------------------------------------------------------------------------------------
!> @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
 use IO

#ifdef PETSc
#include <petsc/finclude/petscsys.h>
   use petscsys
#endif
!$ use OMP_LIB

 implicit none
 private

 integer, protected, public :: &
   iJacoStiffness             =  1, &                                                              !< frequency of stiffness update
   nMPstate                   = 10, &                                                              !< materialpoint state loop limit
   randomSeed                 =  0, &                                                              !< fixed seeding for pseudo-random number generator, Default 0: use random seed
   worldrank                  =  0, &                                                              !< MPI worldrank (/=0 for MPI simulations only)
   worldsize                  =  1, &                                                              !< MPI worldsize (/=1 for MPI simulations only)
   numerics_integrator        =  1                                                                 !< method used for state integration Default 1: fix-point iteration
 integer(4), protected, public :: &
   DAMASK_NumThreadsInt       =  0                                                                  !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive
 real(pReal), protected, public :: &
   defgradTolerance           =  1.0e-7_pReal, &                                                    !< deviation of deformation gradient that is still allowed (used by CPFEM to determine outdated ffn1)
   subStepMinHomog            =  1.0e-3_pReal, &                                                    !< minimum (relative) size of sub-step allowed during cutback in homogenization
   subStepSizeHomog           =  0.25_pReal, &                                                      !< size of first substep when cutback in homogenization
   stepIncreaseHomog          =  1.5_pReal, &                                                       !< increase of next substep size when previous substep converged in homogenization
   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
   charLength                 =  1.0_pReal, &                                                       !< characteristic length scale for gradient problems
   residualStiffness          =  1.0e-6_pReal                                                       !< non-zero residual damage   
 logical, protected, public :: &                                                   
   usePingPong                = .true.

!--------------------------------------------------------------------------------------------------
! field parameters:
 real(pReal), protected, public :: &
   err_struct_tolAbs          =  1.0e-10_pReal, &                                                   !< absolute tolerance for mechanical equilibrium
   err_struct_tolRel          =  1.0e-4_pReal, &                                                    !< relative tolerance for mechanical equilibrium
   err_thermal_tolAbs         =  1.0e-2_pReal, &                                                    !< absolute tolerance for thermal equilibrium
   err_thermal_tolRel         =  1.0e-6_pReal, &                                                    !< relative tolerance for thermal equilibrium
   err_damage_tolAbs          =  1.0e-2_pReal, &                                                    !< absolute tolerance for damage evolution
   err_damage_tolRel          =  1.0e-6_pReal                                                       !< relative tolerance for damage evolution
 integer, protected, public :: &
   itmax                      =  250, &                                                             !< maximum number of iterations
   itmin                      =  1, &                                                               !< minimum number of iterations
   stagItMax                  =  10, &                                                              !< max number of field level staggered iterations
   maxCutBack                 =  3                                                                  !< max number of cut backs

!--------------------------------------------------------------------------------------------------
! spectral parameters:
#ifdef Grid
 real(pReal), protected, public :: &
   err_div_tolAbs             =  1.0e-4_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
   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=pStringLen), protected, public :: &
   petsc_options              = ''
#endif

!--------------------------------------------------------------------------------------------------
! FEM parameters:
#ifdef FEM
 integer, protected, public :: &
   integrationOrder           =  2, &                                                              !< order of quadrature rule required
   structOrder                =  2                                                                 !< order of displacement shape functions
 logical, protected, public :: & 
   BBarStabilisation          = .false.                                                  
 character(len=*), parameter, public :: &
   petsc_defaultOptions    = '-mech_snes_type newtonls &
                             &-mech_snes_linesearch_type cp &
                             &-mech_snes_ksp_ew &
                             &-mech_snes_ksp_ew_rtol0 0.01 &
                             &-mech_snes_ksp_ew_rtolmax 0.01 &
                             &-mech_ksp_type fgmres &
                             &-mech_ksp_max_it 25 &
                             &-mech_pc_type ml &
                             &-mech_mg_levels_ksp_type chebyshev &
                             &-mech_mg_levels_pc_type sor &
                             &-mech_pc_ml_nullspace user'
 character(len=pStringLen), protected, public :: &
   petsc_options           = ''
#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
!$ integer ::                                gotDAMASK_NUM_THREADS = 1
 integer :: i,j, ierr
 integer, allocatable, dimension(:) :: chunkPos
 character(len=pStringLen), dimension(:), allocatable :: fileContent
 character(len=pStringLen) :: &
   tag ,&
   line
 logical :: fexist
!$ character(len=6) DAMASK_NumThreadsString                                                         ! environment variable DAMASK_NUM_THREADS

#ifdef PETSc
 call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,ierr);CHKERRQ(ierr)
 call MPI_Comm_size(PETSC_COMM_WORLD,worldsize,ierr);CHKERRQ(ierr)
#endif
 write(6,'(/,a)') ' <<<+-  numerics init  -+>>>'
 
!$ 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) then                                                              ! could not get number of threads, set it to 1
!$   call IO_warning(35,ext_msg='BEGIN:'//DAMASK_NumThreadsString//':END')
!$   DAMASK_NumThreadsInt = 1_4
!$ else
!$   read(DAMASK_NumThreadsString,'(i6)') DAMASK_NumThreadsInt                                      ! read as integer
!$   if (DAMASK_NumThreadsInt < 1_4) DAMASK_NumThreadsInt = 1_4                                     ! in case of string conversion fails, set it to one
!$ endif
!$ call omp_set_num_threads(DAMASK_NumThreadsInt)                                                   ! set number of threads for parallel execution
 
 inquire(file='numerics.config', exist=fexist)
 
 fileExists: if (fexist) then
   write(6,'(a,/)') ' using values from config file'
   flush(6)
   fileContent = IO_read_ASCII('numerics.config') 
   do j=1, size(fileContent)
   
!--------------------------------------------------------------------------------------------------
! read variables from config file and overwrite default parameters if keyword is present
     line = fileContent(j)
     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
     chunkPos = IO_stringPos(line)
     tag = IO_lc(IO_stringValue(line,chunkPos,1))                                                   ! extract key

     select case(tag)
       case ('defgradtolerance')
         defgradTolerance = IO_floatValue(line,chunkPos,2)
       case ('ijacostiffness')
         iJacoStiffness = IO_intValue(line,chunkPos,2)
       case ('nmpstate')
         nMPstate = IO_intValue(line,chunkPos,2)
       case ('substepminhomog')
         subStepMinHomog = IO_floatValue(line,chunkPos,2)
       case ('substepsizehomog')
         subStepSizeHomog = IO_floatValue(line,chunkPos,2)
       case ('stepincreasehomog')
         stepIncreaseHomog = IO_floatValue(line,chunkPos,2)
       case ('integrator')
         numerics_integrator = IO_intValue(line,chunkPos,2)
       case ('usepingpong')
         usepingpong = IO_intValue(line,chunkPos,2) > 0
       case ('unitlength')
         numerics_unitlength = IO_floatValue(line,chunkPos,2)

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

!--------------------------------------------------------------------------------------------------
! random seeding parameter
       case ('random_seed','fixed_seed')
         randomSeed = IO_intValue(line,chunkPos,2)

!--------------------------------------------------------------------------------------------------
! gradient parameter
       case ('charlength')
         charLength = IO_floatValue(line,chunkPos,2)
       case ('residualstiffness')
         residualStiffness = IO_floatValue(line,chunkPos,2)

!--------------------------------------------------------------------------------------------------
! field parameters
       case ('err_struct_tolabs')
         err_struct_tolAbs = IO_floatValue(line,chunkPos,2)
       case ('err_struct_tolrel')
         err_struct_tolRel = IO_floatValue(line,chunkPos,2)
       case ('err_thermal_tolabs')
         err_thermal_tolabs = IO_floatValue(line,chunkPos,2)
       case ('err_thermal_tolrel')
         err_thermal_tolrel = IO_floatValue(line,chunkPos,2)
       case ('err_damage_tolabs')
         err_damage_tolabs = IO_floatValue(line,chunkPos,2)
       case ('err_damage_tolrel')
         err_damage_tolrel = IO_floatValue(line,chunkPos,2)
       case ('itmax')
         itmax = IO_intValue(line,chunkPos,2)
       case ('itmin')
         itmin = IO_intValue(line,chunkPos,2)
       case ('maxcutback')
         maxCutBack = IO_intValue(line,chunkPos,2)
       case ('maxstaggerediter')
         stagItMax = IO_intValue(line,chunkPos,2)

!--------------------------------------------------------------------------------------------------
! spectral parameters
#ifdef Grid
       case ('err_div_tolabs')
         err_div_tolAbs = IO_floatValue(line,chunkPos,2)
       case ('err_div_tolrel')
         err_div_tolRel = IO_floatValue(line,chunkPos,2)
       case ('err_stress_tolrel')
         err_stress_tolrel = IO_floatValue(line,chunkPos,2)
       case ('err_stress_tolabs')
         err_stress_tolabs = IO_floatValue(line,chunkPos,2)
       case ('petsc_options')
         petsc_options = trim(line(chunkPos(4):))
       case ('err_curl_tolabs')
         err_curl_tolAbs = IO_floatValue(line,chunkPos,2)
       case ('err_curl_tolrel')
         err_curl_tolRel = IO_floatValue(line,chunkPos,2)
       case ('polaralpha')
         polarAlpha = IO_floatValue(line,chunkPos,2)
       case ('polarbeta')
         polarBeta = IO_floatValue(line,chunkPos,2)
#endif

!--------------------------------------------------------------------------------------------------
! FEM parameters
#ifdef FEM
       case ('integrationorder')
         integrationorder = IO_intValue(line,chunkPos,2)
       case ('structorder')
         structorder = IO_intValue(line,chunkPos,2)
       case ('petsc_options')
         petsc_options = trim(line(chunkPos(4):))
       case ('bbarstabilisation')
         BBarStabilisation = IO_intValue(line,chunkPos,2) > 0
#endif
     end select
   enddo


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


!--------------------------------------------------------------------------------------------------
! writing parameters to output
 write(6,'(a24,1x,es8.1)')  ' defgradTolerance:       ',defgradTolerance
 write(6,'(a24,1x,i8)')     ' iJacoStiffness:         ',iJacoStiffness
 write(6,'(a24,1x,i8)')     ' integrator:             ',numerics_integrator
 write(6,'(a24,1x,L8)')     ' use ping pong scheme:   ',usepingpong
 write(6,'(a24,1x,es8.1,/)')' unitlength:             ',numerics_unitlength

 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,'(a16,1x,i16,/)')    ' random_seed:    ',randomSeed
 if (randomSeed <= 0) &
   write(6,'(a,/)')           ' random seed will be generated!'

!--------------------------------------------------------------------------------------------------
! gradient parameter
 write(6,'(a24,1x,es8.1)')   ' charLength:             ',charLength
 write(6,'(a24,1x,es8.1)')   ' residualStiffness:      ',residualStiffness

!--------------------------------------------------------------------------------------------------
! openMP parameter
 !$  write(6,'(a24,1x,i8,/)')   ' number of threads:      ',DAMASK_NumThreadsInt

!--------------------------------------------------------------------------------------------------
! field parameters
 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)')      ' maxStaggeredIter:       ',stagItMax
 write(6,'(a24,1x,es8.1)')   ' err_struct_tolAbs:      ',err_struct_tolAbs
 write(6,'(a24,1x,es8.1)')   ' err_struct_tolRel:      ',err_struct_tolRel
 write(6,'(a24,1x,es8.1)')   ' err_thermal_tolabs:     ',err_thermal_tolabs
 write(6,'(a24,1x,es8.1)')   ' err_thermal_tolrel:     ',err_thermal_tolrel
 write(6,'(a24,1x,es8.1)')   ' err_damage_tolabs:      ',err_damage_tolabs
 write(6,'(a24,1x,es8.1)')   ' err_damage_tolrel:      ',err_damage_tolrel

!--------------------------------------------------------------------------------------------------
! spectral parameters
#ifdef Grid
 write(6,'(a24,1x,es8.1)')   ' err_stress_tolAbs:      ',err_stress_tolAbs
 write(6,'(a24,1x,es8.1)')   ' err_stress_tolRel:      ',err_stress_tolRel
 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_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)')       ' PETSc_options:          ',trim(petsc_options)
#endif

!--------------------------------------------------------------------------------------------------
! spectral parameters
#ifdef FEM
 write(6,'(a24,1x,i8)')      ' integrationOrder:       ',integrationOrder
 write(6,'(a24,1x,i8)')      ' structOrder:            ',structOrder
 write(6,'(a24,1x,a)')       ' PETSc_options:          ',trim(petsc_defaultOptions)//' '//trim(petsc_options)
 write(6,'(a24,1x,L8)')      ' B-Bar stabilisation:    ',BBarStabilisation
#endif
 

!--------------------------------------------------------------------------------------------------
! sanity checks
 if (defgradTolerance <= 0.0_pReal)        call IO_error(301,ext_msg='defgradTolerance')
 if (iJacoStiffness < 1)                   call IO_error(301,ext_msg='iJacoStiffness')
 if (nMPstate < 1)                         call IO_error(301,ext_msg='nMPstate')
 if (subStepMinHomog <= 0.0_pReal)         call IO_error(301,ext_msg='subStepMinHomog')
 if (subStepSizeHomog <= 0.0_pReal)        call IO_error(301,ext_msg='subStepSizeHomog')
 if (stepIncreaseHomog <= 0.0_pReal)       call IO_error(301,ext_msg='stepIncreaseHomog')
 if (numerics_integrator <= 0 .or. numerics_integrator >= 6) &
                                           call IO_error(301,ext_msg='integrator')
 if (numerics_unitlength <= 0.0_pReal)     call IO_error(301,ext_msg='unitlength')
 if (absTol_RGC <= 0.0_pReal)              call IO_error(301,ext_msg='absTol_RGC')
 if (relTol_RGC <= 0.0_pReal)              call IO_error(301,ext_msg='relTol_RGC')
 if (absMax_RGC <= 0.0_pReal)              call IO_error(301,ext_msg='absMax_RGC')
 if (relMax_RGC <= 0.0_pReal)              call IO_error(301,ext_msg='relMax_RGC')
 if (pPert_RGC <= 0.0_pReal)               call IO_error(301,ext_msg='pPert_RGC')
 if (xSmoo_RGC <= 0.0_pReal)               call IO_error(301,ext_msg='xSmoo_RGC')
 if (viscPower_RGC < 0.0_pReal)            call IO_error(301,ext_msg='viscPower_RGC')
 if (viscModus_RGC < 0.0_pReal)            call IO_error(301,ext_msg='viscModus_RGC')
 if (refRelaxRate_RGC <= 0.0_pReal)        call IO_error(301,ext_msg='refRelaxRate_RGC')
 if (maxdRelax_RGC <= 0.0_pReal)           call IO_error(301,ext_msg='maxdRelax_RGC')
 if (maxVolDiscr_RGC <= 0.0_pReal)         call IO_error(301,ext_msg='maxVolDiscr_RGC')
 if (volDiscrMod_RGC < 0.0_pReal)          call IO_error(301,ext_msg='volDiscrMod_RGC')
 if (volDiscrPow_RGC <= 0.0_pReal)         call IO_error(301,ext_msg='volDiscrPw_RGC')
 if (residualStiffness < 0.0_pReal)        call IO_error(301,ext_msg='residualStiffness')
 if (itmax <= 1)                           call IO_error(301,ext_msg='itmax')
 if (itmin > itmax .or. itmin < 1)         call IO_error(301,ext_msg='itmin')
 if (maxCutBack < 0)                       call IO_error(301,ext_msg='maxCutBack')
 if (stagItMax < 0)                        call IO_error(301,ext_msg='maxStaggeredIter')
 if (err_struct_tolRel <= 0.0_pReal)       call IO_error(301,ext_msg='err_struct_tolRel')
 if (err_struct_tolAbs <= 0.0_pReal)       call IO_error(301,ext_msg='err_struct_tolAbs')
 if (err_thermal_tolabs <= 0.0_pReal)      call IO_error(301,ext_msg='err_thermal_tolabs')
 if (err_thermal_tolrel <= 0.0_pReal)      call IO_error(301,ext_msg='err_thermal_tolrel')
 if (err_damage_tolabs <= 0.0_pReal)       call IO_error(301,ext_msg='err_damage_tolabs')
 if (err_damage_tolrel <= 0.0_pReal)       call IO_error(301,ext_msg='err_damage_tolrel')
#ifdef Grid
 if (err_stress_tolrel <= 0.0_pReal)       call IO_error(301,ext_msg='err_stress_tolRel')
 if (err_stress_tolabs <= 0.0_pReal)       call IO_error(301,ext_msg='err_stress_tolAbs')
 if (err_div_tolRel < 0.0_pReal)           call IO_error(301,ext_msg='err_div_tolRel')
 if (err_div_tolAbs <= 0.0_pReal)          call IO_error(301,ext_msg='err_div_tolAbs')
 if (err_curl_tolRel < 0.0_pReal)          call IO_error(301,ext_msg='err_curl_tolRel')
 if (err_curl_tolAbs <= 0.0_pReal)         call IO_error(301,ext_msg='err_curl_tolAbs')
 if (polarAlpha <= 0.0_pReal .or. &
     polarAlpha >  2.0_pReal)              call IO_error(301,ext_msg='polarAlpha')
 if (polarBeta < 0.0_pReal .or. &
     polarBeta > 2.0_pReal)                call IO_error(301,ext_msg='polarBeta')
#endif

end subroutine numerics_init

end module numerics