Additional feature to the RGC scheme: The volumetric consistency constraint. This is implemented by penalizing any volumetric discrepancy which occurs due to severe relaxation. The penalty for volume inconsistency is described using a power-law model. The setting for this function can be done in the numerics.config via three new parameters, i.e., maxvoldiscrepancy_RGC, voldiscrepancymod_RGC and discrepancypower_RGC. For monintoring, an extra output variable in prescribed in the material.config in the RGC section, namely "(output) volumediscrepancy".

IO.f90 and numerics.f90 has been modified accordingly to accommodate these changes

code/IO.f90

@@ -1054,9 +1054,11 @@ endfunction
  case (277)
    msg = 'Non-positive relevant mismatch in RGC'
  case (278)
-   msg = 'Non-positive positive definite viscosity model in RGC'
+   msg = 'Non-positive definite viscosity model in RGC'
  case (288)
    msg = 'Non-positive maximum threshold of relaxation change in RGC'
+ case (289)
+   msg = 'Non-positive definite volume discrepancy penalty in RGC'
 
 !* Error message when selected perturbation method is not defined
  case (299)

code/homogenization_RGC.f90

@@ -23,6 +23,8 @@ MODULE homogenization_RGC
  integer(pInt),     dimension(:,:),   allocatable :: homogenization_RGC_Ngrains
  real(pReal),       dimension(:,:),   allocatable :: homogenization_RGC_xiAlpha, &
                                                      homogenization_RGC_ciAlpha
+ real(pReal),       dimension(:),     allocatable :: homogenization_RGC_maxVol0, &
+                                                     homogenization_RGC_vPower0
  character(len=64), dimension(:,:),   allocatable :: homogenization_RGC_output
 
 CONTAINS
@@ -68,6 +70,8 @@ subroutine homogenization_RGC_init(&
  allocate(homogenization_RGC_Ngrains(3,maxNinstance));       homogenization_RGC_Ngrains = 0_pInt
  allocate(homogenization_RGC_ciAlpha(3,maxNinstance));       homogenization_RGC_ciAlpha = 0.0_pReal
  allocate(homogenization_RGC_xiAlpha(3,maxNinstance));       homogenization_RGC_xiAlpha = 0.0_pReal
+ allocate(homogenization_RGC_maxVol0(maxNinstance));         homogenization_RGC_maxVol0 = 0.0_pReal
+ allocate(homogenization_RGC_vPower0(maxNinstance));         homogenization_RGC_vPower0 = 0.0_pReal
  allocate(homogenization_RGC_output(maxval(homogenization_Noutput),maxNinstance)); homogenization_RGC_output = ''
  
  rewind(file)
@@ -106,6 +110,10 @@ subroutine homogenization_RGC_init(&
               homogenization_RGC_ciAlpha(1,i) = IO_floatValue(line,positions,2)
               homogenization_RGC_ciAlpha(2,i) = IO_floatValue(line,positions,3)
               homogenization_RGC_ciAlpha(3,i) = IO_floatValue(line,positions,4)
+       case ('maxvoldiscrepancy')
+              homogenization_RGC_maxVol0(i) = IO_floatValue(line,positions,2)
+       case ('discrepancypower')
+              homogenization_RGC_vPower0(i) = IO_floatValue(line,positions,2)
      end select
    endif
  enddo
@@ -119,10 +127,13 @@ subroutine homogenization_RGC_init(&
        case('constitutivework')
          homogenization_RGC_sizePostResults(i) = &
          homogenization_RGC_sizePostResults(i) + 1
+       case('magnitudemismatch')
+         homogenization_RGC_sizePostResults(i) = &
+         homogenization_RGC_sizePostResults(i) + 1
        case('penaltyenergy')
          homogenization_RGC_sizePostResults(i) = &
          homogenization_RGC_sizePostResults(i) + 1
-       case('magnitudemismatch')
+       case('volumediscrepancy')
          homogenization_RGC_sizePostResults(i) = &
          homogenization_RGC_sizePostResults(i) + 1
      end select
@@ -132,7 +143,8 @@ subroutine homogenization_RGC_init(&
        = 3*(homogenization_RGC_Ngrains(1,i)-1)*homogenization_RGC_Ngrains(2,i)*homogenization_RGC_Ngrains(3,i) &
          + 3*homogenization_RGC_Ngrains(1,i)*(homogenization_RGC_Ngrains(2,i)-1)*homogenization_RGC_Ngrains(3,i) &
          + 3*homogenization_RGC_Ngrains(1,i)*homogenization_RGC_Ngrains(2,i)*(homogenization_RGC_Ngrains(3,i)-1) &
-         + 3_pInt   ! (1) Average constitutive work, (2) Average penalty energy, (3) Overall mismatch
+         + 4_pInt   ! (1) Average constitutive work, (2) Overall mismatch, (3) Average penalty energy, 
+                    ! (4) Volume discrepancy
  enddo
 
  return
@@ -277,11 +289,11 @@ function homogenization_RGC_updateState(&
  integer(pInt), dimension (4) :: intFaceN,intFaceP,faceID
  integer(pInt), dimension (3) :: nGDim,iGr3N,iGr3P,stresLoc
  integer(pInt), dimension (2) :: residLoc
- integer(pInt) homID,i1,i2,i3,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ival,ipert,iGrain
+ integer(pInt) homID,i1,i2,i3,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ival,ipert,iGrain,nGrain
- real(pReal), dimension (3,3,homogenization_maxNgrains) :: R,pF,pR
+ real(pReal), dimension (3,3,homogenization_maxNgrains) :: R,pF,pR,D,pD
  real(pReal), dimension (homogenization_maxNgrains)     :: NN,pNN
  real(pReal), dimension (3)   :: normP,normN,mornP,mornN
- real(pReal) residMax,stresMax,constitutiveWork,penaltyEnergy
+ real(pReal) residMax,stresMax,constitutiveWork,penaltyEnergy,volDiscrep,penDiscrep
  logical error,RGCdebug,RGCdebugJacobi,RGCcheck
 !
  integer(pInt), parameter :: nFace = 6
@@ -296,6 +308,7 @@ function homogenization_RGC_updateState(&
 !* Get the dimension of the cluster (grains and interfaces)
  homID  = homogenization_typeInstance(mesh_element(3,el))
  nGDim  = homogenization_RGC_Ngrains(:,homID)
+ nGrain = homogenization_Ngrains(mesh_element(3,el))
  nIntFaceTot = (nGDim(1)-1)*nGDim(2)*nGDim(3) + nGDim(1)*(nGDim(2)-1)*nGDim(3) &
                + nGDim(1)*nGDim(2)*(nGDim(3)-1)
 
@@ -316,19 +329,26 @@ function homogenization_RGC_updateState(&
 
 !* Stress-like penalty related to mismatch or incompatibility at interfaces
  call homogenization_RGC_stressPenalty(R,NN,F,ip,el,homID)
-!* Debugging the mismatch, stress and penalty of grains
+
- if (RGCdebug) then
+!* Stress-like penalty related to overall volume discrepancy 
-   do iGrain = 1,homogenization_Ngrains(mesh_element(3,el))
+ call homogenization_RGC_volumePenalty(D,volDiscrep,F,avgF,ip,el,homID)
+
+!* Debugging the mismatch, stress and penalties of grains
+ if (el == 1 .and. ip == 1) then
+   do iGrain = 1,nGrain
      write(6,'(x,a30,x,i3,x,a4,x,e14.8)')'Mismatch magnitude of grain(',iGrain,') :',NN(iGrain)
      write(6,*)' '
-     write(6,'(x,a30,x,i3)')'Stress and penalty of grain: ',iGrain
+     write(6,'(x,a30,x,i3)')'Stress and penalties of grain: ',iGrain
      do i = 1,3
-       write(6,'(x,3(e14.8,x),x,3(e14.8,x))')(P(i,j,iGrain), j = 1,3),(R(i,j,iGrain), j = 1,3)
+       write(6,'(x,3(e14.8,x),x,3(e14.8,x),x,3(e14.8,x))')(P(i,j,iGrain), j = 1,3), &
+                                                          (R(i,j,iGrain), j = 1,3), &
+                                                          (D(i,j,iGrain), j = 1,3)
      enddo
      write(6,*)' '
    enddo
  endif
 
+
 !* Compute the residual stress from the balance of traction at all (interior) interfaces
  do iNum = 1,nIntFaceTot
    call homogenization_RGC_interface1to4(faceID,iNum,homID)
@@ -347,8 +367,8 @@ function homogenization_RGC_updateState(&
      tract(iNum,i) = sign(viscModus_RGC*(abs(drelax(i+3*(iNum-1))/dt))**ratePower_RGC, &
                           drelax(i+3*(iNum-1)))
      do j = 1,3
-       tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP))*normP(j) &
+       tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP) + D(i,j,iGrP))*normP(j) &
-                                     + (P(i,j,iGrN) + R(i,j,iGrN))*normN(j)
+                                     + (P(i,j,iGrN) + R(i,j,iGrN) + D(i,j,iGrN))*normN(j)
        resid(i+3*(iNum-1)) = tract(iNum,i)                              ! map into 1D residual array
      enddo
    enddo
@@ -386,25 +406,29 @@ function homogenization_RGC_updateState(&
    endif
 !    write(6,'(x,a,x,i3,x,a6,x,i3,x,a12)')'RGC_updateState: ip',ip,'| el',el,'converged :)'
 !* Then compute/update the state for postResult, i.e., ...
-!* ... the (bulk) constitutive work and penalty energy
+!* ... all energy densities computed by time-integration
    constitutiveWork = state%p(3*nIntFaceTot+1)
-   penaltyEnergy    = state%p(3*nIntFaceTot+2)
+   penaltyEnergy    = state%p(3*nIntFaceTot+3)
    do iGrain = 1,homogenization_Ngrains(mesh_element(3,el))
      do i = 1,3
      do j = 1,3
-       constitutiveWork = constitutiveWork + P(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))
+       constitutiveWork = constitutiveWork + P(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/dble(nGrain)
-       penaltyEnergy    = penaltyEnergy    + R(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))
+       penaltyEnergy    = penaltyEnergy    + R(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/dble(nGrain)
      enddo
      enddo
    enddo
+!* ... the bulk mechanical/constitutive work
    state%p(3*nIntFaceTot+1) = constitutiveWork
-   state%p(3*nIntFaceTot+2) = penaltyEnergy
 !* ... the overall mismatch
-   state%p(3*nIntFaceTot+3) = sum(NN)
+   state%p(3*nIntFaceTot+2) = sum(NN)/dble(nGrain)
-   if (RGCdebug) then
+   state%p(3*nIntFaceTot+3) = penaltyEnergy
+!* ... the volume discrepancy
+   state%p(3*nIntFaceTot+4) = volDiscrep
+   if (el == 1 .and. ip == 1) then
      write(6,'(x,a30,x,e14.8)')'Constitutive work: ',constitutiveWork
-     write(6,'(x,a30,x,e14.8)')'Penalty energy: ',penaltyEnergy
      write(6,'(x,a30,x,e14.8)')'Magnitude mismatch: ',sum(NN)
+     write(6,'(x,a30,x,e14.8)')'Penalty energy: ',penaltyEnergy
+     write(6,'(x,a30,x,e14.8)')'Volume discrepancy: ',volDiscrep
      write(6,*)''
      call flush(6)
    endif
@@ -489,6 +513,7 @@ function homogenization_RGC_updateState(&
    state%p(1:3*nIntFaceTot) = p_relax
    call homogenization_RGC_grainDeformation(pF,F0,avgF,state,el)
    call homogenization_RGC_stressPenalty(pR,pNN,pF,ip,el,homID)
+   call homogenization_RGC_volumePenalty(pD,volDiscrep,pF,avgF,ip,el,homID)
    p_resid = 0.0_pReal
    do iNum = 1,nIntFaceTot
      call homogenization_RGC_interface1to4(faceID,iNum,homID)
@@ -507,7 +532,9 @@ function homogenization_RGC_updateState(&
      do i = 1,3
      do j = 1,3
        p_resid(i+3*(iNum-1)) = p_resid(i+3*(iNum-1)) + (pR(i,j,iGrP) - R(i,j,iGrP))*normP(j) &
-                                                     + (pR(i,j,iGrN) - R(i,j,iGrN))*normN(j)
+                                                     + (pR(i,j,iGrN) - R(i,j,iGrN))*normN(j) &
+                                                     + (pD(i,j,iGrP) - D(i,j,iGrP))*normP(j) &
+                                                     + (pD(i,j,iGrN) - D(i,j,iGrN))*normN(j)
      enddo
      enddo
    enddo
@@ -708,11 +735,14 @@ pure function homogenization_RGC_postResults(&
      case('constitutivework')
        homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+1)
        c = c + 1
+     case('magnitudemismatch')
+       homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+3)
+       c = c + 1
      case('penaltyenergy')
        homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+2)
        c = c + 1
-     case('magnitudemismatch')
+     case('volumediscrepancy')
-       homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+3)
+       homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+4)
        c = c + 1
    end select
  enddo
@@ -722,13 +752,13 @@ pure function homogenization_RGC_postResults(&
 endfunction
 
 !********************************************************************
-! subroutine to calculate stress-like penalty due to mismatch
+! subroutine to calculate stress-like penalty due to deformation mismatch
 !********************************************************************
 subroutine homogenization_RGC_stressPenalty(&
    rPen, &          ! stress-like penalty
    nMis, &          ! total amount of mismatch
 !
-   fDef, &          ! relaxation vectors
+   fDef, &          ! deformation gradients
    ip, &            ! integration point
    el, &            ! element
    homID &          ! homogenization ID
@@ -761,7 +791,7 @@ subroutine homogenization_RGC_stressPenalty(&
  rPen = 0.0_pReal
  nMis = 0.0_pReal
  do iGrain = 1,homogenization_Ngrains(mesh_element(3,el))
-   call homogenization_RGC_equivalentShearMod(muGrain,constitutive_homogenizedC(iGrain,ip,el))
+   call homogenization_RGC_equivalentModuli(muGrain,constitutive_homogenizedC(iGrain,ip,el))
    call homogenization_RGC_grain1to3(iGrain3,iGrain,homID)
 !* Compute the mismatch tensor at all six interfaces
    do iFace = 1,nFace
@@ -779,7 +809,7 @@ subroutine homogenization_RGC_stressPenalty(&
      if (iGNghb3(3) < 1)        iGNghb3(3) = nGDim(3)
      if (iGNghb3(3) > nGDim(3)) iGNghb3(3) = 1
      call homogenization_RGC_grain3to1(iGNghb,iGNghb3,homID)            ! get the grain neighbor
-     call homogenization_RGC_equivalentShearMod(muGNghb,constitutive_homogenizedC(iGNghb,ip,el))
+     call homogenization_RGC_equivalentModuli(muGNghb,constitutive_homogenizedC(iGNghb,ip,el))
      gDef = 0.5_pReal*(fDef(:,:,iGNghb) - fDef(:,:,iGrain))             ! difference in F with the neighbor
      nDefNorm = 0.0_pReal
      nDef = 0.0_pReal
@@ -832,9 +862,69 @@ subroutine homogenization_RGC_stressPenalty(&
 endsubroutine
 
 !********************************************************************
-! subroutine to compute the equivalent shear modulus from the elasticity tensor
+! subroutine to calculate stress-like penalty due to volume discrepancy
 !********************************************************************
-subroutine homogenization_RGC_equivalentShearMod(&
+subroutine homogenization_RGC_volumePenalty(&
+   vPen, &          ! stress-like penalty due to volume
+   vDiscrep, &      ! total volume discrepancy
+!
+   fDef, &          ! deformation gradients
+   fAvg, &          ! overall deformation gradient
+   ip, &            ! integration point
+   el, &            ! element
+   homID &          ! homogenization ID
+  )
+ use prec, only: pReal,pInt,p_vec
+ use mesh, only: mesh_element
+ use constitutive, only: constitutive_homogenizedC
+ use math, only: math_det3x3,math_inv3x3
+ use material, only: homogenization_maxNgrains,homogenization_Ngrains
+ use numerics, only: maxVolDiscr_RGC,volDiscrMod_RGC,volDiscrPow_RGC
+
+ implicit none
+
+!* Definition of variables
+ real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: vPen
+ real(pReal), intent(out)                  :: vDiscrep
+ real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in)  :: fDef
+ real(pReal), dimension (3,3), intent(in)  :: fAvg
+ integer(pInt), intent(in)                 :: ip,el
+ real(pReal), dimension (homogenization_maxNgrains) :: gVol
+ integer(pInt) homID,iGrain,nGrain,i,j
+! 
+ nGrain = homogenization_Ngrains(mesh_element(3,el))
+
+!* Compute volumes of grain and the effective volume and the total volume discrepancy
+ vDiscrep   = math_det3x3(fAvg)
+ do iGrain = 1,nGrain
+   gVol(iGrain) = math_det3x3(fDef(:,:,iGrain))
+   vDiscrep     = vDiscrep - gVol(iGrain)/dble(nGrain)
+ enddo
+
+!* Calculate the stress and penalty due to volume discrepancy
+ vPen       = 0.0_pReal
+ do iGrain = 1,nGrain
+   vPen(:,:,iGrain) = -1.0_pReal*sign(volDiscrMod_RGC*volDiscrPow_RGC/maxVolDiscr_RGC* &
+                      (abs(vDiscrep)/maxVolDiscr_RGC)**(volDiscrPow_RGC - 1.0)*gVol(iGrain)/dble(nGrain),vDiscrep)* &
+                      transpose(math_inv3x3(fDef(:,:,iGrain)))
+
+!* Debugging the stress-like penalty of volume discrepancy
+!    if (ip == 1 .and. el == 1) then
+!      write(6,'(x,a30,i2)')'Volume penalty of grain: ',iGrain
+!      do i = 1,3
+!        write(6,'(x,3(e10.4,x))')(vPen(i,j,iGrain), j = 1,3)
+!      enddo
+!    endif
+ enddo
+
+ return
+
+endsubroutine
+
+!********************************************************************
+! subroutine to compute the equivalent shear and bulk moduli from the elasticity tensor
+!********************************************************************
+subroutine homogenization_RGC_equivalentModuli(&
    shearMod, &       ! equivalent (isotropic) shear modulus
 !
    elasTens &        ! elasticity tensor in Mandel notation

code/material.config

@@ -17,6 +17,7 @@ overproportionality	1.0e+1	1.0e+1	1.0e+1   # typical range between 0.1 (very lar
 (output)		constitutivework
 (output)		penaltyenergy
 (output)		magnitudemismatch
+(output)		volumediscrepancy
 
 [Taylor2]
 type	isostrain

code/numerics.config

@@ -35,5 +35,8 @@ relevantMismatch_RGC      1.0e-5      # minimum threshold of mismatch
 viscosityRate_RGC         1.0e+0      # power (sensitivity rate) of numerical viscosity in RGC scheme
 viscosityModulus_RGC      0.0e+0      # stress modulus of RGC numerical viscosity (zero = without numerical viscosity)
 maxRelaxation_RGC         1.0e+0      # threshold of maximum relaxation vector increment (if exceed this then cutback)
+maxvoldiscrepancy_RGC     1.0e-5      # threshold of maximum volume discrepancy allowed
+voldiscrepancymod_RGC     0.0e+8      # energy modulus of RGC volume discrepancy (zero = without volume discrepancy constraint)
+discrepancypower_RGC      5.0         # powerlaw penalty for volume discrepancy
 
 fixed_seed                1234        # put any number larger than zero, integer, if you want to have a pseudo random distribution

code/numerics.f90

@@ -28,7 +28,7 @@ real(pReal)                     relevantStrain, &                       ! strain
                                 rTol_crystalliteStress, &               ! relative tolerance in crystallite stress loop
                                 aTol_crystalliteStress, &               ! absolute tolerance in crystallite stress loop
 
-!* RGC parameters: added <<<updated 17.11.2009>>>
+!* RGC parameters: added <<<updated 17.12.2009>>>
                                 absTol_RGC, &                           ! absolute tolerance of RGC residuum
                                 relTol_RGC, &                           ! relative tolerance of RGC residuum
                                 absMax_RGC, &                           ! absolute maximum of RGC residuum
@@ -37,7 +37,10 @@ real(pReal)                     relevantStrain, &                       ! strain
                                 xSmoo_RGC, &                            ! RGC penalty smoothing parameter (hyperbolic tangent)
                                 ratePower_RGC, &                        ! power (sensitivity rate) of numerical viscosity in RGC scheme
                                 viscModus_RGC, &                        ! stress modulus of RGC numerical viscosity
-                                maxdRelax_RGC                           ! threshold of maximum relaxation vector increment (if exceed this then cutback)
+                                maxdRelax_RGC, &                        ! threshold of maximum relaxation vector increment (if exceed this then cutback)
+                                maxVolDiscr_RGC, &                      ! threshold of maximum volume discrepancy allowed
+                                volDiscrMod_RGC, &                      ! stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint)
+                                volDiscrPow_RGC                         ! powerlaw penalty for volume discrepancy
 
 !* Random seeding parameters: added <<<updated 27.08.2009>>>
 integer(pInt)                   fixedSeed                            ! fixed seeding for pseudo-random number generator
@@ -101,7 +104,7 @@ subroutine numerics_init()
   rTol_crystalliteStress  = 1.0e-6_pReal
   aTol_crystalliteStress  = 1.0e-8_pReal            ! residuum is in Lp (hence strain on the order of 1e-8 here)
 
-!* RGC parameters: added <<<updated 17.11.2009>>> with moderate setting
+!* RGC parameters: added <<<updated 17.12.2009>>> with moderate setting
   absTol_RGC              = 1.0e+4
   relTol_RGC              = 1.0e-3
   absMax_RGC              = 1.0e+10
@@ -111,6 +114,9 @@ subroutine numerics_init()
   ratePower_RGC           = 1.0e+0  ! Newton viscosity (linear model)
   viscModus_RGC           = 0.0e+0  ! No viscosity is applied
   maxdRelax_RGC           = 1.0e+0
+  maxVolDiscr_RGC         = 1.0e-5  ! tolerance for volume discrepancy allowed
+  volDiscrMod_RGC         = 1.0e+12
+  volDiscrPow_RGC         = 5.0
   
 !* Random seeding parameters: added <<<updated 27.08.2009>>>
   fixedSeed               = 0_pInt
@@ -170,7 +176,7 @@ subroutine numerics_init()
         case ('atol_crystallitestress')
               aTol_crystalliteStress = IO_floatValue(line,positions,2)
 
-!* RGC parameters: added <<<updated 17.11.2009>>>
+!* RGC parameters: added <<<updated 17.12.2009>>>
         case ('atol_rgc')
               absTol_RGC = IO_floatValue(line,positions,2)
         case ('rtol_rgc')
@@ -189,6 +195,12 @@ subroutine numerics_init()
               viscModus_RGC = IO_floatValue(line,positions,2)
         case ('maxrelaxation_rgc')
               maxdRelax_RGC = IO_floatValue(line,positions,2)
+        case ('maxvoldiscrepancy_rgc')
+              maxVolDiscr_RGC = IO_floatValue(line,positions,2)
+        case ('voldiscrepancymod_rgc')
+              volDiscrMod_RGC = IO_floatValue(line,positions,2)
+        case ('discrepancypower_rgc')
+              volDiscrPow_RGC = IO_floatValue(line,positions,2)
 
 !* Random seeding parameters: added <<<updated 27.08.2009>>>
         case ('fixed_seed')
@@ -230,7 +242,7 @@ subroutine numerics_init()
   write(6,'(a24,x,i8)')   'nMPstate:               ',nMPstate
   write(6,*)
 
-!* RGC parameters: added <<<updated 17.11.2009>>>
+!* RGC parameters: added <<<updated 17.12.2009>>>
   write(6,'(a24,x,e8.1)') 'aTol_RGC:             ',absTol_RGC
   write(6,'(a24,x,e8.1)') 'rTol_RGC:             ',relTol_RGC
   write(6,'(a24,x,e8.1)') 'aMax_RGC:             ',absMax_RGC
@@ -240,6 +252,10 @@ subroutine numerics_init()
   write(6,'(a24,x,e8.1)') 'viscosityrate_RGC:    ',ratePower_RGC
   write(6,'(a24,x,e8.1)') 'viscositymodulus_RGC: ',viscModus_RGC
   write(6,'(a24,x,e8.1)') 'maxrelaxation_RGC:    ',maxdRelax_RGC
+  write(6,'(a24,x,e8.1)') 'maxVolDiscrepancy_RGC:',maxVolDiscr_RGC
+  write(6,'(a24,x,e8.1)') 'volDiscrepancyMod_RGC:',volDiscrMod_RGC
+  write(6,'(a24,x,e8.1)') 'discrepancyPower_RGC: ',volDiscrPow_RGC
+
   write(6,*)
 
 !* Random seeding parameters: added <<<updated 27.08.2009>>>
@@ -279,6 +295,9 @@ subroutine numerics_init()
   if (ratePower_RGC < 0.0_pReal)            call IO_error(278)
   if (viscModus_RGC < 0.0_pReal)            call IO_error(278)
   if (maxdRelax_RGC <= 0.0_pReal)           call IO_error(288)
+  if (maxVolDiscr_RGC <= 0.0_pReal)         call IO_error(289)
+  if (volDiscrMod_RGC < 0.0_pReal)          call IO_error(289)
+  if (volDiscrPow_RGC <= 0.0_pReal)         call IO_error(289)
  
   if (fixedSeed <= 0_pInt)                  write(6,'(a)') 'Random is random!'
 endsubroutine