diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 06f4356b0..99daef5e5 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -158,12 +158,12 @@ Post_AverageDown: - master - release -Post_General: - stage: postprocessing - script: PostProcessing/test.py - except: - - master - - release +#Post_General: +# stage: postprocessing +# script: PostProcessing/test.py +# except: +# - master +# - release Post_GeometryReconstruction: stage: postprocessing @@ -364,12 +364,12 @@ Phenopowerlaw_singleSlip: - master - release -TextureComponents: - stage: spectral - script: TextureComponents/test.py - except: - - master - - release +#TextureComponents: +# stage: spectral +# script: TextureComponents/test.py +# except: +# - master +# - release ################################################################################################### diff --git a/DAMASK_prerequisites.sh b/DAMASK_prerequisites.sh index c9e165c3c..f0a525c34 100755 --- a/DAMASK_prerequisites.sh +++ b/DAMASK_prerequisites.sh @@ -119,6 +119,9 @@ for executable in mpirun mpiexec; do getDetails $executable '--version' done +firstLevel "CMake" +getDetails cmake --version + firstLevel "Abaqus" cd installation/mods_Abaqus # to have the right environment file for executable in abaqus abq2017 abq2018; do diff --git a/VERSION b/VERSION index 939d675ec..8817d8cc0 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.2-983-g0a648459 +v2.0.2-1124-g0ff9e9fa diff --git a/env/DAMASK.csh b/env/DAMASK.csh index 07b4b6817..2e1dc979c 100644 --- a/env/DAMASK.csh +++ b/env/DAMASK.csh @@ -44,7 +44,6 @@ if ( $?prompt ) then echo "DAMASK $DAMASK_ROOT" echo "Spectral Solver $SOLVER" echo "Post Processing $PROCESSING" - echo "Multithreading DAMASK_NUM_THREADS=$DAMASK_NUM_THREADS" if ( $?PETSC_DIR) then echo "PETSc location $PETSC_DIR" endif @@ -52,8 +51,10 @@ if ( $?prompt ) then echo "MSC.Marc/Mentat $MSC_ROOT" endif echo + echo "Multithreading DAMASK_NUM_THREADS=$DAMASK_NUM_THREADS" echo `limit datasize` echo `limit stacksize` + echo endif setenv DAMASK_NUM_THREADS $DAMASK_NUM_THREADS diff --git a/env/DAMASK.sh b/env/DAMASK.sh index 663e9a4b3..e790ae3cc 100644 --- a/env/DAMASK.sh +++ b/env/DAMASK.sh @@ -88,6 +88,7 @@ if [ ! -z "$PS1" ]; then size=$(( 1024*$(ulimit -s) )); \ print('{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), \ ['bytes','KiB','MiB','GiB','TiB','EiB','ZiB'][int(math.log(size,2) / 10) if size else 0]))") + echo fi export DAMASK_NUM_THREADS diff --git a/env/DAMASK.zsh b/env/DAMASK.zsh index 43f682865..dbfde767d 100644 --- a/env/DAMASK.zsh +++ b/env/DAMASK.zsh @@ -81,6 +81,7 @@ if [ ! -z "$PS1" ]; then size=$(( 1024*$(ulimit -s) )); \ print('{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), \ ['bytes','KiB','MiB','GiB','TiB','EiB','ZiB'][int(math.log(size,2) / 10) if size else 0]))") + echo fi export DAMASK_NUM_THREADS diff --git a/lib/damask/orientation.py b/lib/damask/orientation.py index d0e549b4b..c1862e13c 100644 --- a/lib/damask/orientation.py +++ b/lib/damask/orientation.py @@ -33,7 +33,7 @@ class Quaternion: All methods and naming conventions based on Rowenhorst_etal2015 Convention 1: coordinate frames are right-handed Convention 2: a rotation angle ω is taken to be positive for a counterclockwise rotation - when viewing from the end point of the rotation axis unit vector towards the origin + when viewing from the end point of the rotation axis towards the origin Convention 3: rotations will be interpreted in the passive sense Convention 4: Euler angle triplets are implemented using the Bunge convention, with the angular ranges as [0, 2π],[0, π],[0, 2π] @@ -48,206 +48,138 @@ class Quaternion: """ def __init__(self, - quatArray = [1.0,0.0,0.0,0.0]): - """Initializes to identity if not given""" - self.w, \ - self.x, \ - self.y, \ - self.z = quatArray + quat = None, + q = 1.0, + p = np.zeros(3,dtype=float)): + """Initializes to identity unless specified""" + self.q = quat[0] if quat is not None else q + self.p = np.array(quat[1:4]) if quat is not None else p self.homomorph() def __iter__(self): """Components""" - return iter([self.w,self.x,self.y,self.z]) + return iter(self.asList()) def __copy__(self): - """Create copy""" - Q = Quaternion([self.w,self.x,self.y,self.z]) - return Q + """Copy""" + return self.__class__(q=self.q,p=self.p.copy()) copy = __copy__ def __repr__(self): - """Readbable string""" - return 'Quaternion(real=%+.6f, imag=<%+.6f, %+.6f, %+.6f>)' % \ - (self.w, self.x, self.y, self.z) + """Readable string""" + return 'Quaternion(real={q:+.6f}, imag=<{p[0]:+.6f}, {p[1]:+.6f}, {p[2]:+.6f}>)'.format(q=self.q,p=self.p) def __pow__(self, exponent): """Power""" - omega = math.acos(self.w) - vRescale = math.sin(exponent*omega)/math.sin(omega) - Q = Quaternion() - Q.w = math.cos(exponent*omega) - Q.x = self.x * vRescale - Q.y = self.y * vRescale - Q.z = self.z * vRescale - return Q + omega = math.acos(self.q) + return self.__class__(q= math.cos(exponent*omega), + p=self.p * math.sin(exponent*omega)/math.sin(omega)) def __ipow__(self, exponent): """In-place power""" - omega = math.acos(self.w) - vRescale = math.sin(exponent*omega)/math.sin(omega) - self.w = np.cos(exponent*omega) - self.x *= vRescale - self.y *= vRescale - self.z *= vRescale + omega = math.acos(self.q) + self.q = math.cos(exponent*omega) + self.p *= math.sin(exponent*omega)/math.sin(omega) return self def __mul__(self, other): """Multiplication""" + # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 + P = -1.0 try: # quaternion - Aw = self.w - Ax = self.x - Ay = self.y - Az = self.z - Bw = other.w - Bx = other.x - By = other.y - Bz = other.z - Q = Quaternion() - Q.w = - Ax * Bx - Ay * By - Az * Bz + Aw * Bw - Q.x = + Ax * Bw + Ay * Bz - Az * By + Aw * Bx - Q.y = - Ax * Bz + Ay * Bw + Az * Bx + Aw * By - Q.z = + Ax * By - Ay * Bx + Az * Bw + Aw * Bz - return Q + return self.__class__(q=self.q*other.q - np.dot(self.p,other.p), + p=self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p)) except: pass - try: # vector (perform active rotation, i.e. q*v*q.conjugated) - w = self.w - x = self.x - y = self.y - z = self.z - Vx = other[0] - Vy = other[1] - Vz = other[2] + try: # vector (perform passive rotation) + ( x, y, z) = self.p + (Vx,Vy,Vz) = other[0:3] + A = self.q*self.q - np.dot(self.p,self.p) + B = 2.0 * (x*Vx + y*Vy + z*Vz) + C = 2.0 * P*self.q - return np.array([\ - w * w * Vx + 2 * y * w * Vz - 2 * z * w * Vy + \ - x * x * Vx + 2 * y * x * Vy + 2 * z * x * Vz - \ - z * z * Vx - y * y * Vx, - 2 * x * y * Vx + y * y * Vy + 2 * z * y * Vz + \ - 2 * w * z * Vx - z * z * Vy + w * w * Vy - \ - 2 * x * w * Vz - x * x * Vy, - 2 * x * z * Vx + 2 * y * z * Vy + \ - z * z * Vz - 2 * w * y * Vx - y * y * Vz + \ - 2 * w * x * Vy - x * x * Vz + w * w * Vz ]) + return np.array([ + A*Vx + B*x + C*(y*Vz - z*Vy), + A*Vy + B*y + C*(z*Vx - x*Vz), + A*Vz + B*z + C*(x*Vy - y*Vx), + ]) except: pass try: # scalar - Q = self.copy() - Q.w *= other - Q.x *= other - Q.y *= other - Q.z *= other - return Q + return self.__class__(q=self.q*other, + p=self.p*other) except: return self.copy() def __imul__(self, other): """In-place multiplication""" + # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 + P = -1.0 try: # Quaternion - Aw = self.w - Ax = self.x - Ay = self.y - Az = self.z - Bw = other.w - Bx = other.x - By = other.y - Bz = other.z - self.w = - Ax * Bx - Ay * By - Az * Bz + Aw * Bw - self.x = + Ax * Bw + Ay * Bz - Az * By + Aw * Bx - self.y = - Ax * Bz + Ay * Bw + Az * Bx + Aw * By - self.z = + Ax * By - Ay * Bx + Az * Bw + Aw * Bz + self.q = self.q*other.q - np.dot(self.p,other.p) + self.p = self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p) except: pass return self def __div__(self, other): """Division""" if isinstance(other, (int,float)): - w = self.w / other - x = self.x / other - y = self.y / other - z = self.z / other - return self.__class__([w,x,y,z]) + return self.__class__(q=self.q / other, + p=self.p / other) else: return NotImplemented def __idiv__(self, other): """In-place division""" if isinstance(other, (int,float)): - self.w /= other - self.x /= other - self.y /= other - self.z /= other + self.q /= other + self.p /= other return self def __add__(self, other): """Addition""" if isinstance(other, Quaternion): - w = self.w + other.w - x = self.x + other.x - y = self.y + other.y - z = self.z + other.z - return self.__class__([w,x,y,z]) + return self.__class__(q=self.q + other.q, + p=self.p + other.p) else: return NotImplemented def __iadd__(self, other): """In-place addition""" if isinstance(other, Quaternion): - self.w += other.w - self.x += other.x - self.y += other.y - self.z += other.z + self.q += other.q + self.p += other.p return self def __sub__(self, other): """Subtraction""" if isinstance(other, Quaternion): - Q = self.copy() - Q.w -= other.w - Q.x -= other.x - Q.y -= other.y - Q.z -= other.z - return Q + return self.__class__(q=self.q - other.q, + p=self.p - other.p) else: - return self.copy() + return NotImplemented def __isub__(self, other): """In-place subtraction""" if isinstance(other, Quaternion): - self.w -= other.w - self.x -= other.x - self.y -= other.y - self.z -= other.z + self.q -= other.q + self.p -= other.p return self def __neg__(self): """Additive inverse""" - self.w = -self.w - self.x = -self.x - self.y = -self.y - self.z = -self.z + self.q = -self.q + self.p = -self.p return self def __abs__(self): """Norm""" - return math.sqrt(self.w ** 2 + \ - self.x ** 2 + \ - self.y ** 2 + \ - self.z ** 2) + return math.sqrt(self.q ** 2 + np.dot(self.p,self.p)) magnitude = __abs__ def __eq__(self,other): """Equal at e-8 precision""" - return (abs(self.w-other.w) < 1e-8 and \ - abs(self.x-other.x) < 1e-8 and \ - abs(self.y-other.y) < 1e-8 and \ - abs(self.z-other.z) < 1e-8) \ - or \ - (abs(-self.w-other.w) < 1e-8 and \ - abs(-self.x-other.x) < 1e-8 and \ - abs(-self.y-other.y) < 1e-8 and \ - abs(-self.z-other.z) < 1e-8) + return (self-other).magnitude() < 1e-8 or (-self-other).magnitude() < 1e-8 def __ne__(self,other): """Not equal at e-8 precision""" @@ -255,31 +187,26 @@ class Quaternion: def __cmp__(self,other): """Linear ordering""" - return (self.Rodrigues()>other.Rodrigues()) - (self.Rodrigues() np.linalg.norm(other.asRodrigues()) else 0) \ + - (1 if np.linalg.norm(self.asRodrigues()) < np.linalg.norm(other.asRodrigues()) else 0) def magnitude_squared(self): - return self.w ** 2 + \ - self.x ** 2 + \ - self.y ** 2 + \ - self.z ** 2 + return self.q ** 2 + np.dot(self.p,self.p) def identity(self): - self.w = 1. - self.x = 0. - self.y = 0. - self.z = 0. + self.q = 1. + self.p = np.zeros(3,dtype=float) return self def normalize(self): d = self.magnitude() if d > 0.0: - self /= d + self.q /= d + self.p /= d return self def conjugate(self): - self.x = -self.x - self.y = -self.y - self.z = -self.z + self.p = -self.p return self def inverse(self): @@ -290,11 +217,9 @@ class Quaternion: return self def homomorph(self): - if self.w < 0.0: - self.w = -self.w - self.x = -self.x - self.y = -self.y - self.z = -self.z + if self.q < 0.0: + self.q = -self.q + self.p = -self.p return self def normalized(self): @@ -310,27 +235,35 @@ class Quaternion: return self.copy().homomorph() def asList(self): - return [i for i in self] + return [self.q]+list(self.p) def asM(self): # to find Averaging Quaternions (see F. Landis Markley et al.) - return np.outer([i for i in self],[i for i in self]) - + return np.outer(self.asList(),self.asList()) + def asMatrix(self): - qbarhalf = 0.5*(self.w**2 - self.x**2 - self.y**2 - self.z**2) + # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 + P = -1.0 + qbarhalf = 0.5*(self.q**2 - np.dot(self.p,self.p)) return 2.0*np.array( - [[ qbarhalf + self.x**2 , self.x*self.y - self.w*self.z, self.x*self.z + self.w*self.y], - [ self.x*self.y + self.w*self.z, qbarhalf + self.y**2 , self.y*self.z - self.w*self.x], - [ self.x*self.z - self.w*self.y, self.y*self.z + self.w*self.x, qbarhalf + self.z**2 ], + [[ qbarhalf + self.p[0]**2 , + self.p[0]*self.p[1] -P* self.q*self.p[2], + self.p[0]*self.p[2] +P* self.q*self.p[1] ], + [ self.p[0]*self.p[1] +P* self.q*self.p[2], + qbarhalf + self.p[1]**2 , + self.p[1]*self.p[2] -P* self.q*self.p[0] ], + [ self.p[0]*self.p[2] -P* self.q*self.p[1], + self.p[1]*self.p[2] +P* self.q*self.p[0], + qbarhalf + self.p[2]**2 ], ]) def asAngleAxis(self, degrees = False): - if self.w > 1: + if self.q > 1.: self.normalize() - s = math.sqrt(1. - self.w**2) - x = 2*self.w**2 - 1. - y = 2*self.w * s + s = math.sqrt(1. - self.q**2) + x = 2*self.q**2 - 1. + y = 2*self.q * s angle = math.atan2(y,x) if angle < 0.0: @@ -338,26 +271,28 @@ class Quaternion: s *= -1. return (np.degrees(angle) if degrees else angle, - np.array([1.0, 0.0, 0.0] if np.abs(angle) < 1e-6 else [self.x / s, self.y / s, self.z / s])) + np.array([1.0, 0.0, 0.0] if np.abs(angle) < 1e-6 else self.p / s)) def asRodrigues(self): - return np.inf*np.ones(3) if self.w == 0.0 else np.array([self.x, self.y, self.z])/self.w + return np.inf*np.ones(3) if self.q == 0.0 else self.p/self.q def asEulers(self, degrees = False): """Orientation as Bunge-Euler angles.""" - q03 = self.w**2+self.z**2 - q12 = self.x**2+self.y**2 + # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 + P = -1.0 + q03 = self.q**2 + self.p[2]**2 + q12 = self.p[0]**2 + self.p[1]**2 chi = np.sqrt(q03*q12) if abs(chi) < 1e-10 and abs(q12) < 1e-10: - eulers = np.array([math.atan2(-2*self.w*self.z,self.w**2-self.z**2),0,0]) + eulers = np.array([math.atan2(-2*P*self.q*self.p[2],self.q**2-self.p[2]**2),0,0]) elif abs(chi) < 1e-10 and abs(q03) < 1e-10: - eulers = np.array([math.atan2( 2*self.x*self.y,self.x**2-self.y**2),np.pi,0]) + eulers = np.array([math.atan2( 2 *self.p[0]*self.p[1],self.p[0]**2-self.p[1]**2),np.pi,0]) else: - eulers = np.array([math.atan2((self.x*self.z-self.w*self.y)/chi,(-self.w*self.x-self.y*self.z)/chi), + eulers = np.array([math.atan2((self.p[0]*self.p[2]-P*self.q*self.p[1])/chi,(-P*self.q*self.p[0]-self.p[1]*self.p[2])/chi), math.atan2(2*chi,q03-q12), - math.atan2((self.w*self.y+self.x*self.z)/chi,( self.y*self.z-self.w*self.x)/chi), + math.atan2((P*self.q*self.p[1]+self.p[0]*self.p[2])/chi,( self.p[1]*self.p[2]-P*self.q*self.p[0])/chi), ]) return np.degrees(eulers) if degrees else eulers @@ -371,25 +306,26 @@ class Quaternion: @classmethod def fromRandom(cls,randomSeed = None): + import binascii if randomSeed is None: - randomSeed = int(os.urandom(4).encode('hex'), 16) + randomSeed = int(binascii.hexlify(os.urandom(4)),16) np.random.seed(randomSeed) r = np.random.random(3) w = math.cos(2.0*math.pi*r[0])*math.sqrt(r[2]) x = math.sin(2.0*math.pi*r[1])*math.sqrt(1.0-r[2]) y = math.cos(2.0*math.pi*r[1])*math.sqrt(1.0-r[2]) z = math.sin(2.0*math.pi*r[0])*math.sqrt(r[2]) - return cls([w,x,y,z]) + return cls(quat=[w,x,y,z]) @classmethod def fromRodrigues(cls, rodrigues): if not isinstance(rodrigues, np.ndarray): rodrigues = np.array(rodrigues) - halfangle = math.atan(np.linalg.norm(rodrigues)) + norm = np.linalg.norm(rodrigues) + halfangle = math.atan(norm) + s = math.sin(halfangle) c = math.cos(halfangle) - w = c - x,y,z = rodrigues/c - return cls([w,x,y,z]) + return cls(q=c,p=s*rodrigues/norm) @classmethod @@ -397,22 +333,19 @@ class Quaternion: angle, axis, degrees = False): - if not isinstance(axis, np.ndarray): axis = np.array(axis,dtype='d') + if not isinstance(axis, np.ndarray): axis = np.array(axis,dtype=float) axis = axis.astype(float)/np.linalg.norm(axis) angle = np.radians(angle) if degrees else angle s = math.sin(0.5 * angle) - w = math.cos(0.5 * angle) - x = axis[0] * s - y = axis[1] * s - z = axis[2] * s - return cls([w,x,y,z]) + c = math.cos(0.5 * angle) + return cls(q=c,p=axis*s) @classmethod def fromEulers(cls, eulers, degrees = False): - if not isinstance(eulers, np.ndarray): eulers = np.array(eulers,dtype='d') + if not isinstance(eulers, np.ndarray): eulers = np.array(eulers,dtype=float) eulers = np.radians(eulers) if degrees else eulers sigma = 0.5*(eulers[0]+eulers[2]) @@ -420,11 +353,13 @@ class Quaternion: c = np.cos(0.5*eulers[1]) s = np.sin(0.5*eulers[1]) - w = c * np.cos(sigma) - x = -s * np.cos(delta) - y = -s * np.sin(delta) - z = -c * np.sin(sigma) - return cls([w,x,y,z]) + # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 + P = -1.0 + w = c * np.cos(sigma) + x = -P * s * np.cos(delta) + y = -P * s * np.sin(delta) + z = -P * c * np.sin(sigma) + return cls(quat=[w,x,y,z]) # Modified Method to calculate Quaternion from Orientation Matrix, @@ -435,16 +370,18 @@ class Quaternion: if m.shape != (3,3) and np.prod(m.shape) == 9: m = m.reshape(3,3) - w = 0.5*math.sqrt(1.+m[0,0]+m[1,1]+m[2,2]) - x = 0.5*math.sqrt(1.+m[0,0]-m[1,1]-m[2,2]) - y = 0.5*math.sqrt(1.-m[0,0]+m[1,1]-m[2,2]) - z = 0.5*math.sqrt(1.-m[0,0]-m[1,1]+m[2,2]) + # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 + P = -1.0 + w = 0.5*math.sqrt(1.+m[0,0]+m[1,1]+m[2,2]) + x = P*0.5*math.sqrt(1.+m[0,0]-m[1,1]-m[2,2]) + y = P*0.5*math.sqrt(1.-m[0,0]+m[1,1]-m[2,2]) + z = P*0.5*math.sqrt(1.-m[0,0]-m[1,1]+m[2,2]) x *= -1 if m[2,1] < m[1,2] else 1 y *= -1 if m[0,2] < m[2,0] else 1 z *= -1 if m[1,0] < m[0,1] else 1 - return cls( np.array([w,x,y,z])/math.sqrt(w**2 + x**2 + y**2 + z**2)) + return cls(quat=np.array([w,x,y,z])/math.sqrt(w**2 + x**2 + y**2 + z**2)) @classmethod @@ -458,36 +395,30 @@ class Quaternion: assert isinstance(q1, Quaternion) and isinstance(q2, Quaternion) Q = cls() - costheta = q1.w * q2.w + q1.x * q2.x + q1.y * q2.y + q1.z * q2.z + costheta = q1.q*q2.q + np.dot(q1.p,q2.p) if costheta < 0.: costheta = -costheta q1 = q1.conjugated() - elif costheta > 1: - costheta = 1 + elif costheta > 1.: + costheta = 1. theta = math.acos(costheta) if abs(theta) < 0.01: - Q.w = q2.w - Q.x = q2.x - Q.y = q2.y - Q.z = q2.z + Q.q = q2.q + Q.p = q2.p return Q sintheta = math.sqrt(1.0 - costheta * costheta) if abs(sintheta) < 0.01: - Q.w = (q1.w + q2.w) * 0.5 - Q.x = (q1.x + q2.x) * 0.5 - Q.y = (q1.y + q2.y) * 0.5 - Q.z = (q1.z + q2.z) * 0.5 + Q.q = (q1.q + q2.q) * 0.5 + Q.p = (q1.p + q2.p) * 0.5 return Q - ratio1 = math.sin((1 - t) * theta) / sintheta - ratio2 = math.sin(t * theta) / sintheta + ratio1 = math.sin((1.0 - t) * theta) / sintheta + ratio2 = math.sin( t * theta) / sintheta - Q.w = q1.w * ratio1 + q2.w * ratio2 - Q.x = q1.x * ratio1 + q2.x * ratio2 - Q.y = q1.y * ratio1 + q2.y * ratio2 - Q.z = q1.z * ratio1 + q2.z * ratio2 + Q.q = q1.q * ratio1 + q2.q * ratio2 + Q.p = q1.p * ratio1 + q2.p * ratio2 return Q @@ -513,7 +444,7 @@ class Symmetry: def __repr__(self): """Readbable string""" - return '%s' % (self.lattice) + return '{}'.format(self.lattice) def __eq__(self, other): @@ -526,7 +457,7 @@ class Symmetry: def __cmp__(self,other): """Linear ordering""" - myOrder = Symmetry.lattices.index(self.lattice) + myOrder = Symmetry.lattices.index(self.lattice) otherOrder = Symmetry.lattices.index(other.lattice) return (myOrder > otherOrder) - (myOrder < otherOrder) @@ -722,7 +653,7 @@ class Symmetry: else: return True - v = np.array(vector,dtype = float) + v = np.array(vector,dtype=float) if proper: # check both improper ... theComponents = np.dot(basis['improper'],v) inSST = np.all(theComponents >= 0.0) @@ -737,10 +668,10 @@ class Symmetry: if color: # have to return color array if inSST: rgb = np.power(theComponents/np.linalg.norm(theComponents),0.5) # smoothen color ramps - rgb = np.minimum(np.ones(3,'d'),rgb) # limit to maximum intensity + rgb = np.minimum(np.ones(3,dtype=float),rgb) # limit to maximum intensity rgb /= max(rgb) # normalize to (HS)V = 1 else: - rgb = np.zeros(3,'d') + rgb = np.zeros(3,dtype=float) return (inSST,rgb) else: return inSST @@ -780,8 +711,9 @@ class Orientation: self.quaternion = Quaternion.fromRodrigues(Rodrigues) elif isinstance(quaternion, Quaternion): # based on given quaternion self.quaternion = quaternion.homomorphed() - elif isinstance(quaternion, np.ndarray) and quaternion.shape == (4,): # based on given quaternion-like array - self.quaternion = Quaternion(quaternion).homomorphed() + elif (isinstance(quaternion, np.ndarray) and quaternion.shape == (4,)) or \ + (isinstance(quaternion, list) and len(quaternion) == 4 ): # based on given quaternion-like array + self.quaternion = Quaternion(quat=quaternion).homomorphed() self.symmetry = Symmetry(symmetry) @@ -794,10 +726,12 @@ class Orientation: def __repr__(self): """Value as all implemented representations""" - return 'Symmetry: %s\n' % (self.symmetry) + \ - 'Quaternion: %s\n' % (self.quaternion) + \ - 'Matrix:\n%s\n' % ( '\n'.join(['\t'.join(map(str,self.asMatrix()[i,:])) for i in range(3)]) ) + \ - 'Bunge Eulers / deg: %s' % ('\t'.join(map(str,self.asEulers(degrees=True))) ) + return '\n'.join([ + 'Symmetry: {}'.format(self.symmetry), + 'Quaternion: {}'.format(self.quaternion), + 'Matrix:\n{}'.format( '\n'.join(['\t'.join(list(map(str,self.asMatrix()[i,:]))) for i in range(3)]) ), + 'Bunge Eulers / deg: {}'.format('\t'.join(list(map(str,self.asEulers(degrees=True)))) ), + ]) def asQuaternion(self): return self.quaternion.asList() @@ -927,7 +861,7 @@ class Orientation: M = closest.quaternion.asM() * n if i == 0 else M + closest.quaternion.asM() * n # noqa add (multiples) of this orientation to average noqa eig, vec = np.linalg.eig(M/N) - return Orientation(quaternion = Quaternion(quatArray = np.real(vec.T[eig.argmax()])), + return Orientation(quaternion = Quaternion(quat = np.real(vec.T[eig.argmax()])), symmetry = reference.symmetry.lattice) diff --git a/processing/post/addGrainID.py b/processing/post/addGrainID.py index 45034034b..c3b98f4e6 100755 --- a/processing/post/addGrainID.py +++ b/processing/post/addGrainID.py @@ -200,9 +200,9 @@ for name in filenames: if gID != -1 and gID not in alreadyChecked: # indexed point belonging to a grain not yet tested? alreadyChecked[gID] = True # remember not to check again disorientation = o.disorientation(orientations[gID],SST = False)[0] # compare against other orientation - if disorientation.quaternion.w > cos_disorientation: # within threshold ... + if disorientation.quaternion.q > cos_disorientation: # within threshold ... candidates.append(gID) # remember as potential candidate - if disorientation.quaternion.w >= bestDisorientation.w: # ... and better than current best? + if disorientation.quaternion.q >= bestDisorientation.q: # ... and better than current best? matched = True matchedID = gID # remember that grain bestDisorientation = disorientation.quaternion diff --git a/processing/pre/geom_addPrimitive.py b/processing/pre/geom_addPrimitive.py index 35000d8bf..54de558f7 100755 --- a/processing/pre/geom_addPrimitive.py +++ b/processing/pre/geom_addPrimitive.py @@ -64,11 +64,11 @@ if options.dimension is None: parser.error('no dimension specified.') if options.angleaxis is not None: options.angleaxis = list(map(float,options.angleaxis)) - rotation = damask.Quaternion().fromAngleAxis(np.radians(options.angleaxis[0]) if options.degrees else options.angleaxis[0], - options.angleaxis[1:4]) + rotation = damask.Quaternion.fromAngleAxis(np.radians(options.angleaxis[0]) if options.degrees else options.angleaxis[0], + options.angleaxis[1:4]) elif options.quaternion is not None: - options.quaternion = map(float,options.quaternion) - rotation = damask.Quaternion(options.quaternion) + options.quaternion = list(map(float,options.quaternion)) + rotation = damask.Quaternion(quat=options.quaternion) else: rotation = damask.Quaternion() diff --git a/processing/pre/geom_fromTable.py b/processing/pre/geom_fromTable.py index a2f71d0c9..e1157d325 100755 --- a/processing/pre/geom_fromTable.py +++ b/processing/pre/geom_fromTable.py @@ -43,7 +43,7 @@ parser.add_option('-e', '--eulers', parser.add_option('-d', '--degrees', dest = 'degrees', action = 'store_true', - help = 'angles are given in degrees [%default]') + help = 'all angles are in degrees') parser.add_option('-m', '--matrix', dest = 'matrix', type = 'string', metavar = 'string', @@ -71,7 +71,7 @@ parser.add_option('--axes', parser.add_option('-s', '--symmetry', dest = 'symmetry', action = 'extend', metavar = '', - help = 'crystal symmetry %default {{{}}} '.format(', '.join(damask.Symmetry.lattices[1:]))) + help = 'crystal symmetry of each phase %default {{{}}} '.format(', '.join(damask.Symmetry.lattices[1:]))) parser.add_option('--homogenization', dest = 'homogenization', type = 'int', metavar = 'int', @@ -234,7 +234,7 @@ for name in filenames: o = damask.Orientation(Eulers = myData[colOri:colOri+3]*toRadians, symmetry = mySym) elif inputtype == 'matrix': - o = damask.Orientation(matrix = myData[colOri:colOri+9].reshape(3,3).transpose(), + o = damask.Orientation(matrix = myData[colOri:colOri+9].reshape(3,3), symmetry = mySym) elif inputtype == 'frame': o = damask.Orientation(matrix = np.hstack((myData[colOri[0]:colOri[0]+3], @@ -246,7 +246,7 @@ for name in filenames: o = damask.Orientation(quaternion = myData[colOri:colOri+4], symmetry = mySym) - cos_disorientations = -np.ones(1,dtype='f') # largest possible disorientation + cos_disorientations = -np.ones(1,dtype=float) # largest possible disorientation closest_grain = -1 # invalid neighbor if options.tolerance > 0.0: # only try to compress orientations if asked to @@ -258,7 +258,7 @@ for name in filenames: if len(grains) > 0: # check immediate neighborhood first cos_disorientations = np.array([o.disorientation(orientations[grainID], - SST = False)[0].quaternion.w \ + SST = False)[0].quaternion.q \ for grainID in grains]) # store disorientation per grainID closest_grain = np.argmax(cos_disorientations) # grain among grains with closest orientation to myself match = 'local' @@ -269,7 +269,7 @@ for name in filenames: if len(grains) > 0: cos_disorientations = np.array([o.disorientation(orientations[grainID], - SST = False)[0].quaternion.w \ + SST = False)[0].quaternion.q \ for grainID in grains]) # store disorientation per grainID closest_grain = np.argmax(cos_disorientations) # grain among grains with closest orientation to myself match = 'global' diff --git a/processing/pre/hybridIA_linODFsampling.py b/processing/pre/hybridIA_linODFsampling.py index acedc3cca..cf1a473cf 100755 --- a/processing/pre/hybridIA_linODFsampling.py +++ b/processing/pre/hybridIA_linODFsampling.py @@ -244,7 +244,7 @@ for name in filenames: continue damask.util.report(scriptName,name) - randomSeed = int(os.urandom(4).encode('hex'), 16) if options.randomSeed is None else options.randomSeed # random seed per file for second phase + randomSeed = int(os.urandom(4).hex(), 16) if options.randomSeed is None else options.randomSeed # random seed per file for second phase random.seed(randomSeed) # ------------------------------------------ read header and data --------------------------------- diff --git a/src/IO.f90 b/src/IO.f90 index a76959a44..af59b11b9 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -191,7 +191,9 @@ recursive function IO_recursiveRead(fileName,cnt) result(fileContent) l,i, & myStat - if (merge(cnt,0_pInt,present(cnt))>10_pInt) call IO_error(106_pInt,ext_msg=trim(fileName)) + if (present(cnt)) then + if (cnt>10_pInt) call IO_error(106_pInt,ext_msg=trim(fileName)) + endif !-------------------------------------------------------------------------------------------------- ! read data as stream @@ -684,7 +686,11 @@ function IO_stringValue(string,chunkPos,myChunk,silent) logical :: warn - warn = merge(silent,.false.,present(silent)) + if (present(silent)) then + warn = silent + else + warn = .false. + endif IO_stringValue = '' valuePresent: if (myChunk > chunkPos(1) .or. myChunk < 1_pInt) then diff --git a/src/config.f90 b/src/config.f90 index 4d5a76432..7ae800f30 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -513,8 +513,12 @@ character(len=65536) function getString(this,key,defaultVal,raw) type(tPartitionedStringList), pointer :: item logical :: found, & whole + if (present(raw)) then + whole = raw + else + whole = .false. + endif - whole = merge(raw,.false.,present(raw)) ! whole string or white space splitting found = present(defaultVal) if (found) then getString = trim(defaultVal) @@ -661,7 +665,11 @@ function getStrings(this,key,defaultVal,requiredShape,raw) cumulative cumulative = (key(1:1) == '(' .and. key(len_trim(key):len_trim(key)) == ')') - whole = merge(raw,.false.,present(raw)) + if (present(raw)) then + whole = raw + else + whole = .false. + endif found = .false. item => this diff --git a/src/constitutive.f90 b/src/constitutive.f90 index c4d2cacbf..eca8af08a 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -513,7 +513,7 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType of = phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(material_phase(ipc,ip,el)) - call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMp, Mp,instance,of) + call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMp,Mp,instance,of) case (PLASTICITY_KINEHARDENING_ID) plasticityType call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el) @@ -525,9 +525,9 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp), & - temperature(ho)%p(tme),ipc,ip,el) - dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + of = phasememberAt(ipc,ip,el) + instance = phase_plasticityInstance(material_phase(ipc,ip,el)) + call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of) case (PLASTICITY_DISLOUCLA_ID) plasticityType call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp), & @@ -922,8 +922,9 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac call plastic_kinehardening_dotState(math_Mandel33to6(Mp),ipc,ip,el) case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_dotState (math_Mandel33to6(Mp),temperature(ho)%p(tme), & - ipc,ip,el) + of = phasememberAt(ipc,ip,el) + instance = phase_plasticityInstance(material_phase(ipc,ip,el)) + call plastic_dislotwin_dotState (Mp,temperature(ho)%p(tme),instance,of) case (PLASTICITY_DISLOUCLA_ID) plasticityType call plastic_disloucla_dotState (math_Mandel33to6(Mp),temperature(ho)%p(tme), & @@ -1135,20 +1136,27 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el) case (PLASTICITY_ISOTROPIC_ID) plasticityType constitutive_postResults(startPos:endPos) = & plastic_isotropic_postResults(S6,ipc,ip,el) + case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType of = phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(material_phase(ipc,ip,el)) constitutive_postResults(startPos:endPos) = & plastic_phenopowerlaw_postResults(Mp,instance,of) + case (PLASTICITY_KINEHARDENING_ID) plasticityType constitutive_postResults(startPos:endPos) = & plastic_kinehardening_postResults(S6,ipc,ip,el) + case (PLASTICITY_DISLOTWIN_ID) plasticityType + of = phasememberAt(ipc,ip,el) + instance = phase_plasticityInstance(material_phase(ipc,ip,el)) constitutive_postResults(startPos:endPos) = & - plastic_dislotwin_postResults(S6,temperature(ho)%p(tme),ipc,ip,el) + plastic_dislotwin_postResults(Mp,temperature(ho)%p(tme),instance,of) + case (PLASTICITY_DISLOUCLA_ID) plasticityType constitutive_postResults(startPos:endPos) = & plastic_disloucla_postResults(S6,temperature(ho)%p(tme),ipc,ip,el) + case (PLASTICITY_NONLOCAL_ID) plasticityType constitutive_postResults(startPos:endPos) = & plastic_nonlocal_postResults (S6,FeArray,ip,el) diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index 3cec1da4c..3d1de3d0a 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -225,7 +225,7 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar3333, ipc, + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**1 & ! linear coefficient + lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**2 & ! quadratic coefficient ) / & - (1.0_pReal \ + (1.0_pReal & + lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**1 / 1. & + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**2 / 2. & + lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**3 / 3. & diff --git a/src/lattice.f90 b/src/lattice.f90 index 996852a79..2b2a5641d 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -72,7 +72,7 @@ module lattice !-------------------------------------------------------------------------------------------------- ! face centered cubic integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: & - LATTICE_fcc_NslipSystem = int([12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],pInt) !< # of slip systems per family for fcc + LATTICE_fcc_NslipSystem = int([12, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],pInt) !< # of slip systems per family for fcc integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: & LATTICE_fcc_NtwinSystem = int([12, 0, 0, 0],pInt) !< # of twin systems per family for fcc @@ -104,11 +104,19 @@ module lattice 1, 1, 0, 1,-1,-1, & ! A6 0, 1, 1, -1, 1,-1, & ! D1 1, 0,-1, -1, 1,-1, & ! D4 - -1,-1, 0, -1, 1,-1 & ! D6 + -1,-1, 0, -1, 1,-1, & ! D6 + ! Slip system <110>{110} + 1, 1, 0, 1,-1, 0, & + 1,-1, 0, 1, 1, 0, & + 1, 0, 1, 1, 0,-1, & + 1, 0,-1, 1, 0, 1, & + 0, 1, 1, 0, 1,-1, & + 0, 1,-1, 0, 1, 1 & ],pReal),shape(LATTICE_FCC_SYSTEMSLIP)) !< Slip system <110>{111} directions. Sorted according to Eisenlohr & Hantcherli - character(len=*), dimension(1), parameter, public :: LATTICE_FCC_SLIPFAMILY_NAME = & - ['<0 1 -1>{1 1 1}'] + character(len=*), dimension(2), parameter, public :: LATTICE_FCC_SLIPFAMILY_NAME = & + ['<0 1 -1>{1 1 1}', & + '<0 1 -1>{0 1 1}'] real(pReal), dimension(3+3,LATTICE_fcc_Ntwin), parameter, private :: & LATTICE_fcc_systemTwin = reshape(real( [& @@ -166,25 +174,38 @@ module lattice integer(pInt), dimension(LATTICE_fcc_Nslip,lattice_fcc_Nslip), parameter, public :: & LATTICE_fcc_interactionSlipSlip = reshape(int( [& - 1,2,2,4,6,5,3,5,5,4,5,6, & ! ---> slip - 2,1,2,6,4,5,5,4,6,5,3,5, & ! | - 2,2,1,5,5,3,5,6,4,6,5,4, & ! | - 4,6,5,1,2,2,4,5,6,3,5,5, & ! v slip - 6,4,5,2,1,2,5,3,5,5,4,6, & - 5,5,3,2,2,1,6,5,4,5,6,4, & - 3,5,5,4,5,6,1,2,2,4,6,5, & - 5,4,6,5,3,5,2,1,2,6,4,5, & - 5,6,4,6,5,4,2,2,1,5,5,3, & - 4,5,6,3,5,5,4,6,5,1,2,2, & - 5,3,5,5,4,6,6,4,5,2,1,2, & - 6,5,4,5,6,4,5,5,3,2,2,1 & - ],pInt),shape(LATTICE_FCC_INTERACTIONSLIPSLIP),order=[2,1]) !< Slip--slip interaction types for fcc + 1, 2, 2, 4, 6, 5, 3, 5, 5, 4, 5, 6, 9,10, 9,10,11,12, & ! ---> slip + 2, 1, 2, 6, 4, 5, 5, 4, 6, 5, 3, 5, 9,10,11,12, 9,10, & ! | + 2, 2, 1, 5, 5, 3, 5, 6, 4, 6, 5, 4, 11,12, 9,10, 9,10, & ! | + 4, 6, 5, 1, 2, 2, 4, 5, 6, 3, 5, 5, 9,10,10, 9,12,11, & ! v slip + 6, 4, 5, 2, 1, 2, 5, 3, 5, 5, 4, 6, 9,10,12,11,10, 9, & + 5, 5, 3, 2, 2, 1, 6, 5, 4, 5, 6, 4, 11,12,10, 9,10, 9, & + 3, 5, 5, 4, 5, 6, 1, 2, 2, 4, 6, 5, 10, 9,10, 9,11,12, & + 5, 4, 6, 5, 3, 5, 2, 1, 2, 6, 4, 5, 10, 9,12,11, 9,10, & + 5, 6, 4, 6, 5, 4, 2, 2, 1, 5, 5, 3, 12,11,10, 9, 9,10, & + 4, 5, 6, 3, 5, 5, 4, 6, 5, 1, 2, 2, 10, 9, 9,10,12,11, & + 5, 3, 5, 5, 4, 6, 6, 4, 5, 2, 1, 2, 10, 9,11,12,10, 9, & + 6, 5, 4, 5, 6, 4, 5, 5, 3, 2, 2, 1, 12,11, 9,10,10, 9, & + + 9, 9,11, 9, 9,11,10,10,12,10,10,12, 1, 7, 8, 8, 8, 8, & + 10,10,12,10,10,12, 9, 9,11, 9, 9,11, 7, 1, 8, 8, 8, 8, & + 9,11, 9,10,12,10,10,12,10, 9,11, 9, 8, 8, 1, 7, 8, 8, & + 10,12,10, 9,11, 9, 9,11, 9,10,12,10, 8, 8, 7, 1, 8, 8, & + 11, 9, 9,12,10,10,11, 9, 9,12,10,10, 8, 8, 8, 8, 1, 7, & + 12,10,10,11, 9, 9,12,10,10,11, 9, 9, 8, 8, 8, 8, 7, 1 & + ],pInt),[LATTICE_fcc_Nslip,LATTICE_fcc_Nslip],order=[2,1]) !< Slip--slip interaction types for fcc !< 1: self interaction !< 2: coplanar interaction !< 3: collinear interaction !< 4: Hirth locks !< 5: glissile junctions !< 6: Lomer locks + !< 7: crossing (similar to Hirth locks in <110>{111} for two {110} planes) + !< 8: similar to Lomer locks in <110>{111} for two {110} planes + !< 9: similar to Lomer locks in <110>{111} btw one {110} and one {111} plane + !<10: similar to glissile junctions in <110>{111} btw one {110} and one {111} plane + !<11: crossing btw one {110} and one {111} plane + !<12: collinear btw one {110} and one {111} plane integer(pInt), dimension(LATTICE_fcc_Nslip,LATTICE_fcc_Ntwin), parameter, public :: & LATTICE_fcc_interactionSlipTwin = reshape(int( [& 1,1,1,3,3,3,2,2,2,3,3,3, & ! ---> twin @@ -198,7 +219,14 @@ module lattice 3,3,3,3,3,3,1,1,1,2,2,2, & 3,3,3,2,2,2,3,3,3,1,1,1, & 2,2,2,3,3,3,3,3,3,1,1,1, & - 3,3,3,3,3,3,2,2,2,1,1,1 & + 3,3,3,3,3,3,2,2,2,1,1,1, & + + 4,4,4,4,4,4,4,4,4,4,4,4, & + 4,4,4,4,4,4,4,4,4,4,4,4, & + 4,4,4,4,4,4,4,4,4,4,4,4, & + 4,4,4,4,4,4,4,4,4,4,4,4, & + 4,4,4,4,4,4,4,4,4,4,4,4, & + 4,4,4,4,4,4,4,4,4,4,4,4 & ],pInt),shape(LATTICE_FCC_INTERACTIONSLIPTWIN),order=[2,1]) !< Slip--twin interaction types for fcc !< 1: coplanar interaction !< 2: screw trace between slip system and twin habit plane (easy cross slip) @@ -235,7 +263,14 @@ module lattice 3,3,3,3,3,3,1,1,1,2,2,2, & 3,3,3,2,2,2,3,3,3,1,1,1, & 2,2,2,3,3,3,3,3,3,1,1,1, & - 3,3,3,3,3,3,2,2,2,1,1,1 & + 3,3,3,3,3,3,2,2,2,1,1,1, & + + 4,4,4,4,4,4,4,4,4,4,4,4, & + 4,4,4,4,4,4,4,4,4,4,4,4, & + 4,4,4,4,4,4,4,4,4,4,4,4, & + 4,4,4,4,4,4,4,4,4,4,4,4, & + 4,4,4,4,4,4,4,4,4,4,4,4, & + 4,4,4,4,4,4,4,4,4,4,4,4 & ],pInt),shape(LATTICE_FCCTOHEX_INTERACTIONSLIPTRANS),order=[2,1]) !< Slip--trans interaction types for fcc integer(pInt), dimension(LATTICE_fcc_Ntrans,LATTICE_fcc_Nslip), parameter, public :: & diff --git a/src/math.f90 b/src/math.f90 index 440ee5303..725c0446e 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -302,7 +302,7 @@ subroutine math_check endif end subroutine math_check - + !-------------------------------------------------------------------------------------------------- !> @brief Quicksort algorithm for two-dimensional integer arrays @@ -2625,12 +2625,9 @@ real(pReal) pure function math_clip(a, left, right) real(pReal), intent(in) :: a real(pReal), intent(in), optional :: left, right - - math_clip = min ( & - max (merge(left, -huge(a), present(left)), a), & - merge(right, huge(a), present(right)) & - ) - + math_clip = a + if (present(left)) math_clip = max(left,math_clip) + if (present(right)) math_clip = min(right,math_clip) if (present(left) .and. present(right)) & math_clip = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_clip, left>right) diff --git a/src/plastic_dislotwin.f90 b/src/plastic_dislotwin.f90 index 2ed8ebfdf..00534d251 100644 --- a/src/plastic_dislotwin.f90 +++ b/src/plastic_dislotwin.f90 @@ -1,3 +1,8 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @author Su Leen Wong, Max-Planck-Institut für Eisenforschung GmbH +!> @author Nan Jia, Max-Planck-Institut für Eisenforschung GmbH +!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief material subroutine incoprorating dislocation and twinning physics !> @details to be done @@ -6,160 +11,147 @@ module plastic_dislotwin use prec, only: & pReal, & pInt - + implicit none private - integer(pInt), dimension(:), allocatable, public, protected :: & - plastic_dislotwin_sizePostResults !< cumulative size of post results - integer(pInt), dimension(:,:), allocatable, target, public :: & plastic_dislotwin_sizePostResult !< size of each post result output - character(len=64), dimension(:,:), allocatable, target, public :: & plastic_dislotwin_output !< name of each post result output - + real(pReal), parameter, private :: & kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin - integer(pInt), dimension(:), allocatable, target, public :: & - plastic_dislotwin_Noutput !< number of outputs per instance of this plasticity - - integer(pInt), dimension(:), allocatable, public, protected :: & - plastic_dislotwin_totalNslip, & !< total number of active slip systems for each instance - plastic_dislotwin_totalNtwin, & !< total number of active twin systems for each instance - plastic_dislotwin_totalNtrans !< number of active transformation systems - - integer(pInt), dimension(:,:), allocatable, private :: & - plastic_dislotwin_Nslip, & !< number of active slip systems for each family and instance - plastic_dislotwin_Ntwin, & !< number of active twin systems for each family and instance - plastic_dislotwin_Ntrans !< number of active transformation systems for each family and instance - - real(pReal), dimension(:), allocatable, private :: & - plastic_dislotwin_CAtomicVolume, & !< atomic volume in Bugers vector unit - plastic_dislotwin_D0, & !< prefactor for self-diffusion coefficient - plastic_dislotwin_Qsd, & !< activation energy for dislocation climb - plastic_dislotwin_GrainSize, & !< grain size - plastic_dislotwin_pShearBand, & !< p-exponent in shearband velocity - plastic_dislotwin_qShearBand, & !< q-exponent in shearband velocity - plastic_dislotwin_MaxTwinFraction, & !< maximum allowed total twin volume fraction - plastic_dislotwin_CEdgeDipMinDistance, & !< - plastic_dislotwin_Cmfptwin, & !< - plastic_dislotwin_Cthresholdtwin, & !< - plastic_dislotwin_SolidSolutionStrength, & !< Strength due to elements in solid solution - plastic_dislotwin_L0_twin, & !< Length of twin nuclei in Burgers vectors - plastic_dislotwin_L0_trans, & !< Length of trans nuclei in Burgers vectors - plastic_dislotwin_xc_twin, & !< critical distance for formation of twin nucleus - plastic_dislotwin_xc_trans, & !< critical distance for formation of trans nucleus - plastic_dislotwin_VcrossSlip, & !< cross slip volume - plastic_dislotwin_sbResistance, & !< value for shearband resistance (might become an internal state variable at some point) - plastic_dislotwin_sbVelocity, & !< value for shearband velocity_0 - plastic_dislotwin_sbQedge, & !< value for shearband systems Qedge - plastic_dislotwin_SFE_0K, & !< stacking fault energy at zero K - plastic_dislotwin_dSFE_dT, & !< temperature dependance of stacking fault energy - plastic_dislotwin_dipoleFormationFactor, & !< scaling factor for dipole formation: 0: off, 1: on. other values not useful - plastic_dislotwin_aTolRho, & !< absolute tolerance for integration of dislocation density - plastic_dislotwin_aTolTwinFrac, & !< absolute tolerance for integration of twin volume fraction - plastic_dislotwin_aTolTransFrac, & !< absolute tolerance for integration of trans volume fraction - plastic_dislotwin_deltaG, & !< Free energy difference between austensite and martensite - plastic_dislotwin_Cmfptrans, & !< - plastic_dislotwin_Cthresholdtrans, & !< - plastic_dislotwin_transStackHeight !< Stack height of hex nucleus - - real(pReal), dimension(:,:,:,:), allocatable, private :: & - plastic_dislotwin_Ctwin66 !< twin elasticity matrix in Mandel notation for each instance - real(pReal), dimension(:,:,:,:,:,:), allocatable, private :: & - plastic_dislotwin_Ctwin3333 !< twin elasticity matrix for each instance - real(pReal), dimension(:,:,:,:), allocatable, private :: & - plastic_dislotwin_Ctrans66 !< trans elasticity matrix in Mandel notation for each instance - real(pReal), dimension(:,:,:,:,:,:), allocatable, private :: & - plastic_dislotwin_Ctrans3333 !< trans elasticity matrix for each instance - real(pReal), dimension(:,:), allocatable, private :: & - plastic_dislotwin_rhoEdge0, & !< initial edge dislocation density per slip system for each family and instance - plastic_dislotwin_rhoEdgeDip0, & !< initial edge dipole density per slip system for each family and instance - plastic_dislotwin_burgersPerSlipFamily, & !< absolute length of burgers vector [m] for each slip family and instance - plastic_dislotwin_burgersPerSlipSystem, & !< absolute length of burgers vector [m] for each slip system and instance - plastic_dislotwin_burgersPerTwinFamily, & !< absolute length of burgers vector [m] for each twin family and instance - plastic_dislotwin_burgersPerTwinSystem, & !< absolute length of burgers vector [m] for each twin system and instance - plastic_dislotwin_burgersPerTransFamily, & !< absolute length of burgers vector [m] for each trans family and instance - plastic_dislotwin_burgersPerTransSystem, & !< absolute length of burgers vector [m] for each trans system and instance - plastic_dislotwin_QedgePerSlipFamily, & !< activation energy for glide [J] for each slip family and instance - plastic_dislotwin_QedgePerSlipSystem, & !< activation energy for glide [J] for each slip system and instance - plastic_dislotwin_v0PerSlipFamily, & !< dislocation velocity prefactor [m/s] for each family and instance - plastic_dislotwin_v0PerSlipSystem, & !< dislocation velocity prefactor [m/s] for each slip system and instance - plastic_dislotwin_tau_peierlsPerSlipFamily, & !< Peierls stress [Pa] for each family and instance - plastic_dislotwin_Ndot0PerTwinFamily, & !< twin nucleation rate [1/m³s] for each twin family and instance - plastic_dislotwin_Ndot0PerTwinSystem, & !< twin nucleation rate [1/m³s] for each twin system and instance - plastic_dislotwin_Ndot0PerTransFamily, & !< trans nucleation rate [1/m³s] for each trans family and instance - plastic_dislotwin_Ndot0PerTransSystem, & !< trans nucleation rate [1/m³s] for each trans system and instance - plastic_dislotwin_tau_r_twin, & !< stress to bring partial close together for each twin system and instance - plastic_dislotwin_tau_r_trans, & !< stress to bring partial close together for each trans system and instance - plastic_dislotwin_twinsizePerTwinFamily, & !< twin thickness [m] for each twin family and instance - plastic_dislotwin_twinsizePerTwinSystem, & !< twin thickness [m] for each twin system and instance - plastic_dislotwin_CLambdaSlipPerSlipFamily, & !< Adj. parameter for distance between 2 forest dislocations for each slip family and instance - plastic_dislotwin_CLambdaSlipPerSlipSystem, & !< Adj. parameter for distance between 2 forest dislocations for each slip system and instance - plastic_dislotwin_lamellarsizePerTransFamily, & !< martensite lamellar thickness [m] for each trans family and instance - plastic_dislotwin_lamellarsizePerTransSystem, & !< martensite lamellar thickness [m] for each trans system and instance - plastic_dislotwin_interaction_SlipSlip, & !< coefficients for slip-slip interaction for each interaction type and instance - plastic_dislotwin_interaction_SlipTwin, & !< coefficients for slip-twin interaction for each interaction type and instance - plastic_dislotwin_interaction_TwinSlip, & !< coefficients for twin-slip interaction for each interaction type and instance - plastic_dislotwin_interaction_TwinTwin, & !< coefficients for twin-twin interaction for each interaction type and instance - plastic_dislotwin_interaction_SlipTrans, & !< coefficients for slip-trans interaction for each interaction type and instance - plastic_dislotwin_interaction_TransSlip, & !< coefficients for trans-slip interaction for each interaction type and instance - plastic_dislotwin_interaction_TransTrans, & !< coefficients for trans-trans interaction for each interaction type and instance - plastic_dislotwin_pPerSlipFamily, & !< p-exponent in glide velocity - plastic_dislotwin_qPerSlipFamily, & !< q-exponent in glide velocity - plastic_dislotwin_rPerTwinFamily, & !< r-exponent in twin nucleation rate - plastic_dislotwin_sPerTransFamily !< s-exponent in trans nucleation rate - real(pReal), dimension(:,:,:), allocatable, private :: & - plastic_dislotwin_interactionMatrix_SlipSlip, & !< interaction matrix of the different slip systems for each instance - plastic_dislotwin_interactionMatrix_SlipTwin, & !< interaction matrix of slip systems with twin systems for each instance - plastic_dislotwin_interactionMatrix_TwinSlip, & !< interaction matrix of twin systems with slip systems for each instance - plastic_dislotwin_interactionMatrix_TwinTwin, & !< interaction matrix of the different twin systems for each instance - plastic_dislotwin_interactionMatrix_SlipTrans, & !< interaction matrix of slip systems with trans systems for each instance - plastic_dislotwin_interactionMatrix_TransSlip, & !< interaction matrix of trans systems with slip systems for each instance - plastic_dislotwin_interactionMatrix_TransTrans, & !< interaction matrix of the different trans systems for each instance - plastic_dislotwin_forestProjectionEdge, & !< matrix of forest projections of edge dislocations for each instance - plastic_dislotwin_projectionMatrix_Trans !< matrix for projection of slip system shear on fault band (twin) systems for each instance - - real(pReal), dimension(:,:,:,:,:), allocatable, private :: & - plastic_dislotwin_sbSv - enum, bind(c) - enumerator :: undefined_ID, & - edge_density_ID, & - dipole_density_ID, & - shear_rate_slip_ID, & - accumulated_shear_slip_ID, & - mfp_slip_ID, & - resolved_stress_slip_ID, & - threshold_stress_slip_ID, & - edge_dipole_distance_ID, & - stress_exponent_ID, & - twin_fraction_ID, & - shear_rate_twin_ID, & - accumulated_shear_twin_ID, & - mfp_twin_ID, & - resolved_stress_twin_ID, & - threshold_stress_twin_ID, & - resolved_stress_shearband_ID, & - shear_rate_shearband_ID, & - sb_eigenvalues_ID, & - sb_eigenvectors_ID, & - stress_trans_fraction_ID, & - strain_trans_fraction_ID, & - trans_fraction_ID + enumerator :: & + undefined_ID, & + edge_density_ID, & + dipole_density_ID, & + shear_rate_slip_ID, & + accumulated_shear_slip_ID, & + mfp_slip_ID, & + resolved_stress_slip_ID, & + threshold_stress_slip_ID, & + edge_dipole_distance_ID, & + stress_exponent_ID, & + twin_fraction_ID, & + shear_rate_twin_ID, & + accumulated_shear_twin_ID, & + mfp_twin_ID, & + resolved_stress_twin_ID, & + threshold_stress_twin_ID, & + resolved_stress_shearband_ID, & + shear_rate_shearband_ID, & + stress_trans_fraction_ID, & + strain_trans_fraction_ID end enum - integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & - plastic_dislotwin_outputID !< ID of each post result output + + type, private :: tParameters + real(pReal) :: & + mu, & + nu, & + CAtomicVolume, & !< atomic volume in Bugers vector unit + D0, & !< prefactor for self-diffusion coefficient + Qsd, & !< activation energy for dislocation climb + GrainSize, & !>>' write(6,'(/,a)') ' A. Ma and F. Roters, Acta Materialia, 52(12):3603–3612, 2004' - write(6,'(/,a)') ' https://doi.org/10.1016/j.actamat.2004.04.012' + write(6,'(a)') ' https://doi.org/10.1016/j.actamat.2004.04.012' write(6,'(/,a)') ' F.Roters et al., Computational Materials Science, 39:91–95, 2007' - write(6,'(/,a)') ' https://doi.org/10.1016/j.commatsci.2006.04.014' + write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2006.04.014' write(6,'(/,a)') ' Wong et al., Acta Materialia, 118:140–151, 2016' - write(6,'(/,a)') ' https://doi.org/10.1016/j.actamat.2016.07.032' + write(6,'(a,/)') ' https://doi.org/10.1016/j.actamat.2016.07.032' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - maxNinstance = int(count(phase_plasticity == PLASTICITY_DISLOTWIN_ID),pInt) - if (maxNinstance == 0_pInt) return + Ninstance = int(count(phase_plasticity == PLASTICITY_DISLOTWIN_ID),pInt) + if (Ninstance == 0_pInt) return if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - allocate(plastic_dislotwin_sizePostResults(maxNinstance), source=0_pInt) - allocate(plastic_dislotwin_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) - allocate(plastic_dislotwin_output(maxval(phase_Noutput),maxNinstance)) + + allocate(plastic_dislotwin_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt) + allocate(plastic_dislotwin_output(maxval(phase_Noutput),Ninstance)) plastic_dislotwin_output = '' - allocate(plastic_dislotwin_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID) - allocate(plastic_dislotwin_Noutput(maxNinstance), source=0_pInt) - allocate(plastic_dislotwin_Nslip(lattice_maxNslipFamily,maxNinstance), source=0_pInt) - allocate(plastic_dislotwin_Ntwin(lattice_maxNtwinFamily,maxNinstance), source=0_pInt) - allocate(plastic_dislotwin_Ntrans(lattice_maxNtransFamily,maxNinstance), source=0_pInt) - allocate(plastic_dislotwin_totalNslip(maxNinstance), source=0_pInt) - allocate(plastic_dislotwin_totalNtwin(maxNinstance), source=0_pInt) - allocate(plastic_dislotwin_totalNtrans(maxNinstance), source=0_pInt) - allocate(plastic_dislotwin_CAtomicVolume(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_D0(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_Qsd(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_GrainSize(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_pShearBand(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_qShearBand(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_MaxTwinFraction(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_CEdgeDipMinDistance(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_Cmfptwin(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_Cthresholdtwin(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_SolidSolutionStrength(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_L0_twin(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_L0_trans(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_xc_twin(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_xc_trans(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_VcrossSlip(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_aTolRho(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_aTolTwinFrac(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_aTolTransFrac(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_sbResistance(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_sbVelocity(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_sbQedge(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_SFE_0K(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_dSFE_dT(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_dipoleFormationFactor(maxNinstance), source=1.0_pReal) !should be on by default - allocate(plastic_dislotwin_deltaG(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_Cmfptrans(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_Cthresholdtrans(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_transStackHeight(maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_rhoEdge0(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_rhoEdgeDip0(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_burgersPerSlipFamily(lattice_maxNslipFamily,maxNinstance), & - source=0.0_pReal) - allocate(plastic_dislotwin_burgersPerTwinFamily(lattice_maxNtwinFamily,maxNinstance), & - source=0.0_pReal) - allocate(plastic_dislotwin_burgersPerTransFamily(lattice_maxNtransFamily,maxNinstance), & - source=0.0_pReal) - allocate(plastic_dislotwin_QedgePerSlipFamily(lattice_maxNslipFamily,maxNinstance), & - source=0.0_pReal) - allocate(plastic_dislotwin_v0PerSlipFamily(lattice_maxNslipFamily,maxNinstance), & - source=0.0_pReal) - allocate(plastic_dislotwin_tau_peierlsPerSlipFamily(lattice_maxNslipFamily,maxNinstance), & - source=0.0_pReal) - allocate(plastic_dislotwin_pPerSlipFamily(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal) - allocate(plastic_dislotwin_qPerSlipFamily(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal) - allocate(plastic_dislotwin_Ndot0PerTwinFamily(lattice_maxNtwinFamily,maxNinstance), & - source=0.0_pReal) - allocate(plastic_dislotwin_Ndot0PerTransFamily(lattice_maxNtransFamily,maxNinstance), & - source=0.0_pReal) - allocate(plastic_dislotwin_twinsizePerTwinFamily(lattice_maxNtwinFamily,maxNinstance), & - source=0.0_pReal) - allocate(plastic_dislotwin_CLambdaSlipPerSlipFamily(lattice_maxNslipFamily,maxNinstance), & - source=0.0_pReal) - allocate(plastic_dislotwin_rPerTwinFamily(lattice_maxNtwinFamily,maxNinstance),source=0.0_pReal) - allocate(plastic_dislotwin_interaction_SlipSlip(lattice_maxNinteraction,maxNinstance), & - source=0.0_pReal) - allocate(plastic_dislotwin_interaction_SlipTwin(lattice_maxNinteraction,maxNinstance), & - source=0.0_pReal) - allocate(plastic_dislotwin_interaction_TwinSlip(lattice_maxNinteraction,maxNinstance), & - source=0.0_pReal) - allocate(plastic_dislotwin_interaction_TwinTwin(lattice_maxNinteraction,maxNinstance), & - source=0.0_pReal) - allocate(plastic_dislotwin_interaction_SlipTrans(lattice_maxNinteraction,maxNinstance), & - source=0.0_pReal) - allocate(plastic_dislotwin_interaction_TransSlip(lattice_maxNinteraction,maxNinstance), & - source=0.0_pReal) - allocate(plastic_dislotwin_interaction_TransTrans(lattice_maxNinteraction,maxNinstance), & - source=0.0_pReal) - allocate(plastic_dislotwin_sbSv(6,6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), & - source=0.0_pReal) - allocate(plastic_dislotwin_lamellarsizePerTransFamily(lattice_maxNtransFamily,maxNinstance), & - source=0.0_pReal) - allocate(plastic_dislotwin_sPerTransFamily(lattice_maxNtransFamily,maxNinstance),source=0.0_pReal) + + allocate(param(Ninstance)) + allocate(state(Ninstance)) + allocate(dotState(Ninstance)) + allocate(microstructure(Ninstance)) + + do p = 1_pInt, size(phase_plasticityInstance) + if (phase_plasticity(p) /= PLASTICITY_DISLOTWIN_ID) cycle + associate(prm => param(phase_plasticityInstance(p)), & + dot => dotState(phase_plasticityInstance(p)), & + stt => state(phase_plasticityInstance(p)), & + mse => microstructure(phase_plasticityInstance(p))) + + ! This data is read in already in lattice + prm%isFCC = merge(.true., .false., lattice_structure(p) == LATTICE_FCC_ID) + prm%mu = lattice_mu(p) + prm%nu = lattice_nu(p) + prm%C66 = lattice_C66(1:6,1:6,p) + + structure = config_phase(p)%getString('lattice_structure') - rewind(fileUnit) - phase = 0_pInt - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to - line = IO_read(fileUnit) - enddo - - parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part - line = IO_read(fileUnit) - if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') then ! stop at next part - line = IO_read(fileUnit, .true.) ! reset IO_read - exit - endif - if (IO_getTag(line,'[',']') /= '') then ! next phase section - phase = phase + 1_pInt ! advance phase section counter - if (phase_plasticity(phase) == PLASTICITY_DISLOTWIN_ID) then - Nchunks_SlipFamilies = count(lattice_NslipSystem(:,phase) > 0_pInt) - Nchunks_TwinFamilies = count(lattice_NtwinSystem(:,phase) > 0_pInt) - Nchunks_TransFamilies = count(lattice_NtransSystem(:,phase)> 0_pInt) - Nchunks_SlipSlip = maxval(lattice_interactionSlipSlip(:,:,phase)) - Nchunks_SlipTwin = maxval(lattice_interactionSlipTwin(:,:,phase)) - Nchunks_TwinSlip = maxval(lattice_interactionTwinSlip(:,:,phase)) - Nchunks_TwinTwin = maxval(lattice_interactionTwinTwin(:,:,phase)) - Nchunks_SlipTrans = maxval(lattice_interactionSlipTrans(:,:,phase)) - Nchunks_TransSlip = maxval(lattice_interactionTransSlip(:,:,phase)) - Nchunks_TransTrans = maxval(lattice_interactionTransTrans(:,:,phase)) - if(allocated(tempPerSlip)) deallocate(tempPerSlip) - if(allocated(tempPerTwin)) deallocate(tempPerTwin) - if(allocated(tempPerTrans)) deallocate(tempPerTrans) - allocate(tempPerSlip(Nchunks_SlipFamilies)) - allocate(tempPerTwin(Nchunks_TwinFamilies)) - allocate(tempPerTrans(Nchunks_TransFamilies)) +!-------------------------------------------------------------------------------------------------- +! slip related parameters + prm%Nslip = config_phase(p)%getInts('nslip',defaultVal=emptyIntArray) + prm%totalNslip = sum(prm%Nslip) + slipActive: if (prm%totalNslip > 0_pInt) then + prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),& + config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal)) + prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, & + config_phase(p)%getFloats('interaction_slipslip'), & + structure(1:3)) + + prm%rho0 = config_phase(p)%getFloats('rhoedge0', requiredShape=shape(prm%Nslip)) !ToDo: rename to rho_0 + prm%rhoDip0 = config_phase(p)%getFloats('rhoedgedip0',requiredShape=shape(prm%Nslip)) !ToDo: rename to rho_dip_0 + prm%v0 = config_phase(p)%getFloats('v0', requiredShape=shape(prm%Nslip)) + prm%burgers_slip = config_phase(p)%getFloats('slipburgers',requiredShape=shape(prm%Nslip)) + prm%Qedge = config_phase(p)%getFloats('qedge', requiredShape=shape(prm%Nslip)) !ToDo: rename (ask Karo) + prm%CLambdaSlip = config_phase(p)%getFloats('clambdaslip',requiredShape=shape(prm%Nslip)) + prm%p = config_phase(p)%getFloats('p_slip', requiredShape=shape(prm%Nslip)) + prm%q = config_phase(p)%getFloats('q_slip', requiredShape=shape(prm%Nslip)) + prm%B = config_phase(p)%getFloats('b', requiredShape=shape(prm%Nslip), & + defaultVal=[(0.0_pReal, i=1,size(prm%Nslip))]) + prm%tau_peierls = config_phase(p)%getFloats('tau_peierls',requiredShape=shape(prm%Nslip), & + defaultVal=[(0.0_pReal, i=1,size(prm%Nslip))]) + + prm%CEdgeDipMinDistance = config_phase(p)%getFloat('cedgedipmindistance') + + ! expand: family => system + prm%rho0 = math_expand(prm%rho0, prm%Nslip) + prm%rhoDip0 = math_expand(prm%rhoDip0, prm%Nslip) + prm%v0 = math_expand(prm%v0, prm%Nslip) + prm%burgers_slip = math_expand(prm%burgers_slip,prm%Nslip) + prm%Qedge = math_expand(prm%Qedge, prm%Nslip) + prm%CLambdaSlip = math_expand(prm%CLambdaSlip, prm%Nslip) + prm%p = math_expand(prm%p, prm%Nslip) + prm%q = math_expand(prm%q, prm%Nslip) + prm%B = math_expand(prm%B, prm%Nslip) + prm%tau_peierls = math_expand(prm%tau_peierls, prm%Nslip) + + ! sanity checks + if (any(prm%rho0 < 0.0_pReal)) extmsg = trim(extmsg)//'rho0 ' + if (any(prm%rhoDip0 < 0.0_pReal)) extmsg = trim(extmsg)//'rhoDip0 ' + if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//'v0 ' + if (any(prm%burgers_slip <= 0.0_pReal)) extmsg = trim(extmsg)//'burgers_slip ' + if (any(prm%Qedge <= 0.0_pReal)) extmsg = trim(extmsg)//'Qedge ' + if (any(prm%CLambdaSlip <= 0.0_pReal)) extmsg = trim(extmsg)//'CLambdaSlip ' + if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//'B ' + if (any(prm%tau_peierls < 0.0_pReal)) extmsg = trim(extmsg)//'tau_peierls ' + if (any(prm%p<=0.0_pReal .or. prm%p>1.0_pReal)) extmsg = trim(extmsg)//'p ' + if (any(prm%q< 1.0_pReal .or. prm%q>2.0_pReal)) extmsg = trim(extmsg)//'q ' + + else slipActive + allocate(prm%burgers_slip(0)) + endif slipActive + +!-------------------------------------------------------------------------------------------------- +! twin related parameters + prm%Ntwin = config_phase(p)%getInts('ntwin', defaultVal=emptyIntArray) + prm%totalNtwin = sum(prm%Ntwin) + if (prm%totalNtwin > 0_pInt) then + prm%Schmid_twin = lattice_SchmidMatrix_twin(prm%Ntwin,structure(1:3),& + config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal)) + prm%interaction_TwinTwin = lattice_interaction_TwinTwin(prm%Ntwin,& + config_phase(p)%getFloats('interaction_twintwin'), & + structure(1:3)) + + prm%burgers_twin = config_phase(p)%getFloats('twinburgers') + prm%twinsize = config_phase(p)%getFloats('twinsize') + prm%r = config_phase(p)%getFloats('r_twin') + + prm%xc_twin = config_phase(p)%getFloat('xc_twin') + prm%L0_twin = config_phase(p)%getFloat('l0_twin') + prm%MaxTwinFraction = config_phase(p)%getFloat('maxtwinfraction') ! ToDo: only used in postResults + prm%Cthresholdtwin = config_phase(p)%getFloat('cthresholdtwin', defaultVal=0.0_pReal) + prm%Cmfptwin = config_phase(p)%getFloat('cmfptwin', defaultVal=0.0_pReal) ! ToDo: How to handle that??? + + + if (.not. prm%isFCC) then + prm%Ndot0_twin = config_phase(p)%getFloats('ndot0_twin') + prm%Ndot0_twin = math_expand(prm%Ndot0_twin,prm%Ntwin) endif - cycle ! skip to next line + + ! expand: family => system + prm%burgers_twin = math_expand(prm%burgers_twin,prm%Ntwin) + prm%twinsize = math_expand(prm%twinsize,prm%Ntwin) + prm%r = math_expand(prm%r,prm%Ntwin) + + else + allocate(prm%twinsize(0)) + allocate(prm%burgers_twin(0)) + allocate(prm%r(0)) + endif + +!-------------------------------------------------------------------------------------------------- +! transformation related parameters + prm%Ntrans = config_phase(p)%getInts('ntrans', defaultVal=emptyIntArray) + prm%totalNtrans = sum(prm%Ntrans) + if (prm%totalNtrans > 0_pInt) then + prm%burgers_trans = config_phase(p)%getFloats('transburgers') + prm%burgers_trans = math_expand(prm%burgers_trans,prm%Ntrans) + + prm%Cthresholdtrans = config_phase(p)%getFloat('cthresholdtrans', defaultVal=0.0_pReal) ! ToDo: How to handle that??? + prm%transStackHeight = config_phase(p)%getFloat('transstackheight', defaultVal=0.0_pReal) ! ToDo: How to handle that??? + prm%Cmfptrans = config_phase(p)%getFloat('cmfptrans', defaultVal=0.0_pReal) ! ToDo: How to handle that??? + prm%deltaG = config_phase(p)%getFloat('deltag') + prm%xc_trans = config_phase(p)%getFloat('xc_trans', defaultVal=0.0_pReal) ! ToDo: How to handle that??? + prm%L0_trans = config_phase(p)%getFloat('l0_trans') + + prm%interaction_TransTrans = spread(config_phase(p)%getFloats('interaction_transtrans'),2,1) + if (lattice_structure(p) /= LATTICE_fcc_ID) then + prm%Ndot0_trans = config_phase(p)%getFloats('ndot0_trans') + prm%Ndot0_trans = math_expand(prm%Ndot0_trans,prm%Ntrans) + endif + prm%lamellarsizePerTransSystem = config_phase(p)%getFloats('lamellarsize') + prm%lamellarsizePerTransSystem = math_expand(prm%lamellarsizePerTransSystem,prm%Ntrans) + prm%s = config_phase(p)%getFloats('s_trans',defaultVal=[0.0_pReal]) + prm%s = math_expand(prm%s,prm%Ntrans) + else + allocate(prm%lamellarsizePerTransSystem(0)) + allocate(prm%burgers_trans(0)) + endif + + if (sum(prm%Ntwin) > 0_pInt .or. prm%totalNtrans > 0_pInt) then + prm%SFE_0K = config_phase(p)%getFloat('sfe_0k') + prm%dSFE_dT = config_phase(p)%getFloat('dsfe_dt') + prm%VcrossSlip = config_phase(p)%getFloat('vcrossslip') + endif + + if (prm%totalNslip > 0_pInt .and. prm%totalNtwin > 0_pInt) then + prm%interaction_SlipTwin = lattice_interaction_SlipTwin(prm%Nslip,prm%Ntwin,& + config_phase(p)%getFloats('interaction_sliptwin'), & + structure(1:3)) + prm%interaction_TwinSlip = lattice_interaction_TwinSlip(prm%Ntwin,prm%Nslip,& + config_phase(p)%getFloats('interaction_twinslip'), & + structure(1:3)) + endif + + if (prm%totalNslip > 0_pInt .and. prm%totalNtrans > 0_pInt) then + prm%interaction_TransSlip = spread(config_phase(p)%getFloats('interaction_transslip'),2,1) + prm%interaction_SlipTrans = spread(config_phase(p)%getFloats('interaction_sliptrans'),2,1) + endif + + + prm%aTolRho = config_phase(p)%getFloat('atol_rho', defaultVal=0.0_pReal) + prm%aTolTwinFrac = config_phase(p)%getFloat('atol_twinfrac', defaultVal=0.0_pReal) + prm%aTolTransFrac = config_phase(p)%getFloat('atol_transfrac', defaultVal=0.0_pReal) + + prm%CAtomicVolume = config_phase(p)%getFloat('catomicvolume') + prm%GrainSize = config_phase(p)%getFloat('grainsize') + + + prm%D0 = config_phase(p)%getFloat('d0') + prm%Qsd = config_phase(p)%getFloat('qsd') + prm%SolidSolutionStrength = config_phase(p)%getFloat('solidsolutionstrength') + if (config_phase(p)%keyExists('dipoleformationfactor')) call IO_error(1,ext_msg='use /nodipoleformation/') + prm%dipoleformation = .not. config_phase(p)%keyExists('/nodipoleformation/') + prm%sbVelocity = config_phase(p)%getFloat('shearbandvelocity',defaultVal=0.0_pReal) + if (prm%sbVelocity > 0.0_pReal) then + prm%sbResistance = config_phase(p)%getFloat('shearbandresistance') + prm%sbQedge = config_phase(p)%getFloat('qedgepersbsystem') + prm%pShearBand = config_phase(p)%getFloat('p_shearband') + prm%qShearBand = config_phase(p)%getFloat('q_shearband') endif - if (phase > 0_pInt ) then; if (phase_plasticity(phase) == PLASTICITY_DISLOTWIN_ID) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran - instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase - chunkPos = IO_stringPos(line) - tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key - select case(tag) - case ('(output)') - select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt))) - case ('edge_density') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = edge_density_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('dipole_density') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = dipole_density_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('shear_rate_slip','shearrate_slip') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = shear_rate_slip_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('accumulated_shear_slip') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = accumulated_shear_slip_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('mfp_slip') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = mfp_slip_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('resolved_stress_slip') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = resolved_stress_slip_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('threshold_stress_slip') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = threshold_stress_slip_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('edge_dipole_distance') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = edge_dipole_distance_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('stress_exponent') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = stress_exponent_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('twin_fraction') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = twin_fraction_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('shear_rate_twin','shearrate_twin') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = shear_rate_twin_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('accumulated_shear_twin') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = accumulated_shear_twin_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('mfp_twin') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = mfp_twin_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('resolved_stress_twin') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = resolved_stress_twin_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('threshold_stress_twin') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = threshold_stress_twin_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('resolved_stress_shearband') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = resolved_stress_shearband_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('shear_rate_shearband','shearrate_shearband') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = shear_rate_shearband_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('sb_eigenvalues') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = sb_eigenvalues_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('sb_eigenvectors') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = sb_eigenvectors_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('stress_trans_fraction') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = stress_trans_fraction_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('strain_trans_fraction') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = strain_trans_fraction_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case ('trans_fraction','total_trans_fraction') - plastic_dislotwin_Noutput(instance) = plastic_dislotwin_Noutput(instance) + 1_pInt - plastic_dislotwin_outputID(plastic_dislotwin_Noutput(instance),instance) = trans_fraction_ID - plastic_dislotwin_output(plastic_dislotwin_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - end select -!-------------------------------------------------------------------------------------------------- -! parameters depending on number of slip system families - case ('nslip') - if (chunkPos(1) < Nchunks_SlipFamilies + 1_pInt) & - call IO_warning(50_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOTWIN_label//')') - if (chunkPos(1) > Nchunks_SlipFamilies + 1_pInt) & - call IO_error(150_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOTWIN_label//')') - Nchunks_SlipFamilies = chunkPos(1) - 1_pInt - do j = 1_pInt, Nchunks_SlipFamilies - plastic_dislotwin_Nslip(j,instance) = IO_intValue(line,chunkPos,1_pInt+j) - enddo - case ('rhoedge0','rhoedgedip0','slipburgers','qedge','v0','clambdaslip','tau_peierls','p_slip','q_slip') - do j = 1_pInt, Nchunks_SlipFamilies - tempPerSlip(j) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo - select case(tag) - case ('rhoedge0') - plastic_dislotwin_rhoEdge0(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) - case ('rhoedgedip0') - plastic_dislotwin_rhoEdgeDip0(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) - case ('slipburgers') - plastic_dislotwin_burgersPerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) - case ('qedge') - plastic_dislotwin_QedgePerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) - case ('v0') - plastic_dislotwin_v0PerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) - case ('clambdaslip') - plastic_dislotwin_CLambdaSlipPerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) - case ('tau_peierls') - if (lattice_structure(phase) /= LATTICE_bcc_ID) & - call IO_warning(42_pInt,ext_msg=trim(tag)//' for non-bcc ('//PLASTICITY_DISLOTWIN_label//')') - plastic_dislotwin_tau_peierlsPerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) - case ('p_slip') - plastic_dislotwin_pPerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) - case ('q_slip') - plastic_dislotwin_qPerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) - end select -!-------------------------------------------------------------------------------------------------- -! parameters depending on slip number of twin families - case ('ntwin') - if (chunkPos(1) < Nchunks_TwinFamilies + 1_pInt) & - call IO_warning(51_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOTWIN_label//')') - if (chunkPos(1) > Nchunks_TwinFamilies + 1_pInt) & - call IO_error(150_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOTWIN_label//')') - Nchunks_TwinFamilies = chunkPos(1) - 1_pInt - do j = 1_pInt, Nchunks_TwinFamilies - plastic_dislotwin_Ntwin(j,instance) = IO_intValue(line,chunkPos,1_pInt+j) - enddo - case ('ndot0_twin','twinsize','twinburgers','r_twin') - do j = 1_pInt, Nchunks_TwinFamilies - tempPerTwin(j) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo - select case(tag) - case ('ndot0_twin') - if (lattice_structure(phase) == LATTICE_fcc_ID) & - call IO_warning(42_pInt,ext_msg=trim(tag)//' for fcc ('//PLASTICITY_DISLOTWIN_label//')') - plastic_dislotwin_Ndot0PerTwinFamily(1:Nchunks_TwinFamilies,instance) = tempPerTwin(1:Nchunks_TwinFamilies) - case ('twinsize') - plastic_dislotwin_twinsizePerTwinFamily(1:Nchunks_TwinFamilies,instance) = tempPerTwin(1:Nchunks_TwinFamilies) - case ('twinburgers') - plastic_dislotwin_burgersPerTwinFamily(1:Nchunks_TwinFamilies,instance) = tempPerTwin(1:Nchunks_TwinFamilies) - case ('r_twin') - plastic_dislotwin_rPerTwinFamily(1:Nchunks_TwinFamilies,instance) = tempPerTwin(1:Nchunks_TwinFamilies) - end select -!-------------------------------------------------------------------------------------------------- -! parameters depending on number of transformation system families - case ('ntrans') - if (chunkPos(1) < Nchunks_TransFamilies + 1_pInt) & - call IO_warning(53_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOTWIN_label//')') - if (chunkPos(1) > Nchunks_TransFamilies + 1_pInt) & - call IO_error(150_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOTWIN_label//')') - Nchunks_TransFamilies = chunkPos(1) - 1_pInt - do j = 1_pInt, Nchunks_TransFamilies - plastic_dislotwin_Ntrans(j,instance) = IO_intValue(line,chunkPos,1_pInt+j) - enddo - case ('ndot0_trans','lamellarsize','transburgers','s_trans') - do j = 1_pInt, Nchunks_TransFamilies - tempPerTrans(j) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo - select case(tag) - case ('ndot0_trans') - if (lattice_structure(phase) == LATTICE_fcc_ID) & - call IO_warning(42_pInt,ext_msg=trim(tag)//' for fcc ('//PLASTICITY_DISLOTWIN_label//')') - plastic_dislotwin_Ndot0PerTransFamily(1:Nchunks_TransFamilies,instance) = tempPerTrans(1:Nchunks_TransFamilies) - case ('lamellarsize') - plastic_dislotwin_lamellarsizePerTransFamily(1:Nchunks_TransFamilies,instance) = tempPerTrans(1:Nchunks_TransFamilies) - case ('transburgers') - plastic_dislotwin_burgersPerTransFamily(1:Nchunks_TransFamilies,instance) = tempPerTrans(1:Nchunks_TransFamilies) - case ('s_trans') - plastic_dislotwin_sPerTransFamily(1:Nchunks_TransFamilies,instance) = tempPerTrans(1:Nchunks_TransFamilies) - end select -!-------------------------------------------------------------------------------------------------- -! parameters depending on number of interactions - case ('interaction_slipslip','interactionslipslip') - if (chunkPos(1) < 1_pInt + Nchunks_SlipSlip) & - call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOTWIN_label//')') - do j = 1_pInt, Nchunks_SlipSlip - plastic_dislotwin_interaction_SlipSlip(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo - case ('interaction_sliptwin','interactionsliptwin') - if (chunkPos(1) < 1_pInt + Nchunks_SlipTwin) & - call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOTWIN_label//')') - do j = 1_pInt, Nchunks_SlipTwin - plastic_dislotwin_interaction_SlipTwin(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo - case ('interaction_twinslip','interactiontwinslip') - if (chunkPos(1) < 1_pInt + Nchunks_TwinSlip) & - call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOTWIN_label//')') - do j = 1_pInt, Nchunks_TwinSlip - plastic_dislotwin_interaction_TwinSlip(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo - case ('interaction_twintwin','interactiontwintwin') - if (chunkPos(1) < 1_pInt + Nchunks_TwinTwin) & - call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOTWIN_label//')') - do j = 1_pInt, Nchunks_TwinTwin - plastic_dislotwin_interaction_TwinTwin(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo - case ('interaction_sliptrans','interactionsliptrans') - if (chunkPos(1) < 1_pInt + Nchunks_SlipTrans) & - call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOTWIN_label//')') - do j = 1_pInt, Nchunks_SlipTrans - plastic_dislotwin_interaction_SlipTrans(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo - case ('interaction_transslip','interactiontransslip') - if (chunkPos(1) < 1_pInt + Nchunks_TransSlip) & - call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOTWIN_label//')') - do j = 1_pInt, Nchunks_TransSlip - plastic_dislotwin_interaction_TransSlip(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo - case ('interaction_transtrans','interactiontranstrans') - if (chunkPos(1) < 1_pInt + Nchunks_TransTrans) & - call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOTWIN_label//')') - do j = 1_pInt, Nchunks_TransTrans - plastic_dislotwin_interaction_TransTrans(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo -!-------------------------------------------------------------------------------------------------- -! parameters independent of number of slip/twin/trans systems - case ('grainsize') - plastic_dislotwin_GrainSize(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('maxtwinfraction') - plastic_dislotwin_MaxTwinFraction(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('p_shearband') - plastic_dislotwin_pShearBand(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('q_shearband') - plastic_dislotwin_qShearBand(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('d0') - plastic_dislotwin_D0(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('qsd') - plastic_dislotwin_Qsd(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('atol_rho') - plastic_dislotwin_aTolRho(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('atol_twinfrac') - plastic_dislotwin_aTolTwinFrac(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('atol_transfrac') - plastic_dislotwin_aTolTransFrac(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('cmfptwin') - plastic_dislotwin_Cmfptwin(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('cthresholdtwin') - plastic_dislotwin_Cthresholdtwin(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('solidsolutionstrength') - plastic_dislotwin_SolidSolutionStrength(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('l0_twin') - plastic_dislotwin_L0_twin(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('l0_trans') - plastic_dislotwin_L0_trans(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('xc_twin') - plastic_dislotwin_xc_twin(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('xc_trans') - plastic_dislotwin_xc_trans(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('vcrossslip') - plastic_dislotwin_VcrossSlip(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('cedgedipmindistance') - plastic_dislotwin_CEdgeDipMinDistance(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('catomicvolume') - plastic_dislotwin_CAtomicVolume(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('sfe_0k') - plastic_dislotwin_SFE_0K(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('dsfe_dt') - plastic_dislotwin_dSFE_dT(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('dipoleformationfactor') - plastic_dislotwin_dipoleFormationFactor(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('shearbandresistance') - plastic_dislotwin_sbResistance(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('shearbandvelocity') - plastic_dislotwin_sbVelocity(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('qedgepersbsystem') - plastic_dislotwin_sbQedge(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('deltag') - plastic_dislotwin_deltaG(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('cmfptrans') - plastic_dislotwin_Cmfptrans(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('cthresholdtrans') - plastic_dislotwin_Cthresholdtrans(instance) = IO_floatValue(line,chunkPos,2_pInt) - case ('transstackheight') - plastic_dislotwin_transStackHeight(instance) = IO_floatValue(line,chunkPos,2_pInt) + !if (Ndot0PerTwinFamily(f,p) < 0.0_pReal) & + ! call IO_error(211_pInt,el=p,ext_msg='ndot0_twin ('//PLASTICITY_DISLOTWIN_label//')') + + if (prm%CAtomicVolume <= 0.0_pReal) & + call IO_error(211_pInt,el=p,ext_msg='cAtomicVolume ('//PLASTICITY_DISLOTWIN_label//')') + if (prm%D0 <= 0.0_pReal) & + call IO_error(211_pInt,el=p,ext_msg='D0 ('//PLASTICITY_DISLOTWIN_label//')') + if (prm%Qsd <= 0.0_pReal) & + call IO_error(211_pInt,el=p,ext_msg='Qsd ('//PLASTICITY_DISLOTWIN_label//')') + if (prm%totalNtwin > 0_pInt) then + if (dEq0(prm%SFE_0K) .and. & + dEq0(prm%dSFE_dT) .and. & + lattice_structure(p) == LATTICE_fcc_ID) & + call IO_error(211_pInt,el=p,ext_msg='SFE0K ('//PLASTICITY_DISLOTWIN_label//')') + if (prm%aTolRho <= 0.0_pReal) & + call IO_error(211_pInt,el=p,ext_msg='aTolRho ('//PLASTICITY_DISLOTWIN_label//')') + if (prm%aTolTwinFrac <= 0.0_pReal) & + call IO_error(211_pInt,el=p,ext_msg='aTolTwinFrac ('//PLASTICITY_DISLOTWIN_label//')') + endif + if (prm%totalNtrans > 0_pInt) then + if (dEq0(prm%SFE_0K) .and. & + dEq0(prm%dSFE_dT) .and. & + lattice_structure(p) == LATTICE_fcc_ID) & + call IO_error(211_pInt,el=p,ext_msg='SFE0K ('//PLASTICITY_DISLOTWIN_label//')') + if (prm%aTolTransFrac <= 0.0_pReal) & + call IO_error(211_pInt,el=p,ext_msg='aTolTransFrac ('//PLASTICITY_DISLOTWIN_label//')') + endif + !if (prm%sbResistance < 0.0_pReal) & + ! call IO_error(211_pInt,el=p,ext_msg='sbResistance ('//PLASTICITY_DISLOTWIN_label//')') + !if (prm%sbVelocity < 0.0_pReal) & + ! call IO_error(211_pInt,el=p,ext_msg='sbVelocity ('//PLASTICITY_DISLOTWIN_label//')') + !if (prm%sbVelocity > 0.0_pReal .and. & + ! prm%pShearBand <= 0.0_pReal) & + ! call IO_error(211_pInt,el=p,ext_msg='pShearBand ('//PLASTICITY_DISLOTWIN_label//')') + if (prm%sbVelocity > 0.0_pReal .and. & + prm%qShearBand <= 0.0_pReal) & + call IO_error(211_pInt,el=p,ext_msg='qShearBand ('//PLASTICITY_DISLOTWIN_label//')') + + outputs = config_phase(p)%getStrings('(output)', defaultVal=emptyStringArray) + allocate(prm%outputID(0)) + do i= 1_pInt, size(outputs) + outputID = undefined_ID + select case(outputs(i)) + case ('edge_density') + outputID = merge(edge_density_ID,undefined_ID,prm%totalNslip > 0_pInt) + outputSize = prm%totalNslip + case ('dipole_density') + outputID = merge(dipole_density_ID,undefined_ID,prm%totalNslip > 0_pInt) + outputSize = prm%totalNslip + case ('shear_rate_slip','shearrate_slip') + outputID = merge(shear_rate_slip_ID,undefined_ID,prm%totalNslip > 0_pInt) + outputSize = prm%totalNslip + case ('accumulated_shear_slip') + outputID = merge(accumulated_shear_slip_ID,undefined_ID,prm%totalNslip > 0_pInt) + outputSize = prm%totalNslip + case ('mfp_slip') + outputID = merge(mfp_slip_ID,undefined_ID,prm%totalNslip > 0_pInt) + outputSize = prm%totalNslip + case ('resolved_stress_slip') + outputID = merge(resolved_stress_slip_ID,undefined_ID,prm%totalNslip > 0_pInt) + outputSize = prm%totalNslip + case ('threshold_stress_slip') + outputID= merge(threshold_stress_slip_ID,undefined_ID,prm%totalNslip > 0_pInt) + outputSize = prm%totalNslip + case ('edge_dipole_distance') + outputID = merge(edge_dipole_distance_ID,undefined_ID,prm%totalNslip > 0_pInt) + outputSize = prm%totalNslip + case ('stress_exponent') + outputID = merge(stress_exponent_ID,undefined_ID,prm%totalNslip > 0_pInt) + outputSize = prm%totalNslip + + case ('twin_fraction') + outputID = merge(twin_fraction_ID,undefined_ID,prm%totalNtwin >0_pInt) + outputSize = prm%totalNtwin + case ('shear_rate_twin','shearrate_twin') + outputID = merge(shear_rate_twin_ID,undefined_ID,prm%totalNtwin >0_pInt) + outputSize = prm%totalNtwin + case ('accumulated_shear_twin') + outputID = merge(accumulated_shear_twin_ID,undefined_ID,prm%totalNtwin >0_pInt) + outputSize = prm%totalNtwin + case ('mfp_twin') + outputID = merge(mfp_twin_ID,undefined_ID,prm%totalNtwin >0_pInt) + outputSize = prm%totalNtwin + case ('resolved_stress_twin') + outputID = merge(resolved_stress_twin_ID,undefined_ID,prm%totalNtwin >0_pInt) + outputSize = prm%totalNtwin + case ('threshold_stress_twin') + outputID = merge(threshold_stress_twin_ID,undefined_ID,prm%totalNtwin >0_pInt) + outputSize = prm%totalNtwin + + case ('resolved_stress_shearband') + outputID = resolved_stress_shearband_ID + outputSize = 6_pInt + case ('shear_rate_shearband','shearrate_shearband') + outputID = shear_rate_shearband_ID + outputSize = 6_pInt + + case ('stress_trans_fraction') + outputID = stress_trans_fraction_ID + outputSize = prm%totalNtrans + case ('strain_trans_fraction') + outputID = strain_trans_fraction_ID + outputSize = prm%totalNtrans + end select - endif; endif - enddo parsingFile - - sanityChecks: do phase = 1_pInt, size(phase_plasticity) - myPhase: if (phase_plasticity(phase) == PLASTICITY_dislotwin_ID) then - instance = phase_plasticityInstance(phase) - - if (sum(plastic_dislotwin_Nslip(:,instance)) < 0_pInt) & - call IO_error(211_pInt,el=instance,ext_msg='Nslip ('//PLASTICITY_DISLOTWIN_label//')') - if (sum(plastic_dislotwin_Ntwin(:,instance)) < 0_pInt) & - call IO_error(211_pInt,el=instance,ext_msg='Ntwin ('//PLASTICITY_DISLOTWIN_label//')') - if (sum(plastic_dislotwin_Ntrans(:,instance)) < 0_pInt) & - call IO_error(211_pInt,el=instance,ext_msg='Ntrans ('//PLASTICITY_DISLOTWIN_label//')') - do f = 1_pInt,lattice_maxNslipFamily - if (plastic_dislotwin_Nslip(f,instance) > 0_pInt) then - if (plastic_dislotwin_rhoEdge0(f,instance) < 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='rhoEdge0 ('//PLASTICITY_DISLOTWIN_label//')') - if (plastic_dislotwin_rhoEdgeDip0(f,instance) < 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='rhoEdgeDip0 ('//PLASTICITY_DISLOTWIN_label//')') - if (plastic_dislotwin_burgersPerSlipFamily(f,instance) <= 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='slipBurgers ('//PLASTICITY_DISLOTWIN_label//')') - if (plastic_dislotwin_v0PerSlipFamily(f,instance) <= 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='v0 ('//PLASTICITY_DISLOTWIN_label//')') - if (plastic_dislotwin_tau_peierlsPerSlipFamily(f,instance) < 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='tau_peierls ('//PLASTICITY_DISLOTWIN_label//')') - endif - enddo - do f = 1_pInt,lattice_maxNtwinFamily - if (plastic_dislotwin_Ntwin(f,instance) > 0_pInt) then - if (plastic_dislotwin_burgersPerTwinFamily(f,instance) <= 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='twinburgers ('//PLASTICITY_DISLOTWIN_label//')') - if (plastic_dislotwin_Ndot0PerTwinFamily(f,instance) < 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='ndot0_twin ('//PLASTICITY_DISLOTWIN_label//')') - endif - enddo - if (plastic_dislotwin_CAtomicVolume(instance) <= 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='cAtomicVolume ('//PLASTICITY_DISLOTWIN_label//')') - if (plastic_dislotwin_D0(instance) <= 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='D0 ('//PLASTICITY_DISLOTWIN_label//')') - if (plastic_dislotwin_Qsd(instance) <= 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='Qsd ('//PLASTICITY_DISLOTWIN_label//')') - if (sum(plastic_dislotwin_Ntwin(:,instance)) > 0_pInt) then - if (dEq0(plastic_dislotwin_SFE_0K(instance)) .and. & - dEq0(plastic_dislotwin_dSFE_dT(instance)) .and. & - lattice_structure(phase) == LATTICE_fcc_ID) & - call IO_error(211_pInt,el=instance,ext_msg='SFE0K ('//PLASTICITY_DISLOTWIN_label//')') - if (plastic_dislotwin_aTolRho(instance) <= 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='aTolRho ('//PLASTICITY_DISLOTWIN_label//')') - if (plastic_dislotwin_aTolTwinFrac(instance) <= 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='aTolTwinFrac ('//PLASTICITY_DISLOTWIN_label//')') - endif - if (sum(plastic_dislotwin_Ntrans(:,instance)) > 0_pInt) then - if (dEq0(plastic_dislotwin_SFE_0K(instance)) .and. & - dEq0(plastic_dislotwin_dSFE_dT(instance)) .and. & - lattice_structure(phase) == LATTICE_fcc_ID) & - call IO_error(211_pInt,el=instance,ext_msg='SFE0K ('//PLASTICITY_DISLOTWIN_label//')') - if (plastic_dislotwin_aTolTransFrac(instance) <= 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='aTolTransFrac ('//PLASTICITY_DISLOTWIN_label//')') - endif - if (plastic_dislotwin_sbResistance(instance) < 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='sbResistance ('//PLASTICITY_DISLOTWIN_label//')') - if (plastic_dislotwin_sbVelocity(instance) < 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='sbVelocity ('//PLASTICITY_DISLOTWIN_label//')') - if (plastic_dislotwin_sbVelocity(instance) > 0.0_pReal .and. & - plastic_dislotwin_pShearBand(instance) <= 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='pShearBand ('//PLASTICITY_DISLOTWIN_label//')') - if (dNeq0(plastic_dislotwin_dipoleFormationFactor(instance)) .and. & - dNeq(plastic_dislotwin_dipoleFormationFactor(instance), 1.0_pReal)) & - call IO_error(211_pInt,el=instance,ext_msg='dipoleFormationFactor ('//PLASTICITY_DISLOTWIN_label//')') - if (plastic_dislotwin_sbVelocity(instance) > 0.0_pReal .and. & - plastic_dislotwin_qShearBand(instance) <= 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='qShearBand ('//PLASTICITY_DISLOTWIN_label//')') - -!-------------------------------------------------------------------------------------------------- -! Determine total number of active slip or twin systems - plastic_dislotwin_Nslip(:,instance) = min(lattice_NslipSystem(:,phase),plastic_dislotwin_Nslip(:,instance)) - plastic_dislotwin_Ntwin(:,instance) = min(lattice_NtwinSystem(:,phase),plastic_dislotwin_Ntwin(:,instance)) - plastic_dislotwin_Ntrans(:,instance)= min(lattice_NtransSystem(:,phase),plastic_dislotwin_Ntrans(:,instance)) - plastic_dislotwin_totalNslip(instance) = sum(plastic_dislotwin_Nslip(:,instance)) - plastic_dislotwin_totalNtwin(instance) = sum(plastic_dislotwin_Ntwin(:,instance)) - plastic_dislotwin_totalNtrans(instance) = sum(plastic_dislotwin_Ntrans(:,instance)) - endif myPhase - enddo sanityChecks - -!-------------------------------------------------------------------------------------------------- -! allocation of variables whose size depends on the total number of active slip systems - maxTotalNslip = maxval(plastic_dislotwin_totalNslip) - maxTotalNtwin = maxval(plastic_dislotwin_totalNtwin) - maxTotalNtrans = maxval(plastic_dislotwin_totalNtrans) - - allocate(plastic_dislotwin_burgersPerSlipSystem(maxTotalNslip, maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_burgersPerTwinSystem(maxTotalNtwin, maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_burgersPerTransSystem(maxTotalNtrans, maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_QedgePerSlipSystem(maxTotalNslip, maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_v0PerSlipSystem(maxTotalNslip, maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_Ndot0PerTwinSystem(maxTotalNtwin, maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_Ndot0PerTransSystem(maxTotalNtrans, maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_tau_r_twin(maxTotalNtwin, maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_tau_r_trans(maxTotalNtrans, maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_twinsizePerTwinSystem(maxTotalNtwin, maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_CLambdaSlipPerSlipSystem(maxTotalNslip, maxNinstance),source=0.0_pReal) - allocate(plastic_dislotwin_lamellarsizePerTransSystem(maxTotalNtrans, maxNinstance),source=0.0_pReal) - - allocate(plastic_dislotwin_interactionMatrix_SlipSlip(maxval(plastic_dislotwin_totalNslip),& ! slip resistance from slip activity - maxval(plastic_dislotwin_totalNslip),& - maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_interactionMatrix_SlipTwin(maxval(plastic_dislotwin_totalNslip),& ! slip resistance from twin activity - maxval(plastic_dislotwin_totalNtwin),& - maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_interactionMatrix_TwinSlip(maxval(plastic_dislotwin_totalNtwin),& ! twin resistance from slip activity - maxval(plastic_dislotwin_totalNslip),& - maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_interactionMatrix_TwinTwin(maxval(plastic_dislotwin_totalNtwin),& ! twin resistance from twin activity - maxval(plastic_dislotwin_totalNtwin),& - maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_interactionMatrix_SlipTrans(maxval(plastic_dislotwin_totalNslip),& ! slip resistance from trans activity - maxval(plastic_dislotwin_totalNtrans),& - maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_interactionMatrix_TransSlip(maxval(plastic_dislotwin_totalNtrans),& ! trans resistance from slip activity - maxval(plastic_dislotwin_totalNslip),& - maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_interactionMatrix_TransTrans(maxval(plastic_dislotwin_totalNtrans),& ! trans resistance from trans activity - maxval(plastic_dislotwin_totalNtrans),& - maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_forestProjectionEdge(maxTotalNslip,maxTotalNslip,maxNinstance), & - source=0.0_pReal) - allocate(plastic_dislotwin_projectionMatrix_Trans(maxTotalNtrans,maxTotalNslip,maxNinstance), & - source=0.0_pReal) - allocate(plastic_dislotwin_Ctwin66(6,6,maxTotalNtwin,maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_Ctwin3333(3,3,3,3,maxTotalNtwin,maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_Ctrans66(6,6,maxTotalNtrans,maxNinstance), source=0.0_pReal) - allocate(plastic_dislotwin_Ctrans3333(3,3,3,3,maxTotalNtrans,maxNinstance), source=0.0_pReal) - - allocate(state(maxNinstance)) - allocate(state0(maxNinstance)) - allocate(dotState(maxNinstance)) - - initializeInstances: do phase = 1_pInt, size(phase_plasticity) - myPhase2: if (phase_plasticity(phase) == PLASTICITY_dislotwin_ID) then - NofMyPhase=count(material_phase==phase) - instance = phase_plasticityInstance(phase) - - ns = plastic_dislotwin_totalNslip(instance) - nt = plastic_dislotwin_totalNtwin(instance) - nr = plastic_dislotwin_totalNtrans(instance) - -!-------------------------------------------------------------------------------------------------- -! Determine size of postResults array - outputsLoop: do o = 1_pInt,plastic_dislotwin_Noutput(instance) - select case(plastic_dislotwin_outputID(o,instance)) - case(edge_density_ID, & - dipole_density_ID, & - shear_rate_slip_ID, & - accumulated_shear_slip_ID, & - mfp_slip_ID, & - resolved_stress_slip_ID, & - threshold_stress_slip_ID, & - edge_dipole_distance_ID, & - stress_exponent_ID & - ) - mySize = ns - case(twin_fraction_ID, & - shear_rate_twin_ID, & - accumulated_shear_twin_ID, & - mfp_twin_ID, & - resolved_stress_twin_ID, & - threshold_stress_twin_ID & - ) - mySize = nt - case(resolved_stress_shearband_ID, & - shear_rate_shearband_ID & - ) - mySize = 6_pInt - case(sb_eigenvalues_ID) - mySize = 3_pInt - case(sb_eigenvectors_ID) - mySize = 9_pInt - case(stress_trans_fraction_ID, & - strain_trans_fraction_ID, & - trans_fraction_ID & - ) - mySize = nr - end select - - if (mySize > 0_pInt) then ! any meaningful output found - plastic_dislotwin_sizePostResult(o,instance) = mySize - plastic_dislotwin_sizePostResults(instance) = plastic_dislotwin_sizePostResults(instance) + mySize - endif - enddo outputsLoop + + if (outputID /= undefined_ID) then + plastic_dislotwin_output(i,phase_plasticityInstance(p)) = outputs(i) + plastic_dislotwin_sizePostResult(i,phase_plasticityInstance(p)) = outputSize + prm%outputID = [prm%outputID, outputID] + endif + enddo + !-------------------------------------------------------------------------------------------------- ! allocate state arrays + NipcMyPhase=count(material_phase==p) + sizeDotState = int(size(['rho ','rhoDip ','accshearslip']),pInt) * prm%totalNslip & + + int(size(['twinFraction','accsheartwin']),pInt) * prm%totalNtwin & + + int(size(['stressTransFraction','strainTransFraction']),pInt) * prm%totalNtrans + sizeState = sizeDotState - sizeDotState = int(size(['rhoEdge ','rhoEdgeDip ','accshearslip']),pInt) * ns & - + int(size(['twinFraction','accsheartwin']),pInt) * nt & - + int(size(['stressTransFraction','strainTransFraction']),pInt) * nr - sizeDeltaState = 0_pInt - sizeState = sizeDotState & - + int(size(['invLambdaSlip ','invLambdaSlipTwin ','invLambdaSlipTrans',& - 'meanFreePathSlip ','tauSlipThreshold ']),pInt) * ns & - + int(size(['invLambdaTwin ','meanFreePathTwin','tauTwinThreshold',& - 'twinVolume ']),pInt) * nt & - + int(size(['invLambdaTrans ','meanFreePathTrans','tauTransThreshold', & - 'martensiteVolume ']),pInt) * nr - - plasticState(phase)%sizeState = sizeState - plasticState(phase)%sizeDotState = sizeDotState - plasticState(phase)%sizeDeltaState = sizeDeltaState - plasticState(phase)%sizePostResults = plastic_dislotwin_sizePostResults(instance) - plasticState(phase)%nSlip = plastic_dislotwin_totalNslip(instance) - plasticState(phase)%nTwin = plastic_dislotwin_totalNtwin(instance) - plasticState(phase)%nTrans= plastic_dislotwin_totalNtrans(instance) - allocate(plasticState(phase)%aTolState (sizeState), source=0.0_pReal) - allocate(plasticState(phase)%state0 (sizeState,NofMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%subState0 (sizeState,NofMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%state (sizeState,NofMyPhase), source=0.0_pReal) + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0_pInt, & + prm%totalNslip,prm%totalNtwin,prm%totalNtrans) + plasticState(p)%sizePostResults = sum(plastic_dislotwin_sizePostResult(:,phase_plasticityInstance(p))) - allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) - if (any(numerics_integrator == 1_pInt)) then - allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) - endif - if (any(numerics_integrator == 4_pInt)) & - allocate(plasticState(phase)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal) - if (any(numerics_integrator == 5_pInt)) & - allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal) - offset_slip = 2_pInt*plasticState(phase)%nslip - plasticState(phase)%slipRate => & - plasticState(phase)%dotState(offset_slip+1:offset_slip+plasticState(phase)%nslip,1:NofMyPhase) - plasticState(phase)%accumulatedSlip => & - plasticState(phase)%state (offset_slip+1:offset_slip+plasticState(phase)%nslip,1:NofMyPhase) + ! ToDo: do later on + offset_slip = 2_pInt*plasticState(p)%nslip + plasticState(p)%slipRate => & + plasticState(p)%dotState(offset_slip+1:offset_slip+plasticState(p)%nslip,1:NipcMyPhase) + plasticState(p)%accumulatedSlip => & + plasticState(p)%state (offset_slip+1:offset_slip+plasticState(p)%nslip,1:NipcMyPhase) + + allocate(temp1(prm%totalNslip,prm%totalNtrans),source =0.0_pReal) + allocate(prm%forestProjectionEdge(prm%totalNslip,prm%totalNslip),source = 0.0_pReal) + i = 0_pInt + mySlipFamilies: do f = 1_pInt,size(prm%Nslip,1) + index_myFamily = sum(prm%Nslip(1:f-1_pInt)) - !* Process slip related parameters ------------------------------------------------ - slipFamiliesLoop: do f = 1_pInt,lattice_maxNslipFamily - index_myFamily = sum(plastic_dislotwin_Nslip(1:f-1_pInt,instance)) ! index in truncated slip system list - slipSystemsLoop: do j = 1_pInt,plastic_dislotwin_Nslip(f,instance) - - !* Burgers vector, - ! dislocation velocity prefactor, - ! mean free path prefactor, - ! and minimum dipole distance - - plastic_dislotwin_burgersPerSlipSystem(index_myFamily+j,instance) = & - plastic_dislotwin_burgersPerSlipFamily(f,instance) - - plastic_dislotwin_QedgePerSlipSystem(index_myFamily+j,instance) = & - plastic_dislotwin_QedgePerSlipFamily(f,instance) - - plastic_dislotwin_v0PerSlipSystem(index_myFamily+j,instance) = & - plastic_dislotwin_v0PerSlipFamily(f,instance) - - plastic_dislotwin_CLambdaSlipPerSlipSystem(index_myFamily+j,instance) = & - plastic_dislotwin_CLambdaSlipPerSlipFamily(f,instance) + slipSystemsLoop: do j = 1_pInt,prm%Nslip(f) + i = i + 1_pInt + do o = 1_pInt, size(prm%Nslip,1) + index_otherFamily = sum(prm%Nslip(1:o-1_pInt)) + do k = 1_pInt,prm%Nslip(o) ! loop over (active) systems in other family (slip) + prm%forestProjectionEdge(index_myFamily+j,index_otherFamily+k) = & + abs(math_mul3x3(lattice_sn(:,sum(lattice_NslipSystem(1:f-1,p))+j,p), & + lattice_st(:,sum(lattice_NslipSystem(1:o-1,p))+k,p))) + enddo; enddo + do o = 1_pInt,size(prm%Ntrans,1) + index_otherFamily = sum(prm%Ntrans(1:o-1_pInt)) + do k = 1_pInt,prm%Ntrans(o) ! loop over (active) systems in other family (trans) + temp1(index_myFamily+j,index_otherFamily+k) = & + prm%interaction_SlipTrans(lattice_interactionSlipTrans( & + sum(lattice_NslipSystem(1:f-1_pInt,p))+j, & + sum(lattice_NtransSystem(1:o-1_pInt,p))+k, & + p),1 ) + enddo; enddo - !* Calculation of forest projections for edge dislocations - !* Interaction matrices - do o = 1_pInt,lattice_maxNslipFamily - index_otherFamily = sum(plastic_dislotwin_Nslip(1:o-1_pInt,instance)) - do k = 1_pInt,plastic_dislotwin_Nslip(o,instance) ! loop over (active) systems in other family (slip) - plastic_dislotwin_forestProjectionEdge(index_myFamily+j,index_otherFamily+k,instance) = & - abs(math_mul3x3(lattice_sn(:,sum(lattice_NslipSystem(1:f-1,phase))+j,phase), & - lattice_st(:,sum(lattice_NslipSystem(1:o-1,phase))+k,phase))) - plastic_dislotwin_interactionMatrix_SlipSlip(index_myFamily+j,index_otherFamily+k,instance) = & - plastic_dislotwin_interaction_SlipSlip(lattice_interactionSlipSlip( & - sum(lattice_NslipSystem(1:f-1,phase))+j, & - sum(lattice_NslipSystem(1:o-1,phase))+k, & - phase), instance ) - enddo; enddo - - do o = 1_pInt,lattice_maxNtwinFamily - index_otherFamily = sum(plastic_dislotwin_Ntwin(1:o-1_pInt,instance)) - do k = 1_pInt,plastic_dislotwin_Ntwin(o,instance) ! loop over (active) systems in other family (twin) - plastic_dislotwin_interactionMatrix_SlipTwin(index_myFamily+j,index_otherFamily+k,instance) = & - plastic_dislotwin_interaction_SlipTwin(lattice_interactionSlipTwin( & - sum(lattice_NslipSystem(1:f-1_pInt,phase))+j, & - sum(lattice_NtwinSystem(1:o-1_pInt,phase))+k, & - phase), instance ) - enddo; enddo + enddo slipSystemsLoop + enddo mySlipFamilies + prm%interaction_SlipTrans = temp1; deallocate(temp1) + + allocate(prm%C66_twin(6,6,prm%totalNtwin), source=0.0_pReal) + if (lattice_structure(p) == LATTICE_fcc_ID) & + allocate(prm%fcc_twinNucleationSlipPair(2,prm%totalNtwin),source = 0_pInt) + allocate(prm%shear_twin(prm%totalNtwin),source = 0.0_pReal) + i = 0_pInt + twinFamiliesLoop: do f = 1_pInt, size(prm%Ntwin,1) + index_myFamily = sum(prm%Ntwin(1:f-1_pInt)) ! index in truncated twin system list + twinSystemsLoop: do j = 1_pInt,prm%Ntwin(f) + i = i + 1_pInt + prm%shear_twin(i) = lattice_shearTwin(sum(lattice_Ntwinsystem(1:f-1,p))+j,p) + if (lattice_structure(p) == LATTICE_fcc_ID) prm%fcc_twinNucleationSlipPair(1:2,i) = & + lattice_fcc_twinNucleationSlipPair(1:2,sum(lattice_Ntwinsystem(1:f-1,p))+j) + !* Rotate twin elasticity matrices + index_otherFamily = sum(lattice_NtwinSystem(1:f-1_pInt,p)) ! index in full lattice twin list + prm%C66_twin(1:6,1:6,index_myFamily+j) = & + math_Mandel3333to66(math_rotate_forward3333(lattice_C3333(1:3,1:3,1:3,1:3,p),& + lattice_Qtwin(1:3,1:3,index_otherFamily+j,p))) + enddo twinSystemsLoop + enddo twinFamiliesLoop - do o = 1_pInt,lattice_maxNtransFamily - index_otherFamily = sum(plastic_dislotwin_Ntrans(1:o-1_pInt,instance)) - do k = 1_pInt,plastic_dislotwin_Ntrans(o,instance) ! loop over (active) systems in other family (trans) - plastic_dislotwin_interactionMatrix_SlipTrans(index_myFamily+j,index_otherFamily+k,instance) = & - plastic_dislotwin_interaction_SlipTrans(lattice_interactionSlipTrans( & - sum(lattice_NslipSystem(1:f-1_pInt,phase))+j, & - sum(lattice_NtransSystem(1:o-1_pInt,phase))+k, & - phase), instance ) - enddo; enddo - - enddo slipSystemsLoop - enddo slipFamiliesLoop - - !* Process twin related parameters ------------------------------------------------ - twinFamiliesLoop: do f = 1_pInt,lattice_maxNtwinFamily - index_myFamily = sum(plastic_dislotwin_Ntwin(1:f-1_pInt,instance)) ! index in truncated twin system list - twinSystemsLoop: do j = 1_pInt,plastic_dislotwin_Ntwin(f,instance) - !* Burgers vector, - ! nucleation rate prefactor, - ! and twin size - - plastic_dislotwin_burgersPerTwinSystem(index_myFamily+j,instance) = & - plastic_dislotwin_burgersPerTwinFamily(f,instance) + allocate(temp1(prm%totalNtrans,prm%totalNslip), source =0.0_pReal) + allocate(temp2(prm%totalNtrans,prm%totalNtrans), source =0.0_pReal) + allocate(prm%C66_trans(6,6,prm%totalNtrans) ,source=0.0_pReal) + allocate(prm%Schmid_trans(3,3,prm%totalNtrans),source = 0.0_pReal) + i = 0_pInt + transFamiliesLoop: do f = 1_pInt,size(prm%Ntrans,1) + index_myFamily = sum(prm%Ntrans(1:f-1_pInt)) ! index in truncated trans system list + transSystemsLoop: do j = 1_pInt,prm%Ntrans(f) + i = i + 1_pInt + prm%Schmid_trans(1:3,1:3,i) = lattice_Strans(1:3,1:3,sum(lattice_Ntranssystem(1:f-1,p))+j,p) + !* Rotate trans elasticity matrices + index_otherFamily = sum(lattice_NtransSystem(1:f-1_pInt,p)) ! index in full lattice trans list + prm%C66_trans(1:6,1:6,index_myFamily+j) = & + math_Mandel3333to66(math_rotate_forward3333(lattice_trans_C3333(1:3,1:3,1:3,1:3,p),& + lattice_Qtrans(1:3,1:3,index_otherFamily+j,p))) + !* Interaction matrices + do o = 1_pInt,size(prm%Nslip,1) + index_otherFamily = sum(prm%Nslip(1:o-1_pInt)) + do k = 1_pInt,prm%Nslip(o) ! loop over (active) systems in other family (slip) + temp1(index_myFamily+j,index_otherFamily+k) = & + prm%interaction_TransSlip(lattice_interactionTransSlip( & + sum(lattice_NtransSystem(1:f-1_pInt,p))+j, & + sum(lattice_NslipSystem(1:o-1_pInt,p))+k, & + p) ,1 ) + enddo; enddo - plastic_dislotwin_Ndot0PerTwinSystem(index_myFamily+j,instance) = & - plastic_dislotwin_Ndot0PerTwinFamily(f,instance) + do o = 1_pInt,size(prm%Ntrans,1) + index_otherFamily = sum(prm%Ntrans(1:o-1_pInt)) + do k = 1_pInt,prm%Ntrans(o) ! loop over (active) systems in other family (trans) + temp2(index_myFamily+j,index_otherFamily+k) = & + prm%interaction_TransTrans(lattice_interactionTransTrans( & + sum(lattice_NtransSystem(1:f-1_pInt,p))+j, & + sum(lattice_NtransSystem(1:o-1_pInt,p))+k, & + p),1 ) + enddo; enddo - plastic_dislotwin_twinsizePerTwinSystem(index_myFamily+j,instance) = & - plastic_dislotwin_twinsizePerTwinFamily(f,instance) - - !* Rotate twin elasticity matrices - index_otherFamily = sum(lattice_NtwinSystem(1:f-1_pInt,phase)) ! index in full lattice twin list - do l = 1_pInt,3_pInt; do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt; do o = 1_pInt,3_pInt - do p = 1_pInt,3_pInt; do q = 1_pInt,3_pInt; do r = 1_pInt,3_pInt; do s = 1_pInt,3_pInt - plastic_dislotwin_Ctwin3333(l,m,n,o,index_myFamily+j,instance) = & - plastic_dislotwin_Ctwin3333(l,m,n,o,index_myFamily+j,instance) + & - lattice_C3333(p,q,r,s,phase) * & - lattice_Qtwin(l,p,index_otherFamily+j,phase) * & - lattice_Qtwin(m,q,index_otherFamily+j,phase) * & - lattice_Qtwin(n,r,index_otherFamily+j,phase) * & - lattice_Qtwin(o,s,index_otherFamily+j,phase) - enddo; enddo; enddo; enddo - enddo; enddo; enddo; enddo - plastic_dislotwin_Ctwin66(1:6,1:6,index_myFamily+j,instance) = & - math_Mandel3333to66(plastic_dislotwin_Ctwin3333(1:3,1:3,1:3,1:3,index_myFamily+j,instance)) - - !* Interaction matrices - do o = 1_pInt,lattice_maxNslipFamily - index_otherFamily = sum(plastic_dislotwin_Nslip(1:o-1_pInt,instance)) - do k = 1_pInt,plastic_dislotwin_Nslip(o,instance) ! loop over (active) systems in other family (slip) - plastic_dislotwin_interactionMatrix_TwinSlip(index_myFamily+j,index_otherFamily+k,instance) = & - plastic_dislotwin_interaction_TwinSlip(lattice_interactionTwinSlip( & - sum(lattice_NtwinSystem(1:f-1_pInt,phase))+j, & - sum(lattice_NslipSystem(1:o-1_pInt,phase))+k, & - phase), instance ) - enddo; enddo - - do o = 1_pInt,lattice_maxNtwinFamily - index_otherFamily = sum(plastic_dislotwin_Ntwin(1:o-1_pInt,instance)) - do k = 1_pInt,plastic_dislotwin_Ntwin(o,instance) ! loop over (active) systems in other family (twin) - plastic_dislotwin_interactionMatrix_TwinTwin(index_myFamily+j,index_otherFamily+k,instance) = & - plastic_dislotwin_interaction_TwinTwin(lattice_interactionTwinTwin( & - sum(lattice_NtwinSystem(1:f-1_pInt,phase))+j, & - sum(lattice_NtwinSystem(1:o-1_pInt,phase))+k, & - phase), instance ) - enddo; enddo - - enddo twinSystemsLoop - enddo twinFamiliesLoop + enddo transSystemsLoop + enddo transFamiliesLoop + prm%interaction_TransSlip = temp1; deallocate(temp1) + prm%interaction_TransTrans = temp2; deallocate(temp2) - !* Process transformation related parameters ------------------------------------------------ - transFamiliesLoop: do f = 1_pInt,lattice_maxNtransFamily - index_myFamily = sum(plastic_dislotwin_Ntrans(1:f-1_pInt,instance)) ! index in truncated trans system list - transSystemsLoop: do j = 1_pInt,plastic_dislotwin_Ntrans(f,instance) + startIndex=1_pInt + endIndex=prm%totalNslip + stt%rhoEdge=>plasticState(p)%state(startIndex:endIndex,:) + stt%rhoEdge= spread(prm%rho0,2,NipcMyPhase) + dot%rhoEdge=>plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolRho - !* Burgers vector, - ! nucleation rate prefactor, - ! and martensite size + startIndex=endIndex+1 + endIndex=endIndex+prm%totalNslip + stt%rhoEdgeDip=>plasticState(p)%state(startIndex:endIndex,:) + stt%rhoEdgeDip= spread(prm%rhoDip0,2,NipcMyPhase) + dot%rhoEdgeDip=>plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolRho + + startIndex=endIndex+1 + endIndex=endIndex+prm%totalNslip + stt%accshear_slip=>plasticState(p)%state(startIndex:endIndex,:) + dot%accshear_slip=>plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = 1.0e6_pReal + + startIndex=endIndex+1 + endIndex=endIndex+prm%totalNtwin + stt%twinFraction=>plasticState(p)%state(startIndex:endIndex,:) + dot%twinFraction=>plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolTwinFrac + + startIndex=endIndex+1 + endIndex=endIndex+prm%totalNtwin + stt%accshear_twin=>plasticState(p)%state(startIndex:endIndex,:) + dot%accshear_twin=>plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = 1.0e6_pReal + + startIndex=endIndex+1 + endIndex=endIndex+prm%totalNtrans + stt%stressTransFraction=>plasticState(p)%state(startIndex:endIndex,:) + dot%stressTransFraction=>plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolTransFrac + + startIndex=endIndex+1 + endIndex=endIndex+prm%totalNtrans + stt%strainTransFraction=>plasticState(p)%state(startIndex:endIndex,:) + dot%strainTransFraction=>plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolTransFrac - plastic_dislotwin_burgersPerTransSystem(index_myFamily+j,instance) = & - plastic_dislotwin_burgersPerTransFamily(f,instance) + plasticState(p)%state0 = plasticState(p)%state + dot%whole => plasticState(p)%dotState - plastic_dislotwin_Ndot0PerTransSystem(index_myFamily+j,instance) = & - plastic_dislotwin_Ndot0PerTransFamily(f,instance) + allocate(mse%invLambdaSlip(prm%totalNslip,NipcMyPhase),source=0.0_pReal) + allocate(mse%invLambdaSlipTwin(prm%totalNslip,NipcMyPhase),source=0.0_pReal) + allocate(mse%invLambdaTwin(prm%totalNtwin,NipcMyPhase),source=0.0_pReal) + allocate(mse%invLambdaSlipTrans(prm%totalNtrans,NipcMyPhase),source=0.0_pReal) + allocate(mse%invLambdaTrans(prm%totalNtrans,NipcMyPhase),source=0.0_pReal) - plastic_dislotwin_lamellarsizePerTransSystem(index_myFamily+j,instance) = & - plastic_dislotwin_lamellarsizePerTransFamily(f,instance) + allocate(mse%mfp_slip(prm%totalNslip,NipcMyPhase), source=0.0_pReal) + allocate(mse%mfp_twin(prm%totalNtwin,NipcMyPhase), source=0.0_pReal) + allocate(mse%mfp_trans(prm%totalNtrans,NipcMyPhase),source=0.0_pReal) - !* Rotate trans elasticity matrices - index_otherFamily = sum(lattice_NtransSystem(1:f-1_pInt,phase)) ! index in full lattice trans list - do l = 1_pInt,3_pInt; do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt; do o = 1_pInt,3_pInt - do p = 1_pInt,3_pInt; do q = 1_pInt,3_pInt; do r = 1_pInt,3_pInt; do s = 1_pInt,3_pInt - plastic_dislotwin_Ctrans3333(l,m,n,o,index_myFamily+j,instance) = & - plastic_dislotwin_Ctrans3333(l,m,n,o,index_myFamily+j,instance) + & - lattice_trans_C3333(p,q,r,s,phase) * & - lattice_Qtrans(l,p,index_otherFamily+j,phase) * & - lattice_Qtrans(m,q,index_otherFamily+j,phase) * & - lattice_Qtrans(n,r,index_otherFamily+j,phase) * & - lattice_Qtrans(o,s,index_otherFamily+j,phase) - enddo; enddo; enddo; enddo - enddo; enddo; enddo; enddo - plastic_dislotwin_Ctrans66(1:6,1:6,index_myFamily+j,instance) = & - math_Mandel3333to66(plastic_dislotwin_Ctrans3333(1:3,1:3,1:3,1:3,index_myFamily+j,instance)) + allocate(mse%threshold_stress_slip(prm%totalNslip,NipcMyPhase), source=0.0_pReal) + allocate(mse%threshold_stress_twin(prm%totalNtwin,NipcMyPhase), source=0.0_pReal) + allocate(mse%threshold_stress_trans(prm%totalNtrans,NipcMyPhase),source=0.0_pReal) - !* Interaction matrices - do o = 1_pInt,lattice_maxNslipFamily - index_otherFamily = sum(plastic_dislotwin_Nslip(1:o-1_pInt,instance)) - do k = 1_pInt,plastic_dislotwin_Nslip(o,instance) ! loop over (active) systems in other family (slip) - plastic_dislotwin_interactionMatrix_TransSlip(index_myFamily+j,index_otherFamily+k,instance) = & - plastic_dislotwin_interaction_TransSlip(lattice_interactionTransSlip( & - sum(lattice_NtransSystem(1:f-1_pInt,phase))+j, & - sum(lattice_NslipSystem(1:o-1_pInt,phase))+k, & - phase), instance ) - enddo; enddo + allocate(mse%tau_r_twin(prm%totalNtwin,NipcMyPhase), source=0.0_pReal) + allocate(mse%tau_r_trans(prm%totalNtrans,NipcMyPhase), source=0.0_pReal) - do o = 1_pInt,lattice_maxNtransFamily - index_otherFamily = sum(plastic_dislotwin_Ntrans(1:o-1_pInt,instance)) - do k = 1_pInt,plastic_dislotwin_Ntrans(o,instance) ! loop over (active) systems in other family (trans) - plastic_dislotwin_interactionMatrix_TransTrans(index_myFamily+j,index_otherFamily+k,instance) = & - plastic_dislotwin_interaction_TransTrans(lattice_interactionTransTrans( & - sum(lattice_NtransSystem(1:f-1_pInt,phase))+j, & - sum(lattice_NtransSystem(1:o-1_pInt,phase))+k, & - phase), instance ) - enddo; enddo + allocate(mse%twinVolume(prm%totalNtwin,NipcMyPhase), source=0.0_pReal) + allocate(mse%martensiteVolume(prm%totalNtrans,NipcMyPhase), source=0.0_pReal) - !* Projection matrices for shear from slip systems to fault-band (twin) systems for strain-induced martensite nucleation - select case(trans_lattice_structure(phase)) - case (LATTICE_bcc_ID) - do o = 1_pInt,lattice_maxNtransFamily - index_otherFamily = sum(plastic_dislotwin_Nslip(1:o-1_pInt,instance)) - do k = 1_pInt,plastic_dislotwin_Nslip(o,instance) ! loop over (active) systems in other family (trans) - plastic_dislotwin_projectionMatrix_Trans(index_myFamily+j,index_otherFamily+k,instance) = & - lattice_projectionTrans( sum(lattice_NtransSystem(1:f-1,phase))+j, & - sum(lattice_NslipSystem(1:o-1,phase))+k, phase) - enddo; enddo - end select - - enddo transSystemsLoop - enddo transFamiliesLoop - - startIndex=1_pInt - endIndex=ns - state(instance)%rhoEdge=>plasticState(phase)%state(startIndex:endIndex,:) - state0(instance)%rhoEdge=>plasticState(phase)%state0(startIndex:endIndex,:) - dotState(instance)%rhoEdge=>plasticState(phase)%dotState(startIndex:endIndex,:) - - startIndex=endIndex+1 - endIndex=endIndex+ns - state(instance)%rhoEdgeDip=>plasticState(phase)%state(startIndex:endIndex,:) - state0(instance)%rhoEdgeDip=>plasticState(phase)%state0(startIndex:endIndex,:) - dotState(instance)%rhoEdgeDip=>plasticState(phase)%dotState(startIndex:endIndex,:) - - startIndex=endIndex+1 - endIndex=endIndex+ns - state(instance)%accshear_slip=>plasticState(phase)%state(startIndex:endIndex,:) - state0(instance)%accshear_slip=>plasticState(phase)%state0(startIndex:endIndex,:) - dotState(instance)%accshear_slip=>plasticState(phase)%dotState(startIndex:endIndex,:) - - startIndex=endIndex+1 - endIndex=endIndex+nt - state(instance)%twinFraction=>plasticState(phase)%state(startIndex:endIndex,:) - state0(instance)%twinFraction=>plasticState(phase)%state0(startIndex:endIndex,:) - dotState(instance)%twinFraction=>plasticState(phase)%dotState(startIndex:endIndex,:) - - startIndex=endIndex+1 - endIndex=endIndex+nt - state(instance)%accshear_twin=>plasticState(phase)%state(startIndex:endIndex,:) - state0(instance)%accshear_twin=>plasticState(phase)%state0(startIndex:endIndex,:) - dotState(instance)%accshear_twin=>plasticState(phase)%dotState(startIndex:endIndex,:) - - startIndex=endIndex+1 - endIndex=endIndex+nr - state(instance)%stressTransFraction=>plasticState(phase)%state(startIndex:endIndex,:) - state0(instance)%stressTransFraction=>plasticState(phase)%state0(startIndex:endIndex,:) - dotState(instance)%stressTransFraction=>plasticState(phase)%dotState(startIndex:endIndex,:) - - startIndex=endIndex+1 - endIndex=endIndex+nr - state(instance)%strainTransFraction=>plasticState(phase)%state(startIndex:endIndex,:) - state0(instance)%strainTransFraction=>plasticState(phase)%state0(startIndex:endIndex,:) - dotState(instance)%strainTransFraction=>plasticState(phase)%dotState(startIndex:endIndex,:) - - startIndex=endIndex+1 - endIndex=endIndex+ns - state(instance)%invLambdaSlip=>plasticState(phase)%state(startIndex:endIndex,:) - state0(instance)%invLambdaSlip=>plasticState(phase)%state0(startIndex:endIndex,:) - - startIndex=endIndex+1 - endIndex=endIndex+ns - state(instance)%invLambdaSlipTwin=>plasticState(phase)%state(startIndex:endIndex,:) - state0(instance)%invLambdaSlipTwin=>plasticState(phase)%state0(startIndex:endIndex,:) - - startIndex=endIndex+1 - endIndex=endIndex+nt - state(instance)%invLambdaTwin=>plasticState(phase)%state(startIndex:endIndex,:) - state0(instance)%invLambdaTwin=>plasticState(phase)%state0(startIndex:endIndex,:) - - startIndex=endIndex+1 - endIndex=endIndex+ns - state(instance)%invLambdaSlipTrans=>plasticState(phase)%state(startIndex:endIndex,:) - state0(instance)%invLambdaSlipTrans=>plasticState(phase)%state0(startIndex:endIndex,:) - - startIndex=endIndex+1 - endIndex=endIndex+nr - state(instance)%invLambdaTrans=>plasticState(phase)%state(startIndex:endIndex,:) - state0(instance)%invLambdaTrans=>plasticState(phase)%state0(startIndex:endIndex,:) - - startIndex=endIndex+1 - endIndex=endIndex+ns - state(instance)%mfp_slip=>plasticState(phase)%state(startIndex:endIndex,:) - state0(instance)%mfp_slip=>plasticState(phase)%state0(startIndex:endIndex,:) - - startIndex=endIndex+1 - endIndex=endIndex+nt - state(instance)%mfp_twin=>plasticState(phase)%state(startIndex:endIndex,:) - state0(instance)%mfp_twin=>plasticState(phase)%state0(startIndex:endIndex,:) - - startIndex=endIndex+1 - endIndex=endIndex+nr - state(instance)%mfp_trans=>plasticState(phase)%state(startIndex:endIndex,:) - state0(instance)%mfp_trans=>plasticState(phase)%state0(startIndex:endIndex,:) - - startIndex=endIndex+1 - endIndex=endIndex+ns - state(instance)%threshold_stress_slip=>plasticState(phase)%state(startIndex:endIndex,:) - state0(instance)%threshold_stress_slip=>plasticState(phase)%state0(startIndex:endIndex,:) - - startIndex=endIndex+1 - endIndex=endIndex+nt - state(instance)%threshold_stress_twin=>plasticState(phase)%state(startIndex:endIndex,:) - state0(instance)%threshold_stress_twin=>plasticState(phase)%state0(startIndex:endIndex,:) - - startIndex=endIndex+1 - endIndex=endIndex+nr - state(instance)%threshold_stress_trans=>plasticState(phase)%state(startIndex:endIndex,:) - state0(instance)%threshold_stress_trans=>plasticState(phase)%state0(startIndex:endIndex,:) - - startIndex=endIndex+1 - endIndex=endIndex+nt - state(instance)%twinVolume=>plasticState(phase)%state(startIndex:endIndex,:) - state0(instance)%twinVolume=>plasticState(phase)%state0(startIndex:endIndex,:) - - startIndex=endIndex+1 - endIndex=endIndex+nr - state(instance)%martensiteVolume=>plasticState(phase)%state(startIndex:endIndex,:) - state0(instance)%martensiteVolume=>plasticState(phase)%state0(startIndex:endIndex,:) - - call plastic_dislotwin_stateInit(phase,instance) - call plastic_dislotwin_aTolState(phase,instance) - endif myPhase2 - - enddo initializeInstances -end subroutine plastic_dislotwin_init - -!-------------------------------------------------------------------------------------------------- -!> @brief sets the relevant state values for a given instance of this plasticity -!-------------------------------------------------------------------------------------------------- -subroutine plastic_dislotwin_stateInit(ph,instance) - use math, only: & - pi - use lattice, only: & - lattice_maxNslipFamily, & - lattice_mu - use material, only: & - plasticState - - implicit none - integer(pInt), intent(in) :: & - instance, & !< number specifying the instance of the plasticity - ph - - real(pReal), dimension(plasticState(ph)%sizeState) :: tempState - - integer(pInt) :: i,j,f,ns,nt,nr, index_myFamily - real(pReal), dimension(plastic_dislotwin_totalNslip(instance)) :: & - rhoEdge0, & - rhoEdgeDip0, & - invLambdaSlip0, & - MeanFreePathSlip0, & - tauSlipThreshold0 - real(pReal), dimension(plastic_dislotwin_totalNtwin(instance)) :: & - MeanFreePathTwin0,TwinVolume0 - real(pReal), dimension(plastic_dislotwin_totalNtrans(instance)) :: & - MeanFreePathTrans0,MartensiteVolume0 - tempState = 0.0_pReal - ns = plastic_dislotwin_totalNslip(instance) - nt = plastic_dislotwin_totalNtwin(instance) - nr = plastic_dislotwin_totalNtrans(instance) - -!-------------------------------------------------------------------------------------------------- -! initialize basic slip state variables - do f = 1_pInt,lattice_maxNslipFamily - index_myFamily = sum(plastic_dislotwin_Nslip(1:f-1_pInt,instance)) ! index in truncated slip system list - rhoEdge0(index_myFamily+1_pInt: & - index_myFamily+plastic_dislotwin_Nslip(f,instance)) = & - plastic_dislotwin_rhoEdge0(f,instance) - rhoEdgeDip0(index_myFamily+1_pInt: & - index_myFamily+plastic_dislotwin_Nslip(f,instance)) = & - plastic_dislotwin_rhoEdgeDip0(f,instance) + end associate enddo - tempState(1_pInt:ns) = rhoEdge0 - tempState(ns+1_pInt:2_pInt*ns) = rhoEdgeDip0 - -!-------------------------------------------------------------------------------------------------- -! initialize dependent slip microstructural variables - forall (i = 1_pInt:ns) & - invLambdaSlip0(i) = sqrt(dot_product((rhoEdge0+rhoEdgeDip0),plastic_dislotwin_forestProjectionEdge(1:ns,i,instance)))/ & - plastic_dislotwin_CLambdaSlipPerSlipSystem(i,instance) - tempState(3_pInt*ns+2_pInt*nt+2_pInt*nr+1:4_pInt*ns+2_pInt*nt+2_pInt*nr) = invLambdaSlip0 - - forall (i = 1_pInt:ns) & - MeanFreePathSlip0(i) = & - plastic_dislotwin_GrainSize(instance)/(1.0_pReal+invLambdaSlip0(i)*plastic_dislotwin_GrainSize(instance)) - tempState(6_pInt*ns+3_pInt*nt+3_pInt*nr+1:7_pInt*ns+3_pInt*nt+3_pInt*nr) = MeanFreePathSlip0 - - forall (i = 1_pInt:ns) & - tauSlipThreshold0(i) = & - lattice_mu(ph)*plastic_dislotwin_burgersPerSlipSystem(i,instance) * & - sqrt(dot_product((rhoEdge0+rhoEdgeDip0),plastic_dislotwin_interactionMatrix_SlipSlip(i,1:ns,instance))) - - tempState(7_pInt*ns+4_pInt*nt+4_pInt*nr+1:8_pInt*ns+4_pInt*nt+4_pInt*nr) = tauSlipThreshold0 - -!-------------------------------------------------------------------------------------------------- -! initialize dependent twin microstructural variables - forall (j = 1_pInt:nt) & - MeanFreePathTwin0(j) = plastic_dislotwin_GrainSize(instance) - tempState(7_pInt*ns+3_pInt*nt+3_pInt*nr+1_pInt:7_pInt*ns+4_pInt*nt+3_pInt*nr) = MeanFreePathTwin0 - - forall (j = 1_pInt:nt) & - TwinVolume0(j) = & - (pi/4.0_pReal)*plastic_dislotwin_twinsizePerTwinSystem(j,instance)*MeanFreePathTwin0(j)**(2.0_pReal) - tempState(8_pInt*ns+5_pInt*nt+5_pInt*nr+1_pInt:8_pInt*ns+6_pInt*nt+5_pInt*nr) = TwinVolume0 - -!-------------------------------------------------------------------------------------------------- -! initialize dependent trans microstructural variables - forall (j = 1_pInt:nr) & - MeanFreePathTrans0(j) = plastic_dislotwin_GrainSize(instance) - tempState(7_pInt*ns+4_pInt*nt+3_pInt*nr+1_pInt:7_pInt*ns+4_pInt*nt+4_pInt*nr) = MeanFreePathTrans0 - - forall (j = 1_pInt:nr) & - MartensiteVolume0(j) = & - (pi/4.0_pReal)*plastic_dislotwin_lamellarsizePerTransSystem(j,instance)*MeanFreePathTrans0(j)**(2.0_pReal) - tempState(8_pInt*ns+6_pInt*nt+5_pInt*nr+1_pInt:8_pInt*ns+6_pInt*nt+6_pInt*nr) = MartensiteVolume0 - -plasticState(ph)%state0 = spread(tempState,2,size(plasticState(ph)%state(1,:))) - -end subroutine plastic_dislotwin_stateInit - -!-------------------------------------------------------------------------------------------------- -!> @brief sets the relevant state values for a given instance of this plasticity -!-------------------------------------------------------------------------------------------------- -subroutine plastic_dislotwin_aTolState(ph,instance) - use material, only: & - plasticState - - implicit none - integer(pInt), intent(in) :: & - ph, & - instance ! number specifying the current instance of the plasticity - - integer(pInt) :: ns, nt, nr - - ns = plastic_dislotwin_totalNslip(instance) - nt = plastic_dislotwin_totalNtwin(instance) - nr = plastic_dislotwin_totalNtrans(instance) - - ! Tolerance state for dislocation densities - plasticState(ph)%aTolState(1_pInt: & - 2_pInt*ns) = plastic_dislotwin_aTolRho(instance) - - ! Tolerance state for accumulated shear due to slip - plasticState(ph)%aTolState(2_pInt*ns+1_pInt: & - 3_pInt*ns)=1.0e6_pReal - - ! Tolerance state for twin volume fraction - plasticState(ph)%aTolState(3_pInt*ns+1_pInt: & - 3_pInt*ns+nt) = plastic_dislotwin_aTolTwinFrac(instance) - - ! Tolerance state for accumulated shear due to twin - plasticState(ph)%aTolState(3_pInt*ns+nt+1_pInt: & - 3_pInt*ns+2_pInt*nt) = 1.0e6_pReal - -! Tolerance state for stress-assisted martensite volume fraction - plasticState(ph)%aTolState(3_pInt*ns+2_pInt*nt+1_pInt: & - 3_pInt*ns+2_pInt*nt+nr) = plastic_dislotwin_aTolTransFrac(instance) - -! Tolerance state for strain-induced martensite volume fraction - plasticState(ph)%aTolState(3_pInt*ns+2_pInt*nt+nr+1_pInt: & - 3_pInt*ns+2_pInt*nt+2_pInt*nr) = plastic_dislotwin_aTolTransFrac(instance) - -end subroutine plastic_dislotwin_aTolState - +end subroutine plastic_dislotwin_init !-------------------------------------------------------------------------------------------------- !> @brief returns the homogenized elasticity matrix !-------------------------------------------------------------------------------------------------- function plastic_dislotwin_homogenizedC(ipc,ip,el) use material, only: & + material_phase, & phase_plasticityInstance, & - phaseAt, phasememberAt - use lattice, only: & - lattice_C66 + phasememberAt - implicit none - real(pReal), dimension(6,6) :: & - plastic_dislotwin_homogenizedC - integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element + implicit none + real(pReal), dimension(6,6) :: & + plastic_dislotwin_homogenizedC + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + type(tParameters) :: prm + type(tDislotwinState) :: stt - integer(pInt) :: instance,ns,nt,nr,i, & - ph, & + integer(pInt) :: i, & of - real(pReal) :: sumf, sumftr + real(pReal) :: f_unrotated - !* Shortened notation of = phasememberAt(ipc,ip,el) - ph = phaseAt(ipc,ip,el) - instance = phase_plasticityInstance(ph) - ns = plastic_dislotwin_totalNslip(instance) - nt = plastic_dislotwin_totalNtwin(instance) - nr = plastic_dislotwin_totalNtrans(instance) + associate(prm => param(phase_plasticityInstance(material_phase(ipc,ip,el))),& + stt => state(phase_plasticityInstance(material_phase(ipc,ip,el)))) - !* Total twin volume fraction - sumf = sum(state(instance)%twinFraction(1_pInt:nt,of)) ! safe for nt == 0 + f_unrotated = 1.0_pReal & + - sum(stt%twinFraction(1_pInt:prm%totalNtwin,of)) & + - sum(stt%stressTransFraction(1_pInt:prm%totalNtrans,of)) & + - sum(stt%strainTransFraction(1_pInt:prm%totalNtrans,of)) - !* Total transformed volume fraction - sumftr = sum(state(instance)%stressTransFraction(1_pInt:nr,of)) + & - sum(state(instance)%strainTransFraction(1_pInt:nr,of)) - - !* Homogenized elasticity matrix - plastic_dislotwin_homogenizedC = (1.0_pReal-sumf-sumftr)*lattice_C66(1:6,1:6,ph) - do i=1_pInt,nt + plastic_dislotwin_homogenizedC = f_unrotated * prm%C66 + do i=1_pInt,prm%totalNtwin plastic_dislotwin_homogenizedC = plastic_dislotwin_homogenizedC & - + state(instance)%twinFraction(i,of)*plastic_dislotwin_Ctwin66(1:6,1:6,i,instance) + + stt%twinFraction(i,of)*prm%C66_twin(1:6,1:6,i) enddo - do i=1_pInt,nr + do i=1_pInt,prm%totalNtrans plastic_dislotwin_homogenizedC = plastic_dislotwin_homogenizedC & - + (state(instance)%stressTransFraction(i,of) + state(instance)%strainTransFraction(i,of))*& - plastic_dislotwin_Ctrans66(1:6,1:6,i,instance) + +(stt%stressTransFraction(i,of)+stt%strainTransFraction(i,of))*& + prm%C66_trans(1:6,1:6,i) enddo - + end associate end function plastic_dislotwin_homogenizedC - + + !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- subroutine plastic_dislotwin_microstructure(temperature,ipc,ip,el) use math, only: & - pi + PI use material, only: & material_phase, & phase_plasticityInstance, & - !plasticState, & !!!!delete - phaseAt, phasememberAt - use lattice, only: & - lattice_mu, & - lattice_nu + phasememberAt implicit none integer(pInt), intent(in) :: & @@ -1470,222 +829,154 @@ subroutine plastic_dislotwin_microstructure(temperature,ipc,ip,el) temperature !< temperature at IP integer(pInt) :: & - instance, & - ns,nt,nr,s,t,r, & - ph, & + i, & of real(pReal) :: & - sumf,sfe,x0,sumftr - real(pReal), dimension(plastic_dislotwin_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: fOverStacksize - real(pReal), dimension(plastic_dislotwin_totalNtrans(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + sumf_twin,SFE,sumf_trans + real(pReal), dimension(:), allocatable :: & + x0, & + fOverStacksize, & ftransOverLamellarSize + + type(tParameters) :: prm !< parameters of present instance + type(tDislotwinState) :: stt !< state of present instance + type(tDislotwinMicrostructure) :: mse - !* Shortened notation of = phasememberAt(ipc,ip,el) - ph = phaseAt(ipc,ip,el) - instance = phase_plasticityInstance(ph) - ns = plastic_dislotwin_totalNslip(instance) - nt = plastic_dislotwin_totalNtwin(instance) - nr = plastic_dislotwin_totalNtrans(instance) - !* Total twin volume fraction - sumf = sum(state(instance)%twinFraction(1_pInt:nt,of)) ! safe for nt == 0 - - !* Total transformed volume fraction - sumftr = sum(state(instance)%stressTransFraction(1_pInt:nr,of)) + & - sum(state(instance)%strainTransFraction(1_pInt:nr,of)) + associate(prm => param(phase_plasticityInstance(material_phase(ipc,ip,el))),& + stt => state(phase_plasticityInstance(material_phase(ipc,ip,el))),& + mse => microstructure(phase_plasticityInstance(material_phase(ipc,ip,el)))) - !* Stacking fault energy - sfe = plastic_dislotwin_SFE_0K(instance) + & - plastic_dislotwin_dSFE_dT(instance) * Temperature - - !* rescaled twin volume fraction for topology - forall (t = 1_pInt:nt) & - fOverStacksize(t) = & - state(instance)%twinFraction(t,of)/plastic_dislotwin_twinsizePerTwinSystem(t,instance) + sumf_twin = sum(stt%twinFraction(1:prm%totalNtwin,of)) + sumf_trans = sum(stt%stressTransFraction(1:prm%totalNtrans,of)) & + + sum(stt%strainTransFraction(1:prm%totalNtrans,of)) - !* rescaled trans volume fraction for topology - forall (r = 1_pInt:nr) & - ftransOverLamellarSize(r) = & - (state(instance)%stressTransFraction(r,of)+state(instance)%strainTransFraction(r,of))/& - plastic_dislotwin_lamellarsizePerTransSystem(r,instance) + sfe = prm%SFE_0K + prm%dSFE_dT * Temperature + !* rescaled volume fraction for topology + fOverStacksize = stt%twinFraction(1_pInt:prm%totalNtwin,of)/prm%twinsize !ToDo: this is per system + ftransOverLamellarSize = sumf_trans/prm%lamellarsizePerTransSystem !ToDo: But this not ... + !Todo: Physically ok, but naming could be adjusted + + !* 1/mean free distance between 2 forest dislocations seen by a moving dislocation - forall (s = 1_pInt:ns) & - state(instance)%invLambdaSlip(s,of) = & - sqrt(dot_product((state(instance)%rhoEdge(1_pInt:ns,of)+state(instance)%rhoEdgeDip(1_pInt:ns,of)),& - plastic_dislotwin_forestProjectionEdge(1:ns,s,instance)))/ & - plastic_dislotwin_CLambdaSlipPerSlipSystem(s,instance) + forall (i = 1_pInt:prm%totalNslip) & + mse%invLambdaSlip(i,of) = & + sqrt(dot_product((stt%rhoEdge(1_pInt:prm%totalNslip,of)+stt%rhoEdgeDip(1_pInt:prm%totalNslip,of)),& + prm%forestProjectionEdge(1:prm%totalNslip,i)))/prm%CLambdaSlip(i) !* 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation !$OMP CRITICAL (evilmatmul) - state(instance)%invLambdaSlipTwin(1_pInt:ns,of) = 0.0_pReal - if (nt > 0_pInt .and. ns > 0_pInt) & - state(instance)%invLambdaSlipTwin(1_pInt:ns,of) = & - matmul(plastic_dislotwin_interactionMatrix_SlipTwin(1:ns,1:nt,instance),fOverStacksize(1:nt))/(1.0_pReal-sumf) - !$OMP END CRITICAL (evilmatmul) - + if (prm%totalNtwin > 0_pInt .and. prm%totalNslip > 0_pInt) & + mse%invLambdaSlipTwin(1_pInt:prm%totalNslip,of) = & + matmul(prm%interaction_SlipTwin,fOverStacksize)/(1.0_pReal-sumf_twin) + !* 1/mean free distance between 2 twin stacks from different systems seen by a growing twin - !$OMP CRITICAL (evilmatmul) - if (nt > 0_pInt) & - state(instance)%invLambdaTwin(1_pInt:nt,of) = & - matmul(plastic_dislotwin_interactionMatrix_TwinTwin(1:nt,1:nt,instance),fOverStacksize(1:nt))/(1.0_pReal-sumf) - !$OMP END CRITICAL (evilmatmul) + + !ToDo: needed? if (prm%totalNtwin > 0_pInt) & + mse%invLambdaTwin(1_pInt:prm%totalNtwin,of) = & + matmul(prm%interaction_TwinTwin,fOverStacksize)/(1.0_pReal-sumf_twin) + !* 1/mean free distance between 2 martensite lamellar from different systems seen by a moving dislocation - state(instance)%invLambdaSlipTrans(1_pInt:ns,of) = 0.0_pReal - if (nr > 0_pInt .and. ns > 0_pInt) & - state(instance)%invLambdaSlipTrans(1_pInt:ns,of) = & - matmul(plastic_dislotwin_interactionMatrix_SlipTrans(1:ns,1:nr,instance),ftransOverLamellarSize(1:nr))/(1.0_pReal-sumftr) + if (prm%totalNtrans > 0_pInt .and. prm%totalNslip > 0_pInt) & + mse%invLambdaSlipTrans(1_pInt:prm%totalNslip,of) = & + matmul(prm%interaction_SlipTrans,ftransOverLamellarSize)/(1.0_pReal-sumf_trans) !* 1/mean free distance between 2 martensite stacks from different systems seen by a growing martensite (1/lambda_trans) - if (nr > 0_pInt) & - state(instance)%invLambdaTrans(1_pInt:nr,of) = & - matmul(plastic_dislotwin_interactionMatrix_TransTrans(1:nr,1:nr,instance),ftransOverLamellarSize(1:nr))/(1.0_pReal-sumftr) + !ToDo: needed? if (prm%totalNtrans > 0_pInt) & + + mse%invLambdaTrans(1_pInt:prm%totalNtrans,of) = & + matmul(prm%interaction_TransTrans,ftransOverLamellarSize)/(1.0_pReal-sumf_trans) + !$OMP END CRITICAL (evilmatmul) !* mean free path between 2 obstacles seen by a moving dislocation - do s = 1_pInt,ns - if ((nt > 0_pInt) .or. (nr > 0_pInt)) then - state(instance)%mfp_slip(s,of) = & - plastic_dislotwin_GrainSize(instance)/(1.0_pReal+plastic_dislotwin_GrainSize(instance)*& - (state(instance)%invLambdaSlip(s,of) + & - state(instance)%invLambdaSlipTwin(s,of) + & - state(instance)%invLambdaSlipTrans(s,of))) + do i = 1_pInt,prm%totalNslip + if ((prm%totalNtwin > 0_pInt) .or. (prm%totalNtrans > 0_pInt)) then ! ToDo: This is too simplified + mse%mfp_slip(i,of) = & + prm%GrainSize/(1.0_pReal+prm%GrainSize*& + (mse%invLambdaSlip(i,of) + mse%invLambdaSlipTwin(i,of) + mse%invLambdaSlipTrans(i,of))) else - state(instance)%mfp_slip(s,of) = & - plastic_dislotwin_GrainSize(instance)/& - (1.0_pReal+plastic_dislotwin_GrainSize(instance)*(state(instance)%invLambdaSlip(s,of))) !!!!!! correct? + mse%mfp_slip(i,of) = & + prm%GrainSize/& + (1.0_pReal+prm%GrainSize*(mse%invLambdaSlip(i,of))) !!!!!! correct? endif enddo - !* mean free path between 2 obstacles seen by a growing twin - forall (t = 1_pInt:nt) & - state(instance)%mfp_twin(t,of) = & - plastic_dislotwin_Cmfptwin(instance)*plastic_dislotwin_GrainSize(instance)/& - (1.0_pReal+plastic_dislotwin_GrainSize(instance)*state(ph)%invLambdaTwin(t,of)) - - !* mean free path between 2 obstacles seen by a growing martensite - forall (r = 1_pInt:nr) & - state(instance)%mfp_trans(r,of) = & - plastic_dislotwin_Cmfptrans(instance)*plastic_dislotwin_GrainSize(instance)/& - (1.0_pReal+plastic_dislotwin_GrainSize(instance)*state(instance)%invLambdaTrans(r,of)) + !* mean free path between 2 obstacles seen by a growing twin/martensite + mse%mfp_twin(:,of) = prm%Cmfptwin*prm%GrainSize/ (1.0_pReal+prm%GrainSize*mse%invLambdaTwin(:,of)) + mse%mfp_trans(:,of) = prm%Cmfptrans*prm%GrainSize/(1.0_pReal+prm%GrainSize*mse%invLambdaTrans(:,of)) !* threshold stress for dislocation motion - forall (s = 1_pInt:ns) & - state(instance)%threshold_stress_slip(s,of) = & - lattice_mu(ph)*plastic_dislotwin_burgersPerSlipSystem(s,instance)*& - sqrt(dot_product((state(instance)%rhoEdge(1_pInt:ns,of)+state(instance)%rhoEdgeDip(1_pInt:ns,of)),& - plastic_dislotwin_interactionMatrix_SlipSlip(s,1:ns,instance))) + forall (i = 1_pInt:prm%totalNslip) mse%threshold_stress_slip(i,of) = & + prm%mu*prm%burgers_slip(i)*& + sqrt(dot_product(stt%rhoEdge(1_pInt:prm%totalNslip,of)+stt%rhoEdgeDip(1_pInt:prm%totalNslip,of),& + prm%interaction_SlipSlip(i,1:prm%totalNslip))) - !* threshold stress for growing twin - forall (t = 1_pInt:nt) & - state(instance)%threshold_stress_twin(t,of) = & - plastic_dislotwin_Cthresholdtwin(instance)* & - (sfe/(3.0_pReal*plastic_dislotwin_burgersPerTwinSystem(t,instance)) & - + 3.0_pReal*plastic_dislotwin_burgersPerTwinSystem(t,instance)*lattice_mu(ph)/& - (plastic_dislotwin_L0_twin(instance)*plastic_dislotwin_burgersPerSlipSystem(t,instance)) & - ) - - !* threshold stress for growing martensite - forall (r = 1_pInt:nr) & - state(instance)%threshold_stress_trans(r,of) = & - plastic_dislotwin_Cthresholdtrans(instance)* & - (sfe/(3.0_pReal*plastic_dislotwin_burgersPerTransSystem(r,instance)) & - + 3.0_pReal*plastic_dislotwin_burgersPerTransSystem(r,instance)*lattice_mu(ph)/& - (plastic_dislotwin_L0_trans(instance)*plastic_dislotwin_burgersPerSlipSystem(r,instance))& - + plastic_dislotwin_transStackHeight(instance)*plastic_dislotwin_deltaG(instance)/ & - (3.0_pReal*plastic_dislotwin_burgersPerTransSystem(r,instance)) & - ) - - !* final twin volume after growth - forall (t = 1_pInt:nt) & - state(instance)%twinVolume(t,of) = & - (pi/4.0_pReal)*plastic_dislotwin_twinsizePerTwinSystem(t,instance)*& - state(instance)%mfp_twin(t,of)**(2.0_pReal) - - !* final martensite volume after growth - forall (r = 1_pInt:nr) & - state(instance)%martensiteVolume(r,of) = & - (pi/4.0_pReal)*plastic_dislotwin_lamellarsizePerTransSystem(r,instance)*& - state(instance)%mfp_trans(r,of)**(2.0_pReal) + !* threshold stress for growing twin/martensite + if(prm%totalNtwin == prm%totalNslip) & + mse%threshold_stress_twin(:,of) = prm%Cthresholdtwin* & + (sfe/(3.0_pReal*prm%burgers_twin)+ 3.0_pReal*prm%burgers_twin*prm%mu/ & + (prm%L0_twin*prm%burgers_slip)) ! slip burgers here correct? + if(prm%totalNtrans == prm%totalNslip) & + mse%threshold_stress_trans(:,of) = prm%Cthresholdtrans* & + (sfe/(3.0_pReal*prm%burgers_trans) + 3.0_pReal*prm%burgers_trans*prm%mu/& + (prm%L0_trans*prm%burgers_slip) + prm%transStackHeight*prm%deltaG/ (3.0_pReal*prm%burgers_trans) ) + + ! final volume after growth + mse%twinVolume(:,of) = (PI/4.0_pReal)*prm%twinsize*mse%mfp_twin(:,of)**2.0_pReal + mse%martensiteVolume(:,of) = (PI/4.0_pReal)*prm%lamellarsizePerTransSystem*mse%mfp_trans(:,of)**2.0_pReal !* equilibrium separation of partial dislocations (twin) - do t = 1_pInt,nt - x0 = lattice_mu(ph)*plastic_dislotwin_burgersPerTwinSystem(t,instance)**(2.0_pReal)/& - (sfe*8.0_pReal*pi)*(2.0_pReal+lattice_nu(ph))/(1.0_pReal-lattice_nu(ph)) - plastic_dislotwin_tau_r_twin(t,instance)= & - lattice_mu(ph)*plastic_dislotwin_burgersPerTwinSystem(t,instance)/(2.0_pReal*pi)*& - (1/(x0+plastic_dislotwin_xc_twin(instance))+cos(pi/3.0_pReal)/x0) - enddo + x0 = prm%mu*prm%burgers_twin**2.0_pReal/(sfe*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) + mse%tau_r_twin(:,of) = prm%mu*prm%burgers_twin/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%xc_twin)+cos(pi/3.0_pReal)/x0) !* equilibrium separation of partial dislocations (trans) - do r = 1_pInt,nr - x0 = lattice_mu(ph)*plastic_dislotwin_burgersPerTransSystem(r,instance)**(2.0_pReal)/& - (sfe*8.0_pReal*pi)*(2.0_pReal+lattice_nu(ph))/(1.0_pReal-lattice_nu(ph)) - plastic_dislotwin_tau_r_trans(r,instance)= & - lattice_mu(ph)*plastic_dislotwin_burgersPerTransSystem(r,instance)/(2.0_pReal*pi)*& - (1/(x0+plastic_dislotwin_xc_trans(instance))+cos(pi/3.0_pReal)/x0) - enddo + x0 = prm%mu*prm%burgers_trans**2.0_pReal/(sfe*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) + mse%tau_r_trans(:,of) = prm%mu*prm%burgers_trans/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%xc_trans)+cos(pi/3.0_pReal)/x0) +end associate end subroutine plastic_dislotwin_microstructure !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,ipc,ip,el) +subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,Temperature,instance,of) use prec, only: & tol_math_check, & dNeq0 use math, only: & - math_Plain3333to99, & - math_Mandel6to33, & - math_Mandel33to6, & math_eigenValuesVectorsSym, & math_tensorproduct33, & math_symmetric33, & + math_mul33xx33, & math_mul33x3 use material, only: & material_phase, & phase_plasticityInstance, & - phaseAt, phasememberAt - use lattice, only: & - lattice_Sslip, & - lattice_Sslip_v, & - lattice_Stwin, & - lattice_Stwin_v, & - lattice_Strans, & - lattice_Strans_v, & - lattice_maxNslipFamily,& - lattice_maxNtwinFamily, & - lattice_maxNtransFamily, & - lattice_NslipSystem, & - lattice_NtwinSystem, & - lattice_NtransSystem, & - lattice_shearTwin, & - lattice_structure, & - lattice_fcc_twinNucleationSlipPair, & - LATTICE_fcc_ID + phasememberAt implicit none - integer(pInt), intent(in) :: ipc,ip,el - real(pReal), intent(in) :: Temperature - real(pReal), dimension(6), intent(in) :: Tstar_v - real(pReal), dimension(3,3), intent(out) :: Lp - real(pReal), dimension(9,9), intent(out) :: dLp_dTstar99 + real(pReal), dimension(3,3), intent(out) :: Lp + real(pReal), dimension(3,3,3,3), intent(out) :: dLp_dMp + real(pReal), dimension(3,3), intent(in) :: Mp + integer(pInt), intent(in) :: instance,of + real(pReal), intent(in) :: Temperature - integer(pInt) :: instance,ph,of,ns,nt,nr,f,i,j,k,l,m,n,index_myFamily,s1,s2 - real(pReal) :: sumf,sumftr,StressRatio_p,StressRatio_pminus1,StressRatio_r,BoltzmannRatio,DotGamma0,Ndot0_twin,stressRatio, & - Ndot0_trans,StressRatio_s - real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 - real(pReal), dimension(plastic_dislotwin_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & - gdot_slip,dgdot_dtauslip,tau_slip - real(pReal), dimension(plastic_dislotwin_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & - gdot_twin,dgdot_dtautwin,tau_twin - real(pReal), dimension(plastic_dislotwin_totalNtrans(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & - gdot_trans,dgdot_dtautrans,tau_trans - real(pReal), dimension(6) :: gdot_sb,dgdot_dtausb,tau_sb - real(pReal), dimension(3,3) :: eigVectors, sb_Smatrix + integer(pInt) :: i,k,l,m,n,s1,s2 + real(pReal) :: f_unrotated,StressRatio_p,& + StressRatio_r,BoltzmannRatio,Ndot0_twin,stressRatio, & + Ndot0_trans,StressRatio_s, & + dgdot_dtau, & + tau + real(pReal), dimension(param(instance)%totalNslip) :: & + gdot_slip,dgdot_dtau_slip + real(pReal), dimension(param(instance)%totalNtwin) :: & + gdot_twin,dgdot_dtau_twin + real(pReal):: gdot_sb,gdot_trans + real(pReal), dimension(3,3) :: eigVectors, Schmid_shearBand real(pReal), dimension(3) :: eigValues, sb_s, sb_m logical :: error real(pReal), dimension(3,6), parameter :: & @@ -1707,239 +998,104 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature 0, 1,-1, & 0, 1, 1 & ],pReal),[ 3,6]) - !* Shortened notation - of = phasememberAt(ipc,ip,el) - ph = phaseAt(ipc,ip,el) - instance = phase_plasticityInstance(ph) - ns = plastic_dislotwin_totalNslip(instance) - nt = plastic_dislotwin_totalNtwin(instance) - nr = plastic_dislotwin_totalNtrans(instance) + + type(tParameters) :: prm !< parameters of present instance + type(tDislotwinState) :: ste !< state of present instance + + associate(prm => param(instance), stt => state(instance), mse => microstructure(instance)) + + f_unrotated = 1.0_pReal & + - sum(stt%twinFraction(1_pInt:prm%totalNtwin,of)) & + - sum(stt%stressTransFraction(1_pInt:prm%totalNtrans,of)) & + - sum(stt%strainTransFraction(1_pInt:prm%totalNtrans,of)) Lp = 0.0_pReal - dLp_dTstar3333 = 0.0_pReal + dLp_dMp = 0.0_pReal -!-------------------------------------------------------------------------------------------------- -! Dislocation glide part - gdot_slip = 0.0_pReal - dgdot_dtauslip = 0.0_pReal - j = 0_pInt - slipFamiliesLoop: do f = 1_pInt,lattice_maxNslipFamily - index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family - slipSystemsLoop: do i = 1_pInt,plastic_dislotwin_Nslip(f,instance) - j = j+1_pInt + call kinetics_slip(prm,stt,mse,of,Mp,temperature,gdot_slip,dgdot_dtau_slip) + slipContribution: do i = 1_pInt, prm%totalNslip + Lp = Lp + gdot_slip(i)*prm%Schmid_slip(1:3,1:3,i) + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + dgdot_dtau_slip(i) * prm%Schmid_slip(k,l,i) * prm%Schmid_slip(m,n,i) + enddo slipContribution - !* Calculation of Lp - !* Resolved shear stress on slip system - tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph)) + !ToDo: Why do this before shear banding? + Lp = Lp * f_unrotated + dLp_dMp = dLp_dMp * f_unrotated + + shearBandingContribution: if(dNeq0(prm%sbVelocity)) then - if((abs(tau_slip(j))-state(instance)%threshold_stress_slip(j,of)) > tol_math_check) then - !* Stress ratios - stressRatio =((abs(tau_slip(j))- state(instance)%threshold_stress_slip(j,of))/& - (plastic_dislotwin_SolidSolutionStrength(instance)+plastic_dislotwin_tau_peierlsPerSlipFamily(f,instance))) - StressRatio_p = stressRatio** plastic_dislotwin_pPerSlipFamily(f,instance) - StressRatio_pminus1 = stressRatio**(plastic_dislotwin_pPerSlipFamily(f,instance)-1.0_pReal) - !* Boltzmann ratio - BoltzmannRatio = plastic_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature) - !* Initial shear rates - DotGamma0 = & - state(instance)%rhoEdge(j,of)*plastic_dislotwin_burgersPerSlipSystem(j,instance)*& - plastic_dislotwin_v0PerSlipSystem(j,instance) - - !* Shear rates due to slip - gdot_slip(j) = DotGamma0 & - * exp(-BoltzmannRatio*(1-StressRatio_p) ** plastic_dislotwin_qPerSlipFamily(f,instance)) & - * sign(1.0_pReal,tau_slip(j)) - - !* Derivatives of shear rates - dgdot_dtauslip(j) = & - abs(gdot_slip(j))*BoltzmannRatio*plastic_dislotwin_pPerSlipFamily(f,instance)& - *plastic_dislotwin_qPerSlipFamily(f,instance)/& - (plastic_dislotwin_SolidSolutionStrength(instance)+plastic_dislotwin_tau_peierlsPerSlipFamily(f,instance))*& - StressRatio_pminus1*(1-StressRatio_p)**(plastic_dislotwin_qPerSlipFamily(f,instance)-1.0_pReal) - endif - - !* Plastic velocity gradient for dislocation glide - Lp = Lp + gdot_slip(j)*lattice_Sslip(:,:,1,index_myFamily+i,ph) - - !* Calculation of the tangent of Lp - forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLp_dTstar3333(k,l,m,n) = & - dLp_dTstar3333(k,l,m,n) + dgdot_dtauslip(j)*& - lattice_Sslip(k,l,1,index_myFamily+i,ph)*& - lattice_Sslip(m,n,1,index_myFamily+i,ph) - enddo slipSystemsLoop - enddo slipFamiliesLoop - -!-------------------------------------------------------------------------------------------------- -! correct Lp and dLp_dTstar3333 for twinned and transformed fraction - !* Total twin volume fraction - sumf = sum(state(instance)%twinFraction(1_pInt:nt,of)) ! safe for nt == 0 + BoltzmannRatio = prm%sbQedge/(kB*Temperature) + call math_eigenValuesVectorsSym(Mp,eigValues,eigVectors,error) - !* Total transformed volume fraction - sumftr = sum(state(instance)%stressTransFraction(1_pInt:nr,of)) + & - sum(state(instance)%strainTransFraction(1_pInt:nr,of)) - Lp = Lp * (1.0_pReal - sumf - sumftr) - dLp_dTstar3333 = dLp_dTstar3333 * (1.0_pReal - sumf - sumftr) - -!-------------------------------------------------------------------------------------------------- -! Shear banding (shearband) part - if(dNeq0(plastic_dislotwin_sbVelocity(instance)) .and. dNeq0(plastic_dislotwin_sbResistance(instance))) then - gdot_sb = 0.0_pReal - dgdot_dtausb = 0.0_pReal - call math_eigenValuesVectorsSym(math_Mandel6to33(Tstar_v),eigValues,eigVectors,error) - do j = 1_pInt,6_pInt - sb_s = 0.5_pReal*sqrt(2.0_pReal)*math_mul33x3(eigVectors,sb_sComposition(1:3,j)) - sb_m = 0.5_pReal*sqrt(2.0_pReal)*math_mul33x3(eigVectors,sb_mComposition(1:3,j)) - sb_Smatrix = math_tensorproduct33(sb_s,sb_m) - plastic_dislotwin_sbSv(1:6,j,ipc,ip,el) = math_Mandel33to6(math_symmetric33(sb_Smatrix)) + do i = 1_pInt,6_pInt + sb_s = 0.5_pReal*sqrt(2.0_pReal)*math_mul33x3(eigVectors,sb_sComposition(1:3,i)) + sb_m = 0.5_pReal*sqrt(2.0_pReal)*math_mul33x3(eigVectors,sb_mComposition(1:3,i)) + Schmid_shearBand = math_tensorproduct33(sb_s,sb_m) + tau = math_mul33xx33(Mp,Schmid_shearBand) - !* Calculation of Lp - !* Resolved shear stress on shear banding system - tau_sb(j) = dot_product(Tstar_v,plastic_dislotwin_sbSv(1:6,j,ipc,ip,el)) - - !* Stress ratios - if (abs(tau_sb(j)) < tol_math_check) then - StressRatio_p = 0.0_pReal - StressRatio_pminus1 = 0.0_pReal - else - StressRatio_p = (abs(tau_sb(j))/plastic_dislotwin_sbResistance(instance))& - **plastic_dislotwin_pShearBand(instance) - StressRatio_pminus1 = (abs(tau_sb(j))/plastic_dislotwin_sbResistance(instance))& - **(plastic_dislotwin_pShearBand(instance)-1.0_pReal) - endif - - !* Boltzmann ratio - BoltzmannRatio = plastic_dislotwin_sbQedge(instance)/(kB*Temperature) - !* Initial shear rates - DotGamma0 = plastic_dislotwin_sbVelocity(instance) + significantShearBandStress: if (abs(tau) > tol_math_check) then + StressRatio_p = (abs(tau)/prm%sbResistance)**prm%pShearBand + gdot_sb = sign(prm%sbVelocity*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**prm%qShearBand), tau) + dgdot_dtau = ((abs(gdot_sb)*BoltzmannRatio* prm%pShearBand*prm%qShearBand)/ prm%sbResistance) & + * (abs(tau)/prm%sbResistance)**(prm%pShearBand-1.0_pReal) & + * (1.0_pReal-StressRatio_p)**(prm%qShearBand-1.0_pReal) - !* Shear rates due to shearband - gdot_sb(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**& - plastic_dislotwin_qShearBand(instance))*sign(1.0_pReal,tau_sb(j)) - - !* Derivatives of shear rates - dgdot_dtausb(j) = & - ((abs(gdot_sb(j))*BoltzmannRatio*& - plastic_dislotwin_pShearBand(instance)*plastic_dislotwin_qShearBand(instance))/& - plastic_dislotwin_sbResistance(instance))*& - StressRatio_pminus1*(1_pInt-StressRatio_p)**(plastic_dislotwin_qShearBand(instance)-1.0_pReal) - - !* Plastic velocity gradient for shear banding - Lp = Lp + gdot_sb(j)*sb_Smatrix - - !* Calculation of the tangent of Lp - forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLp_dTstar3333(k,l,m,n) = & - dLp_dTstar3333(k,l,m,n) + dgdot_dtausb(j)*& - sb_Smatrix(k,l)*& - sb_Smatrix(m,n) + Lp = Lp + gdot_sb * Schmid_shearBand + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + dgdot_dtau * Schmid_shearBand(k,l) * Schmid_shearBand(m,n) + endif significantShearBandStress enddo - end if + + endif shearBandingContribution -!-------------------------------------------------------------------------------------------------- -! Mechanical twinning part - gdot_twin = 0.0_pReal - dgdot_dtautwin = 0.0_pReal - j = 0_pInt - twinFamiliesLoop: do f = 1_pInt,lattice_maxNtwinFamily - index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family - twinSystemsLoop: do i = 1_pInt,plastic_dislotwin_Ntwin(f,instance) - j = j+1_pInt - - !* Calculation of Lp - !* Resolved shear stress on twin system - tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,ph)) + call kinetics_twin(prm,stt,mse,of,Mp,temperature,gdot_slip,gdot_twin,dgdot_dtau_twin) + gdot_twin = f_unrotated * gdot_twin + dgdot_dtau_twin = f_unrotated * dgdot_dtau_twin + twinContibution: do i = 1_pInt, prm%totalNtwin + Lp = Lp + gdot_twin(i)*prm%Schmid_twin(1:3,1:3,i) + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + dgdot_dtau_twin(i)* prm%Schmid_twin(k,l,i)*prm%Schmid_twin(m,n,i) + enddo twinContibution - !* Stress ratios - if (tau_twin(j) > tol_math_check) then - StressRatio_r = (state(instance)%threshold_stress_twin(j,of)/tau_twin(j))**plastic_dislotwin_rPerTwinFamily(f,instance) - !* Shear rates and their derivatives due to twin - select case(lattice_structure(ph)) - case (LATTICE_fcc_ID) - s1=lattice_fcc_twinNucleationSlipPair(1,index_myFamily+i) - s2=lattice_fcc_twinNucleationSlipPair(2,index_myFamily+i) - if (tau_twin(j) < plastic_dislotwin_tau_r_twin(j,instance)) then - Ndot0_twin=(abs(gdot_slip(s1))*(state(instance)%rhoEdge(s2,of)+state(ph)%rhoEdgeDip(s2,of))+& !!!!! correct? - abs(gdot_slip(s2))*(state(instance)%rhoEdge(s1,of)+state(instance)%rhoEdgeDip(s1,of)))/& - (plastic_dislotwin_L0_twin(instance)*plastic_dislotwin_burgersPerSlipSystem(j,instance))*& - (1.0_pReal-exp(-plastic_dislotwin_VcrossSlip(instance)/(kB*Temperature)*& - (plastic_dislotwin_tau_r_twin(j,instance)-tau_twin(j)))) - else - Ndot0_twin=0.0_pReal - end if - case default - Ndot0_twin=plastic_dislotwin_Ndot0PerTwinSystem(j,instance) - end select - gdot_twin(j) = & - (1.0_pReal-sumf-sumftr)*lattice_shearTwin(index_myFamily+i,ph)*& - state(instance)%twinVolume(j,of)*Ndot0_twin*exp(-StressRatio_r) - dgdot_dtautwin(j) = ((gdot_twin(j)*plastic_dislotwin_rPerTwinFamily(f,instance))/tau_twin(j))*StressRatio_r - endif - - !* Plastic velocity gradient for mechanical twinning - Lp = Lp + gdot_twin(j)*lattice_Stwin(:,:,index_myFamily+i,ph) - - !* Calculation of the tangent of Lp - forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLp_dTstar3333(k,l,m,n) = & - dLp_dTstar3333(k,l,m,n) + dgdot_dtautwin(j)*& - lattice_Stwin(k,l,index_myFamily+i,ph)*& - lattice_Stwin(m,n,index_myFamily+i,ph) - enddo twinSystemsLoop - enddo twinFamiliesLoop + transConstribution: do i = 1_pInt, prm%totalNtrans - !* Phase transformation part - gdot_trans = 0.0_pReal - dgdot_dtautrans = 0.0_pReal - j = 0_pInt - transFamiliesLoop: do f = 1_pInt,lattice_maxNtransFamily - index_myFamily = sum(lattice_NtransSystem(1:f-1_pInt,ph)) ! at which index starts my family - transSystemsLoop: do i = 1_pInt,plastic_dislotwin_Ntrans(f,instance) - j = j+1_pInt + tau = math_mul33xx33(Mp,prm%Schmid_trans(1:3,1:3,i)) - !* Resolved shear stress on transformation system - tau_trans(j) = dot_product(Tstar_v,lattice_Strans_v(:,index_myFamily+i,ph)) + significantTransStress: if (tau > tol_math_check) then + StressRatio_s = (mse%threshold_stress_trans(i,of)/tau)**prm%s(i) - !* Stress ratios - if (tau_trans(j) > tol_math_check) then - StressRatio_s = (state(instance)%threshold_stress_trans(j,of)/tau_trans(j))**plastic_dislotwin_sPerTransFamily(f,instance) - !* Shear rates and their derivatives due to transformation - select case(lattice_structure(ph)) - case (LATTICE_fcc_ID) - s1=lattice_fcc_twinNucleationSlipPair(1,index_myFamily+i) - s2=lattice_fcc_twinNucleationSlipPair(2,index_myFamily+i) - if (tau_trans(j) < plastic_dislotwin_tau_r_trans(j,instance)) then - Ndot0_trans=(abs(gdot_slip(s1))*(state(instance)%rhoEdge(s2,of)+state(instance)%rhoEdgeDip(s2,of))+& !!!!! correct? - abs(gdot_slip(s2))*(state(instance)%rhoEdge(s1,of)+state(instance)%rhoEdgeDip(s1,of)))/& - (plastic_dislotwin_L0_trans(instance)*plastic_dislotwin_burgersPerSlipSystem(j,instance))*& - (1.0_pReal-exp(-plastic_dislotwin_VcrossSlip(instance)/(kB*Temperature)*& - (plastic_dislotwin_tau_r_trans(j,instance)-tau_trans(j)))) - else - Ndot0_trans=0.0_pReal - end if - case default - Ndot0_trans=plastic_dislotwin_Ndot0PerTransSystem(j,instance) - end select - gdot_trans(j) = & - (1.0_pReal-sumf-sumftr)*& - state(instance)%martensiteVolume(j,of)*Ndot0_trans*exp(-StressRatio_s) - dgdot_dtautrans(j) = ((gdot_trans(j)*plastic_dislotwin_sPerTransFamily(f,instance))/tau_trans(j))*StressRatio_s - endif + isFCCtrans: if (prm%isFCC) then + s1=prm%fcc_twinNucleationSlipPair(1,i) + s2=prm%fcc_twinNucleationSlipPair(2,i) + if (tau < mse%tau_r_trans(i,of)) then + Ndot0_trans=(abs(gdot_slip(s1))*(stt%rhoEdge(s2,of)+stt%rhoEdgeDip(s2,of))+& !!!!! correct? + abs(gdot_slip(s2))*(stt%rhoEdge(s1,of)+stt%rhoEdgeDip(s1,of)))/& + (prm%L0_trans*prm%burgers_slip(i))*& + (1.0_pReal-exp(-prm%VcrossSlip/(kB*Temperature)*(mse%tau_r_trans(i,of)-tau))) + else + Ndot0_trans=0.0_pReal + end if + else isFCCtrans + Ndot0_trans=prm%Ndot0_trans(i) + endif isFCCtrans - !* Plastic velocity gradient for phase transformation - Lp = Lp + gdot_trans(j)*lattice_Strans(:,:,index_myFamily+i,ph) + gdot_trans = mse%martensiteVolume(i,of) * Ndot0_trans*exp(-StressRatio_s) + gdot_trans = f_unrotated * gdot_trans + dgdot_dtau = ((gdot_trans*prm%s(i))/tau)*StressRatio_s + Lp = Lp + gdot_trans*prm%Schmid_trans(1:3,1:3,i) + + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + dgdot_dtau * prm%Schmid_trans(k,l,i)* prm%Schmid_trans(m,n,i) + endif significantTransStress - !* Calculation of the tangent of Lp - forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLp_dTstar3333(k,l,m,n) = & - dLp_dTstar3333(k,l,m,n) + dgdot_dtautrans(j)*& - lattice_Strans(k,l,index_myFamily+i,ph)*& - lattice_Strans(m,n,index_myFamily+i,ph) + enddo transConstribution - enddo transSystemsLoop - enddo transFamiliesLoop - - dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) + end associate end subroutine plastic_dislotwin_LpAndItsTangent @@ -1947,598 +1103,596 @@ end subroutine plastic_dislotwin_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el) +subroutine plastic_dislotwin_dotState(Mp,Temperature,instance,of) use prec, only: & tol_math_check, & dEq0 use math, only: & + math_mul33xx33, & + math_Mandel6to33, & pi use material, only: & material_phase, & phase_plasticityInstance, & plasticState, & - phaseAt, phasememberAt - use lattice, only: & - lattice_Sslip_v, & - lattice_Stwin_v, & - lattice_Strans_v, & - lattice_maxNslipFamily, & - lattice_maxNtwinFamily, & - lattice_maxNtransFamily, & - lattice_NslipSystem, & - lattice_NtwinSystem, & - lattice_NtransSystem, & - lattice_sheartwin, & - lattice_mu, & - lattice_structure, & - lattice_fcc_twinNucleationSlipPair, & - lattice_fccTobcc_transNucleationTwinPair, & - lattice_fccTobcc_shearCritTrans, & - LATTICE_fcc_ID + phasememberAt implicit none - real(pReal), dimension(6), intent(in):: & - Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal), intent(in) :: & + real(pReal), dimension(3,3), intent(in):: & + Mp !< Mandel stress + real(pReal), intent(in) :: & temperature !< temperature at integration point - integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element + integer(pInt), intent(in) :: & + instance, & + of - integer(pInt) :: instance,ns,nt,nr,f,i,j,index_myFamily,s1,s2, & - ph, & - of - real(pReal) :: sumf,sumftr,StressRatio_p,StressRatio_pminus1,BoltzmannRatio,DotGamma0,& + integer(pInt) :: i,s1,s2 + real(pReal) :: f_unrotated,StressRatio_p,BoltzmannRatio,& EdgeDipMinDistance,AtomicVolume,VacancyDiffusion,StressRatio_r,Ndot0_twin,stressRatio,& Ndot0_trans,StressRatio_s,EdgeDipDistance, ClimbVelocity,DotRhoEdgeDipClimb,DotRhoEdgeDipAnnihilation, & - DotRhoDipFormation,DotRhoMultiplication,DotRhoEdgeEdgeAnnihilation - real(pReal), dimension(plastic_dislotwin_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & - gdot_slip,tau_slip + DotRhoDipFormation,DotRhoMultiplication,DotRhoEdgeEdgeAnnihilation, & + tau + real(pReal), dimension(plasticState(instance)%Nslip) :: & + gdot_slip - real(pReal), dimension(plastic_dislotwin_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & - tau_twin - real(pReal), dimension(plastic_dislotwin_totalNtrans(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & - tau_trans - !* Shortened notation - of = phasememberAt(ipc,ip,el) - ph = phaseAt(ipc,ip,el) - instance = phase_plasticityInstance(ph) - ns = plastic_dislotwin_totalNslip(instance) - nt = plastic_dislotwin_totalNtwin(instance) - nr = plastic_dislotwin_totalNtrans(instance) + type(tParameters) :: prm + type(tDislotwinState) :: stt, dot + type(tDislotwinMicrostructure) :: mse - !* Total twin volume fraction - sumf = sum(state(instance)%twinFraction(1_pInt:nt,of)) ! safe for nt == 0 - plasticState(ph)%dotState(:,of) = 0.0_pReal - - !* Total transformed volume fraction - sumftr = sum(state(instance)%stressTransFraction(1_pInt:nr,of)) + & - sum(state(instance)%strainTransFraction(1_pInt:nr,of)) - - !* Dislocation density evolution - gdot_slip = 0.0_pReal - j = 0_pInt - do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families - index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family - do i = 1_pInt,plastic_dislotwin_Nslip(f,instance) ! process each (active) slip system in family - j = j+1_pInt - - !* Resolved shear stress on slip system - tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph)) - if((abs(tau_slip(j))-state(instance)%threshold_stress_slip(j,of)) > tol_math_check) then - !* Stress ratios - stressRatio =((abs(tau_slip(j))- state(instance)%threshold_stress_slip(j,of))/& - (plastic_dislotwin_SolidSolutionStrength(instance)+plastic_dislotwin_tau_peierlsPerSlipFamily(f,instance))) - StressRatio_p = stressRatio** plastic_dislotwin_pPerSlipFamily(f,instance) - StressRatio_pminus1 = stressRatio**(plastic_dislotwin_pPerSlipFamily(f,instance)-1.0_pReal) - !* Boltzmann ratio - BoltzmannRatio = plastic_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature) - !* Initial shear rates - DotGamma0 = & - plasticState(ph)%state(j, of)*plastic_dislotwin_burgersPerSlipSystem(j,instance)*& - plastic_dislotwin_v0PerSlipSystem(j,instance) - - !* Shear rates due to slip - gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)** & - plastic_dislotwin_qPerSlipFamily(f,instance))*sign(1.0_pReal,tau_slip(j)) - endif - !* Multiplication - DotRhoMultiplication = abs(gdot_slip(j))/& - (plastic_dislotwin_burgersPerSlipSystem(j,instance)*state(instance)%mfp_slip(j,of)) - !* Dipole formation - EdgeDipMinDistance = & - plastic_dislotwin_CEdgeDipMinDistance(instance)*plastic_dislotwin_burgersPerSlipSystem(j,instance) - if (dEq0(tau_slip(j))) then - DotRhoDipFormation = 0.0_pReal - else - EdgeDipDistance = & - (3.0_pReal*lattice_mu(ph)*plastic_dislotwin_burgersPerSlipSystem(j,instance))/& - (16.0_pReal*pi*abs(tau_slip(j))) - if (EdgeDipDistance>state(instance)%mfp_slip(j,of)) EdgeDipDistance=state(instance)%mfp_slip(j,of) - if (EdgeDipDistance tol_math_check) then - StressRatio_r = (state(instance)%threshold_stress_twin(j,of)/& - tau_twin(j))**plastic_dislotwin_rPerTwinFamily(f,instance) - !* Shear rates and their derivatives due to twin - select case(lattice_structure(ph)) - case (LATTICE_fcc_ID) - s1=lattice_fcc_twinNucleationSlipPair(1,index_myFamily+i) - s2=lattice_fcc_twinNucleationSlipPair(2,index_myFamily+i) - if (tau_twin(j) < plastic_dislotwin_tau_r_twin(j,instance)) then - Ndot0_twin=(abs(gdot_slip(s1))*(state(instance)%rhoEdge(s2,of)+state(instance)%rhoEdgeDip(s2,of))+& - abs(gdot_slip(s2))*(state(instance)%rhoEdge(s1,of)+state(instance)%rhoEdgeDip(s1,of)))/& - (plastic_dislotwin_L0_twin(instance)*plastic_dislotwin_burgersPerSlipSystem(j,instance))*& - (1.0_pReal-exp(-plastic_dislotwin_VcrossSlip(instance)/(kB*Temperature)*& - (plastic_dislotwin_tau_r_twin(j,instance)-tau_twin(j)))) - else - Ndot0_twin=0.0_pReal - end if - case default - Ndot0_twin=plastic_dislotwin_Ndot0PerTwinSystem(j,instance) - end select - dotState(instance)%twinFraction(j,of) = & - (1.0_pReal-sumf-sumftr)*& - state(instance)%twinVolume(j,of)*Ndot0_twin*exp(-StressRatio_r) - !* Dotstate for accumulated shear due to twin - dotState(instance)%accshear_twin(j,of) = dotState(instance)%twinFraction(j,of) * & - lattice_sheartwin(index_myfamily+i,ph) - endif - enddo - enddo + associate(prm => param(instance), stt => state(instance), & + dot => dotstate(instance), mse => microstructure(instance)) - !* Transformation volume fraction evolution - j = 0_pInt - do f = 1_pInt,lattice_maxNtransFamily ! loop over all trans families - index_myFamily = sum(lattice_NtransSystem(1:f-1_pInt,ph)) ! at which index starts my family - do i = 1_pInt,plastic_dislotwin_Ntrans(f,instance) ! process each (active) trans system in family - j = j+1_pInt + dot%whole(:,of) = 0.0_pReal - !* Resolved shear stress on transformation system - tau_trans(j) = dot_product(Tstar_v,lattice_Strans_v(:,index_myFamily+i,ph)) + f_unrotated = 1.0_pReal & + - sum(stt%twinFraction(1_pInt:prm%totalNtwin,of)) & + - sum(stt%stressTransFraction(1_pInt:prm%totalNtrans,of)) & + - sum(stt%strainTransFraction(1_pInt:prm%totalNtrans,of)) - !* Stress ratios - if (tau_trans(j) > tol_math_check) then - StressRatio_s = (state(instance)%threshold_stress_trans(j,of)/& - tau_trans(j))**plastic_dislotwin_sPerTransFamily(f,instance) - !* Shear rates and their derivatives due to transformation - select case(lattice_structure(ph)) - case (LATTICE_fcc_ID) - s1=lattice_fcc_twinNucleationSlipPair(1,index_myFamily+i) - s2=lattice_fcc_twinNucleationSlipPair(2,index_myFamily+i) - if (tau_trans(j) < plastic_dislotwin_tau_r_trans(j,instance)) then - Ndot0_trans=(abs(gdot_slip(s1))*(state(instance)%rhoEdge(s2,of)+state(instance)%rhoEdgeDip(s2,of))+& - abs(gdot_slip(s2))*(state(instance)%rhoEdge(s1,of)+state(instance)%rhoEdgeDip(s1,of)))/& - (plastic_dislotwin_L0_trans(instance)*plastic_dislotwin_burgersPerSlipSystem(j,instance))*& - (1.0_pReal-exp(-plastic_dislotwin_VcrossSlip(instance)/(kB*Temperature)*& - (plastic_dislotwin_tau_r_trans(j,instance)-tau_trans(j)))) - else - Ndot0_trans=0.0_pReal - end if - case default - Ndot0_trans=plastic_dislotwin_Ndot0PerTransSystem(j,instance) - end select - dotState(instance)%strainTransFraction(j,of) = & - (1.0_pReal-sumf-sumftr)*& - state(instance)%martensiteVolume(j,of)*Ndot0_trans*exp(-StressRatio_s) + call kinetics_slip(prm,stt,mse,of,Mp,temperature,gdot_slip) + slipState: do i = 1_pInt, prm%totalNslip + tau = math_mul33xx33(Mp,prm%Schmid_slip(1:3,1:3,i)) + + DotRhoMultiplication = abs(gdot_slip(i))/(prm%burgers_slip(i)*mse%mfp_slip(i,of)) + EdgeDipMinDistance = prm%CEdgeDipMinDistance*prm%burgers_slip(i) + + significantSlipStress2: if (dEq0(tau)) then + DotRhoDipFormation = 0.0_pReal + else significantSlipStress2 + EdgeDipDistance = (3.0_pReal*prm%mu*prm%burgers_slip(i))/(16.0_pReal*PI*abs(tau)) + if (EdgeDipDistance>mse%mfp_slip(i,of)) EdgeDipDistance = mse%mfp_slip(i,of) + if (EdgeDipDistance tol_math_check) then + StressRatio_r = (mse%threshold_stress_twin(i,of)/tau)**prm%r(i) + isFCCtwin: if (prm%isFCC) then + s1=prm%fcc_twinNucleationSlipPair(1,i) + s2=prm%fcc_twinNucleationSlipPair(2,i) + if (tau < mse%tau_r_twin(i,of)) then + Ndot0_twin=(abs(gdot_slip(s1))*(stt%rhoEdge(s2,of)+stt%rhoEdgeDip(s2,of))+& + abs(gdot_slip(s2))*(stt%rhoEdge(s1,of)+stt%rhoEdgeDip(s1,of)))/& + (prm%L0_twin*prm%burgers_slip(i))*(1.0_pReal-exp(-prm%VcrossSlip/(kB*Temperature)*& + (mse%tau_r_twin(i,of)-tau))) + else + Ndot0_twin=0.0_pReal + end if + else isFCCtwin + Ndot0_twin=prm%Ndot0_twin(i) + endif isFCCtwin + dot%twinFraction(i,of) = f_unrotated * mse%twinVolume(i,of)*Ndot0_twin*exp(-StressRatio_r) + dot%accshear_twin(i,of) = dot%twinFraction(i,of) * prm%shear_twin(i) + endif significantTwinStress + + enddo twinState + + transState: do i = 1_pInt, prm%totalNtrans + + tau = math_mul33xx33(Mp,prm%Schmid_trans(1:3,1:3,i)) + + significantTransStress: if (tau > tol_math_check) then + StressRatio_s = (mse%threshold_stress_trans(i,of)/tau)**prm%s(i) + isFCCtrans: if (prm%isFCC) then + s1=prm%fcc_twinNucleationSlipPair(1,i) + s2=prm%fcc_twinNucleationSlipPair(2,i) + if (tau < mse%tau_r_trans(i,of)) then + Ndot0_trans=(abs(gdot_slip(s1))*(stt%rhoEdge(s2,of)+stt%rhoEdgeDip(s2,of))+& + abs(gdot_slip(s2))*(stt%rhoEdge(s1,of)+stt%rhoEdgeDip(s1,of)))/& + (prm%L0_trans*prm%burgers_slip(i))*(1.0_pReal-exp(-prm%VcrossSlip/(kB*Temperature)*& + (mse%tau_r_trans(i,of)-tau))) + else + Ndot0_trans=0.0_pReal + end if + else isFCCtrans + Ndot0_trans=prm%Ndot0_trans(i) + endif isFCCtrans + dot%strainTransFraction(i,of) = f_unrotated * & + mse%martensiteVolume(i,of)*Ndot0_trans*exp(-StressRatio_s) !* Dotstate for accumulated shear due to transformation - !dotState(instance)%accshear_trans(j,of) = dotState(instance)%strainTransFraction(j,of) * & + !dot%accshear_trans(i,of) = dot%strainTransFraction(i,of) * & ! lattice_sheartrans(index_myfamily+i,ph) - endif - - enddo - enddo + endif significantTransStress + + enddo transState + end associate end subroutine plastic_dislotwin_dotState + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates shear rates on slip systems +!-------------------------------------------------------------------------------------------------- +pure subroutine kinetics_slip(prm,stt,mse,of,Mp,temperature,gdot_slip,dgdot_dtau_slip) + use prec, only: & + tol_math_check, & + dNeq0 + use math, only: & + math_mul33xx33 + + implicit none + type(tParameters), intent(in) :: & + prm + type(tDislotwinState), intent(in) :: & + stt + integer(pInt), intent(in) :: & + of + type(tDislotwinMicrostructure), intent(in) :: & + mse + real(pReal), dimension(prm%totalNslip), intent(out) :: & + gdot_slip + real(pReal), dimension(prm%totalNslip), optional, intent(out) :: & + dgdot_dtau_slip + real(pReal), dimension(prm%totalNslip) :: & + dgdot_dtau + real(pReal), dimension(3,3), intent(in) :: & + Mp + real(pReal), intent(in) :: & + temperature + + real, dimension(prm%totalNslip) :: & + tau, & + stressRatio, & + StressRatio_p, & + BoltzmannRatio, & + v_wait_inverse, & !< inverse of the effective velocity of a dislocation waiting at obstacles (unsigned) + v_run_inverse, & !< inverse of the velocity of a free moving dislocation (unsigned) + dV_wait_inverse_dTau, & + dV_run_inverse_dTau, & + dV_dTau, & + tau_eff !< effective resolved stress + integer(pInt) :: i + + do i = 1_pInt, prm%totalNslip + tau(i) = math_mul33xx33(Mp,prm%Schmid_slip(1:3,1:3,i)) + enddo + tau_eff = abs(tau)-mse%threshold_stress_slip(:,of) + + significantStress: where(tau_eff > tol_math_check) + stressRatio = tau_eff/(prm%SolidSolutionStrength+prm%tau_peierls) + StressRatio_p = stressRatio** prm%p + BoltzmannRatio = prm%Qedge/(kB*Temperature) + v_wait_inverse = prm%v0**(-1.0_pReal) * exp(BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q) + v_run_inverse = prm%B/(tau_eff*prm%burgers_slip) + + gdot_slip = sign(stt%rhoEdge(:,of)*prm%burgers_slip/(v_wait_inverse+v_run_inverse),tau) + + dV_wait_inverse_dTau = v_wait_inverse * prm%p * prm%q * BoltzmannRatio & + * (stressRatio**(prm%p-1.0_pReal)) & + * (1.0_pReal-StressRatio_p)**(prm%q-1.0_pReal) & + / (prm%SolidSolutionStrength+prm%tau_peierls) + dV_run_inverse_dTau = v_run_inverse/tau_eff + dV_dTau = (dV_wait_inverse_dTau+dV_run_inverse_dTau) & + / (v_wait_inverse+v_run_inverse)**2.0_pReal + dgdot_dtau = dV_dTau*stt%rhoEdge(:,of)*prm%burgers_slip + else where significantStress + gdot_slip = 0.0_pReal + dgdot_dtau = 0.0_pReal + end where significantStress + if(present(dgdot_dtau_slip)) dgdot_dtau_slip = dgdot_dtau + +end subroutine kinetics_slip + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates shear rates on twin systems +!-------------------------------------------------------------------------------------------------- +pure subroutine kinetics_twin(prm,stt,mse,of,Mp,temperature,gdot_slip,gdot_twin,dgdot_dtau_twin) + use prec, only: & + tol_math_check, & + dNeq0 + use math, only: & + math_mul33xx33 + + implicit none + type(tParameters), intent(in) :: & + prm + type(tDislotwinState), intent(in) :: & + stt + integer(pInt), intent(in) :: & + of + type(tDislotwinMicrostructure), intent(in) :: & + mse + real(pReal), dimension(prm%totalNslip), intent(out) :: & + gdot_slip + real(pReal), dimension(prm%totalNtwin), intent(out) :: & + gdot_twin + real(pReal), dimension(prm%totalNtwin), optional, intent(out) :: & + dgdot_dtau_twin + real(pReal), dimension(3,3), intent(in) :: & + Mp + real(pReal), intent(in) :: & + temperature + + real, dimension(prm%totalNtwin) :: & + tau, & + Ndot0_twin, & + stressRatio_r, & + dgdot_dtau + + integer(pInt) :: i,s1,s2 + + do i = 1_pInt, prm%totalNtwin + tau(i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i)) + isFCC: if (prm%isFCC) then + s1=prm%fcc_twinNucleationSlipPair(1,i) + s2=prm%fcc_twinNucleationSlipPair(2,i) + if (tau(i) < mse%tau_r_twin(i,of)) then + Ndot0_twin=(abs(gdot_slip(s1))*(stt%rhoEdge(s2,of)+stt%rhoEdgeDip(s2,of))+& + abs(gdot_slip(s2))*(stt%rhoEdge(s1,of)+stt%rhoEdgeDip(s1,of)))/& + (prm%L0_twin*prm%burgers_slip(i))*& + (1.0_pReal-exp(-prm%VcrossSlip/(kB*Temperature)*& + (mse%tau_r_twin(i,of)-tau))) + else + Ndot0_twin=0.0_pReal + end if + else isFCC + Ndot0_twin=prm%Ndot0_twin(i) + endif isFCC + enddo + + + significantStress: where(tau > tol_math_check) + StressRatio_r = (mse%threshold_stress_twin(:,of)/tau)**prm%r + gdot_twin = prm%shear_twin * mse%twinVolume(:,of) * Ndot0_twin*exp(-StressRatio_r) + dgdot_dtau = ((gdot_twin*prm%r)/tau)*StressRatio_r + else where significantStress + gdot_twin = 0.0_pReal + dgdot_dtau = 0.0_pReal + end where significantStress + + if(present(dgdot_dtau_twin)) dgdot_dtau_twin = dgdot_dtau + +end subroutine kinetics_twin + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates shear rates on transformation systems +!-------------------------------------------------------------------------------------------------- +pure subroutine kinetics_trans(prm,stt,mse,of,Mp,temperature,gdot_slip,gdot_trans,dgdot_dtau_trans) + use prec, only: & + tol_math_check, & + dNeq0 + use math, only: & + math_mul33xx33 + + implicit none + type(tParameters), intent(in) :: & + prm + type(tDislotwinState), intent(in) :: & + stt + integer(pInt), intent(in) :: & + of + type(tDislotwinMicrostructure), intent(in) :: & + mse + real(pReal), dimension(prm%totalNslip), intent(out) :: & + gdot_slip + real(pReal), dimension(prm%totalNtrans), intent(out) :: & + gdot_trans + real(pReal), dimension(prm%totalNtrans), optional, intent(out) :: & + dgdot_dtau_trans + real(pReal), dimension(3,3), intent(in) :: & + Mp + real(pReal), intent(in) :: & + temperature + + real, dimension(prm%totalNtrans) :: & + tau, & + Ndot0_trans, & + stressRatio_r, & + dgdot_dtau + + integer(pInt) :: i,s1,s2 + + do i = 1_pInt, prm%totalNtrans + tau(i) = math_mul33xx33(Mp,prm%Schmid_trans(1:3,1:3,i)) + isFCC: if (prm%isFCC) then + s1=prm%fcc_twinNucleationSlipPair(1,i) + s2=prm%fcc_twinNucleationSlipPair(2,i) + if (tau(i) < mse%tau_r_trans(i,of)) then + Ndot0_trans=(abs(gdot_slip(s1))*(stt%rhoEdge(s2,of)+stt%rhoEdgeDip(s2,of))+& + abs(gdot_slip(s2))*(stt%rhoEdge(s1,of)+stt%rhoEdgeDip(s1,of)))/& + (prm%L0_trans*prm%burgers_slip(i))*& ! burgers_slip correct? + (1.0_pReal-exp(-prm%VcrossSlip/(kB*Temperature)*& + (mse%tau_r_trans(i,of)-tau))) + else + Ndot0_trans=0.0_pReal + end if + else isFCC + Ndot0_trans=prm%Ndot0_trans(i) + endif isFCC + enddo +! +! +! endif isFCCtrans +! dot%strainTransFraction(i,of) = f_unrotated * & +! mse%martensiteVolume(i,of)*Ndot0_trans*exp(-StressRatio_s) +! !* Dotstate for accumulated shear due to transformation +! !dot%accshear_trans(i,of) = dot%strainTransFraction(i,of) * & +! ! lattice_sheartrans(index_myfamily+i,ph) +! endif significantTransStress +! +! enddo transState +! +! +! significantStress: where(tau > tol_math_check) +! StressRatio_r = (mse%threshold_stress_twin(:,of)/tau)**prm%r +! gdot_twin = prm%shear_twin * mse%twinVolume(:,of) * Ndot0_twin*exp(-StressRatio_r) +! dgdot_dtau = ((gdot_twin*prm%r)/tau)*StressRatio_r +! else where significantStress +! gdot_twin = 0.0_pReal +! dgdot_dtau = 0.0_pReal +! end where significantStress +! +! if(present(dgdot_dtau_twin)) dgdot_dtau_twin = dgdot_dtau +! +end subroutine kinetics_trans + !-------------------------------------------------------------------------------------------------- !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- -function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el) +function plastic_dislotwin_postResults(Mp,Temperature,instance,of) result(postResults) use prec, only: & tol_math_check, & dEq0 use math, only: & - pi, & - math_Mandel6to33, & - math_eigenValuesSym33, & - math_eigenValuesVectorsSym33 + PI, & + math_mul33xx33, & + math_Mandel6to33 use material, only: & material_phase, & + plasticState, & phase_plasticityInstance,& - phaseAt, phasememberAt - use lattice, only: & - lattice_Sslip_v, & - lattice_Stwin_v, & - lattice_maxNslipFamily, & - lattice_maxNtwinFamily, & - lattice_NslipSystem, & - lattice_NtwinSystem, & - lattice_shearTwin, & - lattice_mu, & - lattice_structure, & - lattice_fcc_twinNucleationSlipPair, & - LATTICE_fcc_ID + phasememberAt implicit none - real(pReal), dimension(6), intent(in) :: & - Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + real(pReal), dimension(3,3),intent(in) :: & + Mp !< 2nd Piola Kirchhoff stress tensor in Mandel notation real(pReal), intent(in) :: & temperature !< temperature at integration point integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - - real(pReal), dimension(plastic_dislotwin_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & - plastic_dislotwin_postResults - integer(pInt) :: & - instance,& - ns,nt,nr,& - f,o,i,c,j,index_myFamily,& - s1,s2, & - ph, & + instance, & of - real(pReal) :: sumf,tau,StressRatio_p,StressRatio_pminus1,BoltzmannRatio,DotGamma0,StressRatio_r,Ndot0_twin,dgdot_dtauslip, & + + real(pReal), dimension(sum(plastic_dislotwin_sizePostResult(:,instance))) :: & + postResults + + integer(pInt) :: & + o,c,j,& + s1,s2 + real(pReal) :: sumf_twin,tau,StressRatio_p,StressRatio_pminus1,BoltzmannRatio,DotGamma0,StressRatio_r,Ndot0_twin,dgdot_dtauslip, & stressRatio - real(preal), dimension(plastic_dislotwin_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + real(pReal), dimension(param(instance)%totalNslip) :: & gdot_slip - real(pReal), dimension(3,3) :: eigVectors - real(pReal), dimension (3) :: eigValues - !* Shortened notation - of = phasememberAt(ipc,ip,el) - ph = phaseAt(ipc,ip,el) - instance = phase_plasticityInstance(ph) - ns = plastic_dislotwin_totalNslip(instance) - nt = plastic_dislotwin_totalNtwin(instance) - nr = plastic_dislotwin_totalNtrans(instance) + type(tParameters) :: prm + type(tDislotwinState) :: stt + type(tDislotwinMicrostructure) :: mse + - !* Total twin volume fraction - sumf = sum(state(instance)%twinFraction(1_pInt:nt,of)) ! safe for nt == 0 + associate(prm => param(instance), stt => state(instance), mse => microstructure(instance)) + + sumf_twin = sum(stt%twinFraction(1_pInt:prm%totalNtwin,of)) - !* Required output c = 0_pInt - plastic_dislotwin_postResults = 0.0_pReal - do o = 1_pInt,plastic_dislotwin_Noutput(instance) - select case(plastic_dislotwin_outputID(o,instance)) + postResults = 0.0_pReal + do o = 1_pInt,size(prm%outputID) + select case(prm%outputID(o)) - case (edge_density_ID) - plastic_dislotwin_postResults(c+1_pInt:c+ns) = state(instance)%rhoEdge(1_pInt:ns,of) - c = c + ns - case (dipole_density_ID) - plastic_dislotwin_postResults(c+1_pInt:c+ns) = state(instance)%rhoEdgeDip(1_pInt:ns,of) - c = c + ns - case (shear_rate_slip_ID) - j = 0_pInt - do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families - index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family - do i = 1_pInt,plastic_dislotwin_Nslip(f,instance) ! process each (active) slip system in family - j = j + 1_pInt ! could be taken from state by now! + case (edge_density_ID) + postResults(c+1_pInt:c+prm%totalNslip) = stt%rhoEdge(1_pInt:prm%totalNslip,of) + c = c + prm%totalNslip + case (dipole_density_ID) + postResults(c+1_pInt:c+prm%totalNslip) = stt%rhoEdgeDip(1_pInt:prm%totalNslip,of) + c = c + prm%totalNslip + case (shear_rate_slip_ID) + call kinetics_slip(prm,stt,mse,of,Mp,temperature,postResults(c+1:c+prm%totalNslip)) + c = c + prm%totalNslip + case (accumulated_shear_slip_ID) + postResults(c+1_pInt:c+prm%totalNslip) = stt%accshear_slip(1_pInt:prm%totalNslip,of) + c = c + prm%totalNslip + case (mfp_slip_ID) + postResults(c+1_pInt:c+prm%totalNslip) = mse%mfp_slip(1_pInt:prm%totalNslip,of) + c = c + prm%totalNslip + case (resolved_stress_slip_ID) + do j = 1_pInt, prm%totalNslip + postResults(c+j) = math_mul33xx33(Mp,prm%Schmid_slip(1:3,1:3,j)) + enddo + c = c + prm%totalNslip + case (threshold_stress_slip_ID) + postResults(c+1_pInt:c+prm%totalNslip) = mse%threshold_stress_slip(1_pInt:prm%totalNslip,of) + c = c + prm%totalNslip + case (edge_dipole_distance_ID) + do j = 1_pInt, prm%totalNslip + postResults(c+j) = (3.0_pReal*prm%mu*prm%burgers_slip(j)) & + / (16.0_pReal*PI*abs(math_mul33xx33(Mp,prm%Schmid_slip(1:3,1:3,j)))) + postResults(c+j)=min(postResults(c+j),mse%mfp_slip(j,of)) + ! postResults(c+j)=max(postResults(c+j),& + ! plasticState(ph)%state(4*ns+2*nt+2*nr+j, of)) + enddo + c = c + prm%totalNslip + ! case (resolved_stress_shearband_ID) + ! do j = 1_pInt,6_pInt ! loop over all shearband families + ! postResults(c+j) = dot_product(Tstar_v,sbSv(1:6,j,ipc,ip,el)) + ! enddo + ! c = c + 6_pInt + ! case (shear_rate_shearband_ID) + ! do j = 1_pInt,6_pInt ! loop over all shearbands + ! tau = dot_product(Tstar_v,sbSv(1:6,j,ipc,ip,el)) + ! if (abs(tau) < tol_math_check) then + ! StressRatio_p = 0.0_pReal + ! StressRatio_pminus1 = 0.0_pReal + ! else + ! StressRatio_p = (abs(tau)/prm%sbResistance)**prm%pShearBand + ! StressRatio_pminus1 = (abs(tau)/prm%sbResistance)**(prm%pShearBand-1.0_pReal) + ! endif + ! BoltzmannRatio = prm%sbQedge/(kB*Temperature) + ! DotGamma0 = prm%sbVelocity + ! postResults(c+j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**prm%qShearBand)*& + ! sign(1.0_pReal,tau) + ! enddo + ! c = c + 6_pInt + case (twin_fraction_ID) + postResults(c+1_pInt:c+prm%totalNtwin) = stt%twinFraction(1_pInt:prm%totalNtwin,of) + c = c + prm%totalNtwin + case (shear_rate_twin_ID) + do j = 1_pInt, prm%totalNslip + tau = math_mul33xx33(Mp,prm%Schmid_slip(1:3,1:3,j)) + if((abs(tau)-mse%threshold_stress_slip(j,of)) > tol_math_check) then + StressRatio_p = ((abs(tau)-mse%threshold_stress_slip(j,of))/& + (prm%SolidSolutionStrength+& + prm%tau_peierls(j)))& + **prm%p(j) + StressRatio_pminus1 = ((abs(tau)-mse%threshold_stress_slip(j,of))/& + (prm%SolidSolutionStrength+& + prm%tau_peierls(j)))& + **(prm%p(j)-1.0_pReal) + BoltzmannRatio = prm%Qedge(j)/(kB*Temperature) + DotGamma0 = stt%rhoEdge(j,of)*prm%burgers_slip(j)* prm%v0(j) - !* Resolved shear stress on slip system - tau = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph)) - !* Stress ratios - if((abs(tau)-state(instance)%threshold_stress_slip(j,of)) > tol_math_check) then - !* Stress ratios - stressRatio = ((abs(tau)-state(ph)%threshold_stress_slip(j,of))/& - (plastic_dislotwin_SolidSolutionStrength(instance)+& - plastic_dislotwin_tau_peierlsPerSlipFamily(f,instance))) - StressRatio_p = stressRatio** plastic_dislotwin_pPerSlipFamily(f,instance) - StressRatio_pminus1 = stressRatio**(plastic_dislotwin_pPerSlipFamily(f,instance)-1.0_pReal) - !* Boltzmann ratio - BoltzmannRatio = plastic_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature) - !* Initial shear rates - DotGamma0 = & - state(instance)%rhoEdge(j,of)*plastic_dislotwin_burgersPerSlipSystem(j,instance)* & - plastic_dislotwin_v0PerSlipSystem(j,instance) - - !* Shear rates due to slip - plastic_dislotwin_postResults(c+j) = & - DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**& - plastic_dislotwin_qPerSlipFamily(f,instance))*sign(1.0_pReal,tau) + gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**& + prm%q(j))*sign(1.0_pReal,tau) + else + gdot_slip(j) = 0.0_pReal + endif + enddo + + do j = 1_pInt, prm%totalNtwin + tau = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,j)) + + if ( tau > 0.0_pReal ) then + isFCCtwin: if (prm%isFCC) then + s1=prm%fcc_twinNucleationSlipPair(1,j) + s2=prm%fcc_twinNucleationSlipPair(2,j) + if (tau < mse%tau_r_twin(j,of)) then + Ndot0_twin=(abs(gdot_slip(s1))*(stt%rhoEdge(s2,of)+stt%rhoEdgeDip(s2,of))+& + abs(gdot_slip(s2))*(stt%rhoEdge(s1,of)+stt%rhoEdgeDip(s1,of)))/& + (prm%L0_twin* prm%burgers_slip(j))*& + (1.0_pReal-exp(-prm%VcrossSlip/(kB*Temperature)* (mse%tau_r_twin(j,of)-tau))) else - plastic_dislotwin_postResults(c+j) = 0.0_pReal - endif - - enddo ; enddo - c = c + ns - case (accumulated_shear_slip_ID) - plastic_dislotwin_postResults(c+1_pInt:c+ns) = & - state(instance)%accshear_slip(1_pInt:ns,of) - c = c + ns - case (mfp_slip_ID) - plastic_dislotwin_postResults(c+1_pInt:c+ns) =& - state(instance)%mfp_slip(1_pInt:ns,of) - c = c + ns - case (resolved_stress_slip_ID) - j = 0_pInt - do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families - index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family - do i = 1_pInt,plastic_dislotwin_Nslip(f,instance) ! process each (active) slip system in family - j = j + 1_pInt - plastic_dislotwin_postResults(c+j) =& - dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph)) - enddo; enddo - c = c + ns - case (threshold_stress_slip_ID) - plastic_dislotwin_postResults(c+1_pInt:c+ns) = & - state(instance)%threshold_stress_slip(1_pInt:ns,of) - c = c + ns - case (edge_dipole_distance_ID) - j = 0_pInt - do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families - index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family - do i = 1_pInt,plastic_dislotwin_Nslip(f,instance) ! process each (active) slip system in family - j = j + 1_pInt - plastic_dislotwin_postResults(c+j) = & - (3.0_pReal*lattice_mu(ph)*plastic_dislotwin_burgersPerSlipSystem(j,instance))/& - (16.0_pReal*pi*abs(dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph)))) - plastic_dislotwin_postResults(c+j)=min(plastic_dislotwin_postResults(c+j),& - state(instance)%mfp_slip(j,of)) - ! plastic_dislotwin_postResults(c+j)=max(plastic_dislotwin_postResults(c+j),& - ! plasticState(ph)%state(4*ns+2*nt+2*nr+j, of)) - enddo; enddo - c = c + ns - case (resolved_stress_shearband_ID) - do j = 1_pInt,6_pInt ! loop over all shearband families - plastic_dislotwin_postResults(c+j) = dot_product(Tstar_v, & - plastic_dislotwin_sbSv(1:6,j,ipc,ip,el)) - enddo - c = c + 6_pInt - case (shear_rate_shearband_ID) - do j = 1_pInt,6_pInt ! loop over all shearbands - !* Resolved shear stress on shearband system - tau = dot_product(Tstar_v,plastic_dislotwin_sbSv(1:6,j,ipc,ip,el)) - !* Stress ratios - if (abs(tau) < tol_math_check) then - StressRatio_p = 0.0_pReal - StressRatio_pminus1 = 0.0_pReal - else - StressRatio_p = (abs(tau)/plastic_dislotwin_sbResistance(instance))**& - plastic_dislotwin_pShearBand(instance) - StressRatio_pminus1 = (abs(tau)/plastic_dislotwin_sbResistance(instance))**& - (plastic_dislotwin_pShearBand(instance)-1.0_pReal) - endif - !* Boltzmann ratio - BoltzmannRatio = plastic_dislotwin_sbQedge(instance)/(kB*Temperature) - !* Initial shear rates - DotGamma0 = plastic_dislotwin_sbVelocity(instance) - ! Shear rate due to shear band - plastic_dislotwin_postResults(c+j) = & - DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**plastic_dislotwin_qShearBand(instance))*& - sign(1.0_pReal,tau) - enddo - c = c + 6_pInt - case (twin_fraction_ID) - plastic_dislotwin_postResults(c+1_pInt:c+nt) = state(instance)%twinFraction(1_pInt:nt,of) - c = c + nt - case (shear_rate_twin_ID) - if (nt > 0_pInt) then - - j = 0_pInt - do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families - index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family - do i = 1_pInt,plastic_dislotwin_Nslip(f,instance) ! process each (active) slip system in family - j = j + 1_pInt - - !* Resolved shear stress on slip system - tau = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph)) - !* Stress ratios - if((abs(tau)-state(instance)%threshold_stress_slip(j,of)) > tol_math_check) then - !* Stress ratios - StressRatio_p = ((abs(tau)-state(instance)%threshold_stress_slip(j,of))/& - (plastic_dislotwin_SolidSolutionStrength(instance)+& - plastic_dislotwin_tau_peierlsPerSlipFamily(f,instance)))& - **plastic_dislotwin_pPerSlipFamily(f,instance) - StressRatio_pminus1 = ((abs(tau)-state(instance)%threshold_stress_slip(j,of))/& - (plastic_dislotwin_SolidSolutionStrength(instance)+& - plastic_dislotwin_tau_peierlsPerSlipFamily(f,instance)))& - **(plastic_dislotwin_pPerSlipFamily(f,instance)-1.0_pReal) - !* Boltzmann ratio - BoltzmannRatio = plastic_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature) - !* Initial shear rates - DotGamma0 = & - state(instance)%rhoEdge(j,of)*plastic_dislotwin_burgersPerSlipSystem(j,instance)* & - plastic_dislotwin_v0PerSlipSystem(j,instance) - - !* Shear rates due to slip - gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**& - plastic_dislotwin_qPerSlipFamily(f,instance))*sign(1.0_pReal,tau) - else - gdot_slip(j) = 0.0_pReal - endif - enddo;enddo + Ndot0_twin=0.0_pReal + end if + else isFCCtwin + Ndot0_twin=prm%Ndot0_twin(j) + endif isFCCtwin + StressRatio_r = (mse%threshold_stress_twin(j,of)/tau) **prm%r(j) + postResults(c+j) = (prm%MaxTwinFraction-sumf_twin)*prm%shear_twin(j) & + * mse%twinVolume(j,of)*Ndot0_twin*exp(-StressRatio_r) + endif + enddo + c = c + prm%totalNtwin + case (accumulated_shear_twin_ID) + postResults(c+1_pInt:c+prm%totalNtwin) = stt%accshear_twin(1_pInt:prm%totalNtwin,of) + c = c + prm%totalNtwin + case (mfp_twin_ID) + postResults(c+1_pInt:c+prm%totalNtwin) = mse%mfp_twin(1_pInt:prm%totalNtwin,of) + c = c + prm%totalNtwin + case (resolved_stress_twin_ID) + do j = 1_pInt, prm%totalNtwin + postResults(c+j) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,j)) + enddo + c = c + prm%totalNtwin + case (threshold_stress_twin_ID) + postResults(c+1_pInt:c+prm%totalNtwin) = mse%threshold_stress_twin(1_pInt:prm%totalNtwin,of) + c = c + prm%totalNtwin + case (stress_exponent_ID) + do j = 1_pInt, prm%totalNslip + tau = math_mul33xx33(Mp,prm%Schmid_slip(1:3,1:3,j)) + if((abs(tau)-mse%threshold_stress_slip(j,of)) > tol_math_check) then + StressRatio_p = ((abs(tau)-mse%threshold_stress_slip(j,of))/& + (prm%SolidSolutionStrength+& + prm%tau_peierls(j)))& + **prm%p(j) + StressRatio_pminus1 = ((abs(tau)-mse%threshold_stress_slip(j,of))/& + (prm%SolidSolutionStrength+& + prm%tau_peierls(j)))& + **(prm%p(j)-1.0_pReal) + BoltzmannRatio = prm%Qedge(j)/(kB*Temperature) + DotGamma0 = stt%rhoEdge(j,of)*prm%burgers_slip(j)* prm%v0(j) - j = 0_pInt - do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families - index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family - do i = 1,plastic_dislotwin_Ntwin(f,instance) ! process each (active) twin system in family - j = j + 1_pInt + gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**& + prm%q(j))*sign(1.0_pReal,tau) - tau = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,ph)) - - - !* Shear rates due to twin - if ( tau > 0.0_pReal ) then - select case(lattice_structure(ph)) - case (LATTICE_fcc_ID) - s1=lattice_fcc_twinNucleationSlipPair(1,index_myFamily+i) - s2=lattice_fcc_twinNucleationSlipPair(2,index_myFamily+i) - if (tau < plastic_dislotwin_tau_r_twin(j,instance)) then - Ndot0_twin=(abs(gdot_slip(s1))*(state(instance)%rhoEdge(s2,of)+state(instance)%rhoEdgeDip(s2,of))+& - abs(gdot_slip(s2))*(state(instance)%rhoEdge(s1,of)+state(instance)%rhoEdgeDip(s1,of)))/& - (plastic_dislotwin_L0_twin(instance)*& - plastic_dislotwin_burgersPerSlipSystem(j,instance))*& - (1.0_pReal-exp(-plastic_dislotwin_VcrossSlip(instance)/(kB*Temperature)*& - (plastic_dislotwin_tau_r_twin(j,instance)-tau))) - else - Ndot0_twin=0.0_pReal - end if - case default - Ndot0_twin=plastic_dislotwin_Ndot0PerTwinSystem(j,instance) - end select - StressRatio_r = (state(instance)%threshold_stress_twin(j,of)/tau) & - **plastic_dislotwin_rPerTwinFamily(f,instance) - plastic_dislotwin_postResults(c+j) = & - (plastic_dislotwin_MaxTwinFraction(instance)-sumf)*lattice_shearTwin(index_myFamily+i,ph)*& - state(instance)%twinVolume(j,of)*Ndot0_twin*exp(-StressRatio_r) - endif - - enddo ; enddo - endif - c = c + nt - case (accumulated_shear_twin_ID) - plastic_dislotwin_postResults(c+1_pInt:c+nt) = state(instance)%accshear_twin(1_pInt:nt,of) - c = c + nt - case (mfp_twin_ID) - plastic_dislotwin_postResults(c+1_pInt:c+nt) = state(instance)%mfp_twin(1_pInt:nt,of) - c = c + nt - case (resolved_stress_twin_ID) - if (nt > 0_pInt) then - j = 0_pInt - do f = 1_pInt,lattice_maxNtwinFamily ! loop over all slip families - index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family - do i = 1_pInt,plastic_dislotwin_Ntwin(f,instance) ! process each (active) slip system in family - j = j + 1_pInt - plastic_dislotwin_postResults(c+j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,ph)) - enddo; enddo - endif - c = c + nt - case (threshold_stress_twin_ID) - plastic_dislotwin_postResults(c+1_pInt:c+nt) = state(instance)%threshold_stress_twin(1_pInt:nt,of) - c = c + nt - case (stress_exponent_ID) - j = 0_pInt - do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families - index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family - do i = 1_pInt,plastic_dislotwin_Nslip(f,instance) ! process each (active) slip system in family - j = j + 1_pInt - - !* Resolved shear stress on slip system - tau = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph)) - if((abs(tau)-state(instance)%threshold_stress_slip(j,of)) > tol_math_check) then - !* Stress ratios - StressRatio_p = ((abs(tau)-state(instance)%threshold_stress_slip(j,of))/& - (plastic_dislotwin_SolidSolutionStrength(instance)+& - plastic_dislotwin_tau_peierlsPerSlipFamily(f,instance)))& - **plastic_dislotwin_pPerSlipFamily(f,instance) - StressRatio_pminus1 = ((abs(tau)-state(instance)%threshold_stress_slip(j,of))/& - (plastic_dislotwin_SolidSolutionStrength(instance)+& - plastic_dislotwin_tau_peierlsPerSlipFamily(f,instance)))& - **(plastic_dislotwin_pPerSlipFamily(f,instance)-1.0_pReal) - !* Boltzmann ratio - BoltzmannRatio = plastic_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature) - !* Initial shear rates - DotGamma0 = & - state(instance)%rhoEdge(j,of)*plastic_dislotwin_burgersPerSlipSystem(j,instance)* & - plastic_dislotwin_v0PerSlipSystem(j,instance) - - !* Shear rates due to slip - gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**& - plastic_dislotwin_qPerSlipFamily(f,instance))*sign(1.0_pReal,tau) - - !* Derivatives of shear rates - dgdot_dtauslip = & - abs(gdot_slip(j))*BoltzmannRatio*plastic_dislotwin_pPerSlipFamily(f,instance)& - *plastic_dislotwin_qPerSlipFamily(f,instance)/& - (plastic_dislotwin_SolidSolutionStrength(instance)+& - plastic_dislotwin_tau_peierlsPerSlipFamily(f,instance))*& - StressRatio_pminus1*(1-StressRatio_p)**(plastic_dislotwin_qPerSlipFamily(f,instance)-1.0_pReal) - - else - gdot_slip(j) = 0.0_pReal - dgdot_dtauslip = 0.0_pReal - endif - - !* Stress exponent - plastic_dislotwin_postResults(c+j) = & - merge(0.0_pReal,(tau/gdot_slip(j))*dgdot_dtauslip,dEq0(gdot_slip(j))) - enddo ; enddo - c = c + ns - case (sb_eigenvalues_ID) - plastic_dislotwin_postResults(c+1_pInt:c+3_pInt) = math_eigenvaluesSym33(math_Mandel6to33(Tstar_v)) - c = c + 3_pInt - case (sb_eigenvectors_ID) - call math_eigenValuesVectorsSym33(math_Mandel6to33(Tstar_v),eigValues,eigVectors) - plastic_dislotwin_postResults(c+1_pInt:c+9_pInt) = reshape(eigVectors,[9]) - c = c + 9_pInt - case (stress_trans_fraction_ID) - plastic_dislotwin_postResults(c+1_pInt:c+nr) = & - state(instance)%stressTransFraction(1_pInt:nr,of) - c = c + nr - case (strain_trans_fraction_ID) - plastic_dislotwin_postResults(c+1_pInt:c+nr) = & - state(instance)%strainTransFraction(1_pInt:nr,of) - c = c + nr - case (trans_fraction_ID) - plastic_dislotwin_postResults(c+1_pInt:c+nr) = & - state(instance)%stressTransFraction(1_pInt:nr,of) + & - state(instance)%strainTransFraction(1_pInt:nr,of) - c = c + nr - end select + dgdot_dtauslip = abs(gdot_slip(j))*BoltzmannRatio*prm%p(j) *prm%q(j)/& + (prm%SolidSolutionStrength+ prm%tau_peierls(j))*& + StressRatio_pminus1*(1-StressRatio_p)**(prm%q(j)-1.0_pReal) + else + gdot_slip(j) = 0.0_pReal + dgdot_dtauslip = 0.0_pReal + endif + postResults(c+j) = merge(0.0_pReal,(tau/gdot_slip(j))*dgdot_dtauslip,dEq0(gdot_slip(j))) + enddo + c = c + prm%totalNslip + case (stress_trans_fraction_ID) + postResults(c+1_pInt:c+prm%totalNtrans) = stt%stressTransFraction(1_pInt:prm%totalNtrans,of) + c = c + prm%totalNtrans + case (strain_trans_fraction_ID) + postResults(c+1_pInt:c+prm%totalNtrans) = stt%strainTransFraction(1_pInt:prm%totalNtrans,of) + c = c + prm%totalNtrans + end select enddo + end associate end function plastic_dislotwin_postResults end module plastic_dislotwin diff --git a/src/spectral_damage.f90 b/src/spectral_damage.f90 index c6f2f3822..13ddf0e74 100644 --- a/src/spectral_damage.f90 +++ b/src/spectral_damage.f90 @@ -120,8 +120,8 @@ subroutine spectral_damage_init() trim(snes_type) == 'vinewtonssls') then call DMGetGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr) call DMGetGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr) - call VecSet(lBound,0.0,ierr); CHKERRQ(ierr) - call VecSet(uBound,1.0,ierr); CHKERRQ(ierr) + call VecSet(lBound,0.0_pReal,ierr); CHKERRQ(ierr) + call VecSet(uBound,1.0_pReal,ierr); CHKERRQ(ierr) call SNESVISetVariableBounds(damage_snes,lBound,uBound,ierr) !< variable bounds for variational inequalities like contact mechanics, damage etc. call DMRestoreGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr) call DMRestoreGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr) @@ -134,7 +134,7 @@ subroutine spectral_damage_init() xend = xstart + xend - 1 yend = ystart + yend - 1 zend = zstart + zend - 1 - call VecSet(solution,1.0,ierr); CHKERRQ(ierr) + call VecSet(solution,1.0_pReal,ierr); CHKERRQ(ierr) allocate(damage_current(grid(1),grid(2),grid3), source=1.0_pReal) allocate(damage_lastInc(grid(1),grid(2),grid3), source=1.0_pReal) allocate(damage_stagInc(grid(1),grid(2),grid3), source=1.0_pReal)