diff --git a/PRIVATE b/PRIVATE index ab64793bb..95f7faea9 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit ab64793bb04c506d815ebc850672ed0f2d013e67 +Subproject commit 95f7faea920dd6956884e4a55f72e5d5b1ffcdc8 diff --git a/VERSION b/VERSION index 7942fa57e..69aed17b7 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha3-40-ga1c46b445 +v3.0.0-alpha3-70-g7774ca721 diff --git a/env/DAMASK.sh b/env/DAMASK.sh index edf019921..c74c447ed 100644 --- a/env/DAMASK.sh +++ b/env/DAMASK.sh @@ -30,7 +30,6 @@ PATH=${DAMASK_ROOT}/bin:$PATH SOLVER=$(type -p DAMASK_grid || true 2>/dev/null) [ "x$SOLVER" == "x" ] && SOLVER=$(blink 'Not found!') -[ "x$OMP_NUM_THREADS" == "x" ] && OMP_NUM_THREADS=1 # currently, there is no information that unlimited stack size causes problems # still, http://software.intel.com/en-us/forums/topic/501500 suggest to fix it @@ -55,9 +54,8 @@ if [ ! -z "$PS1" ]; then [[ $(canonicalPath "$PETSC_DIR") == $PETSC_DIR ]] \ || echo " ~~> "$(canonicalPath "$PETSC_DIR") fi - echo -n "MSC.Marc/Mentat " - [ -d $MSC_ROOT ] && echo $MSC_ROOT || blink $MSC_ROOT - echo + [ "x$PETSC_ARCH" != "x" ] && echo "PETSc architecture $PETSC_ARCH" + [ "x$OMP_NUM_THREADS" == "x" ] && export OMP_NUM_THREADS=4 echo "Multithreading OMP_NUM_THREADS=$OMP_NUM_THREADS" echo -n "heap size " [[ "$(ulimit -d)" == "unlimited" ]] \ diff --git a/env/DAMASK.zsh b/env/DAMASK.zsh index 07064eb7a..5813d110a 100644 --- a/env/DAMASK.zsh +++ b/env/DAMASK.zsh @@ -20,7 +20,6 @@ PATH=${DAMASK_ROOT}/bin:$PATH SOLVER=$(which DAMASK_grid || true 2>/dev/null) [[ "x$SOLVER" == "x" ]] && SOLVER=$(blink 'Not found!') -[[ "x$OMP_NUM_THREADS" == "x" ]] && OMP_NUM_THREADS=1 # currently, there is no information that unlimited stack size causes problems # still, http://software.intel.com/en-us/forums/topic/501500 suggest to fix it @@ -45,11 +44,8 @@ if [ ! -z "$PS1" ]; then [[ $(canonicalPath "$PETSC_DIR") == $PETSC_DIR ]] \ || echo " ~~> "$(canonicalPath "$PETSC_DIR") fi - [[ "x$PETSC_ARCH" == "x" ]] \ - || echo "PETSc architecture $PETSC_ARCH" - echo -n "MSC.Marc/Mentat " - [ -d $MSC_ROOT ] && echo $MSC_ROOT || blink $MSC_ROOT - echo + [[ "x$PETSC_ARCH" != "x" ]] && echo "PETSc architecture $PETSC_ARCH" + [[ "x$OMP_NUM_THREADS" == "x" ]] && export OMP_NUM_THREADS=4 echo "Multithreading OMP_NUM_THREADS=$OMP_NUM_THREADS" echo -n "heap size " [[ "$(ulimit -d)" == "unlimited" ]] \ diff --git a/examples/config/phase/thermal/Al.yaml b/examples/config/phase/thermal/Al.yaml new file mode 100644 index 000000000..8837f30af --- /dev/null +++ b/examples/config/phase/thermal/Al.yaml @@ -0,0 +1,7 @@ +references: + - www.engineeringtoolbox.com/thermal-conductivity-metals-d_858.html + - www.engineeringtoolbox.com/specific-heat-metals-d_152.html +c_p: 0.91e3 +K_11: 236.0 +K_22: 236.0 +K_33: 236.0 diff --git a/examples/config/phase/thermal/Steel-0.5C.yaml b/examples/config/phase/thermal/Steel-0.5C.yaml new file mode 100644 index 000000000..861cbbee5 --- /dev/null +++ b/examples/config/phase/thermal/Steel-0.5C.yaml @@ -0,0 +1,7 @@ +references: + - www.engineeringtoolbox.com/thermal-conductivity-metals-d_858.html + - www.engineeringtoolbox.com/specific-heat-metals-d_152.html +c_p: 0.49e3 +K_11: 54.0 +K_22: 54.0 +K_33: 54.0 diff --git a/installation/mods_MarcMentat/apply_DAMASK_modifications.py b/installation/mods_MarcMentat/apply_DAMASK_modifications.py index 5b6b7a2ce..38323b3f7 100755 --- a/installation/mods_MarcMentat/apply_DAMASK_modifications.py +++ b/installation/mods_MarcMentat/apply_DAMASK_modifications.py @@ -7,6 +7,11 @@ from pathlib import Path import damask +def copy_and_replace(in_file,dst,msc_root,editor): + with open(in_file) as f_in, open(dst/Path(in_file).name,'w') as f_out: + f_out.write(f_in.read().replace('%INSTALLDIR%',msc_root).replace('%EDITOR%',editor)) + + parser = argparse.ArgumentParser( description='Apply DAMASK modification to MSC.Marc/Mentat', prog = Path(__file__).name, @@ -17,34 +22,25 @@ parser.add_argument('--editor', dest='editor', metavar='string', default='vi', parser.add_argument('--msc-root', dest='msc_root', metavar='string', default=damask.solver._marc._msc_root, help='MSC.Marc/Mentat root directory') -parser.add_argument('--msc-version', dest='msc_version', metavar='string', +parser.add_argument('--msc-version', dest='msc_version', type=float, metavar='float', default=damask.solver._marc._msc_version, help='MSC.Marc/Mentat version') -parser.add_argument('--damask-root', dest='damask_root', type=float, metavar = 'float', +parser.add_argument('--damask-root', dest='damask_root', metavar = 'string', default=damask.solver._marc._damask_root, help='DAMASK root directory') - args = parser.parse_args() msc_root = Path(args.msc_root) damask_root = Path(args.damask_root) msc_version = args.msc_version -def copy_and_replace(in_file,dst): - with open(in_file) as f: - content = f.read() - content = content.replace('%INSTALLDIR%',str(msc_root)) - content = content.replace('%EDITOR%', args.editor) - with open(dst/Path(in_file).name,'w') as f: - f.write(content) - print('adapting Marc tools...\n') src = damask_root/f'installation/mods_MarcMentat/{msc_version}/Marc_tools' dst = msc_root/f'marc{msc_version}/tools' for in_file in glob.glob(str(src/'*damask*')) + [str(src/'include_linux64')]: - copy_and_replace(in_file,dst) + copy_and_replace(in_file,dst,args.msc_root,args.editor) print('adapting Mentat scripts and menus...\n') @@ -52,12 +48,12 @@ print('adapting Mentat scripts and menus...\n') src = damask_root/f'installation/mods_MarcMentat/{msc_version}/Mentat_bin' dst = msc_root/f'mentat{msc_version}/bin' for in_file in glob.glob(str(src/'*[!.original]')): - copy_and_replace(in_file,dst) + copy_and_replace(in_file,dst,args.msc_root,args.editor) src = damask_root/f'installation/mods_MarcMentat/{msc_version}/Mentat_menus' dst = msc_root/f'mentat{msc_version}/menus' for in_file in glob.glob(str(src/'job_run.ms')): - copy_and_replace(in_file,dst) + copy_and_replace(in_file,dst,args.msc_root,args.editor) print('compiling Mentat menu binaries...') diff --git a/python/damask/_colormap.py b/python/damask/_colormap.py index 02e2b922c..7e82853b4 100644 --- a/python/damask/_colormap.py +++ b/python/damask/_colormap.py @@ -110,13 +110,13 @@ class Colormap(mpl.colors.ListedColormap): low_,high_ = map(Colormap._rgb2msh,low_high) elif model.lower() == 'hsv': - if np.any(low_high<0) or np.any(low_high[:,1:3]>1) or np.any(low_high[:,0]>360): + if np.any(low_high<0) or np.any(low_high>[360,1,1]): raise ValueError(f'HSV color {low} | {high} are out of range.') low_,high_ = map(Colormap._hsv2msh,low_high) elif model.lower() == 'hsl': - if np.any(low_high<0) or np.any(low_high[:,1:3]>1) or np.any(low_high[:,0]>360): + if np.any(low_high<0) or np.any(low_high>[360,1,1]): raise ValueError(f'HSL color {low} | {high} are out of range.') low_,high_ = map(Colormap._hsl2msh,low_high) diff --git a/python/damask/_grid.py b/python/damask/_grid.py index 1b58870e4..237d56ba1 100644 --- a/python/damask/_grid.py +++ b/python/damask/_grid.py @@ -392,7 +392,7 @@ class Grid: seeds_p = seeds coords = grid_filters.coordinates0_point(cells,size).reshape(-1,3) - pool = mp.Pool(int(os.environ.get('OMP_NUM_THREADS',1))) + pool = mp.Pool(int(os.environ.get('OMP_NUM_THREADS',4))) result = pool.map_async(partial(Grid._find_closest_seed,seeds_p,weights_p), [coord for coord in coords]) pool.close() pool.join() @@ -436,8 +436,12 @@ class Grid: """ coords = grid_filters.coordinates0_point(cells,size).reshape(-1,3) - KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds) - devNull,material_ = KDTree.query(coords) + tree = spatial.cKDTree(seeds,boxsize=size) if periodic else \ + spatial.cKDTree(seeds) + try: + material_ = tree.query(coords, workers = int(os.environ.get('OMP_NUM_THREADS',4)))[1] + except TypeError: + material_ = tree.query(coords, n_jobs = int(os.environ.get('OMP_NUM_THREADS',4)))[1] # scipy <1.6 return Grid(material = (material_ if material is None else material[material_]).reshape(cells), size = size, diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 6d4b09d5b..919d9e518 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -134,30 +134,15 @@ class Orientation(Rotation): """ Rotation.__init__(self) if rotation is None else Rotation.__init__(self,rotation=rotation) - if ( lattice not in _lattice_symmetries - and lattice not in _crystal_families): - raise KeyError(f'Lattice "{lattice}" is unknown') - - self.family = None - self.lattice = None - self.a = None - self.b = None - self.c = None - self.alpha = None - self.beta = None - self.gamma = None self.kinematics = None if lattice in _lattice_symmetries: self.family = _lattice_symmetries[lattice] self.lattice = lattice + self.a = 1 if a is None else a self.b = b self.c = c - self.alpha = (np.radians(alpha) if degrees else alpha) if alpha is not None else None - self.beta = (np.radians(beta) if degrees else beta) if beta is not None else None - self.gamma = (np.radians(gamma) if degrees else gamma) if gamma is not None else None - self.a = float(self.a) if self.a is not None else \ (self.b / self.ratio['b'] if self.b is not None and self.ratio['b'] is not None else self.c / self.ratio['c'] if self.c is not None and self.ratio['c'] is not None else None) @@ -169,9 +154,13 @@ class Orientation(Rotation): (self.a * self.ratio['c'] if self.a is not None and self.ratio['c'] is not None else self.b / self.ratio['b'] * self.ratio['c'] if self.c is not None and self.ratio['b'] is not None and self.ratio['c'] is not None else None) - self.alpha = self.alpha if self.alpha is not None else self.immutable['alpha'] if 'alpha' in self.immutable else None - self.beta = self.beta if self.beta is not None else self.immutable['beta'] if 'beta' in self.immutable else None - self.gamma = self.gamma if self.gamma is not None else self.immutable['gamma'] if 'gamma' in self.immutable else None + + self.alpha = np.radians(alpha) if degrees and alpha is not None else alpha + self.beta = np.radians(beta) if degrees and beta is not None else beta + self.gamma = np.radians(gamma) if degrees and gamma is not None else gamma + if self.alpha is None and 'alpha' in self.immutable: self.alpha = self.immutable['alpha'] + if self.beta is None and 'beta' in self.immutable: self.beta = self.immutable['beta'] + if self.gamma is None and 'gamma' in self.immutable: self.gamma = self.immutable['gamma'] if \ (self.a is None) \ @@ -197,7 +186,13 @@ class Orientation(Rotation): {'direction':self.Bravais_to_Miller(uvtw=master[m][:,0:4]), 'plane': self.Bravais_to_Miller(hkil=master[m][:,4:8])} elif lattice in _crystal_families: - self.family = lattice + self.family = lattice + self.lattice = None + + self.a = self.b = self.c = None + self.alpha = self.beta = self.gamma = None + else: + raise KeyError(f'Lattice "{lattice}" is unknown') def __repr__(self): diff --git a/python/damask/_result.py b/python/damask/_result.py index 0520f08f2..8dd915d1b 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -1179,7 +1179,7 @@ class Result: """ chunk_size = 1024**2//8 - pool = mp.Pool(int(os.environ.get('OMP_NUM_THREADS',1))) + pool = mp.Pool(int(os.environ.get('OMP_NUM_THREADS',4))) lock = mp.Manager().Lock() groups = [] diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py index 029abab19..a40a5990f 100644 --- a/python/damask/_vtk.py +++ b/python/damask/_vtk.py @@ -150,9 +150,8 @@ class VTK: ---------- fname : str or pathlib.Path Filename for reading. Valid extensions are .vtr, .vtu, .vtp, and .vtk. - dataset_type : str, optional + dataset_type : {'vtkRectilinearGrid', 'vtkUnstructuredGrid', 'vtkPolyData'}, optional Name of the vtk.vtkDataSet subclass when opening a .vtk file. - Valid types are vtkRectilinearGrid, vtkUnstructuredGrid, and vtkPolyData. Returns ------- diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index 6ba785df1..c8519abad 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -124,9 +124,6 @@ def strain(F,t,m): """ Calculate strain tensor (Seth–Hill family). - For details refer to https://en.wikipedia.org/wiki/Finite_strain_theory and - https://de.wikipedia.org/wiki/Verzerrungstensor - Parameters ---------- F : numpy.ndarray of shape (...,3,3) @@ -142,6 +139,11 @@ def strain(F,t,m): epsilon : numpy.ndarray of shape (...,3,3) Strain of F. + References + ---------- + https://en.wikipedia.org/wiki/Finite_strain_theory + https://de.wikipedia.org/wiki/Verzerrungstensor + """ if t == 'V': w,n = _np.linalg.eigh(deformation_Cauchy_Green_left(F)) @@ -150,7 +152,6 @@ def strain(F,t,m): if m > 0.0: eps = 1.0/(2.0*abs(m)) * (+ _np.einsum('...j,...kj,...lj',w**m,n,n) - _np.eye(3)) - elif m < 0.0: eps = 1.0/(2.0*abs(m)) * (- _np.einsum('...j,...kj,...lj',w**m,n,n) + _np.eye(3)) else: diff --git a/python/damask/seeds.py b/python/damask/seeds.py index 02a050e66..b2e91a352 100644 --- a/python/damask/seeds.py +++ b/python/damask/seeds.py @@ -3,8 +3,8 @@ from scipy import spatial as _spatial import numpy as _np -from . import util -from . import grid_filters +from . import util as _util +from . import grid_filters as _grid_filters def from_random(size,N_seeds,cells=None,rng_seed=None): @@ -34,7 +34,7 @@ def from_random(size,N_seeds,cells=None,rng_seed=None): if cells is None: coords = rng.random((N_seeds,3)) * size else: - grid_coords = grid_filters.coordinates0_point(cells,size).reshape(-1,3,order='F') + grid_coords = _grid_filters.coordinates0_point(cells,size).reshape(-1,3,order='F') coords = grid_coords[rng.choice(_np.prod(cells),N_seeds, replace=False)] \ + _np.broadcast_to(size/cells,(N_seeds,3))*(rng.random((N_seeds,3))*.5-.25) # wobble without leaving cells @@ -73,12 +73,12 @@ def from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=True,rng_seed= s = 1 i = 0 - progress = util._ProgressBar(N_seeds+1,'',50) + progress = _util._ProgressBar(N_seeds+1,'',50) while s < N_seeds: candidates = rng.random((N_candidates,3))*_np.broadcast_to(size,(N_candidates,3)) tree = _spatial.cKDTree(coords[:s],boxsize=size) if periodic else \ _spatial.cKDTree(coords[:s]) - distances, dev_null = tree.query(candidates) + distances = tree.query(candidates)[0] best = distances.argmax() if distances[best] > distance: # require minimum separation coords[s] = candidates[best] # maximum separation to existing point cloud @@ -120,7 +120,7 @@ def from_grid(grid,selection=None,invert=False,average=False,periodic=True): material = grid.material.reshape((-1,1),order='F') mask = _np.full(grid.cells.prod(),True,dtype=bool) if selection is None else \ _np.isin(material,selection,invert=invert).flatten() - coords = grid_filters.coordinates0_point(grid.cells,grid.size).reshape(-1,3,order='F') + coords = _grid_filters.coordinates0_point(grid.cells,grid.size).reshape(-1,3,order='F') if not average: return (coords[mask],material[mask]) diff --git a/python/damask/util.py b/python/damask/util.py index fcf41c251..efd493f28 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -1,3 +1,5 @@ +"""Miscellaneous helper functionality.""" + import sys import datetime import os @@ -177,6 +179,16 @@ def execute(cmd,wd='./',env=None): def natural_sort(key): + """ + Natural sort. + + For use in python's 'sorted'. + + References + ---------- + https://en.wikipedia.org/wiki/Natural_sort_order + + """ convert = lambda text: int(text) if text.isdigit() else text return [ convert(c) for c in re.split('([0-9]+)', key) ] @@ -191,9 +203,9 @@ def show_progress(iterable,N_iter=None,prefix='',bar_length=50): ---------- iterable : iterable/function with yield statement Iterable (or function with yield statement) to be decorated. - N_iter : int + N_iter : int, optional Total # of iterations. Needed if number of iterations can not be obtained as len(iterable). - prefix : str, optional. + prefix : str, optional Prefix string. bar_length : int, optional Character length of bar. Defaults to 50. @@ -509,6 +521,7 @@ def dict_prune(d): v = dict_prune(v) if not isinstance(v,dict) or v != {}: new[k] = v + return new diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90 index ce00c4913..4ec7b8ea9 100644 --- a/src/HDF5_utilities.f90 +++ b/src/HDF5_utilities.f90 @@ -1775,7 +1775,7 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, & integer, dimension(worldsize) :: writeSize !< contribution of all processes integer(HID_T) :: dcpl - integer :: ierr, hdferr + integer :: ierr, hdferr, HDF5_major, HDF5_minor, HDF5_release integer(HSIZE_T), parameter :: chunkSize = 1024_HSIZE_T**2/8_HSIZE_T !------------------------------------------------------------------------------------------------- @@ -1808,14 +1808,18 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, & call h5pcreate_f(H5P_DATASET_CREATE_F, dcpl, hdferr) if(hdferr < 0) error stop 'HDF5 error' if(product(totalShape) >= chunkSize*2_HSIZE_T) then - call h5pset_chunk_f(dcpl, size(totalShape), getChunks(totalShape,chunkSize), hdferr) - if(hdferr < 0) error stop 'HDF5 error' - call h5pset_shuffle_f(dcpl, hdferr) - if(hdferr < 0) error stop 'HDF5 error' - call h5pset_deflate_f(dcpl, 6, hdferr) - if(hdferr < 0) error stop 'HDF5 error' - call h5pset_Fletcher32_f(dcpl,hdferr) - if(hdferr < 0) error stop 'HDF5 error' + call H5get_libversion_f(HDF5_major,HDF5_minor,HDF5_release,hdferr) + if (hdferr < 0) error stop 'HDF5 error' + if (HDF5_major == 1 .and. HDF5_minor >= 12) then ! https://forum.hdfgroup.org/t/6186 + call h5pset_chunk_f(dcpl, size(totalShape), getChunks(totalShape,chunkSize), hdferr) + if (hdferr < 0) error stop 'HDF5 error' + call h5pset_shuffle_f(dcpl, hdferr) + if (hdferr < 0) error stop 'HDF5 error' + call h5pset_deflate_f(dcpl, 6, hdferr) + if (hdferr < 0) error stop 'HDF5 error' + call h5pset_Fletcher32_f(dcpl,hdferr) + if (hdferr < 0) error stop 'HDF5 error' + endif endif !-------------------------------------------------------------------------------------------------- diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 5ddc7219e..b0bb96402 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -21,6 +21,7 @@ #include "lattice.f90" #include "phase.f90" #include "phase_mechanical.f90" +#include "phase_mechanical_elastic.f90" #include "phase_mechanical_plastic.f90" #include "phase_mechanical_plastic_none.f90" #include "phase_mechanical_plastic_isotropic.f90" diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 1ccb3b901..5a1745668 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -90,7 +90,8 @@ subroutine grid_thermal_spectral_init(T_0) !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc - call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-thermal_snes_type ngmres',ierr) + call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-thermal_snes_type newtonls -thermal_snes_mf & + &-thermal_snes_ksp_ew -thermal_ksp_type fgmres',ierr) CHKERRQ(ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr) CHKERRQ(ierr) diff --git a/src/phase.f90 b/src/phase.f90 index 93743d64a..927ba6b2b 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -58,10 +58,6 @@ module phase grain end type tDebugOptions - integer, dimension(:), allocatable, public :: & !< ToDo: should be protected (bug in Intel compiler) - phase_elasticityInstance, & - phase_NstiffnessDegradations - logical, dimension(:), allocatable, public :: & ! ToDo: should be protected (bug in Intel Compiler) phase_localPlasticity !< flags phases with local constitutive law @@ -298,7 +294,6 @@ module phase end interface - type(tDebugOptions) :: debugConstitutive #if __INTEL_COMPILER >= 1900 public :: & diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 298ff57f3..29e9e8ada 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -5,10 +5,6 @@ submodule(phase) mechanical enum, bind(c); enumerator :: & - ELASTICITY_UNDEFINED_ID, & - ELASTICITY_HOOKE_ID, & - STIFFNESS_DEGRADATION_UNDEFINED_ID, & - STIFFNESS_DEGRADATION_DAMAGE_ID, & PLASTICITY_UNDEFINED_ID, & PLASTICITY_NONE_ID, & PLASTICITY_ISOTROPIC_ID, & @@ -23,11 +19,6 @@ submodule(phase) mechanical KINEMATICS_THERMAL_EXPANSION_ID end enum - integer(kind(ELASTICITY_UNDEFINED_ID)), dimension(:), allocatable :: & - phase_elasticity !< elasticity of each phase - integer(kind(STIFFNESS_DEGRADATION_UNDEFINED_ID)), dimension(:,:), allocatable :: & - phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase - type(tTensorContainer), dimension(:), allocatable :: & ! current value phase_mechanical_Fe, & @@ -57,9 +48,27 @@ submodule(phase) mechanical class(tNode), pointer :: phases end subroutine eigendeformation_init + module subroutine elastic_init(phases) + class(tNode), pointer :: phases + end subroutine elastic_init + module subroutine plastic_init end subroutine plastic_init + module subroutine phase_hooke_SandItsTangents(S,dS_dFe,dS_dFi,Fe,Fi,ph,en) + integer, intent(in) :: & + ph, & + en + real(pReal), intent(in), dimension(3,3) :: & + Fe, & !< elastic deformation gradient + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + S !< 2nd Piola-Kirchhoff stress tensor in lattice configuration + real(pReal), intent(out), dimension(3,3,3,3) :: & + dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient + dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient + end subroutine phase_hooke_SandItsTangents + module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en) real(pReal), dimension(3,3), intent(out) :: & Li !< inleastic velocity gradient @@ -73,7 +82,6 @@ submodule(phase) mechanical end subroutine plastic_isotropic_LiAndItsTangent module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken) - integer, intent(in) :: & co, & !< component-ID of integration point ip, & !< integration point @@ -198,17 +206,11 @@ module subroutine mechanical_init(materials,phases) constituents, & constituent, & phase, & - mech, & - elastic, & - stiffDegradation + mech print'(/,a)', ' <<<+- phase:mechanical init -+>>>' !------------------------------------------------------------------------------------------------- -! initialize elasticity (hooke) !ToDO: Maybe move to elastic submodule along with function homogenizedC? - allocate(phase_elasticity(phases%length), source = ELASTICITY_undefined_ID) - allocate(phase_elasticityInstance(phases%length), source = 0) - allocate(phase_NstiffnessDegradations(phases%length),source=0) allocate(output_constituent(phases%length)) allocate(phase_mechanical_Fe(phases%length)) @@ -253,32 +255,8 @@ module subroutine mechanical_init(materials,phases) #else output_constituent(ph)%label = mech%get_as1dString('output',defaultVal=emptyStringArray) #endif - elastic => mech%get('elastic') - if (IO_lc(elastic%get_asString('type')) == 'hooke') then ! accept small letter h for the moment - phase_elasticity(ph) = ELASTICITY_HOOKE_ID - else - call IO_error(200,ext_msg=elastic%get_asString('type')) - endif - stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList) ! check for stiffness degradation mechanisms - phase_NstiffnessDegradations(ph) = stiffDegradation%length enddo - allocate(phase_stiffnessDegradation(maxval(phase_NstiffnessDegradations),phases%length), & - source=STIFFNESS_DEGRADATION_undefined_ID) - - if(maxVal(phase_NstiffnessDegradations)/=0) then - do ph = 1, phases%length - phase => phases%get(ph) - mech => phase%get('mechanical') - stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList) - do stiffDegradationCtr = 1, stiffDegradation%length - if(stiffDegradation%get_asString(stiffDegradationCtr) == 'damage') & - phase_stiffnessDegradation(stiffDegradationCtr,ph) = STIFFNESS_DEGRADATION_damage_ID - enddo - enddo - endif - - do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2) do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) material => materials%get(discretization_materialAt(el)) @@ -306,6 +284,9 @@ module subroutine mechanical_init(materials,phases) enddo; enddo +! initialize elasticity + call elastic_init(phases) + ! initialize plasticity allocate(plasticState(phases%length)) allocate(phase_plasticity(phases%length),source = PLASTICITY_undefined_ID) @@ -313,9 +294,6 @@ module subroutine mechanical_init(materials,phases) call plastic_init() - do ph = 1, phases%length - phase_elasticityInstance(ph) = count(phase_elasticity(1:ph) == phase_elasticity(ph)) - enddo num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) @@ -348,51 +326,6 @@ module subroutine mechanical_init(materials,phases) end subroutine mechanical_init -!-------------------------------------------------------------------------------------------------- -!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to -!> the elastic and intermediate deformation gradients using Hooke's law -!-------------------------------------------------------------------------------------------------- -subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & - Fe, Fi, ph, en) - - integer, intent(in) :: & - ph, & - en - real(pReal), intent(in), dimension(3,3) :: & - Fe, & !< elastic deformation gradient - Fi !< intermediate deformation gradient - real(pReal), intent(out), dimension(3,3) :: & - S !< 2nd Piola-Kirchhoff stress tensor in lattice configuration - real(pReal), intent(out), dimension(3,3,3,3) :: & - dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient - dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient - - real(pReal), dimension(3,3) :: E - real(pReal), dimension(3,3,3,3) :: C - integer :: & - d, & !< counter in degradation loop - i, j - - C = math_66toSym3333(phase_homogenizedC(ph,en)) - - DegradationLoop: do d = 1, phase_NstiffnessDegradations(ph) - degradationType: select case(phase_stiffnessDegradation(d,ph)) - case (STIFFNESS_DEGRADATION_damage_ID) degradationType - C = C * damage_phi(ph,en)**2 - end select degradationType - enddo DegradationLoop - - E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration - S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration - - do i =1, 3;do j=1,3 - dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko - dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn - enddo; enddo - -end subroutine phase_hooke_SandItsTangents - - module subroutine mechanical_results(group,ph) character(len=*), intent(in) :: group @@ -1082,26 +1015,6 @@ module subroutine mechanical_forward() end subroutine mechanical_forward - -!-------------------------------------------------------------------------------------------------- -!> @brief returns the homogenize elasticity matrix -!> ToDo: homogenizedC66 would be more consistent -!-------------------------------------------------------------------------------------------------- -module function phase_homogenizedC(ph,en) result(C) - - real(pReal), dimension(6,6) :: C - integer, intent(in) :: ph, en - - plasticType: select case (phase_plasticity(ph)) - case (PLASTICITY_DISLOTWIN_ID) plasticType - C = plastic_dislotwin_homogenizedC(ph,en) - case default plasticType - C = lattice_C66(1:6,1:6,ph) - end select plasticType - -end function phase_homogenizedC - - !-------------------------------------------------------------------------------------------------- !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 new file mode 100644 index 000000000..9c18bfa7f --- /dev/null +++ b/src/phase_mechanical_elastic.f90 @@ -0,0 +1,135 @@ +submodule(phase:mechanical) elastic + + enum, bind(c); enumerator :: & + ELASTICITY_UNDEFINED_ID, & + ELASTICITY_HOOKE_ID, & + STIFFNESS_DEGRADATION_UNDEFINED_ID, & + STIFFNESS_DEGRADATION_DAMAGE_ID + end enum + + integer, dimension(:), allocatable :: & + phase_NstiffnessDegradations + + integer(kind(ELASTICITY_UNDEFINED_ID)), dimension(:), allocatable :: & + phase_elasticity !< elasticity of each phase + integer(kind(STIFFNESS_DEGRADATION_UNDEFINED_ID)), dimension(:,:), allocatable :: & + phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase + +contains + + +module subroutine elastic_init(phases) + + class(tNode), pointer :: & + phases + + integer :: & + ph, & + stiffDegradationCtr + class(tNode), pointer :: & + phase, & + mech, & + elastic, & + stiffDegradation + + print'(/,a)', ' <<<+- phase:mechanical:elastic init -+>>>' + + allocate(phase_elasticity(phases%length), source = ELASTICITY_undefined_ID) + allocate(phase_NstiffnessDegradations(phases%length),source=0) + + do ph = 1, phases%length + phase => phases%get(ph) + mech => phase%get('mechanical') + elastic => mech%get('elastic') + if(IO_lc(elastic%get_asString('type')) == 'hooke') then ! accept small letter h for the moment + phase_elasticity(ph) = ELASTICITY_HOOKE_ID + else + call IO_error(200,ext_msg=elastic%get_asString('type')) + endif + stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList) ! check for stiffness degradation mechanisms + phase_NstiffnessDegradations(ph) = stiffDegradation%length + enddo + + allocate(phase_stiffnessDegradation(maxval(phase_NstiffnessDegradations),phases%length), & + source=STIFFNESS_DEGRADATION_undefined_ID) + + if(maxVal(phase_NstiffnessDegradations)/=0) then + do ph = 1, phases%length + phase => phases%get(ph) + mech => phase%get('mechanical') + stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList) + do stiffDegradationCtr = 1, stiffDegradation%length + if(stiffDegradation%get_asString(stiffDegradationCtr) == 'damage') & + phase_stiffnessDegradation(stiffDegradationCtr,ph) = STIFFNESS_DEGRADATION_damage_ID + enddo + enddo + endif + +end subroutine elastic_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to +!> the elastic and intermediate deformation gradients using Hooke's law +!-------------------------------------------------------------------------------------------------- +module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & + Fe, Fi, ph, en) + + integer, intent(in) :: & + ph, & + en + real(pReal), intent(in), dimension(3,3) :: & + Fe, & !< elastic deformation gradient + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + S !< 2nd Piola-Kirchhoff stress tensor in lattice configuration + real(pReal), intent(out), dimension(3,3,3,3) :: & + dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient + dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient + + real(pReal), dimension(3,3) :: E + real(pReal), dimension(3,3,3,3) :: C + integer :: & + d, & !< counter in degradation loop + i, j + + C = math_66toSym3333(phase_homogenizedC(ph,en)) + + DegradationLoop: do d = 1, phase_NstiffnessDegradations(ph) + degradationType: select case(phase_stiffnessDegradation(d,ph)) + case (STIFFNESS_DEGRADATION_damage_ID) degradationType + C = C * damage_phi(ph,en)**2 + end select degradationType + enddo DegradationLoop + + E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration + S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration + + do i =1, 3;do j=1,3 + dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko + dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn + enddo; enddo + +end subroutine phase_hooke_SandItsTangents + + +!-------------------------------------------------------------------------------------------------- +!> @brief returns the homogenized elasticity matrix +!> ToDo: homogenizedC66 would be more consistent +!-------------------------------------------------------------------------------------------------- +module function phase_homogenizedC(ph,en) result(C) + + real(pReal), dimension(6,6) :: C + integer, intent(in) :: ph, en + + plasticType: select case (phase_plasticity(ph)) + case (PLASTICITY_DISLOTWIN_ID) plasticType + C = plastic_dislotwin_homogenizedC(ph,en) + case default plasticType + C = lattice_C66(1:6,1:6,ph) + end select plasticType + +end function phase_homogenizedC + + +end submodule elastic