From 3ab5cdc7709551e177c3bafc764ca633e35eeefb Mon Sep 17 00:00:00 2001 From: Denny Tjahjanto Date: Wed, 16 Dec 2009 16:20:53 +0000 Subject: [PATCH] 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 | 4 +- code/homogenization_RGC.f90 | 152 ++++++++++++++++++++++++++++-------- code/material.config | 1 + code/numerics.config | 3 + code/numerics.f90 | 29 +++++-- 5 files changed, 152 insertions(+), 37 deletions(-) diff --git a/code/IO.f90 b/code/IO.f90 index f502f04a5..23699c898 100644 --- a/code/IO.f90 +++ b/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) diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index ff5a12c51..d1ff1bd25 100644 --- a/code/homogenization_RGC.f90 +++ b/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 - real(pReal), dimension (3,3,homogenization_maxNgrains) :: R,pF,pR + 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,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 @@ -294,8 +306,9 @@ function homogenization_RGC_updateState(& RGCdebugJacobi = .false. !* Get the dimension of the cluster (grains and interfaces) - homID = homogenization_typeInstance(mesh_element(3,el)) - nGDim = homogenization_RGC_Ngrains(:,homID) + 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 - do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) + +!* Stress-like penalty related to overall volume discrepancy + 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) & - + (P(i,j,iGrN) + R(i,j,iGrN))*normN(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) + 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)) - penaltyEnergy = penaltyEnergy + R(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))/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) - if (RGCdebug) then + state%p(3*nIntFaceTot+2) = sum(NN)/dble(nGrain) + 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') - homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+3) + case('volumediscrepancy') + 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 diff --git a/code/material.config b/code/material.config index 44ceb53f3..08476d104 100644 --- a/code/material.config +++ b/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 diff --git a/code/numerics.config b/code/numerics.config index a45b68a1c..95a16bcc8 100644 --- a/code/numerics.config +++ b/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 diff --git a/code/numerics.f90 b/code/numerics.f90 index fb3303478..0cd9dd949 100644 --- a/code/numerics.f90 +++ b/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 <<>> +!* RGC parameters: added <<>> 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 <<>> 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 <<>> with moderate setting +!* RGC parameters: added <<>> 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 <<>> fixedSeed = 0_pInt @@ -170,7 +176,7 @@ subroutine numerics_init() case ('atol_crystallitestress') aTol_crystalliteStress = IO_floatValue(line,positions,2) -!* RGC parameters: added <<>> +!* RGC parameters: added <<>> 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 <<>> case ('fixed_seed') @@ -230,7 +242,7 @@ subroutine numerics_init() write(6,'(a24,x,i8)') 'nMPstate: ',nMPstate write(6,*) -!* RGC parameters: added <<>> +!* RGC parameters: added <<>> 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 <<>> @@ -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