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

Merge branch 'development' into NoCoreModule

This commit is contained in:
Martin Diehl 2016-06-27 15:35:46 +02:00
parent 91753d5322 b40cbd68bb
commit 259ee072a6
215 changed files with 188793 additions and 192202 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.gitattributes vendored
View File

@ -3,6 +3,8 @@
# always use LF, even if the files are edited on windows, they need to be compiled/used on unix
* text eol=lf
installation/mods_Abaqus/abaqus_v6_windows.env eol=crlf
# Denote all files that are truly binary and should not be modified.
*.png binary
*.jpg binary
*.cae binary

19
DAMASK_env.bat
View File

@ -1,19 +0,0 @@
:: sets up an environment for DAMASK on Windows
:: usage: call DAMASK_env.bat
@echo off
set LOCATION=%~dp0
set DAMASK_ROOT=%LOCATION%\DAMASK
set DAMASK_NUM_THREADS=2
chcp 1252
Title Düsseldorf Advanced Materials Simulation Kit - DAMASK, MPIE Düsseldorf
echo.
echo Düsseldorf Advanced Materials Simulation Kit - DAMASK
echo Max-Planck-Institut für Eisenforschung, Düsseldorf
echo http://damask.mpie.de
echo.
echo Preparing environment ...
echo DAMASK_ROOT=%DAMASK_ROOT%
echo DAMASK_NUM_THREADS=%DAMASK_NUM_THREADS%
set DAMASK_BIN=%DAMASK_ROOT%\bin
set PATH=%PATH%;%DAMASK_BIN%
set PYTHONPATH=%PYTHONPATH%;%DAMASK_ROOT%\lib

2
VERSION
View File

@ -1 +1 @@
v2.0.0-97-g8b27de7
v2.0.0-290-g12ac5e3

6
code/CPFEM2.f90
View File

@ -43,7 +43,7 @@ subroutine CPFEM_initAll(el,ip)
crystallite_init
use homogenization, only: &
homogenization_init, &
materialpoint_postResults
materialpoint_postResults
use IO, only: &
IO_init
use DAMASK_interface
@ -73,7 +73,7 @@ materialpoint_postResults
call crystallite_init
call homogenization_init
call materialpoint_postResults
call CPFEM_init
call CPFEM_init
end subroutine CPFEM_initAll
@ -251,8 +251,6 @@ subroutine CPFEM_general(age, dt)
crystallite_Tstar0_v, &
crystallite_Tstar_v
use homogenization, only: &
materialpoint_F, &
materialpoint_F0, &
materialpoint_stressAndItsTangent, &
materialpoint_postResults
use IO, only: &

28
code/C_routines.c Normal file
View File

@ -0,0 +1,28 @@
/* Unix */
#include <stdio.h>
#include <unistd.h>
#include <dirent.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <stdio.h>
#include <string.h>
/* http://stackoverflow.com/questions/30279228/is-there-an-alternative-to-getcwd-in-fortran-2003-2008 */
void getcurrentworkdir_c(char cwd[], int *stat ){
char cwd_tmp[1024];
if(getcwd(cwd_tmp, sizeof(cwd_tmp)) == cwd_tmp){
strcpy(cwd,cwd_tmp);
*stat = 0;
}
else{
*stat = 1;
}
}
int isdirectory_c(const char *dir){
struct stat statbuf;
if(stat(dir, &statbuf) != 0)
return 0;
return S_ISDIR(statbuf.st_mode);
}

549
code/IO.f90
View File

@ -6,10 +6,6 @@
!> @brief input/output functions, partly depending on chosen solver
!--------------------------------------------------------------------------------------------------
module IO
#ifdef HDF
use hdf5, only: &
HID_T
#endif
use prec, only: &
pInt, &
pReal
@ -18,22 +14,8 @@ module IO
private
character(len=5), parameter, public :: &
IO_EOF = '#EOF#' !< end of file string
#ifdef HDF
integer(HID_T), public, protected :: tempCoordinates, tempResults
integer(HID_T), private :: resultsFile, tempFile
integer(pInt), private :: currentInc
#endif
public :: &
#ifdef HDF
HDF5_mappingConstitutive, &
HDF5_mappingHomogenization, &
HDF5_mappingCells, &
HDF5_addGroup ,&
HDF5_forwardResults, &
HDF5_addScalarDataset, &
IO_formatIntToString ,&
#endif
IO_init, &
IO_read, &
IO_checkAndRewind, &
@ -117,9 +99,6 @@ subroutine IO_init
#include "compilation_info.f90"
endif mainProcess
#ifdef HDF
call HDF5_createJobFile
#endif
end subroutine IO_init
@ -567,7 +546,7 @@ function IO_hybridIA(Nast,ODFfileName)
!--------------------------------------------------------------------------------------------------
! math module is not available
real(pReal), parameter :: PI = 3.14159265358979323846264338327950288419716939937510_pReal
real(pReal), parameter :: PI = 3.141592653589793_pReal
real(pReal), parameter :: INRAD = PI/180.0_pReal
integer(pInt) :: i,j,bin,NnonZero,Nset,Nreps,reps,phi1,Phi,phi2
@ -687,7 +666,7 @@ function IO_hybridIA(Nast,ODFfileName)
else
prob = 0.0_pReal
endif
dV_V(phi2,Phi,phi1) = prob*dg_0*sin((Phi-1.0_pReal+center)*deltas(2))
dV_V(phi2,Phi,phi1) = prob*dg_0*sin((real(Phi-1_pInt,pReal)+center)*deltas(2))
enddo; enddo; enddo
close(FILEUNIT)
dV_V = dV_V/sum_dV_V ! normalize to 1
@ -734,7 +713,7 @@ function IO_hybridIA(Nast,ODFfileName)
do i=1_pInt,Nast
if (i < Nast) then
call random_number(rnd)
j = nint(rnd*(Nreps-i)+i+0.5_pReal,pInt)
j = nint(rnd*real(Nreps-i,pReal)+real(i,pReal)+0.5_pReal,pInt)
else
j = i
endif
@ -1944,526 +1923,4 @@ recursive function abaqus_assembleInputFile(unit1,unit2) result(createSuccess)
end function abaqus_assembleInputFile
#endif
#ifdef HDF
!--------------------------------------------------------------------------------------------------
!> @brief creates and initializes HDF5 output files
!--------------------------------------------------------------------------------------------------
subroutine HDF5_createJobFile
use hdf5
use DAMASK_interface, only: &
getSolverWorkingDirectoryName, &
getSolverJobName
implicit none
integer :: hdferr
integer(SIZE_T) :: typeSize
character(len=1024) :: path
integer(HID_T) :: prp_id
integer(SIZE_T), parameter :: increment = 104857600 ! increase temp file in memory in 100MB steps
!--------------------------------------------------------------------------------------------------
! initialize HDF5 library and check if integer and float type size match
call h5open_f(hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_createJobFile: h5open_f')
call h5tget_size_f(H5T_NATIVE_INTEGER,typeSize, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_createJobFile: h5tget_size_f (int)')
if (int(pInt,SIZE_T)/=typeSize) call IO_error(0_pInt,ext_msg='pInt does not match H5T_NATIVE_INTEGER')
call h5tget_size_f(H5T_NATIVE_DOUBLE,typeSize, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_createJobFile: h5tget_size_f (double)')
if (int(pReal,SIZE_T)/=typeSize) call IO_error(0_pInt,ext_msg='pReal does not match H5T_NATIVE_DOUBLE')
!--------------------------------------------------------------------------------------------------
! open file
path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//'DAMASKout'
call h5fcreate_f(path,H5F_ACC_TRUNC_F,resultsFile,hdferr)
if (hdferr < 0) call IO_error(100_pInt,ext_msg=path)
call HDF5_addStringAttribute(resultsFile,'createdBy','$Id$')
!--------------------------------------------------------------------------------------------------
! open temp file
path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//'DAMASKoutTemp'
call h5pcreate_f(H5P_FILE_ACCESS_F, prp_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_createJobFile: h5pcreate_f')
call h5pset_fapl_core_f(prp_id, increment, .false., hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_createJobFile: h5pset_fapl_core_f')
call h5fcreate_f(path,H5F_ACC_TRUNC_F,tempFile,hdferr)
if (hdferr < 0) call IO_error(100_pInt,ext_msg=path)
!--------------------------------------------------------------------------------------------------
! create mapping groups in out file
call HDF5_closeGroup(HDF5_addGroup("mapping"))
call HDF5_closeGroup(HDF5_addGroup("results"))
call HDF5_closeGroup(HDF5_addGroup("coordinates"))
!--------------------------------------------------------------------------------------------------
! create results group in temp file
tempResults = HDF5_addGroup("results",tempFile)
tempCoordinates = HDF5_addGroup("coordinates",tempFile)
end subroutine HDF5_createJobFile
!--------------------------------------------------------------------------------------------------
!> @brief creates and initializes HDF5 output file
!--------------------------------------------------------------------------------------------------
subroutine HDF5_closeJobFile()
use hdf5
implicit none
integer :: hdferr
call h5fclose_f(resultsFile,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_closeJobFile: h5fclose_f')
end subroutine HDF5_closeJobFile
!--------------------------------------------------------------------------------------------------
!> @brief adds a new group to the results file, or if loc is present at the given location
!--------------------------------------------------------------------------------------------------
integer(HID_T) function HDF5_addGroup(path,loc)
use hdf5
implicit none
character(len=*), intent(in) :: path
integer(HID_T), intent(in),optional :: loc
integer :: hdferr
if (present(loc)) then
call h5gcreate_f(loc, trim(path), HDF5_addGroup, hdferr)
else
call h5gcreate_f(resultsFile, trim(path), HDF5_addGroup, hdferr)
endif
if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_addGroup: h5gcreate_f ('//trim(path)//' )')
end function HDF5_addGroup
!--------------------------------------------------------------------------------------------------
!> @brief adds a new group to the results file
!--------------------------------------------------------------------------------------------------
integer(HID_T) function HDF5_openGroup(path)
use hdf5
implicit none
character(len=*), intent(in) :: path
integer :: hdferr
call h5gopen_f(resultsFile, trim(path), HDF5_openGroup, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_openGroup: h5gopen_f ('//trim(path)//' )')
end function HDF5_openGroup
!--------------------------------------------------------------------------------------------------
!> @brief closes a group
!--------------------------------------------------------------------------------------------------
subroutine HDF5_closeGroup(ID)
use hdf5
implicit none
integer(HID_T), intent(in) :: ID
integer :: hdferr
call h5gclose_f(ID, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_closeGroup: h5gclose_f')
end subroutine HDF5_closeGroup
!--------------------------------------------------------------------------------------------------
!> @brief adds a new group to the results file
!--------------------------------------------------------------------------------------------------
subroutine HDF5_addStringAttribute(entity,attrLabel,attrValue)
use hdf5
implicit none
integer(HID_T), intent(in) :: entity
character(len=*), intent(in) :: attrLabel, attrValue
integer :: hdferr
integer(HID_T) :: attr_id, space_id, type_id
call h5screate_f(H5S_SCALAR_F,space_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5screate_f')
call h5tcopy_f(H5T_NATIVE_CHARACTER, type_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5tcopy_f')
call h5tset_size_f(type_id, int(len(trim(attrValue)),HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5tset_size_f')
call h5acreate_f(entity, trim(attrLabel),type_id,space_id,attr_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5acreate_f')
call h5awrite_f(attr_id, type_id, trim(attrValue), int([1],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5awrite_f')
call h5aclose_f(attr_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5aclose_f')
call h5sclose_f(space_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5sclose_f')
end subroutine HDF5_addStringAttribute
!--------------------------------------------------------------------------------------------------
!> @brief adds the unique mapping from spatial position and constituent ID to results
!--------------------------------------------------------------------------------------------------
subroutine HDF5_mappingConstitutive(mapping)
use hdf5
implicit none
integer(pInt), intent(in), dimension(:,:,:) :: mapping
integer :: hdferr, NmatPoints,Nconstituents
integer(HID_T) :: mapping_id, dtype_id, dset_id, space_id,instance_id,position_id
Nconstituents=size(mapping,1)
NmatPoints=size(mapping,2)
mapping_ID = HDF5_openGroup("mapping")
!--------------------------------------------------------------------------------------------------
! create dataspace
call h5screate_simple_f(2, int([Nconstituents,NmatPoints],HSIZE_T), space_id, hdferr, &
int([Nconstituents,NmatPoints],HSIZE_T))
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive')
!--------------------------------------------------------------------------------------------------
! compound type
call h5tcreate_f(H5T_COMPOUND_F, 6_SIZE_T, dtype_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tcreate_f dtype_id')
call h5tinsert_f(dtype_id, "Constitutive Instance", 0_SIZE_T, H5T_STD_U16LE, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tinsert_f 0')
call h5tinsert_f(dtype_id, "Position in Instance Results", 2_SIZE_T, H5T_STD_U32LE, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tinsert_f 2')
!--------------------------------------------------------------------------------------------------
! create Dataset
call h5dcreate_f(mapping_id, "Constitutive", dtype_id, space_id, dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive')
!--------------------------------------------------------------------------------------------------
! Create memory types (one compound datatype for each member)
call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), instance_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tcreate_f instance_id')
call h5tinsert_f(instance_id, "Constitutive Instance", 0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tinsert_f instance_id')
call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tcreate_f position_id')
call h5tinsert_f(position_id, "Position in Instance Results", 0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tinsert_f position_id')
!--------------------------------------------------------------------------------------------------
! write data by fields in the datatype. Fields order is not important.
call h5dwrite_f(dset_id, position_id, mapping(1:Nconstituents,1:NmatPoints,1), &
int([Nconstituents, NmatPoints],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5dwrite_f position_id')
call h5dwrite_f(dset_id, instance_id, mapping(1:Nconstituents,1:NmatPoints,2), &
int([Nconstituents, NmatPoints],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5dwrite_f instance_id')
!--------------------------------------------------------------------------------------------------
!close types, dataspaces
call h5tclose_f(dtype_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tclose_f dtype_id')
call h5tclose_f(position_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tclose_f position_id')
call h5tclose_f(instance_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tclose_f instance_id')
call h5dclose_f(dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5dclose_f')
call h5sclose_f(space_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5sclose_f')
call HDF5_closeGroup(mapping_ID)
end subroutine HDF5_mappingConstitutive
!--------------------------------------------------------------------------------------------------
!> @brief adds the unique mapping from spatial position and constituent ID to results
!--------------------------------------------------------------------------------------------------
subroutine HDF5_mappingCrystallite(mapping)
use hdf5
implicit none
integer(pInt), intent(in), dimension(:,:,:) :: mapping
integer :: hdferr, NmatPoints,Nconstituents
integer(HID_T) :: mapping_id, dtype_id, dset_id, space_id,instance_id,position_id
Nconstituents=size(mapping,1)
NmatPoints=size(mapping,2)
mapping_ID = HDF5_openGroup("mapping")
!--------------------------------------------------------------------------------------------------
! create dataspace
call h5screate_simple_f(2, int([Nconstituents,NmatPoints],HSIZE_T), space_id, hdferr, &
int([Nconstituents,NmatPoints],HSIZE_T))
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite')
!--------------------------------------------------------------------------------------------------
! compound type
call h5tcreate_f(H5T_COMPOUND_F, 6_SIZE_T, dtype_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tcreate_f dtype_id')
call h5tinsert_f(dtype_id, "Crystallite Instance", 0_SIZE_T, H5T_STD_U16LE, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tinsert_f 0')
call h5tinsert_f(dtype_id, "Position in Instance Results", 2_SIZE_T, H5T_STD_U32LE, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tinsert_f 2')
!--------------------------------------------------------------------------------------------------
! create Dataset
call h5dcreate_f(mapping_id, "Crystallite", dtype_id, space_id, dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite')
!--------------------------------------------------------------------------------------------------
! Create memory types (one compound datatype for each member)
call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), instance_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tcreate_f instance_id')
call h5tinsert_f(instance_id, "Crystallite Instance", 0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tinsert_f instance_id')
call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tcreate_f position_id')
call h5tinsert_f(position_id, "Position in Instance Results", 0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tinsert_f position_id')
!--------------------------------------------------------------------------------------------------
! write data by fields in the datatype. Fields order is not important.
call h5dwrite_f(dset_id, position_id, mapping(1:Nconstituents,1:NmatPoints,1), &
int([Nconstituents, NmatPoints],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5dwrite_f position_id')
call h5dwrite_f(dset_id, instance_id, mapping(1:Nconstituents,1:NmatPoints,2), &
int([Nconstituents, NmatPoints],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5dwrite_f instance_id')
!--------------------------------------------------------------------------------------------------
!close types, dataspaces
call h5tclose_f(dtype_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tclose_f dtype_id')
call h5tclose_f(position_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tclose_f position_id')
call h5tclose_f(instance_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tclose_f instance_id')
call h5dclose_f(dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5dclose_f')
call h5sclose_f(space_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5sclose_f')
call HDF5_closeGroup(mapping_ID)
end subroutine HDF5_mappingCrystallite
!--------------------------------------------------------------------------------------------------
!> @brief adds the unique mapping from spatial position to results
!--------------------------------------------------------------------------------------------------
subroutine HDF5_mappingHomogenization(mapping)
use hdf5
implicit none
integer(pInt), intent(in), dimension(:,:) :: mapping
integer :: hdferr, NmatPoints
integer(HID_T) :: mapping_id, dtype_id, dset_id, space_id,instance_id,position_id,elem_id,ip_id
NmatPoints=size(mapping,1)
mapping_ID = HDF5_openGroup("mapping")
!--------------------------------------------------------------------------------------------------
! create dataspace
call h5screate_simple_f(1, int([NmatPoints],HSIZE_T), space_id, hdferr, &
int([NmatPoints],HSIZE_T))
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization')
!--------------------------------------------------------------------------------------------------
! compound type
call h5tcreate_f(H5T_COMPOUND_F, 11_SIZE_T, dtype_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tcreate_f dtype_id')
call h5tinsert_f(dtype_id, "Homogenization Instance", 0_SIZE_T, H5T_STD_U16LE, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f 0')
call h5tinsert_f(dtype_id, "Position in Instance Results", 2_SIZE_T, H5T_STD_U32LE, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f 2')
call h5tinsert_f(dtype_id, "Element Number", 6_SIZE_T, H5T_STD_U32LE, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f 6')
call h5tinsert_f(dtype_id, "Material Point Number", 10_SIZE_T, H5T_STD_U8LE, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f 10')
!--------------------------------------------------------------------------------------------------
! create Dataset
call h5dcreate_f(mapping_id, "Homogenization", dtype_id, space_id, dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization')
!--------------------------------------------------------------------------------------------------
! Create memory types (one compound datatype for each member)
call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), instance_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tcreate_f instance_id')
call h5tinsert_f(instance_id, "Homogenization Instance", 0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f instance_id')
call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tcreate_f position_id')
call h5tinsert_f(position_id, "Position in Instance Results", 0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f position_id')
call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), elem_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tcreate_f elem_id')
call h5tinsert_f(elem_id, "Element Number", 0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f elem_id')
call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), ip_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tcreate_f ip_id')
call h5tinsert_f(ip_id, "Material Point Number", 0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f ip_id')
!--------------------------------------------------------------------------------------------------
! write data by fields in the datatype. Fields order is not important.
call h5dwrite_f(dset_id, position_id, mapping(1:NmatPoints,1), &
int([NmatPoints],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5dwrite_f position_id')
call h5dwrite_f(dset_id, instance_id, mapping(1:NmatPoints,2), &
int([NmatPoints],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5dwrite_f position_id')
call h5dwrite_f(dset_id, elem_id, mapping(1:NmatPoints,3), &
int([NmatPoints],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5dwrite_f elem_id')
call h5dwrite_f(dset_id, ip_id, mapping(1:NmatPoints,4), &
int([NmatPoints],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5dwrite_f ip_id')
!--------------------------------------------------------------------------------------------------
!close types, dataspaces
call h5tclose_f(dtype_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tclose_f dtype_id')
call h5tclose_f(position_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tclose_f position_id')
call h5tclose_f(instance_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tclose_f instance_id')
call h5tclose_f(ip_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tclose_f ip_id')
call h5tclose_f(elem_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tclose_f elem_id')
call h5dclose_f(dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5dclose_f')
call h5sclose_f(space_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5sclose_f')
call HDF5_closeGroup(mapping_ID)
end subroutine HDF5_mappingHomogenization
!--------------------------------------------------------------------------------------------------
!> @brief adds the unique cell to node mapping
!--------------------------------------------------------------------------------------------------
subroutine HDF5_mappingCells(mapping)
use hdf5
implicit none
integer(pInt), intent(in), dimension(:) :: mapping
integer :: hdferr, Nnodes
integer(HID_T) :: mapping_id, dset_id, space_id
Nnodes=size(mapping)
mapping_ID = HDF5_openGroup("mapping")
!--------------------------------------------------------------------------------------------------
! create dataspace
call h5screate_simple_f(1, int([Nnodes],HSIZE_T), space_id, hdferr, &
int([Nnodes],HSIZE_T))
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCells: h5screate_simple_f')
!--------------------------------------------------------------------------------------------------
! create Dataset
call h5dcreate_f(mapping_id, "Cell",H5T_NATIVE_INTEGER, space_id, dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCells')
!--------------------------------------------------------------------------------------------------
! write data by fields in the datatype. Fields order is not important.
call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER, mapping, int([Nnodes],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCells: h5dwrite_f instance_id')
!--------------------------------------------------------------------------------------------------
!close types, dataspaces
call h5dclose_f(dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5dclose_f')
call h5sclose_f(space_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5sclose_f')
call HDF5_closeGroup(mapping_ID)
end subroutine HDF5_mappingCells
!--------------------------------------------------------------------------------------------------
!> @brief creates a new scalar dataset in the given group location
!--------------------------------------------------------------------------------------------------
subroutine HDF5_addScalarDataset(group,nnodes,label,SIunit)
use hdf5
implicit none
integer(HID_T), intent(in) :: group
integer(pInt), intent(in) :: nnodes
character(len=*), intent(in) :: SIunit,label
integer :: hdferr
integer(HID_T) :: dset_id, space_id
!--------------------------------------------------------------------------------------------------
! create dataspace
call h5screate_simple_f(1, int([Nnodes],HSIZE_T), space_id, hdferr, &
int([Nnodes],HSIZE_T))
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addScalarDataset: h5screate_simple_f')
!--------------------------------------------------------------------------------------------------
! create Dataset
call h5dcreate_f(group, trim(label),H5T_NATIVE_DOUBLE, space_id, dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addScalarDataset: h5dcreate_f')
call HDF5_addStringAttribute(dset_id,'unit',trim(SIunit))
!--------------------------------------------------------------------------------------------------
!close types, dataspaces
call h5dclose_f(dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addScalarDataset: h5dclose_f')
call h5sclose_f(space_id, hdferr)
end subroutine HDF5_addScalarDataset
!--------------------------------------------------------------------------------------------------
!> @brief returns nicely formatted string of integer value
!--------------------------------------------------------------------------------------------------
function IO_formatIntToString(myInt)
implicit none
integer(pInt), intent(in) :: myInt
character(len=1_pInt + int(log10(real(myInt)),pInt)) :: IO_formatIntToString
write(IO_formatIntToString,'('//IO_intOut(myInt)//')') myInt
end function
!--------------------------------------------------------------------------------------------------
!> @brief copies the current temp results to the actual results file
!--------------------------------------------------------------------------------------------------
subroutine HDF5_forwardResults
use hdf5
implicit none
integer :: hdferr
integer(HID_T) :: new_loc_id
new_loc_id = HDF5_openGroup("results")
currentInc = currentInc + 1_pInt
call h5ocopy_f(tempFile, 'results', new_loc_id,dst_name=IO_formatIntToString(currentInc), hdferr=hdferr)
if (hdferr < 0_pInt) call IO_error(1_pInt,ext_msg='HDF5_forwardResults: h5ocopy_f')
call HDF5_closeGroup(new_loc_id)
end subroutine HDF5_forwardResults
#endif
end module IO

65
code/Makefile
View File

@ -7,11 +7,11 @@ SHELL = /bin/sh
# OPTIONS = standard (alternative): meaning
#-------------------------------------------------------------
# F90 = ifort (gfortran): compiler type, choose Intel or GNU
# COMPILERNAME = name of the compiler executable (if not the same as the ype), e.g. using mpich-g90 instead of ifort
# FCOMPILERNAME = name of the compiler executable (if not the same as the ype), e.g. using mpich-g90 instead of ifort
# PORTABLE = TRUE (FALSE): decision, if executable is optimized for the machine on which it was built.
# OPTIMIZATION = DEFENSIVE (OFF,AGGRESSIVE,ULTRA): Optimization mode: O2, O0, O3 + further options for most files, O3 + further options for all files
# OPENMP = TRUE (FALSE): OpenMP multiprocessor support
# PREFIX = arbitrary prefix (before compilername)
# PREFIX = arbitrary prefix (before FCOMPILERNAME)
# OPTION = arbitrary option (just before file to compile)
# SUFFIX = arbitrary suffix (after file to compile)
# STANDARD_CHECK = checking for Fortran 2008, compiler dependend
@ -24,7 +24,8 @@ include ${PETSC_DIR}/lib/petsc/conf/rules
INCLUDE_DIRS := $(PETSC_FC_INCLUDES) -DPETSc -I../lib
LIBRARIES := $(PETSC_WITH_EXTERNAL_LIB)
COMPILERNAME ?= $(FC)
FCOMPILERNAME ?= $(FC)
CCOMPILERNAME ?= $(CC)
LINKERNAME ?= $(FLINKER)
# MPI compiler wrappers will tell if they are pointing to ifort or gfortran
@ -309,7 +310,7 @@ KINEMATICS_FILES = \
kinematics_vacancy_strain.o kinematics_hydrogen_strain.o
PLASTIC_FILES = \
plastic_dislotwin.o plastic_disloUCLA.o plastic_isotropic.o plastic_j2.o \
plastic_dislotwin.o plastic_disloUCLA.o plastic_isotropic.o \
plastic_phenopowerlaw.o plastic_titanmod.o plastic_nonlocal.o plastic_none.o \
plastic_phenoplus.o
@ -350,7 +351,7 @@ DAMASK_spectral.o: INTERFACENAME := spectral_interface.f90
SPECTRAL_SOLVER_FILES = spectral_mech_AL.o spectral_mech_Basic.o spectral_mech_Polarisation.o \
spectral_thermal.o spectral_damage.o
SPECTRAL_FILES = prec.o DAMASK_interface.o IO.o numerics.o debug.o math.o \
SPECTRAL_FILES = C_routines.o system_routines.o prec.o DAMASK_interface.o IO.o libs.o numerics.o debug.o math.o \
FEsolving.o mesh.o material.o lattice.o \
$(SOURCE_FILES) $(KINEMATICS_FILES) $(PLASTIC_FILES) constitutive.o \
crystallite.o \
@ -368,7 +369,7 @@ DAMASK_spectral.exe: DAMASK_spectral.o
DAMASK_spectral.o: DAMASK_spectral.f90 \
$(SPECTRAL_SOLVER_FILES)
$(PREFIX) $(COMPILERNAME) $(COMPILE_MAXOPTI) -c DAMASK_spectral.f90 $(SUFFIX)
$(PREFIX) $(FCOMPILERNAME) $(COMPILE_MAXOPTI) -c DAMASK_spectral.f90 $(SUFFIX)
spectral_mech_AL.o: spectral_mech_AL.f90 \
spectral_utilities.o
@ -400,7 +401,7 @@ DAMASK_FEM.exe: INCLUDE_DIRS += -I./
FEM_SOLVER_FILES = FEM_mech.o FEM_thermal.o FEM_damage.o FEM_vacancyflux.o FEM_porosity.o FEM_hydrogenflux.o
FEM_FILES = prec.o DAMASK_interface.o FEZoo.o IO.o numerics.o debug.o math.o \
FEM_FILES = prec.o DAMASK_interface.o FEZoo.o IO.o libs.o numerics.o debug.o math.o \
FEsolving.o mesh.o material.o lattice.o \
$(SOURCE_FILES) $(KINEMATICS_FILES) $(PLASTIC_FILES) constitutive.o \
crystallite.o \
@ -415,7 +416,7 @@ DAMASK_FEM.exe: DAMASK_FEM_driver.o
$(FEM_FILES) $(LIBRARIES) $(SUFFIX)
DAMASK_FEM_driver.o: DAMASK_FEM_driver.f90 $(FEM_SOLVER_FILES)
$(PREFIX) $(COMPILERNAME) $(COMPILE_MAXOPTI) -c ../private/FEM/code/DAMASK_FEM_driver.f90 $(SUFFIX)
$(PREFIX) $(FCOMPILERNAME) $(COMPILE_MAXOPTI) -c ../private/FEM/code/DAMASK_FEM_driver.f90 $(SUFFIX)
FEM_mech.o: FEM_mech.f90 \
FEM_utilities.o
@ -440,7 +441,7 @@ FEM_utilities.o: FEM_utilities.f90 \
FEZoo.o: $(wildcard FEZoo.f90) \
IO.o
$(IGNORE) $(PREFIX) $(COMPILERNAME) $(COMPILE) -c ../private/FEM/code/FEZoo.f90 $(SUFFIX)
$(IGNORE) $(PREFIX) $(FCOMPILERNAME) $(COMPILE) -c ../private/FEM/code/FEZoo.f90 $(SUFFIX)
touch FEZoo.o
CPFEM.o: CPFEM.f90 \
@ -579,15 +580,12 @@ plastic_phenoplus.o: plastic_phenoplus.f90 \
plastic_isotropic.o: plastic_isotropic.f90 \
lattice.o
plastic_j2.o: plastic_j2.f90 \
lattice.o
plastic_none.o: plastic_none.f90 \
lattice.o
ifeq "$(F90)" "gfortran"
lattice.o: lattice.f90 \
material.o
$(PREFIX) $(COMPILERNAME) $(COMPILE) -ffree-line-length-240 -c lattice.f90 $(SUFFIX)
$(PREFIX) $(FCOMPILERNAME) $(COMPILE) -ffree-line-length-240 -c lattice.f90 $(SUFFIX)
# long lines for interaction matrix
else
lattice.o: lattice.f90 \
@ -602,7 +600,7 @@ mesh.o: mesh.f90 \
FEsolving.o \
math.o \
FEZoo.o
$(PREFIX) $(COMPILERNAME) $(COMPILE) -c $(MESHNAME) -o mesh.o $(SUFFIX)
$(PREFIX) $(FCOMPILERNAME) $(COMPILE) -c $(MESHNAME) -o mesh.o $(SUFFIX)
FEsolving.o: FEsolving.f90 \
debug.o
@ -614,21 +612,23 @@ debug.o: debug.f90 \
numerics.o
numerics.o: numerics.f90 \
libs.o
libs.o: libs.f90 \
IO.o
IO.o: IO.f90 \
DAMASK_interface.o
ifeq "$(F90)" "gfortran"
DAMASK_interface.o: spectral_interface.f90 \
$(wildcard DAMASK_FEM_interface.f90) \
prec.o
$(PREFIX) $(COMPILERNAME) $(COMPILE) -c $(INTERFACENAME) -fall-intrinsics -o DAMASK_interface.o $(SUFFIX)
#-fall-intrinsics: all intrinsic procedures (including the GNU-specific extensions) are accepted. -Wintrinsics-std will be ignored
# and no user-defined procedure with the same name as any intrinsic will be called except when it is explicitly declared external
# --> allows the use of 'getcwd'
prec.o: prec.f90
$(PREFIX) $(COMPILERNAME) $(COMPILE) -c prec.f90 -fno-range-check -fall-intrinsics -fno-fast-math $(SUFFIX)
$(PREFIX) $(FCOMPILERNAME) $(COMPILE) -c $(INTERFACENAME) -o DAMASK_interface.o $(SUFFIX)
ifeq "$(F90)" "gfortran"
prec.o: prec.f90 \
system_routines.o
$(PREFIX) $(FCOMPILERNAME) $(COMPILE) -c prec.f90 -fno-range-check -fall-intrinsics -fno-fast-math $(SUFFIX)
# fno-range-check: Disable range checking on results of simplification of constant expressions during compilation
# --> allows the definition of DAMASK_NaN
#-fall-intrinsics: all intrinsic procedures (including the GNU-specific extensions) are accepted. -Wintrinsics-std will be ignored
@ -637,17 +637,22 @@ prec.o: prec.f90
#-fno-fast-math:
# --> otherwise, when setting -ffast-math, isnan always evaluates to false (I would call it a bug)
else
DAMASK_interface.o: spectral_interface.f90 \
$(wildcard DAMASK_FEM_interface.f90) \
prec.o
$(PREFIX) $(COMPILERNAME) $(COMPILE) -c $(INTERFACENAME) -diag-remark 7410 -stand none -warn nostderrors -o DAMASK_interface.o $(SUFFIX)
# -diag-disable 7410 should disable warning about directory statement in inquire function, but does not work. hence the other 2 statements
prec.o: prec.f90
$(PREFIX) $(COMPILERNAME) $(COMPILE) -c prec.f90 $(SUFFIX)
prec.o: prec.f90 \
system_routines.o
endif
system_routines.o: system_routines.f90 \
C_routines.o
C_routines.o: C_routines.c
%.o : %.f90
$(PREFIX) $(COMPILERNAME) $(COMPILE) -c $< $(SUFFIX)
$(PREFIX) $(FCOMPILERNAME) $(COMPILE) -c $< $(SUFFIX)
%.o : %.c
$(CCOMPILERNAME) -c $<
.PHONY: tidy
tidy:
@ -674,6 +679,6 @@ cleanDAMASK:
.PHONY: help
help:
F90="$(F90)"
COMPILERNAME="$(COMPILERNAME)"
FCOMPILERNAME="$(FCOMPILERNAME)"
COMPILEROUT="$(COMPILEROUT)"

1
code/commercialFEM_fileList.f90
View File

@ -27,7 +27,6 @@
#include "kinematics_hydrogen_strain.f90"
#include "plastic_none.f90"
#include "plastic_isotropic.f90"
#include "plastic_j2.f90"
#include "plastic_phenopowerlaw.f90"
#include "plastic_phenoplus.f90"
#include "plastic_titanmod.f90"

30
code/constitutive.f90
View File

@ -69,7 +69,6 @@ subroutine constitutive_init()
ELASTICITY_hooke_ID, &
PLASTICITY_none_ID, &
PLASTICITY_isotropic_ID, &
PLASTICITY_j2_ID, &
PLASTICITY_phenopowerlaw_ID, &
PLASTICITY_phenoplus_ID, &
PLASTICITY_dislotwin_ID, &
@ -93,7 +92,6 @@ subroutine constitutive_init()
ELASTICITY_HOOKE_label, &
PLASTICITY_NONE_label, &
PLASTICITY_ISOTROPIC_label, &
PLASTICITY_J2_label, &
PLASTICITY_PHENOPOWERLAW_label, &
PLASTICITY_PHENOPLUS_label, &
PLASTICITY_DISLOTWIN_label, &
@ -114,7 +112,6 @@ subroutine constitutive_init()
use plastic_none
use plastic_isotropic
use plastic_j2
use plastic_phenopowerlaw
use plastic_phenoplus
use plastic_dislotwin
@ -160,7 +157,6 @@ subroutine constitutive_init()
! parse plasticities from config file
if (any(phase_plasticity == PLASTICITY_NONE_ID)) call plastic_none_init
if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_J2_ID)) call plastic_j2_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_PHENOPLUS_ID)) call plastic_phenoplus_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init(FILEUNIT)
@ -217,11 +213,6 @@ subroutine constitutive_init()
thisNoutput => plastic_isotropic_Noutput
thisOutput => plastic_isotropic_output
thisSize => plastic_isotropic_sizePostResult
case (PLASTICITY_J2_ID) plasticityType
outputName = PLASTICITY_J2_label
thisNoutput => plastic_j2_Noutput
thisOutput => plastic_j2_output
thisSize => plastic_j2_sizePostResult
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
outputName = PLASTICITY_PHENOPOWERLAW_label
thisNoutput => plastic_phenopowerlaw_Noutput
@ -408,8 +399,6 @@ function constitutive_homogenizedC(ipc,ip,el)
plastic_titanmod_homogenizedC
use plastic_dislotwin, only: &
plastic_dislotwin_homogenizedC
use plastic_disloucla, only: &
plastic_disloucla_homogenizedC
use lattice, only: &
lattice_C66
@ -423,8 +412,6 @@ function constitutive_homogenizedC(ipc,ip,el)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_DISLOTWIN_ID) plasticityType
constitutive_homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el)
case (PLASTICITY_DISLOUCLA_ID) plasticityType
constitutive_homogenizedC = plastic_disloucla_homogenizedC(ipc,ip,el)
case (PLASTICITY_TITANMOD_ID) plasticityType
constitutive_homogenizedC = plastic_titanmod_homogenizedC (ipc,ip,el)
case default plasticityType
@ -513,7 +500,6 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
thermalMapping, &
PLASTICITY_NONE_ID, &
PLASTICITY_ISOTROPIC_ID, &
PLASTICITY_J2_ID, &
PLASTICITY_PHENOPOWERLAW_ID, &
PLASTICITY_PHENOPLUS_ID, &
PLASTICITY_DISLOTWIN_ID, &
@ -522,8 +508,6 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
PLASTICITY_NONLOCAL_ID
use plastic_isotropic, only: &
plastic_isotropic_LpAndItsTangent
use plastic_j2, only: &
plastic_j2_LpAndItsTangent
use plastic_phenopowerlaw, only: &
plastic_phenopowerlaw_LpAndItsTangent
use plastic_phenoplus, only: &
@ -574,8 +558,6 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
dLp_dMstar = 0.0_pReal
case (PLASTICITY_ISOTROPIC_ID) plasticityType
call plastic_isotropic_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
case (PLASTICITY_J2_ID) plasticityType
call plastic_j2_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
case (PLASTICITY_PHENOPLUS_ID) plasticityType
@ -622,8 +604,6 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v
use material, only: &
phase_plasticity, &
material_phase, &
material_homog, &
phaseAt, phasememberAt, &
phase_kinematics, &
phase_Nkinematics, &
PLASTICITY_isotropic_ID, &
@ -903,7 +883,6 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
homogenization_maxNgrains, &
PLASTICITY_none_ID, &
PLASTICITY_isotropic_ID, &
PLASTICITY_j2_ID, &
PLASTICITY_phenopowerlaw_ID, &
PLASTICITY_phenoplus_ID, &
PLASTICITY_dislotwin_ID, &
@ -916,8 +895,6 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
SOURCE_thermal_externalheat_ID
use plastic_isotropic, only: &
plastic_isotropic_dotState
use plastic_j2, only: &
plastic_j2_dotState
use plastic_phenopowerlaw, only: &
plastic_phenopowerlaw_dotState
use plastic_phenoplus, only: &
@ -971,8 +948,6 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_ISOTROPIC_ID) plasticityType
call plastic_isotropic_dotState (Tstar_v,ipc,ip,el)
case (PLASTICITY_J2_ID) plasticityType
call plastic_j2_dotState (Tstar_v,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
case (PLASTICITY_PHENOPLUS_ID) plasticityType
@ -1117,7 +1092,6 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
homogenization_maxNgrains, &
PLASTICITY_NONE_ID, &
PLASTICITY_ISOTROPIC_ID, &
PLASTICITY_J2_ID, &
PLASTICITY_PHENOPOWERLAW_ID, &
PLASTICITY_PHENOPLUS_ID, &
PLASTICITY_DISLOTWIN_ID, &
@ -1130,8 +1104,6 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
SOURCE_damage_anisoDuctile_ID
use plastic_isotropic, only: &
plastic_isotropic_postResults
use plastic_j2, only: &
plastic_j2_postResults
use plastic_phenopowerlaw, only: &
plastic_phenopowerlaw_postResults
use plastic_phenoplus, only: &
@ -1185,8 +1157,6 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
constitutive_postResults(startPos:endPos) = plastic_titanmod_postResults(ipc,ip,el)
case (PLASTICITY_ISOTROPIC_ID) plasticityType
constitutive_postResults(startPos:endPos) = plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
case (PLASTICITY_J2_ID) plasticityType
constitutive_postResults(startPos:endPos) = plastic_j2_postResults(Tstar_v,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
constitutive_postResults(startPos:endPos) = &
plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el)

99
code/crystallite.f90
View File

@ -197,8 +197,6 @@ subroutine crystallite_init
nMax, & !< maximum number of ip neighbors
myNcomponents, & !< number of components at current IP
section = 0_pInt, &
j, &
p, &
mySize
character(len=65536) :: &
@ -258,7 +256,8 @@ subroutine crystallite_init
allocate(crystallite_orientation(4,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_orientation0(4,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_rotation(4,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_disorientation(4,nMax,cMax,iMax,eMax), source=0.0_pReal)
if (any(plasticState%nonLocal)) &
allocate(crystallite_disorientation(4,nMax,cMax,iMax,eMax),source=0.0_pReal)
allocate(crystallite_localPlasticity(cMax,iMax,eMax), source=.true.)
allocate(crystallite_requested(cMax,iMax,eMax), source=.false.)
allocate(crystallite_todo(cMax,iMax,eMax), source=.false.)
@ -510,7 +509,8 @@ end subroutine crystallite_init
!--------------------------------------------------------------------------------------------------
subroutine crystallite_stressAndItsTangent(updateJaco)
use prec, only: &
tol_math_check
tol_math_check, &
dNeq
use numerics, only: &
subStepMinCryst, &
subStepSizeCryst, &
@ -802,7 +802,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
endif
else
subFracIntermediate = maxval(crystallite_subFrac, mask=.not.crystallite_localPlasticity)
if (abs(subFracIntermediate) > tiny(0.0_pReal)) then
if (dNeq(subFracIntermediate,0.0_pReal)) then
crystallite_neighborEnforcedCutback = .false. ! look for ips that require a cutback because of a nonconverged neighbor
!$OMP PARALLEL
!$OMP DO PRIVATE(neighboring_e,neighboring_i)
@ -843,7 +843,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
!$OMP DO PRIVATE(neighboring_e,neighboring_i)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
if (.not. crystallite_localPlasticity(1,i,e) .and. abs(crystallite_subFrac(1,i,e)) > tiny(0.0_pReal)) then
if (.not. crystallite_localPlasticity(1,i,e) .and. dNeq(crystallite_subFrac(1,i,e),0.0_pReal)) then
do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e))))
neighboring_e = mesh_ipNeighborhood(1,n,i,e)
neighboring_i = mesh_ipNeighborhood(2,n,i,e)
@ -3351,7 +3351,8 @@ end subroutine crystallite_integrateStateFPI
!--------------------------------------------------------------------------------------------------
logical function crystallite_stateJump(ipc,ip,el)
use prec, only: &
prec_isNaN
prec_isNaN, &
dNeq
use debug, only: &
debug_level, &
debug_crystallite, &
@ -3403,7 +3404,7 @@ logical function crystallite_stateJump(ipc,ip,el)
enddo
#ifndef _OPENMP
if (any(plasticState(p)%deltaState(1:mySizePlasticDeltaState,c) /= 0.0_pReal) &
if (any(dNeq(plasticState(p)%deltaState(1:mySizePlasticDeltaState,c),0.0_pReal)) &
.and. iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
@ -3458,7 +3459,8 @@ logical function crystallite_integrateStress(&
)
use prec, only: pLongInt, &
tol_math_check, &
prec_isNaN
prec_isNaN, &
dEq
use numerics, only: nStress, &
aTol_crystalliteStress, &
rTol_crystalliteStress, &
@ -3606,7 +3608,7 @@ logical function crystallite_integrateStress(&
!* inversion of Fp_current...
invFp_current = math_inv33(Fp_current)
if (all(abs(invFp_current) <= tiny(0.0_pReal))) then ! math_inv33 returns zero when failed, avoid floating point comparison
failedInversionFp: if (all(dEq(invFp_current,0.0_pReal))) then
#ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fp_current at el (elFE) ip g ',&
@ -3616,13 +3618,13 @@ logical function crystallite_integrateStress(&
endif
#endif
return
endif
endif failedInversionFp
A = math_mul33x33(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp
!* inversion of Fi_current...
invFi_current = math_inv33(Fi_current)
if (all(abs(invFi_current) <= tiny(0.0_pReal))) then ! math_inv33 returns zero when failed, avoid floating point comparison
failedInversionFi: if (all(dEq(invFi_current,0.0_pReal))) then
#ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fi_current at el (elFE) ip ipc ',&
@ -3632,7 +3634,7 @@ logical function crystallite_integrateStress(&
endif
#endif
return
endif
endif failedInversionFi
!* start LpLoop with normal step length
@ -3882,7 +3884,7 @@ logical function crystallite_integrateStress(&
invFp_new = math_mul33x33(invFp_current,B)
invFp_new = invFp_new / math_det33(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize by det
Fp_new = math_inv33(invFp_new)
if (all(abs(Fp_new)<= tiny(0.0_pReal))) then ! math_inv33 returns zero when failed, avoid floating point comparison
failedInversionInvFp: if (all(dEq(Fp_new,0.0_pReal))) then
#ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip ipc ',&
@ -3894,7 +3896,7 @@ logical function crystallite_integrateStress(&
endif
#endif
return
endif
endif failedInversionInvFp
Fe_new = math_mul33x33(math_mul33x33(Fg_new,invFp_new),invFi_new) ! calc resulting Fe
!* calculate 1st Piola-Kirchhoff stress
@ -3961,7 +3963,6 @@ subroutine crystallite_orientations
use plastic_nonlocal, only: &
plastic_nonlocal_updateCompatibility
implicit none
integer(pInt) &
c, & !< counter in integration point component loop
@ -3977,50 +3978,51 @@ subroutine crystallite_orientations
! --- CALCULATE ORIENTATION AND LATTICE ROTATION ---
!$OMP PARALLEL DO PRIVATE(orientation)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do c = 1_pInt,homogenization_Ngrains(mesh_element(3,e))
!$OMP PARALLEL DO PRIVATE(orientation)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do c = 1_pInt,homogenization_Ngrains(mesh_element(3,e))
! somehow this subroutine is not threadsafe, so need critical statement here; not clear, what exactly the problem is
!$OMP CRITICAL (polarDecomp)
orientation = math_RtoQ(transpose(math_rotationalPart33(crystallite_Fe(1:3,1:3,c,i,e)))) ! rotational part from polar decomposition as quaternion
!$OMP END CRITICAL (polarDecomp)
crystallite_rotation(1:4,c,i,e) = lattice_qDisorientation(crystallite_orientation0(1:4,c,i,e), & ! active rotation from ori0
orientation) ! to current orientation (with no symmetry)
crystallite_orientation(1:4,c,i,e) = orientation
enddo; enddo; enddo
!$OMP END PARALLEL DO
!$OMP CRITICAL (polarDecomp)
orientation = math_RtoQ(transpose(math_rotationalPart33(crystallite_Fe(1:3,1:3,c,i,e))))
!$OMP END CRITICAL (polarDecomp)
crystallite_rotation(1:4,c,i,e) = lattice_qDisorientation(crystallite_orientation0(1:4,c,i,e), &! active rotation from initial
orientation) ! to current orientation (with no symmetry)
crystallite_orientation(1:4,c,i,e) = orientation
enddo; enddo; enddo
!$OMP END PARALLEL DO
! --- UPDATE SOME ADDITIONAL VARIABLES THAT ARE NEEDED FOR NONLOCAL MATERIAL ---
! --- we use crystallite_orientation from above, so need a separate loop
!$OMP PARALLEL DO PRIVATE(myPhase,neighboring_e,neighboring_i,neighboringPhase)
nonlocalPresent: if (any(plasticState%nonLocal)) then
!$OMP PARALLEL DO PRIVATE(myPhase,neighboring_e,neighboring_i,neighboringPhase)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
myPhase = material_phase(1,i,e) ! get my phase (non-local models make no sense with more than one grain per material point)
if (plasticState(myPhase)%nonLocal) then ! if nonlocal model
myPhase = material_phase(1,i,e) ! get my phase (non-local models make no sense with more than one grain per material point)
if (plasticState(myPhase)%nonLocal) then ! if nonlocal model
! --- calculate disorientation between me and my neighbor ---
do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e)))) ! loop through my neighbors
do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e)))) ! loop through my neighbors
neighboring_e = mesh_ipNeighborhood(1,n,i,e)
neighboring_i = mesh_ipNeighborhood(2,n,i,e)
if (neighboring_e > 0 .and. neighboring_i > 0) then ! if neighbor exists
neighboringPhase = material_phase(1,neighboring_i,neighboring_e) ! get my neighbor's phase
if (plasticState(neighboringPhase)%nonLocal) then ! neighbor got also nonlocal plasticity
if (lattice_structure(myPhase) == lattice_structure(neighboringPhase)) then ! if my neighbor has same crystal structure like me
if (neighboring_e > 0 .and. neighboring_i > 0) then ! if neighbor exists
neighboringPhase = material_phase(1,neighboring_i,neighboring_e) ! get my neighbor's phase
if (plasticState(neighboringPhase)%nonLocal) then ! neighbor got also nonlocal plasticity
if (lattice_structure(myPhase) == lattice_structure(neighboringPhase)) then ! if my neighbor has same crystal structure like me
crystallite_disorientation(:,n,1,i,e) = &
lattice_qDisorientation( crystallite_orientation(1:4,1,i,e), &
crystallite_orientation(1:4,1,neighboring_i,neighboring_e), &
lattice_structure(myPhase)) ! calculate disorientation for given symmetry
else ! for neighbor with different phase
crystallite_disorientation(:,n,1,i,e) = [0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal] ! 180 degree rotation about 100 axis
lattice_structure(myPhase)) ! calculate disorientation for given symmetry
else ! for neighbor with different phase
crystallite_disorientation(:,n,1,i,e) = [0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal]! 180 degree rotation about 100 axis
endif
else ! for neighbor with local plasticity
crystallite_disorientation(:,n,1,i,e) = [-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal] ! homomorphic identity
else ! for neighbor with local plasticity
crystallite_disorientation(:,n,1,i,e) = [-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal]! homomorphic identity
endif
else ! no existing neighbor
crystallite_disorientation(:,n,1,i,e) = [-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal] ! homomorphic identity
else ! no existing neighbor
crystallite_disorientation(:,n,1,i,e) = [-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal] ! homomorphic identity
endif
enddo
@ -4031,7 +4033,8 @@ subroutine crystallite_orientations
endif
enddo; enddo
!$OMP END PARALLEL DO
!$OMP END PARALLEL DO
endif nonlocalPresent
end subroutine crystallite_orientations
@ -4114,7 +4117,7 @@ function crystallite_postResults(ipc, ip, el)
mySize = 1_pInt
detF = math_det33(crystallite_partionedF(1:3,1:3,ipc,ip,el)) ! V_current = det(F) * V_reference
crystallite_postResults(c+1) = detF * mesh_ipVolume(ip,el) &
/ homogenization_Ngrains(mesh_element(3,el)) ! grain volume (not fraction but absolute)
/ real(homogenization_Ngrains(mesh_element(3,el)),pReal) ! grain volume (not fraction but absolute)
case (orientation_ID)
mySize = 4_pInt
crystallite_postResults(c+1:c+mySize) = crystallite_orientation(1:4,ipc,ip,el) ! grain orientation as quaternion

6
code/damage_local.f90
View File

@ -223,7 +223,7 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el
use material, only: &
homogenization_Ngrains, &
mappingHomogenization, &
phaseAt, phasememberAt, &
phaseAt, &
phase_source, &
phase_Nsources, &
SOURCE_damage_isoBrittle_ID, &
@ -280,8 +280,8 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el
enddo
enddo
phiDot = phiDot/homogenization_Ngrains(mappingHomogenization(2,ip,el))
dPhiDot_dPhi = dPhiDot_dPhi/homogenization_Ngrains(mappingHomogenization(2,ip,el))
phiDot = phiDot/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal)
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal)
end subroutine damage_local_getSourceAndItsTangent

13
code/damage_nonlocal.f90
View File

@ -184,7 +184,7 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip,
use material, only: &
homogenization_Ngrains, &
mappingHomogenization, &
phaseAt, phasememberAt, &
phaseAt, &
phase_source, &
phase_Nsources, &
SOURCE_damage_isoBrittle_ID, &
@ -241,8 +241,8 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip,
enddo
enddo
phiDot = phiDot/homogenization_Ngrains(mappingHomogenization(2,ip,el))
dPhiDot_dPhi = dPhiDot_dPhi/homogenization_Ngrains(mappingHomogenization(2,ip,el))
phiDot = phiDot/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal)
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal)
end subroutine damage_nonlocal_getSourceAndItsTangent
@ -279,9 +279,7 @@ function damage_nonlocal_getDiffusion33(ip,el)
enddo
damage_nonlocal_getDiffusion33 = &
charLength*charLength* &
damage_nonlocal_getDiffusion33/ &
homogenization_Ngrains(homog)
charLength**2_pInt*damage_nonlocal_getDiffusion33/real(homogenization_Ngrains(homog),pReal)
end function damage_nonlocal_getDiffusion33
@ -310,7 +308,8 @@ real(pReal) function damage_nonlocal_getMobility(ip,el)
damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_DamageMobility(material_phase(ipc,ip,el))
enddo
damage_nonlocal_getMobility = damage_nonlocal_getMobility /homogenization_Ngrains(mesh_element(3,el))
damage_nonlocal_getMobility = damage_nonlocal_getMobility/&
real(homogenization_Ngrains(mesh_element(3,el)),pReal)
end function damage_nonlocal_getMobility

9
code/debug.f90
View File

@ -46,10 +46,10 @@ module debug
integer(pInt),protected, dimension(debug_maxNtype+2_pInt), public :: & ! specific ones, and 2 for "all" and "other"
debug_level = 0_pInt
integer(pInt), public :: &
debug_cumLpCalls = 0_pInt, & !< total number of calls to LpAndItsTangent
debug_cumDeltaStateCalls = 0_pInt, & !< total number of calls to deltaState
debug_cumDotStateCalls = 0_pInt !< total number of calls to dotState
integer(pLongInt), public :: &
debug_cumLpCalls = 0_pLongInt, & !< total number of calls to LpAndItsTangent
debug_cumDeltaStateCalls = 0_pLongInt, & !< total number of calls to deltaState
debug_cumDotStateCalls = 0_pLongInt !< total number of calls to dotState
integer(pInt), protected, public :: &
debug_e = 1_pInt, &
@ -67,6 +67,7 @@ module debug
debug_jacobianMaxLocation = 0_pInt, &
debug_jacobianMinLocation = 0_pInt
integer(pInt), dimension(:), allocatable, public :: &
debug_CrystalliteLoopDistribution, & !< distribution of crystallite cutbacks
debug_MaterialpointStateLoopDistribution, &

23
code/homogenization.f90
View File

@ -71,12 +71,6 @@ contains
!> @brief module initialization
!--------------------------------------------------------------------------------------------------
subroutine homogenization_init
#ifdef HDF
use hdf5, only: &
HID_T
use IO, only : &
HDF5_mappingHomogenization
#endif
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use math, only: &
math_I3
@ -131,12 +125,6 @@ subroutine homogenization_init
character(len=64), dimension(:,:), pointer :: thisOutput
character(len=32) :: outputName !< name of output, intermediate fix until HDF5 output is ready
logical :: knownHomogenization, knownThermal, knownDamage, knownVacancyflux, knownPorosity, knownHydrogenflux
#ifdef HDF
integer(pInt), dimension(:,:), allocatable :: mapping
integer(pInt), dimension(:), allocatable :: InstancePosition
allocate(mapping(mesh_ncpelems,4),source=0_pInt)
allocate(InstancePosition(material_Nhomogenization),source=0_pInt)
#endif
!--------------------------------------------------------------------------------------------------
@ -396,17 +384,6 @@ subroutine homogenization_init
!--------------------------------------------------------------------------------------------------
! allocate and initialize global state and postresutls variables
#ifdef HDF
elementLooping: do e = 1,mesh_NcpElems
myInstance = homogenization_typeInstance(mesh_element(3,e))
IpLooping: do i = 1,FE_Nips(FE_geomtype(mesh_element(2,e)))
InstancePosition(myInstance) = InstancePosition(myInstance)+1_pInt
mapping(e,1:4) = [instancePosition(myinstance),myinstance,e,i]
enddo IpLooping
enddo elementLooping
call HDF5_mappingHomogenization(mapping)
#endif
homogenization_maxSizePostResults = 0_pInt
thermal_maxSizePostResults = 0_pInt
damage_maxSizePostResults = 0_pInt

8
code/homogenization_RGC.f90
View File

@ -386,6 +386,8 @@ end subroutine homogenization_RGC_partitionDeformation
! "happy" with result
!--------------------------------------------------------------------------------------------------
function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
use prec, only: &
dEq
use debug, only: &
debug_level, &
debug_homogenization,&
@ -441,10 +443,10 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix
real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax
if(abs(dt) < tiny(0.0_pReal)) then ! zero time step
homogenization_RGC_updateState = .true. ! pretend everything is fine and return
zeroTimeStep: if(dEq(dt,0.0_pReal)) then
homogenization_RGC_updateState = .true. ! pretend everything is fine and return
return
endif
endif zeroTimeStep
!--------------------------------------------------------------------------------------------------
! get the dimension of the cluster (grains and interfaces)

1
code/kinematics_thermal_expansion.f90
View File

@ -77,7 +77,6 @@ subroutine kinematics_thermal_expansion_init(fileUnit)
integer(pInt) :: maxNinstance,phase,instance,kinematics
character(len=65536) :: &
tag = '', &
output = '', &
line = ''
mainProcess: if (worldrank == 0) then

163
code/lattice.f90
View File

@ -17,13 +17,7 @@ module lattice
LATTICE_maxNslipFamily = 13_pInt, & !< max # of slip system families over lattice structures
LATTICE_maxNtwinFamily = 4_pInt, & !< max # of twin system families over lattice structures
LATTICE_maxNtransFamily = 2_pInt, & !< max # of transformation system families over lattice structures
LATTICE_maxNcleavageFamily = 3_pInt, & !< max # of transformation system families over lattice structures
LATTICE_maxNslip = 52_pInt, & !< max # of slip systems over lattice structures
LATTICE_maxNtwin = 24_pInt, & !< max # of twin systems over lattice structures
LATTICE_maxNinteraction = 182_pInt, & !< max # of interaction types (in hardening matrix part)
LATTICE_maxNnonSchmid = 6_pInt, & !< max # of non schmid contributions over lattice structures
LATTICE_maxNtrans = 12_pInt, & !< max # of transformations over lattice structures
LATTICE_maxNcleavage = 9_pInt !< max # of cleavage over lattice structures
LATTICE_maxNcleavageFamily = 3_pInt !< max # of transformation system families over lattice structures
integer(pInt), allocatable, dimension(:,:), protected, public :: &
lattice_NslipSystem, & !< total # of slip systems in each family
@ -80,25 +74,25 @@ module lattice
lattice_NnonSchmid !< total # of non-Schmid contributions for each structure
!--------------------------------------------------------------------------------------------------
! fcc
! face centered cubic
integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: &
LATTICE_fcc_NslipSystem = int([12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],pInt) !< total # of slip systems per family for fcc
LATTICE_fcc_NslipSystem = int([12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],pInt) !< # of slip systems per family for fcc
integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: &
LATTICE_fcc_NtwinSystem = int([12, 0, 0, 0],pInt) !< total # of twin systems per family for fcc
LATTICE_fcc_NtwinSystem = int([12, 0, 0, 0],pInt) !< # of twin systems per family for fcc
integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: &
LATTICE_fcc_NtransSystem = int([12, 0],pInt) !< total # of transformation systems per family for fcc
LATTICE_fcc_NtransSystem = int([12, 0],pInt) !< # of transformation systems per family for fcc
integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: &
LATTICE_fcc_NcleavageSystem = int([3, 4, 0],pInt) !< total # of cleavage systems per family for fcc
LATTICE_fcc_NcleavageSystem = int([3, 4, 0],pInt) !< # of cleavage systems per family for fcc
integer(pInt), parameter, private :: &
LATTICE_fcc_Nslip = 12_pInt, & ! sum(lattice_fcc_NslipSystem), & !< total # of slip systems for fcc
LATTICE_fcc_Ntwin = 12_pInt, & ! sum(lattice_fcc_NtwinSystem) !< total # of twin systems for fcc
LATTICE_fcc_Nslip = 12_pInt, & !sum(lattice_fcc_NslipSystem), & !< total # of slip systems for fcc
LATTICE_fcc_Ntwin = 12_pInt, & !sum(lattice_fcc_NtwinSystem), & !< total # of twin systems for fcc
LATTICE_fcc_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for fcc
LATTICE_fcc_Ntrans = 12_pInt, & !< total # of transformations for fcc
LATTICE_fcc_Ncleavage = 7_pInt !< total # of cleavage systems for fcc
LATTICE_fcc_Ntrans = 12_pInt, & !sum(lattice_fcc_NtransSystem), & !< total # of transformation systems for fcc
LATTICE_fcc_Ncleavage = 7_pInt !sum(lattice_fcc_NcleavageSystem) !< total # of cleavage systems for fcc
real(pReal), dimension(3+3,LATTICE_fcc_Nslip), parameter, private :: &
LATTICE_fcc_systemSlip = reshape(real([&
@ -312,8 +306,8 @@ module lattice
0.0, 0.0, 1.0, 45.0 &
],[ 4_pInt,LATTICE_fcc_Ntrans])
real(pReal), dimension(LATTICE_fcc_Ntrans,LATTICE_fcc_Ntrans), parameter, private :: & ! Matrix for projection of shear from slip system to fault-band (twin) systems
LATTICE_fccTobcc_projectionTrans = reshape(real([& ! For ns = nt = nr
real(pReal), dimension(LATTICE_fcc_Ntrans,LATTICE_fcc_Ntrans), parameter, private :: & ! Matrix for projection of shear from slip system to fault-band (twin) systems
LATTICE_fccTobcc_projectionTrans = reshape(real([& ! For ns = nt = nr
0, 1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
-1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
@ -363,27 +357,26 @@ module lattice
],pReal),[ 3_pInt + 3_pInt,LATTICE_fcc_Ncleavage])
!--------------------------------------------------------------------------------------------------
! bcc
! body centered cubic
integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: &
LATTICE_bcc_NslipSystem = int([ 12, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], pInt) !< total # of slip systems per family for bcc
LATTICE_bcc_NslipSystem = int([ 12, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], pInt) !< # of slip systems per family for bcc
integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: &
LATTICE_bcc_NtwinSystem = int([ 12, 0, 0, 0], pInt) !< total # of twin systems per family for bcc
LATTICE_bcc_NtwinSystem = int([ 12, 0, 0, 0], pInt) !< # of twin systems per family for bcc
integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: &
LATTICE_bcc_NtransSystem = int([0,0],pInt) !< total # of transformation systems per family for bcc
LATTICE_bcc_NtransSystem = int([0,0],pInt) !< # of transformation systems per family for bcc
integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: &
LATTICE_bcc_NcleavageSystem = int([3,6,0],pInt) !< total # of cleavage systems per family for bcc
LATTICE_bcc_NcleavageSystem = int([3,6,0],pInt) !< # of cleavage systems per family for bcc
integer(pInt), parameter, private :: &
LATTICE_bcc_Nslip = 24_pInt, & ! sum(lattice_bcc_NslipSystem), & !< total # of slip systems for bcc
LATTICE_bcc_Ntwin = 12_pInt, & ! sum(lattice_bcc_NtwinSystem) !< total # of twin systems for bcc
LATTICE_bcc_NnonSchmid = 6_pInt, & !< # of non-Schmid contributions for bcc. 6 known non schmid contributions for BCC (A. Koester, A. Ma, A. Hartmaier 2012)
LATTICE_bcc_Ntrans = 0_pInt, & !< total # of transformations for bcc
LATTICE_bcc_Ncleavage = 9_pInt !< total # of cleavage systems for bcc
LATTICE_bcc_Nslip = 24_pInt, & !sum(lattice_bcc_NslipSystem), & !< total # of slip systems for bcc
LATTICE_bcc_Ntwin = 12_pInt, & !sum(lattice_bcc_NtwinSystem), & !< total # of twin systems for bcc
LATTICE_bcc_NnonSchmid = 6_pInt, & !< total # of non-Schmid contributions for bcc (A. Koester, A. Ma, A. Hartmaier 2012)
LATTICE_bcc_Ntrans = 0_pInt, & !sum(lattice_bcc_NtransSystem), & !< total # of transformation systems for bcc
LATTICE_bcc_Ncleavage = 9_pInt !sum(lattice_bcc_NcleavageSystem) !< total # of cleavage systems for bcc
real(pReal), dimension(3+3,LATTICE_bcc_Nslip), parameter, private :: &
LATTICE_bcc_systemSlip = reshape(real([&
! Slip direction Plane normal
@ -561,25 +554,25 @@ module lattice
],pReal),[ 3_pInt + 3_pInt,LATTICE_bcc_Ncleavage])
!--------------------------------------------------------------------------------------------------
! hex
! hexagonal
integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: &
lattice_hex_NslipSystem = int([ 3, 3, 3, 6, 12, 6, 0, 0, 0, 0, 0, 0, 0],pInt) !< # of slip systems per family for hex
lattice_hex_NslipSystem = int([ 3, 3, 3, 6, 12, 6, 0, 0, 0, 0, 0, 0, 0],pInt) !< # of slip systems per family for hex
integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: &
lattice_hex_NtwinSystem = int([ 6, 6, 6, 6],pInt) !< # of slip systems per family for hex
integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: &
LATTICE_hex_NtransSystem = int([0,0],pInt) !< total # of transformation systems per family for hex
LATTICE_hex_NtransSystem = int([0,0],pInt) !< # of transformation systems per family for hex
integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: &
LATTICE_hex_NcleavageSystem = int([3,0,0],pInt) !< total # of cleavage systems per family for hex
LATTICE_hex_NcleavageSystem = int([3,0,0],pInt) !< # of cleavage systems per family for hex
integer(pInt), parameter , private :: &
LATTICE_hex_Nslip = 33_pInt, & ! sum(lattice_hex_NslipSystem), !< total # of slip systems for hex
LATTICE_hex_Ntwin = 24_pInt, & ! sum(lattice_hex_NtwinSystem) !< total # of twin systems for hex
LATTICE_hex_NnonSchmid = 0_pInt, & !< # of non-Schmid contributions for hex
LATTICE_hex_Ntrans = 0_pInt, & !< total # of transformations for hex
LATTICE_hex_Ncleavage = 3_pInt !< total # of transformations for hex
integer(pInt), parameter, private :: &
LATTICE_hex_Nslip = 33_pInt, & !sum(lattice_hex_NslipSystem), & !< total # of slip systems for hex
LATTICE_hex_Ntwin = 24_pInt, & !sum(lattice_hex_NtwinSystem), & !< total # of twin systems for hex
LATTICE_hex_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for hex
LATTICE_hex_Ntrans = 0_pInt, & !sum(lattice_hex_NtransSystem), & !< total # of transformation systems for hex
LATTICE_hex_Ncleavage = 3_pInt !sum(lattice_hex_NcleavageSystem) !< total # of cleavage systems for hex
real(pReal), dimension(4+4,LATTICE_hex_Nslip), parameter, private :: &
LATTICE_hex_systemSlip = reshape(real([&
@ -842,28 +835,26 @@ module lattice
],pReal),[ 4_pInt + 4_pInt,LATTICE_hex_Ncleavage])
!--------------------------------------------------------------------------------------------------
! bct
! body centered tetragonal
integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: &
LATTICE_bct_NslipSystem = int([2, 2, 2, 4, 2, 4, 2, 2, 4, 8, 4, 8, 8 ],pInt) !< # of slip systems per family for bct (Sn) Bieler J. Electr Mater 2009
integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: &
LATTICE_bct_NtwinSystem = int([0, 0, 0, 0], pInt) !< total # of twin systems per family for bct-example
LATTICE_bct_NtwinSystem = int([0, 0, 0, 0], pInt) !< # of twin systems per family for bct
integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: &
LATTICE_bct_NtransSystem = int([0,0],pInt) !< total # of transformation systems per family for bct
LATTICE_bct_NtransSystem = int([0,0],pInt) !< # of transformation systems per family for bct
integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: &
LATTICE_bct_NcleavageSystem = int([0,0,0],pInt) !< total # of cleavage systems per family for bct
LATTICE_bct_NcleavageSystem = int([0,0,0],pInt) !< # of cleavage systems per family for bct
integer(pInt), parameter , private :: &
LATTICE_bct_Nslip = 52_pInt, & ! sum(lattice_bct_NslipSystem), !< total # of slip systems for bct
LATTICE_bct_Ntwin = 0_pInt, & ! sum(lattice_bcc_NtwinSystem) !< total # of twin systems for bct
LATTICE_bct_NnonSchmid = 0_pInt, & !< # of non-Schmid contributions for bct
LATTICE_bct_Ntrans = 0_pInt, & !< total # of transformations for bct
LATTICE_bct_Ncleavage = 0_pInt !< total # of transformations for bct
integer(pInt), parameter, private :: &
LATTICE_bct_Nslip = 52_pInt, & !sum(lattice_bct_NslipSystem), & !< total # of slip systems for bct
LATTICE_bct_Ntwin = 0_pInt, & !sum(lattice_bct_NtwinSystem), & !< total # of twin systems for bct
LATTICE_bct_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for bct
LATTICE_bct_Ntrans = 0_pInt, & !sum(lattice_bct_NtransSystem), & !< total # of transformation systems for bct
LATTICE_bct_Ncleavage = 0_pInt !sum(lattice_bct_NcleavageSystem) !< total # of cleavage systems for bct
real(pReal), dimension(3+3,LATTICE_bct_Nslip), parameter, private :: &
LATTICE_bct_systemSlip = reshape(real([&
@ -1006,12 +997,25 @@ module lattice
!--------------------------------------------------------------------------------------------------
! isotropic
integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: &
LATTICE_iso_NslipSystem = int([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ],pInt) !< # of slip systems per family for iso
integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: &
LATTICE_iso_NtwinSystem = int([0, 0, 0, 0], pInt) !< # of twin systems per family for iso
integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: &
LATTICE_iso_NtransSystem = int([0, 0],pInt) !< # of transformation systems per family for iso
integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: &
LATTICE_iso_NcleavageSystem = int([3,0,0],pInt) !< total # of cleavage systems per family for isotropic
LATTICE_iso_NcleavageSystem = int([3,0,0],pInt) !< # of cleavage systems per family for iso
integer(pInt), parameter, private :: &
LATTICE_iso_Ncleavage = 3_pInt !< total # of cleavage systems for bcc
LATTICE_iso_Nslip = 0_pInt, & !sum(lattice_iso_NslipSystem), & !< total # of slip systems for iso
LATTICE_iso_Ntwin = 0_pInt, & !sum(lattice_iso_NtwinSystem), & !< total # of twin systems for iso
LATTICE_iso_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for iso
LATTICE_iso_Ntrans = 0_pInt, & !sum(lattice_iso_NtransSystem), & !< total # of transformation systems for iso
LATTICE_iso_Ncleavage = 3_pInt !sum(lattice_iso_NcleavageSystem) !< total # of cleavage systems for iso
real(pReal), dimension(3+3,LATTICE_iso_Ncleavage), parameter, private :: &
LATTICE_iso_systemCleavage = reshape(real([&
! Cleavage direction Plane normal
@ -1022,12 +1026,25 @@ module lattice
!--------------------------------------------------------------------------------------------------
! orthorhombic
integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: &
LATTICE_ortho_NslipSystem = int([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ],pInt) !< # of slip systems per family for ortho
integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: &
LATTICE_ortho_NtwinSystem = int([0, 0, 0, 0], pInt) !< # of twin systems per family for ortho
integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: &
LATTICE_ortho_NtransSystem = int([0, 0],pInt) !< # of transformation systems per family for ortho
integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: &
LATTICE_ortho_NcleavageSystem = int([1,1,1],pInt) !< total # of cleavage systems per family for orthotropic
LATTICE_ortho_NcleavageSystem = int([1,1,1],pInt) !< # of cleavage systems per family for ortho
integer(pInt), parameter, private :: &
LATTICE_ortho_Ncleavage = 3_pInt !< total # of cleavage systems for bcc
LATTICE_ortho_Nslip = 0_pInt, & !sum(lattice_ortho_NslipSystem), & !< total # of slip systems for ortho
LATTICE_ortho_Ntwin = 0_pInt, & !sum(lattice_ortho_NtwinSystem), & !< total # of twin systems for ortho
LATTICE_ortho_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for ortho
LATTICE_ortho_Ntrans = 0_pInt, & !sum(lattice_ortho_NtransSystem), & !< total # of transformation systems for ortho
LATTICE_ortho_Ncleavage = 3_pInt !sum(lattice_ortho_NcleavageSystem) !< total # of cleavage systems for ortho
real(pReal), dimension(3+3,LATTICE_ortho_Ncleavage), parameter, private :: &
LATTICE_ortho_systemCleavage = reshape(real([&
! Cleavage direction Plane normal
@ -1036,16 +1053,36 @@ module lattice
1, 0, 0, 0, 0, 1 &
],pReal),[ 3_pInt + 3_pInt,LATTICE_ortho_Ncleavage])
real(pReal), dimension(:,:,:), allocatable, public, protected :: &
integer(pInt), parameter, public :: &
LATTICE_maxNslip = 52_pInt, &
!LATTICE_maxNslip = maxval([LATTICE_fcc_Nslip,LATTICE_bcc_Nslip,LATTICE_hex_Nslip,\
! LATTICE_bct_Nslip,LATTICE_iso_Nslip,LATTICE_ortho_Nslip]), & !< max # of slip systems over lattice structures
LATTICE_maxNtwin = 24_pInt, &
!LATTICE_maxNtwin = maxval([LATTICE_fcc_Ntwin,LATTICE_bcc_Ntwin,LATTICE_hex_Ntwin,\
! LATTICE_bct_Ntwin,LATTICE_iso_Ntwin,LATTICE_ortho_Ntwin]), & !< max # of twin systems over lattice structures
LATTICE_maxNnonSchmid = 6_pInt, &
!LATTICE_maxNtwin = maxval([LATTICE_fcc_NnonSchmid,LATTICE_bcc_NnonSchmid,\
! LATTICE_hex_NnonSchmid,LATTICE_bct_NnonSchmid,\
! LATTICE_iso_NnonSchmid,LATTICE_ortho_NnonSchmid]), & !< max # of non-Schmid contributions over lattice structures
LATTICE_maxNtrans = 12_pInt, &
!LATTICE_maxNtrans = maxval([LATTICE_fcc_Ntrans,LATTICE_bcc_Ntrans,LATTICE_hex_Ntrans,\
! LATTICE_bct_Ntrans,LATTICE_iso_Ntrans,LATTICE_ortho_Ntrans]),&!< max # of transformation systems over lattice structures
LATTICE_maxNcleavage = 9_pInt, &
!LATTICE_maxNcleavage = maxval([LATTICE_fcc_Ncleavage,LATTICE_bcc_Ncleavage,\
! LATTICE_hex_Ncleavage,LATTICE_bct_Ncleavage,\
! LATTICE_iso_Ncleavage,LATTICE_ortho_Ncleavage]) !< max # of cleavage systems over lattice structures
LATTICE_maxNinteraction = 182_pInt !< max # of interaction types (in hardening matrix part)
real(pReal), dimension(:,:,:), allocatable, public, protected :: &
lattice_C66, lattice_trans_C66
real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: &
real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: &
lattice_C3333, lattice_trans_C3333
real(pReal), dimension(:), allocatable, public, protected :: &
real(pReal), dimension(:), allocatable, public, protected :: &
lattice_mu, &
lattice_nu, &
lattice_trans_mu, &
lattice_trans_nu
real(pReal), dimension(:,:,:), allocatable, public, protected :: &
real(pReal), dimension(:,:,:), allocatable, public, protected :: &
lattice_thermalConductivity33, &
lattice_thermalExpansion33, &
lattice_damageDiffusion33, &
@ -1054,7 +1091,7 @@ module lattice
lattice_porosityDiffusion33, &
lattice_hydrogenfluxDiffusion33, &
lattice_hydrogenfluxMobility33
real(pReal), dimension(:), allocatable, public, protected :: &
real(pReal), dimension(:), allocatable, public, protected :: &
lattice_damageMobility, &
lattice_porosityMobility, &
lattice_massDensity, &
@ -1253,7 +1290,7 @@ subroutine lattice_init
endif mainProcess
!--------------------------------------------------------------------------------------------------
! consistency checks
! consistency checks (required since ifort 15.0 does not support sum/maxval in parameter definition)
if (LATTICE_maxNslip /= maxval([LATTICE_fcc_Nslip,LATTICE_bcc_Nslip,LATTICE_hex_Nslip,LATTICE_bct_Nslip])) &
call IO_error(0_pInt,ext_msg = 'LATTICE_maxNslip')

38
code/material.f90
View File

@ -24,7 +24,6 @@ module material
ELASTICITY_hooke_label = 'hooke', &
PLASTICITY_none_label = 'none', &
PLASTICITY_isotropic_label = 'isotropic', &
PLASTICITY_j2_label = 'j2', &
PLASTICITY_phenopowerlaw_label = 'phenopowerlaw', &
PLASTICITY_phenoplus_label = 'phenoplus', &
PLASTICITY_dislotwin_label = 'dislotwin', &
@ -74,7 +73,6 @@ module material
enumerator :: PLASTICITY_undefined_ID, &
PLASTICITY_none_ID, &
PLASTICITY_isotropic_ID, &
PLASTICITY_j2_ID, &
PLASTICITY_phenopowerlaw_ID, &
PLASTICITY_phenoplus_ID, &
PLASTICITY_dislotwin_ID, &
@ -313,7 +311,6 @@ module material
ELASTICITY_hooke_ID ,&
PLASTICITY_none_ID, &
PLASTICITY_isotropic_ID, &
PLASTICITY_J2_ID, &
PLASTICITY_phenopowerlaw_ID, &
PLASTICITY_phenoplus_ID, &
PLASTICITY_dislotwin_ID, &
@ -351,9 +348,6 @@ module material
HYDROGENFLUX_cahnhilliard_ID, &
HOMOGENIZATION_none_ID, &
HOMOGENIZATION_isostrain_ID, &
#ifdef HDF
material_NconstituentsPhase, &
#endif
HOMOGENIZATION_RGC_ID
private :: &
@ -982,8 +976,6 @@ subroutine material_parsePhase(fileUnit,myPart)
phase_plasticity(section) = PLASTICITY_NONE_ID
case (PLASTICITY_ISOTROPIC_label)
phase_plasticity(section) = PLASTICITY_ISOTROPIC_ID
case (PLASTICITY_J2_label)
phase_plasticity(section) = PLASTICITY_J2_ID
case (PLASTICITY_PHENOPOWERLAW_label)
phase_plasticity(section) = PLASTICITY_PHENOPOWERLAW_ID
case (PLASTICITY_PHENOPLUS_label)
@ -1240,6 +1232,8 @@ end subroutine material_parseTexture
!! calculates the volume of the grains and deals with texture components and hybridIA
!--------------------------------------------------------------------------------------------------
subroutine material_populateGrains
use prec, only: &
dEq
use math, only: &
math_RtoEuler, &
math_EulerToR, &
@ -1379,7 +1373,7 @@ subroutine material_populateGrains
else
forall (i = 1_pInt:FE_Nips(t)) & ! loop over IPs
volumeOfGrain(grain+(i-1)*dGrains+1_pInt:grain+i*dGrains) = &
mesh_ipVolume(i,e)/dGrains ! assign IPvolume/Ngrains@IP to all grains of IP
mesh_ipVolume(i,e)/real(dGrains,pReal) ! assign IPvolume/Ngrains@IP to all grains of IP
grain = grain + FE_Nips(t) * dGrains ! wind forward by Nips*Ngrains@IP
endif
enddo
@ -1401,7 +1395,7 @@ subroutine material_populateGrains
NgrainsOfConstituent = 0_pInt ! reset counter of grains per constituent
forall (i = 1_pInt:myNconstituents) &
NgrainsOfConstituent(i) = nint(microstructure_fraction(i,micro) * myNgrains, pInt) ! do rounding integer conversion
NgrainsOfConstituent(i) = nint(microstructure_fraction(i,micro)*real(myNgrains,pReal),pInt)! do rounding integer conversion
do while (sum(NgrainsOfConstituent) /= myNgrains) ! total grain count over constituents wrong?
sgn = sign(1_pInt, myNgrains - sum(NgrainsOfConstituent)) ! direction of required change
extreme = 0.0_pReal
@ -1442,24 +1436,24 @@ subroutine material_populateGrains
! ...has texture components
if (texture_ODFfile(textureID) == '') then
gauss: do t = 1_pInt,texture_Ngauss(textureID) ! loop over Gauss components
do g = 1_pInt,int(myNorientations*texture_Gauss(5,t,textureID),pInt) ! loop over required grain count
do g = 1_pInt,int(real(myNorientations,pReal)*texture_Gauss(5,t,textureID),pInt) ! loop over required grain count
orientationOfGrain(:,grain+constituentGrain+g) = &
math_sampleGaussOri(texture_Gauss(1:3,t,textureID),&
texture_Gauss( 4,t,textureID))
enddo
constituentGrain = &
constituentGrain + int(myNorientations*texture_Gauss(5,t,textureID)) ! advance counter for grains of current constituent
constituentGrain + int(real(myNorientations,pReal)*texture_Gauss(5,t,textureID)) ! advance counter for grains of current constituent
enddo gauss
fiber: do t = 1_pInt,texture_Nfiber(textureID) ! loop over fiber components
do g = 1_pInt,int(myNorientations*texture_Fiber(6,t,textureID),pInt) ! loop over required grain count
do g = 1_pInt,int(real(myNorientations,pReal)*texture_Fiber(6,t,textureID),pInt) ! loop over required grain count
orientationOfGrain(:,grain+constituentGrain+g) = &
math_sampleFiberOri(texture_Fiber(1:2,t,textureID),&
texture_Fiber(3:4,t,textureID),&
texture_Fiber( 5,t,textureID))
enddo
constituentGrain = &
constituentGrain + int(myNorientations*texture_fiber(6,t,textureID),pInt) ! advance counter for grains of current constituent
constituentGrain + int(real(myNorientations,pReal)*texture_fiber(6,t,textureID),pInt) ! advance counter for grains of current constituent
enddo fiber
random: do constituentGrain = constituentGrain+1_pInt,myNorientations ! fill remainder with random
@ -1470,7 +1464,7 @@ subroutine material_populateGrains
else
orientationOfGrain(1:3,grain+1_pInt:grain+myNorientations) = &
IO_hybridIA(myNorientations,texture_ODFfile(textureID))
if (all(orientationOfGrain(1:3,grain+1_pInt) == -1.0_pReal)) call IO_error(156_pInt)
if (all(dEq(orientationOfGrain(1:3,grain+1_pInt),-1.0_pReal))) call IO_error(156_pInt)
endif
!--------------------------------------------------------------------------------------------------
@ -1509,7 +1503,7 @@ subroutine material_populateGrains
do j = 1_pInt,NgrainsOfConstituent(i)-1_pInt ! walk thru grains of current constituent
call random_number(rnd)
t = nint(rnd*(NgrainsOfConstituent(i)-j)+j+0.5_pReal,pInt) ! select a grain in remaining list
t = nint(rnd*real(NgrainsOfConstituent(i)-j,pReal)+real(j,pReal)+0.5_pReal,pInt) ! select a grain in remaining list
m = phaseOfGrain(grain+t) ! exchange current with random
phaseOfGrain(grain+t) = phaseOfGrain(grain+j)
phaseOfGrain(grain+j) = m
@ -1547,7 +1541,7 @@ subroutine material_populateGrains
randomOrder = math_range(microstructure_maxNconstituents) ! start out with ordered sequence of constituents
call random_number(rndArray) ! as many rnd numbers as (max) constituents
do j = 1_pInt, myNconstituents - 1_pInt ! loop over constituents ...
r = nint(rndArray(j)*(myNconstituents-j)+j+0.5_pReal,pInt) ! ... select one in remaining list
r = nint(rndArray(j)*real(myNconstituents-j,pReal)+real(j,pReal)+0.5_pReal,pInt) ! ... select one in remaining list
c = randomOrder(r) ! ... call it "c"
randomOrder(r) = randomOrder(j) ! ... and exchange with present position in constituent list
grain = sum(NgrainsOfConstituent(1:c-1_pInt)) ! figure out actual starting index in overall/consecutive grain population
@ -1603,14 +1597,4 @@ subroutine material_populateGrains
end subroutine material_populateGrains
#ifdef HDF
integer(pInt) pure function material_NconstituentsPhase(matID)
implicit none
integer(pInt), intent(in) :: matID
material_NconstituentsPhase = count(microstructure_phase == matID)
end function
#endif
end module material

36
code/math.f90
View File

@ -13,10 +13,10 @@ module math
implicit none
private
real(pReal), parameter, public :: PI = 3.14159265358979323846264338327950288419716939937510_pReal !< ratio of a circle's circumference to its diameter
real(pReal), parameter, public :: PI = 3.141592653589793_pReal !< ratio of a circle's circumference to its diameter
real(pReal), parameter, public :: INDEG = 180.0_pReal/PI !< conversion from radian into degree
real(pReal), parameter, public :: INRAD = PI/180.0_pReal !< conversion from degree into radian
complex(pReal), parameter, public :: TWOPIIMG = (0.0_pReal,2.0_pReal)* PI !< Re(0.0), Im(2xPi)
complex(pReal), parameter, public :: TWOPIIMG = (0.0_pReal,2.0_pReal)*(PI,0.0_pReal) !< Re(0.0), Im(2xPi)
real(pReal), dimension(3,3), parameter, public :: &
MATH_I3 = reshape([&
@ -704,6 +704,8 @@ end function math_transpose33
! returns all zeroes if not possible, i.e. if det close to zero
!--------------------------------------------------------------------------------------------------
pure function math_inv33(A)
use prec, only: &
dNeq
implicit none
real(pReal),dimension(3,3),intent(in) :: A
@ -716,7 +718,7 @@ pure function math_inv33(A)
DetA = A(1,1) * math_inv33(1,1) + A(1,2) * math_inv33(2,1) + A(1,3) * math_inv33(3,1)
if (abs(DetA) > tiny(DetA)) then ! use a real threshold here
if (dNeq(DetA,0.0_pReal)) then
math_inv33(1,2) = -A(1,2) * A(3,3) + A(1,3) * A(3,2)
math_inv33(2,2) = A(1,1) * A(3,3) - A(1,3) * A(3,1)
math_inv33(3,2) = -A(1,1) * A(3,2) + A(1,2) * A(3,1)
@ -740,6 +742,8 @@ end function math_inv33
! returns error if not possible, i.e. if det close to zero
!--------------------------------------------------------------------------------------------------
pure subroutine math_invert33(A, InvA, DetA, error)
use prec, only: &
dEq
implicit none
logical, intent(out) :: error
@ -753,7 +757,7 @@ pure subroutine math_invert33(A, InvA, DetA, error)
DetA = A(1,1) * InvA(1,1) + A(1,2) * InvA(2,1) + A(1,3) * InvA(3,1)
if (abs(DetA) <= tiny(DetA)) then
if (dEq(DetA,0.0_pReal)) then
InvA = 0.0_pReal
error = .true.
else
@ -820,6 +824,9 @@ subroutine math_invert(myDim,A, InvA, error)
real(pReal), dimension(myDim,myDim), intent(out) :: invA
logical, intent(out) :: error
external :: &
dgetrf, &
dgetri
invA = A
call dgetrf(myDim,myDim,invA,myDim,ipiv,ierr)
@ -1255,6 +1262,8 @@ end function math_qNorm
!> @brief quaternion inversion
!--------------------------------------------------------------------------------------------------
pure function math_qInv(Q)
use prec, only: &
dNeq
implicit none
real(pReal), dimension(4), intent(in) :: Q
@ -1264,8 +1273,7 @@ pure function math_qInv(Q)
math_qInv = 0.0_pReal
squareNorm = math_qDot(Q,Q)
if (abs(squareNorm) > tiny(squareNorm)) &
math_qInv = math_qConj(Q) / squareNorm
if (dNeq(squareNorm,0.0_pReal)) math_qInv = math_qConj(Q) / squareNorm
end function math_qInv
@ -1626,14 +1634,14 @@ pure function math_qToAxisAngle(Q)
real(pReal) :: halfAngle, sinHalfAngle
real(pReal), dimension(4) :: math_qToAxisAngle
halfAngle = acos(max(-1.0_pReal, min(1.0_pReal, Q(1)))) ! limit to [-1,1] --> 0 to 180 deg
halfAngle = acos(math_limit(Q(1),-1.0_pReal,1.0_pReal))
sinHalfAngle = sin(halfAngle)
if (sinHalfAngle <= 1.0e-4_pReal) then ! very small rotation angle?
smallRotation: if (sinHalfAngle <= 1.0e-4_pReal) then
math_qToAxisAngle = 0.0_pReal
else
else smallRotation
math_qToAxisAngle= [ Q(2:4)/sinHalfAngle, halfAngle*2.0_pReal]
endif
endif smallRotation
end function math_qToAxisAngle
@ -2069,6 +2077,8 @@ end function math_eigenvectorBasisSym33
!> @brief rotational part from polar decomposition of 33 tensor m
!--------------------------------------------------------------------------------------------------
function math_rotationalPart33(m)
use prec, only: &
dEq
use IO, only: &
IO_warning
@ -2080,12 +2090,12 @@ function math_rotationalPart33(m)
U = math_eigenvectorBasisSym33(math_mul33x33(transpose(m),m))
Uinv = math_inv33(U)
if (all(abs(Uinv) <= tiny(Uinv))) then ! math_inv33 returns zero when failed, avoid floating point equality comparison
inversionFailed: if (all(dEq(Uinv,0.0_pReal))) then
math_rotationalPart33 = math_I3
call IO_warning(650_pInt)
else
else inversionFailed
math_rotationalPart33 = math_mul33x33(m,Uinv)
endif
endif inversionFailed
end function math_rotationalPart33

835
code/plastic_disloUCLA.f90
View File

File diff suppressed because it is too large Load Diff

41
code/plastic_dislotwin.f90
View File

@ -199,6 +199,9 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine plastic_dislotwin_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: &
dEq, &
dNeq
use debug, only: &
debug_level,&
debug_constitutive,&
@ -749,8 +752,8 @@ subroutine plastic_dislotwin_init(fileUnit)
if (plastic_dislotwin_Qsd(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='Qsd ('//PLASTICITY_DISLOTWIN_label//')')
if (sum(plastic_dislotwin_Ntwin(:,instance)) > 0_pInt) then
if (abs(plastic_dislotwin_SFE_0K(instance)) <= tiny(0.0_pReal) .and. &
abs(plastic_dislotwin_dSFE_dT(instance)) <= tiny(0.0_pReal) .and. &
if (dEq(plastic_dislotwin_SFE_0K(instance), 0.0_pReal) .and. &
dEq(plastic_dislotwin_dSFE_dT(instance),0.0_pReal) .and. &
lattice_structure(phase) == LATTICE_fcc_ID) &
call IO_error(211_pInt,el=instance,ext_msg='SFE0K ('//PLASTICITY_DISLOTWIN_label//')')
if (plastic_dislotwin_aTolRho(instance) <= 0.0_pReal) &
@ -759,8 +762,8 @@ subroutine plastic_dislotwin_init(fileUnit)
call IO_error(211_pInt,el=instance,ext_msg='aTolTwinFrac ('//PLASTICITY_DISLOTWIN_label//')')
endif
if (sum(plastic_dislotwin_Ntrans(:,instance)) > 0_pInt) then
if (abs(plastic_dislotwin_SFE_0K(instance)) <= tiny(0.0_pReal) .and. &
abs(plastic_dislotwin_dSFE_dT(instance)) <= tiny(0.0_pReal) .and. &
if (dEq(plastic_dislotwin_SFE_0K(instance), 0.0_pReal) .and. &
dEq(plastic_dislotwin_dSFE_dT(instance),0.0_pReal) .and. &
lattice_structure(phase) == LATTICE_fcc_ID) &
call IO_error(211_pInt,el=instance,ext_msg='SFE0K ('//PLASTICITY_DISLOTWIN_label//')')
if (plastic_dislotwin_aTolTransFrac(instance) <= 0.0_pReal) &
@ -773,8 +776,8 @@ subroutine plastic_dislotwin_init(fileUnit)
if (plastic_dislotwin_sbVelocity(instance) > 0.0_pReal .and. &
plastic_dislotwin_pShearBand(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='pShearBand ('//PLASTICITY_DISLOTWIN_label//')')
if (abs(plastic_dislotwin_dipoleFormationFactor(instance)) > tiny(0.0_pReal) .and. &
plastic_dislotwin_dipoleFormationFactor(instance) /= 1.0_pReal) &
if (dNeq(plastic_dislotwin_dipoleFormationFactor(instance), 0.0_pReal) .and. &
dNeq(plastic_dislotwin_dipoleFormationFactor(instance), 1.0_pReal)) &
call IO_error(211_pInt,el=instance,ext_msg='dipoleFormationFactor ('//PLASTICITY_DISLOTWIN_label//')')
if (plastic_dislotwin_sbVelocity(instance) > 0.0_pReal .and. &
plastic_dislotwin_qShearBand(instance) <= 0.0_pReal) &
@ -1628,7 +1631,8 @@ end subroutine plastic_dislotwin_microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,ipc,ip,el)
use prec, only: &
tol_math_check
tol_math_check, &
dNeq
use math, only: &
math_Plain3333to99, &
math_Mandel6to33, &
@ -1775,8 +1779,8 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature
!--------------------------------------------------------------------------------------------------
! Shear banding (shearband) part
if(abs(plastic_dislotwin_sbVelocity(instance)) > tiny(0.0_pReal) .and. &
abs(plastic_dislotwin_sbResistance(instance)) > tiny(0.0_pReal)) then
if(dNeq(plastic_dislotwin_sbVelocity(instance), 0.0_pReal) .and. &
dNeq(plastic_dislotwin_sbResistance(instance),0.0_pReal)) then
gdot_sb = 0.0_pReal
dgdot_dtausb = 0.0_pReal
call math_eigenValuesVectorsSym(math_Mandel6to33(Tstar_v),eigValues,eigVectors,error)
@ -1942,7 +1946,8 @@ end subroutine plastic_dislotwin_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el)
use prec, only: &
tol_math_check
tol_math_check, &
dEq
use math, only: &
pi
use material, only: &
@ -2043,7 +2048,7 @@ subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el)
!* Dipole formation
EdgeDipMinDistance = &
plastic_dislotwin_CEdgeDipMinDistance(instance)*plastic_dislotwin_burgersPerSlipSystem(j,instance)
if (abs(tau_slip(j)) <= tiny(0.0_pReal)) then
if (dEq(tau_slip(j),0.0_pReal)) then
DotRhoDipFormation = 0.0_pReal
else
EdgeDipDistance = &
@ -2071,10 +2076,10 @@ subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el)
plastic_dislotwin_CAtomicVolume(instance)*plastic_dislotwin_burgersPerSlipSystem(j,instance)**(3.0_pReal)
VacancyDiffusion = &
plastic_dislotwin_D0(instance)*exp(-plastic_dislotwin_Qsd(instance)/(kB*Temperature))
if (abs(tau_slip(j)) <= tiny(0.0_pReal)) then
if (dEq(tau_slip(j),0.0_pReal)) then
DotRhoEdgeDipClimb = 0.0_pReal
else
if (EdgeDipDistance-EdgeDipMinDistance <= tiny(0.0_pReal)) then
if (dEq(EdgeDipDistance-EdgeDipMinDistance,0.0_pReal)) then
DotRhoEdgeDipClimb = 0.0_pReal
else
ClimbVelocity = 3.0_pReal*lattice_mu(ph)*VacancyDiffusion*AtomicVolume/ &
@ -2189,7 +2194,8 @@ end subroutine plastic_dislotwin_dotState
!--------------------------------------------------------------------------------------------------
function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
use prec, only: &
tol_math_check
tol_math_check, &
dEq
use math, only: &
pi, &
math_Mandel6to33, &
@ -2504,11 +2510,8 @@ function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
endif
!* Stress exponent
if (abs(gdot_slip(j))<=tiny(0.0_pReal)) then
plastic_dislotwin_postResults(c+j) = 0.0_pReal
else
plastic_dislotwin_postResults(c+j) = (tau/gdot_slip(j))*dgdot_dtauslip
endif
plastic_dislotwin_postResults(c+j) = &
merge(0.0_pReal,(tau/gdot_slip(j))*dgdot_dtauslip,dEq(gdot_slip(j),0.0_pReal))
enddo ; enddo
c = c + ns
case (sb_eigenvalues_ID)

61
code/plastic_isotropic.f90
View File

@ -7,7 +7,6 @@
!! untextured polycrystal
!--------------------------------------------------------------------------------------------------
module plastic_isotropic
use prec, only: &
pReal,&
pInt, &
@ -140,9 +139,10 @@ subroutine plastic_isotropic_init(fileUnit)
sizeDeltaState
character(len=65536) :: &
tag = '', &
outputtag = '', &
line = '', &
extmsg = ''
character(len=64) :: &
outputtag = ''
integer(pInt) :: NipcMyPhase
mainProcess: if (worldrank == 0) then
@ -382,8 +382,7 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
math_mul33xx33, &
math_transpose33
use material, only: &
phaseAt, phasememberAt, &
plasticState, &
phasememberAt, &
material_phase, &
phase_plasticityInstance
@ -413,7 +412,7 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
k, l, m, n
of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember
instance = phase_plasticityInstance(phaseAt(ipc,ip,el)) ! "phaseAt" equivalent to "material_phase" !!
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
Tstar_dev_33 = math_deviatoric33(math_Mandel6to33(Tstar_v)) ! deviatoric part of 2nd Piola-Kirchhoff stress
squarenorm_Tstar_dev = math_mul33xx33(Tstar_dev_33,Tstar_dev_33)
@ -463,8 +462,7 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar_3333,Tstar_v,ipc,ip,e
math_spherical33, &
math_mul33xx33
use material, only: &
phaseAt, phasememberAt, &
plasticState, &
phasememberAt, &
material_phase, &
phase_plasticityInstance
@ -491,34 +489,29 @@ real(pReal) :: &
k, l, m, n
of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember
instance = phase_plasticityInstance(phaseAt(ipc,ip,el)) ! "phaseAt" equivalent to "material_phase" !!
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
Tstar_sph_33 = math_spherical33(math_Mandel6to33(Tstar_v)) ! spherical part of 2nd Piola-Kirchhoff stress
squarenorm_Tstar_sph = math_mul33xx33(Tstar_sph_33,Tstar_sph_33)
norm_Tstar_sph = sqrt(squarenorm_Tstar_sph)
if (param(instance)%dilatation) then
if (norm_Tstar_sph <= 0.0_pReal) then ! Tstar == 0 --> both Li and dLi_dTstar are zero
Li = 0.0_pReal
dLi_dTstar_3333 = 0.0_pReal
else
gamma_dot = param(instance)%gdot0 &
* (sqrt(1.5_pReal) * norm_Tstar_sph / param(instance)%fTaylor / state(instance)%flowstress(of) ) &
**param(instance)%n
if (param(instance)%dilatation .and. norm_Tstar_sph > 0.0_pReal) then ! Tstar == 0 or J2 plascitiy --> both Li and dLi_dTstar are zero
gamma_dot = param(instance)%gdot0 &
* (sqrt(1.5_pReal) * norm_Tstar_sph / param(instance)%fTaylor / state(instance)%flowstress(of) ) &
**param(instance)%n
Li = Tstar_sph_33/norm_Tstar_sph * gamma_dot/param(instance)%fTaylor
Li = Tstar_sph_33/norm_Tstar_sph * gamma_dot/param(instance)%fTaylor
!--------------------------------------------------------------------------------------------------
! Calculation of the tangent of Li
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLi_dTstar_3333(k,l,m,n) = (param(instance)%n-1.0_pReal) * &
Tstar_sph_33(k,l)*Tstar_sph_33(m,n) / squarenorm_Tstar_sph
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
dLi_dTstar_3333(k,l,k,l) = dLi_dTstar_3333(k,l,k,l) + 1.0_pReal
!--------------------------------------------------------------------------------------------------
! Calculation of the tangent of Li
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLi_dTstar_3333(k,l,m,n) = (param(instance)%n-1.0_pReal) * &
Tstar_sph_33(k,l)*Tstar_sph_33(m,n) / squarenorm_Tstar_sph
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
dLi_dTstar_3333(k,l,k,l) = dLi_dTstar_3333(k,l,k,l) + 1.0_pReal
dLi_dTstar_3333 = gamma_dot / param(instance)%fTaylor * &
dLi_dTstar_3333 / norm_Tstar_sph
endif
dLi_dTstar_3333 = gamma_dot / param(instance)%fTaylor * &
dLi_dTstar_3333 / norm_Tstar_sph
else
Li = 0.0_pReal
dLi_dTstar_3333 = 0.0_pReal
@ -531,11 +524,12 @@ end subroutine plastic_isotropic_LiAndItsTangent
!> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
use prec, only: &
dEq
use math, only: &
math_mul6x6
use material, only: &
phaseAt, phasememberAt, &
plasticState, &
phasememberAt, &
material_phase, &
phase_plasticityInstance
@ -558,7 +552,7 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
of !< shortcut notation for offset position in state array
of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember
instance = phase_plasticityInstance(phaseAt(ipc,ip,el)) ! "phaseAt" equivalent to "material_phase" !!
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
!--------------------------------------------------------------------------------------------------
! norm of (deviatoric) 2nd Piola-Kirchhoff stress
@ -578,7 +572,7 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
!--------------------------------------------------------------------------------------------------
! hardening coefficient
if (abs(gamma_dot) > 1e-12_pReal) then
if (abs(param(instance)%tausat_SinhFitA) <= tiny(0.0_pReal)) then
if (dEq(param(instance)%tausat_SinhFitA,0.0_pReal)) then
saturation = param(instance)%tausat
else
saturation = ( param(instance)%tausat &
@ -614,8 +608,7 @@ function plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
math_mul6x6
use material, only: &
material_phase, &
plasticState, &
phaseAt, phasememberAt, &
phasememberAt, &
phase_plasticityInstance
implicit none
@ -639,7 +632,7 @@ function plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
o
of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember
instance = phase_plasticityInstance(phaseAt(ipc,ip,el)) ! "phaseAt" equivalent to "material_phase" !!
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
!--------------------------------------------------------------------------------------------------
! norm of (deviatoric) 2nd Piola-Kirchhoff stress

564
code/plastic_j2.f90
View File

@ -1,564 +0,0 @@
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for isotropic (J2) plasticity
!> @details Isotropic (J2) Plasticity which resembles the phenopowerlaw plasticity without
!! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an
!! untextured polycrystal
!--------------------------------------------------------------------------------------------------
module plastic_j2
#ifdef HDF
use hdf5, only: &
HID_T
#endif
use prec, only: &
pReal,&
pInt
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
plastic_j2_sizePostResults !< cumulative size of post results
integer(pInt), dimension(:,:), allocatable, target, public :: &
plastic_j2_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
plastic_j2_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
plastic_j2_Noutput !< number of outputs per instance
real(pReal), dimension(:), allocatable, private :: &
plastic_j2_fTaylor, & !< Taylor factor
plastic_j2_tau0, & !< initial plastic stress
plastic_j2_gdot0, & !< reference velocity
plastic_j2_n, & !< Visco-plastic parameter
!--------------------------------------------------------------------------------------------------
! h0 as function of h0 = A + B log (gammadot)
plastic_j2_h0, &
plastic_j2_h0_slopeLnRate, &
plastic_j2_tausat, & !< final plastic stress
plastic_j2_a, &
plastic_j2_aTolResistance, &
plastic_j2_aTolShear, &
!--------------------------------------------------------------------------------------------------
! tausat += (asinh((gammadot / SinhFitA)**(1 / SinhFitD)))**(1 / SinhFitC) / (SinhFitB * (gammadot / gammadot0)**(1/n))
plastic_j2_tausat_SinhFitA, & !< fitting parameter for normalized strain rate vs. stress function
plastic_j2_tausat_SinhFitB, & !< fitting parameter for normalized strain rate vs. stress function
plastic_j2_tausat_SinhFitC, & !< fitting parameter for normalized strain rate vs. stress function
plastic_j2_tausat_SinhFitD !< fitting parameter for normalized strain rate vs. stress function
enum, bind(c)
enumerator :: undefined_ID, &
flowstress_ID, &
strainrate_ID
end enum
integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: &
plastic_j2_outputID !< ID of each post result output
#ifdef HDF
type plastic_j2_tOutput
real(pReal), dimension(:), allocatable, private :: &
flowstress, &
strainrate
logical :: flowstressActive = .false., strainrateActive = .false. ! if we can write the output block wise, this is not needed anymore because we can do an if(allocated(xxx))
end type plastic_j2_tOutput
type(plastic_j2_tOutput), allocatable, dimension(:) :: plastic_j2_Output2
integer(HID_T), allocatable, dimension(:) :: outID
#endif
public :: &
plastic_j2_init, &
plastic_j2_LpAndItsTangent, &
plastic_j2_dotState, &
plastic_j2_postResults
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine plastic_j2_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
#ifdef HDF
use hdf5
#endif
use debug, only: &
debug_level, &
debug_constitutive, &
debug_levelBasic
use numerics, only: &
analyticJaco, &
worldrank, &
numerics_integrator
use math, only: &
math_Mandel3333to66, &
math_Voigt66to3333
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_error, &
IO_timeStamp, &
#ifdef HDF
tempResults, &
HDF5_addGroup, &
HDF5_addScalarDataset,&
#endif
IO_EOF
use material, only: &
phase_plasticity, &
phase_plasticityInstance, &
phase_Noutput, &
PLASTICITY_J2_label, &
PLASTICITY_J2_ID, &
material_phase, &
plasticState, &
MATERIAL_partPhase
use lattice
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: &
o, &
phase, &
maxNinstance, &
instance, &
mySize, &
sizeDotState, &
sizeState, &
sizeDeltaState
character(len=65536) :: &
tag = '', &
line = ''
integer(pInt) :: NofMyPhase
#ifdef HDF
character(len=5) :: &
str1
integer(HID_T) :: ID,ID2,ID4
#endif
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_J2_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_plasticity == PLASTICITY_J2_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
#ifdef HDF
allocate(plastic_j2_Output2(maxNinstance))
allocate(outID(maxNinstance))
#endif
allocate(plastic_j2_sizePostResults(maxNinstance), source=0_pInt)
allocate(plastic_j2_sizePostResult(maxval(phase_Noutput), maxNinstance),source=0_pInt)
allocate(plastic_j2_output(maxval(phase_Noutput), maxNinstance))
plastic_j2_output = ''
allocate(plastic_j2_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID)
allocate(plastic_j2_Noutput(maxNinstance), source=0_pInt)
allocate(plastic_j2_fTaylor(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_tau0(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_gdot0(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_n(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_h0(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_h0_slopeLnRate(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_tausat(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_a(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_aTolResistance(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_aTolShear (maxNinstance), source=0.0_pReal)
allocate(plastic_j2_tausat_SinhFitA(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_tausat_SinhFitB(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_tausat_SinhFitC(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_tausat_SinhFitD(maxNinstance), source=0.0_pReal)
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next section
phase = phase + 1_pInt ! advance section counter
if (phase_plasticity(phase) == PLASTICITY_J2_ID) then
instance = phase_plasticityInstance(phase)
endif
cycle ! skip to next line
endif
if (phase > 0_pInt ) then; if (phase_plasticity(phase) == PLASTICITY_J2_ID) then ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran
instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('(output)')
select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
case ('flowstress')
plastic_j2_Noutput(instance) = plastic_j2_Noutput(instance) + 1_pInt
plastic_j2_outputID(plastic_j2_Noutput(instance),instance) = flowstress_ID
plastic_j2_output(plastic_j2_Noutput(instance),instance) = &
IO_lc(IO_stringValue(line,chunkPos,2_pInt))
case ('strainrate')
plastic_j2_Noutput(instance) = plastic_j2_Noutput(instance) + 1_pInt
plastic_j2_outputID(plastic_j2_Noutput(instance),instance) = strainrate_ID
plastic_j2_output(plastic_j2_Noutput(instance),instance) = &
IO_lc(IO_stringValue(line,chunkPos,2_pInt))
case default
end select
case ('tau0')
plastic_j2_tau0(instance) = IO_floatValue(line,chunkPos,2_pInt)
if (plastic_j2_tau0(instance) < 0.0_pReal) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
case ('gdot0')
plastic_j2_gdot0(instance) = IO_floatValue(line,chunkPos,2_pInt)
if (plastic_j2_gdot0(instance) <= 0.0_pReal) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
case ('n')
plastic_j2_n(instance) = IO_floatValue(line,chunkPos,2_pInt)
if (plastic_j2_n(instance) <= 0.0_pReal) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
case ('h0')
plastic_j2_h0(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('h0_slope','slopelnrate')
plastic_j2_h0_slopeLnRate(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat')
plastic_j2_tausat(instance) = IO_floatValue(line,chunkPos,2_pInt)
if (plastic_j2_tausat(instance) <= 0.0_pReal) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
case ('tausat_sinhfita')
plastic_j2_tausat_SinhFitA(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfitb')
plastic_j2_tausat_SinhFitB(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfitc')
plastic_j2_tausat_SinhFitC(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfitd')
plastic_j2_tausat_SinhFitD(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('a', 'w0')
plastic_j2_a(instance) = IO_floatValue(line,chunkPos,2_pInt)
if (plastic_j2_a(instance) <= 0.0_pReal) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
case ('taylorfactor')
plastic_j2_fTaylor(instance) = IO_floatValue(line,chunkPos,2_pInt)
if (plastic_j2_fTaylor(instance) <= 0.0_pReal) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
case ('atol_resistance')
plastic_j2_aTolResistance(instance) = IO_floatValue(line,chunkPos,2_pInt)
if (plastic_j2_aTolResistance(instance) <= 0.0_pReal) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
case ('atol_shear')
plastic_j2_aTolShear(instance) = IO_floatValue(line,chunkPos,2_pInt)
case default
end select
endif; endif
enddo parsingFile
initializeInstances: do phase = 1_pInt, size(phase_plasticity)
myPhase: if (phase_plasticity(phase) == PLASTICITY_j2_ID) then
NofMyPhase=count(material_phase==phase)
instance = phase_plasticityInstance(phase)
!--------------------------------------------------------------------------------------------------
! sanity checks
if (plastic_j2_aTolShear(instance) <= 0.0_pReal) &
plastic_j2_aTolShear(instance) = 1.0e-6_pReal ! default absolute tolerance 1e-6
!--------------------------------------------------------------------------------------------------
! Determine size of postResults array
outputsLoop: do o = 1_pInt,plastic_j2_Noutput(instance)
select case(plastic_j2_outputID(o,instance))
case(flowstress_ID,strainrate_ID)
mySize = 1_pInt
case default
end select
outputFound: if (mySize > 0_pInt) then
plastic_j2_sizePostResult(o,instance) = mySize
plastic_j2_sizePostResults(instance) = &
plastic_j2_sizePostResults(instance) + mySize
endif outputFound
enddo outputsLoop
!--------------------------------------------------------------------------------------------------
! allocate state arrays
sizeState = 2_pInt
sizeDotState = sizeState
sizeDeltaState = 0_pInt
plasticState(phase)%sizeState = sizeState
plasticState(phase)%sizeDotState = sizeDotState
plasticState(phase)%sizeDeltaState = sizeDeltaState
plasticState(phase)%sizePostResults = plastic_j2_sizePostResults(instance)
plasticState(phase)%nSlip = 1
plasticState(phase)%nTwin = 0
plasticState(phase)%nTrans= 0
allocate(plasticState(phase)%aTolState ( sizeState))
plasticState(phase)%aTolState(1) = plastic_j2_aTolResistance(instance)
plasticState(phase)%aTolState(2) = plastic_j2_aTolShear(instance)
allocate(plasticState(phase)%state0 ( sizeState,NofMyPhase))
plasticState(phase)%state0(1,1:NofMyPhase) = plastic_j2_tau0(instance)
plasticState(phase)%state0(2,1:NofMyPhase) = 0.0_pReal
allocate(plasticState(phase)%partionedState0 ( sizeState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%subState0 ( sizeState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%state ( sizeState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase),source=0.0_pReal)
if (.not. analyticJaco) then
allocate(plasticState(phase)%state_backup ( sizeState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase),source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%previousDotState2(sizeDotState,NofMyPhase),source=0.0_pReal)
endif
if (any(numerics_integrator == 4_pInt)) &
allocate(plasticState(phase)%RK4dotState (sizeDotState,NofMyPhase),source=0.0_pReal)
if (any(numerics_integrator == 5_pInt)) &
allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal)
plasticState(phase)%slipRate => plasticState(phase)%dotState(2:2,1:NofMyPhase)
plasticState(phase)%accumulatedSlip => plasticState(phase)%state (2:2,1:NofMyPhase)
endif myPhase
enddo initializeInstances
end subroutine plastic_j2_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
use math, only: &
math_mul6x6, &
math_Mandel6to33, &
math_Plain3333to99, &
math_deviatoric33, &
math_mul33xx33
use material, only: &
phaseAt, phasememberAt, &
plasticState, &
material_phase, &
phase_plasticityInstance
implicit none
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(9,9), intent(out) :: &
dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(3,3) :: &
Tstar_dev_33 !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
real(pReal), dimension(3,3,3,3) :: &
dLp_dTstar_3333 !< derivative of Lp with respect to Tstar as 4th order tensor
real(pReal) :: &
gamma_dot, & !< strainrate
norm_Tstar_dev, & !< euclidean norm of Tstar_dev
squarenorm_Tstar_dev !< square of the euclidean norm of Tstar_dev
integer(pInt) :: &
instance, &
k, l, m, n
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
Tstar_dev_33 = math_deviatoric33(math_Mandel6to33(Tstar_v)) ! deviatoric part of 2nd Piola-Kirchhoff stress
squarenorm_Tstar_dev = math_mul33xx33(Tstar_dev_33,Tstar_dev_33)
norm_Tstar_dev = sqrt(squarenorm_Tstar_dev)
if (norm_Tstar_dev <= 0.0_pReal) then ! Tstar == 0 --> both Lp and dLp_dTstar are zero
Lp = 0.0_pReal
dLp_dTstar99 = 0.0_pReal
else
gamma_dot = plastic_j2_gdot0(instance) &
* (sqrt(1.5_pReal) * norm_Tstar_dev / (plastic_j2_fTaylor(instance) * &
plasticState(phaseAt(ipc,ip,el))%state(1,phasememberAt(ipc,ip,el)))) &
**plastic_j2_n(instance)
Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/plastic_j2_fTaylor(instance)
!--------------------------------------------------------------------------------------------------
! Calculation of the tangent of Lp
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar_3333(k,l,m,n) = (plastic_j2_n(instance)-1.0_pReal) * &
Tstar_dev_33(k,l)*Tstar_dev_33(m,n) / squarenorm_Tstar_dev
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
dLp_dTstar_3333(k,l,k,l) = dLp_dTstar_3333(k,l,k,l) + 1.0_pReal
forall (k=1_pInt:3_pInt,m=1_pInt:3_pInt) &
dLp_dTstar_3333(k,k,m,m) = dLp_dTstar_3333(k,k,m,m) - 1.0_pReal/3.0_pReal
dLp_dTstar99 = math_Plain3333to99(gamma_dot / plastic_j2_fTaylor(instance) * &
dLp_dTstar_3333 / norm_Tstar_dev)
end if
end subroutine plastic_j2_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_j2_dotState(Tstar_v,ipc,ip,el)
use math, only: &
math_mul6x6
use material, only: &
phaseAt, phasememberAt, &
plasticState, &
material_phase, &
phase_plasticityInstance
implicit none
real(pReal), dimension(6), intent(in):: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(6) :: &
Tstar_dev_v !< deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal) :: &
gamma_dot, & !< strainrate
hardening, & !< hardening coefficient
saturation, & !< saturation resistance
norm_Tstar_dev !< euclidean norm of Tstar_dev
integer(pInt) :: &
instance, & !< instance of my instance (unique number of my constitutive model)
of, & !< shortcut notation for offset position in state array
ph !< shortcut notation for phase ID (unique number of all phases, regardless of constitutive model)
of = phasememberAt(ipc,ip,el)
ph = phaseAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
!--------------------------------------------------------------------------------------------------
! norm of deviatoric part of 2nd Piola-Kirchhoff stress
Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal
Tstar_dev_v(4:6) = Tstar_v(4:6)
norm_Tstar_dev = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v))
!--------------------------------------------------------------------------------------------------
! strain rate
gamma_dot = plastic_j2_gdot0(instance) * ( sqrt(1.5_pReal) * norm_Tstar_dev &
/ &!-----------------------------------------------------------------------------------
(plastic_j2_fTaylor(instance)*plasticState(ph)%state(1,of)) )**plastic_j2_n(instance)
!--------------------------------------------------------------------------------------------------
! hardening coefficient
if (abs(gamma_dot) > 1e-12_pReal) then
if (abs(plastic_j2_tausat_SinhFitA(instance)) <= tiny(0.0_pReal)) then
saturation = plastic_j2_tausat(instance)
else
saturation = ( plastic_j2_tausat(instance) &
+ ( log( ( gamma_dot / plastic_j2_tausat_SinhFitA(instance)&
)**(1.0_pReal / plastic_j2_tausat_SinhFitD(instance))&
+ sqrt( ( gamma_dot / plastic_j2_tausat_SinhFitA(instance) &
)**(2.0_pReal / plastic_j2_tausat_SinhFitD(instance)) &
+ 1.0_pReal ) &
) & ! asinh(K) = ln(K + sqrt(K^2 +1))
)**(1.0_pReal / plastic_j2_tausat_SinhFitC(instance)) &
/ ( plastic_j2_tausat_SinhFitB(instance) &
* (gamma_dot / plastic_j2_gdot0(instance))**(1.0_pReal / plastic_j2_n(instance)) &
) &
)
endif
hardening = ( plastic_j2_h0(instance) + plastic_j2_h0_slopeLnRate(instance) * log(gamma_dot) ) &
* abs( 1.0_pReal - plasticState(ph)%state(1,of)/saturation )**plastic_j2_a(instance) &
* sign(1.0_pReal, 1.0_pReal - plasticState(ph)%state(1,of)/saturation)
else
hardening = 0.0_pReal
endif
plasticState(ph)%dotState(1,of) = hardening * gamma_dot
plasticState(ph)%dotState(2,of) = gamma_dot
end subroutine plastic_j2_dotState
!--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results
!--------------------------------------------------------------------------------------------------
function plastic_j2_postResults(Tstar_v,ipc,ip,el)
use math, only: &
math_mul6x6
use material, only: &
material_phase, &
plasticState, &
phaseAt, phasememberAt, &
phase_plasticityInstance
implicit none
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(plastic_j2_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
plastic_j2_postResults
real(pReal), dimension(6) :: &
Tstar_dev_v ! deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal) :: &
norm_Tstar_dev ! euclidean norm of Tstar_dev
integer(pInt) :: &
instance, & !< instance of my instance (unique number of my constitutive model)
of, & !< shortcut notation for offset position in state array
ph, & !< shortcut notation for phase ID (unique number of all phases, regardless of constitutive model)
c, &
o
of = phasememberAt(ipc,ip,el)
ph = phaseAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
!--------------------------------------------------------------------------------------------------
! calculate deviatoric part of 2nd Piola-Kirchhoff stress and its norm
Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal
Tstar_dev_v(4:6) = Tstar_v(4:6)
norm_Tstar_dev = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v))
c = 0_pInt
plastic_j2_postResults = 0.0_pReal
outputsLoop: do o = 1_pInt,plastic_j2_Noutput(instance)
select case(plastic_j2_outputID(o,instance))
case (flowstress_ID)
plastic_j2_postResults(c+1_pInt) = plasticState(ph)%state(1,of)
c = c + 1_pInt
case (strainrate_ID)
plastic_j2_postResults(c+1_pInt) = &
plastic_j2_gdot0(instance) * ( sqrt(1.5_pReal) * norm_Tstar_dev &
/ &!----------------------------------------------------------------------------------
(plastic_j2_fTaylor(instance) * plasticState(ph)%state(1,of)) ) ** plastic_j2_n(instance)
c = c + 1_pInt
end select
enddo outputsLoop
end function plastic_j2_postResults
end module plastic_j2

68
code/plastic_nonlocal.f90
View File

@ -1549,6 +1549,8 @@ end subroutine plastic_nonlocal_aTolState
!> @brief calculates quantities characterizing the microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_nonlocal_microstructure(Fe, Fp, ip, el)
use prec, only: &
dEq
use IO, only: &
IO_error
use math, only: &
@ -1792,7 +1794,7 @@ if (.not. phase_localPlasticity(ph) .and. shortRangeStressCorrection(instance))
- neighbor_rhoExcess(c,s,neighbors(2))
enddo
invConnections = math_inv33(connections)
if (all(abs(invConnections) <= tiny(0.0_pReal))) & ! check for failed in version (math_inv33 returns 0) and avoid floating point equality comparison
if (all(dEq(invConnections,0.0_pReal))) &
call IO_error(-1_pInt,ext_msg='back stress calculation: inversion error')
rhoExcessGradient(c) = math_mul3x3(m(1:3,s,c), &
math_mul33x3(invConnections,rhoExcessDifferences))
@ -2200,6 +2202,8 @@ end subroutine plastic_nonlocal_LpAndItsTangent
!> @brief (instantaneous) incremental change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_nonlocal_deltaState(Tstar_v,ip,el)
use prec, only: &
dNeq
use debug, only: debug_level, &
debug_constitutive, &
debug_levelBasic, &
@ -2322,8 +2326,8 @@ dUpper(1:ns,2) = lattice_mu(ph) * burgers(1:ns,instance) / (4.0_pReal * pi * abs
forall (c = 1_pInt:2_pInt)
where(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)) >= tiny(0.0_pReal)) &
where(dNeq(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)),0.0_pReal)) &
dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) &
+ abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
dUpper(1:ns,c))
@ -2335,7 +2339,7 @@ deltaDUpper = dUpper - dUpperOld
!*** dissociation by stress increase
deltaRhoDipole2SingleStress = 0.0_pReal
forall (c=1_pInt:2_pInt, s=1_pInt:ns, deltaDUpper(s,c) < 0.0_pReal .and. &
abs(dUpperOld(s,c) - dLower(s,c)) > tiny(0.0_pReal)) &
dNeq(dUpperOld(s,c) - dLower(s,c),0.0_pReal)) &
deltaRhoDipole2SingleStress(s,8_pInt+c) = rhoDip(s,c) * deltaDUpper(s,c) &
/ (dUpperOld(s,c) - dLower(s,c))
@ -2382,7 +2386,9 @@ end subroutine plastic_nonlocal_deltaState
subroutine plastic_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, &
timestep,subfrac, ip,el)
use prec, only: DAMASK_NaN
use prec, only: DAMASK_NaN, &
dNeq, &
dEq
use numerics, only: numerics_integrationMode, &
numerics_timeSyncing
use IO, only: IO_error
@ -2616,8 +2622,8 @@ dUpper(1:ns,1) = lattice_mu(ph) * burgers(1:ns,instance) &
dUpper(1:ns,2) = lattice_mu(ph) * burgers(1:ns,instance) &
/ (4.0_pReal * pi * abs(tau))
forall (c = 1_pInt:2_pInt)
where(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)) >= tiny(0.0_pReal)) &
where(dNeq(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)),0.0_pReal)) &
dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) &
+ abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
dUpper(1:ns,c))
@ -2760,8 +2766,7 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then
endif
if (considerEnteringFlux) then
if(numerics_timeSyncing .and. (subfrac(1_pInt,neighbor_ip,neighbor_el) /= subfrac(1_pInt,ip,el))) &
then ! for timesyncing: in case of a timestep at the interface we have to use "state0" to make sure that fluxes n both sides are equal
if(numerics_timeSyncing .and. (dNeq(subfrac(1,neighbor_ip,neighbor_el),subfrac(1,ip,el)))) then ! for timesyncing: in case of a timestep at the interface we have to use "state0" to make sure that fluxes n both sides are equal
forall (s = 1:ns, t = 1_pInt:4_pInt)
neighbor_v(s,t) = plasticState(np)%state0(iV (s,t,neighbor_instance),no)
@ -2830,11 +2835,11 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then
my_rhoSgl = rhoSgl
my_v = v
if(numerics_timeSyncing) then
if (abs(subfrac(1_pInt,ip,el))<= tiny(0.0_pReal)) then
if (dEq(subfrac(1_pInt,ip,el),0.0_pReal)) then
my_rhoSgl = rhoSgl0
my_v = v0
elseif (neighbor_n > 0_pInt) then
if (abs(subfrac(1_pInt,neighbor_ip,neighbor_el))<= tiny(0.0_pReal)) then
if (dEq(subfrac(1_pInt,neighbor_ip,neighbor_el),0.0_pReal)) then
my_rhoSgl = rhoSgl0
my_v = v0
endif
@ -3078,13 +3083,11 @@ slipDirection(1:3,1:ns) = lattice_sd(1:3, slipSystemLattice(1:ns,instance), ph)
!*** start out fully compatible
my_compatibility = 0.0_pReal
forall(s1 = 1_pInt:ns) &
my_compatibility(1:2,s1,s1,1:Nneighbors) = 1.0_pReal
forall(s1 = 1_pInt:ns) my_compatibility(1:2,s1,s1,1:Nneighbors) = 1.0_pReal
!*** Loop thrugh neighbors and check whether there is any my_compatibility.
do n = 1_pInt,Nneighbors
neighbors: do n = 1_pInt,Nneighbors
neighbor_e = mesh_ipNeighborhood(1,n,i,e)
neighbor_i = mesh_ipNeighborhood(2,n,i,e)
@ -3093,8 +3096,7 @@ do n = 1_pInt,Nneighbors
!* Set surface transmissivity to the value specified in the material.config
if (neighbor_e <= 0_pInt .or. neighbor_i <= 0_pInt) then
forall(s1 = 1_pInt:ns) &
my_compatibility(1:2,s1,s1,n) = sqrt(surfaceTransmissivity(instance))
forall(s1 = 1_pInt:ns) my_compatibility(1:2,s1,s1,n) = sqrt(surfaceTransmissivity(instance))
cycle
endif
@ -3107,10 +3109,8 @@ do n = 1_pInt,Nneighbors
neighbor_phase = material_phase(1,neighbor_i,neighbor_e)
if (neighbor_phase /= ph) then
if (.not. phase_localPlasticity(neighbor_phase) .and. .not. phase_localPlasticity(ph)) then
forall(s1 = 1_pInt:ns) &
my_compatibility(1:2,s1,s1,n) = 0.0_pReal ! = sqrt(0.0)
endif
if (.not. phase_localPlasticity(neighbor_phase) .and. .not. phase_localPlasticity(ph))&
forall(s1 = 1_pInt:ns) my_compatibility(1:2,s1,s1,n) = 0.0_pReal
cycle
endif
@ -3141,33 +3141,33 @@ do n = 1_pInt,Nneighbors
else
absoluteMisorientation = lattice_qDisorientation(orientation(1:4,1,i,e), &
orientation(1:4,1,neighbor_i,neighbor_e)) ! no symmetry
do s1 = 1_pInt,ns ! my slip systems
do s2 = 1_pInt,ns ! my neighbor's slip systems
mySlipSystems: do s1 = 1_pInt,ns
neighborSlipSystems: do s2 = 1_pInt,ns
my_compatibility(1,s2,s1,n) = math_mul3x3(slipNormal(1:3,s1), math_qRot(absoluteMisorientation, slipNormal(1:3,s2))) &
* abs(math_mul3x3(slipDirection(1:3,s1), math_qRot(absoluteMisorientation, slipDirection(1:3,s2))))
my_compatibility(2,s2,s1,n) = abs(math_mul3x3(slipNormal(1:3,s1), math_qRot(absoluteMisorientation, slipNormal(1:3,s2)))) &
* abs(math_mul3x3(slipDirection(1:3,s1), math_qRot(absoluteMisorientation, slipDirection(1:3,s2))))
enddo
enddo neighborSlipSystems
my_compatibilitySum = 0.0_pReal
belowThreshold = .true.
do while (my_compatibilitySum < 1.0_pReal .and. any(belowThreshold(1:ns)))
thresholdValue = maxval(my_compatibility(2,1:ns,s1,n), belowThreshold(1:ns)) ! screws always positive
nThresholdValues = real(count(my_compatibility(2,1:ns,s1,n) == thresholdValue),pReal)
nThresholdValues = real(count(my_compatibility(2,1:ns,s1,n) >= thresholdValue),pReal)
where (my_compatibility(2,1:ns,s1,n) >= thresholdValue) &
belowThreshold(1:ns) = .false.
if (my_compatibilitySum + thresholdValue * nThresholdValues > 1.0_pReal) &
where (abs(my_compatibility(1:2,1:ns,s1,n)) == thresholdValue) & ! MD: rather check below threshold?
where (abs(my_compatibility(1:2,1:ns,s1,n)) >= thresholdValue) & ! MD: rather check below threshold?
my_compatibility(1:2,1:ns,s1,n) = sign((1.0_pReal - my_compatibilitySum) &
/ nThresholdValues, my_compatibility(1:2,1:ns,s1,n))
my_compatibilitySum = my_compatibilitySum + nThresholdValues * thresholdValue
enddo
where (belowThreshold(1:ns)) my_compatibility(1,1:ns,s1,n) = 0.0_pReal
where (belowThreshold(1:ns)) my_compatibility(2,1:ns,s1,n) = 0.0_pReal
enddo ! my slip systems cycle
enddo mySlipSystems
endif
enddo ! neighbor cycle
enddo neighbors
compatibility(1:2,1:ns,1:ns,1:Nneighbors,i,e) = my_compatibility
@ -3177,6 +3177,8 @@ end subroutine plastic_nonlocal_updateCompatibility
!* calculates quantities characterizing the microstructure *
!*********************************************************************
function plastic_nonlocal_dislocationstress(Fe, ip, el)
use prec, only: &
dEq
use math, only: math_mul33x33, &
math_mul33x3, &
math_inv33, &
@ -3389,7 +3391,7 @@ if (.not. phase_localPlasticity(ph)) then
Rsquare = R * R
Rcube = Rsquare * R
denominator = R * (R + flipSign * lambda)
if (abs(denominator)<= tiny(0.0_pReal)) exit ipLoop
if (dEq(denominator,0.0_pReal)) exit ipLoop
sigma(1,1) = sigma(1,1) - real(side,pReal) &
* flipSign * z / denominator &
@ -3434,7 +3436,7 @@ if (.not. phase_localPlasticity(ph)) then
Rsquare = R * R
Rcube = Rsquare * R
denominator = R * (R + flipSign * lambda)
if (abs(denominator)<= tiny(0.0_pReal)) exit ipLoop
if (dEq(denominator,0.0_pReal)) exit ipLoop
sigma(1,2) = sigma(1,2) - real(side,pReal) * flipSign * z &
* (1.0_pReal - lattice_nu(ph)) / denominator &
@ -3523,6 +3525,8 @@ end function plastic_nonlocal_dislocationstress
!> @brief return array of constitutive results
!--------------------------------------------------------------------------------------------------
function plastic_nonlocal_postResults(Tstar_v,Fe,ip,el)
use prec, only: &
dNeq
use math, only: &
math_mul6x6, &
math_mul33x3, &
@ -3639,8 +3643,8 @@ dUpper(1:ns,1) = lattice_mu(ph) * burgers(1:ns,instance) &
dUpper(1:ns,2) = lattice_mu(ph) * burgers(1:ns,instance) &
/ (4.0_pReal * pi * abs(tau))
forall (c = 1_pInt:2_pInt)
where(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)) >= tiny(0.0_pReal)) &
where(dNeq(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)),0.0_pReal)) &
dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) &
+ abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
dUpper(1:ns,c))

14
code/plastic_phenoplus.f90
View File

@ -112,6 +112,8 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine plastic_phenoplus_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: &
dEq
use debug, only: &
debug_level, &
debug_constitutive,&
@ -481,7 +483,7 @@ subroutine plastic_phenoplus_init(fileUnit)
if (any(plastic_phenoplus_tausat_slip(:,instance) <= 0.0_pReal .and. &
plastic_phenoplus_Nslip(:,instance) > 0)) &
call IO_error(211_pInt,el=instance,ext_msg='tausat_slip ('//PLASTICITY_PHENOPLUS_label//')')
if (any(abs(plastic_phenoplus_a_slip(instance)) <= tiny(0.0_pReal) .and. &
if (any(dEq(plastic_phenoplus_a_slip(instance),0.0_pReal) .and. &
plastic_phenoplus_Nslip(:,instance) > 0)) &
call IO_error(211_pInt,el=instance,ext_msg='a_slip ('//PLASTICITY_PHENOPLUS_label//')')
if (any(plastic_phenoplus_tau0_twin(:,instance) < 0.0_pReal .and. &
@ -873,7 +875,7 @@ subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el)
ENDDO LOOPFINDNEISHEAR
!***calculate the average accumulative shear and use it as cutoff
avg_acshear_ne = SUM(ne_acshear)/ns
avg_acshear_ne = sum(ne_acshear)/real(ns,pReal)
!***
IF (ph==neighbor_ph) THEN
@ -923,6 +925,8 @@ end subroutine plastic_phenoplus_microstructure
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
use prec, only: &
dNeq
use math, only: &
math_Plain3333to99, &
math_Mandel6to33
@ -1038,7 +1042,7 @@ subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
(gdot_slip_pos+gdot_slip_neg)*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph)
! Calculation of the tangent of Lp
if (abs(gdot_slip_pos) > tiny(0.0_pReal)) then
if (dNeq(gdot_slip_pos,0.0_pReal)) then
dgdot_dtauslip_pos = gdot_slip_pos*plastic_phenoplus_n_slip(instance)/tau_slip_pos
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
@ -1046,7 +1050,7 @@ subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
nonSchmid_tensor(m,n,1)
endif
if (abs(gdot_slip_neg) > tiny(0.0_pReal)) then
if (dNeq(gdot_slip_neg,0.0_pReal)) then
dgdot_dtauslip_neg = gdot_slip_neg*plastic_phenoplus_n_slip(instance)/tau_slip_neg
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
@ -1073,7 +1077,7 @@ subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
Lp = Lp + gdot_twin*lattice_Stwin(1:3,1:3,index_myFamily+i,ph)
! Calculation of the tangent of Lp
if (abs(gdot_twin) > tiny(0.0_pReal)) then
if (dNeq(gdot_twin,0.0_pReal)) then
dgdot_dtautwin = gdot_twin*plastic_phenoplus_n_twin(instance)/tau_twin
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &

80
code/plastic_phenopowerlaw.f90
View File

@ -1,4 +1,3 @@
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for phenomenological crystal plasticity formulation using a powerlaw
@ -60,6 +59,7 @@ module plastic_phenopowerlaw
plastic_phenopowerlaw_tau0_slip, & !< initial critical shear stress for slip (input parameter, per family)
plastic_phenopowerlaw_tau0_twin, & !< initial critical shear stress for twin (input parameter, per family)
plastic_phenopowerlaw_tausat_slip, & !< maximum critical shear stress for slip (input parameter, per family)
plastic_phenopowerlaw_H_int, & !< per family hardening activity(input parameter(optional), per family)
plastic_phenopowerlaw_nonSchmidCoeff, &
plastic_phenopowerlaw_interaction_SlipSlip, & !< interaction factors slip - slip (input parameter)
@ -124,6 +124,8 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine plastic_phenopowerlaw_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: &
dEq
use debug, only: &
debug_level, &
debug_constitutive,&
@ -197,29 +199,30 @@ subroutine plastic_phenopowerlaw_init(fileUnit)
source=0_pInt)
allocate(plastic_phenopowerlaw_output(maxval(phase_Noutput),maxNinstance))
plastic_phenopowerlaw_output = ''
allocate(plastic_phenopowerlaw_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID)
allocate(plastic_phenopowerlaw_Noutput(maxNinstance), source=0_pInt)
allocate(plastic_phenopowerlaw_Nslip(lattice_maxNslipFamily,maxNinstance), source=0_pInt)
allocate(plastic_phenopowerlaw_Ntwin(lattice_maxNtwinFamily,maxNinstance), source=0_pInt)
allocate(plastic_phenopowerlaw_Ntrans(lattice_maxNtransFamily,maxNinstance), source=0_pInt)
allocate(plastic_phenopowerlaw_totalNslip(maxNinstance), source=0_pInt)
allocate(plastic_phenopowerlaw_totalNtwin(maxNinstance), source=0_pInt)
allocate(plastic_phenopowerlaw_totalNtrans(maxNinstance), source=0_pInt)
allocate(plastic_phenopowerlaw_gdot0_slip(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_n_slip(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_tau0_slip(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_outputID(maxval(phase_Noutput),maxNinstance),source=undefined_ID)
allocate(plastic_phenopowerlaw_Noutput(maxNinstance), source=0_pInt)
allocate(plastic_phenopowerlaw_Nslip(lattice_maxNslipFamily,maxNinstance), source=0_pInt)
allocate(plastic_phenopowerlaw_Ntwin(lattice_maxNtwinFamily,maxNinstance), source=0_pInt)
allocate(plastic_phenopowerlaw_Ntrans(lattice_maxNtransFamily,maxNinstance),source=0_pInt)
allocate(plastic_phenopowerlaw_totalNslip(maxNinstance), source=0_pInt)
allocate(plastic_phenopowerlaw_totalNtwin(maxNinstance), source=0_pInt)
allocate(plastic_phenopowerlaw_totalNtrans(maxNinstance), source=0_pInt)
allocate(plastic_phenopowerlaw_gdot0_slip(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_n_slip(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_tau0_slip(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal)
allocate(plastic_phenopowerlaw_tausat_slip(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal)
allocate(plastic_phenopowerlaw_gdot0_twin(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_n_twin(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_tau0_twin(lattice_maxNtwinFamily,maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_spr(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_twinB(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_twinC(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_twinD(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_twinE(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_h0_SlipSlip(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_h0_TwinSlip(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_h0_TwinTwin(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_H_int(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal)
allocate(plastic_phenopowerlaw_gdot0_twin(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_n_twin(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_tau0_twin(lattice_maxNtwinFamily,maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_spr(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_twinB(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_twinC(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_twinD(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_twinE(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_h0_SlipSlip(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_h0_TwinSlip(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_h0_TwinTwin(maxNinstance), source=0.0_pReal)
allocate(plastic_phenopowerlaw_interaction_SlipSlip(lattice_maxNinteraction,maxNinstance), &
source=0.0_pReal)
allocate(plastic_phenopowerlaw_interaction_SlipTwin(lattice_maxNinteraction,maxNinstance), &
@ -340,7 +343,7 @@ subroutine plastic_phenopowerlaw_init(fileUnit)
do j = 1_pInt, Nchunks_SlipFamilies
plastic_phenopowerlaw_Nslip(j,instance) = IO_intValue(line,chunkPos,1_pInt+j)
enddo
case ('tausat_slip','tau0_slip')
case ('tausat_slip','tau0_slip','H_int')
tempPerSlip = 0.0_pReal
do j = 1_pInt, Nchunks_SlipFamilies
if (plastic_phenopowerlaw_Nslip(j,instance) > 0_pInt) &
@ -351,6 +354,8 @@ subroutine plastic_phenopowerlaw_init(fileUnit)
plastic_phenopowerlaw_tausat_slip(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies)
case ('tau0_slip')
plastic_phenopowerlaw_tau0_slip(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies)
case ('H_int')
plastic_phenopowerlaw_H_int(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies)
end select
!--------------------------------------------------------------------------------------------------
! parameters depending on number of twin families
@ -372,7 +377,7 @@ subroutine plastic_phenopowerlaw_init(fileUnit)
! parameters depending on number of transformation families
case ('ntrans')
if (chunkPos(1) < Nchunks_TransFamilies + 1_pInt) &
call IO_warning(51_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPOWERLAW_label//')')
call IO_warning(53_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPOWERLAW_label//')')
if (chunkPos(1) > Nchunks_TransFamilies + 1_pInt) &
call IO_error(150_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPOWERLAW_label//')')
Nchunks_TransFamilies = chunkPos(1) - 1_pInt
@ -484,7 +489,7 @@ subroutine plastic_phenopowerlaw_init(fileUnit)
if (any(plastic_phenopowerlaw_tausat_slip(:,instance) <= 0.0_pReal .and. &
plastic_phenopowerlaw_Nslip(:,instance) > 0)) &
call IO_error(211_pInt,el=instance,ext_msg='tausat_slip ('//PLASTICITY_PHENOPOWERLAW_label//')')
if (any(abs(plastic_phenopowerlaw_a_slip(instance)) <= tiny(0.0_pReal) .and. &
if (any(dEq(plastic_phenopowerlaw_a_slip(instance),0.0_pReal) .and. &
plastic_phenopowerlaw_Nslip(:,instance) > 0)) &
call IO_error(211_pInt,el=instance,ext_msg='a_slip ('//PLASTICITY_PHENOPOWERLAW_label//')')
if (any(plastic_phenopowerlaw_tau0_twin(:,instance) < 0.0_pReal .and. &
@ -683,9 +688,9 @@ subroutine plastic_phenopowerlaw_init(fileUnit)
startIndex = endIndex + 1_pInt
endIndex = endIndex +plastic_phenopowerlaw_totalNtwin(instance)
state (instance)%accshear_slip=>plasticState(phase)%state (startIndex:endIndex,:)
state0 (instance)%accshear_slip=>plasticState(phase)%state0 (startIndex:endIndex,:)
dotState(instance)%accshear_slip=>plasticState(phase)%dotState(startIndex:endIndex,:)
state (instance)%accshear_twin=>plasticState(phase)%state (startIndex:endIndex,:)
state0 (instance)%accshear_twin=>plasticState(phase)%state0 (startIndex:endIndex,:)
dotState(instance)%accshear_twin=>plasticState(phase)%dotState(startIndex:endIndex,:)
call plastic_phenopowerlaw_stateInit(phase,instance)
@ -771,6 +776,8 @@ end subroutine plastic_phenopowerlaw_aTolState
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
use prec, only: &
dNeq
use math, only: &
math_Plain3333to99, &
math_Mandel6to33
@ -860,7 +867,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,
(gdot_slip_pos+gdot_slip_neg)*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph)
! Calculation of the tangent of Lp
if (abs(gdot_slip_pos) > tiny(0.0_pReal)) then
if (dNeq(gdot_slip_pos,0.0_pReal)) then
dgdot_dtauslip_pos = gdot_slip_pos*plastic_phenopowerlaw_n_slip(instance)/tau_slip_pos
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
@ -868,7 +875,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,
nonSchmid_tensor(m,n,1)
endif
if (abs(gdot_slip_neg) > tiny(0.0_pReal)) then
if (dNeq(gdot_slip_neg,0.0_pReal)) then
dgdot_dtauslip_neg = gdot_slip_neg*plastic_phenopowerlaw_n_slip(instance)/tau_slip_neg
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
@ -895,7 +902,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,
Lp = Lp + gdot_twin*lattice_Stwin(1:3,1:3,index_myFamily+i,ph)
! Calculation of the tangent of Lp
if (abs(gdot_twin) > tiny(0.0_pReal)) then
if (dNeq(gdot_twin,0.0_pReal)) then
dgdot_dtautwin = gdot_twin*plastic_phenopowerlaw_n_twin(instance)/tau_twin
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
@ -967,7 +974,6 @@ subroutine plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
offset_accshear_twin = nSlip + nTwin + 2_pInt + nSlip
plasticState(ph)%dotState(:,of) = 0.0_pReal
!--------------------------------------------------------------------------------------------------
! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices
c_SlipSlip = plastic_phenopowerlaw_h0_SlipSlip(instance)*&
@ -986,7 +992,7 @@ subroutine plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family
slipSystems1: do i = 1_pInt,plastic_phenopowerlaw_Nslip(f,instance)
j = j+1_pInt
left_SlipSlip(j) = 1.0_pReal ! no system-dependent left part
left_SlipSlip(j) = 1.0_pReal + plastic_phenopowerlaw_H_int(f,instance) ! modified no system-dependent left part
left_SlipTwin(j) = 1.0_pReal ! no system-dependent left part
right_SlipSlip(j) = abs(1.0_pReal-plasticState(ph)%state(j,of) / &
(plastic_phenopowerlaw_tausat_slip(f,instance)+ssat_offset)) &
@ -1007,12 +1013,14 @@ subroutine plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
enddo nonSchmidSystems
gdot_slip(j) = plastic_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* &
((abs(tau_slip_pos)/(plasticState(ph)%state(j,of)))**plastic_phenopowerlaw_n_slip(instance) &
+(abs(tau_slip_neg)/(plasticState(ph)%state(j,of)))**plastic_phenopowerlaw_n_slip(instance))&
*sign(1.0_pReal,tau_slip_pos)
*sign(1.0_pReal,tau_slip_pos) &
+(abs(tau_slip_neg)/(plasticState(ph)%state(j,of)))**plastic_phenopowerlaw_n_slip(instance) &
*sign(1.0_pReal,tau_slip_neg))
enddo slipSystems1
enddo slipFamilies1
j = 0_pInt
twinFamilies1: do f = 1_pInt,lattice_maxNtwinFamily
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family

17
code/porosity_phasefield.f90
View File

@ -210,8 +210,7 @@ function porosity_phasefield_getFormationEnergy(ip,el)
enddo
porosity_phasefield_getFormationEnergy = &
porosity_phasefield_getFormationEnergy/ &
homogenization_Ngrains(mesh_element(3,el))
porosity_phasefield_getFormationEnergy/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
end function porosity_phasefield_getFormationEnergy
@ -243,8 +242,7 @@ function porosity_phasefield_getSurfaceEnergy(ip,el)
enddo
porosity_phasefield_getSurfaceEnergy = &
porosity_phasefield_getSurfaceEnergy/ &
homogenization_Ngrains(mesh_element(3,el))
porosity_phasefield_getSurfaceEnergy/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
end function porosity_phasefield_getSurfaceEnergy
@ -308,7 +306,7 @@ subroutine porosity_phasefield_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi,
enddo
W_e = W_e + sum(abs(strain*math_mul66x6(C,strain)))
enddo
W_e = W_e/homogenization_Ngrains(homog)
W_e = W_e/real(homogenization_Ngrains(homog),pReal)
phiDot = 2.0_pReal*(1.0_pReal - phi)*(1.0_pReal - Cv)*(1.0_pReal - Cv) - &
2.0_pReal*phi*(W_e + Cv*porosity_phasefield_getFormationEnergy(ip,el))/ &
@ -350,8 +348,7 @@ function porosity_phasefield_getDiffusion33(ip,el)
enddo
porosity_phasefield_getDiffusion33 = &
porosity_phasefield_getDiffusion33/ &
homogenization_Ngrains(homog)
porosity_phasefield_getDiffusion33/real(homogenization_Ngrains(homog),pReal)
end function porosity_phasefield_getDiffusion33
@ -377,10 +374,12 @@ real(pReal) function porosity_phasefield_getMobility(ip,el)
porosity_phasefield_getMobility = 0.0_pReal
do ipc = 1, homogenization_Ngrains(mesh_element(3,el))
porosity_phasefield_getMobility = porosity_phasefield_getMobility + lattice_PorosityMobility(material_phase(ipc,ip,el))
porosity_phasefield_getMobility = porosity_phasefield_getMobility &
+ lattice_PorosityMobility(material_phase(ipc,ip,el))
enddo
porosity_phasefield_getMobility = porosity_phasefield_getMobility/homogenization_Ngrains(mesh_element(3,el))
porosity_phasefield_getMobility = &
porosity_phasefield_getMobility/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
end function porosity_phasefield_getMobility

69
code/prec.f90
View File

@ -21,10 +21,10 @@ module prec
#if (FLOAT==8)
integer, parameter, public :: pReal = 8 !< floating point double precision (was selected_real_kind(15,300), number with 15 significant digits, up to 1e+-300)
#ifdef __INTEL_COMPILER
real(pReal), parameter, public :: DAMASK_NaN = Z'7FF8000000000000' !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran)
real(pReal), parameter, public :: DAMASK_NaN = Z'7FF8000000000000' !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html)
#endif
#ifdef __GFORTRAN__
real(pReal), parameter, public :: DAMASK_NaN = real(Z'7FF8000000000000',pReal) !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran)
real(pReal), parameter, public :: DAMASK_NaN = real(Z'7FF8000000000000',pReal) !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html)
#endif
#else
NO SUITABLE PRECISION FOR REAL SELECTED, STOPPING COMPILATION
@ -115,7 +115,9 @@ module prec
prec_init, &
prec_isNaN, &
dEq, &
dNeq
cEq, &
dNeq, &
cNeq
contains
@ -180,28 +182,69 @@ end function prec_isNaN
!--------------------------------------------------------------------------------------------------
!> @brief equality comparison for double precision
!> @brief equality comparison for float with double precision
! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq
! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
!--------------------------------------------------------------------------------------------------
logical elemental pure function dEq(a,b,tol)
real(pReal), intent(in) :: a,b
real(pReal), intent(in), optional :: tol
real(pReal), parameter :: eps = 2.2204460492503131E-16 ! DBL_EPSILON in C
dEq = merge(.True., .False.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b])))
implicit none
real(pReal), intent(in) :: a,b
real(pReal), intent(in), optional :: tol
real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C
dEq = merge(.True., .False.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b])))
end function dEq
!--------------------------------------------------------------------------------------------------
!> @brief inequality comparison for double precision
!> @brief inequality comparison for float with double precision
! replaces "!=" but for certain (relative) tolerance. Counterpart to dEq
! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
!--------------------------------------------------------------------------------------------------
logical elemental pure function dNeq(a,b,tol)
real(pReal), intent(in) :: a,b
real(pReal), intent(in), optional :: tol
real(pReal), parameter :: eps = 2.2204460492503131E-16 ! DBL_EPSILON in C
dNeq = merge(.False., .True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b])))
implicit none
real(pReal), intent(in) :: a,b
real(pReal), intent(in), optional :: tol
real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C
dNeq = merge(.False., .True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b])))
end function dNeq
!--------------------------------------------------------------------------------------------------
!> @brief equality comparison for complex with double precision
! replaces "==" but for certain (relative) tolerance. Counterpart to cNeq
! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
! probably a component wise comparison would be more accurate than the comparsion of the absolute
! value
!--------------------------------------------------------------------------------------------------
logical elemental pure function cEq(a,b,tol)
implicit none
complex(pReal), intent(in) :: a,b
real(pReal), intent(in), optional :: tol
real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C
cEq = merge(.True., .False.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b])))
end function cEq
!--------------------------------------------------------------------------------------------------
!> @brief inequality comparison for complex with double precision
! replaces "!=" but for certain (relative) tolerance. Counterpart to cEq
! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
! probably a component wise comparison would be more accurate than the comparsion of the absolute
! value
!--------------------------------------------------------------------------------------------------
logical elemental pure function cNeq(a,b,tol)
implicit none
complex(pReal), intent(in) :: a,b
real(pReal), intent(in), optional :: tol
real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C
cNeq = merge(.False., .True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b])))
end function cNeq
end module prec

4
code/spectral_damage.f90
View File

@ -176,10 +176,6 @@ type(tSolutionState) function spectral_damage_solution(guess,timeinc,timeinc_old
itmax, &
err_damage_tolAbs, &
err_damage_tolRel
use spectral_utilities, only: &
tBoundaryCondition, &
Utilities_maskedCompliance, &
Utilities_updateGamma
use mesh, only: &
grid, &
grid3

88
code/spectral_interface.f90
View File

@ -11,7 +11,6 @@
module DAMASK_interface
use prec, only: &
pInt
implicit none
private
#ifdef PETSc
@ -39,7 +38,6 @@ module DAMASK_interface
IIO_intValue, &
IIO_lc, &
IIO_stringPos
contains
!--------------------------------------------------------------------------------------------------
@ -63,6 +61,7 @@ subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn)
tag
integer :: &
i, &
threadLevel, &
worldrank = 0
integer, allocatable, dimension(:) :: &
chunkPos
@ -75,15 +74,22 @@ subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn)
quit,&
MPI_Comm_rank,&
PETScInitialize, &
MPI_Init_Thread, &
MPI_abort
!--------------------------------------------------------------------------------------------------
! PETSc Init
#ifdef PETSc
call PetscInitialize(PETSC_NULL_CHARACTER,ierr) ! according to PETSc manual, that should be the first line in the code
#ifdef _OPENMP
call MPI_Init_Thread(MPI_THREAD_FUNNELED,threadLevel,ierr);CHKERRQ(ierr) ! in case of OpenMP, don't rely on PETScInitialize doing MPI init
if (threadLevel<MPI_THREAD_FUNNELED) then
write(6,'(a)') 'MPI library does not support OpenMP'
call quit(1_pInt)
endif
#endif
call PetscInitialize(PETSC_NULL_CHARACTER,ierr) ! according to PETSc manual, that should be the first line in the code
CHKERRQ(ierr) ! this is a macro definition, it is case sensitive
open(6, encoding='UTF-8') ! modern fortran compilers (gfortran >4.4, ifort >11 support it)
open(6, encoding='UTF-8')
call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,ierr);CHKERRQ(ierr)
#endif
mainProcess: if (worldrank == 0) then
@ -99,7 +105,6 @@ subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn)
write(6,'(/,a)') ' <<<+- DAMASK_interface init -+>>>'
#include "compilation_info.f90"
endif mainProcess
if ( present(loadcaseParameterIn) .and. present(geometryParameterIn)) then ! both mandatory parameters given in function call
geometryArg = geometryParameterIn
loadcaseArg = loadcaseParameterIn
@ -188,7 +193,7 @@ subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn)
call quit(1_pInt)
endif
workingDirectory = storeWorkingDirectory(trim(workingDirArg),trim(geometryArg))
workingDirectory = trim(storeWorkingDirectory(trim(workingDirArg),trim(geometryArg)))
geometryFile = getGeometryFile(geometryArg)
loadCaseFile = getLoadCaseFile(loadCaseArg)
@ -221,49 +226,45 @@ end subroutine DAMASK_interface_init
!> @todo change working directory with call chdir(storeWorkingDirectory)?
!--------------------------------------------------------------------------------------------------
character(len=1024) function storeWorkingDirectory(workingDirectoryArg,geometryArg)
#ifdef __INTEL_COMPILER
use IFPORT
#endif
use system_routines, only: &
isDirectory, &
getCWD
implicit none
character(len=*), intent(in) :: workingDirectoryArg !< working directory argument
character(len=*), intent(in) :: geometryArg !< geometry argument
character(len=1024) :: cwd
character :: pathSep
logical :: dirExists
logical :: error
external :: quit
integer :: error
pathSep = getPathSep()
if (len(workingDirectoryArg)>0) then ! got working directory as input
if (workingDirectoryArg(1:1) == pathSep) then ! absolute path given as command line argument
wdGiven: if (len(workingDirectoryArg)>0) then
absolutePath: if (workingDirectoryArg(1:1) == pathSep) then
storeWorkingDirectory = workingDirectoryArg
else
error = getcwd(cwd) ! relative path given as command line argument
else absolutePath
error = getCWD(cwd)
if (error) call quit(1_pInt)
storeWorkingDirectory = trim(cwd)//pathSep//workingDirectoryArg
endif
if (storeWorkingDirectory(len(trim(storeWorkingDirectory)):len(trim(storeWorkingDirectory))) & ! if path seperator is not given, append it
/= pathSep) storeWorkingDirectory = trim(storeWorkingDirectory)//pathSep
#ifdef __INTEL_COMPILER
inquire(directory = trim(storeWorkingDirectory)//'.', exist=dirExists)
#else
inquire(file = trim(storeWorkingDirectory), exist=dirExists)
#endif
if(.not. dirExists) then ! check if the directory exists
write(6,'(a20,a,a16)') ' working directory "',trim(storeWorkingDirectory),'" does not exist'
call quit(1_pInt)
endif
else ! using path to geometry file as working dir
endif absolutePath
if (storeWorkingDirectory(len(trim(storeWorkingDirectory)):len(trim(storeWorkingDirectory))) /= pathSep) &
storeWorkingDirectory = trim(storeWorkingDirectory)//pathSep ! if path seperator is not given, append it
else wdGiven
if (geometryArg(1:1) == pathSep) then ! absolute path given as command line argument
storeWorkingDirectory = geometryArg(1:scan(geometryArg,pathSep,back=.true.))
else
error = getcwd(cwd) ! relative path given as command line argument
storeWorkingDirectory = trim(cwd)//pathSep//&
geometryArg(1:scan(geometryArg,pathSep,back=.true.))
error = getCWD(cwd) ! relative path given as command line argument
if (error) call quit(1_pInt)
storeWorkingDirectory = trim(cwd)//pathSep//geometryArg(1:scan(geometryArg,pathSep,back=.true.))
endif
endif wdGiven
storeWorkingDirectory = trim(rectifyPath(storeWorkingDirectory))
if(.not. isDirectory(trim(storeWorkingDirectory))) then ! check if the directory exists
write(6,'(a20,a,a16)') ' working directory "',trim(storeWorkingDirectory),'" does not exist'
call quit(1_pInt)
endif
storeWorkingDirectory = rectifyPath(storeWorkingDirectory)
end function storeWorkingDirectory
@ -309,9 +310,8 @@ end function getSolverJobName
!> @brief basename of geometry file with extension from command line arguments
!--------------------------------------------------------------------------------------------------
character(len=1024) function getGeometryFile(geometryParameter)
#ifdef __INTEL_COMPILER
use IFPORT
#endif
use system_routines, only: &
getCWD
implicit none
character(len=1024), intent(in) :: &
@ -319,8 +319,9 @@ character(len=1024) function getGeometryFile(geometryParameter)
character(len=1024) :: &
cwd
integer :: posExt, posSep
logical :: error
character :: pathSep
integer :: error
external :: quit
getGeometryFile = geometryParameter
pathSep = getPathSep()
@ -330,6 +331,7 @@ character(len=1024) function getGeometryFile(geometryParameter)
if (posExt <= posSep) getGeometryFile = trim(getGeometryFile)//('.geom') ! no extension present
if (scan(getGeometryFile,pathSep) /= 1) then ! relative path given as command line argument
error = getcwd(cwd)
if (error) call quit(1_pInt)
getGeometryFile = rectifyPath(trim(cwd)//pathSep//getGeometryFile)
else
getGeometryFile = rectifyPath(getGeometryFile)
@ -344,17 +346,18 @@ end function getGeometryFile
!> @brief relative path of loadcase from command line arguments
!--------------------------------------------------------------------------------------------------
character(len=1024) function getLoadCaseFile(loadCaseParameter)
#ifdef __INTEL_COMPILER
use IFPORT
#endif
use system_routines, only: &
getCWD
implicit none
character(len=1024), intent(in) :: &
loadCaseParameter
character(len=1024) :: &
cwd
integer :: posExt, posSep, error
integer :: posExt, posSep
logical :: error
character :: pathSep
external :: quit
getLoadCaseFile = loadcaseParameter
pathSep = getPathSep()
@ -364,6 +367,7 @@ character(len=1024) function getLoadCaseFile(loadCaseParameter)
if (posExt <= posSep) getLoadCaseFile = trim(getLoadCaseFile)//('.load') ! no extension present
if (scan(getLoadCaseFile,pathSep) /= 1) then ! relative path given as command line argument
error = getcwd(cwd)
if (error) call quit(1_pInt)
getLoadCaseFile = rectifyPath(trim(cwd)//pathSep//getLoadCaseFile)
else
getLoadCaseFile = rectifyPath(getLoadCaseFile)

4
code/spectral_thermal.f90
View File

@ -178,10 +178,6 @@ type(tSolutionState) function spectral_thermal_solution(guess,timeinc,timeinc_ol
itmax, &
err_thermal_tolAbs, &
err_thermal_tolRel
use spectral_utilities, only: &
tBoundaryCondition, &
Utilities_maskedCompliance, &
Utilities_updateGamma
use mesh, only: &
grid, &
grid3

109
code/spectral_utilities.f90
View File

@ -102,8 +102,6 @@ module spectral_utilities
real(pReal) :: density
end type tSolutionParams
type(tSolutionParams), private :: params
type, public :: phaseFieldDataBin !< set of parameters defining a phase field
real(pReal) :: diffusion = 0.0_pReal, & !< thermal conductivity
mobility = 0.0_pReal, & !< thermal mobility
@ -265,8 +263,9 @@ subroutine utilities_init()
enddo
elseif (divergence_correction == 2_pInt) then
do j = 1_pInt, 3_pInt
if (j /= minloc(geomSize/grid,1) .and. j /= maxloc(geomSize/grid,1)) &
scaledGeomSize = geomSize/geomSize(j)*grid(j)
if ( j /= int(minloc(geomSize/real(grid,pReal),1),pInt) &
.and. j /= int(maxloc(geomSize/real(grid,pReal),1),pInt)) &
scaledGeomSize = geomSize/geomSize(j)*real(grid(j),pReal)
enddo
else
scaledGeomSize = geomSize
@ -403,7 +402,7 @@ subroutine utilities_updateGamma(C,saveReference)
integer(pInt) :: &
i, j, k, &
l, m, n, o
logical :: ierr
logical :: err
C_ref = C
if (saveReference) then
@ -427,7 +426,7 @@ subroutine utilities_updateGamma(C,saveReference)
matA(1:3,1:3) = real(temp33_complex); matA(4:6,4:6) = real(temp33_complex)
matA(1:3,4:6) = aimag(temp33_complex); matA(4:6,1:3) = -aimag(temp33_complex)
if (abs(math_det33(matA(1:3,1:3))) > 1e-16) then
call math_invert(6_pInt, matA, matInvA, ierr)
call math_invert(6_pInt, matA, matInvA, err)
temp33_complex = cmplx(matInvA(1:3,1:3),matInvA(1:3,4:6),pReal)
forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt) &
gamma_hat(l,m,n,o,i,j,k-grid3Offset) = temp33_complex(l,n)* &
@ -543,7 +542,7 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
integer(pInt) :: &
i, j, k, &
l, m, n, o
logical :: ierr
logical :: err
if (worldrank == 0_pInt) then
@ -563,7 +562,7 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
matA(1:3,1:3) = real(temp33_complex); matA(4:6,4:6) = real(temp33_complex)
matA(1:3,4:6) = aimag(temp33_complex); matA(4:6,1:3) = -aimag(temp33_complex)
if (abs(math_det33(matA(1:3,1:3))) > 1e-16) then
call math_invert(6_pInt, matA, matInvA, ierr)
call math_invert(6_pInt, matA, matInvA, err)
temp33_complex = cmplx(matInvA(1:3,1:3),matInvA(1:3,4:6),pReal)
forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt) &
gamma_hat(l,m,n,o,1,1,1) = temp33_complex(l,n)*conjg(-xi1st(o,i,j,k))*xi1st(m,i,j,k)
@ -623,6 +622,8 @@ end subroutine utilities_fourierGreenConvolution
!> @brief calculate root mean square of divergence of field_fourier
!--------------------------------------------------------------------------------------------------
real(pReal) function utilities_divergenceRMS()
use IO, only: &
IO_error
use numerics, only: &
worldrank
use mesh, only: &
@ -631,8 +632,8 @@ real(pReal) function utilities_divergenceRMS()
grid3
implicit none
integer(pInt) :: i, j, k
PetscErrorCode :: ierr
integer(pInt) :: i, j, k, ierr
complex(pReal), dimension(3) :: rescaledGeom
external :: &
MPI_Allreduce
@ -641,6 +642,7 @@ real(pReal) function utilities_divergenceRMS()
write(6,'(/,a)') ' ... calculating divergence ................................................'
flush(6)
endif
rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal)
!--------------------------------------------------------------------------------------------------
! calculating RMS divergence criterion in Fourier space
@ -648,23 +650,24 @@ real(pReal) function utilities_divergenceRMS()
do k = 1_pInt, grid3; do j = 1_pInt, grid(2)
do i = 2_pInt, grid1Red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice.
utilities_divergenceRMS = utilities_divergenceRMS &
+ 2.0_pReal*(sum (real(matmul(tensorField_fourier(1:3,1:3,i,j,k),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again
conjg(-xi1st(1:3,i,j,k))*geomSize/scaledGeomSize))**2.0_pReal)& ! --> sum squared L_2 norm of vector
+ 2.0_pReal*(sum (real(matmul(tensorField_fourier(1:3,1:3,i,j,k),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again
conjg(-xi1st(1:3,i,j,k))*rescaledGeom))**2.0_pReal)& ! --> sum squared L_2 norm of vector
+sum(aimag(matmul(tensorField_fourier(1:3,1:3,i,j,k),&
conjg(-xi1st(1:3,i,j,k))*geomSize/scaledGeomSize))**2.0_pReal))
conjg(-xi1st(1:3,i,j,k))*rescaledGeom))**2.0_pReal))
enddo
utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if grid(1) /= 1)
+ sum( real(matmul(tensorField_fourier(1:3,1:3,1 ,j,k), &
conjg(-xi1st(1:3,1,j,k))*geomSize/scaledGeomSize))**2.0_pReal) &
conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2.0_pReal) &
+ sum(aimag(matmul(tensorField_fourier(1:3,1:3,1 ,j,k), &
conjg(-xi1st(1:3,1,j,k))*geomSize/scaledGeomSize))**2.0_pReal) &
conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2.0_pReal) &
+ sum( real(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), &
conjg(-xi1st(1:3,grid1Red,j,k))*geomSize/scaledGeomSize))**2.0_pReal) &
conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2.0_pReal) &
+ sum(aimag(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), &
conjg(-xi1st(1:3,grid1Red,j,k))*geomSize/scaledGeomSize))**2.0_pReal)
conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2.0_pReal)
enddo; enddo
if(grid(1) == 1_pInt) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of grid(1) == 1
call MPI_Allreduce(MPI_IN_PLACE,utilities_divergenceRMS,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='utilities_divergenceRMS')
utilities_divergenceRMS = sqrt(utilities_divergenceRMS) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space
@ -675,6 +678,8 @@ end function utilities_divergenceRMS
!> @brief calculate max of curl of field_fourier
!--------------------------------------------------------------------------------------------------
real(pReal) function utilities_curlRMS()
use IO, only: &
IO_error
use numerics, only: &
worldrank
use mesh, only: &
@ -683,9 +688,9 @@ real(pReal) function utilities_curlRMS()
grid3
implicit none
integer(pInt) :: i, j, k, l
complex(pReal), dimension(3,3) :: curl_fourier
PetscErrorCode :: ierr
integer(pInt) :: i, j, k, l, ierr
complex(pReal), dimension(3,3) :: curl_fourier
complex(pReal), dimension(3) :: rescaledGeom
external :: &
MPI_Reduce, &
@ -695,47 +700,49 @@ real(pReal) function utilities_curlRMS()
write(6,'(/,a)') ' ... calculating curl ......................................................'
flush(6)
endif
rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal)
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
! calculating max curl criterion in Fourier space
utilities_curlRMS = 0.0_pReal
do k = 1_pInt, grid3; do j = 1_pInt, grid(2);
do i = 2_pInt, grid1Red - 1_pInt
do l = 1_pInt, 3_pInt
curl_fourier(l,1) = (+tensorField_fourier(l,3,i,j,k)*xi1st(2,i,j,k)*geomSize(2)/scaledGeomSize(2) &
-tensorField_fourier(l,2,i,j,k)*xi1st(3,i,j,k)*geomSize(3)/scaledGeomSize(3))
curl_fourier(l,2) = (+tensorField_fourier(l,1,i,j,k)*xi1st(3,i,j,k)*geomSize(3)/scaledGeomSize(3) &
-tensorField_fourier(l,3,i,j,k)*xi1st(1,i,j,k)*geomSize(1)/scaledGeomSize(1))
curl_fourier(l,3) = (+tensorField_fourier(l,2,i,j,k)*xi1st(1,i,j,k)*geomSize(1)/scaledGeomSize(1) &
-tensorField_fourier(l,1,i,j,k)*xi1st(2,i,j,k)*geomSize(2)/scaledGeomSize(2))
curl_fourier(l,1) = (+tensorField_fourier(l,3,i,j,k)*xi1st(2,i,j,k)*rescaledGeom(2) &
-tensorField_fourier(l,2,i,j,k)*xi1st(3,i,j,k)*rescaledGeom(3))
curl_fourier(l,2) = (+tensorField_fourier(l,1,i,j,k)*xi1st(3,i,j,k)*rescaledGeom(3) &
-tensorField_fourier(l,3,i,j,k)*xi1st(1,i,j,k)*rescaledGeom(1))
curl_fourier(l,3) = (+tensorField_fourier(l,2,i,j,k)*xi1st(1,i,j,k)*rescaledGeom(1) &
-tensorField_fourier(l,1,i,j,k)*xi1st(2,i,j,k)*rescaledGeom(2))
enddo
utilities_curlRMS = utilities_curlRMS + &
2.0_pReal*sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal)! Has somewhere a conj. complex counterpart. Therefore count it twice.
enddo
do l = 1_pInt, 3_pInt
curl_fourier = (+tensorField_fourier(l,3,1,j,k)*xi1st(2,1,j,k)*geomSize(2)/scaledGeomSize(2) &
-tensorField_fourier(l,2,1,j,k)*xi1st(3,1,j,k)*geomSize(3)/scaledGeomSize(3))
curl_fourier = (+tensorField_fourier(l,1,1,j,k)*xi1st(3,1,j,k)*geomSize(3)/scaledGeomSize(3) &
-tensorField_fourier(l,3,1,j,k)*xi1st(1,1,j,k)*geomSize(1)/scaledGeomSize(1))
curl_fourier = (+tensorField_fourier(l,2,1,j,k)*xi1st(1,1,j,k)*geomSize(1)/scaledGeomSize(1) &
-tensorField_fourier(l,1,1,j,k)*xi1st(2,1,j,k)*geomSize(2)/scaledGeomSize(2))
curl_fourier = (+tensorField_fourier(l,3,1,j,k)*xi1st(2,1,j,k)*rescaledGeom(2) &
-tensorField_fourier(l,2,1,j,k)*xi1st(3,1,j,k)*rescaledGeom(3))
curl_fourier = (+tensorField_fourier(l,1,1,j,k)*xi1st(3,1,j,k)*rescaledGeom(3) &
-tensorField_fourier(l,3,1,j,k)*xi1st(1,1,j,k)*rescaledGeom(1))
curl_fourier = (+tensorField_fourier(l,2,1,j,k)*xi1st(1,1,j,k)*rescaledGeom(1) &
-tensorField_fourier(l,1,1,j,k)*xi1st(2,1,j,k)*rescaledGeom(2))
enddo
utilities_curlRMS = utilities_curlRMS + &
sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal)! this layer (DC) does not have a conjugate complex counterpart (if grid(1) /= 1)
do l = 1_pInt, 3_pInt
curl_fourier = (+tensorField_fourier(l,3,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*geomSize(2)/scaledGeomSize(2) &
-tensorField_fourier(l,2,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*geomSize(3)/scaledGeomSize(3))
curl_fourier = (+tensorField_fourier(l,1,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*geomSize(3)/scaledGeomSize(3) &
-tensorField_fourier(l,3,grid1Red,j,k)*xi1st(1,grid1Red,j,k)*geomSize(1)/scaledGeomSize(1))
curl_fourier = (+tensorField_fourier(l,2,grid1Red,j,k)*xi1st(1,grid1Red,j,k)*geomSize(1)/scaledGeomSize(1) &
-tensorField_fourier(l,1,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*geomSize(2)/scaledGeomSize(2))
curl_fourier = (+tensorField_fourier(l,3,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2) &
-tensorField_fourier(l,2,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*rescaledGeom(3))
curl_fourier = (+tensorField_fourier(l,1,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*rescaledGeom(3) &
-tensorField_fourier(l,3,grid1Red,j,k)*xi1st(1,grid1Red,j,k)*rescaledGeom(1))
curl_fourier = (+tensorField_fourier(l,2,grid1Red,j,k)*xi1st(1,grid1Red,j,k)*rescaledGeom(1) &
-tensorField_fourier(l,1,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2))
enddo
utilities_curlRMS = utilities_curlRMS + &
sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal)! this layer (Nyquist) does not have a conjugate complex counterpart (if grid(1) /= 1)
enddo; enddo
call MPI_Allreduce(MPI_IN_PLACE,utilities_curlRMS,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='utilities_curlRMS')
utilities_curlRMS = sqrt(utilities_curlRMS) * wgt
if(grid(1) == 1_pInt) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of grid(1) == 1
@ -931,6 +938,10 @@ end subroutine utilities_fourierTensorDivergence
!--------------------------------------------------------------------------------------------------
subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
P,C_volAvg,C_minmaxAvg,P_av,forwardData,rotation_BC)
use prec, only: &
dNeq
use IO, only: &
IO_error
use debug, only: &
debug_reset, &
debug_info
@ -969,10 +980,9 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
age
integer(pInt) :: &
j,k
j,k,ierr
real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF
real(pReal) :: max_dPdF_norm, min_dPdF_norm, defgradDetMin, defgradDetMax, defgradDet
PetscErrorCode :: ierr
external :: &
MPI_Reduce, &
@ -1006,7 +1016,9 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
defgradDetMin = min(defgradDetMin,defgradDet)
end do
call MPI_reduce(MPI_IN_PLACE,defgradDetMax,1,MPI_DOUBLE,MPI_MAX,0,PETSC_COMM_WORLD,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce max')
call MPI_reduce(MPI_IN_PLACE,defgradDetMin,1,MPI_DOUBLE,MPI_MIN,0,PETSC_COMM_WORLD,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce min')
if (worldrank == 0_pInt) then
write(6,'(a,1x,es11.4)') ' max determinant of deformation =', defgradDetMax
write(6,'(a,1x,es11.4)') ' min determinant of deformation =', defgradDetMin
@ -1032,7 +1044,9 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
end do
call MPI_Allreduce(MPI_IN_PLACE,max_dPdF,81,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce max')
call MPI_Allreduce(MPI_IN_PLACE,min_dPdF,81,MPI_DOUBLE,MPI_MIN,PETSC_COMM_WORLD,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce min')
C_minmaxAvg = 0.5_pReal*(max_dPdF + min_dPdF)
C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt
@ -1187,6 +1201,10 @@ end function utilities_getFreqDerivative
! convolution
!--------------------------------------------------------------------------------------------------
subroutine utilities_updateIPcoords(F)
use prec, only: &
cNeq
use IO, only: &
IO_error
use math, only: &
math_mul33x3
use mesh, only: &
@ -1198,10 +1216,9 @@ subroutine utilities_updateIPcoords(F)
implicit none
real(pReal), dimension(3,3,grid(1),grid(2),grid3), intent(in) :: F
integer(pInt) :: i, j, k, m
integer(pInt) :: i, j, k, m, ierr
real(pReal), dimension(3) :: step, offset_coords
real(pReal), dimension(3,3) :: Favg
PetscErrorCode :: ierr
external &
MPI_Bcast
@ -1212,8 +1229,8 @@ subroutine utilities_updateIPcoords(F)
call utilities_FFTtensorForward()
call utilities_fourierTensorDivergence()
do k = 1_pInt, grid3; do j = 1_pInt, grid(2) ;do i = 1_pInt, grid1Red
if (any(abs(xi1st(1:3,i,j,k)) > tiny(0.0_pReal))) &
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid1Red
if (any(cNeq(xi1st(1:3,i,j,k),cmplx(0.0_pReal,0.0_pReal)))) &
vectorField_fourier(1:3,i,j,k) = vectorField_fourier(1:3,i,j,k)/ &
sum(conjg(-xi1st(1:3,i,j,k))*xi1st(1:3,i,j,k))
enddo; enddo; enddo
@ -1223,12 +1240,14 @@ subroutine utilities_updateIPcoords(F)
! average F
if (grid3Offset == 0_pInt) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt
call MPI_Bcast(Favg,9,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='update_IPcoords')
!--------------------------------------------------------------------------------------------------
! add average to fluctuation and put (0,0,0) on (0,0,0)
step = geomSize/real(grid, pReal)
if (grid3Offset == 0_pInt) offset_coords = vectorField_real(1:3,1,1,1)
call MPI_Bcast(offset_coords,3,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='update_IPcoords')
offset_coords = math_mul33x3(Favg,step/2.0_pReal) - offset_coords
m = 1_pInt
do k = 1_pInt,grid3; do j = 1_pInt,grid(2); do i = 1_pInt,grid(1)

89
code/system_routines.f90 Normal file
View File

@ -0,0 +1,89 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief provides wrappers to C routines
!--------------------------------------------------------------------------------------------------
module system_routines
implicit none
private
public :: &
isDirectory, &
getCWD
interface
function isDirectory_C(path) BIND(C)
use, intrinsic :: ISO_C_Binding, only: &
C_INT, &
C_CHAR
integer(C_INT) :: isDirectory_C
character(kind=C_CHAR), dimension(1024), intent(in) :: path ! C string is an array
end function isDirectory_C
subroutine getCurrentWorkDir_C(str, stat) bind(C)
use, intrinsic :: ISO_C_Binding, only: &
C_INT, &
C_CHAR
character(kind=C_CHAR), dimension(1024), intent(out) :: str ! C string is an array
integer(C_INT),intent(out) :: stat
end subroutine getCurrentWorkDir_C
end interface
contains
!--------------------------------------------------------------------------------------------------
!> @brief figures out if a given path is a directory (and not an ordinary file)
!--------------------------------------------------------------------------------------------------
logical function isDirectory(path)
use, intrinsic :: ISO_C_Binding, only: &
C_INT, &
C_CHAR, &
C_NULL_CHAR
implicit none
character(len=*), intent(in) :: path
character(kind=C_CHAR), dimension(1024) :: strFixedLength
integer :: i
strFixedLength = repeat(C_NULL_CHAR,len(strFixedLength))
do i=1,len(path) ! copy array components
strFixedLength(i)=path(i:i)
enddo
isDirectory=merge(.True.,.False.,isDirectory_C(strFixedLength) /= 0_C_INT)
end function isDirectory
!--------------------------------------------------------------------------------------------------
!> @brief gets the current working directory
!--------------------------------------------------------------------------------------------------
logical function getCWD(str)
use, intrinsic :: ISO_C_Binding, only: &
C_INT, &
C_CHAR, &
C_NULL_CHAR
implicit none
character(len=*), intent(out) :: str
character(kind=C_CHAR), dimension(1024) :: strFixedLength ! C string is an array
integer(C_INT) :: stat
integer :: i
str = repeat('',len(str))
call getCurrentWorkDir_C(strFixedLength,stat)
do i=1,1024 ! copy array components until Null string is found
if (strFixedLength(i) /= C_NULL_CHAR) then
str(i:i)=strFixedLength(i)
else
exit
endif
enddo
getCWD=merge(.True.,.False.,stat /= 0_C_INT)
end function getCWD
end module system_routines

12
code/thermal_adiabatic.f90
View File

@ -236,7 +236,7 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
use material, only: &
homogenization_Ngrains, &
mappingHomogenization, &
phaseAt, phasememberAt, &
phaseAt, &
thermal_typeInstance, &
phase_Nsources, &
phase_source, &
@ -297,8 +297,8 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
enddo
enddo
Tdot = Tdot/homogenization_Ngrains(homog)
dTdot_dT = dTdot_dT/homogenization_Ngrains(homog)
Tdot = Tdot/real(homogenization_Ngrains(homog),pReal)
dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal)
end subroutine thermal_adiabatic_getSourceAndItsTangent
@ -336,8 +336,7 @@ function thermal_adiabatic_getSpecificHeat(ip,el)
enddo
thermal_adiabatic_getSpecificHeat = &
thermal_adiabatic_getSpecificHeat/ &
homogenization_Ngrains(mesh_element(3,el))
thermal_adiabatic_getSpecificHeat/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
end function thermal_adiabatic_getSpecificHeat
@ -375,8 +374,7 @@ function thermal_adiabatic_getMassDensity(ip,el)
enddo
thermal_adiabatic_getMassDensity = &
thermal_adiabatic_getMassDensity/ &
homogenization_Ngrains(mesh_element(3,el))
thermal_adiabatic_getMassDensity/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
end function thermal_adiabatic_getMassDensity

15
code/thermal_conduction.f90
View File

@ -190,7 +190,7 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
use material, only: &
homogenization_Ngrains, &
mappingHomogenization, &
phaseAt, phasememberAt, &
phaseAt, &
thermal_typeInstance, &
phase_Nsources, &
phase_source, &
@ -252,8 +252,8 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
enddo
enddo
Tdot = Tdot/homogenization_Ngrains(homog)
dTdot_dT = dTdot_dT/homogenization_Ngrains(homog)
Tdot = Tdot/real(homogenization_Ngrains(homog),pReal)
dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal)
end subroutine thermal_conduction_getSourceAndItsTangent
@ -291,8 +291,7 @@ function thermal_conduction_getConductivity33(ip,el)
enddo
thermal_conduction_getConductivity33 = &
thermal_conduction_getConductivity33/ &
homogenization_Ngrains(mesh_element(3,el))
thermal_conduction_getConductivity33/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
end function thermal_conduction_getConductivity33
@ -330,8 +329,7 @@ function thermal_conduction_getSpecificHeat(ip,el)
enddo
thermal_conduction_getSpecificHeat = &
thermal_conduction_getSpecificHeat/ &
homogenization_Ngrains(mesh_element(3,el))
thermal_conduction_getSpecificHeat/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
end function thermal_conduction_getSpecificHeat
@ -369,8 +367,7 @@ function thermal_conduction_getMassDensity(ip,el)
enddo
thermal_conduction_getMassDensity = &
thermal_conduction_getMassDensity/ &
homogenization_Ngrains(mesh_element(3,el))
thermal_conduction_getMassDensity/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
end function thermal_conduction_getMassDensity

21
code/vacancyflux_cahnhilliard.f90
View File

@ -219,7 +219,7 @@ subroutine vacancyflux_cahnhilliard_getSourceAndItsTangent(CvDot, dCvDot_dCv, Cv
use material, only: &
homogenization_Ngrains, &
mappingHomogenization, &
phaseAt, phasememberAt, &
phaseAt, &
phase_source, &
phase_Nsources, &
SOURCE_vacancy_phenoplasticity_ID, &
@ -266,8 +266,8 @@ subroutine vacancyflux_cahnhilliard_getSourceAndItsTangent(CvDot, dCvDot_dCv, Cv
enddo
enddo
CvDot = CvDot/homogenization_Ngrains(mappingHomogenization(2,ip,el))
dCvDot_dCv = dCvDot_dCv/homogenization_Ngrains(mappingHomogenization(2,ip,el))
CvDot = CvDot/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal)
dCvDot_dCv = dCvDot_dCv/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal)
end subroutine vacancyflux_cahnhilliard_getSourceAndItsTangent
@ -301,8 +301,7 @@ function vacancyflux_cahnhilliard_getMobility33(ip,el)
enddo
vacancyflux_cahnhilliard_getMobility33 = &
vacancyflux_cahnhilliard_getMobility33/ &
homogenization_Ngrains(mesh_element(3,el))
vacancyflux_cahnhilliard_getMobility33/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
end function vacancyflux_cahnhilliard_getMobility33
@ -336,8 +335,7 @@ function vacancyflux_cahnhilliard_getDiffusion33(ip,el)
enddo
vacancyflux_cahnhilliard_getDiffusion33 = &
vacancyflux_cahnhilliard_getDiffusion33/ &
homogenization_Ngrains(mesh_element(3,el))
vacancyflux_cahnhilliard_getDiffusion33/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
end function vacancyflux_cahnhilliard_getDiffusion33
@ -371,8 +369,7 @@ real(pReal) function vacancyflux_cahnhilliard_getFormationEnergy(ip,el)
enddo
vacancyflux_cahnhilliard_getFormationEnergy = &
vacancyflux_cahnhilliard_getFormationEnergy/ &
homogenization_Ngrains(mesh_element(3,el))
vacancyflux_cahnhilliard_getFormationEnergy/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
end function vacancyflux_cahnhilliard_getFormationEnergy
@ -408,7 +405,7 @@ real(pReal) function vacancyflux_cahnhilliard_getEntropicCoeff(ip,el)
vacancyflux_cahnhilliard_getEntropicCoeff = &
vacancyflux_cahnhilliard_getEntropicCoeff* &
temperature(material_homog(ip,el))%p(thermalMapping(material_homog(ip,el))%p(ip,el))/ &
homogenization_Ngrains(material_homog(ip,el))
real(homogenization_Ngrains(material_homog(ip,el)),pReal)
end function vacancyflux_cahnhilliard_getEntropicCoeff
@ -467,8 +464,8 @@ subroutine vacancyflux_cahnhilliard_KinematicChemPotAndItsTangent(KPot, dKPot_dC
enddo
enddo
KPot = KPot/homogenization_Ngrains(material_homog(ip,el))
dKPot_dCv = dKPot_dCv/homogenization_Ngrains(material_homog(ip,el))
KPot = KPot/real(homogenization_Ngrains(material_homog(ip,el)),pReal)
dKPot_dCv = dKPot_dCv/real(homogenization_Ngrains(material_homog(ip,el)),pReal)
end subroutine vacancyflux_cahnhilliard_KinematicChemPotAndItsTangent

6
code/vacancyflux_isochempot.f90
View File

@ -235,7 +235,7 @@ subroutine vacancyflux_isochempot_getSourceAndItsTangent(CvDot, dCvDot_dCv, Cv,
use material, only: &
homogenization_Ngrains, &
mappingHomogenization, &
phaseAt, phasememberAt, &
phaseAt, &
phase_source, &
phase_Nsources, &
SOURCE_vacancy_phenoplasticity_ID, &
@ -282,8 +282,8 @@ subroutine vacancyflux_isochempot_getSourceAndItsTangent(CvDot, dCvDot_dCv, Cv,
enddo
enddo
CvDot = CvDot/homogenization_Ngrains(mappingHomogenization(2,ip,el))
dCvDot_dCv = dCvDot_dCv/homogenization_Ngrains(mappingHomogenization(2,ip,el))
CvDot = CvDot/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal)
dCvDot_dCv = dCvDot_dCv/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal)
end subroutine vacancyflux_isochempot_getSourceAndItsTangent

1
config/Crystallite_All.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[all]
(output) phase
(output) texture

1
config/Crystallite_None.config
View File

@ -1,2 +1 @@
### $Id$ ###
[none]

1
config/Crystallite_aLittleSomething.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[aLittleSomething]
(output) f
(output) p

1
config/Homogenization_Damage_NonLocal.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
damage nonlocal
initialDamage 1.0
(output) damage

1
config/Homogenization_HydrogenFlux_CahnHilliard.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
hydrogenflux cahnhilliard
initialHydrogenConc 0.0
(output) hydrogenconc

1
config/Homogenization_Isostrain_Parallel3.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[Parallel3]
type isostrain
Ngrains 3

1
config/Homogenization_Isostrain_SX.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[SX]
type isostrain
Ngrains 1

1
config/Homogenization_Isostrain_Taylor2.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[Taylor2]
type isostrain
Ngrains 2

1
config/Homogenization_None_Dummy.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[directSX]
type none

1
config/Homogenization_Porosity_PhaseField.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
porosity phasefield
(output) porosity

1
config/Homogenization_RGC_8Grains.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[8Grains]
type RGC
Ngrains 8

1
config/Homogenization_Thermal_Conduction.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
thermal conduction
initialT 300.0
(output) temperature

1
config/Homogenization_VacancyFlux_CahnHilliard.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
vacancyflux cahnhilliard
initialVacancyConc 1e-6
(output) vacancyconc

1
config/Homogenization_multiField.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[SX]
type isostrain
Ngrains 1

1
config/Kinematics_Hydrogen_Strain.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
(kinematics) vacancy_strain
vacancy_strain_coeff 0.006

1
config/Kinematics_Thermal_Expansion.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
(kinematics) thermal_expansion
thermal_expansion11 0.00231

1
config/Kinematics_Vacancy_Strain.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
(kinematics) hydrogen_strain
hydrogen_strain_coeff 0.06

1
config/Microstructure_DP_Steel.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[DP_Steel]
/elementhomogeneous/
crystallite 1

1
config/Microstructure_ElementHomogeneous.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[ElementHomogeneous]
/elementhomogeneous/ # put this flag to set ips identical in one element (something like reduced integration)
crystallite 1

1
config/Phase_Damage.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
damage_diffusion11 1.0
damage_mobility 0.001

1
config/Phase_Dislotwin_TWIP-Steel-FeMnC.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[TWIP_Steel_FeMnC]
elasticity hooke

1
config/Phase_Dislotwin_Tungsten.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[Tungsten]
elasticity hooke

1
config/Phase_Hydrogen.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
hydrogenflux_diffusion11 1.0
hydrogenflux_mobility11 1.0
hydrogenVolume 1e-28

3
config/Phase_J2_AluminumIsotropic.config → config/Phase_Isotropic_AluminumIsotropic.config
View File

@ -1,9 +1,8 @@
### $Id$ ###
[Aluminum_Isotropic]
# Kuo, J. C., Mikrostrukturmechanik von Bikristallen mit Kippkorngrenzen. Shaker-Verlag 2004. http://edoc.mpg.de/204079
elasticity hooke
plasticity j2
plasticity isotropic
(output) flowstress
(output) strainrate

1
config/Phase_None_IsotropicVolumePreservation.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[IsotropicVolumePreservation]
elasticity hooke

1
config/Phase_None_Orthorombic.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[Orthorombic]
elasticity hooke

1
config/Phase_Nonlocal_Aluminum.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[Aluminum]
elasticity hooke

1
config/Phase_Phenopowerlaw_Aluminum.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[Aluminum]
elasticity hooke
plasticity phenopowerlaw

1
config/Phase_Phenopowerlaw_BCC-Ferrite.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
# Tasan et.al. 2015 Acta Materalia
# Tasan et.al. 2015 International Journal of Plasticity
# Diehl et.al. 2015 Meccanica

1
config/Phase_Phenopowerlaw_BCC-Martensite.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
# Tasan et.al. 2015 Acta Materalia
# Tasan et.al. 2015 International Journal of Plasticity
# Diehl et.al. 2015 Meccanica

1
config/Phase_Phenopowerlaw_Gold.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
# parameters fitted by D. Ma to:

1
config/Phase_Phenopowerlaw_cpTi-alpha.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[cpTi-alpha]
plasticity phenopowerlaw
elasticity hooke

1
config/Phase_Phenopowerlaw_multiField.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[Aluminum]
elasticity hooke
plasticity phenopowerlaw

1
config/Phase_Porosity.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
porosity_diffusion11 1.0
porosity_mobility 0.001

1
config/Phase_Thermal.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
thermal_conductivity11 237.0
specific_heat 910.0
mass_density 2700.0

1
config/Phase_Vacancy.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
vacancyflux_diffusion11 1.0
vacancyflux_mobility11 1.0
vacancyFormationEnergy 1e-19

1
config/Source_Damage_IsoBrittle.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
(source) damage_isoBrittle
isobrittle_criticalStrainEnergy 1400000.0
isobrittle_atol 0.01

1
config/Source_Thermal_Dissipation.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
(source) thermal_dissipation
dissipation_ColdWorkCoeff 0.95

1
config/Source_Vacancy_Irradiation.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
(source) vacancy_irradiation
irradiation_cascadeprobability 0.00001
irradiation_cascadevolume 1000.0

1
config/Source_Vacancy_PhenoPlasticity.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
(source) vacancy_phenoplasticity
phenoplasticity_ratecoeff 0.01

1
config/Texture_FiberExample.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[FiberExample]
axes x y -z # model coordinate x-, y-, z-axes correspond to which axes during texture measurement? this was a left handed coordinate system!
# fiber axis in spherical coordinates: alpha crystal system, beta sample system

1
config/Texture_Gauss_001.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
[001]
(gauss) phi1 0.000 Phi 0.000 phi2 0.000 scatter 0.000 fraction 1.000

1
config/Texture_Gauss_101.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
[101]
(gauss) phi1 0.000 Phi 45.000 phi2 90.000 scatter 0.000 fraction 1.000

1
config/Texture_Gauss_111.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
[111]
(gauss) phi1 0.000 Phi 54.7356 phi2 45.000 scatter 0.000 fraction 1.000

1
config/Texture_Gauss_123.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
[123]
(gauss) phi1 209.805 Phi 29.206 phi2 63.435 scatter 0.000 fraction 1.000

1
config/Texture_RandomSingleCrystals.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
[RandomSingleCrystals]
(random) scatter 0.000 fraction 1.000

1
config/Texture_Rolling.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[Rolling]
hybridIA rollingTexture.linearODF
symmetry orthotropic # or monoclinic

1
config/debug.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
### debugging parameters ###
# example:

1
config/material.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
### numerical parameters ###
# The material.config file needs to specify five parts:

1
config/numerics.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
### numerical parameters ###
relevantStrain 1.0e-7 # strain increment considered significant (used by crystallite to determine whether strain inc is considered significant)

373258
config/rollingTexture.linearODF
View File

File diff suppressed because it is too large Load Diff

1
examples/.gitignore vendored Normal file
View File

@ -0,0 +1 @@
postProc

4
examples/AbaqusStandard/material.config
View File

@ -1,7 +1,3 @@
#####################
# $Id$
#####################
#-------------------#
<homogenization>
#-------------------#

4
examples/AbaqusStandard/numerics.config
View File

@ -1,2 +1,2 @@
fixed_seed 144921823
analyticJaco 0
fixed_seed 1697667030
analyticJaco 1

4
examples/MSC.Marc/material.config
View File

@ -1,7 +1,3 @@
#####################
# $Id$
#####################
#-------------------#
<homogenization>
#-------------------#

33
examples/SpectralMethod/example_postprocess.py
View File

@ -1,33 +0,0 @@
#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import os
import glob
from subprocess import call
geom_name = '20grains16x16x16_tensionX'
postResults = 'postResults --cr f,p --split --separation x,y,z '+geom_name+'.spectralOut'
sts = call(postResults, shell=True)
os.chdir('./postProc/')
ascii_files = glob.glob(geom_name+'_inc*.txt')
print ascii_files
showTable = "showTable -a "
addCauchy = 'addCauchy '
addMises = 'addMises -s Cauchy '
addStrainTensors = "addStrainTensors -0 -v "
visualize3D = "3Dvisualize -s 'Mises(Cauchy)',1_p Cauchy "
postProc = [addCauchy, addMises, addStrainTensors, visualize3D]
for f in ascii_files:
print f
for p in postProc:
p = p+f
print p
sts = call(p,shell=True)

17
examples/SpectralMethod/example_run.py
View File

@ -1,17 +0,0 @@
#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import sys
resolutions = [16,32,64]
resolution = resolutions[0]
try:
resolution = int(sys.argv[1])
except:
pass
if resolution not in resolutions:
resolution = resolutions[0]
from subprocess import call
call('make run%s'%('x'.join([str(resolution)]*3)), shell=True)

14
installation/compile_CoreModule.py
View File

@ -1,7 +1,7 @@
#!/usr/bin/env python
#!/usr/bin/env python2
# -*- coding: UTF-8 no BOM -*-
import os,sys,glob,string,subprocess,shlex
import os,sys,glob,subprocess,shlex
from damask import Environment
from damask import version as DAMASKVERSION
@ -47,9 +47,13 @@ compileOptions = ' -DSpectral -DFLOAT=8 -DINT=4 -I%s/lib -DDAMASKVERSION=\\\\\"\
#--- this saves the path of libraries to core.so, hence it is known during runtime ----------------
if options['F90'] == 'gfortran':
LDFLAGS = '-shared -Wl,-undefined,dynamic_lookup' # solved error: Undefined symbols for architecture x86_64: "_PyArg_ParseTupleAndKeywords" as found on https://lists.macosforge.org/pipermail/macports-dev/2013-May/022735.html
# solved error: Undefined symbols for architecture x86_64: "_PyArg_ParseTupleAndKeywords"
# as found on https://lists.macosforge.org/pipermail/macports-dev/2013-May/022735.html
LDFLAGS = '-shared -Wl,-undefined,dynamic_lookup'
else:
LDFLAGS = ' -openmp -Wl' # some f2py versions/configurations compile with openMP, so linking against openMP is needed to prevent errors during loading of core module
# some f2py versions/configurations compile with openMP, so linking against openMP is needed
# to prevent errors during loading of core module
LDFLAGS = ' -openmp -Wl'
#--- run path of for fftw during runtime ----------------------------------------------------------
LDFLAGS += ',-rpath,%s/lib,-rpath,%s/lib64'%(options['FFTW_ROOT'],options['FFTW_ROOT'])
@ -96,6 +100,8 @@ os.chdir(codeDir)
cmd = 'f2py damask.core.pyf' +\
' -c --no-lower %s'%(compiler) +\
compileOptions+\
' C_routines.c'+\
' system_routines.f90'+\
' prec.f90'+\
' spectral_interface.f90'+\
' IO.f90'+\

1
installation/mods_Abaqus/abaqus_v6.env
View File

@ -1,7 +1,6 @@
#
# DAMASK Abaqus Environment File
#
# $Id$
# ------------------------------------
# originally taken from Abaqus ver. 6.11.1
#

1
installation/mods_Abaqus/abaqus_v6_serial.env
View File

@ -1,7 +1,6 @@
#
# DAMASK Abaqus Environment File
#
# $Id$
# ------------------------------------
# originally taken from Abaqus ver. 6.11.1
#

354
installation/mods_Abaqus/abaqus_v6_windows.env
View File

@ -1,177 +1,177 @@
#
# System-Wide Abaqus Environment File
# -------------------------------------
standard_parallel = ALL
mp_mode = MPI
mp_file_system = (DETECT,DETECT)
mp_num_parallel_ftps = (4, 4)
mp_environment_export = ('MPI_PROPAGATE_TSTP',
'ABA_CM_BUFFERING',
'ABA_CM_BUFFERING_LIMIT',
'ABA_ITERATIVE_SOLVER_VERBOSE',
'ABA_DMPSOLVER_BWDPARALLELOFF',
'ABA_ELP_SURFACE_SPLIT',
'ABA_ELP_SUSPEND',
'ABA_HOME',
'ABA_MEMORY_MODE',
'ABA_MPI_MESSAGE_TRACKING',
'ABA_MPI_VERBOSE_LEVEL',
'ABA_PATH',
'ABAQUS_CSE_RELTIMETOLERANCE',
'ABA_RESOURCE_MONITOR',
'ABA_RESOURCE_USEMALLINFO',
'ABAQUS_LANG',
'ABAQUS_CSE_CURRCONFIGMAPPING',
'ABAQUS_MPF_DIAGNOSTIC_LEVEL',
'ABAQUSLM_LICENSE_FILE',
'ABQ_CRTMALLOC',
'ABQ_DATACHECK',
'ABQ_RECOVER',
'ABQ_RESTART',
'ABQ_SPLITFILE',
'ABQ_XPL_WINDOWDUMP',
'ABQ_XPL_PARTITIONSIZE',
'ABQLMHANGLIMIT',
'ABQLMQUEUE',
'ABQLMUSER',
'CCI_RENDEZVOUS',
'DOMAIN',
'DOMAIN_CPUS',
'DOUBLE_PRECISION',
'FLEXLM_DIAGNOSTICS',
'FOR0006',
'FOR0064',
'FOR_IGNORE_EXCEPTIONS',
'FOR_DISABLE_DIAGNOSTIC_DISPLAY',
'LD_PRELOAD',
'MP_NUMBER_OF_THREADS',
'MPC_GANG',
'MPI_FLAGS',
'MPI_FLUSH_FCACHE',
'MPI_RDMA_NENVELOPE',
'MPI_SOCKBUFSIZE',
'MPI_USE_MALLOPT_MMAP_MAX',
'MPI_USE_MALLOPT_MMAP_THRESHOLD',
'MPI_USE_MALLOPT_SBRK_PROTECTION',
'MPI_WORKDIR',
'MPCCI_DEBUG',
'MPCCI_CODEID',
'MPCCI_JOBID',
'MPCCI_NETDEVICE',
'MPCCI_TINFO',
'MPCCI_SERVER',
'ABAQUS_CCI_DEBUG',
'NCPUS',
'OMP_DYNAMIC',
'OMP_NUM_THREADS',
'OUTDIR',
'PAIDUP',
'PARALLEL_METHOD',
'RAIDEV_NDREG_LAZYMEM',
'ABA_SYMBOLIC_GENERALCOLLAPSE',
'ABA_SYMBOLIC_GENERAL_MAXCLIQUERANK',
'ABA_ADM_MINIMUMINCREASE',
'ABA_ADM_MINIMUMDECREASE',
'IPATH_NO_CPUAFFINITY',
'MALLOC_MMAP_THRESHOLD_',
'ABA_EXT_SIMOUTPUT',
'SMA_WS',
'SMA_PARENT',
'SMA_PLATFORM',
'ABA_PRE_DECOMPOSITION',
'ACML_FAST_MALLOC',
'ACML_FAST_MALLOC_CHUNK_SIZE',
'ACML_FAST_MALLOC_MAX_CHUNKS',
'ACML_FAST_MALLOC_DEBUG')
import driverUtils, os
#-*- mode: python -*-
# #
# Compile and Link command settings for the Windows 64 Platform #
# ( AMD Opteron / Intel EM64T ) #
# #
compile_fortran=['ifort',
'/c','/DABQ_WIN86_64', '/u',
'/iface:cref', '/recursive', '/Qauto-scalar',
'/QxSSE3', '/QaxAVX',
'/heap-arrays:1',
# '/Od', '/Ob0' # <-- Optimization
# '/Zi', # <-- Debugging
'/include:%I', '/free', '/O1', '/fpp', '/openmp', '/Qmkl']
link_sl=['LINK',
'/nologo', '/NOENTRY', '/INCREMENTAL:NO', '/subsystem:console', '/machine:AMD64',
'/NODEFAULTLIB:LIBC.LIB', '/NODEFAULTLIB:LIBCMT.LIB',
'/DEFAULTLIB:OLDNAMES.LIB', '/DEFAULTLIB:LIBIFCOREMD.LIB', '/DEFAULTLIB:LIBIFPORTMD', '/DEFAULTLIB:LIBMMD.LIB',
'/DEFAULTLIB:kernel32.lib', '/DEFAULTLIB:user32.lib', '/DEFAULTLIB:advapi32.lib',
'/FIXED:NO', '/dll',
'/def:%E', '/out:%U', '%F', '%A', '%L', '%B',
'oldnames.lib', 'user32.lib', 'ws2_32.lib', 'netapi32.lib', 'advapi32.lib']
link_exe=['LINK',
'/nologo', '/INCREMENTAL:NO', '/subsystem:console', '/machine:AMD64', '/STACK:20000000',
'/NODEFAULTLIB:LIBC.LIB', '/NODEFAULTLIB:LIBCMT.LIB', '/DEFAULTLIB:OLDNAMES.LIB', '/DEFAULTLIB:LIBIFCOREMD.LIB',
'/DEFAULTLIB:LIBIFPORTMD', '/DEFAULTLIB:LIBMMD.LIB', '/DEFAULTLIB:kernel32.lib',
'/DEFAULTLIB:user32.lib', '/DEFAULTLIB:advapi32.lib',
'/FIXED:NO', '/LARGEADDRESSAWARE',
'/out:%J', '%F', '%M', '%L', '%B', '%O',
'oldnames.lib', 'user32.lib', 'ws2_32.lib', 'netapi32.lib', 'advapi32.lib']
# Link command to be used for MAKE w/o fortran compiler.
# remove the pound signs in order to remove the comments and have the file take effect.
#
#link_exe=['LINK', '/nologo', 'INCREMENTAL:NO', '/subsystem:console', '/machine:AMD64', '/NODEFAULTLIB:LIBC.LIB', '/NODEFAULTLIB:LIBCMT.LIB',
# '/DEFAULTLIB:OLDNAMES.LIB', '/DEFAULTLIB:MSVCRT.LIB', '/DEFAULTLIB:kernel32.lib', 'DEFAULTLIB:user32.lib', '/DEFAULTLIB:advapi32.lib',
# '/FIXED:NO', '/LARGEADDRESSAWARE', '/DEBUG', '/out:%J', '%F', '%M', '%L', '%B', '%O', 'oldnames.lib', 'user32.lib', 'ws2_32.lib',
# 'netapi32.lib', 'advapi32.lib]
# MPI Configuration
mp_mode = THREADS
mp_mpi_implementation = NATIVE
mp_rsh_command = 'dummy %H -l %U -n %C'
mp_mpirun_path = {}
mpirun = ''
progDir = os.environ.get('ProgramFiles','C:\\Program Files')
for mpiDir in ('Microsoft HPC Pack', 'Microsoft HPC Pack 2008 R2', 'Microsoft HPC Pack 2008', 'Microsoft HPC Pack 2008 SDK'):
mpirun = progDir + os.sep + mpiDir + os.sep + 'bin' + os.sep + 'mpiexec.exe'
if os.path.exists(mpirun):
mp_mpirun_path[NATIVE] = mpirun
mp_mpirun_path[MSSDK] = os.path.join(progDir, mpiDir)
break
if os.environ.has_key('CCP_HOME'):
from queueCCS import QueueCCS
queues['default'] = QueueCCS(queueName='share')
queues['share'] = QueueCCS(queueName='share')
queues['local'] = QueueCCS(queueName='local')
queues['genxmlshare'] = QueueCCS(queueName='genxmlshare')
queues['genxmllocal'] = QueueCCS(queueName='genxmllocal')
del QueueCCS
mpirun = os.path.join(os.environ['CCP_HOME'], 'bin', 'mpiexec.exe')
if os.path.exists(mpirun):
mp_mpirun_path[NATIVE] = mpirun
run_mode=BATCH
if mp_mpirun_path:
mp_mode=MPI
del progDir, mpiDir, mpirun
graphicsEnv = driverUtils.locateFile(os.environ['ABA_PATH'],'site','graphicsConfig','env')
if graphicsEnv:
execfile(graphicsEnv)
else:
raise 'Cannot find the graphics configuration environment file (graphicsConfig.env)'
del driverUtils, os, graphicsEnv
license_server_type=FLEXNET
abaquslm_license_file=""
doc_root="
doc_root_type="html"
academic=RESEARCH
#
# System-Wide Abaqus Environment File
# -------------------------------------
standard_parallel = ALL
mp_mode = MPI
mp_file_system = (DETECT,DETECT)
mp_num_parallel_ftps = (4, 4)
mp_environment_export = ('MPI_PROPAGATE_TSTP',
'ABA_CM_BUFFERING',
'ABA_CM_BUFFERING_LIMIT',
'ABA_ITERATIVE_SOLVER_VERBOSE',
'ABA_DMPSOLVER_BWDPARALLELOFF',
'ABA_ELP_SURFACE_SPLIT',
'ABA_ELP_SUSPEND',
'ABA_HOME',
'ABA_MEMORY_MODE',
'ABA_MPI_MESSAGE_TRACKING',
'ABA_MPI_VERBOSE_LEVEL',
'ABA_PATH',
'ABAQUS_CSE_RELTIMETOLERANCE',
'ABA_RESOURCE_MONITOR',
'ABA_RESOURCE_USEMALLINFO',
'ABAQUS_LANG',
'ABAQUS_CSE_CURRCONFIGMAPPING',
'ABAQUS_MPF_DIAGNOSTIC_LEVEL',
'ABAQUSLM_LICENSE_FILE',
'ABQ_CRTMALLOC',
'ABQ_DATACHECK',
'ABQ_RECOVER',
'ABQ_RESTART',
'ABQ_SPLITFILE',
'ABQ_XPL_WINDOWDUMP',
'ABQ_XPL_PARTITIONSIZE',
'ABQLMHANGLIMIT',
'ABQLMQUEUE',
'ABQLMUSER',
'CCI_RENDEZVOUS',
'DOMAIN',
'DOMAIN_CPUS',
'DOUBLE_PRECISION',
'FLEXLM_DIAGNOSTICS',
'FOR0006',
'FOR0064',
'FOR_IGNORE_EXCEPTIONS',
'FOR_DISABLE_DIAGNOSTIC_DISPLAY',
'LD_PRELOAD',
'MP_NUMBER_OF_THREADS',
'MPC_GANG',
'MPI_FLAGS',
'MPI_FLUSH_FCACHE',
'MPI_RDMA_NENVELOPE',
'MPI_SOCKBUFSIZE',
'MPI_USE_MALLOPT_MMAP_MAX',
'MPI_USE_MALLOPT_MMAP_THRESHOLD',
'MPI_USE_MALLOPT_SBRK_PROTECTION',
'MPI_WORKDIR',
'MPCCI_DEBUG',
'MPCCI_CODEID',
'MPCCI_JOBID',
'MPCCI_NETDEVICE',
'MPCCI_TINFO',
'MPCCI_SERVER',
'ABAQUS_CCI_DEBUG',
'NCPUS',
'OMP_DYNAMIC',
'OMP_NUM_THREADS',
'OUTDIR',
'PAIDUP',
'PARALLEL_METHOD',
'RAIDEV_NDREG_LAZYMEM',
'ABA_SYMBOLIC_GENERALCOLLAPSE',
'ABA_SYMBOLIC_GENERAL_MAXCLIQUERANK',
'ABA_ADM_MINIMUMINCREASE',
'ABA_ADM_MINIMUMDECREASE',
'IPATH_NO_CPUAFFINITY',
'MALLOC_MMAP_THRESHOLD_',
'ABA_EXT_SIMOUTPUT',
'SMA_WS',
'SMA_PARENT',
'SMA_PLATFORM',
'ABA_PRE_DECOMPOSITION',
'ACML_FAST_MALLOC',
'ACML_FAST_MALLOC_CHUNK_SIZE',
'ACML_FAST_MALLOC_MAX_CHUNKS',
'ACML_FAST_MALLOC_DEBUG')
import driverUtils, os
#-*- mode: python -*-
# #
# Compile and Link command settings for the Windows 64 Platform #
# ( AMD Opteron / Intel EM64T ) #
# #
compile_fortran=['ifort',
'/c','/DABQ_WIN86_64', '/u',
'/iface:cref', '/recursive', '/Qauto-scalar',
'/QxSSE3', '/QaxAVX',
'/heap-arrays:1',
# '/Od', '/Ob0' # <-- Optimization
# '/Zi', # <-- Debugging
'/include:%I', '/free', '/O1', '/fpp', '/openmp', '/Qmkl']
link_sl=['LINK',
'/nologo', '/NOENTRY', '/INCREMENTAL:NO', '/subsystem:console', '/machine:AMD64',
'/NODEFAULTLIB:LIBC.LIB', '/NODEFAULTLIB:LIBCMT.LIB',
'/DEFAULTLIB:OLDNAMES.LIB', '/DEFAULTLIB:LIBIFCOREMD.LIB', '/DEFAULTLIB:LIBIFPORTMD', '/DEFAULTLIB:LIBMMD.LIB',
'/DEFAULTLIB:kernel32.lib', '/DEFAULTLIB:user32.lib', '/DEFAULTLIB:advapi32.lib',
'/FIXED:NO', '/dll',
'/def:%E', '/out:%U', '%F', '%A', '%L', '%B',
'oldnames.lib', 'user32.lib', 'ws2_32.lib', 'netapi32.lib', 'advapi32.lib']
link_exe=['LINK',
'/nologo', '/INCREMENTAL:NO', '/subsystem:console', '/machine:AMD64', '/STACK:20000000',
'/NODEFAULTLIB:LIBC.LIB', '/NODEFAULTLIB:LIBCMT.LIB', '/DEFAULTLIB:OLDNAMES.LIB', '/DEFAULTLIB:LIBIFCOREMD.LIB',
'/DEFAULTLIB:LIBIFPORTMD', '/DEFAULTLIB:LIBMMD.LIB', '/DEFAULTLIB:kernel32.lib',
'/DEFAULTLIB:user32.lib', '/DEFAULTLIB:advapi32.lib',
'/FIXED:NO', '/LARGEADDRESSAWARE',
'/out:%J', '%F', '%M', '%L', '%B', '%O',
'oldnames.lib', 'user32.lib', 'ws2_32.lib', 'netapi32.lib', 'advapi32.lib']
# Link command to be used for MAKE w/o fortran compiler.
# remove the pound signs in order to remove the comments and have the file take effect.
#
#link_exe=['LINK', '/nologo', 'INCREMENTAL:NO', '/subsystem:console', '/machine:AMD64', '/NODEFAULTLIB:LIBC.LIB', '/NODEFAULTLIB:LIBCMT.LIB',
# '/DEFAULTLIB:OLDNAMES.LIB', '/DEFAULTLIB:MSVCRT.LIB', '/DEFAULTLIB:kernel32.lib', 'DEFAULTLIB:user32.lib', '/DEFAULTLIB:advapi32.lib',
# '/FIXED:NO', '/LARGEADDRESSAWARE', '/DEBUG', '/out:%J', '%F', '%M', '%L', '%B', '%O', 'oldnames.lib', 'user32.lib', 'ws2_32.lib',
# 'netapi32.lib', 'advapi32.lib]
# MPI Configuration
mp_mode = THREADS
mp_mpi_implementation = NATIVE
mp_rsh_command = 'dummy %H -l %U -n %C'
mp_mpirun_path = {}
mpirun = ''
progDir = os.environ.get('ProgramFiles','C:\\Program Files')
for mpiDir in ('Microsoft HPC Pack', 'Microsoft HPC Pack 2008 R2', 'Microsoft HPC Pack 2008', 'Microsoft HPC Pack 2008 SDK'):
mpirun = progDir + os.sep + mpiDir + os.sep + 'bin' + os.sep + 'mpiexec.exe'
if os.path.exists(mpirun):
mp_mpirun_path[NATIVE] = mpirun
mp_mpirun_path[MSSDK] = os.path.join(progDir, mpiDir)
break
if os.environ.has_key('CCP_HOME'):
from queueCCS import QueueCCS
queues['default'] = QueueCCS(queueName='share')
queues['share'] = QueueCCS(queueName='share')
queues['local'] = QueueCCS(queueName='local')
queues['genxmlshare'] = QueueCCS(queueName='genxmlshare')
queues['genxmllocal'] = QueueCCS(queueName='genxmllocal')
del QueueCCS
mpirun = os.path.join(os.environ['CCP_HOME'], 'bin', 'mpiexec.exe')
if os.path.exists(mpirun):
mp_mpirun_path[NATIVE] = mpirun
run_mode=BATCH
if mp_mpirun_path:
mp_mode=MPI
del progDir, mpiDir, mpirun
graphicsEnv = driverUtils.locateFile(os.environ['ABA_PATH'],'site','graphicsConfig','env')
if graphicsEnv:
execfile(graphicsEnv)
else:
raise 'Cannot find the graphics configuration environment file (graphicsConfig.env)'
del driverUtils, os, graphicsEnv
license_server_type=FLEXNET
abaquslm_license_file=""
doc_root="
doc_root_type="html"
academic=RESEARCH

2
installation/symlink_Code.py
View File

@ -1,4 +1,4 @@
#!/usr/bin/env python
#!/usr/bin/env python2
# -*- coding: UTF-8 no BOM -*-
import os,sys

6
installation/symlink_Processing.py
View File

@ -1,4 +1,4 @@
#!/usr/bin/env python
#!/usr/bin/env python2
# -*- coding: UTF-8 no BOM -*-
# Makes postprocessing routines acessible from everywhere.
@ -21,13 +21,13 @@ if not os.path.isdir(binDir):
#define ToDo list
processing_subDirs = ['pre','post','misc',]
processing_extensions = ['.py',]
processing_extensions = ['.py','.sh',]
for subDir in processing_subDirs:
theDir = os.path.abspath(os.path.join(baseDir,subDir))
for theFile in os.listdir(theDir):
if os.path.splitext(theFile)[1] in processing_extensions: # omit anything not fitting our script extensions (skip .py.bak, .py~, and the like)
if os.path.splitext(theFile)[1] in processing_extensions: # only consider files with proper extensions
src = os.path.abspath(os.path.join(theDir,theFile))
sym_link = os.path.abspath(os.path.join(binDir,os.path.splitext(theFile)[0]))

Some files were not shown because too many files have changed in this diff Show More