diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 8680b2873..165cdc68f 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -81,7 +81,7 @@ checkout: - release ################################################################################################### -processing: +pytest: stage: python script: - cd $DAMASKROOT/python @@ -91,6 +91,15 @@ processing: - master - release +mypy: + stage: python + script: + - cd $DAMASKROOT/python + - mypy -m damask + except: + - master + - release + ################################################################################################### compile_grid_Intel: diff --git a/PRIVATE b/PRIVATE index fabe69749..5a769ec75 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit fabe69749425e8a7aceb3b7c2758b40d97d8b809 +Subproject commit 5a769ec759d9dacc1866c35c6663cd0001e198c5 diff --git a/examples/config/Phase_Dislotungsten_W.yaml b/examples/config/Phase_Dislotungsten_W.yaml deleted file mode 100644 index bf8796cfa..000000000 --- a/examples/config/Phase_Dislotungsten_W.yaml +++ /dev/null @@ -1,26 +0,0 @@ -type: dislotungsten - -N_sl: [12] - -rho_mob_0: [1.0e+9] -rho_dip_0: [1.0] - -nu_a: [9.1e+11] -b_sl: [2.72e-10] -Delta_H_kp,0: [2.61154e-19] # 1.63 eV, Delta_H0 - -tau_Peierls: [2.03e+9] -p_sl: [0.86] -q_sl: [1.69] -h: [2.566e-10] -w: [2.992e-09] -B: [8.3e-5] -D_a: 1.0 # d_edge - -# climb (disabled) -D_0: 0.0 -Q_cl: 0.0 -V_cl: [0.0] - -h_sl-sl: [0.009, 0.72, 0.009, 0.05, 0.05, 0.06, 0.09] -a_nonSchmid: [0.938, 0.71, 4.43] diff --git a/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml b/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml new file mode 100644 index 000000000..55814c3f8 --- /dev/null +++ b/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml @@ -0,0 +1,35 @@ +type: dislotungsten +references: + - D. Cereceda et al., + International Journal of Plasticity 78:242-265, 2016, + http://dx.doi.org/10.1016/j.ijplas.2015.09.002 + - R. Gröger et al., + Acta Materialia 56(19):5412-5425, 2008, + https://doi.org/10.1016/j.actamat.2008.07.037 +output: [Lambda_sl] +N_sl: [12] +b_sl: [2.72e-10] +rho_mob_0: [1.0e+9] # estimated from section 3.2 +rho_dip_0: [1.0] # not given +Q_s: [2.61154e-19] # 1.63 eV, Delta_H0 +B: [8.3e-5] +omega: [9.1e+11] # nu_0 +p_sl: [0.86] +q_sl: [1.69] +tau_Peierls: [2.03e+9] +h: [2.566e-10] +h_sl-sl: [0.009, 0.72, 0.009, 0.05, 0.05, 0.06, 0.09] +w: [2.992e-09] # 11b + +# values in Cereceda et al. are high, using parameters from Gröger et al. +a_nonSchmid: [0.0, 0.56, 0.75] # Table 2 + +# (almost) no annhilation, adjustment needed for simulations beyond the yield point +i_sl: [1] # c, eq. (25) +D: 1.0e+20 # d_g, eq. (25) +D_a: 1.0 # d_edge = D_a*b + +# disable climb (not discussed in Cereceda et al.) +D_0: 0.0 +f_at: 1 +Q_cl: 1.0 diff --git a/python/damask/VERSION b/python/damask/VERSION index 98903601b..bd451b119 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-31-gddb25ad0e +v3.0.0-alpha5-64-g8e08af31e diff --git a/python/damask/__init__.py b/python/damask/__init__.py index 001e46276..231fa8b30 100644 --- a/python/damask/__init__.py +++ b/python/damask/__init__.py @@ -14,8 +14,8 @@ from . import tensor # noqa from . import mechanics # noqa from . import solver # noqa from . import grid_filters # noqa -#Modules that contain only one class (of the same name), are prefixed by a '_'. -#For example, '_colormap' containsa class called 'Colormap' which is imported as 'damask.Colormap'. +# Modules that contain only one class (of the same name), are prefixed by a '_'. +# For example, '_colormap' contains a class called 'Colormap' which is imported as 'damask.Colormap'. from ._rotation import Rotation # noqa from ._crystal import Crystal # noqa from ._orientation import Orientation # noqa diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 3d4d259ff..bb14fd38b 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -125,7 +125,7 @@ class Orientation(Rotation,Crystal): """Create deep copy.""" dup = copy.deepcopy(self) if rotation is not None: - dup.quaternion = Orientation(rotation,family='cubic').quaternion + dup.quaternion = Rotation(rotation).quaternion return dup copy = __copy__ diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 74a3f7419..ac921d70a 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -1,3 +1,5 @@ +import copy + import numpy as np from . import tensor @@ -85,9 +87,12 @@ class Rotation: + str(self.quaternion) - def __copy__(self,**kwargs): + def __copy__(self,rotation=None): """Create deep copy.""" - return self.__class__(rotation=kwargs['rotation'] if 'rotation' in kwargs else self.quaternion) + dup = copy.deepcopy(self) + if rotation is not None: + dup.quaternion = Rotation(rotation).quaternion + return dup copy = __copy__ diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index 816c727cd..42b5a16c4 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -11,10 +11,14 @@ the following operations are required for tensorial data: - D1 = D3.reshape(cells+(-1,)).reshape(-1,9,order='F') """ + +from typing import Sequence, Tuple, Union + from scipy import spatial as _spatial import numpy as _np -def _ks(size,cells,first_order=False): + +def _ks(size: _np.ndarray, cells: Union[_np.ndarray,Sequence[int]], first_order: bool = False) -> _np.ndarray: """ Get wave numbers operator. @@ -41,7 +45,7 @@ def _ks(size,cells,first_order=False): return _np.stack(_np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij'), axis=-1) -def curl(size,f): +def curl(size: _np.ndarray, f: _np.ndarray) -> _np.ndarray: u""" Calculate curl of a vector or tensor field in Fourier space. @@ -72,7 +76,7 @@ def curl(size,f): return _np.fft.irfftn(curl_,axes=(0,1,2),s=f.shape[:3]) -def divergence(size,f): +def divergence(size: _np.ndarray, f: _np.ndarray) -> _np.ndarray: u""" Calculate divergence of a vector or tensor field in Fourier space. @@ -99,7 +103,7 @@ def divergence(size,f): return _np.fft.irfftn(div_,axes=(0,1,2),s=f.shape[:3]) -def gradient(size,f): +def gradient(size: _np.ndarray, f: _np.ndarray) -> _np.ndarray: u""" Calculate gradient of a scalar or vector fieldin Fourier space. @@ -126,7 +130,9 @@ def gradient(size,f): return _np.fft.irfftn(grad_,axes=(0,1,2),s=f.shape[:3]) -def coordinates0_point(cells,size,origin=_np.zeros(3)): +def coordinates0_point(cells: Union[ _np.ndarray,Sequence[int]], + size: _np.ndarray, + origin: _np.ndarray = _np.zeros(3)) -> _np.ndarray: """ Cell center positions (undeformed). @@ -145,8 +151,8 @@ def coordinates0_point(cells,size,origin=_np.zeros(3)): Undeformed cell center coordinates. """ - start = origin + size/cells*.5 - end = origin + size - size/cells*.5 + start = origin + size/_np.array(cells)*.5 + end = origin + size - size/_np.array(cells)*.5 return _np.stack(_np.meshgrid(_np.linspace(start[0],end[0],cells[0]), _np.linspace(start[1],end[1],cells[1]), @@ -154,7 +160,7 @@ def coordinates0_point(cells,size,origin=_np.zeros(3)): axis = -1) -def displacement_fluct_point(size,F): +def displacement_fluct_point(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray: """ Cell center displacement field from fluctuation part of the deformation gradient field. @@ -186,7 +192,7 @@ def displacement_fluct_point(size,F): return _np.fft.irfftn(displacement,axes=(0,1,2),s=F.shape[:3]) -def displacement_avg_point(size,F): +def displacement_avg_point(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray: """ Cell center displacement field from average part of the deformation gradient field. @@ -207,7 +213,7 @@ def displacement_avg_point(size,F): return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),coordinates0_point(F.shape[:3],size)) -def displacement_point(size,F): +def displacement_point(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray: """ Cell center displacement field from deformation gradient field. @@ -227,7 +233,7 @@ def displacement_point(size,F): return displacement_avg_point(size,F) + displacement_fluct_point(size,F) -def coordinates_point(size,F,origin=_np.zeros(3)): +def coordinates_point(size: _np.ndarray, F: _np.ndarray, origin: _np.ndarray = _np.zeros(3)) -> _np.ndarray: """ Cell center positions. @@ -249,7 +255,8 @@ def coordinates_point(size,F,origin=_np.zeros(3)): return coordinates0_point(F.shape[:3],size,origin) + displacement_point(size,F) -def cellsSizeOrigin_coordinates0_point(coordinates0,ordered=True): +def cellsSizeOrigin_coordinates0_point(coordinates0: _np.ndarray, + ordered: bool = True) -> Tuple[_np.ndarray,_np.ndarray,_np.ndarray]: """ Return grid 'DNA', i.e. cells, size, and origin from 1D array of point positions. @@ -292,13 +299,15 @@ def cellsSizeOrigin_coordinates0_point(coordinates0,ordered=True): raise ValueError('Regular cell spacing violated.') if ordered and not _np.allclose(coordinates0.reshape(tuple(cells)+(3,),order='F'), - coordinates0_point(cells,size,origin),atol=atol): + coordinates0_point(list(cells),size,origin),atol=atol): raise ValueError('Input data is not ordered (x fast, z slow).') return (cells,size,origin) -def coordinates0_node(cells,size,origin=_np.zeros(3)): +def coordinates0_node(cells: Union[_np.ndarray,Sequence[int]], + size: _np.ndarray, + origin: _np.ndarray = _np.zeros(3)) -> _np.ndarray: """ Nodal positions (undeformed). @@ -323,7 +332,7 @@ def coordinates0_node(cells,size,origin=_np.zeros(3)): axis = -1) -def displacement_fluct_node(size,F): +def displacement_fluct_node(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray: """ Nodal displacement field from fluctuation part of the deformation gradient field. @@ -343,7 +352,7 @@ def displacement_fluct_node(size,F): return point_to_node(displacement_fluct_point(size,F)) -def displacement_avg_node(size,F): +def displacement_avg_node(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray: """ Nodal displacement field from average part of the deformation gradient field. @@ -364,7 +373,7 @@ def displacement_avg_node(size,F): return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),coordinates0_node(F.shape[:3],size)) -def displacement_node(size,F): +def displacement_node(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray: """ Nodal displacement field from deformation gradient field. @@ -384,7 +393,7 @@ def displacement_node(size,F): return displacement_avg_node(size,F) + displacement_fluct_node(size,F) -def coordinates_node(size,F,origin=_np.zeros(3)): +def coordinates_node(size: _np.ndarray, F: _np.ndarray, origin: _np.ndarray = _np.zeros(3)) -> _np.ndarray: """ Nodal positions. @@ -406,7 +415,8 @@ def coordinates_node(size,F,origin=_np.zeros(3)): return coordinates0_node(F.shape[:3],size,origin) + displacement_node(size,F) -def cellsSizeOrigin_coordinates0_node(coordinates0,ordered=True): +def cellsSizeOrigin_coordinates0_node(coordinates0: _np.ndarray, + ordered: bool = True) -> Tuple[_np.ndarray,_np.ndarray,_np.ndarray]: """ Return grid 'DNA', i.e. cells, size, and origin from 1D array of nodal positions. @@ -441,13 +451,13 @@ def cellsSizeOrigin_coordinates0_node(coordinates0,ordered=True): raise ValueError('Regular cell spacing violated.') if ordered and not _np.allclose(coordinates0.reshape(tuple(cells+1)+(3,),order='F'), - coordinates0_node(cells,size,origin),atol=atol): + coordinates0_node(list(cells),size,origin),atol=atol): raise ValueError('Input data is not ordered (x fast, z slow).') return (cells,size,origin) -def point_to_node(cell_data): +def point_to_node(cell_data: _np.ndarray) -> _np.ndarray: """ Interpolate periodic point data to nodal data. @@ -469,7 +479,7 @@ def point_to_node(cell_data): return _np.pad(n,((0,1),(0,1),(0,1))+((0,0),)*len(cell_data.shape[3:]),mode='wrap') -def node_to_point(node_data): +def node_to_point(node_data: _np.ndarray) -> _np.ndarray: """ Interpolate periodic nodal data to point data. @@ -491,7 +501,7 @@ def node_to_point(node_data): return c[1:,1:,1:] -def coordinates0_valid(coordinates0): +def coordinates0_valid(coordinates0: _np.ndarray) -> bool: """ Check whether coordinates form a regular grid. @@ -513,7 +523,7 @@ def coordinates0_valid(coordinates0): return False -def regrid(size,F,cells): +def regrid(size: _np.ndarray, F: _np.ndarray, cells: Union[_np.ndarray,Sequence[int]]) -> _np.ndarray: """ Return mapping from coordinates in deformed configuration to a regular grid. diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index 0e160523b..1a03f390b 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -5,13 +5,15 @@ All routines operate on numpy.ndarrays of shape (...,3,3). """ -from . import tensor as _tensor -from . import _rotation +from typing import Sequence import numpy as _np +from . import tensor as _tensor +from . import _rotation -def deformation_Cauchy_Green_left(F): + +def deformation_Cauchy_Green_left(F: _np.ndarray) -> _np.ndarray: """ Calculate left Cauchy-Green deformation tensor (Finger deformation tensor). @@ -29,7 +31,7 @@ def deformation_Cauchy_Green_left(F): return _np.matmul(F,_tensor.transpose(F)) -def deformation_Cauchy_Green_right(F): +def deformation_Cauchy_Green_right(F: _np.ndarray) -> _np.ndarray: """ Calculate right Cauchy-Green deformation tensor. @@ -47,7 +49,7 @@ def deformation_Cauchy_Green_right(F): return _np.matmul(_tensor.transpose(F),F) -def equivalent_strain_Mises(epsilon): +def equivalent_strain_Mises(epsilon: _np.ndarray) -> _np.ndarray: """ Calculate the Mises equivalent of a strain tensor. @@ -65,7 +67,7 @@ def equivalent_strain_Mises(epsilon): return _equivalent_Mises(epsilon,2.0/3.0) -def equivalent_stress_Mises(sigma): +def equivalent_stress_Mises(sigma: _np.ndarray) -> _np.ndarray: """ Calculate the Mises equivalent of a stress tensor. @@ -83,7 +85,7 @@ def equivalent_stress_Mises(sigma): return _equivalent_Mises(sigma,3.0/2.0) -def maximum_shear(T_sym): +def maximum_shear(T_sym: _np.ndarray) -> _np.ndarray: """ Calculate the maximum shear component of a symmetric tensor. @@ -102,7 +104,7 @@ def maximum_shear(T_sym): return (w[...,0] - w[...,2])*0.5 -def rotation(T): +def rotation(T: _np.ndarray) -> _rotation.Rotation: """ Calculate the rotational part of a tensor. @@ -120,7 +122,7 @@ def rotation(T): return _rotation.Rotation.from_matrix(_polar_decomposition(T,'R')[0]) -def strain(F,t,m): +def strain(F: _np.ndarray, t: str, m: float) -> _np.ndarray: """ Calculate strain tensor (Seth–Hill family). @@ -160,7 +162,7 @@ def strain(F,t,m): return eps -def stress_Cauchy(P,F): +def stress_Cauchy(P: _np.ndarray, F: _np.ndarray) -> _np.ndarray: """ Calculate the Cauchy stress (true stress). @@ -182,7 +184,7 @@ def stress_Cauchy(P,F): return _tensor.symmetric(_np.einsum('...,...ij,...kj',1.0/_np.linalg.det(F),P,F)) -def stress_second_Piola_Kirchhoff(P,F): +def stress_second_Piola_Kirchhoff(P: _np.ndarray, F: _np.ndarray) -> _np.ndarray: """ Calculate the second Piola-Kirchhoff stress. @@ -205,7 +207,7 @@ def stress_second_Piola_Kirchhoff(P,F): return _tensor.symmetric(_np.einsum('...ij,...jk',_np.linalg.inv(F),P)) -def stretch_left(T): +def stretch_left(T: _np.ndarray) -> _np.ndarray: """ Calculate left stretch of a tensor. @@ -223,7 +225,7 @@ def stretch_left(T): return _polar_decomposition(T,'V')[0] -def stretch_right(T): +def stretch_right(T: _np.ndarray) -> _np.ndarray: """ Calculate right stretch of a tensor. @@ -241,7 +243,7 @@ def stretch_right(T): return _polar_decomposition(T,'U')[0] -def _polar_decomposition(T,requested): +def _polar_decomposition(T: _np.ndarray, requested: Sequence[str]) -> tuple: """ Perform singular value decomposition. @@ -257,21 +259,21 @@ def _polar_decomposition(T,requested): u, _, vh = _np.linalg.svd(T) R = _np.einsum('...ij,...jk',u,vh) - output = () + output = [] if 'R' in requested: - output+=(R,) + output+=[R] if 'V' in requested: - output+=(_np.einsum('...ij,...kj',T,R),) + output+=[_np.einsum('...ij,...kj',T,R)] if 'U' in requested: - output+=(_np.einsum('...ji,...jk',R,T),) + output+=[_np.einsum('...ji,...jk',R,T)] if len(output) == 0: raise ValueError('output needs to be out of V, R, U') - return output + return tuple(output) -def _equivalent_Mises(T_sym,s): +def _equivalent_Mises(T_sym: _np.ndarray, s: float) -> _np.ndarray: """ Base equation for Mises equivalent of a stress or strain tensor. diff --git a/python/damask/seeds.py b/python/damask/seeds.py index 26aa3084b..4d5a8c624 100644 --- a/python/damask/seeds.py +++ b/python/damask/seeds.py @@ -1,5 +1,7 @@ """Functionality for generation of seed points for Voronoi or Laguerre tessellation.""" +from typing import Sequence,Tuple + from scipy import spatial as _spatial import numpy as _np @@ -7,7 +9,7 @@ from . import util as _util from . import grid_filters as _grid_filters -def from_random(size,N_seeds,cells=None,rng_seed=None): +def from_random(size: _np.ndarray, N_seeds: int, cells: _np.ndarray = None, rng_seed=None) -> _np.ndarray: """ Place seeds randomly in space. @@ -41,7 +43,8 @@ def from_random(size,N_seeds,cells=None,rng_seed=None): return coords -def from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=True,rng_seed=None): +def from_Poisson_disc(size: _np.ndarray, N_seeds: int, N_candidates: int, distance: float, + periodic: bool = True, rng_seed=None) -> _np.ndarray: """ Place seeds according to a Poisson disc distribution. @@ -75,18 +78,17 @@ def from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=True,rng_seed= i = 0 progress = _util._ProgressBar(N_seeds+1,'',50) while s < N_seeds: + i += 1 candidates = rng.random((N_candidates,3))*_np.broadcast_to(size,(N_candidates,3)) tree = _spatial.cKDTree(coords[:s],boxsize=size) if periodic else \ _spatial.cKDTree(coords[:s]) distances = tree.query(candidates)[0] best = distances.argmax() if distances[best] > distance: # require minimum separation + i = 0 coords[s] = candidates[best] # maximum separation to existing point cloud s += 1 progress.update(s) - i = 0 - else: - i += 1 if i == 100: raise ValueError('Seeding not possible') @@ -94,22 +96,23 @@ def from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=True,rng_seed= return coords -def from_grid(grid,selection=None,invert=False,average=False,periodic=True): +def from_grid(grid, selection: Sequence[int] = None, + invert: bool = False, average: bool = False, periodic: bool = True) -> Tuple[_np.ndarray, _np.ndarray]: """ Create seeds from grid description. Parameters ---------- grid : damask.Grid - Grid, from which the material IDs are used as seeds. + Grid from which the material IDs are used as seeds. selection : iterable of integers, optional Material IDs to consider. invert : boolean, false - Do not consider the material IDs given in selection. Defaults to False. + Consider all material IDs except those in selection. Defaults to False. average : boolean, optional Seed corresponds to center of gravity of material ID cloud. periodic : boolean, optional - Center of gravity with periodic boundaries. + Center of gravity accounts for periodic boundaries. Returns ------- diff --git a/python/damask/tensor.py b/python/damask/tensor.py index cf5d94020..a735b355e 100644 --- a/python/damask/tensor.py +++ b/python/damask/tensor.py @@ -8,7 +8,7 @@ All routines operate on numpy.ndarrays of shape (...,3,3). import numpy as _np -def deviatoric(T): +def deviatoric(T: _np.ndarray) -> _np.ndarray: """ Calculate deviatoric part of a tensor. @@ -26,7 +26,7 @@ def deviatoric(T): return T - spherical(T,tensor=True) -def eigenvalues(T_sym): +def eigenvalues(T_sym: _np.ndarray) -> _np.ndarray: """ Eigenvalues, i.e. principal components, of a symmetric tensor. @@ -45,7 +45,7 @@ def eigenvalues(T_sym): return _np.linalg.eigvalsh(symmetric(T_sym)) -def eigenvectors(T_sym,RHS=False): +def eigenvectors(T_sym: _np.ndarray, RHS: bool = False) -> _np.ndarray: """ Eigenvectors of a symmetric tensor. @@ -70,7 +70,7 @@ def eigenvectors(T_sym,RHS=False): return v -def spherical(T,tensor=True): +def spherical(T: _np.ndarray, tensor: bool = True) -> _np.ndarray: """ Calculate spherical part of a tensor. @@ -92,7 +92,7 @@ def spherical(T,tensor=True): return _np.einsum('...jk,...',_np.eye(3),sph) if tensor else sph -def symmetric(T): +def symmetric(T: _np.ndarray) -> _np.ndarray: """ Symmetrize tensor. @@ -110,7 +110,7 @@ def symmetric(T): return (T+transpose(T))*0.5 -def transpose(T): +def transpose(T: _np.ndarray) -> _np.ndarray: """ Transpose tensor. diff --git a/python/mypy.ini b/python/mypy.ini new file mode 100644 index 000000000..01001daa6 --- /dev/null +++ b/python/mypy.ini @@ -0,0 +1,14 @@ +[mypy-scipy.*] +ignore_missing_imports = True +[mypy-h5py.*] +ignore_missing_imports = True +[mypy-vtk.*] +ignore_missing_imports = True +[mypy-PIL.*] +ignore_missing_imports = True +[mypy-matplotlib.*] +ignore_missing_imports = True +[mypy-pandas.*] +ignore_missing_imports = True +[mypy-wx.*] +ignore_missing_imports = True diff --git a/src/IO.f90 b/src/IO.f90 index cd7c09c75..717493006 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -432,7 +432,7 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) msg = 'Nconstituents mismatch between homogenization and material' !-------------------------------------------------------------------------------------------------- -! material error messages and related messages in mesh +! material error messages and related messages in geometry case (150) msg = 'index out of bounds' case (153) @@ -499,6 +499,11 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) case (710) msg = 'Closing quotation mark missing in string' +!------------------------------------------------------------------------------------------------- +! errors related to the mesh solver + case (821) + msg = 'order not supported' + !------------------------------------------------------------------------------------------------- ! errors related to the grid solver case (831) diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 64bcc3896..a3856ccaa 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -19,6 +19,7 @@ module FEM_utilities use IO use discretization_mesh use homogenization + use FEM_quadrature implicit none private @@ -29,8 +30,8 @@ module FEM_utilities !-------------------------------------------------------------------------------------------------- ! field labels information - character(len=*), parameter, public :: & - FIELD_MECH_label = 'mechanical' + character(len=*), parameter, public :: & + FIELD_MECH_label = 'mechanical' enum, bind(c); enumerator :: & FIELD_UNDEFINED_ID, & @@ -86,7 +87,9 @@ subroutine FEM_utilities_init class(tNode), pointer :: & num_mesh, & debug_mesh ! pointer to mesh debug options - integer :: structOrder !< order of displacement shape functions + integer :: & + p_s, & !< order of shape functions + p_i !< integration order (quadrature rule) character(len=*), parameter :: & PETSCDEBUG = ' -snes_view -snes_monitor ' PetscErrorCode :: ierr @@ -96,7 +99,14 @@ subroutine FEM_utilities_init print'(/,a)', ' <<<+- FEM_utilities init -+>>>' num_mesh => config_numerics%get('mesh',defaultVal=emptyDict) - structOrder = num_mesh%get_asInt('structOrder', defaultVal = 2) + + p_s = num_mesh%get_asInt('p_s',defaultVal = 2) + p_i = num_mesh%get_asInt('p_i',defaultVal = p_s) + + if (p_s < 1_pInt .or. p_s > size(FEM_nQuadrature,2)) & + call IO_error(821,ext_msg='shape function order (p_s) out of bounds') + if (p_i < max(1_pInt,p_s-1_pInt) .or. p_i > p_s) & + call IO_error(821,ext_msg='integration order (p_i) out of bounds') debug_mesh => config_debug%get('mesh',defaultVal=emptyList) debugPETSc = debug_mesh%contains('PETSc') @@ -119,7 +129,7 @@ subroutine FEM_utilities_init CHKERRQ(ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_mesh%get_asString('PETSc_options',defaultVal=''),ierr) CHKERRQ(ierr) - write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', structOrder + write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', p_s call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),ierr) CHKERRQ(ierr) diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index e5f989484..f1d38760d 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -85,7 +85,7 @@ subroutine discretization_mesh_init(restart) materialAt class(tNode), pointer :: & num_mesh - integer :: integrationOrder !< order of quadrature rule required + integer :: p_i !< integration order (quadrature rule) type(tvec) :: coords_node0 print'(/,a)', ' <<<+- discretization_mesh init -+>>>' @@ -93,7 +93,7 @@ subroutine discretization_mesh_init(restart) !-------------------------------------------------------------------------------- ! read numerics parameter num_mesh => config_numerics%get('mesh',defaultVal=emptyDict) - integrationOrder = num_mesh%get_asInt('integrationorder',defaultVal = 2) + p_i = num_mesh%get_asInt('p_i',defaultVal = 2) !--------------------------------------------------------------------------------- ! read debug parameters @@ -150,9 +150,9 @@ subroutine discretization_mesh_init(restart) call VecGetArrayF90(coords_node0, mesh_node0_temp,ierr) CHKERRQ(ierr) - mesh_maxNips = FEM_nQuadrature(dimPlex,integrationOrder) + mesh_maxNips = FEM_nQuadrature(dimPlex,p_i) - call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,integrationOrder)%p) + call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,p_i)%p) call mesh_FEM_build_ipVolumes(dimPlex) allocate(materialAt(mesh_NcpElems)) diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index d6d314a42..496a82fc5 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -41,7 +41,7 @@ module mesh_mechanical_FEM type, private :: tNumerics integer :: & - integrationOrder, & !< order of quadrature rule required + p_i, & !< integration order (quadrature rule) itmax logical :: & BBarStabilisation @@ -118,7 +118,7 @@ subroutine FEM_mechanical_init(fieldBC) !----------------------------------------------------------------------------- ! read numerical parametes and do sanity checks num_mesh => config_numerics%get('mesh',defaultVal=emptyDict) - num%integrationOrder = num_mesh%get_asInt('integrationorder',defaultVal = 2) + num%p_i = num_mesh%get_asInt('p_i',defaultVal = 2) num%itmax = num_mesh%get_asInt('itmax',defaultVal=250) num%BBarStabilisation = num_mesh%get_asBool('bbarstabilisation',defaultVal = .false.) num%eps_struct_atol = num_mesh%get_asFloat('eps_struct_atol', defaultVal = 1.0e-10_pReal) @@ -135,9 +135,9 @@ subroutine FEM_mechanical_init(fieldBC) !-------------------------------------------------------------------------------------------------- ! Setup FEM mech discretization - qPoints = FEM_quadrature_points( dimPlex,num%integrationOrder)%p - qWeights = FEM_quadrature_weights(dimPlex,num%integrationOrder)%p - nQuadrature = FEM_nQuadrature( dimPlex,num%integrationOrder) + qPoints = FEM_quadrature_points( dimPlex,num%p_i)%p + qWeights = FEM_quadrature_weights(dimPlex,num%p_i)%p + nQuadrature = FEM_nQuadrature( dimPlex,num%p_i) qPointsP => qPoints qWeightsP => qWeights call PetscQuadratureCreate(PETSC_COMM_SELF,mechQuad,ierr); CHKERRQ(ierr) @@ -146,7 +146,7 @@ subroutine FEM_mechanical_init(fieldBC) call PetscQuadratureSetData(mechQuad,dimPlex,nc,nQuadrature,qPointsP,qWeightsP,ierr) CHKERRQ(ierr) call PetscFECreateDefault(PETSC_COMM_SELF,dimPlex,nc,PETSC_TRUE,prefix, & - num%integrationOrder,mechFE,ierr); CHKERRQ(ierr) + num%p_i,mechFE,ierr); CHKERRQ(ierr) call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr) call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr) nBasis = nBasis/nc diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 1e5e120e5..c759cdaad 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -24,7 +24,6 @@ submodule(phase:plastic) dislotungsten tau_Peierls, & !< Peierls stress !* mobility law parameters Q_s, & !< activation energy for glide [J] - v_0, & !< dislocation velocity prefactor [m/s] p, & !< p-exponent in glide velocity q, & !< q-exponent in glide velocity B, & !< friction coefficient @@ -148,7 +147,7 @@ module function plastic_dislotungsten_init() result(myPlasticity) else prm%P_nS_pos = prm%P_sl prm%P_nS_neg = prm%P_sl - endif + end if prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'), & phase_lattice(ph)) @@ -158,7 +157,6 @@ module function plastic_dislotungsten_init() result(myPlasticity) rho_mob_0 = pl%get_as1dFloat('rho_mob_0', requiredSize=size(N_sl)) rho_dip_0 = pl%get_as1dFloat('rho_dip_0', requiredSize=size(N_sl)) - prm%v_0 = pl%get_as1dFloat('v_0', requiredSize=size(N_sl)) prm%b_sl = pl%get_as1dFloat('b_sl', requiredSize=size(N_sl)) prm%Q_s = pl%get_as1dFloat('Q_s', requiredSize=size(N_sl)) @@ -189,18 +187,16 @@ module function plastic_dislotungsten_init() result(myPlasticity) prm%w = math_expand(prm%w, N_sl) prm%omega = math_expand(prm%omega, N_sl) prm%tau_Peierls = math_expand(prm%tau_Peierls, N_sl) - prm%v_0 = math_expand(prm%v_0, N_sl) prm%B = math_expand(prm%B, N_sl) prm%i_sl = math_expand(prm%i_sl, N_sl) prm%f_at = math_expand(prm%f_at, N_sl) prm%d_caron = pl%get_asFloat('D_a') * prm%b_sl ! sanity checks - if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0' + if ( prm%D_0 < 0.0_pReal) extmsg = trim(extmsg)//' D_0' if ( prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' if (any(rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0' if (any(rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0' - if (any(prm%v_0 < 0.0_pReal)) extmsg = trim(extmsg)//' v_0' if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl' if (any(prm%Q_s <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s' if (any(prm%tau_Peierls < 0.0_pReal)) extmsg = trim(extmsg)//' tau_Peierls' @@ -209,13 +205,13 @@ module function plastic_dislotungsten_init() result(myPlasticity) if (any(prm%f_at <= 0.0_pReal)) extmsg = trim(extmsg)//' f_at or b_sl' else slipActive - rho_mob_0= emptyRealArray; rho_dip_0 = emptyRealArray + rho_mob_0 = emptyRealArray; rho_dip_0 = emptyRealArray allocate(prm%b_sl,prm%d_caron,prm%i_sl,prm%f_at,prm%tau_Peierls, & - prm%Q_s,prm%v_0,prm%p,prm%q,prm%B,prm%h,prm%w,prm%omega, & + prm%Q_s,prm%p,prm%q,prm%B,prm%h,prm%w,prm%omega, & source = emptyRealArray) allocate(prm%forestProjection(0,0)) allocate(prm%h_sl_sl (0,0)) - endif slipActive + end if slipActive !-------------------------------------------------------------------------------------------------- ! allocate state arrays @@ -258,7 +254,7 @@ module function plastic_dislotungsten_init() result(myPlasticity) ! exit if any parameter is out of range if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(dislotungsten)') - enddo + end do end function plastic_dislotungsten_init @@ -267,7 +263,7 @@ end function plastic_dislotungsten_init !> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, & - Mp,T,ph,en) + Mp,T,ph,en) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -287,19 +283,20 @@ pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, & dot_gamma_pos,dot_gamma_neg, & ddot_gamma_dtau_pos,ddot_gamma_dtau_neg + Lp = 0.0_pReal dLp_dMp = 0.0_pReal associate(prm => param(ph)) - call kinetics(Mp,T,ph,en,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg) - do i = 1, prm%sum_N_sl - Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P_sl(1:3,1:3,i) - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + ddot_gamma_dtau_pos(i) * prm%P_sl(k,l,i) * prm%P_nS_pos(m,n,i) & - + ddot_gamma_dtau_neg(i) * prm%P_sl(k,l,i) * prm%P_nS_neg(m,n,i) - enddo + call kinetics(Mp,T,ph,en,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg) + do i = 1, prm%sum_N_sl + Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P_sl(1:3,1:3,i) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + ddot_gamma_dtau_pos(i) * prm%P_sl(k,l,i) * prm%P_nS_pos(m,n,i) & + + ddot_gamma_dtau_neg(i) * prm%P_sl(k,l,i) * prm%P_nS_neg(m,n,i) + end do end associate @@ -328,35 +325,36 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en) dot_rho_dip_climb, & d_hat + associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph)) - call kinetics(Mp,T,ph,en,& - dot_gamma_pos,dot_gamma_neg, & - tau_pos_out = tau_pos,tau_neg_out = tau_neg) + call kinetics(Mp,T,ph,en,& + dot_gamma_pos,dot_gamma_neg, & + tau_pos_out = tau_pos,tau_neg_out = tau_neg) - dot%gamma_sl(:,en) = abs(dot_gamma_pos+dot_gamma_neg) + dot%gamma_sl(:,en) = abs(dot_gamma_pos+dot_gamma_neg) - where(dEq0(tau_pos)) ! ToDo: use avg of +/- - dot_rho_dip_formation = 0.0_pReal - dot_rho_dip_climb = 0.0_pReal - else where - d_hat = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), & ! ToDo: use avg of +/- - prm%d_caron, & ! lower limit - dst%Lambda_sl(:,en)) ! upper limit - dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot%gamma_sl(:,en)/prm%b_sl, & - 0.0_pReal, & - prm%dipoleformation) - v_cl = (3.0_pReal*prm%mu*prm%D_0*exp(-prm%Q_cl/(kB*T))*prm%f_at/(2.0_pReal*PI*kB*T)) & - * (1.0_pReal/(d_hat+prm%d_caron)) - dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(d_hat-prm%d_caron) ! ToDo: Discuss with Franz: Stress dependency? - end where + where(dEq0((tau_pos+tau_neg)*0.5_pReal)) + dot_rho_dip_formation = 0.0_pReal + dot_rho_dip_climb = 0.0_pReal + else where + d_hat = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos+tau_neg)*0.5_pReal), & + prm%d_caron, & ! lower limit + dst%Lambda_sl(:,en)) ! upper limit + dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot%gamma_sl(:,en)/prm%b_sl, & + 0.0_pReal, & + prm%dipoleformation) + v_cl = (3.0_pReal*prm%mu*prm%D_0*exp(-prm%Q_cl/(kB*T))*prm%f_at/(2.0_pReal*PI*kB*T)) & + * (1.0_pReal/(d_hat+prm%d_caron)) + dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(d_hat-prm%d_caron) ! ToDo: Discuss with Franz: Stress dependency? + end where - dot%rho_mob(:,en) = dot%gamma_sl(:,en)/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication - - dot_rho_dip_formation & - - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*dot%gamma_sl(:,en) ! Spontaneous annihilation of 2 edges - dot%rho_dip(:,en) = dot_rho_dip_formation & - - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*dot%gamma_sl(:,en) & ! Spontaneous annihilation of an edge with a dipole - - dot_rho_dip_climb + dot%rho_mob(:,en) = dot%gamma_sl(:,en)/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication + - dot_rho_dip_formation & + - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*dot%gamma_sl(:,en) ! Spontaneous annihilation of 2 edges + dot%rho_dip(:,en) = dot_rho_dip_formation & + - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*dot%gamma_sl(:,en) & ! Spontaneous annihilation of an edge with a dipole + - dot_rho_dip_climb end associate @@ -368,21 +366,22 @@ end subroutine dislotungsten_dotState !-------------------------------------------------------------------------------------------------- module subroutine dislotungsten_dependentState(ph,en) - integer, intent(in) :: & + integer, intent(in) :: & ph, & en real(pReal), dimension(param(ph)%sum_N_sl) :: & - dislocationSpacing + Lambda_sl_inv associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) - dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en))) dst%tau_pass(:,en) = prm%mu*prm%b_sl & * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) - dst%Lambda_sl(:,en) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl) + Lambda_sl_inv = 1.0_pReal/prm%D & + + sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))/prm%i_sl + dst%Lambda_sl(:,en) = Lambda_sl_inv**(-1.0_pReal) end associate @@ -423,7 +422,7 @@ module subroutine plastic_dislotungsten_results(ph,group) 'threshold stress for slip','Pa',prm%systems_sl) end select - enddo + end do end associate @@ -456,88 +455,91 @@ pure subroutine kinetics(Mp,T,ph,en, & ddot_gamma_dtau_neg, & tau_pos_out, & tau_neg_out + real(pReal), dimension(param(ph)%sum_N_sl) :: & StressRatio, & StressRatio_p,StressRatio_pminus1, & - dvel, vel, & - tau_pos,tau_neg, & + dvel, & + tau_pos, tau_neg, tau_eff, & t_n, t_k, dtk,dtn integer :: j + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) - do j = 1, prm%sum_N_sl - tau_pos(j) = math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,j)) - tau_neg(j) = math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,j)) - enddo + do j = 1, prm%sum_N_sl + tau_pos(j) = math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,j)) + tau_neg(j) = math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,j)) + end do + if (present(tau_pos_out)) tau_pos_out = tau_pos + if (present(tau_neg_out)) tau_neg_out = tau_neg - if (present(tau_pos_out)) tau_pos_out = tau_pos - if (present(tau_neg_out)) tau_neg_out = tau_neg + associate(BoltzmannRatio => prm%Q_s/(kB*T), & + b_rho_half => stt%rho_mob(:,en) * prm%b_sl * 0.5_pReal, & + effectiveLength => dst%Lambda_sl(:,en) - prm%w) - associate(BoltzmannRatio => prm%Q_s/(kB*T), & - dot_gamma_0 => stt%rho_mob(:,en)*prm%b_sl*prm%v_0, & - effectiveLength => dst%Lambda_sl(:,en) - prm%w) + tau_eff = abs(tau_pos)-dst%tau_pass(:,en) - significantPositiveTau: where(abs(tau_pos)-dst%tau_pass(:,en) > tol_math_check) - StressRatio = (abs(tau_pos)-dst%tau_pass(:,en))/prm%tau_Peierls - StressRatio_p = StressRatio** prm%p - StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) + significantPositiveTau: where(tau_eff > tol_math_check) + StressRatio = tau_eff/prm%tau_Peierls + StressRatio_p = StressRatio** prm%p + StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) - t_n = prm%b_sl/(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)*prm%omega*effectiveLength) - t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) + t_n = prm%b_sl*exp(BoltzmannRatio*(1.0_pReal-StressRatio_p) ** prm%q) & + / (prm%omega*effectiveLength) + t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_eff) ! corrected eq. (14) - vel = prm%h/(t_n + t_k) + dot_gamma_pos = b_rho_half * sign(prm%h/(t_n + t_k),tau_pos) + else where significantPositiveTau + dot_gamma_pos = 0.0_pReal + end where significantPositiveTau - dot_gamma_pos = dot_gamma_0 * sign(vel,tau_pos) * 0.5_pReal - else where significantPositiveTau - dot_gamma_pos = 0.0_pReal - end where significantPositiveTau + if (present(ddot_gamma_dtau_pos)) then + significantPositiveTau2: where(abs(tau_pos)-dst%tau_pass(:,en) > tol_math_check) + dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & + * StressRatio_pminus1 / prm%tau_Peierls + dtk = -1.0_pReal * t_k / tau_pos - if (present(ddot_gamma_dtau_pos)) then - significantPositiveTau2: where(abs(tau_pos)-dst%tau_pass(:,en) > tol_math_check) - dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & - * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls - dtk = -1.0_pReal * t_k / tau_pos + dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal - dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal + ddot_gamma_dtau_pos = b_rho_half * dvel + else where significantPositiveTau2 + ddot_gamma_dtau_pos = 0.0_pReal + end where significantPositiveTau2 + end if - ddot_gamma_dtau_pos = dot_gamma_0 * dvel* 0.5_pReal - else where significantPositiveTau2 - ddot_gamma_dtau_pos = 0.0_pReal - end where significantPositiveTau2 - endif + tau_eff = abs(tau_neg)-dst%tau_pass(:,en) - significantNegativeTau: where(abs(tau_neg)-dst%tau_pass(:,en) > tol_math_check) - StressRatio = (abs(tau_neg)-dst%tau_pass(:,en))/prm%tau_Peierls - StressRatio_p = StressRatio** prm%p - StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) + significantNegativeTau: where(tau_eff > tol_math_check) + StressRatio = tau_eff/prm%tau_Peierls + StressRatio_p = StressRatio** prm%p + StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) - t_n = prm%b_sl/(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)*prm%omega*effectiveLength) - t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) + t_n = prm%b_sl*exp(BoltzmannRatio*(1.0_pReal-StressRatio_p) ** prm%q) & + / (prm%omega*effectiveLength) + t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_eff) ! corrected eq. (14) - vel = prm%h/(t_n + t_k) + dot_gamma_neg = b_rho_half * sign(prm%h/(t_n + t_k),tau_neg) + else where significantNegativeTau + dot_gamma_neg = 0.0_pReal + end where significantNegativeTau - dot_gamma_neg = dot_gamma_0 * sign(vel,tau_neg) * 0.5_pReal - else where significantNegativeTau - dot_gamma_neg = 0.0_pReal - end where significantNegativeTau + if (present(ddot_gamma_dtau_neg)) then + significantNegativeTau2: where(abs(tau_neg)-dst%tau_pass(:,en) > tol_math_check) + dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & + * StressRatio_pminus1 / prm%tau_Peierls + dtk = -1.0_pReal * t_k / tau_neg - if (present(ddot_gamma_dtau_neg)) then - significantNegativeTau2: where(abs(tau_neg)-dst%tau_pass(:,en) > tol_math_check) - dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & - * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls - dtk = -1.0_pReal * t_k / tau_neg + dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal - dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal + ddot_gamma_dtau_neg = b_rho_half * dvel + else where significantNegativeTau2 + ddot_gamma_dtau_neg = 0.0_pReal + end where significantNegativeTau2 + end if - ddot_gamma_dtau_neg = dot_gamma_0 * dvel * 0.5_pReal - else where significantNegativeTau2 - ddot_gamma_dtau_neg = 0.0_pReal - end where significantNegativeTau2 - end if - - end associate + end associate end associate end subroutine kinetics diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index ac179d775..de73cee04 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -890,7 +890,8 @@ pure subroutine kinetics_sl(Mp,T,ph,en, & stressRatio = tau_eff/prm%tau_0 StressRatio_p = stressRatio** prm%p Q_kB_T = prm%Q_sl/(kB*T) - v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(Q_kB_T*(1.0_pReal-StressRatio_p)** prm%q) + v_wait_inverse = exp(Q_kB_T*(1.0_pReal-StressRatio_p)** prm%q) & + / prm%v_0 v_run_inverse = prm%B/(tau_eff*prm%b_sl) dot_gamma_sl = sign(stt%rho_mob(:,en)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau) diff --git a/src/prec.f90 b/src/prec.f90 index 2e67ae76a..d6d161a94 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -101,7 +101,7 @@ logical elemental pure function dEq(a,b,tol) dEq = abs(a-b) <= tol else dEq = abs(a-b) <= PREAL_EPSILON * maxval(abs([a,b])) - endif + end if end function dEq @@ -139,7 +139,7 @@ logical elemental pure function dEq0(a,tol) dEq0 = abs(a) <= tol else dEq0 = abs(a) <= PREAL_MIN * 10.0_pReal - endif + end if end function dEq0 @@ -178,7 +178,7 @@ logical elemental pure function cEq(a,b,tol) cEq = abs(a-b) <= tol else cEq = abs(a-b) <= PREAL_EPSILON * maxval(abs([a,b])) - endif + end if end function cEq