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

Merge branch 'HDF5-out-3' into 'development' 

Hdf5 out 3

See merge request damask/DAMASK!70
This commit is contained in:
Franz Roters 2019-04-10 20:39:37 +02:00
parent 21f308deec c2cb8e2033
commit da034f971c
23 changed files with 2474 additions and 2138 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
installation/mods_Abaqus/abaqus_v6.env
View File

@ -16,7 +16,7 @@ if False:
# use hdf5 compiler wrapper in $PATH
fortCmd = os.popen('h5fc -shlib -show').read().replace('\n','') # complicated way needed to pass in DAMASKVERSION string
link_sl += fortCmd.split()[1:]
fortCmd +=" -DDAMASKHDF5"
fortCmd +=" -DDAMASK_HDF5"
else:
# Use the version in $PATH
fortCmd = "ifort"

2
installation/mods_Abaqus/abaqus_v6_debug.env
View File

@ -16,7 +16,7 @@ if False:
# use hdf5 compiler wrapper in $PATH
fortCmd = os.popen('h5fc -shlib -show').read().replace('\n','') # complicated way needed to pass in DAMASKVERSION string
link_sl += fortCmd.split()[1:]
fortCmd +=" -DDAMASKHDF5"
fortCmd +=" -DDAMASK_HDF5"
else:
# Use the version in $PATH
fortCmd = "ifort"

2
installation/mods_MarcMentat/2018.1/Marc_tools/include_linux64
View File

@ -102,7 +102,7 @@ fi
if test "$DAMASK_HDF5" = "ON";then
H5FC="$(h5fc -shlib -show)"
HDF5_LIB=${H5FC//ifort/}
FCOMP="$H5FC -DDAMASKHDF5"
FCOMP="$H5FC -DDAMASK_HDF5"
echo $FCOMP
else
FCOMP=ifort

14
installation/mods_MarcMentat/2018/Marc_tools/include_linux64
View File

@ -63,7 +63,6 @@ else
INTEGER_PATH=/$MARC_INTEGER_SIZE
fi
FCOMP=ifort
INTELPATH="/opt/intel/compilers_and_libraries_2017/linux"
# find the root directory of the compiler installation:
@ -99,6 +98,16 @@ else
FCOMPROOT=
fi
# DAMASK uses the HDF5 compiler wrapper around the Intel compiler
if test "$DAMASK_HDF5" = "ON";then
H5FC="$(h5fc -shlib -show)"
HDF5_LIB=${H5FC//ifort/}
FCOMP="$H5FC -DDAMASK_HDF5"
echo $FCOMP
else
FCOMP=ifort
fi
# AEM
if test "$MARCDLLOUTDIR" = ""; then
DLLOUTDIR="$MARC_LIB"
@ -535,6 +544,7 @@ else
DAMASKVERSION="'N/A'"
fi
# DAMASK compiler calls: additional flags are in line 2 OpenMP flags in line 3
DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -O0 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018 -DDAMASKVERSION=$DAMASKVERSION \
@ -737,7 +747,7 @@ SECLIBS="-L$MARC_LIB -llapi"
SOLVERLIBS="${BCSSOLVERLIBS} ${VKISOLVERLIBS} ${CASISOLVERLIBS} ${MF2SOLVERLIBS} \
$MKLLIB -L$MARC_MKL -liomp5 \
$MARC_LIB/blas_src.a ${ACSI_LIB}/ACSI_MarcLib.a $KDTREE2_LIB/kdtree2.a "
$MARC_LIB/blas_src.a ${ACSI_LIB}/ACSI_MarcLib.a $KDTREE2_LIB/kdtree2.a $HDF5_LIB "
SOLVERLIBS_DLL=${SOLVERLIBS}
if test "$AEM_DLL" -eq 1

10
src/CPFEM.f90
View File

@ -72,6 +72,12 @@ subroutine CPFEM_initAll(el,ip)
mesh_init
use material, only: &
material_init
#ifdef DAMASK_HDF5
use HDF5_utilities, only: &
HDF5_utilities_init
use results, only: &
results_init
#endif
use lattice, only: &
lattice_init
use constitutive, only: &
@ -100,6 +106,10 @@ subroutine CPFEM_initAll(el,ip)
call FE_init
call mesh_init(ip, el)
call lattice_init
#ifdef DAMASK_HDF5
call HDF5_utilities_init
call results_init
#endif
call material_init
call constitutive_init
call crystallite_init

7
src/CPFEM2.f90
View File

@ -72,9 +72,9 @@ subroutine CPFEM_initAll()
call FE_init
call mesh_init
call lattice_init
call material_init
call HDF5_utilities_init
call results_init
call material_init
call constitutive_init
call crystallite_init
call homogenization_init
@ -300,6 +300,8 @@ subroutine CPFEM_results(inc,time)
use HDF5_utilities
use constitutive, only: &
constitutive_results
use crystallite, only: &
crystallite_results
implicit none
integer(pInt), intent(in) :: inc
@ -307,7 +309,8 @@ subroutine CPFEM_results(inc,time)
call results_openJobFile
call results_addIncrement(inc,time)
call constitutive_results()
call constitutive_results
call crystallite_results
call results_removeLink('current') ! ToDo: put this into closeJobFile
call results_closeJobFile

2289
src/HDF5_utilities.f90
View File

File diff suppressed because it is too large Load Diff

8
src/commercialFEM_fileList.f90
View File

@ -9,10 +9,6 @@
#include "list.f90"
#include "future.f90"
#include "config.f90"
#ifdef DAMASKHDF5
#include "HDF5_utilities.f90"
#include "results.f90"
#endif
#include "math.f90"
#include "quaternions.f90"
#include "Lambert.f90"
@ -26,6 +22,10 @@
#ifdef Marc4DAMASK
#include "mesh_marc.f90"
#endif
#ifdef DAMASK_HDF5
#include "HDF5_utilities.f90"
#include "results.f90"
#endif
#include "material.f90"
#include "lattice.f90"
#include "source_thermal_dissipation.f90"

15
src/config.f90
View File

@ -231,15 +231,21 @@ end function read_materialConfig
!--------------------------------------------------------------------------------------------------
subroutine parse_materialConfig(sectionNames,part,line, &
fileContent)
use prec, only: &
pStringLen
use IO, only: &
IO_intOut
implicit none
character(len=64), allocatable, dimension(:), intent(out) :: sectionNames
type(tPartitionedStringList), allocatable, dimension(:), intent(inout) :: part
character(len=pStringLen), intent(inout) :: line
character(len=pStringLen), dimension(:), intent(in) :: fileContent
integer, allocatable, dimension(:) :: partPosition ! position of [] tags + last line in section
integer :: i, j
logical :: echo
integer, allocatable, dimension(:) :: partPosition !< position of [] tags + last line in section
integer :: i, j
logical :: echo
character(len=pStringLen) :: section_ID
echo = .false.
@ -263,7 +269,8 @@ subroutine parse_materialConfig(sectionNames,part,line, &
partPosition = [partPosition, i] ! needed when actually storing content
do i = 1, size(partPosition) -1
sectionNames(i) = trim(adjustl(IO_getTag(fileContent(partPosition(i)),'[',']')))
write(section_ID,'('//IO_intOut(size(partPosition))//')') i
sectionNames(i) = trim(section_ID)//'_'//trim(adjustl(IO_getTag(fileContent(partPosition(i)),'[',']')))
do j = partPosition(i) + 1, partPosition(i+1) -1
call part(i)%add(trim(adjustl(fileContent(j))))
enddo

149
src/constitutive.f90
View File

@ -9,7 +9,7 @@ module constitutive
implicit none
private
integer(pInt), public, protected :: &
integer, public, protected :: &
constitutive_plasticity_maxSizePostResults, &
constitutive_plasticity_maxSizeDotState, &
constitutive_source_maxSizePostResults, &
@ -37,7 +37,7 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates arrays pointing to array of the various constitutive modules
!--------------------------------------------------------------------------------------------------
subroutine constitutive_init()
subroutine constitutive_init
use prec, only: &
pReal
use debug, only: &
@ -50,8 +50,7 @@ subroutine constitutive_init()
IO_write_jobFile
use config, only: &
material_Nphase, &
phase_name, &
config_deallocate
phase_name
use material, only: &
material_phase, &
phase_plasticity, &
@ -111,14 +110,14 @@ subroutine constitutive_init()
use kinematics_thermal_expansion
implicit none
integer(pInt), parameter :: FILEUNIT = 204_pInt
integer(pInt) :: &
integer, parameter :: FILEUNIT = 204
integer :: &
o, & !< counter in output loop
ph, & !< counter in phase loop
s, & !< counter in source loop
ins !< instance of plasticity/source
integer(pInt), dimension(:,:), pointer :: thisSize
integer, dimension(:,:), pointer :: thisSize
character(len=64), dimension(:,:), pointer :: thisOutput
character(len=32) :: outputName !< name of output, intermediate fix until HDF5 output is ready
logical :: knownPlasticity, knownSource, nonlocalConstitutionPresent
@ -149,15 +148,13 @@ subroutine constitutive_init()
if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init
if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init
call config_deallocate('material.config/phase')
write(6,'(/,a)') ' <<<+- constitutive init -+>>>'
mainProcess: if (worldrank == 0) then
!--------------------------------------------------------------------------------------------------
! write description file for constitutive output
call IO_write_jobFile(FILEUNIT,'outputConstitutive')
PhaseLoop: do ph = 1_pInt,material_Nphase
PhaseLoop: do ph = 1,material_Nphase
activePhase: if (any(material_phase == ph)) then
ins = phase_plasticityInstance(ph)
knownPlasticity = .true. ! assume valid
@ -197,14 +194,14 @@ subroutine constitutive_init()
if (knownPlasticity) then
write(FILEUNIT,'(a)') '(plasticity)'//char(9)//trim(outputName)
if (phase_plasticity(ph) /= PLASTICITY_NONE_ID) then
OutputPlasticityLoop: do o = 1_pInt,size(thisOutput(:,ins))
if(len(trim(thisOutput(o,ins))) > 0_pInt) &
OutputPlasticityLoop: do o = 1,size(thisOutput(:,ins))
if(len(trim(thisOutput(o,ins))) > 0) &
write(FILEUNIT,'(a,i4)') trim(thisOutput(o,ins))//char(9),thisSize(o,ins)
enddo OutputPlasticityLoop
endif
endif
SourceLoop: do s = 1_pInt, phase_Nsources(ph)
SourceLoop: do s = 1, phase_Nsources(ph)
knownSource = .true. ! assume valid
sourceType: select case (phase_source(s,ph))
case (SOURCE_thermal_dissipation_ID) sourceType
@ -242,8 +239,8 @@ subroutine constitutive_init()
end select sourceType
if (knownSource) then
write(FILEUNIT,'(a)') '(source)'//char(9)//trim(outputName)
OutputSourceLoop: do o = 1_pInt,size(thisOutput(:,ins))
if(len(trim(thisOutput(o,ins))) > 0_pInt) &
OutputSourceLoop: do o = 1,size(thisOutput(:,ins))
if(len(trim(thisOutput(o,ins))) > 0) &
write(FILEUNIT,'(a,i4)') trim(thisOutput(o,ins))//char(9),thisSize(o,ins)
enddo OutputSourceLoop
endif
@ -253,17 +250,17 @@ subroutine constitutive_init()
close(FILEUNIT)
endif mainProcess
constitutive_plasticity_maxSizeDotState = 0_pInt
constitutive_plasticity_maxSizePostResults = 0_pInt
constitutive_source_maxSizeDotState = 0_pInt
constitutive_source_maxSizePostResults = 0_pInt
constitutive_plasticity_maxSizeDotState = 0
constitutive_plasticity_maxSizePostResults = 0
constitutive_source_maxSizeDotState = 0
constitutive_source_maxSizePostResults = 0
PhaseLoop2:do ph = 1_pInt,material_Nphase
PhaseLoop2:do ph = 1,material_Nphase
!--------------------------------------------------------------------------------------------------
! partition and inititalize state
plasticState(ph)%partionedState0 = plasticState(ph)%state0
plasticState(ph)%state = plasticState(ph)%partionedState0
forall(s = 1_pInt:phase_Nsources(ph))
forall(s = 1:phase_Nsources(ph))
sourceState(ph)%p(s)%partionedState0 = sourceState(ph)%p(s)%state0
sourceState(ph)%p(s)%state = sourceState(ph)%p(s)%partionedState0
end forall
@ -302,7 +299,7 @@ function constitutive_homogenizedC(ipc,ip,el)
implicit none
real(pReal), dimension(6,6) :: constitutive_homogenizedC
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
@ -341,14 +338,14 @@ subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el)
plastic_disloUCLA_dependentState
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: &
Fe, & !< elastic deformation gradient
Fp !< plastic deformation gradient
integer(pInt) :: &
integer :: &
ho, & !< homogenization
tme, & !< thermal member position
instance, of
@ -412,7 +409,7 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
plastic_nonlocal_LpAndItsTangent
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
@ -428,10 +425,10 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
dLp_dMp !< derivative of Lp with respect to Mandel stress
real(pReal), dimension(3,3) :: &
Mp !< Mandel stress work conjugate with Lp
integer(pInt) :: &
integer :: &
ho, & !< homogenization
tme !< thermal member position
integer(pInt) :: &
integer :: &
i, j, instance, of
ho = material_homogenizationAt(el)
@ -519,7 +516,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
kinematics_thermal_expansion_LiAndItsTangent
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
@ -541,7 +538,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
my_dLi_dS
real(pReal) :: &
detFi
integer(pInt) :: &
integer :: &
k, i, j, &
instance, of
@ -562,7 +559,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS
KinematicsLoop: do k = 1_pInt, phase_Nkinematics(material_phase(ipc,ip,el))
KinematicsLoop: do k = 1, phase_Nkinematics(material_phase(ipc,ip,el))
kinematicsType: select case (phase_kinematics(k,material_phase(ipc,ip,el)))
case (KINEMATICS_cleavage_opening_ID) kinematicsType
call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ipc, ip, el)
@ -583,7 +580,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration
temp_33 = matmul(FiInv,Li)
do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt
do i = 1,3; do j = 1,3
dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi
dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i)
dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i)
@ -612,22 +609,22 @@ pure function constitutive_initialFi(ipc, ip, el)
kinematics_thermal_expansion_initialStrain
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(3,3) :: &
constitutive_initialFi !< composite initial intermediate deformation gradient
integer(pInt) :: &
integer :: &
k !< counter in kinematics loop
integer(pInt) :: &
integer :: &
phase, &
homog, offset
constitutive_initialFi = math_I3
phase = material_phase(ipc,ip,el)
KinematicsLoop: do k = 1_pInt, phase_Nkinematics(phase) !< Warning: small initial strain assumption
KinematicsLoop: do k = 1, phase_Nkinematics(phase) !< Warning: small initial strain assumption
kinematicsType: select case (phase_kinematics(k,phase))
case (KINEMATICS_thermal_expansion_ID) kinematicsType
homog = material_homogenizationAt(el)
@ -650,7 +647,7 @@ subroutine constitutive_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el)
pReal
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
@ -691,7 +688,7 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
STIFFNESS_DEGRADATION_damage_ID
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
@ -705,19 +702,19 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
real(pReal), dimension(3,3) :: E
real(pReal), dimension(3,3,3,3) :: C
integer(pInt) :: &
integer :: &
ho, & !< homogenization
d !< counter in degradation loop
integer(pInt) :: &
integer :: &
i, j
ho = material_homogenizationAt(el)
C = math_66toSym3333(constitutive_homogenizedC(ipc,ip,el))
DegradationLoop: do d = 1_pInt, phase_NstiffnessDegradations(material_phase(ipc,ip,el))
DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phase(ipc,ip,el))
degradationType: select case(phase_stiffnessDegradation(d,material_phase(ipc,ip,el)))
case (STIFFNESS_DEGRADATION_damage_ID) degradationType
C = C * damage(ho)%p(damageMapping(ho)%p(ip,el))**2_pInt
C = C * damage(ho)%p(damageMapping(ho)%p(ip,el))**2
end select degradationType
enddo DegradationLoop
@ -725,7 +722,7 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
dS_dFe = 0.0_pReal
forall (i=1_pInt:3_pInt, j=1_pInt:3_pInt)
forall (i=1:3, j=1:3)
dS_dFe(i,j,1:3,1:3) = &
matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn
@ -790,7 +787,7 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip,
source_thermal_externalheat_dotState
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
@ -805,7 +802,7 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip,
S !< 2nd Piola Kirchhoff stress (vector notation)
real(pReal), dimension(3,3) :: &
Mp
integer(pInt) :: &
integer :: &
ho, & !< homogenization
tme, & !< thermal member position
i, & !< counter in source loop
@ -848,7 +845,7 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip,
subdt,ip,el)
end select plasticityType
SourceLoop: do i = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
SourceLoop: do i = 1, phase_Nsources(material_phase(ipc,ip,el))
sourceType: select case (phase_source(i,material_phase(ipc,ip,el)))
@ -900,7 +897,7 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
source_damage_isoBrittle_deltaState
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
@ -910,7 +907,7 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
Fi !< intermediate deformation gradient
real(pReal), dimension(3,3) :: &
Mp
integer(pInt) :: &
integer :: &
i, &
instance, of
@ -928,7 +925,7 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
end select plasticityType
sourceLoop: do i = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
sourceLoop: do i = 1, phase_Nsources(material_phase(ipc,ip,el))
sourceType: select case (phase_source(i,material_phase(ipc,ip,el)))
@ -994,7 +991,7 @@ function constitutive_postResults(S, Fi, ipc, ip, el)
source_damage_anisoDuctile_postResults
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
@ -1007,9 +1004,9 @@ function constitutive_postResults(S, Fi, ipc, ip, el)
S !< 2nd Piola Kirchhoff stress
real(pReal), dimension(3,3) :: &
Mp !< Mandel stress
integer(pInt) :: &
integer :: &
startPos, endPos
integer(pInt) :: &
integer :: &
ho, & !< homogenization
tme, & !< thermal member position
i, of, instance !< counter in source loop
@ -1021,7 +1018,7 @@ function constitutive_postResults(S, Fi, ipc, ip, el)
ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el)
startPos = 1_pInt
startPos = 1
endPos = plasticState(material_phase(ipc,ip,el))%sizePostResults
of = phasememberAt(ipc,ip,el)
@ -1054,8 +1051,8 @@ function constitutive_postResults(S, Fi, ipc, ip, el)
end select plasticityType
SourceLoop: do i = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
startPos = endPos + 1_pInt
SourceLoop: do i = 1, phase_Nsources(material_phase(ipc,ip,el))
startPos = endPos + 1
endPos = endPos + sourceState(material_phase(ipc,ip,el))%p(i)%sizePostResults
of = phasememberAt(ipc,ip,el)
sourceType: select case (phase_source(i,material_phase(ipc,ip,el)))
@ -1077,7 +1074,7 @@ end function constitutive_postResults
!--------------------------------------------------------------------------------------------------
!> @brief writes constitutive results to HDF5 output file
!--------------------------------------------------------------------------------------------------
subroutine constitutive_results()
subroutine constitutive_results
use material, only: &
PLASTICITY_ISOTROPIC_ID, &
PLASTICITY_PHENOPOWERLAW_ID, &
@ -1085,7 +1082,7 @@ subroutine constitutive_results()
PLASTICITY_DISLOTWIN_ID, &
PLASTICITY_DISLOUCLA_ID, &
PLASTICITY_NONLOCAL_ID
#if defined(PETSc) || defined(DAMASKHDF5)
#if defined(PETSc) || defined(DAMASK_HDF5)
use results
use HDF5_utilities
use config, only: &
@ -1108,36 +1105,40 @@ subroutine constitutive_results()
use plastic_nonlocal, only: &
plastic_nonlocal_results
implicit none
integer :: p
call HDF5_closeGroup(results_addGroup('current/phase'))
do p=1,size(config_name_phase)
call HDF5_closeGroup(results_addGroup('current/phase/'//trim(config_name_phase(p))))
integer :: p
character(len=256) :: group
call HDF5_closeGroup(results_addGroup('current/constitutive'))
do p=1,size(config_name_phase)
group = trim('current/constitutive')//'/'//trim(config_name_phase(p))
call HDF5_closeGroup(results_addGroup(group))
group = trim(group)//'/'//'plastic'
call HDF5_closeGroup(results_addGroup(group))
select case(material_phase_plasticity_type(p))
case(PLASTICITY_ISOTROPIC_ID)
call plastic_isotropic_results(phase_plasticityInstance(p),'current/phase/'//trim(config_name_phase(p)))
call plastic_isotropic_results(phase_plasticityInstance(p),group)
case(PLASTICITY_PHENOPOWERLAW_ID)
call plastic_phenopowerlaw_results(phase_plasticityInstance(p),'current/phase/'//trim(config_name_phase(p)))
call plastic_phenopowerlaw_results(phase_plasticityInstance(p),group)
case(PLASTICITY_KINEHARDENING_ID)
call plastic_kinehardening_results(phase_plasticityInstance(p),'current/phase/'//trim(config_name_phase(p)))
case(PLASTICITY_KINEHARDENING_ID)
call plastic_kinehardening_results(phase_plasticityInstance(p),group)
case(PLASTICITY_DISLOTWIN_ID)
call plastic_dislotwin_results(phase_plasticityInstance(p),'current/phase/'//trim(config_name_phase(p)))
case(PLASTICITY_DISLOTWIN_ID)
call plastic_dislotwin_results(phase_plasticityInstance(p),group)
case(PLASTICITY_DISLOUCLA_ID)
call plastic_disloUCLA_results(phase_plasticityInstance(p),'current/phase/'//trim(config_name_phase(p)))
case(PLASTICITY_DISLOUCLA_ID)
call plastic_disloUCLA_results(phase_plasticityInstance(p),group)
case(PLASTICITY_NONLOCAL_ID)
call plastic_nonlocal_results(phase_plasticityInstance(p),'current/phase/'//trim(config_name_phase(p)))
case(PLASTICITY_NONLOCAL_ID)
call plastic_nonlocal_results(phase_plasticityInstance(p),group)
end select
enddo
enddo
#endif

179
src/crystallite.f90
View File

@ -10,7 +10,8 @@
module crystallite
use prec, only: &
pReal
pReal, &
pStringLen
use rotations, only: &
rotation
use FEsolving, only: &
@ -103,6 +104,13 @@ module crystallite
end enum
integer(kind(undefined_ID)),dimension(:,:), allocatable, private :: &
crystallite_outputID !< ID of each post result output
type, private :: tOutput !< new requested output (per phase)
character(len=65536), allocatable, dimension(:) :: &
label
end type tOutput
type(tOutput), allocatable, dimension(:), private :: output_constituent
procedure(), pointer :: integrateState
public :: &
@ -111,7 +119,8 @@ module crystallite
crystallite_stressTangent, &
crystallite_orientations, &
crystallite_push33ToRef, &
crystallite_postResults
crystallite_postResults, &
crystallite_results
private :: &
integrateStress, &
integrateState, &
@ -156,6 +165,7 @@ subroutine crystallite_init
use config, only: &
config_deallocate, &
config_crystallite, &
config_phase, &
crystallite_name
use constitutive, only: &
constitutive_initialFi, &
@ -296,6 +306,18 @@ subroutine crystallite_init
end select outputName
enddo
enddo
allocate(output_constituent(size(config_phase)))
do c = 1, size(config_phase)
#if defined(__GFORTRAN__)
allocate(output_constituent(c)%label(1))
output_constituent(c)%label(1)= 'GfortranBug86277'
output_constituent(c)%label = config_phase(c)%getStrings('(output)',defaultVal=output_constituent(c)%label )
if (output_constituent(c)%label (1) == 'GfortranBug86277') output_constituent(c)%label = [character(len=pStringLen)::]
#else
output_constituent(c)%label = config_phase(c)%getStrings('(output)',defaultVal=[character(len=pStringLen)::])
#endif
enddo
do r = 1,size(config_crystallite)
@ -340,6 +362,7 @@ subroutine crystallite_init
close(FILEUNIT)
endif
call config_deallocate('material.config/phase')
call config_deallocate('material.config/crystallite')
!--------------------------------------------------------------------------------------------------
@ -1053,6 +1076,158 @@ function crystallite_postResults(ipc, ip, el)
end function crystallite_postResults
!--------------------------------------------------------------------------------------------------
!> @brief writes constitutive results to HDF5 output file
!--------------------------------------------------------------------------------------------------
subroutine crystallite_results
#if defined(PETSc) || defined(DAMASK_HDF5)
use lattice
use results
use HDF5_utilities
use rotations
use config, only: &
config_name_phase => phase_name ! anticipate logical name
use material, only: &
material_phase_plasticity_type => phase_plasticity
implicit none
integer :: p,o
real(pReal), allocatable, dimension(:,:,:) :: selected_tensors
type(rotation), allocatable, dimension(:) :: selected_rotations
character(len=256) :: group,lattice_label
call HDF5_closeGroup(results_addGroup('current/constituent'))
do p=1,size(config_name_phase)
group = trim('current/constituent')//'/'//trim(config_name_phase(p))
call HDF5_closeGroup(results_addGroup(group))
do o = 1, size(output_constituent(p)%label)
select case (output_constituent(p)%label(o))
case('f')
selected_tensors = select_tensors(crystallite_partionedF,p)
call results_writeDataset(group,selected_tensors,'F',&
'deformation gradient','1')
case('fe')
selected_tensors = select_tensors(crystallite_Fe,p)
call results_writeDataset(group,selected_tensors,'Fe',&
'elastic deformation gradient','1')
case('fp')
selected_tensors = select_tensors(crystallite_Fp,p)
call results_writeDataset(group,selected_tensors,'Fp',&
'plastic deformation gradient','1')
case('fi')
selected_tensors = select_tensors(crystallite_Fi,p)
call results_writeDataset(group,selected_tensors,'Fi',&
'inelastic deformation gradient','1')
case('lp')
selected_tensors = select_tensors(crystallite_Lp,p)
call results_writeDataset(group,selected_tensors,'Lp',&
'plastic velocity gradient','1/s')
case('li')
selected_tensors = select_tensors(crystallite_Li,p)
call results_writeDataset(group,selected_tensors,'Li',&
'inelastic velocity gradient','1/s')
case('p')
selected_tensors = select_tensors(crystallite_P,p)
call results_writeDataset(group,selected_tensors,'P',&
'1st Piola-Kirchoff stress','Pa')
case('s')
selected_tensors = select_tensors(crystallite_S,p)
call results_writeDataset(group,selected_tensors,'S',&
'2nd Piola-Kirchoff stress','Pa')
case('orientation')
select case(lattice_structure(p))
case(LATTICE_iso_ID)
lattice_label = 'iso'
case(LATTICE_fcc_ID)
lattice_label = 'fcc'
case(LATTICE_bcc_ID)
lattice_label = 'bcc'
case(LATTICE_bct_ID)
lattice_label = 'bct'
case(LATTICE_hex_ID)
lattice_label = 'hex'
case(LATTICE_ort_ID)
lattice_label = 'ort'
end select
selected_rotations = select_rotations(crystallite_orientation,p)
call results_writeDataset(group,selected_rotations,'orientation',&
'crystal orientation as quaternion',lattice_label)
end select
enddo
enddo
contains
!--------------------------------------------------------------------------------------------------
!> @brief select tensors for output
!--------------------------------------------------------------------------------------------------
function select_tensors(dataset,instance)
use material, only: &
homogenization_maxNgrains, &
material_phaseAt
integer, intent(in) :: instance
real(pReal), dimension(:,:,:,:,:), intent(in) :: dataset
real(pReal), allocatable, dimension(:,:,:) :: select_tensors
integer :: e,i,c,j
allocate(select_tensors(3,3,count(material_phaseAt==instance)*homogenization_maxNgrains))
j=1
do e = 1, size(material_phaseAt,2)
do i = 1, homogenization_maxNgrains !ToDo: this needs to be changed for varying Ngrains
do c = 1, size(material_phaseAt,1)
if (material_phaseAt(c,e) == instance) then
select_tensors(1:3,1:3,j) = dataset(1:3,1:3,c,i,e)
j = j + 1
endif
enddo
enddo
enddo
end function select_tensors
!--------------------------------------------------------------------------------------------------
!> @brief select rotations for output
!--------------------------------------------------------------------------------------------------
function select_rotations(dataset,instance)
use material, only: &
homogenization_maxNgrains, &
material_phaseAt
integer, intent(in) :: instance
type(rotation), dimension(:,:,:), intent(in) :: dataset
type(rotation), allocatable, dimension(:) :: select_rotations
integer :: e,i,c,j
allocate(select_rotations(count(material_phaseAt==instance)*homogenization_maxNgrains))
j=1
do e = 1, size(material_phaseAt,2)
do i = 1, homogenization_maxNgrains !ToDo: this needs to be changed for varying Ngrains
do c = 1, size(material_phaseAt,1)
if (material_phaseAt(c,e) == instance) then
select_rotations(j) = dataset(c,i,e)
j = j + 1
endif
enddo
enddo
enddo
end function select_rotations
#endif
end subroutine crystallite_results
!--------------------------------------------------------------------------------------------------
!> @brief calculation of stress (P) with time integration based on a residuum in Lp and
!> intermediate acceleration of the Newton-Raphson correction

10
src/lattice.f90
View File

@ -507,10 +507,12 @@ module lattice
public :: &
lattice_init, &
lattice_qDisorientation, &
LATTICE_iso_ID, &
LATTICE_fcc_ID, &
LATTICE_bcc_ID, &
LATTICE_bct_ID, &
LATTICE_hex_ID, &
LATTICE_ort_ID, &
lattice_SchmidMatrix_slip, &
lattice_SchmidMatrix_twin, &
lattice_SchmidMatrix_trans, &
@ -581,18 +583,18 @@ subroutine lattice_init
do p = 1, size(config_phase)
tag = config_phase(p)%getString('lattice_structure')
select case(trim(tag))
case('iso','isotropic')
select case(trim(tag(1:3)))
case('iso')
lattice_structure(p) = LATTICE_iso_ID
case('fcc')
lattice_structure(p) = LATTICE_fcc_ID
case('bcc')
lattice_structure(p) = LATTICE_bcc_ID
case('hex','hexagonal')
case('hex')
lattice_structure(p) = LATTICE_hex_ID
case('bct')
lattice_structure(p) = LATTICE_bct_ID
case('ort','orthorhombic')
case('ort')
lattice_structure(p) = LATTICE_ort_ID
end select

73
src/material.f90
View File

@ -147,16 +147,14 @@ module material
damage_initialPhi !< initial damage per each homogenization
! NEW MAPPINGS
integer(pInt), dimension(:), allocatable, public, protected :: &
material_homogenizationAt, & !< homogenization ID of each element (copy of mesh_homogenizationAt)
material_homogenizationMemberAt, & !< position of the element within its homogenization instance
material_aggregateAt, & !< aggregate ID of each element FUTURE USE FOR OUTPUT
material_aggregatMemberAt !< position of the element within its aggregate instance FUTURE USE FOR OUTPUT
integer(pInt), dimension(:,:), allocatable, public, protected :: &
material_phaseAt, & !< phase ID of each element
material_phaseMemberAt, & !< position of the element within its phase instance
material_crystalliteAt, & !< crystallite ID of each element CURRENTLY NOT PER CONSTITUTENT
material_crystalliteMemberAt !< position of the element within its crystallite instance CURRENTLY NOT PER CONSTITUTENT
integer, dimension(:), allocatable, public, protected :: & ! (elem)
material_homogenizationAt !< homogenization ID of each element (copy of mesh_homogenizationAt)
integer, dimension(:,:), allocatable, public, protected :: & ! (ip,elem)
material_homogenizationMemberAt !< position of the element within its homogenization instance
integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem)
material_phaseAt !< phase ID of each element
integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,ip,elem)
material_phaseMemberAt !< position of the element within its phase instance
! END NEW MAPPINGS
! DEPRECATED: use material_phaseAt
@ -275,7 +273,10 @@ contains
!> @details figures out if solverJobName.materialConfig is present, if not looks for
!> material.config
!--------------------------------------------------------------------------------------------------
subroutine material_init()
subroutine material_init
#if defined(PETSc) || defined(DAMASK_HDF5)
use results
#endif
use IO, only: &
IO_error
use debug, only: &
@ -382,21 +383,57 @@ subroutine material_init()
endif debugOut
call material_populateGrains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! new mappings
allocate(material_homogenizationAt,source=theMesh%homogenizationAt)
allocate(material_homogenizationMemberAt(theMesh%elem%nIPs,theMesh%Nelems),source=0)
allocate(CounterHomogenization(size(config_homogenization)),source=0)
do e = 1, theMesh%Nelems
do i = 1, theMesh%elem%nIPs
CounterHomogenization(material_homogenizationAt(e)) = &
CounterHomogenization(material_homogenizationAt(e)) + 1
material_homogenizationMemberAt(i,e) = CounterHomogenization(material_homogenizationAt(e))
enddo
enddo
allocate(material_phaseAt(homogenization_maxNgrains,theMesh%Nelems), source=material_phase(:,1,:))
allocate(material_phaseMemberAt(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0)
allocate(CounterPhase(size(config_phase)),source=0)
do e = 1, theMesh%Nelems
do i = 1, theMesh%elem%nIPs
do c = 1, homogenization_maxNgrains
CounterPhase(material_phaseAt(c,e)) = &
CounterPhase(material_phaseAt(c,e)) + 1
material_phaseMemberAt(c,i,e) = CounterPhase(material_phaseAt(c,e))
enddo
enddo
enddo
#if defined(PETSc) || defined(DAMASK_HDF5)
call results_openJobFile
call results_mapping_constituent(material_phaseAt,material_phaseMemberAt,phase_name)
call results_mapping_materialpoint(material_homogenizationAt,material_homogenizationMemberAt,homogenization_name)
call results_closeJobFile
#endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BEGIN DEPRECATED
allocate(phaseAt ( homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0_pInt)
allocate(phasememberAt ( homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0_pInt)
allocate(mappingHomogenization (2, theMesh%elem%nIPs,theMesh%Nelems),source=0_pInt)
allocate(mappingHomogenizationConst( theMesh%elem%nIPs,theMesh%Nelems),source=1_pInt)
! END DEPRECATED
allocate(material_homogenizationAt,source=theMesh%homogenizationAt)
allocate(material_AggregateAt, source=theMesh%homogenizationAt)
allocate(CounterPhase (size(config_phase)), source=0_pInt)
allocate(CounterHomogenization(size(config_homogenization)),source=0_pInt)
CounterHomogenization=0
CounterPhase =0
! BEGIN DEPRECATED
do e = 1_pInt,theMesh%Nelems
myHomog = theMesh%homogenizationAt(e)
do i = 1_pInt, theMesh%elem%nIPs

2
src/numerics.f90
View File

@ -21,7 +21,7 @@ module numerics
pert_method = 1_pInt, & !< method used in perturbation technique for tangent
randomSeed = 0_pInt, & !< fixed seeding for pseudo-random number generator, Default 0: use random seed
worldrank = 0_pInt, & !< MPI worldrank (/=0 for MPI simulations only)
worldsize = 0_pInt, & !< MPI worldsize (/=0 for MPI simulations only)
worldsize = 1_pInt, & !< MPI worldsize (/=1 for MPI simulations only)
numerics_integrator = 1_pInt !< method used for state integration Default 1: fix-point iteration
integer(4), protected, public :: &
DAMASK_NumThreadsInt = 0 !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive

32
src/plastic_disloUCLA.f90
View File

@ -560,23 +560,41 @@ end function plastic_disloUCLA_postResults
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
subroutine plastic_disloUCLA_results(instance,group)
#if defined(PETSc) || defined(DAMASKHDF5)
use results
#if defined(PETSc) || defined(DAMASK_HDF5)
use results, only: &
results_writeDataset
implicit none
integer, intent(in) :: instance
character(len=*) :: group
integer, intent(in) :: instance
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(instance), stt => state(instance))
associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (rho_mob_ID)
call results_writeDataset(group,stt%rho_mob,'rho_mob',&
'mobile dislocation density','1/m²')
case (rho_dip_ID)
call results_writeDataset(group,stt%rho_dip,'rho_dip',&
'dislocation dipole density''1/m²')
case (dot_gamma_sl_ID)
call results_writeDataset(group,stt%gamma_sl,'dot_gamma_sl',&
'plastic shear','1')
case (Lambda_sl_ID)
call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',&
'mean free path for slip','m')
case (thresholdstress_ID)
call results_writeDataset(group,dst%threshold_stress,'threshold_stress',&
'threshold stress for slip','Pa')
end select
enddo outputsLoop
end associate
#else
integer, intent(in) :: instance
character(len=*) :: group
integer, intent(in) :: instance
character(len=*), intent(in) :: group
#endif
end subroutine plastic_disloUCLA_results

57
src/plastic_dislotwin.f90
View File

@ -148,7 +148,7 @@ module plastic_dislotwin
type(tDislotwinState), allocatable, dimension(:), private :: &
dotState, &
state
type(tDislotwinMicrostructure), allocatable, dimension(:), private :: microstructure
type(tDislotwinMicrostructure), allocatable, dimension(:), private :: dependentState
public :: &
plastic_dislotwin_init, &
@ -241,14 +241,14 @@ subroutine plastic_dislotwin_init
allocate(param(Ninstance))
allocate(state(Ninstance))
allocate(dotState(Ninstance))
allocate(microstructure(Ninstance))
allocate(dependentState(Ninstance))
do p = 1, size(phase_plasticity)
if (phase_plasticity(p) /= PLASTICITY_DISLOTWIN_ID) cycle
associate(prm => param(phase_plasticityInstance(p)), &
dot => dotState(phase_plasticityInstance(p)), &
stt => state(phase_plasticityInstance(p)), &
dst => microstructure(phase_plasticityInstance(p)), &
dst => dependentState(phase_plasticityInstance(p)), &
config => config_phase(p))
prm%aTol_rho = config%getFloat('atol_rho', defaultVal=0.0_pReal)
@ -801,7 +801,7 @@ subroutine plastic_dislotwin_dotState(Mp,T,instance,of)
dot_gamma_tr
associate(prm => param(instance), stt => state(instance), &
dot => dotstate(instance), dst => microstructure(instance))
dot => dotstate(instance), dst => dependentState(instance))
f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,of)) &
@ -897,7 +897,7 @@ subroutine plastic_dislotwin_dependentState(T,instance,of)
associate(prm => param(instance),&
stt => state(instance),&
dst => microstructure(instance))
dst => dependentState(instance))
sumf_twin = sum(stt%f_tw(1:prm%sum_N_tw,of))
sumf_trans = sum(stt%f_tr(1:prm%sum_N_tr,of))
@ -1002,7 +1002,7 @@ function plastic_dislotwin_postResults(Mp,T,instance,of) result(postResults)
integer :: &
o,c,j
associate(prm => param(instance), stt => state(instance), dst => microstructure(instance))
associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
c = 0
@ -1063,20 +1063,53 @@ end function plastic_dislotwin_postResults
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
subroutine plastic_dislotwin_results(instance,group)
#if defined(PETSc) || defined(DAMASKHDF5)
use results
#if defined(PETSc) || defined(DAMASK_HDF5)
use results, only: &
results_writeDataset
implicit none
integer, intent(in) :: instance
character(len=*) :: group
integer :: o
associate(prm => param(instance), stt => state(instance))
associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (rho_mob_ID)
call results_writeDataset(group,stt%rho_mob,'rho_mob',&
'mobile dislocation density','1/m²')
case (rho_dip_ID)
call results_writeDataset(group,stt%rho_dip,'rho_dip',&
'dislocation dipole density''1/m²')
case (dot_gamma_sl_ID)
call results_writeDataset(group,stt%gamma_sl,'dot_gamma_sl',&
'plastic shear','1')
case (Lambda_sl_ID)
call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',&
'mean free path for slip','m')
case (threshold_stress_slip_ID)
call results_writeDataset(group,dst%tau_pass,'tau_pass',&
'passing stress for slip','Pa')
case (f_tw_ID)
call results_writeDataset(group,stt%f_tw,'f_tw',&
'twinned volume fraction','m³/m³')
case (Lambda_tw_ID)
call results_writeDataset(group,dst%Lambda_tw,'Lambda_tw',&
'mean free path for twinning','m')
case (tau_hat_tw_ID)
call results_writeDataset(group,dst%tau_hat_tw,'tau_hat_tw',&
'threshold stress for twinning','Pa')
case (f_tr_ID)
call results_writeDataset(group,stt%f_tr,'f_tr',&
'martensite volume fraction','m³/m³')
end select
enddo outputsLoop
end associate
#else
integer, intent(in) :: instance
character(len=*) :: group
@ -1130,7 +1163,7 @@ pure subroutine kinetics_slip(Mp,T,instance,of, &
tau_eff !< effective resolved stress
integer :: i
associate(prm => param(instance), stt => state(instance), dst => microstructure(instance))
associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
do i = 1, prm%sum_N_sl
tau(i) = math_mul33xx33(Mp,prm%P_sl(1:3,1:3,i))
@ -1203,7 +1236,7 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,&
integer :: i,s1,s2
associate(prm => param(instance), stt => state(instance), dst => microstructure(instance))
associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
do i = 1, prm%sum_N_tw
tau(i) = math_mul33xx33(Mp,prm%P_tw(1:3,1:3,i))
@ -1275,7 +1308,7 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,&
integer :: i,s1,s2
associate(prm => param(instance), stt => state(instance), dst => microstructure(instance))
associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
do i = 1, prm%sum_N_tr
tau(i) = math_mul33xx33(Mp,prm%P_tr(1:3,1:3,i))

22
src/plastic_isotropic.f90
View File

@ -410,8 +410,7 @@ subroutine plastic_isotropic_dotState(Mp,instance,of)
xi_inf_star = prm%xi_inf
else
xi_inf_star = prm%xi_inf &
+ asinh( (dot_gamma / prm%c_1)**(1.0_pReal / prm%c_2) &
)**(1.0_pReal / prm%c_3) &
+ asinh( (dot_gamma / prm%c_1)**(1.0_pReal / prm%c_2))**(1.0_pReal / prm%c_3) &
/ prm%c_4 * (dot_gamma / prm%dot_gamma_0)**(1.0_pReal / prm%n)
endif
dot%xi(of) = dot_gamma &
@ -470,7 +469,7 @@ function plastic_isotropic_postResults(Mp,instance,of) result(postResults)
c = c + 1
case (dot_gamma_ID)
postResults(c+1) = prm%dot_gamma_0 &
* (sqrt(1.5_pReal) * norm_Mp /(prm%M * stt%xi(of)))**prm%n
* (sqrt(1.5_pReal) * norm_Mp /(prm%M * stt%xi(of)))**prm%n
c = c + 1
end select
@ -485,23 +484,28 @@ end function plastic_isotropic_postResults
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_results(instance,group)
#if defined(PETSc) || defined(DAMASKHDF5)
use results
#if defined(PETSc) || defined(DAMASK_HDF5)
use results, only: &
results_writeDataset
implicit none
integer, intent(in) :: instance
character(len=*) :: group
integer, intent(in) :: instance
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(instance), stt => state(instance))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (xi_ID)
call results_writeDataset(group,stt%xi,'xi','resistance against plastic flow','Pa')
end select
enddo outputsLoop
end associate
#else
integer, intent(in) :: instance
character(len=*) :: group
integer, intent(in) :: instance
character(len=*), intent(in) :: group
#endif
end subroutine plastic_isotropic_results

26
src/plastic_kinematichardening.f90
View File

@ -551,8 +551,9 @@ end function plastic_kinehardening_postResults
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
subroutine plastic_kinehardening_results(instance,group)
#if defined(PETSc) || defined(DAMASKHDF5)
use results
#if defined(PETSc) || defined(DAMASK_HDF5)
use results, only: &
results_writeDataset
implicit none
integer, intent(in) :: instance
@ -562,6 +563,27 @@ subroutine plastic_kinehardening_results(instance,group)
associate(prm => param(instance), stt => state(instance))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (crss_ID)
call results_writeDataset(group,stt%crss,'xi_sl', &
'resistance against plastic slip','Pa')
case(crss_back_ID)
call results_writeDataset(group,stt%crss_back,'tau_back', &
'back stress against plastic slip','Pa')
case (sense_ID)
call results_writeDataset(group,stt%sense,'sense_of_shear','tbd','1')
case (chi0_ID)
call results_writeDataset(group,stt%chi0,'chi0','tbd','Pa')
case (gamma0_ID)
call results_writeDataset(group,stt%gamma0,'gamma0','tbd','1')
case (accshear_ID)
call results_writeDataset(group,stt%accshear,'gamma_sl', &
'plastic shear','1')
end select
enddo outputsLoop
end associate

82
src/plastic_none.f90
View File

@ -5,13 +5,13 @@
!> @brief Dummy plasticity for purely elastic material
!--------------------------------------------------------------------------------------------------
module plastic_none
implicit none
private
public :: &
plastic_none_init
implicit none
private
public :: &
plastic_none_init
contains
!--------------------------------------------------------------------------------------------------
@ -19,41 +19,39 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine plastic_none_init
use debug, only: &
debug_level, &
debug_constitutive, &
debug_levelBasic
use material, only: &
phase_plasticity, &
material_allocatePlasticState, &
PLASTICITY_NONE_label, &
PLASTICITY_NONE_ID, &
material_phase, &
plasticState
implicit none
integer :: &
Ninstance, &
p, &
NipcMyPhase
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>'
Ninstance = count(phase_plasticity == PLASTICITY_NONE_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
do p = 1, size(phase_plasticity)
if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NipcMyPhase = count(material_phase == p)
call material_allocatePlasticState(p,NipcMyPhase,0,0,0, &
0,0,0)
plasticState(p)%sizePostResults = 0
enddo
use debug, only: &
debug_level, &
debug_constitutive, &
debug_levelBasic
use material, only: &
phase_plasticity, &
material_allocatePlasticState, &
PLASTICITY_NONE_label, &
PLASTICITY_NONE_ID, &
material_phase, &
plasticState
implicit none
integer :: &
Ninstance, &
p, &
NipcMyPhase
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>'
Ninstance = count(phase_plasticity == PLASTICITY_NONE_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
do p = 1, size(phase_plasticity)
if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle
NipcMyPhase = count(material_phase == p)
call material_allocatePlasticState(p,NipcMyPhase,0,0,0, &
0,0,0)
plasticState(p)%sizePostResults = 0
enddo
end subroutine plastic_none_init

40
src/plastic_nonlocal.f90
View File

@ -2316,7 +2316,7 @@ outputsLoop: do o = 1,size(param(instance)%outputID)
case (rho_dot_ann_ath_ID)
postResults(cs+1:cs+ns) = results(instance)%rhoDotAthermalAnnihilation(1:ns,1,of) &
+ results(instance)%rhoDotAthermalAnnihilation(1:ns,2,of)
+ results(instance)%rhoDotAthermalAnnihilation(1:ns,2,of)
cs = cs + ns
case (rho_dot_ann_the_edge_ID)
@ -2402,8 +2402,9 @@ end function getRho
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
subroutine plastic_nonlocal_results(instance,group)
#if defined(PETSc) || defined(DAMASKHDF5)
use results
#if defined(PETSc) || defined(DAMASK_HDF5)
use results, only: &
results_writeDataset
implicit none
integer, intent(in) :: instance
@ -2413,6 +2414,39 @@ subroutine plastic_nonlocal_results(instance,group)
associate(prm => param(instance), stt => state(instance))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (rho_sgl_mob_edg_pos_ID)
call results_writeDataset(group,stt%rho_sgl_mob_edg_pos, 'rho_sgl_mob_edg_pos', &
'positive mobile edge density','1/m²')
case (rho_sgl_imm_edg_pos_ID)
call results_writeDataset(group,stt%rho_sgl_imm_edg_pos, 'rho_sgl_imm_edg_pos',&
'positive immobile edge density','1/m²')
case (rho_sgl_mob_edg_neg_ID)
call results_writeDataset(group,stt%rho_sgl_mob_edg_neg, 'rho_sgl_mob_edg_neg',&
'negative mobile edge density','1/m²')
case (rho_sgl_imm_edg_neg_ID)
call results_writeDataset(group,stt%rho_sgl_imm_edg_neg, 'rho_sgl_imm_edg_neg',&
'negative immobile edge density','1/m²')
case (rho_dip_edg_ID)
call results_writeDataset(group,stt%rho_dip_edg, 'rho_dip_edg',&
'edge dipole density','1/m²')
case (rho_sgl_mob_scr_pos_ID)
call results_writeDataset(group,stt%rho_sgl_mob_scr_pos, 'rho_sgl_mob_scr_pos',&
'positive mobile screw density','1/m²')
case (rho_sgl_imm_scr_pos_ID)
call results_writeDataset(group,stt%rho_sgl_imm_scr_pos, 'rho_sgl_imm_scr_pos',&
'positive immobile screw density','1/m²')
case (rho_sgl_mob_scr_neg_ID)
call results_writeDataset(group,stt%rho_sgl_mob_scr_neg, 'rho_sgl_mob_scr_neg',&
'negative mobile screw density','1/m²')
case (rho_sgl_imm_scr_neg_ID)
call results_writeDataset(group,stt%rho_sgl_imm_scr_neg, 'rho_sgl_imm_scr_neg',&
'negative immobile screw density','1/m²')
case (rho_dip_scr_ID)
call results_writeDataset(group,stt%rho_dip_scr, 'rho_dip_scr',&
'screw dipole density','1/m²')
case (rho_forest_ID)
call results_writeDataset(group,stt%rho_forest, 'rho_forest',&
'forest density','1/m²')
end select
enddo outputsLoop
end associate

51
src/plastic_phenopowerlaw.f90
View File

@ -563,28 +563,43 @@ end function plastic_phenopowerlaw_postResults
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
subroutine plastic_phenopowerlaw_results(instance,group)
#if defined(PETSc) || defined(DAMASKHDF5)
use results
#if defined(PETSc) || defined(DAMASK_HDF5)
use results, only: &
results_writeDataset
implicit none
integer, intent(in) :: instance
character(len=*) :: group
integer :: o
implicit none
integer, intent(in) :: instance
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(instance), stt => state(instance))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (resistance_slip_ID)
call results_writeVectorDataset(group,stt%xi_slip,'xi_slip','Pa')
case (accumulatedshear_slip_ID)
call results_writeVectorDataset(group,stt%gamma_slip,'gamma_slip','-')
end select
enddo outputsLoop
end associate
associate(prm => param(instance), stt => state(instance))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (resistance_slip_ID)
call results_writeDataset(group,stt%xi_slip, 'xi_sl', &
'resistance against plastic slip','Pa')
case (accumulatedshear_slip_ID)
call results_writeDataset(group,stt%gamma_slip,'gamma_sl', &
'plastic shear','1')
case (resistance_twin_ID)
call results_writeDataset(group,stt%xi_twin, 'xi_tw', &
'resistance against twinning','Pa')
case (accumulatedshear_twin_ID)
call results_writeDataset(group,stt%gamma_twin,'gamma_tw', &
'twinning shear','1')
end select
enddo outputsLoop
end associate
#else
integer, intent(in) :: instance
character(len=*) :: group
integer, intent(in) :: instance
character(len=*), intent(in) :: group
#endif
end subroutine plastic_phenopowerlaw_results

1534
src/results.f90
View File

File diff suppressed because it is too large Load Diff

6
src/rotations.f90
View File

@ -777,12 +777,12 @@ pure function qu2ax(qu) result(ax)
real(pReal) :: omega, s
omega = 2.0 * acos(math_clip(qu%w,-1.0_pReal,1.0_pReal))
! if the angle equals zero, then we return the rotation axis as [001]
if (dEq0(omega)) then
ax = [ 0.0, 0.0, 1.0, 0.0 ]
if (dEq0(sqrt(qu%x**2+qu%y**2+qu%z**2))) then
ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ]
elseif (dNeq0(qu%w)) then
s = sign(1.0_pReal,qu%w)/sqrt(qu%x**2+qu%y**2+qu%z**2)
omega = 2.0_pReal * acos(math_clip(qu%w,-1.0_pReal,1.0_pReal))
ax = [ qu%x*s, qu%y*s, qu%z*s, omega ]
else
ax = [ qu%x, qu%y, qu%z, PI ]