From bd0eb0f154327b9773a38f539313848fe663cc55 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Mon, 18 Dec 2023 17:21:21 +0100 Subject: [PATCH 01/33] load -> yaml --- examples/mesh/tensionY_mono.load | 14 -------------- examples/mesh/tensionY_mono.yaml | 22 ++++++++++++++++++++++ examples/mesh/tensionZ_3g.load | 14 -------------- examples/mesh/tensionZ_3g.yaml | 22 ++++++++++++++++++++++ examples/mesh/tensionZ_5g.load | 18 ------------------ examples/mesh/tensionZ_5g.yaml | 22 ++++++++++++++++++++++ 6 files changed, 66 insertions(+), 46 deletions(-) delete mode 100644 examples/mesh/tensionY_mono.load create mode 100644 examples/mesh/tensionY_mono.yaml delete mode 100644 examples/mesh/tensionZ_3g.load create mode 100644 examples/mesh/tensionZ_3g.yaml delete mode 100644 examples/mesh/tensionZ_5g.load create mode 100644 examples/mesh/tensionZ_5g.yaml diff --git a/examples/mesh/tensionY_mono.load b/examples/mesh/tensionY_mono.load deleted file mode 100644 index e7ab69fc4..000000000 --- a/examples/mesh/tensionY_mono.load +++ /dev/null @@ -1,14 +0,0 @@ -# initial elastic step -$Loadcase 1 t 0.0005 N 1 f_out 1 - Face 3 Y -0.025 - Face 4 X 0.0 - Face 4 Y 0.0 - Face 4 Z 0.0 -$EndLoadcase -# plastic step -$Loadcase 2 t 1.0 N 10 f_out 2 - Face 3 Y -0.025 - Face 4 X 0.0 - Face 4 Y 0.0 - Face 4 Z 0.0 -$EndLoadcase diff --git a/examples/mesh/tensionY_mono.yaml b/examples/mesh/tensionY_mono.yaml new file mode 100644 index 000000000..fb7da1b86 --- /dev/null +++ b/examples/mesh/tensionY_mono.yaml @@ -0,0 +1,22 @@ +--- +loadstep: + - boundary_conditions: + mechanical: + - dot_u: ['x', -0.025, 'x'] + tag: 3 + - dot_u: [0.0, 0.0, 0.0] + tag: 4 + discretization: + t: 0.0005 + N: 1 + f_out: 1 + - boundary_conditions: + mechanical: + - dot_u: ['x', -0.025, 'x'] + tag: 3 + - dot_u: [0.0, 0.0, 0.0] + tag: 4 + discretization: + t: 1.0 + N: 10 + f_out: 2 diff --git a/examples/mesh/tensionZ_3g.load b/examples/mesh/tensionZ_3g.load deleted file mode 100644 index b888873ea..000000000 --- a/examples/mesh/tensionZ_3g.load +++ /dev/null @@ -1,14 +0,0 @@ -# initial elastic step -$Loadcase 1 t 0.0005 N 1 f_out 1 - Face 1 Z 0.01 - Face 2 X 0.0 - Face 2 Y 0.0 - Face 2 Z 0.0 -$EndLoadcase -# plastic step -$Loadcase 2 t 1.0 N 10 f_out 2 - Face 1 Z 0.01 - Face 2 X 0.0 - Face 2 Y 0.0 - Face 2 Z 0.0 -$EndLoadcase diff --git a/examples/mesh/tensionZ_3g.yaml b/examples/mesh/tensionZ_3g.yaml new file mode 100644 index 000000000..e10f80229 --- /dev/null +++ b/examples/mesh/tensionZ_3g.yaml @@ -0,0 +1,22 @@ +--- +loadstep: + - boundary_conditions: + mechanical: + - dot_u: ['x', 'x', 0.01] + tag: 1 + - dot_u: [0.0, 0.0, 0.0] + tag: 2 + discretization: + t: 0.0005 + N: 1 + f_out: 1 + - boundary_conditions: + mechanical: + - dot_u: ['x', 'x', 0.01] + tag: 1 + - dot_u: [0.0, 0.0, 0.0] + tag: 2 + discretization: + t: 1.0 + N: 10 + f_out: 2 diff --git a/examples/mesh/tensionZ_5g.load b/examples/mesh/tensionZ_5g.load deleted file mode 100644 index 475931523..000000000 --- a/examples/mesh/tensionZ_5g.load +++ /dev/null @@ -1,18 +0,0 @@ -# initial elastic step -$Loadcase 1 t 0.0005 N 1 f_out 1 - Face 1 X 0.0 - Face 1 Y 0.0 - Face 1 Z 0.0 - Face 2 X 0.0 - Face 2 Y 0.0 - Face 2 Z 0.0025 -$EndLoadcase -# plastic step -$Loadcase 2 t 1.0 N 10 f_out 2 - Face 1 X 0.0 - Face 1 Y 0.0 - Face 1 Z 0.0 - Face 2 X 0.0 - Face 2 Y 0.0 - Face 2 Z 0.0025 -$EndLoadcase diff --git a/examples/mesh/tensionZ_5g.yaml b/examples/mesh/tensionZ_5g.yaml new file mode 100644 index 000000000..f43e1e326 --- /dev/null +++ b/examples/mesh/tensionZ_5g.yaml @@ -0,0 +1,22 @@ +--- +loadstep: + - boundary_conditions: + mechanical: + - dot_u: [0.0, 0.0, 0.0] + tag: 1 + - dot_u: [0.0, 0.0, 0.0025] + tag: 2 + discretization: + t: 0.0005 + N: 1 + f_out: 1 + - boundary_conditions: + mechanical: + - dot_u: [0.0, 0.0, 0.0] + tag: 1 + - dot_u: [0.0, 0.0, 0.0025] + tag: 2 + discretization: + t: 1.0 + N: 10 + f_out: 2 From 7245aa24d4ed99a8f6b3625893ea2f0575abc1b6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 18 Dec 2023 23:55:38 +0100 Subject: [PATCH 02/33] use NotImplementedError to indicate non-implemented combinations --- python/damask/_result.py | 35 +++++++++++++++++++++-------------- python/damask/util.py | 2 +- python/tests/test_Result.py | 6 +++++- 3 files changed, 27 insertions(+), 16 deletions(-) diff --git a/python/damask/_result.py b/python/damask/_result.py index 6ffbd0352..5b6ac71d5 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -115,8 +115,6 @@ class Result: self.cells = f['geometry'].attrs['cells'] self.size = f['geometry'].attrs['size'] self.origin = f['geometry'].attrs['origin'] - else: - self.add_curl = self.add_divergence = self.add_gradient = None # type: ignore r = re.compile(rf'{prefix_inc}([0-9]+)') self.increments = sorted([i for i in f.keys() if r.match(i)],key=util.natural_sort) @@ -1313,8 +1311,8 @@ class Result: Notes ----- - This function is only available for structured grids, - i.e. fields resulting from the grid solver. + This function is implemented only for structured grids + with one constituent and a single phase. """ def curl(f: DADF5Dataset, size: np.ndarray) -> DADF5Dataset: @@ -1342,8 +1340,8 @@ class Result: Notes ----- - This function is only available for structured grids, - i.e. fields resulting from the grid solver. + This function is implemented only for structured grids + with one constituent and a single phase. """ def divergence(f: DADF5Dataset, size: np.ndarray) -> DADF5Dataset: @@ -1371,8 +1369,8 @@ class Result: Notes ----- - This function is only available for structured grids, - i.e. fields resulting from the grid solver. + This function is implemented only for structured grids + with one constituent and a single phase. """ def gradient(f: DADF5Dataset, size: np.ndarray) -> DADF5Dataset: @@ -1410,13 +1408,13 @@ class Result: Arguments parsed to func. """ - if len(datasets) != 1 or self.N_constituents != 1: - raise NotImplementedError + if self.N_constituents != 1 or len(datasets) != 1 or not self.structured: + raise NotImplementedError('not a structured grid with one constituent and a single phase') at_cell_ph,in_data_ph,at_cell_ho,in_data_ho = self._mappings() increments = self.place(list(datasets.values()),False) - if not increments: raise RuntimeError("received invalid dataset") + if not increments: raise RuntimeError('received invalid dataset') with h5py.File(self.fname, 'a') as f: for increment in increments.items(): for ty in increment[1].items(): @@ -1722,9 +1720,14 @@ class Result: Defaults to False, i.e. the XDMF file expects the DADF5 file at a stable relative path. + Notes + ----- + This function is implemented only for structured grids with + one constituent and a single phase. + """ if self.N_constituents != 1 or len(self.phases) != 1 or not self.structured: - raise TypeError('XDMF output requires structured grid with single phase and single constituent.') + raise NotImplementedError('not a structured grid with one constituent and a single phase') attribute_type_map = defaultdict(lambda:'Matrix', ( ((),'Scalar'), ((3,),'Vector'), ((3,3),'Tensor')) ) @@ -1949,6 +1952,11 @@ class Result: target_dir : str or pathlib.Path, optional Directory to save DREAM3D files. Will be created if non-existent. + Notes + ----- + This function is implemented only for structured grids with + one constituent. + """ def add_attribute(obj,name,data): """DREAM.3D requires fixed length string.""" @@ -1964,11 +1972,10 @@ class Result: return obj[name] if self.N_constituents != 1 or not self.structured: - raise TypeError('DREAM3D output requires structured grid with single constituent.') + raise NotImplementedError('not a structured grid with one constituent') N_digits = int(np.floor(np.log10(max(1,self.incs[-1]))))+1 - at_cell_ph,in_data_ph,_,_ = self._mappings() out_dir = Path.cwd() if target_dir is None else Path(target_dir) diff --git a/python/damask/util.py b/python/damask/util.py index eda668564..0688fca80 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -618,7 +618,7 @@ def _docstringer(docstring: _Union[str, _Callable], adopted_, flags=_re.MULTILINE|_re.DOTALL).group('content') except AttributeError: - raise RuntimeError(f"Function docstring passed for docstring section '{key}' is invalid:\n{docstring}") + raise RuntimeError(f"function docstring passed for docstring section '{key}' is invalid:\n{docstring}") docstring_indent, adopted_indent = (min([len(line)-len(line.lstrip()) for line in section.split('\n') if line.strip()]) for section in [docstring_, adopted_]) diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index bad55c98d..d6bbd63a8 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -478,6 +478,10 @@ class TestResult: for attr in dset.attrs: assert np.array_equal(dset.attrs[attr],cur[path].attrs[attr]) + def test_export_DREAM3D_invalid(self,res_path): + with pytest.raises(NotImplementedError): + Result(res_path/'4grains2x4x3_compressionY.hdf5').export_DREAM3D() + def test_XDMF_datatypes(self,tmp_path,single_phase,update,res_path): for what,shape in {'scalar':(),'vector':(3,),'tensor':(3,3),'matrix':(12,)}.items(): @@ -509,7 +513,7 @@ class TestResult: assert dim_vti == dim_xdmf and bounds_vti == bounds_xdmf def test_XDMF_invalid(self,default): - with pytest.raises(TypeError): + with pytest.raises(NotImplementedError): default.export_XDMF() def test_XDMF_custom_path(self,single_phase,tmp_path): From c4d061ba0a645b9080d03888925bbbddc49f6e0c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 19 Dec 2023 20:05:52 +0100 Subject: [PATCH 03/33] F2023 tokenize (second form) --- src/IO.f90 | 78 +++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 77 insertions(+), 1 deletion(-) diff --git a/src/IO.f90 b/src/IO.f90 index a08c624c3..74f23a678 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -48,7 +48,8 @@ implicit none(type,external) IO_color, & IO_error, & IO_warning, & - IO_STDOUT + IO_STDOUT, & + tokenize contains @@ -742,6 +743,33 @@ pure function CRLF2LF(str) end function CRLF2LF +!-------------------------------------------------------------------------------------------------- +!> @brief Fortran 2023 tokenize (first form). +!-------------------------------------------------------------------------------------------------- +pure subroutine tokenize(string,set,tokens) + + character(len=*), intent(in) :: string, set + character(len=:), dimension(:), allocatable, intent(out) :: tokens + + integer, allocatable, dimension(:,:) :: pos + integer :: i, s, e + + + allocate(pos(2,0)) + e = 0 + do while (e < verify(string,set,back=.true.)) + s = e + merge(verify(string(e+1:),set),1,scan(string(e+1:),set)/=0) + e = s + merge(scan(string(s:),set)-2,len(string(s:))-1,scan(string(s:),set)/=0) + pos = reshape([pos,[s,e]],[2,size(pos)/2+1]) + end do + allocate(character(len=merge(maxval(pos(2,:)-pos(1,:))+1,0,size(pos)>0))::tokens(size(pos,2))) + do i = 1, size(pos,2) + tokens(i) = string(pos(1,i):pos(2,i)) + end do + +end subroutine tokenize + + !-------------------------------------------------------------------------------------------------- !> @brief Write statements to standard error. !-------------------------------------------------------------------------------------------------- @@ -808,6 +836,7 @@ subroutine IO_selfTest() integer, dimension(:), allocatable :: chunkPos character(len=:), allocatable :: str,out + character(len=:), dimension(:), allocatable :: tokens if (dNeq(1.0_pREAL, IO_strAsReal('1.0'))) error stop 'IO_strAsReal' @@ -887,6 +916,53 @@ subroutine IO_selfTest() if ('abc,'//IO_EOL//'xxdefg,'//IO_EOL//'xxhij' /= IO_wrapLines('abc,defg, hij',filler='xx',length=4)) & error stop 'IO_wrapLines/7' + call tokenize('','$',tokens) + if (size(tokens) /= 0 .or. len(tokens) /=0) error stop 'tokenize empty' + call tokenize('abcd','dcba',tokens) + if (size(tokens) /= 0 .or. len(tokens) /=0) error stop 'tokenize empty' + + tokens=['a'] + call test_tokenize('a','#',tokens) + call test_tokenize('#a','#',tokens) + call test_tokenize('a#','#',tokens) + + tokens=['aa'] + call test_tokenize('aa','#',tokens) + call test_tokenize('$aa','$',tokens) + call test_tokenize('aa$','$',tokens) + + tokens=['a','b'] + call test_tokenize('a$b','$',tokens) + call test_tokenize('@a@$b@','$@',tokens) + + tokens=['aa','bb'] + call test_tokenize('aa$bb','$',tokens) + call test_tokenize('aa$$bb','$',tokens) + call test_tokenize('aa$bb$','$',tokens) + + tokens=['aa ','bbb ','cccc'] + call test_tokenize('aa$bbb$cccc','$',tokens) + call test_tokenize('$aa$bbb$cccc$','$',tokens) + call tokenize('#aa@@bbb!!!cccc#','#@!',tokens) + + + contains + subroutine test_tokenize(input,delimiter,solution) + + character(len=*), intent(in) :: input, delimiter + character(len=*), dimension(:), intent(in) :: solution + + character(len=:), dimension(:), allocatable :: tok + integer :: i + + + call tokenize(input,delimiter,tok) + do i = 1,size(tok) + if (solution(i) /= tok(i)) error stop 'tokenize "'//solution(i)//'" vs. "'//tok(i)//'"' + end do + + end subroutine test_tokenize + end subroutine IO_selfTest end module IO From ace7d8f00340a98a09bc9126e24868c93f80129e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 19 Dec 2023 21:19:57 +0100 Subject: [PATCH 04/33] use tokenize (F2023) --- src/grid/VTI.f90 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/grid/VTI.f90 b/src/grid/VTI.f90 index 2749c1bb6..976d11067 100644 --- a/src/grid/VTI.f90 +++ b/src/grid/VTI.f90 @@ -201,26 +201,26 @@ subroutine cellsSizeOrigin(c,s,o,header) real(pREAL), dimension(3), intent(out) :: s,o character(len=*), intent(in) :: header - character(len=:), allocatable :: temp + character(len=:), allocatable, dimension(:) :: temp real(pREAL), dimension(3) :: delta integer :: i - temp = getXMLValue(header,'Direction') - if (temp /= '1 0 0 0 1 0 0 0 1' .and. temp /= '') & ! https://discourse.vtk.org/t/vti-specification/6526 + temp = [getXMLValue(header,'Direction')] + if (temp(1) /= '1 0 0 0 1 0 0 0 1' .and. temp(1) /= '') & ! https://discourse.vtk.org/t/vti-specification/6526 call IO_error(error_ID = 844, ext_msg = 'coordinate order') - temp = getXMLValue(header,'WholeExtent') - if (any([(IO_intValue(temp,IO_strPos(temp),i),i=1,5,2)] /= 0)) & + call tokenize(getXMLValue(header,'WholeExtent'),' ',temp) + if (any([(IO_strAsInt(temp(i)),i=1,5,2)] /= 0)) & call IO_error(error_ID = 844, ext_msg = 'coordinate start') - c = [(IO_intValue(temp,IO_strPos(temp),i),i=2,6,2)] + c = [(IO_strAsInt(temp(i)),i=2,6,2)] - temp = getXMLValue(header,'Spacing') - delta = [(IO_realValue(temp,IO_strPos(temp),i),i=1,3)] + call tokenize(getXMLValue(header,'Spacing'),' ',temp) + delta = [(IO_strAsReal(temp(i)),i=1,3)] s = delta * real(c,pREAL) - temp = getXMLValue(header,'Origin') - o = [(IO_realValue(temp,IO_strPos(temp),i),i=1,3)] + call tokenize(getXMLValue(header,'Origin'),' ',temp) + o = [(IO_strAsReal(temp(i)),i=1,3)] end subroutine cellsSizeOrigin From 51bccd98d535b4d706c72f6b59e75ca6ea584780 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 20 Dec 2023 11:27:57 +0100 Subject: [PATCH 05/33] check is already done centrally --- cmake/Compiler-GNU.cmake | 4 ---- cmake/Compiler-Intel.cmake | 4 ---- cmake/Compiler-IntelLLVM.cmake | 6 +----- 3 files changed, 1 insertion(+), 13 deletions(-) diff --git a/cmake/Compiler-GNU.cmake b/cmake/Compiler-GNU.cmake index a62875d05..ee58f8cc7 100644 --- a/cmake/Compiler-GNU.cmake +++ b/cmake/Compiler-GNU.cmake @@ -1,10 +1,6 @@ ################################################################################################### # GNU Compiler ################################################################################################### -if (CMAKE_Fortran_COMPILER_VERSION VERSION_LESS 9.0) - message (FATAL_ERROR "GCC Compiler version: ${CMAKE_Fortran_COMPILER_VERSION} not supported") -endif () - if (OPENMP) set (OPENMP_FLAGS "-fopenmp") endif () diff --git a/cmake/Compiler-Intel.cmake b/cmake/Compiler-Intel.cmake index 19e75a9fa..2182459ed 100644 --- a/cmake/Compiler-Intel.cmake +++ b/cmake/Compiler-Intel.cmake @@ -1,10 +1,6 @@ ################################################################################################### # 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 () diff --git a/cmake/Compiler-IntelLLVM.cmake b/cmake/Compiler-IntelLLVM.cmake index 76fc0386a..4cdd2588d 100644 --- a/cmake/Compiler-IntelLLVM.cmake +++ b/cmake/Compiler-IntelLLVM.cmake @@ -1,10 +1,6 @@ ################################################################################################### -# Intel Compiler +# IntelLLVM 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 "-fiopenmp") endif () From 4e3b9e658621f6aa693db4ee7263e649aa35ce51 Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 20 Dec 2023 16:58:04 +0100 Subject: [PATCH 06/33] [skip ci] updated version information after successful test of v3.0.0-alpha8-172-gec624a86a --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 4f15e26cc..a45b5eb13 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -3.0.0-alpha8-169-g85191cd02 +3.0.0-alpha8-172-gec624a86a From c0b928021875b73b9c79d35494b1766074d8620d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 19 Dec 2023 21:26:31 +0100 Subject: [PATCH 07/33] variables in stop code not supported by older compilers --- src/IO.f90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index 74f23a678..4c0676ec8 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -919,7 +919,7 @@ subroutine IO_selfTest() call tokenize('','$',tokens) if (size(tokens) /= 0 .or. len(tokens) /=0) error stop 'tokenize empty' call tokenize('abcd','dcba',tokens) - if (size(tokens) /= 0 .or. len(tokens) /=0) error stop 'tokenize empty' + if (size(tokens) /= 0 .or. len(tokens) /=0) error stop 'tokenize only separators' tokens=['a'] call test_tokenize('a','#',tokens) @@ -958,7 +958,8 @@ subroutine IO_selfTest() call tokenize(input,delimiter,tok) do i = 1,size(tok) - if (solution(i) /= tok(i)) error stop 'tokenize "'//solution(i)//'" vs. "'//tok(i)//'"' + !if (solution(i) /= tok(i)) error stop 'tokenize "'//solution(i)//'" vs. "'//tok(i)//'"' ! requires 2018 standard + if (solution(i) /= tok(i)) error stop 'tokenize' end do end subroutine test_tokenize From 8458ea5ecf85266586a9b636791314e02337a344 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 20 Dec 2023 06:30:58 +0100 Subject: [PATCH 08/33] simplified string conversion and related error handling checking for valid characters is not sufficient because signs and exponents are allowed to appear only at specific locations. Since Fortrans internal file read will anyways complain about invalid characters, the check doesn't even have a value. reduced number of error codes by grouping all conversion related errors --- src/IO.f90 | 37 +++++++++++-------------------------- 1 file changed, 11 insertions(+), 26 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index a08c624c3..b8e4d0022 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -380,17 +380,11 @@ integer function IO_strAsInt(str) character(len=*), intent(in) :: str !< string for conversion to int value - integer :: readStatus - character(len=*), parameter :: VALIDCHARS = '0123456789+- ' + integer :: readStatus - valid: if (verify(str,VALIDCHARS) == 0) then - read(str,*,iostat=readStatus) IO_strAsInt - if (readStatus /= 0) call IO_error(111,str) - else valid - IO_strAsInt = 0 - call IO_error(111,str) - end if valid + read(str,*,iostat=readStatus) IO_strAsInt + if (readStatus /= 0) call IO_error(111,'cannot represent "'//str//'" as integer') end function IO_strAsInt @@ -402,27 +396,23 @@ real(pREAL) function IO_strAsReal(str) character(len=*), intent(in) :: str !< string for conversion to real value - integer :: readStatus - character(len=*), parameter :: VALIDCHARS = '0123456789eE.+- ' + integer :: readStatus - valid: if (verify(str,VALIDCHARS) == 0) then - read(str,*,iostat=readStatus) IO_strAsReal - if (readStatus /= 0) call IO_error(112,str) - else valid - IO_strAsReal = 0.0_pREAL - call IO_error(112,str) - end if valid + read(str,*,iostat=readStatus) IO_strAsReal + if (readStatus /= 0) call IO_error(111,'cannot represent "'//str//'" as real') end function IO_strAsReal !-------------------------------------------------------------------------------------------------- !> @brief Return logical value from given string. +!> @details: 'True' and 'true' are converted to .true. +!> @details: 'False' and 'false' are converted to .false. !-------------------------------------------------------------------------------------------------- logical function IO_strAsBool(str) - character(len=*), intent(in) :: str !< string for conversion to int value + character(len=*), intent(in) :: str !< string for conversion to boolean if (trim(adjustl(str)) == 'True' .or. trim(adjustl(str)) == 'true') then @@ -430,8 +420,7 @@ logical function IO_strAsBool(str) elseif (trim(adjustl(str)) == 'False' .or. trim(adjustl(str)) == 'false') then IO_strAsBool = .false. else - IO_strAsBool = .false. - call IO_error(113,str) + call IO_error(111,'cannot represent "'//str//'" as boolean') end if end function IO_strAsBool @@ -498,11 +487,7 @@ subroutine IO_error(error_ID,ext_msg,label1,ID1,label2,ID2) case (110) msg = 'invalid chunk selected' case (111) - msg = 'invalid character for int:' - case (112) - msg = 'invalid character for real:' - case (113) - msg = 'invalid character for logical:' + msg = 'invalid string for conversion' case (114) msg = 'cannot decode base64 string:' From 808ab090530740b9016d7089819e382d33965862 Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 21 Dec 2023 00:13:10 +0100 Subject: [PATCH 09/33] [skip ci] updated version information after successful test of v3.0.0-alpha8-175-g5c71238bc --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index a45b5eb13..a9ff51081 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -3.0.0-alpha8-172-gec624a86a +3.0.0-alpha8-175-g5c71238bc From bdb30118fb2aa6757a42dcd4648cd56b013dc1e9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 20 Dec 2023 17:39:28 +0100 Subject: [PATCH 10/33] starting to support new LLVM compiler does not work yet, but good to have for testing once in a while. Current version is 17, require minimum 19 to make clear that it does not work --- CMakeLists.txt | 3 +++ cmake/Compiler-LLVMFlang.cmake | 12 ++++++++++++ 2 files changed, 15 insertions(+) create mode 100644 cmake/Compiler-LLVMFlang.cmake diff --git a/CMakeLists.txt b/CMakeLists.txt index a501717ed..b630f7e46 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -101,6 +101,9 @@ elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "Intel") elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "IntelLLVM") include(Compiler-IntelLLVM) set(Fortran_COMPILER_VERSION_MIN 19) +elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "LLVMFlang") + include(Compiler-LLVMFlang) + set(Fortran_COMPILER_VERSION_MIN 19) else() message(FATAL_ERROR "Compiler '${CMAKE_Fortran_COMPILER_ID}' not supported") endif() diff --git a/cmake/Compiler-LLVMFlang.cmake b/cmake/Compiler-LLVMFlang.cmake new file mode 100644 index 000000000..c6b0405bc --- /dev/null +++ b/cmake/Compiler-LLVMFlang.cmake @@ -0,0 +1,12 @@ +################################################################################################### +# LLVM Compiler +################################################################################################### +if (OPENMP) + set (OPENMP_FLAGS "-fopenmp") +endif () + +#------------------------------------------------------------------------------------------------ +# Fine tuning compilation options +set (COMPILE_FLAGS "${COMPILE_FLAGS} -cpp") +# preprocessor + From 0518cc67308296b14050497921aea8747b340d41 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 21 Dec 2023 07:54:20 +0100 Subject: [PATCH 11/33] slighly polished yaml files --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 117b65d85..62df7f24f 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 117b65d852659158c0f4ca3bf8ce8db51a7a1961 +Subproject commit 62df7f24f2a95fda255f7d20b130afcfeecb1b4a From 1933c896b28b86c99966de35c0ce26cc36218b51 Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 21 Dec 2023 12:41:31 +0100 Subject: [PATCH 12/33] [skip ci] updated version information after successful test of v3.0.0-alpha8-179-gb9d166652 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index a9ff51081..fa7f3d18a 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -3.0.0-alpha8-175-g5c71238bc +3.0.0-alpha8-179-gb9d166652 From 0e4c20a43042a20b02c247f7581cd3b8d3026c6b Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 21 Dec 2023 22:10:08 +0100 Subject: [PATCH 13/33] [skip ci] updated version information after successful test of v3.0.0-alpha8-186-gaf5bbed00 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index fa7f3d18a..007956b90 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -3.0.0-alpha8-179-gb9d166652 +3.0.0-alpha8-186-gaf5bbed00 From 0c56711032b893e36f0df6dbac067222b48aa9cd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 22 Dec 2023 19:06:02 +0100 Subject: [PATCH 14/33] use OS-independent separators fix tests on Windows --- python/tests/test_Result.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index d6bbd63a8..c7e4187ab 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -63,7 +63,7 @@ def h5py_dataset_iterator(): """Iterate over all datasets in an HDF5 file.""" def _h5py_dataset_iterator(g, prefix=''): for key,item in g.items(): - path = os.path.join(prefix, key) + path = '/'.join([prefix, key]) if isinstance(item, h5py.Dataset): # test for dataset yield (path, item) elif isinstance(item, h5py.Group): # test for group (go down) @@ -472,7 +472,7 @@ class TestResult: c = [_.decode() for _ in cur[path]] r = ['Unknown Phase Type'] + result.phases assert c == r - grp = os.path.split(path)[0] + grp = str(path).rpartition('/')[0] for attr in ref[grp].attrs: assert np.array_equal(ref[grp].attrs[attr],cur[grp].attrs[attr]) for attr in dset.attrs: From 696cccaa6620743bb502e2cae29864f3738213cc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 22 Dec 2023 22:37:47 +0100 Subject: [PATCH 15/33] rho_forest is a dependent state --- src/phase_mechanical_plastic_nonlocal.f90 | 37 +++++++++++------------ 1 file changed, 18 insertions(+), 19 deletions(-) diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 545dec4e6..361ecc98d 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -124,7 +124,8 @@ submodule(phase:plastic) nonlocal type :: tNonlocalDependentState real(pREAL), allocatable, dimension(:,:) :: & tau_pass, & - tau_back + tau_back, & + rho_forest real(pREAL), allocatable, dimension(:,:,:,:,:) :: & compatibility end type tNonlocalDependentState @@ -146,7 +147,6 @@ submodule(phase:plastic) nonlocal rhoDip, & rho_dip_edg, & rho_dip_scr, & - rho_forest, & gamma, & v, & v_edg_pos, & @@ -177,7 +177,7 @@ module function plastic_nonlocal_init() result(myPlasticity) integer :: & ph, & Nmembers, & - sizeState, sizeDotState, sizeDependentState, sizeDeltaState, & + sizeState, sizeDotState, sizeDeltaState, & s1, s2, & s, t, l real(pREAL), dimension(:,:), allocatable :: & @@ -389,8 +389,7 @@ module function plastic_nonlocal_init() result(myPlasticity) 'rhoSglScrewPosImmobile','rhoSglScrewNegImmobile', & 'rhoDipEdge ','rhoDipScrew ', & 'gamma ' ]) * prm%sum_N_sl !< "basic" microstructural state variables that are independent from other state variables - sizeDependentState = size([ 'rhoForest ']) * prm%sum_N_sl !< microstructural state variables that depend on other state variables - sizeState = sizeDotState + sizeDependentState & + sizeState = sizeDotState & + size([ 'velocityEdgePos ','velocityEdgeNeg ', & 'velocityScrewPos ','velocityScrewNeg ', & 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure @@ -477,15 +476,15 @@ module function plastic_nonlocal_init() result(myPlasticity) if (any(plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_pREAL)) & extmsg = trim(extmsg)//' atol_gamma' - stt%rho_forest => plasticState(ph)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nmembers) - stt%v => plasticState(ph)%state (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nmembers) - stt%v_edg_pos => plasticState(ph)%state (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nmembers) - stt%v_edg_neg => plasticState(ph)%state (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nmembers) - stt%v_scr_pos => plasticState(ph)%state (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nmembers) - stt%v_scr_neg => plasticState(ph)%state (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nmembers) + stt%v => plasticState(ph)%state (11*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nmembers) + stt%v_edg_pos => plasticState(ph)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nmembers) + stt%v_edg_neg => plasticState(ph)%state (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nmembers) + stt%v_scr_pos => plasticState(ph)%state (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nmembers) + stt%v_scr_neg => plasticState(ph)%state (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nmembers) allocate(dst%tau_pass(prm%sum_N_sl,Nmembers),source=0.0_pREAL) allocate(dst%tau_back(prm%sum_N_sl,Nmembers),source=0.0_pREAL) + allocate(dst%rho_forest(prm%sum_N_sl,Nmembers),source=0.0_pREAL) allocate(dst%compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,Nmembers),source=0.0_pREAL) end associate @@ -518,7 +517,7 @@ module function plastic_nonlocal_init() result(myPlasticity) iRhoU(s,t,ph) = l end do end do - l = l + (4+2+1+1)*param(ph)%sum_N_sl ! immobile(4), dipole(2), shear, forest + l = l + (4+2+1)*param(ph)%sum_N_sl ! immobile(4), dipole(2), shear do t = 1,4 do s = 1,param(ph)%sum_N_sl l = l + 1 @@ -602,7 +601,7 @@ module subroutine nonlocal_dependentState(ph, en) nu = elastic_nu(ph,en,prm%isotropic_bound) rho = getRho(ph,en) - stt%rho_forest(:,en) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) & + dst%rho_forest(:,en) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) & + matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2)) @@ -612,7 +611,7 @@ module subroutine nonlocal_dependentState(ph, en) myInteractionMatrix = prm%h_sl_sl & * spread(( 1.0_pREAL - prm%f_F & + prm%f_F & - * log(0.35_pREAL * prm%b_sl * sqrt(max(stt%rho_forest(:,en),prm%rho_significant))) & + * log(0.35_pREAL * prm%b_sl * sqrt(max(dst%rho_forest(:,en),prm%rho_significant))) & / log(0.35_pREAL * prm%b_sl * 1e6_pREAL))**2,2,prm%sum_N_sl) else myInteractionMatrix = prm%h_sl_sl @@ -1018,17 +1017,17 @@ module subroutine nonlocal_dotState(Mp,timestep, & isBCC: if (phase_lattice(ph) == 'cI') then forall (s = 1:prm%sum_N_sl, sum(abs(v(s,1:4))) > 0.0_pREAL) rhoDotMultiplication(s,1:2) = sum(abs(dot_gamma(s,3:4))) / prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication - * sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! & ! mean free path + * sqrt(dst%rho_forest(s,en)) / prm%i_sl(s) ! & ! mean free path ! * 2.0_pREAL * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4))) ! ratio of screw to overall velocity determines edge generation rhoDotMultiplication(s,3:4) = sum(abs(dot_gamma(s,3:4))) /prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication - * sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! & ! mean free path + * sqrt(dst%rho_forest(s,en)) / prm%i_sl(s) ! & ! mean free path ! * 2.0_pREAL * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation endforall else isBCC rhoDotMultiplication(:,1:4) = spread( & (sum(abs(dot_gamma(:,1:2)),2) * prm%f_ed_mult + sum(abs(dot_gamma(:,3:4)),2)) & - * sqrt(stt%rho_forest(:,en)) / prm%i_sl / prm%b_sl, 2, 4) ! eq. 3.26 + * sqrt(dst%rho_forest(:,en)) / prm%i_sl / prm%b_sl, 2, 4) ! eq. 3.26 end if isBCC forall (s = 1:prm%sum_N_sl, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en) @@ -1074,7 +1073,7 @@ module subroutine nonlocal_dotState(Mp,timestep, & if (phase_lattice(ph) == 'cF') & forall (s = 1:prm%sum_N_sl, prm%colinearSystem(s) > 0) & rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) & - * 0.25_pREAL * sqrt(stt%rho_forest(s,en)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed + * 0.25_pREAL * sqrt(dst%rho_forest(s,en)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed ! thermally activated annihilation of edge dipoles by climb @@ -1486,7 +1485,7 @@ module subroutine plastic_nonlocal_result(ph,group) call result_writeDataset(stt%rho_dip_scr,group,trim(prm%output(ou)), & 'screw dipole density','1/m²', prm%systems_sl) case('rho_f') - call result_writeDataset(stt%rho_forest,group,trim(prm%output(ou)), & + call result_writeDataset(dst%rho_forest,group,trim(prm%output(ou)), & 'forest density','1/m²', prm%systems_sl) case('v_ed_pos') call result_writeDataset(stt%v_edg_pos,group,trim(prm%output(ou)), & From a603e153dbc25f87977e24eb14d35b2b47cdad8c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 22 Dec 2023 23:29:46 +0100 Subject: [PATCH 16/33] simplified --- src/phase_mechanical_plastic_nonlocal.f90 | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 361ecc98d..6f7be8039 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -477,6 +477,7 @@ module function plastic_nonlocal_init() result(myPlasticity) extmsg = trim(extmsg)//' atol_gamma' stt%v => plasticState(ph)%state (11*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nmembers) + st0%v => plasticState(ph)%state0 (11*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nmembers) stt%v_edg_pos => plasticState(ph)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nmembers) stt%v_edg_neg => plasticState(ph)%state (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nmembers) stt%v_scr_pos => plasticState(ph)%state (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nmembers) @@ -860,13 +861,13 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en) deltaDUpper ! change in maximum stable dipole distance for edges and screws - associate(prm => param(ph),dst => dependentState(ph),del => deltaState(ph)) + associate(prm => param(ph),dst => dependentState(ph),del => deltaState(ph), stt=>state(ph)) mu = elastic_mu(ph,en,prm%isotropic_bound) nu = elastic_nu(ph,en,prm%isotropic_bound) !*** shortcut to state variables - forall (s = 1:prm%sum_N_sl, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) + v = reshape(stt%v(:,en),[prm%sum_N_sl,4]) forall (s = 1:prm%sum_N_sl, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ph),en) rho = getRho(ph,en) @@ -974,7 +975,8 @@ module subroutine nonlocal_dotState(Mp,timestep, & return end if - associate(prm => param(ph), dst => dependentState(ph), dot => dotState(ph), stt => state(ph)) + associate(prm => param(ph), dst => dependentState(ph), dot => dotState(ph), & + stt => state(ph), st0 => state0(ph)) mu = elastic_mu(ph,en,prm%isotropic_bound) nu = elastic_nu(ph,en,prm%isotropic_bound) @@ -989,11 +991,10 @@ module subroutine nonlocal_dotState(Mp,timestep, & rho0 = getRho0(ph,en) my_rhoSgl0 = rho0(:,sgl) - forall (s = 1:prm%sum_N_sl, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) + v = reshape(stt%v(:,en),[prm%sum_N_sl,4]) dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) - ! limits for stable dipole height do s = 1,prm%sum_N_sl tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) + dst%tau_back(s,en) @@ -1030,7 +1031,7 @@ module subroutine nonlocal_dotState(Mp,timestep, & * sqrt(dst%rho_forest(:,en)) / prm%i_sl / prm%b_sl, 2, 4) ! eq. 3.26 end if isBCC - forall (s = 1:prm%sum_N_sl, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en) + v0 = reshape(st0%v(:,en),[prm%sum_N_sl,4]) !**************************************************************************** @@ -1170,7 +1171,8 @@ function rhoDotFlux(timestep,ph,en) associate(prm => param(ph), & dst => dependentState(ph), & - stt => state(ph)) + stt => state(ph), & + st0 => state0(ph)) ns = prm%sum_N_sl dot_gamma = 0.0_pREAL @@ -1180,10 +1182,10 @@ function rhoDotFlux(timestep,ph,en) rho0 = getRho0(ph,en) my_rhoSgl0 = rho0(:,sgl) - forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) !ToDo: MD: I think we should use state0 here + v = reshape(stt%v(:,en),[prm%sum_N_sl,4]) !ToDo: MD: I think we should use state0 here dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) - forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en) + v0 = reshape(st0%v(:,en),[prm%sum_N_sl,4]) !**************************************************************************** !*** calculate dislocation fluxes (only for nonlocal plasticity) From e7543fd71586b7499b5f0fc7fbcf37042e380d59 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 22 Dec 2023 23:45:34 +0100 Subject: [PATCH 17/33] can be stored as dependent state --- src/phase_mechanical_plastic_nonlocal.f90 | 23 ++++++++--------------- 1 file changed, 8 insertions(+), 15 deletions(-) diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 6f7be8039..af1108405 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -44,8 +44,7 @@ submodule(phase:plastic) nonlocal ! BEGIN DEPRECATED integer, dimension(:,:,:), allocatable :: & iRhoU, & !< state indices for unblocked density - iV, & !< state indices for dislocation velocities - iD !< state indices for stable dipole height + iV !< state indices for dislocation velocities !END DEPRECATED real(pREAL), dimension(:,:,:,:,:,:), allocatable :: & @@ -125,7 +124,8 @@ submodule(phase:plastic) nonlocal real(pREAL), allocatable, dimension(:,:) :: & tau_pass, & tau_back, & - rho_forest + rho_forest, & + max_dipole_height real(pREAL), allocatable, dimension(:,:,:,:,:) :: & compatibility end type tNonlocalDependentState @@ -391,8 +391,7 @@ module function plastic_nonlocal_init() result(myPlasticity) 'gamma ' ]) * prm%sum_N_sl !< "basic" microstructural state variables that are independent from other state variables sizeState = sizeDotState & + size([ 'velocityEdgePos ','velocityEdgeNeg ', & - 'velocityScrewPos ','velocityScrewNeg ', & - 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure + 'velocityScrewPos ','velocityScrewNeg ']) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure sizeDeltaState = sizeDotState call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState,0) ! ToDo: state structure does not follow convention @@ -486,6 +485,7 @@ module function plastic_nonlocal_init() result(myPlasticity) allocate(dst%tau_pass(prm%sum_N_sl,Nmembers),source=0.0_pREAL) allocate(dst%tau_back(prm%sum_N_sl,Nmembers),source=0.0_pREAL) allocate(dst%rho_forest(prm%sum_N_sl,Nmembers),source=0.0_pREAL) + allocate(dst%max_dipole_height(2*prm%sum_N_sl,Nmembers),source=0.0_pREAL) ! edge and screw allocate(dst%compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,Nmembers),source=0.0_pREAL) end associate @@ -503,7 +503,6 @@ module function plastic_nonlocal_init() result(myPlasticity) ! BEGIN DEPRECATED---------------------------------------------------------------------------------- allocate(iRhoU(maxval(param%sum_N_sl),4,phases%length), source=0) allocate(iV(maxval(param%sum_N_sl),4,phases%length), source=0) - allocate(iD(maxval(param%sum_N_sl),2,phases%length), source=0) do ph = 1, phases%length @@ -525,13 +524,7 @@ module function plastic_nonlocal_init() result(myPlasticity) iV(s,t,ph) = l end do end do - do t = 1,2 - do s = 1,param(ph)%sum_N_sl - l = l + 1 - iD(s,t,ph) = l - end do - end do - if (iD(param(ph)%sum_N_sl,2,ph) /= plasticState(ph)%sizeState) & + if (iV(param(ph)%sum_N_sl,4,ph) /= plasticState(ph)%sizeState) & error stop 'state indices not properly set (nonlocal)' end do @@ -868,7 +861,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en) !*** shortcut to state variables v = reshape(stt%v(:,en),[prm%sum_N_sl,4]) - forall (s = 1:prm%sum_N_sl, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ph),en) + dUpperOld = reshape(dst%max_dipole_height(:,en),[prm%sum_N_sl,2]) rho = getRho(ph,en) rhoDip = rho(:,dip) @@ -915,7 +908,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en) / (dUpperOld(s,c) - prm%minDipoleHeight(s,c)) forall (t=1:4) deltaRhoDipole2SingleStress(:,t) = -0.5_pREAL * deltaRhoDipole2SingleStress(:,(t-1)/2+9) - forall (s = 1:prm%sum_N_sl, c = 1:2) plasticState(ph)%state(iD(s,c,ph),en) = dUpper(s,c) + dst%max_dipole_height(:,en) = pack(dUpper,.true.) plasticState(ph)%deltaState(:,en) = 0.0_pREAL del%rho(:,en) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*prm%sum_N_sl]) From a8e979e904396f44780f2917126d16010820b5cf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 23 Dec 2023 07:04:04 +0100 Subject: [PATCH 18/33] ensure that files are closed automatically reported by Karo Sedighiani --- python/damask/_colormap.py | 43 ++++++++------------------ python/damask/_loadcasegrid.py | 12 ++++---- python/damask/_table.py | 56 ++++++++++++++++------------------ python/damask/_yaml.py | 16 +++++----- python/damask/util.py | 21 +++++++++---- 5 files changed, 69 insertions(+), 79 deletions(-) diff --git a/python/damask/_colormap.py b/python/damask/_colormap.py index f4834da78..8eb20e9be 100644 --- a/python/damask/_colormap.py +++ b/python/damask/_colormap.py @@ -2,7 +2,7 @@ import os import json import functools import colorsys -from typing import Optional, Union, TextIO +from typing import Optional, Union from itertools import chain import numpy as np @@ -344,30 +344,6 @@ class Colormap(mpl.colors.ListedColormap): return Colormap(np.array(rev.colors),rev.name[:-4] if rev.name.endswith('_r_r') else rev.name) - def _get_file_handle(self, - fname: Union[FileHandle, None], - suffix: str = '') -> TextIO: - """ - Provide file handle. - - Parameters - ---------- - fname : file, str, pathlib.Path, or None - Name or handle of file. - If None, colormap name + suffix. - suffix: str, optional - Extension to use for colormap file. - Defaults to empty. - - Returns - ------- - f : file object - File handle with write access. - - """ - return util.open_text(self.name.replace(' ','_')+suffix if fname is None else fname, 'w') - - def save_paraview(self, fname: Optional[FileHandle] = None): """ @@ -387,9 +363,9 @@ class Colormap(mpl.colors.ListedColormap): 'RGBPoints':list(chain.from_iterable([(i,*c) for i,c in enumerate(self.colors.round(6))])) }] - fhandle = self._get_file_handle(fname,'.json') - json.dump(out,fhandle,indent=4) - fhandle.write('\n') + with util.open_text(self.name.replace(' ','_')+'.json' if fname is None else fname, 'w') as fhandle: + json.dump(out,fhandle,indent=4) + fhandle.write('\n') def save_ASCII(self, @@ -405,7 +381,9 @@ class Colormap(mpl.colors.ListedColormap): """ labels = {'RGBA':4} if self.colors.shape[1] == 4 else {'RGB': 3} t = Table(labels,self.colors,[f'Creator: {util.execution_stamp("Colormap")}']) - t.save(self._get_file_handle(fname,'.txt')) + + with util.open_text(self.name.replace(' ','_')+'.txt' if fname is None else fname, 'w') as fhandle: + t.save(fhandle) def save_GOM(self, fname: Optional[FileHandle] = None): @@ -425,7 +403,8 @@ class Colormap(mpl.colors.ListedColormap): + ' '.join([f' 0 {c[0]} {c[1]} {c[2]} 255 1' for c in reversed((self.colors*255).astype(np.int64))]) \ + '\n' - self._get_file_handle(fname,'.legend').write(GOM_str) + with util.open_text(self.name.replace(' ','_')+'.legend' if fname is None else fname, 'w') as fhandle: + fhandle.write(GOM_str) def save_gmsh(self, @@ -443,7 +422,9 @@ class Colormap(mpl.colors.ListedColormap): gmsh_str = 'View.ColorTable = {\n' \ +'\n'.join([f'{c[0]},{c[1]},{c[2]},' for c in self.colors[:,:3]*255]) \ +'\n}\n' - self._get_file_handle(fname,'.msh').write(gmsh_str) + + with util.open_text(self.name.replace(' ','_')+'.msh' if fname is None else fname, 'w') as fhandle: + fhandle.write(gmsh_str) @staticmethod diff --git a/python/damask/_loadcasegrid.py b/python/damask/_loadcasegrid.py index 3dcffefcb..6b247f97c 100644 --- a/python/damask/_loadcasegrid.py +++ b/python/damask/_loadcasegrid.py @@ -70,9 +70,9 @@ class LoadcaseGrid(YAML): if key not in kwargs: kwargs[key] = default - fhandle = util.open_text(fname,'w') - try: - fhandle.write(yaml.dump(self,Dumper=MaskedMatrixDumper,**kwargs)) - except TypeError: # compatibility with old pyyaml - del kwargs['sort_keys'] - fhandle.write(yaml.dump(self,Dumper=MaskedMatrixDumper,**kwargs)) + with util.open_text(fname,'w') as fhandle: + try: + fhandle.write(yaml.dump(self,Dumper=MaskedMatrixDumper,**kwargs)) + except TypeError: # compatibility with old pyyaml + del kwargs['sort_keys'] + fhandle.write(yaml.dump(self,Dumper=MaskedMatrixDumper,**kwargs)) diff --git a/python/damask/_table.py b/python/damask/_table.py index 9608f0d13..f22ebd9e3 100644 --- a/python/damask/_table.py +++ b/python/damask/_table.py @@ -277,28 +277,28 @@ class Table: Table data from file. """ - f = util.open_text(fname) - f.seek(0) + with util.open_text(fname) as f: + f.seek(0) - comments = [] - while (line := f.readline().strip()).startswith('#'): - comments.append(line.lstrip('#').strip()) - labels = line.split() + comments = [] + while (line := f.readline().strip()).startswith('#'): + comments.append(line.lstrip('#').strip()) + labels = line.split() - shapes = {} - for label in labels: - tensor_column = re.search(r'[0-9,x]*?:[0-9]*?_',label) - if tensor_column: - my_shape = tensor_column.group().split(':',1)[0].split('x') - shapes[label.split('_',1)[1]] = tuple([int(d) for d in my_shape]) - else: - vector_column = re.match(r'[0-9]*?_',label) - if vector_column: - shapes[label.split('_',1)[1]] = (int(label.split('_',1)[0]),) + shapes = {} + for label in labels: + tensor_column = re.search(r'[0-9,x]*?:[0-9]*?_',label) + if tensor_column: + my_shape = tensor_column.group().split(':',1)[0].split('x') + shapes[label.split('_',1)[1]] = tuple([int(d) for d in my_shape]) else: - shapes[label] = (1,) + vector_column = re.match(r'[0-9]*?_',label) + if vector_column: + shapes[label.split('_',1)[1]] = (int(label.split('_',1)[0]),) + else: + shapes[label] = (1,) - data = pd.read_csv(f,names=list(range(len(labels))),sep=r'\s+') + data = pd.read_csv(f,names=list(range(len(labels))),sep=r'\s+') return Table(shapes,data,comments) @@ -339,10 +339,9 @@ class Table: Table data from file. """ - f = util.open_text(fname) - f.seek(0) - - content = f.readlines() + with util.open_text(fname) as f: + f.seek(0) + content = f.readlines() comments = [util.execution_stamp('Table','from_ang')] for line in content: @@ -605,10 +604,9 @@ class Table: labels += [f'{util.srepr(self.shapes[l],"x")}:{i+1}_{l}' \ for i in range(np.prod(self.shapes[l]))] - f = util.open_text(fname,'w') - - f.write('\n'.join([f'# {c}' for c in self.comments] + [' '.join(labels)])+('\n' if labels else '')) - try: # backward compatibility - self.data.to_csv(f,sep=' ',na_rep='nan',index=False,header=False,lineterminator='\n') - except TypeError: - self.data.to_csv(f,sep=' ',na_rep='nan',index=False,header=False,line_terminator='\n') + with util.open_text(fname,'w') as f: + f.write('\n'.join([f'# {c}' for c in self.comments] + [' '.join(labels)])+('\n' if labels else '')) + try: # backward compatibility + self.data.to_csv(f,sep=' ',na_rep='nan',index=False,header=False,lineterminator='\n') + except TypeError: + self.data.to_csv(f,sep=' ',na_rep='nan',index=False,header=False,line_terminator='\n') diff --git a/python/damask/_yaml.py b/python/damask/_yaml.py index 077f8738f..b9d10ce36 100644 --- a/python/damask/_yaml.py +++ b/python/damask/_yaml.py @@ -197,7 +197,9 @@ class YAML(dict): YAML from file. """ - return cls(yaml.load(util.open_text(fname), Loader=SafeLoader)) + with util.open_text(fname) as fhandle: + return cls(yaml.load(fhandle, Loader=SafeLoader)) + def save(self, fname: FileHandle, @@ -220,12 +222,12 @@ class YAML(dict): if 'sort_keys' not in kwargs: kwargs['sort_keys'] = False - fhandle = util.open_text(fname,'w') - try: - fhandle.write(yaml.dump(self,Dumper=NiceDumper,**kwargs)) - except TypeError: # compatibility with old pyyaml - del kwargs['sort_keys'] - fhandle.write(yaml.dump(self,Dumper=NiceDumper,**kwargs)) + with util.open_text(fname,'w') as fhandle: + try: + fhandle.write(yaml.dump(self,Dumper=NiceDumper,**kwargs)) + except TypeError: # compatibility with old pyyaml + del kwargs['sort_keys'] + fhandle.write(yaml.dump(self,Dumper=NiceDumper,**kwargs)) @property diff --git a/python/damask/util.py b/python/damask/util.py index 0688fca80..9b606de7a 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -8,12 +8,13 @@ import shlex as _shlex import re as _re import signal as _signal import fractions as _fractions +import contextlib as _contextlib from collections import abc as _abc, OrderedDict as _OrderedDict from functools import reduce as _reduce, partial as _partial, wraps as _wraps import inspect from typing import Optional as _Optional, Callable as _Callable, Union as _Union, Iterable as _Iterable, \ Dict as _Dict, List as _List, Tuple as _Tuple, Literal as _Literal, \ - Any as _Any, TextIO as _TextIO + Any as _Any, TextIO as _TextIO, Generator as _Generator from pathlib import Path as _Path import numpy as _np @@ -193,11 +194,15 @@ def run(cmd: str, return stdout, stderr - +@_contextlib.contextmanager def open_text(fname: _FileHandle, - mode: _Literal['r','w'] = 'r') -> _TextIO: # noqa + mode: _Literal['r','w'] = 'r') -> _Generator[_TextIO, None, None]: # noqa """ - Open a text file. + Open a text file with Unix line endings + + If a path or string is given, a context manager ensures that + the file handle is closed. + If a file handle is given, it remains unmodified. Parameters ---------- @@ -211,8 +216,12 @@ def open_text(fname: _FileHandle, f : file handle """ - return fname if not isinstance(fname, (str,_Path)) else \ - open(_Path(fname).expanduser(),mode,newline=('\n' if mode == 'w' else None)) + if isinstance(fname, (str,_Path)): + fhandle = open(_Path(fname).expanduser(),mode,newline=('\n' if mode == 'w' else None)) + yield fhandle + fhandle.close() + else: + yield fname def execution_stamp(class_name: str, From 14d73fc528d05113399e9197f3a41e891a58c053 Mon Sep 17 00:00:00 2001 From: Test User Date: Sat, 23 Dec 2023 18:20:44 +0100 Subject: [PATCH 19/33] [skip ci] updated version information after successful test of v3.0.0-alpha8-189-g598b0b238 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 007956b90..118dd277a 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -3.0.0-alpha8-186-gaf5bbed00 +3.0.0-alpha8-189-g598b0b238 From 1bf16951ee46e4f4e81ddd81976bdb43e9d26384 Mon Sep 17 00:00:00 2001 From: Test User Date: Sat, 23 Dec 2023 21:38:33 +0100 Subject: [PATCH 20/33] [skip ci] updated version information after successful test of v3.0.0-alpha8-196-gaed2643af --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 118dd277a..7c220c457 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -3.0.0-alpha8-189-g598b0b238 +3.0.0-alpha8-196-gaed2643af From 169bc977d7061b75b4ec30a082cb7c50164b3dda Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 23 Dec 2023 22:39:23 +0100 Subject: [PATCH 21/33] was a 'write-only' variable --- src/Marc/materialpoint_Marc.f90 | 1 - src/grid/grid_mech_FEM.f90 | 2 -- src/grid/grid_mech_spectral_basic.f90 | 3 --- src/grid/grid_mech_spectral_polarization.f90 | 2 -- src/grid/spectral_utilities.f90 | 2 +- src/homogenization.f90 | 1 - src/homogenization_mechanical.f90 | 3 +-- src/mesh/mesh_mech_FEM.f90 | 1 - 8 files changed, 2 insertions(+), 13 deletions(-) diff --git a/src/Marc/materialpoint_Marc.f90 b/src/Marc/materialpoint_Marc.f90 index 0e3835924..bce0e6db6 100644 --- a/src/Marc/materialpoint_Marc.f90 +++ b/src/Marc/materialpoint_Marc.f90 @@ -142,7 +142,6 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, if (iand(mode, materialpoint_AGERESULTS) /= 0) call materialpoint_forward - homogenization_F0(1:3,1:3,ce) = ffn homogenization_F(1:3,1:3,ce) = ffn1 if (iand(mode, materialpoint_CALCRESULTS) /= 0) then diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 754942ae1..f7346055f 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -269,7 +269,6 @@ subroutine grid_mechanical_FEM_init(num_grid) F = spread(spread(spread(math_I3,3,cells(1)),4,cells(2)),5,cells3) end if restartRead - homogenization_F0 = reshape(F_lastInc, [3,3,product(cells(1:2))*cells3]) ! set starting condition for homogenization_mechanical_response call utilities_updateCoords(F) call utilities_constitutiveResponse(P_current,P_av,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 F, & ! target F @@ -390,7 +389,6 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai F_lastInc = F - homogenization_F0 = reshape(F, [3,3,product(cells(1:2))*cells3]) end if !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 5bcfec438..3fad9d59c 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -226,7 +226,6 @@ subroutine grid_mechanical_spectral_basic_init(num_grid) F = reshape(F_lastInc,[9,cells(1),cells(2),cells3]) end if restartRead - homogenization_F0 = reshape(F_lastInc, [3,3,product(cells(1:2))*cells3]) ! set starting condition for homogenization_mechanical_response call utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F @@ -347,8 +346,6 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_ F_lastInc,reshape(F,[3,3,cells(1),cells(2),cells3]),Delta_t_old, & rotation_BC%rotate(F_aimDot,active=.true.)) F_lastInc = reshape(F,[3,3,cells(1),cells(2),cells3]) - - homogenization_F0 = reshape(F,[3,3,product(cells(1:2))*cells3]) end if !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/grid_mech_spectral_polarization.f90 b/src/grid/grid_mech_spectral_polarization.f90 index efa1a64a9..3b7cbda14 100644 --- a/src/grid/grid_mech_spectral_polarization.f90 +++ b/src/grid/grid_mech_spectral_polarization.f90 @@ -255,7 +255,6 @@ subroutine grid_mechanical_spectral_polarization_init(num_grid) F_tau_lastInc = 2.0_pREAL*F_lastInc end if restartRead - homogenization_F0 = reshape(F_lastInc, [3,3,product(cells(1:2))*cells3]) ! set starting condition for homogenization_mechanical_response call utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F @@ -391,7 +390,6 @@ subroutine grid_mechanical_spectral_polarization_forward(cutBack,guess,Delta_t,D F_lastInc = reshape(F, [3,3,cells(1),cells(2),cells3]) F_tau_lastInc = reshape(F_tau,[3,3,cells(1),cells(2),cells3]) - homogenization_F0 = reshape(F,[3,3,product(cells(1:2))*cells3]) end if !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 9ed5cbe17..19f497d57 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -756,7 +756,7 @@ end function utilities_vectorDivergence !-------------------------------------------------------------------------------------------------- -!> @brief calculate constitutive response from homogenization_F0 to F during Delta_t +!> @brief Calculate constitutive response. !-------------------------------------------------------------------------------------------------- subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& F,Delta_t,rotation_BC) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index ebb18b232..366e38053 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -52,7 +52,6 @@ module homogenization !-------------------------------------------------------------------------------------------------- ! General variables for the homogenization at a material point real(pREAL), dimension(:,:,:), allocatable, public :: & - homogenization_F0, & !< def grad of IP at start of FE increment homogenization_F !< def grad of IP to be reached at end of FE increment real(pREAL), dimension(:,:,:), allocatable, public :: & !, protected :: & Issue with ifort homogenization_P !< first P--K stress of IP diff --git a/src/homogenization_mechanical.f90 b/src/homogenization_mechanical.f90 index 31bd42aa5..6a0d916fb 100644 --- a/src/homogenization_mechanical.f90 +++ b/src/homogenization_mechanical.f90 @@ -77,8 +77,7 @@ module subroutine mechanical_init() call parseMechanical() allocate(homogenization_dPdF(3,3,3,3,discretization_Ncells), source=0.0_pREAL) - homogenization_F0 = spread(math_I3,3,discretization_Ncells) - homogenization_F = homogenization_F0 + homogenization_F = spread(math_I3,3,discretization_Ncells) allocate(homogenization_P(3,3,discretization_Ncells),source=0.0_pREAL) if (any(mechanical_type == MECHANICAL_PASS_ID)) call pass_init() diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 824ebecdd..f3c08c16c 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -686,7 +686,6 @@ subroutine FEM_mechanical_forward(guess,Delta_t,Delta_t_prev,mechBC) ! forward last inc if (guess .and. .not. cutBack) then ForwardData = .True. - homogenization_F0 = homogenization_F call SNESGetDM(mechanical_snes,dm_local,err_PETSc) !< retrieve mesh info from mechanical_snes into dm_local CHKERRQ(err_PETSc) call DMGetSection(dm_local,section,err_PETSc) From e1ca6bbed61a42bbd9ba0f1af02d7044eb569bad Mon Sep 17 00:00:00 2001 From: Test User Date: Sun, 24 Dec 2023 12:39:38 +0100 Subject: [PATCH 22/33] [skip ci] updated version information after successful test of v3.0.0-alpha8-199-g0f69eff16 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 7c220c457..bb39723c5 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -3.0.0-alpha8-196-gaed2643af +3.0.0-alpha8-199-g0f69eff16 From 395b6b543801a4faa36e5e40e25defb83d713ab0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 24 Dec 2023 06:37:30 +0100 Subject: [PATCH 23/33] don't rely on data stored in homogenization --- src/Marc/materialpoint_Marc.f90 | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/Marc/materialpoint_Marc.f90 b/src/Marc/materialpoint_Marc.f90 index bce0e6db6..c8db7d8df 100644 --- a/src/Marc/materialpoint_Marc.f90 +++ b/src/Marc/materialpoint_Marc.f90 @@ -30,7 +30,8 @@ module materialpoint_Marc real(pREAL), dimension (:,:,:), allocatable, private :: & materialpoint_cs !< Cauchy stress real(pREAL), dimension (:,:,:,:), allocatable, private :: & - materialpoint_dcsdE !< Cauchy stress tangent + materialpoint_dcsdE, & !< Cauchy stress tangent + materialpoint_F !< deformation gradient real(pREAL), dimension (:,:,:,:), allocatable, private :: & materialpoint_dcsdE_knownGood !< known good tangent @@ -95,6 +96,7 @@ subroutine materialpoint_init() print'(/,1x,a)', '<<<+- materialpoint init -+>>>'; flush(IO_STDOUT) + allocate(materialpoint_F( 3,3,discretization_nIPs,discretization_Nelems), source= 0.0_pREAL) allocate(materialpoint_cs( 6,discretization_nIPs,discretization_Nelems), source= 0.0_pREAL) allocate(materialpoint_dcsdE( 6,6,discretization_nIPs,discretization_Nelems), source= 0.0_pREAL) allocate(materialpoint_dcsdE_knownGood(6,6,discretization_nIPs,discretization_Nelems), source= 0.0_pREAL) @@ -140,9 +142,10 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, if (iand(mode, materialpoint_RESTOREJACOBIAN) /= 0) & materialpoint_dcsde = materialpoint_dcsde_knownGood - if (iand(mode, materialpoint_AGERESULTS) /= 0) call materialpoint_forward + if (iand(mode, materialpoint_AGERESULTS) /= 0) call materialpoint_forward() homogenization_F(1:3,1:3,ce) = ffn1 + materialpoint_F(1:3,1:3,ip,elCP) = ffn1 if (iand(mode, materialpoint_CALCRESULTS) /= 0) then @@ -167,17 +170,17 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, else terminalIllness ! translate from P to sigma - Kirchhoff = matmul(homogenization_P(1:3,1:3,ce), transpose(homogenization_F(1:3,1:3,ce))) - J_inverse = 1.0_pREAL / math_det33(homogenization_F(1:3,1:3,ce)) + Kirchhoff = matmul(homogenization_P(1:3,1:3,ce), transpose(materialpoint_F(1:3,1:3,ip,elCP))) + J_inverse = 1.0_pREAL / math_det33(materialpoint_F(1:3,1:3,ip,elCP)) materialpoint_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.) ! translate from dP/dF to dCS/dE H = 0.0_pREAL do i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3 H(i,j,k,l) = H(i,j,k,l) & - + homogenization_F(j,m,ce) * homogenization_F(l,n,ce) & - * homogenization_dPdF(i,m,k,n,ce) & - - math_delta(j,l) * homogenization_F(i,m,ce) * homogenization_P(k,m,ce) & + + materialpoint_F(j,m,ip,elCP) * materialpoint_F(l,n,ip,elCP) & + * homogenization_dPdF(i,m,k,n,ce) & + - math_delta(j,l) * materialpoint_F(i,m,ip,elCP) * homogenization_P(k,m,ce) & + 0.5_pREAL * ( Kirchhoff(j,l)*math_delta(i,k) + Kirchhoff(i,k)*math_delta(j,l) & + Kirchhoff(j,k)*math_delta(i,l) + Kirchhoff(i,l)*math_delta(j,k)) end do; end do; end do; end do; end do; end do From 0d747ae5e4aa9685bebcea76256b78aef569f07d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 24 Dec 2023 07:31:02 +0100 Subject: [PATCH 24/33] (incomplete) split: spectral/FFT vs mech functionality --- src/grid/DAMASK_grid.f90 | 1 + src/grid/grid_damage_spectral.f90 | 12 +- src/grid/grid_mech_FEM.f90 | 1 + src/grid/grid_mech_spectral_basic.f90 | 1 + src/grid/grid_mech_spectral_polarization.f90 | 1 + src/grid/grid_mech_utilities.f90 | 253 +++++++++++++++++++ src/grid/grid_thermal_spectral.f90 | 16 +- src/grid/spectral_utilities.f90 | 217 ---------------- 8 files changed, 269 insertions(+), 233 deletions(-) create mode 100644 src/grid/grid_mech_utilities.f90 diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index cd6a60abe..3fd220ce5 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -23,6 +23,7 @@ program DAMASK_grid use materialpoint use material use spectral_utilities + use grid_mech_utilities use grid_mechanical_spectral_basic use grid_mechanical_spectral_polarization use grid_mechanical_FEM diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 648f9ebbb..cce869653 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -44,7 +44,6 @@ module grid_damage_spectral type(tNumerics) :: num - type(tSolutionParams) :: params !-------------------------------------------------------------------------------------------------- ! PETSc data SNES :: SNES_damage @@ -57,7 +56,7 @@ module grid_damage_spectral ! reference diffusion tensor, mobility etc. integer :: totalIter = 0 !< total iteration in current increment real(pREAL), dimension(3,3) :: K_ref - real(pREAL) :: mu_ref + real(pREAL) :: mu_ref, Delta_t_ public :: & grid_damage_spectral_init, & @@ -207,8 +206,7 @@ end subroutine grid_damage_spectral_init !-------------------------------------------------------------------------------------------------- function grid_damage_spectral_solution(Delta_t) result(solution) - real(pREAL), intent(in) :: & - Delta_t !< increment in time for current solution + real(pREAL), intent(in) :: Delta_t !< increment in time for current solution type(tSolutionState) :: solution PetscInt :: devNull @@ -222,7 +220,7 @@ function grid_damage_spectral_solution(Delta_t) result(solution) !-------------------------------------------------------------------------------------------------- ! set module wide availabe data - params%Delta_t = Delta_t + Delta_t_ = Delta_t call SNESSolve(SNES_damage,PETSC_NULL_VEC,phi_PETSc,err_PETSc) CHKERRQ(err_PETSc) @@ -350,12 +348,12 @@ subroutine formResidual(residual_subdomain,x_scal,r,dummy,err_PETSc) ce = 0 do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 - r(i,j,k) = params%Delta_t*(r(i,j,k) + homogenization_f_phi(phi(i,j,k),ce)) & + r(i,j,k) = Delta_t_*(r(i,j,k) + homogenization_f_phi(phi(i,j,k),ce)) & + homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi(i,j,k)) & + mu_ref*phi(i,j,k) end do; end do; end do - r = max(min(utilities_GreenConvolution(r, K_ref, mu_ref, params%Delta_t),phi_lastInc),num%phi_min) & + r = max(min(utilities_GreenConvolution(r, K_ref, mu_ref, Delta_t_),phi_lastInc),num%phi_min) & - phi end associate err_PETSc = 0 diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index f7346055f..59835f250 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -23,6 +23,7 @@ module grid_mechanical_FEM use math use rotations use spectral_utilities + use grid_mech_utilities use config use homogenization use discretization diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 3fad9d59c..0e8ba4841 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -23,6 +23,7 @@ module grid_mechanical_spectral_basic use math use rotations use spectral_utilities + use grid_mech_utilities use homogenization use discretization_grid diff --git a/src/grid/grid_mech_spectral_polarization.f90 b/src/grid/grid_mech_spectral_polarization.f90 index 3b7cbda14..b5cc0b967 100644 --- a/src/grid/grid_mech_spectral_polarization.f90 +++ b/src/grid/grid_mech_spectral_polarization.f90 @@ -23,6 +23,7 @@ module grid_mechanical_spectral_polarization use math use rotations use spectral_utilities + use grid_mech_utilities use config use homogenization use discretization_grid diff --git a/src/grid/grid_mech_utilities.f90 b/src/grid/grid_mech_utilities.f90 new file mode 100644 index 000000000..0ea269319 --- /dev/null +++ b/src/grid/grid_mech_utilities.f90 @@ -0,0 +1,253 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, KU Leuven +!> @brief Utilities used by the mech grid solver variants +!-------------------------------------------------------------------------------------------------- +module grid_mech_utilities +#include + use PETScSys +#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) + use MPI_f08 +#endif + + use prec + use parallelization + use math + use rotations + use IO + use discretization_grid + use discretization + use spectral_utilities + use homogenization + + +#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) + implicit none(type,external) +#else + implicit none +#endif + private + +!-------------------------------------------------------------------------------------------------- +! derived types + type, public :: tBoundaryCondition !< set of parameters defining a boundary condition + real(pREAL), dimension(3,3) :: values = 0.0_pREAL + logical, dimension(3,3) :: mask = .true. + character(len=:), allocatable :: myType + end type tBoundaryCondition + + type, public :: tSolutionParams + real(pREAL), dimension(3,3) :: stress_BC + logical, dimension(3,3) :: stress_mask + type(tRotation) :: rotation_BC + real(pREAL) :: Delta_t + end type tSolutionParams + + public :: & + utilities_maskedCompliance, & + utilities_constitutiveResponse, & + utilities_calculateRate, & + utilities_forwardTensorField + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief Calculate masked compliance tensor used to adjust F to fullfill stress BC. +!-------------------------------------------------------------------------------------------------- +function utilities_maskedCompliance(rot_BC,mask_stress,C) + + real(pREAL), dimension(3,3,3,3) :: utilities_maskedCompliance !< masked compliance + real(pREAL), intent(in), dimension(3,3,3,3) :: C !< current average stiffness + type(tRotation), intent(in) :: rot_BC !< rotation of load frame + logical, intent(in), dimension(3,3) :: mask_stress !< mask of stress BC + + integer :: i, j + logical, dimension(9) :: mask_stressVector + logical, dimension(9,9) :: mask + real(pREAL), dimension(9,9) :: temp99_real + integer :: size_reduced = 0 + real(pREAL), dimension(:,:), allocatable :: & + s_reduced, & !< reduced compliance matrix (depending on number of stress BC) + c_reduced, & !< reduced stiffness (depending on number of stress BC) + sTimesC !< temp variable to check inversion + logical :: errmatinv + character(len=pSTRLEN):: formatString + + + mask_stressVector = .not. reshape(transpose(mask_stress), [9]) + size_reduced = count(mask_stressVector) + if (size_reduced > 0) then + temp99_real = math_3333to99(rot_BC%rotate(C)) + + do i = 1,9; do j = 1,9 + mask(i,j) = mask_stressVector(i) .and. mask_stressVector(j) + end do; end do + c_reduced = reshape(pack(temp99_Real,mask),[size_reduced,size_reduced]) + + allocate(s_reduced,mold = c_reduced) + call math_invert(s_reduced, errmatinv, c_reduced) ! invert reduced stiffness + if (any(IEEE_is_NaN(s_reduced))) errmatinv = .true. + +!-------------------------------------------------------------------------------------------------- +! check if inversion was successful + sTimesC = matmul(c_reduced,s_reduced) + errmatinv = errmatinv .or. any(dNeq(sTimesC,math_eye(size_reduced),1.0e-12_pREAL)) + if (errmatinv) then + write(formatString, '(i2)') size_reduced + formatString = '(/,1x,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))' + print trim(formatString), 'C * S (load) ', transpose(matmul(c_reduced,s_reduced)) + print trim(formatString), 'C (load) ', transpose(c_reduced) + print trim(formatString), 'S (load) ', transpose(s_reduced) + if (errmatinv) error stop 'matrix inversion error' + end if + temp99_real = reshape(unpack(reshape(s_reduced,[size_reduced**2]),reshape(mask,[81]),0.0_pREAL),[9,9]) + else + temp99_real = 0.0_pREAL + end if + + utilities_maskedCompliance = math_99to3333(temp99_Real) + +end function utilities_maskedCompliance + + +!-------------------------------------------------------------------------------------------------- +!> @brief Calculate constitutive response. +!-------------------------------------------------------------------------------------------------- +subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& + F,Delta_t,rotation_BC) + + real(pREAL), intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness + real(pREAL), intent(out), dimension(3,3) :: P_av !< average PK stress + real(pREAL), intent(out), dimension(3,3,cells(1),cells(2),cells3) :: P !< PK stress + real(pREAL), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: F !< deformation gradient target + real(pREAL), intent(in) :: Delta_t !< loading time + type(tRotation), intent(in), optional :: rotation_BC !< rotation of load frame + + + integer :: i + integer(MPI_INTEGER_KIND) :: err_MPI + real(pREAL), dimension(3,3,3,3) :: dPdF_max, dPdF_min + real(pREAL) :: dPdF_norm_max, dPdF_norm_min + real(pREAL), dimension(2) :: valueAndRank !< pair of min/max norm of dPdF to synchronize min/max of dPdF + + print'(/,1x,a)', '... evaluating constitutive response ......................................' + flush(IO_STDOUT) + + homogenization_F = reshape(F,[3,3,product(cells(1:2))*cells3]) ! set materialpoint target F to estimated field + + call homogenization_mechanical_response(Delta_t,1,product(cells(1:2))*cells3) ! calculate P field + if (.not. terminallyIll) & + call homogenization_thermal_response(Delta_t,1,product(cells(1:2))*cells3) + if (.not. terminallyIll) & + call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(cells(1:2))*cells3]) + + P = reshape(homogenization_P, [3,3,cells(1),cells(2),cells3]) + P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt + call MPI_Allreduce(MPI_IN_PLACE,P_av,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) + if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + if (present(rotation_BC)) then + if (any(dNeq(rotation_BC%asQuaternion(), real([1.0, 0.0, 0.0, 0.0],pREAL)))) & + print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & + 'Piola--Kirchhoff stress (lab) / MPa =', transpose(P_av)*1.e-6_pREAL + P_av = rotation_BC%rotate(P_av) + end if + print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & + 'Piola--Kirchhoff stress / MPa =', transpose(P_av)*1.e-6_pREAL + flush(IO_STDOUT) + + dPdF_max = 0.0_pREAL + dPdF_norm_max = 0.0_pREAL + dPdF_min = huge(1.0_pREAL) + dPdF_norm_min = huge(1.0_pREAL) + do i = 1, product(cells(1:2))*cells3 + if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then + dPdF_max = homogenization_dPdF(1:3,1:3,1:3,1:3,i) + dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2) + end if + if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then + dPdF_min = homogenization_dPdF(1:3,1:3,1:3,1:3,i) + dPdF_norm_min = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2) + end if + end do + + valueAndRank = [dPdF_norm_max,real(worldrank,pREAL)] + call MPI_Allreduce(MPI_IN_PLACE,valueAndRank,1_MPI_INTEGER_KIND,MPI_2DOUBLE_PRECISION,MPI_MAXLOC,MPI_COMM_WORLD,err_MPI) + if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call MPI_Bcast(dPdF_max,81_MPI_INTEGER_KIND,MPI_DOUBLE,int(valueAndRank(2),MPI_INTEGER_KIND),MPI_COMM_WORLD,err_MPI) + if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + + valueAndRank = [dPdF_norm_min,real(worldrank,pREAL)] + call MPI_Allreduce(MPI_IN_PLACE,valueAndRank,1_MPI_INTEGER_KIND,MPI_2DOUBLE_PRECISION,MPI_MINLOC,MPI_COMM_WORLD,err_MPI) + if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call MPI_Bcast(dPdF_min,81_MPI_INTEGER_KIND,MPI_DOUBLE,int(valueAndRank(2),MPI_INTEGER_KIND),MPI_COMM_WORLD,err_MPI) + if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + + C_minmaxAvg = 0.5_pREAL*(dPdF_max + dPdF_min) + + C_volAvg = sum(homogenization_dPdF,dim=5) + call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) + if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + C_volAvg = C_volAvg * wgt + + +end subroutine utilities_constitutiveResponse + + +!-------------------------------------------------------------------------------------------------- +!> @brief Calculate forward rate, either as local guess or as homogeneous add on. +!-------------------------------------------------------------------------------------------------- +pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate) + + real(pREAL), intent(in), dimension(3,3) :: & + avRate !< homogeneous addon + real(pREAL), intent(in) :: & + dt !< Delta_t between field0 and field + logical, intent(in) :: & + heterogeneous !< calculate field of rates + real(pREAL), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: & + field0, & !< data of previous step + field !< data of current step + real(pREAL), dimension(3,3,cells(1),cells(2),cells3) :: & + utilities_calculateRate + + + utilities_calculateRate = merge((field-field0) / dt, & + spread(spread(spread(avRate,3,cells(1)),4,cells(2)),5,cells3), & + heterogeneous) + +end function utilities_calculateRate + + +!-------------------------------------------------------------------------------------------------- +!> @brief forwards a field with a pointwise given rate, if aim is given, +!> ensures that the average matches the aim +!-------------------------------------------------------------------------------------------------- +function utilities_forwardTensorField(Delta_t,field_lastInc,rate,aim) + + real(pREAL), intent(in) :: & + Delta_t !< Delta_t of current step + real(pREAL), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: & + field_lastInc, & !< initial field + rate !< rate by which to forward + real(pREAL), intent(in), optional, dimension(3,3) :: & + aim !< average field value aim + + real(pREAL), dimension(3,3,cells(1),cells(2),cells3) :: & + utilities_forwardTensorField + real(pREAL), dimension(3,3) :: fieldDiff !< - aim + integer(MPI_INTEGER_KIND) :: err_MPI + + + utilities_forwardTensorField = field_lastInc + rate*Delta_t + if (present(aim)) then !< correct to match average + fieldDiff = sum(sum(sum(utilities_forwardTensorField,dim=5),dim=4),dim=3)*wgt + call MPI_Allreduce(MPI_IN_PLACE,fieldDiff,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) + if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + fieldDiff = fieldDiff - aim + utilities_forwardTensorField = utilities_forwardTensorField & + - spread(spread(spread(fieldDiff,3,cells(1)),4,cells(2)),5,cells3) + end if + +end function utilities_forwardTensorField + +end module grid_mech_utilities diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index c178c5810..c8e638207 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -43,7 +43,6 @@ module grid_thermal_spectral type(tNumerics) :: num - type(tSolutionParams) :: params !-------------------------------------------------------------------------------------------------- ! PETSc data SNES :: SNES_thermal @@ -56,7 +55,7 @@ module grid_thermal_spectral ! reference diffusion tensor, mobility etc. integer :: totalIter = 0 !< total iteration in current increment real(pREAL), dimension(3,3) :: K_ref - real(pREAL) :: mu_ref + real(pREAL) :: mu_ref, Delta_t_ public :: & grid_thermal_spectral_init, & @@ -186,8 +185,7 @@ end subroutine grid_thermal_spectral_init !-------------------------------------------------------------------------------------------------- function grid_thermal_spectral_solution(Delta_t) result(solution) - real(pREAL), intent(in) :: & - Delta_t !< increment in time for current solution + real(pREAL), intent(in) :: Delta_t !< increment in time for current solution type(tSolutionState) :: solution PetscInt :: devNull @@ -201,7 +199,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution) !-------------------------------------------------------------------------------------------------- ! set module wide availabe data - params%Delta_t = Delta_t + Delta_t_ = Delta_t call SNESSolve(SNES_thermal,PETSC_NULL_VEC,T_PETSc,err_PETSc) CHKERRQ(err_PETSc) @@ -227,7 +225,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution) T_stagInc = T call homogenization_thermal_setField(reshape(T,[product(cells(1:2))*cells3]), & - reshape(T-T_lastInc,[product(cells(1:2))*cells3])/params%Delta_t) + reshape(T-T_lastInc,[product(cells(1:2))*cells3])/Delta_t_) call DMDAVecRestoreArrayF90(DM_thermal,T_PETSc,T,err_PETSc) CHKERRQ(err_PETSc) @@ -264,7 +262,7 @@ subroutine grid_thermal_spectral_forward(cutBack) T = T_lastInc T_stagInc = T_lastInc else - dotT_lastInc = (T - T_lastInc)/params%Delta_t + dotT_lastInc = (T - T_lastInc)/Delta_t_ T_lastInc = T call updateReference() end if @@ -336,13 +334,13 @@ subroutine formResidual(residual_subdomain,x_scal,r,dummy,err_PETSc) ce = 0 do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 - r(i,j,k) = params%Delta_t*(r(i,j,k) + homogenization_f_T(ce)) & + r(i,j,k) = Delta_t_*(r(i,j,k) + homogenization_f_T(ce)) & + homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T(i,j,k)) & + mu_ref*T(i,j,k) end do; end do; end do r = T & - - utilities_GreenConvolution(r, K_ref, mu_ref, params%Delta_t) + - utilities_GreenConvolution(r, K_ref, mu_ref, Delta_t_) end associate err_PETSc = 0 diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 19f497d57..4ea53d038 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -75,19 +75,6 @@ module spectral_utilities termIll = .false. end type tSolutionState - type, public :: tBoundaryCondition !< set of parameters defining a boundary condition - real(pREAL), dimension(3,3) :: values = 0.0_pREAL - logical, dimension(3,3) :: mask = .true. - character(len=:), allocatable :: myType - end type tBoundaryCondition - - type, public :: tSolutionParams - real(pREAL), dimension(3,3) :: stress_BC - logical, dimension(3,3) :: stress_mask - type(tRotation) :: rotation_BC - real(pREAL) :: Delta_t - end type tSolutionParams - type :: tNumerics integer :: & divergence_correction !< scale divergence/curl calculation: [0: no correction, 1: size scaled to 1, 2: size scaled to Npoints] @@ -121,10 +108,6 @@ module spectral_utilities utilities_curlRMS, & utilities_scalarGradient, & utilities_vectorDivergence, & - utilities_maskedCompliance, & - utilities_constitutiveResponse, & - utilities_calculateRate, & - utilities_forwardTensorField, & utilities_updateCoords contains @@ -653,65 +636,6 @@ real(pREAL) function utilities_curlRMS(tensorField) end function utilities_curlRMS -!-------------------------------------------------------------------------------------------------- -!> @brief Calculate masked compliance tensor used to adjust F to fullfill stress BC. -!-------------------------------------------------------------------------------------------------- -function utilities_maskedCompliance(rot_BC,mask_stress,C) - - real(pREAL), dimension(3,3,3,3) :: utilities_maskedCompliance !< masked compliance - real(pREAL), intent(in), dimension(3,3,3,3) :: C !< current average stiffness - type(tRotation), intent(in) :: rot_BC !< rotation of load frame - logical, intent(in), dimension(3,3) :: mask_stress !< mask of stress BC - - integer :: i, j - logical, dimension(9) :: mask_stressVector - logical, dimension(9,9) :: mask - real(pREAL), dimension(9,9) :: temp99_real - integer :: size_reduced = 0 - real(pREAL), dimension(:,:), allocatable :: & - s_reduced, & !< reduced compliance matrix (depending on number of stress BC) - c_reduced, & !< reduced stiffness (depending on number of stress BC) - sTimesC !< temp variable to check inversion - logical :: errmatinv - character(len=pSTRLEN):: formatString - - - mask_stressVector = .not. reshape(transpose(mask_stress), [9]) - size_reduced = count(mask_stressVector) - if (size_reduced > 0) then - temp99_real = math_3333to99(rot_BC%rotate(C)) - - do i = 1,9; do j = 1,9 - mask(i,j) = mask_stressVector(i) .and. mask_stressVector(j) - end do; end do - c_reduced = reshape(pack(temp99_Real,mask),[size_reduced,size_reduced]) - - allocate(s_reduced,mold = c_reduced) - call math_invert(s_reduced, errmatinv, c_reduced) ! invert reduced stiffness - if (any(IEEE_is_NaN(s_reduced))) errmatinv = .true. - -!-------------------------------------------------------------------------------------------------- -! check if inversion was successful - sTimesC = matmul(c_reduced,s_reduced) - errmatinv = errmatinv .or. any(dNeq(sTimesC,math_eye(size_reduced),1.0e-12_pREAL)) - if (errmatinv) then - write(formatString, '(i2)') size_reduced - formatString = '(/,1x,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))' - print trim(formatString), 'C * S (load) ', transpose(matmul(c_reduced,s_reduced)) - print trim(formatString), 'C (load) ', transpose(c_reduced) - print trim(formatString), 'S (load) ', transpose(s_reduced) - if (errmatinv) error stop 'matrix inversion error' - end if - temp99_real = reshape(unpack(reshape(s_reduced,[size_reduced**2]),reshape(mask,[81]),0.0_pREAL),[9,9]) - else - temp99_real = 0.0_pREAL - end if - - utilities_maskedCompliance = math_99to3333(temp99_Real) - -end function utilities_maskedCompliance - - !-------------------------------------------------------------------------------------------------- !> @brief Calculate gradient of scalar field. !-------------------------------------------------------------------------------------------------- @@ -755,147 +679,6 @@ function utilities_vectorDivergence(field) result(div) end function utilities_vectorDivergence -!-------------------------------------------------------------------------------------------------- -!> @brief Calculate constitutive response. -!-------------------------------------------------------------------------------------------------- -subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& - F,Delta_t,rotation_BC) - - real(pREAL), intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness - real(pREAL), intent(out), dimension(3,3) :: P_av !< average PK stress - real(pREAL), intent(out), dimension(3,3,cells(1),cells(2),cells3) :: P !< PK stress - real(pREAL), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: F !< deformation gradient target - real(pREAL), intent(in) :: Delta_t !< loading time - type(tRotation), intent(in), optional :: rotation_BC !< rotation of load frame - - - integer :: i - integer(MPI_INTEGER_KIND) :: err_MPI - real(pREAL), dimension(3,3,3,3) :: dPdF_max, dPdF_min - real(pREAL) :: dPdF_norm_max, dPdF_norm_min - real(pREAL), dimension(2) :: valueAndRank !< pair of min/max norm of dPdF to synchronize min/max of dPdF - - print'(/,1x,a)', '... evaluating constitutive response ......................................' - flush(IO_STDOUT) - - homogenization_F = reshape(F,[3,3,product(cells(1:2))*cells3]) ! set materialpoint target F to estimated field - - call homogenization_mechanical_response(Delta_t,1,product(cells(1:2))*cells3) ! calculate P field - if (.not. terminallyIll) & - call homogenization_thermal_response(Delta_t,1,product(cells(1:2))*cells3) - if (.not. terminallyIll) & - call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(cells(1:2))*cells3]) - - P = reshape(homogenization_P, [3,3,cells(1),cells(2),cells3]) - P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt - call MPI_Allreduce(MPI_IN_PLACE,P_av,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - if (present(rotation_BC)) then - if (any(dNeq(rotation_BC%asQuaternion(), real([1.0, 0.0, 0.0, 0.0],pREAL)))) & - print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & - 'Piola--Kirchhoff stress (lab) / MPa =', transpose(P_av)*1.e-6_pREAL - P_av = rotation_BC%rotate(P_av) - end if - print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & - 'Piola--Kirchhoff stress / MPa =', transpose(P_av)*1.e-6_pREAL - flush(IO_STDOUT) - - dPdF_max = 0.0_pREAL - dPdF_norm_max = 0.0_pREAL - dPdF_min = huge(1.0_pREAL) - dPdF_norm_min = huge(1.0_pREAL) - do i = 1, product(cells(1:2))*cells3 - if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then - dPdF_max = homogenization_dPdF(1:3,1:3,1:3,1:3,i) - dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2) - end if - if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then - dPdF_min = homogenization_dPdF(1:3,1:3,1:3,1:3,i) - dPdF_norm_min = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2) - end if - end do - - valueAndRank = [dPdF_norm_max,real(worldrank,pREAL)] - call MPI_Allreduce(MPI_IN_PLACE,valueAndRank,1_MPI_INTEGER_KIND,MPI_2DOUBLE_PRECISION,MPI_MAXLOC,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - call MPI_Bcast(dPdF_max,81_MPI_INTEGER_KIND,MPI_DOUBLE,int(valueAndRank(2),MPI_INTEGER_KIND),MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - - valueAndRank = [dPdF_norm_min,real(worldrank,pREAL)] - call MPI_Allreduce(MPI_IN_PLACE,valueAndRank,1_MPI_INTEGER_KIND,MPI_2DOUBLE_PRECISION,MPI_MINLOC,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - call MPI_Bcast(dPdF_min,81_MPI_INTEGER_KIND,MPI_DOUBLE,int(valueAndRank(2),MPI_INTEGER_KIND),MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - - C_minmaxAvg = 0.5_pREAL*(dPdF_max + dPdF_min) - - C_volAvg = sum(homogenization_dPdF,dim=5) - call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - C_volAvg = C_volAvg * wgt - - -end subroutine utilities_constitutiveResponse - - -!-------------------------------------------------------------------------------------------------- -!> @brief Calculate forward rate, either as local guess or as homogeneous add on. -!-------------------------------------------------------------------------------------------------- -pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate) - - real(pREAL), intent(in), dimension(3,3) :: & - avRate !< homogeneous addon - real(pREAL), intent(in) :: & - dt !< Delta_t between field0 and field - logical, intent(in) :: & - heterogeneous !< calculate field of rates - real(pREAL), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: & - field0, & !< data of previous step - field !< data of current step - real(pREAL), dimension(3,3,cells(1),cells(2),cells3) :: & - utilities_calculateRate - - - utilities_calculateRate = merge((field-field0) / dt, & - spread(spread(spread(avRate,3,cells(1)),4,cells(2)),5,cells3), & - heterogeneous) - -end function utilities_calculateRate - - -!-------------------------------------------------------------------------------------------------- -!> @brief forwards a field with a pointwise given rate, if aim is given, -!> ensures that the average matches the aim -!-------------------------------------------------------------------------------------------------- -function utilities_forwardTensorField(Delta_t,field_lastInc,rate,aim) - - real(pREAL), intent(in) :: & - Delta_t !< Delta_t of current step - real(pREAL), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: & - field_lastInc, & !< initial field - rate !< rate by which to forward - real(pREAL), intent(in), optional, dimension(3,3) :: & - aim !< average field value aim - - real(pREAL), dimension(3,3,cells(1),cells(2),cells3) :: & - utilities_forwardTensorField - real(pREAL), dimension(3,3) :: fieldDiff !< - aim - integer(MPI_INTEGER_KIND) :: err_MPI - - - utilities_forwardTensorField = field_lastInc + rate*Delta_t - if (present(aim)) then !< correct to match average - fieldDiff = sum(sum(sum(utilities_forwardTensorField,dim=5),dim=4),dim=3)*wgt - call MPI_Allreduce(MPI_IN_PLACE,fieldDiff,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - fieldDiff = fieldDiff - aim - utilities_forwardTensorField = utilities_forwardTensorField & - - spread(spread(spread(fieldDiff,3,cells(1)),4,cells(2)),5,cells3) - end if - -end function utilities_forwardTensorField - - !-------------------------------------------------------------------------------------------------- !> @brief Calculate Filter for Fourier convolution. !> @details this is the full operator to calculate derivatives, i.e. 2 \pi i k for the From af7f73f3796356b987b41f14a2fed426913bbd4b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 27 Dec 2023 15:07:04 +0100 Subject: [PATCH 25/33] CMake knows compiler flags --- CMakeLists.txt | 2 ++ cmake/Compiler-GNU.cmake | 3 --- cmake/Compiler-Intel.cmake | 3 --- cmake/Compiler-IntelLLVM.cmake | 3 --- cmake/Compiler-LLVMFlang.cmake | 2 -- 5 files changed, 2 insertions(+), 11 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index b630f7e46..d0f50affb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -29,6 +29,8 @@ else() endif() add_definitions("-D${DAMASK_SOLVER}") +set(CMAKE_Fortran_PREPROCESS "ON") + # EXPERIMENTAL: This might help to detect HDF5 and FFTW3 in the future if PETSc is not aware of them set(ENV{PKG_CONFIG_PATH} "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/externalpackages:$ENV{PKG_CONFIG_PATH}") pkg_check_modules(HDF5 hdf5) diff --git a/cmake/Compiler-GNU.cmake b/cmake/Compiler-GNU.cmake index ee58f8cc7..61ec2143b 100644 --- a/cmake/Compiler-GNU.cmake +++ b/cmake/Compiler-GNU.cmake @@ -19,9 +19,6 @@ set (STANDARD_CHECK "-std=f2018 -pedantic-errors" ) #------------------------------------------------------------------------------------------------ # Fine tuning compilation options -set (COMPILE_FLAGS "${COMPILE_FLAGS} -cpp") -# preprocessor - set (COMPILE_FLAGS "${COMPILE_FLAGS} -fPIE") # position independent code diff --git a/cmake/Compiler-Intel.cmake b/cmake/Compiler-Intel.cmake index 2182459ed..59ba6d74d 100644 --- a/cmake/Compiler-Intel.cmake +++ b/cmake/Compiler-Intel.cmake @@ -22,9 +22,6 @@ set (LINKER_FLAGS "${LINKER_FLAGS} -shared-intel") #------------------------------------------------------------------------------------------------ # Fine tuning compilation options -set (COMPILE_FLAGS "${COMPILE_FLAGS} -fpp") -# preprocessor - set (COMPILE_FLAGS "${COMPILE_FLAGS} -no-ftz") # disable flush underflow to zero, will be set if -O[1,2,3] diff --git a/cmake/Compiler-IntelLLVM.cmake b/cmake/Compiler-IntelLLVM.cmake index 4cdd2588d..3749b925f 100644 --- a/cmake/Compiler-IntelLLVM.cmake +++ b/cmake/Compiler-IntelLLVM.cmake @@ -24,9 +24,6 @@ set (LINKER_FLAGS "${LINKER_FLAGS} -shared-intel -fc=ifx") #------------------------------------------------------------------------------------------------ # Fine tuning compilation options -set (COMPILE_FLAGS "${COMPILE_FLAGS} -fpp") -# preprocessor - set (COMPILE_FLAGS "${COMPILE_FLAGS} -no-ftz") # disable flush underflow to zero, will be set if -O[1,2,3] diff --git a/cmake/Compiler-LLVMFlang.cmake b/cmake/Compiler-LLVMFlang.cmake index c6b0405bc..ff0bdf6e8 100644 --- a/cmake/Compiler-LLVMFlang.cmake +++ b/cmake/Compiler-LLVMFlang.cmake @@ -7,6 +7,4 @@ endif () #------------------------------------------------------------------------------------------------ # Fine tuning compilation options -set (COMPILE_FLAGS "${COMPILE_FLAGS} -cpp") -# preprocessor From e4e0840479a7d31c47a09c6d01e04405c65ccbef Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 27 Dec 2023 16:09:40 +0100 Subject: [PATCH 26/33] simplified --- CMakeLists.txt | 9 +++------ cmake/Compiler-LLVMFlang.cmake | 2 ++ 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index d0f50affb..a8e551373 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -93,27 +93,24 @@ if(CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY") endif() -list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake) if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") - include(Compiler-GNU) set(Fortran_COMPILER_VERSION_MIN 9.1) elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "Intel") - include(Compiler-Intel) set(Fortran_COMPILER_VERSION_MIN 19) elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "IntelLLVM") - include(Compiler-IntelLLVM) set(Fortran_COMPILER_VERSION_MIN 19) elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "LLVMFlang") - include(Compiler-LLVMFlang) set(Fortran_COMPILER_VERSION_MIN 19) else() message(FATAL_ERROR "Compiler '${CMAKE_Fortran_COMPILER_ID}' not supported") endif() - if(CMAKE_Fortran_COMPILER_VERSION VERSION_LESS Fortran_COMPILER_VERSION_MIN) message(FATAL_ERROR "Version '${CMAKE_Fortran_COMPILER_VERSION}' of '${CMAKE_Fortran_COMPILER_ID}' is not supported") endif() +list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake) +include("Compiler-${CMAKE_Fortran_COMPILER_ID}") + file(STRINGS "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/petsc/conf/petscvariables" PETSC_EXTERNAL_LIB REGEX "PETSC_EXTERNAL_LIB_BASIC = .*$?") string(REPLACE "PETSC_EXTERNAL_LIB_BASIC = " "" PETSC_EXTERNAL_LIB "${PETSC_EXTERNAL_LIB}") message("PETSC_EXTERNAL_LIB:\n${PETSC_EXTERNAL_LIB}\n") diff --git a/cmake/Compiler-LLVMFlang.cmake b/cmake/Compiler-LLVMFlang.cmake index ff0bdf6e8..b28df7027 100644 --- a/cmake/Compiler-LLVMFlang.cmake +++ b/cmake/Compiler-LLVMFlang.cmake @@ -5,6 +5,8 @@ if (OPENMP) set (OPENMP_FLAGS "-fopenmp") endif () +set (STANDARD_CHECK "-std=f2018 -pedantic" ) + #------------------------------------------------------------------------------------------------ # Fine tuning compilation options From 77e1c5a8a7b81d46c09c825a6f2b002366fbbaff Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 28 Dec 2023 19:15:26 +0100 Subject: [PATCH 27/33] standard indexing: ce(ll) instead of element and integration point --- src/phase.f90 | 9 +++------ src/phase_mechanical_plastic_nonlocal.f90 | 17 +++++++++++------ 2 files changed, 14 insertions(+), 12 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index ffa606cc2..148b028a7 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -326,11 +326,8 @@ module phase real(pREAL) :: f end function phase_f_T - module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el) - integer, intent(in) :: & - ph, & - ip, & - el + module subroutine plastic_nonlocal_updateCompatibility(orientation,ce) + integer, intent(in) :: ce type(tRotationContainer), dimension(:), intent(in) :: orientation end subroutine plastic_nonlocal_updateCompatibility @@ -590,7 +587,7 @@ subroutine crystallite_orientations(co,ip,el) call phase_O(ph)%data(en)%fromMatrix(transpose(math_rotationalPart(mechanical_F_e(ph,en)))) if (plasticState(material_ID_phase(1,(el-1)*discretization_nIPs + ip))%nonlocal) & - call plastic_nonlocal_updateCompatibility(phase_O,material_ID_phase(1,(el-1)*discretization_nIPs + ip),ip,el) + call plastic_nonlocal_updateCompatibility(phase_O,(el-1)*discretization_nIPs + ip) end subroutine crystallite_orientations diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 545dec4e6..166c30ea1 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -1331,18 +1331,19 @@ end function rhoDotFlux ! plane normals and signed cosine of the angle between the slip directions. Only the largest values ! that sum up to a total of 1 are considered, all others are set to zero. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el) +module subroutine plastic_nonlocal_updateCompatibility(orientation,ce) type(tRotationContainer), dimension(:), intent(in) :: & orientation ! crystal orientation integer, intent(in) :: & - ph, & - ip, & - el + ce integer :: & n, & ! neighbor index + ph, & en, & + ip, & + el, & neighbor_e, & ! element index of my neighbor neighbor_i, & ! integration point index of my neighbor neighbor_me, & @@ -1350,17 +1351,21 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el) ns, & ! number of active slip systems s1, & ! slip system index (en) s2 ! slip system index (my neighbor) - real(pREAL), dimension(2,param(ph)%sum_N_sl,param(ph)%sum_N_sl,nIPneighbors) :: & + real(pREAL), dimension(2,param(material_ID_phase(1,ce))%sum_N_sl,param(material_ID_phase(1,ce))%sum_N_sl,nIPneighbors) :: & my_compatibility ! my_compatibility for current element and ip real(pREAL) :: & my_compatibilitySum, & thresholdValue, & nThresholdValues - logical, dimension(param(ph)%sum_N_sl) :: & + logical, dimension(param(material_ID_phase(1,ce))%sum_N_sl) :: & belowThreshold type(tRotation) :: mis + ph = material_ID_phase(1,ce) + el = (ce-1)/discretization_nIPs + 1 + ip = modulo(ce-1,discretization_nIPs) + 1 + associate(prm => param(ph)) ns = prm%sum_N_sl From aae99156b4eafdbcc5b460250a6c2d401ba86b96 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 28 Dec 2023 19:18:25 +0100 Subject: [PATCH 28/33] not needed --- src/phase.f90 | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index 148b028a7..2191ca791 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -384,7 +384,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief Initialize constitutive models for individual physics !-------------------------------------------------------------------------------------------------- -subroutine phase_init +subroutine phase_init() integer :: & ph, ce, co, ma @@ -544,11 +544,7 @@ subroutine crystallite_init() ip, & !< counter in integration point loop el, & !< counter in element loop en, ph - type(tDict), pointer :: & - num_phase, & - phases - phases => config_material%get_dict('phase') !$OMP PARALLEL DO PRIVATE(ce,ph,en) do el = 1, discretization_Nelems From 04fdac6556498e53c3d02ba3abb62624301fa859 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 28 Dec 2023 19:51:22 +0100 Subject: [PATCH 29/33] standard access pattern --- src/homogenization.f90 | 2 +- src/phase.f90 | 41 ++++++++++++++++++----------------------- 2 files changed, 19 insertions(+), 24 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 366e38053..fbc21cef2 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -312,7 +312,7 @@ subroutine homogenization_mechanical_response2(Delta_t,FEsolving_execIP,FEsolvin ce = (el-1)*discretization_nIPs + ip ho = material_ID_homogenization(ce) do co = 1, homogenization_Nconstituents(ho) - call crystallite_orientations(co,ip,el) + call crystallite_orientations(co,ce) end do call mechanical_homogenize(Delta_t,ce) end do IpLooping3 diff --git a/src/phase.f90 b/src/phase.f90 index 2191ca791..07f9bd317 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -541,21 +541,16 @@ subroutine crystallite_init() integer :: & ce, & co, & !< counter in integration point component loop - ip, & !< counter in integration point loop - el, & !< counter in element loop en, ph - !$OMP PARALLEL DO PRIVATE(ce,ph,en) - do el = 1, discretization_Nelems - do ip = 1, discretization_nIPs - ce = (el-1)*discretization_nIPs + ip - do co = 1,homogenization_Nconstituents(material_ID_homogenization(ce)) - en = material_entry_phase(co,ce) - ph = material_ID_phase(co,ce) - call crystallite_orientations(co,ip,el) - call plastic_dependentState(ph,en) ! update dependent state variables to be consistent with basic states - end do + !$OMP PARALLEL DO PRIVATE(ph,en) + do ce = 1, size(material_ID_homogenization) + do co = 1,homogenization_Nconstituents(material_ID_homogenization(ce)) + ph = material_ID_phase(co,ce) + en = material_entry_phase(co,ce) + call crystallite_orientations(co,ce) + call plastic_dependentState(ph,en) ! update dependent state variables to be consistent with basic states end do end do !$OMP END PARALLEL DO @@ -565,32 +560,30 @@ end subroutine crystallite_init !-------------------------------------------------------------------------------------------------- -!> @brief calculates orientations +!> @brief Update orientations and, if needed, compatibility. !-------------------------------------------------------------------------------------------------- -subroutine crystallite_orientations(co,ip,el) +subroutine crystallite_orientations(co,ce) integer, intent(in) :: & - co, & !< counter in integration point component loop - ip, & !< counter in integration point loop - el !< counter in element loop + co, & + ce integer :: ph, en - ph = material_ID_phase(co,(el-1)*discretization_nIPs + ip) - en = material_entry_phase(co,(el-1)*discretization_nIPs + ip) + ph = material_ID_phase(co,ce) + en = material_entry_phase(co,ce) call phase_O(ph)%data(en)%fromMatrix(transpose(math_rotationalPart(mechanical_F_e(ph,en)))) - if (plasticState(material_ID_phase(1,(el-1)*discretization_nIPs + ip))%nonlocal) & - call plastic_nonlocal_updateCompatibility(phase_O,(el-1)*discretization_nIPs + ip) + if (plasticState(material_ID_phase(1,ce))%nonlocal) call plastic_nonlocal_updateCompatibility(phase_O,ce) end subroutine crystallite_orientations !-------------------------------------------------------------------------------------------------- -!> @brief Map 2nd order tensor to reference config +!> @brief Map 2nd order tensor to reference configuration. !-------------------------------------------------------------------------------------------------- function crystallite_push33ToRef(co,ce, tensor33) @@ -618,11 +611,13 @@ end function crystallite_push33ToRef !-------------------------------------------------------------------------------------------------- logical pure function converged(residuum,state,atol) - real(pREAL), intent(in), dimension(:) ::& + real(pREAL), intent(in), dimension(:) :: & residuum, state, atol + real(pREAL) :: & rTol + rTol = num%rTol_crystalliteState converged = all(abs(residuum) <= max(atol, rtol*abs(state))) From b2b3a3f3f0e5840ff02bc5fb9f10fd1b8d5bcbc6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 28 Dec 2023 21:28:57 +0100 Subject: [PATCH 30/33] loop over cells, not (element, IP) --- src/Marc/materialpoint_Marc.f90 | 6 ++++-- src/grid/grid_mech_utilities.f90 | 2 +- src/homogenization.f90 | 23 ++++++++++------------- src/mesh/FEM_utilities.f90 | 2 +- src/phase.f90 | 2 +- 5 files changed, 17 insertions(+), 18 deletions(-) diff --git a/src/Marc/materialpoint_Marc.f90 b/src/Marc/materialpoint_Marc.f90 index c8db7d8df..3be9032f3 100644 --- a/src/Marc/materialpoint_Marc.f90 +++ b/src/Marc/materialpoint_Marc.f90 @@ -156,9 +156,11 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, materialpoint_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6) else validCalculation - call homogenization_mechanical_response(dt,(elCP-1)*discretization_nIPs + ip,(elCP-1)*discretization_nIPs + ip) + call homogenization_mechanical_response(dt,(elCP-1)*discretization_nIPs + ip, & + (elCP-1)*discretization_nIPs + ip) if (.not. terminallyIll) & - call homogenization_mechanical_response2(dt,[ip,ip],[elCP,elCP]) + call homogenization_mechanical_response2(dt,(elCP-1)*discretization_nIPs + ip, & + (elCP-1)*discretization_nIPs + ip) terminalIllness: if (terminallyIll) then diff --git a/src/grid/grid_mech_utilities.f90 b/src/grid/grid_mech_utilities.f90 index 0ea269319..dae4cd4ee 100644 --- a/src/grid/grid_mech_utilities.f90 +++ b/src/grid/grid_mech_utilities.f90 @@ -139,7 +139,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& if (.not. terminallyIll) & call homogenization_thermal_response(Delta_t,1,product(cells(1:2))*cells3) if (.not. terminallyIll) & - call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(cells(1:2))*cells3]) + call homogenization_mechanical_response2(Delta_t,1,product(cells(1:2))*cells3) P = reshape(homogenization_P, [3,3,cells(1),cells(2),cells3]) P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt diff --git a/src/homogenization.f90 b/src/homogenization.f90 index fbc21cef2..726c2fc2e 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -273,6 +273,7 @@ subroutine homogenization_thermal_response(Delta_t,cell_start,cell_end) real(pREAL), intent(in) :: Delta_t !< time increment integer, intent(in) :: & cell_start, cell_end + integer :: & co, ce, ho @@ -296,37 +297,33 @@ end subroutine homogenization_thermal_response !-------------------------------------------------------------------------------------------------- !> @brief !-------------------------------------------------------------------------------------------------- -subroutine homogenization_mechanical_response2(Delta_t,FEsolving_execIP,FEsolving_execElem) +subroutine homogenization_mechanical_response2(Delta_t,cell_start,cell_end) real(pREAL), intent(in) :: Delta_t !< time increment - integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP + integer, intent(in) :: & + cell_start, cell_end + integer :: & - ip, & !< integration point number - el, & !< element number co, ce, ho - !$OMP PARALLEL DO PRIVATE(ho,ce) - elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2) - IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2) - ce = (el-1)*discretization_nIPs + ip + !$OMP PARALLEL DO PRIVATE(ho) + do ce = cell_start, cell_end ho = material_ID_homogenization(ce) do co = 1, homogenization_Nconstituents(ho) call crystallite_orientations(co,ce) end do call mechanical_homogenize(Delta_t,ce) - end do IpLooping3 - end do elementLooping3 + end do !$OMP END PARALLEL DO - end subroutine homogenization_mechanical_response2 !-------------------------------------------------------------------------------------------------- !> @brief writes homogenization results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine homogenization_result +subroutine homogenization_result() integer :: ho character(len=:), allocatable :: group_base,group @@ -361,7 +358,7 @@ end subroutine homogenization_result !> @brief Forward data after successful increment. ! ToDo: Any guessing for the current states possible? !-------------------------------------------------------------------------------------------------- -subroutine homogenization_forward +subroutine homogenization_forward() integer :: ho diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 15a2168a8..0764f3443 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -148,7 +148,7 @@ subroutine utilities_constitutiveResponse(Delta_t,P_av,forwardData) call homogenization_mechanical_response(Delta_t,1,mesh_maxNips*mesh_NcpElems) ! calculate P field if (.not. terminallyIll) & - call homogenization_mechanical_response2(Delta_t,[1,mesh_maxNips],[1,mesh_NcpElems]) + call homogenization_mechanical_response2(Delta_t,1,mesh_maxNips*mesh_NcpElems) cutBack = .false. P_av = sum(homogenization_P,dim=3) * wgt diff --git a/src/phase.f90 b/src/phase.f90 index 07f9bd317..7cebda5b7 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -607,7 +607,7 @@ end function crystallite_push33ToRef !-------------------------------------------------------------------------------------------------- -!> @brief determines whether a point is converged +!> @brief Determine whether a point is converged. !-------------------------------------------------------------------------------------------------- logical pure function converged(residuum,state,atol) From 9ebb2d4ddd5fcfb59f8dc59d8757cee9d75491b3 Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 29 Dec 2023 00:04:44 +0100 Subject: [PATCH 31/33] [skip ci] updated version information after successful test of v3.0.0-alpha8-207-ga80fce410 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index bb39723c5..fdb577ca7 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -3.0.0-alpha8-199-g0f69eff16 +3.0.0-alpha8-207-ga80fce410 From 936328ae8bfb7a52e5d540a1fed99115a8131827 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 29 Dec 2023 08:32:52 +0100 Subject: [PATCH 32/33] re-enabled support for older CMake --- CMakeLists.txt | 2 +- cmake/Compiler-GNU.cmake | 2 ++ cmake/Compiler-Intel.cmake | 2 ++ cmake/Compiler-IntelLLVM.cmake | 2 ++ cmake/Compiler-LLVMFlang.cmake | 2 +- 5 files changed, 8 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index a8e551373..787117524 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -29,7 +29,7 @@ else() endif() add_definitions("-D${DAMASK_SOLVER}") -set(CMAKE_Fortran_PREPROCESS "ON") +set(CMAKE_Fortran_PREPROCESS "ON") # works only for CMake >= 3.18 # EXPERIMENTAL: This might help to detect HDF5 and FFTW3 in the future if PETSc is not aware of them set(ENV{PKG_CONFIG_PATH} "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/externalpackages:$ENV{PKG_CONFIG_PATH}") diff --git a/cmake/Compiler-GNU.cmake b/cmake/Compiler-GNU.cmake index 61ec2143b..43ed64af2 100644 --- a/cmake/Compiler-GNU.cmake +++ b/cmake/Compiler-GNU.cmake @@ -19,6 +19,8 @@ set (STANDARD_CHECK "-std=f2018 -pedantic-errors" ) #------------------------------------------------------------------------------------------------ # Fine tuning compilation options +set (COMPILE_FLAGS "${COMPILE_FLAGS} -cpp") # preprocessor, needed for CMake < 3.18 + set (COMPILE_FLAGS "${COMPILE_FLAGS} -fPIE") # position independent code diff --git a/cmake/Compiler-Intel.cmake b/cmake/Compiler-Intel.cmake index 59ba6d74d..d5e2dff7d 100644 --- a/cmake/Compiler-Intel.cmake +++ b/cmake/Compiler-Intel.cmake @@ -22,6 +22,8 @@ set (LINKER_FLAGS "${LINKER_FLAGS} -shared-intel") #------------------------------------------------------------------------------------------------ # Fine tuning compilation options +set (COMPILE_FLAGS "${COMPILE_FLAGS} -fpp") # preprocessor, needed for CMake < 3.18 + set (COMPILE_FLAGS "${COMPILE_FLAGS} -no-ftz") # disable flush underflow to zero, will be set if -O[1,2,3] diff --git a/cmake/Compiler-IntelLLVM.cmake b/cmake/Compiler-IntelLLVM.cmake index 3749b925f..079574c4b 100644 --- a/cmake/Compiler-IntelLLVM.cmake +++ b/cmake/Compiler-IntelLLVM.cmake @@ -24,6 +24,8 @@ set (LINKER_FLAGS "${LINKER_FLAGS} -shared-intel -fc=ifx") #------------------------------------------------------------------------------------------------ # Fine tuning compilation options +set (COMPILE_FLAGS "${COMPILE_FLAGS} -fpp") # preprocessor, needed for CMake < 3.18 + set (COMPILE_FLAGS "${COMPILE_FLAGS} -no-ftz") # disable flush underflow to zero, will be set if -O[1,2,3] diff --git a/cmake/Compiler-LLVMFlang.cmake b/cmake/Compiler-LLVMFlang.cmake index b28df7027..f6be61b45 100644 --- a/cmake/Compiler-LLVMFlang.cmake +++ b/cmake/Compiler-LLVMFlang.cmake @@ -9,4 +9,4 @@ set (STANDARD_CHECK "-std=f2018 -pedantic" ) #------------------------------------------------------------------------------------------------ # Fine tuning compilation options - +set (COMPILE_FLAGS "${COMPILE_FLAGS} -cpp") # preprocessor, needed for CMake < 3.18 From 8797af49ae7f1fdff8235265af901b800bfd10f8 Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 29 Dec 2023 11:19:04 +0100 Subject: [PATCH 33/33] [skip ci] updated version information after successful test of v3.0.0-alpha8-211-gccf4867e0 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index fdb577ca7..360fe5135 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -3.0.0-alpha8-207-ga80fce410 +3.0.0-alpha8-211-gccf4867e0