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

Merge remote-tracking branch 'origin/development' into python-improvements

This commit is contained in:
Sharan Roongta 2021-02-13 16:51:29 +01:00
parent fe63adcca4 2676837578
commit 099cd90890
52 changed files with 859 additions and 827 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit d90babadfb0a33afa1b793044cc4efd4d7430731 Subproject commit e74cf00628285a587ced1e551cc9837c1011ca1c

2
VERSION
View File

@ -1 +1 @@
v3.0.0-alpha2-415-g09c2d7f3f v3.0.0-alpha2-428-gb7b764ab5

3
examples/SpectralMethod/Polycrystal/shearXY.yaml
View File

@ -3,9 +3,6 @@
solver: solver:
mechanical: spectral_basic mechanical: spectral_basic
initial_conditions:
T: 300 #in Kelvin
loadstep: loadstep:
- boundary_conditions: - boundary_conditions:
mechanical: mechanical:

3
examples/SpectralMethod/Polycrystal/shearZX.yaml
View File

@ -3,9 +3,6 @@
solver: solver:
mechanical: spectral_basic mechanical: spectral_basic
initial_conditions:
T: 300 #in Kelvin
loadstep: loadstep:
- boundary_conditions: - boundary_conditions:
mechanical: mechanical:

3
examples/SpectralMethod/Polycrystal/tensionX.yaml
View File

@ -3,9 +3,6 @@
solver: solver:
mechanical: spectral_basic mechanical: spectral_basic
initial_conditions:
T: 300 #in Kelvin
loadstep: loadstep:
- boundary_conditions: - boundary_conditions:
mechanical: mechanical:

81
python/damask/_colormap.py
View File

@ -223,25 +223,46 @@ class Colormap(mpl.colors.ListedColormap):
return Colormap(np.array(rev.colors),rev.name[:-4] if rev.name.endswith('_r_r') else rev.name) return Colormap(np.array(rev.colors),rev.name[:-4] if rev.name.endswith('_r_r') else rev.name)
def _get_file_handle(self,fname,extension):
"""
Provide file handle.
Parameters
----------
fname : file, str, pathlib.Path, or None
Filename or filehandle, will be name of the colormap+extension if None.
extension: str
Extension of the filename.
Returns
-------
f
File handle
"""
if fname is None:
fhandle = open(self.name.replace(' ','_')+'.'+extension,'w',newline='\n')
else:
try:
fhandle = open(fname,'w',newline='\n')
except TypeError:
fhandle = fname
return fhandle
def save_paraview(self,fname=None): def save_paraview(self,fname=None):
""" """
Save as JSON file for use in Paraview. Save as JSON file for use in Paraview.
Parameters Parameters
---------- ----------
fname : file, str, or pathlib.Path, optional. fname : file, str, or pathlib.Path, optional
Filename to store results. If not given, the filename will Filename to store results. If not given, the filename will
consist of the name of the colormap and extension '.json'. consist of the name of the colormap and extension '.json'.
""" """
if fname is None:
fhandle = None
else:
try:
fhandle = open(fname,'w')
except TypeError:
fhandle = fname
colors = [] colors = []
for i,c in enumerate(np.round(self.colors,6).tolist()): for i,c in enumerate(np.round(self.colors,6).tolist()):
colors+=[i]+c colors+=[i]+c
@ -254,8 +275,7 @@ class Colormap(mpl.colors.ListedColormap):
'RGBPoints':colors 'RGBPoints':colors
}] }]
with open(self.name.replace(' ','_')+'.json', 'w') if fhandle is None else fhandle as f: json.dump(out,self._get_file_handle(fname,'json'),indent=4)
json.dump(out, f,indent=4)
def save_ASCII(self,fname=None): def save_ASCII(self,fname=None):
@ -264,24 +284,14 @@ class Colormap(mpl.colors.ListedColormap):
Parameters Parameters
---------- ----------
fname : file, str, or pathlib.Path, optional. fname : file, str, or pathlib.Path, optional
Filename to store results. If not given, the filename will Filename to store results. If not given, the filename will
consist of the name of the colormap and extension '.txt'. consist of the name of the colormap and extension '.txt'.
""" """
if fname is None:
fhandle = None
else:
try:
fhandle = open(fname,'w')
except TypeError:
fhandle = fname
labels = {'RGBA':4} if self.colors.shape[1] == 4 else {'RGB': 3} labels = {'RGBA':4} if self.colors.shape[1] == 4 else {'RGB': 3}
t = Table(self.colors,labels,f'Creator: {util.execution_stamp("Colormap")}') t = Table(self.colors,labels,f'Creator: {util.execution_stamp("Colormap")}')
t.save(self._get_file_handle(fname,'txt'))
with open(self.name.replace(' ','_')+'.txt', 'w') if fhandle is None else fhandle as f:
t.save(f)
def save_GOM(self,fname=None): def save_GOM(self,fname=None):
@ -290,26 +300,19 @@ class Colormap(mpl.colors.ListedColormap):
Parameters Parameters
---------- ----------
fname : file, str, or pathlib.Path, optional. fname : file, str, or pathlib.Path, optional
Filename to store results. If not given, the filename will Filename to store results. If not given, the filename will
consist of the name of the colormap and extension '.legend'. consist of the name of the colormap and extension '.legend'.
""" """
if fname is None:
fhandle = None
else:
try:
fhandle = open(fname,'w')
except TypeError:
fhandle = fname
# ToDo: test in GOM # ToDo: test in GOM
GOM_str = '1 1 {name} 9 {name} '.format(name=self.name.replace(" ","_")) \ GOM_str = '1 1 {name} 9 {name} '.format(name=self.name.replace(" ","_")) \
+ '0 1 0 3 0 0 -1 9 \\ 0 0 0 255 255 255 0 0 255 ' \ + '0 1 0 3 0 0 -1 9 \\ 0 0 0 255 255 255 0 0 255 ' \
+ f'30 NO_UNIT 1 1 64 64 64 255 1 0 0 0 0 0 0 3 0 {len(self.colors)}' \ + f'30 NO_UNIT 1 1 64 64 64 255 1 0 0 0 0 0 0 3 0 {len(self.colors)}' \
+ ' '.join([f' 0 {c[0]} {c[1]} {c[2]} 255 1' for c in reversed((self.colors*255).astype(int))]) \ + ' '.join([f' 0 {c[0]} {c[1]} {c[2]} 255 1' for c in reversed((self.colors*255).astype(int))]) \
+ '\n' + '\n'
with open(self.name.replace(' ','_')+'.legend', 'w') if fhandle is None else fhandle as f:
f.write(GOM_str) self._get_file_handle(fname,'legend').write(GOM_str)
def save_gmsh(self,fname=None): def save_gmsh(self,fname=None):
@ -318,24 +321,16 @@ class Colormap(mpl.colors.ListedColormap):
Parameters Parameters
---------- ----------
fname : file, str, or pathlib.Path, optional. fname : file, str, or pathlib.Path, optional
Filename to store results. If not given, the filename will Filename to store results. If not given, the filename will
consist of the name of the colormap and extension '.msh'. consist of the name of the colormap and extension '.msh'.
""" """
if fname is None:
fhandle = None
else:
try:
fhandle = open(fname,'w')
except TypeError:
fhandle = fname
# ToDo: test in gmsh # ToDo: test in gmsh
gmsh_str = 'View.ColorTable = {\n' \ gmsh_str = 'View.ColorTable = {\n' \
+'\n'.join([f'{c[0]},{c[1]},{c[2]},' for c in self.colors[:,:3]*255]) \ +'\n'.join([f'{c[0]},{c[1]},{c[2]},' for c in self.colors[:,:3]*255]) \
+'\n}\n' +'\n}\n'
with open(self.name.replace(' ','_')+'.msh', 'w') if fhandle is None else fhandle as f: self._get_file_handle(fname,'msh').write(gmsh_str)
f.write(gmsh_str)
@staticmethod @staticmethod

2
python/damask/_config.py
View File

@ -75,7 +75,7 @@ class Config(dict):
""" """
try: try:
fhandle = open(fname,'w') fhandle = open(fname,'w',newline='\n')
except TypeError: except TypeError:
fhandle = fname fhandle = fname

14
python/damask/_configmaterial.py
View File

@ -168,11 +168,11 @@ class ConfigMaterial(Config):
ok = False ok = False
if 'material' in self: if 'material' in self:
for i,v in enumerate(self['material']): for i,m in enumerate(self['material']):
if 'constituents' in v: if 'constituents' in m:
f = 0.0 v = 0.0
for c in v['constituents']: for c in m['constituents']:
f+= float(c['v']) v+= float(c['v'])
if 'O' in c: if 'O' in c:
try: try:
Rotation.from_quaternion(c['O']) Rotation.from_quaternion(c['O'])
@ -180,8 +180,8 @@ class ConfigMaterial(Config):
o = c['O'] o = c['O']
print(f"Invalid orientation: '{o}' in material '{i}'") print(f"Invalid orientation: '{o}' in material '{i}'")
ok = False ok = False
if not np.isclose(f,1.0): if not np.isclose(v,1.0):
print(f"Invalid total fraction '{f}' in material '{i}'") print(f"Invalid total fraction (v) '{v}' in material '{i}'")
ok = False ok = False
return ok return ok

2
python/damask/_result.py
View File

@ -1287,7 +1287,7 @@ class Result:
np.prod(shape))} np.prod(shape))}
data_items[-1].text=f'{os.path.split(self.fname)[1]}:{name}' data_items[-1].text=f'{os.path.split(self.fname)[1]}:{name}'
with open(self.fname.with_suffix('.xdmf').name,'w') as f: with open(self.fname.with_suffix('.xdmf').name,'w',newline='\n') as f:
f.write(xml.dom.minidom.parseString(ET.tostring(xdmf).decode()).toprettyxml()) f.write(xml.dom.minidom.parseString(ET.tostring(xdmf).decode()).toprettyxml())

4
python/damask/_table.py
View File

@ -26,7 +26,7 @@ class Table:
comments_ = [comments] if isinstance(comments,str) else comments comments_ = [comments] if isinstance(comments,str) else comments
self.comments = [] if comments_ is None else [c for c in comments_] self.comments = [] if comments_ is None else [c for c in comments_]
self.data = pd.DataFrame(data=data) self.data = pd.DataFrame(data=data)
self.shapes = { k:(v,) if isinstance(v,(np.int,int)) else v for k,v in shapes.items() } self.shapes = { k:(v,) if isinstance(v,(np.int64,np.int32,int)) else v for k,v in shapes.items() }
self._label_uniform() self._label_uniform()
def __repr__(self): def __repr__(self):
@ -380,7 +380,7 @@ class Table:
[f'# {comment}' for comment in self.comments] [f'# {comment}' for comment in self.comments]
try: try:
fhandle = open(fname,'w') fhandle = open(fname,'w',newline='\n')
except TypeError: except TypeError:
fhandle = fname fhandle = fname

4
python/damask/_vtk.py
View File

@ -347,7 +347,6 @@ class VTK:
See http://compilatrix.com/article/vtk-1 for further ideas. See http://compilatrix.com/article/vtk-1 for further ideas.
""" """
def screen_size():
try: try:
import wx import wx
_ = wx.App(False) # noqa _ = wx.App(False) # noqa
@ -363,7 +362,6 @@ class VTK:
width = 1024 width = 1024
height = 768 height = 768
return (width,height)
mapper = vtk.vtkDataSetMapper() mapper = vtk.vtkDataSetMapper()
mapper.SetInputData(self.vtk_data) mapper.SetInputData(self.vtk_data)
actor = vtk.vtkActor() actor = vtk.vtkActor()
@ -377,7 +375,7 @@ class VTK:
ren.AddActor(actor) ren.AddActor(actor)
ren.SetBackground(0.2,0.2,0.2) ren.SetBackground(0.2,0.2,0.2)
window.SetSize(screen_size[0],screen_size[1]) window.SetSize(width,height)
iren = vtk.vtkRenderWindowInteractor() iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(window) iren.SetRenderWindow(window)

2
python/damask/solver/_marc.py
View File

@ -66,7 +66,7 @@ class Marc:
if logfile is not None: if logfile is not None:
try: try:
f = open(logfile,'w+') f = open(logfile,'w+',newline='\n')
except TypeError: except TypeError:
f = logfile f = logfile
else: else:

6
src/CPFEM.f90
View File

@ -85,7 +85,7 @@ subroutine CPFEM_initAll
call discretization_marc_init call discretization_marc_init
call lattice_init call lattice_init
call material_init(.false.) call material_init(.false.)
call constitutive_init call phase_init
call homogenization_init call homogenization_init
call crystallite_init call crystallite_init
call CPFEM_init call CPFEM_init
@ -257,7 +257,7 @@ end subroutine CPFEM_general
subroutine CPFEM_forward subroutine CPFEM_forward
call homogenization_forward call homogenization_forward
call constitutive_forward call phase_forward
end subroutine CPFEM_forward end subroutine CPFEM_forward
@ -272,7 +272,7 @@ subroutine CPFEM_results(inc,time)
call results_openJobFile call results_openJobFile
call results_addIncrement(inc,time) call results_addIncrement(inc,time)
call constitutive_results call phase_results
call homogenization_results call homogenization_results
call discretization_results call discretization_results
call results_finalizeIncrement call results_finalizeIncrement

10
src/CPFEM2.f90
View File

@ -60,7 +60,7 @@ subroutine CPFEM_initAll
call discretization_grid_init(restart=interface_restartInc>0) call discretization_grid_init(restart=interface_restartInc>0)
#endif #endif
call material_init(restart=interface_restartInc>0) call material_init(restart=interface_restartInc>0)
call constitutive_init call phase_init
call homogenization_init call homogenization_init
call crystallite_init call crystallite_init
call CPFEM_init call CPFEM_init
@ -87,7 +87,7 @@ subroutine CPFEM_init
fileHandle = HDF5_openFile(fileName) fileHandle = HDF5_openFile(fileName)
call homogenization_restartRead(fileHandle) call homogenization_restartRead(fileHandle)
call constitutive_restartRead(fileHandle) call phase_restartRead(fileHandle)
call HDF5_closeFile(fileHandle) call HDF5_closeFile(fileHandle)
endif endif
@ -110,7 +110,7 @@ subroutine CPFEM_restartWrite
fileHandle = HDF5_openFile(fileName,'a') fileHandle = HDF5_openFile(fileName,'a')
call homogenization_restartWrite(fileHandle) call homogenization_restartWrite(fileHandle)
call constitutive_restartWrite(fileHandle) call phase_restartWrite(fileHandle)
call HDF5_closeFile(fileHandle) call HDF5_closeFile(fileHandle)
@ -123,7 +123,7 @@ end subroutine CPFEM_restartWrite
subroutine CPFEM_forward subroutine CPFEM_forward
call homogenization_forward call homogenization_forward
call constitutive_forward call phase_forward
end subroutine CPFEM_forward end subroutine CPFEM_forward
@ -138,7 +138,7 @@ subroutine CPFEM_results(inc,time)
call results_openJobFile call results_openJobFile
call results_addIncrement(inc,time) call results_addIncrement(inc,time)
call constitutive_results call phase_results
call homogenization_results call homogenization_results
call discretization_results call discretization_results
call results_finalizeIncrement call results_finalizeIncrement

27
src/IO.f90
View File

@ -65,8 +65,8 @@ end subroutine IO_init
function IO_readlines(fileName) result(fileContent) function IO_readlines(fileName) result(fileContent)
character(len=*), intent(in) :: fileName character(len=*), intent(in) :: fileName
character(len=pStringLen), dimension(:), allocatable :: fileContent !< file content, separated per lines character(len=pStringLen), dimension(:), allocatable :: fileContent !< file content, separated per lines
character(len=pStringLen) :: line character(len=pStringLen) :: line
character(len=:), allocatable :: rawData character(len=:), allocatable :: rawData
integer :: & integer :: &
@ -75,6 +75,7 @@ function IO_readlines(fileName) result(fileContent)
l l
logical :: warned logical :: warned
rawData = IO_read(fileName) rawData = IO_read(fileName)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -112,16 +113,21 @@ end function IO_readlines
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Read whole file. !> @brief Read whole file.
!> @details ensures that the string ends with a new line (expected UNIX behavior) !> @details ensures that the string ends with a new line (expected UNIX behavior) and rejects
! windows (CRLF) line endings
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function IO_read(fileName) result(fileContent) function IO_read(fileName) result(fileContent)
character(len=*), intent(in) :: fileName character(len=*), intent(in) :: fileName
character(len=:), allocatable :: fileContent character(len=:), allocatable :: fileContent
integer :: & integer :: &
fileLength, & fileLength, &
fileUnit, & fileUnit, &
myStat myStat, &
firstEOL
character, parameter :: CR = achar(13)
inquire(file = fileName, size=fileLength) inquire(file = fileName, size=fileLength)
open(newunit=fileUnit, file=fileName, access='stream',& open(newunit=fileUnit, file=fileName, access='stream',&
@ -137,8 +143,12 @@ function IO_read(fileName) result(fileContent)
if(myStat /= 0) call IO_error(102,ext_msg=trim(fileName)) if(myStat /= 0) call IO_error(102,ext_msg=trim(fileName))
close(fileUnit) close(fileUnit)
if(fileContent(fileLength:fileLength) /= IO_EOL) fileContent = fileContent//IO_EOL ! ensure EOL@EOF if(fileContent(fileLength:fileLength) /= IO_EOL) fileContent = fileContent//IO_EOL ! ensure EOL@EOF
firstEOL = index(fileContent,IO_EOL)
if(scan(fileContent(firstEOL:firstEOL),CR) /= 0) call IO_error(115)
end function IO_read end function IO_read
@ -151,6 +161,7 @@ logical pure function IO_isBlank(string)
integer :: posNonBlank integer :: posNonBlank
posNonBlank = verify(string,IO_WHITESPACE) posNonBlank = verify(string,IO_WHITESPACE)
IO_isBlank = posNonBlank == 0 .or. posNonBlank == scan(string,IO_COMMENT) IO_isBlank = posNonBlank == 0 .or. posNonBlank == scan(string,IO_COMMENT)
@ -170,6 +181,7 @@ pure function IO_stringPos(string)
integer :: left, right integer :: left, right
allocate(IO_stringPos(1), source=0) allocate(IO_stringPos(1), source=0)
right = 0 right = 0
@ -249,6 +261,7 @@ pure function IO_lc(string)
integer :: i,n integer :: i,n
do i=1,len(string) do i=1,len(string)
n = index(UPPER,string(i:i)) n = index(UPPER,string(i:i))
if(n/=0) then if(n/=0) then
@ -271,6 +284,7 @@ function IO_rmComment(line)
character(len=:), allocatable :: IO_rmComment character(len=:), allocatable :: IO_rmComment
integer :: split integer :: split
split = index(line,IO_COMMENT) split = index(line,IO_COMMENT)
if (split == 0) then if (split == 0) then
@ -292,6 +306,7 @@ integer function IO_stringAsInt(string)
integer :: readStatus integer :: readStatus
character(len=*), parameter :: VALIDCHARS = '0123456789+- ' character(len=*), parameter :: VALIDCHARS = '0123456789+- '
valid: if (verify(string,VALIDCHARS) == 0) then valid: if (verify(string,VALIDCHARS) == 0) then
read(string,*,iostat=readStatus) IO_stringAsInt read(string,*,iostat=readStatus) IO_stringAsInt
if (readStatus /= 0) call IO_error(111,ext_msg=string) if (readStatus /= 0) call IO_error(111,ext_msg=string)
@ -313,6 +328,7 @@ real(pReal) function IO_stringAsFloat(string)
integer :: readStatus integer :: readStatus
character(len=*), parameter :: VALIDCHARS = '0123456789eE.+- ' character(len=*), parameter :: VALIDCHARS = '0123456789eE.+- '
valid: if (verify(string,VALIDCHARS) == 0) then valid: if (verify(string,VALIDCHARS) == 0) then
read(string,*,iostat=readStatus) IO_stringAsFloat read(string,*,iostat=readStatus) IO_stringAsFloat
if (readStatus /= 0) call IO_error(112,ext_msg=string) if (readStatus /= 0) call IO_error(112,ext_msg=string)
@ -331,6 +347,7 @@ logical function IO_stringAsBool(string)
character(len=*), intent(in) :: string !< string for conversion to int value character(len=*), intent(in) :: string !< string for conversion to int value
if (trim(adjustl(string)) == 'True' .or. trim(adjustl(string)) == 'true') then if (trim(adjustl(string)) == 'True' .or. trim(adjustl(string)) == 'true') then
IO_stringAsBool = .true. IO_stringAsBool = .true.
elseif (trim(adjustl(string)) == 'False' .or. trim(adjustl(string)) == 'false') then elseif (trim(adjustl(string)) == 'False' .or. trim(adjustl(string)) == 'false') then
@ -356,6 +373,7 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
character(len=:), allocatable :: msg character(len=:), allocatable :: msg
character(len=pStringLen) :: formatString character(len=pStringLen) :: formatString
select case (error_ID) select case (error_ID)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -382,6 +400,9 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
msg = 'invalid character for logical:' msg = 'invalid character for logical:'
case (114) case (114)
msg = 'cannot decode base64 string:' msg = 'cannot decode base64 string:'
case (115)
msg = 'found CR. Windows file endings (CRLF) are not supported.'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! lattice error messages ! lattice error messages

34
src/commercialFEM_fileList.f90
View File

@ -20,19 +20,19 @@
#include "material.f90" #include "material.f90"
#include "lattice.f90" #include "lattice.f90"
#include "phase.f90" #include "phase.f90"
#include "phase_mechanics.f90" #include "phase_mechanical.f90"
#include "phase_mechanics_plastic.f90" #include "phase_mechanical_plastic.f90"
#include "phase_mechanics_plastic_none.f90" #include "phase_mechanical_plastic_none.f90"
#include "phase_mechanics_plastic_isotropic.f90" #include "phase_mechanical_plastic_isotropic.f90"
#include "phase_mechanics_plastic_phenopowerlaw.f90" #include "phase_mechanical_plastic_phenopowerlaw.f90"
#include "phase_mechanics_plastic_kinehardening.f90" #include "phase_mechanical_plastic_kinehardening.f90"
#include "phase_mechanics_plastic_dislotwin.f90" #include "phase_mechanical_plastic_dislotwin.f90"
#include "phase_mechanics_plastic_dislotungsten.f90" #include "phase_mechanical_plastic_dislotungsten.f90"
#include "phase_mechanics_plastic_nonlocal.f90" #include "phase_mechanical_plastic_nonlocal.f90"
#include "phase_mechanics_eigendeformation.f90" #include "phase_mechanical_eigen.f90"
#include "phase_mechanics_eigendeformation_cleavageopening.f90" #include "phase_mechanical_eigen_cleavageopening.f90"
#include "phase_mechanics_eigendeformation_slipplaneopening.f90" #include "phase_mechanical_eigen_slipplaneopening.f90"
#include "phase_mechanics_eigendeformation_thermalexpansion.f90" #include "phase_mechanical_eigen_thermalexpansion.f90"
#include "phase_thermal.f90" #include "phase_thermal.f90"
#include "phase_thermal_dissipation.f90" #include "phase_thermal_dissipation.f90"
#include "phase_thermal_externalheat.f90" #include "phase_thermal_externalheat.f90"
@ -44,10 +44,10 @@
#include "damage_none.f90" #include "damage_none.f90"
#include "damage_nonlocal.f90" #include "damage_nonlocal.f90"
#include "homogenization.f90" #include "homogenization.f90"
#include "homogenization_mechanics.f90" #include "homogenization_mechanical.f90"
#include "homogenization_mechanics_none.f90" #include "homogenization_mechanical_pass.f90"
#include "homogenization_mechanics_isostrain.f90" #include "homogenization_mechanical_isostrain.f90"
#include "homogenization_mechanics_RGC.f90" #include "homogenization_mechanical_RGC.f90"
#include "homogenization_thermal.f90" #include "homogenization_thermal.f90"
#include "homogenization_damage.f90" #include "homogenization_damage.f90"
#include "CPFEM.f90" #include "CPFEM.f90"

66
src/grid/DAMASK_grid.f90
View File

@ -18,9 +18,9 @@ program DAMASK_grid
use CPFEM2 use CPFEM2
use material use material
use spectral_utilities use spectral_utilities
use grid_mech_spectral_basic use grid_mechanical_spectral_basic
use grid_mech_spectral_polarisation use grid_mechanical_spectral_polarisation
use grid_mech_FEM use grid_mechanical_FEM
use grid_damage_spectral use grid_damage_spectral
use grid_thermal_spectral use grid_thermal_spectral
use results use results
@ -83,16 +83,16 @@ program DAMASK_grid
type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases
type(tSolutionState), allocatable, dimension(:) :: solres type(tSolutionState), allocatable, dimension(:) :: solres
procedure(grid_mech_spectral_basic_init), pointer :: & procedure(grid_mechanical_spectral_basic_init), pointer :: &
mech_init mechanical_init
procedure(grid_mech_spectral_basic_forward), pointer :: & procedure(grid_mechanical_spectral_basic_forward), pointer :: &
mech_forward mechanical_forward
procedure(grid_mech_spectral_basic_solution), pointer :: & procedure(grid_mechanical_spectral_basic_solution), pointer :: &
mech_solution mechanical_solution
procedure(grid_mech_spectral_basic_updateCoords), pointer :: & procedure(grid_mechanical_spectral_basic_updateCoords), pointer :: &
mech_updateCoords mechanical_updateCoords
procedure(grid_mech_spectral_basic_restartWrite), pointer :: & procedure(grid_mechanical_spectral_basic_restartWrite), pointer :: &
mech_restartWrite mechanical_restartWrite
external :: & external :: &
quit quit
@ -139,25 +139,25 @@ program DAMASK_grid
debug_grid => config_debug%get('grid',defaultVal=emptyList) debug_grid => config_debug%get('grid',defaultVal=emptyList)
select case (trim(num_grid%get_asString('solver', defaultVal = 'Basic'))) select case (trim(num_grid%get_asString('solver', defaultVal = 'Basic')))
case ('Basic') case ('Basic')
mech_init => grid_mech_spectral_basic_init mechanical_init => grid_mechanical_spectral_basic_init
mech_forward => grid_mech_spectral_basic_forward mechanical_forward => grid_mechanical_spectral_basic_forward
mech_solution => grid_mech_spectral_basic_solution mechanical_solution => grid_mechanical_spectral_basic_solution
mech_updateCoords => grid_mech_spectral_basic_updateCoords mechanical_updateCoords => grid_mechanical_spectral_basic_updateCoords
mech_restartWrite => grid_mech_spectral_basic_restartWrite mechanical_restartWrite => grid_mechanical_spectral_basic_restartWrite
case ('Polarisation') case ('Polarisation')
mech_init => grid_mech_spectral_polarisation_init mechanical_init => grid_mechanical_spectral_polarisation_init
mech_forward => grid_mech_spectral_polarisation_forward mechanical_forward => grid_mechanical_spectral_polarisation_forward
mech_solution => grid_mech_spectral_polarisation_solution mechanical_solution => grid_mechanical_spectral_polarisation_solution
mech_updateCoords => grid_mech_spectral_polarisation_updateCoords mechanical_updateCoords => grid_mechanical_spectral_polarisation_updateCoords
mech_restartWrite => grid_mech_spectral_polarisation_restartWrite mechanical_restartWrite => grid_mechanical_spectral_polarisation_restartWrite
case ('FEM') case ('FEM')
mech_init => grid_mech_FEM_init mechanical_init => grid_mechanical_FEM_init
mech_forward => grid_mech_FEM_forward mechanical_forward => grid_mechanical_FEM_forward
mech_solution => grid_mech_FEM_solution mechanical_solution => grid_mechanical_FEM_solution
mech_updateCoords => grid_mech_FEM_updateCoords mechanical_updateCoords => grid_mechanical_FEM_updateCoords
mech_restartWrite => grid_mech_FEM_restartWrite mechanical_restartWrite => grid_mechanical_FEM_restartWrite
case default case default
call IO_error(error_ID = 891, ext_msg = trim(num_grid%get_asString('solver'))) call IO_error(error_ID = 891, ext_msg = trim(num_grid%get_asString('solver')))
@ -303,7 +303,7 @@ program DAMASK_grid
do field = 1, nActiveFields do field = 1, nActiveFields
select case (loadCases(1)%ID(field)) select case (loadCases(1)%ID(field))
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
call mech_init call mechanical_init
case(FIELD_THERMAL_ID) case(FIELD_THERMAL_ID)
call grid_thermal_spectral_init call grid_thermal_spectral_init
@ -379,7 +379,7 @@ program DAMASK_grid
do field = 1, nActiveFields do field = 1, nActiveFields
select case(loadCases(l)%ID(field)) select case(loadCases(l)%ID(field))
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
call mech_forward (& call mechanical_forward (&
cutBack,guess,timeinc,timeIncOld,remainingLoadCaseTime, & cutBack,guess,timeinc,timeIncOld,remainingLoadCaseTime, &
deformation_BC = loadCases(l)%deformation, & deformation_BC = loadCases(l)%deformation, &
stress_BC = loadCases(l)%stress, & stress_BC = loadCases(l)%stress, &
@ -399,7 +399,7 @@ program DAMASK_grid
do field = 1, nActiveFields do field = 1, nActiveFields
select case(loadCases(l)%ID(field)) select case(loadCases(l)%ID(field))
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
solres(field) = mech_solution(incInfo) solres(field) = mechanical_solution(incInfo)
case(FIELD_THERMAL_ID) case(FIELD_THERMAL_ID)
solres(field) = grid_thermal_spectral_solution(timeinc) solres(field) = grid_thermal_spectral_solution(timeinc)
case(FIELD_DAMAGE_ID) case(FIELD_DAMAGE_ID)
@ -420,7 +420,7 @@ program DAMASK_grid
if ( (all(solres(:)%converged .and. solres(:)%stagConverged)) & ! converged if ( (all(solres(:)%converged .and. solres(:)%stagConverged)) & ! converged
.and. .not. solres(1)%termIll) then ! and acceptable solution found .and. .not. solres(1)%termIll) then ! and acceptable solution found
call mech_updateCoords call mechanical_updateCoords
timeIncOld = timeinc timeIncOld = timeinc
cutBack = .false. cutBack = .false.
guess = .true. ! start guessing after first converged (sub)inc guess = .true. ! start guessing after first converged (sub)inc
@ -463,7 +463,7 @@ program DAMASK_grid
call MPI_Allreduce(interface_SIGUSR2,signal,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(interface_SIGUSR2,signal,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
if (ierr /= 0) error stop 'MPI error' if (ierr /= 0) error stop 'MPI error'
if (mod(inc,loadCases(l)%f_restart) == 0 .or. signal) then if (mod(inc,loadCases(l)%f_restart) == 0 .or. signal) then
call mech_restartWrite call mechanical_restartWrite
call CPFEM_restartWrite call CPFEM_restartWrite
endif endif
if(signal) call interface_setSIGUSR2(.false.) if(signal) call interface_setSIGUSR2(.false.)

135
src/grid/grid_mech_FEM.f90
View File

@ -4,7 +4,7 @@
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Grid solver for mechanics: FEM !> @brief Grid solver for mechanics: FEM
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module grid_mech_FEM module grid_mechanical_FEM
#include <petsc/finclude/petscsnes.h> #include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h> #include <petsc/finclude/petscdmda.h>
use PETScdmda use PETScdmda
@ -45,8 +45,8 @@ module grid_mech_FEM
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
DM :: mech_grid DM :: mechanical_grid
SNES :: mech_snes SNES :: mechanical_snes
Vec :: solution_current, solution_lastInc, solution_rate Vec :: solution_current, solution_lastInc, solution_rate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -79,18 +79,18 @@ module grid_mech_FEM
totalIter = 0 !< total iteration in current increment totalIter = 0 !< total iteration in current increment
public :: & public :: &
grid_mech_FEM_init, & grid_mechanical_FEM_init, &
grid_mech_FEM_solution, & grid_mechanical_FEM_solution, &
grid_mech_FEM_forward, & grid_mechanical_FEM_forward, &
grid_mech_FEM_updateCoords, & grid_mechanical_FEM_updateCoords, &
grid_mech_FEM_restartWrite grid_mechanical_FEM_restartWrite
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields and fills them with data, potentially from restart info !> @brief allocates all necessary fields and fills them with data, potentially from restart info
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_FEM_init subroutine grid_mechanical_FEM_init
real(pReal), parameter :: HGCoeff = 0.0e-2_pReal real(pReal), parameter :: HGCoeff = 0.0e-2_pReal
real(pReal), parameter, dimension(4,8) :: & real(pReal), parameter, dimension(4,8) :: &
@ -114,7 +114,7 @@ subroutine grid_mech_FEM_init
num_grid, & num_grid, &
debug_grid debug_grid
print'(/,a)', ' <<<+- grid_mech_FEM init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- grid_mechanical_FEM init -+>>>'; flush(IO_STDOUT)
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! debugging options ! debugging options
@ -141,8 +141,11 @@ subroutine grid_mech_FEM_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc ! set default and user defined options for PETSc
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type newtonls -mech_ksp_type fgmres & call PetscOptionsInsertString(PETSC_NULL_OPTIONS, &
&-mech_ksp_max_it 25 -mech_pc_type ml -mech_mg_levels_ksp_type chebyshev',ierr) '-mechanical_snes_type newtonls -mechanical_ksp_type fgmres &
&-mechanical_ksp_max_it 25 -mechanical_pc_type ml &
&-mechanical_mg_levels_ksp_type chebyshev', &
ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
@ -155,8 +158,10 @@ subroutine grid_mech_FEM_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,mech_snes,ierr); CHKERRQ(ierr) call SNESCreate(PETSC_COMM_WORLD,mechanical_snes,ierr)
call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetOptionsPrefix(mechanical_snes,'mechanical_',ierr)
CHKERRQ(ierr)
localK = 0 localK = 0
localK(worldrank) = grid3 localK(worldrank) = grid3
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
@ -167,34 +172,44 @@ subroutine grid_mech_FEM_init
1, 1, worldsize, & 1, 1, worldsize, &
3, 1, & 3, 1, &
[grid(1)],[grid(2)],localK, & [grid(1)],[grid(2)],localK, &
mech_grid,ierr) mechanical_grid,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetDM(mech_snes,mech_grid,ierr); CHKERRQ(ierr) call SNESSetDM(mechanical_snes,mechanical_grid,ierr)
call DMsetFromOptions(mech_grid,ierr); CHKERRQ(ierr)
call DMsetUp(mech_grid,ierr); CHKERRQ(ierr)
call DMDASetUniformCoordinates(mech_grid,0.0_pReal,geomSize(1),0.0_pReal,geomSize(2),0.0_pReal,geomSize(3),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call DMCreateGlobalVector(mech_grid,solution_current,ierr); CHKERRQ(ierr) call DMsetFromOptions(mechanical_grid,ierr)
call DMCreateGlobalVector(mech_grid,solution_lastInc,ierr); CHKERRQ(ierr)
call DMCreateGlobalVector(mech_grid,solution_rate ,ierr); CHKERRQ(ierr)
call DMSNESSetFunctionLocal(mech_grid,formResidual,PETSC_NULL_SNES,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call DMSNESSetJacobianLocal(mech_grid,formJacobian,PETSC_NULL_SNES,ierr) call DMsetUp(mechanical_grid,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetConvergenceTest(mech_snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged" call DMDASetUniformCoordinates(mechanical_grid,0.0_pReal,geomSize(1),0.0_pReal,geomSize(2),0.0_pReal,geomSize(3),ierr)
CHKERRQ(ierr)
call DMCreateGlobalVector(mechanical_grid,solution_current,ierr)
CHKERRQ(ierr)
call DMCreateGlobalVector(mechanical_grid,solution_lastInc,ierr)
CHKERRQ(ierr)
call DMCreateGlobalVector(mechanical_grid,solution_rate ,ierr)
CHKERRQ(ierr)
call DMSNESSetFunctionLocal(mechanical_grid,formResidual,PETSC_NULL_SNES,ierr)
CHKERRQ(ierr)
call DMSNESSetJacobianLocal(mechanical_grid,formJacobian,PETSC_NULL_SNES,ierr)
CHKERRQ(ierr)
call SNESSetConvergenceTest(mechanical_snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged"
CHKERRQ(ierr)
call SNESSetMaxLinearSolveFailures(mechanical_snes, huge(1), ierr) ! ignore linear solve failures
CHKERRQ(ierr)
call SNESSetFromOptions(mechanical_snes,ierr) ! pull it all together with additional cli arguments
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr) ! ignore linear solve failures
call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
call VecSet(solution_current,0.0_pReal,ierr);CHKERRQ(ierr) call VecSet(solution_current,0.0_pReal,ierr);CHKERRQ(ierr)
call VecSet(solution_lastInc,0.0_pReal,ierr);CHKERRQ(ierr) call VecSet(solution_lastInc,0.0_pReal,ierr);CHKERRQ(ierr)
call VecSet(solution_rate ,0.0_pReal,ierr);CHKERRQ(ierr) call VecSet(solution_rate ,0.0_pReal,ierr);CHKERRQ(ierr)
call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(mechanical_grid,solution_current,u_current,ierr)
call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) CHKERRQ(ierr)
call DMDAVecGetArrayF90(mechanical_grid,solution_lastInc,u_lastInc,ierr)
CHKERRQ(ierr)
call DMDAGetCorners(mech_grid,xstart,ystart,zstart,xend,yend,zend,ierr) ! local grid extent call DMDAGetCorners(mechanical_grid,xstart,ystart,zstart,xend,yend,zend,ierr) ! local grid extent
CHKERRQ(ierr) CHKERRQ(ierr)
xend = xstart+xend-1 xend = xstart+xend-1
yend = ystart+yend-1 yend = ystart+yend-1
@ -242,9 +257,9 @@ subroutine grid_mech_FEM_init
call utilities_constitutiveResponse(P_current,P_av,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 call utilities_constitutiveResponse(P_current,P_av,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2
F, & ! target F F, & ! target F
0.0_pReal) ! time increment 0.0_pReal) ! time increment
call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr) call DMDAVecRestoreArrayF90(mechanical_grid,solution_current,u_current,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr) call DMDAVecRestoreArrayF90(mechanical_grid,solution_lastInc,u_lastInc,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
restartRead2: if (interface_restartInc > 0) then restartRead2: if (interface_restartInc > 0) then
@ -257,13 +272,13 @@ subroutine grid_mech_FEM_init
endif restartRead2 endif restartRead2
end subroutine grid_mech_FEM_init end subroutine grid_mechanical_FEM_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the FEM scheme with internal iterations !> @brief solution for the FEM scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function grid_mech_FEM_solution(incInfoIn) result(solution) function grid_mechanical_FEM_solution(incInfoIn) result(solution)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! input data for solution ! input data for solution
@ -284,11 +299,13 @@ function grid_mech_FEM_solution(incInfoIn) result(solution)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! solve BVP ! solve BVP
call SNESsolve(mech_snes,PETSC_NULL_VEC,solution_current,ierr); CHKERRQ(ierr) call SNESsolve(mechanical_snes,PETSC_NULL_VEC,solution_current,ierr)
CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! check convergence ! check convergence
call SNESGetConvergedReason(mech_snes,reason,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(mechanical_snes,reason,ierr)
CHKERRQ(ierr)
solution%converged = reason > 0 solution%converged = reason > 0
solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
@ -296,14 +313,14 @@ function grid_mech_FEM_solution(incInfoIn) result(solution)
terminallyIll = .false. terminallyIll = .false.
P_aim = merge(P_aim,P_av,params%stress_mask) P_aim = merge(P_aim,P_av,params%stress_mask)
end function grid_mech_FEM_solution end function grid_mechanical_FEM_solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forwarding routine !> @brief forwarding routine
!> @details find new boundary conditions and best F estimate for end of current timestep !> @details find new boundary conditions and best F estimate for end of current timestep
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,& subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,&
deformation_BC,stress_BC,rotation_BC) deformation_BC,stress_BC,rotation_BC)
logical, intent(in) :: & logical, intent(in) :: &
@ -323,8 +340,10 @@ subroutine grid_mech_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,&
u_current,u_lastInc u_current,u_lastInc
call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(mechanical_grid,solution_current,u_current,ierr)
call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) CHKERRQ(ierr)
call DMDAVecGetArrayF90(mechanical_grid,solution_lastInc,u_lastInc,ierr)
CHKERRQ(ierr)
if (cutBack) then if (cutBack) then
C_volAvg = C_volAvgLastInc C_volAvg = C_volAvgLastInc
@ -371,8 +390,10 @@ subroutine grid_mech_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,&
call VecAXPY(solution_current,Delta_t,solution_rate,ierr); CHKERRQ(ierr) call VecAXPY(solution_current,Delta_t,solution_rate,ierr); CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(mechanical_grid,solution_current,u_current,ierr)
call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(mechanical_grid,solution_lastInc,u_lastInc,ierr)
CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide available data ! set module wide available data
@ -380,31 +401,33 @@ subroutine grid_mech_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,&
params%rotation_BC = rotation_BC params%rotation_BC = rotation_BC
params%timeinc = Delta_t params%timeinc = Delta_t
end subroutine grid_mech_FEM_forward end subroutine grid_mechanical_FEM_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Update coordinates !> @brief Update coordinates
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_FEM_updateCoords subroutine grid_mechanical_FEM_updateCoords
call utilities_updateCoords(F) call utilities_updateCoords(F)
end subroutine grid_mech_FEM_updateCoords end subroutine grid_mechanical_FEM_updateCoords
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Write current solver and constitutive data for restart to file !> @brief Write current solver and constitutive data for restart to file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_FEM_restartWrite subroutine grid_mechanical_FEM_restartWrite
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
PetscScalar, dimension(:,:,:,:), pointer :: u_current,u_lastInc PetscScalar, dimension(:,:,:,:), pointer :: u_current,u_lastInc
character(len=pStringLen) :: fileName character(len=pStringLen) :: fileName
call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(mechanical_grid,solution_current,u_current,ierr)
call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) CHKERRQ(ierr)
call DMDAVecGetArrayF90(mechanical_grid,solution_lastInc,u_lastInc,ierr)
CHKERRQ(ierr)
print*, 'writing solver data required for restart to file'; flush(IO_STDOUT) print*, 'writing solver data required for restart to file'; flush(IO_STDOUT)
@ -427,10 +450,12 @@ subroutine grid_mech_FEM_restartWrite
call HDF5_closeGroup(groupHandle) call HDF5_closeGroup(groupHandle)
call HDF5_closeFile(fileHandle) call HDF5_closeFile(fileHandle)
call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr);CHKERRQ(ierr) call DMDAVecRestoreArrayF90(mechanical_grid,solution_current,u_current,ierr)
call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr);CHKERRQ(ierr) CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(mechanical_grid,solution_lastInc,u_lastInc,ierr)
CHKERRQ(ierr)
end subroutine grid_mech_FEM_restartWrite end subroutine grid_mechanical_FEM_restartWrite
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -498,8 +523,10 @@ subroutine formResidual(da_local,x_local, &
PetscErrorCode :: ierr PetscErrorCode :: ierr
real(pReal), dimension(3,3,3,3) :: devNull real(pReal), dimension(3,3,3,3) :: devNull
call SNESGetNumberFunctionEvals(mech_snes,nfuncs,ierr); CHKERRQ(ierr) call SNESGetNumberFunctionEvals(mechanical_snes,nfuncs,ierr)
call SNESGetIterationNumber(mech_snes,PETScIter,ierr); CHKERRQ(ierr) CHKERRQ(ierr)
call SNESGetIterationNumber(mechanical_snes,PETScIter,ierr)
CHKERRQ(ierr)
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
@ -679,4 +706,4 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr)
end subroutine formJacobian end subroutine formJacobian
end module grid_mech_FEM end module grid_mechanical_FEM

40
src/grid/grid_mech_spectral_basic.f90
View File

@ -4,7 +4,7 @@
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Grid solver for mechanics: Spectral basic !> @brief Grid solver for mechanics: Spectral basic
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module grid_mech_spectral_basic module grid_mechanical_spectral_basic
#include <petsc/finclude/petscsnes.h> #include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h> #include <petsc/finclude/petscdmda.h>
use PETScdmda use PETScdmda
@ -79,18 +79,18 @@ module grid_mech_spectral_basic
totalIter = 0 !< total iteration in current increment totalIter = 0 !< total iteration in current increment
public :: & public :: &
grid_mech_spectral_basic_init, & grid_mechanical_spectral_basic_init, &
grid_mech_spectral_basic_solution, & grid_mechanical_spectral_basic_solution, &
grid_mech_spectral_basic_forward, & grid_mechanical_spectral_basic_forward, &
grid_mech_spectral_basic_updateCoords, & grid_mechanical_spectral_basic_updateCoords, &
grid_mech_spectral_basic_restartWrite grid_mechanical_spectral_basic_restartWrite
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields and fills them with data, potentially from restart info !> @brief allocates all necessary fields and fills them with data, potentially from restart info
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_init subroutine grid_mechanical_spectral_basic_init
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
PetscErrorCode :: ierr PetscErrorCode :: ierr
@ -105,7 +105,7 @@ subroutine grid_mech_spectral_basic_init
num_grid, & num_grid, &
debug_grid debug_grid
print'(/,a)', ' <<<+- grid_mech_spectral_basic init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- grid_mechanical_spectral_basic init -+>>>'; flush(IO_STDOUT)
print*, 'Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013' print*, 'Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013'
print*, 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL print*, 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL
@ -139,7 +139,7 @@ subroutine grid_mech_spectral_basic_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc ! set default and user defined options for PETSc
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mechanical_snes_type ngmres',ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
@ -152,7 +152,7 @@ subroutine grid_mech_spectral_basic_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) call SNESSetOptionsPrefix(snes,'mechanical_',ierr);CHKERRQ(ierr)
localK = 0 localK = 0
localK(worldrank) = grid3 localK(worldrank) = grid3
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
@ -222,13 +222,13 @@ subroutine grid_mech_spectral_basic_init
call utilities_updateGamma(C_minMaxAvg) call utilities_updateGamma(C_minMaxAvg)
call utilities_saveReferenceStiffness call utilities_saveReferenceStiffness
end subroutine grid_mech_spectral_basic_init end subroutine grid_mechanical_spectral_basic_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the basic scheme with internal iterations !> @brief solution for the basic scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function grid_mech_spectral_basic_solution(incInfoIn) result(solution) function grid_mechanical_spectral_basic_solution(incInfoIn) result(solution)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! input data for solution ! input data for solution
@ -262,14 +262,14 @@ function grid_mech_spectral_basic_solution(incInfoIn) result(solution)
terminallyIll = .false. terminallyIll = .false.
P_aim = merge(P_aim,P_av,params%stress_mask) P_aim = merge(P_aim,P_av,params%stress_mask)
end function grid_mech_spectral_basic_solution end function grid_mechanical_spectral_basic_solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forwarding routine !> @brief forwarding routine
!> @details find new boundary conditions and best F estimate for end of current timestep !> @details find new boundary conditions and best F estimate for end of current timestep
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,& subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,&
deformation_BC,stress_BC,rotation_BC) deformation_BC,stress_BC,rotation_BC)
logical, intent(in) :: & logical, intent(in) :: &
@ -339,13 +339,13 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_old,t_
params%rotation_BC = rotation_BC params%rotation_BC = rotation_BC
params%timeinc = Delta_t params%timeinc = Delta_t
end subroutine grid_mech_spectral_basic_forward end subroutine grid_mechanical_spectral_basic_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Update coordinates !> @brief Update coordinates
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_updateCoords subroutine grid_mechanical_spectral_basic_updateCoords
PetscErrorCode :: ierr PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: F PetscScalar, dimension(:,:,:,:), pointer :: F
@ -354,13 +354,13 @@ subroutine grid_mech_spectral_basic_updateCoords
call utilities_updateCoords(F) call utilities_updateCoords(F)
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
end subroutine grid_mech_spectral_basic_updateCoords end subroutine grid_mechanical_spectral_basic_updateCoords
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Write current solver and constitutive data for restart to file !> @brief Write current solver and constitutive data for restart to file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_restartWrite subroutine grid_mechanical_spectral_basic_restartWrite
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
@ -393,7 +393,7 @@ subroutine grid_mech_spectral_basic_restartWrite
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
end subroutine grid_mech_spectral_basic_restartWrite end subroutine grid_mechanical_spectral_basic_restartWrite
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -506,4 +506,4 @@ subroutine formResidual(in, F, &
end subroutine formResidual end subroutine formResidual
end module grid_mech_spectral_basic end module grid_mechanical_spectral_basic

40
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -4,7 +4,7 @@
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Grid solver for mechanics: Spectral Polarisation !> @brief Grid solver for mechanics: Spectral Polarisation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module grid_mech_spectral_polarisation module grid_mechanical_spectral_polarisation
#include <petsc/finclude/petscsnes.h> #include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h> #include <petsc/finclude/petscdmda.h>
use PETScdmda use PETScdmda
@ -90,18 +90,18 @@ module grid_mech_spectral_polarisation
totalIter = 0 !< total iteration in current increment totalIter = 0 !< total iteration in current increment
public :: & public :: &
grid_mech_spectral_polarisation_init, & grid_mechanical_spectral_polarisation_init, &
grid_mech_spectral_polarisation_solution, & grid_mechanical_spectral_polarisation_solution, &
grid_mech_spectral_polarisation_forward, & grid_mechanical_spectral_polarisation_forward, &
grid_mech_spectral_polarisation_updateCoords, & grid_mechanical_spectral_polarisation_updateCoords, &
grid_mech_spectral_polarisation_restartWrite grid_mechanical_spectral_polarisation_restartWrite
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields and fills them with data, potentially from restart info !> @brief allocates all necessary fields and fills them with data, potentially from restart info
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_polarisation_init subroutine grid_mechanical_spectral_polarisation_init
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
PetscErrorCode :: ierr PetscErrorCode :: ierr
@ -118,7 +118,7 @@ subroutine grid_mech_spectral_polarisation_init
num_grid, & num_grid, &
debug_grid debug_grid
print'(/,a)', ' <<<+- grid_mech_spectral_polarisation init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- grid_mechanical_spectral_polarisation init -+>>>'; flush(IO_STDOUT)
print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006' print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006'
@ -157,7 +157,7 @@ subroutine grid_mech_spectral_polarisation_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc ! set default and user defined options for PETSc
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mechanical_snes_type ngmres',ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
@ -172,7 +172,7 @@ subroutine grid_mech_spectral_polarisation_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) call SNESSetOptionsPrefix(snes,'mechanical_',ierr);CHKERRQ(ierr)
localK = 0 localK = 0
localK(worldrank) = grid3 localK(worldrank) = grid3
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
@ -250,13 +250,13 @@ subroutine grid_mech_spectral_polarisation_init
C_scale = C_minMaxAvg C_scale = C_minMaxAvg
S_scale = math_invSym3333(C_minMaxAvg) S_scale = math_invSym3333(C_minMaxAvg)
end subroutine grid_mech_spectral_polarisation_init end subroutine grid_mechanical_spectral_polarisation_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the Polarisation scheme with internal iterations !> @brief solution for the Polarisation scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function grid_mech_spectral_polarisation_solution(incInfoIn) result(solution) function grid_mechanical_spectral_polarisation_solution(incInfoIn) result(solution)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! input data for solution ! input data for solution
@ -294,14 +294,14 @@ function grid_mech_spectral_polarisation_solution(incInfoIn) result(solution)
terminallyIll = .false. terminallyIll = .false.
P_aim = merge(P_aim,P_av,params%stress_mask) P_aim = merge(P_aim,P_av,params%stress_mask)
end function grid_mech_spectral_polarisation_solution end function grid_mechanical_spectral_polarisation_solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forwarding routine !> @brief forwarding routine
!> @details find new boundary conditions and best F estimate for end of current timestep !> @details find new boundary conditions and best F estimate for end of current timestep
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,& subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,&
deformation_BC,stress_BC,rotation_BC) deformation_BC,stress_BC,rotation_BC)
logical, intent(in) :: & logical, intent(in) :: &
@ -393,13 +393,13 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,Delta_t,Delta_t
params%rotation_BC = rotation_BC params%rotation_BC = rotation_BC
params%timeinc = Delta_t params%timeinc = Delta_t
end subroutine grid_mech_spectral_polarisation_forward end subroutine grid_mechanical_spectral_polarisation_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Update coordinates !> @brief Update coordinates
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_polarisation_updateCoords subroutine grid_mechanical_spectral_polarisation_updateCoords
PetscErrorCode :: ierr PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau
@ -408,13 +408,13 @@ subroutine grid_mech_spectral_polarisation_updateCoords
call utilities_updateCoords(FandF_tau(0:8,:,:,:)) call utilities_updateCoords(FandF_tau(0:8,:,:,:))
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
end subroutine grid_mech_spectral_polarisation_updateCoords end subroutine grid_mechanical_spectral_polarisation_updateCoords
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Write current solver and constitutive data for restart to file !> @brief Write current solver and constitutive data for restart to file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_polarisation_restartWrite subroutine grid_mechanical_spectral_polarisation_restartWrite
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
@ -450,7 +450,7 @@ subroutine grid_mech_spectral_polarisation_restartWrite
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
end subroutine grid_mech_spectral_polarisation_restartWrite end subroutine grid_mechanical_spectral_polarisation_restartWrite
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -618,4 +618,4 @@ subroutine formResidual(in, FandF_tau, &
end subroutine formResidual end subroutine formResidual
end module grid_mech_spectral_polarisation end module grid_mechanical_spectral_polarisation

32
src/homogenization.f90
View File

@ -45,10 +45,10 @@ module homogenization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
interface interface
module subroutine mech_init(num_homog) module subroutine mechanical_init(num_homog)
class(tNode), pointer, intent(in) :: & class(tNode), pointer, intent(in) :: &
num_homog !< pointer to mechanical homogenization numerics data num_homog !< pointer to mechanical homogenization numerics data
end subroutine mech_init end subroutine mechanical_init
module subroutine thermal_init module subroutine thermal_init
end subroutine thermal_init end subroutine thermal_init
@ -56,13 +56,13 @@ module homogenization
module subroutine damage_init module subroutine damage_init
end subroutine damage_init end subroutine damage_init
module subroutine mech_partition(subF,ip,el) module subroutine mechanical_partition(subF,ip,el)
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
subF subF
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point ip, & !< integration point
el !< element number el !< element number
end subroutine mech_partition end subroutine mechanical_partition
module subroutine thermal_partition(ce) module subroutine thermal_partition(ce)
integer, intent(in) :: ce integer, intent(in) :: ce
@ -76,19 +76,19 @@ module homogenization
integer, intent(in) :: ip,el integer, intent(in) :: ip,el
end subroutine thermal_homogenize end subroutine thermal_homogenize
module subroutine mech_homogenize(dt,ip,el) module subroutine mechanical_homogenize(dt,ip,el)
real(pReal), intent(in) :: dt real(pReal), intent(in) :: dt
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point ip, & !< integration point
el !< element number el !< element number
end subroutine mech_homogenize end subroutine mechanical_homogenize
module subroutine mech_results(group_base,h) module subroutine mechanical_results(group_base,h)
character(len=*), intent(in) :: group_base character(len=*), intent(in) :: group_base
integer, intent(in) :: h integer, intent(in) :: h
end subroutine mech_results end subroutine mechanical_results
module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy) module function mechanical_updateState(subdt,subF,ip,el) result(doneAndHappy)
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
subdt !< current time step subdt !< current time step
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
@ -97,7 +97,7 @@ module homogenization
ip, & !< integration point ip, & !< integration point
el !< element number el !< element number
logical, dimension(2) :: doneAndHappy logical, dimension(2) :: doneAndHappy
end function mech_updateState end function mechanical_updateState
module function thermal_conduction_getConductivity(ip,el) result(K) module function thermal_conduction_getConductivity(ip,el) result(K)
@ -216,7 +216,7 @@ subroutine homogenization_init()
if (num%nMPstate < 1) call IO_error(301,ext_msg='nMPstate') if (num%nMPstate < 1) call IO_error(301,ext_msg='nMPstate')
call mech_init(num_homog) call mechanical_init(num_homog)
call thermal_init() call thermal_init()
call damage_init() call damage_init()
@ -253,7 +253,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
ce = (el-1)*discretization_nIPs + ip ce = (el-1)*discretization_nIPs + ip
me = material_homogenizationMemberAt2(ce) me = material_homogenizationMemberAt2(ce)
call constitutive_restore(ce,.false.) ! wrong name (is more a forward function) call phase_restore(ce,.false.) ! wrong name (is more a forward function)
if(homogState(ho)%sizeState > 0) homogState(ho)%State(:,me) = homogState(ho)%State0(:,me) if(homogState(ho)%sizeState > 0) homogState(ho)%State(:,me) = homogState(ho)%State0(:,me)
if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%State(:,me) = damageState_h(ho)%State0(:,me) if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%State(:,me) = damageState_h(ho)%State0(:,me)
@ -267,7 +267,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
if (.not. doneAndHappy(1)) then if (.not. doneAndHappy(1)) then
call mech_partition(homogenization_F(1:3,1:3,ce),ip,el) call mechanical_partition(homogenization_F(1:3,1:3,ce),ip,el)
converged = .true. converged = .true.
do co = 1, myNgrains do co = 1, myNgrains
converged = converged .and. crystallite_stress(dt,co,ip,el) converged = converged .and. crystallite_stress(dt,co,ip,el)
@ -276,7 +276,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
if (.not. converged) then if (.not. converged) then
doneAndHappy = [.true.,.false.] doneAndHappy = [.true.,.false.]
else else
doneAndHappy = mech_updateState(dt,homogenization_F(1:3,1:3,ce),ip,el) doneAndHappy = mechanical_updateState(dt,homogenization_F(1:3,1:3,ce),ip,el)
converged = all(doneAndHappy) converged = all(doneAndHappy)
endif endif
endif endif
@ -338,7 +338,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
do co = 1, homogenization_Nconstituents(ho) do co = 1, homogenization_Nconstituents(ho)
call crystallite_orientations(co,ip,el) call crystallite_orientations(co,ip,el)
enddo enddo
call mech_homogenize(dt,ip,el) call mechanical_homogenize(dt,ip,el)
enddo IpLooping3 enddo IpLooping3
enddo elementLooping3 enddo elementLooping3
!$OMP END DO !$OMP END DO
@ -365,7 +365,7 @@ subroutine homogenization_results
group_base = 'current/homogenization/'//trim(material_name_homogenization(ho)) group_base = 'current/homogenization/'//trim(material_name_homogenization(ho))
call results_closeGroup(results_addGroup(group_base)) call results_closeGroup(results_addGroup(group_base))
call mech_results(group_base,ho) call mechanical_results(group_base,ho)
group = trim(group_base)//'/damage' group = trim(group_base)//'/damage'
call results_closeGroup(results_addGroup(group)) call results_closeGroup(results_addGroup(group))

4
src/homogenization_damage.f90
View File

@ -74,7 +74,7 @@ module subroutine damage_partition(ce)
phi = current(material_homogenizationAt2(ce))%phi(material_homogenizationMemberAt2(ce)) phi = current(material_homogenizationAt2(ce))%phi(material_homogenizationMemberAt2(ce))
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
call constitutive_damage_set_phi(phi,co,ce) call phase_damage_set_phi(phi,co,ce)
enddo enddo
end subroutine damage_partition end subroutine damage_partition
@ -120,7 +120,7 @@ module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, p
phiDot = 0.0_pReal phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal dPhiDot_dPhi = 0.0_pReal
call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) call phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)

94
src/homogenization_mechanics.f90 → src/homogenization_mechanical.f90
View File

@ -7,52 +7,52 @@ submodule(homogenization) mechanics
interface interface
module subroutine mech_none_init module subroutine mechanical_pass_init
end subroutine mech_none_init end subroutine mechanical_pass_init
module subroutine mech_isostrain_init module subroutine mechanical_isostrain_init
end subroutine mech_isostrain_init end subroutine mechanical_isostrain_init
module subroutine mech_RGC_init(num_homogMech) module subroutine mechanical_RGC_init(num_homogMech)
class(tNode), pointer, intent(in) :: & class(tNode), pointer, intent(in) :: &
num_homogMech !< pointer to mechanical homogenization numerics data num_homogMech !< pointer to mechanical homogenization numerics data
end subroutine mech_RGC_init end subroutine mechanical_RGC_init
module subroutine mech_isostrain_partitionDeformation(F,avgF) module subroutine mechanical_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 mechanical_isostrain_partitionDeformation
module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of) module subroutine mechanical_RGC_partitionDeformation(F,avgF,instance,of)
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
integer, intent(in) :: & integer, intent(in) :: &
instance, & instance, &
of of
end subroutine mech_RGC_partitionDeformation end subroutine mechanical_RGC_partitionDeformation
module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) module subroutine mechanical_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
real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses
real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
integer, intent(in) :: instance integer, intent(in) :: instance
end subroutine mech_isostrain_averageStressAndItsTangent end subroutine mechanical_isostrain_averageStressAndItsTangent
module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) module subroutine mechanical_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
real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses
real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
integer, intent(in) :: instance integer, intent(in) :: instance
end subroutine mech_RGC_averageStressAndItsTangent end subroutine mechanical_RGC_averageStressAndItsTangent
module function mech_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy) module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy)
logical, dimension(2) :: doneAndHappy logical, dimension(2) :: doneAndHappy
real(pReal), dimension(:,:,:), intent(in) :: & real(pReal), dimension(:,:,:), intent(in) :: &
P,& !< partitioned stresses P,& !< partitioned stresses
@ -63,13 +63,13 @@ submodule(homogenization) mechanics
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
end function mech_RGC_updateState end function mechanical_RGC_updateState
module subroutine mech_RGC_results(instance,group) module subroutine mechanical_RGC_results(instance,group)
integer, intent(in) :: instance !< homogenization instance integer, intent(in) :: instance !< homogenization instance
character(len=*), intent(in) :: group !< group name in HDF5 file character(len=*), intent(in) :: group !< group name in HDF5 file
end subroutine mech_RGC_results end subroutine mechanical_RGC_results
end interface end interface
@ -78,7 +78,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Allocate variables and set parameters. !> @brief Allocate variables and set parameters.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_init(num_homog) module subroutine mechanical_init(num_homog)
class(tNode), pointer, intent(in) :: & class(tNode), pointer, intent(in) :: &
num_homog num_homog
@ -94,17 +94,17 @@ module subroutine mech_init(num_homog)
allocate(homogenization_P(3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal) allocate(homogenization_P(3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal)
num_homogMech => num_homog%get('mech',defaultVal=emptyDict) num_homogMech => num_homog%get('mech',defaultVal=emptyDict)
if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mech_none_init if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mechanical_pass_init
if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mech_isostrain_init if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mechanical_isostrain_init
if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call mech_RGC_init(num_homogMech) if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call mechanical_RGC_init(num_homogMech)
end subroutine mech_init end subroutine mechanical_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Partition F onto the individual constituents. !> @brief Partition F onto the individual constituents.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_partition(subF,ip,el) module subroutine mechanical_partition(subF,ip,el)
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
subF subF
@ -122,25 +122,25 @@ module subroutine mech_partition(subF,ip,el)
Fs(1:3,1:3,1) = subF Fs(1:3,1:3,1) = subF
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
call mech_isostrain_partitionDeformation(Fs,subF) call mechanical_isostrain_partitionDeformation(Fs,subF)
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
call mech_RGC_partitionDeformation(Fs,subF,ip,el) call mechanical_RGC_partitionDeformation(Fs,subF,ip,el)
end select chosenHomogenization end select chosenHomogenization
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
call constitutive_mech_setF(Fs(1:3,1:3,co),co,ip,el) call phase_mechanical_setF(Fs(1:3,1:3,co),co,ip,el)
enddo enddo
end subroutine mech_partition end subroutine mechanical_partition
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Average P and dPdF from the individual constituents. !> @brief Average P and dPdF from the individual constituents.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_homogenize(dt,ip,el) module subroutine mechanical_homogenize(dt,ip,el)
real(pReal), intent(in) :: dt real(pReal), intent(in) :: dt
integer, intent(in) :: & integer, intent(in) :: &
@ -156,15 +156,15 @@ module subroutine mech_homogenize(dt,ip,el)
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization case (HOMOGENIZATION_NONE_ID) chosenHomogenization
homogenization_P(1:3,1:3,ce) = constitutive_mech_getP(1,ip,el) homogenization_P(1:3,1:3,ce) = phase_mechanical_getP(1,ip,el)
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = constitutive_mech_dPdF(dt,1,ip,el) homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ip,el)
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(dt,co,ip,el) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ip,el)
Ps(:,:,co) = constitutive_mech_getP(co,ip,el) Ps(:,:,co) = phase_mechanical_getP(co,ip,el)
enddo enddo
call mech_isostrain_averageStressAndItsTangent(& call mechanical_isostrain_averageStressAndItsTangent(&
homogenization_P(1:3,1:3,ce), & homogenization_P(1:3,1:3,ce), &
homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
Ps,dPdFs, & Ps,dPdFs, &
@ -172,10 +172,10 @@ module subroutine mech_homogenize(dt,ip,el)
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(dt,co,ip,el) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ip,el)
Ps(:,:,co) = constitutive_mech_getP(co,ip,el) Ps(:,:,co) = phase_mechanical_getP(co,ip,el)
enddo enddo
call mech_RGC_averageStressAndItsTangent(& call mechanical_RGC_averageStressAndItsTangent(&
homogenization_P(1:3,1:3,ce), & homogenization_P(1:3,1:3,ce), &
homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
Ps,dPdFs, & Ps,dPdFs, &
@ -183,14 +183,14 @@ module subroutine mech_homogenize(dt,ip,el)
end select chosenHomogenization end select chosenHomogenization
end subroutine mech_homogenize end subroutine mechanical_homogenize
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @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
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy) module function mechanical_updateState(subdt,subF,ip,el) result(doneAndHappy)
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
subdt !< current time step subdt !< current time step
@ -209,22 +209,22 @@ module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy)
if (homogenization_type(material_homogenizationAt(el)) == HOMOGENIZATION_RGC_ID) then if (homogenization_type(material_homogenizationAt(el)) == HOMOGENIZATION_RGC_ID) then
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(subdt,co,ip,el) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(subdt,co,ip,el)
Fs(:,:,co) = constitutive_mech_getF(co,ip,el) Fs(:,:,co) = phase_mechanical_getF(co,ip,el)
Ps(:,:,co) = constitutive_mech_getP(co,ip,el) Ps(:,:,co) = phase_mechanical_getP(co,ip,el)
enddo enddo
doneAndHappy = mech_RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ip,el) doneAndHappy = mechanical_RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ip,el)
else else
doneAndHappy = .true. doneAndHappy = .true.
endif endif
end function mech_updateState end function mechanical_updateState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Write results to file. !> @brief Write results to file.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_results(group_base,h) module subroutine mechanical_results(group_base,h)
use material, only: & use material, only: &
material_homogenization_type => homogenization_type material_homogenization_type => homogenization_type
@ -239,7 +239,7 @@ module subroutine mech_results(group_base,h)
select case(material_homogenization_type(h)) select case(material_homogenization_type(h))
case(HOMOGENIZATION_rgc_ID) case(HOMOGENIZATION_rgc_ID)
call mech_RGC_results(homogenization_typeInstance(h),group) call mechanical_RGC_results(homogenization_typeInstance(h),group)
end select end select
@ -250,7 +250,7 @@ module subroutine mech_results(group_base,h)
!call results_writeDataset(group,temp,'P',& !call results_writeDataset(group,temp,'P',&
! '1st Piola-Kirchhoff stress','Pa') ! '1st Piola-Kirchhoff stress','Pa')
end subroutine mech_results end subroutine mechanical_results
end submodule mechanics end submodule mechanics

24
src/homogenization_mechanics_RGC.f90 → src/homogenization_mechanical_RGC.f90
View File

@ -71,7 +71,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields, reads information from material configuration file !> @brief allocates all necessary fields, reads information from material configuration file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_RGC_init(num_homogMech) module subroutine mechanical_RGC_init(num_homogMech)
class(tNode), pointer, intent(in) :: & class(tNode), pointer, intent(in) :: &
num_homogMech !< pointer to mechanical homogenization numerics data num_homogMech !< pointer to mechanical homogenization numerics data
@ -155,7 +155,7 @@ module subroutine mech_RGC_init(num_homogMech)
prm%N_constituents = homogMech%get_asInts('cluster_size',requiredSize=3) prm%N_constituents = homogMech%get_asInts('cluster_size',requiredSize=3)
if (homogenization_Nconstituents(h) /= product(prm%N_constituents)) & if (homogenization_Nconstituents(h) /= product(prm%N_constituents)) &
call IO_error(211,ext_msg='N_constituents (mech_RGC)') call IO_error(211,ext_msg='N_constituents (mechanical_RGC)')
prm%xi_alpha = homogMech%get_asFloat('xi_alpha') prm%xi_alpha = homogMech%get_asFloat('xi_alpha')
prm%c_alpha = homogMech%get_asFloat('c_alpha') prm%c_alpha = homogMech%get_asFloat('c_alpha')
@ -190,13 +190,13 @@ module subroutine mech_RGC_init(num_homogMech)
enddo enddo
end subroutine mech_RGC_init end subroutine mechanical_RGC_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief partitions the deformation gradient onto the constituents !> @brief partitions the deformation gradient onto the constituents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of) module subroutine mechanical_RGC_partitionDeformation(F,avgF,instance,of)
real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned F per grain real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned F per grain
@ -229,14 +229,14 @@ module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
end associate end associate
end subroutine mech_RGC_partitionDeformation end subroutine mechanical_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
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function mech_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy) module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy)
logical, dimension(2) :: doneAndHappy logical, dimension(2) :: doneAndHappy
real(pReal), dimension(:,:,:), intent(in) :: & real(pReal), dimension(:,:,:), intent(in) :: &
P,& !< partitioned stresses P,& !< partitioned stresses
@ -658,7 +658,7 @@ module function mech_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy
real(pReal), dimension(6,6) :: C real(pReal), dimension(6,6) :: C
C = constitutive_homogenizedC(material_phaseAt(grainID,el),material_phaseMemberAt(grainID,ip,el)) C = phase_homogenizedC(material_phaseAt(grainID,el),material_phaseMemberAt(grainID,ip,el))
equivalentMu = lattice_equivalent_mu(C,'voigt') equivalentMu = lattice_equivalent_mu(C,'voigt')
end function equivalentMu end function equivalentMu
@ -704,13 +704,13 @@ module function mech_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy
end subroutine grainDeformation end subroutine grainDeformation
end function mech_RGC_updateState end function mechanical_RGC_updateState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief derive average stress and stiffness from constituent quantities !> @brief derive average stress and stiffness from constituent quantities
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) module subroutine mechanical_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
@ -722,13 +722,13 @@ module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ins
avgP = sum(P,3) /real(product(param(instance)%N_constituents),pReal) avgP = sum(P,3) /real(product(param(instance)%N_constituents),pReal)
dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%N_constituents),pReal) dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%N_constituents),pReal)
end subroutine mech_RGC_averageStressAndItsTangent end subroutine mechanical_RGC_averageStressAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_RGC_results(instance,group) module subroutine mechanical_RGC_results(instance,group)
integer, intent(in) :: instance integer, intent(in) :: instance
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
@ -754,7 +754,7 @@ module subroutine mech_RGC_results(instance,group)
enddo outputsLoop enddo outputsLoop
end associate end associate
end subroutine mech_RGC_results end subroutine mechanical_RGC_results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

14
src/homogenization_mechanics_isostrain.f90 → src/homogenization_mechanical_isostrain.f90
View File

@ -26,7 +26,7 @@ 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 mechanical_isostrain_init
integer :: & integer :: &
Ninstances, & Ninstances, &
@ -58,7 +58,7 @@ module subroutine mech_isostrain_init
case ('avg') case ('avg')
prm%mapping = average_ID prm%mapping = average_ID
case default case default
call IO_error(211,ext_msg='sum'//' (mech_isostrain)') call IO_error(211,ext_msg='sum'//' (mechanical_isostrain)')
end select end select
Nmaterialpoints = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == h)
@ -70,13 +70,13 @@ module subroutine mech_isostrain_init
enddo enddo
end subroutine mech_isostrain_init end subroutine mechanical_isostrain_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief partitions the deformation gradient onto the constituents !> @brief partitions the deformation gradient onto the constituents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_isostrain_partitionDeformation(F,avgF) module subroutine mechanical_isostrain_partitionDeformation(F,avgF)
real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient
@ -84,13 +84,13 @@ module subroutine mech_isostrain_partitionDeformation(F,avgF)
F = spread(avgF,3,size(F,3)) F = spread(avgF,3,size(F,3))
end subroutine mech_isostrain_partitionDeformation end subroutine mechanical_isostrain_partitionDeformation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief derive average stress and stiffness from constituent quantities !> @brief derive average stress and stiffness from constituent quantities
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) module subroutine mechanical_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
@ -112,6 +112,6 @@ module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dP
end associate end associate
end subroutine mech_isostrain_averageStressAndItsTangent end subroutine mechanical_isostrain_averageStressAndItsTangent
end submodule isostrain end submodule isostrain

8
src/homogenization_mechanics_none.f90 → src/homogenization_mechanical_pass.f90
View File

@ -11,14 +11,14 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields, reads information from material configuration file !> @brief allocates all necessary fields, reads information from material configuration file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_none_init module subroutine mechanical_pass_init
integer :: & integer :: &
Ninstances, & Ninstances, &
h, & h, &
Nmaterialpoints Nmaterialpoints
print'(/,a)', ' <<<+- homogenization:mechanics:none init -+>>>' print'(/,a)', ' <<<+- homogenization:mechanics:pass init -+>>>'
Ninstances = count(homogenization_type == HOMOGENIZATION_NONE_ID) Ninstances = count(homogenization_type == HOMOGENIZATION_NONE_ID)
print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
@ -27,7 +27,7 @@ module subroutine mech_none_init
if(homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle if(homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle
if(homogenization_Nconstituents(h) /= 1) & if(homogenization_Nconstituents(h) /= 1) &
call IO_error(211,ext_msg='N_constituents (mech_none)') call IO_error(211,ext_msg='N_constituents (mechanical_pass)')
Nmaterialpoints = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == h)
homogState(h)%sizeState = 0 homogState(h)%sizeState = 0
@ -36,6 +36,6 @@ module subroutine mech_none_init
enddo enddo
end subroutine mech_none_init end subroutine mechanical_pass_init
end submodule none end submodule none

6
src/homogenization_thermal.f90
View File

@ -78,7 +78,7 @@ module subroutine thermal_partition(ce)
T = current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce)) T = current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce))
dot_T = current(material_homogenizationAt2(ce))%dot_T(material_homogenizationMemberAt2(ce)) dot_T = current(material_homogenizationAt2(ce))%dot_T(material_homogenizationMemberAt2(ce))
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
call constitutive_thermal_setField(T,dot_T,co,ce) call phase_thermal_setField(T,dot_T,co,ce)
enddo enddo
end subroutine thermal_partition end subroutine thermal_partition
@ -91,7 +91,7 @@ module subroutine thermal_homogenize(ip,el)
integer, intent(in) :: ip,el integer, intent(in) :: ip,el
!call constitutive_thermal_getRate(homogenization_dot_T((el-1)*discretization_nIPs+ip), ip,el) !call phase_thermal_getRate(homogenization_dot_T((el-1)*discretization_nIPs+ip), ip,el)
end subroutine thermal_homogenize end subroutine thermal_homogenize
@ -235,7 +235,7 @@ module subroutine thermal_conduction_getSource(Tdot, ip,el)
do co = 1, homogenization_Nconstituents(ho) do co = 1, homogenization_Nconstituents(ho)
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phasememberAt(co,ip,el) me = material_phasememberAt(co,ip,el)
call constitutive_thermal_getRate(dot_T_temp, ph,me) call phase_thermal_getRate(dot_T_temp, ph,me)
Tdot = Tdot + dot_T_temp Tdot = Tdot + dot_T_temp
enddo enddo

8
src/mesh/DAMASK_mesh.f90
View File

@ -18,7 +18,7 @@ program DAMASK_mesh
use config use config
use discretization_mesh use discretization_mesh
use FEM_Utilities use FEM_Utilities
use mesh_mech_FEM use mesh_mechanical_FEM
implicit none implicit none
@ -242,7 +242,7 @@ program DAMASK_mesh
do field = 1, nActiveFields do field = 1, nActiveFields
select case (loadCases(1)%fieldBC(field)%ID) select case (loadCases(1)%fieldBC(field)%ID)
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
call FEM_mech_init(loadCases(1)%fieldBC(field)) call FEM_mechanical_init(loadCases(1)%fieldBC(field))
end select end select
enddo enddo
@ -306,7 +306,7 @@ program DAMASK_mesh
do field = 1, nActiveFields do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID) select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
call FEM_mech_forward (& call FEM_mechanical_forward (&
guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field)) guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
end select end select
@ -320,7 +320,7 @@ program DAMASK_mesh
do field = 1, nActiveFields do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID) select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
solres(field) = FEM_mech_solution (& solres(field) = FEM_mechanical_solution (&
incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field)) incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
end select end select

12
src/mesh/FEM_utilities.f90
View File

@ -127,12 +127,12 @@ subroutine FEM_utilities_init
CHKERRQ(ierr) CHKERRQ(ierr)
if(debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr) if(debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type newtonls & call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mechanical_snes_type newtonls &
&-mech_snes_linesearch_type cp -mech_snes_ksp_ew & &-mechanical_snes_linesearch_type cp -mechanical_snes_ksp_ew &
&-mech_snes_ksp_ew_rtol0 0.01 -mech_snes_ksp_ew_rtolmax 0.01 & &-mechanical_snes_ksp_ew_rtol0 0.01 -mechanical_snes_ksp_ew_rtolmax 0.01 &
&-mech_ksp_type fgmres -mech_ksp_max_it 25 & &-mechanical_ksp_type fgmres -mechanical_ksp_max_it 25 &
&-mech_pc_type ml -mech_mg_levels_ksp_type chebyshev & &-mechanical_pc_type ml -mechanical_mg_levels_ksp_type chebyshev &
&-mech_mg_levels_pc_type sor -mech_pc_ml_nullspace user',ierr) &-mechanical_mg_levels_pc_type sor -mechanical_pc_ml_nullspace user',ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_mesh%get_asString('petsc_options',defaultVal=''),ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_mesh%get_asString('petsc_options',defaultVal=''),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)

110
src/mesh/mesh_mech_FEM.f90
View File

@ -4,7 +4,7 @@
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief FEM PETSc solver !> @brief FEM PETSc solver
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module mesh_mech_FEM module mesh_mechanical_FEM
#include <petsc/finclude/petscdmplex.h> #include <petsc/finclude/petscdmplex.h>
#include <petsc/finclude/petscdm.h> #include <petsc/finclude/petscdm.h>
#include <petsc/finclude/petsc.h> #include <petsc/finclude/petsc.h>
@ -50,7 +50,7 @@ module mesh_mech_FEM
type(tNumerics), private :: num type(tNumerics), private :: num
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
SNES :: mech_snes SNES :: mechanical_snes
Vec :: solution, solution_rate, solution_local Vec :: solution, solution_rate, solution_local
PetscInt :: dimPlex, cellDof, nQuadrature, nBasis PetscInt :: dimPlex, cellDof, nQuadrature, nBasis
PetscReal, allocatable, target :: qPoints(:), qWeights(:) PetscReal, allocatable, target :: qPoints(:), qWeights(:)
@ -65,20 +65,20 @@ module mesh_mech_FEM
real(pReal), parameter :: eps = 1.0e-18_pReal real(pReal), parameter :: eps = 1.0e-18_pReal
public :: & public :: &
FEM_mech_init, & FEM_mechanical_init, &
FEM_mech_solution, & FEM_mechanical_solution, &
FEM_mech_forward FEM_mechanical_forward
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields and fills them with data !> @brief allocates all neccessary fields and fills them with data
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_mech_init(fieldBC) subroutine FEM_mechanical_init(fieldBC)
type(tFieldBC), intent(in) :: fieldBC type(tFieldBC), intent(in) :: fieldBC
DM :: mech_mesh DM :: mechanical_mesh
PetscFE :: mechFE PetscFE :: mechFE
PetscQuadrature :: mechQuad, functional PetscQuadrature :: mechQuad, functional
PetscDS :: mechDS PetscDS :: mechDS
@ -126,8 +126,8 @@ subroutine FEM_mech_init(fieldBC)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Setup FEM mech mesh ! Setup FEM mech mesh
call DMClone(geomMesh,mech_mesh,ierr); CHKERRQ(ierr) call DMClone(geomMesh,mechanical_mesh,ierr); CHKERRQ(ierr)
call DMGetDimension(mech_mesh,dimPlex,ierr); CHKERRQ(ierr) call DMGetDimension(mechanical_mesh,dimPlex,ierr); CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Setup FEM mech discretization ! Setup FEM mech discretization
@ -146,22 +146,22 @@ subroutine FEM_mech_init(fieldBC)
call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr) call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr)
call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr) call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr)
nBasis = nBasis/nc nBasis = nBasis/nc
call DMAddField(mech_mesh,PETSC_NULL_DMLABEL,mechFE,ierr); CHKERRQ(ierr) call DMAddField(mechanical_mesh,PETSC_NULL_DMLABEL,mechFE,ierr); CHKERRQ(ierr)
call DMCreateDS(mech_mesh,ierr); CHKERRQ(ierr) call DMCreateDS(mechanical_mesh,ierr); CHKERRQ(ierr)
call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr) call DMGetDS(mechanical_mesh,mechDS,ierr); CHKERRQ(ierr)
call PetscDSGetTotalDimension(mechDS,cellDof,ierr); CHKERRQ(ierr) call PetscDSGetTotalDimension(mechDS,cellDof,ierr); CHKERRQ(ierr)
call PetscFEDestroy(mechFE,ierr); CHKERRQ(ierr) call PetscFEDestroy(mechFE,ierr); CHKERRQ(ierr)
call PetscQuadratureDestroy(mechQuad,ierr); CHKERRQ(ierr) call PetscQuadratureDestroy(mechQuad,ierr); CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Setup FEM mech boundary conditions ! Setup FEM mech boundary conditions
call DMGetLabel(mech_mesh,'Face Sets',BCLabel,ierr); CHKERRQ(ierr) call DMGetLabel(mechanical_mesh,'Face Sets',BCLabel,ierr); CHKERRQ(ierr)
call DMPlexLabelComplete(mech_mesh,BCLabel,ierr); CHKERRQ(ierr) call DMPlexLabelComplete(mechanical_mesh,BCLabel,ierr); CHKERRQ(ierr)
call DMGetLocalSection(mech_mesh,section,ierr); CHKERRQ(ierr) call DMGetLocalSection(mechanical_mesh,section,ierr); CHKERRQ(ierr)
allocate(pnumComp(1), source=dimPlex) allocate(pnumComp(1), source=dimPlex)
allocate(pnumDof(0:dimPlex), source = 0) allocate(pnumDof(0:dimPlex), source = 0)
do topologDim = 0, dimPlex do topologDim = 0, dimPlex
call DMPlexGetDepthStratum(mech_mesh,topologDim,cellStart,cellEnd,ierr) call DMPlexGetDepthStratum(mechanical_mesh,topologDim,cellStart,cellEnd,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscSectionGetDof(section,cellStart,pnumDof(topologDim),ierr) call PetscSectionGetDof(section,cellStart,pnumDof(topologDim),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
@ -179,10 +179,10 @@ subroutine FEM_mech_init(fieldBC)
numBC = numBC + 1 numBC = numBC + 1
call ISCreateGeneral(PETSC_COMM_WORLD,1,[field-1],PETSC_COPY_VALUES,pbcComps(numBC),ierr) call ISCreateGeneral(PETSC_COMM_WORLD,1,[field-1],PETSC_COPY_VALUES,pbcComps(numBC),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call DMGetStratumSize(mech_mesh,'Face Sets',mesh_boundaries(faceSet),bcSize,ierr) call DMGetStratumSize(mechanical_mesh,'Face Sets',mesh_boundaries(faceSet),bcSize,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
if (bcSize > 0) then if (bcSize > 0) then
call DMGetStratumIS(mech_mesh,'Face Sets',mesh_boundaries(faceSet),bcPoint,ierr) call DMGetStratumIS(mechanical_mesh,'Face Sets',mesh_boundaries(faceSet),bcPoint,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call ISGetIndicesF90(bcPoint,pBcPoint,ierr); CHKERRQ(ierr) call ISGetIndicesF90(bcPoint,pBcPoint,ierr); CHKERRQ(ierr)
call ISCreateGeneral(PETSC_COMM_WORLD,bcSize,pBcPoint,PETSC_COPY_VALUES,pbcPoints(numBC),ierr) call ISCreateGeneral(PETSC_COMM_WORLD,bcSize,pBcPoint,PETSC_COPY_VALUES,pbcPoints(numBC),ierr)
@ -195,32 +195,32 @@ subroutine FEM_mech_init(fieldBC)
endif endif
endif endif
enddo; enddo enddo; enddo
call DMPlexCreateSection(mech_mesh,nolabel,pNumComp,pNumDof, & call DMPlexCreateSection(mechanical_mesh,nolabel,pNumComp,pNumDof, &
numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS,section,ierr) numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS,section,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call DMSetSection(mech_mesh,section,ierr); CHKERRQ(ierr) call DMSetSection(mechanical_mesh,section,ierr); CHKERRQ(ierr)
do faceSet = 1, numBC do faceSet = 1, numBC
call ISDestroy(pbcPoints(faceSet),ierr); CHKERRQ(ierr) call ISDestroy(pbcPoints(faceSet),ierr); CHKERRQ(ierr)
enddo enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,mech_snes,ierr);CHKERRQ(ierr) call SNESCreate(PETSC_COMM_WORLD,mechanical_snes,ierr);CHKERRQ(ierr)
call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr) call SNESSetOptionsPrefix(mechanical_snes,'mechanical_',ierr);CHKERRQ(ierr)
call SNESSetDM(mech_snes,mech_mesh,ierr); CHKERRQ(ierr) !< set the mesh for non-linear solver call SNESSetDM(mechanical_snes,mechanical_mesh,ierr); CHKERRQ(ierr) !< set the mesh for non-linear solver
call DMCreateGlobalVector(mech_mesh,solution ,ierr); CHKERRQ(ierr) !< locally owned displacement Dofs call DMCreateGlobalVector(mechanical_mesh,solution ,ierr); CHKERRQ(ierr) !< locally owned displacement Dofs
call DMCreateGlobalVector(mech_mesh,solution_rate ,ierr); CHKERRQ(ierr) !< locally owned velocity Dofs to guess solution at next load step call DMCreateGlobalVector(mechanical_mesh,solution_rate ,ierr); CHKERRQ(ierr) !< locally owned velocity Dofs to guess solution at next load step
call DMCreateLocalVector (mech_mesh,solution_local ,ierr); CHKERRQ(ierr) !< locally owned velocity Dofs to guess solution at next load step call DMCreateLocalVector (mechanical_mesh,solution_local ,ierr); CHKERRQ(ierr) !< locally owned velocity Dofs to guess solution at next load step
call DMSNESSetFunctionLocal(mech_mesh,FEM_mech_formResidual,PETSC_NULL_VEC,ierr) !< function to evaluate residual forces call DMSNESSetFunctionLocal(mechanical_mesh,FEM_mechanical_formResidual,PETSC_NULL_VEC,ierr) !< function to evaluate residual forces
CHKERRQ(ierr) CHKERRQ(ierr)
call DMSNESSetJacobianLocal(mech_mesh,FEM_mech_formJacobian,PETSC_NULL_VEC,ierr) !< function to evaluate stiffness matrix call DMSNESSetJacobianLocal(mechanical_mesh,FEM_mechanical_formJacobian,PETSC_NULL_VEC,ierr) !< function to evaluate stiffness matrix
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr) !< ignore linear solve failures call SNESSetMaxLinearSolveFailures(mechanical_snes, huge(1), ierr); CHKERRQ(ierr) !< ignore linear solve failures
call SNESSetConvergenceTest(mech_snes,FEM_mech_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,ierr) call SNESSetConvergenceTest(mechanical_snes,FEM_mechanical_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetTolerances(mech_snes,1.0,0.0,0.0,num%itmax,num%itmax,ierr) call SNESSetTolerances(mechanical_snes,1.0,0.0,0.0,num%itmax,num%itmax,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr) call SNESSetFromOptions(mechanical_snes,ierr); CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
@ -236,11 +236,11 @@ subroutine FEM_mech_init(fieldBC)
call PetscDSGetDiscretization(mechDS,0,mechFE,ierr) call PetscDSGetDiscretization(mechDS,0,mechFE,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscFEGetDualSpace(mechFE,mechDualSpace,ierr); CHKERRQ(ierr) call PetscFEGetDualSpace(mechFE,mechDualSpace,ierr); CHKERRQ(ierr)
call DMPlexGetHeightStratum(mech_mesh,0,cellStart,cellEnd,ierr) call DMPlexGetHeightStratum(mechanical_mesh,0,cellStart,cellEnd,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements do cell = cellStart, cellEnd-1 !< loop over all elements
x_scal = 0.0_pReal x_scal = 0.0_pReal
call DMPlexComputeCellGeometryAffineFEM(mech_mesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) call DMPlexComputeCellGeometryAffineFEM(mechanical_mesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex]) cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex])
do basis = 0, nBasis*dimPlex-1, dimPlex do basis = 0, nBasis*dimPlex-1, dimPlex
@ -251,17 +251,17 @@ subroutine FEM_mech_init(fieldBC)
x_scal(basis+1:basis+dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0_pReal) x_scal(basis+1:basis+dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0_pReal)
enddo enddo
px_scal => x_scal px_scal => x_scal
call DMPlexVecSetClosure(mech_mesh,section,solution_local,cell,px_scal,5,ierr) call DMPlexVecSetClosure(mechanical_mesh,section,solution_local,cell,px_scal,5,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
enddo enddo
end subroutine FEM_mech_init end subroutine FEM_mechanical_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the FEM load step !> @brief solution for the FEM load step
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
type(tSolutionState) function FEM_mech_solution( & type(tSolutionState) function FEM_mechanical_solution( &
incInfoIn,timeinc,timeinc_old,fieldBC) incInfoIn,timeinc,timeinc_old,fieldBC)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -278,35 +278,35 @@ type(tSolutionState) function FEM_mech_solution( &
SNESConvergedReason :: reason SNESConvergedReason :: reason
incInfo = incInfoIn incInfo = incInfoIn
FEM_mech_solution%converged =.false. FEM_mechanical_solution%converged =.false.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide availabe data ! set module wide availabe data
params%timeinc = timeinc params%timeinc = timeinc
params%fieldBC = fieldBC params%fieldBC = fieldBC
call SNESSolve(mech_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) ! solve mech_snes based on solution guess (result in solution) call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) ! solve mechanical_snes based on solution guess (result in solution)
call SNESGetConvergedReason(mech_snes,reason,ierr); CHKERRQ(ierr) ! solution converged? call SNESGetConvergedReason(mechanical_snes,reason,ierr); CHKERRQ(ierr) ! solution converged?
terminallyIll = .false. terminallyIll = .false.
if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error
FEM_mech_solution%converged = .false. FEM_mechanical_solution%converged = .false.
FEM_mech_solution%iterationsNeeded = num%itmax FEM_mechanical_solution%iterationsNeeded = num%itmax
else ! >= 1 proper convergence (or terminally ill) else ! >= 1 proper convergence (or terminally ill)
FEM_mech_solution%converged = .true. FEM_mechanical_solution%converged = .true.
call SNESGetIterationNumber(mech_snes,FEM_mech_solution%iterationsNeeded,ierr) call SNESGetIterationNumber(mechanical_snes,FEM_mechanical_solution%iterationsNeeded,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
endif endif
print'(/,a)', ' ===========================================================================' print'(/,a)', ' ==========================================================================='
flush(IO_STDOUT) flush(IO_STDOUT)
end function FEM_mech_solution end function FEM_mechanical_solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forms the FEM residual vector !> @brief forms the FEM residual vector
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,ierr)
DM :: dm_local DM :: dm_local
PetscObject,intent(in) :: dummy PetscObject,intent(in) :: dummy
@ -431,13 +431,13 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
enddo enddo
call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
end subroutine FEM_mech_formResidual end subroutine FEM_mechanical_formResidual
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forms the FEM stiffness matrix !> @brief forms the FEM stiffness matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
DM :: dm_local DM :: dm_local
@ -575,13 +575,13 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr) call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
call MatNullSpaceDestroy(matnull,ierr); CHKERRQ(ierr) call MatNullSpaceDestroy(matnull,ierr); CHKERRQ(ierr)
end subroutine FEM_mech_formJacobian end subroutine FEM_mechanical_formJacobian
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forwarding routine !> @brief forwarding routine
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC) subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC)
type(tFieldBC), intent(in) :: & type(tFieldBC), intent(in) :: &
fieldBC fieldBC
@ -603,7 +603,7 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
if (guess .and. .not. cutBack) then if (guess .and. .not. cutBack) then
ForwardData = .True. ForwardData = .True.
homogenization_F0 = homogenization_F homogenization_F0 = homogenization_F
call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local call SNESGetDM(mechanical_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mechanical_snes into dm_local
call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr) call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call VecSet(x_local,0.0_pReal,ierr); CHKERRQ(ierr) call VecSet(x_local,0.0_pReal,ierr); CHKERRQ(ierr)
@ -634,13 +634,13 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
call VecCopy(solution_rate,solution,ierr); CHKERRQ(ierr) call VecCopy(solution_rate,solution,ierr); CHKERRQ(ierr)
call VecScale(solution,timeinc,ierr); CHKERRQ(ierr) call VecScale(solution,timeinc,ierr); CHKERRQ(ierr)
end subroutine FEM_mech_forward end subroutine FEM_mechanical_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief reporting !> @brief reporting
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) subroutine FEM_mechanical_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
SNES :: snes_local SNES :: snes_local
PetscInt :: PETScIter PetscInt :: PETScIter
@ -662,6 +662,6 @@ subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dumm
' Piola--Kirchhoff stress / MPa =',transpose(P_av)*1.e-6_pReal ' Piola--Kirchhoff stress / MPa =',transpose(P_av)*1.e-6_pReal
flush(IO_STDOUT) flush(IO_STDOUT)
end subroutine FEM_mech_converged end subroutine FEM_mechanical_converged
end module mesh_mech_FEM end module mesh_mechanical_FEM

180
src/phase.f90
View File

@ -62,7 +62,7 @@ module phase
phase_Nsources, & !< number of source mechanisms active in each phase phase_Nsources, & !< number of source mechanisms active in each phase
phase_Nkinematics, & !< number of kinematic mechanisms active in each phase phase_Nkinematics, & !< number of kinematic mechanisms active in each phase
phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase
phase_plasticityInstance, & !< instance of particular plasticity of each phase phase_plasticInstance, & !< instance of particular plasticity of each phase
phase_elasticityInstance !< instance of particular elasticity of each phase phase_elasticityInstance !< instance of particular elasticity of each phase
logical, dimension(:), allocatable, public :: & ! ToDo: should be protected (bug in Intel Compiler) logical, dimension(:), allocatable, public :: & ! ToDo: should be protected (bug in Intel Compiler)
@ -75,15 +75,15 @@ module phase
integer, public, protected :: & integer, public, protected :: &
constitutive_plasticity_maxSizeDotState, & phase_plasticity_maxSizeDotState, &
constitutive_source_maxSizeDotState phase_source_maxSizeDotState
interface interface
! == cleaned:begin ================================================================================= ! == cleaned:begin =================================================================================
module subroutine mech_init(phases) module subroutine mechanical_init(phases)
class(tNode), pointer :: phases class(tNode), pointer :: phases
end subroutine mech_init end subroutine mechanical_init
module subroutine damage_init module subroutine damage_init
end subroutine damage_init end subroutine damage_init
@ -93,83 +93,83 @@ module phase
end subroutine thermal_init end subroutine thermal_init
module subroutine mech_results(group,ph) module subroutine mechanical_results(group,ph)
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
integer, intent(in) :: ph integer, intent(in) :: ph
end subroutine mech_results end subroutine mechanical_results
module subroutine damage_results(group,ph) module subroutine damage_results(group,ph)
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
integer, intent(in) :: ph integer, intent(in) :: ph
end subroutine damage_results end subroutine damage_results
module subroutine mech_windForward(ph,me) module subroutine mechanical_windForward(ph,me)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
end subroutine mech_windForward end subroutine mechanical_windForward
module subroutine mech_forward() module subroutine mechanical_forward()
end subroutine mech_forward end subroutine mechanical_forward
module subroutine thermal_forward() module subroutine thermal_forward()
end subroutine thermal_forward end subroutine thermal_forward
module subroutine mech_restore(ce,includeL) module subroutine mechanical_restore(ce,includeL)
integer, intent(in) :: ce integer, intent(in) :: ce
logical, intent(in) :: includeL logical, intent(in) :: includeL
end subroutine mech_restore end subroutine mechanical_restore
module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF) module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF)
real(pReal), intent(in) :: dt real(pReal), intent(in) :: dt
integer, intent(in) :: & integer, intent(in) :: &
co, & !< counter in constituent loop co, & !< counter in constituent loop
ip, & !< counter in integration point loop ip, & !< counter in integration point loop
el !< counter in element loop el !< counter in element loop
real(pReal), dimension(3,3,3,3) :: dPdF real(pReal), dimension(3,3,3,3) :: dPdF
end function constitutive_mech_dPdF end function phase_mechanical_dPdF
module subroutine mech_restartWrite(groupHandle,ph) module subroutine mechanical_restartWrite(groupHandle,ph)
integer(HID_T), intent(in) :: groupHandle integer(HID_T), intent(in) :: groupHandle
integer, intent(in) :: ph integer, intent(in) :: ph
end subroutine mech_restartWrite end subroutine mechanical_restartWrite
module subroutine mech_restartRead(groupHandle,ph) module subroutine mechanical_restartRead(groupHandle,ph)
integer(HID_T), intent(in) :: groupHandle integer(HID_T), intent(in) :: groupHandle
integer, intent(in) :: ph integer, intent(in) :: ph
end subroutine mech_restartRead end subroutine mechanical_restartRead
module function mech_S(ph,me) result(S) module function mechanical_S(ph,me) result(S)
integer, intent(in) :: ph,me integer, intent(in) :: ph,me
real(pReal), dimension(3,3) :: S real(pReal), dimension(3,3) :: S
end function mech_S end function mechanical_S
module function mech_L_p(ph,me) result(L_p) module function mechanical_L_p(ph,me) result(L_p)
integer, intent(in) :: ph,me integer, intent(in) :: ph,me
real(pReal), dimension(3,3) :: L_p real(pReal), dimension(3,3) :: L_p
end function mech_L_p end function mechanical_L_p
module function constitutive_mech_getF(co,ip,el) result(F) module function phase_mechanical_getF(co,ip,el) result(F)
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ip, el
real(pReal), dimension(3,3) :: F real(pReal), dimension(3,3) :: F
end function constitutive_mech_getF end function phase_mechanical_getF
module function mech_F_e(ph,me) result(F_e) module function mechanical_F_e(ph,me) result(F_e)
integer, intent(in) :: ph,me integer, intent(in) :: ph,me
real(pReal), dimension(3,3) :: F_e real(pReal), dimension(3,3) :: F_e
end function mech_F_e end function mechanical_F_e
module function constitutive_mech_getP(co,ip,el) result(P) module function phase_mechanical_getP(co,ip,el) result(P)
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ip, el
real(pReal), dimension(3,3) :: P real(pReal), dimension(3,3) :: P
end function constitutive_mech_getP end function phase_mechanical_getP
module function constitutive_damage_get_phi(co,ip,el) result(phi) module function phase_damage_get_phi(co,ip,el) result(phi)
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ip, el
real(pReal) :: phi real(pReal) :: phi
end function constitutive_damage_get_phi end function phase_damage_get_phi
module function thermal_T(ph,me) result(T) module function thermal_T(ph,me) result(T)
integer, intent(in) :: ph,me integer, intent(in) :: ph,me
@ -182,20 +182,20 @@ module phase
end function thermal_dot_T end function thermal_dot_T
module subroutine constitutive_mech_setF(F,co,ip,el) module subroutine phase_mechanical_setF(F,co,ip,el)
real(pReal), dimension(3,3), intent(in) :: F real(pReal), dimension(3,3), intent(in) :: F
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ip, el
end subroutine constitutive_mech_setF end subroutine phase_mechanical_setF
module subroutine constitutive_thermal_setField(T,dot_T, co,ce) module subroutine phase_thermal_setField(T,dot_T, co,ce)
real(pReal), intent(in) :: T, dot_T real(pReal), intent(in) :: T, dot_T
integer, intent(in) :: ce, co integer, intent(in) :: ce, co
end subroutine constitutive_thermal_setField end subroutine phase_thermal_setField
module subroutine constitutive_damage_set_phi(phi,co,ce) module subroutine phase_damage_set_phi(phi,co,ce)
real(pReal), intent(in) :: phi real(pReal), intent(in) :: phi
integer, intent(in) :: co, ce integer, intent(in) :: co, ce
end subroutine constitutive_damage_set_phi end subroutine phase_damage_set_phi
! == cleaned:end =================================================================================== ! == cleaned:end ===================================================================================
@ -222,13 +222,13 @@ module phase
logical :: converged_ logical :: converged_
end function crystallite_stress end function crystallite_stress
module function constitutive_homogenizedC(ph,me) result(C) module function phase_homogenizedC(ph,me) result(C)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
real(pReal), dimension(6,6) :: C real(pReal), dimension(6,6) :: C
end function constitutive_homogenizedC end function phase_homogenizedC
module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
@ -237,13 +237,13 @@ module phase
real(pReal), intent(inout) :: & real(pReal), intent(inout) :: &
phiDot, & phiDot, &
dPhiDot_dPhi dPhiDot_dPhi
end subroutine constitutive_damage_getRateAndItsTangents end subroutine phase_damage_getRateAndItsTangents
module subroutine constitutive_thermal_getRate(TDot, ph,me) module subroutine phase_thermal_getRate(TDot, ph,me)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
TDot TDot
end subroutine constitutive_thermal_getRate end subroutine phase_thermal_getRate
module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e)
integer, intent(in) :: & integer, intent(in) :: &
@ -281,39 +281,39 @@ module phase
#endif #endif
public :: & public :: &
constitutive_init, & phase_init, &
constitutive_homogenizedC, & phase_homogenizedC, &
constitutive_damage_getRateAndItsTangents, & phase_damage_getRateAndItsTangents, &
constitutive_thermal_getRate, & phase_thermal_getRate, &
constitutive_results, & phase_results, &
constitutive_allocateState, & phase_allocateState, &
constitutive_forward, & phase_forward, &
constitutive_restore, & phase_restore, &
plastic_nonlocal_updateCompatibility, & plastic_nonlocal_updateCompatibility, &
converged, & converged, &
crystallite_init, & crystallite_init, &
crystallite_stress, & crystallite_stress, &
thermal_stress, & thermal_stress, &
constitutive_mech_dPdF, & phase_mechanical_dPdF, &
crystallite_orientations, & crystallite_orientations, &
crystallite_push33ToRef, & crystallite_push33ToRef, &
constitutive_restartWrite, & phase_restartWrite, &
constitutive_restartRead, & phase_restartRead, &
integrateDamageState, & integrateDamageState, &
constitutive_thermal_setField, & phase_thermal_setField, &
constitutive_damage_set_phi, & phase_damage_set_phi, &
constitutive_damage_get_phi, & phase_damage_get_phi, &
constitutive_mech_getP, & phase_mechanical_getP, &
constitutive_mech_setF, & phase_mechanical_setF, &
constitutive_mech_getF, & phase_mechanical_getF, &
constitutive_windForward phase_windForward
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Initialze constitutive models for individual physics !> @brief Initialze constitutive models for individual physics
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_init subroutine phase_init
integer :: & integer :: &
ph, & !< counter in phase loop ph, & !< counter in phase loop
@ -336,12 +336,12 @@ subroutine constitutive_init
phases => config_material%get('phase') phases => config_material%get('phase')
call mech_init(phases) call mechanical_init(phases)
call damage_init call damage_init
call thermal_init(phases) call thermal_init(phases)
constitutive_source_maxSizeDotState = 0 phase_source_maxSizeDotState = 0
PhaseLoop2:do ph = 1,phases%length PhaseLoop2:do ph = 1,phases%length
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! partition and initialize state ! partition and initialize state
@ -350,18 +350,18 @@ subroutine constitutive_init
damageState(ph)%p(so)%state = damageState(ph)%p(so)%state0 damageState(ph)%p(so)%state = damageState(ph)%p(so)%state0
end forall end forall
constitutive_source_maxSizeDotState = max(constitutive_source_maxSizeDotState, & phase_source_maxSizeDotState = max(phase_source_maxSizeDotState, &
maxval(damageState(ph)%p%sizeDotState)) maxval(damageState(ph)%p%sizeDotState))
enddo PhaseLoop2 enddo PhaseLoop2
constitutive_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState) phase_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState)
end subroutine constitutive_init end subroutine phase_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Allocate the components of the state structure for a given phase !> @brief Allocate the components of the state structure for a given phase
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_allocateState(state, & subroutine phase_allocateState(state, &
Nconstituents,sizeState,sizeDotState,sizeDeltaState) Nconstituents,sizeState,sizeDotState,sizeDeltaState)
class(tState), intent(out) :: & class(tState), intent(out) :: &
@ -387,13 +387,13 @@ subroutine constitutive_allocateState(state, &
allocate(state%deltaState (sizeDeltaState,Nconstituents), source=0.0_pReal) allocate(state%deltaState (sizeDeltaState,Nconstituents), source=0.0_pReal)
end subroutine constitutive_allocateState end subroutine phase_allocateState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Restore data after homog cutback. !> @brief Restore data after homog cutback.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_restore(ce,includeL) subroutine phase_restore(ce,includeL)
logical, intent(in) :: includeL logical, intent(in) :: includeL
integer, intent(in) :: ce integer, intent(in) :: ce
@ -410,21 +410,21 @@ subroutine constitutive_restore(ce,includeL)
enddo enddo
enddo enddo
call mech_restore(ce,includeL) call mechanical_restore(ce,includeL)
end subroutine constitutive_restore end subroutine phase_restore
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Forward data after successful increment. !> @brief Forward data after successful increment.
! ToDo: Any guessing for the current states possible? ! ToDo: Any guessing for the current states possible?
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_forward() subroutine phase_forward()
integer :: ph, so integer :: ph, so
call mech_forward() call mechanical_forward()
call thermal_forward() call thermal_forward()
do ph = 1, size(damageState) do ph = 1, size(damageState)
@ -432,13 +432,13 @@ subroutine constitutive_forward()
damageState(ph)%p(so)%state0 = damageState(ph)%p(so)%state damageState(ph)%p(so)%state0 = damageState(ph)%p(so)%state
enddo; enddo enddo; enddo
end subroutine constitutive_forward end subroutine phase_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes constitutive results to HDF5 output file !> @brief writes constitutive results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_results() subroutine phase_results()
integer :: ph integer :: ph
character(len=:), allocatable :: group character(len=:), allocatable :: group
@ -451,12 +451,12 @@ subroutine constitutive_results()
group = '/current/phase/'//trim(material_name_phase(ph))//'/' group = '/current/phase/'//trim(material_name_phase(ph))//'/'
call results_closeGroup(results_addGroup(group)) call results_closeGroup(results_addGroup(group))
call mech_results(group,ph) call mechanical_results(group,ph)
call damage_results(group,ph) call damage_results(group,ph)
enddo enddo
end subroutine constitutive_results end subroutine phase_results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -557,7 +557,7 @@ end subroutine crystallite_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Wind homog inc forward. !> @brief Wind homog inc forward.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_windForward(ip,el) subroutine phase_windForward(ip,el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
@ -572,7 +572,7 @@ subroutine constitutive_windForward(ip,el)
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el) me = material_phaseMemberAt(co,ip,el)
call mech_windForward(ph,me) call mechanical_windForward(ph,me)
do so = 1, phase_Nsources(material_phaseAt(co,el)) do so = 1, phase_Nsources(material_phaseAt(co,el))
damageState(ph)%p(so)%state0(:,me) = damageState(ph)%p(so)%state(:,me) damageState(ph)%p(so)%state0(:,me) = damageState(ph)%p(so)%state(:,me)
@ -580,7 +580,7 @@ subroutine constitutive_windForward(ip,el)
enddo enddo
end subroutine constitutive_windForward end subroutine phase_windForward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -595,11 +595,11 @@ subroutine crystallite_orientations(co,ip,el)
call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(& call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(&
mech_F_e(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el))))) mechanical_F_e(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)))))
if (plasticState(material_phaseAt(1,el))%nonlocal) & if (plasticState(material_phaseAt(1,el))%nonlocal) &
call plastic_nonlocal_updateCompatibility(crystallite_orientation, & call plastic_nonlocal_updateCompatibility(crystallite_orientation, &
phase_plasticityInstance(material_phaseAt(1,el)),ip,el) phase_plasticInstance(material_phaseAt(1,el)),ip,el)
end subroutine crystallite_orientations end subroutine crystallite_orientations
@ -620,7 +620,7 @@ function crystallite_push33ToRef(co,ip,el, tensor33)
real(pReal), dimension(3,3) :: T real(pReal), dimension(3,3) :: T
T = matmul(material_orientation0(co,ip,el)%asMatrix(),transpose(math_inv33(constitutive_mech_getF(co,ip,el)))) ! ToDo: initial orientation correct? T = matmul(material_orientation0(co,ip,el)%asMatrix(),transpose(math_inv33(phase_mechanical_getF(co,ip,el)))) ! ToDo: initial orientation correct?
crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T))
@ -648,7 +648,7 @@ end function converged
!> @brief Write current restart information (Field and constitutive data) to file. !> @brief Write current restart information (Field and constitutive data) to file.
! ToDo: Merge data into one file for MPI ! ToDo: Merge data into one file for MPI
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_restartWrite(fileHandle) subroutine phase_restartWrite(fileHandle)
integer(HID_T), intent(in) :: fileHandle integer(HID_T), intent(in) :: fileHandle
@ -662,7 +662,7 @@ subroutine constitutive_restartWrite(fileHandle)
groupHandle(2) = HDF5_addGroup(groupHandle(1),material_name_phase(ph)) groupHandle(2) = HDF5_addGroup(groupHandle(1),material_name_phase(ph))
call mech_restartWrite(groupHandle(2),ph) call mechanical_restartWrite(groupHandle(2),ph)
call HDF5_closeGroup(groupHandle(2)) call HDF5_closeGroup(groupHandle(2))
@ -670,14 +670,14 @@ subroutine constitutive_restartWrite(fileHandle)
call HDF5_closeGroup(groupHandle(1)) call HDF5_closeGroup(groupHandle(1))
end subroutine constitutive_restartWrite end subroutine phase_restartWrite
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Read data for restart !> @brief Read data for restart
! ToDo: Merge data into one file for MPI ! ToDo: Merge data into one file for MPI
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_restartRead(fileHandle) subroutine phase_restartRead(fileHandle)
integer(HID_T), intent(in) :: fileHandle integer(HID_T), intent(in) :: fileHandle
@ -691,7 +691,7 @@ subroutine constitutive_restartRead(fileHandle)
groupHandle(2) = HDF5_openGroup(groupHandle(1),material_name_phase(ph)) groupHandle(2) = HDF5_openGroup(groupHandle(1),material_name_phase(ph))
call mech_restartRead(groupHandle(2),ph) call mechanical_restartRead(groupHandle(2),ph)
call HDF5_closeGroup(groupHandle(2)) call HDF5_closeGroup(groupHandle(2))
@ -699,7 +699,7 @@ subroutine constitutive_restartRead(fileHandle)
call HDF5_closeGroup(groupHandle(1)) call HDF5_closeGroup(groupHandle(1))
end subroutine constitutive_restartRead end subroutine phase_restartRead
end module phase end module phase

34
src/phase_damage.f90
View File

@ -196,7 +196,7 @@ end subroutine damage_init
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief returns local part of nonlocal damage driving force !< @brief returns local part of nonlocal damage driving force
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
@ -246,7 +246,7 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi
enddo enddo
enddo enddo
end subroutine constitutive_damage_getRateAndItsTangents end subroutine phase_damage_getRateAndItsTangents
@ -272,9 +272,9 @@ module function integrateDamageState(dt,co,ip,el) result(broken)
size_so size_so
real(pReal) :: & real(pReal) :: &
zeta zeta
real(pReal), dimension(constitutive_source_maxSizeDotState) :: & real(pReal), dimension(phase_source_maxSizeDotState) :: &
r ! state residuum r ! state residuum
real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState real(pReal), dimension(phase_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState
logical :: & logical :: &
converged_ converged_
@ -283,7 +283,7 @@ module function integrateDamageState(dt,co,ip,el) result(broken)
me = material_phaseMemberAt(co,ip,el) me = material_phaseMemberAt(co,ip,el)
converged_ = .true. converged_ = .true.
broken = constitutive_damage_collectDotState(co,ip,el,ph,me) broken = phase_damage_collectDotState(co,ip,el,ph,me)
if(broken) return if(broken) return
do so = 1, phase_Nsources(ph) do so = 1, phase_Nsources(ph)
@ -300,7 +300,7 @@ module function integrateDamageState(dt,co,ip,el) result(broken)
source_dotState(1:size_so(so),1,so) = damageState(ph)%p(so)%dotState(:,me) source_dotState(1:size_so(so),1,so) = damageState(ph)%p(so)%dotState(:,me)
enddo enddo
broken = constitutive_damage_collectDotState(co,ip,el,ph,me) broken = phase_damage_collectDotState(co,ip,el,ph,me)
if(broken) exit iteration if(broken) exit iteration
do so = 1, phase_Nsources(ph) do so = 1, phase_Nsources(ph)
@ -320,7 +320,7 @@ module function integrateDamageState(dt,co,ip,el) result(broken)
enddo enddo
if(converged_) then if(converged_) then
broken = constitutive_damage_deltaState(mech_F_e(ph,me),ph,me) broken = phase_damage_deltaState(mechanical_F_e(ph,me),ph,me)
exit iteration exit iteration
endif endif
@ -393,7 +393,7 @@ end subroutine damage_results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure !> @brief contains the constitutive equation for calculating the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_damage_collectDotState(co,ip,el,ph,me) result(broken) function phase_damage_collectDotState(co,ip,el,ph,me) result(broken)
integer, intent(in) :: & integer, intent(in) :: &
co, & !< component-ID me integration point co, & !< component-ID me integration point
@ -419,7 +419,7 @@ function constitutive_damage_collectDotState(co,ip,el,ph,me) result(broken)
call anisoductile_dotState(co, ip, el) call anisoductile_dotState(co, ip, el)
case (DAMAGE_ANISOBRITTLE_ID) sourceType case (DAMAGE_ANISOBRITTLE_ID) sourceType
call anisobrittle_dotState(mech_S(ph,me),co, ip, el) ! correct stress? call anisobrittle_dotState(mechanical_S(ph,me),co, ip, el) ! correct stress?
end select sourceType end select sourceType
@ -427,7 +427,7 @@ function constitutive_damage_collectDotState(co,ip,el,ph,me) result(broken)
enddo SourceLoop enddo SourceLoop
end function constitutive_damage_collectDotState end function phase_damage_collectDotState
@ -435,7 +435,7 @@ end function constitutive_damage_collectDotState
!> @brief for constitutive models having an instantaneous change of state !> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model !> will return false if delta state is not needed/supported by the constitutive model
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_damage_deltaState(Fe, ph, me) result(broken) function phase_damage_deltaState(Fe, ph, me) result(broken)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
@ -457,7 +457,7 @@ function constitutive_damage_deltaState(Fe, ph, me) result(broken)
sourceType: select case (phase_source(so,ph)) sourceType: select case (phase_source(so,ph))
case (DAMAGE_ISOBRITTLE_ID) sourceType case (DAMAGE_ISOBRITTLE_ID) sourceType
call source_damage_isoBrittle_deltaState(constitutive_homogenizedC(ph,me), Fe, ph,me) call source_damage_isoBrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me)
broken = any(IEEE_is_NaN(damageState(ph)%p(so)%deltaState(:,me))) broken = any(IEEE_is_NaN(damageState(ph)%p(so)%deltaState(:,me)))
if(.not. broken) then if(.not. broken) then
myOffset = damageState(ph)%p(so)%offsetDeltaState myOffset = damageState(ph)%p(so)%offsetDeltaState
@ -470,7 +470,7 @@ function constitutive_damage_deltaState(Fe, ph, me) result(broken)
enddo SourceLoop enddo SourceLoop
end function constitutive_damage_deltaState end function phase_damage_deltaState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -507,7 +507,7 @@ end function source_active
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Set damage parameter !< @brief Set damage parameter
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module subroutine constitutive_damage_set_phi(phi,co,ce) module subroutine phase_damage_set_phi(phi,co,ce)
real(pReal), intent(in) :: phi real(pReal), intent(in) :: phi
integer, intent(in) :: ce, co integer, intent(in) :: ce, co
@ -515,17 +515,17 @@ module subroutine constitutive_damage_set_phi(phi,co,ce)
current(material_phaseAt2(co,ce))%phi(material_phaseMemberAt2(co,ce)) = phi current(material_phaseAt2(co,ce))%phi(material_phaseMemberAt2(co,ce)) = phi
end subroutine constitutive_damage_set_phi end subroutine phase_damage_set_phi
module function constitutive_damage_get_phi(co,ip,el) result(phi) module function phase_damage_get_phi(co,ip,el) result(phi)
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ip, el
real(pReal) :: phi real(pReal) :: phi
phi = current(material_phaseAt(co,el))%phi(material_phaseMemberAt(co,ip,el)) phi = current(material_phaseAt(co,el))%phi(material_phaseMemberAt(co,ip,el))
end function constitutive_damage_get_phi end function phase_damage_get_phi
end submodule damagee end submodule damagee

2
src/phase_damage_anisobrittle.f90
View File

@ -101,7 +101,7 @@ module function anisobrittle_init(source_length) result(mySources)
if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit' if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit'
Nconstituents = count(material_phaseAt==p) * discretization_nIPs Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0)
damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol'

2
src/phase_damage_anisoductile.f90
View File

@ -87,7 +87,7 @@ module function anisoductile_init(source_length) result(mySources)
if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit' if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit'
Nconstituents=count(material_phaseAt==p) * discretization_nIPs Nconstituents=count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0)
damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol'

2
src/phase_damage_isobrittle.f90
View File

@ -74,7 +74,7 @@ module function isobrittle_init(source_length) result(mySources)
if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit' if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit'
Nconstituents = count(material_phaseAt==p) * discretization_nIPs Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,1) call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,1)
damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'

2
src/phase_damage_isoductile.f90
View File

@ -78,7 +78,7 @@ module function isoductile_init(source_length) result(mySources)
if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit'
Nconstituents=count(material_phaseAt==p) * discretization_nIPs Nconstituents=count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0)
damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'

384
src/phase_mechanics.f90 → src/phase_mechanical.f90
View File

@ -31,21 +31,21 @@ submodule(phase) mechanics
type(tTensorContainer), dimension(:), allocatable :: & type(tTensorContainer), dimension(:), allocatable :: &
! current value ! current value
constitutive_mech_Fe, & phase_mechanical_Fe, &
constitutive_mech_Fi, & phase_mechanical_Fi, &
constitutive_mech_Fp, & phase_mechanical_Fp, &
constitutive_mech_F, & phase_mechanical_F, &
constitutive_mech_Li, & phase_mechanical_Li, &
constitutive_mech_Lp, & phase_mechanical_Lp, &
constitutive_mech_S, & phase_mechanical_S, &
constitutive_mech_P, & phase_mechanical_P, &
! converged value at end of last solver increment ! converged value at end of last solver increment
constitutive_mech_Fi0, & phase_mechanical_Fi0, &
constitutive_mech_Fp0, & phase_mechanical_Fp0, &
constitutive_mech_F0, & phase_mechanical_F0, &
constitutive_mech_Li0, & phase_mechanical_Li0, &
constitutive_mech_Lp0, & phase_mechanical_Lp0, &
constitutive_mech_S0 phase_mechanical_S0
integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: & integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: &
@ -97,7 +97,7 @@ submodule(phase) mechanics
broken broken
end function plastic_deltaState end function plastic_deltaState
module subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
S, Fi, co, ip, el) S, Fi, co, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
@ -114,7 +114,7 @@ submodule(phase) mechanics
dLi_dS, & !< derivative of Li with respect to S dLi_dS, & !< derivative of Li with respect to S
dLi_dFi dLi_dFi
end subroutine constitutive_LiAndItsTangents end subroutine phase_LiAndItsTangents
module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
@ -186,7 +186,7 @@ contains
!> @brief Initialize mechanical field related constitutive models !> @brief Initialize mechanical field related constitutive models
!> @details Initialize elasticity, plasticity and stiffness degradation models. !> @details Initialize elasticity, plasticity and stiffness degradation models.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_init(phases) module subroutine mechanical_init(phases)
class(tNode), pointer :: & class(tNode), pointer :: &
phases phases
@ -215,38 +215,38 @@ module subroutine mech_init(phases)
allocate(phase_NstiffnessDegradations(phases%length),source=0) allocate(phase_NstiffnessDegradations(phases%length),source=0)
allocate(output_constituent(phases%length)) allocate(output_constituent(phases%length))
allocate(constitutive_mech_Fe(phases%length)) allocate(phase_mechanical_Fe(phases%length))
allocate(constitutive_mech_Fi(phases%length)) allocate(phase_mechanical_Fi(phases%length))
allocate(constitutive_mech_Fi0(phases%length)) allocate(phase_mechanical_Fi0(phases%length))
allocate(constitutive_mech_Fp(phases%length)) allocate(phase_mechanical_Fp(phases%length))
allocate(constitutive_mech_Fp0(phases%length)) allocate(phase_mechanical_Fp0(phases%length))
allocate(constitutive_mech_F(phases%length)) allocate(phase_mechanical_F(phases%length))
allocate(constitutive_mech_F0(phases%length)) allocate(phase_mechanical_F0(phases%length))
allocate(constitutive_mech_Li(phases%length)) allocate(phase_mechanical_Li(phases%length))
allocate(constitutive_mech_Li0(phases%length)) allocate(phase_mechanical_Li0(phases%length))
allocate(constitutive_mech_Lp0(phases%length)) allocate(phase_mechanical_Lp0(phases%length))
allocate(constitutive_mech_Lp(phases%length)) allocate(phase_mechanical_Lp(phases%length))
allocate(constitutive_mech_S(phases%length)) allocate(phase_mechanical_S(phases%length))
allocate(constitutive_mech_P(phases%length)) allocate(phase_mechanical_P(phases%length))
allocate(constitutive_mech_S0(phases%length)) allocate(phase_mechanical_S0(phases%length))
do ph = 1, phases%length do ph = 1, phases%length
Nconstituents = count(material_phaseAt == ph) * discretization_nIPs Nconstituents = count(material_phaseAt == ph) * discretization_nIPs
allocate(constitutive_mech_Fi(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Fi(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Fe(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Fe(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Fi0(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Fi0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Fp(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Fp(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Fp0(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Fp0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Li(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Li(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Li0(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Li0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Lp0(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Lp0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Lp(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Lp(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_S(ph)%data(3,3,Nconstituents),source=0.0_pReal) allocate(phase_mechanical_S(ph)%data(3,3,Nconstituents),source=0.0_pReal)
allocate(constitutive_mech_P(ph)%data(3,3,Nconstituents),source=0.0_pReal) allocate(phase_mechanical_P(ph)%data(3,3,Nconstituents),source=0.0_pReal)
allocate(constitutive_mech_S0(ph)%data(3,3,Nconstituents),source=0.0_pReal) allocate(phase_mechanical_S0(ph)%data(3,3,Nconstituents),source=0.0_pReal)
allocate(constitutive_mech_F(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_F(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_F0(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_F0(ph)%data(3,3,Nconstituents))
phase => phases%get(ph) phase => phases%get(ph)
mech => phase%get('mechanics') mech => phase%get('mechanics')
@ -287,17 +287,17 @@ module subroutine mech_init(phases)
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el) me = material_phaseMemberAt(co,ip,el)
constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = material_orientation0(co,ip,el)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = material_orientation0(co,ip,el)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) & phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = phase_mechanical_Fp0(ph)%data(1:3,1:3,me) &
/ math_det33(constitutive_mech_Fp0(ph)%data(1:3,1:3,me))**(1.0_pReal/3.0_pReal) / math_det33(phase_mechanical_Fp0(ph)%data(1:3,1:3,me))**(1.0_pReal/3.0_pReal)
constitutive_mech_Fi0(ph)%data(1:3,1:3,me) = math_I3 phase_mechanical_Fi0(ph)%data(1:3,1:3,me) = math_I3
constitutive_mech_F0(ph)%data(1:3,1:3,me) = math_I3 phase_mechanical_F0(ph)%data(1:3,1:3,me) = math_I3
constitutive_mech_Fe(ph)%data(1:3,1:3,me) = math_inv33(matmul(constitutive_mech_Fi0(ph)%data(1:3,1:3,me), & phase_mechanical_Fe(ph)%data(1:3,1:3,me) = math_inv33(matmul(phase_mechanical_Fi0(ph)%data(1:3,1:3,me), &
constitutive_mech_Fp0(ph)%data(1:3,1:3,me))) ! assuming that euler angles are given in internal strain free configuration phase_mechanical_Fp0(ph)%data(1:3,1:3,me))) ! assuming that euler angles are given in internal strain free configuration
constitutive_mech_Fp(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) phase_mechanical_Fp(ph)%data(1:3,1:3,me) = phase_mechanical_Fp0(ph)%data(1:3,1:3,me)
constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) phase_mechanical_Fi(ph)%data(1:3,1:3,me) = phase_mechanical_Fi0(ph)%data(1:3,1:3,me)
constitutive_mech_F(ph)%data(1:3,1:3,me) = constitutive_mech_F0(ph)%data(1:3,1:3,me) phase_mechanical_F(ph)%data(1:3,1:3,me) = phase_mechanical_F0(ph)%data(1:3,1:3,me)
enddo enddo
enddo; enddo enddo; enddo
@ -307,14 +307,14 @@ module subroutine mech_init(phases)
! initialize plasticity ! initialize plasticity
allocate(plasticState(phases%length)) allocate(plasticState(phases%length))
allocate(phase_plasticity(phases%length),source = PLASTICITY_undefined_ID) allocate(phase_plasticity(phases%length),source = PLASTICITY_undefined_ID)
allocate(phase_plasticityInstance(phases%length),source = 0) allocate(phase_plasticInstance(phases%length),source = 0)
allocate(phase_localPlasticity(phases%length), source=.true.) allocate(phase_localPlasticity(phases%length), source=.true.)
call plastic_init() call plastic_init()
do ph = 1, phases%length do ph = 1, phases%length
phase_elasticityInstance(ph) = count(phase_elasticity(1:ph) == phase_elasticity(ph)) phase_elasticityInstance(ph) = count(phase_elasticity(1:ph) == phase_elasticity(ph))
phase_plasticityInstance(ph) = count(phase_plasticity(1:ph) == phase_plasticity(ph)) phase_plasticInstance(ph) = count(phase_plasticity(1:ph) == phase_plasticity(ph))
enddo enddo
num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
@ -345,14 +345,14 @@ module subroutine mech_init(phases)
call eigendeformation_init(phases) call eigendeformation_init(phases)
end subroutine mech_init end subroutine mechanical_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to !> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
!> the elastic and intermediate deformation gradients using Hooke's law !> the elastic and intermediate deformation gradients using Hooke's law
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
Fe, Fi, co, ip, el) Fe, Fi, co, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
@ -376,7 +376,7 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
i, j, ph, me i, j, ph, me
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
C = math_66toSym3333(constitutive_homogenizedC(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el))) C = math_66toSym3333(phase_homogenizedC(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)))
DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(co,el)) DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(co,el))
degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(co,el))) degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(co,el)))
@ -393,10 +393,10 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
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 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 enddo; enddo
end subroutine constitutive_hooke_SandItsTangents end subroutine phase_hooke_SandItsTangents
module subroutine mech_results(group,ph) module subroutine mechanical_results(group,ph)
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
integer, intent(in) :: ph integer, intent(in) :: ph
@ -407,28 +407,28 @@ module subroutine mech_results(group,ph)
select case(phase_plasticity(ph)) select case(phase_plasticity(ph))
case(PLASTICITY_ISOTROPIC_ID) case(PLASTICITY_ISOTROPIC_ID)
call plastic_isotropic_results(phase_plasticityInstance(ph),group//'plastic/') call plastic_isotropic_results(phase_plasticInstance(ph),group//'plastic/')
case(PLASTICITY_PHENOPOWERLAW_ID) case(PLASTICITY_PHENOPOWERLAW_ID)
call plastic_phenopowerlaw_results(phase_plasticityInstance(ph),group//'plastic/') call plastic_phenopowerlaw_results(phase_plasticInstance(ph),group//'plastic/')
case(PLASTICITY_KINEHARDENING_ID) case(PLASTICITY_KINEHARDENING_ID)
call plastic_kinehardening_results(phase_plasticityInstance(ph),group//'plastic/') call plastic_kinehardening_results(phase_plasticInstance(ph),group//'plastic/')
case(PLASTICITY_DISLOTWIN_ID) case(PLASTICITY_DISLOTWIN_ID)
call plastic_dislotwin_results(phase_plasticityInstance(ph),group//'plastic/') call plastic_dislotwin_results(phase_plasticInstance(ph),group//'plastic/')
case(PLASTICITY_DISLOTUNGSTEN_ID) case(PLASTICITY_DISLOTUNGSTEN_ID)
call plastic_dislotungsten_results(phase_plasticityInstance(ph),group//'plastic/') call plastic_dislotungsten_results(phase_plasticInstance(ph),group//'plastic/')
case(PLASTICITY_NONLOCAL_ID) case(PLASTICITY_NONLOCAL_ID)
call plastic_nonlocal_results(phase_plasticityInstance(ph),group//'plastic/') call plastic_nonlocal_results(phase_plasticInstance(ph),group//'plastic/')
end select end select
call crystallite_results(group,ph) call crystallite_results(group,ph)
end subroutine mech_results end subroutine mechanical_results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -503,8 +503,8 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
call plastic_dependentState(co,ip,el) call plastic_dependentState(co,ip,el)
Lpguess = constitutive_mech_Lp(ph)%data(1:3,1:3,me) ! take as first guess Lpguess = phase_mechanical_Lp(ph)%data(1:3,1:3,me) ! take as first guess
Liguess = constitutive_mech_Li(ph)%data(1:3,1:3,me) ! take as first guess Liguess = phase_mechanical_Li(ph)%data(1:3,1:3,me) ! take as first guess
call math_invert33(invFp_current,devNull,error,subFp0) call math_invert33(invFp_current,devNull,error,subFp0)
if (error) return ! error if (error) return ! error
@ -538,7 +538,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
B = math_I3 - Delta_t*Lpguess B = math_I3 - Delta_t*Lpguess
Fe = matmul(matmul(A,B), invFi_new) Fe = matmul(matmul(A,B), invFi_new)
call constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & call phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
Fe, Fi_new, co, ip, el) Fe, Fi_new, co, ip, el)
call plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & call plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, &
@ -582,7 +582,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
+ deltaLp * steplengthLp + deltaLp * steplengthLp
enddo LpLoop enddo LpLoop
call constitutive_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, & call phase_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, &
S, Fi_new, co, ip, el) S, Fi_new, co, ip, el)
!* update current residuum and check for convergence of loop !* update current residuum and check for convergence of loop
@ -633,13 +633,13 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
call math_invert33(Fp_new,devNull,error,invFp_new) call math_invert33(Fp_new,devNull,error,invFp_new)
if (error) return ! error if (error) return ! error
constitutive_mech_P(ph)%data(1:3,1:3,me) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) phase_mechanical_P(ph)%data(1:3,1:3,me) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new)))
constitutive_mech_S(ph)%data(1:3,1:3,me) = S phase_mechanical_S(ph)%data(1:3,1:3,me) = S
constitutive_mech_Lp(ph)%data(1:3,1:3,me) = Lpguess phase_mechanical_Lp(ph)%data(1:3,1:3,me) = Lpguess
constitutive_mech_Li(ph)%data(1:3,1:3,me) = Liguess phase_mechanical_Li(ph)%data(1:3,1:3,me) = Liguess
constitutive_mech_Fp(ph)%data(1:3,1:3,me) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize phase_mechanical_Fp(ph)%data(1:3,1:3,me) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize
constitutive_mech_Fi(ph)%data(1:3,1:3,me) = Fi_new phase_mechanical_Fi(ph)%data(1:3,1:3,me) = Fi_new
constitutive_mech_Fe(ph)%data(1:3,1:3,me) = matmul(matmul(F,invFp_new),invFi_new) phase_mechanical_Fe(ph)%data(1:3,1:3,me) = matmul(matmul(F,invFp_new),invFi_new)
broken = .false. broken = .false.
end function integrateStress end function integrateStress
@ -668,9 +668,9 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
sizeDotState sizeDotState
real(pReal) :: & real(pReal) :: &
zeta zeta
real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: & real(pReal), dimension(phase_plasticity_maxSizeDotState) :: &
r ! state residuum r ! state residuum
real(pReal), dimension(constitutive_plasticity_maxSizeDotState,2) :: & real(pReal), dimension(phase_plasticity_maxSizeDotState,2) :: &
dotState dotState
@ -796,7 +796,7 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip
ph, & ph, &
me, & me, &
sizeDotState sizeDotState
real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic real(pReal), dimension(phase_plasticity_maxSizeDotState) :: residuum_plastic
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
@ -914,7 +914,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D
ph, & ph, &
me, & me, &
sizeDotState sizeDotState
real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState real(pReal), dimension(phase_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
@ -987,28 +987,28 @@ subroutine crystallite_results(group,ph)
select case (output_constituent(ph)%label(ou)) select case (output_constituent(ph)%label(ou))
case('F') case('F')
call results_writeDataset(group//'/mechanics/',constitutive_mech_F(ph)%data,'F',& call results_writeDataset(group//'/mechanics/',phase_mechanical_F(ph)%data,'F',&
'deformation gradient','1') 'deformation gradient','1')
case('F_e') case('F_e')
call results_writeDataset(group//'/mechanics/',constitutive_mech_Fe(ph)%data,'F_e',& call results_writeDataset(group//'/mechanics/',phase_mechanical_Fe(ph)%data,'F_e',&
'elastic deformation gradient','1') 'elastic deformation gradient','1')
case('F_p') case('F_p')
call results_writeDataset(group//'/mechanics/',constitutive_mech_Fp(ph)%data,'F_p', & call results_writeDataset(group//'/mechanics/',phase_mechanical_Fp(ph)%data,'F_p', &
'plastic deformation gradient','1') 'plastic deformation gradient','1')
case('F_i') case('F_i')
call results_writeDataset(group//'/mechanics/',constitutive_mech_Fi(ph)%data,'F_i', & call results_writeDataset(group//'/mechanics/',phase_mechanical_Fi(ph)%data,'F_i', &
'inelastic deformation gradient','1') 'inelastic deformation gradient','1')
case('L_p') case('L_p')
call results_writeDataset(group//'/mechanics/',constitutive_mech_Lp(ph)%data,'L_p', & call results_writeDataset(group//'/mechanics/',phase_mechanical_Lp(ph)%data,'L_p', &
'plastic velocity gradient','1/s') 'plastic velocity gradient','1/s')
case('L_i') case('L_i')
call results_writeDataset(group//'/mechanics/',constitutive_mech_Li(ph)%data,'L_i', & call results_writeDataset(group//'/mechanics/',phase_mechanical_Li(ph)%data,'L_i', &
'inelastic velocity gradient','1/s') 'inelastic velocity gradient','1/s')
case('P') case('P')
call results_writeDataset(group//'/mechanics/',constitutive_mech_P(ph)%data,'P', & call results_writeDataset(group//'/mechanics/',phase_mechanical_P(ph)%data,'P', &
'First Piola-Kirchhoff stress','Pa') 'First Piola-Kirchhoff stress','Pa')
case('S') case('S')
call results_writeDataset(group//'/mechanics/',constitutive_mech_S(ph)%data,'S', & call results_writeDataset(group//'/mechanics/',phase_mechanical_S(ph)%data,'S', &
'Second Piola-Kirchhoff stress','Pa') 'Second Piola-Kirchhoff stress','Pa')
case('O') case('O')
select case(lattice_structure(ph)) select case(lattice_structure(ph))
@ -1067,43 +1067,43 @@ end subroutine crystallite_results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Wind homog inc forward. !> @brief Wind homog inc forward.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_windForward(ph,me) module subroutine mechanical_windForward(ph,me)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp(ph)%data(1:3,1:3,me) phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = phase_mechanical_Fp(ph)%data(1:3,1:3,me)
constitutive_mech_Fi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi(ph)%data(1:3,1:3,me) phase_mechanical_Fi0(ph)%data(1:3,1:3,me) = phase_mechanical_Fi(ph)%data(1:3,1:3,me)
constitutive_mech_F0(ph)%data(1:3,1:3,me) = constitutive_mech_F(ph)%data(1:3,1:3,me) phase_mechanical_F0(ph)%data(1:3,1:3,me) = phase_mechanical_F(ph)%data(1:3,1:3,me)
constitutive_mech_Li0(ph)%data(1:3,1:3,me) = constitutive_mech_Li(ph)%data(1:3,1:3,me) phase_mechanical_Li0(ph)%data(1:3,1:3,me) = phase_mechanical_Li(ph)%data(1:3,1:3,me)
constitutive_mech_Lp0(ph)%data(1:3,1:3,me) = constitutive_mech_Lp(ph)%data(1:3,1:3,me) phase_mechanical_Lp0(ph)%data(1:3,1:3,me) = phase_mechanical_Lp(ph)%data(1:3,1:3,me)
constitutive_mech_S0(ph)%data(1:3,1:3,me) = constitutive_mech_S(ph)%data(1:3,1:3,me) phase_mechanical_S0(ph)%data(1:3,1:3,me) = phase_mechanical_S(ph)%data(1:3,1:3,me)
plasticState(ph)%State0(:,me) = plasticState(ph)%state(:,me) plasticState(ph)%State0(:,me) = plasticState(ph)%state(:,me)
end subroutine mech_windForward end subroutine mechanical_windForward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Forward data after successful increment. !> @brief Forward data after successful increment.
! ToDo: Any guessing for the current states possible? ! ToDo: Any guessing for the current states possible?
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_forward() module subroutine mechanical_forward()
integer :: ph integer :: ph
do ph = 1, size(plasticState) do ph = 1, size(plasticState)
constitutive_mech_Fi0(ph) = constitutive_mech_Fi(ph) phase_mechanical_Fi0(ph) = phase_mechanical_Fi(ph)
constitutive_mech_Fp0(ph) = constitutive_mech_Fp(ph) phase_mechanical_Fp0(ph) = phase_mechanical_Fp(ph)
constitutive_mech_F0(ph) = constitutive_mech_F(ph) phase_mechanical_F0(ph) = phase_mechanical_F(ph)
constitutive_mech_Li0(ph) = constitutive_mech_Li(ph) phase_mechanical_Li0(ph) = phase_mechanical_Li(ph)
constitutive_mech_Lp0(ph) = constitutive_mech_Lp(ph) phase_mechanical_Lp0(ph) = phase_mechanical_Lp(ph)
constitutive_mech_S0(ph) = constitutive_mech_S(ph) phase_mechanical_S0(ph) = phase_mechanical_S(ph)
plasticState(ph)%state0 = plasticState(ph)%state plasticState(ph)%state0 = plasticState(ph)%state
enddo enddo
end subroutine mech_forward end subroutine mechanical_forward
@ -1111,19 +1111,19 @@ end subroutine mech_forward
!> @brief returns the homogenize elasticity matrix !> @brief returns the homogenize elasticity matrix
!> ToDo: homogenizedC66 would be more consistent !> ToDo: homogenizedC66 would be more consistent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function constitutive_homogenizedC(ph,me) result(C) module function phase_homogenizedC(ph,me) result(C)
real(pReal), dimension(6,6) :: C real(pReal), dimension(6,6) :: C
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
plasticityType: select case (phase_plasticity(ph)) plasticType: select case (phase_plasticity(ph))
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticType
C = plastic_dislotwin_homogenizedC(ph,me) C = plastic_dislotwin_homogenizedC(ph,me)
case default plasticityType case default plasticType
C = lattice_C66(1:6,1:6,ph) C = lattice_C66(1:6,1:6,ph)
end select plasticityType end select plasticType
end function constitutive_homogenizedC end function phase_homogenizedC
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -1158,17 +1158,17 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
me = material_phaseMemberAt(co,ip,el) me = material_phaseMemberAt(co,ip,el)
sizeDotState = plasticState(ph)%sizeDotState sizeDotState = plasticState(ph)%sizeDotState
subLi0 = constitutive_mech_Li0(ph)%data(1:3,1:3,me) subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,me)
subLp0 = constitutive_mech_Lp0(ph)%data(1:3,1:3,me) subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,me)
subState0 = plasticState(ph)%State0(:,me) subState0 = plasticState(ph)%State0(:,me)
do so = 1, phase_Nsources(ph) do so = 1, phase_Nsources(ph)
damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state0(:,me) damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state0(:,me)
enddo enddo
subFp0 = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,me)
subFi0 = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,me)
subF0 = constitutive_mech_F0(ph)%data(1:3,1:3,me) subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,me)
subFrac = 0.0_pReal subFrac = 0.0_pReal
subStep = 1.0_pReal/num%subStepSizeCryst subStep = 1.0_pReal/num%subStepSizeCryst
todo = .true. todo = .true.
@ -1186,10 +1186,10 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
if (todo) then if (todo) then
subF0 = subF subF0 = subF
subLp0 = constitutive_mech_Lp(ph)%data(1:3,1:3,me) subLp0 = phase_mechanical_Lp(ph)%data(1:3,1:3,me)
subLi0 = constitutive_mech_Li(ph)%data(1:3,1:3,me) subLi0 = phase_mechanical_Li(ph)%data(1:3,1:3,me)
subFp0 = constitutive_mech_Fp(ph)%data(1:3,1:3,me) subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,me)
subFi0 = constitutive_mech_Fi(ph)%data(1:3,1:3,me) subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,me)
subState0 = plasticState(ph)%state(:,me) subState0 = plasticState(ph)%state(:,me)
do so = 1, phase_Nsources(ph) do so = 1, phase_Nsources(ph)
damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state(:,me) damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state(:,me)
@ -1199,12 +1199,12 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
! cut back (reduced time and restore) ! cut back (reduced time and restore)
else else
subStep = num%subStepSizeCryst * subStep subStep = num%subStepSizeCryst * subStep
constitutive_mech_Fp(ph)%data(1:3,1:3,me) = subFp0 phase_mechanical_Fp(ph)%data(1:3,1:3,me) = subFp0
constitutive_mech_Fi(ph)%data(1:3,1:3,me) = subFi0 phase_mechanical_Fi(ph)%data(1:3,1:3,me) = subFi0
constitutive_mech_S(ph)%data(1:3,1:3,me) = constitutive_mech_S0(ph)%data(1:3,1:3,me) ! why no subS0 ? is S0 of any use? phase_mechanical_S(ph)%data(1:3,1:3,me) = phase_mechanical_S0(ph)%data(1:3,1:3,me) ! why no subS0 ? is S0 of any use?
if (subStep < 1.0_pReal) then ! actual (not initial) cutback if (subStep < 1.0_pReal) then ! actual (not initial) cutback
constitutive_mech_Lp(ph)%data(1:3,1:3,me) = subLp0 phase_mechanical_Lp(ph)%data(1:3,1:3,me) = subLp0
constitutive_mech_Li(ph)%data(1:3,1:3,me) = subLi0 phase_mechanical_Li(ph)%data(1:3,1:3,me) = subLi0
endif endif
plasticState(ph)%state(:,me) = subState0 plasticState(ph)%state(:,me) = subState0
do so = 1, phase_Nsources(ph) do so = 1, phase_Nsources(ph)
@ -1218,9 +1218,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
! prepare for integration ! prepare for integration
if (todo) then if (todo) then
subF = subF0 & subF = subF0 &
+ subStep * (constitutive_mech_F(ph)%data(1:3,1:3,me) - constitutive_mech_F0(ph)%data(1:3,1:3,me)) + subStep * (phase_mechanical_F(ph)%data(1:3,1:3,me) - phase_mechanical_F0(ph)%data(1:3,1:3,me))
constitutive_mech_Fe(ph)%data(1:3,1:3,me) = matmul(subF,math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), & phase_mechanical_Fe(ph)%data(1:3,1:3,me) = matmul(subF,math_inv33(matmul(phase_mechanical_Fi(ph)%data(1:3,1:3,me), &
constitutive_mech_Fp(ph)%data(1:3,1:3,me)))) phase_mechanical_Fp(ph)%data(1:3,1:3,me))))
converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el) converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el)
converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,ip,el) converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,ip,el)
endif endif
@ -1233,7 +1233,7 @@ end function crystallite_stress
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Restore data after homog cutback. !> @brief Restore data after homog cutback.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_restore(ce,includeL) module subroutine mechanical_restore(ce,includeL)
integer, intent(in) :: ce integer, intent(in) :: ce
logical, intent(in) :: & logical, intent(in) :: &
@ -1247,23 +1247,23 @@ module subroutine mech_restore(ce,includeL)
ph = material_phaseAt2(co,ce) ph = material_phaseAt2(co,ce)
me = material_phaseMemberAt2(co,ce) me = material_phaseMemberAt2(co,ce)
if (includeL) then if (includeL) then
constitutive_mech_Lp(ph)%data(1:3,1:3,me) = constitutive_mech_Lp0(ph)%data(1:3,1:3,me) phase_mechanical_Lp(ph)%data(1:3,1:3,me) = phase_mechanical_Lp0(ph)%data(1:3,1:3,me)
constitutive_mech_Li(ph)%data(1:3,1:3,me) = constitutive_mech_Li0(ph)%data(1:3,1:3,me) phase_mechanical_Li(ph)%data(1:3,1:3,me) = phase_mechanical_Li0(ph)%data(1:3,1:3,me)
endif ! maybe protecting everything from overwriting makes more sense endif ! maybe protecting everything from overwriting makes more sense
constitutive_mech_Fp(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) phase_mechanical_Fp(ph)%data(1:3,1:3,me) = phase_mechanical_Fp0(ph)%data(1:3,1:3,me)
constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) phase_mechanical_Fi(ph)%data(1:3,1:3,me) = phase_mechanical_Fi0(ph)%data(1:3,1:3,me)
constitutive_mech_S(ph)%data(1:3,1:3,me) = constitutive_mech_S0(ph)%data(1:3,1:3,me) phase_mechanical_S(ph)%data(1:3,1:3,me) = phase_mechanical_S0(ph)%data(1:3,1:3,me)
plasticState(ph)%state(:,me) = plasticState(ph)%State0(:,me) plasticState(ph)%state(:,me) = plasticState(ph)%State0(:,me)
enddo enddo
end subroutine mech_restore end subroutine mechanical_restore
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculate tangent (dPdF). !> @brief Calculate tangent (dPdF).
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF) module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF)
real(pReal), intent(in) :: dt real(pReal), intent(in) :: dt
integer, intent(in) :: & integer, intent(in) :: &
@ -1297,18 +1297,18 @@ module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF)
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el) me = material_phaseMemberAt(co,ip,el)
call constitutive_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, &
constitutive_mech_Fe(ph)%data(1:3,1:3,me), & phase_mechanical_Fe(ph)%data(1:3,1:3,me), &
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el) phase_mechanical_Fi(ph)%data(1:3,1:3,me),co,ip,el)
call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, & call phase_LiAndItsTangents(devNull,dLidS,dLidFi, &
constitutive_mech_S(ph)%data(1:3,1:3,me), & phase_mechanical_S(ph)%data(1:3,1:3,me), &
constitutive_mech_Fi(ph)%data(1:3,1:3,me), & phase_mechanical_Fi(ph)%data(1:3,1:3,me), &
co,ip,el) co,ip,el)
invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me)) invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me))
invFi = math_inv33(constitutive_mech_Fi(ph)%data(1:3,1:3,me)) invFi = math_inv33(phase_mechanical_Fi(ph)%data(1:3,1:3,me))
invSubFp0 = math_inv33(constitutive_mech_Fp0(ph)%data(1:3,1:3,me)) invSubFp0 = math_inv33(phase_mechanical_Fp0(ph)%data(1:3,1:3,me))
invSubFi0 = math_inv33(constitutive_mech_Fi0(ph)%data(1:3,1:3,me)) invSubFi0 = math_inv33(phase_mechanical_Fi0(ph)%data(1:3,1:3,me))
if (sum(abs(dLidS)) < tol_math_check) then if (sum(abs(dLidS)) < tol_math_check) then
dFidS = 0.0_pReal dFidS = 0.0_pReal
@ -1334,15 +1334,15 @@ module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF)
endif endif
call plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & call plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, &
constitutive_mech_S(ph)%data(1:3,1:3,me), & phase_mechanical_S(ph)%data(1:3,1:3,me), &
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el) phase_mechanical_Fi(ph)%data(1:3,1:3,me),co,ip,el)
dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate dSdF ! calculate dSdF
temp_33_1 = transpose(matmul(invFp,invFi)) temp_33_1 = transpose(matmul(invFp,invFi))
temp_33_2 = matmul(constitutive_mech_F(ph)%data(1:3,1:3,me),invSubFp0) temp_33_2 = matmul(phase_mechanical_F(ph)%data(1:3,1:3,me),invSubFp0)
temp_33_3 = matmul(matmul(constitutive_mech_F(ph)%data(1:3,1:3,me),invFp), invSubFi0) temp_33_3 = matmul(matmul(phase_mechanical_F(ph)%data(1:3,1:3,me),invFp), invSubFi0)
do o=1,3; do p=1,3 do o=1,3; do p=1,3
rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1) rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1)
@ -1370,9 +1370,9 @@ module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! assemble dPdF ! assemble dPdF
temp_33_1 = matmul(constitutive_mech_S(ph)%data(1:3,1:3,me),transpose(invFp)) temp_33_1 = matmul(phase_mechanical_S(ph)%data(1:3,1:3,me),transpose(invFp))
temp_33_2 = matmul(constitutive_mech_F(ph)%data(1:3,1:3,me),invFp) temp_33_2 = matmul(phase_mechanical_F(ph)%data(1:3,1:3,me),invFp)
temp_33_3 = matmul(temp_33_2,constitutive_mech_S(ph)%data(1:3,1:3,me)) temp_33_3 = matmul(temp_33_2,phase_mechanical_S(ph)%data(1:3,1:3,me))
dPdF = 0.0_pReal dPdF = 0.0_pReal
do p=1,3 do p=1,3
@ -1380,129 +1380,129 @@ module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF)
enddo enddo
do o=1,3; do p=1,3 do o=1,3; do p=1,3
dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) & dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) &
+ matmul(matmul(constitutive_mech_F(ph)%data(1:3,1:3,me),dFpinvdF(1:3,1:3,p,o)),temp_33_1) & + matmul(matmul(phase_mechanical_F(ph)%data(1:3,1:3,me),dFpinvdF(1:3,1:3,p,o)),temp_33_1) &
+ matmul(matmul(temp_33_2,dSdF(1:3,1:3,p,o)),transpose(invFp)) & + matmul(matmul(temp_33_2,dSdF(1:3,1:3,p,o)),transpose(invFp)) &
+ matmul(temp_33_3,transpose(dFpinvdF(1:3,1:3,p,o))) + matmul(temp_33_3,transpose(dFpinvdF(1:3,1:3,p,o)))
enddo; enddo enddo; enddo
end function constitutive_mech_dPdF end function phase_mechanical_dPdF
module subroutine mech_restartWrite(groupHandle,ph) module subroutine mechanical_restartWrite(groupHandle,ph)
integer(HID_T), intent(in) :: groupHandle integer(HID_T), intent(in) :: groupHandle
integer, intent(in) :: ph integer, intent(in) :: ph
call HDF5_write(groupHandle,plasticState(ph)%state,'omega') call HDF5_write(groupHandle,plasticState(ph)%state,'omega')
call HDF5_write(groupHandle,constitutive_mech_Fi(ph)%data,'F_i') call HDF5_write(groupHandle,phase_mechanical_Fi(ph)%data,'F_i')
call HDF5_write(groupHandle,constitutive_mech_Li(ph)%data,'L_i') call HDF5_write(groupHandle,phase_mechanical_Li(ph)%data,'L_i')
call HDF5_write(groupHandle,constitutive_mech_Lp(ph)%data,'L_p') call HDF5_write(groupHandle,phase_mechanical_Lp(ph)%data,'L_p')
call HDF5_write(groupHandle,constitutive_mech_Fp(ph)%data,'F_p') call HDF5_write(groupHandle,phase_mechanical_Fp(ph)%data,'F_p')
call HDF5_write(groupHandle,constitutive_mech_S(ph)%data,'S') call HDF5_write(groupHandle,phase_mechanical_S(ph)%data,'S')
call HDF5_write(groupHandle,constitutive_mech_F(ph)%data,'F') call HDF5_write(groupHandle,phase_mechanical_F(ph)%data,'F')
end subroutine mech_restartWrite end subroutine mechanical_restartWrite
module subroutine mech_restartRead(groupHandle,ph) module subroutine mechanical_restartRead(groupHandle,ph)
integer(HID_T), intent(in) :: groupHandle integer(HID_T), intent(in) :: groupHandle
integer, intent(in) :: ph integer, intent(in) :: ph
call HDF5_read(groupHandle,plasticState(ph)%state0,'omega') call HDF5_read(groupHandle,plasticState(ph)%state0,'omega')
call HDF5_read(groupHandle,constitutive_mech_Fi0(ph)%data,'F_i') call HDF5_read(groupHandle,phase_mechanical_Fi0(ph)%data,'F_i')
call HDF5_read(groupHandle,constitutive_mech_Li0(ph)%data,'L_i') call HDF5_read(groupHandle,phase_mechanical_Li0(ph)%data,'L_i')
call HDF5_read(groupHandle,constitutive_mech_Lp0(ph)%data,'L_p') call HDF5_read(groupHandle,phase_mechanical_Lp0(ph)%data,'L_p')
call HDF5_read(groupHandle,constitutive_mech_Fp0(ph)%data,'F_p') call HDF5_read(groupHandle,phase_mechanical_Fp0(ph)%data,'F_p')
call HDF5_read(groupHandle,constitutive_mech_S0(ph)%data,'S') call HDF5_read(groupHandle,phase_mechanical_S0(ph)%data,'S')
call HDF5_read(groupHandle,constitutive_mech_F0(ph)%data,'F') call HDF5_read(groupHandle,phase_mechanical_F0(ph)%data,'F')
end subroutine mech_restartRead end subroutine mechanical_restartRead
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Get first Piola-Kichhoff stress (for use by non-mech physics) !< @brief Get first Piola-Kichhoff stress (for use by non-mech physics)
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module function mech_S(ph,me) result(S) module function mechanical_S(ph,me) result(S)
integer, intent(in) :: ph,me integer, intent(in) :: ph,me
real(pReal), dimension(3,3) :: S real(pReal), dimension(3,3) :: S
S = constitutive_mech_S(ph)%data(1:3,1:3,me) S = phase_mechanical_S(ph)%data(1:3,1:3,me)
end function mech_S end function mechanical_S
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Get plastic velocity gradient (for use by non-mech physics) !< @brief Get plastic velocity gradient (for use by non-mech physics)
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module function mech_L_p(ph,me) result(L_p) module function mechanical_L_p(ph,me) result(L_p)
integer, intent(in) :: ph,me integer, intent(in) :: ph,me
real(pReal), dimension(3,3) :: L_p real(pReal), dimension(3,3) :: L_p
L_p = constitutive_mech_Lp(ph)%data(1:3,1:3,me) L_p = phase_mechanical_Lp(ph)%data(1:3,1:3,me)
end function mech_L_p end function mechanical_L_p
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Get deformation gradient (for use by homogenization) !< @brief Get deformation gradient (for use by homogenization)
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module function constitutive_mech_getF(co,ip,el) result(F) module function phase_mechanical_getF(co,ip,el) result(F)
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ip, el
real(pReal), dimension(3,3) :: F real(pReal), dimension(3,3) :: F
F = constitutive_mech_F(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) F = phase_mechanical_F(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el))
end function constitutive_mech_getF end function phase_mechanical_getF
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Get elastic deformation gradient (for use by non-mech physics) !< @brief Get elastic deformation gradient (for use by non-mech physics)
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module function mech_F_e(ph,me) result(F_e) module function mechanical_F_e(ph,me) result(F_e)
integer, intent(in) :: ph,me integer, intent(in) :: ph,me
real(pReal), dimension(3,3) :: F_e real(pReal), dimension(3,3) :: F_e
F_e = constitutive_mech_Fe(ph)%data(1:3,1:3,me) F_e = phase_mechanical_Fe(ph)%data(1:3,1:3,me)
end function mech_F_e end function mechanical_F_e
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Get second Piola-Kichhoff stress (for use by homogenization) !< @brief Get second Piola-Kichhoff stress (for use by homogenization)
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module function constitutive_mech_getP(co,ip,el) result(P) module function phase_mechanical_getP(co,ip,el) result(P)
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ip, el
real(pReal), dimension(3,3) :: P real(pReal), dimension(3,3) :: P
P = constitutive_mech_P(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) P = phase_mechanical_P(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el))
end function constitutive_mech_getP end function phase_mechanical_getP
! setter for homogenization ! setter for homogenization
module subroutine constitutive_mech_setF(F,co,ip,el) module subroutine phase_mechanical_setF(F,co,ip,el)
real(pReal), dimension(3,3), intent(in) :: F real(pReal), dimension(3,3), intent(in) :: F
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ip, el
constitutive_mech_F(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) = F phase_mechanical_F(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) = F
end subroutine constitutive_mech_setF end subroutine phase_mechanical_setF
end submodule mechanics end submodule mechanics

14
src/phase_mechanics_eigendeformation.f90 → src/phase_mechanical_eigen.f90
View File

@ -126,7 +126,7 @@ end function kinematics_active
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
! ToDo: MD: S is Mi? ! ToDo: MD: S is Mi?
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
S, Fi, co, ip, el) S, Fi, co, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
@ -159,15 +159,15 @@ module subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
dLi_dS = 0.0_pReal dLi_dS = 0.0_pReal
dLi_dFi = 0.0_pReal dLi_dFi = 0.0_pReal
plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) plasticType: select case (phase_plasticity(material_phaseAt(co,el)))
case (PLASTICITY_isotropic_ID) plasticityType case (PLASTICITY_isotropic_ID) plasticType
of = material_phasememberAt(co,ip,el) of = material_phasememberAt(co,ip,el)
instance = phase_plasticityInstance(material_phaseAt(co,el)) instance = phase_plasticInstance(material_phaseAt(co,el))
call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of) call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of)
case default plasticityType case default plasticType
my_Li = 0.0_pReal my_Li = 0.0_pReal
my_dLi_dS = 0.0_pReal my_dLi_dS = 0.0_pReal
end select plasticityType end select plasticType
Li = Li + my_Li Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS dLi_dS = dLi_dS + my_dLi_dS
@ -201,7 +201,7 @@ module subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i)
enddo; enddo enddo; enddo
end subroutine constitutive_LiAndItsTangents end subroutine phase_LiAndItsTangents
end submodule eigendeformation end submodule eigendeformation

0
src/phase_mechanics_eigendeformation_cleavageopening.f90 → src/phase_mechanical_eigen_cleavageopening.f90
View File

0
src/phase_mechanics_eigendeformation_slipplaneopening.f90 → src/phase_mechanical_eigen_slipplaneopening.f90
View File

0
src/phase_mechanics_eigendeformation_thermalexpansion.f90 → src/phase_mechanical_eigen_thermalexpansion.f90
View File

64
src/phase_mechanics_plastic.f90 → src/phase_mechanical_plastic.f90
View File

@ -270,31 +270,31 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
me = material_phasememberAt(co,ip,el) me = material_phasememberAt(co,ip,el)
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) plasticType: select case (phase_plasticity(material_phaseAt(co,el)))
case (PLASTICITY_NONE_ID) plasticityType case (PLASTICITY_NONE_ID) plasticType
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dMp = 0.0_pReal dLp_dMp = 0.0_pReal
case (PLASTICITY_ISOTROPIC_ID) plasticityType case (PLASTICITY_ISOTROPIC_ID) plasticType
call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType case (PLASTICITY_PHENOPOWERLAW_ID) plasticType
call phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) call phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticType
call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticType
call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me,ip,el) call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me,ip,el)
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticType
call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me) call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me)
case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType
call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me) call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me)
end select plasticityType end select plasticType
do i=1,3; do j=1,3 do i=1,3; do j=1,3
dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + & dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + &
@ -323,29 +323,29 @@ module function plastic_dotState(subdt,co,ip,el,ph,me) result(broken)
logical :: broken logical :: broken
Mp = matmul(matmul(transpose(constitutive_mech_Fi(ph)%data(1:3,1:3,me)),& Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,me)),&
constitutive_mech_Fi(ph)%data(1:3,1:3,me)),constitutive_mech_S(ph)%data(1:3,1:3,me)) phase_mechanical_Fi(ph)%data(1:3,1:3,me)),phase_mechanical_S(ph)%data(1:3,1:3,me))
plasticityType: select case (phase_plasticity(ph)) plasticType: select case (phase_plasticity(ph))
case (PLASTICITY_ISOTROPIC_ID) plasticityType case (PLASTICITY_ISOTROPIC_ID) plasticType
call isotropic_dotState(Mp,ph,me) call isotropic_dotState(Mp,ph,me)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType case (PLASTICITY_PHENOPOWERLAW_ID) plasticType
call phenopowerlaw_dotState(Mp,ph,me) call phenopowerlaw_dotState(Mp,ph,me)
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticType
call plastic_kinehardening_dotState(Mp,ph,me) call plastic_kinehardening_dotState(Mp,ph,me)
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticType
call dislotwin_dotState(Mp,thermal_T(ph,me),ph,me) call dislotwin_dotState(Mp,thermal_T(ph,me),ph,me)
case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType
call dislotungsten_dotState(Mp,thermal_T(ph,me),ph,me) call dislotungsten_dotState(Mp,thermal_T(ph,me),ph,me)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticType
call nonlocal_dotState(Mp,thermal_T(ph,me),subdt,ph,me,ip,el) call nonlocal_dotState(Mp,thermal_T(ph,me),subdt,ph,me,ip,el)
end select plasticityType end select plasticType
broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,me))) broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,me)))
@ -369,20 +369,20 @@ module subroutine plastic_dependentState(co, ip, el)
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phasememberAt(co,ip,el) me = material_phasememberAt(co,ip,el)
instance = phase_plasticityInstance(ph) instance = phase_plasticInstance(ph)
plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) plasticType: select case (phase_plasticity(material_phaseAt(co,el)))
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticType
call dislotwin_dependentState(thermal_T(ph,me),instance,me) call dislotwin_dependentState(thermal_T(ph,me),instance,me)
case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType
call dislotungsten_dependentState(instance,me) call dislotungsten_dependentState(instance,me)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticType
call nonlocal_dependentState(instance,me,ip,el) call nonlocal_dependentState(instance,me,ip,el)
end select plasticityType end select plasticType
end subroutine plastic_dependentState end subroutine plastic_dependentState
@ -410,24 +410,24 @@ module function plastic_deltaState(co, ip, el, ph, me) result(broken)
mySize mySize
Mp = matmul(matmul(transpose(constitutive_mech_Fi(ph)%data(1:3,1:3,me)),& Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,me)),&
constitutive_mech_Fi(ph)%data(1:3,1:3,me)),constitutive_mech_S(ph)%data(1:3,1:3,me)) phase_mechanical_Fi(ph)%data(1:3,1:3,me)),phase_mechanical_S(ph)%data(1:3,1:3,me))
instance = phase_plasticityInstance(ph) instance = phase_plasticInstance(ph)
plasticityType: select case (phase_plasticity(ph)) plasticType: select case (phase_plasticity(ph))
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticType
call plastic_kinehardening_deltaState(Mp,instance,me) call plastic_kinehardening_deltaState(Mp,instance,me)
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me))) broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me)))
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticType
call plastic_nonlocal_deltaState(Mp,instance,me,ip,el) call plastic_nonlocal_deltaState(Mp,instance,me,ip,el)
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me))) broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me)))
case default case default
broken = .false. broken = .false.
end select plasticityType end select plasticType
if(.not. broken) then if(.not. broken) then
select case(phase_plasticity(ph)) select case(phase_plasticity(ph))

16
src/phase_mechanics_plastic_dislotungsten.f90 → src/phase_mechanical_plastic_dislotungsten.f90
View File

@ -226,7 +226,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl
sizeState = sizeDotState sizeState = sizeDotState
call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0) call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization
@ -289,16 +289,16 @@ pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, &
integer :: & integer :: &
i,k,l,m,n i,k,l,m,n
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
dot_gamma_pos,dot_gamma_neg, & dot_gamma_pos,dot_gamma_neg, &
ddot_gamma_dtau_pos,ddot_gamma_dtau_neg ddot_gamma_dtau_pos,ddot_gamma_dtau_neg
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dMp = 0.0_pReal dLp_dMp = 0.0_pReal
associate(prm => param(phase_plasticityInstance(ph))) associate(prm => param(phase_plasticInstance(ph)))
call kinetics(Mp,T,phase_plasticityInstance(ph),me,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg) call kinetics(Mp,T,phase_plasticInstance(ph),me,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg)
do i = 1, prm%sum_N_sl do i = 1, prm%sum_N_sl
Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P_sl(1:3,1:3,i) Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P_sl(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -327,7 +327,7 @@ module subroutine dislotungsten_dotState(Mp,T,ph,me)
real(pReal) :: & real(pReal) :: &
VacancyDiffusion VacancyDiffusion
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
gdot_pos, gdot_neg,& gdot_pos, gdot_neg,&
tau_pos,& tau_pos,&
tau_neg, & tau_neg, &
@ -336,10 +336,10 @@ module subroutine dislotungsten_dotState(Mp,T,ph,me)
dot_rho_dip_climb, & dot_rho_dip_climb, &
dip_distance dip_distance
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)),& associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)),&
dot => dotState(phase_plasticityInstance(ph)), dst => dependentState(phase_plasticityInstance(ph))) dot => dotState(phase_plasticInstance(ph)), dst => dependentState(phase_plasticInstance(ph)))
call kinetics(Mp,T,phase_plasticityInstance(ph),me,& call kinetics(Mp,T,phase_plasticInstance(ph),me,&
gdot_pos,gdot_neg, & gdot_pos,gdot_neg, &
tau_pos_out = tau_pos,tau_neg_out = tau_neg) tau_pos_out = tau_pos,tau_neg_out = tau_neg)

36
src/phase_mechanics_plastic_dislotwin.f90 → src/phase_mechanical_plastic_dislotwin.f90
View File

@ -415,7 +415,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
sizeState = sizeDotState sizeState = sizeDotState
call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0) call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and atol ! locally defined state aliases and initialization of state0 and atol
@ -496,8 +496,8 @@ module function plastic_dislotwin_homogenizedC(ph,me) result(homogenizedC)
real(pReal) :: f_unrotated real(pReal) :: f_unrotated
associate(prm => param(phase_plasticityInstance(ph)),& associate(prm => param(phase_plasticInstance(ph)),&
stt => state(phase_plasticityInstance(ph))) stt => state(phase_plasticInstance(ph)))
f_unrotated = 1.0_pReal & f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,me)) & - sum(stt%f_tw(1:prm%sum_N_tw,me)) &
@ -535,11 +535,11 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
BoltzmannRatio, & BoltzmannRatio, &
ddot_gamma_dtau, & ddot_gamma_dtau, &
tau tau
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
dot_gamma_sl,ddot_gamma_dtau_slip dot_gamma_sl,ddot_gamma_dtau_slip
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tw) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_tw) :: &
dot_gamma_twin,ddot_gamma_dtau_twin dot_gamma_twin,ddot_gamma_dtau_twin
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tr) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_tr) :: &
dot_gamma_tr,ddot_gamma_dtau_trans dot_gamma_tr,ddot_gamma_dtau_trans
real(pReal):: dot_gamma_sb real(pReal):: dot_gamma_sb
real(pReal), dimension(3,3) :: eigVectors, P_sb real(pReal), dimension(3,3) :: eigVectors, P_sb
@ -564,7 +564,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
0, 1, 1 & 0, 1, 1 &
],pReal),[ 3,6]) ],pReal),[ 3,6])
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph))) associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)))
f_unrotated = 1.0_pReal & f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,me)) & - sum(stt%f_tw(1:prm%sum_N_tw,me)) &
@ -573,7 +573,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dMp = 0.0_pReal dLp_dMp = 0.0_pReal
call kinetics_slip(Mp,T,phase_plasticityInstance(ph),me,dot_gamma_sl,ddot_gamma_dtau_slip) call kinetics_slip(Mp,T,phase_plasticInstance(ph),me,dot_gamma_sl,ddot_gamma_dtau_slip)
slipContribution: do i = 1, prm%sum_N_sl slipContribution: do i = 1, prm%sum_N_sl
Lp = Lp + dot_gamma_sl(i)*prm%P_sl(1:3,1:3,i) Lp = Lp + dot_gamma_sl(i)*prm%P_sl(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -581,7 +581,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
+ ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) + ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i)
enddo slipContribution enddo slipContribution
call kinetics_twin(Mp,T,dot_gamma_sl,phase_plasticityInstance(ph),me,dot_gamma_twin,ddot_gamma_dtau_twin) call kinetics_twin(Mp,T,dot_gamma_sl,phase_plasticInstance(ph),me,dot_gamma_twin,ddot_gamma_dtau_twin)
twinContibution: do i = 1, prm%sum_N_tw twinContibution: do i = 1, prm%sum_N_tw
Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i) Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -589,7 +589,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
+ ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) + ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
enddo twinContibution enddo twinContibution
call kinetics_trans(Mp,T,dot_gamma_sl,phase_plasticityInstance(ph),me,dot_gamma_tr,ddot_gamma_dtau_trans) call kinetics_trans(Mp,T,dot_gamma_sl,phase_plasticInstance(ph),me,dot_gamma_tr,ddot_gamma_dtau_trans)
transContibution: do i = 1, prm%sum_N_tr transContibution: do i = 1, prm%sum_N_tr
Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i) Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -653,24 +653,24 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
tau, & tau, &
sigma_cl, & !< climb stress sigma_cl, & !< climb stress
b_d !< ratio of Burgers vector to stacking fault width b_d !< ratio of Burgers vector to stacking fault width
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
dot_rho_dip_formation, & dot_rho_dip_formation, &
dot_rho_dip_climb, & dot_rho_dip_climb, &
rho_dip_distance_min, & rho_dip_distance_min, &
dot_gamma_sl dot_gamma_sl
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tw) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_tw) :: &
dot_gamma_twin dot_gamma_twin
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tr) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_tr) :: &
dot_gamma_tr dot_gamma_tr
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)), & associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)), &
dot => dotState(phase_plasticityInstance(ph)), dst => dependentState(phase_plasticityInstance(ph))) dot => dotState(phase_plasticInstance(ph)), dst => dependentState(phase_plasticInstance(ph)))
f_unrotated = 1.0_pReal & f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,me)) & - sum(stt%f_tw(1:prm%sum_N_tw,me)) &
- sum(stt%f_tr(1:prm%sum_N_tr,me)) - sum(stt%f_tr(1:prm%sum_N_tr,me))
call kinetics_slip(Mp,T,phase_plasticityInstance(ph),me,dot_gamma_sl) call kinetics_slip(Mp,T,phase_plasticInstance(ph),me,dot_gamma_sl)
dot%gamma_sl(:,me) = abs(dot_gamma_sl) dot%gamma_sl(:,me) = abs(dot_gamma_sl)
rho_dip_distance_min = prm%D_a*prm%b_sl rho_dip_distance_min = prm%D_a*prm%b_sl
@ -721,10 +721,10 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
- 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,me)*abs(dot_gamma_sl) & - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,me)*abs(dot_gamma_sl) &
- dot_rho_dip_climb - dot_rho_dip_climb
call kinetics_twin(Mp,T,dot_gamma_sl,phase_plasticityInstance(ph),me,dot_gamma_twin) call kinetics_twin(Mp,T,dot_gamma_sl,phase_plasticInstance(ph),me,dot_gamma_twin)
dot%f_tw(:,me) = f_unrotated*dot_gamma_twin/prm%gamma_char dot%f_tw(:,me) = f_unrotated*dot_gamma_twin/prm%gamma_char
call kinetics_trans(Mp,T,dot_gamma_sl,phase_plasticityInstance(ph),me,dot_gamma_tr) call kinetics_trans(Mp,T,dot_gamma_sl,phase_plasticInstance(ph),me,dot_gamma_tr)
dot%f_tr(:,me) = f_unrotated*dot_gamma_tr dot%f_tr(:,me) = f_unrotated*dot_gamma_tr
end associate end associate

8
src/phase_mechanics_plastic_isotropic.f90 → src/phase_mechanical_plastic_isotropic.f90
View File

@ -135,7 +135,7 @@ module function plastic_isotropic_init() result(myPlasticity)
sizeDotState = size(['xi ','gamma']) sizeDotState = size(['xi ','gamma'])
sizeState = sizeDotState sizeState = sizeDotState
call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0) call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization
@ -190,7 +190,7 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
integer :: & integer :: &
k, l, m, n k, l, m, n
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph))) associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)))
Mp_dev = math_deviatoric33(Mp) Mp_dev = math_deviatoric33(Mp)
squarenorm_Mp_dev = math_tensordot(Mp_dev,Mp_dev) squarenorm_Mp_dev = math_tensordot(Mp_dev,Mp_dev)
@ -275,8 +275,8 @@ module subroutine isotropic_dotState(Mp,ph,me)
xi_inf_star, & !< saturation xi xi_inf_star, & !< saturation xi
norm_Mp !< norm of the (deviatoric) Mandel stress norm_Mp !< norm of the (deviatoric) Mandel stress
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)), & associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)), &
dot => dotState(phase_plasticityInstance(ph))) dot => dotState(phase_plasticInstance(ph)))
if (prm%dilatation) then if (prm%dilatation) then
norm_Mp = sqrt(math_tensordot(Mp,Mp)) norm_Mp = sqrt(math_tensordot(Mp,Mp))

16
src/phase_mechanics_plastic_kinehardening.f90 → src/phase_mechanical_plastic_kinehardening.f90
View File

@ -180,7 +180,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names
sizeState = sizeDotState + sizeDeltaState sizeState = sizeDotState + sizeDeltaState
call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,sizeDeltaState) call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,sizeDeltaState)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization
@ -255,16 +255,16 @@ pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
integer :: & integer :: &
i,k,l,m,n i,k,l,m,n
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
gdot_pos,gdot_neg, & gdot_pos,gdot_neg, &
dgdot_dtau_pos,dgdot_dtau_neg dgdot_dtau_pos,dgdot_dtau_neg
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dMp = 0.0_pReal dLp_dMp = 0.0_pReal
associate(prm => param(phase_plasticityInstance(ph))) associate(prm => param(phase_plasticInstance(ph)))
call kinetics(Mp,phase_plasticityInstance(ph),me,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) call kinetics(Mp,phase_plasticInstance(ph),me,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
do i = 1, prm%sum_N_sl do i = 1, prm%sum_N_sl
Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%P(1:3,1:3,i) Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%P(1:3,1:3,i)
@ -292,14 +292,14 @@ module subroutine plastic_kinehardening_dotState(Mp,ph,me)
real(pReal) :: & real(pReal) :: &
sumGamma sumGamma
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
gdot_pos,gdot_neg gdot_pos,gdot_neg
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)),& associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)),&
dot => dotState(phase_plasticityInstance(ph))) dot => dotState(phase_plasticInstance(ph)))
call kinetics(Mp,phase_plasticityInstance(ph),me,gdot_pos,gdot_neg) call kinetics(Mp,phase_plasticInstance(ph),me,gdot_pos,gdot_neg)
dot%accshear(:,me) = abs(gdot_pos+gdot_neg) dot%accshear(:,me) = abs(gdot_pos+gdot_neg)
sumGamma = sum(stt%accshear(:,me)) sumGamma = sum(stt%accshear(:,me))

2
src/phase_mechanics_plastic_none.f90 → src/phase_mechanical_plastic_none.f90
View File

@ -44,7 +44,7 @@ module function plastic_none_init() result(myPlasticity)
phase => phases%get(p) phase => phases%get(p)
if(.not. myPlasticity(p)) cycle if(.not. myPlasticity(p)) cycle
Nconstituents = count(material_phaseAt2 == p) Nconstituents = count(material_phaseAt2 == p)
call constitutive_allocateState(plasticState(p),Nconstituents,0,0,0) call phase_allocateState(plasticState(p),Nconstituents,0,0,0)
enddo enddo
end function plastic_none_init end function plastic_none_init

64
src/phase_mechanics_plastic_nonlocal.f90 → src/phase_mechanical_plastic_nonlocal.f90
View File

@ -407,7 +407,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure
sizeDeltaState = sizeDotState sizeDeltaState = sizeDotState
call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,sizeDeltaState) call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,sizeDeltaState)
plasticState(p)%nonlocal = pl%get_asBool('nonlocal') plasticState(p)%nonlocal = pl%get_asBool('nonlocal')
if(plasticState(p)%nonlocal .and. .not. allocated(IPneighborhood)) & if(plasticState(p)%nonlocal .and. .not. allocated(IPneighborhood)) &
@ -642,8 +642,8 @@ module subroutine nonlocal_dependentState(instance, me, ip, el)
rho0 = getRho0(instance,me,ip,el) rho0 = getRho0(instance,me,ip,el)
if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then
ph = material_phaseAt(1,el) ph = material_phaseAt(1,el)
invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me)) invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me))
invFe = math_inv33(constitutive_mech_Fe(ph)%data(1:3,1:3,me)) invFe = math_inv33(phase_mechanical_Fe(ph)%data(1:3,1:3,me))
rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg) rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg)
rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg) rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg)
@ -662,7 +662,7 @@ module subroutine nonlocal_dependentState(instance, me, ip, el)
neighbor_ip = IPneighborhood(2,n,ip,el) neighbor_ip = IPneighborhood(2,n,ip,el)
no = material_phasememberAt(1,neighbor_ip,neighbor_el) no = material_phasememberAt(1,neighbor_ip,neighbor_el)
if (neighbor_el > 0 .and. neighbor_ip > 0) then if (neighbor_el > 0 .and. neighbor_ip > 0) then
neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el)) neighbor_instance = phase_plasticInstance(material_phaseAt(1,neighbor_el))
if (neighbor_instance == instance) then if (neighbor_instance == instance) then
nRealNeighbors = nRealNeighbors + 1.0_pReal nRealNeighbors = nRealNeighbors + 1.0_pReal
@ -782,25 +782,25 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
l, & l, &
t, & !< dislocation type t, & !< dislocation type
s !< index of my current slip system s !< index of my current slip system
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,8) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,8) :: &
rhoSgl !< single dislocation densities (including blocked) rhoSgl !< single dislocation densities (including blocked)
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,10) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,10) :: &
rho rho
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,4) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,4) :: &
v, & !< velocity v, & !< velocity
tauNS, & !< resolved shear stress including non Schmid and backstress terms tauNS, & !< resolved shear stress including non Schmid and backstress terms
dv_dtau, & !< velocity derivative with respect to the shear stress dv_dtau, & !< velocity derivative with respect to the shear stress
dv_dtauNS !< velocity derivative with respect to the shear stress dv_dtauNS !< velocity derivative with respect to the shear stress
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
tau, & !< resolved shear stress including backstress terms tau, & !< resolved shear stress including backstress terms
gdotTotal !< shear rate gdotTotal !< shear rate
associate(prm => param(phase_plasticityInstance(ph)),dst=>microstructure(phase_plasticityInstance(ph)),& associate(prm => param(phase_plasticInstance(ph)),dst=>microstructure(phase_plasticInstance(ph)),&
stt=>state(phase_plasticityInstance(ph))) stt=>state(phase_plasticInstance(ph)))
ns = prm%sum_N_sl ns = prm%sum_N_sl
!*** shortcut to state variables !*** shortcut to state variables
rho = getRho(phase_plasticityInstance(ph),me,ip,el) rho = getRho(phase_plasticInstance(ph),me,ip,el)
rhoSgl = rho(:,sgl) rhoSgl = rho(:,sgl)
do s = 1,ns do s = 1,ns
@ -820,7 +820,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
! edges ! edges
call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), & call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), &
tau, tauNS(:,1), dst%tau_pass(:,me),1,Temperature, phase_plasticityInstance(ph)) tau, tauNS(:,1), dst%tau_pass(:,me),1,Temperature, phase_plasticInstance(ph))
v(:,2) = v(:,1) v(:,2) = v(:,1)
dv_dtau(:,2) = dv_dtau(:,1) dv_dtau(:,2) = dv_dtau(:,1)
dv_dtauNS(:,2) = dv_dtauNS(:,1) dv_dtauNS(:,2) = dv_dtauNS(:,1)
@ -833,7 +833,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
else else
do t = 3,4 do t = 3,4
call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), & call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), &
tau, tauNS(:,t), dst%tau_pass(:,me),2,Temperature, phase_plasticityInstance(ph)) tau, tauNS(:,t), dst%tau_pass(:,me),2,Temperature, phase_plasticInstance(ph))
enddo enddo
endif endif
@ -992,7 +992,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
c, & !< character of dislocation c, & !< character of dislocation
t, & !< type of dislocation t, & !< type of dislocation
s !< index of my current slip system s !< index of my current slip system
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,10) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,10) :: &
rho, & rho, &
rho0, & !< dislocation density at beginning of time step rho0, & !< dislocation density at beginning of time step
rhoDot, & !< density evolution rhoDot, & !< density evolution
@ -1000,17 +1000,17 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide) rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide)
rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation
rhoDotThermalAnnihilation !< density evolution by thermal annihilation rhoDotThermalAnnihilation !< density evolution by thermal annihilation
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,8) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,8) :: &
rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles)
my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles)
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,4) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,4) :: &
v, & !< current dislocation glide velocity v, & !< current dislocation glide velocity
v0, & v0, &
gdot !< shear rates gdot !< shear rates
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
tau, & !< current resolved shear stress tau, & !< current resolved shear stress
vClimb !< climb velocity of edge dipoles vClimb !< climb velocity of edge dipoles
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,2) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,2) :: &
rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) rhoDip, & !< current dipole dislocation densities (screw and edge dipoles)
dLower, & !< minimum stable dipole distance for edges and screws dLower, & !< minimum stable dipole distance for edges and screws
dUpper !< current maximum stable dipole distance for edges and screws dUpper !< current maximum stable dipole distance for edges and screws
@ -1022,22 +1022,22 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
return return
endif endif
associate(prm => param(phase_plasticityInstance(ph)), & associate(prm => param(phase_plasticInstance(ph)), &
dst => microstructure(phase_plasticityInstance(ph)), & dst => microstructure(phase_plasticInstance(ph)), &
dot => dotState(phase_plasticityInstance(ph)), & dot => dotState(phase_plasticInstance(ph)), &
stt => state(phase_plasticityInstance(ph))) stt => state(phase_plasticInstance(ph)))
ns = prm%sum_N_sl ns = prm%sum_N_sl
tau = 0.0_pReal tau = 0.0_pReal
gdot = 0.0_pReal gdot = 0.0_pReal
rho = getRho(phase_plasticityInstance(ph),me,ip,el) rho = getRho(phase_plasticInstance(ph),me,ip,el)
rhoSgl = rho(:,sgl) rhoSgl = rho(:,sgl)
rhoDip = rho(:,dip) rhoDip = rho(:,dip)
rho0 = getRho0(phase_plasticityInstance(ph),me,ip,el) rho0 = getRho0(phase_plasticInstance(ph),me,ip,el)
my_rhoSgl0 = rho0(:,sgl) my_rhoSgl0 = rho0(:,sgl)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,phase_plasticityInstance(ph)),me) forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,phase_plasticInstance(ph)),me)
gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
#ifdef DEBUG #ifdef DEBUG
@ -1086,7 +1086,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
* sqrt(stt%rho_forest(:,me)) / prm%i_sl / prm%b_sl, 2, 4) * sqrt(stt%rho_forest(:,me)) / prm%i_sl / prm%b_sl, 2, 4)
endif isBCC endif isBCC
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,phase_plasticityInstance(ph)),me) forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,phase_plasticInstance(ph)),me)
!**************************************************************************** !****************************************************************************
@ -1142,7 +1142,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
- rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
- rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have
rhoDot = rhoDotFlux(timestep, phase_plasticityInstance(ph),me,ip,el) & rhoDot = rhoDotFlux(timestep, phase_plasticInstance(ph),me,ip,el) &
+ rhoDotMultiplication & + rhoDotMultiplication &
+ rhoDotSingle2DipoleGlide & + rhoDotSingle2DipoleGlide &
+ rhoDotAthermalAnnihilation & + rhoDotAthermalAnnihilation &
@ -1284,8 +1284,8 @@ function rhoDotFlux(timestep,instance,me,ip,el)
m(1:3,:,3) = -prm%slip_transverse m(1:3,:,3) = -prm%slip_transverse
m(1:3,:,4) = prm%slip_transverse m(1:3,:,4) = prm%slip_transverse
my_F = constitutive_mech_F(ph)%data(1:3,1:3,me) my_F = phase_mechanical_F(ph)%data(1:3,1:3,me)
my_Fe = matmul(my_F, math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me))) my_Fe = matmul(my_F, math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me)))
neighbors: do n = 1,nIPneighbors neighbors: do n = 1,nIPneighbors
@ -1301,9 +1301,9 @@ function rhoDotFlux(timestep,instance,me,ip,el)
opposite_n = IPneighborhood(3,opposite_neighbor,ip,el) opposite_n = IPneighborhood(3,opposite_neighbor,ip,el)
if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient
neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el)) neighbor_instance = phase_plasticInstance(material_phaseAt(1,neighbor_el))
neighbor_F = constitutive_mech_F(np)%data(1:3,1:3,no) neighbor_F = phase_mechanical_F(np)%data(1:3,1:3,no)
neighbor_Fe = matmul(neighbor_F, math_inv33(constitutive_mech_Fp(np)%data(1:3,1:3,no))) neighbor_Fe = matmul(neighbor_F, math_inv33(phase_mechanical_Fp(np)%data(1:3,1:3,no)))
Favg = 0.5_pReal * (my_F + neighbor_F) Favg = 0.5_pReal * (my_F + neighbor_F)
else ! if no neighbor, take my value as average else ! if no neighbor, take my value as average
Favg = my_F Favg = my_F

22
src/phase_mechanics_plastic_phenopowerlaw.f90 → src/phase_mechanical_plastic_phenopowerlaw.f90
View File

@ -231,7 +231,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
sizeState = sizeDotState sizeState = sizeDotState
call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0) call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization
@ -300,18 +300,18 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
integer :: & integer :: &
i,k,l,m,n i,k,l,m,n
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
gdot_slip_pos,gdot_slip_neg, & gdot_slip_pos,gdot_slip_neg, &
dgdot_dtauslip_pos,dgdot_dtauslip_neg dgdot_dtauslip_pos,dgdot_dtauslip_neg
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tw) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_tw) :: &
gdot_twin,dgdot_dtautwin gdot_twin,dgdot_dtautwin
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dMp = 0.0_pReal dLp_dMp = 0.0_pReal
associate(prm => param(phase_plasticityInstance(ph))) associate(prm => param(phase_plasticInstance(ph)))
call kinetics_slip(Mp,phase_plasticityInstance(ph),me,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) call kinetics_slip(Mp,phase_plasticInstance(ph),me,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
slipSystems: do i = 1, prm%sum_N_sl slipSystems: do i = 1, prm%sum_N_sl
Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%P_sl(1:3,1:3,i) Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%P_sl(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -320,7 +320,7 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
+ dgdot_dtauslip_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(m,n,i) + dgdot_dtauslip_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(m,n,i)
enddo slipSystems enddo slipSystems
call kinetics_twin(Mp,phase_plasticityInstance(ph),me,gdot_twin,dgdot_dtautwin) call kinetics_twin(Mp,phase_plasticInstance(ph),me,gdot_twin,dgdot_dtautwin)
twinSystems: do i = 1, prm%sum_N_tw twinSystems: do i = 1, prm%sum_N_tw
Lp = Lp + gdot_twin(i)*prm%P_tw(1:3,1:3,i) Lp = Lp + gdot_twin(i)*prm%P_tw(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -348,12 +348,12 @@ module subroutine phenopowerlaw_dotState(Mp,ph,me)
c_SlipSlip,c_TwinSlip,c_TwinTwin, & c_SlipSlip,c_TwinSlip,c_TwinTwin, &
xi_slip_sat_offset,& xi_slip_sat_offset,&
sumGamma,sumF sumGamma,sumF
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
left_SlipSlip,right_SlipSlip, & left_SlipSlip,right_SlipSlip, &
gdot_slip_pos,gdot_slip_neg gdot_slip_pos,gdot_slip_neg
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)), & associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)), &
dot => dotState(phase_plasticityInstance(ph))) dot => dotState(phase_plasticInstance(ph)))
sumGamma = sum(stt%gamma_slip(:,me)) sumGamma = sum(stt%gamma_slip(:,me))
sumF = sum(stt%gamma_twin(:,me)/prm%gamma_tw_char) sumF = sum(stt%gamma_twin(:,me)/prm%gamma_tw_char)
@ -373,9 +373,9 @@ module subroutine phenopowerlaw_dotState(Mp,ph,me)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! shear rates ! shear rates
call kinetics_slip(Mp,phase_plasticityInstance(ph),me,gdot_slip_pos,gdot_slip_neg) call kinetics_slip(Mp,phase_plasticInstance(ph),me,gdot_slip_pos,gdot_slip_neg)
dot%gamma_slip(:,me) = abs(gdot_slip_pos+gdot_slip_neg) dot%gamma_slip(:,me) = abs(gdot_slip_pos+gdot_slip_neg)
call kinetics_twin(Mp,phase_plasticityInstance(ph),me,dot%gamma_twin(:,me)) call kinetics_twin(Mp,phase_plasticInstance(ph),me,dot%gamma_twin(:,me))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! hardening ! hardening

14
src/phase_thermal.f90
View File

@ -124,7 +124,7 @@ end subroutine thermal_init
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief calculates thermal dissipation rate !< @brief calculates thermal dissipation rate
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module subroutine constitutive_thermal_getRate(TDot, ph,me) module subroutine phase_thermal_getRate(TDot, ph,me)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
@ -153,13 +153,13 @@ module subroutine constitutive_thermal_getRate(TDot, ph,me)
enddo enddo
end subroutine constitutive_thermal_getRate end subroutine phase_thermal_getRate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure !> @brief contains the constitutive equation for calculating the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_thermal_collectDotState(ph,me) result(broken) function phase_thermal_collectDotState(ph,me) result(broken)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
logical :: broken logical :: broken
@ -178,7 +178,7 @@ function constitutive_thermal_collectDotState(ph,me) result(broken)
enddo SourceLoop enddo SourceLoop
end function constitutive_thermal_collectDotState end function phase_thermal_collectDotState
module function thermal_stress(Delta_t,ph,me) result(converged_) module function thermal_stress(Delta_t,ph,me) result(converged_)
@ -207,7 +207,7 @@ function integrateThermalState(Delta_t, ph,me) result(broken)
so, & so, &
sizeDotState sizeDotState
broken = constitutive_thermal_collectDotState(ph,me) broken = phase_thermal_collectDotState(ph,me)
if(broken) return if(broken) return
do so = 1, thermal_Nsources(ph) do so = 1, thermal_Nsources(ph)
@ -264,7 +264,7 @@ end function thermal_dot_T
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Set temperature !< @brief Set temperature
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module subroutine constitutive_thermal_setField(T,dot_T, co,ce) module subroutine phase_thermal_setField(T,dot_T, co,ce)
real(pReal), intent(in) :: T, dot_T real(pReal), intent(in) :: T, dot_T
integer, intent(in) :: ce, co integer, intent(in) :: ce, co
@ -273,7 +273,7 @@ module subroutine constitutive_thermal_setField(T,dot_T, co,ce)
current(material_phaseAt2(co,ce))%T(material_phaseMemberAt2(co,ce)) = T current(material_phaseAt2(co,ce))%T(material_phaseMemberAt2(co,ce)) = T
current(material_phaseAt2(co,ce))%dot_T(material_phaseMemberAt2(co,ce)) = dot_T current(material_phaseAt2(co,ce))%dot_T(material_phaseMemberAt2(co,ce)) = dot_T
end subroutine constitutive_thermal_setField end subroutine phase_thermal_setField

4
src/phase_thermal_dissipation.f90
View File

@ -56,7 +56,7 @@ module function dissipation_init(source_length) result(mySources)
prm%kappa = src%get_asFloat('kappa') prm%kappa = src%get_asFloat('kappa')
Nconstituents = count(material_phaseAt2 == ph) Nconstituents = count(material_phaseAt2 == ph)
call constitutive_allocateState(thermalState(ph)%p(so),Nconstituents,0,0,0) call phase_allocateState(thermalState(ph)%p(so),Nconstituents,0,0,0)
end associate end associate
endif endif
@ -78,7 +78,7 @@ module subroutine dissipation_getRate(TDot, ph,me)
associate(prm => param(ph)) associate(prm => param(ph))
TDot = prm%kappa*sum(abs(mech_S(ph,me)*mech_L_p(ph,me))) TDot = prm%kappa*sum(abs(mechanical_S(ph,me)*mechanical_L_p(ph,me)))
end associate end associate
end subroutine dissipation_getRate end subroutine dissipation_getRate

2
src/phase_thermal_externalheat.f90
View File

@ -69,7 +69,7 @@ module function externalheat_init(source_length) result(mySources)
prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n)) prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n))
Nconstituents = count(material_phaseAt2 == ph) Nconstituents = count(material_phaseAt2 == ph)
call constitutive_allocateState(thermalState(ph)%p(so),Nconstituents,1,1,0) call phase_allocateState(thermalState(ph)%p(so),Nconstituents,1,1,0)
end associate end associate
endif endif
enddo enddo