Merge branch 'development' into fix-Intel2021

PRIVATE

@@ -1 +1 @@
-Subproject commit d10a94a36056d99ea0892d1c6c69aa0e1dbe9d46
+Subproject commit c0716e46b72439a49212c5e9803a5886c03b4739

VERSION

@@ -1 +1 @@
-v3.0.0-alpha3-101-g0f4e7e96d
+v3.0.0-alpha3-129-g7440cc32a

examples/config/phase/damage/anisobrittle_cubic.yaml

@@ -1,8 +1,11 @@
-N_cl: [3]
-dot_o: 1e-3
-g_crit: [0.50e7]
-q: 20
-s_crit: [0.006666]
 type: anisobrittle
+N_cl: [3]
+g_crit: [0.50e7]
+s_crit: [0.006666]
+dot_o: 1e-3
+q: 20
+
+output: [f_phi]
+
 D_11: 1.0
 M: 0.001

examples/config/phase/damage/isobrittle_generic.yaml

@@ -0,0 +1,9 @@
+type: isobrittle
+W_crit: 1400000.0
+m: 1.0
+isoBrittle_atol: 0.01
+
+output: [f_phi]
+
+D_11: 1.0
+M: 0.001

examples/mesh/material.yaml

@@ -1,15 +1,15 @@
 homogenization:
   SX:
     N_constituents: 1
-    mechanics: {type: pass}
+    mechanical: {type: pass}
 
 phase:
   Aluminum:
     lattice: cF
-    mechanics:
+    mechanical:
       output: [F, P, F_e, F_p, L_p]
-      elasticity: {type: Hooke, C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9}
+      elastic: {type: Hooke, C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9}
-      plasticity:
+      plastic:
         type: phenopowerlaw
         N_sl: [12]
         a_sl: 2.25

python/damask/_rotation.py

@@ -629,7 +629,7 @@ class Rotation:
         else:
             if np.any(qu[...,0] < 0.0):
                 raise ValueError('Quaternion with negative first (real) component.')
-        if not np.all(np.isclose(np.linalg.norm(qu,axis=-1), 1.0)):
+        if not np.all(np.isclose(np.linalg.norm(qu,axis=-1), 1.0,rtol=0.0)):
             raise ValueError('Quaternion is not of unit length.')
 
         return Rotation(qu)

python/damask/_table.py

@@ -384,7 +384,7 @@ class Table:
 
         Returns
         -------
-        udated : damask.Table
+        updated : damask.Table
             Updated table.
 
         """
@@ -412,7 +412,7 @@ class Table:
 
         Returns
         -------
-        udated : damask.Table
+        updated : damask.Table
             Updated table.
 
         """
@@ -435,7 +435,7 @@ class Table:
 
         Returns
         -------
-        udated : damask.Table
+        updated : damask.Table
             Updated table.
 
         """
@@ -461,7 +461,7 @@ class Table:
 
         Returns
         -------
-        udated : damask.Table
+        updated : damask.Table
             Updated table.
 
         """
@@ -494,7 +494,7 @@ class Table:
 
         Returns
         -------
-        udated : damask.Table
+        updated : damask.Table
             Updated table.
 
         """
@@ -519,7 +519,7 @@ class Table:
 
         Returns
         -------
-        udated : damask.Table
+        updated : damask.Table
             Updated table.
 
         """

src/IO.f90

@@ -113,8 +113,7 @@ end function IO_readlines
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Read whole file.
-!> @details ensures that the string ends with a new line (expected UNIX behavior) and rejects
+!> @details ensures that the string ends with a new line (expected UNIX behavior)
-! windows (CRLF) line endings
 !--------------------------------------------------------------------------------------------------
 function IO_read(fileName) result(fileContent)
 
@@ -124,8 +123,7 @@ function IO_read(fileName) result(fileContent)
   integer ::  &
     fileLength, &
     fileUnit, &
-    myStat, &
+    myStat
-    firstEOL
   character, parameter :: CR = achar(13)
 
 
@@ -143,12 +141,17 @@ function IO_read(fileName) result(fileContent)
   if(myStat /= 0) call IO_error(102,ext_msg=trim(fileName))
   close(fileUnit)
 
+  foundCRLF: if (scan(fileContent(:index(fileContent,IO_EOL)),CR) /= 0) then
+    CRLF2LF: block
+      integer :: c
+      do c=1, len(fileContent)
+        if (fileContent(c:c) == CR) fileContent(c:c) = ' '
+      enddo
+    end block CRLF2LF
+  endif foundCRLF
 
   if(fileContent(fileLength:fileLength) /= IO_EOL) fileContent = fileContent//IO_EOL                ! ensure EOL@EOF
 
-  firstEOL = index(fileContent,IO_EOL)
-  if(scan(fileContent(firstEOL:firstEOL),CR) /= 0) call IO_error(115)
-
 end function IO_read
 
 
@@ -400,9 +403,6 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
       msg = 'invalid character for logical:'
     case (114)
       msg = 'cannot decode base64 string:'
-    case (115)
-      msg = 'found CR. Windows file endings (CRLF) are not supported.'
-
 
 !--------------------------------------------------------------------------------------------------
 ! lattice error messages
@@ -644,11 +644,16 @@ subroutine selfTest
   if(dNeq(1.0_pReal, IO_stringAsFloat('1.0')))      error stop 'IO_stringAsFloat'
   if(dNeq(1.0_pReal, IO_stringAsFloat('1e0')))      error stop 'IO_stringAsFloat'
   if(dNeq(0.1_pReal, IO_stringAsFloat('1e-1')))     error stop 'IO_stringAsFloat'
+  if(dNeq(0.1_pReal, IO_stringAsFloat('1.0e-1')))   error stop 'IO_stringAsFloat'
+  if(dNeq(0.1_pReal, IO_stringAsFloat('1.00e-1')))  error stop 'IO_stringAsFloat'
+  if(dNeq(10._pReal, IO_stringAsFloat(' 1.0e+1 '))) error stop 'IO_stringAsFloat'
 
   if(3112019  /= IO_stringAsInt( '3112019'))        error stop 'IO_stringAsInt'
   if(3112019  /= IO_stringAsInt(' 3112019'))        error stop 'IO_stringAsInt'
   if(-3112019 /= IO_stringAsInt('-3112019'))        error stop 'IO_stringAsInt'
   if(3112019  /= IO_stringAsInt('+3112019 '))       error stop 'IO_stringAsInt'
+  if(3112019  /= IO_stringAsInt('03112019 '))       error stop 'IO_stringAsInt'
+  if(3112019  /= IO_stringAsInt('+03112019'))       error stop 'IO_stringAsInt'
 
   if(.not. IO_stringAsBool(' true'))                error stop 'IO_stringAsBool'
   if(.not. IO_stringAsBool(' True '))               error stop 'IO_stringAsBool'

src/phase.f90

@@ -67,10 +67,6 @@ module phase
     damageState
 
 
-  integer, public, protected :: &
-    phase_plasticity_maxSizeDotState, &
-    phase_source_maxSizeDotState
-
   interface
 
 ! == cleaned:begin =================================================================================
@@ -235,11 +231,16 @@ module phase
       logical :: converged_
     end function crystallite_stress
 
+    !ToDo: Try to merge the all stiffness functions
     module function phase_homogenizedC(ph,en) result(C)
       integer, intent(in) :: ph, en
       real(pReal), dimension(6,6) :: C
     end function phase_homogenizedC
-
+    module function phase_damage_C(C_homogenized,ph,en) result(C)
+      real(pReal), dimension(3,3,3,3), intent(in)  :: C_homogenized
+      integer,                         intent(in)  :: ph,en
+      real(pReal), dimension(3,3,3,3) :: C
+    end function phase_damage_C
 
     module function phase_f_phi(phi,co,ce) result(f)
       integer, intent(in) :: ce,co
@@ -372,19 +373,6 @@ subroutine phase_init
   call damage_init
   call thermal_init(phases)
 
-
-  phase_source_maxSizeDotState = 0
-  PhaseLoop2:do ph = 1,phases%length
-!--------------------------------------------------------------------------------------------------
-! partition and initialize state
-    plasticState(ph)%state = plasticState(ph)%state0
-    if(damageState(ph)%sizeState > 0) &
-      damageState(ph)%state  = damageState(ph)%state0
-  enddo PhaseLoop2
-
-  phase_source_maxSizeDotState     = maxval(damageState%sizeDotState)
-  phase_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState)
-
 end subroutine phase_init
 
 
@@ -540,8 +528,7 @@ subroutine crystallite_init()
   phases => config_material%get('phase')
 
   do ph = 1, phases%length
-    if (damageState(ph)%sizeState > 0) &
+    if (damageState(ph)%sizeState > 0) allocate(damageState(ph)%subState0,source=damageState(ph)%state0)  ! ToDo: hack
-      allocate(damageState(ph)%subState0,source=damageState(ph)%state0)                 ! ToDo: hack
   enddo
 
   print'(a42,1x,i10)', '    # of elements:                       ', eMax

src/phase_damage.f90

@@ -15,13 +15,15 @@ submodule(phase) damage
     DAMAGE_ANISOBRITTLE_ID
   end enum
 
+  integer :: phase_damage_maxSizeDotState
+
 
   type :: tDataContainer
     real(pReal), dimension(:), allocatable :: phi, d_phi_d_dot_phi
   end type tDataContainer
 
   integer(kind(DAMAGE_UNDEFINED_ID)),     dimension(:), allocatable :: &
-    phase_source                                                                                    !< active sources mechanisms of each phase
+    phase_damage                                                                                    !< active sources mechanisms of each phase
 
   type(tDataContainer), dimension(:), allocatable :: current
 
@@ -126,17 +128,38 @@ module subroutine damage_init
 
   enddo
 
-  allocate(phase_source(phases%length), source = DAMAGE_UNDEFINED_ID)
+  allocate(phase_damage(phases%length), source = DAMAGE_UNDEFINED_ID)
 
   if (damage_active) then
-    where(isobrittle_init()  ) phase_source = DAMAGE_ISOBRITTLE_ID
+    where(isobrittle_init()  ) phase_damage = DAMAGE_ISOBRITTLE_ID
-    where(isoductile_init()  ) phase_source = DAMAGE_ISODUCTILE_ID
+    where(isoductile_init()  ) phase_damage = DAMAGE_ISODUCTILE_ID
-    where(anisobrittle_init()) phase_source = DAMAGE_ANISOBRITTLE_ID
+    where(anisobrittle_init()) phase_damage = DAMAGE_ANISOBRITTLE_ID
   endif
 
+  phase_damage_maxSizeDotState     = maxval(damageState%sizeDotState)
+
 end subroutine damage_init
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief returns the degraded/modified elasticity matrix
+!--------------------------------------------------------------------------------------------------
+module function phase_damage_C(C_homogenized,ph,en) result(C)
+
+  real(pReal), dimension(3,3,3,3), intent(in)  :: C_homogenized
+  integer,                         intent(in)  :: ph,en
+  real(pReal), dimension(3,3,3,3) :: C
+
+  damageType: select case (phase_damage(ph))
+    case (DAMAGE_ISOBRITTLE_ID) damageType
+     C = C_homogenized * damage_phi(ph,en)**2
+    case default damageType
+     C = C_homogenized
+  end select damageType
+
+end function phase_damage_C
+
+
 !----------------------------------------------------------------------------------------------
 !< @brief returns local part of nonlocal damage driving force
 !----------------------------------------------------------------------------------------------
@@ -155,7 +178,7 @@ module function phase_f_phi(phi,co,ce) result(f)
   ph = material_phaseID(co,ce)
   en = material_phaseEntry(co,ce)
 
-  select case(phase_source(ph))
+  select case(phase_damage(ph))
     case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ISODUCTILE_ID,DAMAGE_ANISOBRITTLE_ID)
       f = 1.0_pReal &
         - phi*damageState(ph)%state(1,en)
@@ -187,9 +210,9 @@ module function integrateDamageState(dt,co,ip,el) result(broken)
     size_so
   real(pReal) :: &
     zeta
-  real(pReal), dimension(phase_source_maxSizeDotState) :: &
+  real(pReal), dimension(phase_damage_maxSizeDotState) :: &
     r                                                                                               ! state residuum
-  real(pReal), dimension(phase_source_maxSizeDotState,2) :: source_dotState
+  real(pReal), dimension(phase_damage_maxSizeDotState,2) :: source_dotState
   logical :: &
     converged_
 
@@ -275,10 +298,10 @@ module subroutine damage_results(group,ph)
   integer,          intent(in) :: ph
 
 
-  if (phase_source(ph) /= DAMAGE_UNDEFINED_ID) &
+  if (phase_damage(ph) /= DAMAGE_UNDEFINED_ID) &
     call results_closeGroup(results_addGroup(group//'damage'))
 
-  sourceType: select case (phase_source(ph))
+  sourceType: select case (phase_damage(ph))
 
     case (DAMAGE_ISOBRITTLE_ID) sourceType
       call isobrittle_results(ph,group//'damage/')
@@ -309,7 +332,7 @@ function phase_damage_collectDotState(ph,me) result(broken)
 
   if (damageState(ph)%sizeState > 0) then
 
-    sourceType: select case (phase_source(ph))
+    sourceType: select case (phase_damage(ph))
 
       case (DAMAGE_ISODUCTILE_ID) sourceType
         call isoductile_dotState(ph,me)
@@ -376,7 +399,7 @@ function phase_damage_deltaState(Fe, ph, me) result(broken)
 
   if (damageState(ph)%sizeState == 0) return
 
-   sourceType: select case (phase_source(ph))
+   sourceType: select case (phase_damage(ph))
 
     case (DAMAGE_ISOBRITTLE_ID) sourceType
       call isobrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me)

src/phase_damage_isobrittle.f90

@@ -66,14 +66,14 @@ module function isobrittle_init() result(mySources)
 
         Nmembers = count(material_phaseID==ph)
         call phase_allocateState(damageState(ph),Nmembers,1,1,1)
-        damageState(ph)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal)
+        damageState(ph)%atol = src%get_asFloat('isobrittle_atol',defaultVal=1.0e-3_pReal)
         if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'
 
         end associate
 
 !--------------------------------------------------------------------------------------------------
 !  exit if any parameter is out of range
-        if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isoBrittle)')
+        if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isobrittle)')
       endif
 
   enddo

src/phase_damage_isoductile.f90

@@ -70,14 +70,14 @@ module function isoductile_init() result(mySources)
 
         Nmembers=count(material_phaseID==ph)
         call phase_allocateState(damageState(ph),Nmembers,1,1,0)
-        damageState(ph)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal)
+        damageState(ph)%atol = src%get_asFloat('isoductile_atol',defaultVal=1.0e-3_pReal)
         if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'
 
         end associate
 
 !--------------------------------------------------------------------------------------------------
 !  exit if any parameter is out of range
-        if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isoDuctile)')
+        if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isoductile)')
       endif
   enddo
 

src/phase_mechanical.f90

@@ -41,6 +41,7 @@ submodule(phase) mechanical
   integer(kind(PLASTICITY_undefined_ID)), dimension(:),   allocatable :: &
     phase_plasticity                                                                                !< plasticity of each phase
 
+  integer :: phase_plasticity_maxSizeDotState
 
   interface
 
@@ -294,6 +295,11 @@ module subroutine mechanical_init(materials,phases)
 
   call plastic_init()
 
+  do ph = 1,phases%length
+    plasticState(ph)%state0 = plasticState(ph)%state
+  enddo
+  phase_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState)
+
 
   num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
 

src/phase_mechanical_elastic.f90

@@ -2,18 +2,11 @@ submodule(phase:mechanical) elastic
 
   enum, bind(c); enumerator :: &
     ELASTICITY_UNDEFINED_ID, &
-    ELASTICITY_HOOKE_ID, &
+    ELASTICITY_HOOKE_ID
-    STIFFNESS_DEGRADATION_UNDEFINED_ID, &
-    STIFFNESS_DEGRADATION_DAMAGE_ID
   end enum
  
-  integer, dimension(:), allocatable :: &
-    phase_NstiffnessDegradations
- 
   integer(kind(ELASTICITY_UNDEFINED_ID)), dimension(:),   allocatable :: &
     phase_elasticity                                                                                !< elasticity of each phase
-  integer(kind(STIFFNESS_DEGRADATION_UNDEFINED_ID)),     dimension(:,:), allocatable :: &
-    phase_stiffnessDegradation                                                                      !< active stiffness degradation mechanisms of each phase
 
 contains
 
@@ -24,18 +17,15 @@ module subroutine elastic_init(phases)
     phases
 
   integer :: &
-    ph, &
+    ph
-    stiffDegradationCtr
   class(tNode), pointer :: &
     phase, &
     mech, &
-    elastic, &
+    elastic
-    stiffDegradation
  
   print'(/,a)', ' <<<+-  phase:mechanical:elastic init  -+>>>'
 
   allocate(phase_elasticity(phases%length), source = ELASTICITY_undefined_ID)
-  allocate(phase_NstiffnessDegradations(phases%length),source=0)
   
   do ph = 1, phases%length
     phase   => phases%get(ph)
@@ -46,25 +36,8 @@ module subroutine elastic_init(phases)
     else
       call IO_error(200,ext_msg=elastic%get_asString('type'))
     endif
-    stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList)                      ! check for stiffness degradation mechanisms
-    phase_NstiffnessDegradations(ph) = stiffDegradation%length
   enddo
 
-  allocate(phase_stiffnessDegradation(maxval(phase_NstiffnessDegradations),phases%length), &
-                        source=STIFFNESS_DEGRADATION_undefined_ID)
-
-  if(maxVal(phase_NstiffnessDegradations)/=0) then
-    do ph = 1, phases%length
-      phase => phases%get(ph)
-      mech    => phase%get('mechanical')
-      stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList)
-      do stiffDegradationCtr = 1, stiffDegradation%length
-        if(stiffDegradation%get_asString(stiffDegradationCtr) == 'damage') &
-            phase_stiffnessDegradation(stiffDegradationCtr,ph) = STIFFNESS_DEGRADATION_damage_ID
-      enddo
-    enddo
-  endif
-
 end subroutine elastic_init
 
 
@@ -94,13 +67,7 @@ module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
     i, j
 
   C = math_66toSym3333(phase_homogenizedC(ph,en))
-
+  C = phase_damage_C(C,ph,en)
-  DegradationLoop: do d = 1, phase_NstiffnessDegradations(ph)
-    degradationType: select case(phase_stiffnessDegradation(d,ph))
-      case (STIFFNESS_DEGRADATION_damage_ID) degradationType
-        C = C * damage_phi(ph,en)**2
-    end select degradationType
-  enddo DegradationLoop
 
   E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3)                                                  !< Green-Lagrange strain in unloaded configuration
   S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi))                                        !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration

src/phase_mechanical_plastic_dislotungsten.f90

@@ -254,8 +254,6 @@ module function plastic_dislotungsten_init() result(myPlasticity)
     allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers),         source=0.0_pReal)
     allocate(dst%threshold_stress(prm%sum_N_sl,Nmembers),  source=0.0_pReal)
 
-    plasticState(ph)%state0 = plasticState(ph)%state                                                ! ToDo: this could be done centrally
-
     end associate
 
 !--------------------------------------------------------------------------------------------------
@@ -414,7 +412,7 @@ module subroutine plastic_dislotungsten_results(ph,group)
                                                      'mobile dislocation density','1/m²')
       case('rho_dip')
         if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_dip,trim(prm%output(o)), &
-                                                     'dislocation dipole density''1/m²')
+                                                     'dislocation dipole density','1/m²')
       case('gamma_sl')
         if(prm%sum_N_sl>0) call results_writeDataset(group,stt%gamma_sl,trim(prm%output(o)), &
                                                      'plastic shear','1')

src/phase_mechanical_plastic_dislotwin.f90

@@ -467,8 +467,6 @@ module function plastic_dislotwin_init() result(myPlasticity)
     allocate(dst%tau_r_tr              (prm%sum_N_tr,Nmembers),source=0.0_pReal)
     allocate(dst%V_tr                  (prm%sum_N_tr,Nmembers),source=0.0_pReal)
 
-    plasticState(ph)%state0 = plasticState(ph)%state                                                ! ToDo: this could be done centrally
-
     end associate
 
 !--------------------------------------------------------------------------------------------------

src/phase_mechanical_plastic_isotropic.f90

@@ -140,8 +140,6 @@ module function plastic_isotropic_init() result(myPlasticity)
     ! global alias
     plasticState(ph)%slipRate => plasticState(ph)%dotState(2:2,:)
 
-    plasticState(ph)%state0 = plasticState(ph)%state                                                ! ToDo: this could be done centrally
-
     end associate
 
 !--------------------------------------------------------------------------------------------------

src/phase_mechanical_plastic_kinehardening.f90

@@ -213,8 +213,6 @@ module function plastic_kinehardening_init() result(myPlasticity)
     stt%gamma0 => plasticState(ph)%state     (startIndex  :endIndex  ,:)
     dlt%gamma0 => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:)
 
-    plasticState(ph)%state0 = plasticState(ph)%state                                                  ! ToDo: this could be done centrally
-
     end associate
 
 !--------------------------------------------------------------------------------------------------

src/phase_mechanical_plastic_phenopowerlaw.f90

@@ -265,8 +265,6 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
     if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
 
-    plasticState(ph)%state0 = plasticState(ph)%state                                                ! ToDo: this could be done centrally
-
     end associate
 
 !--------------------------------------------------------------------------------------------------

src/rotations.f90

@@ -178,8 +178,7 @@ subroutine fromQuaternion(self,qu)
   class(rotation), intent(out)          :: self
   real(pReal), dimension(4), intent(in) :: qu
 
-  if (dNeq(norm2(qu),1.0_pReal)) &
+  if (dNeq(norm2(qu),1.0_pReal,1.0e-8_pReal)) call IO_error(402,ext_msg='fromQuaternion')
-    call IO_error(402,ext_msg='fromQuaternion')
 
   self%q = qu
 
@@ -1394,7 +1393,7 @@ end function conjugate_quaternion
 !--------------------------------------------------------------------------------------------------
 !> @brief Check correctness of some rotations functions.
 !--------------------------------------------------------------------------------------------------
-subroutine selfTest
+subroutine selfTest()
 
   type(rotation)                  :: R
   real(pReal), dimension(4)       :: qu, ax, ro
@@ -1405,7 +1404,7 @@ subroutine selfTest
   integer :: i
 
 
-  do i = 1, 10
+  do i = 1, 20
 
 #if defined(__GFORTRAN__) && __GNUC__<9
     if(i<7) cycle
@@ -1433,6 +1432,7 @@ subroutine selfTest
             sin(2.0_pReal*PI*x(1))*A]
       if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal)
     endif
+
     if(.not. quaternion_equal(om2qu(qu2om(qu)),qu)) error stop 'om2qu/qu2om'
     if(.not. quaternion_equal(eu2qu(qu2eu(qu)),qu)) error stop 'eu2qu/qu2eu'
     if(.not. quaternion_equal(ax2qu(qu2ax(qu)),qu)) error stop 'ax2qu/qu2ax'
@@ -1479,6 +1479,8 @@ subroutine selfTest
     if(all(dNeq(R%rotTensor4(R%rotTensor4(t3333),active=.true.),t3333,1.0e-12_pReal))) &
       error stop 'rotTensor4'
 
+    call R%fromQuaternion(qu * (1.0_pReal + merge(+5.e-9_pReal,-5.e-9_pReal, mod(i,2) == 0)))       ! allow reasonable tolerance for ASCII/YAML
+
   enddo
 
   contains