diff --git a/CMakeLists.txt b/CMakeLists.txt index a501717ed..787117524 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -29,6 +29,8 @@ else() endif() add_definitions("-D${DAMASK_SOLVER}") +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}") pkg_check_modules(HDF5 hdf5) @@ -91,24 +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") 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/PRIVATE b/PRIVATE index 117b65d85..62df7f24f 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 117b65d852659158c0f4ca3bf8ce8db51a7a1961 +Subproject commit 62df7f24f2a95fda255f7d20b130afcfeecb1b4a diff --git a/VERSION b/VERSION index 4f15e26cc..360fe5135 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -3.0.0-alpha8-169-g85191cd02 +3.0.0-alpha8-211-gccf4867e0 diff --git a/cmake/Compiler-GNU.cmake b/cmake/Compiler-GNU.cmake index a62875d05..43ed64af2 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 () @@ -23,8 +19,7 @@ set (STANDARD_CHECK "-std=f2018 -pedantic-errors" ) #------------------------------------------------------------------------------------------------ # Fine tuning compilation options -set (COMPILE_FLAGS "${COMPILE_FLAGS} -cpp") -# preprocessor +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 19e75a9fa..d5e2dff7d 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 () @@ -26,8 +22,7 @@ set (LINKER_FLAGS "${LINKER_FLAGS} -shared-intel") #------------------------------------------------------------------------------------------------ # Fine tuning compilation options -set (COMPILE_FLAGS "${COMPILE_FLAGS} -fpp") -# preprocessor +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 76fc0386a..079574c4b 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 () @@ -28,8 +24,7 @@ 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} -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 new file mode 100644 index 000000000..f6be61b45 --- /dev/null +++ b/cmake/Compiler-LLVMFlang.cmake @@ -0,0 +1,12 @@ +################################################################################################### +# LLVM Compiler +################################################################################################### +if (OPENMP) + set (OPENMP_FLAGS "-fopenmp") +endif () + +set (STANDARD_CHECK "-std=f2018 -pedantic" ) + +#------------------------------------------------------------------------------------------------ +# Fine tuning compilation options +set (COMPILE_FLAGS "${COMPILE_FLAGS} -cpp") # preprocessor, needed for CMake < 3.18 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 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/_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/_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 eda668564..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, @@ -618,7 +627,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..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,12 +472,16 @@ 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: 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): diff --git a/src/IO.f90 b/src/IO.f90 index a08c624c3..89f86ec63 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 @@ -380,17 +381,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 +397,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 +421,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 +488,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:' @@ -742,6 +728,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 +821,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 +901,54 @@ 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 only separators' + + 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)//'"' ! requires 2018 standard + if (solution(i) /= tok(i)) error stop 'tokenize' + end do + + end subroutine test_tokenize + end subroutine IO_selfTest end module IO diff --git a/src/Marc/materialpoint_Marc.f90 b/src/Marc/materialpoint_Marc.f90 index 0e3835924..3be9032f3 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,10 +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_F0(1:3,1:3,ce) = ffn homogenization_F(1:3,1:3,ce) = ffn1 + materialpoint_F(1:3,1:3,ip,elCP) = ffn1 if (iand(mode, materialpoint_CALCRESULTS) /= 0) then @@ -154,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 @@ -168,17 +172,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 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/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 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 754942ae1..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 @@ -269,7 +270,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 +390,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..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 @@ -226,7 +227,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 +347,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..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 @@ -255,7 +256,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 +391,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/grid_mech_utilities.f90 b/src/grid/grid_mech_utilities.f90 new file mode 100644 index 000000000..dae4cd4ee --- /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,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 9ed5cbe17..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 from homogenization_F0 to F during Delta_t -!-------------------------------------------------------------------------------------------------- -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 diff --git a/src/homogenization.f90 b/src/homogenization.f90 index ebb18b232..726c2fc2e 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 @@ -274,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 @@ -297,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,ip,el) + 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 @@ -362,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/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/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index cc906a92b..f0bafbac3 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -131,7 +131,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/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 16cf94ae5..f040f6b50 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) diff --git a/src/phase.f90 b/src/phase.f90 index ffa606cc2..7cebda5b7 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 @@ -387,7 +384,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief Initialize constitutive models for individual physics !-------------------------------------------------------------------------------------------------- -subroutine phase_init +subroutine phase_init() integer :: & ph, ce, co, ma @@ -544,25 +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 - type(tDict), pointer :: & - num_phase, & - phases - phases => config_material%get_dict('phase') - !$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 @@ -572,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,material_ID_phase(1,(el-1)*discretization_nIPs + ip),ip,el) + 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) @@ -621,15 +607,17 @@ 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) - 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))) diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 545dec4e6..bd855e46f 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 :: & @@ -124,7 +123,9 @@ submodule(phase:plastic) nonlocal type :: tNonlocalDependentState real(pREAL), allocatable, dimension(:,:) :: & tau_pass, & - tau_back + tau_back, & + rho_forest, & + max_dipole_height 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,11 +389,9 @@ 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 + '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 @@ -477,15 +475,17 @@ 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) + 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) + 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%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 @@ -518,20 +517,14 @@ 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 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 @@ -602,7 +595,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 +605,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 @@ -861,14 +854,14 @@ 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) - forall (s = 1:prm%sum_N_sl, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ph),en) + v = reshape(stt%v(:,en),[prm%sum_N_sl,4]) + 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]) @@ -975,7 +968,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) @@ -990,11 +984,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) @@ -1018,20 +1011,20 @@ 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) + v0 = reshape(st0%v(:,en),[prm%sum_N_sl,4]) !**************************************************************************** @@ -1074,7 +1067,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 @@ -1171,7 +1164,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 @@ -1181,10 +1175,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) @@ -1331,18 +1325,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 +1345,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 @@ -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)), &