From 596a8ecd119fded5604e02043bfabfa7bea0ccfe Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 6 Dec 2021 18:34:14 +0100 Subject: [PATCH 01/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-251-g65c4417a2 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 656df5ec5..2b61790a8 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-244-gb57351045 +v3.0.0-alpha5-251-g65c4417a2 From be95f9a98ab0419dd5553f15049673acd8af880d Mon Sep 17 00:00:00 2001 From: Daniel Otto de Mentock Date: Mon, 6 Dec 2021 23:49:12 +0000 Subject: [PATCH 02/74] Typehints for crystal --- python/damask/_crystal.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/python/damask/_crystal.py b/python/damask/_crystal.py index fd87a90d4..ea9effffe 100644 --- a/python/damask/_crystal.py +++ b/python/damask/_crystal.py @@ -1,3 +1,5 @@ +from typing import Union, Dict, List, Tuple + import numpy as np from . import util @@ -177,7 +179,7 @@ class Crystal(): @property - def standard_triangle(self): + def standard_triangle(self) -> Union[Dict[str, np.ndarray], None]: """ Corners of the standard triangle. @@ -238,7 +240,7 @@ class Crystal(): [-1., 0., 0.], [ 0., 1., 0.] ]), }} - return _basis.get(self.family,None) + return _basis.get(self.family, None) @property @@ -256,7 +258,7 @@ class Crystal(): @property - def basis_real(self): + def basis_real(self) -> np.ndarray: """ Return orthogonal real space crystal basis. @@ -280,7 +282,7 @@ class Crystal(): @property - def basis_reciprocal(self): + def basis_reciprocal(self) -> np.ndarray: """Return reciprocal (dual) crystal basis.""" return np.linalg.inv(self.basis_real.T) @@ -312,7 +314,7 @@ class Crystal(): + _lattice_points.get(self.lattice if self.lattice == 'hP' else \ self.lattice[-1],None),dtype=float) - def to_lattice(self,*,direction=None,plane=None): + def to_lattice(self, *, direction: np.ndarray = None, plane: np.ndarray = None) -> np.ndarray: """ Calculate lattice vector corresponding to crystal frame direction or plane normal. @@ -336,7 +338,7 @@ class Crystal(): return np.einsum('il,...l',basis,axis) - def to_frame(self,*,uvw=None,hkl=None): + def to_frame(self, *, uvw: np.ndarray = None, hkl: np.ndarray = None) -> np.ndarray: """ Calculate crystal frame vector along lattice direction [uvw] or plane normal (hkl). @@ -359,7 +361,7 @@ class Crystal(): return np.einsum('il,...l',basis,axis) - def kinematics(self,mode): + def kinematics(self, mode: str) -> Dict[str, List[np.ndarray]]: """ Return crystal kinematics systems. @@ -555,7 +557,7 @@ class Crystal(): 'plane': [m[:,3:6] for m in master]} - def relation_operations(self,model): + def relation_operations(self, model: str) -> Tuple[str, Rotation]: """ Crystallographic orientation relationships for phase transformations. From c56642f64a84e51e0d99652c4b55c9d6bc1f1fd5 Mon Sep 17 00:00:00 2001 From: Test User Date: Tue, 7 Dec 2021 08:48:55 +0100 Subject: [PATCH 03/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-253-ga4cd663fc --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 656df5ec5..caa7f1748 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-244-gb57351045 +v3.0.0-alpha5-253-ga4cd663fc From 722c9828a03d101186b74789e06c6ecba6c8809c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 7 Dec 2021 09:27:24 +0100 Subject: [PATCH 04/74] trying to fix errors like cannot lock ref 'refs/heads/master': is at 2d39da4d8d387b7dbb5b63fae67bff41a57543ce but expected 596a8ecd119fded5604e02043bfabfa7bea0ccfe --- .gitlab-ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 756237b3f..8ceb7575a 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -208,7 +208,7 @@ update_revision: - git checkout master - git pull - git merge $UPDATEDREV -s recursive -X ours # conflicts occur only for inconsistent state - - git push origin master # master is now tested version and has updated VERSION file + - git push - git checkout development - git pull - git merge master -s recursive -X ours -m "[skip ci] Merge branch 'master' into development" # only possible conflict is in VERSION file From 811879788af1aa856bd89e279e3d1d0fd664bd94 Mon Sep 17 00:00:00 2001 From: Test User Date: Tue, 7 Dec 2021 11:37:02 +0100 Subject: [PATCH 05/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-261-g722c9828a --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 2b61790a8..cd9e9a549 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-251-g65c4417a2 +v3.0.0-alpha5-261-g722c9828a From fbbe5ec006c302ef3456abb84d6f8e86c6175513 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 8 Dec 2021 21:04:22 +0000 Subject: [PATCH 06/74] fixed typo --- python/damask/_orientation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index bb14fd38b..e727c54ae 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -13,7 +13,7 @@ _parameter_doc = \ """ family : {'triclinic', 'monoclinic', 'orthorhombic', 'tetragonal', 'hexagonal', 'cubic'}, optional. Name of the crystal family. - Will be infered if 'lattice' is given. + Family will be inferred if 'lattice' is given. lattice : {'aP', 'mP', 'mS', 'oP', 'oS', 'oI', 'oF', 'tP', 'tI', 'hP', 'cP', 'cI', 'cF'}, optional. Name of the Bravais lattice in Pearson notation. a : float, optional From 63c76abed870fc3c9fdf7461b35a4160e77820e4 Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 9 Dec 2021 00:09:39 +0100 Subject: [PATCH 07/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-263-gfbbe5ec00 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index cd9e9a549..b22bd41df 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-261-g722c9828a +v3.0.0-alpha5-263-gfbbe5ec00 From 947bf946e1dfae6f278d61513cc4b7e092987ed7 Mon Sep 17 00:00:00 2001 From: Eureka Pai Date: Fri, 10 Dec 2021 11:31:26 -0500 Subject: [PATCH 08/74] added Colormap.at(fraction) to interpolate --- python/damask/_colormap.py | 28 +++++++++++++++++++++++++--- python/tests/test_Colormap.py | 7 +++++++ 2 files changed, 32 insertions(+), 3 deletions(-) diff --git a/python/damask/_colormap.py b/python/damask/_colormap.py index 233be569a..36bcf45c0 100644 --- a/python/damask/_colormap.py +++ b/python/damask/_colormap.py @@ -6,6 +6,7 @@ from pathlib import Path from typing import Sequence, Union, TextIO import numpy as np +import scipy.interpolate as interp import matplotlib as mpl if os.name == 'posix' and 'DISPLAY' not in os.environ: mpl.use('Agg') @@ -191,6 +192,28 @@ class Colormap(mpl.colors.ListedColormap): return Colormap.from_range(definition['low'],definition['high'],name,N) + def at(self, + fraction : float) -> np.ndarray: + """ + Interpolate color at fraction. + + Parameters + ---------- + fraction : float + Fractional coordinate to evaluate Colormap at. + + Returns + ------- + color : np.ndarray, shape(3) + RGB values of interpolated color. + + """ + return interp.interp1d(np.linspace(0,1,self.N), + self.colors, + axis=0, + assume_sorted=True)(fraction) + + def shade(self, field: np.ndarray, bounds: Sequence[float] = None, @@ -213,7 +236,6 @@ class Colormap(mpl.colors.ListedColormap): RGBA image of shaded data. """ - N = len(self.colors) mask = np.logical_not(np.isnan(field) if gap is None else \ np.logical_or (np.isnan(field), field == gap)) # mask NaN (and gap if present) @@ -227,7 +249,7 @@ class Colormap(mpl.colors.ListedColormap): return Image.fromarray( (np.dstack(( - self.colors[(np.round(np.clip((field-lo)/(hi-lo),0.0,1.0)*(N-1))).astype(np.uint16),:3], + self.colors[(np.round(np.clip((field-lo)/(hi-lo),0.0,1.0)*(self.N-1))).astype(np.uint16),:3], mask.astype(float) ) )*255 @@ -343,7 +365,7 @@ class Colormap(mpl.colors.ListedColormap): # ToDo: test in GOM GOM_str = '1 1 {name} 9 {name} '.format(name=self.name.replace(" ","_")) \ + '0 1 0 3 0 0 -1 9 \\ 0 0 0 255 255 255 0 0 255 ' \ - + f'30 NO_UNIT 1 1 64 64 64 255 1 0 0 0 0 0 0 3 0 {len(self.colors)}' \ + + f'30 NO_UNIT 1 1 64 64 64 255 1 0 0 0 0 0 0 3 0 {self.N}' \ + ' '.join([f' 0 {c[0]} {c[1]} {c[2]} 255 1' for c in reversed((self.colors*255).astype(int))]) \ + '\n' diff --git a/python/tests/test_Colormap.py b/python/tests/test_Colormap.py index ab9bcf92f..2c2c78560 100644 --- a/python/tests/test_Colormap.py +++ b/python/tests/test_Colormap.py @@ -139,6 +139,13 @@ class TestColormap: c += c assert (np.allclose(c.colors[:len(c.colors)//2],c.colors[len(c.colors)//2:])) + @pytest.mark.parametrize('N,cmap',[ + (8,'gray'), + (17,'gray'), + ]) + def test_at(self, N, cmap): + assert np.allclose(Colormap.from_predefined(cmap,N=N).at(0.5)[:3],0.5,rtol=0.005) + @pytest.mark.parametrize('bounds',[None,[2,10]]) def test_shade(self,ref_path,update,bounds): data = np.add(*np.indices((10, 11))) From 7bb6707b0e8af00f44a97cc4d7d9a2012b71ead9 Mon Sep 17 00:00:00 2001 From: Franz Roters Date: Fri, 10 Dec 2021 18:07:04 +0100 Subject: [PATCH 09/74] [skip ci] Web page was moved --- python/damask/_colormap.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/_colormap.py b/python/damask/_colormap.py index 36bcf45c0..6b8da4909 100644 --- a/python/damask/_colormap.py +++ b/python/damask/_colormap.py @@ -42,7 +42,7 @@ class Colormap(mpl.colors.ListedColormap): https://doi.org/10.1016/j.ijplas.2012.09.012 Matplotlib colormaps overview - https://matplotlib.org/tutorials/colors/colormaps.html + https://matplotlib.org/stable/tutorials/colors/colormaps.html """ From 2b6283bbe301048b140a966deed708e4e964c78f Mon Sep 17 00:00:00 2001 From: Eureka Pai Date: Fri, 10 Dec 2021 12:52:44 -0500 Subject: [PATCH 10/74] flexibility to return with alpha channel and/or as string --- python/damask/_colormap.py | 32 +++++++++++++++++++++++--------- python/tests/test_Colormap.py | 23 +++++++++++++++++------ 2 files changed, 40 insertions(+), 15 deletions(-) diff --git a/python/damask/_colormap.py b/python/damask/_colormap.py index 6b8da4909..97330baf3 100644 --- a/python/damask/_colormap.py +++ b/python/damask/_colormap.py @@ -193,25 +193,39 @@ class Colormap(mpl.colors.ListedColormap): def at(self, - fraction : float) -> np.ndarray: + fraction : Union[float,Sequence[float]], + with_alpha : bool = False, + as_string : bool = False) -> np.ndarray: """ Interpolate color at fraction. Parameters ---------- - fraction : float - Fractional coordinate to evaluate Colormap at. + fraction : float or sequence of float + Fractional coordinate(s) to evaluate Colormap at. + with_alpha : bool + Provide opacity as fourth value. + as_string : bool + Encode color as 'rgb(r,g,b)' or 'rgba(r,g,b,a)'. Returns ------- - color : np.ndarray, shape(3) - RGB values of interpolated color. + color : np.ndarray, shape(3 or 4), or string + RGB(A) values of interpolated color. """ - return interp.interp1d(np.linspace(0,1,self.N), - self.colors, - axis=0, - assume_sorted=True)(fraction) + def _stringify(color): + d = color.shape[-1] + c = np.array([f'rgba({c[0]},{c[1]},{c[2]},{c[3]})' if d==4 + else f'rgb({c[0]},{c[1]},{c[2]})' for c in np.atleast_2d(color)]) + return c if len(c)>1 else c[0] + + + color = interp.interp1d(np.linspace(0,1,self.N), + self.colors, + axis=0, + assume_sorted=True)(fraction)[...,:4 if with_alpha else 3] + return _stringify(color) if as_string else color def shade(self, diff --git a/python/tests/test_Colormap.py b/python/tests/test_Colormap.py index 2c2c78560..d411ca39d 100644 --- a/python/tests/test_Colormap.py +++ b/python/tests/test_Colormap.py @@ -139,12 +139,23 @@ class TestColormap: c += c assert (np.allclose(c.colors[:len(c.colors)//2],c.colors[len(c.colors)//2:])) - @pytest.mark.parametrize('N,cmap',[ - (8,'gray'), - (17,'gray'), - ]) - def test_at(self, N, cmap): - assert np.allclose(Colormap.from_predefined(cmap,N=N).at(0.5)[:3],0.5,rtol=0.005) + @pytest.mark.parametrize('N,cmap,at,with_alpha,result',[ + (8,'gray',0.5,False,[0.5,0.5,0.5]), + (17,'gray',0.5,False,[0.5,0.5,0.5]), + (17,'gray',[0.5,0.75],False,[[0.5,0.5,0.5],[0.75,0.75,0.75]]), + ]) + def test_at_value(self, N, cmap, at, with_alpha, result): + assert np.allclose(Colormap.from_predefined(cmap,N=N).at(at,with_alpha=with_alpha), + result, + rtol=0.005) + + @pytest.mark.parametrize('N,cmap,at,with_alpha,result',[ + (8,'gray',0.5,False,'rgb(0.5,0.5,0.5)'), + (8,'gray',0.5,True,'rgba(0.5,0.5,0.5,1.0)'), + (8,'gray',[0.5,0.25],True,['rgba(0.5,0.5,0.5,1.0)','rgba(0.25,0.25,0.25,1.0)']), + ]) + def test_at_string(self, N, cmap, at, with_alpha, result): + assert np.all(Colormap.from_predefined(cmap,N=N).at(at,with_alpha=with_alpha,as_string=True) == result) @pytest.mark.parametrize('bounds',[None,[2,10]]) def test_shade(self,ref_path,update,bounds): From 7d7d0c26591ac067fb2e385e9468d83907b6abdc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Dec 2021 08:49:30 +0100 Subject: [PATCH 11/74] only local variable are good variables --- src/grid/DAMASK_grid.f90 | 9 +++++++++ src/grid/spectral_utilities.f90 | 15 +-------------- 2 files changed, 10 insertions(+), 14 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index de324e0f5..b1da0c2a2 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -43,6 +43,15 @@ program DAMASK_grid logical :: estimate_rate !< follow trajectory of former loadcase end type tLoadCase +!-------------------------------------------------------------------------------------------------- +! field labels information + enum, bind(c); enumerator :: & + FIELD_UNDEFINED_ID, & + FIELD_MECH_ID, & + FIELD_THERMAL_ID, & + FIELD_DAMAGE_ID + end enum + integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:) !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 29777db3b..69b1d8fda 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -28,15 +28,6 @@ module spectral_utilities include 'fftw3-mpi.f03' -!-------------------------------------------------------------------------------------------------- -! field labels information - enum, bind(c); enumerator :: & - FIELD_UNDEFINED_ID, & - FIELD_MECH_ID, & - FIELD_THERMAL_ID, & - FIELD_DAMAGE_ID - end enum - !-------------------------------------------------------------------------------------------------- ! grid related information information real(pReal), protected, public :: wgt !< weighting factor 1/Nelems @@ -139,11 +130,7 @@ module spectral_utilities utilities_calculateRate, & utilities_forwardField, & utilities_updateCoords, & - utilities_saveReferenceStiffness, & - FIELD_UNDEFINED_ID, & - FIELD_MECH_ID, & - FIELD_THERMAL_ID, & - FIELD_DAMAGE_ID + utilities_saveReferenceStiffness contains From f51633d43a32a98b4117009971d7eb37a6783692 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Dec 2021 09:01:42 +0100 Subject: [PATCH 12/74] forall is deprecated do concurrent is the successor but ifort had problems and generated faulty code --- src/grid/spectral_utilities.f90 | 40 ++++++++++++++++++++------------- 1 file changed, 24 insertions(+), 16 deletions(-) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 69b1d8fda..f40025039 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -375,21 +375,24 @@ subroutine utilities_updateGamma(C) gamma_hat = cmplx(0.0_pReal,0.0_pReal,pReal) ! for the singular point and any non invertible A do k = grid3Offset+1, grid3Offset+grid3; do j = 1, grid(2); do i = 1, grid1Red if (any([i,j,k] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 - forall(l = 1:3, m = 1:3) & + do concurrent (l = 1:3, m = 1:3) xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-grid3Offset))*xi1st(m,i,j,k-grid3Offset) - forall(l = 1:3, m = 1:3) & + end do + do concurrent(l = 1:3, m = 1:3) temp33_complex(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx) + end do A(1:3,1:3) = temp33_complex%re; A(4:6,4:6) = temp33_complex%re A(1:3,4:6) = temp33_complex%im; A(4:6,1:3) = -temp33_complex%im if (abs(math_det33(A(1:3,1:3))) > 1e-16) then call math_invert(A_inv, err, A) temp33_complex = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal) - forall(l=1:3, m=1:3, n=1:3, o=1:3) & + do concurrent(l=1:3, m=1:3, n=1:3, o=1:3) gamma_hat(l,m,n,o,i,j,k-grid3Offset) = temp33_complex(l,n)* & conjg(-xi1st(o,i,j,k-grid3Offset))*xi1st(m,i,j,k-grid3Offset) - endif - endif - enddo; enddo; enddo + end do + end if + end if + end do; end do; end do endif end subroutine utilities_updateGamma @@ -492,32 +495,37 @@ subroutine utilities_fourierGammaConvolution(fieldAim) memoryEfficient: if (num%memory_efficient) then do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red if (any([i,j,k+grid3Offset] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 - forall(l = 1:3, m = 1:3) & + do concurrent(l = 1:3, m = 1:3) xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k) - forall(l = 1:3, m = 1:3) & + end do + do concurrent(l = 1:3, m = 1:3) temp33_complex(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx) + end do A(1:3,1:3) = temp33_complex%re; A(4:6,4:6) = temp33_complex%re A(1:3,4:6) = temp33_complex%im; A(4:6,1:3) = -temp33_complex%im if (abs(math_det33(A(1:3,1:3))) > 1e-16) then call math_invert(A_inv, err, A) temp33_complex = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal) - forall(l=1:3, m=1:3, n=1:3, o=1:3) & + do concurrent(l=1:3, m=1:3, n=1:3, o=1:3) gamma_hat(l,m,n,o,1,1,1) = temp33_complex(l,n)*conjg(-xi1st(o,i,j,k))*xi1st(m,i,j,k) + end do else gamma_hat(1:3,1:3,1:3,1:3,1,1,1) = cmplx(0.0_pReal,0.0_pReal,pReal) - endif - forall(l = 1:3, m = 1:3) & + end if + do concurrent(l = 1:3, m = 1:3) temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3,1,1,1)*tensorField_fourier(1:3,1:3,i,j,k)) + end do tensorField_fourier(1:3,1:3,i,j,k) = temp33_Complex - endif - enddo; enddo; enddo + end if + end do; end do; end do else memoryEfficient do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red - forall(l = 1:3, m = 1:3) & + do concurrent(l = 1:3, m = 1:3) temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k) * tensorField_fourier(1:3,1:3,i,j,k)) + end do tensorField_fourier(1:3,1:3,i,j,k) = temp33_Complex - enddo; enddo; enddo - endif memoryEfficient + end do; end do; end do + end if memoryEfficient if (grid3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal) From d6ba73d9e28589f9fb4a4e860e679edb19bd7dd2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Dec 2021 09:54:46 +0100 Subject: [PATCH 13/74] consistent names --- src/phase_mechanical.f90 | 38 +++++++-------- src/phase_mechanical_eigen.f90 | 18 ++++---- src/phase_mechanical_elastic.f90 | 2 +- src/phase_mechanical_plastic.f90 | 56 +++++++++++------------ src/phase_mechanical_plastic_nonlocal.f90 | 4 +- 5 files changed, 59 insertions(+), 59 deletions(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index db2faf914..37e6fe7bf 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -5,17 +5,17 @@ submodule(phase) mechanical enum, bind(c); enumerator :: & - PLASTICITY_UNDEFINED_ID, & - PLASTICITY_NONE_ID, & - PLASTICITY_ISOTROPIC_ID, & - PLASTICITY_PHENOPOWERLAW_ID, & - PLASTICITY_KINEHARDENING_ID, & - PLASTICITY_DISLOTWIN_ID, & - PLASTICITY_DISLOTUNGSTEN_ID, & - PLASTICITY_NONLOCAL_ID, & - KINEMATICS_UNDEFINED_ID, & - KINEMATICS_CLEAVAGE_OPENING_ID, & - KINEMATICS_THERMAL_EXPANSION_ID + PLASTIC_UNDEFINED_ID, & + PLASTIC_NONE_ID, & + PLASTIC_ISOTROPIC_ID, & + PLASTIC_PHENOPOWERLAW_ID, & + PLASTIC_KINEHARDENING_ID, & + PLASTIC_DISLOTWIN_ID, & + PLASTIC_DISLOTUNGSTEN_ID, & + PLASTIC_NONLOCAL_ID, & + EIGEN_UNDEFINED_ID, & + EIGEN_CLEAVAGE_OPENING_ID, & + EIGEN_THERMAL_EXPANSION_ID end enum type(tTensorContainer), dimension(:), allocatable :: & @@ -37,7 +37,7 @@ submodule(phase) mechanical phase_mechanical_S0 - integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: & + integer(kind(PLASTIC_undefined_ID)), dimension(:), allocatable :: & phase_plasticity !< plasticity of each phase integer :: phase_plasticity_maxSizeDotState @@ -291,7 +291,7 @@ module subroutine mechanical_init(phases) call elastic_init(phases) allocate(plasticState(phases%length)) - allocate(phase_plasticity(phases%length),source = PLASTICITY_undefined_ID) + allocate(phase_plasticity(phases%length),source = PLASTIC_UNDEFINED_ID) call plastic_init() do ph = 1,phases%length plasticState(ph)%state0 = plasticState(ph)%state @@ -340,22 +340,22 @@ module subroutine mechanical_results(group,ph) select case(phase_plasticity(ph)) - case(PLASTICITY_ISOTROPIC_ID) + case(PLASTIC_ISOTROPIC_ID) call plastic_isotropic_results(ph,group//'mechanical/') - case(PLASTICITY_PHENOPOWERLAW_ID) + case(PLASTIC_PHENOPOWERLAW_ID) call plastic_phenopowerlaw_results(ph,group//'mechanical/') - case(PLASTICITY_KINEHARDENING_ID) + case(PLASTIC_KINEHARDENING_ID) call plastic_kinehardening_results(ph,group//'mechanical/') - case(PLASTICITY_DISLOTWIN_ID) + case(PLASTIC_DISLOTWIN_ID) call plastic_dislotwin_results(ph,group//'mechanical/') - case(PLASTICITY_DISLOTUNGSTEN_ID) + case(PLASTIC_DISLOTUNGSTEN_ID) call plastic_dislotungsten_results(ph,group//'mechanical/') - case(PLASTICITY_NONLOCAL_ID) + case(PLASTIC_NONLOCAL_ID) call plastic_nonlocal_results(ph,group//'mechanical/') end select diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index c73665d17..f0b86319c 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -3,9 +3,9 @@ submodule(phase:mechanical) eigen integer, dimension(:), allocatable :: & Nmodels - integer(kind(KINEMATICS_UNDEFINED_ID)), dimension(:,:), allocatable :: & + integer(kind(EIGEN_UNDEFINED_ID)), dimension(:,:), allocatable :: & model - integer(kind(KINEMATICS_UNDEFINED_ID)), dimension(:), allocatable :: & + integer(kind(EIGEN_UNDEFINED_ID)), dimension(:), allocatable :: & model_damage interface @@ -57,15 +57,15 @@ module subroutine eigen_init(phases) Nmodels(ph) = kinematics%length end do - allocate(model(maxval(Nmodels),phases%length), source = KINEMATICS_undefined_ID) + allocate(model(maxval(Nmodels),phases%length), source = EIGEN_undefined_ID) if (maxval(Nmodels) /= 0) then - where(thermalexpansion_init(maxval(Nmodels))) model = KINEMATICS_thermal_expansion_ID + where(thermalexpansion_init(maxval(Nmodels))) model = EIGEN_thermal_expansion_ID endif - allocate(model_damage(phases%length), source = KINEMATICS_UNDEFINED_ID) + allocate(model_damage(phases%length), source = EIGEN_UNDEFINED_ID) - where(damage_anisobrittle_init()) model_damage = KINEMATICS_cleavage_opening_ID + where(damage_anisobrittle_init()) model_damage = EIGEN_cleavage_opening_ID end subroutine eigen_init @@ -173,7 +173,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & plasticType: select case (phase_plasticity(ph)) - case (PLASTICITY_isotropic_ID) plasticType + case (PLASTIC_isotropic_ID) plasticType call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,ph,en) Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS @@ -182,7 +182,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & KinematicsLoop: do k = 1, Nmodels(ph) kinematicsType: select case (model(k,ph)) - case (KINEMATICS_thermal_expansion_ID) kinematicsType + case (EIGEN_thermal_expansion_ID) kinematicsType call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,en) Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS @@ -191,7 +191,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & end do KinematicsLoop select case (model_damage(ph)) - case (KINEMATICS_cleavage_opening_ID) + case (EIGEN_cleavage_opening_ID) call damage_anisobrittle_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, en) Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index ea113ddfb..acbbac236 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -223,7 +223,7 @@ module function phase_homogenizedC66(ph,en) result(C) plasticType: select case (phase_plasticity(ph)) - case (PLASTICITY_DISLOTWIN_ID) plasticType + case (PLASTIC_DISLOTWIN_ID) plasticType C = plastic_dislotwin_homogenizedC(ph,en) case default plasticType C = elastic_C66(ph,en) diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index b2a9bdcc4..4951712fd 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -224,15 +224,15 @@ module subroutine plastic_init print'(/,1x,a)', '<<<+- phase:mechanical:plastic init -+>>>' - where(plastic_none_init()) phase_plasticity = PLASTICITY_NONE_ID - where(plastic_isotropic_init()) phase_plasticity = PLASTICITY_ISOTROPIC_ID - where(plastic_phenopowerlaw_init()) phase_plasticity = PLASTICITY_PHENOPOWERLAW_ID - where(plastic_kinehardening_init()) phase_plasticity = PLASTICITY_KINEHARDENING_ID - where(plastic_dislotwin_init()) phase_plasticity = PLASTICITY_DISLOTWIN_ID - where(plastic_dislotungsten_init()) phase_plasticity = PLASTICITY_DISLOTUNGSTEN_ID - where(plastic_nonlocal_init()) phase_plasticity = PLASTICITY_NONLOCAL_ID + where(plastic_none_init()) phase_plasticity = PLASTIC_NONE_ID + where(plastic_isotropic_init()) phase_plasticity = PLASTIC_ISOTROPIC_ID + where(plastic_phenopowerlaw_init()) phase_plasticity = PLASTIC_PHENOPOWERLAW_ID + where(plastic_kinehardening_init()) phase_plasticity = PLASTIC_KINEHARDENING_ID + where(plastic_dislotwin_init()) phase_plasticity = PLASTIC_DISLOTWIN_ID + where(plastic_dislotungsten_init()) phase_plasticity = PLASTIC_DISLOTUNGSTEN_ID + where(plastic_nonlocal_init()) phase_plasticity = PLASTIC_NONLOCAL_ID - if (any(phase_plasticity == PLASTICITY_undefined_ID)) call IO_error(201) + if (any(phase_plasticity == PLASTIC_undefined_ID)) call IO_error(201) end subroutine plastic_init @@ -262,7 +262,7 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & i, j - if (phase_plasticity(ph) == PLASTICITY_NONE_ID) then + if (phase_plasticity(ph) == PLASTIC_NONE_ID) then Lp = 0.0_pReal dLp_dFi = 0.0_pReal dLp_dS = 0.0_pReal @@ -272,22 +272,22 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & plasticType: select case (phase_plasticity(ph)) - case (PLASTICITY_ISOTROPIC_ID) plasticType + case (PLASTIC_ISOTROPIC_ID) plasticType call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) - case (PLASTICITY_PHENOPOWERLAW_ID) plasticType + case (PLASTIC_PHENOPOWERLAW_ID) plasticType call phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) - case (PLASTICITY_KINEHARDENING_ID) plasticType + case (PLASTIC_KINEHARDENING_ID) plasticType call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) - case (PLASTICITY_NONLOCAL_ID) plasticType + case (PLASTIC_NONLOCAL_ID) plasticType call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) - case (PLASTICITY_DISLOTWIN_ID) plasticType + case (PLASTIC_DISLOTWIN_ID) plasticType call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) - case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType + case (PLASTIC_DISLOTUNGSTEN_ID) plasticType call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) end select plasticType @@ -321,28 +321,28 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken) logical :: broken - if (phase_plasticity(ph) /= PLASTICITY_NONE_ID) then + if (phase_plasticity(ph) /= PLASTIC_NONE_ID) then Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),& phase_mechanical_Fi(ph)%data(1:3,1:3,en)),phase_mechanical_S(ph)%data(1:3,1:3,en)) plasticType: select case (phase_plasticity(ph)) - case (PLASTICITY_ISOTROPIC_ID) plasticType + case (PLASTIC_ISOTROPIC_ID) plasticType call isotropic_dotState(Mp,ph,en) - case (PLASTICITY_PHENOPOWERLAW_ID) plasticType + case (PLASTIC_PHENOPOWERLAW_ID) plasticType call phenopowerlaw_dotState(Mp,ph,en) - case (PLASTICITY_KINEHARDENING_ID) plasticType + case (PLASTIC_KINEHARDENING_ID) plasticType call plastic_kinehardening_dotState(Mp,ph,en) - case (PLASTICITY_DISLOTWIN_ID) plasticType + case (PLASTIC_DISLOTWIN_ID) plasticType call dislotwin_dotState(Mp,thermal_T(ph,en),ph,en) - case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType + case (PLASTIC_DISLOTUNGSTEN_ID) plasticType call dislotungsten_dotState(Mp,thermal_T(ph,en),ph,en) - case (PLASTICITY_NONLOCAL_ID) plasticType + case (PLASTIC_NONLOCAL_ID) plasticType call nonlocal_dotState(Mp,thermal_T(ph,en),subdt,ph,en,ip,el) end select plasticType end if @@ -372,13 +372,13 @@ module subroutine plastic_dependentState(co, ip, el) plasticType: select case (phase_plasticity(ph)) - case (PLASTICITY_DISLOTWIN_ID) plasticType + case (PLASTIC_DISLOTWIN_ID) plasticType call dislotwin_dependentState(thermal_T(ph,en),ph,en) - case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType + case (PLASTIC_DISLOTUNGSTEN_ID) plasticType call dislotungsten_dependentState(ph,en) - case (PLASTICITY_NONLOCAL_ID) plasticType + case (PLASTIC_NONLOCAL_ID) plasticType call nonlocal_dependentState(ph,en,ip,el) end select plasticType @@ -406,7 +406,7 @@ module function plastic_deltaState(ph, en) result(broken) broken = .false. select case (phase_plasticity(ph)) - case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID) + case (PLASTIC_NONLOCAL_ID,PLASTIC_KINEHARDENING_ID) Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),& phase_mechanical_Fi(ph)%data(1:3,1:3,en)),& @@ -414,10 +414,10 @@ module function plastic_deltaState(ph, en) result(broken) plasticType: select case (phase_plasticity(ph)) - case (PLASTICITY_KINEHARDENING_ID) plasticType + case (PLASTIC_KINEHARDENING_ID) plasticType call plastic_kinehardening_deltaState(Mp,ph,en) - case (PLASTICITY_NONLOCAL_ID) plasticType + case (PLASTIC_NONLOCAL_ID) plasticType call plastic_nonlocal_deltaState(Mp,ph,en) end select plasticType diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 6bd648eaa..c762039d4 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -1277,7 +1277,7 @@ function rhoDotFlux(timestep,ph,en,ip,el) !* The entering flux from my neighbor will be distributed on my slip systems according to the !* compatibility if (neighbor_n > 0) then - if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTICITY_NONLOCAL_ID .and. & + if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTIC_NONLOCAL_ID .and. & any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) then forall (s = 1:ns, t = 1:4) @@ -1323,7 +1323,7 @@ function rhoDotFlux(timestep,ph,en,ip,el) !* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density. !* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations. if (opposite_n > 0) then - if (phase_plasticity(material_phaseAt(1,opposite_el)) == PLASTICITY_NONLOCAL_ID) then + if (phase_plasticity(material_phaseAt(1,opposite_el)) == PLASTIC_NONLOCAL_ID) then normal_me2neighbor_defConf = math_det33(Favg) & * matmul(math_inv33(transpose(Favg)),IPareaNormal(1:3,n,ip,el)) ! normal of the interface in (average) deformed configuration (pointing en => neighbor) From 6ba2a08e5a5fea121c1bb8693213f9bad249d4ca Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Dec 2021 11:37:17 +0100 Subject: [PATCH 14/74] easier to read --- src/grid/grid_damage_spectral.f90 | 5 ++--- src/grid/grid_thermal_spectral.f90 | 5 ++--- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 187a6f552..ef229368e 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -167,7 +167,7 @@ function grid_damage_spectral_solution(Delta_t) result(solution) integer :: i, j, k, ce type(tSolutionState) :: solution PetscInt :: devNull - PetscReal :: phi_min, phi_max, stagNorm, solnNorm + PetscReal :: phi_min, phi_max, stagNorm PetscErrorCode :: ierr SNESConvergedReason :: reason @@ -189,9 +189,8 @@ function grid_damage_spectral_solution(Delta_t) result(solution) solution%iterationsNeeded = totalIter end if stagNorm = maxval(abs(phi_current - phi_stagInc)) - solnNorm = maxval(abs(phi_current)) call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD,ierr) - solution%stagConverged = stagNorm < max(num%eps_damage_atol, num%eps_damage_rtol*solnNorm) + solution%stagConverged = stagNorm < max(num%eps_damage_atol, num%eps_damage_rtol*maxval(phi_current)) call MPI_Allreduce(MPI_IN_PLACE,solution%stagConverged,1,MPI_LOGICAL,MPI_LAND,MPI_COMM_WORLD,ierr) phi_stagInc = phi_current diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 868d3101e..67c7ba1c3 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -162,7 +162,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution) integer :: i, j, k, ce type(tSolutionState) :: solution PetscInt :: devNull - PetscReal :: T_min, T_max, stagNorm, solnNorm + PetscReal :: T_min, T_max, stagNorm PetscErrorCode :: ierr SNESConvergedReason :: reason @@ -184,9 +184,8 @@ function grid_thermal_spectral_solution(Delta_t) result(solution) solution%iterationsNeeded = totalIter end if stagNorm = maxval(abs(T_current - T_stagInc)) - solnNorm = maxval(abs(T_current)) call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD,ierr) - solution%stagConverged = stagNorm < max(num%eps_thermal_atol, num%eps_thermal_rtol*solnNorm) + solution%stagConverged = stagNorm < max(num%eps_thermal_atol, num%eps_thermal_rtol*maxval(T_current)) call MPI_Allreduce(MPI_IN_PLACE,solution%stagConverged,1,MPI_LOGICAL,MPI_LAND,MPI_COMM_WORLD,ierr) T_stagInc = T_current From 5ef918e75b1e0d3cdfb8805e5a5389fa512f040b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Dec 2021 14:55:56 +0100 Subject: [PATCH 15/74] PETSc 3.16.2 is out --- .github/workflows/Fortran.yml | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/.github/workflows/Fortran.yml b/.github/workflows/Fortran.yml index 0d11918bc..f852893ea 100644 --- a/.github/workflows/Fortran.yml +++ b/.github/workflows/Fortran.yml @@ -48,17 +48,17 @@ jobs: uses: actions/cache@v2 with: path: download - key: petsc-3.16.0.tar.gz + key: petsc-3.16.2.tar.gz - name: PETSc - Download if: steps.petsc-download.outputs.cache-hit != 'true' run: | - wget -q https://ftp.mcs.anl.gov/pub/petsc/release-snapshots/petsc-3.16.0.tar.gz -P download + wget -q https://ftp.mcs.anl.gov/pub/petsc/release-snapshots/petsc-3.16.2.tar.gz -P download - name: PETSc - Prepare run: | - tar -xf download/petsc-3.16.0.tar.gz -C . - export PETSC_DIR=${PWD}/petsc-3.16.0 + tar -xf download/petsc-3.16.2.tar.gz -C . + export PETSC_DIR=${PWD}/petsc-3.16.2 export PETSC_ARCH=gcc${GCC_V} printenv >> $GITHUB_ENV @@ -66,13 +66,13 @@ jobs: id: petsc-install uses: actions/cache@v2 with: - path: petsc-3.16.0 - key: petsc-3.16.0-${{ matrix.os }}-gcc${{ matrix.gcc_v }}-${{ hashFiles('**/petscversion.h') }} + path: petsc-3.16.2 + key: petsc-3.16.2-${{ matrix.os }}-gcc${{ matrix.gcc_v }}-${{ hashFiles('**/petscversion.h') }} - name: PETSc - Install (Linux) if: contains( matrix.os, 'ubuntu') run: | - cd petsc-3.16.0 + cd petsc-3.16.2 ./configure --with-fc=gfortran --with-cc=gcc --with-cxx=g++ \ --download-mpich --download-fftw --download-hdf5 --download-hdf5-fortran-bindings=1 --download-zlib \ --with-mpi-f90module-visibility=0 @@ -81,7 +81,7 @@ jobs: - name: PETSc - Install (macOS) if: contains( matrix.os, 'macos') run: | - cd petsc-3.16.0 + cd petsc-3.16.2 ./configure --with-fc=gfortran-${GCC_V} --with-cc=gcc-${GCC_V} --with-cxx=g++-${GCC_V} \ --download-openmpi --download-fftw --download-hdf5 --download-hdf5-fortran-bindings=1 --download-zlib make all @@ -132,17 +132,17 @@ jobs: uses: actions/cache@v2 with: path: download - key: petsc-3.16.0.tar.gz + key: petsc-3.16.2.tar.gz - name: PETSc - Download if: steps.petsc-download.outputs.cache-hit != 'true' run: | - wget -q https://ftp.mcs.anl.gov/pub/petsc/release-snapshots/petsc-3.16.0.tar.gz -P download + wget -q https://ftp.mcs.anl.gov/pub/petsc/release-snapshots/petsc-3.16.2.tar.gz -P download - name: PETSc - Prepare run: | - tar -xf download/petsc-3.16.0.tar.gz -C . - export PETSC_DIR=${PWD}/petsc-3.16.0 + tar -xf download/petsc-3.16.2.tar.gz -C . + export PETSC_DIR=${PWD}/petsc-3.16.2 export PETSC_ARCH=intel-${INTEL_V} printenv >> $GITHUB_ENV @@ -150,13 +150,13 @@ jobs: id: petsc-install uses: actions/cache@v2 with: - path: petsc-3.16.0 - key: petsc-3.16.0-intel-${{ matrix.intel_v }}-${{ hashFiles('**/petscversion.h') }} + path: petsc-3.16.2 + key: petsc-3.16.2-intel-${{ matrix.intel_v }}-${{ hashFiles('**/petscversion.h') }} - name: PETSc - Install (classic) if: contains( matrix.intel_v, 'classic') run: | - cd petsc-3.16.0 + cd petsc-3.16.2 ./configure --with-fc=mpiifort --with-cc=mpiicc --with-cxx=mpiicpc \ --download-fftw --download-hdf5 --download-hdf5-fortran-bindings=1 --download-zlib make all @@ -164,7 +164,7 @@ jobs: - name: PETSc - Install (LLVM) if: contains( matrix.intel_v, 'llvm') run: | - cd petsc-3.16.0 + cd petsc-3.16.2 ./configure --with-fc=mpiifort --with-cc="mpiicc -cc=icx" --with-cxx="mpiicpc -cxx=icpx" \ --download-fftw --download-hdf5 --download-hdf5-fortran-bindings=1 --download-zlib make all From ad5262a7a8698e00ec60c4c142032b31fa186c2c Mon Sep 17 00:00:00 2001 From: Test User Date: Sat, 11 Dec 2021 22:21:22 +0100 Subject: [PATCH 16/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-270-g6e36c4c30 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index b22bd41df..1a4f13343 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-263-gfbbe5ec00 +v3.0.0-alpha5-270-g6e36c4c30 From 3192a31e1ea14c716b23a5b325c7cc2980bb22c9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Dec 2021 23:08:18 +0100 Subject: [PATCH 17/74] single source of truth --- .github/workflows/Fortran.yml | 43 ++++++++++++++++++----------------- 1 file changed, 22 insertions(+), 21 deletions(-) diff --git a/.github/workflows/Fortran.yml b/.github/workflows/Fortran.yml index f852893ea..5cc241e00 100644 --- a/.github/workflows/Fortran.yml +++ b/.github/workflows/Fortran.yml @@ -2,11 +2,12 @@ name: Grid and Mesh Solver on: [push] env: - HOMEBREW_NO_ANALYTICS: "ON" # Make Homebrew installation a little quicker - HOMEBREW_NO_AUTO_UPDATE: "ON" - HOMEBREW_NO_BOTTLE_SOURCE_FALLBACK: "ON" - HOMEBREW_NO_GITHUB_API: "ON" - HOMEBREW_NO_INSTALL_CLEANUP: "ON" + PETSC_VERSION: '3.16.2' + HOMEBREW_NO_ANALYTICS: 'ON' # Make Homebrew installation a little quicker + HOMEBREW_NO_AUTO_UPDATE: 'ON' + HOMEBREW_NO_BOTTLE_SOURCE_FALLBACK: 'ON' + HOMEBREW_NO_GITHUB_API: 'ON' + HOMEBREW_NO_INSTALL_CLEANUP: 'ON' jobs: @@ -48,17 +49,17 @@ jobs: uses: actions/cache@v2 with: path: download - key: petsc-3.16.2.tar.gz + key: petsc-${{ env.PETSC_VERSION }}.tar.gz - name: PETSc - Download if: steps.petsc-download.outputs.cache-hit != 'true' run: | - wget -q https://ftp.mcs.anl.gov/pub/petsc/release-snapshots/petsc-3.16.2.tar.gz -P download + wget -q https://ftp.mcs.anl.gov/pub/petsc/release-snapshots/petsc-${PETSC_VERSION}.tar.gz -P download - name: PETSc - Prepare run: | - tar -xf download/petsc-3.16.2.tar.gz -C . - export PETSC_DIR=${PWD}/petsc-3.16.2 + tar -xf download/petsc-${PETSC_VERSION}.tar.gz -C . + export PETSC_DIR=${PWD}/petsc-${PETSC_VERSION} export PETSC_ARCH=gcc${GCC_V} printenv >> $GITHUB_ENV @@ -66,13 +67,13 @@ jobs: id: petsc-install uses: actions/cache@v2 with: - path: petsc-3.16.2 - key: petsc-3.16.2-${{ matrix.os }}-gcc${{ matrix.gcc_v }}-${{ hashFiles('**/petscversion.h') }} + path: petsc-${{ env.PETSC_VERSION }} + key: petsc-${{ env.PETSC_VERSION }}-${{ matrix.os }}-gcc${{ matrix.gcc_v }}-${{ hashFiles('**/petscversion.h') }} - name: PETSc - Install (Linux) if: contains( matrix.os, 'ubuntu') run: | - cd petsc-3.16.2 + cd petsc-${PETSC_VERSION} ./configure --with-fc=gfortran --with-cc=gcc --with-cxx=g++ \ --download-mpich --download-fftw --download-hdf5 --download-hdf5-fortran-bindings=1 --download-zlib \ --with-mpi-f90module-visibility=0 @@ -81,7 +82,7 @@ jobs: - name: PETSc - Install (macOS) if: contains( matrix.os, 'macos') run: | - cd petsc-3.16.2 + cd petsc-${PETSC_VERSION} ./configure --with-fc=gfortran-${GCC_V} --with-cc=gcc-${GCC_V} --with-cxx=g++-${GCC_V} \ --download-openmpi --download-fftw --download-hdf5 --download-hdf5-fortran-bindings=1 --download-zlib make all @@ -132,17 +133,17 @@ jobs: uses: actions/cache@v2 with: path: download - key: petsc-3.16.2.tar.gz + key: petsc-${{ env.PETSC_VERSION }}.tar.gz - name: PETSc - Download if: steps.petsc-download.outputs.cache-hit != 'true' run: | - wget -q https://ftp.mcs.anl.gov/pub/petsc/release-snapshots/petsc-3.16.2.tar.gz -P download + wget -q https://ftp.mcs.anl.gov/pub/petsc/release-snapshots/petsc-${PETSC_VERSION}.tar.gz -P download - name: PETSc - Prepare run: | - tar -xf download/petsc-3.16.2.tar.gz -C . - export PETSC_DIR=${PWD}/petsc-3.16.2 + tar -xf download/petsc-${PETSC_VERSION}.tar.gz -C . + export PETSC_DIR=${PWD}/petsc-${PETSC_VERSION} export PETSC_ARCH=intel-${INTEL_V} printenv >> $GITHUB_ENV @@ -150,13 +151,13 @@ jobs: id: petsc-install uses: actions/cache@v2 with: - path: petsc-3.16.2 - key: petsc-3.16.2-intel-${{ matrix.intel_v }}-${{ hashFiles('**/petscversion.h') }} + path: petsc-${{ env.PETSC_VERSION }} + key: petsc-${{ env.PETSC_VERSION }}-intel-${{ matrix.intel_v }}-${{ hashFiles('**/petscversion.h') }} - name: PETSc - Install (classic) if: contains( matrix.intel_v, 'classic') run: | - cd petsc-3.16.2 + cd petsc-${PETSC_VERSION} ./configure --with-fc=mpiifort --with-cc=mpiicc --with-cxx=mpiicpc \ --download-fftw --download-hdf5 --download-hdf5-fortran-bindings=1 --download-zlib make all @@ -164,7 +165,7 @@ jobs: - name: PETSc - Install (LLVM) if: contains( matrix.intel_v, 'llvm') run: | - cd petsc-3.16.2 + cd petsc-${PETSC_VERSION} ./configure --with-fc=mpiifort --with-cc="mpiicc -cc=icx" --with-cxx="mpiicpc -cxx=icpx" \ --download-fftw --download-hdf5 --download-hdf5-fortran-bindings=1 --download-zlib make all From a89eff2991a16b116010fe0a9e7297d83ba2d5fd Mon Sep 17 00:00:00 2001 From: Test User Date: Sun, 12 Dec 2021 02:03:10 +0100 Subject: [PATCH 18/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-272-g3192a31e1 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 1a4f13343..27dd17c2c 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-270-g6e36c4c30 +v3.0.0-alpha5-272-g3192a31e1 From b41d9411a7d8257017bdd46cf709c0d8edc598e3 Mon Sep 17 00:00:00 2001 From: Nikhil Prabhu Date: Mon, 13 Dec 2021 16:30:46 +0100 Subject: [PATCH 19/74] Temperature-dependent Elastic Constants for Silver --- .../phase/mechanical/elastic/Hooke_Ag.yaml | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 examples/config/phase/mechanical/elastic/Hooke_Ag.yaml diff --git a/examples/config/phase/mechanical/elastic/Hooke_Ag.yaml b/examples/config/phase/mechanical/elastic/Hooke_Ag.yaml new file mode 100644 index 000000000..e968b567e --- /dev/null +++ b/examples/config/phase/mechanical/elastic/Hooke_Ag.yaml @@ -0,0 +1,18 @@ +type: Hooke +references: + - Neighbours, J. R. & Alers, G. A. + Phys. Rev., American Physical Society, 1958, 111, 707-712 + https://journals.aps.org/pr/abstract/10.1103/PhysRev.111.707 + - Chang, Y. A. & Himmel, L. + Journal of Applied Physics, American Institute of Physics, 1966, 37, 3567-3572 + https://aip.scitation.org/doi/10.1063/1.1708903 +T_ref: 300 +C_11: 122.9e+9 +C_11, T: -313.5e+5 +C_11, T^2: -107.3e+2 +C_12: 915.5e+8 +C_12, T: -164.1e+5 +C_12, T^2: -681.6e+1 +C_44: 426.3e+8 +C_44, T: -180.5e+5 +C_44, T^2: -353.8e+1 From e330348c984dca7f409d14cd37a32e89369acdb8 Mon Sep 17 00:00:00 2001 From: Nikhil Prabhu Date: Mon, 13 Dec 2021 16:38:31 +0100 Subject: [PATCH 20/74] Temperature-dependent Elastic Constants for Aluminum --- .../phase/mechanical/elastic/Hooke_Al.yaml | 22 ++++++++++++++----- 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/examples/config/phase/mechanical/elastic/Hooke_Al.yaml b/examples/config/phase/mechanical/elastic/Hooke_Al.yaml index 6007f7695..fe2702c79 100644 --- a/examples/config/phase/mechanical/elastic/Hooke_Al.yaml +++ b/examples/config/phase/mechanical/elastic/Hooke_Al.yaml @@ -1,8 +1,18 @@ type: Hooke references: - - J. Vallin et al., - Journal of Applied Physics 35(6):1825-1826, 1964, - https://doi.org/10.1063/1.1713749 -C_11: 107.3e+9 -C_12: 60.8e+9 -C_44: 28.3e+9 + - Kamm, G. N. & Alers, G. A. + Journal of Applied Physics, American Institute of Physics, 1964, 35, 327-330 + https://doi.org/10.1063/1.1713309 + - Gerlich, D. & Fisher, E. S. + Journal of Physics and Chemistry of Solids, 1969, 30, 1197-1205 + https://doi.org/10.1016/0022-3697(69)90377-1 +T_ref: 300 +C_11: 106.1e+9 +C_11, T: -359.3e+5 +C_11, T^2: -152.7e+2 +C_12: 578.3e+8 +C_12, T: -781.6e+4 +C_12, T^2: -551.3e+1 +C_44: 243.1e+8 +C_44, T: -142.9e+5 +C_44, T^2: -404.6e+1 From d3a8b3d3d1922e917a54ce98425b13ca7381c427 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 13 Dec 2021 16:46:01 +0100 Subject: [PATCH 21/74] polishing --- .../phase/mechanical/elastic/Hooke_Al.yaml | 28 +++++++++++-------- 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/examples/config/phase/mechanical/elastic/Hooke_Al.yaml b/examples/config/phase/mechanical/elastic/Hooke_Al.yaml index fe2702c79..f4c266d54 100644 --- a/examples/config/phase/mechanical/elastic/Hooke_Al.yaml +++ b/examples/config/phase/mechanical/elastic/Hooke_Al.yaml @@ -1,18 +1,22 @@ type: Hooke references: - - Kamm, G. N. & Alers, G. A. - Journal of Applied Physics, American Institute of Physics, 1964, 35, 327-330 + - G.N. Kamm and G.A. Alers, + Journal of Applied Physics 35:327-330, 1964, https://doi.org/10.1063/1.1713309 - - Gerlich, D. & Fisher, E. S. - Journal of Physics and Chemistry of Solids, 1969, 30, 1197-1205 + - D. Gerlich and E.S. Fisher, + Journal of Physics and Chemistry of Solids 30:1197-1205, 1969 https://doi.org/10.1016/0022-3697(69)90377-1 + T_ref: 300 + C_11: 106.1e+9 -C_11, T: -359.3e+5 -C_11, T^2: -152.7e+2 -C_12: 578.3e+8 -C_12, T: -781.6e+4 -C_12, T^2: -551.3e+1 -C_44: 243.1e+8 -C_44, T: -142.9e+5 -C_44, T^2: -404.6e+1 +C_11,T: -359.3e+5 +C_11,T^2: -152.7e+2 + +C_12: 57.83e+9 +C_12,T: -781.6e+4 +C_12,T^2: -551.3e+1 + +C_44: 24.31e+9 +C_44,T: -142.9e+5 +C_44,T^2: -404.6e+1 From fd3e8532ee0379d6d4e72f400cc07ef0c3553a66 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 13 Dec 2021 16:50:47 +0100 Subject: [PATCH 22/74] polishing --- .../phase/mechanical/elastic/Hooke_Ag.yaml | 32 +++++++++++-------- 1 file changed, 18 insertions(+), 14 deletions(-) diff --git a/examples/config/phase/mechanical/elastic/Hooke_Ag.yaml b/examples/config/phase/mechanical/elastic/Hooke_Ag.yaml index e968b567e..9485a4bb2 100644 --- a/examples/config/phase/mechanical/elastic/Hooke_Ag.yaml +++ b/examples/config/phase/mechanical/elastic/Hooke_Ag.yaml @@ -1,18 +1,22 @@ type: Hooke references: - - Neighbours, J. R. & Alers, G. A. - Phys. Rev., American Physical Society, 1958, 111, 707-712 - https://journals.aps.org/pr/abstract/10.1103/PhysRev.111.707 - - Chang, Y. A. & Himmel, L. - Journal of Applied Physics, American Institute of Physics, 1966, 37, 3567-3572 - https://aip.scitation.org/doi/10.1063/1.1708903 + - J.R. Neighbours and G.A. Alers, + Physical Review 111:707-712, 1958, + https://doi.org/10.1103/PhysRev.111.707 + - Y.A. Chang and L. Himmel, + Journal of Applied Physics 37:3567-3572, 1966, + https://doi.org/10.1063/1.1708903 + T_ref: 300 + C_11: 122.9e+9 -C_11, T: -313.5e+5 -C_11, T^2: -107.3e+2 -C_12: 915.5e+8 -C_12, T: -164.1e+5 -C_12, T^2: -681.6e+1 -C_44: 426.3e+8 -C_44, T: -180.5e+5 -C_44, T^2: -353.8e+1 +C_11,T: -313.5e+5 +C_11,T^2: -107.3e+2 + +C_12: 91.55e+9 +C_12,T: -164.1e+5 +C_12,T^2: -681.6e+1 + +C_44: 42.63e+9 +C_44,T: -180.5e+5 +C_44,T^2: -353.8e+1 From 32595299c9e779d2c321a81d34c18f15858472af Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Mon, 13 Dec 2021 17:44:24 -0500 Subject: [PATCH 23/74] minimal clean API for Colormap.at --- python/damask/_colormap.py | 37 ++++++++++++++++--------------------- 1 file changed, 16 insertions(+), 21 deletions(-) diff --git a/python/damask/_colormap.py b/python/damask/_colormap.py index 97330baf3..5a13f440d 100644 --- a/python/damask/_colormap.py +++ b/python/damask/_colormap.py @@ -193,9 +193,7 @@ class Colormap(mpl.colors.ListedColormap): def at(self, - fraction : Union[float,Sequence[float]], - with_alpha : bool = False, - as_string : bool = False) -> np.ndarray: + fraction : Union[float,Sequence[float]]) -> np.ndarray: """ Interpolate color at fraction. @@ -203,29 +201,26 @@ class Colormap(mpl.colors.ListedColormap): ---------- fraction : float or sequence of float Fractional coordinate(s) to evaluate Colormap at. - with_alpha : bool - Provide opacity as fourth value. - as_string : bool - Encode color as 'rgb(r,g,b)' or 'rgba(r,g,b,a)'. Returns ------- - color : np.ndarray, shape(3 or 4), or string - RGB(A) values of interpolated color. + color : np.ndarray, shape(...,4) + RGBA values of interpolated color(s). + + Examples + -------- + >>> import damask + >>> cmap = damask.Colormap.from_predefined('gray') + >>> cmap.at(0.5) + array([0.5, 0.5, 0.5, 1. ]) + >>> 'rgb({},{},{})'.format(*cmap.at(0.5)) + 'rgb(0.5,0.5,0.5)' """ - def _stringify(color): - d = color.shape[-1] - c = np.array([f'rgba({c[0]},{c[1]},{c[2]},{c[3]})' if d==4 - else f'rgb({c[0]},{c[1]},{c[2]})' for c in np.atleast_2d(color)]) - return c if len(c)>1 else c[0] - - - color = interp.interp1d(np.linspace(0,1,self.N), - self.colors, - axis=0, - assume_sorted=True)(fraction)[...,:4 if with_alpha else 3] - return _stringify(color) if as_string else color + return interp.interp1d(np.linspace(0,1,self.N), + self.colors, + axis=0, + assume_sorted=True)(fraction) def shade(self, From 2b4487f42667c6572d6f86b56cac10690ce718da Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Mon, 13 Dec 2021 19:07:18 -0500 Subject: [PATCH 24/74] forgot to update test --- python/tests/test_Colormap.py | 20 ++++++-------------- 1 file changed, 6 insertions(+), 14 deletions(-) diff --git a/python/tests/test_Colormap.py b/python/tests/test_Colormap.py index d411ca39d..156907670 100644 --- a/python/tests/test_Colormap.py +++ b/python/tests/test_Colormap.py @@ -139,24 +139,16 @@ class TestColormap: c += c assert (np.allclose(c.colors[:len(c.colors)//2],c.colors[len(c.colors)//2:])) - @pytest.mark.parametrize('N,cmap,at,with_alpha,result',[ - (8,'gray',0.5,False,[0.5,0.5,0.5]), - (17,'gray',0.5,False,[0.5,0.5,0.5]), - (17,'gray',[0.5,0.75],False,[[0.5,0.5,0.5],[0.75,0.75,0.75]]), + @pytest.mark.parametrize('N,cmap,at,result',[ + (8,'gray',0.5,[0.5,0.5,0.5]), + (17,'gray',0.5,[0.5,0.5,0.5]), + (17,'gray',[0.5,0.75],[[0.5,0.5,0.5],[0.75,0.75,0.75]]), ]) - def test_at_value(self, N, cmap, at, with_alpha, result): - assert np.allclose(Colormap.from_predefined(cmap,N=N).at(at,with_alpha=with_alpha), + def test_at_value(self, N, cmap, at, result): + assert np.allclose(Colormap.from_predefined(cmap,N=N).at(at)[...,:3], result, rtol=0.005) - @pytest.mark.parametrize('N,cmap,at,with_alpha,result',[ - (8,'gray',0.5,False,'rgb(0.5,0.5,0.5)'), - (8,'gray',0.5,True,'rgba(0.5,0.5,0.5,1.0)'), - (8,'gray',[0.5,0.25],True,['rgba(0.5,0.5,0.5,1.0)','rgba(0.25,0.25,0.25,1.0)']), - ]) - def test_at_string(self, N, cmap, at, with_alpha, result): - assert np.all(Colormap.from_predefined(cmap,N=N).at(at,with_alpha=with_alpha,as_string=True) == result) - @pytest.mark.parametrize('bounds',[None,[2,10]]) def test_shade(self,ref_path,update,bounds): data = np.add(*np.indices((10, 11))) From 4235c70aea054bff3cce9531dee683a6b4847755 Mon Sep 17 00:00:00 2001 From: Daniel Otto de Mentock Date: Tue, 14 Dec 2021 17:01:40 +0100 Subject: [PATCH 25/74] added first typehints for rotation and orientation modules --- python/damask/_orientation.py | 136 ++++++++++-------- python/damask/_rotation.py | 251 ++++++++++++++++++---------------- 2 files changed, 211 insertions(+), 176 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index e727c54ae..6c0dd2fa3 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -1,5 +1,6 @@ import inspect import copy +from typing import Union, Callable, Sequence, Dict, Any, Tuple, List import numpy as np @@ -93,12 +94,12 @@ class Orientation(Rotation,Crystal): @util.extend_docstring(_parameter_doc) def __init__(self, - rotation = np.array([1.0,0.0,0.0,0.0]), *, - family = None, - lattice = None, - a = None,b = None,c = None, - alpha = None,beta = None,gamma = None, - degrees = False): + rotation: Union[Sequence[float], np.ndarray, Rotation] = np.array([1.,0.,0.,0.]), *, + family: str = None, + lattice: str = None, + a: float = None, b: float = None, c: float = None, + alpha: float = None, beta: float = None, gamma: float = None, + degrees: bool = False): """ New orientation. @@ -115,13 +116,12 @@ class Orientation(Rotation,Crystal): a=a,b=b,c=c, alpha=alpha,beta=beta,gamma=gamma, degrees=degrees) - def __repr__(self): + def __repr__(self) -> str: """Represent.""" return '\n'.join([Crystal.__repr__(self), Rotation.__repr__(self)]) - - def __copy__(self,rotation=None): + def __copy__(self,rotation: Union[Sequence[float], np.ndarray, Rotation] = None) -> "Orientation": """Create deep copy.""" dup = copy.deepcopy(self) if rotation is not None: @@ -131,7 +131,8 @@ class Orientation(Rotation,Crystal): copy = __copy__ - def __eq__(self,other): + + def __eq__(self, other: object) -> bool: """ Equal to other. @@ -141,12 +142,14 @@ class Orientation(Rotation,Crystal): Orientation to check for equality. """ + if not isinstance(other, Orientation): + raise TypeError matching_type = self.family == other.family and \ self.lattice == other.lattice and \ self.parameters == other.parameters return np.logical_and(matching_type,super(self.__class__,self.reduced).__eq__(other.reduced)) - def __ne__(self,other): + def __ne__(self, other: object) -> bool: """ Not equal to other. @@ -156,10 +159,16 @@ class Orientation(Rotation,Crystal): Orientation to check for equality. """ + if not isinstance(other, Orientation): + raise TypeError return np.logical_not(self==other) - def isclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True): + def isclose(self, + other: object, + rtol: float = 1e-5, + atol: float = 1e-8, + equal_nan: bool = True) -> bool: """ Report where values are approximately equal to corresponding ones of other Orientation. @@ -180,6 +189,8 @@ class Orientation(Rotation,Crystal): Mask indicating where corresponding orientations are close. """ + if not isinstance(other, Orientation): + raise TypeError matching_type = self.family == other.family and \ self.lattice == other.lattice and \ self.parameters == other.parameters @@ -187,7 +198,11 @@ class Orientation(Rotation,Crystal): - def allclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True): + def allclose(self, + other: object, + rtol: float = 1e-5, + atol: float = 1e-8, + equal_nan: bool = True) -> Union[bool, np.bool_]: """ Test whether all values are approximately equal to corresponding ones of other Orientation. @@ -208,10 +223,12 @@ class Orientation(Rotation,Crystal): Whether all values are close between both orientations. """ + if not isinstance(other, Orientation): + raise TypeError return np.all(self.isclose(other,rtol,atol,equal_nan)) - def __mul__(self,other): + def __mul__(self, other: Union[Rotation, "Orientation"]) -> "Orientation": """ Compose this orientation with other. @@ -233,7 +250,7 @@ class Orientation(Rotation,Crystal): @staticmethod - def _split_kwargs(kwargs,target): + def _split_kwargs(kwargs: Dict[str, Any], target: Callable) -> Tuple[Dict[str, Any], ...]: """ Separate keyword arguments in 'kwargs' targeted at 'target' from general keyword arguments of Orientation objects. @@ -252,7 +269,7 @@ class Orientation(Rotation,Crystal): Valid keyword arguments of Orientation object. """ - kws = () + kws: Tuple[Dict[str, Any], ...] = () for t in (target,Orientation.__init__): kws += ({key: kwargs[key] for key in set(inspect.signature(t).parameters) & set(kwargs)},) @@ -264,85 +281,88 @@ class Orientation(Rotation,Crystal): @classmethod - @util.extended_docstring(Rotation.from_random,_parameter_doc) - def from_random(cls,**kwargs): + @util.extended_docstring(Rotation.from_random, _parameter_doc) + def from_random(cls, **kwargs) -> "Orientation": kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_random) return cls(rotation=Rotation.from_random(**kwargs_rot),**kwargs_ori) @classmethod @util.extended_docstring(Rotation.from_quaternion,_parameter_doc) - def from_quaternion(cls,**kwargs): + def from_quaternion(cls, **kwargs) -> "Orientation": kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_quaternion) return cls(rotation=Rotation.from_quaternion(**kwargs_rot),**kwargs_ori) @classmethod @util.extended_docstring(Rotation.from_Euler_angles,_parameter_doc) - def from_Euler_angles(cls,**kwargs): + def from_Euler_angles(cls, **kwargs) -> "Orientation": kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_Euler_angles) return cls(rotation=Rotation.from_Euler_angles(**kwargs_rot),**kwargs_ori) @classmethod @util.extended_docstring(Rotation.from_axis_angle,_parameter_doc) - def from_axis_angle(cls,**kwargs): + def from_axis_angle(cls, **kwargs) -> "Orientation": kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_axis_angle) return cls(rotation=Rotation.from_axis_angle(**kwargs_rot),**kwargs_ori) @classmethod @util.extended_docstring(Rotation.from_basis,_parameter_doc) - def from_basis(cls,**kwargs): + def from_basis(cls, **kwargs) -> "Orientation": kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_basis) return cls(rotation=Rotation.from_basis(**kwargs_rot),**kwargs_ori) @classmethod @util.extended_docstring(Rotation.from_matrix,_parameter_doc) - def from_matrix(cls,**kwargs): + def from_matrix(cls, **kwargs) -> "Orientation": kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_matrix) return cls(rotation=Rotation.from_matrix(**kwargs_rot),**kwargs_ori) @classmethod @util.extended_docstring(Rotation.from_Rodrigues_vector,_parameter_doc) - def from_Rodrigues_vector(cls,**kwargs): + def from_Rodrigues_vector(cls, **kwargs) -> "Orientation": kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_Rodrigues_vector) return cls(rotation=Rotation.from_Rodrigues_vector(**kwargs_rot),**kwargs_ori) @classmethod @util.extended_docstring(Rotation.from_homochoric,_parameter_doc) - def from_homochoric(cls,**kwargs): + def from_homochoric(cls, **kwargs) -> "Orientation": kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_homochoric) return cls(rotation=Rotation.from_homochoric(**kwargs_rot),**kwargs_ori) @classmethod @util.extended_docstring(Rotation.from_cubochoric,_parameter_doc) - def from_cubochoric(cls,**kwargs): + def from_cubochoric(cls, **kwargs) -> "Orientation": kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_cubochoric) return cls(rotation=Rotation.from_cubochoric(**kwargs_rot),**kwargs_ori) @classmethod @util.extended_docstring(Rotation.from_spherical_component,_parameter_doc) - def from_spherical_component(cls,**kwargs): + def from_spherical_component(cls, **kwargs) -> "Orientation": kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_spherical_component) return cls(rotation=Rotation.from_spherical_component(**kwargs_rot),**kwargs_ori) @classmethod @util.extended_docstring(Rotation.from_fiber_component,_parameter_doc) - def from_fiber_component(cls,**kwargs): + def from_fiber_component(cls, **kwargs) -> "Orientation": kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_fiber_component) return cls(rotation=Rotation.from_fiber_component(**kwargs_rot),**kwargs_ori) @classmethod @util.extend_docstring(_parameter_doc) - def from_directions(cls,uvw,hkl,**kwargs): + def from_directions(cls, + uvw: Union[Sequence[float], np.ndarray], + hkl: Union[Sequence[float], np.ndarray], + **kwargs) -> "Orientation": """ Initialize orientation object from two crystallographic directions. @@ -362,7 +382,7 @@ class Orientation(Rotation,Crystal): @property - def equivalent(self): + def equivalent(self) -> "Orientation": """ Orientations that are symmetrically equivalent. @@ -376,7 +396,7 @@ class Orientation(Rotation,Crystal): @property - def reduced(self): + def reduced(self) -> "Orientation": """Select symmetrically equivalent orientation that falls into fundamental zone according to symmetry.""" eq = self.equivalent ok = eq.in_FZ @@ -385,9 +405,8 @@ class Orientation(Rotation,Crystal): sort = 0 if len(loc) == 1 else np.lexsort(loc[:0:-1]) return eq[ok][sort].reshape(self.shape) - @property - def in_FZ(self): + def in_FZ(self) -> Union[np.bool_, np.ndarray]: """ Check whether orientation falls into fundamental zone of own symmetry. @@ -431,7 +450,7 @@ class Orientation(Rotation,Crystal): @property - def in_disorientation_FZ(self): + def in_disorientation_FZ(self) -> np.ndarray: """ Check whether orientation falls into fundamental zone of disorientations. @@ -471,8 +490,7 @@ class Orientation(Rotation,Crystal): else: return np.ones_like(rho[...,0],dtype=bool) - - def disorientation(self,other,return_operators=False): + def disorientation(self, other, return_operators = False): """ Calculate disorientation between myself and given other orientation. @@ -518,8 +536,8 @@ class Orientation(Rotation,Crystal): if self.family != other.family: raise NotImplementedError('disorientation between different crystal families') - blend = util.shapeblender(self.shape,other.shape) - s = self.equivalent + blend: Tuple[int, ...] = util.shapeblender(self.shape,other.shape) + s = self.equivalent o = other.equivalent s_ = s.reshape((s.shape[0],1)+ self.shape).broadcast_to((s.shape[0],o.shape[0])+blend,mode='right') @@ -534,10 +552,9 @@ class Orientation(Rotation,Crystal): r = np.where(np.any(forward[...,np.newaxis],axis=(0,1),keepdims=True), r_.quaternion, _r.quaternion) - loc = np.where(ok) - sort = 0 if len(loc) == 2 else np.lexsort(loc[:1:-1]) - quat = r[ok][sort].reshape(blend+(4,)) - + loc: Tuple[float] = np.where(ok) + sort: np.ndarray = 0 if len(loc) == 2 else np.lexsort(loc[:1:-1]) + quat: np.ndarray = r[ok][sort].reshape(blend+(4,)) return ( (self.copy(rotation=quat), (np.vstack(loc[:2]).T)[sort].reshape(blend+(2,))) @@ -546,7 +563,7 @@ class Orientation(Rotation,Crystal): ) - def average(self,weights=None,return_cloud=False): + def average(self, weights = None, return_cloud = False): """ Return orientation average over last dimension. @@ -587,7 +604,10 @@ class Orientation(Rotation,Crystal): ) - def to_SST(self,vector,proper=False,return_operators=False): + def to_SST(self, + vector: np.ndarray, + proper: bool = False, + return_operators: bool = False) -> np.ndarray: """ Rotate vector to ensure it falls into (improper or proper) standard stereographic triangle of crystal symmetry. @@ -626,7 +646,7 @@ class Orientation(Rotation,Crystal): ) - def in_SST(self,vector,proper=False): + def in_SST(self, vector: np.ndarray, proper: bool = False) -> Union[np.bool_, np.ndarray]: """ Check whether given crystal frame vector falls into standard stereographic triangle of own symmetry. @@ -667,7 +687,7 @@ class Orientation(Rotation,Crystal): return np.all(components >= 0.0,axis=-1) - def IPF_color(self,vector,in_SST=True,proper=False): + def IPF_color(self, vector: np.ndarray, in_SST: bool = True, proper: bool = False) -> np.ndarray: """ Map vector to RGB color within standard stereographic triangle of own symmetry. @@ -715,30 +735,30 @@ class Orientation(Rotation,Crystal): components_improper = np.around(np.einsum('...ji,...i', np.broadcast_to(self.standard_triangle['improper'], vector_.shape+(3,)), vector_), 12) - in_SST = np.all(components_proper >= 0.0,axis=-1) \ + in_SST_ = np.all(components_proper >= 0.0,axis=-1) \ | np.all(components_improper >= 0.0,axis=-1) - components = np.where((in_SST & np.all(components_proper >= 0.0,axis=-1))[...,np.newaxis], + components = np.where((in_SST_ & np.all(components_proper >= 0.0,axis=-1))[...,np.newaxis], components_proper,components_improper) else: components = np.around(np.einsum('...ji,...i', np.broadcast_to(self .standard_triangle['improper'], vector_.shape+(3,)), np.block([vector_[...,:2],np.abs(vector_[...,2:3])])), 12) - in_SST = np.all(components >= 0.0,axis=-1) + in_SST_ = np.all(components >= 0.0,axis=-1) with np.errstate(invalid='ignore',divide='ignore'): rgb = (components/np.linalg.norm(components,axis=-1,keepdims=True))**0.5 # smoothen color ramps rgb = np.clip(rgb,0.,1.) # clip intensity rgb /= np.max(rgb,axis=-1,keepdims=True) # normalize to (HS)V = 1 - rgb[np.broadcast_to(~in_SST[...,np.newaxis],rgb.shape)] = 0.0 + rgb[np.broadcast_to(~in_SST_[...,np.newaxis],rgb.shape)] = 0.0 return rgb @property - def symmetry_operations(self): + def symmetry_operations(self) -> Rotation: """Symmetry operations as Rotations.""" - _symmetry_operations = { + _symmetry_operations: Dict[str, List[List]] = { 'cubic': [ [ 1.0, 0.0, 0.0, 0.0 ], [ 0.0, 1.0, 0.0, 0.0 ], @@ -808,7 +828,10 @@ class Orientation(Rotation,Crystal): #################################################################################################### # functions that require lattice, not just family - def to_pole(self,*,uvw=None,hkl=None,with_symmetry=False): + def to_pole(self, *, + uvw: np.ndarray = None, + hkl: np.ndarray = None, + with_symmetry: bool = False) -> np.ndarray: """ Calculate lab frame vector along lattice direction [uvw] or plane normal (hkl). @@ -839,7 +862,9 @@ class Orientation(Rotation,Crystal): @ np.broadcast_to(v,blend+(3,)) - def Schmid(self,*,N_slip=None,N_twin=None): + def Schmid(self, *, + N_slip: Sequence[int] = None, + N_twin: Sequence[int] = None) -> np.ndarray: u""" Calculate Schmid matrix P = d ⨂ n in the lab frame for selected deformation systems. @@ -876,6 +901,7 @@ class Orientation(Rotation,Crystal): (self.kinematics('twin'),N_twin) if active == '*': active = [len(a) for a in kinematics['direction']] + assert active d = self.to_frame(uvw=np.vstack([kinematics['direction'][i][:n] for i,n in enumerate(active)])) p = self.to_frame(hkl=np.vstack([kinematics['plane'][i][:n] for i,n in enumerate(active)])) P = np.einsum('...i,...j',d/np.linalg.norm(d,axis=1,keepdims=True), @@ -886,7 +912,7 @@ class Orientation(Rotation,Crystal): @ np.broadcast_to(P.reshape(util.shapeshifter(P.shape,shape)),shape) - def related(self,model): + def related(self, model: str) -> "Orientation": """ Orientations derived from the given relationship. diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index ac921d70a..0a4421090 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -6,6 +6,8 @@ from . import tensor from . import util from . import grid_filters +from typing import Union, Sequence, Tuple, Literal, List + _P = -1 # parameters for conversion from/to cubochoric @@ -61,7 +63,7 @@ class Rotation: __slots__ = ['quaternion'] - def __init__(self,rotation = np.array([1.0,0.0,0.0,0.0])): + def __init__(self, rotation: Union[Sequence[float], np.ndarray, "Rotation"] = np.array([1.0,0.0,0.0,0.0])): """ New rotation. @@ -73,6 +75,7 @@ class Rotation: Defaults to no rotation. """ + self.quaternion: np.ndarray if isinstance(rotation,Rotation): self.quaternion = rotation.quaternion.copy() elif np.array(rotation).shape[-1] == 4: @@ -81,13 +84,13 @@ class Rotation: raise TypeError('"rotation" is neither a Rotation nor a quaternion') - def __repr__(self): + def __repr__(self) -> str: """Represent rotation as unit quaternion(s).""" return f'Quaternion{" " if self.quaternion.shape == (4,) else "s of shape "+str(self.quaternion.shape[:-1])+chr(10)}'\ + str(self.quaternion) - def __copy__(self,rotation=None): + def __copy__(self, rotation: Union[Sequence[float], np.ndarray, "Rotation"] = None) -> "Rotation": """Create deep copy.""" dup = copy.deepcopy(self) if rotation is not None: @@ -97,14 +100,14 @@ class Rotation: copy = __copy__ - def __getitem__(self,item): + def __getitem__(self, item: Union[Tuple[int], int, bool, np.bool_, np.ndarray]): """Return slice according to item.""" return self.copy() \ if self.shape == () else \ self.copy(rotation=self.quaternion[item+(slice(None),)] if isinstance(item,tuple) else self.quaternion[item]) - def __eq__(self,other): + def __eq__(self, other: object) -> bool: """ Equal to other. @@ -114,11 +117,13 @@ class Rotation: Rotation to check for equality. """ + if not isinstance(other, Rotation): + raise TypeError return np.logical_or(np.all(self.quaternion == other.quaternion,axis=-1), np.all(self.quaternion == -1.0*other.quaternion,axis=-1)) - def __ne__(self,other): + def __ne__(self, other: object) -> bool: """ Not equal to other. @@ -128,10 +133,12 @@ class Rotation: Rotation to check for inequality. """ + if not isinstance(other, Rotation): + raise TypeError return np.logical_not(self==other) - def isclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True): + def isclose(self, other: "Rotation", rtol: float = 1e-5, atol: float = 1e-8, equal_nan: bool = True) -> bool: """ Report where values are approximately equal to corresponding ones of other Rotation. @@ -158,7 +165,7 @@ class Rotation: np.all(np.isclose(s,-1.0*o,rtol,atol,equal_nan),axis=-1)) - def allclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True): + def allclose(self, other: "Rotation", rtol: float = 1e-5, atol: float = 1e-8, equal_nan: bool = True) -> Union[np.bool_, bool]: """ Test whether all values are approximately equal to corresponding ones of other Rotation. @@ -188,27 +195,27 @@ class Rotation: @property - def size(self): + def size(self) -> int: return self.quaternion[...,0].size @property - def shape(self): + def shape(self) -> Tuple[int, ...]: return self.quaternion[...,0].shape - def __len__(self): + def __len__(self) -> int: """Length of leading/leftmost dimension of array.""" return 0 if self.shape == () else self.shape[0] - def __invert__(self): + def __invert__(self) -> "Rotation": """Inverse rotation (backward rotation).""" dup = self.copy() dup.quaternion[...,1:] *= -1 return dup - def __pow__(self,exp): + def __pow__(self, exp: int) -> "Rotation": """ Perform the rotation 'exp' times. @@ -222,7 +229,7 @@ class Rotation: p = self.quaternion[...,1:]/np.linalg.norm(self.quaternion[...,1:],axis=-1,keepdims=True) return self.copy(rotation=Rotation(np.block([np.cos(exp*phi),np.sin(exp*phi)*p]))._standardize()) - def __ipow__(self,exp): + def __ipow__(self, exp: int) -> "Rotation": """ Perform the rotation 'exp' times (in-place). @@ -235,7 +242,7 @@ class Rotation: return self**exp - def __mul__(self,other): + def __mul__(self, other: "Rotation") -> "Rotation": """ Compose with other. @@ -261,7 +268,7 @@ class Rotation: else: raise TypeError('Use "R@b", i.e. matmul, to apply rotation "R" to object "b"') - def __imul__(self,other): + def __imul__(self, other: "Rotation") -> "Rotation": """ Compose with other (in-place). @@ -274,7 +281,7 @@ class Rotation: return self*other - def __truediv__(self,other): + def __truediv__(self, other: "Rotation") -> "Rotation": """ Compose with inverse of other. @@ -294,7 +301,7 @@ class Rotation: else: raise TypeError('Use "R@b", i.e. matmul, to apply rotation "R" to object "b"') - def __itruediv__(self,other): + def __itruediv__(self, other: "Rotation") -> "Rotation": """ Compose with inverse of other (in-place). @@ -307,7 +314,7 @@ class Rotation: return self/other - def __matmul__(self,other): + def __matmul__(self, other: np.ndarray) -> np.ndarray: """ Rotate vector, second order tensor, or fourth order tensor. @@ -350,13 +357,13 @@ class Rotation: apply = __matmul__ - def _standardize(self): + def _standardize(self) -> "Rotation": """Standardize quaternion (ensure positive real hemisphere).""" self.quaternion[self.quaternion[...,0] < 0.0] *= -1 return self - def append(self,other): + def append(self, other: Union["Rotation", List["Rotation"]]) -> "Rotation": """ Extend array along first dimension with other array(s). @@ -369,7 +376,7 @@ class Rotation: [self]+other if isinstance(other,list) else [self,other])))) - def flatten(self,order = 'C'): + def flatten(self, order: Literal['C','F','A'] = 'C') -> "Rotation": """ Flatten array. @@ -382,7 +389,7 @@ class Rotation: return self.copy(rotation=self.quaternion.reshape((-1,4),order=order)) - def reshape(self,shape,order = 'C'): + def reshape(self,shape: Union[int, Tuple[int, ...]], order: Literal['C','F','A'] = 'C') -> "Rotation": """ Reshape array. @@ -396,7 +403,7 @@ class Rotation: return self.copy(rotation=self.quaternion.reshape(tuple(shape)+(4,),order=order)) - def broadcast_to(self,shape,mode = 'right'): + def broadcast_to(self, shape: Union[int, Tuple[int, ...]], mode: Literal["left", "right"] = 'right') -> "Rotation": """ Broadcast array. @@ -458,7 +465,7 @@ class Rotation: accept_homomorph = True) - def misorientation(self,other): + def misorientation(self, other: "Rotation") -> "Rotation": """ Calculate misorientation to other Rotation. @@ -479,7 +486,7 @@ class Rotation: ################################################################################################ # convert to different orientation representations (numpy arrays) - def as_quaternion(self): + def as_quaternion(self) -> np.ndarray: """ Represent as unit quaternion. @@ -492,7 +499,7 @@ class Rotation: return self.quaternion.copy() def as_Euler_angles(self, - degrees = False): + degrees: bool = False) -> np.ndarray: """ Represent as Bunge Euler angles. @@ -526,8 +533,8 @@ class Rotation: return eu def as_axis_angle(self, - degrees = False, - pair = False): + degrees: bool = False, + pair: bool = False) -> Union[Tuple[float, ...], np.ndarray]: """ Represent as axis–angle pair. @@ -554,11 +561,11 @@ class Rotation: (array([0., 0., 1.]), array(0.)) """ - ax = Rotation._qu2ax(self.quaternion) + ax: np.ndarray = Rotation._qu2ax(self.quaternion) if degrees: ax[...,3] = np.degrees(ax[...,3]) return (ax[...,:3],ax[...,3]) if pair else ax - def as_matrix(self): + def as_matrix(self) -> np.ndarray: """ Represent as rotation matrix. @@ -582,7 +589,7 @@ class Rotation: return Rotation._qu2om(self.quaternion) def as_Rodrigues_vector(self, - compact = False): + compact: bool = False) -> np.ndarray: """ Represent as Rodrigues–Frank vector with separate axis and angle argument. @@ -615,7 +622,7 @@ class Rotation: else: return ro - def as_homochoric(self): + def as_homochoric(self) -> np.ndarray: """ Represent as homochoric vector. @@ -636,7 +643,7 @@ class Rotation: """ return Rotation._qu2ho(self.quaternion) - def as_cubochoric(self): + def as_cubochoric(self) -> np.ndarray: """ Represent as cubochoric vector. @@ -661,9 +668,9 @@ class Rotation: # Static constructors. The input data needs to follow the conventions, options allow to # relax the conventions. @staticmethod - def from_quaternion(q, - accept_homomorph = False, - P = -1): + def from_quaternion(q: Union[Sequence[Sequence[float]], np.ndarray], + accept_homomorph: bool = False, + P: Literal[1, -1] = -1) -> "Rotation": """ Initialize from quaternion. @@ -678,7 +685,7 @@ class Rotation: Sign convention. Defaults to -1. """ - qu = np.array(q,dtype=float) + qu: np.ndarray = np.array(q,dtype=float) if qu.shape[:-2:-1] != (4,): raise ValueError('Invalid shape.') if abs(P) != 1: @@ -696,8 +703,8 @@ class Rotation: return Rotation(qu) @staticmethod - def from_Euler_angles(phi, - degrees = False): + def from_Euler_angles(phi: np.ndarray, + degrees: bool = False) -> "Rotation": """ Initialize from Bunge Euler angles. @@ -725,10 +732,10 @@ class Rotation: return Rotation(Rotation._eu2qu(eu)) @staticmethod - def from_axis_angle(axis_angle, - degrees = False, - normalize = False, - P = -1): + def from_axis_angle(axis_angle: np.ndarray, + degrees:bool = False, + normalize: bool = False, + P: Literal[1, -1] = -1) -> "Rotation": """ Initialize from Axis angle pair. @@ -763,9 +770,9 @@ class Rotation: return Rotation(Rotation._ax2qu(ax)) @staticmethod - def from_basis(basis, - orthonormal = True, - reciprocal = False): + def from_basis(basis: np.ndarray, + orthonormal: bool = True, + reciprocal: bool = False) -> "Rotation": """ Initialize from lattice basis vectors. @@ -799,7 +806,7 @@ class Rotation: return Rotation(Rotation._om2qu(om)) @staticmethod - def from_matrix(R): + def from_matrix(R: np.ndarray) -> "Rotation": """ Initialize from rotation matrix. @@ -812,8 +819,9 @@ class Rotation: return Rotation.from_basis(R) @staticmethod - def from_parallel(a,b, - **kwargs): + def from_parallel(a: np.ndarray, + b: np.ndarray, + **kwargs) -> "Rotation": """ Initialize from pairs of two orthogonal lattice basis vectors. @@ -841,9 +849,9 @@ class Rotation: @staticmethod - def from_Rodrigues_vector(rho, - normalize = False, - P = -1): + def from_Rodrigues_vector(rho: np.ndarray, + normalize: bool = False, + P: Literal[1, -1] = -1) -> "Rotation": """ Initialize from Rodrigues–Frank vector (angle separated from axis). @@ -873,8 +881,8 @@ class Rotation: return Rotation(Rotation._ro2qu(ro)) @staticmethod - def from_homochoric(h, - P = -1): + def from_homochoric(h: np.ndarray, + P: Literal[1, -1] = -1) -> "Rotation": """ Initialize from homochoric vector. @@ -900,8 +908,8 @@ class Rotation: return Rotation(Rotation._ho2qu(ho)) @staticmethod - def from_cubochoric(x, - P = -1): + def from_cubochoric(x: np.ndarray, + P: Literal[1, -1] = -1) -> "Rotation": """ Initialize from cubochoric vector. @@ -928,8 +936,8 @@ class Rotation: @staticmethod - def from_random(shape = None, - rng_seed = None): + def from_random(shape: Tuple[int, ...] = None, + rng_seed: Union[None, int, Sequence[int], np.ndarray] = None) -> "Rotation": """ Initialize with random rotation. @@ -958,13 +966,13 @@ class Rotation: @staticmethod - def from_ODF(weights, - phi, - N = 500, - degrees = True, - fractions = True, - rng_seed = None, - **kwargs): + def from_ODF(weights: np.ndarray, + phi: np.ndarray, + N: int = 500, + degrees: bool = True, + fractions: bool = True, + rng_seed: Union[None, int, Sequence[int], np.ndarray] = None, + **kwargs) -> "Rotation": """ Sample discrete values from a binned orientation distribution function (ODF). @@ -1017,11 +1025,11 @@ class Rotation: @staticmethod - def from_spherical_component(center, - sigma, - N = 500, - degrees = True, - rng_seed = None): + def from_spherical_component(center: "Rotation", + sigma: float, + N: int = 500, + degrees: bool = True, + rng_seed: Union[None, int, Sequence[float], np.ndarray] = None) -> "Rotation": """ Calculate set of rotations with Gaussian distribution around center. @@ -1052,12 +1060,12 @@ class Rotation: @staticmethod - def from_fiber_component(alpha, - beta, - sigma = 0.0, - N = 500, - degrees = True, - rng_seed = None): + def from_fiber_component(alpha: np.ndarray, + beta: np.ndarray, + sigma: float = 0.0, + N: int = 500, + degrees: bool = True, + rng_seed: bool = None): """ Calculate set of rotations with Gaussian distribution around direction. @@ -1080,7 +1088,8 @@ class Rotation: """ rng = np.random.default_rng(rng_seed) - sigma_,alpha_,beta_ = map(np.radians,(sigma,alpha,beta)) if degrees else (sigma,alpha,beta) + sigma_: np.ndarray; alpha_: np.ndarray; beta_: np.ndarray + sigma_,alpha_,beta_ = map(np.radians,(sigma,alpha,beta)) if degrees else (sigma,alpha,beta) # type: ignore d_cr = np.array([np.sin(alpha_[0])*np.cos(alpha_[1]), np.sin(alpha_[0])*np.sin(alpha_[1]), np.cos(alpha_[0])]) d_lab = np.array([np.sin( beta_[0])*np.cos( beta_[1]), np.sin( beta_[0])*np.sin( beta_[1]), np.cos( beta_[0])]) @@ -1134,7 +1143,7 @@ class Rotation: #################################################################################################### #---------- Quaternion ---------- @staticmethod - def _qu2om(qu): + def _qu2om(qu: np.ndarray) -> np.ndarray: qq = qu[...,0:1]**2-(qu[...,1:2]**2 + qu[...,2:3]**2 + qu[...,3:4]**2) om = np.block([qq + 2.0*qu[...,1:2]**2, 2.0*(qu[...,2:3]*qu[...,1:2]-_P*qu[...,0:1]*qu[...,3:4]), @@ -1149,7 +1158,7 @@ class Rotation: return om @staticmethod - def _qu2eu(qu): + def _qu2eu(qu: np.ndarray) -> np.ndarray: """Quaternion to Bunge Euler angles.""" q02 = qu[...,0:1]*qu[...,2:3] q13 = qu[...,1:2]*qu[...,3:4] @@ -1177,7 +1186,7 @@ class Rotation: return eu @staticmethod - def _qu2ax(qu): + def _qu2ax(qu: np.ndarray) -> np.ndarray: """ Quaternion to axis–angle pair. @@ -1193,7 +1202,7 @@ class Rotation: return ax @staticmethod - def _qu2ro(qu): + def _qu2ro(qu: np.ndarray) -> np.ndarray: """Quaternion to Rodrigues–Frank vector.""" with np.errstate(invalid='ignore',divide='ignore'): s = np.linalg.norm(qu[...,1:4],axis=-1,keepdims=True) @@ -1207,7 +1216,7 @@ class Rotation: return ro @staticmethod - def _qu2ho(qu): + def _qu2ho(qu: np.ndarray) -> np.ndarray: """Quaternion to homochoric vector.""" with np.errstate(invalid='ignore'): omega = 2.0 * np.arccos(np.clip(qu[...,0:1],-1.0,1.0)) @@ -1218,14 +1227,14 @@ class Rotation: return ho @staticmethod - def _qu2cu(qu): + def _qu2cu(qu: np.ndarray) -> np.ndarray: """Quaternion to cubochoric vector.""" return Rotation._ho2cu(Rotation._qu2ho(qu)) #---------- Rotation matrix ---------- @staticmethod - def _om2qu(om): + def _om2qu(om: np.ndarray) -> np.ndarray: """ Rotation matrix to quaternion. @@ -1267,7 +1276,7 @@ class Rotation: return qu @staticmethod - def _om2eu(om): + def _om2eu(om: np.ndarray) -> np.ndarray: """Rotation matrix to Bunge Euler angles.""" with np.errstate(invalid='ignore',divide='ignore'): zeta = 1.0/np.sqrt(1.0-om[...,2,2:3]**2) @@ -1286,7 +1295,7 @@ class Rotation: return eu @staticmethod - def _om2ax(om): + def _om2ax(om: np.ndarray) -> np.ndarray: """Rotation matrix to axis–angle pair.""" diag_delta = -_P*np.block([om[...,1,2:3]-om[...,2,1:2], om[...,2,0:1]-om[...,0,2:3], @@ -1307,24 +1316,24 @@ class Rotation: return ax @staticmethod - def _om2ro(om): + def _om2ro(om: np.ndarray) -> np.ndarray: """Rotation matrix to Rodrigues–Frank vector.""" return Rotation._eu2ro(Rotation._om2eu(om)) @staticmethod - def _om2ho(om): + def _om2ho(om: np.ndarray) -> np.ndarray: """Rotation matrix to homochoric vector.""" return Rotation._ax2ho(Rotation._om2ax(om)) @staticmethod - def _om2cu(om): + def _om2cu(om: np.ndarray) -> np.ndarray: """Rotation matrix to cubochoric vector.""" return Rotation._ho2cu(Rotation._om2ho(om)) #---------- Bunge Euler angles ---------- @staticmethod - def _eu2qu(eu): + def _eu2qu(eu: np.ndarray) -> np.ndarray: """Bunge Euler angles to quaternion.""" ee = 0.5*eu cPhi = np.cos(ee[...,1:2]) @@ -1337,7 +1346,7 @@ class Rotation: return qu @staticmethod - def _eu2om(eu): + def _eu2om(eu: np.ndarray) -> np.ndarray: """Bunge Euler angles to rotation matrix.""" c = np.cos(eu) s = np.sin(eu) @@ -1355,7 +1364,7 @@ class Rotation: return om @staticmethod - def _eu2ax(eu): + def _eu2ax(eu: np.ndarray) -> np.ndarray: """Bunge Euler angles to axis–angle pair.""" t = np.tan(eu[...,1:2]*0.5) sigma = 0.5*(eu[...,0:1]+eu[...,2:3]) @@ -1374,7 +1383,7 @@ class Rotation: return ax @staticmethod - def _eu2ro(eu): + def _eu2ro(eu: np.ndarray) -> np.ndarray: """Bunge Euler angles to Rodrigues–Frank vector.""" ax = Rotation._eu2ax(eu) ro = np.block([ax[...,:3],np.tan(ax[...,3:4]*.5)]) @@ -1383,19 +1392,19 @@ class Rotation: return ro @staticmethod - def _eu2ho(eu): + def _eu2ho(eu: np.ndarray) -> np.ndarray: """Bunge Euler angles to homochoric vector.""" return Rotation._ax2ho(Rotation._eu2ax(eu)) @staticmethod - def _eu2cu(eu): + def _eu2cu(eu: np.ndarray) -> np.ndarray: """Bunge Euler angles to cubochoric vector.""" return Rotation._ho2cu(Rotation._eu2ho(eu)) #---------- Axis angle pair ---------- @staticmethod - def _ax2qu(ax): + def _ax2qu(ax: np.ndarray) -> np.ndarray: """Axis–angle pair to quaternion.""" c = np.cos(ax[...,3:4]*.5) s = np.sin(ax[...,3:4]*.5) @@ -1403,7 +1412,7 @@ class Rotation: return qu @staticmethod - def _ax2om(ax): + def _ax2om(ax: np.ndarray) -> np.ndarray: """Axis-angle pair to rotation matrix.""" c = np.cos(ax[...,3:4]) s = np.sin(ax[...,3:4]) @@ -1420,12 +1429,12 @@ class Rotation: return om if _P < 0.0 else np.swapaxes(om,-1,-2) @staticmethod - def _ax2eu(ax): + def _ax2eu(ax: np.ndarray) -> np.ndarray: """Rotation matrix to Bunge Euler angles.""" return Rotation._om2eu(Rotation._ax2om(ax)) @staticmethod - def _ax2ro(ax): + def _ax2ro(ax: np.ndarray) -> np.ndarray: """Axis–angle pair to Rodrigues–Frank vector.""" ro = np.block([ax[...,:3], np.where(np.isclose(ax[...,3:4],np.pi,atol=1.e-15,rtol=.0), @@ -1436,36 +1445,36 @@ class Rotation: return ro @staticmethod - def _ax2ho(ax): + def _ax2ho(ax: np.ndarray) -> np.ndarray: """Axis–angle pair to homochoric vector.""" f = (0.75 * ( ax[...,3:4] - np.sin(ax[...,3:4]) ))**(1.0/3.0) ho = ax[...,:3] * f return ho @staticmethod - def _ax2cu(ax): + def _ax2cu(ax: np.ndarray) -> np.ndarray: """Axis–angle pair to cubochoric vector.""" return Rotation._ho2cu(Rotation._ax2ho(ax)) #---------- Rodrigues-Frank vector ---------- @staticmethod - def _ro2qu(ro): + def _ro2qu(ro: np.ndarray) -> np.ndarray: """Rodrigues–Frank vector to quaternion.""" return Rotation._ax2qu(Rotation._ro2ax(ro)) @staticmethod - def _ro2om(ro): + def _ro2om(ro: np.ndarray) -> np.ndarray: """Rodgrigues–Frank vector to rotation matrix.""" return Rotation._ax2om(Rotation._ro2ax(ro)) @staticmethod - def _ro2eu(ro): + def _ro2eu(ro: np.ndarray) -> np.ndarray: """Rodrigues–Frank vector to Bunge Euler angles.""" return Rotation._om2eu(Rotation._ro2om(ro)) @staticmethod - def _ro2ax(ro): + def _ro2ax(ro: np.ndarray) -> np.ndarray: """Rodrigues–Frank vector to axis–angle pair.""" with np.errstate(invalid='ignore',divide='ignore'): ax = np.where(np.isfinite(ro[...,3:4]), @@ -1475,7 +1484,7 @@ class Rotation: return ax @staticmethod - def _ro2ho(ro): + def _ro2ho(ro: np.ndarray) -> np.ndarray: """Rodrigues–Frank vector to homochoric vector.""" f = np.where(np.isfinite(ro[...,3:4]),2.0*np.arctan(ro[...,3:4]) -np.sin(2.0*np.arctan(ro[...,3:4])),np.pi) ho = np.where(np.broadcast_to(np.sum(ro[...,0:3]**2.0,axis=-1,keepdims=True) < 1.e-8,ro[...,0:3].shape), @@ -1483,29 +1492,29 @@ class Rotation: return ho @staticmethod - def _ro2cu(ro): + def _ro2cu(ro: np.ndarray) -> np.ndarray: """Rodrigues–Frank vector to cubochoric vector.""" return Rotation._ho2cu(Rotation._ro2ho(ro)) #---------- Homochoric vector---------- @staticmethod - def _ho2qu(ho): + def _ho2qu(ho: np.ndarray) -> np.ndarray: """Homochoric vector to quaternion.""" return Rotation._ax2qu(Rotation._ho2ax(ho)) @staticmethod - def _ho2om(ho): + def _ho2om(ho: np.ndarray) -> np.ndarray: """Homochoric vector to rotation matrix.""" return Rotation._ax2om(Rotation._ho2ax(ho)) @staticmethod - def _ho2eu(ho): + def _ho2eu(ho: np.ndarray) -> np.ndarray: """Homochoric vector to Bunge Euler angles.""" return Rotation._ax2eu(Rotation._ho2ax(ho)) @staticmethod - def _ho2ax(ho): + def _ho2ax(ho: np.ndarray) -> np.ndarray: """Homochoric vector to axis–angle pair.""" tfit = np.array([+1.0000000000018852, -0.5000000002194847, -0.024999992127593126, -0.003928701544781374, @@ -1528,12 +1537,12 @@ class Rotation: return ax @staticmethod - def _ho2ro(ho): + def _ho2ro(ho: np.ndarray) -> np.ndarray: """Axis–angle pair to Rodrigues–Frank vector.""" return Rotation._ax2ro(Rotation._ho2ax(ho)) @staticmethod - def _ho2cu(ho): + def _ho2cu(ho: np.ndarray) -> np.ndarray: """ Homochoric vector to cubochoric vector. @@ -1573,32 +1582,32 @@ class Rotation: #---------- Cubochoric ---------- @staticmethod - def _cu2qu(cu): + def _cu2qu(cu: np.ndarray) -> np.ndarray: """Cubochoric vector to quaternion.""" return Rotation._ho2qu(Rotation._cu2ho(cu)) @staticmethod - def _cu2om(cu): + def _cu2om(cu: np.ndarray) -> np.ndarray: """Cubochoric vector to rotation matrix.""" return Rotation._ho2om(Rotation._cu2ho(cu)) @staticmethod - def _cu2eu(cu): + def _cu2eu(cu: np.ndarray) -> np.ndarray: """Cubochoric vector to Bunge Euler angles.""" return Rotation._ho2eu(Rotation._cu2ho(cu)) @staticmethod - def _cu2ax(cu): + def _cu2ax(cu: np.ndarray) -> np.ndarray: """Cubochoric vector to axis–angle pair.""" return Rotation._ho2ax(Rotation._cu2ho(cu)) @staticmethod - def _cu2ro(cu): + def _cu2ro(cu: np.ndarray) -> np.ndarray: """Cubochoric vector to Rodrigues–Frank vector.""" return Rotation._ho2ro(Rotation._cu2ho(cu)) @staticmethod - def _cu2ho(cu): + def _cu2ho(cu: np.ndarray) -> np.ndarray: """ Cubochoric vector to homochoric vector. @@ -1641,7 +1650,7 @@ class Rotation: @staticmethod - def _get_pyramid_order(xyz,direction=None): + def _get_pyramid_order(xyz: np.ndarray, direction: Literal['forward', 'backward']) -> np.ndarray: """ Get order of the coordinates. From 3b9822ad6750cb1ffc37e51db2144d20c09aff9d Mon Sep 17 00:00:00 2001 From: Daniel Otto de Mentock Date: Tue, 14 Dec 2021 16:04:18 +0000 Subject: [PATCH 26/74] Revert "added first typehints for rotation and orientation modules" This reverts commit 4235c70aea054bff3cce9531dee683a6b4847755 --- python/damask/_orientation.py | 136 ++++++++---------- python/damask/_rotation.py | 251 ++++++++++++++++------------------ 2 files changed, 176 insertions(+), 211 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 6c0dd2fa3..e727c54ae 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -1,6 +1,5 @@ import inspect import copy -from typing import Union, Callable, Sequence, Dict, Any, Tuple, List import numpy as np @@ -94,12 +93,12 @@ class Orientation(Rotation,Crystal): @util.extend_docstring(_parameter_doc) def __init__(self, - rotation: Union[Sequence[float], np.ndarray, Rotation] = np.array([1.,0.,0.,0.]), *, - family: str = None, - lattice: str = None, - a: float = None, b: float = None, c: float = None, - alpha: float = None, beta: float = None, gamma: float = None, - degrees: bool = False): + rotation = np.array([1.0,0.0,0.0,0.0]), *, + family = None, + lattice = None, + a = None,b = None,c = None, + alpha = None,beta = None,gamma = None, + degrees = False): """ New orientation. @@ -116,12 +115,13 @@ class Orientation(Rotation,Crystal): a=a,b=b,c=c, alpha=alpha,beta=beta,gamma=gamma, degrees=degrees) - def __repr__(self) -> str: + def __repr__(self): """Represent.""" return '\n'.join([Crystal.__repr__(self), Rotation.__repr__(self)]) - def __copy__(self,rotation: Union[Sequence[float], np.ndarray, Rotation] = None) -> "Orientation": + + def __copy__(self,rotation=None): """Create deep copy.""" dup = copy.deepcopy(self) if rotation is not None: @@ -131,8 +131,7 @@ class Orientation(Rotation,Crystal): copy = __copy__ - - def __eq__(self, other: object) -> bool: + def __eq__(self,other): """ Equal to other. @@ -142,14 +141,12 @@ class Orientation(Rotation,Crystal): Orientation to check for equality. """ - if not isinstance(other, Orientation): - raise TypeError matching_type = self.family == other.family and \ self.lattice == other.lattice and \ self.parameters == other.parameters return np.logical_and(matching_type,super(self.__class__,self.reduced).__eq__(other.reduced)) - def __ne__(self, other: object) -> bool: + def __ne__(self,other): """ Not equal to other. @@ -159,16 +156,10 @@ class Orientation(Rotation,Crystal): Orientation to check for equality. """ - if not isinstance(other, Orientation): - raise TypeError return np.logical_not(self==other) - def isclose(self, - other: object, - rtol: float = 1e-5, - atol: float = 1e-8, - equal_nan: bool = True) -> bool: + def isclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True): """ Report where values are approximately equal to corresponding ones of other Orientation. @@ -189,8 +180,6 @@ class Orientation(Rotation,Crystal): Mask indicating where corresponding orientations are close. """ - if not isinstance(other, Orientation): - raise TypeError matching_type = self.family == other.family and \ self.lattice == other.lattice and \ self.parameters == other.parameters @@ -198,11 +187,7 @@ class Orientation(Rotation,Crystal): - def allclose(self, - other: object, - rtol: float = 1e-5, - atol: float = 1e-8, - equal_nan: bool = True) -> Union[bool, np.bool_]: + def allclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True): """ Test whether all values are approximately equal to corresponding ones of other Orientation. @@ -223,12 +208,10 @@ class Orientation(Rotation,Crystal): Whether all values are close between both orientations. """ - if not isinstance(other, Orientation): - raise TypeError return np.all(self.isclose(other,rtol,atol,equal_nan)) - def __mul__(self, other: Union[Rotation, "Orientation"]) -> "Orientation": + def __mul__(self,other): """ Compose this orientation with other. @@ -250,7 +233,7 @@ class Orientation(Rotation,Crystal): @staticmethod - def _split_kwargs(kwargs: Dict[str, Any], target: Callable) -> Tuple[Dict[str, Any], ...]: + def _split_kwargs(kwargs,target): """ Separate keyword arguments in 'kwargs' targeted at 'target' from general keyword arguments of Orientation objects. @@ -269,7 +252,7 @@ class Orientation(Rotation,Crystal): Valid keyword arguments of Orientation object. """ - kws: Tuple[Dict[str, Any], ...] = () + kws = () for t in (target,Orientation.__init__): kws += ({key: kwargs[key] for key in set(inspect.signature(t).parameters) & set(kwargs)},) @@ -281,88 +264,85 @@ class Orientation(Rotation,Crystal): @classmethod - @util.extended_docstring(Rotation.from_random, _parameter_doc) - def from_random(cls, **kwargs) -> "Orientation": + @util.extended_docstring(Rotation.from_random,_parameter_doc) + def from_random(cls,**kwargs): kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_random) return cls(rotation=Rotation.from_random(**kwargs_rot),**kwargs_ori) @classmethod @util.extended_docstring(Rotation.from_quaternion,_parameter_doc) - def from_quaternion(cls, **kwargs) -> "Orientation": + def from_quaternion(cls,**kwargs): kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_quaternion) return cls(rotation=Rotation.from_quaternion(**kwargs_rot),**kwargs_ori) @classmethod @util.extended_docstring(Rotation.from_Euler_angles,_parameter_doc) - def from_Euler_angles(cls, **kwargs) -> "Orientation": + def from_Euler_angles(cls,**kwargs): kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_Euler_angles) return cls(rotation=Rotation.from_Euler_angles(**kwargs_rot),**kwargs_ori) @classmethod @util.extended_docstring(Rotation.from_axis_angle,_parameter_doc) - def from_axis_angle(cls, **kwargs) -> "Orientation": + def from_axis_angle(cls,**kwargs): kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_axis_angle) return cls(rotation=Rotation.from_axis_angle(**kwargs_rot),**kwargs_ori) @classmethod @util.extended_docstring(Rotation.from_basis,_parameter_doc) - def from_basis(cls, **kwargs) -> "Orientation": + def from_basis(cls,**kwargs): kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_basis) return cls(rotation=Rotation.from_basis(**kwargs_rot),**kwargs_ori) @classmethod @util.extended_docstring(Rotation.from_matrix,_parameter_doc) - def from_matrix(cls, **kwargs) -> "Orientation": + def from_matrix(cls,**kwargs): kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_matrix) return cls(rotation=Rotation.from_matrix(**kwargs_rot),**kwargs_ori) @classmethod @util.extended_docstring(Rotation.from_Rodrigues_vector,_parameter_doc) - def from_Rodrigues_vector(cls, **kwargs) -> "Orientation": + def from_Rodrigues_vector(cls,**kwargs): kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_Rodrigues_vector) return cls(rotation=Rotation.from_Rodrigues_vector(**kwargs_rot),**kwargs_ori) @classmethod @util.extended_docstring(Rotation.from_homochoric,_parameter_doc) - def from_homochoric(cls, **kwargs) -> "Orientation": + def from_homochoric(cls,**kwargs): kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_homochoric) return cls(rotation=Rotation.from_homochoric(**kwargs_rot),**kwargs_ori) @classmethod @util.extended_docstring(Rotation.from_cubochoric,_parameter_doc) - def from_cubochoric(cls, **kwargs) -> "Orientation": + def from_cubochoric(cls,**kwargs): kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_cubochoric) return cls(rotation=Rotation.from_cubochoric(**kwargs_rot),**kwargs_ori) @classmethod @util.extended_docstring(Rotation.from_spherical_component,_parameter_doc) - def from_spherical_component(cls, **kwargs) -> "Orientation": + def from_spherical_component(cls,**kwargs): kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_spherical_component) return cls(rotation=Rotation.from_spherical_component(**kwargs_rot),**kwargs_ori) @classmethod @util.extended_docstring(Rotation.from_fiber_component,_parameter_doc) - def from_fiber_component(cls, **kwargs) -> "Orientation": + def from_fiber_component(cls,**kwargs): kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_fiber_component) return cls(rotation=Rotation.from_fiber_component(**kwargs_rot),**kwargs_ori) @classmethod @util.extend_docstring(_parameter_doc) - def from_directions(cls, - uvw: Union[Sequence[float], np.ndarray], - hkl: Union[Sequence[float], np.ndarray], - **kwargs) -> "Orientation": + def from_directions(cls,uvw,hkl,**kwargs): """ Initialize orientation object from two crystallographic directions. @@ -382,7 +362,7 @@ class Orientation(Rotation,Crystal): @property - def equivalent(self) -> "Orientation": + def equivalent(self): """ Orientations that are symmetrically equivalent. @@ -396,7 +376,7 @@ class Orientation(Rotation,Crystal): @property - def reduced(self) -> "Orientation": + def reduced(self): """Select symmetrically equivalent orientation that falls into fundamental zone according to symmetry.""" eq = self.equivalent ok = eq.in_FZ @@ -405,8 +385,9 @@ class Orientation(Rotation,Crystal): sort = 0 if len(loc) == 1 else np.lexsort(loc[:0:-1]) return eq[ok][sort].reshape(self.shape) + @property - def in_FZ(self) -> Union[np.bool_, np.ndarray]: + def in_FZ(self): """ Check whether orientation falls into fundamental zone of own symmetry. @@ -450,7 +431,7 @@ class Orientation(Rotation,Crystal): @property - def in_disorientation_FZ(self) -> np.ndarray: + def in_disorientation_FZ(self): """ Check whether orientation falls into fundamental zone of disorientations. @@ -490,7 +471,8 @@ class Orientation(Rotation,Crystal): else: return np.ones_like(rho[...,0],dtype=bool) - def disorientation(self, other, return_operators = False): + + def disorientation(self,other,return_operators=False): """ Calculate disorientation between myself and given other orientation. @@ -536,8 +518,8 @@ class Orientation(Rotation,Crystal): if self.family != other.family: raise NotImplementedError('disorientation between different crystal families') - blend: Tuple[int, ...] = util.shapeblender(self.shape,other.shape) - s = self.equivalent + blend = util.shapeblender(self.shape,other.shape) + s = self.equivalent o = other.equivalent s_ = s.reshape((s.shape[0],1)+ self.shape).broadcast_to((s.shape[0],o.shape[0])+blend,mode='right') @@ -552,9 +534,10 @@ class Orientation(Rotation,Crystal): r = np.where(np.any(forward[...,np.newaxis],axis=(0,1),keepdims=True), r_.quaternion, _r.quaternion) - loc: Tuple[float] = np.where(ok) - sort: np.ndarray = 0 if len(loc) == 2 else np.lexsort(loc[:1:-1]) - quat: np.ndarray = r[ok][sort].reshape(blend+(4,)) + loc = np.where(ok) + sort = 0 if len(loc) == 2 else np.lexsort(loc[:1:-1]) + quat = r[ok][sort].reshape(blend+(4,)) + return ( (self.copy(rotation=quat), (np.vstack(loc[:2]).T)[sort].reshape(blend+(2,))) @@ -563,7 +546,7 @@ class Orientation(Rotation,Crystal): ) - def average(self, weights = None, return_cloud = False): + def average(self,weights=None,return_cloud=False): """ Return orientation average over last dimension. @@ -604,10 +587,7 @@ class Orientation(Rotation,Crystal): ) - def to_SST(self, - vector: np.ndarray, - proper: bool = False, - return_operators: bool = False) -> np.ndarray: + def to_SST(self,vector,proper=False,return_operators=False): """ Rotate vector to ensure it falls into (improper or proper) standard stereographic triangle of crystal symmetry. @@ -646,7 +626,7 @@ class Orientation(Rotation,Crystal): ) - def in_SST(self, vector: np.ndarray, proper: bool = False) -> Union[np.bool_, np.ndarray]: + def in_SST(self,vector,proper=False): """ Check whether given crystal frame vector falls into standard stereographic triangle of own symmetry. @@ -687,7 +667,7 @@ class Orientation(Rotation,Crystal): return np.all(components >= 0.0,axis=-1) - def IPF_color(self, vector: np.ndarray, in_SST: bool = True, proper: bool = False) -> np.ndarray: + def IPF_color(self,vector,in_SST=True,proper=False): """ Map vector to RGB color within standard stereographic triangle of own symmetry. @@ -735,30 +715,30 @@ class Orientation(Rotation,Crystal): components_improper = np.around(np.einsum('...ji,...i', np.broadcast_to(self.standard_triangle['improper'], vector_.shape+(3,)), vector_), 12) - in_SST_ = np.all(components_proper >= 0.0,axis=-1) \ + in_SST = np.all(components_proper >= 0.0,axis=-1) \ | np.all(components_improper >= 0.0,axis=-1) - components = np.where((in_SST_ & np.all(components_proper >= 0.0,axis=-1))[...,np.newaxis], + components = np.where((in_SST & np.all(components_proper >= 0.0,axis=-1))[...,np.newaxis], components_proper,components_improper) else: components = np.around(np.einsum('...ji,...i', np.broadcast_to(self .standard_triangle['improper'], vector_.shape+(3,)), np.block([vector_[...,:2],np.abs(vector_[...,2:3])])), 12) - in_SST_ = np.all(components >= 0.0,axis=-1) + in_SST = np.all(components >= 0.0,axis=-1) with np.errstate(invalid='ignore',divide='ignore'): rgb = (components/np.linalg.norm(components,axis=-1,keepdims=True))**0.5 # smoothen color ramps rgb = np.clip(rgb,0.,1.) # clip intensity rgb /= np.max(rgb,axis=-1,keepdims=True) # normalize to (HS)V = 1 - rgb[np.broadcast_to(~in_SST_[...,np.newaxis],rgb.shape)] = 0.0 + rgb[np.broadcast_to(~in_SST[...,np.newaxis],rgb.shape)] = 0.0 return rgb @property - def symmetry_operations(self) -> Rotation: + def symmetry_operations(self): """Symmetry operations as Rotations.""" - _symmetry_operations: Dict[str, List[List]] = { + _symmetry_operations = { 'cubic': [ [ 1.0, 0.0, 0.0, 0.0 ], [ 0.0, 1.0, 0.0, 0.0 ], @@ -828,10 +808,7 @@ class Orientation(Rotation,Crystal): #################################################################################################### # functions that require lattice, not just family - def to_pole(self, *, - uvw: np.ndarray = None, - hkl: np.ndarray = None, - with_symmetry: bool = False) -> np.ndarray: + def to_pole(self,*,uvw=None,hkl=None,with_symmetry=False): """ Calculate lab frame vector along lattice direction [uvw] or plane normal (hkl). @@ -862,9 +839,7 @@ class Orientation(Rotation,Crystal): @ np.broadcast_to(v,blend+(3,)) - def Schmid(self, *, - N_slip: Sequence[int] = None, - N_twin: Sequence[int] = None) -> np.ndarray: + def Schmid(self,*,N_slip=None,N_twin=None): u""" Calculate Schmid matrix P = d ⨂ n in the lab frame for selected deformation systems. @@ -901,7 +876,6 @@ class Orientation(Rotation,Crystal): (self.kinematics('twin'),N_twin) if active == '*': active = [len(a) for a in kinematics['direction']] - assert active d = self.to_frame(uvw=np.vstack([kinematics['direction'][i][:n] for i,n in enumerate(active)])) p = self.to_frame(hkl=np.vstack([kinematics['plane'][i][:n] for i,n in enumerate(active)])) P = np.einsum('...i,...j',d/np.linalg.norm(d,axis=1,keepdims=True), @@ -912,7 +886,7 @@ class Orientation(Rotation,Crystal): @ np.broadcast_to(P.reshape(util.shapeshifter(P.shape,shape)),shape) - def related(self, model: str) -> "Orientation": + def related(self,model): """ Orientations derived from the given relationship. diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 0a4421090..ac921d70a 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -6,8 +6,6 @@ from . import tensor from . import util from . import grid_filters -from typing import Union, Sequence, Tuple, Literal, List - _P = -1 # parameters for conversion from/to cubochoric @@ -63,7 +61,7 @@ class Rotation: __slots__ = ['quaternion'] - def __init__(self, rotation: Union[Sequence[float], np.ndarray, "Rotation"] = np.array([1.0,0.0,0.0,0.0])): + def __init__(self,rotation = np.array([1.0,0.0,0.0,0.0])): """ New rotation. @@ -75,7 +73,6 @@ class Rotation: Defaults to no rotation. """ - self.quaternion: np.ndarray if isinstance(rotation,Rotation): self.quaternion = rotation.quaternion.copy() elif np.array(rotation).shape[-1] == 4: @@ -84,13 +81,13 @@ class Rotation: raise TypeError('"rotation" is neither a Rotation nor a quaternion') - def __repr__(self) -> str: + def __repr__(self): """Represent rotation as unit quaternion(s).""" return f'Quaternion{" " if self.quaternion.shape == (4,) else "s of shape "+str(self.quaternion.shape[:-1])+chr(10)}'\ + str(self.quaternion) - def __copy__(self, rotation: Union[Sequence[float], np.ndarray, "Rotation"] = None) -> "Rotation": + def __copy__(self,rotation=None): """Create deep copy.""" dup = copy.deepcopy(self) if rotation is not None: @@ -100,14 +97,14 @@ class Rotation: copy = __copy__ - def __getitem__(self, item: Union[Tuple[int], int, bool, np.bool_, np.ndarray]): + def __getitem__(self,item): """Return slice according to item.""" return self.copy() \ if self.shape == () else \ self.copy(rotation=self.quaternion[item+(slice(None),)] if isinstance(item,tuple) else self.quaternion[item]) - def __eq__(self, other: object) -> bool: + def __eq__(self,other): """ Equal to other. @@ -117,13 +114,11 @@ class Rotation: Rotation to check for equality. """ - if not isinstance(other, Rotation): - raise TypeError return np.logical_or(np.all(self.quaternion == other.quaternion,axis=-1), np.all(self.quaternion == -1.0*other.quaternion,axis=-1)) - def __ne__(self, other: object) -> bool: + def __ne__(self,other): """ Not equal to other. @@ -133,12 +128,10 @@ class Rotation: Rotation to check for inequality. """ - if not isinstance(other, Rotation): - raise TypeError return np.logical_not(self==other) - def isclose(self, other: "Rotation", rtol: float = 1e-5, atol: float = 1e-8, equal_nan: bool = True) -> bool: + def isclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True): """ Report where values are approximately equal to corresponding ones of other Rotation. @@ -165,7 +158,7 @@ class Rotation: np.all(np.isclose(s,-1.0*o,rtol,atol,equal_nan),axis=-1)) - def allclose(self, other: "Rotation", rtol: float = 1e-5, atol: float = 1e-8, equal_nan: bool = True) -> Union[np.bool_, bool]: + def allclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True): """ Test whether all values are approximately equal to corresponding ones of other Rotation. @@ -195,27 +188,27 @@ class Rotation: @property - def size(self) -> int: + def size(self): return self.quaternion[...,0].size @property - def shape(self) -> Tuple[int, ...]: + def shape(self): return self.quaternion[...,0].shape - def __len__(self) -> int: + def __len__(self): """Length of leading/leftmost dimension of array.""" return 0 if self.shape == () else self.shape[0] - def __invert__(self) -> "Rotation": + def __invert__(self): """Inverse rotation (backward rotation).""" dup = self.copy() dup.quaternion[...,1:] *= -1 return dup - def __pow__(self, exp: int) -> "Rotation": + def __pow__(self,exp): """ Perform the rotation 'exp' times. @@ -229,7 +222,7 @@ class Rotation: p = self.quaternion[...,1:]/np.linalg.norm(self.quaternion[...,1:],axis=-1,keepdims=True) return self.copy(rotation=Rotation(np.block([np.cos(exp*phi),np.sin(exp*phi)*p]))._standardize()) - def __ipow__(self, exp: int) -> "Rotation": + def __ipow__(self,exp): """ Perform the rotation 'exp' times (in-place). @@ -242,7 +235,7 @@ class Rotation: return self**exp - def __mul__(self, other: "Rotation") -> "Rotation": + def __mul__(self,other): """ Compose with other. @@ -268,7 +261,7 @@ class Rotation: else: raise TypeError('Use "R@b", i.e. matmul, to apply rotation "R" to object "b"') - def __imul__(self, other: "Rotation") -> "Rotation": + def __imul__(self,other): """ Compose with other (in-place). @@ -281,7 +274,7 @@ class Rotation: return self*other - def __truediv__(self, other: "Rotation") -> "Rotation": + def __truediv__(self,other): """ Compose with inverse of other. @@ -301,7 +294,7 @@ class Rotation: else: raise TypeError('Use "R@b", i.e. matmul, to apply rotation "R" to object "b"') - def __itruediv__(self, other: "Rotation") -> "Rotation": + def __itruediv__(self,other): """ Compose with inverse of other (in-place). @@ -314,7 +307,7 @@ class Rotation: return self/other - def __matmul__(self, other: np.ndarray) -> np.ndarray: + def __matmul__(self,other): """ Rotate vector, second order tensor, or fourth order tensor. @@ -357,13 +350,13 @@ class Rotation: apply = __matmul__ - def _standardize(self) -> "Rotation": + def _standardize(self): """Standardize quaternion (ensure positive real hemisphere).""" self.quaternion[self.quaternion[...,0] < 0.0] *= -1 return self - def append(self, other: Union["Rotation", List["Rotation"]]) -> "Rotation": + def append(self,other): """ Extend array along first dimension with other array(s). @@ -376,7 +369,7 @@ class Rotation: [self]+other if isinstance(other,list) else [self,other])))) - def flatten(self, order: Literal['C','F','A'] = 'C') -> "Rotation": + def flatten(self,order = 'C'): """ Flatten array. @@ -389,7 +382,7 @@ class Rotation: return self.copy(rotation=self.quaternion.reshape((-1,4),order=order)) - def reshape(self,shape: Union[int, Tuple[int, ...]], order: Literal['C','F','A'] = 'C') -> "Rotation": + def reshape(self,shape,order = 'C'): """ Reshape array. @@ -403,7 +396,7 @@ class Rotation: return self.copy(rotation=self.quaternion.reshape(tuple(shape)+(4,),order=order)) - def broadcast_to(self, shape: Union[int, Tuple[int, ...]], mode: Literal["left", "right"] = 'right') -> "Rotation": + def broadcast_to(self,shape,mode = 'right'): """ Broadcast array. @@ -465,7 +458,7 @@ class Rotation: accept_homomorph = True) - def misorientation(self, other: "Rotation") -> "Rotation": + def misorientation(self,other): """ Calculate misorientation to other Rotation. @@ -486,7 +479,7 @@ class Rotation: ################################################################################################ # convert to different orientation representations (numpy arrays) - def as_quaternion(self) -> np.ndarray: + def as_quaternion(self): """ Represent as unit quaternion. @@ -499,7 +492,7 @@ class Rotation: return self.quaternion.copy() def as_Euler_angles(self, - degrees: bool = False) -> np.ndarray: + degrees = False): """ Represent as Bunge Euler angles. @@ -533,8 +526,8 @@ class Rotation: return eu def as_axis_angle(self, - degrees: bool = False, - pair: bool = False) -> Union[Tuple[float, ...], np.ndarray]: + degrees = False, + pair = False): """ Represent as axis–angle pair. @@ -561,11 +554,11 @@ class Rotation: (array([0., 0., 1.]), array(0.)) """ - ax: np.ndarray = Rotation._qu2ax(self.quaternion) + ax = Rotation._qu2ax(self.quaternion) if degrees: ax[...,3] = np.degrees(ax[...,3]) return (ax[...,:3],ax[...,3]) if pair else ax - def as_matrix(self) -> np.ndarray: + def as_matrix(self): """ Represent as rotation matrix. @@ -589,7 +582,7 @@ class Rotation: return Rotation._qu2om(self.quaternion) def as_Rodrigues_vector(self, - compact: bool = False) -> np.ndarray: + compact = False): """ Represent as Rodrigues–Frank vector with separate axis and angle argument. @@ -622,7 +615,7 @@ class Rotation: else: return ro - def as_homochoric(self) -> np.ndarray: + def as_homochoric(self): """ Represent as homochoric vector. @@ -643,7 +636,7 @@ class Rotation: """ return Rotation._qu2ho(self.quaternion) - def as_cubochoric(self) -> np.ndarray: + def as_cubochoric(self): """ Represent as cubochoric vector. @@ -668,9 +661,9 @@ class Rotation: # Static constructors. The input data needs to follow the conventions, options allow to # relax the conventions. @staticmethod - def from_quaternion(q: Union[Sequence[Sequence[float]], np.ndarray], - accept_homomorph: bool = False, - P: Literal[1, -1] = -1) -> "Rotation": + def from_quaternion(q, + accept_homomorph = False, + P = -1): """ Initialize from quaternion. @@ -685,7 +678,7 @@ class Rotation: Sign convention. Defaults to -1. """ - qu: np.ndarray = np.array(q,dtype=float) + qu = np.array(q,dtype=float) if qu.shape[:-2:-1] != (4,): raise ValueError('Invalid shape.') if abs(P) != 1: @@ -703,8 +696,8 @@ class Rotation: return Rotation(qu) @staticmethod - def from_Euler_angles(phi: np.ndarray, - degrees: bool = False) -> "Rotation": + def from_Euler_angles(phi, + degrees = False): """ Initialize from Bunge Euler angles. @@ -732,10 +725,10 @@ class Rotation: return Rotation(Rotation._eu2qu(eu)) @staticmethod - def from_axis_angle(axis_angle: np.ndarray, - degrees:bool = False, - normalize: bool = False, - P: Literal[1, -1] = -1) -> "Rotation": + def from_axis_angle(axis_angle, + degrees = False, + normalize = False, + P = -1): """ Initialize from Axis angle pair. @@ -770,9 +763,9 @@ class Rotation: return Rotation(Rotation._ax2qu(ax)) @staticmethod - def from_basis(basis: np.ndarray, - orthonormal: bool = True, - reciprocal: bool = False) -> "Rotation": + def from_basis(basis, + orthonormal = True, + reciprocal = False): """ Initialize from lattice basis vectors. @@ -806,7 +799,7 @@ class Rotation: return Rotation(Rotation._om2qu(om)) @staticmethod - def from_matrix(R: np.ndarray) -> "Rotation": + def from_matrix(R): """ Initialize from rotation matrix. @@ -819,9 +812,8 @@ class Rotation: return Rotation.from_basis(R) @staticmethod - def from_parallel(a: np.ndarray, - b: np.ndarray, - **kwargs) -> "Rotation": + def from_parallel(a,b, + **kwargs): """ Initialize from pairs of two orthogonal lattice basis vectors. @@ -849,9 +841,9 @@ class Rotation: @staticmethod - def from_Rodrigues_vector(rho: np.ndarray, - normalize: bool = False, - P: Literal[1, -1] = -1) -> "Rotation": + def from_Rodrigues_vector(rho, + normalize = False, + P = -1): """ Initialize from Rodrigues–Frank vector (angle separated from axis). @@ -881,8 +873,8 @@ class Rotation: return Rotation(Rotation._ro2qu(ro)) @staticmethod - def from_homochoric(h: np.ndarray, - P: Literal[1, -1] = -1) -> "Rotation": + def from_homochoric(h, + P = -1): """ Initialize from homochoric vector. @@ -908,8 +900,8 @@ class Rotation: return Rotation(Rotation._ho2qu(ho)) @staticmethod - def from_cubochoric(x: np.ndarray, - P: Literal[1, -1] = -1) -> "Rotation": + def from_cubochoric(x, + P = -1): """ Initialize from cubochoric vector. @@ -936,8 +928,8 @@ class Rotation: @staticmethod - def from_random(shape: Tuple[int, ...] = None, - rng_seed: Union[None, int, Sequence[int], np.ndarray] = None) -> "Rotation": + def from_random(shape = None, + rng_seed = None): """ Initialize with random rotation. @@ -966,13 +958,13 @@ class Rotation: @staticmethod - def from_ODF(weights: np.ndarray, - phi: np.ndarray, - N: int = 500, - degrees: bool = True, - fractions: bool = True, - rng_seed: Union[None, int, Sequence[int], np.ndarray] = None, - **kwargs) -> "Rotation": + def from_ODF(weights, + phi, + N = 500, + degrees = True, + fractions = True, + rng_seed = None, + **kwargs): """ Sample discrete values from a binned orientation distribution function (ODF). @@ -1025,11 +1017,11 @@ class Rotation: @staticmethod - def from_spherical_component(center: "Rotation", - sigma: float, - N: int = 500, - degrees: bool = True, - rng_seed: Union[None, int, Sequence[float], np.ndarray] = None) -> "Rotation": + def from_spherical_component(center, + sigma, + N = 500, + degrees = True, + rng_seed = None): """ Calculate set of rotations with Gaussian distribution around center. @@ -1060,12 +1052,12 @@ class Rotation: @staticmethod - def from_fiber_component(alpha: np.ndarray, - beta: np.ndarray, - sigma: float = 0.0, - N: int = 500, - degrees: bool = True, - rng_seed: bool = None): + def from_fiber_component(alpha, + beta, + sigma = 0.0, + N = 500, + degrees = True, + rng_seed = None): """ Calculate set of rotations with Gaussian distribution around direction. @@ -1088,8 +1080,7 @@ class Rotation: """ rng = np.random.default_rng(rng_seed) - sigma_: np.ndarray; alpha_: np.ndarray; beta_: np.ndarray - sigma_,alpha_,beta_ = map(np.radians,(sigma,alpha,beta)) if degrees else (sigma,alpha,beta) # type: ignore + sigma_,alpha_,beta_ = map(np.radians,(sigma,alpha,beta)) if degrees else (sigma,alpha,beta) d_cr = np.array([np.sin(alpha_[0])*np.cos(alpha_[1]), np.sin(alpha_[0])*np.sin(alpha_[1]), np.cos(alpha_[0])]) d_lab = np.array([np.sin( beta_[0])*np.cos( beta_[1]), np.sin( beta_[0])*np.sin( beta_[1]), np.cos( beta_[0])]) @@ -1143,7 +1134,7 @@ class Rotation: #################################################################################################### #---------- Quaternion ---------- @staticmethod - def _qu2om(qu: np.ndarray) -> np.ndarray: + def _qu2om(qu): qq = qu[...,0:1]**2-(qu[...,1:2]**2 + qu[...,2:3]**2 + qu[...,3:4]**2) om = np.block([qq + 2.0*qu[...,1:2]**2, 2.0*(qu[...,2:3]*qu[...,1:2]-_P*qu[...,0:1]*qu[...,3:4]), @@ -1158,7 +1149,7 @@ class Rotation: return om @staticmethod - def _qu2eu(qu: np.ndarray) -> np.ndarray: + def _qu2eu(qu): """Quaternion to Bunge Euler angles.""" q02 = qu[...,0:1]*qu[...,2:3] q13 = qu[...,1:2]*qu[...,3:4] @@ -1186,7 +1177,7 @@ class Rotation: return eu @staticmethod - def _qu2ax(qu: np.ndarray) -> np.ndarray: + def _qu2ax(qu): """ Quaternion to axis–angle pair. @@ -1202,7 +1193,7 @@ class Rotation: return ax @staticmethod - def _qu2ro(qu: np.ndarray) -> np.ndarray: + def _qu2ro(qu): """Quaternion to Rodrigues–Frank vector.""" with np.errstate(invalid='ignore',divide='ignore'): s = np.linalg.norm(qu[...,1:4],axis=-1,keepdims=True) @@ -1216,7 +1207,7 @@ class Rotation: return ro @staticmethod - def _qu2ho(qu: np.ndarray) -> np.ndarray: + def _qu2ho(qu): """Quaternion to homochoric vector.""" with np.errstate(invalid='ignore'): omega = 2.0 * np.arccos(np.clip(qu[...,0:1],-1.0,1.0)) @@ -1227,14 +1218,14 @@ class Rotation: return ho @staticmethod - def _qu2cu(qu: np.ndarray) -> np.ndarray: + def _qu2cu(qu): """Quaternion to cubochoric vector.""" return Rotation._ho2cu(Rotation._qu2ho(qu)) #---------- Rotation matrix ---------- @staticmethod - def _om2qu(om: np.ndarray) -> np.ndarray: + def _om2qu(om): """ Rotation matrix to quaternion. @@ -1276,7 +1267,7 @@ class Rotation: return qu @staticmethod - def _om2eu(om: np.ndarray) -> np.ndarray: + def _om2eu(om): """Rotation matrix to Bunge Euler angles.""" with np.errstate(invalid='ignore',divide='ignore'): zeta = 1.0/np.sqrt(1.0-om[...,2,2:3]**2) @@ -1295,7 +1286,7 @@ class Rotation: return eu @staticmethod - def _om2ax(om: np.ndarray) -> np.ndarray: + def _om2ax(om): """Rotation matrix to axis–angle pair.""" diag_delta = -_P*np.block([om[...,1,2:3]-om[...,2,1:2], om[...,2,0:1]-om[...,0,2:3], @@ -1316,24 +1307,24 @@ class Rotation: return ax @staticmethod - def _om2ro(om: np.ndarray) -> np.ndarray: + def _om2ro(om): """Rotation matrix to Rodrigues–Frank vector.""" return Rotation._eu2ro(Rotation._om2eu(om)) @staticmethod - def _om2ho(om: np.ndarray) -> np.ndarray: + def _om2ho(om): """Rotation matrix to homochoric vector.""" return Rotation._ax2ho(Rotation._om2ax(om)) @staticmethod - def _om2cu(om: np.ndarray) -> np.ndarray: + def _om2cu(om): """Rotation matrix to cubochoric vector.""" return Rotation._ho2cu(Rotation._om2ho(om)) #---------- Bunge Euler angles ---------- @staticmethod - def _eu2qu(eu: np.ndarray) -> np.ndarray: + def _eu2qu(eu): """Bunge Euler angles to quaternion.""" ee = 0.5*eu cPhi = np.cos(ee[...,1:2]) @@ -1346,7 +1337,7 @@ class Rotation: return qu @staticmethod - def _eu2om(eu: np.ndarray) -> np.ndarray: + def _eu2om(eu): """Bunge Euler angles to rotation matrix.""" c = np.cos(eu) s = np.sin(eu) @@ -1364,7 +1355,7 @@ class Rotation: return om @staticmethod - def _eu2ax(eu: np.ndarray) -> np.ndarray: + def _eu2ax(eu): """Bunge Euler angles to axis–angle pair.""" t = np.tan(eu[...,1:2]*0.5) sigma = 0.5*(eu[...,0:1]+eu[...,2:3]) @@ -1383,7 +1374,7 @@ class Rotation: return ax @staticmethod - def _eu2ro(eu: np.ndarray) -> np.ndarray: + def _eu2ro(eu): """Bunge Euler angles to Rodrigues–Frank vector.""" ax = Rotation._eu2ax(eu) ro = np.block([ax[...,:3],np.tan(ax[...,3:4]*.5)]) @@ -1392,19 +1383,19 @@ class Rotation: return ro @staticmethod - def _eu2ho(eu: np.ndarray) -> np.ndarray: + def _eu2ho(eu): """Bunge Euler angles to homochoric vector.""" return Rotation._ax2ho(Rotation._eu2ax(eu)) @staticmethod - def _eu2cu(eu: np.ndarray) -> np.ndarray: + def _eu2cu(eu): """Bunge Euler angles to cubochoric vector.""" return Rotation._ho2cu(Rotation._eu2ho(eu)) #---------- Axis angle pair ---------- @staticmethod - def _ax2qu(ax: np.ndarray) -> np.ndarray: + def _ax2qu(ax): """Axis–angle pair to quaternion.""" c = np.cos(ax[...,3:4]*.5) s = np.sin(ax[...,3:4]*.5) @@ -1412,7 +1403,7 @@ class Rotation: return qu @staticmethod - def _ax2om(ax: np.ndarray) -> np.ndarray: + def _ax2om(ax): """Axis-angle pair to rotation matrix.""" c = np.cos(ax[...,3:4]) s = np.sin(ax[...,3:4]) @@ -1429,12 +1420,12 @@ class Rotation: return om if _P < 0.0 else np.swapaxes(om,-1,-2) @staticmethod - def _ax2eu(ax: np.ndarray) -> np.ndarray: + def _ax2eu(ax): """Rotation matrix to Bunge Euler angles.""" return Rotation._om2eu(Rotation._ax2om(ax)) @staticmethod - def _ax2ro(ax: np.ndarray) -> np.ndarray: + def _ax2ro(ax): """Axis–angle pair to Rodrigues–Frank vector.""" ro = np.block([ax[...,:3], np.where(np.isclose(ax[...,3:4],np.pi,atol=1.e-15,rtol=.0), @@ -1445,36 +1436,36 @@ class Rotation: return ro @staticmethod - def _ax2ho(ax: np.ndarray) -> np.ndarray: + def _ax2ho(ax): """Axis–angle pair to homochoric vector.""" f = (0.75 * ( ax[...,3:4] - np.sin(ax[...,3:4]) ))**(1.0/3.0) ho = ax[...,:3] * f return ho @staticmethod - def _ax2cu(ax: np.ndarray) -> np.ndarray: + def _ax2cu(ax): """Axis–angle pair to cubochoric vector.""" return Rotation._ho2cu(Rotation._ax2ho(ax)) #---------- Rodrigues-Frank vector ---------- @staticmethod - def _ro2qu(ro: np.ndarray) -> np.ndarray: + def _ro2qu(ro): """Rodrigues–Frank vector to quaternion.""" return Rotation._ax2qu(Rotation._ro2ax(ro)) @staticmethod - def _ro2om(ro: np.ndarray) -> np.ndarray: + def _ro2om(ro): """Rodgrigues–Frank vector to rotation matrix.""" return Rotation._ax2om(Rotation._ro2ax(ro)) @staticmethod - def _ro2eu(ro: np.ndarray) -> np.ndarray: + def _ro2eu(ro): """Rodrigues–Frank vector to Bunge Euler angles.""" return Rotation._om2eu(Rotation._ro2om(ro)) @staticmethod - def _ro2ax(ro: np.ndarray) -> np.ndarray: + def _ro2ax(ro): """Rodrigues–Frank vector to axis–angle pair.""" with np.errstate(invalid='ignore',divide='ignore'): ax = np.where(np.isfinite(ro[...,3:4]), @@ -1484,7 +1475,7 @@ class Rotation: return ax @staticmethod - def _ro2ho(ro: np.ndarray) -> np.ndarray: + def _ro2ho(ro): """Rodrigues–Frank vector to homochoric vector.""" f = np.where(np.isfinite(ro[...,3:4]),2.0*np.arctan(ro[...,3:4]) -np.sin(2.0*np.arctan(ro[...,3:4])),np.pi) ho = np.where(np.broadcast_to(np.sum(ro[...,0:3]**2.0,axis=-1,keepdims=True) < 1.e-8,ro[...,0:3].shape), @@ -1492,29 +1483,29 @@ class Rotation: return ho @staticmethod - def _ro2cu(ro: np.ndarray) -> np.ndarray: + def _ro2cu(ro): """Rodrigues–Frank vector to cubochoric vector.""" return Rotation._ho2cu(Rotation._ro2ho(ro)) #---------- Homochoric vector---------- @staticmethod - def _ho2qu(ho: np.ndarray) -> np.ndarray: + def _ho2qu(ho): """Homochoric vector to quaternion.""" return Rotation._ax2qu(Rotation._ho2ax(ho)) @staticmethod - def _ho2om(ho: np.ndarray) -> np.ndarray: + def _ho2om(ho): """Homochoric vector to rotation matrix.""" return Rotation._ax2om(Rotation._ho2ax(ho)) @staticmethod - def _ho2eu(ho: np.ndarray) -> np.ndarray: + def _ho2eu(ho): """Homochoric vector to Bunge Euler angles.""" return Rotation._ax2eu(Rotation._ho2ax(ho)) @staticmethod - def _ho2ax(ho: np.ndarray) -> np.ndarray: + def _ho2ax(ho): """Homochoric vector to axis–angle pair.""" tfit = np.array([+1.0000000000018852, -0.5000000002194847, -0.024999992127593126, -0.003928701544781374, @@ -1537,12 +1528,12 @@ class Rotation: return ax @staticmethod - def _ho2ro(ho: np.ndarray) -> np.ndarray: + def _ho2ro(ho): """Axis–angle pair to Rodrigues–Frank vector.""" return Rotation._ax2ro(Rotation._ho2ax(ho)) @staticmethod - def _ho2cu(ho: np.ndarray) -> np.ndarray: + def _ho2cu(ho): """ Homochoric vector to cubochoric vector. @@ -1582,32 +1573,32 @@ class Rotation: #---------- Cubochoric ---------- @staticmethod - def _cu2qu(cu: np.ndarray) -> np.ndarray: + def _cu2qu(cu): """Cubochoric vector to quaternion.""" return Rotation._ho2qu(Rotation._cu2ho(cu)) @staticmethod - def _cu2om(cu: np.ndarray) -> np.ndarray: + def _cu2om(cu): """Cubochoric vector to rotation matrix.""" return Rotation._ho2om(Rotation._cu2ho(cu)) @staticmethod - def _cu2eu(cu: np.ndarray) -> np.ndarray: + def _cu2eu(cu): """Cubochoric vector to Bunge Euler angles.""" return Rotation._ho2eu(Rotation._cu2ho(cu)) @staticmethod - def _cu2ax(cu: np.ndarray) -> np.ndarray: + def _cu2ax(cu): """Cubochoric vector to axis–angle pair.""" return Rotation._ho2ax(Rotation._cu2ho(cu)) @staticmethod - def _cu2ro(cu: np.ndarray) -> np.ndarray: + def _cu2ro(cu): """Cubochoric vector to Rodrigues–Frank vector.""" return Rotation._ho2ro(Rotation._cu2ho(cu)) @staticmethod - def _cu2ho(cu: np.ndarray) -> np.ndarray: + def _cu2ho(cu): """ Cubochoric vector to homochoric vector. @@ -1650,7 +1641,7 @@ class Rotation: @staticmethod - def _get_pyramid_order(xyz: np.ndarray, direction: Literal['forward', 'backward']) -> np.ndarray: + def _get_pyramid_order(xyz,direction=None): """ Get order of the coordinates. From 9c79c29288a4e49aa85535909e8dd7249e374823 Mon Sep 17 00:00:00 2001 From: Nikhil Prabhu Date: Tue, 14 Dec 2021 18:10:11 +0100 Subject: [PATCH 27/74] Lattice, density and reference --- PRIVATE | 2 +- examples/config/phase/Ag.yaml | 4 ++++ 2 files changed, 5 insertions(+), 1 deletion(-) create mode 100644 examples/config/phase/Ag.yaml diff --git a/PRIVATE b/PRIVATE index e6e1f93a3..2ad27552c 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit e6e1f93a36d63348359a81d7c373083a39977694 +Subproject commit 2ad27552c43316735b6ef425737fe3c8a5231598 diff --git a/examples/config/phase/Ag.yaml b/examples/config/phase/Ag.yaml new file mode 100644 index 000000000..ebcd21082 --- /dev/null +++ b/examples/config/phase/Ag.yaml @@ -0,0 +1,4 @@ +references: + - https://en.wikipedia.org/wiki/Silver +lattice: cF +rho: 10490.0 From d30d3b21283708b8d84acac7dfa289dbe1eff324 Mon Sep 17 00:00:00 2001 From: Test User Date: Tue, 14 Dec 2021 20:21:48 +0100 Subject: [PATCH 28/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-275-g3b9822ad6 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 27dd17c2c..959348642 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-272-g3192a31e1 +v3.0.0-alpha5-275-g3b9822ad6 From 67cc36daf33a4b0de6b4b75ba359c9b666a22358 Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 15 Dec 2021 12:04:02 +0100 Subject: [PATCH 29/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-283-gdacd08f39 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 959348642..7ba8f4262 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-275-g3b9822ad6 +v3.0.0-alpha5-283-gdacd08f39 From 99267afda3d01203758959cd7e45c6a1f0093335 Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 15 Dec 2021 14:21:37 +0100 Subject: [PATCH 30/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-289-g6eccc6a4a --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 959348642..91d183b4a 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-275-g3b9822ad6 +v3.0.0-alpha5-289-g6eccc6a4a From f40d731fe1fa4fcd62b99de1f0e1b8ecd2501cf9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 16 Dec 2021 18:22:42 +0100 Subject: [PATCH 31/74] use the Box-Muller transform instead of random sampling still needs testing. --- src/math.f90 | 47 +++++++++++------------ src/phase_mechanical_plastic_nonlocal.f90 | 24 +++++------- 2 files changed, 31 insertions(+), 40 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 207a27e78..3303f3b4c 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -961,39 +961,36 @@ pure function math_3333toVoigt66(m3333) end function math_3333toVoigt66 - !-------------------------------------------------------------------------------------------------- -!> @brief draw a random sample from Gauss variable +!> @brief Draw a sample from a normal distribution. +!> @details https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform +!> https://masuday.github.io/fortran_tutorial/random.html !-------------------------------------------------------------------------------------------------- -real(pReal) function math_sampleGaussVar(mu, sigma, width) +impure elemental subroutine math_normal(x,mu,sigma) - real(pReal), intent(in) :: mu, & !< mean - sigma !< standard deviation - real(pReal), intent(in), optional :: width !< cut off as multiples of standard deviation + real(pReal), intent(out) :: x + real(pReal), intent(in), optional :: mu, sigma + + real(pReal) :: sigma_, mu_ + real(pReal), dimension(2) :: rnd - real(pReal), dimension(2) :: rnd ! random numbers - real(pReal) :: scatter, & ! normalized scatter around mean - width_ - - if (abs(sigma) < tol_math_check) then - math_sampleGaussVar = mu + + if (present(mu)) then + mu_ = mu else - if (present(width)) then - width_ = width - else - width_ = 3.0_pReal ! use +-3*sigma as default scatter - endif + mu_ = 0.0_pReal + end if - do - call random_number(rnd) - scatter = width_ * (2.0_pReal * rnd(1) - 1.0_pReal) - if (rnd(2) <= exp(-0.5_pReal * scatter**2)) exit ! test if scattered value is drawn - enddo + if (present(sigma)) then + sigma_ = sigma + else + sigma_ = 1.0_pReal + end if - math_sampleGaussVar = scatter * sigma - endif + call random_number(rnd) + x = mu_ + sigma_ * sqrt(-2.0_pReal*log(1.0_pReal-rnd(1)))*cos(2.0_pReal*PI*(1.0_pReal - rnd(2))) -end function math_sampleGaussVar +end subroutine math_normal !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index c762039d4..bf0dc3b19 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -1592,21 +1592,15 @@ subroutine stateInit(ini,phase,Nentries) stt%rhoSglMobile(s,e) = densityBinning end do else ! homogeneous distribution with noise - do e = 1, Nentries - do f = 1,size(ini%N_sl,1) - from = 1 + sum(ini%N_sl(1:f-1)) - upto = sum(ini%N_sl(1:f)) - do s = from,upto - noise = [math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u), & - math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u)] - stt%rho_sgl_mob_edg_pos(s,e) = ini%rho_u_ed_pos_0(f) + noise(1) - stt%rho_sgl_mob_edg_neg(s,e) = ini%rho_u_ed_neg_0(f) + noise(1) - stt%rho_sgl_mob_scr_pos(s,e) = ini%rho_u_sc_pos_0(f) + noise(2) - stt%rho_sgl_mob_scr_neg(s,e) = ini%rho_u_sc_neg_0(f) + noise(2) - end do - stt%rho_dip_edg(from:upto,e) = ini%rho_d_ed_0(f) - stt%rho_dip_scr(from:upto,e) = ini%rho_d_sc_0(f) - end do + do f = 1,size(ini%N_sl,1) + from = 1 + sum(ini%N_sl(1:f-1)) + upto = sum(ini%N_sl(1:f)) + call math_normal(stt%rho_sgl_mob_edg_pos(from:upto,:),ini%rho_u_ed_pos_0(f),ini%sigma_rho_u) + call math_normal(stt%rho_sgl_mob_edg_neg(from:upto,:),ini%rho_u_ed_neg_0(f),ini%sigma_rho_u) + call math_normal(stt%rho_sgl_mob_scr_pos(from:upto,:),ini%rho_u_sc_pos_0(f),ini%sigma_rho_u) + call math_normal(stt%rho_sgl_mob_scr_neg(from:upto,:),ini%rho_u_sc_neg_0(f),ini%sigma_rho_u) + stt%rho_dip_edg(from:upto,:) = ini%rho_d_ed_0(f) + stt%rho_dip_scr(from:upto,:) = ini%rho_d_sc_0(f) end do end if From f833d348e0337947c3852b591f4af98b5e2c7cdc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 17 Dec 2021 08:01:15 +0100 Subject: [PATCH 32/74] testing random sampling --- PRIVATE | 2 +- src/math.f90 | 28 ++++++++++++++++++++++++---- 2 files changed, 25 insertions(+), 5 deletions(-) diff --git a/PRIVATE b/PRIVATE index 2ad27552c..96c32ba42 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 2ad27552c43316735b6ef425737fe3c8a5231598 +Subproject commit 96c32ba4237a51eaad92cd139e1a716ee5b32493 diff --git a/src/math.f90 b/src/math.f90 index 3303f3b4c..2d8f8fb3b 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -895,7 +895,7 @@ pure function math_33toVoigt6_stress(sigma) result(sigma_tilde) sigma_tilde = [sigma(1,1), sigma(2,2), sigma(3,3), & - sigma(3,2), sigma(3,1), sigma(1,2)] + sigma(3,2), sigma(3,1), sigma(1,2)] end function math_33toVoigt6_stress @@ -910,7 +910,7 @@ pure function math_33toVoigt6_strain(epsilon) result(epsilon_tilde) epsilon_tilde = [ epsilon(1,1), epsilon(2,2), epsilon(3,3), & - 2.0_pReal*epsilon(3,2), 2.0_pReal*epsilon(3,1), 2.0_pReal*epsilon(1,2)] + 2.0_pReal*epsilon(3,2), 2.0_pReal*epsilon(3,1), 2.0_pReal*epsilon(1,2)] end function math_33toVoigt6_strain @@ -970,11 +970,11 @@ impure elemental subroutine math_normal(x,mu,sigma) real(pReal), intent(out) :: x real(pReal), intent(in), optional :: mu, sigma - + real(pReal) :: sigma_, mu_ real(pReal), dimension(2) :: rnd - + if (present(mu)) then mu_ = mu else @@ -1431,6 +1431,26 @@ subroutine selfTest if (dNeq0(math_LeviCivita(ijk(1),ijk(2),ijk(3)))) & error stop 'math_LeviCivita' + normal_distribution: block + real(pReal), dimension(500000) :: r + real(pReal) :: mu, sigma + + call random_number(mu) + call random_number(sigma) + + sigma = 1.0_pReal + sigma*5.0_pReal + mu = (mu-0.5_pReal)*10_pReal + + call math_normal(r,mu,sigma) + + if (abs(mu -sum(r)/real(size(r),pReal))>5.0e-2_pReal) & + error stop 'math_normal(mu)' + + mu = sum(r)/real(size(r),pReal) + if (abs(sigma**2 -1.0_pReal/real(size(r)-1,pReal) * sum((r-mu)**2))/sigma > 5.0e-2_pReal) & + error stop 'math_normal(sigma)' + end block normal_distribution + end subroutine selfTest end module math From 501465dfd1d6d456faed9c49032297dd474537d4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 17 Dec 2021 09:42:10 +0100 Subject: [PATCH 33/74] more alpha releases than expected hopefully, all DAMASK 2 users have been migrated once 3.0 is released --- python/damask/_grid.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/damask/_grid.py b/python/damask/_grid.py index 135cc6b66..d852f0642 100644 --- a/python/damask/_grid.py +++ b/python/damask/_grid.py @@ -197,7 +197,7 @@ class Grid: Grid-based geometry from file. """ - warnings.warn('Support for ASCII-based geom format will be removed in DAMASK 3.1.0', DeprecationWarning,2) + warnings.warn('Support for ASCII-based geom format will be removed in DAMASK 3.0.0', DeprecationWarning,2) try: f = open(fname) except TypeError: @@ -629,7 +629,7 @@ class Grid: Compress geometry with 'x of y' and 'a to b'. """ - warnings.warn('Support for ASCII-based geom format will be removed in DAMASK 3.1.0', DeprecationWarning,2) + warnings.warn('Support for ASCII-based geom format will be removed in DAMASK 3.0.0', DeprecationWarning,2) header = [f'{len(self.comments)+4} header'] + self.comments \ + ['grid a {} b {} c {}'.format(*self.cells), 'size x {} y {} z {}'.format(*self.size), From 0b6af19e5468d4867baba1e7b512363eb71fbffb Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 17 Dec 2021 09:47:58 +0100 Subject: [PATCH 34/74] use of keywords is more intuitive code needs cleaning after revoming of 'what' and 'datasets'. For the moment, keep the old style for compatibility with existing evaluation scripts --- python/damask/_result.py | 115 +++++++++++++++++++++++++++++------- python/tests/test_Result.py | 22 +++---- 2 files changed, 104 insertions(+), 33 deletions(-) diff --git a/python/damask/_result.py b/python/damask/_result.py index 141b795c2..981b0dcb3 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -4,6 +4,7 @@ import fnmatch import os import copy import datetime +import warnings import xml.etree.ElementTree as ET import xml.dom.minidom from pathlib import Path @@ -27,6 +28,20 @@ h5py3 = h5py.__version__[0] == '3' chunk_size = 1024**2//8 # for compression in HDF5 +def _view_transition(what,datasets,increments,times,phases,homogenizations,fields): + if (datasets is not None and what is None) or (what is not None and datasets is None): + raise ValueError('"what" and "datasets" need to be used as a pair') + if datasets is not None or what is not None: + warnings.warn('Arguments "what" and "datasets" will be removed in DAMASK v3.0.0-alpha7', DeprecationWarning,2) + return what,datasets + if sum(1 for _ in filter(None.__ne__, [increments,times,phases,homogenizations,fields])) > 1: + raise ValueError('Only one out of "increments", "times", "phases", "homogenizations", and "fields" can be used') + else: + if increments is not None: return "increments", increments + if times is not None: return "times", times + if phases is not None: return "phases", phases + if homogenizations is not None: return "homogenizations", homogenizations + if fields is not None: return "fields", fields def _read(dataset): """Read a dataset and its metadata into a numpy.ndarray.""" @@ -79,7 +94,7 @@ class Result: >>> r.add_Cauchy() >>> r.add_equivalent_Mises('sigma') >>> r.export_VTK() - >>> r_last = r.view('increments',-1) + >>> r_last = r.view(increments=-1) >>> sigma_vM_last = r_last.get('sigma_vM') """ @@ -155,10 +170,10 @@ class Result: """Show summary of file content.""" visible_increments = self.visible['increments'] - first = self.view('increments',visible_increments[0:1]).list_data() + first = self.view(increments=visible_increments[0:1]).list_data() last = '' if len(visible_increments) < 2 else \ - self.view('increments',visible_increments[-1:]).list_data() + self.view(increments=visible_increments[-1:]).list_data() in_between = '' if len(visible_increments) < 3 else \ ''.join([f'\n{inc}\n ...\n' for inc in visible_increments[1:-1]]) @@ -285,7 +300,6 @@ class Result: selected.append(self.increments[i]) return selected - def times_in_range(self,start,end): """ Get all increments within a given time range. @@ -310,17 +324,35 @@ class Result: return selected - def view(self,what,datasets): + def view(self,what=None,datasets=None,*, + increments=None, + times=None, + phases=None, + homogenizations=None, + fields=None): """ Set view. + Wildcard matching with '?' and '*' is supported. + True is equivalent to '*', False is equivalent to []. + Parameters ---------- what : {'increments', 'times', 'phases', 'homogenizations', 'fields'} - Attribute to change. + Attribute to change. DEPRECATED. datasets : (list of) int (for increments), (list of) float (for times), (list of) str, or bool - Name of datasets; supports '?' and '*' wildcards. + Name of datasets; supports '?' and '*' wildcards. DEPRECATED. True is equivalent to '*', False is equivalent to []. + increments: (list of) int, (list of) str, or bool, optional. + Number(s) of increments to select. + times: (list of) float, (list of) str, or bool, optional. + Simulation time(s) of increments to select. + phases: (list of) str, or bool, optional. + Name(s) of phases to select. + homogenizations: (list of) str, or bool, optional. + Name(s) of homogenizations to select. + fields: (list of) str, or bool, optional. + Name(s) of fields to select. Returns ------- @@ -333,29 +365,48 @@ class Result: >>> import damask >>> r = damask.Result('my_file.hdf5') - >>> r_first = r.view('increment',0) + >>> r_first = r.view(increment=0) Get a view that shows all results between simulation times of 10 to 40: >>> import damask >>> r = damask.Result('my_file.hdf5') - >>> r_t10to40 = r.view('times',r.times_in_range(10.0,40.0)) + >>> r_t10to40 = r.view(times=r.times_in_range(10.0,40.0)) """ - return self._manage_view('set',what,datasets) + what_, datasets_ = _view_transition(what,datasets,increments,times,phases,homogenizations,fields) + return self._manage_view('set',what_,datasets_) - def view_more(self,what,datasets): + def view_more(self,what=None,datasets=None,*, + increments=None, + times=None, + phases=None, + homogenizations=None, + fields=None): """ Add to view. + Wildcard matching with '?' and '*' is supported. + True is equivalent to '*', False is equivalent to []. + Parameters ---------- what : {'increments', 'times', 'phases', 'homogenizations', 'fields'} - Attribute to change. + Attribute to change. DEPRECATED. datasets : (list of) int (for increments), (list of) float (for times), (list of) str, or bool - Name of datasets; supports '?' and '*' wildcards. + Name of datasets; supports '?' and '*' wildcards. DEPRECATED. True is equivalent to '*', False is equivalent to []. + increments: (list of) int, (list of) str, or bool, optional. + Number(s) of increments to select. + times: (list of) float, (list of) str, or bool, optional. + Simulation time(s) of increments to select. + phases: (list of) str, or bool, optional. + Name(s) of phases to select. + homogenizations: (list of) str, or bool, optional. + Name(s) of homogenizations to select. + fields: (list of) str, or bool, optional. + Name(s) of fields to select. Returns ------- @@ -367,25 +418,44 @@ class Result: Get a view that shows only results from first and last increment: >>> import damask - >>> r_empty = damask.Result('my_file.hdf5').view('increments',False) - >>> r_first = r_empty.view_more('increments',0) - >>> r_first_and_last = r.first.view_more('increments',-1) + >>> r_empty = damask.Result('my_file.hdf5').view(increments=False) + >>> r_first = r_empty.view_more(increments=0) + >>> r_first_and_last = r.first.view_more(increments=-1) """ - return self._manage_view('add',what,datasets) + what_, datasets_ = _view_transition(what,datasets,increments,times,phases,homogenizations,fields) + return self._manage_view('add',what_,datasets_) - def view_less(self,what,datasets): + def view_less(self,what=None,datasets=None,*, + increments=None, + times=None, + phases=None, + homogenizations=None, + fields=None): """ Remove from view. + Wildcard matching with '?' and '*' is supported. + True is equivalent to '*', False is equivalent to []. + Parameters ---------- what : {'increments', 'times', 'phases', 'homogenizations', 'fields'} - Attribute to change. + Attribute to change. DEPRECATED. datasets : (list of) int (for increments), (list of) float (for times), (list of) str, or bool - Name of datasets; supports '?' and '*' wildcards. + Name of datasets; supports '?' and '*' wildcards. DEPRECATED. True is equivalent to '*', False is equivalent to []. + increments: (list of) int, (list of) str, or bool, optional. + Number(s) of increments to select. + times: (list of) float, (list of) str, or bool, optional. + Simulation time(s) of increments to select. + phases: (list of) str, or bool, optional. + Name(s) of phases to select. + homogenizations: (list of) str, or bool, optional. + Name(s) of homogenizations to select. + fields: (list of) str, or bool, optional. + Name(s) of fields to select. Returns ------- @@ -398,10 +468,11 @@ class Result: >>> import damask >>> r_all = damask.Result('my_file.hdf5') - >>> r_deformed = r_all.view_less('increments',0) + >>> r_deformed = r_all.view_less(increments=0) """ - return self._manage_view('del',what,datasets) + what_, datasets_ = _view_transition(what,datasets,increments,times,phases,homogenizations,fields) + return self._manage_view('del',what_,datasets_) def rename(self,name_src,name_dst): diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index 3f2555aec..3bc363ebc 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -25,7 +25,7 @@ def default(tmp_path,ref_path): fname = '12grains6x7x8_tensionY.hdf5' shutil.copy(ref_path/fname,tmp_path) f = Result(tmp_path/fname) - return f.view('times',20.0) + return f.view(times=20.0) @pytest.fixture def single_phase(tmp_path,ref_path): @@ -58,14 +58,14 @@ class TestResult: def test_view_all(self,default): - a = default.view('increments',True).get('F') + a = default.view(increments=True).get('F') - assert dict_equal(a,default.view('increments','*').get('F')) - assert dict_equal(a,default.view('increments',default.increments_in_range(0,np.iinfo(int).max)).get('F')) + assert dict_equal(a,default.view(increments='*').get('F')) + assert dict_equal(a,default.view(increments=default.increments_in_range(0,np.iinfo(int).max)).get('F')) - assert dict_equal(a,default.view('times',True).get('F')) - assert dict_equal(a,default.view('times','*').get('F')) - assert dict_equal(a,default.view('times',default.times_in_range(0.0,np.inf)).get('F')) + assert dict_equal(a,default.view(times=True).get('F')) + assert dict_equal(a,default.view(times='*').get('F')) + assert dict_equal(a,default.view(times=default.times_in_range(0.0,np.inf)).get('F')) @pytest.mark.parametrize('what',['increments','times','phases','fields']) # ToDo: discuss homogenizations def test_view_none(self,default,what): @@ -314,7 +314,7 @@ class TestResult: @pytest.mark.parametrize('overwrite',['off','on']) def test_add_overwrite(self,default,overwrite): - last = default.view('increments',-1) + last = default.view(increments=-1) last.add_stress_Cauchy() @@ -377,7 +377,7 @@ class TestResult: @pytest.mark.parametrize('inc',[4,0],ids=range(2)) @pytest.mark.xfail(int(vtk.vtkVersion.GetVTKVersion().split('.')[0])<9, reason='missing "Direction" attribute') def test_vtk(self,request,tmp_path,ref_path,update,patch_execution_stamp,patch_datetime_now,output,fname,inc): - result = Result(ref_path/fname).view('increments',inc) + result = Result(ref_path/fname).view(increments=inc) os.chdir(tmp_path) result.export_VTK(output,parallel=False) fname = fname.split('.')[0]+f'_inc{(inc if type(inc) == int else inc[0]):0>2}.vti' @@ -400,7 +400,7 @@ class TestResult: result.export_VTK(output,mode) def test_marc_coordinates(self,ref_path): - result = Result(ref_path/'check_compile_job1.hdf5').view('increments',-1) + result = Result(ref_path/'check_compile_job1.hdf5').view(increments=-1) c_n = result.coordinates0_node + result.get('u_n') c_p = result.coordinates0_point + result.get('u_p') assert len(c_n) > len(c_p) @@ -440,7 +440,7 @@ class TestResult: dim_xdmf = reader_xdmf.GetOutput().GetDimensions() bounds_xdmf = reader_xdmf.GetOutput().GetBounds() - single_phase.view('increments',0).export_VTK(parallel=False) + single_phase.view(increments=0).export_VTK(parallel=False) fname = os.path.splitext(os.path.basename(single_phase.fname))[0]+'_inc00.vti' reader_vti = vtk.vtkXMLImageDataReader() reader_vti.SetFileName(fname) From c1eabc33846df7a1e77e366fe443b1df29576cb3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 17 Dec 2021 10:31:41 +0100 Subject: [PATCH 35/74] simplified interface --- python/damask/_result.py | 58 +++++++++++++++---------------------- python/tests/test_Result.py | 10 +++---- 2 files changed, 28 insertions(+), 40 deletions(-) diff --git a/python/damask/_result.py b/python/damask/_result.py index 981b0dcb3..2453001a6 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -246,36 +246,6 @@ class Result: return dup - def modification_enable(self): - """ - Allow modification of existing data. - - Returns - ------- - modified_view : damask.Result - View without write-protection of existing data. - - """ - print(util.warn('Warning: Modification of existing datasets allowed!')) - dup = self.copy() - dup._allow_modification = True - return dup - - def modification_disable(self): - """ - Prevent modification of existing data (default case). - - Returns - ------- - modified_view : damask.Result - View with write-protection of existing data. - - """ - dup = self.copy() - dup._allow_modification = False - return dup - - def increments_in_range(self,start,end): """ Get all increments within a given range. @@ -329,7 +299,8 @@ class Result: times=None, phases=None, homogenizations=None, - fields=None): + fields=None, + protected=None): """ Set view. @@ -353,6 +324,8 @@ class Result: Name(s) of homogenizations to select. fields: (list of) str, or bool, optional. Name(s) of fields to select. + protected: bool, optional. + Protection status of existing data. Returns ------- @@ -374,8 +347,23 @@ class Result: >>> r_t10to40 = r.view(times=r.times_in_range(10.0,40.0)) """ - what_, datasets_ = _view_transition(what,datasets,increments,times,phases,homogenizations,fields) - return self._manage_view('set',what_,datasets_) + v = _view_transition(what,datasets,increments,times,phases,homogenizations,fields) + if protected is not None: + if v == None: + dup = self.copy() + else: + what_,datasets_ = v + dup = self._manage_view('set',what_,datasets_) + else: + what_,datasets_ = v + dup = self._manage_view('set',what_,datasets_) + + if protected is False: + dup._allow_modification = True + if protected is True: + print(util.warn('Warning: Modification of existing datasets allowed!')) + dup._allow_modification = False + return dup def view_more(self,what=None,datasets=None,*, @@ -495,7 +483,7 @@ class Result: >>> import damask >>> r = damask.Result('my_file.hdf5') - >>> r_unprotected = r.modification_enable() + >>> r_unprotected = r.view(protected=False) >>> r_unprotected.rename('F','def_grad') """ @@ -534,7 +522,7 @@ class Result: >>> import damask >>> r = damask.Result('my_file.hdf5') - >>> r_unprotected = r.modification_enable() + >>> r_unprotected = r.view(protected=False) >>> r_unprotected.remove('F') """ diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index 3bc363ebc..4fea08eb3 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -322,9 +322,9 @@ class TestResult: created_first = datetime.strptime(created_first,'%Y-%m-%d %H:%M:%S%z') if overwrite == 'on': - last = last.modification_enable() + last = last.view(protected=False) else: - last = last.modification_disable() + last = last.view(protected=True) time.sleep(2.) try: @@ -344,10 +344,10 @@ class TestResult: def test_rename(self,default,allowed): if allowed == 'on': F = default.place('F') - default = default.modification_enable() + default = default.view(protected=False) default.rename('F','new_name') assert np.all(F == default.place('new_name')) - default = default.modification_disable() + default = default.view(protected=True) with pytest.raises(PermissionError): default.rename('P','another_new_name') @@ -355,7 +355,7 @@ class TestResult: @pytest.mark.parametrize('allowed',['off','on']) def test_remove(self,default,allowed): if allowed == 'on': - unsafe = default.modification_enable() + unsafe = default.view(protected=False) unsafe.remove('F') assert unsafe.get('F') is None else: From a1e42af8602a63f7c5d8fbaa1b71cdd5028f117a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 17 Dec 2021 10:37:45 +0100 Subject: [PATCH 36/74] easier to understand --- python/damask/_result.py | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/python/damask/_result.py b/python/damask/_result.py index 2453001a6..019f9c87d 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -156,7 +156,7 @@ class Result: self.fname = Path(fname).absolute() - self._allow_modification = False + self._protected = True def __copy__(self): @@ -349,20 +349,18 @@ class Result: """ v = _view_transition(what,datasets,increments,times,phases,homogenizations,fields) if protected is not None: - if v == None: + if v is None: dup = self.copy() else: what_,datasets_ = v dup = self._manage_view('set',what_,datasets_) + if protected is True: + print(util.warn('Warning: Modification of existing datasets allowed!')) + dup._protected = protected else: what_,datasets_ = v dup = self._manage_view('set',what_,datasets_) - if protected is False: - dup._allow_modification = True - if protected is True: - print(util.warn('Warning: Modification of existing datasets allowed!')) - dup._allow_modification = False return dup @@ -487,7 +485,7 @@ class Result: >>> r_unprotected.rename('F','def_grad') """ - if not self._allow_modification: + if self._protected: raise PermissionError('Renaming datasets not permitted') with h5py.File(self.fname,'a') as f: @@ -526,7 +524,7 @@ class Result: >>> r_unprotected.remove('F') """ - if not self._allow_modification: + if self._protected: raise PermissionError('Removing datasets not permitted') with h5py.File(self.fname,'a') as f: @@ -1417,7 +1415,7 @@ class Result: lock.acquire() with h5py.File(self.fname, 'a') as f: try: - if self._allow_modification and '/'.join([group,result['label']]) in f: + if not self._protected and '/'.join([group,result['label']]) in f: dataset = f['/'.join([group,result['label']])] dataset[...] = result['data'] dataset.attrs['overwritten'] = True From 2fc25c6fa9d81783a15b2ab243e706954f3259ca Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 17 Dec 2021 14:28:12 +0100 Subject: [PATCH 37/74] correct reporting --- python/damask/_result.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/_result.py b/python/damask/_result.py index 019f9c87d..f47a7da80 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -354,7 +354,7 @@ class Result: else: what_,datasets_ = v dup = self._manage_view('set',what_,datasets_) - if protected is True: + if not protected: print(util.warn('Warning: Modification of existing datasets allowed!')) dup._protected = protected else: From 0468bfd3e1e09d9f51afb552cd808791b990bf65 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 18 Dec 2021 17:21:46 +0000 Subject: [PATCH 38/74] Use ArrayLike for numpy >= 1.20 --- python/damask/_colormap.py | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/python/damask/_colormap.py b/python/damask/_colormap.py index 5a13f440d..ba7617fbb 100644 --- a/python/damask/_colormap.py +++ b/python/damask/_colormap.py @@ -6,6 +6,10 @@ from pathlib import Path from typing import Sequence, Union, TextIO import numpy as np +try: + from numpy.typing import ArrayLike +except ImportError: + ArrayLike = Union[np.ndarray,Sequence[float]] # type: ignore import scipy.interpolate as interp import matplotlib as mpl if os.name == 'posix' and 'DISPLAY' not in os.environ: @@ -78,8 +82,8 @@ class Colormap(mpl.colors.ListedColormap): @staticmethod - def from_range(low: Sequence[float], - high: Sequence[float], + def from_range(low: ArrayLike, + high: ArrayLike, name: str = 'DAMASK colormap', N: int = 256, model: str = 'rgb') -> 'Colormap': @@ -129,7 +133,7 @@ class Colormap(mpl.colors.ListedColormap): if model.lower() not in toMsh: raise ValueError(f'Invalid color model: {model}.') - low_high = np.vstack((low,high)) + low_high = np.vstack((low,high)).astype(float) out_of_bounds = np.bool_(False) if model.lower() == 'rgb': @@ -142,7 +146,7 @@ class Colormap(mpl.colors.ListedColormap): out_of_bounds = np.any(low_high[:,0]<0) if out_of_bounds: - raise ValueError(f'{model.upper()} colors {low} | {high} are out of bounds.') + raise ValueError(f'{model.upper()} colors {low_high[0]} | {low_high[1]} are out of bounds.') low_,high_ = map(toMsh[model.lower()],low_high) msh = map(functools.partial(Colormap._interpolate_msh,low=low_,high=high_),np.linspace(0,1,N)) @@ -225,7 +229,7 @@ class Colormap(mpl.colors.ListedColormap): def shade(self, field: np.ndarray, - bounds: Sequence[float] = None, + bounds: ArrayLike = None, gap: float = None) -> Image: """ Generate PIL image of 2D field using colormap. @@ -235,7 +239,7 @@ class Colormap(mpl.colors.ListedColormap): field : numpy.array, shape (:,:) Data to be shaded. bounds : sequence of float, len (2), optional - Value range (low,high) spanned by colormap. + Value range (left,right) spanned by colormap. gap : field.dtype, optional Transparent value. NaN will always be rendered transparent. @@ -248,17 +252,17 @@ class Colormap(mpl.colors.ListedColormap): mask = np.logical_not(np.isnan(field) if gap is None else \ np.logical_or (np.isnan(field), field == gap)) # mask NaN (and gap if present) - lo,hi = (field[mask].min(),field[mask].max()) if bounds is None else \ - (min(bounds[:2]),max(bounds[:2])) + l,r = (field[mask].min(),field[mask].max()) if bounds is None else \ + np.array(bounds,float)[:2] - delta,avg = hi-lo,0.5*(hi+lo) + delta,avg = r-l,0.5*abs(r+l) - if delta * 1e8 <= avg: # delta is similar to numerical noise - hi,lo = hi+0.5*avg,lo-0.5*avg # extend range to have actual data centered within + if abs(delta) * 1e8 <= avg: # delta is similar to numerical noise + l,r = l-0.5*avg*np.sign(delta),r+0.5*avg*np.sign(delta), # extend range to have actual data centered within return Image.fromarray( (np.dstack(( - self.colors[(np.round(np.clip((field-lo)/(hi-lo),0.0,1.0)*(self.N-1))).astype(np.uint16),:3], + self.colors[(np.round(np.clip((field-l)/delta,0.0,1.0)*(self.N-1))).astype(np.uint16),:3], mask.astype(float) ) )*255 From 28449393bfb851fe6155c2e9fa2d005566b533ae Mon Sep 17 00:00:00 2001 From: Test User Date: Sat, 18 Dec 2021 22:44:25 +0100 Subject: [PATCH 39/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-297-g5ecfba1e5 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 7ba8f4262..42ae7a5db 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-283-gdacd08f39 +v3.0.0-alpha5-297-g5ecfba1e5 From d060eabb459844f2e10056a21c020c722efc737a Mon Sep 17 00:00:00 2001 From: Test User Date: Sun, 19 Dec 2021 00:24:12 +0100 Subject: [PATCH 40/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-299-g29122ef41 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 7ba8f4262..3485580fd 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-283-gdacd08f39 +v3.0.0-alpha5-299-g29122ef41 From 5af6cc288b6e892f594c611689d1fc42d5f52181 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 19 Dec 2021 21:46:10 +0100 Subject: [PATCH 41/74] whitespace adjustments --- src/phase_mechanical_plastic_nonlocal.f90 | 313 +++++++++++----------- 1 file changed, 158 insertions(+), 155 deletions(-) diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index bf0dc3b19..7ff016514 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -772,66 +772,67 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & tau, & !< resolved shear stress including backstress terms dot_gamma !< shear rate + associate(prm => param(ph),dst=>dependentState(ph),stt=>state(ph)) - !*** shortcut to state variables - rho = getRho(ph,en) - rhoSgl = rho(:,sgl) + !*** shortcut to state variables + rho = getRho(ph,en) + rhoSgl = rho(:,sgl) - do s = 1,prm%sum_N_sl - tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) - tauNS(s,1) = tau(s) - tauNS(s,2) = tau(s) - if (tau(s) > 0.0_pReal) then - tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_pos(1:3,1:3,s)) - tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_neg(1:3,1:3,s)) - else - tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_neg(1:3,1:3,s)) - tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_pos(1:3,1:3,s)) - end if - end do - tauNS = tauNS + spread(dst%tau_back(:,en),2,4) - tau = tau + dst%tau_back(:,en) - - ! edges - call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), & - tau, tauNS(:,1), dst%tau_pass(:,en),1,Temperature, ph) - v(:,2) = v(:,1) - dv_dtau(:,2) = dv_dtau(:,1) - dv_dtauNS(:,2) = dv_dtauNS(:,1) - - !screws - if (prm%nonSchmidActive) then - do t = 3,4 - call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), & - tau, tauNS(:,t), dst%tau_pass(:,en),2,Temperature, ph) + do s = 1,prm%sum_N_sl + tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) + tauNS(s,1) = tau(s) + tauNS(s,2) = tau(s) + if (tau(s) > 0.0_pReal) then + tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_pos(1:3,1:3,s)) + tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_neg(1:3,1:3,s)) + else + tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_neg(1:3,1:3,s)) + tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_pos(1:3,1:3,s)) + end if end do - else - v(:,3:4) = spread(v(:,1),2,2) - dv_dtau(:,3:4) = spread(dv_dtau(:,1),2,2) - dv_dtauNS(:,3:4) = spread(dv_dtauNS(:,1),2,2) - end if + tauNS = tauNS + spread(dst%tau_back(:,en),2,4) + tau = tau + dst%tau_back(:,en) - stt%v(:,en) = pack(v,.true.) + ! edges + call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), & + tau, tauNS(:,1), dst%tau_pass(:,en),1,Temperature, ph) + v(:,2) = v(:,1) + dv_dtau(:,2) = dv_dtau(:,1) + dv_dtauNS(:,2) = dv_dtauNS(:,1) - !*** Bauschinger effect - forall (s = 1:prm%sum_N_sl, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) & - rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t)) + !screws + if (prm%nonSchmidActive) then + do t = 3,4 + call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), & + tau, tauNS(:,t), dst%tau_pass(:,en),2,Temperature, ph) + end do + else + v(:,3:4) = spread(v(:,1),2,2) + dv_dtau(:,3:4) = spread(dv_dtau(:,1),2,2) + dv_dtauNS(:,3:4) = spread(dv_dtauNS(:,1),2,2) + end if - dot_gamma = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl + stt%v(:,en) = pack(v,.true.) - Lp = 0.0_pReal - dLp_dMp = 0.0_pReal - do s = 1,prm%sum_N_sl - Lp = Lp + dot_gamma(s) * prm%P_sl(1:3,1:3,s) - forall (i=1:3,j=1:3,k=1:3,l=1:3) & - dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) & - + prm%P_sl(i,j,s) * prm%P_sl(k,l,s) & - * sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%b_sl(s) & - + prm%P_sl(i,j,s) & - * (+ prm%P_nS_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) & - - prm%P_nS_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s) - end do + !*** Bauschinger effect + forall (s = 1:prm%sum_N_sl, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) & + rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t)) + + dot_gamma = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl + + Lp = 0.0_pReal + dLp_dMp = 0.0_pReal + do s = 1,prm%sum_N_sl + Lp = Lp + dot_gamma(s) * prm%P_sl(1:3,1:3,s) + forall (i=1:3,j=1:3,k=1:3,l=1:3) & + dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) & + + prm%P_sl(i,j,s) * prm%P_sl(k,l,s) & + * sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%b_sl(s) & + + prm%P_sl(i,j,s) & + * (+ prm%P_nS_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) & + - prm%P_nS_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s) + end do end associate @@ -870,7 +871,8 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en) dUpper, & ! current maximum stable dipole distance for edges and screws dUpperOld, & ! old maximum stable dipole distance for edges and screws deltaDUpper ! change in maximum stable dipole distance for edges and screws - + + associate(prm => param(ph),dst => dependentState(ph),del => deltaState(ph)) mu = elastic_mu(ph,en) @@ -1394,78 +1396,79 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) belowThreshold type(rotation) :: mis + associate(prm => param(ph)) - ns = prm%sum_N_sl + ns = prm%sum_N_sl - en = material_phaseMemberAt(1,i,e) - !*** start out fully compatible - my_compatibility = 0.0_pReal - forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pReal + en = material_phaseMemberAt(1,i,e) + !*** start out fully compatible + my_compatibility = 0.0_pReal + forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pReal - neighbors: do n = 1,nIPneighbors - neighbor_e = IPneighborhood(1,n,i,e) - neighbor_i = IPneighborhood(2,n,i,e) - neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e) - neighbor_phase = material_phaseAt(1,neighbor_e) + neighbors: do n = 1,nIPneighbors + neighbor_e = IPneighborhood(1,n,i,e) + neighbor_i = IPneighborhood(2,n,i,e) + neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e) + neighbor_phase = material_phaseAt(1,neighbor_e) - if (neighbor_e <= 0 .or. neighbor_i <= 0) then - !* FREE SURFACE - forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_surface) - elseif (neighbor_phase /= ph) then - !* PHASE BOUNDARY - if (plasticState(neighbor_phase)%nonlocal .and. plasticState(ph)%nonlocal) & - forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_pReal - elseif (prm%chi_GB >= 0.0_pReal) then - !* GRAIN BOUNDARY - if (any(dNeq(phase_O_0(ph)%data(en)%asQuaternion(), & - phase_O_0(neighbor_phase)%data(neighbor_me)%asQuaternion())) .and. & - plasticState(neighbor_phase)%nonlocal) & - forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB) - else - !* GRAIN BOUNDARY ? - !* Compatibility defined by relative orientation of slip systems: - !* The my_compatibility value is defined as the product of the slip normal projection and the slip direction projection. - !* Its sign is always positive for screws, for edges it has the same sign as the slip normal projection. - !* Since the sum for each slip system can easily exceed one (which would result in a transmissivity larger than one), - !* only values above or equal to a certain threshold value are considered. This threshold value is chosen, such that - !* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one. - !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one. - !* All values below the threshold are set to zero. - mis = orientation(ph)%data(en)%misorientation(orientation(neighbor_phase)%data(neighbor_me)) - mySlipSystems: do s1 = 1,ns - neighborSlipSystems: do s2 = 1,ns - my_compatibility(1,s2,s1,n) = math_inner(prm%slip_normal(1:3,s1), & - mis%rotate(prm%slip_normal(1:3,s2))) & - * abs(math_inner(prm%slip_direction(1:3,s1), & - mis%rotate(prm%slip_direction(1:3,s2)))) - my_compatibility(2,s2,s1,n) = abs(math_inner(prm%slip_normal(1:3,s1), & - mis%rotate(prm%slip_normal(1:3,s2)))) & - * abs(math_inner(prm%slip_direction(1:3,s1), & - mis%rotate(prm%slip_direction(1:3,s2)))) - end do neighborSlipSystems + if (neighbor_e <= 0 .or. neighbor_i <= 0) then + !* FREE SURFACE + forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_surface) + elseif (neighbor_phase /= ph) then + !* PHASE BOUNDARY + if (plasticState(neighbor_phase)%nonlocal .and. plasticState(ph)%nonlocal) & + forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_pReal + elseif (prm%chi_GB >= 0.0_pReal) then + !* GRAIN BOUNDARY + if (any(dNeq(phase_O_0(ph)%data(en)%asQuaternion(), & + phase_O_0(neighbor_phase)%data(neighbor_me)%asQuaternion())) .and. & + plasticState(neighbor_phase)%nonlocal) & + forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB) + else + !* GRAIN BOUNDARY ? + !* Compatibility defined by relative orientation of slip systems: + !* The my_compatibility value is defined as the product of the slip normal projection and the slip direction projection. + !* Its sign is always positive for screws, for edges it has the same sign as the slip normal projection. + !* Since the sum for each slip system can easily exceed one (which would result in a transmissivity larger than one), + !* only values above or equal to a certain threshold value are considered. This threshold value is chosen, such that + !* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one. + !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one. + !* All values below the threshold are set to zero. + mis = orientation(ph)%data(en)%misorientation(orientation(neighbor_phase)%data(neighbor_me)) + mySlipSystems: do s1 = 1,ns + neighborSlipSystems: do s2 = 1,ns + my_compatibility(1,s2,s1,n) = math_inner(prm%slip_normal(1:3,s1), & + mis%rotate(prm%slip_normal(1:3,s2))) & + * abs(math_inner(prm%slip_direction(1:3,s1), & + mis%rotate(prm%slip_direction(1:3,s2)))) + my_compatibility(2,s2,s1,n) = abs(math_inner(prm%slip_normal(1:3,s1), & + mis%rotate(prm%slip_normal(1:3,s2)))) & + * abs(math_inner(prm%slip_direction(1:3,s1), & + mis%rotate(prm%slip_direction(1:3,s2)))) + end do neighborSlipSystems - my_compatibilitySum = 0.0_pReal - belowThreshold = .true. - do while (my_compatibilitySum < 1.0_pReal .and. any(belowThreshold)) - thresholdValue = maxval(my_compatibility(2,:,s1,n), belowThreshold) ! screws always positive - nThresholdValues = real(count(my_compatibility(2,:,s1,n) >= thresholdValue),pReal) - where (my_compatibility(2,:,s1,n) >= thresholdValue) belowThreshold = .false. - if (my_compatibilitySum + thresholdValue * nThresholdValues > 1.0_pReal) & - where (abs(my_compatibility(:,:,s1,n)) >= thresholdValue) & - my_compatibility(:,:,s1,n) = sign((1.0_pReal - my_compatibilitySum)/nThresholdValues,& - my_compatibility(:,:,s1,n)) - my_compatibilitySum = my_compatibilitySum + nThresholdValues * thresholdValue - end do + my_compatibilitySum = 0.0_pReal + belowThreshold = .true. + do while (my_compatibilitySum < 1.0_pReal .and. any(belowThreshold)) + thresholdValue = maxval(my_compatibility(2,:,s1,n), belowThreshold) ! screws always positive + nThresholdValues = real(count(my_compatibility(2,:,s1,n) >= thresholdValue),pReal) + where (my_compatibility(2,:,s1,n) >= thresholdValue) belowThreshold = .false. + if (my_compatibilitySum + thresholdValue * nThresholdValues > 1.0_pReal) & + where (abs(my_compatibility(:,:,s1,n)) >= thresholdValue) & + my_compatibility(:,:,s1,n) = sign((1.0_pReal - my_compatibilitySum)/nThresholdValues,& + my_compatibility(:,:,s1,n)) + my_compatibilitySum = my_compatibilitySum + nThresholdValues * thresholdValue + end do - where(belowThreshold) my_compatibility(1,:,s1,n) = 0.0_pReal - where(belowThreshold) my_compatibility(2,:,s1,n) = 0.0_pReal + where(belowThreshold) my_compatibility(1,:,s1,n) = 0.0_pReal + where(belowThreshold) my_compatibility(2,:,s1,n) = 0.0_pReal - end do mySlipSystems - end if + end do mySlipSystems + end if - end do neighbors + end do neighbors - compatibility(:,:,:,:,i,e) = my_compatibility + compatibility(:,:,:,:,i,e) = my_compatibility end associate @@ -1647,52 +1650,52 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T, criticalStress_S !< maximum obstacle strength associate(prm => param(ph)) - v = 0.0_pReal - dv_dtau = 0.0_pReal - dv_dtauNS = 0.0_pReal + v = 0.0_pReal + dv_dtau = 0.0_pReal + dv_dtauNS = 0.0_pReal - do s = 1,prm%sum_N_sl - if (abs(tau(s)) > tauThreshold(s)) then + do s = 1,prm%sum_N_sl + if (abs(tau(s)) > tauThreshold(s)) then - !* Peierls contribution - tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s)) - lambda_P = prm%b_sl(s) - activationVolume_P = prm%w *prm%b_sl(s)**3 - criticalStress_P = prm%peierlsStress(s,c) - activationEnergy_P = criticalStress_P * activationVolume_P - tauRel_P = min(1.0_pReal, tauEff / criticalStress_P) - tPeierls = 1.0_pReal / prm%nu_a & - * exp(activationEnergy_P / (K_B * T) & - * (1.0_pReal - tauRel_P**prm%p)**prm%q) - dtPeierls_dtau = merge(tPeierls * prm%p * prm%q * activationVolume_P / (K_B * T) & - * (1.0_pReal - tauRel_P**prm%p)**(prm%q-1.0_pReal) * tauRel_P**(prm%p-1.0_pReal), & - 0.0_pReal, & - tauEff < criticalStress_P) + !* Peierls contribution + tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s)) + lambda_P = prm%b_sl(s) + activationVolume_P = prm%w *prm%b_sl(s)**3 + criticalStress_P = prm%peierlsStress(s,c) + activationEnergy_P = criticalStress_P * activationVolume_P + tauRel_P = min(1.0_pReal, tauEff / criticalStress_P) + tPeierls = 1.0_pReal / prm%nu_a & + * exp(activationEnergy_P / (K_B * T) & + * (1.0_pReal - tauRel_P**prm%p)**prm%q) + dtPeierls_dtau = merge(tPeierls * prm%p * prm%q * activationVolume_P / (K_B * T) & + * (1.0_pReal - tauRel_P**prm%p)**(prm%q-1.0_pReal) * tauRel_P**(prm%p-1.0_pReal), & + 0.0_pReal, & + tauEff < criticalStress_P) - ! Contribution from solid solution strengthening - tauEff = abs(tau(s)) - tauThreshold(s) - lambda_S = prm%b_sl(s) / sqrt(prm%c_sol) - activationVolume_S = prm%f_sol * prm%b_sl(s)**3 / sqrt(prm%c_sol) - criticalStress_S = prm%Q_sol / activationVolume_S - tauRel_S = min(1.0_pReal, tauEff / criticalStress_S) - tSolidSolution = 1.0_pReal / prm%nu_a & - * exp(prm%Q_sol / (K_B * T)* (1.0_pReal - tauRel_S**prm%p)**prm%q) - dtSolidSolution_dtau = merge(tSolidSolution * prm%p * prm%q * activationVolume_S / (K_B * T) & - * (1.0_pReal - tauRel_S**prm%p)**(prm%q-1.0_pReal)* tauRel_S**(prm%p-1.0_pReal), & - 0.0_pReal, & - tauEff < criticalStress_S) + ! Contribution from solid solution strengthening + tauEff = abs(tau(s)) - tauThreshold(s) + lambda_S = prm%b_sl(s) / sqrt(prm%c_sol) + activationVolume_S = prm%f_sol * prm%b_sl(s)**3 / sqrt(prm%c_sol) + criticalStress_S = prm%Q_sol / activationVolume_S + tauRel_S = min(1.0_pReal, tauEff / criticalStress_S) + tSolidSolution = 1.0_pReal / prm%nu_a & + * exp(prm%Q_sol / (K_B * T)* (1.0_pReal - tauRel_S**prm%p)**prm%q) + dtSolidSolution_dtau = merge(tSolidSolution * prm%p * prm%q * activationVolume_S / (K_B * T) & + * (1.0_pReal - tauRel_S**prm%p)**(prm%q-1.0_pReal)* tauRel_S**(prm%p-1.0_pReal), & + 0.0_pReal, & + tauEff < criticalStress_S) - !* viscous glide velocity - tauEff = abs(tau(s)) - tauThreshold(s) + !* viscous glide velocity + tauEff = abs(tau(s)) - tauThreshold(s) - v(s) = sign(1.0_pReal,tau(s)) & - / (tPeierls / lambda_P + tSolidSolution / lambda_S + prm%B /(prm%b_sl(s) * tauEff)) - dv_dtau(s) = v(s)**2 * (dtSolidSolution_dtau / lambda_S + prm%B / (prm%b_sl(s) * tauEff**2)) - dv_dtauNS(s) = v(s)**2 * dtPeierls_dtau / lambda_P + v(s) = sign(1.0_pReal,tau(s)) & + / (tPeierls / lambda_P + tSolidSolution / lambda_S + prm%B /(prm%b_sl(s) * tauEff)) + dv_dtau(s) = v(s)**2 * (dtSolidSolution_dtau / lambda_S + prm%B / (prm%b_sl(s) * tauEff**2)) + dv_dtauNS(s) = v(s)**2 * dtPeierls_dtau / lambda_P - end if - end do + end if + end do end associate From 00230d482f30ce8155852d1e201613bbf6b37a06 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 19 Dec 2021 22:07:23 +0100 Subject: [PATCH 42/74] use data from other physics directly more clear code, simplified interfaces --- src/phase.f90 | 2 +- src/phase_mechanical_plastic.f90 | 22 +--- ...phase_mechanical_plastic_dislotungsten.f90 | 8 +- src/phase_mechanical_plastic_dislotwin.f90 | 112 +++++++++--------- src/phase_mechanical_plastic_nonlocal.f90 | 21 ++-- src/phase_thermal.f90 | 2 +- 6 files changed, 81 insertions(+), 86 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index 214b0f5fa..6035b4491 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -155,7 +155,7 @@ module phase real(pReal), dimension(3,3) :: P end function phase_P - module function thermal_T(ph,en) result(T) + pure module function thermal_T(ph,en) result(T) integer, intent(in) :: ph,en real(pReal) :: T end function thermal_T diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 4951712fd..0a24b0cbf 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -73,47 +73,37 @@ submodule(phase:mechanical) plastic en end subroutine kinehardening_LpAndItsTangent - module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) + module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) 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 - real(pReal), intent(in) :: & - T integer, intent(in) :: & ph, & en end subroutine dislotwin_LpAndItsTangent - pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) + pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) 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 - real(pReal), intent(in) :: & - T integer, intent(in) :: & ph, & en end subroutine dislotungsten_LpAndItsTangent - module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & - Mp,Temperature,ph,en) + module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) 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 !< Mandel stress - real(pReal), intent(in) :: & - Temperature integer, intent(in) :: & ph, & en @@ -282,13 +272,13 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) case (PLASTIC_NONLOCAL_ID) plasticType - call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) + call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) case (PLASTIC_DISLOTWIN_ID) plasticType - call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) + call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) case (PLASTIC_DISLOTUNGSTEN_ID) plasticType - call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) + call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) end select plasticType diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index c2a584ea1..cd71a7fd7 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -257,27 +257,27 @@ 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,ph,en) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress - real(pReal), intent(in) :: & - T !< temperature integer, intent(in) :: & ph, & en integer :: & i,k,l,m,n + real(pReal) :: & + T !< temperature real(pReal), dimension(param(ph)%sum_N_sl) :: & dot_gamma_pos,dot_gamma_neg, & ddot_gamma_dtau_pos,ddot_gamma_dtau_neg + T = thermal_T(ph,en) Lp = 0.0_pReal dLp_dMp = 0.0_pReal diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index c5e043e3f..cfc9a002f 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -513,20 +513,20 @@ end function plastic_dislotwin_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- -module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) +module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) 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, intent(in) :: ph,en - real(pReal), intent(in) :: T integer :: i,k,l,m,n real(pReal) :: & - f_unrotated,StressRatio_p,& - E_kB_T, & - ddot_gamma_dtau, & - tau + f_unrotated,StressRatio_p,& + E_kB_T, & + ddot_gamma_dtau, & + tau, & + T real(pReal), dimension(param(ph)%sum_N_sl) :: & dot_gamma_sl,ddot_gamma_dtau_sl real(pReal), dimension(param(ph)%sum_N_tw) :: & @@ -556,69 +556,71 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) 0, 1, 1 & ],pReal),[ 3,6]) - associate(prm => param(ph), stt => state(ph)) - - f_unrotated = 1.0_pReal & - - sum(stt%f_tw(1:prm%sum_N_tw,en)) & - - sum(stt%f_tr(1:prm%sum_N_tr,en)) + T = thermal_T(ph,en) Lp = 0.0_pReal dLp_dMp = 0.0_pReal - call kinetics_sl(Mp,T,ph,en,dot_gamma_sl,ddot_gamma_dtau_sl) - slipContribution: do i = 1, prm%sum_N_sl - Lp = Lp + dot_gamma_sl(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_sl(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) - end do slipContribution + associate(prm => param(ph), stt => state(ph)) - call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw) - twinContibution: do i = 1, prm%sum_N_tw - Lp = Lp + dot_gamma_tw(i)*prm%P_tw(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_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) - end do twinContibution + f_unrotated = 1.0_pReal & + - sum(stt%f_tw(1:prm%sum_N_tw,en)) & + - sum(stt%f_tr(1:prm%sum_N_tr,en)) - call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr) - transContibution: do i = 1, prm%sum_N_tr - Lp = Lp + dot_gamma_tr(i)*prm%P_tr(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_tr(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) - end do transContibution + call kinetics_sl(Mp,T,ph,en,dot_gamma_sl,ddot_gamma_dtau_sl) + slipContribution: do i = 1, prm%sum_N_sl + Lp = Lp + dot_gamma_sl(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_sl(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) + end do slipContribution - Lp = Lp * f_unrotated - dLp_dMp = dLp_dMp * f_unrotated + call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw) + twinContibution: do i = 1, prm%sum_N_tw + Lp = Lp + dot_gamma_tw(i)*prm%P_tw(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_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) + end do twinContibution - shearBandingContribution: if (dNeq0(prm%v_sb)) then + call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr) + transContibution: do i = 1, prm%sum_N_tr + Lp = Lp + dot_gamma_tr(i)*prm%P_tr(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_tr(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) + end do transContibution - E_kB_T = prm%E_sb/(K_B*T) - call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? + Lp = Lp * f_unrotated + dLp_dMp = dLp_dMp * f_unrotated - do i = 1,6 - P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),& - matmul(eigVectors,sb_mComposition(1:3,i))) - tau = math_tensordot(Mp,P_sb) + shearBandingContribution: if (dNeq0(prm%v_sb)) then - significantShearBandStress: if (abs(tau) > tol_math_check) then - StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb - dot_gamma_sb = sign(prm%v_sb*exp(-E_kB_T*(1-StressRatio_p)**prm%q_sb), tau) - ddot_gamma_dtau = abs(dot_gamma_sb)*E_kB_T*prm%p_sb*prm%q_sb/prm%xi_sb & - * (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) & - * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal) + E_kB_T = prm%E_sb/(K_B*T) + call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? - Lp = Lp + dot_gamma_sb * P_sb - 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 * P_sb(k,l) * P_sb(m,n) - end if significantShearBandStress - end do + do i = 1,6 + P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),& + matmul(eigVectors,sb_mComposition(1:3,i))) + tau = math_tensordot(Mp,P_sb) - end if shearBandingContribution + significantShearBandStress: if (abs(tau) > tol_math_check) then + StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb + dot_gamma_sb = sign(prm%v_sb*exp(-E_kB_T*(1-StressRatio_p)**prm%q_sb), tau) + ddot_gamma_dtau = abs(dot_gamma_sb)*E_kB_T*prm%p_sb*prm%q_sb/prm%xi_sb & + * (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) & + * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal) - end associate + Lp = Lp + dot_gamma_sb * P_sb + 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 * P_sb(k,l) * P_sb(m,n) + end if significantShearBandStress + end do + + end if shearBandingContribution + + end associate end subroutine dislotwin_LpAndItsTangent diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 7ff016514..7118055fe 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -741,7 +741,7 @@ end subroutine nonlocal_dependentState !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & - Mp,Temperature,ph,en) + Mp,ph,en) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -749,9 +749,6 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & integer, intent(in) :: & ph, & en - real(pReal), intent(in) :: & - Temperature !< temperature - real(pReal), dimension(3,3), intent(in) :: & Mp !< derivative of Lp with respect to Mp @@ -771,8 +768,14 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & real(pReal), dimension(param(ph)%sum_N_sl) :: & tau, & !< resolved shear stress including backstress terms dot_gamma !< shear rate + real(pReal) :: & + Temperature !< temperature + Temperature = thermal_T(ph,en) + Lp = 0.0_pReal + dLp_dMp = 0.0_pReal + associate(prm => param(ph),dst=>dependentState(ph),stt=>state(ph)) !*** shortcut to state variables @@ -821,8 +824,6 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & dot_gamma = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl - Lp = 0.0_pReal - dLp_dMp = 0.0_pReal do s = 1,prm%sum_N_sl Lp = Lp + dot_gamma(s) * prm%P_sl(1:3,1:3,s) forall (i=1:3,j=1:3,k=1:3,l=1:3) & @@ -1649,10 +1650,12 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T, criticalStress_P, & !< maximum obstacle strength criticalStress_S !< maximum obstacle strength + + v = 0.0_pReal + dv_dtau = 0.0_pReal + dv_dtauNS = 0.0_pReal + associate(prm => param(ph)) - v = 0.0_pReal - dv_dtau = 0.0_pReal - dv_dtauNS = 0.0_pReal do s = 1,prm%sum_N_sl if (abs(tau(s)) > tauThreshold(s)) then diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 4ea21d039..a3e4dd628 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -271,7 +271,7 @@ end subroutine thermal_forward !---------------------------------------------------------------------------------------------- !< @brief Get temperature (for use by non-thermal physics) !---------------------------------------------------------------------------------------------- -module function thermal_T(ph,en) result(T) +pure module function thermal_T(ph,en) result(T) integer, intent(in) :: ph, en real(pReal) :: T From e10dea5b6cefeaa668a3e6f608cd1b64306c5391 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 19 Dec 2021 22:53:48 +0100 Subject: [PATCH 43/74] easier to understand --- src/phase_mechanical_plastic_dislotwin.f90 | 34 +++++++++++----------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index cfc9a002f..78681fa24 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -476,18 +476,18 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC) C66_tw, & C66_tr integer :: i - real(pReal) :: f_unrotated + real(pReal) :: f_matrix C = elastic_C66(ph,en) associate(prm => param(ph), stt => state(ph)) - f_unrotated = 1.0_pReal & - - sum(stt%f_tw(1:prm%sum_N_tw,en)) & - - sum(stt%f_tr(1:prm%sum_N_tr,en)) + f_matrix = 1.0_pReal & + - sum(stt%f_tw(1:prm%sum_N_tw,en)) & + - sum(stt%f_tr(1:prm%sum_N_tr,en)) - homogenizedC = f_unrotated * C + homogenizedC = f_matrix * C twinActive: if (prm%sum_N_tw > 0) then C66_tw = lattice_C66_twin(prm%N_tw,C,phase_lattice(ph),phase_cOverA(ph)) @@ -522,7 +522,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) integer :: i,k,l,m,n real(pReal) :: & - f_unrotated,StressRatio_p,& + f_matrix,StressRatio_p,& E_kB_T, & ddot_gamma_dtau, & tau, & @@ -563,9 +563,9 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) associate(prm => param(ph), stt => state(ph)) - f_unrotated = 1.0_pReal & - - sum(stt%f_tw(1:prm%sum_N_tw,en)) & - - sum(stt%f_tr(1:prm%sum_N_tr,en)) + f_matrix = 1.0_pReal & + - sum(stt%f_tw(1:prm%sum_N_tw,en)) & + - sum(stt%f_tr(1:prm%sum_N_tr,en)) call kinetics_sl(Mp,T,ph,en,dot_gamma_sl,ddot_gamma_dtau_sl) slipContribution: do i = 1, prm%sum_N_sl @@ -591,8 +591,8 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) + ddot_gamma_dtau_tr(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) end do transContibution - Lp = Lp * f_unrotated - dLp_dMp = dLp_dMp * f_unrotated + Lp = Lp * f_matrix + dLp_dMp = dLp_dMp * f_matrix shearBandingContribution: if (dNeq0(prm%v_sb)) then @@ -640,7 +640,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) integer :: i real(pReal) :: & - f_unrotated, & + f_matrix, & d_hat, & v_cl, & !< climb velocity tau, & @@ -663,9 +663,9 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) mu = elastic_mu(ph,en) nu = elastic_nu(ph,en) - f_unrotated = 1.0_pReal & - - sum(stt%f_tw(1:prm%sum_N_tw,en)) & - - sum(stt%f_tr(1:prm%sum_N_tr,en)) + f_matrix = 1.0_pReal & + - sum(stt%f_tw(1:prm%sum_N_tw,en)) & + - sum(stt%f_tr(1:prm%sum_N_tr,en)) call kinetics_sl(Mp,T,ph,en,dot_gamma_sl) dot%gamma_sl(:,en) = abs(dot_gamma_sl) @@ -711,10 +711,10 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) - dot_rho_dip_climb call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw) - dot%f_tw(:,en) = f_unrotated*dot_gamma_tw/prm%gamma_char + dot%f_tw(:,en) = f_matrix*dot_gamma_tw/prm%gamma_char call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr) - dot%f_tr(:,en) = f_unrotated*dot_gamma_tr + dot%f_tr(:,en) = f_matrix*dot_gamma_tr end associate From 32e741ae2cc0f2e9c144ace41789d998d75aedba Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Tue, 21 Dec 2021 15:54:46 +0100 Subject: [PATCH 44/74] elastic constants for bcc iron --- .../phase/mechanical/elastic/Hooke_Fe.yaml | 24 +++++++++++++------ 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/examples/config/phase/mechanical/elastic/Hooke_Fe.yaml b/examples/config/phase/mechanical/elastic/Hooke_Fe.yaml index e736829a6..94e761aaf 100644 --- a/examples/config/phase/mechanical/elastic/Hooke_Fe.yaml +++ b/examples/config/phase/mechanical/elastic/Hooke_Fe.yaml @@ -1,9 +1,19 @@ type: Hooke references: - - J.P. Hirth and J. Lothe, - Theory of Dislocations, 1982, - John Wiley & Sons, - page 837 -C_11: 242.e9 -C_12: 146.5e+9 -C_44: 112.e9 + - D.J. Dever, + Journal of Applied Physics 43(8):3293-3301, 1972, + https://doi.org/10.1063/1.1661710 + +T_ref: 300 + +C_11: 231.7e+9 +C_11,T: -47.6e+6 +C_11,T^2: -54.4e+3 + +C_12: 135.8e+9 +C_12,T: -12.9e+6 +C_12,T^2: -7.3e+3 + +C_44: 116.8e+9 +C_44,T: -19.4e+6 +C_44,T^2: -2.5e+3 From 3e01f72daa09e48d9cf6402216ea40296ae35013 Mon Sep 17 00:00:00 2001 From: Test User Date: Tue, 21 Dec 2021 22:19:29 +0100 Subject: [PATCH 45/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-308-gb4f3ac457 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 42ae7a5db..290dfd544 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-297-g5ecfba1e5 +v3.0.0-alpha5-308-gb4f3ac457 From 508083082d35c4dc11bd5421f50e52066f88aa90 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 21 Dec 2021 23:50:43 +0100 Subject: [PATCH 46/74] simplified --- python/damask/_colormap.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/python/damask/_colormap.py b/python/damask/_colormap.py index ba7617fbb..c2721e0fa 100644 --- a/python/damask/_colormap.py +++ b/python/damask/_colormap.py @@ -236,7 +236,7 @@ class Colormap(mpl.colors.ListedColormap): Parameters ---------- - field : numpy.array, shape (:,:) + field : numpy.ndarray, shape (:,:) Data to be shaded. bounds : sequence of float, len (2), optional Value range (left,right) spanned by colormap. @@ -616,7 +616,7 @@ class Colormap(mpl.colors.ListedColormap): @staticmethod - def _lab2xyz(lab: np.ndarray, ref_white: np.ndarray = None) -> np.ndarray: + def _lab2xyz(lab: np.ndarray, ref_white: np.ndarray = _REF_WHITE) -> np.ndarray: """ CIE Lab to CIE Xyz. @@ -624,6 +624,8 @@ class Colormap(mpl.colors.ListedColormap): ---------- lab : numpy.ndarray, shape (3) CIE lab values. + ref_white : numpy.ndarray, shape (3) + Reference white, default value is the standard 2° observer for D65. Returns ------- @@ -642,10 +644,10 @@ class Colormap(mpl.colors.ListedColormap): f_x**3. if f_x**3. > _EPS else (116.*f_x-16.)/_KAPPA, ((lab[0]+16.)/116.)**3 if lab[0]>_KAPPA*_EPS else lab[0]/_KAPPA, f_z**3. if f_z**3. > _EPS else (116.*f_z-16.)/_KAPPA - ])*(ref_white if ref_white is not None else _REF_WHITE) + ])*ref_white @staticmethod - def _xyz2lab(xyz: np.ndarray, ref_white: np.ndarray = None) -> np.ndarray: + def _xyz2lab(xyz: np.ndarray, ref_white: np.ndarray = _REF_WHITE) -> np.ndarray: """ CIE Xyz to CIE Lab. @@ -653,6 +655,8 @@ class Colormap(mpl.colors.ListedColormap): ---------- xyz : numpy.ndarray, shape (3) CIE Xyz values. + ref_white : numpy.ndarray, shape (3) + Reference white, default value is the standard 2° observer for D65. Returns ------- @@ -664,7 +668,6 @@ class Colormap(mpl.colors.ListedColormap): http://www.brucelindbloom.com/index.html?Eqn_Lab_to_XYZ.html """ - ref_white = ref_white if ref_white is not None else _REF_WHITE f = np.where(xyz/ref_white > _EPS,(xyz/ref_white)**(1./3.),(_KAPPA*xyz/ref_white+16.)/116.) return np.array([ From dc5172ce6351450f5d220ad65a946875ca921ebc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 22 Dec 2021 12:41:16 +0100 Subject: [PATCH 47/74] don't use deprecated functions --- python/tests/test_VTK.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/tests/test_VTK.py b/python/tests/test_VTK.py index 26d1c4a53..68224ff33 100644 --- a/python/tests/test_VTK.py +++ b/python/tests/test_VTK.py @@ -20,7 +20,7 @@ def default(): """Simple VTK.""" cells = np.array([5,6,7],int) size = np.array([.6,1.,.5]) - return VTK.from_rectilinear_grid(cells,size) + return VTK.from_image_data(cells,size) class TestVTK: @@ -116,7 +116,7 @@ class TestVTK: def test_add_extension(self,tmp_path,default): default.save(tmp_path/'default.txt',parallel=False) - assert os.path.isfile(tmp_path/'default.txt.vtr') + assert os.path.isfile(tmp_path/'default.txt.vti') def test_invalid_get(self,default): @@ -160,7 +160,7 @@ class TestVTK: def test_comments(self,tmp_path,default): default.add_comments(['this is a comment']) default.save(tmp_path/'with_comments',parallel=False) - new = VTK.load(tmp_path/'with_comments.vtr') + new = VTK.load(tmp_path/'with_comments.vti') assert new.get_comments() == ['this is a comment'] @pytest.mark.xfail(int(vtk.vtkVersion.GetVTKVersion().split('.')[0])<8, reason='missing METADATA') From 6e833d2327b1bc4ea21cd95add449ccdde89fb28 Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 22 Dec 2021 12:56:07 +0100 Subject: [PATCH 48/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-311-g38d497819 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 290dfd544..cd791b3f1 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-308-gb4f3ac457 +v3.0.0-alpha5-311-g38d497819 From 8c2266bc50bf57c2ca406dd5ce5f9926de829bf4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 23 Dec 2021 19:23:02 +0100 Subject: [PATCH 49/74] bugfix/workaround for failing pipelines on GitHub credits to @miketimofeev and @mikhailkoliada for the fix --- .github/workflows/Python.yml | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/.github/workflows/Python.yml b/.github/workflows/Python.yml index c4fb15bdc..4c199b81c 100644 --- a/.github/workflows/Python.yml +++ b/.github/workflows/Python.yml @@ -43,12 +43,13 @@ jobs: pip install pytest - name: Install dependencies + # https://github.com/actions/virtual-environments/issues/4790 run: > sudo apt-get update && - sudo apt-get install python3-pip python3-pytest python3-pandas python3-scipy - python3-h5py python3-vtk7 python3-matplotlib python3-yaml -y + sudo apt-get remove mysql* && + sudo apt-get install python3-pandas python3-scipy python3-h5py python3-vtk7 python3-matplotlib python3-yaml -y - name: Run unit tests run: | export PYTHONPATH=${PWD}/python - COLUMNS=256 python -m pytest python + COLUMNS=256 pytest python From 368f4ac9854c22e292754f9f4bc73ebbcc4e077b Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 23 Dec 2021 22:35:44 +0100 Subject: [PATCH 50/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-313-g8c2266bc5 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index cd791b3f1..852dd348b 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-311-g38d497819 +v3.0.0-alpha5-313-g8c2266bc5 From 68b34f6022a4bb5751a654cf6fbd69f3ce12cf52 Mon Sep 17 00:00:00 2001 From: Test User Date: Sat, 25 Dec 2021 07:15:40 +0100 Subject: [PATCH 51/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-316-gd8ef23fc9 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 852dd348b..4815547f4 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-313-g8c2266bc5 +v3.0.0-alpha5-316-gd8ef23fc9 From ac7caa0a643337403f45bf212442d404129f759c Mon Sep 17 00:00:00 2001 From: Test User Date: Sat, 25 Dec 2021 08:39:13 +0100 Subject: [PATCH 52/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-323-g2ea293063 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 852dd348b..56ff18f14 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-313-g8c2266bc5 +v3.0.0-alpha5-323-g2ea293063 From e34b8a632c616531b2ba84c53a8681c2ef18af6d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 25 Dec 2021 13:11:45 +0100 Subject: [PATCH 53/74] new Intel compilers are available --- .gitlab-ci.yml | 6 +++--- cmake/Compiler-Intel.cmake | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 8ceb7575a..d60921b8c 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -36,13 +36,13 @@ variables: # Names of module files to load # =============================================================================================== # ++++++++++++ Compiler +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - COMPILER_INTEL: "Compiler/Intel/19.1.2 Libraries/IMKL/2020" + COMPILER_INTEL: "Compiler/Intel/2022.0.1 Libraries/IMKL/2022.0.1" COMPILER_GNU: "Compiler/GNU/10" # ++++++++++++ MPI ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - MPI_INTEL: "MPI/Intel/19.1.2/IntelMPI/2019" + MPI_INTEL: "MPI/Intel/2022.0.1/IntelMPI/2021.5.0" MPI_GNU: "MPI/GNU/10/OpenMPI/4.1.1" # ++++++++++++ PETSc ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - PETSC_INTEL: "Libraries/PETSc/3.16.1/Intel-19.1.2-IntelMPI-2019" + PETSC_INTEL: "Libraries/PETSc/3.16.2/Intel-2022.0.1-IntelMPI-2021.5.0" PETSC_GNU: "Libraries/PETSc/3.16.1/GNU-10-OpenMPI-4.1.1" # ++++++++++++ MSC Marc +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ MSC: "FEM/MSC/2021.3.1" diff --git a/cmake/Compiler-Intel.cmake b/cmake/Compiler-Intel.cmake index 5b551069e..b964a31e6 100644 --- a/cmake/Compiler-Intel.cmake +++ b/cmake/Compiler-Intel.cmake @@ -106,7 +106,7 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -fpe-all=0") #set (DEBUG_FLAGS "${DEBUG_FLAGS},stderrors") # ... warnings about Fortran standard violations are changed to errors -set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug-parameters all") +#set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug-parameters all") # generate debug information for parameters # Additional options From a9293db1714d5c412848cf7ef6d67270eb503cc4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 25 Dec 2021 18:33:54 +0100 Subject: [PATCH 54/74] new LLVM-based Intel compilers work quick test shows significant performance improvement --- CMakeLists.txt | 2 + cmake/Compiler-Intel.cmake | 3 +- cmake/Compiler-IntelLLVM.cmake | 121 +++++++++++++++++++++++++++++++++ 3 files changed, 125 insertions(+), 1 deletion(-) create mode 100644 cmake/Compiler-IntelLLVM.cmake diff --git a/CMakeLists.txt b/CMakeLists.txt index 56416d0fd..e8b774cfc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -82,6 +82,8 @@ if (CMAKE_Fortran_COMPILER_ID STREQUAL "Intel") include(Compiler-Intel) elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") include(Compiler-GNU) +elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "IntelLLVM") + include(Compiler-IntelLLVM) else() message(FATAL_ERROR "Compiler type(CMAKE_Fortran_COMPILER_ID) not recognized") endif() diff --git a/cmake/Compiler-Intel.cmake b/cmake/Compiler-Intel.cmake index b964a31e6..3afffd2be 100644 --- a/cmake/Compiler-Intel.cmake +++ b/cmake/Compiler-Intel.cmake @@ -106,8 +106,9 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -fpe-all=0") #set (DEBUG_FLAGS "${DEBUG_FLAGS},stderrors") # ... warnings about Fortran standard violations are changed to errors -#set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug-parameters all") +#set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug-parameters all") # generate debug information for parameters +# Disabled due to ICE when compiling phase_damage.f90 (not understandable, there is no parameter in there) # Additional options # -heap-arrays: Should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Probably it helps also to unlimit other limits diff --git a/cmake/Compiler-IntelLLVM.cmake b/cmake/Compiler-IntelLLVM.cmake new file mode 100644 index 000000000..5cd6a0e63 --- /dev/null +++ b/cmake/Compiler-IntelLLVM.cmake @@ -0,0 +1,121 @@ +################################################################################################### +# Intel Compiler +################################################################################################### +if (CMAKE_Fortran_COMPILER_VERSION VERSION_LESS 18.0) + message (FATAL_ERROR "Intel Compiler version: ${CMAKE_Fortran_COMPILER_VERSION} not supported") +endif () + +if (OPENMP) + set (OPENMP_FLAGS "-qopenmp -parallel") +endif () + +if (OPTIMIZATION STREQUAL "OFF") + set (OPTIMIZATION_FLAGS "-O0 -no-ip") +elseif (OPTIMIZATION STREQUAL "DEFENSIVE") + set (OPTIMIZATION_FLAGS "-O2") +elseif (OPTIMIZATION STREQUAL "AGGRESSIVE") + set (OPTIMIZATION_FLAGS "-ipo -O3 -no-prec-div -fp-model fast=2 -xHost") + # -fast = -ipo, -O3, -no-prec-div, -static, -fp-model fast=2, and -xHost" +endif () + +# -assume std_mod_proc_name (included in -standard-semantics) causes problems if other modules +# (PETSc, HDF5) are not compiled with this option (https://software.intel.com/en-us/forums/intel-fortran-compiler-for-linux-and-mac-os-x/topic/62172) +set (STANDARD_CHECK "-stand f18 -assume nostd_mod_proc_name") +set (LINKER_FLAGS "${LINKER_FLAGS} -shared-intel") +# Link against shared Intel libraries instead of static ones + +#------------------------------------------------------------------------------------------------ +# Fine tuning compilation options +set (COMPILE_FLAGS "${COMPILE_FLAGS} -fpp") +# preprocessor + +set (COMPILE_FLAGS "${COMPILE_FLAGS} -ftz") +# flush underflow to zero, automatically set if -O[1,2,3] + +set (COMPILE_FLAGS "${COMPILE_FLAGS} -diag-disable") +# disables warnings ... +set (COMPILE_FLAGS "${COMPILE_FLAGS} 5268") +# ... the text exceeds right hand column allowed on the line (we have only comments there) +set (COMPILE_FLAGS "${COMPILE_FLAGS},7624") +# ... about deprecated forall (has nice syntax and most likely a performance advantage) + +set (COMPILE_FLAGS "${COMPILE_FLAGS} -warn") +# enables warnings ... +set (COMPILE_FLAGS "${COMPILE_FLAGS} declarations") +# ... any undeclared names (alternative name: -implicitnone) +set (COMPILE_FLAGS "${COMPILE_FLAGS},general") +# ... warning messages and informational messages are issued by the compiler +set (COMPILE_FLAGS "${COMPILE_FLAGS},usage") +# ... questionable programming practices +set (COMPILE_FLAGS "${COMPILE_FLAGS},interfaces") +# ... checks the interfaces of all SUBROUTINEs called and FUNCTIONs invoked in your compilation against an external set of interface blocks +set (COMPILE_FLAGS "${COMPILE_FLAGS},ignore_loc") +# ... %LOC is stripped from an actual argument +set (COMPILE_FLAGS "${COMPILE_FLAGS},alignments") +# ... data that is not naturally aligned +set (COMPILE_FLAGS "${COMPILE_FLAGS},unused") +# ... declared variables that are never used + +# Additional options +# -warn: enables warnings, where +# truncated_source: Determines whether warnings occur when source exceeds the maximum column width in fixed-format files. +# (too many warnings because we have comments beyond character 132) +# uncalled: Determines whether warnings occur when a statement function is never called +# all: +# -name as_is: case sensitive Fortran! + +#------------------------------------------------------------------------------------------------ +# Runtime debugging +set (DEBUG_FLAGS "${DEBUG_FLAGS} -g") +# Generate symbolic debugging information in the object file + +set (DEBUG_FLAGS "${DEBUG_FLAGS} -traceback") +# Generate extra information in the object file to provide source file traceback information when a severe error occurs at run time + +set (DEBUG_FLAGS "${DEBUG_FLAGS} -gen-interfaces") +# Generate an interface block for each routine. http://software.intel.com/en-us/blogs/2012/01/05/doctor-fortran-gets-explicit-again/ + +set (DEBUG_FLAGS "${DEBUG_FLAGS} -fp-stack-check") +# Generate extra code after every function call to ensure that the floating-point (FP) stack is in the expected state + +set (DEBUG_FLAGS "${DEBUG_FLAGS} -fp-model strict") +# Trap uninitalized variables + +set (DEBUG_FLAGS "${DEBUG_FLAGS} -check" ) +# Checks at runtime ... +set (DEBUG_FLAGS "${DEBUG_FLAGS} bounds") +# ... if an array index is too small (<1) or too large! +set (DEBUG_FLAGS "${DEBUG_FLAGS},format") +# ... for the data type of an item being formatted for output. +set (DEBUG_FLAGS "${DEBUG_FLAGS},output_conversion") +# ... for the fit of data items within a designated format descriptor field. +set (DEBUG_FLAGS "${DEBUG_FLAGS},pointers") +# ... for certain disassociated or uninitialized pointers or unallocated allocatable objects. +set (DEBUG_FLAGS "${DEBUG_FLAGS},uninit") +# ... for uninitialized variables. +set (DEBUG_FLAGS "${DEBUG_FLAGS} -ftrapuv") +# ... initializes stack local variables to an unusual value to aid error detection +set (DEBUG_FLAGS "${DEBUG_FLAGS} -fpe-all=0") +# ... capture all floating-point exceptions, sets -ftz automatically + +# disable due to compiler bug https://community.intel.com/t5/Intel-Fortran-Compiler/false-positive-stand-f18-and-IEEE-SELECTED-REAL-KIND/m-p/1227336 +#set (DEBUG_FLAGS "${DEBUG_FLAGS} -warn") +# enables warnings ... +#set (DEBUG_FLAGS "${DEBUG_FLAGS} errors") +# ... warnings are changed to errors +#set (DEBUG_FLAGS "${DEBUG_FLAGS},stderrors") +# ... warnings about Fortran standard violations are changed to errors + +set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug-parameters all") +# generate debug information for parameters + +# Additional options +# -heap-arrays: Should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Probably it helps also to unlimit other limits +# -check: Checks at runtime, where +# arg_temp_created: will cause a lot of warnings because we create a bunch of temporary arrays (performance?) +# stack: + +#------------------------------------------------------------------------------------------------ +# precision settings +set (PRECISION_FLAGS "${PRECISION_FLAGS} -real-size 64") +# set precision for standard real to 32 | 64 | 128 (= 4 | 8 | 16 bytes, type pReal is always 8 bytes) From a7440da26a8e72927147107ac5f4a907a92294cc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 26 Dec 2021 13:06:41 +0100 Subject: [PATCH 55/74] options not supported by LLVM-based Intel compilers --- cmake/Compiler-IntelLLVM.cmake | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cmake/Compiler-IntelLLVM.cmake b/cmake/Compiler-IntelLLVM.cmake index 5cd6a0e63..326cfe319 100644 --- a/cmake/Compiler-IntelLLVM.cmake +++ b/cmake/Compiler-IntelLLVM.cmake @@ -6,15 +6,15 @@ if (CMAKE_Fortran_COMPILER_VERSION VERSION_LESS 18.0) endif () if (OPENMP) - set (OPENMP_FLAGS "-qopenmp -parallel") + set (OPENMP_FLAGS "-qopenmp") endif () if (OPTIMIZATION STREQUAL "OFF") - set (OPTIMIZATION_FLAGS "-O0 -no-ip") + set (OPTIMIZATION_FLAGS "-O0") elseif (OPTIMIZATION STREQUAL "DEFENSIVE") set (OPTIMIZATION_FLAGS "-O2") elseif (OPTIMIZATION STREQUAL "AGGRESSIVE") - set (OPTIMIZATION_FLAGS "-ipo -O3 -no-prec-div -fp-model fast=2 -xHost") + set (OPTIMIZATION_FLAGS "-ipo -O3 -fp-model fast=2 -xHost") # -fast = -ipo, -O3, -no-prec-div, -static, -fp-model fast=2, and -xHost" endif () From 4168b101f4f775c4392758c870d70adacc647d1d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 26 Dec 2021 13:07:58 +0100 Subject: [PATCH 56/74] testing LLVM-based compilers does not work for grid due to missing ML --- .gitlab-ci.yml | 44 +++++++++++++++++++++++++++----------------- 1 file changed, 27 insertions(+), 17 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index d60921b8c..6f18e3f1e 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -36,14 +36,17 @@ variables: # Names of module files to load # =============================================================================================== # ++++++++++++ Compiler +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + COMPILER_GNU: "Compiler/GNU/10" + COMPILER_INTELLLVM: "Compiler/oneAPI/2022.0.1 Libraries/IMKL/2022.0.1" COMPILER_INTEL: "Compiler/Intel/2022.0.1 Libraries/IMKL/2022.0.1" - COMPILER_GNU: "Compiler/GNU/10" # ++++++++++++ MPI ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - MPI_INTEL: "MPI/Intel/2022.0.1/IntelMPI/2021.5.0" MPI_GNU: "MPI/GNU/10/OpenMPI/4.1.1" + MPI_INTELLLVM: "MPI/oneAPI/2022.0.1/IntelMPI/2021.5.0" + MPI_INTEL: "MPI/Intel/2022.0.1/IntelMPI/2021.5.0" # ++++++++++++ PETSc ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + PETSC_GNU: "Libraries/PETSc/3.16.1/GNU-10-OpenMPI-4.1.1" + PETSC_INTELLLVM: "Libraries/PETSc/3.16.2/oneAPI-2022.0.1-IntelMPI-2021.5.0" PETSC_INTEL: "Libraries/PETSc/3.16.2/Intel-2022.0.1-IntelMPI-2021.5.0" - PETSC_GNU: "Libraries/PETSc/3.16.1/GNU-10-OpenMPI-4.1.1" # ++++++++++++ MSC Marc +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ MSC: "FEM/MSC/2021.3.1" IntelMarc: "Compiler/Intel/19.1.2 Libraries/IMKL/2020" @@ -76,20 +79,6 @@ mypy: ################################################################################################### -test_grid_Intel: - stage: compile - script: - - module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL} - - cd PRIVATE/testing/pytest - - pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_Intel - -test_mesh_Intel: - stage: compile - script: - - module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL} - - cd PRIVATE/testing/pytest - - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_Intel - test_grid_GNU: stage: compile script: @@ -104,6 +93,27 @@ test_mesh_GNU: - cd PRIVATE/testing/pytest - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_GNU +test_mesh_IntelLLVM: + stage: compile + script: + - module load ${COMPILER_INTELLLVM} ${MPI_INTELLLVM} ${PETSC_INTELLLVM} + - cd PRIVATE/testing/pytest + - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_IntelLLVM + +test_grid_Intel: + stage: compile + script: + - module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL} + - cd PRIVATE/testing/pytest + - pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_Intel + +test_mesh_Intel: + stage: compile + script: + - module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL} + - cd PRIVATE/testing/pytest + - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_Intel + test_Marc: stage: compile script: From 2d347e1be18c16aefc92043d7a2c37a240137a2b Mon Sep 17 00:00:00 2001 From: Test User Date: Sun, 26 Dec 2021 21:05:58 +0100 Subject: [PATCH 57/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-333-g01cd92755 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 4815547f4..e657e6a62 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-316-gd8ef23fc9 +v3.0.0-alpha5-333-g01cd92755 From 7e4a374c9e047d030a5f4211f9aa11ee7ba00fc3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 27 Dec 2021 09:20:15 +0100 Subject: [PATCH 58/74] new material parameters - austenitic steel AISI304 - SAE1050 after quenching to martensite I don't know how the old set of elastic constants for martenstite (originally coming from 10.1088/0965-0393/23/4/045005) was derived, but it is hard to believe that the elastic behavior can be that strongly modified by a heat treatment. --- PRIVATE | 2 +- .../Phase_Phenopowerlaw_BCC-Martensite.yaml | 23 ++++++++----------- examples/config/phase/AISI304.yaml | 6 +++++ .../eigen/thermalexpansion_AISI304.yaml | 7 ++++++ .../mechanical/elastic/Hooke_AISI304.yaml | 8 +++++++ .../elastic/Hooke_SAE1050-martensite.yaml | 8 +++++++ examples/config/phase/thermal/AISI304.yaml | 9 ++++++++ .../{Steel-0.5C.yaml => steel-0.5C.yaml} | 0 8 files changed, 48 insertions(+), 15 deletions(-) create mode 100644 examples/config/phase/AISI304.yaml create mode 100644 examples/config/phase/mechanical/eigen/thermalexpansion_AISI304.yaml create mode 100644 examples/config/phase/mechanical/elastic/Hooke_AISI304.yaml create mode 100644 examples/config/phase/mechanical/elastic/Hooke_SAE1050-martensite.yaml create mode 100644 examples/config/phase/thermal/AISI304.yaml rename examples/config/phase/thermal/{Steel-0.5C.yaml => steel-0.5C.yaml} (100%) diff --git a/PRIVATE b/PRIVATE index 96c32ba42..b898a8b55 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 96c32ba4237a51eaad92cd139e1a716ee5b32493 +Subproject commit b898a8b5552bd9d1c555edc3d8134564dd32fe53 diff --git a/examples/config/Phase_Phenopowerlaw_BCC-Martensite.yaml b/examples/config/Phase_Phenopowerlaw_BCC-Martensite.yaml index 3ad6779eb..f2066b373 100644 --- a/examples/config/Phase_Phenopowerlaw_BCC-Martensite.yaml +++ b/examples/config/Phase_Phenopowerlaw_BCC-Martensite.yaml @@ -1,17 +1,12 @@ # Tasan et.al. 2015 Acta Materalia # Tasan et.al. 2015 International Journal of Plasticity # Diehl et.al. 2015 Meccanica -Martensite: - lattice: cI - mechanical: - elastic: {C_11: 417.4e+9, C_12: 242.4e+9, C_44: 211.1e+9, type: Hooke} - plastic: - N_sl: [12, 12] - a_sl: 2.0 - dot_gamma_0_sl: 0.001 - h_0_sl-sl: 563.0e+9 - h_sl-sl: [1, 1.4, 1, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4] - n_sl: 20 - type: phenopowerlaw - xi_0_sl: [405.8e+6, 456.7e+6] - xi_inf_sl: [872.9e+6, 971.2e+6] +N_sl: [12, 12] +a_sl: 2.0 +dot_gamma_0_sl: 0.001 +h_0_sl-sl: 563.0e+9 +h_sl-sl: [1, 1.4, 1, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4] +n_sl: 20 +type: phenopowerlaw +xi_0_sl: [405.8e+6, 456.7e+6] +xi_inf_sl: [872.9e+6, 971.2e+6] diff --git a/examples/config/phase/AISI304.yaml b/examples/config/phase/AISI304.yaml new file mode 100644 index 000000000..0ca55635b --- /dev/null +++ b/examples/config/phase/AISI304.yaml @@ -0,0 +1,6 @@ +references: + - H.M. Ledbetter + physica status solidi (a) 85(1):89-96, 1984 + https://doi.org/10.1002/pssa.2210850111 +lattice: cF +rho: 7937.0 diff --git a/examples/config/phase/mechanical/eigen/thermalexpansion_AISI304.yaml b/examples/config/phase/mechanical/eigen/thermalexpansion_AISI304.yaml new file mode 100644 index 000000000..6a3b339bb --- /dev/null +++ b/examples/config/phase/mechanical/eigen/thermalexpansion_AISI304.yaml @@ -0,0 +1,7 @@ +type: thermalexpansion +references: + - R.H. Bogaard et al. + Thermochimica Acta 218:373-393, 1993 + https://doi.org/10.1016/0040-6031(93)80437-F +A_11: 15.0e-6 +T_ref: 300.0 diff --git a/examples/config/phase/mechanical/elastic/Hooke_AISI304.yaml b/examples/config/phase/mechanical/elastic/Hooke_AISI304.yaml new file mode 100644 index 000000000..70fbd25c7 --- /dev/null +++ b/examples/config/phase/mechanical/elastic/Hooke_AISI304.yaml @@ -0,0 +1,8 @@ +type: Hooke +references: + - H.M. Ledbetter + physica status solidi (a) 85(1):89-96, 1984 + https://doi.org/10.1002/pssa.2210850111 +C_11: 204.6e+9 +C_12: 137.7e+9 +C_44: 126.2e+9 diff --git a/examples/config/phase/mechanical/elastic/Hooke_SAE1050-martensite.yaml b/examples/config/phase/mechanical/elastic/Hooke_SAE1050-martensite.yaml new file mode 100644 index 000000000..7ba941066 --- /dev/null +++ b/examples/config/phase/mechanical/elastic/Hooke_SAE1050-martensite.yaml @@ -0,0 +1,8 @@ +type: Hooke +references: + - S.A. Kim and W.L. Johnson, + Materials Science & Engineering A 452-453:633-639, 2007, + https://doi.org/10.1016/j.msea.2006.11.147 +C_11: 268.1e+9 +C_12: 111.2e+9 +C_44: 79.06e+9 diff --git a/examples/config/phase/thermal/AISI304.yaml b/examples/config/phase/thermal/AISI304.yaml new file mode 100644 index 000000000..ef09fd850 --- /dev/null +++ b/examples/config/phase/thermal/AISI304.yaml @@ -0,0 +1,9 @@ +references: + - B.F. Blackwell et al. + Proceedings of 34th National Heat Transfer Conference 2000 + https://www.osti.gov/servlets/purl/760791 + - R.H. Bogaard et al. + Thermochimica Acta 218:373-393, 1993 + https://doi.org/10.1016/0040-6031(93)80437-F +C_p: 470.0 +K_11: 14.34 diff --git a/examples/config/phase/thermal/Steel-0.5C.yaml b/examples/config/phase/thermal/steel-0.5C.yaml similarity index 100% rename from examples/config/phase/thermal/Steel-0.5C.yaml rename to examples/config/phase/thermal/steel-0.5C.yaml From b2d0fd8ff5c0c2675af9ce8c1fdd449aef223969 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 28 Dec 2021 10:19:17 +0000 Subject: [PATCH 59/74] Lambert azimuthal equal-area (laea) projection --- python/damask/util.py | 70 +++++++++++++++++++++++++++++++-------- python/tests/test_util.py | 19 +++++++++-- 2 files changed, 72 insertions(+), 17 deletions(-) diff --git a/python/damask/util.py b/python/damask/util.py index a39a7b8c7..0581302db 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -18,11 +18,11 @@ from . import version __all__=[ 'srepr', 'emph', 'deemph', 'warn', 'strikeout', - 'run', + 'run', 'natural_sort', 'show_progress', 'scale_to_coprime', - 'project_stereographic', + 'project_equal_angle', 'project_equal_area', 'hybrid_IA', 'execution_stamp', 'shapeshifter', 'shapeblender', @@ -267,13 +267,13 @@ def scale_to_coprime(v): return m -def project_stereographic(vector,direction='z',normalize=True,keepdims=False): +def project_equal_angle(vector,direction='z',normalize=True,keepdims=False): """ - Apply stereographic projection to vector. + Apply equal-angle projection to vector. Parameters ---------- - vector : numpy.ndarray of shape (...,3) + vector : numpy.ndarray, shape (...,3) Vector coordinates to be projected. direction : str Projection direction 'x', 'y', or 'z'. @@ -281,32 +281,74 @@ def project_stereographic(vector,direction='z',normalize=True,keepdims=False): normalize : bool Ensure unit length of input vector. Defaults to True. keepdims : bool - Maintain three-dimensional output coordinates. - Default two-dimensional output uses right-handed frame spanned by + Maintain three-dimensional output coordinates. Defaults to False. + Two-dimensional output uses right-handed frame spanned by the next and next-next axis relative to the projection direction, e.g. x-y when projecting along z and z-x when projecting along y. Returns ------- - coordinates : numpy.ndarray of shape (...,2 | 3) + coordinates : numpy.ndarray, shape (...,2 | 3) Projected coordinates. Examples -------- >>> import damask >>> import numpy as np - >>> project_stereographic(np.ones(3)) + >>> project_equal_angle(np.ones(3)) [0.3660254, 0.3660254] - >>> project_stereographic(np.ones(3),direction='x',normalize=False,keepdims=True) + >>> project_equal_angle(np.ones(3),direction='x',normalize=False,keepdims=True) [0, 0.5, 0.5] - >>> project_stereographic([0,1,1],direction='y',normalize=True,keepdims=False) + >>> project_equal_angle([0,1,1],direction='y',normalize=True,keepdims=False) [0.41421356, 0] """ shift = 'zyx'.index(direction) - v_ = np.roll(vector/np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector, - shift,axis=-1) - return np.roll(np.block([v_[...,:2]/(1+np.abs(v_[...,2:3])),np.zeros_like(v_[...,2:3])]), + v = np.roll(vector/np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector, + shift,axis=-1) + return np.roll(np.block([v[...,:2]/(1.0+np.abs(v[...,2:3])),np.zeros_like(v[...,2:3])]), + -shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2] + +def project_equal_area(vector,direction='z',normalize=True,keepdims=False): + """ + Apply equal-area projection to vector. + + Parameters + ---------- + vector : numpy.ndarray, shape (...,3) + Vector coordinates to be projected. + direction : str + Projection direction 'x', 'y', or 'z'. + Defaults to 'z'. + normalize : bool + Ensure unit length of input vector. Defaults to True. + keepdims : bool + Maintain three-dimensional output coordinates. Defaults to False. + Two-dimensional output uses right-handed frame spanned by + the next and next-next axis relative to the projection direction, + e.g. x-y when projecting along z and z-x when projecting along y. + + Returns + ------- + coordinates : numpy.ndarray, shape (...,2 | 3) + Projected coordinates. + + Examples + -------- + >>> import damask + >>> import numpy as np + >>> project_equal_area(np.ones(3)) + [0.45970084, 0.45970084] + >>> project_equal_area(np.ones(3),direction='x',normalize=False,keepdims=True) + [0.0, 0.70710678, 0.70710678] + >>> project_equal_area([0,1,1],direction='y',normalize=True,keepdims=False) + [0.5411961, 0.0] + + """ + shift = 'zyx'.index(direction) + v = np.roll(vector/np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector, + shift,axis=-1) + return np.roll(np.block([v[...,:2]/np.sqrt(1.0+np.abs(v[...,2:3])),np.zeros_like(v[...,2:3])]), -shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2] diff --git a/python/tests/test_util.py b/python/tests/test_util.py index 13731eda2..ee345605f 100644 --- a/python/tests/test_util.py +++ b/python/tests/test_util.py @@ -59,9 +59,22 @@ class TestUtil: ([1,1,0],'x',False,False,[0.5,0]), ([1,1,1],'y',True, True, [0.3660254, 0,0.3660254]), ]) - def test_project_stereographic(self,point,direction,normalize,keepdims,answer): - assert np.allclose(util.project_stereographic(np.array(point),direction=direction, - normalize=normalize,keepdims=keepdims),answer) + def test_project_equal_angle(self,point,direction,normalize,keepdims,answer): + assert np.allclose(util.project_equal_angle(np.array(point),direction=direction, + normalize=normalize,keepdims=keepdims),answer) + + @pytest.mark.parametrize('point,direction,normalize,keepdims,answer', + [ + ([1,0,0],'z',False,True, [1,0,0]), + ([1,0,0],'z',True, False,[1,0]), + ([0,1,1],'z',False,True, [0,0.70710678,0]), + ([0,1,1],'y',True, False,[0.5411961,0]), + ([1,1,0],'x',False,False,[0.70710678,0]), + ([1,1,1],'y',True, True, [0.45970084,0,0.45970084]), + ]) + def test_project_equal_area(self,point,direction,normalize,keepdims,answer): + assert np.allclose(util.project_equal_area(np.array(point),direction=direction, + normalize=normalize,keepdims=keepdims),answer) @pytest.mark.parametrize('fro,to,mode,answer', [ From e19c17c2bccba2287a7ad30a01605e95dcb58bb1 Mon Sep 17 00:00:00 2001 From: Test User Date: Tue, 28 Dec 2021 14:19:37 +0100 Subject: [PATCH 60/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-336-g6871eb302 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index e657e6a62..7a82d0bd4 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-333-g01cd92755 +v3.0.0-alpha5-336-g6871eb302 From 4583c17080e1ce956dc5778632b1ab09d7c22436 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 29 Dec 2021 07:00:28 +0100 Subject: [PATCH 61/74] corrent 'intent' specification - http://www.netlib.org/lapack/explore-html/d7/d3b/group__double_g_esolve_ga5ee879032a8365897c3ba91e3dc8d512.html - http://www.netlib.org/lapack/explore-html/dd/d9a/group__double_g_ecomputational_ga56d9c860ce4ce42ded7f914fdb0683ff.html --- src/LAPACK_interface.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/LAPACK_interface.f90 b/src/LAPACK_interface.f90 index 7d3043ed0..91d2451a4 100644 --- a/src/LAPACK_interface.f90 +++ b/src/LAPACK_interface.f90 @@ -23,7 +23,7 @@ module LAPACK_interface integer, intent(in) :: n,nrhs,lda,ldb real(pReal), intent(inout), dimension(lda,n) :: a integer, intent(out), dimension(n) :: ipiv - real(pReal), intent(out), dimension(ldb,nrhs) :: b + real(pReal), intent(inout), dimension(ldb,nrhs) :: b integer, intent(out) :: info end subroutine dgesv @@ -39,7 +39,7 @@ module LAPACK_interface use prec integer, intent(in) :: n,lda,lwork real(pReal), intent(inout), dimension(lda,n) :: a - integer, intent(out), dimension(n) :: ipiv + integer, intent(in), dimension(n) :: ipiv real(pReal), intent(out), dimension(max(1,lwork)) :: work integer, intent(out) :: info end subroutine dgetri From 59bb264b5fd9d3d7dc7d87b8aedba66727506d78 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 29 Dec 2021 07:09:52 +0100 Subject: [PATCH 62/74] LAPACK routines can be considered pure all arguments have 'intent' specification and don't access any global variables. output to screen only occurs in the case that someting goes wrong --- src/LAPACK_interface.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/LAPACK_interface.f90 b/src/LAPACK_interface.f90 index 91d2451a4..e11f74875 100644 --- a/src/LAPACK_interface.f90 +++ b/src/LAPACK_interface.f90 @@ -6,7 +6,7 @@ module LAPACK_interface interface - subroutine dgeev(jobvl,jobvr,n,a,lda,wr,wi,vl,ldvl,vr,ldvr,work,lwork,info) + pure subroutine dgeev(jobvl,jobvr,n,a,lda,wr,wi,vl,ldvl,vr,ldvr,work,lwork,info) use prec character, intent(in) :: jobvl,jobvr integer, intent(in) :: n,lda,ldvl,ldvr,lwork @@ -18,7 +18,7 @@ module LAPACK_interface integer, intent(out) :: info end subroutine dgeev - subroutine dgesv(n,nrhs,a,lda,ipiv,b,ldb,info) + pure subroutine dgesv(n,nrhs,a,lda,ipiv,b,ldb,info) use prec integer, intent(in) :: n,nrhs,lda,ldb real(pReal), intent(inout), dimension(lda,n) :: a @@ -27,7 +27,7 @@ module LAPACK_interface integer, intent(out) :: info end subroutine dgesv - subroutine dgetrf(m,n,a,lda,ipiv,info) + pure subroutine dgetrf(m,n,a,lda,ipiv,info) use prec integer, intent(in) :: m,n,lda real(pReal), intent(inout), dimension(lda,n) :: a @@ -35,7 +35,7 @@ module LAPACK_interface integer, intent(out) :: info end subroutine dgetrf - subroutine dgetri(n,a,lda,ipiv,work,lwork,info) + pure subroutine dgetri(n,a,lda,ipiv,work,lwork,info) use prec integer, intent(in) :: n,lda,lwork real(pReal), intent(inout), dimension(lda,n) :: a @@ -44,7 +44,7 @@ module LAPACK_interface integer, intent(out) :: info end subroutine dgetri - subroutine dsyev(jobz,uplo,n,a,lda,w,work,lwork,info) + pure subroutine dsyev(jobz,uplo,n,a,lda,w,work,lwork,info) use prec character, intent(in) :: jobz,uplo integer, intent(in) :: n,lda,lwork From fb51e3c4cd45372b394f5963c88063472a3a776a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 29 Dec 2021 07:19:26 +0100 Subject: [PATCH 63/74] functions have no side-effects, hence 'pure' --- src/lattice.f90 | 4 ++-- src/math.f90 | 12 ++++++------ 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 50d4b6c51..3089aa30b 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2070,7 +2070,7 @@ end function getlabels !> @brief Equivalent Poisson's ratio (ν) !> @details https://doi.org/10.1143/JPSJ.20.635 !-------------------------------------------------------------------------------------------------- -function lattice_equivalent_nu(C,assumption) result(nu) +pure function lattice_equivalent_nu(C,assumption) result(nu) real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation) character(len=5), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress) @@ -2103,7 +2103,7 @@ end function lattice_equivalent_nu !> @brief Equivalent shear modulus (μ) !> @details https://doi.org/10.1143/JPSJ.20.635 !-------------------------------------------------------------------------------------------------- -function lattice_equivalent_mu(C,assumption) result(mu) +pure function lattice_equivalent_mu(C,assumption) result(mu) real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation) character(len=5), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress) diff --git a/src/math.f90 b/src/math.f90 index 2d8f8fb3b..34cb6eecb 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -512,7 +512,7 @@ end subroutine math_invert33 !-------------------------------------------------------------------------------------------------- !> @brief Inversion of symmetriced 3x3x3x3 matrix !-------------------------------------------------------------------------------------------------- -function math_invSym3333(A) +pure function math_invSym3333(A) real(pReal),dimension(3,3,3,3) :: math_invSym3333 @@ -538,7 +538,7 @@ end function math_invSym3333 !-------------------------------------------------------------------------------------------------- !> @brief invert quadratic matrix of arbitrary dimension !-------------------------------------------------------------------------------------------------- -subroutine math_invert(InvA, error, A) +pure subroutine math_invert(InvA, error, A) real(pReal), dimension(:,:), intent(in) :: A real(pReal), dimension(size(A,1),size(A,1)), intent(out) :: invA @@ -996,7 +996,7 @@ end subroutine math_normal !-------------------------------------------------------------------------------------------------- !> @brief eigenvalues and eigenvectors of symmetric matrix !-------------------------------------------------------------------------------------------------- -subroutine math_eigh(w,v,error,m) +pure subroutine math_eigh(w,v,error,m) real(pReal), dimension(:,:), intent(in) :: m !< quadratic matrix to compute eigenvectors and values of real(pReal), dimension(size(m,1)), intent(out) :: w !< eigenvalues @@ -1021,7 +1021,7 @@ end subroutine math_eigh !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @details See http://arxiv.org/abs/physics/0610206 (DSYEVH3) !-------------------------------------------------------------------------------------------------- -subroutine math_eigh33(w,v,m) +pure subroutine math_eigh33(w,v,m) real(pReal), dimension(3,3),intent(in) :: m !< 3x3 matrix to compute eigenvectors and values of real(pReal), dimension(3), intent(out) :: w !< eigenvalues @@ -1114,7 +1114,7 @@ end function math_rotationalPart !> @brief Eigenvalues of symmetric matrix ! will return NaN on error !-------------------------------------------------------------------------------------------------- -function math_eigvalsh(m) +pure function math_eigvalsh(m) real(pReal), dimension(:,:), intent(in) :: m !< symmetric matrix to compute eigenvalues of real(pReal), dimension(size(m,1)) :: math_eigvalsh @@ -1137,7 +1137,7 @@ end function math_eigvalsh !> but apparently more stable solution and has general LAPACK powered version for arbritrary sized !> matrices as fallback !-------------------------------------------------------------------------------------------------- -function math_eigvalsh33(m) +pure function math_eigvalsh33(m) real(pReal), intent(in), dimension(3,3) :: m !< 3x3 symmetric matrix to compute eigenvalues of real(pReal), dimension(3) :: math_eigvalsh33,I From 2f74e0d0700b1f581226c17daceeb4d2a3082f42 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 30 Dec 2021 21:53:00 +0100 Subject: [PATCH 64/74] avoid failing self test increase number of samples to have less corner cases. Needs to be allocatable to avoid stack/heap issue on ifort --- src/math.f90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 34cb6eecb..3ea4d6a08 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1432,9 +1432,11 @@ subroutine selfTest error stop 'math_LeviCivita' normal_distribution: block - real(pReal), dimension(500000) :: r + integer, parameter :: N = 1000000 + real(pReal), dimension(:), allocatable :: r real(pReal) :: mu, sigma + allocate(r(N)) call random_number(mu) call random_number(sigma) @@ -1443,11 +1445,11 @@ subroutine selfTest call math_normal(r,mu,sigma) - if (abs(mu -sum(r)/real(size(r),pReal))>5.0e-2_pReal) & + if (abs(mu -sum(r)/real(N,pReal))>5.0e-2_pReal) & error stop 'math_normal(mu)' - mu = sum(r)/real(size(r),pReal) - if (abs(sigma**2 -1.0_pReal/real(size(r)-1,pReal) * sum((r-mu)**2))/sigma > 5.0e-2_pReal) & + mu = sum(r)/real(N,pReal) + if (abs(sigma**2 -1.0_pReal/real(N-1,pReal) * sum((r-mu)**2))/sigma > 5.0e-2_pReal) & error stop 'math_normal(sigma)' end block normal_distribution From b34655b7fc3cc2bc999af482a6e1abd828d7e812 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 30 Dec 2021 21:55:49 +0100 Subject: [PATCH 65/74] functions without side-effects are 'pure' basically all 'getter' functions should be pure --- src/phase_mechanical.f90 | 6 +++--- src/phase_mechanical_elastic.f90 | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 37e6fe7bf..1917d81c9 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -168,17 +168,17 @@ submodule(phase) mechanical integer, intent(in) :: ph,en end function plastic_dislotwin_homogenizedC - module function elastic_C66(ph,en) result(C66) + pure module function elastic_C66(ph,en) result(C66) real(pReal), dimension(6,6) :: C66 integer, intent(in) :: ph, en end function elastic_C66 - module function elastic_mu(ph,en) result(mu) + pure module function elastic_mu(ph,en) result(mu) real(pReal) :: mu integer, intent(in) :: ph, en end function elastic_mu - module function elastic_nu(ph,en) result(nu) + pure module function elastic_nu(ph,en) result(nu) real(pReal) :: nu integer, intent(in) :: ph, en end function elastic_nu diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index acbbac236..24797a880 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -86,7 +86,7 @@ end subroutine elastic_init !-------------------------------------------------------------------------------------------------- !> @brief return 6x6 elasticity tensor !-------------------------------------------------------------------------------------------------- -module function elastic_C66(ph,en) result(C66) +pure module function elastic_C66(ph,en) result(C66) integer, intent(in) :: & ph, & @@ -140,7 +140,7 @@ end function elastic_C66 !-------------------------------------------------------------------------------------------------- !> @brief return shear modulus !-------------------------------------------------------------------------------------------------- -module function elastic_mu(ph,en) result(mu) +pure module function elastic_mu(ph,en) result(mu) integer, intent(in) :: & ph, & @@ -157,7 +157,7 @@ end function elastic_mu !-------------------------------------------------------------------------------------------------- !> @brief return Poisson ratio !-------------------------------------------------------------------------------------------------- -module function elastic_nu(ph,en) result(nu) +pure module function elastic_nu(ph,en) result(nu) integer, intent(in) :: & ph, & From 4a0a1f7ac997fedc430e61abf317879b980a1bd7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 2 Jan 2022 22:29:48 +0100 Subject: [PATCH 66/74] paper is online --- .../config/phase/mechanical/plastic/dislotwin_IF-steel.yaml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml b/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml index 362797d0b..4e4ff1a9f 100644 --- a/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml +++ b/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml @@ -4,7 +4,8 @@ references: International Journal of Plasticity 134:102779, 2020, https://doi.org/10.1016/j.ijplas.2020.102779 - K. Sedighiani et al., - Mechanics of Materials, submitted + Mechanics of Materials, 164:104117, 2022, + https://doi.org/10.1016/j.mechmat.2021.104117 output: [rho_dip, rho_mob] N_sl: [12, 12] b_sl: [2.49e-10, 2.49e-10] From 3dd37cdf229cec09b73efe182fc9b74d166a9987 Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 3 Jan 2022 02:19:10 +0100 Subject: [PATCH 67/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-338-g4a0a1f7ac --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 7a82d0bd4..0a2f00913 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-336-g6871eb302 +v3.0.0-alpha5-338-g4a0a1f7ac From 510a26ded90dfebee28c9b5a05fb362e13471013 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 3 Jan 2022 10:36:23 +0100 Subject: [PATCH 68/74] is used in dislotwin according to original paper --- src/constants.f90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/constants.f90 b/src/constants.f90 index dd26ce78c..43bb60953 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -9,7 +9,8 @@ module constants public real(pReal), parameter :: & - T_ROOM = 300.0_pReal, & !< Room temperature in K - K_B = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin + T_ROOM = 300.0_pReal, & !< Room temperature in K. ToDo: IUPAC: 298.15 + K_B = 1.38e-23_pReal, & !< Boltzmann constant in J/Kelvin + N_A = 6.02214076e-23_pReal !< Avogadro constant in 1/mol end module constants From c51976e2a9b44598149b200196f22244562080eb Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 3 Jan 2022 15:47:34 +0100 Subject: [PATCH 69/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-340-g510a26ded --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 0a2f00913..b273e2008 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-338-g4a0a1f7ac +v3.0.0-alpha5-340-g510a26ded From 09d68b1fff9fcf20e2255e786d9d3e15de297abd Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 3 Jan 2022 21:46:38 +0100 Subject: [PATCH 70/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-342-ga177e32ff --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 0a2f00913..701bcd98b 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-338-g4a0a1f7ac +v3.0.0-alpha5-342-ga177e32ff From b12a18097488ba4e8a7794099791eb86ac8d4845 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Thu, 6 Jan 2022 09:54:33 -0500 Subject: [PATCH 71/74] added literature reference for constitutive law --- src/phase_mechanical_plastic_kinehardening.f90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index 6c0e1e0cc..03eb27f31 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -86,6 +86,8 @@ module function plastic_kinehardening_init() result(myPlasticity) print'(/,1x,a)', '<<<+- phase:mechanical:plastic:kinehardening init -+>>>' print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) + print'(/,1x,a)', 'J.A. Wollmershauser et al., International Journal of Fatigue 36:181–193, 2012' + print'( 1x,a)', 'https://doi.org/10.1016/j.ijfatigue.2011.07.008' phases => config_material%get('phase') allocate(param(phases%length)) From a2f6475baa6e06d6550aa183aeaa013ff98eeeaf Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 7 Jan 2022 07:23:43 +0100 Subject: [PATCH 72/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-355-gc29428a60 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index b273e2008..a467f728c 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-340-g510a26ded +v3.0.0-alpha5-355-gc29428a60 From e3a233a681ea7c827063806777ff7ce90573e867 Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 7 Jan 2022 21:38:12 +0100 Subject: [PATCH 73/74] [skip ci] updated version information after successful test of v3.0.0-alpha5-358-g81a7c32a5 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index a467f728c..09d546001 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-355-gc29428a60 +v3.0.0-alpha5-358-g81a7c32a5 From dd6e9a016efe823fa06c250d6f6f2d59e81226c6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 8 Jan 2022 12:07:29 +0000 Subject: [PATCH 74/74] just off by 46 orders of magnitude ;) --- src/constants.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/constants.f90 b/src/constants.f90 index 43bb60953..3b99763de 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -11,6 +11,6 @@ module constants real(pReal), parameter :: & T_ROOM = 300.0_pReal, & !< Room temperature in K. ToDo: IUPAC: 298.15 K_B = 1.38e-23_pReal, & !< Boltzmann constant in J/Kelvin - N_A = 6.02214076e-23_pReal !< Avogadro constant in 1/mol + N_A = 6.02214076e23_pReal !< Avogadro constant in 1/mol end module constants