223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

keep connected data together 

- avoids dependencies
- easier to read and modify
This commit is contained in:
Martin Diehl 2019-04-11 07:24:04 +02:00
parent 8d90cfb600
commit 67eb39255a
2 changed files with 49 additions and 56 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

73
src/crystallite.f90
View File

@ -111,6 +111,20 @@ module crystallite
end type tOutput end type tOutput
type(tOutput), allocatable, dimension(:), private :: output_constituent type(tOutput), allocatable, dimension(:), private :: output_constituent
type, private :: tNumerics
real(pReal) :: &
subStepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback
subStepSizeCryst, & !< size of first substep when cutback
subStepSizeLp, & !< size of first substep when cutback in Lp calculation
subStepSizeLi, & !< size of first substep when cutback in Li calculation
stepIncreaseCryst, & !< increase of next substep size when previous substep converged
rTol_crystalliteState, & !< relative tolerance in state loop
rTol_crystalliteStress, & !< relative tolerance in stress loop
aTol_crystalliteStress !< absolute tolerance in stress loop
end type tNumerics
type(tNumerics) :: num ! numerics parameters. Better name?
procedure(), pointer :: integrateState procedure(), pointer :: integrateState
public :: & public :: &
@ -165,6 +179,7 @@ subroutine crystallite_init
use config, only: & use config, only: &
config_deallocate, & config_deallocate, &
config_crystallite, & config_crystallite, &
config_numerics, &
config_phase, & config_phase, &
crystallite_name crystallite_name
use constitutive, only: & use constitutive, only: &
@ -243,6 +258,24 @@ subroutine crystallite_init
allocate(crystallite_sizePostResult(maxval(crystallite_Noutput), & allocate(crystallite_sizePostResult(maxval(crystallite_Noutput), &
size(config_crystallite)), source=0) size(config_crystallite)), source=0)
num%subStepMinCryst = config_numerics%getFloat('subStepMinCryst', defaultVal=1.0e-3_pReal)
num%subStepSizeCryst = config_numerics%getFloat('subStepSizeCryst', defaultVal=0.25_pReal)
num%subStepSizeLp = config_numerics%getFloat('subStepSizeLp', defaultVal=0.5_pReal)
num%subStepSizeLi = config_numerics%getFloat('subStepSizeLi', defaultVal=0.5_pReal)
num%stepIncreaseCryst = config_numerics%getFloat('stepIncreaseCryst', defaultVal=1.5_pReal)
num%rTol_crystalliteState = config_numerics%getFloat('rTol_crystalliteState', defaultVal=1.0e-6_pReal)
num%rTol_crystalliteStress = config_numerics%getFloat('rTol_crystalliteStress',defaultVal=1.0e-6_pReal)
num%aTol_crystalliteStress = config_numerics%getFloat('aTol_crystalliteStress',defaultVal=1.0e-8_pReal)
if(num%subStepMinCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepMinCryst')
if(num%subStepSizeCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeCryst')
if(num%stepIncreaseCryst <= 0.0_pReal) call IO_error(301,ext_msg='stepIncreaseCryst')
if(num%subStepSizeLp <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeLp')
if(num%subStepSizeLi <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeLi')
if(num%rTol_crystalliteState <= 0.0_pReal) call IO_error(301,ext_msg='rTol_crystalliteState')
if(num%rTol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='rTol_crystalliteStress')
if(num%aTol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='aTol_crystalliteStress')
select case(numerics_integrator) select case(numerics_integrator)
case(1) case(1)
integrateState => integrateStateFPI integrateState => integrateStateFPI
@ -433,10 +466,6 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
use prec, only: & use prec, only: &
tol_math_check, & tol_math_check, &
dNeq0 dNeq0
use numerics, only: &
subStepMinCryst, &
subStepSizeCryst, &
stepIncreaseCryst
#ifdef DEBUG #ifdef DEBUG
use debug, only: & use debug, only: &
debug_level, & debug_level, &
@ -519,7 +548,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partionedF0(1:3,1:3,c,i,e) crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partionedF0(1:3,1:3,c,i,e)
crystallite_subS0(1:3,1:3,c,i,e) = crystallite_partionedS0(1:3,1:3,c,i,e) crystallite_subS0(1:3,1:3,c,i,e) = crystallite_partionedS0(1:3,1:3,c,i,e)
crystallite_subFrac(c,i,e) = 0.0_pReal crystallite_subFrac(c,i,e) = 0.0_pReal
crystallite_subStep(c,i,e) = 1.0_pReal/subStepSizeCryst crystallite_subStep(c,i,e) = 1.0_pReal/num%subStepSizeCryst
crystallite_todo(c,i,e) = .true. crystallite_todo(c,i,e) = .true.
crystallite_converged(c,i,e) = .false. ! pretend failed step of 1/subStepSizeCryst crystallite_converged(c,i,e) = .false. ! pretend failed step of 1/subStepSizeCryst
endif homogenizationRequestsCalculation endif homogenizationRequestsCalculation
@ -554,7 +583,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
formerSubStep = crystallite_subStep(c,i,e) formerSubStep = crystallite_subStep(c,i,e)
crystallite_subFrac(c,i,e) = crystallite_subFrac(c,i,e) + crystallite_subStep(c,i,e) crystallite_subFrac(c,i,e) = crystallite_subFrac(c,i,e) + crystallite_subStep(c,i,e)
crystallite_subStep(c,i,e) = min(1.0_pReal - crystallite_subFrac(c,i,e), & crystallite_subStep(c,i,e) = min(1.0_pReal - crystallite_subFrac(c,i,e), &
stepIncreaseCryst * crystallite_subStep(c,i,e)) num%stepIncreaseCryst * crystallite_subStep(c,i,e))
crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on? crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on?
if (crystallite_todo(c,i,e)) then if (crystallite_todo(c,i,e)) then
@ -584,7 +613,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! cut back (reduced time and restore) ! cut back (reduced time and restore)
else else
crystallite_subStep(c,i,e) = subStepSizeCryst * crystallite_subStep(c,i,e) crystallite_subStep(c,i,e) = num%subStepSizeCryst * crystallite_subStep(c,i,e)
crystallite_Fp (1:3,1:3,c,i,e) = crystallite_subFp0(1:3,1:3,c,i,e) crystallite_Fp (1:3,1:3,c,i,e) = crystallite_subFp0(1:3,1:3,c,i,e)
crystallite_invFp(1:3,1:3,c,i,e) = math_inv33(crystallite_Fp (1:3,1:3,c,i,e)) crystallite_invFp(1:3,1:3,c,i,e) = math_inv33(crystallite_Fp (1:3,1:3,c,i,e))
crystallite_Fi (1:3,1:3,c,i,e) = crystallite_subFi0(1:3,1:3,c,i,e) crystallite_Fi (1:3,1:3,c,i,e) = crystallite_subFi0(1:3,1:3,c,i,e)
@ -602,7 +631,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
enddo enddo
! cant restore dotState here, since not yet calculated in first cutback after initialization ! cant restore dotState here, since not yet calculated in first cutback after initialization
crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > subStepMinCryst ! still on track or already done (beyond repair) crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > num%subStepMinCryst ! still on track or already done (beyond repair)
#ifdef DEBUG #ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 &
.and. ((e == debug_e .and. i == debug_i .and. c == debug_g) & .and. ((e == debug_e .and. i == debug_i .and. c == debug_g) &
@ -652,7 +681,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! integrate --- requires fully defined state array (basic + dependent state) ! integrate --- requires fully defined state array (basic + dependent state)
if (any(crystallite_todo)) call integrateState ! TODO: unroll into proper elementloop to avoid N^2 for single point evaluation if (any(crystallite_todo)) call integrateState ! TODO: unroll into proper elementloop to avoid N^2 for single point evaluation
where(.not. crystallite_converged .and. crystallite_subStep > subStepMinCryst) & ! do not try non-converged but fully cutbacked any further where(.not. crystallite_converged .and. crystallite_subStep > num%subStepMinCryst) & ! do not try non-converged but fully cutbacked any further
crystallite_todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation crystallite_todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation
@ -1238,11 +1267,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
use prec, only: tol_math_check, & use prec, only: tol_math_check, &
dEq0 dEq0
use numerics, only: nStress, & use numerics, only: nStress, &
aTol_crystalliteStress, & iJacoLpresiduum
rTol_crystalliteStress, &
iJacoLpresiduum, &
subStepSizeLp, &
subStepSizeLi
#ifdef DEBUG #ifdef DEBUG
use debug, only: debug_level, & use debug, only: debug_level, &
debug_e, & debug_e, &
@ -1442,8 +1467,8 @@ logical function integrateStress(ipc,ip,el,timeFraction)
#endif #endif
!* update current residuum and check for convergence of loop !* update current residuum and check for convergence of loop
aTolLp = max(rTol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error aTolLp = max(num%rTol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error
aTol_crystalliteStress) ! minimum lower cutoff num%aTol_crystalliteStress) ! minimum lower cutoff
residuumLp = Lpguess - Lp_constitutive residuumLp = Lpguess - Lp_constitutive
if (any(IEEE_is_NaN(residuumLp))) then if (any(IEEE_is_NaN(residuumLp))) then
@ -1463,7 +1488,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
Lpguess_old = Lpguess Lpguess_old = Lpguess
steplengthLp = 1.0_pReal ! ...proceed with normal step length (calculate new search direction) steplengthLp = 1.0_pReal ! ...proceed with normal step length (calculate new search direction)
else ! not converged and residuum not improved... else ! not converged and residuum not improved...
steplengthLp = subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction steplengthLp = num%subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction
Lpguess = Lpguess_old + steplengthLp * deltaLp Lpguess = Lpguess_old + steplengthLp * deltaLp
#ifdef DEBUG #ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 &
@ -1541,8 +1566,8 @@ logical function integrateStress(ipc,ip,el,timeFraction)
#endif #endif
!* update current residuum and check for convergence of loop !* update current residuum and check for convergence of loop
aTolLi = max(rTol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error aTolLi = max(num%rTol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error
aTol_crystalliteStress) ! minimum lower cutoff num%aTol_crystalliteStress) ! minimum lower cutoff
residuumLi = Liguess - Li_constitutive residuumLi = Liguess - Li_constitutive
if (any(IEEE_is_NaN(residuumLi))) then ! NaN in residuum... if (any(IEEE_is_NaN(residuumLi))) then ! NaN in residuum...
#ifdef DEBUG #ifdef DEBUG
@ -1561,7 +1586,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
Liguess_old = Liguess Liguess_old = Liguess
steplengthLi = 1.0_pReal ! ...proceed with normal step length (calculate new search direction) steplengthLi = 1.0_pReal ! ...proceed with normal step length (calculate new search direction)
else ! not converged and residuum not improved... else ! not converged and residuum not improved...
steplengthLi = subStepSizeLi * steplengthLi ! ...try with smaller step length in same direction steplengthLi = num%subStepSizeLi * steplengthLi ! ...try with smaller step length in same direction
Liguess = Liguess_old + steplengthLi * deltaLi Liguess = Liguess_old + steplengthLi * deltaLi
cycle LiLoop cycle LiLoop
endif endif
@ -2295,14 +2320,14 @@ end subroutine setConvergenceFlag
!> @brief determines whether a point is converged !> @brief determines whether a point is converged
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
logical pure function converged(residuum,state,aTol) logical pure function converged(residuum,state,aTol)
use prec, only: &
dEq0
use numerics, only: &
rTol => rTol_crystalliteState
implicit none implicit none
real(pReal), intent(in), dimension(:) ::& real(pReal), intent(in), dimension(:) ::&
residuum, state, aTol residuum, state, aTol
real(pReal) :: &
rTol
rTol = num%rTol_crystalliteState
converged = all(abs(residuum) <= max(aTol, rTol*abs(state))) converged = all(abs(residuum) <= max(aTol, rTol*abs(state)))

32
src/numerics.f90
View File

@ -26,17 +26,9 @@ module numerics
DAMASK_NumThreadsInt = 0 !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive DAMASK_NumThreadsInt = 0 !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive
real(pReal), protected, public :: & real(pReal), protected, public :: &
defgradTolerance = 1.0e-7_pReal, & !< deviation of deformation gradient that is still allowed (used by CPFEM to determine outdated ffn1) defgradTolerance = 1.0e-7_pReal, & !< deviation of deformation gradient that is still allowed (used by CPFEM to determine outdated ffn1)
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 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 subStepSizeHomog = 0.25_pReal, & !< size of first substep when cutback in homogenization
subStepSizeLp = 0.5_pReal, & !< size of first substep when cutback in Lp calculation
subStepSizeLi = 0.5_pReal, & !< size of first substep when cutback in Li calculation
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 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_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 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 absTol_RGC = 1.0e+4_pReal, & !< absolute tolerance of RGC residuum
relTol_RGC = 1.0e-3_pReal, & !< relative tolerance of RGC residuum relTol_RGC = 1.0e-3_pReal, & !< relative tolerance of RGC residuum
@ -205,15 +197,10 @@ subroutine numerics_init
case ('nstress') case ('nstress')
nStress = IO_intValue(line,chunkPos,2_pInt) nStress = IO_intValue(line,chunkPos,2_pInt)
case ('substepmincryst') case ('substepmincryst')
subStepMinCryst = IO_floatValue(line,chunkPos,2_pInt)
case ('substepsizecryst') case ('substepsizecryst')
subStepSizeCryst = IO_floatValue(line,chunkPos,2_pInt)
case ('stepincreasecryst') case ('stepincreasecryst')
stepIncreaseCryst = IO_floatValue(line,chunkPos,2_pInt)
case ('substepsizelp') case ('substepsizelp')
subStepSizeLp = IO_floatValue(line,chunkPos,2_pInt)
case ('substepsizeli') case ('substepsizeli')
subStepSizeLi = IO_floatValue(line,chunkPos,2_pInt)
case ('substepminhomog') case ('substepminhomog')
subStepMinHomog = IO_floatValue(line,chunkPos,2_pInt) subStepMinHomog = IO_floatValue(line,chunkPos,2_pInt)
case ('substepsizehomog') case ('substepsizehomog')
@ -221,11 +208,8 @@ subroutine numerics_init
case ('stepincreasehomog') case ('stepincreasehomog')
stepIncreaseHomog = IO_floatValue(line,chunkPos,2_pInt) stepIncreaseHomog = IO_floatValue(line,chunkPos,2_pInt)
case ('rtol_crystallitestate') case ('rtol_crystallitestate')
rTol_crystalliteState = IO_floatValue(line,chunkPos,2_pInt)
case ('rtol_crystallitestress') case ('rtol_crystallitestress')
rTol_crystalliteStress = IO_floatValue(line,chunkPos,2_pInt)
case ('atol_crystallitestress') case ('atol_crystallitestress')
aTol_crystalliteStress = IO_floatValue(line,chunkPos,2_pInt)
case ('integrator') case ('integrator')
numerics_integrator = IO_intValue(line,chunkPos,2_pInt) numerics_integrator = IO_intValue(line,chunkPos,2_pInt)
case ('usepingpong') case ('usepingpong')
@ -350,16 +334,8 @@ subroutine numerics_init
write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness
write(6,'(a24,1x,i8)') ' iJacoLpresiduum: ',iJacoLpresiduum write(6,'(a24,1x,i8)') ' iJacoLpresiduum: ',iJacoLpresiduum
write(6,'(a24,1x,i8)') ' nCryst: ',nCryst 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,es8.1)') ' subStepSizeLp: ',subStepSizeLp
write(6,'(a24,1x,es8.1)') ' subStepSizeLi: ',subStepSizeLi
write(6,'(a24,1x,i8)') ' nState: ',nState write(6,'(a24,1x,i8)') ' nState: ',nState
write(6,'(a24,1x,i8)') ' nStress: ',nStress write(6,'(a24,1x,i8)') ' nStress: ',nStress
write(6,'(a24,1x,es8.1)') ' rTol_crystalliteState: ',rTol_crystalliteState
write(6,'(a24,1x,es8.1)') ' rTol_crystalliteStress: ',rTol_crystalliteStress
write(6,'(a24,1x,es8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress
write(6,'(a24,1x,i8)') ' integrator: ',numerics_integrator write(6,'(a24,1x,i8)') ' integrator: ',numerics_integrator
write(6,'(a24,1x,L8)') ' use ping pong scheme: ',usepingpong 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,/)')' unitlength: ',numerics_unitlength
@ -446,17 +422,9 @@ subroutine numerics_init
if (nCryst < 1_pInt) call IO_error(301_pInt,ext_msg='nCryst') 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 (nState < 1_pInt) call IO_error(301_pInt,ext_msg='nState')
if (nStress < 1_pInt) call IO_error(301_pInt,ext_msg='nStress') 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 (subStepSizeLp <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeLp')
if (subStepSizeLi <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeLi')
if (subStepMinHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepMinHomog') 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 (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 (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_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 (numerics_integrator <= 0_pInt .or. numerics_integrator >= 6_pInt) & if (numerics_integrator <= 0_pInt .or. numerics_integrator >= 6_pInt) &
call IO_error(301_pInt,ext_msg='integrator') call IO_error(301_pInt,ext_msg='integrator')
if (numerics_unitlength <= 0.0_pReal) call IO_error(301_pInt,ext_msg='unitlength') if (numerics_unitlength <= 0.0_pReal) call IO_error(301_pInt,ext_msg='unitlength')