diff --git a/CMakeLists.txt b/CMakeLists.txt index a501717ed..0f2cf2763 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -17,17 +17,16 @@ pkg_get_variable(CMAKE_Fortran_COMPILER PETSc fcompiler) pkg_get_variable(CMAKE_C_COMPILER PETSc ccompiler) # Solver determines name of project -string(TOUPPER "${DAMASK_SOLVER}" DAMASK_SOLVER) -if(DAMASK_SOLVER STREQUAL "GRID") - project(damask-grid HOMEPAGE_URL https://damask.mpie.de LANGUAGES Fortran C) -elseif(DAMASK_SOLVER STREQUAL "MESH") - project(damask-mesh HOMEPAGE_URL https://damask.mpie.de LANGUAGES Fortran C) -elseif(DAMASK_SOLVER STREQUAL "TEST") - project(damask-test HOMEPAGE_URL https://damask.mpie.de LANGUAGES Fortran C) +string(TOUPPER "${DAMASK_SOLVER}" DAMASK_SOLVER_UPPER) +string(TOLOWER "${DAMASK_SOLVER}" DAMASK_SOLVER_LOWER) +if("${DAMASK_SOLVER_UPPER}" MATCHES "^(GRID|MESH|TEST)$") + project("damask-${DAMASK_SOLVER_LOWER}" HOMEPAGE_URL https://damask.mpie.de LANGUAGES Fortran C) else() message(FATAL_ERROR "Invalid solver: DAMASK_SOLVER=${DAMASK_SOLVER}") endif() -add_definitions("-D${DAMASK_SOLVER}") +add_definitions("-D${DAMASK_SOLVER_UPPER}") + +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}") @@ -91,24 +90,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 5a715996e..62df7f24f 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 5a715996e6e8418a59fbcaf3715a2516ad05ed51 +Subproject commit 62df7f24f2a95fda255f7d20b130afcfeecb1b4a diff --git a/VERSION b/VERSION index a45b5eb13..cdbc225bb 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -3.0.0-alpha8-172-gec624a86a +3.0.0-alpha8-246-g112b5be1a 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/_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 74c731ae0..b4ee64283 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 time_stamp() -> str: """Provide current time as formatted string.""" diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index e9e054171..b03b0249a 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) @@ -489,7 +489,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: diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 4d0d8cef0..1aef09c3f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,16 +1,7 @@ file(GLOB damask-sources CONFIGURE_DEPENDS *.f90 *.c) -if(PROJECT_NAME STREQUAL "damask-grid") - set(executable-name "DAMASK_grid") - file(GLOB solver-sources CONFIGURE_DEPENDS grid/*.f90) -elseif(PROJECT_NAME STREQUAL "damask-mesh") - set(executable-name "DAMASK_mesh") - file(GLOB solver-sources CONFIGURE_DEPENDS mesh/*.f90) -elseif(PROJECT_NAME STREQUAL "damask-test") - set(executable-name "DAMASK_test") - file(GLOB solver-sources CONFIGURE_DEPENDS test/*.f90) -endif() - +set(executable-name "DAMASK_${DAMASK_SOLVER_LOWER}") +file(GLOB solver-sources CONFIGURE_DEPENDS ${DAMASK_SOLVER_LOWER}/*.f90) if(NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY") add_executable(${executable-name} ${damask-sources} ${solver-sources}) diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90 index 65ee66af9..413446399 100644 --- a/src/HDF5_utilities.f90 +++ b/src/HDF5_utilities.f90 @@ -1560,7 +1560,7 @@ subroutine initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_ call HDF5_chkerr(hdferr) call MPI_Allgather(int(localShape(ubound(localShape,1)),MPI_INTEGER_KIND),1_MPI_INTEGER_KIND,MPI_INTEGER,& readSize,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) end if #endif myStart = int(0,HSIZE_T) @@ -1667,7 +1667,7 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, & if (parallel) then call MPI_Allgather(int(localShape(ubound(localShape,1)),MPI_INTEGER_KIND),1_MPI_INTEGER_KIND,MPI_INTEGER,& writeSize,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) end if #endif myStart = int(0,HSIZE_T) 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..4d8f333ec 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,8 @@ 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) - if (.not. terminallyIll) & - call homogenization_mechanical_response2(dt,[ip,ip],[elCP,elCP]) + call homogenization_mechanical_response(dt,(elCP-1)*discretization_nIPs + ip, & + (elCP-1)*discretization_nIPs + ip) terminalIllness: if (terminallyIll) then @@ -168,17 +169,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..68ed40cf3 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 @@ -362,7 +363,7 @@ program DAMASK_grid end if; flush(IO_STDOUT) call MPI_Allreduce(signal_SIGUSR1,sig,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) if (mod(inc,loadCases(l)%f_out) == 0 .or. sig) then print'(/,1x,a)', '... saving results ........................................................' flush(IO_STDOUT) @@ -370,7 +371,7 @@ program DAMASK_grid end if if (sig) call signal_setSIGUSR1(.false.) call MPI_Allreduce(signal_SIGUSR2,sig,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) if (mod(inc,loadCases(l)%f_restart) == 0 .or. sig) then do field = 1, nActiveFields select case (ID(field)) @@ -386,7 +387,7 @@ program DAMASK_grid end if if (sig) call signal_setSIGUSR2(.false.) call MPI_Allreduce(signal_SIGINT,sig,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) if (sig) exit loadCaseLooping end if skipping 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/discretization_grid.f90 b/src/grid/discretization_grid.f90 index 2cb5dbf9f..0c18b3317 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -97,12 +97,12 @@ subroutine discretization_grid_init(restart) call MPI_Bcast(cells,3_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) if (cells(1) < 2) call IO_error(844, ext_msg='cells(1) must be larger than 1') call MPI_Bcast(geomSize,3_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call MPI_Bcast(origin,3_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) print'(/,1x,a,i0,a,i0,a,i0)', 'cells: ', cells(1), ' × ', cells(2), ' × ', cells(3) print '(1x,a,es8.2,a,es8.2,a,es8.2,a)', 'size: ', geomSize(1), ' × ', geomSize(2), ' × ', geomSize(3), ' m³' @@ -126,15 +126,15 @@ subroutine discretization_grid_init(restart) call MPI_Gather(product(cells(1:2))*cells3Offset,1_MPI_INTEGER_KIND,MPI_INTEGER,displs,& 1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call MPI_Gather(product(myGrid), 1_MPI_INTEGER_KIND,MPI_INTEGER,sendcounts,& 1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) allocate(materialAt(product(myGrid))) call MPI_Scatterv(materialAt_global,sendcounts,displs,MPI_INTEGER,materialAt,size(materialAt),& MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call discretization_init(materialAt, & IPcoordinates0(myGrid,mySize,cells3Offset), & @@ -318,10 +318,10 @@ function discretization_grid_getInitialCondition(label) result(ic) real(pREAL), dimension(:), allocatable :: ic_global, ic_local integer(MPI_INTEGER_KIND) :: err_MPI - integer, dimension(worldsize) :: & displs, sendcounts + if (worldrank == 0) then ic_global = VTI_readDataset_real(IO_read(CLI_geomFile),label) else @@ -330,15 +330,15 @@ function discretization_grid_getInitialCondition(label) result(ic) call MPI_Gather(product(cells(1:2))*cells3Offset, 1_MPI_INTEGER_KIND,MPI_INTEGER,displs,& 1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call MPI_Gather(product(cells(1:2))*cells3, 1_MPI_INTEGER_KIND,MPI_INTEGER,sendcounts,& 1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) allocate(ic_local(product(cells(1:2))*cells3)) call MPI_Scatterv(ic_global,sendcounts,displs,MPI_DOUBLE,ic_local,size(ic_local),& MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) ic = reshape(ic_local,[cells(1),cells(2),cells3]) diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 648f9ebbb..6b53c0a75 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, & @@ -130,7 +129,7 @@ subroutine grid_damage_spectral_init(num_grid) CHKERRQ(err_PETSc) call MPI_Allgather(int(cells3,MPI_INTEGER_KIND),1_MPI_INTEGER_KIND,MPI_INTEGER,& cells3_global,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call DMDACreate3D(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point @@ -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) @@ -241,10 +239,10 @@ function grid_damage_spectral_solution(Delta_t) result(solution) phi_max = maxval(phi) stagNorm = maxval(abs(phi - phi_stagInc)) call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) solution%stagConverged = stagNorm < max(num%eps_damage_atol, num%eps_damage_rtol*phi_max) call MPI_Allreduce(MPI_IN_PLACE,solution%stagConverged,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LAND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) phi_stagInc = phi call homogenization_set_phi(reshape(phi,[product(cells(1:2))*cells3])) @@ -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 @@ -381,10 +379,10 @@ subroutine updateReference() K_ref = K_ref*wgt call MPI_Allreduce(MPI_IN_PLACE,K_ref,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) mu_ref = mu_ref*wgt call MPI_Allreduce(MPI_IN_PLACE,mu_ref,1_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) end subroutine updateReference diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 754942ae1..2b856ff1c 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 @@ -172,7 +173,7 @@ subroutine grid_mechanical_FEM_init(num_grid) CHKERRQ(err_PETSc) call MPI_Allgather(int(cells3,MPI_INTEGER_KIND),1_MPI_INTEGER_KIND,MPI_INTEGER,& cells3_global,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call DMDACreate3d(PETSC_COMM_WORLD, & DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, & DMDA_STENCIL_BOX, & @@ -245,16 +246,16 @@ subroutine grid_mechanical_FEM_init(num_grid) call HDF5_read(P_aim,groupHandle,'P_aim',.false.) call MPI_Bcast(P_aim,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call HDF5_read(F_aim,groupHandle,'F_aim',.false.) call MPI_Bcast(F_aim,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call HDF5_read(F_aim_lastInc,groupHandle,'F_aim_lastInc',.false.) call MPI_Bcast(F_aim_lastInc,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call HDF5_read(F_aimDot,groupHandle,'F_aimDot',.false.) call MPI_Bcast(F_aimDot,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call HDF5_read(temp33n,groupHandle,'F') F = reshape(temp33n,[3,3,cells(1),cells(2),cells3]) call HDF5_read(temp33n,groupHandle,'F_lastInc') @@ -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 @@ -283,10 +283,10 @@ subroutine grid_mechanical_FEM_init(num_grid) print'(1x,a,1x,i0)', 'loading additional restart data of increment', CLI_restartInc call HDF5_read(C_volAvg,groupHandle,'C_volAvg',.false.) call MPI_Bcast(C_volAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call HDF5_read(C_volAvgLastInc,groupHandle,'C_volAvgLastInc',.false.) call MPI_Bcast(C_volAvgLastInc,81_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call HDF5_closeGroup(groupHandle) call HDF5_closeFile(fileHandle) @@ -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 !-------------------------------------------------------------------------------------------------- @@ -576,7 +575,7 @@ subroutine formResidual(da_local,x_local, & P_av,C_volAvg,devNull, & F,params%Delta_t,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) !-------------------------------------------------------------------------------------------------- ! stress BC handling diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 5bcfec438..22113a6fd 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 @@ -168,7 +169,7 @@ subroutine grid_mechanical_spectral_basic_init(num_grid) CHKERRQ(err_PETSc) call MPI_Allgather(int(cells3,MPI_INTEGER_KIND),1_MPI_INTEGER_KIND,MPI_INTEGER,& cells3_global,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call DMDACreate3d(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point @@ -206,16 +207,16 @@ subroutine grid_mechanical_spectral_basic_init(num_grid) call HDF5_read(P_aim,groupHandle,'P_aim',.false.) call MPI_Bcast(P_aim,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call HDF5_read(F_aim,groupHandle,'F_aim',.false.) call MPI_Bcast(F_aim,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call HDF5_read(F_aim_lastInc,groupHandle,'F_aim_lastInc',.false.) call MPI_Bcast(F_aim_lastInc,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call HDF5_read(F_aimDot,groupHandle,'F_aimDot',.false.) call MPI_Bcast(F_aimDot,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call HDF5_read(temp33n,groupHandle,'F') F = reshape(temp33n,[9,cells(1),cells(2),cells3]) call HDF5_read(temp33n,groupHandle,'F_lastInc') @@ -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 @@ -238,13 +238,13 @@ subroutine grid_mechanical_spectral_basic_init(num_grid) print'(1x,a,1x,i0)', 'loading additional restart data of increment', CLI_restartInc call HDF5_read(C_volAvg,groupHandle,'C_volAvg',.false.) call MPI_Bcast(C_volAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call HDF5_read(C_volAvgLastInc,groupHandle,'C_volAvgLastInc',.false.) call MPI_Bcast(C_volAvgLastInc,81_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call HDF5_read(C_minMaxAvg,groupHandle,'C_minMaxAvg',.false.) call MPI_Bcast(C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call HDF5_closeGroup(groupHandle) call HDF5_closeFile(fileHandle) @@ -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 !-------------------------------------------------------------------------------------------------- @@ -521,7 +519,7 @@ subroutine formResidual(residual_subdomain, F, & P_av,C_volAvg,C_minMaxAvg, & F,params%Delta_t,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) err_div = utilities_divergenceRMS(P) end associate diff --git a/src/grid/grid_mech_spectral_polarization.f90 b/src/grid/grid_mech_spectral_polarization.f90 index efa1a64a9..a4da7452d 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 @@ -189,7 +190,7 @@ subroutine grid_mechanical_spectral_polarization_init(num_grid) CHKERRQ(err_PETSc) call MPI_Allgather(int(cells3,pPetscInt),1_MPI_INTEGER_KIND,MPI_INTEGER,& cells3_global,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call DMDACreate3d(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point @@ -229,16 +230,16 @@ subroutine grid_mechanical_spectral_polarization_init(num_grid) call HDF5_read(P_aim,groupHandle,'P_aim',.false.) call MPI_Bcast(P_aim,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call HDF5_read(F_aim,groupHandle,'F_aim',.false.) call MPI_Bcast(F_aim,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call HDF5_read(F_aim_lastInc,groupHandle,'F_aim_lastInc',.false.) call MPI_Bcast(F_aim_lastInc,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call HDF5_read(F_aimDot,groupHandle,'F_aimDot',.false.) call MPI_Bcast(F_aimDot,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call HDF5_read(temp33n,groupHandle,'F') F = reshape(temp33n,[9,cells(1),cells(2),cells3]) call HDF5_read(temp33n,groupHandle,'F_lastInc') @@ -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 @@ -267,13 +267,13 @@ subroutine grid_mechanical_spectral_polarization_init(num_grid) print '(1x,a,1x,i0)', 'loading additional restart data of increment', CLI_restartInc call HDF5_read(C_volAvg,groupHandle,'C_volAvg',.false.) call MPI_Bcast(C_volAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call HDF5_read(C_volAvgLastInc,groupHandle,'C_volAvgLastInc',.false.) call MPI_Bcast(C_volAvgLastInc,81_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call HDF5_read(C_minMaxAvg,groupHandle,'C_minMaxAvg',.false.) call MPI_Bcast(C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call HDF5_closeGroup(groupHandle) call HDF5_closeFile(fileHandle) @@ -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 !-------------------------------------------------------------------------------------------------- @@ -574,7 +573,7 @@ subroutine formResidual(residual_subdomain, FandF_tau, & F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt call MPI_Allreduce(MPI_IN_PLACE,F_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' + call parallelization_chkerr(err_MPI) call SNESGetNumberFunctionEvals(SNES_mech,nfuncs,err_PETSc) CHKERRQ(err_PETSc) diff --git a/src/grid/grid_mech_utilities.f90 b/src/grid/grid_mech_utilities.f90 new file mode 100644 index 000000000..bae5c309b --- /dev/null +++ b/src/grid/grid_mech_utilities.f90 @@ -0,0 +1,249 @@ +!-------------------------------------------------------------------------------------------------- +!> @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 + + 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) + call parallelization_chkerr(err_MPI) + 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) + call parallelization_chkerr(err_MPI) + call MPI_Bcast(dPdF_max,81_MPI_INTEGER_KIND,MPI_DOUBLE,int(valueAndRank(2),MPI_INTEGER_KIND),MPI_COMM_WORLD,err_MPI) + call parallelization_chkerr(err_MPI) + + 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) + call parallelization_chkerr(err_MPI) + call MPI_Bcast(dPdF_min,81_MPI_INTEGER_KIND,MPI_DOUBLE,int(valueAndRank(2),MPI_INTEGER_KIND),MPI_COMM_WORLD,err_MPI) + call parallelization_chkerr(err_MPI) + + 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) + call parallelization_chkerr(err_MPI) + 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) + call parallelization_chkerr(err_MPI) + 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..1c46050ff 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, & @@ -124,7 +123,7 @@ subroutine grid_thermal_spectral_init(num_grid) CHKERRQ(err_PETSc) call MPI_Allgather(int(cells3,pPETSCINT),1_MPI_INTEGER_KIND,MPI_INTEGER,& cells3_global,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call DMDACreate3D(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point @@ -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) @@ -220,14 +218,14 @@ function grid_thermal_spectral_solution(Delta_t) result(solution) T_max = maxval(T) stagNorm = maxval(abs(T - T_stagInc)) call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) solution%stagConverged = stagNorm < max(num%eps_thermal_atol, num%eps_thermal_rtol*T_max) call MPI_Allreduce(MPI_IN_PLACE,solution%stagConverged,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LAND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) 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 @@ -325,6 +323,8 @@ subroutine formResidual(residual_subdomain,x_scal,r,dummy,err_PETSc) real(pREAL), dimension(3,cells(1),cells(2),cells3) :: vectorField + call homogenization_thermal_response(Delta_t_,1,product(cells(1:2))*cells3) + associate(T => x_scal) vectorField = utilities_ScalarGradient(T) ce = 0 @@ -336,13 +336,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 @@ -367,10 +367,10 @@ subroutine updateReference() K_ref = K_ref*wgt call MPI_Allreduce(MPI_IN_PLACE,K_ref,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) mu_ref = mu_ref*wgt call MPI_Allreduce(MPI_IN_PLACE,mu_ref,1_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) end subroutine updateReference diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 9ed5cbe17..90544b2f8 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -75,22 +75,9 @@ 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] + divergence_correction !< scale divergence/curl calculation logical :: & memory_efficient !< calculate gamma operator on the fly end type tNumerics @@ -121,10 +108,6 @@ module spectral_utilities utilities_curlRMS, & utilities_scalarGradient, & utilities_vectorDivergence, & - utilities_maskedCompliance, & - utilities_constitutiveResponse, & - utilities_calculateRate, & - utilities_forwardTensorField, & utilities_updateCoords contains @@ -580,7 +563,7 @@ real(pREAL) function utilities_divergenceRMS(tensorField) conjg(-xi1st(1:3,cells1Red,k,j))*rescaledGeom))**2) end do; end do call MPI_Allreduce(MPI_IN_PLACE,utilities_divergenceRMS,1_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) utilities_divergenceRMS = sqrt(utilities_divergenceRMS) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space if (cells(1) == 1) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pREAL ! counted twice in case of cells(1) == 1 @@ -646,72 +629,13 @@ real(pREAL) function utilities_curlRMS(tensorField) end do; end do call MPI_Allreduce(MPI_IN_PLACE,utilities_curlRMS,1_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) utilities_curlRMS = sqrt(utilities_curlRMS) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space if (cells(1) == 1) utilities_curlRMS = utilities_curlRMS * 0.5_pREAL ! counted twice in case of cells(1) == 1 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 @@ -995,7 +778,7 @@ subroutine utilities_updateCoords(F) ! average F if (cells3Offset == 0) Favg = tensorField_fourier(1:3,1:3,1,1,1)%re*wgt call MPI_Bcast(Favg,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) !-------------------------------------------------------------------------------------------------- ! integration in Fourier space to get fluctuations of cell center displacements @@ -1015,24 +798,24 @@ subroutine utilities_updateCoords(F) !-------------------------------------------------------------------------------------------------- ! pad cell center fluctuations along z-direction (needed when running MPI simulation) - c = product(shape(u_tilde_p_padded(:,:,:,1))) !< amount of data to transfer + c = product(shape(u_tilde_p_padded(:,:,:,1))) !< amount of data to transfer rank_t = modulo(worldrank+1_MPI_INTEGER_KIND,worldsize) rank_b = modulo(worldrank-1_MPI_INTEGER_KIND,worldsize) ! send bottom layer to process below call MPI_Isend(u_tilde_p_padded(:,:,:,1), c,MPI_DOUBLE,rank_b,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(1),err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call MPI_Irecv(u_tilde_p_padded(:,:,:,cells3+1),c,MPI_DOUBLE,rank_t,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(2),err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) ! send top layer to process above call MPI_Isend(u_tilde_p_padded(:,:,:,cells3) ,c,MPI_DOUBLE,rank_t,1_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(3),err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call MPI_Irecv(u_tilde_p_padded(:,:,:,0), c,MPI_DOUBLE,rank_b,1_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(4),err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call MPI_Waitall(4,request,status,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) ! ToDo #else @@ -1085,7 +868,7 @@ subroutine selfTest() call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier) call MPI_Allreduce(sum(sum(sum(tensorField_real_,dim=5),dim=4),dim=3),tensorSum,9_MPI_INTEGER_KIND, & MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) if (worldrank==0) then if (any(dNeq(tensorSum/tensorField_fourier(:,:,1,1,1)%re,1.0_pREAL,1.0e-12_pREAL))) & error stop 'mismatch avg tensorField FFT <-> real' @@ -1101,7 +884,7 @@ subroutine selfTest() call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier) call MPI_Allreduce(sum(sum(sum(vectorField_real_,dim=4),dim=3),dim=2),vectorSum,3_MPI_INTEGER_KIND, & MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) if (worldrank==0) then if (any(dNeq(vectorSum/vectorField_fourier(:,1,1,1)%re,1.0_pREAL,1.0e-12_pREAL))) & error stop 'mismatch avg vectorField FFT <-> real' @@ -1117,7 +900,7 @@ subroutine selfTest() call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier) call MPI_Allreduce(sum(sum(sum(scalarField_real_,dim=3),dim=2),dim=1),scalarSum,1_MPI_INTEGER_KIND, & MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) if (worldrank==0) then if (dNeq(scalarSum/scalarField_fourier(1,1,1)%re,1.0_pREAL,1.0e-12_pREAL)) & error stop 'mismatch avg scalarField FFT <-> real' @@ -1129,7 +912,7 @@ subroutine selfTest() call random_number(r) call MPI_Bcast(r,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) scalarField_real_ = r(1,1) if (maxval(abs(utilities_scalarGradient(scalarField_real_)))>5.0e-9_pREAL) error stop 'non-zero grad(const)' diff --git a/src/homogenization.f90 b/src/homogenization.f90 index ebb18b232..2da1d73b2 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 @@ -169,7 +168,6 @@ module homogenization public :: & homogenization_init, & homogenization_mechanical_response, & - homogenization_mechanical_response2, & homogenization_thermal_response, & homogenization_thermal_active, & homogenization_mu_T, & @@ -228,7 +226,8 @@ subroutine homogenization_mechanical_response(Delta_t,cell_start,cell_end) doneAndHappy - !$OMP PARALLEL DO PRIVATE(en,ho,co,converged,doneAndHappy) + !$OMP PARALLEL + !$OMP DO PRIVATE(en,ho,co,converged,doneAndHappy) do ce = cell_start, cell_end en = material_entry_homogenization(ce) @@ -261,7 +260,18 @@ subroutine homogenization_mechanical_response(Delta_t,cell_start,cell_end) terminallyIll = .true. end if end do - !$OMP END PARALLEL DO + !$OMP END DO + + !$OMP 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 + !$OMP END DO + !$OMP END PARALLEL end subroutine homogenization_mechanical_response @@ -274,6 +284,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 @@ -294,40 +305,10 @@ subroutine homogenization_thermal_response(Delta_t,cell_start,cell_end) end subroutine homogenization_thermal_response -!-------------------------------------------------------------------------------------------------- -!> @brief -!-------------------------------------------------------------------------------------------------- -subroutine homogenization_mechanical_response2(Delta_t,FEsolving_execIP,FEsolving_execElem) - - real(pREAL), intent(in) :: Delta_t !< time increment - integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP - 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 - ho = material_ID_homogenization(ce) - do co = 1, homogenization_Nconstituents(ho) - call crystallite_orientations(co,ip,el) - end do - call mechanical_homogenize(Delta_t,ce) - end do IpLooping3 - end do elementLooping3 - !$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 +343,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/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index 370dc3535..6b6400d62 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -23,43 +23,32 @@ program DAMASK_mesh implicit none(type,external) type :: tLoadCase - real(pREAL) :: time = 0.0_pREAL !< length of increment - integer :: incs = 0, & !< number of increments - outputfrequency = 1 !< frequency of result writes - logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase - integer, allocatable, dimension(:) :: faceID - type(tFieldBC), allocatable, dimension(:) :: fieldBC + real(pREAL) :: t = 0.0_pREAL !< length of increment + integer :: N = 0, & !< number of increments + f_out = 1 !< frequency of result writes + logical :: estimate_rate = .true. !< follow trajectory of former loadcase + type(tMechBC), allocatable, dimension(:) :: mechBC end type tLoadCase -!-------------------------------------------------------------------------------------------------- -! variables related to information from load case and geom file - integer, allocatable, dimension(:) :: chunkPos ! this is longer than needed for geometry parsing - integer :: & - N_def = 0 !< # of rate of deformation specifiers found in load case file - character(len=:), allocatable :: & - line !-------------------------------------------------------------------------------------------------- ! loop variables, convergence etc. integer, parameter :: & subStepFactor = 2 !< for each substep, divide the last time increment by 2.0 real(pREAL) :: & - time = 0.0_pREAL, & !< elapsed time - time0 = 0.0_pREAL, & !< begin of interval - timeinc = 0.0_pREAL, & !< current time interval - timeIncOld = 0.0_pREAL, & !< previous time interval - remainingLoadCaseTime = 0.0_pREAL !< remaining time of current load case + t = 0.0_pREAL, & !< elapsed time + t_0 = 0.0_pREAL, & !< begin of interval + Delta_t = 0.0_pREAL, & !< current time interval + Delta_t_prev = 0.0_pREAL !< previous time interval logical :: & guess, & !< guess along former trajectory stagIterate integer :: & l, & - i, & + m, & errorID, & cutBackLevel = 0, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$ stepFraction = 0, & !< fraction of current time interval - currentLoadcase = 0, & !< current load case - currentFace = 0, & inc, & !< current increment in current load case totalIncsCounter = 0, & !< total # of increments statUnit = 0, & !< file unit for statistics output @@ -67,8 +56,16 @@ program DAMASK_mesh component type(tDict), pointer :: & num_solver, & - num_mesh - character(len=pSTRLEN), dimension(:), allocatable :: fileContent + num_mesh, & + load, & + load_step, & + step_bc, & + mech_BC, & + step_discretization + type(tList), pointer :: & + load_steps, & + mech_u, & + step_mech character(len=pSTRLEN) :: & incInfo, & loadcase_string @@ -80,9 +77,11 @@ program DAMASK_mesh type(tSolutionState), allocatable, dimension(:) :: solres PetscInt :: faceSet, currentFaceSet, dimPlex PetscErrorCode :: err_PETSc - integer(kind(COMPONENT_UNDEFINED_ID)) :: ID external :: & quit + character(len=:), allocatable :: & + fileContent, fname + !-------------------------------------------------------------------------------------------------- ! init DAMASK (all modules) @@ -104,135 +103,92 @@ program DAMASK_mesh CHKERRA(err_PETSc) allocate(solres(1)) -!-------------------------------------------------------------------------------------------------- -! reading basic information from load case file and allocate data structure containing load cases - fileContent = IO_readlines(trim(CLI_loadFile)) - do l = 1, size(fileContent) - line = fileContent(l) - if (IO_isBlank(line)) cycle ! skip empty lines + if (worldrank == 0) then + fileContent = IO_read(CLI_loadFile) + fname = CLI_loadFile + if (scan(fname,'/') /= 0) fname = fname(scan(fname,'/',.true.)+1:) + call result_openJobFile(parallel=.false.) + call result_addSetupFile(fileContent,fname,'load case definition (mesh solver)') + call result_closeJobFile() + end if - chunkPos = IO_strPos(line) - do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase - select case (IO_strValue(line,chunkPos,i)) - case('$Loadcase') - N_def = N_def + 1 - end select - end do ! count all identifiers to allocate memory and do sanity check - end do + call parallelization_bcast_str(fileContent) + load => YAML_parse_str_asDict(fileContent) + load_steps => load%get_list('loadstep') - if (N_def < 1) call IO_error(error_ID = 837) - allocate(loadCases(N_def)) + allocate(loadCases(load_steps%length)) - do i = 1, size(loadCases) - allocate(loadCases(i)%fieldBC(1)) - loadCases(i)%fieldBC(1)%ID = FIELD_MECH_ID - end do - do i = 1, size(loadCases) - loadCases(i)%fieldBC(1)%nComponents = dimPlex !< X, Y (, Z) displacements - allocate(loadCases(i)%fieldBC(1)%componentBC(loadCases(i)%fieldBC(1)%nComponents)) - do component = 1, loadCases(i)%fieldBC(1)%nComponents - select case (component) - case (1) - loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_X_ID - case (2) - loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Y_ID - case (3) - loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Z_ID - end select + do l = 1, load_steps%length + load_step => load_steps%get_dict(l) + step_bc => load_step%get_dict('boundary_conditions') + step_mech => step_bc%get_list('mechanical') + allocate(loadCases(l)%mechBC(mesh_Nboundaries)) + loadCases(l)%mechBC(:)%nComponents = dimPlex !< X, Y (, Z) displacements + do faceSet = 1, mesh_Nboundaries + allocate(loadCases(l)%mechBC(faceSet)%Value(dimPlex), source = 0.0_pREAL) + allocate(loadCases(l)%mechBC(faceSet)%Mask(dimPlex), source = .false.) end do - do component = 1, loadCases(i)%fieldBC(1)%nComponents - allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pREAL) - allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.) + + do m = 1, step_mech%length + mech_BC => step_mech%get_dict(m) + currentFaceSet = -1 + do faceSet = 1, mesh_Nboundaries + if (mesh_boundaries(faceSet) == mech_BC%get_asInt('tag')) currentFaceSet = faceSet + end do + if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC') + mech_u => mech_BC%get_list('dot_u') + do component = 1, dimPlex + if (mech_u%get_asStr(component) /= 'x') then + loadCases(l)%mechBC(currentFaceSet)%Mask(component) = .true. + loadCases(l)%mechBC(currentFaceSet)%Value(component) = mech_u%get_asReal(component) + end if + end do end do - end do + step_discretization => load_step%get_dict('discretization') + loadCases(l)%t = step_discretization%get_asReal('t') + loadCases(l)%N = step_discretization%get_asInt ('N') -!-------------------------------------------------------------------------------------------------- -! reading the load case and assign values to the allocated data structure - do l = 1, size(fileContent) - line = fileContent(l) - if (IO_isBlank(line)) cycle ! skip empty lines - - chunkPos = IO_strPos(line) - do i = 1, chunkPos(1) - select case (IO_strValue(line,chunkPos,i)) -!-------------------------------------------------------------------------------------------------- -! loadcase information - case('$Loadcase') - currentLoadCase = IO_intValue(line,chunkPos,i+1) - case('Face') - currentFace = IO_intValue(line,chunkPos,i+1) - currentFaceSet = -1 - do faceSet = 1, mesh_Nboundaries - if (mesh_boundaries(faceSet) == currentFace) currentFaceSet = faceSet - end do - if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC') - case('t') - loadCases(currentLoadCase)%time = IO_realValue(line,chunkPos,i+1) - case('N') - loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1) - case('f_out') - loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,chunkPos,i+1) - case('estimate_rate') - loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory - -!-------------------------------------------------------------------------------------------------- -! boundary condition information - case('X','Y','Z') - select case(IO_strValue(line,chunkPos,i)) - case('X') - ID = COMPONENT_MECH_X_ID - case('Y') - ID = COMPONENT_MECH_Y_ID - case('Z') - ID = COMPONENT_MECH_Z_ID - end select - - do component = 1, loadcases(currentLoadCase)%fieldBC(1)%nComponents - if (loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%ID == ID) then - loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Mask (currentFaceSet) = & - .true. - loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Value(currentFaceSet) = & - IO_realValue(line,chunkPos,i+1) - end if - end do - end select - end do + if (load_step%get_asStr('f_out',defaultVal='n/a') == 'none') then + loadCases(l)%f_out = huge(0) + else + loadCases(l)%f_out = load_step%get_asInt('f_out', defaultVal=1) + end if + loadCases(l)%estimate_rate = (load_step%get_asBool('estimate_rate',defaultVal=.true.) .and. l>1) end do !-------------------------------------------------------------------------------------------------- ! consistency checks and output of load case - loadCases(1)%followFormerTrajectory = .false. ! cannot guess along trajectory for first inc of first currentLoadCase errorID = 0 - checkLoadcases: do currentLoadCase = 1, size(loadCases) - write (loadcase_string, '(i0)' ) currentLoadCase - print'(/,1x,a,1x,i0)', 'load case:', currentLoadCase - if (.not. loadCases(currentLoadCase)%followFormerTrajectory) & + checkLoadcases: do l = 1, load_steps%length + write (loadcase_string, '(i0)' ) l + print'(/,1x,a,1x,i0)', 'load case:', l + if (.not. loadCases(l)%estimate_rate) & print'(2x,a)', 'drop guessing along trajectory' print'(2x,a)', 'Field '//trim(FIELD_MECH_label) do faceSet = 1, mesh_Nboundaries - do component = 1, loadCases(currentLoadCase)%fieldBC(1)%nComponents - if (loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Mask(faceSet)) & + do component = 1, dimPlex + if (loadCases(l)%mechBC(faceSet)%Mask(component)) & print'(a,i2,a,i2,a,f12.7)', & ' Face ', mesh_boundaries(faceSet), & ' Component ', component, & - ' Value ', loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Value(faceSet) + ' Value ', loadCases(l)%mechBC(faceSet)%Value(component) end do end do - print'(2x,a,f12.6)', 'time: ', loadCases(currentLoadCase)%time - if (loadCases(currentLoadCase)%incs < 1) errorID = 835 ! non-positive incs count - print'(2x,a,i5)', 'increments: ', loadCases(currentLoadCase)%incs - if (loadCases(currentLoadCase)%outputfrequency < 1) errorID = 836 ! non-positive result frequency + print'(2x,a,f12.6)', 'time: ', loadCases(l)%t + if (loadCases(l)%N < 1) errorID = 835 ! non-positive incs count + print'(2x,a,i5)', 'increments: ', loadCases(l)%N + if (loadCases(l)%f_out < 1) errorID = 836 ! non-positive result frequency print'(2x,a,i5)', 'output frequency: ', & - loadCases(currentLoadCase)%outputfrequency - if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message + loadCases(l)%f_out + if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message end do checkLoadcases !-------------------------------------------------------------------------------------------------- ! doing initialization depending on active solvers call FEM_Utilities_init(num_mesh) - call FEM_mechanical_init(loadCases(1)%fieldBC(1),num_mesh) + call FEM_mechanical_init(loadCases(1)%mechBC,num_mesh) call config_numerics_deallocate() if (worldrank == 0) then @@ -244,46 +200,45 @@ program DAMASK_mesh flush(IO_STDOUT) call materialpoint_result(0,0.0_pREAL) - loadCaseLooping: do currentLoadCase = 1, size(loadCases) - time0 = time ! load case start time - guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc + loadCaseLooping: do l = 1, load_steps%length + t_0 = t ! load case start time + guess = loadCases(l)%estimate_rate ! change of load case? homogeneous guess for the first inc - incLooping: do inc = 1, loadCases(currentLoadCase)%incs + incLooping: do inc = 1, loadCases(l)%N totalIncsCounter = totalIncsCounter + 1 !-------------------------------------------------------------------------------------------------- ! forwarding time - timeIncOld = timeinc ! last timeinc that brought former inc to an end - timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pREAL) - timeinc = timeinc * real(subStepFactor,pREAL)**real(-cutBackLevel,pREAL) ! depending on cut back level, decrease time step + Delta_t_prev = Delta_t ! last timeinc that brought former inc to an end + Delta_t = loadCases(l)%t/real(loadCases(l)%N,pREAL) + Delta_t = Delta_t * real(subStepFactor,pREAL)**real(-cutBackLevel,pREAL) ! depending on cut back level, decrease time step stepFraction = 0 ! fraction scaled by stepFactor**cutLevel subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel) - remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time - time = time + timeinc ! forward target time + t = t + Delta_t ! forward target time stepFraction = stepFraction + 1 ! count step !-------------------------------------------------------------------------------------------------- ! report begin of new step print'(/,1x,a)', '###########################################################################' print'(1x,a,es12.5,6(a,i0))',& - 'Time', time, & - 's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,& + 'Time', t, & + 's: Increment ', inc, '/', loadCases(l)%N,& '-', stepFraction, '/', subStepFactor**cutBackLevel,& - ' of load case ', currentLoadCase,'/',size(loadCases) + ' of load case ', l,'/',load_steps%length write(incInfo,'(4(a,i0))') & - 'Increment ',totalIncsCounter,'/',sum(loadCases%incs),& + 'Increment ',totalIncsCounter,'/',sum(loadCases%N),& '-',stepFraction, '/', subStepFactor**cutBackLevel flush(IO_STDOUT) - call FEM_mechanical_forward(guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1)) + call FEM_mechanical_forward(guess,Delta_t,Delta_t_prev,loadCases(l)%mechBC) !-------------------------------------------------------------------------------------------------- ! solve fields stagIter = 0 stagIterate = .true. do while (stagIterate) - solres(1) = FEM_mechanical_solution(incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1)) + solres(1) = FEM_mechanical_solution(incInfo,Delta_t,Delta_t_prev,loadCases(l)%mechBC) if (.not. solres(1)%converged) exit stagIter = stagIter + 1 @@ -294,13 +249,13 @@ program DAMASK_mesh ! check solution cutBack = .False. - if (.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found + if (.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found if (cutBackLevel < maxCutBack) then ! do cut back cutBack = .True. stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator cutBackLevel = cutBackLevel + 1 - time = time - timeinc ! rewind time - timeinc = timeinc/2.0_pREAL + t = t - Delta_t ! rewind time + Delta_t = Delta_t/2.0_pREAL print'(/,1x,a)', 'cutting back' else ! default behavior, exit if spectral solver does not converge if (worldrank == 0) close(statUnit) @@ -308,10 +263,10 @@ program DAMASK_mesh end if else guess = .true. ! start guessing after first converged (sub)inc - timeIncOld = timeinc + Delta_t_prev = Delta_t end if if (.not. cutBack .and. worldrank == 0) then - write(statUnit,*) totalIncsCounter, time, cutBackLevel, & + write(statUnit,*) totalIncsCounter, t, cutBackLevel, & solres%converged, solres%iterationsNeeded ! write statistics about accepted solution flush(statUnit) end if @@ -325,10 +280,10 @@ program DAMASK_mesh print'(/,1x,a,1x,i0,1x,a)', 'increment', totalIncsCounter, 'NOT converged' end if; flush(IO_STDOUT) - if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency + if (mod(inc,loadCases(l)%f_out) == 0) then ! at output frequency print'(/,1x,a)', '... saving results ........................................................' call FEM_mechanical_updateCoords() - call materialpoint_result(totalIncsCounter,time) + call materialpoint_result(totalIncsCounter,t) end if @@ -344,4 +299,5 @@ program DAMASK_mesh call quit(0) ! no complains ;) + end program DAMASK_mesh diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index b1c218172..220852b89 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -18,6 +18,7 @@ module FEM_utilities use math use misc use IO + use parallelization use discretization_mesh use homogenization use FEM_quadrature @@ -38,18 +39,6 @@ module FEM_utilities character(len=*), parameter, public :: & FIELD_MECH_label = 'mechanical' - enum, bind(c); enumerator :: & - FIELD_UNDEFINED_ID, & - FIELD_MECH_ID - end enum - enum, bind(c); enumerator :: & - COMPONENT_UNDEFINED_ID, & - COMPONENT_MECH_X_ID, & - COMPONENT_MECH_Y_ID, & - COMPONENT_MECH_Z_ID - end enum - - !-------------------------------------------------------------------------------------------------- ! derived types type, public :: tSolutionState !< return type of solution from FEM solver variants @@ -58,17 +47,11 @@ module FEM_utilities PetscInt :: iterationsNeeded = 0_pPETSCINT end type tSolutionState - type, public :: tComponentBC - integer(kind(COMPONENT_UNDEFINED_ID)) :: ID + type, public :: tMechBC + integer :: nComponents = 0 real(pREAL), allocatable, dimension(:) :: Value logical, allocatable, dimension(:) :: Mask - end type tComponentBC - - type, public :: tFieldBC - integer(kind(FIELD_UNDEFINED_ID)) :: ID - integer :: nComponents = 0 - type(tComponentBC), allocatable, dimension(:) :: componentBC - end type tFieldBC + end type tMechBC external :: & ! ToDo: write interfaces PetscSectionGetFieldComponents, & @@ -78,12 +61,7 @@ module FEM_utilities public :: & FEM_utilities_init, & utilities_constitutiveResponse, & - utilities_projectBCValues, & - FIELD_MECH_ID, & - COMPONENT_UNDEFINED_ID, & - COMPONENT_MECH_X_ID, & - COMPONENT_MECH_Y_ID, & - COMPONENT_MECH_Z_ID + utilities_projectBCValues contains @@ -142,24 +120,23 @@ end subroutine FEM_utilities_init !-------------------------------------------------------------------------------------------------- !> @brief calculates constitutive response !-------------------------------------------------------------------------------------------------- -subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) +subroutine utilities_constitutiveResponse(Delta_t,P_av,forwardData) - real(pREAL), intent(in) :: timeinc !< loading time + real(pREAL), intent(in) :: Delta_t !< loading time logical, intent(in) :: forwardData !< age results real(pREAL),intent(out), dimension(3,3) :: P_av !< average PK stress integer(MPI_INTEGER_KIND) :: err_MPI + print'(/,1x,a)', '... evaluating constitutive response ......................................' - call homogenization_mechanical_response(timeinc,1,mesh_maxNips*mesh_NcpElems) ! calculate P field - if (.not. terminallyIll) & - call homogenization_mechanical_response2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) + call homogenization_mechanical_response(Delta_t,1,mesh_maxNips*mesh_NcpElems) ! calculate P field cutBack = .false. P_av = sum(homogenization_P,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' + call parallelization_chkerr(err_MPI) end subroutine utilities_constitutiveResponse @@ -168,7 +145,7 @@ end subroutine utilities_constitutiveResponse !-------------------------------------------------------------------------------------------------- !> @brief Project BC values to local vector !-------------------------------------------------------------------------------------------------- -subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCValue,BCDotValue,timeinc) +subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCValue,BCDotValue,Delta_t) Vec :: localVec PetscInt :: field, comp, nBcPoints, point, dof, numDof, numComp, offset @@ -176,7 +153,7 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa IS :: bcPointsIS PetscInt, pointer :: bcPoints(:) real(pREAL), pointer :: localArray(:) - real(pREAL) :: BCValue,BCDotValue,timeinc + real(pREAL) :: BCValue,BCDotValue,Delta_t PetscErrorCode :: err_PETSc @@ -193,7 +170,7 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa call PetscSectionGetFieldOffset(section,bcPoints(point),field,offset,err_PETSc) CHKERRQ(err_PETSc) do dof = offset+comp+1, offset+numDof, numComp - localArray(dof) = localArray(dof) + BCValue + BCDotValue*timeinc + localArray(dof) = localArray(dof) + BCValue + BCDotValue*Delta_t end do end do call VecRestoreArrayF90(localVec,localArray,err_PETSc) diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index be9be3b19..22eb6147f 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -121,13 +121,13 @@ subroutine discretization_mesh_init(restart) CHKERRQ(err_PETSc) mesh_Nboundaries = int(Nboundaries) call MPI_Bcast(mesh_Nboundaries,1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call MPI_Bcast(mesh_NcpElemsGlobal,1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) dim = int(dimPlex) call MPI_Bcast(dim,1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) dimPlex = int(dim,pPETSCINT) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) if (worldsize == 1) then call DMClone(globalMesh,geomMesh,err_PETSc) @@ -149,7 +149,7 @@ subroutine discretization_mesh_init(restart) call ISRestoreIndicesF90(faceSetIS,pFaceSets,err_PETSc) end if call MPI_Bcast(mesh_boundaries,mesh_Nboundaries,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) call DMDestroy(globalMesh,err_PETSc) CHKERRQ(err_PETSc) diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 2d5556e63..13ae970ae 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -36,8 +36,8 @@ module mesh_mechanical_FEM !-------------------------------------------------------------------------------------------------- ! derived types type tSolutionParams - type(tFieldBC) :: fieldBC - real(pREAL) :: timeinc + type(tMechBC),allocatable, dimension(:) :: mechBC + real(pREAL) :: Delta_t end type tSolutionParams type(tSolutionParams) :: params @@ -97,9 +97,9 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields and fills them with data !-------------------------------------------------------------------------------------------------- -subroutine FEM_mechanical_init(fieldBC,num_mesh) +subroutine FEM_mechanical_init(mechBC,num_mesh) - type(tFieldBC), intent(in) :: fieldBC + type(tMechBC), dimension(:), intent(in):: mechBC type(tDict), pointer, intent(in) :: num_mesh DM :: mechanical_mesh @@ -112,7 +112,7 @@ subroutine FEM_mechanical_init(fieldBC,num_mesh) PetscInt, dimension(:), pointer :: pNumComp, pNumDof, pBcField, pBcPoint PetscInt :: numBC, bcSize, nc, & - field, faceSet, topologDim, nNodalPoints, & + component, faceSet, topologDim, nNodalPoints, & cellStart, cellEnd, cell, basis IS :: bcPoint @@ -208,17 +208,17 @@ subroutine FEM_mechanical_init(fieldBC,num_mesh) CHKERRQ(err_PETSc) end do numBC = 0 - do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries - if (fieldBC%componentBC(field)%Mask(faceSet)) numBC = numBC + 1 + do faceSet = 1, mesh_Nboundaries; do component = 1, dimPlex + if (mechBC(faceSet)%Mask(component)) numBC = numBC + 1 end do; end do allocate(pbcField(numBC), source=0_pPETSCINT) allocate(pbcComps(numBC)) allocate(pbcPoints(numBC)) numBC = 0 - do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries - if (fieldBC%componentBC(field)%Mask(faceSet)) then + do faceSet = 1, mesh_Nboundaries; do component = 1, dimPlex + if (mechBC(faceSet)%Mask(component)) then numBC = numBC + 1 - call ISCreateGeneral(PETSC_COMM_WORLD,1_pPETSCINT,[field-1],PETSC_COPY_VALUES,pbcComps(numBC),err_PETSc) + call ISCreateGeneral(PETSC_COMM_WORLD,1_pPETSCINT,[component-1],PETSC_COPY_VALUES,pbcComps(numBC),err_PETSc) CHKERRQ(err_PETSc) call DMGetStratumSize(mechanical_mesh,'Face Sets',mesh_boundaries(faceSet),bcSize,err_PETSc) CHKERRQ(err_PETSc) @@ -320,15 +320,15 @@ end subroutine FEM_mechanical_init !> @brief solution for the FEM load step !-------------------------------------------------------------------------------------------------- type(tSolutionState) function FEM_mechanical_solution( & - incInfoIn,timeinc,timeinc_old,fieldBC) + incInfoIn,Delta_t,Delta_t_prev,mechBC) !-------------------------------------------------------------------------------------------------- ! input data for solution real(pREAL), intent(in) :: & - timeinc, & !< increment in time for current solution - timeinc_old !< increment in time of last increment - type(tFieldBC), intent(in) :: & - fieldBC + Delta_t, & !< increment in time for current solution + Delta_t_prev !< increment in time of last increment + type(tMechBC), dimension(:),intent(in) :: & + mechBC character(len=*), intent(in) :: & incInfoIn @@ -339,8 +339,8 @@ type(tSolutionState) function FEM_mechanical_solution( & FEM_mechanical_solution%converged = .false. !-------------------------------------------------------------------------------------------------- ! set module wide availabe data - params%timeinc = timeinc - params%fieldBC = fieldBC + params%Delta_t = Delta_t + params%mechBC = mechBC call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,err_PETSc) ! solve mechanical_snes based on solution guess (result in solution) CHKERRQ(err_PETSc) @@ -380,7 +380,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc real(pREAL), dimension(cellDof), target :: f_scal PetscReal :: IcellJMat(dimPlex,dimPlex) PetscReal, dimension(:),pointer :: pV0, pCellJ, pInvcellJ, basisField, basisFieldDer - PetscInt :: cellStart, cellEnd, cell, field, face, & + PetscInt :: cellStart, cellEnd, cell, component, face, & qPt, basis, comp, cidx, & numFields, & bcSize,m,i @@ -406,14 +406,14 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc CHKERRQ(err_PETSc) call VecWAXPY(x_local,1.0_pREAL,xx_local,solution_local,err_PETSc) CHKERRQ(err_PETSc) - do field = 1_pPETSCINT, dimPlex; do face = 1_pPETSCINT, mesh_Nboundaries - if (params%fieldBC%componentBC(field)%Mask(face)) then + do face = 1_pPETSCINT, mesh_Nboundaries; do component = 1_pPETSCINT, dimPlex + if (params%mechBC(face)%Mask(component)) then call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc) if (bcSize > 0) then call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc) CHKERRQ(err_PETSc) - call utilities_projectBCValues(x_local,section,0_pPETSCINT,field-1,bcPoints, & - 0.0_pREAL,params%fieldBC%componentBC(field)%Value(face),params%timeinc) + call utilities_projectBCValues(x_local,section,0_pPETSCINT,component-1,bcPoints, & + 0.0_pREAL,params%mechBC(face)%Value(component),params%Delta_t) call ISDestroy(bcPoints,err_PETSc) CHKERRQ(err_PETSc) end if @@ -459,9 +459,9 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call utilities_constitutiveResponse(params%timeinc,P_av,ForwardData) + call utilities_constitutiveResponse(params%Delta_t,P_av,ForwardData) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call parallelization_chkerr(err_MPI) ForwardData = .false. !-------------------------------------------------------------------------------------------------- @@ -529,7 +529,7 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P real(pREAL),dimension(cellDOF,cellDOF), target :: K_e real(pREAL),dimension(cellDOF,cellDOF) :: K_eA, K_eB - PetscInt :: cellStart, cellEnd, cell, field, face, & + PetscInt :: cellStart, cellEnd, cell, component, face, & qPt, basis, comp, cidx,bcSize, m, i IS :: bcPoints @@ -556,14 +556,14 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P CHKERRQ(err_PETSc) call VecWAXPY(x_local,1.0_pREAL,xx_local,solution_local,err_PETSc) CHKERRQ(err_PETSc) - do field = 1, dimPlex; do face = 1, mesh_Nboundaries - if (params%fieldBC%componentBC(field)%Mask(face)) then + do face = 1, mesh_Nboundaries; do component = 1, dimPlex + if (params%mechBC(face)%Mask(component)) then call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc) if (bcSize > 0) then call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc) CHKERRQ(err_PETSc) - call utilities_projectBCValues(x_local,section,0_pPETSCINT,field-1,bcPoints, & - 0.0_pREAL,params%fieldBC%componentBC(field)%Value(face),params%timeinc) + call utilities_projectBCValues(x_local,section,0_pPETSCINT,component-1,bcPoints, & + 0.0_pREAL,params%mechBC(face)%Value(component),params%Delta_t) call ISDestroy(bcPoints,err_PETSc) CHKERRQ(err_PETSc) end if @@ -665,17 +665,17 @@ end subroutine FEM_mechanical_formJacobian !-------------------------------------------------------------------------------------------------- !> @brief forwarding routine !-------------------------------------------------------------------------------------------------- -subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC) +subroutine FEM_mechanical_forward(guess,Delta_t,Delta_t_prev,mechBC) - type(tFieldBC), intent(in) :: & - fieldBC + type(tMechBC), dimension(:), intent(in) :: & + mechBC real(pREAL), intent(in) :: & - timeinc_old, & - timeinc + Delta_t_prev, & + Delta_t logical, intent(in) :: & guess - PetscInt :: field, face, bcSize + PetscInt :: component, face, bcSize DM :: dm_local Vec :: x_local PetscSection :: section @@ -686,7 +686,6 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC) ! 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) @@ -701,14 +700,14 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC) CHKERRQ(err_PETSc) call VecAXPY(solution_local,1.0_pREAL,x_local,err_PETSc) CHKERRQ(err_PETSc) - do field = 1, dimPlex; do face = 1, mesh_Nboundaries - if (fieldBC%componentBC(field)%Mask(face)) then + do face = 1, mesh_Nboundaries; do component = 1, dimPlex + if (mechBC(face)%Mask(component)) then call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc) if (bcSize > 0) then call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc) CHKERRQ(err_PETSc) - call utilities_projectBCValues(solution_local,section,0_pPETSCINT,field-1,bcPoints, & - 0.0_pREAL,fieldBC%componentBC(field)%Value(face),timeinc_old) + call utilities_projectBCValues(solution_local,section,0_pPETSCINT,component-1,bcPoints, & + 0.0_pREAL,mechBC(face)%Value(component),Delta_t_prev) call ISDestroy(bcPoints,err_PETSc) CHKERRQ(err_PETSc) end if @@ -721,12 +720,12 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC) ! update rate and forward last inc call VecCopy(solution,solution_rate,err_PETSc) CHKERRQ(err_PETSc) - call VecScale(solution_rate,timeinc_old**(-1),err_PETSc) + call VecScale(solution_rate,Delta_t_prev**(-1),err_PETSc) CHKERRQ(err_PETSc) end if call VecCopy(solution_rate,solution,err_PETSc) CHKERRQ(err_PETSc) - call VecScale(solution,timeinc,err_PETSc) + call VecScale(solution,Delta_t,err_PETSc) CHKERRQ(err_PETSc) end subroutine FEM_mechanical_forward diff --git a/src/phase.f90 b/src/phase.f90 index ffa606cc2..39a8679be 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -326,14 +326,6 @@ 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 - type(tRotationContainer), dimension(:), intent(in) :: orientation - end subroutine plastic_nonlocal_updateCompatibility - module subroutine plastic_dependentState(ph,en) integer, intent(in) :: & ph, & @@ -366,7 +358,6 @@ module phase phase_allocateState, & phase_forward, & phase_restore, & - plastic_nonlocal_updateCompatibility, & converged, & phase_mechanical_constitutive, & phase_thermal_constitutive, & @@ -387,7 +378,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief Initialize constitutive models for individual physics !-------------------------------------------------------------------------------------------------- -subroutine phase_init +subroutine phase_init() integer :: & ph, ce, co, ma @@ -544,25 +535,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 +554,27 @@ 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) - - 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 +598,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.f90 b/src/phase_mechanical_plastic.f90 index 5d1a462e7..2d0324742 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -204,6 +204,11 @@ submodule(phase:mechanical) plastic en end subroutine plastic_nonlocal_deltaState + module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,en) + integer, intent(in) :: ph, en + type(tRotationContainer), dimension(:), intent(in) :: orientation + end subroutine plastic_nonlocal_updateCompatibility + end interface contains @@ -359,6 +364,7 @@ module subroutine plastic_dependentState(ph,en) case (MECHANICAL_PLASTICITY_NONLOCAL) plasticType call nonlocal_dependentState(ph,en) + if (plasticState(ph)%nonlocal) call plastic_nonlocal_updateCompatibility(phase_O,ph,en) end select plasticType diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 545dec4e6..3af650f1d 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -6,18 +6,18 @@ !-------------------------------------------------------------------------------------------------- submodule(phase:plastic) nonlocal use geometry_plastic_nonlocal, only: & - nIPneighbors => geometry_plastic_nonlocal_nIPneighbors, & + nCellNeighbors => geometry_plastic_nonlocal_nIPneighbors, & IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, & - IPvolume => geometry_plastic_nonlocal_IPvolume0, & - IParea => geometry_plastic_nonlocal_IParea0, & - IPareaNormal => geometry_plastic_nonlocal_IPareaNormal0, & + IPvolume0 => geometry_plastic_nonlocal_IPvolume0, & + IParea0 => geometry_plastic_nonlocal_IParea0, & + IPareaNormal0 => geometry_plastic_nonlocal_IPareaNormal0, & geometry_plastic_nonlocal_disable type :: tGeometry - real(pREAL), dimension(:), allocatable :: V_0 + real(pREAL), dimension(:), allocatable :: v_0 + real(pREAL), dimension(:,:), allocatable :: a_0, x_0 + real(pREAL), dimension(:,:,:), allocatable :: n_0 integer, dimension(:,:,:), allocatable :: IPneighborhood - real(pREAL), dimension(:,:), allocatable :: IParea, IPcoordinates - real(pREAL), dimension(:,:,:), allocatable :: IPareaNormal end type tGeometry type(tGeometry), dimension(:), allocatable :: geom @@ -41,13 +41,6 @@ submodule(phase:plastic) nonlocal mob_scr_pos = 3, & !< mobile screw positive mob_scr_neg = 4 !< mobile screw positive - ! BEGIN DEPRECATED - integer, dimension(:,:,:), allocatable :: & - iRhoU, & !< state indices for unblocked density - iV, & !< state indices for dislocation velocities - iD !< state indices for stable dipole height - !END DEPRECATED - real(pREAL), dimension(:,:,:,:,:,:), allocatable :: & compatibility !< slip system compatibility between en and my neighbors @@ -124,7 +117,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 @@ -132,21 +127,20 @@ submodule(phase:plastic) nonlocal type :: tNonlocalState real(pREAL), pointer, dimension(:,:) :: & rho, & ! < all dislocations - rhoSgl, & - rhoSglMobile, & ! iRhoU + rho_sgl, & + rho_sgl_mob, & rho_sgl_mob_edg_pos, & rho_sgl_mob_edg_neg, & rho_sgl_mob_scr_pos, & rho_sgl_mob_scr_neg, & - rhoSglImmobile, & + rho_sgl_imm, & rho_sgl_imm_edg_pos, & rho_sgl_imm_edg_neg, & rho_sgl_imm_scr_pos, & rho_sgl_imm_scr_neg, & - rhoDip, & + rho_dip, & rho_dip_edg, & rho_dip_scr, & - rho_forest, & gamma, & v, & v_edg_pos, & @@ -177,9 +171,8 @@ module function plastic_nonlocal_init() result(myPlasticity) integer :: & ph, & Nmembers, & - sizeState, sizeDotState, sizeDependentState, sizeDeltaState, & - s1, s2, & - s, t, l + sizeState, sizeDotState, sizeDeltaState, & + s1, s2 real(pREAL), dimension(:,:), allocatable :: & a_nS !< non-Schmid coefficients character(len=:), allocatable :: & @@ -383,26 +376,24 @@ module function plastic_nonlocal_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays Nmembers = count(material_ID_phase == ph) - sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', & - 'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', & - 'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', & - '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 & + sizeDotState = size(['rho_sgl_mob_edg_pos', 'rho_sgl_mob_edg_neg', & + 'rho_sgl_mob_scr_pos', 'rho_sgl_mob_scr_neg', & + 'rho_sgl_imm_edg_pos', 'rho_sgl_imm_edg_neg', & + 'rho_sgl_imm_scr_pos', 'rho_sgl_imm_scr_neg', & + 'rho_dip_edg ', 'rho_dip_scr ', & + '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 - allocate(geom(ph)%V_0(Nmembers)) - allocate(geom(ph)%IPneighborhood(3,nIPneighbors,Nmembers)) - allocate(geom(ph)%IPareaNormal(3,nIPneighbors,Nmembers)) - allocate(geom(ph)%IParea(nIPneighbors,Nmembers)) - allocate(geom(ph)%IPcoordinates(3,Nmembers)) + allocate(geom(ph)%v_0(Nmembers)) + allocate(geom(ph)%a_0(nCellNeighbors,Nmembers)) + allocate(geom(ph)%x_0(3,Nmembers)) + allocate(geom(ph)%n_0(3,nCellNeighbors,Nmembers)) + allocate(geom(ph)%IPneighborhood(3,nCellNeighbors,Nmembers)) call storeGeometry(ph) if (plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) & @@ -414,13 +405,14 @@ module function plastic_nonlocal_init() result(myPlasticity) del%rho => plasticState(ph)%deltaState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) plasticState(ph)%atol(1:10*prm%sum_N_sl) = prm%atol_rho - stt%rhoSgl => plasticState(ph)%state (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) - dot%rhoSgl => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) - del%rhoSgl => plasticState(ph)%deltaState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) + stt%rho_sgl => plasticState(ph)%state (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) + dot%rho_sgl => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) + del%rho_sgl => plasticState(ph)%deltaState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) - stt%rhoSglMobile => plasticState(ph)%state (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) - dot%rhoSglMobile => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) - del%rhoSglMobile => plasticState(ph)%deltaState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) + st0%rho_sgl_mob => plasticState(ph)%state0 (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) + stt%rho_sgl_mob => plasticState(ph)%state (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) + dot%rho_sgl_mob => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) + del%rho_sgl_mob => plasticState(ph)%deltaState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) stt%rho_sgl_mob_edg_pos => plasticState(ph)%state (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:) dot%rho_sgl_mob_edg_pos => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:) @@ -438,9 +430,9 @@ module function plastic_nonlocal_init() result(myPlasticity) dot%rho_sgl_mob_scr_neg => plasticState(ph)%dotState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) del%rho_sgl_mob_scr_neg => plasticState(ph)%deltaState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) - stt%rhoSglImmobile => plasticState(ph)%state (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) - dot%rhoSglImmobile => plasticState(ph)%dotState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) - del%rhoSglImmobile => plasticState(ph)%deltaState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) + stt%rho_sgl_imm => plasticState(ph)%state (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) + dot%rho_sgl_imm => plasticState(ph)%dotState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) + del%rho_sgl_imm => plasticState(ph)%deltaState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) stt%rho_sgl_imm_edg_pos => plasticState(ph)%state (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:) dot%rho_sgl_imm_edg_pos => plasticState(ph)%dotState (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:) @@ -458,9 +450,9 @@ module function plastic_nonlocal_init() result(myPlasticity) dot%rho_sgl_imm_scr_neg => plasticState(ph)%dotState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) del%rho_sgl_imm_scr_neg => plasticState(ph)%deltaState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) - stt%rhoDip => plasticState(ph)%state (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:) - dot%rhoDip => plasticState(ph)%dotState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:) - del%rhoDip => plasticState(ph)%deltaState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:) + stt%rho_dip => plasticState(ph)%state (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:) + dot%rho_dip => plasticState(ph)%dotState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:) + del%rho_dip => plasticState(ph)%deltaState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:) stt%rho_dip_edg => plasticState(ph)%state (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:) dot%rho_dip_edg => plasticState(ph)%dotState (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:) @@ -477,16 +469,18 @@ 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%compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,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),nCellNeighbors,Nmembers),source=0.0_pREAL) end associate if (Nmembers > 0) call stateInit(ini,ph,Nmembers) @@ -497,44 +491,9 @@ module function plastic_nonlocal_init() result(myPlasticity) end do - allocate(compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,& + allocate(compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nCellNeighbors,& discretization_nIPs,discretization_Nelems), source=0.0_pREAL) -! 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 - - if (.not. myPlasticity(ph)) cycle - - phase => phases%get_dict(ph) - Nmembers = count(material_ID_phase == ph) - l = 0 - do t = 1,4 - do s = 1,param(ph)%sum_N_sl - l = l + 1 - 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 - 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) & - error stop 'state indices not properly set (nonlocal)' - end do - end function plastic_nonlocal_init @@ -548,9 +507,9 @@ module subroutine nonlocal_dependentState(ph, en) en integer :: & - no, & !< neighbor offset - neighbor_el, & ! element number of neighboring material point - neighbor_ip, & ! integration point of neighboring material point + en_nbr, & ! neighbor phase entry + el_nbr, & ! element number of neighboring material point + ip_nbr, & ! integration point of neighboring material point c, & ! index of dilsocation character (edge, screw) s, & ! slip system index dir, & @@ -565,7 +524,7 @@ module subroutine nonlocal_dependentState(ph, en) real(pREAL), dimension(2) :: & rhoExcessGradient, & rhoExcessGradient_over_rho, & - rhoTotal + rho_sum real(pREAL), dimension(3) :: & rhoExcessDifferences, & normal_latticeConf @@ -574,25 +533,23 @@ module subroutine nonlocal_dependentState(ph, en) invFp, & !< inverse of plastic deformation gradient connections, & invConnections - real(pREAL), dimension(3,nIPneighbors) :: & + real(pREAL), dimension(3,nCellNeighbors) :: & connection_latticeConf - real(pREAL), dimension(2,param(ph)%sum_N_sl) :: & - rhoExcess real(pREAL), dimension(param(ph)%sum_N_sl) :: & rho_edg_delta, & rho_scr_delta real(pREAL), dimension(param(ph)%sum_N_sl,10) :: & rho, & - rho0, & - rho_neighbor0 + rho_0, & + rho_0_nbr real(pREAL), dimension(param(ph)%sum_N_sl,param(ph)%sum_N_sl) :: & myInteractionMatrix ! corrected slip interaction matrix - real(pREAL), dimension(param(ph)%sum_N_sl,nIPneighbors) :: & - rho_edg_delta_neighbor, & - rho_scr_delta_neighbor - real(pREAL), dimension(2,maxval(param%sum_N_sl),nIPneighbors) :: & - neighbor_rhoExcess, & ! excess density at neighboring material point - neighbor_rhoTotal ! total density at neighboring material point + real(pREAL), dimension(param(ph)%sum_N_sl,nCellNeighbors) :: & + rho_edg_delta_nbr, & + rho_scr_delta_nbr + real(pREAL), dimension(2,maxval(param%sum_N_sl),nCellNeighbors) :: & + rho_delta_nbr, & ! excess density at neighboring material point + rho_sum_nbr ! total density at neighboring material point real(pREAL), dimension(3,param(ph)%sum_N_sl,2) :: & m ! direction of dislocation motion @@ -602,7 +559,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 +569,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 @@ -628,60 +585,57 @@ module subroutine nonlocal_dependentState(ph, en) ! ToDo: MD: this is most likely only correct for F_i = I !################################################################################################# - rho0 = getRho0(ph,en) + rho_0 = getRho0(ph,en) if (plasticState(ph)%nonlocal .and. prm%shortRangeStressCorrection) then invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en)) invFe = math_inv33(phase_mechanical_Fe(ph)%data(1:3,1:3,en)) - rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg) - rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg) + rho_edg_delta = rho_0(:,mob_edg_pos) - rho_0(:,mob_edg_neg) + rho_scr_delta = rho_0(:,mob_scr_pos) - rho_0(:,mob_scr_neg) - rhoExcess(1,:) = rho_edg_delta - rhoExcess(2,:) = rho_scr_delta - - FVsize = geom(ph)%V_0(en)**(1.0_pREAL/3.0_pREAL) + FVsize = geom(ph)%v_0(en)**(1.0_pREAL/3.0_pREAL) !* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities nRealNeighbors = 0.0_pREAL - neighbor_rhoTotal = 0.0_pREAL - do n = 1,nIPneighbors - neighbor_el = geom(ph)%IPneighborhood(1,n,en) - neighbor_ip = geom(ph)%IPneighborhood(2,n,en) + rho_sum_nbr = 0.0_pREAL + do n = 1,nCellNeighbors + el_nbr = geom(ph)%IPneighborhood(1,n,en) + ip_nbr = geom(ph)%IPneighborhood(2,n,en) - if (neighbor_el > 0 .and. neighbor_ip > 0) then - if (material_ID_phase(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip) == ph) then - no = material_entry_phase(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip) + if (el_nbr > 0 .and. ip_nbr > 0) then + if (material_ID_phase(1,(el_nbr-1)*discretization_nIPs + ip_nbr) == ph) then + en_nbr = material_entry_phase(1,(el_nbr-1)*discretization_nIPs + ip_nbr) nRealNeighbors = nRealNeighbors + 1.0_pREAL - rho_neighbor0 = getRho0(ph,no) + rho_0_nbr = getRho0(ph,en_nbr) - rho_edg_delta_neighbor(:,n) = rho_neighbor0(:,mob_edg_pos) - rho_neighbor0(:,mob_edg_neg) - rho_scr_delta_neighbor(:,n) = rho_neighbor0(:,mob_scr_pos) - rho_neighbor0(:,mob_scr_neg) + rho_edg_delta_nbr(:,n) = rho_0_nbr(:,mob_edg_pos) - rho_0_nbr(:,mob_edg_neg) + rho_scr_delta_nbr(:,n) = rho_0_nbr(:,mob_scr_pos) - rho_0_nbr(:,mob_scr_neg) - neighbor_rhoTotal(1,:,n) = sum(abs(rho_neighbor0(:,edg)),2) - neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor0(:,scr)),2) + rho_sum_nbr(1,:,n) = sum(abs(rho_0_nbr(:,edg)),2) + rho_sum_nbr(2,:,n) = sum(abs(rho_0_nbr(:,scr)),2) - connection_latticeConf(1:3,n) = matmul(invFe, geom(ph)%IPcoordinates(1:3,no) & - - geom(ph)%IPcoordinates(1:3,en)) - normal_latticeConf = matmul(transpose(invFp), geom(ph)%IPareaNormal(1:3,n,en)) + connection_latticeConf(1:3,n) = matmul(invFe, geom(ph)%x_0(1:3,en_nbr) & + - geom(ph)%x_0(1:3,en)) + normal_latticeConf = matmul(transpose(invFp), geom(ph)%n_0(1:3,n,en)) if (math_inner(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pREAL) & ! neighboring connection points in opposite direction to face normal: must be periodic image - connection_latticeConf(1:3,n) = normal_latticeConf * geom(ph)%V_0(en)/geom(ph)%IParea(n,en) ! instead take the surface normal scaled with the diameter of the cell + connection_latticeConf(1:3,n) = normal_latticeConf * geom(ph)%v_0(en)/geom(ph)%a_0(n,en) ! instead take the surface normal scaled with the diameter of the cell else ! local neighbor or different lattice structure or different constitution instance -> use central values instead connection_latticeConf(1:3,n) = 0.0_pREAL - rho_edg_delta_neighbor(:,n) = rho_edg_delta - rho_scr_delta_neighbor(:,n) = rho_scr_delta + rho_edg_delta_nbr(:,n) = rho_edg_delta + rho_scr_delta_nbr(:,n) = rho_scr_delta end if else ! free surface -> use central values instead connection_latticeConf(1:3,n) = 0.0_pREAL - rho_edg_delta_neighbor(:,n) = rho_edg_delta - rho_scr_delta_neighbor(:,n) = rho_scr_delta + rho_edg_delta_nbr(:,n) = rho_edg_delta + rho_scr_delta_nbr(:,n) = rho_scr_delta end if end do - neighbor_rhoExcess(1,:,:) = rho_edg_delta_neighbor - neighbor_rhoExcess(2,:,:) = rho_scr_delta_neighbor + rho_delta_nbr(1,:,:) = rho_edg_delta_nbr + rho_delta_nbr(2,:,:) = rho_scr_delta_nbr !* loop through the slip systems and calculate the dislocation gradient by !* 1. interpolation of the excess density in the neighorhood @@ -698,8 +652,8 @@ module subroutine nonlocal_dependentState(ph, en) neighbors(2) = 2 * dir connections(dir,1:3) = connection_latticeConf(1:3,neighbors(1)) & - connection_latticeConf(1:3,neighbors(2)) - rhoExcessDifferences(dir) = neighbor_rhoExcess(c,s,neighbors(1)) & - - neighbor_rhoExcess(c,s,neighbors(2)) + rhoExcessDifferences(dir) = rho_delta_nbr(c,s,neighbors(1)) & + - rho_delta_nbr(c,s,neighbors(2)) end do invConnections = math_inv33(connections) if (all(dEq0(invConnections))) call IO_error(-1,ext_msg='back stress calculation: inversion error') @@ -712,11 +666,11 @@ module subroutine nonlocal_dependentState(ph, en) rhoExcessGradient(2) = rhoExcessGradient(2) + sum(rho(s,imm_scr)) / FVsize ! ... normalized with the total density ... - rhoTotal(1) = (sum(abs(rho(s,edg))) + sum(neighbor_rhoTotal(1,s,:))) / (1.0_pREAL + nRealNeighbors) - rhoTotal(2) = (sum(abs(rho(s,scr))) + sum(neighbor_rhoTotal(2,s,:))) / (1.0_pREAL + nRealNeighbors) + rho_sum(1) = (sum(abs(rho(s,edg))) + sum(rho_sum_nbr(1,s,:))) / (1.0_pREAL + nRealNeighbors) + rho_sum(2) = (sum(abs(rho(s,scr))) + sum(rho_sum_nbr(2,s,:))) / (1.0_pREAL + nRealNeighbors) rhoExcessGradient_over_rho = 0.0_pREAL - where(rhoTotal > 0.0_pREAL) rhoExcessGradient_over_rho = rhoExcessGradient / rhoTotal + where(rho_sum > 0.0_pREAL) rhoExcessGradient_over_rho = rhoExcessGradient / rho_sum ! ... gives the local stress correction when multiplied with a factor dst%tau_back(s,en) = - mu * prm%b_sl(s) / (2.0_pREAL * PI) & @@ -750,7 +704,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & t, & !< dislocation type s !< index of my current slip system real(pREAL), dimension(param(ph)%sum_N_sl,8) :: & - rhoSgl !< single dislocation densities (including blocked) + rho_sgl !< single dislocation densities (including blocked) real(pREAL), dimension(param(ph)%sum_N_sl,10) :: & rho real(pREAL), dimension(param(ph)%sum_N_sl,4) :: & @@ -773,7 +727,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & !*** shortcut to state variables rho = getRho(ph,en) - rhoSgl = rho(:,sgl) + rho_sgl = rho(:,sgl) do s = 1,prm%sum_N_sl tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) @@ -806,20 +760,20 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & stt%v(:,en) = pack(v,.true.) !*** Bauschinger effect - forall (s = 1:prm%sum_N_sl, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pREAL) & - rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t)) + forall (s = 1:prm%sum_N_sl, t = 5:8, rho_sgl(s,t) * v(s,t-4) < 0.0_pREAL) & + rho_sgl(s,t-4) = rho_sgl(s,t-4) + abs(rho_sgl(s,t)) - dot_gamma = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl + dot_gamma = sum(rho_sgl(:,1:4) * v, 2) * prm%b_sl do s = 1,prm%sum_N_sl Lp = Lp + dot_gamma(s) * prm%P_sl(1:3,1:3,s) forall (i=1:3,j=1:3,k=1:3,l=1:3) & dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) & + prm%P_sl(i,j,s) * prm%P_sl(k,l,s) & - * sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%b_sl(s) & + * sum(rho_sgl(s,1:4) * dv_dtau(s,1:4)) * prm%b_sl(s) & + prm%P_sl(i,j,s) & - * (+ prm%P_nS_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) & - - prm%P_nS_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s) + * (+ prm%P_nS_pos(k,l,s) * rho_sgl(s,3) * dv_dtauNS(s,3) & + - prm%P_nS_neg(k,l,s) * rho_sgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s) end do end associate @@ -855,23 +809,23 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en) real(pREAL), dimension(param(ph)%sum_N_sl) :: & tau ! current resolved shear stress real(pREAL), dimension(param(ph)%sum_N_sl,2) :: & - rhoDip, & ! current dipole dislocation densities (screw and edge dipoles) + rho_dip, & ! current dipole dislocation densities (screw and edge dipoles) dUpper, & ! current maximum stable dipole distance for edges and screws dUpperOld, & ! old maximum stable dipole distance for edges and screws deltaDUpper ! change in maximum stable dipole distance for edges and screws - associate(prm => param(ph),dst => dependentState(ph),del => deltaState(ph)) + 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) + rho_dip = rho(:,dip) !**************************************************************************** !*** dislocation remobilization (bauschinger effect) @@ -911,11 +865,11 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en) deltaRhoDipole2SingleStress = 0.0_pREAL forall (c=1:2, s=1:prm%sum_N_sl, deltaDUpper(s,c) < 0.0_pREAL .and. & dNeq0(dUpperOld(s,c) - prm%minDipoleHeight(s,c))) & - deltaRhoDipole2SingleStress(s,8+c) = rhoDip(s,c) * deltaDUpper(s,c) & + deltaRhoDipole2SingleStress(s,8+c) = rho_dip(s,c) * deltaDUpper(s,c) & / (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]) @@ -941,28 +895,24 @@ module subroutine nonlocal_dotState(Mp,timestep, & integer :: & c, & !< character of dislocation - t, & !< type of dislocation s !< index of my current slip system real(pREAL), dimension(param(ph)%sum_N_sl,10) :: & rho, & - rho0, & !< dislocation density at beginning of time step rhoDot, & !< density evolution rhoDotMultiplication, & !< density evolution by multiplication rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide) rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation rhoDotThermalAnnihilation !< density evolution by thermal annihilation real(pREAL), dimension(param(ph)%sum_N_sl,8) :: & - rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) - my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) + rho_sgl !< current single dislocation densities (positive/negative screw and edge without dipoles) real(pREAL), dimension(param(ph)%sum_N_sl,4) :: & v, & !< current dislocation glide velocity - v0, & dot_gamma !< shear rates real(pREAL), dimension(param(ph)%sum_N_sl) :: & tau, & !< current resolved shear stress v_climb !< climb velocity of edge dipoles real(pREAL), dimension(param(ph)%sum_N_sl,2) :: & - rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) + rho_dip, & !< current dipole dislocation densities (screw and edge dipoles) dLower, & !< minimum stable dipole distance for edges and screws dUpper !< current maximum stable dipole distance for edges and screws real(pREAL) :: & @@ -975,7 +925,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) @@ -984,15 +935,12 @@ module subroutine nonlocal_dotState(Mp,timestep, & tau = 0.0_pREAL dot_gamma = 0.0_pREAL - rho = getRho(ph,en) - rhoSgl = rho(:,sgl) - rhoDip = rho(:,dip) - 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) - dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) + rho = getRho(ph,en) + rho_sgl = rho(:,sgl) + rho_dip = rho(:,dip) + v = reshape(stt%v(:,en),[prm%sum_N_sl,4]) + dot_gamma = rho_sgl(:,1:4) * v * spread(prm%b_sl,2,4) ! limits for stable dipole height @@ -1018,21 +966,19 @@ 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) - !**************************************************************************** ! dipole formation and annihilation @@ -1040,20 +986,20 @@ module subroutine nonlocal_dotState(Mp,timestep, & ! formation by glide do c = 1,2 rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pREAL * dUpper(:,c) / prm%b_sl & - * ( rhoSgl(:,2*c-1) * abs(dot_gamma(:,2*c)) & ! negative mobile --> positive mobile - + rhoSgl(:,2*c) * abs(dot_gamma(:,2*c-1)) & ! positive mobile --> negative mobile - + abs(rhoSgl(:,2*c+4)) * abs(dot_gamma(:,2*c-1))) ! positive mobile --> negative immobile + * ( rho_sgl(:,2*c-1) * abs(dot_gamma(:,2*c)) & ! negative mobile --> positive mobile + + rho_sgl(:,2*c) * abs(dot_gamma(:,2*c-1)) & ! positive mobile --> negative mobile + + abs(rho_sgl(:,2*c+4)) * abs(dot_gamma(:,2*c-1))) ! positive mobile --> negative immobile rhoDotSingle2DipoleGlide(:,2*c) = -2.0_pREAL * dUpper(:,c) / prm%b_sl & - * ( rhoSgl(:,2*c-1) * abs(dot_gamma(:,2*c)) & ! negative mobile --> positive mobile - + rhoSgl(:,2*c) * abs(dot_gamma(:,2*c-1)) & ! positive mobile --> negative mobile - + abs(rhoSgl(:,2*c+3)) * abs(dot_gamma(:,2*c))) ! negative mobile --> positive immobile + * ( rho_sgl(:,2*c-1) * abs(dot_gamma(:,2*c)) & ! negative mobile --> positive mobile + + rho_sgl(:,2*c) * abs(dot_gamma(:,2*c-1)) & ! positive mobile --> negative mobile + + abs(rho_sgl(:,2*c+3)) * abs(dot_gamma(:,2*c))) ! negative mobile --> positive immobile rhoDotSingle2DipoleGlide(:,2*c+3) = -2.0_pREAL * dUpper(:,c) / prm%b_sl & - * rhoSgl(:,2*c+3) * abs(dot_gamma(:,2*c)) ! negative mobile --> positive immobile + * rho_sgl(:,2*c+3) * abs(dot_gamma(:,2*c)) ! negative mobile --> positive immobile rhoDotSingle2DipoleGlide(:,2*c+4) = -2.0_pREAL * dUpper(:,c) / prm%b_sl & - * rhoSgl(:,2*c+4) * abs(dot_gamma(:,2*c-1)) ! positive mobile --> negative immobile + * rho_sgl(:,2*c+4) * abs(dot_gamma(:,2*c-1)) ! positive mobile --> negative immobile rhoDotSingle2DipoleGlide(:,c+8) = abs(rhoDotSingle2DipoleGlide(:,2*c+3)) & + abs(rhoDotSingle2DipoleGlide(:,2*c+4)) & @@ -1066,15 +1012,15 @@ module subroutine nonlocal_dotState(Mp,timestep, & rhoDotAthermalAnnihilation = 0.0_pREAL forall (c=1:2) & rhoDotAthermalAnnihilation(:,c+8) = -2.0_pREAL * dLower(:,c) / prm%b_sl & - * ( 2.0_pREAL * (rhoSgl(:,2*c-1) * abs(dot_gamma(:,2*c)) + rhoSgl(:,2*c) * abs(dot_gamma(:,2*c-1))) & ! was single hitting single - + 2.0_pREAL * (abs(rhoSgl(:,2*c+3)) * abs(dot_gamma(:,2*c)) + abs(rhoSgl(:,2*c+4)) * abs(dot_gamma(:,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single - + rhoDip(:,c) * (abs(dot_gamma(:,2*c-1)) + abs(dot_gamma(:,2*c)))) ! single knocks dipole constituent + * ( 2.0_pREAL * (rho_sgl(:,2*c-1) * abs(dot_gamma(:,2*c)) + rho_sgl(:,2*c) * abs(dot_gamma(:,2*c-1))) & ! was single hitting single + + 2.0_pREAL * (abs(rho_sgl(:,2*c+3)) * abs(dot_gamma(:,2*c)) + abs(rho_sgl(:,2*c+4)) * abs(dot_gamma(:,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single + + rho_dip(:,c) * (abs(dot_gamma(:,2*c-1)) + abs(dot_gamma(:,2*c)))) ! single knocks dipole constituent ! annihilated screw dipoles leave edge jogs behind on the colinear system 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 @@ -1083,8 +1029,8 @@ module subroutine nonlocal_dotState(Mp,timestep, & v_climb = D_SD * mu * prm%V_at & / (PI * (1.0_pREAL-nu) * (dUpper(:,1) + dLower(:,1)) * K_B * Temperature) ! eq. 3.54 forall (s = 1:prm%sum_N_sl, dUpper(s,1) > dLower(s,1)) & - rhoDotThermalAnnihilation(s,9) = max(- 4.0_pREAL * rhoDip(s,1) * v_climb(s) / (dUpper(s,1) - dLower(s,1)), & - - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & + rhoDotThermalAnnihilation(s,9) = max(- 4.0_pREAL * rho_dip(s,1) * v_climb(s) / (dUpper(s,1) - dLower(s,1)), & + - rho_dip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have rhoDot = rhoDotFlux(timestep, ph,en) & @@ -1125,38 +1071,36 @@ function rhoDotFlux(timestep,ph,en) ns, & !< short notation for the total number of active slip systems c, & !< character of dislocation n, & !< index of my current neighbor - neighbor_el, & !< element number of my neighbor - neighbor_ip, & !< integration point of my neighbor - neighbor_n, & !< neighbor index pointing to en when looking from my neighbor + el_nbr, & !< element number of my neighbor + ip_nbr, & !< integration point of my neighbor + n_nbr, & !< neighbor index pointing to en when looking from my neighbor opposite_neighbor, & !< index of my opposite neighbor opposite_ip, & !< ip of my opposite neighbor opposite_el, & !< element index of my opposite neighbor opposite_n, & !< neighbor index pointing to en when looking from my opposite neighbor t, & !< type of dislocation - no,& !< neighbor offset shortcut - np,& !< neighbor phase shortcut + en_nbr,& !< neighbor phase entry + ph_nbr,& !< neighbor phase ID topp, & !< type of dislocation with opposite sign to t s !< index of my current slip system real(pREAL), dimension(param(ph)%sum_N_sl,10) :: & rho, & - rho0, & !< dislocation density at beginning of time step + rho_0, & !< dislocation density at beginning of time step rhoDotFlux !< density evolution by flux - real(pREAL), dimension(param(ph)%sum_N_sl,8) :: & - rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) - neighbor_rhoSgl0, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) - my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) real(pREAL), dimension(param(ph)%sum_N_sl,4) :: & - v, & !< current dislocation glide velocity - v0, & - neighbor_v0, & !< dislocation glide velocity of enighboring ip - dot_gamma !< shear rates + rho_0_sgl_mob, & !< mobile dislocation densities of neighboring ip (positive/negative screw and edge) + rho_0_sgl_mob_nbr, & !< mobile dislocation densities of neighboring ip (positive/negative screw and edge) + v, & !< dislocation glide velocity + v_0, & + v_0_nbr, & !< dislocation glide velocity of enighboring ip + dot_gamma !< shear rates real(pREAL), dimension(3,param(ph)%sum_N_sl,4) :: & m !< direction of dislocation motion real(pREAL), dimension(3,3) :: & my_F, & !< my total deformation gradient - neighbor_F, & !< total deformation gradient of my neighbor + F_nbr, & !< total deformation gradient of my neighbor my_Fe, & !< my elastic deformation gradient - neighbor_Fe, & !< elastic deformation gradient of my neighbor + F_e_nbr, & !< elastic deformation gradient of my neighbor Favg !< average total deformation gradient of en and my neighbor real(pREAL), dimension(3) :: & normal_neighbor2me, & !< interface normal pointing from my neighbor to en in neighbor's lattice configuration @@ -1164,27 +1108,28 @@ function rhoDotFlux(timestep,ph,en) normal_me2neighbor, & !< interface normal pointing from en to my neighbor in my lattice configuration normal_me2neighbor_defConf !< interface normal pointing from en to my neighbor in shared deformed configuration real(pREAL) :: & - area, & !< area of the current interface + a, & !< area of the current interface transmissivity, & !< overall transmissivity of dislocation flux to neighboring material point lineLength !< dislocation line length leaving the current interface 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 - rho = getRho(ph,en) - rhoSgl = rho(:,sgl) - rho0 = getRho0(ph,en) - my_rhoSgl0 = rho0(:,sgl) + rho = getRho(ph,en) + rho_0 = getRho0(ph,en) + rho_0_sgl_mob = rho_0(:,mob) - 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 - dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) + v = reshape(stt%v(:,en),[prm%sum_N_sl,4]) !ToDo: MD: I think we should use state0 here + dot_gamma = rho(:,mob) * 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) + v_0 = reshape(st0%v(:,en),[prm%sum_N_sl,4]) !**************************************************************************** !*** calculate dislocation fluxes (only for nonlocal plasticity) @@ -1193,8 +1138,8 @@ function rhoDotFlux(timestep,ph,en) !*** check CFL (Courant-Friedrichs-Lewy) condition for flux if (any( abs(dot_gamma) > 0.0_pREAL & ! any active slip system ... - .and. prm%C_CFL * abs(v0) * timestep & - > geom(ph)%V_0(en)/ maxval(geom(ph)%IParea(:,en)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) + .and. prm%C_CFL * abs(v_0) * timestep & + > geom(ph)%v_0(en)/ maxval(geom(ph)%a_0(:,en)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) rhoDotFlux = IEEE_value(1.0_pREAL,IEEE_quiet_NaN) ! enforce cutback return end if @@ -1211,28 +1156,28 @@ function rhoDotFlux(timestep,ph,en) my_F = phase_mechanical_F(ph)%data(1:3,1:3,en) my_Fe = matmul(my_F, math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en))) - neighbors: do n = 1,nIPneighbors + neighbors: do n = 1,nCellNeighbors - neighbor_el = geom(ph)%IPneighborhood(1,n,en) - neighbor_ip = geom(ph)%IPneighborhood(2,n,en) - neighbor_n = geom(ph)%IPneighborhood(3,n,en) - np = material_ID_phase(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip) - no = material_entry_phase(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip) + el_nbr = geom(ph)%IPneighborhood(1,n,en) + ip_nbr = geom(ph)%IPneighborhood(2,n,en) + n_nbr = geom(ph)%IPneighborhood(3,n,en) + ph_nbr = material_ID_phase(1,(el_nbr-1)*discretization_nIPs + ip_nbr) + en_nbr = material_entry_phase(1,(el_nbr-1)*discretization_nIPs + ip_nbr) opposite_neighbor = n + mod(n,2) - mod(n+1,2) opposite_el = geom(ph)%IPneighborhood(1,opposite_neighbor,en) opposite_ip = geom(ph)%IPneighborhood(2,opposite_neighbor,en) opposite_n = geom(ph)%IPneighborhood(3,opposite_neighbor,en) - if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient - neighbor_F = phase_mechanical_F(np)%data(1:3,1:3,no) - neighbor_Fe = matmul(neighbor_F, math_inv33(phase_mechanical_Fp(np)%data(1:3,1:3,no))) - Favg = 0.5_pREAL * (my_F + neighbor_F) + if (n_nbr > 0) then ! if neighbor exists, average deformation gradient + F_nbr = phase_mechanical_F(ph_nbr)%data(1:3,1:3,en_nbr) + F_e_nbr = matmul(F_nbr, math_inv33(phase_mechanical_Fp(ph_nbr)%data(1:3,1:3,en_nbr))) + Favg = 0.5_pREAL * (my_F + F_nbr) else ! if no neighbor, take my value as average Favg = my_F end if - neighbor_v0 = 0.0_pREAL ! needed for check of sign change in flux density below + v_0_nbr = 0.0_pREAL ! needed for check of sign change in flux density below !* FLUX FROM MY NEIGHBOR TO ME !* This is only considered, if I have a neighbor of nonlocal plasticity @@ -1242,38 +1187,37 @@ function rhoDotFlux(timestep,ph,en) !* my neighbor's interface. !* The entering flux from my neighbor will be distributed on my slip systems according to the !* compatibility - if (neighbor_n > 0) then - if (mechanical_plasticity_type(np) == MECHANICAL_PLASTICITY_NONLOCAL .and. & + if (n_nbr > 0) then + if (mechanical_plasticity_type(ph_nbr) == MECHANICAL_PLASTICITY_NONLOCAL .and. & any(dependentState(ph)%compatibility(:,:,:,n,en) > 0.0_pREAL)) then - forall (s = 1:ns, t = 1:4) - neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,np),no) - neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,np),no),0.0_pREAL) - endforall + ! ToDo MD: Not sure if ns is correct here, but I think that compatibility is 0 if different phase + v_0_nbr = reshape(state0(ph_nbr)%v(:,en_nbr),[ns,4]) + rho_0_sgl_mob_nbr = max(reshape(state0(ph_nbr)%rho_sgl_mob(:,en_nbr),[ns,4]),0.0_pREAL) - where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pREAL < prm%rho_min & - .or. neighbor_rhoSgl0 < prm%rho_significant) & - neighbor_rhoSgl0 = 0.0_pREAL + where (rho_0_sgl_mob_nbr * geom(ph_nbr)%v_0(en_nbr) ** 0.667_pREAL < prm%rho_min & + .or. rho_0_sgl_mob_nbr < prm%rho_significant) & + rho_0_sgl_mob_nbr = 0.0_pREAL normal_neighbor2me_defConf = math_det33(Favg) * matmul(math_inv33(transpose(Favg)), & - IPareaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! normal of the interface in (average) deformed configuration (pointing neighbor => en) - normal_neighbor2me = matmul(transpose(neighbor_Fe), normal_neighbor2me_defConf) & - / math_det33(neighbor_Fe) ! interface normal in the lattice configuration of my neighbor - area = IParea(neighbor_n,neighbor_ip,neighbor_el) * norm2(normal_neighbor2me) + geom(ph_nbr)%n_0(1:3,n_nbr,en_nbr)) ! normal of the interface in (average) deformed configuration (pointing neighbor => en) + normal_neighbor2me = matmul(transpose(F_e_nbr), normal_neighbor2me_defConf) & + / math_det33(F_e_nbr) ! interface normal in the lattice configuration of my neighbor + a = geom(ph_nbr)%a_0(n_nbr,en_nbr) * norm2(normal_neighbor2me) normal_neighbor2me = normal_neighbor2me / norm2(normal_neighbor2me) ! normalize the surface normal to unit length do s = 1,ns do t = 1,4 c = (t + 1) / 2 topp = t + mod(t,2) - mod(t+1,2) - if (neighbor_v0(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_pREAL & ! flux from my neighbor to en == entering flux for en - .and. v0(s,t) * neighbor_v0(s,t) >= 0.0_pREAL ) then ! ... only if no sign change in flux density - lineLength = neighbor_rhoSgl0(s,t) * neighbor_v0(s,t) & - * math_inner(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface + if (v_0_nbr(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_pREAL & ! flux from my neighbor to en == entering flux for en + .and. v_0(s,t) * v_0_nbr(s,t) >= 0.0_pREAL ) then ! ... only if no sign change in flux density + lineLength = rho_0_sgl_mob_nbr(s,t) * v_0_nbr(s,t) & + * math_inner(m(1:3,s,t), normal_neighbor2me) * a ! positive line length that wants to enter through this interface where (dependentState(ph)%compatibility(c,:,s,n,en) > 0.0_pREAL) & rhoDotFlux(:,t) = rhoDotFlux(1:ns,t) & - + lineLength/geom(ph)%V_0(en)*dependentState(ph)%compatibility(c,:,s,n,en)**2 ! transferring to equally signed mobile dislocation type + + lineLength/geom(ph)%v_0(en)*dependentState(ph)%compatibility(c,:,s,n,en)**2 ! transferring to equally signed mobile dislocation type where (dependentState(ph)%compatibility(c,:,s,n,en) < 0.0_pREAL) & rhoDotFlux(:,topp) = rhoDotFlux(:,topp) & - + lineLength/geom(ph)%V_0(en)*dependentState(ph)%compatibility(c,:,s,n,en)**2 ! transferring to opposite signed mobile dislocation type + + lineLength/geom(ph)%v_0(en)*dependentState(ph)%compatibility(c,:,s,n,en)**2 ! transferring to opposite signed mobile dislocation type end if end do @@ -1289,29 +1233,29 @@ function rhoDotFlux(timestep,ph,en) !* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density. !* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations. if (opposite_n > 0) then - if (mechanical_plasticity_type(np) == MECHANICAL_PLASTICITY_NONLOCAL) then + if (mechanical_plasticity_type(ph_nbr) == MECHANICAL_PLASTICITY_NONLOCAL) then normal_me2neighbor_defConf = math_det33(Favg) & - * matmul(math_inv33(transpose(Favg)),geom(ph)%IPareaNormal(1:3,n,en)) ! normal of the interface in (average) deformed configuration (pointing en => neighbor) + * matmul(math_inv33(transpose(Favg)),geom(ph)%n_0(1:3,n,en)) ! normal of the interface in (average) deformed configuration (pointing en => neighbor) normal_me2neighbor = matmul(transpose(my_Fe), normal_me2neighbor_defConf) & / math_det33(my_Fe) ! interface normal in my lattice configuration - area = geom(ph)%IParea(n,en) * norm2(normal_me2neighbor) + a = geom(ph)%a_0(n,en) * norm2(normal_me2neighbor) normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor) ! normalize the surface normal to unit length do s = 1,ns do t = 1,4 c = (t + 1) / 2 - if (v0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pREAL ) then ! flux from en to my neighbor == leaving flux for en (might also be a pure flux from my mobile density to dead density if interface not at all transmissive) - if (v0(s,t) * neighbor_v0(s,t) >= 0.0_pREAL) then ! no sign change in flux density + if (v_0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pREAL ) then ! flux from en to my neighbor == leaving flux for en (might also be a pure flux from my mobile density to dead density if interface not at all transmissive) + if (v_0(s,t) * v_0_nbr(s,t) >= 0.0_pREAL) then ! no sign change in flux density transmissivity = sum(dependentState(ph)%compatibility(c,:,s,n,en)**2) ! overall transmissivity from this slip system to my neighbor else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor transmissivity = 0.0_pREAL end if - lineLength = my_rhoSgl0(s,t) * v0(s,t) & - * math_inner(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface - rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / geom(ph)%V_0(en) ! subtract dislocation flux from current type + lineLength = rho_0_sgl_mob(s,t) * v_0(s,t) & + * math_inner(m(1:3,s,t), normal_me2neighbor) * a ! positive line length of mobiles that wants to leave through this interface + rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / geom(ph)%v_0(en) ! subtract dislocation flux from current type rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) & - + lineLength / geom(ph)%V_0(en) * (1.0_pREAL - transmissivity) & - * sign(1.0_pREAL, v0(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point + + lineLength / geom(ph)%v_0(en) * (1.0_pREAL - transmissivity) & + * sign(1.0_pREAL, v_0(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point end if end do end do @@ -1331,26 +1275,24 @@ 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,ph,en) type(tRotationContainer), dimension(:), intent(in) :: & orientation ! crystal orientation integer, intent(in) :: & - ph, & - ip, & - el + ph, en integer :: & n, & ! neighbor index - en, & - neighbor_e, & ! element index of my neighbor - neighbor_i, & ! integration point index of my neighbor - neighbor_me, & - neighbor_phase, & + el_nbr, & ! element index of my neighbor + ip_nbr, & ! integration point index of my neighbor + ce_nbr, & + en_nbr, & + ph_nbr, & 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(ph)%sum_N_sl,param(ph)%sum_N_sl,nCellNeighbors) :: & my_compatibility ! my_compatibility for current element and ip real(pREAL) :: & my_compatibilitySum, & @@ -1364,29 +1306,26 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el) associate(prm => param(ph)) ns = prm%sum_N_sl - en = material_entry_phase(1,(el-1)*discretization_nIPs + ip) !*** start out fully compatible my_compatibility = 0.0_pREAL forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pREAL - neighbors: do n = 1,nIPneighbors - neighbor_e = IPneighborhood(1,n,ip,el) - neighbor_i = IPneighborhood(2,n,ip,el) - neighbor_me = material_entry_phase(1,(neighbor_e-1)*discretization_nIPs + neighbor_i) - neighbor_phase = material_ID_phase(1,(neighbor_e-1)*discretization_nIPs + neighbor_i) + neighbors: do n = 1,nCellNeighbors + el_nbr = geom(ph)%IPneighborhood(1,n,en) + ip_nbr = geom(ph)%IPneighborhood(2,n,en) + ce_nbr = (el_nbr-1)*discretization_nIPs + ip_nbr + en_nbr = material_entry_phase(1,ce_nbr) + ph_nbr = material_ID_phase(1,ce_nbr) - if (neighbor_e <= 0 .or. neighbor_i <= 0) then - !* FREE SURFACE + if (ce_nbr <= 0) then ! free surface forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_surface) - elseif (neighbor_phase /= ph) then - !* PHASE BOUNDARY - if (plasticState(neighbor_phase)%nonlocal .and. plasticState(ph)%nonlocal) & + elseif (ph_nbr /= ph) then ! phase boundary + if (plasticState(ph_nbr)%nonlocal .and. plasticState(ph)%nonlocal) & forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_pREAL - elseif (prm%chi_GB >= 0.0_pREAL) then - !* GRAIN BOUNDARY + elseif (prm%chi_GB >= 0.0_pREAL) then ! grain boundary if (any(dNeq(phase_O_0(ph)%data(en)%asQuaternion(), & - phase_O_0(neighbor_phase)%data(neighbor_me)%asQuaternion())) .and. & - plasticState(neighbor_phase)%nonlocal) & + phase_O_0(ph_nbr)%data(en_nbr)%asQuaternion())) .and. & + plasticState(ph_nbr)%nonlocal) & forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB) else !* GRAIN BOUNDARY ? @@ -1398,7 +1337,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el) !* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one. !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one. !* All values below the threshold are set to zero. - mis = orientation(ph)%data(en)%misorientation(orientation(neighbor_phase)%data(neighbor_me)) + mis = orientation(ph)%data(en)%misorientation(orientation(ph_nbr)%data(en_nbr)) mySlipSystems: do s1 = 1,ns neighborSlipSystems: do s2 = 1,ns my_compatibility(1,s2,s1,n) = math_inner(prm%slip_normal(1:3,s1), & @@ -1432,7 +1371,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el) end do neighbors - dependentState(ph)%compatibility(:,:,:,:,material_entry_phase(1,(el-1)*discretization_nIPs + ip)) = my_compatibility + dependentState(ph)%compatibility(:,:,:,:,en) = my_compatibility end associate @@ -1486,7 +1425,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)), & @@ -1544,8 +1483,8 @@ subroutine stateInit(ini,phase,Nentries) associate(stt => state(phase)) if (ini%random_rho_u > 0.0_pREAL) then ! randomly distribute dislocation segments on random slip system and of random type in the volume - totalVolume = sum(geom(phase)%V_0) - minimumIPVolume = minval(geom(phase)%V_0) + totalVolume = sum(geom(phase)%v_0) + minimumIPVolume = minval(geom(phase)%v_0) densityBinning = ini%random_rho_u_binning / minimumIpVolume ** (2.0_pREAL / 3.0_pREAL) ! fill random material points with dislocation segments until the desired overall density is reached @@ -1554,8 +1493,8 @@ subroutine stateInit(ini,phase,Nentries) call random_number(rnd) e = nint(rnd(1)*real(Nentries,pREAL) + 0.5_pREAL) s = nint(rnd(2)*real(sum(ini%N_sl),pREAL)*4.0_pREAL + 0.5_pREAL) - meanDensity = meanDensity + densityBinning * geom(phase)%V_0(e) / totalVolume - stt%rhoSglMobile(s,e) = densityBinning + meanDensity = meanDensity + densityBinning * geom(phase)%v_0(e) / totalVolume + stt%rho_sgl_mob(s,e) = densityBinning end do else ! homogeneous distribution with noise do f = 1,size(ini%N_sl,1) @@ -1685,7 +1624,7 @@ pure function getRho(ph,en) result(rho) rho(:,mob) = max(rho(:,mob),0.0_pREAL) rho(:,dip) = max(rho(:,dip),0.0_pREAL) - where(abs(rho) < max(prm%rho_min/geom(ph)%V_0(en)**(2.0_pREAL/3.0_pREAL),prm%rho_significant)) & + where(abs(rho) < max(prm%rho_min/geom(ph)%v_0(en)**(2.0_pREAL/3.0_pREAL),prm%rho_significant)) & rho = 0.0_pREAL end associate @@ -1697,56 +1636,57 @@ end function getRho !> @brief returns copy of current dislocation densities from state !> @details raw values is rectified !-------------------------------------------------------------------------------------------------- -pure function getRho0(ph,en) result(rho0) +pure function getRho0(ph,en) result(rho_0) integer, intent(in) :: ph, en - real(pREAL), dimension(param(ph)%sum_N_sl,10) :: rho0 + real(pREAL), dimension(param(ph)%sum_N_sl,10) :: rho_0 associate(prm => param(ph)) - rho0 = reshape(state0(ph)%rho(:,en),[prm%sum_N_sl,10]) + rho_0 = reshape(state0(ph)%rho(:,en),[prm%sum_N_sl,10]) ! ensure positive densities (not for imm, they have a sign) - rho0(:,mob) = max(rho0(:,mob),0.0_pREAL) - rho0(:,dip) = max(rho0(:,dip),0.0_pREAL) + rho_0(:,mob) = max(rho_0(:,mob),0.0_pREAL) + rho_0(:,dip) = max(rho_0(:,dip),0.0_pREAL) - where (abs(rho0) < max(prm%rho_min/geom(ph)%V_0(en)**(2.0_pREAL/3.0_pREAL),prm%rho_significant)) & - rho0 = 0.0_pREAL + where (abs(rho_0) < max(prm%rho_min/geom(ph)%v_0(en)**(2.0_pREAL/3.0_pREAL),prm%rho_significant)) & + rho_0 = 0.0_pREAL end associate end function getRho0 +!-------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- subroutine storeGeometry(ph) integer, intent(in) :: ph - integer :: ce, co, nCell - real(pREAL), dimension(:), allocatable :: V + integer :: ce, nCell + real(pREAL), dimension(:), allocatable :: v_0 + real(pREAL), dimension(:,:), allocatable :: a_0, x_0 + real(pREAL), dimension(:,:,:), allocatable :: n_0 integer, dimension(:,:,:), allocatable :: neighborhood - real(pREAL), dimension(:,:), allocatable :: area, coords - real(pREAL), dimension(:,:,:), allocatable :: areaNormal - nCell = product(shape(IPvolume)) - V = reshape(IPvolume,[nCell]) - neighborhood = reshape(IPneighborhood,[3,nIPneighbors,nCell]) - area = reshape(IParea,[nIPneighbors,nCell]) - areaNormal = reshape(IPareaNormal,[3,nIPneighbors,nCell]) - coords = reshape(discretization_IPcoords,[3,nCell]) + nCell = product(shape(IPVolume0)) + + v_0 = reshape(IPVolume0,[nCell]) + a_0 = reshape(IPArea0,[nCellNeighbors,nCell]) + x_0 = reshape(discretization_IPcoords,[3,nCell]) + n_0 = reshape(IPAreaNormal0,[3,nCellNeighbors,nCell]) + neighborhood = reshape(IPneighborhood,[3,nCellNeighbors,nCell]) do ce = 1, size(material_entry_homogenization,1) - do co = 1, homogenization_maxNconstituents - if (material_ID_phase(co,ce) == ph) then - geom(ph)%V_0(material_entry_phase(co,ce)) = V(ce) - geom(ph)%IPneighborhood(:,:,material_entry_phase(co,ce)) = neighborhood(:,:,ce) - geom(ph)%IParea(:,material_entry_phase(co,ce)) = area(:,ce) - geom(ph)%IPareaNormal(:,:,material_entry_phase(co,ce)) = areaNormal(:,:,ce) - geom(ph)%IPcoordinates(:,material_entry_phase(co,ce)) = coords(:,ce) - end if - end do + if (material_ID_phase(1,ce) == ph) then + geom(ph)%v_0(material_entry_phase(1,ce)) = v_0(ce) + geom(ph)%a_0(:,material_entry_phase(1,ce)) = a_0(:,ce) + geom(ph)%x_0(:,material_entry_phase(1,ce)) = x_0(:,ce) + geom(ph)%n_0(:,:,material_entry_phase(1,ce)) = n_0(:,:,ce) + geom(ph)%IPneighborhood(:,:,material_entry_phase(1,ce)) = neighborhood(:,:,ce) + end if end do end subroutine storeGeometry