223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch 'more-HDF5-postprocessing-cleaned' into development

This commit is contained in:
Vitesh Shah 2019-05-27 12:17:12 +02:00
parent e743a55ff2 b35465b591
commit 5abcca50ee
11 changed files with 231 additions and 418 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit d96bfb32920c96a8a43958f76a209d34c6bd841a Subproject commit 3a2f89547c264044a7bfab9d33aee78eec495a76

37
processing/post/DADF5_postResults.py
View File

@ -22,10 +22,15 @@ parser.add_argument('filenames', nargs='+',
help='DADF5 files') help='DADF5 files')
parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string', parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string',
help='name of subdirectory to hold output') help='name of subdirectory to hold output')
parser.add_argument('--mat', nargs='+',
help='labels for materialpoint/homogenization',dest='mat')
parser.add_argument('--con', nargs='+',
help='labels for constituent/crystallite/constitutive',dest='con')
options = parser.parse_args() options = parser.parse_args()
options.labels = ['Fe','Fp','xi_sl'] if options.mat is None: options.mat=[]
if options.con is None: options.con=[]
# --- loop over input files ------------------------------------------------------------------------ # --- loop over input files ------------------------------------------------------------------------
@ -49,15 +54,12 @@ for filename in options.filenames:
data = np.array([inc['inc'] for j in range(np.product(results.grid))]).reshape([np.product(results.grid),1]) data = np.array([inc['inc'] for j in range(np.product(results.grid))]).reshape([np.product(results.grid),1])
header+= 'inc' header+= 'inc'
data = np.concatenate((data,np.array([j+1 for j in range(np.product(results.grid))]).reshape([np.product(results.grid),1])),1)
header+=' node'
coords = coords.reshape([np.product(results.grid),3]) coords = coords.reshape([np.product(results.grid),3])
data = np.concatenate((data,coords),1) data = np.concatenate((data,coords),1)
header+=' 1_pos 2_pos 3_pos' header+=' 1_pos 2_pos 3_pos'
results.active['increments'] = [inc] results.active['increments'] = [inc]
for label in options.labels: for label in options.con:
for o in results.c_output_types: for o in results.c_output_types:
results.active['c_output_types'] = [o] results.active['c_output_types'] = [o]
for c in results.constituents: for c in results.constituents:
@ -67,12 +69,33 @@ for filename in options.filenames:
continue continue
label = x[0].split('/')[-1] label = x[0].split('/')[-1]
array = results.read_dataset(x,0) array = results.read_dataset(x,0)
d = np.product(np.shape(array)[1:]) d = int(np.product(np.shape(array)[1:]))
array = np.reshape(array,[np.product(results.grid),d]) array = np.reshape(array,[np.product(results.grid),d])
data = np.concatenate((data,array),1) data = np.concatenate((data,array),1)
header+= ''.join([' {}_{}'.format(j+1,label) for j in range(d)]) if d>1:
header+= ''.join([' {}_{}'.format(j+1,label) for j in range(d)])
else:
header+=' '+label
for label in options.mat:
for o in results.m_output_types:
results.active['m_output_types'] = [o]
for m in results.materialpoints:
results.active['materialpoints'] = [m]
x = results.get_dataset_location(label)
if len(x) == 0:
continue
label = x[0].split('/')[-1]
array = results.read_dataset(x,0)
d = int(np.product(np.shape(array)[1:]))
array = np.reshape(array,[np.product(results.grid),d])
data = np.concatenate((data,array),1)
if d>1:
header+= ''.join([' {}_{}'.format(j+1,label) for j in range(d)])
else:
header+=' '+label
dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir)) dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir))
try: try:

11
processing/post/DADF5_vtk_cells.py
View File

@ -23,10 +23,15 @@ parser.add_argument('filenames', nargs='+',
help='DADF5 files') help='DADF5 files')
parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string', parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string',
help='name of subdirectory to hold output') help='name of subdirectory to hold output')
parser.add_argument('--mat', nargs='+',
help='labels for materialpoint/homogenization',dest='mat')
parser.add_argument('--con', nargs='+',
help='labels for constituent/crystallite/constitutive',dest='con')
options = parser.parse_args() options = parser.parse_args()
options.labels = ['Fe','Fp','xi_sl'] if options.mat is None: options.mat=[]
if options.con is None: options.con=[]
# --- loop over input files ------------------------------------------------------------------------ # --- loop over input files ------------------------------------------------------------------------
@ -54,7 +59,9 @@ for filename in options.filenames:
print('Output step {}/{}'.format(i+1,len(results.increments))) print('Output step {}/{}'.format(i+1,len(results.increments)))
vtk_data = [] vtk_data = []
results.active['increments'] = [inc] results.active['increments'] = [inc]
for label in options.labels:
for label in options.con:
for o in results.c_output_types: for o in results.c_output_types:
results.active['c_output_types'] = [o] results.active['c_output_types'] = [o]
if o != 'generic': if o != 'generic':

57
python/damask/dadf5.py
View File

@ -49,15 +49,23 @@ class DADF5():
self.c_output_types.append(o) self.c_output_types.append(o)
self.c_output_types = list(set(self.c_output_types)) # make unique self.c_output_types = list(set(self.c_output_types)) # make unique
self.m_output_types = []
for m in self.materialpoints:
for o in f['inc{:05}/materialpoint/{}'.format(self.increments[0]['inc'],m)].keys():
self.m_output_types.append(o)
self.m_output_types = list(set(self.m_output_types)) # make unique
self.active= {'increments': self.increments, self.active= {'increments': self.increments,
'constituents': self.constituents, 'constituents': self.constituents,
'materialpoints': self.materialpoints, 'materialpoints': self.materialpoints,
'constituent': self.Nconstituents, 'constituent': self.Nconstituents,
'c_output_types': self.c_output_types} 'c_output_types': self.c_output_types,
'm_output_types': self.m_output_types}
self.filename = filename self.filename = filename
self.mode = mode self.mode = mode
def list_data(self): def list_data(self):
"""Shows information on all datasets in the file""" """Shows information on all datasets in the file"""
with h5py.File(self.filename,'r') as f: with h5py.File(self.filename,'r') as f:
@ -73,6 +81,16 @@ class DADF5():
print(' {} ({})'.format(x,f[group_output_types+'/'+x].attrs['Description'].decode())) print(' {} ({})'.format(x,f[group_output_types+'/'+x].attrs['Description'].decode()))
except: except:
pass pass
for m in self.active['materialpoints']:
group_materialpoint = group_inc+'/materialpoint/'+m
for t in self.active['m_output_types']:
print(' {}'.format(t))
group_output_types = group_materialpoint+'/'+t
try:
for x in f[group_output_types].keys():
print(' {} ({})'.format(x,f[group_output_types+'/'+x].attrs['Description'].decode()))
except:
pass
def get_dataset_location(self,label): def get_dataset_location(self,label):
@ -81,14 +99,25 @@ class DADF5():
with h5py.File(self.filename,'r') as f: with h5py.File(self.filename,'r') as f:
for i in self.active['increments']: for i in self.active['increments']:
group_inc = 'inc{:05}'.format(i['inc']) group_inc = 'inc{:05}'.format(i['inc'])
for c in self.active['constituents']: for c in self.active['constituents']:
group_constituent = group_inc+'/constituent/'+c group_constituent = group_inc+'/constituent/'+c
for t in self.active['c_output_types']: for t in self.active['c_output_types']:
try: try:
f[group_constituent+'/'+t+'/'+label] f[group_constituent+'/'+t+'/'+label]
path.append(group_constituent+'/'+t+'/'+label) path.append(group_constituent+'/'+t+'/'+label)
except: except Exception as e:
pass print('unable to locate constituents dataset: '+ str(e))
for m in self.active['materialpoints']:
group_materialpoint = group_inc+'/materialpoint/'+m
for t in self.active['m_output_types']:
try:
f[group_materialpoint+'/'+t+'/'+label]
path.append(group_materialpoint+'/'+t+'/'+label)
except Exception as e:
print('unable to locate materialpoints dataset: '+ str(e))
return path return path
@ -100,12 +129,28 @@ class DADF5():
""" """
with h5py.File(self.filename,'r') as f: with h5py.File(self.filename,'r') as f:
shape = (self.Nmaterialpoints,) + np.shape(f[path[0]])[1:] shape = (self.Nmaterialpoints,) + np.shape(f[path[0]])[1:]
if len(shape) == 1: shape = shape +(1,)
dataset = np.full(shape,np.nan) dataset = np.full(shape,np.nan)
for pa in path: for pa in path:
label = pa.split('/')[2] label = pa.split('/')[2]
p = np.where(f['mapping/cellResults/constituent'][:,c]['Name'] == str.encode(label))[0] try:
u = (f['mapping/cellResults/constituent'][p,c]['Position']) p = np.where(f['mapping/cellResults/constituent'][:,c]['Name'] == str.encode(label))[0]
dataset[p,:] = f[pa][u,:] u = (f['mapping/cellResults/constituent'][p,c]['Position'])
a = np.array(f[pa])
if len(a.shape) == 1:
a=a.reshape([a.shape[0],1])
dataset[p,:] = a[u,:]
except Exception as e:
print('unable to read constituent: '+ str(e))
try:
p = np.where(f['mapping/cellResults/materialpoint']['Name'] == str.encode(label))[0]
u = (f['mapping/cellResults/materialpoint'][p.tolist()]['Position'])
a = np.array(f[pa])
if len(a.shape) == 1:
a=a.reshape([a.shape[0],1])
dataset[p,:] = a[u,:]
except Exception as e:
print('unable to read materialpoint: '+ str(e))
return dataset return dataset

2
src/commercialFEM_fileList.f90
View File

@ -46,7 +46,6 @@
#include "plastic_nonlocal.f90" #include "plastic_nonlocal.f90"
#include "constitutive.f90" #include "constitutive.f90"
#include "crystallite.f90" #include "crystallite.f90"
#include "homogenization_mech_RGC.f90"
#include "thermal_isothermal.f90" #include "thermal_isothermal.f90"
#include "thermal_adiabatic.f90" #include "thermal_adiabatic.f90"
#include "thermal_conduction.f90" #include "thermal_conduction.f90"
@ -56,4 +55,5 @@
#include "homogenization.f90" #include "homogenization.f90"
#include "homogenization_mech_none.f90" #include "homogenization_mech_none.f90"
#include "homogenization_mech_isostrain.f90" #include "homogenization_mech_isostrain.f90"
#include "homogenization_mech_RGC.f90"
#include "CPFEM.f90" #include "CPFEM.f90"

255
src/homogenization.f90
View File

@ -6,7 +6,20 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module homogenization module homogenization
use prec use prec
use IO
use config
use debug
use math
use material use material
use numerics
use constitutive
use crystallite
use mesh
use FEsolving
#if defined(PETSc) || defined(DAMASK_HDF5)
use results
use HDF5_utilities
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! General variables for the homogenization at a material point ! General variables for the homogenization at a material point
@ -22,7 +35,6 @@ module homogenization
materialpoint_results !< results array of material point materialpoint_results !< results array of material point
integer, public, protected :: & integer, public, protected :: &
materialpoint_sizeResults, & materialpoint_sizeResults, &
homogenization_maxSizePostResults, &
thermal_maxSizePostResults, & thermal_maxSizePostResults, &
damage_maxSizePostResults damage_maxSizePostResults
@ -47,11 +59,24 @@ module homogenization
module subroutine mech_isostrain_init module subroutine mech_isostrain_init
end subroutine mech_isostrain_init end subroutine mech_isostrain_init
module subroutine mech_RGC_init
end subroutine mech_RGC_init
module subroutine mech_isostrain_partitionDeformation(F,avgF) module subroutine mech_isostrain_partitionDeformation(F,avgF)
real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient
real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point
end subroutine mech_isostrain_partitionDeformation end subroutine mech_isostrain_partitionDeformation
module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient
real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point
integer, intent(in) :: &
instance, &
of
end subroutine mech_RGC_partitionDeformation
module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
@ -61,6 +86,36 @@ module homogenization
integer, intent(in) :: instance integer, intent(in) :: instance
end subroutine mech_isostrain_averageStressAndItsTangent end subroutine mech_isostrain_averageStressAndItsTangent
module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses
real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
integer, intent(in) :: instance
end subroutine mech_RGC_averageStressAndItsTangent
module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
logical, dimension(2) :: mech_RGC_updateState
real(pReal), dimension(:,:,:), intent(in) :: &
P,& !< partitioned stresses
F,& !< partitioned deformation gradients
F0 !< partitioned initial deformation gradients
real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
real(pReal), dimension(3,3), intent(in) :: avgF !< average F
real(pReal), intent(in) :: dt !< time increment
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
end function mech_RGC_updateState
module subroutine mech_RGC_results(instance,group)
integer, intent(in) :: instance !< homogenization instance
character(len=*), intent(in) :: group !< group name in HDF5 file
end subroutine mech_RGC_results
end interface end interface
public :: & public :: &
@ -68,11 +123,6 @@ module homogenization
materialpoint_stressAndItsTangent, & materialpoint_stressAndItsTangent, &
materialpoint_postResults, & materialpoint_postResults, &
homogenization_results homogenization_results
private :: &
partitionDeformation, &
updateState, &
averageStressAndItsTangent, &
postResults
contains contains
@ -81,36 +131,12 @@ contains
!> @brief module initialization !> @brief module initialization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_init subroutine homogenization_init
use math, only: &
math_I3
use debug, only: &
debug_level, &
debug_homogenization, &
debug_levelBasic, &
debug_e, &
debug_g
use mesh, only: &
theMesh, &
mesh_element
use constitutive, only: &
constitutive_plasticity_maxSizePostResults, &
constitutive_source_maxSizePostResults
use crystallite, only: &
crystallite_maxSizePostResults
use config, only: &
config_deallocate, &
config_homogenization, &
homogenization_name
use homogenization_mech_RGC
use thermal_isothermal use thermal_isothermal
use thermal_adiabatic use thermal_adiabatic
use thermal_conduction use thermal_conduction
use damage_none use damage_none
use damage_local use damage_local
use damage_nonlocal use damage_nonlocal
use IO
use numerics, only: &
worldrank
integer, parameter :: FILEUNIT = 200 integer, parameter :: FILEUNIT = 200
integer :: e,i,p integer :: e,i,p
@ -120,10 +146,9 @@ subroutine homogenization_init
character(len=32) :: outputName !< name of output, intermediate fix until HDF5 output is ready character(len=32) :: outputName !< name of output, intermediate fix until HDF5 output is ready
logical :: valid logical :: valid
if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mech_none_init if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mech_none_init
if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mech_isostrain_init if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mech_isostrain_init
if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call homogenization_RGC_init if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call mech_RGC_init
if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init
if (any(thermal_type == THERMAL_adiabatic_ID)) call thermal_adiabatic_init if (any(thermal_type == THERMAL_adiabatic_ID)) call thermal_adiabatic_init
@ -135,39 +160,14 @@ subroutine homogenization_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! write description file for homogenization output ! write description file for homogenization output
mainProcess2: if (worldrank == 0) then mainProcess: if (worldrank == 0) then
call IO_write_jobFile(FILEUNIT,'outputHomogenization') call IO_write_jobFile(FILEUNIT,'outputHomogenization')
do p = 1,size(config_homogenization) do p = 1,size(config_homogenization)
if (any(material_homogenizationAt == p)) then if (any(material_homogenizationAt == p)) then
i = homogenization_typeInstance(p) ! which instance of this homogenization type
valid = .true. ! assume valid
select case(homogenization_type(p)) ! split per homogenization type
case (HOMOGENIZATION_NONE_ID)
outputName = HOMOGENIZATION_NONE_label
thisOutput => null()
thisSize => null()
case (HOMOGENIZATION_ISOSTRAIN_ID)
outputName = HOMOGENIZATION_ISOSTRAIN_label
thisOutput => null()
thisSize => null()
case (HOMOGENIZATION_RGC_ID)
outputName = HOMOGENIZATION_RGC_label
thisOutput => homogenization_RGC_output
thisSize => homogenization_RGC_sizePostResult
case default
valid = .false.
end select
write(FILEUNIT,'(/,a,/)') '['//trim(homogenization_name(p))//']' write(FILEUNIT,'(/,a,/)') '['//trim(homogenization_name(p))//']'
if (valid) then write(FILEUNIT,'(a)') '(type) n/a'
write(FILEUNIT,'(a)') '(type)'//char(9)//trim(outputName) write(FILEUNIT,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p)
write(FILEUNIT,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p)
if (homogenization_type(p) /= HOMOGENIZATION_NONE_ID .and. &
homogenization_type(p) /= HOMOGENIZATION_ISOSTRAIN_ID) then
do e = 1,size(thisOutput(:,i))
write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
enddo
endif
endif
i = thermal_typeInstance(p) ! which instance of this thermal type i = thermal_typeInstance(p) ! which instance of this thermal type
valid = .true. ! assume valid valid = .true. ! assume valid
select case(thermal_type(p)) ! split per thermal type select case(thermal_type(p)) ! split per thermal type
@ -197,6 +197,7 @@ subroutine homogenization_init
enddo enddo
endif endif
endif endif
i = damage_typeInstance(p) ! which instance of this damage type i = damage_typeInstance(p) ! which instance of this damage type
valid = .true. ! assume valid valid = .true. ! assume valid
select case(damage_type(p)) ! split per damage type select case(damage_type(p)) ! split per damage type
@ -229,7 +230,7 @@ subroutine homogenization_init
endif endif
enddo enddo
close(FILEUNIT) close(FILEUNIT)
endif mainProcess2 endif mainProcess
call config_deallocate('material.config/homogenization') call config_deallocate('material.config/homogenization')
@ -237,7 +238,7 @@ subroutine homogenization_init
! allocate and initialize global variables ! allocate and initialize global variables
allocate(materialpoint_dPdF(3,3,3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) allocate(materialpoint_dPdF(3,3,3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal)
allocate(materialpoint_F0(3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) allocate(materialpoint_F0(3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal)
materialpoint_F0 = spread(spread(math_I3,3,theMesh%elem%nIPs),4,theMesh%nElems) ! initialize to identity materialpoint_F0 = spread(spread(math_I3,3,theMesh%elem%nIPs),4,theMesh%nElems) ! initialize to identity
allocate(materialpoint_F(3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) allocate(materialpoint_F(3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal)
materialpoint_F = materialpoint_F0 ! initialize to identity materialpoint_F = materialpoint_F0 ! initialize to identity
allocate(materialpoint_subF0(3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) allocate(materialpoint_subF0(3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal)
@ -252,18 +253,15 @@ subroutine homogenization_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate and initialize global state and postresutls variables ! allocate and initialize global state and postresutls variables
homogenization_maxSizePostResults = 0
thermal_maxSizePostResults = 0 thermal_maxSizePostResults = 0
damage_maxSizePostResults = 0 damage_maxSizePostResults = 0
do p = 1,size(config_homogenization) do p = 1,size(config_homogenization)
homogenization_maxSizePostResults = max(homogenization_maxSizePostResults,homogState (p)%sizePostResults)
thermal_maxSizePostResults = max(thermal_maxSizePostResults, thermalState (p)%sizePostResults) thermal_maxSizePostResults = max(thermal_maxSizePostResults, thermalState (p)%sizePostResults)
damage_maxSizePostResults = max(damage_maxSizePostResults ,damageState (p)%sizePostResults) damage_maxSizePostResults = max(damage_maxSizePostResults ,damageState (p)%sizePostResults)
enddo enddo
materialpoint_sizeResults = 1 & ! grain count materialpoint_sizeResults = 1 & ! grain count
+ 1 + homogenization_maxSizePostResults & ! homogSize & homogResult + 1 + thermal_maxSizePostResults &
+ thermal_maxSizePostResults &
+ damage_maxSizePostResults & + damage_maxSizePostResults &
+ homogenization_maxNgrains * (1 + crystallite_maxSizePostResults & ! crystallite size & crystallite results + homogenization_maxNgrains * (1 + crystallite_maxSizePostResults & ! crystallite size & crystallite results
+ 1 + constitutive_plasticity_maxSizePostResults & ! constitutive size & constitutive results + 1 + constitutive_plasticity_maxSizePostResults & ! constitutive size & constitutive results
@ -273,11 +271,6 @@ subroutine homogenization_init
write(6,'(/,a)') ' <<<+- homogenization init -+>>>' write(6,'(/,a)') ' <<<+- homogenization init -+>>>'
if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0) then if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0) then
#ifdef TODO
write(6,'(a32,1x,7(i8,1x))') 'homogenization_state0: ', shape(homogenization_state0)
write(6,'(a32,1x,7(i8,1x))') 'homogenization_subState0: ', shape(homogenization_subState0)
write(6,'(a32,1x,7(i8,1x))') 'homogenization_state: ', shape(homogenization_state)
#endif
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F0: ', shape(materialpoint_F0) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F0: ', shape(materialpoint_F0)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F: ', shape(materialpoint_F) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F: ', shape(materialpoint_F)
@ -303,51 +296,6 @@ end subroutine homogenization_init
!> @brief parallelized calculation of stress and corresponding tangent at material points !> @brief parallelized calculation of stress and corresponding tangent at material points
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine materialpoint_stressAndItsTangent(updateJaco,dt) subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
use numerics, only: &
subStepMinHomog, &
subStepSizeHomog, &
stepIncreaseHomog, &
nMPstate
use FEsolving, only: &
FEsolving_execElem, &
FEsolving_execIP, &
terminallyIll
use mesh, only: &
mesh_element
use crystallite, only: &
crystallite_F0, &
crystallite_Fp0, &
crystallite_Fp, &
crystallite_Fi0, &
crystallite_Fi, &
crystallite_Lp0, &
crystallite_Lp, &
crystallite_Li0, &
crystallite_Li, &
crystallite_S0, &
crystallite_S, &
crystallite_partionedF0, &
crystallite_partionedF, &
crystallite_partionedFp0, &
crystallite_partionedLp0, &
crystallite_partionedFi0, &
crystallite_partionedLi0, &
crystallite_partionedS0, &
crystallite_dt, &
crystallite_requested, &
crystallite_stress, &
crystallite_stressTangent, &
crystallite_orientations
#ifdef DEBUG
use debug, only: &
debug_level, &
debug_homogenization, &
debug_levelBasic, &
debug_levelExtensive, &
debug_levelSelective, &
debug_e, &
debug_i
#endif
real(pReal), intent(in) :: dt !< time increment real(pReal), intent(in) :: dt !< time increment
logical, intent(in) :: updateJaco !< initiating Jacobian update logical, intent(in) :: updateJaco !< initiating Jacobian update
@ -640,14 +588,6 @@ end subroutine materialpoint_stressAndItsTangent
!> @brief parallelized calculation of result array at material points !> @brief parallelized calculation of result array at material points
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine materialpoint_postResults subroutine materialpoint_postResults
use FEsolving, only: &
FEsolving_execElem, &
FEsolving_execIP
use mesh, only: &
mesh_element
use crystallite, only: &
crystallite_sizePostResults, &
crystallite_postResults
integer :: & integer :: &
thePos, & thePos, &
@ -697,12 +637,6 @@ end subroutine materialpoint_postResults
!> @brief partition material point def grad onto constituents !> @brief partition material point def grad onto constituents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine partitionDeformation(ip,el) subroutine partitionDeformation(ip,el)
use mesh, only: &
mesh_element
use crystallite, only: &
crystallite_partionedF
use homogenization_mech_RGC, only: &
homogenization_RGC_partitionDeformation
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point ip, & !< integration point
@ -719,7 +653,7 @@ subroutine partitionDeformation(ip,el)
materialpoint_subF(1:3,1:3,ip,el)) materialpoint_subF(1:3,1:3,ip,el))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
call homogenization_RGC_partitionDeformation(& call mech_RGC_partitionDeformation(&
crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
materialpoint_subF(1:3,1:3,ip,el),& materialpoint_subF(1:3,1:3,ip,el),&
ip, & ip, &
@ -734,15 +668,6 @@ end subroutine partitionDeformation
!> "happy" with result !> "happy" with result
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function updateState(ip,el) function updateState(ip,el)
use mesh, only: &
mesh_element
use crystallite, only: &
crystallite_P, &
crystallite_dPdF, &
crystallite_partionedF,&
crystallite_partionedF0
use homogenization_mech_RGC, only: &
homogenization_RGC_updateState
use thermal_adiabatic, only: & use thermal_adiabatic, only: &
thermal_adiabatic_updateState thermal_adiabatic_updateState
use damage_local, only: & use damage_local, only: &
@ -758,14 +683,14 @@ function updateState(ip,el)
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
updateState = & updateState = &
updateState .and. & updateState .and. &
homogenization_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
crystallite_partionedF0(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el),& crystallite_partionedF0(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el),&
materialpoint_subF(1:3,1:3,ip,el),& materialpoint_subF(1:3,1:3,ip,el),&
materialpoint_subdt(ip,el), & materialpoint_subdt(ip,el), &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
ip, & ip, &
el) el)
end select chosenHomogenization end select chosenHomogenization
chosenThermal: select case (thermal_type(mesh_element(3,el))) chosenThermal: select case (thermal_type(mesh_element(3,el)))
@ -793,12 +718,6 @@ end function updateState
!> @brief derive average stress and stiffness from constituent quantities !> @brief derive average stress and stiffness from constituent quantities
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine averageStressAndItsTangent(ip,el) subroutine averageStressAndItsTangent(ip,el)
use mesh, only: &
mesh_element
use crystallite, only: &
crystallite_P,crystallite_dPdF
use homogenization_mech_RGC, only: &
homogenization_RGC_averageStressAndItsTangent
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point ip, & !< integration point
@ -818,7 +737,7 @@ subroutine averageStressAndItsTangent(ip,el)
homogenization_typeInstance(mesh_element(3,el))) homogenization_typeInstance(mesh_element(3,el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
call homogenization_RGC_averageStressAndItsTangent(& call mech_RGC_averageStressAndItsTangent(&
materialpoint_P(1:3,1:3,ip,el), & materialpoint_P(1:3,1:3,ip,el), &
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),&
crystallite_P(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & crystallite_P(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
@ -834,10 +753,6 @@ end subroutine averageStressAndItsTangent
!> if homogenization_sizePostResults(i,e) > 0 !! !> if homogenization_sizePostResults(i,e) > 0 !!
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function postResults(ip,el) function postResults(ip,el)
use mesh, only: &
mesh_element
use homogenization_mech_RGC, only: &
homogenization_RGC_postResults
use thermal_adiabatic, only: & use thermal_adiabatic, only: &
thermal_adiabatic_postResults thermal_adiabatic_postResults
use thermal_conduction, only: & use thermal_conduction, only: &
@ -856,23 +771,12 @@ function postResults(ip,el)
postResults postResults
integer :: & integer :: &
startPos, endPos ,& startPos, endPos ,&
of, instance, homog homog
postResults = 0.0_pReal postResults = 0.0_pReal
startPos = 1 startPos = 1
endPos = homogState(material_homogenizationAt(el))%sizePostResults endPos = thermalState(material_homogenizationAt(el))%sizePostResults
chosenHomogenization: select case (homogenization_type(mesh_element(3,el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
instance = homogenization_typeInstance(material_homogenizationAt(el))
of = mappingHomogenization(1,ip,el)
postResults(startPos:endPos) = homogenization_RGC_postResults(instance,of)
end select chosenHomogenization
startPos = endPos + 1
endPos = endPos + thermalState(material_homogenizationAt(el))%sizePostResults
chosenThermal: select case (thermal_type(mesh_element(3,el))) chosenThermal: select case (thermal_type(mesh_element(3,el)))
case (THERMAL_adiabatic_ID) chosenThermal case (THERMAL_adiabatic_ID) chosenThermal
@ -905,14 +809,9 @@ end function postResults
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_results subroutine homogenization_results
#if defined(PETSc) || defined(DAMASK_HDF5) #if defined(PETSc) || defined(DAMASK_HDF5)
use results
use homogenization_mech_RGC
use HDF5_utilities
use config, only: & use config, only: &
config_name_homogenization => homogenization_name ! anticipate logical name config_name_homogenization => homogenization_name ! anticipate logical name
use material, only: & use material, only: &
homogenization_typeInstance, &
material_homogenization_type => homogenization_type material_homogenization_type => homogenization_type
integer :: p integer :: p

235
src/homogenization_mech_RGC.f90
View File

@ -6,16 +6,7 @@
!> @brief Relaxed grain cluster (RGC) homogenization scheme !> @brief Relaxed grain cluster (RGC) homogenization scheme
!> Nconstituents is defined as p x q x r (cluster) !> Nconstituents is defined as p x q x r (cluster)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module homogenization_mech_RGC submodule(homogenization) homogenization_mech_RGC
use prec
use material
implicit none
private
integer, dimension(:,:), allocatable,target, public :: &
homogenization_RGC_sizePostResult
character(len=64), dimension(:,:), allocatable,target, public :: &
homogenization_RGC_output ! name of each post result output
enum, bind(c) enum, bind(c)
enumerator :: & enumerator :: &
@ -28,7 +19,7 @@ module homogenization_mech_RGC
magnitudemismatch_ID magnitudemismatch_ID
end enum end enum
type, private :: tParameters type :: tParameters
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
Nconstituents Nconstituents
real(pReal) :: & real(pReal) :: &
@ -39,11 +30,11 @@ module homogenization_mech_RGC
angles angles
integer :: & integer :: &
of_debug = 0 of_debug = 0
integer(kind(undefined_ID)), dimension(:), allocatable :: & integer(kind(undefined_ID)), dimension(:), allocatable :: &
outputID outputID
end type tParameters end type tParameters
type, private :: tRGCstate type :: tRGCstate
real(pReal), pointer, dimension(:) :: & real(pReal), pointer, dimension(:) :: &
work, & work, &
penaltyEnergy penaltyEnergy
@ -51,7 +42,7 @@ module homogenization_mech_RGC
relaxationVector relaxationVector
end type tRGCstate end type tRGCstate
type, private :: tRGCdependentState type :: tRGCdependentState
real(pReal), allocatable, dimension(:) :: & real(pReal), allocatable, dimension(:) :: &
volumeDiscrepancy, & volumeDiscrepancy, &
relaxationRate_avg, & relaxationRate_avg, &
@ -62,56 +53,25 @@ module homogenization_mech_RGC
orientation orientation
end type tRGCdependentState end type tRGCdependentState
type(tparameters), dimension(:), allocatable, private :: & type(tparameters), dimension(:), allocatable :: &
param param
type(tRGCstate), dimension(:), allocatable, private :: & type(tRGCstate), dimension(:), allocatable :: &
state, & state, &
state0 state0
type(tRGCdependentState), dimension(:), allocatable, private :: & type(tRGCdependentState), dimension(:), allocatable :: &
dependentState dependentState
public :: &
homogenization_RGC_init, &
homogenization_RGC_partitionDeformation, &
homogenization_RGC_averageStressAndItsTangent, &
homogenization_RGC_updateState, &
homogenization_RGC_postResults, &
mech_RGC_results ! name suited for planned submodule situation
private :: &
relaxationVector, &
interfaceNormal, &
getInterface, &
grain1to3, &
grain3to1, &
interface4to1, &
interface1to4
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields, reads information from material configuration file !> @brief allocates all necessary fields, reads information from material configuration file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_RGC_init() module subroutine mech_RGC_init
use debug, only: &
#ifdef DEBUG
debug_i, &
debug_e, &
#endif
debug_level, &
debug_homogenization, &
debug_levelBasic
use math, only: &
math_EulerToR, &
INRAD
use IO, only: &
IO_error
use config, only: &
config_homogenization
integer :: & integer :: &
Ninstance, & Ninstance, &
h, i, & h, i, &
NofMyHomog, outputSize, & NofMyHomog, &
sizeState, nIntFaceTot sizeState, nIntFaceTot
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
@ -139,9 +99,6 @@ subroutine homogenization_RGC_init()
allocate(state0(Ninstance)) allocate(state0(Ninstance))
allocate(dependentState(Ninstance)) allocate(dependentState(Ninstance))
allocate(homogenization_RGC_sizePostResult(maxval(homogenization_Noutput),Ninstance),source=0)
allocate(homogenization_RGC_output(maxval(homogenization_Noutput),Ninstance))
homogenization_RGC_output=''
do h = 1, size(homogenization_type) do h = 1, size(homogenization_type)
if (homogenization_type(h) /= HOMOGENIZATION_RGC_ID) cycle if (homogenization_type(h) /= HOMOGENIZATION_RGC_ID) cycle
@ -176,28 +133,20 @@ subroutine homogenization_RGC_init()
case('constitutivework') case('constitutivework')
outputID = constitutivework_ID outputID = constitutivework_ID
outputSize = 1
case('penaltyenergy') case('penaltyenergy')
outputID = penaltyenergy_ID outputID = penaltyenergy_ID
outputSize = 1
case('volumediscrepancy') case('volumediscrepancy')
outputID = volumediscrepancy_ID outputID = volumediscrepancy_ID
outputSize = 1
case('averagerelaxrate') case('averagerelaxrate')
outputID = averagerelaxrate_ID outputID = averagerelaxrate_ID
outputSize = 1
case('maximumrelaxrate') case('maximumrelaxrate')
outputID = maximumrelaxrate_ID outputID = maximumrelaxrate_ID
outputSize = 1
case('magnitudemismatch') case('magnitudemismatch')
outputID = magnitudemismatch_ID outputID = magnitudemismatch_ID
outputSize = 3
end select end select
if (outputID /= undefined_ID) then if (outputID /= undefined_ID) then
homogenization_RGC_output(i,homogenization_typeInstance(h)) = outputs(i)
homogenization_RGC_sizePostResult(i,homogenization_typeInstance(h)) = outputSize
prm%outputID = [prm%outputID , outputID] prm%outputID = [prm%outputID , outputID]
endif endif
@ -211,7 +160,7 @@ subroutine homogenization_RGC_init()
+ size(['avg constitutive work ','average penalty energy']) + size(['avg constitutive work ','average penalty energy'])
homogState(h)%sizeState = sizeState homogState(h)%sizeState = sizeState
homogState(h)%sizePostResults = sum(homogenization_RGC_sizePostResult(:,homogenization_typeInstance(h))) homogState(h)%sizePostResults = 0
allocate(homogState(h)%state0 (sizeState,NofMyHomog), source=0.0_pReal) allocate(homogState(h)%state0 (sizeState,NofMyHomog), source=0.0_pReal)
allocate(homogState(h)%subState0(sizeState,NofMyHomog), source=0.0_pReal) allocate(homogState(h)%subState0(sizeState,NofMyHomog), source=0.0_pReal)
allocate(homogState(h)%state (sizeState,NofMyHomog), source=0.0_pReal) allocate(homogState(h)%state (sizeState,NofMyHomog), source=0.0_pReal)
@ -235,23 +184,17 @@ subroutine homogenization_RGC_init()
enddo enddo
end subroutine homogenization_RGC_init end subroutine mech_RGC_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief partitions the deformation gradient onto the constituents !> @brief partitions the deformation gradient onto the constituents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_RGC_partitionDeformation(F,avgF,instance,of) module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
#ifdef DEBUG
use debug, only: &
debug_level, &
debug_homogenization, &
debug_levelExtensive
#endif
real(pReal), dimension (:,:,:), intent(out) :: F !< partioned F per grain real(pReal), dimension (:,:,:), intent(out) :: F !< partioned F per grain
real(pReal), dimension (:,:), intent(in) :: avgF !< averaged F real(pReal), dimension (3,3), intent(in) :: avgF !< averaged F
integer, intent(in) :: & integer, intent(in) :: &
instance, & instance, &
of of
@ -291,45 +234,14 @@ subroutine homogenization_RGC_partitionDeformation(F,avgF,instance,of)
end associate end associate
end subroutine homogenization_RGC_partitionDeformation end subroutine mech_RGC_partitionDeformation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief update the internal state of the homogenization scheme and tell whether "done" and !> @brief update the internal state of the homogenization scheme and tell whether "done" and
! "happy" with result ! "happy" with result
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) module procedure mech_RGC_updateState
#ifdef DEBUG
use debug, only: &
debug_level, &
debug_homogenization,&
debug_levelExtensive
#endif
use math, only: &
math_invert2
use numerics, only: &
absTol_RGC, &
relTol_RGC, &
absMax_RGC, &
relMax_RGC, &
pPert_RGC, &
maxdRelax_RGC, &
viscPower_RGC, &
viscModus_RGC, &
refRelaxRate_RGC
real(pReal), dimension(:,:,:), intent(in) :: &
P,& !< array of P
F,& !< array of F
F0 !< array of initial F
real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< array of current grain stiffness
real(pReal), dimension(3,3), intent(in) :: avgF !< average F
real(pReal), intent(in) :: dt !< time increment
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
logical, dimension(2) :: homogenization_RGC_updateState
integer, dimension(4) :: intFaceN,intFaceP,faceID integer, dimension(4) :: intFaceN,intFaceP,faceID
integer, dimension(3) :: nGDim,iGr3N,iGr3P integer, dimension(3) :: nGDim,iGr3N,iGr3P
@ -347,7 +259,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
#endif #endif
zeroTimeStep: if(dEq0(dt)) then zeroTimeStep: if(dEq0(dt)) then
homogenization_RGC_updateState = .true. ! pretend everything is fine and return mech_RGC_updateState = .true. ! pretend everything is fine and return
return return
endif zeroTimeStep endif zeroTimeStep
@ -466,12 +378,12 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
endif endif
#endif #endif
homogenization_RGC_updateState = .false. mech_RGC_updateState = .false.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! If convergence reached => done and happy ! If convergence reached => done and happy
if (residMax < relTol_RGC*stresMax .or. residMax < absTol_RGC) then if (residMax < relTol_RGC*stresMax .or. residMax < absTol_RGC) then
homogenization_RGC_updateState = .true. mech_RGC_updateState = .true.
#ifdef DEBUG #ifdef DEBUG
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 & if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 &
.and. prm%of_debug == of) write(6,'(1x,a55,/)')'... done and happy' .and. prm%of_debug == of) write(6,'(1x,a55,/)')'... done and happy'
@ -513,7 +425,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! if residual blows-up => done but unhappy ! if residual blows-up => done but unhappy
elseif (residMax > relMax_RGC*stresMax .or. residMax > absMax_RGC) then ! try to restart when residual blows up exceeding maximum bound elseif (residMax > relMax_RGC*stresMax .or. residMax > absMax_RGC) then ! try to restart when residual blows up exceeding maximum bound
homogenization_RGC_updateState = [.true.,.false.] ! with direct cut-back mech_RGC_updateState = [.true.,.false.] ! with direct cut-back
#ifdef DEBUG #ifdef DEBUG
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 & if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 &
@ -656,9 +568,10 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! ... of the numerical viscosity traction "rmatrix" ! ... of the numerical viscosity traction "rmatrix"
allocate(rmatrix(3*nIntFaceTot,3*nIntFaceTot),source=0.0_pReal) allocate(rmatrix(3*nIntFaceTot,3*nIntFaceTot),source=0.0_pReal)
forall (i=1:3*nIntFaceTot) & do i=1,3*nIntFaceTot
rmatrix(i,i) = viscModus_RGC*viscPower_RGC/(refRelaxRate_RGC*dt)* & ! tangent due to numerical viscosity traction appears rmatrix(i,i) = viscModus_RGC*viscPower_RGC/(refRelaxRate_RGC*dt)* & ! tangent due to numerical viscosity traction appears
(abs(drelax(i))/(refRelaxRate_RGC*dt))**(viscPower_RGC - 1.0_pReal) ! only in the main diagonal term (abs(drelax(i))/(refRelaxRate_RGC*dt))**(viscPower_RGC - 1.0_pReal) ! only in the main diagonal term
enddo
#ifdef DEBUG #ifdef DEBUG
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0) then if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0) then
@ -710,7 +623,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
enddo; enddo enddo; enddo
stt%relaxationVector(:,of) = relax + drelax ! Updateing the state variable for the next iteration stt%relaxationVector(:,of) = relax + drelax ! Updateing the state variable for the next iteration
if (any(abs(drelax) > maxdRelax_RGC)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large if (any(abs(drelax) > maxdRelax_RGC)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large
homogenization_RGC_updateState = [.true.,.false.] mech_RGC_updateState = [.true.,.false.]
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(1x,a,1x,i3,1x,a,1x,i3,1x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback' write(6,'(1x,a,1x,i3,1x,a,1x,i3,1x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback'
write(6,'(1x,a,1x,e15.8)')'due to large relaxation change =',maxval(abs(drelax)) write(6,'(1x,a,1x,e15.8)')'due to large relaxation change =',maxval(abs(drelax))
@ -736,10 +649,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
!> @brief calculate stress-like penalty due to deformation mismatch !> @brief calculate stress-like penalty due to deformation mismatch
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine stressPenalty(rPen,nMis,avgF,fDef,ip,el,instance,of) subroutine stressPenalty(rPen,nMis,avgF,fDef,ip,el,instance,of)
use math, only: &
math_civita
use numerics, only: &
xSmoo_RGC
real(pReal), dimension (:,:,:), intent(out) :: rPen !< stress-like penalty real(pReal), dimension (:,:,:), intent(out) :: rPen !< stress-like penalty
real(pReal), dimension (:,:), intent(out) :: nMis !< total amount of mismatch real(pReal), dimension (:,:), intent(out) :: nMis !< total amount of mismatch
@ -852,13 +761,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
!> @brief calculate stress-like penalty due to volume discrepancy !> @brief calculate stress-like penalty due to volume discrepancy
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine volumePenalty(vPen,vDiscrep,fAvg,fDef,nGrain,instance,of) subroutine volumePenalty(vPen,vDiscrep,fAvg,fDef,nGrain,instance,of)
use math, only: &
math_det33, &
math_inv33
use numerics, only: &
maxVolDiscr_RGC,&
volDiscrMod_RGC,&
volDiscrPow_RGC
real(pReal), dimension (:,:,:), intent(out) :: vPen ! stress-like penalty due to volume real(pReal), dimension (:,:,:), intent(out) :: vPen ! stress-like penalty due to volume
real(pReal), intent(out) :: vDiscrep ! total volume discrepancy real(pReal), intent(out) :: vDiscrep ! total volume discrepancy
@ -907,8 +809,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
! deformation ! deformation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function surfaceCorrection(avgF,instance,of) function surfaceCorrection(avgF,instance,of)
use math, only: &
math_invert33
real(pReal), dimension(3) :: surfaceCorrection real(pReal), dimension(3) :: surfaceCorrection
@ -940,8 +840,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
!> @brief compute the equivalent shear and bulk moduli from the elasticity tensor !> @brief compute the equivalent shear and bulk moduli from the elasticity tensor
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function equivalentModuli(grainID,ip,el) function equivalentModuli(grainID,ip,el)
use constitutive, only: &
constitutive_homogenizedC
real(pReal), dimension(2) :: equivalentModuli real(pReal), dimension(2) :: equivalentModuli
@ -1012,13 +910,13 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
end subroutine grainDeformation end subroutine grainDeformation
end function homogenization_RGC_updateState end procedure mech_RGC_updateState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief derive average stress and stiffness from constituent quantities !> @brief derive average stress and stiffness from constituent quantities
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
@ -1030,68 +928,19 @@ subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,
avgP = sum(P,3) /real(product(param(instance)%Nconstituents),pReal) avgP = sum(P,3) /real(product(param(instance)%Nconstituents),pReal)
dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%Nconstituents),pReal) dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%Nconstituents),pReal)
end subroutine homogenization_RGC_averageStressAndItsTangent end subroutine mech_RGC_averageStressAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief return array of homogenization results for post file inclusion
!--------------------------------------------------------------------------------------------------
pure function homogenization_RGC_postResults(instance,of) result(postResults)
integer, intent(in) :: &
instance, &
of
integer :: &
o,c
real(pReal), dimension(sum(homogenization_RGC_sizePostResult(:,instance))) :: &
postResults
associate(stt => state(instance), dst => dependentState(instance), prm => param(instance))
c = 0
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (constitutivework_ID)
postResults(c+1) = stt%work(of)
c = c + 1
case (magnitudemismatch_ID)
postResults(c+1:c+3) = dst%mismatch(1:3,of)
c = c + 3
case (penaltyenergy_ID)
postResults(c+1) = stt%penaltyEnergy(of)
c = c + 1
case (volumediscrepancy_ID)
postResults(c+1) = dst%volumeDiscrepancy(of)
c = c + 1
case (averagerelaxrate_ID)
postResults(c+1) = dst%relaxationrate_avg(of)
c = c + 1
case (maximumrelaxrate_ID)
postResults(c+1) = dst%relaxationrate_max(of)
c = c + 1
end select
enddo outputsLoop
end associate
end function homogenization_RGC_postResults
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file !> @brief writes results to HDF5 output file
! ToDo: check wheter units are correct ! ToDo: check wheter units are correct
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mech_RGC_results(instance,group) module subroutine mech_RGC_results(instance,group)
#if defined(PETSc) || defined(DAMASK_HDF5) #if defined(PETSc) || defined(DAMASK_HDF5)
use results, only: &
results_writeDataset
integer, intent(in) :: instance integer, intent(in) :: instance
character(len=*) :: group character(len=*), intent(in) :: group
integer :: o integer :: o
associate(stt => state(instance), dst => dependentState(instance), prm => param(instance)) associate(stt => state(instance), dst => dependentState(instance), prm => param(instance))
@ -1122,8 +971,8 @@ subroutine mech_RGC_results(instance,group)
end associate end associate
#else #else
integer, intent(in) :: instance integer, intent(in) :: instance
character(len=*) :: group character(len=*), intent(in) :: group
#endif #endif
end subroutine mech_RGC_results end subroutine mech_RGC_results
@ -1262,7 +1111,7 @@ integer pure function interface4to1(iFace4D, nGDim)
else else
interface4to1 = iFace4D(4) + nGDim(3)*(iFace4D(2)-1) & interface4to1 = iFace4D(4) + nGDim(3)*(iFace4D(2)-1) &
+ nGDim(3)*nGDim(1)*(iFace4D(3)-1) & + nGDim(3)*nGDim(1)*(iFace4D(3)-1) &
+ (nGDim(1)-1)*nGDim(2)*nGDim(3) ! total number of interfaces normal //e1 + (nGDim(1)-1)*nGDim(2)*nGDim(3) ! total # of interfaces normal || e1
endif endif
case(3) case(3)
@ -1271,8 +1120,8 @@ integer pure function interface4to1(iFace4D, nGDim)
else else
interface4to1 = iFace4D(2) + nGDim(1)*(iFace4D(3)-1) & interface4to1 = iFace4D(2) + nGDim(1)*(iFace4D(3)-1) &
+ nGDim(1)*nGDim(2)*(iFace4D(4)-1) & + nGDim(1)*nGDim(2)*(iFace4D(4)-1) &
+ (nGDim(1)-1)*nGDim(2)*nGDim(3) & ! total number of interfaces normal //e1 + (nGDim(1)-1)*nGDim(2)*nGDim(3) & ! total # of interfaces normal || e1
+ nGDim(1)*(nGDim(2)-1)*nGDim(3) ! total number of interfaces normal //e2 + nGDim(1)*(nGDim(2)-1)*nGDim(3) ! total # of interfaces normal || e2
endif endif
case default case default
@ -1296,23 +1145,23 @@ pure function interface1to4(iFace1D, nGDim)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! compute the total number of interfaces, which ... ! compute the total number of interfaces, which ...
nIntFace = [(nGDim(1)-1)*nGDim(2)*nGDim(3), & ! ... normal //e1 nIntFace = [(nGDim(1)-1)*nGDim(2)*nGDim(3), & ! ... normal || e1
nGDim(1)*(nGDim(2)-1)*nGDim(3), & ! ... normal //e2 nGDim(1)*(nGDim(2)-1)*nGDim(3), & ! ... normal || e2
nGDim(1)*nGDim(2)*(nGDim(3)-1)] ! ... normal //e3 nGDim(1)*nGDim(2)*(nGDim(3)-1)] ! ... normal || e3
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! get the corresponding interface ID in 4D (normal and local position) ! get the corresponding interface ID in 4D (normal and local position)
if (iFace1D > 0 .and. iFace1D <= nIntFace(1)) then ! interface with normal //e1 if (iFace1D > 0 .and. iFace1D <= nIntFace(1)) then ! interface with normal || e1
interface1to4(1) = 1 interface1to4(1) = 1
interface1to4(3) = mod((iFace1D-1),nGDim(2))+1 interface1to4(3) = mod((iFace1D-1),nGDim(2))+1
interface1to4(4) = mod(int(real(iFace1D-1,pReal)/real(nGDim(2),pReal)),nGDim(3))+1 interface1to4(4) = mod(int(real(iFace1D-1,pReal)/real(nGDim(2),pReal)),nGDim(3))+1
interface1to4(2) = int(real(iFace1D-1,pReal)/real(nGDim(2),pReal)/real(nGDim(3),pReal))+1 interface1to4(2) = int(real(iFace1D-1,pReal)/real(nGDim(2),pReal)/real(nGDim(3),pReal))+1
elseif (iFace1D > nIntFace(1) .and. iFace1D <= (nIntFace(2) + nIntFace(1))) then ! interface with normal //e2 elseif (iFace1D > nIntFace(1) .and. iFace1D <= (nIntFace(2) + nIntFace(1))) then ! interface with normal || e2
interface1to4(1) = 2 interface1to4(1) = 2
interface1to4(4) = mod((iFace1D-nIntFace(1)-1),nGDim(3))+1 interface1to4(4) = mod((iFace1D-nIntFace(1)-1),nGDim(3))+1
interface1to4(2) = mod(int(real(iFace1D-nIntFace(1)-1,pReal)/real(nGDim(3),pReal)),nGDim(1))+1 interface1to4(2) = mod(int(real(iFace1D-nIntFace(1)-1,pReal)/real(nGDim(3),pReal)),nGDim(1))+1
interface1to4(3) = int(real(iFace1D-nIntFace(1)-1,pReal)/real(nGDim(3),pReal)/real(nGDim(1),pReal))+1 interface1to4(3) = int(real(iFace1D-nIntFace(1)-1,pReal)/real(nGDim(3),pReal)/real(nGDim(1),pReal))+1
elseif (iFace1D > nIntFace(2) + nIntFace(1) .and. iFace1D <= (nIntFace(3) + nIntFace(2) + nIntFace(1))) then ! interface with normal //e3 elseif (iFace1D > nIntFace(2) + nIntFace(1) .and. iFace1D <= (nIntFace(3) + nIntFace(2) + nIntFace(1))) then ! interface with normal || e3
interface1to4(1) = 3 interface1to4(1) = 3
interface1to4(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1),nGDim(1))+1 interface1to4(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1),nGDim(1))+1
interface1to4(3) = mod(int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/real(nGDim(1),pReal)),nGDim(2))+1 interface1to4(3) = mod(int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/real(nGDim(1),pReal)),nGDim(2))+1
@ -1322,4 +1171,4 @@ pure function interface1to4(iFace1D, nGDim)
end function interface1to4 end function interface1to4
end module homogenization_mech_RGC end submodule homogenization_mech_RGC

10
src/homogenization_mech_isostrain.f90
View File

@ -6,8 +6,6 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(homogenization) homogenization_mech_isostrain submodule(homogenization) homogenization_mech_isostrain
implicit none
enum, bind(c) enum, bind(c)
enumerator :: & enumerator :: &
parallel_ID, & parallel_ID, &
@ -30,14 +28,6 @@ contains
!> @brief allocates all neccessary fields, reads information from material configuration file !> @brief allocates all neccessary fields, reads information from material configuration file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_isostrain_init module subroutine mech_isostrain_init
use debug, only: &
debug_HOMOGENIZATION, &
debug_level, &
debug_levelBasic
use IO, only: &
IO_error
use config, only: &
config_homogenization
integer :: & integer :: &
Ninstance, & Ninstance, &

8
src/homogenization_mech_none.f90
View File

@ -6,20 +6,12 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(homogenization) homogenization_mech_none submodule(homogenization) homogenization_mech_none
implicit none
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields, reads information from material configuration file !> @brief allocates all neccessary fields, reads information from material configuration file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_none_init module subroutine mech_none_init
use debug, only: &
debug_HOMOGENIZATION, &
debug_level, &
debug_levelBasic
use config, only: &
config_homogenization
integer :: & integer :: &
Ninstance, & Ninstance, &

4
src/plastic_dislotwin.f90
View File

@ -1067,8 +1067,8 @@ subroutine plastic_dislotwin_results(instance,group)
case (rho_dip_ID) case (rho_dip_ID)
call results_writeDataset(group,stt%rho_dip,'rho_dip',& call results_writeDataset(group,stt%rho_dip,'rho_dip',&
'dislocation dipole density''1/m²') 'dislocation dipole density''1/m²')
case (dot_gamma_sl_ID) case (gamma_sl_ID)
call results_writeDataset(group,stt%gamma_sl,'dot_gamma_sl',& call results_writeDataset(group,stt%gamma_sl,'gamma_sl',&
'plastic shear','1') 'plastic shear','1')
case (Lambda_sl_ID) case (Lambda_sl_ID)
call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',& call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',&

18
src/results.f90
View File

@ -299,18 +299,26 @@ end subroutine results_writeVectorDataset_real
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine results_writeTensorDataset_real(group,dataset,label,description,SIunit) subroutine results_writeTensorDataset_real(group,dataset,label,description,SIunit)
character(len=*), intent(in) :: label,group,description character(len=*), intent(in) :: label,group,description
character(len=*), intent(in), optional :: SIunit character(len=*), intent(in), optional :: SIunit
real(pReal), intent(inout), dimension(:,:,:) :: dataset real(pReal), intent(in), dimension(:,:,:) :: dataset
integer :: i
integer(HID_T) :: groupHandle integer(HID_T) :: groupHandle
real(pReal), dimension(:,:,:), allocatable :: dataset_transposed
allocate(dataset_transposed,mold=dataset)
do i=1,size(dataset,3)
dataset_transposed(1:3,1:3,i) = transpose(dataset(1:3,1:3,i))
enddo
groupHandle = results_openGroup(group) groupHandle = results_openGroup(group)
#ifdef PETSc #ifdef PETSc
call HDF5_write(groupHandle,dataset,label,.true.) call HDF5_write(groupHandle,dataset_transposed,label,.true.)
#else #else
call HDF5_write(groupHandle,dataset,label,.false.) call HDF5_write(groupHandle,dataset_transposed,label,.false.)
#endif #endif
if (HDF5_objectExists(groupHandle,label)) & if (HDF5_objectExists(groupHandle,label)) &